Tucker Modeling based Kronecker Constrained Block Sparse Algorithm

Abstract

This paper studies synthetic aperture radar (SAR) imaging problem which the scatterers are often distributed in block sparse pattern. To exploiting the sparse geometrical feature, a Kronecker constrained SAR imaging algorithm is proposed by combining the block sparse characteristics with the multiway sparse reconstruction framework with Tucker modeling. We validate the proposed algorithm via real data and it shows that the our algorithm can achieve better accuracy and convergence than the reference methods even in the demanding environment. Meanwhile, the complexity is smaller than that of the existing methods. The simulation experiments confirmed the effectiveness of the algorithm as well.

1. Introduction

 The synthetic aperture radar (SAR) is a promising radar that is employed to reconstruct a two-dimensional or three-dimensional images of desired objects. SAR has been widely applied in many fields in recent decades [1], [2]. The quality of the SAR image is often influenced by many factors. For instance, it is restricted by the Nyquist sampling limit, the receiver needs a high sampling rate for the wideband waveform, which causes to the large resource consumption for data storage and processing. Besides, imaging quality of the traditional methods, such as polar format algorithm (PFA) [3], [4], may be degraded significantly by noises and sidelobe disturbance. Recently, compressed sensing (CS) technology has been utilized for SAR imaging, by which the necessary information can be obtained by the sampling rate which is much lower than Nyquist sampling limit [5]. The SAR imaging methods in the CS framework have the robustness to noise and sidelobe disturbance [5]-[7]. However, these studies merely consider the spatial sparsity of the observation, which don’t exploit the structural characteristics of scatterers.It is the fact that the scatterers can be aggregated together as nonzero blocks [2], [7]. In 2016, a CS based SAR imaging algorithm is proposed using the Bayesian learning procedure with consideration of the structured sparsity constraint [8]. But this kind of procedure requires the interesting signal signal data rigidly to match with the model since it introduces the multi-level random model.

 

Fig. 1. Figure illustration of SAR geometry for point scatterers.

 In order to solve the existing problems and obtain a better SAR imaging performance, this paper investigates a tensor decomposition-based CS method which can make use of the block sparsity of scatterers, where tensor refers to the multidimensional array [9]. Usually, the signals in practice are multidimensional and highly structured [10]. Tensor decomposition is appropriate to handle such higher-order data, so it has been studied extensively in recent years [10]-[12]. Because there exist curse of dimensionality in tensor signal processing [10], [13], [14], researches are mainly focused on how to take advantage of the sparsity of signals [15], [16]. In the SAR imaging area, there has few work combining sparse reconstruction with the tensor decomposition until now. Inspired by the research of higher-order CS in hyperspectral imaging [17], [18], this paper formulizes the echo signal of SAR as a tucker model and then develops an imaging algorithm by using the block sparsity of target scene, as shown in Fig. 1. Experiments are conducted to validate the our proposed algorithm.

 The remainder of our paper can be organized as follows. In the Section 2, we introduces the received signal model of point scattering target; Section 3 describes the tensor modeling based on the tucker decomposition; The proposed algorithm is introduced in Section 4, and the Section 5 will demonstrate several typical simulations to evaluate our proposed method, and makes detailed comparisons with conventional schemes in performance; Finally, we will conclude this paper.

2. Signal Model

 Some necessary notions are explained as follows. The italic letter denotes the scalar, e.g., x ; the bold italic letters indicate vector and matrix, respectively, e.g., x and X ; the calligraphic capital letter indicate the tenso, e.g., \(\mathcal{X}\) ; the script letter presents the index of a set, e.g., \(\mathfrak{X}\) ; the \(\ell_{0}\)-norm \(\|\boldsymbol{x}\|_{0}\) equals to the nonzeros number in x ; \(\|\mathcal{X}\|_{F}\) is the Frobenius-norm of tensor \(\mathcal{X}\) ; the operator ⊗ is denoted as Kronecker product operation between two matrices; the ×_n presents the mode-n product operation, where the mode denotes the order of a tensor [9].

 Consider the spotlight SAR and the geometry is depicted in Fig. 1. The parameters θ and \(\phi\) denote respectively, the azimuth and pitch angle of the radio wave. \(R\left(x_{0}, y_{0}, z_{0}\right)\) represents the distance between the target and antenna at \(\left(x_{0}, y_{0}, z_{0}\right)\); \(R_c\) indicates the distance between the center of the scene and the antenna. Then the received signal can be formulated as

\(y(t, \theta, \phi)=\sum_{x, y, z} g(x, y, z) s(t-\tau(x, y, z))\)       (1)

where \(g(x, y, z)\) indicates the scattering function of the target located at \((x, y, z)\) ; \(s(t)\) is the autocorrelation function for transmitted signal. After transformation to the frequency domain, the Eq. (1) is transformed as

\(Y(f, \theta, \phi)=C(f) \sum_{x, y, z} g(x, y, z) \exp \left(H_{x, y, z}(f, \theta, \phi)\right)\)       (2)

where \(H_{x, y, z}(f, \theta, \phi)=(4 \pi j / c)(x f \cos \theta \cos \phi+y f \cos \theta \sin \phi+z f \sin \theta)\) is the phase factor [19]. The constant term \(C(f)=S(f) \exp \left(-j 2 \pi f \tau_{c}\right)\) can be ignored, where \(S (f)\) is the Fourier transform of the autocorrelation function of transmitted signal; \(\tau_{c}=2 R_{c} / c\) denotes the propagation time, where c is the speed of light.

 Because of the coupled parameters in (2), it is necessary to implement an interpolation method for decoupling [2]. Thus polar coordinates are transformed into Cartesian coordinate system [20]. The numbers of grid points in \(x, y, z\)directions are \(M_{1}, M_{2}, M_{3}\) and the grid points in \(f, \theta, \phi\) directions are \(P_{1}, P_{2}, P_ 3\) . Defining the steering vector of the scatterer at \(x_{j}, y_{k}, z_{l}\) as

\(d_{j, k, l}=\left[b_{j, k, l}(1,1,1), \ldots, b_{j, k, l}\left(P_{1}, P_{2}, 1\right), \ldots,\left(P_{1}, P_{2}, P_{3}\right)\right]^{T}\)

where \(b_{j, k, l}\left(p_{1}, p_{2}, p_{3}\right)=\exp \left((4 \pi j / c)\left(x_{j} u_{p 1}+y_{j} v_{p 2}+z_{l} w_{p 3}\right)\right)\), the steering matrix \(\boldsymbol{D} \in \mathbb{C}^{P \times M}\) is \(\boldsymbol{D}=\left[d_{1,1,1}, \ldots, d_{M_{1}, M_{2}, 1}, \ldots, d_{1,1, M_{3}}, \ldots, d_{M_{1}, M_{2}, M_{3}}\right]\), where \(M=M_{1} \times M_{2} \times M_{3}, P=P_{1} \times P_{2} \times P_{3}\). Then the received signal model can be represented as

\(y=D s\)       (3)

where \(\boldsymbol{s}=\left[g\left(x_{1}, y_{1}, z_{1}\right), \ldots, g\left(x_{M_{1}}, y_{M_{2}}, z_{1}\right), \ldots, g\left(x_{1}, y_{1}, z_{M_{3}}\right), \ldots, g\left(x_{M_{1}}, y_{M_{2}}, z_{M_{3}}\right)\right]^{T}\). Considering the \(\left(p_{1} p_{2} p_{3}, m_{1} m_{2} m_{3}\right)\) -th element of \(\boldsymbol D\) is \(\exp \left((4 \pi j / c)\left(x_{m_{1}} u_{p_{1}}+y_{m_{2}} v_{p_{2}}+z_{m_ 3} w_{p_{3}}\right)\right)\) , the steering matrix can be expressed as \(\boldsymbol{D}=\boldsymbol{D}_{3} \otimes \boldsymbol{D}_{2} \otimes \boldsymbol{D}_{1}\) . Hence, the received signal model can be formulated as

\(y=\left(D_{3} \otimes D_{2} \otimes D_{1}\right) s\)       (4)

Algorithm 1. The proposed algorithm using tensor-based block sparse constraint.

Input: Observation data \(\mathcal{Y}\) ; the number of maximum nonzero entries \(k_{\max }\) ; dictionaries \(\left\{\boldsymbol{D}_{1}, \boldsymbol{D}_{2}, \cdots, \boldsymbol{D}_{N}\right\}\) ; the threshold of error \(\mathcal{E}\) ; Output: Estimation of the scattering coefficients \(\hat{S}\) ; Initialize the index set \(\mathcal{I}_{1}=\varnothing, n=1,2, \ldots, N\) ; Let \(\mathcal{R}=\mathcal{Y}, \hat{S}=0, k=1\); While \(\left|\mathcal{I}_{1}\right|\left|\mathcal{I}_{2}\right| \cdots\left|\mathcal{I}_{N}\right| \leq k_{\max }\) and \(\|\boldsymbol{R}\|_{F}>\varepsilon\) do 1. \(\left\{i_{1}^{k}, \cdots, i_{N}^{k}\right\}=\arg \max _{\left\{i_{1}^{k}, \cdots, i_{N}^{k}\right\}}\left|\boldsymbol{R}_{\times 1} \boldsymbol{D}_{1}^{T}\left(:, i_{1}\right)_{\times 2 \cdots \times N} \boldsymbol{D}_{N}^{T}\left(:, i_{N}\right)\right|\) 2. Update the index set \(\left|\mathcal{I}_{n}\right|=\left|\mathcal{I}_{n}\right| \cup\left\{i_{n}^{k}\right\}, i=1,2, \cdots, N\) 3. \(\boldsymbol{B}_{n}=\boldsymbol{D}_{n}\left(:, \mathcal{I}_{n}\right)\) 4. Substitute \(\boldsymbol{B}_{n}(n=1,2, \cdots, N)\) into (8), combining Cholesky decomposition, obtain the estimates of nonzero coefficients \(\boldsymbol{s}_{n z}\) ; 5. Reshape \(\boldsymbol{s}_{n z}\) to the tensor form \(\boldsymbol{S}_{n z}\) ; 6. Computing the residual \(\boldsymbol{R}=\mathcal{Y}-\boldsymbol{S}_{n z} \times_{1} \boldsymbol{B}_{1} \times_{2} \boldsymbol{B}_{2} \cdots \times_{N} \boldsymbol{B}_{N}\) 7. \(k=k+1\) End While Reconstruct the scattering coefficients \(\hat{\boldsymbol{S}}\) from \(\boldsymbol{S}_{n z}\) and index set \(\mathcal{I}=\left\{\mathcal{I}_{1}, \mathcal{I}_{2}, \cdots, \mathcal{I}_{N}\right\}\), satisfying that \(\hat{S}\left(\mathcal{I}_{1}, \mathcal{I}_{2}, \cdots, \mathcal{I}_{N}\right)=S_{n z}\).

 

3. Tucker decomposition for SAR model

 The tensor signal is a multidimensional data with more than two modes. The tensor can be expanded into a one-dimension vector, i.e., \(y=\operatorname{vec}(\mathcal{Y})=\operatorname{vec}\left(Y_{(1)}\right)\). Then the Tucker decomposition of \(\mathcal{Y}\) is \(\mathcal{S} \times_{1} \boldsymbol{D}_{1} \times_{2} \boldsymbol{D}_{2} \times_{3} \ldots \times_{N} \boldsymbol{D}_{N}\) [10], which can be written as

\(y=\left(D_{N} \otimes D_{N-1} \otimes \cdots \otimes D_{1}\right) s\)       (5)

Because a K-sparse representation of signal y exists a dictionary \(\boldsymbol{D}\), we have

\({y}={D} {s}, \text { s.t. }\|{s}\|_{0} \leq K\)       (6)

where K is referred to the sparsity. Since the columns of the dictionary is more than the rows, the system is undetermined. Because of the sparsity feature of the desired signal, then a solution can be achieved under the condition [21],

\(K<\frac{1}{2}\left(1+\frac{1}{\mu(D)}\right)\)       (7)

in which \(\mu(\boldsymbol{D})=\max _{i \neq j}\left|\left\langle\boldsymbol{d}_{i}, \boldsymbol{d}_{j}\right\rangle\right|\) denotes the coherence coefficient of the matrix \(\boldsymbol{D}\). When the dictionary possesses Kronecker structure, i.e., \(\boldsymbol{D}_{N} \otimes \boldsymbol{D}_{N-1} \otimes \cdots \otimes \boldsymbol{D}_{1}\), the coherence coefficient should be \(\mu(\boldsymbol{D})=\max \left\{\mu_{1}, \mu_{1}, \cdots, \mu_{N}\right\}\), with \(\mu_{n}=\mu\left(\boldsymbol{D}_{n}\right), n=1,2, \cdots, N\) [21].

4. The proposed algorithm for SAR imaging

 According to the Kronecker structure of SAR imaging model and the blocks structural characteristics of scatterers, in this section, we will propose our SAR imaging algorithm based on the block sparsity feature with Kronecker dictionary.

Definition 1. Multiway block sparsity.

 A tensor data \(\mathcal{Y}\) is called as \(\left(K_{1}, K_{2}, \cdots, K_{N}\right)\)-block sparse, if only \(K_{n}(n=1,2, \cdots, N)\) columns in the factor matrices of its tucker decomposition need to be computed, i.e.,

\(\mathcal{Y}=\mathcal{S} \times_{1} \boldsymbol{D}_{1} \times_{2} \boldsymbol{D}_{2} \times_{3} \ldots \times_{N} \boldsymbol{D}_{N}, \text { s.t. }\)

\(\boldsymbol{s}_{i, i_{2}, \cdots, i_{N}}=0 \forall\left(i_{1}, i_{2}, \cdots, i_{N}\right) \notin\left\{\mathcal{I}_{1}, \mathcal{I}_{2}, \cdots, \mathcal{I}_{N}\right\}\)

in which \(\mathcal{I}_{n}=\left\{i_{n}^{1}, i_{n}^{2}, i_{n}^{K_{n}}\right\}\) is a subset of the index set for mode-n (\(n=1,2, \cdots, N\))

 The Definition 1 indicates that nonzero elements in core tensor are concentrated within a subtensor \(\mathcal{S}\left(\mathcal{I}_{1}, \mathcal{I}_{2}, \cdots, \mathcal{I}_{N}\right)\) . Therefore, \(y=\operatorname{vec}(\mathcal{Y})\) is K-sparse (\(K=\prod_{i=1}^{K} K_{i}\)) with respect to the Kronecker dictionary \(\boldsymbol{D}=\boldsymbol{D}_{N} \otimes \boldsymbol{D}_{N-1} \otimes \cdots \otimes \boldsymbol{D}_{1}\). For the typical underdetermined system\(y=\boldsymbol D x\), where the coefficient vector x contains a series data segments with \(\left\{d_{m}\right\}_{m=1}^{M}\), and the m-th segment (block) can be defined as \(\boldsymbol{x}[m]\). Define the \(\ell_{0}\)-norm of x as \(\|\boldsymbol{x}\|_{0, b}=\sum_{m=1}^{M} f_{b}\left(\|\boldsymbol{x}[m]\|_{2}>0\right)\), where \(\|\boldsymbol{x}[m]\|_{2}\) is the \(ℓ_2\)-norm of the m-th data block; \(f_{b}\left(\|\boldsymbol{x}[m]\|_{2}>0\right)\) represents an indicator function which returns 1 if \(\|\boldsymbol{x}[m]\|_{2}>0\), otherwise returns 0. If \(\|x\|_{0, b}, x is block k-sparse.. Defining \(\boldsymbol{B}_{n}=\boldsymbol{D}_{n}\left(:, \mathcal{I}_{n}\right)\), \(n=1,2, \cdots, N\), the tucker decomposition for the original signal (5) can be expressed as

\(\hat{\boldsymbol{y}}=\left(\boldsymbol{B}_{N} \otimes \boldsymbol{B}_{N-1} \otimes \cdots \otimes \boldsymbol{B}_{1}\right) \boldsymbol{s}_{n z}\)       (8)

where \(\boldsymbol{s}_{n z} \in \mathbb{R}^{K}\) is the vectorization of all nonzero scattering coefficients. As a result, the solution of the problem can be given as

\(s_{n z}=\underset{u}{\operatorname{argmin}}\left\|\left(B_{N} \otimes B_{N-1} \otimes \cdots \otimes B_{1}\right) u-y\right\|_{2}^{2}\)       (9)

Denoting \(\boldsymbol{B}_{N} \otimes \boldsymbol{B}_{N-1} \otimes \cdots \otimes \boldsymbol{B}_{1} \) as \( \boldsymbol{B}\) , then the solution of equation (9) can be derived as \(\boldsymbol{s}_{n z}=\left(\boldsymbol{B}^{T} \boldsymbol{B}\right)^{1}\). A more efficient calculation step can be added by means of Cholesky decomposition [22]. It is necessary to ensure that the nonzero elements of the core tensor must be included in a small subset of index during iterations. As the output, all nonzero elements and the index set can be obtained; the expected scattering coefficients will be recovered as well. The algorithm flow is shown in Algorithm 1.

 The complexity of the proposed algorithm is reduced to \(\mathcal{O}\left(M^{N} P+k P^{N}\right)\), comparing with OMP algorithm which has the complexity of \(\mathcal{O}\left((M P)^{N}+(k P)^{N}\right)\)[21]. For a \(\left(K_{0}, K_{0}, \cdots, K_{0}\right)\)-block sparse multiway decomposition problem described as Section 3, a sparse representation can be guaranteed, if \(\left(K_{0} \mu\right)^{N}<2-\left(1+\left(K_{0}-1\right) \mu\right)^{N}\), where \(\mu=\max \left\{\mu_{1}, \mu_{2}, \cdots, \mu_{N}\right\}\), and the iterations of the proposed method is less than \(N K_{0}\), while the OMP is \(K_{0}^{N}\) . Here, it is necessary to note that the iteration number of our proposed method is between \(K_0\) and \(NK_0\). After the maximum correlated atom is choiced in step.1 in Algorithm 1, the position of the \((M \times M \times \ldots \times M)\) multiway array determines the indices \(i_{1}^{k}, \cdots, i_{N}^{k}\) to be added to the current \(\mathcal{I}_{n}\) subsets. There are two conditions about the computational cost. For the better condition, a new index is incremented to each mode-n dictionary at every iteration. For the worse condition, only one new index is incremented to each mode at every iteration.As a consequence, the minimum number of iterations is \(K_0\) and the maximum number of iterations is \(NK_0\). This indicates that the proposed algorithm can estimate the scattering coefficients with much fewer computation burdens than the classic matching pursuit algorithm.

5. Simulation experiments

 In this section will evluate the performance of our proposed algorithm through several numerical experiments. The experiments are classified as two parts: 1) SAR imaging implementations, including ideal point-scatterer simulation and real data experiment; 2) performance comparisons of the proposed algorithm with the existing methods, including two indicators of root mean square error (RMSE) and computational cost. Three algorithms are chosen as reference methods: PFA [3], [4], OMP [21] and CoSaMP [23]. All the experiments are conducted in the personal computer equipped with Windows 10 OS, Intel Core i7-5500U 2.4GHz CPU, and 8GB RAM.

Table 1. Simulation parameters.

 

5.1. Implementation for SAR Imaging

 To test the effectiveness of the proposed method, two groups of experiments will be presented in the following, including the simulated ideal point targets and the practical measured data in SAR imaging. All the parameters are set as follows: the center frequency is set to 9 GHz; the bandwidth is 1 GHz and the frequency resolution is 0.01 GHz. What’s more, the angle aperture and the angle resolution are all set to 5o; and the subsampling rate is 50%.

 1) Simulations of ideal-point target: The simulation scenario presents 20 ideal point scatterers distributed randomly. When the SNR is set as 3 dB and 30 dB, the SAR imaging implementations are shown in Fig. 2. These two experiments are provided to show the imaging quality of the proposed method. For the high SNR scenarios, the imaging results show no distinct difference among all algorithms. But for the low SNR case, the imaging results of OMP and CoSaMP display many outliers, meanwhile PFA shows severe noisy ambiguity. In contrast, the imaging results of the proposed algorithm show less noisy points and output the correct images.

 2) Practical measured data applications: This group of experiments is carried out based on the real data of spotlight SAR. The observation target is a crawler-type engineering vehicle. The observation direction is ( , ) (90 ,50 ) o o φ θ = . We set the observation angle samples number and the frequency samples number as both 101, thus the total pixels of scene image are 10201. In this work, the sparsity is set to 200. The SAR imaging implementations are shown in Fig. 3. These experiments demonstrate that the result of PFA shows obvious imaging blur. Due to the sidelobe effect, PFA cannot perform well for adjacent targets or weak scattering targets. The classic matching pursuit methods, such as OMP and CoSaMP, also generate the imaging ambiguity, despite they perform much better than PFA. By contrast, the proposed algorithm eliminates these drawbacks and generates the correct images under the same scenario, which displays the effectiveness of our proposed algorithm in practical application.

 

Fig. 2. SAR imaging results for ideal point scatterer model of different SNRs using four algorithms of PFA, OMP, CoSaMP and our proposed algorithm.

5.2. Performance Comparison

The definition of RMSE for SAR imaging model is

\(\operatorname{RMSE}\left(\hat{s}_{l}\right)=\sqrt{\frac{1}{L} \sum_{l=1}^{L}\left(\left\|\boldsymbol{y}-\boldsymbol{A} \hat{s}_{l}\right\|_{2}^{2} /\|\boldsymbol{y}\|_{2}^{2}\right)}\)       (10)

where ˆ l s denotes the estimate values for scattering coefficients. The SNR range is 3~30dB, and 𝐿𝐿 is set to 500. The RMSE curves of the simulations for 10 scatterers and 50 scatterers, are shown in Fig. 4.

 The RMSE curves demonstrate the performance of the algorithms. The reconstruction error of PFA is obviously big at low SNR since its sensitivity to noise. From the results, we can find that matching pursuit methods of OMP and CoSamp are more robust to noise. The proposed algorithm outperforms the reference methods at the same SNR in both the cases of small number and large number of scattering points. Since the Kronecker structure information of signal is fully exploited, the RMSE of the proposed algorithm is still quite small even under the serious noise interference.

 

Fig. 3. SAR imaging for the practical measured data.

 

Fig. 4. Performance comparison for two kinds of scatterers versus SNRs.

 2) Computational cost: This simulation experiment surveys the requirements for computing resources. CPU operation time is considered as the evaluation indicator. The configuration of this experiment is the same as Section 5-1. The results of computational cost are shown in Table 1.

Table 2. Average operation time of different SNRs.

 

 The statistical results present that the PFA and the proposed algorithm are faster than OMP and CoSamp. However, the accuracy of PFA is the worst among all reference methods as demonstrated in Section 5-B(1). By contrast, OMP and CoSaMP perform better, but they require more computing resources. Especially OMP algorithm runs much time-consumingly since it must match every atom at each iteration step, which makes it unsuitable for the real-time applications. The proposed algorithm could not only run faster, but also can achieve the better imaging performance.

6. Conclusion

 In this paper, we investigate a spotlight SAR imaging algorithm for point scattering targets based on Tucker tensor decomposition. In this work, the geometry feature, i.e., the block sparse distribution of scatterers for a typical SAR target scene has been exploited by the tensor modeling. For this purpose, the scattering parameters estimation is formulized as a multiway sparse reconstruction problem. Since this SAR imaging model has Kronecker structure to depict the block sparsity of scatterer distribution reasonably, the iterative convergence of the algorithm is obviously accelerated and the robustness against the interference is enhanced according to the theoretical analysis and numerical experiments, compared with the reference methods. These features imply that the proposed algorithm can be well applied in practical applications with real-time requirements.