Temporal rate up-conversion of synthetic aperture radar via low-rank matrix recovery

TEMPORAL RATE UP-CONVERSION OF SYNTHETIC APERTURE RADAR VIA LOW-RANK MATRIX RECOVERY Minh Dao1 , Lam Nguyen2 , and Trac D. Tran1 1

Department of Electrical and Computer Engineering, The Johns Hopkins University, MD 21218 2 U.S. Army Research Laboratory, MD 20783 ABSTRACT

The radar data to form synthetic aperture radar (SAR) imagery is normally transmitted and received by moving platforms like aircraft or vehicles. In many situations, the platforms move at high speed; which reduces the number of sampling records collected to the synthetic aperture, hence degrades the quality of the reconstructed SAR images. Therefore, it is necessary to develop an algorithm that is capable of increasing the temporal frequency rate of the received data. In this paper, we propose a novel technique to generate intermediate records from the existing ones by a locally-adaptive low-rank matrix recovery framework. The system first fills in the blank records using a bi-directional motion estimation scheme. The initialized aperture records are then refined by a robust low-rank matrix completion algorithm using the reference from neighborhood clean records. Experiments demonstrate that the proposed method provides comparative results when up-converting the aperture rate by a factor of two or four, both in mean square error of the raw SAR signal and PSNR performance of the recovered SAR images. Index Terms— Synthetic aperture radar (SAR), rate upconversion, low-rank matrix completion, robust PCA. 1. INTRODUCTION A synthetic aperture radar transmits and receives electromagnetic waves by a sensor attached to a platform moving along the cross-range direction to generate the synthetic aperture [1]. The system shoots and collects signals at a constant pulse repetition interval (PRI) along the radar path to produce equally spaced aperture records. The collected data is then processed by a backprojection algorithm to form high resolution SAR imagery. The more number of sensing records the system receives, the sharper and higher quality the formed image is. An aperture rate up-conversion algorithm is therefore important. Generally, if we can double the number of sensing records, the vehicle carrying the radar system can boost the travelling speed by a factor of two. The problem of inserting new sensing records between the existing ones can be considered as the problem of completing full missing aperture records along the temporal domain and

978-1-4799-2341-0/13/$31.00 ©2013 IEEE

has not been thoroughly researched for SAR data. Most of the existing completion techniques solve for the cases of recovering missing or corrupted samples [2, 3] or interpolating the signal in the spectral domain. However, almost no methods have been successful in tackle the challenging problem of recovering whole missing sensing records without damaging image structures and details. Recent advances on low-rank matrix recovery techniques have opened a new trend in approaching completion problems. Matrix recovery problem is a robust framework to recover incomplete or corrupted data by extracting lowdimensional structures from high-dimensional observations. Matrix completion [4] and robust principal component analysis (robust PCA or RPCA) [5] are two highly applicable lowrank matrix recovery techniques in which matrix completion retrieves missing elements while RPCA recovers an underlying low-rank matrix from its sparse but grossly corrupted entries. These two problems have been beneficial in solving a wide range of applications including background modeling, target tracking [5], image alignment [6] or video completion [7] problems. Nevertheless, these techniques cannot be directly applied for matrices with fully missing rows or columns. Furthermore, most of the existing applications assume global low-rank property of the input data which is not commonly true for the case. With consideration of these observations, we propose a new method to take full advantage of matrix completion and RPCA methods to effectively interpolate full aperture records of SAR data. The way to construct dictionary of the proposed method is different from existing techniques in the sense that it uses the reference from the neighboring signal, while for most of other SAR sparsity-driven completion approaches [2, 3], the dictionary is formed based on the transmitting signal of the SAR system. The advantage of the local-based dictionary formation is that it can make use of the high temporal correlation of SAR data. The algorithm is then summarized by two main steps: 1. the initialization step fills in new records by a simple motion-estimated scheme; and 2. the refinement step updates the records by a proposed robust matrix completion using confident weights obtained from the first step. For simple instruction, the detailed algorithm in this paper only describes the case of temporal upsampling by a factor of two.

2358

ICIP 2013

The rest of the paper is organized as follows. Section 2 briefly gives an overview of matrix completion and RPCA techniques and introduce the proposed robust matrix completion algorithm. Section 3 describes the proposed temporal rate up-conversion for SAR data in a two-stage framework. The experimental results for the simulated SAR data are evaluated in section 4 and we conclude this paper in section 5. 2. LOW-RANK MATRIX RECOVERY In the most general form, low-rank matrix recovery problem consists of recovering a low-rank matrix X ∈RN1 ×N2 from a set of M linear measurements: y = A(X) where A : RN1 ×N2 → RM is a linear map. In order to reconstruct X, one would like to find the simplest model that fits the low-rank observations: M in Rank(X) s.t. y = A(X) X

(1)

The above rank-minimization problem is a non-convex and intractable problem. Under some mild conditions [4], however, the rank-minimization problem can be relaxed to the convex problem of nuclear norm minimization: M in kXk∗ s.t. y = A(X) X

(2)

where the nuclear matrix norm kXk∗ is defined as sum of all singular values of the matrix X. A highly applicable subset of low-rank matrix recovery problems in (2) is the matrix completion problem where the goal is to recover an unknown matrix from a subset of its entries. Typically, given an incomplete observation matrix Y = X|Ω , where Ω is the index set of available entries of X, we want to recover back the original matrix X with the prior knowledge that X is low-rank. Here, A is the linear operator that sets unobserved entries to zero. Again, to achieve X, nuclear norm minimization is proposed as follows [4]: M in kXk∗ s.t. Y = X|Ω X

(3)

Another recent landmark result introduced by Candes et al investigates the problem of robust principal component analysis (RPCA) [5]. The question is how to efficiently recover a low-rank structure X from the corrupted observed version Y = X + E, where E is a sparse matrix which contains most of zero or near-zero entries but the other magnitudes can be arbitrarily large. To separate X and E, a principle component pursuit strategy is proposed [5]: M in kXk∗ + λ kEk1 s.t. Y = X + E X,E

Algorithm 1 Robust matrix completion (RMC) Input: Data input Y, weighting matrix WE 1. Z0 = 0, E0 = 0, µ0 > 0, ρ > 1, k = 1 2. While not converged, do 3. (U, Σ, V) = svd(Y − Ek + µ−1 k Zk ); 4. Xk+1 = US µ−1 [Σ]VT ; k

5. Ek+1 = S µ−1 WE (Y − Xk+1 + µ−1 k Zk )|Ω k 6. Zk+1 = Zk + µk (Y − Xk+1 − Ek+1 ) 7. µk+1 = ρµk 8. k = k + 1 9. end while where S  (a) = max(|a| − , 0) sgn(a) is the element-wise soft-thresholding operator. Outputs: (Xk , Ek ). cases, we have no prior information of the support locations of the outlier entries. Therefore, we propose a robust matrix completion (RMC) algorithm which can solve for the case the input data contains both missing elements and sparse noise components. Mathematically, the RMC problem recover the low-rank matrix X and the sparse matrix E from the incomplete observation Y = (X + E)|Ω : M in kXk∗ + λ kEk1 s.t. Y = (X + E)|Ω X,E

(5)

Furthermore, if the prior knowledge of where the sparse components in E are more likely to appear is provided, i.e. each element in E is appointed with a weighting parameter representing how much assurance that element is not an outlier. The global weighting parameter λ is then replaced by a component-wise weighting matrix WE in the cost function: M in kXk∗ + kWE ◦Ek1 s.t. Y = (X + E)|Ω X,E

(6)

where ◦ denotes the Hadamard (pointwise) product of two matrices. The weighted RMC framework (6) can be efficiently solved by a number of methods in which we promote to use Augmented Lagrange Multipliers (ALM) [8] because of its fast convergence property. The corresponding augmented Lagrangian function of (6) is: L(X, E, Z) = kXk∗ + kWE ◦Ek1

(7) µ + h(Y − X − E)|Ω , Zi + k(Y − X − E)|Ω kF 2 where Z is the Lagrangian multiplier, µ > 0 is a penalty parameter, h·i is the inner product and k.kF denotes the Frobenius norm. The unconstrained minimization of L(X, E, Z) can be solved by updating one variable at a time, iteratively. The detailed algorithm is summarized in Algorithm 1.

(4)

where λ is a positive weighting parameter, and k.k1 is the l1 norm of a matrix. This problem in some sense can be viewed as a generalization of matrix completion. In fact, if ¯ , we set entries of E such that Eij = −Xij with (i, j) ∈ Ω then RPCA turns into the matrix completion problem. However, this assumption generally does not hold since in many

3. TEMPORAL RATE UP-CONVERSION VIA LOW-RANK MATRIX RECOVERY Let denote S = [s1 , s2 , ..., sN ] ∈ RK×N as the full SAR data where sn ∈ RK (n = 1, 2, ..., N ) is the received aperture record at a given time n along the radar travelling path.

2359

Without the loss of generality, we assume that all the odd numbered records (n = 1, 3, ...) are received by the SAR system, and all the even numbered records (n = 2, 4, ...) are to be interpolated. The aperture rate up-conversion problem can be formulated as the problem of estimating S from the low temporal resolution observation SL = [s1 , s3 , ...]. The proposed algorithm to upscale the number of records can be defined by two main steps: the initialization step and the refinement step. In the initialization step, the interpolated record is simply determined by bi-directional motion estimations from the two nearby records. This step also assigns to each samples in the record a confident weight based on both the estimation similarity and local consistency. In the refinement step, every sample in the initialized record is updated by a locally-adaptive low-rank matrix recovery algorithm. 3.1. Sensing record initialization step. Let sn (n < N ) be an unknown sensing record that we want to interpolate; sn−1 and sn+1 be the previous and future observed records. For each patch pn (k) of size p centered at the sample k (k = 1, 2, ..., K) in sn , we search for the best match trajectory passing through the patch pn (k) in term of smallest sum of absolute difference (SAD) from a patch pn−1 (k − d) in the aperture sn−1 to a patch pn+1 (k + d) in the aperture sn+1 , within a search region d ∈ [−∆d, ∆d]. The initial value of sample k is then calculated as the average ˆ n−1 (k − of the two center samples in the best match patches p ˆ ˆ This bi-directional motion estimation ˆ n+1 (k + d). d)and p process can be illustrated as figure 1(a). The initialization step certainly provides many inaccurate motion estimations and need to be updated. Intuitively, the errors are normally caused by estimations with large SAD as well as those with inconsistent trajectory motion with the neighborhoods. Hence, the motion estimations should be further classified and processed rather than directly accepted for the result. A sample can be categorized as totally unreliable if its SAD of the motion estimation is larger than a threshold or the trajectory going through its location is very different with those passing its nearby pixels. These samples are not used for the refinement process and can be considered as incomplete samples. The remaining set can contains both reliable and unreliable samples but cannot be directly classified. However, each sample k (k = 1, 2, ..., K) is assigned with a weighting value w(k) which is the summation of the similarity of the best motion estimation going through that sample w(M ) (k) (defined as one over SAD) and the consistency of the trajectory motion at the sample with its surrounding neighborhoods w(C) (k): w(k) = w(M ) (k) + w(C) (k). The larger a component of w is, the more certainty it is a correct motion estimation. 3.2. Low-rank matrix recovery refinement step.

(a)

(b)

Fig. 1. (a) Bi-directional motion estimation, and (b) low-rank matrix construction ple values in the record sn , but also a support set of unreliable samples and a vector of confident weights w. In the refinement step, the initialized record sn is first divided into distinct patches of the same length. For each patch pi of length L (say L = 16) in sn , a matrix Di ∈ Rk is constructed such that the first column of Di is the vector pi and the remaining columns are the patches of the same length sliding all neighboring area of the previous and future records. The reference records in this process are not constrained in the right previous or next aperture record, but can be extended to a number r of available records in both temporal directions. Because of the temporal correlation within neighborhood aperture, each of the reference record should contain one patch that has similar underlying structure as the patch pi , and each patch in any other record should also have similar patches from the neighbors. Thus, the columns of Di should stay in very low-dimensional subspace, making the matrix Di be low-rank. The low-rank matrix construction is depicted in figure 1(b). The first column of the low-rank structure Di is a section pi of the initialization of sn , hence it can be incomplete. Let Ωi be the set of missing elements of Di corresponding to the totally unreliable motion estimations in pi determined by the initialization step. In the available index set of Di , i.e. the complement set Ωi , it is expected that some of the elements in Di are sparse noise from the low-rank representation of Di , especially the atoms from the first column. Di then can be formulated as Di = (Ai + Ei )|Ωi , where Ai is a lowrank matrix and Ei is a sparse noise. The outliers in Ei come from the uncorrected motion estimations processed in step one and are automatically detected and updated in the RMC algorithm. The other columns of Di are the estimations from the clean apertures. Therefore, the probability that they behave as sparse noise is low. With this analysis, we define the weighting matrix Wi with the first column corresponding to the weights in w of the elements from the first column of Di , and the remaining columns having some small nonnegative value. Then the RMA algorithm for record interpolation is formulated as the following optimization:



Mi ini Ai ∗ + Wi ◦Ei 1 s.t. Di = (Ai + Ei )|Ωi (8) A ,E

The outputs of the initialization step are not only the sam-

The recovery of (8) via algorithm 1 both recovers the missing

2360

100

50

50

50

100

100

100

150

150

150

200

200

250

250

300

300

350

350

400

400

300 Down range

Sample amplitude

200

400

200 250 300

500 350 600

400

700

450 450 50

100

150 Aperture records

200

250

300

50

100

150

200

(a)

250 300 Cross range

350

400

450

500

500

−3

−3

−3

−4

x 10

3.5

10

1.5

100

150

200

−3

300

350

400

450

500

50

100

150

200

250

300

350

400

450

500

(b)

Fig. 4. SAR image formulations for double rate upconversion by (a) bi-cubic interpolation (PSNR = 22.50 dB) and (b) the proposed method (PSNR = 28.74dB).

−4

x 10 x 10

x 10 10

1.5 3

1

250

(a)

(b)

x 10 x 10

3

500 50

Fig. 2. (a) Original SAR data in 2D form, and (b) Formed SAR image. 3.5

450

500

800

1

5

5

50

50

0

100

100

150

150

200

200

250

250

300

300

350

350

400

400

0.5

0.5

2.5

2.5 0

0

0

2 −0.5

2 −0.5 −1 1.5

−1

−5 320

340

360

440

460

1.5

480

1

1

0.5

0.5

0

0

−0.5 Original data Bi−cubic method −1

−5 320

340

360

440

460

480

−0.5 Original data Proposed meth

−1 0

100

200

300

400

500

(a)

600

700

800

0

100

200

300

400

500

600

700

800

450

(b)

450

500

500 50

Fig. 3. Reconstructions of the aperture record at n = 40 by (a) bi-cubic interpolation (RM SE = 4.32e − 4) and (b) the proposed method (RM SE = 1.37e − 4) in comparison with the baseline signal for the double rate up-conversion case. elements of pi and updates the other samples by replacing pi with the first column of the output low-rank matrix Ai . This refinement process is then applied for all patches in sn to complete the full interpolation of one aperture record.

4. EXPERIMENTAL RESULTS In this section we test the performance of the algorithm by a simulation SAR data. The data is configured to represent many point targets randomly locating throughout the scene and having different amplitudes. The number of records in the aperture is set to N = 300. This SAR data is processed through a backprojection image formation [9] to form the SAR image. Both of the SAR data and the formed image are used as the ground truth for performance comparison purpose. Figure 2(a) and 2(b) represent the simulated raw SAR data in 2D form and the reconstructed SAR image, respectively. In the first experiment, the SAR data is temporally reduced by a factor of two, i.e. all the records at even number times are eliminated. The missing records are then interpolated by both bi-cubic interpolation along the temporal direction and the proposed algorithm, and compared with the original data. The patch length used in the experiment is fixed at L = 16. Three nearest records in both previous and future directions are used as the references in the low-rank construction step. Figure 3 shows the reconstructions and root mean square errors (RMSE) from the record at time n = 40 by both bicubic interpolation and the proposed method in comparison with the original signal. While the bi-cubic interpolation result does not follow the ground truth in almost all targets, the

100

150

200

250

300

350

400

450

500

50

100

150

(a)

200

250

300

350

400

450

500

(b)

Fig. 5. SAR image formulations for four times rate upconversion by (a) bi-cubic interpolation (PSNR = 16.73 dB) and (b) the proposed method (PSNR = 23.24dB). output of the proposed method captures most peaks as well as side lobes. In some cases, the lobes are one sample away from the origins, but the magnitudes still reserve. This only results a little blurring of the corresponding targets in the formed SAR image. Figure 4 shows the recovered SAR images by two methods which demonstrate the outperformance in both visual quality and PSNR of the proposed algorithm. In the second experiment, the data is down-sampled by four time in the temporal direction and recovered back. The SAR image formation of the proposed technique (figure 5(b)) still outperforms bi-cubic interpolation technique (figure 5(a)). However, more noise are introduced and the targets are detected at a higher false alarm rate.

5. CONCLUSIONS In this paper, we propose a novel robust approach to incorporate low-rank matrix recovery technique to temporally rate up-convert synthetic aperture radar data by a two-step framework. In the first step, the record to-be-interpolated is initialized by a bi-directional motion estimation scheme. In the second step, a novel framework to construct low-rank structures based on the intrinsic high correlation along the temporal direction of SAR data is proposed. A robust matrix completion algorithm to solve for the case both missing elements and sparse outliers with weighting prior are present is also promoted to update the initialized record. The proposed method much outperforms bi-cubic interpolation and shows comparative performance with the original data.

2361

6. REFERENCES [1] M. Ressler, L. Happ, L. Nguyen, T. Ton, & M. Bennett, “The Army Research Laboratory ultra-wide band testbed radars”, Radar Conference (RADAR), 1995. [2] L. Nguyen, and T. D. Tran, “Robust recovery of synthetic aperture radar data from uniformly under-sampled measurements”, IEEE Conf. Geoscience and Remote Sensing Symposium (IGARSS), 2011. [3] L. Nguyen, and T. Do, “Recovery of missing spectral information in ultra-wideband synthetic aperture radar (SAR) data”, Radar Conference (RADAR), 2012. [4] E. J. Candes and B. Recht, “Exact matrix completion via convex optimization”, Foundations of Computational Mathematics, vol. 9, pp. 717-772, 2008.

[6] Y. Peng, A. Ganesh, J. Wright, and Y. Ma, “Robust alignment by sparse and low-rank decomposition for linearly correlated images,” IEEE Conf. Comp. Vision and Pattern Recog. (CVPR), 2010. [7] M. Dao, D. Nguyen, Y. Cao, and T. D. Tran, “Video concealment via matrix completion at high missing rates,” Proc. IEEE 44th Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, Oct. 2010. [8] Z. Lin, M. Chen, & Y. Ma, “The augmented lagrange multiplier method for exact recovery of corrupted lowrank matrices,” arXiv preprint arXiv:1009.5055, 2010. [9] L. Nguyen, “Signal and Image Processing Algorithms for the Army Research Lab Ultra-Wideband Synchronous Impulse Reconstruction (UWB SIRE) Radar,” ARL Technical Report, ARL-TR-4784, April 2009.

[5] E. J. Candes, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis”, Journal of the ACM, 58(3), 2011.

2362