Fracture blocking rate inversion method based on GNSS earth crust deformation observation and elastic dislocation theory constraint

Abstract

The invention provides a fracture blocking rate inversion method based on GNSS earth crust deformation observation and elastic dislocation theory constraint, and the method comprises the steps: constructing a fracture three-dimensional model and an elastic dislocation constraint model, and calculating a fracture blocking rate. Through multi-period result comparison, strain field verification and earthquake activity verification, a high-risk region of pregnancy and earthquake can be identified. According to the method, the GNSS earth crust deformation observation elastic dislocation theory is adopted to constrain the negative dislocation model, and the precision and stability of fracture blocking rate inversion are remarkably improved. According to the invention, multiple periods of GNSS data are integrated to realize dynamic inversion, and the evolution trend of the locking state can be captured. According to the method, secondary blocks are finely divided, strain rate field constraints are introduced, locking distribution of a local strain enhancement area is accurately identified, the reliability of judgment of a large earthquake risk area is improved, a fracture high-locking-rate area is defined, through theoretical constraint and multi-source data fusion, the fracture locking rate inversion precision and the dynamic capture capability are improved, and the method is suitable for large-scale popularization and application. And reliable technical support is provided for earthquake risk assessment.

Description

Technical Field

[0001] This invention relates to the field of earthquake monitoring and crustal deformation analysis technology, and more specifically, to a method for inverting fracture locking rate under the constraints of GNSS crustal deformation observation and elastic dislocation theory. Background Technology

[0002] Fault locking status is a core indicator for assessing regional seismic risk. The accumulation of slip loss and stress concentration in locked fault segments are directly related to the likelihood of major earthquakes. Monitoring the locking characteristics of major active faults in the North-South Seismic Belt is of great significance for earthquake disaster early warning.

[0003] Existing methods for inverting fault locking rates largely rely on single-period GNSS velocity field data constraints, lacking the rigorous constraints of elastic dislocation theory, resulting in insufficient stability and accuracy of the inversion results. Furthermore, they lack sufficient comparative analysis of slip parameters across different time periods, making it difficult to capture the dynamic evolution characteristics of the locking state. In addition, traditional methods do not provide sufficiently fine division of secondary blocks, failing to accurately reflect the locking distribution patterns in areas of increased local strain, thus limiting the reliability of seismic risk assessment.

[0004] The China Mainland Tectonic Environment Monitoring Network has accumulated long-term, continuous GNSS observation data, providing a high-precision data source for crustal deformation monitoring. How to construct a scientific inversion model based on this data and elastic dislocation theory to achieve dynamic and high-precision inversion of fault locking rates has become a pressing technical problem in the field of earthquake monitoring. Therefore, this paper proposes a fault locking rate inversion method based on GNSS crustal deformation observations and elastic dislocation theory constraints. Summary of the Invention

[0005] The purpose of this invention is to achieve high-precision inversion of fracture slip rate, slip deficit, and lockout rate by integrating multi-period GNSS velocity field data constraints and utilizing negative dislocation models and elastic dislocation theory [Note: In a sense, negative dislocation models and elastic dislocation theory are two expressions of the same thing!]. Simultaneously, it reveals the dynamic evolution characteristics of the lockout state, providing a reliable basis for earthquake risk assessment. To achieve the above-mentioned objective, this invention provides the following technical solution: a fracture lockout rate inversion method based on GNSS crustal deformation observations and elastic dislocation theory constraints, comprising the following steps: Step (1) GNSS observation data acquisition: Collect multi-period observation data from GNSS reference stations and rover stations covering the target area of ​​the China Continental Tectonic Environment Monitoring Network. The target area is the northern segment of the North-South Seismic Belt, which is 32°~43°N and 95°~112°E. Step (2) Data preprocessing: The GNSS observation data is subjected to gross error removal based on spatiotemporal correlation, and high-precision atmospheric delay correction is performed by combining the global atmospheric grid model and regional water vapor radiometer data; the influence of coseismic deformation is stripped off by InSAR coseismic deformation field constraint, and a clean GNSS deformation time series with both time resolution and accuracy is generated. Step (3) GNSS velocity field calculation: Clean observation data is calculated using the MIT open-source software GAMIT / GLOBK to construct a multi-period GNSS crustal movement velocity field under a unified reference frame (Eurasian plate stability frame); Step (4) Three-dimensional geometric modeling of the fault: Based on the geological survey results and previous research data, establish a three-dimensional geometric model of the main active faults in the target area, and clarify the fault strike, dip, dip angle and spatial distribution of 0-30km depth; Step (5) Construction of elastic dislocation constraint model: The fault-locked segment is equivalent to a rectangular negative dislocation unit embedded in an elastic half-space medium, and the Green's function response matrix between the fault negative dislocation slip deficit and the surface GNSS deformation rate is established; on this basis, a dual physical constraint mechanism is introduced: the strain rate tensor inside the secondary block calculated based on the block dynamics model is used as a priori constraint to ensure that the inversion results are consistent with the regional tectonic stress field; the historical earthquake recurrence interval, paleoseismic slip and earthquake gap distribution are integrated to construct a probability weight function and give higher inversion priority to high seismic hazard segments; Step (6) Regularized inversion of fracture parameters: The adaptive Tikhonov regularization method is adopted to automatically adjust the regularization factor according to the geometric complexity of the fracture segment, solve the constrained inversion equation, and simultaneously calculate the long-term slip rate, coseismic slip amount and post-earthquake viscoelastic relaxation rate of each fracture segment to obtain the accurate fracture slip loss rate. Step (7) Quantitative calculation of fracture blocking rate: Calculate the blocking rate of each fracture segment based on the ratio of fracture slip loss rate to long-term slip rate, and then calculate the blocking rate based on the fracture slip loss rate obtained from the inversion. Long-term fault slip rate estimated in conjunction with geological / geometry According to the formula Calculate the blocking rate of each segment and introduce an uncertainty propagation model to quantify the confidence interval of the blocking rate; Step (8) Result verification and risk identification: The reliability of the inversion results is confirmed by multi-period GNSS velocity field time series collaborative verification, InSAR deformation rate field cross-verification and historical strong earthquake recurrence cycle matching analysis; combined with the three-dimensional distribution of fault lockout rate and Coulomb stress loading rate, high-risk earthquake-prone areas are identified.

[0006] As a preferred technical solution of the present invention, the multi-period observation data in step (1) includes continuous observation data from 1998 to 2018, and data from 10 subdivided periods from 1999 to 2001, 2001 to 2004, ..., 2013 to 2015, and 2015 to 2020.

[0007] As a preferred technical solution of the present invention, the accuracy of the GNSS velocity field obtained in step (3) is better than 0.5 mm / a, and the velocity projection data of 18 baseline length time series and 5 profiles across the active fault are extracted.

[0008] As a preferred technical solution of the present invention, the main active faults mentioned in step (4) include the northern margin fault of the Altyn Tagh Mountains, the eastern foothills of the Helan Mountains, the Zhuanglanghe Fault, the northern margin fault of the western Qinling Mountains, the eastern Kunlun Fault, the Haiyuan Fault, the Liupanshan Fault, and the northern margin fault of the Qilian Mountains.

[0009] As a preferred technical solution of the present invention, the response relationship in step (5) is characterized by the elastic dislocation Green's function matrix, and the formula is: ,in This is the velocity field vector observed by GNSS. Let be the Green's function matrix, and d be the sliding parameter vector.

[0010] As a preferred technical solution of the present invention, the regularized inversion equation in step (6) is: ,in The regularization parameter is determined by the L-curve method and has a value of 0.01. L is the smoothing matrix.

[0011] As a preferred technical solution of the present invention, the formula for calculating the breakage lock-in rate in step (7) is as follows: ,in Where is the fracture lock-up rate, and D is the fracture slip deficit rate (mm / a). The long-term sliding rate of the fracture is (mm / a). The value range is 0 to 100%.

[0012] As a preferred technical solution of the present invention, the multi-period result comparison in step (8) includes calculating the change in fracture slip loss rate between the period 2015-2020 and the period 2013-2015, and the strain rate field co-verification includes the correlation analysis of the principal strain rate field, the maximum shear strain rate field and the surface strain field.

[0013] As a preferred technical solution of the present invention, when performing three-dimensional geometric modeling of the fracture in step (4), each fracture is divided into sub-segments with a length of 5 to 10 km, and the sub-segment division accuracy matches the distribution density of GNSS stations.

[0014] As a preferred technical solution of the present invention, the high-risk area for major earthquakes mentioned in step (8) includes 12 areas with a closure rate higher than 70%, such as the Subei Kuantanshan section of the northern margin fault of the Altyn Tagh, the Yongdeng County-Lanzhou border section of the Zhuanglanghe fault, the Tianshui section and the Luanfeng section of the northern margin fault of the West Qinling Mountains.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention employs a GNSS crustal deformation observation-constrained negative dislocation model, which significantly improves the accuracy and stability of fault locking rate inversion and solves the problem of lack of theoretical constraints in traditional methods.

[0016] This invention integrates multi-period GNSS data to achieve dynamic inversion, which can capture the evolution trend of the closure state over time and provide dynamic data support for medium- and long-term earthquake risk assessment.

[0017] This invention finely divides secondary blocks and introduces strain rate field constraints to accurately identify the locked distribution of local strain-enhanced regions, thereby improving the reliability of the determination of high earthquake risk areas.

[0018] The inversion results of this invention cover the major active faults of the North-South Seismic Belt and identify areas with high fault closure rates, providing clear target areas for earthquake monitoring and early warning. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the method flow provided by the present invention; Figure 2 A schematic diagram of the GNSS velocity field calculation results provided by this invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention.

[0021] Therefore, the following detailed description of the embodiments of the present invention is not intended to limit the scope of the claimed invention, but merely illustrates some embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention. It should be noted that, in the absence of conflict, the embodiments and features and technical solutions in the embodiments of the present invention can be combined with each other. It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0022] Example 1: The method for inverting the fault locking rate under the constraint of elastic dislocation theory based on GNSS includes the following steps: Step (1) GNSS observation data acquisition: Collect multi-period observation data of GNSS reference stations and rover stations covering the target area of ​​the China Mainland Tectonic Environment Monitoring Network. The target area is the northern segment of the North-South Seismic Belt of 32°~43°N and 95°~112°E. Step (2) Data preprocessing: The GNSS observation data is subjected to gross error removal based on spatiotemporal correlation, and high-precision atmospheric delay correction is performed by combining the global atmospheric grid model and regional water vapor radiometer data; the coseismic influence is removed by InSAR coseismic deformation field constraint, and a clean GNSS deformation time series with both time resolution and accuracy is generated. Step (3) GNSS velocity field calculation: Clean observation data is calculated using the MIT open-source software GAMIT / GLOBK to construct a multi-period GNSS crustal movement velocity field under a unified reference frame (Eurasian plate stability frame); Step (4) Three-dimensional geometric modeling of the fault: Based on the geological survey results and previous research data, establish a three-dimensional geometric model of the main active faults in the target area, and clarify the fault strike, dip, dip angle and spatial distribution of 0-30km depth; Step (5) Construction of elastic dislocation constraint model: The fault-locked segment is equivalent to a rectangular negative dislocation unit embedded in an elastic half-space medium, and the Green's function response matrix between the negative dislocation slip deficit and the surface GNSS deformation rate is established; on this basis, a dual physical constraint mechanism is introduced: the strain rate tensor inside the secondary block calculated based on the block dynamics model is used as a priori constraint to ensure that the inversion results are consistent with the regional tectonic stress field; the historical earthquake recurrence interval, paleoseismic slip and earthquake gap distribution are integrated to construct a probability weight function and give higher inversion priority to high seismic hazard segments; Step (6) Regularized inversion of fracture parameters: The adaptive Tikhonov regularization method is adopted to automatically adjust the regularization factor according to the geometric complexity of the fracture segment, solve the constrained inversion equation, and simultaneously calculate the long-term slip rate, coseismic slip amount and post-earthquake viscoelastic relaxation rate of each fracture segment to obtain the accurate fracture slip loss rate. Step (7) Quantitative calculation of fracture blocking rate: Calculate the blocking rate of each fracture segment based on the ratio of fracture slip loss rate to long-term slip rate, and then calculate the blocking rate based on the fracture slip loss rate obtained from the inversion. Long-term fault slip rate estimated in conjunction with geological / geometry According to the formula Calculate the blocking rate of each segment and introduce an uncertainty propagation model to quantify the confidence interval of the blocking rate; Step (8) Result verification and risk identification: The reliability of the inversion results is confirmed by multi-period GNSS velocity field time series collaborative verification, InSAR deformation rate field cross-verification and historical strong earthquake recurrence cycle matching analysis; combined with the three-dimensional distribution of fault lockout rate and Coulomb stress loading rate, high-risk earthquake-prone areas are identified.

[0023] The multi-period observation data mentioned in step (1) includes continuous observation data from 1998 to 2018, and data from 10 subdivided periods: 1999 to 2001, 2001 to 2004, ..., 2013 to 2015, and 2015 to 2020.

[0024] The accuracy of the GNSS velocity field obtained in step (3) is better than 0.5 mm / a. At the same time, the time series of 18 baseline lengths and the velocity projection data of 5 profiles across the active fault are extracted.

[0025] The main active faults mentioned in step (4) include the northern margin fault of the Altyn Tagh Mountains, the eastern foothills of the Helan Mountains, the Zhuanglanghe Fault, the northern margin fault of the western Qinling Mountains, the eastern Kunlun Mountains, the Haiyuan Fault, the Liupanshan Fault, and the northern margin fault of the Qilian Mountains.

[0026] The response relationship described in step (5) is characterized by the Green's function matrix of the elastic dislocation, and the formula is as follows: ,in This is the velocity field vector observed by GNSS. Let be the Green's function matrix, and d be the sliding parameter vector.

[0027] The regularized inversion equation mentioned in step (6) is: ,in The regularization parameter is determined by the L-curve method and has a value of 0.01. L is the smoothing matrix.

[0028] The formula for calculating the breakage lockout rate in step (7) is as follows: ,in Where is the fracture lock-up rate, and D is the fracture slip deficit rate (mm / a). The long-term sliding rate of the fracture is (mm / a). The value range is 0 to 100%.

[0029] The multi-period result comparison in step (8) includes calculating the change in fracture slip loss rate between the period 2015-2020 and the period 2013-2015. The strain rate field co-verification includes the correlation analysis of the principal strain rate field, the maximum shear strain rate field and the surface strain field.

[0030] When performing three-dimensional geometric modeling of the fracture in step (4), each fracture is divided into segments with a length of 5 to 10 km. The segment division accuracy matches the distribution density of GNSS stations.

[0031] The high-risk areas for major earthquakes mentioned in step (8) include 12 areas with a closure rate higher than 70%, such as the Subei Kuantanshan section of the northern margin fault of the Altyn Tagh, the Yongdeng County-Lanzhou border section of the Zhuanglanghe fault, and the Tianshui and Luanfeng sections of the northern margin fault of the West Qinling Mountains.

[0032] Example 2: Inversion Method for Fracture Lockout Rate under Elastic Dislocation Theory Constraints Based on GNSS (a) Data preprocessing steps GNSS data acquisition: Collect observation data from GNSS reference stations and rover stations of the China Mainland Tectonic Environment Monitoring Network, covering multiple periods of observation data in the target area (32°~43°N, 95°~112°E), including key periods such as 1998–2018 and 2015–2020, as well as subdivided periods such as 1999–2001 and 2001–2004.

[0033] Data processing: The MIT open-source software GAMIT / GLOBK was used to calculate the daily relaxation solutions for each site, constructing a multi-period GNSS velocity field under the Eurasian framework. At the same time, the Ordos block was used as a reference to construct velocity field data under a local framework, eliminating interference factors such as coseismic effects.

[0034] Data validation: The reliability of the velocity field data was verified by calculating the time series of the baseline length across the active fault. Five profiles covering the main faults were selected, and the velocity projections of each station within a range of 50 to 350 kilometers on both sides of the fault were extracted, which are perpendicular or parallel to the fault strike, to provide basic data for subsequent inversion.

[0035] (II) Construction of the Elastic Dislocation Theory Constraint Model Fault geometry modeling: Based on previous research, a three-dimensional geometric model of the main active faults in the North-South Seismic Belt is established, including the northern margin fault of the Altun Mountains, the eastern foothill fault of the Helan Mountains, the Zhuanglanghe Fault, the northern margin fault of the western Qinling Mountains, and the eastern Kunlun Fault, etc., clarifying the fault strike, dip, dip angle and spatial distribution within the 0-30km depth range.

[0036] Negative dislocation model setup: The fracture-locked segment is equivalent to a negative dislocation. Under the assumption of elastic half-space, the response relationship between negative dislocation and GNSS-observed deformation is constructed, as shown in the following formula:

[0037] in, This is the velocity field vector observed by GNSS. is the Green's function matrix of elastic dislocations, calculated based on elastic dislocation theory, and d is the fracture slip parameter vector, including fracture slip rate and slip deficit.

[0038] Constraint setting: Based on the results of secondary block division, the strain rate change within the block is used as a constraint to limit the spatial continuity of fracture sliding; historical seismic activity data is introduced to constrain the spatial distribution range of the locked segment.

[0039] (III) Fracture Parameter Inversion Steps Initial parameter estimation: Based on the fracture slip rate results of previous studies, the initial value range of the slip parameters is set, including the rate range corresponding to motion types such as left-handed strike-slip and reverse compression.

[0040] Regularized inversion calculation: The inversion equation is solved using the Tikhonov regularization method, balancing the data fit and model smoothness, to calculate the slip rate and slip deficit at depths of 0–30k for each fracture segment, as shown in the following formulas:

[0041] in, Here, L is the regularization parameter, and L is the smoothing matrix.

[0042] Fracture lockout ratio calculation: The fracture lockout ratio is defined as the ratio of fracture slip deficit to long-term slip rate. :

[0043] Where D is the fracture slip loss rate (mm / a). The long-term sliding rate of the fracture is (mm / a). The value ranges from 0 to 100%, with a higher value indicating a higher degree of locking.

[0044] (iv) Dynamic evolution analysis and verification Multi-period result comparison: Quantitative analysis was conducted on the results of fault slip rate, slip deficit and closure rate in different periods (such as 2013-2015 and 2015-2020), and the changes were calculated (such as the change in slip deficit rate of the Alak Lake-Toso Lake segment of the East Kunlun Fault) to reveal the trend evolution characteristics of the closure state.

[0045] Synergistic verification of strain rate field: By combining the data of principal strain rate field, maximum shear strain rate field and surface strain field, the spatial correspondence between the locked section and the strain enhancement region is verified; the stress accumulation in the locked section is corroborated by the trend of cross-fracture baseline length variation.

[0046] Seismic activity verification: Integrate historical earthquake data within the study area, analyze the correlation between areas with high fault lock-in rates and earthquake occurrence, and verify the reliability of the inversion results.

[0047] Example 3: Inversion Method for Fracture Lockout Rate under Elastic Dislocation Theory Constraints Based on GNSS I. Implementation Area The study area of ​​this embodiment is the northern segment of the North-South Seismic Belt (32°~43°N, 95°~112°E), covering major active faults such as the northern margin fault of the Altun Mountains, the eastern foothill fault of the Helan Mountains, the Zhuanglanghe Fault, the northern margin fault of the western Qinling Mountains, the eastern Kunlun Fault, the Haiyuan Fault, the Liupanshan Fault, and the northern margin fault of the Qilian Mountains.

[0048] II. Data Sources and Processing Data source: GPS reference station and rover observation data from the China Mainland Tectonic Environment Monitoring Network were selected, including continuous observation data from 1998 to 2018, as well as data from 10 subdivided periods, namely 1999–2001, 2001–2004, …, 2013–2015, and 2015–2020.

[0049] Data processing: The GAMIT / GLOBK software was used for processing. The Eurasian Plate Stability Framework was adopted, with the Ordos Block as the local reference benchmark. Coseismic effects and atmospheric delays were eliminated to obtain velocity field data for each period. The velocity accuracy was better than 0.5 mm / a.

[0050] Profile and baseline extraction: Five profiles covering the Qilian Mountain Fault and Haiyuan Fault were selected, and the velocity projections of stations within 50-350 km on both sides of the faults, perpendicular or parallel to the fault strike, were extracted; the length time series of 18 baselines of 14 reference stations from September 2010 to March 2016 were calculated.

[0051] III. Model Parameter Settings Fault geometry parameters: Based on the geological survey results, the strike, dip, and dip angle of each fault are set. For example, the Zhuanglanghe Fault strikes NNW, dips NE, and has a dip angle of 75°; the North Margin Fault of the West Qinling Mountains strikes NWW, dips NNE, and has a dip angle of 60° to 80°, with a fault depth range of 0 to 30k.

[0052] Negative dislocation model parameters: Each fracture is divided into segments with a length of 5 to 10 km. Each segment is set as a rectangular negative dislocation. The dislocation direction is determined according to the fracture motion type (left-hand strike-slip, thrust compression, etc.).

[0053] Regularization parameter: The regularization parameter λ = 0.01 was determined using the L-curve method to balance data fit and model smoothness.

[0054] IV. Inversion Process and Results Fault slip rate inversion: Through regularized inversion calculation, the slip rate of each fault segment at depths of 0–30k was obtained. Among them, the Zhuanglanghe Fault exhibits left-lateral strike-slip motion with a slip rate of 3.2–4.5 mm / a; the slip rate of the Xiugou-Alakehu segment of the East Kunlun Fault is stable at 4.0–5.2 mm / a.

[0055] Calculation of fault slip deficit: The slip deficit of the Zhuanglanghe Fault at the Yongdeng County-Lanzhou border reaches 2.69 mm / a, the slip deficit of the Luanfeng segment of the northern margin fault of the West Qinling Mountains ranges from 0 to 3.22 mm / a, and the overall slip deficit of the Liupanshan Fault increases to 1.39 to 6.34 mm / a.

[0056] Inversion of fault closure rates: The closure rates of the Subei Kuantanshan segment of the northern margin fault of the Altun Mountains and the southern segment of the eastern foothills fault of the Helan Mountains are higher than 80%; the closure rates of the East Kunlun Fault are evenly distributed in other parts except for Tuosuohu and Maqu, ranging from 70% to 85%; the closure rate of the Hasishan-Machangshan segment of the Haiyuan Fault reaches 82%.

[0057] Dynamic comparison: The sliding loss rate of the Alakhu-Tosohu segment of the East Kunlun Fault increased by 5.12 mm / a from 2015 to 2020 compared with that from 2013 to 2015; the thrust compression and blocking of the Qilian Fault began to be obvious around 2009, and the thrust compression continued to intensify from 2013 to 2015.

[0058] V. Result Validation and Application Verification results: The high-value regions of the locking rate are highly consistent with the regions of enhanced strain rate field, such as the high-value region of the maximum shear strain rate corresponding to the high-locking segment of the Zhuanglanghe Fault; the change in the length of the fracture baseline shows that the baseline shortening rate of the high-locking segment is significant, indicating obvious stress accumulation.

[0059] Earthquake Risk Assessment: Based on the inversion results and seismic activity analysis, 12 areas, including the Subei Kuantanshan segment of the northern margin of the Altyn Tagh fault and the Yongdeng County-Lanzhou border segment of the Zhuanglanghe fault, were identified as high-risk areas for major earthquakes. There is a possibility of an 8-magnitude earthquake occurring in the northern segment of the North-South Seismic Belt in the next 50 years. The Qilian Fault Zone and Haiyuan Fault have priority conditions for the occurrence of 6-7 magnitude earthquakes.

[0060] The above embodiments are only used to illustrate the present invention and are not intended to limit the technical solutions described herein. Although the present invention has been described in detail with reference to the above embodiments, the present invention is not limited to the specific embodiments described above. Therefore, any modifications or equivalent substitutions to the present invention, as well as all technical solutions and improvements that do not depart from the spirit and scope of the invention, are covered within the scope of the claims of the present invention.

Claims

1. A method for inverting fracture locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints, characterized in that, Includes the following steps: Step (1) GNSS observation data acquisition: Collect multi-period observation data from GNSS reference stations and rover stations covering the target area of ​​the China Continental Tectonic Environment Monitoring Network. The target area is the northern segment of the North-South Seismic Belt, which is 32°~43°N and 95°~112°E. Step (2) Data preprocessing: The GNSS observation data is subjected to gross error removal based on spatiotemporal correlation, and high-precision atmospheric delay correction is performed by combining the global atmospheric grid model and regional water vapor radiometer data; the coseismic influence is removed by InSAR coseismic deformation field constraint, and a clean GNSS deformation sequence with both time resolution and accuracy is generated. Step (3) GNSS velocity field calculation: Clean observation data is calculated using the MIT open-source software GAMIT / GLOBK to construct a multi-period GNSS velocity field under a unified reference frame (Eurasian plate stability frame); Step (4) Three-dimensional geometric modeling of the fault: Based on the geological survey results and previous research data, establish a three-dimensional geometric model of the main active faults in the target area, and clarify the fault strike, dip, dip angle and spatial distribution of 0-30km depth; Step (5) Construction of elastic dislocation constraint model: The fault-locked segment is equivalent to a rectangular negative dislocation unit embedded in an elastic half-space medium, and the Green's function response matrix between the fault negative dislocation slip deficit and the surface GNSS deformation rate is established; on this basis, a dual physical constraint mechanism is introduced: the strain rate tensor inside the secondary block calculated based on the block dynamics model is used as a priori constraint to ensure that the inversion results are consistent with the regional tectonic stress field; the historical earthquake recurrence interval, paleoseismic slip and earthquake gap distribution are integrated to construct a probability weight function and give higher inversion priority to high seismic hazard segments; Step (6) Regularized inversion of fracture parameters: The adaptive Tikhonov regularization method is adopted to automatically adjust the regularization factor according to the geometric complexity of the fracture segment, solve the constrained inversion equation, and simultaneously calculate the long-term slip rate, coseismic slip amount and post-earthquake viscoelastic relaxation rate of each fracture segment to obtain the accurate fracture slip loss rate. Step (7) Quantitative calculation of fracture blocking rate: Calculate the blocking rate of each fracture segment based on the ratio of fracture slip loss rate to long-term slip rate, and then calculate the blocking rate based on the fracture slip loss rate obtained from the inversion. Long-term fault slip rate estimated in conjunction with geological / geometry According to the formula Calculate the blocking rate of each segment and introduce an uncertainty propagation model to quantify the confidence interval of the blocking rate; Step (8) Result verification and risk identification: The reliability of the inversion results is confirmed by multi-period GNSS velocity field time series collaborative verification, InSAR deformation rate field cross-verification and historical strong earthquake recurrence cycle matching analysis. By combining the three-dimensional distribution of fracture locking rate and Coulomb stress loading rate, high-risk seismogenic areas for major earthquakes can be identified.

2. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The multi-period observation data mentioned in step (1) includes continuous observation data from 1998 to 2018, and data from 10 subdivided periods: 1999 to 2001, 2001 to 2004, ..., 2013 to 2015, and 2015 to 2020.

3. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The accuracy of the GNSS velocity field obtained in step (3) is better than 0.5 mm / a. At the same time, the time series of 18 baseline lengths and the velocity projection data of 5 profiles across the active fault are extracted.

4. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The main active faults mentioned in step (4) include the northern margin fault of the Altyn Tagh Mountains, the eastern foothills of the Helan Mountains, the Zhuanglanghe Fault, the northern margin fault of the western Qinling Mountains, the eastern Kunlun Mountains, the Haiyuan Fault, the Liupanshan Fault, and the northern margin fault of the Qilian Mountains.

5. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The response relationship described in step (5) is characterized by the elastic dislocation Green's function matrix, and the formula is as follows: ,in This is the velocity field vector observed by GNSS. Let be the Green's function matrix, and d be the sliding parameter vector.

6. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The regularized inversion equation mentioned in step (6) is: ,in The regularization parameter is determined by the L-curve method and has a value of 0.

01. L is the smoothing matrix.

7. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The formula for calculating the locking rate in step (7) is as follows: ,in Where is the fracture lock-up rate, and D is the fracture slip deficit rate (mm / a). The long-term sliding rate of the fracture is (mm / a). The value range is 0 to 100%.

8. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The multi-period result comparison in step (8) includes calculating the change in fracture slip loss rate between the period 2015-2020 and the period 2013-2015. The strain rate field co-verification includes the correlation analysis of the principal strain rate field, the maximum shear strain rate field and the surface strain field.

9. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, In step (4), when performing three-dimensional geometric modeling of the fracture, each fracture is divided into segments with a length of 5 to 10 km. The segment division accuracy matches the distribution density of GNSS stations.

10. The method for inverting fault locking rate based on GNSS crustal deformation observation and elastic dislocation theory constraints as described in claim 1, characterized in that, The high-risk areas for major earthquakes mentioned in step (8) include 12 areas with a closure rate higher than 70%, such as the Subei Kuantanshan section of the northern margin fault of the Altyn Tagh, the Yongdeng County-Lanzhou border section of the Zhuanglanghe fault, and the Tianshui and Luanfeng sections of the northern margin fault of the West Qinling Mountains.

Images