Method, computer device and readable storage medium for composite seismic source parameter inversion

Abstract

The present invention provides a method for inverting composite seismic source parameters, a computer device, and a readable storage medium. The method includes the steps of: obtaining the original seismic wave band data of the target area, picking the first arrivals from the original seismic wave band data to obtain the arrival time of the first arrivals; obtaining the initial seismic wave velocity model of the area according to the well logging data or past geological data of the target area; using the first arrival information and / or seismic wave waveform of the seismic waves in the target area to achieve the initial positioning of the seismic source; calculating the dip polarization angle and azimuth polarization angle of the three components of the original seismic wave band data according to the initial positioning result of the seismic source and the coordinates of the observation points, and performing vector synthesis on the original seismic waves to obtain new seismic wave band data; performing joint inversion according to the new seismic wave band data and the initial velocity model. The computer device includes: at least one processor for executing the instructions of the method. The present invention has the advantages of being able to achieve multi-mode vibration positioning, accurate confirmation of the vibration source position, and high accuracy.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data monitoring, and specifically, to a method for inverting composite seismic source parameters, a computer device, and a readable storage medium. Background Art

[0002] Earthquakes can be divided into active earthquakes and passive earthquakes, mainly including three-dimensional seismic exploration active source positioning, surface vibration source positioning, small seismic activities caused by natural earthquakes or small-scale natural faults / fractures, vibrations caused by geological disasters, etc. These earthquakes or vibrations are collectively referred to as vibration sources, and the seismic data of these vibration sources are collectively referred to as seismic data. With the rapid development of the national economy, the consumption of oil and gas resources continues to grow. Seismic exploration and engineering construction mainly based on oil and gas resources will all cause a large amount of vibrations to generate seismic wave signals. Since the vibration sources are unknown and have a greater impact on engineering construction and the environment, therefore, how to quickly and accurately locate the vibration sources is of great significance for industrial production, engineering operations, security, etc.

[0003] The present invention uses the seismic wave signals observed and recorded by seismic equipment, and utilizes information such as seismic wave amplitude, polarity, and first arrival to realize the automatic detection and positioning of vibration sources, thereby providing important technical support for geological exploration resource exploration, engineering construction risk assessment, geological disaster early warning, etc., and solving technical problems such as difficult confirmation of vibration source positions, uncertain determination of source types, and low source positioning accuracy in active and passive seismic monitoring. Summary of the Invention

[0004] The purpose of the present invention is to solve at least one of the above-mentioned deficiencies existing in the prior art. For example, one of the purposes of the present invention is to provide a method for inverting composite seismic source parameters based on seismic wave dynamics and kinematic parameters. Another purpose of the present invention is to provide a computer program for executing the method for inverting composite seismic source parameters based on seismic wave dynamics and kinematic parameters.

[0005] To achieve the above purpose, on the one hand, the present invention provides a method for inverting composite seismic source parameters, and the method includes the steps of:

[0006] Obtain the original seismic wave band data of the target area, pick the first arrival of the original seismic wave band data, and obtain the arrival first arrival time;

[0007] Obtain the initial seismic wave velocity model of the area according to the well logging data or past geological data of the target area;

[0008] Utilize the first arrival information of the seismic waves in the target area and / or the seismic wave waveform to achieve the initial positioning of the seismic source;

[0009] Calculate the dip polarization angle and azimuth polarization angle of the three components of the original seismic wave band data based on the initial seismic source location result and the observation point coordinates, and use the dip polarization angle and azimuth polarization angle to perform vector synthesis on the original seismic wave to obtain new seismic wave band data;

[0010] Perform joint inversion according to the new seismic wave band data and the initial velocity model through Equation 1,

[0011]

[0012] where q is an exponent, dimensionless; is the P-wave end time of the i-th observation point, with the unit of s, ms or sampling points; is the SV-wave end time of the i-th observation point, with the unit of s, ms or sampling points; is the SH-wave end time of the i-th observation point, with the unit of s, ms or sampling points; represents the possible seismic source χ S and the relative propagation time of the SV wave to the observation point χ R，i ; represents the possible seismic source χ S and the relative propagation time of the SH wave to the observation point χ R，i ; represents the possible seismic source χ S and the relative propagation time of the P wave to the observation point χ R，i ; a is the weighting coefficient for location inversion based on scalar / vector superposition of P-wave seismic band data, dimensionless; b is the weighting coefficient for location inversion based on scalar / vector superposition of SV-wave seismic band data, dimensionless; c is the weighting coefficient for location inversion based on scalar / vector superposition of SH-wave seismic band data, dimensionless; M is the seismic source mechanism, which is a 3×3 matrix; V(x, y, z) is the new velocity model, with the unit of m / s or km / s; represents the P-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V P (x, y, z) according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ; represents the SV-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V SV (x, y, z) according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ; represents the SH-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V SH (x, y, z) according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ;

[0013] In an exemplary embodiment of one aspect of the present invention, the original seismic wave band data may include P-wave band data, S-wave band data, or both P-wave band data and S-wave band data, where the S-wave includes SV-wave and / or SH-wave.

[0014] In an exemplary embodiment of one aspect of the present invention, the P-wave band data can be expressed as Equation 2,

[0015] S P，k = [S P，1，k (χ R，1 , t1, t0),..., S P，i，k (χ R，i , t i , t0),..., S P，N，k (χ R，N , t N , t0)] Equation 2

[0016] Wherein, S p,k is the P-wave of the seismic wave recorded by the k-th component, with the unit of m, cm, mm or μm; S P，1，k (χ R，1 , t1, t0) is the seismic P-wave waveform of the k-th component obtained when the observation point is located at χ R，， and the starting time of the band is t ， , with the unit of m, cm, mm or μm; S P，i，k (χ R，i , t i , t0) is the seismic P-wave waveform of the k-th component obtained when the observation point is located at χ R，i and the starting time of the band is t i , with the unit of m, cm, mm or μm; S P，N，k (χ R，N , t N , t0) is the seismic P-wave waveform of the k-th component obtained when the observation point is located at χ R，N and the starting time of the band is t N , with the unit of m, cm, mm or μm; t0 is the earthquake origin time, with the unit of s, ms or time sample points;

[0017] The S-wave band data can be expressed as Equation 3,

[0018] S S，k = [S S，1，k (χ R，1 , t1, t0),..., S S，i，k (χ R，i , t i , t0),..., S S，N，k (χ R，N , t N, t0)] Equation 3

[0019] Wherein, S S,k is the S-wave of the seismic wave recorded for the k-th component, with the unit of m, cm, mm or μm; S S，1，k (χ R，1 , t1, t0) is the waveform of the seismic S-wave of the k-th component obtained when the observation point is located at χ R，1 and the starting time of the wave band is t1, with the unit of m, cm, mm or μm; S S，i，k (χ R，i , t i , t0) is the waveform of the seismic S-wave of the k-th component obtained when the observation point is located at χ R，i and the starting time of the wave band is t i obtained, with the unit of m, cm, mm or μm; S S，N，k (χ R，N , t N , t0) is the waveform of the seismic S-wave of the k-th component obtained when the observation point is located at χ R，N and the starting time of the wave band is t N obtained, with the unit of m, cm, mm or μm; t0 is the earthquake origin time, with the unit of s, ms or the number of time samples.

[0020] In an exemplary embodiment of one aspect of the present invention, the first arrival of the P-wave band data is picked up, and its arrival first arrival time can be expressed as Equation 4,

[0021] T P =(t P，1 ,..., t P，i ,..., t P，N , t0) Equation 4

[0022] Wherein, T P is the first arrival of the P-wave picked up from N trace gathers, with the unit of s or ms; t P，1 is the first arrival of the P-wave of the seismic wave data recorded at the first observation point, with the unit of s, ms or the number of time samples; t P，i is the first arrival of the P-wave of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or the number of time samples; t P，N is the first arrival of the P-wave of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or the number of time samples, and t0 is the earthquake origin time, with the unit of s, ms or the number of time samples;

[0023] The first arrival of the S-band data is picked up, and its arrival first arrival time can be expressed as Equation 5,

[0024] T S =(t S，1 ,..., t S，i ,..., t S，N , t0) Equation 5

[0025] Among them, T S is the S-wave first arrival picked from N gather, with the unit of s, ms or time sample points; t S，1 is the S-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t S，i is the S-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t S，N is the S-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points.

[0026] In an exemplary embodiment of one aspect of the present invention, the realizing the initial positioning of the earthquake source by using the first arrival information of seismic waves in the target area may include:

[0027] For the seismic wave data that can obtain both P-wave and S-wave, use Equation 6 to calculate the initial coordinates of the earthquake source to realize the initial positioning of the earthquake source;

[0028]

[0029] Among them, S(χ S ) is the time residual, with the unit of s, ms or sampling points; χ S is the initial coordinates of the earthquake source, with the unit of m or km; t S，i，c is the calculated travel time of the S-wave from the possible earthquake source to the i-th observation point, with the unit of s, ms or sampling points; t P，i，c is the calculated travel time of the P-wave from the possible earthquake source to the i-th observation point, with the unit of s, ms or sampling points; t S，i is the actual first arrival time of the S-wave from the possible earthquake source to the i-th observation point, with the unit of s, ms or sampling points; t P，i is the actual first arrival time of the P-wave from the possible earthquake source to the i-th observation point, with the unit of s, ms or sampling points;

[0030] For the vibration source waveform with only P-wave or S-wave, use Equation 7 to calculate the initial coordinates of the earthquake source to realize the initial positioning of the earthquake source;

[0031]

[0032] Among them, S(χ S ) is the time residual, with the unit of s, ms or sampling points; χ S is the initial coordinates of the earthquake source, with the unit of m or km; t m，i，c is the calculated arrival time of the m-wave from the possible earthquake source to the i-th observation point, with the unit of s, ms or sampling points; t m，i is the actual first arrival time of the m-wave from the possible earthquake source to the i-th observation point, with the unit of s, ms or sampling points; tm，j，c The arrival time of the m-wave from a possible seismic source to the j-th observation point, in units of s, ms, or number of sampling points; t m，j The actual first arrival time of the m-wave from a possible seismic source to the j-th observation point, in units of s, ms, or number of sampling points.

[0033] In an exemplary embodiment of one aspect of the present invention, the initial source location using seismic wave waveforms may include:

[0034] Based on seismic wave dynamics and kinematic parameters, the initial source location is achieved through Equation 8,

[0035]

[0036] where S E (χ S ) is the waveform vector / scalar energy parameter, which can be expressed as the vector relative energy value or the scalar relative energy value, dimensionless; q is an exponent, dimensionless; t P，i is the start time of the P-wave at the i-th observation point, in units of s, ms, or number of sampling points; is the end time of the P-wave at the i-th observation point, in units of s, ms, or number of sampling points; t S，i is the start time of the S-wave at the i-th observation point, in units of s, ms, or number of sampling points; is the end time of the S-wave at the i-th observation point, in units of s, ms, or number of sampling points; represents the relative S-wave propagation time between the possible seismic source χ S and the observation point χ R，i ; represents the relative P-wave propagation time between the possible seismic source χ S and the observation point χ R，i ; A i and B i are the seismic wave band amplitude compensation factors due to the seismic wave propagation path from the vibration point to the observation point, dimensionless; a is the weight coefficient for location inversion based on the vector / scalar superposition of P-wave seismic wave band data, dimensionless; b is the weight coefficient for location inversion based on the vector / scalar superposition of S-wave seismic wave band data, dimensionless.

[0037] In an exemplary embodiment of one aspect of the present invention, the initial source location using seismic wave waveforms may further include the step of performing a cross-correlation once:

[0038] During the inversion process, according to the relationship between each possible source coordinate and the observation point, the relative arrival time is calculated to obtain reliable band data of the original seismic wave. W observation points are selected from N observation points for cross-correlation. When the proportion of the pairwise cross-correlation values that reach the threshold value among the W observation points exceeds a certain ratio, it indicates that the source coordinate is a reliable source coordinate. Among them, the first cross-correlation value is calculated by Equation 9,

[0039]

[0040] where val i，j is the correlation coefficient of the first cross-correlation, dimensionless; represents the full waveform of the χ R，i th observation of the m-wave, with the unit of m, cm, mm, or μm;

[0041] represents the full waveform of the χ R，j th observation of the m-wave, with the unit of m, cm, mm, or μm.

[0042] In an exemplary embodiment of one aspect of the present invention, the method may further include the step of performing a second cross-correlation:

[0043] During the inversion process, according to the relationship between each possible source coordinate and the observation point, the relative arrival time is calculated to obtain reliable band data of the original seismic wave. W observation points are selected from N observation points for cross-correlation. When the proportion of the pairwise cross-correlation values that reach the new threshold value among the W observation points exceeds a certain ratio, it indicates that the source coordinate χ S , the source mechanism M, and V m (x, y, z) are mutually restricted and are better source parameters that satisfy the current source characteristics. V m (x, y, z) is a more suitable m-wave velocity model for the current source. Among them, the second cross-correlation value is calculated by Equation 10,

[0044]

[0045] where is the correlation coefficient of the second cross-correlation, dimensionless; represents the new full waveform of the m-wave observation point at χ R，i , with the unit of m, cm, mm, or μm;

[0046] represents the new full waveform of the m-wave observation point at χ R，j , with the unit of m, cm, mm, or μm.

[0047] On the other hand, the present invention provides a computer device, including:

[0048] At least one processor;

[0049] A memory storing program instructions, wherein the program instructions are configured to be executed by the at least one processor, and the program instructions include instructions for executing the method according to any one of the above.

[0050] In another aspect of the present invention, there is provided a computer-readable storage medium storing computer program instructions, and when the computer program instructions are executed by a processor, the method according to any one of the above is implemented.

[0051] Compared with the prior art, the beneficial effects of the present invention include at least one of the following:

[0052] (1) The method for inverting composite seismic source parameters of the present invention uses the seismic wave signals observed and recorded by seismic equipment, and utilizes information such as seismic wave amplitude, polarity, and first arrival to realize the automatic detection and positioning of the vibration source, thereby providing important technical support for geological exploration resource exploration, engineering construction risk assessment, geological disaster early warning, etc.;

[0053] (2) The method for inverting composite seismic source parameters of the present invention can be applied to the response of underground fracture activities caused by artificial injection and production such as hydraulic fracturing, deep geothermal exploitation, mine exploitation, carbon dioxide geological storage, gas storage injection and production, and wastewater reinjection, or the fault characterization of natural earthquakes;

[0054] (3) The method for inverting composite seismic source parameters of the present invention can also be applied to the detection and positioning of active seismic sources in oil and gas seismic exploration, the detection and positioning of noise sources, and the detection and positioning of related vibration sources such as geological disasters, etc., and has broad prospects for engineering and industrial technology applications and scientific research. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Through the following description in conjunction with the drawings, the above and other objects and / or features of the present invention will become clearer, wherein:

[0056] Figure 1 A flowchart showing the method for inverting composite seismic source parameters according to an exemplary embodiment of the present invention is shown. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0057] In the following, the method for inverting composite seismic source parameters, computer equipment, and readable storage medium of the present invention will be described in detail in conjunction with exemplary embodiments.

[0058] Figure 1 A flowchart showing the method for inverting composite seismic source parameters according to an exemplary embodiment of the present invention is shown.

[0059] In the first exemplary embodiment of the present invention, the method for inverting composite seismic source parameters includes the steps:

[0060] Obtain the original seismic wave band data of the target area, pick the first arrivals of the original seismic wave band data, and obtain the arrival first arrival time. Here, by performing preprocessing such as noise attenuation and event recognition on the seismic data of the target area, the band data of the seismic waves of the target discrimination element is obtained. Among them, the band data m wave of the seismic wave may include P wave and / or S wave.

[0061] Obtain the initial velocity model of the seismic wave in this area according to the well logging data or past geological data of the target area. Here, the initial velocity model of the seismic wave can be a horizontal layered velocity model or a three-dimensional grid velocity model. For example, the migration stacking velocity model, etc., which can be expressed as V m，0 (x, y, z).

[0062] Use the first arrival information of the seismic wave in the target area and / or the seismic wave waveform to achieve the initial source location. Here, the initial source location can be achieved by using the first arrival information of the seismic wave in the target area, or by using the seismic wave waveform in the target area. Of course, both can also be used together for the initial source location. However, the present invention is not limited to this, and other source location methods can also be used.

[0063] Calculate the dip polarization angle and azimuth polarization angle of the three components of the original seismic wave band data according to the initial source location result and the observation point coordinates, and use the dip polarization angle and azimuth polarization angle to perform vector synthesis on the original seismic wave to obtain new seismic wave band data. Here, using the dip polarization angle and azimuth polarization angle to perform vector synthesis on the original seismic wave to obtain new seismic wave band data can adopt the conventional synthesis methods in the art.

[0064] Finally, perform joint inversion according to the above-obtained new seismic wave band data and the initial velocity model through Equation 1,

[0065]

[0066] where q is an exponent, dimensionless; is the P-wave end time of the i-th observation point, with the unit of s, ms or sampling points; is the SV-wave end time of the i-th observation point, with the unit of s, ms or sampling points; is the SH-wave end time of the i-th observation point, with the unit of s, ms or sampling points; represents the possible source χ S and the observation point χ R，i the relative propagation time of the SV wave, represents the possible source χ S and the observation point χ R，i the relative propagation time of the SH wave, represents the possible source χ S and the observation point χ R，iThe relative arrival time of P-waves; a is the weighting coefficient for location inversion based on scalar / vector superposition of P-wave seismic band data, dimensionless; b is the weighting coefficient for location inversion based on scalar / vector superposition of SV-wave seismic band data, dimensionless; c is the weighting coefficient for location inversion based on scalar / vector superposition of SH-wave seismic band data, dimensionless; M is the seismic source mechanism, which is a 3×3 matrix and can be decomposed into an isotropic component ISO part, a double couple component DC part, and a compensated linear vector dipole CLVD through moment tensor. Among them, the DC part can be represented by azimuth Azi, dip angle Dip, and rake angle Rake, and the units of these three angles are degrees. V(x, y, z) is the new velocity model, with the unit of m / s or km / s; it includes P-waves, SV-waves, and SH-waves, that is, it can also be expressed as V m represents P-wave, SV-wave, and SH-wave respectively; represents the P-wave amplitude and polarity compensation factor at the observation point χ S , under the action of the velocity model V P (x, y, z), according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ; represents the SV-wave amplitude and polarity compensation factor at the observation point χ S , under the action of the velocity model V SV (x, y, z), according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ; represents the SH-wave amplitude and polarity compensation factor at the observation point χ S , under the action of the velocity model V SH (x, y, z), according to the seismic source mechanism M and the possible seismic source coordinates χ R，i .

[0067] In this exemplary embodiment, the original seismic wave band data may include P-wave band data, S-wave band data, or both P-wave band data and S-wave band data exist, where the S-wave includes SV-wave and / or SH-wave.

[0068] In this exemplary embodiment, the P-wave band data can be expressed as Equation 2,

[0069] S P，k =[S P， 1，k (χ R，1 , t1, t0), …, S P，i，k (χ R，i , t i , t0),..., S P，N，k (χ R，N , t N , t0)] Equation 2

[0070] where S p,kThe P-wave of the seismic wave recorded for the k-th component, with the unit of m, cm, mm or μm; S P，1，k (χ R，1 , t1, t0) is the waveform of the P-wave of the seismic wave of the k-th component obtained when the starting time of the χ R，1 band is t1, with the unit of m, cm, mm or μm; S P，i，k (χ R，i , t i , t0) is the waveform of the P-wave of the seismic wave of the k-th component obtained when the starting time of the χ R，i band is t i obtained, with the unit of m, cm, mm or μm; S P，N，k (χ R，N , t N , t0) is the waveform of the P-wave of the seismic wave of the k-th component obtained when the starting time of the χ R，N band is t N obtained, with the unit of m, cm, mm or μm; t0 is the earthquake origin time, with the unit of s, ms or the number of time samples;

[0071] The S-wave band data can be expressed as Equation 3,

[0072] S S，k =[S S，1，k (χ R，1 , t1, t0), …, S S，i，k (χ R，i , t i , t0),..., S S，N，k (χ R，N , t N , t0)] Equation 3

[0073] where, S S,k is the S-wave of the seismic wave recorded for the k-th component, with the unit of m, cm, mm or μm; S S，1，k (χ R，1 , t1, t0) is the waveform of the S-wave of the seismic wave of the k-th component obtained when the starting time of the χ R，1 band is t1, with the unit of m, cm, mm or μm; S S，i，k (χ R，i , t i , t0) is the waveform of the S-wave of the seismic wave of the k-th component obtained when the starting time of the χ R，i band is t i obtained, with the unit of m, cm, mm or μm; S S，N，k (χ R，N , t N , t0) is the waveform of the S-wave of the seismic wave of the k-th component obtained when the starting time of the χ R，N band is t NThe seismic S-wave waveform of the obtained k-th component, with the unit of m, cm, mm or μm; t0 is the earthquake origin time, with the unit of s, ms or time sample points.

[0074] In this exemplary embodiment, for the P-wave band data, pick the first arrival, and its arrival first arrival time can be expressed as Equation 4,

[0075] T P =(T P，1 ,..., t P，i ,..., t P，N , t0 ) Equation 4

[0076] where T P is the P-wave first arrival picked from N trace gathers, with the unit of s or ms; t P，1 is the P-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t P，i is the P-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t P，N is the P-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points;

[0077] For the S-wave band data, pick the first arrival, and its arrival first arrival time can be expressed as Equation 5,

[0078] T S =(t S，1 ,..., t S，i ,..., t S，N , t0) Equation 5

[0079] where T S is the S-wave first arrival picked from N trace gathers, with the unit of s, ms or time sample points; t S，1 is the S-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t S，i is the S-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t S，N is the S-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points.

[0080] In this exemplary embodiment, the realizing the initial positioning of the earthquake origin by using the first arrival information of the seismic waves in the target area may include:

[0081] For the seismic wave data that can obtain both P-waves and S-waves, use Equation 6 to calculate the initial coordinates of the earthquake origin and realize the initial positioning of the earthquake origin;

[0082]

[0083] Among them, S(χ S ) is the time residual, with the unit of s, ms, or number of sampling points; χ S is the initial coordinate of the seismic source, with the unit of m or km; t S，i，c is the calculated travel time of the S-wave from the possible seismic source to the i-th observation point, with the unit of s, ms, or number of sampling points; t P，i，c is the calculated travel time of the P-wave from the possible seismic source to the i-th observation point, with the unit of s, ms, or number of sampling points; t S，i is the actual first arrival time of the S-wave from the possible seismic source to the i-th observation point, with the unit of s, ms, or number of sampling points; t P，i is the actual first arrival time of the P-wave from the possible seismic source to the i-th observation point, with the unit of s, ms, or number of sampling points;

[0084] For the vibration source waveform with only P-waves or S-waves, the initial coordinate of the seismic source is calculated using Equation 7 to achieve the initial positioning of the seismic source;

[0085]

[0086] Among them, S(χ S ) is the time residual, with the unit of s, ms, or number of sampling points; χ S is the initial coordinate of the seismic source, with the unit of m or km; t m，i，c is the calculated arrival time of the m-wave from the possible seismic source to the i-th observation point, with the unit of s, ms, or number of sampling points; t m，i is the actual first arrival time of the m-wave from the possible seismic source to the i-th observation point, with the unit of s, ms, or number of sampling points; t m，j，c is the calculated arrival time of the m-wave from the possible seismic source to the j-th observation point, with the unit of s, ms, or number of sampling points; t m，j is the actual first arrival time of the m-wave from the possible seismic source to the j-th observation point, with the unit of s, ms, or number of sampling points.

[0087] In this exemplary embodiment, the initial positioning of the seismic source using the seismic wave waveform may include:

[0088] Based on the dynamic and kinematic parameters of the seismic wave, the initial positioning of the seismic source is achieved through Equation 8,

[0089]

[0090] Among them, S E (χ S ) is the waveform vector / scalar energy parameter, which can be expressed as a vector relative energy value or a scalar relative energy value, dimensionless; q is an exponent, dimensionless; t P，i is the start time of the P-wave at the i-th observation point, with the unit of s, ms, or number of sampling points; is the P-wave end time of the i-th observation point, in units of s, ms, or number of sampling points; t S，i is the S-wave start time of the i-th observation point, in units of s, ms, or number of sampling points; is the S-wave end time of the i-th observation point, in units of s, ms, or number of sampling points; S(χ S ) is the residual shown in Equation 5 or Equation 6, in units of s, ms, or number of sampling points; represents the possible source χ S and the relative S-wave propagation time between the observation point χ R，i ; represents the relative P-wave propagation time between the possible source χ S and the observation point χ R，i ; A i and B i are seismic wave band amplitude compensation factors caused by the seismic wave propagation path of the vibration point from the observation point, dimensionless; a is the weight coefficient for location inversion based on scalar / vector superposition of P-wave seismic wave band data, dimensionless; b is the weight coefficient for location inversion based on scalar / vector superposition of S-wave seismic wave band data, dimensionless.

[0091] In this exemplary embodiment, the method for realizing the initial location of the source by using the seismic wave waveform may further include the step of performing a first cross-correlation:

[0092] During the inversion process, according to the relationship between each possible source coordinate and the observation point, calculate the relative arrival time, so as to obtain reliable band data of the original seismic wave. Select W observation points from N observation points for cross-correlation. When the proportion of the pairwise cross-correlation values reaching the threshold value among the W observation points exceeds a certain ratio, it indicates that the source coordinate is a reliable source coordinate. Among them, the first cross-correlation value is calculated by Equation 9,

[0093]

[0094] where val i，j is the correlation coefficient of the first cross-correlation, dimensionless; represents the full waveform of the m-th wave at the χ R，i -th observation, in units of m, cm, mm, or μm;

[0095] represents the full waveform of the m-th wave at the χ R，j -th observation, in units of m, cm, mm, or μm.

[0096] In this exemplary embodiment, the method may further include the step of performing a second cross-correlation:

[0097] During the inversion process, according to the relationship between each possible source coordinate and the observation point, the relative arrival time is calculated to obtain reliable band data of the original seismic wave. W observation points are selected from N observation points for cross-correlation. When the proportion of the pairwise cross-correlation values reaching a new threshold value among the W observation points exceeds a certain ratio, it indicates the source coordinate χ S , the source mechanism M and V m (x, y, z) are mutually restricted and are better source parameters that meet the current source characteristics. V m (x, y, z) is an appropriate m-wave velocity model for the current source. Among them, the quadratic cross-correlation value is calculated by formula 10,

[0098]

[0099] Among them, is the correlation coefficient of the quadratic cross-correlation, dimensionless; represents the new full waveform of the m-wave observation point as χ R，i , with the unit of m, cm, mm or μm;

[0100] represents the new full waveform of the m-wave observation point as χ R，j , with the unit of m, cm, mm or μm.

[0101] In the second exemplary embodiment of the present invention, as shown in Figure 1 , the method for inverting the composite source parameters includes the steps:

[0102] (1) Acquisition of seismic exploration band data

[0103] Through preprocessing such as noise attenuation and event recognition of the seismic data in the target area, the band data of the seismic wave is obtained. Among them, the band data of the seismic wave m-wave includes P-wave, S-wave or both P-wave and S-wave. Here, the P-wave band data is expressed as:

[0104] S P，k =[S P，1，k (χ R，1 , t1, t0),..., S P，i，k (χ R，i , t i , t0),…, S P，N，k (χ R，N , t N , t0)

[0105] Among them, S p,k is the P-wave of the seismic wave recorded by the k-th component, with the unit of m, cm, mm or μm; S P，1，k (χ R，1 , t1, t0) is the observation point located at χ R，1The seismic P-wave waveform of the k-th component obtained at the starting time t1 of the frequency band, with the unit of m, cm, mm, or μm; S P，i，k (χ R，i , t i , t0) represents the observation point located at χ R，i The seismic P-wave waveform of the k-th component obtained at the starting time t i of the frequency band, with the unit of m, cm, mm, or μm; S P，N，k (χ R，N , t N , t0) represents the observation point located at χ R，N The seismic P-wave waveform of the k-th component obtained at the starting time t N of the frequency band, with the unit of m, cm, mm, or μm; t0 is the earthquake origin time, with the unit of s, ms, or time sample points.

[0106] Similarly, the S-wave frequency band data can be expressed as:

[0107] S S，k = [S S，1，k (χ R，1 , t1, t0),..., S S，i，k (χ R，i , t i , t0),..., S s，N，k (χ R，N , t N , t0)]

[0108] Among them, S S,k is the seismic S-wave of the k-th component record, with the unit of m, cm, mm, or μm; S S，1，k (χ R，1 , t1, t0) represents the seismic S-wave waveform of the k-th component obtained at the starting time t1 when the observation point is located at χ R，1 , with the unit of m, cm, mm, or μm; S S，i，k (χ R，i , t i , t0) represents the seismic S-wave waveform of the k-th component obtained at the starting time t R，i when the observation point is located at χ i , with the unit of m, cm, mm, or μm; S S，N，k (χ R，N , t N , t0) represents the seismic S-wave waveform of the k-th component obtained at the starting time t R，N when the observation point is located at χ N , with the unit of m, cm, mm, or μm; t0 is the earthquake origin time, with the unit of s, ms, or time sample points.

[0109] Pick the first arrival of the P-wave, and its arrival first arrival time is expressed as:

[0110] T P =(t P，1 ,...,t P，i ,...,t P，N ,t0)

[0111] Among them, T P is the P-wave first arrival picked from N gather, with the unit of s or ms; t P，1 is the P-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t P，i is the P-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t P，N is the P-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points.

[0112] For the earthquake source location based on the dynamic and kinematic parameters of seismic waves, pick the first arrival of S waves, and its arrival first arrival time is expressed as:

[0113] T S =(t S，1 ,...,t S，i ,...,t S，N ,t0)

[0114] Among them, T S is the S-wave first arrival picked from N gather, with the unit of s, ms or time sample points; t S，1 is the S-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t S，i is the S-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t S，N is the S-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points.

[0115] (2) Establishment of the initial velocity model

[0116] Establish an initial velocity model according to well logging data or previous geological data. Here, the initial velocity model can be a horizontal layered velocity model or a three-dimensional grid velocity model. For example, an offset stacking velocity model, etc., which can be expressed as v m，0 (x, y, z), where m represents P-wave or S-wave. If m is P-wave, the P-wave velocity model is expressed as v P，0 (x, y, z), then the S-wave velocity model is expressed as v S，0 (x, y, z).

[0117] (3) Inversion of source parameters

[0118] a. Initial source location using first arrival information

[0119] For seismic data where both P - waves and S - waves can be obtained, the initial source coordinates are calculated using the following formula to achieve initial source location:

[0120]

[0121] where S(χ S ) is the time residual, with the unit of s, ms or sampling points; χ S is the initial source coordinate, with the unit of m or km; t S，i，c is the calculated travel time of the S - wave from the possible source to the i - th observation point, with the unit of s, ms or sampling points; t P，i，c is the calculated travel time of the P - wave from the possible source to the i - th observation point, with the unit of s, ms or sampling points; t S，i is the actual first arrival time of the S - wave from the possible source to the i - th observation point, with the unit of s, ms or sampling points; t P，i is the actual first arrival time of the P - wave from the possible source to the i - th observation point, with the unit of s, ms or sampling points.

[0122] For the vibration source waveform with only P - waves or S - waves, the initial source coordinates are calculated using the following formula to achieve initial source location:

[0123]

[0124] where S(χ S ) is the time residual, with the unit of s, ms or sampling points; χ S is the initial source coordinate, with the unit of m or km; t m，i，c is the calculated arrival time of the m - wave from the possible source to the i - th observation point, with the unit of s, ms or sampling points; t m，i is the actual first arrival time of the m - wave from the possible source to the i - th observation point, with the unit of s, ms or sampling points; t m，j，c is the calculated arrival time of the m - wave from the possible source to the j - th observation point, with the unit of s, ms or sampling points; t m，j is the actual first arrival time of the m - wave from the possible source to the j - th observation point, with the unit of s, ms or sampling points.

[0125] b. Initial source location based on seismic wave waveform characteristics

[0126] For the source location method based on the dynamic and kinematic parameters of seismic waves, the initial source location equations are as follows:

[0127]

[0128] Among them, S E (χ S ) is the waveform vector / scalar energy parameter, which can be expressed as the vector relative energy value or the scalar relative energy value, dimensionless; q is the exponent, dimensionless; t P，i is the P-wave start time of the i-th observation point, with the unit of s, ms or the number of sampling points; is the P-wave end time of the i-th observation point, with the unit of s, ms or the number of sampling points; t S，i is the S-wave start time of the i-th observation point, with the unit of s, ms or the number of sampling points; is the S-wave end time of the i-th observation point, with the unit of s, ms or the number of sampling points; represents the relative S-wave propagation time between the possible seismic source χ S and the observation point χ R，i ; represents the relative P-wave propagation time between the possible seismic source χ S and the observation point χ R，i ; A i and B i are the differences in seismic wave band amplitudes caused by the seismic wave propagation paths of the vibration points from the observation points, dimensionless; a is the weight coefficient for location inversion based on the superposition of P-wave seismic wave band data, dimensionless; b is the weight coefficient for location inversion based on the superposition of S-wave seismic wave band data, dimensionless.

[0129] During the inversion process, according to the relationship between each possible seismic source coordinate χ S and the observation points, calculate the relative arrival time For the S-wave, it is For the P-wave, it is Thus, reliable waveform data of the m-wave can be obtained For the P-wave, the reliable waveform data is For the S-wave, the reliable waveform data is Select W observation points from N observation points and perform one cross-correlation through the following formula:

[0130]

[0131] Among them, val i，j is the correlation coefficient of one cross-correlation, dimensionless; represents the full waveform of the χ R，i -th observation of the m-wave, with the unit of m, cm, mm or μm;

[0132] represents the full waveform of the χ R，j -th observation of the m-wave, with the unit of m, cm, mm or μm. When val i，jThe correlation coefficient reaches the threshold value, and the pairwise cross-correlation quantity of P-waves or S-waves, or S-waves and S-waves reaches a certain proportion. For example, 70%, and this value can be set according to data quality and the number of observation points, then it indicates that χ S is a relatively reliable hypocenter coordinate.

[0133] c. Calculate the inversion first arrival according to the hypocenter position of the seismic wave, and accurately obtain the seismic wave band data according to the first arrival information

[0134] Calculate the source parameters according to the preliminary hypocenter positioning result. According to the preliminary positioning χ S and the observation point coordinates χ R，j calculate the seismic wave S P，i，k (χ R，i , t i， t0) of the dip polarization angle θ i and the azimuth polarization angle φ i .

[0135] Use the dip polarization angle θ i and the azimuth polarization angle φ i , and perform vector synthesis on the original seismic wave S m，i，k (χ R，i , t i , t0) so that the P-wave radial components of these three components point along the line connecting the vibration source χ S and the observation point χ R，j . Obtain the new three-component seismic wave When m is a P-wave, the waveform of the vector-synthesized P-wave is When m is an SV-wave, the waveform of the vector-synthesized SV-wave is When m is an SH-wave, the waveform of the vector-synthesized SV-wave is d. Generate new seismic wave data according to the preliminary hypocenter positioning result Establish an inversion operator according to the following formula, and carry out joint inversion of the source mechanism, source velocity model and source location.

[0136]

[0137] Among them, q is an exponent, dimensionless; is the end time of the P-wave at the i-th observation point, with the unit of s, ms or the number of sampling points; is the end time of the SV-wave at the i-th observation point, with the unit of s, ms or the number of sampling points; is the end time of the SH-wave at the i-th observation point, with the unit of s, ms or the number of sampling points; represents the relative propagation time of the SV-wave between the possible source χ S and the observation point χ R，i . represents the possible source χ SThe relative time of SH wave propagation with respect to the observation point χ R，i ; represents the relative time of P wave propagation between the possible seismic source χ S and the observation point χ R，i ; a is the weighting coefficient for location inversion based on scalar / vector superposition of P-wave seismic band data, dimensionless; b is the weighting coefficient for location inversion based on scalar / vector superposition of SV-wave seismic band data, dimensionless; c is the weighting coefficient for location inversion based on scalar / vector superposition of SH-wave seismic band data, dimensionless; M is the seismic source mechanism, which is a 3×3 matrix and can be decomposed into an isotropic component ISO part, a double couple component DC part, and a compensated linear vector dipole CLVD through moment tensor decomposition. Among them, the DC part can be represented by the azimuth Azi, dip angle Dip, and rake angle Rake, and the units of these three angles are degrees; V(x, y, z) is the new velocity model, with the unit of m / s or km / s; it includes P-wave, SV-wave, and SH-wave, that is, it can also be expressed as V m (x, y, z); m represents P-wave, SV-wave, and SH-wave respectively; represents the P-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V P (x, y, z) according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ; represents the SV-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V SV (x, y, z) according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ; represents the SH-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V SH (x, y, z) according to the seismic source mechanism M and the possible seismic source coordinates χ R，i ;

[0138] Similarly, in the inversion process, according to the relationship between each possible seismic source coordinate χ S and the observation point, the relative arrival time is calculated to obtain reliable waveform data of the m-wave For P-wave, the reliable waveform data is For SV-wave, the reliable waveform data is For SH-wave, the reliable waveform data is Select W observation points from N observation points for cross-correlation. Select W observation points from N observation points for secondary cross-correlation according to the following formula:

[0139]

[0140] Among them, is the correlation coefficient of the second-order cross-correlation, dimensionless; represents the new full waveform at the m-wave observation point χ R，i with the unit of m, cm, mm or μm; represents the new full waveform at the m-wave observation point χR,j, with the unit of m, cm, mm or μm. When the vali,j correlation coefficient reaches a new threshold value, determined according to the data quality, such as 0.7, 0.8, 0.9, etc.; and the pairwise cross-correlation quantity of P-wave or S-wave, or S-wave and S-wave reaches a certain proportion, it indicates that the source coordinates χ S , the source mechanism M and V m (x, y, z) are mutually restricted and are better source parameters that satisfy the current source characteristics, that is, the current source coordinates χ S , the source mechanism M is a relatively reliable source coordinate and source mechanism solution, V m (x, y, z is a suitable m-wave velocity model for the current source.

[0141] In the third exemplary embodiment of the present invention, a computer device is provided, including:

[0142] At least one processor; a memory storing program instructions. Among them, the program instructions are configured to be executed by the at least one processor, and the program instructions include instructions for executing the method according to the first or second exemplary embodiment as above.

[0143] In the fourth exemplary embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method according to the first or second exemplary embodiment as above is implemented.

[0144] In summary, the beneficial effects of the present invention include at least one of the following:

[0145] (1) The method for inverting composite source parameters of the present invention uses the seismic wave signals recorded by seismic equipment, and utilizes information such as seismic wave amplitude, polarity, and first arrival to realize the automatic detection and positioning of the vibration source, thereby providing important technical support for geological exploration resource exploration, engineering construction risk assessment, geological disaster early warning, etc.;

[0146] (2) The method for inverting composite source parameters of the present invention can be applied to the response of underground fracture activities caused by artificial injection and production such as hydraulic fracturing, deep geothermal exploitation, mine exploitation, carbon dioxide geological storage, gas storage injection and production, wastewater reinjection, etc., or the fault characterization of natural earthquakes;

[0147] (3) The method for inverting the composite seismic source parameters of the present invention can also be applied to active seismic source detection in oil and gas seismic exploration, noise source detection and positioning, detection and positioning of related vibration sources such as geological disasters, etc., and has broad prospects for engineering and industrial technology applications and scientific research.

[0148] Although the present invention has been described above in connection with exemplary embodiments and the accompanying drawings, those of ordinary skill in the art should understand that various modifications can be made to the above embodiments without departing from the spirit and scope of the claims.

Claims

1. A method for inverting composite seismic source parameters, characterized in that The method includes the steps of: Obtaining the original seismic wave band data of the target area, picking the first arrivals of the original seismic wave band data, and obtaining the arrival time of the first arrivals; Obtaining the initial seismic wave velocity model of the area according to the logging data or past geological data of the target area; Using the first arrival information of the seismic waves in the target area and / or the seismic wave waveforms to achieve the initial positioning of the seismic source; Calculating the dip polarization angle and azimuth polarization angle of the three components of the original seismic wave band data according to the initial positioning result of the seismic source and the coordinates of the observation points, and using the dip polarization angle and azimuth polarization angle to perform vector synthesis on the original seismic waves to obtain new seismic wave band data; Performing joint inversion on the new seismic wave band data and the initial velocity model through Equation 1; where q is an exponent and dimensionless; is the end time of the P-wave at the i-th observation point, with the unit of s, ms or the number of sampling points; is the end time of the SV-wave at the i-th observation point, with the unit of s, ms or the number of sampling points; is the end time of the SH-wave at the i-th observation point, with the unit of s, ms or the number of sampling points; represents the possible source χ S and the relative propagation time of the SV-wave to the observation point χ R，i ; represents the possible source χ S and the relative propagation time of the SH-wave to the observation point χ R，i ; represents the relative propagation time of the P-wave between the possible source χ S and the observation point χ R，i ; a is the weighting coefficient for location inversion based on scalar / vector superposition of P-wave seismic band data, dimensionless; b is the weighting coefficient for location inversion based on scalar / vector superposition of SV-wave seismic band data, dimensionless; c is the weighting coefficient for location inversion based on scalar / vector superposition of SH-wave seismic band data, dimensionless; M is the seismic source mechanism, which is a 3×3 matrix; V(x, y, z) is the new velocity model, with the unit of m / s or km / s; represents the P-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V P (x, y, z) according to the source mechanism M and the possible source coordinates χ R，i ; represents the SV-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V SV (x, y, z) according to the source mechanism M and the possible source coordinates χ R，i ; represents the SH-wave amplitude and polarity compensation factor of the observation point χ S under the action of the velocity model V SH (x, y, z) according to the source mechanism M and the possible source coordinates χ R，i .

2. The method for inverting composite seismic source parameters according to claim 1, wherein The original seismic wave band data includes P-wave band data, S-wave band data, or both P-wave band data and S-wave band data exist, where the S-wave includes SV-wave and / or SH-wave.

3. The method for inverting composite seismic source parameters according to claim 2, characterized in that The P-wave band data is expressed as Equation 2; S P，k = [S P，1，k (χ R，1 , t1, t0),..., S P，i，k (χ R，i , t i , t0),..., S P，N，k (χ R，N , t N , t0)] Equation 2 Among them, S p，k is the P-wave of the seismic wave recorded for the k-th component, with the unit of m, cm, mm or μm; S P，1，k (χ R，1 , t1, t0) is the waveform of the seismic P-wave of the k-th component obtained when the observation point is located at χ R，1 and the starting time of the wave band is t1, with the unit of m, cm, mm or μm; S P，i，k (χ R，i , t i , t0) is the waveform of the seismic P-wave of the k-th component obtained when the observation point is located at χ R，i and the starting time of the wave band is t i obtained, with the unit of m, cm, mm or μm; S P，N，k (χ R，N , t N , t0) is the waveform of the seismic P-wave of the k-th component obtained when the observation point is located at χ R，N and the starting time of the wave band is t N obtained, with the unit of m, cm, mm or μm; t0 is the earthquake origin time, with the unit of s, ms or the number of time samples; The S-wave band data is expressed as Equation 3; S S，k = [S S，1，k (χ R，1 , t1, t0),..., S S，i，k (χ R，i , t i , t0),..., S S，N，k (χ R，N , t N , t0)] Equation 3 Among them, S S，k is the S-wave of the seismic wave recorded for the k-th component, with the unit of m, cm, mm or μm; S S，1，k (χ R，1 , t1, t0) is the waveform of the seismic S-wave of the k-th component obtained when the starting time of the χ R，1 band is t1, with the unit of m, cm, mm or μm; S S，i，k (χ R，i , t i , t0) is the waveform of the seismic S-wave of the k-th component obtained when the starting time of the χ R，i band is t i obtained, with the unit of m, cm, mm or μm; S S，N，k (χ R，N , t N , t0) is the waveform of the seismic S-wave of the k-th component obtained when the starting time of the χ R，N band is t N obtained, with the unit of m, cm, mm or μm; t0 is the earthquake origin time, with the unit of s, ms or the number of time samples.

4. The method for inverting composite seismic source parameters according to claim 2 or 3, characterized in that, Picking the first arrivals of the P-wave band data, and its arrival time of the first arrivals is expressed as Equation 4; T P =(t P，1 ,..., t P，i ,..., L P，N , t0) Equation 4 Among them, T P is the P-wave first arrival picked from N gather, with the unit of s or ms; t P，1 is the P-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t P，i is the P-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t P，N is the P-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points; Picking the first arrivals of the S-band data, and its arrival time of the first arrivals is expressed as Equation 5; T S =(t S，1 ,...,t S，i ,...,t S，N ,t0) Equation 5 Among them, T S is the S-wave first arrival picked up from N gather, with the unit of s, ms or time sample points; t S，1 is the S-wave first arrival of the seismic wave data recorded at the first observation point, with the unit of s, ms or time sample points; t S，i is the S-wave first arrival of the seismic wave data recorded at the i-th observation point, with the unit of s, ms or time sample points; t S，N is the S-wave first arrival of the seismic wave data recorded at the N-th observation point, with the unit of s, ms or time sample points, and t0 is the earthquake origin time, with the unit of s, ms or time sample points.

5. The method for inverting composite seismic source parameters according to claim 1, characterized in that The realizing the initial positioning of the seismic source by using the first arrival information of the seismic waves in the target area includes: For the seismic wave data for which both P-waves and S-waves can be obtained, calculating the initial coordinates of the seismic source by using Equation 6 to achieve the initial positioning of the seismic source; where S(χ S ) is the time residual, with the unit of s, ms or the number of sampling points; χ S is the initial coordinates of the seismic source, with the unit of m or km; t S，i，c is the calculated travel time of the S-wave from the possible seismic source to the i-th observation point, with the unit of s, ms or the number of sampling points; t P，i，c is the calculated travel time of the P-wave from the possible seismic source to the i-th observation point, with the unit of s, ms or the number of sampling points; t S，i is the actual first arrival time of the S-wave from the possible seismic source to the i-th observation point, with the unit of s, ms or the number of sampling points; t P，i is the actual first arrival time of the P-wave from the possible seismic source to the i-th observation point, with the unit of s, ms or the number of sampling points; For the vibration source waveforms with only P-waves or S-waves, calculating the initial coordinates of the seismic source by using Equation 7 to achieve the initial positioning of the seismic source; Among them, S(χ S ) is the time residual, with the unit of s, ms or the number of sampling points; χ S is the initial coordinate of the seismic source, with the unit of m or km; t m，i，c is the calculated arrival time of the m-wave from the possible seismic source to the i-th observation point, with the unit of s, ms or the number of sampling points; t m，i is the actual first arrival time of the m-wave from the possible seismic source to the i-th observation point, with the unit of s, ms or the number of sampling points; t m，j，c is the calculated arrival time of the m-wave from the possible seismic source to the j-th observation point, with the unit of s, ms or the number of sampling points; t m，j is the actual first arrival time of the m-wave from the possible seismic source to the j-th observation point, with the unit of s, ms or the number of sampling points.

6. The method for inverting composite seismic source parameters according to claim 1, characterized in that The realizing the initial positioning of the seismic source by using the seismic wave waveforms includes: Based on the dynamic and kinematic parameters of the seismic waves, achieving the initial positioning of the seismic source through Equation 8; Among them, S E (χ S ) is the waveform vector / scalar energy parameter, which can be expressed as the vector relative energy value or the scalar relative energy value, dimensionless; q is the exponent, dimensionless; t P，i is the P-wave start time of the i-th observation point, with the unit of s, ms or sampling points; is the P-wave end time of the i-th observation point, with the unit of s, ms or sampling points; t S，i is the S-wave start time of the i-th observation point, with the unit of s, ms or sampling points; is the S-wave end time of the i-th observation point, with the unit of s, ms or sampling points; represents the possible seismic source χ S and the relative S-wave propagation time between the possible seismic source χ R，i and the observation point χ represents the possible seismic source χ S and the relative P-wave propagation time between the possible seismic source χ R，i and the observation point χ; A i and B i are the seismic wave band amplitude compensation factors caused by the seismic wave propagation path of the vibration point to the observation point, dimensionless; a is the weight coefficient for location inversion based on the vector / scalar superposition of P-wave seismic wave band data, dimensionless; b is the weight coefficient for location inversion based on the vector / scalar superposition of S-wave seismic wave band data, dimensionless.

7. The method for inverting composite seismic source parameters according to claim 6, characterized in that, The realizing the initial positioning of the seismic source by using the seismic wave waveforms further includes the step of performing one cross-correlation: In the inversion process, according to the relationship between each possible seismic source coordinate and the observation point, calculating the relative arrival time, so as to obtain reliable band data of the original seismic waves, selecting W observation points from N observation points for cross-correlation. When the proportion of the cross-correlation values reaching the threshold value among the pairwise cross-correlations of the W observation points exceeds a certain proportion, it indicates that the seismic source coordinate is a reliable seismic source coordinate, where the one cross-correlation value is calculated by Equation 9; Among them, val i，j is the correlation coefficient of the first cross-correlation, dimensionless; represents the full waveform of the χ R，i th observation of the m-wave, with the unit of m, cm, mm or μm; Indicates the full waveform of the χth observation of the m-wave R，j in units of m, cm, mm, or μm.

8. The method for inverting composite seismic source parameters according to claim 1, characterized in that The method further includes the step of performing two cross-correlations: During the inversion process, according to the relationship between each possible source coordinate and the observation point, the relative arrival time is calculated to obtain reliable band data of the original seismic wave. W observation points are selected from N observation points for cross-correlation. When the proportion of the pairwise cross-correlation values among the W observation points that reach a new threshold value exceeds a certain ratio, it indicates the source coordinate χ S , the source mechanism M and V m (x, y, z) are mutually restricted and are better source parameters that satisfy the current source characteristics. V m (x, y, z) is an m-wave velocity model suitable for the current source. Among them, the secondary cross-correlation is calculated by Equation 10 Among them, is the correlation coefficient of the second-order cross-correlation, dimensionless; indicates that the observation point of the m-wave is χ R，i new full waveform, with units of m, cm, mm or μm; Indicates that the m-wave observation point is χ R，j The new full waveform, with units of m, cm, mm, or μm.

9. A computer device, characterized in that, Including: At least one processor; A memory storing program instructions, where the program instructions are configured to be executed by the at least one processor, and the program instructions include instructions for executing the method according to any one of claims 1-8.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, the method according to any one of claims 1-8 is implemented.

Images