A Tensor Decomposition Based Multiway Structured Sparse SAR Imaging Algorithm with Kronecker Constraint

Yu-Fei Gao 1, Xun-Chao Cong 1, Yue Yang 1, Qun Wan 1 and Guan Gui 2,* 1 School of Electronic Engineering, University of Electronic Science and Technology of China, No. 2006, Xiyuan Ave., West Hi-Tech Zone, Chengdu 611731, China; yufeiee@gmail.com (Y.-F.G.); XunchaoCong@std.uestc.edu.cn (X.-C.C.); 201611020223@std.uestc.edu.cn (Y.Y.); wanqun@uestc.edu.cn (Q.W.) 2 College of Telecommunication and Information Engineering, Nanjing University of Posts and Telecommunications, 66 Xinmofan Road, Nanjing 210003, China * Correspondence: guiguan@njupt.edu.cn


Introduction
Synthetic aperture radar (SAR) is a kind of radar system with the features of all-weather, multi-polarization and penetrability, which played a key role in remote sensing, detection and mapping fields in the past decades [1][2][3].There are four common observation modes for SAR: Stripmap, Scan, Inverse and Spotlight [4].Because the high resolution is conducive to capturing the structure information of the SAR targets, this paper focuses on Spotlight SAR.It is a beamforming-based observation technology, for which the radio beam will always point to the observation scene.The SAR image reconstruction is to acquire the scatter coefficients and to recover spatial distribution of the scatterers from the reflected waves.To obtain a better imaging resolution, a larger aperture will be synthesised by the coherent processing for the motions of targets relative to the antenna, as the pulse compression is usually applied on the range direction.However, in practical applications, traditional SAR imaging methods encounter many problems.Firstly, the receiver needs a higher sampling rate for the wideband waveform because of the Nyquist sampling limit, which leads to a considerable cost for data storing and processing.Secondly, the conventional techniques often suffer from the sidelobe disturbance which may lead to the performance deterioration of reconstruction, such as fuzzy in the SAR target image, etc.
Compressed sensing (CS) [5][6][7], as a novel theoretical framework appearing in literature recently, can obtain necessary information via sub-Nyquist sampling rate for sparse signal, which can reduce storage requirement and simplify the design of hardware system.That is why CS has attracted greater attention in SAR imaging.In 2001, Cetin et al. introduced the sparse constraint into spotlight SAR imaging, and proposed a method to improve the target feature [8].CS was first applied in super-resolution radar imaging by Baraniuk et al. in 2007, and a new type of radar system was proposed based on sparse reconstruction [9].In 2009, Gurbuz et al. analyzed the application of sparse reconstruction in ground-penetrating radar imaging, and proposed a sparse imaging algorithm based on l 1 -norm minimization [10].In the next year, a SAR imaging method based on sparse reconstruction theory was introduced by Potter et al. [11], which can achieve a higher reconstruction resolution than the traditional imaging methods.In 2011, the effect of parameter selection of SAR imaging based on regularized sparse reconstruction was studied by Batu et al. [12].In order to accelerate the speed of SAR imaging algorithm, a sparse regularization technique based on iterative threshold method was proposed by Fang et al. in 2014 [13].For moving targets, Zhang et al. proposed an efficient SAR imaging method in 2015 [14] combining fractional Fourier transform technology with CS framework.However, these studies dealt only with the spatial sparsity of the observation without exploiting the structural characteristics of scatterers in the target scene.In practice, the imaging results show that the scatterers tend to be aggregated together as clumps.The sparse reconstruction schemes have ability to accomplish SAR imaging based on object-level information of scatterers and possess the robustness to noise [15].In 2016, Zhao proposed a CS-based SAR imaging algorithm with structured sparsity constraint, in which a hierarchical Bayesian learning framework was applied [16].Since hierarchical model is designed, the requirement for model matching of such kind of algorithm is strict, thus the imaging quality may deteriorate once observation data is mismatched with the model.
To solve the above problems, we introduce a tensor decomposition algorithm in this paper both to reap the model robustness and to lower resource requirement.The real-world signal is usually multidimensional and structured.In order to deal with such higher-order data processing problem, tensor decomposition is studied widely.For instance, tensor modeling is utilized in chemometrics to analyse spectroscopy [17], and is also applied to acquire two-dimensional direction of arrival (DOA) estimation by taking advantage of the Vandermonde structure in array manifold [18,19].For SAR imaging problems, tensor modeling is capable of exploiting the structured sparse distribution of the scatterers to gain more robustness.Because of the curse of dimensionality in higher-order data processing [20][21][22], making use of the sparsity is becoming a research hotspot in recent years.In 2010, Lim and Comon applied tensor decomposition to propose a sparse reconstruction method, and derived the uniqueness condition of tensor decomposition [23].Duarte et al. first introduced the Kronecker constraint into the construction of dictionary in 2011, and proposed the Kronecker-CS model [24].In low rank tensor signal processing, Sidiropoulos et al. introduced the CS theorem into a two-step sparse signal reconstruction algorithm in 2012 [25].The following year, Caiafa and Cichocki proposed two types of matching pursuit algorithms based on tensor model incorporating with CS method [26].Although there has been little public research about combining sparse representation with tensor decomposition in SAR imaging field, potential advantages of this method with robustness and efficiency inspire us to apply it in solving the open problems.The application of higher-order CS in image processing, such as hyperspectral imaging [24,27] provides the inspiration of this work, but this is still in the stage of preliminary study.We formulate the SAR echo signal as a Tucker model based on its Kronecker information and then construct a sparse reconstruction algorithm by utilizing the structured sparsity of target scene.Both the theoretical analysis and numerical experiment verifications reveal that the proposed algorithm is superior to the existing methods, such as spatially variant apodization (SVA) [28,29] and matching pursuit algorithms [30][31][32].
The following is the organization of other sections in this paper: Section 2 describes the received signal model of SAR imaging and its tensor modeling; Section 3 constructs the algorithm based on Kronecker constraint and multiway structured sparsity of the signal model; Section 4 provides several groups of numerical experiments to verify the proposed method, and compare with the existing methods in performance; Section 5 concludes this paper.

Received Signal Model
For convenience of statement hereafter, several notions of this paper are explained as follows.The scalar is indicated by italic letter, e.g., s; the vector and matrix are indicated by bold italic letters, e.g., v, M; the tensor is indicated by calligraphic capital letter, e.g., T ; the index of a set is indicated by script letter, e.g., I .
Here we consider the point-scattering models [3].The spotlight SAR geometry is illustrated in Figure 1.The target scene is represented by O(x, y, z) coordinate system.The azimuth and pitch angles of antenna relating to the center of the scene are φ and θ.The distance between target and antenna is R (x, y, z) = R c − x cos (θ) cos (φ) + y cos (θ) sin (φ) + z sin (θ), where R c denotes the distance between the antenna and the center of target area.The time delay for the target point (x, y, z) in SAR echo signal is τ (x, y, z) = 2R (x, y, z) /c, where c is the speed of light.After matched filtering, the reflection signal is where g (x, y, z) is the scattering function of scatterer placed at (x, y, z); s (t) is the autocorrelation function of transmitted signal; e (t) is the additive noise.In the following discussion the noise is neglected for notational simplicity.Then the phase history can be obtained through Fourier transform: where can be ignored as it is a constant for all scatterers, and the phase factor is H x,y,z ( f , θ, φ) = j 4π c (x f cos θ cos φ + y f cos θ sin φ + z f sin θ).For decoupling between the parameters in the phase factor, an interpolation approach is applied and the polarization coordinates is transformed to a Cartesian coordinates [33].Assume there are M 1 , M 2 , M 3 grid points along the directions of the axis of x, y, z, respectively; and the number of grid points for f , θ, φ are P 1 , P 2 , P 3 , respectively.Then for the scatterer at x j , y k , z l , the steering vector is a j,k,l = b j,k,l (1, 1, 1) , . . . ,b j,k,l (P 1 , P 2 , 1) , . . . ,b j,k,l (1, 1, P 3 ) , . . . ,b j,k,l (P 1 , P 2 , P 3 ) T , where then the received signal model ( 2) can be rewritten as where Note that each element in the steering matrix can be rewritten as where ⊗ denotes the Kronecker product.
, A 2 and A 3 have the similar formats.It is demonstrated from (5) that the steering matrix possesses the Kronecker constraint.Therefore, (4) can be written as

Kronecker Dictionary
For a vector signal y, we consider that a K-sparse representation is existed on a dictionary A, if the following relationship is satisfied: where s 0 is the l 0 -norm of s, indicating the number of nonzero elements in s; K is defined as sparsity.
Usually, K M. Because the number of rows is less than the columns, this is a underdetermined system which means there is no unique solution of s.However, a sparse solution can be obtained if the signal has the feature of sparsity and the following condition is satisfied [30]: where µ (A) = max i =j a T i a j , indicating the coherence coefficient of dictionary A.

If the dictionary has Kronecker constraint, i.e., A = A
where It is shown from ( 9) that for a Kronecker dictionary, the overall coherence depends merely on the factor matrix with the greatest coherence [34].
The previous studies [35,36] indicated that, for the signal recovery problem with the same sparsity level, the reconstruction procedure based on Kronecker dictionary has much less strict requirements for coherence than the classic matching pursuit method.In other words, the former has the higher recovery bound [26].In addition, by using the structural information of signal, the reconstruction algorithm based on Kronecker dictionary performs robustly under the condition of high coherence.

Multiway Sparsity
One-dimensional expansion of a tensor is defined as stacking all elements along the mode-1, e.g., y = vec (Y ) = vec Y (1) , where mode is considered as the order of a tensor.Then the Tucker decomposition of Y is expressed as where × n denotes the mode-n product of tensor with matrix [37], which can be written as It is indicated from ( 10) that a Tucker decomposition of a multiway data is equivalent to a linear representation with Kronecker dictionary.Definition 1 (Multiway sparsity).If y in the Tucker model ( 10) is K-sparse with respect to The columns of A n are the atoms along each mode.
For an N-dimensional signal, if it has a K-sparse representation, the core tensor of Tucker decomposition has K nonzero elements, i.e., where "•" is the outer product operator.Because the nonzero entries of scatter coefficients are non-uniformly distributed and always gather together in the practical application of SAR imaging, this feature results in the highly structured echo signal, by which more rapid algorithm may be constructed.
By Definition 2, the nonzero elements of core tensor are concentrated in a subtensor S (I 1 , I 2 , . . . ,I N ).Then the relationship between multiway structured sparsity and the vector signal sparsity can be deduced: It should be noted that the definition of structured sparsity here is different from the definition of block sparsity used for one-dimensional signal.For the latter, a coefficient vector is called block K-sparse if it consists of a series of concatenate segments, named blocks, and the number of segments without nonzero Euclidean norm is less than K [38,39].

Sparse Reconstruction Based on Kronecker Dictionary
Compressed sensing can reconstruct the source signal from the undersampled measurements.Applying the technique to SAR super-resolution imaging still faces many challenges.The sparse imaging method is mainly to solve an ill-posed linear inverse problem by using the sparsity of radar cross section (RCS).For the large scale data, the classic l 0 procedure will be rigorously challenged by many computation problems.For instance, considering a scene with 1024 × 1024 observation points, the classic reconstruction approach will deal with data in length up to 1048576, which leads to a tremendous linear system, consuming computation resource significantly.From the analysis above, it is known that the target scene has the sparsity and the scatterers usually clustered together to blocky areas.Therefore, considering the structured sparse feature, the SAR imaging can be treated as a sparse reconstruction problem with Kronecker dictionary.A multiway structured sparse matching pursuit algorithm can be constructed to estimate the scatter coefficients with much fewer iterations.
Firstly, the SAR echo signal must be decoupled, for which the tensor data Y ∈ R P 1 ×P 2 ×•••×P N , and dictionaries along each mode A n ∈ R P n ×M n , n = 1, 2, . . ., N have to be constructed.Our purpose is reconstructing the nonzero scatter coefficients from Y, i.e., finding a (K 1 , K 2 , . . . ,K N )-structured sparse representation corresponding to A n (n = 1, 2, . . ., N).Because of the sparsity of target scene, only a fewer rows need to be randomly extracted from the dictionary.
Let B n = A n (:, I n ) , n = 1, 2, . . ., N, and then the Tucker decomposition of the original signal can be rewritten as where s nz ∈ R K is the vectorization of all nonzero elements.Consequently, the solution of the problem can be described as , the solution of ( 14) can be written as s nz = B T B −1 By.
Efficient calculation steps can be deduced by Cholesky decomposition [40].
Algorithm 1 Tucker decomposition based SAR imaging algorithm.
Input: decoupled echo signal Y; the number of maximum nonzero entries k max ; dictionaries along each mode {A 1 , A 2 , ..., A N }; the threshold of error .Output: scattering coefficients estimate Ŝ.
In iterative process, make sure that the nonzero elements of core tensor are always contained within a small index set.Finally all the nonzero elements and the corresponding index set are obtained and the expected scatter coefficients can be recovered as well.The algorithm flow is represented in Algorithm 1.

Algorithm Complexity Analysis
The first step of iteration needs to calculate R × 1 A T 1 (: . In addition, the absolute maximum optimization should be calculated, which needs 2M N operations.Comparing with this, the classic matching pursuit method, such as orthogonal matching pursuit (OMP) [30,32], requires the computation cost 2M N P N + 1 in this step.
With the application of the Cholesky decomposition, the computational complexity of estimating s nz can be reduced efficiently in the iteration.Even if every mode need to be updated, the total computational cost can not exceed 2Nk (P + 1) + Nk 2 + 2Nk N+1 .By comparison, the computational cost of OMP algorithm is k N P N + 3k N in this step.
In the residual calculation step, it is needed to compute and subtract it from Y.The corresponding calculations are P N (2k + 1).The ones of the OMP algorithm in this step, by contrast, are P N 2k N + 1 .
With regard to iterations, for a (K 0 , K 0 , . . . ,K 0 )-structured sparse problem, the maximum iterations is NK 0 .But OMP algorithm needs K N 0 times iteration for the same problem.

Numerical Simulations and Discussion
This section will demonstrate several experiments to verify the proposed method.The simulations are classified into two parts, one is spotlight SAR imaging realization, including ideal point-scatterer target setting and real data experiments; the other is performance comparisons between the proposed algorithm and the existing methods, including error curves and computational cost.We choose three algorithms as references: SVA [28], OMP [30] and compressive sampling MP (CoSaMP) [31].
The SVA algorithm is a kind of apodization filtering SAR imaging method, which is designed to achieve sidelobe control.However, the performance of such classic algorithm is susceptible by noise, and the resolution is limited, especially for adjacent targets.Both OMP and CoSaMP belong to sparse reconstruction algorithms, which have the advantage of super-resolution due to the robustness to sidelobe disturbance and noise.The hardware configuration for the simulation system is: Intel Core i7-5500U 2.4 GHz, 8 GB RAM, Windows 10.

Realization of SAR Imaging
In this section, two groups of simulation experiments are presented to verify the effectiveness of the proposed method.One investigates the SAR imaging of ideal point targets and the other considers the practical measured data.The simulation scenario is stated as Table 1.This section includes two settings of simulations: The first group investigates the imaging results via modifying the signal-to-noise ratio (SNR) under the condition of fixed number of scatterers.The second one inspects the imaging resolution via changing the the number of scatterers with fixed SNR.
(1) Different SNRs We investigate the effectiveness of the proposed method according to different noise conditions, from the low to the high SNR configurations.Five ideal point scatterers are placed around the centre of target scene in this simulation scenario.
The Figures 2 and 3 show the SAR imaging realizations at the situation of the 3 dB and 30 dB SNR, respectively.These two simulations demonstrate the imaging quality of the proposed method.For the low SNR, the imaging result of SVA shows severe ambiguity owing to the noise instability of this procedure.In comparison, the proposed algorithm performs much better, and the two matching pursuit algorithms also show good results.For the highest SNR, the imaging quality of all algorithms approach to the same.
(2) Different number of scatterers In this simulation scenario, the SNR is fixed at 5 dB.All the scatterers are gathered around three target clumps and such scatter centers are distributed along the diagonal line of the scene.
At the situation of scatterer number being 30, the SAR imaging realizations are shown in Figure 4.When the number of scatterers is set to 150, the SAR imaging realizations are shown in Figure 5. From the two groups of simulations above, it is indicated that for the small number of scatterers, all the methods are able to reconstruct target scene well.Compared with SVA, the sparse-reconstruction-based methods obtain more explicit images.It is worth noting that simulation results of the proposed method contain few noisy points.As the scatterer number increases, the scatterers imaging results appear fuzzy.Under the same scenario, the reconstructed scene of the proposed method achieve the most accurate images.The matrix A is constructed as the indication of the Equation ( 3), and Q = γ • M rows are extracted from it to form a projection matrix.Here, γ = 0.5 is taken.Note that the number of observations Q = O (K log M), where K is the sparsity.Therefore, the sparsity can be expressed as λQ/ log M. K is designed as 200 in this experiment.
In the following, we investigate the results at different azimuth angles under the condition at a fixed pitch angle firstly, and observe the results at different pitch angles under the condition at a fixed azimuth angle secondly.The above experiments demonstrate the effectiveness of the proposed method in practical implementation.With the comparison between the selected algorithms, the proposed method obtains the better imaging results.When the azimuth central angle is 100 • , the imaging result of SVA is dim and blurred, but the vehicle contour is clearly visible with the proposed method.For the experiments of different azimuth angles, although the classic matching pursuit approaches, i.e., OMP and CoSaMP, perform much better than SVA, the problem of imaging target ambiguity still exists.For instance, the CoSaMP imaging result has a remarkable bright streak at 90 • of the azimuth central angle, as shown in Figure 7. Besides, as indicated by the arrows in Figures 6 and 8, CoSaMP shows raster lines over the target image at the central region, which will disturb the recognition of targets significantly.Under the same scenario, there exist incorrect scattering points on the sides of the target in the OMP imaging results, which are circled in white for each images.By contrast, the proposed method does not have such problems in all the experimental tests.With the comparison between the experiments of different pitch angles, the superiority of the proposed method is demonstrated.As shown in Figure 11, the imaging result of SVA contains many noisy points appearing blurred; OMP outputs several outliers which are marked by white circles; for CoSaMP, a part of the scatter points in the central are lost, which is indicated by arrow C; but the proposed method obtains the clear and accurate imaging result.As shown in Figures 9 and 10, the imaging results of CoSaMP contain raster lines over the targets which are indicated by arrow A and B, respectively; the other three methods achieve the similar performance as demonstrated above.Although the two classic matching pursuit methods can obtain clearer images than SVA, but both of them result in errors of the estimation beyond tolerance.By contrast, the proposed method can achieve the better imaging results.This is true especially as the pitch angle increases.

Performance Comparison
This section focuses on two performance indicators of the proposed method and the reference methods, one is the root mean square error (RMSE) curve, the other is the computational resource requirements.

RMSE
The RMSE of SAR imaging model is defined as follows: where ŝl is the estimate of scatter coefficients in the lth Monte-Carlo independent trial.
(1) Different SNRs The configuration of this experiment is the same as Section 4.1.1.The SNR range is 3∼30 dB, and L, the number of independent trials, is set to 500.The RMSE curves are shown in Figure 12.
It is demonstrated from the RMSE curves that the performance of the proposed method is better than the others under the condition with low SNR.As the SNR increasing, the RMSE curves of the reference methods are approximately approaching to the one of the proposed method.Because of the sensitivity to noise, the reconstruction error of SVA algorithm is considerable in low SNR.In this respect, the sparse-reconstruction-based methods are robust to noise.The simulation results show that the performance of CoSaMP is slightly worse than that of OMP.The reason is that CoSaMP only chooses 2k max atoms at each iteration, which can help to improve the efficiency of the algorithm, but may reduce the success rate of reconstruction.In comparison, the proposed method performs much better, even if the noise interference is serious.This is because the Kronecker information of the echo signal is fully utilized.(2) Different number of scatterers The number of scattering points in the target area will directly affect the imaging quality.In the following experiment, the performance of the proposed method versus different number of scatterers will be investigated.The configuration of this simulation is the same as Section 4.1.1.The SNR is set to 5 dB, and the RMSE curves are shown in Figure 13.The simulation result indicates that the performance of proposed method is superior to the reference methods with all the test scenarios.When the numbers of scatterers is less than 30, the performance of all reference methods seems to be close.With the number of scatterers increasing, the performance of SVA, limited by the resolution of itself, becomes worse.The sparse-reconstruction-based methods with super-resolution maintain better robustness even as the scatterers become dense.Exploiting the structured feature of target, the proposed method provides the advantage for imaging, and achieves effective sidelobe control.Of all the curves of RMSE versus scatterer quantity, the one of the proposed method grows most slowly and always has the lowest RMSE value, according to the curves shown in Figure 13.Even the number of scatterers is larger than 90, the RMSE of the proposed method is still lower 3 dB than that of the SVA method.

Computational Resource Requirement
This group of experiments investigate the requirements for computing resources of different methods.The evaluation indicator is CPU time.The configuration of this experiment is the same as in Section 4.2.1.The statistical results for the case of different SNRs are shown in Table 2.The last experiment considers the case of different number of scatter points.Setting the configuration be the same as in Section 4.2.1, the statistical results are shown in Table 3. From the statistical results in Tables 2 and 3, it is shown that the fastest algorithm is SVA, the second fastest one is the proposed method, followed by CoSaMP and OMP in sequence.Although SVA is the fastest algorithm, its performance is the worst among all these methods, especially in the low SNR case, which is shown in Figure 12.Compared with this, OMP and CoSaMP perform better, but both the two algorithms require a large amount of system resources.Since OMP algorithm needs to match every atom at each iteration step, it runs extremely time-consuming, which means it unsuitable for the real-time applications.So far as that is concerned, the proposed method of this paper could not only run faster one order of magnitude than CoSaMP, but also ensure the better performance than the latter, which determines its prominent advantages in practical applications.

Conclusions
In this paper, a tensor decomposition algorithm for the spotlight SAR imaging is developed based on Tucker model which formulizes SAR echo signal with Kronecker constraint.Under the common scenario in SAR imaging, the scatterers usually cluster together and the targets are sparsely located.This structured sparse feature provides a prerequisite for the proposed method.Not only has this method better robustness for signal reconstruction, but also produces less algorithm complexity and requires smaller computation resource as well, due to exploiting the Kronecker information in each mode sufficiently.
In comparison with the existing SAR imaging methods, such as the SVA and classic matching pursuit algorithms, the proposed method performs much better with respect to robustness about noise and sidelobe disturbance.Numerical analysis has demonstrated that the recovery upper bound of the multiway sparse reconstruction based on structured sparsity is larger than that of the classic matching pursuit methods.Both the theoretical analysis and simulation results also indicate that the proposed algorithm requires much less computational resources than the reference methods, which is by far the most important in practical applications.

Figure 1 .
Figure 1.The geometry of spotlight SAR for point scatterers.

4. 1 . 2 .
SAR Imaging for Practical Measured Data This group of experiments deal with real data, measured by a spotlight SAR.The observed target is a crawler-type vehicle.The radar data are obtained from different azimuth and pitch angles.The angle aperture is 5 • .The center frequency is 9 GHz and the bandwidth is 1 GHz.The number of observation angle samples and frequency samples are both set to 101, thus the number of total pixels of image is M = 10201.

Table 1 .
System parameters setting for spotlight SAR.

Table 2 .
Average running time of each algorithm versus different SNRs.

Table 3 .
Average running time of each algorithm versus different number of scatterers.