Journal of Harbin Institute of Technology (New Series)  2018, Vol. 25 Issue (6): 46-58  DOI: 10.11916/j.issn.1005-9113.18024 0

### Citation

Yong Wang, Shuai Feng, Qingxiang Zhang. Inverse Synthetic Aperture Radar Imaging of Multi-Targets Based on Image Processing Technique[J]. Journal of Harbin Institute of Technology (New Series), 2018, 25(6): 46-58. DOI: 10.11916/j.issn.1005-9113.18024.

### Fund

Sponsored by the National Natural Science Foundation of China(Grant No. 61622107)

### Corresponding author

Winner of the National Science Fund for Excellent Youth Scholars and the Program for New Century Excellent Talents in University of Ministry of Education. E-mail: wangyong6012@hit.edu.cn

### Article history

Inverse Synthetic Aperture Radar Imaging of Multi-Targets Based on Image Processing Technique
Yong Wang, Shuai Feng, Qingxiang Zhang
Research Institute of Electronic Engineering Technology, Harbin Institute of Technology, Harbin 150001, China
Abstract: This paper presents a new multi-targets inverse synthetic aperture radar (ISAR) imaging approach via the image segmentation processing. This method can separate multi-targets with similar velocities, and there is no strict limit on the rotational state of the targets. Firstly, the motion compensation for the completely multi-targets echo is carried out and the coarse image can be achieved with the Range-Doppler (RD) technique. Then a series of image processing methods and image segmentation processing are used to separate the echo data of each mono-target. At last, the image with high quality of each target can be achieved with the RD technique and the Range-Instantaneous-Doppler (RID) technique. ISAR imaging results of simulated and measured data validate the validity of the proposed approach.
Key words: ISAR     multi-targets     image segmentation
1 Introduction

Here, we present a new technique which can be applied to the situation that all the targets share the analogous translations and different oscillations. This hypothesis is consistent with the motion state of the multi-ship-targets on the sea surface [23]. The proposed method can obtain the coarse imaging of multi-targets based on the traditional mono-target imaging approaches. With the aim being to reduce the background noise and enable the sparse image to be processed by the traditional image segmentation methods, a series of image processing algorithms are implemented, including noise suppression, image interception based on the center of mass, interference fringe suppression, median filtering, and pixel expansion. After the processing in image domain, the signals of each target in the image domain will gather in the independent connected regions, and the image segmentation based on the geometric clustering can be used to extract different targets. Finally, we can use the echo signals of each target to achieve the high-resolution focused imaging results by the traditional mono-target imaging processing, including the RD and RID technique.

2 Imaging Geometry and Signal Model of Multi-Ship-Targets

The geometric model of multi-ship-targets is shown in Fig. 1. At the coherent processing interval (CPI), it is assumed that there are Q targets can be observed in the same antenna beam, different targets possess different velocities and three-dimensional oscillations. The blue lines O-UVW represent the target observation coordinate system, where the green lines represent the instantaneous radar antenna beam.

 Fig.1 Imaging geometry of multi-targets

It is assumed that the target q, q∈{1, 2, ..., Q} can be depicted as the summation of Kq scattering points, the instantaneous distance between the scattering point k, k∈{1, 2, ..., Kq} of target q to the radar at instant t is Rq, k(t). Generally, the transmitted signal is Linear-Frequency-Modulation (LFM) signal, which can be written as

 $s\left( {\hat t,{t_m}} \right) = rect\left( {\frac{{\hat t}}{{{T_p}}}} \right)\exp \left( {{\rm{j}}2{\rm{ \mathit{ π} }}\left( {{f_c}t + \frac{1}{2}\gamma {{\hat t}^2}} \right)} \right)$ (1)

where $rect\left[ x \right] = \left\{ \begin{array}{l} 1, \;\;\left| x \right| \le 0.5\\ 0, \;\;\left| x \right| > 0.5 \end{array} \right.$ denotes the rectangle window. Tp denotes the pulse width. fc and γ represent the central frequency and chirp rate, respectively. ${\hat t}$=t-m·T and tm=m·T represent the fast and slow time, while T denotes the pulse repetition interval (PRI). Then the received echoes of multi-targets can be illustrated as

 $\begin{array}{l} {s_r}\left( {\hat t,{t_m}} \right) = \sum\limits_{q = 1}^Q {\sum\limits_{k = 1}^{{K_q}} {\left\{ {{\sigma _{k,q}}rect\left[ {\left( {\hat t - \frac{{2{R_{q,k}}\left( t \right)}}{c}} \right)/{T_p}} \right] \cdot } \right.} } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left( {{\rm{j}}2{\rm{ \mathit{ π} }}\left[ {{f_c}\left( {t - \frac{{2{R_{q,k}}\left( t \right)}}{c}} \right) + } \right.} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {\frac{1}{2}\gamma {{\left( {\hat t - \frac{{2{R_{q,k}}\left( t \right)}}{c}} \right)}^2}} \right]} \right\} \end{array}$ (2)

where σk, q denotes the backscattering coefficient for the scattering k in the target q.

For the special characters of LFM signals, the Dechirping method is used to implement the pulse compression in this paper. The reference distance is assumed to be Rref, then the reference signal is written as

 $\begin{array}{l} {s_{{\rm{ref}}}}\left( {\hat t,{t_m}} \right) = rect\left[ {\left( {\hat t - \frac{{2{R_{{\rm{ref}}}}}}{c}} \right)/{T_{{\rm{ref}}}}} \right] \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left( {{\rm{j}}2{\rm{ \mathit{ π} }}\left[ {{f_c}\left( {t - \frac{{2{R_{{\rm{ref}}}}}}{c}} \right) + \frac{1}{2}\gamma {{\left( {\hat t - \frac{{2{R_{{\rm{ref}}}}}}{c}} \right)}^2}} \right]} \right) \end{array}$ (3)

where Tref is the pulse width for the reference signal, which satisfies Tref > Tp. Then, we can obtain the base band signal by

 ${s_{if}}\left( {\hat t,{t_m}} \right) = {s_r}\left( {\hat t,{t_m}} \right) \cdot s_{{\rm{ref}}}^ * \left( {\hat t,{t_m}} \right)$ (4)

Then, the base band signal is shown as

 $\begin{array}{l} {s_{if}}\left( {\hat t,{t_m}} \right) = \sum\limits_{q = 1}^Q {\sum\limits_{k = 1}^{{K_q}} {\left\{ {{\sigma _{k,q}}rect\left[ {\left( {\hat t - \frac{{2{R_{q,k}}\left( t \right)}}{c}} \right)/{T_p}} \right] \cdot } \right.} } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{{\rm{e}}^{ - {\rm{j}}\frac{{4{\rm{ \mathit{ π} }}}}{c}{R_\mathit{\Delta }}\gamma \left( {\hat t - \frac{{2{R_{{\rm{ref}}}}}}{c}} \right)}}{{\rm{e}}^{ - {\rm{j}}\frac{{4{\rm{ \mathit{ π} }}}}{c}{f_c}{R_\mathit{\Delta }}}}{{\rm{e}}^{{\rm{j}}\frac{{4{\rm{ \mathit{ π} }}}}{{{c^2}}}\gamma R_\mathit{\Delta }^2}}} \right\} \end{array}$ (5)

where RΔ=Rq, k(t)－Rref. Perform the Fourier transform to Eq.(5), we can get

 $\begin{array}{l} {S_{if}}\left( {{f_i},{t_m}} \right) = \sum\limits_{q = 1}^Q {\sum\limits_{k = 1}^{{K_q}} {\left\{ {{\sigma _{k,q}}{T_p}{\mathop{\rm sinc}\nolimits} \left[ {{T_p}\left( {{f_i} + 2\frac{\gamma }{c}{R_\mathit{\Delta }}} \right)} \right] \cdot } \right.} } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\exp \left( { - {\rm{j}}\frac{{4{\rm{ \mathit{ π} }}}}{c}{f_c}{R_\mathit{\Delta }}} \right)\exp \left( {{\rm{j}}\frac{{4{\rm{ \mathit{ π} }}}}{{{c^2}}}\gamma R_\mathit{\Delta }^2} \right)} \right\} \end{array}$ (6)

where σk, qTpsinc[Tp(fi+$2\;\frac{\gamma }{c}{R_\mathit{\Delta }}$)] donates the envelope of the range profile, while the RΔ shows the translation of the envelope.

Due to the diversity of Rq(t) of different targets, most of the motion compensation method that manifest as the whole move of the range profile, can hardly accomplish the translational compensation of each target simultaneously. This is the essential reason why the traditional mono-target motion compensation method performs poorly in a multi-target imaging situation.

In addition, the ships sailing on sea surface also possess complex oscillations in three-dimensions[24-25]: yaw, pitch and roll, as shown in Fig. 1. According to the results in Ref.[26], the three-dimensional oscillation for the ship target can be approximated as the periodic cosine function. The angular velocities in the three directions can be shown as

 ${\theta _l}\left( t \right) = {A_l}\cos \left( {{w_l}t + {\chi _l}} \right),\left( {l = yaw,pitch,roll} \right)$ (7)

where Al represents the maximum angular motion amplitude. wl=2π/Tl represents the angular velocity of oscillation while Tl donates cycle of oscillation. χl represents the initial angle. Then the three-dimensional oscillation matrix of the target can be expressed as

 $\begin{array}{l} \mathit{\boldsymbol{roll}}\left( {{\theta _r}\left( t \right)} \right) = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos {\theta _r}\left( t \right)}&{ - \sin {\theta _r}\left( t \right)}\\ 0&{\sin {\theta _r}\left( t \right)}&{\cos {\theta _r}\left( t \right)} \end{array}} \right]\\ \mathit{\boldsymbol{yaw}}\left( {{\theta _y}\left( t \right)} \right) = \left[ {\begin{array}{*{20}{c}} {\cos {\theta _y}\left( t \right)}&{ - \sin {\theta _y}\left( t \right)}&0\\ {\sin {\theta _y}\left( t \right)}&{\cos {\theta _y}\left( t \right)}&0\\ 0&0&1 \end{array}} \right]\\ \mathit{\boldsymbol{pitch}}\left( {{\theta _p}\left( t \right)} \right) = \left[ {\begin{array}{*{20}{c}} {\cos {\theta _p}\left( t \right)}&0&{\sin {\theta _p}\left( t \right)}\\ 0&1&0\\ { - \sin {\theta _p}\left( t \right)}&0&{\cos {\theta _p}\left( t \right)} \end{array}} \right] \end{array}$ (8)

Moreover, the synthetic oscillation matrix is shown as

 $\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Rot}}\left( {{\theta _r}\left( t \right),{\theta _p}\left( t \right),{\theta _y}\left( t \right)} \right) = \mathit{\boldsymbol{roll}}\left( {{\theta _r}\left( t \right)} \right) \cdot }\\ {\mathit{\boldsymbol{pitch}}\left( {{\theta _p}\left( t \right)} \right)\mathit{\boldsymbol{yaw}}\left( {{\theta _y}\left( t \right)} \right)} \end{array}$ (9)

The ship coordinate system (O, ξ, η, ζ) is established at the origin of the ship oscillation center as shown in Fig. 1. We suppose that the initial coordinate for the point k is (ξ0k, η0k, ζ0k) in the coordinate system (O′, ξ, η, ζ). Then, the coordinate of the scattering point k in the coordinate system O-UVW at the instant t can be written as

 $\begin{array}{l} {\left[ {{U^k}\left( t \right),{V^k}\left( t \right),{W^k}\left( t \right)} \right]^{\rm{T}}} = \mathit{\boldsymbol{Rot}}\left( {{\theta _r}\left( t \right),{\theta _p}\left( t \right),} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{\theta _y}\left( t \right)} \right){\left[ {\xi _0^k,\eta _0^k,\zeta _0^k} \right]^{\rm{T}}} \end{array}$ (10)

Assume ${{\overrightarrow{\mathit{\boldsymbol{i}}}}_{U}},{{\overrightarrow{\mathit{\boldsymbol{i}}}}_{V}},{{\overrightarrow{\mathit{\boldsymbol{i}}}}_{W}}$ represent the unit vector of the three direction of coordinate axis in O-UVW, (U0, V0, W0) is the initial coordinate of the target center in O-UVW, then ${{\overrightarrow{\mathit{\boldsymbol{R}}}}_{k,q}}\left( t \right)$ can be expressed as

 $\begin{array}{l} {{\mathit{\boldsymbol{\vec R}}}_{k,q}}\left( t \right) = U\left( t \right){{\mathit{\boldsymbol{\vec i}}}_U} + V\left( t \right){{\mathit{\boldsymbol{\vec i}}}_V} + W\left( t \right){{\mathit{\boldsymbol{\vec i}}}_W} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{U_0} + \int_0^t {{v_U}\left( \tau \right){\rm{d}}\tau } + {U^k}\left( t \right)} \right]{{\mathit{\boldsymbol{\vec i}}}_U} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{V_0} + \int_0^t {{v_v}\left( \tau \right){\rm{d}}\tau } + {V^k}\left( t \right)} \right]{{\mathit{\boldsymbol{\vec i}}}_V} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{W_0} + \int_0^t {{v_W}\left( \tau \right){\rm{d}}\tau } + {W^k}\left( t \right)} \right]{{\mathit{\boldsymbol{\vec i}}}_W} \end{array}$ (11)

where vU(t), vV(t), vW(t) represent the velocities in the three direction of coordinate axis U, V, W, respectively.

As we can see in Eq.(11), due to the time-varying velocity and the three-dimensional oscillation, the rotation of the target could be very complex. In this situation, the classical RD algorithm could cause serious defocus on the azimuth direction [27-28]. As a result, the RID technique is widely used to solve this problem. The most commonly used RID algorithms are based on the LFM signal model[29-35], such as the Wigner-Hough transform, the match Fourier transform, the Dechirp method and the signal decomposition technique. However, if the target has complicated rotational movement, just like the ship model above, the high-order phase term will come into being. In this situation, the RID technique based on the cubic phase signal model will behave better[36-37], such as the PHMT algorithm and the PGCPF algorithm.

3 Processing in Image Domain

The detailed steps of the imaging process are introduced in this section.

(a) Coarse imaging of multi-targets.

For the separation in image domains, the first step in common is to obtain the coarse imaging results. In this paper, the cross-correlation method and constant phase-error elimination method[9] are used to the envelope alignment and phase correction, respectively.

1) Cross-correlation method.

In ISAR imaging, the frequency of pulse repetition is generally high, which makes the envelopes of two adjacent range profiles possess strong correlation. The cross-correlation method calculates the cross-correlation function of the adjacent envelopes and find the peak position, and then utilizes the delay corresponding to the peak to process the envelope alignment. However, due to small error accumulation and single abnormal echo, the traditional cross-correlation method may cause envelope drift and jump error. To avoid this problem, a novel method based on accumulation has been proposed. This method replaces the single echo with the accumulation of N time's echoes, and achieves the optimal alignment in the whole. The cross-correlation function in this method is

 $R\left( s \right) = \int {\left( {\sum\limits_{i = 1}^{N - 1} {e_i^ * \left( r \right)} } \right){e_N}\left( {r + s} \right){\rm{d}}r}$ (12)

where ei(r) represents the real envelope of the ith echo after pulse compression. The jump error can be avoided and the envelope drift can be effectively reduced after appropriate accumulation.

2) Constant phase-error elimination method.

After the envelope alignment, the amplitudes of the echoes are corrected accurately, but the phases are still chaotic. The constant phase-error elimination method is a simplified algorithm to process phase correction for practical engineering applications. Assume the sub-echo envelope of the pth range unit can be simplified as

 ${s_p}\left( m \right) = {\sigma _{1p}}\left( m \right){{\rm{e}}^{{\rm{j}}\left( {{\varphi _{1p0}} + \frac{{4{\rm{ \mathit{ π} }}}}{\lambda }m{\chi _{1p}} + {\psi _{1p}}\left( m \right) + {\gamma _m}} \right)}}$ (13)

where σ1p(m) and ψ1p(m) represent the modulation of amplitude and phase caused by small clutter and noise, respectively. φ1p0 represents the initial phase. (4π/λ)1p represents the phase component due to the Doppler frequency which changes over the slow time. γm represents the initial phase error.

Multiply the mth echo shown in Eq.(13) with the conjugate of the (m-1)th echo, we can get

 $\begin{array}{*{20}{c}} {{s_p}\left( m \right)s_p^ * \left( {m - 1} \right) = {\sigma _{1p}}\left( m \right){\sigma _{1p}}\left( {m - 1} \right) \cdot }\\ {{{\rm{e}}^{{\rm{j}}\left( {\frac{{4{\rm{ \mathit{ π} }}}}{\lambda }{\chi _{1p}} + \mathit{\Delta }{\psi _{1p}}\left( m \right) + \mathit{\Delta }{\gamma _m}} \right)}}} \end{array}$ (14)

where

 $\mathit{\Delta }{\psi _{1p}}\left( m \right) = {\psi _{1p}}\left( m \right) - {\psi _{1p}}\left( {m - 1} \right)$

and Δγm=γmγm－1 represent the difference between adjacent phase fluctuations and the difference between adjacent initial phase error, respectively. Obviously, the initial phase φ1p0 has been dropped and the Doppler phase has been converted to the constant (4π/λ)χ1p. Then the phase required for correction can be calculated by

 $\mathit{\Delta }{{\tilde \gamma }_m} = \arg \left[ {\sum\limits_{p = 1}^L {{s_p}\left( m \right)s_p^ * \left( {m - 1} \right)} } \right]$ (15)

3) RD algorithm for coarse imaging.

After the holistic motion compensation, the range profiles have been aligned roughly. The range profiles can hardly align strictly, because each target has different translational motion. Hence, it is meaningless to use complex accurate imaging method to obtain the coarse image. The RD algorithm is simple, easy to realize by computer and reversible. In addition, due to the complex three-dimensional oscillation of each target, the azimuth direction will create severe defocusing.

In a word, the coarse image obtained by the above method would not focus accurately. However, the signals of different targets will assemble into independent regions. Besides, due to the character of ISAR imaging and complex electromagnetic environment, the pixels in the target area are sparse and the strong background noise will spread throughout the image. The traditional image segmentation methods cannot divide images in this case. A series of image processing methods are used to solve this problem, the specific steps are shown as follows.

(b) Global noise suppression.

The first step of the image processing methods is the global noise suppression. Due to the sea-clutter and the interference scattering points that randomly distribute in the environment, the noise distribution of ISAR image is often complex and random. To suppress the background noise adaptively, this paper propose a novel global noise suppression method by referring to the two-channel CFAR method in the signal detection. The statistical property of the noise is obtained by the elimination of the signal, and the threshold of adaptive detection is calculated by the mean and variance of noise. The detailed steps are:

Step1   Compute the mean μ0 and variance δ02 of the normalized image, then get the initial threshold V0=μ0+g0·δ0.

Step 2   Remove the points with the energy higher than the threshold, and calculate the mean μ0 and variance δ02 of the new image.

Step 3   Calculate the relative errors

 $\alpha = \frac{{{\mu _i} - {\mu _{i - 1}}}}{{{\mu _{i - 1}}}},\beta = \frac{{\delta _i^2 - \delta _{i - 1}^2}}{{\delta _{i - 1}^2}}$

where μi and δi2 represent the mean and the variance of the image after the ith iteration, respectively.

Step 4   Terminate the iteration and enter Step 5 when α < ε1 and β < ε2, otherwise calculate the new threshold by Vi=μi+g0·δi and return to Step 2.

Step 5   Calculate the final threshold as VT=μi+gT·δi and set the value of the pixel in image below that threshold to be zero.

In the above steps, ε1 and ε2 are the default error values, which are set to 0.02 and 0.35 respectively by experience.g0 and gT are the threshold factors, which are set to 3 and 3.7 respectively by experience.

(c) Adaptive image intercept.

After the global noise suppression, a large number of useless white backgrounds will appear in the image, these backgrounds will add unnecessary computational burden to the image processing. An effective way to remove the blank background is to intercept the target area. This paper applies a method of adaptive interception based on the calculation of mass center. The specific steps are:

Step 1   Calculate the zero and first order moment of the image

 $\left\{ \begin{array}{l} {m_{00}} = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {I\left( {i,j} \right)} } \\ {m_{10}} = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {i \cdot I\left( {i,j} \right)} } \\ {m_{01}} = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {j \cdot I\left( {i,j} \right)} } \end{array} \right.$ (16)

where I(i, j) denotes the pixel intensity, M×N is the size of image.

Step 2   Obtain the coordinate of mass center by the zero order moment and the first moment

 ${i_c} = \frac{{{m_{10}}}}{{{m_{00}}}},\;\;\;\;{j_c} = \frac{{{m_{01}}}}{{{m_{00}}}}$ (17)

Step 3   Take the coordinate as the center, and take an image of the size K×K in the original image. In this paper, K is set to 400.

(d) Fringe and spot noise suppression.

Fringe and spot noise may arise from the strong interference point in the background, or the overlap of the secondary flaps of strong scattering spots. The intensity of these noises is higher than some weak scattering points of the targets, so the global threshold cannot be used to suppress them. Strong spot noise is represented as a small set of pixels in the ISAR image, and the pixel mean of the two-dimensional neighborhood is generally low. Strong stripe noise is generally not continuous in the range unit or the Doppler unit, so the pixel mean in the one-dimensional neighborhood in this direction is weak. In contrast, the scattering points of target tend to be continuous, and the pixel mean in each neighborhood is large. According to the difference, this paper uses the following algorithm to suppress the fringe interference and spot noise:

Step 1   Set reasonable center pixel threshold Thc and 2-D neighborhood threshold Th2. The values of two thresholds will directly affect the noise suppression effect, and the thresholds should be increased appropriately when the noise suppression effect is not obvious. In this paper, the thresholds Thc and Th2 are set to 0.8 and 0.4 respectively by experience.

Step 2   Iterate through all the pixels, find the pixels that satisfy

 $\left\{ \begin{array}{l} I\left( {i,j} \right) < T{h_c}\\ \left[ {\sum\limits_{a = i - 1}^{i + 1} {\sum\limits_{b = j - 1}^{j + 1} {I\left( {a,b} \right) - I\left( {i,j} \right)} } } \right]/8 < T{h_2} \end{array} \right.$ (18)

and set I(i, j) to zero.

After the fringe and spot noise suppression, the edges of the target area are clearer, and the probability of misconnecting between different target areas is also reduced.

(e) Image cavity filling.

After the noise suppression, there are still plenty of holes in the target area. In order to focus the signal of the mono-target scattering points in a connected independent region, these holes need to be filled in.

The main steps of cavity filling are median filtering, pixel expansion, image smoothing, and the binary processing. The median filtering is used to suppress isolated strong scattering points, while protecting the edge and shape of the targets as completely as possible. The pixel expansion is used to fill the holes in the area of targets. The filling effect is directly related to the expansion ratio. The greater the expansion ratio, the better the effect, but too much expansion will connect the different target regions that are originally independent. Therefore, it is very important to choose appropriate expansion ratio. Image smoothing and binary processing can further fill the holes in the target area without undue influence on the outer edge of the target. Since these steps are commonly used in image processing, the detailed principles are not covered in this paper.

(f) Image segmentation.

The image segmentation method based on geometric clustering is used in this paper. The principle of this method is to detect all connected independent regions in the image and then label them consequently. The number of pixels in the connected region determines the order of the label. First, we can select the area of the largest number of pixels to be the region of the target, and set the point outside of the area to be zero. Then, the outline of the selected target area can be obtained. Due to the small expansion radio in the image cavity filling, in the target area there may still exist some holes. The ISAR imaging process shows that the signal of mono-target will gather in a contiguous area in the image domain. The holes can only indicate that the signal energy is low and does not mean that the signals do not exist. To ensure the integrity of the extracted signal, the holes in the target contour should be all filled.

The reasonable filling methods should ensure that the external contour of the target is unchanged, so the image expansion and other algorithms are not applicable here. The endpoint connection method is used to solve the internal cavity problem. The specific steps are:

Step 1   Go through each row in the image, find the first pixel and the last pixel of the image in each row, and set the pixel values between the two pixels to one. Then the row-connection image A has been obtained.

Step 2   Go through each column in the image, find the first pixel and the last pixel of the image in each column, and set the pixel values between the two pixels to one. Then the column-connection image B has been obtained.

Step 3   Multiply the corresponding pixels of A and B, and get the composite connection image C.

This method is similar to the strategy in WeiQi(or go chess). Only the blank areas surrounded by pixel points will be filled, and the blank areas connected to the external blank background will remain unchanged.

4 Imaging Process Based on Image Segmentation

This section will describe the overall process of the proposed methods systematically:

Step 1   Obtain the range profiles by the pulse compression of the range direction. Then, use the method (a) in Section 3 to get the coarse image of the multi-targets.

Step 2   Obtain the external contours of the separating targets by using the processing of image domain in (b) to (f) in Section 3.

Step 3   Extract the separating signals of the targets by the external contours. The extraction method is similar to the band-pass filtering. Set the coordinates of the target area to the passband, and the outer coordinates of the zone are the stopband, the filter in image domain can be obtained. By multiplying the filter and the coarse image, the signals of separating targets can be obtained.

To enable the signals of the selected targets to be extracted as completely as possible, the filter can be expanded before filtering, which is equivalent to amplify the passband slightly.

Step 4   Return the extracted signals to the data domain and then perform the motion compensation and imaging processing of mono-target. Because of the irregular 3-D oscillation of the ship, the RD algorithm may not be able to get clear imaging results. To solve this problem, a RID imaging method based on PGCPF is used in this paper. The received signals are modeled as cubic phase signals, and the parameters of the signals can be estimated by PGCPF algorithm. The PGCPF of s(n) is defined in Ref.[37] as:

 $\begin{array}{l} PGCPF\left( \beta \right) = \prod\limits_{l = 1}^L {\sum\limits_{m = 0}^{\left( {N - 1} \right)/2} {s\left( {n + m} \right)s\left( {n - m} \right){s^ * }} } \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( { - n + m} \right){s^ * }\left( { - n - m} \right){s^ * }{{\rm{e}}^{ - {\rm{j}}\beta n{m^2}}} \end{array}$ (19)

After image processing, we can obtain the RD imaging results and the instantaneous image of each moment.

Fig. 2 is the flow diagram for the proposed method.

 Fig.2 Flow diagram of multi-targets ISAR imaging

5 ISAR Imaging Results

In this section, the effectiveness of the proposed method are testified by the simulated and measured data.

(a) Simulated data

The carrier frequency of the LFM signal is 5.52 GHz with the bandwidth is 400 MHz. The Pulse Repetition Interval (PRI) is 2.5 ms. The SNR of the received echo is 10 dB. The motional parameters of the targets are shown in Table 1. Fig. 3 shows the simulated target model, and each ship consists of 386 scattering points.

Table 1 Motion parameters of targets

 Fig.3 Simulated target model

Fig. 4 is the normalized coarse image via the method proposed in Section 3. The targets are faintly visible and blurred severely in the azimuth direction. In addition, the influence of environmental noise is obvious in the coarse image. Fig. 5 is the image after the global noise suppression. The background noise has been suppressed well after the processing. At the same time, the target areas are conserved integrally.

 Fig.4 Normalized coarse image

 Fig.5 Image after global noise suppression

Fig. 6 (a)-(f) are the results of a series of processing in the image domain. Fig. 6(a) is the image after adaptive image intercept, the target area is entirely intercepted. Fig. 6(b) is the image after fringe and spot noise suppression, the isolated strong spot noise and fringe interference are suppressed, and the target is retained perfectly. Fig. 6(c)-(f) are the images after the image cavity filling. Fig. 6(c) is the image after median filtering, the holes in the target area are initially filled while the edge and shape of the targets are protected as completely as possible. Fig. 6(d) is the image after pixel expansion, the holes in the target area are nearly filled, and the areas of different targets still stay independent because of the small expansion ratio. Fig. 6(e) is the image after the image smoothing, the holes in the target areas are further filled without undue influence on the outer edge. Fig. 6(f) is the image after the binary processing, the different targets become independently connected regions. In addition, the boundaries are clear enough to be segmented.

 Fig.6 Results of processing in image domain

Fig. 7 (a)-(d) show the results of image segmentation based on geometric clustering and endpoint connection. The separating targets are fully extracted, and there are no holes in the target area.

 Fig.7 Results of image segmentation

Fig. 8(a)-(d) are the radar images of the Target 1-4 via the RD algorithm. Compared with the coarse image shown in Fig. 4 and the results of processing in image domain shown in Fig. 6, the imaging quality has been greatly improved. This indicates that the important reason for the low quality of ISAR image is the difficulty of accurate motion compensation of each target. Although, the image results can be recognized as a ship, there is still some defocus. In order to estimate the parameters as accurately as possible, this paper considers the multi-component cubic phase signal as the signal model of the echo signals in same range cell. The imaging results in Fig. 9 are obtained by the PGCPF algorithm at selected time positions. Compared with the images in Fig. 8, the defocusing phenomenon no longer appears, and results in Fig. 9 can accurately reflect the external characteristics of the targets.

 Fig.8 Radar image via RD algorithm

 Fig.9 Radar image via RID algorithm

In order to evaluate quantitatively, the image entropy is used as the numerical evaluation index. The image entropy En(I) is defined as follows:

 $\begin{array}{l} En\left( I \right) = - \sum\limits_M {\sum\limits_N {P\left( {i,j} \right)\ln \left( {P\left( {i,j} \right)} \right)} } \\ P\left( {i,j} \right) = I{\left( {i,j} \right)^2}/\sum\limits_M {\sum\limits_N {I{{\left( {i,j} \right)}^2}} } \end{array}$ (20)

where I (i, j) denotes the image intensity, M×N is the size of image. The image entropy of different image results are shown in Table 2. The image entropies of the RD imaging results are much lower than the coarse image, which means the imaging quality has been improved after the target separation and refocusing. In addition, the lowest entropy of the RID image proves the effectiveness of the RID imaging algorithm.

Table 2 The image entropy of different image results

From the processing results of simulation data, the proposed separation method in the image domain is practical and effective, and has certain adaptability. Furthermore, the multi-ship-targets can be accurately focused by the RID imaging algorithm.

(b) Measured data.

To testify the effectiveness of the proposed method in the actual application scenario, a set of measured double-ship ISAR data is selected. The radar works at X band with bandwidth is 400 MHz and the PRI is 2 ms.

Fig. 10 is the normalized coarse image of measured data, and the background noise is very complicated. Fig. 11 is the image after global noise suppression. Most of the background noise has been suppressed well after the processing except for a small amount of strong spot noise. The results of a series of processing in the image domain before image segmentation are shown in Fig. 12.

 Fig.10 Normalized coarse image

 Fig.11 Image after global noise suppression

 Fig.12 Results of processing in image domain

Similar to the processing effect of the simulated data, the different targets become independently connected regions and the boundaries are clear enough to be segmented. Fig. 13 (a)-(b) shows the results of image segmentation based on geometric clustering and endpoint connection. The main bodies of the separating targets are extracted completely, but there is still a small fractured part of Target 1, which is not recognized as the target.

 Fig.13 Results of image segmentation

Fig. 14 (a)-(b) are the imaging results of the targets via the RD algorithm. For the Target 2, the imaging quality is relatively high due to the suitable imaging angle and the small three-dimensional oscillations. In the meantime, the Target 1 is almost unrecognizable because of the poor imaging angle and the serve defocusing phenomenon in the azimuth direction. Fig. 15(a)-(b) are the imaging results by the PGCPF algorithm at selected time positions. For Target 2, the image quality is slightly improved. As for Target 1, due to the bad imaging angle of Target 2, no suitable instantaneous imaging instance can be found in the CPI. The image entropy of different image results are shown in Table 3. By the entropy value, we can obtain the conclusion which is same as the visual judgment.

 Fig.14 Radar image via RD algorithm

 Fig.15 Radar image via RID algorithm

Table 3 The image entropy of different image results

6 Conclusions

A novel method for multi-targets ISAR imaging based on image segmentation is proposed in this paper. This method can be applied to multi-targets ISAR imaging with similar motion. After coarse imaging, processing in the image domain, image segmentation and the target extraction, the signal of each target can be obtained. Then, the well-focused radar image can be obtained by utilizing the imaging algorithm of mono-target. The ISAR imaging results of the simulated and measured data have testified the effectiveness of the proposed method. Nevertheless, the proposed method relies on the quality of coarse image. When the image quality of the coarse image is poor, for example, the azimuth exists a large number of false targets; the performance of the algorithm will be poor. This issue is the focus of future research.

References