A method for earthquake rupture inversion based on multi-source observation data

Abstract

The present invention relates to an earthquake rupture inversion method based on multi-source observation data, including: obtaining a focal mechanism solution of a specified earthquake event; obtaining geometric parameters of the fault of the earthquake event according to the focal mechanism solution; generating a three-dimensional grid model of the fault plane to which the fault belongs according to a pre-constructed geometric shape model of the fault and the geometric parameters of the fault; obtaining the slip of the sub-fault to which the earthquake rupture belongs according to the limiting conditions of each sub-fault in the three-dimensional grid model; performing inversion verification with the aid of an inversion model according to the observation data to which the earthquake event belongs and the slip of the sub-fault, and if the verification passes, using the slip of all sub-faults when the verification passes as the earthquake rupture inversion result. The above method improves the stability and reliability of the inversion result, and can better estimate the slip vector at each sub-fault, thereby more comprehensively understanding the slip process of the entire earthquake fault.

Description

Technical Field

[0001] The present invention relates to the field of computer technology, and in particular to an earthquake rupture inversion method based on multi-source observation data. Background Art

[0002] Tsunamis are often referred to as "earthquake tsunamis". my country has vast sea areas, and its coastal areas are within the Pacific Rim seismic belt. Since tsunamis can travel across oceans over long distances and lose less energy during the propagation process, my country has long faced the threat of tsunamis from transoceanic tsunamis around the Pacific Rim and seismic belts such as the Ryukyu Trench and the Manila Trench. The prediction and mitigation of submarine earthquake sources has always been one of the focuses of tsunami warning and mitigation work.

[0003] Earthquake rupture is the process of releasing energy through the rapid sliding of faults after the internal stress of the earth's crust accumulates to a certain level. Earthquakes are essentially the process of rupture and dislocation starting, spreading and ending on the fault plane. The rupture process of a large earthquake usually includes the successive rupture of several asperities or obstacles on the fault plane, and the rupture behavior is relatively complex. In a complex fault system, there may even be a rupture mode in which the earthquake jumps between faults. The range of rupture will vary with the magnitude and focal mechanism, covering an area ranging from a dozen kilometers to hundreds of kilometers; the duration of the rupture can range from a dozen seconds to hundreds of seconds.

[0004] The development of modern geophysical space and imaging geodesy and ocean observation technology has provided unprecedented multi-dimensional observation data for the study of large earthquakes. These data provide important constraints for people to comprehensively study the rupture process of large earthquakes. By using observation data of different dimensions to invert the rupture characteristics of earthquake ruptures, it is possible to provide the tsunami warning system with earthquake sources that are more in line with reality, thereby improving the accuracy of tsunami warning forecasts. In this way, it is possible to effectively reduce casualties and property losses caused by tsunamis, and provide strong guarantees for the safety and stability of human society.

[0005] Traditional earthquake rupture inversion methods mainly rely on seismic waveform inversion technology. Due to the long distance and long travel time of teleseismic waves, it takes about 6 to 13 minutes for P waves to propagate to the station 30 degrees away. In order to invert teleseismic data, it is necessary to try multiple parameter combinations to determine the optimal solution, which greatly affects the inversion efficiency and the offshore tsunami warning with higher time requirements. The multi-source observation data method integrates a variety of geophysical methods, provides constraints for near-field observations, and improves the efficiency, accuracy and resolution of the inversion results.

[0006] Traditional methods mainly use seismic waveform data recorded by seismic stations to determine earthquake rupture parameters through inversion algorithms. This method is relatively single in data source and may be affected by factors such as uneven distribution of stations and seismic wave attenuation. Summary of the invention

[0007] 1. Technical issues to be resolved

[0008] In view of the above-mentioned shortcomings and deficiencies of the prior art, the present invention provides an earthquake rupture inversion method based on multi-source observation data, which improves the stability and reliability of the inversion results and enhances the accuracy of the source parameters.

[0009] (II) Technical solution

[0010] In order to achieve the above object, the main technical solutions adopted by the present invention include:

[0011] In a first aspect, an embodiment of the present invention provides an earthquake rupture inversion method based on multi-source observation data, comprising:

[0012] S100, obtaining a focal mechanism solution of a designated earthquake event;

[0013] S200, obtaining geometric parameters of the fault of the earthquake event according to the focal mechanism solution;

[0014] S300, generating a three-dimensional grid model of the fault plane to which the fault belongs according to the pre-constructed geometric shape model of the fault and the geometric parameters of the fault;

[0015] S400, obtaining the slip amount of the sub-fault to which the earthquake rupture belongs according to the limiting conditions of each sub-fault in the three-dimensional grid model;

[0016] S500: Perform inversion verification based on the observation data of the earthquake event and the slip amount of the sub-faults with the help of the inversion model. If the verification is successful, the slip amount of all sub-faults when the verification is successful is used as the earthquake rupture inversion result.

[0017] Optionally, the geometric parameters of the fault of the specified earthquake event include: strike angle α, dip angle β, and slip angle;

[0018] The focal mechanism solution of the earthquake event is estimated based on the database of historical earthquake events at the same location, including: epicenter location, magnitude, strike angle, dip angle, and slip angle.

[0019] Optionally, the S300 includes:

[0020] The fault geometry model is constructed based on the data of Slab 2.0, a global comprehensive subduction zone geometry model library, and is modified based on the geometric parameters of the fault.

[0021] Specifically, the normal vector of the fault plane is calculated; the strike vector S is expressed by the strike angle α as:

[0022] S=(sinα,cosa,0)

[0023] The dip vector d is expressed by the fault dip angle β as: d = (cosβ, 0, -sinβ)

[0024] Then the normal vector n of the fault plane is: n = d × s = (cosαsinβ, -sinαsinβ, cosαcosβ)

[0025] Assume that there is a given point A on the fault plane: A(x0,y0,z0), and the coordinates of any point B on the fault plane are B(x,y,z), then the fault plane equation is:

[0026] That is, the three-dimensional grid model is z=z0+tanβ(x0-x)-tanαtanβ(y0-y),

[0027] Point A is the point selected as the epicenter.

[0028] Optionally, the S400 includes:

[0029] Assuming the dislocation plane and parameterizing the shape of the slip rate function, the slip inversion is transformed into the following formula (1) with double summation over the nth row and mth column:

[0030]

[0031] Among them, D jk is the slip amplitude, i.e. the slip amount of the sub-fault to be solved, λ jk is the sliding angle, i.e. the parameter of focal mechanism solution, S jk (t) is a given rise time function,

[0032] V jk is the average rupture velocity between the earthquake source and the sub-fault jk, V jk The limiting condition is that the maximum allowable rupture velocity is 80% of the deepest shear wave velocity in the crust velocity model of the fault Crust1.0;

[0033] is the sub-fault Green function; i is 1 or 2, t is time, j and k represent the number of rows and columns of the sub-fault respectively; u(t) is the total slip.

[0034] Optionally, based on the focal mechanism solution and known parameters, obtain the sub-fault Green function The calculation formula of Green's function G(r,θ,z,t) is:

[0035]

[0036] Among them (U z ,U r ,U θ ) is the known surface harmonic vector basis coordinates, k is the known horizontal wave number set in advance, ω is the known circular frequency set in advance, r, θ, z are the cylindrical coordinate system, and t is the time.

[0037] Optionally, before S500, the method further includes:

[0038] Get multi-source observation data belonging to the specified earthquake event;

[0039] The multi-source observation data includes: InSAR data before and after the earthquake, GPS data before and after the earthquake, water level data corresponding to the earthquake event, and seismic wave data corresponding to the earthquake time;

[0040] Calibrate each observation data in the multi-source observation data, and achieve time synchronization and spatial registration to obtain pre-processed multi-source observation data;

[0041] Accordingly, S500 includes:

[0042] According to the preprocessed observation data and the slip amount of the sub-faults, the inversion verification is carried out with the help of the inversion model. If the verification passes, the slip amount of all sub-faults when the verification passes will be used as the earthquake rupture inversion result.

[0043] Optionally, the S500 includes:

[0044] Based on the pre-given weight of each data type in the multi-source observation data, a data vector is formed

[0045] And according to the data vector Generate diagonal weight matrix

[0046] According to the diagonal weight matrix Constructing the diagonal covariance matrix

[0047] Based on the following formula (3) and constraint condition (4), the least square method is used to solve until verification is passed;

[0048]

[0049] in, is the high-frequency GPS-GNSS data vector, is the seismic p-wave data vector, is the earthquake strong motion data vector; the norm symbol |||| is used to represent the norm of the vector,

[0050] is the high-frequency GPS data-GNSS diagonal covariance matrix, is the diagonal covariance matrix of seismic p-wave data, is the diagonal covariance matrix of earthquake strong motion data, G is the Green's function, m is the sub-fault slip, L is the Laplace operator, and β is the smoothing factor.

[0051] Optionally, multi-source observation data belonging to a designated earthquake event is obtained, each observation data in the multi-source observation data is calibrated, and time synchronization and spatial registration are achieved to obtain pre-processed multi-source observation data, including:

[0052] Select the data of seismic stations at 30°-90° from the epicenter of the specified earthquake event, i.e. acceleration.

[0053] Integrate the acceleration to velocity with a 100 to 5 Hz filter,

[0054] A bandpass filter between 0.05 and 0.4 Hz is performed, the velocity waveform is integrated into the displacement waveform, and the time window data of the specified time period from the arrival of the P wave is intercepted as the preprocessed seismic wave data

[0055] and / or,

[0056] Obtain static GPS data of the first duration through GPS stations, where the first duration includes the time period from before the earthquake to after the earthquake;

[0057] The time series of position coordinates is obtained by PPP processing of static GPS data in precise point positioning mode.

[0058] The co-seismic deformation is obtained by the difference between the average value of the position coordinates a few days before the earthquake and the average value of the position coordinates a few days after the earthquake.

[0059] Optionally, obtaining multi-source observation data belonging to a designated earthquake event, calibrating each observation data in the multi-source observation data, and achieving time synchronization and spatial registration to obtain pre-processed multi-source observation data, further includes:

[0060] Acquire raw high-frequency GPS data of a second duration of the GNSS site, where the second duration is less than or equal to the first duration;

[0061] Using the clock error and orbit information provided by the international GNSS service, the XYZ component displacement of the GNSS station is estimated, and high-frequency GPS data is calculated to obtain the time series of displacement at a sampling rate of 5 Hz;

[0062] The displacement time series is windowed, and the co-seismic displacement is calculated from the difference between the post-seismic position and the pre-seismic position in the windowed displacement time series. The co-seismic displacement is low-pass filtered at 0.5 Hz to obtain the final co-seismic displacement.

[0063] and / or,

[0064] Acquire InSAR data, i.e., interferograms spanning a long period before and after an earthquake, and use software tools to process the interferograms spanning a long period before and after an earthquake;

[0065] Use time series analysis method to analyze the interferogram, remove the atmospheric influence, and obtain the corrected interferogram;

[0066] For the corrected interferogram, 90m resolution SRTM DEM data is used to remove the terrain phase effect, and the interferogram without terrain phase effect is obtained; and the co-seismic deformation field is obtained by unwrapping.

[0067] The co-seismic deformation field is downscaled and resampled using the quadtree method to obtain the pre-processed co-seismic deformation field;

[0068] and / or, obtaining water level data from buoys or tide stations and removing tidal information from the water level data.

[0069] In a second aspect, an embodiment of the present invention further provides a computing device, comprising: a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program in the memory to specifically perform the steps of a seismic rupture inversion method based on multi-source observation data as described in any one of the first aspects above.

[0070] (III) Beneficial effects

[0071] The present invention realizes comprehensive monitoring of crustal deformation and motion state by fusing multi-source observation data, including InSAR, GPS, water level, seismic waveform and other data, and provides more abundant data support for the analysis of earthquake rupture process; a concise inversion algorithm is proposed, which can effectively process observation data and dynamically simulate earthquake rupture process, and can reflect the spatiotemporal evolution characteristics of earthquake rupture in real time, improve the stability and reliability of inversion results, and improve the accuracy of source parameters. The method of the present invention can provide a more reliable source for tsunami warning and disaster assessment, and is also of great significance for improving the level of earthquake science research and enhancing earthquake disaster defense capabilities, and has broad application prospects.

[0072] The method of the present invention combines the modern advantages of earth observation technology and ocean observation technology to fully capture the complexity of the earthquake rupture process. This method not only inherits the in-depth understanding of seismic wave propagation and earthquake source characteristics of traditional seismology, but also incorporates the high-precision surface deformation monitoring of earth observation technology and the detailed observation of seabed seismic activity by ocean observation technology. Through this multidisciplinary and multi-technical integration, it is possible to examine the complex process of earthquake rupture from multiple dimensions and multiple levels, and then construct a more detailed and accurate earthquake rupture model. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 A schematic flow chart of a method for earthquake rupture inversion based on multi-source observation data provided by one embodiment of the present invention;

[0074] Figure 2 A schematic diagram of a three-dimensional grid model for generating a fault plane according to the present invention;

[0075] Figure 3 for Figure 2 A schematic diagram of the model from another angle;

[0076] Figure 4A is a schematic diagram of the horizontal component 1 of the first inversion case,

[0077] Figure 4B is a schematic diagram of the horizontal component 2 of the first inversion case;

[0078] Figure 4C is a schematic diagram of the vertical component of the first inversion case;

[0079] Figure 5A For Figure 4A The schematic diagram after the waveform is clipped and band-pass filtered by the horizontal component 1;

[0080] Figure 5B For Figure 4B The schematic diagram after the waveform is clipped and band-pass filtered by the horizontal component 2;

[0081] Figure 5C For Figure 4C Schematic diagram of the vertical component of the waveform clipping and bandpass filtering;

[0082] Figure 6 This is the fault model diagram in the first inversion case;

[0083] Figure 7 An example diagram of the sub-fault Green's function in the first inversion case;

[0084] Figure 8 Schematic diagram of the comparison between the predicted earthquake P-wave record and the synthetic seismogram in the first inversion case;

[0085] Fig. 9 This is a schematic diagram of the 3D results of the limited fault inversion of the 7.6 magnitude Mindanao earthquake in the first inversion case;

[0086] Fig.10 Schematic diagram of the slip plane results of the finite fault inversion of the 7.6 magnitude Mindanao earthquake in the first inversion case. DETAILED DESCRIPTION

[0087] In order to better explain the present invention and facilitate understanding, the present invention is described in detail below through specific implementation modes in conjunction with the accompanying drawings.

[0088] The multi-source observation data joint inversion method proposed in the embodiment of the present invention combines observation technologies such as InSAR (synthetic aperture radar interferometry), GNSS (global positioning system), and water level observation to provide richer, direct and indirect observation information. It can well reveal the spatiotemporal process of earthquake rupture, especially in remote areas or areas with sparse stations.

[0089] Observational data usually need to be used in conjunction with seismic waveform data. Although observational data provides additional information, it also brings technical challenges to data processing and fusion, such as data accuracy calibration and time synchronization.

[0090] By combining traditional seismological methods with the modern advantages of land observation technology and ocean observation technology, this inversion method fully captures the complexity of the earthquake rupture process. This method not only inherits the in-depth understanding of seismic wave propagation and earthquake source characteristics of traditional seismology, but also incorporates the high-precision surface deformation monitoring of land observation technology and the detailed observation of seabed seismic activity by ocean observation technology. Through this multidisciplinary and multi-technical integration, it is possible to examine the complex process of earthquake rupture from multiple dimensions and levels, and then construct a more detailed and accurate earthquake rupture model.

[0091] In order to better understand the above technical solution, 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 described herein. On the contrary, these embodiments are provided to enable a clearer and more thorough understanding of the present invention and to fully convey the scope of the present invention to those skilled in the art.

[0092] Embodiment 1

[0093] like Figure 1 As shown, this embodiment provides a flowchart of an earthquake rupture inversion method based on multi-source observation data. The method of this embodiment includes the following steps:

[0094] S100, obtaining a focal mechanism solution of a designated earthquake event;

[0095] S200, obtaining geometric parameters of the fault of the earthquake event according to the focal mechanism solution; for example, the geometric parameters of the fault of the designated earthquake event include: strike angle α, dip angle β, slip angle, etc.;

[0096] Estimate the focal mechanism solution of the earthquake event based on the historical earthquake event database at the same location, including: epicenter location, magnitude, strike angle, dip angle, slip angle, etc.;

[0097] In another embodiment, based on the seismic wave data of a designated earthquake event, a focal mechanism solution is obtained by fast calculation using W-phase.

[0098] S300: Generate a three-dimensional grid model of the fault plane to which the fault belongs according to the pre-constructed geometric shape model of the fault and the geometric parameters of the fault.

[0099] In this embodiment, the geometric shape model of the fault is constructed based on the data of Slab 2.0, a global comprehensive subduction zone geometric model library, and is modified based on the geometric parameters of the fault;

[0100] Specifically, the normal vector of the fault plane is calculated; the strike vector S is expressed by the strike angle α as:

[0101] S=(sinα,cosa,0)

[0102] The dip vector d is expressed by the fault dip angle β as: d = (cosβ, 0, -sinβ)

[0103] Then the normal vector n of the fault plane is: n = d × s = (cosαsinβ, -sinαsinβ, cosαcosβ)

[0104] Assume that there is a given point A on the fault plane: A(x0,y0,z0), and the coordinates of any point B on the fault plane are B(x,y,z), then the fault plane equation is:

[0105] That is, the three-dimensional grid model is z=z0+tanβ(x0-x)-tanαtanβ(y0-y),

[0106] Point A is the point selected as the epicenter.

[0107] S400: Obtain the slip amount of the sub-fault to which the earthquake rupture belongs according to the limiting conditions of each sub-fault in the three-dimensional grid model.

[0108] For example, assuming the dislocation plane and parameterizing the shape of the slip rate function, the slip inversion is transformed into the following formula (1) with double summation over the nth row and mth column:

[0109]

[0110] Among them, D jk is the slip amplitude, i.e. the slip amount of the sub-fault to be solved, λ jk is the sliding angle, i.e. the parameter of focal mechanism solution, S jk (t) is the given rise time function, u(t) is the total slip;

[0111] V jk is the average rupture velocity between the earthquake source and the sub-fault jk, V jk The limiting condition is that the maximum allowable rupture velocity is 80% of the deepest shear wave velocity in the crust velocity model of the fault Crust1.0;

[0112] is the sub-fault Green function. i is 1 or 2, t is the time, j and k represent the number of rows and columns of the sub-fault respectively, and n and m are natural numbers greater than 1.

[0113] In this embodiment, based on the focal mechanism solution and known parameters, the sub-fault Green function is obtained. The calculation formula of Green's function G(r,θ,z,t) is:

[0114]

[0115] Among them (U z ,U r ,U θ ) is the known surface harmonic vector basis coordinates, k is the known horizontal wave number set in advance, ω is the known circular frequency set in advance, r, θ, z are the cylindrical coordinate system, and t is the time.

[0116] S500: Perform inversion verification based on the observation data of the earthquake event and the slip amount of the sub-faults with the help of the inversion model. If the verification is successful, the slip amount of all sub-faults when the verification is successful is used as the earthquake rupture inversion result.

[0117] For example, based on the pre-given weight of each data type in the multi-source observation data, the data vector And according to the data vector Generate diagonal weight matrix According to the diagonal weight matrix Constructing the diagonal covariance matrix

[0118] Based on the following formula (3) and constraint condition (4), the least square method is used to solve until verification is passed;

[0119]

[0120] in, is the high-frequency GPS-GNSS data vector, is the seismic p-wave data vector, is the earthquake strong motion data vector; the norm symbol |||| is used to represent the norm of the vector,

[0121] is the high-frequency GPS data-GNSS diagonal covariance matrix, is the diagonal covariance matrix of seismic p-wave data, is the diagonal covariance matrix of earthquake strong motion data, G is the Green's function, m is the sub-fault slip, L is the Laplace operator, and β is the smoothing factor.

[0122] In this embodiment, before S500, the method further includes:

[0123] S001. Obtain multi-source observation data belonging to a specified earthquake event;

[0124] The multi-source observation data includes: InSAR data before and after the earthquake, GPS data before and after the earthquake, water level data corresponding to the earthquake event, and seismic wave data corresponding to the earthquake time;

[0125] Calibrate each observation data in the multi-source observation data, and achieve time synchronization and spatial registration to obtain pre-processed multi-source observation data;

[0126] Specifically, the data of seismic stations at a distance of 30°-90° from the epicenter of the specified earthquake event, i.e., acceleration, are selected.

[0127] Integrate the acceleration to velocity with a 100 to 5 Hz filter,

[0128] And perform bandpass filtering between 0.05 and 0.4 Hz, integrate the velocity waveform into the displacement waveform, and intercept the time window data of the specified time period from the arrival of the P wave as the pre-processed seismic wave data;

[0129] and / or, obtaining static GPS data of a first duration through a GPS station, wherein the first duration includes a time period from before the earthquake to after the earthquake,

[0130] The static GPS data were processed in precise point positioning mode PPP to obtain the time series of position coordinates. The co-seismic deformation was obtained by the difference between the average value of the time series of position coordinates a few days before the earthquake and the average value of the time series of position coordinates a few days after the earthquake.

[0131] and / or, obtaining original high-frequency GPS data of a second duration of the GNSS site, the second duration being less than or equal to the first duration;

[0132] Using the clock error and orbit information provided by the international GNSS service, the XYZ component displacement of the GNSS station is estimated, and high-frequency GPS data is calculated to obtain the time series of displacement at a sampling rate of 5 Hz;

[0133] The displacement time series is windowed, and the co-seismic displacement is calculated from the difference between the post-seismic position and the pre-seismic position in the windowed displacement time series. The co-seismic displacement is low-pass filtered at 0.5 Hz to obtain the final co-seismic displacement.

[0134] and / or, obtaining InSAR data, i.e., an interferogram spanning a long period of time before and after an earthquake, and processing the interferogram spanning a long period of time before and after an earthquake using software tools;

[0135] Use time series analysis method to analyze the interferogram, remove the atmospheric influence, and obtain the corrected interferogram;

[0136] For the corrected interferogram, 90m resolution SRTM DEM data is used to remove the terrain phase effect, and the interferogram without terrain phase effect is obtained; and the co-seismic deformation field is obtained by unwrapping processing; the co-seismic deformation field is downscaled and resampled using the quadtree method; and the pre-processed co-seismic deformation field is obtained;

[0137] and / or, obtaining water level data from buoys or tide stations and removing tidal information from the water level data.

[0138] The method of this embodiment realizes comprehensive monitoring of crustal deformation and movement state by fusing multi-source observation data, including InSAR, GPS, water level, seismic waveform and other data, and provides richer data support for the analysis of earthquake rupture process; a concise inversion algorithm is proposed, which can effectively process observation data and dynamically simulate earthquake rupture process, and can reflect the spatiotemporal evolution characteristics of earthquake rupture in real time, improve the stability and reliability of inversion results, and improve the accuracy of source parameters. It can provide a more reliable source for tsunami warning and disaster assessment. At the same time, it is also of great significance to improve the level of earthquake science research and enhance earthquake disaster defense capabilities, and has broad application prospects.

[0139] Embodiment 2

[0140] This embodiment provides a complete real-time process of an earthquake rupture inversion method based on multi-source observation data, which is as follows:

[0141] Step 200: Acquire multi-source observation data.

[0142] In this embodiment, relevant observation data such as InSAR, GPS, water level and seismic waveform are collected before and after the earthquake event.

[0143] In this embodiment, InSAR includes: pre-earthquake data within one month before the earthquake and data 5-10 days after the earthquake; GPS data includes: pre-earthquake data 3-15 days before the earthquake and post-earthquake data 1-15 days after the earthquake; water level data includes: data from 6 hours before the earthquake event to the end of the tsunami wave; earthquake waveform data includes: data from 30 minutes before the earthquake event to 30 minutes after the earthquake event.

[0144] The acquisition of geodetic observation data mainly relies on a variety of space technologies, including but not limited to seismic waveform data, satellite remote sensing, global navigation and positioning system (GNSS), synthetic aperture radar interferometry (InSAR), deep sea water pressure and buoy water level data, etc. These technologies can provide key data such as surface deformation and fault displacement before and after earthquakes.

[0145] Seismic wave data: Teleseismic body waves have high horizontal apparent velocity. Teleseismic body waves can provide accurate rupture occurrence time, and teleseismic surface waves can improve the resolution of seismic moments.

[0146] Satellite remote sensing: High-resolution remote sensing images can be used to monitor surface changes caused by earthquakes, such as surface cracks and landslides.

[0147] GNSS data: provides high-precision surface displacement data that can capture small changes before and after an earthquake, which is important for determining the amount of slip on the fault;

[0148] InSAR data: Through radar interferometry technology, the surface deformation field before and after the earthquake can be monitored, providing intuitive information on the spatial distribution of the earthquake rupture process;

[0149] Ocean observation data: Deep sea water pressure and buoy observations use deep-sea tsunami waves to better interpret the displacement of shallow layers in the subduction zone.

[0150] Step 201: preprocessing the multi-source observation data, that is, implementing data quality control to obtain preprocessed multi-source observation data.

[0151] After obtaining the original observation data, a series of calibration, time synchronization and spatial alignment are required to ensure the quality and reliability of the data.

[0152] Error analysis: Evaluate data uncertainty and error sources, including instrument errors, atmospheric effects, surface conditions, etc., to ensure the accuracy of the inversion results.

[0153] Data calibration: Calibrate the original data, remove system errors and noise such as instrument response through PZ method, atmospheric effect through time series analysis method, and terrain effect through DEM to improve data accuracy.

[0154] Time synchronization: Ensuring the temporal consistency of data from different sources is critical for dynamic monitoring of earthquake rupture processes.

[0155] Spatial registration: Spatial registration of data of different resolutions and sources, projecting all data into the same spatial coordinate system to facilitate subsequent data input and analysis.

[0156] Through careful processing and analysis of geodetic observation data, high-quality input data is provided for earthquake rupture inversion.

[0157] Different data have different advantages and disadvantages. Strong earthquake data are unreliable for a long time, while GPS data faithfully records low-frequency ground motions, has low sensitivity and low sampling rate, resulting in increased noise levels at higher frequencies. The single interferogram of InSAR data has a large amount of atmospheric noise, which partially masks the co-seismic signal and causes large errors in estimating co-seismic deformation. Different data have different magnitudes, units, and error ranges. When using any of the different types of data for slip inversion, they must be pre-processed in a pattern to maximize their strengths and avoid their weaknesses while ensuring the efficiency of the inversion calculation.

[0158] Seismic wave data (i.e., seismic waveform data): The time series information of the elastic wave field recorded by independent instruments (i.e., seismic wave data, which is acceleration) and the data of stations at a distance of 30°-90° from the epicenter are used to minimize the errors caused by the medium model in waveform inversion. The acceleration is integrated into velocity with a filter with a frequency of 100 to 5 Hz, and a bandpass filter between 0.05 and 0.4 Hz is performed. The velocity waveform is further integrated into the displacement waveform, and the time window of 30 to 100 s starting from the arrival of the P wave is selected as the teleseismic data (referring to the seismic wave data recorded by the more distant seismic stations).

[0159] GPS data: Static GPS data is processed in Precise Point Positioning (PPP) mode to obtain a time series of position coordinates. The co-seismic deformation is obtained by the difference between the average value of the time series a few days before the earthquake and the average value of the time series a few days after the earthquake.

[0160] High-frequency GPS data: Get the original high-frequency GPS data of the station (such as GNSS site), use the clock error and orbit information provided by the international GNSS service to estimate the three-component (XYZ) displacement of the GNSS site, obtain high-frequency GPS data, and calculate the time series of displacement at a sampling rate of 5Hz. The time length of high-frequency GPS data is shorter than the time length of the corresponding static GPS data.

[0161] The time series of displacement is windowed to avoid the influence of ground motion of the reference station. The co-seismic displacement is calculated by taking the difference between the post-seismic position and the pre-seismic position, and then the co-seismic displacement is low-pass filtered at 0.5 Hz to obtain the co-seismic displacement for subsequent calculations. The pre-processing of this embodiment can assign certain parameters to the deformation data in the vertical direction and reduce its inversion weight to eliminate the higher noise level in this direction of motion.

[0162] InSAR data: InSAR data, i.e., interferograms spanning a long period before and after the earthquake, are obtained. Gamma, MintPy and other software tools are used to process a series of interferograms spanning a long period before and after the earthquake. At the same time, the time series analysis method is used to improve the coseismic information mainly from atmospheric masking through the PyAPS atmospheric correction of the European Center for Medium-Term Weather Forecasts Reanalysis (ERA5) weather model; 90m resolution SRTM DEM data is further used to remove the terrain phase effect, and finally the interferogram is unwrapped to obtain the coseismic deformation field. Before inversion, the coseismic deformation field is downscaled and resampled using the quadtree method to obtain the preprocessed coseismic deformation field to ensure the efficiency of the inversion calculation; each sampling point is regarded as a separate station, similar to the GPS observation value of static offset.

[0163] Water level data: Tidal information was removed from the recorded water level data, which recorded short-term changes in sea surface height, and resampled at a sampling rate of 15 seconds. Tide gauge data existed near the source area and were used for static inversion. The background noise level of the sensor was calculated using 6 hours of pre-event noise, while the waveform envelope was estimated using the Hilbert transform, and the duration of the tsunami was estimated by measuring the time required for the record to decay to the pre-event state through comparative analysis.

[0164] Step 202: construct a fault model.

[0165] Through the subdivision of the three-dimensional fault geometry, the entire structure is divided into multiple sub-faults to more accurately describe the slip process. At the same time, with reference to the global comprehensive subduction zone geometry model library Slab 2.0, this data constructs a detailed fault geometry model.

[0166] While obtaining information on the focal mechanism, one can also refer to historical stratigraphic research results published by other scholars and research reports and information published by relevant geological or geoscience institutions to ensure the accuracy of the geometric structure.

[0167] When the geometric parameters of this fault (focal mechanism solution) are known, including the strike angle α and dip angle β, such as Figure 2 As shown, these parameters are used to generate a three-dimensional grid model of the fault plane.

[0168] First, calculate the normal vector of the fault plane. The strike vector S is expressed by the strike angle α as:

[0169] S=(sinα,cosa,0)

[0170] The dip vector d is expressed by the fault dip angle β as:

[0171] d = (cosβ, 0, -sinβ)

[0172] Then the normal vector n of the fault plane is:

[0173] n=d×s=(cosαsinβ,-sinαsinβ,cosαcosβ)

[0174] Assuming there is a given point A (select the epicenter position) A(x0, y0, z0) on the fault plane, and the coordinates of any point B on the fault plane are B(x, y, z), then the fault plane equation can be expressed as:

[0175] That is, the three-dimensional grid model z: z = z0 + tanβ (x0-x) - tanαtanβ (y0-y);

[0176] Each small grid is called a sub-fault of the large fault plane. The response of each small sub-fault to the earthquake finally constitutes the rupture surface of the fault plane.

[0177] Step 203: Joint inversion.

[0178] If the time, magnitude, direction and geographical location of an earthquake are unknown, then inverting the temporal and spatial variation characteristics of the earthquake source is a complex problem.

[0179] A. Sub-fault response.

[0180] Based on the rupture kinematics theory, the slip inversion is transformed into a steady-state problem by assuming the dislocation plane and parameterizing the shape of the slip rate function. The four parameters (slip amount, slip direction, rise time function, and rupture velocity function) used to describe the response of the sub-fault are used. The motion (displacement) can be obtained by double summing the nth row and the mth column:

[0181]

[0182] Among them, D jk is the slip amplitude (offset, slip of the sub-fault to be solved), λ jk is the sliding angle (a parameter in the focal mechanism solution), S jk (t) is a given rise time function, V jk is the average rupture velocity between the earthquake source and the sub-fault jk, is the sub-fault Green function.

[0183] If the rise time is further assumed, the inversion becomes a linear problem. The linear multi-time window method is used to allow the sub-fault to slip in multiple intervals. This method can adapt to different rupture velocities and can produce complex source time functions.

[0184] These sub-fault slip rate functions allow slip to occur in multiple overlapping regions, with an overlap of 50% in each region. The rise time of each slip rate function lasts for several seconds, ensuring that the dynamic changes of the slip process can be accurately captured. The rise time is based on the expected value of the magnitude scaling law obtained from earthquake events occurring worldwide.

[0185] log(s)=-0.5323+0.293log 10(M)

[0186] Where: s is the average rise time, M is the magnitude. In this way, the reliability and accuracy of the model in practical applications are ensured.

[0187] In addition, multiple possible rupture velocity tests were performed to evaluate the fault rupture behavior under different scenarios. In these tests, the maximum allowable rupture velocity was set to 80% of the deepest shear wave velocity in the crust velocity model of the fault Crust1.0.

[0188] B. Green's function calculation.

[0189] A method for solving the point source elastodynamic and elastostatic problems in multilayer half-spaces using the frequency-wavenumber method of Zhu and Rivera. Using the propagator matrices in the dynamic and static cases, when the dynamic propagator matrices and solutions converge to their static frequency ω→0, the dynamic solution close to zero frequency is a satisfactory static deformation.

[0190] In the calculation process, the Fourier-Hankel transformation is first used to transform the wave equation in the cylindrical coordinate system (r, θ, z) in the time-space domain into a wave equation in the frequency-wavenumber domain. The equation is simplified and converted into a matrix ordinary differential equation about the depth z. According to the stress and displacement boundary conditions, the coefficient matrices of the P-SV and SH systems are solved respectively by the Thompson-Haskell transfer coefficient matrix method to obtain the response caused by the earthquake source. Then, the response in the time-space domain (Ur, Uz, Uz) is obtained by the Fourier-Hankel inverse transformation. θ ) T , that is, the surface harmonic vector displacement solution, the final Green's function expression in the time-space domain is:

[0191]

[0192] Among them, (U z ,U r ,U θ ) is the surface harmonic vector basis coordinates, k is the horizontal wave number, and ω is the circular frequency. In the calculation, the Thompson-Haskell transfer coefficient matrix was decomposed into Jordan standard form to simplify the frequency-containing expression. The decoupling of the wave number and frequency calculation process and the integral operation in the wave number domain greatly improve the calculation speed. The calculation of the Green's function is independent of the source description. The crustal velocity structure model expresses the changes of the density, velocity, quality factor and other characteristics of the crustal medium with depth, controls the amplitude, frequency component and duration of the Green's function, and expresses the attenuation characteristics of the earthquake motion. Therefore, after the existing regional velocity structure model is determined, the Green's function can be calculated.

[0193] The Green's function of each sub-fault and the static and elastic dynamic Green's functions of InSAR, GPS, GNSS and SM data are calculated in the range of 0 to 0.5 Hz.

[0194] The Green's function of the water level data is calculated by setting the seafloor deformation field of the high-resolution grid within the model range obtained by the slip calculation on each sub-fault and inputting this deformation field into the Geoclaw tsunami model.

[0195] C. Inversion model.

[0196] For joint kinematic inversion using different geophysical data types, defining the relative weight of each type of data is an ambiguous process. First, by normalizing the vector norm, each data type is given equal importance. Consider a data vector consisting of InSAR, SM, tidal measurement, and HR-GNSS data vectors.

[0197] Before inversion, a diagonal weight matrix is ​​defined

[0198]

[0199] Among them, ||.|| is the norm of the vector of each data type. This method completely excludes any preferentially suitable data type and ensures that no specific data type will be given priority in the selection process. At the same time, in order to elaborate on the "credibility" of each measurement value, the concept of diagonal covariance matrix is ​​further introduced. The diagonal covariance matrix is ​​a special matrix whose non-diagonal elements are all zero, and the elements on the diagonal represent the variance of each measurement value. In this way, the uncertainty of each measurement value can be quantified and its reliability can be evaluated. Specifically, the diagonal covariance matrix Each diagonal element in represents the variance of a measurement value. The smaller the variance, the higher the credibility of the measurement value, and vice versa. This definition can more accurately analyze and compare the credibility of different measurement values, thereby providing a more reliable basis in data analysis and decision-making.

[0200]

[0201] where the standard deviation of the pre-event noise for each record is used to create and Similarly, It is also obtained from the standard deviation of the pre-event portion of the tide gauge record. According to the estimation from the far-field noise, where the pixels are correlated as a function of their distance from each other. During the inversion process, two weight matrices are multiplied on both sides of the equation to form the Green function matrix G,

[0202]

[0203] The slip vector at each sub-fault is represented by the variable m. In order to solve this problem more accurately, the least square method is used.

[0204]

[0205] in, is the deformation observation; is the weight matrix of the observations, G is the Green's function, m is the sub-fault slip, L is the Laplace operator, and β is the smoothing factor. This method can effectively deal with the noise and uncertainty in the data by introducing a constraint condition, thereby improving the stability and reliability of the solution. In this way, the slip vector at each sub-fault can be better estimated, and the slip process of the entire earthquake fault can be more fully understood.

[0206] Embodiment 3

[0207] The fault rupture of the 2024 Mindanao earthquake with a magnitude of 7.6 is used as an inversion case.

[0208] Step 1: Observation data and processing.

[0209] The original three-component P-wave data of the nine seismic stations used for calculation and inversion are calibrated and quality controlled as follows: Figure 4A , Figure 4B and Figure 4C shown.

[0210] The waveform was trimmed to 75 seconds in length and downsampled to 5 Hz. The velocity was integrated into the time series of displacement changes and bandpass filtered between 0.005 and 0.2 Hz for inversion. The results are shown in Figure 5A , Figure 5B and Figure 5C shown.

[0211] Step 2: Joint inversion.

[0212] A. Source model of the Mindanao earthquake.

[0213] After the Mindanao earthquake, referring to the focal mechanism parameters (epicenter, magnitude, strike angle, strike angle, etc.) published by multiple agencies including the United States Geological Survey (USGS), the USGS measurement results were selected (epicenter 126.449E, 8.527N, focal depth 32.8km; strike 167°, dip 17°, slip angle 63°), and referring to the expected values ​​of the magnitude scaling law during rise obtained from earthquake events occurring worldwide. Given an average rise time of 2s, the rupture velocity range is 2.0 to 3.0km / s, the maximum allowable rupture velocity is equivalent to 80% of the deepest layered shear wave velocity spanned by the fault model, and the possible slip angle range is between 60° and 66°.

[0214] B. Tectonic fault model.

[0215] At the same time, the fault model was constructed by referring to the source parameters of multiple institutions and databases such as Slab2.0, and 361 sub-faults were divided, such as Figure 6 shown.

[0216] C. Crustal velocity structure model.

[0217] The determination of a detailed three-dimensional velocity structure model requires the comprehensive use of seismic exploration results, regional geological data, drilling data, and formation wave velocity test data. The underground undulating structure and topographic factors are not considered here, and it is assumed that the wave velocity, density, and damping ratio of the same horizontal rock layer in the study area are constant.

[0218] The one-dimensional crustal velocity structure information of the Mindanao region is constructed with reference to the Crust1.0 model (as shown in the table below). This model describes the main stratification of the crust, and the area below the Moho surface is simplified to half space. The frequency-wavenumber method uses a one-dimensional crustal model that matches the current level of understanding of the crustal structure, which expresses the stratified characteristics of the crustal medium and ensures computational efficiency.

[0219] Crustal velocity structure model in Mindanao

[0220]

[0221] D. Calculate the Green's function.

[0222] The Green's function is calculated using the frequency-wavenumber method (FK). The Fourier-Hankel transform is used to transform the wave equation in the cylindrical coordinate system (r, θ, z) in the time-space domain into a wave equation in the frequency-wavenumber domain. The equation is simplified and converted into a matrix ordinary differential equation about the depth z. According to the stress and displacement boundary conditions, the coefficient matrices of the P-SV and SH systems are solved respectively by the Thompson-Haskell transfer coefficient matrix method to obtain the response caused by the earthquake source. The response in the time-space domain is then obtained by the Fourier-Hankel inverse transform, and finally the Green's function in the time-space domain is obtained. Figure 7 shown.

[0223] E. Inversion results.

[0224] The optimal parameters are determined by minimizing the residual between the observed and simulated values. Overall, the inversion results show that the observed seismic motions are consistent with the simulated seismic motions. The details are better reflected; in the inversion of the fault slip surface, the focal center area is well simulated; the slip on the fault plane is mainly distributed in the range of 0km to -30km from the epicenter and -20km to -60km in dip. The fault mainly slips in the southeast direction near the epicenter. The maximum slip of this earthquake is located between the focal center and the northeast direction of about 20 to 40km; in terms of geographical coordinates, it is located at 127° west longitude, 9.25° north latitude, and a depth of 20km. The maximum slip is 1.5m. Figure 8 , Fig. 9 and Fig.10 shown.

[0225] The present invention relates to a method for joint inversion of earthquake rupture process by combining multi-source observation data, which belongs to the intersection of seismology, geodesy and oceanography. The method dynamically simulates the earthquake source rupture process by integrating InSAR, GPS, water level and other observation data with seismic wave data, thereby improving the accuracy and efficiency of earthquake rupture process inversion.

[0226] The method of the present invention integrates multiple observation data, including but not limited to InSAR, GPS and seismic waveform data, to overcome the limitation of a single data source. The method of the present invention provides a parameter optimization method for the joint inversion of earthquake rupture process using multi-source data, thereby improving the accuracy and speed of inversion model parameter estimation. Based on the constraint optimization method, the stability of the inversion model and the rationality of the solution are enhanced.

[0227] The method of the present invention effectively simplifies the dynamic rupture problem into a linear problem. First, the assumption of uniform source propagation velocity is not applicable to the rupture mode of all earthquakes. Therefore, in the inversion process, the rupture velocity is adjusted by fitting and comparing the inversion results with the data. Secondly, multi-time window linearization usually requires more parameters to simulate the rupture process, and the inversion results are stabilized by introducing stronger constraints. At the same time, the method uses two parameters to describe the source time function (displacement, sliding direction) of the sub-fault, thereby reducing the number of inversion parameters.

[0228] In joint inversion, the weights and relative errors between different data are the core challenges. Since different data have different magnitudes, units, and error ranges, it is challenging to integrate them into an inversion system. We first normalize the different data and then perform the inversion. In order to consider the credibility of the information between different data, we use the covariance matrix of the data for weighting, so that the relative weight between the data is inversely proportional to the size of each error, thereby increasing the weight of the credible data. The weights of each observation are assigned using different attributes of the observations and repeated testing methods.

[0229] In addition, an embodiment of the present invention further provides a computing device, comprising: a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program in the memory to specifically perform the steps of a seismic rupture inversion method based on multi-source observation data described in any of the above embodiments.

[0230] It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as methods, systems or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0231] The present invention is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems) and computer program products according to the embodiments of the present invention. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions.

[0232] It should be noted that in the claims, any reference numerals placed between brackets shall not be construed as limiting the claims. The word "comprising" does not exclude the presence of components or steps not listed in the claims. The word "a" or "an" preceding a component does not exclude the presence of a plurality of such components. The invention may be implemented by means of hardware comprising several different components and by means of a suitably programmed computer. In the claims enumerating several means, several of these means may be embodied by the same hardware. The use of the words first, second, third, etc., is for convenience of expression only and does not indicate any order. These words may be understood as part of the component name.

[0233] In addition, it should be noted that, in the description of this specification, the description of the terms "one embodiment", "some embodiments", "embodiment", "example", "specific example" or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, unless they are contradictory.

[0234] Although the preferred embodiments of the present invention have been described, those skilled in the art may make other changes and modifications to these embodiments after knowing the basic creative concept. Therefore, the claims should be interpreted as including the preferred embodiments and all changes and modifications falling within the scope of the present invention.

[0235] Obviously, those skilled in the art can make various modifications and variations to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention should also include these modifications and variations.

Claims

1. A method for earthquake rupture inversion based on multi-source observation data, characterized in that: include: S100, obtaining a focal mechanism solution of a designated earthquake event; S200, obtaining geometric parameters of the fault of the earthquake event according to the focal mechanism solution; S300, generating a three-dimensional grid model of the fault plane to which the fault belongs according to the pre-constructed geometric shape model of the fault and the geometric parameters of the fault; S400, obtaining the slip amount of the sub-fault to which the earthquake rupture belongs according to the limiting conditions of each sub-fault in the three-dimensional grid model; S500: Perform inversion verification based on the observation data of the earthquake event and the slip amount of the sub-faults with the help of the inversion model. If the verification is successful, the slip amount of all sub-faults when the verification is successful is used as the earthquake rupture inversion result.

2. The method according to claim 1, characterized in that The geometric parameters of the fault that specifies the earthquake event include: strike angle α, dip angle β, and slip angle; The focal mechanism solution of the earthquake event is estimated based on the database of historical earthquake events at the same location, including: epicenter location, magnitude, strike angle, dip angle, and slip angle.

3. The method according to claim 1, characterized in that The S300 includes: The fault geometry model is constructed based on the data of Slab 2.0, a global comprehensive subduction zone geometry model library, and is modified based on the geometric parameters of the fault. Specifically, the normal vector of the fault plane is calculated; the strike vector S is expressed by the strike angle α as: S=(sinα,cosa,0) The dip vector d is expressed by the fault dip angle β as: d = (cosβ, 0, -sinβ) Then the normal vector n of the fault plane is: n = d × s = (cosαsinβ, -sinαsinβ, cosαcosβ) Assume that there is a given point A on the fault plane: A(x0,y0,z0), and the coordinates of any point B on the fault plane are B(x,y,z), then the fault plane equation is: That is, the three-dimensional grid model is z=z0+tanβ(x0-x)-tanαtanβ(y0-y), Point A is the point selected as the epicenter.

4. The method according to claim 3, characterized in that The S400 includes: Assuming the dislocation plane and parameterizing the shape of the slip rate function, the slip inversion is transformed into the following formula (1) with double summation over the nth row and mth column: Among them, D jk is the slip amplitude, i.e. the slip amount of the sub-fault to be solved, λ jk is the sliding angle, i.e. the parameter of focal mechanism solution, S jk (t) is a given rise time function, V jk is the average rupture velocity between the earthquake source and the sub-fault jk, V jk The limiting condition is that the maximum allowable rupture velocity is 80% of the deepest shear wave velocity in the crust velocity model of the fault Crust1.0; is the sub-fault Green function; i is 1 or 2, t is time, j and k represent the row and column numbers of the sub-fault respectively; u(t) is the total slip.

5. The method according to claim 4, characterized in that Based on the focal mechanism solution and known parameters, obtain the sub-fault Green function The calculation formula of Green's function G(r,θ,z,t) is: Among them (U z ,U r ,U θ ) is the known surface harmonic vector basis coordinates, k is the known horizontal wave number set in advance, ω is the known circular frequency set in advance, r, θ, z are the cylindrical coordinate system, and t is the time.

6. The method according to claim 4, characterized in that Before S500, the method further includes: Get multi-source observation data belonging to the specified earthquake event; The multi-source observation data includes: InSAR data before and after the earthquake, GPS data before and after the earthquake, water level data corresponding to the earthquake event, and seismic wave data corresponding to the earthquake time; Calibrate each observation data in the multi-source observation data, and achieve time synchronization and spatial registration to obtain pre-processed multi-source observation data; Accordingly, S500 includes: According to the preprocessed observation data and the slip amount of the sub-faults, the inversion verification is carried out with the help of the inversion model. If the verification passes, the slip amount of all sub-faults when the verification passes will be used as the earthquake rupture inversion result.

7. The method according to claim 6, characterized in that The S500 includes: Based on the pre-given weight of each data type in the multi-source observation data, a data vector is formed And according to the data vector Generate diagonal weight matrix According to the diagonal weight matrix Constructing the diagonal covariance matrix Based on the following formula (3) and constraint condition (4), the least square method is used to solve until verification is passed; in, is the high-frequency GPS-GNSS data vector, is the seismic p-wave data vector, is the earthquake strong motion data vector; the norm symbol |||| is used to represent the norm of the vector, is the high-frequency GPS data-GNSS diagonal covariance matrix, is the diagonal covariance matrix of seismic p-wave data, is the diagonal covariance matrix of earthquake strong motion data, G is the Green's function, m is the sub-fault slip, L is the Laplace operator, and β is the smoothing factor.

8. The method according to claim 6, characterized in that Acquire multi-source observation data belonging to a specified earthquake event, calibrate each observation data in the multi-source observation data, and achieve time synchronization and spatial registration to obtain pre-processed multi-source observation data, including: Select the data of seismic stations at 30°-90° from the epicenter of the specified earthquake event, i.e. acceleration. Integrate the acceleration to velocity with a 100 to 5 Hz filter, A bandpass filter between 0.05 and 0.4 Hz is performed, the velocity waveform is integrated into the displacement waveform, and the time window data of the specified time period from the arrival of the P wave is intercepted as the preprocessed seismic wave data and / or, Obtain static GPS data of the first duration through GPS stations, where the first duration includes the time period from before the earthquake to after the earthquake; The time series of position coordinates is obtained by PPP processing of static GPS data in precise point positioning mode. The co-seismic deformation is obtained by the difference between the average value of the position coordinates a few days before the earthquake and the average value of the position coordinates a few days after the earthquake.

9. The method according to claim 6, characterized in that Acquire multi-source observation data belonging to a designated earthquake event, calibrate each observation data in the multi-source observation data, and achieve time synchronization and spatial registration to obtain pre-processed multi-source observation data, and also include: Acquire raw high-frequency GPS data of a second duration of the GNSS site, where the second duration is less than or equal to the first duration; Using the clock error and orbit information provided by the international GNSS service, the XYZ component displacement of the GNSS station is estimated, and high-frequency GPS data is calculated to obtain the time series of displacement at a sampling rate of 5 Hz; The displacement time series is windowed, and the co-seismic displacement is calculated from the difference between the post-seismic position and the pre-seismic position in the windowed displacement time series. The co-seismic displacement is low-pass filtered at 0.5 Hz to obtain the final co-seismic displacement. and / or, Acquire InSAR data, i.e., interferograms spanning a long period before and after an earthquake, and use software tools to process the interferograms spanning a long period before and after an earthquake; Use time series analysis method to analyze the interferogram, remove the atmospheric influence, and obtain the corrected interferogram; For the corrected interferogram, 90m resolution SRTM DEM data is used to remove the terrain phase effect, and the interferogram without terrain phase effect is obtained; and the co-seismic deformation field is obtained by unwrapping. The co-seismic deformation field is downscaled and resampled using the quadtree method to obtain the pre-processed co-seismic deformation field; and / or, obtaining water level data from buoys or tide stations and removing tidal information from the water level data.

10. A computing device, characterized in that include: A memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program in the memory to specifically perform the steps of the earthquake rupture inversion method based on multi-source observation data as described in any one of claims 1 to 9.

Images