High-precision, fast geolocation method for spaceborne synthetic aperture radar

Geolocation using spaceborne synthetic aperture radar (SAR) is essential for imagery applications, and high performance geolocation methods need to be developed to promote SAR imagery applications. Starting from the SAR imaging principle, this paper reveals and analyzes two basic characteristics of SAR imaging geometry, and demonstrates the rationality of the two characteristics. On this basis, a high-precision and fast geolocation method is proposed. We conducted a precision analysis on four SAR satellites (Germany’s TerraSAR-X, Italy’s COSMO-SkyMed, Japan’s ALOS-PalSAR and Canada’s Radarsat-2 satellites), and the results show that the precision of the proposed method meets practical needs. We then used TerraSAR-X SpotLight SAR real data to implement the fast geolocation, and found from performance evaluation that the computation cost is greatly reduced while high geolocation accuracy is maintained. We thus verified the efficiency and accuracy of the proposed method.

With the technological advancement of spaceborne synthetic aperture radar (SAR), SAR image data are being increasingly used [1][2][3] in military and civil applications. Because of the effects of the side-looking operation of SAR and topography on the SAR imaging geometry, an SAR image is more geometrically distorted than an optical image. Therefore, geometric correction or geocoding is needed for high-precision SAR geolocation.
Currently, the most widely used geolocation model applied by spaceborne SAR is the Range Doppler (RD) model proposed by Brown [4] and Curlander [5]. The model has been widely investigated in terms of model parameter extraction and optimization, coordinate transformation, model solving, and geolocation error source analysis. The three equations applied by the RD model are the range, Doppler and Earth model equations. Governed by the three equations, the exact solution for target geolocation can only be obtained using a complicated iteration method. Li et al. [6] was the first to assume that the Earth is a sphere with a local radius at the subsatellite point. Under this assumption, the location of the target on the Earth's surface can be determined analytically and without iteration. This represents a significant saving in computation. Yuan [3,7] conducted a systematic study on the geolocation analytic solution for this approximation. Several algorithms [8][9][10][11][12] have been proposed for the solutions of the RD model with a variety of approximations. The US Jet Propulsion Laboratory, German Aerospace Center, Italian Space Agency and other international research institutions have developed corresponding SAR image geolocation processing modules based on the RD model.
Initially, because of poor measurement accuracy for the satellite orbit, the loss of geolocation accuracy resulting from making a variety of approximations was acceptable for near Earth remote sensing missions. However, with the development of measurement technology for spaceborne radar, especially with the improvement in accuracy of the orbit measurement, SAR geolocation has achieved meter level accuracy or better. The geolocation errors resulting from employing conventional approximation methods are in the same order like the geolocation accuracy. Taking the German TerraSAR-X satellite for example, the measurement accuracy of its science orbit is 3 cm [13], and the geolocation of the SAR imagery product has accuracy better than 10 m or even submeter accuracy [14,15] in certain working modes. Here, the loss of geolocation accuracy due to employing fast geolocation methods based on approximations cannot be ignored. In addition, the use of arbitrary geolocation algorithms affects the geolocation accuracy and more importantly, may affect the results of scientific investigations. There is a recent important example in astronomy where the algorithm used for the Wilkinson Microwave Anisotropy Probe [16] produces an illusion [17] of a larger quadrupole for cosmic background radiation. Thus, it is important to develop new methods for highly accurate and fast geolocation, either by improving the geolocation accuracy of SAR imaging products or investigating other issues.
Considering the basic characteristics of the RD model, this paper presents a high-precision and fast geolocation method for spaceborne SAR. The method solves the problem of the computational cost of an iterative algorithm and greatly improves the operation speed, with a loss of geolocation accuracy of only 10 4 m, which is insignificant when comparing with the exact iterative solution. The effectiveness of the method is demonstrated in processing Ter-raSAR-X SpotLight data.

Analysis and characterization of the RD model
Spaceborne SAR geolocation involves solving the threedimensional coordinates of a ground target point for the image pixel by combining SAR image information and radar satellite orbit data obtained by imaging processing with a range equation, Doppler equation and Earth model equation. The range equation is the Doppler equation is usually, an oblate ellipsoid is used to model the Earth's shape where R is the slant range of the SAR imagery; P S and V S are the position and velocity vectors of the SAR antenna phase center respectively; P T = (x T , y T , z T ) and V T are position and velocity vectors of the ground target point, respectively; V T is equal to the zero vector in the geographic frame coordinate system and  is the wavelength; f DC is the Doppler center frequency of SAR imaging; R e and R p are the equatorial radius and polar radius of the Earth, respectively; h is the corresponding height of the assumed Earth model for the target. Note that the validity of the Earth model affects the SAR geolocation accuracy, and some scholars [18] have thus modified the Earth model equation. The Earth model we apply here is widely used, and the proposed method is applicable to other Earth models.
Eqs. (1)-(3) form a nonlinear system, and an iterative approach is thus required for target geolocation. An important task for a spaceborne SAR system is to make global observations for the creation of global SAR products. However, it is difficult for a conventional point-by-point iterative algorithm to meet needs in the processing of a massive quantity of data. Scholars are thus devoting great efforts to develop fast solution methods for SAR geolocation.
Geometric distortions in SAR imagery are mainly due to side-looking radar imaging and variation in the terrain elevation. We found from analysis that the space geometry model ( Figure 1) of spaceborne SAR imaging has two basic characteristics. (i) The geometry change in SAR imaging is locally small and successive; that is, the imaging geometry for point M in Figure 1 is little different from that of the four surrounding points (M i , i = 1,2,3,4). In spaceborne SAR image processing, the parameters of a two-dimensional matched filter are determined by the SAR imaging geometry. Additionally, the spatial variance of the filter can be ignored over a range of tens to hundreds of kilometers along the azimuth considering the resolution of typical spaceborne SAR at present. Usually, the depth of focus of tens to hundreds of meters along the range can also be ignored. Even the two-dimensional spatial variance curve of the parameters of the matched filter is flat, and can be considered a loworder curve. These characteristics as the basic premise for the frequency-domain fast imaging method have been fully verified. The spatial variance of a two-dimensional matched filter is a direct reflection of the spatial variance of the imaging geometry, which explains the gradual variance of SAR space geometry. (ii) Restricted by eqs. (1)-(3), the relationship between elevation h of point M and the three-dimensional coordinates (x, y, z) can be described by polynomials p x (h), p y (h) and p z (h).
We assume that second-order polynomials p x (h), p y (h) and p z (h) are used to approximate the relationship between three-dimensional coordinates (x, y, z) of a target point and the corresponding elevation h. According to the Taylor expansion theorem, when a function where  (h a , h b ) and  (h) = (hh 0 )(hh 1 )(hh 2 ). The loss of accuracy due to the approximation of the polynomials in SAR geometry imaging can be obtained by taking the approach of Wang et al. [19]. For the sake of brevity, only the accuracy results are presented here. When the system parameters of the four most used SAR satellites (Germany's TerraSAR-X, Italy's COSMO-SkyMed, Japan's ALOS-PalSAR and Canada's Radarsat-2) are substituted into the error model, the loss of accuracy arising from using the second-order polynomials is as shown in Figure 2.
In Figure 2, the horizontal ordinates correspond to the actual Earth elevation range (0-9000 m) and the vertical ordinates correspond to the working range of the viewing angle for the SAR satellite system. Figure 2 shows that for the four international SAR satellites, the accuracy loss in geolocation due to the use of a second-order polynomial is of the order of 10 4 m, which is acceptable in terms of actual needs in geolocation. The above analysis verifies the two basic characteristics of the SAR imaging geometry and that the use of a second-order polynomial results in sufficient geolocation accuracy.

Realization of the proposed method
Section 1.1 introduced the basic principles of the RD geolocation model and analyzed and verified the two basic characteristics of SAR imaging geometry. On this basis, a highly precise and fast geolocation method is presented.
Here the steps taken in implementing the method are described.
Step 1: Constructing the derivative Earth model. We assume that n Earth model equations were derived when substituting n elevations (h i , i=1,2,···n) into the Earth model.   (5) n is greater than the order of the polynomial to be fitted. h i is well distributed in the elevation range of the observation area of the SAR image. In the case of an unknown observation area, the values of h i are well distributed within the actual Earth elevation variation range (0-9000 m). Taking the second-order polynomial for example, we have n=3, h 1 =0 m, h 2 =4500 m, and h 3 =9000 m.
Step 2: Coarse sampling and geolocation of an SAR image. A coarsely sampled SAR image is obtained after coarse sampling of the original SAR image at certain intervals (k, l respectively) along the azimuth and range. The geolocation of the coarsely sampled SAR image is determined with the corresponding range equation, Doppler equation and all the derivative Earth model equations. Additionally, an iterative algorithm is employed for the geolocation solution. We thus obtain n groups of geolocation results P Ti =(x Ti , y Ti , z Ti ), i=1,2,···n governed by n Earth model equations for the coarsely sampled SAR image.
Step 3: Polynomial fitting and interpolation. From the n groups of geolocation results P Ti = (x Ti , y Ti , z Ti ), i=1,2,···n obtained in the second step with the corresponding Earth model elevation (h i , i=1,2,···n), three m(m< n)-order polynomials p x T (a s , r s ) (h), p y T (a s , r s ) (h) and p z T (a s , r s ) (h) with elevation h as the independent variable are fitted for each pixel position (a s , r s ) of the SAR image after coarse sampling. The three polynomials independently describe the relationship between the geolocation solution results x T (a s , r s ), y T (a s , r s ) and z T (a s , r s ) for pixel (a s , r s ) and elevation h. The polynomials p′ x T (a, r) (h), p′ y T (a, r) (h) and p′ z T (a, r) (h) for each pixel position (a, r) of the original SAR image are then obtained by bilinear interpolation of the polynomials for the pixel position in the SAR image after coarse sampling where where bilinear() denotes bilinear interpolation and floor() returns the closest integer less than or equal to the argument.
Step 4: Geolocation solution for the original SAR image. The actual Earth model elevation h e is substituted into polynomials p′ x T (a, r) (h), p′ y T (a, r) (h) and p′ z T (a, r) (h) to determine the geolocation result P T (a, r) for each pixel of the original SAR image.
The above geolocation solution only requires an iterative algorithm for one pixel of the SAR image after coarse sampling, and for other pixels, the corresponding geolocation result can be obtained immediately by directly evaluating the polynomials. In this process, the number of repeated iterations is greatly reduced by polynomial approximation, and the calculation cost decreases. Here we prove the efficiency of the method in terms of computation complexity, which is an important indicator of the efficiency of algorithms. The most accurate representation of complexity is the numbers of addition and multiplication operations carried out by the algorithm. For a complicated hybrid algorithm, however, it is difficult to accurately count the numbers of addition and multiplication operations, and asymptotic complexity is thus used to describe the complexity of the algorithm. Because the proposed method applies a hybrid algorithm, iteration is employed in solving nonlinear equations, polynomial fitting and interpolation and other mathematical processing. Therefore, we used asymptotic complexity as the indicator of complexity. When using iteration to solve nonlinear equations, we take the solution for a current pixel as the initial value for the next pixel, and the average number of iterations required for each pixel is then about 5. According to [20], the computational complexity of conventional point-by-point iteration is 5MNO(3 2 ), where M, N are the numbers of rows and columns in an SAR image. The computational complexity of the proposed method is where n is the number of elevations, m the order of polynomials to be fitted, and k, l the coarse sampling intervals along azimuth and range. Taking the real data in section 2 as an example, M = 6207, N = 10912, m = 2 and n = 3, and the computational complexity of the proposed method becomes less than that of conventional point-by-point iteration when the product kl > 3.4 for coarse sampling intervals along the azimuth and range. Thus, the high efficiency of the proposed method is demonstrated in terms of computational complexity.

Experimental results
The real data used in this section were obtained from Germany's TerraSAR-X satellite operating in SpotLight mode on February 23, 2009. The image had dimensions of 6207 × 10912 pixels, was centered on 25.34° latitude and 131.02° longitude, and had azimuth resolution of 1.1 m and range resolution of 0.6 m. The observation area was of Ayers Rock, a famous landmark in Australia; the circumference of the rock base is about 8.5 km and the height of the rock is 348 m, making Ayers Rock the world's largest rock. The proposed method is employed to determine the geolocation of the real data, and the results are presented in Figure 3. Figure 3 shows that the SAR image after geolocation processing fits well with the texture information of roads and other features in a Google Earth image. Figure 4 shows the variation curves of each order coefficient for secondorder polynomials that correspond to data lines running along the range (red dotted line) and azimuth (blue dotted line) in Figure 3(a). Figure 4 shows that each polynomial can be approximated as a linear variation along the azimuth and range, which once again demonstrates the first characteristic of SAR imaging geometry. Geolocation errors of SAR products can be divided into two types: geolocation errors arising from various measurement sources (e.g., orbit measurement error and atmospheric delay error) and geolocation errors arising from approximation when solving the geolocation model. This paper presents a new geolocation solution, and taking the same measurement sources of error, we need to evaluate the geolocation errors resulting from employing the proposed method. By referencing to the geolocation solution results of the point-by-point iteration, we evaluated the performance of the proposed method. We used a 2.33 GHz Intel Core 2 Quad CPU with 2 GB RAM and Interactive Data Language in the experiment. Table 1 presents the processing time for the point-bypoint iteration solution, and the processing performance with different parameters for the proposed method. The root-mean-square error is taken as a measure of coordinate errors. Thus, as the sampling interval increases, the computation cost significantly reduces while high accuracy is maintained when applying the proposed method. However, when the coarse sampling interval is greater than 30 × 30 pixels, the running time appears to be stable on the whole. This is because the geolocation time consists of two parts: the time to reach the iteration solution and the time for polynomial interpolation. The SpotLight SAR image has relatively few pixels, and with the coarse sampling interval increasing, the time required for polynomial interpolation as a

Figure 4
Curves of each order coefficient for the second-order polynomial. The first column are coefficients of each order for the X coordinate polynomial, the second column are coefficients of each order for the Y coordinate polynomial, and the third column are coefficients of each order for the Z coordinate polynomial. growing share of the total time required gradually becomes constant. When processing real data recorded by Ter-raSAR-X in StripMap mode (32000 × 15000 pixels) for containing more pixels, the processing speed is about 200 times faster for a sampling interval of 100 × 100 pixels, and the accuracy loss of geolocation is of the order of 10 4 m. Therefore, in applications, we choose appropriate solution parameters that strike the best balance between efficiency and accuracy. Note that the proposed method with a slight change can be applied to SAR image geocoding, including ellipsoid-corrected and terrain-corrected geocoding. In summary, the proposed method has an obvious performance advantage of combining efficiency with accuracy, and will be highly useful when processing massive amounts of data recorded in SAR global observations.

Conclusion
The high-precision geolocation of spaceborne SAR data is a prerequisite for the application of such data. In handling massive amounts of data recorded in global observations, it is important to consider how to improve the data processing speed while maintaining geolocation accuracy. Existing efficient geolocation algorithms have low precision and hardly meet the growing demand for geolocation accuracy.
On the basis of analysis and verification of the two fundamental characteristics of SAR imaging geometry, a highprecision and fast SAR geolocation solution was presented. The proposed method was used to process SpotLight SAR data obtained from the TerraSAR-X satellite on February 23, 2009, and the performance evaluation demonstrated the effectiveness of this method, which had greatly reduced computation costs while maintaining high accuracy. The implementation of the proposed method paves the way to realizing the fast mass production of global high-precision SAR imagery products.
The authors thank Germany Aerospace Center (DLR) and Infoterra for providing TerraSAR-X data, Michael Eineder (DLR) for many communications that formed and developed the core idea for this article, and reviewers for their valuable comments.