A Wide-Swath Spaceborne TOPS SAR Image Formation Algorithm Based on Chirp Scaling and Chirp-Z Transform

Based on the terrain observation by progressive scans (TOPS) mode, an efficient full-aperture image formation algorithm for focusing wide-swath spaceborne TOPS data is proposed. First, to overcome the Doppler frequency spectrum aliasing caused by azimuth antenna steering, the range-independent derotation operation is adopted, and the signal properties after derotation are derived in detail. Then, the azimuth deramp operation is performed to resolve image folding in azimuth. The traditional dermap function will introduce a time shift, resulting in appearance of ghost targets and azimuth resolution reduction at the scene edge, especially in the wide-swath coverage case. To avoid this, a novel solution is provided using a modified range-dependent deramp function combined with the chirp-z transform. Moreover, range scaling and azimuth scaling are performed to provide the same azimuth and range sampling interval for all sub-swaths, instead of the interpolation operation for the sub-swath image mosaic. Simulation results are provided to validate the proposed algorithm.


Introduction
State-of-the-art spaceborne SAR systems are capable of operating in several novel imaging modes by adopting the electronic beam steering technique in both elevation and along-track directions, and one such an example is the terrain observation with progressive scan (TOPS) mode [1][2][3]. The TOPS mode was presented by Zan and Guarnieri in [1] for wide-swath observation, similar to the scanSAR mode [3]. Compared with the scanSAR mode, the TOPS mode can avoid scalloping and the azimuth-dependent distributed target ambiguity ratio effect on SAR images, improving the image radiometric quality significantly [3]. Therefore, the TOPS mode is not only widely used in current spaceborne SAR systems, such as Sentinel-1 and TerraSAR-X/TanDem-X [4][5][6], but also adopted in future SAR systems, such as TerraSAR-X2 [7].
The TOPS mode scans the scene with very long bursts and rotates the antenna beam throughout the acquisition from backward to forward in azimuth, resulting in a wider total Doppler bandwidth compared with the pulse repetition frequency (PRF) and a longer strip compared with the standard strip mode within the same time span [1], so classic strip algorithms are not suitable for processing the TOPS data due to Doppler frequency spectrum aliasing and image folding in azimuth.
Several algorithms were developed to overcome these problems. The wavenumber domain algorithm (WDA) kernel was first adopted for TOPS data processing by adding pre-and post-  (1) where D is the antenna length, and    r is: According to Equations (1) and (2), the azimuth resolution is dependent on r , which means that targets located at the near slant range can have a higher azimuth resolution than those located at the far slant range.
As shown in Figure 2, assuming a linear frequency modulated (LFM) pulse is transmitted by the radar, after demodulation to the baseband, the received signal for point target A can be described as [13]: where  [ ] rect represents the rectangular window, b is the signal frequency modulation rate, c is the speed of light,  is the carrier frequency,  and  are range time and azimuth time, respectively, The acquisition geometry of the TOPS mode in the slant range plane is shown in Figure 2, where r s and r represent the minimum distance from the virtual rotation point to the SAR sensor and to the ground, respectively, v e is the effective velocity, T B is the burst time, A is the point target located at (x a , r), X s corresponds to the valid area with full illumination time, and X f corresponds to the area with insufficient illumination time. The TOPS mode employs a counter-clockwise antenna beam steering at a constant rotation rate of ω ϕ , leading to reduction of target illumination time. By adopting the hybrid factor γ(r) [1], the azimuth resolution ρ a can be approximately given by: where D is the antenna length, and γ(r) is: According to Equations (1) and (2), the azimuth resolution is dependent on r, which means that targets located at the near slant range can have a higher azimuth resolution than those located at the far slant range.
As shown in Figure 2, assuming a linear frequency modulated (LFM) pulse is transmitted by the radar, after demodulation to the baseband, the received signal for point target A can be described as [13]: where rect[·] represents the rectangular window, b is the signal frequency modulation rate, c is the speed of light, λ is the carrier frequency, τ and η are range time and azimuth time, respectively, X f = λr/D is the length of azimuth antenna beam footprint, and R(η; x a , r) is the distance between the satellite and the point target A(x a , r). Defining x a = v e η a , R(η; x a , r) can be expressed as follows:

Signal Model
TOPS is proposed as a wide-swath imaging mode with low/medium resolutions. So, the range cell migration (RCM) is relatively small, which can be corrected by classic algorithms, such as the chirp scaling algorithm (CSA) [22], and the Range-Doppler algorithm (RDA) [23]. After range cell migration correction (RCMC), the azimuth signal expression is given by: Applying the azimuth fast Fourier transform (FFT) to Equation (5), the azimuth signal Doppler spectrum is: The properties of azimuth signal Doppler spectrum have been analyzed in detail in [1]. As the time-frequency diagram (TFD) shown in Figure 3 where rotation k is the sweep rate [6]:

Signal Model
TOPS is proposed as a wide-swath imaging mode with low/medium resolutions. So, the range cell migration (RCM) is relatively small, which can be corrected by classic algorithms, such as the chirp scaling algorithm (CSA) [22], and the Range-Doppler algorithm (RDA) [23]. After range cell migration correction (RCMC), the azimuth signal expression is given by: Applying the azimuth fast Fourier transform (FFT) to Equation (5), the azimuth signal Doppler spectrum is: where k s (r) is the Doppler rate, and B ∆θ is the Doppler bandwidth corresponding to the instantaneous antenna beamwidth, as given below: The properties of azimuth signal Doppler spectrum have been analyzed in detail in [1]. As the time-frequency diagram (TFD) shown in Figure 3, the total azimuth bandwidth B T , consisting of B ∆θ and B steer , may span over serval PRF intervals that indicated by f pr f , which results in Doppler frequency spectrum aliasing. B steer denotes the Doppler bandwidth resulting from azimuth antenna-beam steering, given by: where k rotation is the sweep rate [6]: In order to overcome spectrum aliasing, some methods have been proposed, including frequency mosaic [1], sub-aperture processing [10], and the derotation operation [13,15]. Among them the derotation operation is the most efficient. Moreover, the related range-dependent and range-independent functions were derived in [10,13], respectively.
Another problem in TOPS imaging is image folding in azimuth and one solution is the deramp technique. However, the time shift caused by the deramp operation requires an increased PRF [14]. On the other hand, the wide-swath coverage needs a low PRF for a large echo-receiving window, which leads to contradiction between efficient processing and wide-swath coverage. As a solution, in this work, a novel image formation algorithm is proposed for wide-swath coverage TOPS data processing, which can overcome the problems of spectrum aliasing, image folding and time shift without sub-aperture combination.

Azimuth De-Rotation
Similar to the case of sliding spotlight [24,25], the range-independent de-rotation function is employed to overcome the Doppler frequency spectrum aliasing problem, with the expression given by: The de-rotation operation involves azimuth signal convolution between the azimuth signal and de-rotation function. Based on the principle of stationary phase (POSP) [26], the convolution result is given by:  In order to overcome spectrum aliasing, some methods have been proposed, including frequency mosaic [1], sub-aperture processing [10], and the derotation operation [13,15]. Among them the derotation operation is the most efficient. Moreover, the related range-dependent and range-independent functions were derived in [10,13], respectively.
Another problem in TOPS imaging is image folding in azimuth and one solution is the deramp technique. However, the time shift caused by the deramp operation requires an increased PRF [14]. On the other hand, the wide-swath coverage needs a low PRF for a large echo-receiving window, which leads to contradiction between efficient processing and wide-swath coverage. As a solution, in this work, a novel image formation algorithm is proposed for wide-swath coverage TOPS data processing, which can overcome the problems of spectrum aliasing, image folding and time shift without sub-aperture combination.

Azimuth De-Rotation
Similar to the case of sliding spotlight [24,25], the range-independent de-rotation function is employed to overcome the Doppler frequency spectrum aliasing problem, with the expression given by: The de-rotation operation involves azimuth signal convolution between the azimuth signal and de-rotation function. Based on the principle of stationary phase (POSP) [26], the convolution result is given by: where: In the next, we analyze the value range of the azimuth signal in detail. According to Equation (12), after de-rotation, the value range of the azimuth signal is determined by three terms. The first term limits the value of η by: With respect to η a , the value range is determined by the second term: Considering the third term, the value of η varies with η a . As a result: Substituting Equation (15) into Equation (16), we can obtain the value of η corresponding to targets located at the azimuth edge, with η a = X s /(2v e ) and η a = −X s /(2v e ), respectively: Comparing Equations (14) and (17), the boundary value difference can be calculated as follows: In the same way, we can calculate the boundary value difference with respect to the target located at η a = −X s /(2v e ), which is also equal to zero. As shown in Figure 4, the range of η 2 is always wider than η 1 for all targets located in X s . So, the first term plays a dominant role in Equation (12), and the azimuth signal after de-rotation can be rewritten as follows: After de-rotation, the equivalent PRF, referred to as f pr f , is given by: where N 1 is the azimuth output point number after de-rotation. Moreover, f pr f should be larger than B T to resolve spectrum aliasing, which determines the value of N 1 by [12]: where N A is the azimuth point number of raw data. With azimuth FFT, S A,η (η; η a , r) is transformed into the azimuth Doppler domain, yielding:  Figure 4. The value range of azimuth signal after a derotation operation, which is limited by Equation (14) (black solid rectangular window), Equation (17) Since the range-independent derotation operation is performed, the azimuth signal is completely overlapped in the time domain, as shown in Equation (20). Moreover, the time domain width of the output signal, referred to as DE T , is wider than the time domain width of the azimuth signal, as shown in Equation (24), which means that the derotation operation will not cause time aliasing: Furthermore, according to Equation (23) Figure 5 shows the azimuth signal time-frequency diagram after the derotation operation.

Range Scaling
The TOPS sub-swaths data have different signal sampling rate s f , resulting in different sampling intervals in the range dimension. Consequently, resampling in range is necessary for subswaths combination, requiring additional computation. In order to avoid the resampling operation, Figure 4. The value range of azimuth signal after a derotation operation, which is limited by Equation (14) (black solid rectangular window), Equation (17) (red dashed rectangular window), and Equation (18) (blue dash-dot rectangular window).
Since the range-independent derotation operation is performed, the azimuth signal is completely overlapped in the time domain, as shown in Equation (20). Moreover, the time domain width of the output signal, referred to as T DE , is wider than the time domain width of the azimuth signal, as shown in Equation (24), which means that the derotation operation will not cause time aliasing: Furthermore, according to Equation (23), the frequency domain width after derotation is limited by the first term, so targets located at different azimuth positions are separated, with the accumulated Doppler bandwidth B tar (r) and Doppler centroid f η,d (r) given by: Note that B tar (r) and f η,d (r) are range-dependent, but the total azimuth bandwidth B T is range-independent. Figure 5 shows the azimuth signal time-frequency diagram after the derotation operation.  Figure 4. The value range of azimuth signal after a derotation operation, which is limited by Equation (14) (black solid rectangular window), Equation (17) Since the range-independent derotation operation is performed, the azimuth signal is completely overlapped in the time domain, as shown in Equation (20). Moreover, the time domain width of the output signal, referred to as DE T , is wider than the time domain width of the azimuth signal, as shown in Equation (24), which means that the derotation operation will not cause time aliasing: Furthermore, according to Equation (23)

Range Scaling
The TOPS sub-swaths data have different signal sampling rate s f , resulting in different sampling intervals in the range dimension. Consequently, resampling in range is necessary for subswaths combination, requiring additional computation. In order to avoid the resampling operation,

Range Scaling
The TOPS sub-swaths data have different signal sampling rate f s , resulting in different sampling intervals in the range dimension. Consequently, resampling in range is necessary for sub-swaths combination, requiring additional computation. In order to avoid the resampling operation, range scaling is performed in the range-Doppler domain, providing the same sampling intervals. According to the range scaling principle proposed in [14], the range scaling function is given by: where r re f is the reference range, b f a ; r re f and 2R f η ; r re f /c are the effective FM chirp rate in range and the time delay in range-Doppler domain, respectively, whose expression can be found in [22], and Cs (i) ( f a ) is the range scaling factor: where i stands for the i-th sub-swath, and α (i) is given by: with ∆τ (i) being the range signal sampling rate of the i-th sub-swath, and ∆τ 0 the chosen output signal sampling rate after range scaling.

Azimuth Deramp
After range scaling, the standard CSA kernel [22] is adopted for range cell migration correction, range compression, second range compression (SRC), and hyperbolic phase removal. Moreover, the residual phase caused by de-rotation is compensated by: However, if the azimuth inverse FFT (IFFT) is applied directly to focus the image after residual phase compensation, image folding in azimuth will likely occur. To solve the problem, a constant linear frequency modulation for the whole sub-swath is introduced, given by: where: k e r re f = 2v 2 e λ r s + r re f Then, after IFFT, the azimuth signal in the time domain is given by: According to Equation (33), the deramp operation leads to a significant range-dependent time shift given by: Choosing r re f as the middle range, the maximum value of ∆T(r, η a ) is approximately given by: where β is the incident angle, X s and X r are the swath-width in azimuth and range, respectively. According to Equation (35), since the time shift is proportional to the sub-swath coverage, time aliasing will occur for wide-swath coverage TOPS data processing, especially at the scene edge. With the parameters listed in Table 1, Figure 6 shows the time aliasing area, which is highlighted in red. where  is the incident angle, s X and r X are the swath-width in azimuth and range, respectively.
According to Equation (35), since the time shift is proportional to the sub-swath coverage, time aliasing will occur for wide-swath coverage TOPS data processing, especially at the scene edge. With the parameters listed in Table 1, Figure 6 shows the time aliasing area, which is highlighted in red.   The time aliasing effects of point targets P 1 , P 2 and P 3 are shown in Figure 7, with their positions indicated in Figure 6. Since P 1 and P 3 are located at the time aliasing area, significant time aliasing has occurred as shown in Figure 7a, which results in ghost appearance and contour plot distortion, with reduced azimuth resolution, as shown in Figure 7b. where  is the incident angle, s X and r X are the swath-width in azimuth and range, respectively.
According to Equation (35), since the time shift is proportional to the sub-swath coverage, time aliasing will occur for wide-swath coverage TOPS data processing, especially at the scene edge. With the parameters listed in Table 1, Figure 6 shows the time aliasing area, which is highlighted in red.

Modified Azimuth Deramp and Chirp-z
In order to overcome the time shift problem, data division before the deramp operation can be applied, and then parameters are updated for different subdata blocks [27]. However, due to the data division, additional subdata combinations are acquired, which will lead to reduction in processing efficiency. Here, a modified method is proposed, using the range-dependent deramp function given by   e k r is range-dependent, azimuth geometry distortion will take place, which will be analyzed in next section. As for geometry distortion, interpolations are usually performed to realize the same Doppler frequency sampling interval, resulting in additional computations. In this paper, the chirp-z transform is adopted to accommodate geometry distortion, by which more accurate geometry correction results can be obtained without interpolation. The chirp-z transform in the discrete domain is given by [28,29]:

Modified Azimuth Deramp and Chirp-z
In order to overcome the time shift problem, data division before the deramp operation can be applied, and then parameters are updated for different subdata blocks [27]. However, due to the data division, additional subdata combinations are acquired, which will lead to reduction in processing efficiency. Here, a modified method is proposed, using the range-dependent deramp function given by H deramp f η , r and: With the range-dependent deramp function, we can get the time domain expression of the azimuth signal by azimuth IFFT: Compared with Equation (33), the time shift is completely compensated w.r.t. the first term. However, a range-dependent quadric phase is introduced, which should be compensated by: By azimuth IFFT, targets are focused in the frequency domain, with position at k e (r)η a . Since k e (r) is range-dependent, azimuth geometry distortion will take place, which will be analyzed in next section.
As for geometry distortion, interpolations are usually performed to realize the same Doppler frequency sampling interval, resulting in additional computations. In this paper, the chirp-z transform is adopted to accommodate geometry distortion, by which more accurate geometry correction results can be obtained without interpolation. The chirp-z transform in the discrete domain is given by [28,29]: where I is the final focused result, CZT represents the chirp-z transform, ⊗ * is the convolution operation which can be implemented by FFT and complex multiplications as shown in Figure 8.
For Q and W, they are given by: where ∆ f η (r) is the output frequency sampling interval, given by: So, the equivalent time sampling interval is: γ r re f f pr f

(43)
Note that ∆η is range-independent, which means the azimuth geometry distortion is corrected. Moreover, we can also select a fixed time sampling interval for all the sub-swathes, which avoids azimuth interpolation for the sub-swath image mosaic.  where I is the final focused result, CZT represents the chirp-z transform,   is the convolution operation which can be implemented by FFT and complex multiplications as shown in Figure 8. For Q and W , they are given by: where    f r is the output frequency sampling interval, given by: So, the equivalent time sampling interval is: Note that   is range-independent, which means the azimuth geometry distortion is corrected. Moreover, we can also select a fixed time sampling interval for all the sub-swathes, which avoids azimuth interpolation for the sub-swath image mosaic. The flowchart of the proposed algorithm is given in Figure 9. First, in the derotation step, the FFT and complex multiplications are applied instead of the convolution operation in real implementation, including the azimuth dechirp, azimuth FFT and azimuth rechirp [12]. Then, after transforming the data to the range-Doppler domain using azimuth FFT, range scaling can be implemented by multiplying the data by scl H r , which should be updated in range. Finally, the chirp-z transform is adopted to accommodate geometry distortion, which has been discussed earlier in this section. The flowchart of the proposed algorithm is given in Figure 9. First, in the derotation step, the FFT and complex multiplications are applied instead of the convolution operation in real implementation, including the azimuth dechirp, azimuth FFT and azimuth rechirp [12]. Then, after transforming the data to the range-Doppler domain using azimuth FFT, range scaling can be implemented by multiplying the data by H scl,r f η , τ . Next, transforming the data to the 2-D frequency domain by range FFT, the standard CSA kernel is employed for range cell migration correction, range compression, second range compression (SRC). After range IFFT, range compression is completed, and the data is transformed back to the range-Doppler domain, in which the hyperbolic phase removal is performed according to the standard CSA kernel. Note that the residual phase compensation and the range-dependent deramp operations are also implemented in this step by multiplying the data by H de_com f η and H deramp f η , r , respectively. After azimuth IFFT, the data is transformed to the time domain, and we can compensate the quadric phase by multiplying the data by H scl (η; r), which should be updated in range. Finally, the chirp-z transform is adopted to accommodate geometry distortion, which has been discussed earlier in this section.

Simulation Results and Discussion
In order to validate the proposed algorithm, the TOPS raw data of nine point targets were simulated using the parameters listed in Table 1. The nine point targets are arranged in a scene of 50 km × 50 km, as shown in Figure 6.

Comparative Experiments and Discussion
To show the performance of the proposed algorithm, time aliasing, geometry distortion and imaging results are disscussed with the traditional range-independent deramp algorithm and our proposed algorithm.
As shown in Figure 7a, using the range-independent deramp, significant time aliasing occurs, which leads to ghost target appearance and azimuth resolution reduction. By implementing the range-dependent deramp operation proposed in our algorithm, time aliasing is completely removed as shown in Figure 10.

Comparative Experiments and Discussion
To show the performance of the proposed algorithm, time aliasing, geometry distortion and imaging results are disscussed with the traditional range-independent deramp algorithm and our proposed algorithm.
As shown in Figure 7a, using the range-independent deramp, significant time aliasing occurs, which leads to ghost target appearance and azimuth resolution reduction. By implementing the range-dependent deramp operation proposed in our algorithm, time aliasing is completely removed as shown in Figure 10. Moreover, Figure 11 illustrates the geometry correction results with the traditional FFT transform and chrip-z transform proposed in our algorithm. Due to the range-dependent deramp, the azimuth sampling rate varies with range bins. By using the FFT transform, it causes clear geometry distortions, especially at the azimuth edge, resulting in an offset of 69 pixels, as shown in Figure 11a. By employing the chirp-z transform, geometry distortions are totally compensated, as shown in Figure 11b.  Moreover, Figure 11 illustrates the geometry correction results with the traditional FFT transform and chrip-z transform proposed in our algorithm. Due to the range-dependent deramp, the azimuth sampling rate varies with range bins. By using the FFT transform, it causes clear geometry distortions, especially at the azimuth edge, resulting in an offset of 69 pixels, as shown in Figure 11a. By employing the chirp-z transform, geometry distortions are totally compensated, as shown in Figure 11b.

Comparative Experiments and Discussion
To show the performance of the proposed algorithm, time aliasing, geometry distortion and imaging results are disscussed with the traditional range-independent deramp algorithm and our proposed algorithm.
As shown in Figure 7a, using the range-independent deramp, significant time aliasing occurs, which leads to ghost target appearance and azimuth resolution reduction. By implementing the range-dependent deramp operation proposed in our algorithm, time aliasing is completely removed as shown in Figure 10. Moreover, Figure 11 illustrates the geometry correction results with the traditional FFT transform and chrip-z transform proposed in our algorithm. Due to the range-dependent deramp, the azimuth sampling rate varies with range bins. By using the FFT transform, it causes clear geometry distortions, especially at the azimuth edge, resulting in an offset of 69 pixels, as shown in Figure 11a. By employing the chirp-z transform, geometry distortions are totally compensated, as shown in Figure 11b.  Finally, Table 2 shows the image results of P 1 , P 2 and P 3 , including the resolution, the peak sidelobe ratio (PSLR) and the integrated sidelobe ratio (ISLR). The corresponding interpolated contour plots are given in Figures 7b and 12a-c, respectively. Comparing the image results, the traditional algorithm and the proposed algorithm have nearly the same resolution, PLSR and ISLR for point target P 2 . However, as for P 1 and P 3 , the traditional algorithm suffers from significant resolution and PSLR degradation at scene edge, while the proposed algorithm has excellent focusing performance, which means the proposed algorithm is more suitable for wide-swath TOPS data processing.

Imaging Simulation Experiments and Analysis
To show the focusing performance in the whole scene shown in Figure 6, the interpolated contour plot of the nine point targets and the corresponding quality measurement results are shown in Figure 12 and Table 3, respectively. All nine point targets are well focused, with no more than 1% resolution deviations compared with the theoretical values, which verifies the accuracy of the proposed algorithm. Moreover, the simulation data are processed without data division and interpolation operations by employing the proposed algorithm, which means the proposed algorithm is more efficient. Furthermore, to meet the requirement of more applications, such as high-resolution image formation, the proposed algorithm for sliding spotlight data processing will be investigated as a topic of our future study. Finally, Table 2 shows the image results of 1 P , 2 P and 3 P , including the resolution, the peak sidelobe ratio (PSLR) and the integrated sidelobe ratio (ISLR). The corresponding interpolated contour plots are given in Figures 7b and 12a-c, respectively. Comparing the image results, the traditional algorithm and the proposed algorithm have nearly the same resolution, PLSR and ISLR for point target 2 P . However, as for 1 P and 3 P , the traditional algorithm suffers from significant resolution and PSLR degradation at scene edge, while the proposed algorithm has excellent focusing performance, which means the proposed algorithm is more suitable for wide-swath TOPS data processing.

Imaging Simulation Experiments and Analysis
To show the focusing performance in the whole scene shown in Figure 6, the interpolated contour plot of the nine point targets and the corresponding quality measurement results are shown in Figure 12 and Table 3, respectively. All nine point targets are well focused, with no more than 1% resolution deviations compared with the theoretical values, which verifies the accuracy of the proposed algorithm. Moreover, the simulation data are processed without data division and interpolation operations by employing the proposed algorithm, which means the proposed algorithm is more efficient. Furthermore, to meet the requirement of more applications, such as high-resolution image formation, the proposed algorithm for sliding spotlight data processing will be investigated as a topic of our future study.

Conclusions
In this paper, a novel imaging algorithm has been proposed for processing wide-swath coverage spaceborne terrain observation by progressive scans (TOPS) SAR data, without raw data division and any interpolations. To overcome the Doppler frequency aliasing problem, a range-independent rotation operation was applied to compensate the Doppler bandwidth resulting from azimuth antenna-beam steering, and the azimuth signal properties after derotation were derived and discussed in detail. Then, range scaling was performed to provide the same range sampling interval for all sub-swath data, avoiding the additional range resampling operation. Moreover, a modified azimuth deramp method was presented to remove the time shift by adopting a range-dependent dermap function, which solves the time aliasing problem and avoids the appearance of ghost targets. Furthermore, to correct geometry distortions caused by the range-dependent dermap, the chirp-z transform is employed, and the same azimuth sampling interval for all sub-swath data can be obtained. As demonstrated by simulations results, the proposed imaging algorithm has performed well for the wide-swath scene. Finally, the propose algorithms also can be adopted for sliding spotlight data processing, which will be studied in future research work.