Seismic wavelet estimation method and system, electronic equipment and storage medium

Abstract

The invention provides a seismic wavelet estimation method and system, electronic equipment and a storage medium. The method comprises the following steps: inputting initial values of an order, an amplitude coefficient, a dominant frequency parameter and a phase angle parameter into a seismic wavelet model to generate a first seismic wavelet; establishing an objective function according to the first seismic wavelet, inputting the obtained amplitude coefficient, the main frequency parameter and the phase angle parameter which enable the value of the objective function to be minimum into a seismic wavelet model to generate a second seismic wavelet, and obtaining a current synthetic seismic record according to the second seismic wavelet; calculating a correlation coefficient between the time domain seismic record fragment and the current synthetic seismic record; changing the initial value of the order, and repeating the steps to obtain a plurality of correlation coefficients; inputting the numerical values of the four parameters corresponding to the maximum correlation coefficient into the seismic wavelet model to obtain a seismic wavelet estimation result; according to the method, a multi-parameter seismic wavelet model is adopted, the spatial dimension is larger, and the precision of an estimation result is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data analysis, and in particular relates to a seismic wavelet estimation method, system, electronic device and storage medium. Background Art

[0002] Currently, seismic wavelets are the basis for processes such as seismic data quality assessment, deconvolution, well-seismic calibration, seismic forward and inverse modeling, and hydrocarbon reservoir prediction. Seismic wavelet estimation is one of the important steps for qualitative and quantitative interpretation and analysis of seismic data. The process of model-based time-domain seismic wavelet estimation is usually based on a seismic wavelet model with an analytical expression. First, the control parameters of the seismic wavelet model are estimated from seismic data and well logging data, and then the seismic wavelet waveform is obtained from the seismic wavelet analytical formula. To obtain reliable seismic wavelet results using the model-based time-domain seismic wavelet estimation method, the key lies in the seismic wavelet model used. If the model space of the seismic wavelet is small, for example, the Ricker wavelet model commonly used in hydrocarbon seismic exploration, which has only two model control parameters: amplitude coefficient and dominant frequency parameter, and the model space dimension is 2, then it is impossible to obtain high-precision seismic wavelet results from the input seismic data and well logging data, which will further affect the result accuracy of subsequent well-seismic calibration, seismic inversion, hydrocarbon reservoir prediction, etc.

[0003] Therefore, there is an urgent need to develop a seismic wavelet estimation method, system, electronic device and storage medium to solve the above technical problems. Summary of the Invention

[0004] To solve the above technical problems, the present invention proposes a seismic wavelet estimation method, system, electronic device and storage medium to solve the above technical problems.

[0005] The first aspect of the present invention discloses a seismic wavelet estimation method, and the method includes:

[0006] Step S1: Calculate the corresponding time-domain reflection coefficient using well logging velocity data and well logging density data;

[0007] Step S2: Extract a time-domain seismic record corresponding to the well logging position coordinates from the three-dimensional work area time-domain seismic data, and intercept a time-domain seismic record segment whose time range corresponds to the time-domain reflection coefficient from the time-domain seismic record;

[0008] Step S3: Input the initial values of the set order, amplitude coefficient, dominant frequency parameter and phase angle parameter into a pre-established seismic wavelet model to generate a first seismic wavelet;

[0009] Step S4: Establish an objective function based on the first seismic wavelet, and obtain the amplitude coefficient, dominant frequency parameter and phase angle parameter when the value of the objective function reaches the minimum.

[0010] Step S5: Input the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum into the seismic wavelet model to generate a second seismic wavelet, and produce the current synthetic seismic record based on the second seismic wavelet;

[0011] Step S6: Calculate the correlation coefficient between the time-domain seismic record segment and the current synthetic seismic record;

[0012] Step S7: Change the initial value of the order, and repeat the above Steps S3 - S6 until multiple correlation coefficients are obtained;

[0013] Step S8: Input the values of the order, amplitude coefficient, main frequency parameter, and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate a third seismic wavelet and use it as the seismic wavelet estimation result.

[0014] According to the method of the first aspect of the present invention, in the said Step S3, the seismic wavelet model is a four-parameter seismic wavelet model, denoted as w, and the expression is:

[0015] w(t|u,a,f,θ)=w(t|u,a,f)cosθ - w(t|u,a,f) h sinθ (1)

[0016] where t is time, t ∈ (-∞, +∞), with the unit of millisecond; u is the order, u ∈ [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dimensionless; a is the amplitude coefficient, a ∈ (0, +∞), dimensionless; f is the main frequency parameter, with the unit of hertz; dt is the time sampling interval of the time-domain seismic data, with the unit of millisecond; θ is the phase angle parameter, θ ∈ (-180°, 180°], with the unit of degree; w(t|u,a,f) h is the Hilbert transform of w(t|u,a,f);

[0017] The specific analytical formula of w(t|u,a,f) is:

[0018]

[0019] where x represents the independent variable, for Equation (2), exp() is the exponential function with the natural constant e as the base; P0(x) = 1, P1(x) = 2x, and for the remaining u orders, P u is calculated recursively according to the following formula

[0020] P u (x) = xP u-1 (x) - (u - 1)P u-2 (x) (3).

[0021] According to the method of the first aspect of the present invention, step S3 specifically includes:

[0022] Step S31: Set the seismic wavelet model;

[0023] Step S32: Set the initial value of the order u to u0, and let u = u0 = 1;

[0024] Step S33: Set the time axis vector of the seismic wavelet as t = [t1, …, t l , …, t L , L is the number of elements in the time axis vector, and L is a positive odd number. l represents the element serial number in the time axis vector;

[0025] Step S34: Set the initial values of the amplitude coefficient a, the main frequency parameter f, and the phase angle parameter θ to a0, f0, and θ0 respectively;

[0026] Step S35: Let a = a0, f = f0, θ = θ0, and substitute them into the seismic wavelet model to generate the first seismic wavelet.

[0027] According to the method of the first aspect of the present invention, in step S3, the initial value a0 of the amplitude coefficient is set by using the ratio of the maximum value of the amplitude envelope of the time-domain reflection coefficient r to the time-domain seismic record segment y, which is expressed as:

[0028]

[0029] where r env and y env are respectively the amplitude envelopes of the time-domain reflection coefficient r and the time-domain seismic record segment y.

[0030] According to the method of the first aspect of the present invention, in step S3, the initial value f0 of the main frequency parameter is determined by using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y, which is expressed as:

[0031] f0 = max(|fft(y)|) (5)

[0032] where fft() represents the Fourier transform operation.

[0033] According to the method of the first aspect of the present invention, in step S3, the initial value θ0 of the phase angle parameter is estimated from the time-domain seismic record segment y by using the kurtosis phase estimation method.

[0034] According to the method of the first aspect of the present invention, in step S4, the objective function is represented as F, and the expression is:

[0035]

[0036] wherein, * represents a convolution operation, and σ is the standard deviation of y - r * w(t|u, a, f, θ), which is the amplitude coefficient to be solved, dominant frequency parameter and phase angle parameter to form a parameter vector μ c = [a f θ] T ,

[0037] The amplitude coefficient, dominant frequency parameter, and phase angle parameter that minimize the value of the objective function F are obtained by using an unconstrained multivariable optimization algorithm. dominant frequency parameter and phase angle parameter That is,

[0038] The second aspect of the present invention discloses a seismic wavelet estimation system, which includes:

[0039] A first processing module configured to calculate the corresponding time - domain reflection coefficient by using logging velocity data and logging density data;

[0040] A second processing module configured to extract a time - domain seismic record corresponding to the logging position coordinates from the three - dimensional work area time - domain seismic data, and intercept a time - domain seismic record segment whose time range corresponds to the time - domain reflection coefficient from the time - domain seismic record;

[0041] A third processing module configured to input the initial values of the set order, amplitude coefficient, dominant frequency parameter, and phase angle parameter into a pre - established seismic wavelet model to generate a first seismic wavelet;

[0042] A fourth processing module configured to establish an objective function based on the first seismic wavelet and obtain the amplitude coefficient, dominant frequency parameter, and phase angle parameter that minimize the value of the objective function;

[0043] A fifth processing module configured to input the amplitude coefficient, dominant frequency parameter, and phase angle parameter that minimize the value of the objective function into the seismic wavelet model to generate a second seismic wavelet, and produce a current synthetic seismic record based on the second seismic wavelet;

[0044] A sixth processing module configured to calculate the correlation coefficient between the time - domain seismic record segment and the current synthetic seismic record;

[0045] A seventh processing module configured to change the initial value of the order and repeat the processing of the above - mentioned third processing module to the sixth processing module until multiple correlation coefficients are obtained;

[0046] The eighth processing module is configured to input the values of the order corresponding to the maximum correlation coefficient, the amplitude coefficient, the main frequency parameter, and the phase angle parameter into the seismic wavelet model to generate a third seismic wavelet and use it as the seismic wavelet estimation result.

[0047] According to the system of the second aspect of the present invention, in the third processing module, the seismic wavelet model is a four-parameter seismic wavelet model, expressed as:

[0048] w(t|u,a,f,θ) = w(t|u,a,f)cosθ - w(t|u,a,f) h sinθ (1)

[0049] where t is time, t ∈ (-∞, +∞), with the unit of millisecond; u is the order, u ∈ [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dimensionless; a is the amplitude coefficient, a ∈ (0, +∞), dimensionless; f is the main frequency parameter, with the unit of hertz; dt is the time sampling interval of the time-domain seismic data, with the unit of millisecond; θ is the phase angle parameter, θ ∈ (-180°, 180°], with the unit of degree; w(t|u,a,f) h is the Hilbert transform of w(t|u,a,f);

[0050] The specific analytical formula of w(t|u,a,f) is:

[0051]

[0052] where x represents the independent variable. For equation (2), exp() is the exponential function with the natural constant e as the base; P0(x) = 1, P1(x) = 2x, and for the remaining u-th order, P u is calculated recursively according to the following formula

[0053] P u (x) = xP u-1 (x) - (u - 1)P u-2 (x) (3).

[0054] According to the system of the second aspect of the present invention, the third processing module is specifically configured to:

[0055] Set the seismic wavelet model;

[0056] Set the initial value of the order u to u0 and let u = u0 = 1;

[0057] Set the time axis vector of the seismic wavelet to t = [t1, …, t l , …, t L , L is the number of elements in the time axis vector, and L is a positive odd number;

[0058] Set the initial values of the amplitude coefficient a, the dominant frequency parameter f, and the phase angle parameter θ to a0, f0, and θ0 respectively;

[0059] Let a = a0, f = f0, θ = θ0, and substitute them into the seismic wavelet model to generate the first seismic wavelet.

[0060] According to the system of the second aspect of the present invention, in the third processing module, the initial value a0 of the amplitude coefficient is set by using the ratio of the time-domain reflection coefficient r to the maximum value of the amplitude envelope of the time-domain seismic record segment y, expressed as:

[0061]

[0062] where r env and y env are respectively the amplitude envelopes of the time-domain reflection coefficient r and the time-domain seismic record segment y.

[0063] According to the system of the second aspect of the present invention, in the third processing module, the initial value f0 of the dominant frequency parameter is determined by using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y, expressed as:

[0064] f0 = max(|fft(y)|) (5)

[0065] where fft() represents the Fourier transform operation.

[0066] According to the system of the second aspect of the present invention, in the third processing module, the initial value θ0 of the phase angle parameter is estimated from the time-domain seismic record segment y by using the kurtosis phase estimation method.

[0067] According to the system of the second aspect of the present invention, in the fourth processing module, the objective function F is expressed as:

[0068]

[0069] where * represents the convolution operation, σ is the standard deviation of y - r * w(t|u, a, f, θ), is the parameter vector composed of the amplitude coefficient to be solved the dominant frequency parameter and the phase angle parameter constitute, μ c = [a f θ] T ,

[0070] Use the unconstrained multivariable optimization algorithm to obtain the amplitude coefficient when the value of the objective function F reaches the minimum the dominant frequency parameter and the phase angle parameter That is,

[0071] A third aspect of the present invention discloses an electronic device. The electronic device includes a memory and a processor. The memory stores a computer program. When the processor executes the computer program, the steps in a seismic wavelet estimation method according to any one of the first aspects of the present disclosure are implemented.

[0072] A fourth aspect of the present invention discloses a computer-readable storage medium. A computer program is stored on the computer-readable storage medium. When the computer program is executed by a processor, the steps in a seismic wavelet estimation method according to any one of the first aspects of the present disclosure are implemented.

[0073] In summary, the solution proposed by the present invention adopts a multi-parameter seismic wavelet model that expands the space of the existing Ricker wavelet model. Its model space has a larger dimension and the waveform is more flexible and variable. At the same time, a corresponding objective function is set based on the multi-parameter seismic wavelet model to estimate the control parameters retroactively, improving the accuracy of the estimation result in time-domain seismic wavelet estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for use in the description of the specific embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0075] Figure 1 FIG. is a flowchart of a seismic wavelet estimation method according to an embodiment of the present invention;

[0076] Figure 2 FIG. is a schematic diagram of the seismic wavelet model changing with different control parameters according to an embodiment of the present invention;

[0077] Figure 3 FIG. is the input data used in the first example according to an embodiment of the present invention; wherein, (a) is the well logging density data that has been converted to the time domain, (b) is the well logging longitudinal wave velocity data that has been converted to the time domain, (c) is the time domain reflection coefficient calculated based on the well logging density data and the well logging longitudinal wave velocity data shown in (a) and (b), and (d) is a time domain seismic record segment corresponding to the time domain reflection coefficient;

[0078] Figure 4Input data used in the second example according to an embodiment of the present invention; wherein, (a) is the well logging density data that has been transformed into the time domain, (b) is the well logging compressional wave velocity data that has been transformed into the time domain, (c) is the reflection coefficient in the time domain calculated based on the well logging density data and the well logging compressional wave velocity data shown in (a) and (b), and (d) is the time domain seismic record segment corresponding to the reflection coefficient in the time domain;

[0079] Figure 5 Schematic diagram for determining the initial value of the main frequency parameter in the first example according to an embodiment of the present invention;

[0080] Figure 6 Schematic diagram for determining the initial value of the phase angle parameter in the first example according to an embodiment of the present invention;

[0081] Figure 7 Schematic diagram for determining the initial value of the main frequency parameter in the second example according to an embodiment of the present invention;

[0082] Figure 8 Schematic diagram for determining the initial value of the phase angle parameter in the second example according to an embodiment of the present invention;

[0083] Figure 9 Schematic diagram of the (b) seismic wavelet estimation result obtained based on the (a) seismic wavelet model of the input data in the first example according to an embodiment of the present invention;

[0084] Figure 10 Schematic diagram of the (b) seismic wavelet estimation result obtained based on the (a) Ricker wavelet model of the input data in the first example according to an embodiment of the present invention;

[0085] Figure 11 Schematic diagram of the (b) seismic wavelet estimation result obtained based on the (a) seismic wavelet model of the input data in the second example according to an embodiment of the present invention;

[0086] Figure 12 Schematic diagram of the (b) seismic wavelet estimation result obtained based on the (a) Ricker wavelet model of the input data in the second example according to an embodiment of the present invention;

[0087] Figure 13 Structural diagram of a seismic wavelet estimation system according to an embodiment of the present invention;

[0088] Figure 14 Structural diagram of an electronic device according to an embodiment of the present invention. Detailed implementation manner

[0089] To make the objectives, technical solutions and advantages of the present invention more clear, the following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without making creative efforts shall fall within the protection scope of the present invention.

[0090] The first aspect of the present invention discloses a seismic wavelet estimation method. Figure 1 As shown in the flowchart of a seismic wavelet estimation method according to an embodiment of the present invention, Figure 1 the method includes:

[0091] Step S1: Calculate the corresponding reflection coefficient in the time domain using well logging velocity data and well logging density data;

[0092] Step S2: Extract a time-domain seismic record corresponding to the well logging position coordinates from the three-dimensional work area time-domain seismic data, and intercept a time-domain seismic record segment with a time range corresponding to the time-domain reflection coefficient from the time-domain seismic record;

[0093] Step S3: Input the initial values of the set order, amplitude coefficient, main frequency parameter, and phase angle parameter into a pre-established seismic wavelet model to generate a first seismic wavelet;

[0094] Step S4: Establish an objective function based on the first seismic wavelet, and obtain the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum;

[0095] Step S5: Input the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum into the seismic wavelet model to generate a second seismic wavelet, and produce the current synthetic seismic record according to the second seismic wavelet;

[0096] Step S6: Calculate the correlation coefficient between the time-domain seismic record segment and the current synthetic seismic record;

[0097] Step S7: Change the initial value of the order, and repeat the above steps S3 - S6 until multiple correlation coefficients are obtained;

[0098] Step S8: Input the values of the order, amplitude coefficient, main frequency parameter, and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate a third seismic wavelet and use it as the seismic wavelet estimation result.

[0099] In step S1, the corresponding reflection coefficient r in the time domain is calculated using well logging velocity data and well logging density data.

[0100] In step S2, extract the seismic record in the time domain corresponding to the well logging position coordinates from the 3D work area time domain seismic data, and further intercept a time domain seismic record segment y from this seismic record, whose time range corresponds to the time domain reflection coefficient.

[0101] In step S3, input the initial values of the set order, amplitude coefficient, main frequency parameter, and phase angle parameter into the pre-established seismic wavelet model to generate the first seismic wavelet.

[0102] In some embodiments, in step S3, the seismic wavelet model is a four-parameter seismic wavelet model, expressed as:

[0103] w(t|u,a,f,θ)=w(t|u,a,f)cosθ - w(t|u,a,f) h sinθ (1)

[0104] where t is time, t ∈ (-∞, +∞), with the unit of millisecond; u is the order, u ∈ [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dimensionless; a is the amplitude coefficient, a ∈ (0, +∞), dimensionless; f is the main frequency parameter, with the unit of hertz; dt is the time sampling interval of the time domain seismic data, with the unit of millisecond; θ is the phase angle parameter, θ ∈ (-180°, 180°], with the unit of degree; w(t|u,a,f) h is the Hilbert transform of w(t|u,a,f);

[0105] The specific analytical formula of w(t|u,a,f) is:

[0106]

[0107] where x represents the independent variable. For formula (2), exp() is the exponential function with the natural constant e as the base; P0(x) = 1, P1(x) = 2x, and for the remaining u orders, P u is calculated recursively according to the following formula

[0108] P u (x) = xP u-1 (x) - (u - 1)P u-2 (x) (3).

[0109] Specifically, Figure 2It is a schematic diagram of the seismic wavelet model of the embodiment of the present invention changing with different control parameters. Among them, in the first column of figures from the left, the amplitude coefficient a of the seismic wavelet model is 1, the main frequency parameter f is 40 Hz, the phase angle parameter θ is 0°, and the orders are u = 1, u = 2, and u = 10 respectively; in the second column of figures from the left, the order u of the seismic wavelet model is 2, the main frequency parameter f is 40 Hz, the phase angle parameter θ is 0°, and the amplitude coefficients are a = 2, a = 8, and a = 16 respectively; in the third column of figures from the left, the amplitude coefficient a of the seismic wavelet model is 1, the order u is 2, the phase angle parameter θ is 0°, and the main frequency parameters are f = 20 Hz, f = 40 Hz, and f = 60 Hz respectively; in the fourth column of figures from the left, the amplitude coefficient a of the seismic wavelet model is 1, the order u is 2, the main frequency parameter f is 40 Hz, and the phase angle parameters are θ = -170°, θ = 90°, and θ = 180° respectively.

[0110] Figure 2 It shows the waveform changes of the seismic wavelet model used in this embodiment with four model control parameters. Among them, when the order u = 2 and the phase angle parameter is θ = 0°, it is the Ricker wavelet. It can be seen that the seismic wavelet model used in this embodiment has a larger model space and more flexible waveform changes.

[0111] In some embodiments, step S3 specifically includes:

[0112] Step S31: Set the seismic wavelet model;

[0113] Step S32: Set the initial value of the order u as u0, and let u = u0 = 1;

[0114] Step S33: Set the time axis vector of the seismic wavelet as t = [t1,…,t l ,…,t L , L is the number of elements in the time axis vector, and L is a positive odd number. l represents the element serial number in the time axis vector;

[0115] Step S34: Set the initial values of the amplitude coefficient a, the main frequency parameter f, and the phase angle parameter θ as a0, f0, and θ0 respectively;

[0116] Step S35: Let a = a0, f = f0, θ = θ0, and substitute them into the seismic wavelet model to generate the first seismic wavelet.

[0117] In some embodiments, the initial value a0 of the amplitude coefficient is set by using the ratio of the time domain reflection coefficient r to the maximum value of the amplitude envelope of the time domain seismic record segment y, which is expressed as:

[0118]

[0119] where r env and y env are the amplitude envelopes of the time-domain reflection coefficient r and the time-domain seismic record segment y, respectively.

[0120] In some embodiments, the initial value f0 of the dominant frequency parameter is determined using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y, expressed as:

[0121] f0 = max(|fft(y)|) (5)

[0122] where fft() represents the Fourier transform operation.

[0123] In some embodiments, the initial value θ0 of the phase angle parameter is estimated from the time-domain seismic record segment y using the kurtosis phase estimation method.

[0124] Step S4, establish an objective function based on the first seismic wavelet, and obtain the amplitude coefficient, dominant frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum.

[0125] In some embodiments, the objective function F is expressed as:

[0126]

[0127] where * represents the convolution operation, σ is the standard deviation of y - r * w(t|u,a,f,θ), is the parameter vector composed of the amplitude coefficient to be solved dominant frequency parameter and phase angle parameter , μ c = [a f θ] T ,

[0128] Use an unconstrained multivariable optimization algorithm to obtain the amplitude coefficient when the value of the objective function F reaches the minimum dominant frequency parameter and phase angle parameter That is,

[0129] Specifically, taking two specific examples for illustration, the method of the embodiments of the present invention is used to obtain the seismic wavelet estimation results from the input data of the first example and the second example, including the following main steps:

[0130] (1) Calculate the corresponding time-domain reflection coefficient r using the logging velocity data and logging density data.

[0131] (2) Extract the time-domain seismic record corresponding to the well logging position coordinates from the 3D work area time-domain seismic data, and further extract the time-domain seismic record segment y whose time range corresponds to the time-domain reflection coefficient from this seismic record.

[0132] Figure 3 is the input data used in the first example; where (a) is the well logging density data that has been converted to the time domain, (b) is the well logging P-wave velocity data that has been converted to the time domain, (c) is the time-domain reflection coefficient calculated based on the well logging density data and well logging P-wave velocity data shown in (a) and (b), and (d) is the time-domain seismic record segment corresponding to the time-domain reflection coefficient.

[0133] Figure 4 is the input data used in the second example; where (a) is the well logging density data that has been converted to the time domain, (b) is the well logging P-wave velocity data that has been converted to the time domain, (c) is the time-domain reflection coefficient calculated based on the well logging density data and well logging P-wave velocity data shown in (a) and (b), and (d) is the time-domain seismic record segment corresponding to the time-domain reflection coefficient.

[0134] (3) Set the seismic wavelet model as

[0135] w(t|u,a,f,θ)=w(t|u,a,f)cosθ - w(t|u,a,f) h sinθ (1)

[0136] In the above formula, t is time, t∈(-∞,+∞), and the unit is millisecond (ms); u is the order, u∈[1,2,3,4,5,6,7,8,9,10], dimensionless; a is the amplitude coefficient, a∈(0,+∞), dimensionless; f is the main frequency parameter, the unit is Hertz (Hz), dt is the time sampling interval of the time-domain seismic data, and the unit is millisecond (ms); θ is the phase angle parameter, θ∈(-180°, 180°], and the unit is degree (°); w(t|u,a,f) h is the Hilbert transform of w(t|u,a,f), and the specific analytical formula of w(t|u,a,f) is:

[0137]

[0138] where, exp() is the exponential function with the natural constant e as the base; P0(x) = 1, P1(x) = x, and for the remaining u orders, P u is calculated recursively according to the following formula

[0139] P u (x)=xP u-1 (x)-(n - 1)Pu-2 (x)(3).

[0140] (4) Set the initial value of the order to u0 = 1, and let u = u0.

[0141] (5) Set the time axis vector of the seismic wavelet to t = [t1, …, t l , …, t L , where L is the number of elements in the time axis vector, and L is a positive odd number.

[0142] (6) Set the initial values of the amplitude coefficient, dominant frequency parameter, and phase angle parameter to a0, f0, θ0. For the initial value a0 of the amplitude coefficient, it can be set using the ratio of the amplitude envelope maximum value of the time-domain reflection coefficient r to the time-domain seismic record segment y, i.e.,

[0143]

[0144] where r env and y env are the amplitude envelopes of the time-domain reflection coefficient r and the time-domain seismic record segment y respectively; for the initial value f0 of the dominant frequency parameter, it can be determined using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y, i.e.,

[0145] f0 = max(|fft(y)|) (5)

[0146] where fft() represents the Fourier transform operation.

[0147] For the initial value θ0 of the phase angle parameter, it can be estimated from the time-domain seismic record segment y using the kurtosis phase estimation method: within a certain angular range, for each angular value, first perform phase rotation on the time-domain seismic record segment y to obtain y θ , and then calculate the corresponding kurtosis kurt θ :

[0148]

[0149] where is the mean of y θ , N is the number of elements in the time-domain seismic record segment y, and then from these kurtosis results, select the angle corresponding to the maximum kurtosis as the initial value θ0 of the phase angle parameter.

[0150] Figure 5 Figure 53 is a schematic diagram of determining the initial value f0 of the dominant frequency parameter using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y in the first example, and the determined initial value of the dominant frequency parameter f0 = 24.1 Hz.

[0151] Figure 6 In the first example, it is a schematic diagram for determining the initial value θ0 of the phase angle parameter from the time-domain seismic record segment y using the kurtosis phase estimation method. The determined initial value of the phase angle parameter θ0 = -1°, and the angular range used by the kurtosis phase estimation method is [-90°, 90°].

[0152] Figure 7 In the second example, it is a schematic diagram for determining the initial value f0 of the dominant frequency parameter using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y. The determined initial value of the dominant frequency parameter f0 = 30.94 Hz.

[0153] Figure 8 In the second example, it is a schematic diagram for determining the initial value θ0 of the phase angle parameter from the time-domain seismic record segment y using the kurtosis phase estimation method. The determined initial value of the phase angle parameter θ0 = 35°, the angular range used by the kurtosis phase estimation method is [-90°, 90°], with an interval of 1°.

[0154] (7) Let a = a0, f = f0, θ = θ0, and generate the seismic wavelet w(t|u, a, f, θ).

[0155] (8) Set the objective function F

[0156]

[0157] where * represents the convolution operation, σ is the standard deviation of y - r * w(t|u, a, f, θ), is the amplitude coefficient to be solved dominant frequency parameter and phase angle parameter constitute the parameter vector μ c = [a f θ] T , Use the Nelder-Mead simplex search method (Lagarias et al., 1998) or the quasi-Newton method to obtain the amplitude coefficient dominant frequency parameter and phase angle parameter when the value of the objective function F reaches the minimum, that is,

[0158] (9) Let Then generate the seismic wavelet w(t|u, a, f, θ) and produce the synthetic seismic record

[0159] (10) Calculate the correlation coefficient k between the time-domain seismic record segment y and the current synthetic seismic record u .

[0160] Let \(u = u + 1\), return to step (6), and repeat steps (6)-(10) until the process of \(u = 10\) is completed.

[0161] (12) From the above correlation coefficient results \([k1,k2,k3,k4,k5,k6,k7,k8,k9,k 10 , select the maximum correlation coefficient \(k opt =\max[k1,k2,k3,k4,k5,k6,k7,k8,k9,k 10 , and the order \(u opt \) corresponding to \(k opt , the amplitude coefficient \(a opt , the main frequency parameter \(f opt \), and the phase angle parameter \(\theta opt , and generate the seismic wavelet \(w(t|u opt ,a opt ,f opt ,\theta opt ) as the final seismic wavelet estimation result.

[0162] Table 1 and Table 2 respectively give the estimation results of the seismic wavelet model control parameters obtained by using the method of the present invention from the input data of the first example and the second example.

[0163] Figure 9 Schematic diagram of (b) the seismic wavelet estimation result obtained based on the (a) seismic wavelet model of the first example input data according to the embodiment of the present invention; Figure 10 Schematic diagram of (b) the seismic wavelet estimation result obtained based on the (a) Ricker wavelet model of the first example input data according to the embodiment of the present invention.

[0164] For the first example, \(k opt = 90.56\%\), \(u opt = 1\), \(a opt = 1.13×10 6 , \(f opt = 18.91Hz\), \(\theta opt = 17.17^{\circ}\), and the corresponding seismic wavelet \(w(t|u opt ,a opt ,f opt ,\theta opt ), as shown in (a) and (b) of Figure 9 , indicates the correlation coefficient between the synthetic seismic record (the dashed line indicated in (b) of Figure 9 ) made from the seismic wavelet estimation result (shown in (a) of Figure 9 ) and the time-domain seismic record segment \(y( Figure 9 the solid line indicated in (b)) is 90.56\%; refer to Figure 10As shown in (a) and (b), the correlation coefficient between the synthetic seismic record (the dashed line shown in (b) of Figure 10 ) made using the seismic wavelet estimation result obtained based on the conventional Ricker wavelet model (shown in (a) of Figure 10 ) and the time-domain seismic record segment y (the solid line shown in (b) of Figure 10 ) is only 84.18%.

[0165] Figure 11 Schematic diagram of the (b) seismic wavelet estimation result obtained based on the (a) seismic wavelet model of the present invention for the second instance input data; Figure 12 Schematic diagram of the (b) seismic wavelet estimation result obtained based on the (a) Ricker wavelet model for the second instance input data of the embodiments of the present invention.

[0166] For the second instance, k opt = 67.34%, u opt = 4, a opt = 3.11×10 4 , f opt = 31.37Hz, θ opt = 46.51°, and the corresponding seismic wavelet w(t|u opt , a opt , f opt , θ opt ), as shown in (a) and (b) of Figure 11 , the correlation coefficient between the synthetic seismic record (the dashed line shown in (b) of Figure 11 ) made from the seismic wavelet estimation result (shown in (a) of Figure 11 ) and the time-domain seismic record segment y (the solid line shown in (b) of Figure 11 ) is 67.34%; referring to (a) and (b) of Figure 12 , the correlation coefficient between the synthetic seismic record (the dashed line shown in (b) of Figure 12 ) made using the seismic wavelet estimation result obtained based on the conventional Ricker wavelet model (shown in (a) of Figure 12 ) and the time-domain seismic record segment y (the solid line shown in (b) of Figure 12 ) is only 43.54%. In the above figures, the time sampling intervals of the seismic wavelet, logging density, logging longitudinal wave velocity, reflection coefficient, seismic record segment, and synthetic seismic record are all dt = 2ms. It can be seen that compared with the existing method, the seismic wavelet estimation method using the embodiments of the present invention obtains a seismic wavelet estimation result with higher accuracy.

[0167] Table 1 Seismic wavelet estimation results of the first instance of the present invention

[0168]

[0169]

[0170] Table 2 Seismic wavelet estimation results of the second example of the present invention

[0171]

[0172] In summary, the solution proposed by the present invention adopts a multi-parameter seismic wavelet model that expands the space of the existing Ricker wavelet model. Its model space has a larger dimension and more flexible waveforms. At the same time, a corresponding objective function is set based on the multi-parameter seismic wavelet model to estimate the control parameters retroactively, improving the accuracy of the estimation results in time-domain seismic wavelet estimation.

[0173] The second aspect of the present invention discloses a seismic wavelet estimation system. Figure 13 Structural diagram of a seismic wavelet estimation system according to an embodiment of the present invention; as Figure 13 shown, the system 100 includes:

[0174] The first processing module 101 is configured to calculate the corresponding time-domain reflection coefficient by using well logging velocity data and well logging density data;

[0175] The second processing module 102 is configured to extract a time-domain seismic record corresponding to the well logging position coordinates from the three-dimensional work area time-domain seismic data, and intercept a time-domain seismic record segment with a time range corresponding to the time-domain reflection coefficient from the time-domain seismic record;

[0176] The third processing module 103 is configured to input the initial values of the set order, amplitude coefficient, main frequency parameter, and phase angle parameter into a pre-established seismic wavelet model to generate a first seismic wavelet;

[0177] The fourth processing module 104 is configured to establish an objective function based on the first seismic wavelet and obtain the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum;

[0178] The fifth processing module 105 is configured to input the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum into the seismic wavelet model to generate a second seismic wavelet, and produce the current synthetic seismic record according to the second seismic wavelet;

[0179] The sixth processing module 106 is configured to calculate the correlation coefficient between the time-domain seismic record segment and the current synthetic seismic record;

[0180] The seventh processing module 107 is configured to change the initial value of the order, and repeat the processing of the above third processing module to the sixth processing module until a plurality of correlation coefficients are obtained;

[0181] The eighth processing module 108 is configured to input the values of the order, amplitude coefficient, main frequency parameter, and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate a third seismic wavelet and use it as the seismic wavelet estimation result.

[0182] According to the system of the second aspect of the present invention, in the third processing module 103, the seismic wavelet model is a four-parameter seismic wavelet model, which is expressed as:

[0183] w(t|u,a,f,θ)=w(t|u,a,f)cosθ - w(t|u,a,f) h sinθ (1)

[0184] Where t is time, t ∈ (-∞, +∞), with the unit of millisecond; u is the order, u ∈ [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dimensionless; a is the amplitude coefficient, a ∈ (0, +∞), dimensionless; f is the main frequency parameter, with the unit of hertz; dt is the time sampling interval of the seismic data in the time domain, with the unit of millisecond; θ is the phase angle parameter, θ ∈ (-180°, 180°], with the unit of degree; w(t|u,a,f) h is the Hilbert transform of w(t|u,a,f);

[0185] The specific analytical formula of w(t|u,a,f) is:

[0186]

[0187] Where x represents the independent variable. For Equation (2), exp() is the exponential function with the natural constant e as the base; P0(x) = 1, P1(x) = 2x, and for the remaining u orders, P u is recursively calculated according to the following formula

[0188] P u (x) = xP u-1 (x) - (u - 1)P u-2 (x) (3).

[0189] According to the system of the second aspect of the present invention, the third processing module 103 is specifically configured to:

[0190] Set the seismic wavelet model;

[0191] Set the initial value of the order u to u0, and let u = u0 = 1;

[0192] Set the time axis vector of the seismic wavelet as \(t = [t_1, \ldots, t l , \ldots, t L \), where \(L\) is the number of elements in the time axis vector and \(L\) is a positive odd number;

[0193] Set the initial values of the amplitude coefficient \(a\), the dominant frequency parameter \(f\), and the phase angle parameter \(\theta\) as \(a_0\), \(f_0\), and \(\theta_0\) respectively;

[0194] Let \(a = a_0\), \(f = f_0\), \(\theta = \theta_0\), and substitute them into the seismic wavelet model to generate the first seismic wavelet.

[0195] According to the system of the second aspect of the present invention, in the third processing module 103, set the initial value \(a_0\) of the amplitude coefficient by using the ratio of the time-domain reflection coefficient \(r\) to the maximum value of the amplitude envelope of the time-domain seismic record segment \(y\), which is expressed as:

[0196]

[0197] where \(r env and \(y env are the amplitude envelopes of the time-domain reflection coefficient \(r\) and the time-domain seismic record segment \(y\) respectively.

[0198] According to the system of the second aspect of the present invention, in the third processing module 103, determine the initial value \(f_0\) of the dominant frequency parameter by using the peak frequency of the amplitude spectrum of the time-domain seismic record segment \(y\), which is expressed as:

[0199] \(f_0=\max(|fft(y)|)\) (5)

[0200] where \(fft()\) represents the Fourier transform operation.

[0201] According to the system of the second aspect of the present invention, in the third processing module 103, estimate the initial value \(\theta_0\) of the phase angle parameter from the time-domain seismic record segment \(y\) by using the kurtosis phase estimation method.

[0202] According to the system of the second aspect of the present invention, in the fourth processing module 104, the objective function \(F\) is expressed as:

[0203]

[0204] where \(*\) represents the convolution operation, \(\sigma\) is the standard deviation of \(y - r*w(t|u,a,f,\theta)\), \(u\) is the parameter vector composed of the amplitude coefficient dominant frequency parameter and phase angle parameter to be solved, \(\mu c =[af\theta] T ,

[0205] The amplitude coefficient when the value of the objective function F reaches the minimum is obtained by using an unconstrained multivariable optimization algorithm. Main frequency parameter and phase angle parameter That is,

[0206] The third aspect of the present invention discloses an electronic device. The electronic device includes a memory and a processor. The memory stores a computer program. When the processor executes the computer program, the steps in a seismic wavelet estimation method according to any one of the first aspects disclosed in the present invention are implemented.

[0207] Figure 14 FIG. is a structural diagram of an electronic device according to an embodiment of the present invention. As Figure 14 shown, the electronic device includes a processor, a memory, a communication interface, a display screen, and an input device connected through a system bus. Among them, the processor of the electronic device is used to provide computing and control capabilities. The memory of the electronic device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a computer program. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The communication interface of the electronic device is used to communicate with an external terminal in a wired or wireless manner. The wireless manner can be implemented through WIFI, a carrier network, near field communication (NFC), or other technologies. The display screen of the electronic device can be a liquid crystal display screen or an electronic ink display screen. The input device of the electronic device can be a touch layer covered on the display screen, or a button, a trackball, or a touchpad provided on the housing of the electronic device, or an external keyboard, a touchpad, or a mouse, etc.

[0208] Those skilled in the art can understand that Figure 14 the structure shown in is only a structural diagram of a part related to the technical solution of the present disclosure, and does not constitute a limitation on the electronic device to which the solution of the present application is applied. The specific electronic device may include more or fewer components than those shown in the figure, or combine certain components, or have a different component layout.

[0209] The fourth aspect of the present invention discloses a computer-readable storage medium. A computer program is stored on the computer-readable storage medium. When the computer program is executed by a processor, the steps in a seismic wavelet estimation method according to any one of the first aspects disclosed in the present invention are implemented.

[0210] Please note that the technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as within the scope described in this specification. The above embodiments only express several implementation manners of the present application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several modifications and improvements can still be made, and these all belong to the protection scope of the present application. Therefore, the protection scope of the patent of the present application shall be subject to the appended claims.

[0211] The above are the preferred implementation manners of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can still be made, and these improvements and refinements should also be regarded as within the protection scope of the present invention.

Claims

1. An earthquake wavelet estimation method, characterized in that, The method includes: Step S1: Calculate the corresponding time-domain reflection coefficient using well logging velocity data and well logging density data; Step S2: Extract a time-domain seismic record corresponding to the well logging position coordinates from the three-dimensional work area time-domain seismic data, and intercept a time-domain seismic record segment with a time range corresponding to the time-domain reflection coefficient from the time-domain seismic record; Step S3: Input the initial values of the set order, amplitude coefficient, main frequency parameter, and phase angle parameter into the pre-established seismic wavelet model to generate a first seismic wavelet; Step S4: Establish an objective function based on the first seismic wavelet, and obtain the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum; Step S5: Input the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum into the seismic wavelet model to generate a second seismic wavelet, and produce the current synthetic seismic record according to the second seismic wavelet; Step S6: Calculate the correlation coefficient between the time-domain seismic record segment and the current synthetic seismic record; Step S7: Change the initial value of the order, and repeat the above steps S3 - S6 until multiple correlation coefficients are obtained; Step S8: Input the values of the order, amplitude coefficient, main frequency parameter, and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate a third seismic wavelet and use it as the seismic wavelet estimation result.

2. The method for estimating a seismic wavelet according to claim 1, characterized in that, In the said step S3, the seismic wavelet model is a four-parameter seismic wavelet model, denoted as w, and the expression is as follows: w(t|u,a,f,θ) = w(t|u,a,f)cosθ - w(t|u,a,f) h sinθ (1) where t is time, t ∈ (-∞, +∞), with the unit of millisecond; u is the order, u ∈ [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dimensionless; a is the amplitude coefficient, a ∈ (0, +∞), dimensionless; f is the main frequency parameter, with the unit of hertz; dt is the time sampling interval of the time-domain seismic data, with the unit of millisecond; θ is the phase angle parameter, θ ∈ (-180°, 180°], with the unit of degree; w(t|u,a,f) h is the Hilbert transform of w(t|u,a,f); The specific analytical formula of w(t|u,a,f) is: Among them, x represents the independent variable; exp() is the exponential function with the natural constant e as the base; P0(x) = 1, P1(x) = 2x, and for the remaining u-th order, P u It is calculated recursively according to the following formula: P u P(x)=x u-1 P(x)-(u - 1) u-2 P(x) (3).

3. A method for estimating seismic wavelets according to claim 2, characterized in that, The said step S3 specifically includes: Step S31: Set the seismic wavelet model; Step S32: Set the initial value of the order u to u0, and let u = u0 = 1; Step S33: Set the time axis vector of the seismic wavelet as \(t = [t_1, \ldots, t l , \ldots, t L \), where \(L\) is the number of elements in the time axis vector, and \(L\) is a positive odd number, and \(l\) represents the element serial number in the time axis vector; Step S34: Set the initial values of the amplitude coefficient a, main frequency parameter f, and phase angle parameter θ to a0, f0, and θ0 respectively; Step S35: Let a = a0, f = f0, θ = θ0, and substitute them into the seismic wavelet model to generate a first seismic wavelet.

4. A method for estimating seismic wavelets according to claim 3, characterized in that, In the said step S3, the initial value a0 of the amplitude coefficient is set using the ratio of the time-domain reflection coefficient r to the maximum value of the amplitude envelope of the time-domain seismic record segment, expressed as: where r env and y env are the amplitude envelopes of the time-domain reflection coefficient r and the time-domain seismic record segment y, respectively.

5. A method for estimating seismic wavelets according to claim 3, characterized in that In the said step S3, the initial value f0 of the main frequency parameter is determined using the peak frequency of the amplitude spectrum of the time-domain seismic record segment, expressed as: f0 = max(|fft(y)|) (5) Where fft() represents the Fourier transform operation.

6. A method for estimating seismic wavelets according to claim 3, characterized in that, In the said step S3, the initial value θ0 of the phase angle parameter is estimated from the time-domain seismic record segment y using the kurtosis phase estimation method.

7. A method for estimating seismic wavelets according to claim 3, characterized in that, In the said step S4, the objective function is denoted as F, and the expression is: where * represents the convolution operation, and σ is the standard deviation of y - r * w(t|u, a, f, θ). is the amplitude coefficient to be solved main frequency parameter and phase angle parameter constitute the parameter vector μ c = [a f θ] T , The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariable optimization algorithm Main frequency parameter and phase angle parameter That is 8. An earthquake wavelet estimation system, characterized in that The system includes: A first processing module, configured to calculate the corresponding time-domain reflection coefficient using well logging velocity data and well logging density data; A second processing module, configured to extract a time-domain seismic record corresponding to the well logging position coordinates from the three-dimensional work area time-domain seismic data, and intercept a time-domain seismic record segment with a time range corresponding to the time-domain reflection coefficient from the time-domain seismic record; A third processing module, configured to input initial values of a set order, amplitude coefficient, main frequency parameter, and phase angle parameter into a pre-established seismic wavelet model to generate a first seismic wavelet; A fourth processing module, configured to establish an objective function based on the first seismic wavelet and obtain the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum; A fifth processing module, configured to input the amplitude coefficient, main frequency parameter, and phase angle parameter when the value of the objective function reaches the minimum into the seismic wavelet model to generate a second seismic wavelet, and produce a current synthetic seismic record according to the second seismic wavelet; A sixth processing module, configured to calculate the correlation coefficient between the time-domain seismic record segment and the current synthetic seismic record; A seventh processing module, configured to change the initial value of the order, and repeat the processing of the above third processing module to sixth processing module until multiple correlation coefficients are obtained; An eighth processing module, configured to input the values of the order, amplitude coefficient, main frequency parameter, and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate a third seismic wavelet and use it as the seismic wavelet estimation result.

9. An electronic device, characterized in that, The electronic device includes a memory and a processor. When the processor executes the computer program stored in the memory, the steps in a seismic wavelet estimation method according to any one of claims 1 to 7 are implemented.

10. A computer-readable storage medium, characterized in that, A computer program is stored on the computer-readable storage medium. When the computer program is executed by a processor, the steps in a seismic wavelet estimation method according to any one of claims 1 to 7 are implemented.

Images