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.

Keywords:

• Direct positioning
• RD model
• SAR
• iteration convergence condition

1. 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, stereoSAR, 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 and Zhang 2012; Cheng et al. 2012, 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.

2. 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. Mathematical deduction is carried out to derive the convergence conditions and the impact factors for the convergence speed are analyzed.

2.1. 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(1) )).(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 coordinates of target T (X_T, Y_T, Z_T) that is H_T = f (X_T, Y_T, Z_T). 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(2) ), R (i, j) means the distance between satellite and ground target at the moment t_ij. Equation (3(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(1) )−(3(3) ) comprise the RD model together (Cheng 2004; Sekhar e tal. 2014).(2) (3)

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 , we can determine the position coordinate of the initial projected point A₁ (, , ) in the object space by Equations (1(1) )–(3(3) ). Next, vertically trace from A₁ to surface point B₁, and interpolate its height from the DEM using the planar position coordinates of A₁. Update H_T with and then use Equations (1(1) )–(3(3) ) to solve the projected point A₂ (, , ). And then repeat the same procedure to determine subsequent projected points A₂, A₃, A₄, … , 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₁, A₂, … , 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₁, A₂, … , A_n approach point T from only one side (as shown in Figure 1(b)).

Figure 1. Iterative positioning method. (a) Profile of back-radar slope; (b) Profile of toward-radar slope.

2.2. 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₁ 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.

Figure 2. Simplification of the RD model. (a) Back-radar slope; (b) Toward-radar slope.

The RD model can be simplified in three aspects. First, MA₁ can be viewed as line segment MA₁ because its length is relatively very short compared with the distance between Earth center O and the first projected point A₁. Second, A₁A₂ can be replaced by line segment A₁A₂ because calculation for actual data shows that length of A₁A₂ is relatively very short compared with the distance between the sensor S and A₁. Third, 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₁/OA₁ = 0.0000386, A₁A₂/SA₁ = 0.000247, . Suppose and θ represents the radar depression angle. For the back-radar slope, there is (as shown in Figure 2(a)); for the toward-radar slope, (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.

For the back-radar slope, , , d₂ = d₁ tan β/tan γ, d₃ = d₂ tan β/tan γ, and so on. In a general form, d_n = d_n–1 tan β/tan γ = d₁(tan β/tan γ)^n–1. Hence, d₁, d₂, d₃, … , d_n is a geometric sequence with ratio of tan β/tan According to the convergence condition of geometric sequence, it is easy to obtain.(4)

If the geometric sequence is convergent, A₁, A₂, … , A_n approaches T gradually. Hence, the sufficient and necessary condition for convergence is β ≤ γ. Also . Finally, the convergence condition for back-radar slope is .

For the toward-radar slope, , , , and so on. In a general form, . Hence, d₁, d₂, d₃, …, d_n is a geometric sequence with a positive ratio of . According to the convergence condition of geometric sequence, we can obtain(5)

If the geometric sequence is convergent, A₁, A₂, …, A_n approaches T gradually and iteration is convergent. Hence, the sufficient and necessary condition for convergence is . Also . Finally, the convergence condition for toward-radar slope is . When the sum of the slope angle and the depression angle is greater than , that is , the layover is generated. Hence, the non-layover areas are convergent.

For back-radar slopes, , we can assume ; for toward-radar slopes, , assume . Then 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 and the actual elevation for target T is . We may assume and derive n by solving this following equation.(6)

where [x] is the truncation function which means taking the smallest integer greater than x. According to Equation (6(6) ), the convergence speed is influenced by slope angle β in zero Doppler plane, the radar depression angle θ, convergence threshold ΔG and initial elevation .

3. Experiments

3.1. 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 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.

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 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.

Table 1. Experimental results of a planar surface.

β(°)	β + θ(°)	Back-radar slope		Toward-radar slope
		Iterations	Converge to (m)	Iterations	Converge to (m)
5	66.038	6	43.624	6	58.559
10	71.038	8	38.622	8	70.877
15	76.038	12	34.538	12	90.509
20	81.038	19	31.093	19	127.510
25	86.038	36	28.104	36	225.739
30	91.038	246	25.450	199	1257.358
40	101.038	Divergent	Divergent	Divergent	Divergent
50	111.038	Divergent	Divergent	Divergent	Divergent
60	121.038	Divergent	Divergent	Divergent	Divergent
70	131.038	Divergent	Divergent	Divergent	Divergent
80	141.038	Divergent	Divergent	Divergent	Divergent
90	151.038	Divergent	Divergent	Divergent	Divergent

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 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.

Figure 3. Relationship between initial elevation and convergence speed.

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.

3.2. 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.

Table 2. Experimental results of back-radar slope points.

Screenshot of Sample points	(Column, Row)	Slope angle β	Look angle θ	β + θ	Iteration result
	(800, 331)	35.962°	63.339°	99.301°	Divergent
	(236, 2153)	41.568°	63.045°	104.613°	Divergent
	(246, 2136)	66.121°	63.047°	129.168°	Divergent
	(994, 2200)	76.530°	63.038°	139.568°	Divergent
	(818, 2706)	8.633°	62.957°	71.59°	Convergent
	(512, 2704)	9.552°	62.957°	72.509°	Convergent
	(1145, 1382)	16.440°	63.169°	79.609°	Convergent
	(748, 508)	22.121°	63.310°	85.431°	Convergent

Table 3. Experimental results of toward-radar slope points.

Screenshot of Sample points	(Column, Row)	Slope angle β	Look angle θ	β + θ	Iteration result
	(832, 262)	36.469°	63.350°	99.819°	Divergent
	(323, 1648)	48.322°	63.126°	111.448°	Divergent
	(1366, 406)	61.025°	63.327°	124.352°	Divergent
	(432, 3072)	66.915°	62.898°	129.813°	Divergent
	(107, 172)	3.096°	63.364°	66.46°	Convergent
	(1939, 675)	7.786°	63.284°	71.07°	Convergent
	(1288, 2517)	13.284°	62.987°	76.271°	Convergent
	(1337, 2573)	19.205°	62.978°	82.183°	Convergent

3.3. Analysis of experimental results

(1)	If the iteration is convergent, for back-radar slope, the projected points A1, A2, … , An 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 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.

4. 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₁, A₂, … , 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.

Funding

This work was supported by the National Natural Science Foundation of China [grant number 41271457] and the Demonstration System of High Resolution Remote Sensing Applications in Urban Fine Management Area [grant number 06-Y30B04–9002-13/15].

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.

Acknowledgments

The authors would like to thank the German Aerospace Center (DLR) for providing the test data-sets via the DLR AO LAN0793 and LAN0634, and Prof. Miaozhong Xu of LIESMARS for providing the photogrammetric DEM with spatial resolution of 3 m.