Interseismic deformation rate of the Haiyuan fault system based on the modified SBAS method

Yng Liu , Sen Hn , Luyun Xiong , Ying Peng ,d, Mengyo Go , Cijun Xu ,b,c

a School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China

b Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China

c Key Laboratory of Geophysical Geodesy, Ministry of Natural Resources, Wuhan University, Wuhan 430079, China

d Changguang Satellite Technology Co., LTD, Changchun 130102, China

Keywords:InSAR Modified SBAS method Haiyuan fault system Interseismic deformation rate

ABSTRACT Radar interferograms are usually influenced by factors such as atmospheric artifacts, orbital errors, and terrain errors.It is difficult to reduce the influence by using the conventional small baseline subset(SBAS)method when determining the deformation rate.This study uses the adjustment model with systematic parameters to improve the conventional SBAS method and employs it to determine the interseismic deformation rate of the Haiyuan fault system,providing a data reference for exploring the locking depth,strain accumulation state, and potential seismic risk assessment of different segments of the Haiyuan fault system. The results are as follows: (1) the simulation experiment verifies the feasibility and robustness of the modified SBAS method. This method can effectively reduce the influence of residual signals such as atmospheric artifacts, orbital errors and terrain errors in the interferograms. The deformation rate map can be significantly improved;(2)the deformation rate field in the radar's Line of Sight(LOS) direction shows that there are obvious differences between the north and south sides of Haiyuan fault system, which is consistent with the characteristics of the left-lateral strike-slip movement of the Haiyuan fault system.The deformation rate field and profiles reflect the complex trends among different segments of Haiyuan fault system in detail. (3) the deformation rate of the Jingtai pull-apart basin is higher than that of the surrounding areas, possibly indicating strong regional activity, which provides a reference for studying the seismic risk of the Jingtai pull-apart basin; and (4) the interseismic deformation rate and profiles across the fault show that the middle section of the Lao Hu Shan(LHS)segment and the western and middle sections of the Haiyuan segment are locked.

1. Introduction

In recent years, many strong earthquakes have occurred in China, like the 1976 MW7.5 Tangshan earthquake, the 2008 MW7.9 Wenchuan earthquake,the 2010 MW6.9 Yushu earthquake,and the 2013 MW6.6 Lushan earthquake, causing numerous casualties and significant economic losses [1,2]. The formation and occurrence of earthquakes are closely related to crustal deformation [1-3]. To a certain extent, the possibility of earthquake occurrence can be determined by analyzing crustal deformation.Therefore,long-term crustal deformation monitoring in fault systems is an essential part of earthquake prediction research [1]. The Haiyuan fault system located on the north-eastern margin of the Qinghai-Tibet Plateau is a large left-lateral strike-slip fault,composed of the Leng Long Ling segment(LLL),Jin Qiang He segment(JQH),Mao Mao Shan segment(MMS), Lao Hu Shan segment (LHS), and Haiyuan segment from west to east [2-5]. Historically, strong earthquakes have occurred in the generalized Haiyuan fault system and surrounding areas.The MS8.5 Haiyuan earthquake occurred in the eastern part of the fault system in 1920,and a 220-km-long rupture occurred along the fault[3,4].The MS8.0 Gulang earthquake occurred in the western part of the fault system in 1927(Fig.1)[1-4].In the 260-km-long Tianzhu seismic gap between the two earthquakes, where no earthquakes above MS6.0 have been recorded,there exists the danger of a large earthquake [3,4]. According to the paleo earthquake research, the earthquake recurrence interval of the Tianzhu gap is approximately 1000 years,and the last major earthquake occurred in 1092,about 1000 years ago [4,6]. Therefore, it is crucial to investigate the current activity of the Haiyuan fault system.

Many studies have analyzed the Haiyuan fault system's activity and seismic risk by using geological and GPS geodetic observations.Zhang et al. [7], Burchfiel et al. [8], and Li et al. [9] used faulted landforms or dating techniques to obtain the slip rates of the Haiyuan fault system as 8 ± 2 mm/yr, 5-10 mm/yr, and 4.5 ± 1.0 mm/yr, respectively. Based on the block or dislocation models,Gan et al.[10],Thatcher[11],Duvall and Clark[12],and Cui et al. [13]derived the strike-slip rates of the Haiyuan fault system using GPS observation as 8.6 mm/yr,5-6 mm/yr,0-4.5 mm/yr,and 3.6-6.5 mm/yr.Although the slip rates of the Haiyuan fault system obtained by different studies differ,they are generally in the range of 2-10 mm/yr [14]. As these methods have some disadvantages,such as high monitoring costs and long monitoring cycle, the derived results usually have shortcomings, such as untimely monitoring and incomplete data. Since the interferometric synthetic aperture radar (InSAR) technology was first proposed by Graham in 1974 [15], this technology has been widely used in deformation monitoring fields, such as volcanoes, earthquakes,urban subsidences,and glacier movements,owing to its advantages of being full-time, large area, and high spatial resolution [16]. Differential InSAR (DInSAR) technology can overcome the limitations of traditional methods and quickly acquire large-scale crustal deformations.However,DInSAR technology is easily affected by orbit errors, atmospheric artifacts, spatio-temporal decoherence, and terrain errors; therefore, it has some limitations in measuring the crust's micro-deformation [2].

Fig.1. Tectonic setting of the study area.Black solid circles represent cities;black lines represent faults; blue line shows the Jingtai pull-apart basin boundary; red stars represent the epicenters of the 1920 MS8.5 Haiyuan and 1927 MS8.0 Gulang earthquakes,respectively;red dashed box shows the coverage of Sentinel-1 SAR data;white circles indicate MS ＜6 earthquakes between 1900 and 2019; red circles represent MS ≥6 earthquakes between 1900 and 2019 (Data source: National Earthquake Data Center, China Earthquake Networks Center).

In recent decades,a series of high-resolution SAR satellites have been successfully operated, and high-resolution images provide a good guarantee for obtaining high-precision deformation. Timeseries techniques, such as permanent scatterers InSAR (PS-InSAR),small baseline subset (SBAS), and stacking, have been gradually applied to the study of interseismic deformation in fault systems[1]. Stacking technology can reduce the influence of atmospheric artifacts.However,it does not consider the effects of terrain errors[10]. PS-InSAR technology can overcome the influence of factors such as decoherence, noise, atmospheric artifacts, and terrain errors to the greatest extent. However, the quality of observation largely depends on the density of PS points and the number of SAR images. In the radar interferograms, the deformation phase is usually mixed with the atmospheric phase and orbit phase[10-18].It is usually to use an external independent atmospheric model to correct atmospheric artifacts,and use orbital correction method to remove orbit errors. However, different methods cannot completely remove the errors, resulting in some residual errors.The residual errors are often coupled with small interseismic deformation signals, and it is difficult to separate them [20].

This study proposes a modified SBAS method-the conventional SBAS method combined with the adjustment model with systematic parameters. This method considers the residual errors in the radar interferogram as unknown parameters,and then reduces the influence of these errors on the extraction of deformation maps.First,this study uses a series of simulation experiments to verify the feasibility and robustness of the modified SBAS method. Then,based on the modified SBAS method, this study extracts the interseismic deformation of the Haiyuan fault system with Sentinel-1 SAR data from 2018 to 2020. The results can provide more data support for assessing seismic risk, predicting strong earthquake locations,and exploring the deformation mechanism of the northeastern margin of the Qinghai-Tibet Plateau.

2. Modified SBAS method

2.1. Conventional SBAS method

The SBAS method was first proposed by Berardino et al.[21].The basic principle of the SBAS method is to select image pairs with a short spatial baseline from many SAR images by appropriate combinations, and then use the least squares (LS) or singular value decomposition(SVD)method to obtain high-precision deformation[8,16,17,21-23].

Assuming that N+1 images covering the same area are acquired according to the time series t0,t1...tn,they are co-registered with any image as the master image; the time baseline threshold and spatial perpendicular baseline threshold (generally less than 300 m) are also set. Image pairs that meet the time and space baseline thresholds form one or several baseline sets, and the image pairs in each set are processed to obtain the corresponding radar interferograms. Then, M interferograms can be obtained[15-22] with the following range:

When there are two or more small baseline sets, the matrix rank is deficient, and the SVD method will be used to solve the system equation [18,19]. From the above formulae, it can be seen that the influences of the residual errors related to terrain errors, atmospheric artifacts, orbit errors, and noise on deformation estimation are not considered in the conventional SBAS method.

2.2. Modified SBAS method

The modified SBAS method is proposed to reduce the influence of residual errors on the deformation extraction. This method is based on the conventional SBAS method and the adjustment model with systematic parameters [24]. This modified SBAS method considers residual errors by introducing the parameters S1,S2...Sk.Thus, equation (7) can be expressed as:

3. Simulation experiment

In the simulation experiment,the fault parameters are based on the results of the Haiyuan fault system from previous studies. The locking depth dris 13 km, and the fault slip rate Srbeneath the locking depth is 6 mm/a[20].The surface deformation rate field is generated by the screw dislocation model[25,26],which is shown as follows:

According to the projection parameters (-0.6136, -0.1092,0.7783) of Sentinel-1A ascending orbit in three directions of eastwest, north-south, and vertical in the Haiyuan area, we converted the deformation rate from the fault-parallel direction to the LOS direction[27],as shown in Fig.2(a).We design 20 LOS deformation maps with a time span of 3-5 years. To account for the residual errors existing in the real interferograms, a series of simulation errors, such as noises, terrain errors, orbit errors and atmospheric artifacts, are added to each deformation map.

Each map is added with random Gaussian white noise with a standard deviation of 0.5 mm. The random noise represents decoherent noises, terrain errors, and errors caused by processing methods [26,28-30]. One example of random noise is shown in Fig. 2 (b). Each map is added with different orbital errors, and the values of the range and azimuth gradients are 0.001-0.002 mm/km and 0.003-0.005 mm/km, respectively[19,28,31]. One example of orbital error is shown in Fig.2(c).Among them,large errors exist in the N deformation maps,and correlation noises are added to them[19,32,33]. One example of correlation noise with a mean value of 0 mm and a standard deviation of 2 mm is shown in Fig. 2 (d).Fig. 2(e) and (f) are examples of synthetic deformation map with correlated noise (t = 3.1) and deformation maps without correlated noise (t = 3.2), respectively.

Fig. 2. Simulation data. Black line represents fault. (a) annual average deformation rate in LOS direction; (b) Gaussian random noise with a standard deviation of 0.5 mm; (c)simulated orbit errors; (d) correlated noise with a standard deviation of 2 mm; (e) and (f) are synthetic deformation maps with correlated noise (t = 3.1) and without correlated noise (t = 3.2), respectively.

Four schemes are designed,as shown in Table 1.According to the deformation rate calculated by the two methods(the conventional and modified SBAS methods), we calculate the root mean squared error (RMSE) and the improved percent of the modified SBAS method compared with the conventional method, and comprehensively discuss the feasibility and robustness of the modified SBAS method.In schemeI,the standard deviation of the correlation noise with 2 mm is added to the N deformation maps. The feasibility and robustness of the conventional and modified SBAS methods are compared with the increase of N. In schemeII -Ⅳ, in which N = 5, N = 10, and N = 15, respectively, the synthetic deformation maps contain correlated noises with different standard deviations.

Table 1 Four schemes in the simulation experiment.

The RMSEs of the two methods in schemeI are shown in Fig.3(a).As N increases, the RMSEs of the conventional SBAS method gradually increase, while those of the modified SBAS method are relatively stable. Correspondingly, the improved percentage gradually increases. The results indicate that the modified SBAS method is superior to the conventional method, which has certain robustness when there are some low-quality deformation maps,and can obtain more accurate deformation rate. Fig. 3(b)-(d) show the RMSEs for the other three schemes,and the results are almost consistent.With the increase of the standard deviation of the correlation noise in N deformation maps, the RMSEs of the conventional SBAS method gradually increase, while those of the modified SBAS method are relatively stable,and the improved percent gradually increases.The results of the three schemes indicate that the larger the noises, the more obvious the advantage of the modified SBAS method. These simulation experiments demonstrate that the modified SBAS method can effectively reduce the influences of the residual errors when obtaining interseismic deformation, and then verify the feasibility and robustness of the proposed method.

4. Interseismic deformation rate of Haiyuan fault system

4.1. SAR data and processing

Fig.3. The RMSEs of two methods,and the improved percent of the modified SBAS method compared with the conventional SBAS method in four schemes.(a)-(d)for scheme I,II,III and Ⅳ, respectively.

Table 2 Basic parameters of Sentinel-1ASARdata.

We collect 15 Sentinel-1ASAR images from January 2018 to May 2020.The image range is shown in Fig.1,and the parameters of SAR data are in Table 2. Interferograms are generated using GAMMA software [34]. To ensure the better interferometric coherence,spatial perpendicular baselines of less than 100 m are used. To improve the magnitude of deformation signal, image pairs with temporal baselines of more than 480 days are employed.During the processing,the Shuttle Radar Topography Mission Digital Elevation Model (DEM) [35] is used to improve the SAR images’ coregistration accuracy and remove the influence of the terrain phase. To improve the interferometric coherence, interferograms are multi-looked with a factor of 20:4 and filtered using an adaptive filter function[4,15-18].Phase unwrapping is conducted based on the minimum-cost flow algorithm [4,15-18]. A total of 77 interferograms are generated. Finally, 36 interferograms with better coherence are selected to participate in the subsequent deformation rate estimation.The spatiotemporal baseline is shown in Fig.4.

The conventional and modified SBAS methods are used to calculate the deformation rate with the 36 interferograms. When using the modified SBAS method, five interferograms(20180121-20180720, 20180521-20180720, 20180521-20190727,20180930-20190925, and 20190528-20200522) are judged to have significant non-deformation signals, such as atmospheric artifacts and terrain errors, and then five system parameters (S1, S2,S3, S4, S5) are added in the equation (10).

Fig. 4. Spatiotemporal baselines of the interferograms.

4.2. Results

Fig. 5 shows the interseismic deformation rate of the Haiyuan fault system based on the modified SBAS method.According to the Fig. 5, it can be found that: (1) Most areas have continuous deformation distribution; there are incoherence in the desert areas in the northern part of the deformation map and in the near-fault area of the MMS segment. (2) The movements of the north and south sides of the Haiyuan fault system are different. The deformation rate is mainly negative on the south side,and mainly positive on the north side.(3)The deformation rate of the Haiyuan fault system is mainly between -4 and 4 mm/yr. From west to east, the deformation rate shows an increasing trend in the north side and a decreasing trend in the south side.(4)The deformation rates in the epicentre area of the 1920 Haiyuan and 1927 Gulang earthquakes are negative, which are different from the whole deformation characteristics on the north side of the fault.

5. Discussion and analysis

Fig. 6 shows the interseismic deformation rate of the Haiyuan fault system obtained by the conventional and modified SBAS methods.The distributions of the deformation rate fields obtained by two methods are generally consistent.The deformation rate is in the range of-4-4 mm/a.The deformation rate on the south side of the fault is mainly negative, indicating that the south side moves away from the satellite. The rate on the north side is positive,indicating that the northern side moves toward the satellite. This spatial distribution pattern is consistent with the characteristics of the left-lateral strike-slip movement of the Haiyuan fault system[2,14]. The difference between the two deformation rate fields mainly lies in the near fault zone.On the south side of the Haiyuan segment, the deformation rate field obtained by the conventional SBAS method shows a positive value in partial areas near the fault.In the deformation rate field obtained by the modified SBAS method, the positive value area is significantly reduced, which is more consistent with the general characteristics of the cross-fault deformation rate. This result is similar to the deformation pattern obtained by geological, GPS, and other observations [8-14,20,36].

Fig. 5. Interseismic deformation rate field (unit: mm/a) of Haiyuan fault system obtained by the modified SBAS method. Yellow line shows the Jingtai pull-apart basin boundary; black lines represent faults; red stars represent the epicenters of the 1920 MS 8.5 Haiyuanand 1927 MS 8.0 Gulang earthquakes, respectively; the profiles A-A′,B-B′ and C-C′ are used to display the cross-fault deformation in Fig. 7.

To further analyze the feasibility of the modified SBAS method,we calculated the standard deviation of the deformation rate in the far-field area. For the interseismic deformation signal, the further away from the fault, the smaller the deformation gradient; the deformation signal in the far-field sample area is generally consistent[19].A rectangular area is selected as the sample area in the north part of the MMS-LHS section (Fig. 6), and the standard deviation of the deformation rate in the area is calculated. The standard deviations of the conventional and modified SBAS methods in the sample area are 0.54 mm and 0.50 mm, respectively. To a certain extent, the result indicates that the modified SBAS method is more feasible and robust than the conventional SBAS method.

To intuitively display the characteristics of the deformation rate across the Haiyuan fault system, three profiles (Fig. 5) are drawn.The profile A-A′is located at the middle section of the LHS segment,and the profiles B-B′and C-C′are located at the western and middle sections of the Haiyuan segment, respectively. The profile is 4 km wide and 80 km long. The moving window is 4 km × 4 km, and the mean value and standard deviation are calculated in the window, respectively. The results are shown in Fig. 7.

Fig. 7(a)-(c) show the deformation rate of the profiles A-A′,B-B′and C-C′. The rate of the profile A-A′shows a decreasing trend from north to south,varying from 2 mm/a to-2.1 mm/a,with certain fluctuations. The rate of the profile B-B′changes from 1.9 mm/a to-1.3mm/a from north to south.The rate of the profile C-C′changes from 1.7 mm/a to-1.6 mm/a from north to south.In three profiles, there are no sudden variations in the near-field region across the fault, and the rate variations in the far-field region are relatively gentle, which are consistent with the interseismic deformation characteristics of locking fault [1,2,36]. This indicates that the middle section of the LHS segment and the middle and western sections of the Haiyuan segment are locked. The area 5-18 km away from the fault in the south part of the profile C-C′has a funnel-shaped change trend (Zone 1 in Fig. 7(c)), with a maximum deformation rate of approximately -2.8 mm/a. This funnel-shaped change may be related to the shallow surface deformation caused by local coal mining [2,3].

To further reveal the deformation rate characteristics along the Haiyuan fault system,the other four profiles are drawn(Fig.8(a)).The profiles X-X′and W-W′are located at the north and south sides of the MMS-LHS segment,and the profiles Y-Y′and Z-Z′are located at the north and south sides of the Haiyuan segment. The deformation ratechanges of the four profiles are shown in Fig.8(b)and (c).

In Fig. 8 (b), the profiles X-X′and W-W′show that the deformation rate of the north side of the MMS-LHS section increases from west to east, varying from 0 mm/a to 1.9 mm/a, with some regional fluctuations.The deformation rate of the south side in the MMS segment shows a decreasing trend, and in the LHS segment increases first and then decreases. Based on the deformation characteristics in the profiles X-X′and W-W′, the rate difference between the north and south sides shows an increasing trend in the MMS segment, and decreases first and then increases in the LHS segment.In Fig.8 (c),the profile Z-Z′shows that the deformation rate in the south side of the Haiyuan segment decreases from west to east. The profile Y-Y′shows that the deformation rate in the western section of the Haiyuan segment generally decreases first and then increases.Based on the deformation characteristics in the profiles Z-Z′and Y-Y′, the deformation difference between the north and south sides of the Haiyuan segment first decreases and then increases from west to east, which is consistent with the slip rate characteristics obtained by other techniques [2,8-14]. The profiles X-X′and Y-Y′show that the deformation rate of the Jingtai pull-apart basin is higher than that of the surrounding area(Zone 2 in Fig. 8 (b), Zone 3 in Fig. 8 (c)), indicating that the Jingtai pullapart basin is more active [4,36]. According to the recurrence interval is about 1000 years in this Tianzhu gap,and about 930 years have passed since the last major earthquake [4,6], therefore, one can speculate that the Jingtai pull-apart basin in the east of the Tianzhu gap has a higher seismic risk. Variations of deformation rate on the profiles reflect the complex movement trends among different segments of the Haiyuan fault system.

Fig. 6. Interseismic deformation rate field(unit:mm/a) of Haiyuan fault system obtained by the conventional(a)and modified SBAS methods (b). Black lines represent faults and the red dashed box shows the sample area.

Fig.7. Interseismic deformation rate of the cross-fault profiles at different locations of the Haiyuan fault system.(a)profile A-A′;(b)profile B-B′,and(c)profile C-C′.Black dashed line in (a)-(c) indicates the fault. Red circle in (c) shows the funnel-shaped deformation.

The standard deviations of the seven profiles (Figs. 7 and 8)indicate that the deformation rates are consistent in the calculation window. The standard deviation of the funnel-shaped area in the south side of the Haiyuan segment is relatively large (Zone 1 in Fig. 7 (c)). This may be due to a funnel-shaped rate change in this area,which leads to a significant difference in the deformation rate field.

Fig. 8. Location and deformation rate of profiles. (a) Location of four profiles (W-W′, X-X′, Y-Y′ and Z-Z′) along the approximately fault-parallel direction; (b) Interseismic deformation rate of the profiles W-W′and X-X′;(c)Interseismic deformation rate of the profiles Y-Y′and Z-Z′.Yellow line in(a)shows the Jingtai pull-apart basin boundary;red starsin (a) represent the epicenters of the 1920 MS8.5 Haiyuanand 1927 MS8.0 Gulang earthquakes, respectively. Black lines in (a) represent faults.

6. Conclusion

We propose a modified SBAS method based on the conventional SBAS method and the adjustment model with systematic parameters. The feasibility and robustness of the method are verified by simulation experiments.This modified SBAS method can reduce the negative influence of the residual signals related to terrain errors,atmospheric artifacts, orbit errors and noise, thereby obtaining a more accurate deformation rate map. We apply this method to extract the interseismic deformation rate of the Haiyuan fault system from 2018 to 2020. The interseismic deformation rate field shows significant differences in the deformation rate of the north and south sides of the Haiyuan fault system;this result is consistent with the characteristics of the left-lateral strike-slip movement.The deformation rate field and profiles reflect the complex movement trends among different segments of Haiyuan fault system in detail.The Jingtai pull-apart basin is more active than the surrounding areas.The middle section of the LHS segment and the western and middle sections of the Haiyuan segment is in a locking state.These results can provide more information to examine the locking depth,strain accumulation state, regional seismic risk of different segments of the Haiyuan fault system,and the deformation mechanism of the north-eastern margin of the Qinghai-Tibet Plateau.

Author statement

Yang Liu: Conceptualization, Methodology, Investigation, Data curation, Formal analysis, Writing - Original Draft.Sen Han:Methodology,Investigation,Data curation,Formal analysis,Writing-Original Draft.Luyun Xiong:Methodology,Data curation,Formal analysis.Ying Peng:Methodology,Formal analysis.Mengyao Gao:Data curation, Formal analysis.Caijun Xu: Methodology, Formal analysis.

Conflicts of interest

The authors declare that there is no conflicts of interest.

Acknowledgements

This research was supported by the National Natural Science Foundation of China(41874011,41861134009),and the National Key Research and Development Program of China (2018YFC1503603).The authors would like to thank the European Space Agency for providing the Sentinel-1A data.