A Modified Range Doppler Algorithm for High-Squint SAR Data Imaging

Abstract

The high-squint airborne Synthetic Aperture Radar (SAR) has the ability to detect the target area flexibly, and the detection swath is significantly increased compared with the side-looking SAR system. Therefore, it is of great significance to carry out research on high-precision imaging methods for high-squint airborne SAR. However, the high-squint SAR echoes have large Range Cell Migration (RCM), resulting in severe range–azimuth coupling and strong spatial variation. In this paper, a Modified Range Doppler Algorithm (MRDA) is proposed to cope with these effects introduced by the significant RCM in high-squint airborne SAR imaging. The bulk compensation preprocessing is first adopted to remove the considerable RCM and severe cross-coupling in a two-dimensional frequency domain. Then, Non-Linear Chirp Scaling (NLCS) in the range direction is utilized to equalize the range-variant chirp rate caused by the residual RCM and coupling and, therefore, the consistent range phase compensation can be fulfilled in range frequency domain. The modified correlation processing is executed to compensate the residual Doppler phase modulation, the residual RCM and the range-variant cubic phase modulation, which guarantees the characteristics of high efficiency and high precision. The simulations have demonstrated that the MRDA can focus the SAR echoes with large squint angles more effectively than the algorithms based on the Linear Range Walk Correction (LRWC) method.

1. Introduction

As a typical active microwave remote sensing technology, different from optical satellite systems, Synthetic Aperture Radar (SAR) has all-day and all-weather working capabilities [1,2]. SAR is an imaging radar that can provide two-dimensional high-resolution images, which are widely used in various fields, such as environmental protection, resource exploration, disaster monitoring, battlefield environment reconnaissance, target precision attack, etc., and play an important role in many civil and military fields. Since the concept of SAR was first proposed by Carl Wiley of Goodyear Aircraft Corporation in the 1950s, this field has developed in the directions of high resolution, multi-mode, multi-polarization and multi-platform [3,4]. The traditional SAR system operates in broadside mode, which means the radar Line of Sight (LOS) is perpendicular to the platform’s flight path. When the system works in this mode, it can only observe and image the broadside area of the platform, and therefore flexible observation and information acquisition for complex environments cannot be achieved. Normally, the angle between the LOS direction and the vertical direction of the flight path is defined as the squint angle [2]. The large squint angle SAR platform allows the LOS direction to vary in a large scope, which can break through the above limitations and realize flexible multi-angle detection. Additionally, the SAR system with large squint angle can also acquire the multi-angle scattering characteristics of the target, so as to realize the precise identification of the key targets [2,5,6]. Therefore, high-squint airborne SAR has significant application value in many fields.

Due to the significant Range Cell Migration (RCM) and severe range–azimuth cross-coupling of high-squint SAR echoes, many existing image-processing algorithms cannot be applied to large squint angle and large scene cases [5,6,7,8,9]. How to realize RCM Correction (RCMC) and cross-coupling removal for high-squint airborne SAR is the premise of high-precision imaging. According to the implement domain of the algorithm, the existing imaging algorithms applied to the squint SAR system can be divided into two categories: time-domain and frequency-domain imaging algorithms. Theoretically, the time-domain imaging algorithm has the ability of high-resolution imaging and can be used in almost any mode and configuration of the SAR system. The typical time-domain imaging algorithm is represented by the Back Projection Algorithm (BPA) [10]. Although the BP algorithm has high precision and excellent adaptability, its computational complexity in the literature [10] is O(N³), which indicates its low computational efficiency. Aiming at the calculation efficiency problem of the BP algorithm, various improved algorithms have been proposed, as shown below. In [11], the Local BP Algorithm (LBPA) is proposed, and it reduces the complexity of the BP algorithm from O(N³) to O(N^2.5). The Factorization Fast Back Projection Algorithm (FFBPA) is presented in [12], the complexity of which is O(N²log₂N). As a classic BP-based efficient algorithm, it is widely used in practice. However, the efficiency improvements of the typical BP algorithm and its extended algorithms above are still rather limited in real-time imaging, especially in large-scale scene imaging.

The frequency-domain imaging algorithms can make full use of the “azimuth-invariant” property and improve the computational efficiency through batch processing. As the primary imaging algorithms at present, the significant RCM and severe coupling of the high-squint SAR echoes have hindered the application of the frequency-domain algorithms. Based on the classic Range Doppler Algorithm (RDA) in [13], an improved RD algorithm with Secondary Range Compression (SRC) operation [14] is proposed to realize SAR image-processing at small squint angles, which is not suitable for SAR imaging at large squint angles. A two-dimensional frequency algorithm, for example, the ωK Algorithm (ωKA), can accomplish squint SAR imaging under the condition of narrow imaging swath or poor resolution [15,16,17,18]. But when the squint angle further increases, the imaging algorithm is still incapable. This limitation is mainly attributed to the fact that the imaging parameters cannot be updated throughout the range cells. In addition, the Stolt interpolation required in the ωKA will have a great side effect on SAR imaging efficiency. The Polar Format Algorithm (PFA) is another effective approach to cope with the large squint angle SAR imaging problem. A PFA-based Two-Stage Focusing Algorithm (TSFA) combining coarse focusing and fine focusing procedures is proposed in [19] to correct the significant RCM and severe range–azimuth decoupling via a “squint minimization” operation. It should be noted that the PFA is generally used in the SAR spotlight mode, and the imaging results are obtained by polar coordinate projection and transformation, which require a large amount of calculation and therefore limit its application. Affected by the wavefront curvature, the traditional PFA cannot meet the imaging requirement of a large scene. Although the new PFA-based algorithm in [20] can expand the imaging size, the two-step PFA processing and the block processing bring a greater computational burden. As another frequency-domain algorithm, the Chirp Scaling Algorithm (CSA) [21,22] plays a key role in the research of squint SAR imaging algorithms. The “squint minimization” preprocessing is first proposed in [23] to realize the echo RCMC and cross-decoupling, after which operation the Azimuth Non-Linear Chirp Scaling (ANLCS) method is used to resolve the azimuth variation problem introduced by the preprocessing. In [24], an extended two-step focusing algorithm is proposed, which solves the two-dimensional spectrum distortion problem caused by the squint angles and removes the two-dimensional coupling by means of a “non-linear shift” method. The “squint minimization” and “non-linear shift” methods are collectively referred to as Linear Range Walk Correction (LRWC) in the literature [25]. An Extended Non-Linear Chirp Scaling Algorithm (ENLCSA) is given to equalize the azimuth-variant frequency modulation rate and the higher-order terms are also considered in [25]. Essentially, the aforementioned LRWC-based algorithms result in range cell shifts for azimuthal non-central targets while eliminating the significant RCM and severe cross-coupling, which destroys the “azimuth-invariant” property. Since this azimuth variation caused by the ANLCS processing cannot be completely removed, the azimuth swath will be significantly reduced. Therefore, the LRWC-based algorithms in [19,23,24,25] cannot adapt to large scene imaging processing, especially for large squint angles cases.

In this paper, a Modified Range Doppler Algorithm (MRDA) is proposed for the high-squint airborne SAR system. In order to eliminate the significant RCM and severe range–azimuth cross-coupling of high-squint SAR echoes, this novel algorithm begins with bulk compensation preprocessing, which is performed in the two-dimensional frequency domain. Unlike the LRWC method in [19,23,24,25], the bulk compensation preprocessing does not introduce range cell shifts for the echo signals, and thus the “azimuth-invariant” property can be preserved, which means the batch focusing processing can be accomplished in the azimuth frequency domain. After the bulk compensation preprocessing, only the residual RCM and the residual cross-coupling can be observed in the echo signals. The Non-Linear Chirp Scaling (NLCS) method is then adopted to equalize the range-variant frequency modulation rate of the range pulses caused by the residual cross-coupling. The consistent range compression can subsequently be realized by the range phase compensation function in range frequency domain. The modified correlation processing is then be applied to precisely remove the residual RCM, the residual Doppler phase modulation and the range-variant cubic phase modulation introduced by the residual cross-coupling. The characteristics of the modified correlation reference function can be summarized as follows: firstly, as the bulk RCM has been corrected by the bulk compensation preprocessing and the variation of the quadratic range frequency term has been eliminated by the NLCS method, the correlation window length is significantly shortened, which ensures processing efficiency. Secondly, the construction of the reference function takes into account the range-variant cubic phase modulation introduced by the residual cross-coupling. Therefore, the modified correlation processing can correct the asymmetry of the range side-lobes while completing the residual RCMC and the residual Doppler phase compensation, which can realize high-precision processing. Ultimately, the final focused SAR image data can be obtained via an Azimuth Inverse Fast Fourier Transform (AIFFT) operation. The feasibility of the proposed MRDA is then verified by comparative simulations.

The rest of this paper is organized as follows: the geometry and signal models are displayed in Section 2. The details of the proposed MRDA for a high-squint airborne SAR system are shown in Section 3. The simulations and results are given in Section 4. A further discussion and conclusion are present in Section 5 and Section 6, respectively.

2. Geometry and Signal Models

The geometry model is provided in Figure 1 for the high-squint airborne SAR system. The platform flies along the y-axis with a fixed velocity v. In this paper, the squint angle α is defined as the complementary angle between the line-of-sight direction and the flight path of the platform. Point P is an arbitrary target in the scene. R₀ is the closest slant range between the sensor and the point target P, and η = 0 represents the corresponding azimuth time when the sensor passes the closest slant range R₀. The range history R(η) of point target P can then be expressed as

$$R\left(\eta \right)=\sqrt{R_0^2+v^2{\eta }^2}$$

(1)

The baseband signal for a single point target after demodulation has the following complex form:

$$s\left(\tau ,\eta \right)=w_a\left(\eta -{\eta }_c\right)w_r\left[\tau -\frac{2R\left(\eta \right)}{c}\right]\exp \left[-j4\pi \frac{R\left(\eta \right)}{\lambda }\right]\exp \left\{j\pi k_r{\left[\tau -\frac{2R\left(\eta \right)}{c}\right]}^2\right\}$$

(2)

where τ and η are the range fast time and azimuth slow time, respectively. w_a(∙) and w_r(∙) are the antenna patterns in azimuth and range directions, respectively. η_c represents the beam center-crossing time in the azimuth direction, λ is the wavelength of the transmitted signal, c is the speed of light in space, k_r is the range frequency modulation rate of the transmitted range pulse, R(η) is the instant slant range in different azimuth times as denoted in the range of Equation (1).

Using the Principle of Stationary Phase (POSP) with Equation (2), we can obtain

$$S_{2df}\left(f_{\tau },f_{\eta }\right)=W_r\left(f_{\tau }\right)W_a\left(f_{\eta }-f_{\eta c}\right)\exp \left[-j\frac{4\pi R_0\left(f_0+f_{\tau }\right)}{c}{D}_{2df}\right]\exp \left(-j\pi \frac{f_{\tau }^2}{k_r}\right)$$

(3)

where f_τ and f_η represent the frequencies in range and azimuth directions, respectively. f_ηc is the Doppler centroid frequency for a certain point target, and f₀ is the carrier frequency of the transmitted signal. W_a(∙) and W_r(∙) are the envelopes of the signal spectrum in azimuth and range directions, respectively. D_2df is the cosine value of the instantaneous squint angle at the azimuth time η in the straight-line geometry shown in Figure 1 and is also considered as the migration factor in the two-dimensional frequency domain. The specific expression of D_2df is

$$D_{2df}=\sqrt{1-\frac{c^2{f}_{\eta }^2}{4v^2{\left(f_0+f_{\tau }\right)}^2}}$$

(4)

3. Modified Range Doppler Algorithm

Based on the geometry and signal models shown in Section 2, the MRDA is proposed in Section 3 for a high-squint airborne SAR system. The block diagram of the proposed MRDA is shown in Figure 2 below. Firstly, the raw echo data are transformed into the two-dimensional frequency domain, and the bulk compensation preprocessing is applied to remove the range–azimuth cross-coupling and the total RCM of the reference range cell without introducing the azimuth-variation effect. Then, the NLCS processing is performed along the range direction in the Range Doppler (RD) domain to equalize the range-variant frequency modulation rate introduced by the residual cross-coupling. Subsequently, the range matched filter can be realized via the uniform range phase compensation, and the compressed range pulses can therefore be obtained. Finally, the modified correlation processing in the RD domain is adopted to accomplish the residual RCMC, the residual Doppler phase modulation compensation and the range-variant cubic phase modulation compensation, and the focused SAR data can be acquired at the end via an AIFFT operation.

3.1. Bulk Compensation Preprocessing

In order to remove the range–azimuth coupling, the bulk compensation filter is constructed in the two-dimensional frequency domain and is given the form of

$$H_{Bulk}\left(f_{\tau },f_{\eta }\right)=\exp \left[j\frac{4\pi R_{ref}\left(f_0+f_{\tau }\right)}{c}{D}_{2df}\right]$$

(5)

where R_ref is the referent closest slant range in the central range gate. By multiplying Equation (3) with (5), the signal after bulk compensation preprocessing in the two-dimensional frequency domain can be written as

$$S_{res}\left(f_{\tau },f_{\eta }\right)=W_r\left(f_{\tau }\right)W_a\left(f_{\eta }-f_{\eta c}\right)\exp \left[j{\theta }_{res}\left(f_{\tau },f_{\eta }\right)\right]$$

(6)

$${\theta }_{res}\left(f_{\tau },f_{\eta }\right)=-\frac{4\pi \left(R_0-R_{ref}\right)\left(f_0+f_{\tau }\right)}{c}{D}_{2df}-\pi \frac{f_{\tau }^2}{k_r}$$

(7)

A simulation using the parameters listed in

Table 1. List of simulation parameters.

Description	Value	Units
Radar height	20	km
Radar velocity	200	m/s
Wavelength	0.03	m
Pulse width	30	μs
Bandwidth	150	MHz
Sample frequency	180	MHz
PRF	300	MHZ
Look angle	60	deg
Squint angle	45	deg
Antenna length	2	m

is conducted to explain the effect of the bulk compensation preprocessing. Figure 3a,b depict the raw echo signals in the two-dimensional time domain after bulk compensation for central and non-central range cells, respectively. Three point targets are positioned at the same range gate with a constant azimuth interval of 1 km. The range position displacement is set to be 3 km for these three point targets, shown in Figure 3b. After bulk compensation preprocessing, the total RCM, Doppler phase and the range–azimuth coupling for the target positioned at the central range cell are removed as is shown in Figure 3a. The signal for the central point target has been zoomed into, shown by the red solid line box. However, the slight skew of the signal along the horizontal direction and the extension along the vertical direction indicate that there are still residual RCM and residual Doppler phase modulation, respectively, for targets located in non-central range cells. Moreover, this extension can also account for the residual cross-coupling. Therefore, it can be concluded from the above simulation results that the bulk compensation preprocessing can eliminate the total RCM of the reference range cell while maintaining the “azimuth-invariant” property and realize the range–azimuth decoupling processing.

3.2. Non-Linear Chirp Scaling

Using Taylor expansion and keeping Equation (7) to the third order, we can obtain

$$\begin{array}{cc}\hfill {\theta }_{res}\left(f_{\tau },f_{\eta }\right)& =-\frac{4\pi \left(R_0-R_{ref}\right)}{\lambda }\left[D_f+\frac{1}{D_f}\frac{f_{\tau }}{f_0}-\frac{{\left(\frac{\lambda f_{\eta }}{2v}\right)}^2}{2D_f^3}{\left(\frac{f_{\tau }}{f_0}\right)}^2+\frac{{\left(\frac{\lambda f_{\eta }}{2v}\right)}^2}{2D_f^5}{\left(\frac{f_{\tau }}{f_0}\right)}^3\right]-\pi \frac{f_{\tau }^2}{k_r}\hfill \\ \hfill & =-\frac{4\pi \left(R_0-R_{ref}\right)D_f}{\lambda }-\frac{4\pi \left(R_0-R_{ref}\right)}{c{D}_f}{f}_{\tau }+\pi \left(z_1-\frac{1}{k_r}\right)f_{\tau }^2+\pi z_2{f}_{\tau }^3\hfill \end{array}$$

(8)

$$z_1=\frac{c{f}_{\eta }^2}{2v^2{f}_0^3{D}_f^3}\left(R_0-R_{ref}\right)$$

(9)

$$z_2=-\frac{c{f}_{\eta }^2}{2v^2{f}_0^4{D}_f^5}\left(R_0-R_{ref}\right)$$

(10)

$$D_f=\sqrt{1-\frac{c^2{f}_{\eta }^2}{4v^2{f}_0^2}}$$

(11)

where D_f is the range cell migration factor in the RD domain, which represents the cosine value of the instantaneous squint angle corresponding to different azimuth positions of a certain range gate.

After applying Taylor expansion to Equation (7), the uniformly compensated echo signal is further transformed into the RD domain using the POSP:

$$S_{rd}\left(\tau ,f_{\eta }\right)={\displaystyle {\int }_{-\infty }^{\infty }{S}_{res}\left(f_{\tau },f_{\eta }\right)\exp \left(j2\pi f_{\tau }\tau \right)d{f}_{\tau }}=w_a\left(f_{\eta }-f_{\eta c}\right){\displaystyle {\int }_{-\infty }^{\infty }{w}_r\left(f_{\tau }\right)\exp \left[j\theta \left(f_{\tau },f_{\eta }\right)\right]d{f}_{\tau }}$$

(12)

$$\theta \left(f_{\tau },f_{\eta }\right)={\theta }_{res}\left(f_{\tau },f_{\eta }\right)+2\pi f_{\tau }\tau$$

(13)

Neglecting the cubic term of f_τ in Equation (8), the stationary phase point can be obtained by the following equation:

$$\frac{d\theta \left(f_{\tau },f_{\eta }\right)}{d{f}_{\tau }}=-\frac{4\pi R_{res}}{c{D}_f}+2\pi \left(z_1-\frac{1}{k_r}\right)f_{\tau }+2\pi \tau$$

(14)

where R_res is the relative closest range and can be defined as

$$R_{res}=R_0-R_{ref}$$

(15)

When Equation (14) is set to be zero, the range frequency f_τ and fast time τ can be related by

$$f_{\tau }=k_m\left(\tau -\frac{2R_{res}}{c{D}_f}\right)$$

(16)

$$k_m=\frac{k_r}{1-k_r{z}_1}$$

(17)

The RD spectrum after bulk compensation preprocessing can be given as

$$\begin{array}{cc}\hfill S_{rd}\left(\tau ,f_{\eta }\right)& =w_r\left[\frac{k_m}{k_r}\left(\tau -\frac{2R_{res}}{c{D}_f}\right)\right]w_a\left(f_{\eta }-f_{\eta c}\right)\exp \left(-j\frac{4\pi R_{res}{D}_f}{\lambda }\right)\hfill \\ \hfill & \cdot \exp \left\{j\pi k_m{\left(\tau -\frac{2R_{res}}{c{D}_f}\right)}^2\right\}\exp \left\{j\pi z_2{\left[k_m\left(\tau -\frac{2R_{res}}{c{D}_f}\right)\right]}^3\right\}\hfill \end{array}$$

(18)

where the range antenna pattern w_r(·) reveals the residual RCM of the echo signal. The magnitude and phase of the azimuth signal are represented by the second and third terms in Equation (18), respectively. The fourth term and the fifth term successively represent the quadratic and cubic range frequency modulation terms of the range signal. It can be determined from the above analysis that the updated range frequency modulation rate and the residual RCM are both variant in the range dimensional. Since Equation (18) satisfies the condition of |k_rz₁| << 1, the following approximation can be obtained:

$$k_m\approx k_r\cdot \left(1+k_r{z}_1\right)=k_r+k_{m1}$$

(19)

$$k_{m1}=k_{m1}\left(f_{\eta },R_{res}\right)=k_r^2{z}_1$$

(20)

According to Equation (20) above, k_m1 is the range-variant part of the updated range frequency modulation rate k_m, which is related to both the azimuth frequency f_η and the relative closest range R_res. Actually, the variation of k_m1 with R_res is the dominant term, so the variation of k_m1 along the azimuth is ignored. The Non-Linear Chirp Scaling (NLCS) filter is then given a cubic form of fast time τ in the range direction to remove the range-variant component of range frequency modulation rate k_m.

$$H_{nlcs}\left(\tau ,f_{\eta }\right)=\exp \left(-j\pi \gamma {\tau }^3\right)$$

(21)

where

$$\gamma ={\frac{k_{m1}c{D}_f}{6R_{res}}|}_{f_{\eta }=f_{dc}}=\frac{c^2{k}_r^2{f}_{dc}^2}{12v^2{f}_0^3{D}_{fref}^2}$$

(22)

$$D_{fref}={D_f|}_{f_{\eta }=f_{dc}}$$

(23)

The influence of the range-variant linear term introduced by the NLCS operation on the matched filter will be taken into account in the modified correlation processing. The quadratic term is useful for eliminating the variation of k_m in the range direction. It should be noted that the cubic term introduced by the NLCS operation is uniform for different range cells and can be compensated consistently by the range matched filtering.

After utilizing the NLCS operation, the echo signal in Equation (18) is changed into

$$\begin{array}{cc}\hfill S_{rd1}\left(\tau ,f_{\eta }\right)& =w_r\left(\tau -\frac{2R_{res}}{c{D}_f}\right)W_a\left(f_{\eta }-f_{\eta c}\right)\exp \left(-j\frac{4\pi R_{res}{D}_f}{\lambda }\right)\hfill \\ \hfill & \cdot \exp \left[j\pi k_r{\left(\tau -\frac{2R_{res}}{c{D}_f}\right)}^2\right]\exp \left[j\pi z_2{k}_r^3{\left(\tau -\frac{2R_{res}}{c{D}_f}\right)}^3\right]\hfill \end{array}$$

(24)

So far, the compression processing in the range direction can be completed by the batching processing in the range frequency domain.

3.3. Range Phase Compensation

In order to realize range compression, the signal in Equation (24) is then transformed into the two-dimensional frequency domain by Range FFT (RFFT):

$$\begin{array}{cc}\hfill S_{2df1}\left(f_{\tau },f_{\eta }\right)& ={\displaystyle {\int }_{-\infty }^{\infty }{S}_{rd1}\left(\tau ,f_{\eta }\right)\exp \left(-j2\pi f_{\tau }\tau \right)d\tau }\hfill \\ \hfill & =W_r\left(f_{\tau }\right)W_a\left(f_{\eta }-f_{\eta c}\right)\exp \left(-j\frac{4\pi R_{res}{D}_f}{\lambda }\right)\hfill \\ \hfill & \cdot \exp \left(-j\pi \frac{f_{\tau }^2}{k_r}\right)\exp \left(j\pi z_2{f}_{\tau }^3\right)\exp \left(-j2\pi f_{\tau }\frac{2R_{res}}{c{D}_f}\right)\hfill \end{array}$$

(25)

where the second exponential term is the consistent quadratic range frequency modulation phase. The third and fourth exponential terms represent the range-variant cubic frequency modulation phase and the residual RCM, respectively, the compensation of which can be realized via the followed modified correlation processing. Therefore, the expression of the range compression filter is

$$H_{rc}\left(f_{\tau },f_{\eta }\right)=\exp \left[j\pi \frac{f_{\tau }^2}{k_r}\right]$$

(26)

Additionally, the cubic phase introduced by the NLCS filter also should be taken into consideration, and this cubic phase compensation function is given the form of

$$H_{cpc}\left(f_{\tau },f_{\eta }\right)=\exp \left\{j\pi \gamma {\left[\tau \left(f_{\tau }\right)-\frac{2R_{res}}{c{D}_f}\right]}^3\right\}=\exp \left\{j\pi \gamma {\left(\frac{f_{\tau }}{k_r}\right)}^3\right\}$$

(27)

where τ(f_τ) is the stationary phase point of RFFT and can be expressed as

$$\tau \left(f_{\tau }\right)=\frac{f_{\tau }}{k_r}+\frac{2R_{res}}{c{D}_f}$$

(28)

The range phase compensation filter can therefore be integrated as

$$H_{rpc}\left(f_{\tau },f_{\eta }\right)=H_{rc}\left(f_{\tau },f_{\eta }\right)\cdot H_{cpc}\left(f_{\tau },f_{\eta }\right)=\exp \left[j\pi \frac{f_{\tau }^2}{k_r}+j\pi \gamma {\left(\frac{f_{\tau }}{k_r}\right)}^3\right]$$

(29)

The two-dimensional frequency signal after range phase compensation via Equation (29) becomes

$$S_{2df2}\left(f_{\tau },f_{\eta }\right)=W_r\left(f_{\tau }\right)W_a\left(f_{\eta }-f_{\eta c}\right)\exp \left(-j\frac{4\pi R_{res}{D}_f}{\lambda }\right)\exp \left(j\pi z_2{f}_{\tau }^3\right)\exp \left(-j2\pi f_{\tau }\frac{2R_{res}}{c{D}_f}\right)$$

(30)

The three exponential terms in the above Equation (30) represent the range-variant azimuth phase modulation, cubic phase modulation of range frequency and residual RCM, respectively.

3.4. Modified Correlation Processing

After the Range Inverse FFT (RIFFT) is applied to Equation (30), the signal is again transformed into the RD domain:

$$S_{rd2}\left(\tau ,f_{\eta };R_0\right)=r\left(\tau -\frac{2R_{res}}{c{D}_f},f_{\eta };R_0\right)W_a\left(f_{\eta }-f_{\eta c}\right)\exp \left(-j\frac{4\pi R_{res}{D}_f}{\lambda }\right)$$

(31)

$$r\left(\tau ,f_{\eta };R_0\right)={\displaystyle {\int }_{-\infty }^{\infty }{W}_r\left(f_{\tau }\right)\exp \left(j\pi z_2{f}_{\tau }^3\right)\exp \left(j2\pi f_{\tau }\tau \right)d{f}_{\tau }}={ℱ}_r^{-1}\left[W_r\left(f_{\tau }\right)\exp \left(j\pi z_2{f}_{\tau }^3\right)\right]$$

(32)

where $ℱ$_r⁻¹(∙) denotes the RIFFT operation. Figure 4a shows the range compressed result before NLCS processing in the range direction, which reflects the range pulse broadening due to the residual cross-coupling for the range marginal target. It can be observed from Section 3.2 that the variation of the range frequency modulation rate caused by the residual cross-coupling has been eliminated, which allows the signal to be consistently compressed by the range phase compensation in Section 3.3. The range compressed result after NLCS processing is depicted in Figure 4b, and the width of the range signal in the fast time domain is significantly reduced, which indicates that the window width of the correlation processing will be remarkably shortened for large squint angle cases. Additionally, Equation (32) shows that there is an additional range-variant cubic phase modulation of the rectangular window function in the range direction, which is also introduced by the residual cross-coupling. Different from the normal sinc function, the compressed range pulse r(τ, f_η; R₀) therefore has an asymmetric form as is shown in the partially enlarged view of Figure 4b, where obvious side-lobe energy asymmetry appears. This is the key point to distinguish the modified correlation processing from the traditional ones.

To compensate the range-variant cubic phase modulation in Equation (32), we first quantify the maximum value of this cubic phase at different azimuth frequencies for different range cells. Assuming that the cubic phase modulation quantization scope is [−φ, φ], the quantization interval is 2φ/(N₁ − 1) when the corresponding quantization point is N₁, which is shown in the height dimensional in Figure 5. For the n₁th quantization interval, the corresponding cubic phase modulation ϕ_n1 is

$${\varphi }_{n1}=-\phi +\frac{2\phi }{N_1-1}\left(n_1-1\right),\begin{array}{cc}& n_1=1,\dots ,N_1\end{array}$$

(33)

The compressed range pulse with phase modulation ϕ_n1 can be expressed as

$$g\left(\tau ,{\varphi }_{n1}\right)={ℱ}_r^{-1}\left[W_r\left(f_{\tau }\right)\exp \left(j{\varphi }_{n1}\right)\right]$$

(34)

where n₁ = 1, …, N₁. Figure 5a,b show the compressed range pulse g(τ, ϕ_n1) with normalized amplitude and dB-scaled normalized amplitude when phase modulation ϕ_n1 = π/4, respectively. Kernel functions with quantization shifts can be obtained by quantizing g(τ, ϕ_n1) at equal intervals, and thus a two-dimensional matrix of kernel function is constructed as is indicated by the yellow boxes in Figure 6. N₂ and N₃ denote the length of the correlation window and the number of quantization shifts for RCMC, respectively. The above analyses illustrate that the kernel function under different cubic phase modulations in the modified correlation processing is a three-dimensional matrix, which can compensate the range-variant cubic phase modulation while completing the residual RCMC.

The modified correlation reference function can then be given as

$$H_a\left(\tau ,f_{\eta };R_0\right)=g\left(\tau ,{\varphi }_{n1}\right)\exp \left(-j\frac{4\pi R_{res}{D}_f}{\lambda }\right)$$

(35)

In order to achieve accurate focusing of non-central range cells in high-squint observation, the modified correlation processing is adopted, which can be described as

$$S_{rd3}\left(\tau ,f_{\eta };R_0\right)={\displaystyle \sum _{n_2=-\frac{N_2}{2}+1}^{\frac{N_2}{2}}{S}_{rd2}\left(\tau -\frac{\lfloor {\tau }_0^{\prime }\cdot f_s\rfloor }{f_s}+n_2,f_{\eta }\right)}\cdot H_a\left(\tau -\frac{\lfloor {\tau }_0^{\prime }\cdot f_s\rfloor }{f_s},f_{\eta };R_0\right)$$

(36)

where

$${\tau }_0=\frac{2R_{res}}{c{D}_f}$$

(37)

$${\tau }_0^{\prime }={\tau }_0+\frac{3\gamma {\tau }_0^2}{2k_r}$$

(38)

τ₀ and τ′₀ are the residual range-variant RCM after bulk compensation preprocessing and after NLCS processing, respectively. The second term in Equation (38) is the extra RCM introduced by the range-variant linear term of the NLCS operation in Section 3.2.

Through the modified correlation processing, the residual RCM, the residual Doppler phase modulation and the range-variant cubic phase modulation introduced by the residual cross-coupling are removed, and the focused SAR image can finally be obtained by

$$s\left(\tau ,\eta ,;R_0\right)={ℱ}_a^{-1}\left[S_{rd3}\left(\tau ,f_{\eta };R_0\right)\right]=p\left(\eta -{\eta }_c\right)p\left(\tau -{\tau }_c\right)$$

(39)

where $ℱ$_a⁻¹(∙) represents the AIFFT operation and p(∙) denotes the sinc function. According to Equations (37) and (38), τ_c can be expressed as

$${\tau }_c={{\tau }_0^{\prime }|}_{f_{\eta }=f_{dc}}$$

(40)

Different from simulated examples, the actual system is affected by various noises. The impact of the Signal-to-Noise Ratio (SNR) on the algorithm is mainly reflected in the estimation of Doppler parameters in real data processing, which limits the estimation accuracy of Doppler parameters. In order to fulfill precise focusing processing, especially in practical real-time imaging systems, it is necessary to ensure the accuracy and efficiency of Doppler parameter estimation.

4. Simulations and Results

In order to validate the feasibility of the proposed MRDA for the high-squint airborne SAR imaging, point target simulations are carried out in this section employing the parameters listed in Table 1. As depicted in Figure 7, a 5 × 5 point matrix is evenly arranged in a square area of 10 km × 10 km in the ground scene, and the distance interval is 2.5 km both in range and azimuth directions. The squint angle α is set to be 45°. Rectangular antenna patterns are assumed in both these two directions.

The simulation results adopting the proposed MRDA are displayed in Figure 8a and Figure 9. The focused SAR image of the point matrix is shown in Figure 8a, and the slant range plane has not been projected onto the ground range plane. The target energy is displayed at the dB scale to facilitate the identification of each point target. The imaging results (two-dimensional contours, azimuth profiles and the range profiles) of the chosen point targets 1, 7, 13, 19 and 25 are shown in Figure 10a–e, respectively. For the point targets 7 and 19 positioned at 2.5 km away from the scene center, the imaging performance is quite good, and it is almost consistent with that of central point target 13 as can be seen in Figure 10b,d. From Figure 10a,e, however, slight distortion can be observed in the two-dimensional contour and range profile of the marginal point target 5 km away from the scene center. The following reasons lead to this distortion in the range dimensional. On the one hand, the variation of range frequency modulation rate k_m1 along the azimuth frequency is ignored during the NLCS processing. The bulk compensation preprocessing eliminates the major range–azimuth cross-coupling of the signal in non-central range cells, but the residual cross-coupling still exists, which will result in spatial variation of the range frequency modulation rate k_m1. Actually, k_m1 is not only related to the relative closest slant range R_res but also to the azimuth frequency f_η. However, k_m1 is considered to be weakly correlated with f_η. On the other hand, additional terms introduced by the NLCS operations also cause slight degradation. While equalizing the comprehensive range frequency modulation rate k_m, the NLCS processing brings about a range-variant linear term and identical quadratic term of range frequency f_τ, which can be compensated in the subsequent processing. The range position displacement caused by the linear term will lead to incomplete equalization, which means the existence of small residual range variation of k_m. To illustrate the performance of this imaging algorithm from a quantitative perspective, the evaluation indicators of the chosen point targets are given in

Table 2. Evaluation results of the chosen point targets using the proposed MRDA.

Target	Azimuth			Range
	ρa (m)	PSLR (dB)	ISLR (dB)	ρr (m)	PSLR (dB)	ISLR (dB)
1	1.01	−13.17	−10.60	1.06	−12.96	−10.56
7	1.00	−13.20	−10.43	1.05	−13.10	−10.38
13	1.01	−13.25	−10.39	1.06	−13.25	−10.35
19	1.00	−13.23	−10.41	1.05	−13.21	−10.43
25	1.00	−13.18	−10.58	1.05	−13.22	−10.58

, including the evaluation resolution, the Peak Side-Lobe Ratio (PSLR) and the Integrated Side-Lobe Ratio (ISLR) of profiles both in range and azimuth dimensions. When the simulation parameters listed in Table 1 are employed and the squint angle α is set to be 0°, the theoretical values of the slant range resolution and azimuth resolution are both 0.89 m. When the rectangle patterns are assumed both in range and azimuth directions, the theoretical values of the PSLR and the ISLR are −13.26 dB and −9.68 dB, respectively. It can be seen from Table 2 that all these point targets are well focused in the azimuth direction within the whole scene arranged in Figure 7—the azimuth PSLR deviates from the ideal value of −13.26 dB by less than 0.09 dB. Slight distortion occurs for the range profiles of the chosen point targets—the range PSLR has a degradation of less than 0.3 dB from the expected value. Overall, all chosen point targets are well focused and show great consistency, which verifies the effectiveness of the proposed MRDA.

By using the LRWC-based algorithm, another point target simulation is carried out to compare the performance of the proposed MRDA. The simulation parameter configuration is the same as that of the MRDA as shown in Table 1, whereas the scene parameter configuration is different, which is reflected in the fact that the point matrix intervals both in range and azimuth directions are 1 km. Therefore, the size of the imaging scene is set to be a square area of 4 km × 4 km in this comparison experiment. The imaging result of the point matrix when the LRWC-based algorithm is adopted is depicted in Figure 8b and the details of the imaging performances for the selected point targets are shown in Figure 10. It is worth noting that the focused SAR image in Figure 8b has been subjected to geometric correction processing, so the geometric distortion caused by the LRWC operation has been removed. Significant degradation can be observed from the two-dimensional contours and the azimuth profiles of point targets 1 and 25 shown in Figure 10. The azimuth phase error has even reached 3/4π, which far exceeds the phase error standard of 1/4π. Significant main-lobe extensions and side-lobe elevations can be seen in Figure 10a,e. In order to describe the degree of distortion caused by the LRWC-based method from a quantitative point of view,

Table 3. Evaluation results of the chosen point targets using the LRWC-based algorithm.

Target	Azimuth			Range
	ρa (m)	PSLR (dB)	ISLR (dB)	ρr (m)	PSLR (dB)	ISLR (dB)
1	1.20	−6.64	−4.97	1.06	−13.24	−10.69
7	1.04	−12.21	−9.91	1.05	−13.26	−10.70
13	1.00	−13.25	−10.50	1.05	−13.25	−10.70
19	1.00	−12.02	−9.80	1.05	−13.25	−10.71
25	1.17	−5.28	−4.02	1.05	−13.18	−10.72

presents the evaluation indicators of the selected point targets. Only the central point target 13 is well focused in the azimuth direction, and no deterioration occurs in azimuth evaluation indicators. Point targets 1 and 25 are the farthest from the center of the scene and hence suffer the severest quality degradation. From the analysis of azimuth resolution, Table 3 shows that the Impulse Response Widths (IRWs) in the azimuth direction of point targets 1 and 25 are 20% and 17% worse than that of the scene center point target 13, respectively, which indicates a clear deterioration in azimuth resolution. The PSLR and ISLR degradations in the azimuth direction for point targets 1 and 25 are approximately 7–8 dB and 5–6 dB, respectively. For the point targets whose displacements in azimuth and range directions are both 1 km, the azimuth PSLR degradation in the azimuth direction is more than 1.05 dB from the theoretical value. Practically, the above distortions are attributed to complicated azimuth variation introduced by the LRWC operation. For point targets located at different azimuth positions of the same range gate, different range cell shifts occur, which leads to the failure of the “azimuth-invariant” property. Since severe spatial variation of Doppler parameters along the azimuth direction, the frequency-domain algorithm can only guarantee target focus within a small azimuth swath. As shown in the simulation results shown in Figure 10 and Table 3, this azimuth variation cannot be completely eliminated even when employing the azimuth NLCS processing after the LRWC operation. By comparison with the proposed MRDA, the MRDA shows obvious superiority in focusing performance.

5. Discussion

5.1. Phase Error Analysis

From the analyses in Section 3.2, it can be seen that k_m1 has a weak correlation with azimuth frequency f_η. Therefore, a very small integrated quadratic modulation term of range frequency remains after the NLCS operation in the range direction. For the range edge point targets, the maximum phase mismatch of a single range pulse at a certain azimuth time can be given by

$$\mathsf{\Delta }{\varphi }_r=\pi \mathsf{\Delta }{k}_r{\left(\frac{{\tau }_p}{2}\right)}^2$$

(41)

where τ_p is the range pulse width and Δk_r is the offset of the range chirp rate related to the azimuth frequency f_η, which can be calculated by

$$\mathsf{\Delta }{k}_r=\left|k_{m1}-{k_{m1}|}_{f_{\eta }=f_{dc}}\right|$$

(42)

where |·| represents the absolute value of an element. The integrated quadratic frequency modulation term resulting from the above NLCS processing is very small, which can be confirmed by the range profiles of point targets 1 and 25 in Figure 9.

5.2. Ground Range Error Analysis

In order to investigate the impact of the proposed algorithm on slant range geometry projection, the following analysis is conducted. The shortest slant range of a focused target in the ground scene is

$$R_0^{\prime }=\frac{c\cdot D_{fref}}{2}\mathsf{\Delta }\tau +R_{ref}$$

(43)

where Δτ represents range time delay of the focused target relative to the central range cell. When the focused target is projected to the ground range plane, the corresponding range displacement relative to the scene center is

$$\mathsf{\Delta }{R}_g=R_0^{\prime }\cdot \sin \left({\beta }_0\right)-R_{ref}\cdot \sin \left(\beta \right)$$

(44)

where β and β₀ are look angles for central range cells and non-central range cells, respectively. Using the simulation parameters in Table 1, the relative range displacements for point targets 1 and 25 shown in Figure 9 can be obtained as −4999.46 m and 5000.32 m, respectively, and the corresponding ground range projection errors are 0.54 m and 0.32 m, respectively. Essentially, these projection errors are caused by the additional small linear phase modulation introduced by the NLCS operation in Section 3.2.

5.3. Computational Complexity Analysis

The computational complexities of the proposed MRDA and the referenced LRWC-based algorithm are also discussed in this section and are summarized in

Table 4. Comparison of computational complexities.

Algorithm	Computational Complexity
The LRWC-based algorithm	$C_1=2N_r{N}_a{\log }_2{N}_a+2N_a{N}_r{\log }_2{N}_r+5N_a{N}_r$
The proposed MRDA	$C_2=N_r{N}_a{\log }_2{N}_a+2N_a{N}_r{\log }_2{N}_r+\left(3+N_k\right)N_a{N}_r$

and Figure 11. Range and azimuth sampling points are denoted by N_r and N_a, respectively. The modified correlation window length in the range direction is represented by N_k, which is set to be 32 in the above simulation. Assuming the sampling points in range and azimuth directions are the same, i.e., N_r = N_a = N, and set to be 2¹⁴, the calculation complexities (C_i, i = 1, 2) in Table 4 are 1.64 × 10¹⁰ and 2.07 × 10¹⁰, respectively. The ratio of the computational complexities between the proposed MRDA and the LRWC-based algorithm can be calculated as

$$\frac{C_2}{C_1}=1.26$$

(45)

Equation (45) implies that the computational complexity of the proposed MRDA is increased by 26% compared to the LRWC-based algorithm using the simulation parameters listed in Table 1. For larger squint angles and better resolutions, the length of this correlation window will increase, which will further lead to a certain degree of decrease in the efficiency of the proposed algorithm.

5.4. Motion Error Compensation

The airborne SAR platform is easily affected by external factors such as airflow and weather in practical applications, which makes the actual motion of the SAR platform deviate from the ideal trajectory, resulting in unknown motion errors. Due to the impact of motion errors on imaging quality, SAR motion compensation also needs to be considered. Motion compensation based on echo data is the key step to achieve accurate imaging of the target area. There are two main approaches for motion compensation: parameter estimation and autofocus processing. In terms of parameter estimation, the research focuses on the estimation of Doppler centroid frequency and Doppler frequency modulation rate. The estimation methods of the former mainly include the energy equalization method, clutter locking method, radon transform method radon-transform-based method and its improved form [26,27,28,29]. The latter is represented by the map-drift method [30]. Phase Gradient Autofocus (PGA) is a classic autofocus method, which can realize high-order motion error compensation directly by estimating high-order phase errors [31]. The modified PGA based on the Weighted Maximum Likelihood (WML) can significantly improve convergence speed and promote the processing efficiency [32]. Accurate imaging of the target area can then be achieved by embedding the motion compensation into the imaging algorithms.

6. Conclusions

In this paper, a Modified Range Doppler Algorithm (MRDA) is proposed for a high-squint airborne SAR system, which can realize high-precision SAR echo imaging. The bulk compensation preprocessing is used to remove the significant RCM and severe range–azimuth cross-coupling of SAR echoes with large squint angles, which can simultaneously maintain the “azimuth-invariant” property. To equalize the range chirp rate, the range NLCS method is introduced in the imaging processing, after which the consistent range compression can be fulfilled by the range phase compensation operation in the range frequency domain. The modified correlation processing is performed in the RD domain to eliminate the residual range-variant RCM, the residual Doppler phase modulation and the range-variant cubic phase modulation introduced by the residual cross-coupling. The final focused high-squint airborne SAR image can be obtained via an AIFFT operation. Simulations are conducted in this paper to verify the effectiveness of the proposed MRDA. Evaluation results of the proposed MRDA show better performance than that of the LRWC-based algorithm for a high-squint airborne SAR system.