sensors Article

Generalized Chirp Scaling Combined with Baseband Azimuth Scaling Algorithm for Large Bandwidth Sliding Spotlight SAR Imaging Tianzhu Yi 1, *, Zhihua He 1 , Feng He 1 , Zhen Dong 1 and Manqing Wu 2 1

2

*

School of Electronic Science and Engineering, National University of Defense Technology, Sanyi Avenue, Changsha 410073, China; [email protected] (Z.H.); [email protected] (F.H.); [email protected] (Z.D.) China Electronics Technology Group Corporation (CETC), China Academy of Electronics and Information Technology, Beijing 100846, China; [email protected] Correspondence: [email protected]; Tel.: +86-187-7318-5317

Academic Editor: Vittorio M.N. Passaro Received: 13 February 2017; Accepted: 26 May 2017; Published: 29 May 2017

Abstract: This paper presents an efficient and precise imaging algorithm for the large bandwidth sliding spotlight synthetic aperture radar (SAR). The existing sub-aperture processing method based on the baseband azimuth scaling (BAS) algorithm cannot cope with the high order phase coupling along the range and azimuth dimensions. This coupling problem causes defocusing along the range and azimuth dimensions. This paper proposes a generalized chirp scaling (GCS)-BAS processing algorithm, which is based on the GCS algorithm. It successfully mitigates the deep focus along the range dimension of a sub-aperture of the large bandwidth sliding spotlight SAR, as well as high order phase coupling along the range and azimuth dimensions. Additionally, the azimuth focusing can be achieved by this azimuth scaling method. Simulation results demonstrate the ability of the GCS-BAS algorithm to process the large bandwidth sliding spotlight SAR data. It is proven that great improvements of the focus depth and imaging accuracy are obtained via the GCS-BAS algorithm. Keywords: generalized chirp scaling; baseband azimuth scaling; large bandwidth; sliding spotlight

1. Introduction High resolution acquisition has been a hot topic in the field of synthetic aperture radar (SAR) systems and signal processing [1]. In addition, considerable studies have been performed on high-resolution SAR involved in both the large bandwidth [2–6] and sliding spotlight modes [1]. The conventional baseband azimuth scaling (BAS) algorithm [7] can work well with a sliding spotlight SAR signal with a carrier frequency that is much larger than the bandwidth. However, this algorithm is not suitable for the large bandwidth sliding spotlight SAR signals, because the approximation of the algorithm for sub-aperture data cannot decouple the sub-aperture signal in the large bandwidth mode. Targets deviating from the reference slant range defocus due to their inability to decouple the two dimensions’ high order phase binding caused by the large bandwidth mode. Additionally, the high resolution enabled by the large bandwidth sliding spotlight SAR is attracting increasing attention. Reference [8] first proposed the processing principle of the chirp scaling (CS) algorithm to solve the SAR signal in the large bandwidth mode and large squint angle mode. Zaugg of American Brigham Young University subsequently developed this principle and proposed a generalized high-order CS algorithm [9] that could compensate for the defocusing caused by the large bandwidth mode. At present, the generalized chirp scaling (GCS) algorithm has been a generalized signal processing method for SAR work in both large bandwidth and large squint modes.

Sensors 2017, 17, 1237; doi:10.3390/s17061237

www.mdpi.com/journal/sensors

Sensors 2017, 17, 1237

2 of 19

The sliding spotlight SAR appropriately settles the conflicts between the high resolution and the large mapping range due to the fact that it can flexibly control the beam direction of the antenna [10]. For this reason, it has attracted many researchers to study sliding spotlight SAR more thoroughly and intensively. The requirements for data with high quality and accuracy have increased in the field of remote sensing in recent years. The higher accuracy imaging algorithms of high-resolution SAR are required by the increasing demands of high-resolution remote sensing data. Among sliding spotlight SAR imaging algorithms, many well-developed methods exist, such as the up-sampling imaging method, full-aperture processing method, and BAS algorithms. The up-sampling imaging method requires a large amount of data due to zero-padding processing [11]. Typically it acquires Ba_total /PRF (where Ba_total is the azimuth total bandwidth of the sliding spotlight SAR, and PRF is the pulse repetition frequency) times the data after zero-padding. The ω − k algorithm, which is introduced in the full-aperture processing method to solve the coupling problems along the range and azimuth dimensions, results in the algorithm time consumption [12]. Alberto Moreira of the German Aerospace Center (DLR) proposed an extended chirp scaling (ECS) algorithm that can solve the imaging problem of the sliding spotlight SAR in the antenna scanning modes [13]. Since then, the ECS algorithm has been the basic framework for multi-mode SAR imaging methods [13]. Pau Prats later developed the ECS algorithm and presented the BAS algorithm [7] to solve the imaging problem in the spotlight, Terrain Observation by Progressive Scans (TOPS), and sliding spotlight modes. In this method, the conventional CS algorithm is used to process the sub-aperture data of sliding spotlight SAR to eliminate the coupling between the range and azimuth dimensions successfully. After decoupling in the BAS algorithm, the Doppler rate of these sub-aperture data is adjusted again. After recombination of all sub-apertures, a de-rotation function and an azimuth compressed function are used to focus the sliding spotlight SAR data. Moreover, the experimental results of Terra SAR in the sliding spotlight and TOPS modes also prove the validity of the BAS algorithm. With the increasing requirements for data accuracy in the remote sensing field, studies on the feasibility of high-resolution spaceborne/airborne SAR systems and imaging algorithms for high-resolution SAR are gaining more attention. However, the conditions of large bandwidth and low carrier frequency are not taken into consideration according to the BAS algorithm [14]. In other words, when Br /(2 f 0 ) > 0.04, the condition: Br (1) f0  2 is invalid, and the second-order approximation model of the BAS algorithm cannot address the high-order term, which can cause large phase errors and defocus the SAR images under the large bandwidth condition. The impact of the large bandwidth is involved in this paper in detail. In addition, taking the ideas of [8,9] about the processing of the large bandwidth SAR into account, this paper presents the method of the GCS-BAS algorithm for large bandwidth sliding spotlight SAR imaging. The proposed algorithm can improve the focus depth and accuracy of the processing results of the BAS algorithm for high-resolution signals from the large bandwidth sliding spotlight SAR. The remainder of this paper is organized as follows. Section 2 analyzes the influence of the second-order approximation model in the conventional BAS algorithm on the phase of the signal from the large bandwidth sliding spotlight SAR. Section 3 presents the algorithm framework of the proposed GCS-BAS imaging algorithm and deduces the phase expressions and phase filters in each step of the proposed algorithm in detail. Finally, Section 4 compares the imaged results of the simulation by the conventional BAS algorithm, the back projection (BP) algorithm, and the GCS-BAS algorithm proposed in this paper. 2. Analysis of Model Error From Figure 1, the echo of the target in ( X 0 , R0 ) demodulated to the baseband can be illustrated as:

Sensors 2017, 17, 1237

3 of 19

  ssliding_spotlight (η, τ ) = Wr (τ − τ0 ) · Wa (η − ηc ) · exp jπKr (τ − τ0 )2 ·   (Vg /Vr ) x − X 0 exp(− j2π f 0 τ0 ) · rect La

(2)

where τ is the fast time and η is the slow time. ηc represents the moment of the Doppler centroid. Wr (·) and Wa (·) are the antenna patterns along the range and azimuth dimensions, respectively. Tr is the pulse width. Kr is the chirp rate of the transmit signal. Vg and Vr represent the velocities of beam 2R(η )

steering on the ground and the radar, respectively. τ0 = c represents the delay of the target. R(η ) q is the range history of the target, expressed as R(η ) = R20 + Vr2 η 2 . R0 is the nearest range between the radar and the target. c is the speed of light. x is the trajectory of the radar. L a is the length of the synthetic aperture. rect(·) represents the illumination area of the beam [15]. A two dimensional (2D) fast Fourier transform (FFT) yields a signal in the 2D frequency domain: θ2D f τ , f η r


where D f η , Vr =

1−




c2 f η2 ,f 4Vr2 f 02 η

4πR0 f 0 =− c

s


fτ fτ π f τ2 D f η , Vr + 2 + 2 − f0 Kr f0

(3)

and f τ are the azimuth frequency and range frequency, respectively.

Br 2 f0

On the condition of > 0.04 [9], the conventional CS is only able to compensate the high-order phase error of the reference range. However, it cannot resolve the high-order phase error caused by the large bandwidth that varies along the range. This problem leads to reduced precision for the secondary range compression (SRC) and range cell migration correction (RCMC). This degradation causes sub-aperture focusing performance that results in the deterioration of azimuth focusing later in the operation of azimuth scaling of the BAS. Therefore, the high-order phase error variation along the range should be taken into account to separate the coupling between the range and azimuth dimensions. To illustrate the influence of the order-2 model approximation, a simulation is performed with the parameters listed in Table 1. Table 1. Parameters for the model phase error. Parameter

Value

Carrier frequency/ f 0 Bandwidth/Br Reference range/Rre f Pulse repetition frequency/PRF Velocity of radar/Vr Pulse width/Tr

8 GHz 1 GHz 30 km 600 Hz 240 m/s 40 µs

Figure 2a–d illustrates the phase error for different order approximations. The horizontal and vertical axes represent the range frequency and range, respectively. The color bar represents the magnitude of the phase error in different colors. Figure 2a illustrates that the phase error of the order-2 approximation is small for bandwidths between 400 and 600 MHz. This phase error can be neglected in the process of the sliding spotlight SAR. When the bandwidth is larger than 600 MHz, the phase error becomes large regardless of the reference range or edge range. The impact of the phase error is discussed in the previous paragraph. Figure 2b–d show the phase error of different order model approximations. The figures demonstrate that the phase error becomes lower as the order of the model increases. In Figure 2d, the phase error is reduced to 1 degree; this phase error is considered to have a small impact on the range compression and RCMC. High-order phase coupling should be considered in the processing of the large bandwidth sliding spotlight SAR. Each processing step for the large bandwidth sliding spotlight SAR is discussed below.

neglected processin ofthe the previous sliding spotlight SAR.Figure When the bandwidth largererror than 600 MHz, phase error in is the discussed paragraph. 2b–d show theisphase of different the phase becomes largeThe regardless the reference range or edgeerror range. The impact thethe order modelerror approximations. figures of demonstrate that the phase becomes lowerofas phase error is discussed in the previous paragraph. Figure 2b–d show the phase error of different order of the model increases. In Figure 2d, the phase error is reduced to 1 degree; this phase error is order model approximations. The figures demonstrate that the phase error becomes lower as the considered to have a small impact on the range compression and RCMC. High-order phase order of the model increases. In Figure 2d, the phase error is reduced to 1 degree; this phase error is should be considered in the processing of the large bandwidth sliding spotlight SAR. Each4 of 19 Sensorscoupling 2017, 17, 1237 considered to have a small impact on the range compression and RCMC. High-order phase processing step for the large bandwidth sliding spotlight SAR is discussed below. coupling should be considered in the processing of the large bandwidth sliding spotlight SAR. Each processing step for the large bandwidth sliding spotlight SAR is discussed below. Vr

fη Vr

fη

ψS ψS

ψE

Binst

ψE

Binst

PRF / 2 PRF / 2

BT

( X ' , R0 )

BT

'

( X , R0 )

η η PRF

− PRF / 2

PRF

− PRF / 2 TB

TB

(a) (a)

(b) (b)

Figure 1. Schematic of the acquisition geometry for sliding spotlight spotlight synthetic aperture

Figure Figure 1. Schematic of theofacquisition geometry for for sliding spotlight spotlight synthetic aperture aperture radar 1. Schematic the acquisition geometry sliding spotlight spotlight radar (SAR) (a) and time frequency distribution (TFD) of sliding spotlight SARsynthetic data (b). ψ S and (SAR) (a) and(SAR) time (a) frequency (TFD) of (TFD) slidingofspotlight SAR data and represent ψ SψEand radar and timedistribution frequency distribution sliding spotlight SAR(b). dataψS(b). ψ E represent the start and end scanning angles, respectively. PRF is the pulse repetition frequency. the startψand end scanning angles, respectively. PRF is the pulse repetition frequency. B represents represent the start and end scanning angles, respectively. PRF is the pulse repetition frequency. inst Binst E represents the instant azimuth bandwidth [16], and BT is the additional band-width the instant and bandwidth BT is the additional band-width resulting from the beam Binst azimuth representsbandwidth the instant [16], azimuth [16], and B is the additional band-width T resulting from the beam steering [16]. steeringresulting [16]. from the beam steering [16].

26 27

Phase error/ [deg] Phase error/ [deg] 25

25 25

1000 1000

26

800 800

27

600 600

2828 Range/ km Range/ km

200 200

3030

00

3131

-200 -200

3232

28 28

6060

29 29

5050

30 30

4040

31 31

-800 -800

3535 -600-400 -400-200 -200 0 0 200 200 400 400 -600 RangeFrequency Frequency/M /MHz Hz Range Sensors 2017, 17, 1237 Phase error/ [deg]

1010

35 35 -600 -600 -400 -400 -200 -200 00 200 200 400 400 Range RangeFrequency Frequency/M /MHzHz

0 0

Phase error/ [deg] 25 Figure 2. Cont. 6 2. Figure Cont.

0.45

26

4

27

0.4

27

0.35

28

2

30

0

31 32

-2

33 -4

34 35 -600 -400 -200 0 200 400 Range Frequency /M Hz

(c)

Range/ km

28 29

5 of 20

(b) (b)

26

Range/ km

2020

34 34

(a) (a)

25

3030

33 33

-600 -600

3434

70 70

32 32

-400 -400

3333

80 80

27 27

400 400

2929

Phase error/ [deg] Phase error/ [deg]

26 26

Range/ Range/km km

25

0.3

29 30

0.25

31

0.2

32

0.15

33

0.1

34

0.05

35 -600 -400 -200 0 200 400 Range Frequency /M Hz

0

(d)

Figure 2. (a–d) shows the phase error of order-2–order-5 approximations, respectively. The color-bar

Figure 2. (a–d) shows the phase error of order-2–order-5 approximations, respectively. The color-bar represents the magnitude of the phase error by different colors and the horizontal and vertical axes represents the magnitude of the phase error by different colors and the horizontal and vertical axes represent the range frequency and range, respectively. represent the range frequency and range, respectively. 3. GCS-BAS 3.1. Procedures of GCS-BAS The blocks of processing procedures for GCS-BAS are shown in Figure 3, while the following text and Appendix A derive the signal of each step for GCS-BAS.

35 -600 -400 -200 0 200 400 Range Frequency /M Hz

35 -600 -400 -200 0 200 400 Range Frequency /M Hz

(c)

0

(d)

Figure 2. (a–d) shows the phase error of order-2–order-5 approximations, respectively. The color-bar 5 of 19 represents the magnitude of the phase error by different colors and the horizontal and vertical axes represent the range frequency and range, respectively.

Sensors 2017, 17, 1237

3. GCS-BAS 3. GCS-BAS 3.1. Procedures of GCS-BAS 3.1. Procedures of GCS-BAS The blocks of processing procedures for GCS-BAS are shown in Figure 3, while the following text The blocks processing procedures for for GCS-BAS are shown in Figure 3, while the following and Appendix A of derive the signal of each step GCS-BAS. text and Appendix A derive the signal of each step for GCS-BAS.

H5

H1

H2

H3

H6

H7

H8

H4

Figure 3. Blocks of proposed generalized chirp scaling-baseband azimuth scaling (GCS-BAS) for the Figure 3. Blocks of proposed generalized chirp scaling-baseband azimuth scaling (GCS-BAS) for the large bandwidth sliding spotlight SAR image. large bandwidth sliding spotlight SAR image.

Reference [7] explains how to divide the raw data into the processed sub-apertures’ data. The steps are as follows: Step 1: Considering the motion error of the platform, we adopt the method presented by [17] to eliminate the impact of motion error (readers can refer to [17]). This method has the advantage of eliminating the motion error on the basis of the platform trajectory measured by a precise device like global position system (GPS) or inertial navigation system (INS). The precise trajectory of the platform is easy to obtain in the simulation. Then the 1st and 2nd motion errors are compensated by the parameters of the 3D motion error which are extracted from the GPS/INS data [17]. Step 2: FFT is implemented to transfer the sliding spotlight SAR sub-aperture data to the 2D frequency domain. Then, the sliding spotlight SAR sub-aperture data in the 2D frequency domain are multiplied by the high-order compensation function; Step 3: After high-order compensation, the sub-aperture data are transformed by inverse FFT (IFFT) along the range, and are then multiplied by a high-order CS function; Step 4: After high-order CS, the sub-aperture data are transferred to the 2D frequency domain by FFT along the range. Then, the sub-aperture data are processed along the range with range compression (RC), SRC, and RCMC; Step 5: The sub-aperture data are transferred into the range-Doppler domain by IFFT along the range. The phase error caused by the motion error is corrected by the residual error compensation

Sensors 2017, 17, 1237

6 of 19

which is calculated from the GPS/INS data [17]. It can be thought that the impact of the motion error is eliminated after the processing of the residual error correction [17]. Then, the data are multiplied by the azimuth compression function to achieve sub-aperture data imaging. Steps 6–9: These steps of the BAS are described in [7] in detail. Functions of H1 ∼ H8 are explained in detail in the following text. 3.2. Theoretical Formulation The following formulation is simplified with the form of the phase for the GCS-BAS algorithm which primarily processes phase functions. Equation (3) is expanded by a high-order Taylor series expansion. θ2D


= − 4πRc 0 f0 D f η − n

4πR0 f cD ( f η ) τ

−

2 2 4πR0 f 0 (−1+ D ( f η )) f τ c 2D3 ( f η ) f 02

−

π f τ2 Kr

− 4πRc 0 f0 ∑ γi f τi

(4)

i =3

where R0 is the slant range vector, and the coefficients γi are discussed in [9]. We filter the sub-aperture raw data of the sliding spotlight and the large bandwidth SAR in the 2D frequency domain with the function: "  # n 
4Rre f f 0 γi H1 f τ , f η = exp jπ ∑ Xi + f τi (5) c i =3 where Xi is a high order turbulent item, which can simplify the GCS algorithm expressions, and its
specific equation is discussed in Appendix A. The second item in H1 f τ , f η is the compensation of a high order phase in the reference slant range that can minimize the range variation of a high order phase [18]. After high order phase compensation, the phase in the 2D frequency domain can be
expressed by θ1 f τ , f η :

θ1 f τ , f η




    n 4 R0 − Rre f f 0
4πR0 f 0 π f τ2 4πR0
fτ − =− D fη − + π ∑  Xi − γi  f τi c Km c cD f η i =3

(6)

Then, IFFT is performed along the range for the sub-aperture SAR data after high order phase compensation. The solution of a stationary phase point can be calculated based on the principle of the stationary phase (POSP), f τ, POSP = Km (τ − τd ). The processed sub-aperture data phase in the Range
Doppler (RD) domain can be expressed by θ2 τ, f η : θ2 τ, f η





= − 4πRc 0 f0 D f η + πKm (τ − τd )2 + n 
 i π ∑ Xi − 2 f 0 D f η γi ∆τ Km (τ − τd )i

(7)

i =3

Differential can be achieved by the method shown in Figure 4. The relationships between  RCMC  τs , τre f , τd , ∆τ f ηre f , and ∆τ are explained in detail in Appendix A. The method for deducing the high order CS function is also discussed in Appendix A. The high order CS function [19] can be expressed
by H2 τ, f η : " #  2  i n
H2 τ, f η = exp jπq2 τ − τre f + jπ ∑ qi τ − τre f (8) i =3

where the expressions of qi (i ≥ 2) are explained in Appendix A. This step mainly uses the high order CS function to compensate for the influences of high order range cell migration (RCM) caused by the large bandwidth [20]. The phase of the sub-aperture data after processing high orders can be expressed
by θ2 τ, f η :

Sensors 2017, 17, 1237

θ2 τ, f η




7 of 19

 2  i n
= − 4πRc 0 f0 D f η + πKm (τ − τd )2 + πq2 τ − τre f + π ∑ qi τ − τre f n

i =3

+π ∑ Xi − 2 f 0 D f η γi ∆τ 




i =3



i (τ Km

− τd )

i

(9)


Also θ2 τ, f η can be expressed by the other form for the relationships between the coefficients qi (i ≥ 2) and Xi (i ≥ 3): θ2 τ, f η




n
= − 4πRc 0 f0 D f η + πC0 (∆τ ) + ∑ Ci (τ − τs )i i =1
K 4πR0 f 0 ≈ − c D f η + πC0 (∆τ ) + π αf (τ − τs )2   n + ∑ qi + K if Xi (τ − τs )i

(10)

i =3

The definitions of C0 (∆τ ), α, Ci , and K f are provided in Appendix A. FFT is performed along the range dimension for the sub-aperture data with the phase of Equation (10) in the RD domain. Neglecting the impact of high order items (3 or higher) on the stationary phase point [12], we can obtain the solution of the stationary phase point by the POSP, τPOSP = Kα f τ + τs . The phase after FFT f
can be expressed by θ3 f τ , f η : θ3 f τ , f η





≈ − 4πRc 0 f0 D f η + πC0 (∆τ ) − π Kαf f τ2 − 2π f τ τs  i n  +π ∑ qi + K if Xi Kαf f τi i =3
= − 4πRc 0 f0 D f η + πC0 (∆τ ) − 2πτre f (1 − α) f τ − 2πατd f τ  i n  −π Kαf f τ2 + π ∑ qi + K if Xi Kαf f τi

(11)

i =3

In Equation (11), the first term is the phase item related to the azimuth compression; the second term is the so-called residual video phase (RVP) [21]; the third term is the bulk RCM; the forth term is the linear phase, which is related to the target’s position along the range; the fifth and sixth terms are
the phase that influence the accuracy of RC. H3 f τ , f η is multiplied to achieve the RC, secondary RC and bulk RCMC.  !i      n 
α α 2 H3 f τ , f η = exp j2πτre f (1 − α) f τ + π f τ − π ∑ qi + K if Xi f τi  (12)   Kf Kf i =3

θ4


Filtered by H3 f τ , f η , the phase of the signal in the 2D frequency domain can be expressed as
fτ , fη :

4πR0 f 0 θ4 f τ , f η = − D f η + πC0 (∆τ ) − 2πατd f τ (13) c

The second term in Equation (13) is the RVP that must be eliminated in the RD domain. A range IFFT using POSP leads to the following phase:

4πR0 f 0 θ4 τ, f η = − D f η + πC0 (∆τ ) c

(14)

Thus, we can obtain the azimuth compression function based on Equation (14): H4 τ, f η









 4πR0 f 0  = exp jπ D f η − 1 − πC0 (∆τ ) c

Then, the signal in the RD domain can be expressed as:

 (15)

Sensors 2017, 17, 1237

8 of 19

     

4πR0 f 0 2R0 2R S4 τ, f η = ωr τ − · Wa f ηc · sinc τ − 0 · exp − j c c c

(16)


After being processed by H4 τ, f η , the phase coupling between the range and azimuth caused by the large bandwidth can be eliminated. If the coupling between the range and azimuth dimensions is not eliminated clearly, it leads to the unmatched Doppler rate in the following BAS algorithm’s azimuth scaling. This also results in defocusing in the azimuth dimension. The following operations are similar to the BAS algorithm’s azimuth scaling and hence are described only briefly here. Sensors 17, 1237 compression for the sub-aperture data, a purely quadratic phase9 of 20 After2017, azimuth shape is Sensors 2017, 17, 1237 9 of 20 added using:   Kscl ( R0 ) is the Doppler rate the range vector of
factor, which is related
to the velocity π Vr and 2 Kscl ( R0 ) is the H V Doppler rate factor, which is related to the velocity and the range vector of (17) τ, f = exp − j2π f η R exp − j f ( ) r η η η v 5 0 Rscl ( R0 ) . Kscl ( R0 ) Rscl ( R0 ) .

Kscl ( R0 ) is the Doppler rate factor, which is related 2Vr2 2 to the velocity Vr and the range vector of K scl ( R0 ) = − (18) 2V Rscl ( R0 ). K scl ( R0 ) = −λ Rscl ( rR0 ) (18) λ Rscl2V ( Rr20 ) =− 0 ) Equation scl ( Rin The factor is different from theKone causes the extension along the (18) λRscl(15) ( R0 )that The factor is different from the one in Equation (15) that causes the extension along the azimuth dimension after azimuth scaling. The time shift ηv ( R0 ) with a linear phase ramp is azimuth dimension after azimuth scaling. The time a linear phase is The factor is different from the one in Equation (15)shift that ηcauses the extension along ramp the azimuth v ( R0 ) with introduced to minimize the extension of azimuth dimension, whose expression is given in [15]. An introduced toazimuth minimizescaling. the extension of azimuth dimension, whose expression is given in [15]. An dimension after The time shift η R with a linear phase ramp is introduced to ( ) 0 azimuth IFFT is performed on the sub-aperturev data after azimuth scaling. The shift and azimuth IFFT is performed on the sub-aperture data after azimuth scaling. The shift and minimize the extension of azimuth whose expression isare given in [15].byAnH azimuth IFFT is recombination operations need todimension, be done while the sub-apertures processed 1 ~ H 5 . The recombination operations need to be done while the sub-apertures areand processed by H1 ~ Hoperations 5 . The performed on the sub-aperture data after azimuth scaling. The shift recombination shift and recombination operations in Figure 5, as well as the meaning of variances, are discussed shift and recombination operations in Figure 5, as well as variances, arerecombination discussed needwell to be done while the sub-apertures are processed bythe H1meaning ∼ H5 . of The shift and in [22]. well in [22]. operations in Figure 5, as well as the meaning of variances, are discussed well in [22]. fη fη

τ ref τ ref

τs τs

τd τd

Δτ Δτ

(( ))

Δτ fηref Δτ fηref

fηref fηref

R0 R0

Rref Rref

r r

Figure 4. 4.Chirp differentialrange range cell migration correction. Figure Chirpscaling scaling for for differential cell migration correction. Figure 4. Chirp scaling for differential range cell migration correction.

fη fη

Δηi'' '' Δηi

Δηi' ' Δηi

ηis s'' ηis ' ηi e'' ηi '' ηis ' ηi '' ηis s Δηi ηie ηi Δηi ηie e

ηiee' ηi '

η η

Figure time-frequencydistribution distribution (TFD). Figure5.5.Variation Variation of of time-frequency (TFD). Figure 5. Variation of time-frequency distribution (TFD).

However, the total azimuth bandwidth still exceeds the PRF after all sub-apertures are However, the total azimuth bandwidth still exceeds the PRF after all sub-apertures are assembled. Therefore a demodulation operation similar to the two-step algorithm [23] is carried out assembled. Therefore a demodulation operation similar to the two-step algorithm [23] is carried out to move the signal to the azimuth baseband using the demodulation function: to move the signal to the azimuth baseband using the demodulation function: H 6 (τ ,η ) = exp ( − jπ K rot ( R0 )η 2 2) (19) H (τ ,η ) = exp ( − jπ K ( R )η ) (19)

Sensors 2017, 17, 1237

9 of 19

However, the total azimuth bandwidth still exceeds the PRF after all sub-apertures are assembled. Therefore a demodulation operation similar to the two-step algorithm [23] is carried out to move the signal to the azimuth baseband using the demodulation function: Sensors 2017, 17, 1237

10 of 20

  H6 (τ, η ) = exp − jπKrot ( R0 )η 2

(19)

2Vr22 K rot ( R0 ) = − 2V Krot ( R0 ) = − λ Rrot (rR0 ) λRrot ( R0 )

(20) (20)

R00)) ininEquation Equation (18) (18) and and RRrot Equation(20) (20)are aredescribed describedin indetail detailin in [7]. [7]. The The TFD TFD ( RR00)) ininEquation scl ( rot((R RRscl after azimuth derotation derotation is is shown shown in in Figure Figure 6. 6. Time-frequency distribution 800

Azimuth frequency/Hz

600 400 200 0 -200 -400 -600 -800

-4

-3

-2

-1 0 1 Azimuth time/s

2

3

4

Figure 6. TFD of SAR data after azimuth derotation.

After filtered filteredbybyH (Hτ,6 (ητ),η sliding spotlight SAR azimuth time sampling and image ) , the After , the sliding spotlight SAR azimuth time sampling and image sampling 6 sampling with the following [13]: vary with vary the following equations equations [13]:   R  1 −Rscl0scl 0  1⋅ − rot0  PRF0 PRF  RR (21) rot R 0 ( R ) (21) ∆x f inal = ∆xorig 1 − R scl ( R0 )  R rot ( R )0  Δx final = Δxorig  1 − scl 0  )  sampling, and ∆xorig is provided rot ( R0image where PRF0 is the new azimuth sampling, ∆x f inal is theRfinal in [13], and the effective chirp rate has been changed after azimuth scaling and azimuth derotation: where PRF0 is the new azimuth sampling, Δx final is the final image sampling, and Δxorig is 1 1 1 =1 PRF0 = PRF ·

provided in [13], and the effectiveKchirp rate has been changed after azimuth scaling and azimuth (22) e f f ( R0 ) = Kscl ( R0 ) − Krot ( R0 ) derotation: An azimuth FFT results in the sliding spotlight data in the RD domain. The azimuth match filter K eff ( R0 ) = K scl ( R0 ) − K rot ( R0 ) (22) is as following: ! data in the RD domain. The azimuth match An azimuth FFT results in the sliding spotlight

π PRF 2 filter is as following: H7 τ, f η = Wa f η · exp j f η , − PRF + f dc (23) + f dc < f η < 2 K e f f ( R0 ) 2   π PRF PRF + f dc < fη < +f
H 7 (τ , fη ) = Wa ( fη ) ⋅ exp  j fη2  , − (23)  of  where f dc is the mean Doppler centroid the whole data acquisition, and Wa f ηdc is a weighting K R 2 2 ( ) eff 0   function which can suppress the side-lobe. An IFFT on the data after being filtered by the function
is thethe mean Doppler of preservation the whole data acquisition, and should Wa ( fη )beismultiplied a weighting where H7 τ, f η f dcyields focused data.centroid For phase purposes, the data by the function: function which can suppress the side-lobe. An IFFT onthe data  after!being filtered by the function rscl0 2 2 exp jπK η the data should be multiplied (24) ) = phase H 7 (τ , fη ) yields the focusedHdata. preservation t ( R0 ) 1 −purposes, 8 ( τ, ηFor rrot0 by the function: 2    r  H 8 (τ ,η ) = exp  jπ K t ( R0 )  1 − scl 0  η 2     rrot 0   

where Kt ( R0 ) = −

(24)

2Vr2 . The proposed function achieves the phase preservation λ ( rrot ( R0 ) − rscl ( R0 ) )

Sensors 2017, 17, 1237 Sensors 2017, 17, 1237

10 of 19 11 of 20 2V 2

Time-frequency distribution where Kt ( R0 ) = − λ(r ( R )−r r ( R )) . The proposed function achieves the phase preservation efficiently rot 0 0 scl 800 Sensors 2017, 17, 1237 11 of 20 without interpolation. The TFD after phase preservation is shown in Figure 7. Time-frequency distribution

400800 600

200 Azimuth frequency/Hz

Azimuth frequency/Hz

600

400

0

200

-200

0

-400

-200

-600

-400

-800-600 -5

0 Azimuth time/s

-800 -5

5

0 Azimuth time/s

5

Figure 7. TFD of SAR data after phase preservation. Figure TFDof ofSAR SAR data after Figure 7. 7.TFD afterphase phasepreservation. preservation.

4. Experimental Results and Discussion

4. Experimental Results and Discussion 4. Experimental Results and Discussion In this section, all experiments are carried out by a PC with an i7 Dual-Core CPU of 3.3 GHz Inmemory. this section, allverify experiments are carried outGCS-BAS by a PC with an i7 Dual-Core CPU of 3.3 GHz and 32 GB To the validity of out the algorithm, a simulation performed In this section, all experiments are carried by a PC with an i7 Dual-Core CPU ofis3.3 GHz and and 32 GB memory. To verify the validity of the GCS-BAS algorithm, a simulation is performed with the parameters listed in Table 2. 32 GB memory. To verify the validity of the GCS-BAS algorithm, a simulation is performed with the with the parameters listed in Table 2.

parameters listed in Table 2.

Table 2. Parameters for sliding spotlight SAR. Table 2. Parameters for sliding spotlight SAR.

Table 2. Parameter Parameters for sliding spotlight SAR. Value Parameter Value Carrier frequency/ GHz Carrier frequency/ ff00 88Value GHz Parameter

Bandwidth/ Bandwidth/BBr r Carrier frequency/ f0 Scaling range/ Scaling range/RRrrefref Bandwidth/B Scaling range/R Pulse repetition frequency/ re f PRF PRF Pulse repetition frequency/ Pulse repetition frequency/PRF V Velocity of platform/ Velocity of platform/ Vrr Velocity of platform/Vr Pulse width/TTr Pulse Pulsewidth/ width/Tr r Imageswath swath Image Image swath Factor of sliding spotlight/ M Factor of sliding spotlight/M Factor of sliding spotlight/ M

GHz 118GHz GHz

30 km 30 km 1 GHz 30Hz km 600 600 Hz 600m/s Hz 240 240 m/s 240 m/s 40 μs 40 40 μs µs 55km × 5 km km × km 5 km × 55 km 0.4 0.4 0.4

The results of the conventional BAS, GCS-BAS, BP, and full-aperture imaging algorithms are

The results of the conventional conventional BAS, GCS-BAS, BP, and full-aperture full-aperture imaging algorithms are The results the GCS-BAS, BP, and imaging algorithms presented toof compare the focusingBAS, performance of the conventional BAS algorithm and GCS-BAS are presented to compare compare theoffocusing focusing performance of the thein conventional BAS algorithm algorithm and and GCS-BAS GCS-BAS presented to the performance of conventional algorithm. The results the four algorithms are shown Figures 8–11.BAS algorithm.The Theresults resultsof ofthe thefour fouralgorithms algorithmsare areshown shownin inFigures Figures 8–11. 8–11. algorithm.

0 -1 -2 -3 -2

Range Slice

0

0 -1 -2 -3 -2

-1

-1

0 Range / m

(a)

0 Range / m

(a)

1

1

2

2

-10

-10

-20

-10

0

-10

-20

-20 -30

-20-30 -30

Azimuth Slice Magnitude/ dB

1

Azimuth Slice 0

Magnitude/ dB

1

Magnitude/ dB

2

Range Slice

0

Magnitude/ dB

2 Azimuth / m

Azimuth / m

3

3

-1 0 1 Range Samples/ m

(b)

-1 0 1 Range Samples/ m

Figure (b)8. Cont. Figure8. 8. Cont. Cont. Figure

-30

-2 0 2 Azimuth Samples/ m

(c)

-2 0 2 Azimuth Samples/ m

(c)

Sensors 2017, 17, 1237

11 of 19

Sensors 2017, 17, 1237 Sensors 2017, 17, 1237

12 of 20 12 of 20

-10 -10 -20 -20 -30 -30

2 2

(d) (d)

0 1 0 /m 1 Range Range / m

-20 -20 -30 -30

-1 0 1 -1Range Samples/ 0 1m Range Samples/ m

-2 0 2 -2 0 2m Azimuth Samples/ Azimuth Samples/ m

(e) Slice Range Range Slice

(f) Azimuth Slice Azimuth Slice

-10 -10 -20 -20 -30 -30

2 2

(f)

0 0 -10 -10 -20 -20 -30 -30

-1 0 1 -1Range Samples/ 0 1 m Range Samples/ m

(g) (g)

Azimuth Slice Azimuth Slice

-10 -10

(e)

0 0

-1 -1

0 0 Magnitude/ dBdB Magnitude/

Magnitude/ dBdB Magnitude/

-1 0 1 -1 Range 0 /m 1 Range / m

Range Slice Range Slice

Magnitude/ dBdB Magnitude/

3 3 2 2 1 1 0 0 -1 -1 -2 -2 -3 -3 -2 -2

0 0

Magnitude/ dBdB Magnitude/

Azimuth / m/ m Azimuth

Azimuth / m/ m Azimuth

3 3 2 2 1 1 0 0 -1 -1 -2 -2 -3 -3 -2 -2

-2 0 2 -2 0 2m Azimuth Samples/ Azimuth (i)Samples/ m

(h) (h)

(i)

-1 -1

0 0 /m Range Range / m

1 1

Magnitude/ dBdB Magnitude/ -1 0 1 -1 Range Samples/ 0 m1 Range Samples/ m

Magnitude/ dB Magnitude/ dB

(d) (d)

-30 -30 -2 0 2 -2 Azimuth Samples/ 0 m 2 Azimuth Samples/ m

(b) (b)

1 1

2 2

Azimuth Slice Azimuth Slice

-20 -20

-30 -30

0 0

0 0 /m Range Range / m

0 0 -10 -10

-20 -20

(a) (a)

-1 -1

Range Slice Range Slice

-10 -10

2 2

Azim uthuth / m/ m Azim

3 3 2 2 1 1 0 0 -1 -1 -2 -2 -3 -3 -2 -2

0 0

(c) (c)

Range Slice Range Slice

-10 -10

0 0 Magnitude/ dB Magnitude/ dB

3 3 2 2 1 1 0 0 -1 -1 -2 -2 -3 -3 -2 -2

Magnitude/ dBdB Magnitude/

Azimuth / m/ m Azimuth

Figure 8. Imaging results of baseband azimuth scaling (BAS). (a–c) show contour, range slice, and Figure 8. 8. Imaging Imagingresults resultsof of baseband basebandazimuth azimuthscaling scaling(BAS). (BAS).(a–c) (a–c) show show contour, contour, range range slice, slice, and and Figure azimuth slice of the target in the far range, respectively. (d–f) show contour, range slice, and azimuth slice of the target in the far range, respectively. (d–f) show contour, range slice, and azimuth azimuth slice of the target in the far range, respectively. (d–f) show contour, range slice, and azimuth slice of the target near the mid-range, respectively. (g–i) show contour, range slice, and slice of the target neartarget the mid-range, (g–i) show contour, rangecontour, slice, andrange azimuth slice of azimuth slice of the near the respectively. mid-range, respectively. (g–i) show slice, and azimuth slice of the target in the near range, respectively. the targetslice in the neartarget range, azimuth of the inrespectively. the near range, respectively.

-20 -20 -30 -30 -1 0 1 m -1Range Samples/ 0 1 Range Samples/ m

(e) (e) Figure 9. Cont. Figure 9. 9. Cont. Cont. Figure

Azimuth Slice Azimuth Slice

-10 -10 -20 -20 -30 -30 -3 -3

-2 -1 0 1 2 Samples/ -2 Azimuth -1 0 1 m 2 Azimuth Samples/ m

(f) (f)

Sensors 2017, 17, 1237 Sensors 2017, 17, 1237 Sensors 2017, 17, 1237

1

2

0

1

0

Range Slice Range Slice 0

-10 -10

0

0 -1 -1 -2 -2 -3 -3 -2 -2

-10 -10

-20 -20

-20 -20

-30 -30

-1 -1

0 0 /m Range Range / m

1

1

2

-30 -30

-1 0 1 -1 Range Samples/ 0 1m Range Samples/ m

2

(g) (g)

Azimuth Slice Azimuth Slice

0

M ag n itu d e/ d B M ag n itu d e/ d B

Azimuth / m Azimuth / m

2

3

Magnitude/ dB Magnitude/ dB

3

12 of 19 13 of 20 13 of 20

-3 -3

(h) (h)

-2 -2

-1 0 1 2 -1 0 1 m 2 Azimuth Samples/ Azimuth Samples/ m

(i) (i)

Figure 9. Imaging results of GCS-BAS. (a–c) show contour, range slice, and azimuth slice of the Figure 9. (a–c) show contour, range slice,slice, and azimuth slice of the of target Figure 9. Imaging Imagingresults resultsofofGCS-BAS. GCS-BAS. (a–c) show contour, range and azimuth slice the target in the far range, respectively. (d–f) show contour, range slice, and azimuth slice of the target in the far range, respectively. (d–f) show contour, range slice, and azimuth slice of the target near the target in the far range, respectively. (d–f) show contour, range slice, and azimuth slice of the target near the mid-range, respectively. (g–i) show contour, range slice, and azimuth slice of the target in mid-range, respectively. (g–i) show contour, range slice, andslice, azimuth slice of the target in target the near near the mid-range, respectively. (g–i) show contour, range and azimuth slice of the in the near range, respectively. range, respectively. the near range, respectively.

1

0 -1 -1 -2 -2 -3 -3 -2 -2

-30 -30

2

2

0

0 -1 -1 -2 -2 -3 -3 -2 -2

1

-1 -1

0 1 0 /m 1 Range Range / m

2

2

(d) (d)

3

0

1

0

0 1 0 /m 1 Range Range / m

(g) (g)

2

2

(f)

Range Slice

0

0

(f) Azimuth Slice Azimuth Slice

-10 -10

-20 -20

-1 -1

-2 0 2 -2Azimuth 0Samples/ 2m Azimuth Samples/ m

-1 0 1 -1 Range Samples/ 0 1 m Range Samples/ m

-10 -10

0

0 -1 -1 -2 -2 -3 -3 -2 -2

-30 -30

(e) (e) Range Slice

2

Azimuth Slice Azimuth Slice

-20 -20

-20 -20 -30 -30

0

-10 -10

Magnitude/ dB Magnitude/ dB

Azim uth / m Azim uth / m

2

0

-10 -10

0

3

0

(c) (c)

Range Slice Range Slice

Magnitude/ dB Magnitude/ dB

Azim uth / m Azim uth / m

1

1

-2 0 2 -2Azimuth 0 Samples/2m Azimuth Samples/ m

(b) (b)

2

2

-30 -30

-1 0 1 -1 Range Samples/ 0 1m Range Samples/ m

(a) (a)

3

3

-20 -20

-20 -20

-1 0 1 -1 0 /m 1 Range Range / m

Azimuth Slice Azimuth Slice

-10 -10

-10 -10

0

0

0

Magnitude/ dB Magnitude/ dB

0

Range Slice Range Slice

Magnitude/ dB Magnitude/ dB

1

2

0

Magnitude/ dB Magnitude/ dB

Az im uth / m Az im uth / m

2

3

Magnitude/ dB Magnitude/ dB

3

-20 -20 -30 -30

-30 -30 -1 0 1 -1 Range Samples/ 0 1m Range Samples/ m

(h) (h)

-2 0 2 -2 Azimuth 0Samples/ 2m Azimuth Samples/ m

(i) (i)

Figure 10. Imaging results of the full-aperture imaging algorithm with ω − k process kernel. (a–c) ω − kk process Figure 10. Imaging results of the full-aperture imaging algorithm with process kernel. kernel. (a–c) Figure 10. Imaging results full-aperture with ω − respectively. show contour, range slice, of andthe azimuth slice ofimaging the targetalgorithm near the far range, (d–f) (a–c) show show contour, range slice, slice, and and azimuth slice of the target near the far range, respectively. (d–f) show show contour, range azimuth slice of the target near the far range, respectively. (d–f) contour, range slice, and azimuth slice of the target in the mid-range, respectively. (g–i)show show contour, rangeslice, slice, and azimuth slice of target the target inmid-range, the mid-range, respectively. (g–i) show contour, ofof the inin the respectively. (g–i) show contour, contour,range range slice,and andazimuth azimuthslice slice the target the near range, respectively. contour, range slice, and azimuth slice of the target in the near range, respectively. range slice, and azimuth slice of the target in the near range, respectively.

Sensors 2017, 17, 1237

13 of 19

Sensors 2017, 17, 1237

14 of 20 Range Slice

3

0 -1 -2

-10

-20

-1

0 Range / m

1

2

-10

-20

-30

-30

-3 -2

-2 0 2 Azimuth Samples/ m

-1 0 1 Range Samples/ m

(a)

(b)

(c) Azimuth Slice

Range Slice

3

0

Magnitude/ dB

1 0 -1 -2

Magnitude/ dB

0

2

-10

-20

-1

0 Range / m

1

2

-20

-1 0 1 Range Samples/ m

(d)

-2 0 2 Azimuth Samples/ m

(e)

(f)

Range Slice

3

Azimuth Slice

0

2 Magnitude/ dB

1 0 -1 -2 -3 -2

-10

-30

-30

-3 -2

0

-10

Magnitude/ dB

Azimuth / m

0

Magnitude/ dB

Magnitude/ dB

Azimuth / m

1

Azimuth / m

Azimuth Slice

0

2

-20

0 Range / m

(g)

1

2

-20

-30

-30

-1

-10

-1 0 1 Range Samples/ m

(h)

-4

-2 0 2 Azimuth Samples/ m

4

(i)

Figure 11. Imaging results of back-projection (BP). (a–c) show contour, range slice, and azimuth slice Figure 11. Imaging results of back-projection (BP). (a–c) show contour, range slice, and azimuth slice of the target near the far range, respectively. (d–f) show contour, range slice, and azimuth slice of of the target near the far range, respectively. (d–f) show contour, range slice, and azimuth slice of the the target in the mid-range, respectively. (g–i) show contour, range slice, and azimuth slice of the target in the mid-range, respectively. (g–i) show contour, range slice, and azimuth slice of the target in target in the near range, respectively. the near range, respectively.

This paper analyses the focusing performance of the three algorithms based on such indices as This analyses focusing performance of the three algorithms based on(Res). such The indices as the peakpaper side-lobe ratiothe (PSLR), integrated side-lobe ratio (ISLR), and resolution indices the peak side-lobe ratio (PSLR), integrated side-lobe ratio (ISLR), and resolution (Res). The indices of of target focusing for the four algorithms are shown in Table 3. target focusing for the four algorithms are shown in Table 3. Figure 8 and Table 3 illustrate that the conventional BAS algorithm causes deterioration of the Figure 8 and Table illustrate that conventional algorithm deterioration the target image when the3system is in thethe large bandwidthBAS mode, such ascauses with the parametersof shown target image when the system is in the large bandwidth mode, such as with the parameters shown in in Table 2. The PSLR, ISLR, and resolution degrade considerably in the two dimensions, Table 2. The PSLR, ISLR, and resolution degrade considerably in the two dimensions, particularly along particularly along the range edges. Degradation is also observed along the azimuth dimension. the range9edges. is also observed along the azimuth dimension. Figure shows thetargets resultsin Figure showsDegradation the results obtained by the GCS-BAS algorithm. The contours for9 different obtained by the GCS-BAS algorithm. The contours for different targets in Figure 9a,d,g, the range Figure 9a,d,g, the range slices in Figure 9b,e,h, and the azimuth slices in Figure 9c,f,i show that the slices in Figure 9b,e,h,isand azimuththe slices inbandwidth Figure 9c,f,isliding show spotlight that the GCS-BAS algorithm is GCS-BAS algorithm ablethe to process large SAR data whether what able process the large sliding SAR data whether the target located the to target is located at bandwidth the edges or in thespotlight center. Figure 10 and Tablewhat 3 illustrate thatis the focus atindexes the edges or in the center. Figure 10 and Table 3 illustrate that the focus indexes of the GCS-BAS of the GCS-BAS algorithm are notably better than the results processed by the conventional algorithm are notably better than the processed byas thethe conventional BAS Table algorithm, which BAS algorithm, which performs inresults a similar manner BP algorithm. 4 shows the performs in a similar manner as the BP algorithm. Table 4 shows the statistical properties of targets statistical properties of targets at the azimuth edge. The standard deviations indicate that at the the azimuth edge. The standard deviations indicate that the full-aperture imaging algorithm with the ω − k full-aperture imaging algorithm with the processing kernel has a focusing capability which ωapproaches − k processing kernel has a focusing capability GCS-BAS. Thethat results Figurestwo 9 GCS-BAS. The results of Figureswhich 9 andapproaches 10 and Table 4 show theofformer imaging algorithms have imaging capabilities which are closer to BP.

Sensors 2017, 17, 1237

14 of 19

and 10 and Table 4 show that the former two imaging algorithms have imaging capabilities which are closer to BP. Table 3. Performance of the three algorithm imaging results. Target Position

Dimension

Index

BAS

GCS-BAS

Full-Aperture Imaging

BP

Azimuth

PSLR/dB ISLR/dB Res/m

−13.6593 −12.4162 0.2490

−13.3273 −10.7198 0.2371

−13.3364 −10.6901 0.2365

−13.4676 −10.6779 0.2359

Range

PSLR/dB ISLR/dB Res/m

−9.7625 −6.2740 0.1379

−12.8918 −10.1871 0.1328

−13.0021 −10.0817 0.1328

−13.2873 −9.9379 0.1328

Azimuth

PSLR/dB ISLR/dB Res/m

−14.0884 −11.7920 0.2321

−13.4163 −10.5783 0.2229

−13.3701 −10.5001 0.2203

−13.6043 −10.1517 0.2191

Range

PSLR/dB ISLR/dB Res/m

−9.5632 −5.6929 0.1359

−13.2104 −9.9811 0.1328

−13.2498 −9.98758 0.1328

−13.2759 −9.9629 0.1328

Azimuth

PSLR/dB ISLR/dB Res/m

−14.8360 −10.9848 0.2371

−13.2982 −10.3992 0.2234

−13.3125 −10.0126 0.2190

−13.6205 −9.9223 0.2148

Range

PSLR/dB ISLR/dB Res/m

−9.4085 −5.6916 0.1395

−12.9537 −9.6545 0.1328

−13.0529 −9.6581 0.1328

−13.2856 −9.9877 0.1328

Far range

Mid-range

Near range

Table 4. The statistical properties of targets at the azimuth edge. Dimension

Range

Azimuth

Index

BAS

GCS-BAS

Full-Aperture Imaging

BP

Res

Mean/m Standard deviation

0.1351 0.0090

0.1327 0.0010

0.1328 0.0012

0.1324 0.0002

PSLR

Mean/dB Standard deviation

−10.8216 7.4778

−13.0939 0.5546

−13.2842 0.1766

−13.2889 0.0558

ISLR

Mean/dB Standard deviation

−7.4356 6.6592

−10.0103 1.2482

−9.9508 1.2809

−9.8435 0.7351

Res

Mean/m Standard deviation

0.2423 0.0238

0.2348 0.0004

0.2354 0.0018

0.2348 0.0004

PSLR

Mean/dB Standard deviation

−13.4884 0.7070

−13.3351 0.1968

−13.3710 0.3371

−13.3012 0.0948

ISLR

Mean/dB Standard deviation

−11.1393 5.0295

−10.6051 3.0250

−10.9491 2.7883

−10.5896 2.6693

However, the four algorithms discussed in this paper has different computational burdens. Suppose that the number of pixels in the range direction and azimuth direction is Nrg and Naz , respectively. Table 5 shows the computational burden of the BAS, GCS-BAS, BP, and full-aperture imaging algorithms with the ω − k processing kernel. Table 5. The comparison of algorithms’ computational burden. Algorithms

Computational Burden

BAS GCS-BAS Full-aperture imaging BP


29Nrg Naz + 20Nrg Naz log2 ( Naz ) + 10Nrg Naz log2 Nrg
64Nrg Naz + 20Nrg Naz log2 ( Naz ) + 20Nrg Naz log2 Nrg
(4Mker + 40) Nrg Naz + 30Nrg Naz log2 ( Naz ) + 10Nrg Naz log2 Nrg 2 Nrg Naz

Sensors 2017, 17, 1237

15 of 19

Mker is the Stolt interpolation length of the ω − k processing kernel. To maintain a high focusing accuracy, Mker is always greater than 16. Figure 12 illustrates the different computational burdens of the four algorithms with different Naz Sensors 2017, 17, 1237 16 of 20 and a constant Nrg . It is easy to conclude that the computational burden of BP is much heavier than that of the other three algorithms. The computational burden of the full-aperture imaging algorithm algorithm with the ω −k kernel processing kernel is slightly heavier than that of BAS-GCS. 5 with the ω − k processing is slightly heavier than that of BAS-GCS. Table 5 illustratesTable that the ω − k illustrates that burden the computational burden of the BPimaging and full-aperture with kernel computational of the BP and full-aperture algorithmimaging with ω algorithm − k processing −k processing kernel are much larger than that of GCS-BAS. The Stolt interpolation in theis ω are much larger than that of GCS-BAS. The Stolt interpolation in the ω − k processing kernel much processing kernel is much more time consuming than the common multiplication operation [24] in more time consuming than the common multiplication operation [24] in BAS or GCS-BAS. Hence, BAS or GCS-BAS. Hence,three compared withdiscussed the other in three discussed has in this paper, the compared with the other algorithms thisalgorithms paper, the GCS-BAS the advantages GCS-BAS has the advantages of computational efficiency and focusing accuracy in the processing of of computational efficiency and focusing accuracy in the processing of the large bandwidth sliding the large bandwidth sliding spotlight SAR. spotlight SAR. 9000

180 BAS GCS-BAS Full-aperture imaging

140

8000 Total computation/[G Times]

Total computation/[G Times]

160

120 100 80 60 40

6000 5000 4000 3000 2000 1000

20 0

BP

7000

0

0.5

1

1.5 2 Azimuth Size

2.5

0

3

0

0.5

1

4

x 10

1.5 2 Azimuth Size

(a)

2.5

3 4

x 10

(b)

Figure12. 12.Computational Computational burden burden of of the the four four algorithms. algorithms. (a) Figure (a) shows shows the the computational computationalburden burdenwith with respect to azimuth size, with Nrg = 8192 , for BAS (solid blue), GCS-BAS (dotted red), and respect to azimuth size, with Nrg = 8192, for BAS (solid blue), GCS-BAS (dotted red), and full-aperture full-aperture imaging black). (b) shows the burden computational burden of to BPazimuth with respect to imaging (dashed black).(dashed (b) shows the computational of BP with respect size, with N = 8192 size, with . Nazimuth = 8192. rg rg

Foraamore moreintuitive intuitive understanding understanding of the computational For computational efficiency efficiencyof ofthe thefour fouralgorithms, algorithms,we we process the the data data on the with an an i7 Dual-Core CPUCPU of process the parameters parametersshown shownininTable Table2 2bybya computer a computer with i7 Dual-Core 3.33.3GHz. raw data dataisis NNrgrg(= . The algorithms’ time 4 0 0 0 0)) × × NNa az = 2 4 0 0 0 24000 of GHz.The Thesize size of of the raw . The algorithms’ time (= ) ( = 40000 ( ) z consumption consumptionisisshown shown in in Table Table6. 6. Table of the the algorithms. algorithms. Table 6. 6. Time Time consumption consumption of

Algorithms Time consumption Algorithms

Time consumption

BAS BAS 1153.1293s 1153.1293 s

GCS-BAS Full-Aperture Imaging Full-Aperture GCS-BAS Imaging 1405.0256s 3685.2732s 1405.0256 s

3685.2732 s

BP BP 160,484.2220s 160,484.2220 s

5. Conclusions

5. Conclusions The large bandwidth sliding spotlight SAR critically couples between the range and azimuth dimensions due to its largesliding bandwidth relative to critically the carriercouples frequency. Its azimuth spectrum is also The large bandwidth spotlight SAR between the range and azimuth aliased for its antenna beam steering in velocity orientation, which causes the total azimuth dimensions due to its large bandwidth relative to the carrier frequency. Its azimuth spectrum is bandwidth to be larger than that of PRF. This paper presents a method called the GCS-BAS also aliased for its antenna beam steering in velocity orientation, which causes the total azimuth algorithm to solve the coupling problem of the large bandwidth sliding spotlight SAR. The bandwidth to be larger than that of PRF. This paper presents a method called the GCS-BAS algorithm algorithm has a flexibility to choose the compensation order for the large bandwidth sliding to solve the coupling problem of the large bandwidth sliding spotlight SAR. The algorithm has a spotlight SAR based on the ratio f 0 / Br . Compared with the up-sampling algorithm and the full-aperture processing algorithm, the GCS-BAS algorithm processes the large bandwidth sliding spotlight SAR data more efficiently without zero-padding or interpolation. The simulation and real data processing indicate that the algorithm is able to solve the phase coupling problem caused by a large bandwidth. Higher precision and accuracy can be achieved for the large bandwidth sliding

Sensors 2017, 17, 1237

16 of 19

flexibility to choose the compensation order for the large bandwidth sliding spotlight SAR based on the ratio f 0 /Br . Compared with the up-sampling algorithm and the full-aperture processing algorithm, the GCS-BAS algorithm processes the large bandwidth sliding spotlight SAR data more efficiently without zero-padding or interpolation. The simulation and real data processing indicate that the algorithm is able to solve the phase coupling problem caused by a large bandwidth. Higher precision and accuracy can be achieved for the large bandwidth sliding spotlight SAR data compared to the conventional BAS algorithm. It is also more efficient than the full-aperture imaging algorithm with the ω − k processing kernel or the BP algorithm at the same focusing depth and accuracy. The simulation of the former four algorithms indicates that a fine focus can be achieved for the large bandwidth sliding spotlight data by the GCS-BAS algorithm, while an appropriate trade-off can be obtained between computational efficiency and imaging accuracy. Acknowledgments: This work was supported by the Chinese National Natural Science Fund under Grant 61401480. The authors would like to thank the anonymous reviewers for their very competent comments and helpful suggestions. Author Contributions: Tianzhu Yi, Zhihua He, and Feng He conceived and designed the experiments; Tianzhu Yi performed the experiments and analyzed the data; Zhen Dong and Manqing Wu contributed materials; and Tianzhu Yi wrote the paper. Conflicts of Interest: The authors declare no conflict of interest.

Appendix A This appendix is intended to help understand the different steps of sub-aperture processing of GCS-BAS. Expanding Equation (3) with the Tyler principle and the sub-aperture data with the phase of Equation (4) in the 2D frequency domain can be described as:

S1 f τ , f η = Wr ( f τ )Wa f η exp( jθ2D )

(A1)


θ2D is described in Equation (4). After the high order phase compensation by H1 f τ , f η , the signal in the 2D domain can be expressed as:



S1 f τ , f η = Wr ( f τ )Wa f η exp jθ1 f τ , f η

(A2)


θ1 f τ , f η is described in Equation (6). A range IFFT results in the signal in the RD domain:



S1 τ, f η = wr (τ )Wa f η exp jθ2 τ, f η

(A3)


θ2 τ, f η is discussed in Equation (7). The SRC factor Km and time in the RD domain are: Kr

Km = 1 − Kr τd =

(A4)

cR0 f η2 2Vr2 f 03 D ( f η )

3

2R0 cD ( f η ) 2Rre f

τre f =

cD ( f η )

(A5)

∆τ = τd − τre f The following factors can be deduced from Equations (A4) and (A5): Ks = Kf =

c2 f η2 4v2 f 03 D2 Kr 1−Kr τre f Ks

The time in the RD domain is represented as the total RCM:

(A6)

Sensors 2017, 17, 1237

17 of 19

=

RCMtotal

2R0 cD ( f η )

2R0  cD f ηre f

−

=

2R0 cD ( f η )

1−

D( fη )   D f ηre f

! (A7)

= τd (1 − α) We define α =

D( fη ) .  D f ηre f

RCMbulk

=

The time delay of the bulk RCM can be expressed as:

2R  re f  cD f η ,Vrre f

−

2Rre f   cD f ηre f ,Vrre f

=

2R  re f  cD f η ,Vrre f

1−

  ! D f η ,Vrre f   D f ηre f ,Vrre f

(A8)

≈ τre f − ατre f The following equation can be obtained from Figure 4: 

∆τ f ηre f



=

2R0  cD f ηre f

−

2Rre f   cD f ηre f

=

D( fη )   D f ηre f



2R0 cD ( f η )

−

2Rre f



cD ( f η )

(A9)

= α∆τ   τs = τre f + ∆τ f ηre f = τre f + α∆τ = τd − ∆τ (1 − α)

(A10)

The differential RCM can be expressed as: RCMdi f f

= RCMtotal  − RCMbulk = τd (1 − α) − τre f (1 − α) = (1 − α) τd − τre f = (1 − α)∆τ

(A11)

Differential RCMC is achieved by the CS operation. From the phase perspective, the CS realizes the addition of Equations (7) and (8). Expand Equation (9) at τ = τs , and the following equation is obtained:

θ2 τ, f η = − 4πRc 0 f0 D f η + πKm ((τ − τs ) − (1 − α)∆τ )2 n 
 i +π ∑ Xi − 2 f 0 D f η γi ∆τ Km ((τ − τs ) − (1 − α)∆τ )i i =3

n

+πq2 ((τ − τs ) + α∆τ )2 + π ∑ qi ((τ − τs ) + α∆τ )i

(A12)

i =3

n
= − 4πRc 0 f0 D f η + πC0 (∆τ ) + ∑ Ci (τ − τs )i i =1

C0 (∆τ ) is the polynomial of ∆τ, and it can be expressed as: C0 (∆τ )

    3 X ∆τ 3 = α2 q2 + (α − 1)2 Km ∆τ 2 + α3 q3 + (α − 1)3 Km 3  n 
i i i −1 i −1 i + ∑ α qi + (α − 1) Km Xi − 2D f η f 0 γi−1 (α − 1) Km ∆τ i

(A13)

i =4

The coefficients of (τ − τs ) are calculated from Equation (A12) (this paper only calculates to C3 due to space restraints).

Sensors 2017, 17, 1237

C1 C2

C3

18 of 19

h i n 3 X ∆τ 2 + = 2[(−1 + α)Km + αq2 ]∆τ + 3α2 q3 + 3(−1 + α)2 Km ∑ C1i ∆τ i 3 i =3   3 X ∆τ = q2 + Km + 3 αq3 + (−1 + α)Km 3 i h n
4 X − (−1 + α ) K 3 D f i 2 +6 α2 q4 + (−1 + α)2 Km f γ η 0 3 ∆τ + ∑ C2i ∆τ 4 m i = 3 
 3 X + 4αq + 4K 4 X (−1 + α ) − 2K 3 D f =hq3 + Km η f 0 γ3i ∆τ 3 4 m 4 m
2 5 X − 8(−1 + α ) K 4 D f + 10α2 q5 + 10(−1 + α)2 Km η f 0 γ4 ∆τ 5 m

(A14)

n

+ ∑ C3i ∆τ i i =3

...... Equation (A4) is expanded at ∆τ = 0 with the Tyler principle. Then Equation (A4) can be expressed as: Kf Km = (A15) ≈ K f + Ks K2f ∆τ + Ks2 K3f ∆τ 2 1 − K f ∆τKs Based on Equations (A15) and (A14), we can obtain the solutions of qi (i ≥ 2) and Xi (i ≥ 3) for setting the first- and second-order coefficients of ∆τ to zero.  K (1− α )  q2 = f α    (1−α)K2f Ks    q3 =  3α     X3 = (−2+α)Ks 3(−1+α)K

f

K3 (2Ks2 +9(−1+α)K f Ks X3 −6D ( f η )(−1+α) f 0 γ3 )   q4 = − f   12α   2Ks2 +9(−2+α)K f Ks X3 −6D ( f η )(−2+α) f 0 γ3   X4 =   12(−1+α)K f   ......

(A16)

After the CS operation, the signal in the RD domain can be expressed as:



S2 τ, f η = wr (τ )Wa f η exp jθ2 τ, f η

(A17)


θ2 τ, f η is described in Equation (10). A range FFT is used for the transformation into the 2D frequency domain. Neglecting the impact of a third order higher-order item, the result of the stationary phase point is τPOSP = Kα f τ + τs . The signal can be expressed as: f




 S3 f τ , f η = Wr ( f τ )Wa f η exp jπθ3 f τ , f η

(A18)


θ3 f τ , f η is described in Equation (11). References 1. 2.

3.

4.

Balz, T.; Caspari, G.; Fu, B.; Liao, M. Discernibility of Burial Mounds in High-Resolution X-Band SAR Images for Archaeological Prospections in the Altai Mountains. Remote Sens. 2016, 8, 817. [CrossRef] Redadaa, S.; Le Caillec, J.M.; Solaiman, B.; Benslama, M. Focusing problems of subsurface imaging by a low-frequency SAR. In Proceedings of the International Geoscience and Remote Sensing Symposium, Barcelona, Spain, 23–27 July 2007; pp. 4101–4104. Ulander, L.M.H.; Blom, M.; Flood, B.; Follo, P.; Frolind, P.O.; Gustavsson, A.; Jonsson, T.; Larsson, B.; Murdin, D.; Pettersson, M.; et al. Development of the largeband LORA SAR operating in the VHF/UHF-band. In Proceedings of the International Geoscience Remote Sensing Symposium, Toulouse, France, 21–25 July 2003; pp. 4268–4270. Ressler, M.A. The Army Research Laboratory ultra wideband Boom-SAR. In Proceedings of the International Geoscience Remote Sensing Symposium, Lincoln, NE, USA, 27–31 May 1996; pp. 1886–1888.

Sensors 2017, 17, 1237

5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17.

18. 19. 20. 21. 22. 23. 24.

19 of 19

Hellsten, H. CARABAS—An UWB low frequency SAR. In Proceedings of the IEEE MTT-S International Microwave Symposium Digest, Albuquerque, NM, USA, 1–5 June 1992; pp. 1495–1498. Rau, R.; McClellan, J.H. Analytic Models and Postprocessing Techniques for UWB SAR. IEEE Trans. Aerosp. Electron. Syst. 2000, 36, 1058–1074. Prats, P.; Scheiber, R.; Mittermayer, J.; Meta, A.; Moreria, A. Processing of Sliding Spotlight and TOPS SAR Data Using Baseband Azimuth Scaling. IEEE Trans. Geosci. Remote Sens. 2010, 48, 770–780. [CrossRef] Davidson, G.W.; Cumming, I.G. A Chirp Scaling Approach for Processing Squint Mode SAR Data. IEEE Trans. Aerosp. Electron. Syst. 1996, 32, 121–133. [CrossRef] Zaugg, E.C.; Long, D.G. Generalized Frequency-Domain SAR Processing. IEEE Trans. Geosci. Remote Sens. 2009, 47, 3761–3773. [CrossRef] Lanari, R.; Tesauro, M.; Sansosti, E.; Fornaro, G. Spotlight SAR data focusing based on a two-step processing approach. IEEE Trans. Geosci. Remote Sens. 2001, 39, 1993–2004. [CrossRef] Tang, Y.; Wang, Y.F.; Zhang, B.C. A Study of Sling Spotlight SAR Imaging Mode. J. Electron. Inf. Technol. 2007, 29, 26–29. Wei, X.; Deng, Y.K. Full-Aperture SAR Data Focusing in the Spaceborne Squinted Sliding spotlight Mode. IEEE Trans. Geosci. Remote Sens. 2014, 52, 4596–4607. [CrossRef] Mittermayer, J.; Moreira, A.; Loffeld, O. Spotlight SAR Data Processing Using the Frequency Scaling Algorithm. IEEE Trans. Geosci. Remote Sens. 1999, 37, 2198–2214. [CrossRef] Sun, G.C.; Xing, M.D.; Liu, Y.; Sun, L.; Bao, Z.; Wu, Y. Extended NCS Based on Method of Series Reversion for Imaging of Highly Squinted SAR. IEEE Geosci. Remote Sens. Lett. 2011, 8, 446–450. [CrossRef] Neo, Y.L.; Wong, F.; Cumming, I.G. A Two-Dimensional Spectrum for Bistatic SAR Processing Using Series Reversion. IEEE Geosci. Remote Sens. Lett. 2007, 4, 93–96. [CrossRef] Luo, X.L.; Deng, Y.K.; Wang, R.; Xu, W.; Luo, Y.H.; Guo, L. Image Formation Processing for Sliding Spotlight SAR With Stepped Frequency Chirps. IEEE Geosci. Remote Sens. Lett. 2014, 11, 1692–1696. Zheng, X.; Yu, W.; Li, Z. A Novel Algorithm for Wide Beam SAR Motion Compensation Based on Frequency Division. In Proceedings of the IEEE International Symposium on Geoscience and Remote Sensing, Beijing, China, 10–15 July 2016; pp. 3160–3163. Wang, K.; Liu, X. Quartic-phase algorithm for highly squinted SAR data processing. IEEE Geosci. Remote Sens. Lett. 2007, 4, 246–250. [CrossRef] Wong, F.H.; Cumming, I.G.; Neo, Y.L. Focusing Bistatic SAR Data Using the Nonlinear Chirp Scaling Algorithm. IEEE Trans. Geosci. Remote Sens. 2008, 46, 2493–2505. [CrossRef] Zaugg, E.C.; Long, D.G. Generalized Frequency Scaling and Backprojection for LFM-CW SAR Processing. IEEE Trans. Geosci. Remote Sens. 2015, 53, 3600–3614. [CrossRef] Raney, R.K.; Runge, H.; Bamler, R.; Cumming, I.G.; Wong, F.H. Precision SAR Processing Using Chirp Scaling. IEEE Trans. Geosci. Remote Sens. 1994, 32, 786–799. [CrossRef] He, F.; Chen, Q.; Dong, Z.; Sun, Z. Processing of Ultrahigh-Resolution Spaceborne Sliding Spotlight SAR Data. IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 1–21. [CrossRef] Lanari, R.; Tesauro, M.; Sansosti, E.; Fornaro, G. Spotlight SAR Data Focusing Based on a Two-Step Processing Approach. IEEE Trans. Geosci. Remote Sens. 2001, 39, 1993–2004. [CrossRef] Cumming, I.G.; Wong, F.H. Digital Processing of Synthetic Aperture Radar Data Algorithms and Implementation; ARTECH HOUSE, Inc.: Norwood, MA, USA, 2005. © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).