A Novel Topography Retrieval Algorithm Based on Single-Pass Polarimetric SAR Data and Terrain Dependent Error Analysis

: Polarimetric synthetic aperture radar (PolSAR) data provide an alternative way for topography retrieval, especially when limited PolSAR data are available. This article proposes a novel topography retrieval algorithm based on the Lambertian backscatter model that further improves the vertical precision of digital elevation model (DEM) generation and requires only one ﬂight. The key idea of the proposed algorithm is to avoid data ﬂuctuations caused by the ratio of the azimuth slope angle to the polarimetric orientation angle (POA). The previous research has conﬁrmed the feasibility of generating a DEM based on single-pass PolSAR data, but its effect on the quality of reference DEM has not been well-explained. To analyze this effect, a large number of experiments on DEM with different resolutions are conducted. In addition, an in-depth analysis of non-linear and terrain-dependent errors is performed. The L-band PolSAR data of NASA/JPL TOPSAR and ALOS-2 PALSAR-2 and interferometric SAR (InSAR) DEM data are used to verify the proposed algorithm. The experimental results show that PolSAR data can be used as an additional reliable information source for DEM fusion under certain conditions to improve the quality of public DEM.


Introduction
Polarimetric synthetic aperture radar (PolSAR) data provides another method of terrain reconstruction, especially when radar interferometry cannot be applied. In general, surface topography retrieval can be achieved using orthogonal two-pass or single-pass PolSAR data [1,2]. The polarimetric orientation angle (POA) is a critical polarimetric parameter in topography retrieval, and several classical POA estimation methods have been proposed [3][4][5][6][7]. Jin [8] proposed to use an adaptive threshold method and an image morphology thinning algorithm for linear texture to calculate the azimuth slope when limited single-pass PolSAR data are available. Chen et al. [9] found that the method in [8] was complex and computationally expensive and proposed another slope estimation approach that combines the topography-POA shift model and the refined Lambertian backscattering model, which is used in radarclinometry [10]. Furthermore, Li et al. [11] proposed a new single-pass PolSAR topography retrieval method based on the polarization-dependent intensity ratio, which proved to be more effective than the method in [9].
Instead of the POA, using the polarimetric signatures for the derivation of sea-ice surface topography was discussed in [12], and an accurate DEM of sea-ice surface topography was generated based on co-polarimetric (coPol) coherence. However, this method requires both PolSAR and interferometric SAR (InSAR) data, and only the sea-ice surface topography was verified.
The algorithm proposed in [11] still has two problems, although it can further improve the accuracy of the generated digital elevation model (DEM). First, when the incidence angle is smaller than the ground range slope angle, the azimuth slope direction cannot be correctly determined. Second, the ground range slope angle is easily affected by data fluctuations, especially when the surface of terrain is almost flat. The first problem can be solved by using a large incidence angle for steep relief variation areas, and the second problem can be solved by revising the slopes by the second-order weighted full multigrid (WFMG) algorithm [11]. However, the incidence angle of a SAR image is fixed, which inevitably limits the application range of the algorithm. This paper proposes a novel topography retrieval algorithm based on single-pass PolSAR data; the main contributions are as follows: (1) New expressions of azimuth slope and range slope are derived, which can effectively avoid the data fluctuation caused by the ratio of azimuth slope to POA; (2) The method proposed in [13] is used to correct the irregular points of the slope. The principle is that when the azimuth and ground range pixel size are the same, the integral of the slope along a closed-loop should be zero. After obtaining the azimuth and ground range slope, the method of solving the Poisson equation is the same as the weighted full multigrid algorithm described in [14]; (3) The terrain dependent error was fully analyzed. The experimental results show that the quality of the generated DEM is related to the terrain, and it is especially suitable for terrain with a small elevation variation range and that is free from the interference of vegetation and buildings; (4) The application potential of single-pass PolSAR data in DEM generation is deeply analyzed. The experimental results show that PolSAR data can be used as an additional reliable information source for DEM fusion under certain conditions to improve the quality of public DEM.

Background of POA
In the monostatic backscattering case, the 2 × 2 complex backscattering matrix can be expressed as [15]: where, for a reciprocal target matrix, S hv = S vh . Lee et al. [5,15,16] have indicated that the POA of terrain surface is geometrically related to the topographical slopes and radar look angle, and has a simple form of where ξ is the POA of the ground patch, ω is the azimuth slope angle, γ is the ground range slope angle, and φ is the radar look angle. The SAR imaging geometry of a tilted surface patch is shown in Figure 1a, and the geometric diagrams of the ground azimuth slope tan ω and the ground range slope tan γ are shown in Figure 1b. The projection of surface normal to the xz and yz planes, and N xz and N yz , which are the corresponding projection components. The blue-marked parts d x , d y , and d z in (b) are the changes along the positive direction of the xyz coordinate axis, respectively. The azimuth slope is negative and the ground range slope is positive.
The circular polarization algorithm (CPA) developed by Lee et al. [5,17,18] based on the terrain reflection symmetry hypothesis is one of the most popular algorithms for POA estimation. Based on this hypothesis, for a reflection symmetrical medium, the circular polarization parameter S RR S † LL is real and the circular polarization estimator is expressed as [17] Re(A) denotes the real part of A, and arctan(·) is computed in the range of [−π, π]. Thus, the POA range is limited to [−π/4, π/4], but the actual POA range is [−π/2, π/2], which indicates a phase wrapping problem in steep terrain. The ambiguous problem of POA estimation on steep terrain can be successfully solved by combining physical scattering characteristics. As for the natural terrain dominated by Bragg scattering, the unambiguous estimation of POA over steep topography can be expressed as [19], which can be called the vertical polarization dominated E = 0 deorientation algorithm (VEDA), where E is the Huynen parameter [20].
C denotes the deoriented Huynen parameter,S hh andS vv are elements after deorientation of complex backscattering S, and ξ is the unambiguous estimation of POA. However, it is not clear whether the VEDA is effective for L-band data, because only the real measurement data of the P-band are verified in [19]. To address this issue, an example with NASA/JPL airborne L-band PolSAR data from Camp Roberts is presented, as shown in Figure 2. The pixel size in the azimuth and ground range directions is 10 m. The DEM generated by C-band InSAR of the corresponding area in Figure 2a is shown in Figure 3a, which shows the roughness of the terrain. The POA estimated by the DEM is used as a reference, shown in Figure 3b. The POAs estimated from the L-band PolSAR data using the classic CPA algorithm and the VEDA algorithm are shown in Figure 3c,d, respectively. The L-band PolSAR data have been filtered using a 5 × 5 boxcar filter to reduce the speckle noise.
The comparison of Figures 2a and 3c show that the POA noise mainly appears in the lower-left corner and right edge of the image, where there are valleys and trees. This is because the volume scattering does not fit the tilted surface model [16]. The image in Figure 3d is smoother than that in Figure 3c, but there is still noise due to the volume scattering. An additional 3 × 3 boxcar filter has been applied in Figure 3d [19]. Next, the root mean square difference (RMSD) between the POA from InSAR and POLSAR data is calculated for the whole study area, with Figure 3b as the ground truth value. The results are shown in Table 1, which proves that the VEDA is effective for L-band PolSAR data.

Two-Dimensional Slope Estimation and Correction
To further solve the two problems in [11], that is, the fact that the azimuth slope direction cannot be uniquely determined when the range slope angle is greater than the incident angle and the ground range slope angle can be easily affected by data fluctuations, this paper proposes a novel two-dimensional slope estimation method.
For a terrain that fits the Lambertian backscatter model, the intensity of a SAR image and the slope of a local terrain can be expressed as [10,21,22], where K is a calibration constant of a SAR system, σ is a backscatter coefficient of the ground surface, and R g and R a indicate the ground range and azimuth resolution of a SAR image, respectively.
The estimate of the ground range slope angle in [11,23] is Obviously, the ground range slope angle is easily affected by data fluctuations, especially when the surface is almost flat. The new estimate of ground range slope angle, which is derived in detail in Appendix A, can be expressed as where P is given by and I(0, 0) can be estimated by the flat region in the image. The term |S hh − S vv | 2 in the coherency matrix shows a consistently increasing trend after compensation, so cos ω can be obtained by [23] cos where |S hh −S vv | 2 denotes the results after compensation. However, ω cannot be uniquely determined only according to (11). Therefore, combining (2) with the estimated γ c in (9), the corrected azimuth slope angle is expressed as where the ω s is Finally, ω c and γ c denote the azimuth slope angle and the ground range slope angle estimated by the proposed algorithm, respectively. The detailed processing flow is shown in Figure 4. Moreover, the improved algorithm also avoids the problem that the ground range slope in [9] cannot be uniquely determined due to the arccosine equation.  However, the intensity I(0, 0) still needs to be estimated. This is because the illposed slopes equation cannot be solved directly in the case of single-pass PolSAR data. Nevertheless, the following experimental data analysis shows that the cost is worthwhile.
In addition, since in areas with speckle noise and shadows, radar information cannot truly reflect the undulating characteristics of the actual terrain, this paper uses the method proposed in [13] to eliminate the irregular points of a slope. The principle is that if the grid spacing of the azimuth direction and the ground range direction is the same, the integral value of the slope along the closed loop should be zero.

Elevation Estimation and POA Correction
After estimating the azimuth and ground range slope, the surface elevation can also be estimated by solving the Poisson equation [1,2].
The Poisson equation can be expressed as, where ∇ 2 is the Laplacian operator ∂ 2 /∂x 2 + ∂ 2 /∂y 2 , φ(x, y) is the height value at (x, y), and the source function δ(x, y) on the right side is the surface curvature value at (x, y). Therefore, the discretized Poisson equation is as follows where δ(x, y) is given by and ∆ g (x, y) = R g tan(ω(x, y)), ∆ a (x, y) = R a tan(γ(x, y)).
The WFMG algorithm [14] is used to solve the Poisson equation, and the weight is obtained by the method proposed in [13]. The key idea of the WFMG algorithm is to convert low-frequency components of error into high-frequency components and then to remove them by the Gauss-Seidel relaxation method [24].
To obtain absolute rather than relative elevation values, at least one elevation tie-point must be known for each integrated slope profile [3]. The tie-points can be obtained by down-sampling the public DEM, such as in the shuttle radar topography mission (SRTM) DEM [11]. Then, the DEM can be obtained according to the terrain retrieval processing map, as shown in Figure 5.
Moreover, according to (2), the DEM retrieved from full PolSAR data can also be used to calculate the POA, which provides an effective method to improve the estimation precision of POA without using high-resolution DEM. This will be explained in Section 3.4.

Azimuth Slope
Ground Range Slope

Relative Elevation
Absolute Elevation Tie Points Figure 5. Single-pass PolSAR topography retrieval processing diagram after the azimuth slope and the ground range slope are estimated.

Algorithm Performance Analysis
First, the method proposed in [11] was used to generate the tie-point matrix and the C-band InSAR DEM presented in Figure 3a was downsampled to the resolution of 640 m, which was the same as the resolution in [11], as shown in Figure 6. Meanwhile, for better performance analysis, the DEM in Figure 3a is repeated in Figure 6a and a 5 × 5 boxcar filter is used to filter the L-band PolSAR data. Then, the DEM inversion was performed using the proposed method; the DEM map derived from L-band PolSAR data is shown in Figure 7. Compared with Figure 6a, it can be found that the main features of the topography can be reconstructed through single-pass PolSAR data.
Further, the differences between the PolSAR-derived DEM and the C-band AIRSAR InSAR-derived DEM were deeply analyzed. The scatter plot of the DEM maps obtained from the overall InSAR, as the ground truth value, and PolSAR data is shown in Figure 8a. The cumulative distribution function of the elevation difference is shown in Figure 8b and we can see that 90.35% of the points have an elevation difference of less than 18.9 m. Two profiles were extracted from C-band InSAR DEM and L-band PolSAR DEM along lines A and B in Figure 6a; the comparison results are presented in Figure 8c,d. The overall trend of elevation change is the same, but it can be found that the error is greater in steep areas, which can be obtained by combining the profile at the intersection of lines A and B. Moreover, for quantitative comparison with previous research results, the data involved are resampled to 36 m × 36 m resolution. The RMSD of slope and height between C-band AIRSAR InSAR DEM, as the ground truth value, and L-band PolSAR DEM are shown in Table 2 and Table 3, respectively. The comparison with the results of Chen's method [9] and Li's method [11] shows that the performance of the proposed method is effective. However, the terrain-dependent errors shown by Figure 8c,d require further analysis.

Non-Linear and Terrain Dependent Error Analysis
To better assess the impacts of terrain on the retrieval of azimuthal and ground slopes, the NASA/JPL airborne L-band PolSAR data from Camp Roberts was analyzed in detail, including catchments in multiple directions, as shown in Figure 9. In particular, the data of Figures 2 and 9 are referred to as PD1 and PD2, respectively.
Likewise, a 5 × 5 boxcar filter is used to filter the L-band PolSAR data and the data involved are resampled to 36 m × 36 m resolution. However, the case where the resolution of the tie-points is 90m is additionally analyzed, as shown in Figure 10. The inversion results corresponding to the two resolutions are shown in Figure 11 and the height and slope RMSD between the PolSAR-derived DEM and C-band InSAR DEM are given in Tables 4 and 5, respectively. The effectiveness of the proposed method is further verified by comparing Tables 4 and 5 with Tables 2 and 3.
However, the accuracy of PD2 is lower than that of PD1, although the processing methods and data sources are the same. Comparing Figure 2 and Figure 9, the terrain of PD2 includes buildings, highways, vegetation, etc., which is more complex than PD1. Therefore, it can be considered that the elevation accuracy of topography retrieval is related to the terrain. In addition, comparing Figure 11a and Figure 11b, it is also related to the resolution of the tie-points. This will be analyzed in detail in the next section.  For further analysis, three regions of interest (ROI) are selected, as shown in Figure 12a, and the tie-points are shown in Figure 12b-d with a resolution of 90 m. The red (ROI1)and green (ROI2)-marked areas correspond to relatively flat terrain in different directions and the blue (ROI3)-marked areas correspond to relatively steep terrain.
The topography retrieval results of the red-, green-, and blue-marked areas in Figure 12a are shown in Figure 13. Comparing the difference images of the three ROIs in the fourth column, it can be found that the error of steep terrain is greater. A careful analysis of the first two rows shows that the high error areas in the difference image are mainly distributed at the boundary positions of different terrains, which is due to the non-homogeneous scattering of the ground cover at the boundary. This indicates that the proposed method is more suitable for terrain dominated by a single scatterer.
Moreover, the elevation profiles marked by line A-line D in the three ROI regions were extracted, as shown in Figure 14. Combining line A and line D in Figure 13a and line A, line B, and line D in Figure 13e, it can be seen that the poor estimation results in a, j, b, e, and k in Figure 14 are due to the fact that these profiles all span multiple terrains.
The height and slope RMSD between the PolSAR-derived DEM and C-band AIRSAR InSAR DEM are given in Tables 6 and 7, respectively. It can be seen that ROI3, the bluemarked area, has the largest estimation error. Combining Figures 13 and 14, the estimation accuracy of DEM generated by PolSAR data is related not only to the terrain, but also to the homogeny of ground covering scattering. Therefore, these factors should be considered in practice.  Figure 13. The first and second columns are the corresponding C-band InSAR DEM and Pauli images, respectively, the third column is the DEM generated from PolSAR data and the fourth column is the difference image between the first and third columns. The first to third rows correspond to the red-, green-, and blue-marked areas in Figure 12a, respectively. Further, combining the results of the first to third rows, it can be seen that the error of flat terrain is smaller, and the error is larger at the boundaries of different terrains.   Figure 14. The first to third columns correspond to Figure 13a, e, and i, respectively, and the first to fourth rows correspond to the elevation profiles marked by lines A and D, respectively. It can be seen that the error depends on the terrain, and the error is larger where the elevation changes sharply. The more complex the terrain, the larger the error dynamic range.

Analysis of Low-Resolution DEM Matching High-Resolution PolSAR Data
To evaluate the effect of low-resolution DEM matching the high-resolution polarization SAR data, the proposed method was used to analyze the situation when tie-points had the resolutions of 90 m and 30 m, as shown in Figure 15. The inversion results corresponding to the two resolutions are shown in Figure 16, where it can be seen that the higher the resolution of the tie-point matrix, the higher the similarity between the PolSAR-derived DEM and InSAR-derived DEM. (c) (d) Figure 16. The first row is the DEM image derived from L-band PolSAR data and the second row is the scatter diagram between C-band InSAR DEM (truths) and PolSAR-derived DEM (estimates). Furthermore, the first and second columns correspond to the cases where the tie-points have 90 m and 30 m resolution, respectively. It can be seen that the higher the resolution of the tie-points, the closer the elevation estimate is to the ground truth.
A quantitative analysis was also conducted, and the height and slope RMSD between the PolSAR-derived DEM and C-band InSAR DEM are given in Tables 8 and 9, respectively. It should be noted that in order to analyze the high-resolution situation, the data involved are not resampled here, but a 5 × 5 boxcar filter is still used. Experimental results show that PolSAR data can integrate terrain textures into lowresolution DEMs and improve the resolution of DEMs without interpolation.

POA Correction Effect Analysis
According to (2), the DEM retrieved from PolSAR data was used to calculate the POA, and the results are shown in Figure 17a. Comparison of the results in Figure 17a and those in Figure 3c,d shows that after correction, the POA in Figure 17a is smoother and more similar to the POA obtained from the C-band InSAR DEM. Combining Figure 17b,d, 99.03% of the absolute values in the difference image are less than 39°, and most are within 20°. Furthermore, it can be clearly seen from Figure 17c that the estimated value and the real value are symmetrically distributed on both sides of the diagonal, and the point density in the middle pink area is the largest, indicating that the corrected POA is highly consistent with the POA derived from DEM.
The experimental results show that the DEM derived from L-band PolSAR data has high application potential for POA correction.

Application Potential of PolSAR Data in DEM Generation
To further verify the reliability of PolSAR data for DEM generation, a field scientific expedition was conducted in Salt Lake, Qinghai Province, China. The L-band fully PolSAR data obtained by ALOS-2 PALSAR-2 in this area in August 2020 was used for analysis. The RGB Pauli color-coded image is shown in Figure 18a and the aerial view of the same area obtained from Google Earth is shown in Figure 18b The ROI marked by the red box in Figure 18 is shown in Figure 19. Similarly, the ROI includes a variety of terrain.
The proposed method was used to analyze the situation when tie-points had resolutions of 90 m, and 30 m, as shown in Figure 20. The inversion results corresponding to the two resolutions are shown in Figure 21, and the height and slope RMSDs between the PolSAR-derived DEM and SRTM DEM are given in Table 10 and Table 11, respectively.  (c) (d) Figure 21. The first row is the DEM image derived from L-band PolSAR data and the second row is the scatter diagram between SRTM DEM (truths) and PolSAR-derived DEM (estimates). The first and second columns correspond to the cases where the tie-points have 90 m and 30 m resolution, respectively. It can be seen that the higher the resolution of the tie-points, the closer the elevation estimate is to the ground truth.
The experimental results of ALOS2-PALSAR2 and AIRSAR data show that the elevation accuracy of DEM generated from POLSAR data depends on terrain, and the elevation accuracy is higher for relatively flat terrain without vegetation and artificial target interference.

Discussion
The experimental results in Section 3.1 show that the error of topography retrieval is terrain-dependent. To better assess the impacts of terrain on the retrieval of azimuthal and ground slopes, the PD2 data from Camp Roberts was analyzed in detail. The experimental results in Section 3.2 show that the error of steep terrain is greater. Moreover, the estimation accuracy of DEM generated by PolSAR data is related not only to the terrain, but also to the homogeny of ground covering scattering. Therefore, these factors should be considered in practice.
To evaluate the effect of low-resolution DEM matching the high-resolution polarization SAR data, tie-point matrices with resolutions of 90 m and 30 m, respectively, are analyzed in depth in Section 3.3. It was found that the higher the resolution of the tie-point matrix, the higher the similarity between the PolSAR-derived DEM and InSAR-derived DEM. Further, it can be concluded that PolSAR data can integrate terrain textures into low-resolution DEMs and improve the resolution of DEMs without interpolation.
Moreover, the experimental results in Section 3.4 show that the DEM derived from L-band PolSAR data has high application potential for POA correction. This will contribute to the application of base POA.
To further verify the reliability of PolSAR data for DEM generation, a field scientific expedition was conducted in Salt Lake, Qinghai Province, China. Combining the experimental results of ALOS2-PALSAR2 in Section 3.5 and AIRSAR data, it can be further concluded that the elevation accuracy of DEM generated from POLSAR data depends on terrain and that the elevation accuracy is higher for relatively flat terrain without vegetation and artificial target interference.
In addition, we know that the surface of the moon is bare and dry, and there is no vegetation interference. Therefore, the technology for generating DEM based on PolSAR data may be applied to the topographic mapping of the lunar surface, especially lunar permanently shadowed regions (PSR). Then, similar to the idea of combining the Lunar Orbiter Laser Altimeter (LOLA) and SELENE Terrain Camera (TC) with data adopted in [25], it may be possible to combine LOLA with PolSAR data for topographic mapping of lunar PSR.

Conclusions
In this paper, a novel topography retrieval algorithm based on single-pass PolSAR data is proposed. First, a new ground range slope estimation equation is derived and the azimuth slope estimation is also improved. Therefore, data fluctuations caused by the ratio of the azimuth slope angle to the POA are avoided.
Furthermore, the L-band PolSAR data from AIRSAR and ALOS2-PALSAR2 containing a variety of terrain are deeply analyzed. We find that the quality of DEM generated is related to terrain, and the error will be greater at the boundary of different terrains. Therefore, in practice, more attention should be paid to ROIs where surface scattering is homogeneous and the elevation dynamic range is small.
In addition, a detailed comparative analysis of the enhancement effect of PolSAR data on the DEM is conducted. The analysis results demonstrate that PolSAR data can not only be used to generate DEM, but also can further improve DEM resolution under certain elevation accuracy. The application potential of DEM generated from PolSAR data for POA estimation is also analyzed.
In the future, different scattering models need to be optimized according to different terrain characteristics to further improve the precision of topography retrieval. Next, we set P = I(ω,γ) sin 2 φ I(0,0) cos(γ) cos(φ) , θ = γ + φ, and obtain sin 2 θ − P 1 − sin 2 θ = 0. (A2) Therefore, according to the root formula of the univariate quadratic equation, we can obtain Finally, the new estimate of the ground range slope angle can be expressed as