Focusing Translational Variant Bistatic Forward-Looking SAR Using Keystone Transform and Extended Nonlinear Chirp Scaling

School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China; threadkite@sina.com (Z.S.); heidy_vic@sina.com (Z.L.); yulinhuang@uestc.edu.cn (Y.H.); jyyang@uestc.edu.cn (J.Y.); sunzhichao1110@sina.com (Z.L.) * Correspondence: junjie_wu@uestc.edu.cn; Tel.: +86-28-6183-1533 † These authors contributed equally to this work. ‡ Current address: School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China.


Introduction
In recent years, synthetic aperture radar (SAR) has become increasingly attractive in both civilian and military fields, because it can provide high-resolution day-and night-time images of the observed area, independent of weather conditions [1][2][3][4].Typical fields of SAR applications include disaster monitoring, resources exploration, geological mapping, military reconnaissance, and so on.However, the traditional monostatic SAR cannot achieve forward-looking imaging because of its inherent characteristic [5].The range and azimuth directions coincide with each other, and a two-dimensional image can not be obtained.This restricts the applications of SAR technology in areas such as navigation, terminal guidance, self-landing, etc.With the choice of proper geometry, bistatic SAR can image the forward-looking terrain, which is called bistatic forward-looking SAR (BFSAR).A bistatic SAR system can break through the limitations of traditional monostatic SAR in forward-looking imaging and achieve high-resolution imaging of the forward-looking terrain of the receiver.This might result in some new applications of bistatic SAR in aviation by imaging obstacles (e.g., mountains) in the flight direction which are covered by fog, rain, and clouds.With such systems, we can also improve the safety of landing airplanes by imaging the runway in undesirable weather conditions [6].Current research and publications on BFSAR include imaging theories [7][8][9], imaging algorithms [10][11][12][13][14][15][16], and experiments [17][18][19][20][21].In Wu J. et al. [7], BFSAR is categorized into three types: translational invariant (TI), translational variant (TV), and stationary transmitter cases.In TV mode, the transmitter and receiver move with different velocities.The spaceborne/airborne BFSAR and the airborne/missleborne BFSAR often work in this mode.In comparison with the TI mode, TV-BFSAR provides higher flexibility and mobility, which makes it suitable for more complicated applications.
Some reconstruction methods for BFSAR, such as the range-Doppler algorithm (RDA) [10][11][12], Omega-k [13,14], and back-projection [15] are ongoing.The RDA neglects the spatial variation of the secondary range compression (SRC), which is very important for BFSAR.As for the Omega-k for BFSAR, the current methods all use a monostatic Omega-k algorithm after a data conversion from bistatic to monostatic.This makes it fail for cases with large bistatic angle.In addition, because of the symmetric effect of the receiver about its flight path, the range cell migration (RCM) is severely nonlinear with cross-track coordinate.Thus, the Omega-k algorithm is actually not suitable for TV-BFSAR [22].A back-projection algorithm is a precise way to focus the raw data of TV-BFSAR.However, because of its low computational efficiency, the applications of BFSAR in real-time navigation and reconnaissance are limited, which are the major application fields of BFSAR.
In TV-BFSAR, the RCM varies with the positions of the targets, and the linear component of RCM is the most significant which should be handled.Traditional range cell migration correction (RCMC) methods [23][24][25] operated in range frequency and azimuth time domains all use correction factor multiplication to remove the linear component of RCM.However, the linear RCMC factor is constructed by using the linear RCM of the reference target, and the spatial variance of linear RCM is neglected.So, residual linear RCM still exists for targets away from the reference target.Keystone Transform (KT) is a data reformatting technique that is generally used to correct the range walk [26] or range curvature [27] of moving targets.Because it is independent of the motion parameters, it can correct the range walk or range curvature of all targets in the scene without knowing their velocities, regardless of the target position.First-order KT has already been applied to the stationary scene data processing of bistatic SAR to remove the linear RCM in [28][29][30], but these methods neglect the residual higher-order RCM, which degrades the azimuth focusing performance.Another problem is the 2-D variation of Doppler coefficients.Different from TI-BFSAR, the Doppler coefficients in the TV case depend on the range and azimuth positions of the targets.Moreover, the target undergoes a position shift after RCMC operated in range frequency and azimuth time domain, and it further aggravates the 2-D variation.Thus, the Doppler coefficients of the targets in the same range gate should first be equalized so that these targets can be compressed with the same azimuth-matched filter.Nonlinear chirp scaling (NLCS) was first adopted by Wong and Yeo in monostatic and bistatic SAR data processing [23].Current publications on bistatic SAR NLCS algorithm can be found in References [31][32][33][34].In Qiu et al. [34], an azimuth NLCS algorithm based on numerical fitting is proposed to process bistatic SAR data with a stationary transmitter.However, this algorithm can only deal with the spatial variance of Doppler frequency modulation rate.In some recent publications on highly squint monostatic SAR data focusing [24,25], an extended azimuth nonlinear chirp scaling algorithm (ENLCS) has been proposed.The ENLCS can simultaneously equalize the Doppler FM rate and third-order coefficient.However, in TV-BFSAR, Doppler centroid is spatial-variant for targets in the same range gate after RCMC, and it induces a variation term into the Doppler FM rate, which is ignored in the algorithm.This deteriorates the performance of azimuth focusing.Moreover, if the spatial variation of the Doppler centroid is not removed, the magnitude of the targets away from the scene center is diminished, and the azimuth resolution is coarsened because the reference function mismatch the echo data in azimuth frequency domain.
In this paper, we first analyze the characteristic of the echo signal of TV-BFSAR.Based on the analysis, an extended azimuth NLCS algorithm is proposed to process the raw data and obtain a well-focused image of the TV-BFSAR.The algorithm aims to deal with two problems of TV-BFSAR: 2D spatial-variation of linear RCM and azimuth-variant Doppler parameters.The processing procedures are as follows.First, KT is performed in the range frequency and azimuth time domain to remove the spatial-variant linear RCM.Second, a correction factor is multiplied to the echo to compensate for the residual RCM.Then, the azimuth dependence of Doppler character and the influence of spatial-variant Doppler centroid on azimuth compression are analyzed.Polynomial approximations of Doppler centroid and FM rate are obtained by using numerical curve fitting.At last, a fourth-order filtering in the azimuth frequency domain together with NLCS is carried out to equalize the Doppler centroid and FM rate based on the above analysis.The proposed algorithm resolves the problem of the variation of linear RCM and Doppler centroid, and thus has a better imaging performance for multiple targets in the scene.This method is computationally efficient compared to time-domain methods such as back projection (BP) and achieves better focusing performance than other frequency domain methods.
The rest of the paper is organized as follows.Section 2 gives the geometric and signal model of TV-BFSAR.Section 3 presents the range processing of echo signal, including KT, range compression, and higher-order RCMC.In Section 4, the azimuth dependence of Doppler coefficients is first analyzed, and the extended NLCS is then derived in detail.The numerical simulations are given in Section 5. Section 6 presents some discussionof the proposed method, including advantages, practical limitations, and future work.Section 7 concludes the paper.

Signal Model
The geometric graph in Figure 1 shows the configuration of a general translational variant bistatic forward-looking SAR.The transmitter works in squint mode and the receiver works in forward-looking mode.The receiver travels along the y axis with a velocity of V R , and the velocity of the transmitter is V T .The transmitter and the receiver both work in stripmap mode.The angle between the flight direction of the transmitter and the receiver is α.The original coordinates of the transmitter and the receiver are (x T , y T , h T ) and (x R , y R , h R ), respectively.The downward-looking angle of the receiver is defined as θ R and the squint angle of the transmitter is defined as θ T , respectively.The origin of the coordinate O is the scene center, and P(x, y) is an arbitrary point target in the scene.
Imaging geometry configuration of translational variant bistatic forward-looking synthetic aperture radar (TV-BFSAR).
Supposing the transmitted signal is an linear frequency modulated (LFM) pulse, the received signal after demodulation is where ω r and ω a are the range and azimuth envelopes, respectively.τ is the fast time and t is the slow time.λ is the wavelength and λ = c/ f c .c is the speed of light and f c is the carrier frequency.K r represents the range FM rate, T a is the synthetic aperture time, and t R is the azimuth position of the target.R(t; x, y) is the bistatic range history, which is given by The range histories with respect to target P(x, y) of the receiver and the transmitter are: where r 0R and r 0T are the slant ranges of target P to the receiver and transmitter at t R and t T , respectively.t R = y/V R and t T = y/V T cos α are the beam center crossing time of the receiver and the transmitter, respectively.From Figure 1, r 0R and r 0T are given by After range Fourier transform, the signal in the range frequency domain is given by: Expanding R R (t; x, y) and R T (t; x, y) (jointly denoted as R R,T (t; x, y) for simplicity) at t R and t T (jointly denoted as t R,T for simplicity) into Taylor series and keeping up to third-order terms yields where A R,T , B R,T , and C R,T are the expanding coefficients.
In bistatic forward-looking SAR, the absolute values of B R,T and C R,T are much smaller than that of A R,T .Thus, the RCM components caused by B and C are much smaller than that of A, and linear range cell migration (LRCM) is more significant than higher-order RCM, which will be illustrated in Section 5.

LRCMC via Keystone Transform and Higher-Order RCMC
In this section, we consider the problem of RCM and range compression.As mentioned above, the main part of range cell migration is linear RCM.Thus, KT is first introduced to deal with the spatial-variant LRCM, and range compression is carried out after KT.Then, to compensate for the residual RCM, a third-order phase term of azimuth time variable is multiplied to the echo.

Linear RCMC
Substituting (Equation ( 8)) into Equation (2), R(t; x, y) is rewritten as follows where R(0; x, y) represents the bistatic slant range at t = 0.After the slant range expansion in Equation ( 12), the phase term of the received echo after range fast Fourier transform (FFT) is given by The first term is range modulation term, while the second term contains azimuth information and range-azimuth coupling.Notice that the second phase term is a linear phase of range frequency variable f τ , so it represents a position shift of echo signal in range direction in the t versus τ plane.We can see that the shift is a third-order polynomial of t.From the derivation above, R(0; x, y), A , B , and C change with the coordinates of the target P(x, y).So, the RCM of the targets in this configuration is 2-D spatially variant.If we construct a compensation factor in the range frequency and azimuth time domain using the linear RCM of the scene center [24], the LRCM of the targets away from the scene center cannot be fully compensated for, because it ignores the spatial variance of linear RCM.Moreover, the further the distance between a target and the scene center is, the larger the residual RCM of the target.Thus, in order to remove the spatially-variant linear RCM, KT is performed where t m is the slow-time variable after KT.The phase term after KT is given by As can be seen in Equation ( 19), the first-order polynomial of t m is −2πA f c t m /c.Thus, the first-order coupling of range frequency variable f τ and azimuth time variable t m is removed.KT is essentially a coordinate transformation, and it changes the original azimuth time variable t to a new azimuth time variable t m .Figure 2 illustrates the 2-D spectrum of simulated SAR data before and after KT.It is obvious that the spectrum before KT is skewed, and this phenomenon is mainly caused by first-order coupling of echo signal.After KT, the skew problem is corrected, and azimuth and range direction become orthogonal to each other.In addition, before KT, the range cell of target P(x, y) is After KT, however, the targets with the same R (0; x, y) are in the same range gate.That is to say, KT gathers those targets having the same R(0; x, y) into the same range bin.
Expanding the above phase term Equation ( 19) into a Taylor series of f τ and keeping up to second-order term yields Inspecting Equation ( 21), the phase term in the first bracket is the azimuth modulation term and is independent of range variable f τ .The second term is the first-order term of range frequency represents the RCM of the target after range compression.It is obvious that the range gate of target P(x, y) is R(0; x, y) and the residual RCM still exists, which is to be compensated for in the next step.The third term is the range modulation term, and the range FM rate is changed.The new range FM rate is Thus, the range compression should be carried out after KT.As analyzed above, B and C are fairly small, and we can neglect their spatial changes in TV-BFSAR.So, here the range compression is performed using the K of the scene center: The range compression factor is given by

Higher-Order RCMC
After KT and range compression, the phase term in the range frequency domain is Compared with Equation ( 13), the first brackets in Equation ( 25) can be written as R (t m ; x, y).
The echo is then transformed to the 2-D time domain by applying the Principle of Stationary Phase As mentioned above, the residual RCM is RCM res is a third-order polynomial of azimuth time variable t m .Thus, we can construct a phase factor in the time domain to compensate for the residual RCM, which is given by where the spatial variations of B and C' are neglected.The coefficient B (0, 0)t 2 m /2 + C (0, 0)t 3 m /3 of the exponential phase term above is equal to the residual RCM of the scene center.
After higher-order RCMC, the echo in the 2-D time domain is Inspecting Equation ( 28), R(0; x, y) is the new slant range of echo signal after KT.The second and third terms are azimuth envelope and phase terms, respectively.We can find that the slow-time modulation in Equation ( 28) for a single target has not been changed after KT.

Azimuth Compression via Extended NLCS
In previous publications on azimuth NLCS processing, only the spatial variation of the Doppler FM rate and third-order terms are taken into consideration, and the variation of the Doppler centroid is neglected.However, in TV-BFSAR, the first-order expanding coefficient of slant range (i.e., A ) is 2-D variant.This variation brings two problems.First, the linear RCM is spatially-variant along both azimuth and slant range directions.This problem is solved in Section 3. Second, the Doppler centroid is essentially 2-D spatially-variant.Furthermore, the targets undergo a position shift along the range direction due to KT, which aggravates the variation.In this section, an azimuth-extended NLCS process is derived to deal with the problem.
As previously explained, those targets having the same R(0; x, y) fall into the same range bin after KT with different x.Therefore, they have different Doppler centroids, FM rates, and slow-time reference functions.Thus, targets in the same range gate cannot be compressed with the same matched filter unless reference functions are preprocessed to be the same.Based on the analysis in [28], the representation of the x coordinate in each range gate can be expressed as where Then, substituting (x(R, y), y) into Equation ( 17), we can compute the Doppler coefficients of every point target whose energy center corresponds to the sampling location in the dataset after KT.

Azimuth-Dependent Doppler Coefficients Analysis
First, we expand the azimuth phase term at t = t R and keep up to the third-order term as follows: where f dc , f dr , and f d3 are the Doppler centroid, FM rate, and third-order coefficient, respectively.Now, we shall evaluate the spatial variance of the Doppler centroid and analyze its influence on azimuth compression.f dc is given by where cos θ R and cos θ T are given by cos θ T = sin αx − cos αy where x = x − x T + V T sin αy/V R and y = y − y T − V T cos αy/V R are the contributions of the transmitter to the Doppler centroid along the cross-track direction and the azimuth direction, respectively.Inspecting Equations ( 31)- (33), it is obvious that the Doppler centroid of a target varies with the cross-range coordinate.First of all, generally speaking, in order to suppress the sidelobes and avoid azimuth ambiguity, the azimuth reference function undergoes a truncation operation so that its width in the time and frequency domains are both limited.If the azimuth variation of the Doppler centroid is not removed, the echo data of these targets are not in the same azimuth position in the frequency domain.Thus, some of the echo information of the targets away from the azimuth center is missed because of the reference function multiplication.The effective Doppler bandwidth is decreased, and the azimuth resolution is coarsened.The magnitudes of these targets are diminished compared with the azimuth center.When the reference function totally mismatches the echo data of the targets in the frequency domain, the targets will vanish from the scene.Moreover, if we construct the azimuth reference function using the Doppler centroid of the reference target, the imaging performance of the targets with different Doppler centroid is deteriorated.The azimuth reference function in the time domain is given by: Transforming Equation ( 34) into the frequency domain by principle of stationary phase (POSP), the phase term is given by The first term is a constant phase term, which should not exceed π/4.The second phase term represents first-order term of azimuth frequency, which is related to target position after azimuth compression.The third term is the second-order term of azimuth frequency, representing azimuth modulation.The last term is related to the sidelobes of the focused image.We can see from the above that the Doppler centroid influences the azimuth FM rate, first-order azimuth frequency, and the constant phase term.The spatial variation of the Doppler centroid induces a variant term into the Doppler FM rate.Thus, the Doppler FM rate is changed because of the difference of Doppler centroid between an arbitrary target and the reference target.The new FM rate is given by Here, we model the azimuth dependence of the Doppler centroid as a first-order polynomial of azimuth position.
where f dc0 is the Doppler centroid of the reference target and a is the coefficient of the first-order term of azimuth position t R .
Similarly, the Doppler FM rate is modeled as a second-order polynomial of azimuth position, as follows: where f dr0 is the Doppler FM rate of the reference target, and b and c are the coefficients of first-order and second-order term of azimuth position t R .The quadratic and cubic phase errors induced by the approximation are given by where ∆ f dr is the second-order approximation error of Doppler FM rate and ∆ f dr is where ∆ f dc is the first-order approximation error of the Doppler centroid.

Derivation of the Extended NLCS
In the previous subsection, the character of the spatial variation of f dc and f dr is obtained.Now, the extended NLCS is derived based on the presented analysis.Let us assume that the azimuth sampling rate is pulse repetition frequency (PRF) and the Doppler centroid of a target in a given range gate is f dc ; then, we have where f dc is the baseband Doppler centroid and n is the ambiguity number.If f dc = 0, the Doppler spectrum of the target is not at the center of PRF, which is different from broadside SAR.Due to the 2-D variation of the Doppler centroid, f dc is different among targets in the same range gate.This difference broadens the Doppler spectrum.If f dc + B a /2 > PRF/2, the echo signal is aliased, where B a is the total Doppler bandwidth.To avoid the Doppler ambiguity, an azimuth phase term is multiplied to Equation (30), which is given by After the phase multiplication, the Doppler spectrum is moved to the center of PRF, and part of the coupling term of t m and f τ is removed.The azimuth phase is then given by First, the data after phase multiplication is transformed into the range Doppler (RD) domain.Then, a fourth-order filtering is performed where Y 3 and Y 4 are the coefficients of H F ( f ).The fourth-order filtering can introduce two degrees of freedom in the azimuth phase term, which are used for the derivation of the ENLCS.The azimuth phase in the frequency domain is Equation ( 46) is then transformed into the time domain: After the fourth-order filtering, an azimuth NLCS factor is introduced to equalize the spatially-variant Doppler centroid as well as the FM rate.The NLCS factor is given as follows: where q 2 , q 3 , and q 4 are the coefficients of the NLCS operator.The azimuth phase in the time domain after fourth-order filtering is then multiplied by the NLCS factor and transformed into the frequency domain: (q 2 + f dr ) 4  (49) Using the approximation presented in Equations ( 37) and (38), we expand the above phase term into a power series of t R and f a where E, F, G, H, and I are the expanding coefficients of the coupling term of t R and f a , which are given in Equation (51).Inspecting Equation (50), the first term D ( f a ) represents the azimuth modulation term, which is independent of azimuth position.The second term is the azimuth position of the targets after compression.The third term causes a nonlinear shift of azimuth position.The fourth and fifth terms represent the spatial variation of Doppler FM rate along the azimuth direction.The sixth term is the spatial variation of the third-order coefficient of Doppler frequency, which influences the sidelobes of the target after azimuth compression.The seventh term is independent of azimuth frequency, and can be neglected.
In order to eliminate the azimuth variation of Doppler coefficients, the coefficient of the first-order coupling of f a and t R is set to −π/β, where β is a constant scaling factor which determines the azimuth position of the targets.Furthermore, the coefficients of other coupling terms is set to zero.Thus, we obtain an equation system with five unknowns: Solving Equation (56), the parameters of fourth-order filtering and NLCS are given by where L, M, and N are given in Equations (62)−(64).
After the fourth-order filtering and NLCS process, the phase of the echo signal in the frequency domain is given by It is obvious that in Equation (65), the spatial variance of the Doppler centroid and the FM rate is removed.
The Doppler coefficients of all the targets in the same range gate are processed to be the same.Notice that the azimuth FM rate after the process is changed to f drnew = q 2 + f dr0 .When β equals 0.5, the target after azimuth compression is at its real position.While β is larger than 0.5, the azimuth distance among targets in the same range gate is enlarged.However, β cannot be too far from 0.5, because the larger β is, the larger is the absolute value of f drnew .Additionally, the Doppler bandwidth is widened, which may cause azimuth aliasing.On the other hand, β should not be too small, because the targets may get too close to each other.Thus, β is better chosen to be 0.5.Figure 3 shows the Doppler spectra of the simulated targets before and after the extended NLCS operation.It is clear that before extended NLCS, the Doppler spectra of different targets are located at different azimuth positions, due to 2-D variance of the Doppler centroid.After the extended NLCS operation, the Doppler spectra of all targets are moved to the center of the azimuth frequency, and the azimuth variation of the Doppler centroid is removed with respect to every range gate.The Doppler bandwidth of the three range gates 3350, 4096, and 4950 are 215.53Hz, 228.85 Hz, and 240.14 Hz, respectively.This phenomenon is due to their differences in azimuth processing.The coefficients of fourth-order filtering, NLCS, and the azimuth compression factor are calculated with respect to every range gate.Thus, the Doppler FM rates of the targets after azimuth NLCS are slightly different from each other.According to Equation (65), the azimuth compression function is given by (66) Figure 4 shows the flowchart of the proposed algorithm.In the derivation of the proposed algorithm, the azimuth dependence of the Doppler centroid and FM rate are approximated by a first-order polynomial function and a second-order polynomial function, respectively.Moreover, the spatial variation of the third-order Doppler coefficient is neglected.Based on the analysis in this paper, the approximations here are valid, and the corresponding errors are small enough that we can ignore them.However, these approximations can cause some errors in specific geometry configurations where the azimuth dependence of the Doppler character is more complicated, and the order of the polynomial function chosen in this paper is not enough.In addition, if we ignore the spatial variance of the third-order Doppler coefficient, the corresponding phase error may be larger than π/4, and it will influence the azimuth compression performance.Thus, the imaging quality can get worse.In such cases, we can increase the order of approximation of the polynomial function and take the spatial variance of third-order Doppler coefficient into consideration, but it will make the derivation much more complex.

Computation Complexity
From the above analysis, the proposed imaging algorithm includes two instances of range FFT and inverse FFT (IFFT) and four instances of azimuth FFT and IFFT.The major operations include a keystone transform and five instances of 2-D phase factor multiplication; i.e., range compression, second order RCMC, fourth order filtering, extended NLCS, and azimuth compression.Suppose that the range sample number is denoted as N r and the azimuth sample number is N a .The total number of real floating-point operations is where M ker is the number of neighbor samples used for the range and azimuth interpolation.Therefore, the computational complexity is of order O N 2 log 2 N , where N is the 1-D size of the data.It can be concluded that the computation complexity of the proposed method is of the same order as traditional range Doppler algorithm (RDA) and ωKA and more efficient than back-projection, whose computation complexity is of order O N 3 .

Numerical Simulation Verification
Numerical simulations are carried out in this subsection to verify the effectiveness of the proposed imaging algorithm.The simulation parameters are listed in Table 1.The receiver can be mounted on a fighter.In the following simulation, the locations of the 12 targets are shown in Figure 5.The azimuth distance between targets P1-P6 and scene center O are 350 m, and the azimuth distance between targets P7-P12 and the scene center are 175 m.The coordinates of targets P1-P6 are (−54.6576,350) m, (−443.4703,350) m, (−872.9143,350) m, (707.1249,−350) m, (381.2807,−350) m, and (32.1523, −350) m, respectively.Targets P7-P12 are the mid-points on the lines of P1-P6 and scene center O.The positions of the targets are chosen in this way so that targets P1 and P4, P2, O, and P5, P3, and P6 are in the same range gates after KT and higher-order RCMC, respectively.Thus, we can evaluate the effectiveness of the algorithm both in the azimuth and range directions.First of all, the expanding coefficients A R , B R , C R , and A T , B T , C T in the simulated case are 249.6151m/s, 3.8402 m/s 2 , 0.3988 m/s 3 , and −49.2518 m/s, 0.7537 m/s 2 , and −0.0111 m/s 3 , respectively.It can be seen that A R,T is much larger than B R,T and C R,T , and LRCM is the most significant in bistatic forward-looking SAR.
Figure 6 shows the residual RCM of the three simulated targets and their corresponding spatial variance.It can be seen that the residual RCMs of the simulated targets are approximately 10 m with 2 s synthetic aperture time.The residual RCMs exceed one resolution cell and should be compensated.However, the spatial variance of residual RCM between the targets is 0.1 m, which is smaller than one resolution cell and can be ignored.Figure 7 shows the echo of the seven simulated targets (targets O and P1-P6) before and after the higher-order RCMC operation.We can find that in Figure 7a,c, the energy of a single target is spread across multiple range gates due to higher-order RCM.In Figure 7b,d, this problem is solved and the energy is gathered in one range gate.
Figure 8 shows the Doppler FM rate deviation induced by the azimuth variation of the Doppler centroid according to Equation (41). Figure 8 illustrates that FM rate varies considerably when the Doppler centroid changes, and it has a great influence on the azimuth focusing.A variance of 500 Hz in the Doppler centroid can cause a FM rate deviation of 30 Hz/s.If the problem of spatial variation of the Doppler centroid is not settled properly, only the reference target can be well-focused, and the imaging quality of the rest of the targets suffer great degradation.Figures 9 and 10 show the first-order fitting of the Doppler centroid and the second-order fitting of FM rate, and their corresponding approximation errors, respectively.It can be seen that ∆φ err2 and ∆φ err3 are less than π/4.Thus, the first-order and second-order approximations of the Doppler centroid and FM rate are precise enough, respectively.The approximation errors are small enough that we can ignore their impacts on azimuth compression.Figure 11 shows the imaging results of the simulated targets in the scene.The results are interpolated by adding zeros in the frequency domain and transforming back to the time domain.All of the targets are well-focused by the proposed imaging algorithm.We can find that the aforementioned targets are gathered in the same range gates, and the azimuth positions of the targets are not changed by the processing because β is chosen to be 0.5 in this simulation.Figures 12 and 13 show the detailed imaging results of targets O, P2, P6, P7, P9, and P11.The results of the proposed method are compared with two state-of-the-art methods-i.e., the traditional NLCS method of Qiu et al. [34] and the fast backprojection (FBP) method.In the traditional azimuth NLCS method, the spatial variation of the Doppler FM rate is first obtained by local curve fitting, and the azimuth perturbation function is then derived by numerical integration.However, this ignores the azimuth variation of the Doppler centroid, and the FM rate deviation induced by it is not removed by the NLCS process.Thus, after the traditional NLCS process, the Doppler FM rates of the targets in the same range gate are still not equalized.Inspecting Figure 12g-l, it is clear that only the azimuth scene center is well-focused, and the azimuth compression performance of other targets suffer great degradation because the Doppler FM rates are still different from each other.In the upper two rows, the Doppler centroid and the FM rates are both equalized by the extended NLCS algorithm.Thus, when multiplied by the same azimuth compression factor, all of the targets get well focused.On the other hand, comparing the imaging results of the proposed method with that of the FBP method in Figure 13, it can be seen that the focusing performance of our method is close to the FBP method.However, the FBP method has more computational complexity than the proposed method, which is the drawback of time-domain imaging methods.The measured parameters peak sidelobe ratio (PSLR) and integrated sidelobe ratio (ISLR) of targets P2, P5, P6, P7, P9, and P11 are listed in Table 2.The PLSR and ISLR of the imaging results of the proposed are around −13 dB and −10 dB, which are nearly theoretical values.However, the parameters of the NLCS algorithm are worse compared to the proposed algorithm in the azimuth direction.

Real Data Processing Results
In order to further verify the effectiveness of the proposed method, real airborne bistatic SAR data processing results are given in Figure 14, and comparison between the traditional NLCS method of Qiu in Figure 15 et al. [34] is made in this subsection.The bistatic SAR system works in the X band with 75 MHz chirp bandwidth.The pulse repetition frequency is 500 Hz, and the pulse width is 20 µs.The velocity is 42 m/s, and the flying height for the transmitter and the receiver are 800 m and 600 m, respectively.

Discussion
Based on the results and analysis in Section 5, it can be concluded that the proposed algorithm can achieve accurate focusing for TV-BFSAR data.The main problems for TV-BFSAR imaging are 2-D spatial variant RCM and azimuth variant Doppler parameters (i.e., Doppler centroid and FM rate).For the first problem, the proposed RCMC method using keystone transform and higher-order RCM correction can remove the spatially-variant RCM of all the targets precisely compared with current methods, such as linear range walk correction (LRWC).On the other hand, the extended azimuth NLCS can equalize the azimuth variant Doppler centroid as well as FM rate, which obtains better azimuth focusing performance than the traditional NLCS method.
It should be noted that the ancillary data (such as Tx and Rx position, velocity, and attitude) should be accurate enough for the application of the method.These data can be either obtained by GPS and attitude instrument or by autofocus techniques [35][36][37].Here we assume that the ancillary data are derived from GPS and attitude instrument onboard the Tx and Rx and are accurate enough for the processing method.Moreover, if motion error exists for the transmitting and receiving platforms, motion compensation should be utilized to obtain satisfactory imaging results.Autofocus techniques can be integrated in the proposed method to estimate the ancillary parameters and further improve the focusing performance, which is our future work.

Conclusions
In this paper, an imaging algorithm for translational-variant bistatic forwarding SAR is proposed.The algorithm can be grouped into two parts: the range RCMC and azimuth extended NLCS.The range RCMC mainly deals with the problem of the spatial variation of first-order RCM and the residual higher-order RCM.The extended NLCS equalizes the Doppler centroid as well as FM rate.The range RCMC procedures are as follows.First, Keystone Transform is performed on the echo signal to remove the spatially-variant linear RCM.Second, the echo signal is compressed using a new range FM rate.Third, a correction factor is multiplied to the echo data to compensate for the higher-order RCM.Then, the azimuth dependence of the Doppler characters is calculated based on the relationship of cross-range coordinate x with azimuth position y and slant range R. The azimuth processing can also be divided into three parts.First, a phase factor is multiplied to the echo to move the Doppler spectra of the reference targets to the center of the azimuth frequency with respect to every range gate.Second, a fourth-order filtering is performed, which can introduce two degrees of freedom in the azimuth phase term.Third, an azimuth NLCS operation is carried out to equalize the azimuth-variant Doppler centroid and FM rate.At last, the targets in the same range gate are compressed by the same reference function.
Verified by the simulated and real TV-BFSAR data, this algorithm resolves the problem of the variation of both Doppler centroid and linear RCM, and thus has higher accuracy for all the targets in the scene.The residual RCM after keystone transform and higher-order RCMC is less than one resolution cell, which is accurate enough for TV-BFSAR imaging.The proposed method can achieve accurate focusing of TV-BFSAR data compared with tradition NLCS methods.The PSLR and ISLR of the simulated target are approximately −13 dB and −10 dB.Last but not least, the proposed method is more efficient than the backprojection method.

Figure 5 .
Figure 5. Target area used in the simulation.

Figure 6 .
Figure 6.Residual range cell migration (RCM) after KT, and its spatial variance.(a) Diagram of residual RCMs of O, P2, and P5 after KT; (b) Spatial variance of residual RCM between O and P2, O, and P5.

Figure 7 .Figure 8 .
Figure 7. Effect of the higher-order range cell migration correction (RCMC).(a) Seven targets before higher-order RCMC; (b) Seven targets after higher-order RCMC; (c) Detailed figure of the three targets in the 4097 range bin before higher-order RCMC; (d) Detailed figure of the three targets in the 4097 range bin after higher-order RCMC.

Figure 9 .Figure 10 .
Figure 9. Diagrams of the first-order fitting of the Doppler centroid and its approximation error.(a) Doppler centroid and its fitting curve; (b) Doppler centroid approximation error.

Figure 11 .
Figure 11.Imaging results of the seven simulated targets.

Figure 12 .Figure 13 Figure 13 .
Figure 12.Imaging results.(a)-(f) are the imaging results processed by the proposed algorithm; (g)-(l) are the imaging results processed by the traditional NLCS algorithm.

Figure 14 .
Figure 14.Real data processing results by the proposed method.

Figure 15 .
Figure 15.Real data processing results by traditional NLCS method.