Focus Improvement of Airborne High-Squint Bistatic SAR Data Using Modiﬁed Azimuth NLCS Algorithm Based on Lagrange Inversion Theorem

: In this paper, a modiﬁed azimuth nonlinear chirp scaling (NLCS) algorithm is derived for high-squint bistatic synthetic aperture radar (BiSAR) imaging to solve its inherent difﬁcult issues, including the large range cell migration (RCM), azimuth-dependent Doppler parameters, and the sensibility of the higher order terms. First, using the Lagrange inversion theorem, an accurate spectrum suitable for processing airborne high-squint BiSAR data is introduced. Different from the spectrum that is based on the method of series reversion (MSR), it is allowed to derive the bistatic stationary phase point while retaining the double square root (DSR) of the slant range history. Based the spectrum, a linear RCM correction is used to remove the most of the linear RCM components and mitigate the range-azimuth coupling, and, then, bulk secondary range compression is implemented to compensate the residual RCM and cross-coupling terms. Following this, a modiﬁed azimuth NLCS operation is applied to eliminate the azimuth-dependence of Doppler parameters and equalize the azimuth frequency modulation for azimuth compression. The experimental results, with better focusing performance, prove the high accuracy and effectiveness of the proposed algorithm. Contributions: Conceptualization, C.L. H.Z.; methodology, software, vali-dation, C.L.; formal analysis, C.L.; investigation, C.L.; resources, Y.D.; data curation, C.L.; writing— preparation, C.L.; writing—review and editing, H.Z.; visualization, C.L.; supervision, Y.D.; project administration, Y.D.; funding acquisition, H.Z.


Introduction
Bistatic synthetic aperture radar (BiSAR) refers to a SAR system where the transmitter and receiver are separately mounted on different platforms, which, hence, increases the appreciable capability, flexibility, and reliability in designing BiSAR configuration [1][2][3][4]. Compared with the monostatic SAR system (MonoSAR), BiSAR system could offer many additional advantages, such as frequent monitoring, reduced vulnerability, single-pass interferometry, capacity of forward-or backward-looking SAR imaging, and saved costs using an existing illuminators of opportunity with several receive-only systems [2][3][4][5][6][7]. Generally, for airborne BiSAR, to obtain a good overlapping beam, the transmitter and receiver antennas both operate in squint mode. However, the squint BiSAR echo brings more complexities for imaging due to its acquisition geometry. To produce well-focused images, the particular problems for the high-squint BiSAR must be taken into consideration when designing the imaging algorithms.
Whereas time-domain algorithms can handle general bistatic cases [8], they are very inefficient; therefore, frequency-domain methods are preferred. There are two main challenges in developping the frequency-domain algorithm for processing high-squint BiSAR data, the first one is an accurate analytical solution of the 2-D spectrum for BiSAR, which is the basis of the high-precision bistatic imaging algorithm. Because of the range history of a bistatic target is the sum of two different hyperbolic range equations, which leads to The rest of this paper is organized, as follows. The echo model of the general BiSAR is introduced in Section 2, and a more accurate spectrum suitable for processing high-squint BiSAR data is introduced in this section. In Section 3, the derivation of the modified azimuth NLCS algorithm is described in detail. The range residual RCM is evaluated, and the azimuth-dependent Doppler parameters are analyzed in this section. The point target simulations of high-squint BiSAR imaging are employed to validate the proposed algorithm in Section 4. Finally, the conclusion is drawn in Section 5. Figure 1 presents the general bistatic geometry, in which the motion trajectories of transmitter and receiver are not parallel. The SAR sensors fly along a given direction with constant velocities V T and V R at constant altitudes H T and H R . t c is defined as the composite beam center cross time of the reference target P, and θ ST and θ SR are the squint angles for target P at t = t c . R 0T and R 0R are the slant ranges of the transmitter and receiver to the target at t c , respectively. Suppose that the transimitted signal is a linear FM pulse, the received echo after demodulation is

BiSAR Signal Model and LIT Spectrum
where w r and w az are the range and azimuth envelopes, respectively. τ is fast time and t is the slow time. f c is the carrier frequency, k r represents the range frequency modulation rate, and c is the speed of light. From Figure 1, we can obtain the expression of the instantaneous bistatic slant range R(t), which is given by (2) Substituting (2) into (1) and preforming the 2-D Fourier transform to (1) yields where W r represents the range frequency envelope, f τ is the range frequency, and f t is the azimuth frequency. The principle of stationary phase is usually used to solve the integral term in the above formula. However, because of the existence of the DSR term in the bistatic slant range history R(t), it is difficult to directly obtain the analytical expression of the bistatic stationary phase pointt. Many approximate methods [2,[9][10][11]14] have been proposed to solve this problem. Among them, the MSR spectrum [14] first performs Taylor expansion of R(t) in (2), which is, , and, then, after removing the linear term, MSR is used to find an approximate solution oft. The accuracy of this spectrum is controlled by the two expansions of the slant range and the series inversion. Different from MSR, LIT, in [11], is allowed to derive the bistatic stationary phase point while retaining the DSR term of R(t). Combining LIT and the principle of stationary phase, the spectrum in (3) can be written as where W az represents the azimuth frequency envelope andt can be expressed as We can see thatt is derived directly by the Lagrange inversion theorem without intermediate steps, such as the linear term removal and the shift in time domain in the case where the MSR is used [11]. The accuracy of LIT spectrum is only controlled by the number of analytical solution of the stationary phase retained in (6). Additionally, according to reference [11] and [21], under the same conditions, the bistatic spectrum accuracy that is based on LIT is higher.

Modified LIT Spectrum for High-Squint BiSAR
Although the idea of LIT spectrum solves the problem of deriving the bistatic stationary phase point while retaining the DSR, and the accuracy of the spectrum is improved, it also limits the derivation of the high-squint BiSAR imaging algorithm that is based on the LIT spectrum. In [21], the RDA based on LIT spectrum is derived for processing the azimuth-invariant synthetic aperture sonar data. The RDA directly deals with the original data, but, for high-squint BiSAR data, the target's 2-D spectrum is skew and causes parts of the spectrum to cross into the adjacent pulse-repetition frequency (PRF) band. Therefore, the 2-D skew spectrum usually requires high PRF or extra computational load in the imaging processing due to the increased samples in azimuth [22,23,32] with a high squint and big bandwidth. Additionally, since the range dependence of the secondary range compression and the azimuth dependence of the Doppler parameters are not considered, the RDA is only suitable for dealing with the BiSAR data with azimuth-invariant or low squint case. To solve the above problems, a new accurate spectrum that is based on LIT suitable for processing high-squint BiSAR data is introduced next.
First of all, the bistatic slant range R(t) can be rewritten as The above bistatic slant range approximation incorporates the quadratic range curvature (ie, quadratic term of slow time) into the square root term, which is similar to the form of RCM in broadside BiSAR mode. (7) has higher accuracy than the traditional third-order approximation [33,34], and Figure 2 shows the relationship between the slant range errors that are induced by the different approximations and the azimuth exposure time, and the imaging geometric parameters are shown in Table 1. Each curve represents the error of a target point under a different range gate, and the range interval between the target points is 400 m. In addition, (7) also separates the linear RCM term to prepare for the subsequent LRCMC operation, and the DSR term in (7) can more intuitively show the link between the high-squint BiSAR spectrum after the LRCMC and the broadside BiSAR spectrum. Substituting (7) into (1) and applying the range Fourier transform to convert the signal to the range frequency domain, we can get the expression of the signal at this time, as where the phase term is From Formula (8), it can be seen that the trajectory of the bistatic target contains both linear RCM and non-linear RCM components. The linear RCM component arises when the composite antenna beam is squinted. Additionally, the range-azimuth coupling mainly comes from the linear RCM component, particularly in the high-squint BiSAR mode [22][23][24]27]. Therefore, in order to alleviate range-azimuth coupling, many highsquint MonoSAR imaging algorithms apply the LRCMC operation in the azimuth time range frequency domain to remove linear RCM component [23,24,[27][28][29][30]. This idea is also applicable for high-squint BiSAR data processing. Additionally, the phase matching function of LRCMC operation for BiSAR is Multiplying the above formula with (8) and observing the residual phase term in (8), we can find that the function of LRCMC is equal to transforming the squinted BiSAR signal to the broadside BiSAR signal, with a change of the moving velocity from V i to V i cos(θ i ) (i = R or T), which is analogous to the role of LRCMC in squinted MonoSAR [23,24,27]. Figure 3 shows the 2-D spectrum of high-squint BiSAR raw data before and after LRCMC operation; the simulation parameters are shown in Table 1, where the squint angles of transmitter and receiver are θ ST = 50 • and θ SR = 45 • , respectively, and the PRF is set to 250, which is approximately twice the Doppler bandwidth. From Figure 3a, due to the the skew of the 2-D spectrum, part of the spectrum is folded into the adjacent PRF band, causing azimuth aliasing. After the LRCMC operation, the spectrum skew phenomenon is corrected, the azimuth aliasing is eliminated, and the orthogonality between the range and the azimuth is significantly increased, which means that the range-azimuth coupling of the 2-D spectrum is significantly alleviated. In [23], the MonoSAR spectrum after LRCMC is called the squint-minimized spectrum, which is also appropriate for the high-squint BiSAR spectrum. However, the LRCMC factor in (9) varies with range, a practical way to deal with this problem is to perform LRCMC within an invariance region to keep the variations of squint angles small, and the size analysis of invariance region can be found in [16,20] (see Appendix C in [16] and Section IV in [20]).  When applying azimuth Fourier transform to (8), the modified spectrum can be derived by means of LIT, which is, the phase term is wheret M is the new bistatic stationary phase, which be solved in the same way as (6), and its analytical expression is It can be seen from the above derivation process that the modified LIT spectrum in (10) can be easily combined with the LRCMC operation, which makes it a reality to derive the high-squint BiSAR imaging algorithm that is based on the idea of LIT spectrum. Next, based on this spectrum, a new bistatic imaging algorithm will be introduced.

A Modified Bistatic Azimuth NLCS Algorithm for the High-Squint BiSAR Imaging
The proposed algorithm in this section mainly includes five steps. First, a LRCMC operation in the azimuth time range frequency domain is to remove most of the linear RCM and mitigate the range-azimuth coupling. Second, a BSRC operation completes the compensation of residual RCM, secondary range compression, and high-order coupling terms in 2-D spectrum. Third, a fourth-order azimuth filtering operation is used to increase the degree of freedom in deriving the azimuth NLCS coefficient expression. Fourth, an azimuth NLCS operation in range Doppler domain is to equalize the azimuth FM rate and remove the azimuth dependence of the quadratic FM rate and cubic phase. Finally, a residual azimuth phase compression operation completes the azimuth focusing. In order to clearly present each step, the flowchart of the proposed algorithm is summarized in Figure 4, where the first two steps belong to range focusing processing, and the last three steps belong to azimuth focusing processing. In what follows, the necessary derivations and explanations of each step are detailed.

Range Focusing via the Bulk Secondary Range Compression Operation
After the LRCMC operation in (9), the RCM mainly includes quadratic and higher order terms, although the higher order terms are usually very small when compared to the quadratic term. Uncorrected range curvature leads to impulse response degradation in both the range and azimuth directions [16,22]. The range focusing of the signal that is processed by LRCMC mainly faces two problems. The first is residual RCM correction.
Combining (8)-(11), the 2-D spectrum after LRCMC operation can be obtained as where the first exponential term is the range FM term, which is used for range compression by multiplying its complex conjugate. The second exponential term contains target azimuth position information, it can be seen that, even though the linear coupling term has been removed, an azimuth-position-dependent range shift is brought about, which further leads to the azimuth spatial variance and brings great challenges for azimuth focusing.The range-azimuth coupling terms are the third and fourth exponential terms. Expanding the phase term in (12) into a power series of the range frequency f τ , the result can be given by where In (13), φ 0 ( f t ) contains all of the information of the azimuth modualation. φ 1 ( f t ) is linearly dependent on the range frequency f τ and it represents the residual linear RCM. φ 2 ( f t ), like the first exponential term in (12), is the quadratic term of f τ , and it can be used to reduce the mainlobe broadening of the range impulse response after range compression, so it is also called the secondary range compression (SRC) term [35]. Additionally, we can see that the terms are range dependent, next it is necessary to analyze the range-dependent characteristics of the residual RCM and the coupling phase.
Consider a 1-m resolution bistatic X-band system with a transmitter range of 14.1 km and a receiver range of 15.2 km from the scene center target T0, as shown in Figure 5a. The transmitter/receiver separation or baseline is 2 km, the squint angles are θ SR = 30 • and θ ST = 45 • . The distance interval between the targets T0, T1, and T2 is 2.5 km. Target T2 has the same minimum slant range as target T4 and the same beam center crossing time as targets T0 and T1, targets T0, T3, and T4 have a similar linear RCM slope and lie in the same range cell after LRCMC, as shown in Figure 5b. Thus, the LRCMC operation results in another problem that is an azimuth-dependent displacement in range envelope. In other words, the azimuth targets thhat are located in the same gate bin will have different RCM trajectories in the RD domain, as shown in Figure 5b, which shakes the fundamental assumption of the traditional CS and RD algorithms and makes the algorithms no longer applicable. We can see that the target points T0, T3, T4 with different RCM trajectories are at the same range gate after the LRCMC operation. At this time, the traditional RD and CS algorithm will no longer be applicable. Figure 6 shows the evaluation results of the residual RCM at each point in the scene. From Figure 6a, it can be seen that, after the LRCMC operation, the residual RCMs are greater than three range resolution units, so the residual RCM correction is indispensable. The residual RCM change between the midswath and edge targets (measured at the end of the azimuth aperture) is less than 20% of a resolution cell, which means that the range dependence of the residual RCMs of the targets in the scene is negligible. For second problem, due to the residual RCMs of targets T0, T3, and T4 coincide quite closely in azimuth frequency domain, it can be removed with the same RCM correction shift. For more precise processing, the ideas in [36] are recommended. Thus, a bulk compensation function is adopted to complete residual RCM correction. It is essentially the complex conjugate of the BiSAR two-dimensional spectrum in (12) at a reference range (R 0Rre f , R 0Tre f ) (that is usually set to the center of the scene), and it can be expressed as The above compensation function is similar to the bulk range processing transfer function in [33], which can perform range compression, residual RCM correction, secondary range compression, and the compensation of all higher-order phase terms for all points located at the reference range. Additionally, it uses the double root form of the 2-D spectrum in (12) instead of its series expansion in (13), which can effectively reduce the phase error that is introduced by the root term expansion, especially for the bistatic imaging with a large fractional bandwidth [37]. Following the reference [24], this compensation function is called a bulk secondary range compression (BSRC) operation.
The BSRC operation can avoid a time-consuming interpolation operation without affecting the quality of range focusing, which is more efficient than the bistatic RDA in [7,15,21] and NLCS algorithm in [16]. While the BSRC that corrects the phase at the scene center is often sufficient for the whole scene, for wider range swaths, it may be necessary to segment the scene into range invariance regions whose width is selected by the QPE allowed [22,24].

Azimuth Focusing via the Modified Bistatic Azimuth NLCS Operation
The trajectories become parallel to the azimuth frequency axis after the LRCMC and BSRC operation, and then only the azimuth compression needs to be completed. However, the azimuth signal in a given range cell now consists of a group of point targets that originated from different ranges and, comsequently, have different azimuth Doppler parameters. For instance, we consider targets T0, T3, and T4, which lie in the same range cell after LRCMC, as shown in Figure 5. The azimuth phase curves of the three targets have different frequency modulation (FM) rates, because this parameter varies with the original range position of the target. This is the problem that needs to be solved in this subsection.

Azimuth-Dependent Characteristic Evaluation
Multiplying (14) with (12) and transforming the result into the range Doppler domain by the range inverse Fourier transform, the range focused BiSAR signal can be expressed as where The first two exponential terms in (15) are constant and linear shift terms, respectively, which contain the position information of the target. The azimuth FM terms are the third and fourth exponential terms, which are essential for azimuth focusing in high-squint case. In traditional imaging algorithms, azimuth compression is accomplished by applying a matched filter in the range Doppler domain, which is the complex conjugate of the azimuth FM terms. However, because LRCMC operation causes these modulation terms to be azimuth-dependent, the direct matched filtering method will no longer be available. Therefore, next we need to evaluate the azimuth-dependent characteristics.
First of all, the azimuth quadratic FM rate k a2 can be expressed as where k a20re f is the azimuth quadratic FM rate at reference position (R 0Rre f , R 0Tre f ), and k a21re f and k a22re f are the cofficients of the spatial variance components of k a2 .
In [16], the first-order spatial variance term k a21re f is corrected. However, in highsquint case, the influence of the second-order spatial variance term k a22re f cannot be ignored. Figure 7a shows the variation range of the azimuth quadratic FM rate with the azimuth position, and Figure 7b shows the relationship between the quadratic phase errors (QPEs) that are induced by the different order approximations of the azimuth quadratic FM rate and azimuth position, the simulation parameters are shown in Table 1. A phase error of π/4 is usually used as an upper limit to obtain a good focusing quality [22]. When the azimuth imaging swath exceeds 1.84 km, the QPE is larger than π/4, and the azimuth focusing result of the bistatic NLCS algorithm in [16] will deteriorate. The quadratic FM rate directly affects the main lobe broadening of the focusing result. An accurate azimuth FM rate for azimuth compression can significantly reduce the azimuth main lobe broadening, thereby improving the azimuth resolution. Next, the azimuth cubic phase k a3 can be expressed as where Under the same simulation parameters, as compared with quadratic FM rate, the azimuth-dependent cubic phase is much smaller. Figure 8 shows the relationship between the azimuth-dependent cubic phase and azimuth position in different squint cases. When the azimuth imaging swath is larger than 3.6 km at the squint angles are θ SR = 50 • and θ ST = 55 • , the phase error that is caused by cubic phase will be larger than π/4, the spatial variance of cubic phase term becomes significant. Theoretically, the cubic phase mismatch will cause the sidelobe asymmetry of the focusing result, which will lead to the deterioration of the peak sidelobe ratio (PSLR) and integral sidelobe ratio (ISLR). Therefore, for high-squint biSAR imaging, the azimuth-dependent cubic phase term correction, which is ignored in [16], needs to be considered.

Derivation of the Modified Bistatic Azimuth NLCS Coeffcients
According to the analysis that is presented in the previous sections, the purpose of azimuth processing at this time is clear, that is, to eliminate the azimuth dependence related to t c in k a2 and k a3 . It can be solved by using the idea of azimuth NLCS in [16,24], and the schematic diagram is shown in Figure 9. After NLCS operation, the different phase curves of different targets in the same range gate are unified. However, the algorithm faces two problems when dealing with high-squint BiSAR data. One is that the first-order approximation of the azimuth quadratic FM rate is used when deriving the perturbation function, which will cause the mainlobe of the azimuth focusing result to broaden and the resolution to decrease, the other is ignoring the azimuth dependence of the cubic phase term shown in (17), which will lead to insufficient azimuth focus depth and asymmetric sidelobes. First of all, it needs to be explained that, due to the LRCMC operation, the spectrum of the reference target has been moved to the center of the azimuth frequency, which is, the interval of Dopper frequency changes from to f dc − B a/ 2 ≤ f t ≤ f dc +B a/ 2 to −B a/ 2 ≤ f t ≤ B a/ 2, where f dc is Doppler center and B a is the Doppler bandwidth.
Combining (15)- (17), the range focused BiSAR signal can be rewritten as Before performing the azimuth NLCS operation, a fourth-order azimuth frequency filter should be first applied, which is analogous to the nonlinear FM filtering method implemented in the range frequency domain in [38] to eliminate the range dependence of the secondary range compression. Additionally, the main purpose of the fourth-order filter is to derive the expressions for the chirp scaling phase function coeffcients. The fourth-order filter factor is given by The azimuth filtered BiSAR signal can be represented by (20) Subsequently, the signal is transformed to the 2-D time domain and multiplied by a new azimuth NLCS factor. The new azimuth NLCS factor is as follows: When compared with the disturbance function in [16], which is the cubic term of the azimuth time, the new azimuth NLCS factor belongs to an extended form.
Converting the scaled signal to the azimuth frequency domain, and then expanding it into a Taylor series of t c and f t and keeping up to the third-order coupling term yields To remove the azimuth dependence of the quadratic FM rate and cubic phase, the coefficient of the first-order coupling term is set to −π/a, where a = 0.5 is an introduced constant factor, and all of the coefficients of higher order coupling terms are set to zero.
Thus, we get an equation system with five equations and five unknowns. The solution to the equation system is , Additionally, after azimuth NLCS operation, the residual azimuth phase term, which is azimuth-independence, is as follows: Thus, after the fourth-order azimuth filter and modified azimuth NLCS operation, the expression of the BiSAR signal can be given by Inpecting (25), it is chear that the azimuth-dependent characteristic of the azimuth modulation terms has been completely removed. The first exponential term represents the azimuth position of target, and it is dependent on the parameter a. Here, a is chosen to be 0.5, so that the targets are in their real positions after azimuth compression. The second exponential term is the residual azimuth phase term, which needs to be compensated to complete the final azimuth compression, so the azimuth compression factor is given by Subsequently, transform the signal after azimuth compression to the azimuth time domain, the final focused SAR image is obtained as From (27), we can find that the target is focused at position τ = (R 0 + k 0 t c ) c, t = t c/ 2α , which is offset from the correct position. This offset can be easily corrected through geometry correction. Related methods can be found [16,24], and they will not be repeated here.

Simulation Results by the Modified Bistatic ANLCS Algorithm
The experiments with simulated data using airborne BiSAR parameters that are shown in Table 1 are carried out to prove the accuracy and imaging performance of the proposed algorithm on processing high-squint BiSAR data.
Following the experiment that was done in [16], the simulation uses an array of 25 targets, which are laied out on a 1.6-km 2 grid in ground range/azimuth, as shown in Figure 10. The point target layout corresponds to the range invariance region that limits the resolution broadening to 3% [16,20]. A rectangular window is used in both the range and azimuth processing, instead of an antenna pattern, i.e. no weighting is used, which enables more quantitative evaluation of imaging results. The oversampling ratio is 1.2 in range; thus, the theoretical range resolution is 0.75 m, which is based on the definition of resolution in [39,40]. According to the PRF and velocity of the senor, we can roughly calculate the azimuth oversampling rate is 3.63 and the theoretical azimuth resolution is 1 m.   Figure 11 provides the results from the main steps of our proposed algorithm. Figure  11a gives the echo data of the 25 point targets in 2-D time domain. Figure 11b,c, respectively, show the RCM trajectories of each target in 2-D time domain after range compression and LRCMC operation. It is observed that, after the LRCMC operation, most of the linear RCM is removed, but this also causes the target's range position to shift. Figure 11d,e present the RCM trajectories of each target before and after the BSRC operation in the range Doppler domain, respectively. From the two figures, the residual RCM is corrected by BSRC operation, and the range energy of each target is gathered into one range gate. After the LRCMC operation, the Doppler parameters of the targets 5, 13, and 21 at the adjacent range gates are different, as shown in Figure 11f. For the convenience of comparison, the phase curves in the figure are obtained after phase unwrapping and removing the constant term. Figure 11g provides the azimuth phase curves of the three targets after the proposed modified azimuth NLCS operation. From the figure, the different azimuth phase curves of different targets at same or adjacent range gates are unified. After the above-mentioned steps are processed, the energy of each target in the 2-D time domain gathered into one range bin in both the range and azimuth directions, which can be seen by comparing Figure 11c,h. Additionally, Figure 12 shows the final resulting image.  In the imaging, the chirp scaling factor a is set to 0.55. Moreover, the geometric correction is unused, which is helpful for computing the measured parameters of the targets. Observing the results, we can find that all of the targets are well focused. Target 13 is the center of the scene (which is also the reference point target), target 25 is the furthest target from target 13, and, therefore, suffers the largest azimuth and range phase degradation. They are the most representative points to evaluate the imaging performance of the algorithm [16,22]. Figure 13 and Figure 14 show the impulse responses of these two representative target points, respectively, with the measured values of the resolution, the peak sidelobe ratio (PSLR), and the integrated sidelobe ratio (ISLR) annotated. For precise comparison, the figures also give the imaging results (red line) that are focused by the backprojection (BP) algorithm in [8], which is a time-domain algorithm with high accuracy.
From Figure 13, it can be seen that the reference target 13 is focused perfectly, it has no broadening in both range and azimuth, and its measurement indexes (blue font) are very close to that obtained by the BP algorithm (red font). According to the contour plot in Figure 14a, the focusing result of target 25 as compared to target 13 is not perfect due to the range and azimuth phase degradation. However, for the edge target, there is negligible broadening in range, whereas azimuth broadening is less than 3%. The range PSLR deviates from the theoretical value of −13.3 dB by less than 0.06 dB, the azimuth one is less than 0.12 dB, and the range and azimuth ISLR is within 1dB from the expected value of −10.0 dB. Moreover, under the given simulation parameters, the proposed algorithm can obtain similar focusing capabilities with the BP algorithm, and the depth of focus is slightly inferior to the time domain algorithm, which prove the excellent imaging performance of the proposed frequency domain algorithm.

Comparison Results and Discussion
The existing algorithms, i.e., the Zhang's indirect RDA based on the LIT spectrum in [21], and Wong's NLCS algorithm based on MSR spectrum in [16], are employed to process the high-squint BiSAR data in same simulation parameters.
Zhang's RDA is derived for processing synthetic aperture sonar data, which is azimuth-invariant, which is, the velocities of the transmitter and the receiver are equal, and the baseline between them is constant. This algorithm cannot be used to deal with BiSAR data with general configuration, so we first need to modify Zhang's RD algorithm. Zhang's RDA mainly contains three transfer functions, which are derived from Bamler's summary of existing imaging algorithms [33]. By combining the LIT spectrum in (6) and Zhang's considerations for deriving RDA, we give the three modifed transfer functions that can process BiSAR data with general configuration.
The first is the azimuth transfer function, which can be expressed as wheret It is obtained after the range frequency of the 2-D spectrum is set to zero, which is, f τ = 0. It is used for the azimuth compression of the algorithm.
The second function is the bulk range processing transfer function (BRTF), and its expression is where conj stands for phase conjugation operation. BRTF is the residual term after removing the azimuth modulation and the range position information of the 2-D spectrum at a reference range (R 0Rre f , R 0Tre f ). It is similar to the aforementioned BSRC operation, and its main purpose is to remove the secondary range compression term and azimuth-range coupling terms at the reference range. However, the difference from the BSRC operation is that, due to the lack of LRCMC operation, the RCM in the BiSAR data processed by BRTF is very large, and the azimuth-range coupling is serious, so further processing of the RCM in the data is required. The last is differential RCMC, which can be expressed as Its purpose is to remove the residual RCM that is related to the range in the RD domain to complete range compression. This process is realized by utilizing interpolation operation. Figures 15 and 16 show the subimages of target 14 and target 15, which are extracted from the BiSAR images obtained by the three imaging algorithms, respectively. Figures 15a and 16a are obtained using the aforementioned modified Zhang's RDA. It can be seen that the algorithm directly processes the original high-squint BiSAR data without "squint minimization", so the focus result is also skew. The skew is caused by the variation of the azimuth FM rate with range in the matched fliter, which interacts with the nonzero Doppler centroid to change the registration in the azimuth direction [22]. Because the algorithm ignores the range dependence of the secondary range compression term and the high-order range-azimuth coupling term, the residual RCM of the targets at the non-reference positions are very large, and the azimuth sidelobes of the targets cross several range bins, which results in a mismatch in the azimuth FM rate. Thus, both the azimuth and the range are defocused. And, as the distance in range direction between the target and the reference position increases, the target-focusing quality degrades quikly. Figures 15b and 16b are obtained using Wong's NLCS algorithm. Because the algorithm uses LRCMC operation to alleviate the range-azimuth coupling, the range dependence of the residual RCM is weakened, so the focus quality of the target is significantly improved. However, the main lobe of the target in azimuth direction is broadened, due to the firstorder approximation in (16) for azimuth quadratic FM rate, and azimuth sidelobes are unsymmetrical due to ignore the influence of the azimuth dependence of the azimuth cubic phase term in (17). The imaging results obtained by the proposed algorithm are shown in Figures 15c and 16c, and all of the targets in range and azimuth direction are well focused. A worse imaging quality is obtained by the RD algorithm and the NLCS algorithm than the proposed algorithm because of the lack of residual RCM compensation and the increased mismatch of azimuthdependent azimuth FM rate. To further quantify the performance of the proposed algorithm, the resolution, PSLR, and ISLR are also used as the performace measure. The corresponding imaging performances of the targets 14, 15, 18, and 23 are calculated, and the results are listed in Table 2. Inspecting the Table 2, it is found that the 2-D measured parameters resolution, PSLR, and ISLR obtained by the proposed algorithm agree with the theoretical values. However, the measured parameters that were obtained by the Zhang's RDA and Wong's NLCS algorithm are not good, especially in the azimuth direction.
According to the algorithm flowcharts, we now compare the computational complexity of these three algorithms. Generally speaking, the 5Nlog 2 (N) floating-point operations (FLOPs) are needed to compute N-point fast Fourier transform or inverse fast Fourier transform and 6N FLOPs are needed for one-time complex multiplication [27]. Suppose that the azimuth sample number is N a , the range sample number is N r , and the interpolator kernel is N k . Therefore, the computational loads of these three algorithms, in terms of FLOPs, are, respectively C Zhang = 10N a N r log 2 N r + 10N a N r log 2 N r + 12N a N r + 2(2N k − 1)N a N r , C Wong = 20N a N r log 2 N r + 20N a N r log 2 N r + 30N a N r + 2(2N k − 1)N a N r , C proposed = 10N a N r log 2 N r + 20N a N r log 2 N r + 30N a N r .
It can be seen that the computational load of the proposed algorithm is better than the Wong's NLCS algorithm, but slightly heavier than the Zhang's RDA. However, since there is no time-consuming interpolation operation, its processing efficiency is higher.  In addition, in order to more intuitively demonstrate the imaging performance of the proposed algorithm, a scene with distributed targets is simulated, in which there are 2000 × 2000 point targets, the backscatter coefficients are replaced by the grayscale values of a existing SAR image, and the interval between the target points is one resolution. Figure  19 shows the imaging results that were processed by the proposed algorithm and the BP algorithm. It can be seen that the scattering characteristics of the SAR images are well displayed, and their focusing results are nearly consistent.

Conclusions
In this paper, a modified bistatic azimuth NLCS algorithm for the high-squint BiSAR data processing has been discussed. The algorithm is based on an accurate modified LIT spectrum that is suitable for processing high-squint BiSAR data. In the algorithm, a LRCMC operation is first performed to remove the most of the linear RCM components and mitigate the range-azimuth coupling. Subsequently, a BSRC operation is applied for compensating the residual RCM and higher order cross-coupling terms to complete the range focusing processing. Finally, by adopting higher order approximation to the azimuth FM rate and cubic phase term, a modified azimuth NLCS operation is used for eliminating the azimuth-dependence of Doppler parameters and equalizing the azimuth FM rate to accomplish azimuth compression, and the more accurate focusing results prove the effectiveness of the proposed algorithm. Data Availability Statement: All data generated or analyzed during this study are included in this article.