GEOLOGICAL STRUCTURE CHARACTERIZATION METHOD BASED ON SEISMIC VELOCITY SIGNAL AND ACCELERATION SIGNAL

Abstract

A geological structure characterization method based on a seismic velocity signal and a seismic acceleration signal is provided. The seismic acceleration signal and the seismic velocity signal are acquired and normalized. Constant-phase wavelet and minimum-phase wavelet extraction or mixed-phase wavelet extraction is performed. In the constant-phase wavelet and minimum-phase wavelet extraction, a first minimum-phase wavelet is extracted based on the seismic velocity signal, and a constant-phase wavelet is extracted based on the seismic acceleration signal, and converted to a second minimum-phase wavelet. A residual between the first minimum-phase wavelet and the second minimum-phase wavelet is calculated. If the residual is greater than a preset threshold, the constant-phase wavelet extraction is performed again, otherwise, the first minimum-phase wavelet and the constant-phase wavelet are output.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of priority from Chinese Patent Application No. 202310807076.4, filed on Jul. 4, 2023. The content of the aforementioned application, including any intervening amendments thereto, is incorporated herein by reference in its entirety.

TECHNICAL FIELD

This application relates to petroleum geophysical exploration, and more specifically to a geological structure characterization method based on a seismic velocity signal and a seismic acceleration signal, which involves wavelet extraction of the seismic velocity signal and the seismic acceleration signal and deconvolution.

BACKGROUND

Petroleum and natural gas are important strategic energy resources, but in recent years, the increasing dependence on petroleum import has greatly aggravated the supply and security problems of petroleum and natural gas. In view of this, extensive attempts have been made to enhance the exploration and production of petroleum and natural gas. Reflective seismic exploration is proposed to obtain the structure and lithology of strata by using reflected wave data, while the recognition accuracy of the structure and lithology of strata is seriously restricted by the unsatisfied characterization accuracy of geological structure.

At present, the deconvolution of the extracted seismic wavelets is commonly used in the fine characterization of geological structure. In a conventional method, only one type of signal (velocity signal or acceleration signal) is acquired and processed to characterize the underground structure. Obviously, this method cannot effectively take characteristics of the two signals into account, and is dependent on the phase correction factor during the wavelet extraction. Therefore, it is urgent to develop a wavelet extraction method with improved accuracy to improve the characterization accuracy of subsurface structures.

SUMMARY

An objective of the present disclosure is to provide a geological structure characterization method based on a seismic velocity signal and an acceleration signal to overcome the above technical problems.

Technical solutions of the present disclosure are described below.

A geological structure characterization method based on a seismic velocity signal and a seismic acceleration signal, comprising:

• • acquiring the seismic acceleration signal and the seismic velocity signal;
  • normalizing the seismic acceleration signal and the seismic velocity signal; and
  • based on wavelet phase requirements, setting parameters to carry out constant-phase wavelet and minimum-phase wavelet extraction or to extract a mixed-phase wavelet.

In some embodiments, the constant-phase wavelet and minimum-phase wavelet extraction is performed through the following steps:

• • (S1) extracting a first minimum-phase wavelet based on a normalized seismic velocity signal;
  • (S2) extracting a constant-phase wavelet based on a normalized seismic acceleration signal; and converting the constant-phase wavelet to obtain a second minimum-phase wavelet; and
  • (S3) calculating a residual between the first minimum-phase wavelet and the second minimum-phase wavelet;
  • if the residual is greater than a preset threshold, returning to step (S2) and performing constant-phase wavelet extraction again;
  • if the residual is less than or equal to the preset threshold, outputting the first minimum-phase wavelet and the constant-phase wavelet.

In some embodiments, in step (S1), the first minimum-phase wavelet is extracted based on the normalized seismic velocity signal through the following formulas:

$a\left(n\right)=u\left(n\right)a_o\left(n\right);$ $u\left(n\right)=\{\begin{array}{cc}0& n<0\\ 1& n=0\\ 2& n>0\end{array};$ $a_o\left(n\right)=\mathrm{IFFT}\left(A_o\left(\omega \right)\right);\mathrm{and}$ $A_o\left(\omega \right)=\ln ❘X\left(\omega \right)❘;$

• • wherein a(n) represents the first minimum-phase wavelet; u(n) represents a coefficient; a_o(n) is an odd-part sequence of a(n); n represents a time sequence; IFFT represents Fourier inverse transform; A_o(ω) represents a Fourier transform result of a_o(n); and |X(ω)| represents an amplitude spectrum of a seismic record.

In some embodiments, in step (S2), the constant-phase wavelet is extracted based on the normalized seismic acceleration signal through the following formulas:

$g\left(a,n\right)=b\left(n\right)\cos \left(a\right)-h_b\left(n\right)\sin \left(a\right);\mathrm{and}$ $b\left(n\right)=\mathrm{IFFT}\left(❘X\left(\omega \right)❘\right);$

• • wherein g(a, n) represents the constant-phase wavelet; b(n) represents a zero-phase wavelet; a represents a phase correction factor; h_b(n) is a Hilbert transform result of b(n); IFFT represents Fourier inverse transform; and |X(ω)| represents an amplitude spectrum of a seismic record.

In some embodiments, the phase correction factor is expressed as:

$a=\frac{\pi }{k},k=1,2,3\dots ;$

• • wherein k is a positive integer.

In some embodiments, the mixed-phase wavelet is extracted through the following formula:

$w=\mathrm{IFFT}\left(W\right);\mathrm{and}$ $W=\sin \left({\Phi }^{\prime }\right)❘W❘+j\cos \left({\Phi }^{\prime }\right)❘W❘;$

• • wherein w represents the mixed-phase wavelet; IFFT represents Fourier inverse transform; W represents a wavelet frequency spectrum; Φ′ represents a wavelet phase spectrum; |W| represents a wavelet amplitude spectrum; and j is an imaginary number.

In some embodiments, the wavelet phase spectrum Φ′ is calculated through the following formulas:

${\Phi }^{\prime }={\left(A_{\mathrm{pha}}^T{A}_{\mathrm{pha}}\right)}^{-1}{A}_{\mathrm{pha}}^T{\phi }^{\prime };$ $A_{\mathrm{pha}}=\left[\begin{array}{ccccccc}2& -1& 0& 0& \dots & 0& 0\\ 1& 1& -1& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 0& 0& 0& \dots & 1& -1\\ 0& 2& 0& -1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \ddots & 0& 1\end{array}\right];\mathrm{and}$ ${\phi }^{\prime }={\left({\phi }^{\prime }\left(1,1\right),{\phi }^{\prime }\left(1,2\right),\dots ,{\phi }^{\prime }\left(1,N-1\right),{\phi }^{\prime }\left(2,2\right),\dots ,{\phi }^{\prime }\left(2,N-2\right),\dots ,{\phi }^{\prime }\left(\frac{N}{2},\frac{N}{2}\right)\right)}^T;$

• • wherein A_pha is a (N²/4)×N matrix; T represents a matrix transpose operation; −1 denotes a matrix inversion operation; q′ is a (N²/4)×1 column vector; and N represents a size of a higher-order cumulant matrix.

In some embodiments, the amplitude spectrum |W| of the mixed-phase wavelet is calculated through the following formulas:

$❘W\left(n\right)❘=\{\begin{array}{c}{e}^{W^{\prime }\left(n\right)},n=1,2,3,\dots ,N/2\\ ❘W\left(N-n\right)❘,n=\frac{N}{2}+1,\frac{N}{2}+2,\dots ,N-1\end{array};$ $W^{\prime }={\left(A_{amp}^T{A}_{amp}\right)}^{-1}{A}_{amp}^T{X}^{\prime };$ $A_{\mathrm{amp}}=\left[\begin{array}{ccccccc}2& 1& 0& 0& \dots & 0& 0\\ 1& 1& 1& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 0& 0& 0& \dots & 1& 0\\ 0& 2& 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \ddots & 0& 1\end{array}\right];\mathrm{and}$ $X^{\prime }={\left(X^{\prime }\left(1,1\right),X^{\prime }\left(1,2\right),\dots ,X^{\prime }\left(1,\frac{N}{2}-1\right),X^{\prime }\left(2,2\right),\dots ,X^{\prime }\left(2,\frac{N}{2}-2\right),\dots ,X^{\prime }\left(\frac{N}{4},\frac{N}{4}\right)\right)}^T;$

• • wherein |W(n)| represents the wavelet amplitude spectrum; e represents a base of a natural logarithm function; W′(n) represents a natural logarithm of the wavelet amplitude spectrum; n represents a time sequence; N represents a size of a higher-order cumulant matrix; |W(N−n)| represents a part of the wavelet amplitude spectrum greater than one half of a wavelet length; W′ represents a (N/2)×1 column vector, A_amp represents a (N²/16)×(N/2) matrix; T represents a matrix transpose operation; −1 represents a matrix inversion operation; and X′ represents a (N²/16)×1 column vector.

The beneficial effects of the present disclosure are described below.

In the present disclosure, wavelets are extracted based on combination of seismic velocity signals and seismic acceleration signal, which can improve the accuracy of wavelet extraction of seismic signals, and solve the problems in the existing technology that the wavelet extraction only involves a single signal, and is excessively dependent on the selected signal and the phase correction factor.

The deconvolution process is performed as follows.

The minimum-phase wavelet a or the constant-phase wavelet g or the mixed-phase wavelet w obtained by the combined wavelet extraction are used to calculate the corresponding inverse wavelets @′, g′ and we′, respectively. The relationship between a wavelet and its corresponding inverse wavelet (taking the mixed-phase wavelet w as an example) is expressed by:

$\begin{array}{cc}w*w^{\prime }=1;& \left(1\right)\end{array}$

• • where * represents a convolution symbol.

Therefore, the inverse wavelet w′ of the mixed-phase wavelet can be expressed as:

$\begin{array}{cc}{w}^{\prime }=\mathrm{IFFT}\left(2\mathrm{\pi \delta }\left(\omega \right)-\mathrm{FFT}\left(\omega \right)\right);& \left(2\right)\end{array}$

• • where δ(ω) represents an impulse function.

The reflection coefficient r can be obtained by convolution of the inverse wavelet w obtained by formula (2) with the seismic record x, expressed as:

$\begin{array}{cc}r=w^{\prime }*x.& \left(3\right)\end{array}$

The reflection coefficient obtained by formula (3) is subjected to convolution with the high-frequency wavelet to obtain the fine characterization result.

On the fine characterization image, the number of traces and time for recognizing pinch-out points are represented by Trc and Tim, respectively. Assuming that the trace spacing is dx (the distance between two neighboring seismic traces), the transverse position X_j of the recognized pinch-out point is:

$\begin{array}{cc}\mathrm{X\_j}=\mathrm{Trc}*\mathrm{dx}.& \left(4\right).\end{array}$

Assuming that the average velocity from the surface to the pinch-out point stratum is v, the longitudinal position of the recognized pinch-out point D_j is:

$\begin{array}{cc}\mathrm{D\_j}=\mathrm{Tim}*v.& \left(5\right).\end{array}$

The pinch-out point location (X_j, D_j) after fine characterization can be obtained through formulas (4) and (5).

The present disclosure can improve the resolution of seismic data through fine characterization of geological structure based on joint extraction of seismic velocity signals and seismic acceleration signals, which can realize more fine and accurate identification of geological formations, and solve the problems of excessive dependence on the phase correction factor during the characterization and wavelet extraction using only a single kind of signal in the existing technology.

BRIEF DESCRIPTION OF THE DRAWINGS

To more clearly illustrate the technical solutions of the embodiments of the present disclosure or technical solutions in the prior art, the accompanying drawings needed in the description of the embodiments or the prior art will be introduced briefly below. Obviously, presented in the accompanying drawings are only some embodiments of the present disclosure, and other accompanying drawings can be obtained by one of ordinary skill in the art based on these drawings without creative labor.

FIG. 1 is a flow chart of a geological structure characterization method based on seismic velocity signals and seismic acceleration signals according to an embodiment of the present disclosure;

FIGS. 2a-2b schematically show seismic signals according to Embodiment 1 of the present disclosure, where (a): velocity signal; and (b) acceleration signal;

FIGS. 3a-3b schematically show a velocity wavelet and an acceleration wavelet extracted in the characterization method according to Embodiment 1 of the present disclosure, where (a): velocity wavelet; and (b): acceleration wavelet;

FIG. 4 schematically show a mixed-phase wavelet extracted in the characterization method according to Embodiment 1 of the present disclosure;

FIGS. 5a-5b schematically show seismic signals according to Embodiment 2 of the present disclosure, where (a): velocity signal; and (b) acceleration signal;

FIGS. 6a-6b schematically show a velocity wavelet and an acceleration wavelet extracted in the characterization method according to Embodiment 2 of the present disclosure, where (a): velocity wavelet; and (b): acceleration wavelet;

FIGS. 7a-7b schematically show a velocity wavelet and an acceleration wavelet respectively extracted by an existing method, where (a): velocity wavelet; and (b): acceleration wavelet;

FIG. 8 schematically shows a mixed-phase wavelet extracted by the characterization method according to Embodiment 2 of the present disclosure;

FIG. 9 schematically shows a mixed-phase wavelet extracted by an existing method;

FIGS. 10a-10b schematically show fine characterization results according to Embodiment 2 of the present disclosure, where (a): fine characterization results acquired by the method of the present disclosure; and (b) fine characterization results acquired by a conventional method;

FIGS. 11a-11b respectively are enlarged views of FIG. 10a and FIG. 10b;

FIG. 12 is a partially enlarged view of initial data in Embodiment 2 of the present disclosure; and

FIGS. 13a-13b schematically show fine characterization results according to Embodiment 2 of the present disclosure, where (a): fine characterization results acquired by the method of the present disclosure; and (b) fine characterization results acquired by a conventional method.

DETAILED DESCRIPTION OF EMBODIMENTS

The technical solutions in the embodiments of the present disclosure will be described clearly and completely with reference to the accompanying drawings. Individual embodiments and the technical features recited therein can be combined in the case of no contradiction. Unless otherwise specified, all technical and scientific terms used in this application have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. As used herein, term “including” or “comprising” is intended to mean that the preceding component or object covers the following components or objects without excluding other components or objects.

As shown in FIG. 1, this application provides a geological structure characterization method based on joint extraction of a seismic velocity signal and a seismic acceleration signal, which includes the following steps. The seismic acceleration signal and the seismic velocity signal are acquired respectively through an acceleration seismometer and a velocity seismometer, and normalized. Parameters are set based on wavelet phase requirements to carry out constant-phase wavelet and minimum-phase wavelet extraction or to extract a mixed-phase wavelet.

The acceleration signal has a huge advantage over the velocity signal in terms of bandwidth and octave range, while the velocity signal has richer effective information in the low-frequency part than the acceleration signal. The present disclosure combines the velocity signal and the acceleration signal for wavelet extraction, which preserves factors, such as bandwidth, and the low-frequency effective signals, thereby improving the accuracy of the wavelet extraction results.

It should be noted that the normalized processing of the signal is known in the prior art, and the specific steps will not be repeated herein. The setting of parameters, including the number of processing channels, time windows, etc., is also known in the prior art, thus the specific details are not repeated herein.

In an embodiment, the extraction of the constant-phase wavelet and the extraction of the minimum-phase wavelets are carried out through the following steps.

(S1) A first minimum-phase wavelet is extracted based on the seismic velocity signal to obtain through the following formulas:

$\begin{array}{cc}a\left(n\right)=u\left(n\right)a_o\left(n\right);& \left(1\right)\end{array}$ $\begin{array}{cc}u\left(n\right)=\{\begin{array}{cc}0& n<0\\ 1& n=0\\ 2& n>0\end{array};& \left(2\right)\end{array}$ $\begin{array}{cc}{a}_o\left(n\right)=\mathrm{IFFT}\left(A_o\left(\omega \right)\right);\mathrm{and}& \left(3\right)\end{array}$ $\begin{array}{cc}{A}_o\left(\omega \right)=\ln ❘X\left(\omega \right)❘;& \left(4\right)\end{array}$

• • where a(n) represents the first minimum-phase wavelet; u(n) represents a coefficient; a_o(n) is an odd part sequence of a(n); n represents a time sequence; IFFT represents Fourier inverse transform; A_o(ω) represents a Fourier transform result of a_o(n); and |X(ω)| represents an amplitude spectrum of a seismic record.

In the above embodiment, the extraction formula for the first minimum-phase wavelet is derived by the following steps.

Assuming that the seismic wavelet is b(n), the logarithmic spectrum of the wavelet A(ω) can be written as:

$\begin{array}{cc}A\left(\omega \right)=\ln B\left(\omega \right);& \left(17\right)\end{array}$

• • where B(ω) represents a frequency spectrum of the seismic wavelet; and ln represents a logarithm function with a base of e.

Assuming that a(n) is the time sequence of the logarithmic spectrum of the seismic wavelet b(n) (i.e., the minimum phase wavelet). It is theoretically proven that if a wavelet is a minimum-phase wavelet, the time sequence of its logarithmic spectrum is a real-valued causal sequence. Any real-valued sequence can be written as the sum of its odd-part sequence and its even-part sequence, so a(n) can be written as:

$\begin{array}{cc}a\left(n\right)=a_o\left(n\right)+a_e\left(n\right);& \left(18\right)\end{array}$

• • where a_e(n) represents an even-part sequence of the time sequence a(n) of the log-spectrum.

Since a(n) has causality, a(n) can be written as the expression shown in formula (1), and u(n) in formula (1) is shown in formula (2). Hence, only the odd-part sequence a_o(n) is obtained, and the minimum-phase wavelet can be obtained by formula (1).

Taking logarithms for the left and right sides of formula (31), the following formula is obtained:

$\begin{array}{cc}\ln ❘B\left(\omega \right)❘=\ln ❘X\left(\omega \right)❘.& \left(19\right)\end{array}$

Meanwhile, the frequency spectrum of the wavelet can be expressed as:

$\begin{array}{cc}B\left(\omega \right)=❘B\left(\omega \right)❘e^{i\phi \left(\omega \right)};& \left(20\right)\end{array}$

• • where φ(ω) represents a phase of the wavelet.

Then, formula (19) is substituted into formula (17) to obtain:

$\begin{array}{cc}A\left(\omega \right)=\ln B\left(\omega \right)=\ln ❘B\left(\omega \right)❘+i\phi \left(\omega \right).& \left(21\right)\end{array}$

At the same time, there are:

$\begin{array}{cc}{A}_o\left(\omega \right)=A_R\left(\omega \right);& \left(22\right)\end{array}$ $\begin{array}{cc}{A}_e\left(\omega \right)=A_I\left(\omega \right);\mathrm{and}& \left(23\right)\end{array}$ $\begin{array}{cc}A\left(\omega \right)=A_o\left(\omega \right)+{\mathrm{iA}}_e\left(\omega \right);& \left(24\right)\end{array}$

• • where A_o(ω) and A_e(ω) are Fourier transform results of a_o(n) and a_e(n), respectively; A_R(ω) and A_I(ω) are the real and imaginary parts of the wavelet logarithmic spectrum, respectively; and i represents an imaginary unit.

Formulas (21)-(24) are subjected to simultaneity to obtain:

$\begin{array}{cc}{A}_o\left(\omega \right)=\ln ❘B\left(\omega \right)❘;& \left(25\right)\end{array}$ That is, A_o(ω)=ln|B(ω)|=ln|X(ω)|  (26);

Thus, the A_o(ω) shown in formula (4) can be computed, and from this, the a_o(n) shown in formula (3) can also be computed. Finally, u(n) is determined from the value of n, and the odd-part sequence a_o(n) can be calculated by formula (3) to extract the minimum-phase wavelet based on the velocity signal through formula (1).

(S2) The constant-phase wavelet is extracted based on the seismic acceleration signal to obtain the constant-phase wavelet, and the constant-phase wavelet is converted to obtain a second minimum-phase wavelet.

In an embodiment, the constant-phase wavelet is extracted based on the seismic acceleration signal through the following formulas:

$\begin{array}{cc}g\left(a,n\right)=b\left(n\right)\cos \left(a\right)-h_b\left(n\right)\sin \left(a\right);\mathrm{and}& \left(5\right)\end{array}$ $\begin{array}{cc}b\left(n\right)=\mathrm{IFFT}\left(❘X\left(\omega \right)❘\right);& \left(6\right)\end{array}$

• • where g(a, n) represents the constant-phase wavelet; b(n) represents a zero-phase wavelet; a represents a phase correction factor; and h_b(n) is a Hilbert transform result of b(n).

Optionally, the phase correction factor is expressed as:

$a=\frac{\pi }{k},k=1,2,3\dots ;$

• • where k is a positive integer.

In the above embodiment, the constant-phase wavelet extraction formula is derived by the following steps.

Assuming that the seismic record is x(n), the reflection coefficient is r(n), and the seismic wavelet is b(n) (the seismic wavelet is a zero-phase wavelet when the wavelet phase is zero). Then, the seismic convolution model can be expressed as:

$\begin{array}{cc}x\left(n\right)=b\left(n\right)*r\left(n\right).& \left(27\right)\end{array}$

Firstly, the Fourier transform is performed on both sides of formula (27) to obtain the expression of the convolution model in the frequency domain, as shown below:

$\begin{array}{cc}X\left(\omega \right)=B\left(\omega \right)*R\left(\omega \right);& \left(28\right)\end{array}$

• • where X(ω) is the Fourier transform of x(n), B(ω) is the Fourier transform of b(n), and R(ω) is the Fourier transform of r(n).

Assuming that the reflection coefficient is a white noise sequence, the amplitude spectrum is expressed as:

$\begin{array}{cc}❘R\left(\omega \right)❘=1.& \left(29\right)\end{array}$

Then, the amplitude spectrum of the seismic record x(n) can be written as:

$\begin{array}{cc}❘X\left(\omega \right)❘=❘B\left(\omega \right)❘❘R\left(\omega \right)❘=❘B\left(\omega \right)❘;& \left(30\right)\end{array}$ that is, |B(ω)|=|X(ω)|  (31).

Formula (31) shows that the amplitude spectrum |B(ω)| of the wavelet can be estimated from the amplitude spectrum |X(ω)| of the seismic record.

Meanwhile, the frequency spectrum of the wavelet can be expressed as:

$\begin{array}{cc}B\left(\omega \right)=❘B\left(\omega \right)❘e^{i\phi \left(\omega \right)};& \left(32\right)\end{array}$

• • where φ(ω) represents the phase of the wavelet.

Secondly, when the phase of the wavelet is zero, i.e., φ(ω)=0 in formula (32), the following formula can be obtained:

$\begin{array}{cc}b\left(n\right)=\mathrm{IFFT}\left(B\left(\omega \right)\right)=\mathrm{IFFT}\left(❘B\left(\omega \right)❘\right)=\mathrm{IFFT}\left(❘X\left(\omega \right)❘\right).& \left(33\right)\end{array}$

From this, the zero-phase wavelet calculation formula shown in formula (6) can be obtained, and the zero-phase wavelet b(n) can be obtained by performing Fourier inverse transform on the amplitude spectrum |X(ω)| of the seismic record.

Again, the phase correction can be performed on the zero-phase wavelet to obtain the expression of the constant-phase wavelet in the frequency domain:

$\begin{array}{cc}G\left(a,\omega \right)=\{\begin{array}{cc}B\left(\omega \right)e^{\mathrm{ia}}=B\left(\omega \right)\left(\cos \left(a\right)+i\sin \left(a\right)\right)& \omega >0\\ B\left(\omega \right)e^{-\mathrm{ia}}=B\left(\omega \right)\left(\cos \left(a\right)-i\sin \left(a\right)\right)& \omega <0\end{array};& \left(34\right)\end{array}$

where G(α, ω) is the frequency spectrum of the constant-phase wavelet.

Next, assuming that the frequency response of the Hilbert transform factor is H(ω), then:

$\begin{array}{cc}H\left(\omega \right)=\{\begin{array}{cc}-i& \omega >0\\ i& \omega <0\end{array}.& \left(35\right)\end{array}$

Formula (35) is substituted into formula (34) to obtain:

$\begin{array}{cc}G\left(a,\omega \right)=B\left(\omega \right)\cos \left(a\right)-H\left(\omega \right)\sin \left(a\right).& \left(36\right)\end{array}$

From this, the constant-phase wavelet extraction formula shown in formula (5) can be obtained, and the constant-phase wavelet of the acceleration signal can be extracted by using this formula.

(S3) A residual between the first minimum-phase wavelet and the second minimum-phase wavelet is calculated. If the residual is greater than a preset threshold, returning to step (S2) and the constant-phase wavelet extraction is performed again. If the residual is less than or equal to the preset threshold, the first minimum-phase wavelet and the constant-phase wavelet are output.

It should be noted that when the residual is greater than a preset threshold to return to step (S2), a new constant-phase factor is extracted by resetting a new phase correction factor. The new phase correction factor can be obtained by formula (7), and the value of k is taken one place backward for each update of the phase correction factor.

It should be noted that the residual threshold is an artificially set threshold, and the smaller the residual threshold is, the higher the accuracy of the obtained results. The present disclosure utilizes that the seismic acceleration signal is equal to the derivative of the seismic velocity signal, with the first minimum-phase wavelet obtained by the extraction of the velocity signal as a benchmark, to convert the constant-phase wavelet obtained by the extraction of the acceleration signal into a second minimum-phase wavelet. By comparing the first minimum-phase wavelet with the second minimum-phase wavelet, the constant-phase factor can be determined. This can avoid the uncertainty brought about by the artificially-set constant phase factor, and ultimately improves the extraction accuracy of the constant phase wavelet.

In an embodiment, the mixed-phase wavelet is extracted through the following formula:

$\begin{array}{cc}w=\mathrm{IFFT}\left(W\right);\mathrm{and}& \left(8\right)\\ W=\sin \left({\Phi }^{\prime }\right)❘W❘+j\cos \left({\Phi }^{\prime }\right)❘W❘;& \left(9\right)\end{array}$

where w represents the mixed-phase wavelet; W represents a wavelet frequency spectrum; Φ′ represents a wavelet phase spectrum; |W| represents a wavelet amplitude spectrum; and j is an imaginary number.

Optionally, the phase spectrum Φ′ of the mixed-phase wavelet is calculated through the following formulas:

$\begin{array}{cc}{\Phi }^{\prime }={\left(A_{\mathrm{pha}}^T{A}_{\mathrm{pha}}\right)}^{-1}{A}_{\mathrm{pha}}^T{\phi }^{\prime };& \left(10\right)\\ A_{\mathrm{pha}}=\left[\begin{array}{ccccccc}2& -1& 0& 0& \cdots & 0& 0\\ 1& 1& -1& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 0& 0& 0& \dots & 1& -1\\ 0& 2& 0& -1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \ddots & 0& 1\end{array}\right];\mathrm{and}& \left(11\right)\\ {\phi }^{\prime }={\left({\phi }^{\prime }\left(1,1\right),{\phi }^{\prime }\left(1,2\right),\dots ,{\phi }^{\prime }\left(1,N-1\right),{\phi }^{\prime }\left(2,2\right),\dots ,{\phi }^{\prime }\left(2,N-2\right),\dots ,{\phi }^{\prime }\left(\frac{N}{2},\frac{N}{2}\right)\right)}^T;& \left(12\right)\end{array}$

where A_pha is a (N²/4)×N matrix; T represents a matrix transpose operation; −1 denotes a matrix inversion operation; φ′ is a (N²/4)×1 column vector; and N represents a size of a higher-order cumulant matrix.

Optionally, the amplitude spectrum |W| of the mixed-phase wavelet is calculated through the following formulas:

$\begin{array}{cc}❘W\left(n\right)❘=\{\begin{array}{c}{e}^{W^{\prime }\left(n\right)},n=1,2,3\cdots ,N/2\\ ❘W\left(N-n\right)❘,n=\frac{N}{2}+1,\frac{N}{2}+2,\cdots ,N-1\end{array};& \left(13\right)\\ W^{\prime }={\left(A_{\mathrm{amp}}^T{A}_{\mathrm{amp}}\right)}^{-1}{A}_{\mathrm{amp}}^T{X}^{\prime };& \left(14\right)\\ A_{\mathrm{amp}}=\left[\begin{array}{ccccccc}2& 1& 0& 0& \cdots & 0& 0\\ 1& 1& 1& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 0& 0& 0& \dots & 1& 0\\ 0& 2& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \ddots & 0& 1\end{array}\right];\mathrm{and}& \left(15\right)\\ X^{\prime }={\left(X^{\prime }\left(1,1\right),X^{\prime }\left(1,2\right),\dots ,X^{\prime }\left(1-\frac{N}{2}-1\right),X^{\prime }\left(2,2\right),\dots ,X^{\prime }\left(2,\frac{N}{2}-2\right),\dots ,X^{\prime }\left(\frac{N}{4},\frac{N}{4}\right)\right)}^T;& \left(16\right)\end{array}$

where |W(n)| represents a wavelet amplitude spectrum; e represents a base of a natural logarithm function; W′(n) represents a natural logarithm of the wavelet amplitude spectrum; n represents a time sequence; N represents a size of a higher-order cumulant matrix; |W(N−n)| represents a part of the wavelet amplitude spectrum greater than one half of a wavelet length; W′ represents a (N/2)×1 column vector; A_amp represents a (N²/16)×(N/2) matrix; T represents a matrix transpose operation; −1 represents a matrix inversion operation; and X′ represents a (N²/16)×1 column vector.

In the above embodiment, the mixed phase wavelet is extracted by utilizing the higher-order accumulation of the acceleration signal and the velocity signal, which can satisfy the characteristics of the attenuated wavelet to render the extraction results more in line with reality.

In the above embodiment, the mixed-phase wavelet extraction formula is derived by the following steps.

Assuming that there exists a zero-mean smooth stochastic process y(n), where n=0, ±1, ±2, . . . .

Then the third-order cumulative quantity of y(n) can be expressed as:

$\begin{array}{cc}{\mathrm{cum}}_{3y}\left({\tau }_1,{\tau }_2\right)=E\left[y\left(n\right)y\left(n+{\tau }_1\right)y\left(n+{\tau }_2\right)\right];& \left(37\right)\end{array}$

where E[·] denotes the calculated mathematical expectation.

The above equation can be rewritten as:

$\begin{array}{cc}{\mathrm{cum}}_{3y}\left({\tau }_1,{\tau }_2\right)={\sum }_{n=-\infty }^{\infty }y\left(n\right)y\left(n+{\tau }_1\right)y\left(n+{\tau }_2\right).& \left(38\right)\end{array}$

The k−1^th Fourier transform of the k^th order cumulative quantity is its k^th order spectrum, also called k−1 spectrum. Then, the third-order spectrum (bispectrum) of the above third-order cumulative quantity can be expressed as:

$\begin{array}{cc}{C}_{3y}\left({\omega }_1,{\omega }_2\right)={\sum }_{{\tau }_1=-\infty }^{\infty }{\sum }_{{\tau }_2=-\infty }^{\infty }{\mathrm{cum}}_{3y}\left({\tau }_1,{\tau }_2\right)e^{-j\left({\omega }_1{\tau }_1+{\omega }_2{\tau }_2\right)}.& \left(39\right)\end{array}$

The formula (38) is substituted into formula (39) to obtain:

$\begin{array}{cc}{C}_{3y}\left({\omega }_1,{\omega }_2\right)={\sum }_{n=-\infty }^{\infty }{\sum }_{{\tau }_1=-\infty }^{\infty }{\sum }_{{\tau }_2=-\infty }^{\infty }y\left(n\right)y\left(n+{\tau }_1\right)y\left(n+{\tau }_2\right)e^{-j\left({\omega }_1{\tau }_1+{\omega }_2{\tau }_2\right)};& \left(40\right)\end{array}$

i.e.:

$\begin{array}{cc}{C}_{3y}\left({\omega }_1,{\omega }_2\right)={\sum }_{n=-\infty }^{\infty }y\left(n\right)e^{jn\left({\omega }_1+{\omega }_2\right)}{\sum }_{{\tau }_1=-\infty }^{\infty }y\left(n+{\tau }_1\right)e^{-j{\omega }_1\left(n+{\tau }_1\right)}{\sum }_{{\tau }_2=-\infty }^{\infty }y\left(n+{\tau }_2\right)e^{-j{\omega }_2\left(n+{\tau }_2\right)}=Y\left({\omega }_1+{\omega }_2\right)Y\left({\omega }_1\right)Y\left({\omega }_2\right)=❘{\mathrm{cum}}_{3y}\left({\omega }_1,{\omega }_2\right)❘e^{j{\phi }_{3x}\left({\omega }_1,{\omega }_2\right)};& \left(41\right)\end{array}$

where Y is the Fourier transform of y(n).

Then, the amplitude spectrum and phase spectrum of the third-order spectrum can be respectively expressed as:

$\begin{array}{cc}❘\mathrm{cu}{m}_{3y}\left({\omega }_1,{\omega }_2\right)❘=❘Y\left({\omega }_1\right)❘❘Y\left({\omega }_2\right)❘❘Y\left({\omega }_1+{\omega }_2\right)❘;\mathrm{and}& \left(42\right)\end{array}$ $\begin{array}{cc}{\phi }_{3y}\left({\omega }_1,{\omega }_2\right)={\Phi }_y\left({\omega }_1\right)+{\Phi }_y\left({\omega }_2\right)-{\Phi }_y\left({\omega }_1+{\omega }_2\right).& \left(43\right)\end{array}$

Assuming that the seismic record is x(t), the convolution model can be expressed as:

$\begin{array}{cc}x\left(t\right)=w\left(t\right)*r\left(t\right)+v\left(t\right);& \left(44\right)\end{array}$

where w(t) represents the seismic wavelet sequence; r(t) represents the reflection coefficient sequence; v(t) represents the Gaussian noise sequence; and * represents convolution.

Formula (44) is subjected to Fourier transform to obtain:

$\begin{array}{cc}X\left(\omega \right)=W\left(\omega \right)*R\left(\omega \right)+V\left(\omega \right);& \left(45\right)\end{array}$

where X(ω), W(ω), R(ω), and V(ω) are the Fourier transforms of x(t), w(t), r(t), and v(t), respectively.

Then, the bi-spectrum of x(t) is:

$\begin{array}{cc}{C}_{3x}\left({\omega }_1,{\omega }_2\right)=C_{3w}\left({\omega }_1,{\omega }_2\right)C_{3r}\left({\omega }_1,{\omega }_2\right)+C_{3v}\left({\omega }_1,{\omega }_2\right).& \left(46\right)\end{array}$

Since the higher-order cumulative quantity of the Gaussian random variable is zero, then:

$\begin{array}{cc}{C}_{3v}\left({\omega }_1,{\omega }_2\right)=0.& \left(47\right)\end{array}$

Hence, formula (47) is substituted into formula (46) to obtain:

$\begin{array}{cc}{C}_{3x}\left({\omega }_1,{\omega }_2\right)=C_{3w}\left({\omega }_1,{\omega }_2\right)C_{3r}\left({\omega }_1,{\omega }_2\right).& \left(48\right)\end{array}$

Typically, the higher-order accumulation of the reflection coefficient is assumed to be a multidimensional impulse function, then formula (48) can be written as:

$\begin{array}{cc}{C}_{3x}\left({\omega }_1,{\omega }_2\right)=C_{3w}\left({\omega }_1,{\omega }_2\right)/{\sigma }^2;& \left(49\right)\end{array}$

where σ represents the variance of the reflection coefficient sequence.

Formula (49) shows that there is only one scale factor difference between the third-order spectrum of the wavelet and the third-order spectrum of the seismic record. Therefore, the third-order spectrum of the seismic record can be used to determine the spectrum of the wavelet.

Taking the logarithm of both sides of formula (42) to obtain:

$\begin{array}{cc}\ln ❘{\mathrm{cum}}_{3y}\left({\omega }_1,{\omega }_2\right)❘=\ln ❘Y\left({\omega }_1\right)❘+\ln ❘Y\left({\omega }_2\right)❘+\ln ❘Y\left({\omega }_1+{\omega }_2\right)❘.& \left(50\right)\end{array}$

Formula (50) is applied to the seismic data to obtain:

$\begin{array}{cc}\ln ❘{\mathrm{cum}}_{3x}\left({\omega }_1,{\omega }_2\right)❘=\ln ❘W\left({\omega }_1\right)❘+\ln ❘W\left({\omega }_2\right)❘+\ln ❘W\left({\omega }_1+{\omega }_2\right)❘.& \left(51\right)\end{array}$

Allowing ln|cum_3x(ω₁, ω₂)|=X′ (k, l) and ln|W|=W′, and formula (51) can be written as:

$\begin{array}{cc}{X}^{\prime }\left(k,l\right)=W^{\prime }\left(k\right)+W^{\prime }\left(l\right)+W^{\prime }\left(k+l\right).& \left(52\right)\end{array}$

Formula (52) shows that X′ is a linear combination of W′, so formula (52) can be written in matrix form as follows:

$\begin{array}{cc}{X}^{\prime }=A_{amp}{W}^{\prime };& \left(53\right)\end{array}$

where the specific forms of X′ and A_amp are shown in formulas (16) and (15), respectively.

A_amp is a column non-singular matrix. Therefore, the generalized inverse matrix can be used to solve formula (53) to solve the expression of W′ shown in formula (14), and the wavelet amplitude spectrum can be solved using formula (13).

Formula (13) is applied to the seismic data to obtain:

$\begin{array}{cc}{\phi }_{3x}\left({\omega }_1,{\omega }_2\right)={\Phi }_w\left({\omega }_1\right)+{\Phi }_w\left({\omega }_2\right)-{\Phi }_w\left({\omega }_1+{\omega }_2\right).& \left(54\right)\end{array}$

Formula (54) indicates that φ_3x is a linear combination of Φ_w, and thus formula (54) can be written in a matrix form as follows:

$\begin{array}{cc}{\phi }^{\prime }=A_{pha}{\Phi }^{\prime };& \left(55\right)\end{array}$

where Φ′ represents N×1 column vector.

The specific expressions of φ′ and A_pha in the formula (55) are shown in formula (12) and formula (11), respectively.

A_pha is a column non-singular matrix. Therefore, the generalized inverse matrix can be used to solve formula (55) to solve the expression for the wavelet phase spectrum shown in formula (10), which can be used to solve the wavelet phase spectrum.

The method using the above-mentioned phase recovery utilizes all possible bi-spectral phase values that can be obtained, and is a non-recursive algorithm. Thus, there is no cumulative error in the calculation process.

Based on the calculation results of formulas (10) and (13), the spectrum of the wavelet can be obtained computationally using formula (9). Then the Fourier inverse transform is performed on formula (9) to obtain the mixed-phase wavelet extraction formula shown in formula (8), and the mixed-phase wavelet of the seismic record can be solved by formula (8).

In an embodiment, a method for extracting wavelets using the above-mentioned method for combined wavelet extraction of seismic velocity signals and seismic acceleration signals is provided, which includes the following steps.

(1) A seismic acceleration signal and a seismic velocity signal are acquired and normalized.

(2) Parameters (number of processing channels, time window, etc.) are set, and the data within the time window are averaged to obtain an averaged seismic record and its amplitude spectrum.

(3) A Fourier inverse transform is performed on the amplitude spectrum obtained from the acceleration signal to obtain a zero-phase wavelet.

(4) A Hilbert transform is performed on the zero-phase wavelet, and a constant-phase factor is set to extract an acceleration wavelet (i.e., the constant-phase wavelet).

(5) A logarithmic operation is performed on the amplitude spectrum obtained from the velocity signal to obtain a logarithmic spectrum sequence.

(6) A Fourier inverse transform is performed on the logarithmic spectrum to obtain a frequency spectrum, and then a velocity wavelet (i.e., the first minimum phase wavelet) is extracted.

(7) The acceleration wavelet is converted into a velocity wavelet (i.e., the second minimum phase wavelet), and the residual between the converted velocity wavelet and the original velocity wavelet is calculated. If the residual satisfies the accuracy requirement, the calculation is stopped, and the jointly extracted velocity wavelet and acceleration wavelet are output. If the accuracy requirement is not satisfied, the phase correction factor is re-selected, and the above operations are repeated until the accuracy requirement is satisfied;

(8) After extracting the constant-phase wavelet and the minimum-phase wavelets, the data after preprocessing (normalization, parameterization, and averaging) are converted to the acceleration domain and spliced with the acceleration signal after preprocessing to obtain fused data.

(9) A third-order cumulative quantity of the fused data and the third-order spectrum (bi-spectrum) of the fused data are calculated. The third-order cumulative quantity is filtered using a Parzen window, and thereafter an amplitude spectrum and a phase spectrum of the wavelet are calculated using the bi-spectrum. Ultimately, the mixed-phase wavelet is reconstructed.

Embodiment 1

In Embodiment 1, the seismic signal of a work area is selected, as shown in FIG. 2, with the time range of 0-2 s, the number of channels of 1-853, and the residual threshold of 0.001. The final outputs of the jointly extracted velocity wavelet and acceleration wavelet are shown in FIGS. 3a-3b, respectively. The mixed-phase wavelet reconstructed by using the bi-spectrum of the higher-order cumulative quantity is shown in FIG. 4.

The mixed-phase wavelet (shown in FIG. 4) extracted by the method provided in the present disclosure, the wavelet obtained by the conventional method and the original seismic acceleration signal (shown in FIG. 2b) are subjected to the deconvolution to obtain the fine characterization results as shown in FIGS. 10a-10b, where (a): fine characterization results obtained by the method proposed in the present disclosure; and (b): fine characterization results obtained by using the conventional method. FIGS. 11a-11b respectively are enlarged results of FIG. 10a and FIG. 10b. Comparing the fine characterization results obtained by different methods in FIGS. 11a-11b, the fine characterization results obtained by the present disclosure have a clearer image at the position pointed out by the arrow, and have a better continuity, i.e., the method in the present disclosure can recognize the continuous structure at the 200^th channel at 0.45 s, but the conventional method can only recognize the 200^th channel at 0.43 s. The channel spacing is 10 m, and the speed is 1000 m/s, that is, the method of the present disclosure identifies the structure with a lateral position of 2000 m and a longitudinal position of 450 m, while the conventional method identifies the structure with a lateral position of 2000 m and a longitudinal position of 430 m. Combined with the actual seismic situation that the structure at a lateral position of 2000 m has a longitudinal depth of 500 m, it indicates the location obtained by the characterization method of the present disclosure is closer to the actual situation than that identified by the conventional method.

Embodiment 2

In Embodiment 2, seismic signals from a work zone as shown in FIGS. 5a-5b for the marmousi2 model are selected, with a time of 0-1.2 s, a channel number of 1-3500, and a residual threshold of 0.001. The extracted velocity wavelet and acceleration wavelet are shown in FIGS. 6a-6b, and the mixed-phase wavelet reconstructed using the bi-spectral reconstruction of the higher-order cumulants is shown in FIG. 8.

The mixed-phase wavelet extracted by the method proposed in Embodiment 2 of the present disclosure (shown in FIG. 8), the wavelet obtained by the conventional method and the original seismic acceleration signal as shown in FIG. 5b are subjected to deconvolution to obtain fine characterization results. To show the relevant details, the fine characterization results are enlarged, as shown in FIGS. 13a-13b, where a (a): fine characterization results acquired by the depiction method of the present disclosure; and (b) fine characterization results acquired by the conventional method. The results shown in FIG. 12 are amplified display results of the original signal. By comparing fine characterization results obtained by different methods in FIGS. 13a-13b, with a channel spacing of 10 m (i.e., the distance between each two channels is 10 m) and a speed of 500 m/s, the pinch-out point location identified by the method proposed by the present disclosure is at 1600^th channel at 1.80 s, which corresponds to the actual situation of the position with a lateral position of 1600 m and a longitudinal depth of the 900 m. However, the conventional method can only identify the pinch-out point location at 1560^th channel at 1.84 s, which corresponds to the actual situation of the position with a lateral position of 1560 m and a longitudinal depth of 920 m. Considering that the position of the actual pinch-out point is at a lateral position of 1650 m and a longitudinal depth of 860 m. Therefore, the fine characterization results obtained by the method proposed in the present disclosure are closer to the actual location of the pinch-out point than the results obtained by the conventional method.

In the wavelet extraction on the above Marmousi2 model, instead of using the combined wavelet extraction method of the present disclosure, the velocity wavelet and the acceleration wavelet are extracted separately by using the existing method with all other parameters (e.g., number of processing channels, time window range) being the same, and the results are shown in FIGS. 7a-7b. The mixed-phase wavelet results extracted by using the existing method (without filter smoothing of third-order cumulative quantities) are shown in FIG. 9.

By comparing FIGS. 6a-6b and FIGS. 7a-7b, it can be seen that the velocity wavelet extracted separately by the existing method is basically the same as the velocity wavelet obtained by the joint extraction of the present disclosure, but the acceleration wavelet extracted separately is weaker than the acceleration wavelet obtained by the joint extraction of the present disclosure in terms of amplitude strength and other aspects.

By comparing FIGS. 8 and 9, it can be seen that the mixed-phase wavelet obtained by the joint extraction method of the present disclosure is closer to the theoretical wavelet than the mixed-phase wavelet separate extracted by the existing method, and has a higher accuracy.

In summary, the present disclosure enables the extraction of theoretical wavelets such as constant-phase wavelet and minimum phase wavelet, and can also extract the mixed-phase wavelet that has the characteristics of attenuated wavelets. The two extraction methods both involve the combination of the seismic velocity signal and the acceleration signal, allowing for improved extraction accuracy. Compared with the prior art, the present disclosure has achieved significant improvement.

The above are only preferred embodiments of the present disclosure, and are not intended to limit the present disclosure in any forms. Though the present disclosure has been illustrated in detail above, one of ordinary skill in the art can still make various changes or modifications to the technical features recited in the embodiments. It should be understood that those modifications or changes made without departing from the spirit of the present disclosure shall still fall within the scope of the present disclosure defined by the appended claims.

Claims

1. A geological structure characterization method based on a seismic velocity signal and a seismic acceleration signal, comprising: acquiring the seismic acceleration signal and the seismic velocity signal;normalizing the seismic acceleration signal and the seismic velocity signal;based on wavelet phase requirements, setting parameters to carry out constant-phase wavelet and minimum-phase wavelet extraction or to extract a mixed-phase wavelet; andwherein the constant-phase wavelet and minimum-phase wavelet extraction is performed through the following steps: (S1) extracting a first minimum-phase wavelet based on a normalized seismic velocity signal;(S2) extracting a constant-phase wavelet based on a normalized seismic acceleration signal; and converting the constant-phase wavelet to obtain a second minimum-phase wavelet; and(S3) calculating a residual between the first minimum-phase wavelet and the second minimum-phase wavelet;if the residual is greater than a preset threshold, returning to step (S2) and performing constant-phase wavelet extraction again; andif the residual is less than or equal to the preset threshold, outputting the first minimum-phase wavelet and the constant-phase wavelet;wherein the geological structure characterization method further comprises:after the constant-phase wavelet and minimum-phase wavelet extraction, converting the normalized seismic velocity signal in an acceleration domain followed by splicing with the normalized seismic acceleration signal to obtain fused data; andcalculating a third-order cumulant of the fused data and a third-order spectrum of the fused data; filtering the third-order cumulant using Parzen window; and calculating an amplitude spectrum and a phase spectrum of the constant-phase wavelet and the first minimum-phase wavelet using bispectrum analysis to reconstruct the mixed-phase wavelet; the mixed-phase wavelet is extracted through the following formulas:

2. The geological structure characterization method of claim 1, wherein in step (S1), the first minimum-phase wavelet is extracted based on the normalized seismic velocity signal through the following formulas:

3. The geological structure characterization method of claim 1, wherein in step (S2), the constant-phase wavelet is extracted based on the normalized seismic acceleration signal through the following formulas:

4. The geological structure characterization method of claim 3, wherein the phase correction factor is expressed as:

5. The geological structure characterization method of claim 1, wherein the wavelet phase spectrum Φ′ is calculated through the following formulas:

6. The geological structure characterization method of claim 1, wherein the wavelet amplitude spectrum |W| is calculated through the following formulas: