Focus Improvement of Spaceborne-Missile Bistatic SAR Data Using the Modiﬁed NLCS Algorithm Based on the Method of Series Reversion

: The speed and direction of a missile shifts sharply in the dive phase, making the azimuth frequency modulation (FM) rate change with the azimuthal position, leading to azimuth ambiguities and image distortion. To solve this problem, a modiﬁed nonlinear chirp scaling (NLCS) algorithm was adopted to compensate for the azimuth FM rate. First, the geometric conﬁguration and echo signal model of the spaceborne missile bistatic synthetic aperture radar (SAR) were built, and then the Doppler frequency correction was performed, and the 2-D spectrum of the signal was derived by the method of series reversion. Next, range migration correction and range compression were ﬁnished in the 2-D frequency domain. Following this, a modiﬁed NLCS algorithm was proposed to solve the space variance of Doppler phase problem. After compensating for the azimuth FM rate, the azimuth compression focusing was completed and the imaging result was obtained. Finally, by comparing the calculation amount, imaging effect, and performance index with the traditional NLCS algorithm, it can be concluded that the algorithm reduced the calculation amount by 1.0128 × 10 8 ﬂoating point operations per second (FLOPs) compared with the traditional NLCS algorithm, and the azimuth focusing effect of the edge point was greatly improved. Its resolution, peak sidelobe ratio (PSLR), and integrated sidelobe ratio (ISLR) were improved by 0.87 m, 3.32 dB, and 1.79 dB, respectively, which proved the effectiveness and feasibility of this method.


Introduction
The synthetic aperture radar (SAR) can penetrate through all weather conditions and can operate all day long [1,2]. It has been widely used in combat and geological surveys [3,4]. Compared with the traditional monostatic SAR, the new bistatic SAR has the advantages of strong concealment, rich information acquisition, and forward-looking imaging [5]. As the transmitting source, the satellite operates at a certain orbital altitude, with strong anti-strike capability and wide beam coverage [6]. As the receiver, the missile has the advantage of high resolution [7]. The new system of spaceborne missile bistatic SAR can combine the advantages of spaceborne SAR and missile-borne SAR. It can combine the large-scale coverage and high security of spaceborne SAR with the high resolution and high mobility of missile-borne SAR, effectively obtaining scene targets including background, improving the ability to strike targets, and realizing the integration of observation and attack [8][9][10].
Although there are many advantages in the special geometry configuration, it also brings about two challenges in applying the algorithm to processing the echo signal data. First, it is difficult to obtain an accurate expression of the 2-D spectrum for spaceborne missile bistatic SAR, which is the basis of the imaging algorithm. In a bistatic system, the range history is the sum of two hyperbolic range equations, and therefore there is a double square root term in the range history [11], which makes it difficult to solve the 2-D spectrum analytically by directly applying the stationary phase principle [12][13][14][15]. The derivations about the spectrum of the bistatic SAR were completed in [16]. However, these derivations only apply to some special geometric configurations (azimuth-invariant or parallel) that do not directly obtain the analytical solution of the bistatic stationary phase leading to a significant decrease in accuracy. For the flexible spaceborne missile bistatic system, these approaches in [16] cannot be applied to handle it. To solve the accuracy issue, the Legendre expansion [17] and orthogonal expansion [18] methods have been preliminarily studied, but solving a higher-order 2-D spectrum is challenging. The 2-D spectrum was derived in [19] using the Lagrange inversion theorem (LIT). The bistatic stationary phase can be derived while maintaining the double square root of the slant range history. Furthermore, its accuracy is solely determined by the number of analytical solutions retained in the stationary phase. In this paper, we derived the 2-D spectrum of spaceborne missile bistatic SAR by using the method of series reversion [20]. The idea is to apply Taylor expansion to the slant range history, then determine the analytical solution of the stationary phase point by series reversion after removing the linear range cell migration and the linear phase term. The precision of this spectrum is controlled by the two slant range expansions as well as the series reversion.
The second challenge comes from the space variance of the Doppler phase. The modified Omega-K, chirp scaling algorithm, and range-Doppler algorithms reported in [16,[21][22][23] are the forward-looking imaging algorithms that are suited for airborne and spaceborne bistatic SAR, but they are not suitable for missile-borne bistatic SAR. Back-projection algorithms are capable of handling general bistatic configurations with topographic dependences and flexible moving platforms [24][25][26][27], but the efficiency of the algorithm is low. A nonlinear chirp scaling algorithm that is more effective for handling data with the bistatic geometry where the transmitter or receiver is stationary has been presented [28,29]. In these papers, the azimuth focusing depth is limited as only the symmetrical azimuth FM rates and the parabolic are considered. Moreover, [30] considers the nonlinear component of azimuth FM rate, but only monostatic squinted missile-borne SAR data can be processed by this algorithm. In [31], a frequency-domain algorithm is proposed for curvilinear flight SAR, which compensates the spatial variation of the Doppler frequency modulation term. However, the residual azimuth spatial variation leads to an affect imaging focus in the azimuth direction for large scenes.
On the basis of the previous research, a modified NLCS algorithm was used to process the data of spaceborne missile bistatic SAR. The algorithm can be divided into the following steps. First, the echo signal and the geometric configuration of spaceborne missile bistatic SAR were established. Then, the Doppler frequency correction operation was conducted in range frequency and the azimuth time domain. Second, the 2-D spectrum of the spaceborne missile bistatic SAR echo signal was successfully derived by the method of series reversion. Following this, the range migration correction and range compression were completed in the 2-D frequency domain, and the imaging position in the range was determined. Finally, a modified NLCS algorithm was applied to compensate for the azimuth FM rate, eliminating the distortion of the image in azimuth, and a great imaging result was presented in the simulation.
The rest of the paper is organized as follows. Section 2 establishes the geometry of the spaceborne missile bistatic SAR and builds the echo signal model. As outlined in Section 3, the Doppler frequency correction was performed in the range frequency and azimuth time domain. What is more, the accurate 2-D spectrum of the echo signal was derived in detail by using the method of series reversion. Moreover, Section 3 shows the performance of the range migration correction and range compression in the 2-D frequency domain, determining the imaging position in range of the point target and describing the derivation of the modified NLCS algorithm in detail. Section 4 outlines the performance of the simulated experiments and presents the simulation results, comparisons between different algorithms, and imaging performance analysis. Finally, Section 5 concludes this paper.

Spaceborne Missile Bistatic SAR Geometry and Signal Model
The spaceborne missile bistatic SAR configuration with a moving satellite transmitter and a moving missile receiver is shown in Figure 1. The satellite operates in the Earth Centered Inertial coordinates, and its coordinate is determined by orbital elements. It is in a different coordinate system from the missile, and its coordinate needs to be converted into the local coordinate system [32]. The coordinate of the satellite in Figure 1 was converted to the local coordinate system. The initial coordinate of the satellite after conversion was in the same coordinate system as the missile, and corresponding calculations can be performed directly. The initial coordinates of the satellite is (x T , y T , z T ), the velocity of the satellite along the y-axis (range) direction is v ty , and the velocity of the satellite along the x-axis (azimuth) direction is v tx . The initial coordinate of the missile is (0, y R , z R ), and v y and v z are the velocity of the missile along the y-axis (range) direction and z-axis (height) direction, respectively. The acceleration of the missile along the y-axis (range) direction and z-axis (height) direction are a y and a z , respectively. For the convenience of calculation, the coordinate origin O was set as the target point at the center of the scene. P(x, y, z) represents an arbitrary point target in the imaging area. R T0 and R R0 represent the initial distances of the satellite and missile to the origin of the coordinates at the original time, respectively. R T (η) and R R (η) represent the instantaneous slant ranges of the satellite and missile to the point target P at slow time η, respectively. The transmitter was a geostationary earth orbit satellite that can be considered to be flying at a constant velocity in synthetic aperture time due to the relatively small amount of velocity change. In the dive-down section, the missile moved along the curve with a constant acceleration. η and τ represent the azimuth time (slow time) and range time (fast time), respectively. Assuming the initial moment of the system, the azimuth time η = 0, and the amplitude of the signal echo at each target point was set as the ideal case, not varying with the bistatic slant range R(η). Ignoring the problem of system synchronization, it was assumed that the operation of the satellite and the missile was synchronized at all times, and the electromagnetic scattering coefficient of the target remained unchanged [33]. Suppose that the transmitted signal was a linear frequency modulated (LEM) signal, and the received echo after demodulation can be written as where w a and w r represent the azimuth and range envelopes, respectively. f c represents the carrier frequency; K r is the range frequency modulation rate; c represents the speed of light; and R(η) represents the bistatic range history, which is the sum of R T (η) and R R (η). R T (η) and R R (η) can be described by and in order to further analyze the slant ranges, R T (η) and R R (η) can be expanded to their Taylor's series at η = 0, and then we have Remote Sens. 2022, 14, 5770  Different from airborne vehicles that fly along a certain direction with constant velocities, the missile is in the dive phase, with two-dimensional acceleration, and its velocity magnitude and direction change quickly, with the first-order to fourth-order term range cell migration being serious. Therefore, it is important to take the efficiency approximation of each order term into account. Some simulations are presented using the parameters listed in Table 1, taking the coordinate origin O as the reference point target. The results are displayed in Figure 2. According to Figure 2, the fourth-order approximation Different from airborne vehicles that fly along a certain direction with constant velocities, the missile is in the dive phase, with two-dimensional acceleration, and its velocity magnitude and direction change quickly, with the first-order to fourth-order term range cell migration being serious. Therefore, it is important to take the efficiency approximation of each order term into account. Some simulations are presented using the parameters listed in Table 1, taking the coordinate origin O as the reference point target. The results are displayed in Figure 2. According to Figure 2, the fourth-order approximation error was only 1.313 m, while the highest approximation error when keeping the terms up to the third-order term was around 7.125 m. The sum of higher order terms above the fourth order was only 0.4687 m, which had little influence on the imaging effect and can be neglected. Therefore, in order to ensure good high precision imaging, the Taylor expansion series was taken as the fourth-order term. error was only 1.313 m, while the highest approximation error when keeping the terms up to the third-order term was around 7.125 m. The sum of higher order terms above the fourth order was only 0.4687 m, which had little influence on the imaging effect and can be neglected. Therefore, in order to ensure good high precision imaging, the Taylor expansion series was taken as the fourth-order term.

The Specific Steps of the Algorithm
The imaging method proposed in the paper was divided into the following four main steps. First, the Doppler frequency correction was performed in the range frequency azimuth time to obtain a corresponding spectrum in the next calculation. Second, the 2-D spectrum of spaceborne missile bistatic SAR echo signal was derived. Third, the operations in range complete range cell migration correction, range compression and determine the range direction position of the imaging scene. Finally, a modified NLCS algorithm was introduced to compensate for the azimuth FM rate, complete azimuth compression, and obtain the final image. According to Figure 3, the flowchart for the proposed algorithm

The Specific Steps of the Algorithm
The imaging method proposed in the paper was divided into the following four main steps. First, the Doppler frequency correction was performed in the range frequency azimuth time to obtain a corresponding spectrum in the next calculation. Second, the 2-D spectrum of spaceborne missile bistatic SAR echo signal was derived. Third, the operations in range complete range cell migration correction, range compression and determine the range direction position of the imaging scene. Finally, a modified NLCS algorithm was introduced to compensate for the azimuth FM rate, complete azimuth compression, and obtain the final image. According to Figure 3, the flowchart for the proposed algorithm was summarized to make each step of the process clear (the steps of algorithm improvement are marked with yellow boxes). The specific operations and derivations for each step are described in detail below.

Doppler Frequency Correction
Applying the range Fourier transform (FT) to (1) yields where f r represents the frequency domain variable corresponding to τ. Due to factors such as the missile's own acceleration and forward-looking imaging that cause Doppler central frequency deviation, for later calculation in the two-dimensional spectrum, the spectrum needs to be shifted to the zero-frequency position. It is assumed that the motion parameters of the satellite and missile can be calculated accurately during operation. According to Formula (6), the Doppler central frequency is calculated as

Doppler Frequency Correction
Applying the range Fourier transform (FT) to (1) yields where r f represents the frequency domain variable corresponding to τ . Due to factors such as the missile's own acceleration and forward-looking imaging that cause Doppler central frequency deviation, for later calculation in the two-dimensional spectrum, the spectrum needs to be shifted to the zero-frequency position. It is assumed that the motion parameters of the satellite and missile can be calculated accurately during operation. According to Formula(6), the Doppler central frequency is calculated as The Doppler frequency deviation term − needs to be corrected, and the Doppler frequency correction factor is Remote Sens. 2022, 14, 5770 7 of 23 multiplying (7) and (5) in the range frequency and azimuth time domain yields

Derivation of the Signal 2-D Spectrum
Application of the azimuth Fourier transform to (8) needs to use the method of series reversion [20]. Assuming where Then, applying the azimuth Fourier transform to (9) yields where f a represents the frequency domain variable correspond to η. Substituting (10) into (11), assuming applying the series reversion to (14) and inverting the power series yield where Substituting (15) into (11), we obtain the 2-D spectrum of S A (τ, η). where According to the shift properties and FFT skew, we have then, we obtain the original 2-D point target spectrum where According to the properties of series convergence [34] substituting (22) into (21) and expanding (21) into its Taylor series of f r at f r = 0 yields where Inspecting (23), ϕ 0 ( f a ) represents the azimuth modulation, which has no correlation with range variables. ϕ 1 ( f a ) is the range frequency linear term coefficient, which is related to large range cell migration and can be matched with the corresponding filter to achieve linear range migration correction. The first component in ϕ 2 ( f a ) represents the range modulation, and the remainder of components are the secondary range compression (SRC) terms that represent the range-azimuth coupling. ϕ 3 ( f a ) is the cubic term of the range frequency, which represents the rate of change of the range chirp rate. Because of the range-dependent characteristic of the above formula, it is difficult to figure out an accurate compensation for all targets at different ranges. To solve this problem, it is necessary to analyze the range-dependent characteristic of the residual range cell migration (RCM) and SRC. Because the effect of higher order terms above, the third order is very small, and only the third term coefficients were considered.
The three points, T 1 , T 6 , and T 16 , are located at the position of range directions 0, 150 m, and 300 m, respectively, as shown in Figure 4. The range-dependent characteristic evaluation of the SRC and the residual RCM are shown in Figure 5. Observing the images, it can be seen that the SRC phase error and the residual RCM error were smaller than π/4 and one range resolution cells, respectively. Thus, the imaging results would not be affected by those effects of the range dependence of the residual SRC and RCM. The three points, T1, T6, and T16, are located at the position of range directions 0, 150 m, and 300 m, respectively, as shown in Figure 4. The range-dependent characteristic evaluation of the SRC and the residual RCM are shown in Figure 5. Observing the images, it can be seen that the SRC phase error and the residual RCM error were smaller than / 4 π and one range resolution cells, respectively. Thus, the imaging results would not be affected by those effects of the range dependence of the residual SRC and RCM.

The Operations in Range
In order to keep the range migration error within 1/4 range sampling, we performed bulk range cell migration correction in the 2-D frequency domain by multiplying with the conjugate of ϕ 1 ( f a ), setting the scene center as the reference point target. The factor was given by The range compression term is exp(−jϕ 2 ( f a ) f 2 r ). Because the missile was in the highfrequency band, only the first component term in ϕ 2 ( f a ) was kept, and the remainder terms in ϕ 2 ( f a ). were ignored, having little impact on range compression. The range compression factor is expressed as The echo signal of any point target in the scene will become a vertical line in the azimuth after range migration correction and range pulse compression. The linear term was eliminated after the range migration correction. In order to determine the range direction position of the imaging scene, the function was multiplied by Multiplying (20) with (28), (29), (30) and transforming the result into range time azimuth frequency domain by the range inverse Fourier transform, the range focused signal can be written as where To illustrate the range processing of the aforementioned algorithm, some simulations on the echo trajectories and 2-D spectrum of the simulated targets are presented before and after the range operation, and the results are shown in Figure 6. The simulation parameters are shown in Table 1. Figure 6a shows the original signal in the time domain, Figure 6b represents the original signal in the two-dimensional frequency domain, Figure 6d is the image of the signal after range compression, Figure 6e is the image of the signal after range compression and migration correction, and Figure 6c shows the two-dimensional spectrum after range migration correction. Combined with the original signal echo in Figure 6a,d, the original signal echo was compressed into nine slanted straight lines after range compression, and then the slanted lines were straightened by range migration correction, as shown in Figure 6e. The spectrum was also corrected accordingly, as shown in Figure 6b

Azimuth Compression by a Modified Bistatic Azimuth NLCS Algorithm
Different from the bistatic configuration with constant speed [28][29][30], after range pulse compression and range migration correction, the azimuth FM rate of the targets located in different azimuth directions of the same range unit changed significantly and could not be regarded as a constant (specific simulation analysis was conducted later). If the azimuthal direction is compressed directly by the azimuthal matching filter function, which is constructed by the azimuth FM rate of the azimuthal center, the edge target points of the imaging scene will be difficult to focus. Therefore, before the azimuth compression, the variation of azimuth FM rate along the azimuth direction must be compensated. The azimuth nonlinear scaling algorithm was introduced here to improve the effect of azimuth focusing.
The linear term coefficient of the azimuth frequency in the signal can be written as ( ) g y , which represents the azimuthal position of the target point.

Azimuth Compression by a Modified Bistatic Azimuth NLCS Algorithm
Different from the bistatic configuration with constant speed [28][29][30], after range pulse compression and range migration correction, the azimuth FM rate of the targets located in different azimuth directions of the same range unit changed significantly and could not be regarded as a constant (specific simulation analysis was conducted later). If the azimuthal direction is compressed directly by the azimuthal matching filter function, which is constructed by the azimuth FM rate of the azimuthal center, the edge target points of the imaging scene will be difficult to focus. Therefore, before the azimuth compres-sion, the variation of azimuth FM rate along the azimuth direction must be compensated. The azimuth nonlinear scaling algorithm was introduced here to improve the effect of azimuth focusing.
The linear term coefficient of the azimuth frequency in the signal can be written as g(y), which represents the azimuthal position of the target point.
Expanding (32) into its Taylor series of y at y = 0 and keeping up to the first-order term yield g(y) ≈ g 0 + ρy, where g 0 represents the whole azimuth displacement, which is easily corrected by multiplying with its conjugate. ρy reflects the compression or stretching of the image along the azimuth direction, which is caused by the acceleration and height speed. Due to the acceleration of the missile, the speed changes rapidly, making ρy = y/v y . In order to correct ρy and make ρy = y/v y , the nonlinear chirp scaling method proposed later was adopted to eliminate the influence on the image. First, the whole azimuth displacement g 0 was corrected by its conjugate factor as follows.
H 3 ( f a ) = exp(−j2πg 0 f a ), (36) multiplying (31) and (36) yields where Due to the influence of missile-borne acceleration, the azimuth FM rate mentioned above will change with the azimuth position y, as shown in Figure 7a. Considering the azimuth FM rate in the center of the scene as the whole azimuthal FM rate, it will dramatically affect the imaging focus on both sides of the scene edges. Therefore, it is very essential to perform a deep analysis in the azimuth-dependent characteristic of the azimuth FM rate. Expanding K a into its Taylor series of ∆η at ∆η = 0 yields where azimuth FM rate in the center of the scene as the whole azimuthal FM rate, it will dramatically affect the imaging focus on both sides of the scene edges. Therefore, it is very essential to perform a deep analysis in the azimuth-dependent characteristic of the azimuth FM rate. Expanding a K into its Taylor series of η Δ at 0 where  The third-and fourth-term coefficients 1 Z and 2 Z of a f in (37) had a small effect with the change of azimuth, as shown in Figure 7b, and for the convenience of calculation, the amplitude of variation can be ignored. Therefore, the coefficient amplitude of the whole azimuthal direction can be considered as the coefficient amplitude of the azimuthal center. Before the azimuth FM rate compensation and nonlinear chirp scaling operation, in order to eliminate the influence of the third-order and fourth-order terms of azimuth frequency a f , it is necessary to introduce the third-order and fourth-order filters in the range time azimuth frequency domain. The third-order and fourth-order were are as follows: where A and B are the uncertain parameters, which were calculated later. Multiplying (37) with (39) yields The third-and fourth-term coefficients Z 1 and Z 2 of f a in (37) had a small effect with the change of azimuth, as shown in Figure 7b, and for the convenience of calculation, the amplitude of variation can be ignored. Therefore, the coefficient amplitude of the whole azimuthal direction can be considered as the coefficient amplitude of the azimuthal center. Before the azimuth FM rate compensation and nonlinear chirp scaling operation, in order to eliminate the influence of the third-order and fourth-order terms of azimuth frequency f a , it is necessary to introduce the third-order and fourth-order filters in the range time azimuth frequency domain. The third-order and fourth-order were are as follows: where A and B are the uncertain parameters, which were calculated later. where (41) Applying the principle of stationary phase (POSP) to (40) and transforming (40) into the azimuth time domain yields The next step is to compensate the azimuth FM rate so that the azimuth FM rate of each target point is equal to the azimuth FM rate of the scene center. According to the algorithm of nonlinear chirp scaling, a fourth-order chirp scaling factor H 5 was constructed to multiply the above formula to meet the above requirements.
where l 2 , l 3 , l 4 are the uncertain coefficients. Multiplying (42) with (43) and transforming the result into the azimuth frequency domain by applying the azimuth Fourier transform yields expressing the φ( f a ) as a power series of ∆η and f a yields where As can be seen from (45), the first term contained only azimuth frequencies, which can be removed directly by constructing a conjugate-matching filter function during azimuth compression, while the second term is a linear term of azimuthal frequency, which contains the position information of the target imaging azimuth and cannot be removed. The remaining items represent the azimuth-dependent terms that degrade the effect of azimuth focusing. To compensate for the azimuth FM rate and eliminate azimuth-dependent modulation, making its value not change with azimuth position, we set c 2 (l 2 ) to −2π/α(α = ρv y ) and c 3 , c 4 , c 5 , c 6 , c 7 to 0. The formulas can be solved as follows.
substituting (47) into (44) yields Inspecting (48), The second, third, and fourth term coefficients of azimuth frequency are constant and no longer vary with azimuthal change, and thus the compensation correction was completed and the azimuth-dependent characteristic of the azimuth FM rate was removed entirely. Therefore, the azimuth direction can be compressed directly, and the azimuth matching filter function can be written as multiplying (48) with (49) yields Transforming (50) into the range time azimuth time domain by azimuth inverse FT yields the final SAR image.

Simulation Results by the Proposed Algorithm
As shown in Figure 4, the simulations chose an array of 5 × 5 point targets at intervals of 100 m and 150 m in the azimuth and range axes, respectively. T 1 is the target point at the center of the scene, and T 5 is the edge point with the most severe azimuthal variation. The system parameters used in the simulation are shown in Table 1. The results of the point target imaging simulation are shown in Figure 8a. Due to the influence of forward-looking imaging, there was a certain geometric distortion in the range and azimuth direction, which was not consistent with the real target position. The results of the distortion correction are shown in Figure 8b.
Transforming (50) into the range time azimuth time domain by azimuth inverse FT yields the final SAR image.

Simulation Results by the Proposed Algorithm
As shown in Figure 4, the simulations chose an array of 5 × 5 point targets at intervals of 100 m and 150 m in the azimuth and range axes, respectively. T1 is the target point at the center of the scene, and T5 is the edge point with the most severe azimuthal variation. The system parameters used in the simulation are shown in Table 1. The results of the point target imaging simulation are shown in Figure 8a. Due to the influence of forwardlooking imaging, there was a certain geometric distortion in the range and azimuth direction, which was not consistent with the real target position. The results of the distortion correction are shown in Figure 8b.

Comparison of Algorithm Efficiency
As can be seen from the flow chart, the increase in the computation of the algorithm mainly came from the azimuth nonlinear scaling operations, including one FFT, one IFFT, and five times phase multiplication (as shown in the azimuth-processing section of the flow chart). According to [6], a complex phase multiplication takes six FLOPs, and an Npoint FFT or IFFT takes

Comparison of Algorithm Efficiency
As can be seen from the flow chart, the increase in the computation of the algorithm mainly came from the azimuth nonlinear scaling operations, including one FFT, one IFFT, and five times phase multiplication (as shown in the azimuth-processing section of the flow chart). According to [6], a complex phase multiplication takes six FLOPs, and an N-point FFT or IFFT takes 5N log 2 (N) FLOPs [35]. In the simulation experiment, the number of sampling points in the range direction was N f ast = 2048, the number of sampling points in the azimuth direction was N slow = 512, and the number of output points selected after range compression was N r = 1856. It can be calculated that the added floating point arithmetic amount after the azimuth nonlinear operation was ∆Q = 2 × 5N r N slow log 2 (N slow ) + 5 × 6N r N slow = 0.950272 × 10 6 (FLOPs) Although the calculation amount was slightly increased compared with that without considering the space variance of azimuth FM rate, the azimuthal focusing was obviously improved and the image effect was better.
The computational complexity of the traditional NLCS algorithm in [36] was (setting the interpolation kernel value to 6) the total computation of the proposed algorithm in this paper is It can be seen that the computational amount of the algorithm proposed in this paper was better than the traditional NLCS algorithm in [36], and the algorithm in this paper had no time-consuming interpolation operation, and thus its processing efficiency was higher. Moreover, the azimuth focusing was significantly improved, and the image effect was better (a detailed comparison is made in the next section).

Comparison of Imaging Effects
In order to evaluate the imaging performance of the proposed method, the traditional azimuth NLCS algorithm in [36] was introduced here to process the same simulated data, which solely considered the linear component of azimuth FM rate and disregarded the effects of residual RCM. Figure 9b shows the imaging result obtained by the traditional azimuth NLCS algorithm. By comparing Figure 9a,b, it can be seen that the traditional azimuth NLCS algorithm failed to achieve the same focus on the azimuth direction as the proposed azimuth NLCS algorithm because of the mismatched azimuth FM rate. The imaging result of point 1, point 3, and point 5 that were positioned at the same range but in different azimuth directions are shown in Figure 9a,b by using the traditional azimuth NLCS algorithm and the proposed azimuth NLCS algorithm, respectively. By comparing Figure 9a,b, it can be seen that in Figure 9a, the imaging effect was good with the change of azimuth distance from the center of the scene to the edge of the scene, while in Figure 9b, the image quality in the azimuth direction decreased rapidly. The specific performance comparison of imaging results is shown in Figure 10. Furthermore, Table 2 shows the calculated PSLR and ISLR for the corresponding targets in order to evaluate the imaging performance of the proposed NLCS algorithm. According to the Table 2, the proposed NLCS algorithm achieved nearly theoretical PSLR and ISLR, whereas the traditional NLCS was only able to achieve focused performance in the range direction, and in azimuth direction, the measured parameters were poorer than the theoretical values.

Conclusions
In this paper, a modified NLCS algorithm for the spaceborne missile bistatic SAR data processing is presented. Firstly, the system configuration and signal model were established, and then the Doppler frequency correction was carried out to eliminate Doppler center frequency deviation, which was caused by the acceleration of the missile itself and the forward-looking imaging. Next, the two-dimensional spectrum expression was derived by the method of series reversion, and range migration correction and range compression were performed in the two-dimensional frequency domain. Then, a modified azimuth NLCS operation was applied to eliminate the azimuth dependence of Doppler parameters through the use of a higher order approximation to the cubic phase term and

Conclusions
In this paper, a modified NLCS algorithm for the spaceborne missile bistatic SAR data processing is presented. Firstly, the system configuration and signal model were established, and then the Doppler frequency correction was carried out to eliminate Doppler center frequency deviation, which was caused by the acceleration of the missile itself and the forward-looking imaging. Next, the two-dimensional spectrum expression was derived by the method of series reversion, and range migration correction and range compression were performed in the two-dimensional frequency domain. Then, a modified azimuth NLCS operation was applied to eliminate the azimuth dependence of Doppler parameters through the use of a higher order approximation to the cubic phase term and azimuth FM rate. The modified azimuth NLCS operation equalized the azimuth FM rate so that azimuth compression can be achieved. Finally, the high efficiency and feasibility of the algorithm were proved by comparing the imaging effect, computational efficiency, and performance analysis between different algorithms.