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

Abstract

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 first presented and then the flight-path constraint, characterized by geometric configurations, 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.

1. 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 missile-borne SAR will inevitably bring some limitations, such as the large power antenna and high probability of being detected. Moreover, the monostatic SAR is incapable of forward-looking imaging, which seriously limits its applicability in the military field because information in the movement direction cannot be obtained directly. Compared with the monostatic missile-borne SAR, the ST/MR SAR has many advantages and they are:

(1)

The large-maneuver and receiving-only features of the missile-borne receiver reduce the detection probability. The corresponding anti-interference and anti-interception are greatly improved in the ST/MR SAR system.

(2)

The radar antenna, installed on the missile-borne platform, is used to receive the echo signals without a transmitting function, which indicates that a large power device can be avoided in the missile system to greatly save space and costs.

(3)

Spaceborne transmitter has the characteristics of wide coverage, unrestricted geography, net flexibility, and a high revisiting-frequency, which could provide some degrees of freedom (DOFs) in target area selection, flight trajectories, and number of missile-borne receivers.

(4)

By setting the transmitter and receiver on different platforms with a special style, the non-complete coupling areas between iso-range and iso-Doppler lines could be formed in front of missile, which means that the missile has the capability of forwarding-looking imaging.

Note that these advantages for the ST/MR SAR also hold for other BiSAR systems, such as the spaceborne/airborne SAR and the airborne BiSAR, which are of special significance for the missiles mainly applied in the military fields. Generally, these common features of the BiSAR can avoid disadvantages inherent in the monostatic missile-borne SAR and make the ST/MR SAR a flexible and effective tool for information retrieval [11,12,13,14,15]. However, there are also a number of problems brought about by the ST/MR SAR, such as the complex geometric model, high cross-couplings, and serious spatial variations. Thus, further studies are still necessary for ST/MR SAR.

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.

2. Modeling

2.1. 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. ${\mathit{v}}_t$ and ${\mathit{v}}_r$ are respectively the transmitter and receiver velocity vectors and ${\mathit{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 $\left|\mathit{r}\left(\eta ,A\right)\right|$ of the point A is expressed as:

$$\begin{array}{ll}\left|\mathit{r}\left(\eta ,A\right)\right|& =\left|{\mathit{r}}_t\left(\eta ,A\right)\right|+\left|{\mathit{r}}_r\left(\eta ,A\right)\right|\\ & =\left|{\mathit{r}}_t\left(0,A\right)-{\mathit{v}}_t\eta \right|+\left|{\mathit{r}}_r\left(0,A\right)-{\mathit{v}}_r\eta -{\mathit{a}}_r{\eta }^2/2\right|\end{array}$$

(1)

where ${\mathit{r}}_t\left(0,A\right)$ and ${\mathit{r}}_r\left(0,A\right)$ are the slant range vectors of the transmitter and receiver, respectively, to point A at ACM; $\eta$ 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:

$$\left|\mathit{r}\left(\eta ,A\right)\right|={\displaystyle \sum _0\frac{1}{n!}{\mu }_n\left(A\right)·{\eta }^n}$$

(2)

where ${\mu }_n\left(A\right)\left(n=0,1,2,3,\cdots \right)$ are expanded coefficients of $\left|\mathit{r}\left(\eta ,A\right)\right|$, which are expressed in Appendix A.

2.2. 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 $\nabla$, the rates of changes of ${\mu }_0$ and ${\mu }_1$ in the directions of range and azimuth are equal to $\nabla {\mu }_0$ and $\nabla {\mu }_1$, respectively. Accordingly, the ground space can be spanned by two unit vectors $\hslash$ and $ƛ$, which are the projections of $\nabla {\mu }_0$ and $\nabla {\mu }_1$ on the ground, respectively. Vectors $\hslash$ and $ƛ$ are normal to the iso-range and iso-Doppler lines, respectively. Thus, the angle $\alpha$ between the iso-range and iso-Doppler lines is derived as:

$$\alpha ={\cos }^{-1}\left\{\left|\langle \hslash ,ƛ\rangle /\left(|\hslash |·\left|ƛ\right|\right)\right|\right\}$$

(3)

where $\alpha$ is the angle between the iso-range and iso-Doppler lines and $0\le \alpha \le \pi /2$. The angle $\alpha$ 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 $\alpha$ is close to $\pi /2$. Generally, ST/MR SAR with a large $\alpha$ has a better geometric structure than that of the ST/MR SAR with a small one. Employing the simulation parameters in

Table 1. Parameter settings.

	Platform	Transmitter	Receiver 1	Receiver 2
Parameter
Radar Position at ACM		(0, 1000, 500) km	(12.5, −4, 25) km	(−6, −4, 30) km
Velocity Vector		(0, 7600, 0) m/s	(0, 1000, 0) m/s	(−100, −1200, 800) m/s
Acceleration Vector		/	(10, −30, −50) m/s2	(50, −60, −35) m/s2
Elevation Angle		84.1°	45.8°	32.3°

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 $\alpha$ 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.

2.3. 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 $\tau$ 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 $\chi$ at $T_{end}$, as shown in Figure 4, and $T_{start}$ is a known initial time for the transmission. Thus, the propagation time $\tau$ can be derived as $\tau =T_{end}-T_{start}$. The basic procedure of the iterative method is listed in Algorithm 1, where $c$ is the speed of light, and ${\mathit{r}}_t\left(0,A\right)$, ${\mathit{r}}_r\left(0,A\right)$, ${\mathit{v}}_t$, ${\mathit{v}}_r$, and ${\mathit{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 $\mathsf{\Delta }$. The iterative algorithm is applicable for different SAR with maneuvers.

Algorithm 1. Basic procedure.

1: Fixed input parameters: ${\mathit{r}}_t\left(0,A\right)$, ${\mathit{r}}_r\left(0,A\right)$, ${\mathit{v}}_t$, ${\mathit{v}}_r$, ${\mathit{a}}_r$, $T_{start}$, $c$ 2: Initialization (i=0): $T_{end\left(0\right)}=T_{start}+\left|\mathit{r}\left(T_{start},A\right)\right|/c$ 3: Iterations (i≥1): 4: Calculate the delay: ${\tau }_{i-1}=T_{end\left(i-1\right)}-T_{start}$ 5: Calculate the slant range of the receiver: $\left|{\mathit{r}}_r\left(T_{end\left(i-1\right)},A\right)\right|=\left|{\mathit{r}}_r\left(0,A\right)-{\mathit{v}}_r{T}_{end\left(i-1\right)}-{\mathit{a}}_r{T}_{end\left(i-1\right)}^2/2\right|$ 6: If $\left|{\mathit{r}}_t\left(T_{start},A\right)\right|+\left|{\mathit{r}}_r\left(T_{end\left(i-1\right)},A\right)\right|-{\tau }_{i-1}·c\le \mathsf{\Delta }$, return $\tau ={\tau }_{i-1}$ and stop. 7: Calculate the time of the receiver: $T_{end\left(i\right)}=T_{start}+\left\{\left|{\mathit{r}}_t\left(T_{start},A\right)\right|+\left|{\mathit{r}}_r\left(T_{end\left(i-1\right)},A\right)\right|\right\}/c$ 8: Set $i=i+1$.

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:

$$\begin{array}{ll}\left|{\mathit{r}}_{new}\left(\eta ,A\right)\right|& =\left|{\mathit{r}}_t\left(0,A\right)-{\mathit{v}}_t\eta \right|+\left|{\mathit{r}}_r\left(0,A\right)-{\mathit{v}}_r\left(\eta +\tau \right)-{\mathit{a}}_r{\left(\eta +\tau \right)}^2/2\right|\\ & ={\displaystyle \sum _0\frac{1}{n!}{\kappa }_n\left(A\right)·{\eta }^n}\end{array}$$

(4)

where $\tau =\left|\mathit{r}\left(\eta ,A\right)\right|/c=\left(\left|{\mathit{r}}_t\left(0,A\right)-{\mathit{v}}_t\eta \right|+\left|{\mathit{r}}_r\left(0,A\right)-{\mathit{v}}_r\eta -{\mathit{a}}_r{\eta }^2/2\right|\right)/c$ and ${\kappa }_n\left(A\right)\left(n=0,1,2,3,\cdots \right)$ are the expansion coefficients which can be found in Appendix B. Employing the simulation parameters in

Table 2. Parameter settings for Case I.

	Platform	Transmitter	Receiver
Parameter
Radar Position at ACM		(0, 0, 754) km	(132, 0, 15) km
Velocity Vector		(0, 7600, 0) m/s	(300, 1200, −200) m/s
Acceleration Vector		/	(30, 54, −26) m/s2
Elevation Angle		74.1°	44.6°
Carrier Frequency		9.65 GHz
Pulse Width		20 μs
Pulse Bandwidth		200 MHz
Sampling Frequency		300 MHz
Azimuth Resolution		0.823 m

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.

3. 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.

3.1. 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 $\gamma$. After range compression, the received echo signal in the range frequency/azimuth time domain can be expressed as

$$\begin{array}{ll}{S}_0\left(f_r,\eta \right)& =A_r\left(f_r\right)A_a\left(\eta \right)\exp \left\{-j\frac{2\pi \left(f_c+f_r\right)}{c}\left|{\mathit{r}}_{new}\left(\eta ,A\right)\right|\right\}\\ & =A_r\left(f_r\right)A_a\left(\eta \right)\exp \left\{-j\frac{2\pi \left(f_c+f_r\right)}{c}{\displaystyle \sum _0{\kappa }_n\left(A\right){\eta }^n}\right\}\end{array}$$

(5)

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.,

$$H_{1st}\left(f_r,\eta \right)=\exp \left[-j\frac{2\pi \left(f_r+f_c\right)}{c}{\displaystyle \sum _{i=2}^4{\varsigma }_i{\eta }^i}\right]$$

(6)

where the coefficients ${\varsigma }_2$, ${\varsigma }_3$, and ${\varsigma }_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

$$\begin{array}{ll}{S}_1\left(f_r,\eta \right)& =\exp \left\{-j\frac{2\pi \left(f_c+f_r\right)}{c}\left[{\displaystyle \sum _{i=2}^4{\varsigma }_i{\eta }^i}+{\displaystyle \sum _0^4{\kappa }_n\left(A\right){\eta }^n}\right]\right\}\\ & ·\exp \left\{-j\frac{2\pi \left(f_c+f_r\right)}{c}{\displaystyle \sum _5{\kappa }_n\left(A\right){\eta }^n}\right\}\end{array}$$

(7)

Note that the spatial variations of the high-order terms corresponding to ${\kappa }_n\left(A\right)\left(n=5,6,\cdots \right)$ 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

$$S_2\left(f_r,\eta \right)=\exp \left\{-j\frac{2\pi \left(f_c+f_r\right)}{c}\left[{\displaystyle \sum _{i=2}^4{\varsigma }_i{\eta }^i}+{\displaystyle \sum _0^4{\kappa }_n\left(A\right){\eta }^n}\right]\right\}$$

(8)

3.2. 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) ${\kappa }_2^C$ is greatly affected by the motion parameter vectors ${\mathit{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 ${\kappa }_2^C+\zeta$, 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 ${\kappa }_2^C+\zeta$ 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

$$S_3\left(f_r,f_{\eta }^{\prime }\right)=\exp \left\{j2\pi ·{\displaystyle \sum _{i=0}^4{\xi }_i\left(A\right)·{f_{\eta }^{\prime }}^{\left(i\right)}·{\left[\frac{\left(f_r+f_c\right)}{c}\right]}^{\left(1-i\right)}}\right\}$$

(9)

where $f_{\eta }^{\prime }$ is the new azimuth frequency corresponding to the new azimuth time ${\eta }^{\prime }$ after distortion minimization, i.e., $f_{\eta }^{\prime }=\left({\kappa }_2\left(C\right)+\zeta \right)·\eta /\lambda$ and $f_{\eta }=\left({\kappa }_2\left(C\right)+\zeta \right)·{\eta }^{\prime }/\lambda$, where $f_{\eta }$ is the azimuth frequency before distortion minimization. ${\xi }_i\left(A\right)\left(i=0,1,2,3,4\right)$ 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_{\eta }^{\prime }$ and $f_{\eta }$.

3.3. 2-D Scaling in Range Frequency/Azimuth Frequency Domain

The coefficients ${\xi }_0\left(A\right)$ and ${\xi }_1\left(A\right)$ are the range and Doppler positions of A, respectively. To obtain the well-focused image, the higher-order terms corresponding to ${\xi }_2\left(A\right)$, ${\xi }_3\left(A\right)$, and ${\xi }_4\left(A\right)$ 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_{\eta }^{\prime }$), is presented, i.e.,

$$H_{2nd}\left(f_r,f_{\eta }^{\prime }\right)=\exp \left\{-j2\pi ·{\displaystyle \sum _{i=2}^4{\omega }_i^{}·{f_{\eta }^{\prime }}^{\left(i\right)}·{\left[\frac{\left(f_r+f_c\right)}{c}\right]}^{\left(1-i\right)}}\right\}$$

(10)

where ${\omega }_2$, ${\omega }_3$, and ${\omega }_4$ are all the scaling factors similar to those of ${\varsigma }_2$, ${\varsigma }_3$, and ${\varsigma }_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

$${\phi }_s\left(f_r,f_{\eta }^{\prime }\right)=2\pi {\xi }_0\left(A\right)\frac{\left(f_r+f_c\right)}{c}+2\pi {\xi }_1\left(A\right)f_{\eta }^{\prime }+2\pi ·{\displaystyle \sum _{i=2}^4\left({\xi }_i\left(A\right)-{\omega }_i^{}\right)·{f_{\eta }^{\prime }}^{\left(i\right)}·{\left[\frac{\left(f_r+f_c\right)}{c}\right]}^{\left(1-i\right)}}$$

(11)

Clearly, the range, azimuth, and cross-coupling spatial variations of A are all included in the third term in Equation (11) and they should be eliminated to guarantee the image qualities. Setting the terms involving the spatially variant coefficients ${\xi }_i\left(A\right)\left(i=2,3,4\right)$ to zero is a crucial operation for the 2-D scaling algorithm, i.e., ${\left[\begin{array}{ccc}{\xi }_2\left(A\right)-{\omega }_2& {\xi }_3\left(A\right)-{\omega }_3& {\xi }_4\left(A\right)-{\omega }_4\end{array}\right]}^T={\left[\begin{array}{ccc}0& 0& 0\end{array}\right]}^T$. According to the detailed expressions ${\xi }_i^{}\left(i=2,3,4\right)$ in Appendix B, the following equations are obtained:

$$\left[\begin{array}{l}{\xi }_2\left(A\right)-{\omega }_2\\ {\xi }_3\left(A\right)-{\omega }_3\\ {\xi }_4\left(A\right)-{\omega }_4\end{array}\right]=\left[\begin{array}{ccc}3{\beta }_4\left(A\right)/2& -3{\beta }_3\left(A\right)/4& {\beta }_2\left(A\right)/4\\ 0& {\beta }_4\left(A\right)/2& -{\beta }_3\left(A\right)/8\\ 0& 0& {\beta }_4\left(A\right)/16\end{array}\right]\left[\begin{array}{c}{\left\{{\kappa }_1\left(A\right)\right\}}^2\\ {\kappa }_1\left(A\right)\\ 1\end{array}\right]-\left[\begin{array}{c}{\omega }_2\\ {\omega }_3\\ {\omega }_4\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$$

(12)

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

$$\{\begin{array}{l}{\beta }_2\left(A\right)=G_0+G_1·\left\{{\kappa }_1\left(A\right)-{\kappa }_1\left(C\right)\right\}+G_2·{\left\{{\kappa }_1\left(A\right)-{\kappa }_1\left(C\right)\right\}}^2\\ {\beta }_3^A=Y_0+Y_1·\left\{{\kappa }_1\left(A\right)-{\kappa }_1\left(C\right)\right\}\\ {\beta }_4^A=M_0\end{array}$$

(13)

where $G_0=4{\omega }_2-24{\kappa }_1\left(C\right){\omega }_3+6{\beta }_4\left(A\right){\left\{{\kappa }_1\left(A\right)\right\}}^2$, $G_1=12\left({\beta }_4\left(A\right){\kappa }_1\left(C\right)-2{\omega }_3\right)$, $G_2=6{\beta }_4\left(A\right)$, $Y_0=4{\beta }_4\left(A\right){\kappa }_1\left(C\right)-8{\omega }_3$, $Y_1=4{\beta }_4\left(A\right)$, and $M_0=16{\omega }_4$. Then, ${\beta }_2\left(A\right)$, ${\beta }_3\left(A\right)$, and ${\beta }_4\left(A\right)$ are respectively expanded in the power series with respect to ${\kappa }_1\left(A\right)$, i.e.,

$$\{\begin{array}{l}{\beta }_2\left(A\right)={\beta }_2\left(C\right)+{\dot{\beta }}_2\left(C\right)·\left\{{\kappa }_1\left(A\right)-{\kappa }_1\left(C\right)\right\}+{\ddot{\beta }}_2\left(C\right)·{\left\{{\kappa }_1\left(A\right)-{\kappa }_1\left(C\right)\right\}}^2\\ {\beta }_3\left(A\right)={\beta }_3\left(C\right)+{\dot{\beta }}_3{\left(C\right)}_1·\left\{{\kappa }_1\left(A\right)-{\kappa }_1\left(C\right)\right\}\\ {\beta }_4\left(A\right)={\beta }_4\left(C\right)\end{array}$$

(14)

where ${\beta }_2\left(C\right)$, ${\dot{\beta }}_2\left(C\right)$, and ${\ddot{\beta }}_2\left(C\right)$ are expansion coefficients of ${\beta }_2\left(A\right)$, ${\beta }_3\left(C\right)$ and ${\dot{\beta }}_3\left(C\right)$ are the expansion coefficients of ${\beta }_3\left(A\right)$, and ${\beta }_4\left(C\right)$ is the expansion coefficient of ${\beta }_4\left(A\right)$, i.e., ${\dot{\beta }}_n\left(C\right)={d{\beta }_n\left(A\right)/d{\kappa }_1\left(A\right)|}_{{\kappa }_1\left(A\right)={\kappa }_1\left(C\right)}$ and ${\dot{\beta }}_n\left(C\right)=(1/2)d^2{\beta }_n\left(A\right)/{{d\left\{{\kappa }_1\left(A\right)\right\}}^2|}_{{\kappa }_1\left(A\right)={\kappa }_1\left(C\right)}$. Note that both Equations (13) and (14) are used for the expressions of β₂(A), β₃(A), and β₄(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.,

$$\begin{array}{l}{G}_0={\beta }_2\left(C\right),G_1={\dot{\beta }}_2\left(C\right),G_2={\ddot{\beta }}_2\left(C\right)\\ Y_0={\beta }_3\left(C\right),Y_1={\dot{\beta }}_3\left(C\right),M_0={\beta }_4\left(C\right)\end{array}$$

(15)

Solving this equation set in Equation (15), we can derive its solutions as

$$\{\begin{array}{l}{\varsigma }_2=\left(-P_1^2-P_3\right)/P_2-{\kappa }_2\left(C\right)\\ {\varsigma }_3=\left(-P_1^3-P_3{P}_1\right)/P_2-{\kappa }_3\left(C\right)\\ {\varsigma }_4=\left(P_3^2-P_1^2{P}_3\right)/P_2-{\kappa }_4\left(C\right)\end{array}$$

(16)

$$\{\begin{array}{l}{\omega }_2=24{\varsigma }_4{\left\{{\kappa }_1\left(C\right)\right\}}^2+6{\varsigma }_3{\kappa }_1\left(C\right)-P_2/\left(8\left(P_1^2+P_3\right)\right)\\ {\omega }_3=8{\varsigma }_4{\kappa }_1\left(C\right)-P_1{P}_2^2/\left(48{\left(P_1^2+P_3\right)}^2\right)\\ {\omega }_4=-\left(3P_1^2+P_3\right)P_2^3/\left(384{\left(P_1^2+P_3\right)}^3\right)\end{array}$$

(17)

where $P_1={d{\kappa }_2\left(A\right)/d{\kappa }_1\left(A\right)|}_{{\kappa }_1\left(A\right)={\kappa }_1\left(C\right)}$, $P_2={\left(1/2\right)d^2{\kappa }_2\left(A\right)/d{\left\{{\kappa }_1\left(A\right)\right\}}^2|}_{{\kappa }_1\left(A\right)={\kappa }_1\left(C\right)}$, and $P_3={d{\kappa }_3\left(A\right)/d{\kappa }_1\left(A\right)|}_{{\kappa }_1\left(A\right)={\kappa }_1\left(C\right)}$. After the 2-D scaling processing in the hybrid domain by using Equations (16) and (17), the echo signal can be expressed as

$$S_4\left(f_r,f_{\eta }^{\prime }\right)=\exp \left[j2\pi \left({\xi }_0\left(A\right)·\frac{\left(f_r+f_c\right)}{c}+{\xi }_1\left(A\right)·f_{\eta }^{\prime }\right)\right]$$

(18)

Clearly, the phase in Equation 18 is a 2-D linear one with respect to the range frequency $f_r$ and azimuth frequency $f_{\eta }^{\prime }$. 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.

4. Simulation Results

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

4.1. 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.

Figure 7a–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.

Figure 8 shows the comparative results of the reference targets which are extracted from the simulated results in Figure 7. Clearly, the 2-D impulse response by the proposed model is well-focused with relative clear separations of the main lobe and first and subsequential sidelobes, whereas the ones from the traditional non-‘Stop-Go’ and ‘Stop-Go’ models are defocused completely, which indicates that the impacts of the accelerations on the ‘Stop-Go’ model must be taken into consideration when design the focusing algorithm.

4.2. 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. Parameter settings for Case II.

	Platform	Transmitter	Receiver
Parameter
Radar Position at ACM		(0, 0, 510) km	(112, −78, 25) km
Velocity Vector		(0, 7600, 0) m/s	(−170, 800, −640) m/s
Acceleration Vector		/	(13, −34, −68) m/s2
Elevation Angle		81.3°	27.7°
Carrier Frequency		9.65 GHz
Pulse Width		20 μs
Pulse Bandwidth		240 MHz
Sampling Frequency		360 MHz
Azimuth Resolution		0.716

. 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.

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. Image quality parameters.

		Range			Azimuth
Method	Target	IRW (m)	PSLR (dB)	ISLR (dB)	IRW (m)	PSLR (dB)	ISLR (dB)
Proposed	PT1	0.626	−13.21	−9.97	0.723	−13.09	−9.91
	PT3	0.625	−13.19	−10.01	0.712	−13.04	−9.84
NCSA	PT1	0.625	−13.20	−9.99	0.983	−7.04	−7.36
	PT3	0.627	−13.18	−9.98	0.965	−6.83	−7.28
GPFA	PT1	0.626	−13.20	−10.02	1.124	−10.12	−8.07
	PT3	0.628	−13.19	−9.99	1.059	−12.25	−6.63

. 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.

5. Discussion

5.1. 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$\pi$ and the extra phase error brought about by the acceleration term is about 0.007$\pi$, 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

$$H_C\left(f_r,f_{\eta }\right)=\exp \left(j2\pi \frac{f_r{f}_{\eta }}{\gamma }+j\pi \frac{f_{\eta }^2}{\gamma }\right)$$

(19)

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

$$H_{2D-rc}\left(f_r,f_{\eta }\right)=\exp \left\{j\pi \frac{{\left(f_r+f_{\eta }\right)}^2}{\gamma }\right\}$$

(20)

5.2. 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+3NM{\log }_{2}M/2+4NM$, 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+NM{\log }_{2}M+3NM$.

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. Computational load analysis.

Data Size in Samples (×103)	1	2	4	8	16
CSA/Proposed	0.793	0.794	0.794	0.795	0.795
BPA/Proposed	35.59	91.25	194.81	387.78	750.61

).

6. 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.