Parameter estimation for SAR micromotion target based on sparse signal representation

In this article, we address the parameter estimation of micromotion targets in synthetic aperture radar (SAR), where scattering parameters and micromotion parameters of targets are coupled resulting in a nonlinear parameter estimation problem. The conventional methods address this nonlinear problem by matched filter, which are computationally expensive and of lower resolutions. In contrast, we address this problem by linearizing the forward model as a linear combination of elements of an over-complete dictionary. The essential idea of sparse signal representation models comes from the fact that SAR micromotion targets are sparsely distributed in the observation scene. Accordingly, we propose to jointly estimate the target micromotion and scattering parameters via a Bayesian approach with sparsity-inducing priors. In addition, we present a variational approximation framework for Bayesian computation. Numerical simulations demonstrate the proposed sparsity-inducing reconstruction method achieves higher resolution and better performance with smaller measures compared to conventional methods.


Introduction
Target micromotion and micro-doppler are attracting an increasingly great interest from the synthetic aperture radar (SAR) community since they can provide additional and favorable information for understanding SAR images. Micromotion is mainly embodied by rotation and vibration, and typical SAR micromotion targets include ground/ship-borne search antennas for air traffic control/surveillance [1], rotor blades of hovering helicopters, and vibrating vehicles as well as their tires/ engines. Micromotion parameters, such as the rotating frequency and radius, record targets' attributed information. Thus their estimation is very important for micromotion compensation and refocusing in SAR imagery, and the estimated results can also be directly used as signatures for target recognition. However, it's a huge challenge for micromotion parameter estimation in SAR, since (1) micromotion target signals are hard to be separated from stationary-clutter ones, (2) they are also distributed over multiple range cells (especially for large rotating radii), i.e., range cell migration (RCM) occurs, which is disadvantageous for target energy integration. Either, it's not practical to estimate them in the SAR gray image domain because of defocusing, ghost images [2] and other energy-spread image characteristics induced by target micromotion [3].
A few algorithms have been proposed for the estimation of SAR micromotion targets [1,[4][5][6]. All of them manipulate a single range cell and take micromotion target azimuthal echoes as sinusoidal frequency modulated (SFM) signals. The cyclic spectral density [4], a time-frequency method [6], and the adaptive optimal kernel one [5], have been used to estimate the vibrating frequency of simulated or real SAR targets. Then in [1], the wavelet or chirplet decomposition is used to separate the signal of a rotating radar dish from that of stationary clutter and then auto correlation is utilized to get its rotating frequency. All these methods, however, haven't addressed the aforementioned two key problems ever-present in SAR, i.e., clutter and RCM, which hinders their application in reality. In effect, unlike uniformly moving targets, RCM correction is very difficult for micromotion ones due to their sinusoidal range history [7].
Matched filter is commonly used for motion or micromotion target imaging [8,9]. It performs the reconstruction at every pixel for every possible velocity of the motion, resulting in a huge space-velocity cube [8]. Worse still, for the fact that each slice of the velocity is estimated independently, it brings in ambiguous results. To improve this, an adaptive matched filtering method, called filtered back projection, was proposed by Cheney [9]. However, all these methods yield high computational cost and ambiguity unavoidably caused by independent estimation. Recently, sparse signal representation and compressive sensing (CS) have become a standing interest for SAR imaging [10][11][12][13]. A joint spatial reflectivity signal inversion method based on an over-complete dictionary of target velocities was applied to SAR moving targets imaging [10]. However, large scaled matrix computation is still treated as an open problem.
Hence we propose to obtain micromotion parameters from the viewpoint of scattering center estimation, which circumvents the tough issues mentioned above via target model priors. The scattering center model, however, must herein consider target micromotion, and thus more parameters, besides target position, and higher dimensions are involved which create adverse effects on fast and global optimization. Fortunately we observe finer target s-parsity due to an increase of the parameter space dimension. Therefore we will exploit target priors and estimate the model based on sparse signal reconstruction. We recast the micro-motion target imaging problem as a problem of signal representation in an over-complete dictionary. To enforce sparsity, we consider two Baysian prior models: generalized Gaussian and Student-t. Then we examine the expression of posterior laws, either the maximum A poseriori (MAP) estimator or the posterior means using the variational Bayes approximation (VBA) [14]. Compared to conventional methods, besides overcoming two difficulties aforementioned, the advantages of our method include: (1) putting the micromotion target imaging and parameter estimation into a unified Bayesian parameter estimation framework, which could also handle the hyperparameter estimation; (2) breaking through the classic Relay resolution's limit, providing the capability of super-resolution; (3) being capable of estimating micromotion parameters from limited observations; (4) being robust to noise.
The rest of the article is organized as follows. Section 2 presents the SAR signal model of micromotion targets. In Section 3 we review the different sparse modeling and optimization criteria. In particular, l 1 regularization approach conducts us to the Bayesian approach which is developed in Section 4. We provide two priors as generalized Gaussian priors and Student-t priors, which enforce the sparsity. Section 5 provides simulation results and performance analysis. Finally, Section 6 summarizes our conclusions.
2 Wavenumber-domain signal model of SAR micromotion targets As illustrated in Figure 1, the radar moves at velocity V a . Then for slowtime t it moves to We could see that θ has the similar meaning as slowtime t. Considering an arbitrarily moving target, let vector ϑ represent the target micromotion parameters, such as the initial position (x, y), velocity, rotation frequency etc. Suppose the target moves to (x ϑ,θ , y ϑ,θ ) when radar is located at y'. f (ϑ) is the scattering coefficient. Thus the distance model of the target is Spotlight SAR echo of the target could be represented in the wavenumber domain as where P (K) is the Fourier transform (FT) of the transmitted signal. Then the total echoes of all targets are After range compression and motion compensation, the first two terms of s (⋅) in Equation (3) disappear, and then the target signal model becomes When the target experiences micromotion, e.g., rotation or vibration, we have where micromotion parameters compose a parameter vector and (x,y) is the position of the micromotion center, r is the micromotion amplitude, i.e., rotating radius or vibrating amplitude, f m is the micromotion frequency, and 0 is the initial micromotion phase. Substituting Equations (6) and (7) into (5) leads to where We can clearly see that, Equation (10) has an additional exponential component representing target micromotion, compared with the stationary scattering center model [5].
We now try to discretize Equation (10). Without loss of generality, suppose there are I rotated targets. Then for the ith one, let f i denotes the scatter coefficient, ϑ i denotes the micromotion parameter, both of which are unknown. The model of Equation (9) could be discretized as where noise has been added via i (K, θ). Note K and θ can also be discretized into M and N values respectively, and therefore Equation (11) can be expressed in a matrix form as where is a vector of size M N representing the data, is a vector of size M N representing the errors (modeling and measurement), is a matrix of dimensions M N × N x N y PQJ representing the forward modeling matrix system and is a vector of size N x N y PQJ of parameters representing targets in the scene. In this expression A(x n x , y n y , r p , f q , ϕ 0 j ) is the coefficient at position (x n x , y n y ) with micromotion frequency f q , micromotion range r p and initial micromotion phase ϕ 0 j To this end, the problem of scattering and micromotion parameter estimation can be reformulated as a linear inversion problem subject to sparsity constraints.

Sparse signal representation and deterministic optimization
The main idea behind sparse signal representation is, to find the most compact representation of a signal as a linear combination of a few elements (or atoms), in an over-complete dictionary [15][16][17][18]. Compared with the conventional orthogonal transform representation, this most parsimonious representation of a signal over a redundant collection of generated basis offers efficient capability of signal modeling. Finding such a sparse representation of a signal involves solving an optimization problem. Mathematically, it can be formulated as follows. For Equation (2), assume g = H f in absence of noise where g ℂ M × 1 is a vector of data, H ℂ M × N a matrix whose elements can be considered as an overcomplete dictionary as its columns and f ℂ N × 1 the corresponding linear coefficients. In particular, M ≪ N leads the null space of F to be non-empty such that there are many different possibilities to represent g with the elements in H. The problem of sparse representation is then to find the coefficients f with the most few nonzero elements, i.e., ||f|| 0 is minimized while g = H f. Formally, where ||f|| 0 is the l 0 norm which is the cardinality of f. However, the combinatorial optimization problem Equation (17) is NP-hard and intractable. A large body of approximation methods are proposed to address this optimization problem, such as greedy pursuit [19] based methods like matching pursuit [20], or convex-relaxation [21] based methods that replace the l 0 with the l 1 norm, Candes et al. [22] show that for K-sparsity signal that only has K non-zero element in f, the reconstruction of f with M ≥ O(K log (N/K)) [15] measures can be achieved with high probability by l 1 norm minimization. Moreover, to efficiently reconstruct f, the mapping matrix H should satisfy the restricted isometry property (RIP) [23] which requires that This RIP of H is connected to the mutual coherence between the atoms of the dictionary which is defined as where the a i is the ith column of H. Large mutual coherence indicates that there are two atoms that are closely related will degrade the reconstruction algorithm. Hence, the dictionary is required to have low coherence so that the submatrix H with K atoms is nearly orthogonal [18].
If the observation g is noisy, the problem of the sparse representation for a noisy signal can be formulated as where δ is a noise allowance. Equivalently, the Equation (21) can be reformulated to minimize the following objective function where l > 0 is the regularization parameter that balances the trade-off between the reconstruction error and the sparsity of f. The formulation Equation (22) can also be interpreted as the MAP estimation in the Bayesian philosophy as we will see in the next section.
To this end, the micromotion parameter estimation is now cast as the sparse reconstruction of f associated with the parameter hypothesis at the position of nonzero elements of f. There are a large number of methods to solve the Equations (21) or (22), such as the method of compressive sampling matching pursuit (CoSaMP) presented in [24] which has been widely used for its simplification and effectiveness. Here, we will compare our proposed method with this method.

Bayesian approach to sparse reconstruction
Even if the sparse representation has originally been introduced as an optimization problem such as Equations (17), (18), (21), or (22), it can also be presented as a Bayesian MAP estimation problem [25,26]: where To understand this, firstly let us assume the error in Equation (12) is centered, Gaussian and white: ε ∼ N (ε|0, v ε I). It brings us to the expression of the likelihood: Secondly, choose the separable double exponential probability density [27] as the prior of f: it is then easy to see that the MAP estimation with this prior becomeŝ which can be compared to Equation (22). The prior information that the targets are sparsely distributed in the observation scene can be modeled by the two following probability density functions (PDF) [14]: • Generalized Gaussian priors: which give the double exponential for b = 1 and Gaussian for b = 2 and are also more useful for sparse representation with 0 <b < 1. With these priors, the MAP estimate can be computed by optimizing the following criterion: which can be done with any gradient based algorithm when 1 <b ≤ 2. There also exist appropriate algorithms for b = 1 and 0 <b < 1. In this article, we used a gradient based algorithm.
• Student-t priors: where These priors are interesting due to its link to l 1 regularization and secondly due to the mixture of Gaussian representation of the Student-t probability density: which gives the possibility of proposing a hierarchical model via the positive hidden variables τ j : Using this hierarchical model, we can write the joint prior of f and τ we obtain: where which is summarized as follows: Joint optimization of this criterion, alternatively with respect to f (with fixed τ)f and with respect to τ (with fixed f) results in the following iterative algorithm: Note that τ j is inverse of a variance and we have 1/τ j = f 2 j + β/α . We can interpret this as an iterative quadratic regularization inversion followed by the estimates of variances τ j which are used in the next iteration to define the variance matrix D(τ). This algorithm is simple to implement. However, we are not sure about its convergency. To obtain a better solution and at the same time to be able to estimate the variance of the noise, we propose to use the VBA [28][29][30] which consists in approximating the joint posterior by a separable one and then using it to do the inference.
Here we summarize this approach: • Model for the noise: • Model for the sparse signal: • Joint posterior: • VBA: p(f,τ,τ |g) is approximated by and the expressions of the needed expectations are: This algorithm can be summarized as follows: • Initialization:τ ε = 0.1,Ṽ = diag[τ j /τ ε ] withτj = 1 • Iterations: computeα ε ,β ε and soτ ε =α ε /β ε computeα j ,β j and soτ j =α j /β j The only difficult and costly part is the estimation of andμ . Due to the fact that we only needμ and˜ jj , we propose the following approximation:μ is computed through the optimization of J(f ) = τ ε g − Hf 2 + 1 2 j τ j f 2 j with respect to f and˜ jj which is the variance of f j is approximated by the empirical variance of f j during the iterations of the optimization algorithm.
This is the method we implemented, tested and compared to other classical methods.

Numerical experiments
In this section, we conduct several numerical experiments to demonstrate our method based on the sparse signal representation. The imaging geometry is shown in Figure 1. The range R 0 from the original to the center of the target is 10 km, and the velocity of the platform V a is 200 m/s. The central frequency f c is 10 GHz with bandwidth B = 400 MHz associated with the Rayleigh resolution along the range direction 0.375 m, and the angular extent of azimuth is 10°with cross-range resolution 0.0861 m.
Based on the compressive sensing principle, the targets can be recovered with a smaller randomly sampled measures. Figure 2 shows that sampling pattern in the wavenumber domain, the uniform sampling in Figure 2a and the randomly sampling in Figure 2b. Two targets are located at (0,0), (5,1), respectively. With the randomly sampled measures, Figure 3 compares the reconstruction results between the traditional method of fast FT (FFT), the CoSaMP and the Bayesian sparse method when no micromotion is present. It is shown that the CoSaMP and the proposed Bayesian method come out with clearer images and are capable to recover the true position of scatters, compared with the traditional method of FFT, even with smaller randomly measures.
When targets experience micromotion, the initial phases are assumed to be both zeros. The micromotion frequencies are 0.5, 1Hz, respectively, and the micromotion range is 1 and 0.5, respectively. In Figure 4b, the range profile appears clearly in the micromotion pattern compared with Figure 3b. The presentation of micromotion blurs the reconstruction images without motion compensation as shown in Figure 4a image in Figure 4d recovering the true parameters (x,ŷ,r,f m ,φ 0 ) = (0, 0, 1, 0.5, 0) , and (5,1,0.5,1,0) for the two scatter points, respectively. Figure 4c illustrates the reconstruction result via the CoSaMP method. We then set the micromotion range 0.5 and initial phase 0 for both targets but the micromotion frequencies are 0.5 and 1 Hz, respectively. We adopt the matched filtering in the 3D range-Azimuth-micromotion frequency space by scanning a large number of possible scatterer positions and micromotion frequencies, resulting in a large space-micromotion frequency cube. Figure  5a shows the 3D data cube. Figure 5b,c illustrate the two slices after matched filtering at micromotion frequencies f m = 1 Hz and f m = 0.5 Hz, respectively. It is computationally expensive and not well focused being of low resolution. In addition, it is rather difficult to perform RCM such that the position cannot be estimated accurately. In contrast, our method can overcome these drawbacks of traditional methods and yield a more precise estimate. Figure 6 shows that our proposed method resolves the two very closely spaced micromo-tion targets localized at positions of (0, 0) and (0.25, 0.25), respectively. The reconstruction image by FFT is illustrated in Figure 6a and the corresponding range profile in Figure 6b. It shows that the range profiles of the two targets are overlapped so that the two targets cannot be discerned.  prove the super-resolution capability of the proposed method. Figure 7 depicts the estimation root mean square (RMS) error varies with SNR which demonstrates our method can recover the targets signature parameters accurately. It can be observed that the RMS decreases sharply as the SNR increases and arrives at high precision estimations after 0dB, indicating the robustness of our method to loss and noise of measurement.

Conclusions
In this article, we proposed a sparsity-inducing method to estimate the scattering and mi-cromotion parameters of SAR targets jointly and further formatted it in the Bayesian framework. It was done by formulating the original nonlinear problem as a sparse representation problem over an over-complete dictionary. In addition, an efficient computation algorithm as VBA estimator was applied to the hierarchical Bayesian models. The proposed method can exactly recover the scattering and micromotion parameters of targets, even for near spacing targets, achieving good performance, as demonstrated by the simulation experiments.

CoSaMP algorithm
The basic idea of CoSaMP [24] algorithm is that: for Ssparse signal f with S non zero elements, the z = H ⋆ H f can serve as a proxy for the signal where H ⋆ is the Hermitian transpose of H, since the energy in each set of S components of z approximates the energy in the corresponding components of f. In particular, the largest S entries of the proxy z point toward the largest S entries of the signal f. a The basic steps are (appendix table): (1) Identification: Compute z H ⋆ y to find a proxy of the residual from the current samples and locate the largest components Ω = supp(z 2S ) of the proxy z; (2) Support merger: The set of newly identified components Ω is united with the set of the components that appear in the previous approximation supp(f k-1 ), i.e., T = Ω ∪ supp(f k-1 ); (3) Estimation: Solve a least-square problem to approximate the target signal on the merged set T of component, b| T ← H † T g (4) Pruning: Obtain a new approximation by retaining only the largest entries in this least-square signal approximation, f k b S ; (5) Sample update: Finally, update the residual g -H f k .