# GPS and InSAR Time Series Analysis: Deformation Monitoring Application in a Hydraulic Engineering Resettlement Zone, Southwest China

Authors: Ruya Xiao; Xiufeng He
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Booming development of hydropower in China has resulted in increasing concerns about the related resettlement issues. Both global positioning system (GPS) and persistent scatterer interferometric synthetic aperture radar (InSAR) time series analysis are applied to measuring the magnitude and monitoring the spatial and temporal variations of land surface displacement in Hanyuan, a hydraulic engineering resettlement zone, southwest China. The results from the GPS monitoring system established in Hanyuan match well the digital inclinometer results, suggesting that the GPS monitoring system can be employed as a complement to the traditional ground movement monitoring methods. The InSAR time series witness various patterns and magnitudes of deformation in the resettlement zone. Combining the two complementary techniques will overcome the limitations of the single method.


---

## Body

## 1. Introduction

Pubugou Hydropower Station with a 186 m high core wall rock-fill dam located in the middle reaches of Dadu River, southwest China, is built for power generation, flood control, and sediment retaining retention. Several towns have been inundated by the backwater, and the population in the whole area needs to relocate, involving more than 100,000 residents, which in resettlement scale is next only to the Three Gorges and Xiaolangdi water control project among modern hydraulic engineering projects in China.The resettlement zone, new Hanyuan County, located on Luobogang hillock, consists of unstable broken limestone, weathered shale, and schist shone with laminated loess. In 2009 and 2010, two landslide failures occurred in Hanyuan, resulting in great loss of lives and properties. Thus, to detect early indications of catastrophic movements, monitoring is essential to predict landslides.However, the current conventional deformation monitoring techniques are always time-consuming and costly. In addition, even if taken regularly, the monitoring points are not usually dense enough to assist in understanding the mechanisms of the landslide. Therefore, in this paper, we address deformation monitoring applications through several geodetic techniques, that is, global positioning system (GPS) and interferometric synthetic aperture radar (InSAR) in the reservoir resettlement county.The remainder of the paper is organized as follows. Section2 briefly introduces the principles of GPS and InSAR time series analysis. Site conditions and SAR data used are described in Section 3. In Section 4, we present the GPS monitoring system established in the study area and comparison between GPS and digital inclinometers in a landslide case. Then, time series results obtained by both GPS and InSAR are also displayed. Finally, conclusions of the paper are drawn in Section 5.


## 2.  Method

Two modern geodetic techniques developed since the late 1980s, global positioning system (GPS) and interferometric synthetic aperture radar (InSAR), have revolutionized the way we observe the Earth surface. GPS can provide 3D information at each GPS observation points while InSAR can image the line of sight (LOS) component of deformation over a large area at spatial resolution of tens of meters [1–4]. The technology using GPS to monitor deformation is mature; thus, we just briefly introduce the multiantenna monitoring system, aiming at decreasing the device investment.


### 2.1. GPS Multiantenna Monitoring System

GPS is considered as one of the best techniques for deformation monitoring thanks to its higher automation and lower labor consuming than that of traditional geodetic survey techniques [5–8]. Depending on density of the GPS antenna, the spatial resolution is low. Therefore, GPS does have disadvantages, with the major drawback being the high cost associated with placing a permanent GPS receiver at each monitoring point.The approach of linking a single GPS receiver with multiple antennas mounted at several monitoring points has been implemented. We use a dedicated electronic timed switching device, the GPS multiantenna switch (GMAS), to connect one receiver with individual antennas. This significantly reduces the number of receivers required, thus decreasing hardware investment. GMAS has been patented in China (Patent no: 00219891.6, Owners: Xiufeng He, Xiaoli Ding, Yongrong Sun, Yongqi Chen et al.). The receiver conducts standard pseudo-range and carrier-phase observations for each antenna [9–12]. It takes efficient advantage of the characteristics of deformation survey, such as the position of a surveyed point, which changes relatively little between two consecutive survey epochs. Figure 1 shows the outline of multiantenna GPS system.Figure 1
Outline of multiantenna GPS system for deformation monitoring.


### 2.2. Persistent Scatterers InSAR Time Series

InSAR can provide measurements of surface movements at a considerably improved spatial resolution with quite high accuracy, that is, centimeter or even millimeter level, over the whole area. Thus, repeat-pass satellite InSAR is a rapidly evolving remote sensing technology for observing the Earth surface and now well developed with many applications [13–18].Due to the slow surface displacements rate, phase ambiguity problems related to temporal and geometrical decorrelations, and strong tropospheric effects, it is difficult to obtain reliable measurements of subsidence in densely vegetated areas using conventional InSAR techniques. In the late 1990s it is noticed that some radar targets maintain stable backscattering characteristics for a period of month or years, and the phase information from these stable targets (hereafter called Persistent Scatterers, PS) can be used, even over a long time period. This led to the development of PS InSAR methodologies, adopting both amplitude stability and coherence stability (i.e., correlation) as pixel selection criteria [19–21]. The persistent scatterers technique is a powerful tool that overcomes the above-said limitation of routine InSAR by simultaneously exploiting all the available SAR images gathered during repeated satellite passes.The deformation was modeled in time and a temporal model of evolution was fitted to the double difference phase during the process as the phase is unwrapped. By taking the phase difference between nearby persistent scatterer candidates, atmospheric and orbit errors are mainly reduced. Besides, signals lacking in consistency in both space and time dimension are always considered as the atmospheric error rather than the deformation signal. By high-pass filtering the unwrapped phase in time and low-pass filtering in space we estimate the spatially correlated errors. Deformation and DEM error contributing to the double difference phase are later modeled for the whole time series. The residuals between the model and the double phase signal were considered as the noise level. Then for each persistent scatterer candidates, network adjustment is applied to estimating the noise.Stanford method for persistent scatterers (StaMPS) is used for analyzing crustal deformation in nonurban environments in this paper. Unlike Ferretti method, StaMPS produces a time series of deformation, with no prior assumptions about the temporal nature of deformation. Therefore, it can get better performance and be more reliable under the circumstances of unknown or not readily parameterized deformation. Here only PS identification and selection steps of the innovative PS processing method are discussed; more detailed information about the algorithm, such as the phase unwrapping and spatially correlated terms removal, can be found in [22].


#### 2.2.1. Amplitude Analysis

Ferretti et al. [19] defined the amplitude dispersion index as
(1)DA≡σAmA,
where σA and mA are, respectively, the standard deviation and the mean of the amplitude values. For a constant signal and high signal to noise ratio (SNR), DA≈σΦ, where σΦ is the phase standard deviation. As lower DA indicates higher phase stability, consideration of amplitude is useful both to reduce the number of pixels for phase analysis and to better estimate the probability of a pixel being a PS. The pixels selected as PS candidates have DA≤0.4 [23].


#### 2.2.2. Phase Analysis

The residual phase (with flat earth and DEM reference phase removal),Φint, of the xth pixel in the ith interferogram can be written as
(2)Φint,x,i=Φdef,x,i+ΔΦε,x,i+Φatm,x,i+ΔΦorb,x,i+Φn,x,i,
where Φdef is the phase change due to the deformation of the pixel in the satellite line-of-sight (LOS) direction, ΔΦε is the residual topographic phase due to DEM error, Φatm is the phase equivalent of the difference in atmospheric retardation between passes, ΔΦorb is the residual phase due to orbit inaccuracies, and Φn is the noise term due to variability in scattering from the pixel, thermal noise, and coregistration errors.The spatially correlated phase is estimated through low-pass adaptive filtering. A 5th order Butterworth filter,L(x,y), with a typical cutoff wavelength of 800 m was used as a narrow low-pass filter response. The adaptive part of the combined filter determines the pass band based on the dominant frequencies present in the phase of the pixels themselves [20]. The response is calculated as
(3)H(x,y)=|Z(x,y)|,
where Z is the smoothed intensity of the 2D FFT. Combine H(x,y) with L(x,y) to form the new filter response,
(4)G(x,y)=L(x,y)+β(H(x,y)H-(x,y)-1)α,
where H-(x,y) is the median value of H(x,y). α (typically 1) and β (typically 0.3) are adjustable weighting parameters.Subtracting the filtered phase valueΦ-int,x,i from Φint,x,i leaves
(5)Φint,x,i-Φ-int,x,i=ΔΦϵ,x,i+Φn,x,i-Φ-n,x,i′.The phase error caused by DEM inaccuracies is proportional to the perpendicular baseline,B⊥, so
(6)ΔΦϵ,x,i=B⊥,x,iKϵ,x,
where Kϵ is a proportionality constant. Substituting this expression into (6) gives
(7)Φint,x,i-Φ-int,x,i=B⊥,x,iKϵ,x+Φn,x,i-Φ-n,x,i′.Using all the available interferograms we can estimateKϵ for pixel x in a least square sense, as this is the only term that correlates with baseline.A measure based on the temporal coherence of pixelx is defined as
(8)γx=1N|∑i=1Nexp{j(Φint,x,i-Φ-int,x,i-ΔΦ^ε,x,i)}|,
where N is the number of available interferograms and ΔΦ^ε,x,i is the estimate of Φ-int,x,i. Assuming that Φ-n,x,i′ values are small,γx is a measure of the phase stability of the pixel and hence an indicator of whether the pixel is a PS.


#### 2.2.3. PS Probability

Theγx values are binned and normalized to get a probability density. All pixels are divided into two sets: PS pixels and non-PS pixels. p(γx) is, then, a weighted sum of the probability density for the PS pixels, pPS(γx), and the probability density for the non-PS pixels, pR(γx); that is,
(9)p(γx)=αpPS(γx)+(1-α)pR(γx),
where 0≤α≤1. To derive pR(γx), 106 pseudo-pixels with random phase are simulated; that is, W{ψx,i-ψ~x,i}=exp(jR), where R is a random variable in the interval [-π,π], and follow the steps described above to arrive at a value of γx for each pseudo-pixel. We bin these values and normalize the distribution to obtain an estimate for pR(γx). For low values of γx, that is, <0.3, pPS(γx)≈0 which implies that
(10)∫00.3p(γx)dγx=(1-α)∫00.3pR(γx)dγx.We use the data to evaluate the integral on the left and the simulation to evaluate the integral on the right. Thus, we are able to estimate a conservative value ofα. For pixel x, the probability that it is a PS is
(11)P(x∈PS)=1-(1-α)pR(γx)p(γx).


#### 2.2.4. PS Selection

After every iteration, the newγx is compared to the old value to see the changes of γx. The root-mean-square of the change are calculated. Once the RMS falls below the RMS threshold, for example, 0.005, the solution should have converged. The selection of PS pixels is based on the PS probability [21]. The probability can be calculated more accurately by considering the amplitude dispersion of the pixels, D^A,x, as well as γx. Pixels are binned by D^A,x, ensuring that there are at least 104 pixels in every bin, resulting in a number of data probability distributions, p(γx,D^A,x). For each distribution,p(γx,D^A,x), we estimate α(D^A,x). If only pixels with γx above a threshold value, γthresh(D^A,x), are selected, the number of those pixels that are non-PS pixels is given by
(12)(1-α(D^A,x))∫γthresh1PR(γx)dγx.We chooseγthresh(D^A,x) such that the fraction of non-PS pixels to the total number of chosen pixels is acceptable for our particular application; that is,
(13)(1-α(D^A,x))∫γthresh1PR(γx)dγx∫γthresh1P(γx,D^A,x)dγx=q,
where q is the maximum fraction of all the selected pixels that we will accept to be non-PS pixels (false positives). We expect γx to decrease with increasing D^A,x. This implies that as D^A,x increases, P(γx,D^A,x) will skew more to lower values of γx. The net effect on γthresh(D^A,x) is to increase with increasing D^A. Empirically, we find that the relationship is approximately linear; that is, γthresh=κD^A, where κis a constant. We find the best-fitting κby least-squares inversion and select pixels for which γx>κD^A,x as PS.


## 2.1. GPS Multiantenna Monitoring System

GPS is considered as one of the best techniques for deformation monitoring thanks to its higher automation and lower labor consuming than that of traditional geodetic survey techniques [5–8]. Depending on density of the GPS antenna, the spatial resolution is low. Therefore, GPS does have disadvantages, with the major drawback being the high cost associated with placing a permanent GPS receiver at each monitoring point.The approach of linking a single GPS receiver with multiple antennas mounted at several monitoring points has been implemented. We use a dedicated electronic timed switching device, the GPS multiantenna switch (GMAS), to connect one receiver with individual antennas. This significantly reduces the number of receivers required, thus decreasing hardware investment. GMAS has been patented in China (Patent no: 00219891.6, Owners: Xiufeng He, Xiaoli Ding, Yongrong Sun, Yongqi Chen et al.). The receiver conducts standard pseudo-range and carrier-phase observations for each antenna [9–12]. It takes efficient advantage of the characteristics of deformation survey, such as the position of a surveyed point, which changes relatively little between two consecutive survey epochs. Figure 1 shows the outline of multiantenna GPS system.Figure 1
Outline of multiantenna GPS system for deformation monitoring.


## 2.2. Persistent Scatterers InSAR Time Series

InSAR can provide measurements of surface movements at a considerably improved spatial resolution with quite high accuracy, that is, centimeter or even millimeter level, over the whole area. Thus, repeat-pass satellite InSAR is a rapidly evolving remote sensing technology for observing the Earth surface and now well developed with many applications [13–18].Due to the slow surface displacements rate, phase ambiguity problems related to temporal and geometrical decorrelations, and strong tropospheric effects, it is difficult to obtain reliable measurements of subsidence in densely vegetated areas using conventional InSAR techniques. In the late 1990s it is noticed that some radar targets maintain stable backscattering characteristics for a period of month or years, and the phase information from these stable targets (hereafter called Persistent Scatterers, PS) can be used, even over a long time period. This led to the development of PS InSAR methodologies, adopting both amplitude stability and coherence stability (i.e., correlation) as pixel selection criteria [19–21]. The persistent scatterers technique is a powerful tool that overcomes the above-said limitation of routine InSAR by simultaneously exploiting all the available SAR images gathered during repeated satellite passes.The deformation was modeled in time and a temporal model of evolution was fitted to the double difference phase during the process as the phase is unwrapped. By taking the phase difference between nearby persistent scatterer candidates, atmospheric and orbit errors are mainly reduced. Besides, signals lacking in consistency in both space and time dimension are always considered as the atmospheric error rather than the deformation signal. By high-pass filtering the unwrapped phase in time and low-pass filtering in space we estimate the spatially correlated errors. Deformation and DEM error contributing to the double difference phase are later modeled for the whole time series. The residuals between the model and the double phase signal were considered as the noise level. Then for each persistent scatterer candidates, network adjustment is applied to estimating the noise.Stanford method for persistent scatterers (StaMPS) is used for analyzing crustal deformation in nonurban environments in this paper. Unlike Ferretti method, StaMPS produces a time series of deformation, with no prior assumptions about the temporal nature of deformation. Therefore, it can get better performance and be more reliable under the circumstances of unknown or not readily parameterized deformation. Here only PS identification and selection steps of the innovative PS processing method are discussed; more detailed information about the algorithm, such as the phase unwrapping and spatially correlated terms removal, can be found in [22].


### 2.2.1. Amplitude Analysis

Ferretti et al. [19] defined the amplitude dispersion index as
(1)DA≡σAmA,
where σA and mA are, respectively, the standard deviation and the mean of the amplitude values. For a constant signal and high signal to noise ratio (SNR), DA≈σΦ, where σΦ is the phase standard deviation. As lower DA indicates higher phase stability, consideration of amplitude is useful both to reduce the number of pixels for phase analysis and to better estimate the probability of a pixel being a PS. The pixels selected as PS candidates have DA≤0.4 [23].


### 2.2.2. Phase Analysis

The residual phase (with flat earth and DEM reference phase removal),Φint, of the xth pixel in the ith interferogram can be written as
(2)Φint,x,i=Φdef,x,i+ΔΦε,x,i+Φatm,x,i+ΔΦorb,x,i+Φn,x,i,
where Φdef is the phase change due to the deformation of the pixel in the satellite line-of-sight (LOS) direction, ΔΦε is the residual topographic phase due to DEM error, Φatm is the phase equivalent of the difference in atmospheric retardation between passes, ΔΦorb is the residual phase due to orbit inaccuracies, and Φn is the noise term due to variability in scattering from the pixel, thermal noise, and coregistration errors.The spatially correlated phase is estimated through low-pass adaptive filtering. A 5th order Butterworth filter,L(x,y), with a typical cutoff wavelength of 800 m was used as a narrow low-pass filter response. The adaptive part of the combined filter determines the pass band based on the dominant frequencies present in the phase of the pixels themselves [20]. The response is calculated as
(3)H(x,y)=|Z(x,y)|,
where Z is the smoothed intensity of the 2D FFT. Combine H(x,y) with L(x,y) to form the new filter response,
(4)G(x,y)=L(x,y)+β(H(x,y)H-(x,y)-1)α,
where H-(x,y) is the median value of H(x,y). α (typically 1) and β (typically 0.3) are adjustable weighting parameters.Subtracting the filtered phase valueΦ-int,x,i from Φint,x,i leaves
(5)Φint,x,i-Φ-int,x,i=ΔΦϵ,x,i+Φn,x,i-Φ-n,x,i′.The phase error caused by DEM inaccuracies is proportional to the perpendicular baseline,B⊥, so
(6)ΔΦϵ,x,i=B⊥,x,iKϵ,x,
where Kϵ is a proportionality constant. Substituting this expression into (6) gives
(7)Φint,x,i-Φ-int,x,i=B⊥,x,iKϵ,x+Φn,x,i-Φ-n,x,i′.Using all the available interferograms we can estimateKϵ for pixel x in a least square sense, as this is the only term that correlates with baseline.A measure based on the temporal coherence of pixelx is defined as
(8)γx=1N|∑i=1Nexp{j(Φint,x,i-Φ-int,x,i-ΔΦ^ε,x,i)}|,
where N is the number of available interferograms and ΔΦ^ε,x,i is the estimate of Φ-int,x,i. Assuming that Φ-n,x,i′ values are small,γx is a measure of the phase stability of the pixel and hence an indicator of whether the pixel is a PS.


### 2.2.3. PS Probability

Theγx values are binned and normalized to get a probability density. All pixels are divided into two sets: PS pixels and non-PS pixels. p(γx) is, then, a weighted sum of the probability density for the PS pixels, pPS(γx), and the probability density for the non-PS pixels, pR(γx); that is,
(9)p(γx)=αpPS(γx)+(1-α)pR(γx),
where 0≤α≤1. To derive pR(γx), 106 pseudo-pixels with random phase are simulated; that is, W{ψx,i-ψ~x,i}=exp(jR), where R is a random variable in the interval [-π,π], and follow the steps described above to arrive at a value of γx for each pseudo-pixel. We bin these values and normalize the distribution to obtain an estimate for pR(γx). For low values of γx, that is, <0.3, pPS(γx)≈0 which implies that
(10)∫00.3p(γx)dγx=(1-α)∫00.3pR(γx)dγx.We use the data to evaluate the integral on the left and the simulation to evaluate the integral on the right. Thus, we are able to estimate a conservative value ofα. For pixel x, the probability that it is a PS is
(11)P(x∈PS)=1-(1-α)pR(γx)p(γx).


### 2.2.4. PS Selection

After every iteration, the newγx is compared to the old value to see the changes of γx. The root-mean-square of the change are calculated. Once the RMS falls below the RMS threshold, for example, 0.005, the solution should have converged. The selection of PS pixels is based on the PS probability [21]. The probability can be calculated more accurately by considering the amplitude dispersion of the pixels, D^A,x, as well as γx. Pixels are binned by D^A,x, ensuring that there are at least 104 pixels in every bin, resulting in a number of data probability distributions, p(γx,D^A,x). For each distribution,p(γx,D^A,x), we estimate α(D^A,x). If only pixels with γx above a threshold value, γthresh(D^A,x), are selected, the number of those pixels that are non-PS pixels is given by
(12)(1-α(D^A,x))∫γthresh1PR(γx)dγx.We chooseγthresh(D^A,x) such that the fraction of non-PS pixels to the total number of chosen pixels is acceptable for our particular application; that is,
(13)(1-α(D^A,x))∫γthresh1PR(γx)dγx∫γthresh1P(γx,D^A,x)dγx=q,
where q is the maximum fraction of all the selected pixels that we will accept to be non-PS pixels (false positives). We expect γx to decrease with increasing D^A,x. This implies that as D^A,x increases, P(γx,D^A,x) will skew more to lower values of γx. The net effect on γthresh(D^A,x) is to increase with increasing D^A. Empirically, we find that the relationship is approximately linear; that is, γthresh=κD^A, where κis a constant. We find the best-fitting κby least-squares inversion and select pixels for which γx>κD^A,x as PS.


## 2.2.1. Amplitude Analysis

Ferretti et al. [19] defined the amplitude dispersion index as
(1)DA≡σAmA,
where σA and mA are, respectively, the standard deviation and the mean of the amplitude values. For a constant signal and high signal to noise ratio (SNR), DA≈σΦ, where σΦ is the phase standard deviation. As lower DA indicates higher phase stability, consideration of amplitude is useful both to reduce the number of pixels for phase analysis and to better estimate the probability of a pixel being a PS. The pixels selected as PS candidates have DA≤0.4 [23].


## 2.2.2. Phase Analysis

The residual phase (with flat earth and DEM reference phase removal),Φint, of the xth pixel in the ith interferogram can be written as
(2)Φint,x,i=Φdef,x,i+ΔΦε,x,i+Φatm,x,i+ΔΦorb,x,i+Φn,x,i,
where Φdef is the phase change due to the deformation of the pixel in the satellite line-of-sight (LOS) direction, ΔΦε is the residual topographic phase due to DEM error, Φatm is the phase equivalent of the difference in atmospheric retardation between passes, ΔΦorb is the residual phase due to orbit inaccuracies, and Φn is the noise term due to variability in scattering from the pixel, thermal noise, and coregistration errors.The spatially correlated phase is estimated through low-pass adaptive filtering. A 5th order Butterworth filter,L(x,y), with a typical cutoff wavelength of 800 m was used as a narrow low-pass filter response. The adaptive part of the combined filter determines the pass band based on the dominant frequencies present in the phase of the pixels themselves [20]. The response is calculated as
(3)H(x,y)=|Z(x,y)|,
where Z is the smoothed intensity of the 2D FFT. Combine H(x,y) with L(x,y) to form the new filter response,
(4)G(x,y)=L(x,y)+β(H(x,y)H-(x,y)-1)α,
where H-(x,y) is the median value of H(x,y). α (typically 1) and β (typically 0.3) are adjustable weighting parameters.Subtracting the filtered phase valueΦ-int,x,i from Φint,x,i leaves
(5)Φint,x,i-Φ-int,x,i=ΔΦϵ,x,i+Φn,x,i-Φ-n,x,i′.The phase error caused by DEM inaccuracies is proportional to the perpendicular baseline,B⊥, so
(6)ΔΦϵ,x,i=B⊥,x,iKϵ,x,
where Kϵ is a proportionality constant. Substituting this expression into (6) gives
(7)Φint,x,i-Φ-int,x,i=B⊥,x,iKϵ,x+Φn,x,i-Φ-n,x,i′.Using all the available interferograms we can estimateKϵ for pixel x in a least square sense, as this is the only term that correlates with baseline.A measure based on the temporal coherence of pixelx is defined as
(8)γx=1N|∑i=1Nexp{j(Φint,x,i-Φ-int,x,i-ΔΦ^ε,x,i)}|,
where N is the number of available interferograms and ΔΦ^ε,x,i is the estimate of Φ-int,x,i. Assuming that Φ-n,x,i′ values are small,γx is a measure of the phase stability of the pixel and hence an indicator of whether the pixel is a PS.


## 2.2.3. PS Probability

Theγx values are binned and normalized to get a probability density. All pixels are divided into two sets: PS pixels and non-PS pixels. p(γx) is, then, a weighted sum of the probability density for the PS pixels, pPS(γx), and the probability density for the non-PS pixels, pR(γx); that is,
(9)p(γx)=αpPS(γx)+(1-α)pR(γx),
where 0≤α≤1. To derive pR(γx), 106 pseudo-pixels with random phase are simulated; that is, W{ψx,i-ψ~x,i}=exp(jR), where R is a random variable in the interval [-π,π], and follow the steps described above to arrive at a value of γx for each pseudo-pixel. We bin these values and normalize the distribution to obtain an estimate for pR(γx). For low values of γx, that is, <0.3, pPS(γx)≈0 which implies that
(10)∫00.3p(γx)dγx=(1-α)∫00.3pR(γx)dγx.We use the data to evaluate the integral on the left and the simulation to evaluate the integral on the right. Thus, we are able to estimate a conservative value ofα. For pixel x, the probability that it is a PS is
(11)P(x∈PS)=1-(1-α)pR(γx)p(γx).


## 2.2.4. PS Selection

After every iteration, the newγx is compared to the old value to see the changes of γx. The root-mean-square of the change are calculated. Once the RMS falls below the RMS threshold, for example, 0.005, the solution should have converged. The selection of PS pixels is based on the PS probability [21]. The probability can be calculated more accurately by considering the amplitude dispersion of the pixels, D^A,x, as well as γx. Pixels are binned by D^A,x, ensuring that there are at least 104 pixels in every bin, resulting in a number of data probability distributions, p(γx,D^A,x). For each distribution,p(γx,D^A,x), we estimate α(D^A,x). If only pixels with γx above a threshold value, γthresh(D^A,x), are selected, the number of those pixels that are non-PS pixels is given by
(12)(1-α(D^A,x))∫γthresh1PR(γx)dγx.We chooseγthresh(D^A,x) such that the fraction of non-PS pixels to the total number of chosen pixels is acceptable for our particular application; that is,
(13)(1-α(D^A,x))∫γthresh1PR(γx)dγx∫γthresh1P(γx,D^A,x)dγx=q,
where q is the maximum fraction of all the selected pixels that we will accept to be non-PS pixels (false positives). We expect γx to decrease with increasing D^A,x. This implies that as D^A,x increases, P(γx,D^A,x) will skew more to lower values of γx. The net effect on γthresh(D^A,x) is to increase with increasing D^A. Empirically, we find that the relationship is approximately linear; that is, γthresh=κD^A, where κis a constant. We find the best-fitting κby least-squares inversion and select pixels for which γx>κD^A,x as PS.


## 3. Study Area

### 3.1. Site Conditions

Three ancient landslides, namely, Luanshigang, Futang, and Kangjiaping were identified by geologists on Luobogang hillock. The stability of the hill slopes is adversely affected by the water surface rising and fluctuation in the reservoir. Large-scale catastrophic movements may be triggered when the shear stress exceeds the shear resistance of the material.Hundreds of retaining walls were built to minimize large movement of the landslides. But cracks were soon observed on the retaining walls (Figure2). The displacements of landslides cause continuous damage to buildings and infrastructures.Figure 2
Cracks found on the retaining walls in Hanyuan.


### 3.2. SAR Data Set

The present study used a data set composed of 34 images produced by ENVISAT/ASAR systems spanning the time from September 20, 2003 to August 14, 2010. The ENVISAT ASAR data (descending orbits track 61) were provided under the European Space Agency (ESA)-National Remote Sensing Center of China (NRSCC) Dragon cooperation project, collected with a nominal radar look angle of 23°. DORIS precise orbits data provided by the ESRIN help desk of ESA were applied to calculating ENVISAT interferometric baselines. The 3 arcsecond (~90 m) shuttle radar topography mission (SRTM) DEM data were used to correct the interferometric phase for topography at the first step and then for geocoding the C-band SAR data (transforming Range-Doppler coordinates into Universal Transverse Mercator map geometry system).The size of the test area is about 12 * 12 km, covering about 29.26° to 29.40°N and 102.53° to 102.72°E (Figure3). The scene acquired on May 31, 2008 is selected as the “master” image. The maximum temporal baseline is 1715 days, and the maximum normal baseline is about 746 m. Figure 4 and Table 1 show more details about the temporal and perpendicular baselines.Table 1
ENVISAT ASAR data set used in the PS time series analysis.
Date
Temporal baseline (day)
Perp. baseline (m)
Date
Temporal baseline (day)
Perp. baseline (m)
20-Sep-2003
−1715
354.5
13-Sep-2008
105
351.5
25-Oct-2003
−1680
644.8
27-Dec-2008
210
−28.4
29-Nov-2003
−1645
172.3
31-Jan-2009
245
230.4
3-Jan-2004
−1610
745.9
7-Mar-2009
280
351.8
13-Mar-2004
−1540
716.9
11-Apr-2009
315
−66.6
17-Apr-2004
−1505
416.8
16-May-2009
350
84.2
26-Jun-2004
−1435
373.2
20-Jun-2009
385
403.7
13-Nov-2004
−1295
682.7
25-Jul-2009
420
−58.7
18-Dec-2004
−1260
184.3
3-Oct-2009
490
−88.5
26-Feb-2005
−1190
118.1
12-Dec-2009
560
398.7
18-Nov-2006
−560
682.8
16-Jan-2010
595
145.8
27-Jan-2007
−490
369.7
20-Feb-2010
630
118.6
21-Jul-2007
−315
−5.6
27-Mar-2010
665
277.2
3-Nov-2007
−210
251.4
1-May-2010
700
193.3
31-May-2008
0
0
5-Jun-2010
735
232.5
5-Jul-2008
35
246.8
10-Jul-2010
770
−43
9-Aug-2008
70
257.4
14-Aug-2010
805
−5.3Figure 3
Location map of Hanyuan County, Sichuan. The bigger black rectangle indicates the ENVISAT scene footprint (descending orbits track 61) used in the PS analysis and the smaller one represents the region of interest in this paper.Figure 4
Perpendicular baseline relative to the “master” scene plot for the interferometric stack.


## 3.1. Site Conditions

Three ancient landslides, namely, Luanshigang, Futang, and Kangjiaping were identified by geologists on Luobogang hillock. The stability of the hill slopes is adversely affected by the water surface rising and fluctuation in the reservoir. Large-scale catastrophic movements may be triggered when the shear stress exceeds the shear resistance of the material.Hundreds of retaining walls were built to minimize large movement of the landslides. But cracks were soon observed on the retaining walls (Figure2). The displacements of landslides cause continuous damage to buildings and infrastructures.Figure 2
Cracks found on the retaining walls in Hanyuan.


## 3.2. SAR Data Set

The present study used a data set composed of 34 images produced by ENVISAT/ASAR systems spanning the time from September 20, 2003 to August 14, 2010. The ENVISAT ASAR data (descending orbits track 61) were provided under the European Space Agency (ESA)-National Remote Sensing Center of China (NRSCC) Dragon cooperation project, collected with a nominal radar look angle of 23°. DORIS precise orbits data provided by the ESRIN help desk of ESA were applied to calculating ENVISAT interferometric baselines. The 3 arcsecond (~90 m) shuttle radar topography mission (SRTM) DEM data were used to correct the interferometric phase for topography at the first step and then for geocoding the C-band SAR data (transforming Range-Doppler coordinates into Universal Transverse Mercator map geometry system).The size of the test area is about 12 * 12 km, covering about 29.26° to 29.40°N and 102.53° to 102.72°E (Figure3). The scene acquired on May 31, 2008 is selected as the “master” image. The maximum temporal baseline is 1715 days, and the maximum normal baseline is about 746 m. Figure 4 and Table 1 show more details about the temporal and perpendicular baselines.Table 1
ENVISAT ASAR data set used in the PS time series analysis.
Date
Temporal baseline (day)
Perp. baseline (m)
Date
Temporal baseline (day)
Perp. baseline (m)
20-Sep-2003
−1715
354.5
13-Sep-2008
105
351.5
25-Oct-2003
−1680
644.8
27-Dec-2008
210
−28.4
29-Nov-2003
−1645
172.3
31-Jan-2009
245
230.4
3-Jan-2004
−1610
745.9
7-Mar-2009
280
351.8
13-Mar-2004
−1540
716.9
11-Apr-2009
315
−66.6
17-Apr-2004
−1505
416.8
16-May-2009
350
84.2
26-Jun-2004
−1435
373.2
20-Jun-2009
385
403.7
13-Nov-2004
−1295
682.7
25-Jul-2009
420
−58.7
18-Dec-2004
−1260
184.3
3-Oct-2009
490
−88.5
26-Feb-2005
−1190
118.1
12-Dec-2009
560
398.7
18-Nov-2006
−560
682.8
16-Jan-2010
595
145.8
27-Jan-2007
−490
369.7
20-Feb-2010
630
118.6
21-Jul-2007
−315
−5.6
27-Mar-2010
665
277.2
3-Nov-2007
−210
251.4
1-May-2010
700
193.3
31-May-2008
0
0
5-Jun-2010
735
232.5
5-Jul-2008
35
246.8
10-Jul-2010
770
−43
9-Aug-2008
70
257.4
14-Aug-2010
805
−5.3Figure 3
Location map of Hanyuan County, Sichuan. The bigger black rectangle indicates the ENVISAT scene footprint (descending orbits track 61) used in the PS analysis and the smaller one represents the region of interest in this paper.Figure 4
Perpendicular baseline relative to the “master” scene plot for the interferometric stack.


## 4. Results and Analysis

### 4.1. GPS Monitoring System

More than 100 monitoring stations are installed in the whole county. The 28 continuous monitoring stations operate 7 * 24 hours per week and engineers will periodically check the other 90 stations. In key areas with more buildings, the density of GPS monitoring stations is also higher. Detailed information about the GPS monitoring network can be found in [24].The power to the GPS receivers and data transmission devices is supplied by the uninterruptible power system (UPS) at the monitoring sites. Although solar panels are all employed at each site, power from the city electricity grid is always the preferred option when it is available.Each site needs to transfer massive amount of data to the control center which is located within a town area. A public wireless phone network or other media may be used. In this case, we chose the wireless transmission and 3rd generation telecommunication technologies.The central servers housed inside the control center are connected with a local area network (LAN). The control center serves two main functions. Firstly, it receives the monitoring data, verifies the integrity of the data, assesses the status of all the devices, and issues control commands to the field stations. Secondly, the control center is responsible for data analysis, storage of processed data, deformation modeling, forecast, and transmission of information.


### 4.2. GPS versus Digital Inclinometer in a Landslide Case

The GPS monitoring system is designed to detect mm level movements. Figure5 illustrates the movements on the three components captured by the system in August 2011. The result for site TP6-1 showed a significant displacement which occurred in August, and it became stable again in September. The site moves by almost 15 mm to the north as well as 5 mm to the west within 20 days.Figure 5
The movements captured by the landslide monitoring system during August 2011.The engineers periodically collect subsurface deformation data using RST MEMS digital inclinometer system at GPS monitoring stations. Measurements were taken along two perpendicular directions, bothA and B axes, from the depth of 0.5 m to 28.5 m at 0.5 m intervals. Figure 6 shows the comparison between the GPS and digital inclinometer results on TP6-1 during the unstable period. The total displacements at the depth of 0.5 m from inclinometer (which is the measurement point closest to the ground surface) were compared with the results obtained from GPS. Generally, the displacement amplitude on the ground surface will be larger than that in subsurface. Therefore, the magnitude from digital inclinometer is smaller than that from GPS monitoring. The displacement trend and time of occurrence captured by both methods match well.Figure 6
GPS versus digital inclinometer results for the displacement on TP6-1 in August 2011.Field investigations were conducted shortly after the displacement was discovered. Construction activities near site TP6-1 influenced the stability of the local ground. A mound of soil heaped in the south results in much more load and pressure to the retaining wall. We believe that this is the main reason for the monitoring site moving north.Due to page limit, detailed information about the GPS monitoring system established in Hanyuan and the result time series is omitted here in this paper. The exhaustive report and analysis can be found in [24].


### 4.3. Time Series InSAR Results

Pubugou Hydropower Station started to impound water for the first time on November 1, 2009. The old county no longer existed from then on. In the mean SAR amplitude image (Figure7), we can clearly see the Dadu River and the Luobogang resettlement zone. A reservoir, named Hanyuan Lake, was formed at the location of the old Hanyuan County.Figure 7
The new resettlement county, Luobogang hillock, in the SAR image. After completing the hydropower project and water storage, a reservoir was formed and the old Hanyuan County is now under water.With StaMPS, the density of PS identified over most of the region was about 20 per km2 (1794 PSs in total, excluding the water area). No PSs were identified in the densely vegetated area due to the seasonal cultivation, flood irrigation, and crop growth. Even in some parts of urban area (on Luobogang), relatively rare PSs were identified because of the construction activities in the resettlement zone which lead to local interferometric phase decorrelation. In addition, due to the side-looking radar signals, the slope and aspect in the mountain areas also have considerable influence on PS detection.Unlike other PS algorithms, StaMPS extracts the LOS deformation via phase unwrapping rather than assuming a model for how displacement rates vary with time. Figure8 demonstrates the LOS deformation time series at Hanyuan during the period from 2003 to 2010. For the scene predating the “master” image, positive phase implies movement towards the satellite, and vice versa.Figure 8
Time series of LOS deformation at Hanyuan spanning from September 2003 to August 2010, referenced temporally to the “master” scene (May 31, 2008) and superimposed on the mean amplitude image.Based on the unwrapped phase time series, mean LOS velocity of PSs was calculated. As shown in Figure9, the density of PS in the southern part of the new county reached more than 40 per km2. Due to the construction activities, rare PSs were detected in downtown. Scatterers appearing on all the images were accepted as PS and the ones existed in only a subset of the data stack were rejected. Therefore, no scatterers were identified in the old county which disappeared after the water storage. Standard deviation of the mean LOS velocity shows the precision of the linear regression.Mean LOS deformation velocity map of the PSs and the standard deviation of mean LOS velocity. The units are mm/year with positive values being towards the satellite along the LOS direction.
(a)
(b)Overall, the mean LOS velocity of the PSs is relatively slow. Several PSs in two typical regions, displayed in Figure10, are selected for further analysis.Figure 10
LOS deformation time series of the selected PSs. The red lines indicate linear trends.The deformation histories of the PSs (4 PSs in Figure10 upper case and 5 PSs in the lower case) in such a small region match well. This, to a certain degree, verified the reliability of the phase unwrapping. The average velocities of the adjacent PSs in the two cases are about 3 mm/year and −4 mm/year, respectively. It should be noted that the sign of the deformation and the velocity implies the movement towards or away from the satellite on the LOS direction. Without the information of the look angle and aspects of the slope, we cannot determine the movement direction (uplift or subsidence) only from the InSAR results. This is also the main limitation of the InSAR technique, which can be improved by integrating ascending/descending orbits data from multiplatform satellite sensors.The real-time GPS landslide monitoring system of new Hanyuan county began trial operation in May 2011. The ENVISAT satellite was brought to a new lower orbit in October 2010. The orbital change ends the PS time series analysis using ASAR data. Unfortunately, no synchronous data can be collected from the two sources; thus, verification between GPS and InSAR is impossible.


## 4.1. GPS Monitoring System

More than 100 monitoring stations are installed in the whole county. The 28 continuous monitoring stations operate 7 * 24 hours per week and engineers will periodically check the other 90 stations. In key areas with more buildings, the density of GPS monitoring stations is also higher. Detailed information about the GPS monitoring network can be found in [24].The power to the GPS receivers and data transmission devices is supplied by the uninterruptible power system (UPS) at the monitoring sites. Although solar panels are all employed at each site, power from the city electricity grid is always the preferred option when it is available.Each site needs to transfer massive amount of data to the control center which is located within a town area. A public wireless phone network or other media may be used. In this case, we chose the wireless transmission and 3rd generation telecommunication technologies.The central servers housed inside the control center are connected with a local area network (LAN). The control center serves two main functions. Firstly, it receives the monitoring data, verifies the integrity of the data, assesses the status of all the devices, and issues control commands to the field stations. Secondly, the control center is responsible for data analysis, storage of processed data, deformation modeling, forecast, and transmission of information.


## 4.2. GPS versus Digital Inclinometer in a Landslide Case

The GPS monitoring system is designed to detect mm level movements. Figure5 illustrates the movements on the three components captured by the system in August 2011. The result for site TP6-1 showed a significant displacement which occurred in August, and it became stable again in September. The site moves by almost 15 mm to the north as well as 5 mm to the west within 20 days.Figure 5
The movements captured by the landslide monitoring system during August 2011.The engineers periodically collect subsurface deformation data using RST MEMS digital inclinometer system at GPS monitoring stations. Measurements were taken along two perpendicular directions, bothA and B axes, from the depth of 0.5 m to 28.5 m at 0.5 m intervals. Figure 6 shows the comparison between the GPS and digital inclinometer results on TP6-1 during the unstable period. The total displacements at the depth of 0.5 m from inclinometer (which is the measurement point closest to the ground surface) were compared with the results obtained from GPS. Generally, the displacement amplitude on the ground surface will be larger than that in subsurface. Therefore, the magnitude from digital inclinometer is smaller than that from GPS monitoring. The displacement trend and time of occurrence captured by both methods match well.Figure 6
GPS versus digital inclinometer results for the displacement on TP6-1 in August 2011.Field investigations were conducted shortly after the displacement was discovered. Construction activities near site TP6-1 influenced the stability of the local ground. A mound of soil heaped in the south results in much more load and pressure to the retaining wall. We believe that this is the main reason for the monitoring site moving north.Due to page limit, detailed information about the GPS monitoring system established in Hanyuan and the result time series is omitted here in this paper. The exhaustive report and analysis can be found in [24].


## 4.3. Time Series InSAR Results

Pubugou Hydropower Station started to impound water for the first time on November 1, 2009. The old county no longer existed from then on. In the mean SAR amplitude image (Figure7), we can clearly see the Dadu River and the Luobogang resettlement zone. A reservoir, named Hanyuan Lake, was formed at the location of the old Hanyuan County.Figure 7
The new resettlement county, Luobogang hillock, in the SAR image. After completing the hydropower project and water storage, a reservoir was formed and the old Hanyuan County is now under water.With StaMPS, the density of PS identified over most of the region was about 20 per km2 (1794 PSs in total, excluding the water area). No PSs were identified in the densely vegetated area due to the seasonal cultivation, flood irrigation, and crop growth. Even in some parts of urban area (on Luobogang), relatively rare PSs were identified because of the construction activities in the resettlement zone which lead to local interferometric phase decorrelation. In addition, due to the side-looking radar signals, the slope and aspect in the mountain areas also have considerable influence on PS detection.Unlike other PS algorithms, StaMPS extracts the LOS deformation via phase unwrapping rather than assuming a model for how displacement rates vary with time. Figure8 demonstrates the LOS deformation time series at Hanyuan during the period from 2003 to 2010. For the scene predating the “master” image, positive phase implies movement towards the satellite, and vice versa.Figure 8
Time series of LOS deformation at Hanyuan spanning from September 2003 to August 2010, referenced temporally to the “master” scene (May 31, 2008) and superimposed on the mean amplitude image.Based on the unwrapped phase time series, mean LOS velocity of PSs was calculated. As shown in Figure9, the density of PS in the southern part of the new county reached more than 40 per km2. Due to the construction activities, rare PSs were detected in downtown. Scatterers appearing on all the images were accepted as PS and the ones existed in only a subset of the data stack were rejected. Therefore, no scatterers were identified in the old county which disappeared after the water storage. Standard deviation of the mean LOS velocity shows the precision of the linear regression.Mean LOS deformation velocity map of the PSs and the standard deviation of mean LOS velocity. The units are mm/year with positive values being towards the satellite along the LOS direction.
(a)
(b)Overall, the mean LOS velocity of the PSs is relatively slow. Several PSs in two typical regions, displayed in Figure10, are selected for further analysis.Figure 10
LOS deformation time series of the selected PSs. The red lines indicate linear trends.The deformation histories of the PSs (4 PSs in Figure10 upper case and 5 PSs in the lower case) in such a small region match well. This, to a certain degree, verified the reliability of the phase unwrapping. The average velocities of the adjacent PSs in the two cases are about 3 mm/year and −4 mm/year, respectively. It should be noted that the sign of the deformation and the velocity implies the movement towards or away from the satellite on the LOS direction. Without the information of the look angle and aspects of the slope, we cannot determine the movement direction (uplift or subsidence) only from the InSAR results. This is also the main limitation of the InSAR technique, which can be improved by integrating ascending/descending orbits data from multiplatform satellite sensors.The real-time GPS landslide monitoring system of new Hanyuan county began trial operation in May 2011. The ENVISAT satellite was brought to a new lower orbit in October 2010. The orbital change ends the PS time series analysis using ASAR data. Unfortunately, no synchronous data can be collected from the two sources; thus, verification between GPS and InSAR is impossible.


## 5. Concluding Remarks

The measurement of ground surface movements over Hanyuan is one of the important field monitoring activities which would provide advance warning against sudden failures. A real-time landslide monitoring system using GPS in new Hanyuan County is now in normal operation. The comparison between GPS and digital inclinometer results shows a high consistency in respect of the displacement trend. This suggests that the GPS monitoring system can be employed as a complement to the traditional ground movement monitoring methods.New Hanyuan County is not an ideal environment for the application of conventional InSAR technique because of its dense vegetation cover, high soil moisture content and humidity, and the active constructions over the past years. However, the persistent scatterer InSAR approach can overcome the problems, and the time series of surface displacements show various patterns and magnitudes of deformation in the resettlement zone.GPS provides 3D information at relatively sparse observation points while InSAR captures LOS deformation on a surface (in effect, as for PS approach, the “surface” refers to denser points) with wide coverage. Therefore, it is obvious that the two techniques are complementary. Unfortunately, neither ENVISAT ASAR nor ALOS PALSAR satellite data can be synchronous with the ground GPS monitoring results. Nevertheless, the historical information obtained via InSAR technique using archived scenes cannot be replicated.



---

*Source: 601209-2013-12-29.xml*