A Novel Ultra − High Resolution Imaging Algorithm Based on the Accurate High − Order 2 − D Spectrum for Space − Borne SAR

: Ultra − high spatial resolution, which can bring more detail to ground observation, is a constant pursuit of the modern space − borne synthetic aperture radar. However, the exact imaging in this case has always been a complex technical problem due to its complicated imaging geometry and signal structure. To achieve those applications’ strict requirements, a novel ultra − high resolution imaging algorithm based on an accurate high − order 2 − D spectrum is presented in this paper. The only ﬁrst two Doppler parameters needed as range models in the defective spectrum are replaced by a polynomial range model, which can derive coefﬁcients from the relative motion between the radar and the targets. Then, the new spectrum is calculated through the Lagrange inversion formula. Based on this, the novel imaging algorithm is elaborated in detail as follows: The range high − order term of the spectrum is compensated completely, and the range chirp rate space variance is eliminated by the cubic phase term. Two steps of range cell migration correct are applied in this algorithm before and after the range compression; one is the traditional linear chirp scaling method, and another is the interpolation to correct the quadratic range cell migration introduced by the range chirp rate equalization. The simulation results illustrate that the proposed algorithm can handle the exact imaging processing with a 0.25 m resolution around the azimuth and range in 2 km × 6 km, which validates the feasibility of the proposed algorithm.


Introduction
Space−borne synthetic aperture radar (SAR) is an ideal approach in earth observation, which has risen in recent years, both in military and civil applications, due to its advantages of capability with all−time and all−weather. To receive more detailed information in SAR images, the requirement of the resolution becomes higher. Currently, many advanced orbiting space−borne SARs such as TerraSAR−X [1][2][3], Cosmo−SkyMed [4,5], and GaoFen−3 are capable of reaching the 1 m resolution, and the TerraSAR NEXT GENERATION (Ter-raSAR NG) [6] plans to enhance this ability to 0.25 m.
In the traditional space−borne SAR system, many classic imaging algorithms possess great ability with the simple air−borne imaging geometry assumptions like the range Doppler algorithm (RDA) [7] and the chirp scaling algorithm (CSA) [8]. The hyperbolic range equation model (HREM) [9] or equivalent squint range model (ESRM) [10] applied in RDA and CSA only require the Doppler centroid and Doppler rate, which supposes the satellite trajectory as a straight line during the integration time. Meanwhile, these methods regard the three− and higher−terms of the signal spectrum as the phase error and ignore them. Conversely, the wavenumber domain algorithm (ωkA) [11] is a classic tion between the satellite and the targets is accomplished by complex multiplication in the RD domain. Furthermore, to get enough focus depth, the range chirp rate's space variance, which will render the range defocused, is non−negligible in this case. A non−linear FM filter method proposed by Davidson [25] is a common approach to compensate the residual phase raised by the squint angle in the quadratic range compression (SRC). It was validated through the simulation based on Seasat and ERS−1 and achieves the well−focusing with squint angle of up to 35 deg for L−band and 50 deg for C−band. Moreover, the article mentions a limitation: the reference Doppler frequency must be selected outside of the Doppler bandwidth, which will introduce the extra phase error. Another efficient method proposed by Wong [26] in 2008 is applied in this paper. The range chirp rate is assumed to vary linearly, so the cubic phase function can equalize the dominant quadratic phase term of each target, which is like the chirp scaling principle. Nevertheless, there is still a problem that needs to be handled. The residual linear phase introduced above would cause the chirp phase center to shift in range direction with the Doppler frequency and be reflected as the quadratic RCM after the range compression due to the corresponding square of range difference. This quadratic RCM is usually neglected because of the slight variation of chirp rate for low resolution. In contrast, it must be corrected for our research case, so the interpolation is applied to fulfill the correction and finish the range processing. Finally, the well−focused targets can be obtained after the azimuth processing.
This paper is organized as follows: The Infinite polynomial range model and the accurate high−order signal spectrum are introduced in Section 2. In Section 3, the proposed imaging algorithm for the ultra−high space−borne SAR is described in detail. In addition, an analysis of the proposed algorithm's performance and an integrated simulation is presented in Section 4. Finally, conclusions are drawn in Section 5.

Infinite Polynomial Range Model and Signal Spectrum
Eldhuset gives a four−order range model to describe the SAR range history in [16], but it only corresponds to the first four−order relative motion parameters between the satellite and the targets. To describe the exact range history of the long illumination time, a general infinite polynomial is presented in this paper. The range history vector → R(η) [16] can be expressed through the Taylor expansion in azimuth time η as while → X i is the ith−order relative motion parameters. Modulo both sides of Equation (1) and processed through the binomial theorem, the scalar value of the range history can be acquired in Equation (2).
Expanding Equation (2) by the Taylor series, the no error polynomial range model is given as while R i represents the ith−order range model coefficient.
The space−borne SAR received echo from a point target, after demodulated to the baseband, can be derived as where the σ 0 represents the target's back scattering coefficient; ω r (·) and ω a (·) denote the antenna pattern functions on the range and azimuth directions, respectively; c is the velocity of light; η c is the Doppler central time; f 0 is the signal carrier frequency; K r is the range chirp rate; τ is the range time.
To acquire the 2−D spectrum of the echo, the POSP is applied. It is easy to achieve this operation on the range direction and obtain the signal in the range−frequency (RD) domain as where the f τ is the range frequency. However, in the azimuth direction, the stationary phase point cannot work out through differentiating η on the integrand phase Equation (7) directly, so we must find other ways to solve this problem.
According to the time−frequency properties of the Fourier transform, the linear term multiplied in the time domain is equivalent to a frequency shift. Therefore, we can extract the linear term in Equation (3) and acquire the new form of the signal as while The new integrand phase of Equation (9) can be expressed as Equation (10), and the stationary phase point can be derived from its differential.
To make the followed steps more intuitive, variable P is used to simplify the differential equation as Obviously, Equation (12) is an infinite series without any constant term. Therefore, we can use the Lagrange inversion [27] to obtain its inverse function as another infinite series by while So that the 2−D spectrum of S ( f τ , η) can be achieved by introducing Equation (14) to Equation (10), and the result is given as Finally, the linear term R 1 η is added as a frequency bias into the Equation (16) and the precise high−order 2−D spectrum of s(τ, η) is elaborated as

Imaging Algorithm
Based on the high−order 2−D spectrum proposed above, a novel ultra−high−resolution imaging algorithm for space−borne SAR is presented in this section. The sever azimuth aliasing caused by the spotlight or sliding spotlight mode is eliminated by the two−step processing method [28,29]. Then, a range high−order phase filter and a linear RCMC filter are applied to fulfill the range preprocessing. Further, a cubic phase function is introduced to eliminate the space−variance of the range chirp rate, and a quadratic RCMC filter (after the range SRC) is employed to correct the new RCM generated by the chirp rate equalization. Finally, the azimuth compression filter is applied to focus the targets, and the azimuth scaling is adopted to avoid the azimuth time domain aliasing caused by the deramp operation. Figure 1 gives the block diagram of the proposed imaging algorithm, and all these parts mentioned above are illustrated in the following:

Azimuth Preprocessing
The total Doppler bandwidth of an ultra−high−resolution space−borne SAR echo is always far larger than the pulse repetition frequency (PRF), which is limited by the range ambiguity and data volume. Therefore, the signal spectrum suffered from the severe aliasing at the azimuth direction. The time−frequency diagram of the sliding spotlight, which is widely developed in many advanced SAR systems and applied in our simulations, is elaborated in Figure 2. The azimuth bandwidth of the whole scene B total is composed of B steer and B 3dB , while B steer is rendered by the azimuth antenna steering and B 3dB is the antenna beam 3 dB bandwidth. B T is the single point target's azimuth bandwidth and T B is the azimuth integration time. As shown in Figure 2, the aliasing is inevitable and must be removed for post−processing.
The two−step processing method is an efficient method and is employed in our algorithm. It utilizes a reference LFM signal with an opposite Doppler chirp rate convoluted with the received signal in the azimuth direction. B steer is removed through this operation and the processed signal can be expressed as Equation (20) s (τ, η) = s(τ, η) * exp −jπK r_rot η 2 (20) while K r_rot is the azimuth chirp rate introduced by the antenna steering and * is the convolution operator. The sketches of the spectrum, whether processed by this method or not, are elaborated in Figure 3, which demonstrate the azimuth wrapping is eliminated clearly in Figure 3b.  In addition, the convolution operation requires many computing resources. To simplify the calculation, it can be divided into three parts [30]: the two complex multiplications and one Fast Fourier Transform (FFT) as Equation (21) while η 1 = K r_rot η, and the PRF is increased equivalently as while N a is the azimuth sampling number. The unwrapped 2−D dimension signal spectrum can be obtained based on the time−frequency properties of the Fourier transform as Equations (23) and (24).

High−Order Phase Compensation
The phase mathematical expression of the 2−D spectrum is portrayed in Equation (24), which embodies significant range−azimuth coupling. To illustrate validly, the phase of spectrum can be decomposed by the Taylor series with f τ and expressed as Equation (25).
while D 0 f η is the term of azimuth Doppler frequency modulation, which is recognized as the LFM term in most instances. D 1 f η includes the information of the RCM for each target, and D 2 f η corresponds to the range modulation. Their approximate expression is given in Appendix A, and the remaining higher order terms donate the high−order range-azimuth coupling.
Many frequency algorithms are based on the first−three order of the spectrum and focus less on the others, such as the RDA and CSA. However, the negligence of the high−order term does not work for the high−resolution case. The range distortion caused by it is elaborated in Figure 4. [24] gives an analysis of phase error introduced by each order term of the resolution of 1 m, and concludes that the sixth order coupling phase term still requires considering. Through this, it is evident that the higher resolution required, the higher phase terms need to be considered. To obtain the well−focused targets, the high−order phase term is compensated in this paper, which can be calculated by (26) and the compensated 2−D spectrum can be achieved by Transforming the signal to the range time domain by a range inverse fast Fourier transform (IFFT), at last, the compensated signal in the RD domain is achieved by

Linear RCMC
As the principle of chirp scaling, the chirp signal can adjust the phase center by multiplying itself with another chirp signal. By these means, the different RCM produced by each target, which is located in different range gates, can be adjusted to the same, and this is named differential RCMC in many papers.
For the expression, according to the Equations (28) and (29), the range chirp rate br f η and range cell migration(RCM) R rd f η in the RD domain can be derived as Therefore, the chirp scaling factor is where R re f represented the referenced slant range of the referenced range gate. The differential RCMC is achieved by multiplying Φ cs τ, f η with S c τ, f η in the RD domain.
After the chirp scaling operation is finished, the remaining RCMC can be achieved together by the RCM of the referenced range gate and is always realized by a complex multiplication in the 2−D frequency domain. According to the time−frequency properties of the Fourier transform, the remaining RCMC filter in the 2−D frequency domain can be expressed as Equation (36) and the signal in RD domain after finished the linear RCMC is while the R is the referenced slant range for each range gate.

Chirp Rate Equalization and Quadratic RCMC
The chirp rate is space−variant in the RD domain. This variance is always regarded as a system error for the traditional imaging algorithm and neglected because of the low resolution. However, for the ultra−high−resolution case, this is sufficient to lead the targets at the edge to defocus seriously. Based on the first exponential term in Equation (37), the range chirp rate in the RD domain can be expressed as K r f η , R while it is assumed to vary linearly with R in this paper like Equation (39) and A f η is the rate of the K r f η , R . So, this quadratic exponential term can be re−expressed as Based on [26], it is feasible to apply a cubic phase perturbation factor, given in Equation (41), to equalize the chirp rate, like the chirp scaling principle the quadratic term after multiplied with Equation (41) becomes Equation (42).
The chirp rate for each target located at each range gate is equalized to the reference chirp rate in Equation (42), and it can be well−compressed with a match filter Equation (43) in the RD domain accurately. The first exponential term is the cubic phase term, which can be compensated in the RD domain by the POSP, and the last exponential term is constant.
The most difficult is the remaining linear term. It will cause the phase center shift in the range time domain and bandwidth shift in the frequency domain. Moreover, a false target will be produced if the bandwidth shifts too much. To explain this phenomenon straightforwardly, Figure 5 shows the false target caused by the spectrum wrapping. The bandwidth of the original spectrum is included in [−F s /2 : F s /2], so its compressed result is well−focused and located at the phase center of the chirp signal. However, if the bandwidth shifts out of the sampling frequency, as illustrated in Figure 5a, the energy leakage will occur, and the false target will emerge, as shown in Figure 5b. Therefore, the sampling frequency must be large enough to ensure the bandwidth does not alias after the shift. Nevertheless, phase center shifting is inevitable in the range time domain. Figure 6a gives a plot of a real part from a chirp signal, adjusted to the chirp rate by the aforementioned cubic phase perturbation factor. The phase center of the original signal is located at the center of the plot. However, it shifted to the left after applying the chirp rate equal-ization, displayed in Figure 6a below. Their compressed results, focused on the reference chirp rate, are displayed in Figure 6b. The defocusing emerges on the original signal due to the mismatch of the match filter, and the equalized signal reflects the well−focusing. The peaks of each pulse appear in the phase center so that there is a deviation of location between them. When this bias varies along with the azimuth frequency, it manifests as the RCM introduced artificially. Furthermore, due to the shifts corresponding to the square of the range difference, which embodies the third exponential of Equation (40), the newly generated RCM is represented as quadratic and toward the same direction as illustrated in Figure 7.  The principle of chirp scaling is no longer applicable for the correction because the quadratic RCM is introduced by itself. To accomplish the RCMC, interpolation is employed in this paper. The shift of each range bin can be calculated in Equation (44) precisely without any high−order term. Therefore, the interpolation kernel is longer, and the RCMC by this method is more accurate. At last, the expression of the signal after RCMC and residual phase compensation is given in Equation (45).
while the A r (·) is the range envelope after range compression and usually as a sinc function.

Azimuth Compression and Resampling
According to the Equation (45), the remaining phase includes two parts. The penultimate term is the azimuth modulation phase mentioned above. The last exponential term is introduced by the chirp scaling operation, both of which need to be compensated accurately. Therefore, the azimuth compression filter can be obtained through the conjugate of these parts.
To eliminate the azimuth time domain aliasing caused by the deramp, the azimuth scaling option is necessary to achieve the azimuth resampling. The scaling function is given in Equation (47) as while the H f is the hybrid factor. After performing the azimuth IFFT, the residual phase is given in Equation (48), and the focused target can be expressed as Equation (49).
while the A a (·) is the azimuth envelope function.

Discussion and Simulation Results
The discussion and experiment results of simulation are presented in this section. To verify the proposed algorithm, the main simulation parameters are given in Table 1.

The Discussion of The Polynomial Range Model and High−Order 2−D Spectrum Analyses
The infinite range model is unachievable in engineering. Therefore, the truncation of it is necessary. There are two infinite series while the range model deduction. First is the range vector illustrated in Equation (1). Each order relative motion parameter vector is required in this expression. However, the higher−order term is non−essential because the illumination time is limited. Figure 8 shows the phase error calculated by the modulus of the low−order truncated range vector and the real range history, and the red line annotates the safe line under the azimuth phase error criterion of 0.25 π. As portrayed in Figure 8, the first four−order relative motion parameters can fit the real range history in 20 s. This is enough to meet the requirement of the ultra−high−resolution imaging for now, and for some time to come. The second infinite series is the scalar expression of the range model acquired by the Taylor expansion. Similarly, it is not necessary and impossible to apply all coefficients to obtain the no error range history. The first several order approximations will offer acceptable performance and reduce the calculation in engineering. Figure 9 portrays the approximation results based on the first−four order relative motion parameters in the same manner as above. Moreover, the only even−order truncated results are elaborated in the plot because the odd−order coefficients are calculated close to zero through the simulation parameters presented in Table 1. The six−order approximated range model is enough for the 20 s integration time, while the requirement of the ultra−high−resolution imaging is already achieved. However, the curve has an apparent upward trend after −6 s, so to satisfy the development of the SAR system in the future, the eight−order range model will be a better choice, and the remaining experiments are based upon it. Furthermore, we aimed to verify the accuracy of the high−order 2−D spectrum generated by the Lagrange inversion. The peak sidelobe ratio (PSLR) and the normalized impulse response width (IRW) measured from the point targets, compressed by the calculated spectrum directly based on the different order truncated Lagrange inversion, are presented in Figure 10. Obviously, the three−order Lagrange inversion is far from reaching the requirement of ultra−high−resolution imaging. Serious defocusing has emerged since the 0.8 m of the azimuth resolution. In addition, the four− and five−order inversions exhibit the near−precision of the derived spectrum because the fifth−order coefficient is too small to affect the precision, which is almost the same as the sixth−order parameter in this experiment. Moreover, they are all invalid since the azimuth resolution is higher than 0.17 m, where the spectrum obtained by the six−order inversion still performs well. Similarly, the accuracy of the 2−D spectrum can be improved by increasing the order of the Lagrange inversion to approach the actual accuracy like the Taylor series, and the six−order inversion is employed in the next simulation experiment to calculate the 2−D spectrum.

Image Quality Evaluation and Analysis
The simulated scene is shown in Figure 11, the nine point targets arranged evenly along with the range, and azimuth in 6 km and 2 km as illuminated in the sliding spotlight mode. The main simulation parameters are given in Table 1. To fit the engineering applications, all simulation results are weighted by the 1/3 simplified Taylor window both in range and azimuth, which can reduce the sidelobe level and avoid flooding the weak targets [31].  Figure 12 illustrates the necessity of the range chirp rate equalization. The diagrams of the PT1, whether it eliminated the chirp rate space variance or not after the range compression in the RD domain, is given in Figures 12a and 12b respectively. As seen in Figure 12a, the mismatch of the match filter is reflected obviously, especially at the edge of the Doppler bandwidth. This leads the results to range defocusing finally. In contrast, Figure 12b shows the influence of the chirp rate equalization achieved by the cubic phase function. The mismatch is eliminated effectively through this method. Conversely, a tiny shift in the range direction corresponding to the Doppler frequency is caused and brings the azimuth distortion for each target in the end. To fulfill the quadratic RCMC, the interpolation is employed, and its result is presented in Figure 12c, in which it can be observed that every shift at each Doppler frequency is removed clearly, and the whole energy of a target is well−compressed. The interpolated slices of PT1, PT5, and PT9 focused by the proposed algorithm are displayed in Figure 13. In contrast, Figure 14 gives the same targets focused by the conventional CSA, and Figure 15 gives the results from the BPA. The range and azimuth profiles of the PT5 are presented in Figure 15.   As seen from the results processed by CSA, due to the tremendous phase error introduced by the ESRM in the ultra−high−resolution case, each point target suffers from severe degradation, especially at the edge of the azimuth extent. The insufficient precision of the 2−D spectrum computed in CSA causes the rise of PSRL and the widening of IRW. In Figure 16, the profile cannot be recognized as the sinc function. Conversely, the results from the proposed algorithm, as presented in Figure 13, indicate that with the more accurate range model and the 2−D spectrum based on it, the focusing performance is dramatically improved. Every point target on the scene is well−focused. The results from the BPA are displayed in Figure 15 as a reference. The well−focused targets validate the high imaging accuracy of the time domain algorithm. Compared with Figure 13, two imaging algorithms have similar performances in image quality, but the proposed algorithm only requires eight times FFT, seven times complex multiplication, and one interpolation, which is far more efficient than the BPA.
To better quantify the image quality, the point target analysis results are listed in Table 2. The theoretical resolution, weighted by the 1/3 simplified Taylor window, is calculated in where L a is the antenna length and H f is the hybrid factor. It can be observed from Table 2 that the deterioration of the IRW is less than 1%, both in range and azimuth direction. Meanwhile, the PSRL degradation in the table is less than 2%. These promising results prove the proposed algorithm has remarkable performance and can meet the requirement of the ultra−high−resolution imaging for the space−borne SAR for the present and future.

Geolocation Analysis
The information of geolocation, which is widely used in the post−processing of SAR images, is another criterion for evaluating the performance of the imaging algorithm. Table 3 gives the deviation of the positioning results for every focused target in the image scene, while the targets are located in the WGS84 coordinate system. These results indicate that the deviation at each coordinate axis is less than 10 −2 m, proving the proposed algorithm's high position accuracy.

The Discussion of the LEO SAR Motion Error and Autofocus Algorithms
As can be observed from the infinite polynomial range model, the exact relative motion parameters are required for every order. Nevertheless, the parameters error, especially for the higher order, is inevitable in engineering and degrades the image quality. Motion compensation is a common method of handling this problem. It relies on the inertial navigation system installed on the SAR load to record the movement of the radar and compensate for the residual phase through the stored data. This method requires very high accuracy of the inertial navigation system, which increases the system's final cost.
Autofocus is another way to solve the problem. This algorithm is driven by the echo or image data without any auxiliary equipment. The phase gradient autofocus (PGA) [32] algorithm is the most widely used autofocus method, and the maximum likelihood estimation is employed in the kernel to calculate the phase gradient based on the strongest scattering point in the defocused SAR image, and the estimated phase error gradient compensates the phase error. It is worth noting that PGA can eliminate any order of phase error so that it has a good performance in the high−resolution space−borne SAR. With the development of deep learning, many autofocus algorithms based on the neural network are proposed now [33,34]. They treat the entire imaging process as a linear transformation of the SAR signal and use the neural network to accomplish the linear mapping. Under certain conditions, they have a similar or even better performance than the classic autofocus algorithm, but they still have some limitations, such as only being applied on the polynomial form phase error [34] and the vast number of samples required while the neural network is training.

Conclusions
In this paper, a novel 2−D spectrum required imaging algorithm for the ultra−high− resolution space−borne SAR is proposed. To achieve the accuracy of ultra−high−resolution, a high−precision high−order 2−D spectrum is illustrated in detail at the beginning. The infinite polynomial range model is introduced to the 2−D spectrum analysis, and the expression of the new spectrum is derived from it through the Taylor expansion and Lagrange inversion. In addition, the novel imaging algorithm kernel is presented. The azimuth spectrum aliasing caused by the antenna rotation of the sliding spotlight mode is eliminated by the two−step processing. Moreover, to keep the symmetry on the range direction of the point target's impulse response, the high−order phase term of the echo's spectrum is compensated through the derived new spectrum expression. The linear RCM is corrected through the principle of chirp scaling. The space variance of the chirp rate in the range direction is non−negligible in the ultra−high−resolution case. To solve this problem, a cubic phase function is introduced to equalize the chirp rate and guarantee the quality of range focusing. Meanwhile, the quadratic RCM introduced by the cubic phase function is corrected by interpolating after the range compression. In the end, the azimuth compression is implemented, and the well−focused targets of the large scene can be obtained.
The accuracy of the new spectrum is discussed with the order of the Taylor and Lagrange series in the experiment in this paper, and it verifies that the new spectrum can achieve the requirement of ultra−high−resolution. Then, the image quality of the algorithm is evaluated by the simulation. As a comparison, the imaging results from the traditional CSA are also presented. From the detailed view, it is obvious that every well−focused target derived from the proposed algorithm is defocused and distorted by the CSA. The table of the PSLR and ISLR, calculated from the well−focused targets, further proved the performance of the proposed algorithm. Additionally, the geolocation analysis is also given in this paper. The high positioning accuracy shown in the table ensures the post−processing of a SAR image. All these results verify the effectiveness of the proposed imaging algorithm.