A Modiﬁed Cartesian Factorized Backprojection Algorithm Integrating with Non-Start-Stop Model for High Resolution SAR Imaging

: High resolution synthetic aperture radar (SAR) imaging has extensive application value especially in military reconnaissance and disaster monitoring. The motion of the satellite during the transmission and reception of the signal introduces notable errors in the high resolution SAR spotlight mode, which will lead to a defocused SAR image if not handled. To address this problem, an accurate correct echo model based on non-start-stop model is derived to describe the property of the SAR signal in the paper. Then, in the imaging processing, an azimuth-time-varying range frequency modulation rate is used for range compression. The range history and compensation phase are also derived based on the correct echo model. Then, combining the correct echo model and Cartesian factorized backprojection (CFBP) algorithm, a modiﬁed CFBP algorithm is proposed for SAR imaging to improve the accuracy and efﬁciency of processing. Besides, the inﬂuence of residual error due to mismatch is analyzed in detail. In the end, the simulation experiment and Gaofen-3 (GF-3) data experiment are carried out to demonstrate the feasibility of the proposed algorithm.


Introduction
Synthetic aperture radar (SAR), which can provide remote sensing images during the day and night regardless of weather conditions, plays an important role in Earth and planetary observation [1][2][3]. In recent years, high resolution SAR imaging attracts growing interest due to its wide applications especially in military reconnaissance and disaster monitoring [4][5][6][7][8][9][10]. Several spaceborne SAR systems have the ability of submeter imaging. For example, the COSMO-SkyMed SAR can achieve a high azimuth resolution better than 0.9 m [6], while TerraSAR-X can achieve 0.16 m azimuth resolution [7]. In addition, the Gaofen-3 (GF-3) satellite, a C-band multi-polarization SAR of China launched in 2016, can achieve 0.4 m azimuth resolution [5,10]. The realization of high resolution SAR imaging not only puts forward new requirements to the hardware of SAR system but also poses new challenges to the imaging algorithm.
In the high resolution spotlight mode, the satellite steers its antenna beam to continuously illuminate the same area of interest. As a result, the raw signal azimuth bandwidth is larger than pulse repetition frequency (PRF) [11], which increases the difficulty of imaging processing compared with the stripmap mode. A standard approach employed when processing the SAR data is based on the start-stop approximation, with the assumption that the radar transmits and receives the signal at the same position [12]. Although it has been widely used in SAR echo model, it is no longer applicable in high resolution spotlight case [7]. Therefore, the non-start-stop model should be used in imaging processing [13]. Additionally, due to the long trajectory, the error introduced by curved orbit should also be corrected [6]. Topography dependence and atmospheric effects should be considered [7,14], which further increases the complexity of imaging processing.
In order to solve the problems mentioned above, numerous algorithms have been proposed for high resolution SAR imaging. The full-aperture algorithm and subaperture algorithm are two typical frequency domain algorithms. The two-step processing approach [11], integrated with the efficiency of SPECAN algorithms and the precision of stripmap focusing techniques, is a typical full aperture approach and has been further studied in [15]. The subaperture algorithm first divides the raw data into subapertures in the azimuth time domain and the high-resolution is obtained by synthesizing the coarse resolution images of each subaperture [7,16]. The subaperture algorithm solves the curved orbit problem by motion compensation. However, there are residual errors due to the variance of echo signal property along range and azimuth directions. Besides, there are some inconveniences for frequency domain algorithms to adapt to the topography variance in the ground. The time domain algorithm is available for different SAR configurations and the imaging result on the ground can be obtained directly by time domain algorithm [14,17,18]. The backprojection algorithm (BPA), a typical time domain algorithm, has been widely used in SAR imaging [19,20]. The disadvantage of BPA is its high computational burden. There are two strategies for improving the efficiency of BPA. One strategy is using parallel computing tool, for example, GPU. In recent years, the development of GPU has made significant progress. The parallel BPA based on the GPU has been attracting more and more attention and there is great progress in the acceleration of the BPA [19][20][21]. Another strategy is developing the fast algorithms. Numerous fast algorithms of BPA have been developed in the last few decades [22][23][24][25][26][27][28][29]. The fast factorized BP (FFBP) algorithm, which splits long synthetic aperture into short subapertures to generate coarse resolution images and then fuses coarse resolution images to form the final image via interpolation, is an effective algorithm to reduce the computational burden [23,30]. However, the massive interpolations limit the efficiency and introduce errors that result in artifacts in the final image. The accelerating fast BP (AFBP) algorithm is proposed for the high-resolution SAR imaging [29]. The AFBP algorithm only contains two stages: subimage formation in unified polar coordinate and subimage fusion in the 2-D wavenumber domain, which can improve efficiency compared with FFBP algorithm [29]. The Cartesian factorized BP (CFBP) algorithm is proposed in [25,31] and has been further studied in [26,32,33]. Two spectrum compressing filters are applied in CFBP algorithm, which can decrease the cross-range Nyquist sampling rate (NSR) enormously, therefore reducing the computational burden. In addition, during image combination, interpolation can be easily achieved by zero-padding in the azimuth frequency domain, leading to a substantial improvement on both the accuracy and the efficiency [33].
In this paper, a modified CFBP algorithm is proposed for high resolution spotlight SAR imaging. Since the traditional echo model is derived based on start-stop approximation, there will be intolerable error in high resolution case, which will cause defocusing of the target. In order to solve this problem, an accurate correct echo model is firstly derived. The non-start-stop model and curved orbit are all taken into consideration in the correct model. Then, based on the correct echo mode, an azimuth-time-varying range frequency modulation (FM) rate is used for range compression. The position of the target after range compression is also derived. In addition, the azimuth phase in the correct echo model is also changed compared with the traditional echo mode. A discussion for residual phase error after range compression is given, which shows that the error is small and has little effect on the focusing of the target. The correct echo model and CFBP algorithm are combined, which can integrate the non-start-stop model and curved orbit, further improving the imaging processing efficiency. The algorithm can be divided into seven steps: aperture division and subimage formation, first compression filter step, second compression filter step, azimuth sampling, second compression filter recovery step, first compression filter recovery step, coherent accumulation of subimage and recursion. After recursion step, the final image with full resolution can be obtained. The detailed processing flow of the proposed algorithm is introduced in detail in the paper. In addition, the simulation experiment with 0.07 m azimuth resolution and 0.21 m ground range resolution is performed to show the performance of the proposed algorithm. GF-3 data with 8.58 s synthetic aperture time are also performed to verify the effectiveness of the proposed algorithm. The azimuth-time-varying FM rates of targets in the center and edge of the scene are calculated for comparison. The focusing quality of targets is also evaluated.
The rest of this paper is organized as follows. The echo model of SAR is introduced in Section 2. The comparison of correct echo model and traditional echo model is discussed. The modified CFBP algorithm is also described in detail and the processing flow of the proposed algorithm is presented. In Section 3, the validity of the proposed method is confirmed by simulation experiment and GF-3 data experiment. In Section 4, a discussion is given to show the performance and advantages of the proposed algorithm. Section 5 concludes this paper.

Materials and Methods
In this section, the traditional echo model and proposed correct echo model are introduced in detail, followed by a comparison between the two models. Then, the modified CFBP algorithm is described and the residual error in the modified CFBP algorithm is analyzed. In addition, the processing flow of the algorithm is presented.

Echo Model of SAR
In this subsection, the traditional echo model is introduced first. Then, the correct echo model based on the non-start-stop model is derived in detail.

Traditional Echo Model
In the spotlight mode, the beam of SAR system is steered to illuminate the common spot on the surface of Earth, as Figure 1 shows. In O-xyz coordinate, original point O is located at the center of the scene. The x-axis is along the direction of mean velocity of the satellite and y-axis is along the ground range direction. Suppose the position of nth target is T n = (x n , y n , z n ) and the satellite position is S(η) = (x s (η), y s (η), z s (η)) at azimuth time η. The start-stop approximation is widely used in conventional SAR algorithms. Under this approximation, range the history of SAR signal can be expressed as R(η) = S(η) − T n , where · denotes the norm. The propagation delay τ 0 (η) = 2R(η)/c is a function of azimuth time η, where c denotes the speed of light. The echo model can be written as [34] where A 0 is a complex constant, f 0 is the carrier frequency, K r is the FM rate of the transmitted signal, w r is range envelope, and w a is azimuth envelope.

Problem
Although the start-stop approximation is widely used in SAR imaging processing, it will introduce intolerable errors in spaceborne high resolution case. As mentioned in [7], the start-stop approximation will introduce two effects, slow-time effects, and and fast-time effects. The error introduced by the slow-time effect will introduce an azimuth shift which will not influence the focusing of the target. In 2D frequency domain, the phase error introduced by fast-time effect can be expressed as [7] where f a is the azimuth frequency. f r is the range frequency with −B r /2 ≤ f r ≤ B r /2, where B r is the range bandwidth. In the edge of the 2D frequency domain (corresponding to the position where azimuth frequency and range frequency reach their maximum), the error can be given by where f a,max is the maximum of azimuth frequency, f r,max is the maximum of range frequency, and T r is the pulse width. For a SAR data with 600 MHz range bandwidth and 20,000 Hz azimuth bandwidth, the phase of error H( f a , f r ) can be seen in Figure 2a (suppose the pulse width T r is 50 µs and the corresponding K r is 1.2 × 10 13 Hz/s). The phase error error edge in the edge of the 2D frequency domain will reach 80 degrees, which should be corrected in the imaging processing. Furthermore, the change of error edge varied with T r and f a,max is shown in Figure 2b. As we can see, the phase error error edge increases with the increase of T r and f a,max . In high resolution case, error edge is large, thus, the start-stop approximation is no longer applicable.

Correct Echo Model
As mentioned above, in the high resolution spotlight mode of SAR, the motion of the satellite during the transmission and reception of the signal will introduce intolerable error [7,35]. Therefore, the non-start-stop model should be used in the high resolution spotlight mode of SAR. The propagation delay τ d in non-start-stop model is a function of the fast time τ and the azimuth time η. A refined model that takes this motion into account is defined as [35,36] where (η + τ s ) denotes the transmitting time of the radar signal, ) denotes the velocity of the satellite at azimuth time η. Assume V(η) is invariant during the transmission and reception of the signal. Therefore, Equation (4) can be expressed as where θ is given as follows Therefore, the propagation delay τ d can be derived as [17] As a result, the correct echo model can be expressed as follows For simplicity, A 0 , w r (τ) and w a (η) in Equation (11) are ignored in the subsequent analysis. Thus, Equation (11) can be derived as where K n (η) denotes the range FM rate which is varying with azimuth time η.

Comparison between the Traditional Echo Model and the Correct Echo Model
As shown in Table 1, the difference between the traditional echo model and the correct echo model can be summarized as follows: • A new range FM rate expression. K n (η) can be regarded as an azimuth-time-varying FM rate in the correct echo model, while the range FM rate K r is constant in the traditional echo model.

•
The position of the target after range compression is shifted. In the traditional echo model, the position of the target after range compression is τ 0 (η). However, in the correct echo model, the peak position is The azimuth phase is different. In the traditional echo model, the azimuth phase is 2π f 0 τ 0 (η).
However, in the correct echo model, the azimuth phase is 2π f 0

Modified Cartesian Factorized Backprojection Algorithm
In this section, the proposed modified CFBP algorithm is introduced in detail. An azimuth-time varying range FM rate is used for range compression. The residual error due to mismatch is also analyzed in detail. Then, the range history and compensation phase are derived for imaging processing. Finally, the processing flow chart of the proposed algorithm is presented.

Processing Flow of the Modified CFBP Algorithm
In the first step of the imaging processing, an azimuth-time-varying FM rate K n (η) is used for range compression. However, it can be seen from Equation (13) that K n (η) is dependent on the relative position between the target and SAR platform. As a result, K n (η) is different for different targets.
In general, in order to improve the efficiency of processing, the target T 0 in the center of the scene is selected as the reference point. Therefore, K 0 (η) calculated from the target T 0 is used for range compression. Thus, the focusing quality of the target T 0 is the best in the scene. The residual error of other targets in the scene introduced by the range compression mismatch should be analyzed. The detailed analysis can be seen in the next subsection.
After range compression, the signal can be expressed as follows where B r is the range frequency bandwidth. Transform the Equation (14) into the range frequency domain In first processing step of CFBP algorithm, the synthetic aperture is split into N subapertures. The impulse response function (IRF) at position (x i , y i , z i ) generated from nth subaperture can be written as where η n,begin and η n,end denotes the start and the end times of the nth subaperture.
To decrease the NSR of subimage, two spectrum compression filters of nth subaperture are used in CFBP algorithm, which can be expressed as follows [31,33] where (X s (η n,cen ) , Y s (η n,cen ) , Z s (η n,cen )) denotes the aperture center position of nth subaperture.
x i,cen is the center value of x i in the scene. f r,g is the subimage range frequency and its relationship with f r can be expressed as where θ is the incidence angle. (X s (η cen ) , Y s (η cen ) , Z s (η cen )) denotes the aperture center position of whole subaperture. The first spectrum compression filter F n (x i , y i ) is used to compress the 2-D wavenumber support region [33]. However, the azimuth bandwidth may still large due to the varying range frequency. Therefore, the second spectrum compression filter is performed in the range frequency domain to correct the azimuth spectrum expand [33]. The detailed derivation can be seen in [33]. After operating with two spectrum compression filters, the azimuth bandwidth of the subperture image is reduced, thus, a small NSR is enough for subaperture imaging. As a result, the computational burden is reduced greatly while operating backprojection algorithm in subaperture for subimage formation. The following step is upsampling by zero-padding in 2D frequency domain to get a high NSR of the subimage. After decompressed by two spectrum compression filters, two subimages can be added up to get a higher resolution image. After step-and-repeat operation, the final image with full resolution can be obtained. In the imaging processing, another problem that needs to be considered is how to determine the number of subaperture. In [25], the formula of length of subaperture l is given as follows where β represents the azimuth width of the beam, γ represents the azimuth oversampling ratio and |x i | R is the ground range width of a subscene. As stated in [25], the whole scene can be divided along the range direction according to Equation (20) in spotlight mode of SAR system.

Derivation of the Error Due to Mismatch
The range FM rate K n (η) is varying with azimuth time η in the correct echo model. Therefore, in the first step of the imaging processing, an azimuth-time-varying FM rate K n (η) is used for range compression. However, it can be seen from Equation (13) that K n (η) is dependent on the relative position between the target and SAR platform. As a result, K n (η) is different for different targets.
In general, in order to improve the efficiency of processing, the target T 0 in the center of the scene is selected as the reference point. Therefore, K 0 (η) calculated from the target T 0 is used for range compression. Thus, the focusing quality of the target T 0 is the best in the scene. The residual error of other targets in the scene introduced by the range compression mismatch should be analyzed.
The range FM rate difference ∆K n,0 (η) between the target T n and T 0 can be written as ∆K n,0 (η) causes a filter mismatch in range direction, which broadens the impulse response width (IRW) and increases the size of sidelobes [34]. The quadratic phase error (QPE), which characterizes the broadening of range signal after range compression, is defined as follows [34] QPE = π∆K n,0 (η) T r 2 2 (22) where T r is the pulse width of range signal. The phase error ∆φ in the peak position caused by ∆K 0 (η) can be expressed as follows [34] In the imaging processing, if the IRW broadening and phase error ∆φ are at a tolerable level, K 0 (η) calculated from the target T 0 in the center of the scene can be used for range compression. However, if the IRW broadening and phase error ∆φ are large, which would degrade the focusing quality seriously, K 0 (η) used for range compression is no longer applicable. From the simulation experiment in Section 3, we can see that the IRW broadening and phase error are small which will not influence the focus of edge targets. In addition, the phase error due to mismatch can be compensated in the processing, which can further improve the accuracy of phase in the final image.

Algorithm Steps
Based on the former investigation, the algorithm steps are summarized as follows for the convenience of implementation. The processing flow can also be seen in Figure 3. 1 Range compression by a varying frequency modulation rate K 0 (η). 2 The peak position in the compressed signal for target T n is R mod (η) = τ d,n (η) 1−α n (η) + f 0 α n (η) K n (η) . 3 The peak position in the compressed signal for target T n is R mod (η) = τ d,n (η) The compensation phase is φ mod (η) = 2π f 0 denotes the residual compensation phase due to mismatch. It should be noted that in the processing of phase compensation in step2, the phase error ∆φ com introduced by range FM rate mismatch is also compensated, which improves the phase accuracy of the imaging result. As a result, the error of IRW broadening exists in the final image. From the simulation analysis in Section 3, it can be seen that the IRW boarding error is small which will cause negligible error in the final image.

Results
To show the effectiveness of the proposed algorithm, the simulation experiment and GF-3 data experiment are performed in this section. The results are compared with the traditional echo model. The azimuth-time-varying FM rates of different targets in the scene are calculated for comparison. In addition, the corresponding phase error is demonstrated to analyze its influence on the focusing of the target.

Simulation Experiment
In this subsection, the simulation experiments are carried out to validate the proposed algorithm. The satellite orbit parameters in Table 2 are used for simulation. The height of the satellite is 600 km. The orbit eccentricity is 0.0011. The parameters of the SAR system can be seen in Table 3. The synthetic aperture time is 20.47 s. Therefore, a resolution higher than 0.1 m in azimuth direction can be reached. The pulse width is 40 µs. The motion of the platform between transmitting and receiving signals cannot be ignored due to the large pulse width. The spotlight mode SAR echo data is simulated with a 3 × 3 point target array being arranged in the simulation scene, which can be seen in Figure 4. The SAR echo is obtained by using the iteration method proposed in [13], which can simulate the continuous motion of the satellite.   In the first step, the azimuth-time-varying range FM rates of target P1, P5 and P9 are calculated, which are shown in Figure 5a. Assume ∆K P1,P5 (η) denotes the range FM rate difference between target P1 and target P5, and ∆K P9,P5 (η) denotes the range FM rate difference between target P9 and target P5. ∆K P1,P5 (η) and ∆K P9,P5 (η) are shown in Figure 5a. If the range FM rate of P5 is used for range compression, the corresponding phase error of P1 and P9 can be written as (a) (b) The change of ∆φ P1 and ∆φ P9 over azimuth time η can be seen in Figure 6a. The maximum ∆φ P1 and ∆φ P9 are less than 0.5 degrees, thus, the residual error due to mismatch in the whole scene is small. In addition, the IRW broadening of targets P1 and P9 can be seen in Figure 6b. The value of IRW broadening is less than 0.06%, which indicates that the influence of mismatch is small and the mismatch can be ignored in the analysis. After processing with the proposed algorithm, the interpolated contour of points P1~P9 can be seen in Figure 7. In addition, quantitative results of the imaging quality for P1∼P9 are listed in Table 4, including the IRW, the peak sidelobe ratio (PSLR), and the integration sidelobe ratio (ISLR) in range and azimuth directions. It can be seen that all targets have a good focusing performance, which verifies the feasibility of the proposed algorithm. In contrast, when the traditional echo model is employed, the imaging result of target P5 is shown in Figure 8 and the quantitative result of target P5 is shown in Table 5. In addition, the range and azimuth profiles of P5 by using the start-stop model and the non-start-stop model are shown in Figures 9 and 10 respectively. The IRWs of range and azimuth direction are 0.4889 m and 0.1279 m respectively, which are much higher than the theoretical values. Compared with Figure 7e, the point P5 in Figure 8b is seriously defocused, which proves that the start-stop approximation is no longer applicable in high resolution case. In addition, the computational times of the modified CFBP algorithm and traditional BPA are also compared. The computational time of modified CFBP algorithm is about 7% computational time of traditional BPA, which denotes that the modified CFBP algorithm has high efficiency.

Real Data Experiment
In this subsection, SAR raw data of spotlight mode acquired from GF-3 satellite are used to verify the effectiveness of the proposed algorithm.
The parameters of the raw data can be seen in Table 6. The carrier frequency is 5.4 GHz. The range bandwidth of the raw data is 240 MHz. The pulse duration of transmitted signal is 45 µs. Although the beam steering ability of GF-3 satellite is −1.9 • to 1.9 • , the steering range of the spotlight mode raw data is from −1.78 • to 1.78 • in this experiment. The corresponding acquisition time is 8.58 s. The theoretical azimuth resolution and ground range resolutions are 0.396 m and 0.889 m, respectively. In the processing of the raw data, the azimuth-time-varying FM rate is calculated for range compression. K 0 (η) of point target in the center of the scene is shown in Figure 11a. The standard FM rate K r is 5.33333 × 10 12 Hz/s. The ratio of K 0 (η) to K r is shown in Figure 11b. The maximum of the ratio is about 1.000003. By using the proposed algorithm, the imaging result of GF-3 spotlight data can be seen in Figure 12. For comparison, the imaging result of local area A based on start-stop model is shown in Figure 13a. Compared with image in Figure 13a, the image in Figure 13b has a better focusing quality. In addition, the contour of one target with high reflectivity in the scene is shown in Figure 14. It can be seen that based on the non-start-stop model, the target has a high focusing quality, which proves the effectiveness of the proposed algorithm.

Discussion
Nowadays, the spaceborne SAR system is characterized by high resolution and swath coverage. The high resolution mode of SAR system brings new difficulties and challenges. In this paper, a modified CFBP algorithm integrating with non-start-stop model is proposed for high resolution SAR imaging. The simulation experiment with 0.21 m range resolution and 0.07 m azimuth resolution demonstrates that the proposed algorithm can achieve high focusing quality. The GF-3 spotlight mode with 0.889 m range resolution and 0.396 m azimuth resolution is also used to verify the effectiveness of proposed algorithm. Compared with the tradition echo model, the proposed correct echo model considers the motion of the satellite during the transmission and reception of the SAR signal. A new range FM rate is derived, which is varied with azimuth time. The position of the target after range compression and the azimuth phase are also derived. The proposed echo model can reflect the true motion of satellite; therefore, it has high accuracy and can be applied to high resolution SAR mode case. The CFBP algorithm proposed in [25,26,[31][32][33] is an efficient algorithm. In order to improve the efficiency of imaging processing, the correct echo model and CFBP algorithm are combined. In addition, the residual error due to mismatch is also derived. The results show that the residual error is small which can be ignored in the analysis.
The advantages of the proposed algorithm are as follows. The curved orbit and non-start-stop mode are all considered in the correct echo model, therefore, it has high accuracy. The whole aperture is divided into subapertures and the final image is obtained by coherent combination of subimages, which can improve the efficiency. Although in the subimage formation, the BPA is used which may increase computational burden; however, this step can be done by using parallel computing tool which can greatly increase efficiency [19,24,37].

Conclusions
The high resolution spotlight mode of SAR plays an increasingly important role in remote sensing. In this paper, a modified CFBP algorithm integrating with non-start-stop model is proposed for very high resolution spotlight mode of SAR data processing. In the high resolution case, the start-stop model will introduce errors, which will cause the target to be defocused in the image. In addition, under the condition of long synthetic aperture time, the error introduced by curved orbit should also be considered, further increasing the difficulty of imaging processing.
In this paper, an accurate correct echo model is used to depict the SAR signal with high accuracy. The non-start-stop model and curved orbit are all considered in the correct echo model. Then, The correct echo model and CFBP algorithm are combined to improve the accuracy and efficiency of imaging processing. Based on the echo model, at first, an azimuth-time-varying range FM rate of the center target is used for range compression. The residual error of the proposed algorithm is also discussed in detail. Then, the range history and compensation phase are derived for imaging processing. The modified CFBP algorithm is divided into seven steps. After all the steps are complete, the final image with full resolution is obtained. In the end, the simulation experiment and GF-3 data experiment demonstrate the feasibility of the proposed algorithm. What is more, the high resolution mode is sensitive to the topographic change. Therefore, combing the imaging algorithm and topographic change is a new challenge. In the future, further research will be carried out on this problem.

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