A method and system for constructing a four-dimensional dynamic evolution model of an active fault

Abstract

The present invention discloses a method and system for constructing a four-dimensional dynamic evolution model of active faults, comprising: dividing the period since the Late Quaternary into a prehistoric stage, a historical stage, and a present stage; then determining the evolutionary characteristics of active faults in the prehistoric stage, including paleoearthquakes and slip rates; determining the evolutionary characteristics of active faults in the historical stage, including the time of historical earthquakes, macroscopic epicenters, rupture lengths, and earthquake magnitudes; determining the evolutionary characteristics of active faults in the present stage, including the regional crustal strain accumulation rate, interseismic motion closure status, modern earthquake slip distribution characteristics, and post-earthquake afterslip distribution characteristics; and constructing a four-dimensional dynamic evolution model of active faults in chronological order based on the evolutionary characteristics of active faults in the prehistoric stage, the historical stage, and the present stage. The present invention reveals the evolution and changes of active blocks at different time scales in three-dimensional space and the kinematic processes of active faults, providing a scientific basis for studying and assessing large earthquake hazard zones.

Description

Technical Field

[0001] The present invention belongs to the field of earthquake disaster prevention, and in particular relates to a method and system for constructing a four-dimensional dynamic evolution model of an active fault. Background Art

[0002] Active faults are those that have been active since the Late Quaternary Period (~100,000 years ago). These faults are still active today, or have been active in recent geological periods and may become active again in the future.

[0003] The four-dimensional dynamic evolution model of active faults refers to the model of changes in the geometric, kinematic, and physical properties of faults since the Late Quaternary (~100,000 years ago), specifically including changes in the 3D geometric shape of the fault, seismic activity, slip direction, and rate.

[0004] Instrumental earthquakes, or modern earthquakes, are earthquakes recorded by seismic instruments. Historical earthquakes are earthquakes that occurred before instrumental records and are confirmed based on historical documents. Paleoearthquakes are earthquakes that have no written records but are discovered through geological methods. The geological slip rate is the relative slip or displacement velocity calculated based on the magnitude of the fault's displacement and the age of the landform surface. Surface deformation refers to surface movement and deformation observed using methods such as GNSS, InSAR, and leveling.

[0005] While there are numerous publications on constructing three-dimensional models of active faults, there are no reports on four-dimensional dynamic evolution models of active faults. There is still a lack of a sound and effective method for building reliable four-dimensional dynamic evolution models of active faults using existing seismic, geological, geophysical, and geodetic data. Summary of the Invention

[0006] In view of the above problems, the present invention proposes a method and system for constructing a four-dimensional dynamic evolution model of an active fault.

[0007] The present invention adopts the following technical solution: a method for constructing a four-dimensional dynamic evolution model of an active fault, comprising the following steps:

[0008] Step S1: Divide the period since the Late Quaternary into three stages: prehistoric stage, historical stage, and present stage;

[0009] Step S2: determining the evolution characteristics of active faults in the prehistoric period, including paleoearthquakes and slip rates;

[0010] Step S3: Determine the evolution characteristics of active faults in the historical stage, including the time of historical earthquakes, macroscopic epicenters, rupture lengths, and earthquake magnitudes;

[0011] Step S4: Determine the evolution characteristics of active faults at the current stage, including the regional crustal strain accumulation rate, interseismic motion locking status, modern earthquake slip distribution characteristics, and post-earthquake afterslip distribution characteristics;

[0012] Step S5: Based on the evolutionary characteristics of active faults in the prehistoric stage, the historical stage, and the present stage, a four-dimensional dynamic evolution model of the active fault is produced using animation software in chronological order.

[0013] Furthermore, the prehistoric stage is a scale of 100,000 to 10,000 years ago, the historical stage is a scale of 1,000 to 100 years ago, and the present stage is the time period since instrumental records were available.

[0014] Furthermore, step S2 includes the following sub-steps:

[0015] Step S2.1: Collect literature and extract active fault paleoseismic data and slip rate data;

[0016] Step S2.2: If paleoseismic data is lacking, trenching is used to excavate paleoseismic trenches, identify paleoseismic events, and use dating techniques to constrain the time of earthquake occurrence and establish a paleoseismic event time series.

[0017] Step S2.3: If the slip rate data is missing, the fault slip rate is calculated using the displacement of the fault surface and the cumulative time of the displacement;

[0018] Step S2.4: For thrust faults and normal faults, calculate the magnitude of the paleoearthquake based on the coseismic displacement revealed by the trenching using the displacement-magnitude statistical relationship;

[0019] Step S2.5: For strike-slip faults, the possible rupture length of the paleoearthquake is inferred based on the geometric segmentation of adjacent exploration trenches or faults, and the magnitude of the paleoearthquake is calculated based on the statistical relationship between rupture length and magnitude.

[0020] Furthermore, step S3 includes the following sub-steps:

[0021] Step S3.1: Collect historical earthquake catalogs and analyze historical earthquake data to determine the time of historical earthquakes;

[0022] Step S3.2: Using the historical earthquake catalog and other relevant data, conduct a survey of coseismic surface ruptures and secondary fractures associated with historical earthquakes and paleoseismic trenching to determine the macroscopic epicenters and seismogenic faults of historical earthquakes.

[0023] Step S3.3: For historical earthquakes with surface ruptures, determine the length of their coseismic surface ruptures by investigating their coseismic surface ruptures. For historical earthquakes without surface ruptures but for which a complete intensity distribution can be determined, determine the rupture length of the historical earthquakes using the empirical relationship between earthquake intensity and rupture extension. For historical earthquakes without surface ruptures and for which a complete intensity distribution cannot be determined, supplement paleoseismic trenching studies based on damage and felt information from historical materials, combined with earthquake intensity attenuation patterns and active faults, to approximately infer the location and extension of the historical earthquake ruptures.

[0024] Step S3.4: Determine the magnitude of historical earthquakes based on the rupture length-magnitude statistical relationship.

[0025] Furthermore, the step S3.3 of determining the rupture length of historical earthquakes using the earthquake intensity-rupture extension relationship specifically includes the following steps:

[0026] (1) Based on the written records of historical earthquakes, combined with the distribution of active faults and modern small earthquakes, the intensity distribution of relevant historical earthquake events that have been published was verified, and any erroneous intensity distribution was corrected. The revised intensity distribution of historical earthquake events was plotted on the active fault distribution map of the study area.

[0027] (2) Analyze the intensity distribution of historical earthquakes and their relationship with active faults in the earthquake zone to determine the earthquake-causing fault;

[0028] (3) For each maximum intensity I h For historical earthquake events with a magnitude of ≥ VIII, based on their intensity distribution, the empirical relationship between earthquake intensity and rupture extension is used to determine the highest intensity of the event. h The corresponding intensity range of rupture extension [I h , I L '], circle the intensity ≥I L ', and the rupture position and possible maximum extension of the event are preliminarily determined by the range; h The highest intensity, I L ' is the lower limit of the average intensity at both ends of the rupture;

[0029] (4) Considering the geometric structure of the earthquake fault and combining the relevant surface rupture information, the intensity ≥I L ' range is further modified to make its extension along the seismogenic fault closer to the intensity range of the actual rupture extension [I h , I L ] is the relatively heavy damage area; L is the average intensity at both ends of the rupture;

[0030] (5) The position and length of the relatively heavy damage zone along the earthquake fault are used to represent the position and extension of the rupture of the corresponding historical earthquake event. At this time, the relatively heavy damage zone is the rupture zone.

[0031] Furthermore, step S4 includes the following sub-steps:

[0032] Step S4.1: Using the crustal movement velocity observed by GNSS and various continuous deformation models, the crustal strain accumulation rate field in the area where the active fault is located is obtained through inversion;

[0033] Step S4.2: Using the crustal movement velocity observed by GNSS and InSAR and the 3D active block model, the 3D interseismic kinematic locking state of the active fault is obtained by inversion;

[0034] Step S4.3: Using surface deformation and seismic waveform data observed by GNSS and InSAR, the slip distribution characteristics of the current earthquake on the active fault are obtained based on the dislocation theory inversion;

[0035] Step S4.4: Using the surface deformation observed by GNSS and InSAR, the afterslip distribution characteristics on the active fault after the earthquake are obtained based on the dislocation theory.

[0036] Furthermore, the step 5 includes:

[0037] Step S5.1: displaying the 3D fault model on a 3D visualization platform;

[0038] Step S5.2: Display the evolution characteristics of the 3D fault in the prehistoric stage on a 3D visualization platform, including the temporal and spatial distribution of paleoearthquakes and slip rate;

[0039] Step S5.3: Display the evolution characteristics of the 3D fault history stage on a 3D visualization platform, including the temporal and spatial distribution of historical earthquakes;

[0040] Step S5.4: Display the evolution characteristics of the 3D fault at the current stage on the 3D visualization platform, including the regional crustal strain accumulation rate, interseismic motion locking status, modern earthquake slip distribution, and post-earthquake afterslip distribution.

[0041] In another aspect, the present invention provides a system for constructing a four-dimensional dynamic evolution model of an active fault, comprising:

[0042] Stage division module: It is used to divide the period since the Late Quaternary into three stages: prehistoric stage, historical stage and present stage;

[0043] Prehistoric stage evolution characteristics determination module: It is used to determine the evolution characteristics of active faults in the prehistoric stage, including paleoearthquakes and slip rates;

[0044] Historical stage evolution characteristics determination module: It is used to determine the historical stage active fault evolution characteristics, including historical earthquake occurrence time, macroscopic epicenter, rupture length and earthquake magnitude;

[0045] Current stage evolution characteristics determination module: It is used to determine the current stage evolution characteristics of active faults, including regional crustal strain accumulation rate, interseismic motion locking state, modern earthquake slip distribution characteristics and post-earthquake afterslip distribution characteristics;

[0046] Four-dimensional dynamic evolution model construction module: It is used to create a four-dimensional dynamic evolution model of active faults using animation software in chronological order based on the evolutionary characteristics of active faults in prehistoric, historical and present stages;

[0047] The active fault four-dimensional dynamic evolution model construction system is used to execute the steps in the active fault four-dimensional dynamic evolution model construction method.

[0048] Compared with the prior art, the present invention has the following beneficial effects:

[0049] The four-dimensional dynamic evolution model of active faults constructed by the present invention can reveal the evolution and changes of active blocks at different time scales in three-dimensional space, as well as the kinematic process of active faults, and can provide a scientific basis for studying and evaluating large earthquake risk zones. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Various other advantages and benefits will become apparent to those skilled in the art upon reading the detailed description of the preferred embodiment below. The accompanying drawings are for illustration purposes only and are not to be considered as limiting the present invention. The same reference symbols are used throughout the drawings to represent the same components. In the drawings:

[0051] Figure 1 A flow chart for determining the evolutionary characteristics of active faults in the prehistoric period according to an embodiment of the present invention;

[0052] Figure 2 A flow chart for determining the evolutionary characteristics of active faults during the historical phase of an embodiment of the present invention;

[0053] Figure 3 A flowchart for determining the evolution characteristics of active faults in the current stage of an embodiment of the present invention;

[0054] Figure 4 The distribution diagram of paleoearthquakes on the Xianshuihe fault in different time periods according to the embodiment of the present invention is shown as follows: (a) distribution of paleoearthquakes on the Xianshuihe fault from 20ka to 4ka years ago, and (b) distribution of paleoearthquakes on the Xianshuihe fault from 4ka to 0.7ka years ago.

[0055] Figure 5This is a three-dimensional view of the slip rate of the Xianshuihe fault at different time periods according to an embodiment of the present invention. (a) and (b) are the slip rates of the Xianshuihe fault north of Kangding since ~12.5 ka and ~22 ka, respectively. (c) is the slip rate of the Xianshuihe fault south of Kangding since ~11 ka.

[0056] Figure 6 This is a schematic diagram of the relationship between the maximum intensity and rupture extension intensity range of 10 earthquakes in and near the main fault zone on the eastern boundary of the Sichuan-Yunnan block in an embodiment of the present invention;

[0057] Figure 7 This is an intensity distribution map of some strong and large earthquakes in the Xianshuihe fault zone according to an embodiment of the present invention;

[0058] Figure 8 This is a schematic diagram of the four-dimensional dynamic evolution of the Xianshuihe Fault Zone during the historical stages according to an embodiment of the present invention;

[0059] Figure 9 This is a diagram showing the GNSS velocity field and main fault distribution in the Sichuan-Yunnan region according to an embodiment of the present invention;

[0060] Figure 10 This is a schematic diagram of covariance function selection according to an embodiment of the present invention;

[0061] Figure 11 This is a distribution diagram of the second invariant of strain rate in the Sichuan-Yunnan region according to an embodiment of the present invention;

[0062] Figure 12 This is a schematic diagram of the interseismic velocity field according to an embodiment of the present invention;

[0063] Figure 13 Schematic diagram of the double-layer crust structure model in western Sichuan according to an embodiment of the present invention;

[0064] Figure 14 This is a schematic diagram of a rupture segmentation model of the Xianshui River-Anning River-Zemu River fault zone according to an embodiment of the present invention;

[0065] Figure 15 Schematic diagram of the sliding loss rate input model of the present invention and the output model obtained by inversion using different weight ratios-smoothing factors;

[0066] Figure 16 Schematic diagram of the interseismic kinematic locking characteristics of the Xianshuihe fault according to an embodiment of the present invention, including (a) long-term slip rate distribution, (b) interseismic coupling ratio distribution, and (c) interseismic creep rate distribution.

[0067] Figure 17 This is a panoramic schematic diagram of InSAR co-seismic deformation according to an embodiment of the present invention;

[0068] Figure 18Figure 2 shows the InSAR interference fringe patterns and fitting results of an embodiment of the present invention. (a) Interference fringe pattern, (b) Line-of-sight deformation map after resampling, (c) and (d) are the fitted InSAR interference fringe pattern and residual fringes, respectively. The color scale used in (a) is the same.

[0069] Figure 19 This is a schematic diagram of the fitting error of the observation data according to an embodiment of the present invention, wherein the inclination angles of the faults of the compromise curve between the roughness of the slip distribution model and the faults are set at 90°, 82°, and 72°, respectively;

[0070] Figure 20 This is a schematic diagram of the variation of the InSAR data fitting residual with the fault dip angle according to an embodiment of the present invention;

[0071] Figure 21 Schematic diagram of the spatial distribution of coseismic slip according to an embodiment of the present invention, (a) half-space elastic crust model, (b) CRUST2.0 model. DETAILED DESCRIPTION

[0072] Exemplary embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although exemplary embodiments of the present invention are shown in the accompanying drawings, it should be understood that the present invention can be implemented in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the present invention and to fully convey the scope of the present invention to those skilled in the art.

[0073] Example 1

[0074] This embodiment provides a method for constructing a four-dimensional dynamic evolution model of an active fault, comprising the following steps:

[0075] Step S1: Divide the period since the Late Quaternary into three stages: prehistoric stage, historical stage, and present stage;

[0076] Step S2: determining the evolution characteristics of active faults in the prehistoric period, including paleoearthquakes and slip rates;

[0077] Step S3: Determine the evolution characteristics of active faults in the historical stage, including the time of historical earthquakes, macroscopic epicenters, rupture lengths, and earthquake magnitudes;

[0078] Step S4: Determine the evolution characteristics of active faults at the current stage, including the regional crustal strain accumulation rate, interseismic motion locking status, modern earthquake slip distribution characteristics, and post-earthquake afterslip distribution characteristics;

[0079] Step S5: Based on the evolutionary characteristics of active faults in the prehistoric stage, the historical stage, and the present stage, a four-dimensional dynamic evolution model of the active fault is produced using animation software in chronological order.

[0080] The specific process is as follows:

[0081] 1. Such as Figure 1 As shown in Figure 2, the specific process of determining the evolutionary characteristics of active faults in the prehistoric period is as follows:

[0082] 1.1 Data Collection

[0083] 1.1.1 Paleoseismic data collection

[0084] Paleoseismic data is an important parameter reflecting the intensity of active fault activity. In addition to extracting paleoseismic data from literature, paleoseismic data can also be obtained by conducting trenching directly on active faults. The specific process is as follows:

[0085] Step 1: Data collection and comprehensive analysis of the current status of active fault paleoearthquake research;

[0086] Step 2: Image interpretation: Comprehensively utilize multi-source satellite image data to identify the micro-geomorphological environment suitable for excavating paleoseismic trenches and preliminarily determine the trench excavation location;

[0087] Step 3: Field survey, on-site analysis of the structural and geomorphological environment of the trenching location, and optimization of the trenching layout and excavation plan;

[0088] Step 4: Topographic mapping. Before excavation, it is best to conduct aerial mapping of the terrain near the preset location of the exploration trench to draw high-resolution topographic data;

[0089] Step 5: Trench excavation. If the road permits, a large excavator is usually rented to excavate the trench. If the road is not accessible, manual excavation can also be used.

[0090] Step 6: Clean the trench. Manually remove the excavator bucket tooth marks on the trench wall and level the trench wall.

[0091] Step 7: Hang lines in the trench and take photos. Generally, hang an engineering grid with a one-meter interval on the wall of the trench and then take photos;

[0092] Step 8: Use stitching software indoors to stitch the trench photos together to create a complete trench stitching photo;

[0093] Step 9: Trench interpretation, which is divided into field interpretation and indoor interpretation, analyzes paleoseismic events and determines event horizons;

[0094] Step 10: Sample collection: Based on the determined event horizon, determine the horizon for collecting age samples. Generally, samples for radiocarbon dating or optically stimulated luminescence dating are collected.

[0095] Step 11: Sample processing and analysis: send the sample to the laboratory, process the sample, and obtain the dating data;

[0096] Step 12: Time constraint of paleo-earthquake events: constrain the time of paleo-earthquake occurrence based on the age results of samples in pre- and post-earthquake sedimentary strata.

[0097] 1.1.2 Sliding rate data collection

[0098] Slip rate is also an important parameter reflecting the activity intensity of active faults. In addition to extracting slip rate from literature, geological surveys can also be conducted on active faults to obtain slip rate data. The specific process is as follows:

[0099] Step 1: Data collection and comprehensive analysis of the current status of research on active fault slip rates;

[0100] Step 2: Image interpretation: Comprehensively utilize multi-source satellite image data to identify fault-fault landforms suitable for slip rate limitation.

[0101] Step 3: Field survey, on-site analysis of the geomorphological characteristics of the dislocation points, and formulation of a dislocation measurement plan and a chronological sample collection plan;

[0102] Step 4: Topographic mapping. In terrain environments with less vegetation cover, drone aerial photography can be used to map the faulted terrain. If vegetation cover is more abundant, ground-based LiDAR or drone LiDAR can be used to scan the faulted terrain.

[0103] Step 5: Terrain data processing: Based on the surveyed data, a digital elevation model (DEM) of the faulted landform is established;

[0104] Step 6: Measure the fault displacement based on DEM data and field survey results;

[0105] Step 7: Sample processing and analysis: send the sample to the laboratory, process the sample, and obtain the dating data;

[0106] Step 8: Slip rate constraint: Calculate the fault slip rate based on the measured displacement values ​​and the age of the relevant landform surface.

[0107] 1.2 Spatial and temporal distribution of paleoearthquakes along the Xianshuihe fault in the prehistoric period

[0108] Taking the Xianshuihe fault as an example, the existing paleo-earthquake data of the fault are summarized in Table 1. Based on the paleo-earthquake data and the three-dimensional fault model, the paleo-earthquake data distribution map of the Xianshuihe fault in different periods was constructed, see Appendix. Figure 4 .

[0109] Table 1. Data of paleoearthquakes on the Xianshuihe fault

[0110]

[0111]

[0112] Note: The latitude and longitude coordinates are the trench locations, which are the epicenter by default. The focal depth is the assigned depth. To distinguish multiple paleoseismic events at the same trench point in the 3D display, the focal depths of adjacent paleoseismic events differ by 2 km.

[0113] 1.3 Slip rate of the Xianshuihe fault in the prehistoric period

[0114] Taking the Xianshuihe fault as an example, we have collected the existing slip rate data of the fault. Based on these data and the 3D fault model, we constructed a 3D view of the slip rate of the Xianshuihe fault at different times. Figure 5 .

[0115] 2. Such as Figure 2 As shown in Figure 2, the specific evolutionary characteristics of active faults in the historical stage are as follows:

[0116] 2.1 Data Collection

[0117] For historical earthquakes (M ≥ 6.5) without instrumental records, the occurrence dates of most historical earthquakes are documented (or can be estimated). The rupture locations and extensions of these historical earthquakes can be reliably determined using the "active tectonic geological survey + intensity distribution" analysis method of Wen et al. (2008). The required data are documented data that can reliably constrain the intensity distribution of historical earthquakes. For earlier historical earthquakes lacking reliable intensity records, their rupture locations and extensions can be roughly inferred based on limited information on damage / felt distribution combined with active tectonic analysis.

[0118] 2.2 Determining the rupture location and extension of historical earthquakes on the Xianshuihe fault

[0119] 2.2.1 Using the “earthquake intensity-rupture extension” relationship to establish the rupture location and extension

[0120] Wen et al. (2008) established an empirical relationship between earthquake intensity and rupture extension on the eastern boundary of the Sichuan-Yunnan block by statistically analyzing the intensity and rupture extension of historical earthquakes on the eastern boundary of the Sichuan-Yunnan block (see Appendix). Figure 6 Based on this empirical relationship between earthquake intensity and rupture extension, the rupture location and extension of historical earthquake events can be comprehensively determined using earthquake intensity distribution (isoseismic maps) and active fault distribution information. The specific methods and steps are:

[0121] (1) Reinterpret historical earthquake damage, casualties, and felt earthquake records. Combined with information on active fault distribution and modern small earthquake distribution, verify the published intensity distribution of relevant historical earthquake events and correct any erroneous intensity distributions (including the highest intensity value). Then, plot the revised intensity distribution of historical earthquake events on the active fault distribution map of the study area.

[0122] (2) Analyze the intensity distribution of historical earthquakes and their relationship with active faults in the earthquake zone to determine the earthquake-causing fault (segment); for earthquakes with small earthquake records, they can be analyzed together with the distribution of modern small earthquakes.

[0123] (3) For each maximum intensity I h For historical earthquake events with a magnitude of ≥ VIII, based on their intensity distribution, the empirical relationship of “earthquake intensity-rupture extension” (Appendix Figure 6 ), determine the highest intensity I h The corresponding intensity range of rupture extension [I h , I L '], circle the intensity ≥I L ' range, and preliminarily determine the rupture location and possible maximum extension of the event based on this range.

[0124] (4) Considering the geometric structure of the earthquake fault and combining the relevant surface rupture information, the intensity ≥I L ' to further limit and modify its range so that its extension along the seismogenic fault is closer to the intensity zone of the actual rupture extension [I h , I L ]——“relatively heavy damage area”.

[0125] (5) The location and length of the “relatively severely damaged zone” along the seismogenic fault are used to represent the location and extension of the rupture of the corresponding historical earthquake event. In this case, the “relatively severely damaged zone” is the “rupture zone”.

[0126] Attachment Figure 6 The relationship between the maximum intensity and rupture extension intensity of the 10 earthquakes in and near the main fault zone on the eastern boundary of the Sichuan-Yunnan block (according to the M7 Special Working Group, 2012)

[0127] The serial numbers indicate earthquake events with good data records; the intensity range covered by the double-headed arrows is [I h , I L ], which corresponds to the actual extension range of the rupture along the seismogenic fault; the intensity range between the two outer envelopes is [I h , I L’ ], representing the maximum possible extension of the rupture.

[0128] 2.2.2 Reanalysis of historical earthquake intensity / impact distribution

[0129] On the basis of systematically sorting out the historical textual data of earthquakes in the Xianshuihe fault zone, referring to the standards of the earthquake intensity scale, analyzing, confirming or revising the published intensity distribution, and finally drawing the earthquake intensity distribution map superimposed with the detailed Xianshuihe fault distribution (see Appendix). Figure 7 ).

[0130] 3. Such as Figure 3 As shown in Figure 2, the specific process of determining the evolutionary characteristics of active faults at the current stage is as follows:

[0131] 3.1 Data Collection

[0132] 3.1.1 GNSS data

[0133] GNSS can accurately monitor crustal deformation and static co-seismic deformation. The GNSS observation network can simultaneously monitor surface deformation caused by tectonic and non-tectonic movements. When using GNSS observations for surface deformation monitoring, it is necessary to measure the position of the GNSS station before and after the deformation occurs. For example, in the monitoring of earthquake deformation measurements, the position of the same GNSS station before and after the earthquake is generally solved to obtain the co-seismic displacement caused by the earthquake. In order to make the results more accurate, the time series of the station is usually fitted using long-term observations before and after the earthquake. At this time, the position where the time series jumps is the displacement caused by the earthquake. Taking GAMIT software as an example, the data processing flow is roughly as follows:

[0134] Step 1: First, use GAMIT software to solve the GNSS observation data and IGS data together to obtain a single-day relaxation solution. The antenna phase center model uses the absolute PCVs model provided by IGS, the ocean tide model uses FES2004, and the tropospheric mapping function uses GMF.

[0135] Step 2: When solving the time series of continuous stations, use QOCA software to combine all single-day solutions for global adjustment. Select stations in the IGS station as frame points. Use these frame points to solve the similarity transformation parameters relative to the global reference frame ITRF2008, thereby obtaining the coordinate time series of each station under ITRF2008.

[0136] Step 3: When solving the velocity field, use GLOBK software to convert each subnet h file into a glx file, then merge multiple glx files for each day into one GLX file, and then use QOCA software to perform whole-network adjustment.

[0137] 3.1.2 InSAR data

[0138] The scientific community now routinely uses InSAR to monitor continental crustal movement and geological hazards. Commonly used InSAR data processing software include commercial software such as GAMMA, SARscape, ENVI, ROI_PAC, and the open-source GMTSAR. Using GAMMA as an example, we describe the data processing workflow for extracting surface deformation using DInSAR technology:

[0139] Step 1: SAR image registration. The process generally includes three steps: coarse registration, fine registration, and transformation from image grid coordinates and resampling.

[0140] Step 2: Generate the interferogram. Before performing conjugate multiplication to obtain the interferogram, the SAR image is generally subjected to noise reduction processing to improve the signal-to-noise ratio of the interferogram.

[0141] Step 3: Estimation of the interferometric baseline.

[0142] Step 4: DEM simulation of SAR image.

[0143] Step 5: Differential processing.

[0144] Step 6: Interferogram phase filtering.

[0145] Step 7: Phase unwrapping.

[0146] Step 8: Geocoding.

[0147] 3.2 Determining the crustal strain rate in the Xianshuihe fault area

[0148] Regional crustal strain rate is very important for understanding the strain accumulation state of active faults and assessing earthquake risks. We calculated the strain rate field in the Sichuan-Yunnan region where the Xianshuihe fault is located based on the least squares collocation method. The observation data used were the GNSS interseismic velocity field from 1999 to 2016. The coordinate reference frame was the Eurasian continent, and a total of 450 stations were included (see Appendix). Figure 9 ).

[0149] When the least squares collocation method is used to establish the velocity and strain fields of crustal movement, the surface observed velocity is composed of the rigid motion of the block (dip parameter) and the random strain field (random signal) at the station. The specific calculation process is as follows:

[0150] Step 1: Calculate the statistics of variance and covariance functions based on the given observation data (Appendix 2). Set the distance interval and the number of intervals, and count the number of stations that fall within different distance ranges. Commonly used candidate functions for covariance functions include Gaussian function, Silver function and exponential function (Appendix 2). Figure 10 ).

[0151] Step 2: Select the function that best fits the data as the covariance function, and use the statistical values ​​of variance and covariance as samples to invert the parameters of the covariance function.

[0152] Step 3: Use the covariance function to build the variance-covariance matrix.

[0153] Step 4: Solve the observation equation and derive the random signal component at the station to obtain the strain field (see Appendix Figure 11 ).

[0154] Appendix 2 Sampling interval and number of intervals selection

[0155]

[0156] The results show that high strain accumulation zones exist in the northern Xianshuihe and Xiaojiang fault zones, reaching ~150 nanostrain / yr. Strain accumulation of ~80 nanostrain / yr also occurs in the northern Honghe and southern Heihe fault zones. The analysis of the continuous deformation model derived from the least squares collocation method provides an important foundation for constructing a four-dimensional dynamic evolution model of active faults.

[0157] 3.3 Determining the interseismic locking state of the Xianshuihe fault

[0158] We further used GNSS interseismic velocity field and InSAR data (Appendix Figure 12 ) to determine the interseismic motion locking state of the Xianshuihe fault. The specific process is as follows:

[0159] Step 1: Based on existing seismological and crustal flow studies, we divide the crust of the study area into a 20 km thick elastic upper crust and a 30 km thick middle and lower crustal flow, and use them as Maxwell bodies (see Appendix Figure 13 ).

[0160] The spatiotemporal distribution of ancient and historical earthquakes shows that the Xianshuihe-Anninghe-Zemuhe fault zone has obvious segmented earthquake rupture characteristics. The magnitude and recurrence interval of characteristic earthquakes in each segment are listed in Appendix 3.

[0161] Step 2: Based on existing research, we divide the Xianshuihe fault zone into 9 sections (see Appendix Figure 14 The surface trace of the Xianshuihe-Anninghe-Zemuhe fault zone was extracted, and a 3D geometric model of the fault was constructed using triangular dislocation elements (see Appendix). Figure 15 ).

[0162] Step 3: By adjusting the optimal weight ratio and smoothing coefficient (see Figure 15 ), using measured data to invert the rotation and strain parameters of the active block, as well as the slip loss rate and locking characteristic distribution of the fault.

[0163] Table 3 Magnitudes and recurrence periods of characteristic earthquakes on the nine segments of the Xianshuihe-Anninghe-Zemuhe fault zone

[0164]

[0165] The inversion results show that there are 6 completely closed potential asperities in the Xianshuihe-Anninghe-Zemuhe fault zone, including one each between Zhuwo-Luhuo, at the bottom of Kangding, between Shimian-Mianning, and near Qiaojia, and two between Mianning and Xichang (see Appendix). Figure 16Additionally, a partially coupled potential asperity exists between Daofu and Kangding, with a length exceeding 140 km, exceeding the minimum length resolvable by existing geodetic data. Seismic moment accumulation calculations indicate that the Songlinkou-Selaha segment (S4), the Shimian-Mianning segment (S7), and the Mianning-Xichang segment (S8) possess very high seismic rupture risks. These results reflect the current kinematic locked nature of the Xianshuihe fault.

[0166] 3.4 Determination of the distribution characteristics of instrumental seismic slip on the Xianshuihe fault

[0167] Between November 22 and 29, 2014, more than 1,000 earthquakes occurred on the Xianshuihe fault northwest of Kangding, Sichuan, China, including the largest event on November 22, 2014, and a secondary event on November 25, 2014 (Appendix Table 4).

[0168] Table 4 Source parameters of the two largest Kangding earthquakes in November 2014

[0169]

[0170] We obtained the surface deformation of the Kangding earthquake by processing radar images from the ALOS-2 PALSAR-2 satellite. The details are as follows: Using the Caltech / JPL software ROI_PAC (version 3.1beta) in dual-track differential interferometry mode, we obtained the coseismic deformation field from a pair of ascending ALOS-2 PALSAR-2 images (see Appendix). Figure 17 The two images were acquired on September 26, 2014, and December 5, 2014, respectively, with a vertical baseline of 12.6 m. The interferograms were corrected for satellite position differences using JAXA's precise orbital parameters, then phase unwrapped using the SNAPHU program and finally converted to the WGS-84 coordinate system (see Appendix). Figure 18 ).

[0171] On this basis, the co-seismic slip on each slider was further inverted using the bounded least squares algorithm and 593 resampled InSAR observation data. The elastic half-space dislocation theory was used to calculate the Green's function between fault slip and surface displacement during the inversion. In order to find the optimal Laplace smoothing operator, the fault dip was first set to 90°. The curve between the root mean square error (RMSE) of the data fitting and the roughness of the slip distribution shows that the optimal Laplace smoothing operator is 0.1 (Appendix Figure 19 Then, the fault dip angle is changed and the inversion is performed under the optimal Laplace smoothing operator. The results show that when the fault dip angle is 82°, the fitting residual is the smallest (see Appendix Figure 20In addition, more inversions show that the optimal smoothing factor selected when the fault dip is 90° is also applicable to non-vertical fault models (Appendix Figure 19 ).

[0172] Attachment Figure 21 The figure shows the optimal deep slip distribution model obtained. The main slip occurred on the S1 fault, with a maximum value of 0.47 m occurring at a depth of ~9 km. However, there was no significant slip on the S2 fault, especially near the epicenter of the November 25 earthquake. Based on the slip distribution model obtained by inversion, it was calculated that the total seismic moment released by the two secondary earthquake sequences reached 2.36×10 18 Nm, equivalent to an earthquake of Mw 6.2 (stiffness 30GPa). Figure 18 (c) and attached Figure 18 (d) The fitted InSAR interferogram and residual fringes, respectively. No significant residual fringes are observed near the seismogenic faults (S1 and S2). The RMS residual between the InSAR observations and the predicted LOS displacement is 0.4 cm.

[0173] Example 2

[0174] This embodiment provides a system for constructing a four-dimensional dynamic evolution model of an active fault, including:

[0175] Stage division module: It is used to divide the period since the Late Quaternary into three stages: prehistoric stage, historical stage and present stage;

[0176] Prehistoric Evolution Characteristic Determination Module: This module is used to determine the evolutionary characteristics of active faults in the prehistoric period, including paleoearthquakes and slip rates. Based on paleoearthquake data and a 3D fault model, paleoearthquake data distribution maps for different time periods of the fault are constructed. Based on slip rate data and a 3D fault model, a 3D view of slip rates for different time periods of the fault is constructed.

[0177] Historical stage evolution characteristics determination module: This module is used to determine the historical stage evolution characteristics of active faults, including the time of historical earthquakes, macro-epicenters, rupture lengths, and earthquake magnitudes. Based on the historical stage evolution characteristics of active faults, the module displays the rupture zones with time attributes on a 3D visualization platform, and obtains a 4D historical stage evolution model of the faults.

[0178] Current stage evolution characteristics determination module: It is used to determine the current stage evolution characteristics of active faults, including the regional crustal strain accumulation rate, interseismic motion closure state, modern earthquake slip distribution characteristics and post-earthquake afterslip distribution characteristics. Based on the current stage evolution characteristics of active faults, a four-dimensional evolution model of the current stage is constructed.

[0179] Four-dimensional dynamic evolution model construction module: It is used to create a four-dimensional dynamic evolution model of active faults using animation software in chronological order based on the evolutionary characteristics of active faults in the prehistoric stage, historical stage and present stage.

[0180] Similarly, it should be understood that in order to streamline the present invention and aid in understanding one or more of the various inventive aspects, in the above description of exemplary embodiments of the invention, various features of the embodiments of the invention are sometimes grouped together into a single embodiment, figure, or description thereof. However, this disclosed method should not be interpreted as reflecting an intention that the claimed invention requires more features than are expressly recited in each claim. Rather, as reflected in the claims, inventive aspects lie in less than all the features of the individual embodiments previously disclosed. Accordingly, the claims that follow the detailed description are hereby expressly incorporated into this detailed description, with each claim standing on its own as a separate embodiment of the invention.

[0181] Furthermore, those skilled in the art will appreciate that although some embodiments herein include certain features included in other embodiments but not other features, combinations of features from different embodiments are intended to be within the scope of the present invention and to form different embodiments. For example, any of the claimed embodiments may be used in any combination.

[0182] The various component embodiments of the present invention can be implemented in hardware, or in software modules running on one or more processors, or in a combination thereof. It will be appreciated by those skilled in the art that a microprocessor or digital signal processor (DSP) can be used in practice to implement some or all of the functions of some or all of the components according to an embodiment of the present invention. The present invention can also be implemented as a device or apparatus program (e.g., a computer program and a computer program product) for executing a part or all of the methods described herein. Such a program implementing the present invention can be stored on a computer-readable medium, or can have the form of one or more signals. Such a signal can be downloaded from an Internet website, or provided on a carrier signal, or provided in any other form.

[0183] It should be noted that the above embodiments illustrate rather than limit the invention, and that alternative embodiments may be devised by a person skilled in the art without departing from the scope of the appended claims. In the claims, any reference signs placed between brackets should not be construed as limiting the claims. The word "comprising" does not exclude the presence of elements or steps not listed in the claims. The word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements. The present invention may be implemented by means of hardware comprising several different elements and by means of appropriately programmed computers. In a unit claim enumerating several means, several of these means may be embodied by the same item of hardware. The use of the words first, second, and third etc. does not indicate any order. These words may be interpreted as names. The steps in the above embodiments should not be understood as limiting the order of execution unless otherwise specified.

Claims

1. A method for constructing a four-dimensional dynamic evolution model of an active fault, characterized in that: The steps include: Step S1: Divide the period since the Late Quaternary into three stages: prehistoric stage, historical stage, and present stage; Step S2: determining the evolution characteristics of active faults in the prehistoric period, including paleoearthquakes and slip rates; Step S3: Determine the evolution characteristics of active faults in the historical stage, including the time of historical earthquakes, macroscopic epicenters, rupture lengths, and earthquake magnitudes; Step S4: Determine the evolution characteristics of active faults at the current stage, including the regional crustal strain accumulation rate, interseismic motion locking status, modern earthquake slip distribution characteristics, and post-earthquake afterslip distribution characteristics; Step S5: Based on the evolutionary characteristics of active faults in the prehistoric stage, the historical stage, and the present stage, a four-dimensional dynamic evolution model of the active fault is produced using animation software in chronological order.

2. The method for constructing a four-dimensional dynamic evolution model of an active fault according to claim 1, characterized in that: The prehistoric stage is a scale of 100,000 to 10,000 years ago, the historical stage is a scale of 1,000 to 100 years ago, and the present stage is the time period since instrumental records.

3. The method for constructing a four-dimensional dynamic evolution model of an active fault according to claim 1, characterized in that: The step S2 includes the following sub-steps: Step S2.1: Collect literature and extract active fault paleoseismic data and slip rate data; Step S2.2: If paleoseismic data is lacking, trenching is used to excavate paleoseismic trenches, identify paleoseismic events, and use dating techniques to constrain the time of earthquake occurrence and establish a paleoseismic event time series. Step S2.3: If the slip rate data is missing, the fault slip rate is calculated using the displacement of the fault surface and the cumulative time of the displacement; Step S2.4: For thrust faults and normal faults, calculate the magnitude of the paleoearthquake based on the coseismic displacement revealed by the trenching using the displacement-magnitude statistical relationship; Step S2.5: For strike-slip faults, the possible rupture length of the paleoearthquake is inferred based on the geometric segmentation of adjacent exploration trenches or faults, and the magnitude of the paleoearthquake is calculated based on the statistical relationship between rupture length and magnitude.

4. The method for constructing a four-dimensional dynamic evolution model of an active fault according to claim 1, characterized in that: The step S3 includes the following sub-steps: Step S3.1: Collect historical earthquake catalogs and analyze historical earthquake data to determine the time of historical earthquakes; Step S3.2: Using the historical earthquake catalog and other relevant data, conduct a survey of coseismic surface ruptures and secondary fractures associated with historical earthquakes and paleoseismic trenching to determine the macroscopic epicenters and seismogenic faults of historical earthquakes. Step S3.3: For historical earthquakes with surface ruptures, determine the length of their coseismic surface ruptures by investigating their coseismic surface ruptures. For historical earthquakes without surface ruptures but for which a complete intensity distribution can be determined, determine the rupture length of the historical earthquakes using the empirical relationship between earthquake intensity and rupture extension. For historical earthquakes without surface ruptures and for which a complete intensity distribution cannot be determined, supplement paleoseismic trenching studies based on damage and felt information from historical materials, combined with earthquake intensity attenuation patterns and active faults, to approximately infer the location and extension of the historical earthquake ruptures. Step S3.4: Determine the magnitude of historical earthquakes based on the rupture length-magnitude statistical relationship.

5. The method for constructing a four-dimensional dynamic evolution model of an active fault according to claim 4, characterized in that: Determining the rupture length of historical earthquakes using the earthquake intensity-rupture extension relationship in step S3.3 specifically includes the following steps: (1) Based on the written records of historical earthquakes, combined with the distribution of active faults and modern small earthquakes, the intensity distribution of relevant historical earthquake events that have been published was verified, and any erroneous intensity distribution was corrected. The revised intensity distribution of historical earthquake events was plotted on the active fault distribution map of the study area. (2) Analyze the intensity distribution of historical earthquakes and their relationship with active faults in the earthquake zone to determine the earthquake-causing fault; (3) For each maximum intensity I h For historical earthquake events with a magnitude of ≥ VIII, based on their intensity distribution, the empirical relationship between earthquake intensity and rupture extension is used to determine the highest intensity of the event. h The corresponding intensity range of rupture extension [I h , I L '], circle the intensity ≥I L ', and the rupture position and possible maximum extension of the event are preliminarily determined by the range; h The highest intensity, I L ' is the lower limit of the average intensity at both ends of the rupture; (4) Considering the geometric structure of the earthquake fault and combining the relevant surface rupture information, the intensity ≥I L ' range is further modified to make its extension along the seismogenic fault closer to the intensity range of the actual rupture extension [I h , I L ] is the relatively heavy damage area; L is the average intensity at both ends of the rupture; (5) The position and length of the relatively heavy damage zone along the earthquake fault are used to represent the position and extension of the rupture of the corresponding historical earthquake event. At this time, the relatively heavy damage zone is the rupture zone.

6. The method for constructing a four-dimensional dynamic evolution model of an active fault according to claim 1, characterized in that: The step S4 includes the following sub-steps: Step S4.1: Using the crustal movement velocity observed by GNSS and various continuous deformation models, the crustal strain accumulation rate field in the area where the active fault is located is obtained through inversion; Step S4.2: Using the crustal movement velocity observed by GNSS and InSAR and the 3D active block model, the 3D interseismic kinematic locking state of the active fault is obtained by inversion; Step S4.3: Using surface deformation and seismic waveform data observed by GNSS and InSAR, the slip distribution characteristics of the current earthquake on the active fault are obtained based on the dislocation theory inversion; Step S4.4: Using the surface deformation observed by GNSS and InSAR, the afterslip distribution characteristics on the active fault after the earthquake are obtained based on the dislocation theory.

7. The method for constructing a four-dimensional dynamic evolution model of an active fault according to claim 1, characterized in that: The step 5 comprises: Step S5.1: displaying the 3D fault model on a 3D visualization platform; Step S5.2: Display the evolution characteristics of the 3D fault in the prehistoric stage on a 3D visualization platform, including the temporal and spatial distribution of paleoearthquakes and slip rate; Step S5.3: Display the evolution characteristics of the 3D fault history stage on a 3D visualization platform, including the temporal and spatial distribution of historical earthquakes; Step S5.4: Display the evolution characteristics of the 3D fault at the current stage on the 3D visualization platform, including the regional crustal strain accumulation rate, interseismic motion locking status, modern earthquake slip distribution, and post-earthquake afterslip distribution.

8. A system for constructing a four-dimensional dynamic evolution model of an active fault, characterized in that: include: Stage division module: It is used to divide the period since the Late Quaternary into three stages: prehistoric stage, historical stage and present stage; Prehistoric stage evolution characteristics determination module: It is used to determine the evolution characteristics of active faults in the prehistoric stage, including paleoearthquakes and slip rates; Historical stage evolution characteristics determination module: It is used to determine the historical stage active fault evolution characteristics, including historical earthquake occurrence time, macroscopic epicenter, rupture length and earthquake magnitude; Current stage evolution characteristics determination module: It is used to determine the current stage evolution characteristics of active faults, including regional crustal strain accumulation rate, interseismic motion locking state, modern earthquake slip distribution characteristics and post-earthquake afterslip distribution characteristics; Four-dimensional dynamic evolution model construction module: It is used to create a four-dimensional dynamic evolution model of active faults based on the evolution characteristics of active faults in prehistoric, historical and current stages in chronological order using animation software; The active fault four-dimensional dynamic evolution model construction system is used to execute the steps in the active fault four-dimensional dynamic evolution model construction method according to any one of claims 1 to 7.

Images