Near-field motion of the Haiyuan fault zone in the northeastern margin of the Tibetan plateau derived from InSAR permanent scatterers analysis

Abstract The Haiyuan fault zone is a major discontinuity in the northeastern margin of the Tibetan plateau. A magnitude 8.5 earthquake occurred there in 1920. Both geological investigations and GPS measurements show that this fault zone is still highly active, with a slip rate of 3 to 10     mm / year , exhibiting a large range of variance in both space and time. We attempt to use the permanent scatterers interferometric (PSI) synthetic aperture radar technique to better detect the near-field motion of this fault zone. We process and analyze 38 scenes of ENVISAT/ASAR images from two neighboring descending orbits using the PSI method. The results show a remarkable velocity gradient of about 5     mm / year across the central segment of the fault zone and a rate of about 5 to 6     mm / year on its eastern segment. The motion senses are consistent with a left-lateral strike slip. The motion histories of most PS points show a stable linear variation trend in time series. In addition to these motion features that agree with those from geological, GPS and other observations, the dense PS analysis also reveals spatially continuous variations of crustal motion around the fault zone.


Introduction
The Haiyuan fault zone, lying in the northeastern margin of the Tibetan plateau, is a major zone of tectonic discontinuity between active blocks and an important seismic belt in western China. It begins south of Xiaokou, Guyuan county, Ningxia Autonomous Region in southeast, and extends northwestward to Xingquanpu, Jingtai county, and Gansu province, belonging to the eastern section of the North Qilian fault system. This zone comprises three en echelon left-slip faults, overall trending in the N280 to 300°W area, with a total length ∼240 km (Fig. 1). Left slip has continued to characterize its motion since the late Quaternary (describes a period in geologic history) with a cumulative displacement up to 10 to 15 km. 1 The Haiyuan M8.5 earthquake of 1920 took place on the eastern segment of this fault zone, producing a left-lateral strike slip of up to 10.5 m, as well as a 237-km long surface rupture belt. [1][2][3] The analysis of long term and short term active behaviors of the Haiyuan fault zone is very important not only for estimating and predicting seismic hazards, but also for understanding the crustal motion in the northeastern margin of the Tibetan plateau. There have been a lot of studies on the slip rate of the Haiyuan fault zone. Using geological methods and investigations of the features of the terrain around the fault, and using geological dating, the Holocene slip rates of various segments of the Haiyuan fault zone are estimated to be 3 to 10 mm∕year. [4][5][6][7] Most of the recent GPS measurements suggest that this fault zone is in a left-slip state at rates of 3 to 8 mm∕year, [8][9][10][11] though individual work reports slip rates of only 1.1 to 1.7 mm∕year. 12 The discrepancy in these geodetic estimations may be associated with differences in measurement sites, station spacing, and analysis methods.
In recent years, InSAR technology has been applied to the measurement of crustal motion around faults and especially to the observation of co-seismic deformation fields. Although in theory, this method can detect displacement of mm orders of magnitude on the surface, interseismic crustal motion or fault slip is so slow and tiny that its signal can be blanketed by errors from decoherence factors such as satellite orbits atmospheric delay, and terrain conditions. To solve this problem, some improved InSAR techniques have been proposed, such as permanent or persistent scatterers InSAR (PSI), small baseline subsets (SBAS), and interferograms stacking. All these methods have been applied to detection of surface deformation.  Of them, the PSI uses highly coherent point targets, called "PS," on the time series synthetic aperture radar (SAR) image with a high temporal stability of the backscattered signal, such as houses, bridges, big stones, and mountain crests, as objects for phase analysis, so that noncoherent image elements with low signal-to-noise ratio are avoided. 13,14 Once a PS data set is determined based on a criterion and algorithm, the analysis of the the phase time series for each PS allows one to detect near-field crustal motion of a fault.
In this work, we attempt to use the PSI method to conduct an experimental study of crustal deformation measurements on the Haiyuan fault zone. Our purpose is to examine whether this approach can reduce decorrelation of phase in repeat SAR images caused by some natural conditions in western China, and help better detect the near-field motion of active faults like the Haiyuan fault zone. A dry climate and rare vegetation characterize this region, most of which is intact mountains. The exposed ridges and slopes are natural point targets for PS selection, so this area is suitable for measurements of fault activity using PSI. At present, due to limitations of the technology, the PSI can be used to detect crustal motion of small-sized areas, so we choose two small experimental areas for this work. One crosses the central section (104.75°E, 36.78°N) of the Haiyuan fault zone, as shown in Fig. 1(b). On the south and north sides of this area are mountainous lands with steep relief and little vegetation. Its central portion, where the Haiyuan fault zone runs, is relatively flat with towns, villages, and croplands. The Yellow River and flood land are present in the southwest corner. The northern side of the eastern section (105.60°E, 36.50°N) with relatively flat terrain, is where Haiyuan county is located as shown in Fig. 1(c). Here, on the south side of the Haiyuan fault zone are steep mountains, while on the north side is flat land with a town. Our tests show that the experimental area covering both sides of the fault zone with large terrain differences can yield bad results, so we have to select an area on the flat northern side.
The SAR data used in this work are ENVISAT/ASAR from the European Space Agency, with image mode, IS2 incidence mode, incidence 19 to 26 deg with a 23 deg average, C band and VV polarization, in a descending orbit from track 2290 to 2018. We have 17 scenes from 20040119 to 20100308 in track 2290, and 21 scenes from 20030813 to 20091209 in track 2018. Considering the spatial perpendicular baseline, temporal baseline, Doppler center frequency difference, and atmospheric errors, we selected one master image and form interferometric pairs (including one automatic pair). The ASAR data and interferometric parameters are shown in Tables 1 and 2 and Fig. 2, respectively. For most pairs, the perpendicular baseline is <400 m, which is helpful for maintaining good coherence in the mountainous area. Here, it should be mentioned that although there should be no spatial decorrelation occurring on the point target due to its dominated echo in the located pixel, in fact its phase is still affected by geometric decorrelation to some extent because of the limited number real point targets.

Data Processing
We use the interferometric point target analysis (IPTA) module of the software Gamma to process the ASAR data. The whole flowchart is shown in Fig. 3. PS candidates (PSCs) are determined based on the co-registered single look complex (SLC) images. For each PSC, the SLC values are extracted and written to a point data stack. Next, the differential interferograms of the point data stack are calculated. From the differential phases of good PSCs, deformation, height correction, and atmospheric delay components are estimated using an iterative procedure. That is to say, the interferograms are only generated and interpreted for the selected points in IPTA. For efficiency and data storage reasons, vector format data structures are used instead of the raster data format used in conventional interferometry. 43 Before interferometric processing, we make a radiometric calibration of each SAR image used in this study acquired at different times, so that their intensity values are comparable, and good to the followed PSCs selection [convenient for selection of permanent scatterer candidates (PSCs) from many SAR images]. We also make perpendicular baseline corrections to the interferograms based on the accurate orbit data from Delft University. The initial topographic phase component is removed using the 3-arc-sec-SRTM DEM. Considering the mountainous terrain of the study area, we adopted a new co-registration algorithm based on a lookup table, which can greatly improve co-registration accuracy and increase coherence of the image pairs in the mountain area. As shown in Tables 1 and 2, most image pairs have a coregistration accuracy <0.2 pixel, which is the precondition for success of the chosen processing methods. The detailed processing procedures are described below.

Generation of PSCs
In processing the time-series InSAR analysis, one usually adopts amplitude stability and coherence stability or either of them (i.e., correlation) as pixel selection criteria. The choice of the selection criterion depends on the individual application and the point target processing method to be used. In this study, we use Gamma software programs to select the PSCs automatically. 43 Both amplitude stability (low temporal variability) and spectral phase diversity are utilized as selection criteria of PSCs. We first identify point target candidates (PSC1) based on the lower temporal variability of the backscattering (ratio of mean to standard deviation) in the average amplitude image of all co-registered SLCs, which is better for large SLC data stacks. Then we   detect point target candidates (PSC2) based on low spectral phase diversity in a single co-registered SLC, which is better for small SLC data stacks. Finally, the two point target candidates are merged into a single one in which only those with common positions (column/pixel) are kept (Fig. 3). In the followed regression analysis, the quality of these PSCs is evaluated and some noisy PSCs are dropped.
In the experimental area of the central fault segment, we select 69143 PSCs from point targets in bare mountainous terrain. These PSCs are distributed nonuniformly, being relatively dense in the south and north mountainous sections, but sparse in the middle and southwest flat parts. Field investigations show sandstone, mudstone, and shale exposures in the northern and southern sections of the study area with steep terrains and linear ridges result in strong and stable backscattering, providing a lot of point targets. In the middle of the experimental area, some villages and towns are present in a flat terrain covered with farm fields and vegetation, thus few PS points can be identified. In the other experimental area, located at the eastern fault segment, we select 1204 PSCs that are concentrated around Haiyuan county. This implies that the distribution and density of PS points depend on ground surface conditions. As the quality of these PSCs would be evaluated and those of poor quality would be eliminated in the subsequent regressive analysis, the threshold for PSC selection set in this stage is not very high.

Geocoding of SAR data
For InSAR data processing, orthorectification means geocoding of SAR data. This is the coordinate transformation between the coordinates of range-Doppler coordinates of the SAR and orthonormal map coordinates. This type of geocoding takes the image in the radar geometry (range azimuth) and re-samples the data into a map projection, such as UTM with an associated datum. Depending on the ground surface considered for geocoding, it can be distinguished between ellipsoid corrected and terrain corrected geocoding. The former does not require a digital elevation model (DEM) and solely considers a constant height, while the latter needs the DEM and takes into account the surface topography which definitely has an effect on the look vector. 43 Due to the slope terrain of this study area, we use the terrain corrected geocoding method. The geocoding procedure includes three steps: First, it builds an initial geocoding lookup table between the SAR and map geometry based on the orbital data and the parameters of the master image to be geocoded as well as that of the DEM. Second, it refines the lookup table through simulating the SAR intensity image based on the DEM and the SAR imaging geometry, and calculates the registration offsets between the simulated SAR image and the true SAR image. Third, it implements the transform from radar geometry into map geometry based on the refined geocoding lookup table, with the advantage of just a single re-sampling and interpolation (Fig. 3).

Calculation of initial differential phases
We extracted SLC data of PSCs and generated initial interferograms by interferometric processing focused on only the point data stack, as well as extracted elevation values of PSCs from 3-arc-sec SRTM DEM and calculated the simulated topographic phases (initial phase model) of the PSCs. Then by subtracting the initial phase model from the point interferograms, the initial differential phases are determined (Fig. 4). Figure 4 shows 17 differential interferograms of PSCs with different perpendicular baselines in the central segment of the Haiyuan fault zone. It is clear that for the interferometric pairs with a short perpendicular baseline <300 m, the initial differential interferograms are relatively smooth, implying good coherence. For those pairs with relatively long perpendicular baselines >300 m, however, they are rough in variation and disordered in color, indicative of high noise level. It suggests that the perpendicular baseline is still a very important influence factor for point interferometric processing, so we adopt an appropriate strategy of regression analysis according to perpendicular baselines of all point interferograms.

Regression analysis of point differential phases
The initial differential phases of PSCs may contain various phase error components, such as orbital error, topographic residual, atmospheric delay, and other noise. These error phases can be estimated and removed through regression analyses iteratively based on a two-dimensional linear phase model, which indicates the topographic phase is linearly dependent on the perpendicular baseline component of the interferometric pairs, and the crustal motion rate is linearly dependent on the time interval between the image pairs. Thus the related slopes of regressions correspond to relative terrain height corrections and relative linear deformation rates. During regression analyses, a Delaunay triangle grid is constructed and phase unwrapping of differential phases of PSCs is completed by the minimum cost flow (MCF) algorithm. The phase standard deviation and residual phases derived from the regression are used as a measure for the quality of phase unwrapping and regressive fitting. Low phase standard deviation and smooth residual phases are indicative of correct unwrapping, whereas large phase standard deviation and uneven residual phases exhibiting block-looking variation mean wrong unwrapping and iterative calculation needs to continue. 43 As the original differential interferograms have distinct coherences, we adopt the following procedures for regression analyses. First, regression processing and phase unwrapping are done on those interferometric pairs (layers) with a perpendicular baseline <300 m in order to search the layers of high quality. Using these good layers, the DEM error model and linear deformation rate model for some PS points are reliably established. Second, more layers and points are unwrapped based on these models, while those points of poor quality are removed. Third, perpendicular baseline correction is performed by the least square algorithm based on stable points with motion rates <2 mm∕year. The regression processing is conducted again using the optimized baseline, and the DEM error model and linear deformation rate model are renewed. Finally, the atmospheric delay phase is estimated and eliminated by spatial filtering to residual phases and a new linear deformation model after atmospheric correction as well as the uncertainty of parameter estimates and phase standard deviation is generated again. These resultant images are shown in Figs. 5-12. It should be mentioned that the reliability of estimation of the linear motion rate is related to the quality of PS points, while the quality of PS points is associated with their density. A trade-off is needed between the quality and the density of PS points.  agrees with the feature of fault motion. On the other hand, some colored spots with low quality can be found in the condition of more PS points, such as Figs. 5(a), especially in the northern mountain land of the study area. If the number of PS points is reduced by setting a lower threshold of phase standard deviation in the regression analysis, these low-quality points could decrease considerably or even disappear [Figs. 5(b)-5(d)]. It implies that both the quality of PS points and reliability of motion rate estimates are enhanced with a decrease of PS points. However, the observation from sparse and few PS points would lose the advantage of revealing the spatial variation trend of the motion rate around the fault zone and motion differences on either sides. Thus, a trade-off between the number of PS points and the quality of deformation rate estimates is required, as in Fig. 5(b), which should be the optimal result in a relative sense. Figure 6 displays phase standard deviations corresponding to PS point amounts 67665, 59584, 39993, and 12333, respectively, which measure the quality of regression analyses and reliability of motion rate estimates. It shows that with a decrease of PS point number, the standard deviations decrease. For example, in Figs. 6(b)-6(d), these deviations are as small as ≤1.0 rad in general, while in Fig. 6(a) with relatively more PS points, the deviation increases in the north and south parts and ≥1.0 rad, especially in the north mountains, which suggests that there are some points of low quality there. The motion rates of PS points are also more heterogeneous and dispersed, indicative of poor reliability [ Fig. 5(a)]. The reason for this may be radar shadows and overlapping in the mountain areas which lead to a decline of point quality. Therefore, in the authors' opinion, it would be better to set different and proper thresholds of criteria for PS point selection according to specific conditions of topography and vegetation coverage, and perform regression analyses separately in distinct regions.     In order to describe the motion history of the study area, we select four PS points on the north side of the central Haiyuan fault zone and plot their cumulative displacement-time curves (Fig. 9). They exhibit linear changes with time, and a largely consistent cumulative displacement of about 4 cm on average with a negative sense, i.e., increase of the distance between the satellite and the Earth. Assuming the ground moves along with the fault zone only horizontally, these changes mean westward motion of PS points on the north side in the mode of the right looking and descending orbit, which is consistent with the left slip of the Haiyuan fault zone.

Results and Analysis
Cumulative displacement images of all 17 interferometric pairs are shown in order of time (Fig. 10). Each of them reveals displacements of PS points over the interval between two acquisition times of one image and the reference image. The colors of these serial displacement images change slowly and stably with time, implying a linear variation of PS point motion. They should reflect a long-term slow interseismic motion of the crust rather than local changes caused by the atmosphere or other factors. It is noted that four PS point interferograms before 2008 (20040119, 20040607, 20041025, and 20070806) have a relatively conspicuous color contrast with respect to those after 2008, perhaps indicating a change of crustal motion around 2008. The other 12 interferograms after 2008 (20080303 to 20100201) have a small color contrast, probably indicating a stable and slow motion process. It is likely that such a phenomenon is associated with the Wenchuan great earthquake of 2008, or implies an accelerated motion on the Haiyuan fault zone even before that event. It is reported that an anomalous deformation was observed at several sites on the Qilian Shan and Haiyuan faults >600 km to the north of Wenchuan, 44,45 consistent with the above results. Unfortunately, there is little ENVISAT/ASAR data from 2006 and 2007 available for the study area, limiting further exploration of the crustal motion anomalies prior to the 2008 Wenchuan great earthquake.  The resultant motion rate maps of PS points on the eastern segment of the Haiyuan fault zone are shown in Figs. 11 and 12, respectively, where the PS points are mainly concentrated in Haiyuan county in the center of the experimental area. In this test area, the fault motion rates revealed by most PS points are about 5 to 6 mm∕year [ Fig. 11(a)], and the corresponding phase standard deviations are mostly <0.7 rad [ Fig. 11(b)]. A time series of cumulative displacement images of PS points (Fig. 12), similar to those in the central segment of the study area, display a gradually stable variation trend, and a total displacement up to 4 cm from August 2003 to December 2009. However, the cumulative displacement images before 2008 show more visible color differences than those after 2008, which is a very interesting phenomenon. We shall study it based on more data in the future.

Discussion
The radar satellite measures range variation in the LOS between two flights over the same target.
In theory, such measurements can be converted into strike slip and dip slip of a fault by a geometric relationship (Fig. 13). For a fault like the Haiyuan fault zone, we assume that Dstrike is a strike slip of the fault (southeastward is positive), Ddip is the dip slip (upward is positive), Dlos is displacement in LOS (toward satellite is positive), θ is the incidence angle of the radar, α is the orbit azimuth of the satellite (clockwise is positive), and β is the strike of the fault (including the angle to the north). Then, we have the following relationship: In this work, we use the images ENVISAT/ASAR of a descending pass with an average incidence of 23 deg and a satellite orbit azimuth 192 deg. The Haiyuan fault zone tends in 280 to 300 deg, so we take 290 deg. If we assume that this fault is of a pure strike-slip motion, the calculation by the relationship above shows that a motion rate of 5 mm∕year in LOS corresponds to 12.9 mm∕year along the fault. If we assume that the fault also has a 1 mm∕year of dip slip in addition to its strike-slip motion, a strike-slip rate of 10.5 mm∕year can be derived from the 5 mm∕year in LOS. GPS measurements (3 to 8 mm∕year) are the average values on a large scale. In this work, the PSI analysis attempts to reveal motions in a small area, likely with a localized feature. Besides, the results from this approach might include the dip slip component of the fault, while GPS is only a horizontal motion. Thus, we think the relatively high value in this work (10.51 mm∕year with dip slip or 12.9 mm∕year pure strike slip, see Table 3) is acceptable, though it does have some error. Of course, this issue should be studied further. The PSI analysis reveals a gradient contrast of motion rates across the central segment of the Haiyuan fault zone. This is not what we expected if the upper crust portion is presumably locked above the seismogenic layer. Nevertheless, such a phenomenon has been observed on many active faults in the world. 20,23 Using the interferogram stacking method that converts average phase change into velocity, Cavalié et al. 22 found a zone of high velocity gradient across the western Haiyuan fault zone which is a few kilometers wide. They suggest a screw dislocation model in an elastic half-space with locking at depth or an alternative model with shallow creep and deep aseismic slip to explain the InSAR observations. The current PSI cannot resolve the discrimination between these two models by itself. We expect to combine PSI, GPS, and seismic observations to yield a better interpretation on the motion behavior of the Haiyuan fault in the future.
This work shows that the PS point cumulative displacements of all interferograms in the two experimental areas exhibit conspicuous changes before and after 2008, likely indicating the effect of the 2008 Wenchuan great earthquake on the Haiyuan fault zone. Because only a few SAR scenes before the shock are available, with data from 2006 and 2007 absent and only three scenes available for 2005 in two tracks, our effort to study the variation of motion rates on the Haiyuan fault zone before and after the Wenchuan event is limited. We hope to acquire more relevant data to fill this gap in our following work.
This work demonstrates that the reliability of the processing result of PSI is associated with the number and quality of PS points that are related to surface coverage of the study area. PSI can work well in mountainous regions with moderate vegetation where PS points are usually dense and stable. But in such settings, the distribution of PS points are sometimes not uniform due to the concentration of partial PS points on ridges and slopes that are facing the incidence of radar wave fronts. In the regions with heavy vegetation, PS points are sparse and unstable, thus the usage of PSI is limited. A possible approach to solve this problem is to combine PS and CR (corner reflector) techniques, i.e., to increase observational point targets by the deployment of corner reflectors in the places lacking PS points.
This work was conducted on a small area of the central and eastern Haiyuan fault zone, thus cannot describe the variation of near-field motion rates along the entire fault zone as well as a wider range across the fault zone. Meanwhile, the SAR data used a span only about 6 year. With the increasing spatial scale of a study area, the phase errors of PS points from baselines, topography, and atmospheric delay would become larger, probably leading to failure of phase unwrapping. This is an issue in PSI that remains to be solved. One possible approach is to perform processing on many small areas separately and then stitch them together so that these small areas are continuous and comparable in motion rates. It will be a focused topic of our future work. On the other hand, we will use other time-series InSAR methods such as SBAS, and stacking to research larger regions.

Conclusions
In this work, we attempt to detect near-field crustal motion of the Haiyuan fault zone based on 38 ENVISAT/ASAR images of two neighbor descending tracks through a PSI analysis. The 3-arcsec SRTM DEM and IPTA module of the software Gamma are used in the data processing. The results show that the motion rates of PS points in the two experimental areas of the fault zone are all smooth and consistent in spatial distribution, exhibiting a feature of steady and linear variance in time. Meanwhile, a conspicuous contrast of motion rates is present on opposite sides of the Haiyuan fault zone, and a relative motion rate in LOS about of about 5 mm∕a is derived, which can be converted into the left-lateral strike slip of 13 mm∕year of the fault under the assumption of pure strike-slip faulting. This value is larger than 3 to 10 and 3 to 8 mm∕year determined by geological and GPS observations, respectively. A possible explanation for this discrepancy is that the estimate from geology is a long-term average (e.g., since Holocene) of crustal motion or fault slip, and estimates from GPS are averaged on a large scale, while the result from PSI represents a near-field motion picture of a fault with a localized character. Besides, the motion rates from PSI analysis likely contain dip-slip displacement of the fault in addition to the horizontal slip. Of course, this issue requires further studies with more data and modeling.
Guohong Zhang received his PhD degree from Institute of Geology, China Earthquake Administration in 2011. His major is in geophysics, his current research interests mainly include inversion of co-seismic fault slip, simulation of crustal deformation field and remote sensing. He has published more than 10 related papers in academic journals and conference proceedings.
Xiaogang Song Song received his PhD degrees in photogrammetry and remote sensing from Wuhan University in 2008, writing his thesis on atmospheric correction in InSAR interferogram using GPS and MODIS data. In June 2008. Now he works for Institute of Geology, Chinese Earthquake Administration. His research interests are in InSAR and earthquake monitoring. He has published more than 10 related papers in academic journals and conference proceedings.
Guifang Zhang received her PhD degree from Institute of Geology, China Earthquake Administration in 2012. Her major is in tectonic geology, her research focuses are mainly on InSAR and Earth's surface deformation associated with the earthquake cycle and other geophysical processes. She has published some related papers in academic journals and conference proceedings.