Iteration convergence condition modeling for single spaceborne SAR image direct positioning

Abstract Single SAR image direct positioning is to determine the ground coordinate for each pixel in the SAR image assisted with a reference DEM. During this procedure, an iterative procedure is essentially needed to solve the uncertainty in elevation of each pixel in the SAR image. However, such an iterative procedure may suffer from the problem of divergence in shaded and serious layover areas. To investigate this problem, we performed a theoretical analysis on the convergence conditions that has not been intensively studied till now. The Range-Doppler (RD) model was simplified and then the general surface is degenerated into a planar surface. Mathematical deduction was then carried out to derive the convergence conditions and the impact factors for the convergence speed were evaluated. The theoretical findings were validated by experiments for both simulated and real scenarios.


Introduction
Single SAR image positioning is to determine the ground coordinate for each pixel in the SAR image assisted with an external DEM, which is widely used in the SAR ortho-image generation, multibaseline InSAR, stereo-SAR, and so on. (Cheng 2004;Jiang 2012;Toutin 2004) SAR image positioning could be realized through either direct or indirect method (Cheng 2004;Jiang 2012). The indirect method starts from the ground coordinates of DEM resolution cells to solve the corresponding image coordinates, and it is currently the most widely used geocoding method as it does not require the iterative procedure. The direct method, which goes in the opposite direction, is to calculate the ground coordinate from the image coordinate of every pixel (Toutin 2004;Zhang 2013). During the direct positioning process, due to the uncertainty of the elevation of each pixel in the SAR image, an iterative procedure is needed, which suffers from the problem of divergence in shaded and serious layover areas (Cheng 2004;Yang et al. 2006). As a good supplement to the indirect method, the direct method has explicit geometric significance which could be helpful to study the geometric meaning of the rigorous imaging model of the SAR image and develop the multi-image positioning method Cheng et al. , 2013. In the optical photogrammetry area, Sheng (2005) has given the convergence condition of direct method through theoretical analysis, while the convergence condition suitable for SAR images have not been intensively studied by now.
Similar to the collinearity equations in photogrammetry (Zhang and Zhang 1996), Range-Doppler (RD) model is the rigorous imaging geometry model for the SAR image relating the image coordinates with the ground coordinates, which was first proposed by Brown (1981) and developed by Curlander (1982). In this paper, we borrow some ideas from optical photogrammetry and derive the convergence condition for single SAR image direct positioning based on a simplified RD model, which can help us evaluate the applicability of the algorithm. Also with the known terrain conditions, the iteration convergence condition can be used to predict the divergent areas. With the theoretical analysis, we can also sort out factors that affect the convergence rate, especially the initial elevation. These theoretical results are tested with various scenarios of simulated or true slopes.

Theoretical analysis of convergence condition
In order to theoretically analyze the convergence conditions, the RD model is simplified and then the general surface is degenerated into a planar surface.

OPEN ACCESS
Mathematical deduction is carried out to derive the convergence conditions and the impact factors for the convergence speed are analyzed.

Basic principle for direct single SAR image positioning
For a particular image pixel, we use the RD model to calculate its ground coordinates with the given initial elevation, interpolate its height in the georeferenced DEM data and then update the initial elevation with the interpolated one. This procedure is repeated until the elevation difference between the two adjacent iterations is less than a certain threshold, and then the three-dimensional ground coordinates of the pixel are determined. However, as an iterative algorithm, it inevitably suffers from the divergence problem, especially in the areas with steep terrain. Suppose the position vector of satellite is R S = (X S , Y S , Z S ) T , the corresponding velocity vector is V SC = (X v , Y v , Z v ) T , and the position vector of the ground target T in the ground orthogonal coordinate system of geocentric space is R TC = (X T , Y T , Z T ) T , which satisfies the earth ellipsoid equation (Equation (1)).
where R p is the earth ellipsoid polar radius; R e is the earth ellipsoid equatorial radius; H T means the elevation of ground target T and it changes with the position coor- (1) If the target point T at the surface of the earth, then H T equals 0 m; and parameters R p and R e are determined by the specific ellipsoid model. Every pixel in the SAR images can be indexed by column j and row i. The azimuth slow time t ij when the radar signal arrives at the ground target can be determined by t ij = i·δ ta , where i is the row number and δ ta means the imaging time of each line in azimuth direction. In Equation (2), R (i, j) means the distance between satellite and ground target at the moment t ij . Equation (3) makes use of the Doppler information at the moment t ij (Cheng 2004;Sekhar, Kumar, and Dadhwal 2014). The received radar signal suffers from Doppler Effect due to the non-orthogonal angle between the flight direction and echo direction. The Doppler frequency f d can be acquired through the meta data of the SAR image. Here in our experiment, we use the TerraSAR data which have been zero-Doppler processed and hence f d = 0. Equations (1)−(3) comprise the RD model together (Cheng 2004;Sekhar e tal. 2014).
There are four unknown variables in above equations which are the position coordinates of target T (X T , Y T , Z T ) and its height H T , but there are only three equations to solve. Therefore, an iterative procedure is needed. As shown in Figure 1, it is a zero Doppler profile and the "S" represents where the SAR satellite is. Using a given initial elevation H (0) T , we can determine the position coordinate of the initial projected point A 1 (X (1) T , Y (1) T , Z (1) T ) in the object space by Equations (1)-(3). Next, vertically trace from A 1 to surface point B 1 , and interpolate its height H (1) T from the DEM using the planar position coordinates of A 1 . Update H T with H (1) T and then use Equations (1)- . And then repeat the same procedure to determine subsequent projected points A 2 , A 3 , A 4 , … , A n until A n is convergent to target point T. Finally, output the position vector of A n as the ground coordinates of T. If the target T loactes at profile of back-radar slope, the projected points A 1 , A 2 , … , A n approach point T from both sides (as shown in Figure 1(a)), while if it locates at Profile of toward-radar slope, the projected points A 1 , A 2 , … , A n approach point T from only one side (as shown in Figure 1(b)).

Analysis of the iterative procedure
It is difficult to analyze the convergence condition on the general surface. In order to make the theoretical analysis easier, in this paper we first simplify the RD model, then discuss the convergence condition and further derive the convergence speed on planar surface, and finally analyze those on general surface. Figure 2 shows a zero Doppler plane orthogonal to the satellite flight direction. S is the sensor, T is the ground target, and the flight direction of sensor is perpendicular to the paper inwardly. MA 1 represents the earth ellipsoid surface and MT represents the simulated slope. In order to facilitate the modeling analysis, we assume S as the origin of the local coordinate system, where x is the direction from sensor S to Earth ellipsoid center O, z is the flight direction of the satellite, and y completes the right hand system S-xyz.
The RD model can be simplified in three aspects. First, MA 1 can be viewed as line segment MA 1 because its length is relatively very short compared with the distance between Earth center O and the first projected point A 1 . Second, A 1 A 2 can be replaced by line segment A 1 A 2 because calculation for actual data shows that length of A 1 A 2 is relatively very short compared with the distance between the sensor S and A 1 . Third, ∠SA 1 A 2 is viewed as a right angle due to the same reason as in the second assumption above.
In order to validate the model simplification, we select the pixel located in the 10th row and 10th column in Spotlight TerraSAR-X image of Mount Song area acquired on 25 September 2011. According to the calculation, we can obtain MA 1 /OA 1 = 0.0000386, A 1 A 2 /SA 1 = 0.000247, ∠SA 1 A 2 = 90.028 • . Suppose ∠MA 1 A 2 = and θ represents the radar depression angle. For the back-radar slope, there is = ∕2 − (as shown in Figure 2(a)); for the toward-radar slope, = ∕2 + (as shown in Figure 2(b)). Convergence on the planar terrain surface is analyzed from the viewpoint of geometry. Suppose the distance between T and the nth projected point A n is d n and β represents the slope angle in zero-Doppler plane.
If the geometric sequence d n is convergent, A 1 , A 2 , … , A n approaches T gradually. Hence, the sufficient and necessary condition for convergence is β ≤ γ. Also = ∕2 − . Finally, the convergence condition for back-radar slope is + ≤ ∕2.
For the toward-radar slope, ∈ (0, ∕2) , ∈ ( ∕2, ), d 2 = d 1 tan ( − )∕ tan , d 3 = d 2 tan ( − )∕ tan and so on. In a general form, Hence, d 1 , d 2 , d 3 , …, d n is a geometric sequence with a positive ratio of tan ( − )∕ tan . According to the convergence condition of geometric sequence, we can obtain If the geometric sequence d n is convergent, A 1 , A 2 , …, A n approaches T gradually and iteration is convergent. Hence, the sufficient and necessary condition for convergence is − ≥ . Also = ∕2 + . Finally, the convergence condition for toward-radar slope is + ≤ ∕2. When the sum of the slope angle and the depression angle is greater than ∕2, that is + > ∕2 , the layover is generated. Hence, the non-layover areas are convergent.
For back-radar slopes, ∈ (0, ∕2), we can assume ∈ (0, ∞); for toward-radar slopes, ∈ ( ∕2, ), assume ∈ (−∞, 0). Then d n = d 1 (tan ∕ tan ) n−1 (n ≥ 2) fits for both back-and toward-radar slopes. When the iteration is convergent, assume the convergence threshold is ΔG and d n <ΔG. The initial elevation is H (0) T and the actual elevation for target T is H real T . We may assume where [x] is the truncation function which means taking the smallest integer greater than x. According to Equation (6), the convergence speed is influenced by slope angle β in zero Doppler plane, the radar depression angle θ, convergence threshold ΔG and initial elevation   As indicated by the results in Table 1, the iteration convergence can be achieved only when the slope angle and depression angle satisfy the deduced convergence condition. The theoretical analysis and experimental results agree with each other perfectly. The iterations required for convergence increase with the slope angle. Generally when the slope angle in the zero Doppler plane is about 30°, the convergence speed will be very slow. The initial elevation in these experiments is set to 0 m, that is H T (0) = 0 m. In fact, initial elevation has a great influence on the convergence speed. Taking the back-radar slope with slope angle of 10° as an example, the pixel selected in the TerraSAR-X image is still (10, 10) and its elevation converges to 38.622 m.
As shown in Figure 3, the height range 0-200 m are divided by interval of 10 m and every number of this sequence is set to be the initial elevation for iteration

Experiments on planar surface
The experimental data used was a TerraSAR-X Spotlight image of Mount Song area acquired on 25 September 2011. We randomly selected a pixel with line and column number (10, 10) and set the initial elevation H (0) T and θ = 61.038°. Assuming a planar surface with different slope angles for the specified pixel (10, 10), an iterative positioning process is implemented. If convergent, the iterations and the final elevation that the iteration converges to were recorded. The sum of the slope angle and the depression angle were calculated and compared with π/2. If it was greater than π/2, the iteration process will be divergent according to the theoretical convergence condition. The toward-and back-radar slopes were tested differently. 90° which does not meet the convergence condition derived above. (3) The general surface is varied and rolling and it is very different from the simulated planar surface. Therefore, the theoretical convergence condition cannot fit for every ground point. Also the iteration process is limited by the resolution of the auxiliary DEM.

Conclusions
When geocoding a single SAR image with a DEM, for simplified planar surface, the convergence of iteration is depended on both the slope angle and the depression angle. The convergence speed is determined by the depression angle, the slope angle, convergence threshold, and initial elevation. The closer the distance between the initial elevation and the true elevation, the faster the iteration converges. When the initial elevation deviates a lot from the true elevation, the convergence speed will be slow and stable. When the true elevation of ground target is unknown, it is reasonable to set the initial elevation as 0 m. The simplification of the RD model is reasonable which is verified by the experiment. Hence, it can be used to deduce the convergence condition. The experiment results of a planar surface fit well with the theoretical convergence condition, which validate the theoretical findings. For a general surface, the divergent points occur mainly in the layover and shaded areas as the theoretical analysis has expected. However, it cannot be guaranteed that the convergence of all the points is consistent with the theoretical prediction because the general surface is complicated and projected points A 1 , A 2 , … , A n fall on different slopes while the convergence condition is deduced based on the planar surface.
As for the general ground surface, the convergence condition can also be used by dividing the complicated general surface into small locally flat grounds, in which the convergence condition should be the product of convergence conditions with different slope angles and depression angles. in this experiment. We can see that when the initial elevation is 40 m, the convergence speed is the fastest. Likewise, 37-40 m are divided by interval of 0.1 m and when the initial elevation equals 38.6 m, the convergence speed is the fastest; and then the height range 38.61-38.63 m are divided by interval of 0.001 m and when the initial height are set between 38.619 and 38.626 m, the convergence speed is the fastest. The experimental results show that the closer the distance between the initial elevation and the final elevation that the iteration converges to, the faster the convergence speed is.
It can also be seen in Figure 3, when the initial elevation is larger than 60 m (until 200 m), the convergence speed is stable and slow. That is to say, if the initial elevation differentiates substantially from the final converge to elevation, the convergence speed will be slow and stable. When we know nothing about the true elevation, we can set the initial elevation as zero m, which would not have much impact on the convergence speed.

Experiments on general surface
The experimental data included a Stripmap mode TerraSAR-X image of Mount Song area acquired on 18 July 2011 and a reference DEM with the resolution of 3 m covering the same area. In this experiment, 20 typical feature points are chosen to test the theoretical convergence condition, including ten points for toward-radar and ten for back-radar slopes, respectively. The experimental results were shown in Tables 2 and 3, where the slope angle and look angle of each sample points were calculated to test if the theoretical prediction meets the real iteration results. Also the screenshot of each sample point is obtained and the red stars represent the sample points, which show the terrain conditions (flat, layover, or shaded place) of the sample points directly.

Analysis of experimental results
(1) If the iteration is convergent, for back-radar slope, the projected points A 1 , A 2 , … , A n approach the ground point T from two sides; for toward-radar slope, the projected points approach T from one side.
(2) For the general surface, iteration is divergent in layover and shaded areas. In layover areas, multiple points on the ground become one point in the SAR image and hence when iterating, the projected points oscillate and cannot converge to the ground target T. In shaded areas, echo signal information is inadequate. From the other side, the slope angle is greater than the depression angle in shaded areas. For TerraSAR-X data, the depression angle is about 60° and the sum of the slope angle and the depression angle is obviously greater than

Notes on contributors
Yuting Dong is a PhD candidate of LIESMARS, Wuhan University. Her research interest is the geolocation of high resolution SAR images and radargrammetry.
Lu Zhang is Professor of LIESMARS, Wuhan University. His research focuses on radar interferometry, radargrammetry, and geolocation of high resolution SAR images.
Mingsheng Liao is Professor of LIESMARS, Wuhan University. He is mainly engaged in time-series InSAR and their applications in topography, geological hazard monitoring, infrastructure security etc.