Next Article in Journal
Interannual Monitoring of Cropland in South China from 1991 to 2020 Based on the Combination of Deep Learning and the LandTrendr Algorithm
Next Article in Special Issue
Estimating Urban Forests Biomass with LiDAR by Using Deep Learning Foundation Models
Previous Article in Journal
A Multi-Sensor Approach to Characterize Winter Water-Level Drawdown Patterns in Lakes
Previous Article in Special Issue
Evaluating Predictive Models of Tree Foliar Moisture Content for Application to Multispectral UAS Data: A Laboratory Study
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Underlying Topography Estimation over Forest Using Maximum a Posteriori Inversion with Spaceborne Polarimetric SAR Interferometry

1
Key Laboratory of Technology in Geo-Spatial Information Processing and Application System, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100190, China
2
Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China
3
School of Electronic, Electrical and Communication Engineering, University of Chinese Academy of Sciences, Beijing 100049, China
*
Author to whom correspondence should be addressed.
Submission received: 16 January 2024 / Revised: 6 March 2024 / Accepted: 6 March 2024 / Published: 8 March 2024

Abstract

:
This paper presents a method for extracting the digital elevation model (DEM) of forested areas from polarimetric interferometric synthetic aperture radar (PolInSAR) data. The method models the ground phase as a Von Mises distribution, with a mean of the topographic phase computed from an external DEM. By combining the prior distribution of the ground phase with the complex Wishart distribution of the observation covariance matrix, we derive the maximum a posterior (MAP) inversion method based on the RVoG model and analyze its Cramer–Rao Lower Bound (CRLB). Furthermore, considering the characteristics of the objective function, this paper introduces a Four-Step Optimization (FSO) method based on gradient optimization, which solves the inefficiency problem caused by exhaustive search in solving ground phase using the MAP method. The method is validated using spaceborne L-band repeat-pass SAOCOM data from a test forest area. The test results for FSO indicate that it is approximately 5.6 times faster than traditional methods without compromising accuracy. Simultaneously, the experimental results demonstrate that the method effectively solves the problem of elevation jumps in DEM inversion when modeling the ground phase with the Gaussian distribution. ICESAT-2 data are used to evaluate the accuracy of the inverted DEM, revealing that our method improves the root mean square error (RMSE) by about 23.6% compared to the traditional methods.

1. Introduction

According to statistics, the world has a total forest aera of 4.06 billion hectares (ha), which is 31% of the total land area [1]. The acquisition of forest understory digital elevation model (DEM) has been a challenging problem. Accurate DEM data can assist decision makers in effectively planning forest land utilization, afforestation, and logging activities, thereby achieving the conservation and optimal utilization of forest resources [2]. Meanwhile, the forest understory DEM is crucial for identifying landslides and debris flows [3]. Through the analysis of accurate DEM data, we are able to pinpoint the location and path of potential landslides and debris flows, aiding in the implementation of effective disaster prevention measures.
Light Detection and Ranging (LiDAR) and Polarimetric Interferometric Synthetic Aperture Radar (PolInSAR) are the primary means in remote sensing for obtaining DEM under the forest. LiDAR can penetrate the forest canopy and has good accuracy in estimating DEM [4,5,6]. However, it is limited by the atmosphere, mist and clouds. Due to high cost and spatial discontinuity, it is difficult to achieve global and long-term coverage [7]. PolInSAR combines the sensitivity of InSAR to the spatial and height distribution of scatterers with the sensitivity of PolSAR to the shape and orientation of scatterers. It has advantages including all-weather and all-season functionality, as well as high temporal and spatial resolution, making it widely used in obtaining DEM [8,9,10,11,12,13,14].
To obtain the required physical parameters, it is necessary to establish a model that correlates the PolInSAR data with these parameters. In recent years, more and more models have emerged for retrieving vegetation parameters from PolInSAR data. Yamada [15] assumed that there were two uncorrelated scattering centers in the forest area, located at the top of the canopy and on the ground, respectively. This assumption transformed the problem into the direction-of-arrival (DOA) estimation. The interference phase of the ground and canopy can be obtained by ESPRIT (Estimation of Signal Parameters via Rotational Invariance Technique) [16,17,18,19]. However, because of the strong volume scattering component, the detection accuracy of this technique becomes worse for dense forest regions [20,21]. Based on the theory of electromagnetic scattering, Treuhaft et al. proposed the Random Volume over Ground (RVoG) model [22,23]. The model simplifies the forest region into a layer of randomly oriented uniform particles on the ground and derives the expression for complex interference coherence. It is a function of four physical parameters: the forest height that determines the thickness of the volume layer; the mean extinction coefficient that represents the attenuation of electromagnetic waves through the canopy; the ground phase related to the underlying topography; and the ground-to-volume amplitude ratio (GVR) that varies with polarization channels.
Solving the RVoG model parameters is a difficult problem, especially for the ground phase, which is the basis for solving the other parameters and an important parameter for inverting the DEM. In the RVOG model, the complex interference coherence can be expressed as a line in the complex plane. This line intersects the unit circle with two points, one of which represents the ground phase and the other has no practical physical meaning. This phenomenon is known as the double-candidate effect of ground phase. Based on this geometric property, Cloude et al. [24] propose a three-stage inversion (TSI) method. This method divides the inversion process into three stages: (1) least squares line fit; (2) vegetation bias removal; (3) height and extinction estimation. From the perspective of estimation theory, Tabb et al. [25,26] proposed using a maximum likelihood estimation (MLE) to retrieve forest parameters. This approach begins with the observation that the covariance matrix of the PolInSAR data follows a complex Wishart distribution. By maximizing the log-likelihood function, an objective function related solely to the ground phase is obtained. MLE provides the ability to completely separate volume scattering from surface scattering. However, the objective function regarding ground phase in MLE usually presents two indistinguishable peaks, corresponding to the two intersections of the coherence line and the unit circle in the TSI, making this method also suffers from the double-candidate effect.
Due to the double-candidate effect, the solution of the ground phase becomes complex. The three-stage method often uses some assumptions based on physical scattering to solve this problem, such as selecting the intersection point that is far away from the complex coherence under HV polarization as the ground phase [24]. However, these methods do not fully meet the assumptions in practice. Therefore, the accuracy of parameter inversion is limited. Some researchers propose solving the ground phase using the relationship between the vertical wavenumber (kz) and two candidate phases [27], which performs well with airborne data but is less effective with certain spaceborne data. For MLE, Huang et al. [28] proposed the maximum a posteriori inversion method for the RVoG model by modeling the ground phase as a Gaussian distribution (MAPG) with the mean of the topographic phase computed from the external DEM. This method effectively suppressed the double-candidate effect. However, the assumption of Gaussian distribution cannot satisfy the circular data, which may lead to incorrect results when the topographic phase appears near the phase jump point.
To make the MAP method more robust, this paper proposes a maximum a posteriori with the von Mises distribution as the prior (MAPV) method for solving the RVoG model. The method models the ground phase as a Von Mises distribution [29,30,31], which effectively addresses the problems caused by modeling the ground phase with the Gaussian distribution, making the estimation of the ground phase more continuous and accurate. Additionally, since the solution method for MAP requires exhaustive searching of the phase, this will significantly decrease computational efficiency. This paper proposes a four-step optimization (FSO) method for solving ground phase based on the characteristics of the MAPV objective function, employing a gradient-based optimization approach [32,33]. This method markedly improves solving efficiency without sacrificing accuracy.
This paper is organized as follows. Section 2 provides a detailed description of the experimental materials. Section 3 introduces the RVoG model and provides a brief overview of the three-stage inversion process. Section 4 presents the proposed MAPV method, and analyzes and simulates its CRLB. Furthermore, a four-step optimization method is proposed. Results and discussions are presented in Section 5. Finally, conclusions are given in Section 6.

2. Materials

The test area is located in the central region of Sardinia ( 39 48 N and 8 45 E), which is in the southwest of Italy. The area is heavily forested, with the predominant vegetation being holm oak woods. The test area consists of mountains and plains, with ground elevations ranging from 0 to 760 m above mean sea level. Figure 1a depicts the Google Earth optical image of the area projected onto the SAR coordinate system, with the dark green area indicating vegetated regions.
Paired quad-polarized data used for interferometry were acquired by the SAOCOM-1A and SAOCOM-1B satellites, respectively. The SAOCOM (Satélite Argentino de Observación COn Microondas) satellite series is developed by the Argentine Space Agency and comprises two satellites: SAOCOM 1A and SAOCOM 1B. This mission provides L-band (approximately 1.275 GHz) full polarimetric data with spatial resolutions ranging from 10 to 100 m in both real-time and stored modes, and an incidence angle between 20 to 50°. Figure 1b displays a Pauli-based PolSAR image of the test area, with the RGB channels representing HH-VV, HV, and HH+VV, respectively.
Height checkpoints are from Ice, Cloud and land Elevation Satellite-2 (ICESat-2) data [34]. This satellite carries the Advanced Topographic Laser Altimeter System (ATLAS), which introduced single photon detection technology for the first time in Earth’s elevation measurement, significantly improving the data acquisition rate for terrain detection. This article employs the land and vegetation height data (ATL08) from ICESat-2, with a vertical uncertainty of 0.2 m for flat terrain and 2.0 m for mountainous terrain [35].
The external DEMs selected for the experiment are the ALOS Global Digital Surface Model “ALOS World 3D - 30 m” (AW3D30) and the Shuttle Radar Topography Mission 30 m (SRTM30), both datasets having an approximate resolution of 30 m. Figure 1c,d show the Alos DEM projected to the SAR coordinate system and the topographic phase, respectively.
Figure 2 shows the flowchart of the proposed DEM inversion method. The entire inversion framework consists of (a) ground phase estimation based on MAPV and (b) DEM generation.

3. Scattering Model

3.1. PolInSAR Data Description

A single-baseline PolInSAR system acquires the complex scattering matrix of each resolution element via dual-polarized antennas ( h , v ) from two slightly different look angles. According to the reciprocity ( S h v = S v h ), the scattering matrix can be expressed in the form of Pauli vector [36], defined as
k i = 1 2 S h h i + S v v i S h h i S v v i 2 S h v i T , i = 1 , 2
where the S p q represents the scattering coefficients, p , q { h , v } , and the superscript ( · ) T denotes the vector transpose operation.
The 6 × 6 covariance matrix can be defined as
R ^ = k 1 k 2 k 1 H k 2 H = T 1 ^ Ω ^ Ω ^ H T 2 ^
where ( · ) H is the conjugate transpose, and · is the average operation in data processing. T 1 ^ and T 2 ^ contain the polarimetric information of the two images, respectively, called the polarimetric coherence matrix, while Ω ^ contains the polarimetric and interferometric information of the two images, called the polarimetric inteferometry matrix.
For any given non-zero scattering mechanism ω , the complex interferometric coherence can be expressed as
γ ( ω ) = ω H Ω ^ ω ω H T ^ ω
where T ^ = ( T 1 ^ + T 2 ^ ) / 2 .

3.2. The Random Volume over Ground (RVoG) Model

Without considering other sources of decorrelation, such as system and temporal factors, and focusing solely on the decorrelation effect induced by scatterers, the complex coherence (3) can be expressed as the integral ratio form of the vertical distribution function F ( z ) of the scatterers [37]
γ ( ω ) = e j k z z 0 0 h v F ( z ) e j k z z d z 0 h v F ( z ) d z
where z 0 is the vertical position of the ground surface, h v reperesents the forest height, and k z denotes the wavenumber.
The RVoG model assumes that the forested area is composed of randomly oriented uniform particles, with an exponential vertical distribution function, and that the ground layer is impenetrable with a delta function for its vertical distribution, defined as
F ( z ) = m v ( ω ) e 2 σ z cos θ 0 + m g ( ω ) δ ( z )
where m v and m g are the volume and the ground scattering amplitudes, respectively. θ 0 is the incidence angle, and σ represents the mean extinction coefficient. By inserting (5) into (4), we can obtain
γ ( ω ) = e j ϕ 0 γ v + μ ω 1 + μ ω = e j ϕ 0 ( γ v + μ ω 1 + μ ω ( 1 γ v ) )
where γ v denotes the volume coherence, ϕ 0 denotes the ground phase, and μ ( ω ) represents the GVR.
Notice that (6) represents the complex interference coherence as a straight line on the complex plane. The three-stage inversion method mentioned above is based on this geometric property. This method fits the complex coherence into a straight line that intersects the unit circle at two points. The ground phase is then selected from these two points based on the feature that the direct ground scattering signal is weak in the L-band HV channel [24].

4. Maximum a Posteriori Estimation of the Ground Phase

Given the observed value R ^ , the maximum a posteriori (MAP) estimate of the parameter ϕ can be expressed as
ϕ ^ MAP = arg   max ϕ p ( ϕ | R ^ ) arg   max ϕ p ( R ^ | ϕ ) p ( ϕ ) arg   max ϕ log p ( R ^ | ϕ ) + log p ( ϕ )
where p ( R ^ | ϕ ) is the likelihood function, which represents the performance of the data set, and p ( ϕ ) represents a prior distribution of ϕ .
It is known that the observed coherence matrix of PolInSAR after multilook (2) follows a complex wishart distribution
p ( R ^ ; R , N ) = c ( R ^ ) | R | N exp ( N · t r ( R 1 R ^ ) )
where N is the number of looks, c ( R ^ ) is a constant for normalization, | · | is the determinant operator, and t r is the trace operator.
A priori knowledge about the ground phase is obtained by introducing an external DEM.

4.1. A Prior Probability Model of the Ground Phase

With InSAR geometry, the topographic phase can be solved by an external DEM [38]
ϕ topo = 4 π B λ R sin θ h .
The ground phase in PolInSAR is the same as the topographic phase, both representing the contribution of the underlying ground to the interferometric phase. Considering the influence of vertical accuracy and resolution, there is an error between the topographic phase calculated from an external DEM and the actual ground phase. The ground phase can be represented as:
ϕ = ϕ topo + ϕ e
where ϕ e denotes the error component. Some researchers modeled the ϕ e as a Gaussian distribution with zero mean and variance dependent on the external DEM. Therefore, the ground phase follows a Gaussian distribution with a mean of the ϕ topo . The maximum a posteriori estimation method was used to solve the RVoG model parameters. This method effectively suppressed the double-candidate effect, making the two local maxima of the objective function clearly distinguishable [28].
However, since the estimated ground phase is the principal value of the phase within the [ π , π ) interval after the periodic wrapping, π and π differ numerically by 2 π but represent the same point in the complex plane. This point is called the phase jump point. The assumption of the Gaussian distribution does not take into account the characteristic of phase wrapping, so it is discontinuous near the phase jump point, which can lead to errors in phase estimation.
Figure 3a,c show the Gaussian distribution curves in the Cartesian coordinate system and on the unit circle of the complex plane, respectively, for one phase period. In the complex plan, it is evident that among the two candidates ( c 1 and c 2 ), c 2 is closer to ϕ topo , so it should have a higher probability of being the ground phase. However, in Figure 3a, the probability of c 1 being the ground phase is greater than that of c 2 , which is obviously incorrect.
The modeling of the ground phase must consider the characteristic of phase wrapping while retaining the selectivity of the Gaussian distribution. At present, probability models have not found much application to circular data. Common ones include: Uniform, Cardioid, Wrapped Normal, Wrapped Cauchy, and Von Mises distributions [29,39]. This paper employs the Von Mises distribution to model ground phase, offering significant advantages in this context. Firstly, its circular nature enables better adaptation to periodic data, akin to Cardioid and Wrapped Normal distributions, but more flexible than Wrapped Cauchy distribution. Secondly, the Von Mises distribution possesses parametric flexibility, enabling adjustment of its parameters to accommodate different degrees of data concentration. Furthermore, the mathematical form of the Von Mises distribution is simple and easy to understand, facilitating straightforward computation and interpretation, while also possessing robust statistical properties such as mean, variance, etc., making it more manageable and analyzable in practical applications [30,40].
The probability density function of the Von Mises distribution is
p ( ϕ ; ϕ topo , κ ) = e κ cos ( ϕ ϕ topo ) 2 π I 0 ( κ )
where κ represents the degree of concentration (similar to the inverse of variance), I 0 ( κ ) is the modified Bessel function of the first kind of order 0, which is used for normalization.
As can be seen in Figure 3b,d, the Von Mises distribution effectively solves the problem of discontinuity at phase jump point in the Gauss distribution, making the ground phase estimation more robust.

4.2. MAP with Von Mises Distribution as Prior

By changing m v ( ω ) and m g ( ω ) in (5) to matrix form
m v ( ω ) = ω H T v ω m g ( ω ) = ω H T g ω
and inserting them into (4), we can obtain the matrix form of the RVoG model in the same form of (3), and
T = I 1 T v + a T g Ω = e j ϕ ( I 2 T v + a T g )
where a , I 1 , I 2 are functions of the forest height h v and extinction coefficient, T v and T g are the volume and ground coherence matrices, respectively. The volume only coherence can be obtained by γ v = I 2 / I 1 .
It is worth noting that R ^ , T ^ , and Ω ^ are the ensemble averages of the data and should be distinguished from the ideal mathematical expectations R , T , and Ω .
Combining the complex Wishart distribution (8) and the prior probability distribution (11), the objective function of the MAPV estimation of the RVoG model is given by (7)
f MAPV = log p ( R ^ ; R , N ) + log p ( ϕ ; ϕ topo , κ )
where the first term is the objective function of the MLE. Maximising it yields [25]
log p ( R ^ ; R , N ) = N · ( 3 log ( 1 cos θ ) log   | A θ + ϕ |     log   | A ϕ | )
where
A α = T ^ 1 2 ( e j α Ω ^ + e j α Ω ^ H )
and θ represents the angle of the intersection of the unit circle with the extension of the line segment connecting 1 and γ v in the complex plane, geometrically.
Substituting (11) and (15) into (14), the objective function with respect to the ϕ and θ is obtained as
f MAPV = 3 log ( 1 cos θ ) log   | A θ + ϕ |     log   | A ϕ | + κ N cos ( ϕ ϕ topo ) .
Differentiating the objective function with respect to θ and ϕ , respectively, and setting them to 0 yields
f θ = 3 sin θ 1 cos θ d d α log   | A α | = 0 f ϕ = d d α log   | A α | | A ϕ | | A ϕ | κ N sin ( ϕ ϕ topo ) = 0 .
Combining the two equations of (18), we can obtain the expression of θ with respect to ϕ
θ = 2 tan 1 3 | A ϕ | | A ϕ | + κ N sin ( ϕ ϕ topo ) .
Figure 4 shows the difference between the objective functions of MAPG and MAPV. By observing the positions of ϕ topo , c 1 , and c 2 in the complex plane in Figure 3, it can be seen that the expected peak of the objective function should appear near c 2 , which is consistent with the MAPV objective function.
According to (14) and (19), the ground phase is given by
ϕ 0 = arg max ϕ ( 1 cos θ ) 3 | A θ + ϕ | | A ϕ | e κ N cos ( ϕ ϕ topo )
which is related to the underlying topography. The DEM can be obtained from part (b) of Figure 2.

4.3. The Cramer–Rao Lower Bound Analysis

To evaluate the validity of the method, it is necessary to analyze its Cramer–Rao Lower Bound (CRLB), which provides the minimum bound of variance if it is unbiased.
Note that the above derivation of MAPV is done under the noise-free assumption. In practice, the PolInSAR system inevitably introduces noise, the noise-affected R can be expressed as [41]
R n = T Ω Ω H T + N 1 0 0 N 2
where N 1 and N 2 are the noise covariance matrices of the receivers at each end of the baseline. Assuming that the noise between channels is uncorrelated, N i ( i = 1 , 2 ) can be expressed as
N i = 1 2 σ h h i 2 + σ v v i 2 σ h h i 2 σ v v i 2 0 σ h h i 2 σ v v i 2 σ h h i 2 + σ v v i 2 0 0 0 σ h v i 2 + σ v h i 2
where σ p q 2 , p , q { h , v } denotes the noise variance of each channel. Note that since the noise of the two receivers is uncorrelated, the noise has no effect on Ω .
The CRLB can be derived from the inverse of the Fisher information. The derivation of the Fisher information is presented in Appendix A for clarity, which can be expressed as [42,43]
I F = N t r R n 1 R n ϕ R n 1 R n ϕ + κ I 1 ( κ ) I 0 ( κ ) .
Therefore, the CRLB is
Var ϕ ^ 1 I F = 1 N t r R n 1 R n ϕ R n 1 R n ϕ + κ I 1 ( κ ) I 0 ( κ ) .
Figure 5 shows the simulation results using the data from [44] with N = 50 and κ = 3.65 (approximately equivalent to 30 ).
It is worth noting that the variance of the TSI method falls below the theoretical CRLB curve when the SNR is reduced below about −14 dB. This is due to the fact that the phase error is calculated in the sense of phase wrapping, whereas the CRLB is calculated using the local curvature of the distribution function [45]. From Figure 5, it can be seen that the inclusion of the prior effectively reduces the CRLB, and also demonstrates the superiority of the MAPV method over the TSI method for ground phase estimation.

4.4. Four-Step Optimized Solution for MAPV

To determine the ground phase, which corresponds to maximizing the objective function (17), the typical approach involves conducting an exhaustive search (ES) over ϕ . Although this ensures accuracy in finding the solution, it leads to high computational complexity and low computational efficiency.
Since the objective function of MAPV with respect to ϕ is a continuous curve, as shown in Figure 4. Based on this property, this paper employs the gradient optimization method to enhance the efficiency of solving ground phase. The gradient of the MAPV objective function is presented in Appendix B.
Given the significant variations in the objective function for each pixel, the method requires addressing two distinct problems:
  • Initial phase: Due to the presence of two peaks in the objective function, a single initial point is prone to fall into the local maximum. Therefore, it is necessary to choose an appropriate strategy to ensure that the results converge to the global maximum.
  • Global learning rate: Because of the large differences between the objective functions of different pixels, it is very important to choose the global learning rate so that the method can adapt to all objective functions.

4.4.1. Initial Phase

Taking into account the characteristics of the MAPV objective function depicted in Figure 6, this paper introduces a four-step solving approach.
  • Step 1: Gradient descent using ϕ topo as initial phase. Get the ϕ vally between the two peaks.
  • Step 2: Give the ϕ vally a Δ ϕ in the opposite direction to ϕ topo as a ϕ seed .
  • Step 3: Perform gradient ascent separately using this ϕ seed and ϕ topo to obtain the two peak values and their respective phases.
  • Step 4: Compare the magnitudes of these two peak values and select the phase corresponding to the larger peak value as the final result.

4.4.2. Global Learning Rate

In the FSO method, we employed both gradient ascent and gradient descent algorithms. For convenience, this paper converts the gradient ascent process into a gradient descent process, with the specific conversion procedure given below.
The gradient descent algorithm is based on the Root Mean Square Propagation (RMSProp) [33]. Although RMSprop adjusts the learning rate during iterations, the global learning rate determines the speed, convergence, and eventual performance of the model during training. If the learning rate is too high, it may lead to an unstable training process and difficulties in convergence. Conversely, a low learning rate might result in slower training speed, requiring more iterations to converge to an optimal solution.
Due to significant differences in the objective functions of various pixels, selecting a uniform global learning rate poses a challenge. Therefore, this paper proposes improvements to RMSprop to reduce the model’s reliance on the global learning rate to accommodate all objective functions, as shown in Algorithm 1.
When the current gradient g has opposite signs with the minimum gradient g m i n , or when the current function value f ( ϕ ) is greater than the minimum function value f m i n , it indicates that the step size might be too large, causing oscillations. In such cases, the minimum phase is necessarily situated between ϕ and ϕ m i n , closer to the side with the lower function value. In the fourth line of Algorithm 1, a weighted computation based on function values is applied for this condition, assigning higher weight to the lower function value when updating the phase. As the updated phase is very close to the minimum phase, it is advisable to moderately decrease the learning rate for better convergence.
Algorithm 1 assumes the objective function to be positive, but neglecting certain constant terms during derivation may result in the objective function potentially becoming negative. This could result in an error on line 4 during execution. This problem can be solved by preprocessing the coherence matrix, let
T = E 1 2 T E 1 2 Ω = E 1 2 Ω E 1 2
where E is the diagonal matrix formed by the eigenvalues of T .
Algorithm 1: Framework of improved RMSprop algorithm.
   Require: Global learning rate ϵ , Decay rate ρ , Initial phase ϕ
   Initialize: r 0 , f m i n inf , ϕ m i n 0 , g m i n 0 , n 0 , δ 10 6
Remotesensing 16 00948 i001
In FSO, the gradient ascent process can be converted into a gradient descent process. In the gradient ascent task, calculate the function value f and gradient g at that point, then let
g = g / f 2 f = 1 / f .
Then, Algorithm 1 can be used for solving.
Although the method requires performing three iterations of gradient descent, the improvement made in RMSprop based on the characteristics of the objective function allows each gradient descent to converge to the extremum within ten iterations. Furthermore, the precision obtained by this algorithm often surpasses that achieved through exhaustive search (<1 degree).

5. Results

5.1. Evaluation Indicator

The evaluation indicators are mean error m, root mean square error (RMSE) r, and accuracy α σ
m = 1 N n = 1 N ( x ˜ n x n ) r = 1 N n = 1 N ( x ˜ n x n ) 2 α σ = 1 N n = 1 N I | x ˜ n x n | σ
where I | x ˜ n x n | σ is the idicative function, defined as
I | x ˜ n x n | σ = 1 , if | x ˜ n x n |     σ 0 , otherwise .

5.2. Efficiency of FSO for MAPV

As mentioned above, the discrepancy in the objective function between pixels poses a challenge in the selection of the global learning rate when using traditional RMSprop. Algorithm 1 proposed in this paper is an improvement specifically designed to address this issue, and it is necessary to discuss its convergence here.
Figure 7 presents the convergence performance of RMSprop and Algorithm 1 under different global learning rates. As depicted in Figure 7a,d, under conditions of a relatively small global learning rate, both methods exhibit robust convergence. However, this comes at the cost of requiring more than 20 iterations. As the global learning rate increases, Algorithm 1 demonstrates a pronounced acceleration in convergence. In contrast, the RMSprop fails to converge to extrema due to oscillations.
While a larger ϵ appears to contribute to faster convergence in this experiment, it is important to note that, due to the potentially steep peaks in the MAPV objective function, excessively high global learning rates may cause the algorithm to jump out of these peaks. Therefore, it is also not advisable to set the global learning rate too high.
Table 1 presents the efficiency of the four-step optimization method for the ground phase of the test area. Here, ’ES’ denotes achieving 1-degree accuracy through exhaustive search. It is evident that the proposed method is 5.6 times faster than exhaustive search while ensuring accuracy. On average, convergence is achieved in only 24.5 steps per pixel.

5.3. Results of DEM Inversion in the Test Area

Given that the Gaussian distribution is not suitable for modeling the ground phase, MAPG encounters discontinuities when solving for ground phase. These discontinuities are propagated to the DEM during phase inversion, resulting in discontinuities within the DEM. This section will employ the Alos DEM as a benchmark to evaluate the extent of this discontinuity. Despite vertical accuracy differences between Alos DEM and actual values, it remains an essential reference in assessing DEM continuity.
As shown in Figure 8, the phase estimated by MAPG and MAPV show similar fringes that are basically consistent with the topographic phase computed by the external DEM. However, the comparison shows that there are significant differences between the two methods in the vicinity of the phase jump point (the blue and red boundary region in the phase diagram). In this region, the MAPG typically transitions with an intermediate value, resulting in a phase discontinuity that causes elevation jumps when inverting the DEM.
The DEM inverted by the two methods are shown in Figure 9a,b. Due to the large elevation range in the test area, it is not easy to observe differences over an entire image. Therefore, we select two areas for comparison, corresponding to the two areas in Figure 8. Different scales are set for the two areas. As shown in Figure 9, the MAPG-inverted DEM shows a strong discontinuity in some regions that are highly correlated with the vicinity of the phase jump point in Figure 8a. From the elevation change curves of the marked line segments in Area 1 and Area 2 shown in Figure 10a,b, it can be seen that the jump between adjacent pixels in the MAPG-inverted DEM reaches 20 m in some areas, which does not correspond to the actual situation. On the other hand, the MAPV-inverted DEM shows better results without significant elevation jumps.
Table 2 evaluates the continuity of Area 1 and Area 2. The criterion used is the accuracy α , where discrepancies greater than the set threshold compared to the Alos DEM are deemed as discontinuities. It is evident that compared to the MAPG-inverted DEM, the MAPV-inverted DEM exhibits better continuity.

5.4. Performance Assessment of the Proposed Method

In this section, IceSat-2 data are used to confirm the effectiveness of MAPV in generating DEM for forested areas.
Figure 11 shows the distribution of IceSat-2 data in the test forest area, where the color of the dots indicates the magnitude of the difference between the calculated elevation at that point and the LiDAR data. Figure 11a,b show the difference between the Alos and SRTM DEM with the LiDAR DEM, respectively. With reference to the colorbar on the right, it can be seen that the Alos and SRTM DEM are generally higher than the LiDAR DEM. Figure 11c shows the difference between the TSI-inverted DEM and the LiDAR DEM. It is worth noting that, as mentioned above, the assumptions made by TSI in solving the double-candidate effect cannot satisfy all conditions, and errors will occur in some low-coherence regions, resulting in elevation jumps in these regions when inverting the DEM. Figure 11d shows the difference between the MAPV-inverted DEM and the LiDAR DEM, and it can be seen that the DEM obtained by this method is closer to the LiDAR DEM.
Table 3 shows the mean error (ME) and the root mean square error (RMSE) of the different source DEM compared to the LiDAR DEM. Compared to the DEM of Alos and SRTM, the RMSE of the MAPV-inverted DEM is improved by 23.1% and 24.1%, respectively. It should be noted that due to the elevation jumps in the TSI-inverted DEM, the elevation in some areas may differ significantly from the LiDAR DEM, resulting in a large RMSE.
Figure 12 shows a histogram of the distribution of elevation differences between different source DEM. Obviously, compared to the DEM of Alos, SRTM, and TSI, the error center of MAPV-inverted DEM is closer to 0.
The method proposed in this paper validates the feasibility of using spaceborne L-band PolInSAR for retrieving forest understory DEM. The rapid development of spaceborne PolInSAR technology has expanded the coverage and increased the update speed of PolInSAR data, enabling the generation of large-scale, accurate DEM. High-precision, rapidly updated DEM plays a crucial role in fields such as environmental monitoring and disaster monitoring, providing precise data support for assessing forest coverage and terrain changes, among other issues. However, solving for the ground phase requires high-quality data, with the interferometric pair needing suitable temporal and spatial baselines, as well as good coherence. Low-quality data can result in inaccurate phase solutions, leading to errors in the inverted DEM. Furthermore, future research needs to further explore the applicability of this method under different climatic conditions or forest types and address potential challenges, such as considering the impact of varying climatic conditions on data quality and investigating interference effects and L-band penetration depth in different types of forests.

6. Conclusions

This paper proposes a MAPV inversion algorithm for RVoG model, and generates high-precision radar DEM data in forested area. The method is based on maximum a posteriori estimation and models the ground phase as a Von Mises distribution, with its mean derived from the topographic phase computed from external DEM. This method effectively overcomes the double-candidate effect, addresses the phase jump issue caused by modeling ground phases as a Gaussian distribution, and ensures the continuity and accuracy of ground phase solutions. Additionally, this paper derives, analyzes, and simulates the CRLB for MAPV. The results indicate that the introduction of prior information significantly reduces the CRLB of this method. In response to the characteristics of the MAPV objective function, this paper proposes a FSO method to enhance the efficiency of ground phase determination. FSO obtains two candidate phases separately based on the gradient descent algorithm and selects the one with the higher function value as the ground phase, as the gradient descent algorithm requires setting a global learning rate and there exist significant differences in the objective functions among different pixels. This paper improves the traditional RMSprop algorithm by reducing its dependency on the global learning rate, enabling it to handle various objective functions.
In this paper, the L-band spaceborne SAOCOM data are used to estimate the PolInSAR understory DEM to evaluate and verify the effectiveness of the method. The results demonstrate that:
  • Compared to the traditional exhaustive search method, FSO significantly improves computational efficiency without compromising accuracy.
  • The MAPV method effectively addresses the issue of elevation jumps in DEM caused by the discontinuity in ground phase solutions by MAPG.
  • Using IceSat-2 data as a benchmark, the DEM of the test forest area is compared with the DEMs of Alos, SRTM, and TSI. The results show that MAPV has better estimation performance, with improvements in both mean error (ME) and root mean square error (RMSE).

Author Contributions

Conceptualization, X.L. (Xiaoshuai Li); methodology, X.L. (Xiaoshuai Li); software, X.L. (Xiaoshuai Li); validation, X.L. (Xiaoshuai Li); formal analysis, X.L. (Xiaoshuai Li); investigation, X.L. (Xiaoshuai Li), X.L. (Xiaolei Lv) and Z.H.; resources, X.L. (Xiaoshuai Li), X.L. (Xiaolei Lv) and Z.H.; data curation, X.L. (Xiaoshuai Li); writing—original draft preparation, X.L. (Xiaoshuai Li); writing—review and editing, X.L. (Xiaoshuai Li), X.L. (Xiaolei Lv) and Z.H.; visualization, X.L. (Xiaoshuai Li); supervision, X.L. (Xiaolei Lv) and Z.H.; project administration, X.L. (Xiaolei Lv); funding acquisition, X.L. (Xiaolei Lv). All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the LuTan-1 L-Band Spaceborne Bistatic SAR data processing program, grant number E0H2080702.

Data Availability Statement

In this study, the remote sensing data were obtained from various sources to support our analyses. We accessed Satélite Argentino de Observación COn Microondas (SAOCOM) at http://saocom.asi.it:8081/#/home (accessed on 5 March 2024), Ice, Cloud and land Elevation Satellite-2 (ICESat-2) at https://nsidc.org/data/icesat-2 (accessed on 5 March 2024), Advanced Land Observing Satellite (ALOS) at https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm (accessed on 5 March 2024), Shuttle Radar Topography Mission 30 m (SRTM30) at https://earthexplorer.usgs.gov/ (accessed on 5 March 2024). These diverse data sources played a crucial role in our research and provided a comprehensive foundation for our remote sensing investigations.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DEMDigital Elevation Model
PolInSARPolarimetric Interferometric Synthetic Aperture Radar
LidarLight Detection and Ranging
RVoGRandom Volume over Ground
TSIThree-Stage Inversion
MAPMaximum a Posteriori
MAPGMaximum a Posteriori with Gaussian distribution as prior
MAPVMaximum a Posteriori with Von Mises distribution as prior
CRLBCramer–Rao Lower Bound
FSOFour-Step Optimization
MEMean Error
RMSERoot Mean Square Error

Appendix A. The Fisher Information of MAPV

It can be inferred from (14) that the Fisher information consists of two parts, with the first part provided by the RVoG model,
I F RVoG = E 2 log P ( R ^ ; R n , N ) ϕ 2 = E 2 ( log C ( R ^ ) N log | R n | N t r ( R n 1 R ^ ) ) ϕ 2 = E N | R n | | R n | ϕ N t r ( R n 1 R ^ ) ϕ = E N t r ( R n 1 R n ϕ ) N t r ( R n 1 R n ϕ R 1 R ^ ) = E N t r R n 1 ϕ R n ϕ + R n 1 2 R n ϕ 2 + N t r 2 R n 1 ϕ R n ϕ R n 1 + R n 1 2 R ϕ 2 R n 1 R ^ = E N t r R n 1 R n ϕ R n 1 R n ϕ + R n 1 2 R n ϕ 2 + N t r 2 R n 1 R n ϕ R n 1 R n ϕ R n 1 + R n 1 2 R n ϕ 2 R n 1 R ^ = N t r R n 1 R n ϕ R n 1 R n ϕ + N t r R n 1 2 R n ϕ 2 + 2 N t r R n 1 R n ϕ R n 1 R n ϕ N t r R n 1 2 R n ϕ 2 = N t r R n 1 R n ϕ R n 1 R n ϕ
and the second part being the information incorporated from the prior,
I F VM = E 2 log p ( ϕ ; ϕ topo , κ ) ϕ 2 = E 2 log 2 π I 0 ( κ ) + κ cos ( ϕ ϕ 0 ) ϕ 2 = E κ cos ( ϕ ϕ 0 ) = π π 1 2 π I 0 ( κ ) κ cos ( ϕ ϕ 0 ) e κ cos ( ϕ ϕ 0 ) d ϕ = κ I 1 ( κ ) I 0 ( κ ) .
Thus, the Fisher information for the ground phase ϕ can be expressed as
I F = I F RVoG + I F VM = N t r R n 1 R n ϕ R n 1 R n ϕ + κ I 1 ( κ ) I 0 ( κ ) .

Appendix B. Gradient of the MAPV Objective Function

From (17), we have
f ϕ = 3 sin θ 1 cos θ · θ ϕ 1 | A θ + ϕ | · | A θ + ϕ | ( θ + ϕ ) · ( θ + ϕ ) ϕ 1 | A ϕ | · | A ϕ | ϕ k sin ( ϕ ϕ 0 ) / N
where [46]
| A α | α = | A α | t r A α 1 A α α
and from (16),
A α α = j 2 e j α Ω ^ e j α Ω ^ H .
From (19),
θ ϕ = 2 arctan 3 | A ϕ | A ϕ + k N sin ( ϕ ϕ topo ) 1 ϕ = 2 1 + 9 | A ϕ | | A ϕ | + k N sin ( ϕ ϕ topo ) 2 · 3 | A ϕ | | A ϕ | + k N sin ( ϕ ϕ topo ) 1 ϕ = 2 | A ϕ | | A ϕ | + k N sin ( ϕ ϕ topo ) 2 | A ϕ | | A ϕ | + k N sin ( ϕ ϕ topo ) 2 + 9 · 3 | A ϕ | | A ϕ |     | A ϕ | 2 | A ϕ | 2 + k N cos ( ϕ ϕ topo ) | A ϕ | | A ϕ | + k N sin ( ϕ ϕ topo ) 2 = 6 | A ϕ | | A ϕ |     | A ϕ | 2 | A ϕ | 2 + k N cos ( ϕ ϕ topo ) | A ϕ | | A ϕ | + k N sin ( ϕ ϕ topo ) 2 + 9
where | A ϕ | represents the second derivative of | A ϕ | with respect to ϕ ,
| A ϕ | =   | A ϕ | t r A ϕ 1 A ϕ ϕ + | A ϕ | t r A ϕ 1 ϕ A ϕ ϕ + A ϕ 1 2 A ϕ ϕ 2 =   | A ϕ | t r 2 A ϕ 1 A ϕ ϕ + | A ϕ | t r A ϕ 1 A ϕ ϕ A ϕ 1 A ϕ ϕ + A ϕ 1 2 A ϕ ϕ 2
and
2 A ϕ ϕ 2 = 1 2 e j ϕ Ω ^ + e j ϕ Ω ^ H .
Substituting (A8) into (A7), we can obtain that
θ ϕ = 6 t r A ϕ 1 A ϕ ϕ A ϕ 1 A ϕ ϕ + A ϕ 1 2 A ϕ ϕ 2 + k N cos ( ϕ ϕ topo ) t r A ϕ 1 A ϕ ϕ + k N sin ( ϕ ϕ topo ) 2 + 9 .
Therefore,
f ϕ = 3 sin θ 1 cos θ θ ϕ t r A θ + ϕ 1 A θ + ϕ ( θ + ϕ ) ( θ + ϕ ) ϕ t r A ϕ 1 A ϕ ϕ k N sin ( ϕ ϕ 0 )
the gradient of the MAPV objective function can be obtained by substituting (A6), (A9), and (A10) into (A11).

References

  1. FAO. Global Forest Resources Assessment 2020; FAO: Rome, Italy, 2020. [Google Scholar] [CrossRef]
  2. Li, F.; Zhao, Y. Development of a GIS-based decision-support system of forest resource management. Sci. China Ser. E Technol. Sci. 2006, 49, 76–85. [Google Scholar] [CrossRef]
  3. Dias, H.C.; Hölbling, D.; Dias, V.C.; Grohmann, C.H. Application of Object-Based Image Analysis for Detecting and Differentiating between Shallow Landslides and Debris Flows. GI_Forum 2023 2023, 11, 34–44. [Google Scholar] [CrossRef]
  4. Akay, A.E.; Oğuz, H.; Karas, I.R.; Aruga, K. Using LiDAR Technology in Forestry Activities. Environ. Monit. Assess. 2009, 151, 117–125. [Google Scholar] [CrossRef] [PubMed]
  5. Luo, W.; Ma, H.; Yuan, J.; Zhang, L.; Ma, H.; Cai, Z.; Zhou, W. High-Accuracy Filtering of Forest Scenes Based on Full-Waveform LiDAR Data and Hyperspectral Images. Remote Sens. 2023, 15, 3499. [Google Scholar] [CrossRef]
  6. Cheng, L.; Hao, R.; Cheng, Z.; Li, T.; Wang, T.; Lu, W.; Ding, Y.; Hu, H. Modeling the Global Relationship via the Point Cloud Transformer for the Terrain Filtering of Airborne LiDAR Data. Remote Sens. 2023, 15, 5434. [Google Scholar] [CrossRef]
  7. Pourshamsi, M.; Xia, J.; Yokoya, N.; Garcia, M.; Lavalle, M.; Pottier, E.; Balzter, H. Tropical Forest Canopy Height Estimation from Combined Polarimetric SAR and LiDAR Using Machine-Learning. ISPRS J. Photogramm. Remote Sens. 2021, 172, 79–94. [Google Scholar] [CrossRef]
  8. Lopez-Martinez, C.; Alonso, A.; Fabregas, X.; Papathannassiou, K.P. Ground Topography Estimation over Forests Considering Polarimetric SAR Interferometry. In Proceedings of the 2010 IEEE International Geoscience and Remote Sensing Symposium, Honolulu, HI, USA, 25–30 July 2010; pp. 3612–3615. [Google Scholar] [CrossRef]
  9. Fu, H.; Wang, C.; Zhu, J.; Xie, Q.; Zhang, B. Estimation of Pine Forest Height and Underlying DEM Using Multi-Baseline P-Band PolInSAR Data. Remote Sens. 2016, 8, 820. [Google Scholar] [CrossRef]
  10. Fu, H.; Zhu, J.; Wang, C.; Wang, H.; Zhao, R. Underlying Topography Estimation over Forest Areas Using High-Resolution P-Band Single-Baseline PolInSAR Data. Remote Sens. 2017, 9, 363. [Google Scholar] [CrossRef]
  11. Fu, H.; Zhu, J.; Wang, C.; Li, Z. Underlying Topography Extraction over Forest Areas from Multi-Baseline PolInSAR Data. J. Geod. 2018, 92, 727–741. [Google Scholar] [CrossRef]
  12. Lu, H.; Suo, Z.; Guo, R.; Bao, Z. S-RVoG Model for Forest Parameters Inversion over Underlying Topography. Electron. Lett. 2013, 49, 618–620. [Google Scholar] [CrossRef]
  13. Shenglong Guo, S.G.; Yurun Tian, Y.T.; Yang Li, Y.L.; Qiang Yin, Q.Y.; Wen Hong, W.H. Estimation of Ground Topography under Forests with Polarimetric SAR Interferometry Data. In Proceedings of the IET International Radar Conference 2015, Hangzhou, China, 14–16 October 2015; p. 6. [Google Scholar] [CrossRef]
  14. Mette, T.; Papathanassiou, K.; Hajnsek, I. Biomass Estimation from Polarimetric SAR Interferometry over Heterogeneous Forest Terrain. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, 2004—IGARSS ’04: Proceedings, 2004, Anchorage, AK, USA, 20–24 September 2004; Volume 1, pp. 511–514. [Google Scholar] [CrossRef]
  15. Yamada, H.; Yamaguchi, Y.; Rodriguez, E.; Kim, Y.; Boerner, W. Polarimetric SAR Interferometry for Forest Canopy Analysis by Using the Super-Resolution Method. In Proceedings of the IGARSS 2001: Scanning the Present and Resolving the Future: Proceedings: IEEE 2001 International Geoscience and Remote Sensing Symposium (Cat. No.01CH37217), Sydney, NSW, Australia, 9–13 July 2001; Volume 3, pp. 1101–1103. [Google Scholar] [CrossRef]
  16. Duan, D. A Modified ESPRIT Algorithm to Estimate Tree Height Using Simulated Dual-Polarization PolInSAR Datasets. In Proceedings of the EUSAR 2018; 12th European Conference on Synthetic Aperture Radar, Aachen, Germany, 4–7 June 2018. [Google Scholar]
  17. Lei, Y.; Jun, Z.Y.; Gang, W.Z. Joint Phase and Power Estimation for Polarimetric Interferometric SAR Based on TEL-ESPRIT Algorithm. In Proceedings of the 2006 CIE International Conference on Radar, Shanghai, China, 16–19 October 2006; pp. 1–4. [Google Scholar] [CrossRef]
  18. Yamada, H.; Sato, K.; Yamaguchi, Y.; Boerner, W.M. Interferometric Phase and Coherence of Forest Estimated by ESPRIT-based Polarimetric SAR Interferometry. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, Toronto, ON, Canada, 24–28 June 2002; Volume 2, pp. 829–831. [Google Scholar] [CrossRef]
  19. Yamada, H.; Yamaguchi, Y.; Boerner, W. Forest Height Feature Extraction in Polarimetric SAR Interferometry by Using Rotational Invariance Property. In Proceedings of the IGARSS 2003: 2003 IEEE International Geoscience and Remote Sensing Symposium: Proceedings (IEEE Cat. No.03CH37477), Toulouse, France, 21–25 July 2003; Volume 3, pp. 1426–1428. [Google Scholar] [CrossRef]
  20. Minh, N.P.; Wang, C.; Zou, B.; Nguyen, Q.T.; Le, V.N. Forest Height Extraction from PolInSAR Image Using a Hybrid Method. Int. J. Signal Process. Image Process. Pattern Recognit. 2014, 7, 257–274. [Google Scholar] [CrossRef]
  21. Guillaso, S.; Reigber, A.; Ferro-Famil, L. Evaluation of the ESPRIT Approach in Polarimetric Interferometric SAR. In Proceedings of the 2005 IEEE International Geoscience and Remote Sensing Symposium, 2005. IGARSS ’05, Seoul, Republic of Korea, 29 July 2005; Volume 1, pp. 32–35. [Google Scholar] [CrossRef]
  22. Treuhaft, R.N.; Madsen, S.N.; Moghaddam, M.; Van Zyl, J.J. Vegetation Characteristics and Underlying Topography from Interferometric Radar. Radio Sci. 1996, 31, 1449–1485. [Google Scholar] [CrossRef]
  23. Treuhaft, R.N.; Siqueira, P.R. Vertical Structure of Vegetated Land Surfaces from Interferometric and Polarimetric Radar. Radio Sci. 2000, 35, 141–177. [Google Scholar] [CrossRef]
  24. Cloude, S.; Papathanassiou, K. Three-Stage Inversion Process for Polarimetric SAR Interferometry. IEE Proc.—Radar Sonar Navig. 2003, 150, 125. [Google Scholar] [CrossRef]
  25. Tabb, M.; Flynn, T.; Carande, R. Full Maximum Likelihood Inversion of Polinsar Scattering Models. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, 2004, IGARSS ’04. Proceedings, Anchorage, AK, USA, 20–24 September 2004; Volume 2, pp. 1232–1235. [Google Scholar] [CrossRef]
  26. Tabb, M.; Carande, R. Robust Inversion of Vegetation Structure Parameters from Low-Frequency, Polarimetric Interferometric SAR. In Proceedings of the IGARSS 2001: Scanning the Present and Resolving the Future: Proceedings: IEEE 2001 International Geoscience and Remote Sensing Symposium (Cat. No.01CH37217), Sydney, NSW, Australia, 9–13 July 2001; Volume 7, pp. 3188–3190. [Google Scholar] [CrossRef]
  27. Kugler, F.; Lee, S.-K.; Hajnsek, I.; Papathanassiou, K.P. Forest Height Estimation by Means of Pol-InSAR Data Inversion: The Role of the Vertical Wavenumber. IEEE Trans. Geosci. Remote Sens. 2015, 53, 5294–5311. [Google Scholar] [CrossRef]
  28. Huang, Z.; Lv, X.; Li, X.; Chai, H. Maximum a Posteriori Inversion for Forest Height Estimation Using Spaceborne Polarimetric SAR Interferometry. IEEE Trans. Geosci. Remote Sens. 2023, 61, 1–14. [Google Scholar] [CrossRef]
  29. Fisher, N.I. Statistical Analysis of Circular Data; Cambridge University Press: Cambridge, UK; New York, NY, USA, 1993. [Google Scholar]
  30. Letzepis, N. On the von Mises Approximation for the Distribution of the Phase Angle between Two Independent Complex Gaussian Vectors. In Proceedings of the 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), South Brisbane, QLD, Australia, 19–24 April 2015; pp. 3247–3251. [Google Scholar] [CrossRef]
  31. Jammalamadaka, S.R.; Sengupta, A. Topics in Circular Statistics; Number v. 5 in Series on Multivariate Analysis; World Scientific: River Edge, NJ, USA, 2001. [Google Scholar]
  32. Haji, S.H.; Abdulazeez, A.M. Comparison of optimization techniques based on gradient descent algorithm: A review. PalArch’s J. Archaeol. Egypt/Egyptol. 2021, 18, 2715–2743. [Google Scholar]
  33. Ruder, S. An Overview of Gradient Descent Optimization Algorithms. arXiv 2016, arXiv:1609.04747. [Google Scholar]
  34. Abdalati, W.; Zwally, H.J.; Bindschadler, R.; Csatho, B.; Farrell, S.L.; Fricker, H.A.; Harding, D.; Kwok, R.; Lefsky, M.; Markus, T.; et al. The ICESat-2 Laser Altimetry Mission. Proc. IEEE 2010, 98, 735–751. [Google Scholar] [CrossRef]
  35. Tian, X.; Shan, J. Comprehensive Evaluation of the ICESat-2 ATL08 Terrain Product. IEEE Trans. Geosci. Remote Sens. 2021, 59, 8195–8209. [Google Scholar] [CrossRef]
  36. Cloude, S.; Papathanassiou, K. Polarimetric SAR Interferometry. IEEE Trans. Geosci. Remote Sens. 1998, 36, 1551–1565. [Google Scholar] [CrossRef]
  37. Cloude, S.R. Polarization Coherence Tomography: Polarization coherence tomography. Radio Sci. 2006, 41, 1–27. [Google Scholar] [CrossRef]
  38. Ferretti, A.; Monti-Guarnieri, A.; Prati, C.; Rocca, F.; Massonet, D. InSAR Principles—Guidelines for SAR Interferometry Processing and Interpretation. ESA Train. Man. 2007, 19, A17–A19. [Google Scholar]
  39. Lee, A. Circular data. Wiley Interdiscip. Rev. Comput. Stat. 2010, 2, 477–486. [Google Scholar] [CrossRef]
  40. Neumann, M.; Ferro-Famil, L.; Reigber, A. Estimation of Forest Structure, Ground, and Canopy Layer Characteristics From Multibaseline Polarimetric Interferometric SAR Data. IEEE Trans. Geosci. Remote Sens. 2010, 48, 1086–1104. [Google Scholar] [CrossRef]
  41. Wang, X.; Xu, F. A PolinSAR Inversion Error Model on Polarimetric System Parameters for Forest Height Mapping. IEEE Trans. Geosci. Remote Sens. 2019, 57, 5669–5685. [Google Scholar] [CrossRef]
  42. Roueff, A.; Arnaubec, A.; Dubois-Fernandez, P.C.; Refregier, P. Cramer–Rao Lower Bound Analysis of Vegetation Height Estimation With Random Volume Over Ground Model and Polarimetric SAR Interferometry. IEEE Geosci. Remote Sens. Lett. 2011, 8, 1115–1119. [Google Scholar] [CrossRef]
  43. Arnaubec, A.; Roueff, A.; Dubois-Fernandez, P.C.; Refregier, P. Vegetation Height Estimation Precision With Compact PolInSAR and Homogeneous Random Volume Over Ground Model. IEEE Trans. Geosci. Remote Sens. 2014, 52, 1879–1891. [Google Scholar] [CrossRef]
  44. Refregier, P.; Roueff, A.; Arnaubec, A.; Dubois-Fernandez, P.C. Invariant Contrast Parameters of PolInSAR Homogenous RVoG Model. IEEE Geosci. Remote Sens. Lett. 2014, 11, 1414–1417. [Google Scholar] [CrossRef]
  45. Seymour, M.; Cumming, I. Maximum Likelihood Estimation for SAR Interferometry. In Proceedings of the IGARSS ’94—1994 IEEE International Geoscience and Remote Sensing Symposium, Pasadena, CA, USA, 8–12 August 1994; Volume 4, pp. 2272–2275. [Google Scholar] [CrossRef]
  46. Petersen, K.B.; Pedersen, M.S. The Matrix Cookbook; Technical University of Denmark: Lyngby, Denmark, 2008; Volume 7, pp. 7–9. [Google Scholar]
Figure 1. (a) Optical image of the test area in Google Earth. (b) L-band SAR image in the Pauli basis (R: HH - VV; G: HV; B: HH + VV) (c) Alos DEM projected to SAR coordinates. (d) The topographic phase estimated by Alos DEM.
Figure 1. (a) Optical image of the test area in Google Earth. (b) L-band SAR image in the Pauli basis (R: HH - VV; G: HV; B: HH + VV) (c) Alos DEM projected to SAR coordinates. (d) The topographic phase estimated by Alos DEM.
Remotesensing 16 00948 g001
Figure 2. The DEM inversion flowchart for spaceborne PolInSAR. (a) Ground phase estimation based on MAPV. (b) DEM generation.
Figure 2. The DEM inversion flowchart for spaceborne PolInSAR. (a) Ground phase estimation based on MAPV. (b) DEM generation.
Remotesensing 16 00948 g002
Figure 3. The difference between Gauss distribution and Von Mises distribution. Up: In the Cartesian coordinate system. (a) Gauss distribution. (b) Von Mises distribution. Down: In polar coordinate. (c) Gauss distribution. (d) Von Mises distribution. The hollow black point represents the position of the topographic phase on the unit circle of the complex plane, and the asterisks correspond to the phases at the two peaks of the MLE.
Figure 3. The difference between Gauss distribution and Von Mises distribution. Up: In the Cartesian coordinate system. (a) Gauss distribution. (b) Von Mises distribution. Down: In polar coordinate. (c) Gauss distribution. (d) Von Mises distribution. The hollow black point represents the position of the topographic phase on the unit circle of the complex plane, and the asterisks correspond to the phases at the two peaks of the MLE.
Remotesensing 16 00948 g003
Figure 4. Illustration of the difference between MAPG and MAPV. The solid line, dotted line, and dashed line represent the normalized objective functions of the MLE method, MAPG method, and MAPV method, respectively.
Figure 4. Illustration of the difference between MAPG and MAPV. The solid line, dotted line, and dashed line represent the normalized objective functions of the MLE method, MAPG method, and MAPV method, respectively.
Remotesensing 16 00948 g004
Figure 5. (a) Evolution of the CRLB with signal-to-noise ratio (SNR). (b) A detail extracted from (a).
Figure 5. (a) Evolution of the CRLB with signal-to-noise ratio (SNR). (b) A detail extracted from (a).
Remotesensing 16 00948 g005
Figure 6. Illustration of the four-step optimization procedure for solving the ground phase.
Figure 6. Illustration of the four-step optimization procedure for solving the ground phase.
Remotesensing 16 00948 g006
Figure 7. The convergence of the improved RMSprop. The top (ac) and bottom (df) rows, respectively, depict the convergence of phase and objective function values when the global learning rate ϵ = 0.1, 0.4, and 0.7.
Figure 7. The convergence of the improved RMSprop. The top (ac) and bottom (df) rows, respectively, depict the convergence of phase and objective function values when the global learning rate ϵ = 0.1, 0.4, and 0.7.
Remotesensing 16 00948 g007
Figure 8. Ground phase estimated by (a) MAPG. (b) MAPV.
Figure 8. Ground phase estimated by (a) MAPG. (b) MAPV.
Remotesensing 16 00948 g008
Figure 9. (a) MAPG-inverted DEM. (b) MAPV-inverted DEM.
Figure 9. (a) MAPG-inverted DEM. (b) MAPV-inverted DEM.
Remotesensing 16 00948 g009
Figure 10. (a) Section elevation change curve of marked line segment in Area 1. (b) Section elevation change curve of marked line segment in Area 2.
Figure 10. (a) Section elevation change curve of marked line segment in Area 1. (b) Section elevation change curve of marked line segment in Area 2.
Remotesensing 16 00948 g010
Figure 11. The geocoded DEM maps and the difference between the LiDAR DEM and the DEM from different sources. (a) Alos-30m. (b) SRTM-30m (c) TSI. (d) MAPV.
Figure 11. The geocoded DEM maps and the difference between the LiDAR DEM and the DEM from different sources. (a) Alos-30m. (b) SRTM-30m (c) TSI. (d) MAPV.
Remotesensing 16 00948 g011
Figure 12. Histogram of the difference between the DEM obtained through different methods and the LiDAR DEM.
Figure 12. Histogram of the difference between the DEM obtained through different methods and the LiDAR DEM.
Remotesensing 16 00948 g012
Table 1. Efficiency assessment of four-step optimization method.
Table 1. Efficiency assessment of four-step optimization method.
SizeTime (s)Iterations α 1 α 2
ESFSOESFSO
7015 × 26736880.01235.136024.599.97%99.99%
Table 2. Differential comparison of DEM continuity in inversion methods.
Table 2. Differential comparison of DEM continuity in inversion methods.
ID α 15
MAPGMAPV
Area 187.48%99.06%
Area 289.13%99.14%
Table 3. Assessment of the elevation accuracy of DEM.
Table 3. Assessment of the elevation accuracy of DEM.
DEMME (m)RMSE (m)
TSI−13.240773.7140
ALOS2.63707.7956
SRTM2.13117.8972
MAPV0.21115.9944
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Li, X.; Lv, X.; Huang, Z. Underlying Topography Estimation over Forest Using Maximum a Posteriori Inversion with Spaceborne Polarimetric SAR Interferometry. Remote Sens. 2024, 16, 948. https://0-doi-org.brum.beds.ac.uk/10.3390/rs16060948

AMA Style

Li X, Lv X, Huang Z. Underlying Topography Estimation over Forest Using Maximum a Posteriori Inversion with Spaceborne Polarimetric SAR Interferometry. Remote Sensing. 2024; 16(6):948. https://0-doi-org.brum.beds.ac.uk/10.3390/rs16060948

Chicago/Turabian Style

Li, Xiaoshuai, Xiaolei Lv, and Zenghui Huang. 2024. "Underlying Topography Estimation over Forest Using Maximum a Posteriori Inversion with Spaceborne Polarimetric SAR Interferometry" Remote Sensing 16, no. 6: 948. https://0-doi-org.brum.beds.ac.uk/10.3390/rs16060948

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop