3y ago

33 Views

2 Downloads

1.51 MB

10 Pages

Transcription

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 56, NO. 7, JULY 20183953A Modified EM Algorithm for ISAR ScattererTrajectory Matrix CompletionLei Liu , Feng Zhou , Member, IEEE, Xueru Bai , Member, IEEE, John Paisley, Member, IEEE,and Hongbing Ji, Senior Member, IEEEAbstract— The anisotropy of radar cross section of scatterersmakes the scatterer trajectory matrix incomplete in sequentialinverse synthetic aperture radar images. As a result, factorizationmethods cannot be directly applied to reconstruct the 3-D geometry of scatterers without additional consideration. We proposea modified expectation-maximization (EM) algorithm to retrievethe complete scatterer trajectory matrix. First, we derive themotion dynamics of the projected scatterer, which approximatesan ellipse. Then, based on the estimated ellipse parameters usingthe known data of each scatterer trajectory, we use the Kalmanfilter to initialize the missing data. To address the limitations of atraditional EM, which only considers the rank-deficient characteristics of the scatterer trajectory matrix, we propose to augmentEM by using both the known rank-deficient and elliptical motioncharacteristics. Experimental results on simulated data verify theeffectiveness of the proposed method.Index Terms— Inverse synthetic aperture radar (ISAR), matrixcompletion, Modified expectation-maximization (EM), scatterertrajectory.I. I NTRODUCTIONINVERSE synthetic aperture radar (ISAR) can acquire2-D high-resolution images of airplanes and space targets and has found wide applications in military and civilareas thanks to its all-day and all-weather observation capacity [1]–[3]. High-resolution 2-D images can be obtainedby range compression techniques in the range dimensionand using rotation movement relative to the radar line-ofsight (LOS). However, 2-D ISAR images are the projectionof the target’s 3-D geometry on the ISAR imaging planes.Therefore, only the 2-D projected geometry of an object canbe obtained from 2-D ISAR images, which creates limitationsManuscript received September 12, 2017; revised January 16, 2018;accepted February 12, 2018. Date of publication April 17, 2018; date ofcurrent version June 22, 2018. This work was supported in part by theProject funded by the China Postdoctoral Science Foundation under Grant2016M602775 and Grant 2017M613076, in part by the National NaturalScience Foundation of China under Grant U1430123, Grant 61631019,Grant 61471284, Grant 61522114, and Grant 61571349, in part by theFoundation for the Author of National Excellent Doctoral Dissertation ofChina under Grant 201448, and in part by the Young Scientist Awardof Shaanxi Province under Grant 2015KJXX-19 and Grant 2016KJXX-82.(Corresponding authors: Lei Liu; Feng Zhou.)L. Liu and F. Zhou are with the Ministry Key Laboratory of Electronic Information Countermeasure and Simulation, Xidian University, Xi’an 710071,China (e-mail: liulei xidian@163.com; fzhou@mail.xidian.edu.cn).X. Bai is with the National Laboratory of Radar Signal Processing, XidianUniversity, Xi’an 710071, China (e-mail: xrbai@xidian.edu.cn).J. Paisley is with the Department of Electrical Engineering and the DataScience Institute, Columbia University, New York, NY 10027 USA (e-mail:jpaisley@columbia.edu).H. Ji is with the School of Electronic Engineering, Xidian University,Xi’an 710071, China (e-mail: hbji@xidian.edu.cn).Color versions of one or more of the figures in this paper are availableonline at http://ieeexplore.ieee.org.Digital Object Identifier 10.1109/TGRS.2018.2817650for target feature extraction when used in target recognitionand classification problems. Hence, methods for performing3-D ISAR imaging have attracted much attention. Traditional3-D ISAR imaging algorithms can generally be divided intotwo classes: interferometric ISAR (InISAR)-based 3-D imaging [4]–[6] and sequential ISAR image-based 3-D geometryreconstruction [7]–[9]. Compared with InISAR techniques,3-D imaging algorithms based on 2-D sequential ISAR imagesare independent of complicated antenna array design andcomplex signal processing. Traditional monostatic ISAR issufficient to obtain sequential ISAR images through wideangle and long-term observation. Therefore, research on3-D geometry reconstruction using sequential ISAR imageshas been a focus in recent ISAR imaging problems.In 2-D sequential ISAR image-based 3-D geometry reconstruction algorithms, matrix factorization is performed on thescatterer trajectory matrix to reconstruct the 3-D geometry ofthe target and the projection matrix. Therefore, the extractionand association of scatterers are the basis for using factorization methods. As is well known, the reflected electromagneticwaves of the observed target can be modeled as the summation of the electromagnetic waves of different scatterersin the microwave band. Hence, the ISAR image is sparseand a certain number of scatterers dominate the total echoenergy. Therefore, scatterers can be extracted using sparsemodeling techniques, e.g., estimation of signal parameters viarotational invariance techniques [10]–[12], the root multiplesignal classification algorithm [13], and amplitude and phaseestimation algorithm [14]. The trajectory of each scatterer canthen be found using the Markov chain Monte Carlo (MCMC)data association [15]–[17]. However, an important issue isthat the associated trajectory matrix may be incomplete. Forwide-angle and long-term observation, the anisotropy of theradar cross section (RCS) of each scatterer can generallynot be ignored. Yet, the RCS has such a small partialobservation angle that the scatterer often cannot be extractedsuccessfully [10]–[12], and so the trajectories of anisotropyscatterers are incomplete. In other words, the existing intervalsof different scatterers are different, with the result that thescatterer trajectory matrix is incomplete. Therefore, matrixfactorization methods [18], [19] cannot be applied directlyto reconstruct the 3-D geometry of the observed target. Thismakes retrieval of the incomplete scatterer trajectory matrixthe key problem for the sequential ISAR image-based 3-Dgeometry reconstruction.Performing 3-D geometric reconstruction using 2-Dsequential ISAR images is a structure-from-motion (SFM)problem [20]. Since the 2-D ISAR image is a projection0196-2892 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications standards/publications/rights/index.html for more information.

3954IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 56, NO. 7, JULY 2018of the 3-D scatterer position on the 2-D imaging plane,the rank of the noiseless, complete scatterer trajectory matrixis known to equal three a priori [2]. Therefore, the estimationof missing data can be modeled as a rank-deficient retrievalproblem. General methods for the missing data retrieval problem in the SFM have been considered [20]. For example,Tomasi and Kanade [18] factorized the largest complete submatrix to obtain a partial shape and motion matrix. Thisinitial solution then grows by one row or one column ata time to approximate the full matrix. In the work of [21]and [22], each column of the incomplete rank-deficient matrixis modeled as an affine subspace. The unknown data arerecovered by least squares regression on that subspace. Theselection strategy of submatrices is studied in [22] to improvethe robustness of the algorithm, and [23]–[25] introduce theexpectation-maximization (EM) and row–column algorithmsto estimate the missing data. The limitation of these methodsis that they only consider the low-rank property while ignoringmotion characteristics in ISAR imaging scenarios, such asthe motion dynamics of the scatterer trajectory. Although theretrieval of the incomplete scatterer trajectory using the abovementioned methods satisfies the rank-deficient property, it maynot coincide with the practical motion of the scatterer.Motivated by the EM algorithm approach proposedin [23]–[25], we propose a novel initialization method anda modified EM method to estimate the missing data of thescatterer trajectory matrix. Compared with the traditional EM,we will make use of the motion characteristics of the scatterertrajectory. We first analyze the 2-D ISAR imaging modeland the projection mechanism. Then, we will show that theprojected trajectory of each scatterer approximates an ellipseunder the stationary motion of the observed target. We thendefine the scatterer state vector as the projected range andcross-range positions at each instant and derive the dynamicsof the scatterer state vector transition. After that, we use thebidirectional Kalman filter to initialize the missing data bymaking use of the motion characteristics of each scatterertrajectory. Finally, we propose a modified EM algorithm toestimate the missing data. Since the motion of each scattereris considered in the E-step of modified EM, the estimatedmissing data has an explicit physical meaning and correspondsto a trajectory.The remainder of this paper is organized as follows.Section II introduces the ISAR imaging model and analyzesthe motion physics of each scatterer trajectory. We review thetraditional EM algorithm in Section III, while also pointingout its limitations for the scatterer trajectory retrieval problem.Then, Section IV presents the proposed bidirectional Kalmanfilter to initialize the missing data and a modification to EMis discussed in detail. Experimental results based on simulateddata are presented in Section V to show the effectiveness ofthe proposed method. Section VI concludes this paper.II. ISAR I MAGING M ODEL AND A NALYSISOF S CATTERER T RAJECTORYThe ISAR can obtain long-term and wide-angle echoes oftargets having steady orbits and stationary motion. SequentialISAR images of the observed target can be acquired byFig. 1.ISAR observation and imaging model.dividing the data into overlapping subapertures and performingISAR imaging algorithms on each [2]. This ISAR observationand imaging model is shown in Fig. 1, in which translationalmotion is assumed to have been compensated for. The targetcan therefore be modeled as rotating around a fixed axis. Letthe 3-D coordinate system be O XY Z and the target rotatearound O Z with rotation angular velocity ω. O correspondsto the center of rotation motion, which is also set to bethe origin of the coordinate system. The azimuth angle andpitch angle of the radar LOS are θl and φl , respectively.Since the observed target has steady motion, we assume thatthe radar LOS remains unchanged relative to the rotationaxis. In other words, θl and φl remain constants for theentire ISAR imaging interval. If the wide-angle data aredivided into K subapertures, then K ISAR images can beacquired using traditional range-Doppler imaging algorithm.The interval of adjacent subapertures is T . In the microwavefrequency band, the ISAR echoes of one target can be viewedas the summation of echoes of the dominant scatterers, whichare usually corners, edges, and flat surfaces [26]. The 3-Dgeometry of the target can be obtained by the 3-D reconstruction of those dominant scatterers. We will assume thatthe observed target consists of N scatterers. For any arbitraryscatterer pn [x n , yn , z n ]T , the projected 2-D positions inthe kth image frame can be expressed as described in thefollowing [2].ttLet u nk and v nk indicate the projected 2-D positions of onescatterer in one image frame. These u and v values denote thescatterer positions projected on range dimension and crossrange dimension, respectively. The subscript tk indicates theobservation instant for the kth ISAR image frame and thesuperscript n denotes the scatterer index. Then tk unv ntk sin(ωtk θl )cos(ωtk θl )0 cos(φl ) cos(ωtk θl ) cos(φl ) sin(ωtk θl ) sin(φl ) xn yn .(1)zn

LIU et al.: MODIFIED EM ALGORITHM FOR ISAR SCATTERER TRAJECTORY MATRIX COMPLETION3955We note that tk tk 1 T is a constant. To better analyzethe characteristics of the scatterer’s motion, we can rewrite (1)in the following equivalent way: u kn v nk z n sin(φl ) sin(ωtk θl ) cos(ωtk θl ) x n . cos(ωtk θl ) sin(ωtk θl )yncos(φl )(2)In the above-mentioned equation, we simplify u tnk and v ntk byremoving t and only show its measurement index k. The lefthand side (LHS) of (2) is equal to the measured values ofthe initial position vector pn rotating around the O Z -axisaccording to ωtk θl . Since the observed target rotatesaround O Z , the value of (x n )2 (yn )2 (rn )2 remainsunchanged. The right-hand side is a rotation transformationof pn . Hence, the LHS satisfies the same constraint(v nk z n sin(φl ))2 (rn )2cos2 (φ L )(u kn )2(v nk z n sin(φl ))2 1.(rn )2(rn cos(φ L ))2(u kn )2 k 1unv nk 1 Illustration of incomplete trajectory matrix.the product of the projection matrix R and the scatterer3-D geometry matrix S W RS(3)As is evident,to an ellipse function, which (3) correspondsTshows that u kn , v nk measured for different image frames alllie on an ellipse. The center of the ellipse is [0, z n sin(φl )]T .The major and minor axes are rn and rn cos(φl ), respectively.Hence, the ellipse centers of different scatterers are eachlocated on the same line parallel to the radar LOS. Theratio of the minor axis to the major axis is cos(φl ), and theprojected scatterer trajectory is a circle when the radar LOS isperpendicular to the rotation axis, i.e., φl 0. Equations (2)and (3) describe the linear relationship between a scatterer’strue 3-D position and its measured 2-D projection.The motion dynamics of a scatterer’s 2-D position can bederived as follows: k sin(ω · T )cos(ω· T)un cos(φl )v nkcos(φl ) sin(ω · T ) cos(ω · T ) Fig. 2.A sin(ω · T ) . cos(φl )cos(ω · T )z n sin(φl ) z n sin(φl )b(4)Equation (4) shows that A and b are independent of thespecific image frame indexed by k. Therefore, once sequentialISAR images have been obtained, the analytical expressions ofA and b are determined. If we can extract the 2-D positions ofall N scatterers at each image frame, then we will be able toconstruct the scatterer trajectory matrix W by using scattererposition association processing [15], where u kn0 k K, 0 n NWkn (5)k KvnK k 2K , 0 n Ndenotes the instantaneous position of each scatterer. ClearlyW R2K N . According to (1), W can be expressed by(6)whereRk· [sin(ωt θ),cos(ωt θl ), 0] ,0 k K klk cos(φl ) cos(ωtk K θl ), cos(φl ) sin(ωtk K θl ), (7) sin(φl )] ,K k 2KS·n [x n , yn , z n ]T ,0 n N.(8)Clearly, R R2K 3 and S R3 N . Rk· denotes the kth rowof the projection matrix R, and S·n represents the nth columnof the scatterer’s 3-D geometry matrix S. The rank of W is3 under the noiseless scenario. As is well known, the RCSof anisotropy scatterers changes with the observation angleand the transmitted signal frequency. Therefore, anisotropyscatterers may not be extracted for all image frames, since thecorresponding RCS is too small. After obtaining sequentialISAR images, one should extract the dominant scatterers first.Then, one can perform the MCMC data association methodof [15] on the scatterer position set to obtain the associatedscatterer trajectory. However, the trajectory matrix of scatterersis partially observed and incomplete, as shown in Fig. 2.3-D geometry reconstruction using sequential ISAR imagesis actually an SFM problem, and so factorization methods canbe applied directly to the complete scatterer trajectory matrixto recover the projection matrix and the 3-D positions of thescatterers. However, since the measured trajectory matrix isincomplete, factorization method cannot be directly applied,but rather more complex algorithms are required. In thiscase, the rank of the scatterer trajectory matrix W is stillapproximately three, even after accounting for noise, and sothe matrix retrieval task is a low-rank problem. Many methodshave been proposed to deal with the rank-deficient retrievalproblem for image processing, for example those based on theEM algorithm [23]–[25]. In Section III, we introduce the EMprocedure and discuss the limitations of existing techniquesfor the scatterer trajectory matrix retrieval problem.

3956IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 56, NO. 7, JULY 2018Fig. 3. Limitation of trajectory estimation using traditional EM. (a) Scatterer model. (b) Retrieved trajectory under 50% missing data ratio. (c) Retrievedtrajectory under 20% missing data ratio.III. E XPECTATION -M AXIMIZATION A LGORITHMAs discussed in Section II, the true rank of the scatterertrajectory W is known to be three and is approximated wellby this rank in moderately noisy settings. Thus, the estimationof W is actually a rank-deficient matrix retrieval problem.The EM algorithm has had previous success for this problem [23]–[25]. In this section, we will first review the EMprocedure and then discuss potential limitations when appliedto the scatterer trajectory retrieval problem. Let M denote themask matrix where m i j is either 1 or 0; m i j 1 indicatesthat Wi j is observed, whereas m i j 0 indicates that Wi jis missing. The scatterer trajectory matrix can be expressed W M, where denotes the Hadamard (elemenas Wtwise) product. Therefore, the missing scatterer positions are The problem is to estimate the missingexpressed by 0 in W. entries of W. One way to approach this problem is to constructthe following function:E( W) min S3W ( W W) M2(9)where W denotes the estimated scatterer trajectory matrixwhose rank is strictly three. The EM algorithm is performed intwo alternating steps: 1) the E-step estimates the missing data and 2) the M-step estimatesof the scatterer trajectory matrix W the rank three matrix W that best matches the “complete” data.The detailed procedure is presented as follows.1) E-Step (Estimation of the Missing Scatterer Positions):Given the previous estimate of Wq 1 where q denotes i j : m i j 0}the iterate number, the missing data {W i j of W.of W are simply the corresponding entries WqWe will then build a trajectory matrix W , whose entry i j of theW i j is equal to the corresponding entry W i j ifmatrix W if Wi j was observed, or its estimate WWi j is unknown i j if m i j 1WW ij (10) i j if m i j 0Wor, in matrix notation qW W M Wq 1 1 M .(11)2) M-Step (Estimation of Rank Three Matrix): We nowqhave obtained the complete trajectory matrix W withthe estimates of the missing scatterer positions from theE-step. The rank three matrix Wq that best matchesqW in the Frobenius norm is then computed from theqsingular value decomposition of W , as in q Wq S3W(12)where S3 denotes the set of rank three matrices. Hence,qW S3 represents estimating the rank three matrixwhich minimizes the cost function (9).The EM algorithm is a common method that can dealwith general rank-deficient matrix retrieval problem withmissing entries. However, the general EM algorithm onlymakes use of the rank-deficient characteristics, whereas additional constraints for specific problems are not considered,e.g., the motion characteristics of a scatterer’s 2-D projectedtrajectory. Hence, although the scatterer trajectory matrix canbe retrieved satisfying rank three, it may not coincide with thereal motion dynamics of the scatterer trajectory. We presenta simple simulation to illustrate this limitation in Fig. 3. The3-D geometry of 16 scatterers is shown in Fig. 3(a). After thestandard EM algorithm, we compare the retrieved trajectorieswith the theoretical ones in Fig. 3(b) and (c). The red solidlines with “*” denote the theoretical trajectory. The black solidand dashed lines with “ ” represent the known data and theestimated data, respectively. Fifty percentage of each scatterertrajectory is missing in Fig. 3(b) and 20% of each scatterertrajectory is missing in Fig. 3(c). The three largest singularvalues of the retrieved trajectory matrix using EM are fa

Index Terms—Inverse synthetic aperture radar (ISAR), matrix completion, Modiﬁed expectation-maximization (EM), scatterer trajectory. I. INTRODUCTION I NVERSE synthetic aperture radar (ISAR) can acquire 2-D high-resolution images of airplanes and space tar-gets and has found wide applications in military and civil

Related Documents: