Multi-Baseline Bistatic SAR Three-Dimensional Imaging Method Based on Phase Error Calibration Combining PGA and EB-ISOA

Abstract

Tomographic synthetic aperture radar (TomoSAR) is an advanced three-dimensional (3D) synthetic aperture radar (SAR) imaging technology that can obtain multiple SAR images through multi-track observations, thereby reconstructing the 3D spatial structure of targets. However, due to system limitations, the multi-baseline (MB) monostatic SAR (MonoSAR) encounters temporal decorrelation issues when observing the scene such as forests, affecting the accuracy of the 3D reconstruction. Additionally, during TomoSAR observations, the platform jitter and inaccurate position measurement will contaminate the MB SAR data, which may result in the multiplicative noise with phase errors, thereby leading to the decrease in the imaging quality. To address the above issues, this paper proposes a MB bistatic SAR (BiSAR) 3D imaging method based on the phase error calibration that combines the phase gradient autofocus (PGA) and energy balance intensity-squared optimization autofocus (EB-ISOA). Firstly, the signal model of the MB one-stationary (OS) BiSAR is established and the 3D imaging principle is presented, and then the phase error caused by platform jitter and inaccurate position measurement is analyzed. Moreover, combining the PGA and EB-ISOA methods, a 3D imaging method based on the phase error calibration is proposed. This method can improve the accuracy of phase error calibration, avoid the vertical displacement, and has the noise robustness, which can obtain the high-precision 3D BiSAR imaging results. The experimental results are shown to verify the effectiveness and practicality of the proposed MB BiSAR 3D imaging method.

1. Introduction

Tomographic synthetic aperture radar (TomoSAR) is an advanced synthetic aperture radar (SAR) technology that can effectively extract the three-dimensional (3D) spatial information of scene targets from multiple two-dimensional (2D) SAR images. It has evolved under the influence of computer tomography techniques and has seen rapid development in recent years, achieving good 3D imaging performance [1,2,3,4,5,6]. The elevation information of the target is retained within the phase of 2D SAR images. By acquiring coherent images from different angles, TomoSAR forms a synthetic aperture in the direction perpendicular to the sight of radar, enabling the resolution of targets in the tomographic direction [7,8,9,10,11,12]. TomoSAR has significant applications in various fields, such as urban building extraction, forest parameter retrieval, and glacier structure analysis.

In 1994, Knaell and Cardillo [13] first applied the tomography theory to the 3D radar imaging, providing a basis for analyzing the radar imaging resolution and performance. Subsequently, multiple configurations of the multi-baseline (MB) monostatic SAR (MonoSAR) were explored, and the experimental data was collected. In 2000, the German Aerospace Center (DLR) performed multiple near-parallel flights using a single aircraft and applied Fourier transform algorithms to focus in the tomographic direction, achieving a vertical resolution of 2.9 m [2]. In 2002, She [14] applied the TomoSAR theory to spaceborne systems for the first time, proposing a signal processing workflow and conducting experiments using the ERS-1 satellite data. In 2006, the Air Force Research Laboratory (AFRL) carried out the airborne MB circular SAR (CSAR) field flight tests [15]. In 2012, DLR used its airborne F-SAR system to conduct MB CSAR field tests in the Vordemwald region of Switzerland [16], acquiring a set of fully polarized MB CSAR data on the circular trajectory. In 2012, based on the compressed sensing method, Zhu [17,18] proposed a novel 3D imaging method and conducted relevant experiments using the TerraSAR-X data to verify the performance of the method.

Due to the characteristic of the MB bistatic SAR (BiSAR) with one transmission and multiple receptions, the time baselines between different SAR images can be shortened, and the data coherence can be better maintained. This has led researchers to explore MB BiSAR imaging configurations that effectively address the issue of the temporal decorrelation. In 2009, Duque [19] extended the two-receiver antenna configuration of the bistatic interferometric SAR (InSAR) to multiple receiving antennas, forming a 3D imaging configuration for the MB BiSAR. In 2019, Hu [20] conducted theoretical analyses of spaceborne MB BiSAR 3D imaging principles and carried out the first experimental verification using the data from Beidou inclined geosynchronous orbit navigation satellites. In the same year, Ge [21] applied a similar bistatic dataset to address the differential TomoSAR problem in urban areas by dividing the original problem into two subproblems, i.e., processing them separately, and merging the results to achieve the building observation.

Currently, the research on the MB BiSAR imaging is limited, mostly remaining in the stages of theoretical analysis and experimental validation. Many new configurations are being proposed to solve specific issues. In this paper, the MB one-stationary (OS) BiSAR imaging configuration is studied, aimed at reducing temporal decorrelation effects. This configuration uses a static radar and multiple airborne radars flying in the formation to collect the echo data, thus minimizing the time baselines between data acquisitions. For MB SAR observations, due to the platform jitter and inaccurate position measurement, multiple 2D images are easily contaminated by phase errors. This contamination can lead to the degraded focusing in the tomographic direction, including reduced and broadened main lobes, as well as elevated side lobes. In the MB OS BiSAR imaging configuration studied in this paper, such issues are similarly unavoidable. To address these problems, two classical phase error estimation methods have been proposed.

One method is the phase gradient autofocus (PGA) approach [22], which employs the principle of the inverse filtering to estimate phase errors, typically providing accurate approximations of phase errors caused by various factors. In 2019, the PGA method is applied to the holographic SAR (HoloSAR) for correcting inter-channel phase errors [23], based on the assumption that phase errors change slowly in space. In 2021, in order to solve the vertical displacement problem of the PGA method, Lu [24] estimated the height of significant elements based on the network construction (NC) method, and proposed the NC-PGA method which effectively alleviated the vertical displacement effect. In 2024, Huang [25] further addressed the issue of the vertical displacement by proposing a beamforming PGA (BF-PGA) method and a cumulative intensity projection registration method, Which further corrected the phase error and achieved better vehicle reconstruction results. The main limitation of this type of method was that it requires the assumption of spatially invariant phase errors, meaning it cannot estimate tomographic phase errors on a per-pixel basis, but can only estimate an overall approximate error for the block in the range-azimuth region.

Another category of autofocus methods is based on the sharpness metric criteria. This type of method takes the phase error to be estimated as a variable, uses the sharpness as the optimization objective function, and estimates the phase that needs to be compensated through search. The entropy and contrast of an image are often used as criteria for sharpness measurement [26,27]. In 2012, based on the minimization of the elevation profile entropy, Pardini [28] proposed a new calibration method and demonstrated its phase error correction capability using simulated and real data. However, there is a phenomenon of vertical displacement in the experimental results. To solve the problem of vertical displacement, Aghababaee [29] proposed a method based on the phase derivative constraint optimization, which can correct phase errors on a per-pixel basis without requiring assumptions about the scene. However, it involves solving a multidimensional optimization problem for each pixel. In [30], based on the intensity-squared optimization autofocus (ISOA), an autofocus method was proposed for 2D back-projected images, which allows a closed-form solution in single-variable cases through the coordinate descent method, avoiding the variable search process and significantly reducing processing time.

Based on the studied MB OS SAR imaging configuration, to solve the issues of the PGA method’s inability to estimate phase errors pixel-by-pixel and its susceptibility to the vertical displacement, this paper proposes a MB BiSAR 3D imaging method based on the phase error calibration combining the PGA and energy balance ISOA (EB-ISOA) methods. The main contributions are as follows.

• This paper explores a novel MB OS BiSAR 3D imaging configuration, introduces the concepts of incident angle and projection angle, presents the specific formulation of the de-ramping processing and validates its 3D imaging capability through both theoretical analysis and experimental evaluation.
• The criterion of maximizing the square of the intensity is introduced to estimate phase errors pixel-by-pixel. Moreover, this paper uses the preliminary results of the PGA method as the initial values for the EB-ISOA method, achieving better phase correction results.
• This paper improves the optimization model by proposing the energy balance processing to maintain the focus on weak targets. At the same time, this paper used the PS height obtained by the EB-ISOA method to compensate for the linear phase term introduced by the PGA method, which effectively alleviates the vertical displacement.

The remainder of this paper is formed as follows. In Section 2, the signal model of the MB OS BiSAR is established and the 3D imaging principle is presented, and then the characteristics of phase errors introduced by the platform jitter and inaccurate position measurements are analyzed. In Section 3, a 3D imaging method based on the phase error calibration combining the PGA and EB-ISOA is proposed. The effectiveness and practicability of the proposed MB BiSAR 3D imaging method are validated in Section 4. In Section 5, the projection relationship of the MB BiSAR imaging is discussed, and the calculation process of the projection angle and the incident angle is presented. In Section 6, the conclusions of the paper are presented.

2. Theoretical Foundation

In this section, the signal model of the MB OS BiSAR is established and the 3D imaging principle is presented. Then, the characteristics of phase errors introduced by the platform jitter and inaccurate position measurements are analyzed.

2.1. Signal Model and 3D Imaging Principle

The observation model of the MB OS BiSAR imaging is illustrated in Figure 1. In this setup, the receiver is mounted on an aircraft, while the transmitter remains stationary at a fixed location. The stationary transmitter and multiple airborne receivers form the MB BiSAR 3D imaging configuration. In Figure 1, $L_n\left(n=1,2,\dots ,N\right)$ represents the trajectory of the aircraft for the n-th flight, and N denotes the number of flight tracks. The trajectory $L_{N/2}$ is chosen as the reference trajectory. Based on the observation model depicted in Figure 1, the 2D imaging is performed separately for the echo data collected from each of the N trajectories, resulting in N 2D BiSAR images of the same target area.

Suppose there is an arbitrary point target P, with the coordinates $(x_P,y_P,z_P)$. The ranges from the receiver and the transmitter on the n-th trajectory to the target P at the slow time $\eta$ are denoted as $r_{R,n}\left(\eta \right)$ and $r_T\left(\eta \right)$, respectively. Assuming the signal emitted by the transmitter is a linear frequency modulation signal, the signal received on the n-th trajectory is:

$$\begin{array}{cc}\hfill s_n(\tau ,\eta )& ={\gamma }_P·w_{\mathrm{r}}\left[\tau -{\displaystyle \frac{r_{R,n}\left(\eta \right)+r_T\left(\eta \right)}{c}}\right]·w_{\mathrm{a}}(\eta -{\eta }_c)\hfill \\ \hfill & ·exp\left[-j2\pi f_c{\displaystyle \frac{r_{R,n}\left(\eta \right)+r_T\left(\eta \right)}{c}}+j\pi K_c{\left(\tau -{\displaystyle \frac{r_{R,n}\left(\eta \right)+r_T\left(\eta \right)}{c}}\right)}^2\right]\hfill \end{array}$$

(1)

where $\tau$ is the fast time, $\eta$ is the slow time, ${\gamma }_P$ is the complex scattering coefficient of the target P, c is the speed of the electromagnetic wave propagation, $r_{R,n}\left(\eta \right)$ is the slant range from the receiving station to the target P on the n-th trajectory at the slow time $\eta$, $r_T\left(\eta \right)$ is the slant range from the transmitting station to the target P at the time $\eta$, $w_r$ is the range envelope, $w_a$ is the azimuth envelope, ${\eta }_c$ is the synthetic aperture center time of the receiving station, $f_c$ is the center frequency, and $K_c$ is the chirp rate.

The back-projection (BP) method is used for the 2D BiSAR imaging, which better preserves the phase information compared to frequency-domain methods. The 2D complex SAR image on the imaging plane at the zero-height can be expressed as [31]:

$$\begin{array}{cc}\hfill g_n(x,y,0)& ={\gamma }_P·h_{\mathrm{Bi}}(x-x_{P^{\ast }},y-y_{P^{\ast }})\hfill \\ \hfill & ·exp\left[j2\pi \left(r_{\mathrm{Bi}}({\eta }_c;x,y,0)-r_{\mathrm{Bi}}({\eta }_c;x_P,y_P,z_P)\right)/\lambda \right]\hfill \end{array}$$

(2)

where $g_n(x,y,0)$ represents the complex value of the n-th SAR image at the grid point $(x,y,0)$, $(x_{P^{\ast }},y_{P^{\ast }},0)$ represents the coordinates of the target P projected onto the imaging plane according to the elliptical geometry, $\lambda =c/f_c$ is the signal wavelength, $h_{\mathrm{Bi}}(x,y)$ is the 2D response function of the point target in the BiSAR imaging, $r_{\mathrm{Bi}}({\eta }_c;x,y,0)$ is the round-trip slant range of the point $(x,y,0)$ at the slow time ${\eta }_c$, $r_{\mathrm{Bi}}({\eta }_c;x_P,y_P,z_P)$ is the round-trip slant range of $(x_P,y_P,z_P)$ at the slow time ${\eta }_c$.

For the sake of simplicity, it is assumed that $h_{\mathrm{Bi}}(x,y)$ is the 2D Dirac delta function [12]. After obtaining the imaging result using the BP algorithm, the phase caused by $r_{\mathrm{Bi}}({\eta }_c;x,y,0)$ can be compensated. In addition, because the transmitting station is fixed, its range to the target remains unchanged in N observations, causing the same phase in N SAR images, which can be incorporated into the complex scattering coefficient. Therefore, by omitting the coordinate index of the 2D image, $g_n(x,y,0)$ can be simplified as Equation (3) after the approximation and phase compensation.

$$g_n={\gamma }_P·exp\left(-{\displaystyle \frac{j2\pi r_{R,n}}{\lambda }}\right)$$

(3)

In Figure 2, the geometric relationship of the target projection in the MB OS BiSAR imaging is depicted. The positions of the receiving stations at the synthetic aperture center time are denoted as $R_n$$(n=1,2,\dots ,N)$. The position of the transmitting station is represented by T. In the BiSAR imaging, the target P is projected to the point $P^{\ast }$ on the imaging plane, so the round-trip slant ranges of P and $P^{\ast }$ are equal.

After obtaining N 2D images, the 2D image from the reference trajectory is selected as the master image, while the images obtained from other trajectories are considered as slave images for the image registration. Then, the registered images undergo de-ramping processing. For the n-th image, the de-ramping processing can be expressed as:

$$g_n={\gamma }_P·exp\left(-{\displaystyle \frac{j2\pi r_{R,n}^P}{\lambda }}\right)·exp\left({\displaystyle \frac{j2\pi r_{R,n}^{P^{\ast }}}{\lambda }}\right)$$

(4)

where $r_{R,n}^P$ represents the slant range from $R_n$ to P, and $r_{R,n}^{P^{\ast }}$ represents the slant range from $R_n$ to $P^{\ast }$. From the geometric relationship and far-field condition in Figure 2, $r_{R,n}^P$ and $r_{R,n}^{P^{\ast }}$ can be expressed as:

$$r_{R,n}^P\approx r+{\displaystyle \frac{{\left(s·sin\left({\alpha }_1+{\alpha }_2\right)-b_{n\perp }\right)}^2}{2r}}$$

(5)

$$r_{R,n}^{P^{\ast }}\approx r+{\displaystyle \frac{{\left(b_{n\perp }\right)}^2}{2r}}$$

(6)

where r represents the range from the grid point to the reference trajectory, ${\alpha }_1$ and ${\alpha }_2$ are the angles as shown in Figure 2, and ${\alpha }_1$ and ${\alpha }_2$ are analyzed in Section 5. s represents the range between P and $P^{\ast }$. $b_{n\perp }$ represents the vertical baseline of the n-th trajectory.

Therefore, Equation (4) can be expressed as:

$$g_n={\gamma }_P·exp\left[-{\displaystyle \frac{j\pi s^2·{sin}^2\left({\alpha }_1+{\alpha }_2\right)}{\lambda r}}\right]·exp\left[{\displaystyle \frac{j2\pi s·sin\left({\alpha }_1+{\alpha }_2\right)·b_{n\perp }}{\lambda r}}\right]$$

(7)

Because this paper is only concerned with the amplitude of the 3D image and not its phase, the first phase term in Equation (7) can be ignored [12]. Additionally, if multiple point targets are projected onto the same grid, Equation (7) can be rewritten as:

$$g_n=\int \gamma \left(s\right)·exp\left[j2\pi ·{\displaystyle \frac{sin({\alpha }_1+{\alpha }_2)·b_{n\perp }}{\lambda r}}·s\right]ds$$

(8)

where $\gamma \left(s\right)$ is the complex scattering coefficient of the target at the elevation s. It can be observed from Equation (8) that, the 2D image and the scattering coefficient form a Fourier transform relationship. The spatial frequency is ${\xi }_n=sin({\alpha }_1+{\alpha }_2){\displaystyle \frac{b_{n\perp }}{\lambda r}}$. Therefore, performing a Fourier transform along the s direction for each pixel position in N images yields the 3D imaging result. The specific form is shown as:

$$\widehat{\gamma }\left(s\right)=\int g_n·exp\left[-j2\pi ·{\displaystyle \frac{sin({\alpha }_1+{\alpha }_2)s}{\lambda r}}·b_{n\perp }\right]d{b}_{n\perp }$$

(9)

After obtaining the result of Equation (9), the imaging results can be obtained in the $x^{\prime }-y^{\prime }-s$ coordinate system established with the ground as the reference, which needs to be converted to the visually intuitive $x-y-z$ coordinate system.

$$\left\{\begin{array}{c}x=x^{\prime }\hfill \\ y=y^{\prime }+s·cos\left({\alpha }_2\right)\hfill \\ z=s·sin\left({\alpha }_2\right)\hfill \end{array}\right.$$

(10)

2.2. Phase Error Analysis

Due to the insufficient accuracy of the measuring equipment, inaccurate position measurement of the radar often occurs when the platform moves, leading to phase errors in the generated 2D image. Figure 3 shows the measurement geometric of the n-th trajectory. $R_n$ represents the actual position of radar on the n-th trajectory, and $R_n^{\prime }$ represents the estimated position of radar. There is a difference $(\Delta y,\Delta z)$ between $R_n$ and $R_n^{\prime }$. ${\theta }_n$ is the downward angle of radar for the n-th trajectory. $r_n$ is the true range, and $r_n^{\prime }$ is the range with deviation. The range error and phase error can be represented as:

$$\Delta r_n=r_n^{\prime }-r_n=\Delta y\sin {\theta }_n-\Delta z\cos {\theta }_n$$

(11)

$${\phi }_n={\displaystyle \frac{2\pi \Delta r_n}{\lambda }}$$

(12)

The phase errors on each image vary with the change of ${\theta }_n$. In addition, due to the different radar position deviations on each trajectory, the phase error of multiple trajectories will vary with spatial position changes. The phase error between different trajectories is independent of each other and needs to be estimated separately. The existence of phase error terms may lead to problems such as the blurring, distortion, and reduced resolution in the tomographic profile. Substituting Equation (12) into Equation (8) yields the result of the phase error contamination as:

$$g_n=exp\left(j{\phi }_n\right)·\int \gamma \left(s\right)·exp\left[j2\pi ·{\displaystyle \frac{sin({\alpha }_1+{\alpha }_2)·b_{n\perp }}{\lambda r}}·s\right]ds$$

(13)

3. Proposed MB BiSAR 3D Imaging Method

In this section, a MB BiSAR 3D imaging method based on the phase error calibration combining the PGA and EB-ISOA is proposed. The overall process of the proposed 3D imaging method is shown in Figure 4, which mainly includes five distinct steps:

• Registering images and de-ramping;
• Selecting PS based on the amplitude dispersion index (ADI) criteria;
• Using the PGA method to estimate the phase error and obtain coarse calibration results;
• Using the EB-ISOA method to process coarse calibration results and obtain precise calibration results;
• Performing the tomographic imaging and coordinate transformation to obtain final 3D imaging results.

Since the first step has already been presented in Section 2, the remaining part of this section will cover the last four steps.

3.1. PS Selection Based on ADI Criteria

The ADI criterion is used to screen PS among candidate pixels [32], and it can be expressed as:

$$\mu =\sqrt{N\sum _{n=1}^N{\left|g_n\right|}^2/{\left(\sum _{n=1}^N\left|g_n\right|\right)}^2-1}$$

(14)

Select pixels with $\mu <Th{r}_{\mu }$ as the PS, and $Th{r}_{\mu }$ is usually set to 0.25 [33]. Due to the fact that the ADI considers the amplitude stability of the selected pixels on different baselines, rather than directly considering the amplitude level. Therefore, when weakly scattered pixels meet the ADI, they can also be selected as PS. After calculating ADI and filtering the pixels, the data of the q-th PS selected from the imaging area on the n-th image can be represented as:

$$g_{q,n}={\gamma }_{q}exp(j{\phi }_n+j2\pi {\xi }_n{s}_q)+w_{q,n}$$

(15)

where ${\gamma }_q$ is the complex scattering coefficient of the q-th PS, $s_q$ is the elevation of the q-th PS, and $w_{q,n}$ is the additive noise [34].

3.2. Coarse Phase Error Estimation Based on PGA

The PGA method is based on the assumption that phase error can be considered invariant within a certain range [35,36]. The PGA method mainly includes the following steps:

Firstly, the linear phase term removal. In this paper, the discrete Fourier transform is used to estimate $s_q$. After compensating for this linear phase term, Equation (15) can be expressed as:

$$g_{q,n}^{\prime }={\gamma }_{q}exp\left(j{\phi }_n\right)+w_{q,n}^{\prime }$$

(16)

where $w_{q,n}^{\prime }=w_{q,n}·exp(-j2\pi {\xi }_n{s}_q)$.

Secondly, the phase gradient estimation. Based on the principle of maximum likelihood estimation, the phase gradient $\Delta \widehat{\phi }$ can be estimated as:

$$\Delta {\widehat{\phi }}_n=arg\left\{\sum _{q=1}^Q[{\left(g_{q,n-1}^{\prime }\right)}^{\ast }·g_{q,n}^{\prime }]\right\}$$

(17)

where $arg(·)$ denotes the phase extraction operation, ${(·)}^{\ast }$ denotes the conjugate operation, and Q is the number of the PS.

Finally, the phase error estimation. By summing up $\Delta {\widehat{\phi }}_n$, ${\phi }_n$ can be obtained as:

$${\widehat{\phi }}_n=\left\{\begin{array}{cc}0,\hfill & n=1\hfill \\ \sum _{i=2}^n\Delta {\widehat{\phi }}_i,\hfill & n=2,3,\dots ,N.\hfill \end{array}\right.$$

(18)

Considering the influence of noise, the estimation result of phase error may not be close to the true phase error. Therefore, in the PGA method, the strategy of cyclic iteration is used to gradually improve the accuracy of phase error estimation. The convergence condition of the iterative process can be expressed as [22]:

$$\sum _{n=1}^N{\left|{\widehat{\phi }}_n^i-{\widehat{\phi }}_n^{i-1}\right|}^2<Th{r}_{iter}$$

(19)

where ${\widehat{\phi }}_n^i$ is the results of phase error estimation at the i-th iteration and the $Th{r}_{iter}$ is usually set to 0.001 [13]. In this stage of the proposed method, it only needs the coarse estimation of the phase error, so the convergence condition can be relaxed. In this paper, it is set to 0.01.

3.3. Precise Phase Error Estimation Based on EB-ISOA

The ISOA method is based on the intensity-squared criterion and does not require assumptions of the phase error invariance or noise specific statistical models. It can correct phase errors pixel-by-pixel. The phase error estimation process of the ISOA method can be represented by:

$$\widehat{\mathit{\phi }}=\underset{\mathit{\phi }}{max}\left\{C[f(s_1,\mathit{\phi }),f(s_2,\mathit{\phi }),\cdots ,f(s_N,\mathit{\phi })]\right\}$$

(20)

where $C[·]$ is the objective function, $\mathit{\phi }=[{\phi }_1,\dots ,{\phi }_n]$ is the calibration phase vector to be calculated, $s_d$ is the elevation, and $f(s_d,\mathit{\phi })={\left|\widehat{\gamma }\left(s_d,\mathit{\phi }\right)\right|}^2,(1\le d\le N)$ represents the intensity at $s_d$ in the tomographic focusing result.

The intensity-square is chosen as the measure of the sharpness, and maximizing the intensity-square can maximize the variance of the intensity relative to its average value, resulting in the optimal sharpness of the focused signal [27]. The objective function is defined as:

$$C\left[f(s,\mathit{\phi })\right]=\sum _{d=1}^N{f}^2(s_d,\mathit{\phi })$$

(21)

According to Equation (9), it can be concluded that:

$$\begin{array}{cc}\hfill \widehat{\gamma }(s_d,\widehat{\mathit{\phi }})& =\sum _{n=1}^N{g}_n·{\mathrm{e}}^{-j2\pi {\xi }_n{s}_d}·{\mathrm{e}}^{-j{\widehat{\phi }}_n}\hfill \\ & =\sum _{n=1}^N{\nu }_{d,n}·{\mathrm{e}}^{-j{\widehat{\phi }}_n}\hfill \end{array}$$

(22)

where ${\nu }_{d,n}=g_n·{\mathrm{e}}^{-j2\pi {\xi }_n{s}_d}$.

The optimization problem of Equation (20) does not have a closed-form solution and is usually solved using numerical iterative methods. The coordinate descent method can be to transform the above problem into the following form:

$${\widehat{\phi }}_n^{i+1}=\underset{\phi }{max}\left\{C\left[f(\mathit{s},{\widehat{\phi }}_1^{i+1},\cdots ,{\widehat{\phi }}_{n-1}^{i+1},\phi ,{\widehat{\phi }}_{n+1}^i,\cdots ,{\widehat{\phi }}_N^i)\right]\right\}$$

(23)

where $\phi$ is the phase error to be estimated in the $(i+1)$-th iteration, Equation (22) can be converted into:

$$\begin{array}{cc}\hfill \widehat{\gamma }(s_d,\widehat{\phi })& =\sum _{l=1}^{n-1}{\nu }_{d,l}{\mathrm{e}}^{-j{\widehat{\phi }}_l^{i+1}}+\sum _{l=n+1}^N{\nu }_{d,l}{\mathrm{e}}^{-j{\widehat{\phi }}_l^i}+{\nu }_{d,n}{\mathrm{e}}^{-j\widehat{\phi }}\hfill \\ & =p_d+{\mathrm{e}}^{-j\widehat{\phi }}{q}_d\hfill \end{array}$$

(24)

where $p_d={\sum }_{l=1}^{n-1}{\nu }_{d,l}{\mathrm{e}}^{-j{\widehat{\phi }}_l^{i+1}}+{\sum }_{l=n+1}^N{\nu }_{d,l}{\mathrm{e}}^{-j{\widehat{\phi }}_l^i}$, and $q_d={\nu }_{d,n}$. Equation (24) conforms to Equation (9) in [30], and subsequent implementation details can be found in [30].

However, the intensity squared criterion tends to focus on strong energy sources that are affected by phase noise interference, where the magnitude of the intensity squared is mainly determined by the strong scattering points in the scene. Therefore, the estimation of phase error tends to prioritize focusing on strong scattering points to obtain larger intensity squared values. This results in the focusing result after phase correction being mainly determined by the strong energy component. To address this issue, we introduce image energy balance to improve the ISOA method, which also named EB-ISOA. We propose a balancing operator $\Omega (·)$ to balance the energy distribution of ${\nu }_{d,n}$. To avoid loss of phase information, the balancing operator can be expressed as:

$$\Omega \left({\nu }_{d,n}\right)=\Psi \left(V_{d,n}\right){\mathrm{e}}^{j{\phi }_{d,n}}$$

(25)

where $V_{d,n}=\left|{\nu }_{d,n}\right|$, ${\phi }_{d,n}=angle\left({\nu }_{d,n}\right)$, and the function of $\Psi (·)$ is to reduce the dynamic range of $V_{d,n}$, that is, to reduce the difference in the contribution of projection energy from different regions to $C\left[f\right(s,\mathit{\phi }\left)\right]$. Common nonlinear operators can be used to reduce the dynamic range of data energy. In our data processing, the logarithmic function $log(·)$ is chosen as the balancing operator.

The ISOA method does not introduce linear phase error terms during the optimization, and therefore does not produce the vertical displacement. However, the PGA-based methods may not be accurate in estimating the elevation of the PS, resulting in the vertical displacement when eliminating the linear phase term, and this vertical displacement may gradually increase during the iteration process. So, before using the EB-ISOA method to process the coarse calibration results obtained by the PGA method, the data of the PS can be calibrated using the EB-ISOA method to obtain the elevation of the PS. Next, the PS elevation obtained by the EB-ISOA and PGA methods is used to perform the elevation calibration and obtain the vertical displacement amount of the PGA calibration result. The specific elevation calibration process is as follows:

$$\Delta s={\displaystyle \frac{1}{Q}}\sum _{q=1}^Q\left(s_q^{A_1}-s_q^{A_2}\right)$$

(26)

where $s_q^{A_1}$ expresses the elevation of the q-th PS obtained by the EB-ISOA method, and $s_q^{A_2}$ expresses the elevation of the q-th PS obtained by the PGA method.

After obtaining the vertical displacement amount, the corresponding linear phase term is compensated to the PGA calibration result, which can complete the elevation calibration of the PGA result. After completing the elevation calibration, the EB-ISOA method can be used to process the coarse calibration data of the PGA method pixel-by-pixel.

3.4. Tomographic Imaging and Coordinate Transformation

After the phase error calibration, it is necessary to perform focusing processing on the tomographic signal. From Equation (9), it can be seen that using discrete Fourier transform can achieve focusing on tomographic signals. However, considering that noise often exists in actual signals, which can affect the focusing performance of discrete Fourier transform, this paper chooses the ABF method with anti noise capability to achieve focusing of tomographic signals. The ABF method is based on filtering signals with different spatial frequencies through a spatial filter $\mathbf{h}\left(s\right)$ of order N, so that the reflected signal of the target at the elevation s is not distorted and the amplitude of signals at other frequencies is reduced. By using the ABF method, the estimation of intensity distribution can be denoted as:

$${\widehat{P}}_{ABF}\left(s\right)={\mathbf{h}}^H\left(s\right)\widehat{\mathbf{R}}\mathbf{h}\left(s\right)$$

(27)

where $\mathbf{h}\left(s\right)$ can be expressed as:

$$\mathbf{h}\left(s\right)={\displaystyle \frac{{\widehat{\mathbf{R}}}^{-1}\mathbf{a}\left(s\right)}{{\mathbf{a}}^H\left(s\right){\widehat{\mathbf{R}}}^{-1}\mathbf{a}\left(s\right)}}$$

(28)

where $\mathbf{a}\left(s\right)={\left[exp\left(j2\pi {\xi }_{1}s\right),exp\left(j2\pi {\xi }_{2}s\right),\cdots ,exp\left(j2\pi {\xi }_{N}s\right)\right]}^{\mathrm{T}}$, ${(·)}^H$ denotes the Hermitian operator, and $\widehat{\mathbf{R}}$ is the covariance matrix which can be expressed as:

$$\widehat{\mathbf{R}}={\displaystyle \frac{1}{M}}\sum _{m=1}^M{\mathbf{g}}_m{{\mathbf{g}}_m}^H$$

(29)

The ${\mathbf{g}}_{m=1,\dots ,M}$ can be chosen from the pixels around a selected pixel. After the tomographic imaging, Equation (10) can be used to obtain the visually intuitive results.

4. Experimental Results and Analysis

In this section, the 3D imaging performance of the proposed MB OS BiSAR 3D imaging method is verified through simulation experiments, which demonstrates the effectiveness and practicality of the proposed method. Based on the imaging geometry configuration shown in Figure 2, the system parameters and baseline configuration in

Table 1. System parameters.

Parameter	Value
Center frequency	9.375 GHz
Bandwidth	100 MHz
Sampling rate	120 MHz
Pulse repetition frequency	300 Hz
Transmitting station height	500 m
Receiving station height (Reference trajectory)	10 km
Receiving station speed	100 m/s
Number of Receiving station trajectories	25
Trajectory interval	20 m

are used to generate the simulation data for three target scenarios. Notably, 15,000 m is the ground range between the scene center and the reference trajectory. Two experiments are conducted using the parameters in Table 1. The first experiment is a point targets imaging experiment in an environment without tomographic phase error, aimed at verifying the 3D imaging capability of the MB OS BiSAR and the effectiveness of the proposed 3D imaging method (eliminating the phase error calibration step). The second experiment is the tomographic phase error calibration experiment, aimed at verifying the effectiveness and practicality of the proposed phase error calibration method.

4.1. Three-Dimensional Imaging Capability Experiment

As shown in Figure 5, five point targets are set with the following coordinates: A $(-10,-10,4)$m, B $(-10,10,-4)$m, C $(10,-10,-4)$m, D $(10,10,4)$m, and E $(0,0,0)$m. In this experiment, the trajectory is uniformly distributed along the height direction, and the total length of baseline is 480 m. $30\mathrm{m}\times 50\mathrm{m}$ (X-axis × Y-axis) is the size of the imaging scene. The observation angle of the reference trajectory towards the scene center is 56.34°. The theoretical elevation resolution is 1.50 m, and the unambiguous elevation is 36.04 m. The scattering coefficients of the point targets are all 1 ${\mathrm{m}}^2$, with the highest target at 4 m and the lowest target at −4 m. The OS BiSAR data for 25 voyages has been obtained through the simulation.

The BP algorithm is used to perform the 2D imaging processing on the echo data received by the receiving radar on all trajectories, and 25 BiSAR images are obtained. The imaging results obtained by performing the de-ramping processing and Fourier transform focusing on each image are shown in Figure 6. The discrete points of the sidelobes in the 3D imaging results of point targets can affect the 3D observation of the imaging results. To better display the imaging results without changing the imaging results, only discrete points with amplitudes greater than 0.05 times of the maximum amplitude value are shown in Figure 6. From the imaging results of the point targets, it can be seen that its 3D shape approximates the discretization of the 3D sinc function, which can visually distinguish the point targets more clearly, verifying the 3D imaging capability of the proposed MB OS BiSAR method. In order to have a clearer view of the coordinate details of the 3D imaging results, Figure 7 shows the results of the 3D imaging results from three perspectives after coordinate transformation. From the results of the three perspectives, it can be seen that the imaging results of the point target exhibit the shape of a 2D sinc function on each projection plane, which is in line with theoretical expectations.

Due to the fact that the BP algorithm is used for the 2D imaging and its focusing ability has been validated by numerous researches, the following analysis will only focus on the focusing effect in the tomographic direction. The tomographic profiles of the imaging results of five point targets before coordinate transformation are shown in Figure 8. From Figure 8, it can be seen that the tomographic profiles of the five targets all follow the form of a one-dimensional sinc function, indicating that the targets have been well-focused in the tomographic direction. Since Figure 8 shows the tomographic profiles of the imaging results without coordinate transformation, the peak positions in each subgraph are not the actual heights of the point targets.

In order to further quantitatively demonstrate the performance of the 3D imaging of the MB OS BiSAR,

Table 2. Imaging quality and focusing position of point targets imaging result.

		Expected Values	Point Targets
			A	B	C	D	E
Range	IRW (m)	1.39	1.40	1.41	1.45	1.45	1.42
	PSLR (dB)	−13.2	−13.32	−13.20	−13.29	−13.48	−13.09
	ISLR (dB)	−10	−10.54	−10.61	−10.86	−10.90	−10.78
Azimuth	IRW (m)	0.65	0.66	0.66	0.66	0.66	0.66
	PSLR (dB)	−13.2	−13.27	−13.26	−13.27	−13.26	−13.27
	ISLR (dB)	−10	−10.13	−10.12	−10.13	−10.13	−10.12
Elevation	IRW (m)	1.50	1.23	1.22	1.23	1.23	1.22
	PSLR (dB)	−13.2	−12.85	−12.89	−12.93	−12.90	−12.93
	ISLR (dB)	−10	−9.79	−9.82	−9.71	−9.77	−9.76
True position (m)			(−10, −10, 4)	(−10, 10, −4)	(10, −10, −4)	(10, 10, 4)	(0, 0, 0)
Focusing position (m)			(−10.00, −9.86, 4.33)	(−10.00, 10.06, −3.90)	(10.06, −10.00, −3.90)	(10.06, 10.20, 4.33)	(0.06, 0.14, 0.26)

presents the focus positions of the five point targets and the three performance parameters of the point targets. From the impulse response width (IRW) indicator in Table 2, it can be seen that point targets have a certain resolution in the three direction, while the peak side lobe ratio (PSLR) and integrated side lobe ratio (ISLR) indicators indicate a significant attenuation of the sidelobes of point targets relative to the main lobe in the three direction. Observing the expected and measured values of the indicators at different point targets in the range and azimuth directions, it can be found that the measured and expected values of the IRW of the five point targets are close to each other, which is in line with theoretical expectations. However, PSLR and ISLR are slightly lower than the expected values, which meets the requirements because the smaller the PSLR and ISLR, the more concentrated the point target energy is in the main lobe, and the less energy leakage. Observing the three indicators in the elevation direction, the IRW measurement values of the five points are slightly lower than the expected values, while the PSLR and ISLR are slightly higher than the expected values. This is mainly because the sidelobe section of the tomographic direction is not perpendicular to the tomographic axis, resulting in the tomographic data sampled along the tomographic axis not being completely in the ideal sinc distribution form, causing slight deviations in the calculation results of the indicators. However, it can also be seen that the level of energy focusing in the tomographic direction is close to the performance of the sinc function. In addition, it can be seen that the 3D coordinates of the point targets in the 3D imaging results are close to their true coordinates, but not equal to their true coordinates. This is due to the approximation operation during the de-ramping process, registration errors between different images and the insufficient up-sampling factor during the position calculation. The above experimental results objectively verify that the proposed MB BiSAR configuration has the ability to obtain resolution and energy focusing ability in three directions: range, azimuth and elevation.

4.2. Phase Error Calibration Experiment

4.2.1. Experiment for Simple Targets and Scenarios

Set the target $Ta{r}_z$ composed of multiple point targets, as shown in Figure 9a. Notably, $20\mathrm{m}\times 40\mathrm{m}$ (X-axis × Y-axis) is the size of the imaging scene. The scattering coefficients of the point targets that make up the target $Ta{r}_z$ are all 1 ${\mathrm{m}}^2$, and the highest height of the point targets is 5 m and the lowest height is −5 m. The distance intervals between adjacent points in the X, Y and Z directions are 2 m, 2 m and 2.5 m, respectively. The 3D imaging results of the target $Ta{r}_z$ without the phase error obtained according to the above settings are shown in Figure 9b.

A set of Gaussian distributed phase errors along the tomographic direction within the range of $[-\pi ,\pi ]$ radians is set, with a mean of 0 and a standard deviation of 0.32$\pi$ [29]. Based on the above data, the phase errors that vary linearly along the range and azimuth directions are added to each registered SAR image [37]. Considering the spatial relationship between the scene and the radar and the scale of the imaging scene, the phase error slowly changes on each image. In addition, the signal-to-noise ratio (SNR) is set to 5 dB.

The 2D imaging result generated from the acquisition signal corresponding to the reference trajectory is shown in Figure 10a, while the PS extracted by the ADI criterion are shown in Figure 10b. After statistics, the number of extracted PS is 468. Comparing the 2D imaging results with the PS extraction results, the position of PS basically corresponds to the position of isolated scattered points on the 2D image, and the energy aliasing part of the dense points located in the middle of the Y-axis of Figure 10a has not been extracted. This indicates that after screening by ADI, the extracted PS are the scattering points with isolated properties.

The 3D imaging results after the phase error contamination and phase error calibration using different methods are shown in Figure 11. From Figure 11a, it can be seen that the due to the influence of phase error, the imaging results have obvious blurring, making it difficult to obtain the shape and structural information of the target. After the phase error calibration, the focusing effect of the target is improved and the profile is clear. From Figure 11b–d, it can be seen that using the PGA, ISOA, and BF-PGA methods individually can each produce imaging results with some improvement, but they all exhibit certain problems in the reconstructed vertical plane. The results obtained by the PGA method showed an overall displacement, while the suppression of sidelobes is insufficient. The results obtained by the BF-PGA method can show the more distinct profile, but there is also the vertical displacement. The PGA method and BF-PGA method both rely on the assumption of phase error invariance and can only compensate for the same phase for all pixels in the same image, so the phase compensation effect for different pixels will vary. For some pixels, the tomographic focusing effect may deteriorate after compensation by PGA or BF-PGA methods. In addition, both PGA and BF-PGA methods require estimation of the height of the PS, and whether using DFT or BF methods, the estimated height of the PS may deviate, leading to vertical displacement in the height of the 3D image. The results obtained by the ISOA method do not produce overall offset, but the sidelobe suppression levels of different pixels are inconsistent, and the profiles are not smooth enough. This is mainly because the original ISOA method overly focuses on strong scattering targets, to some extent plundering the energy of weak scattering targets, leading to the phenomenon of over focusing. From Figure 11e, the proposed phase error calibration method can achieve the higher calibration effects and better sidelobe suppression effect. At the same time, there is no vertical displacement in the results, and a clear target profile can be seen. Comparing Figure 11c,e, it can be observed that after energy balance processing, targets with weaker energy can achieve better focusing without obvious over-focusing phenomena.

Figure 12 shows the tomographic profile of a selected PS in the results obtained by the above methods. It can be seen that the results of the PGA and BF-PGA methods both have the vertical displacement, while the elevation of the ISOA method and the proposed method remained stable without obvious vertical displacement. In addition, the proposed method has the good sidelobe suppression effect, which is superior to other comparison methods. For images with phase errors, even after the phase calibration processing, their phase is still close to the error free phase, and it is difficult to fully recover to the clean phase. In this case, the sidelobes in the case of ’without phase errors’ are weaker than those in the proposed method, indicating that there is less energy leakage and better focusing effect in the absence of the error. For the tomographic profile of the BF-PGA method, it can be found that its sidelobe level is not stably attenuated. This phenomenon is mainly due to the assumption that the BF-PGA method can only compensate for one phase value for each image based on the phase error space invariance, resulting in inconsistent phase error correction effects for different PS points. By comparing the results of the proposed method, it can be found that the sidelobes gradually attenuate, and the attenuation trend is consistent with the trend without the phase error, indicating that the signal energy is mainly concentrated in the main lobe with the less energy leakage.

The entropy and contrast values of different phase methods on the 3D BiSAR images mentioned above are given in

Table 3. Comparison of the image quality evaluation indicators for imaging results of the target $Ta{r}_z$.

Methods	Without Calibration	PGA Method	ISOA Method	BF-PGA Method	Proposed Method
IE	10.62	9.27	9.15	9.29	8.63
IC	4.09	8.30	8.74	7.71	9.21

, which are used to comprehensively evaluate the imaging quality. The image entropy can be obtained by [38]:

$$\mathrm{IE}=-\sum _{a=1}^{N_a}\sum _{b=1}^{N_b}\sum _{c=1}^{N_c}p(a,b,c)ln\left[p(a,b,c)\right]$$

(30)

where $N_a$, $N_b$, and $N_c$ are the size of the 3D BiSAR image, $p(a,b,c)$ is the intensity probability density of the selected pixel. $p(a,b,c)$ can be expressed as:

$$p(a,b,c)={\displaystyle \frac{{\left|I(a,b,c)\right|}^2}{{\sum }_{a=1}^{N_a}{\sum }_{b=1}^{N_b}{\sum }_{c=1}^{N_c}{\left|I(a,b,c)\right|}^2}}$$

(31)

where $I(a,b,c)$ is the intensity of pixel. The image contrast is given by [25]:

$$\mathrm{IC}={\displaystyle \frac{\sqrt{M{(I^2(a,b,c)-M\left(I^2(a,b,c)\right))}^2}}{M\left(I^2(a,b,c)\right)}}$$

(32)

where $M(·)$ denotes the expectancy factor.

From Table 3, it can be seen that the proposed method has the highest contrast and lowest entropy on 3D BiSAR images. Overall, compared with other methods, the proposed method has demonstrated its higher focusing degree, smaller sidelobes, more detailed inversion results and superior calibration performance.

This paper sets up a Monte Carlo experiment with the SNR as the variable. When examining the effectiveness of different algorithms under different SNR, Gaussian phase noise and uniform distributed phase noise are used to contaminate the data. The mean of Gaussian phase noise is 0, and the standard deviation is 0.32$\pi$. Uniform distributed phase noise follows a uniform distribution of $[-\pi ,\pi ]$. The number of repeated experiments for each variable’s discrete value is set to 100. Figure 13 and Figure 14 show the variation curves and standard deviations of entropy and contrast without phase error calibration and after phase error calibration, respectively, using four methods in Gaussian phase noise and uniform phase noise environments.

From Figure 13 and Figure 14, it can be seen that as the SNR increases, the entropy and contrast of the processing results of these methods decrease and increase, respectively. Compared to the other three methods, the proposed method has the smallest entropy and the highest contrast in both noisy environments. This indicates that the proposed method further corrects pixel-by-pixel based on the initial values obtained by PGA, achieving finer phase error correction of the data and obtaining higher quality 3D focused images. In addition, the PGA-based methods are prone to failure in the low SNR, while the ISOA method and the proposed method still have certain correction effects. This is mainly because at low SNR, the assumption that the phase error remains unchanged is difficult to hold, making it difficult for PGA based methods to be effective. This further demonstrates the advantage of the proposed method having a wider range of applicable scenarios. In terms of standard deviation in Figure 13 and Figure 14, the proposed method has slightly higher overall standard deviation for IE and IC than PGA and BF-PGA methods, and the standard deviation variation range of the three methods is not large and relatively stable.

4.2.2. Experiment for Complex Targets and Natural Scenes

To further demonstrate the effectiveness of the proposed algorithm, an experiment with complex targets and natural scenes is conducted. Notably, $40\mathrm{m}\times 40\mathrm{m}$ (X-axis × Y-axis) is the size of the imaging scene. Figure 15 shows the 3D perspective view of the complex target model and the projection plane view graphics. In addition, in order to be closer to the actual situation of phase error generation, unlike the phase error generation process in Section 4.2.1, this experiment directly adds perturbed position errors to the set aircraft trajectory, thereby causing distance errors and further simulating the phase error generation process in natural environments. The relevant system parameters of the experiment are consistent with the parameters in Table 1.

Due to disturbance errors, the generated 2D SAR images may appear defocused. In this paper, the autofocus BP method is used for 2D autofocus processing. The 2D SAR image under trajectory disturbance is shown in Figure 16a, which exhibits obvious defocusing and blurring phenomena, making it difficult to see the approximate outline and structure of the target. The 2D SAR image after 2D autofocus processing is shown in Figure 16b. The clear target structure and contour can be seen in the image, and the defocusing phenomenon of the 2D image is well suppressed.

Perform image registration on the 2D autofocus processed image stack using the reference trajectory image as the reference image, and perform oblique processing to obtain preliminary processed image stack data. The proposed phase calibration method and the other three phase calibration methods are used to perform phase error calibration and unified ABF focusing processing on the image stack data, and the results obtained are shown in Figure 17 and Figure 18.

Figure 17a shows the result of direct focusing without phase error correction, indicating that under the influence of phase error, the tomographic signal was not successfully compressed and the focusing basically failed. Figure 17b,d demonstrate that PGA and BF-PGA can effectively correct phase, but the correction effect varies significantly among different pixel points. The pixel points located in the middle area of the X-axis have better correction effect and more concentrated energy, while the pixel points located on both sides of the X-axis have more dispersed energy. In addition, the processing results of PGA and BF-PGA have a slight degree of vertical displacement. The processing result of Figure 17c is poor, as the weak point energy is plundered by the strong point, making it difficult to focus on the weak point and the overall outline unclear. Observing Figure 17e, the focused image corrected by the proposed phase correction method has a clear contours, good focusing of both strong and weak scattered points and no obvious vertical displacement. Observe the three areas circled by red lines in the subgraph of Figure 17, which are partially weak scattering areas in the scene. Comparing the results of four methods for processing regions A, B, and C, region A in Figure 17b has a few weak targets lost, while regions A and B in Figure 17c have lost most weak targets. However, regions A and B in Figure 17d and Figure 17e have achieved good focusing on weak targets, but the focusing depth in Figure 17e is deeper. For region C, four methods can achieve the focus on the target, but the results of the proposed method have better differentiation of the target. Figure 18 shows the observation images of the 3D imaging results from three perspectives. Comparing the three views of the target in Figure 15, it can be seen that these four methods have basically achieved the inversion of the contour of the target, but the contour of the result processed by the ISOA method is relatively blurry, losing the construction of detailed contours such as weak targets. In addition, the imaging results processed by PGA and BF-PGA methods both exhibit obvious vertical displacement. These results reflect the good focusing ability of the proposed method on weak targets, which can better invert the structure of targets in strong and weak target environments. The IE and IC of the image after phase error correction are given in

Table 4. Comparison of the image quality evaluation indicators.

Methods	Without Calibration	PGA Method	ISOA Method	BF-PGA Method	Proposed Method
IE	11.04	9.72	9.73	9.67	9.26
IC	5.45	10.94	10.32	11.03	11.76

. It can also be seen from the table that the focused image of the proposed method has the lowest image entropy and the highest contrast, resulting in the best overall phase calibration effect on the target.

5. Discussion

5.1. Analysis of BiSAR Projection

In the MonoSAR imaging, the projection relationship from the target to the imaging plane conforms to a circular projection relationship. While in the OS BiSAR imaging, the projection relationship from the target to the imaging plane conforms to an elliptical projection relationship. Under far-field conditions, due to the observation scene being much smaller than the range between the radar and the target, both circular and elliptical projection relationships can be approximated as the linear projection relationships. However, although they both can be approximated as a linear projection relationship, the angle between the approximated line (i.e., the tomographic direction) and the sight direction of the radar is not the same. In the MonoSAR imaging, the sight direction of the radar and the tomographic direction can be considered perpendicular. While in the BiSAR imaging, the sight direction of radar and the tomographic direction generally do not satisfy a vertical relationship, and projection angles need to be introduced for the geometric analysis.

The projection relationship under MonoSAR and BiSAR configurations is shown in Figure 19. As mentioned above, the projection relationship of both configurations can be approximated as a linear projection. The direction of this straight line is actually the tangent direction at the intersection of the imaging plane and the circle or ellipse. In Figure 19a, it can be seen that in the MonoSAR configuration, the tangent of the circle at $P^{\ast }$ is perpendicular to the $R_{OS}{P}^{\ast }$, i.e., perpendicular to the sight direction of the radar towards $P^{\ast }$.

Assuming that $P^{\ast }$ is at the center of the scene, $R_{OS}{P}^{\ast }$ can be regarded as the slant range direction, while the tangent direction is the tomographic direction. Under the MonoSAR configuration, the tomographic direction is perpendicular to the slant range direction. In Figure 19b, it can be seen that in the BiSAR configuration, the tangent direction of the ellipse at $P^{\ast }$ may not be perpendicular to the sight $R{P}^{\ast }$ of the radar. Therefore, it is necessary to use elliptical geometric relationships to calculate the angle between the tangent direction s and the line of sight direction $R{P}^{\ast }$, and this angle is set as $\theta$.

5.2. Calculation of Projection Angle and Incident Angle

According to the geometric relationship in Figure 19b, it can be inferred that $\theta =\pi -{\alpha }_1-{\alpha }_2$. ${\alpha }_1$ is the incident angle from the radar to the center of the scene, and ${\alpha }_2$ can be regarded as the projection angle from the target P to the center of the scene $P^{\ast }$. Assuming at the center moment of the synthetic aperture, the coordinates of R are $(x_R,y_R,z_R)$, and the coordinates of $P^{\ast }$ are $(x_{P^{\ast }},y_{P^{\ast }},z_{P^{\ast }})$. Then, we can obtain $x_R=x_{P^{\ast }}$ and ${\alpha }_1$ can be given by:

$${\alpha }_1=arctan({\displaystyle \frac{z_R-z_{P^{\ast }}}{y_R-y_{P^{\ast }}}})$$

(33)

Compared to the angle ${\alpha }_1$, the calculation of the projection angle ${\alpha }_2$ is more complex. As shown in Figure 19b, based on the original O-XYZ coordinate system, a Cartesian coordinate system O-UV (with the U axis parallel to the Y axis and the V axis parallel to the Z axis) is established with the center of $RT$ as the origin. Assuming the coordinates of T are $(x_T,y_T,z_T)$, and considering that R, T, and $P^{\ast }$ are all in the O-UV plane, then $x_T=x_{P^{\ast }}=x_R$. Furthermore, the coordinates $x_O,y_O,z_O$ of the midpoint of $RT$ in the O-XYZ coordinate system can be calculated as:

$$\left\{\begin{array}{c}{x}_O=x_R\hfill \\ y_O={\displaystyle \frac{y_T+y_R}{2}}\hfill \\ z_O={\displaystyle \frac{z_T+z_R}{2}}\hfill \end{array}\right.$$

(34)

In the O-UV coordinate system, the coordinates of $P^{\ast }$ are $(u_{P^{\ast }},v_{P^{\ast }})$, and $u_{P^{\ast }}=y_{P^{\ast }}-y_O$, $v_{P^{\ast }}=z_{P^{\ast }}-z_O$. According to the properties of ellipses, the length of the major axis of the ellipse $2a$ can be calculated as

$$2a=r_R^{P^{\ast }}+r_T^{P^{\ast }}$$

(35)

where $r_R^{P^{\ast }}$ represents the range from $P^{\ast }$ to R, and $r_T^{P^{\ast }}$ represents the range from $P^{\ast }$ to T.

If the endpoint of the ellipse with the major axis closer to R is $C_1$, then $|C_{1}R|$ can be expressed as:

$$|C_{1}R|=(2a-\left|RT\right|)/2$$

(36)

where $\left|RT\right|=\sqrt{{(x_T-x_R)}^2+{(y_T-y_R)}^2+{(z_T-z_R)}^2}$.

After that, we can obtain $|C_{1}O|=|C_{1}R|+1/2|RT|$, and the slope $k_{RT}$ of the $RT$ in O-UV coordinates $k_{RT}=(z_R-z_T)/(y_R-y_T)$. According to $k_{RT}$, the coordinates of $C_1$ under the O-UV coordinate system can be calculated as:

$$(u_{C_1},v_{C_1})=\left\{\begin{array}{c}(0,|C_{1}O\left|\right),k_{RT}=Inf\hfill \\ \left(\right|C_{1}O|\ast cos(arctan\left(k_{RT}\right)),|C_{1}O|\ast \sin \left(\arctan \left(k_{RT}\right)\right)),other\hfill \end{array}\right.$$

(37)

The projection ellipse in the O-UV coordinate system takes the form of a general ellipse, which can be expressed as ${\left(mu\right)}^2+{\left(nv\right)}^2=1$. Given that $C_1$$(u_{C1},v_{C1})$ and $P^{\ast }$$(u_{P^{\ast }},v_{P^{\ast }})$ are two points on an ellipse, the values of m and n can be obtained by substituting them into the above equation.

Under far-field observation conditions, it can be considered that the slope at $P^{\ast }$ is positive in the O-UV coordinate system. From the general elliptic equation, $dv/du$ can be calculated as $dv/du=-mu/nv$. So, the slope $k_{P^{\ast }}$ and ${\alpha }_2$ can be expressed as:

$$k_{P^{\ast }}=-{\displaystyle \frac{m{u}_{P^{\ast }}}{n{v}_{P^{\ast }}}}$$

(38)

$${\alpha }_2=arctan{k}_{P^{\ast }}$$

(39)

5.3. Computational Complexity

The main difference between the PGA method and the BF-PGA method is the method used to calculate the PS height. The PGA method uses the discrete Fourier transform method to calculate the PS height, usually implemented using the fast Fourier transform. Its overall computational complexity can be denoted as:

$$O\left(QNlogN\right)$$

(40)

where Q is the number of PS points that have been filtered, and N is the number of flight trajectories.

For the BF-PGA method, the beamforming algorithm is used to calculate the PS height, and the specific calculation formula can be expressed as:

$${\widehat{P}}_{BF}\left(s\right)={\displaystyle \frac{{\mathbf{a}}^H\left(s\right)\widehat{\mathbf{R}}\mathbf{a}\left(s\right)}{{\mathbf{a}}^H\left(s\right)\mathbf{a}\left(s\right)}}$$

(41)

Considering the same number of PS points, the beamforming method performs spatial spectral search using Equation (41), with a corresponding computational complexity of

$$O\left(Q\right(N+1\left)NS\right)$$

(42)

where S denotes the number of samples alone elevation and S usually takes a value greater than N.

For the ISOA method, the main calculation is the optimization processing step for estimating the phase error of each pixel point, and its computational complexity is

$$O\left(XY{N}^3\right)$$

(43)

where X and Y are the pixel sizes of the 2D image, respectively.

For the proposed method, the main calculation steps are the height calculation process of PGA and the phase error calculation process of each pixel point of EB-ISOA. The time consumption of PGA is almost negligible, so the computational complexity of the proposed method is consistent with that of Equation (43).

The average computation time of Monte Carlo experiments at SNR = 5 dB is selected to investigate the time consumption of several algorithms. The Q size in Equation (42) is about 450, N is 25, S is 25. X and Y are 81 and 161 in Equation (43), respectively. The average processing time of the PGA method is 0.029 s, the average processing time of the BF-PGA method is 1.567 s, the average processing time of the ISOA method is 108.722 s, and the average processing time of the proposed method is 104.458 s. The time consumption ratio of several phase correction methods conforms to the computational complexity given earlier. The method proposed in this paper, which performs phase error correction pixel-by-pixel, has a longer computation time compared to PGA based methods, but the phase correction effect is better.

6. Conclusions

To address the issue of the temporal decorrelation in the MB MonoSAR imaging, this paper investigates a MB OS BiSAR imaging configuration, which utilizes the imaging mode with one-transmit-multiple-receive to reduce the temporal baseline. First, the mathematical derivation of the signal model for the proposed MB OS BiSAR is conducted, yielding the signal form for de-ramping. Then, the 3D imaging process for this configuration is also presented. Furthermore, to address the phase error contamination caused by the platform jitter and inaccurate position measurements, a 3D imaging method based on the phase error calibration combining the PGA and EB-ISOA is proposed. The method first uses the ADI criterion to select PS, then applies the PGA method for the coarse phase error calibration, followed by the EB-ISOA for the precise phase error calibration of the coarse result. Finally, high-precision imaging results are obtained through the tomographic focusing and coordinate transformation.

To verify the practicality of the studied MB BiSAR configuration and the effectiveness of the proposed method, two simulation experiments are conducted. The first simulation experiment uses point target imaging results and three metrics to demonstrate the 3D imaging capability of the studied configuration. The second simulation experiment compares the proposed method with three phase error calibration methods. Qualitative evaluation of the imaging results from the simulation data is performed, along with the quantitative assessment based on entropy and contrast metrics. The results show that, compared to other comparison methods, the proposed method not only solves the issue of vertical displacement, but also enables the pixel-by-pixel calibration. The proposed method overcomes the limitation of the PGA method and the problem of weak target focus failure in ISOA, yielding clearer 3D BiSAR imaging results.

In the future work, we will reduce the computational complexity of the proposed algorithm and explore dividing the scene into blocks for the phase error estimation. Additionally, we will use the polynomial fitting to obtain approximate estimates to reduce processing time, enabling the method to be applied to more scenarios, such as large-scale urban and mountainous areas.