An Efficient Maximum Likelihood Estimation Approach of Multi-Baseline SAR Interferometry for Refined Topographic Mapping in Mountainous Areas

For InSAR topographic mapping, multi-baseline InSAR height estimation is known to be an effective way to facilitate phase unwrapping by significantly increasing the ambiguity intervals and maintaining good height measurement sensitivity, especially in mountainous areas. In this paper, an efficient multi-baseline SAR interferometry approach based on maximum likelihood estimation is developed for refined topographic mapping in mountainous areas. In the algorithm, maximum likelihood (ML) height estimation is used to measure the topographic details and avoid the complicated phase unwrapping process. In order to be well-adapted to the mountainous terrain conditions, the prior height probability is re-defined to take the local terrain conditions and neighboring height constraint into consideration in the algorithm. In addition, three strategies are used to optimize the maximum likelihood height estimation process to obtain higher computational efficiency, so that this method is more suitable for spaceborne InSAR data. The strategies include substituting a rational function model into the complicated conversion process from candidate height to interferometric phase, discretizing the continuous height likelihood probability, and searching for the maximum likelihood height with a flexible step length. The experiment with simulated data is designed to verify the improvement of the ML height estimation accuracy with the re-defined prior height distribution. Then the optimized processing procedure is tested with the multi-baseline L-band ALOS/PALSAR data covering the Mount Tai area in China. The height accuracy of the generated multi-baseline InSAR DEM can meet both standards of American DTED-2 and Chinese national 1:50,000 DEM (mountain) Level 2.


Introduction
SAR interferometry (InSAR) is an effective tool for large-area topographic mapping due to its all-weather imaging and high sensitivity to terrain relief [1,2].The InSAR height measurement accuracy is greatly influenced by the phase unwrapping accuracy and the length of normal baselines [3,4].Longer normal baselines allow more accurate height estimation but also generate higher frequency of the interferometric fringes, which increases the complexity of phase unwrapping.On the other hand, shorter normal baselines reduce the complexity of phase unwrapping but suffer from poorer phase-to-height sensitivity [5,6].Therefore, the contradiction between the sensitivity of height measurement and the reliability of phase unwrapping caused by the length of normal baseline over rough terrain is inevitable for single-baseline InSAR.
To combine the advantages of large and short baselines in topographic mapping, a multi-baseline InSAR principle has been proposed to estimate terrain height by joint analysis of multiple interferometric pairs with diverse normal baselines [7,8].Researchers have proposed a variety of methods for multi-baseline InSAR processing.The standard method to increase the accuracy of interferometric DEMs is to stack multiple geocoded layers, weighted with the individual error estimates [9][10][11].However, all DEM stacking approaches assume correctly unwrapped interferograms.Moreover, it is difficult to co-register different geocoded DEM layers accurately and the horizontal displacements can introduce height errors.Other multi-baseline estimation methods have been published to facilitate the phase unwrapping process by taking advantage of baseline diversity, such as the Least Square estimation method [5], the iterative multi-baseline method [12], the Chinese Remainder Theorem (CRT) method [13,14], and so on.These multi-baseline phase unwrapping methods can significantly increase the ambiguity intervals of interferometric phases and keep the topographic details as well; however, these methods still have to solve the phase unwrapping problem correctly.
In order to avoid the phase unwrapping process and determine the target height directly, the statistical method using the criterion of maximum likelihood [9,[15][16][17][18] is exploited to combine the multi-baseline information for target height estimation.However, in actual data processing, atmospheric effects, orbital errors, and decorrelation will introduce phase noise that cannot be ignored.Therefore, the maximum likelihood estimation method is not robust enough to search for the target height within all elevation values.In order to realize more robust and reliable height estimation, the prior height information from the reference DEM is incorporated into the ML estimation to restrict the height searching range [18], but the problem is that the local terrain conditions and the neighboring height constraint are not considered in the prior height distribution.The maximum a posteriori (MAP) estimation tries to introduce the neighboring height constraints by using Markov random fields to model the prior distribution of the unknown images [19,20].This method allows recovering topographic profiles affected by strong height discontinuities and performing efficient noise rejections.However, the MAP methods also have limits concerning the computational time and the optimization step because there is no guarantee of finding the global optimum.
In this article, we apply the maximum likelihood estimation for multi-baseline InSAR DEM generation.The prior height probability is re-defined to take the local terrain conditions and neighboring height constraint into consideration.An experiment with simulated data is designed to verify the improvement of the ML height estimation accuracy with the well-defined prior height distribution.Furthermore, to make the maximum likelihood estimation method well adapted to spaceborne SAR data, the processing flow is optimized for higher computational efficiency with the following innovative points: (1) Replacing the rigorous height-to-phase conversion with the rational function model (RFM); (2) substituting the complicated height likelihood probability function with two-dimension lookup table ; (3) searching for the maximum likelihood height with flexible search step length instead of fixed search step length.This processing flow is testified with the L-band ALOS/PALSAR data, which can be less influenced by the temporal and volume decorrelation and have a longer critical baseline than InSAR data with a shorter wavelength (such as X band or C band) [21].Since the ALOS/PALSAR data were acquired in the repeat-pass mode, the above processing flow integrates atmospheric effect correction to improve the reliability of multi-baseline estimation.The rest of this article is structured as follows: Section 2 introduces the principle of ML estimation with prior DEM; Section 3 presents the improved proposed processing flow for spaceborne datasets; detailed descriptions of the experiments and results are given in Section 4; Section 5 is a discussion of the experimental results; finally, the conclusions are drawn.

Basic Principle
The maximum likelihood (ML) height estimation method of multi-baseline InSAR considers the target height h as a parameter of the probability distribution of the interferometric phase φ, denoted as pd f (φ|h) , and combines the probability distributions of all the interferometric phase observations to estimate the target height with the maximum likelihood criterion [22].Aiming at reducing the undesired variation in the maximum likelihood height estimation caused by non-negligible phase noise, the prior height distribution from a reference DEM is integrated into the maximum likelihood estimation, which not only greatly improves the estimation reliability but also narrows the search range and increases the computational efficiency.Therefore, the maximum likelihood height estimation assisted by the reference DEM is shown in Equation ( 1): where pd f (h) is the a priori distribution function of height provided by the reference DEM.pd f (φ i |h) is the height likelihood function for the ith interferogram.According to [23], pd f (φ|h) can be estimated as the edge probability density function (PDF) of the interferometric phase φ, as shown in Equation ( 2): where β =|ρ|cos(φ − φ 0 ) ; φ 0 is the mathematical expectation of interferometric phase; ρ is the complex coherence coefficient; and L is the effective number of looks (ENL).φ 0 can be represented by the target height h through height-to-phase conversion function F h2φ (h).When φ 0 is set as zero, pd f (φ|0, ρ, L) describes the probability distribution of the interferometric phase noise.Hence the standard deviation of phase noise σ φ can be derived by Equation (3).Given L, σ φ is determined by the coherence coefficient ρ: The standard deviation of the height errors σ h can be approximated as: where λ is the wavelength; R is the slant range; θ is the incidence angle; and B ⊥ is the normal baseline.

Definition of the Prior Height Probability
Suppose the height acquired from the prior DEM corresponding to a resolution unit of the interferogram is h prior .Generally, we think that the system error of the prior DEM has been corrected.Then the height error is the accidental error and obeys Gaussian distribution.The standard deviation of the height errors for each cell in the interferogram is σ h .Hence, the prior height probability distribution for each cell can be defined as in Equation ( 5) [18]: In this article, the SRTM DEM is used as the reference DEM with standard deviation of the height errors σ SRTM .The σ SRTM of the SRTM DEM in the Eurasian continent is 3.8 m [24].For each cell of the interferogram, if we suppose that σ h = σ SRTM directly, Equation (5) defines the prior height probability distribution over all the geographic coverage of SRTM DEM, which can be much larger than the SAR image coverage, causing the following problems: (1) the local terrain conditions such as plains or mountains, are not taken into account, while σ h under different terrain conditions is not consistent; (2) the neighboring height constraints are not considered that the height probability distribution defined by Equation ( 5) is only related to the height of the resolution unit itself.Aimed at these two problems, Equation ( 5) is modified to take the local terrain conditions in the neighborhood into consideration.The prior height probability distribution is re-defined as in Equation ( 6): where N is composed of the resolution units and its adjacent units, which are usually four-neighbor, eight-neighbor, or 24-neighbor.h SRTM,i , i ∈ N represent the heights of the resolution units, T represents the number of the pixels in the neighborhood, and σ local is the standard deviation of height errors in the neighborhood.Both h SRTM and σ local are obtained from the SRTM DEM.It can be seen from Equation ( 6) that when σ local is larger than σ SRTM , σ h is set to σ local .Under this condition, σ h is no longer a constant.In the undulating terrain areas, σ h depends on the local terrain conditions, while in the flat areas it is still conservatively set as σ SRTM .
The size of neighborhood used to define the prior height distribution probability is determined based on the spatial resolution of the prior DEM and the coherence level of the interferogram.When the spatial resolution of the prior DEM is much lower than that of the interferogram and the coherence level of the interferogram is high, a smaller neighborhood such as four-neighbor is preferred to reduce the influence of the prior height and neighboring heights in the height probability distribution so that the resulting DEM will not be too smooth.Otherwise, when the spatial resolution of the prior DEM is close to that of the interferogram or the coherence level of the interferogram is not high, a larger neighborhood such as eight-neighbor or even 24-neighbor is selected to enhance the constraints of the prior height and neighboring heights on the height probability distribution so that the impact of the phase noise is suppressed and the height estimate is more robust.
In this article, SRTM DEM with a cell size of 30 m is used as the prior DEM and the interferogram has a comparable cell size, about 22 m after 3 × 7 multi-look processing.The mean coherence of the interferogram is about 0.5, which is a moderate coherence level.Therefore, eight-neighbor is chosen in the experiment.
The correspondence between SRTM DEM and the resolution cell of the interferogram needs to be established by radar coding.The height error introduced by the radar coding procedure will increase the value of σ SRTM , hence, in practical calculations, σ SRTM will be adjusted empirically to make the height distribution curve more reasonable.

Optimized Processing Flow for Spaceborne Multi-Baseline InSAR Datasets
When applying the multi-baseline InSAR height estimation with maximum likelihood criterion in spaceborne InSAR datasets, the processing flow can be divided into three major stages.First is the interferometric processing including interferometric pairs combination and differential interferogram generation, as in Section 3.1.Second is the maximum likelihood height estimation process based on the principle introduced in Section 2. In order to make the ML estimation method well adapted to the spaceborne data, the processing flow is optimized for higher computational efficiency with the following new points: (1) Replacing the rigorous height-to-phase conversion with the rational function model (RFM), as in Section 3.2.1;(2) substituting the complicated height likelihood probability function with two-dimension lookup table, as in Section 3.2.2;(3) searching for the maximum likelihood height with flexible search step length instead of the fixed search step length, as in Section 3.2.3.Thirdly, the estimated height map in SAR image coordinate can be geocoded into the geographic coordinate system or the universal transverse mercator (UTM) system, as in Section 3.3.The flowchart of this approach is outlined in Figure 1.A brief description and rationale for each step are given as follows.
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 20 interferogram generation, as in Section 3.1.Second is the maximum likelihood height estimation process based on the principle introduced in Section 2. In order to make the ML estimation method well adapted to the spaceborne data, the processing flow is optimized for higher computational efficiency with the following new points: (1) Replacing the rigorous height-to-phase conversion with the rational function model (RFM), as in Section 3.2.1;(2) substituting the complicated height likelihood probability function with two-dimension lookup table, as in Section 3.2.2;(3) searching for the maximum likelihood height with flexible search step length instead of the fixed search step length, as in Section 3.2.3.Thirdly, the estimated height map in SAR image coordinate can be geocoded into the geographic coordinate system or the universal transverse mercator (UTM) system, as in Section 3.3.The flowchart of this approach is outlined in Figure 1.A brief description and rationale for each step are given as follows.

Interferometric Processing
For the repeat-pass interferometry, we first need to select the suitable SAR images in the given dataset to constitute interferometric pairs and then choose the master interferogram for other interferograms to register with, while for single-pass interferometry we merely need to choose the master interferogram.To keep good coherence, the interferometric pairs should have temporal baselines as short as possible under the premise that the normal baseline is shorter than the critical normal baseline.In order to select the proper master interferogram, the principle is to select a master interferogram at around the center of the time axis among all the interferograms to make it easier for the other interferograms to register to it.
Next is the interferometric processing step for each interferometric pair, which is the basic processing unit for the maximum likelihood height estimation, including complex image coregistration, interferogram generation, flattening, and filtering.The SRTM DEM is projected into the azimuth and slant-range coordinate system of the master SAR image in the interferometric pair.Then the interferometric phase is flattened by the radar-coded DEM, i.e., major 2 phase jumps due to π

Interferometric Processing
For the repeat-pass interferometry, we first need to select the suitable SAR images in the given dataset to constitute interferometric pairs and then choose the master interferogram for other interferograms to register with, while for single-pass interferometry we merely need to choose the master interferogram.To keep good coherence, the interferometric pairs should have temporal baselines as short as possible under the premise that the normal baseline is shorter than the critical normal baseline.In order to select the proper master interferogram, the principle is to select a master interferogram at around the center of the time axis among all the interferograms to make it easier for the other interferograms to register to it.
Next is the interferometric processing step for each interferometric pair, which is the basic processing unit for the maximum likelihood height estimation, including complex image co-registration, interferogram generation, flattening, and filtering.The SRTM DEM is projected into the azimuth and slant-range coordinate system of the master SAR image in the interferometric pair.Then the interferometric phase is flattened by the radar-coded DEM, i.e., major 2π phase jumps due to topography are removed.For repeat-pass interferometry, the atmospheric phase screen (APS) introduces non-negligible height error into the InSAR height [25].The APS consists of a vertically stratified component and a turbulent mixing one [6].Based on spatial pattern analyses of these two Remote Sens. 2018, 10, 454 6 of 20 APS components, a SRTM elevation-to-phase regression model and a low-pass plus adaptive combined filter are employed to estimate and remove them from the differential interferogram sequentially [25].
After interferometric processing, the differential interferograms are registered to the image space of the chosen master interferogram.The phase shift between the interferograms is a constant.In the repeat-pass interferometry mode, the phase shift can be removed through the combined filter used for the APS correction [25], while in the single-pass interferometry mode it can be obtained through statistical average value of the differential phase maps.

Rational Function Model (RFM) for Height-to-Phase Conversion
The candidate height must be converted to the interferometric phase to calculate the corresponding height likelihood probability as in Equation ( 2).This height-to-phase conversion starts from a pixel in the master image with the image coordinates (i M , j M ), and its slant range R M can be determined with the orbital information.The candidate height for pixel (i M , j M ) is represented by h DEM .The Range-Doppler (RD) model can be iteratively solved to calculate the geographic coordinates for pixel (i M , j M ), which is the direct positioning process.With the calculated geographic coordinates and orbital information of the slave image, the RD model can be iteratively solved again to calculate its image coordinates (i S , j S ) in the slave image, which is the indirect positioning process.The slant range for the pixel (i S , j S ) can be determined for the slave image as R S .The corresponding interferometric phase for the chosen pixel with the candidate height h DEM can be determined through Equation ( 7), where λ is the wavelength.For repeat-pass interferometry, k = 2 and for single-pass interferometry, k = 1.
As we can see, this height-to-phase conversion process requires iteratively solving the RD model twice for every candidate height in the height search range; this process has to be performed for every pixel in the SAR image, which is extremely time-consuming.
In order to improve the efficiency of height-to-phase conversion, we try to no longer care about the specific analytical form of height-to-phase conversion functions, but rather write them directly as a rational function of the image coordinates (L_a,P_r), the height h, as shown in Equation ( 8): They are called the rational function model of the height-to-phase conversion function F h2φ (h). a i,j,n , b i,j,n are the unknown parameters to be solved.The value of b 0,0,0 is set to 1.The way to solve the unknown parameters is the same as in [26].
The rational function model for the height-to-phase conversion has been established in this article with the spaceborne InSAR data and tested to see whether it can replace the rigorous method.The number of the control points is 50 × 50 × 10, that is, there are 50 × 50 regular grid points in the plane and 10 layers in the height range.The height interval is between −500 m and 10,000 m.We randomly generate 10,000 checkpoints at which the rational function model and the rigorous method are used for the height-to-phase conversion, respectively.The conversion error of the rational function model is calculated using the results of the rigorous method as a reference.Table 1 lists the conversion error for different spaceborne SAR data, indicating that the conversion error for the rational function model is completely negligible in practical applications.Therefore, the rational function model can replace the rigorous method for height-to-phase conversion at each resolution unit.2) is used to calculate the height likelihood probability, which is a very complicated expression.In order to improve the computational efficiency, we propose using the two-dimensional look-up table as a substitution of Equation ( 2) to calculate the height likelihood probability.It is a regular sampling of the continuous likelihood probability that is to replace the continuous function surface with a discrete numerical table.The height likelihood probability is related to the phase, coherence, and effective number of looks.For a multi-baseline InSAR dataset, the effective number of looks of each interferogram is the same and hence the look-up table to be established is indexed by phase and coherence.The value range of the phase is [−π, π], and the sampling interval is π/180 rad.The value range of the coherence coefficient is [0,1], and the sampling interval is 1/100.Figure 2a shows a three-dimensional view of the look-up table with an effective number of looks ( 16) and Figure 2b shows the profile perpendicular to the coherence axis.
are used for the height-to-phase conversion, respectively.The conversion error of the rational function model is calculated using the results of the rigorous method as a reference.Table 1 lists the conversion error for different spaceborne SAR data, indicating that the conversion error for the rational function model is completely negligible in practical applications.Therefore, the rational function model can replace the rigorous method for height-to-phase conversion at each resolution unit.2) is used to calculate the height likelihood probability, which is a very complicated expression.In order to improve the computational efficiency, we propose using the two-dimensional look-up table as a substitution of Equation ( 2) to calculate the height likelihood probability.It is a regular sampling of the continuous likelihood probability that is to replace the continuous function surface with a discrete numerical table.The height likelihood probability is related to the phase, coherence, and effective number of looks.For a multi-baseline InSAR dataset, the effective number of looks of each interferogram is the same and hence the look-up table to be established is indexed by phase and coherence.The value range of the phase is [ , ] , and the sampling interval is / 180 π rad.The value range of the coherence coefficient is [0,1], and the sampling interval is 1/100.Figure 2a shows a three-dimensional view of the look-up table with an effective number of looks ( 16) and Figure 2b shows the profile perpendicular to the coherence axis.First, larger steps are used to narrow the candidate height range.Then a smaller step is applied to search for the ML height.Repeat the process until the target height changes less than the given threshold.The specific implementation flow is as follows:

ML Height Estimation with Flexible Search Step Length
From Section 3.1, the ML height satisfying Equation ( 1) can be estimated by searching the candidate height range.The search step length determines the accuracy of the ML height and the amount of computation time.Instead of the fixed step length, we propose using a flexible search step length to ensure the efficiency and accuracy of the ML height searching.
First, larger steps are used to narrow the candidate height range.Then a smaller step is applied to search for the ML height.Repeat the process until the target height changes less than the given threshold.The specific implementation flow is as follows: (1) Obtain the initial height value h 0 from the prior DEM, which is radar-coded to the SAR image coordinates; (2) Set the height search range to [h 0 − ∆H i , h 0 + ∆H i ] and the search step to ∆h i .∆H i can be set to an integral multiple of σ h .The optimal height obtained by the maximum likelihood estimation is The optimal height obtained by the maximum likelihood estimation is h i+1 ; (4) Test whether (h i+1 − h i ) is less than the given threshold.If yes, then stop the search and return h i+1 as the optimal height.If no, then repeat Step (3).
It should be noted that the step length ∆h must be less than half of the minimum height ambiguity to satisfy the Nyquist sampling law and ensure that the correct position of the height likelihood probability peak can be detected.

Geocoding
The geographic coordinates of each cell with estimated height are calculated through the Range Doppler (RD) model with the World Geodetic System 1984 (WGS84) as the reference spheroid.Afterward, the multi-baseline InSAR DEM could be geo-referenced and gridded in the geographic coordinate system (latitude and longitude) or in the projection coordinates such as the universal transverse mercator (UTM) system with regular spacing in the East and North dimensions [11].

Experiments and Results
The purpose of the simulated experiment in Section 4.1 is to testify the improvement of the ML height accuracy with the re-defined prior height probability distribution by comparing the height accuracy of ML height with or without the prior DEM.Then, the proposed processing flow was applied in the L-band ALOS/PALSAR data covering the Mount Tai area of China, as in Section 4.2.

Simulation of SAR Interferograms
An American NED DEM of 10 m resolution is used to simulate three SAR interferograms with different normal baselines and evaluate height accuracy of the generated DEMs.The wavelength of the simulated radar system is 0.031 m (X band).The orbit height is about 600 km.The incidence angle is 35 • .The effective number of looks (ENL) of the interferogram is 16.For the three simulated interferograms, their respective normal baseline, height ambiguity (H amb ), coherence coefficient, and standard deviation (Std.) of phase noise are shown in Table 2.The geometric decorrelation caused by terrain changes is not considered and the temporal decorrelation for the three interferograms is assumed to be the same.Therefore, the coherence differences among the three interferograms are induced only by the normal baseline differences.The standard deviation of the phase noise induced by decorrelation is calculated by Equation (3). Figure 3a shows the hillshade of NED DEM and Figure 3b is obtained by 15 × 15 smoothing filtering of the NED DEM, which is used for calculating the prior height probabilities and lots of detailed topographic information is lost, as in Figure 3c. Figure 3d-f represents the simulated interferograms I, II, and II, which are calculated from W{2πh/H amb }, (W{•} is the wrapping operator) and h is the height value acquired from the NED DEM (Figure 3a). Figure 3g-i shows atmospheric turbulence phase generated by statistical simulation.Figure 3j-l is the simulated interferogram I, II, and III superimposed by the atmospheric phase (g-i), respectively.The simulated interferograms superimposed by the phase noise (Figure 3d-f) are used to generate a multi-baseline InSAR DEM with ML estimation method with or without the prior DEM.The height search range is set to between 500 m and 2000 m and the search step is 1 m.For ML estimation without the prior DEM, the prior height probability is evenly distributed within this range, while for ML estimation with the prior DEM it is calculated through Equation ( 6).The standard deviation of the prior DEM height error DEM σ is calculated from the height error map generated by differentiating the NED DEM (Figure 3a) and the smoothed NED DEM (Figure 3b).DEM σ is 4.7 m by calculation and empirically enlarged to 6 m.
Figure 4 shows the generated DEMs and their error maps for ML estimation with the prior DEM, (Figure 4a,b) and without the prior DEM (Figure 4c,d).We can see the topographic information is

Test of the Impact of the Prior Height on ML Estimation
The simulated interferograms superimposed by the phase noise (Figure 3d-f) are used to generate a multi-baseline InSAR DEM with ML estimation method with or without the prior DEM.The height search range is set to between 500 m and 2000 m and the search step is 1 m.For ML estimation without the prior DEM, the prior height probability is evenly distributed within this range, while for ML estimation with the prior DEM it is calculated through Equation ( 6).The standard deviation of the prior DEM height error σ DEM is calculated from the height error map generated by differentiating the NED DEM (Figure 3a) and the smoothed NED DEM (Figure 3b).σ DEM is 4.7 m by calculation and empirically enlarged to 6 m.
Figure 4 shows the generated DEMs and their error maps for ML estimation with the prior DEM, (Figure 4a,b) and without the prior DEM (Figure 4c,d).We can see the topographic information is almost all covered by the noise value in Figure 4a and the height error map in Figure 4b ranges between −960 m and 960 m.Compare Figure 4b to Figure 3d-f, it can be seen that the spatial distribution of the height errors is still related to the wrapped terrain phase, and hence the ML estimation without the prior DEM does not solve the height ambiguity problem and is very sensitive to phase noise.The DEM obtained with ML estimation with the prior DEM depicts the terrain condition well in Figure 4c.The height error map ranges between −8 m and 8 m and is almost randomly spatially distributed except that it is slightly larger in the fluctuating area.From all the above comparisons, it is obvious that the height estimation accuracy and anti-noise ability of ML estimation with the prior DEM are much better than those of ML estimation without the prior DEM.The statistical values of the height errors in the simulation experiment are calculated as shown in Table 3.The theoretical error of the single-baseline interferogram is caused by the decorrelation phase noise without the phase unwrapping error.From Table 3, the Std. of the height errors are up to 408.8 m for ML estimated height without the prior DEM while shrink to 1.6 m for ML estimated height with the prior DEM.The height accuracy of ML estimation with the prior DEM is better than that of the single-baseline InSAR DEMs in terms of the standard deviation of height errors.
In Figure 3b, a small area in the black rectangle is selected and enlarged as shown in Figure 3c.The upper part of Figure 3c is the hillshade of the 10 m NED DEM and the lower part of Figure 3c is the hillshade of the 15 × 15 smoothing filtered NED DEM.Most of the terrain details are lost in the 15 × 15 smoothing filtered NED DEM.The same area of the ML-estimated InSAR DEM is also selected and enlarged in Figure 4e.By comparative analysis of Figures 3c and 4e, it can be seen that although there is high-frequency random noise (caused by decorrelation noise) in the ML-estimated DEM in Figure 4e, the ML-estimated DEM well reconstructs the topographic details lost in the prior DEM as in the lower part of Figure 3c and obviously improves the spatial resolution of the prior DEM.The Std. of height error of ML estimation with a prior DEM (1.6 m) is less than that of prior DEM (4.7 m) in Table 3.The statistical values of the height errors in the simulation experiment are calculated as shown in Table 3.The theoretical error of the single-baseline interferogram is caused by the decorrelation phase noise without the phase unwrapping error.From Table 3, the Std. of the height errors are up to 408.8 m for ML estimated height without the prior DEM while shrink to 1.6 m for ML estimated height with the prior DEM.The height accuracy of ML estimation with the prior DEM is better than that of the single-baseline InSAR DEMs in terms of the standard deviation of height errors.

Mean
Std.In Figure 3b, a small area in the black rectangle is selected and enlarged as shown in Figure 3c.The upper part of Figure 3c is the hillshade of the 10 m NED DEM and the lower part of Figure 3c is the hillshade of the 15 × 15 smoothing filtered NED DEM.Most of the terrain details are lost in the 15 × 15 smoothing filtered NED DEM.The same area of the ML-estimated InSAR DEM is also selected and enlarged in Figure 4e.By comparative analysis of Figures 3c and 4e, it can be seen that although there is high-frequency random noise (caused by decorrelation noise) in the ML-estimated DEM in Figure 4e, the ML-estimated DEM well reconstructs the topographic details lost in the prior DEM as in the lower part of Figure 3c and obviously improves the spatial resolution of the prior DEM.The Std. of height error of ML estimation with a prior DEM (1.6 m) is less than that of prior DEM (4.7 m) in Table 3.

Test of the Impact of the Atmospheric Effects on ML Estimation
Figure 5 shows the multi-baseline InSAR DEM generated from interferograms with the atmospheric effects shown in Figure 3j-l and the corresponding height error map.Comparing the height error map Figure 5b to Figure 4d, it can be seen that there is a trend towards height error in Figure 5b, caused by the atmospheric effects.Figure 5 shows the multi-baseline InSAR DEM generated from interferograms with the atmospheric effects shown in Figure 3j-l and the corresponding height error map.Comparing the height error map Figures 5b to 4d, it can be seen that there is a trend towards height error in Figure 5b, caused by the atmospheric effects.Furthermore, Table 4 shows the statistical values of height errors of the single and multi-baseline InSAR DEMs generated from interferograms with the atmospheric effects in Figure 3j-l.For both single-baseline interferometry and multi-baseline InSAR height estimation, the atmospheric effects can evidently increase the height errors according to Table 4.The multi-baseline InSAR DEM has better height accuracy than the single-baseline InSAR DEMs, indicating that the ML estimation with prior DEM can effectively suppress the atmospheric effects to some extent.However, comparing Tables 3 with 4, it can be seen the standard height error of the ML estimation with prior DEM has increased from 1.6 m to 4.1 m, revealing that the atmospheric effects can cause non-negligible height errors and should be removed from each interferogram.

Experimental Area
The experimental area covers the central and southern parts of Mount Tai (including its main peak, Yuhuangding), and the the central plains and hills (including Taian city at the southern foot of Mount Tai) of Shandong Province.A large portion of this area is covered by vegetation.From the optical image acquired from Google Earth in Figure 6a, we can see that the northern part of the experimental area, including Mount Tai and hills, and Mount Culai in the southeast corner, are Furthermore, Table 4 shows the statistical values of height errors of the single and multi-baseline InSAR DEMs generated from interferograms with the atmospheric effects in Figure 3j-l.For both single-baseline interferometry and multi-baseline InSAR height estimation, the atmospheric effects can evidently increase the height errors according to Table 4.The multi-baseline InSAR DEM has better height accuracy than the single-baseline InSAR DEMs, indicating that the ML estimation with prior DEM can effectively suppress the atmospheric effects to some extent.However, comparing Table 3 with Table 4, it can be seen the standard height error of the ML estimation with prior DEM has increased from 1.6 m to 4.1 m, revealing that the atmospheric effects can cause non-negligible height errors and should be removed from each interferogram.The experimental area covers the central and southern parts of Mount Tai (including its main peak, Yuhuangding), and the the central plains and hills (including Taian city at the southern foot of Mount Tai) of Shandong Province.A large portion of this area is covered by vegetation.From the optical image acquired from Google Earth in Figure 6a, we can see that the northern part of the experimental area, including Mount Tai and hills, and Mount Culai in the southeast corner, are basically covered by woods.Moreover, there is a large amount of farmland in the central and southern parts of the experimental area, which is covered by crops.In this vegetation-covered area, the L-band ALOS/PALSAR data can maintain better temporal coherence and penetrate the vegetation cover to some extent.The terrain types of the experimental area are diverse, including plains, hills, and mountains.The height varies greatly such that the height of the plains is about 100 m above the sea level while the height of Mount Tai peak is about 1533 m.Therefore, we can say it is a suitable experimental area for multi-baseline InSAR DEM generation experiment.

ALOS/PALSAR Data
The SAR interferometric data used in this experiment are ALOS/PALSAR data obtained in the FBS mode and have a total of six images.Table 5 lists the detailed image parameters.The amplitude image shown in Figure 6b is acquired on 6 February 2008 and has been through multi-look processing with a ratio of 7:3 (azimuth:range).Since the data are acquired from ascending orbital direction with right-looking imaging, the amplitude image in Figure 6b is flipped over vertically compared with the image coverage in Figure 6a.The ALOS/PALSAR multi-baseline dataset used in the experiment has six images and Figure 7 shows the spatial-temporal baseline distribution of the interferometric dataset.We selected four interferometric pairs based on the criteria given in Section 3.1, connected by the blue lines as shown in Figure 7.It can be seen that the temporal baselines of interferometric pairs are all about 46 d.For

ALOS/PALSAR Data
The SAR interferometric data used in this experiment are ALOS/PALSAR data obtained in the FBS mode and have a total of six images.Table 5 lists the detailed image parameters.The amplitude image shown in Figure 6b is acquired on 6 February 2008 and has been through multi-look processing with a ratio of 7:3 (azimuth:range).Since the data are acquired from ascending orbital direction with right-looking imaging, the amplitude image in Figure 6b is flipped over vertically compared with the image coverage in Figure 6a.The ALOS/PALSAR multi-baseline dataset used in the experiment has six images and Figure 7 shows the spatial-temporal baseline distribution of the interferometric dataset.We selected four interferometric pairs based on the criteria given in Section 3.1, connected by the blue lines as shown in Figure 7.It can be seen that the temporal baselines of interferometric pairs are all about 46 d.For interferometric pairs 1-2 and 2-3, image 2 is chosen as the master image and for interferometric pairs 4-5 and 5-6, image 5 is chosen as the master image.Then we register interferograms 4-5 and 5-6 to the image space of image 2, as shown by the red line in Figure 7. Table 6 shows the parameters of the four interferograms; the distribution of the normal baselines ranges from −784 m to 185 m and the corresponding height ambiguity ranges from 82 m to 833 m.These are the suitable normal baseline combinations for multi-baseline InSAR processing.In order to improve the signal-to-noise ratio, we perform 7 × 3 multi-look processing of the SAR images in the azimuth and range direction, and then the spatial resolution of the interferogram is about 22 m × 22 m.

Elevation Datasets
In the experiment, the SRTM DEM with a spatial resolution of 90 m × 90 m is used as the prior DEM and the standard deviation of the global height error of SRTM DEM is about 5 m [21].In order to validate the accuracy of the single/multi-baseline InSAR DEMs, the 1:25,000 aerial photogrammetric DEM provided by the Shandong Provincial Land and Surveying and Mapping Institute is used as the reference DEM with the root mean square error (RMSE) less than 4 m [27].The spatial resolution of 1:25,000 DEM is 10 m × 10 m.The area of the 1:25,000 DEM coverage is about 12 km × 14 km marked as the green rectangle in Figure 6a, including the southern part of Mount Tai and Tai'an City.

Experimental Results
Figure 8 shows the flattened interferograms of the ALOS/PALSAR data.The interferometric fringes are sparse in the interferogram with short baseline in Figure 8b,d, and become very dense in Table 6 shows the parameters of the four interferograms; the distribution of the normal baselines ranges from −784 m to 185 m and the corresponding height ambiguity ranges from 82 m to 833 m.These are the suitable normal baseline combinations for multi-baseline InSAR processing.In order to improve the signal-to-noise ratio, we perform 7 × 3 multi-look processing of the SAR images in the azimuth and range direction, and then the spatial resolution of the interferogram is about 22 m × 22 m.

Elevation Datasets
In the experiment, the SRTM DEM with a spatial resolution of 90 m × 90 m is used as the prior DEM and the standard deviation of the global height error of SRTM DEM is about 5 m [21].In order to validate the accuracy of the single/multi-baseline InSAR DEMs, the 1:25,000 aerial photogrammetric DEM provided by the Shandong Provincial Land and Surveying and Mapping Institute is used as the reference DEM with the root mean square error (RMSE) less than 4 m [27].The spatial resolution of 1:25,000 DEM is 10 m × 10 m.The area of the 1:25,000 DEM coverage is about 12 km × 14 km marked as the green rectangle in Figure 6a, including the southern part of Mount Tai and Tai'an City.

Experimental Results
Figure 8 shows the flattened interferograms of the ALOS/PALSAR data.The interferometric fringes are sparse in the interferogram with short baseline in Figure 8b,d, and become very dense in the interferogram with long baseline in Figure 8a,c.The interferograms are clearly visible even in the mountainous and vegetation-covered areas, indicating that the L-band SAR data have good capability of keeping geometric and temporal coherence.Figure 9 shows the hillshades of the single-/multi-baseline DEMs and SRTM DEM corresponding to the black rectangle in Figure 8a.Besides the multi-baseline InSAR DEM (Figure 8a), the radar-coded SRTM DEM (Figure 8b) and single-baseline InSAR DEMs (Figure 8c-f) of the same area are also presented as comparisons.
As to quantitatively evaluate the height accuracy with the reference 1:25,000 DEM, all the InSAR DEMs and are geocoded into the geographic coordinate system and then projected into the UTM coordinate system with spatial resolution of 20 m.Since the spatial resolution of 1:25,000 DEM is about 10 m, it needs to be down-resampled to 20 m when calculating the height error map.Moreover, in order to verify the height accuracy improvement of the multi-baseline estimation, the SRTM DEM and single-baseline InSAR DEMs are also evaluated as comparison.Figure 10 shows the hillshade of the generated multi-baseline InSAR DEM with height ranging between 20 m and 1535 m, which clearly shows the topographic conditions of Mount Tai, Mount Culai (located at the southeast corner) and a large number of hills.In the calculation of the prior height probability, σ SRTM is empirically set to 20 m (considering the spatial resolution difference and the interpolation error of the radar coding) with the neighborhood system N = 8.When calculating the height likelihood probability, the initial height value is obtained from the radar-coded SRTM DEM and the initial search height range and the step size are set at ±300 m and 20 m, respectively.
Figure 9 shows the hillshades of the single-/multi-baseline DEMs and SRTM DEM corresponding to the black rectangle in Figure 8a.Besides the multi-baseline InSAR DEM (Figure 8a), the radar-coded SRTM DEM (Figure 8b) and single-baseline InSAR DEMs (Figure 8c-f) of the same area are also presented as comparisons.
As to quantitatively evaluate the height accuracy with the reference 1:25,000 DEM, all the InSAR DEMs and are geocoded into the geographic coordinate system and then projected into the UTM coordinate system with spatial resolution of 20 m.Since the spatial resolution of 1:25,000 DEM is about 10 m, it needs to be down-resampled to 20 m when calculating the height error map.Moreover, in order to verify the height accuracy improvement of the multi-baseline estimation, the SRTM DEM and single-baseline InSAR DEMs are also evaluated as comparison.Figure 10 shows the hillshade of the generated multi-baseline InSAR DEM with height ranging between 20 m and 1535 m, which clearly shows the topographic conditions of Mount Tai, Mount Culai (located at the southeast corner) and a large number of hills.

Comparative Analysis of the Single-and Multi-Baseline InSAR DEMs
For the single-baseline InSAR, the interferograms in Figure 7 with longer normal baselines (Figure 7a,c) have denser fringes than the interferograms with shorter normal baselines (Figure 7b,d).After the phase unwrapping and phase-to-height conversion, the single-baseline InSAR DEMs with longer normal baselines (Figure 8c,e) have more topographic detail than InSAR DEMs with shorter normal baselines (Figure 8d,f).Quantitatively, the height error of single-baseline InSAR DEMs with longer normal baselines (Figure 10c,e) is evidently smaller than the InSAR DEMs with shorter normal baselines (Figure 10d,f), which can also be reflected in Table 6-the interferograms I, III with longer normal baselines a have smaller standard deviation of height errors and a larger percentage of the absolute height error less than 10 m than the interferograms II and IV with shorter normal baselines.Then, can it be concluded that the longer the normal baseline, the better height measurement accuracy for single-baseline InSAR?The answer is obviously no.Interferogram I (Figure 7a) has a longer normal baseline than interferogram III (Figure 7c); however, the height accuracy of the InSAR DEM generated by interferogram II, as in Figures 8e and 10e, is better than that of the InSAR DEM generated by interferogram I, as in Figures 8c and 10c.According to Equation (4), the height accuracy of single-baseline InSAR is determined by both the length of the normal baseline and the phase noise, while too long a normal baseline can also introduce non-negligible phase noise during the delicate phase unwrapping.
In this article, this contradiction is solved by a multi-baseline InSAR with maximum likelihood estimation criterion.In the ALOS/PALSAR data experiment, the multi-baseline InSAR DEM (Figure 9a) have more topographic detail than the single-baseline InSAR DEMs (Figure 9c-f).The height error

Comparative Analysis of the Single-and Multi-Baseline InSAR DEMs
For the single-baseline InSAR, the interferograms in Figure 7 with longer normal baselines (Figure 7a,c) have denser fringes than the interferograms with shorter normal baselines (Figure 7b,d).After the phase unwrapping and phase-to-height conversion, the single-baseline InSAR DEMs with longer normal baselines (Figure 8c,e) have more topographic detail than InSAR DEMs with shorter normal baselines (Figure 8d,f).Quantitatively, the height error of single-baseline InSAR DEMs with longer normal baselines (Figure 10c,e) is evidently smaller than the InSAR DEMs with shorter normal baselines (Figure 10d,f), which can also be reflected in Table 6-the interferograms I, III with longer normal baselines a have smaller standard deviation of height errors and a larger percentage of the absolute height error less than 10 m than the interferograms II and IV with shorter normal baselines.Then, can it be concluded that the longer the normal baseline, the better height measurement accuracy for single-baseline InSAR?The answer is obviously no.Interferogram I (Figure 7a) has a longer normal baseline than interferogram III (Figure 7c); however, the height accuracy of the InSAR DEM generated by interferogram II, as in Figures 8e and 10e, is better than that of the InSAR DEM generated by interferogram I, as in Figures 8c and 10c.According to Equation (4), the height accuracy of single-baseline InSAR is determined by both the length of the normal baseline and the phase noise, while too long a normal baseline can also introduce non-negligible phase noise during the delicate phase unwrapping.
In this article, this contradiction is solved by a multi-baseline InSAR with maximum likelihood estimation criterion.In the ALOS/PALSAR data experiment, the multi-baseline InSAR DEM (Figure 9a) more topographic detail than the single-baseline InSAR DEMs (Figure 9c-f).The height error of the single-baseline InSAR DEMs (Figure 10c-f) is larger than that of the multi-baseline InSAR DEM (Figure 10a).In Table 6, the standard deviation of height error of the multi-baseline InSAR DEM is the smallest and the percentage of the absolute height error less than 10 m of the multi-baseline InSAR DEM is the largest among all the DEMs.Therefore, the height accuracy of the multi-baseline DEM is better than that of the single-baseline DEMs on the whole, indicating that the proposed method is effective.

Comparative Analysis of the Prior Height's Impact on ML Estimation
From the simulated data experiment in Section 4.1.2,the ML estimated DEM with the prior height probability distribution in Figure 4c,d has higher height accuracy than the ML estimated DEM without the prior height probability distribution in Figure 4a,b.Here, we analyze the influence of the prior height on ML estimation.A resolution unit in the simulated interferogram with true height of 1320.4 m is used as an example.Figure 12 shows the height probability distribution of the chosen resolution unit.Figure 12a shows the joint height probability distribution of three simulated interferograms without phase noise; when h is 1320 m, the joint probability density function (PDF) reaches its peak.Figure 12b shows the joint PDF with the decorrelation noise and, although there is a local peak at h = 1320 m, the highest peak is located at h = 1469 m with height error up to 149 m. Figure 12c shows the prior height probability distribution, with the highest peak located at h = 1332 m. Figure 12d shows the joint height PDF by multiplying the height PDF in Figure 12b and the prior height distribution in Figure 12c.There is only one peak in Figure 12d of the single-baseline InSAR DEMs (Figure 10c-f) is larger than that of the multi-baseline InSAR DEM (Figure 10a).In Table 6, the standard deviation of height error of the multi-baseline InSAR DEM is the smallest and the percentage of the absolute height error less than 10 m of the multi-baseline InSAR DEM is the largest among all the DEMs.Therefore, the height accuracy of the multi-baseline DEM is better than that of the single-baseline DEMs on the whole, indicating that the proposed method is effective.

Comparative Analysis of the Prior Height's Impact on ML Estimation
From the simulated data experiment in Section 4.1.2,the ML estimated DEM with the prior height probability distribution in Figure 4c,d has higher height accuracy than the ML estimated DEM without the prior height probability distribution in Figure 4a,b.Here, we analyze the influence of the prior height on ML estimation.A resolution unit in the simulated interferogram with true height of 1320.4 m is used as an example.Figure 12 shows the height probability distribution of the chosen resolution unit.Figure 12a shows the joint height probability distribution of three simulated interferograms without phase noise; when h is 1320 m, the joint probability density function (PDF) reaches its peak.Figure 12b shows the joint PDF with the decorrelation noise and, although there is a local peak at 1320 h = m, the highest peak is located at 1469 h = m with height error up to 149 m. Figure 12c shows the prior height probability distribution, with the highest peak located at 1332 h = m.Figure 12d shows the joint height PDF by multiplying the height PDF in Figure 12b and the prior height distribution in Figure 12c.There is only one peak in Figure 12d

Comparative Analysis of the Multi-Baseline InSAR DEM and SRTM DEM
The height accuracy of the multi-baseline DEM and SRTM DEM is comparatively evaluated.The multi-baseline InSAR DEM (Figure 8a) has more topographic details than the radar-coded SRTM DEM (Figure 8b).In Figure 10, the height error of multi-baseline InSAR DEM (Figure 10a) is obviously smaller than that of SRTM DEM (Figure 10b), especially in the mountainous areas.From the statistical Compared with SRTM DEM, the multi-baseline InSAR DEM has obvious advantages in terms of resolution and precision.Hence the multi-baseline InSAR estimation can be viewed as a topographical information update of the historical low-resolution DEMs.(4) The multi-baseline InSAR DEM generated from ALOS/PALSAR datasets meets the American DTED-2 standard and Chinese 1:50,000 DEM (mountain) Level 2 in the case of spatial resolution and height accuracy.
In the future, the proposed multi-baseline InSAR topographic mapping flow can be tested and improved further with more spaceborne datasets.The maximum likelihood height estimation method with the prior DEM provides a promising solution for DEM mass production in mountainous areas.

Figure 2 .
Figure 2. (a) The fitted surface by the look-up table with ENL L = 16 and (b) the profile perpendicular to the coherence axis with 0.3, 0.6, 0.9 ρ =

Figure 2 .
Figure 2. (a) The fitted surface by the look-up table with ENL L = 16 and (b) the profile perpendicular to the coherence axis with ρ = 0.3, 0.6, 0.9.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 20 topographic information is lost, as in Figure3c.Figure3d-f represents the simulated interferograms I, II, and II, which are calculated from {2 / } W ⋅ is the wrapping operator) and h is the height value acquired from the NED DEM (Figure3a).Figure3g-i shows atmospheric turbulence phase generated by statistical simulation.Figure3j-l is the simulated interferogram I, II, and III superimposed by the atmospheric phase (g-i), respectively.

Figure 3 .
Figure 3. (a) Hillshade of 10 m NED DEM; (b) hillshade of 15 × 15 smoothing filtered NED DEM; (c) NED DEM (upper) and smoothing filtered NED DEM (lower); (d-f) are the simulated interferograms I, II, and III, respectively, and the phase noise introduced by decorrelation is superimposed; (g-i) are the simulated atmospheric phase for interferograms I, II, and III, respectively; (j-l) are the simulated interferograms I, II, and III superimposed by the atmospheric phase (g-i), respectively.

Figure 3 .
Figure 3. (a) Hillshade of 10 m NED DEM; (b) hillshade of 15 × 15 smoothing filtered NED DEM; (c) NED DEM (upper) and smoothing filtered NED DEM (lower); (d-f) are the simulated interferograms I, II, and III, respectively, and the phase noise introduced by decorrelation is superimposed; (g-i) are the simulated atmospheric phase for interferograms I, II, and III, respectively; (j-l) are the simulated interferograms I, II, and III superimposed by the atmospheric phase (g-i), respectively.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 20almost all covered by the noise value in Figure4aand the height error map in Figure4branges between −960 m and 960 m.Compare Figure4bto Figure3d-f, it can be seen that the spatial distribution of the height errors is still related to the wrapped terrain phase, and hence the ML estimation without the prior DEM does not solve the height ambiguity problem and is very sensitive to phase noise.The DEM obtained with ML estimation with the prior DEM depicts the terrain condition well in Figure4c.The height error map ranges between −8 m and 8 m and is almost randomly spatially distributed except that it is slightly larger in the fluctuating area.From all the above comparisons, it is obvious that the height estimation accuracy and anti-noise ability of ML estimation with the prior DEM are much better than those of ML estimation without the prior DEM.

Figure 4 .
Figure 4. Multi-baseline InSAR DEMs (without atmospheric effects) and height error maps, without the prior DEM (a,b) and with the prior DEM (c-e) is the hillshade of (c) with the enlarged topographic details of a small area.

Figure 4 .
Figure 4. Multi-baseline InSAR DEMs (without atmospheric effects) and height error maps, without the prior DEM (a,b) and with the prior DEM (c-e) is the hillshade of (c) with the enlarged topographic details of a small area.
Interferogram III (Figure 3f) −0.001 m 2.0 m ML without prior DEM 70.072 m 408.8 m ML with prior DEM −0.003 m 1.6 m

3 .
Test of the Impact of the Atmospheric Effects on ML Estimation

Figure 5 .
Figure 5. (a) Multi-baseline InSAR DEM with prior DEM (with atmospheric effects); (b) height error map; (c) is the hillshade map of (a).
ML with prior DEM −0.006 m 4.1 m

Figure 5 .
Figure 5. (a) Multi-baseline InSAR DEM with prior DEM (with atmospheric effects); (b) height error map; (c) is the hillshade map of (a).
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 20 the L-band ALOS/PALSAR data can maintain better temporal coherence and penetrate the vegetation cover to some extent.The terrain types of the experimental area are diverse, including plains, hills, and mountains.The height varies greatly such that the height of the plains is about 100 m above the sea level while the height of Mount Tai peak is about 1533 m.Therefore, we can say it is a suitable experimental area for multi-baseline InSAR DEM generation experiment.

Figure 6 .
Figure 6.(a) The coverage of ALOS/PALSAR images marked by the blue rectangle and 1:25,000 DEM marked by the green rectangle shown in Google Earth; (b) the amplitude image of ALOS/PALSAR data (acquisition time: 6 February 2008).

Figure 6 .
Figure 6.(a) The coverage of ALOS/PALSAR images marked by the blue rectangle and 1:25,000 DEM marked by the green rectangle shown in Google Earth; (b) the amplitude image of ALOS/PALSAR data (acquisition time: 6 February 2008).

20 Figure 7 .
Figure 7.The temporal and spatial baseline distribution of ALOS/PALSAR interferometric pairs.The images connected by the blue line constitute interferometric pairs and the interferograms connected by the red line are co-registered to the same image space.

Figure 7 .
Figure 7.The temporal and spatial baseline distribution of ALOS/PALSAR interferometric pairs.The images connected by the blue line constitute interferometric pairs and the interferograms connected by the red line are co-registered to the same image space.

Figure 10 .
Figure 10.Hillshade of multi-baseline InSAR DEM.The black rectangle marks the coverage of the 1:25,000 DEM used for accuracy validation.

Figure 11
Figure11shows the height error maps of single/multi-baseline InSAR DEMs and SRTM DEM, which directly reflect the distribution and amount of the height error for different DEMs.Table7shows the statistical values of the height error maps of the corresponding InSAR DEMs and SRTM DEM shown in Figure10, which can clearly reflect the statistical characteristics of the height error distribution.

Figure 10 .
Figure 10.Hillshade of multi-baseline InSAR DEM.The black rectangle marks the coverage of the 1:25,000 DEM used for accuracy validation.

Figure 11
Figure11shows the height error maps of single/multi-baseline InSAR DEMs and SRTM DEM, which directly reflect the distribution and amount of the height error for different DEMs.Table7shows the statistical values of the height error maps of the corresponding InSAR DEMs and SRTM DEM shown in Figure10, which can clearly reflect the statistical characteristics of the height error distribution.

Figure 10 .
Figure 10.Hillshade of multi-baseline InSAR DEM.The black rectangle marks the coverage of the 1:25,000 DEM used for accuracy validation.

Figure 11
Figure11shows the height error maps of single/multi-baseline InSAR DEMs and SRTM DEM, which directly reflect the distribution and amount of the height error for different DEMs.Table7shows the statistical values of the height error maps of the corresponding InSAR DEMs and SRTM DEM shown in Figure10, which can clearly reflect the statistical characteristics of the height error distribution.

Figure 11 .Table 7 .
Figure 11.The height error map of InSAR DEMs, (a) the DEM generated by ML estimation with prior DEM; (b) SRTM DEM; (c-f) are the DEMs generated by single-baseline interferograms I, II, III, and IV.

Figure 11 .
Figure 11.The height error map of InSAR DEMs, (a) the DEM generated by ML estimation with prior DEM; (b) SRTM DEM; (c-f) are the DEMs generated by single-baseline interferograms I, II, III, and IV.
, located at h = 1320 m, very close to the true height of h = 1320.4m.Therefore, we can say that the prior DEM forms an external height constraint on the multi-baseline InSAR observations and can effectively remove the erroneous estimate and solve the height ambiguity problem.Remote Sens. 2018, 10, x FOR PEER REVIEW 17 of 20 , located at 1320 h = m, very close to the true height of 1320.4 h = m.Therefore, we can say that the prior DEM forms an external height constraint on the multi-baseline InSAR observations and can effectively remove the erroneous estimate and solve the height ambiguity problem.

Figure 12 .
Figure 12.The height probability distribution of a cell in the simulated interferogram with a true height of 1320.4 h = m.(a,b) The joint height probability distribution with and without decorrelation phase noise, respectively; (c) the prior height probability distribution acquired from the prior DEM; (d) the joint height probability distribution with the prior height probability distribution, obtained by multiplying the two probability distributions in (b,c).

Figure 12 .
Figure 12.The height probability distribution of a cell in the simulated interferogram with a true height of h = 1320.4m. (a,b) The joint height probability distribution with and without decorrelation phase noise, respectively; (c) the prior height probability distribution acquired from the prior DEM; (d) the joint height probability distribution with the prior height probability distribution, obtained by multiplying the two probability distributions in (b,c).

Table 1 .
Height-to-phase conversion errors of rational function model.

Table 1 .
Height-to-phase conversion errors of rational function model.

Table 2 .
The simulation parameters of the interferograms.

Table 3 .
Statistical values of height errors without atmospheric effects.

Table 3 .
Statistical values of height errors without atmospheric effects.

Table 4 .
Statistical values of height errors with atmospheric effects.

Table 4 .
Statistical values of height errors with atmospheric effects.

Table 7 .
Statistical values of height errors of single/multi-baseline InSAR DEMs and SRTM DEM.