PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of sensorsMDPI Open Access JournalsMDPI Open Access JournalsThis articleThis JournalInstructions for authorssubscribe
 
Sensors (Basel). 2017 February; 17(2): 366.
Published online 2017 February 14. doi:  10.3390/s17020366
PMCID: PMC5336065

3D Imaging of Rapidly Spinning Space Targets Based on a Factorization Method

Jonathan Li, Academic Editor

Abstract

Three-dimensional (3D) imaging of space targets can provide crucial information about the target shape and size, which are significant supports for the application of automatic target classification and recognition. In this paper, a new 3D imaging of space spinning targets via a factorization method is proposed. Firstly, after the translational compensation, the scattering centers two-dimensional (2D) range and range-rate sequence induced by the target spinning is extracted using a high resolution spectral estimation technique. Secondly, measurement data association is implemented to obtain the scattering center trajectory matrix by using a range-Doppler tracker. Then, we use an initial coarse angular velocity to generate the projection matrix, which consists of the scattering centers range and cross-range, and a factorization method is applied iteratively to the projection matrix to estimate the accurate angular velocity. Finally, we use the accurate estimate spinning angular velocity to rescale the projection matrix and the well-scaled target 3D geometry is reconstructed. Compared to the previous literature methods, ambiguity in the spatial axes can be removed by this method. Simulation results have demonstrated the effectiveness and robustness of the proposed method.

Keywords: 3D imaging, spinning target, factorization, scaling

1. Introduction

Because of its abundant information about the target’s structure and size, three-dimensional (3D) radar imaging of space targets plays a vital role in the space target recognition and classification field. Recently, the research on 3D imaging of space targets, such as space debris, ballistic targets, satellites and so on, has drawn intensive attention [1,2,3,4]. Here by 3D imaging we mean estimating the 3D coordinates of the target scattering centers from the target echoes. Up to now, various techniques have been proposed to generate high resolution 3D images of space targets using wideband inverse synthetic aperture radar (ISAR).

According to the number of sensors, the 3D imaging technique can be roughly classified into two categories [5]. The first category is interferometric ISAR imaging [2,6,7,8,9]. By using conventional two-dimensional (2D) ISAR imaging algorithms, the target echoes received by the different antennas are processed to form 2D range-Doppler (RD) images. Then the height of the scattering centers is extracted via the phase difference between the 2D images. These conventional 2D imaging algorithms assume the scattering center echo signal energy can be focused in one range and the Doppler resolution cell during the imaging time [8,9]. However, for rapidly spinning targets, there may be migration in the range and Doppler resolution cells, which results in a smeared 2D ISAR image. Therefore, how to achieve well-focused images by the conventional algorithms due to the time-varying range and Doppler information is a challenge [10].

The second one is based on monostatic ISAR, which exploits one-dimensional (1D) range series or 2D range-Doppler image sequences [11,12,13,14,15,16,17] to reconstruct the 3D distribution of the scattering centers. In [1,11,12,13,14,15] the singular value decomposition (SVD) is applied to the 1D range-only matrix to extract the 3D coordinates of scattering centers. This method requires the target to have a 3D rotational motion after the translational motion compensation, which is not applicable to a spinning target. In [16], GRT-CLEAN is proposed to form the 3D image of the spinning target. It is worth mentioning that the 3D images achieved via this algorithm are modified by a scaling factor in the cylindrical coordinate system. Mayhan et al. [17] assumed the target motion parameters are known and proposed the “snapshot” 3D imaging method which employs a sequence of range and range-rate data to extract the scattering centers’ 3D coordinates to obtain unambiguous scattering center coordinates. In a practical scenario, the target motion parameters are usually unknown and need to be estimated before using this method. With the development of the new reconstruction algorithm in the computer vision community [18,19], McFadden [20] and Li et al. [21] applied the factorization method to 2D RD image sequences to obtain the 3D geometry reconstructions based on the 3D-to-2D orthographic projection between the 3D geometry and 2D image sequence. Because of the unknown spinning angular velocity, the reconstructed points also have a certain scale ambiguity in the spatial axes.

The abovementioned problems in 3D imaging of spinning targets motivate our research. In this paper, a novel 3D imaging algorithm for rapidly spinning space targets based on the factorization method is presented. Compared to currently available methods, it can remove the 3D image scaling ambiguity in the shape matrix. Firstly, after the translational compensation of the echo data, the sequence of the independent 2D high resolution scattering center range and range-rate is extracted using 2D spectral estimation techniques [22]. Compared to the conventional Fourier Transform (FT) method, this technique offers two advantages: (1) it provides enhanced resolution both in range and range-rate; (2) it produces a scattering center trajectory matrix consisting of range and range-rate after implementing the scattering centers correlation developed later in this work. Secondly, the scattering centers association is implemented based on the Kalman filter and minimum Euclidean distance criterion [23]. Thirdly, by using a coarse angular velocity, the scattering centers projection matrix consisting of scattering centers’ range and cross-range is obtained. Then, a more accurate target spinning angular velocity estimate is achieved by applying the factorization method to the projection matrix. After that, the estimate angular velocity is used to rescale the projection matrix and the scaled projection matrix is factorized to obtain the newly accurate angular velocity estimate iteratively. Finally, the accurate angular velocity estimate is used to rescale the projection matrix and then factorize the scaled projection matrix. The advantage of the proposed method is that it achieves a well-scaled 3D image of the target.

The remainder of this paper is organized as follows: After an introduction, the 3D imaging geometry of the spinning target is introduced in Section 2. Subsequently, Section 3 presents the range and range-rate acquisition process. Then, the scattering center association process is described in Section 4. Section 5 describes the 3D imaging and scaling procedure in detail. The simulation results to verify the effectiveness and robustness of the proposed method are presented in Section 6. Finally, Section 7 draws the conclusion.

2. Imaging Geometry

In this section, the 3D imaging geometry of space spinning targets will be derived and analyzed in detail. For the applications considered here, we assume that the observed space target is spinning around one major axis during the observation interval [16,24]. Figure 1 shows the 3D imaging geometry after the coherent pre-processing and translational compensation [24]. Assume the target local coordinates system is O−XYZ and O is the origin. The space target spins around the Z axis and the angular velocity ω=(0,0,Ω). After the translational compensation, the unit vector of the RLOS (radar line of sight) is a constant [4]. Suppose that the azimuth and elevation angle of the RLOS are [var phi] and θ, respectively. The unit vector of the RLOS can be written as n^=(cosθcosϕ,cosϕsinϕ,sinθ).

Figure 1
3D imaging geometry of a spinning target.

When the radar transmits a high frequency wideband signal, the target can be represented by a reasonable number of fiducial non-coplanar dominant scattering centers. Suppose that the target is rigid and consists of S point scattering centers, whose positions are denoted as Ps=(xs,ys,zs)T,(s=1,,S), So time instant t, after the coherent pre-processing and translational compensation [25], the scattering center instantaneous radial range rs(t) and range-rate r˙s(t) induced by the spinning motion can be derived as follows:

rs(t)=n^tPs
(1)

r˙s(t)=n^(ω×(tPs))=(n^×ω)tPs
(2)

where t3×3 is called the 3D rotation matrix [26] and here, t=[cos(Ωt)sin(Ωt)0sin(Ωt)cos(Ωt)0001].

From Equations (1) and (2), the range and range-rate (meter-meter/s) sequence can be obtained. To obtain the 2D coordinates in the range and cross-range domain (meter-meter):

u(t)=rs(t)=[cosθcos(Ωtφ)cosθsin(Ωtφ)sinθ]Ps
(3)

v(t)=r˙s(t)Ωcosθ=[sin(Ωtφ)cos(Ωtφ)0]Ps
(4)

One can note that Equations (3) and (4) are the projection equations in range and cross-range dimensions, respectively. For the convenience in the following discussion, the corresponding range projection vector and cross-range projection vector are defined as follows:

{it=[cosθcos(Ωtφ)cosθsin(Ωtφ)sinθ]jt=[sin(Ωtφ)cos(Ωtφ)0]
(5)

Furthermore, note that the two vectors meet the following conditions:

{it=1jt=1itjt
(6)

According to Equations (3)–(5), the projection can be written in matrix form as follows:

W=MP
(7)

where the motion matrix M=[i1,,iM;j1jm]T, shape matrix P=[P1,,PS], and the projection matrix or the measurement matrix is:

W=[u11u12u1SuM1uM2uMSv11v12v1SvM1vM2vMS]
(8)

According to Equation (6), Equation (7) is an orthographic projection and the measurement matrix is the product of the motion matrix M2M×3 and the shape matrix P3×S. Here, this paper focuses on recovering the space spinning target shape matrix I from the measurement matrix W via the factorization method [18]. To begin the 3D imaging procedure, the first step is to extract the scattering centers range and range-rate data set, which will be discussed in the next section.

3. Feature Extraction

In this section, the feature extraction process will be discussed. Suppose the radar transmits wideband linear frequency modulation (LFM) signal and the echo signal can be written as:

E(tm,f)=s=1SAsrect(fB)rect(tmTa)exp[j4πc(fc+f)Rs(tm)]
(9)

where Rs(tm)=R0(tm)+rs(tm), R0(tm) is the radar-target range, rect(u)={1|u|120|u|>0,  B=γt˜ is the signal bandwidth, γ is the frequency modulation rate, t˜ is the fast time, fc is the carrier frequency, Ta is the observation time window, tm is the slow time, c is the light speed and As is the well-known reflectivity coefficient. Generally, the scattering coefficient is a complex function of frequency and radar look-angle. However, considering the narrow-angle measurement here, assume the backscattered intensity of each scatterer does not vary with the frequency and the view angle.

For the targets in high speed motion, the variation of radar-target range R0(tm) in the fast time, together with the residual video phase (RVP) terms, results in distorted and widened high resolution range profiles (HRRPs). Therefore, coherent pre-processing is carried out to eliminate their effects. After the coherent pre-processing method [25] is applied to the ISAR raw data:

E1(tm,f)=k=1SAsrect(fB)rect(tmTa)exp[j4πc(fc+f)rs(tm)]
(10)

The instantaneous range induced by spinning motion could then be achieved using FT on Equation (10). However, to realize the enhanced resolution, the 2D spectral estimation algorithm presented in [22] is used to extract the range and the range-rate data. Therefore, Equation (10) can be rewritten in discrete form as follows:

E1(m,n)=s=1SAsexp(j4π(f0+nΔf)(rs+r˙smΔt)/c)
(11)

where Δt is the pulse repetition time, rs and r˙s are the range and range rate of the scattering center Ps at time zero, respectively. n = 0,…, N − 1, N is the sample number of the LFM pulse in the fast-time domain, Δf = B/N, f0 = fcB/2, m = 0,…, Mhit− 1, Mhit is the total number of coherent process hits. One can notice that the cross product j4πmnΔfΔtr˙s/c of the exponent in Equation (11) keeps it from conforming to the state space signal structure sAsexp(jmαs+jnβs). To remove the cross product j4πmnΔfΔtr˙s/c, a new time-sampling interval Δt’ satisfying (f0 + nΔft = f0Δt’ is chosen, then the data matrix can be expressed as:

E2(m,n)=sAsexp(jmαs)exp(jnβs)
(12)

where αs=4πf0r˙sΔt/c, βs=4πrsΔf/c, As=Asexp(j4πf0rs/c). Now Δt is the original time sampling interval and the ranges of both m and n over the rows and columns of E2 are centered on zero. After transforming the measured data array E2 in a form of the output of two coupled eigenvalue problems, the estimation result α^s and β^s of the parameters can be straightforwardly determined. The amplitude estimation can be evaluated by a “least squares” procedure applied to (10), the α^s and β^s now being known quantities. Then, the pairs provide a direct extraction of the range rate and the range:

r˙^s=c4πf0Δtα^s
(13)

r^s=c4πΔfβ^s
(14)

4. Scattering Center Association

After the acquisition of the r^sr˙^s sequence from Equations (13) and (14), it is imperative that the sequence of r^s and r˙^s estimates are associated with the same physical scattering centres. Note that the measurement matrix is the input of the factorization method. Therefore, a good scattering centres association is very important for the whole procedure, especially in low signal-to-noise ratio (SNR) environments. In [27], the scattering centres association is realised by exploiting minimum Euclidean distance criterion based on the Kalman filter (KF) [23] together with motion features. For the results presented here, the scattering centres association is realised based on the Kalman filter (KF) [23] and together with the scattering centres’ amplitude and motion features.

For a discrete linear Gaussian kinematic and measurement system, the KF provides a closed- form solution to estimate the measurement system state vector and covariance from the measurement histories. Without the loss of generality, the measurement models at time tm can be expressed as:

x(m+1|m)=Fmx(m)+nxm
(15)

z(m+1|m)=Hmx(m)+nzm
(16)

where x(m) is the state vector, Fm is the transition matrix of the system, z(m) is the measurement vector, Hm is the measurement matrix, nxm and nzm are the process noise and measurement noise, respectively. The index m means the index of the time sequence, so the estimated state vector can be given as:

x^(m)=Fmx^(m1)+Gm(z(m)HmFmx^(m1))
(17)

where Gm is the Kalman gain.

Here, the KF is used to estimate the scattering centres’ motion state vector and covariance from measurement histories. In the absence of noise:

x^(m+1|m)=Fmx^(m)
(18)

From Equations (17) and (18), the prediction can be written as:

x^(m+1|m)=Fmx^(m|m1)+FmGm(z(m)HmFmx^(m1))
(19)

where:

FmGm=FmΨm/m1HmT[FmΨm/m1HmT+Γm]1
(20)

and the estimation error covariance matrix is:

Ψm/m1=FmΨm1FmT+Πm1
(21)

where F and Π are the covariance matrixes of the measurement noise and the process noise, respectively.

Considering the application here, it’s sufficient to describe the scattering center motion feature by the state vector which consists of the scattering centers range, range-rate and acceleration, i.e.,

x(m)=[rr˙r¨]T
(22)

We assume that the transition matrix is time-invariant, i.e., Fm=F=[1T0.5T201T001], Hm = H = [1 1 0], nzx~N(0,Γ), Γ=σz2I3, I3 is the identity matrix. nxm~ N(0,Π) with Π=[T4/4T3/20.5T2T3/2T2TT2/2T1]σx2, where T is the pulse repetition time, σx2 and σz2 are the variance of the Gaussian noise.

Now, the approach for scattering centers trajectory tracking is given in detail. Let the number of the slow time samples be Ma, so the range index m [set membership] [1, Ma]:

  • Step 1: Initialization. Let m = 1.
  • Step 2: For the m-th range and range-rate, we have the measurements {Aι,m,rι,m,r˙ι,m}, where ι is the scattering center track index, Aι,m is the amplitude estimate. Define the initial measurements as {Aι,0,rι,0,r˙ι,0}.

We assume that the transition matrix is time-invariant, i.e.: Fm=F=[1T0.5T201T001], Hm = H = [1 1 0], nzx~N(0,Γ), Γ=σz2I3, I3 is the identity matrix. nxm ~ N(0,Π) with Π=[T4/4T3/20.5T2T3/2T2TT2/2T1]σx2, where T is the pulse repetition time, σx2 and σz2 are the variance of the Gaussian noise.

Now, the approach for scattering centers trajectory tracking is given in detail. Let the number of the slow time samples be Ma, so the range index m [set membership] [1, Ma]:

  • Step 1: Initialization. Let m = 1.
  • Step 2: For the m-th range and range-rate, we have the measurements {Ai,m,ri,m,r˙i,m}, where i is the scattering center track index, Ai,m is the amplitude estimate. We define the initial measurements as {Ai,0,ri,0,r˙i,0}.
  • Step 3: Let m = m + 1. For the i-th scattering center track, search for the candidate scattering centers within the search window centered at ri,m1, then record the candidate set {Ai,m,ri,mi,r˙i,m}. Here, for a small observation interval, we assume the amplitude difference between the adjacent observation times is very small, so the optimal candidate for the i-th track can be determined by:
    {Aι,m,rι,m,r˙ι,m}=argmin(A^ι,m,r^ι,m,r˙^ι,m)Aι,mAι,m12
    (23)

    The other tracks are similarly associated with the scattering centers according to Equation (23). The acceleration of the scattering centers can be calculated as follows:

    r¨ι,0=1T(r˙ι,mr˙ι,m1)
    (24)

    Now, the initial state vector can be denoted as x(0)=[ri,0 r˙i,0 r¨i,0]T. Then, the next state {x^i,m+1/m} can be obtained according to Equations (19) to (21).

  • Step 4: Move to the next time index and let m=m+1. According to the minimum Euclidean distance criterion [23], the optimal candidate of the ιth scattering center track for the m-th range can be obtained by:
    {Aι,m,rι,m,r˙ι,m}=argmin(Aι,m,rι,m,r˙ι,m)ξι,mξι,m/m12
    (25)
    where {ξi,m} is the feature set of the candidate scattering centers within the gating area and ξi,m/m1 are the features of the predictions from the previous observables. ξι,m can be expressed as:
    ξι,m=[μ1Aι,m,μ2rι,m,μ3r˙ι,m]
    (26)
    where μ1, μ2 and μ3 are weight factors and μ1+μ2+μ3=1. According to Equations (23) and (26), Equation (25) can be rewritten as:
    {Aι,m,rι,m,r˙ι,m}=argmin(Aι,m,rι,m,r˙ι,m)(μ1Aι,mAι,m1Aι,m1)2+(μ2rι,mx^ι,m/m1(1)x^ι,m/m1(1))2+(μ3r˙ι,mx^ι,m/m1(2)x^ι,m/m1(2))2
    (27)

    Finally, by using Equation (30), the scattering centres association is accomplished.

  • Step 5: Repeat Step 4 to finish the whole measurements association.

It should be noted that the difference among scattering centres near the intersection points is so small that incorrect associations may appear. In this case, the multistep prediction in Steps 3 and 4 can be performed to increase the separability. Within the gating area, the combinations of features for all the candidate scattering centres are used to calculate the Euclidean distance. For avoiding the jump at the intersection points, the predicted amplitude will be replaced by the average amplitude of the existing associated tracks. And then, the minimum one in Equation (25) is the optimal association. Here, the three-step prediction is taken as an example. Suppose the feature set of all the candidates is {ξi,[k,k+2]}, then Equation (27) can be rewritten as:

{Aι,k,rι,k,r˙ι,k}=argmin(Aι,k,rι,k,r˙ι,k)ξι,[k,k+2]ξι,[k,k+2]2
(28)

where:

ξι,[k,k+2]=[μ1A¯ι,k1,μ2[x^ι,k/k1(1),x^ι,k+1/k(1),x^ι,k+2/k+1(1)],μ3[x^ι,k/k1(2),x^ι,k+1/k(2),x^ι,k+2/k+1(2)]]
(29)

where A¯i,k1 is the average amplitude of the existing associated i-th track. x^i,ς+1/ς (ς=k1,k,k+1) are the predictions.

5. Factorization-Based 3D Imaging and Scaling

Based on the above development, the general flow chart for processing a sequence of pulse echoes to form a 3D image is depicted in Figure 2. In this section, the 3D imaging procedure of the proposed method will be discussed in detail.

Figure 2
3D space spinning target imaging methodology.

5.1. 3D Imaging Based on Factorization Method

Because the spinning angular velocity is unknown, according to Equation (4), an initial coarse angular velocity Ω0 is chosen to scale the cross-range to obtain the measurement matrix. The rank of the measurement matrix W is highly rank-deficient [19]. In a real scenario, the rank of W is not exactly three, but approximately three. Therefore, orthogonal matrices through SVD decomposition can be achieved as:

W=[(U1)2M×3(U2)2M×(S3)][(Σ1)3×3000][(V1)3×S(V2)(S3)×S]
(30)

So we can factorize W into:

W=M^P^
(31)

where M^=U1(Σ1)1/2 and P^=(Σ1)1/2V1.

Since the decomposition is not unique, the following step is to find an invertible matrix Δ3×3 and then the true solutions M˜ and P˜ can be obtained as follows:

M˜=M^Δ
(32)

P˜=Δ1P^
(33)

According to the corresponding constraints in Equation (6), we obtain the system of 3M overdetermined equations, such that:

i^mTLi^m=1j^mTLj^m=1i^mTLj^m=0
(34)

where i^m and j^m are the mth row and the (M + m)-throw of M^, respectively. L = ΔΔT is the symmetric matrix:

L=[l1l2l3l2l4l5l3l5l6]
(35)

Equation (34) can be rewritten as:

Θl=χ
(36)

where Θ3M×6, l6×1, and χ3M×1 are defined by:

Θ=[gT(i1,i1)gT(iM,iM)gT(j1,j1)gT(jM,jM)gT(i1,j1)gT(iM,jM)],l=[l1l6],χ=[11}2M00}M]
(37)

And for two arbitrary vectors a = [a1 a2 a3] and b = [b1 b2 b3], we define:

gT(a,b)=[a1b1,a1b2+a2b1,a1b3+a3b1,a2b2,a2b3+a3b2,a3b3]
(38)

So l^ can be solved by the pseudo-inverse method, that is:

l^=(ΘTΘ)1ΘTχ
(39)

Then, the symmetric matrix L can be constructed using l^. Applying eigenvalue decomposition to L:

L=BΛBT
(40)

where B and Λ are the eigenvector matrix and the diagonal eigenvalue matrix, respectively. Then, the solutions for the invertible matrix can be obtained as:

Δ=BΛ1/2
(41)

Finally, M˜ and P˜ can be obtained by substituting Equation (41) into Equations (32) and (33), respectively.

5.2. 3D Image Scaling

Although the 3D image of the target can be obtained according to Equation (33), the image is modified by a scaling factor. This is because the cross-range is scaled by the coarse spinning angular velocity Ω0 and the scale ambiguities in cross-range will result in ambiguities in the spatial axes of the reconstructed target geometry. In this section, the image scaling process is introduced to remove these ambiguities.

Comparing to the real 3D geometry and projection matrix, it is worth noting that the estimated result M˜ and P˜ obtained via Equations (32) and (33) are the rotated projection matrix and 3D geometry matrix, respectively. Because the real 3D geometry and motion matrix is expressed in O-XYZ, suppose the rotated coordinate system is O-X’Y’Z’, as shown in Figure 3.

Figure 3
The new coordinate system O-X’Y’Z’ after 3D rotation.

The real 3D geometry estimate P˜r and motion matrix estimate M˜r can be obtained via the rotation operation, i.e.,

M˜r=M˜QT
(42)

P˜r=QP˜
(43)

where Q is the matrix with QTQ = QQT = I3×3, M˜r=[i˜1r,,i˜Mr;j˜1r,,j˜Mr]T, P˜r=[P˜1r,,P˜sr]. According to [24], the matrix Q can be represented by three continuous rotations around the three coordinate axes. We define the rotational angles around the three axes as [var phi]x, [var phi]y and [var phi]z, respectively, that is:

Q=RzRyRx
(44)

where:

Rz=[cosψzsinψz0sinψzcosψz0001], Ry=[cosψy0sinψy010sinψy0cosψy], Rx=[1000cosψxsinψx0sinψxcosψx]

According to Equation (6), j˜m is perpendicular to the OZ axis, as shown in Figure 3. After the rotation, j˜mr should be still perpendicular to the OZ axis, i.e., j˜mr(3)=0, so Q can be estimated via solving the optimization problem below:

(ψx,ψy,ψz)=argmin{[(jr)col,3]T(jr)col,3}
(45)

where jr=[j˜1r,,j˜Mr], and (jr)col,3 represents the third column of jr.

The optimization problem in Equation (45) is nonlinear and can be solved by exhaustive search. After Q is formed, M˜r and P˜r can be obtained via Equations (42) and (43). A new vector can be formulated as follows:

h(m)=j˜mr(2)+jj˜mr(1)=exp(j(Ωtm+φc))
(46)

where [var phi]c is the constant angle, which is an unimportant constant factor. The phase of Equation (46) can be written as:

h(m)=mΩT+φc
(47)

Due to the phase [var phi]c or the fast spinning angular velocity, during the observation time, phase ambiguity may occur. The ambiguity can be removed by judging the continuity of h(m)=mΩT+ϕc. Then by extracting the polynomial coefficients η^ of the phase, the estimate of the angular velocity can be obtained:

Ω^=|η^|T
(48)

Based on the estimated cross-range resolution, scaling of the cross-range coordinates in W is achieved. After that, the factorization method is performed to W again. The scaling and factorization steps will be carried out iteratively until the value of Equation (45) is smaller than a threshold ε0. The steps of the proposed method are described in Algorithm 1.

Algorithm 1: Processing steps of the proposed method
Input: Raw data
Pre-processing:
-Apply the pre-coherency processing [25] to the data
-Extract the range and range-rate using the spectral estimate method [22]
-Correlate the range and range-rate sequenceusing method developed in Section 4
Initialization: initialize k = 0, and choose an initial coarse angluar velocity Ω^0 and the error threshold ε0. Perform scaling to the sequential range-rate to form the projection matrix W0 according to Equations (3), (4) and (8)
Main iteration: increase k by 1 and perform the following procedure:
Step 1: Submit Ω^0 to Equation (4) to rescale the cross-range and obtain the scaled matrix Wk
Step 2: Calculate the motion matrix M^k and shape matrix P^k by substituting Wk into Equations (32) and (33)
Step 3: According to Equations (42)−(45) and the result in Step 2, calculate the rotational transform matrix Qk
Step 4: Estimate the angular velocity Ω^k by substituting Qk into Equations (46)−(48)
Step 5: If ([(jr)col,3]T(jr)col,3)<ε0, then Ω^k=Ω^true, and stop the iteration. Otherwise, repeat the main iteration.
Output: Estimate Ω^true. Consequently, the well-scaled shape matrix P^r can be generated using the optimal Ω^true.

6. Simulation Results

6.1. Simulation Results

In this section, the performance of the proposed method is verified using the point scattering centre model. As shown in Figure 4, we suppose the 3D target consists of eight scattering centres. The coordinates and the scattering coefficients are listed in Table 1.

Figure 4
The distribution of the scattering centers.
Table 1
Parameters of the scattering centers.

The radar signal centre frequency is 10 GHz and the bandwidth is 1 GHz, giving a range resolution of 0.15 m. The pulse repetition frequency is 100 Hz and the whole observation time is 1 s. Suppose the initial Euler angles are 30°, 20° and 50°, respectively. The elevation and azimuth angle of the RLOS are both 45°. The target spinning frequency is 2 Hz. The pulse time width is 0.3 ms and the dechirping signal sampling is rate 40 ZMHz. In this simulation, Gaussian noise is added to the simulated echo signal and the SNR is 10 dB.

Before using the spectrum analysis method [18] to extract the range and range-rate data, the coherent pre-processing [25] is performed on the echoes in advance. As a result, Figure 5a,b illustrate the sequential range and range-rate estimates, respectively. The colors in Figure 5 depict the normalized scattering coefficient of the scattering centers. Note that the scattering centers are clearly resolved in range and range-rate.

Figure 5
1D range and range rate estimates: (a) range estimate; (b) range rate estimate.

To obtain the measurement matrix W in Equation (8), the scattering center range and range-rate data set needs to be associated. Figure 6 presents the association results of the range and range-rate sequence using the method developed in Section 4. Here, the weight factors μ1, μ2, μ3 are set as 0.4, 0.4 and 0.2, respectively. The correlated scattering centers correspond to the same color.

Figure 6
1D range and range rate data association result: (a) 1D range data association result; (b) range rate data association result.

After finishing the extraction and tracking of the eight scattering centers, the initial coarse value of the angular velocity is set to 10 rad/s to form the projection matrix W0 according to Equations (3), (4) and (8).

The reconstructed 3D geometry using the McFadden’s method [20] is presented in Figure 7. The blue circles represent the reconstructed scattering centers. From Figure 7a, one can find there are transformation and scale ambiguities between the reconstructed target and the true one. This is because the cross-range in the projection matrix is not accurate. By using McFadden’s method [20], the ambiguity of the cross-range will result in the scale ambiguities of the reconstructed scattering centers position.

Figure 7
The reconstructed result using McFadden’s method [20]. (a) The distribution of the reconstructed 3D scattering centers before transformation; (b) The distribution of the 3D scattering centers after transformation; (c) the distribution on XY plane; ...

After the transformation, Figure 7b–e shows the distribution of the reconstructed scattering centers in 3D space, XY plane, XZ plane and YZ plane, respectively. From Figure 7b–e, one can find the apparent scale ambiguities in the 3D dimensions.

To remove the ambiguities shown in the Figure 7, the proposed scaling algorithm is implemented. Figure 8a shows the final reconstructed results using the proposed scaling algorithm and McFadden’s method, respectively. Here, the error threshold ε0 or the iteration procedure is set to 5e-2. After the transformation, the coordinates of the scattering centers are listed in Table 2. Figure 8b–e show the distribution of the reconstructed scattering centers in 3D space, XY plane, XZ plane and YZ plane. As shown in Figure 8b–e, the reconstructed target coincides well with the true one. Comparing to the reconstructed result via McFadden’s method, one can find the ambiguity has been removed by using the proposed scaling algorithm. The reconstruction result in Figure 8 proves the effectiveness of the proposed method.

Figure 8Figure 8
The reconstructed result using the proposed method. (a) The distribution of the reconstructed 3D scattering centers before transformation; (b) The distribution of the 3D scattering centers after transformation; (c) the distribution on XY plane; (d) the ...
Table 2
Coordinates of reconstructed scattering centers.

6.2. Performance Analysis

To evaluate the performance of the proposed algorithm quantitatively, the root mean square error (RMSE) of the recovered 3D target geometry is calculated in the terms of Euler distance error:

RMSE(P)=E[(P^P)T(P^P)]
(49)

The RMSE of the estimate of the spinning angular velocity is defined as:

RMSE(Ω)=E[|Ω^Ω|2]
(50)

where E[X] denotes the average of the X.

6.2.1. Effect of the SNR Level and Pulse Quantity

We note that the proposed algorithm performance is mainly affected by two factors: the noise level and the quantity of the pulses. Therefore, the experiments are designed to analyze the two factors. The first experiment is to analyze the effects of different SNR level on the algorithm performance. In this experiment, the pulse number is fixed at 100 and the SNR level varies from 0 to 20 dB. For each SNR, the experiment is carried out with 500 Monte-Carlo simulations, and the two RMSE curves against SNR are presented in Figure 9. As shown in Figure 9, with the increasing of SNR, the RMSEs of the two parameters decrease and low SNR level has an obvious impact on the 3D reconstruction and angular velocity estimation performance.

Figure 9
RMSE with respect to SNR. (a) 3D geometry; (b) spinning angular velocity.

In the second experiment, we test the algorithm performance by varying the quantity of the pulses. The SNR is set to 10 dB and the echo pulse number varies from 10 to 100 in steps of 10. The experiment is also carried out with 500 Monte-Carlo simulations for each pulse number, and the RMSE curves against pulse quantity is described in Figure 10. As shown in Figure 10, it can be found that when the pulse quantity is less than 30, both of the RMSEs of the two parameters decrease quickly with the increment of the pulse number. Therefore, more pulses will benefit the robustness of the proposed method. Meanwhile, when the pulse number is more than 50, there’s little difference in performance. This is because the target spinning angle is over 360° when the pulse number is over 50 and the information for reconstruction is abundant. From Figure 9 and Figure 10, we can draw the conclusion that the large number of utilized pulses and accurate extraction and tracking will bring good performance of the proposed method.

Figure 10
RMSE against pulse number. (a) 3D geometry; (b) spinning angular velocity.

6.2.2. Effect of Coarse Initial Angular Velocity

To analyze the impact of the different initial angular velocity values on the 3D imaging and angular velocity estimation, here another experiment is designed. In this experiment, the coarse angular velocity used to start the reconstruction procedure varies from 10.0 rad/s to 13.5 rad/s. The number of the pulse used for imaging is 100. The SNR level varies from 0 to 20 dB with a step of 5 dB. For each SNR level and initial angular velocity, 500 Monte-Carlo simulations are carried out. Figure 11 presents the RMSEs calculation result of the 3D reconstruction and angular velocity estimate. The largest RMSEs of 3D imaging and angular velocity estimation are 0.34 m and 0.16 rad/s, respectively which are acceptable. From Figure 9a,b, one can find that RMSEs decrease with the initial coarse value getting closer to the true value. And The RMSEs of the two estimates are minimized when the initial angular velocity is equal to the true value. Therefore, the precise initial angular velocity will benefit the target 3D imaging.

Figure 11
RMSEs of 3D reconstruction and angular velocity estimation. (a) RMSE of the 3D reconstruction against different initial angular velocity and SNR; (b) RMSE of the angular velocity estimate against different initial angular velocity and SNR.

7. Conclusions

In this paper, a 3D imaging and scaling algorithm for rapidly spinning target based on factorization method is proposed. Due to the lack of freedom, the recovered 3D imaging of the spinning target via traditional factorization method is a rotated and scaled version of the true one. The proposed method provides a new solution for 3D imaging of spinning targets and removes the scale ambiguities in the recovered 3D image. Simulation results show the effectiveness and robustness of the proposed algorithm. In the future, we will test the algorithm with real measured data.

Acknowledgments

This study was co-supported by the National Natural Science Foundation of China (Nos. 61671035, 61501012, and 61501011).

Author Contributions

Author Contributions

In this study, Yanxian Bi and Jun Wang conceived and design this novel algorithm and simulation experiment; Shaoming Wei and Shiyi Mao provide the advices in the details; Yanxian Bi wrote this paper.

Conflicts of Interest

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Xue R.B., Feng Z., Zheng B. High-resolution three-dimensional imaging of space targets in micromotion. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 2015;8:3428–3440.
2. Wei Q., Martorella M., Jian X.Z., Hong Z.Z., Qiang F. Three-dimensional inverse synthetic aperture radar imaging based on compressive sensing. IET Radar Sonar Navig. 2015;9:411–420.
3. Bai X., Bao Z. High-resolution 3D imaging of procession cone-shaped targets. IEEE Trans. Antennas Propag. 2014;3:1680–1689.
4. Zhou X.P., Wei G.H., Wu S.L., Wang D.W. Three-dimensional ISAR imaging method for high-speed targets in short-range using impulse radar based on SIMO array. Sensors. 2016;16:364 doi: 10.3390/s16030364. [PMC free article] [PubMed] [Cross Ref]
5. Lei L., Feng Z., Xue R.B., Ming L.T., Zi J.Z. Joint cross-range scaling and 3D geometry reconstruction of ISAR targets based on factorization method. IEEE Trans. Image Process. 2016;25:1740–1750. [PubMed]
6. Xu X., Narayanan R.M. Three-dimensional interferometric ISAR imaging for target scattering diagnosis and modelling. IEEE Trans. Image Process. 2001;10:1094–1102. [PubMed]
7. Suwa K., Wakayama T., Iwamoto M. Three-dimensional target geometry and target motion estimation method using multistatic ISAR movies and its performance. IEEE Trans. Geosci. Remote Sens. 2011;49:2361–2373. doi: 10.1109/TGRS.2010.2095423. [Cross Ref]
8. Martorella M., Stagliano D., Salvetti F., Battisti N. 3D interferometric ISAR imaging of noncooperative targets. IEEE Trans. Aerosp. Electron. Syst. 2014;50:3102–3114. doi: 10.1109/TAES.2014.130210. [Cross Ref]
9. Qian Q.C., Gang X., Lei Z., Meng D.X., Zheng B. Three-dimensional interferometric inverse synthetic aperture radar imaging with limited pulses by exploiting joint sparsity. IET Radar Sonar Navig. 2015;9:692–701.
10. Wang Y., Abdelkader A.C., Zhao B., Wang J. ISAR imaging of maneuvering targets based on the modified discrete polynomial-phase transform. Sensors. 2015;15:22401–22418. doi: 10.3390/s150922401. [PMC free article] [PubMed] [Cross Ref]
11. Stuff M.A. Three-dimensional analysis of moving target radar signals: Methods and implications for ATR and feature-aided tracking. Proc. SPIE. 1999;3721:485–496.
12. Stuff M., Sanchez P., Biancalana M. Extraction of three-dimensional motion and geometric invariants from range dependent signals. Multidimen. Syst. Sig. Process. 2003;14:161–181. doi: 10.1023/A:1022229326812. [Cross Ref]
13. Ferrara M., Arnolad G., Stuff M. Shape and Motion Reconstruction from 3D-to-1D Orthographically Projected Data via Object-Image Relations. IEEE Trans. Pattern Anal. Mach. Intell. 2009;31:1906–1912. doi: 10.1109/TPAMI.2008.294. [PubMed] [Cross Ref]
14. Wu S., Hong L. Modelling 3D rigid-body object motion and structure estimation with HRR/GMTI measurements. IET Control Theory Appl. 2007;4:1023–1032. doi: 10.1049/iet-cta:20060034. [Cross Ref]
15. Zhou J.X., Zhao H.Z., Shi Z.G. Performance analysis of 1D scattering center extraction from wideband radar measurements; Proceedings of the IEEE International Conference on Acoustics Speech and Signal Processing Proceedings; Toulouse, France. 14–19 May 2006; pp. 1096–1099.
16. Wang Q., Xing M.D., Lu G. High-resolution three-dimensional radar imaging for rapidly spinning targets. IEEE Trans. Geosci. Remote Sens. 2008;46:22–30. doi: 10.1109/TGRS.2007.909086. [Cross Ref]
17. Mayhan J.T., Burrows K.M., Piou J.E. High resolution 3D “Snapshot” ISAR imaging and feature extraction. IEEE Trans. Aerosp. Electron. Syst. 2001;2:630–642. doi: 10.1109/7.937474. [Cross Ref]
18. Tomasi C., Kanade T. Shape and motion from image streams under orthography: A factorization method. Int. J. Comput. Vis. 1992;9:137–154. doi: 10.1007/BF00129684. [Cross Ref]
19. Morita T., Kanade T. A sequential factorization method for recovering shape and motion from image streams. IEEE Trans. Pattern Anal. Mach. Intell. 1997;19:858–867. doi: 10.1109/34.608289. [Cross Ref]
20. McFadden F.E. Three-dimensional reconstruction from ISAR sequences. Proc. SPIE. 2002;4744:58–67.
21. Li G., Liu Y., Wu S. Three-dimensional reconstruction using ISAR sequences. Proc. SPIE. 2013;8919:891–908.
22. Michael L.B. Two-dimensional ESPRIT with tracking for radar imaging and feature extraction. IEEE Trans. Antennas Propag. 2004;52:1680–1689.
23. Blackman S., Popoli R. Design and Analysis of Modern Tracking Systems. Artech House; Norwood, MA, USA: 1999.
24. Chen V.C., Li F.Y., Ho S.-S., Wechsler H. Micro-doppler effect in radar phenomenon, model and simulation study. IEEE Trans. Aerosp. Electron. Syst. 2006;42:2–21. doi: 10.1109/TAES.2006.1603402. [Cross Ref]
25. Xing M., Wu R., Bao Z. High resolution ISAR imaging of high speed moving targets. IET Radar Sonar Navig. 2005;152:58–67. doi: 10.1049/ip-rsn:20045084. [Cross Ref]
26. Diebel J. Representing Attitude: Euler Angles, Unit Quaternions, and Rotation Vectors. [(accessed on 3 February 2017)]. Available online: https://www.astro.rug.nl/software/kapteyn/_downloads/attitude.pdf.
27. Bai X., Zhou F., Bao Z. High-resolution radar imaging of space targets based on HRRP series. IEEE Trans. Geosci. Remote Sens. 2014;52:2369–2381. doi: 10.1109/TGRS.2013.2260342. [Cross Ref]

Articles from Sensors (Basel, Switzerland) are provided here courtesy of Multidisciplinary Digital Publishing Institute (MDPI)