Feature Preserving Autofocus Algorithm for Phase Error Correction of SAR Images

Autofocus is an essential technique for airborne synthetic aperture radar (SAR) imaging to correct phase errors mainly due to unexpected motion error. There are several well-known conventional autofocus methods such as phase gradient autofocus (PGA) and minimum entropy (ME). Although these methods are still widely used for various SAR applications, each method has drawbacks such as limited bandwidth of estimation, low convergence rate, huge computation burden, etc. In this paper, feature preserving autofocus (FPA) algorithm is newly proposed. The algorithm is based on the minimization of the cost function containing a regularization term. The algorithm is designed for postprocessing purpose, which is different from the existing regularization-based algorithms such as sparsity-driven autofocus (SDA). This difference makes the proposed method far more straightforward and efficient than those existing algorithms. The experimental results show that the proposed algorithm achieves better performance, convergence, and robustness than the existing postprocessing autofocus algorithms.


Introduction
Synthetic aperture radar (SAR) can create high cross-range resolution images through coherent processing of the returned pulses received at different antenna position due to a moving platform, which enables provision for the effect of a large virtual aperture size [1][2][3]. To achieve full performance of SAR, the exact flight trajectory of the platform should be provided to the signal processor for proper motion compensation. Especially in the case of airborne SAR, the effect of the unexpected platform motion due to the factors such as wind gusts and aircraft vibrations should be compensated before SAR processing. These compensations for airborne SAR are usually performed by using the measured navigation data such as the output of global positioning system (GPS), inertial navigation system (INS), or embedded GPS/INS (EGI) [4]. However, these data also contain the measurement errors due to the inaccuracy of the navigation sensors, and these residual errors cause the phase errors of SAR data, which degrades the quality of the SAR image. There are several methods to solve this kind of image quality degradation, which are usually called autofocus.
Phase gradient autofocus (PGA) is one of the most widely used algorithms to estimate the phase error [5]. The method assumes all the complex reflectivity in the image windowed with appropriate window size, except the center-shifted point target at each range bin, are distributed as zero-mean Gaussian random noises. These assumptions limit the performance of the algorithm in spite of its robustness and fast convergence. The point target and random noise assumption are inappropriate to the scene that contains dominant targets close to each other in the azimuth direction. Quality phase gradient autofocus (QPGA) and generalized phase gradient autofocus (GPGA) are proposed to alleviate this problem [6,7]. However, both QPGA and GPGA use the limited window size, which also limits the bandwidth of the estimated phase error. Therefore, these PGA-based methods

Phase Error Correction in SAR Images
The typical SAR image formations are achieved through the two-dimensional inverse discrete Fourier transform (DFT) of the motion compensated echo signal formatted to Cartesian coordinates. If we define the range-compressed, motion compensated echo signal as G(n, m), where n and m are the indices of range cells of the scene and received pulses, respectively, the complex image g(n, k) is formulated as where k is the indices of azimuth cells of the scene, M is the number of received pulses, and j 2 = −1. The residual phase errors after the motion compensation in the SAR image are mainly due to the errors in the navigation data, and these errors corrupt each returned pulse signal phase. Therefore, most of the phase distortions are in azimuth direction, and the phase-corrected complex image can be formulated aŝ where ϕ(m) is the estimated residual phase error. If we can obtain an accurate estimation of the phase error, the corrected image is also obtained through (2), which can be carried out through FFT.

Phase Adjustment for Desired Image
If we have the desired complex image and want to obtain the phase adjustment for minimizing the error between the desired and phase adjusted image, we can formulate the following cost function: where h(n, k) is the desired complex image and N is the number of range cells of the scene; (·) * and Re(·) denotes complex conjugate and the real part of the complex number, respectively. We assumed that the number of azimuth cells of the scene is equal to M. Let the derivative with respect to ϕ(m) be zero to minimize the cost function in (3): The second and third additive terms in (3) are ignored because the phase does not have an effect to the variation of the total intensity. The derivative ofĝ(n, k) with respect to ϕ(m) is derived from (2) as given in the following equation: Then, (4) becomes where Im(·) denotes the imaginary part of the complex number and H(n, m) is defined as H(n, m) in (7) can be computed through FFT. From (6), the sufficient solution for ϕ(m) is obtained from the following equation: where ∠(·) denotes the phase of the complex number. The result in (8) will be used in Section 3 to explain the proposed method.

Brief Review of Regularized Reconstruction
SAR image formation is basically solving an inverse problem to obtain the image from the returned signal. If we know the exact SAR observation model, SAR imaging is identical to solving the following equation with respect to x [25]: where x, y, and v comprise an MN × 1 vector, which is the column stacked version of the complex image, returned echo, and measurement noise, respectively, and C is the discretized observation kernel. The following cost function is usually used to solve the inverse problem in (9) [25]: where · p denotes the l p -norm and λ is a weighting parameter for regularization term. The first term represents the fidelity term for preserving the dependence on the observation model and the second term is the regularization term for enhancement of the feature of the image. The typical choice for p is 1 [26]. The optimization can be achieved through two ways, one is the gradient descent algorithm [26] and the other is the iterative shrinkage thresholding algorithm (ISTA) [28]. For p = 1, ISTA is represented as where α is a sufficiently small positive number, x i is estimate of x in i-th iteration, and S λ (·) is a soft-thresholding function defined as Iteration in (11) achieves the optimal image for (10) if α is smaller than the largest eigenvalue of C H C. Unlike the inverse problem of (9) and (10), if we already have the complex image formatted by (1) and want to carry out the feature enhancement in postprocessing, x and y become the column stacked version of the complex image before and after the postprocessing, respectively, and kernel C becomes the identity matrix. Then, (10) is rewritten as where g(n, k) is the postprocessed complex image. Unlike the iterative solution in (11), the optimization of (13) can be achieved by the following the closed-form solution [27]: If we approximate |g(n, k)| to be [g * (n, k)g(n, k) + ε] 1 2 for a sufficiently small positive scalar ε, then it satisfies ∂J ∂g(n, k) g(n,k)=S λ (g(n, k)) ≈ 0 From (14), each element of g whose absolute value is smaller than λ is removed, and the absolute value of the other elements are decrease by λ. This so-called "denoising" [28] is an important procedure for the proposed algorithm described in the next section.

Proposed Method
In this section, we define a cost function and derive the equation for its optimality in an indirect optimization manner. We then propose an algorithm to obtain the solution for the equation.

Cost Function and Its Minimization
The phase-corrupted image, owing to the residual phase error, shows distorted features such as degraded resolution and unexpected side-lobes, and hence the goal for phase estimation can be interpreted as to remove these distortions. Therefore, the objective of our method is a phase adjustment to achieve the feature enhancement, which can be realized by concurrent optimization of (3) and (13), i.e., the following cost function: whereĝ is the function of ϕ according to (2). The minimization of (16) enables obtaining the feature-enhanced image g and the phase adjustment ϕ to make the phase-corrected image close to g. We can achieve the minimization through the following indirect optimization. Let g be the solution of the necessary condition for the minimum of (16), i.e., ∂J/∂g = 0; then, the solution is obtained from (14) as If (17) holds for every ϕ, then we can replace the cost function in (16) by the following indirect cost function: The derivative of (18) with respect to ϕ is because ∂J/∂g ≈ 0 for g in (17). Let (19) to be zero for minimization of V( ϕ), then the condition for optimality can be derived by using the result in (6) as given in the following equation: where G(n, m) is defined as The solution of (20) is the final phase error estimate of our proposed method.

Algorithm
The sufficient condition for (20) can be written in a manner similar to that in (8): Equation (22) is not a closed-form solution of ϕ(m) unlike (8), because G(n, m) is still the function of ϕ(m). This implicit equation can be solved by the fixed-point iteration method, which has been applied to solve the similar problems for ME [13,14]. The flowchart of the proposed algorithm is represented in Figure 1. The computation of (2) and (21) in the algorithm can be carried out through FFT, which enables a fast iteration.
In every iteration, g i takes the same role of the desired complex image, i.e., h in (7), which is the denoised image ofĝ i . Because of the soft-thresholding function in (17), a few elements ofĝ i whose absolute value is larger than λ are preserved in g i and the others are set to be zero. We call these preserved elements "features" in this paper. These remained features are the reference for the phase estimation in every iteration. In other words, the phase adjustment is achieved to fit the given image to the feature-preserved image. Therefore, we refer to the proposed algorithm as "feature preserving autofocus (FPA)" in this paper. Since the existing regularization-based autofocus algorithms such as SDA deal with the returned pulse data before processing, image formation algorithms such as the gradient descent algorithm and ISTA described in Section 2.3 are required for every iteration. Therefore, these algorithms require huge computation for a large scene size. Meanwhile, the proposed FPA requires only one soft-thresholding to obtain the reference image for phase adjustment. This simple iteration sequence makes the proposed algorithm more suitable for on-board SAR processing. In every iteration, i g takes the same role of the desired complex image, i.e., h in (7), which is the denoised image of i ĝ . Because of the soft-thresholding function in (17), a few elements of i ĝ whose absolute value is larger than λ are preserved in i g and the others are set to be zero. We call these preserved elements "features" in this paper. These remained features are the reference for the phase estimation in every iteration. In other words, the phase adjustment is achieved to fit the given image to the feature-preserved image. Therefore, we refer to the proposed algorithm as "feature preserving autofocus (FPA)" in this paper. Since the existing regularization-based autofocus algorithms such as SDA deal with the returned pulse data before processing, image formation algorithms such as the gradient descent algorithm and ISTA described in Section 2.3 are required for every iteration. Therefore, these algorithms require huge computation for a large scene size. Meanwhile, the proposed FPA requires only one soft-thresholding to obtain the reference image for phase adjustment. This simple iteration sequence makes the proposed algorithm more suitable for on-board SAR processing.

Selection of the Threshold
It is worth noting that a careful choice of the threshold, i.e., λ , is recommended for the proposed FPA. A choice of small λ enables the features to contain a sufficiently large number of scatterers, which achieves accurate estimation. However, if λ is too small, the remaining features would also contain the artifacts of the SAR image such as side-lobes even though their magnitude are reduced. Thus, it may cause low speed of convergence, and insufficient estimation accuracy if the artifacts still remain in the features at the final iteration. On the other hand, if λ is too large, some dominant scatterers for phase estimation would be removed from the features. In this case, the estimation is too concentrated in a few remained strong scatterers, which leads to local optimal estimation even if it enables to achieve faster convergence rate.
To solve these tradeoffs, we propose varying the threshold to satisfy both fast convergence and optimal performance. The fast convergence can be achieved by setting a large λ at the initial iteration loop, and the accurate estimation can be carried out through gradually decreasing the value of λ for the rest of the iterations. For example, we suggest a simple asymptotically decreasing model for λ :

Selection of the Threshold
It is worth noting that a careful choice of the threshold, i.e., λ, is recommended for the proposed FPA. A choice of small λ enables the features to contain a sufficiently large number of scatterers, which achieves accurate estimation. However, if λ is too small, the remaining features would also contain the artifacts of the SAR image such as side-lobes even though their magnitude are reduced. Thus, it may cause low speed of convergence, and insufficient estimation accuracy if the artifacts still remain in the features at the final iteration. On the other hand, if λ is too large, some dominant scatterers for phase estimation would be removed from the features. In this case, the estimation is too concentrated in a few remained strong scatterers, which leads to local optimal estimation even if it enables to achieve faster convergence rate.
To solve these tradeoffs, we propose varying the threshold to satisfy both fast convergence and optimal performance. The fast convergence can be achieved by setting a large λ at the initial iteration loop, and the accurate estimation can be carried out through gradually decreasing the value of λ for the rest of the iterations. For example, we suggest a simple asymptotically decreasing model for λ: where α and λ i represent the positive forgetting factor less than 1 and the threshold for i-th iteration, respectively. The threshold for the initial loop, i.e., λ 0 , would be a user-defined sufficiently large positive value, but small enough to achieve an appropriate number of features containing the main scatterers. There are various ways to select the main scatterers, such as the methods in [6,7,13]. Any of these methods would work to select the appropriate initial threshold. Alternatively, we observed that the proposed algorithm with an initial threshold slightly smaller than the maximum amplitude of the image, which enables the features to contain at least one scatterer, performs sufficient convergence and performance for every case in this paper. The performance comparison for constant and varying threshold is described in the next section.

Experimental Results
In this section, some experimental results are demonstrated to verify the benefits of our proposed algorithm. In Section 4.1, we verify the performance and characteristics of the proposed method with a different threshold λ, and explain the tradeoff between the convergence and accuracy. Then, we compare the results for the constant threshold with those of the varying threshold to compromise the tradeoff, which is introduced in the previous section. In Section 4.2, we demonstrate the convergence and the performance of the proposed method with various types of the phase error, and compare them with those of the existing autofocus algorithms. PFA is utilized for SAR imaging to all the experimental results. The quantitative measures to verify the performance of the autofocus algorithm are defined as where E(·) represents the spatial mean operator, and IC(ĝ ) and IE(ĝ ) represent the contrast and entropy of the imageĝ , respectively. The initial image g is scaled to have a maximum magnitude of 1 for all the cases of the experiment, which makes the threshold λ for the soft-threshold function within 0 to 1. The constant µ for stop criteria is set to be 10 -4 for all cases.

Performance and Convergence of the Proposed Method
The SAR image used for the proposed FPA is shown in Figure 2a. The size of the image is 4000 × 4000 pixels. The vertical and horizontal coordinates of the image represent the range and azimuth, respectively. The definitions of coordinates are the same for all of the other SAR images in this paper. Figure 2b,c shows the image enlarged at the point near the center of the image both without and with the phase error of Figure 3, respectively. The images corrected by FPA are shown in Figure 4, and the variation of contrast and entropy at each iteration are represented in Figure 5. The images in Figure 4a-d are the results for a constant threshold λ fixed by 0.01, 0.1, 0.3, and 0.9, respectively. Although the images look similar to each other, their contrast and entropy are slightly different, as shown in Figure 5 and Table 1. For a small value of λ, the convergence rate is low, even though the measures tend to converge to its optimal value. On the contrary, a large value of λ enables fast convergence, whereas the final performances are degraded. The reason for these characteristics can be explained by the number of features in Figure 5c. It shows that the smaller value of the threshold applied, the larger the number of features the algorithm uses for each iteration. As mentioned in the previous section, the features may contain the artifacts if the number of features is too large, which causes an adverse effect on the convergence rate. In the case of a small number of features, on the other hand, some dominant scatterers are omitted and it may interfere with the global optimality of the estimation, even if it achieves fast convergence.  As explained in the previous section, we applied a varying threshold to the proposed FPA to satisfy both fast convergence and optimal performance. We set the threshold λ for initial iteration to a relatively large value, and gradually decrease the value at each iteration as in (23). The initial threshold 0 λ and the forgetting factor α are set to 0.9 and 0.5, respectively, and the corrected image from this method is shown in Figure 4e. As shown in Figure 5a,b, the quantitative measures reach the optimal value with an appropriate iteration number. Unlike the fixed threshold cases, the number of features increases  As explained in the previous section, we applied a varying threshold to the proposed FPA to satisfy both fast convergence and optimal performance. We set the threshold λ for initial iteration to a relatively large value, and gradually decrease the value at each iteration as in (23). The initial threshold 0 λ and the forgetting factor α are set to 0.9 and 0.5, respectively, and the corrected image from this method is shown in Figure 4e. As shown in Figure 5a,b, the quantitative measures reach the optimal value with an appropriate iteration number. Unlike the fixed threshold cases, the number of features increases As explained in the previous section, we applied a varying threshold to the proposed FPA to satisfy both fast convergence and optimal performance. We set the threshold λ for initial iteration to a relatively large value, and gradually decrease the value at each iteration as in (23). The initial threshold λ 0 and the forgetting factor α are set to 0.9 and 0.5, respectively, and the corrected image from this method is shown in Figure 4e. As shown in Figure 5a,b, the quantitative measures reach the optimal value with an appropriate iteration number. Unlike the fixed threshold cases, the number of features increases for each iteration, as shown in Figure 5c, which enables the features of the current iteration to contain the dominant scatterers omitted at the previous iteration. for each iteration, as shown in Figure 5c, which enables the features of the current iteration to contain the dominant scatterers omitted at the previous iteration.  for each iteration, as shown in Figure 5c, which enables the features of the current iteration to contain the dominant scatterers omitted at the previous iteration.

Comparison with Existing Autofocus Algorithms
We have compared the proposed method with existing postprocessing autofocus methods of PGA [5], GPGA [7], and ME [12]. Constant false alarm rate (CFAR) detection was used to select the strongest scatterers for GPGA, which is a modified algorithm of PGA. The stop conditions for PGA, GPGA, and ME were the same as that of proposed FPA, and we used the varying threshold for FPA as described in Section 3.3. The value of λ 0 and α for FPA were the same as the previous experiment.
We added various types of phase errors as shown in Figure 6 to the scene in Figure 2a. The images corrected by PGA, GPGA, ME, and the proposed FPA are shown in Figure 7. The contrast and entropy variation for each iteration are represented in Figures 8 and 9, respectively. The performance measures of the images for each method and phase error are represented in Table 2, and the number of iteration and total computing time are shown in Table 3. The computation time was measured with a workstation equipped with Intel ® Xeon ® Gold 6140 CPU. MATLAB was used as the programming language.

Comparison with Existing Autofocus Algorithms
We have compared the proposed method with existing postprocessing autofocus methods of PGA [5], GPGA [7], and ME [12]. Constant false alarm rate (CFAR) detection was used to select the strongest scatterers for GPGA, which is a modified algorithm of PGA. The stop conditions for PGA, GPGA, and ME were the same as that of proposed FPA, and we used the varying threshold for FPA as described in Section 3.3. The value of 0 λ and α for FPA were the same as the previous experiment.
We added various types of phase errors as shown in Figure 6 to the scene in Figure  2a. The images corrected by PGA, GPGA, ME, and the proposed FPA are shown in Figure  7. The contrast and entropy variation for each iteration are represented in Figures 8 and  9, respectively. The performance measures of the images for each method and phase error are represented in Table 2, and the number of iteration and total computing time are shown in Table 3. The computation time was measured with a workstation equipped with Intel ® Xeon ® Gold 6140 CPU. MATLAB was used as the programming language.        It is well-known fact that PGA shows fast convergence and sufficient performance for low-frequency errors, such as the quadratic error shown in Figure 6a, which can be observed from the results in Figures 7a, 8a and 9a, and Tables 2 and 3. GPGA shows almost the same result with PGA with less iteration, but requires more computation time due to the additional process such as the strongest peak selection through CFAR for every iteration. The proposed method shows a similar trend with less computation time, because both PGA and GPGA require additional procedure for center-shifting, windowing, elimination of linear phase, etc. Meanwhile, ME requires more iterations due to the large phase error, even though it slowly converges to the optimal performance as shown in those figures. Unlike the quadratic error case, the image with random phase error cannot be corrected by PGA and GPGA, as shown in Figures 7b, 8b and 9b, and Tables 2 and 3. It is a natural result because of the assumptions and limited window size of PGA, described in Section 1. Meanwhile, ME and FPA show nearly optimal performance with few iterations. Therefore, it can be inferred from these results that the performance and convergence rate of the proposed FPA are not limited by the bandwidth of the phase error, unlike PGA and ME.
In a practical SAR system, Wiener process and discontinuous phase error can occurr because of the navigation systems for motion compensation. If we use an INS system, the Wiener process errors are generated due to the integration of IMU measurements that contain the Gaussian white noises. Furthermore, the navigation data experience some discontinuities if the system uses GPS measurement updates. The phase errors in Figure 6c,d are this kind of error, and the autofocus results for these errors are represented in (c) and (d) of Figures 7-9, and Tables 2 and 3. These results verify that FPA also shows the best performance with appropriate iteration number and computation time. Therefore, we verified the performance and convergence of the proposed FPA as well as its robustness for most types of phase errors through these results. We performed the same procedure for two more scenes in Figure 10 to verify the reliability of the algorithm. Scenes A and B are selected to have a higher and lower value of entropy, i.e., lower and higher contrast, than those of the image in Figure 2a. The results for these two scenes are represented in Tables 4-7, and similar trends of performance and convergence are observed when compared to that of the results in Tables 2 and 3. FPA shows best performance with sufficiently small iterations and computation time for all cases, and it verifies that FPA produces reliable performance.              (a)-(d) denote the types of phase errors, same as Figure 6.

Conclusions
In this paper, we proposed and demonstrated a new autofocus method for postprocessing of a phase-corrupted SAR image based on minimization of the cost function that consists of fidelity and regularization terms. The equation to achieve the optimality is derived by indirect optimization, and the algorithm to solve the equation is proposed. Each iteration in the proposed algorithm requires only one soft-thresholding for the reference image formation, and it enables more efficient processing than the existing regularizationbased autofocus such as SDA. The tradeoff between the performance and convergence for the proposed FPA can be compromised by selecting proper constant threshold or using an asymptotically decreasing threshold with appropriate initial value and forgetting factor. The experimental results verified its better performance, convergence, and robustness when compared to the existing methods of PGA and ME. We verified the reliability of the proposed method by performing additional two experiments with a different scene.
Although the proposed FPA shows sufficient performance and convergence for these experiments, there are still remaining factors to improve, such as selection of features and determining the threshold for each iteration. These factors would depend on the scene. Hence, the modifications through adaptive methods would improve the performance and convergence of the proposed algorithm, which are the future works for this study.