Progress In Electromagnetics Research, Vol. 129, 287–313, 2012 CONVOLUTION BACK-PROJECTION IMAGING ALGO- RITHM FOR DOWNWARD-LOOKING SPARSE LINEAR ARRAY THREE DIMENSIONAL SYNTHETIC APER- TURE RADAR X.M.Peng*, W.X.Tan, Y.P.Wang, W.Hong, andY.R.Wu Science and Technology on Microwave Imaging Laboratory, Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China Abstract—General side-looking synthetic aperture radar (SAR) cannot obtain scattering information about the observed scenes which are constrained by lay over and shading effects. Downward-looking sparse linear array three-dimensional SAR (DLSLA 3D SAR) can be placed on small and mobile platform, allows for the acquisition of full 3Dmicrowaveimagesandovercomestherestrictionsofshadingandlay over effects in side-looking SAR. DLSLA 3D SAR can be developed for various applications, such as city planning, environmental monitoring, DigitalElevationModel(DEM)generation, disasterrelief, surveillance and reconnaissance, etc. In this paper, we give the imaging geometry and dechirp echo signal model of DLSLA 3D SAR. The sparse linear array is composed of multiple transmitting and receiving array elements placed sparsely along cross-track dimension. The radar works on time-divided transmitting-receiving mode. Particularly, the platform motion impact on the echo signal during the time-divided transmitting-receiving procedure is considered. Then we analyse the wave propagation, along-track and cross-track dimensional echo signal bandwidth before and after dechrip processing. In the following we extend the projection-slice theorem which is widely used in computerizedaxialtomography(CAT)toDLSLA3DSARimaging. In consideration of the flying platform motion compensation during time- divided transmitting-receiving procedure and parallel implementation on multi-core CPU or Graphics processing units (GPU) processor, the convolution back-projection (CBP) imaging algorithm is proposed for DLSLA3DSARimagereconstruction.Atlast,thefocusingcapabilities Received 13 May 2012, Accepted 7 June 2012, Scheduled 26 June 2012 * Correspondingauthor:XuemingPeng([email protected]). 288 Peng et al. of our proposed imaging algorithm are investigated and verified by numerical simulations and theoretical analysis. 1. INTRODUCTION Owing to side-looking geometry in general SAR systems, shading and lay over effects by trees, buildings, special terrain shapes will hide essential information of the observed areas, particularly in urban areas and deep valleys in mountain areas [1–3]. DLSLA 3D SAR acquires full 3D microwave images by wave propagation dimensional pulse compression, along-track dimensional aperture synthesis with flying platform movement and cross-track dimensional aperture synthesis with a sparse linear array [4]. And DLSLA 3D SAR observes nadir areas which means it overcomes restrictions of shading and lay over effects in side-looking SAR. Downward-Looking Imaging Radar (DLIR) [5] was first intro- duced by Gierull in 1999. He took advantage of aperture synthesis by platform movement along along-track dimension and linear array aperture synthesis along cross-track dimension to form a two dimen- sional aperture and transmitted a single frequency signal to obtain the two dimensional image. Nouvel et al. in ONERA [6–8] made use of chirp signal instead of single frequency signal and developed the concept of downward-looking 3D-SAR and developed a Ka band lin- ear array 3D-SAR system. Researchers in FGAN-FHR designed an Airborne Radar for Three-dimensional Imaging and Nadir Observa- tion (ARTINO) [4,9]. The system places a thinned linear array along cross-track dimension, it obtains three dimensional resolution by wave propagation pulse compression, aperture synthesis along the fly di- mension and beam-forming along cross-track dimension. The imaging algorithms of the two systems are not illustrated in detail and the 3D imaging results of the two systems have not been publicly reported yet [10]. DLSLA 3D SAR requires pretty fast A/D sampling devices and pretty high capacity data collection devices for echo data recording as there are more transmitting-receiving channels in DLSLA 3D SAR than in side-looking SAR. Sparse linear array with the time- divided transmitting-receiving is often applied to the 3D SAR which helps to reduce the complexity of the system. The motion of the platform during the time-divided transmitting-receiving procedure degrades the reconstructed image and the motion effect should be compensated in the imaging algorithm. Most of the published imaging algorithms suppose the echo signal without carrier frequency and be digitized directly, and take no consideration of A/D sampling Progress In Electromagnetics Research, Vol. 129, 2012 289 rate limitation, echo data voxel volume and flying platform motion effectduringthetime-dividedtransmitting-receivingprocedure. These algorithmsmainlygetthewavepropagationandalong-trackdimension compression first with Range Doppler algorithm, Chirp Scaling algorithmandRangeMigrationalgorithmlikegeneraltwodimensional SAR image reconstruction, then obtain the cross-track dimension compression via SPECAN, beam-forming, CS (Compressed Sensing), etc. [11–16]. These algorithms may be suited for multi-baseline 3D SAR whose echo signal is composed of two dimensional general SAR echo signal along multi-baseline. Multi-baseline 3D SAR systems put less burden on A/D sampling devices and data collection devices while require quite a long time to obtain the 3D echo signal, and the samplingintervalbetweenthemulti-baselineisnotuniform[13,15,16], which restricts the imaging algorithms on the basis of uniform Nyquist sampling theorem. What’s more, any change of the imaging scene destroys the coherence of the imaging scene, while the coherence of the imaging scene is very important for microwave imaging [2]. Thedechirpsignalprocessingtechniqueismainlyusedinspotlight mode SAR systems which helps to reduce the echo signal bandwidth. Dechirp signal processing technique mixes the echo signal with a reference function which is the echo signal from APC (Antenna Phase Center) to SC (Scene Center) [2]. The reduced bandwidth is only related to the range extension of the imaging scene. DLSLA 3D SAR operatesinnadirobservation. Andinwavepropagationdimension,the range gate is the undulation range of the imaging scene. This range gate is several hundred meters in urban areas or plane areas, and the corresponding range gate time delay is very small. While, the pulse of thetransmittingsignalisoftenlargerthantherangegatetimedelayin ordertoobtainenoughechopowerandechoSNR.Sothedechirpsignal processing technique is very suitable for DLSLA 3D SAR for which helps to reduce the A/D sampling bandwidth and echo data volume. Our work below is based on dechirp processing technique. Projection- slice theorem was introduced and used in CAT first [17]. Later, some researchers extended projection-slice theorem to spotlight mode SAR image reconstruction [18,19]. Knaell and Cardillo [20] discussed the radar tomography for the generation of three-dimensional images, but theyneglectedthecenterfrequencyshiftoftheconvolutionkernel, and theydidnotgivethecross-trackdimensionalimaginggeometryandthe three dimensional point spread function. What’s more, the algorithm they introduced takes no consideration of the platform motion effect duringthetime-dividedtransmitting-receivingprocedurewhichmeans it cannot be used for DLSLA 3D SAR imaging. In this paper, we offer a new point of view to DLSLA 3D SAR echo signal acquisition on the 290 Peng et al. basis of projection-slice theorem. And, the CBP imaging algorithm which compensates the flying platform motion effect during the time- divided transmitting-receiving procedure and contains the advantage of parallel implementation [21,22], is introduced for DLSLA 3D SAR. The structure of this paper is organized as follows. Downward- looking sparse linear array 3D SAR imaging geometry, echo signal model, platform motion effect during time-divided transmitting- receivingprocedureandechosignalbandwidthbeforeandafterdechirp are established in Section 2; then, in Section 3, projection-slice theorem and its application in DLSLA 3D SAR are described in detail; CBP imaging algorithm with platform motion compensation during time-divided transmitting-receiving procedure and the parallel implementation of the algorithm is discussed in Section 4; Section 5 gives some computer simulations and result analysis. Finally, conclusion is provided in Section 6. 2. DOWNWARD-LOOKING SPARSE LINEAR ARRAY THREE DIMENSIONAL SAR IMAGING GEOMETRY AND ECHO MODEL Thissectionissplittedintofoursubsections. First,wegivetheDLSLA 3D SAR imaging geometry and illustrate why DLSLA 3D SAR can gather three dimensional scatter information that avoided shading and lay over effects. Then the echo signal model with dechirp signal processing technique is introduced. Next the platform motion effect during the time-divided transmitting-receiving procedure is analyzed. At last, the echo signal bandwidth in every dimension before and after dechirp processing is analysed. 2.1. Downward-looking Sparse Linear Array Three Dimensional SAR Imaging Geometry As is shown in Fig. 1. X-axis is parallel to the flight path of the platform, Y-axis is parallel to the cross-track dimensional sparse linear array, Z-axis is perpendicular to the XY plane, the origin of the coordinateO isthecenterofthethreedimensionalimagingscene(O is also used as the reference point for dechirp signal processing). Q is the Antenna Phase Center (APC), Q(cid:48) is the projection point of Q on XY −−→ plane, P is the target in the imaging scene. APC flight path QX(cid:48) is −−→ parallel to X-axis. QO is the instantaneous reference range from APC −−→ to SC with the distance of R , OP is the instantaneous range from 0 −−→ SC to target with the distance of R, QP is the instantaneous slant Progress In Electromagnetics Research, Vol. 129, 2012 291 (a) (b) Figure 1. DLSLA 3D SAR imaging mode and imaging geometry. (a) DLSLA 3D SAR imaging mode. (b) DLSLA 3D SAR imaging geometry. −−→ −−→ −−→ range from APC to target with the distance of R . QP = QO +OP t (R = R +R) on the hypothesis of planar wave front. ψ is the angle t 0 −−→ between OQ and XY plane, called as grazing angle. ϕ is the angle −−→ betweenOQandZ-axis,calledasincidentangle. θistheanglebetween −−→ −−→ OQ(cid:48) and X-axis, called as slant angle. γ is the angle between QX(cid:48) 1 −−→ and QP, called as along-track dimension Doppler cone angle. γ is the 2 −−→ −−→ angle between QQ(cid:48) and QP, called as cross-track dimension Doppler cone angle. TheSARprocessorneedsthreekindsofinformationtoreconstruct animagefromsensorsignals. First, itneedspositioninginformationto locate the source of the return echoes in the imaging scene. Second, it needs spatial or angular resolution information to differentiate among return echoes from separate scatters or scene areas. Finally, it needs the intensity information of the signal associated with each other or scene areas in the target field. In the case of radar, the intensity is the radar cross section of individual scatter or the radar backscatter coefficient of distributed areas [2]. As shown in Fig. 1, general side- looking SAR measures range from APC to target and along-track Doppler cone angle γ . It can not distinguish targets on the yellow 1 292 Peng et al. region as they have the same range and the same along-track Doppler cone angle [3]. And grazing angle in side-looking SAR is small, especially in the far range gate unit, the red region is shaded by other targets, and SAR sensor can not obtain scattering information of this region[3]. DLSLA3DSARmeasuresrangefromAPCtotarget,along- track dimension Doppler cone angle γ , cross-track dimension Doppler 1 cone angle γ and the radar cross section of the imaging scene which 2 means it acquires full 3D microwave image of the observed areas. And, DLSLA 3D SAR operates nadir observation, therefore, layover and shading effects can be overcome. 2.2. Downward-looking Sparse Linear Array Three Dimensional SAR Imaging Echo Model The sparse linear array is composed of multiple transmitting array elements and multiple receiving array elements placed sparsely along the cross-track dimension and works on the time-divided transmitting- receiving mode. According to the equivalent phase center theorem, the sparse linear array with multiple transmitting array elements and multiple receiving array elements on the basis of time division can be equivalenttoauniformlineararraythateveryequivalentarrayelement transmits and receives signal by itself [12]. All the equations below are based on the equivalent phase center theorem as the range from APC to target is large enough that the phase error generated by equivalent phase center theorem can be ignored. The sparse linear array we design is based on ARTINO array distribution [3] as the length of equivalentarraywiththismethodislongerthanothermethods[23,24]. What’s more, ARTINO array distribution is quite easy for hardware implementation [25]. Suppose that there are L transmitting array 1 elements, L receiving array elements. The equivalent phase center 2 number is L L according to ARTINO array distribution [3]. A chirp 1 2 signal with carrier frequency f , chirp rate K , pulse width T is c r p transmitted and received. The transmitted signal can be written as (cid:181) (cid:182) tˆ (cid:169) (cid:163) (cid:164)(cid:170) S(m,n,t) = rect exp j 2πf t+πK tˆ2 , (1) c r T p where m is the cross-track dimensional equivalent array element number, n is the along-track dimensional pulse number and tˆis fast time. Here we are going to develop an equation for the signal phase received by DLSLA 3D SAR system from a single scatter object P at scene coordinate (x, y, z). This development assumes an ideal point scatter object with radar cross section σ whose amplitude and phase characteristics do not vary with frequency and aspect angle. For simplicity, the receiving signal model ignores antenna gain, amplitude Progress In Electromagnetics Research, Vol. 129, 2012 293 effects of propagation on the signal and any additional time delays due to atmospheric effects. The signal received from target P at array element number m and pulse number n is (cid:181) (cid:182) tˆ−t (cid:110) (cid:104) (cid:161) (cid:162) (cid:105)(cid:111) S (m,n,t) = a rect d exp j 2πf (t−t )+πK tˆ−t 2 , (2) r t c d r d T p √ where a = σ, t = 2Rt is the dual time delay from APC to target t d c P. It is appropriate to view the received signal as a three dimensional signal in the coordinates m, n and tˆ. The reference receiving signal from APC to reference point at SC, O, is (cid:110) (cid:104) (cid:105)(cid:111) (cid:161) (cid:162) S (m,n,t) = exp j 2πf (t−t )+πK tˆ−t 2 , (3) ref c 0 r 0 where t = 2R0 is the dual time delay from APC to reference point at 0 c scenecenterO. Theradarreceivingsignalisthevideofrequencysignal generated by mixing the received signal from target P with reference received signal from reference point O. It is convenient to write the video frequency echo signal in the form S (m,n,t) = S (m,n,t)×S∗ (m,n,t) if r ref (cid:181) (cid:182) tˆ−t (cid:169) (cid:161) (cid:162)(cid:170) = a rect d exp jΦ m,n,tˆ . (4) t T p The phase term Φ(m,n,tˆ) in Eq. (4) can be written as, (cid:161) (cid:162) 4πK Φ m,n,tˆ = −KR+ rR2, (5) c2 where R = R −R . t 0 (cid:169) (cid:170) InEq.(5), K = 4πKr(fc +tˆ−2R0)andexp j4πKrR2 istheResidual c Kr c c2 VideoPhase(RVP)term. TheRVPtermistheconsequenceofdechirp- on-receiveapproach. Itcanbecompletelyremovedfromtheradarecho signalbyapreprocessingoperationwhichisillustratedinAppendixA. 2.3. Platform Motion Effects During Time-divided Transmitting-receiving Procedure Now, we run back over the echo data collection procedure. The sparselineararraymaintainsmultipletransmittingarrayelementsand multiple receiving array elements and the sparse array works in time- divided mode. As shown in Fig. 2, There are L transmitting array 1 elements and L receiving array elements. The transmitting array 2 elements work sequentially, array element T transmits signal first, 1 294 Peng et al. all the receiving array elements receive echo signal simultaneously and L equivalent array elements are obtained, the time interval of 2 this transmitting and receiving procedure is ∆T. Then array element T transmits signal, and the procedure loops until array element 2 T finishes transmitting signal. Finally the radar obtains L L L1 1 2 equivalent array elements. The cross-track dimensional equivalent element number m ∈ [1,L L ] in the video frequency echo signal 1 2 S (m,n,t). The echo data collection procedure is shown in Fig. 2. if As the radar works in the time division mode, the whole transmitting and receiving procedure of the sparse linear array is not short enough to be treated under stop-and-go approximation. The platform motion during the time-divided transmitting-receiving procedure should be considered in the video frequency echo signal S (m,n,t). As shown in Fig. 3, O is the SC with coordinates if (0,0,0). The platform movement is ∆x = l×V×∆T, where l means array element T transmits signal, V is flying platform velocity. The l coordinates of APC is (x +∆x,y ,z ). The coordinates of target P is i i i (x, y, z). The range between scene center and APC is R , the range 0 between target and APC is R , where t (cid:113) R = (x +∆x)2+r2 0 i 0 (cid:113) R = [(x +∆x)−x]2+r2, (6) t i t (cid:113) (cid:112) and where r = y2+z2, r = (y−y )2+(z−z )2, then making 0 i i t i i the substitution (without considering RVP) and taking along-track Figure 2. Time divided MIMO echo data collection. Progress In Electromagnetics Research, Vol. 129, 2012 295 (a) (b) Figure 3. Reference and instantaneous range with platform motion. (a) Reference range. (b) Instantaneous range. dimensional Fourier transform, we obtain (cid:90) F {S (m,n,t)} = S (m,n,t)exp{−jK x }dx x if if x i i (cid:90) x = A exp{−jKR}exp{−jK x }dx 0 x i i x (cid:110) (cid:112) (cid:111) = A exp{jK ∆x}exp j K2−K2(r −r )−jK x , (7) 1 x x t 0 x where (cid:113) (cid:113) R = R −R = r2+[(x +∆x)−x]2− r2+(x +∆x). t 0 t i 0 i and A and A are the corresponding amplitudes before and after 0 1 Fourier transform. Equation (7) describes the platform motion during the time- divided transmitting-receiving procedure causes an along-track dimensional phase term, this phase term should be compensated in image reconstruction procedure. 296 Peng et al. 2.4. Echo Signal Bandwidth Before and after Dechirp Processing As shown in Fig. 4, the AD sampling begins at the same time as the dechirp starts. The video signal is digitized by the AD converter. Corresponding to the return of the center of the pulse from the quantized reference range R , the AD trigger for the first sample steps 0 back by T /2 to the start of the pulse and steps back T /2 to the start p s of the range window. The total sampling time is T +T . p s TheADconvertersamplesthevideosignalatfasttimetˆ(m,n,k): s k 2R T T tˆ(m,n,k) = + 0 − p − s. (8) s F c 2 2 s Eq. (8) describes the value of a sample of the array element m, along- track pulse n and wave-propagation dimension sample k. F is the s sampling frequency, T = 2R = 4∆R is the time interval associated s c c with the wave-propagation dimensional range gate R(R∈[−∆R,∆R]). The phase of the signal from APC sample cell (m, n, k) to target P is (cid:181) (cid:182) 4πK f k T T 4πK Φ(m,n,k) = − r c + − p − s R+ rR2. (9) c K F 2 2 c2 r s Figure 4. Timing diagram for AD samples.