Method for high resolution p and s wave velocity ratio time series inversion based on seismic travel times

Abstract

The application relates to a method for high-resolution P and S wave velocity ratio time series inversion based on seismic wave arrival time, which comprises the following steps: S1, obtaining P wave and S wave arrival time difference data of a seismic pair; S2, performing linear relationship fitting on the P wave and S wave arrival time difference data to extract average wave velocity ratio of the seismic pair; S3, constructing an initial wave velocity ratio time series according to a seismic time of a seismic sequence; S4, constructing a time series inversion matrix and performing inversion solving on the time series inversion matrix to obtain a wave velocity ratio time series; S5, judging whether the wave velocity ratio time series is stable, if the wave velocity ratio time series is stable, then returning to S3 with time series unknown number plus one, and if the wave velocity ratio time series is not stable, then outputting the inverted wave velocity ratio time series. The application can more accurately reveal the dynamic change rule of the wave velocity ratio of a natural seismic zone or an artificial blasting area.

Description

Technical Field

[0001] This invention relates to the field of underground structure inversion technology, and in particular to a method for inverting high-resolution P- and S-wave velocity ratio time series based on seismic wave arrival time. Background Technology

[0002] In geophysical exploration and seismological research, inverting the velocity structure of subsurface media is crucial for understanding crustal deformation, fault activity, and earthquake gestation mechanisms. The velocity ratio (Vp / Vs; hereinafter referred to as wave velocity ratio) between P-waves (P-waves) and S-waves (S-waves) is particularly sensitive to processes such as subsurface fluid transport and pore pressure changes. Currently, commonly used Vp / Vs wave velocity ratio tomography techniques have achieved significant results in static spatial structure inversion, but still have considerable limitations in temporal resolution and dynamic monitoring. Existing dynamic wave velocity ratio monitoring techniques mainly include the Hilda method and the double-difference wave velocity ratio method. The Hilda method utilizes the arrival time t of P-waves from multiple seismic events recorded by stations. P and the time difference between P-wave and S-wave (t) S -t P Add one to the fitted slope to obtain the wave speed ratio, i.e., Vp / Vs = (t S -t P ) / t P +1. This method primarily reflects the average wave velocity ratio of the subsurface medium between the earthquake and the station, but it cannot reflect the wave velocity ratio of a specific local area underground, thus its spatial resolution is insufficient. The double-difference wave velocity ratio method, on the other hand, mainly utilizes the P-wave arrival time difference Δt recorded by multiple stations for the same pair of earthquakes. P and transverse wave to time difference Δt S The wave speed ratio, Vp / Vs = Δt, can be obtained from the slope of the fitted straight line. S / Δt P This method primarily reflects the wave velocity ratio in local areas between earthquake pairs, thus improving spatial resolution. However, because the occurrence times of earthquake pairs may be far apart, and similar earthquake events cannot be received simultaneously by many stations, the temporal resolution and fitting error of this method are significantly affected by the original data. In summary, traditional methods for obtaining underground structures based on wave velocity ratios suffer from drawbacks such as low temporal resolution, reliance on earthquake location accuracy, and large errors in wave velocity ratio structure calculations.

[0003] The prior art discloses a method for real-time acquisition of microseismic wave velocity in hard rock tunnels constructed by drill-and-blast method. This method is mainly based on the Heda method and is used for real-time acquisition of wave velocity. However, it relies on the accuracy of seismic positioning and has no ability to distinguish the changes in wave velocity over time.

[0004] Therefore, there is an urgent need to develop a new wave velocity ratio calculation method based on inversion technology, which can maximize the use of P-wave and S-wave arrival time data of earthquakes to achieve high-resolution acquisition of crustal wave velocity ratio time series, thereby providing a more reliable technical means for dynamic monitoring of underground media. Summary of the Invention

[0005] The purpose of this invention is to provide a method for retrieving high-resolution P- and S-wave velocity ratio time series based on seismic wave arrival time inversion, focusing on refined modeling and inversion in the time dimension, thereby revealing the dynamic variation law of the wave velocity ratio more accurately.

[0006] To achieve the above objectives, the present invention provides the following solution:

[0007] Methods based on seismic wave arrival time inversion to retrieve high-resolution P- and S-wave velocity ratio time series include:

[0008] S1. Obtain the P-wave and S-wave time difference data of the earthquake;

[0009] S2. Perform linear relationship fitting on the P-wave and S-wave time difference data, and extract the average wave velocity ratio of the earthquake pairs;

[0010] S3. Construct an initial PS wave velocity ratio time series based on the earthquake sequence;

[0011] S4. Using the average wave velocity ratio as the observation data and the initial PS wave velocity ratio time series as the unknown, construct a time series inversion matrix, and solve the time series inversion matrix to obtain the wave velocity ratio time series.

[0012] S5. Determine whether the wave speed ratio time series is stable. If the wave speed ratio time series is stable, increment the time series unknowns and return to S3. If the wave speed ratio time series is unstable, output the inverted wave speed ratio time series.

[0013] Optionally, acquiring P-wave and S-wave time difference data for earthquakes includes:

[0014] Acquire the arrival time data of the P-wave and S-wave of the earthquake and pair the earthquakes together;

[0015] Acquire P-wave to time difference data and S-wave to time difference data for earthquake pairs.

[0016] Optionally, the average wave velocity ratio of the seismic pairs can be extracted, including:

[0017] The least squares method was used to fit the P-wave and S-wave time difference data of the earthquake pair to obtain the linear relationship between the P-wave and S-wave time difference data.

[0018] The average wave velocity ratio of the earthquake pairs is extracted based on the slope of the linear relationship.

[0019] Optionally, constructing the initial PS wave velocity ratio time series includes: constructing the initial PS wave velocity ratio time series at equal time intervals based on the initial and final times of the earthquake sequence, or constructing the initial PS wave velocity ratio time series based on the time density of the earthquake sequence.

[0020] Optionally, the time series inversion matrix includes:

[0021] d = Gm

[0022]

[0023] in, Let m be the average wave velocity ratio measured for all m earthquakes with respect to slope, where m = [r1, r2, r3, ... r n Let ] be n time series Vp / Vs to be determined, G be an m*n sensitivity kernel coefficient matrix, and w ij Let t represent the sensitivity of the i-th earthquake pair to the Vp / Vs at time j. j and t j+1 Let t represent the j-th time and the (j+1)-th time respectively. e Let represent any one of the two occurrence times of the i-th earthquake pair.

[0024] Optionally, determining whether the wave velocity ratio time series is stable includes: comparing the number of unknowns in the time series with the model's second norm ||m|| to determine whether the wave velocity ratio time series is stable.

[0025] Optionally, outputting the inverted wave velocity ratio time series includes: obtaining wave velocity ratio time series with a number of time points less than the number of unknowns in the time series, performing an envelope operation, and outputting the inverted wave velocity ratio time series.

[0026] The beneficial effects of the present invention are: (1) Improved time resolution: By directly utilizing the time difference data of P-wave and S-wave arrival of earthquake pairs, the path calculation and earthquake occurrence time calculation errors in the traditional tomographic travel time calculation process are avoided, and the independent inversion of wave velocity ratio time series is achieved.

[0027] (2) Reduce inversion uncertainty: By gradually increasing the number of unknowns through iterative algorithms, the fit is steadily improved while maintaining inversion stability, eliminating the influence of seismic location error on time-varying wave velocity ratio structure, and improving the reliability of dynamic monitoring.

[0028] (3) Realize dynamic monitoring: Provide a high-precision inversion method suitable for long-term evolution analysis of crustal medium parameters (such as fluid migration, stress adjustment, etc.), and provide new technical means for earthquake precursor monitoring and time-varying research of underground structures.

[0029] This invention does not rely on the global inversion framework of traditional tomography, but focuses on refined modeling and inversion in the time dimension, thereby revealing the dynamic change law of wave velocity ratio more accurately. Attached Figure Description

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

[0031] Figure 1 This is a flowchart of a method for retrieving high-resolution P- and S-wave velocity ratio time series based on seismic wave arrival time, according to an embodiment of the present invention.

[0032] Figure 2 This is a schematic diagram of the method for retrieving high-resolution P-wave and S-wave velocity ratio time series based on seismic wave arrival time in an embodiment of the present invention.

[0033] Figure 3 This is a line graph of the P and S wave arrival time difference based on the seismic pair, according to an embodiment of the present invention.

[0034] Figure 4 The figures show the growth curves of the L2 norm and Vp / Vs time series parameters of the time series model in this embodiment of the invention.

[0035] Figure 5 The time series inversion results and envelope results are selected at different times for different parameters in the embodiments of the present invention;

[0036] Figure 6 This is a comparison chart of seismic activity and Vp / Vs time series of the Noto Peninsula according to an embodiment of the present invention;

[0037] Figure 7 The wave velocity ratio time series obtained by the method inversion in this embodiment of the invention and the time series directly obtained by the traditional double-difference wave velocity ratio method are shown. Detailed Implementation

[0038] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0039] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0040] like Figure 1 As shown, this embodiment provides a method for retrieving high-resolution P- and S-wave velocity ratio time series based on seismic wave arrival time inversion, including:

[0041] S1. Obtain the P-wave and S-wave time difference data of the earthquake;

[0042] S2. For each pair of earthquakes, the P-wave arrival time difference of the same station is used as the abscissa and the S-wave arrival time difference is used as the ordinate. A linear relationship is fitted to the P-wave and S-wave arrival time difference data of a series of stations, and the slope of the fitted line is extracted as the average wave velocity ratio of the earthquake pair.

[0043] S3. Based on the earthquake sequence, construct an initial PS wave velocity ratio time series containing n unknowns;

[0044] S4. Using the average wave velocity ratio of each pair of earthquakes as the observation data, the initial wave velocity ratio time series is constructed as the unknown. A time series inversion matrix is ​​constructed, and the time series inversion matrix is ​​inverted and solved to obtain the wave velocity ratio time series.

[0045] S5. Determine whether the wave speed ratio time series is stable. If the wave speed ratio time series is stable, save the series and increment the number of unknowns n in the time series, then return to S3. If the wave speed ratio time series is unstable, output the wave speed ratio time series inverted through the above steps.

[0046] Specifically, the method in this embodiment has high temporal resolution: it directly inverts time difference parameters, avoiding the grid discretization error of traditional tomography. It also has strong anti-interference capabilities: by utilizing seismic pulse time difference data, it reduces the impact of seismic location errors and path uncertainties. Furthermore, it has wide applicability: it is suitable for dynamic detection scenarios of underground structures, such as volcano monitoring, fluid migration analysis, earthquake precursor research, and monitoring of artificially blasted underground structures.

[0047] Furthermore, acquiring the P-wave and S-wave time difference data of the earthquake includes:

[0048] Acquire the arrival time data of the P-wave and S-wave of the earthquake and pair the earthquakes in pairs; during the acquisition, ensure that the distance between earthquakes is small, and that the distance between the station and the earthquake is greater than the distance between the earthquake pairs;

[0049] Acquire P-wave to time difference data and S-wave to time difference data for earthquake pairs.

[0050] Specifically, dense array observation requires that the stations be distributed densely enough to ensure that the P-wave and S-wave of each seismic event are recorded by at least 3 stations.

[0051] Seismic phase catalog: The origin time of each earthquake and the accurate arrival times of P-waves and S-waves recorded by each station need to be provided. Two earthquakes recorded by the same seismic station are collected to obtain earthquake pairs.

[0052] Furthermore, extracting the average wave velocity ratio of earthquake pairs includes:

[0053] Using the least squares method to fit the linear relationship of the arrival time difference data of P-waves and S-waves to obtain the linear relationship of the arrival time difference data of P-waves and S-waves of earthquake pairs;

[0054] According to the slope of the linear relationship, the average wave velocity ratio of earthquake pairs is extracted.

[0055] Furthermore, constructing the initial wave velocity ratio time series includes: constructing the initial P-S wave velocity ratio time series at equal time intervals according to the initial and end times of the earthquake sequence, or constructing the initial P-S wave velocity ratio time series according to the time density of the earthquake sequence.

[0056] Furthermore, the time series inversion matrix includes:

[0057] d = Gm

[0058] where, is the average wave velocity ratio measured for all m earthquake pairs, m = [r1, r2, r3,... r n is the Vp / Vs time series of n unknowns to be solved, and G is the sensitivity kernel coefficient matrix of m*n.

[0059] Specifically, for the Vp / Vs measured for an earthquake pair, an above equation can be obtained. For a long-term earthquake sequence composed of a large number of earthquakes, a large number of wave velocity ratios can be measured, and a system of linear equations can be formed from the above equations: d = Gm. G is the coefficient matrix of size m*n, representing the sensitivity of the mth earthquake pair to the Vp / Vs at the nth moment. To maintain the stability of the solution process, usually the number of parameters of the Vp / Vs time series is much smaller than the number of measured average wave velocity ratios, that is, n < m. Then, using the existing overdetermined equation inversion tool, the solution of this system of equations can be realized, and the time series of Vp / Vs can be obtained.

[0060] Furthermore, judging whether the wave velocity ratio time series is stable includes: comparing the number of unknowns in the time series with the model two-norm ||m|| to judge whether the wave velocity ratio time series is stable.

[0061] Furthermore, outputting the inverted wave velocity ratio time series includes: obtaining the wave velocity ratio time series with the number of moments less than the number of unknowns in the time series and performing envelope processing to output the inverted wave velocity ratio time series.

[0062] The following combines with the appendix Figure 2The method of this embodiment will be further explained as follows:

[0063] This embodiment provides a time series inversion method based on the arrival time differences of P-waves and S-waves in earthquake pairs. By analyzing the arrival time differences of P-waves and S-waves in earthquake pairs, a linear regression model is established, and a time series inversion matrix is ​​constructed. Finally, a high-resolution VP / VS time evolution sequence is obtained. This method overcomes the problems of static inversion results and dependence on the accuracy of earthquake location in traditional tomography techniques, and can more accurately characterize the dynamic changes of the subsurface medium.

[0064] 1. Data premise:

[0065] Dense array observation: requires stations to be distributed densely enough to ensure that the P-wave and S-wave of each seismic event are recorded by at least 3 stations.

[0066] Earthquake phase catalog: The time of occurrence of each earthquake and the precise hourly data of P-waves and S-waves recorded by each station must be provided.

[0067] 2. Technical Principles:

[0068] The arrival time differences of P-waves and S-waves in an earthquake pair satisfy the following linear relationship:

[0069]

[0070] In the formula T represents the arrival times of the P and S waves emitted by earthquakes A and B, respectively. AB This represents the time difference between the occurrence times of earthquakes A and B. The slope of the linear regression is the average Vp / Vs of the earthquake on the corresponding medium, reflecting the mean wave velocity ratio between the occurrence times of the two earthquakes.

[0071] This formula shows that by obtaining the arrival times of the P-wave and S-wave between two earthquakes, the wave velocity ratio can be fitted using the slope, such as... Figure 2 As shown.

[0072] When the wave velocity ratio in the subsurface medium does not change over time, the slope is a constant. However, when the wave velocity ratio in the subsurface medium changes over time and the stations are located on opposite sides of the earthquake pair, the slope of the fitted straight line represents the average wave velocity ratio at the time of the two earthquakes. Therefore, this characteristic can be used to invert the time series of Vp / Vs from long-term earthquake sequences, such as... Figure 2 As shown.

[0073] 3. Calculate the average Vp / Vs using linear regression:

[0074] For each earthquake pair, the least squares method is used for fitting. and The linear relationship, the slope of which is the ratio of the observed average wave velocity r of the earthquake pair.AB .

[0075] 4. Inversion stage:

[0076] The inversion principle is as follows:

[0077] Based on the earthquake occurrence time, the Vp / Vs time series m = (r1, r2, ... r) is generated. n There are two methods for constructing the time series. The first method constructs n time points based on the equal time intervals between the initial and final moments of the earthquake sequence: t1, t1+Δt, t1+2Δt, ..., t1+nΔt. The second method constructs the time points of the time series based on the time density of the earthquake sequence. Although the time intervals constructed in this way are not equal, it is beneficial to the stability inversion process. First, the occurrence times of all earthquakes are arranged, and the occurrence times of every m / n earthquakes are used as the time parameters of the Vp / Vs time series: t1, t2, ..., t3, ..., t4, ..., t5, ..., t6, ..., t7, ..., t8, ..., t9, ..., t1 ... 1+m / n ,t 1+2m / n ,t 1+3m / n ...t 1+m .

[0078] Using the average Vp / Vs as the observed data, and the constructed Vp / Vs time series as the unknown, a coefficient matrix is ​​established between the two.

[0079] Let the difference parameters of the time series be r1, r2, ..., r n , representing the change in Vp / Vs at different times.

[0080] The average wave velocity ratio r of earthquake to AB AB This can be expressed as a weighted sum of the difference parameters at adjacent time points:

[0081]

[0082] For a single earthquake pair with measured Vp / Vs, the above equation can be obtained. For a long-term earthquake sequence consisting of a large number of earthquakes, a large number of wave velocity ratios can be measured. The above equations can be combined to form a system of linear equations, i.e., the time series inversion matrix:

[0083] d = Gm

[0084] in Let m be the average wave velocity ratio measured for all m earthquakes with respect to slope, where m = [r1, r2, r3, ... r n Let G be n time series Vp / Vs to be determined. G is a coefficient matrix of size m*n, where the coefficient w in the i-th row and j-th column is... ij Let Vp / Vs represent the sensitivity of the i-th earthquake pair to the j-th time interval, where i = 1, 2, 3, ..., m, j = 1, 2, 3, ..., n. The expression is:

[0085]

[0086] where t j and t j+1 represent the j-th moment and the (j + 1)-th moment respectively, and t e represents any one of the two earthquake occurrence times of the i-th pair of earthquakes.

[0087] To maintain the stability of the solution process, usually the number of parameters of the Vp / Vs time series is much smaller than the number of measured average wave velocity ratios, that is, n < m. Then, using the existing overdetermined equation inversion tool, the solution of this system of equations can be realized, and the time series of Vp / Vs can be obtained.

[0088] 5. Inversion to solve the time series:

[0089] Adopt optimization algorithms such as LSQR, LSMR or BFGS to solve the linear equation system d = Gm, and obtain high-precision Vp / Vs time series r1, r2,..., r n .

[0090] Optimization algorithms such as LSQR need to set damping parameters. However, when the number of unknowns is much smaller than the number of equations, the damping can be set to 0. Therefore, the inversion result of this technology does not depend on parameter adjustment.

[0091] 6. Evaluate the uncertainty of the time series:

[0092] Since the set time interval and the number of inversion unknowns will greatly affect the time resolution of the result, in this stage, the number of time series unknowns n is gradually increased from 1 to m, and then according to each different n, the 4th and 5th stages are repeated in turn. Finally, m different time series are obtained. Then analyze the growth curve between the number of time series unknowns n and the model two-norm. If the growth is slow and smooth, it is judged that the wave velocity ratio time series is stable and no overfitting occurs; if the growth is steep or extreme values appear, it is judged that the time series is unstable (i.e., overfitting). Then randomly select the decrease values of the model two-norm and the inversion residual of these time series in turn. As n increases, the originally smooth and gentle curve will have inflection points or even extreme points. The n value corresponding to the first inflection point is the optimal number of unknowns, that is, the situation with the best time resolution and stable inversion, as Figure 4 shown. The envelope of the n time series obtained by the above inversion can be used as the uncertainty of the Vp / Vs time series, as Figure 5 shown, the time series inversion results and envelope results for different parameter selections.

[0093] ]>The following further illustrates the method of this embodiment with the Noto Peninsula of a certain country as the research area:

[0094] Methods for inverting time series of wave velocity ratios from seismic P-waves and S-waves based on time differences include:

[0095] 1. Data collection and preprocessing:

[0096] Study area: Noto Peninsula, a country (August 2021 to October 2024), during which a series of moderate to strong earthquakes occurred, including the 7.6 magnitude mainshock on January 1, 2024 and subsequent aftershocks.

[0097] Data source: Earthquake catalog provided by the Japan Meteorological Agency (JMA), including: earthquake longitude, latitude, depth, and time of occurrence t0 (UTC time); arrival times t0 of P-waves and S-waves recorded by each station. P , t S Each earthquake must be recorded by at least three stations to ensure the stability of subsequent regression calculations; station coordinate information.

[0098] Data filtering: Events with low signal-to-noise ratios or unreliable arrival times were removed. Earthquakes less than 10 km from the mainshock were selected to ensure that the distance between earthquakes was small enough that the spatial variation of the Vp / Vs wave velocity ratio was negligible. Finally, 2,843 earthquakes and 34,482 sets of P-wave and S-wave arrival times were retained to form the basic dataset for inversion.

[0099] 2. Earthquake impact on compilation and wave velocity ratio calculation:

[0100] Earthquake generation rules:

[0101] Two earthquakes must be recorded by at least 10 identical stations;

[0102] The arrival times of the P-waves and S-waves of both earthquakes were recorded by these stations.

[0103] Calculation process:

[0104] (1) Select two earthquakes recorded by the same station to form an earthquake pair (e.g., earthquake A and earthquake B). For each earthquake pair (e.g., earthquake A and earthquake B), calculate the P-wave arrival time difference (Δt) for each station. p =t_ PB -t_ PA ) and S-wave arrival time difference (Δt) S =t_ SB -t_ SA ).

[0105] (2) Fitting Δt using the least squares method S With △t p The linear relationship is represented by the slope, which is the average Vp / Vs wave velocity ratio (r_AB) of the earthquake pair. Figure 3This paper presents the linear regression and slope fitting process for two earthquake pairs on the Noto Peninsula. Earthquake A occurred on January 4, 2024, with a magnitude of 2.5, and Earthquake B occurred on January 3, 2024, with a magnitude of 2.6. Their P and S waves were received by 16 surrounding stations, and their locations are shown in the upper left subfigure. The linear fitting of the arrival time differences of the P and S waves of earthquakes AB recorded by these 16 stations is shown in the upper right subfigure. Based on the least squares method for slope fitting, the average wave velocity ratio of this earthquake pair is 1.68. The locations of earthquakes C and D are shown in the lower left subfigure. The linear fitting of the arrival time differences of the P and S waves of earthquake CD is shown in the lower right subfigure.

[0106] (3) To improve the accuracy of linear regression, if the root mean square residual of the fit is greater than 0.2 seconds or the slope uncertainty of the fit is greater than 0.1, the earthquake pair should be removed.

[0107] Results: A total of 12,517 valid earthquake pairs were generated, covering the entire study period.

[0108] 3. Construction and solution of time series inversion matrix:

[0109] Inversion equation establishment:

[0110] Discretize the time axis into daily intervals, and let the daily changes in Vp / Vs be unknowns r1, r2, ... r n .

[0111] Vp / Vs observations for each earthquake pair r AB The parameter relationship with adjacent time nodes is as follows:

[0112]

[0113] By combining the equations of all earthquake pairs, an overdetermined linear equation system d = Gm is formed (matrix G dimension: 12,517 x 1,095, corresponding to 3 years of data).

[0114] Inversion algorithm: The LSQR algorithm is used to solve the Vp / Vs time series, and the regularization damping parameter is set to 0.

[0115] By continuously increasing the number of unknowns, new Vp / Vs time series are solved again until the model overfits. All stable Vp / Vs time series models are output and their envelopes are calculated to obtain the final Vp / Vs diurnal variation curves for the Noto Peninsula.

[0116] 4. Inversion Results and Verification:

[0117] Vp / Vs time series:

[0118] Inversion results show that the Vp / Vs ratio in the Noto Peninsula source area exhibited multiple short-term anomalous fluctuations (with a 2-5 fold increase in power spectral density in bands with a period of 10-80 days), corresponding to deep underground fluid upwelling events. Figure 6 The graph shows the seismic activity and Vp / Vs time series comparison of the Noto Peninsula. The first row represents the number of earthquakes recorded every 10 days. The second row represents the focal depth, with E1-E5 representing 5 lower crustal seismic events. The third row represents the Vp / Vs time series, with the shaded areas above and below the curve representing the envelope. The fourth row represents the fluctuation of the Vp / Vs time series at different periods, with the vertical axis in days.

[0119] The timing of the abnormal fluctuations was highly consistent with the period of increased activity of deep crustal earthquakes (depth > 25 km) (correlation coefficient 0.89), indicating that the earthquakes in the lower crust opened up fluid migration channels and promoted the migration of fluids in the seismic zone.

[0120] Key event relevance:

[0121] April 2023: Vp / Vs fluctuated significantly, and a 6.5 magnitude foreshock occurred one month later;

[0122] December 2023: Vp / Vs fluctuated significantly, and a 7.6 magnitude main shock and a 6.1 magnitude aftershock occurred a month later.

[0123] 5. Method reliability verification:

[0124] Synthetic data testing: In synthetic datasets with known Vp / Vs variation patterns, the time series recovered by this method has an error of less than 2% compared to the true model.

[0125] Compared with traditional methods: such as Figure 7 As shown, the wave velocity ratio time series obtained by inversion using this embodiment is compared with the time series directly obtained by the traditional double-difference wave velocity ratio method. The traditional double-difference wave velocity ratio method requires that the occurrence time of each pair of earthquakes be relatively close, so a large number of earthquake pairs are ignored in the results. However, this embodiment uses earthquake pairs with a longer time interval, so the time resolution is significantly improved, and the diurnal variation dynamic process of fluid activity is successfully revealed.

[0126] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for inverting high-resolution P and S wave velocity ratio time series based on seismic travel times, characterized in that, The method comprises the following steps: S1, obtaining P-wave and S-wave travel time difference data of earthquakes; S2, performing linear relationship fitting on the P-wave and S-wave travel time difference data to extract average wave velocity ratio of earthquake pairs; The step of extracting average wave velocity ratio of earthquake pairs comprises: performing linear relationship fitting on the P-wave and S-wave travel time difference data of the earthquake pairs by using least square method to obtain linear relationship of P-wave and S-wave travel time difference data; extracting average wave velocity ratio of the earthquake pairs according to the slope of the linear relationship; S3, constructing an initial PS wave velocity ratio time sequence according to a seismic sequence; The step of constructing an initial PS wave velocity ratio time sequence comprises: constructing the initial PS wave velocity ratio time sequence at equal time intervals according to initial and end time of the seismic sequence, or constructing the initial PS wave velocity ratio time sequence according to time density of the seismic sequence; S4, constructing a time sequence inversion matrix by taking the average wave velocity ratio as observation data and taking the initial PS wave velocity ratio time sequence as unknown quantity, and performing inversion solving on the time sequence inversion matrix to obtain a wave velocity ratio time sequence; The time sequence inversion matrix comprises: d=Gm; ; Where d=[ , [, ...] represents the average wave velocity ratio measured for all m earthquake pairs with slope, where m = [r1, r2, r3, ... r n Let ] be n time series Vp / Vs to be determined, G be an m*n sensitivity kernel coefficient matrix, and w ij Let t represent the sensitivity of the i-th earthquake pair to the Vp / Vs at time j. j and t j+1 Let t represent the j-th time and the (j+1)-th time respectively. e Represents any one of the two occurrence times of the i-th earthquake; S5, judging whether the wave velocity ratio time sequence is stable, if the wave velocity ratio time sequence is stable, time sequence unknown quantity plus one returns to S3, if the wave velocity ratio time sequence is not stable, the inversion wave velocity ratio time sequence is outputted; The step of outputting the inversion wave velocity ratio time sequence comprises: obtaining a wave velocity ratio time sequence with a time point number less than the number of time sequence unknown quantities and performing envelope on the wave velocity ratio time sequence, and outputting the inversion wave velocity ratio time sequence.

2. The method for high resolution P and S wave velocity ratio time series inversion based on seismic travel times of claim 1, wherein, The step of obtaining P-wave and S-wave travel time difference data of earthquakes comprises: obtaining P-wave and S-wave travel time data of earthquakes and pairing the earthquakes two by two; obtaining P-wave travel time difference data and S-wave travel time difference data of earthquake pairs.

3. The method for high resolution P and S wave velocity ratio time series inversion based on seismic travel times of claim 1, wherein, The judging whether the wave velocity ratio time series is stable comprises: comparing the number of time series unknowns with the model two norm The judging whether the wave velocity ratio time series is stable comprises: comparing the number of time series unknowns with the model two norm

Images