Airborne Downward Looking Sparse Linear Array 3-D SAR Heterogeneous Parallel Simulation

Abstract

The airborne downward looking sparse linear array three dimensional synthetic aperture radar (DLSLA 3-D SAR) operates nadir observation with the along-track synthetic aperture formulated by platform movement and the cross-track synthetic aperture formulated by physical sparse linear array. Considering the lack of DLSLA 3-D SAR data in the current preliminary study stage, it is very important and essential to develop DLSLA 3-D SAR simulation (echo generation simulation and image reconstruction simulation, including point targets simulation and 3-D distributed scene simulation). In this paper, DLSLA 3-D SAR imaging geometry, the echo signal model and the heterogeneous parallel technique are discussed first. Then, heterogeneous parallel echo generation simulation with time domain correlation and the frequency domain correlation method is described. In the following, heterogeneous parallel image reconstruction simulation with two imaging algorithms, e.g., 3-D polar format algorithm, polar formatting and L1 regularization algorithm is discussed. Finally, the point targets and the 3-D distributed scene simulation are demonstrated to validate the effectiveness and performance of our proposed heterogeneous parallel simulation technique. The 3-D distributed scene employs airborne X-band DEM and P-band Circular SAR image of the same area as simulation scene input.

1. Introduction

The airborne downward looking sparse linear array three dimensional synthetic aperture radar (DLSLA 3-D SAR) can be placed on small and mobile platforms to acquire high resolution full 3-D microwave images. The 3-D resolution is acquired by wave propagation dimensional pulse compression with wide band chirp signal, along-track aperture synthesis with flying platform movement, and cross-track aperture synthesis with physical sparse linear array [1–5]. DLSLA 3-D SAR observes the nadir areas of the platform, which means it overcomes the restrictions of shading and lay over effects by trees, buildings, and special terrain shapes in general side-looking SAR [2,3]. DLSLA 3-D SAR can be developed for various applications, such as city planning, environmental monitoring, Digital Surface Model (DSM) generation, disaster relief, surveillance and reconnaissance, etc. [6–12].

Downward Looking Imaging Radar (DLIR) was first introduced by Gierull in 1999 [13]. He took advantage of aperture synthesis by along-track platform movement and cross-track physical linear array to generate along-track and cross-track two dimensional images with a single frequency transmitted signal. Nouvel et al. in ONERA made use of the chirp transmitted signal instead of the single frequency signal, and developed the concept of downward looking three dimensional synthetic aperture radar (DL 3-D SAR) [14–17]. Klare et al., in FGAN FHR, designed an Airborne Radar for Three-dimensional Imaging and Nadir Observation (ARTINO) [18–20]. Recently, DLSLA 3-D SAR image reconstruction algorithms are widely researched [21–30], especially for cross-track sparse reconstruction. However, the DLSLA 3-D SAR simulation technique is insufficient, especially the DLSLA 3-D SAR echo generation simulation and 3-D distributed scene image reconstruction simulation.

In order to assess the key technology and validate the imaging algorithms, it is very important and essential to develop DLSLA 3-D SAR simulation. DLSLA 3-D SAR echo generation simulation with airborne platform motion error can be used for DLSLA 3-D SAR motion compensation research. DLSLA 3-D SAR echo generation simulation with array channel amplitude and phase inconsistency error can be used for DLSLA 3-D SAR array channel amplitude and phase inconsistency compensation research. Furthermore, the effectiveness and performance of the DLSLA 3-D SAR imaging algorithms can be demonstrated with image reconstruction simulation.

Considering the heavy computation complexity and computation load in 3-D circumstances, heterogeneous parallel simulation technique for DLSLA 3-D SAR simulation is demonstrated in this paper [31]. The heterogeneous parallel time domain correlation method and frequency domain correlation method for the DLSLA 3-D SAR echo generation is illustrated. The time domain correlation method for echo generation is precise and easy, while being quite time consuming, which is suitable for point targets simulation. The frequency domain correlation method for DLSLA 3-D SAR echo generation is computation efficient, which is very suitable for 3-D distributed scene simulation [32,33]. The heterogeneous parallel DLSLA 3-D SAR image reconstruction simulation, with the 3-D polar format algorithm [34], polar formatting and an L1 regularization algorithm [35], is performed with wave front curvature phase error compensation. The heterogeneous parallel DLSLA 3-D SAR image reconstruction simulation can be used for uniform and non-uniform sparse cross-track imaging which possesses the advantage of low computation (such as the polar formatting operation) and high precision (wave front curvature phase error compensation) [34,35].

The remainder of the paper is organized as follows. In Section 2, the DLSLA 3-D SAR imaging geometry, the echo signal model and the heterogeneous parallel technique are illustrated. In Section 3, DLSLA 3-D SAR heterogeneous parallel echo generation simulation with time domain correlation and frequency domain correlation method is discussed. In Section 4, DLSLA 3-D SAR heterogeneous parallel image reconstruction simulation with 3-D polar format algorithm, polar formatting and L1 regularization algorithm is described. The effectiveness and performance of the proposed DLSLA 3-D SAR heterogeneous simulation is demonstrated with point targets and 3-D distributed scene simulation in Section 5. Finally, the conclusion and the current focus of our research are outlined in Section 6.

2. DLSLA 3-D SAR Imaging Geometry, Echo Signal Model and Heterogeneous Parallel Technique

2.1. DLSLA 3-D SAR Imaging Geometry

Airborne DLSLA 3-D SAR operates nadir observation and its working geometry is shown in Figure 1. X axis is the along-track axis, Y axis is the cross-track axis, Z axis is from top to bottom perpendicular to the XY plane, O is the origin. Q is the sparse linear array equivalent virtual Antenna Phase Center (APC). P is the target in the imaging scene, P′ is the projection point of target P onto YZ plane. APC flight path QX′ is parallel to X axis. QY′ is parallel to Y axis. OP is the reference range from the coordinate center to the target P with the distance of ρ. QP is the instantaneous range from the virtual APC to the target P with the distance of ρ′. γ₁ is the angle between QP and QX′ called the along-track Doppler cone angle. γ₂ is the angle between QP and QY′ called the cross-track Doppler cone angle. ϕ is the angle between OP and OP′, θ is the angle between OP′ and OZ. The coordinates of the target P in Cartesian coordinate can be written as (ρsinϕ, ρcosϕsinθ, ρcosϕcosθ).

For conventional side-looking 2-D SAR echo signal acquisition, only the along-track Doppler cone angle γ₁ is obtained with platform movement, which means only along-track resolution is obtained. The targets with different cross-track Doppler cone angle located in the same range and along-track bin cannot be distinguished in 2-D SAR. For airborne DLSLA 3-D SAR echo signal acquisition, along-track Doppler cone angle γ₁ is obtained with platform movement and cross-track Doppler cone angle γ₂ is obtained with sparse array, which means the scatters can be differentiated in the along-track and the cross-track dimension. Shadowing and lay over effects cannot be avoided in the side-looking 2-D SAR as the grazing angle is small, especially in the far range gate units. Airborne DLSLA 3-D SAR owns large grazing angle as it operates nadir observation, therefore, shadowing and lay over effects can be overcome [2].

2.2. DLSLA 3-D SAR Echo Signal Model

A chirp signal with carrier frequency f_c, chirp rate K_r, pulse width T_P is transmitted and received. The transmitted signal can be written as [1]

$$S(m,n,t)=\text{rect}(\frac{t}{T_p})\times \text{exp}(\text{j}2\pi f_{c}t+\text{j}\pi K_t{t}^2)$$

(1)

where m is the along-track sample number, n is the cross-track sample number and t is the fast time. Here we are going to develop an equation for the signal received by the airborne DLSLA 3-D SAR for the scatter object P at scene coordinate (ρsinϕ, ρcosϕsinθ, ρcosϕcosθ). This development assumes an ideal point scatter object with radar cross section σ whose amplitude and phase characters do not vary with frequency and aspect angle [1].

The APC coordinate is (x_m, y_n, 0). The radar recorded signal is the video frequency signal generated by mixing the received signal with the carrier frequency signal. For simplicity, the receiving signal model ignores the antenna gain, amplitude effects of propagation on the signal and any additional time delays due to the atmospheric effects [1]. It is convenient to write the video frequency echo signal received from target P in the form

$$S(m,n,t)=a_t\times \text{rect}(\frac{t}{T_p})\times \text{exp}(-\text{j}2\pi f_c{t}_d+\text{j}\pi K_r{(t-t_d)}^2)$$

(2)

where $a_t=\sqrt{\sigma }$, $t_d=\frac{2{\rho }^{\prime }}{c}$ is the dual time delay from APC to target P.

2.3. Heterogeneous Parallel Technique

In this paper, we apply a heterogeneous parallel technique on the Central Processing Unit (CPU) [36] and Graphical Processing Unit (GPU) [37] platform for the DLSLA 3-D SAR simulation. The CPU is mainly used for logical control, while the GPU is mainly used for large-scale computation. We choose GPU for large-scale computation for two reasons: (1) GPU is specialized for compute intensive, highly parallel computation as more transistors are devoted to data processing rather than data caching and flow control; (2) GPU is especially well suited to address problems that can be expressed as data parallel computations (the same program is executed on many data elements in parallel) with high arithmetic intensity. Because the same program is executed for each data element, there is a lower requirement for sophisticated flow control, and because it is executed on many data elements and has high arithmetic intensity, the memory access latency can be hidden with calculations instead of big data caches [38]. The comparison of floating-point operations per second and memory bandwidth for the CPU and GPU [38] is shown in Figure 2.

Heterogeneous parallel model is shown in Figure 3, the serial code, which cannot run in parallel runs on the CPU. The parallel code runs on the GPU. The whole program contains several serial code parts and parallel code parts. The serial code and parallel code are connected with I/O. In the parallel code, computation grid is set up with multiple computation blocks. Every computation block contains multiple threads. All the threads in the computation block run in parallel. All the computation blocks in the computation grid run in parallel. The optimized setup of computation grid, block and thread can be referenced according to NVIDIA CUDA C Programming Guide [38] and NVIDIA CUDA C Best Practices Guide [39].

3. DLSLA 3-D SAR Heterogeneous Parallel Echo Generation Simulation

3.1. Heterogeneous Parallel Echo Generation Simulation with Time Domain Correlation Method

For an imaging scene, the echo signal can be generated according to Equation (2) called time domain correlation method. For every target, we compute the range from the target to the APC, then the echo signal of the target is recorded by the A/D device for sample depth f_s · T_P (f_s is A/D device sampling rate) with dual time delay t_d from the target P to the APC. The time domain correlation echo generation method is precise, while quite time consuming, as the sample depth needs to be judged for every target.

Heterogeneous parallel echo generation simulation with the time domain correlation method is shown in Figure 4. POS data, array amplitude and phase error data, 3-D DEM data and radar cross Section (RCS) data of the imaging scene are the input of the simulation. The echo generation procedure is operated with the heterogeneous parallel technique on the CPU and GPU platform. We offer every target in the 3-D distributed imaging scene with a computation block in the GPU, which means the block number in GPU computation grid is the target number. The echo generation time for all the targets with heterogeneous parallel computation is almost the same as echo generation time for a single target as all the computation blocks run in parallel. The computation block can be offered with parallel threads to accelerate the echo generation of the target (e.g., the parallel threads can be set up to 1024 for a GPU device with computation capability 2.x or 3.x). The threads in the computation block compute the echo signal from APC to the target according to Equation (2). After the echo signal from all the targets in the imaging scene to the APC is computed, we perform a synchronization among all the computation blocks to ensure all the echo data is effective. Then the echo signal simulation is finished after sending the echo data from CPU to GPU and writing the echo data to the hard disk.

3.2. Heterogeneous Parallel Echo Generation Simulation with Frequency Domain Correlation Method

For 3-D distributed imaging scene, it will be very time consuming to generate the DLSLA 3-D SAR echo signal with time domain correlation method. This is caused by the sample depth computation and judgment for so many targets in the imaging scene. We perform wave propagation dimensional FFT to Equation (2) with variable t. Then the echo signal is modeled in wave propagation frequency domain [32,33].

$$S(m,n,f_k)=a_t\times \text{exp}(-\text{j}\frac{4\pi \left(f_c+f_k\right)}{c}{\rho }^{\prime })\times S(f_k)$$

(3)

where, f_c is the carrier frequency, f_k is the baseband frequency, $S(f_k)=\text{exp}(-\text{j}\frac{\pi f_k^2}{{\kappa }_r})\times \text{exp}(\text{j}\frac{\pi }{4})$ is the FFT of the chirp transmitted signal. For echo generation, we can compute the frequency domain phase delay $\text{exp}(-\text{j}\frac{4\pi (f_c+f_k)}{c}{\rho }^{\prime })$ of every target, then multiply the phase delay with the frequency chirp transmitted signal S(f_k), the echo signal can be transformed from S(m, n, f_k) to S(m, n, t) with wave propagation dimensional IFFT, finally, rectangular windowing function is applied to S(m, n, t) to finish range gating for all targets. This is frequency domain correlation echo generation method. This method is computation efficient as no sample depth computation and judgment is required during echo generation for every target.

The heterogeneous parallel echo generation simulation with frequency domain correlation method is shown in Figure 5. POS data, array amplitude and phase error data, 3-D DEM data and radar cross section data of the imaging scene is the input of the simulation. The echo generation procedure is operated with the heterogeneous parallel technique on the CPU and GPU platform. We offer every target in the 3-D distributed imaging scene a computation block in GPU, which means the block number in GPU computation grid is the target number. The echo generation time for all the targets with heterogeneous parallel computation is almost the same as echo generation time for a single target as all the computation blocks run in parallel. The computation block can be offered parallel threads to accelerate the echo generation of the target (e.g., the parallel threads can be set up to 1024 for a GPU device with computation capability 2.x or 3.x). The FFT of the chirp transmitted signal is operated with the GPU first. Then the multiple threads in the computation block compute the phase delay from the target to the APC. The multiplication between the phase delay and the frequency domain chirp transmitted signal is operated with multiple threads in computation block. IFFT is applied to the product of the phase delay and the frequency domain chirp transmitted signal. Finally, a wave propagation windowing operation is performed to keep the time domain echo signal in the wave propagation sample range gate. Then the echo signal simulation is finished after sending the echo data from CPU to GPU and writing the echo data to the hard disk.

3.3. Heterogeneous Parallel Echo Generation Simulation Applicability

Heterogeneous parallel echo generation simulation with time domain correlation method is time consuming while quite precise. Heterogeneous parallel echo generation simulation with frequency domain correlation method eliminates sample depth computation for every target. There exists truncated error in the frequency domain correlation method caused by use of the windowing operation with finite window length. We recommend point targets simulation to be used for imaging index (e.g., PSLR and ISLR) computation, with the heterogeneous parallel time domain correlation echo generation method, although the truncated error with heterogeneous parallel frequency domain correlation is very small and can be ignored in most instances. For large scale 3-D distributed imaging scene simulation, heterogeneous parallel frequency domain correlation echo generation method is the best choice [32,33].

4. DLSLA 3-D SAR Heterogeneous Parallel Image Reconstruction Simulation

After wave propagation dimensional frequency domain matched filtering, the signal can be written as

$$S(m,n,f_k)=a_t\times \text{exp}(-\text{j}\frac{4\pi (f_c+f_k)}{c}\sqrt{{\left(x_m-\rho \hspace{0.17em}\text{sin}\hspace{0.17em}\varphi \right)}^2+{\left(y_n-\rho \hspace{0.17em}\text{cos}\hspace{0.17em}\varphi \hspace{0.17em}\text{sin}\hspace{0.17em}\theta \right)}^2+{\left(\rho \hspace{0.17em}\text{cos}\hspace{0.17em}\varphi \hspace{0.17em}\text{cos}\hspace{0.17em}\theta \right)}^2})$$

(4)

The signal in Equation (4) can be rewritten in the form

$$S(m,n,f_k)=S_P(m,n,f_k)\times S_E(m,n,f_k)$$

(5)

where $S_P(m,n,f_k)=a_t\times \text{exp}(-j\frac{4\pi (f_c+f_k)}{c}(\rho -\mathit{sin}\varphi x_m-\mathit{cos}\varphi \mathit{sin}\theta y_n))$ is the basic imaging term which is used for DLSLA 3-D SAR imaging, $S_E(m,n,f_k)=\text{exp}(-j\frac{4\pi (f_c+f_k)}{c}({\sum }_i{\sum }_{j}O(x_m^i\cdot y_n^j)))$ is the wave front curvature phase error term that causes geometry distortion and defocus the reconstructed DLSLA 3-D SAR image. According to Equation (5), the basic imaging term can be written as

$$S_P(m,n,f_k)=\frac{S(m,n,f_k)}{S_E(m,n,f_k)}=S(m,n,f_k)\times S_E^H\left(m,n,f_k\right)$$

(6)

where is $S_E^H(m,n,f_k)$ the complex conjugate of the wave front curvature phase error term, which can be generated with POS data. The radar reflectivity a_t at target P can be obtained from the basic imaging term S_P (m, n, f_k) via

$$a_t=\sum _{m=0}^{\mathit{Na}-1}\sum _{n=0}^{\mathit{Ne}-1}\sum _{k=0}^{\mathit{Nr}-1}{S}_P(m,n,f_k)\times \text{exp}(\text{j}\frac{4\pi \left(f_c+f_k\right)}{c}\left(\rho -\text{sin}\hspace{0.17em}\varphi \hspace{0.17em}{x}_m-\text{cos}\hspace{0.17em}\varphi \hspace{0.17em}\text{sin}\hspace{0.17em}\theta \hspace{0.17em}{y}_n\right))$$

(7)

4.1. DLSLA 3-D SAR Heterogeneous Parallel Image Reconstruction Simulation with 3-D Polar Format Algorithm

The basic imaging term can be polar formatted in the along-track dimension and cross-track dimension from sample x_m, y_n to x′_m, y′_n with the relationship x_m(f_c + f_k) = x′_mf_c and y_n (f_c + f_k) = y′_nf_c. Then the radar reflectivity at target P in Equation (7) can be obtained with the endomorphism mapping principle of the polar formatting [34,40,41]

$$\begin{array}{ll}{a}_t\hfill & =\sum _{m=0}^{\mathit{Na}-1}\sum _{n=0}^{\mathit{Ne}-1}\sum _{k=0}^{\mathit{Nr}-1}\mathbf{P}[S_P(m,n,f_k)]\times \text{exp}(\text{j}2\pi ((f_c+f_k)\alpha -x_m^{\prime }\beta -y_n^{\prime }\gamma ))\hfill \\ \hfill & =\sum _{m=0}^{\mathit{Na}-1}\sum _{n=0}^{\mathit{Ne}-1}\sum _{k=0}^{\mathit{Nr}-1}\mathbf{P}[S(m,n,f_k)\times S_E^H(m,n,f_k)]\times \text{exp}(\text{j}2\pi ((f_c+f_k)\alpha -x_m^{\prime }\beta -y_n^{\prime }\gamma ))\hfill \\ \hfill & =\sum _{m=0}^{\mathit{Na}-1}\sum _{n=0}^{\mathit{Ne}-1}\sum _{k=0}^{\mathit{Nr}-1}\mathbf{P}[S(m,n,f_k)]\times \mathbf{P}[S_E^H(m,n,f_k)]\times \text{exp}(\text{j}2\pi ((f_c+f_k)\alpha -x_m^{\prime }\beta -y_n^{\prime }\gamma ))\hfill \end{array}$$

(8)

where P means along-track dimension and cross-track dimension polar formatting. $\alpha =\frac{2\rho }{c}$ is directly proportional to range ρ, $\beta =\frac{2\mathit{sin}\varphi }{{\lambda }_c}$ is a sinusoidal function of the polar angle ϕ, $\gamma =\frac{2\mathit{cos}\varphi \mathit{sin}\theta }{{\lambda }_c}$ is a cosine and sinusoidal function of polar angle ϕ and θ, N_a is the along-track sample number, N_e is the cross-track sample number. The radar reflectivity at target P can be obtained in (β, γ, α) coordinate from Equation (8) with wave preoperational IFFT, along-track and cross-track FFT. The reconstructed image can be written as

$$a_t=\text{sinc}(B_r\alpha )\times \text{sinc}(L_a\beta )\times \text{sinc}(L_e\gamma )$$

(9)

where L_a is the length of the along-track aperture, L_e is the length of the cross-track aperture. The (β, γ, α) coordinate image can be interpolated to Cartesian coordinate with the relationship [34]

$$x=\frac{c{\lambda }_c}{4}\alpha \beta ,y=\frac{c{\lambda }_c}{4}\alpha \gamma ,z=\frac{c{\lambda }_c}{4}\alpha \sqrt{1-{\beta }^2-{\gamma }^2}$$

(10)

This algorithm should be applied to uniform along-track and cross-track sample circumstances. The flow diagram of the heterogeneous parallel image reconstruction simulation with 3-D polar format algorithm is shown in Figure 6. The wave front curvature phase error compensation term is generated with POS data on the GPU. After wave propagation matched filtering, the echo signal and wave front curvature phase error compensation terms are polar formatted in along-track and cross-track dimensions with the GPU. Then, the multiplication between the polar formatted echo signal and the wave front curvature phase error compensation term is operated on the GPU. Finally, the FFT/IFFT of the product is operated on the GPU. After transferred from the GPU to CPU, the imaging result can be written on the hard disk or displayed.

4.2. DLSLA 3-D SAR Heterogeneous Parallel Image Reconstruction Simulation with Polar Formatting and L1 Regularization Algorithm

In the real system, the cross-track phase centers are usually non-uniformly and sparsely distributed due to the array element installation positions restricted by the airborne platform and the airborne wing tremor effect. Then the cross-track FFT in Equation (8) cannot be operated. In this circumstance, we perform wave propagation IFFT and along-track FFT, and obtain the wave propagation and along-track imaging result in the form [35]

$$a_t=\text{sinc}(B_r\alpha )\times \text{sinc}(L_a\beta )\times \sum _{n=0}^{\mathit{Ne}-1}{a}_t\times \text{exp}(\text{j}2\pi y_n^{\prime }\gamma )$$

(11)

Equation (11) can be written in the matrix form

$${\mathbf{U}}_{\mathbf{n}}=\mathbf{R}\times \mathbf{A}+\mathbf{\xi }$$

(12)

where U_n is the measurement vector of length N_e, R with R_n×l = exp (j2πy′_nγ_l) is an N_e × L mapping (Fourier transform) matrix, and A is uniformly sampled in γ dimension at γ_l(l=1,2 ... L.). ξ is the noise, e.g., non-scene-center residual wave front curvature error and thermal noise. In ordinary circumstances, N_e ≪ L, which means the system model in Equation (12) is very undetermined. Hence, the solutions of A are infinite. As described in the introduction and as outlined in detail in [24], for DLSLA 3-D SAR, A is sparse and assumed to contain no more than 10 non-zero elements (positions, amplitudes and phases) in cross-track domain with the wave propagation and along-track pixel (α, β) [25,26]. The α and β dimension focused image with the model in Equation (12) satisfies the RIP feature. The Compressed Sensing (CS) theory tells that, if the mapping matrix R fulfills the RIP and incoherence properties and $N_e=O(K\hspace{0.17em}\mathit{log}(\frac{L}{K}))$, the convex L1 minimization give the resolution of the undetermined model in Equation (12) [25,26]. In case there is no prior sparsity about K and in the presence of measurement noise, it can be approximated by the following convex L1 minimization [26].

$$\widehat{\gamma }=\text{arg}\hspace{0.17em}\underset{\gamma }{\text{min}}\{{\Vert {\mathbf{U}}_{\mathbf{n}}-\mathbf{R}\times \mathbf{A}\Vert }_2^2+{\lambda }_k{\Vert \gamma \Vert }_1\}$$

(13)

where λ_k is the regularization parameter as a function of N_e and ξ. The convex L1 regularization is robust to phase noise and does not require model selection to provide the prior knowledge about the number of scatters K [26].

The polar formatting and L1 regularization algorithm can be applied to non-uniform and sparse cross-track sample circumstances. The flow diagram of the heterogeneous parallel image reconstruction simulation with polar formatting and L1 regularization algorithm is shown in Figure 7. The wave front curvature phase error compensation term can be generated with POS data. After wave propagation matched filtering, the echo signal and wave front curvature phase error compensation term are polar formatted in along-track and cross-track dimension with GPU. Then, the multiplication between the polar formatted echo signal and the wave front curvature phase error compensation term is operated on the GPU. Wave propagation and along-track dimension imaging is operated with IFFT/FFT on GPU. Finally, the cross-track imaging with L1 regularization is operated on the CPU (we are working on the GPU version). The imaging result is processed with CPU, which can be written on the hard disk or displayed.

5. Simulation Results

In this section, we present two DLSLA 3-D SAR numerical simulation experiments to illustrate the performance of our proposed heterogeneous parallel simulation method. The sparse linear array we used is composed of 8 transmitting array elements and 32 receiving array elements [34,35]. The configuration of the sparse linear array is shown in Figure 8. The realistic array configuration is on the top of Figure 8, while the equivalent virtual APC distribution is at the bottom of Figure 8. The coordinates of transmitting array elements, receiving elements, equivalent APC are shown in

Table 1. Cross-track coordinates of array elements.

Coordinates	Value (m)
Transmitted Array	Ti = −1.34 + i × 0.02, i = 1, 2,…,4.
Transmitted Array	Ti = −1.14 + i × 0.02, i = 5, 6,…,8.
Receiving Array	Ri = −1.32 + i × 0.08, i = 1, 2,…,32.
Equivalent Array	TRi = −1.28 + (i − 1) × 0.01, i = 1, 2,…,256.

.

5.1. Point Targets Heterogeneous Parallel Simulation

The point targets simulation is used to illustrate the capability of the reconstructed image on different height with the proposed imaging algorithm. The measurement parameters used in point targets simulation are shown in

Table 2. Measurement parameters used in the simulation.

Parameters	Value
Center Frequency	37.5 GHz
Transmitting Signal Bandwidth	300 MHz
A/D Sampling Frequency	360 MHz
Platform Fly Height	1,000 m
Platform Fly Velocity	50 m/s
Transmitting Signal Pulse Width	1.0 us
A/D Sampling Range Gate	[750.0, 1026.7 m]
Range Sample Number	1024
PRF (Pulse Repetition Frequency)	5,000 Hz
Along-track Dimension Sampling Interval	0.01 m
Along-track Dimension Sampling Number	256
Transmitting Array Elements	8
Receiving Array Elements	32
Beam width of T/R Array	14 × 14
Equivalent Phase Center Number	256
Cross-track Dimension Sampling Interval	0.01 m

. 3 × 8 point targets with identical radar cross Section are distributed on three circles at different height. A sketch of the point targets distribution in the 3-D space is shown in Figure 9a. The radius of the three circles from bottom to top is 150 m, 100 m, and 50 m, respectively. The Z axis coordinate value of the three circles from bottom to top is 980 m, 880 m, and 780 m, respectively. The projection of the three circles onto XY plane is shown in Figure 9b. The projection of the three circles onto XZ plane is shown in Figure 9c. The projection of three circles onto YZ plane is shown in Figure 9d.

The heterogeneous parallel echo generation simulation with frequency domain correlation method possesses ringing effect after wave propagation dimensional matched filtering. This may reduce the echo generation precision. We compare the wave propagation dimensional matched filtering profile and the imaging result profile with the two echo generation methods in Figure 10. The ringing effect caused by wave propagation windowing during echo generation is obvious in Figure 10a. However, the wave propagation imaging result with the echo signal generated with the two methods seems almost the same, which is obvious in Figure 10b.

The image reconstructed by 3-D PFA without wave front curvature error compensation is shown in Figure 11. Figure 11a is the 3-D imaging result in (β, γ, α) coordinates without wave front curvature error compensation. Figure 11b is the projection of the 3-D imaging result onto βγ plane. Figure 11c is the projection of the 3-D imaging result onto βα plane. Figure 11d is the projection of the 3-D imaging result onto γα plane.

By comparing Figure 11 with Figure 12, the reconstructed image with wave front curvature error compensation shows better focusing performance. In order to illustrate the imaging results with and without wave front curvature error compensation specifically, the β, γ, α dimensional profile of the imaging result is shown in Figure 13a, Figure 13b, Figure 13c, respectively. It is obvious that a sharper main lobe and lower side lobe is obtained after wave front curvature error compensation. A sharper main lobe brings better resolution in the imaging result, while a lower side lobe brings high contrast ratio in the reconstructed image.

After wave front curvature error compensation, the reconstructed image is in polar coordinate grid, we need to transform it to Cartesian coordinate grid according to Equation (10). The 3-D imaging results in Cartesian coordinate grid, as shown in Figure 14a. The projection of the 3-D imaging result onto the XY plane (orthographic projection) is shown in Figure 14b. The projection of the 3-D imaging result onto the XZ plane is shown in Figure 14c. The projection of the 3-D imaging result onto the YZ plane is shown in Figure 14d. In order to demonstrate the performance of the proposed imaging algorithm, the Peak Side Lobe Rates (PSLR) and the Integrated Side Lobe Rates (ISLR) are all evaluated, as indicated in

Table 3. The measured Peak Side Lobe Rates (PSLR) and Integrated Side Lobe Rates (ISLR) without wave front curvature error compensation.

Measured Parameter	Along-Track	Cross-Track	Wave-Propagation
PSLR (dB)	−9.3	−9.4	−16.7
ISLR (dB)	−10.5	−10.3	−12.7

and

Table 4. The measured PSLR and ISLR with wave front curvature error compensation.

Measured Parameter	Along-Track	Cross-Track	Wave-Propagation
PSLR (dB)	−13.3	−13.3	−13.4
ISLR (dB)	−14.9	−14.3	−12.7

. The PSLR is the ratio between the height of the largest side lobe and the height of the main lobe. The ISLR is often used to analyze the side lobe power of point spread function (PSF).

5.2. 3-D Distributed Scene Heterogeneous Parallel Simulation

In order to make the numerical simulation more precise and consistent with the real 3-D distributed imaging scene, we choose an airborne X-band DEM surface and airborne P-band Circular SAR image (the image is formulated on the DEM surface with Back-projection algorithm) of the same area for distributed scene simulation [42]. The DEM distributes in a region of 400 m (X dimension) × 400 m (Y dimension) × 40 m (Z dimension). The coordinate of scatters for the simulation is (x, y, z): x is uniformly distributed from [−200 m, 200 m] with 2 m interval, y is uniformly distributed from [−200 m, 200 m] with 2 m interval, z is the elevation value of the DEM at(x,y). The P-band Circular SAR image is used as a radar cross section of the distributed imaging scene. The surface figure of the DEM region is shown in Figure 15a. The airborne P-band Circular SAR image is shown in Figure 15b. The measurement parameters used in three dimensional distributed surface targets simulation is shown in

Table 5. Measurement parameters used in the simulation.

Parameters	Value
Center Frequency	37.5 GHz
Transmitting Signal Bandwidth	300 MHz
A/D Sampling Frequency	360 MHz
Platform Fly Height	1,000 m
Platform Fly Velocity	50 m/s
Transmitting Signal Pulse Width	3.85 us
A/D Sampling Range Gate	[974.0, 1026.5 m]
Range Sample Number	1024
PRF (Pulse Repetition Frequency)	5,000 Hz
Along-track Dimension Sampling Interval	0.01 m
Along-track Dimension Sampling Number	256
Transmitting Array Elements	8
Receiving Array Elements	32
Beam width of T/R Array	14 × 14
Equivalent Phase Center Number	256
Cross-track Dimension Sampling Interval	0.01 m

. In order to accelerate the echo generation procedure, frequency correlation method is used for three dimensional echo generation [32,33]. The elapsed-time of echo generation with time domain correlation method with the parameters in Table 5 by using heterogeneous parallel technique is about 190 h. The elapsed-time of echo generation with frequency domain correlation method with the parameters in Table 5 by using heterogeneous parallel technique is about 5 h. It is obvious that the frequency domain correlation method possesses much higher computation efficiency.

The 3-D reconstructed image in polar coordinates without wave front curvature error compensation is shown in Figure 16a. The projection of the 3-D reconstructed image in polar coordinates without wave front curvature error compensation onto βγ plane is shown in Figure 16b. The 3-D reconstructed image in polar coordinates with wave front curvature error compensation is shown in Figure 16c. The projection of the 3-D reconstructed image in polar coordinates without wave front curvature error compensation onto βγ plane is shown in Figure 16d. The reconstructed image possesses better focus performance after wave front curvature error compensation. The distortion effect before and after wave front curvature error compensation is not very clear in visual.

The 3-D reconstructed image in Cartesian coordinates after wave front curvature error compensation is shown in Figure 17a. By comparison with DEM, three dimensional topography is well reconstructed with PFA. The projection of the 3-D reconstructed image in Cartesian coordinates onto XY plane is shown in Figure 17b. The Circular SAR image coincides with the orthographic projection of the 3-D reconstructed image in Cartesian coordinates.

According to the ARTINO sparse array configuration principle, the equivalent phase centers are uniformly distributed. In this simulation, we pick out 50% phase centers randomly in order to simulate non-uniform and sparse cross-track sampling circumstance. The reconstructed image in polar coordinates with 50% random cross-track virtual phase centers by using the proposed polar formatting and L1 regularization method is shown in Figure 18. Figure 18a is the 3-D reconstructed image in polar coordinate. Figure 18b is the orthographic projection of the 3-D polar image onto βγ plane. The image reconstructed with the proposed polar formatting and L1 regularization method possesses good focus performance. The 3-D reconstructed image in the polar coordinate grid is interpolated into a Cartesian coordinate grid with the relationship in Equation (9). Figure 18c is the 3-D reconstructed image in a Cartesian coordinate grid. Figure 18d is the orthographic projection of the 3-D Cartesian image.

We compare the reconstructed image with the proposed polar formatting and L1 regularization method tailored for non-uniform and sparse cross-track virtual phase centers and the reconstructed image with 3-D polar format algorithm tailored for uniform cross-track virtual phase centers. The comparison of the 3-D Cartesian image is shown in Figure 19. Figure 19a is the 3-D Cartesian image reconstructed with 50% random cross-track virtual phase centers by using the proposed polar formatting and the L1 regularization method. Figure 19b is the 3-D Cartesian image reconstructed with uniform cross-track virtual phase centers by using 3-D polar format algorithm. The measurement parameters for the two simulations are the same as is shown in Table 3. It can be seen in the box in Figure 19, the bridge pier and houses look much better in Figure 19a. The comparison of the orthographic projection of the 3-D Cartesian image is shown in Figure 20. Figure 20a is the orthographic projection of the 3-D Cartesian image reconstructed with 50% random cross-track virtual phase centers by using the proposed reconstruction method. Figure 20b is the orthographic projection of the 3-D Cartesian image reconstructed with uniform cross-track virtual phase centers by using 3-D polar format algorithm. As is shown in the box in Figure 20, the bridge pier, housing and farm areas look much better in Figure 20a. What's more, Figure 20a owns lower side lobe and better contrast than Figure 20b. The dynamic range of the orthographic projection reconstructed with the 3-D polar format algorithm is smaller than that of the polar formatting and L1 regularization method. A well reconstructed 3-D image can be developed for city planning [43,44], environmental monitoring [45], Digital Surface Model (DSM) generation [46], disaster relief [47], surveillance and reconnaissance [48], etc.

6. Conclusions

Downward looking sparse linear array three dimensional synthetic aperture radar (DLSLA 3-D SAR) imaging geometry and echo signal model are illustrated as the foundation of the simulation. Heterogeneous parallel technique is employed to accelerate the heavy computations in the echo generation simulation and the image reconstruction simulation. The time domain correlation method and frequency domain correlation method can be selected for heterogeneous parallel echo generation simulation, according to a computation and precision balance (time domain correlation method for point targets echo generation, frequency domain correlation method for large scale 3-D distributed imaging scene echo generation). Space domain image reconstruction methods, e.g., 3-D polar format algorithm, polar formatting and L1 regularization algorithm, are preferred for cross-track imaging to reduce the memory cost as the cross-track aperture is much shorter than the cross-track swath (3-D polar format algorithm for uniform cross-track virtual antenna phase center circumstance, polar formatting and L1 regularization algorithm for non-uniform cross-track virtual antenna phase center circumstance). The effectiveness and performance of our proposed heterogeneous parallel simulation technique are validated with the point targets scene and the 3-D distributed scene simulation.