ATMOSPHERIC PHASE DELAY CORRECTION OF D-INSAR BASED ON SENTINEL-1 A

In this paper, we used the Generic Atmospheric Correction Online Service for InSAR (GACOS) tropospheric delay maps to correct the atmospheric phase delay of the differential interferometric synthetic aperture radar (D-InSAR) monitoring, and we improved the accuracy of subsidence monitoring using D-InSAR technology. Atmospheric phase delay, as one of the most important errors that limit the monitoring accuracy of InSAR, would lead to the masking of true phase in subsidence monitoring. For the problem, this paper used the Sentinel-1A images and the tropospheric delay maps got from GACOS to monitor the subsidence of the Yellow River Delta in Shandong Province. The conventional D-InSAR processing was performed using the GAMMA software. The MATLAB codes were used to correct the atmospheric delay of the D-InSAR results. The results before and after the atmospheric phase delay correction were verified and analyzed in the main subsidence area. The experimental results show that atmospheric phase influences the deformation results to a certain extent. After the correction, the measurement error of vertical deformation is reduced by about 18 mm, which proves that the removal of atmospheric effects can improve the accuracy of the D-InSAR monitoring.


INTRODUCTION
Atmospheric effects are one of the main errors in InSAR applications, and they are important factors that limit the acquisition of high-precision terrain and deformation information.The difference in atmospheric conditions during two SAR imaging moments will bring about atmospheric delay in the propagation process of the electromagnetic wave signal, and then the atmospheric delay phase is produced in the results.These effects are called the atmospheric effects of the InSAR process.In order to improve the accuracy of the results, atmospheric effects need to be reduced as much as possible.Therefore, atmospheric correction becomes an indispensable part of InSAR processing. 1 Current tropospheric delay correction methods are divided into two categories: atmospheric correction based on SAR data itself; atmospheric correction based on external data.Chengsheng Yang analyzed the relationship between the stratified signals and topography of atmospheric wet delays, and he studied the atmospheric delay estimation of SAR interferograms based on terrain and GPS observations [Chengsheng Yang et al., 2011].Wenjun Zhan proposed a strategy for modeling and estimating atomospheric phase of SAR interferogram, which estimated the stratified and turbulent signals from atmosphere.The proposed method was validated with ASAR pair over the Yima area in Henan province [Wenjun Zhan et al., 2015].Liang Chang proposed the use of a differential linear calibration model (DLCM) to calibrate the MODIS infrared water vapor product.Experiments have shown that the DLCM model can effectively 1,2 Corresponding author Xi Li, E-mail address: lix_2018@163.comimprove the estimation accuracy of atmospheric water vapor in eastern China [Liang Chang et al., 2016].
In this paper, GACOS tropospheric delay maps [Yu et al., 2017] are used to correct the atmospheric phase.The key is to propose an iterative tropospheric decomposition interpolation model that decouples the elevation and turbulent tropospheric delay components [Yu et al., 2017].GACOS tropospheric delay maps has the following key features：globally available; operational in a near real time mode.By using the tropospheric delay maps in the Yellow River Delta of Shandong Province, some tentative studies have been made on atmospheric corrections in Sentinel-1A interferograms separated by 84 and 72 days.
The chapters are organized as follows: In section 2, we explain the data used for the experiment and the necessity of atmospheric phase correction.In section 3, we introduce the effect of atmospheric delay on the phase and deformation accuracy of interferograms.In section 4, we use GAMMA software to deal with D-InSAR and correct the atmospheric phase.In section 5, we compare and analyze the deformation results before and after atmospheric correction in Shandong Yellow Triangle.In section 6, we draw some conclusions.

EXPERIMENTAL DATA
We used the Sentinel-1A image of the Yellow River Delta in Shandong as the interference data.The data details are shown in the table 1. Coverage is shown in the red box in Figure 1 (a).The terrain is provided by the SRTM DEM with 30m resolution.The vertical baseline of the interference pair 1 is 26m, and the interference pair 2 is 16m.The elevation difference is greater than 1km at most in the Yellow River Delta.The topography is undulating, and the water vapor distribution is likely to be uneven.So the atmospheric delay seriously influences the interferometric measurement in this area.Therefore, it is necessary to correct the atmosphere of the D-InSAR results.In atmospheric phase correction, high-resolution tropospheric delay maps were used.

Tropospheric Delay Maps
GACOS utilises the Iterative Tropospheric Decomposition (ITD) model (Yu et al., 2017) to separate stratified and turbulent signals from tropospheric total delays, and generate high spatial resolution zenith total delay maps to be used for correcting InSAR measurements and other applications [Yu et al., 2017].Figures 2(a

INFLUENCE OF ATMOPHERIC DELAY ON THE ACCURACY OF INTERFEROGRAM AND DEFORMATION
The atmospheric delay can be divided into the zenith static delay (dry delay) and wet delay, and the dry delay mid-latitude can reach 2.3m, but the surface temperature, pressure, and moisture partial pressure observation values can be used to calculate the model, and the accuracy is up to 1mm [Sasstamoinen et all., 1972].While the wet delay is about 0.3-0.6mand drastic changes, the model estimation accuracy is low.
For side-looking imaging radars, the double-pass phase delay caused by atmospheric changes in single-view complex (SLC) images can be simply expressed as: The influence on the differential interference phase is: Studies have shown that there is irrelevant between atmospheric vapors with time greater than one day [Hanssen et all., 1998].
Therefore, the effect of ZTD error on the interferograms accuracy of can be expressed as: For Sentinel-1A satellites, λ=5.6 cm, inc  =33°.From equation (4) we can see that when the deformation accuracy is better than 1cm, the ZTD accuracy should be better than 6mm.ZTD  including interference delay error and wet delay error, the dry delay accuracy can be estimated to 1mm, so ZTD  can be considered as the wet delay error ZWD  .

Concatenate Sentinel1 SLCs
Due to the large monitoring area in the study area, it is covered by the Sentinel1 image of the upper and lower scenes of the same track.In order to obtain the overall monitoring effect in the monitoring area and save the data processing time.Two scenes of the same track were concatenated.Mosaic image is shown in 4 (a).

Radiometric calibration, Co-registration and Multi-look Processing
SAR image radiometric calibration and co-registration are the basic and necessary steps for interferometric processing.(The geometric position and radiation intensity of the data is unified, and image co-registration and radiation correction are performed before the target point is extracted).
Due to the special imaging mode of the Sentinel1 TOPS mode, the registration accuracy of all images with respect to the master image reaches within 0.009 pixels in the distance.The distance direction: azimuth direction is adopted with 10:2 in multi-look processing.

Simulate Terrain Phase Generation
In the D-InSAR process, DEM needs to be acquired to generate an elevation map in the radar coordinate system, to generate simulated terrain phase and geocoding.In this paper, we use the SRTM DEM data of 30m resolution: after registration with the radar image, The DEM of the study area are obtained in the radar coordinate system; and the simulated topographic phase is generated with the master image for D-InSAR.

Differential Processing and Phase Filtering
Differential processing is performed by using registered master and slave images and simulate topographic phase generated by DEM.We get a differential interferogram, e.g. Figure 4(b).In order to remove the phase noise of differential interferograms and reduce the phase unwrapping error image, it is necessary to filter the differential interferograms, e.g. Figure 4(c).

Differential Interferogram Unwrapping and Geocoding
To mask out the low coherence point and reduce the unwrapping error, phase unwrapping is performed by using a mask file with a coherence threshold of 0.4.We geocode the differential interferometric unwrapping phase image, e.g. Figure 5.

Atmospheric Phase Correction
This article uses the research of real-time mode high-resolution water vapor fields from GPS observations [Yu et al., 2017].The atmospheric phase correction was performed in the MATLAB codes using this tropospheric delay maps.The theory [Yu et al., 2017] is as follows: Tropospheric delays, especially the part due to atmospheric water vapor, vary both vertically and laterally over short distances and are often considered as the sum of (i) a stratified component highly correlated with topography which therefore delineates the vertical tropospheric profile and (ii) a turbulent component resulting from disturbance processes (e.g., severe weather) in the troposphere which trigger uncertain patterns in space and time.We present an iterative tropospheric decomposition (ITD) model to effectively separate the turbulent and elevation-dependent ZTD components.The ITD model decouples the ZTD into stratified and turbulent delays, which enable the more accurate interpolation of dense ZTD fields from pointwise values from a set of GPS reference stations across a region.It is defined as where, for the ZTD at location k, T represents the turbulent component and k x is the station coordinate vector in the local topocentric coordinate system; the stratified component is represented with an exponential function with coefficient β in which 0 L is, for the region considered, the stratified component delay at sea level and  , 1998], so they are not used.
In order to decompose the ZTD into stratified and turbulent components, which can account for substantial amounts of the ZTD but behave very differently, an iterative separation procedure was used: 1).The ZTDs from all GPS stations within the user to reference station decorrelated range limit are used to estimate initial values for the exponential coefficients β and 0 L , assuming that the turbulent component values in equation ( 5) are zero.
2).The residuals k  , which are the summation of the unmodeled errors and the turbulent component, are computed by subtracting per station the stratified delay (as modeled by the estimated exponential coefficients) from the ZTD.
3).The turbulent component, T in equation ( 5), is computed per reference station from the residuals k  by using the IDW function w ui given in equation ( 6): 4).The updated values for the turbulent component per reference station are subtracted from the ZTD per reference station in equation ( 5), and a new set of exponential coefficients are obtained.

5).
Steps 2-4 are repeated until the exponential coefficients β and 0 L converge.The final outputs are the exponential coefficients for the decorrelated range limit considered, plus the turbulent delay component and residuals per reference station.
6).The ZTD at the location of interest is obtained by interpolation of the reference station turbulent component and residuals and added to the stratified delay computed by using the final values of the exponential coefficients.

COMPARISON AND ANALYSIS OF RESULTS
First, the whole deformation maps before and after atmospheric correction are analyzed.Figure 6 is a comparison of the vertical deformation maps before and after atmospheric correction.First, we compare the main deformation areas shown in Figures 6 (a It can be seen that the deformation area is larger before the atmospheric phase correction.In order to better characterize the deformation differences of the settlement areas before and after the atmospheric correction, four north-south section lines were selected in the subsidence area [Cui Xie et all., 2013 ]. Figure 8 D1-D4 shows the north-south profile of the vertical deformation results before and after atmospheric correction.It can be seen that the spatial distribution of the surface deformation of nonatmospheric correction is large, and the value of the deformation variable is generally overestimated by about 11mm.

CONCLUSIONS
In this paper, the process of D-InSAR atmospheric phase correction is realized based on GACOS tropospheric delay maps.We qualitatively and quantitatively compare the results before and after atmospheric correction.The results of atmospheric phase correction experiments in this study indicate: 1) Compared with the results after atmospheric correction, the spatial distribution of deformation monitored by the results without considering atmospheric effects is large, and the deformation value is high (reachable 11mm).There are even false deformation information; 2) Through the atmospheric phase correction in this study, the error of the monitored vertical deformation is reduced by about 18mm.The above conclusions indicate that when using Sentinel-1A data to monitor surface deformation, atmospheric effects will have a greater impact on the final measurement results.In order to obtain more accurate measurement results, InSAR atmospheric correction processing is very necessary.
Figure 1 (b).The monitoring extent is Figure 1 (a).(a) Study area (b) SRTM DEM(unit: m) )(b) and (c)(d) are the tropospheric delay maps used for two interferograms.It can be seen from the figure that there are significant atmospheric effects.(a) 20161011.ztd(b) 20170103.ztd(c) 20161116.ztd(d) 20170127.ztdFigure 2, Tropospheric delay maps(unit: m) change of the radar echo signal λ= radar wavelength ZTD = total zenith delay ZTD  = zenith total delay difference inc  = radar wave incident angle accuracy.Therefore, the relationship between r

Figure 3 ,
Figure 3, Technical process of this study (a) Intensity image after Concatenating (b) Differential interferogram (c) Filtered differential interferogram


represents the remaining unmodeled residual errors, including unmodeled stratified and turbulent signals.The turbulent component usually consists of medium-to-long wavelength signals that can be interpolated by an IDW method.If n GPS stations are used in the region considered, then the IDW model reads as interpolation coefficient; u and i are indices for the user and reference stations, respectively; and ui d represents the horizontal distance from the user to reference station.Reference stations at distances more than 100km from the user station show limited correlation [Emardson and Johansson ) and (b).From experience, the interference fringe in the figure should be a reflection of the deformation of the surface during master-slave image acquisition.In addition, the spatial distribution of surface deformation reflected in the uncorrected results is larger and the deformation value is higher.Comparing Figure6(c) with (d), we can also draw the same conclusion.Moreover, the large red deformed areas in the southeast of 6(b) and (d) disappeared after atmospheric correction, indicating that false deformation information was present in the results without atmospheric correction.(Positive value indicates subsidence): (a) Vertical deformation map before atmospheric correction; (b) Vertical deformation map after atmospheric correction (c) (d) 20161116-20170127(Positive value indicates settlement): (c) Vertical deformation map before atmospheric correction; (d) Vertical deformation map after atmospheric correction Figure 6, Influence Analysis of Atmospheric Effect: Comparison of Vertical Deformation Diagrams Before and After Atmospheric Correction (unit: m) Then select a subsidence area (20161011-20170103) to analyze the deformation results before and after atmospheric correction.Figure 7(b) is the difference map of vertical deformation before and after atmospheric correction.Most of the absolute values is concentrated in the range of 11-25mm.In areas other than the black box in Figure 7(c), except for the two small deformation areas (red color), the value of the vertical deformation after atmospheric correction basically fluctuates 0m, indicating that there is basically no deformation.Only considering the regions with no deformations nearly, the measurement error of vertical deformation after atmospheric correction was reduced by about 18 mm In this experiment.Figure 7 (a) and (c) are the subsidence areas before and after the atmospheric correction.
Figure 8, North-south profile of D1-D4 before and after atmospheric correction

Table 1 ,
Basic parameters of Sentinel-