Full moment tensor seismic source mechanism inversion method, system, device and storage medium

Abstract

The present invention discloses a full moment tensor focal mechanism inversion method, system, device and storage medium. The method includes: Step S1, obtaining control parameters, where the control parameters include: velocity model, focal depth and station location; Step S2, calculating Green's functions with different focal depths by using the fast generalized inverse transmission coefficient algorithm with the control parameters obtained in Step S1; Step S3, substituting each Green's function into the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion to obtain the final focal depth and focal mechanism solution. By using the fast generalized inverse transmission coefficient algorithm to calculate Green's functions, not only the calculation efficiency is improved, but also the storage amount is reduced, and seismic records with any station-source distribution can be calculated. By using the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion, the efficiency of focal mechanism inversion is improved compared with grid search.

Description

Technical Field

[0001] The present invention relates to the field of seismic source mechanism inversion, and particularly to a full moment tensor seismic source mechanism inversion method, system, device and storage medium. Background Art

[0002] Seismic source mechanism inversion is an important means to reveal the activities of seismic faults or fractures. The obtained seismic moment tensor can reflect the type of rupture, estimate the scale of rupture and the magnitude of energy release, etc. For most seismic faults, it can be explained by a simple dislocation model (as shown in Figure 1 ), that is, an earthquake occurs due to relative sliding on the fault plane. There are three control parameters in the simple dislocation model. Among them, the fault strike φ and dip angle δ control the direction of the fault plane, and the slip angle λ controls the relative sliding on both sides of the fault plane. However, when an earthquake is caused by volcanic eruption, hydraulic fracturing, collapse or explosion, etc., the simple dislocation model is not applicable. At this time, two dimensionless quantities β and γ need to be introduced to determine the proportion of non-simple dislocation in the seismic moment tensor. Therefore, accurately obtaining the seismic moment tensor can give relevant parameters such as the type of seismic event. Therefore, seismic source mechanism inversion is widely used not only in natural earthquake observation, but also in the fields of hydraulic fracturing for oil and gas exploitation, nuclear explosion monitoring, etc. For example, if it is a hydraulic fracturing earthquake, the scale of fracturing can be estimated and the fracturing effect can be evaluated; if it is an explosion, the type of explosion and the equivalent can be judged and estimated.

[0003] There are two existing full waveform seismic source mechanism inversion algorithms. One is the time-domain seismic source mechanism inversion algorithm (TDMT_INVC, Time-Domain Moment Tensor INVerse Code) proposed by Dreger and Helmberger in 1993, and the algorithm flow is as shown in Figure 2 . The other is the GCAP method (Cut And Paste) improved by Zhu and Helmberger in 1996 (the algorithm flow is as shown in Figure 3 ). By comparing the two, it is found that the calculation idea of the Green's function is the same, which is to calculate the three-component displacement or velocity records caused by three simple dislocation sources and one explosion source at the station respectively and save these 12 quantities. In the subsequent inversion process, TDMT_INVC adopts the idea of linearizing the nonlinear problem and then performing iterative inversion, while the GCAP method uses grid search to directly search for 5 parameters such as strike, dip angle, and slip angle.

[0004] When existing algorithms for focal mechanism inversion using waveforms generate Green's functions, they generally calculate the displacements caused by four specific faults proposed by Minson and Dreger (2008) as Green's functions. These four specific faults are: vertical strike-slip fault, vertical dip-slip fault, dip-slip fault with a 45° dip angle, and explosion source. However, this calculation method requires calculating the displacements caused by four specific faults separately, and this calculation process may require two or four times for some algorithms. For example, FK (Wand & Herrmann, 1980; Zhu & Rivera, 2002) needs to calculate twice and store the records of three components of four seismic sources.

[0005] The existing methods for calculating Green's functions have at least the following disadvantages:

[0006] (1) It is necessary to perform two or four forward calculations to obtain the Green's functions required for full moment tensor inversion, and the calculation efficiency is low.

[0007] (2) It is necessary to store the records of 12 components, and the storage capacity is large.

[0008] (3) The calculation speed is slow, especially when the focal depth and the station burial depth are close.

[0009] (4) It cannot calculate the cases where the focal point and the station are at the same depth or the station has a certain burial depth. Such situations mainly occur in the observation of hydraulic fracturing earthquakes, the observation of earthquakes caused by collapses or explosions.

[0010] The existing grid search algorithm used in inversion has the disadvantage that the inversion speed and inversion accuracy are controlled by the grid size. The smaller the grid, the slower the inversion speed.

[0011] In view of this, the present invention is specifically proposed. Summary of the Invention

[0012] The purpose of the present invention is to provide a full moment tensor focal mechanism inversion method, system, device and storage medium, which can improve the calculation efficiency of Green's functions and reduce storage, as well as improve the efficiency of focal mechanism inversion, thereby solving the above technical problems existing in the prior art.

[0013] The purpose of the present invention is achieved through the following technical solutions:

[0014] A full moment tensor focal mechanism inversion method includes:

[0015] Step S1, obtaining control parameters from each seismic data center or from the obtained original waveform data. The control parameters include: velocity model, focal position, and station position;

[0016] Step S2: Calculate Green's functions for different source depths by using the control parameters obtained in Step S1 through the fast generalized inverse transmission coefficient algorithm;

[0017] Step S3: Substitute the calculated Green's functions for each source depth into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion, obtain the source mechanism inversion results and error functions for different source depths, and select the source depth and source mechanism solution with the minimum error function as the final source depth and source mechanism solution of the source, thus completing the inversion of the source mechanism.

[0018] A full moment tensor source mechanism inversion system, comprising:

[0019] A control parameter acquisition unit, a Green's function calculation unit, and a source mechanism inversion unit; wherein,

[0020] The control parameter acquisition unit is configured to acquire control parameters, and the control parameters include: a velocity model, a source location, and a station location;

[0021] The Green's function calculation unit is communicatively connected to the control parameter acquisition unit, and is configured to calculate Green's functions for different source depths by using the control parameters acquired by the control parameter acquisition unit through the fast generalized inverse transmission coefficient algorithm;

[0022] The source mechanism inversion unit is communicatively connected to the Green's function calculation unit, and is configured to substitute the Green's functions for each source depth calculated by the Green's function calculation unit into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion, obtain the source mechanism inversion results and error functions for different source depths, and select the source depth and source mechanism solution with the minimum error function as the final source depth and source mechanism solution of the source.

[0023] A processing device, comprising:

[0024] At least one memory for storing one or more programs;

[0025] At least one processor capable of executing the one or more programs stored in the memory, and when the one or more programs are executed by the processor, enabling the processor to implement the method of the present invention.

[0026] A readable storage medium storing a computer program, which when executed by a processor can implement the method of the present invention.

[0027] Compared with the prior art, the full moment tensor source mechanism inversion method, system, device, and storage medium provided by the present invention have the following beneficial effects:

[0028] By using the fast generalized reflection-transmission coefficient algorithm to calculate the Green's function, not only the calculation efficiency is improved, but also the storage requirement is reduced, and seismic records for any station-source distribution can be calculated. By using the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion, the efficiency of focal mechanism inversion is improved compared to grid search, and the source and stations can be located at any position, with a wider scope of application. Description of the Drawings

[0029] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the description of the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0030] Figure 1 Schematic diagram of a simple dislocation model representing most earthquakes.

[0031] Figure 2 Flowchart of the existing time-domain source parameter inversion method.

[0032] Figure 3 Flowchart of the existing GCAP inversion method.

[0033] Figure 4 Flowchart of the full moment tensor focal mechanism inversion method provided by the embodiments of the present invention.

[0034] Figure 5 Specific processing flowchart of the full moment tensor focal mechanism inversion method provided by the embodiments of the present invention.

[0035] Figure 6 Flowchart of the fast generalized reflection-transmission coefficient algorithm of the full moment tensor focal mechanism inversion method provided by the embodiments of the present invention for calculating Green's functions at different source depths.

[0036] Figure 7 Flowchart of the neighborhood algorithm with the moment tensor distance as the stopping criterion of the full moment tensor focal mechanism inversion method provided by the embodiments of the present invention.

[0037] Figure 8 Schematic diagram of the comparison of calculation results between the fast GRTM algorithm (i.e., Green's function calculation method) provided by the embodiments of the present invention and the existing GRTM algorithm; among them, (a) is the waveform comparison diagram calculated by the original GRTM (solid line) and the fast GRTM of the present invention (dashed line) when both the station and the source are located on the surface; (b) is the comparison of the time required for the original GRTM and the fast GRTM of the present invention to calculate the displacements caused by sources with different source depths at the same epicentral distance for the same station (GZ.BJT), where the triangles represent the original GRTM and the circles represent the FGRTM.

[0038] Figure 9 Schematic diagram for comparing the calculation results of the fast GRTM algorithm provided by the embodiments of the present invention with the existing FK algorithm; wherein, (a) is a waveform comparison diagram of the FK and FGRTM calculations for the same station (GZ.BJT) with the same epicentral distance but different source depths (the source depth range is from 0.1 km to 9 km), where the solid black line is FK and the dashed black line is FGRTM; (b) is a comparison of the calculation times of the two algorithms, where the triangles represent FK and the circles represent FGRTM.

[0039] Figure 10 Block diagram of the structure of the full moment tensor source mechanism inversion system provided by the embodiments of the present invention. Detailed implementation manners

[0040] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the specific content of the present invention; obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments, which does not constitute a limitation to the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0041] First, the following explanations are made for the terms that may be used in this article:

[0042] The term "and / or" means that either one or both of the two can be realized. For example, X and / or Y means that it includes both the case of "X" or "Y" and the three cases of "X and Y".

[0043] The description of terms such as "comprising", "including", "containing", "having" or other similar semantics should be interpreted as non-exclusive inclusion. For example: including a certain technical feature element (such as raw materials, components, ingredients, carriers, dosage forms, materials, dimensions, parts, components, mechanisms, devices, steps, processes, methods, reaction conditions, processing conditions, parameters, algorithms, signals, data, products or articles, etc.) should be interpreted as including not only the clearly listed certain technical feature element, but also other well-known technical feature elements in the art that are not clearly listed.

[0044] The term "consisting of" means excluding any technical feature element that is not clearly listed. If this term is used in a claim, this term will make the claim a closed type, so that it does not include technical feature elements other than the clearly listed ones, except for related conventional impurities. If this term only appears in a certain clause of a claim, then it only limits the elements clearly listed in that clause, and the elements recorded in other clauses are not excluded from the overall claim.

[0045] Unless otherwise clearly specified or defined, the terms "installation", "connection", "attachment", "fixation", etc. shall be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium, and it can be the communication inside two components. For those of ordinary skill in the art, the specific meanings of the above terms in this text can be understood according to specific circumstances.

[0046] The orientation or positional relationship indicated by the terms "center", "longitudinal", "transverse", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", etc. is based on the orientation or positional relationship shown in the drawings. It is only for the convenience of description and simplification of description, rather than explicitly or implicitly indicating that the device or component referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation to this text.

[0047] The full moment tensor source mechanism inversion method, system, device and storage medium provided by the present invention will be described in detail below. The content not described in detail in the embodiments of the present invention belongs to the prior art well-known to those of ordinary skill in the art. In the embodiments of the present invention, those not specified in specific conditions are carried out according to the conventional conditions in the art or the conditions recommended by the manufacturer. For the reagents or instruments not specified in the production manufacturer in the embodiments of the present invention, they are all conventional products that can be obtained through commercial purchase.

[0048] As Figure 4 shown, the embodiments of the present invention provide a full moment tensor source mechanism inversion method, including:

[0049] Step S1, obtaining control parameters from each seismic data center or from the obtained original waveform data, where the control parameters include: velocity model, source location, and station location;

[0050] Step S2, calculating Green's functions at different source depths by using the control parameters obtained in step S1 through the fast generalized anti-transmission coefficient algorithm;

[0051] Step S3, substituting the obtained Green's functions at each source depth into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion, obtaining the source mechanism inversion results and error functions at different source depths, and selecting the source depth and source mechanism solution with the minimum error function as the final source depth and source mechanism solution of the source.

[0052] See Figure 6, in step S2 of the above method, the Green's functions with different source depths are calculated by the fast generalized inverse transmission coefficient algorithm using the obtained control parameters in the following manner, including:

[0053] Step S21, for each pair of source stations, calculate the generalized inverse transmission coefficients with different frequencies and wavenumbers between the source and the station;

[0054] Step S22, calculate the Bessel function values of Bessel functions of different orders at different frequencies and wavenumbers;

[0055] Step S23, combining the generalized inverse transmission coefficients and Bessel function values calculated in steps S21 and S22, calculate the integrand values of the Green's function I i at different frequencies and wavenumbers, i = 1…10;

[0056] Step S24, use the valley-peak averaging method to perform wavenumber integration on the Green's function I i values obtained in step S23, and then perform inverse Fourier transform to obtain its value in the time domain, that is, the Green's functions with different source depths are obtained.

[0057] See Figure 7 , in step S3 of the above method, substitute the obtained Green's functions with each source depth into the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion in the following manner, including:

[0058] Step S31, initialize, generate n s initial models uniformly randomly distributed in the model space and calculate the error function of each initial model. The error function e is: where, is the observation record of a certain source station; is the synthetic record calculated according to the formula for the displacement caused by the earthquake at this source station; cc is and the correlation coefficient of; the model space is composed of 5 parameters such as fault strike, dip angle, and slip angle, and each initial model is composed of 5 random parameters;

[0059] Step S32, start iterative solution, let loop = 1;

[0060] Step S33, select the first n s models with the smallest error function from the loop×n r generated in the previous loop steps;

[0061] Step S34, regenerate n r in the Voronoi body composed of the model parameters of n sA new model is generated and the error function of each new model is calculated. Let loop = loop + 1;

[0062] Step S35: Determine whether the number of iterations has reached the preset maximum number of iterations or the average moment tensor distance is less than the preset value. If not, return to step S33 to continue the iteration; if so, execute step S36;

[0063] Step S36: Output the model parameters with the minimum error function, which are the source depth and source mechanism solution with the minimum error function, as the final source depth and source mechanism solution of the source, that is, the inversion of the source mechanism is completed.

[0064] See Figure 10 , the embodiment of the present invention also provides a full moment tensor source mechanism inversion system, including:

[0065] A control parameter acquisition unit, a Green's function calculation unit, and a source mechanism inversion unit; where

[0066] The control parameter acquisition unit is used to obtain control parameters from each seismic data center or from the obtained original waveform data (i.e., seismic original waveform data). The control parameters include: velocity model, source location (latitude and longitude), and station location;

[0067] The Green's function calculation unit is communicatively connected to the control parameter acquisition unit and is used to calculate Green's functions of different source depths through the fast generalized anti-transmission coefficient algorithm using the control parameters obtained by the control parameter acquisition unit;

[0068] The source mechanism inversion unit is communicatively connected to the Green's function calculation unit and is used to substitute the Green's functions of each source depth calculated by the Green's function calculation unit into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion, obtain the source mechanism inversion results and error functions of different source depths, and select the source depth and source mechanism solution with the minimum error function as the final source depth and source mechanism solution of the source.

[0069] See Figure 6 , in the Green's function calculation unit of the above system, the Green's functions of different source depths are calculated through the fast generalized anti-transmission coefficient algorithm using the obtained control parameters in the following manner, including:

[0070] Step S21: For each source-station pair, calculate the generalized anti-transmission coefficients of different frequencies and wavenumbers between the source and the source station;

[0071] Step S22: Calculate the Bessel function values of different orders of the Bessel function at different frequencies and wavenumbers;

[0072] Step S23: Calculate the integrand values of the Green's function I at different frequencies and wavenumbers by combining the generalized anti-transmission coefficients and Bessel function values calculated in Steps S21 and S22, where i = 1…10; i

[0073] Step S24: Perform wavenumber integration on the Green's function I values obtained in Step S23 using the valley-peak averaging method, and then perform an inverse Fourier transform to obtain its value in the time domain, that is, obtain the Green's function at different source depths. i

[0074] See Figure 7 In the source mechanism inversion unit of the above system, substitute the Green's function of each obtained source depth into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion in the following manner, including:

[0075] Step S31: Initialize, generate n initial models uniformly randomly distributed in the model space and calculate the error function of each initial model. The error function e is: s where is the observation record of a certain source station; is the synthetic record calculated according to the formula for the displacement caused by the earthquake at this source station; cc is and the correlation coefficient of;

[0076] Step S32: Start iterative solution, let loop = 1;

[0077] Step S33: Select the first n models with the smallest error function from the loop×n generated in the previous loop steps s ; r

[0078] Step S34: Regenerate n new models in the Voronoi body composed of the n model parameters and calculate the error function of each new model, let loop = loop + 1; r s

[0079] Step S35: Determine whether the number of iterations reaches the preset maximum number of iterations or the average moment tensor distance is less than the preset value. If not, return to Step S33 to continue the iteration; if so, execute Step S36;

[0080] Step S36: Output the model parameters with the smallest error function, that is, the source depth and source mechanism solution with the smallest error function, as the final source depth and source mechanism solution of this source, that is, complete the inversion of the source mechanism.

[0081] The embodiment of the present invention also provides a processing device, including:​​​​​​

[0082] At least one memory for storing one or more programs;

[0083] At least one processor capable of executing the one or more programs stored in the memory, and when the one or more programs are executed by the processor, enabling the processor to implement the above method.

[0084] An embodiment of the present invention further provides a readable storage medium storing a computer program, which can implement the above method when executed by a processor.

[0085] In summary, in the method and system of the embodiment of the present invention, by using the fast generalized back-transmission coefficient algorithm to calculate the Green's function, since the calculation steps of correcting the back-transmission coefficient are saved, not only the calculation efficiency is improved, but also the storage amount is reduced. Moreover, the peak-valley averaging method is introduced to calculate the wave velocity integral, and the seismic records of any station-source distribution can be calculated, even if the station and the source are at the same depth (for this case, other FK algorithms cannot calculate). By using the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion, the efficiency of focal mechanism inversion is improved compared with grid search.

[0086] In order to more clearly show the technical solutions provided by the present invention and the technical effects produced, the full moment tensor focal mechanism inversion method and system provided by the embodiments of the present invention will be described in detail below with specific embodiments.

[0087] Embodiment 1

[0088] As Figure 4 、 5 shown, this embodiment provides a full moment tensor focal mechanism inversion method, including:

[0089] Step S1: Obtain control parameters from each seismic data center (or from the obtained original waveform data (i.e., seismic original waveform data)), and the control parameters include: velocity model, source location, and source-station location;

[0090] Step S2: Calculate the Green's function at different source depths by using the control parameters obtained in Step S1 through the fast generalized back-transmission coefficient algorithm;

[0091] Step S3: Substitute the Green's function at each source depth obtained into the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion to obtain the focal mechanism inversion results and error functions at different source depths;

[0092] Step S4: Select the source depth and focal mechanism solution with the minimum error function as the final source depth and focal mechanism solution of the source, that is, complete the focal mechanism inversion.

[0093] See Figure 6 In the above step S2, the fast generalized back-transmission coefficient algorithm calculates the Green's function I i The process is as follows Figure 5 As shown, its calculation principle is based on the following formula (1). According to the principle of the generalized back-transmission algorithm (Chen, 1999, Eq. 77a-c and 78a-c), for an earthquake, the displacement caused at any station can be rewritten as:

[0094]

[0095] where θ is the station azimuth; the Green's function I i , i = 1…10 are wavenumber integral functions related to the Bessel function and the generalized back-transmission coefficient, and the formula is as follows.

[0096]

[0097] where J i and J i ′, i = 0, 1, 2 are the Bessel function and its derivative; ω is the frequency of the calculated waveform; k is the wavenumber; r is the distance from the earthquake source center to the station; z is the earthquake source depth; and is related to the generalized back-transmission coefficient of the j-th layer (for the specific formula, see the content of Chen, 1999, Eq. 74a-c); M pq , p, q = x, y, z are the normalized seismic moment tensors, which are functions of the fault strike, dip, slip angle, and β and γ.

[0098] Therefore, based on the above formula (1), the earthquake source mechanism inversion can be carried out, and I i can be regarded as an equivalent Green's function, which is only related to the epicentral distance and earthquake source depth of the station and has nothing to do with the station azimuth. In addition, in the original GRTM algorithm, before calculating the generalized back-transmission coefficient, a modified back-transmission coefficient needs to be calculated. And the fast generalized back-transmission coefficient algorithm used in the present invention removes the calculation part of the modified back-transmission coefficient and directly calculates the generalized back-transmission coefficient. In this way, the calculation efficiency of the fast generalized back-transmission coefficient algorithm of the present invention is 40% higher than that of the original GRTM (see Figure 8 b); in addition, due to the introduction of the peak-valley averaging method to calculate the wave velocity integral, the fast generalized back-transmission coefficient algorithm can calculate the seismic records of any station-earthquake source distribution, and can also calculate even if the station and the earthquake source are at the same depth (for this case, other FK algorithms cannot calculate) ( Figure 8 a).

[0099] Because the source mechanism, i.e., the seismic moment tensor, is a function of the fault strike 0 ≤ φ < 360, dip 0 ≤ δ ≤ 90, slip angle -180 < λ ≤ 180, and the colatitude 0 ≤ β ≤ 180 and longitude -30 ≤ γ ≤ 30 on the source sphere, when performing source mechanism inversion, the method of the present invention does not directly invert m pq Instead, it inverses the above five angular values. An error function e will be set during the inversion process: wherein, is the observation record of the source station; is the synthetic record calculated according to formula (1); cc is and the correlation coefficient of.

[0100] See Figure 7 , in the above step S3, the Green's function of each source depth obtained is substituted into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion. The specific process is as follows:

[0101] Step S31, initialization, generate n s initial models uniformly randomly distributed in the model space and calculate the error function of each model;

[0102] Step S32, start iterative solution, let loop = 1;

[0103] Step S33, select the n s models with the smallest error function from the loop × n r generated in the previous loop steps;

[0104] Step S34, regenerate n r new models in the Voronoi body composed of the n s model parameters and calculate the corresponding error function, let loop = loop + 1;

[0105] Step S35, determine whether the number of iterations reaches the preset maximum number of iterations or the average moment tensor distance is less than the preset value. If not, return to step S33 to continue the iteration; if so, execute step 36;

[0106] Step S36, output the model parameters with the smallest error function, which are the source depth and source mechanism solution with the smallest error function, as the final source depth and source mechanism solution of the source, that is, complete the source mechanism inversion.

[0107] In the neighborhood algorithm of the present invention with the moment tensor distance as the stopping criterion, the results of the previous iterative steps generate new models near the optimal solution. Compared with the model space with non-optimal solutions, this algorithm is more inclined to the model space near the optimal solution. This neighborhood algorithm has two control parameters: the number of models to be generated in each iteration ns and the number n of Voronoi bodies to be resampled r . In this neighborhood algorithm, assume n s and n r are equal. At this time, the neighborhood algorithm is closer to a random sampling algorithm rather than an optimization algorithm. In addition, the concept of two moment tensor distances is introduced in this neighborhood algorithm. The distance between two tensors X = (x ij ) and Y = (y ij ) is defined as where, ‖X‖ = X·X. The value range of the distance between two tensors is [0, 180].

[0108] In summary, the method and system of the embodiment of the present invention calculate the Green's function by using the fast generalized anti-transmission coefficient algorithm. Since the calculation step of correcting the anti-transmission coefficient is omitted, not only the calculation efficiency is improved, but also the storage amount is reduced. The peak-valley averaging method is introduced to calculate the wave velocity integral, and the Green's function I of any station-source distribution can be calculated i . Moreover, under the same calculation accuracy, the calculation efficiency can be greatly improved ( Figure 9 (a), Figure 9 (b)). In addition, for each station, the traditional FK algorithm needs to store 12 components as the Green's function, while the fast generalized anti-transmission coefficient algorithm of the present invention only needs to store 10 I i , reducing the storage amount by about 17%; the neighborhood algorithm with the moment tensor distance as the stopping criterion in step 3 of the present invention greatly reduces the number of error function calculations compared with the traditional grid search algorithm. And the efficiency of inversion is proportional to the number of error function calculations. The more times the error function is calculated, the slower the inversion speed. If the grid size of the grid search is 10°, then the accuracy of the grid search is 10°. Therefore, to complete one inversion, the grid search needs to perform 1,723,680 error function calculations. For this neighborhood algorithm with the moment tensor distance as the stopping criterion, after testing, when n s is 1000, both the inversion efficiency and accuracy can be taken into account. Assuming the maximum number of iterations is 50, then the number of error function calculations required by the neighborhood algorithm is 51,000. Therefore, the neighborhood algorithm is at least 33 times faster than the grid search, greatly improving the inversion efficiency.

[0109] Those of ordinary skill in the art can understand that all or part of the processes in the above-described embodiment methods can be completed by instructing relevant hardware through a program. The program can be stored in a computer-readable storage medium. When the program is executed, it can include the processes of the above-described method embodiments. Among them, the storage medium can be a magnetic disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM), etc.

[0110] As described above, the above are only the preferred specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed by the present invention should be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims. The information disclosed in the background art part of this article is only intended to deepen the understanding of the overall background art of the present invention, and should not be regarded as an admission or any form of implication that this information constitutes the prior art already known to those skilled in the art.

Claims

1. A full moment tensor seismic source mechanism inversion method, characterized in that, Including: Step S1: Obtain control parameters from each seismic data center or from the acquired original waveform data. The control parameters include: velocity model, source location, and station location. Step S2: Calculate Green's functions for different source depths using the control parameters obtained in Step S1 through the fast generalized reflection-transmission coefficient algorithm. Specifically: Step S21: For each source-station pair, calculate the generalized reflection-transmission coefficients at different frequencies and wavenumbers between the source and the station. Step S22: Calculate the Bessel function values of different orders of the Bessel function at different frequencies and wavenumbers. Step S23: Calculate the Green's function I at different frequencies and wavenumbers by combining the generalized anti-transmission coefficient and the Bessel function value calculated in steps S21 and S22 i Calculate the value of the integrand for i = 1…10; Step S24: Use the valley-peak averaging method to perform wavenumber integration on the Green's function values obtained in Step S23, and then perform inverse Fourier transform to obtain their values in the time domain. Step S3: Substitute the calculated Green's functions for each source depth into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion, obtain the source mechanism inversion results and error functions for different source depths, and select the source depth and source mechanism solution with the minimum error function as the final source depth and source mechanism solution for this source, thus completing the inversion of the source mechanism.

2. The full moment tensor seismic source mechanism inversion method according to claim 1, wherein In Step S3, substitute the obtained Green's functions for each source depth into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion in the following manner, including: Step S31, initialization, generate n initial models uniformly randomly distributed in the model space and calculate the error function of each initial model. The error function e is as follows: s where is the observation record of a certain seismic source station; is the synthetic record calculated according to the formula of the displacement caused by the earthquake at this seismic source station; cc is the correlation coefficient between and . Step S32: Start iterative solution, and set loop = 1. Step S33, select n s models with the smallest error function from the loop×n r models generated in the previous loop step; Step S34, regenerate n r new models in the Voronoi body composed of n s model parameters, calculate the error function of each new model, and set loop = loop + 1; Step S35: Determine whether the number of iterations has reached the preset maximum number of iterations or the average moment tensor distance is less than the preset value. If not, return to Step S33 to continue the iteration; if so, execute Step S36. Step S36: Output the model parameters with the minimum error function, which are the source depth and source mechanism solution with the minimum error function, as the final source depth and source mechanism solution for this source, thus completing the inversion of the source mechanism.

3. A full moment tensor seismic source mechanism inversion system, characterized in that, Including: A control parameter acquisition unit, a Green's function calculation unit, and a source mechanism inversion unit; where The control parameter acquisition unit is used to acquire control parameters, and the control parameters include: velocity model, source location, and station location. The Green's function calculation unit is communicatively connected to the control parameter acquisition unit and is used to calculate Green's functions for different source depths using the control parameters acquired by the control parameter acquisition unit through the fast generalized reflection-transmission coefficient algorithm. Specifically: Step S21: For each source-station pair, calculate the generalized reflection-transmission coefficients at different frequencies and wavenumbers between the source and the station. Step S22: Calculate the Bessel function values of different orders of the Bessel function at different frequencies and wavenumbers. Step S23: Combine the generalized anti-transmission coefficients and Bessel function values calculated in Steps S21 and S22 to calculate the integrand values of the Green's function I at different frequencies and wavenumbers, where i = 1...10; i ​ Step S24: Use the valley-peak averaging method to perform wavenumber integration on the Green's function I obtained in step S23, and then perform an inverse Fourier transform to obtain its value in the time domain, that is, obtain the Green's function for different source depths; i value, and then perform an inverse Fourier transform to obtain its value in the time domain, that is, obtain the Green's function for different source depths; The source mechanism inversion unit is communicatively connected to the Green's function calculation unit and is used to substitute the Green's functions for each source depth calculated by the Green's function calculation unit into the neighborhood algorithm with the moment tensor distance as the stopping criterion for source mechanism inversion, obtain the source mechanism inversion results and error functions for different source depths, and select the source depth and source mechanism solution with the minimum error function as the final source depth and source mechanism solution for this source.

4. The full moment tensor seismic source mechanism inversion system according to claim 3, characterized in that In the focal mechanism inversion unit, the Green's function of each focal depth obtained is substituted into the neighborhood algorithm with the moment tensor distance as the stopping criterion for focal mechanism inversion in the following manner, including: Step S31, initialization: Generate n initial models uniformly and randomly distributed in the model space and calculate the error function for each initial model. The error function e is as follows: s where is the observation record of a certain seismic source station; is the synthetic record calculated according to the formula for the displacement caused by the earthquake at this seismic source station; cc is the correlation coefficient between and . Step S32, start iterative solution, and set loop = 1; Step S33, select the top n models with the smallest error function from the loop×n generated in the previous loop step s ; r ​ Step S34, regenerate n r new models in the Voronoi body composed of n s model parameters, calculate the error function of each new model, and let loop = loop + 1; Step S35, determine whether the number of iterations reaches the preset maximum number of iterations or the average moment tensor distance is less than the preset value. If not, return to Step S33 to continue the iteration; if so, execute Step S36; Step S36, output the model parameters with the minimum error function, that is, the focal depth and focal mechanism solution with the minimum error function, as the final focal depth and focal mechanism solution of the earthquake source, namely, complete the inversion of the focal mechanism.

5. A processing device, characterized in that, Including: At least one memory for storing one or more programs; At least one processor capable of executing the one or more programs stored in the memory. When the one or more programs are executed by the processor, the processor can implement the method according to any one of claims 1-2.

6. A readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it can implement the method according to any one of claims 1-2.

Images