CN112462418A - Seismic signal processing method based on pseudo WVD and CZT - Google Patents

Seismic signal processing method based on pseudo WVD and CZT Download PDF

Info

Publication number
CN112462418A
CN112462418A CN202011221950.9A CN202011221950A CN112462418A CN 112462418 A CN112462418 A CN 112462418A CN 202011221950 A CN202011221950 A CN 202011221950A CN 112462418 A CN112462418 A CN 112462418A
Authority
CN
China
Prior art keywords
signal
czt
time
frequency
seismic
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202011221950.9A
Other languages
Chinese (zh)
Inventor
徐天吉
李思源
秦正晔
胡光岷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202011221950.9A priority Critical patent/CN112462418A/en
Publication of CN112462418A publication Critical patent/CN112462418A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Complex Calculations (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a seismic signal processing method based on pseudo WVD and CZT, which comprises the following steps: s1, Hilbert transformation is carried out on the seismic signals to obtain frequency domain signals, and the seismic signals and the frequency domain signals are added to obtain analytic signals; s2, obtaining a PWVD time frequency distribution result; s3, obtaining a CZT time-frequency distribution result of the PWVD; s4, rearranging the PWVD time frequency distribution result according to frequency to obtain a high-resolution single-frequency profile of the seismic signal; and (4) performing inverse transformation on the CZT time frequency distribution result, dividing the inverse transformed first line by the signal phasor and taking the real part of the division result to obtain a reconstructed time domain signal. The method comprises the steps of firstly transforming an original signal to a frequency domain, then carrying out high-resolution processing in the frequency domain, and finally returning to a time domain to obtain a high-resolution time-frequency spectrogram and a reconstructed seismic signal. The reconstructed signal is smoother, and the energy distribution condition of the thin layer can be better reflected; the description capacity of the seismic signal to the details of the edge of the reservoir is improved.

Description

Seismic signal processing method based on pseudo WVD and CZT
Technical Field
The invention belongs to the field of geoscience, and particularly relates to a seismic signal processing method based on pseudo WVD and CZT.
Background
Seismic exploration techniques are the most common means of oil and gas exploration. Seismic exploration is a geophysical exploration method which utilizes the difference of elasticity and density of underground media and deduces the properties and forms of underground rock layers by observing and analyzing the propagation rule of seismic waves generated by artificial earthquake in the underground. With the continuous improvement of oil exploration, the key points of oil and gas exploration are changed into types of complex structures, hidden lithologic oil and gas reservoirs and the like, and the explored geological structure is more complex and has the characteristics of small trap scale, thin reservoir thickness and the like. The resolution of seismic signals is required to be higher and higher in the whole exploration process, and seismic waves are typical non-stationary signals, and the energy distribution of the seismic waves mainly exists in a low-frequency area; because the traditional analysis method like Fourier transform can only perform global transformation in a time domain or a frequency domain, the time-frequency domain localization analysis on the seismic signals cannot be performed, and a good effect cannot be obtained in processing non-stationary signals. Due to a series of reasons such as deep gas reservoir burial, thin reservoir, rapid seismic energy and frequency band attenuation and the like, exploration personnel are difficult to accurately identify exploration targets, and deepening of oil-gas exploration and discovery of new targets are limited to a great extent. The resolution of the seismic signals is a key factor for acquiring reservoir information in oil and gas exploration work, and has important significance for researching a thin-layer geological structure; therefore, it is important to reasonably and effectively improve the resolution of the seismic signal according to the characteristics of the seismic signal.
Currently, common seismic signal processing methods include short-time fourier transform (STFT), wavelet transform, Wigner-Ville distribution, Cohen-like distribution, and the like. Both short-time Fourier transform and wavelet transform are in linear distribution, the distribution relation between the time domain and the frequency domain of seismic signals cannot be completely expressed, and the defects are generated to a certain extent when the seismic signals are analyzed. The Wigner-Ville distribution (WVD) and Cohen type distribution belong to bilinear distribution, and the Wigner-Ville distribution can visually describe the distribution of the energy of a signal in time domain and frequency domain, has excellent time-frequency focusing capability, but also has fatal defects: a large number of cross-interference terms are generated, and the underground reservoir structure is not easy to identify by researchers. A pseudo-Wigner-Ville distribution (PWVD) belongs to one of the Cohen-like distributions, which is based on the Wigner-Ville distribution by changing the kernel function by adding a moving window function, thereby suppressing the generation of cross terms to a large extent. Although the pseudo-Wigner-Ville distribution achieves a better time-frequency distribution result by suppressing the generation of cross terms, it cannot further improve the resolution of the seismic signal.
Chirp-Z transform (CZT) is a signal analysis method proposed by Lawrence Rabiner, originally used to analyze speech signals. The CZT is sampled at equal intervals on a unit circle of a Z plane by using a spiral line, sampling points are distributed at equal angles on the spiral line, and spectral analysis of signals can be realized on the spiral line on the Z plane. The method has the advantages of being effective in calculating the time-frequency distribution of the power spectrum, fast in calculation, free of limitation of the length of a data sequence and the like, and is very suitable for processing complex seismic signals. By adjusting the sampling interval and the number of sampling points, the time-frequency spectrogram can be refined, and the resolution of signals is improved.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provides a seismic signal processing method based on pseudo WVD and CZT, wherein the method comprises the steps of firstly converting an original signal into a frequency domain, then carrying out high-resolution processing in the frequency domain, and finally returning to a time domain to obtain a high-resolution time-frequency spectrogram and a reconstructed seismic signal, and can improve the description capacity of the seismic signal on reservoir edge details.
The purpose of the invention is realized by the following technical scheme: the seismic signal processing method based on the pseudo WVD and the CZT comprises the following steps:
s1, carrying out Hilbert transformation on the seismic signals to obtain frequency domain signals of the seismic signals, and adding the seismic signals and the frequency domain signals to obtain analytic signals;
s2, taking the time difference and the time sequence as independent variables to obtain the instantaneous autocorrelation function of the analysis signal; multiplying the instantaneous autocorrelation function by a rectangular window with the length of 2L-1 to obtain a kernel function of pseudo Wigner-Ville distribution; performing cyclic displacement on the kernel function to obtain a new kernel function, and performing Fourier transform on the new kernel function after the cyclic displacement to obtain a PWVD time-frequency distribution result, wherein 2L is a time sampling point of the seismic time-domain signal;
s3, multiplying the new kernel function after the cyclic shift by the weighting coefficient to obtain a sequence f*(l) The sequence f*(l) Circularly shifting to obtain a new sequence f (l) and constructing a unit response sequence g (l); then, respectively carrying out Fourier transform on a new sequence f (l) obtained by cyclic displacement and a unit response sequence g (l), multiplying the Fourier transform result in a frequency domain to obtain Y (r), further carrying out inverse transformation to obtain y (l), and multiplying the y (l) by a weighting coefficient to obtain a CZT time-frequency distribution result of the PWVD;
s4, rearranging the PWVD time frequency distribution result according to frequency to obtain a high-resolution single-frequency profile of the seismic signal; and performing inverse transformation on the CZT time-frequency distribution result, then dividing the first line of the inversely transformed data by the signal phasor and taking the real part of the division result to obtain a reconstructed time-domain signal.
Further, the step S1 includes the following sub-steps:
s11, acquiring a seismic section, and acquiring a single seismic signal x (n) by using the seismic section;
s12, carrying out Hilbert transform on a single-channel seismic signal x (n) to obtain a frequency domain signal H [ x (n) ], and then adding the frequency domain signal H [ x (n) ], the time domain seismic signal to obtain an analytic signal:
Figure BDA0002762351050000021
z(n)＝x(n)+jH[x(n)] (2)。
further, the step S2 includes the following sub-steps:
s21, with the time difference L and the time series n as arguments, wherein L ═ L +1, … 0,1, … L; n is 0,1, … 2L-1; determining the instantaneous autocorrelation function r (n, l) of the analytic signal z (n):
r(n,l)＝z(n+l)z(n-l)* (3)
wherein z (n + l) is an analysis signal at the time of n + l, and z (n-l)*Resolving the conjugate of the signal for n-l time;
s22, multiplying the instantaneous autocorrelation function r (n, L) by a rectangular window h (L) of length 2L-1, and when | L | > L, h (L) is 0; after being windowedDiscrete instantaneous autocorrelation function of Rz(n,l)
Rz(n,l)＝z(n+l)z(n-l)*h(l)h(-l)* (4)
In the formula, h (-l)*Is the conjugate of the window function h (-l); h (-l) replacing l in h (l) with-l;
s23, let Z (n, l) ═ Z (n + l) h (l), obtain kernel function K of pseudo-Wigner-Ville distributionn(l)：
Figure BDA0002762351050000031
In the formula, Z (n, -l)*Is the conjugate of Z (n, -l); the range of L is [ -min (n,2L-n), min (n,2L-n)]；
S24 Kernel function Kn(l) Making cyclic shift to make L be positioned in 0-2L-1 to obtain rearranged kernel function KN(n,l)：
Figure BDA0002762351050000032
In formula (6): z (n, -L +2L)*Is the conjugate of Z (n, -L + 2L);
s25 Kernel function KNFourier transform is carried out on (n, l) to obtain a discrete PWVD time frequency distribution resultz(n,σ)：
Figure BDA0002762351050000033
In formula (7): e.g. of the type-j2lσσ ═ n/2L is the twiddle factor of the fourier transform.
Further, the step S3 includes the following sub-steps:
s31, the kernel function K of the PWVD obtained in the step S2NMultiplying (n, l) by a weighting factor
Figure BDA0002762351050000034
Obtaining the sequence f*(l)：
Figure BDA0002762351050000035
In the formula, A-lIs the length of the sample amplitude of the CZT,
Figure BDA0002762351050000041
is a sample rotation factor of CZT;
s32, pairing sequence f*(l) Circularly shifting to obtain a new sequence f (l) after the displacement, and constructing a unit response sequence g (l):
Figure BDA0002762351050000042
Figure BDA0002762351050000043
in the expressions (9) and (10), Q is a sequence length constructed for CZT calculation, and is a quadratic power closest to N + M-1, N is a time sequence length of a discrete signal, and N is 2L; m is the number of sampling points of CZT;
s33, Fourier transform is carried out on g (l), f (l), and G (r), F (r) obtained by Fourier transform are multiplied in the frequency domain to obtain Y (r):
G(r)＝fft[g(l)],F(r)＝fft[f(l)] (11)
Y(r)＝F(r)G(r) (12)；
s34, carrying out inverse transformation on the Y (r) to obtain y (l); finally multiplying the coefficient to y (l) to obtain CZT distribution result CZT (z) of PWVDl)：
y(l)＝ifft[Y(r)] (13)
Figure BDA0002762351050000044
In the formula (I), the compound is shown in the specification,
Figure BDA0002762351050000045
is CZT weighting coefficient; because the repetition period in the PWVD frequency domain is pi, the number of sampling points of CZT is twice of the selected frequency band; in order to avoid the generation of the frequency spectrum aliasing phenomenon, the sampling factor of CZT is also increased by two times; the final CZT distribution should range from 0 to M/4.
Further, the step S4 includes the following sub-steps:
s41, the PWVD time frequency distribution result obtained in the step S2 is PWVDz(n, sigma) rearranging according to the frequency according to the channel number of the signal in sequence to obtain a single-frequency profile with high resolution for seismic interpretation;
s42, the CZT time frequency distribution result CZT (z) obtained in the step S3l) Inverse transformation to ifKN(n,l)：
ifKN(n,l)＝ifft[CZT(n,l)] (15)
In the formula, ifft is inverse Fourier transform;
s43, inverse transformation result ifKNDividing (n, l) with the signal vector and taking the real part of the signal vector to obtain a reconstructed time domain signal X after the refined frequency spectrum inverse transformation*(n)：
X*(n)＝re[ifKN(1,l)/x(n)] (16)
Where x (n) is the original time domain signal and re represents the real part of the signal.
The invention has the beneficial effects that: the method processes seismic signals based on pseudo Wigner-Ville distribution and Chrip-Z transformation, firstly transforms original signals to a frequency domain, then carries out high-resolution processing in the frequency domain, and finally returns to a time domain to obtain a high-resolution time-frequency spectrogram and reconstructed seismic signals. The CZT and the PWVD are combined for the first time and are used for analyzing the seismic signals; on the basis of PWVD, the seismic signal time-frequency spectrogram is refined, and the signal resolution is further improved. Compared with the original time domain signal, the reconstructed time domain signal obtained by the CZT inverse transformation is smoother, and can better reflect the energy distribution condition of the thin layer; meanwhile, the interlayer relation of the stratum is more obvious, the description capacity of the seismic signal to the details of the edge of the reservoir is improved, the oil-gas mineral deposit can be identified more conveniently, and more accurate seismic data are provided for high-precision identification of the underground thin interbed medium.
Drawings
FIG. 1 is a flow chart of a pseudo WVD and CZT based seismic signal processing method of the present invention;
FIG. 2 is a schematic diagram of a transformation path of CZT;
FIG. 3 is a schematic diagram of a CZT calculation process;
FIG. 4 is a time domain waveform diagram of an emulated signal;
FIGS. 5(a), (b), (c), (d), (e), and (f) are respectively a short-time Fourier transform result diagram, a three-dimensional diagram of a short-time Fourier transform result, a pseudo Wigner-Ville distribution result diagram, a three-dimensional diagram of a pseudo Wigner-Ville distribution result diagram, a Chrip-Z transform result diagram, and a three-dimensional diagram of a Chrip-Z transform result of the simulation signal of FIG. 4;
FIG. 6 is a reconstructed signal obtained by inverse transformation of CZT after the simulation signal is processed by a CZT algorithm based on PWVD;
FIG. 7 is a time domain waveform, STFT, PWVD, and CZT distribution plot of a numerically simulated seismic signal;
FIGS. 8(a) and (b) are a cross-sectional view of a digital analog signal PWVD and a cross-sectional view of CZT, respectively;
FIGS. 9(a) and (b) are respectively an original time domain signal diagram of a numerical analog signal and a reconstructed signal diagram obtained by CZT inverse transformation after being processed by a CZT algorithm based on PWVD;
FIG. 10 is a short-time Fourier transform profile of Crossline data;
FIG. 11 is a single frequency cross section of Crossline data processed by a CZT algorithm based on PWVD;
FIG. 12 is a short-time Fourier transform profile of Inline data;
FIG. 13 is a single frequency cross section of Inline data processed by a PWVD-based CZT algorithm;
fig. 14(a), (b), (c), and (d) are time domain signals of original Crossline data, Crossline time domain signals after CZT inverse transformation, time domain signals of original Inline data, and Inline time domain signals after CZT inverse transformation, respectively.
Detailed Description
The technical scheme of the invention is further explained by combining the attached drawings.
As shown in fig. 1, the seismic signal processing method based on pseudo WVD and CZT of the present invention includes the following steps:
s1, Hilbert transforming the seismic signals x (n) to obtain frequency domain signals H [ x (n) ], and adding the seismic signals x (n) and the frequency domain signals H [ x (n) ], so as to obtain analytic signals z (n); the method specifically comprises the following substeps:
s11, acquiring a seismic section, and acquiring a single seismic signal x (n) by using the seismic section;
s12, carrying out Hilbert transform on a single channel seismic signal x (n) to obtain a frequency domain signal H [ x (n) ], and then adding the frequency domain signal H [ x (n) ], the time domain seismic signal and the time domain seismic signal to obtain an analytic signal containing time domain and frequency domain information:
Figure BDA0002762351050000061
z(n)＝x(n)+jH[x(n)] (2)。
s2, using the time difference l and the time series n as independent variables, obtaining the instantaneous autocorrelation function r (n, l) of the analysis signal z (n); multiplying the instantaneous autocorrelation function r (n, L) by a rectangular window h (L) with the length of 2L-1 to obtain a kernel function K of pseudo Wigner-Ville distributionn(l) (ii) a Kernel function Kn(l) Performing cyclic shift to obtain new kernel function KN(n, l) such that the computational domain of the Fourier transform is non-negative. For new kernel function K after cyclic shiftN(n, l) Fourier transform to obtain PWVD time frequency distribution resultz(n, σ), wherein 2L is a time sampling point of the seismic time domain signal.
The method specifically comprises the following substeps:
s21, with the time difference L and the time series n as arguments, wherein L ═ L +1, … 0,1, … L; n is 0,1, … 2L-1; determining the instantaneous autocorrelation function r (n, l) of the analytic signal z (n):
r(n,l)＝z(n+l)z(n-l)* (3)
wherein z (n + l) is an analysis signal at the time of n + l, and z (n-l)*Resolving the conjugate of the signal for n-l time;
s22, multiplying the instantaneous autocorrelation function r (n, L) by a rectangular window h (L) of length 2L-1, and when | L | > L, h (L) is 0; obtaining a windowed discrete instantaneous autocorrelation function as Rz(n,l)
Rz(n,l)＝z(n+l)z(n-l)*h(l)h(-l)* (4)
In the formula, h (-l)*Is the conjugate of the window function h (-l); h (-l) replacing l in h (l) with-l;
s23, let Z (n, l) ═ Z (n + l) h (l), obtain kernel function K of pseudo-Wigner-Ville distributionn(l)：
Figure BDA0002762351050000071
In the formula, Z (n, -l)*Is the conjugate of Z (n, -l); the range of L is [ -min (n,2L-n), min (n,2L-n)]；
S24, since the computation domain of the fourier transform (FFT) is a non-negative number (L ═ 0,1,2 · · · 2L-1), the kernel function K needs to be paired withn(l) Making cyclic shift to make L be positioned in 0-2L-1 to obtain rearranged kernel function KN(n,l)：
Figure BDA0002762351050000072
In formula (6): z (n, -L +2L)*Is the conjugate of Z (n, -L + 2L);
s25 Kernel function KNFourier transform is carried out on (n, l) to obtain a discrete PWVD time frequency distribution resultz(n,σ)：
Figure BDA0002762351050000073
In formula (7): e.g. of the type-j2lσIs a rotation factor of Fourier transform, and is pi n/2L; since the period in the frequency domain of WVD is pi, not2 pi, the twiddle factor should also be increased by a factor of two.
S3, shifting the loop to a new kernel function KNMultiplying (n, l) by a weighting factor
Figure BDA0002762351050000074
Obtaining the sequence f*(l) The sequence f*(l) Circularly shifting to obtain a new sequence f (l) and constructing a unit response sequence g (l); then respectively carrying out Fourier transform on a new sequence f (l) obtained by cyclic displacement and a unit response sequence g (l), then multiplying Fourier transform results F (r), G (r) in a frequency domain to obtain Y (r), further carrying out inverse transformation on Y (r) to obtain y (l), and multiplying y (l) by a weighting coefficient to obtain a CZT time-frequency distribution result CZT (z) of the PWVDl)；
The method specifically comprises the following substeps:
s31, the kernel function K of the PWVD obtained in the step S2NMultiplying (n, l) by a weighting factor
Figure BDA0002762351050000081
Obtaining the sequence f*(l)：
Figure BDA0002762351050000082
In the formula, A-lIs the length of the sample amplitude of the CZT,
Figure BDA0002762351050000083
is a sample rotation factor of the CZT,
Figure BDA0002762351050000084
W0is the length of the sample
Figure BDA0002762351050000085
Is the sampling interval;
s32, K is derived from the derivation formula of PWVDN(n, l) changing original negative number domain data into positive number domain data through cyclic displacement, wherein data sections between the new positive number domain data and the original positive number domain data are all zero; therefore, in the CZT algorithm, tooZero padding is carried out in the middle of the data segment, and the consistency of the convolution interval is ensured. For sequence f*(l) Circularly shifting to obtain a new sequence f (l) after the displacement, and constructing a unit response sequence g (l):
Figure BDA0002762351050000086
Figure BDA0002762351050000087
in the expressions (9) and (10), Q is a sequence length constructed for CZT calculation, and is a quadratic power closest to N + M-1, N is a time sequence length of a discrete signal, and N is 2L; m is the number of sampling points of CZT;
s33, Fourier transform is carried out on g (l), f (l), and G (r), F (r) obtained by Fourier transform are multiplied in the frequency domain to obtain Y (r):
G(r)＝fft[g(l)],F(r)＝fft[f(l)] (11)
Y(r)＝F(r)G(r) (12)；
s34, carrying out inverse transformation on the Y (r) to obtain y (l); finally multiplying the coefficient to y (l) to obtain CZT distribution result CZT (z) of PWVDl)：
y(l)＝ifft[Y(r)] (13)
Figure BDA0002762351050000091
In the formula (I), the compound is shown in the specification,
Figure BDA0002762351050000092
is CZT weighting coefficient; because the repetition period in the PWVD frequency domain is pi, the number of sampling points of CZT is twice of the selected frequency band; in order to avoid the generation of the frequency spectrum aliasing phenomenon, the sampling factor of CZT is also increased by two times; the final CZT distribution should range from 0 to M/4.
The principle of CZT involved in the step is as follows:
the Z transformation of a discrete signal sequence of known finite length x (n) into X (Z) includes
Figure BDA0002762351050000093
Let Zr＝AW-r，
Figure BDA0002762351050000094
Then there are:
Figure BDA0002762351050000095
in the formula A0，W0Given A as any positive real number0，W0，θ0，
Figure BDA0002762351050000096
When r is 0,1, …, ∞, a point z on the z-plane can be obtained0,z1,…,z∞Taking the Z transform at these points:
Figure BDA0002762351050000097
as can be seen from formula (16), when r is 0,
Figure BDA0002762351050000098
the amplitude of the point in the z-plane is A0Argument is theta0Is the starting point of CZT; when r is equal to 1, the compound is,
Figure BDA0002762351050000099
z1amplitude of the dot becomes
Figure BDA00027623510500000910
Angle at theta0On the basis of (1) has an increment
Figure BDA00027623510500000911
It is not difficult to imagine that point z varies with r0,z1,…zM-1Constituting the path of the CZT transform.
The CZT conversion path is shown in fig. 2.
When performing spectrum analysis on a signal, CZT should be realized on a unit circle. A. the0And W0All should be taken to be 1, let the length of the discrete signal x (N) be 0,1, …, N-1, and the length of the transform r be 0,1, …, M-1, i.e.
Figure BDA00027623510500000912
From Bluestein formula:
Figure BDA00027623510500000913
equation (18) can in turn be written as:
Figure BDA0002762351050000101
order to
Figure BDA0002762351050000102
Can obtain the product
Figure BDA0002762351050000103
Equation (20) is the definition of CZT, and it can be seen that CZT is obtained by first multiplying a discrete signal by a coefficient
Figure BDA0002762351050000104
Carrying out weighting treatment once to obtain f (n); then, by a linear system with unit response g (n); finally, the convolution result is multiplied by a coefficient
Figure BDA0002762351050000105
Once again, weighting achieves helical sampling in three-dimensional space. FIG. 3 is a calculation process of CZT.
S4, PWVD time frequency distribution result PWVDz(n, σ) rearranging according to the frequency to obtain a high-resolution single-frequency profile of the seismic signal; then, the CZT time frequency distribution result CZT (z)l) Inverse transform is performed, and then inverse transformed data ifK is takenNThe first row of (n, l) is divided by the signal phasor and the real part of the division result is taken to obtain the reconstructed time domain signal X*(n)。
The method specifically comprises the following substeps:
s41, the PWVD time frequency distribution result obtained in the step S2 is PWVDz(n, sigma) rearranging according to the frequency according to the channel number of the signal in sequence to obtain a single-frequency profile with high resolution for seismic interpretation;
s42, the CZT time frequency distribution result CZT (z) obtained in the step S3l) Inverse transformation to ifKN(n,l)：
ifKN(n,l)＝ifft[CZT(n,l)] (21)
In the formula, ifft is inverse Fourier transform;
s43, inverse transformation result ifKN(n, l) is similar to the instantaneous autocorrelation matrix K in WVDN(n, l), and the first row (K) of the instantaneous autocorrelation matrixN(1,1),KN(1,2)……KN(1, N)), both are the time instants of the signal sequence multiplied by the conjugate matrix of the signal itself at the same time instant. Therefore, the result ifK of inverse transformation is obtainedNDividing (n, l) with the signal vector and taking the real part of the signal vector to obtain a reconstructed time domain signal X after the refined frequency spectrum inverse transformation*(n)：
X*(n)＝re[ifKN(1,l)/x(n)] (22)
Where x (n) is the original time domain signal and re represents the real part of the signal.
The signal processing effect of the present invention is verified by specific signals.
Example 1 simulation Signal
To prove the feasibility and effectiveness of the fractional algorithm for seismic signal processing based on pseudo Wigner-Ville distribution and Chrip-Z transformation. Firstly, adopting the following simulation signals, wherein the peak values of the signals are at 0.2s, 0.6s and 1s, and the corresponding frequency components are 40Hz, 20Hz and 10Hz respectively, which is the result of simulating a single-channel seismic signal; the sampling frequency of the signal is 1000Hz, and the number of sampling points is 1200. The time domain waveform of the signal is shown in fig. 4.
Figure BDA0002762351050000111
FIG. 5(a) is the result of a short-time Fourier transform of a simulated signal; FIG. 5(b) is a three-dimensional plot of the results of a short-time Fourier transform; FIG. 5(c) shows the pseudo Wigner-Ville distribution of the simulated signal; FIG. 5(d) is a three-dimensional plot of the results of a pseudo Wigner-Ville distribution; FIG. 5(e) is the result of the Chrip-Z transformation of the simulated signal; FIG. 5(f) is a three-dimensional plot of the results of the Chrip-Z transformation. As can be seen from the time-frequency spectrogram after simulation signal processing, the PWVD added with the window function not only has excellent time-frequency focusing performance, but also keeps the mathematical characteristics of the WVD, better inhibits the generation of interference terms, and obtains a better time-frequency analysis result. By combining CZT and PWVD and setting different argument parameters by CZT, not only is the time-frequency spectrogram of the simulation signal refined, but also the time-frequency focusing capacity of the simulation signal is further improved, and therefore the time-frequency resolution of the signal is improved.
Fig. 6 shows a reconstructed signal obtained by inverse transformation of CZT. It can be seen that the time domain waveform of the original signal is substantially restored by the inverse transformation of CZT; because the simulation signal is simple and smooth and free of noise, and the amplitude and phase information of the signal after the signal passes through the WVD and the CZT are changed, the peak values of the reconstructed signal near 0.2s and 0.6s are slightly different from those of the original signal.
Example 2 forward modeling of a numerical analog signal
After the simulated signals verify the effectiveness of the algorithm, the algorithm is then applied to the numerically simulated seismic signals. The adopted theoretical seismic signals are stored in a two-dimensional CDP mode, the theoretical seismic signals totally comprise 101 channels, the sampling period is 1ms, each row of channels are provided with 231 time sampling points, and the sampling frequency is 1000 Hz.
Fig. 7(a) - (c) are time domain waveforms of the channel in column 1, 60 and 101, respectively; fig. 7(d) to (f) are the STFT of the 1 st column channel, the STFT of the 60 th column channel, and the STFT of the 101 st column channel, respectively; (g) (ii) PWVD of the 1 st column channel, PWVD of the 60 th column channel, and PWVD of the 101 st column channel, respectively; (j) the CZT distribution of the 1 st channel, the CZT distribution of the 60 th channel and the CZT distribution of the 101 th channel are respectively shown in the (l) to (l). By processing the time-frequency spectrogram obtained by the analog signal of the forward modeling value, on the basis of PWVD of a single-channel seismic signal, CZT carries out interpolation processing on time-frequency distribution in a sampling mode by setting a sampling frequency band and refining points on a Z plane through spiral sampling, and a corresponding sampling interval is set, so that the resolution of the original PWVD is further improved, and the effect of frequency spectrum refining is achieved.
FIG. 8(a) is a cross-sectional view of a numerical analog signal PWVD; fig. 8(b) is a cross-sectional view of the numerical analog signal CZT. Comparing the two steps, it can be seen that the technical solutions of steps S2 and S3 in the present invention realize the spectrum refinement of the profile, and improve the time-frequency resolution of the signal.
FIG. 9(a) is an original time domain signal of a digital analog signal, which contains 101 channels of signals; fig. 9(b) shows a reconstructed signal obtained by inverse CZT transformation. The comparison shows that through the technical scheme of the step S4 in the embodiment of the invention, the reconstructed signal obtained by CZT is smoother, not only has better resolution and energy focusing, but also the local edge details of the reservoir trend are clearer.
Example 3 actual seismic data
Based on the theoretical model research, the method is used for actual seismic data of 6142.0-6152.4 meters of an upper reservoir and a lower reservoir located at 6171.8-6228.0 meters of a well 1 of a test area with a western depression in the west of the Sichuan basin. The data contains Crossline data in the horizontal direction and Inline data in the vertical direction, and the total number of 1011 channels and the number of time sampling points is 1001.
FIG. 10 is a cross-sectional view of a short-time Fourier transform of Crossline; FIGS. 11(a) to (f) are single-frequency cross-sectional diagrams of Crossline data processed by CZT algorithm time-frequency based on PWVD, and the single-frequency cross-sectional diagrams correspond to frequencies of 35Hz, 45Hz, 55Hz, 65Hz, 75Hz, and 85Hz, respectively; fig. 12 is a short-time fourier transform cross-section of Inline. Fig. 13(a) to (f) are single-frequency cross-sectional views of Inline data processed by the CZT algorithm based on PWVD, and correspond to frequencies of 35Hz, 45Hz, 55Hz, 65Hz, 75Hz, and 85Hz, respectively. The comparison shows that after the PWVD and the CZT, the energy distribution of the signals is more obvious, the frequency spectrum bands at different heights are obviously thinned, and the time-frequency resolution is improved, so that the frequency spectrum refinement of Crossline and Inline data is realized. The refined time-frequency spectrogram can better distinguish the geological structure of the horizontal trend of the reservoir, distinguish the edge details of the reservoir and be more beneficial to recognizing the underground thin interbed medium.
Fig. 14(a) is a time domain signal of original Crossline data; FIG. 14(b) is a Crossline time domain signal after CZT inverse transformation; fig. 14(c) is a time domain signal of the original Inline data; fig. 14(d) shows the Inline time domain signal after the CZT inverse transform. Compared with the original signal, the inverse-transformed time domain signal not only increases the edge details of the reservoir, so that the intervals between the stratums with different levels are more obvious; meanwhile, areas with similar energy distribution of the reservoir are distinguished, and the reservoir identification is facilitated. And the time-frequency graph after thinning is inversely transformed, so that the precision of the time-domain signal is improved.
It will be appreciated by those of ordinary skill in the art that the embodiments described herein are intended to assist the reader in understanding the principles of the invention and are to be construed as being without limitation to such specifically recited embodiments and examples. Those skilled in the art can make various other specific changes and combinations based on the teachings of the present invention without departing from the spirit of the invention, and these changes and combinations are within the scope of the invention.

Claims (5)

1. The seismic signal processing method based on the pseudo WVD and the CZT is characterized by comprising the following steps of:
s1, carrying out Hilbert transformation on the seismic signals to obtain frequency domain signals of the seismic signals, and adding the seismic signals and the frequency domain signals to obtain analytic signals;
s2, taking the time difference and the time sequence as independent variables to obtain the instantaneous autocorrelation function of the analysis signal; multiplying the instantaneous autocorrelation function by a rectangular window with the length of 2L-1 to obtain a kernel function of pseudo Wigner-Ville distribution; performing cyclic displacement on the kernel function to obtain a new kernel function, and performing Fourier transform on the new kernel function after the cyclic displacement to obtain a PWVD time-frequency distribution result, wherein 2L is a time sampling point of the seismic time-domain signal;
s3, multiplying the new kernel function after the cyclic shift by the weighting coefficient to obtain a sequence f*(l) The sequence f*(l) Circularly shifting to obtain a new sequence f (l) and constructing a unit response sequence g (l); then, respectively carrying out Fourier transform on a new sequence f (l) obtained by cyclic displacement and a unit response sequence g (l), multiplying the Fourier transform result in a frequency domain to obtain Y (r), further carrying out inverse transformation to obtain y (l), and multiplying the y (l) by a weighting coefficient to obtain a CZT time-frequency distribution result of the PWVD;
s4, rearranging the PWVD time frequency distribution result according to frequency to obtain a high-resolution single-frequency profile of the seismic signal; and performing inverse transformation on the CZT time-frequency distribution result, then dividing the first line of the inversely transformed data by the signal phasor and taking the real part of the division result to obtain a reconstructed time-domain signal.
2. The pseudo-WVD and CZT based seismic signal resolution enhancement processing method according to claim 1, wherein the step S1 includes the sub-steps of:
s11, acquiring a seismic section, and acquiring a single seismic signal x (n) by using the seismic section;
s12, carrying out Hilbert transform on a single-channel seismic signal x (n) to obtain a frequency domain signal H [ x (n) ], and then adding the frequency domain signal H [ x (n) ], the time domain seismic signal to obtain an analytic signal:
Figure FDA0002762351040000011
z(n)＝x(n)+jH[x(n)] (2)。
3. the pseudo-WVD and CZT based seismic signal resolution enhancement processing method according to claim 1, wherein the step S2 includes the sub-steps of:
s21, with the time difference L and the time series n as arguments, wherein L ═ L +1, … 0,1, … L; n is 0,1, … 2L-1; determining the instantaneous autocorrelation function r (n, l) of the analytic signal z (n):
r(n,l)＝z(n+l)z(n-l)* (3)
wherein z (n + l) is an analysis signal at the time of n + l, and z (n-l)*Resolving the conjugate of the signal for n-l time;
s22, multiplying the instantaneous autocorrelation function r (n, L) by a rectangular window h (L) of length 2L-1, and when | L | > L, h (L) is 0; obtaining a windowed discrete instantaneous autocorrelation function as Rz(n,l)
Rz(n,l)＝z(n+l)z(n-l)*h(l)h(-l)* (4)
In the formula, h (-l)*Is the conjugate of the window function h (-l); h (-l) replacing l in h (l) with-l;
s23, let Z (n, l) ═ Z (n + l) h (l), obtain kernel function K of pseudo-Wigner-Ville distributionn(l)：
Figure FDA0002762351040000021
In the formula, Z (n, -l)*Is the conjugate of Z (n, -l); the range of L is [ -min (n,2L-n), min (n,2L-n)]；
S24 Kernel function Kn(l) Making cyclic shift to make L be positioned in 0-2L-1 to obtain rearranged kernel function KN(n,l)：
Figure FDA0002762351040000022
In formula (6): z (n, -L +2L)*Is the conjugate of Z (n, -L + 2L);
s25 Kernel function KNFourier transform is carried out on (n, l) to obtain a discrete PWVD time frequency distribution resultz(n,σ)：
Figure FDA0002762351040000023
In formula (7):
Figure FDA0002762351040000024
σ ═ n/2L is the twiddle factor of the fourier transform.
4. The pseudo WVD and CZT based seismic signal resolution enhancement processing method according to claim 2, wherein the step S3 includes the sub-steps of:
s31, the kernel function K of the PWVD obtained in the step S2NMultiplying (n, l) by a weighting factor
Figure FDA0002762351040000025
Obtaining the sequence f*(l)：
Figure FDA0002762351040000026
In the formula, A-lIs the length of the sample amplitude of the CZT,
Figure FDA0002762351040000027
is a sample rotation factor of CZT;
s32, pairing sequence f*(l) Circularly shifting to obtain a new sequence f (l) after the displacement, and constructing a unit response sequence g (l):
Figure FDA0002762351040000031
Figure FDA0002762351040000032
in the expressions (9) and (10), Q is a sequence length constructed for CZT calculation, and is a quadratic power closest to N + M-1, N is a time sequence length of a discrete signal, and N is 2L; m is the number of sampling points of CZT;
s33, Fourier transform is carried out on g (l), f (l), and G (r), F (r) obtained by Fourier transform are multiplied in the frequency domain to obtain Y (r):
G(r)＝fft[g(l)],F(r)＝fft[f(l)] (11)
Y(r)＝F(r)G(r) (12)；
s34, carrying out inverse transformation on the Y (r) to obtain y (l); finally multiplying the coefficient to y (l) to obtain CZT distribution result CZT (z) of PWVDl)：
y(l)＝ifft[Y(r)] (13)
Figure FDA0002762351040000033
5. The pseudo-WVD and CZT based seismic signal resolution enhancement processing method according to claim 4, wherein the step S4 includes the sub-steps of:
s41, the PWVD time frequency distribution result obtained in the step S2 is PWVDz(n, sigma) rearranging according to the frequency according to the channel number of the signal in sequence to obtain a single-frequency profile with high resolution for seismic interpretation;
s42, the CZT time frequency distribution result CZT (z) obtained in the step S3l) Inverse transformation to ifKN(n,l)：
ifKN(n,l)＝ifft[CZT(n,l)] (15)
In the formula, ifft is inverse Fourier transform;
s43, inverse transformation result ifKNDividing (n, l) with the signal vector and taking the real part of the signal vector to obtain a reconstructed time domain signal X after the refined frequency spectrum inverse transformation*(n)：
X*(n)＝re[ifKN(1,l)/x(n)] (16)
Where x (n) is the original time domain signal and re represents the real part of the signal.
CN202011221950.9A 2020-11-05 2020-11-05 Seismic signal processing method based on pseudo WVD and CZT Pending CN112462418A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011221950.9A CN112462418A (en) 2020-11-05 2020-11-05 Seismic signal processing method based on pseudo WVD and CZT

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011221950.9A CN112462418A (en) 2020-11-05 2020-11-05 Seismic signal processing method based on pseudo WVD and CZT

Publications (1)

Publication Number Publication Date
CN112462418A true CN112462418A (en) 2021-03-09

Family

ID=74826240

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011221950.9A Pending CN112462418A (en) 2020-11-05 2020-11-05 Seismic signal processing method based on pseudo WVD and CZT

Country Status (1)

Country Link
CN (1) CN112462418A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114188930A (en) * 2022-02-16 2022-03-15 国网江西省电力有限公司电力科学研究院 Ferromagnetic resonance overvoltage intelligent identification method and device

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102466819A (en) * 2010-11-03 2012-05-23 中国石油天然气集团公司 Spectrum analysis method of seismic signal and apparatus thereof
CN103293515A (en) * 2012-11-19 2013-09-11 西北工业大学 Ship and warship line spectrum noise source longitudinal distribution characteristic measuring method
US20140297188A1 (en) * 2013-03-29 2014-10-02 Cgg Services Sa Time-frequency representations of seismic traces using wigner-ville distributions
CN105388527A (en) * 2015-11-30 2016-03-09 中国石油大学(北京) Hydrocarbon detection method based on complex domain matching pursuit algorithm

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102466819A (en) * 2010-11-03 2012-05-23 中国石油天然气集团公司 Spectrum analysis method of seismic signal and apparatus thereof
CN103293515A (en) * 2012-11-19 2013-09-11 西北工业大学 Ship and warship line spectrum noise source longitudinal distribution characteristic measuring method
US20140297188A1 (en) * 2013-03-29 2014-10-02 Cgg Services Sa Time-frequency representations of seismic traces using wigner-ville distributions
CN105388527A (en) * 2015-11-30 2016-03-09 中国石油大学(北京) Hydrocarbon detection method based on complex domain matching pursuit algorithm

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
田丰 等: "窄带细化Wigner—Ville分布分析的快速实现方法", 《电子与信息学报》 *
田琳 等: "平滑伪Wigner-Ville分布在地震信号处理中的应用", 《新疆师范大学学报（自然科学版）》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114188930A (en) * 2022-02-16 2022-03-15 国网江西省电力有限公司电力科学研究院 Ferromagnetic resonance overvoltage intelligent identification method and device

Similar Documents

Publication Publication Date Title
Båth Spectral analysis in geophysics
Poiata et al. Multiband array detection and location of seismic sources recorded by dense seismic networks
Lacoss Data adaptive spectral analysis methods
Wang Multichannel matching pursuit for seismic trace decomposition
US11333783B2 (en) Integrated method for estimation of seismic wavelets and synthesis of seismic records in depth domain
CN107678062B (en) The integrated forecasting deconvolution of the domain hyperbolic Radon and feedback loop methodology multiple suppression model building method
Herrera et al. Applications of the synchrosqueezing transform in seismic time-frequency analysis
Pinnegar et al. The S-transform with windows of arbitrary and varying shape
CN102221708B (en) Fractional-Fourier-transform-based random noise suppression method
Baig et al. Denoising seismic noise cross correlations
CN107608935B (en) Time rearrangement compression transformation-based impact signal time-frequency analysis and reconstruction method
US20140297188A1 (en) Time-frequency representations of seismic traces using wigner-ville distributions
CN109375265B (en) Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit
WO2008152364A1 (en) Method of representing seismic signals
CN111797552A (en) Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum
Cai et al. Seismic data denoising based on mixed time-frequency methods
CN113608259B (en) Seismic thin layer detection method based on ICEEMDAN constraint generalized S transformation
Zhong et al. Simulation of seismic-prospecting random noise in the desert by a Brownian-motion-based parametric modeling algorithm
CN109459789A (en) Time-domain full waveform inversion method based on amplitude decaying and linear interpolation
CN112462418A (en) Seismic signal processing method based on pseudo WVD and CZT
Liu et al. Estimating correlations of neighbouring frequencies in ambient seismic noise
Tianji et al. A microscopic ancient river channel identification method based on maximum entropy principle and Wigner-Ville Distribution and its application
CN111474582A (en) Precise S transformation method for generating high-precision time frequency spectrum
Likkason Spectral analysis of geophysical data
Tribolet Applications of short-time homomorphic signal analysis to seismic wavelet estimation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination