Modeling and Precise Processing for Spaceborne Transmitter/Missile-Borne Receiver SAR Signals

: The spaceborne transmitter/missile-borne receiver (ST/MR) synthetic aperture radar (SAR) could provide several unique advantages, such as wide coverage, unrestricted geography, a small detection probability of the missile, and forward-looking imaging. However, it is also accompanied by problems in imaging, including geometric model establishment and focusing algorithm design. In this paper, an ST/MR SAR model is ﬁrst presented and then the ﬂight-path constraint, characterized by geometric conﬁgurations, is derived. Considering the impacts brought about by the maneuvers of the missile, a non-‘Stop-Go’ mathematical model is devised and it can avoid the large errors introduced by the acceleration, which is neglected in the traditional model. Finally, a two-dimensional (2-D) scaling algorithm is developed to focus the ST/MR data. Without introducing any extra operations, it can greatly remove the spatial variations of the range, azimuth, and cross-coupling phases simultaneously and entirely in the 2-D hybrid domain. Simulation results verify the effectiveness of the proposed model and focusing approach. transmitter and the receiver, the HRH would introduce significant phase errors for the satellite, missile, or aircraft platforms with a curved path. In [11,16,17], several accurate models, based on the Taylor series expansion, are established for the air-missile borne or missile BiSARs. Compared with the HRH, the presented models introduce accelerations into the range histories. However, for ST/MR SAR, the variation of the instantaneous slant range introduced by the continuous motion is no longer negligible since the conventional ‘stop - go’ approximation does not hold. In [18 – 20] , the effects of the ‘Stop - Go’ approximations of the monostatic SAR have been analyzed and the non- ‘Stop - Go’ models have been presented. However, by neglecting the maneuvers of the missile in the BiSAR system, the errors of these non- ‘Stop - Go’ models would greatly deteriorate the final image for the ST/MR SAR.


Introduction
The bistatic synthetic aperture radar (BiSAR) is characterized by different locations for the transmitter and receiver and hence offers a considerable capability, reliability, and flexibility in designing BiSAR missions. Compared with the monostatic synthetic aperture radar (SAR), BiSAR systems can achieve many benefits, like the exploitation of additional information contained in the bistatic reflectivity of targets, an increased radar cross section, and bistatic along-track interferometry [1,2]. Unlike the conventional BiSAR, installed on aircrafts [3,4], satellites [5,6], or satellite/aircraft hybrid platforms [7,8], the spaceborne transmitter/missile-borne receiver SAR (ST/MR SAR) transmits and receives the radar signals by a transmitter on the satellite and receiver on the missile, respectively, as shown in Figure 1. As a transition platform between the satellite and the aircraft corresponding to the altitude and the velocity, the missile with great maneuvers in the ST/MR SAR system is able to provide a wider observation area than that of the aircraft and greater flexibility than that of the satellite. It has broad application prospects in the military fields, namely, comprehensive detection and precision guidance, according to its strong survival ability, fast response, and wide detection range features [9].
Considering the functions and characteristics of the missile, it is better to improve the concealment, reduce the power consumption, and increase the detection range in the SAR system design [10]. Because the transmitter and receiver devices are both mounted on the missile platform, the monostatic In the traditional BiSAR system, the platforms are required to fly at constant velocities (magnitudes and directions), and the hyperbolic range history (HRH) hypotheses can be wellrepresented for both the transmitter and receiver motions during the synthetic aperture [1,4]. Considering the complex geometry and flight characteristics of both the transmitter and the receiver, the HRH would introduce significant phase errors for the satellite, missile, or aircraft platforms with a curved path. In [11,16,17], several accurate models, based on the Taylor series expansion, are established for the air-missile borne or missile BiSARs. Compared with the HRH, the presented models introduce accelerations into the range histories. However, for ST/MR SAR, the variation of the instantaneous slant range introduced by the continuous motion is no longer negligible since the conventional 'stop-go' approximation does not hold. In [18][19][20], the effects of the 'Stop-Go' approximations of the monostatic SAR have been analyzed and the non-'Stop-Go' models have been presented. However, by neglecting the maneuvers of the missile in the BiSAR system, the errors of these non-'Stop-Go' models would greatly deteriorate the final image for the ST/MR SAR.
Several studies on BiSAR focusing approaches have been published in recent years. Some popular monostatic SAR signal processing algorithms can be used for the BiSAR based on some In the traditional BiSAR system, the platforms are required to fly at constant velocities (magnitudes and directions), and the hyperbolic range history (HRH) hypotheses can be well-represented for both the transmitter and receiver motions during the synthetic aperture [1,4]. Considering the complex geometry and flight characteristics of both the transmitter and the receiver, the HRH would introduce significant phase errors for the satellite, missile, or aircraft platforms with a curved path. In [11,16,17], several accurate models, based on the Taylor series expansion, are established for the air-missile borne or missile BiSARs. Compared with the HRH, the presented models introduce accelerations into the range histories. However, for ST/MR SAR, the variation of the instantaneous slant range introduced by the continuous motion is no longer negligible since the conventional 'stop-go' approximation does not hold. In [18][19][20], the effects of the 'Stop-Go' approximations of the monostatic SAR have been analyzed and the non-'Stop-Go' models have been presented. However, by neglecting the maneuvers of the missile in the BiSAR system, the errors of these non-'Stop-Go' models would greatly deteriorate the final image for the ST/MR SAR. Several studies on BiSAR focusing approaches have been published in recent years. Some popular monostatic SAR signal processing algorithms can be used for the BiSAR based on some approximations and simplifications of the dual square root (DSR) [21,22]. The common shortcoming of these algorithms is their rigorous applying condition. For example, the range between the transmitter and receiver is always required to be relatively short. Neo etc. [23] proposed a way to get the formula of a two-dimensional (2-D) spectrum using the method of series reversion (MSR). Bamler and Boerner [24] and Bamler etc. [25] suggested the use of numerically computed transfer functions for bistatic SAR processing. In the case of the 2-D spectrum with MSR, the range Doppler algorithm (RDA) [26], the chirp scaling algorithms (CSA) [27], the frequency scaling algorithm (FSA) [28], the omega-K algorithm (OKA) [29], and their extended forms [30,31] have been accommodated. All the aforementioned algorithms, however, directly deal with the original data and focus on the analysis of the characteristics of the BiSAR data spectrum without taking the maneuvers into consideration. Considering the maneuvers of the platforms, [17] proposed a frequency domain algorithm (FDA) to eliminate the cross-coupling phases for the missile-borne BiSAR. However, neglecting the spatial variations would greatly defocus the image. Zhang etc. [32] proposed a nonlinear CSA (NCSA) for the maneuvering-platform BiSAR which could eliminate the spatial variations of the azimuth modulation phases, but could not do so for the range modulation and cross-coupling phases. Considering the acceleration of the missile, a generalized polar format algorithm (GPFA) [16] has been proposed to focus the air-missile borne BiSAR data. However, the wavefront curvature would greatly affect the depth of field and it should be compensated for by sub-aperture techniques [33,34]. Moreover, the BiSAR data with non-uniform curves, caused by the acceleration, can be focused by the back projection algorithm (BPA) [35][36][37]. However, in terms of computation burden, the BPA is not always a best choice compared with the frequency domain algorithms.
In this paper, the geometric model is established firstly, and the flight-path constraint and the mathematic model are analyzed and presented accordingly. This new mathematic model can well-present the motion of the satellite and the missile by taking several factors into consideration. Then, a 2-D scaling algorithm is proposed, which could greatly remove the range, azimuth, and cross-coupling spatially variant terms simultaneously in the 2-D hybrid domain by using two perturbation functions. Moreover, to avoid possible spectra distortion and aliasing, the distortion minimization is operated between these two perturbations.
This paper starts with a description of the geometric model by means of vector analysis (Section 2). Then, considering the complex geometric configurations, a 2-D scaling algorithm is proposed (Section 3). Section 4 presents the discussions of the proposed approach and Section 5 gives the numerical simulated results to validate the effectiveness of the presented mathematic model and imaging approach. In Section 6, the conclusion is drawn.

Geometric Model
Assume that the frequency modulation chirp pulse is used as the transmitted signal for the ST/MR SAR. Figure 2 shows the ST/MR SAR geometric model. In the coordinate of the imaged scene, point C is the central reference point (CRP) and point A is an arbitrary point on the ground scene. Points S and M are the positions of the transmitter and receiver at aperture center moment (ACM), respectively. v t and v r are respectively the transmitter and receiver velocity vectors and a r is the acceleration vector of the missile-borne receiver. In this paper, 't' and 'r' denote the transmitter and receiver, respectively. After reviewing the studies in the available literature [9,10,16,[38][39][40][41][42][43], we integrate the important results and then show that the flight path of spaceborne SAR can be equivalent to a curved path, which is mainly caused by the curved orbit and the rotating planetary surface and the flight path of the missile can be described by using the motion parameters, including the velocity, acceleration, and higher-order parameters, obtained from the information provided by the inertial navigation system (INS). According to the imaging geometry, the instantaneous slant range |r(η, A)| of the point A is expressed as: where r t (0, A) and r r (0, A) are the slant range vectors of the transmitter and receiver, respectively, to point A at ACM; η is the azimuth time; and |·| is the absolute operator. In this work, the symbols 'A' and 'C' denote the arbitrary target A and the reference target C, respectively. Particularly, the higher-order motion parameters of the satellite and the missile are not included in this work because they have no impact on the image qualities, which have been analyzed in Appendix A. One way to treat the complex range history is to expand it in a power series in azimuth time as: where µ n (A)(n = 0, 1, 2, 3, · · ·) are expanded coefficients of |r(η, A)|, which are expressed in Appendix A.
A are the slant range vectors of the transmitter and receiver, respectively, to point A at ACM;  is the azimuth time; and  is the absolute operator. In this work, the symbols 'A' and 'C' denote the arbitrary target A and the reference target C, respectively. Particularly, the higher-order motion parameters of the satellite and the missile are not included in this work because they have no impact on the image qualities, which have been analyzed in Appendix A. One way to treat the complex range history is to expand it in a power series in azimuth time as:

Flight-Path Constraint
A key factor in determining the imaging performance of the BiSAR is the spatial resolution, which we can also say is the acquisition geometry, including the positions and velocities of both the transmitter and the receiver. This is because it determines the spatial resolution [44]. In this part, the flight-path constraint of the ST/MR SAR characterized by arbitrary configurations is analyzed. The application of vector methods to physical problems most frequently takes the form of differential operations. According to the definition of the gradient operator  , the rates of changes of  0 and  1 in the directions of range and azimuth are equal to   0 and   1 , respectively. Accordingly, the ground space can be spanned by two unit vectors and , which are the projections of   0 and   1 on the ground, respectively. Vectors and are normal to the iso-range and iso-Doppler lines, respectively. Thus, the angle  between the iso-range and iso-Doppler lines is derived as: where  is the angle between the iso-range and iso-Doppler lines and   0 2 . The angle  in Equation (3) can be used to validate the rationality of the BiSAR geometry. The geometric model is more like that of the broadside monostatic SAR when  is close to  2 . Generally, ST/MR SAR with a large  has a better geometric structure than that of the ST/MR SAR with a small one.
Employing the simulation parameters in Table 1

Flight-Path Constraint
A key factor in determining the imaging performance of the BiSAR is the spatial resolution, which we can also say is the acquisition geometry, including the positions and velocities of both the transmitter and the receiver. This is because it determines the spatial resolution [44]. In this part, the flight-path constraint of the ST/MR SAR characterized by arbitrary configurations is analyzed. The application of vector methods to physical problems most frequently takes the form of differential operations. According to the definition of the gradient operator ∇, the rates of changes of µ 0 and µ 1 in the directions of range and azimuth are equal to ∇µ 0 and ∇µ 1 , respectively. Accordingly, the ground space can be spanned by two unit vectors and l, which are the projections of ∇µ 0 and ∇µ 1 on the ground, respectively. Vectors and l are normal to the iso-range and iso-Doppler lines, respectively. Thus, the angle α between the iso-range and iso-Doppler lines is derived as: where α is the angle between the iso-range and iso-Doppler lines and 0 ≤ α ≤ π/2. The angle α in Equation (3) can be used to validate the rationality of the BiSAR geometry. The geometric model is more like that of the broadside monostatic SAR when α is close to π/2. Generally, ST/MR SAR with a large α has a better geometric structure than that of the ST/MR SAR with a small one. Employing the simulation parameters in Table 1 and the simulation scene shown in Figure 3a, Figure 3b shows the iso-range/iso-Doppler lines by using the parameters of Receiver 1 in Table 1 and Figure 3c is the simulation result by the parameters of Receiver 2. In these figures, the angle α between iso-range and iso-Doppler lines of each point on the ground scene can be derived by using Equation (3). Clearly, in the blue area in these two figures, the angles between iso-range and iso-Doppler lines are very small, even close to zero, which means that the cross-couplings between the range and azimuth are extremely serious. Thus, it is difficult to focus the echo data received from these areas.  Figure 3b shows the iso-range/iso-Doppler lines by using the parameters of Receiver 1 in Table 1 and Figure 3c is the simulation result by the parameters of Receiver 2. In these figures, the angle α between iso-range and iso-Doppler lines of each point on the ground scene can be derived by using Equation (3). Clearly, in the blue area in these two figures, the angles between iso-range and iso-Doppler lines are very small, even close to zero, which means that the cross-couplings between the range and azimuth are extremely serious. Thus, it is difficult to focus the echo data received from these areas.

Flight-Path Constraint
Traditional methods considering the 'Stop-Go' assumption have been widely studied and can be well-applied for the SAR on the high-speed platform [18][19][20]. However, in the case of large maneuvers of the missile, the impacts brought about by the accelerations are beyond consideration. In this part, a new mathematical model is established to describe the ST/MR SAR with a high accuracy.

Flight-Path Constraint
Traditional methods considering the 'Stop-Go' assumption have been widely studied and can be well-applied for the SAR on the high-speed platform [18][19][20]. However, in the case of large maneuvers of the missile, the impacts brought about by the accelerations are beyond consideration. In this part, a new mathematical model is established to describe the ST/MR SAR with a high accuracy.
It is generally known that the exact solution of the time delay τ brought about by the 'Stop-Go' assumption is difficult to obtain because of the higher-order transcendental function. To solve this problem, an iterative method is presented to explore an efficient solution. Figure 4 shows the transmitted and received instants, taking the platform movement into account. The radar signal is transmitted by the radar on the satellite at slow time T start and received by the missile at the position χ at T end , as shown in Figure 4, and T start is a known initial time for the transmission. Thus, the propagation time τ can be derived as τ = T end − T start . The basic procedure of the iterative method is listed in Algorithm 1, where c is the speed of light, and r t (0, A), r r (0, A), v t , v r , and a r have been defined in Section 2.1. The iterative range history can be regarded as a satisfactory result if the iterative error is less than the error tolerance ∆. The iterative algorithm is applicable for different SAR with maneuvers.

Flight-Path Constraint
Traditional methods considering the 'Stop-Go' assumption have been widely studied and can be well-applied for the SAR on the high-speed platform [18][19][20]. However, in the case of large maneuvers of the missile, the impacts brought about by the accelerations are beyond consideration. In this part, a new mathematical model is established to describe the ST/MR SAR with a high accuracy.  1: Fixed input parameters: Clearly, the complicated procedure would lead to a complex result and make the focusing algorithm difficult to design. To avoid the cumbersome procedure, a 'fixed-value' iterative method is developed. By analyzing the effects brought about by different iterative numbers, we find that the delay time errors tend to the extremely small fixed-values after one time iteration. Thus, using one-time iteration, the new range history can be derived as: − a r η 2 /2 /c and κ n (A)(n = 0, 1, 2, 3, · · ·) are the expansion coefficients which can be found in Appendix B. Employing the simulation parameters in Table 2 and the scene in Figure 3a, Figure 5 shows the comparative results of the residual phase errors. It is noted that the proposed model can be well-applied for the ST/MR SAR, whereas the traditional method cannot because the acceleration is not considered.

Imaging Approach
The geometric model is established, and the mathematic model of the ST/MR SAR is derived in Section 2. Based on the accurate model, we concentrate on the corresponding imaging algorithm in this section. A 2-D scaling algorithm is proposed to focus the ST/MR SAR data. The flowchart is shown in Figure 6. Because of the complicated mathematic model, we design the imaging approach performed in the 2-D hybrid domain. Two short remarks will be helpful to understanding the 2-D scaling algorithm:  Cross-couplings. Traditionally, the cross-coupling phases in the 2-D spectrum derived by

Imaging Approach
The geometric model is established, and the mathematic model of the ST/MR SAR is derived in Section 2. Based on the accurate model, we concentrate on the corresponding imaging algorithm in this section. A 2-D scaling algorithm is proposed to focus the ST/MR SAR data. The flowchart is shown in Figure 6. Because of the complicated mathematic model, we design the imaging approach performed in the 2-D hybrid domain. Two short remarks will be helpful to understanding the 2-D scaling algorithm: Cross-couplings. Traditionally, the cross-coupling phases in the 2-D spectrum derived by using the principle of stationary point (POSP) and MSR are divided into two parts, i.e., the range-and azimuth-dependent ones, and they are eliminated separately, which would introduce errors in the final image. The 2-D scaling algorithm in this work can avoid this division and make the processing more accurate than that of the traditional methods. Spatial variations. Generally, the range, azimuth, and cross-coupling spatial variations in the 2-D spectrum are difficult to be removed simultaneously and entirely because of the complex phase expressions. In the 2-D scaling algorithm, two perturbation functions performed in the hybrid domains can eliminate the range, azimuth, and cross-coupling spatially variant terms simultaneously and entirely, which can avoid the defect of the traditional methods that the range and cross-coupling spatial variations cannot be removed.

Imaging Approach
The geometric model is established, and the mathematic model of the ST/MR SAR is derived in Section 2. Based on the accurate model, we concentrate on the corresponding imaging algorithm in this section. A 2-D scaling algorithm is proposed to focus the ST/MR SAR data. The flowchart is shown in Figure 6. Because of the complicated mathematic model, we design the imaging approach performed in the 2-D hybrid domain. Two short remarks will be helpful to understanding the 2-D scaling algorithm:  Cross-couplings. Traditionally, the cross-coupling phases in the 2-D spectrum derived by using the principle of stationary point (POSP) and MSR are divided into two parts, i.e., the range-and azimuth-dependent ones, and they are eliminated separately, which would introduce errors in the final image. The 2-D scaling algorithm in this work can avoid this division and make the processing more accurate than that of the traditional methods.  Spatial variations. Generally, the range, azimuth, and cross-coupling spatial variations in the 2-D spectrum are difficult to be removed simultaneously and entirely because of the complex phase expressions. In the 2-D scaling algorithm, two perturbation functions performed in the hybrid domains can eliminate the range, azimuth, and cross-coupling spatially variant terms simultaneously and entirely, which can avoid the defect of the traditional methods that the range and cross-coupling spatial variations cannot be removed.

2-D Scaling in Range Frequency/Azimuth Time Domain
Assume that a linear frequency modulation (LFM) signal is transmitted with the frequency modulation (FM) rate being γ. After range compression, the received echo signal in the range frequency/azimuth time domain can be expressed as where A r (·) and A a (·) are the range and azimuth windows, respectively; and f c and f r are the carrier and range frequencies, respectively. Note that the traditional methods are difficult to apply for the ST/MR SAR case directly because of the complex phase in Equation 5. In the 2-D scaling algorithm, a perturbation function, performed in the 2-D hybrid domain (hybrid domain between range frequency f r and azimuth time η), is firstly designed to weaken the spatially variant phases and also provides many more DOFs for the echo signal, i.e., where the coefficients ς 2 , ς 3 , and ς 4 are all the scaling factors and can be derived at the next step. Multiplying Equation (6) with Equation (5) and ignoring the unimportant amplitudes, one can obtain that Note that the spatial variations of the high-order terms corresponding to κ n (A)(n = 5, 6, · · ·) in Equation (7) are generally very small and can be ignored, which indicates that the second phase term can be compensated for by using the bulk compensation function corresponding to the reference target. Then, the echo signal can be expressed as

Distortion Minimization
After perturbation processing, the echo signal is transformed into the 2-D spectrum. Of particular note is the spectral distortion and aliasing after azimuth Fourier transform (FT). Three main uncertainty factors should be faced and they could greatly deteriorate the final image if we do not pay attention to them: Generally, by steering the beam to different fixed centers, the spatial resolution can be controlled flexibly during the data acquisition interval. Consequently, the corresponding azimuth signal bandwidth might span over several pulse repetition frequency (PRF) intervals [45,46]. The expression of the azimuth bandwidth differs from the usual representation based on the linear path because the rate of the Doppler frequency modulation (FM) κ C 2 is greatly affected by the motion parameter vectors a r . Thus, possible azimuth bandwidth widening should be considered in this work [47,48]. After multiplying Equation (6) with echoes, the rate of the Doppler FM becomes κ C 2 + ζ, which means that the azimuth bandwidth may display a big change. Thus, the possible spectral distortion should be eliminated before the azimuth FT of the signal.
The spectral aliasing elimination method has been discussed in [48] in detail and one can use it for efficient data processing. The processes include the deramping function H 1 , compensation function H 2 , and quadratic function H 3 [48]. The only difference is that the original Doppler FM should be substituted by κ C 2 + ζ so that the spectral distortion and aliasing partly introduced by the perturbation function can be eliminated entirely. Particularly, the spectral aliasing elimination method is an extension of the traditional two-step approach in [49,50], whose essence is the time-frequency rotation. Thus, two azimuth FTs are used. The first one is used as a filter to change the sampling interval and the second one is a normal transformation.
After distortion minimization, the 2-D spectrum of the echo signal is given by where f η is the new azimuth frequency corresponding to the new azimuth time η after distortion minimization, i.e., f η = (κ 2 (C) + ζ)·η/λ and f η = (κ 2 (C) + ζ)·η /λ, where f η is the azimuth frequency before distortion minimization. ξ i (A)(i = 0, 1, 2, 3, 4) represents Doppler parameters of the range equation corresponding to A, all of which are spatially variant terms and their expressions can be found in Appendix C. It is important to note that the expressions of the 2-D spectrum before and after distortion minimization are the same, except for f η and f η .

2-D Scaling in Range Frequency/Azimuth Frequency Domain
The coefficients ξ 0 (A) and ξ 1 (A) are the range and Doppler positions of A, respectively. To obtain the well-focused image, the higher-order terms corresponding to ξ 2 (A), ξ 3 (A), and ξ 4 (A) should be removed entirely. In this sub-section, another perturbation function, performed in the 2-D hybrid domain (hybrid domain between range frequency f r and azimuth frequency f η ), is presented, i.e., where ω 2 , ω 3 , and ω 4 are all the scaling factors similar to those of ς 2 , ς 3 , and ς 4 in Equation (6). After the second perturbation processing by multiplying Equation (10) with Equation (9), the phase terms of the echo signal are derived as Clearly, the range, azimuth, and cross-coupling spatial variations of A are all included in the third term in Equation (11) Clearly, three equations with six unknown scaling factors make it difficult to derive the solutions for Equation (12). To solve this problem, we firstly rewrite Equation 12 as , and M 0 = 16ω 4 . Then, β 2 (A), β 3 (A), and β 4 (A) are respectively expanded in the power series with respect to κ 1 (A), i.e., (14) where β 2 (C),  . Note that both Equations (13) and (14) are used for the expressions of β 2 (A), β 3 (A), and β 4 (A), but with different formulae. Making use of the similarities and differences between Equations (13) and (14), we can obtain the following six equations with six unknowns ξ i , ω i (i = 2, 3, 4), i.e., (15) Solving this equation set in Equation (15), we can derive its solutions as where , and P 3 = . After the 2-D scaling processing in the hybrid domain by using Equations (16) and (17), the echo signal can be expressed as Clearly, the phase in Equation 18 is a 2-D linear one with respect to the range frequency f r and azimuth frequency f η . By performing a 2-D inverse FT (IFT) on Equation (18), one can obtain a focused image for ST/MR SAR data. Particularly, without introducing any extra operations, the range, azimuth, and cross-coupling spatially variant terms in Equation (5) can be removed simultaneously and entirely by using the 2-D scaling algorithm, which indicates that the proposed approach could be easy and efficient to implement for the ST/MR SAR data. Next, some discussions are given to facilitate a better understanding of the proposed method.

Simulation Results
In this section, some simulation results from the two experiments are shown to validate the proposed method and analysis.

Case I
In this experiment, the effectiveness of the mathematical model presented in Equation (4) is analyzed and verified. The traditional 'Stop-Go' and the non-'Stop-Go' models in [20] are compared with the proposed 2-D scaling algorithm. The simulation parameters are listed in Table 2. The ST/MR SAR data are simulated with one reference target PT0 being arranged on the ground scene.  with the proposed 2-D scaling algorithm. The simulation parameters are listed in Table 2. The ST/MR SAR data are simulated with one reference target PT0 being arranged on the ground scene. Figure 7a,b,c show the focused results of the ground scene by the proposed model, the traditional non-'Stop-Go' model, and the traditional 'Stop-Go' model, respectively. Clearly, considering the negative impacts brought about by the missile's accelerations, the focused result of the reference target PT0 in Figure 7a is visually well-focused, with the position being the scene center. However, using the traditional non-'Stop-Go' and 'Stop-Go' models, the imaging results exhibit great deterioration in the azimuth directions with wrong target positions marked on the ground scenes, as shown in Figure 7b,c.

Case II
In order to verify the effectiveness of the proposed approach, some simulated results are given in this case. The simulated parameters are listed in Table 3. A 3 × 3 dot-matrix is arranged on the simulated ground scene, with center target PT2 being the reference target. The geometry of the scene is presented in Figure 9, with the scene size being 2 km × 2 km.    Table 3. A 3 × 3 dot-matrix is arranged on the simulated ground scene, with center target PT2 being the reference target. The geometry of the scene is presented in Figure 9, with the scene size being 2km × 2km.  By performing the proposed algorithm, simulated results of targets PT1, PT2, and PT3 are shown in Figure 10a. It can be seen that the main lobe and the side lobes are generally separated, and the center point PT2 and the edge points PT1 and PT3 are all well-focused. For comparisons, the 2-D impulse responses of these three targets obtained by the NCSA [32] and the GPFA [16] are respectively shown in Figure 10b,c. Experimental results indicate that the center points on the ground scene are focused well, whereas the ones on the edges remain highly defocused in the azimuth direction. For the NCSA, the range and cross-coupling spatial variations are not considered and they would lead to great deterioration in the final image. Moreover, the residual phases after NCSA, eliminated by the deramping function, would also introduce errors in the imaging result. The GPFA can be used for the ST/MR SAR; however, its depth of field is greatly affected by the wavefront curvature.

Case II
The profiles of azimuth spread functions of targets PT1 and PT3 are presented in Figure 11. To quantify the precision of the proposed approach, the impulse-response width (IRW), peak sidelobe ratio (PSLR), and integrated sidelobe ratio (ISLR) are also used as criteria. The results are listed in Table 4   By performing the proposed algorithm, simulated results of targets PT1, PT2, and PT3 are shown in Figure 10a. It can be seen that the main lobe and the side lobes are generally separated, and the center point PT2 and the edge points PT1 and PT3 are all well-focused. For comparisons, the 2-D impulse responses of these three targets obtained by the NCSA [32] and the GPFA [16] are respectively shown in Figure 10b,c. Experimental results indicate that the center points on the ground scene are focused well, whereas the ones on the edges remain highly defocused in the azimuth direction. For the NCSA, the range and cross-coupling spatial variations are not considered and they would lead to great deterioration in the final image. Moreover, the residual phases after NCSA, eliminated by the deramping function, would also introduce errors in the imaging result. The GPFA can be used for the ST/MR SAR; however, its depth of field is greatly affected by the wavefront curvature.
The profiles of azimuth spread functions of targets PT1 and PT3 are presented in Figure 11. To quantify the precision of the proposed approach, the impulse-response width (IRW), peak sidelobe ratio (PSLR), and integrated sidelobe ratio (ISLR) are also used as criteria. The results are listed in Table 4. The proposed method is noted to be generally acceptable. Moreover, to further elevate the performance of the 2-D scaling algorithm, the position errors are computed in this part. The deviations of point targets PT1, PT2, and PT3 from their ideal values are 0.912 m, 0.475 m, and 1.016 m, respectively. Clearly, the deviations were generally within or around one IRW. Consequently, with the incorporation of 2-D perturbations and distortion minimization, promising results are obtained for the ST/MR SAR.

2-D Range Compression
Antenna motion during the transmission and reception of a pulse causes pulse distortion and is analogous to the familiar Doppler frequency shift associated with sound waves [1]. Typically, the phase effect of motion during a pulse is small, but this effect can be become significant as pulse lengths increase and motion becomes faster. Employing the simulation parameters in Table 2, the phase error brought about by the motion during transmission and reception of a pulse is about 3.73π and the extra phase error brought about by the acceleration term is about 0.007π, which indicates that the effect of the motion during transmission and reception should be considered, whereas the impact by the acceleration can be ignored. Actually, the effect has been already discussed in [18][19][20]. Based on the analytical derivations in [19,20], the effect can be compensated by Multiplying Equation (19) by the range compression function in the range frequency domain, the range compression together with the effect compensation (in this work, we call it the 2-D range compression operation) can be expressed as

Computational Load
The computational load is analyzed in this section. In the 2-D scaling algorithm, two range FTs, three azimuth FTs, and four multiplications are utilized. The computational load can be derived as MN log 2 N + 3N M log 2 M/2 + 4N M, where N and M are the numbers of range and azimuth samples, respectively. In comparison, for the BPA (e.g., [35]), for each pixel in the final image, a signal vector with a length M is extracted from range compressed data, multiplied with a phase function, and summed; thus, the total computational load can be expressed as O B = M 2 N. In a standard CSA (e.g., [51]), two range FTs, one azimuth FT, and three multiplications are used. Therefore, the computational load can be expressed as MN log 2 N + N M log 2 M + 3N M.
Clearly, the BPA has the highest computational load. The load of the 2-D scaling algorithm is slightly heavier than that of the CSA. Ratios at five sample sizes of both range and azimuth are computed to quantify the comparison (Table 5).

Conclusions
In this paper, a new BiSAR imaging configuration is presented, with the satellite being the transmitter and missile being the receiver. This system could provide many advantages, including wide coverage, unrestricted geography, low detection probability of the missile, and forward-looking imaging. A geometric model for the ST/MR SAR is first established by using motion parameter vectors and the flight-path constraint is derived, and the corresponding focusing approach is then proposed. In this approach, two perturbation functions are designed in the 2-D hybrid domain to eliminate the range, azimuth, and cross-coupling spatially variant terms in the ST/MR data and distortion minimization is applied to avoid the possible spectral distortion and aliasing in the 2-D spectrum. Furthermore, discussions including the 2-D range compression and the computational load are studied. Validity and applicability are studied through theoretical analysis and numerical experiments.
Author Contributions: S.T. conceived the main idea; P.G. and L.Z. conceived and designed the experiments; S.T. and P.G. analyzed the data and wrote the paper; S.T. revised this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this Appendix, the coefficients of the Taylor series are listed as Table A1. Coefficients of Taylor series expansion.

Coefficient
Transmitter Coefficients Receiver σ 0 r t (0, A), r t (0, A) ν 0 r r (0, A), r r (0, A) where σ i (i = 0, 1, 2, 3, 4) and ν i (i = 0, 1, 2, 3, 4) in Table A1 are the Taylor series expansion coefficients of the range histories corresponding to the transmitter and the receiver, respectively. , is the inner operator. Thus, the coefficients µ i (i = 0, 1, 2, 3, 4) are derived as µ i = σ i + ν i . Considering the high-order accelerations, the range history of A can be expressed as r high (η, A) = |r t (η, A)| + |r r (η, A)| = r t (0, A) − v t η − a t η 2 /2 − b t η 3 /6 − c t η 4 /24 + r r (0, A) − v r η − a r η 2 /2 − b r η 3 /6 − c r η 4 /24 , where a t , b t , and c t are the second-, third-, and fourth-order accelerations of the transmitter, respectively. b r and c r are the third-and fourth-order accelerations of the receiver, respectively. Employing the simulated parameters in Table 2 and Equation (A1), Figure A1 shows the residual errors brought about by the high-order motion parameters. The unit of the contour maps is π. Clearly, the bulk phase errors are all larger than π/4, whereas the spatially variant ones are all far less than π/4, which indicates that the high-order motion parameters can be compensated by r high (η, C) − |r(η, C)| without considering the spatial variations and they have no impact on the focusing algorithm design.

Appendix B
By using the Taylor series expansion in Equation (2), the time delay τ can be expressed as