High-Precision Spotlight SAR Imaging with Joint Modified Fast Factorized Backprojection and Autofocus

Fast factorized backprojection (FFBP) takes advantage of high accuracy of time-domain algorithms while also possessing high efficiency comparable with conventional frequency domain algorithms. When phase errors need to be compensated for high-resolution synthetic aperture radar (SAR) imaging, however, neither polar formatted subimages within FFBP flow nor the final Cartesian image formed by FFBP is suitable for phase gradient autofocus (PGA). +is is because these kinds of images are not capable of providing PGA with a clear Fourier transform relationship (FTR) between image domain and range-compressed phase history domain. In this paper, we make some essential modifications to the original FFBP and present a scheme to incorporate overlapped-subaperture frame for an accurate PGA processing. +e raw data collected by an airborne high-resolution spotlight SAR are used to demonstrate the performance of this algorithm.


Introduction
Backprojection (BP) is essentially a method of phasedarray beamforming in the time domain. It is capable of handling an arbitrary flight path and producing synthetic aperture radar (SAR) images free of geometric distortions and defocus effects [1]. Despite of these advantages, the direct BP algorithm is still difficult to put in practice due to its low efficiency. Some studies about polar formatting method for spotlight-mode SAR imaging are conducted to improve the computational efficiency [2]. During the past decade, a number of fast BP algorithms [3][4][5][6][7] for accelerating BP integral have been proposed, and fast factorized BP (FFBP) is the most representative. Based on recursive partitioning, FFBP consists of a series of stages. A major characteristic of FFBP is that subaperture length increases successively, whereas the beam spacing decreases. Finer angular resolution subimages are formed in local polar coordinates by combining coarse subimages from the previous stage. With high accuracy and efficiency, FFBP takes the best from both time and frequency domain algorithms. erefore, SAR imaging via FFBP has been attracting more and more attention in recent years.
During airborne SAR imaging, air turbulence inevitability induces platform deviation from ideal trajectory. However, it is noteworthy that even the most sophisticated motion measurement system may be not sufficient to ensure diffraction-limited SAR images. erefore, autofocus techniques are required for further processing, which is studied in [8], but research studies about autofocus combined with BP or FFBP are not comprehensive. Different from the methods based on varying antenna path parameters [9] and maximizing image sharpness [10][11][12][13], our primary aim is to achieve an accurate autofocus processing within FFBP by using current autofocus techniques, such as mapdrift (MD) autofocus, minimum-entropy autofocus, and phase gradient autofocus (PGA). In [14], basic theories about the use of autofocus were developed and a scheme of incorporating multiple aperture mapdrift (MAMD) autofocus into FFBP was presented. As a promising start, we firmly believe FFBP is capable of working with other autofocus techniques and can be extended to satisfy our growing need in much more complex scenarios. Generally, phase errors are predominantly low-order phase errors introduced by deviation from ideal trajectory, while mechanical vibration and random propagation will inevitably induce highorder phase errors [15]. However, MAMD and other parametric model-based algorithms all predefine a loworder explicit expression for phase errors and neglect the high-order components. In contrast, PGA has been verified to be a robust technique and is being widely used throughout the SAR community [16]. As a nonparametric algorithm, PGA imposes no restriction on the form of phase errors, so it can achieve excellent restorations under more general conditions. However, directly incorporating PGA into the time-domain algorithms is not easy. To our best knowledge, [17] is the only work considering that PGA is used to refocus SAR images formed by the direct BP. As revealed in [18], there are some significant factors which should be taken into account. Firstly, the spread of azimuth impulse response function (IRF) do not lie in the horizontal dimension of a backprojected image. Secondly, the Fourier transform relationship (FTR) between the image and its corresponding range-compressed phase history is not sufficiently explicit for a backprojected image.
In this paper, we propose a scheme to incorporate PGA into FFBP and achieve accurate autofocus processing for high-resolution spotlight SAR imaging. e main contributions of this work are twofold: (1) Pseudopolar coordinate system is employed as the imaging plane, and an analytic expression for azimuth IRF in a backprojected image is derived. Actually, a line-of-sight (LOS) pseudopolar coordinate system not only ensures the spread of azimuth IRF along the horizontal dimension but also provides an approximate FTR for the use of PGA. (2) Subaperture phase errors estimated by PGA cannot be directly applied to correct the range-compressed phase history data, otherwise they will cause relative shifts among subimages. To achieve an accurate phase correction, overlapped-subaperture frame (OSF) is introduced into FFBP to obtain a full-aperture phase error by coherently combining subaperture phase errors. As phase correction is recursively performed, phase estimation tends to converge and the degraded image becomes refocused.

Modified FFBP for Focusing Spotlight SAR Data
In the original FFBP, both polar formatted subimages and the final Cartesian image are formed, but neither of them can provide PGA with a FTR between image domain and range-compressed phase history domain. In order to establish a FTR, we focus on theoretical analysis and make some significant modifications to the original FFBP in this section.

Azimuth IRF in the Pseudopolar Coordinate
System. An ideal imaging geometry of spotlight SAR is shown in Figure 1; the model is based on the assumption of flat land. e platform moves along the x-axis with velocity v and generates a synthetic aperture of length L. During the data acquisition, the sensor continuously steers its antenna beam onto the scene center C. Taking the aperture center O as the origin of the polar coordinate system, the coordinates of point target P is (r p , θ p ), where r p is the range measured from O to P and θ p is the bearing angle measured away from the broadside direction. For the sake of clarity, we define a pseudopolar coordinate system described by r and Θ, where Θ is the sine value of θ, namely, Θ � sin θ. e slant range from the sensor to target P is given by where X is aperture position and Θ p is equal to sin θ p . e sensor transmits linear frequency-modulated (LFM) signal with bandwidth B and carrier frequency f c . Considering target P with a unit reflectivity coefficient, the baseband signal after range compression is expressed as where τ is the fast time, W a represents a common antenna pattern, and τ p � 2R(r p , Θ p ; X)/c corresponds to the roundway time delay. Equation (2) ignores antenna gain and 1/R 2 propagation decay factor in the amplitude term.
BP is essentially an integral process performed in the time domain. For a given aperture position X 0 , s MF (τ; X 0 ) needs to be backprojected to every pixel in the desired imaging grid, and each pixel is assigned a value according to the time delay corresponding to the range between the sensor and the pixel; for all available aperture positions, BP integral is achieved by coherent accumulation along the range history of each pixel. Importantly, an enormous advantage of BP is that it allows for image reconstruction onto an arbitrary imaging grid. In [10,17], a range and sine of the bearing angle grid is employed as an imaging plane. In this paper, we also adopt this imaging grid and name it pseudopolar coordinate system.
Due to the limited Doppler bandwidth of spotlight SAR, azimuth IRF cannot be a Dirac function which is just located at one angular cell. For target P, the peak of its azimuth IRF is concentrated at the Θ p angular cell and the side lobes spread out in the azimuth direction at the r p range cell. Without loss of generality, the contribution of target P to a pixel located at (r p , Θ) in the pseudopolar coordinate system is given by 2 International Journal of Antennas and Propagation where ΔR(r p , Θ; X) � R(r p , Θ p ; X) − R(r p , Θ; X) is the differential range between R(r p , Θ p ; X) and R(r p , Θ; X). During the derivation of equation (3), we introduce a general approximation that antenna pattern W a (X) is modeled as a rectangular window function, i.e., rect(X/L) � 1 − (L/2) ≤ X < (L/2) 0 elsewhere . It should be emphasized that h(r p , Θ) is not an image as normally formed by BP integral, but an impulse response for the pixels located at the r p range bin. When Θ � Θ p , the main lobe of the azimuth IRF is located at the Θ p angular cell. Despite all that it is not hard to explain the relationship between a reconstructed image and its pixels. Given impulse responses for all the pixels, theoretically, the image could be reached by To derive an analytic expression of azimuth IRF, ΔR(r p , Θ; X) is extended into the third-order Taylor series, i.e., By neglecting the X 2 and other higher-order terms in ΔR(r p , Θ; X), we can rewrite equation (3) and obtain a compact expression, which is given by where K Θ represents angular wavenumber. In equation (5), a generalized Fourier transform is approximately developed between K Θ and X. It states that the azimuth impulse response along X in the r p range bin is achieved by the Fourier transform of the range-compressed data, which is critical to the use of autofocus techniques.
In the following, we analyze the constraint which needs to be fulfilled for the validity of equation (5). In [14], it assumes that the majority of energy of azimuth IRF is concentrated with a small region size corresponding to the angular resolution ΔΘ � λ/2L. e contributions of higherorder terms take minimal effect on the BP integral if quadric phase error (QPE) induced by X 2 term is less than π/16 in the neighborhood of Θ p , namely, Noting that ΔΘ keeps positive, we can determine QPE as follows: By simplifying equation (6)fd6, we can obtain As equation (7) is the constraint of equation (5), if the pixel close to the edge of the imaging grid can satisfy equation (7), all the pixels on this grid can meet this constraint, too. erefore, the final constraint is given by where Θ max is the maximum value of Θ and r min is the minimum value of r. For spotlight SAR, the angular range of pseudopolar coordinates satisfies |Θ| ≤ sin(θ BW /2), where θ BW is antenna azimuth beamwidth. In this case, equation (8) can be straightforwardly written as For the sake of clarity, let us take the X-band spotlight SAR used in the experiments in Section 4 as an example, where L � 1.17 km, r min � 10.1 km, and θ BW ≈ 6.78 degrees. By substituting these parameters into equation (9), 8L · sin(θ BW /2) + 2λ ≈ 553.72 m is less than r min , which satisfies the constraint of equation (9). erefore, the approximation of equation (5) keeps valid and is precise enough to ensure the validity of FTR. Besides, wavelength λ is relatively small compared with L and r min , which can be removed from equations (8) and (9).

e Modified FFBP Algorithm.
Original FFBP splits up BP integral into an infinite sum over finite subapertures, so reconstructing each subimage can be considered as an individual direct BP. Motivated by Section 2.1, we can introduce the pseudopolar coordinate system into the original FFBP instead of the polar coordinate system and refer to this new algorithm as the modified FFBP.
Based on aperture factorization and recursive fusion, FFBP consists of a series of stages. Assume that there are K subapertures in the gth stage, where g � 1, 2, . . . , G and G is the total number of stages. K subapertures generate K subimages. Each subimage is created onto a local polar grid with the origin located at its corresponding subaperture center. Noting the differences among these subimages, we use (r, Θ (g) k ) to describe the kth polar grid in the gth stage. According to Section 2.1, the kth subimage I (g) k can be expressed as where K Θ (g) k represents the angular wavenumber corresponding to the kth subaperture.
Based on the definition of Doppler f a , K Θ (k) is in proportion to f a with a scale factor of 2π/v. In this sense, azimuth focusing of I (g) k is performed in the Doppler domain, which is equivalent to the azimuth focusing domain of SPECAN imaging [18]. Although the modified FFBP implements BP integral in the time domain, by using pseudopolar coordinates, subimages could be considered as range-Doppler images. Moreover, an implied demand of autofocus is that the spread of azimuth IRF should be along the horizontal dimension of an image. e solution to this problem is simple, but it is not categorically stated in [4,5]. By learning experience from polar format algorithm (PFA), we find that LOS polar interpolation can rotate the orientation of the imaged scene while keeping the spread of azimuth IRF along the horizontal dimension [15]. In this paper, all local pseudopolar coordinates in the modified FFBP are established along the direction of the LOS from each subaperture center to the scene center. In this way, we can obtain nonsquinted data spectrum and ensure the orthogonality between its range and azimuth dimension so that we can keep the azimuth IRF spreading out in the azimuth direction. Figure 2 shows the last two subimages of the modified FFBP. Figures 2(a) and 2(b) represent exactly the same scenes (five point targets arranged in a cross line) and correspond to the last two subapertures. Because of LOS coordinates, the imaged scenes are rotated through a slight angle in opposite directions as the red dashed lines outline. However, it will be found that the spread of azimuth impulse response for all targets lies in the horizontal dimension.
It should be noted that the LOS polar grid used in the modified FFBP is a polar grid in the image domain, while PFA uses a polar grid in the frequency domain. e modified FFBP inherently accommodates wavefront curvature and elevation-induced defocus effects arising from out-of-plane motion of platform [17]. erefore, the modified FFBP is capable of producing SAR images without geometric distortions and defocus effects. By contrast, PFA has a limited depth of focus. If we attempt to produce distortion-free SAR image by using PFA, several postprocessing corrections need to be applied [15].

Overlapped-Subaperture Frames and PGA Processing
FTR is what we emphasize in this paper, but it is not an inherence for time-domain algorithms after all. In Section 2, we make some necessary modifications to the original FFBP by introducing LOS pseudopolar coordinates. Noting that all the subimages created by the modified FFBP can satisfy the FTR requirement of PGA, we are able to employ PGA to extract subaperture phase errors from subimages. However, here emerges another issue: subaperture phase errors cannot be directly used to refocus the degraded subimages, otherwise they will cause relative shifts among the subimages. is fact arises for the reason that the linear component of a subaperture phase error has been removed by "circular shifting" operation [16]. erefore, it is necessary to retrieve the linear component belonging to each subaperture phase error before phase correction. How to achieve this goal is the main topic of this section.
In [19], overlapped-subaperture frame (OSF) is introduced to build a bridge connecting original subaperture phase errors and a full-aperture phase error. Each desired subaperture phase error can be updated with its linear component according to the full-aperture phase error function. Enlightened by this thought, this section includes two parts: building OSF and implementing PGA processing within the recursion of the modified FFBP.

Overlapped-Subaperture Frames.
Overlapped-subaperture technology was first combined with PFA to correct spatiallyvariant phase errors in spotlight SAR imaging [20]. Since then, it has been widely used to recover continuous fullaperture phase errors or estimate range-dependent phase errors in wide-swath stripmap SAR. Compared with the manner of building overlapped subapertures in [20], the modified FFBP should take a little more specific details into consideration. For a clear interpretation, a counter example that the once-through designed OSF fails our need is first given as follows.
In the first stage of the modified FFBP, the whole aperture is divided into several overlapped subapertures. Assuming that each subaperture contains 2N samples and the overlap ratio is (1/2), there should be N overlapped samples between two adjacent subapertures. In the second stage, each new subaperture will contain 2 × 2N − N � 3N samples after aperture fusion. Because the number of the samples shared by two adjacent subapertures remains N, and the overlap ratio is 1/3 at this time. As the processing goes on, the length of subaperture increases, but the overlap ratio keeps descending, which is unfavorable to a coherent combination of subaperture phase errors.
To recover a reliable full-aperture phase error function, obviously, we need to guarantee that the overlap ratio of two adjacent subapertures is big enough. In the following, we will present an innovative method of founding OSF within the modified FFBP. To illustrate this method, here we take a four-step processing as an example for analysis, and its critical steps are as follows: Step 1: like the original FFBP, the modified FFBP performs the identical processing until 8 nonoverlapped subapertures are left, and each one contains N samples. is condition is described as the first row (from top to bottom) in Figure 3(a), where the blue line segment represents subaperture and the red circle denotes subaperture center.
Step 2: perform aperture fusion among any two adjacent subapertures, and then 7 subapertures can be formed from 8 previous subapertures. As a consequence, the subaperture length is extended to 2N and the overlap ratio is 1/2. Compared with the second row in Figure 3(a), the red line segments in the second row in Figure 3(b) are newly formed by OSF.
Step 3: employ PGA to extract subaperture phase errors from the 7 subimages formed by Step 2, and the detailed description will be given in the next section. It should be emphasized that, of the 7 new subapertures, only the odd ones, namely, 4 nonoverlapped subapertures will be transported to the next stage for aperture fusion, which are shown as the blue line segments in the second row of Figure 3(b). e other 3 subapertures shown as the red ones will be directly deleted after phase estimation. International Journal of Antennas and Propagation Step 4: repeat Steps 2∼3 and then the last two subapertures are generated. In this situation, there is no need to build OSF, but should carry out the same processing as the original FFBP does in the last stage.
An advantage of this method is that the overlap ratio being 1/2 keeps unchanged in spite of different subaperture lengths at different stages. OSF for phase correction is not required in a stage-by-stage manner for the reason that PGA is instinctively a recursive process. erefore, building OSF in two or three successive stages is enough for practical use.

PGA Processing. "
Original subaperture phase errors-Full-aperture phase error-Desired subaperture phase errors" is a direct way to achieve an effective phase correction. Designed to obtain a reliable full-aperture phase error, OSF is an essential link in this chain. With LOS pseudopolar coordinates and OSF, preparations for PGA processing are all complete.
Assume that phase correction is not performed in the modified FFBP until the gth stage. Assuming the influence of phase errors becomes obvious on degrading image quality in this stage, we need to build OSF and perform an effective PGA processing. Assume that there are K subapertures in the beginning of gth stage. According to Section 3.1, K subapertures can generate K − 1 overlapped subapertures with the overlap ratio being 1/2 and each subaperture corresponds a single subimage. Based on the FTR derived in Section 2, the range-compressed phase history of I (g) k is given by where K Θ (g) k and X e value range limitation of K Θ (g) k is given by Supposing the bulk of motion errors have been effectively corrected by the high-precision motion measurement system, such as global positioning system (GPS) and/or inertial navigation unit (IMU), migration caused by the residual errors is less than one range cell. In this section, the sign ϕ k denotes the phase error retained in the real range-compressed phase history H ′ (g) k , and ϕ (g) k is a function of subaperture position X (g) k . erefore, the relationship between the theoretical value H (g) k and the real value H ′ (g) k could be described as Based on equation (13), we can employ PGA [16] or weight PGA (WPGA) [21] to extract subaperture phase error from H ′ (k) . As a derivative of PGA, WPGA utilizes the signal-toclutter ratio (SCR) to compute the weights of range cells. Large weights are applied to range cells with high SCRs. Assume that WPGA has been employed to each subimage, and the estimated subaperture phase errors are denoted by ϕ (g) K− 1 , where "∧" represents the estimated value of a variable. Compared with the real values, the most difference of these K − 1 estimated values is that the linear components reflecting the interrelation of subapertures have been removed by "circular shifting" operation in WPGA. According to the previous analysis, we need to update each subaperture phase error with its linear phase from a full-aperture phase error function.
To recover a full-aperture phase error function, subaperture phase estimates are coherently combined by eliminating the unknown linear difference between the phase errors extracted from the overlapped samples of the neighboring subapertures, which are illustrated in Figure 4. In practice, the subaperture phase estimates fluctuate sharply and the precision degrades with the increase of clutter and high frequency noise. For a smooth full-aperture phase error function, the phase filter is usually required. After a coherent combination and a smooth filtering, the full-aperture phase error function ϕ (g) (X) is recovered. Noting that each subaperture has a definite subaperture span X (k) , k � 1, 2, . . . , K − 1, we can recalculate subaperture phase errors from ϕ (g) (X). Of these K − 1 subimages, the even ones (totally (K/2) − 1) are no longer needed after phase estimation and deleted from memory, and the other K/2 subimages should be refocused before they are transported to the g + 1th stage. So, we perform motion compensation to these K/2 subimages with desired subaperture phase errors, including residual range cell migration (RCM) correction and phase error removal. Until now, we complete all PGA processing in the gth stage. e detailed illustration is shown in Figure 5. If it is necessary to perform phase error correction in the g + 1th stage, we just repeat the operations mentioned above.
Profit from OSF scheme, PGA processing can be performed in a recursive manner and compatibly blended into the modified FFBP. Undoubtedly, OSF introduces computational burden, but only two or three successive stages are enough for PGA to refocus a full-resolution image up to an acceptable level. It can be regarded as a tradeoff between the image quality and the processing efficiency. Moreover, the "windowing" step of PGA in this paper should be emphasized that the maximum window width becomes small after the first applying PGA, which is of benefit to the rapid convergence of phase estimation.

Experiments
In this section, the raw data collected by an airborne highresolution spotlight SAR are used to demonstrate the performance of this algorithm. e main parameters of the SAR system are tabulated in Table 1. In the experiment, the synthetic aperture time is 7.8 seconds and each pulse contains 8192 range bins. e image is 0.8 km×1.2 km in size with range and azimuth resolutions of 0.15 m×0.15 m (angular resolution varies in the range direction and 0.15 m is the resolution corresponding to the closest slant range). Empirically, when the number of samples is 2048 in each subaperture, angular resolution is 1.2 m and the influence of phase errors emerges gradually. en, we begin to build OSF and perform PGA processing to coarse subimages. Although phase correction has been performed in the current stage, fractional phase errors may still exist in subimages. Due to the fact that angular resolution increases to 0.6 m in the next stage, the influence of residual errors will be magnified. erefore, we also need to perform PGA processing to correct residual errors in the next stage. On the contrary, the signal-to-noise ratio of the subimages will be higher as the subaperture becomes longer, which also enhances the performance of PGA. e processed results and estimated motion error are presented in Figure 6. Due to the existence of motion errors, Figure 6(a) without autofocus involves in serious image blurs, while Figure 6(b) processed by the proposed method is well-focused over the whole scene. Radial motion error estimated by WGPA is observed in Figure 6(c). In order to assess the performance of the algorithm, subscenes A-C from Figure 6 are magnified in Figure 7. For each pair of subscene images in Figure 6, the upper blurred image is extracted from Figure 6(a), while the lower one is obtained by the proposed method. Obviously, all subscene images in Figure 7 generated by the proposed method are well-focused, such as the shadow of the tree and the corner reflector array aligned in "T" form. In Scene C, two close-located corner reflectors outlined by circle "1" are arranged with 0.15 m distance in azimuth, and another two outlined by circle "2" are placed with 0.15 m distance in range. To demonstrate the performance of the proposed method, we give imaging Range-compressed phase history data Recover full aperture phase error with linear phase remove Overlapped subaperture 0 X (1) X (2) X (3) X (4) X (K−1) X ϕ (1) ϕ (2) ϕ (3) ϕ (4) ϕ (K-1) ϕ (X) Figure 4: Full-aperture phase error combination.    Figure 8, where one can clearly find that the depth of notch between the two peaks is lower than − 30 dB and the side lobes of IRFs are almost below − 20 dB. is fact provides a clear evidence to show the great discrimination of our proposed method. To evaluate the image quality, both range    and azimuth IRFs of a single reflector outlined by circle "3" are given in Figure 9. e level of side lobes below − 20 dB also gives us a compelling reason to believe the effectiveness of our proposed method. From Figure 6(c), the radial error is larger than one range resolution cell, which is beyond the assumption that migration caused by the residual errors should be usually less than one range cell. Fortunately, RCM can be corrected in the range frequency domain during FFBP imaging. rough the results above, this skill has succeeded in image processing and RCM correction. It makes us to believe that some other skills or ideas in frequency domain algorithms may be adopted in the modified FFBP, such as the correction of range-dependent phase errors in wide-swath SAR images [22]. Another set of results is presented below. e main parameters of the SAR system are tabulated in Table 2. e processed results and estimated motion error are presented in Figure 10. Figure 10(a) is blurred with motion errors and shown in Figure 10(c), while Figure 10(b) shows the SAR image restored by PGA. Magnified images of subscenes A-C in Figure 10 are given in Figure 11 such that the results processed by the proposed method are well-focused. To evaluate the image quality, both range and azimuth IRFs of a single reflector in Scene C is shown in Figure 12, and the peak side lobes of the selected reflector are less than 20 dB.

Conclusion
In this paper, we incorporate PGA into the modified FFBP to achieve an accurate autofocus for high-resolution spotlight SAR imaging and investigate its performance in detail. In this method, LOS pseudopolar subimages are reconstructed to provide PGA with a FTR between image domain and its corresponding range-compressed phase history domain. Based on OSF, full-aperture phase error function can be obtained to retrieve the linear phase belonging to each subaperture phase error. Within the recursion of the modified FFBP, phase estimation and correction can be achieved at a high precision and efficiency. It should be noted that the unknown DEM will cause the degradation of the focusing and positioning performance of the method in this paper so that we will study in the future work.
Data Availability e data used in this article are obtained from a private data source.

Conflicts of Interest
e authors declare that they have no conflicts of interest.