A Fast, Three-Dimensional, Indirect Geolocation Method Using IAGM and DSM Data without GCPs for Spaceborne SAR Images

To determine the geolocation of a pixel for spaceborne synthetic aperture radar (SAR) images, traditional indirect geolocation methods can cause great computational complexity. In this paper, a fast, three-dimensional, indirect geolocation method without ground control points (GCPs) is presented. First, the Range-Doppler (RD) geolocation model with all the equations in the Earth-centered rotating (ECR) coordinate system is introduced. By using an iterative analytical geolocation method (IAGM), the corner point locations of a quadrangle SAR image on the Earth’s surface are obtained. Then, a three-dimensional (3D) grid can be built by utilizing the digital surface model (DSM) data in this quadrangle. Through the proportional relationship for every pixel in the 3D grid, the azimuth time can be estimated, which is the key to decreasing the calculation time of the Doppler centroid. The results show that the proposed method is about 12 times faster than the traditional method, and that it maintains geolocation accuracy. After acquiring the precise azimuth time, it is easy to obtain the range location. Therefore, the spaceborne SAR image can be geolocated to the Earth surface precisely based on the high-resolution DSM data.


Introduction
Spaceborne synthetic aperture radar (SAR) images are widely used in a variety of fields [1]. With the paramount need to target surveillance and state security, i.e., target geolocation, the accuracy of the absolute locations of sensitive targets on the ground in spaceborne SAR images is very critical. The absolute accuracy of SAR geolocation depends on multiple factors, such as the orbital precision of the satellite, sensor stability, radar accuracy, knowledge of the propagation medium, SAR processing accuracy, timing accuracy, target detection accuracy, coordinate transformation accuracy, etc. [2][3][4][5][6][7].
The traditional indirect geolocation method establishes the mapping association between image coordinates and geodetic coordinates by repetitively calculating the relation between the digital elevation model (DEM) grid and the satellite position, which requires several processing iterations [8][9][10]. A good deal of work has been done using the traditional indirect geolocation method to achieve good accuracy. Geolocation accuracy can be effectively improved by using the precise latitude, longitude, and elevation provided by ground control points (GCPs), such as Wettzell corner reflector experiment. As reported in [11,12], TerraSAR-X has the capability of submeter-level geolocation accuracy in range. An in-depth analysis of TSX-1 carried out as part of the Wettzell corner reflector experiment, the pixel geolocation accuracy was improved to approximately 1 cm in range [11]. However, especially in Sensors 2019, 19, 5062 2 of 23 overseas, mountainous areas, canyons, deserts, etc., it is difficult to obtain GCPs, while the geolocation of these areas is very important. Therefore, a method of geolocation without GCPs is necessary in such cases. The RD model is the most widely-used physical sensor model for spaceborne SAR, which can be used to geolocate without GCPs [13]. In [14], a fast geocoding method of spaceborne SAR images using graphics processing units was proposed. This method accelerated the algorithm processing by performing block processing and parallel processing using the graphic processing unit or the central processing unit. However, this method is only an engineering implementation, and does not improve the calculation model. In [15], a hybrid algorithm based on the RD location model was proposed, which could achieve a higher level of accuracy than the analytic geolocation method (AGM) and relative geolocation method. This method took 10.9 s to process the first, middle, and last rows of a SAR image [15].
In this paper, a method to improve the Doppler centroid calculation model is proposed by estimating the geometric relationship between a target point inside a SAR image and the four corner points of the quadrangle SAR image. (For a N a × N r pixels SAR image, the four corner points are located at (1, 1), (1, N r ), (N a , 1) and (N a , N r ). The target point means any pixel except for the four corner points in the N a × N r SAR image.) This geolocation method for SAR images does not use GCPs, and since the iterative calculations of the Doppler centroid are omitted, it is more efficient to geolocate targets using this method while ensuring the geolocation accuracy. Furthermore, the proposed geolocation method is a general method that can be applied for all current spaceborne SAR systems. In this paper, we choose TanDEM-X images as the experimental data; the experimental results demonstrate the efficiency and reliability of the proposed method.
The rest of the paper is structured as follows. In Section 2 geolocation theory is introduced and divided into four subsections: geometric analyses for SAR images, the RD model for SAR images, the details about the IAGM, and the atmospheric propagation delay for microwaves that contains ionospheric and tropospheric delays. A fast, three-dimensional, indirect geolocation method is then proposed in Section 3. The experimental results in Section 4 validate the performance of the proposed method. Section 5 discusses the results and future works. Finally, Section 6 concludes the findings and contributions of this paper. Figure 1 illustrates the local coordinate system S − xyz based on the ECR coordinate system. The ECR coordinate system is widely used for RD model geolocation [16]. As shown in Figure 1, O − XYZ is a geocentric coordinate system, N represents the North Pole, and the arc NG is the Greenwich meridian. O is set as the origin, which coincides with the Earth's center of mass. The Z axis and the mean rotational axis of the Earth coincide; the X axis is pointing to the mean Greenwich meridian, and the Y axis is directed to complete a right-handed Cartesian coordinate system. E is the point where the Y axis passes through the surface of the Earth. S − xyz is the local coordinate system based on the ECR coordinate system. S represents the position of the satellite, and S is the sub-satellite point. T is the target on the Earth's surface, T is the projected point of T on the ellipsoid surface, T T is the elevation of the target T, and R is the slant range between the satellite and the ground target. S N is an arc passing S . S is set as the origin of the local coordinate system. The OS vector is regarded as the z axis; the direction of the OS vector is set as the positive direction. The x axis is in the OS N plane, which passes the subsatellite point S and is perpendicular to line S S. The y axis is then defined according to the right-handed rule of the Cartesian coordinate system. In addition, α is the angle for ∠TOS, and β is the angle between the plane OS N and the plane TOS . χ is the geocentric longitude of T, and φ is the geocentric latitude of T [16][17][18][19][20][21]. Sensors 2019, 19, 5062 3 of 23 Figure 1. The local coordinate system '− S xyz based on ECR coordinate system.

The Coordinate System and Transformation Formula
In Figure 1, the target T 's location is ( ) , , x y z . Its geocentric longitude and latitude in the spherical coordinate system are given by [22]: x y (2) in which χ T is the geocentric longitude of T , and φ T is the geocentric latitude of T . In the ellipsoidal coordinate system, the target T 's geodetic longitude and latitude are γ T and ϕ T , respectively. The target T 's geodetic longitude and geodetic latitude are the same [19]: The relation between geocentric latitude and geodetic latitude could be given by [22]: where c e is the first eccentricities of Earth, T h is the ellipsoid height of T . N is the radius of curvature in the prime vertical [19,22]: x y z are the parameters in the space rectangular coordinate system. The transformations from the geodetic coordinate to the space rectangular coordinate are given by [19,22]: In Figure 1, the target T's location is (x T , y T , z T ). Its geocentric longitude and latitude in the spherical coordinate system are given by [22]: in which χ T is the geocentric longitude of T, and φ T is the geocentric latitude of T.
In the ellipsoidal coordinate system, the target T's geodetic longitude and latitude are γ T and ϕ T , respectively. The target T's geodetic longitude and geodetic latitude are the same [19]: The relation between geocentric latitude and geodetic latitude could be given by [22]: where e c is the first eccentricities of Earth, h T is the ellipsoid height of T. N is the radius of curvature in the prime vertical [19,22]: (γ T , ϕ T , h T ) are the parameters in the geodetic coordinate system, and (x T , y T , z T ) are the parameters in the space rectangular coordinate system. The transformations from the geodetic coordinate to the space rectangular coordinate are given by [19,22]: 2.2. The RD Model for SAR Images

The Earth Model Equation
Since the Earth's shape is technically an oblate ellipsoid, the difference between its semimajor and semiminor axes is small. But for precise measurements, we must consider the model of the Earth as an oblate ellipsoid. Therefore, the Earth model equation is given by [23,24]: where a is the semimajor axis and b is the semiminor axis.
Considering the ellipsoid height, the Earth model equation can be rewritten as [25]:

The SAR Doppler Equation
The Doppler frequency in the ECR coordinate system for a generic point target is given by [23]: where f D is the Doppler frequency, λ is the radar wavelength, R S = (x S , y S , z S ) is the satellite position vector, R T = (x T , y T , z T ) is the ground target position vector, V S = .
x S , . y S , . z S is the satellite velocity vector, and V T is the ground target velocity vector. The velocity vector V T in the ECR coordinate system is given by [23]:

The SAR Range Equation
The slant range between the sensor and the ground target at given time t i is written as the following expression [23][24][25]: where R(i, j) = (x S − x T ) 2 + (y S − y T ) 2 + (z S − z T ) 2 , i represents the row number for the SAR image, and j represents the column number for the SAR image. j can be written as [25]: where R min is the nearest slant range for the SAR image, f s is the sampling frequency, and c is the speed of light.

The RD Model
The RD model is established between a satellite and a ground target on the Earth's surface [24], which includes Equations (8), (9), and (11). At every azimuth time, the position and velocity vectors of the satellite can be fitted though GPS ephemeris. Therefore, R S , V S , and V T are known, and the ground target R T is unknown, which can be calculated by solving the three equations using AGM.

Iterative Analytical Geolocation Method (IAGM)
As we know, the shape of the Earth is an ellipsoid. Point T is a small area around the subsatellite point S , so OS = OT = R L . R L is called the local radius, which is given by the following expression [16]: where φ S is the geocentric latitude of S . According to the cosine theorem, the angle α is given by [16]: where . R S is the magnitude of the sensor position vector.
The detailed steps of the IAGM are as follows: (a) Calculate the geocentric latitude and longitude of the sub-satellite point S : (b) In the local coordinate system, the vectors R S , V S , and R T in the ECR coordinate system can be transformed to R Sl , V Sl , and R Tl [16]: where z Sl = R S − R L . Equations (17)- (19) are the parameters for the Doppler equation in the local coordinate system S − xyz.
(c) Calculate the angle β from the Doppler equation.
In the local coordinate system S − xyz, the Doppler equation is determined by: Equations (17)- (19) can be directly substituted into Equation (22): where Suppose [16]: Therefore, Equation (23) can be rewritten by [16]: The above equation can be solved as [16]: The sign "±" depends on the SAR observation mode. If the look direction is right, "+" is chosen. Otherwise, "−" is chosen.
(d) Calculate the geocentric latitude and longitude for the ground point T by using the angle α and the angle β [16].
Transform geocentric latitude and longitude to geodetic latitude and longitude [19,22]: Therefore, the vector of the ground target is R T = (γ T , ϕ T , h T ). The iterative steps for this AGM are executed in order to obtain more precise result. The process steps are described as follows: 1.
Get R L , R S , and R(i, j) according to the above method.

3.
According to the steps from (a) to (d), the position vector R T = (γ T , ϕ T , h T ) can be obtained.

4.
Calculate the value R L = R T − h T . If R L − R L < 0.01, then stop the iteration and get the final result (γ T , ϕ T , h T ). Otherwise, let R L = R L and re-execute step 2.

Atmospheric Propagation Delay for Microwaves
The SAR signals travel through the atmosphere slower than they travel in a vacuum due to air refractivity. The distance ∆L is related to the atmospheric refractive index n(s), and can be written as [26]: where ∆Z is the zenith delay, and MF(θ) is the projection function. Generally, the incidence angle of θ for a SAR system is between 20 • and 55 • . A simple projection function is given by [26]: The zenith delay includes two delay parts: the tropospheric delay and the ionospheric delay [26,27]:

The Ionospheric Delay
The ionospheric refractive index is given by [28]: where N e is the electron density, B 0 is the magnitude of the magnetic field vector B 0 , f is the radar frequency in Hz, and the two constants C X and C Y are given by [28]: where e is the electron charge (e = 1.602 × 10 19 C), ε 0 is the permittivity of free space ε 0 = 8.854 × 10 −12 F/m , m e is the electron mass (m e = 9.109 × 10 −31 kg), µ 0 is the permeability of free space µ 0 = 4π × 10 −7 H/m . The first two terms in Equation (33) are regarded as the first order refractive index. Due to the much smaller magnitude of the third-and fourth-order terms, only the electron density within the ionosphere is considered [28]. So, Equation (33) can be rewritten as [28]: Therefore, the group refractive index is expressed as [28]: The ionospheric delay can be calculated by the total number of electrons on the propagation path [26]: where TEC is the total electron content. Therefore: The radar signal ionospheric delay ∆L iono converts the path delay from nadir to the path at a constant incidence angle θ.
The real TEC value can be downloaded from the ftp server of the international Global Navigation Satellite Systems (GNSS) service (IGS) [29].

The Tropospheric Delay
Since most SAR systems operate at frequencies below 40 GHz, and the refractivity using the millimeter-wave propagation model (MPM) is more or less constant, the refractivity can be considered frequency independent [5].
In the integral zenith model, the tropospheric delay is divided into hydrostatic, wet, and liquid components [30,31]: where s is the actual propagation path.
where k 1 = 0.776KPa −1 , k 2 = 0.716KPa −1 , and k 3 = 3.75 × 10 3 K 2 Pa −1 are the constants, P d is the partial pressure of dry air in Pa, e w is the partial pressure of water vapor in Pa, T is the absolute temperature in Kelvin, and W cl is the cloud water content in g/m 3 . The liquid delay is very small, so the last term of liquid water in the above equation can be neglected [7]. The tropospheric delay can be rewritten as [7]: where z 0 is the height of the Earth surface, z t is the height of the upper limit of the troposphere.

Real Data for the Ionosphere and Troposphere
The international GNSS service provides orbits, ionospheric delay, tropospheric delay, and other high-quality GNSS data products online in near real time for the whole world [32][33][34][35].
In order to define a global ionospheric model, the vertical total electron content is represented as a function of geocentric longitude, latitude, and time in the form of a raster grid. The global ionosphere maps are provided in the IONospheric EXchange (IONEX) format [29]. The time resolution of the maps is two hours, and the spatial resolution of each gird is 5 degrees in longitude and 2.5 degrees in latitude. Figure 2 shows the global TEC maps at 10:00:00 UTC on 20 December 2018. Since 1997, IGS has provided zenith path delay (ZPD) products, which are archived and can be accessed through the FTP server [36]. As shown in Figure 3, the ZPD products for 20 December 2018 are available at Beijing Fangshan (BJFS) station. Therefore, the BJFS data can be used to represent the ZPD of Beijing.

A Fast, Three-Dimensional, Indirect Geolocation Method
Geodetic coordinates, i.e., latitude, longitude, and height, are required by the traditional indirect geolocation method. Therefore, high precision DEM data or DSM data is the basic requirement for precise geolocation. In short, the Doppler centroid and slant range can be obtained by calculating the relationship between the target geodetic coordinates and the satellite position. Then, the pixel location in the azimuth and range can be obtained according to the Doppler centroid and slant range, respectively. The traditional indirect geolocation method requires iterative calculations to obtain the azimuth time, which is very time consuming. Therefore, in this paper, we introduce a fast indirect geolocation method to speed up the calculation of the azimuth time.
In order to clearly express the proposed algorithm and the process steps, a flowchart is presented in Figure 4. As shown in Figure 5, the geometric illustration of the estimated azimuth time is given. Since 1997, IGS has provided zenith path delay (ZPD) products, which are archived and can be accessed through the FTP server [36]. As shown in Figure 3, the ZPD products for 20 December 2018 are available at Beijing Fangshan (BJFS) station. Therefore, the BJFS data can be used to represent the ZPD of Beijing. Since 1997, IGS has provided zenith path delay (ZPD) products, which are archived and can be accessed through the FTP server [36]. As shown in Figure 3, the ZPD products for 20 December 2018 are available at Beijing Fangshan (BJFS) station. Therefore, the BJFS data can be used to represent the ZPD of Beijing.

A Fast, Three-Dimensional, Indirect Geolocation Method
Geodetic coordinates, i.e., latitude, longitude, and height, are required by the traditional indirect geolocation method. Therefore, high precision DEM data or DSM data is the basic requirement for precise geolocation. In short, the Doppler centroid and slant range can be obtained by calculating the relationship between the target geodetic coordinates and the satellite position. Then, the pixel location in the azimuth and range can be obtained according to the Doppler centroid and slant range, respectively. The traditional indirect geolocation method requires iterative calculations to obtain the azimuth time, which is very time consuming. Therefore, in this paper, we introduce a fast indirect geolocation method to speed up the calculation of the azimuth time.
In order to clearly express the proposed algorithm and the process steps, a flowchart is presented in Figure 4. As shown in Figure 5, the geometric illustration of the estimated azimuth time is given.

A Fast, Three-Dimensional, Indirect Geolocation Method
Geodetic coordinates, i.e., latitude, longitude, and height, are required by the traditional indirect geolocation method. Therefore, high precision DEM data or DSM data is the basic requirement for precise geolocation. In short, the Doppler centroid and slant range can be obtained by calculating the relationship between the target geodetic coordinates and the satellite position. Then, the pixel location in the azimuth and range can be obtained according to the Doppler centroid and slant range, respectively. The traditional indirect geolocation method requires iterative calculations to obtain the azimuth time, which is very time consuming. Therefore, in this paper, we introduce a fast indirect geolocation method to speed up the calculation of the azimuth time.
In order to clearly express the proposed algorithm and the process steps, a flowchart is presented in Figure 4. As shown in Figure 5, the geometric illustration of the estimated azimuth time is given. . The rectangular space coordinate for the four corner points according to Equations (5) and (6) are recorded as ( ) , , , , , , The detailed steps of the proposed fast, three-dimensional, indirect geolocation method are as follows: (a) Calculate the slopes of the latitude and the longitude between the first azimuth time and the last azimuth time at the nearest and furthest slant ranges, respectively.   . The rectangular space coordinate for the four corner points according to Equations (5) and (6) are recorded as ( ) , , , , , , The detailed steps of the proposed fast, three-dimensional, indirect geolocation method are as follows: (a) Calculate the slopes of the latitude and the longitude between the first azimuth time and the last azimuth time at the nearest and furthest slant ranges, respectively. Suppose the elevations are zero for the four corner points, which are marked C1, C2, C3, and C4. The corner point C1 refers to the pixel (1, 1) in the SAR image, whose geolocation point on the Earth's surface is (γ C1 , ϕ C1 , 0). The corner point C2 refers to the pixel (1, N r ) in the SAR image, whose geolocation point on the Earth's surface is (γ C2 , ϕ C2 , 0). The corner point C3 refers to the pixel (N a , 1) in the SAR image, whose geolocation point on the Earth's surface is (γ C3 , ϕ C3 , 0). The corner point C4 refers to the pixel (N a , N r ) in the SAR image, whose geolocation point on the Earth's surface is (γ C4 , ϕ C4 , 0). The rectangular space coordinate for the four corner points according to Equations (5) and (6) The detailed steps of the proposed fast, three-dimensional, indirect geolocation method are as follows: (a) Calculate the slopes of the latitude and the longitude between the first azimuth time and the last azimuth time at the nearest and furthest slant ranges, respectively.
Then, the average values for the slopes of latitude and the longitude are given by: (b) Calculate the distance between C1 and C3: (c) As point T's projection point is T , the positions of T and T are (γ T , ϕ T , h T ) and (γ T , ϕ T , 0), respectively. The rectangular space coordinates of point T and point T are (x T , y T , z T ) and (x T , y T , z T ), respectively. Then calculate the distance between C1 and T .
(d) Estimate point C g on the straight line C1C3, whose geodetic coordinate γ C g , ϕ C g , 0 is given by: The rectangular space coordinate for point C g is x C g , y C g , z C g . Therefore, R T C g and R C1C g can be expressed as: According to the cosine theorem, the angle ζ can be given: Estimate the projection point C v for point T on the straight line C1C g .
(e) From the ratio relationship, the azimuth time t i_est of point T can be derived: where t 1 is the SAR image generation starting time, t line is the imaging time interval of the SAR image, and int(·) returns the integer part of a decimal number.
According to the geometry estimation, the estimated azimuth time t i_est of point T can be obtained, though it has an error that needs to be corrected.
(f) Obtain point T's precise azimuth time t i . Calculate the Doppler frequency at the estimated azimuth time t i_est by Equations (9)- (11).
Similarly, we could obtain the Doppler frequencies at the first and last azimuth times for the point at the scene center. The Doppler frequency slope can be written as: Therefore, the precise azimuth time t i can be obtained by the second estimation.
(g) Calculate the slant range of point T.
where x S_t i , y S_t i , z S_t i is the satellite position at the azimuth time t i , ρ r is the range spacing.
The proposed fast, three-dimensional, indirect geolocation method contains seven steps (a~g). By calculating the spatial geometric relationship between the four corner points and point T, the estimated azimuth time t i_est can be obtained. Then, the precise azimuth time t i can be easily derived by Equations (58)-(60), and the pixel location in azimuth i is acquired. Finally, by calculating the slant range of point T, the pixel location in range j can be obtained.

Geolocation Results
The proposed algorithm is tested with the TanDEM-X's spotlight mode SAR image. All the parameters are listed in Table 1. The semimajor and semiminor axes are 6378137 m and 6356752.315 m, respectively. The speed of light is 2.99792458E + 5 km/s. First, we find the relevant DSM data according to the TanDEM-X SAR image, and ensure the gridded DSM is suitable to the resolution of the image. Then, all the steps are executed according to the algorithm process flowchart in Figure 4. Figure 6 shows the TanDEM-X SAR image geolocation results using the proposed algorithm in this paper. This figure gives the geolocation results in geodetic coordinates. Figure 6a provides the 0.69 m × 0.53 m gridded DSM, whose accuracy is the same as elevation 1 [37]. It should be noted that the height of DSM data is orthometric height [38]. The proposed algorithm is deduced based on an ellipsoidal Earth model, which means that the height is ellipsoid height. The geoid height is the difference between the ellipsoid and geoid [38]. Therefore, we assume the ellipsoid height, the orthometric height, and the geoid height are recorded as 'h', 'H', and 'N', respectively. Considering that this experimental area is small, a geoid height at the scene center can replace the other geoid heights in the scene. The geoid height can be easily obtained through the EGM2008 Geopotential Model [39]. The three-and the two-dimensional images are given in Figure 6b,c, respectively. In Figure 6b, the geodetic latitude, geodetic longitude, and ellipsoid height are included for every pixel. First, we find the relevant DSM data according to the TanDEM-X SAR image, and ensure the gridded DSM is suitable to the resolution of the image. Then, all the steps are executed according to the algorithm process flowchart in Figure 4. Figure 6 shows the TanDEM-X SAR image geolocation results using the proposed algorithm in this paper. This figure gives the geolocation results in geodetic coordinates. Figure 6a provides the 0.69 m × 0.53 m gridded DSM, whose accuracy is the same as elevation 1 [37]. It should be noted that the height of DSM data is orthometric height [38]. The proposed algorithm is deduced based on an ellipsoidal Earth model, which means that the height is ellipsoid height. The geoid height is the difference between the ellipsoid and geoid [38]. Therefore, we assume the ellipsoid height, the orthometric height, and the geoid height are recorded as 'h', 'H', and 'N', respectively. Considering that this experimental area is small, a geoid height at the scene center can replace the other geoid heights in the scene. The geoid height can be easily obtained through the EGM2008 Geopotential Model [39]. The three-and the two-dimensional images are given in Figure 6b,c, respectively. In Figure 6b, the geodetic latitude, geodetic longitude, and ellipsoid height are included for every pixel. Since the Qianxun Spatial Intelligence Inc. (Qianxun SI, Shanghai, China) offers a positioning service with centimeter-level accuracy (FindCM) for its customers, in this paper, the QianXun (QX) FindCM positioning results are used as the true positions. Two methods are used to compare the algorithm accuracy. Method 1 (M1) has no atmospheric propagation delay correction. Method 2 (M2) has atmospheric propagation delay correction, which includes tropospheric delay and ionospheric delay corrections. From Figure 2, it can be observed that the TEC map in Beijing is about 7.8 TECU. Since the Qianxun Spatial Intelligence Inc. (Qianxun SI, Shanghai, China) offers a positioning service with centimeter-level accuracy (FindCM) for its customers, in this paper, the QianXun (QX) FindCM positioning results are used as the true positions. Two methods are used to compare the algorithm accuracy. Method 1 (M1) has no atmospheric propagation delay correction. Method 2 (M2) has atmospheric propagation delay correction, which includes tropospheric delay and ionospheric delay corrections. From Figure 2, it can be observed that the TEC map in Beijing is about 7.8 TECU.
The value of 7.8 TECU is essentially extracted from the IONEX file [29]. Using TEC and other parameters in Equation (39), the ionospheric delay of the TanDEM-X SAR image can be obtained. Figure 3 shows that the ZPD of Beijing area is 2.368 m. Using Equation (42), the tropospheric delay can be acquired. The proposed method can be easily proved in M1 and M2, with M2 obviously being a better choice than M1. Figure 7 presents the geolocation results for the Banshan Pavilion at Western Hills in Beijing which is picked up from Figure 6. The image size of Figure 7 is 1024 × 1024 pixels. Figure 7a,b are the two-dimensional images for M1 and M2, respectively. In order to see the eight check points (P1~P8) at the Banshan Pavilion clearly, the red windows in Figure 7a,b are enlarged and shown in Figure 7c,d, respectively. The red windows in Figure 7a,b have the same longitude and latitude ranges. Compared with Figure 7c, the eight check points (P1~P8) at Banshan Pavilion in Figure 7d obviously shift several pixels longitudinally. The value of 7.8 TECU is essentially extracted from the IONEX file [29]. Using TEC and other parameters in Equation (39), the ionospheric delay of the TanDEM-X SAR image can be obtained. Figure 3 shows that the ZPD of Beijing area is 2.368 m. Using Equation (42), the tropospheric delay can be acquired. The proposed method can be easily proved in M1 and M2, with M2 obviously being a better choice than M1. Figure 7 presents the geolocation results for the Banshan Pavilion at Western Hills in Beijing which is picked up from Figure 6. The image size of Figure 7 is 1024 × 1024 pixels. Figure 7a,b are the two-dimensional images for M1 and M2, respectively. In order to see the eight check points (P1~P8) at the Banshan Pavilion clearly, the red windows in Figure 7a,b are enlarged and shown in Figure  7c,d, respectively. The red windows in Figure 7a,b have the same longitude and latitude ranges. Compared with Figure 7c, the eight check points (P1~P8) at Banshan Pavilion in Figure 7d obviously shift several pixels longitudinally.   Figure 6) in Figure 7a. (d) The enlarged image for the red window (Area 1 in Figure 6) in Figure 7b. Figure 8 illustrates the geolocation errors for the Banshan Pavilion with these two methods. Figure 8a,c illustrate the geolocation errors in the north, east, and vertical displacements for M1 and M2, respectively. The red line with circle markers represents the north displacement, the green line with rectangle markers represents the east displacement, and the blue line with asterisks represents the vertical displacement. In Figure 8b,d, the red lines with blue triangle markers represent the three-dimensional errors for M1 and M2, respectively. window (Area 1 in Figure 6) in Figure 7a. (d) The enlarged image for the red window (Area 1 in Figure  6) in Figure 7b. Figure 8 illustrates the geolocation errors for the Banshan Pavilion with these two methods. Figure 8a,c illustrate the geolocation errors in the north, east, and vertical displacements for M1 and M2, respectively. The red line with circle markers represents the north displacement, the green line with rectangle markers represents the east displacement, and the blue line with asterisks represents the vertical displacement. In Figure 8b,d, the red lines with blue triangle markers represent the threedimensional errors for M1 and M2, respectively.  Table 2 presents the geolocation accuracy for the eight corner points. The QX positioning results are taken as the true positions for the eight check points. By comparing the geolocation results of M1 (or M2) with QX, we can get the geolocation accuracy of the check points. As shown in Table 2, there are large displacements in east for M1, which are mainly due to atmospheric propagation delay in range. The displacements of the eight corner points in northing ( ΔN ), easting ( ΔE ), orthometric height ( ΔH ), and spatial ( Δ Spatial) are summarized in Table 2. The M1 geolocation errors are very large; the three-dimensional errors range from 3.76 m to 5.85 m. The M2 geolocation errors are smaller than those of M1, with three-dimensional errors ranging from 0.72 m to 1.60 m.  Table 2 presents the geolocation accuracy for the eight corner points. The QX positioning results are taken as the true positions for the eight check points. By comparing the geolocation results of M1 (or M2) with QX, we can get the geolocation accuracy of the check points. As shown in Table 2, there are large displacements in east for M1, which are mainly due to atmospheric propagation delay in range. The displacements of the eight corner points in northing (∆N), easting (∆E), orthometric height (∆H), and spatial (∆Spatial) are summarized in Table 2. The M1 geolocation errors are very large; the three-dimensional errors range from 3.76 m to 5.85 m. The M2 geolocation errors are smaller than those of M1, with three-dimensional errors ranging from 0.72 m to 1.60 m.  Figure 9 shows the geolocation results for the street lights and barriers located on Xingshikou Road in Beijing, taken from Figure 6. The image size of Figure 9 is 2048 × 2048 pixels. Figure 9a,b are the two-dimensional images for M1 and M2, respectively. The red windows in Figure 9a,b are enlarged, as shown in Figure 9c,d, respectively. So, the eleven check points for the street lights (P1~P5 and P8~P11) and the street barriers (P6, P7) located on Xingshikou Road in Beijing can be clearly seen. The red windows in Figure 9a,b have the same longitude and latitude ranges. Compared with Figure 9c, the eleven check points for the street lights (P1~P5 and P8~P11) and the street barriers (P6, P7) located on Xingshikou Road in Figure 9d have obviously shifted longitudinally by several pixels.  Figure 9 shows the geolocation results for the street lights and barriers located on Xingshikou Road in Beijing, taken from Figure 6. The image size of Figure 9 is 2048 × 2048 pixels. Figure 9a,b are the two-dimensional images for M1 and M2, respectively. The red windows in Figure 9a,b are enlarged, as shown in Figure 9c,d, respectively. So, the eleven check points for the street lights (P1~P5 and P8~P11) and the street barriers (P6, P7) located on Xingshikou Road in Beijing can be clearly seen. The red windows in Figure 9a,b have the same longitude and latitude ranges. Compared with Figure  9c, the eleven check points for the street lights (P1~P5 and P8~P11) and the street barriers (P6, P7) located on Xingshikou Road in Figure 9d have obviously shifted longitudinally by several pixels. Similar to Figure 8, Figure 10 shows the geolocation errors for the street lights and the street barriers located on Xingshikou Road with the two methods mentioned above. Figure 10a,c illustrates the geolocation errors in the north, east, and vertical displacements for M1 and M2, respectively. The red line with circle markers represents the north displacement, the green line with rectangle markers represents the east displacement, and the blue line with asterisks represents the vertical displacement. In Figure 10b,d, the red lines with blue triangle markers represent three-dimensional errors for M1 and M2, respectively. Similar to Figure 8, Figure 10 shows the geolocation errors for the street lights and the street barriers located on Xingshikou Road with the two methods mentioned above. Figure 10a,c illustrates the geolocation errors in the north, east, and vertical displacements for M1 and M2, respectively. The red line with circle markers represents the north displacement, the green line with rectangle markers represents the east displacement, and the blue line with asterisks represents the vertical displacement. In Figure 10b,d, the red lines with blue triangle markers represent three-dimensional errors for M1 and M2, respectively. Similar to Figure 8, Figure 10 shows the geolocation errors for the street lights and the street barriers located on Xingshikou Road with the two methods mentioned above. Figure 10a,c illustrates the geolocation errors in the north, east, and vertical displacements for M1 and M2, respectively. The red line with circle markers represents the north displacement, the green line with rectangle markers represents the east displacement, and the blue line with asterisks represents the vertical displacement. In Figure 10b,d, the red lines with blue triangle markers represent three-dimensional errors for M1 and M2, respectively.  Table 3 presents the geolocation accuracy for the street lights and street barriers. The QX positioning results are taken as the true positions for the eleven check points. By comparing the geolocation results of M1 (or M2) with QX, we can get the geolocation accuracy of the check points. As shown in Table 3, there are large displacements in east for M1, which are mainly due to atmospheric propagation delay in range. The displacements of the eight corner points in northing ( ΔN ), easting ( ΔE ), orthometric height ( ΔH ), and spatial ( Δ Spatial) are summarized in Table 3. The M1 geolocation errors are very large; the three-dimensional errors range from 3.43 m to 4.80 m. The M2 geolocation errors are smaller than those of M1, with three-dimensional errors ranging from 0.47 m to 1.69 m.   Table 3 presents the geolocation accuracy for the street lights and street barriers. The QX positioning results are taken as the true positions for the eleven check points. By comparing the geolocation results of M1 (or M2) with QX, we can get the geolocation accuracy of the check points. As shown in Table 3, there are large displacements in east for M1, which are mainly due to atmospheric propagation delay in range. The displacements of the eight corner points in northing (∆N), easting (∆E), orthometric height (∆H), and spatial (∆Spatial) are summarized in Table 3. The M1 geolocation errors are very large; the three-dimensional errors range from 3.43 m to 4.80 m. The M2 geolocation errors are smaller than those of M1, with three-dimensional errors ranging from 0.47 m to 1.69 m.

Computational Efficiency and Accuracy
The computing platform used was a Thinkpad E570 (Lenovo group co., Ltd., Beijing, China) with Core i5-7200U at 2.5GHz, 8GB memory, solid state drive and running the Windows 7 operating system. In addition, the processing software used was MATLAB R2016b (The MathWorks, Inc., Natick, MA, USA). There are two geolocation experiments: (i) The gridded DSM with 1024 × 1024 pixels, (ii) The gridded DSM with 2048 × 2048 pixels. All the geolocation experiments are tested with the TanDEM-X's spotlight mode SAR image. The main parameters of the SAR image are shown in Table 1.
The scale and processing time are shown in Table 4 to compare the proposed method (M2) to the traditional three-dimensional, indirect geolocation method. The traditional method needs to iterate different times according to different thresholds. In this experiment, the threshold of time difference is set to 0.0001 s. The proposed method is not an iterative method, so it is not necessary to set a threshold. It is apparent that the proposed fast, three-dimensional, indirect geolocation method is about 12 times faster than the traditional three-dimensional, indirect geolocation method.
The mean errors of the three-dimensional geolocation for check points are shown in Table 5. The geolocation accuracy of the traditional method is the same as that of the proposed method within the scale of 1024 × 1024 and 2048 × 2048.

Discussion
This paper proposes a fast, three-dimensional, indirect geolocation method using IAGM and DSM data for spaceborne SAR images. The traditional indirect geolocation method based on DSM or DEM data is computationally complex to obtain the accurate azimuth time and slant range, because several iterations must be executed [8][9][10]. The proposed fast method includes two steps to obtain the accurate azimuth time, while keeping the geolocation accuracy. By utilizing the geometric relation between the four corner points and the target, it is easier to estimate the azimuth time. It has to be mentioned that the proposed method can be widely used in any kind of terrain including flat ground, hills, mountains etc. DSM refers to the height of ground surface, buildings, bridges and trees. But DEM contains only elevation information of the terrain, and does not contain other surface information. If there is any need for geolocation regarding buildings, bridges, and other ground surface information, DSM data must be used.
A TanDEM-X SAR image is used to prove this fast, three-dimensional, indirect geolocation method without GCPs. Generally, the geolocation results have no atmospheric propagation delay correction; however, for some special requirements, such as accurate positioning, atmospheric propagation delay correction must be considered. Therefore, a fast, indirect geolocation method is used in such situations. In addition, M1 and M2 have the same computational complexity, and M2 only increases the correction for slant range. The three-dimensional errors in range for the eight corner points at the Banshan Pavilion are 3.76 m to 5.85 m and 0.72 m to 1.60 m for M1 and M2, respectively. For the street lights and barriers located on Xingshikou Road, their three-dimensional errors in range are 3.43 m to 4.80 m and 0.47 m to 1.69 m for M1 and M2, respectively. M1 is a general geolocation method that does not require additional data support. M2 includes atmospheric propagation delay correction, which needs the support of atmospheric data. It has been shown that this fast method can be generally used in spaceborne SAR images for indirect geolocation. It is about 12 times faster than the traditional three-dimensional, indirect geolocation method.
In the pursuit of accurate positioning, future research will consider the solid Earth tides, continental drift, and so on.

Conclusions
In this paper, a fast, three-dimensional, indirect geolocation method for spaceborne SAR images is proposed. It decreases the calculation time by avoiding several iterations, as used in the traditional indirect geolocating method. For precise geolocation of a SAR image, atmospheric propagation delays have been compensated, and the orthometric height of the DSM data has also been transformed to ellipsoidal height. Therefore, the geolocation could achieve accurate results within 2 m from the true positions for the check points. It has been shown that the proposed method is reliable and efficient for geolocation using TanDEM-X's spotlight mode SAR images. Moreover, this paper provides a method which can also be used in other spaceborne SAR systems, such as COSMO-SkyMed, Radarsat, Sentinel etc. The accurate azimuth time can be calculated with less than one pulse repetition time (PRT) with very few errors compared to the traditional three-dimensional, indirect geolocation method.