Method for inverting high-resolution P and S wave velocity ratio time sequence based on seismic wave arrival time

Abstract

The invention relates to a method for inverting a high-resolution P and S wave velocity ratio time sequence based on seismic wave arrival time. The method comprises the following steps: S1, acquiring longitudinal wave and transverse wave arrival time difference data of a seismic pair; s2, carrying out linear relation fitting on longitudinal wave and transverse wave arrival time difference data, and extracting an average wave velocity ratio of an earthquake pair; s3, constructing an initial wave velocity ratio time sequence according to the earthquake generating moment of the earthquake sequence; s4, constructing a time sequence inversion matrix, and performing inversion solution on the time sequence inversion matrix to obtain a wave velocity ratio time sequence; and S5, judging whether the wave velocity ratio time sequence is stable or not, if the wave velocity ratio time sequence is stable, adding one to the unknown number of the time sequence, returning to S3, and if the wave velocity ratio time sequence is unstable, outputting the inverted wave velocity ratio time sequence. According to the method, the dynamic change rule of the wave velocity ratio of a natural seismic zone or an artificial blasting area can be more accurately revealed.

Description

Technical Field

[0001] The invention relates to the technical field of underground structure inversion, 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 Art

[0002] In geophysical exploration and seismological research, inverting the velocity structure of underground media is of great significance for understanding crustal deformation, fault activity and earthquake generation mechanisms. Among them, the velocity ratio of longitudinal waves (P waves) and transverse waves (S waves) (Vp / Vs; hereinafter referred to as the velocity ratio) is particularly sensitive to processes such as underground fluid migration and pore pressure changes. At present, the commonly used Vp / Vs velocity ratio tomography technology has achieved remarkable results in static spatial structure inversion, but it still has great limitations in time resolution and dynamic monitoring. There are two main existing wave velocity ratio dynamic monitoring technologies, including the Heda method and the double difference wave velocity ratio method. The Heda method uses the longitudinal wave arrival times t of multiple seismic events recorded by the station. P and the time difference between longitudinal and transverse waves (t S -t P ), add one to the slope of the fit to get the wave velocity ratio, that is, Vp / Vs=(t S -t P ) / t P +1. This method mainly reflects the average wave velocity ratio of the underground medium between the earthquake and the station, and cannot reflect the wave velocity ratio of a local area underground, so the spatial resolution is insufficient. The double difference wave velocity ratio method mainly uses the longitudinal wave arrival time difference Δt of the same pair of earthquakes recorded by multiple stations. P and the time difference Δt S , the wave velocity ratio is obtained from the slope of the fitting line, that is, Vp / Vs=Δt S / Δt P . This method mainly reflects the wave velocity ratio of the local area between earthquake pairs, so the spatial resolution is improved. However, since the occurrence time of earthquake pairs may be far apart, and similar earthquake events cannot be received by many stations at the same time, the time resolution and fitting error of this method are greatly affected by the original data. In short, the traditional wave velocity ratio underground structure acquisition method has the defects of low time resolution, dependence on earthquake location accuracy, and large error in wave velocity ratio structure calculation.

[0003] The prior art discloses a method for real-time acquisition of microseismic wave velocity in hard rock tunnels constructed by drilling and blasting. The method is mainly based on the Heda method, which is used for real-time acquisition of wave velocity and relies on earthquake positioning accuracy. However, it does not have the ability to resolve changes in wave velocity over time.

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

[0005] The object of the present invention is to provide a method for inverting high-resolution P- and S-wave velocity ratio time series based on the arrival times of seismic waves, focusing on refined modeling inversion in the time dimension, so as to more accurately reveal the dynamic variation law of the wave velocity ratio.

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

[0007] A method for inverting high-resolution P- and S-wave velocity ratio time series based on the arrival times of seismic waves includes:

[0008] S1. Obtain the arrival time difference data of P-waves and S-waves of earthquakes;

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

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

[0011] S4. Use the average wave velocity ratio as the observation data and the initial PS wave velocity ratio time series as the unknowns, construct a time series inversion matrix, and perform inversion solution on the time series inversion matrix to obtain the wave velocity ratio time series;

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

[0013] Optionally, obtaining the arrival time difference data of P-waves and S-waves of earthquakes includes:

[0014] Obtain the arrival time data of P-waves and S-waves of earthquakes and pair the earthquakes in pairs;

[0015] Obtain the arrival time difference data of P-waves and the arrival time difference data of S-waves of earthquake pairs.

[0016] Optionally, extracting the average wave velocity ratio of earthquake pairs includes:

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

[0018] Extract the average wave velocity ratio of the earthquake pairs according to 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 according to the initial and end moments of the earthquake sequence, or constructing the initial PS wave velocity ratio time series according to the time density of the earthquake sequence.

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

[0021] d = Gm

[0022]

[0023] where is the average wave velocity ratio measured for the slopes of all m earthquake pairs, m = [r1, r2, r3,... r n is the Vp / Vs time series of n unknowns to be solved, G is the sensitivity kernel coefficient matrix of m * n, w ij is the sensitivity of Vp / Vs at the j-th moment of the i-th earthquake pair, t j and t j+1 represent the j-th moment and the j + 1-th moment respectively, t e represents any one of the two earthquake 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 two-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 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.

[0026] The beneficial effects of the present invention are: (1) improving time resolution: by directly using the arrival time difference data of P waves and S waves of earthquake pairs, avoiding the path calculation and calculation errors of earthquake occurrence times in the traditional tomography travel time calculation process, and realizing the independent inversion of the wave velocity ratio time series.

[0027] (2) Reducing inversion uncertainty: by gradually increasing the number of unknowns through an iterative algorithm, steadily improving the fitting degree while maintaining inversion stability, eliminating the influence of earthquake location errors on the time-varying structure of the wave velocity ratio, and improving the reliability of dynamic monitoring.

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

[0029] Rather than relying on the global inversion framework of traditional tomography, the present invention focuses on refined modeling inversion in the time dimension, thereby more accurately revealing the dynamic variation law of the wave velocity ratio. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

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

[0032] Figure 2 It is a schematic diagram of the method for inverting the high-resolution P and S wave velocity ratio time series based on seismic wave arrival times according to an embodiment of the present invention;

[0033] Figure 3 It is a fitting straight line diagram of the P and S wave arrival time differences of the seismic pairs according to an embodiment of the present invention;

[0034] Figure 4 It is a growth curve of the two-norm of the time series model and the parameters of the Vp / Vs time series according to an embodiment of the present invention;

[0035] Figure 5 It is the time series inversion result and envelope result with different parameter selections according to an embodiment of the present invention;

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

[0037] Figure 7 It is the wave velocity ratio time series obtained by inverting the method according to an embodiment of the present invention and the time series directly obtained by using the traditional double-difference wave velocity ratio method. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0038] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the scope of protection of the present invention.

[0039] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below in conjunction with the drawings and specific embodiments.

[0040] As shown Figure 1 in the figure, this embodiment provides a method for inverting the high-resolution P and S wave velocity ratio time series based on the arrival times of seismic waves, including:

[0041] S1. Obtain the arrival time difference data of the longitudinal wave and the transverse wave of the earthquake;

[0042] S2. For each pair of earthquakes, take the longitudinal wave arrival time difference at the same station as the abscissa and the transverse wave arrival time difference as the ordinate, perform a linear relationship fitting on the longitudinal and transverse wave arrival time difference data of a series of stations, and extract the slope of the fitting line as the average wave velocity ratio of the earthquake pair;

[0043] S3. According to 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 observed data and the constructed initial wave velocity ratio time series as the unknowns, construct a time series inversion matrix, and perform an inversion solution on the time series inversion matrix to obtain the wave velocity ratio time series;

[0045] S5. Judge whether the wave velocity ratio time series is stable. If the wave velocity ratio time series is stable, save the series, and at the same time increment the number n of the time series unknowns by one and return to S3. If the wave velocity ratio time series is unstable, output the wave velocity ratio time series inverted through the above steps.

[0046] Specifically, the method of this embodiment has high time resolution: directly invert the time difference parameters to avoid the grid discretization error of traditional tomography. It has strong anti-interference ability: utilize the arrival time difference data of earthquake pairs to reduce the influence of earthquake location error and path uncertainty. It has wide applicability: applicable to underground structure dynamic detection scenarios such as volcanic monitoring, fluid migration analysis, earthquake precursor research, and underground structure monitoring of artificial blasting.

[0047] Furthermore, obtaining the arrival time difference data of the longitudinal wave and the transverse wave of the earthquake includes:

[0048] Obtain the arrival time data of the longitudinal wave and the transverse wave of the earthquake and pair the earthquakes two by two; ensure that the distance between earthquakes is small during acquisition, and at the same time the distance between the station and the earthquake is greater than the distance between earthquake pairs;

[0049] Obtain the longitudinal wave arrival time difference data and the transverse wave arrival time difference data of the earthquake pair.

[0050] Specifically, for dense array observation: it is required that the stations are distributed densely enough to ensure that the P wave and S wave of each earthquake 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 the earthquake pair;

[0054] According to the slope of the linear relationship, the average wave velocity ratio of the earthquake pair 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 from the slopes of 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 linear equation system can be formed by 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 equation system 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 the appendix Figure 2A further description of the method of this embodiment is as follows:

[0063] This embodiment provides a method for inverting the wave velocity ratio time series based on the time difference between the arrival times of P-waves and S-waves in earthquakes. By analyzing the time difference between the arrival times of P-waves and S-waves of earthquake pairs, a linear regression model is established, and a time series inversion matrix is constructed. Finally, a high-resolution VP / VS time evolution series is solved. This method overcomes the problems that the inversion results of traditional tomography techniques are static and rely on the accuracy of earthquake location, and can more accurately characterize the dynamic changes of underground media.

[0064] 1. Data prerequisites:

[0065] Dense array observation: It is required that the station distribution is dense enough to ensure that the P-waves and S-waves of each earthquake event are recorded by at least 3 stations.

[0066] Earthquake phase catalog: The origin time of each earthquake and the accurate arrival time data of P-waves and S-waves recorded by each station need to be provided.

[0067] 2. Technical principle:

[0068] The time difference between the arrival times of P-waves and S-waves of earthquake pairs satisfies the following linear relationship:

[0069]

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

[0071] This formula shows that as long as the time differences between the arrival times of P-waves and S-waves between two earthquakes are obtained, the wave velocity ratio can be fitted by the slope, as Figure 2 shown.

[0072] When the wave velocity ratio of the underground medium does not change with time, then the slope is a constant value. However, when the wave velocity ratio of the underground medium changes with time and the stations are distributed on both sides of the earthquake pair, the slope of the fitted straight line represents the average wave velocity ratio between the occurrence times of the two earthquakes. Therefore, this feature can be utilized to invert the time series of Vp / Vs using a long-term earthquake sequence, as Figure 2 shown.

[0073] 3. Obtaining the average Vp / Vs by linear regression:

[0074] For each earthquake pair, use the least squares method to fit the linear relationship between and The slope is the observed average wave velocity ratio r of this earthquake pairAB 。

[0075] 4. Inversion stage:

[0076] The inversion principle is as follows:

[0077] Construct the Vp / Vs time series m = (r1, r2,... r n ) according to the occurrence times of the earthquake sequence. There are two methods. The first is to construct n time points at equal time intervals based on the initial and end times of the earthquake sequence: t1, t1+Δt, t1+2Δt,... t1+nΔt. The second is to construct the time points of the time series according to the time density of the earthquake sequence. Although the time intervals constructed in this way are not equal, it is beneficial to stabilize the inversion process. First, arrange the occurrence times of all earthquakes, and use the occurrence times of every m / n earthquakes as the time parameters of the Vp / Vs time series: t1, t 1+m / n , t 1+2m / n , t 1+3m / n ... t 1+m 。

[0078] Take the average Vp / Vs as the observed data, the constructed Vp / Vs time series as the unknowns, and establish the coefficient matrix 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 the earthquake to A - B AB can be expressed as the weighted sum of the difference parameters at adjacent times:

[0081]

[0082] For the Vp / Vs measured for an earthquake pair, an above - mentioned 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 linear equation system can be formed from the above - mentioned equations, that is, the time series inversion matrix:

[0083] d = Gm

[0084] where is the average wave velocity ratio measured from the slopes of all m earthquake pairs, m = [r1, r2, r3,... r n is the n Vp / Vs time series to be solved. G is the coefficient matrix of size m*n, where the coefficient w ij in the i - th row and j - th column represents the sensitivity of the Vp / Vs of the i - th earthquake pair at the j - th time, i = 1, 2, 3... m, j = 1, 2, 3... n, and 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 and solution of the time series:

[0089] Adopt optimization algorithms such as LSQR, LSMR or BFGS to solve the linear equation system d = Gm to 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 there is no overfitting; if the growth is steep or there are extreme values, 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 case with the best time resolution and relatively 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 by taking the Noto Peninsula of a certain country as the research area:

[0094] A method for inverting the time series of the wave velocity ratio based on the time difference between the arrival times of P-waves and S-waves in an earthquake, including:

[0095] 1. Data collection and preprocessing:

[0096] Study area: The Noto Peninsula in a certain country (from August 2021 to October 2024). A series of moderate to strong earthquakes occurred in this area during this period, including the magnitude 7.6 main earthquake 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 origin time t0 (UTC time); the arrival times of P-waves and S-waves recorded at each station t P t S , t

[0098] Data screening: Exclude events with low signal-to-noise ratio or unreliable arrival time markings, and screen out earthquakes with a distance less than 10 kilometers from the main earthquake to ensure that the distance between earthquakes is small enough and the spatial variation of the Vp / Vs wave velocity ratio can be ignored. Finally, 2,843 earthquakes and 34,482 groups of P-wave and S-wave arrival times are retained to form the inversion basic data set.

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

[0100] Rules for generating earthquake pairs:

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

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

[0103] Calculation process:

[0104] (1) Screen two earthquakes recorded by the same stations to form an earthquake pair (such as earthquake A and earthquake B). For each earthquake pair (such as earthquake A and earthquake B), calculate the P-wave arrival time difference (△t p =t_ PB -t_ PA ) and S-wave arrival time difference (△t S =t_ SB -t_ SA ) at each station.

[0105] (2) Use the least squares method to fit the linear relationship between △t S and △t p . The slope is the average Vp / Vs wave velocity ratio (r_AB) of this earthquake pair. Figure 3The linear regression and slope fitting process of two pairs of earthquakes on the Noto Peninsula are shown. 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 P and S wave arrival time differences of earthquake AB recorded by the 16 stations is shown in the upper right subfigure. According to the least squares fitting slope, the average wave velocity ratio of the earthquake pair is 1.68. The locations of earthquakes C and D are shown in the lower left subfigure. The linear fitting of the P and S wave arrival time differences of the recorded earthquake CD is shown in the lower right subfigure.

[0106] (3) To improve the accuracy of linear regression, if the RMS residual of the fit is greater than 0.2 seconds or the uncertainty of the slope of the fit is greater than 0.1, the earthquake pair is 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] The inversion equation is established:

[0110] Discretize the time axis into daily intervals, and assume that the daily Vp / Vs change is the unknown number r1, r2, ...r n .

[0111] The observed value of Vp / Vs for each earthquake pair r AB The parameter relationship with adjacent time nodes is:

[0112]

[0113] The equations of all earthquake pairs are combined to form an overdetermined linear equation system d=Gm (matrix G dimension: 12,517x1,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] The number of unknowns is continuously increased, and the new Vp / Vs time series is re-solved until the model is overfitted. All stable Vp / Vs time series models are output and enveloped to obtain the final Vp / Vs daily variation curve of the Noto Peninsula.

[0116] 4. Inversion results and verification:

[0117] Vp / Vs time series:

[0118] The inversion results show that there are multiple short-term abnormal fluctuation events in the Vp / Vs of the source area of the Noto Peninsula (the power spectral density of the band with a period of 10 - 80 days increases by 2 - 5 times), corresponding to the upwelling events of deep underground fluids. As Figure 6 shown in the comparison chart of the seismicity and Vp / Vs time series of the Noto Peninsula. The first row shows the number of earthquakes recorded every 10 days. The second row represents the focal depth, and E1 - E5 represent 5 crustal earthquakes. The third row represents the Vp / Vs time series, and the shaded areas above and below the curve represent the envelope range. The fourth row represents the fluctuations of the Vp / Vs time series with different periods, and the unit of the vertical axis is days.

[0119] The time of abnormal fluctuations is highly consistent with the enhanced period of deep crustal earthquakes (depth > 25 km) (correlation coefficient 0.89), indicating that the deep crustal earthquakes have opened up fluid migration channels and promoted the migration of fluids in the seismic zone.

[0120] Relevance of key events:

[0121] April 2023: Significant fluctuations in Vp / Vs, and a 6.5 - magnitude foreshock occurred one month later;

[0122] December 2023: Significant fluctuations in Vp / Vs, and a 7.6 - magnitude main shock and a 6.1 - magnitude aftershock occurred one month later.

[0123] 5. Verification of method reliability:

[0124] Synthetic data test: In the synthetic dataset with known Vp / Vs variation patterns, the error between the time series recovered by this method and the true model is < 2%.

[0125] Comparison with traditional methods: As Figure 7 shown, the time series of wave velocity ratio obtained by inversion in this embodiment is compared with the time series directly obtained by the traditional double - difference wave velocity ratio method. Since the traditional double - difference wave velocity ratio method requires the earthquake origin times of each pair of earthquakes to be close, a large number of earthquake pairs are ignored in the results. However, in this embodiment, due to the use of earthquake pairs with a longer time interval, the time resolution is significantly improved, and the daily variation dynamic process of fluid activity is successfully revealed.

[0126] The embodiments described above are only descriptions of the preferred embodiments of the present invention, and do not limit the scope of the present invention. Without departing from the design spirit of the present invention, various deformations and improvements made by those of ordinary skill in the art to the technical solutions of the present invention shall fall within the protection scope determined 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 wave arrival times, characterized in that Including: S1. Obtain the arrival time difference data of the P-wave and S-wave of the earthquake; S2. Perform linear relationship fitting on the arrival time difference data of the P-wave and S-wave, and extract the average wave velocity ratio of the earthquake pair; S3. Construct an initial PS wave velocity ratio time series according to the earthquake sequence; S4. Take the average wave velocity ratio as the observed data and the initial PS wave velocity ratio time series as the unknowns, construct a time series inversion matrix, and perform inversion solution on the time series inversion matrix to obtain the wave velocity ratio time series; S5. Judge whether the wave velocity ratio time series is stable. If the wave velocity ratio time series is stable, increment the time series unknown by one and return to S3. If the wave velocity ratio time series is unstable, output the inverted wave velocity ratio time series.

2. The method for inverting the high-resolution P- and S-wave velocity ratio time series based on the seismic wave arrival time according to claim 1, characterized in that, Obtaining the arrival time difference data of the P-wave and S-wave of the earthquake includes: Obtain the arrival time data of the P-wave and S-wave of the earthquake and pair the earthquakes in pairs; Obtain the arrival time difference data of the P-wave and the arrival time difference data of the S-wave of the earthquake pair.

3. The method for inverting the time series of the high-resolution P- and S-wave velocity ratio based on the seismic wave arrival time according to claim 2, wherein, Extracting the average wave velocity ratio of the earthquake pair includes: Use the least squares method to perform linear relationship fitting on the arrival time difference data of the P-wave and S-wave of the earthquake pair to obtain the linear relationship of the arrival time difference data of the P-wave and S-wave; Extract the average wave velocity ratio of the earthquake pair according to the slope of the linear relationship.

4. The method for inverting the high-resolution P- and S-wave velocity ratio time series based on seismic wave arrival times according to claim 1, characterized in that, Constructing the initial PS wave velocity ratio time series includes: constructing the initial PS wave velocity ratio time series at equal time intervals according to the initial and end times of the earthquake sequence, or constructing the initial PS wave velocity ratio time series according to the time density of the earthquake sequence.

5. The method for inverting the high-resolution P- and S-wave velocity ratio time series based on the arrival time of seismic waves according to claim 1, wherein The time series inversion matrix includes: d = Gm wherein, is the average wave velocity ratio measured for all m earthquake pairs, m = [r1, r2, r3,... r n is the n Vp / Vs time series to be determined, G is the sensitivity kernel coefficient matrix of m*n, w ij is the sensitivity of the Vp / Vs of the i-th earthquake pair at the j-th moment, t j and t j+1 respectively represent the j-th moment and the j+1-th moment, t e represents any one of the two earthquake occurrence times of the i-th earthquake pair.

6. The method for inverting the high-resolution P- and S-wave velocity ratio time series based on seismic wave arrival times according to claim 5, characterized in that, Judging whether the wave velocity ratio time series is stable includes: comparing the number of time series unknowns with the model two-norm ||m|| to judge whether the wave velocity ratio time series is stable.

7. The method for inverting the time series of high-resolution P- and S-wave velocity ratios based on seismic wave arrival times according to claim 1, characterized in that, Outputting the inverted wave velocity ratio time series includes: obtaining the wave velocity ratio time series with the number of time instants less than the number of time series unknowns and performing envelope processing, and outputting the inverted wave velocity ratio time series.

Images