Using Dual-Polarization Interferograms to Correct Atmospheric Effects for InSAR Topographic Mapping

Atmospheric effect represents one of the major error sources for interferometric synthetic aperture radar (InSAR), particularly for the repeat-pass InSAR data. In order to further improve the practicability of InSAR technology, it is essential to study how to estimate and eliminate the undesired impact of atmospheric effects. In this paper, we propose the multi-resolution weighted correlation analysis (MRWCA) method between the dual-polarization InSAR data to estimate and correct atmospheric effects for InSAR topographic mapping. The study is based on the a priori knowledge that atmospheric effects is independent of the polarization. To find the identical atmospheric phase (ATP) signals of interferograms in different polarizations, we need to remove the other same or similar phase components. Using two different topographic data, differential interferometry was firstly performed so that the obtained differential interferograms (D-Infs) have different topographic error phases. A polynomial fitting method is then used to remove the orbit error phases. Thus, the ATP signals are the only identical components in the final obtained D-Infs. By using a forward wavelet transform, we break down the obtained D-Infs into building blocks based on their frequency properties. We then applied weighted correlation analysis to estimate the wavelet coefficients attributed to the atmospheric effects. Thus, the ATP signals can be obtained by the refined wavelet coefficients during inverse wavelet transform (IWT). Lastly, we tested the proposed method by the L-band Advanced Land Observing Satellite (ALOS)-1 PALSAR dual-polarization SAR data pairs covering the San Francisco (USA) and Moron (Mongolia) regions. By using Ice, Cloud, and land Elevation Satellite (ICESat) data as the reference data, we evaluated the vertical accuracy of the InSAR digital elevation models (DEMs) with and without atmospheric effects correction, which shows that, for the San Francisco test site, the corrected interferogram could provide a DEM with a root-mean-square error (RMSE) of 7.79 m, which is an improvement of 40.5% with respect to the DEM without atmospheric effects correction. For the Moron test site, the corrected interferogram could provide a DEM with an RMSE of 10.74 m, which is an improvement of 30.2% with respect to the DEM without atmospheric effects correction.


Introduction
Interferometric synthetic aperture radar (InSAR) is a powerful technology for Earth observations, which is noted for its all-weather and day-and-night working capability, wide spatial coverage, and fine resolution [1].In particular, the repeat-pass interferometric mode has been widely employed to acquire InSAR data [2].For repeat-pass synthetic aperture radar (SAR) interferograms, different image acquisition times lead to different atmospheric delays for different SAR images, which result in significant atmospheric error in the interferometric phases.As an inevitable error source, the presence of atmospheric effects is a major problem in repeat-pass InSAR topographic mapping measurement [3].
In recent years, to estimate and correct atmospheric effects, three main types of methods have been proposed [4][5][6][7][8][9][10]: (1) Correlation analysis between interferometric phases and elevation [4,5]: Correlation analysis can be applied to remove the atmospheric phase (ATP) signals strongly correlated with topography.However, it is difficult to achieve satisfactory results with the interferograms acquired over plain areas.
(2) Calibration with external data (e.g., GPS data [6], or Moderate Resolution Imaging Spectroradiometer (MODIS) [7] and Medium Resolution Imaging Spectrometer (MERIS) water vapor data [8]): This method based on external data has attracted increasing attention in recent years since the accuracy and resolution of the water vapor data are improved.However, we need to solve the problems of the inconsistent resolution and nonsynchronous acquisition times between the external data and SAR images.
(3) Spatio-temporal analysis by multi-baseline InSAR data [9,10]: spatio-temporal analysis can be used to mitigate the atmospheric effects without relying on the relationship between the ATP signals and the topography or auxiliary data.However, to perform reliable spatio-temporal analysis, we need a large volume of repeat-pass SAR data covering the same area.
Owing to the requirements for external auxiliary data or large volume of SAR data, it is difficult to use the above-mentioned methods to correct the atmospheric effects for accurate InSAR topographic mapping.As the ATP signals are independent of the polarization, multi-polarization SAR images contain the same atmospheric delay [11].As a result, this provides a new idea for us to detect and remove the identical ATP components by use of multi-polarization interferograms.To achieve this goal, we propose the multi-resolution weighted correlation analysis (MRWCA) method between dual-polarization interferograms to estimate and remove the atmospheric effects for InSAR digital elevation model (DEM) generation.The main advantage of the MRWCA method is that it can remove the ATP signals for single-baseline InSAR data without depending on any a prior knowledge or external water vapor data.
This paper is organized as follows.We begin by describing the phase components of differential interferograms (D-Infs) and the MRWCA method in Section 2. We then test the proposed method with real SAR data acquired over different topographic conditions in Section 3. In Section 4, we discuss the characteristics of the MRWCA method and some conclusions are drawn in the last section.

Methodology
When two interferograms acquired in different polarizations are available, it provides us with the possibility to estimate the common ATP signals by correlation analysis.Before performing this, it is necessary to analyze their phase components and any other similar or identical phase components should be removed.

Phase in the Differential Interferograms
The phase component of an interferogram is composed of flat-earth phase, topographic phase, ATP and noise phase (NP) [12].In such a case, it is difficult to identify the ATP signals directly from the interferogram.In order to detect the ATP signals in the dual-polarization interferograms, differential interferometry is adopted to generate D-Infs whose phases are decomposed as: where 'p1' and 'p2' are the two selected polarizations (e.g., HH, HV, VV [13]).φ p εtopo is the topographic error phase (TEP), which is caused by the height error of auxiliary DEM used to simulate the topographic phase.φ atm is the ATP, and φ orbit is the orbit error phase caused by the uncertainty in the orbit parameters.φ p noise is the NP that is mainly attributed to thermal noise and de-correlation effect.φ p de f and φ p iono are the ground deformation and ionosphere phases, respectively, and are not considered in our study.
Among the four phase components, except the ATP, the orbit error phase is also polarization-independent, which should be removed [14].The orbit error phase is usually shown as a systematic phase trend and can be removed by a polynomial model incorporating elevation information [10,15].In addition, since the SAR signals in different polarizations are sensitive to different scattering mechanisms, they are characterized by different phase centers and NP signals.In particular, the different phase centers will result in different TEPs φ p 1 εtopo and φ p 2 εtopo .However, if common external DEM is used to simulate the topographic phase for the dual-polarization interferograms, TEPs φ p 1 εtopo and φ p 2 εtopo will be similar because they have the phase components caused by the identical ground elevation errors.In order to solve this problem, two different external DEMs are used to simulate the topographic phases for the two interferograms, respectively.Thus, the only common component in the two obtained D-Infs is the ATP φ atm .

Estimation of the Atmospheric Effects
Although previous studies have demonstrated that the TEP and the NP signals are composed by short-wavelength components, while the wavelengths of the ATP signals can range from tens of meters to tens of kilometers [5].As a result, it is difficult to directly separate the ATP and non-ATP (TEP and NP) signals from the D-Inf only by filtering components characterized by different wavelengths.However, for the phase components characterized by different wavelengths, which are composed by different ATP and non-ATP signals.In such a case, it provides us the possibility to detect the common ATP signals mixed with different non-ATP signals.To achieve this goal, multi-resolution analysis performed by wavelet can be considered.
To begin with, in the wavelet domain, a D-Inf can be expressed as multiple components characterized by different wavelengths.Consequently, we can investigate the spatial distribution of the D-Inf at different decomposition scales.In this paper, we introduce the multi-resolution analysis [16] to the two D-Infs.Assuming φ p i di f f (x, y) to be a 2D image with a size of m × n, the two D-Infs can be decomposed in the wavelet domain as follows [16]: where Φ and Ψ are the smoothing and the mother wavelet functions, respectively; p 1 w and p 2 w are the wavelet coefficients; J is the number of decomposition scales, and ε = 1, 2, 3 represents the horizontal (H), vertical (V), and diagonal (D) detail information, respectively.The high-frequency components of the two D-Infs, which contain the common ATP signals, have the same or similar values.To evaluate this similarity, we calculate the linear relationship between the wavelet coefficients in different polarizations for each scale j as follows: where f is the slope and ranges from 0 to 1, which is called the correlation coefficient since it indicates the similarity of the ATP signals within the two polarimetric interferograms; and c j represents a constant bias term.The larger the slope, the more similar the ATP components.In this paper, least squares (LS) estimation is used to calculate the slope f and constant bias c j .Although f can determine the ratio of the ATP signals in the obtained D-Infs, it only expresses the global correlation relationship.In such a case, the correlation coefficient f has low sensitivity to the ATP variation since the ATP ratios are different for each pixel.To solve this problem, the perpendicular distance from the wavelet coefficient point to the fitting line is then used to formulate a weighted indicator [17] and the correlation coefficient f is refined by: where d i is the perpendicular distance from the ith wavelet coefficient point to the fitting line, and d max is the maximum distance.λ i is the weighted indicator that indicates that when a wavelet coefficient point has a lager perpendicular distance to the fitting line, it makes a smaller contribution to the global correlation coefficient f and has less ATP signals.Thus, the wavelet coefficients atm w ε j of the ATP signals can be obtained as follows: Following the refined wavelet coefficients, the ATP signals can be reconstructed by IWT.Finally, the corrected D-Infs can be obtained by removing the ATP signals from the original D-Infs.
Another problem should be considered is that how to identify an appropriate number of decomposition scales for the wavelet decomposition to better separate the high-frequency and low-frequency components.This problem can be solved either by visually or based on a priori information.However, owing to the difficulties for either visually or the acquisition of priori information of ATP signals, it is difficult to use the above-mentioned methods to determine an appropriate number of scales for the wavelet decomposition.Thus, in our study a relatively straightforward method was selected, which assumes the ATP signals have a variety of wavelengths from pixel size to broader scales as large as the whole interferogram and was described in [5].
To better understand the MRWCA method, its workflow is presented in Figure 1, which can be divided into the D-InSAR processing, wavelet transform, weight correlation analysis and the estimation of ATP phases.Firstly, two independent external DEMs are adopted to simulate the topographic phases for different polarimetric interferograms so that the obtained D-Infs contain different TEP signals.We then use two-dimensional discrete wavelet transform (2D-DWT) to decompose the D-Infs into different phase components characterized by different wavelengths.However, before doing this, we need to adopt suitable strategy to remove other similar or identical phase components (e.g., orbit error phases).After this, the correlations between the two D-Infs are calculated, which are mainly attributed to the identical ATP signals.Lastly, we can estimate the ATP signals in the D-Infs by IWT and obtain the corrected D-Infs.

Experiments and Analysis
To demonstrate the validity of the proposed approach, we applied it to two interferometric pairs of L-band Advanced Land Observing Satellite (ALOS)-1 PALSAR images covering the San Francisco (USA) and Moron (Mongolia) regions.The first test area is near the Pacific Ocean, while the second test site is located in inland, and they are characterized by different weather and topographic conditions.This results in different atmospheric distributions, allowing us to investigate the effectiveness of the proposed method under different conditions.The images were acquired in fine-beam dual-polarization (FBD) mode, and the spatial resolution is azimuth × range = 3.1 m × 9.4 m.Moreover, since the two test sites are at high geographic latitudes, the total electron content (TEC) is relatively low [18], and the ionospheric effects can be ignored.
For the accuracy assessment for the InSAR DEM, due to the lack of a high-accuracy DEM, the Geoscience Laser Altimeter System (GLAS)/ Ice, Cloud, and land Elevation Satellite (ICESat) L2 Global Land Surface Altimetry Data Version 34 (GLAH 14) [19] are used as the reference data.The GLAS instrument produced laser spots on the Earth's surface with a 70-m diameter.Its vertical accuracy is better than 1 m [20].

San Francisco Test Site
The first test site is located in San Francisco (USA).The dual-polarization (HH and HV polarizations) pair of L-band ALOS-1 PALSAR images have a 46-day interval, and the coverage is roughly 80 km × 75 km, as shown in Figure 2. The test site is close to the Pacific Ocean and features a high water vapor content.

Experiments and Analysis
To demonstrate the validity of the proposed approach, we applied it to two interferometric pairs of L-band Advanced Land Observing Satellite (ALOS)-1 PALSAR images covering the San Francisco (USA) and Moron (Mongolia) regions.The first test area is near the Pacific Ocean, while the second test site is located in inland, and they are characterized by different weather and topographic conditions.This results in different atmospheric distributions, allowing us to investigate the effectiveness of the proposed method under different conditions.The images were acquired in fine-beam dual-polarization (FBD) mode, and the spatial resolution is azimuth × range = 3.1 m × 9.4 m.Moreover, since the two test sites are at high geographic latitudes, the total electron content (TEC) is relatively low [18], and the ionospheric effects can be ignored.
For the accuracy assessment for the InSAR DEM, due to the lack of a high-accuracy DEM, the Geoscience Laser Altimeter System (GLAS)/ Ice, Cloud, and land Elevation Satellite (ICESat) L2 Global Land Surface Altimetry Data Version 34 (GLAH 14) [19] are used as the reference data.The GLAS instrument produced laser spots on the Earth's surface with a 70-m diameter.Its vertical accuracy is better than 1 m [20].

San Francisco Test Site
The first test site is located in San Francisco (USA).The dual-polarization (HH and HV polarizations) pair of L-band ALOS-1 PALSAR images have a 46-day interval, and the coverage is roughly 80 km × 75 km, as shown in Figure 2. The test site is close to the Pacific Ocean and features a high water vapor content.Figure 3a displays the 30-m resolution DEM of the test site acquired by the Shuttle Radar Topography Mission (SRTM) [21], in which the elevation varies between 10 m and 1440 m above mean sea level.Moreover, another 30-m DEM measured by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [22] was also collected and is shown in Figure 3b.The elevation errors of the two DEMs are different and independent, because they were measured by two different techniques characterized by different error sources.The difference map of the two DEMs is shown in Figure 3c, which has a global mean of 0.07 m with a standard deviation (STD) of 7.07 m. (1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 300 m and a temporal baseline of 46 days.The acquisition time interval was only 46 days, so we could exclude the phase caused by topography change in this test site.We firstly adopted two-pass differential interferometry to produce the polarimetric D-Infs.
The interferograms were multi-looked with 10 pixels in the azimuth direction and 2 pixels in the range direction to obtain a final resolution of 30 m × 30 m.The two DEMs (ASTER DEM and SRTM DEM) were then used to remove the topographic phase for the HH-polarization and HV-polarization interferograms, respectively.Moreover, to suppress noise, the interferograms were filtered with the improved Goldstein filter [23], and then the minimum cost flow (MCF) method [24] was applied for the unwrapping process to those pixels that presented coherence of at least 0.3 [25].However, while the interferogram is seriously affected because of decorrelation, which will result in some uncertainty problems in noise filtering and phase unwrapping [26,27].To solve these problems, we should consider the improved methods for D-InSAR processing in the future work.Then, to evaluate the wavelet transform, the data should be continuous on a regular grid, so the Figure 3a displays the 30-m resolution DEM of the test site acquired by the Shuttle Radar Topography Mission (SRTM) [21], in which the elevation varies between 10 m and 1440 m above mean sea level.Moreover, another 30-m DEM measured by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [22] was also collected and is shown in Figure 3b.The elevation errors of the two DEMs are different and independent, because they were measured by two different techniques characterized by different error sources.The difference map of the two DEMs is shown in Figure 3c, which has a global mean of 0.07 m with a standard deviation (STD) of 7.07 m. Figure 3a displays the 30-m resolution DEM of the test site acquired by the Shuttle Radar Topography Mission (SRTM) [21], in which the elevation varies between 10 m and 1440 m above mean sea level.Moreover, another 30-m DEM measured by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [22] was also collected and is shown in Figure 3b.The elevation errors of the two DEMs are different and independent, because they were measured by two different techniques characterized by different error sources.The difference map of the two DEMs is shown in Figure 3c, which has a global mean of 0.07 m with a standard deviation (STD) of 7.07 m. (1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 300 m and a temporal baseline of 46 days.The acquisition time interval was only 46 days, so we could exclude the phase caused by topography change in this test site.We firstly adopted two-pass differential interferometry to produce the polarimetric D-Infs.
The interferograms were multi-looked with 10 pixels in the azimuth direction and 2 pixels in the range direction to obtain a final resolution of 30 m × 30 m.The two DEMs (ASTER DEM and SRTM DEM) were then used to remove the topographic phase for the HH-polarization and HV-polarization interferograms, respectively.Moreover, to suppress noise, the interferograms were filtered with the improved Goldstein filter [23], and then the minimum cost flow (MCF) method [24] was applied for the unwrapping process to those pixels that presented coherence of at least 0.3 [25].However, while the interferogram is seriously affected because of decorrelation, which will result in some uncertainty problems in noise filtering and phase unwrapping [26,27].To solve these problems, we should consider the improved methods for D-InSAR processing in the future work.Then, to evaluate the wavelet transform, the data should be continuous on a regular grid, so the (1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 300 m and a temporal baseline of 46 days.The acquisition time interval was only 46 days, so we could exclude the phase caused by topography change in this test site.We firstly adopted two-pass differential interferometry to produce the polarimetric D-Infs.
The interferograms were multi-looked with 10 pixels in the azimuth direction and 2 pixels in the range direction to obtain a final resolution of 30 m × 30 m.The two DEMs (ASTER DEM and SRTM DEM) were then used to remove the topographic phase for the HH-polarization and HV-polarization interferograms, respectively.Moreover, to suppress noise, the interferograms were filtered with the improved Goldstein filter [23], and then the minimum cost flow (MCF) method [24] was applied for the unwrapping process to those pixels that presented coherence of at least 0.3 [25].However, while the interferogram is seriously affected because of decorrelation, which will result in some uncertainty problems in noise filtering and phase unwrapping [26,27].To solve these problems, we should consider the improved methods for D-InSAR processing in the future work.Then, to evaluate the wavelet transform, the data should be continuous on a regular grid, so the voids with the low-coherence pixels were filled using linear interpolation.Then, after the removal of the orbit error phases by a polynomial model incorporating elevation information [10,15], the obtained D-Infs were superimposed on the SRTM DEM and are shown in Figure 4a  (2) Estimation and correction of the ATP: To begin with, we decomposed the obtained dual-polarization D-Infs into different scales.Before doing this, we needed to choose an optimal wavelet decomposition scale J , which could ensure the separation of the ATP signals and the other non-ATP signals.In this study, the optimal decomposition scale was set to 11, for which we refer the reader to [5].We could then establish the linear relationship between the wavelet coefficients for each level of decomposition with Equation (3).To better understand the correlation analysis, the red solid lines in Figure 5 show scatter plots of the wavelet coefficients in the horizontal direction.It is clear to see that the correlation of the wavelet coefficients increases with the decomposition scale.The high correlation indicates that the wavelet coefficients are almost all contributed by the same ATP signals.In addition, the scatter plots have similar trends in the vertical and diagonal directions.We can see that the correlation between the coefficients expresses the global trend; however, the ratios of atmosphere composition are different for each pixel, so the values of the ATP ratio should also have some differences.In this study, we append an appropriate weighting indicator to refine the correlation coefficients, which can be used to extract the ATP signals accurately.Thus, the (2) Estimation and correction of the ATP: To begin with, we decomposed the obtained dual-polarization D-Infs into different scales.Before doing this, we needed to choose an optimal wavelet decomposition scale J, which could ensure the separation of the ATP signals and the other non-ATP signals.In this study, the optimal decomposition scale was set to 11, for which we refer the reader to [5].We could then establish the linear relationship between the wavelet coefficients for each level of decomposition with Equation (3).To better understand the correlation analysis, the red solid lines in Figure 5 show scatter plots of the wavelet coefficients in the horizontal direction.It is clear to see that the correlation of the wavelet coefficients increases with the decomposition scale.The high correlation indicates that the wavelet coefficients are almost all contributed by the same ATP signals.
In addition, the scatter plots have similar trends in the vertical and diagonal directions.(2) Estimation and correction of the ATP: To begin with, we decomposed the obtained dual-polarization D-Infs into different scales.Before doing this, we needed to choose an optimal wavelet decomposition scale J , which could ensure the separation of the ATP signals and the other non-ATP signals.In this study, the optimal decomposition scale was set to 11, for which we refer the reader to [5].We could then establish the linear relationship between the wavelet coefficients for each level of decomposition with Equation (3).To better understand the correlation analysis, the red solid lines in Figure 5 show scatter plots of the wavelet coefficients in the horizontal direction.It is clear to see that the correlation of the wavelet coefficients increases with the decomposition scale.The high correlation indicates that the wavelet coefficients are almost all contributed by the same ATP signals.In addition, the scatter plots have similar trends in the vertical and diagonal directions.We can see that the correlation between the coefficients expresses the global trend; however, the ratios of atmosphere composition are different for each pixel, so the values of the ATP ratio should also have some differences.In this study, we append an appropriate weighting indicator to refine the correlation coefficients, which can be used to extract the ATP signals accurately.Thus, the We can see that the correlation between the coefficients expresses the global trend; however, the ratios of atmosphere composition are different for each pixel, so the values of the ATP ratio should also have some differences.In this study, we append an appropriate weighting indicator to refine the correlation coefficients, which can be used to extract the ATP signals accurately.Thus, the wavelet coefficients of ATP signals can be calculated by Equation (5), and the ATP signals can be obtained by computing the IWT using the refined wavelet coefficients.
The final corrected HH and HV D-Infs are presented in Figure 6a,b, respectively.Clearly, most of the positive and negative signals in Figure 4a,b are removed by the use of the MRWCA method.Figure 6c shows the difference map of the corrected HV and HH D-Infs.(3) Assessment of the derived DEMs: After mitigating the ATP signals, the corrected HH interferogram was used to generate the DEM since it has higher phase quality, which is of benefit to extract DEM.To verify the effectiveness of the proposed method, we evaluated the vertical accuracy of the InSAR DEMs that were generated without atmospheric effects correction, with atmospheric effects correction using Shirzaei's atmospheric correction method [5], and with atmospheric effects correction using the MRWCA method, as shown in Figure 7. Firstly, owing to the lack of actual ground measurement elevation data, while the ASTER DEM was used as input data for the HH-polarization interferogram, we selected the SRTM DEM covering the experimental region as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the elevation difference distribution, as shown in Figure 8. From Figure 8, we can see the elevation difference without atmospheric correction is significantly larger than that after the atmosphere correction process.In addition, the InSAR DEM with Shirzaei's atmospheric correction method also contains large errors when compared with the result of the MRWCA method, which can be attributed to the fact that Shirzaei's method is used to evaluate the ATP signals that are correlated with topography.However, the MRWCA method can remove both topography-dependent ATP signals (TD-ATP) and the turbulent ATP (T-ATP) signals [12], without depending on the relationship between ATP signals and topography.(3) Assessment of the derived DEMs: After mitigating the ATP signals, the corrected HH interferogram was used to generate the DEM since it has higher phase quality, which is of benefit to extract DEM.To verify the effectiveness of the proposed method, we evaluated the vertical accuracy of the InSAR DEMs that were generated without atmospheric effects correction, with atmospheric effects correction using Shirzaei's atmospheric correction method [5], and with atmospheric effects correction using the MRWCA method, as shown in Figure 7.   (3) Assessment of the derived DEMs: After mitigating the ATP signals, the corrected HH interferogram was used to generate the DEM since it has higher phase quality, which is of benefit to extract DEM.To verify the effectiveness of the proposed method, we evaluated the vertical accuracy of the InSAR DEMs that were generated without atmospheric effects correction, with atmospheric effects correction using Shirzaei's atmospheric correction method [5], and with atmospheric effects correction using the MRWCA method, as shown in Figure 7. Firstly, owing to the lack of actual ground measurement elevation data, while the ASTER DEM was used as input data for the HH-polarization interferogram, we selected the SRTM DEM covering the experimental region as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the elevation difference distribution, as shown in Figure 8. From Figure 8, we can see the elevation difference without atmospheric correction is significantly larger than that after the atmosphere correction process.In addition, the InSAR DEM with Shirzaei's atmospheric correction method also contains large errors when compared with the result of the MRWCA method, which can be attributed to the fact that Shirzaei's method is used to evaluate the ATP signals that are correlated with topography.However, the MRWCA method can remove both topography-dependent ATP signals (TD-ATP) and the turbulent ATP (T-ATP) signals [12], without depending on the relationship between ATP signals and topography.Firstly, owing to the lack of actual ground measurement elevation data, while the ASTER DEM was used as input data for the HH-polarization interferogram, we selected the SRTM DEM covering the experimental region as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the elevation difference distribution, as shown in Figure 8. From Figure 8, we can see the elevation difference without atmospheric correction is significantly larger than that after the atmosphere correction process.In addition, the InSAR DEM with Shirzaei's atmospheric correction method also contains large errors when compared with the result of the MRWCA method, which can be attributed to the fact that Shirzaei's method is used to evaluate the ATP signals that are correlated with topography.However, the MRWCA method can remove both topography-dependent ATP signals (TD-ATP) and the turbulent ATP (T-ATP) signals [12], without depending on the relationship between ATP signals and topography.In order to further verify the accuracy, we used ICESat altimetry data as a reference for the quantitative analysis.The GLAS instrument produces laser spots on the Earth's surface with a 70-m diameter in each footprint location of ICESat, so we calculated an average value of the DEMs within 70 m to evaluate the difference between the ICESat and InSAR-derived DEMs.To quantify the difference, Figure 9 shows the histograms of the differences.Quantitatively, we evaluate the absolute vertical accuracy of the InSAR-derived DEMs in Figure 7, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 13.10 m to 10.16 m, while the RMSE has decreased to 7.79 m with the MRWCA correction method.

Moron Test Site
The second test area covers approximately 80 km × 75 km of the Moron area, located in the northwestern part of Mongolia.This area has a landform of mountains with steep slopes, where the elevation changes substantially, and features complex terrain conditions.Two L-band ALOS-1 PALSAR images acquired in dual-polarization (HH and HV) were collected, as shown in Figure 10.In order to further verify the accuracy, we used ICESat altimetry data as a reference for the quantitative analysis.The GLAS instrument produces laser spots on the Earth's surface with a 70-m diameter in each footprint location of ICESat, so we calculated an average value of the DEMs within 70 m to evaluate the difference between the ICESat and InSAR-derived DEMs.To quantify the difference, Figure 9 shows the histograms of the differences.Quantitatively, we evaluate the absolute vertical accuracy of the InSAR-derived DEMs in Figure 7, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 13.10 m to 10.16 m, while the RMSE has decreased to 7.79 m with the MRWCA correction method.In order to further verify the accuracy, we used ICESat altimetry data as a reference for the quantitative analysis.The GLAS instrument produces laser spots on the Earth's surface with a 70-m diameter in each footprint location of ICESat, so we calculated an average value of the DEMs within 70 m to evaluate the difference between the ICESat and InSAR-derived DEMs.To quantify the difference, Figure 9 shows the histograms of the differences.Quantitatively, we evaluate the absolute vertical accuracy of the InSAR-derived DEMs in Figure 7, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 13.10 m to 10.16 m, while the RMSE has decreased to 7.79 m with the MRWCA correction method.

Moron Test Site
The second test area covers approximately 80 km × 75 km of the Moron area, located in the northwestern part of Mongolia.This area has a landform of mountains with steep slopes, where the elevation changes substantially, and features complex terrain conditions.Two L-band ALOS-1 PALSAR images acquired in dual-polarization (HH and HV) were collected, as shown in Figure 10.

Moron Test Site
The second test area covers approximately 80 km × 75 km of the Moron area, located in the northwestern part of Mongolia.This area has a landform of mountains with steep slopes, where the elevation changes substantially, and features complex terrain conditions.Two L-band ALOS-1 PALSAR images acquired in dual-polarization (HH and HV) were collected, as shown in Figure 10.(1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 190 m and a temporal baseline of 46 days.We then applied the same two-pass differential interferometry strategy to obtain the two polarimetric D-Infs, which were superimposed on the SRTM DEM and are shown in Figure 12a    (1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 190 m and a temporal baseline of 46 days.We then applied the same two-pass differential interferometry strategy to obtain the two polarimetric D-Infs, which were superimposed on the SRTM DEM and are shown in Figure 12a   (1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 190 m and a temporal baseline of 46 days.We then applied the same two-pass differential interferometry strategy to obtain the two polarimetric D-Infs, which were superimposed on the SRTM DEM and are shown in Figure 12a,b  (1) Generation of differential interferograms: The dual-polarization SAR images were used to constitute two polarimetric interferograms with a perpendicular baseline [12] of 190 m and a temporal baseline of 46 days.We then applied the same two-pass differential interferometry strategy to obtain the two polarimetric D-Infs, which were superimposed on the SRTM DEM and are shown in Figure 12a,b   (2) Estimation and correction of the ATP: To begin with, we also decompose the obtained dual-polarization D-Infs into different scales.For this test site, the optimal decomposition scale was set as 11.We could then establish the linear relationship between the high-frequency wavelet coefficients for each level of decomposition with Equation (3).The red solid lines in Figure 13 show scatter plots of the wavelet coefficients in the vertical direction, where it can be seen that the scatter plots have similar trends in the horizontal and diagonal directions.
(2) Estimation and correction of the ATP: To begin with, we also decompose the obtained dual-polarization D-Infs into different scales.For this test site, the optimal decomposition scale was set as 11.We could then establish the linear relationship between the high-frequency wavelet coefficients for each level of decomposition with Equation (3).The red solid lines in Figure 13 show scatter plots of the wavelet coefficients in the vertical direction, where it can be seen that the scatter plots have similar trends in the horizontal and diagonal directions.(3) Assessment of the derived DEMs: Lastly, after eliminating the effects of ATP signals, the corrected HH interferogram was used to generate the InSAR DEMs.We evaluated the vertical accuracy of the InSAR DEMs extracted without atmospheric effects correction, with atmospheric effects correction using Shirzaei's correction method [5], and with atmospheric effects correction using the MRWCA method, as shown in Figure 15.(2) Estimation and correction of the ATP: To begin with, we also decompose the obtained dual-polarization D-Infs into different scales.For this test site, the optimal decomposition scale was set as 11.We could then establish the linear relationship between the high-frequency wavelet coefficients for each level of decomposition with Equation (3).The red solid lines in Figure 13 show scatter plots of the wavelet coefficients in the vertical direction, where it can be seen that the scatter plots have similar trends in the horizontal and diagonal directions.(3) Assessment of the derived DEMs: Lastly, after eliminating the effects of ATP signals, the corrected HH interferogram was used to generate the InSAR DEMs.We evaluated the vertical accuracy of the InSAR DEMs extracted without atmospheric effects correction, with atmospheric effects correction using Shirzaei's correction method [5], and with atmospheric effects correction using the MRWCA method, as shown in Figure 15.(3) Assessment of the derived DEMs: Lastly, after eliminating the effects of ATP signals, the corrected HH interferogram was used to generate the InSAR DEMs.We evaluated the vertical accuracy of the InSAR DEMs extracted without atmospheric effects correction, with atmospheric effects correction using Shirzaei's correction method [5], and with atmospheric effects correction using the MRWCA method, as shown in Figure 15.We also used the SRTM DEM as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the elevation difference distribution, as shown in Figure 16.We can also see that the elevation difference without atmospheric correction is significantly larger than that after the correction process, and the MRWCA correction method shows the best elevation accuracy.We then used the ICESat altimetry data as the reference for the quantitative analysis, and the elevation differences between the InSAR-derived DEMs and the ICESat data were analyzed.The histograms of the differences are displayed in Figure 17.Quantitatively, we evaluate the vertical accuracy of the InSAR-derived DEMs in Figure 15, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 15.40 m to 12 m, and after MRWCA atmospheric correction, the RMSE has decreased to 10.74 m.We also used the SRTM DEM as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the difference distribution, as shown in Figure 16.We can also see that the elevation difference without atmospheric correction is significantly larger than that after the correction process, and the MRWCA correction method shows the best elevation accuracy.We also used the SRTM DEM as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the elevation difference distribution, as shown in Figure 16.We can also see that the elevation difference without atmospheric correction is significantly larger than that after the correction process, and the MRWCA correction method shows the best elevation accuracy.We then used the ICESat altimetry data as the reference for the quantitative analysis, and the elevation differences between the InSAR-derived DEMs and the ICESat data were analyzed.The histograms of the differences are displayed in Figure 17.Quantitatively, we evaluate the vertical accuracy of the InSAR-derived DEMs in Figure 15, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 15.40 m to 12 m, and after MRWCA atmospheric correction, the RMSE has decreased to 10.74 m.We then used the ICESat altimetry data as the reference for the quantitative analysis, and the elevation differences between the InSAR-derived DEMs and the ICESat data were analyzed.The histograms of the differences are displayed in Figure 17.Quantitatively, we evaluate the vertical accuracy of the InSAR-derived DEMs in Figure 15, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 15.40 m to 12 m, and after MRWCA atmospheric correction, the RMSE has decreased to 10.74 m.We also used the SRTM DEM as the validation data for the qualitative analysis.We calculated the vertical difference of the InSAR-derived DEMs with the validation SRTM DEM, and obtained the elevation difference distribution, as shown in Figure 16.We can also see that the elevation difference without atmospheric correction is significantly larger than that after the correction process, and the MRWCA correction method shows the best elevation accuracy.We then used the ICESat altimetry data as the reference for the quantitative analysis, and the elevation differences between the InSAR-derived DEMs and the ICESat data were analyzed.The histograms of the differences are displayed in Figure 17.Quantitatively, we evaluate the vertical accuracy of the InSAR-derived DEMs in Figure 15, which shows that after the use of Shirzaei's atmospheric correction method, the RMSE has decreased from 15.40 m to 12 m, and after MRWCA atmospheric correction, the RMSE has decreased to 10.74 m.

Summary of the MRWCA Method
This paper has presented an effective atmosphere correction method for InSAR topographic mapping, which integrates wavelet multi-resolution analysis and weighted correlation analysis.The main idea of the MRWCA method is based on a priori knowledge that the dual-polarization interferograms have identical ATP signals.To detect the same ATP signals, we firstly apply two-pass differential interferometry technology to generate D-Infs and remove any other similar or identical phase components, the only common component in the two obtained D-Infs is the ATP signals.Thus, the obtained wavelet coefficients have the same or similar values over a range of scales, which represents the key rationale of correlation analysis.Following the decomposition of different polarimetric D-Infs, the global correlation coefficient between the corresponding wavelet coefficients are most likely to be the ratio of the ATP signals in the obtained D-Infs.To improve the sensitivity of the global correlation coefficient to the ATP variations, a weighted indicator is formulated to refine the global correlation coefficient.In most applications of wavelet transforms, the choice of wavelet family is critical.However, in this paper, this choice is rather straightforward.We choose a wavelet that can separate the non-ATP signals (e.g., TEP and NP) from the D-Inf [28].In other words, the selected wavelet should be able to better separate the ATP and non-ATP signals into different phase components characterized by different wavelengths, which will reduce the effects of confounding signals (e.g.non-ATP signals) and improve the robustness of the correlation analysis.One should note that continuous wavelets may have better capacity to support the correlation analysis since they have higher spectral resolution, which will be investigated in the future works.

Self-Consistency Checking between Differences Maps before and after MRWCA Method Correction
In order to quantitatively check the phase quality after MRWCA method correction, the correlation of the differences between the HV and HH D-Infs before and after atmospheric effects correction were investigated, the joint distributions were shown in Figure 18a,b.We can see that, in the San Francisco experiment, the corrected and original difference maps, as shown in Figures 6c and 4c, have a correlation coefficient of 0.94 and an RMSE of 0.07 rad.In the Moron test site, a correlation coefficient of 0.84 with an RMSE of 0.05 rad was obtained, as presented in Figure 18b.In both of the cases, the high correlations and the small RMSEs indicate that the MRWCA method can correct ATP signals, without significantly reducing the phase quality.With the consideration of phase-to-height conversion parameters associated with the two interferometric pairs, these phase losses cannot be the major factor mitigating the DEM quality.However, for both of the cases, the corresponding correlation coefficients less than 1 indicate that some of the non-ATP signals are regarded as ATP signals by the MRWCA method.This can be explained by two possible reasons: (1) the global correlation analysis cannot accurately describe the correlations of some pixels which attributed to ATP signals.As a result, the ATP signals can be over-or under-estimated.(2) If the phase centers of different polarimetric interferograms are similar, some of the non-ATP signals can also be regarded as the ATP signals.

Effects of Different Decomposition Scales on ATP Correction
To better understand the effect of decomposition scale on ATP correction, as an example, we chose four different decomposition scales (3,6,9,11) for ATP correction for Moron test site as shown in Figure 19.We firstly used the SRTM DEM as the validation data for the qualitative analysis, the elevation difference distribution maps between InSAR DEMs and SRTM DEM are shown in Figure 19.Clearly, although a part of ATP signals can be removed by using the small decomposition scales, the residual error attributed to ATP signals are still significant.In order to investigate the corrected effects more clearly, the ICESat altimetry data was used for the quantitative analysis, the corresponding RMSEs of elevation difference are 15 m, 14.84 m, 12.43 m and 10.74 m, respectively.
Based on the above analysis, we can find that, the smaller decomposition scales can only correct the ATP signals with short wavelength, however, the ATP signals of the Moron test site are dominated by longer wavelength that can be confirmed by Figure 13.Therefore, it is necessary to select an appropriate number of decomposition scales for the wavelet decomposition which can separate almost all the ATP signals mixed with different non-ATP signals from the D-Inf.

Comparison with the Existing Method
We extracted the ATP maps that are shown in Figure 20a,d for the two experiment regions.For the two HH D-Infs, the ATP signals make up 63.6% and 44.8% of the total D-Infs signals.For the San Francisco test site, the corrected interferogram could provide a DEM with a RMSE of 7.79 m (Figure

Effects of Different Decomposition Scales on ATP Correction
To better understand the effect of decomposition scale on ATP correction, as an example, we chose four different decomposition scales (3,6,9,11) for ATP correction for Moron test site as shown in Figure 19.

Effects of Different Decomposition Scales on ATP Correction
To better understand the effect of decomposition scale on ATP correction, as an example, we chose four different decomposition scales (3,6,9,11) for ATP correction for Moron test site as shown in Figure 19.We firstly used the SRTM DEM as the validation data for the qualitative analysis, the elevation difference distribution maps between InSAR DEMs and SRTM DEM are shown in Figure 19.Clearly, although a part of ATP signals can be removed by using the small decomposition scales, the residual error attributed to ATP signals are still significant.In order to investigate the corrected effects more clearly, the ICESat altimetry data was used for the quantitative analysis, the corresponding RMSEs of elevation difference are 15 m, 14.84 m, 12.43 m and 10.74 m, respectively.
Based on the above analysis, we can find that, the smaller decomposition scales can only correct the ATP signals with short wavelength, however, the ATP signals of the Moron test site are dominated by longer wavelength that can be confirmed by Figure 13.Therefore, it is necessary to select an appropriate number of decomposition scales for the wavelet decomposition which can separate almost all the ATP signals mixed with different non-ATP signals from the D-Inf.

Comparison with the Existing Method
We extracted the ATP maps that are shown in Figure 20a,d   We firstly used the SRTM DEM as the validation data for the qualitative analysis, the elevation difference distribution maps between InSAR DEMs and SRTM DEM are shown in Figure 19.Clearly, although a part of ATP signals can be removed by using the small decomposition scales, the residual error attributed to ATP signals are still significant.In order to investigate the corrected effects more clearly, the ICESat altimetry data was used for the quantitative analysis, the corresponding RMSEs of elevation difference are 15 m, 14.84 m, 12.43 m and 10.74 m, respectively.
Based on the above analysis, we can find that, the smaller decomposition scales can only correct the ATP signals with short wavelength, however, the ATP signals of the Moron test site are dominated by longer wavelength that can be confirmed by Figure 13.Therefore, it is necessary to select an appropriate number of decomposition scales for the wavelet decomposition which can separate almost all the ATP signals mixed with different non-ATP signals from the D-Inf.

Comparison with the Existing Method
We extracted the ATP maps that are shown in Figure 20a,d for the two experiment regions.For the two HH D-Infs, the ATP signals make up 63.6% and 44.8% of the total D-Infs signals.For the San Francisco test site, the corrected interferogram could provide a DEM with a RMSE of 7.79 m (Figure 9), which is an improvement of 40.5% with respect to the DEM before atmospheric effects correction; for the Moron test site, the corrected interferogram could provide a DEM with an RMSE of 10.74 m, which is an improvement of 30.2% with respect to the DEM without atmospheric effects correction (Figure 17), which demonstrates the MRWCA method is effective in detecting the ATP signals with different ATP content conditions.
From the extracted ATP maps (Figure 20a,d), we can also find that the San Francisco experiment has serious atmospheric effects, which are attributed to this test site because it is close to the Pacific Ocean and features a high water vapor content.The effect of ATP signals can be divided into TD-ATP and T-ATP [12].Although the polynomial model incorporating elevation information can remove the TD-ATP signals in the orbit error removing [10,15], only the ATP signals linearly correlated with topography can be deal estimated.In order to solve this problem, Shirzaei proposed a wavelet decomposition based method to remove the TD-ATP signals [5].With Shirzaei's atmospheric correction method, the extracted TD-ATP maps are shown in Figure 20b,e From the extracted ATP maps (Figure 20a,d), we can also find that the San Francisco experiment has serious atmospheric effects, which are attributed to this test site because it is close to the Pacific Ocean and features a high water vapor content.The effect of ATP signals can be divided into TD-ATP and T-ATP [12].Although the polynomial model incorporating elevation information can remove the TD-ATP signals in the orbit error removing [10,15], only the ATP signals linearly correlated with topography can be deal estimated.In order to solve this problem, Shirzaei proposed a wavelet decomposition based method to remove the TD-ATP signals [5].With Shirzaei's atmospheric correction method, the extracted TD-ATP maps are shown in Figure 20b,e

Potential and Limits of MRWCA Method in Atmospheric Correction
Using the proposed method, which does not depend on any assumption or external meteorological data, we only need single-baseline dual-polarization SAR data and two different external DEMs, thus it is relatively easier to collect and process these data.So we believe this study will provide an effective method to correct the atmospheric effects by using single-baseline SAR data (e.g., Sentinel-1 [29], ALOS-2 [30], and future Tandem-L [31]), which contain multi-polarization SAR data, have high resolution, and a short revisit period.As a result, we can use the MRWCA method to correct ATP signals for InSAR DEM generation and invert high resolution water vapor data.In this study, the selected interferometric pairs had not underwent a significant ionospheric effect, which can be divided into a Faraday angle and ionospheric delay.Although the ionosphere induced Faraday angle is polarization-dependent [32], the ionospheric delay is identical for all polarizations [33].In such a case, we can also remove the ionospheric phase from the D-Inf by the MRWCA method and have no impact on DEM inversion.But when the D-Infs contain large proportion ionosphere components, which will make the polarimetric D-Infs more similar, in particular for

Potential and Limits of MRWCA Method in Atmospheric Correction
Using the proposed method, which does not depend on any assumption or external meteorological data, we only need single-baseline dual-polarization SAR data and two different external DEMs, thus it is relatively easier to collect and process these data.So we believe this study will provide an effective method to correct the atmospheric effects by using single-baseline SAR data (e.g., Sentinel-1 [29], ALOS-2 [30], and future Tandem-L [31]), which contain multi-polarization SAR data, have high resolution, and a short revisit period.As a result, we can use the MRWCA method to correct ATP signals for InSAR DEM generation and invert high resolution water vapor data.In this study, the selected interferometric pairs had not underwent a significant ionospheric effect, which can be divided into a Faraday angle and ionospheric delay.Although the ionosphere induced Faraday angle is polarization-dependent [32], the ionospheric delay is identical for all polarizations [33].In such a case, we can also remove the ionospheric phase from the D-Inf by the MRWCA method and have no impact on DEM inversion.But when the D-Infs contain large proportion ionosphere components, which will make the polarimetric D-Infs more similar, in particular for interferograms acquired at low-frequency, and make it difficult to estimate the ionospheric and atmospheric phases without losing the TEP signals by the correlation analysis.To solve this problem, before performing the MRWCA method, the ionospheric phases should be removed by the proper methods [34,35].However, after correcting the ionospheric effect, the interferograms still contain residual ionospheric phases [36,37], and the effect of residual ionospheric phases on the correlation analysis is dependent on how much of the ionospheric phase is left.In addition, the effect of residual ionospheric phases on the correlation analysis still needs to be further investigated.
The major limit of this method is that the unwrapped phase is required.If the unwrapping step fails due to very low coherence, this can affect the application of the proposed approach.In this study, although the phase unwrapping method can affect the atmospheric phase removal, it is not the major limited factor.Thus, we have not done further research about the effect of phase unwrapping on the atmospheric phase removal.Moreover, the proposed approach is not to distinguish the atmospheric effects and the ground surface deformation signals which are similar or same in the interferograms acquired in different polarizations.Therefore, if the deformation signals are identical in different polarimetric D-Infs, the MRWCA method can remove the deformation signals during the ATP removal and have no impact on DEM production.However, when the polarimetric D-Infs have similar deformation signals, some of the deformation signals will be miss-interpreted as TEP signals, which will be converted into the wrong topographic information.

Conclusions
In this paper, a novel approach to correct the atmospheric effects in a repeat-pass interferogram has been presented, which combines wavelet multi-resolution analysis and weighted correlation analysis to identify the same ATP signals from dual-polarization InSAR interferograms.We successfully tested the method on two test sites over San Francisco (USA) and Moron (Mongolia), which have different topographic characteristics and different water vapor distributions.The results demonstrate that the atmospheric effects can be effectively eliminated.The estimated DEMs have improvements of 40.5% and 30.2% in term of RMSE for the two test sites, respectively.In addition, compared with the existing method [5], the proposed method cannot only deal with the TD-ATP signals, but also the T-ATP signals.The proposed method does not depend on any assumption or external meteorological data and we only need single-baseline dual-polarization SAR data and two different external DEMs, which allow us to extract the DEM with large-scale and high-resolution.
In our future work, we will further test the application scope of the method by selecting SAR data acquired from different topographic conditions.In addition, we will investigate how to use full-polarization data, which contain more scattering information being helpful to conduct the correlation analysis, to identify more accurate ATP signals.Furthermore, how to use these ATP signals to estimate high resolution water vapor data may also help us to better understand the atmosphere.Lastly, in this study, the effects of the ground deformation and ionospheric delay have been neglected, how to reduce these effects on DEM extraction needs to be further investigated.

Figure 1 .
Figure 1.The workflow of the multi-resolution weighted correlation analysis (MRWCA) method.

Figure 1 .
Figure 1.The workflow of the multi-resolution weighted correlation analysis (MRWCA) method.

Figure 2 .
Figure 2. The footprint of the ALOS-1 PALSAR images over the San Francisco region and the spatial coverage of the Ice, Cloud, and land Elevation Satellite/Geoscience Laser Altimeter System (ICESat/GLAS) data over the study area.

Figure 3 .
Figure 3. (a) Shuttle Radar Topography Mission digital elevation model (SRTM DEM); (b) Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM; (c) Difference map between the SRTM DEM and the ASTER DEM.

Figure 2 .
Figure 2. The footprint of the ALOS-1 PALSAR images over the San Francisco region and the spatial coverage of the Ice, Cloud, and land Elevation Satellite/Geoscience Laser Altimeter System (ICESat/GLAS) data over the study area.

Figure 2 .
Figure 2. The footprint of the ALOS-1 PALSAR images over the San Francisco region and the spatial coverage of the Ice, Cloud, and land Elevation Satellite/Geoscience Laser Altimeter System (ICESat/GLAS) data over the study area.

Figure 3 .
Figure 3. (a) Shuttle Radar Topography Mission digital elevation model (SRTM DEM); (b) Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM; (c) Difference map between the SRTM DEM and the ASTER DEM.

Figure 3 .
Figure 3. (a) Shuttle Radar Topography Mission digital elevation model (SRTM DEM); (b) Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM; (c) Difference map between the SRTM DEM and the ASTER DEM.
,b.It is clear that the obtained D-Infs contain obvious positive and negative signals, which are contributed by the identical ATP signals.In addition, the difference map of the HV and HH D-Infs is shown in Figure 4c.Since the same ATP signals have been counteracted, the difference map is made up of the topography-relevant phase signals caused by the difference of the external elevation data and the difference of the polarimetric phase centers.Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 18 voids with the low-coherence pixels were filled using linear interpolation.Then, after the removal of the orbit error phases by a polynomial model incorporating elevation information [10,15], the obtained D-Infs were superimposed on the SRTM DEM and are shown in Figure 4a,b.It is clear that the obtained D-Infs contain obvious positive and negative signals, which are contributed by the identical ATP signals.In addition, the difference map of the HV and HH D-Infs is shown in Figure 4c.Since the same ATP signals have been counteracted, the difference map is made up of the topography-relevant phase signals caused by the difference of the external elevation data and the difference of the polarimetric phase centers.

Figure 4 .
Figure 4. Original (a) HH and (b) HV D-Infs; (c) the difference map between the HV and HH D-Infs.

Figure 5 .
Figure 5. Scatter plots (HV vs. HH) of the wavelet coefficients in horizontal direction (a-k indicates the wavelet decomposition scales ranging from 1 to 11).

Figure 4 .
Figure 4. Original (a) HH and (b) HV D-Infs; (c) the difference map between the HV and HH D-Infs.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 18 voids with the low-coherence pixels were filled using linear interpolation.Then, after the removal of the orbit error phases by a polynomial model incorporating elevation information [10,15], the obtained D-Infs were superimposed on the SRTM DEM and are shown in Figure 4a,b.It is clear that the obtained D-Infs contain obvious positive and negative signals, which are contributed by the identical ATP signals.In addition, the difference map of the HV and HH D-Infs is shown in Figure 4c.Since the same ATP signals have been counteracted, the difference map is made up of the topography-relevant phase signals caused by the difference of the external elevation data and the difference of the polarimetric phase centers.

Figure 4 .
Figure 4. Original (a) HH and (b) HV D-Infs; (c) the difference map between the HV and HH D-Infs.

Figure 5 .
Figure 5. Scatter plots (HV vs. HH) of the wavelet coefficients in horizontal direction (a-k indicates the wavelet decomposition scales ranging from 1 to 11).

Figure 5 .
Figure 5. Scatter plots (HV vs. HH) of the wavelet coefficients in horizontal direction (a-k indicates the wavelet decomposition scales ranging from 1 to 11).
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 18 wavelet coefficients of ATP signals can be calculated by Equation (5), and the ATP signals can be obtained by computing the IWT using the refined wavelet coefficients.The final corrected HH and HV D-Infs are presented in Figure 6a,b, respectively.Clearly, most of the positive and negative signals in Figure 4a,b are removed by the use of the MRWCA method.
Figure 6c shows the difference map of the corrected HV and HH D-Infs.

Figure 6 .
Figure 6.(a) HH and (b) HV D-Infs corrected by the MRWCA method.(c) The difference map between the corrected HV and HH D-Infs.

Figure 6 .
Figure 6.(a) HH and (b) HV D-Infs corrected by the MRWCA method.(c) The difference map between the corrected HV and HH D-Infs.
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 18 wavelet coefficients of ATP signals can be calculated by Equation (5), and the ATP signals can be obtained by computing the IWT using the refined wavelet coefficients.The final corrected HH and HV D-Infs are presented in Figure 6a,b, respectively.Clearly, most of the positive and negative signals in Figure 4a,b are removed by the use of the MRWCA method.
Figure 6c shows the difference map of the corrected HV and HH D-Infs.

Figure 6 .
Figure 6.(a) HH and (b) HV D-Infs corrected by the MRWCA method.(c) The difference map between the corrected HV and HH D-Infs.

Figure 9 .
Figure 9. Histograms of the elevation differences between the ICESat data and the InSAR-derived DEMs generated without atmospheric correction, with Shirzaei's atmospheric correction and MRWCA atmospheric correction, respectively.

Figure 8 .
Figure 8.The elevation difference distribution (a) without atmospheric correction; (b) with Shirzaei's atmospheric correction method; and (c) with MRWCA atmospheric correction.

Figure 9 .
Figure 9. Histograms of the elevation differences between the ICESat data and the InSAR-derived DEMs generated without atmospheric correction, with Shirzaei's atmospheric correction and MRWCA atmospheric correction, respectively.

Figure 9 .
Figure 9. Histograms of the elevation differences between the ICESat data and the InSAR-derived DEMs generated without atmospheric correction, with Shirzaei's atmospheric correction and MRWCA atmospheric correction, respectively.

Figure 10 .
Figure 10.The footprint of the ALOS-1 PALSAR images over the Moron region and the spatial coverage of the ICESat/GLAS data over the study area.The reference DEMs were the SRTM DEM and the ASTER DEM with a spatial resolution of 30 m, which are shown in Figure 11a,b.The elevation in this area varies from 1220 m to 2560 m above mean sea level.The difference map of the two DEMs is shown in Figure 11c, which has a global mean of 0.05 m with a standard deviation (STD) of 10.71 m.
,b.In addition, the difference map is shown in Figure 12c.Clearly, the obvious phase signals in Figure 12a,b have disappeared, which confirms that the original HH and HV-polarization D-Infs contain common ATP signals.Moreover, the time interval of the InSAR pairs was 46 days, so we could also exclude surface deformation, so the D-Infs could be considered to be composed of ATP signals, TEP signals, and NP signals.

Figure 12 .
Figure 12.Original (a) HH and (b) HV D-Infs.(c) The difference map between the HV and HH D-Infs.

Figure 10 . 18 Figure 10 .
Figure 10.The footprint of the ALOS-1 PALSAR images over the Moron region and the spatial coverage of the ICESat/GLAS data over the study area.
,b.In addition, the difference map is shown in Figure 12c.Clearly, the obvious phase signals in Figure 12a,b have disappeared, which confirms that the original HH and HV-polarization D-Infs contain common ATP signals.Moreover, the time interval of the InSAR pairs was 46 days, so we could also exclude surface deformation, so the D-Infs could be considered to be composed of ATP signals, TEP signals, and NP signals.

Figure 12 .
Figure 12.Original (a) HH and (b) HV D-Infs.(c) The difference map between the HV and HH D-Infs.

18 Figure 10 .
Figure 10.The footprint of the ALOS-1 PALSAR images over the Moron region and the spatial coverage of the ICESat/GLAS data over the study area.The reference DEMs were the SRTM DEM and the ASTER DEM with a spatial resolution of 30 m, which are shown in Figure 11a,b.The elevation in this area varies from 1220 m to 2560 m above mean sea level.The difference map of the two DEMs is shown in Figure 11c, which has a global mean of 0.05 m with a standard deviation (STD) of 10.71 m.
. In addition, the difference map is shown in Figure 12c.Clearly, the obvious phase signals in Figure 12a,b have disappeared, which confirms that the original HH and HV-polarization D-Infs contain common ATP signals.Moreover, the time interval of the InSAR pairs was 46 days, so we could also exclude surface deformation, so the D-Infs could be considered to be composed of ATP signals, TEP signals, and NP signals.

Figure 12 .
Figure 12.Original (a) HH and (b) HV D-Infs.(c) The difference map between the HV and HH D-Infs.Figure 12. Original (a) HH and (b) HV D-Infs.(c) The difference map between the HV and HH D-Infs.

Figure 12 .
Figure 12.Original (a) HH and (b) HV D-Infs.(c) The difference map between the HV and HH D-Infs.Figure 12. Original (a) HH and (b) HV D-Infs.(c) The difference map between the HV and HH D-Infs.

Figure 13 .
Figure 13.Scatter plots (HV vs. HH) of the wavelet coefficients in the vertical direction (a-k indicates the wavelet decomposition scales ranging from 1 to 11).

Figure 14 .
Figure 14.(a) HH and (b) HV D-Infs corrected by the MRWCA method.(c) The difference map between the corrected HV and HH D-Infs.

Figure 13 .
Figure 13.Scatter plots (HV vs. HH) of the wavelet coefficients in the vertical direction (a-k indicates the wavelet decomposition scales ranging from 1 to 11).

Figure 13 .
Figure 13.Scatter plots (HV vs. HH) of the wavelet coefficients in the vertical direction (a-k indicates the wavelet decomposition scales ranging from 1 to 11).

Figure 14 .
Figure 14.(a) HH and (b) HV D-Infs corrected by the MRWCA method.(c) The difference map between the corrected HV and HH D-Infs.

Figure 14 .
Figure 14.(a) HH and (b) HV D-Infs corrected by the MRWCA method.(c) The difference map between the corrected HV and HH D-Infs.

Figure 16 .
Figure 16.The elevation difference distribution (a) without atmospheric correction; (b) with Shirzaei's atmospheric correction method; and (c) with MRWCA atmospheric correction.

Figure 17 .
Figure17.Histograms of the elevation differences between the InSAR-derived DEMs and the ICESat data, generated without ATP correction, with Shirzaei's atmospheric correction method, and with MRWCA atmospheric correction.

Figure 16 .
Figure 16.The elevation difference distribution (a) without atmospheric correction; (b) with Shirzaei's atmospheric correction method; and (c) with MRWCA atmospheric correction.

Figure 17 .
Figure 17.Histograms of the elevation differences between the InSAR-derived DEMs and the ICESat data, generated without ATP correction, with Shirzaei's atmospheric correction method, and with MRWCA atmospheric correction.

Figure 16 .
Figure 16.The elevation difference distribution (a) without atmospheric correction; (b) with Shirzaei's atmospheric correction method; and (c) with MRWCA atmospheric correction.

Figure 16 .
Figure 16.The elevation difference distribution (a) without atmospheric correction; (b) with Shirzaei's atmospheric correction method; and (c) with MRWCA atmospheric correction.

Figure 17 .
Figure17.Histograms of the elevation differences between the InSAR-derived DEMs and the ICESat data, generated without ATP correction, with Shirzaei's atmospheric correction method, and with MRWCA atmospheric correction.

Figure 17 .
Figure17.Histograms of the elevation differences between the InSAR-derived DEMs and the ICESat data, generated without ATP correction, with Shirzaei's atmospheric correction method, and with MRWCA atmospheric correction.

Figure 18 .
Figure 18.Joint distribution of the differences of the HV and HH D-Infs before and after atmospheric effects correction (a) San Francisco test site; (b) Moron test site.

Figure 19 .
Figure 19.The elevation difference maps between InSAR DEMs and SRTM DEM.(a-d) indicates the wavelet decomposition scale with 3, 6, 9 and 11 (the selection of this paper), respectively.

Figure 18 .
Figure 18.Joint distribution of the differences of the HV and HH D-Infs before and after atmospheric effects correction (a) San Francisco test site; (b) Moron test site.

Figure 18 .
Figure 18.Joint distribution of the differences of the HV and HH D-Infs before and after atmospheric effects correction (a) San Francisco test site; (b) Moron test site.

Figure 19 .
Figure 19.The elevation difference maps between InSAR DEMs and SRTM DEM.(a-d) indicates the wavelet decomposition scale with 3, 6, 9 and 11 (the selection of this paper), respectively.
for the two experiment regions.For the two HH D-Infs, the ATP signals make up 63.6% and 44.8% of the total D-Infs signals.For the San Francisco test site, the corrected interferogram could provide a DEM with a RMSE of 7.79 m (Figure

Figure 19 .
Figure 19.The elevation difference maps between InSAR DEMs and SRTM DEM.(a-d) indicates the wavelet decomposition scale with 3, 6, 9 and 11 (the selection of this paper), respectively.
. The residual ATP signals can thus be regarded as the T-ATP signals, as shown in Figure 20c,f.Quantitatively, only 47.2% and 55.7% of the total ATP signals could be detected by Shirzaei's method in these two experiments.However, the MRWCA method can remove both of the TD-ATP and T-ATP signals.Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 18 9), which is an improvement of 40.5% with respect to the DEM before atmospheric effects correction; for the Moron test site, the corrected interferogram could provide a DEM with an RMSE of 10.74 m, which is an improvement of 30.2% with respect to the DEM without atmospheric effects correction (Figure 17), which demonstrates the MRWCA method is effective in detecting the ATP signals with different ATP content conditions.
. The residual ATP signals can thus be regarded as the T-ATP signals, as shown in Figure 20c,f.Quantitatively, only 47.2% and 55.7% of the total ATP signals could be detected by Shirzaei's method in these two experiments.However, the MRWCA method can remove both of the TD-ATP and T-ATP signals.

Figure 20 .
Figure 20.(a) ATP map, (b) TD-ATP map, and (c) T-ATP map for San Francisco test site.(d) ATP map, (e) TD-ATP map, and (f) T-ATP map for Moron Test site.

Figure 20 .
Figure 20.(a) ATP map, (b) TD-ATP map, and (c) T-ATP map for San Francisco test site.(d) ATP map, (e) TD-ATP map, and (f) T-ATP map for Moron Test site.