Medium and small earthquake focus mechanism solution automatic inversion method

Abstract

The invention discloses an automatic inversion method for medium and small seismic focus mechanism solutions. The method comprises the following steps: preprocessing original seismic data; screening seismic data; calculating the Green function of the whole seismic source area; defining a weight file parameter, and generating a weight file; according to the weight file, the seismic data and the Green function, seismic source mechanism inversion calculation is carried out, and the optimal seismic source mechanism and the optimal seismic source depth are obtained through further searching. According to the invention, an accurate seismic source mechanism can be automatically obtained within a short time after an earthquake occurs; according to the method, waveform fitting condition judgment is introduced, so that the accuracy of an inversion result can be effectively improved.

Description

Technical Field

[0001] The invention belongs to the technical field of focal seismology, and in particular relates to an automatic inversion method for focal mechanism solutions of small and medium earthquakes. Background Art

[0002] The focal mechanism solution is the basic element for understanding the mechanism of earthquake occurrence. Rapidly and accurately determining the focal mechanism solution of the earthquake after a large earthquake occurs can help us understand the possible seismogenic faults and seismogenic mechanisms of the earthquake, so as to provide some necessary help for earthquake relief and earthquake risk assessment in a timely manner. In addition, through the focal mechanism solutions of a large number of small and medium earthquakes, the temporal and spatial variation characteristics of the regional stress field can be obtained, and the gestation and rupture process of strong earthquakes can be understood. The focal mechanism of small and medium aftershocks of strong earthquakes can provide important information for the seismogenic mechanism of the main earthquake.

[0003] Some earthquake-related institutions abroad at present, relying on for example IRIS (Incorporated Research Institutions for Seismology) global seismic network or other regional seismic network, can provide focal mechanism solutions of global 5.5 magnitude or above earthquakes in real time, such as GCMT (The Global Centroid Moment Tensor) institutions of the United States. Also only part of the research team releases the focal mechanism information of larger earthquakes in addition, and there is also a longer time delay (usually the second day after the earthquake). Because the focal mechanism solution of small earthquakes is greatly affected by velocity model and data signal-to-noise ratio, therefore invention is only for 3.5-5.5 magnitude earthquakes, and the goal is to automatically provide reliable focal mechanism solutions in a short time after the earthquake. Summary of the invention

[0004] The purpose of the present invention is to address the problem that seismic institutions lack real-time focal mechanism information and to provide an automatic inversion method for focal mechanism solutions of small and medium-sized earthquakes. This method does not require manual data preprocessing and focal mechanism inversion, shortens calculation time and reduces manual workload, and provides stable and reliable earthquake focal mechanism solutions.

[0005] The above-mentioned object of the present invention is achieved by the following technical solutions:

[0006] An automatic inversion method for focal mechanism solutions of small and medium earthquakes comprises the following steps:

[0007] Step 1: Preprocess the original seismic data to obtain the seismic data in SAC format and the corresponding station information;

[0008] Step 2: Delete the earthquake data whose amplitude ratio of the waveform before and after the P wave is less than the amplitude ratio threshold of the waveform before and after the P wave;

[0009] Step 3: Calculate the Green's function of the entire source area based on the station information and earthquake event information;

[0010] Step 4: Define weight file parameters and generate weight file;

[0011] Step 5: Perform focal mechanism inversion calculation according to the weight file, the seismic data processed in step 2, and the Green's function obtained in step 3 to obtain the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data, and obtain the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data;

[0012] Step 6: Based on the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data of the station, as well as the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data, by searching the focal mechanism and focal depth parameters, the error between the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data and the error between the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data are minimized, and the optimal focal mechanism and focal depth are obtained.

[0013] The amplitude ratio of the waveform before and after the P wave in step 2 as described above is obtained based on the following steps:

[0014] The theoretical P-wave arrival time of the seismic data is calculated based on the seismic data in the SAC format and the theoretical arrival time software. The ratio of the seismic data in the set time window before the theoretical P-wave to the seismic data in the set time window after the theoretical P-wave is calculated to obtain the waveform amplitude ratio before and after the P-wave.

[0015] The calculation of the Green's function of the entire source area in step 3 as described above includes the following steps:

[0016] Step 3.1: Based on the earthquake location carried in the seismic data, the velocity structure model of the source area is obtained using the global earth shell model;

[0017] Step 3.2: Calculate the Green's function based on the generated velocity structure model of the source area and the station information output in step 1.

[0018] As mentioned above, step 5 also includes the step of judging based on the actual observed waveforms of the body wave component and the surface wave component of the seismic data and the fitting results of the theoretical synthetic waveforms:

[0019] If the cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform of the body wave component of the seismic data of the station is less than the cross-correlation coefficient threshold and the amplitude ratio between the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data of the station is greater than the ratio threshold, the body wave component of the seismic data of the station is removed, a new weight file is regenerated and step 5 is re-implemented;

[0020] If the cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform of the surface wave component of the seismic data of the station is less than the cross-correlation coefficient threshold and the amplitude ratio between the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data of the station is greater than the ratio threshold, then remove the surface wave component of the seismic data of the station, regenerate a new weight file and re-implement step 5.

[0021] A computer device comprises a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps of the above-mentioned inversion method when executing the computer program.

[0022] A computer-readable storage medium stores a computer program, which implements the steps of the above-mentioned inversion method when executed by a processor.

[0023] A computer program product comprises a computer program, which implements the steps of the above-mentioned inversion method when executed by a processor.

[0024] Compared with the prior art, the present invention has the following advantages:

[0025] 1. The present invention can automatically obtain accurate focal mechanisms within a short period of time after an earthquake occurs;

[0026] 2. The present invention divides the entire inversion process into blocks, and can replace the corresponding contents of step 1, step 2, step 3 and step 4 in the implementation scheme as needed;

[0027] 3. The present invention introduces waveform fitting situation judgment, which can effectively improve the accuracy of inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 It is a flow chart of the present invention.

[0029] Figure 2 The waveform fitting results of the outputs of some stations for the 4.9-magnitude earthquake in Xingwen, Sichuan on May 4, 2023 are calculated for the present invention.

[0030] Figure 3 The distribution of earthquake stations with a magnitude of 4.9 in Xingwen, Sichuan on May 4, 2023 is calculated for the present invention;

[0031] Figure 4 The diagram is a comparison between the focal mechanism results calculated by the present invention and the focal mechanism solution results calculated manually; wherein (a) is a diagram of the focal mechanism results calculated by the present invention; and (b) is a diagram of the focal mechanism results calculated manually. DETAILED DESCRIPTION

[0032] In order to facilitate those skilled in the art to understand and implement the present invention, the present invention is further described in detail below in conjunction with embodiments. It should be understood that the implementation examples described herein are only used to illustrate and explain the present invention and are not used to limit the present invention.

[0033] Embodiment 1:

[0034] An automatic inversion method for focal mechanism solutions of small and medium earthquakes, such as Figure 1 As shown, the following steps are included:

[0035] Step 1: Preprocess the original seismic data of each station to obtain the seismic data in SAC format and the corresponding station information;

[0036] Step 1 specifically includes the following steps:

[0037] The original seismic data are automatically de-instrumented, de-meaned, de-linearized, de-wrinkled, cut to 300 seconds, rotated the three-component coordinates, resampled to 10 Hz, and renamed, and the station information and the corresponding SAC format seismic data are output.

[0038] Step 2: Remove seismic data with low signal-to-noise ratio based on the amplitude ratio threshold of the waveform before and after the P wave.

[0039] Step 2 specifically includes the following steps:

[0040] The theoretical P wave arrival time of the seismic data is calculated based on the seismic data in SAC format and the theoretical arrival time software (such as the open source Taup software), and then the ratio of the seismic data in the set window (5s) before the theoretical P wave to the seismic data in the set window (5s) after the theoretical P wave is calculated. If the ratio is less than the amplitude ratio threshold of the waveform before and after the P wave, which is 10 in this embodiment, the seismic data is deleted.

[0041] Step 3: Calculate the Green's function of the entire source area based on the station information and earthquake event information.

[0042] Step 3 specifically includes the following steps:

[0043] Step 3.1: According to the earthquake position carried in the initial seismic data, the velocity structure model of the source area is automatically obtained using the global earth crust model (CRUST2.0), or the velocity structure model of the source area prepared in advance is retrieved.

[0044] Step 3.2: Based on the generated velocity structure model of the source area (or the velocity structure model of the source area prepared in advance) and the station information output in step 1, call the open source frequency-wavenumber Green's function (FK) program to calculate the Green's function. The FK program is fast and accurate, and is a commonly used method for calculating the Green's function. In addition, the Green's function of the entire source area can be calculated in advance to form a Green's function library.

[0045] Step 4: Define the weight file parameters and generate the weight file. The weight file is a control file for the input parameters of seismic data in the inversion calculation. The weight file determines which stations' seismic data participate in the inversion, and determines the proportion of the body wave component and the surface wave component of the seismic data recorded by each station in the inversion. In the present invention, the weight of the body wave component and the surface wave component of the seismic data is initially set to 2:1.

[0046] Step 5: Based on the weight file, the seismic data processed in step 2, and the Green's function (or Green's function library) automatically calculated in step 3, the open source general gCAP program (general Cut and Past) is used to perform focal mechanism inversion calculations. The gCAP program divides the seismic data into body wave components and surface wave components, and then calculates the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data, and calculates the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data.

[0047] Step 5 also includes determining whether to perform secondary calculation based on the fitting results of the actual observed waveforms of the body wave component and the surface wave component of the seismic data with the theoretical synthetic waveforms.

[0048] The actual observed waveforms of body wave components and surface wave components are fitted with the theoretical synthetic waveforms:

[0049] If the cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform of the body wave component of the station's seismic data is less than the cross-correlation coefficient threshold of 50 and the amplitude ratio between the theoretical synthetic waveform and the actual observed waveform of the body wave component of the station's seismic data is greater than the ratio threshold of 10, then remove the body wave component of the seismic data of the station, regenerate a new weight file, and re-implement step 5 to perform focal mechanism inversion calculation.

[0050] If the cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform of the surface wave component of the seismic data of the station is less than the cross-correlation coefficient threshold of 50 and the amplitude ratio between the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data of the station is greater than the ratio threshold of 10, then the surface wave component of the seismic data of the station is removed, a new weight file is regenerated, and step 5 is re-implemented to perform the focal mechanism inversion calculation.

[0051] Step 6: Based on the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data of the station, as well as the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data, by searching the focal mechanism and focal depth parameters, the error between the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data and the error between the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data are minimized, and the optimal focal mechanism and focal depth are obtained.

[0052] A person of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing related hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods.

[0053] Embodiment 2:

[0054] Taking the 4.9 magnitude earthquake that occurred in Xingwen, Sichuan on May 4, 2023 as an example, the present invention is applied to the actual earthquake event for inversion calculation, and steps 1-6 of Example 1 are implemented. The waveform fitting results of the actual observed waveforms of the output body wave components and surface wave components of some stations and the theoretical synthetic waveforms are as follows: Figure 2 The stations are distributed as shown in Figure 3 shown. Figure 2 In the graph, Pz and Pr represent body wave components, Sz, Sr and Sh represent surface wave components; the two rows of letters and numbers on the left side of the waveform graph represent the station name and the station epicentral moment, respectively; the two rows of numbers below the waveform graph represent the waveform cross-correlation time shift and waveform cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform, respectively. Figure 3 In the figure, the five-pointed star indicates the location of the 4.9-magnitude Xingwen earthquake in 2023, the triangle indicates the distribution location of the stations, and the letters under the triangle indicate the names of the stations. The mutual correlation coefficients of the waveform fitting results of the actual observed waveforms of the body wave components and surface wave components of most stations and the theoretical synthetic waveforms are greater than 60%, especially for the surface wave component fitting results, the waveform mutual correlation coefficients of most stations are greater than 70%. The variance reduction is 84.3%, indicating that the overall fitting effect is good (the closer the variance reduction is to 100%, the smaller the fitting error). In order to further explore the reliability of the results of the present invention, the manual inversion focal mechanism was used to compare the results, as shown in the figure. Figure 4 As shown. The focal mechanism solution calculated manually is similar to the focal mechanism solution automatically calculated by the present invention, indicating that the inversion result of the present invention is robust and reliable. In addition, the entire focal mechanism automatic inversion process of the present invention takes 46 minutes, of which the calculation of the Green's function in step 3 takes 29 minutes, steps 5-6 take 16 minutes, and the rest takes less than 1 minute in total. Therefore, from the perspective of calculation time, the present invention can give the focal mechanism result in a relatively short time after the earthquake occurs. If the Green's function library is prepared in advance, the calculation time can be greatly reduced.

[0055] The present invention can effectively and automatically calculate the focal mechanism solution of 3.5-5.5 magnitude earthquakes, which is reflected in: (1) data preprocessing, using the general processing flow of seismic data, the preprocessing control parameters are few and the influence is small, so the automatic data preprocessing results in the present invention are basically consistent with the manual processing results. (2) The velocity structure model of the focal region will affect the accuracy of the focal mechanism solution inversion results to a certain extent, especially in areas where the velocity structure changes dramatically, such as the southeastern edge of the Qinghai-Tibet Plateau. The present invention mainly reduces the influence of the velocity structure model of the focal region from the following three aspects: first, the velocity model of the area with drastic velocity changes is built-in according to the published research results, and then the velocity structure model of the focal region of the corresponding area is directly called according to the focal position; second, the waveform fitting situation of the actual observed waveform and the theoretical synthetic waveform is introduced in step 5 to determine the waveform fitting situation, and the seismic data of the stations with poor fitting in the inversion process are removed and refitted; third, the selected focal mechanism inversion method is highly robust and less dependent on the velocity structure model of the focal region.

[0056] Embodiment 3:

[0057] This embodiment provides a computer device, including a memory and a processor. The memory stores a computer program, and the processor implements the steps in the above-mentioned embodiment 1 when executing the computer program.

[0058] Embodiment 4:

[0059] This embodiment provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the steps in the above-mentioned embodiment 1 are implemented.

[0060] Embodiment 5:

[0061] This embodiment provides a computer program product, including a computer program, which implements the steps in the above-mentioned embodiment 1 when executed by a processor.

[0062] The above are only preferred embodiments of the present invention, and the protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, some improvements and modifications without departing from the principle of the present invention should be regarded as the protection scope of the present invention.

Claims

1. A method for automatic inversion of focal mechanism solutions for small and medium earthquakes, characterized in that: The following steps are involved: Step 1: Preprocess the original seismic data to obtain the seismic data in SAC format and the corresponding station information; Step 2: Delete the earthquake data whose amplitude ratio of the waveform before and after the P wave is less than the amplitude ratio threshold of the waveform before and after the P wave; Step 3: Calculate the Green's function of the entire source area based on the station information and earthquake event information; Step 4: Define weight file parameters and generate weight file; Step 5: Perform focal mechanism inversion calculation according to the weight file, the seismic data processed in step 2, and the Green's function obtained in step 3 to obtain the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data, and obtain the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data; Step 6: Based on the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data of the station, as well as the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data, by searching the focal mechanism and focal depth parameters, the error between the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data and the error between the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data are minimized, and the optimal focal mechanism and focal depth are obtained.

2. According to claim 1, a method for automatic inversion of focal mechanism solutions of small and medium earthquakes is characterized by: The amplitude ratio of the P wave before and after in step 2 is obtained based on the following steps: The theoretical P-wave arrival time of the seismic data is calculated based on the seismic data in the SAC format and the theoretical arrival time software. The ratio of the seismic data in the set time window before the theoretical P-wave to the seismic data in the set time window after the theoretical P-wave is calculated to obtain the waveform amplitude ratio before and after the P-wave.

3. The automatic inversion method for focal mechanism solutions of small and medium earthquakes according to claim 1, characterized in that: Calculating the Green's function of the entire source area in step 3 includes the following steps: Step 3.1: Based on the earthquake location carried in the seismic data, the velocity structure model of the source area is obtained using the global earth shell model; Step 3.2: Calculate the Green's function based on the generated velocity structure model of the source area and the station information output in step 1.

4. The method for automatic inversion of focal mechanism solutions of small and medium earthquakes according to claim 1, characterized in that: The step 5 also includes the step of judging based on the actual observed waveforms of the body wave component and the surface wave component of the seismic data and the fitting results of the theoretical synthetic waveforms: If the cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform of the body wave component of the seismic data of the station is less than the cross-correlation coefficient threshold and the amplitude ratio between the theoretical synthetic waveform and the actual observed waveform of the body wave component of the seismic data of the station is greater than the ratio threshold, the body wave component of the seismic data of the station is removed, a new weight file is regenerated and step 5 is re-implemented; If the cross-correlation coefficient between the actual observed waveform and the theoretical synthetic waveform of the surface wave component of the seismic data of the station is less than the cross-correlation coefficient threshold and the amplitude ratio between the theoretical synthetic waveform and the actual observed waveform of the surface wave component of the seismic data of the station is greater than the ratio threshold, then remove the surface wave component of the seismic data of the station, regenerate a new weight file and re-implement step 5.

5. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the inversion method according to any one of claims 1 to 4 are implemented.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the inversion method according to any one of claims 1 to 4 are implemented.

7. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the inversion method according to any one of claims 1 to 4 are implemented.

Images