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

Abstract

The invention provides a full waveform inversion method based on time-frequency domain amplitude phase combination, and belongs to the technical field of mine geologic structure detection, and the method comprises the steps: S1, building an initial smooth velocity model; s2, preprocessing the seismic data; s3, calculating a time-frequency domain amplitude regulation factor, and constructing a time-frequency domain amplitude phase objective function; s4, performing time-frequency domain amplitude phase inversion of the seismic data envelope line based on the objective function and the initial smoothing velocity model to obtain an updated smoothing velocity model; s5, performing optimization solution on the updated smooth velocity model by using an improved simulated annealing method to obtain an initial high-precision velocity model; and S6, optimizing a seismic inversion algorithm by using the initial high-precision speed model and the time-frequency domain amplitude phase information to obtain a final high-precision speed model. The method can make full use of the amplitude and phase information of seismic data, improves the inversion precision and resolution, and achieves the precise detection of the deep complex structure of a mine.

Description

Technical Field

[0001] The present invention relates to the technical field of mine geological structure detection, and particularly to a full waveform inversion method based on the joint amplitude and phase in the time-frequency domain. Background Technique

[0002] Mine geological structures are key geological factors affecting coal mine production and construction, which are commonly present in each coal mine and play a controlling role in coal mine mining activities. The complexity of these geological structures not only increases the difficulty of coal mine mining, but also poses a serious threat to the safe production of the mine. Traditional mine geological structure detection methods, such as geological mapping, drilling, etc., although can provide information on geological structures to a certain extent, often have problems such as low detection accuracy and insufficient resolution, and are difficult to meet the detection requirements under complex geological conditions.

[0003] With the development of seismic exploration technology, the full waveform inversion method has gradually become an important means for mine geological structure detection due to its characteristics of high precision and high resolution. The full waveform inversion utilizes the frequency, amplitude, and phase information in mine seismic data, and obtains the physical property parameter model of the underground medium through wave equation inversion, which can accurately depict the physical property parameters of the underground. However, most of the existing full waveform inversion methods only utilize a single attribute (such as amplitude or phase) of seismic data, which results in the inversion results being easily interfered by noise and difficult to comprehensively reflect the geological characteristics of the underground medium. In addition, for the detection of complex deep structures in coal fields, 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 structures in coal fields 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 the present invention is to provide a full waveform inversion method based on the joint amplitude and phase in the time-frequency domain, which can make full use of the amplitude and phase information of seismic data, improve the accuracy and resolution of inversion, and achieve accurate detection of complex deep structures in coal fields.

[0006] To achieve the above purpose, the present invention provides a full waveform inversion method based on the joint amplitude and phase in the time-frequency domain, including the following steps:

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

[0008] S2. Preprocess the seismic data to obtain the 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 the seismic data envelope to obtain an updated smooth velocity model;

[0011] S5. Use the improved simulated annealing method to optimize and solve the updated smooth velocity model, considering the probability parameter, 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 the time-frequency domain amplitude-phase information to obtain the final high-precision velocity model.

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

[0014] V init = V0 + α(D - E) (1);

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

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

[0017] Among them, the Wiener filter scale decomposition is as follows:

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

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

[0020] The Hilbert transform is as follows:

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

[0022] Among them, 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] Among them, S represents the short-time Fourier transform result; ω represents the window function; t represents the time translation 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 represents the time-frequency domain amplitude adjustment factor.

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

[0030]

[0031] Among them, F represents the time-frequency domain amplitude-phase objective function; 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 smooth velocity model, the time-frequency domain amplitude-phase inversion of the seismic data envelope is performed to obtain the updated smooth velocity model, and the specific operation is as follows:

[0033]

[0034] Among them, V update represents the updated smooth velocity model; represents the phase adjustment amount, and the calculation formula is as follows:

[0035]

[0036] Among them, η 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 represents the optimized high-precision velocity model; SA represents the simulated annealing algorithm function; P represents the probability parameter; T init represents the initial temperature; T final represents 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, set the initial temperature T init and the 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: Randomly generate a new solution V current in the neighborhood of the current solution V new , and 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 a probability of exp(-ΔE / T), where T represents the current temperature;

[0048] S54. Cooling strategy: Reduce the temperature according to a 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] where V final represents the final high-precision velocity model; β represents the final adjustment factor.

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

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

[0054] The present 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 underground media, thereby improving the accuracy and resolution of inversion.

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

[0056] By combining the amplitude and phase information in the time-frequency domain, the present invention can more accurately characterize the geological characteristics of complex deep structures in coalfields, providing a more reliable geological basis for coal mining.

[0057] (3) Enhance the ability to resist noise interference:

[0058] In the preprocessing stage of seismic data, the present invention adopts technical means such as Wiener filter scale decomposition, Hilbert transform, and short-time Fourier transform, effectively reducing the interference of noise on the inversion result and improving the stability and reliability of inversion.

[0059] (4) Optimize the solution process:

[0060] The present invention uses an improved simulated annealing method to optimize and solve the updated smooth velocity model. Through steps such as iterative search and acceptance criteria, it gradually approaches the optimal solution and obtains an initial high-precision velocity model.

[0061] (5) Algorithm optimizability:

[0062] After obtaining the initial high-precision velocity model, the present invention further optimizes the seismic inversion algorithm using the initial high-precision velocity model and the amplitude-phase information in the time-frequency domain, thereby obtaining the final high-precision velocity model. This optimizability enables the method to continuously adapt to new data and geological conditions in practical applications and improve the inversion effect. Brief Description of the Drawings

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

[0064] Figure 2 is a flowchart of optimizing and solving using an improved simulated annealing method. Detailed Embodiments

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

[0066] Unless otherwise defined, the technical terms or scientific terms used in the present invention shall have the ordinary meanings understood by those of ordinary skill in the field to which the present invention belongs.

[0067] Example 1

[0068] AsFigure 1 As shown in the figure, it is a flowchart of a full waveform inversion method based on the joint amplitude and phase in the time-frequency domain according to the present invention, which specifically includes the following steps:

[0069] S1. Input the initial velocity model, source wavelet and seismic data, establish the initial smooth velocity model, and provide a basis for the subsequent inversion process.

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

[0071] V init = V0 + α(D - E) (1);

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

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

[0074] The preprocessing includes Wiener filter scale decomposition, Hilbert transform and short-time Fourier transform;

[0075] Among them, the Wiener filter scale decomposition is as follows:

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

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

[0078] The Hilbert transform is as follows:

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

[0080] Among them, H represents the Hilbert transform result;

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

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

[0083] Among them, S represents the short-time Fourier transform result; ω represents the window function; t represents the time translation parameter.

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

[0085]

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

[0087]

[0088] Among them, F represents the time-frequency domain amplitude-phase objective function; 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 smooth velocity model, perform time-frequency domain amplitude-phase inversion of the seismic envelope to obtain an updated smooth velocity model:

[0090]

[0091] Among them, V update represents the updated smooth velocity model; represents the phase adjustment amount, and the calculation formula is as follows:

[0092]

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

[0094] S5. Use the improved simulated annealing method to optimize and solve the updated smooth velocity model, considering the probability parameter, 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 represents the optimized high-precision velocity model; SA represents the simulated annealing algorithm function; P represents the probability parameter; T init represents the initial temperature; T final represents the final temperature.

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

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

[0100] The energy function is as follows:

[0101]

[0102] Among them, 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: Randomly generate a new solution V current within the neighborhood of the current solution V new , and 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 a probability of exp(-ΔE / T), where T represents the current temperature;

[0106] S54. Cooling strategy: Reduce the temperature according to the predetermined cooling rate δ, that is, 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 the amplitude-phase information in the time-frequency domain 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 velocity model; β represents the final adjustment factor.

[0110] It should be noted that the content not elaborated in detail in the present invention is all prior art and is well known to those skilled in the art.

[0111] Therefore, by adopting the above-mentioned full waveform inversion method based on the joint amplitude-phase in the time-frequency domain, the present invention can make full use of the amplitude and phase information of seismic data, improve the accuracy and resolution of inversion, and achieve precise detection of the complex deep structure of coal fields.

[0112] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions of the present invention or make equivalent replacements, and these modifications or equivalent replacements do not enable 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 the combination of amplitude and phase in time and frequency domains, characterized in that: The following steps are involved: S1, input the initial velocity model, source wavelet and seismic data to establish the initial smooth velocity model; S2, preprocessing the seismic data to obtain preprocessed data; S3, using the preprocessed data to calculate the time-frequency domain amplitude adjustment factor, 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, performing time-frequency domain amplitude phase inversion of the seismic data envelope to obtain an updated smooth velocity model; S5. Optimizing and solving the updated smooth velocity model using an 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 the amplitude and phase information in the time-frequency domain to obtain the final high-precision velocity model.

2. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 1 is characterized in that: In step S1, the calculation formula of the initial smoothing velocity model is as follows: V init =V0+α(DE) (1); Among them, V init represents the initial smooth velocity model; V0 represents the initial velocity model; D represents the source wavelet; E represents the seismic data; α represents the initial adjustment factor.

3. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 2 is characterized in that: In step S2, preprocessing includes Wiener filter scale decomposition, Hilbert transform and short-time Fourier transform; Among them, the Wiener filter scale is decomposed as follows: E wiener =Wiener(E,σ 2 ) (2); Among them, E wiener Represents the data after decomposition by Wiener filter; σ 2 represents the noise power spectral density; The Hilbert transform is as follows: H=Hilbert(E Wiener ) (3); Where H represents the Hilbert transform result; The short-time Fourier transform is as follows: S=STFT(E Wiener ,ω,t) (4); Among them, S represents the short-time Fourier transform result; ω represents the window function; and t represents the time shift parameter.

4. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 3 is characterized in that: In step S3, the time-frequency domain amplitude adjustment factor is calculated using the preprocessed data. The specific operation is as follows: Among them, A reg Represents the amplitude adjustment factor in the time-frequency domain.

5. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 4 is characterized in that: In step S3, the amplitude-phase objective function in the time-frequency domain is constructed, and the specific operation is as follows: Where F represents the amplitude phase objective function in the time-frequency domain; represents the phase of the short-time Fourier transform; represents the phase of the Hilbert transform; i represents the imaginary unit.

6. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 5 is characterized in that: In step S4, based on the time-frequency domain amplitude phase objective function and the initial smooth velocity model, the time-frequency domain amplitude phase inversion of the seismic data envelope is performed to obtain an updated smooth velocity model. The specific operations are: Among them, V update represents the updated smoothed velocity model; Represents the phase adjustment amount, and the calculation formula is as follows: Among them, η represents the learning rate; Represents the gradient of the objective function F with respect to the phase.

7. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 6 is characterized in that: In step S5, the initial high-precision velocity model is as follows: In high =SA(V update ,P,T init ,T final ) (9); Among them, V high represents the optimized high-precision speed model; SA represents the simulated annealing algorithm function; P represents the probability parameter; T init represents the initial temperature; T final Indicates the final temperature.

8. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 7 is characterized in that: In step S5, the improved simulated annealing method includes the following steps: S51, parameter initialization, select an initial velocity model, set the initial temperature T init and the final temperature T final And define the energy function E(V); The energy function is as follows: Wherein, V represents the velocity model; 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; S52, iterative search, in the current solution V current A new solution V is randomly generated in the neighborhood of new , calculate the energy difference ΔE between the current solution and the new solution: ΔE=E(V new )-E(V current )(11); 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; S54, cooling strategy, reducing the temperature according to the predetermined cooling rate δ, when the termination condition is reached, outputting the initial high-precision velocity model V high .

9. The full waveform inversion method based on the combination of amplitude and phase in the time-frequency domain according to claim 8 is characterized in that: In step S6, the final high-precision velocity model is as follows: V final =V high +β(V update -V high )(12); Among them, V final represents the final high-precision velocity model; β represents the final adjustment factor.

Images