An integrative synchronization and imaging approach for bistatic spaceborne/stratospheric SAR with a fixed receiver

Bistatic spaceborne/stratospheric synthetic aperture radar (SAR) with a fixed receiver is a novel hybrid bistatic SAR system, in which a spaceborne SAR serves as the transmitter of opportunity, while a fixed receiver is mounted on a stratospheric platform. This paper presents an integrative synchronization and imaging approach for this particular system. Firstly, a novel synchronization method using the direct-path signal, which can be collected by a dedicated antenna, is proposed and applied. The synchronization error can be completely removed using the proposed method. However, as the cost of synchronization, the characteristic of synchronized echo’s range history becomes quite different from that of general bistatic SAR data. To focus this particular synchronized data, its 2-D spectrum is derived under linear approximations and then a frequency-domain imaging algorithm using two-dimensional inverse scaled Fourier transform (2-DISFT) is proposed. At last, the proposed integrative synchronization and imaging algorithm is verified by simulations.

Bistatic synthetic aperture radar (SAR) has been an active research area in the last decade, where in particular, the method based on spaceborne illuminator appears to be more attractive. Several experiments of bistatic SAR with spaceborne illuminator have been conducted by numerous organizations. In these experiments, the receiver is either spaceborne, aircraft, or stationary on the ground. The promising imagery results show the great potentials and capabilities of bistatic SAR as an innovative imaging system [1–6].

This paper discusses a particular sub-class of spaceborne hybrid configuration, where a spaceborne SAR, e.g. TerraSAR-X, serves as the transmitter of opportunity, while a fixed receiver is mounted on a stratospheric platform, e.g. a stratospheric aerostat. Its advantages include low-cost, reduced vulnerability to counter measurement, high operational flexibility and wide observation scene [7, 8]. Such a system is a good tool for high-resolution imaging and high-precision height information extraction, which has the potential for future mission both in civil and military applications.

Time and phase synchronization is essential and foremost for such a kind of bistatic system, since independent local oscillators (LO) are used, and the common time reference is missed for the transmitter and the receiver. The synchronization errors result not only in range cell migration (RCM) error but also in the distortion of the azimuth dependent phase history [9–12]. This implies that synchronization errors will degrade the quality of the bistatic SAR images, and therefore, corresponding compensation must be implemented.

The synchronization scheme has been well studied for cooperative bistatic SAR systems, in which the dedicated synchronization link can be constructed [13–17]. An echo-domain phase synchronization approach by using correlation of bistatic echo is proposed for bistatic SAR in alternating bistatic/ping-pong mode in [18]. For un-cooperative bistatic SAR systems, it is a common method to achieve synchronization by using direct-path signal [19–23].

Assuming that the direct-path channel and the reflected-path channel are previous balanced and the common LO is used for both channels, the reflected signal and the direct-path signal will be contaminated by the same synchronization errors. However, without any other auxiliary data, it is quite difficult to estimate the phase synchronization error with high precision [22]. On the one hand, the accuracy of the satellite's trajectory is not sufficient for separating the phase synchronization error from the nominal phase caused by direct-path range history. On the other hand, the phase synchronization error caused by the LO phase noise is a random component, which is difficult to estimate [24].

In this paper, we propose a new synchronization method using direct-path signal's time delays and peak phases, without the isolation of synchronization errors, to compensate corresponding components in the reflected signal directly. In this way, the time and synchronization errors can be completely removed, and this method is quite simple and fast.

However, due to the fact that direct-path signal’s time delays and peak phases constitute both the synchronization error component and the nominal range component, the range history of direct-path signal will be removed in the synchronized echo data as well. This means that the system impulse response of the synchronized data will be quite different from that of the general bistatic SAR data. It is most straightforward to use time-domain algorithms such as back-projection algorithm (BPA) to focus the synchronized data [25–27]. However, it requires high-precision position measurement and suffers from severe computational load. Therefore, a frequency-domain algorithm is the preferred choice, but has to be re-engineered accordingly.

The calculation of an analytical bistatic point-target reference spectrum (BPTRS) is the key to develop frequency-domain imaging algorithm, since the bistatic range history loses its hyperbolic form [28]. Based on the approximate BPTRS, including Loffeld bistatic formula (LBF) [29, 30] and series reversion [31], several image formations have been proposed. The 2-DISFT [32] algorithm and chirp scaling (CS) algorithm [33] were developed from LBF for corresponding bistatic configurations. Based on the series reversion method, a non-linear CS (NLCS) process was applied to equalize the azimuth chirp rate in bistatic SAR focusing for general configuration [34], and a range Doppler (RD) algorithm was proposed to handle the azimuth-invariant bistatic case [35]. In bistatic SAR with a fixed receiver configuration, a two-dimensional NLCS processor was applied in [36] to deal with the large bistatic angle case. In [37], a highly accurate bistatic range migration algorithm (RMA) was proposed for asymmetric bistatic SAR system with a fixed receiver. A modified range Doppler algorithm was presented for space-surface bistatic SAR (SS-BSAR) [38]. In the spaceborne/airborne configuration, a frequency-domain processing method based on 2-DISFT was proposed in [39, 40].

According to the property of the synchronized data, a frequency-domain imaging algorithm based on 2-DISFT is proposed in this paper. Firstly, for the sake of triple square-root terms in the range history, the Taylor series is applied to obtain the stationary phase point when deriving the BPTRS. Further, under the properly designed linear approximation of BPTRS, the 2-D spectrum of the synchronized scene data is derived. It is found in the 2-D spectrum that the dominant component is the 2-D scaled Fourier transform of target's bistatic backscattering coefficient. At last, a frequency-domain imaging algorithm based on 2-DISFT is proposed.

This paper is organized as follows. In Section 2, the geometry of the bistatic spaceborne/stratospheric SAR with a fixed receiver is given, and the signal model of the echo data with synchronization errors is derived. In Section 3, the synchronization implementation using direct-path signal is presented. Section 4 derives the BPTRS and 2-D spectrum for synchronized data. The analysis of approximation errors for the derivation of 2-D spectrum is presented in Section 5. Then, a frequency-domain imaging algorithm using 2-DISFT is proposed in Section 6. To validate the proposed algorithm, simulation experiments are carried out in Section 7. Finally, in Section 8, some conclusions are presented.

Figure 1 shows the geometry of the bistatic spaceborne/stratospheric SAR with a fixed receiver considered, where the spaceborne transmitter T moves along straight line with velocity v, while the stratospheric receiver R keeps stationary. It is worth pointing out that the receiver has been designed with two antennas. One is the radar channel antenna which is used to collect reflected signal, while the other is the direct-path channel antenna which is used to collect direct-path signal. As can be seen in the Figure 1, for the point target P, the reflected signal travels along r _T + r _R to be received by the radar channel antenna, while the direct-path signal travels along r _D to be received by the dedicated direct-path channel antenna.

Geometry of the bistatic spaceborne/stratospheric SAR with a fixed receiver.

Full size image

As shown in Figure 1, the coordinate of the receiver is (0,0,z _R), and the coordinate of the transmitter is (x _T,y _T,z _T). According to the definition of geometry, z _R, x _T, and z _T are constant, and y _T(t) satisfies

where t is the azimuth (slow) time. Assuming that the coordinate of the point target P is (x, y, 0), then the transmitter-to-target range and the target-to-receiver range for P can be expressed as

To make the following derivation easier, we use the transmitter-referenced coordinate to define the coordinate of target space. Assuming that

where t _0T is the zero Doppler time for P, r _0T is the slant range from transmitter to target at t _0T. Therefore, (2) and (3) can be rewritten as

where $r_{0\mathrm{R}}=\sqrt{x^2+z_{\mathrm{R}}^2}$. From (5), we can see that x is the function of the variable r _0T; hence, r _0R is expressed as r _0T(r _0T).

Meanwhile, the slant range history of direct-path signal can be written as

where $r_{0\mathrm{d}}=\sqrt{x_{\mathrm{T}}^2+{\left(z_{\mathrm{T}}-z_{\mathrm{R}}\right)}^2}$ is the closest distance from the transmitter to the receiver.

Supposing that the transmitted signal is

where rect[·] is the window function with rectangular shape, τ is the range (fast) time, T _P is the duration time of the transmitted pulse, f ₀ is the radar carrier frequency, and k is the chirp rate. Therefore, the ideal reflected raw signal from the point target P at the receiver can be given as

where σ(·) is the bistatic backscattering coefficient of the point target P, t _delay = (r _T(t,t _0T,r _0T) + r _R(t _0T,r _0T))/c is the time delay corresponding to the time it takes the signal to travel the transmitter-target-receiver distance, c is the speed of light, and T _S is the synthetic aperture time for the point target P.

In bistatic satellite/stratosphere SAR system, there is no common time reference between the transmitter and the receiver. Moreover, independent oscillators are used in the transmitter and the receiver, thus the phase noise of the oscillator cannot be cancelled out as in monostatic SAR. Therefore, time and phase synchronization errors must be considered in this system.

Firstly, only considering that the time synchronization error is e(t), the practical reflected signal can be rewritten as

Next, assuming that the oscillator in the spaceborne transmitter has an ideal performance, the combined phase synchronization error ∅ _e(t) is assigned to the oscillator in the receiver [14–22]. Thus, omitting the amplitude and the initial phase which are inessential, the output signal of the receiver's oscillator can be modelled as

After the quadrature demodulation process, the practical signal reflected from the point target P, considering time and phase synchronization errors, can be given as

As can be seen from (13), the existing time and phase synchronization errors result not only in a drift of the echo sampling windows which will cause RCM errors but also in the distortion of the azimuth-dependent phase history. From the knowledge of SAR imaging, we can conclude that the time and phase synchronization errors will degrade the quality of bistatic SAR image [9].Therefore, the time and phase synchronization compensation must be implemented.

Assuming that the direct-path channel and the radar channel are previously balanced and the common LOS is used for both channels, the radar channel signal and the direct-path channel signal will be contaminated by the same synchronization errors. Thus, the direct-path signal after demodulation can be written as

where t _{D_delay} = r _D(t)/C.

As stated before, extracting the synchronization errors from direct-path signal is a straightforward idea for the bistatic SAR synchronization. It is possible to extract time synchronization error from direct-path signal [22]. However, it seems to be difficult to estimate and extract the phase synchronization error with high precision if no other auxiliary data are available. The reason for this is double sided. On one hand, the accuracy of the satellite's trajectory is not sufficient to separate the phase synchronization error from the nominal phase caused by direct-path range history. On the other hand, the phase synchronization error caused by the LO phase noise is a random component, which is difficult to be estimated [24]. Therefore, a different approach has to be applied here.

Taking into account the advantages of high SNR and clean phases, after range compression, time delays and peak phases of the direct-path signal can be precisely extracted. Ignoring trivial constant offsets, the extracted time delays and peak phases can be formulated as

As can be seen from (15) and (16), the extracted time delays and peak phases constitute both the nominal range component and the synchronization error component. Instead of separating the synchronization error component from the nominal range component, we use the whole extracted time delay $\stackrel{\wedge }{t_{\mathrm{D}\_\mathrm{delay}}\left(t\right)}$ and peak phase $\stackrel{\wedge }{{\varphi }_{\mathrm{D}}\left(t\right)}$ to compensate the corresponding terms in the reflected signal. This can be done by applying an opposite time-delay shift and by multiplying an opposite phase term on each pulse.

At first, performing the range FT with respect to the variable τ, we can transform (13) into the range-frequency/azimuth-time domain:

where f is the range frequency. Then, according to (15), establishing the reference function of time synchronization yields

Thus, the reflected signal after time synchronization, in the range-frequency/azimuth-time domain, can be given as

and then transforming (19) into 2-D time domain gives

Similarly, the reference function of phase synchronization according to (16) can be given as

Finally, multiplying (20) by the reference function (21) yields

where

The combination of extraction and compensation can be viewed as the whole synchronization process. Its block diagram is shown in Figure 2.

The block diagram of synchronization process.

Full size image

From (22), we can see that after the implementation of synchronization, synchronization errors are completely removed. However, as the cost of synchronization, the range history of the point target P r _T + r _R is replaced by r _T + r _R - r _D. We note here that the range history is quite different from that of the general bistatic SAR. Thus, some special processing approaches have to be implemented, which will be discussed in the following sections.

To understand the features of the synchronized data, this section will describe the processing steps performed to obtain its 2-D spectrum, which is based on the signal model (22) and the principle of stationary phase (POSP).

4.1 Derivation of BPTRS

The BPTRS can be obtained by performing the Fourier transform (FT) with respect to the variable τ and the variable t, respectively. Firstly, performing the range FT with respect to the variable τ, we can transform (22) into the range-frequency/azimuth-time domain. Using POSP, the result can be given as

Then we can transform (24) into the range-frequency/azimuth-frequency domain by performing the azimuth FT with respect to the variable t.

where f _a is the azimuth frequency, and the phase term considered in the integral is

It can be seen from (6) to (8) and (23) that there are three square root terms in φ _a(f _a,f,t,t _0T,r _0T). This makes it difficult to obtain BPTRS by applying POSP, since it seems to be impossible to obtain the analytical expression of the stationary point for (26). To overcome this limitation, the Taylor series expansion is used to simplify the phase terms in (26), with a sufficient precision of approximation which has been verified in [39]. Expanding the three square root terms in Taylor series and keeping terms up to the second order gives

Based on (27) to (29), we can calculate some key parameters of the Doppler signal. Assuming that the transmitted signal is the narrowband signal, i.e. $\frac{f+f_0}{c}\approx \frac{1}{\lambda }$, the Doppler centroid and the Doppler bandwidth can be derived by

Therefore, the Doppler time-bandwidth production (TBP) is ${\mathrm{TBP}}_{\mathrm{D}}=\frac{v^2{T}_{\mathrm{s}}^2}{\lambda r_{0\mathrm{T}}{r}_{0\mathrm{d}}}·\left|r_{0\mathrm{T}}-r_{0\mathrm{d}}\right|$. For the point target located in the scene centre, we get v = 7,600 m/s, λ = 0.031 m, r _0T = 726.9 km, r _0d = 645.8 km, and T _s = 0.484susing the parameters listed in the Table 1 and then TBP_D ≈ 78. This means that the accuracy of POSP is sufficient to obtain the analytical expression of the spectrum. Thus, the phase term shown as (26) can be rewritten as

Setting $\frac{\partial {\phi }_{\mathrm{a}}\left(f_{\mathrm{a}},f,t,t_{0\mathrm{T}},r_{0\mathrm{T}}\right)}{\partial t}=0$, we can obtain the stationary point

Substituting t _P for t in the integral of (25) yields

where the phase term ψ (f _a ,f,t _0T,r _0T), which represents the significant features of BPTRS, can be given as

The first and the second terms in (35) represent the RCM which is not only azimuth variant but also range-variant. The third term is a linear function of azimuth frequency f _a, expressing the azimuth position of the point target in the focused image. The forth term stands for the cross coupling between azimuth and range, and the last one is responsible for the range modulation. We can see from (35) that BPTRS is dependent on both the azimuth coordinate t _0T and the range coordinate r _0T of the point target P.

4.2 Derivation of 2-D spectrum of the scene data

The 2-D spectrum of a complete scene data can be obtained by integrating the BPTRS over all the point target spectra [41, 42]:

As shown in (35), the phase term ψ (f _a,f,t _0T,r _0T) is a complicated function with respect to the position of target, i.e., variables t _0T and r _0T.To obtain the analytical solution of the integral, the linearity approximation of ψ (f _a,f,t _0T,r _0T) with regard to variables t _0T and r _0T has to be carried out.

First of all, even though both the first term and the second term of ψ (f _a,f,t _0T,r _0T) represent the RCM shift, the second term is a quadratic function of t _0T. Therefore, the second term will be automatically omitted in the linearity approximation expression. The impact of this approximation will be given in the following section.

Then, by expanding (35) in Taylor series at t _0T = 0 and r _0T = r ₀, and only keeping terms up to the first order gives

where $r_0=\sqrt{{\left(x_0-x_{\mathrm{T}}\right)}^2+z_{\mathrm{T}}^2}$, (x ₀,0) is the coordinate of the target located at the scene centre, and

where (t _0T = 0, r _0T = r ₀) expresses the scene centre point. In (37), the first term is the space-invariant phase component, representing the space-invariant range modulation, RCM, and azimuth modulation. The last two terms in (37) represent the space-variant component of ψ _L (f _a,f,t _0T,r _0T). After some tedious algebra manipulations, the corresponding result can be given as

where

We can see that there is across coupling between azimuth and range in (43). To remove the coupling, another linear approximation can be made to (43). By expanding ψ _rL (f _a,f,0,r ₀) in Taylor series of f and keeping terms up to the first order, ψ _rL (f _a,f,0,r ₀) becomes

where

where λ = c/f ₀ is the wavelength. In (45), the first term is a pure Doppler phase term which represents the range-variant azimuth modulation, and the second term is the scaled range-frequency phase term which expresses the residual space-variant RCM component.

At last, letting r = r _0T - r ₀ and substituting (37) and (45) into (36) yields the2-D spectrum of the synchronized scene data

where

The exponential term in the integral of (48), as it is linear with respect to t _0T and r _0T, can be interpreted as a Fourier kernel and thus we can obtain

However, due to the expressions of ψ _tL(f _a, f, 0, r ₀) and ψ _rL(f _a, f, 0, r ₀), σ(ψ _tL, ψ _rL) should be viewed as the 2-D scaled Fourier transform of the bistatic backscattering coefficient σ(t _0T, r _0T).

Examining the derivation of (50), it can be seen that two linear approximations are made to the phase term ψ (f _a,f,t _0T,r _0T). To validate these operations, these approximation errors should be analysed in detail.

Firstly, as shown in (45), ψ _rL (f _a,f,0,r ₀) is approximated by a linear function with respect to the variable f. This implies that the second-order term and higher-order terms are neglected with respect to the variable f. According to (43) and (45), the approximation error can be given as

For the case of narrowband signal (f ≪ f ₀), it is clear that the quadratic term is the most dominant component. Therefore, only the quadratic term needs to be considered in the following analysis. Considering the worst case when the range frequency f reaches its maximum value, the approximation phase error is calculated and shown in Figure 3.

Approximation phase error, denoted by ( 51 ).

Full size image

Using the parameters listed in Table 1, we get r ₀ = 726.9 km, r _0d = 645.8 km, f = 25 MHz. The azimuth frequency f _a and the slant range deviation r = r _0T - r ₀ are used as independent variables in the simulation. From Figure 3, we can see that |Δψ _E1| ≪ π/8 is satisfied. This implies that the linear approximation operation of (45) is accurate enough.

Secondly, as shown in (37), ψ (f _a,f,t _0T,r _0T) is approximated by ψ _L(f _a, f, t _0T, r _0T) which is a linear function with respect to variables t _0T and r _0T. Therefore, the approximation error can be expressed as

Substituting (35) and (37) into (52) gives

In (53), the first term is due to the azimuth displacement error. The second term and the third term represent the residual RCM error and may result in a range displacement error. The last term is the residual coupling phase which may result in a defocus.

From (53), the azimuth displacement error can be given as

From (35), the correct range position of the point target in the focused image can be given as

where vt _0T represents the azimuth coordinate of the point target, indicating the azimuth-variant RCM shift. However, after this approximation, the range position of the point target in the focused image will be

Thus, the range displacement error can be expressed as

We can see that the geometric distortion will be introduced by this linear approximation. However, this geometric distortion can be corrected by interpolating after image focusing [43]. Moreover, the geometric correction is a necessary step for SAR imaging in frequency domain, and thus this process will not increase the computational load.

Therefore, only the last phase error term in (53) is considered here. Under the assumption of narrowband signal, e.g. c/f + f ₀ ≈ λ, expanding it in Taylor series at r _0T = r ₀ gives

In (58), the ratio of the second term against to the first term is given as $\frac{1}{2}·\left|\frac{r_{0\mathrm{T}}-r_0}{r_{0\mathrm{d}}-r_0}\right|$, which is smaller than 0.04 for the same parameters as before. Furthermore, according to the principle of Taylor expansion, the higher the order of the term, the smaller the error value. Now, we can see that the dominant phase error is the quadratic function of the f _a and the r = r _0T - r ₀. Thus, the phase error increases quickly when both f _a and r increase. If we use π/8 as the tolerable threshold, the constraint for this approximation can be given as

Using the same parameters as before, we get |f _a r| ≤ 4.87 × 10⁵. For example, if the scene width is assumed to be 10 km in the slant-range direction, i.e. max(r) = 5,000 m, the azimuth frequency needs to satisfy |f _a| ≤ 97.4 Hz.

However, from (34), we can see that the f _a is restricted by $\mathrm{rect}\left[\frac{f_{\mathrm{a}}-f_{\mathrm{DC}}}{B_{\mathrm{a}}}\right]$. Moreover, from (30), we note that the Doppler centroid f _DC is dependent on the azimuth coordinate t _0T of the point target. Assuming that the scene width is W _a in the azimuth direction, and then the azimuth frequency of the synchronized scene data meets

Using the parameters listed in Table 1, we get B _a = 158.99 Hz. If we assume that W _a = 2,000 m; thus, max(|f _a|) = 458.02 Hz. However, it is beyond the restriction given above. This means that a compromise between f _a and r is needed to meet the requirement.

According to the discussion above, Figure 4 gives an example of compromise, in which max(|f _a r|) = 4.5 × 10⁵. The size of the corresponding scene block is about 1 km × 3 km (azimuth direction × slant-range direction).

Approximation phase error, denoted by ( 58 ).

Full size image

If we intend to focus a larger scene, the synchronized data can be processed in small blocks both on range direction and on azimuth direction. For different range blocks, the corresponding reference slant range r ₀ is used. For different azimuth blocks, Doppler centroid correction (DCC) is applied by multiplying the synchronized echo with the DCC function in the 2-D time domain. The DCC function can be given as

where

In (62), t _0TN is the azimuth coordinate of the centre point in the block N. According to the shifting theorem of FT, the corresponding 2-D spectrum can be given as

In this way, the 2-D spectrum of the block N will satisfy the constriction denoted by (59). However, blocks have to overlap, and then the efficiency of the algorithm decreases with increasing the number of blocks.

6.1 Analysis of 2-D spectrum

As is pointed out, σ(ψ _tL,ψ _rL) can be interpreted as a 2-D scaled Fourier transform of the bistatic backscattering coefficient σ(t _0T,r _0T); thus, we rewrite σ(ψ _tL,ψ _rL) as

It can be seen that both the azimuth frequency and range frequency are scaled by coefficients. In particular, the azimuth frequency is scaled by |r _0d/(r _0d - r ₀)|. Meanwhile, from (31), we can see that the Doppler bandwidth can be rewritten as

where $B_{\mathrm{a}}^{\text{'}}=\frac{v^2}{\lambda r_{0\mathrm{T}}}·T_{\mathrm{s}}$ is the Doppler bandwidth for the traditional bistatic SAR with a fixed receiver configuration. In this case, if we apply inverse Fourier transform along azimuth direction directly, the azimuth resolution of inferred image will be ${\rho }_{\mathrm{a}}=\left|\frac{r_{0\mathrm{d}}}{r_{0\mathrm{d}}-r_0}\right|·\frac{v}{B_{\mathrm{a}}^{\text{'}}}$. Compared with the azimuth resolution of the traditional bistatic SAR with a fixed receiver configuration, which is ${\rho }_{\mathrm{a}0}=v/B_{\mathrm{a}}^{\text{'}}$, the inferred azimuth resolution is scaled by |r _0d/(r _0d - r ₀)|. When the value of |r _0d/(r _0d - r ₀)| is many times larger than 1, the inferred image would suffer a very poor azimuth resolution. Using the same parameters as before, we get ρ _a0 = 5.31 m, |r _0d/(r _0d - r ₀)| = 7.96, and then ρ _a0 = 4.82 m. Similarly, the range frequency suffers the same problem, but not as serious as the azimuth one. To circumvent this limitation, we propose to apply ISFT both on range and on azimuth direction, which will be presented in the next section.

6.2 Principle of ISFT

The inverse scaled Fourier transform, abbreviated as ISFT, was introduced to focus monostatic SAR data in [41] and was compared with chirp scaling method in [42]. The ISFT was also applied for bistatic SAR data imaging [32, 40].

The principle of ISFT is shown below [41]. As shown in Figure 5, assuming s(t)↔S(f) is a Fourier transform pair, if the input frequency-domain signal is S(a · f), after ISFT implementation, the output time-domain signal will be $\frac{1}{\left|a\right|}·s\left(t\right)$. The digital implementation of ISFT can be found in [42].

The principle of ISFT.

Full size image

6.3 Imaging process

From (50), we can see that the 2-D spectrum of the synchronized scene data can be viewed as the consequence of a multiplication of the space-invariant phase term H ₀(f _a,f) and the 2-D scaled Fourier transform of the 'brightness’ of the point target σ(ψ _tL(f _a, f, 0, r ₀), ψ _rL(f _a, f, 0, r ₀)). Meanwhile, the imaging can be interpreted as the process of obtaining the backscattering coefficient σ(t _0T,r _0T) from the 2-D spectrum H(f _a,f). Therefore, in this paper, a 2-DISFT imaging algorithm is proposed to obtain the focused image. The implementation steps are listed below, and the block diagram of proposed algorithm is shown in Figure 6.

The block diagram of 2-DISFT imaging process.

Full size image

1. 1.

   According to the constriction (59), the whole synchronized scene data is divided into smaller blocks.

2. 2.

   DCC was performed by multiplying the DCC function in 2-D time-domain for different azimuth blocks.

   $$s_{\mathrm{DCC}}\left(t,\tau ,t_{0\mathrm{T}},r_{0\mathrm{T}}\right)=s\left(t,\tau ,t_{0\mathrm{T}},r_{0\mathrm{T}}\right)·H_{\mathrm{DCC}}\left(t\right).$$

   (66)

1. 3.

   The 2-D time-domain data was transformed into 2-D frequency domain using Fourier transform.

2. 4.

   Bulk RCM correction (RCMC) and 2-D compression by multiplying the reference function (RF). As is pointed out, H ₀(f _a, f) represents the space-invariant range modulation, RCM and azimuth modulation. Therefore, RF should be given as the conjugate of H ₀(f _a,f). However, to keep the phase which can be used in the interferometric application, RF can be expressed as

   $$\begin{array}{ll}{H}_{\mathrm{RF}}\left(f_{\mathrm{a}},f\right)& =H_0^{\ast }(f_{\mathrm{a}},f)·exp\left[-j2\pi f_0·\frac{r_0+r_{0\mathrm{R}}\left(r_0\right)-r_{0\mathrm{d}}}{c}\right]\\ =\mathrm{rect}\left[\frac{f}{k·T_{\mathrm{p}}}\right]·exp\left[\mathit{j\pi }\frac{f^2}{k}\right]\\ ·exp\left[j2\mathit{\pi f}·\frac{r_0+r_{0\mathrm{R}}\left(r_0\right)-r_{0\mathrm{d}}}{c}\right]\\ ·exp\left[-\mathit{j\pi }·\frac{f_{\mathrm{a}}^{2}c{r}_0{r}_{0\mathrm{d}}}{v^2\left(r_{0\mathrm{d}}-r_0\right)\left(f+f_0\right)}\right].\end{array}$$

   (67)

   After this implementation, the remaining signal can be given by

   $$\begin{array}{ll}{H}_1\left(f_{\mathrm{a}},f\right)& =H(f_{\mathrm{a}},f)·H_{\mathrm{RF}}(f_{\mathrm{a}},f)\\ =exp\left[-j2\pi f_0·\frac{r_0+r_{0\mathrm{R}}\left(r_0\right)-r_{0\mathrm{d}}}{c}\right]\\ ·\sigma \left[\frac{r_{0\mathrm{d}}}{r_{0\mathrm{d}}-r_0}·f_{\mathrm{a}},{\psi }_{\mathrm{rL}1}\left(f_{\mathrm{a}},0,r_0\right)+{\psi }_{\mathrm{rL}2}\left(f_{\mathrm{a}},0,r_0\right)·f\right].\end{array}$$

   (68)
3. 5.

   ISFT with regard to the range frequency, which can be expressed as

   $$\begin{array}{ll}{H}_2\left(f_{\mathrm{a}},r\right)=& {\displaystyle \int H_1\left(f_{\mathrm{a}},f\right)}exp\left[j2\pi E_{\mathrm{RFS}}\left(f_{\mathrm{a}},0,r_0\right)·f_r·r\right]\\ \times d\left(E_{\mathrm{RFS}}\right(f_{\mathrm{a}},0,r_0)·f_{\mathrm{r}})=\sigma \left[\frac{r_{0\mathrm{d}}}{r_{0\mathrm{d}}-r_0}·f_{\mathrm{a}},r\right]\\ ·exp\left[-j2\pi {\psi }_{\mathrm{rL}1}\left(f_{\mathrm{a}},0,r_0\right)r\right]\\ ·exp\left[-j2\pi f_0·\frac{r_0+r_{0\mathrm{R}}\left(r_0\right)-r_{0\mathrm{d}}}{c}\right],\end{array}$$

   (69)

   where

   $$\begin{array}{ll}{E}_{\mathrm{RFS}}\left(f_{\mathrm{a}},r_0\right)& =\frac{c}{1+M}·{\psi }_{\mathrm{rL}2}\left(f_{\mathrm{a}},0,r_0\right)\\ =1+\frac{c}{1+M}·\frac{f_{\mathrm{a}}^2\lambda r_{0\mathrm{d}}^2}{2v^2{\left(r_{0\mathrm{d}}-r_0\right)}^2{f}_0}.\end{array}$$

   (70)

   In the implementation of (69), it is worth noting that f _r is the frequency variable with respect to the variable r, while f is the frequency variable with respect to the fast-time variable τ. The relation between f _r and f is $f_{\mathrm{r}}=f·\frac{1+M}{c}$.

4. 6.

   Residual azimuth compression by multiplying the range-variant phase functions in range-time/azimuth-frequency domain.

   $$\begin{array}{ll}{H}_3\left(f_{\mathrm{a}},r\right)& =H_2(f_{\mathrm{a}},r)·exp\left[j2\pi {\psi }_{\mathrm{rL}1}\left(f_{\mathrm{a}},0,r_0\right)r\right]\\ =\sigma \left[\frac{r_{0\mathrm{d}}}{r_{0\mathrm{d}}-r_0}·f_{\mathrm{a}},r\right]\\ ·exp\left[-j2\pi f_0·\frac{r_0+r_{0\mathrm{R}}\left(r_0\right)-r_{0\mathrm{d}}}{c}\right].\end{array}$$

   (71)
5. 7.

   ISFT with regard to the azimuth frequency f _a to remove azimuth scaling.

   $$\begin{array}{ll}{s}_{\mathrm{out}}\left(t_{0\mathrm{T}},r\right)& ={\displaystyle \int H_3\left(f_{\mathrm{a}},r\right)}exp\left[j2\pi \frac{r_{0\mathrm{d}}}{r_{0\mathrm{d}}-r_0}·f_{\mathrm{a}}·t_{0\mathrm{T}}\right]d\left(\frac{r_{0\mathrm{d}}}{r_{0\mathrm{d}}-r_0}·f_{\mathrm{a}}\right)\\ =\sigma \left(t_{0\mathrm{T}},r\right)·exp\left[-j2\pi f_0·\frac{r_0+r_{0\mathrm{R}}\left(r_0\right)-r_{0\mathrm{d}}}{c}\right]\end{array}$$

   (72)
6. 8.

   According to (54) and (57), applying geometric correction by interpolation

7. 9.

   Stitching the divided blocks into a whole SAR image.

To validate the proposed algorithm, the simulations are carried out. The simulation parameters are listed in Table 1. The transmitter's parameters are referred to TerraSAR-X, and the receiver's parameters are referred to a stratospheric aerostat.

From the knowledge of synchronization, the linear component is the dominant component in the time synchronization error [9]. For the phase synchronization error, both the fixed carrier frequency offset and the phase noise are considered [17, 22]. Therefore, the slope of the time synchronization error is given, and the corresponding parameters of phase synchronization error are also listed in Table 1, where 1 ppm means that the fixed carrier frequency offset is ∆f = f _c  · 10^-6, and the Alan variance σ(τ = 1 s) = 1 × 10^-11 can be regarded as a representative example for the ultra stable oscillators (USO) of current spaceborne SAR system [10].

7.1 Echo characteristics before and after synchronization

At first, the point target located at the scene centre is used to show the characteristic of echo data before and after performing the proposed synchronization implementation algorithm. As shown in Figure 7, before synchronization, the echo support domain is skew both in the time domain and in the frequency domain, and 2-D spectrum appears an azimuth displacement in frequency domain. After synchronization, the echo support domain becomes straight. This will benefit the following focus process. However, it can be seen from Figure 7d that the bandwidth of azimuth frequency is |r _0d/(r _0d - r ₀)| times narrower than the one before synchronization, which verifies the analysis in the former section.

Echocharacteristics of the point target located at scene centre before and after synchronization implementation. (a) Time domain echo before synchronization. (b) Frequency domain echo before synchronization. (c) Time domain echo after synchronization. (d) Frequency domain echo after synchronization.

Full size image

7.2 Imaging simulation for small scene

When the scene size satisfies the constraint condition (59), the proposed algorithm can be applied without block dividing. As is shown earlier, a reasonable scene size is 1 km × 3 km (azimuth direction × slant-range direction). Because the incidence angle of the transmitter is 45°, the scene width along the X-axis (the ground-range direction) is 3 km/cos(45°) ≃ 4.24 km. Here, an imaged scene with size 4 km × 1 km (X-axis × Y-axis) is assumed. According to the geometry depicted in Figure 1 and the parameters listed in Table 1, the coordinate of the point target (Target 5) located in the scene centre is (97.9796 km,0,0). To make the following discussion clearer, as shown in Figure 8, a new X-Y coordinate frame is created with the origin located at (97.9796 km,0,0) in the old frame. Nine point targets are located in the scene, with a 3 × 3 matrix. Then, the proposed algorithm is used to process the simulated echo. The focused image after geometric correction is shown in Figure 9.

The scene with nine point targets for imaging simulation.

Full size image

Imaging simulation result using the proposed algorithm.

Full size image

It can be seen from Figure 9 that all targets are focused precisely and located in the right positions. To further demonstrate the performance of the proposed algorithm, the profiles of targets 1, 5 and 9 are shown in Figure 10, and the quality parameters of imaging result are listed in Table 2.

Two-dimensional profiles of point targets 1, 5 and 9(from left to right). (a) Azimuth profile of point target 1, (b) azimuth profile of point target 5, (c) azimuth profile of point target 9, (d) range profile of point target 1, (e) range profile of point target 5, (f) range profile of point target 9.

Full size image

As shown in Figure 10, the profiles of targets 1, 5 and 9 are very close to those of the ideal results. This implies that these three targets are focused well with the proposed algorithm. It can also be found from Figure 10 (a),(c) that the first side lobes are slightly asymmetric for the two edge points. The reason for this is the approximation error shown in (58).

The resolutions of bistatic SAR are 2-D dependent [44]. Therefore, to demonstrate the performance of the proposed algorithm, the resolution before geometric correction is presented in Table 2. According to parameters listed in Table 1, the theoretical value of azimuth resolution and the slant-range resolution is 5.31 m. The theoretical value of peak side lobe ratio (PSLR) is -13.26 dB, and the theoretical value of integrated sidelobe ratio (ISLR) is -9.72 dB [40].

In Table 2, we can find that both the azimuth resolution and the range resolution have a deviation less than 0.08 m with respect to the theoretical value. The maximum deviation of the measured azimuth PSLR is less than 0.49 dB of the theoretical value. The maximum deviation of the measured range PSLR is less than 0.14 dB of the theoretical value. For ISLR, the maximum deviation is less than 0.48 and 0.65 dB along the azimuth and the range directions, respectively.

7.3 Imaging simulation for large scene

According to (60), to avoid the Doppler ambiguity, the Doppler frequency should meet max(|f _a|) < PRF/2. In other words, the azimuth width of the imaged scene should satisfy $W_{\mathrm{a}}<\frac{\lambda r_{0\mathrm{d}}}{v}·\left|\mathrm{PRF}-B_{\mathrm{a}}\right|$. Based on the parameters listed in Table 1, we get W _a < 4.86 km. If we want to extend the azimuth width of the target scene, higher PRF should be applied.

If the scene size extends beyond the constraint condition (59), small blocks should be used in the imaging process. The block process along the range direction is simple, since only the reference slant range r ₀ is replaced for different range blocks. However, as stated before, DCC has to be applied for different azimuth blocks. Assuming that the azimuth coordinate of a block is 1,000 m, the Doppler centroid of this block is 378.53Hz. Then, before and after DCC, the 2-D spectrum of the point target located at (0,1,000 m) is shown in Figure 11.

The 2-D spectrum of the point target located at the block centre. (a) Before DCC and (b) after DCC.

Full size image

It can be seen from Figure 11 that the Doppler frequency is shifted to baseband using DCC. Therefore, the same imaging process can be applied for this block.

Then, an imaged scene with size 8 km × 4 km (X-axis × Y-axis) is assumed. According to the constraint condition (59), we divide the whole scene into 3 × 5 (X-axis × Y-axis) blocks with the consideration of small overlaps between blocks. Three point targets with coordinates are located in this scene. The coordinates can be listed as (-4,000,-2,000), (0,0) and (4,000,2,000). Using the proposed algorithm to process the data, the final focused image and zoomed contour graphs of these three targets are shown in Figure 12.

Simulation result for large scene using the proposed algorithm.

Full size image

It can be seen from Figure 12 that three targets are well focused and located at correct position using the proposed algorithm. This demonstrates the capability of the proposed algorithm for large scene imaging.

Bistatic spaceborne/stratospheric SAR has the potential to play an important role in the future missions. The crucial steps for such a system are synchronization and imaging. There have been quite a number of published studies in these two fields. However, the combined research of synchronization and imaging has not been developed. This paper proposes an integrative synchronization and imaging approach for a particular configuration-bistatic spaceborne/stratospheric SAR with a fixed receiver.

In published methods, the direct-path signal is a common choice for the synchronization of bistatic SAR systems, since the direct-path signal takes the advantage of high SNR and clean phases. However, it is still very difficult to extract synchronization errors with high precision from the direct-path signal, if no other auxiliary data is variable. In this paper, as a novel idea, time delays and peak phases of direct-path signal are utilized to compensate corresponding components of reflected signal directly. Time and phase synchronization errors can be completely removed using the presented method. Meanwhile, as the cost of synchronization, the system impulse response of the synchronized echo becomes quite different from that of the general bistatic SAR. To focus the particular synchronized data, the 2-D spectrum of synchronized data is derived under linear approximations, and then a frequency-domain imaging algorithm is proposed. The theoretical analysis and simulation results show that the proposed approach can provide accurate time and phase synchronization and well-focused images.

Based on the presented method, we continue our research on bistatic spaceborne/stratospheric SAR. The next steps of research are multiple. One is the study of synchronization and imaging for general bistatic spaceborne/stratospheric SAR configuration. The other is the study of the interferometric application of such a bistatic system.

The authors express their appreciation to Dr. Zeng Zhanfan. Thanks for his proofreading and helpful discussions. We also would like to appreciate the reviewers for their very patient work and kind suggestions, which helped a lot in correcting many mistakes, increasing the readability, and improving the quality of this paper.

Authors and Affiliations

1. College of Electronic Science and Engineering, National University of Defence Technology, Changsha, Hunan, 410073, People’s Republic of China

   Qilei Zhang, Wenge Chang & Xiangyang Li

Corresponding author

Correspondence to Qilei Zhang.

Competing interests

The authors declare that they have no competing interests.

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions