A seismic wavelet estimation method, system, electronic device and storage medium

Abstract

The application provides a seismic wavelet estimation method, system, electronic equipment and storage medium. The method comprises the following steps: inputting initial values of an order, an amplitude coefficient, a main frequency parameter and a phase angle parameter into a seismic wavelet model to generate a first seismic wavelet; establishing a target function according to the first seismic wavelet, inputting the amplitude coefficient, the main frequency parameter and the phase angle parameter obtained when the value of the target function reaches the minimum into the 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 a time domain seismic record segment and the current synthetic seismic record; changing the initial value of the order, and repeating the above steps to obtain a plurality of correlation coefficients; and inputting the values of the four parameters corresponding to the maximum correlation coefficient into the seismic wavelet model to obtain a seismic wavelet estimation result. The method adopts a multi-parameter seismic wavelet model, the spatial dimension is larger, and the accuracy of the estimation result is improved.

Description

Technical Field

[0001] This invention belongs to the field of earthquake data analysis technology, and particularly relates to a seismic wavelet estimation method, system, electronic device and storage medium. Background Technology

[0002] Currently, seismic wavelets are fundamental to seismic data quality assessment, deconvolution, well-seismic calibration, seismic forward and inversion retrieval, and oil and gas reservoir prediction. Seismic wavelet estimation is a crucial step in the qualitative and quantitative interpretation and analysis of seismic data. The model-based time-domain seismic wavelet estimation process typically uses a seismic wavelet model with analytical expressions as its basis. First, the control parameters of the seismic wavelet model are estimated from seismic and well logging data. Then, the seismic wavelet waveform is obtained from the analytical expression. The key to obtaining reliable seismic wavelet results using model-based time-domain seismic wavelet estimation methods 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 oil and gas seismic exploration has only two model control parameters: amplitude coefficient and dominant frequency parameter, and a model space dimension of 2—then it is impossible to obtain high-precision seismic wavelet results from the input seismic and well logging data, which will affect the accuracy of subsequent well-seismic calibration, seismic inversion, and oil and gas reservoir prediction results.

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

[0004] To address the aforementioned technical problems, this invention proposes a seismic wavelet estimation method, system, electronic device, and storage medium.

[0005] The first aspect of this invention discloses a seismic wavelet estimation method, the method comprising:

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

[0007] Step S2: Extract a time-domain seismic record corresponding to the well location coordinates from the three-dimensional seismic data of the work area, and extract a time-domain seismic record segment corresponding to the time range and 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 the pre-established seismic wavelet model to generate the first seismic wavelet;

[0009] Step S4: Establish the 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;

[0010] Step S5: 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 the second seismic wavelet, and then use the second seismic wavelet to create the current synthetic seismic record.

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

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

[0013] Step S8: Input the values ​​of the order, amplitude coefficient, dominant frequency parameter and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate the 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 step S3, the seismic wavelet model is a four-parameter seismic wavelet model, denoted as w, and expressed as:

[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∈(-∞,+∞), in milliseconds; u is the order, u∈[1,2,3,4,5,6,7,8,9,10], dimensionless; a is the amplitude coefficient, a∈(0,+∞), dimensionless; and f is the dominant frequency parameter. The unit is Hertz; dt is the time sampling interval of the time-domain seismic data, in milliseconds; θ is the phase angle parameter, θ∈(-180°, 180°], in degrees; w(t|u,a,f) h Let w(t|u,a,f) be the Hilbert transform of w(t|u,a,f).

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

[0018]

[0019] Where x represents the independent variable, for equation (2), exp() is an exponential function with the natural constant e as its base; P0(x) = 1, P1(x) = 2x, and P for the remaining u-th order is... u Calculate by recursion using 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 up 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 to t = [t1, ..., t2]. l ,…,t L ], L is the number of elements in the time axis vector, and L is a positive odd number, where l represents the element index in the time axis vector;

[0025] Step S34: Set the initial values ​​of amplitude coefficient a, main frequency parameter f, and 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, an 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 y, expressed as:

[0028]

[0029] Where, r env and y env , respectively, represent the time-domain reflection coefficient r and the amplitude envelope of 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 dominant frequency parameter is determined using the peak frequency of the amplitude spectrum of the time-domain seismic record segment y, expressed as:

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

[0032] Here, 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 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 denoted as F, and its expression is:

[0035]

[0036] Where * denotes the convolution operation, and σ is the standard deviation of yr*w(t|u,a,f,θ). The amplitude coefficient to be solved clock speed parameters and phase angle parameters The parameter vector formed, μ c =[afθ] T ,

[0037] The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariate optimization algorithm. clock speed parameters and phase angle parameters Right now,

[0038] A second aspect of the present invention discloses a seismic wavelet estimation system, the system comprising:

[0039] The first processing module is configured to calculate the corresponding time-domain reflection coefficient using logging velocity data and logging density data.

[0040] The second processing module is configured to extract a time-domain seismic record corresponding to the well location coordinates from the three-dimensional seismic data of the work area, and extract a time-domain seismic record segment corresponding to the time range and time-domain reflection coefficient from the time-domain seismic record.

[0041] The third processing module is 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 the first seismic wavelet;

[0042] The fourth processing module is 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] The fifth processing module is 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 the second seismic wavelet, and to produce the current synthetic seismic record based on the second seismic wavelet.

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

[0045] The seventh processing module is configured to change the initial value of the order and repeat the processing of the third to sixth processing modules until multiple correlation coefficients are obtained.

[0046] The eighth processing module is configured to input the values ​​of the order, amplitude coefficient, dominant frequency parameter, and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate the 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, the seismic wavelet model in the third processing module 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∈(-∞,+∞), in milliseconds; u is the order, u∈[1,2,3,4,5,6,7,8,9,10], dimensionless; a is the amplitude coefficient, a∈(0,+∞), dimensionless; and f is the dominant frequency parameter. The unit is Hertz; dt is the time sampling interval of the time-domain seismic data, in milliseconds; θ is the phase angle parameter, θ∈(-180°, 180°], in degrees; w(t|u,a,f) h Let w(t|u,a,f) be the Hilbert transform of w(t|u,a,f).

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

[0051]

[0052] Where x represents the independent variable, for equation (2), exp() is an exponential function with the natural constant e as its base; P0(x) = 1, P1(x) = 2x, and P for the remaining u-th order is... u Calculate by recursion using 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 as follows:

[0055] Set up 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, ..., t2]. 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 amplitude coefficient a, dominant frequency parameter f, and 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 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 , respectively, represent the time-domain reflection coefficient r and the amplitude envelope of 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 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] Here, 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 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 * denotes the convolution operation, and σ is the standard deviation of yr*w(t|u,a,f,θ). The amplitude coefficient to be solved clock speed parameters and phase angle parameters The parameter vector formed, μ c =[afθ] T ,

[0070] The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariate optimization algorithm. clock speed parameters and phase angle parameters Right now,

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

[0072] A fourth aspect of this invention discloses a computer-readable storage medium. The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of a seismic wavelet estimation method according to any one of the first aspects of this disclosure.

[0073] In summary, the proposed solution adopts a multi-parameter seismic wavelet model that expands the space of the existing Ricker wavelet model. This model has a larger spatial dimension and more flexible and varied waveforms. At the same time, based on the multi-parameter seismic wavelet model, a corresponding objective function is set to estimate the control parameters and inverse motion, thereby improving the accuracy of the estimation results in the time domain seismic wavelet estimation. Attached Figure Description

[0074] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

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

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

[0077] Figure 3 The input data used in the first example of the present invention includes: (a) well logging density data converted to the time domain, (b) well logging P-wave velocity data converted to the time domain, (c) 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) time domain seismic record segment corresponding to the time domain reflection coefficient.

[0078] Figure 4The input data used in the second example of the present invention includes: (a) well logging density data converted to the time domain, (b) well logging P-wave velocity data converted to the time domain, (c) 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) time domain seismic record segment corresponding to the time domain reflection coefficient.

[0079] Figure 5 This is a schematic diagram illustrating the determination of the initial value of the main frequency parameter in a first embodiment of the present invention;

[0080] Figure 6 This is a schematic diagram illustrating the determination of the initial value of the phase angle parameter in a first embodiment of the present invention;

[0081] Figure 7 This is a schematic diagram illustrating the determination of the initial value of the main frequency parameter in a second embodiment of the present invention;

[0082] Figure 8 This is a schematic diagram illustrating the determination of the initial value of the phase angle parameter in a second embodiment of the present invention;

[0083] Figure 9 A schematic diagram of the input data for a first instance according to an embodiment of the present invention, based on the seismic wavelet model of the present invention, and the seismic wavelet estimation result (b).

[0084] Figure 10 A schematic diagram of the input data for a first instance according to an embodiment of the present invention, based on (a) the seismic wavelet estimation result obtained from the Ricker wavelet model;

[0085] Figure 11 A schematic diagram of (b) seismic wavelet estimation results obtained based on (a) the seismic wavelet model of the present invention as input data for a second instance according to an embodiment of the present invention;

[0086] Figure 12 A schematic diagram of the input data for a second instance according to an embodiment of the present invention, based on (a) the seismic wavelet estimation result obtained from the Ricker wavelet model;

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

[0088] Figure 14 This is a structural diagram of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0089] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of the embodiments of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

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

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

[0092] Step S2: Extract a time-domain seismic record corresponding to the well location coordinates from the three-dimensional seismic data of the work area, and extract a time-domain seismic record segment corresponding to the time range and time-domain reflection coefficient from the time-domain seismic record.

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

[0094] Step S4: Establish the 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;

[0095] Step S5: 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 the second seismic wavelet, and then use the second seismic wavelet to create the current synthetic seismic record.

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

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

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

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

[0100] In step S2, the time-domain seismic record corresponding to the well location coordinates is extracted from the three-dimensional seismic data of the work area, and the time-domain seismic record segment y corresponding to the time range and time-domain reflection coefficient is further extracted from the seismic record.

[0101] In step S3, the initial values ​​of the set order, amplitude coefficient, dominant frequency parameter and phase angle parameter are input 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, represented 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∈(-∞,+∞), in milliseconds; u is the order, u∈[1,2,3,4,5,6,7,8,9,10], dimensionless; a is the amplitude coefficient, a∈(0,+∞), dimensionless; and f is the dominant frequency parameter. The unit is Hertz; dt is the time sampling interval of the time-domain seismic data, in milliseconds; θ is the phase angle parameter, θ∈(-180°, 180°], in degrees; w(t|u,a,f) h Let w(t|u,a,f) be the Hilbert transform of w(t|u,a,f).

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

[0106]

[0107] Where x represents the independent variable, for equation (2), exp() is an exponential function with the natural constant e as its base; P0(x) = 1, P1(x) = 2x, and P for the remaining u-th order is... u Calculate by recursion using the following formula

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

[0109] Specifically, Figure 2This is a schematic diagram illustrating the changes of the seismic wavelet model with different control parameters according to an embodiment of the present invention. In the first column of diagrams from the left, the amplitude coefficient a = 1, the dominant frequency parameter f = 40Hz, and the phase angle parameter θ = 0° remain constant, with orders u = 1, u = 2, and u = 10 respectively. In the second column of diagrams from the left, the order u = 2, the dominant frequency parameter f = 40Hz, and the phase angle parameter θ = 0° remain constant, with amplitude coefficients a = 2, a = 8, and a = 16 respectively. In the third column of diagrams from the left, the amplitude coefficient a = 1, the order u = 2, and the phase angle parameter θ = 0° remain constant, with dominant frequencies f = 20Hz, f = 40Hz, and f = 60Hz respectively. In the fourth column of diagrams from the left, the amplitude coefficient a = 1, the order u = 2, the dominant frequency parameter f = 40Hz remain constant, and the phase angle parameters θ = -170°, θ = 90°, and θ = 180° respectively.

[0110] Figure 2 The waveform changes of the seismic wavelet model used in this embodiment with four model control parameters are shown. When the order u = 2 and the phase angle parameter θ = 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 up the seismic wavelet model;

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

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

[0115] Step S34: Set the initial values ​​of amplitude coefficient a, main frequency parameter f, and phase angle parameter θ to 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 the ratio of the time-domain reflection coefficient r to the maximum amplitude envelope of the time-domain seismic record segment y, expressed as:

[0118]

[0119] Where, r env and y env , respectively, represent the time-domain reflection coefficient r and the amplitude envelope of the time-domain seismic record segment y.

[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] Here, 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 the 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.

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

[0126]

[0127] Where * denotes the convolution operation, and σ is the standard deviation of yr*w(t|u,a,f,θ). The amplitude coefficient to be solved clock speed parameters and phase angle parameters The parameter vector formed, μ c =[afθ] T ,

[0128] The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariate optimization algorithm. clock speed parameters and phase angle parameters Right now,

[0129] Specifically, using two concrete examples, the method of this embodiment of the invention to obtain seismic wavelet estimation results from the input data of the first and second examples includes the following main steps:

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

[0131] (2) Extract the time-domain seismic record corresponding to the well location coordinate from the three-dimensional seismic data of the work area, and further extract the time-domain seismic record segment y corresponding to the time range and time-domain reflection coefficient from the seismic record.

[0132] Figure 3 The input data used in the first instance are as follows: (a) is the well logging density data converted to the time domain, (b) is the well logging P-wave velocity data 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 The input data used in the second example are: (a) well logging density data converted to the time domain, (b) well logging P-wave velocity data converted to the time domain, (c) time domain reflection coefficients calculated based on the well logging density data and well logging P-wave velocity data shown in (a) and (b), and (d) time domain seismic record segments corresponding to the time domain reflection coefficients.

[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 represents time, t∈(-∞,+∞), and the unit is milliseconds (ms); u represents the order, u∈[1,2,3,4,5,6,7,8,9,10], which is dimensionless; a represents the amplitude coefficient, a∈(0,+∞), which is dimensionless; and f represents the dominant frequency parameter. The unit is Hertz (Hz); dt is the time sampling interval of the time-domain seismic data, in milliseconds (ms); θ is the phase angle parameter, θ∈(-180°, 180°], in degrees (°); w(t|u,a,f) h The Hilbert transform of w(t|u,a,f) is given by the following analytical expression:

[0137]

[0138] Where exp() is an exponential function with the natural constant e as its base; P0(x) = 1, P1(x) = x, and P for the remaining u-th order are also exponential functions. u Calculate by recursion using 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 as t = [t1, ..., t2]. l ,…,t L ], 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 amplitude coefficient, dominant frequency parameter, and phase angle parameter to a0, f0, and θ0, respectively. The initial value of amplitude coefficient a0 can be set using the ratio of the time-domain reflection coefficient r to the maximum amplitude envelope of the time-domain seismic record segment y, i.e.

[0143]

[0144] In the formula, r env and y env Let r be the time-domain reflection coefficient and y be the amplitude envelope of the time-domain seismic record segment y, respectively. The initial value of the dominant frequency parameter f0 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] In the formula, fft() represents the Fourier transform operation.

[0147] For the initial value θ0 of the phase angle parameter, the kurtosis phase estimation method can be used to estimate it from the time-domain seismic record segment y: within a certain angle range, using each angle value, the time-domain seismic record segment y is first phase-rotated to obtain y. θ Then calculate the corresponding kurtosis. θ :

[0148]

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

[0150] Figure 5 This 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. The determined initial value of the dominant frequency parameter f0 = 24.1 Hz.

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

[0152] Figure 7 This 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 second example. The determined initial value of the dominant frequency parameter f0 = 30.94 Hz.

[0153] Figure 8 This is a schematic diagram of determining the initial value θ0 of the phase angle parameter from the time-domain seismic record segment y using the kurtosis phase estimation method in the second example. The determined initial value of the phase angle parameter θ0 = 35°. The angle 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 * denotes the convolution operation, and σ is the standard deviation of yr*w(t|u,a,f,θ). The amplitude coefficient to be solved clock speed parameters and phase angle parameters The parameter vector formed, μ c =[afθ] T , The amplitude coefficients that minimize the objective function F can be obtained using the Nelder-Mead simplex search method (Lagarias et al., 1998) or the quasi-Newton method. clock speed parameters and phase angle parameters Right now,

[0158] (9) Order Then, a seismic wavelet w(t|u,a,f,θ) is generated, and a synthetic seismic record is created.

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

[0160] (11) Let u = u + 1, return to step (6), 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 In the [ ], select the largest correlation coefficient k opt =max[k1,k2,k3,k4,k5,k6,k7,k8,k9,k 10 ], and k opt The corresponding order u opt Amplitude coefficient a opt , main frequency parameter f opt and phase angle parameter θ opt And generate seismic wavelet w(t|u opt ,a opt ,f opt ,θ opt This serves as the final seismic wavelet estimation result.

[0162] Tables 1 and 2 present the estimation results of the seismic wavelet model control parameters obtained from the input data of the first and second instances using the method of the present invention.

[0163] Figure 9 A schematic diagram of the input data for a first instance according to an embodiment of the present invention, based on the seismic wavelet model of the present invention, and the seismic wavelet estimation result (b). Figure 10 A schematic diagram of the input data for a first instance according to an embodiment of the present invention, based on (a) the seismic wavelet estimation result obtained from the Ricker wavelet model and (b) the seismic wavelet estimation result.

[0164] For the first instance, k opt =90.56%, u opt =1, a opt =1.13×10 6 f opt =18.91Hz, θ opt =17.17°, corresponding to the seismic wavelet w(t|u opt ,a opt ,f opt ,θ opt ),like Figure 9 As shown in (a) and (b), the seismic wavelet estimation results are illustrated. Figure 9 Synthetic seismic records (as shown in (a)) were created. Figure 9 (b) shows the dashed line and the time-domain seismic record segment y( Figure 9 The correlation coefficient between the solid lines shown in (b) is 90.56%; (Refer to...) Figure 10As shown in (a) and (b), the seismic wavelet estimation results obtained using the conventional Ricker wavelet model ( Figure 10 Synthetic seismic records (as shown in (a)) were created. Figure 10 (b) shows the dashed line and the time-domain seismic record segment y( Figure 10 The correlation coefficient between the solid lines shown in (b) is only 84.18%.

[0165] Figure 11 A schematic diagram of (b) seismic wavelet estimation results obtained based on (a) the seismic wavelet model of the present invention as input data for a second instance according to an embodiment of the present invention; Figure 12 A schematic diagram of (b) seismic wavelet estimation results obtained based on (a) the Ricker wavelet model, for the second instance input data according to an embodiment 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°, corresponding to the seismic wavelet w(t|u opt ,a opt ,f opt ,θ opt ),like Figure 11 As shown in (a) and (b), the results of seismic wavelet estimation ( Figure 11 Synthetic seismic records (as shown in (a)) were created. Figure 11 (b) shows the dashed line and the time-domain seismic record segment y( Figure 11 The correlation coefficient between the solid lines shown in (b) is 67.34%; (Refer to...) Figure 12 As shown in (a) and (b), the seismic wavelet estimation results obtained using the conventional Ricker wavelet model ( Figure 12 Synthetic seismic records (as shown in (a)) were created. Figure 12 (b) shows the dashed line and the time-domain seismic record segment y( Figure 12 The correlation coefficient between the solid lines shown in (b) is only 43.54%. In the above figures, the time sampling interval for seismic wavelet, well logging density, well logging P-wave velocity, reflection coefficient, seismic record fragments, and synthetic seismic records is dt = 2 ms. It can be seen that, compared with existing methods, the seismic wavelet estimation method using this embodiment of the invention obtains more accurate seismic wavelet estimation results.

[0167] Table 1. Seismic wavelet estimation results of the first embodiment 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 proposed solution adopts a multi-parameter seismic wavelet model that expands the space of the existing Ricker wavelet model. This model has a larger spatial dimension and more flexible and varied waveforms. At the same time, based on the multi-parameter seismic wavelet model, a corresponding objective function is set to estimate the control parameters and inverse motion, thereby improving the accuracy of the estimation results in the time domain seismic wavelet estimation.

[0173] A second aspect of the present invention discloses a seismic wavelet estimation system. Figure 13 This is a structural diagram of a seismic wavelet estimation system according to an embodiment of the present invention; as shown below. Figure 13 As shown, the system 100 includes:

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

[0175] The second processing module 102 is configured to extract a time-domain seismic record corresponding to the well location coordinates from the three-dimensional seismic data of the work area in the time domain, and to extract a time-domain seismic record segment corresponding to the time range and 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, dominant frequency parameter and phase angle parameter into a pre-established seismic wavelet model to generate the 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, dominant frequency parameter and phase angle parameter that minimize the value of the objective function;

[0178] The fifth processing module 105 is 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 the second seismic wavelet, and to produce the current synthetic seismic record based on the second seismic wavelet.

[0179] The sixth processing module 106 is configured to calculate the correlation coefficient between time-domain seismic record segments 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 third to sixth processing modules until multiple correlation coefficients are obtained.

[0181] The eighth processing module 108 is configured to input the values ​​of the order, amplitude coefficient, dominant frequency parameter and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate the 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, the seismic wavelet model in the third processing module 103 is a four-parameter seismic wavelet model, 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∈(-∞,+∞), in milliseconds; u is the order, u∈[1,2,3,4,5,6,7,8,9,10], dimensionless; a is the amplitude coefficient, a∈(0,+∞), dimensionless; and f is the dominant frequency parameter. The unit is Hertz; dt is the time sampling interval of the time-domain seismic data, in milliseconds; θ is the phase angle parameter, θ∈(-180°, 180°], in degrees; w(t|u,a,f) h Let w(t|u,a,f) be the Hilbert transform of w(t|u,a,f).

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

[0186]

[0187] Where x represents the independent variable, for equation (2), exp() is an exponential function with the natural constant e as its base; P0(x) = 1, P1(x) = 2x, and P for the remaining u-th order is... u Calculate by recursion using 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 as follows:

[0190] Set up 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 to t = [t1, ..., t2]. l ,…,t L ], L is the number of elements in the time axis vector, and L is a positive odd number;

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

[0194] Let a = a0, f = f0, θ = θ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, 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 y, expressed as:

[0196]

[0197] Where, r env and y env , respectively, represent the time-domain reflection coefficient r and the amplitude envelope of the time-domain seismic record segment y.

[0198] According to the system of the second aspect of the present invention, in the third processing module 103, 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:

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

[0200] Here, 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, 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.

[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 * denotes the convolution operation, and σ is the standard deviation of yr*w(t|u,a,f,θ). The amplitude coefficient to be solved clock speed parameters and phase angle parameters The parameter vector formed, μ c =[afθ] T ,

[0205] The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariate optimization algorithm. clock speed parameters and phase angle parameters Right now,

[0206] A third aspect of this invention discloses an electronic device. The electronic device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements the steps of a seismic wavelet estimation method according to any one of the first aspects of this invention.

[0207] Figure 14 This is a structural diagram of an electronic device according to an embodiment of the present invention, such as... Figure 14 As shown, the electronic device includes a processor, memory, communication interface, display screen, and input device connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, carrier networks, Near Field Communication (NFC), or other technologies. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, buttons, a trackball, or a touchpad mounted on the device's casing, or an external keyboard, touchpad, or mouse.

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

[0209] A fourth aspect of this invention discloses a computer-readable storage medium. The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of a seismic wavelet estimation method according to any one of the first aspects of this invention.

[0210] Please note that the technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments have been described. However, as long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification. The above embodiments only illustrate several implementation methods of this 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 pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims.

[0211] The above are preferred embodiments of the present invention. It should be noted that, for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A seismic wavelet estimation method, characterized in that, The method includes: Step S1: Calculate the corresponding time-domain reflection coefficient using logging velocity data and logging density data; Step S2: Extract a time-domain seismic record corresponding to the well location coordinates from the three-dimensional seismic data of the work area, and extract a time-domain seismic record segment corresponding to the time range and time-domain reflection coefficient from the time-domain seismic record. Step S3: Input the initial values ​​of the set order, amplitude coefficient, dominant frequency parameter and phase angle parameter into the pre-established seismic wavelet model to generate the first seismic wavelet; Step S4: Establish the 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; Step S5: 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 the second seismic wavelet, and then use the second seismic wavelet to create the current synthetic seismic record. Step S6: Calculate the correlation coefficient between the time-domain seismic record fragment and the current synthetic seismic record; Step S7: Change the initial value of the order and repeat steps S3-S6 above until multiple correlation coefficients are obtained. Step S8: Input the values ​​of the order, amplitude coefficient, dominant frequency parameter and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate the third seismic wavelet and use it as the seismic wavelet estimation result; The seismic wavelet model is a four-parameter seismic wavelet model, expressed as follows: w The expression is as follows: Where t is time, The unit is milliseconds; u is the order. , dimensionless; a is the amplitude coefficient, Dimensionless; Main frequency parameters, The unit is Hertz; dt is the time sampling interval of the time-domain seismic data, in milliseconds; θ is the phase angle parameter. The unit is degrees; for Hilbert transform; The specific analytical expression is: in, x represents the independent variable; exp() is an exponential function with the natural constant e as its base. , For the remaining u-order times Calculate recursively using the following formula: ； The objective function is represented by F, and its expression is: in, This represents the convolution operation. for The standard deviation of , r is the time-domain reflection coefficient, and y is the time-domain seismic record segment. The amplitude coefficient to be solved Main frequency parameters and phase angle parameters The parameter vector formed , ; The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariate optimization algorithm. Main frequency parameters and phase angle parameters ,Right now, .

2. The seismic wavelet estimation method according to claim 1, characterized in that, Step S3 specifically includes: Step S31: Set up 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... , , L The number of elements in the time axis vector, and L It is a positive odd number. l Indicates the element index in the time axis vector; Step S34: Set the initial values ​​for the amplitude coefficient a, the dominant frequency parameter f, and the phase angle parameter θ as follows: a 0、 f 0、 θ 0; Step S35: Let a=a0, f=f0, θ=θ0, and substitute them into the seismic wavelet model to generate the first seismic wavelet.

3. The seismic wavelet estimation method according to claim 2, characterized in that, In step S3, the initial value of the amplitude coefficient is set using the ratio of the time-domain reflection coefficient r to the maximum amplitude envelope of the time-domain seismic record segment y. a 0 is represented as: in, and , respectively, represent the time-domain reflection coefficient r and the amplitude envelope of the time-domain seismic record segment y.

4. The seismic wavelet estimation method according to claim 2, characterized in that, In step S3, Initial values ​​of the dominant frequency parameter are determined using the peak frequency of the amplitude spectrum of a time-domain seismic record segment y. , is represented as: Here, fft() represents the Fourier transform operation.

5. The seismic wavelet estimation method according to claim 2, characterized in that, In 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.

6. A seismic wavelet estimation system, characterized in that, The system includes: The first processing module is configured to calculate the corresponding time-domain reflection coefficient using logging velocity data and logging density data. The second processing module is configured to extract a time-domain seismic record corresponding to the well location coordinates from the three-dimensional seismic data of the work area, and extract a time-domain seismic record segment corresponding to the time range and time-domain reflection coefficient from the time-domain seismic record. The third processing module is 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 the first seismic wavelet; The fourth processing module is 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. The fifth processing module is 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 the second seismic wavelet, and to produce the current synthetic seismic record based on the second seismic wavelet. The sixth processing module is configured to calculate the correlation coefficient between time-domain seismic record segments and the current synthetic seismic record; The seventh processing module is configured to change the initial value of the order and repeat the processing of the third to sixth processing modules until multiple correlation coefficients are obtained. The eighth processing module is configured to input the values ​​of the order, amplitude coefficient, dominant frequency parameter and phase angle parameter corresponding to the maximum correlation coefficient into the seismic wavelet model to generate the third seismic wavelet and use it as the seismic wavelet estimation result. The seismic wavelet model is a four-parameter seismic wavelet model, expressed as follows: w The expression is as follows: Where t is time, The unit is milliseconds; u is the order. , dimensionless; a is the amplitude coefficient, Dimensionless; Main frequency parameters, The unit is Hertz; dt is the time sampling interval of the time-domain seismic data, in milliseconds; θ is the phase angle parameter. The unit is degrees; for Hilbert transform; The specific analytical expression is: in, x represents the independent variable; exp() is an exponential function with the natural constant e as its base. , For the remaining u-order times Calculate recursively using the following formula: ； The objective function is represented by F, and its expression is: in, This represents the convolution operation. for The standard deviation of , r is the time-domain reflection coefficient, and y is the time-domain seismic record segment. The amplitude coefficient to be solved Main frequency parameters and phase angle parameters The parameter vector formed , ; The amplitude coefficient that minimizes the value of the objective function F is obtained using an unconstrained multivariate optimization algorithm. Main frequency parameters and phase angle parameters ,Right now, .

7. An electronic device, characterized in that, The electronic device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements the steps of the seismic wavelet estimation method according to any one of claims 1 to 5.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of a seismic wavelet estimation method according to any one of claims 1 to 5.

Images