Method for identifying spatial position and blocking strength of blocking concave-convex body of active fault

Abstract

The invention discloses an active fault locking concave-convex body space position and locking intensity identification method, and relates to the geophysical exploration and earthquake monitoring and forecasting technology field, and the method comprises the steps: processing InSAR data to obtain an InSAR earth surface deformation rate field; obtaining a horizontal velocity field of the GNSS observation station based on the data of the GNSS observation station; eliminating aftershock data in the seismic directory data to obtain a seismic fine positioning result; constructing a target fault three-dimensional geometric model; obtaining the spatial distribution of the locking coefficient of the target fault region; determining the position of a concave-convex body in the target fault plane; determining spatial distribution of blocking strength values of the target fault plane; and generating an active fault atresia concavo-convex body distribution and atresia intensity identification graph based on the position of the concavo-convex body and the spatial distribution of the atresia intensity value. According to the invention, the spatial position of the locking concave-convex body and the locking strength identification precision can be improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical exploration and earthquake monitoring and prediction, in particular to a method for identifying the spatial position and locking strength of an active fault locking concave-convex body. BACKGROUND

[0002] The locking concave-convex body on an active fault is a region with high intensity and strong locking degree on the fault plane, and accumulates strain to cause an earthquake when ruptured. Accurate identification of the concave-convex body is crucial for earthquake risk assessment. The current main identification methods include the following three types: 1. Geophysical observation method (1) Seismological imaging and aftershock positioning method: This method deduces the position of the concave-convex body through the spatial clustering of aftershocks after the main shock (aftershocks are often distributed around the boundary of the concave-convex body), and analyzes the stress concentration area combined with the focal mechanism solution. The disadvantages are: the spatial resolution of the identified concave-convex body is low, usually greater than 1 km; the aftershock distribution is often affected by the main shock rupture process, so it cannot fully reflect the complete state of the locking area.

[0003] (2) Geodetic monitoring inversion method: This method inverts the fault locking degree through surface deformation, and the locking area is represented as a high value area of strain accumulation. The disadvantages are: the inversion result is not unique, and the boundary between the locking area and the creep area is ambiguous; the surface deformation has weak response to deep locking, and it is difficult to detect the concave-convex body at depth.

[0004] (3) Geophysical field anomaly method: This method uses the characteristics that the basement property mutation zone (magnetic anomaly zone) is easy to form a stress concentration area, thereby indicating the position of the potential concave-convex body. The disadvantages are: geophysics may be caused by multiple geological factors, and there is no direct mathematical relationship with the locking state; the resolution is low, and it is only suitable for large tectonic plate boundaries.

[0005] 2. Geological and geomorphic analysis method (1) Geological and geomorphic landmark identification method: This method uses field geological and geomorphic survey to determine the position of the concave-convex body by fault scarps, triangular faces, and water system faults. The disadvantages are: geological and geomorphic features are easily modified by erosion, and the characteristics are not obvious; the locking strength cannot be quantified, and only the activity can be qualitatively judged.

[0006] (2) Fault rock and structural feature analysis: This method is based on the principle that high-strength fault rock (such as granite) is more likely to form a concave-convex body, and believes that the structural lens, the folding zone, and the fault gouge thickness reflect the historical sliding behavior. The disadvantages are: the scale of the fault outcrop is very limited, and it is difficult to represent the deep fault state; the ancient structure may be superimposed by later activities.

[0007] (3) Quaternary geology and chronology method: This method determines the latest activity age by optically stimulated luminescence (OSL) and carbon 14 (14C) dating, and identifies the recurrence period by combining paleoseismic trenching. Its shortcomings are: large dating error (OSL error is often > 10%), lack of dating materials for young sediments (< 200 years), and ignoring the dynamic process of stress accumulation in the estimation of recurrence period.

[0008] 3. Experimental simulation and numerical simulation technology (1) Rock friction experiment: This method simulates the stick-slip behavior of faults in the laboratory, and uses acoustic emission to locate the high micro-fracture area (equivalent to asperities). Its shortcomings are: scale effect (laboratory sample < 1 m, natural fault > km) leads to unreliable behavior extrapolation; it is difficult to reproduce the underground temperature and pressure and fluid conditions in the experiment.

[0009] (2) Numerical simulation and computational mechanics method: This method simulates the role of asperities in interseismic locking and coseismic rupture based on the rate-state friction law. Its shortcomings are: friction parameters (a-b values) are difficult to constrain, and the model simplifies the geometric complexity (such as ignoring branch faults).

[0010] (3) Fault roughness modeling: This method obtains the micron-level topography by laser scanning, and considers that low-value areas of roughness correspond to locked asperities. Its shortcomings are: weathering of fault outcrop surface affects the authenticity of the topography; deep fault morphology cannot be directly measured and needs to rely on seismic data for speculation, which is uncertain.

[0011] In summary, the existing technology lacks a comprehensive method that can accurately locate the spatial geometric position of the locked asperities of active faults and simultaneously quantitatively evaluate the locking strength with clear physical meaning. This limits our understanding of the current stress state of the fault and the seismic risk area, and thus affects the scientificity and accuracy of seismic risk assessment. SUMMARY

[0012] The purpose of the present application is to provide a method for identifying the spatial position and locking strength of locked asperities of active faults, which can improve the accuracy of identifying the spatial position and locking strength of locked asperities, can be directly applied to fine seismic risk assessment, and thus increase the scientificity and accuracy of seismic risk assessment.

[0013] To achieve the above purpose, the present application provides the following solutions: In a first aspect, the present application provides a method for identifying the spatial position and locking strength of locked asperities of active faults, comprising: obtaining InSAR data, GNSS observation station data and earthquake catalog data of a target fault area; processing the InSAR data to obtain an InSAR ground deformation rate field; obtaining a GNSS observation station horizontal velocity field based on GNSS observation station data; obtaining aftershock data in the earthquake catalog data, and obtaining earthquake precise positioning results; constructing a target fault three-dimensional geometric model based on the earthquake precise positioning results; obtaining a spatial distribution of a target fault area locking coefficient based on the InSAR ground surface deformation rate field, the GNSS observation station horizontal velocity field and the target fault three-dimensional geometric model; determining a position of a concave-convex body in a target fault plane based on the spatial distribution of the locking coefficient and the earthquake precise positioning results; determining a spatial distribution of a locking strength value of the target fault plane based on the spatial distribution of the locking coefficient, the target fault three-dimensional geometric model, the InSAR ground surface deformation rate field and the GNSS observation station horizontal velocity field; generating an active fault locking concave-convex body distribution and a locking strength identification map based on the position of the concave-convex body and the spatial distribution of the locking strength value.

[0014] In an embodiment, the InSAR data is InSAR data of multiple periods, multiple orbits and a set spatial resolution; processing the InSAR data to obtain an InSAR ground surface deformation rate field, comprising: using ISCE software to sequentially perform interference pair combination selection, precise orbit correction, phase unwrapping, atmospheric delay correction, terrain phase removal and geocoding processing on the InSAR data of multiple periods to generate a ground surface deformation time series; using StaMPS software to process the ground surface deformation time series to obtain the InSAR ground surface deformation rate field.

[0015] In an embodiment, the GNSS observation station horizontal velocity field is obtained based on GNSS observation station data, comprising: obtaining a coordinate time series of each observation station in an international terrestrial reference frame based on GNSS observation station data; obtaining a GNSS observation station horizontal velocity field based on the coordinate time series of each observation station in the international terrestrial reference frame.

[0016] In an embodiment, the coordinate time series of each observation station in the international terrestrial reference frame is obtained based on the GNSS observation station data, comprising: using GAMIT / GLOBK software to process the GNSS observation station data to obtain the coordinate time series of each observation station in the international terrestrial reference frame.

[0017] In an embodiment, the GNSS observation station horizontal velocity field is obtained based on the coordinate time series of each observation station in the international terrestrial reference frame, comprising: screening the coordinate time series of each observation station under the international terrestrial reference frame based on the set error, to obtain a screened coordinate time series; fitting the screened coordinate time series by using a function containing tectonic movement and periodic non-tectonic movement, to obtain a horizontal velocity field of the GNSS observation station.

[0018] In an embodiment, a three-dimensional geometric model of a target fault is constructed based on the earthquake precise positioning result, including: acquiring a historical active fault geometric structure; determining a target fault geometric structure based on the earthquake precise positioning result and the historical active fault geometric structure; the target fault geometric structure includes a strike, a dip angle, segmentation information, an upper and lower wall boundary, and a deep extension range of a target fault surface; based on the target fault geometric structure, discretizing the target fault surface to obtain the three-dimensional geometric model of the target fault.

[0019] In an embodiment, a spatial distribution of a target fault region locking coefficient is obtained based on the InSAR ground surface deformation rate field, the GNSS observation station horizontal velocity field, and the three-dimensional geometric model of the target fault, including: taking the InSAR ground surface deformation rate field and the GNSS observation station horizontal velocity field as observation values; constructing a relationship model of ground surface deformation and fault surface locking coefficient based on the three-dimensional geometric model of the target fault by using a negative dislocation model; inverting the relationship model of ground surface deformation and fault surface locking coefficient based on the observation values by using a constrained least squares inversion method, to obtain the spatial distribution of the target fault region locking coefficient.

[0020] In an embodiment, a position of a target fault surface concave-convex body is determined based on the spatial distribution of the target fault region locking coefficient and the earthquake precise positioning result, including: based on the earthquake precise positioning result, acquiring a spatial distribution of earthquakes reaching a set magnitude in a set time range in a neighboring region within a set range of a target fault; determining a potential concave-convex body candidate region in the spatial distribution of earthquakes based on a set condition; performing b value spatial scanning in a target fault region to determine a potential target fault surface high locking region; determining the position of the target fault surface concave-convex body based on the potential concave-convex body candidate region, the potential target fault surface high locking region, and the spatial distribution of the target fault region locking coefficient.

[0021] In an embodiment, the spatial distribution of the locking strength value of the target fault plane is determined based on the spatial distribution of the target fault region locking coefficient, the three-dimensional geometric model of the target fault, the InSAR ground surface deformation rate field and the GNSS observation station horizontal velocity field, and comprises: the spatial distribution of the target fault region locking coefficient is taken as a fault constraint condition; a three-dimensional viscoelastic half-space finite element model is constructed based on the three-dimensional geometric model of the target fault; a ground surface deformation rate boundary condition consistent with the relative motion of the target fault region block is determined based on the InSAR ground surface deformation rate field and the GNSS observation station horizontal velocity field; the three-dimensional viscoelastic half-space finite element model is subjected to quasi-static simulation based on the fault constraint condition and the ground surface deformation rate boundary condition, so as to obtain the spatial distribution of the locking strength value of the target fault plane.

[0022] In an embodiment, the active fault locking relief body distribution and locking strength identification map are generated based on the position of the target fault plane relief body and the spatial distribution of the locking strength value of the target fault plane, and comprise: a locking relief body of the target fault region is determined based on the position of the target fault plane relief body and the spatial distribution of the locking strength value of the target fault plane by adopting a set identification condition, and information data of the locking relief body is recorded; the information data of the locking relief body comprises: spatial position, area, average locking strength value and maximum locking strength value; the locking strength grade of the locking relief body is determined based on the average locking strength value of the locking relief body and a set strength grade condition; the active fault locking relief body distribution and locking strength identification map are generated based on the spatial distribution of the locking strength value of the target fault plane, the information data of the locking relief body and the locking strength grade of the locking relief body.

[0023] According to the specific embodiments provided in the present application, the present application has the following technical effects: The application provides a method for identifying the spatial position and locking strength of an active fault locking concave-convex body. BRIEF DESCRIPTION OF DRAWINGS

[0024] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without any creative effort on the basis of these drawings.

[0025] Figure 1 A flow chart of a method for identifying the spatial position and locking strength of an active fault locking concave-convex body according to an embodiment of the present application is shown in FIG. 1. Figure 2 A schematic diagram of the overall process of a method for identifying the spatial position and locking strength of an active fault locking concave-convex body according to an embodiment of the present application is shown in FIG. 2. Figure 3 A schematic diagram of an InSAR ground deformation rate field according to an embodiment of the present application is shown in FIG. 3. Figure 4 A schematic diagram of a GNSS observation station horizontal velocity field according to an embodiment of the present application is shown in FIG. 4. Figure 5 A schematic diagram of a target fault three-dimensional geometric model according to an embodiment of the present application is shown in FIG. 5. Figure 6 A schematic diagram of the spatial distribution of the locking coefficient of a target fault region according to an embodiment of the present application is shown in FIG. 6. Figure 7 A schematic diagram of the analysis result of the seismic activity of a target fault region according to an embodiment of the present application is shown in FIG. 7. Figure 8 A schematic diagram of a three-dimensional viscoelastic half-space finite element model provided in one embodiment of the present application; Figure 9 A schematic diagram of a target fault closure asperity and its closure strength provided in one embodiment of the present application; Figure 10 A schematic diagram of the structure of a computer device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0026] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0027] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application is further described in detail below with reference to the accompanying drawings and specific implementation methods.

[0028] In an exemplary embodiment, Figure 1 As shown, a method for identifying the spatial position and locking strength of a locked concave-convex body of an active fault is provided, comprising: Step 100: Acquire InSAR data, GNSS observation station data, and earthquake catalog data of the target fault area.

[0029] Step 200: Process the InSAR data to obtain the InSAR surface deformation rate field. The GNSS station horizontal velocity field is obtained based on the GNSS station data. Aftershock data from the earthquake catalog is removed to obtain the earthquake precise location result.

[0030] Step 300: construct a three-dimensional geometric model of the target fault based on the earthquake precise positioning results.

[0031] Step 400: obtaining the spatial distribution of the locking coefficient of the target fault region based on the InSAR surface deformation rate field, the GNSS observation station horizontal velocity field, and the three-dimensional geometric model of the target fault.

[0032] Step 500: Determine the position of the asperity in the target fault plane based on the spatial distribution of the locking coefficient and the seismic precise positioning result.

[0033] Step 600: Determine the spatial distribution of the blocking strength value of the target fault plane based on the spatial distribution of the blocking coefficient, the three-dimensional geometric model of the target fault, the InSAR surface deformation rate field, and the horizontal velocity field of the GNSS observation station.

[0034] Step 700, generating the active fault locking asperity distribution and locking strength identification map based on the position of the asperity and the spatial distribution of the locking strength value.

[0035] By implementing steps 100-700, the active fault locking asperity distribution and locking strength identification map can be generated, which intuitively reflects the spatial position of the locking asperity and the locking strength, and improves the accuracy of identifying the spatial position of the locking asperity and the locking strength. The generated active fault locking asperity distribution and locking strength identification map can be directly applied to aspects such as fine seismic risk assessment, thereby increasing the scientificity and accuracy of seismic risk assessment.

[0036] As an optional implementation, in order to improve the accuracy of identifying the spatial position and boundary of the locking asperity and overcome the problem of insufficient resolution of a single data source, the InSAR data is InSAR data of multiple periods, multiple orbits, and a set spatial resolution. Based on this, the process of obtaining the InSAR ground surface deformation rate field in step 200 can include: using ISCE software to sequentially perform interference pair combination selection, precise orbit correction, phase unwrapping, atmospheric delay correction, terrain phase removal, and geographic coding processing on the InSAR data of multiple periods, to generate a ground surface deformation time series. Using StaMPS software to process the ground surface deformation time series to obtain the InSAR ground surface deformation rate field. Based on this, the GNSS observation station horizontal velocity field is obtained based on the GNSS observation station data in step 200, which specifically includes: Step 210, obtaining the coordinate time series of each observation station in the International Terrestrial Reference Frame based on the GNSS observation station data. Step 210 includes: using GAMIT / GLOBK software to process the GNSS observation station data to obtain the coordinate time series of each observation station in the International Terrestrial Reference Frame.

[0037] Step 220, obtaining the GNSS observation station horizontal velocity field based on the coordinate time series of each observation station in the International Terrestrial Reference Frame. Step 220 includes: based on a set error, filtering the coordinate time series of each observation station in the International Terrestrial Reference Frame to obtain a filtered coordinate time series. Using a function containing tectonic movement and periodic non-tectonic movement to fit the filtered coordinate time series to obtain the GNSS observation station horizontal velocity field.

[0038] For example, the InSAR (Interferometric Synthetic Aperture Radar) data of a target fault region is obtained, which is a synthetic aperture radar interferometry data, and the InSAR data is obtained by acquiring a plurality of periods (200 periods in the embodiment, each period is 12 days, and the time interval is 12 days), a plurality of tracks (5 tracks in the embodiment) and a set spatial resolution (the spatial resolution is less than 100 m in the embodiment) of the target fault region. The InSAR data of the plurality of periods is sequentially subjected to interference pair combination selection, precise orbit correction, phase unwrapping, atmospheric delay correction, terrain phase removal and geocoding processing by using ISCE (InSAR Scientific Computing Environment) software, so as to generate a ground surface deformation time sequence. The ground surface deformation time sequence is processed by using StaMPS (Stanford Method for Persistent Scatterers) software, so as to obtain an InSAR ground surface deformation rate field.

[0039] GNSS (Global Navigation Satellite System) observation station data of a target fault region is obtained, that is, GNSS observation station data. The GNSS observation station data is processed by using GAMIT / GLOBK software, so as to obtain a coordinate time sequence of each observation station in the ITRF (International Terrestrial Reference Frame). Data with an observation error greater than three times the mean error (that is, a set error) in the coordinate time sequence of each observation station in the ITRF is removed, so as to obtain a screened coordinate time sequence. The screened coordinate time sequence is fitted by using a function containing tectonic movement and periodic non-tectonic movement, so as to obtain a GNSS observation station horizontal velocity field. The GAMIT / GLOBK software is high-precision GNSS data processing software developed by MIT (Massachusetts Institute of Technology) and SIO (Scripps Institution of Oceanography).

[0040] Seismic catalog data of a target fault region is obtained, which contains information such as time, latitude and longitude, depth, magnitude and focal mechanism of earthquake occurrence. The aftershock data in the seismic catalog data is removed, so as to obtain a seismic precise positioning result.

[0041] As an optional implementation, in order to provide a basis for obtaining the spatial distribution of the locking coefficient, step 300 comprises: obtaining historical active fault geometry. The target fault geometry is determined based on the earthquake precise positioning result and the historical active fault geometry. The target fault geometry comprises the strike, dip angle, segmentation information, upper and lower wall boundary, and deep extension range of the target fault plane. Based on the target fault geometry, the target fault plane is discretized to obtain a three-dimensional geometric model of the target fault.

[0042] For example, the historical active fault geometry is obtained. The three-dimensional geometry of the target fault (i.e., the target fault geometry) is determined based on the earthquake precise positioning result and the historical active fault geometry. The target fault geometry comprises the strike, tendency, dip angle, segmentation information, upper and lower wall boundary, and deep extension range of the target fault plane. According to the obtained three-dimensional geometry of the target fault, the target fault plane is discretized into regular rectangular (or triangular or other shapes) grid cells to form a three-dimensional geometric model of the target fault. The grid size is set according to the actual research requirements and resolution requirements.

[0043] As an optional implementation, in order to improve the accuracy of the spatial distribution of the locking coefficient of the target fault region, step 400 comprises: taking the InSAR ground deformation rate field and the GNSS observation station horizontal velocity field as observation values. A relationship model between ground deformation and fault plane locking coefficient is constructed based on the three-dimensional geometric model of the target fault using a negative dislocation model. The relationship model between ground deformation and fault plane locking coefficient is inverted based on the observation values using a constrained least squares inversion method to obtain the spatial distribution of the locking coefficient of the target fault region.

[0044] For example, the InSAR ground deformation rate field and the GNSS observation station horizontal velocity field obtained in step 200 are taken as observation values. According to the three-dimensional geometric model of the target fault constructed in step 300, a relationship model between ground deformation and fault plane locking coefficient is constructed using a negative dislocation model. The relationship model between ground deformation and fault plane locking coefficient is inverted based on the observation values using a constrained least squares inversion method to invert the locking coefficient on each grid cell of the target fault plane. Laplace smoothing and optimization of the corresponding regularization parameters are considered in the inversion process, so that the predicted deformation and the observation values can be optimally fitted, and the rationality of the relationship model between ground deformation and fault plane locking coefficient is ensured, thereby obtaining the spatial distribution of the locking coefficient of the target fault region.

[0045] As an optional implementation, in order to improve the accuracy of the position identification of the target fault plane concave-convex body, step 500 comprises: based on the seismic fine positioning result, obtaining the spatial distribution of earthquakes reaching a set magnitude in a set time range in the adjacent area within the set range of the target fault. Determine the potential concave-convex body candidate area in the spatial distribution based on the set condition. Perform b value spatial scanning in the target fault area to determine the potential target fault plane high locking area. Determine the position of the target fault plane concave-convex body based on the spatial distribution of the potential concave-convex body candidate area, the potential target fault plane high locking area, and the locking coefficient of the target fault area.

[0046] For example, according to the seismic fine positioning result obtained in step 200, the spatial distribution of M≥3.0 earthquakes in the adjacent area (for example, within 20 km on both sides of the fault) within the set range of the target fault in the past 30 years (that is, the set time range) is counted. M represents the magnitude of the earthquake. In the spatial distribution of earthquakes, the earthquake empty area meeting the conditions of an area greater than 50 km 2 , the duration of no earthquake (specifically, no M≥3.0 earthquake) is greater than 3 years, and the surrounding is surrounded by significant seismic activity (that is, the set condition) is identified. Mark the position of these earthquake empty areas as potential concave-convex body candidate areas.

[0047] Divide the grid cells in the target fault area (the grid size can be 0.1x0.1), and calculate the b value of M≥3.0 earthquakes in the target fault area and each grid cell in the target fault area based on the maximum likelihood method. The b value of the target fault area is defined as the background b value, and the b value map of each grid cell is drawn based thereon, and the low b value area of each grid cell significantly lower than the background b value of the target fault area (that is, the area with a b value less than the background b value minus three times the standard error of the background b value) is identified. The low b value area is used as the potential target fault plane high locking area.

[0048] The position of the area meeting the conditions of earthquake empty area, low b value area, and locking coefficient of the area being large (that is, the locking coefficient is greater than 0.7) is marked as the concave-convex body, and the position of the target fault plane concave-convex body is identified.

[0049] As an optional implementation, in order to quantitatively calculate the locking strength value of the target fault plane, step 600 comprises: taking the spatial distribution of the locking coefficient of the target fault area as the fault constraint condition. Based on the three-dimensional geometric model of the target fault, a three-dimensional viscoelastic half-space finite element model is constructed. Based on the InSAR ground deformation rate field and the GNSS observation station horizontal velocity field, the ground deformation rate boundary condition consistent with the relative motion of the target fault area block is determined. Based on the fault constraint condition and the ground deformation rate boundary condition, the three-dimensional viscoelastic half-space finite element model is simulated in a quasi-static state to obtain the spatial distribution of the locking strength value of the target fault plane.

[0050] For example, the locking strength of a grid cell on the target fault plane can be defined as the rate of shear stress accumulation (unit: MPa / yr) at this point caused by fault locking.

[0051] The spatial distribution of the target fault zone locking coefficient obtained in step 400 is taken as the fault constraint condition. According to the three-dimensional geometric model of the target fault constructed in step 300, the rock mechanics parameters (such as Young's modulus, Poisson's ratio) are set, and then a three-dimensional viscoelastic half-space finite element model is constructed. According to the InSAR ground deformation rate field and the GNSS observation station horizontal velocity field obtained in step 200, the boundary condition of the surface deformation rate consistent with the relative motion of the target fault zone block (a geological unit that behaves as a relatively rigid whole motion in tectonic movement) is determined, that is, the boundary condition of the three-dimensional viscoelastic half-space finite element model.

[0052] On the target fault plane, for a grid cell with a locking coefficient of , its slip rate is constrained to be , where is the long-term slip rate of the grid cell (which can be given by the block relative motion rate determined by geodetic surveying).

[0053] According to the fault constraint condition and the ground deformation rate boundary condition, the three-dimensional viscoelastic half-space finite element model is subjected to quasi-static simulation, and the change of shear stress accumulated on each grid cell of the target fault plane with time is calculated under the condition of applying the fault constraint condition and the ground deformation rate boundary condition. The change value of shear stress on each grid cell of the target fault plane after the quasi-static simulation reaches a steady state (for example, after 1 year of simulation) is . The locking strength of each grid cell is represented by the rate of shear stress accumulation at this point, and the specific formula is: . Wherein, is the time to reach a steady state in the quasi-static simulation. The locking strength value (unit: MPa / yr) of each grid cell on the target fault plane is calculated by this formula, and the spatial distribution of the locking strength value of the target fault plane is obtained.

[0054] As an optional embodiment, in order to visually represent the spatial distribution of the position of the asperity and the locking strength value, the step 700 comprises: adopting a set identification condition, determining the locking asperity of the target fault region based on the position of the asperity of the target fault surface and the spatial distribution of the locking strength value of the target fault surface, and recording the information data of the locking asperity. The information data of the locking asperity comprises: spatial position, area, average locking strength value and maximum locking strength value. The locking strength grade of the locking asperity is determined based on the average locking strength value of the locking asperity and the set strength grade condition. The active fault locking asperity distribution and the locking strength identification map are generated based on the spatial distribution of the locking strength value of the target fault surface, the information data of the locking asperity and the locking strength grade of the locking asperity.

[0055] For example, the asperity identification is performed according to the position of the asperity in the target fault surface obtained in the step 500 and the spatial distribution of the locking strength value of the target fault surface obtained in the step 600, which comprises the following steps: (1) Set a locking strength threshold S_threshold (for example, S_threshold = 0.04 MPa / yr) and a minimum area threshold A_min (for example, A_min = 50 km 2 ) as the set identification condition.

[0056] (2) According to the spatial distribution of the locking strength value of the target fault surface, find all the regions (i.e. continuous regions) on the target fault surface which are S≥S_threshold and spatially continuous. Wherein S represents the average locking strength value.

[0057] (3) For the continuous regions in (2), determine whether the area is greater than or equal to A_min. Identify each continuous region with an area ≥ A_min as an independent locking asperity. Record the information data of the identified locking asperity. The information data of the locking asperity comprises: spatial position (including the boundary range of the center coordinates), area, average locking strength, maximum locking strength.

[0058] (4) For the identified locking asperity, perform strength grading according to the average locking strength value to determine the locking strength grade of the locking asperity. In this embodiment, the locking strength grade is classified as: ① high locking strength: S≥0.1 MPa / yr; ② medium locking strength: 0.05 MPa / yr≤S<0.1 MPa / yr; ③ low locking strength: S<0.05 MPa / yr (but still greater than S_threshold).

[0059] (5) Generate an active fault locking asperity distribution and locking strength identification map based on the spatial distribution of the locking strength values ​​of the target fault plane, the information data of the locking asperities, and the locking strength levels of the locking asperities. The active fault locking asperity distribution and locking strength identification map specifically includes: ① the background is the spatial distribution of the locking strength values ​​obtained in step 600 on the target fault plane; ② the boundary contours of each identified locking asperity are clearly marked; ③ the strength level of each locking asperity is marked within or next to the contour of the asperity.

[0060] In an exemplary embodiment, in combination with the steps in the above embodiments, the Laohushan-Haiyuan Fault Zone on the eastern edge of the Qinghai-Tibet Plateau is used as an example for explanation. Figure 2 shown.

[0061] S1, data acquisition and data processing.

[0062] InSAR data: Ascending and descending orbit data from the European Space Agency's Sentinel-1 A / B satellites covering the region from 2015 to 2022 were acquired. ISCE software was used to perform interferometric pair selection (maximum baseline <150 m), precise orbit correction, phase unwrapping, atmospheric delay correction, terrain phase removal, and geocoding on the ascending and descending orbit data to generate a surface deformation time series. StaMPS software was used to process the surface deformation time series to generate the average LOS InSAR surface deformation rate field, as shown in the following example: Figure 3 As shown ( Figure 3 Part (a) is the surface deformation rate field of the orbit-raising data. Figure 3 Part (b) in the figure shows the surface deformation rate field of the descending orbit data).

[0063] GNSS observation station data: The original observation data of the China Crustal Movement Observation Network and regional encrypted stations during this period were obtained. GAMIT / GLOBK (version 7.01) was used for processing, and the reference frame was ITRF2014. After time series fitting, the horizontal velocity field of the GNSS observation station was extracted, as shown in the following example: Figure 4 shown.

[0064] Earthquake catalog data: The precise location catalog of earthquakes with magnitude ≥ 3 from 2000 to 2023 provided by the China Earthquake Networks Center is used.

[0065] S2, construct a three-dimensional geometric model of the target fault.

[0066] Based on the results of geological surveys and deep exploration, the geometric structure of the fault area is determined, and the fault surface (i.e., the target fault surface) is discretized into a rectangular grid of 2 km (along the strike) × 2 km (along the dip). The three-dimensional geometric model of the fault surface (i.e., the three-dimensional geometric model of the target fault) is obtained, as shown in the figure. Figure 5 shown.

[0067] S3, spatial distribution inversion of target fault zone closure coefficient.

[0068] The average LOS to InSAR ground deformation rate field and GNSS observation station horizontal velocity field (N sites) are taken as observation values.

[0069] The relationship model between surface deformation and fault plane closure coefficient is constructed based on the target fault three-dimensional geometric model using the elastic half-space dislocation model. DEFNODE (negative dislocation model inversion program) is used for slip behavior inversion. The objective function is the least squares difference between the observation value and the predicted value of the relationship model between surface deformation and fault plane closure coefficient, with Laplace second-order smoothing regularization. The optimal regularization parameter is determined by the L curve method. The spatial distribution of the preliminary closure coefficient (i.e. the spatial distribution of the target fault zone closure coefficient) is obtained, i.e. the spatial distribution of the closure coefficient of the LHSF, HYF, LPSF and GGBJF regions, as shown in Figure 6 .

[0070] S4, concave-convex body position determination fused with seismic activity information.

[0071] Seismic gap identification: a significant M≥3.0 seismic gap (area ~50km 2 , duration >3 years) is identified in the middle of the fault. Waveform cross-correlation identifies several clustered repeating earthquakes, mainly distributed at both ends of the fault, and few in the gap, as shown in Figure 7 part (a).

[0072] b-value spatial scanning: b-value calculation (grid size 0.1x0.1, minimum earthquake number 50) finds a significant low b-value zone (b-value 0.6, background b-value ~1.0) under the above seismic gap, as shown in Figure 7 part (b).

[0073] Information fusion and position constraint: fusion method: taking the identified seismic gap / low b-value zone as the center, the closure coefficient is modified to above 0.95 in a circular area with a radius of 7km, and the edge is Gaussian smoothed. The modified closure distribution shows that the closure range in this area is more concentrated and the closure degree is higher.

[0074] S5, quantitative calculation of closure strength.

[0075] A three-dimensional viscoelastic half-space finite element model (size 800 800 km³) is established using Pylith. The Young's modulus is set to 30 GPa and the Poisson's ratio is 0.25. As shown in Figure 8 . The regional background strain rate field is given by the geodetic velocity field as the displacement boundary condition.

[0076] The modified locking coefficient distribution is mapped to the grid cells in the fault region, and the slip constraint of each grid cell is calculated .

[0077] The shear stress accumulation under the time scale of 1 year is simulated. The shear stress increment of each grid cell on the fault surface is extracted . The locking strength value , unit: MPa / yr, the spatial distribution of the locking strength value is obtained, which is displayed in the center of the modified high locking region (potential asperity).

[0078] S6, asperity identification and strength grading map.

[0079] Set S_threshold=0.04MPa / yr, A_min=40km². A main continuous high S region (S>0.04MPa / yr) is identified, with an area of about 6 km² and an average S=0.08MPa / yr. It is identified as a locking asperity. Strength grading: average S=0.08MPa / yr (<0.1MPa / yr), medium locking strength.

[0080] Finally, the locking asperity distribution and locking strength identification map is generated: the background is the S distribution (color bar), which clearly outlines the boundary of the locking asperity (black closed curve), and marks "M" inside it, as shown in Figure 9 .

[0081] In combination with the steps in the above embodiments, the locking asperity spatial position and locking strength identification method provided by the application has the following advantages: 1. High-precision spatial positioning: The geodetic deformation data (providing overall locking pattern constraints) and seismic activity analysis (providing local high-resolution stress state indicators) are combined, which significantly improves the accuracy of identifying the spatial position and boundary of the locking asperity, and overcomes the problem of insufficient resolution of a single data source.

[0082] 2. Locking strength quantification: The direct and quantitative calculation of fault locking strength (shear stress accumulation rate) based on the physical and mechanical model (finite element simulation) is clearly proposed and implemented. The final strength value has a clear physical meaning and can be directly used to evaluate the ability and rate of fault segment to accumulate strain energy, which is a core quantitative index for evaluating seismic potential and risk.

[0083] 3. Clear physical mechanism: The entire technical solution is based on dislocation theory and solid mechanics, and the physical mechanism from data inversion to strength calculation process is clear, and the result is highly reliable.

[0084] 4. Multi-source data collaboration: Make full use of the advantages of modern intensive geodetic network (GPS, InSAR) and seismic network observation, realize the effective complement and collaborative constraint of multi-source information.

[0085] 5. High practicability: The final output contains the distribution of active fault locking concave-convex body and the locking strength identification map of the concave-convex body with accurate position and quantitative intensity level, which is intuitive in form and can be directly used for fine seismic risk assessment, major engineering site safety evaluation and earthquake monitoring and prediction practice, and has important application value.

[0086] 6. Strong generalizability: The above method is clear and can be applied to the research of major active faults (strike-slip, thrust, normal fault) and subduction zones under different tectonic backgrounds around the world.

[0087] In an exemplary embodiment, a computer device is provided, which can be a server or a terminal, and its internal structure diagram can be as shown in Figure 10 The computer device includes a processor, a memory, an input / output interface (I / O) and a communication interface. The processor, the memory and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used to store the spatial position of the active fault locking concave-convex body and the related data of the locking strength identification method. The input / output interface of the computer device is used to exchange information between the processor and external devices. The communication interface of the computer device is used to communicate with external terminals through network connection. The computer program is executed by the processor to implement an active fault locking concave-convex body spatial position and locking strength identification method.

[0088] Those skilled in the art can understand that Figure 10 the structure shown in the figure is only a block diagram of part of the structure related to the scheme of the present application, and does not constitute a limitation on the computer device to which the scheme of the present application is applied. The specific computer device can include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components. In an exemplary embodiment, a computer device is provided, including a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps in the above method embodiments.

[0089] In an exemplary embodiment, a computer readable storage medium storing a computer program is provided, the computer program, when executed by a processor, implements the steps of any of the above method embodiments.

[0090] In an exemplary embodiment, a computer program product is provided, comprising a computer program which, when executed by a processor, implements the steps of any of the above method embodiments.

[0091] It should be noted that the user information (including but not limited to user equipment information, user personal information, etc.) and data (including but not limited to data for analysis, stored data, displayed data, etc.) involved in the present application are all information and data authorized by the user or authorized by all parties, and the collection, use and processing of related data need to comply with relevant regulations.

[0092] It can be understood by those skilled in the art that all or part of the processes in the above embodiments can be completed by a computer program instructing related hardware, and the computer program can be stored in a non-volatile computer readable storage medium. When the computer program is executed, it can include the processes of the above embodiments. Any reference to memory, database or other medium used in the embodiments provided by the present application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical storage, high-density embedded non-volatile memory, resistive memory (ReRAM), magnetoresistive random access memory (MRAM), ferroelectric memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. As an illustration but not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0093] The database involved in each of the embodiments provided in the present application can include at least one of a relational database and a non-relational database. The non-relational database can include a distributed database based on a blockchain, and the like, without being limited thereto. The processor involved in each of the embodiments provided in the present application can be a general-purpose processor, a central processing unit, a graphics processing unit, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, and the like, without being limited thereto.

[0094] The technical features of the above embodiments can be combined in any manner. To make the description concise, all possible combinations of the technical features in the above embodiments are not described, but it should be considered that any combination of the technical features is within the scope of the present disclosure, as long as there is no contradiction.

[0095] The principles and implementation modes of the present application are described by applying specific examples herein, and the above embodiments are only used to help understand the method of the present application and its core idea. Meanwhile, for those skilled in the art, the specific implementation modes and application ranges will be changed according to the idea of the present application. In summary, the content of the present description should not be understood as a limitation of the present application.

Claims

1. A method for identifying the spatial position and locking strength of an active fault locking concave-convex body, characterized in that: include: Obtain InSAR data, GNSS observation station data, and earthquake catalog data for the target fault area; The InSAR data are processed to obtain the InSAR surface deformation rate field; The horizontal velocity field of the GNSS observation station is obtained based on the GNSS observation station data; Eliminate aftershock data from earthquake catalog data to obtain precise earthquake location results; Construct a 3D geometric model of the target fault based on the earthquake precise positioning results; The spatial distribution of the locking coefficient in the target fault area is obtained based on the InSAR surface deformation rate field, the horizontal velocity field of the GNSS observation station and the three-dimensional geometric model of the target fault. Determine the location of the asperity in the target fault plane based on the spatial distribution of the locking coefficient and the seismic precise positioning results; The spatial distribution of the locking intensity value of the target fault plane is determined based on the spatial distribution of the locking coefficient, the three-dimensional geometric model of the target fault, the InSAR surface deformation rate field and the horizontal velocity field of the GNSS observation station; Based on the spatial distribution of the positions of the asperities and the locking strength values, the active fault locking asperity distribution and locking strength identification map are generated.

2. The method for identifying the spatial position and locking strength of an active fault locking concave-convex body according to claim 1, characterized in that: The InSAR data is InSAR data of multiple periods, multiple orbits and set spatial resolution; The InSAR data is processed to obtain an InSAR surface deformation rate field, including: Using ISCE software, the InSAR data of multiple periods are processed in sequence through interferometric pair selection, precise orbit correction, phase unwrapping, atmospheric delay correction, terrain phase removal and geocoding to generate surface deformation time series. The surface deformation time series is processed using StaMPS software to obtain the InSAR surface deformation rate field.

3. The method for identifying the spatial position and locking strength of an active fault locking concavo-convex body according to claim 1, characterized in that: The horizontal velocity field of the GNSS observation station is obtained based on the GNSS observation station data, including: Based on the GNSS observation station data, the coordinate time series of each observation station in the international terrestrial reference frame is obtained; The horizontal velocity field of the GNSS observation station is obtained based on the coordinate time series of each observation station in the international terrestrial reference frame.

4. The method for identifying the spatial position and locking strength of an active fault locking concave-convex body according to claim 3, characterized in that: Based on the GNSS observation station data, the coordinate time series of each observation station in the international terrestrial reference frame is obtained, including: The GNSS observation station data were processed using GAMIT / GLOBK software to obtain the coordinate time series of each observation station in the international terrestrial reference frame.

5. The method for identifying the spatial position and locking strength of an active fault locking asperity according to claim 3, characterized in that: The horizontal velocity field of the GNSS observation station is obtained based on the coordinate time series of each observation station in the international terrestrial reference frame, including: Based on the set error, the coordinate time series of each observation station in the international terrestrial reference frame is filtered to obtain the filtered coordinate time series; The filtered coordinate time series is fitted using a function including tectonic motion and periodic non-tectonic motion to obtain the horizontal velocity field of the GNSS observation station.

6. The method for identifying the spatial position and locking strength of an active fault locking asperity according to claim 1, characterized in that: Constructing a three-dimensional geometric model of the target fault based on the earthquake precise positioning results, including: Obtain historical active fault geometry; Determining a target fault geometry based on the earthquake precise positioning result and the historical active fault geometry; the target fault geometry includes the strike, dip, segmentation information, upper and lower wall boundaries, and deep extension of the fault plane; Based on the target fault geometric structure, the target fault plane is discretized to obtain the three-dimensional geometric model of the target fault.

7. The method for identifying the spatial position and locking strength of an active fault locking concave-convex body according to claim 1, characterized in that: The spatial distribution of the locking coefficient of the target fault region is obtained based on the InSAR surface deformation rate field, the GNSS observation station horizontal velocity field, and the target fault three-dimensional geometric model, including: Taking the InSAR surface deformation rate field and the GNSS observation station horizontal velocity field as observation values; A negative dislocation model is used to construct a relationship model between surface deformation and fault plane locking coefficient based on the three-dimensional geometric model of the target fault; The constrained least squares inversion method is used to invert the relationship model between the surface deformation and the fault plane locking coefficient based on the observation values ​​to obtain the spatial distribution of the locking coefficient of the target fault area.

8. The method for identifying the spatial position and locking strength of an active fault locking asperity according to claim 1, characterized in that: Determining the position of the target fault surface asperity based on the spatial distribution of the locking coefficient of the target fault region and the earthquake precise positioning result includes: Based on the earthquake precise positioning results, obtaining the spatial distribution of earthquakes reaching a set magnitude within a set time range in an adjacent area within a set range of the target fault; Determining a potential asperity candidate area in the earthquake spatial distribution based on set conditions; Perform b-value spatial scanning in the target fault area to determine the high closure area of ​​the potential target fault surface; The position of the target fault plane asperity is determined based on the spatial distribution of the potential asperity candidate area, the potential target fault plane high closure area and the target fault region closure coefficient.

9. The method for identifying the spatial position and locking strength of an active fault locking asperity according to claim 1, characterized in that: Determining the spatial distribution of the blocking strength value of the target fault plane based on the spatial distribution of the blocking coefficient of the target fault region, the three-dimensional geometric model of the target fault, the InSAR surface deformation rate field, and the horizontal velocity field of the GNSS observation station includes: Taking the spatial distribution of the locking coefficient of the target fault area as a fault constraint condition; Constructing a three-dimensional viscoelastic half-space finite element model based on the three-dimensional geometric model of the target fault; Determining a surface deformation rate boundary condition consistent with the relative motion of blocks in the target fault region based on the InSAR surface deformation rate field and the GNSS observation station horizontal velocity field; Based on the fault constraint condition and the surface deformation rate boundary condition, a quasi-static simulation is performed on the three-dimensional viscoelastic half-space finite element model to obtain the spatial distribution of the locking strength value of the target fault plane.

10. The method for identifying the spatial position and locking strength of an active fault locking asperity according to claim 1, characterized in that: Generating an active fault locked concavo-convex body distribution and a locked strength identification map based on the position of the target fault plane concavo-convex body and the spatial distribution of the locked strength value of the target fault plane includes: Using set identification conditions, determining the closed concavo-convex body of the target fault region based on the position of the concavo-convex body of the target fault plane and the spatial distribution of the closed strength value of the target fault plane, and recording information data of the closed concavo-convex body; the information data of the closed concavo-convex body includes: spatial position, area, average closed strength value and maximum closed strength value; determining a locking strength level of the locking concavo-convex body based on an average locking strength value of the locking concavo-convex body and a set strength level condition; The active fault locking asperity distribution and locking strength identification map is generated based on the spatial distribution of the locking strength value of the target fault plane, the information data of the locking asperity and the locking strength level of the locking asperity.

Images