# Applied Research on InSAR and GPS Data Fusion in Deformation Monitoring

Authors: Ziwen Zhang; Xuelian Wang; Yongdong Wu; Zengpeng Zhao; Yang E
Journal: Scientific Programming
 (2021)

---

## Abstract

With the enrichment of land subsidence monitoring means, data fusion of multisource land subsidence data has gradually become a research hotspot. The Interferometry Synthetic Aperture Radar (InSAR) is a potential Earth observation approach, and it has been verified to have a variety of applications in measuring ground movement, urban subsidence, and landslides but similar to Global Positioning System (GPS). The InSAR observation accuracy and measurements are affected by the tropospheric delay error as well as by the Earth’s ionospheric and tropospheric layers. In order to rectify the InSAR result, there is a need to interpolate the GPS-derived tropospheric delay. Keeping in view of the above, this research study has presented an improved Inverse Distance Weighting (IIDW) interpolation method based on Inverse Distance Weighting (IDW) interpolation by using Sentinel-1 radar satellite image provided by European Space Agency (ESA) and the measured data from the Continuously Operating Reference Stations (CORS) provided by the Survey and Mapping Office of the Lands Department of Hong Kong. Furthermore, the corrected differential tropospheric delay correction is used to correct the InSAR image. The experimental results show that the correction of tropospheric delay by IIDW interpolation not only improves the accuracy of Differential Interferometry Synthetic Aperture Radar (D-InSAR) but also provides a new idea for the solution of InSAR and GPS data fusion.


---

## Body

## 1. Introduction

Synthetic Aperture Radar (SAR) is coherent active microwave remotely sensing equipment [1–3] that has been well recognized for its capacity to effectively record the scattering properties of the ground atmosphere [4–6]. A SAR sensor, which has a side-looking light orientation and can precisely extract the position of an item related to the location of the platform in which the sensor is installed, can be deployed on board an aircraft or satellite. The acquired geometry and the physical features of the imaged section lead to the creation of the obtained scattering radar signal, which, when correctly processed, results in the rebuilding of a complicated-valued high-resolution microwave picture of the surrounding area [2]. Interferometry Synthetic Aperture Radar (InSAR) is a wonderful platform for topographic and ground surface deformation imaging because of its entire-weather and day-to-night imaging ability, wide-ranging geographical cover, high resolution, and measurement accuracy [7–11]. The phase latency that occurs when radio signals travel through the atmosphere is a key source of error in repeat-permit InSAR [12–14], and it can be decreased using GPS techniques. But the three-dimensional resolution of GPS, on the other hand, is poor.In order to fix the InSAR tropospheric delay inaccuracy on a pixel-wise basis, the GPS-derived tropospheric delay should be adjusted. In most cases, tropospheric delay correction is affected not just by alterations in horizontal distance but also by variations in height. The influence of height on tropospheric delay has not been completely studied in standard Inverse Distance Weighting (IDW) interpolation. In this study, an effective aspect for the IDW is created, which is referred to as improved IDW (IIDW). For this experiment, we have used Sentinel-1 satellite images from August 20, 2016, to September 25, 2016, given by ESA, as well as measured data from the CORS provided by the Survey and Mapping Office of the Hong Kong Lands Department. We have also investigated D-InSAR results qualitatively as well as quantitatively in our experiments. After correcting for atmospheric delay, the deformation values in most areas of Hong Kong during the studied time period were between −1.75 mm and 1.71 mm, which were more consistent with actual subsidence in the study area and provided accurate data support for ground subsidence monitoring and forecasting.The remainder of this paper is organized as follows: In Section2, we discuss the related works of researchers; in Section 3, we explain our proposed work; in Section 4, the experimental work of our proposed system is discussed; and, finally, conclusion is presented in Section 5.


## 2. Related Work

The Global Positioning System (GPS) and Interferometry Synthetic Aperture Radar (InSAR) are effective methods for obtaining topography and Earth’s surface movement for Earth’s crust deformation research [12, 13]. According to [14], each of these methods is based on the principle of properly calculating the transit time of an electromagnetic signal generated by satellites orbiting high above the Earth’s surface. The authors in [15] investigated that electromagnetic signals are slowed and propagation speed is slowed as they pass through the troposphere. The degree of latency is principally determined by the geographically and temporally varying temperature and pressure of the atmosphere’s constituent gases, particularly water vapour, and causes major mistakes in repeat-pass InSAR observations. According to [16], significant research has been conducted in the last two centuries in order to better understand and reduce tropospheric influences incorporated into InSAR observations.According to the researchers of [12], a 20% variation in relative humidity in space or time could cause 10 cm of inaccuracies in deformation results. Furthermore, adjusting for low stratification tropospheric latencies can improve unwrapping dramatically over rugged terrain where the interferometric phase may be aliased [17]. Tropospheric impacts in InSAR images can be reduced using two different approaches, according to the authors in [18]: empirical, which is focused on direct projections of the raw InSAR phase delay information, and predictive, which is premised on exterior sets of data such as GPS, climate models, or satellite photos. Tropospheric impacts in InSAR images are decreased using the first method, which involves averaging many separate interferograms. The tropospheric delays are handled as bright lights in [18], which is acceptable if there is no association between topography and tropospheric delay or if the InSAR collection equally samples the time dynamics of the deformation. Using the reverse wavelet transform, the authors in [19] rectified the topographical linked parts by downweighting detected related coefficients, which reflected the impact of the atmospheric delay. The ionospheric impact on InSAR signals impacts both azimuth subpixel offsets and differential interferograms, according to the authors in [20]. The first bias is created by directional changes in azimuth pixels’ relative subpixel location, while the second bias is caused by comparative lengthening of wave pathways across two InSAR collections, which affects the interferogram image. Inspired from the above works, this research study aims to develop an improved IDW interpolation method based on IDW interpolation by using Sentinel-1 radar satellite image provided by ESA and the measured data from the CORS provided by the Survey and Mapping Office of the Lands Department of Hong Kong.


## 3. Proposed System

In this section, we will explain our proposed research work on InSAR and GPS data fusion in deformation monitoring; at the beginning, we will give an overview of Global Positioning System (GPS) and Interferometry Synthetic Aperture Radar (InSAR) along with their connection to each other. After that, we will explain the InSAR and GPS data fusion process using a diagram (Figure1). Finally, we will compare two-difference algorithms, Inverse Distance Weighting (IDW) and our proposed improved Inverse Distance Weighting (IIDW), for tropospheric delay correction.Figure 1
Flow chart of InSAR and GPS data fusion.Level measurement is one of the most important means to monitor ground subsidence and it has attracted great attention these days. With the continuous progress of science and technology, in recent years, the means of ground settlement monitoring has been enriched, with the emergence of Global Positioning System (GPS), Interferometry Synthetic Aperture Radar (InSAR), etc. Due to the inherent characteristics of radar measurement, the InSAR technology is incomparable to other methods by virtue of its advantages of all-day, all-weather ability, spatial continuity, high accuracy, and contactlessness [21]. On other hand, the GPS technology can accurately determine tropospheric and ionospheric delays and achieve high accuracy positioning at the same time. Therefore, data fusion of InSAR and GPS data can not only correct the satellite orbit errors and atmospheric delay errors that are difficult to be eliminated by InSAR data itself but also realize the effective unification of high spatial resolution and high-range deformation accuracy of InSAR technology and high temporal resolution and high plane position accuracy of GPS technology [22–26].GPS can obtain single-point continuous high-precision plane position data, while InSAR can provide higher-precision wide-range monitoring results for vertical ground deformation. Therefore, the fusion of GPS data and InSAR data can realize the effective combination of high spatial resolution and high elevation deformation accuracy of InSAR technology and high temporal resolution and high planar position accuracy of GPS and correct the errors that are difficult to be eliminated by InSAR data itself. With the purpose of studying the fusion of InSAR and multisource data, GPS is used to correct the atmospheric delay error of InSAR to better reflect the subsidence pattern of the ground subsidence region [27].


### 3.1. InSAR and GPS Data Fusion

InSAR (Interferometry with Synthetic Aperture Radar) is a technology that derives ground height and distortion information using phase’s information extracted from Synthetic Aperture Radar multiple data. Differential Interferometry Synthetic Aperture Radar (D-InSAR) is based on InSAR, and its theoretical precision when used in subsurface deformation monitoring can reach subcentimeter levels. The rest of this section will debate on the process of InSAR along with GPS data fusion and tropospheric delay correction techniques.


### 3.2. InSAR and GPS Data Fusion Process

The study of InSAR and GPS data fusion for ground settlement monitoring is proposed based on a thorough collection, review, and analysis of related literature. During the process of InSAR and GPS data fusion, in our proposed system, we have used GPS data in order to correct the atmospheric delay error of D-InSAR results; after that, this GPS data is used to accurately determine the satellite orbit parameters, and, finally, the fusion effect is evaluated, and further ground subsidence is analyzed and predicted. The main technical process is shown in Figure1.


### 3.3. Two-Difference Algorithm for Tropospheric Delay Correction

For tropospheric delay correction, assume that point A on the SAR image is fixed and B is another point on top of the SAR image. Assuming that the atmospheric delay of point A at SAR imagej estimated from GPS is DAj and the atmospheric delay of point B at SAR image j is DBj, the delay difference between the two stations can be obtained as(1)DABj=DBj−DAj.With point A as the fixed reference point, the single differential delay between other GPS points and point A can be solved according to the above equation, and then the obtained single differential delay is interpolated to obtain the atmospheric delay correction image map.Consider two sites A and B while assuming two timesj (primary image) and k (secondary image); two single-difference scores can be calculated from equation (1):(2)DABj=DBj−DAj,DABk=DBk−DAk.The two single-difference results are further differenced once again, that is, double-differenced.(3)DABjk=DABk−DABj=DBk−DAk−DBj−DAj=DBk−DBj−DAj−DAj.


### 3.4. Interpolation Model for Tropospheric Delay Correction

The station spacing of GPS networks is generally tens or even hundreds of kilometers, which causes the spatial resolution of the GPS-acquired atmospheric delay correction to be unable to meet the interpolation requirements of InSAR imagery. The atmospheric delay correction obtained by processing GPS needs to be encrypted and interpolated to further provide effective atmospheric delay correction for InSAR images [28–30]. There are two main types of interpolation methods: random and deterministic.


#### 3.4.1. Kriging Interpolation

Kriging, commonly known as Gaussian process regression, is a geostatistical interpolation technique. Interpolated values are represented using a Gaussian process guided by previous covariance matrix in Kriging. Kriging is a technique for predicting values in a particular region.Assuming that there aren known points in the desired interpolation region, a linear combination can be used to obtain an estimate of any point in the desired interpolation region:(4)Zv∗=∑i=1nλiZxi.In equation (5), Zxi is the value of the known points, and the weighting factor is denoted by the weighting factor, which shows the influence of the known points on the points to be estimated.The Kriging interpolation method stipulates that the estimated value of any point in the interpolation area and its true value have the following relationship under unbiased conditions:(5)EZv∗−Zv=0,EZv∗=E∑i=1nλiZxi=∑i=1nλiEZxi=EZv=m,where m is the mathematical expectation.Because(6)EZxi=m,it results in(7)∑i=1nλi=1.Kriging’s interpolation method specifies that the estimation variance minimization condition should be satisfied in the solution process as follows:(8)EZv∗−Zv2=min.The equation for the variance in the solution is(9)σE2=EZv∗−Zv2=∑i=0n∑j=0nλiλjCxi,xj=Cx0,x0−2∑j=1nλiCx0,xj+∑i=0n∑j=0nλiλjCxi,xj.In equation (9), Cxi,xj is the covariance function of xi and xj.The computational process of Kriging interpolation is very simple and the implementation process mainly includes statistical analysis of data, simulation of variorums, creation of surfaces, and detection of surface changes [31].


#### 3.4.2. Inverse Distance Weighting (IDW) Interpolation

Inverse Distance Weighting (IDW) is based on the principle that each known point in the region to be interpolated has a local influence, which becomes smaller as the distance increases. Observations closer to the unmeasured point have more influence on the unmeasured point than observations farther away from the unmeasured point, and the magnitude of the influence is expressed in terms of weights, with the weights of points closer to the expected point being larger than the weights of points farther away from the expected point.IDW equation for the atmospheric delay obtained from GPS observations is(10)D^λ0,ϕ0=∑i=1NwiDλi,ϕi.In equation (10), the interpolated atmospheric delay at the point with coordinates East λ0 and North ϕ0 is D^λ0,ϕ0.The GPS delay correction at the point with coordinates Eastλ0 and North ϕ0 is Dλi,ϕi.The weights of IDW are(11)wi=di0−P∑i=1Ndi0−P,where wi denotes the weight related to the delay of GPS acquisition. The number of GPS points near the expected point used for interpolation is denoted by N. The interpolated weights become smaller as the distance increases.(12)∑i=1Nwi=1,where di0=x0−x22+y0−y22+z0−z22.In equation (12), the distance from the sample point i to the point to be interpolated is di0=i=1,2,…,N and P is the power parameter.


#### 3.4.3. Improved Inverse Distance Weighting (IIDW) Interpolation

From the principle of IDW, it is clear that the IDW interpolation method treats the horizontal distance and elevation difference between the data points and the points to be interpolated to the same degree of influence on the points to be interpolated. This is an ideal state, which does not correspond to the actual situation. In order to better simulate the most realistic situation for interpolation, it is necessary to differentiate the influence of the horizontal distance and elevation difference on the interpolation points. Based on the above IDW principle, an influence factorα is introduced; that is, equation (12) is rewritten as(13)di0=αdi0s2+1−αdi0h2α,where di0s=x0−x22+y0−y22 is the horizontal distance from the data point to the point to be interpolated; the difference in elevation between data point i and the interpolated point is denoted by di0h=hi−h0i=1,2,…,N.Therefore,(14)D^λ0,ϕ0=∑i=1NwiDλi,ϕi,wi=di0−P∑i=1Ndi0−P,∑i=1Nwi=1,where di0=αdi0s2+1−αdi0h2/α.These three equations constitute the improved Inverse Distance Weighting (IIDW) method. The influence factor can be determined by achieving the optimal interpolation accuracy through multiple trials.


## 3.1. InSAR and GPS Data Fusion

InSAR (Interferometry with Synthetic Aperture Radar) is a technology that derives ground height and distortion information using phase’s information extracted from Synthetic Aperture Radar multiple data. Differential Interferometry Synthetic Aperture Radar (D-InSAR) is based on InSAR, and its theoretical precision when used in subsurface deformation monitoring can reach subcentimeter levels. The rest of this section will debate on the process of InSAR along with GPS data fusion and tropospheric delay correction techniques.


## 3.2. InSAR and GPS Data Fusion Process

The study of InSAR and GPS data fusion for ground settlement monitoring is proposed based on a thorough collection, review, and analysis of related literature. During the process of InSAR and GPS data fusion, in our proposed system, we have used GPS data in order to correct the atmospheric delay error of D-InSAR results; after that, this GPS data is used to accurately determine the satellite orbit parameters, and, finally, the fusion effect is evaluated, and further ground subsidence is analyzed and predicted. The main technical process is shown in Figure1.


## 3.3. Two-Difference Algorithm for Tropospheric Delay Correction

For tropospheric delay correction, assume that point A on the SAR image is fixed and B is another point on top of the SAR image. Assuming that the atmospheric delay of point A at SAR imagej estimated from GPS is DAj and the atmospheric delay of point B at SAR image j is DBj, the delay difference between the two stations can be obtained as(1)DABj=DBj−DAj.With point A as the fixed reference point, the single differential delay between other GPS points and point A can be solved according to the above equation, and then the obtained single differential delay is interpolated to obtain the atmospheric delay correction image map.Consider two sites A and B while assuming two timesj (primary image) and k (secondary image); two single-difference scores can be calculated from equation (1):(2)DABj=DBj−DAj,DABk=DBk−DAk.The two single-difference results are further differenced once again, that is, double-differenced.(3)DABjk=DABk−DABj=DBk−DAk−DBj−DAj=DBk−DBj−DAj−DAj.


## 3.4. Interpolation Model for Tropospheric Delay Correction

The station spacing of GPS networks is generally tens or even hundreds of kilometers, which causes the spatial resolution of the GPS-acquired atmospheric delay correction to be unable to meet the interpolation requirements of InSAR imagery. The atmospheric delay correction obtained by processing GPS needs to be encrypted and interpolated to further provide effective atmospheric delay correction for InSAR images [28–30]. There are two main types of interpolation methods: random and deterministic.


### 3.4.1. Kriging Interpolation

Kriging, commonly known as Gaussian process regression, is a geostatistical interpolation technique. Interpolated values are represented using a Gaussian process guided by previous covariance matrix in Kriging. Kriging is a technique for predicting values in a particular region.Assuming that there aren known points in the desired interpolation region, a linear combination can be used to obtain an estimate of any point in the desired interpolation region:(4)Zv∗=∑i=1nλiZxi.In equation (5), Zxi is the value of the known points, and the weighting factor is denoted by the weighting factor, which shows the influence of the known points on the points to be estimated.The Kriging interpolation method stipulates that the estimated value of any point in the interpolation area and its true value have the following relationship under unbiased conditions:(5)EZv∗−Zv=0,EZv∗=E∑i=1nλiZxi=∑i=1nλiEZxi=EZv=m,where m is the mathematical expectation.Because(6)EZxi=m,it results in(7)∑i=1nλi=1.Kriging’s interpolation method specifies that the estimation variance minimization condition should be satisfied in the solution process as follows:(8)EZv∗−Zv2=min.The equation for the variance in the solution is(9)σE2=EZv∗−Zv2=∑i=0n∑j=0nλiλjCxi,xj=Cx0,x0−2∑j=1nλiCx0,xj+∑i=0n∑j=0nλiλjCxi,xj.In equation (9), Cxi,xj is the covariance function of xi and xj.The computational process of Kriging interpolation is very simple and the implementation process mainly includes statistical analysis of data, simulation of variorums, creation of surfaces, and detection of surface changes [31].


### 3.4.2. Inverse Distance Weighting (IDW) Interpolation

Inverse Distance Weighting (IDW) is based on the principle that each known point in the region to be interpolated has a local influence, which becomes smaller as the distance increases. Observations closer to the unmeasured point have more influence on the unmeasured point than observations farther away from the unmeasured point, and the magnitude of the influence is expressed in terms of weights, with the weights of points closer to the expected point being larger than the weights of points farther away from the expected point.IDW equation for the atmospheric delay obtained from GPS observations is(10)D^λ0,ϕ0=∑i=1NwiDλi,ϕi.In equation (10), the interpolated atmospheric delay at the point with coordinates East λ0 and North ϕ0 is D^λ0,ϕ0.The GPS delay correction at the point with coordinates Eastλ0 and North ϕ0 is Dλi,ϕi.The weights of IDW are(11)wi=di0−P∑i=1Ndi0−P,where wi denotes the weight related to the delay of GPS acquisition. The number of GPS points near the expected point used for interpolation is denoted by N. The interpolated weights become smaller as the distance increases.(12)∑i=1Nwi=1,where di0=x0−x22+y0−y22+z0−z22.In equation (12), the distance from the sample point i to the point to be interpolated is di0=i=1,2,…,N and P is the power parameter.


### 3.4.3. Improved Inverse Distance Weighting (IIDW) Interpolation

From the principle of IDW, it is clear that the IDW interpolation method treats the horizontal distance and elevation difference between the data points and the points to be interpolated to the same degree of influence on the points to be interpolated. This is an ideal state, which does not correspond to the actual situation. In order to better simulate the most realistic situation for interpolation, it is necessary to differentiate the influence of the horizontal distance and elevation difference on the interpolation points. Based on the above IDW principle, an influence factorα is introduced; that is, equation (12) is rewritten as(13)di0=αdi0s2+1−αdi0h2α,where di0s=x0−x22+y0−y22 is the horizontal distance from the data point to the point to be interpolated; the difference in elevation between data point i and the interpolated point is denoted by di0h=hi−h0i=1,2,…,N.Therefore,(14)D^λ0,ϕ0=∑i=1NwiDλi,ϕi,wi=di0−P∑i=1Ndi0−P,∑i=1Nwi=1,where di0=αdi0s2+1−αdi0h2/α.These three equations constitute the improved Inverse Distance Weighting (IIDW) method. The influence factor can be determined by achieving the optimal interpolation accuracy through multiple trials.


## 3.4.1. Kriging Interpolation

Kriging, commonly known as Gaussian process regression, is a geostatistical interpolation technique. Interpolated values are represented using a Gaussian process guided by previous covariance matrix in Kriging. Kriging is a technique for predicting values in a particular region.Assuming that there aren known points in the desired interpolation region, a linear combination can be used to obtain an estimate of any point in the desired interpolation region:(4)Zv∗=∑i=1nλiZxi.In equation (5), Zxi is the value of the known points, and the weighting factor is denoted by the weighting factor, which shows the influence of the known points on the points to be estimated.The Kriging interpolation method stipulates that the estimated value of any point in the interpolation area and its true value have the following relationship under unbiased conditions:(5)EZv∗−Zv=0,EZv∗=E∑i=1nλiZxi=∑i=1nλiEZxi=EZv=m,where m is the mathematical expectation.Because(6)EZxi=m,it results in(7)∑i=1nλi=1.Kriging’s interpolation method specifies that the estimation variance minimization condition should be satisfied in the solution process as follows:(8)EZv∗−Zv2=min.The equation for the variance in the solution is(9)σE2=EZv∗−Zv2=∑i=0n∑j=0nλiλjCxi,xj=Cx0,x0−2∑j=1nλiCx0,xj+∑i=0n∑j=0nλiλjCxi,xj.In equation (9), Cxi,xj is the covariance function of xi and xj.The computational process of Kriging interpolation is very simple and the implementation process mainly includes statistical analysis of data, simulation of variorums, creation of surfaces, and detection of surface changes [31].


## 3.4.2. Inverse Distance Weighting (IDW) Interpolation

Inverse Distance Weighting (IDW) is based on the principle that each known point in the region to be interpolated has a local influence, which becomes smaller as the distance increases. Observations closer to the unmeasured point have more influence on the unmeasured point than observations farther away from the unmeasured point, and the magnitude of the influence is expressed in terms of weights, with the weights of points closer to the expected point being larger than the weights of points farther away from the expected point.IDW equation for the atmospheric delay obtained from GPS observations is(10)D^λ0,ϕ0=∑i=1NwiDλi,ϕi.In equation (10), the interpolated atmospheric delay at the point with coordinates East λ0 and North ϕ0 is D^λ0,ϕ0.The GPS delay correction at the point with coordinates Eastλ0 and North ϕ0 is Dλi,ϕi.The weights of IDW are(11)wi=di0−P∑i=1Ndi0−P,where wi denotes the weight related to the delay of GPS acquisition. The number of GPS points near the expected point used for interpolation is denoted by N. The interpolated weights become smaller as the distance increases.(12)∑i=1Nwi=1,where di0=x0−x22+y0−y22+z0−z22.In equation (12), the distance from the sample point i to the point to be interpolated is di0=i=1,2,…,N and P is the power parameter.


## 3.4.3. Improved Inverse Distance Weighting (IIDW) Interpolation

From the principle of IDW, it is clear that the IDW interpolation method treats the horizontal distance and elevation difference between the data points and the points to be interpolated to the same degree of influence on the points to be interpolated. This is an ideal state, which does not correspond to the actual situation. In order to better simulate the most realistic situation for interpolation, it is necessary to differentiate the influence of the horizontal distance and elevation difference on the interpolation points. Based on the above IDW principle, an influence factorα is introduced; that is, equation (12) is rewritten as(13)di0=αdi0s2+1−αdi0h2α,where di0s=x0−x22+y0−y22 is the horizontal distance from the data point to the point to be interpolated; the difference in elevation between data point i and the interpolated point is denoted by di0h=hi−h0i=1,2,…,N.Therefore,(14)D^λ0,ϕ0=∑i=1NwiDλi,ϕi,wi=di0−P∑i=1Ndi0−P,∑i=1Nwi=1,where di0=αdi0s2+1−αdi0h2/α.These three equations constitute the improved Inverse Distance Weighting (IIDW) method. The influence factor can be determined by achieving the optimal interpolation accuracy through multiple trials.


## 4. Experimental Work

For experimental work, we have used cross validation technique, wherein 36 days’ session data observed from August 20, 2016 to September 25, 2016 by the Survey and Mapping Office of the Lands Department of Hong Kong were used to investigate the efficiency of the proposed IIDW. Out of 18 stations, 17 stations were considered as measured locations and the remaining one was used as prediction location for determining the tropospheric delay correction and comparison with its GPS-derived delay. This process is repeated until all tropospheric delays of 18 stations were predicted and compared with their GPS-derived delays. Figure2 shows the heights and locations of the GPS sites for this experiment.Figure 2
Heights of the GPS site.


### 4.1. GPS Actual Measurement Data Processing

The “Hong Kong Satellite Positioning Reference Station Network” of the Survey and Mapping Office of the Hong Kong Lands Department, SatRef, consists of 18 continuously operating CORS evenly distributed throughout Hong Kong, as shown in Figure3.Figure 3
Distribution map of the measured stations of the CORS.Among the 18 sites in the SatRef system, the monitoring data of the sites with gaps were removed because there were gaps in the data of one or two sites. Ten of these sites were selected as GPS observation sites (blue rectangular points in Figure3) and five sites were selected as predicted sites for tropospheric delay correction (red circular points in Figure 3). The significance of the predicted sites is that the predicted values can be compared with the measured values obtained by solving these sites.GPS observations were made during one hour before and after the radar satellite pass in the study area with a sampling interval of 30 s. Two days of observations from the SatRef system were used for data processing with the GAMIT software. During the data processing, the tropospheric delay correction of the station was set to a 30 min interval, so that there were a total of four tropospheric correction parameters in the data processing. In this experiment, the Sha Tin (HKST) station is set as the reference point, and it is assumed that it remains unchanged. Based on this station as the reference, the single-difference tropospheric correction of other stations with this point is further calculated. This requires another differential on top of the single differential in the previous step to satisfy the tropospheric delay correction of the InSAR image, that is, double differential. This results in the double-differential tropospheric delay correction of InSAR images. The double-difference results for most stations are smaller than the single-difference results, and the results of double-difference range from −5.02 cm to 5.34 cm, and this correction is of great importance for D-InSAR to obtain subcentimeter-level accuracy.Three different interpolation methods, namely, Kriging, IDW, and IIDW interpolation methods, are selected to interpolate the double-difference results, and the IIDW interpolation method yields a plot of the double-difference results, as shown in Figure4, and the obtained results can be directly corrected for the convective layer delay of the D-InSAR data.Figure 4
Interpolation of IIDW. (a) 30 min. (b) 60 min. (c) 90 min. (d) 120 min.
(a)(b)(c)(d)


### 4.2. Sentinel-1 Satellite Measured Data Processing

Sentinel-1 satellite images from August 20, 2016, to September 25, 2016, were selected for the experimental data. The deformation image pair time interval was 36 days and interferometry was performed on the two images using the D-InSAR method.The baseline estimation is performed first, followed by the interferogram generation, the filtering and correlation calculations are performed after getting the processing results of the previous step, the filtered interferogram and coherence coefficient map are obtained, the phase deconvolution processing is performed in the next step, followed by the phase transformation as well as the geocoding, and the final deformation data can be obtained through the above steps, as shown in Figures5 and 6.Figure 5
D-InSAR results before atmospheric delay corrections of the test area.Figure 6
D-InSAR results after atmospheric delay corrections of the test area.The subsidence map of the test area was obtained by D-InSAR data processing technique and combined with GIS technology. From Figures5 and 6, it can be seen that the maximum subsidence value of the ground surface before the atmospheric delay correction is 3.85 mm and the maximum uplift value is 11.98 mm, and the maximum subsidence value of the ground surface after the atmospheric delay correction is 4.02 mm and the maximum uplift value is 10.07 mm. From Figure 6, it can be seen that the subsidence changes in most areas of the study area are small after the atmospheric delay correction, and it can be combined with the approximate location; and extent of ground subsidence in the studied area can be obtained by combining the geographic location map of the study area. In the study area, some areas show some minor subsidence, some areas form “funnel-shaped subsidence,” and some areas experience uplift. In the study area covered by the radar images, there are mainly two obvious subsidence targets. One of them is located in the northeast of the image in Tai Po, Luohu, and Yantian districts, reaching a maximum subsidence value of 4 mm; the other is located in the southwest of the image in the island, reaching a maximum uplift value of 10 mm, which may be caused by land reclamation. D-InSAR results before atmospheric delay corrections of test area are shown in Figure 5. D-InSAR results after atmospheric delay corrections of test area are shown in Figure 6.


## 4.1. GPS Actual Measurement Data Processing

The “Hong Kong Satellite Positioning Reference Station Network” of the Survey and Mapping Office of the Hong Kong Lands Department, SatRef, consists of 18 continuously operating CORS evenly distributed throughout Hong Kong, as shown in Figure3.Figure 3
Distribution map of the measured stations of the CORS.Among the 18 sites in the SatRef system, the monitoring data of the sites with gaps were removed because there were gaps in the data of one or two sites. Ten of these sites were selected as GPS observation sites (blue rectangular points in Figure3) and five sites were selected as predicted sites for tropospheric delay correction (red circular points in Figure 3). The significance of the predicted sites is that the predicted values can be compared with the measured values obtained by solving these sites.GPS observations were made during one hour before and after the radar satellite pass in the study area with a sampling interval of 30 s. Two days of observations from the SatRef system were used for data processing with the GAMIT software. During the data processing, the tropospheric delay correction of the station was set to a 30 min interval, so that there were a total of four tropospheric correction parameters in the data processing. In this experiment, the Sha Tin (HKST) station is set as the reference point, and it is assumed that it remains unchanged. Based on this station as the reference, the single-difference tropospheric correction of other stations with this point is further calculated. This requires another differential on top of the single differential in the previous step to satisfy the tropospheric delay correction of the InSAR image, that is, double differential. This results in the double-differential tropospheric delay correction of InSAR images. The double-difference results for most stations are smaller than the single-difference results, and the results of double-difference range from −5.02 cm to 5.34 cm, and this correction is of great importance for D-InSAR to obtain subcentimeter-level accuracy.Three different interpolation methods, namely, Kriging, IDW, and IIDW interpolation methods, are selected to interpolate the double-difference results, and the IIDW interpolation method yields a plot of the double-difference results, as shown in Figure4, and the obtained results can be directly corrected for the convective layer delay of the D-InSAR data.Figure 4
Interpolation of IIDW. (a) 30 min. (b) 60 min. (c) 90 min. (d) 120 min.
(a)(b)(c)(d)


## 4.2. Sentinel-1 Satellite Measured Data Processing

Sentinel-1 satellite images from August 20, 2016, to September 25, 2016, were selected for the experimental data. The deformation image pair time interval was 36 days and interferometry was performed on the two images using the D-InSAR method.The baseline estimation is performed first, followed by the interferogram generation, the filtering and correlation calculations are performed after getting the processing results of the previous step, the filtered interferogram and coherence coefficient map are obtained, the phase deconvolution processing is performed in the next step, followed by the phase transformation as well as the geocoding, and the final deformation data can be obtained through the above steps, as shown in Figures5 and 6.Figure 5
D-InSAR results before atmospheric delay corrections of the test area.Figure 6
D-InSAR results after atmospheric delay corrections of the test area.The subsidence map of the test area was obtained by D-InSAR data processing technique and combined with GIS technology. From Figures5 and 6, it can be seen that the maximum subsidence value of the ground surface before the atmospheric delay correction is 3.85 mm and the maximum uplift value is 11.98 mm, and the maximum subsidence value of the ground surface after the atmospheric delay correction is 4.02 mm and the maximum uplift value is 10.07 mm. From Figure 6, it can be seen that the subsidence changes in most areas of the study area are small after the atmospheric delay correction, and it can be combined with the approximate location; and extent of ground subsidence in the studied area can be obtained by combining the geographic location map of the study area. In the study area, some areas show some minor subsidence, some areas form “funnel-shaped subsidence,” and some areas experience uplift. In the study area covered by the radar images, there are mainly two obvious subsidence targets. One of them is located in the northeast of the image in Tai Po, Luohu, and Yantian districts, reaching a maximum subsidence value of 4 mm; the other is located in the southwest of the image in the island, reaching a maximum uplift value of 10 mm, which may be caused by land reclamation. D-InSAR results before atmospheric delay corrections of test area are shown in Figure 5. D-InSAR results after atmospheric delay corrections of test area are shown in Figure 6.


## 5. Conclusion

This research work puts forward an improved Inverse Distance Weighting (IIDW) interpolation method based on traditional Inverse Distance Weighting (IDW) interpolation method for InSAR and GPS data fusion in deformation monitoring. The Inverse Distance Weighting (IDW) interpolation method does not fully consider the effect of height on tropospheric delay, so we have analyzed the influence of atmospheric delay on InSAR, where the differential model of atmospheric delay correction is studied, and the InSAR data processing process based on GPS atmospheric delay correction is designed. Furthermore, GPS observation data and D-InSAR data were fused and studied by using GPS monitoring data from Hong Kong SatRef network and Sentinel-1 radar satellite monitoring data. Finally, the obtained D-InSAR measurements were analyzed qualitatively as well as quantitatively, where total tropospheric delay over the Hong Kong SatRef network stations was solved by using GAMIT software, and the difference maps of atmospheric delay were generated according to Kriging interpolation, IDW interpolation, and IIDW interpolation methods. The experimental results after atmospheric delay correction showed that the deformation values in most areas of Hong Kong were between −1.75 mm and 1.71 mm during the studied time period, which were more consistent with the actual subsidence in the study area and provided accurate data support for ground subsidence monitoring and forecasting.



---

*Source: 3888975-2021-07-23.xml*