Monitoring of Land-Surface Deformation in the Karamay Oilﬁeld, Xinjiang, China, Using SAR Interferometry

: Synthetic Aperture Radar (SAR) interferometry is a technique that provides high-resolution measurements of the ground displacement associated with various geophysical processes. To investigate the land-surface deformation in Karamay, a typical oil-producing city in the Xinjiang Uyghur Autonomous Region, China, Advanced Land Observing Satellite (ALOS) Phased Array L-band Synthetic Aperture Radar (PALSAR) data were acquired for the period from 2007 to 2009, and a two-pass differential SAR interferometry (D-InSAR) process was applied. The experimental results showed that two sites in the north-eastern part of the city exhibit a clear indication of land deformation. For a further evaluation of the D-InSAR result, the Persistent Scatterer (PS) and Small Baseline Subset (SBAS)-InSAR techniques were applied for 21 time series Environmental Satellite (ENVISAT) C-band Advanced Synthetic Aperture Radar (ASAR) data from 2003 to 2010. The comparison between the D-InSAR and SBAS-InSAR measurements had better agreement than that from the PS-InSAR measurement. The maximum deformation rate attributed to subsurface water injection for the period from 2003 to 2010 was up to approximately 33 mm/year in the line of sight (LOS) direction. The interferometric phase change from November 2007 to June 2010 showed a clear deformation pattern, and the rebound center has been expanding in scale and increasing in quantity.


Introduction
Interferometric Synthetic Aperture Radar (InSAR) is a proven remote sensing technique that uses the phase information of SAR images to measure ground surface movements. However, the conventional InSAR method is limited by the so-called temporal and geometrical decorrelation as well as by atmospheric artifacts [1]. Over the years, new advanced InSAR techniques have been proposed to improve the precision of deformation measurement by the joint analysis of a set of SAR images [2,3], for example, Permanent Scatterer InSAR (PS-InSAR™)-use of time series SAR images to detect coherent radar signals over a sequence of interferograms [4,5]; Small Baseline Subset (SBAS)-use of all available SAR images with small baselines to achieve a high degree of spatial coverage of distributed scatters [6]; Persistent Scatterer Pairs (PSP)-works with pairs of points to identify and analyze persistent scatterer [7]; Quasi Persistent Scatterers (QPS)-utilizes partially coherent targets to increase the spatial density of the observations [8]; Stable Points Network (SPN), which has three key features (i.e., pixel selection, use of multi-master images and modeling capability) and less sensitive to geometric decorrelation [9]; SqueeSAR™-a second generation of images were as follows: reference system for planning (RSP) number 94; path number 501; and acquired in the ascending orbit with an off-nadir angle of 34.3°. We also used twenty-one C-band ENVISAT-ASAR images provided by the European Space Agency (ESA) and acquired in the period from 30 September 2003 to 15 June 2010. These data sets were acquired in the descending orbit with an incidence angle of 22.9°. A Shuttle Radar Topography Mission-3 (SRTM-3) version-4 (90 m resolution) was used to eliminate the topographic phase, which was downloaded from (http://www.cgiar-csi.org/). Figure 1 shows the cover ranges of the ALOS and ENVISAT data.

D-InSAR
D-InSAR is a technique capable of detecting land-surface deformation by analyzing a single interferogram that is derived from a pair of SAR images with the addition of a DEM [38]. In this study, the two-pass interferometry method was implemented by using two ALOS-PALSAR Single Look Complex (SLC) images for interferogram generation, and then, the topographic phase was removed using the SRTM DEM data. To remove noise and to smooth the interferogram, the Goldstein-Werner filtering process was applied [39], and the coefficient in the filtering process was (c)

Data Sets
The L-band ALOS satellite data sets covering the region of interest were provided by the Japan Aerospace Exploration Agency (JAXA). The SAR interferograms were computed from PALSAR fine-beam single-polarization (FBS) data taken on three different dates (20 January 2007, 10 December 2008, and 25 January 2009), and predominantly over the winter season to minimize the adverse impact that vegetation has on the accuracy of SAR interferometry. Observation parameters for all the images were as follows: reference system for planning (RSP) number 94; path number 501; and acquired in the ascending orbit with an off-nadir angle of 34.3 • . We also used twenty-one C-band ENVISAT-ASAR images provided by the European Space Agency (ESA) and acquired in the period from 30 September 2003 to 15 June 2010. These data sets were acquired in the descending orbit with an incidence angle of 22.9 • . A Shuttle Radar Topography Mission-3 (SRTM-3) version-4 (90 m resolution) was used to eliminate the topographic phase, which was downloaded from (http://www.cgiar-csi.org/). Figure 1 shows the cover ranges of the ALOS and ENVISAT data.

D-InSAR
D-InSAR is a technique capable of detecting land-surface deformation by analyzing a single interferogram that is derived from a pair of SAR images with the addition of a DEM [38]. In this study, the two-pass interferometry method was implemented by using two ALOS-PALSAR Single Look Complex (SLC) images for interferogram generation, and then, the topographic phase was removed using the SRTM DEM data. To remove noise and to smooth the interferogram, the Goldstein-Werner filtering process was applied [39], and the coefficient in the filtering process was set to 0.2. Finally, the InSAR products were geocoded from the Range-Doppler coordinates to the map geometry with a pixel resolution of 25 m. The SARscape ® Modules (5.1) for ENVI (5.3) software suit was employed to process the level-1.1 data and perform interferometric analyses.

PS and SBAS-InSAR
PS-InSAR is one of the promising approaches that improves the precision of conventional InSAR displacement measurements. The PS-InSAR algorithm utilizes a time series of radar images to detect coherent radar signals from PS points to derive information of the terrain motion [4,5,40]. Another algorithm called SBAS-InSAR was also employed, which utilizes all available SAR images with small baselines to get high degree of spatial coverage of distributed scatters [6]. Both algorithms can compensate the disadvantages of the conventional D-InSAR, namely, phase errors due to geometrical and temporal decorrelation as well as atmospheric disturbance [31].
The ENVISAT-ASAR data were processed using both the PS and SBAS-InSAR methods. For PS-InSAR, we selected the slant-range image (28 November 2006) as a master image and generated 20 interferograms. The perpendicular baseline ranged from 43 m to 951 m, the temporal baseline ranged from 68 days to 1277 days, and selected the potential PS candidates with a coherence threshold of 0.67. For SBAS-InSAR, we selected the slant-range image (28 October 2008) as a super master image. To reduce the geometrical and temporal decorrelation, the threshold criteria with a maximum temporal baseline of 735 days and a maximum perpendicular baseline of 483 m were used, and then, 68 interferograms were generated. Figure 2 shows the time-baseline plot for the ENVISAT-ASAR images used for PS and SBAS-InSAR processing. A combination of the minimum cost flow (MCF) network [41] and Delauney 3D method [42] was employed for phase unwrapping with an unwrapping coherence threshold of 0.35. After that, a screening process of flattened, filtered interferograms and unwrapped phases for checking some unwanted behaviors and data problems, which were caused by strong orbit inaccuracy, non-coherent pairs, atmospheric artefacts, residual topography etc., and 12 interferometric pairs were discarded from further processing. For refinement and reflattening, the reference points where the unwrapped phase value close to zero and flat areas identified from a topographic map projected in the line of sight (LOS) direction were chosen [43]. A custom atmospheric filtering was performed with a low-pass spatial filter with a 1.2 km × 1.2 km window on each single acquisition, and a high-pass filter at 365 days on these preliminary results to recover the final and cleaned displacement time series. Finally, geocoding was done in the original satellite LOS direction with a pixel resolution of 25 m. The same software package used for D-InSAR was also used in the PS and SBAS-InSAR processing chain.

PS and SBAS-InSAR
PS-InSAR is one of the promising approaches that improves the precision of conventional InSAR displacement measurements. The PS-InSAR algorithm utilizes a time series of radar images to detect coherent radar signals from PS points to derive information of the terrain motion [4,5,40]. Another algorithm called SBAS-InSAR was also employed, which utilizes all available SAR images with small baselines to get high degree of spatial coverage of distributed scatters [6]. Both algorithms can compensate the disadvantages of the conventional D-InSAR, namely, phase errors due to geometrical and temporal decorrelation as well as atmospheric disturbance [31].
The ENVISAT-ASAR data were processed using both the PS and SBAS-InSAR methods. For PS-InSAR, we selected the slant-range image (28 November 2006) as a master image and generated 20 interferograms. The perpendicular baseline ranged from 43 m to 951 m, the temporal baseline ranged from 68 days to 1277 days, and selected the potential PS candidates with a coherence threshold of 0.67. For SBAS-InSAR, we selected the slant-range image (28 October 2008) as a super master image. To reduce the geometrical and temporal decorrelation, the threshold criteria with a maximum temporal baseline of 735 days and a maximum perpendicular baseline of 483 m were used, and then, 68 interferograms were generated. Figure 2 shows the time-baseline plot for the ENVISAT-ASAR images used for PS and SBAS-InSAR processing. A combination of the minimum cost flow (MCF) network [41] and Delauney 3D method [42] was employed for phase unwrapping with an unwrapping coherence threshold of 0.35. After that, a screening process of flattened, filtered interferograms and unwrapped phases for checking some unwanted behaviors and data problems, which were caused by strong orbit inaccuracy, non-coherent pairs, atmospheric artefacts, residual topography etc., and 12 interferometric pairs were discarded from further processing. For refinement and reflattening, the reference points where the unwrapped phase value close to zero and flat areas identified from a topographic map projected in the line of sight (LOS) direction were chosen [43]. A custom atmospheric filtering was performed with a low-pass spatial filter with a 1.2 km × 1.2 km window on each single acquisition, and a high-pass filter at 365 days on these preliminary results to recover the final and cleaned displacement time series. Finally, geocoding was done in the original satellite LOS direction with a pixel resolution of 25 m. The same software package used for D-InSAR was also used in the PS and SBAS-InSAR processing chain.

Pearson Correlation Coefficient and Root Mean Square Error
The deformation rate derived from different methods with respect to each pixel in major deformation areas were compared using the Pearson correlation coefficient (r) and the root mean square error (RMSE). During the estimation, we only considered the common valid pixels taken from the deformation map of three InSAR methods. The Pearson's r is often measured as a correlation coefficient with 0 indicating no linear relationship between the variables, +1 indicating a perfect increasing linear relationship and vice versa [44].
where n is the total number of common valid pixels of the selected reference points, x and y (x and y mean values) are the values of the pixels taken from the selected reference points.
The RMSE values can be used to compare the individual model performance to that of other predictive models; it is often used in the field of geosciences, where researchers use the RMSE for model errors and comparison [45][46][47].
where n is the total number of common valid pixels of the selected reference points, X m,i and X s,i are the deformation rate values of the reference points taken from the deformation map of two InSAR methods, respectively.

Results of D-InSAR Interferometry
The coherence is an indicator of InSAR data quality and significantly influences the accuracy of the phase in interferograms [48]. As shown in Figure 3, the mean coherence is lower (0.39) for the data pairs with long time intervals and relatively large perpendicular baselines (20 January 2007 and 10 December 2008, B p = 643 m). Meanwhile, the mean coherence is higher (0.65) for the data pairs with short time intervals and relatively small perpendicular baselines (10 December 2008 and 25 January 2009, B p = 472 m). Good quality interferograms were generated from these two data pairs ( Figure 4); the interferogram with a relatively small perpendicular baseline shows good coherence (Figure 3b), but no significant land deformation was found (Figure 4b). It is due to the short interval of time, which was 46 days. Therefore, we used the differential interferogram generated from the large perpendicular baseline pair, which showed clear deformation at two sites in the study area. These areas are indicated by blue dashed circles in Figure 4a. Although the small time baseline image pair was not used for final displacement map generation, it served as a good example of showing the difference of different perpendicular baseline image pairs, and, provided useful clues for further data selection and analyzes.  In general, the sequence of color fringes in the interferogram can be used to determine whether uplift or subsidence has taken place. We chose two close-up views of typical fringe maps (Figures 4c,d) from these two deformation sites for further analysis. The ALOS-PALSAR data were in the ascending orbit in which observations were made from the west. The sequence of color fringes in these two figure shows a decrease in range (yellow-purple-blue), which denotes that the ground objective moves closer to the satellite. In this study, we did not consider the horizontal movement, and thus, it corresponded to an uplift. During the process of extracting the ground displacement, we unwrapped the interferogram in order to solve the 2π ambiguity, and corrected the satellite orbit   In general, the sequence of color fringes in the interferogram can be used to determine whether uplift or subsidence has taken place. We chose two close-up views of typical fringe maps (Figures 4c,d) from these two deformation sites for further analysis. The ALOS-PALSAR data were in the ascending orbit in which observations were made from the west. The sequence of color fringes in these two figure shows a decrease in range (yellow-purple-blue), which denotes that the ground objective moves closer to the satellite. In this study, we did not consider the horizontal movement, and thus, it corresponded to an uplift. During the process of extracting the ground displacement, we unwrapped the interferogram in order to solve the 2π ambiguity, and corrected the satellite orbit In general, the sequence of color fringes in the interferogram can be used to determine whether uplift or subsidence has taken place. We chose two close-up views of typical fringe maps (Figure 4c,d) from these two deformation sites for further analysis. The ALOS-PALSAR data were in the ascending orbit in which observations were made from the west. The sequence of color fringes in these two figure shows a decrease in range (yellow-purple-blue), which denotes that the ground objective moves closer to the satellite. In this study, we did not consider the horizontal movement, and thus, it corresponded to an uplift. During the process of extracting the ground displacement, we unwrapped the interferogram in order to solve the 2π ambiguity, and corrected the satellite orbit inaccuracy and phase offset using the collected external ground control points (GCPs). The final geocoded ground displacement map in the LOS direction is shown in Figure 5. The figure shows that two sites in the north-eastern part of the city exhibit clear indications of land uplift, and the maximum LOS displacement was approximately 13 cm in the study period.
Appl. Sci. 2017, 7, 72 7 of 14 inaccuracy and phase offset using the collected external ground control points (GCPs). The final geocoded ground displacement map in the LOS direction is shown in Figure 5. The figure shows that two sites in the north-eastern part of the city exhibit clear indications of land uplift, and the maximum LOS displacement was approximately 13 cm in the study period.

Interferometric Phase Change over Time
The interferometric phase is useful for analyzing the ground deformation trends over time in the LOS direction [49]. The time series interferograms generated from each data pair from (a) to (k) are shown in Figure 6. Generally, we could roughly identify the deformation areas in the large scale by analyzing the fringe density in differential interferograms. In addition, the sequence of color fringes in the interferogram can be a good indicator of whether uplift or subsidence has taken place. Figure 6a shows a typical interferometric phase image, in which the fringes clearly show the areas with a high deformation rate. The ENVISAT-ASAR data were in the descending orbit and observations were made from the east. As shown in Figure 6b, the sequence of color fringes shows a decrease in range (yellow-purple-blue), which denotes that the ground objective moves closer to the satellite. In this study, we did not consider the horizontal movement, and thus, it  (Figure 6g-k), and a pattern of rebound, expanding in area and increasing in quantity, can be observed during this period.

Interferometric Phase Change over Time
The interferometric phase is useful for analyzing the ground deformation trends over time in the LOS direction [49]. The time series interferograms generated from each data pair from (a) to (k) are shown in Figure 6. Generally, we could roughly identify the deformation areas in the large scale by analyzing the fringe density in differential interferograms. In addition, the sequence of color fringes in the interferogram can be a good indicator of whether uplift or subsidence has taken place. Figure 6a shows a typical interferometric phase image, in which the fringes clearly show the areas with a high deformation rate. The ENVISAT-ASAR data were in the descending orbit and observations were made from the east. As shown in Figure 6b, the sequence of color fringes shows a decrease in range (yellow-purple-blue), which denotes that the ground objective moves closer to the satellite. In this study, we did not consider the horizontal movement, and thus, it corresponded to an uplift.

Mean Deformation Velocities and Temporal Evolutions
Analyses of the PS and SBAS-InSAR results from the ENVISAT-ASAR dataset show that significant surface uplift occurred over or near the oilfields over the study period (Figure 7). The spatial distributions and mean displacement velocity for the period from 30

Mean Deformation Velocities and Temporal Evolutions
Analyses of the PS and SBAS-InSAR results from the ENVISAT-ASAR dataset show that significant surface uplift occurred over or near the oilfields over the study period (Figure 7). The spatial distributions and mean displacement velocity for the period from 30 September 2003 to 15 June 2010 are shown in Figure 7a (for PS-InSAR), Figure 7b (for SBAS-InSAR). Temporal decorrelation of bare land and agricultural areas were excluded by the coherence threshold. Major uplifted areas were highlighted by both PS and SBAS-InSAR methods and indicated in red dashed circles in Figure 7a,b. Dense PS points were detected around the reservoir area, whereas the distributions of PS points were sparse in other areas. The maximum deformation velocity was approximately 24 mm/year. The locations and the mean deformation velocity detected by the SBAS-InSAR technique were much clearer than the PS-InSAR results (Figure 7b). SBAS-InSAR also revealed a major uplift area that was not detected using the PS-InSAR, as indicated by the black dashed circles in Figure 7a InSAR technique were much clearer than the PS-InSAR results (Figure 7b). SBAS-InSAR also revealed a major uplift area that was not detected using the PS-InSAR, as indicated by the black dashed circles in Figure 7a,b. The estimated maximum deformation velocity was approximately 33 mm/year in the LOS direction.   (Figure 7). Three typical reference points were selected with the induced deformation as T1-T3 and one stable reference point as T4, which are shown in Figure 7a,b. Figure 8 illustrates that both PS-InSAR and SBAS-InSAR methods have good agreement in the deformation trend, with the time evolutions of each selected pixel sharing similar characteristics, such as increasing uplift. Despite a notable

Comparison and Discussion
Quantitative accuracy assessment is essential for evaluating the usability and reliability of InSAR-derived deformation [49]. However, except for the one GPS station, no other ground measurements are available, and thus, an alternative approach was used to perform cross-validation among the results from the three InSAR methods.
To compare the deformation results obtained from the D-InSAR, PS and SBAS-InSAR methods, the mean deformation velocity of D-InSAR processing results was calculated by averaging the displacement values. As explained in previous Section 4.

Comparison and Discussion
Quantitative accuracy assessment is essential for evaluating the usability and reliability of InSAR-derived deformation [49]. However, except for the one GPS station, no other ground measurements are available, and thus, an alternative approach was used to perform cross-validation among the results from the three InSAR methods.
To compare the deformation results obtained from the D-InSAR, PS and SBAS-InSAR methods, the mean deformation velocity of D-InSAR processing results was calculated by averaging the displacement values. As explained in previous Section 4. We selected 500 reference points from the PS-InSAR results in the study area, including both the stable and deformed areas (Figure 9d). The LOS displacement velocity for these points was plotted compared with the SBAS and D-InSAR results, and the correlation coefficients were calculated. The results show good agreement between the SBAS and D-InSAR methods with a correlation of 0.75, and the RMSE was approximately 4.918 mm/year. However, the results from the D-InSAR and PS-InSAR were not as comparable, with a low correlation of 0.54 and an RMSE of 6.043 mm/year. However, all the three methods could identify the main deformation areas, except for one of the main deformation areas in PS-InSAR indicated in black circles in Figure 7a. This may be caused by low backscatter in the wide non-urban land use, where the sufficient number of PS points could not be detected in the PS-InSAR processing. Although the XJKL is the only GPS station located within the study area and is insufficient to measure and represent the entire deformation area, it provides a means of assessing the accuracy of the InSAR data and methods used in this study. The deformation rate at the XJKL continuous GPS station was 9.6 mm/year, while the D-InSAR and SBAS-InSAR derived deformation rates were 1.2 and 7.4 mm/year, respectively. No displacement was found in the PS-InSAR measurement results. In addition, the identified deformation area and the temporal evolution pattern in this study agree well with the results obtained by Ji et al. [36]. In our study, the maximum deformations from 2003 to 2010 for SBAS and PS-InSAR were approximately 22.7 cm and 17.8 cm, respectively, and the results of D-InSAR was 13 cm for the period 2007 to 2009, which are comparable to the previous research results in which the maximum uplift magnitude was reported as 20 cm over the years 2006 to 2010 [36]. The comparison result shows slight differences between the two results, which might be caused by the differences in the deformation trend in the different periods, rather than any deficiencies in one method over the other. As shown in Figure 6, the deformation was evident and more accelerated from 2006 to 2010 than during the previous time of 2003 to 2006.

Conclusions
To assess the degree of deformation in the oilfields of Karamay, Xinjiang, China, the D-InSAR technique was applied to ALOS-PALSAR data and the PS-InSAR and SBAS-InSAR to the ENVISAT-ASAR dataset. The results showed that all the three methods could provide useful information for Although the XJKL is the only GPS station located within the study area and is insufficient to measure and represent the entire deformation area, it provides a means of assessing the accuracy of the InSAR data and methods used in this study. The deformation rate at the XJKL continuous GPS station was 9.6 mm/year, while the D-InSAR and SBAS-InSAR derived deformation rates were 1.2 and 7.4 mm/year, respectively. No displacement was found in the PS-InSAR measurement results. In addition, the identified deformation area and the temporal evolution pattern in this study agree well with the results obtained by Ji et al. [36]. In our study, the maximum deformations from 2003 to 2010 for SBAS and PS-InSAR were approximately 22.7 cm and 17.8 cm, respectively, and the results of D-InSAR was 13 cm for the period 2007 to 2009, which are comparable to the previous research results in which the maximum uplift magnitude was reported as 20 cm over the years 2006 to 2010 [36]. The comparison result shows slight differences between the two results, which might be caused by the differences in the deformation trend in the different periods, rather than any deficiencies in one method over the other. As shown in Figure 6, the deformation was evident and more accelerated from 2006 to 2010 than during the previous time of 2003 to 2006.

Conclusions
To assess the degree of deformation in the oilfields of Karamay, Xinjiang, China, the D-InSAR technique was applied to ALOS-PALSAR data and the PS-InSAR and SBAS-InSAR to the ENVISAT-ASAR dataset. The results showed that all the three methods could provide useful information for identifying the boundaries of deformation and revealed two areas of land uplift near the oilfield. The maximum deformation velocity was estimated to be 33 mm/year. As previous research suggests, the subsurface water injection used to enhance oil recovery is the likely cause of this land uplift. The comparison of results from each of the three methods indicates that the correlation of the SBAS-InSAR and D-InSAR results is higher than that of the PS-InSAR and D-InSAR results. For this study area, the SBAS-InSAR method was found to be more robust than the PS-InSAR method. Future investigation can improve the findings of this study by utilizing GPS base stations in the high deformation areas and monitoring long-time land deformation in combination with more accurate and improved satellite data.