Full waveform inversion method based on time-frequency domain amplitude and phase joint

Abstract

This invention provides a full-waveform inversion method based on time-frequency domain amplitude and phase joint analysis, belonging to the field of mine geological structure detection technology. The method includes: S1, establishing an initial smoothed velocity model; S2, preprocessing seismic data; S3, calculating a time-frequency domain amplitude adjustment factor and constructing a time-frequency domain amplitude and phase objective function; S4, performing time-frequency domain amplitude and phase inversion on seismic data envelopment lines based on the objective function and the initial smoothed velocity model to obtain an updated smoothed velocity model; S5, optimizing the updated smoothed velocity model using an improved simulated annealing method to obtain an initial high-precision velocity model; S6, optimizing the seismic inversion algorithm using the initial high-precision velocity model and time-frequency domain amplitude and phase information to obtain the final high-precision velocity model. This invention can fully utilize the amplitude and phase information of seismic data, improve the accuracy and resolution of the inversion, and achieve precise detection of complex deep structures in mines.

Description

Technical Field

[0001] This invention relates to the field of mine geological structure detection technology, and in particular to a full waveform inversion method based on time-frequency domain amplitude-phase joint. Background Technology

[0002] Mine geological structures are key geological factors affecting coal mine production and construction. They are ubiquitous in all coal mines and play a controlling role in mining activities. The complexity of these geological structures not only increases the difficulty of coal mining but also poses a serious threat to mine safety. Traditional methods for detecting mine geological structures, such as geological mapping and drilling, while providing information on geological structures to some extent, often suffer from low detection accuracy and insufficient resolution, making it difficult to meet the detection needs under complex geological conditions.

[0003] With the development of seismic exploration technology, the full-waveform inversion method has gradually become an important tool for mine geological structure detection due to its high accuracy and high resolution. Full-waveform inversion utilizes the frequency, amplitude, and phase information from mine seismic data to obtain a physical property model of the subsurface medium through wave equation inversion, enabling high-precision characterization of subsurface physical parameters. However, most existing full-waveform inversion methods only utilize a single attribute of seismic data (such as amplitude or phase), which makes the inversion results susceptible to noise interference and difficult to comprehensively reflect the geological characteristics of the subsurface medium. Furthermore, for the detection of complex deep coalfield structures, the full-waveform inversion method also faces challenges such as strong dependence on the initial model and lack of low-frequency information.

[0004] Therefore, how to make full use of the amplitude and phase information of seismic data, improve the accuracy and resolution of full waveform inversion, and achieve accurate detection of complex deep coalfield structures has become a hot and difficult issue in the current research on mine geological structure detection technology. Summary of the Invention

[0005] The purpose of this invention is to provide a full waveform inversion method based on the joint amplitude and phase in the time and frequency domains, which can make full use of the amplitude and phase information of seismic data, improve the accuracy and resolution of the inversion, and achieve accurate detection of complex deep coalfield structures.

[0006] To achieve the above objectives, this invention provides a full waveform inversion method based on time-frequency domain amplitude-phase joint analysis, comprising the following steps:

[0007] S1. Input the initial velocity model, source wavelet, and seismic data to establish the initial smooth velocity model;

[0008] S2. Preprocess the seismic data to obtain preprocessed data;

[0009] S3. Calculate the time-frequency domain amplitude adjustment factor using the preprocessed data, and construct the time-frequency domain amplitude-phase objective function;

[0010] S4. Based on the time-frequency domain amplitude-phase objective function and the initial smooth velocity model, perform time-frequency domain amplitude-phase inversion of seismic data envelopment lines to obtain the updated smooth velocity model.

[0011] S5. The updated smooth velocity model is optimized and solved using the improved simulated annealing method, taking into account probability parameters, initial temperature, and final temperature, to obtain an initial high-precision velocity model.

[0012] S6. Optimize the seismic inversion algorithm using the initial high-precision velocity model and time-frequency domain amplitude and phase information to obtain the final high-precision velocity model.

[0013] Preferably, in step S1, the calculation formula for the initial smooth velocity model is as follows:

[0014] V init =V0+α(DE) (1);

[0015] Among them, V init V0 represents the initial smoothed velocity model; D represents the source wavelet; E represents the seismic data; and α represents the initial adjustment factor.

[0016] Preferably, in step S2, the preprocessing includes Wiener filtering scaling, Hilbert transform, and short-time Fourier transform;

[0017] The Wiener filter scale decomposition is as follows:

[0018] E wiener =Wiener(E,σ 2 (2);

[0019] Among them, E wiener This represents the data after Wiener filtering decomposition; σ 2 Indicates the noise power spectral density;

[0020] The Hilbert transform is as follows:

[0021] H = Hilbert(E) Wiener (3);

[0022] Where H represents the result of the Hilbert transform;

[0023] The short-time Fourier transform is as follows:

[0024] S = STFT(E) Wiener ,ω,t) (4);

[0025] Where S represents the short-time Fourier transform result; ω represents the window function; and t represents the time shift parameter.

[0026] Preferably, in step S3, the time-frequency domain amplitude adjustment factor is calculated using the preprocessed data, and the specific operation is as follows:

[0027]

[0028] Among them, A reg This represents the time-frequency domain amplitude adjustment factor.

[0029] Preferably, in step S3, the objective function for amplitude and phase in the time-frequency domain is constructed, specifically as follows:

[0030]

[0031] Where F represents the time-frequency domain amplitude-phase objective function; This represents the phase of the short-time Fourier transform; represents the phase of the Hilbert transform; i represents the imaginary unit.

[0032] Preferably, in step S4, based on the time-frequency domain amplitude-phase objective function and the initial smoothed velocity model, time-frequency domain amplitude-phase inversion of the seismic data envelopment line is performed to obtain the updated smoothed velocity model. Specifically, the operation is as follows:

[0033]

[0034] Among them, V update This represents the updated smooth velocity model; The phase adjustment amount is expressed by the following formula:

[0035]

[0036] Where η represents the learning rate; ▽F represents the gradient of the objective function F with respect to the phase.

[0037] Preferably, in step S5, the initial high-precision velocity model is as follows:

[0038] V high =SA(V update ,P,T init ,T final (9);

[0039] Among them, V high This represents the optimized high-precision velocity model; SA represents the simulated annealing algorithm function; P represents the probability parameter; T init Indicates the initial temperature; T final This indicates the final temperature.

[0040] Preferably, in step S5, the improved simulated annealing method includes the following steps:

[0041] S51. Parameter initialization: Select an initial velocity model and set the initial temperature T. init and final temperature T final And define the energy function E(V);

[0042] The energy function is as follows:

[0043]

[0044] Where V represents the velocity model; ▽V represents the gradient of the velocity model; f represents the frequency; E obs (τ,f) represents the time-frequency representation of the observed seismic data; λ represents the regularization parameter;

[0045] S52, Iterative search, in the current solution V current A new solution V is randomly generated within the neighborhood of V. new Calculate the energy difference ΔE between the current solution and the new solution:

[0046] ΔE=E(V new )-E(V current (11);

[0047] S53. Acceptance Criterion: When ΔE is less than 0, accept the new solution as the current solution; when ΔE is greater than or equal to 0, accept the new solution with probability exp(-ΔE / T), where T represents the current temperature.

[0048] S54. Cooling strategy: Reduce the temperature according to the predetermined cooling rate δ. When the termination condition is reached, output the initial high-precision velocity model V. high .

[0049] Preferably, in step S6, the final high-precision velocity model is as follows:

[0050] V final =V high +β(V update -V high (12);

[0051] Among them, V final β represents the final high-precision speed model; β represents the final adjustment factor.

[0052] Therefore, the present invention employs the above-mentioned full waveform inversion method based on the joint amplitude and phase in the time and frequency domains, and the beneficial technical effects are as follows:

[0053] (1) Improve inversion accuracy and resolution:

[0054] This invention makes full use of the amplitude and phase information of seismic data. Compared with the traditional full waveform inversion method that only uses a single attribute, it can more comprehensively reflect the geological characteristics of the underground medium, thereby improving the accuracy and resolution of the inversion.

[0055] (2) Achieve precise detection of complex structures deep within coalfields:

[0056] This invention, by combining time-frequency domain amplitude and phase information, can more accurately characterize the geological features of complex deep coalfield structures, providing a more reliable geological basis for coal mining.

[0057] (3) Enhanced noise interference resistance:

[0058] In the seismic data preprocessing stage, this invention employs techniques such as Wiener filtering scale decomposition, Hilbert transform, and short-time Fourier transform to effectively reduce noise interference with the inversion results and improve the stability and reliability of the inversion.

[0059] (4) Optimize the solution process:

[0060] This invention utilizes an improved simulated annealing method to optimize and solve the updated smooth velocity model. Through iterative search and acceptance criteria, the optimal solution is gradually approximated, and an initial high-precision velocity model is obtained.

[0061] (5) Algorithm Optimizability:

[0062] After obtaining an initial high-precision velocity model, this invention further optimizes the seismic inversion algorithm using the initial high-precision velocity model and time-frequency domain amplitude and phase information, thereby obtaining a final high-precision velocity model. This optimizability allows the method to continuously adapt to new data and geological conditions in practical applications, improving the inversion results. Attached Figure Description

[0063] Figure 1 This is a flowchart of a full waveform inversion method based on the joint amplitude and phase in the time and frequency domain according to the present invention;

[0064] Figure 2 The flowchart for optimizing the solution using the improved simulated annealing method is shown. Detailed Implementation

[0065] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0066] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning as understood by one of ordinary skill in the art to which this invention pertains.

[0067] Example 1

[0068] like Figure 1 The diagram shows a flowchart of a full waveform inversion method based on time-frequency domain amplitude-phase joint operation according to the present invention, which specifically includes the following steps:

[0069] S1. Input the initial velocity model, source wavelet, and seismic data to establish an initial smooth velocity model, which provides the basis for subsequent inversion processes.

[0070] The calculation formula for the initial smooth velocity model is as follows:

[0071] V init =V0+α(DE) (1);

[0072] Among them, V init V0 represents the initial smoothed velocity model; D represents the source wavelet; E represents the seismic data; and α represents the initial adjustment factor.

[0073] S2. Preprocess the seismic data to obtain preprocessed data.

[0074] Preprocessing includes Wiener filtering scaling, Hilbert transform, and short-time Fourier transform;

[0075] The Wiener filter scale decomposition is as follows:

[0076] E wiener =Wiener(E,σ 2 (2);

[0077] Among them, E wiener This represents the data after Wiener filtering decomposition; σ 2 Indicates the noise power spectral density;

[0078] The Hilbert transform is as follows:

[0079] H = Hilbert(E) Wiener (3);

[0080] Where H represents the result of the Hilbert transform;

[0081] The short-time Fourier transform is as follows:

[0082] S = STFT(E) Wiener ,ω,t) (4);

[0083] Where S represents the short-time Fourier transform result; ω represents the window function; and t represents the time shift parameter.

[0084] S3. Calculate the time-frequency domain amplitude adjustment factor using the preprocessed data, and construct the time-frequency domain amplitude-phase objective function. The specific steps are as follows:

[0085]

[0086] Among them, A reg This represents the time-frequency domain amplitude adjustment factor.

[0087]

[0088] Where F represents the time-frequency domain amplitude-phase objective function; This represents the phase of the short-time Fourier transform; represents the phase of the Hilbert transform; i represents the imaginary unit.

[0089] S4. Based on the time-frequency domain amplitude-phase objective function and the initial smoothed velocity model, perform time-frequency domain amplitude-phase inversion of the seismic data envelopment line to obtain the updated smoothed velocity model:

[0090]

[0091] Among them, V update This represents the updated smooth velocity model; The phase adjustment amount is expressed by the following formula:

[0092]

[0093] Where η represents the learning rate; ▽F represents the gradient of the objective function F with respect to the phase.

[0094] S5. The updated smooth velocity model is optimized and solved using the improved simulated annealing method, taking into account probability parameters, initial temperature, and final temperature, to obtain an initial high-precision velocity model.

[0095] The initial high-precision velocity model is as follows:

[0096] V high =SA(V update ,P,T init ,T final (9);

[0097] Among them, V high This represents the optimized high-precision velocity model; SA represents the simulated annealing algorithm function; P represents the probability parameter; T init Indicates the initial temperature; T final This indicates the final temperature.

[0098] like Figure 2 As shown, the improved simulated annealing method includes the following steps:

[0099] S51. Parameter initialization: Select an initial velocity model and set the initial temperature T. init and final temperature T final And define the energy function E(V);

[0100] The energy function is as follows:

[0101]

[0102] Where V represents the velocity model; ▽V represents the gradient of the velocity model; f represents the frequency; E obs (τ,f) represents the time-frequency representation of the observed seismic data; λ represents the regularization parameter;

[0103] S52, Iterative search, in the current solution V current A new solution V is randomly generated within the neighborhood of V. new Calculate the energy difference ΔE between the current solution and the new solution:

[0104] ΔE=E(V new )-E(V current (11);

[0105] S53. Acceptance Criterion: When ΔE is less than 0, accept the new solution as the current solution; when ΔE is greater than or equal to 0, accept the new solution with probability exp(-ΔE / T), where T represents the current temperature.

[0106] S54. Cooling strategy: Reduce the temperature according to the predetermined cooling rate δ, i.e., T1 = δT. When the termination condition is reached, output the initial high-precision velocity model V. high .

[0107] S6. Optimize the seismic inversion algorithm using the initial high-precision velocity model and time-frequency domain amplitude and phase information to obtain the final high-precision velocity model:

[0108] V final =V high +β(V update -V high (12);

[0109] Among them, V final β represents the final high-precision speed model; β represents the final adjustment factor.

[0110] It is worth noting that all contents not described in detail in this invention are existing technologies and are well known to those skilled in the art.

[0111] Therefore, the present invention employs the above-mentioned full waveform inversion method based on the joint amplitude and phase in the time and frequency domain, which can make full use of the amplitude and phase information of seismic data, improve the accuracy and resolution of the inversion, and achieve accurate detection of complex deep coalfield structures.

[0112] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A full waveform inversion method based on time-frequency domain amplitude-phase joint method, characterized in that, Includes the following steps: S1. Input the initial velocity model, source wavelet, and seismic data to establish the initial smooth velocity model; S2. Preprocess the seismic data to obtain preprocessed data; S3. Calculate the time-frequency domain amplitude adjustment factor using the preprocessed data, and construct the time-frequency domain amplitude-phase objective function; S4. Based on the time-frequency domain amplitude-phase objective function and the initial smooth velocity model, perform time-frequency domain amplitude-phase inversion of seismic data envelopment lines to obtain the updated smooth velocity model. S5. The updated smooth velocity model is optimized and solved using the improved simulated annealing method, taking into account probability parameters, initial temperature, and final temperature, to obtain an initial high-precision velocity model. S6. Optimize the seismic inversion algorithm using the initial high-precision velocity model and time-frequency domain amplitude and phase information to obtain the final high-precision velocity model; In step S1, the calculation formula for the initial smooth velocity model is as follows: （1）； in, This represents the initial smooth velocity model; Indicates the initial velocity model; Indicates the source wavelet; Represents earthquake data; Indicates the initial adjustment factor; In step S2, the preprocessing includes Wiener filtering scaling, Hilbert transform, and short-time Fourier transform; The Wiener filter scale decomposition is as follows: （2）； in, This represents the data after Wiener filtering decomposition; Indicates the noise power spectral density; The Hilbert transform is as follows: （3）； in, Indicate the result of the Hilbert transform; The short-time Fourier transform is as follows: （4）； in, This represents the result of the short-time Fourier transform; Indicates the window function; Indicates the time shift parameter; In step S3, the time-frequency domain amplitude adjustment factor is calculated using the preprocessed data. The specific operation is as follows: （5）； in, This represents the time-frequency domain amplitude adjustment factor; In step S3, the objective function of amplitude and phase in the time-frequency domain is constructed, specifically as follows: （6）； in, This represents the objective function of amplitude and phase in the time-frequency domain. This represents the phase of the short-time Fourier transform; This represents the phase of the Hilbert transform; Represents the imaginary unit; In step S4, based on the time-frequency domain amplitude-phase objective function and the initial smoothed velocity model, time-frequency domain amplitude-phase inversion of the seismic data envelopment line is performed to obtain the updated smoothed velocity model. The specific operation is as follows: （7）； in, This represents the updated smooth velocity model; The phase adjustment amount is expressed by the following formula: （8）； in, Indicates the learning rate; Describe the objective function Gradient with respect to phase.

2. The full waveform inversion method based on time-frequency domain amplitude-phase joint as described in claim 1, characterized in that, In step S5, the initial high-precision velocity model is as follows: （9）； in, This represents the optimized high-precision velocity model; This represents the simulated annealing algorithm function; Represents probability parameters; Indicates the initial temperature; This indicates the final temperature.

3. The full waveform inversion method based on time-frequency domain amplitude-phase joint as described in claim 2, characterized in that, In step S5, the improved simulated annealing method includes the following steps: S51. Parameter initialization: Select an initial velocity model and set the initial temperature. and final temperature And define the energy function ; The energy function is as follows: （10）； in, Represents the velocity model; The gradient represents the velocity model; Indicates frequency; The time-frequency representation of observed seismic data; Represents the regularization parameter; S52, Iterative search, in the current solution A new solution is randomly generated within the neighborhood of . Calculate the energy difference between the current solution and the new solution. : （11）； S53, Acceptance Criteria, when When the value is less than 0, accept the new solution as the current solution. When greater than or equal to 0, with probability Accept the new solution, in which, Indicates the current temperature; S54, Cooling strategy, according to the predetermined cooling rate. Lower the temperature, and when the termination condition is met, output the initial high-precision velocity model. .

4. The full waveform inversion method based on time-frequency domain amplitude-phase joint as described in claim 3, characterized in that, In step S6, the final high-precision velocity model is as follows: （12）； in, This represents the final high-precision velocity model; This represents the final adjustment factor.

Images