A method and system for inverting source rupture process based on multi-source observation data

Abstract

The present application relates to the technical field of seismic data processing, and more particularly to a source rupture process inversion method and system based on multi-source observation data; the present application jointly uses multi-source observation data of earthquakes and geodetic survey, combines multi-dimensional advantage resolution capability, and obtains more accurate inversion results; meanwhile, the present application constructs a fault geometry model adaptive to different complexity sources, which is more suitable for the real three-dimensional form of the seismogenic fault; in addition, the present application proposes a space-time scale adaptive smoothing constraint based on the driving of the rupture physical process, realizes the synchronous regulation of the space-time constraint strength by a single hyperparameter, and obviously reduces the difficulty of hyperparameter selection; the present application also realizes the compatibility of at least two fault motion mechanisms by the non-negative constraint of the half-space inequality of the main sliding direction, while prohibiting non-physical reverse sliding, and improves the physical rationality of the inversion model.

Description

Technical Field

[0001] This invention relates to the field of earthquake data processing technology, and in particular to a method and system for inverting the source rupture process based on multi-source observation data. Background Technology

[0002] For large earthquakes with a moment magnitude exceeding 6, the focal rupture scale can reach tens of kilometers or more. Reconstructing the focal rupture process is the core approach to accurately obtain detailed rupture characteristics and a key means to deeply understand the causes and mechanisms of earthquakes. The focal rupture process characterizes the spatiotemporal evolution of the slippage of the seismogenic fault due to rupture, which is of great significance for understanding the earthquake mechanism and assessing earthquake damage.

[0003] Using observational data to infer the spatiotemporal distribution of fault slip is the primary method for establishing models of earthquake source rupture processes. Modern earthquake observations, including seismograph and strong-motion meter records of ground motion waveforms, can be used to invert source rupture processes. With the rapid development of geodetic observation technologies based on satellite systems, such as Synthetic Aperture Radar (SAR) and Global Navigation Satellite System (GNSS), surface displacement data acquired through these observations can also be used to invert fault slip distributions or source rupture processes. Typical source rupture process inversions use only one or more sets of data of a certain type from seismic or geodetic observation data. Because different observation data originate from different observation modes, inversions using only one type of data inevitably suffer from bias, leading to significant uncertainties in the inversion results and limiting the ability to resolve source information.

[0004] The aspect ratio of the fault plane of a large earthquake is generally greater than 1. Especially for terrestrial earthquakes, the width of the ruptured fault is limited to tens of kilometers by the thickness of the seismogenic layer, while the rupture length along the strike can reach hundreds of kilometers or more. The larger rupture scale makes the geometric characteristics of the fault along the strike more complex, which is even more pronounced in cascading large earthquakes with continuous multi-fault ruptures. However, conventional source rupture process inversion uses single or multiple planar faults to parameterize the source model, resulting in significant differences between the source model characteristics and the actual fault geometry, leading to inaccurate source rupture process inversion results.

[0005] Furthermore, existing smoothing constraints in earthquake source rupture process inversion have shortcomings. Traditional Laplace smoothing constraints are purely mathematical, equal-weighted difference forms that do not fully incorporate the physical laws of earthquake rupture. Spatial and temporal constraints require separate hyperparameter settings, leading to highly subjective and computationally expensive parameter tuning processes, and making them prone to imbalances in spatiotemporal constraint strength. Simultaneously, the non-negative constraints used in conventional inversion strictly limit the fault slip direction to a narrow range on both sides of the main slip direction, forcing the entire fault to exhibit only a single fault motion mechanism. This makes it difficult to adapt to the complex scenarios of multiple source mechanisms coexisting in terrestrial earthquakes, resulting in systematic deviations between the inversion results and the actual earthquake occurrence process. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention provides a method and system for inverting the source rupture process based on multi-source observation data, which solves the problems existing in the prior art.

[0007] This invention provides a method for inverting the source rupture process based on multi-source observation data, the method comprising the following steps: S100: Acquire seismic and geodetic observation data from multiple sources; S200: Construct an observation value vector based on the earthquake observation data and the geodetic observation data; S300: Establish a fault geometry model, adaptively discretize the fault plane, and obtain multiple rectangular sub-fault dislocation elements; S400: Based on the seismic observation vector and the fault geometric model, calculate the Green's function and establish the joint inversion equation; S500: Solve the joint inversion equations to obtain the inverted source rupture process model; S500 specifically refers to: S510: Apply spatial smoothing constraints to the inversion of the source rupture process to make the slip or seismic moment release of adjacent sub-faults equal. At the same time, introduce temporal smoothing constraints to make the slip or seismic moment release of the same sub-fault equal in adjacent time windows. S520: In the inversion of the source rupture process, the reverse slip of the fault is not considered. The slip of each sub-fault is set to vary freely in half-space in the main slip direction. By solving the joint inversion equation with non-negative inequality constraints, the spatiotemporal distribution of fault slip characterizing the source rupture process is obtained. S530: The Jackknife test is used to quantitatively obtain the uncertainty of fault slip distribution, reflecting the impact of observational data on the inversion results of the source rupture process.

[0008] Preferably, in S100, the earthquake observation data includes: strong motion three-component waveform recording data obtained from strong motion meter stations near the epicenter, vertical component P-wave recording data and tangential component SH-wave recording data obtained from far-field seismograph stations near the epicenter; the geodetic observation data includes: surface deformation field data obtained from synthetic aperture radar satellite interferometry covering the earthquake rupture range, and surface displacement data obtained from global navigation satellite system stations near the epicenter.

[0009] Preferably, in step S200, different datasets are processed into vector form and combined together to form an observation vector, specifically expressed as follows: , in, Indicates the first A vector of observed values ​​from an observation dataset. Indicates the first The first observation dataset There are [number] observations, and the total number of observation datasets is [number]. The number of observations in each dataset is , Observation dataset Together they form the observation vector.

[0010] Preferably, S300 specifically comprises: S310: Use planar Cartesian coordinates that are continuously distributed along the fault strike to describe the location of the fault crest edge, i.e., the fault trace, and specify the fault dip angle for each trace point; S320: The fault trace is discretized using the principle of linearization, and the coordinates and dip parameters of the control points of the fault trace are calculated by piecewise linear interpolation; S330: To determine the coordinates of the remaining sub-fault control points, first, the center coordinates of the fault trace control points need to be calculated. Then, the arithmetic mean of the plane coordinates of the fault trace control points is calculated, along with the average strike. Finally, based on the average strike and a specified dip angle, the three-dimensional coordinates of each sub-fault control point are calculated. S340: Traverse the sub-fault mesh along the strike and dip directions of the fault, and calculate the coordinates of the calculation center, strike, normal vector, area, dip angle, and size parameters of each sub-fault.

[0011] Preferably, S320 specifically includes: S321: Determine the number of rectangular sub-faults along the strike and dip of the fault; S322: First, a zero matrix is ​​constructed to store the three-dimensional coordinates of the control points of each sub-fault's spatial location, and a zero vector is constructed to store the dip angle parameters that vary along the strike. The coordinates of the top control point of the shallowest sub-fault are used as the fault trace control points. For the intermediate discrete points among the fault trace control points, their depths are first fixed. ,in, For a fixed depth value, traverse each control point segment, and then determine whether the cumulative distance along the direction of the current discrete point falls within the range of the first segment. If the cumulative length interval of the description points satisfies the condition, then the interpolation coefficients are calculated. ;based on For the first and the Linear interpolation is performed on the plane coordinates and dip angle of each description point to obtain the plane coordinates and dip angle of the current fault trace control point. After completing the interpolation calculation of the current control point, the traversal of the control point segment is terminated, and the next fault trace control point is processed.

[0012] Preferably, in step S321, the number of rectangular sub-faults along the fault strike and dip is determined, for a length of... Width is Sub-faults, strike direction, number of sub-faults The total length of the fault trace The length of the sub-fault is Calculation, through Ensure that at least two sub-faults are defined to maintain geometric rationality. , Number of dip direction sub-faults Based on the fault dip width with sub-fault width Similarly, the minimum constraint value is 2. , A slight correction was made to the spatial distribution step size of the sub-faults: , .

[0013] Preferably, in step S322, the discretization of the fault trace is achieved through linear interpolation, first constructing a dimension of The zero matrix stores the three-dimensional coordinates of the control points of each sub-fault's spatial location. ,in East-west coordinates, North-south coordinates For depth coordinates, construct a structure with length... The zero vector stores the tilt angle parameters that vary along the azimuth. The coordinates of the top control point of the shallowest sub-fault are used as the fault trace control points, where the first and last points are directly taken from the coordinates of the fault trace description points. and specify the depth of the fault crest as , , , , , , , , , For intermediate discrete points in the fault trace control points First, fix its depth. Iterate through each control point segment and then determine the cumulative distance along the direction of the current discrete point. Whether it falls into the first The cumulative length interval of each description point If the conditions are met, then calculate the interpolation coefficients. : , based on For the first and the Linear interpolation is performed on the planar coordinates and inclination angles of each description point to obtain the planar coordinates and inclination angles of the current trace control point. , , , After completing the interpolation calculation for the current control point, terminate the traversal of the control point segments and continue processing the next fault trace control point.

[0014] Preferably, S400 specifically comprises: S410: Establish linear equations for inverting the source rupture process based on the source representation theorem; S420: Calculate the Green's function based on the three-dimensional coordinates of each sub-fault in the fault geometry model and the geographic coordinates of seismic observations; S430: Construct joint inversion equations for multi-source data using inversion equations for various types of observation data.

[0015] According to another aspect of the present invention, a source rupture process inversion system based on multi-source observation data is provided. The system employs the aforementioned source rupture process inversion method based on multi-source observation data. The system includes: The observation data acquisition module is used to acquire seismic observation data and geodetic observation data from various sources. The multi-source data processing module is used to process the earthquake observation data and geodetic observation data to construct an earthquake observation value vector; The fault geometry generation module is used to establish the fault geometry model and perform adaptive discretization sampling to obtain multiple rectangular sub-fault dislocation elements. The inversion equation construction module is used to calculate the Green's function and establish joint inversion equations based on the seismic observation vector and fault geometry model. The joint inversion solution module is used to solve the joint inversion equations to obtain the inverted source rupture process model.

[0016] The embodiments of the present invention have the following technical effects: The inversion method proposed in this invention combines multi-source observation data from seismic and geodetic surveys. The strength of the observation data's ability to resolve source information is reflected in three dimensions: spatial resolution, temporal resolution, and slip / seismic moment resolution. The superior resolution of observation data is usually concentrated in one or two of the above three resolution dimensions. The combined use of multi-source data combines the superior resolution of multiple dimensions, achieving the best resolution effect allowed by each type of observation data in the spatial, temporal, and slip / seismic moment dimensions, thereby obtaining more accurate inversion results. Meanwhile, the inversion method of this invention constructs a complex fault geometry model based on the fault trace and a specified dip angle and adaptively achieves discrete parameterization. By continuously changing the strike and dip angle in space, it reflects the complex geometric features of the fault plane, which is closer to the true three-dimensional shape of the seismogenic fault. It can also maintain the spatiotemporal continuity of the rupture propagation process and enhance the reliability of the inversion results.

[0017] Meanwhile, the temporal and spatial smoothing constraints used in the inversion method of this invention link the sub-fault rupture scale, the time window rupture scale, and the overall rupture scale of the seismic source. By utilizing the relative proportional relationship, the smoothing constraint forms in the temporal and spatial domains are effectively unified, enabling the strength of temporal and spatial smoothing constraints to be controlled simultaneously using only one hyperparameter, which significantly reduces the difficulty of selecting hyperparameters.

[0018] Meanwhile, the inversion method of the present invention, by applying non-negative constraints of the half-space inequality in the main slip direction, allows the sub-fault slip to change freely in the half-space of the main slip direction, which can naturally accommodate at least two types of fault motion mechanisms. This makes the inversion model no longer limited to the preset source mechanism, effectively restores the complex source mechanism of the real seismogenic fault, and significantly improves the physical rationality of the inversion model. Attached Figure Description

[0019] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0020] Figure 1This is a schematic diagram of the source rupture process inversion method based on multi-source observation data provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the specific process of S300 provided in an embodiment of the present invention; Figure 3 A schematic diagram of a fault geometric model established using the Mw7.8 and Mw7.7 double earthquakes in a certain area in 2023 as an example, provided for an embodiment of the present invention; Figure 4 A schematic diagram of spatial and temporal smoothing constraints provided for embodiments of the present invention; Figure 5 A schematic diagram of InSAR data fitting for an embodiment of the present invention, using the Mw7.8 and Mw7.7 double earthquakes in a certain area in 2023 as an example; Figure 6 A schematic diagram of GNSS horizontal data fitting provided for an embodiment of the present invention, taking the double earthquakes of Mw7.8 and Mw7.7 in a certain area in 2023 as an example; Figure 7 A schematic diagram of fitting strong ground motion three-component waveform data, taking the Mw7.8 earthquake in a certain area in 2023 as an example, as provided in an embodiment of the present invention; Figure 8 A schematic diagram of fitting strong ground motion three-component waveform data, taking the Mw7.7 earthquake in a certain area in 2023 as an example, as provided in an embodiment of the present invention; Figure 9 A schematic diagram of far-field seismometric P-wave data fitting, using the Mw7.8 earthquake in a certain area in 2023 as an example, provided as an embodiment of the present invention; Figure 10 A schematic diagram of far-field seismometric P-wave data fitting, using the Mw7.7 earthquake in a certain area in 2023 as an example, provided as an embodiment of the present invention; Figure 11 A schematic diagram illustrating the distribution of fault slip, slip standard deviation, and slip variation coefficient, using the Mw7.8 earthquake in a certain area in 2023 as an example, provided for embodiments of the present invention; Figure 12 A schematic diagram illustrating the distribution of fault slip, slip standard deviation, and slip variation coefficient, using the Mw7.7 earthquake in a certain area in 2023 as an example, provided for embodiments of the present invention; Figure 13 This is a schematic diagram of the structure of the source rupture process inversion system based on multi-source observation data provided in this embodiment of the invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0022] Example 1, Figure 1 A flowchart of a method for inverting the source rupture process based on multi-source observation data is shown, as follows: Figure 1 As shown, a method for inverting the source rupture process based on multi-source observation data includes the following steps: S100: Acquire seismic and geodetic observation data from multiple sources; Specifically, the earthquake observation data includes: strong motion three-component waveform records obtained from strong-motion stations near the epicenter; vertical P-wave records and tangential SH-wave records obtained from far-field (30° to 90°) seismograph stations near the epicenter. Geodetic observation data includes: surface deformation field data obtained from interferometric synthetic aperture radar (InSAR) covering the earthquake rupture zone; and surface displacement data obtained from Global Navigation Satellite System (GNSS) stations near the epicenter.

[0023] S200: Construct an observation value vector based on the earthquake observation data and the geodetic observation data; Specifically, observation datasets from different sources have different observation values. For example, strong ground motion record data is a time series of the three components of ground motion velocity, far-field seismic record data is a time series of ground motion displacement, InSAR deformation field data is the satellite line-of-sight component of coseismic surface displacement, and GNSS displacement data is the three components of static surface displacement.

[0024] Different datasets are processed into vector form separately and then combined to form an observation vector. The specific expression is as follows: , in, Indicates the first A vector of observed values ​​from an observation dataset. Indicates the first The first observation dataset There are [number] observations, and the total number of observation datasets is [number]. The number of observations in each dataset is , Observation dataset Together they form the observation vector.

[0025] S300: Establish a fault geometry model, adaptively discretize the fault plane, and obtain multiple rectangular sub-fault dislocation elements; The seismogenic fault has a relatively complex geometric shape, and its strike and dip angle may change in space. First, the position of the top edge of the fault is specified as the fault trace, and then the fault dip angle is specified sequentially along the top trace. If the strike or dip angle of the fault changes in space, a fault geometric model with a curved surface shape can be established. like Figure 2 As shown, S300 specifically includes: S310: Use A continuous distribution of planar Cartesian coordinates along the fault strike. To describe the location of the top edge of the fault, i.e., the fault trace, where East-west coordinates, Use north-south coordinates and specify the fault dip angle for each trace description point. ; The cumulative length of each descriptive point in the fault trace Obtained through the following formula: , S320: The fault trace is discretized using the principle of linearization, and the coordinates and dip parameters of the control points of the fault trace are calculated by piecewise linear interpolation; The specific process is as follows: S321: Determine the number of rectangular sub-faults along the strike and dip of the fault, for a length of... (Step length along the direction), width is Sub-faults (along the dip step). Number of sub-faults in the strike direction. The total length of the fault trace The length of the sub-fault is Calculation, through Ensure that at least two sub-faults are defined to maintain geometric rationality. , Number of dip direction sub-faults Based on the fault dip width with sub-fault width Similarly, the minimum constraint value is 2. , A slight correction was made to the spatial distribution step size of the sub-faults: , .

[0026] S322: Discretization of fault traces is achieved through linear interpolation. First, a dimension of... The zero matrix stores the three-dimensional coordinates of the control points of each sub-fault's spatial location. ,in East-west coordinates, North-south coordinates For depth coordinates, construct a structure with length... The zero vector stores the tilt angle parameters that vary along the azimuth. The coordinates of the top control point of the shallowest sub-fault are used as the fault trace control points, where the first and last points are directly taken from the coordinates of the fault trace description points. and specify the depth of the fault crest as , , , , , , , , , For intermediate discrete points in the fault trace control points First, fix its depth. ,in, For a fixed depth value, traverse each control point segment, and then determine the cumulative distance along the direction of the current discrete point. Whether it falls into the first The cumulative length interval of each description point If the conditions are met, then calculate the interpolation coefficients. : , based on For the first and the Linear interpolation is performed on the planar coordinates and inclination angles of each description point to obtain the planar coordinates and inclination angles of the current trace control point. , , , After completing the interpolation calculation for the current control point, terminate the traversal of the control point segments and continue processing the next fault trace control point.

[0027] S330: Determine the coordinates of the remaining sub-fault control points; first, calculate the center coordinates of the fault trace control points as the central benchmark for subsequent average strike fitting. Then, perform an arithmetic mean of the plane coordinates of the fault trace control points to obtain the center coordinates. , , .

[0028] The control points of the fault trace are centered, and the coordinates of the control points relative to the center are calculated. The sum of squares and the sum of cross products: , , , in , They are respectively , Dispersion of direction for Cumulative covariance in direction.

[0029] like ,Right now Large directional dispersion, regression slope average strike of fault Obtained from the following formula: , in , It is the arctangent function in the four quadrants; like ,Right now Large directional dispersion, regression slope average strike of fault Represented as: , in .

[0030] average trend Corrected to :like ,but ;like Then take the modulus .

[0031] For each sub-fault control point below the fault trace, first traverse along the strike ( ), and then traverse along the tendency ( Calculate the three-dimensional coordinates of each control point. Based on the dip fault step size. With tilt angle Calculate along the direction of the first One, along the tendency The horizontal projection distance of each control point relative to the upper edge of the fault With vertical depth increment : , , Combined with the average strike of the fault The three-dimensional coordinates of the control points are obtained by combining the azimuth correction of the horizontal projection distance with the vertical depth increment. , , .

[0032] S340: Calculate the geometric parameters of each sub-fault based on the three-dimensional coordinates of the sub-fault control points, and traverse the sub-fault mesh along the strike and dip directions of the faults. , Each sub-fault is calculated, including its center coordinates, strike, normal vector, area, dip angle, and dimensional parameters.

[0033] First, accumulate the number of sub-faults ( And by taking the arithmetic mean of the coordinates of the four control points of the sub-fault, the eastward direction of the sub-fault center was obtained. ), North ( ), vertical ( )coordinate, , , .

[0034] sub-fault strike Take the average of the azimuth angles of the two sidelines along its direction, and correct it to... interval, , like Then let .

[0035] The normal vector of the sub-fault is obtained by taking the cross product of the two diagonal vectors of the sub-fault, and then the area of ​​the sub-fault is calculated. The two diagonal vectors of the sub-fault are respectively... , ,in, , , , , , , The cross product of the diagonal vectors is the normal vector of the sub-fault. , , , , The area of ​​the sub-fault is half the magnitude of the normal vector. .

[0036] The dip angle of a sub-fault is the angle between the normal vector and the vertical direction. , The sub-fault length is corrected to the average length of the two edge segments along the strike direction. , The width of the sub-fault is corrected to the ratio of its area to its length. .

[0037] S400: Based on the seismic observation vector and the fault geometric model, calculate the Green's function and establish the joint inversion equation; Specifically, S400 is as follows: S410: Establish linear equations for inverting the source rupture process based on the source representation theorem; In the focal representation theorem, each point on the fault plane is displaced at a spatial point. The resulting displacement is: , in For being in station displacement Directional components, This represents a fault dislocation discontinuity surface. Let be the elastic modulus of the medium at the fracture surface. for The normal cosine of the fault plane direction. for Constantly acting on the fault plane place Unit concentrated force in the direction The location generates the Green's function in Displacement generated at time Directional component.

[0038] After discretization in time and space, the motion time history (displacement, velocity, or acceleration) of a point in space can be expressed as the weighted superposition of the rupture effects of a series of discrete sub-faults on the fault plane within multiple time windows. The source representation theorem then changes to the following form: , in The number of time windows. The number of sliding components is used; generally, two mutually orthogonal sliding components are used to represent the sliding vector. For time window functions, , This is the distance between the sub-fault and the rupture initiation point. The interval of the time window. Indicates the first The upper edge of the fault The first sliding direction The coefficients of each time window are the parameters of the source model to be solved in the inversion process. Time window function It can be taken as a smooth slope function, representing the slip time function or seismic moment time function, with a rise time of... , .

[0039] Based on this, the linear equations required for the inversion are constructed, and the coefficient matrix is ​​composed of the convolution of the theoretical Green's function and the time window function from the rupture of each sub-fault point source to the response of the observation station. The time window functions of each sub-fault have coefficients that need to be solved to form a vector. Earthquake observations form a vector Taking the use of observational data as a record of ground motion velocity as an example, the inversion equation is as follows: , in , , , , The number of earthquake observations. .

[0040] S420: Calculate the Green's function based on the three-dimensional coordinates of each sub-fault in the fault geometry model and the geographic coordinates of seismic observations; The Green's function is a function used to solve non-homogeneous differential equations with initial or boundary conditions. In the source representation theorem, it represents the response at the observation point produced by a unit concentrated force applied by a point source. After discretizing the fault geometry model, where each sub-fault satisfies the point source assumption, the Green's function at a specified observation location can be calculated based on its spatial coordinates.

[0041] Different methods were used to calculate the Green's function for observational data with varying epicentral distances. For strong ground motion records and GNSS displacement data near the epicenter, as well as InSAR deformation field data covering the earthquake rupture zone, the FK software was used to calculate the Green's function. FK software is a high-precision surface response calculation method suitable for layered half-space media. By performing global integration in the frequency-wavenumber domain, it rigorously solves for the propagation, reflection, and conversion effects of various body waves and surface waves. It can output a complete surface response, including near-field strong ground motion, static displacement field, and long-term deformation field, in a single operation. It achieves a natural transition and unified expression from high-frequency ground motion to quasi-static surface deformation without distinguishing between dynamic and static fields, making it particularly suitable for simulating near-field strong ground motion records, GNSS static displacement, and InSAR surface deformation observations.

[0042] For seismic records at far-field epicentral distances (30° to 90°), the Green's function was calculated using Multitel3 software. Multitel3 is a dedicated Green's function calculation method for teleseismic body wave data. Based on layered medium wave theory and ray propagation path construction, it can automatically distinguish between direct P-waves / SH waves and core-mantle boundary reflection phases under far-field conditions, and independently calculate the travel time, amplitude, and geometric diffusion attenuation of different phases. It can output the Green's function for specific phases without calculating the full waveform, effectively suppressing phase aliasing and high-frequency noise, and significantly improving the stability and resolution of teleseismic waveform simulation.

[0043] In calculating the Green's function for different types of observation data, a subsurface velocity structure model matching the epicentral distance needs to be adopted. For strong ground motion records, GNSS static displacement, and InSAR surface deformation observations from the near-field of the epicenter to the regional area, the subsurface velocity structure model is constructed with reference to seismic tomography results or regional geological data of the earthquake area, and typically includes stratum thickness (km) and P-wave velocity (…). (km / s), S-wave velocity ( (km / s) and formation density ( g / cm 3 A velocity model for layered media is presented. For velocity models containing only partial parameters, the missing parameters are supplemented using empirical conversion formulas to obtain a complete velocity model. The conversion relationships between P-wave velocity, S-wave velocity, and formation density are as follows: , , , .

[0044] For far-field seismic records at epicentral distance, a layered medium velocity model is constructed using the CRUST1.0 model. The CRUST1.0 model sets nodes at 1° intervals and covers the entire globe, providing complete parameters for formation thickness, P-wave velocity, S-wave velocity, and formation density, which can meet the needs of widely distributed far-field stations. In actual calculations, model parameters at several surrounding nodes are typically interpolated based on the latitude and longitude of the station location to establish a layered medium velocity model at the corresponding location.

[0045] S430: Construct joint inversion equations for multi-source data using inversion equations for various types of observation data; For any single type of observation data, a linear inversion equation can be established based on the source representation theorem. In the joint inversion, the joint inversion equation is constructed by combining the linear inversion equations corresponding to each type of observation data.

[0046] The various types of observation data can be further categorized into dynamic data (e.g., strong ground motion waveform records, seismic waveform records, GNSS high-frequency displacement records) and static data (e.g., InSAR deformation field data, GNSS static displacement records) based on whether the observed values ​​are dynamic sequences that change over time. For dynamic data, multiple time windows are set as sampling of the source time function in the time domain to invert the source rupture process; while static data is typically used to invert fault slip distribution and constrain the zero-frequency characteristics of the source. The joint inversion equation obtained by linearly superimposing the inversion equations corresponding to dynamic and static data is as follows: , in, This represents the Green's function coefficient matrix and adjustment coefficients corresponding to different types of observation datasets. Simultaneously assigned the Green's function matrix and observation data vector , The inversion parameters to be solved to describe the source rupture process are indicated by subscripts. and These represent dynamic and static data, respectively. Expanding the above formula further yields... , in, and These represent the Green's function matrix elements corresponding to dynamic and static data, respectively. Since the Green's function for dynamic data is itself a dynamic time series, to simultaneously utilize both dynamic and static data in joint inversion, the Green's function matrix corresponding to the static data needs to be extended in the time domain to meet the requirement of setting multiple time windows for dynamic data. .

[0047] Adjustment coefficient Includes normalization coefficients and data weight coefficients Two parts of contribution, .

[0048] Because the dimensions and orders of magnitude of observations differ significantly among different types of data, normalization is used to map all data to a uniform scale, eliminating these differences and avoiding the neglect of the constraints imposed by smaller observations on the inversion results. This ensures that data from different dimensions have equal weighting for the dataset. Normalization coefficient Use the inverse of the L2 norm of the corresponding dataset or time series. , To maximize the resolving power of different data types while suppressing the influence of errors and noise in observational data to enhance the stability of the inversion calculation, different datasets are assigned corresponding data weighting coefficients. After normalization, the appropriate range of the weight coefficients for each observation dataset is between 0 and 1. By using grid search, the final value of the weight coefficients is determined when the fitting residuals of each dataset tend to have no significant change.

[0049] S500: Solve the joint inversion equations to obtain the inverted source rupture process model.

[0050] Specifically, S500 is as follows: S510: An improved spatial smoothing constraint is applied to the inversion of the source rupture process to make the slip or seismic moment release of adjacent sub-faults equal. At the same time, an improved temporal smoothing constraint is introduced to make the slip or seismic moment release of the same sub-fault equal in adjacent time windows. Since the source rupture originates in the fault region that first meets the critical conditions, and the surrounding medium gradually reaches its strength limit and ruptures through stress transmission, this indicates a strong correlation between the rupture behaviors of adjacent media in the source region. Therefore, spatial smoothing constraints need to be applied to the source rupture process inversion to ensure that the slip or moment release of adjacent sub-faults is approximately equal. Simultaneously, rupture propagation exhibits a delay effect dominated by rupture velocity, and its temporal evolution shows continuous and gradual characteristics. Therefore, temporal smoothing constraints also need to be introduced to ensure that the slip or moment release of the same sub-fault is approximately equal in adjacent time windows. Smoothing Constraints By suppressing abrupt changes in the inversion solution in the spatial and temporal domains, the stability of the inversion problem is enhanced, and this approach conforms to the physical laws of source rupture. Traditional Laplace smoothing constraints are purely mathematical, equal-weighted difference forms that only enforce the continuity of adjacent spatiotemporal parameter values, failing to fully incorporate the physical laws of rupture propagation and stress transfer. The constraint weights are disconnected from the physical correlation of the rupture process. Furthermore, traditional Laplace-style spatial and temporal constraints are two independent weighting systems, requiring separate hyperparameter settings to adjust constraint strength. Parameter selection is highly subjective and computationally expensive, and the uniformity of the physical scale of spatiotemporal constraints cannot be guaranteed, easily leading to imbalances between spatiotemporal constraint strengths. Therefore, the proposed spatial smoothing constraint... With time smoothing constraints An improvement is made to the difference form of the Laplace second derivative, which differs from a simple spatiotemporal unification based solely on mathematical form. Instead, it transforms into a spatiotemporal scale adaptive constraint driven by a fractured physical process, with only a single smooth constraint hyperparameter set globally. To achieve synchronous control of spatial and temporal smooth constraint strength. , The spatial smoothing constraint is typically applied to one sub-fault (model parameters are...). )and Adjacent sub-faults (model parameters are) ), specifically in the form of: , in, , , In the formula As a physical reference parameter, it represents the maximum moment when the first time window of each sub-fault is triggered, directly corresponding to the overall rupture spatial scale of the earthquake source. It is the total time boundary of the entire rupture propagation and stress transfer, providing a unified physical reference for spatial constraints. and These are the rupture trigger times of the constrained sub-fault and its adjacent sub-fault, respectively. The physical correlation parameter represents the rupture triggering time difference between the adjacent sub-fault and the current sub-fault, directly quantifying the degree of rupture propagation correlation between two adjacent sub-faults during the rupture process; the proportional term... Based on the overall rupture scale of the earthquake source, the degree of rupture correlation between adjacent sub-faults is normalized and quantified. The larger the ratio, the higher the degree of rupture correlation between adjacent sub-faults, and stronger smoothing constraints need to be applied to conform to the continuity of rupture propagation. The smaller the ratio, the lower the degree of rupture correlation between adjacent sub-faults, and the constraint intensity needs to be appropriately reduced to preserve the local characteristics of the rupture process.

[0051] The time smoothing constraint is typically applied within one time window (model parameters are...). ) and 2 intervals Adjacent time windows (model parameters are) and ), specifically in the form of: , in, , .

[0052] In the formula It still represents the maximum moment triggered by the first time window of each sub-fault, and uses the same physical benchmark as the spatial constraints. The interval between adjacent time windows directly quantifies the degree of temporal correlation in the rupture process within those windows. Because the weights of the spatial and temporal smoothing constraints employ completely homogeneous physical benchmarks and consistent quantization logic, scale unification of the spatial and temporal smoothing constraints is achieved at a physical level, rather than being a purely mathematical concatenation. Therefore, only a single hyperparameter is required. This enables synchronous global control of the spatial and temporal constraint strength, overcoming the limitation of traditional methods where spatial and temporal constraints need to be adjusted separately. Furthermore, the physical meaning of the selected parameters is clear, avoiding the problem of constraint strength imbalance caused by human subjectivity.

[0053] S520: In the inversion of the source rupture process, the reverse slip of the fault is not considered. The slip of each sub-fault is set to vary freely in half-space in the main slip direction. By solving the joint inversion equation with non-negative inequality constraints, the spatiotemporal distribution of fault slip characterizing the source rupture process is obtained. According to the elastic rebound theory, fault slippage occurs when accumulated tectonic stress exceeds the strength of the medium. When slippage occurs between the hanging wall and footwall of a fault, strong friction is generated, causing the fault slippage to resemble damped vibration, gradually decaying until it eventually heals. Its direction of motion is essentially consistent with the direction of regional tectonic stress, making it difficult for slippage to occur completely opposite to the stress direction. Therefore, reverse slippage of the fault is not considered in the inversion of the focal rupture process. The slip angle specified by the focal mechanism solution can be used as the main slip direction of the fault. Generally, non-negative constraints are used to constrain the slip direction of sub-faults within ±45° of this main slip direction, forcing the entire fault to have only a single type of fault movement mechanism. This is difficult to adapt to the complex situations commonly found in terrestrial earthquakes, such as "main rupture accompanied by other types of local slip components," "multi-segmental change mechanisms of cascade earthquakes," and "mixed strike-slip and dip-slip rupture," resulting in a systematic deviation between the inversion results and the actual rupture process. Therefore, under the premise of strictly prohibiting reverse slip and ensuring inversion stability, taking the main slip direction driven by tectonic stress as the benchmark, and by widening the range of variation of the sub-fault slip direction to within ±90° of the main slip direction, allowing the sub-fault slip to freely vary within half-space of the main slip direction, it can naturally accommodate at least two types of fault motion mechanisms. This makes the inversion model no longer limited by the preset source mechanism, effectively restoring the differences in source mechanism caused by spatial variations in geometry and stress of the real seismogenic fault, and significantly improving the physical rationality and authenticity of the inversion model. Its specific form is as follows: , in, For the pre-defined main sliding direction, For the first fault plane The slip vector of a sub-fault. Written in matrix form as: ,in For diagonal matrices: .

[0054] At this point, for the form of The joint inversion equations, under the application of smoothing constraints Regularization (regularization parameter is) ) and inequality constraints In this case, the inversion problem is: .

[0055] Due to the Hessian matrix of this inversion problem It is a symmetric positive definite matrix, and It is a standard linear inequality constraint and naturally belongs to a convex set. Therefore, this inversion problem belongs to a convex quadratic programming problem, and there exists an optimal solution under the condition of satisfying the (Karush-Kuhn-Tucker, KKT) condition.

[0056] First, regarding inequality constraints Introducing Lagrange multiplier vectors Construct the Lagrange function: .

[0057] The optimal solution must satisfy the following KKT conditions: , , .

[0058] The optimal solution can be obtained from the above equation. ,in , which is the least squares benchmark solution to the inversion problem without inequality constraints. Further, the coefficients... Absorption into Lagrange multipliers To obtain the optimal solution The form is: .

[0059] Then, slack variables are introduced. Transform inequality constraints into equality constraints: , .

[0060] The optimal solution Substituting into the above equation, we get: , .

[0061] At this point, the inversion problem with inequality constraints is transformed into a standard linear complementarity problem (LCP), and the Lemke algorithm is used to solve for the vectors. Then the optimal solution can be obtained. Ultimately, the spatiotemporal distribution of fault slip, characterizing the earthquake source rupture process, was obtained.

[0062] Following the above inversion solution process, the fitting results of various types of observation datasets in the joint inversion, using the Mw7.8 and Mw7.7 double earthquakes in a certain area in 2023 as an example, are as follows: Figure 3 - Figure 10 As shown, these correspond to the InSAR data, GNSS horizontal direction data, strong ground motion waveform data, and far-field seismometry P-wave data used for the inversion, respectively. Figures 9-10In the waveform recording, COLA, KDAK, TIXI, etc., on the left are the station names. It can be seen that the observation and fitted data of each type achieve good consistency, meaning that the inversion model accurately fits the observation data, indicating that the inversion method proposed in this invention can effectively achieve joint inversion of multi-source data.

[0063] S530: The Jackknife test is used to quantitatively obtain the uncertainty of fault slip distribution, reflecting the impact of observational data on the inversion results of the source rupture process; If the slip of the sub-fault has a small dispersion in all Jackknife inversion results, it indicates that the parameter is weakly affected by specific observation data and has low uncertainty; a large dispersion indicates that the parameter is sensitive to specific observation data and has high uncertainty.

[0064] First, the results of the source rupture process inversion based on the complete observation dataset are used as the benchmark solution. ,in Number of sub-faults Indicates the first The slip of each sub-fault. The parameter settings in the fixed inversion, such as sub-fault division, smoothing constraints, and time window function form, are kept consistent in all subsequent Jackknife inversions to ensure comparability of results.

[0065] Set the number of repeated inversions in the Jackknife test to [number]. (Generally 100 to 200 times), the proportion of observation data removed in each inversion is... (Generally 10% to 20%), thus constructing The Jackknife sample dataset was used. Considering the spatial and path correlations in InSAR deformation, GNSS displacement, and seismic waveform data, sampling on a per-observation basis would disrupt the inherent correlations between data points, leading to systematic biases in uncertainty assessments. Therefore, highly correlated spatially or path-dependent observations were divided into data blocks, and sampling and removal were performed on a block-by-block basis, preserving the original correlation structure of the data. For InSAR deformation and GNSS displacement data, observations from stations with adjacent pixels or small spatial distances showed strong correlations. For example, in InSAR data, blocks with a spatial size not exceeding 10 km × 10 km and containing no more than 2% of the total pixels were grouped into one data block; in GNSS data, blocks with station distances within 5 to 10 km and not exceeding 10% of the total number of stations were grouped into one data block. For strong-motion and seismic records, stations with similar azimuth angles and epicentral distances exhibit a strong correlation between waveform path effects and observation errors. Data can be divided into blocks according to specific azimuth and epicentral distance intervals. For example, strong-motion data can be divided into two-dimensional blocks at 30° azimuth intervals and three epicentral distance ranges (0~50 km / 50~100 km / >100 km); seismic data can be divided into two-dimensional blocks at 45° azimuth intervals and two epicentral distance ranges (30°~60° / 60°~90°). After data block division, all observation data sampling uses these blocks as the basic sampling unit. Considering that Green's function calculation is a time-consuming step in the inversion process, a complete Green's function matrix containing all observation points corresponding to all sub-faults needs to be calculated before performing the Jackknife test. This way, during each Jackknife sampling inversion, only the corresponding sub-matrices need to be extracted from the complete Green's function matrix based on the retained observation data, eliminating the need for repeated calculations of the Green's function and significantly improving computational efficiency.

[0066] After that, Each of the Jackknife sample datasets was subjected to source rupture process inversion with the same parameter settings as the original inversion, resulting in... Jackknife inversion solution: .based on For a set of Jackknife solutions, calculate the statistical discrete index of the slip volume of each sub-fault. Commonly used statistics characterizing its uncertainty include: standard deviation. 95% confidence interval, coefficient of variation.

[0067] Standard deviation It is the most intuitive indicator of absolute uncertainty, with dimensions consistent with slip volume, making it easy to directly determine the fluctuation range of sub-fault slip volume. The larger the value, the higher the uncertainty. Its expression is: , in , for the first Jackknife mean of individual fault slip.

[0068] The confidence interval reflects the boundary characteristics of the uncertainty in the sub-fault slip. Based on the assumption that the observed data satisfy a normal distribution, the 95% confidence interval... Quantile calculation based on Jackknife solution .

[0069] coefficient of variation It reflects the proportion of the uncertainty of sub-fault slip relative to the slip amplitude, that is, the fluctuation amplitude corresponding to a unit slip, and has comparability across regions and amplitudes. A low value indicates that the inversion results of the slip in this region are highly reliable and less affected by observation data. Its expression is: .

[0070] The fault slip distribution benchmark solution, slip standard deviation distribution, and slip variation coefficient distribution obtained using the Mw7.8 and Mw7.7 double earthquakes in a certain area in 2023 as an example are as follows: Figure 11 and Figure 12 As shown, the maximum slip in the Mw7.8 earthquake was 6.6 m, the maximum slip standard deviation was 1.0 m, and the maximum slip coefficient of variation was 1.0; in the Mw7.7 earthquake, the maximum slip was 7.6 m, the maximum slip standard deviation was 1.4 m, and the maximum slip coefficient of variation was 1.8. It can be seen that although the maximum slip standard deviation in both earthquakes occurred in areas of significant fault slip, it was much smaller than the baseline slip at the corresponding locations, and the coefficient of variation in areas of significant slip was generally within 0.2; while fault areas with larger coefficients of variation did not show significant slip, meaning they did not have a practical impact on the stability of the inversion results. The inversion method proposed in this invention can intuitively assess the uncertainty of the obtained source parameters, and joint inversion can significantly enhance the accuracy of the inversion results.

[0071] The inversion method proposed in this embodiment can reconstruct the source rupture process by combining multi-source observation data from seismic and geodetic surveys. By combining the resolution capabilities of different observation data, it can obtain source parameters more accurately. It can establish a fault geometric model with continuously varying strike and dip angles, and on this basis, invert the source rupture process with complex propagation modes, providing a new and effective means for studying the source characteristics of complex seismogenic faults.

[0072] Example 2: This invention also provides a source rupture process inversion system based on multi-source observation data. The system employs a source rupture process inversion method based on multi-source observation data from Example 1, such as... Figure 13As shown, the system includes: The observation data acquisition module 600 is used to acquire seismic observation data and geodetic observation data from various sources. The multi-source data processing module 610 is used to process the seismic observation data and geodetic observation data to construct a seismic observation value vector. The fault geometry generation module 620 is used to establish the fault geometry model and perform adaptive discretization sampling to obtain multiple rectangular sub-fault dislocation elements. The inversion equation construction module 630 is used to calculate the Green's function and establish a joint inversion equation based on the earthquake observation vector and fault geometry model. The joint inversion solution module 640 is used to solve the joint inversion equations to obtain the inverted source rupture process model.

[0073] Example 3: The present invention also provides an electronic device, including one or more processors and a memory.

[0074] A processor can be a central processing unit (CPU) or other form of processing unit with data processing and / or instruction execution capabilities, and can control other components in an electronic device to perform desired functions.

[0075] The memory may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. The volatile memory may include, for example, random access memory (RAM) and / or cache memory. The non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc. One or more computer program instructions may be stored on the computer-readable storage medium, and a processor may execute the program instructions to implement the source rupture process inversion method based on multi-source observation data described in any embodiment of this application, and / or other desired functions. Various contents such as initial extrinsic parameters and thresholds may also be stored in the computer-readable storage medium.

[0076] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A method for inverting the source rupture process based on multi-source observation data, characterized in that, The method includes the following steps: S100: Acquire seismic and geodetic observation data from multiple sources; S200: Construct an observation value vector based on the earthquake observation data and the geodetic observation data; S300: Establish a fault geometry model, adaptively discretize the fault plane, and obtain multiple rectangular sub-fault dislocation elements; S400: Based on the observed value vector and the fault geometric model, calculate the Green's function and establish the joint inversion equation; S500: Solve the joint inversion equations to obtain the inverted source rupture process model; S500 specifically refers to: S510: Apply spatial smoothing constraints to the inversion of the source rupture process to make the slip or seismic moment release of adjacent sub-faults equal. At the same time, introduce temporal smoothing constraints to make the slip or seismic moment release of the same sub-fault equal in adjacent time windows. S520: In the inversion of the source rupture process, the reverse slip of the fault is not considered. The slip of each sub-fault is set to vary freely in half-space in the main slip direction. By solving the joint inversion equation with non-negative inequality constraints, the spatiotemporal distribution of fault slip characterizing the source rupture process is obtained. S530: The Jackknife test is used to quantitatively obtain the uncertainty of fault slip distribution, reflecting the impact of observational data on the inversion results of the source rupture process.

2. The method for inverting the source rupture process based on multi-source observation data according to claim 1, characterized in that, In S100, the earthquake observation data includes: strong motion three-component waveform recording data obtained from strong motion meter stations near the epicenter, vertical component P-wave recording data and tangential component SH-wave recording data obtained from far-field seismograph stations near the epicenter; the geodetic observation data includes: surface deformation field data obtained from synthetic aperture radar satellite interferometry covering the earthquake rupture range, and surface displacement data obtained from global navigation satellite system stations near the epicenter.

3. The method for inverting the source rupture process based on multi-source observation data according to claim 1, characterized in that, In step S200, different datasets are processed into vector form and then combined to form an observation vector, as specifically expressed as: ， in, Indicates the first A vector of observed values ​​from an observation dataset. Indicates the first The first observation dataset There are [number] observations, and the total number of observation datasets is [number]. The number of observations in each dataset is , Observation dataset Together they form the observation vector.

4. The method for inverting the source rupture process based on multi-source observation data according to claim 1, characterized in that, Specifically, S300 is: S310: Use planar Cartesian coordinates that are continuously distributed along the fault strike to describe the location of the fault crest edge, i.e., the fault trace, and specify the fault dip angle for each trace point; S320: The fault trace is discretized using the principle of linearization, and the coordinates and dip parameters of the control points of the fault trace are calculated by piecewise linear interpolation; S330: To determine the coordinates of the remaining sub-fault control points, the center coordinates of the fault trace control points must first be calculated. The plane coordinates of the fault trace control points are then arithmetically averaged and the average strike is calculated. Finally, the three-dimensional coordinates of each sub-fault control point are calculated by traversing the fault trace control points according to the average strike and the specified dip angle. S340: Traverse the sub-fault mesh along the strike and dip directions of the fault, and calculate the coordinates of the calculation center, strike, normal vector, area, dip angle, and size parameters of each sub-fault.

5. The method for inverting the source rupture process based on multi-source observation data according to claim 4, characterized in that, Specifically, S320 is: S321: Determine the number of rectangular sub-faults along the strike and dip of the fault; S322: First, a zero matrix is ​​constructed to store the three-dimensional coordinates of the control points of each sub-fault's spatial location, and a zero vector is constructed to store the dip angle parameters that vary along the strike. The coordinates of the top control point of the shallowest sub-fault are used as the fault trace control points. For the intermediate discrete points among the fault trace control points, their depths are first fixed. ,in, For a fixed depth value, traverse each control point segment, and then determine whether the cumulative distance along the direction of the current discrete point falls within the range of the first segment. If the cumulative length interval of the description points satisfies the condition, then the interpolation coefficients are calculated. ;based on For the and the Linear interpolation is performed on the plane coordinates and dip angle of each description point to obtain the plane coordinates and dip angle of the current fault trace control point. After completing the interpolation calculation of the current control point, the traversal of the control point segment is terminated, and the next fault trace control point is processed.

6. The method for inverting the source rupture process based on multi-source observation data according to claim 5, characterized in that, In step S321, the number of rectangular sub-faults along the fault strike and dip is determined, for a length of... Width is Sub-faults, strike direction, number of sub-faults The total length of the fault trace The length of the sub-fault is Calculation, through Ensure that at least two sub-faults are defined to maintain geometric rationality. ， Number of dip direction sub-faults Based on the fault dip width with sub-fault width Similarly, the minimum constraint value is 2. ， A slight correction was made to the spatial distribution step size of the sub-faults: ， 。 7. The method for inverting the source rupture process based on multi-source observation data according to claim 5, characterized in that, In step S322, the discretization of the fault traces is achieved through linear interpolation. First, a dimension of [missing information] is constructed. The zero matrix stores the three-dimensional coordinates of the control points of each sub-fault's spatial location. ,in East-west coordinates, North-south coordinates For depth coordinates, construct a structure with length... The zero vector stores the tilt angle parameters that vary along the azimuth. The coordinates of the top control point of the shallowest sub-fault are used as the control points of the fault trace, where the first and last points are directly taken from the coordinates of the fault trace description points. and specify the depth of the fault crest as , ， ， ， ， ， ， ， ； For intermediate discrete points in the fault trace control points First, fix its depth. Iterate through each control point segment and then determine the cumulative distance along the direction of the current discrete point. Whether it falls into the first The cumulative length interval of each description point If the conditions are met, then calculate the interpolation coefficients. : ， based on For the and the Linear interpolation is performed on the planar coordinates and inclination angles of each description point to obtain the planar coordinates and inclination angles of the current trace control point. ， ， ， After completing the interpolation calculation for the current control point, terminate the traversal of the control point segments and continue processing the next fault trace control point.

8. The method for inverting the source rupture process based on multi-source observation data according to claim 1, characterized in that, Specifically, S400 is: S410: Establish linear equations for inverting the source rupture process based on the source representation theorem; S420: Calculate the Green's function based on the three-dimensional coordinates of each sub-fault in the fault geometry model and the geographic coordinates of seismic observations; S430: Construct joint inversion equations for multi-source data using inversion equations for various types of observation data.

9. A source rupture process inversion system based on multi-source observation data, characterized in that, The system employs the source rupture process inversion method based on multi-source observation data as described in any one of claims 1-8, and the system comprises: The observation data acquisition module is used to acquire seismic observation data and geodetic observation data from various sources. The multi-source data processing module is used to process the earthquake observation data and geodetic observation data to construct an earthquake observation value vector; The fault geometry generation module is used to establish the fault geometry model and perform adaptive discretization sampling to obtain multiple rectangular sub-fault dislocation elements. The inversion equation construction module is used to calculate the Green's function and establish joint inversion equations based on the seismic observation vector and fault geometry model. The joint inversion solution module is used to solve the joint inversion equations to obtain the inverted source rupture process model.

Images