Three-Dimensional Surface deformation of the 2022 Mw 6.6 Menyuan, China, Earthquake from InSAR and GF-7 Stereo Satellite Images

1. Introduction

On January 8, 2022, an Mw 6.6 earthquake occurred in Menyuan County, Qinghai Province, China, where has suffered from a series of moderate to strong earthquakes, such as the 1986 Mw 6.0 and 2016 Mw 5.9 Menyuan earthquake, Figure 1 [1]. The epicenter of the 2022 Menyuan earthquake was reported to (101.260°E, 37.770°N) by the China Earthquake Networks Center (CENC), and the depth was of 10km. Focal mechanism solutions of this earthquake indicated it was predominantly a left-lateral slip event. This earthquake ruptured the Lenglongling fault (LLLF) and the Tuolaishan fault (TLSF), Figure 1 [1], which belonged to the Qilian-Haiyuan fault system on the northeastern margin of the Tibetan Plateau [2,3]. The movement patterns in this region were mainly exhibited by strike-slip faults (Qilian-Haiyuan fault) and thrust faults (North Tuolaishan fault, Sunan-Qilian fault, North Qilian thrust faults) with strong earthquake activities, such as the 1920 Mw 7.9 Haiyuan earthquake [4] and 1927 Mw 8.0 Gulang earthquake [5].

Various studies have been carried out by researchers after the 2022 Mw 6.6 Menyuan earthquake. The field work found that the earthquake caused severe surface ruptures, which mainly produced surface cracks, as well as the collapse of the Lanzhou-Xinjiang Railway bridge and tunnel. The maximum coseismic horizontal displacement measured in the field were 2.6-3.7 m [1,6,7,8]. In addition, newly published detailed field survey indicated that this earthquake broke the junctions of LLLF with Sunan-Qilian fault and North Tuolaishan fault respectively [6]. And a secondary fault developed on the north of the middle rupture segment of the LLLF and produced 0.47m vertical displacement [8]. The above abundant field investigations demonstrated that the ruptures of the Menyuan earthquake with moderate magnitude were extremely complex and the amount of surface displacement were prominent. [9,10,11] processed the Sentinel-1 and Advanced Land Observation Satellite-2 (ALOS-2) data to obtain the line-of-sight (LOS) and azimuth (AZI) directions deformation based on the InSAR technology. The surface three-dimensional deformation fields have been solved with the method of SM-VCE just based on the SAR data, which indicated the maximum horizontal and vertical displacements were 1.9 m and 0.6 m, respectively [9].

Quantitative measurements of coseismic surface deformation have always been the focus of earthquake researches. Traditional field survey could accurately identify surface rupture traces and measure the surface displacements recorded by offset landforms. Especially, the developments of high-resolution optical satellites, Light Detection and Ranging (LiDAR), and other remote sensing technologies have greatly alleviated and enriched field measurements [12,13,14,15]. In addition, the geodetic technology of InSAR allowed for rapid acquisition of surface deformation fields after an earthquake. However, InSAR was bounded by incoherence, especially in the areas of coseismic rupture fault zone. The absent of near-fault displacement information directly induced to the discrepancy between field measurements and InSAR observations. Global Navigation Satellite System (GNSS) observations were so sparse to capture the complete deformation field across the full fault rupture zone. Recently, the development of optical pixel correlation technique based on high-resolution optical remote sensing data has enabled the acquisition of near-fault deformation fields [16,17]. A number of studies based on the near-field deformation field obtained by optical pixel correlation technique have revealed that the displacements occurring on the seismic rupture zone contained not only the displacements which localized to within tens of meters of the surface traces of the faults, but also the distributed off-fault deformation (OFD) that extended tens to hundreds of meters away from the surface traces and were difficult to be acquired by the InSAR or field measurements [18,19,20,21,22,23,24]. Our previous research has acquired the horizontal deformation field of the 2022 Mw 6.6 Menyuan earthquake based on GaoFen-2/7 stereo satellite images, and performed quantitative measurements for the OFD [24].

In this study, we comprehensively solved the complete surface three-dimensional deformation fields of the 2022 Mw 6.6 Menyuan earthquake by incorporating the published observations of InSAR [9], optical pixel correlation [24], and differencing DEMs calculated by this study. Firstly, we extracted the pre- and post- 2022 Mw 6.6 Menyuan earthquake DEMs using the GF-7 stereo satellite images, and obtained the vertical deformation field by DEMs differential calculation. Then, taking into account the InSAR LOS and AZI directions deformation obtained from Sentinel-1 and ALOS-2, horizontal deformation calculated based on optical pixel correlating to pre- and post-earthquake GF-2/7 images, and the vertical deformation computed from DEMs differencing, we established the geometric relationships between the above observation directions and the surface three-dimensional directions, and solved the complete three-dimensional deformation fields of the 2022 Mw 6.6 Menyuan earthquake. In addition, the 2D strain field were derived from the surface horizontal displacement maps calculated based on optical pixel correlation method. Finally, the kinematic characteristics of the 2022 Menyuan earthquake were obtained by analyzing the three-dimensional deformation fields and strains.

2. Data and Methods

2.1. Data

We collected the pre- and post- earthquake GF-7 stereo images data from High-resolution Remote Sensing Data Center, China Earthquake Administration (Table 1) to extract the pre- and post-earthquake DEMs for the 2022 Mw 6.6 Menyuan earthquake. In addition, the LOS deformation of DInSAR, and AZI deformation of pixel offset-tracking (POT), multiple aperture InSAR (MAI) and burst overlap InSAR (BOI), respectively, were used too [9]. At the same time, we incorporated the near-field horizontal east-west deformation from [24] and the vertical deformation processed by this study. The detail data were listed in Table 2.

2.2. Calculation of Vertical Displacement

The processes of vertical deformation field calculation of the 2022 Menyuan earthquake based on the pre- and post-earthquake DEMs were described below. Firstly, the high-density point clouds data were obtained by processing the GF-7 stereo satellite image pairs of before and after earthquake, respectively. Then, the point clouds of post-earthquake were aligned to that of pre-earthquake using the method of Iterative Closest Point (IPC), which was designed to process LiDAR point cloud [25] and could be achieved to by the open-source software Cloud Compare (https://www.cloudcompare.org/). Next, the point clouds were meshed to 1x1 m spatial resolution DEMs. It should be noted that the terrain changes were attributed to both horizontal and vertical displacements caused by earthquake. Therefore, in order to get the real vertical displacement, the post-earthquake DEM needed to be corrected in terms of the coseismic horizontal displacement, i.e., the horizontal displacements result of the optical pixel correlation of GF-2/7. Finally, the coseismic vertical deformation filed was obtained by subtracting the pre-earthquake DEM from the rectified post-earthquake DEM.

2.3. Calculation of Coseismic Surface Three-Dimensional Deformation Fields

According to the geometry of SAR satellite observation, the relationships between the LOS and AZI direction displacements observed by InSAR with the three-dimensional deformation of the ground surface were presented in Equations (1) and (2), respectively. In addition, it was set that the horizontal E-W displacement calculated by method of optical pixel correlation based on the GF-2/7 images was the coseismic surface E-W displacement. And the vertical displacement calculated by the DEMs difference was the coseismic surface vertical displacement. Based on the deformation data listed in Table 2, we constructed the matrix relationships between the different directional deformation obtained by each observation means and the ground surface three-dimensional deformation. Finally, by solving the matrix Equation (3) via the least-squares method, the complete surface three-dimensional deformation fields with incorporating both the near- and far-field were derived.

$$d_{LOS}=d_U\cos \theta -\sin \theta \left[d_N\cos \left(\alpha -3\pi /2\right)+d_E\sin \left(\alpha -3\pi /2\right)\right]$$

(1)

$$d_{AZI}=d_N\sin \left(\alpha -3\pi /2\right)-d_E\cos (\alpha -3\pi /2)$$

(2)

$d_{LOS}$ and $d_{AZI}$ were line-of-sight and azimuth direction displacement by InSAR observation, respectively. $\theta$ and $\alpha$ were the angels of incidence and azimuth of the SAR satellite observation, respectively. $\left[d_E{d}_N{d}_U\right]$ were ground surface east-west, north-west and vertical direction displacement, respectively.

$$\left[\begin{array}{ccc}-{\sin \left(\alpha -3\pi /2\right)}_{LOS}& -{\sin \theta \cos (\alpha -3\pi /2)}_{LOS}& {\cos \theta }_{LOS}\\ \dots & \dots & \dots \\ {-\cos (\alpha -3\pi /2)}_{AZI}& {\sin (\alpha -3\pi /2)}_{AZI}& 0\\ \dots & \dots & 0\\ 1_{GF-EW}& 0& 0\\ 0& 0& 1_{GF-UN}\end{array}\right]\left[\begin{array}{c}{d}_E\\ d_N\\ d_U\end{array}\right]=\left[\begin{array}{c}{d}_{LOS}\\ \dots \\ d_{AZI}\\ \dots \\ d_{GF-EW}\\ d_{GF-UN}\end{array}\right]$$

(3)

3. Results

3.1. The Vertical Displacement from DEMs Difference

Figure 2a illustrated the vertical deformation field of the 2022 Menyuan earthquake obtained by differencing the pre- and post-earthquake DEMs. This map explicitly demonstrated that the southern sides of the rupture faults were uplifted relative to the northern sides of that. Along the main coseismic ruptures of the LLLF, there were two regions with significant vertical displacements. One region was at the fault bend on the western segment of main ruptures, where the notable vertical displacements were attributed to the local compressive stress induced by the right-step left-lateral strike slip faults (Figure 2b). Another region with pronounced vertical displacements was at the epicenter on the middle segment of the main ruptured LLLF, where suffered serious ground damages and collapse of the Lanzhou-Xinjiang Railway tunnel and bridge.

We extracted the displacement profiles across both sides of the rupture fault. Then the relative vertical displacements of the southern side versus to the northern side of the rupture fault were calculated. The calculation result shown that the vertical displacement at the fault bend region on the western segment of main rupture fault (LLLF) were 1.25-1.43 m (Figure 2c-e). At the region of the epicenter of the ruptured middle segment, the cross-fault vertical profiles revealed relative vertical displacement occurred on both the main rupture fault and its northern secondary fault simultaneously. Specifically, the relative vertical displacement along the main rupture on the LLLF were 1.08-2.30 m (Figure 2f, g). The relative vertical displacement along the secondary fault to the northern of the main rupture were 0.36-0.76 m (Figure 2f, g), which were consistent with the field measurements [8]. In addition, on the eastern end of the secondary rupture, we still measured a relative vertical displacement of 1.01 m.

3.2. The Surface Three-Dimensional Deformation Fields

Figure 3 demonstrated the resolved results of the coseismic surface three-dimensional deformation fields of the 2022 Mw 6.6 Menyuan earthquake by incorporating multiple remote sensing data, in which the east-west (E-W) deformation was the largest reaching to 2 m, and the north-south and vertical deformations were approximately 1 m. The east-west deformation presented in Figure 3a revealed a homogeneous distribution pattern. In contrast, the north-south deformation presented in Figure 3b shown an inhomogeneous distribution pattern. In the western segment of the coseismic ruptures along the ruptured LLLF and eastern end of the ruptured TLSF, both sides of the rupture fault appeared to move in a northward direction. While, in the middle and eastern segments of the ruptured LLLF, the southern side of the rupture fault displayed obvious southward movement. Figure 3c showed the vertical deformation feature with two contrasting centers of uplift-subsidence. Specifically, the southern side along the coseismic rupture of the TLSF exhibited subsidence. While the southern side along the main rupture of the LLLF exhibited uplift, which was consistent with the result of the coseismic vertical displacement from the DEMs differential calculation in Section 3.1. The arrows in the Figure 3d demonstrated the horizontal displacements which were composed by the east-west and north-south displacements, showing a clockwise rotational movement by the relative left-lateral slipping of the two sides of the rupture fault. For the quantitative analysis of the three-dimensional deformation fields of the 2022 Mw 6.6 Menyuan earthquake which incorporated the near- and far-field displacements, as well as their comparing with other results [9,24], could refer to the Section 4.1.

4. Discussion

4.1. Characteristics of the Near-Field Deformation

In this paper, we obtained the complete surface three-dimensional fields of the 2022 Menyuan earthquake by incorporating the near- and far-field deformation data (Figure 3). We expected to learn a more comprehensive understanding from the complete deformation fields than traditional InSAR observation. Therefore, we selected the E-W direction displacement with the most significant deformation information from the complete three-dimensional fields, and compared it with the displacement results from InSAR and optical pixel correlation, respectively (Figure 4). Firstly, the displacement profiles across both sides of the rupture fault within a distance of 1-1.5km (Figure 4d-i) demonstrated that the results from optical pixel correlation were able to better constrain the near-fault deformation that occurred within tens to hundreds of meters away from the surface rupture traces, which were depicted by black displacement points in the Figure 4d-i. However, the InSAR technique cannot constrain the near-fault deformation owing to the incoherence, and thereby led to the inability of the InSAR technique in determining the location of the coseismic rupture traces and the width of the fault rupture zone. Secondly, in terms of the trend of displacement distributions, the displacement results from optical pixel correlation, especially on the near-fault regions, were larger than the measurement results from InSAR. On one hand, the cause was that InSAR was insufficient for near-field deformation measurements. On the other hand, the deformation on the rupture fault zone of the coseismic earthquake contained not only the offset localized on the primary fault (on-fault displacement), but also diffuse displacements distributed on both sides of the primary fault (off-fault displacement) [24], e.g., penetrative continuous shear, block rotation accompanied by internal antithetic shear, and small block rotation [26,27]. The above factors well explained that the near-fault displacement measurements from optical pixel correlation were larger than the InSAR, and even the field measurements [24]. Finally, the solution results by jointing the InSAR and optical pixel correlation deformation data, i.e., the red displacement points in Figure 4d-i, lied between the InSAR and optical pixel correlation deformation data, which indicated that both data provided certain constraints on the joint calculation results. However, since in this study, only one optical deformation data was used, while the InSAR deformation results were abundant, the joint calculation results were more influenced by the InSAR data. In order to better constrain the near-field deformation, sufficient optical remote sensing data should be used to extract the coseismic near-field deformation.

4.2. Coseismic Strain

Compared to the displacement maps (Figure 3), the strain fields, which were derived from the displacement maps, could provide more insights into the crustal deformation. Generally, the geodetic techniques (GNSS and InSAR) were utilized to calculate the interseismic deformation, and then it was used to generate the interseismic strain rate field. Based on the accumulation and difference in strain rate distribution on the seismological faults, the degree of seismic urgency in different segments could thus be identified [28,29]. In order to quantify the coseismic strain, the strain tensor can be derived from the surface displacement gradient, and then the dilatation, rotation and shear strain fields could be obtained by decomposing the strain tensor. The coseismic strain invariants could reflect more specific kinematic characteristics and rupture deformation patterns. The dilatation strain indicated the contraction and extension deformation. The rotation demonstrated the rotational deformation. The maximum shear strain characterized the region with significant relative displacement mostly concentrated on the rupture zone [30].

In this study, the horizontal E-W and N-S displacements calculated from optical pixel correlation method, which were down-sampled to a spatial resolution of 3 x 3 m, were employed to calculate the 2D strain fields of the 2022 Mw 6.6 Menyuan earthquake. The software written by [31] was used, and it was set to be a calculation window of 2 x 2 pixel, and a mimimum-moment regularization of 0.25. Figure 5 illustrated the strain maps. It should be noted that the InSAR deformation data with very low spatial resolution was not be used by this study for strain calculation. By contrast, [9] just used the InSAR deformation data to calculate the coseismic strain of the 2022 Mw 6.6 Menyuan earthquake. However, the strain results [9] cannot reflect the characteristics of the strain distribution on the coseismic rupture fault zone.

According to the Figure 5, the strain maps clearly demonstrated the concentration of high strains on the rupture fault zone. Specifically, the dilatation strain on the rupture zone shown negative, indicating that the coseismic deformation had extrusion component (Figure 5b). The negative values of the rotation strain (Figure 5c) indicated a relative clockwise motion of the two sides of the rupture fault, which was consistent with the motion trend illustrated in Figure 3d. As for shear strain (Figure 5d), it depicted the ground traces with prominent relative displacements on the rupture zone. That was, for large strike-slip earthquakes, the location of surface rupture zones and rupture traces could be clearly marked on the basis of the shear strain map. In addition, there were some special characteristics for the strain distribution of the Menyuan earthquake coseismic rupture zone. Firstly, in terms of the strain distribution along the rupture fault, the signs of positive and negative values of strain component along the main rupture zone did not change, suggesting that the stress state was not subjected to reverse changes due to variations in stress environment or local geometry structure. Secondly, along the eastern segment of the main rupture of the LLLF, the strains distributed dispersedly, especially where the ruptures passed through the Liuhuang River. Thirdly, the coseismic shear strain calculated in this study were much larger than the commonly assumed upper limit of 0.5% yield strain for rocks in the laboratory [32]. In fact, the factors influencing the ruptures of the fault zone were very complicated. With the research development, the concept of fault structural maturity has been proposed. When strong earthquakes occurred on fault zones with different maturity levels, their surface deformation and strain distributions would behave different characteristics. For example, the 2021 Mw 7.4 Maduo earthquake were considered to occurred on immature faults with coseismic shear strain less than 0.5% of the laboratory threshold [33,34]. Further, [34] proposed a coseismic rupture strain threshold of 0.5-1.5% by measuring the shear strain derived from the deformation field calculated from the SPOT optical satellite images for the 2021 Maduo earthquake. According to our early study [24], the seismogenic fault (i.e., the LLLF) of the 2022 Menyuan earthquake was a mature fault, from which it could be assumed that the large coseismic shear strain occurred on the 2022 Menyuan earthquake surface rupture zone seemed reasonable. As for the assessment of coseismic strain thresholds for mature and immature ruptures, the supports of more seismic cases were definitely needed.

5. Conclusions

In this study, we comprehensively solved the surface three-dimensional deformation fields of the 2022 Mw 6.6 Menyuan earthquake by combining the surface deformation field solved by multi-source remote sensing data and technological means, i.e., the LOS and AZI directions deformation calculated by InSAR based on the Sentinel-1 and ALOS-2 data, the surface horizontal deformation field solved by optical pixel correlation based on the GF-2/7 satellite data, and the vertical deformation field calculated by DEMs difference based on the pre- and post-earthquake GF-7 stereo image pairs data. Firstly, the pre- and post-earthquake topographic point cloud data were extracted based on the GF-7 stereo image pair data. Then the pre- and post-earthquake point clouds were aligned and gridded into DEMs with a spatial resolution of 1 x 1 m. After removing the impact of coseismic horizontal displacement, the coseismic vertical deformation field was obtained by DEMs difference. The coseismic vertical deformation field revealed two regions with significant vertical displacements. One of them exhibited vertical displacement due to the local compressive stress at a bend of right-step left-lateral strike slip fault. The other was the epicenter area, where showed significant uplift to the south of both the main rupture fault and the secondary fault. Secondly, this study demonstrated that the optical pixel correlation method based on the very-high spatial resolution satellite images could provide exhaustive constraint on the near-field deformation than the InSAR technique. The main reasons were that InSAR was subject to incoherence and the low resolution was insufficient to characterize the deformation distribution on the rupture zone. Therefore, to obtain a complete surface three-dimensional deformation field, abundant high spatial resolution optical images were needed to provide near-field deformation constraints. Finally, the coseismic strain fields were derived based on the high spatial resolution deformation field calculated by the optical pixel correlation method. The coseismic strain field showed high strain concentration on the coseismic rupture zone, especially on the shear strain map, which could clearly recognize the surface rupture trace. In conclusion, it can be concluded from our research that a single technical means or data source introduced limitations in understanding the characteristics and patterns of coseismic surface rupture deformation. Therefore, the application of diverse research methods and data sources should been applied and could bring breakthroughs in the study of earthquake rupture deformation.