CN110764147A - Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition - Google Patents

Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition Download PDF

Info

Publication number
CN110764147A
CN110764147A CN201911084413.1A CN201911084413A CN110764147A CN 110764147 A CN110764147 A CN 110764147A CN 201911084413 A CN201911084413 A CN 201911084413A CN 110764147 A CN110764147 A CN 110764147A
Authority
CN
China
Prior art keywords
data
sub
signal
noise
spectrum
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.)
Granted
Application number
CN201911084413.1A
Other languages
Chinese (zh)
Other versions
CN110764147B (en
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.)
Jilin University
Original Assignee
Jilin University
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 Jilin University filed Critical Jilin University
Priority to CN201911084413.1A priority Critical patent/CN110764147B/en
Publication of CN110764147A publication Critical patent/CN110764147A/en
Application granted granted Critical
Publication of CN110764147B publication Critical patent/CN110764147B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G01V2210/324—Filtering

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)
  • Radar Systems Or Details Thereof (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention relates to a desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition, and belongs to the technical field of geophysical. Transforming the noisy seismic data into a frequency-displacement domain by local windowed Fourier transform to obtain data blocks and F-X spectrums thereof; performing modal decomposition on the real part and the imaginary part of each line of data in the F-X spectrum along the displacement direction by utilizing a VMD algorithm, and performing summation reconstruction on the effective modes after respective decomposition to obtain the F-X spectrum after processing each data block; and carrying out Fourier inversion and windowing inversion on the processed F-X spectrum of each data block to obtain a recovery result. The method avoids the problem of partial frequency component loss caused by sub-mode function reconstruction, considers the dual factors of signal and noise forms and frequency, has better integrity of the recovered signal, has simple calculation process, higher speed and high accuracy, and has important significance for desert seismic data processing and subsequent parameter inversion problems.

Description

Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition
Technical Field
The invention belongs to the technical field of geophysical, in particular to a method for recovering weak signals in desert exploration,
background
As is well known, the signal-to-noise ratio and resolution of seismic data directly affect subsequent structural prediction and parametric inversion results. Therefore, how to improve the seismic data quality and reconstruct the weak seismic effective signals becomes the key point in the seismic data processing. The main factor influencing the recovery of weak effective signals in seismic exploration data in desert areas is the unique geological condition. Because the earth surface of the desert area is covered by sand dunes, the propagation energy of effective seismic waves is seriously absorbed and attenuated, the generated random noise has the characteristics of non-stability, non-Gaussian, non-linearity, low frequency and the like, and weak effective seismic signals are submerged in the noise and are difficult to identify and recover.
At present, there are many methods related to seismic exploration random noise suppression and effective signal recovery at home and abroad, and the methods include the following developed methods such as Curvelet transformation, Contourlet transformation, S transformation, time-frequency peak value filtering (TFPF) and the like besides the traditional median filtering, wiener filtering, adaptive filtering, wavelet transformation, Empirical Mode Decomposition (EMD), f-x filtering, Empirical Mode Decomposition (EMD) and the like. However, the recovery results of effective signals in desert exploration data are not ideal enough in the existing methods, so that the search of a new method suitable for recovering weak signals in desert seismic exploration data becomes a problem of close attention of researchers.
Disclosure of Invention
The invention provides a desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition, which aims to solve the problem of weak signal recovery under a complex noise background in a desert area.
The technical scheme adopted by the invention is that the method comprises the following steps:
1) based on the diversity of the same-phase axes in seismic exploration data and the time-varying property of signals, firstly, two-dimensional windowing localization processing is adopted for noisy seismic data, and the similarity of data characteristics in each window is ensured as much as possible; for a sheet n1×n2Contains n2Each seismic trace of n1Selecting a proper time window length n with the seismic channel direction as the displacement directiontAnd displacement window length nlSliding the two-dimensional window on the noisy seismic data, with a certain coincidence between two adjacent windows, using a coincidence coefficient ξ (0)<ξ<1) Means that the data D is divided into blocks, and the divided ith data block is used as D(i)(i is 1,2, …, N), where N is the number of divided data blocks:
D(i)＝[D(i)(a,b)]，a＝1～nt，b＝1～nl(1)
thus D(i)Column b in (1) for data D(i)(:,b)The a-th data is represented by D(i)(a,:)To represent;
then, for the noisy data block D(i)Each column of data (i.e., in the time direction) in (i-1, 2, …, N) is fourier transformed by D(i)Column b data D in (1)(i)(:,b)For example, equation (2) gives the equation for which the fourier transform is performed:
Figure BDA0002263123090000021
where k is the frequency domain sampling point number, n is the time domain sampling point number, j is the imaginary unit, F(i),(:,b)(k) As time domain data D(i),(:,b)(n) corresponding frequency domain data; d can be calculated by the formula (2)(i)In the frequency domain data corresponding to each column of data, then D(i)F-X spectrum F of(i)It can be obtained by combining the frequency domain data of each column of data, i.e. F(i)＝[F(i)(:,1)(k),F(i),(:,2)(k),…,F(i),(:,nl)(k)]So as to analogize the situation of other data blocks;
2) VMD is an adaptive and non-recursive decomposition algorithm that can be viewed as a constrained variational optimization problem, whose continuous form is shown in equation (3):
Figure BDA0002263123090000022
wherein f (t) is data to be decomposed, t is time variable, j is imaginary unit, um(t) is the mth sub-modal function, ω, obtained after decompositionmIs the center frequency of the mth sub-modal function, delta (t) is the dirac function,is um(t) a hilbert transform, operative to convert the hilbert transform into an analytic signal having a single-sided frequency spectrum;
applying a VMD algorithm to the F-X spectrum data of each data block obtained in the step 1) to realize the separation of effective signals and noise components in the F-X spectrum data, and utilizing the VMD algorithmThe row vectors (i.e. along the direction of displacement) of the F-X spectral data are decomposed, F(i)May be denoted as F(i)(a,:)(ii) a Due to F(i)Is complex, the direct application of VMD algorithm can not realize the decomposition of complex number, so we adopt the way of respectively decomposing and recombining real part and imaginary part to realize the decomposition process of F-X spectrum data, the formula (4) gives F(i)A row data F(i)(a,:)The complex form of (a):
F(i)(a,:)＝F(i)(a,:)real+jF(i)(a,:)img，a＝1～nt(4)
wherein j is an imaginary unit, F(i)(a,:)realAnd F(i)(a,:)imgRespectively real part data and imaginary part data;
f is decomposed by using VMD decomposition algorithm in formula (3)(i)(a,:)realAnd F(i)(a,:)imgDecomposing into multiple sub-mode functions, each sub-mode function after decomposition and F(i)(a,:)realAnd F(i)(a,:)imgThe relationship between them is:
Figure BDA0002263123090000024
Figure BDA0002263123090000025
where K is the number of sub-mode functions, usually between 4 and 7, and up,realAnd up,imgRespectively representing real part data F(i)(a,:)realAnd imaginary data F(i)(a,:)imgThe p-th sub-modal function of (1);
3) function of submode up,realAnd up,imgThe method comprises two main classes of a sub-mode function taking signal characteristics as a leading part and a sub-mode function taking noise as a leading part, signal components are mainly collected in the signal leading part sub-mode function, the noise is collected by the noise leading part sub-mode function, and in order to accurately judge the class of the sub-mode function, a difference factor epsilon is respectively defined by an equation (7) and an equation (8)realAnd εimgTo measure whether the sub-modal function is a signal-dominant or noise-dominant modal function;
Figure BDA0002263123090000031
Figure BDA0002263123090000032
where s (k) denotes the desired signal, noi (k) denotes noise, up,real(k) And up,img(k) P-th sub-modal functions, ε, representing real and imaginary data, respectivelyrealAnd εimgThe larger the value of the sub-mode function is, the larger the difference between the sub-mode function and the effective signal is, and when the difference value is larger than a certain threshold value, the sub-mode function is judged to be a sub-mode function taking noise as the leading factor; otherwise, the signal is judged to be a sub-mode function taking the signal as the leading factor;
assuming that the first M mode functions in each decomposed sub-mode function are signal-dominant sub-mode functions, then equations (5) and (6) can be rewritten as:
Figure BDA0002263123090000033
Figure BDA00022631230900000311
wherein u isp,real signalAnd up,img signalRepresenting signal-dominant sub-mode functions, up,real noiseAnd up,img noiseRepresenting a noise-dominant sub-mode function, wherein M represents the number of signal-dominant sub-mode functions, discarding the noise-dominant sub-mode functions, and summing the M signal-dominant sub-mode functions to obtain:
Figure BDA0002263123090000034
Figure BDA0002263123090000035
then will beAnd
Figure BDA0002263123090000037
combining into complex form to obtain the processed F-X spectrum as follows:
Figure BDA0002263123090000038
wherein
Figure BDA0002263123090000039
Is F(i)A row data F(i)(a,:)Processed F-X spectral data, F(i)The other rows of data are analogized, and F can be obtained after the data of each row are processed in the same way(i)F-X spectra after treatment
Figure BDA00022631230900000310
Finally, toThe data in each column in (1) is subjected to inverse Fourier transform, and the data in the b-th columnThe result of performing the inverse fourier transform is:
Figure BDA0002263123090000043
wherein k is the serial number of the sampling point in the frequency domain, n is the serial number of the sampling point in the time domain, j is the unit of imaginary number,
Figure BDA0002263123090000044
the data blocks are corresponding time domain sequences, namely the recovered effective signals, and the rest of the column data are analogized, so that the recovered multi-channel effective signal data blocks are obtained
Figure BDA0002263123090000045
Combining each effective signal data block to obtain a whole effective signal data matrix
Figure BDA0002263123090000046
Seismic recordings can be viewed as two-dimensional data containing different frequency, different directional components. In seismic recording, the effective signal is composed of a plurality of reflected waves of fixed apparent velocity, i.e., the signal is predictable, while the noise is unpredictable. And each frequency in an F-X spectrum obtained after Fourier transformation of each reflected wave with a single apparent velocity has predictability, namely frequency division processing can be carried out in an F-X domain.
The VMD algorithm is an adaptive and non-recursive decomposition algorithm, can adaptively decompose the seismic data, and utilizes a wiener filtering process embedded in the process to obtain decomposition with slight modal aliasing, and the sub-band of the decomposition is finer and more compact. The process can be seen as a constrained variational optimization problem described in the following equation:
Figure BDA0002263123090000047
wherein f (t) is data to be decomposed, t is time variable, j is imaginary unit, um(t) is the m-th decomposition mode function, omega, obtained after decompositionmIs the center frequency of the mth mode function, delta (t) is the dirac function,
Figure BDA0002263123090000048
is umThe Hilbert transform of (t) functions to convert it into an analytic signal having a single-sided spectrum.
The method carries out two-dimensional windowing and blocking processing on noisy data, applies a Variational Modal Decomposition (VMD) algorithm to F-X spectrum decomposition of seismic exploration data in desert regions, and decomposes data of each row in an F-X spectrum of a local data block along a displacement direction. Because the VMD algorithm can not decompose complex data, a real part and an imaginary part are respectively decomposed and then combined, so that the separation of effective signals and noise components in a spectrum is realized, and the effective recovery of weak signals is realized by effectively decomposing the form and the frequency of seismic data in an F-X spectrum. Because the decomposition process is carried out in the same frequency, the problem of partial frequency component loss caused by time domain decomposition and sub-mode function reconstruction is avoided, and the dual factors of the similarity between the frequency and adjacent channels, the form and the frequency of signals and noise are considered, so that the integrity of the recovered effective signals is better.
The method has the beneficial effects that according to the complex and low-frequency characteristics of noise in the desert seismic exploration data, a VMD decomposition algorithm in time-frequency analysis is introduced to provide a new desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition. The method considers the similarity of the same phase axes between adjacent seismic channels and the diversity and time-varying property of effective signals in seismic exploration data, carries out two-dimensional windowing block processing on noisy data, applies the VMD algorithm to the decomposition of a noisy data F-X spectrum along the displacement direction, and adopts a mode of respectively decomposing a real part and an imaginary part, thereby solving the problem that the VMD algorithm can not decompose complex data. The decomposition process is carried out in the same frequency, so that the problem of partial frequency component loss caused by sub-mode function reconstruction is avoided, dual factors of signal and noise form and frequency are considered, the recovered signal has better integrity, the calculation process is simple, the speed is higher, the accuracy is high, and the method has important significance for desert seismic data processing and subsequent parameter inversion problems.
Drawings
FIG. 1 is a schematic diagram of a two-dimensional windowing implementation of block processing of noisy data according to the present invention, where ntIs the time window length and nlThe two-dimensional window slides rightwards along the horizontal direction for the length of the displacement window, and a certain overlap is formed between adjacent windows; then moving downwards to a new position, continuing sliding rightwards along the horizontal direction, and so on to realize the block processing of the seismic data;
FIG. 2(a) is a diagram of a multi-channel waveform of a data block obtained by selecting a suitable two-dimensional window from a noisy seismic data block, wherein the lower right side of the multi-channel waveform comprises an effective event axis;
fig. 2(b) is an F-X amplitude spectrum obtained by performing fourier transform on each column of data in fig. 2(a), and it can be seen that effective signals and noise components are overlapped together, and the noise frequency is low;
FIG. 3(a) is a F-X magnitude spectrum of a clean signal within a certain two-dimensional window;
FIG. 3(b) is a F-X magnitude spectrum of noisy data over a certain two-dimensional window, where it can be seen that the frequency of the noise component is low;
FIG. 3(c) is a F-X amplitude spectrum of noisy data contained within a certain two-dimensional window, with significant signal components superimposed on the noise components;
FIG. 3(d) is the real part data of the spectral data after decomposition and effective sub-modal reconstruction along the displacement direction using VMD algorithm;
FIG. 3(e) is the imaginary data of the spectral data after decomposition and effective sub-modal reconstruction along the displacement direction using VMD algorithm;
FIG. 3(F) is a plot of the processed F-X amplitude spectrum, where the effective signal components remain intact and the noise is effectively suppressed;
FIG. 4(a) clean synthetic seismic data. The device comprises 5 same-phase axes, 2 more complete same-phase axes and 3 shorter same-phase axes, wherein the inclination directions of the same-phase axes are different;
fig. 4(b) desert noise data. The noise is intercepted from the actual seismic data of the desert area, so that the noise has larger fluctuation;
FIG. 4(c) synthesized noisy seismic data. The record is obtained by adding pure seismic data and desert noise data together, wherein a certain distortion occurs after a homophase axis is interfered by noise;
fig. 4(d) processing results of the EMD method. It can be seen that part of the noise is effectively eliminated, but more noise remains, especially near the in-phase axis, and the in-phase axis energy is weakened to a certain extent, and the amplitude becomes smaller;
FIG. 4(e) processing results of the VMD method. Most of the noise is effectively eliminated, but the waveform of the effective in-phase axis is distorted, and the energy is greatly lost;
fig. 4(f) the processing results of the method in this patent show that most of the noise in the processed seismic records is effectively eliminated, the effective event can be clearly recovered, and the recovery effect of the weak effective event is more continuous and complete compared with the processing results of the first two methods;
FIG. 5 is a comparison graph of processing results of noisy seismic data and three methods, and it can be seen from waveforms that there are large fluctuations in the noisy seismic data, and effective signals are completely submerged in the noisy seismic data and cannot be identified; after the three methods are used for processing, the noise is eliminated to a certain extent, and the effective signal position and the waveform can be displayed; comparing the results of the three methods for recovering the waveform, the EMD method can recover the effective signal position and partial waveform, and simultaneously still has larger waveform noise; the waveform form of the effective signal recovered by the VMD method is good, the noise suppression is also good, but the energy of the recovered effective signal is seriously lost, and the amplitude is very small; compared with the recovery result of the VMD method, the effective signal recovered by the method has better amplitude preservation and better noise suppression effect;
fig. 6(a) shows that in the noisy seismic data in the actual desert area, the noise before first arrival has large fluctuation, the effective event is submerged in the noise in the data after first arrival, waveform distortion occurs, and part of the original event forms cannot be identified, as shown by the oval circled area in the figure;
fig. 6(b) shows the recovery result of the EMD method, and it can be seen from the figure that the noise before the first arrival is effectively suppressed, the noise in the data after the first arrival is also reduced to some extent, and the in-phase axis becomes clearer, but the effective in-phase axis is removed in the area enclosed by the rectangular frame in the figure, and the shaft fracture occurs; the coherence of the same phase axis in the elliptical area is not good;
as can be seen from the recovery result of the VMD method in fig. 6(c), the noise suppression effect before first arrival is good, and the homophase axis recovery in the data after first arrival clearly appears, but the coherence of the homophase axis is still poor, as shown by the oval area in the figure;
fig. 6(d) the present invention recovers the results, the post-processing pre-first arrival noise and the post-first arrival noise are both better eliminated, the in-phase axis becomes clearly visible while the continuity is better, and the energy of the effective signal is not too much lost.
Detailed Description
Comprises the following steps:
1) based on the diversity of the same-phase axes in seismic exploration data and the time-varying property of signals, firstly, two-dimensional windowing localization processing is adopted for noisy seismic data, and the similarity of data characteristics in each window is ensured as much as possible; for a sheet n1×n2Contains n2Each seismic trace of n1Selecting a proper time window length n with the seismic channel direction as the displacement directiontAnd displacement window length nlSliding the two-dimensional window on the noisy seismic data along the horizontal direction to the right, then moving the two-dimensional window downwards, and then continuously sliding the two-dimensional window along the horizontal direction to the right, wherein a certain coincidence exists between two adjacent windows, and the coincidence coefficient ξ (0) is used<ξ<1) To illustrate, fig. 1 shows a schematic diagram of a windowing blocking process, and finally, a blocking process is performed on the whole noisy seismic data D. Each data block is nt×nlUsing the divided ith data block as D(i)(i is 1,2, …, N), where N is the number of divided data blocks:
D(i)＝[D(i)(a,b)]，a＝1～nt，b＝1～nl is connected with(1)
Where a and b are sample point numbers, then D(i)Column b in (1) for data D(i)(:,b)The a-th data is represented by D(i)(a,:)To represent;
then, for the noisy data block D(i)Each column of data (i.e., in the time direction) in (i-1, 2, …, N) is fourier transformed by D(i)Column b data D in (1)(i)(:,b)For example, equation (2) gives the equation for which the fourier transform is performed:
Figure BDA0002263123090000071
where k is the frequency domain sampling point number, n is the time domain sampling point number, j is the imaginary unit, F(i),(:,b)(k) As time domain data D(i),(:,b)(n) corresponding to the frequency domain data, D can be calculated by the equation (2)(i)In the frequency domain data corresponding to each column of data, then D(i)F-X spectrum F of(i)It can be obtained by combining the frequency domain data of each column of data, i.e. F(i)＝[F(i)(:,1)(k),F(i),(:,2)(k),…,F(i),(:,nl)(k)]Analogizing the situation of other data blocks, wherein fig. 2(a) is a multi-channel waveform of a certain data block intercepted from certain synthesized noisy seismic data, and it can be seen that the right lower side of the data block comprises an effective homophase axis, and each line of data is respectively subjected to fourier transform and then combined to obtain an F-X spectrum of the data block, and fig. 2(b) provides an F-X amplitude spectrum of the data block, wherein the distribution situation of each single frequency in different channels is relatively similar and has strong distribution stability with the distribution situation of different frequencies in the same channel, so that the frequency component loss situation is obviously reduced when the data block is processed along a seismic channel (displacement direction);
2) VMD is an adaptive and non-recursive decomposition algorithm that can be viewed as a constrained variational optimization problem, whose continuous form is shown in equation (3):
Figure BDA0002263123090000072
wherein f (t) is data to be decomposed, t is time variable, j is imaginary unit, um(t) is the mth sub-modal function, ω, obtained after decompositionmIs the center frequency of the mth sub-modal function, delta (t) is the dirac function,
Figure BDA0002263123090000073
is um(t) a hilbert transform, operative to convert the hilbert transform into an analytic signal having a single-sided frequency spectrum; the constraint variation problem in equation (3) can be solved by the lagrangian equation in the following equation:
Figure BDA0002263123090000081
wherein α is Lagrangian, um，ωmAnd λ may be updated gradually according to equation (3):
where c is the number of iterations and,
Figure BDA0002263123090000083
and
Figure BDA0002263123090000084
respectively, f (t), ui(t), λ (t) and
Figure BDA0002263123090000085
a Fourier transform of (1);
according to the formula (3), a complex function (or signal) can be decomposed into a plurality of sub-mode functions with different frequencies through a VMD algorithm, the sum of the sub-mode functions is an original decomposition function (or signal), and the decomposition meaning is that information in different frequency bands in complex data can be separated, so that effective information can be conveniently extracted;
then, applying the VMD algorithm to the F-X spectrum data of each data block obtained in the step 1) to realize the separation of effective signals and noise components in the F-X spectrum data, and obtaining the data block D from the step 1)(i)Containing nlColumn, ntThe line number of data points, the corresponding F-X spectrum F(i)Also contains nlColumn, ntLine data points, F(i)May be denoted as F(i)(a,:)The VMD algorithm is used to decompose the row vector (i.e., along the direction of displacement) of the F-X spectral data, due to F(i)Is complex, the direct application of VMD algorithm can not realize the decomposition of complex number, so we adopt the way of respectively decomposing and recombining real part and imaginary part to realize the decomposition process of F-X spectrum data, the formula (4) gives F(i)A row data F(i)(a,:)The complex form of (a):
F(i)(a,:)＝F(i)(a,:)real+jF(i)(a,:)img，a＝1～nt(4)
wherein j is an imaginary unit, F(i)(a,:)realAnd F(i)(a,:)imgF is divided into real part data and imaginary part data by using VMD decomposition algorithm in formula (3)(i)(a,:)realAnd F(i)(a,:)imgDecomposing the data into a plurality of sub-mode functions, each sub-mode function after decomposition and F(i)(a,:)realAnd F(i)(a,:)imgThe relationship between is
Figure BDA0002263123090000086
Where K is the number of sub-mode functions, usually between 4 and 7, and up,realAnd up,imgRespectively representing real part data F(i)(a,:)realAnd imaginary data F(i)(a,:)imgThe p-th sub-modal function of (1).
3) To F(i)(a,:)realAnd F(i)(a,:)imgDecomposed sub-mode functions up,realAnd up,imgThe method comprises two main categories of a sub-mode function taking signal characteristics as a main characteristic and a sub-mode function taking noise as a main characteristic, wherein signal components are mainly collected in the signal-dominant sub-mode function, and the noise is collected by the noise-dominant sub-mode function. In order to accurately discriminate the class of the sub-mode function, the difference factor ε is defined by the equations (7) and (8), respectivelyrealAnd εimgTo measure whether the sub-modal function is a signal-dominant or noise-dominant modal function:
Figure BDA0002263123090000091
Figure BDA0002263123090000092
where s (k) denotes the desired signal, noi (k) denotes noise, up,real(k) And up,img(k) P-th sub-modal functions, ε, representing real and imaginary data, respectivelyrealAnd εimgThe larger the value of the sub-mode function is, the larger the difference between the sub-mode function and the effective signal is, and when the difference value is larger than a certain threshold value, the sub-mode function is judged to be a sub-mode function taking noise as the leading factor; otherwise, the signal is judged to be a sub-mode function taking the signal as the leading factor;
assuming that the first M mode functions in each decomposed sub-mode function are signal-dominant sub-mode functions, then equations (5) and (6) can be rewritten as:
Figure BDA0002263123090000094
wherein u isp,real signalAnd up,img signalRepresenting signal-dominant sub-mode functions, u, of sub-mode functionsp,real noiseAnd up,img noiseRepresenting a noise-dominant sub-mode function in the sub-mode functions, wherein M represents the number of the signal-dominant sub-mode functions, discarding the noise-dominant sub-mode functions, and summing the M signal-dominant sub-mode functions to obtain:
Figure BDA0002263123090000096
then will be
Figure BDA0002263123090000097
And
Figure BDA0002263123090000098
combining into complex form to obtain the processed F-X spectrum as follows:
Figure BDA0002263123090000101
wherein
Figure BDA0002263123090000102
Is F(i)(a,:)Processed F-X spectral data, F(i)The other rows of data are analogized, and F can be obtained after the data of each row are processed in the same way(i)F-X spectra after treatmentFIG. 3 shows a comparison analysis of clean signal, noisy data, and real part data and imaginary part data of VMD decomposition of noisy data in F-X spectrum displacement direction and processed F-X spectrum in a certain two-dimensional window; fig. 3(a), (b) and (c) are F-X amplitude spectra of clean signals, noisy data and noisy data, and fig. 3(d), (e) and (F) are processed F-X amplitude spectra obtained after VMD decomposition along the displacement direction, and summation reconstruction of real part data, imaginary part data and effective sub-mode functions. Comparing the F-X amplitude spectrum of the original noisy data with the F-X amplitude spectrum after decomposition processing shows that the noise components are effectively suppressed, the effective signal components are kept more complete, and only less noise is left;
finally, the processed F-X spectrum data is processed
Figure BDA0002263123090000104
The data in each column in (1) is subjected to inverse Fourier transform, and the data in the b-th column
Figure BDA0002263123090000105
The result of performing the inverse fourier transform is:
Figure BDA0002263123090000106
wherein k is the serial number of the sampling point in the frequency domain, n is the serial number of the sampling point in the time domain, j is the unit of imaginary number,
Figure BDA0002263123090000107
the data blocks are corresponding time domain sequences, namely the recovered effective signals, and the rest of the column data are analogized, so that the recovered multi-channel effective signal data blocks are obtainedFinally, all the effective signal data blocks are combined to obtain a whole effective signal data matrix
Figure BDA0002263123090000109
The effects of the present invention will be described below by two specific experimental examples.
Experimental example 1 Synthesis record
Fig. 4(a) shows a piece of 160 synthetic pure seismic data, each of 2048 sampling points, which contains 1 complete long event and 4 short event in order to consider the distribution of different event, and the tilt directions of the event are different from each other. Fig. 4(b) shows desert noise, which is obtained by intercepting the first arrival noise in the actual desert seismic data, and the fluctuation range is large and the dominant frequency is low. The noise-containing seismic data in fig. 4(c) can be obtained by adding the clean seismic data and the noise, and it can be seen that the effective signal is submerged in the large-amplitude noise, and part of the waveform is distorted. And processing the noisy data by respectively utilizing the EMD, the VMD and the method provided by the invention. Setting the mode number as 5 in the EMD method, and selecting a mode 2 and a mode 3 to recover effective signals; in the VMD method, the number of modes is set to be 4, and the modes 2 and 3 are selected to reconstruct effective signals. In the method of the invention, the length of the time window n is selectedt64, displacement window length nlThe window overlap factor is 0.6 at 50, and the three methods recover valid signals as shown in fig. 4(d) -4 (f).
As can be seen from the recovery results of the three methods, firstly, part of the noise after the EMD method in fig. 4(d) is effectively suppressed, but more noise remains, and the in-phase axis can be displayed, but the strength is weakened, which indicates that the effective components are reduced; then, after the VMD method in fig. 4(e) processes, noise is effectively reduced, the in-phase axis becomes clear, but the wavelet waveform of the in-phase axis is severely distorted, which is caused by that part of the effective components contained in the noise mode is directly discarded; finally, most of the noise is effectively suppressed after the processing of the method in fig. 4(f), the in-phase axis is clearly shown, and the coherence and the smoothness are good.
The 80 th data is extracted from the noisy data, the EMD method processing result and the VMD method processing result respectively, and the single-channel waveforms are drawn for comparison, as shown in FIG. 5. It can be seen from the figure that the noisy data contains large-amplitude noise, and the effective signal is submerged in the noisy data and cannot be identified; after the EMD method is used for processing, partial noise is suppressed, but the noise with increased amplitude is still mistaken for a valid signal and is reserved; after the VMD method is processed, most of noise is effectively suppressed, but effective components are seriously lost while the noise is suppressed, and the amplitude of an effective signal is obviously reduced and more effective components are discarded as noise components as can be seen from the waveform after the VMD processing; the method of the invention clearly shows the effective signals in the processed waveform, and effectively suppresses the noise.
Experimental example 2 practical record
Finally, the actual seismic data from the Tarim oil field is selected to further verify the performance of the method. FIG. 6(a) is a 135 trace of actual noisy seismic data at 3800 sample points per trace, at a sampling frequency of 500 Hz. The processing results of the noise suppression and the effective signal recovery by using the EMD, the VMD and the method of the invention are shown in FIGS. 6(b) -4 (d). It can be seen from observation that the effective signal in the original noisy data is covered by strong random noise, and part of the waveform is distorted. After the processing by the EMD method, the first arrival noise in fig. 6(b) is basically eliminated, most of strong noise and part of the surface wave are effectively suppressed, but the middle effective signal is missing, as shown by the circled part of the rectangle in the figure. After the VMD method in fig. 6(c), strong noise is effectively suppressed in addition to the first arrival noise and the surface wave, but the effective signal with both sides submerged in the noise is not effectively recovered. Finally, after the processing by the method of the present invention, it can be seen from the result in fig. 6(d) that the first arrival noise is effectively eliminated, part of the surface wave is suppressed, the effective signal in the middle position is completely retained, and the weak effective signal submerged in the strong noise at the two side positions is clearly and continuously recovered, as shown by the oval circled part in the figure.

Claims (2)

1. A desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition is characterized by comprising the following steps:
1) based on the diversity of the same-phase axes in seismic exploration data and the time-varying property of signals, firstly, two-dimensional windowing localization processing is adopted for noisy seismic data, and the similarity of data characteristics in each window is ensured; for a sheet n1×n2Selecting proper time window length n from the noisy seismic data DtAnd displacement window length nlSliding the two-dimensional window on the noisy seismic data, with a certain coincidence between two adjacent windows, using a coincidence coefficient ξ (0)<ξ<1) Representing to realize the block processing of the noisy seismic data D; dividing the ith data block into D(i)(i is 1,2, …, N), where N is the number of divided data blocks:
D(i)＝[D(i)(a,b)]，a＝1～nt，b＝1～nl(1)
thus D(i)Column b in (1) for data D(i)(:,b)The a-th data is represented by D(i)(a,:)To represent;
then, for the noisy data block D(i)Each column of data (i: 1,2, …, N), that is, data in the time direction is fourier-transformed by D(i)Column b data D in (1)(i)(:,b)For example, equation (2) gives the equation for which the fourier transform is performed:
Figure FDA0002263123080000011
where k is the frequency domain sampling point number, n is the time domain sampling point number, j is the imaginary unit, F(i),(:,b)(k) As time domain data D(i),(:,b)(n) corresponding frequency domain data; d can be calculated by the formula (2)(i)In the frequency domain data corresponding to each column of data, then D(i)F-X spectrum F of(i)It can be obtained by combining the frequency domain data of each column of data, i.e. F(i)＝[F(i)(:,1)(k),F(i),(:,2)(k),…,F(i),(:,nl)(k)]So as to analogize the situation of other data blocks;
2) applying VMD algorithm to the F-X spectrum data of each data block obtained in the step 1) to realize the separation of effective signals and noise components in the F-X spectrum data, decomposing the row vector of the F-X spectrum data, namely along the displacement direction by utilizing the VMD algorithm, and performing F-X spectrum analysis on the F-X spectrum data(i)Is denoted as F(i)(a,:)Due to F(i)The F-X spectrum data decomposition process is realized by respectively decomposing and recombining a real part and an imaginary part in a mode of directly applying a VMD algorithm, wherein the complex number cannot be decomposed; formula (4) gives F(i)(a,:)The complex form of (a):
F(i)(a,:)＝F(i)(a,:)real+jF(i)(a,:)img，a＝1～nt(4)
wherein j is an imaginary unit, F(i)(a,:)realAnd F(i)(a,:)imgRespectively real part data and imaginary part data;
then, F is converted by the formula (3)(i)(a,:)realAnd F(i)(a,:)imgRespectively decomposing the data into a form of weighted sum of K sub-modal functions, and decomposing each sub-modal function and F(i)(a,:)realAnd F(i)(a,:)imgThe relationship between them is:
Figure FDA0002263123080000012
Figure FDA0002263123080000013
where K is the number of sub-mode functions, usually between 4 and 7, and up,realAnd up,imgRespectively represent F(i)(a,:)realAnd F(i)(a,:)imgThe p-th sub-modal function of (1);
3) function of submode up,realAnd up,imgThe method comprises two categories of a sub-mode function taking signal characteristics as a leading factor and a sub-mode function taking noise as a leading factor, wherein signal components are mainly collected in the signal leading sub-mode function, and the noise is collected by the noise leading sub-mode function;
in order to accurately discriminate the class of the sub-mode function, the difference factor ε is defined by the equations (7) and (8), respectivelyrealAnd εimgTo measure whether the sub-modal function is a signal-dominant or noise-dominant modal function;
where s (k) denotes the desired signal, noi (k) denotes the noise,. epsilonrealAnd εimgThe larger the value of the sub-mode function is, the larger the difference between the sub-mode function and the effective signal is, and when the difference value is larger than a certain threshold value, the sub-mode function is judged to be a sub-mode function taking noise as the leading factor; otherwise, the signal is judged to be a sub-mode function taking the signal as the leading factor;
when the first M mode functions in each decomposed sub-mode function are signal-dominant sub-mode functions, then the equations (5) and (6) can be rewritten into
Figure FDA0002263123080000023
Figure FDA0002263123080000024
Wherein u isp,real signalAnd up,img signalRepresenting signal-dominant sub-mode functions, up,real noiseAnd up,img noiseRepresenting by noiseThe number of the sub-mode functions taking the signal as the leading part is M; discarding the noise-dominant sub-mode functions, and summing M signal-dominant sub-mode functions to obtain:
Figure FDA0002263123080000026
then will be
Figure FDA0002263123080000027
And
Figure FDA0002263123080000028
combining into complex form to obtain the processed F-X spectrum as follows:
Figure FDA0002263123080000029
wherein
Figure FDA0002263123080000031
Is F(i)A row data F(i)(a,:)Processed F-X spectral data, F(i)The other rows of data are analogized, and F can be obtained after the data of each row are processed in the same way(i)F-X spectra after treatment
Figure FDA0002263123080000032
Finally, to
Figure FDA0002263123080000033
The data in each column in (1) is subjected to inverse Fourier transform, and the data in the b-th column
Figure FDA0002263123080000034
The result of inverse fourier transform on b-1-nl is:
Figure FDA0002263123080000035
wherein k is the serial number of the sampling point in the frequency domain, n is the serial number of the sampling point in the time domain, j is the unit of imaginary number,the data is a corresponding time domain sequence after inverse transformation, namely a recovered effective signal, and other column data are analogized, so that a plurality of recovered effective signal data blocks are obtained
Figure FDA0002263123080000037
Combining each effective signal data block to obtain a whole effective signal data matrix
Figure FDA0002263123080000038
2. The method for recovering weak signals in desert exploration based on VMD local F-X spectrum decomposition as claimed in claim 1, wherein: in step 2), VMD is an adaptive and non-recursive decomposition algorithm, and its continuous form is shown in formula (3):
Figure FDA0002263123080000039
wherein f (t) is data to be decomposed, t is time variable, j is imaginary unit, um(t) is the mth sub-modal function, ω, obtained after decompositionmIs the center frequency of the mth sub-modal function, delta (t) is the dirac function,
Figure FDA00022631230800000310
is um(t) Hilbert transform, which is operative to convert the Hilbert transform into an analytic signal having a single-sided spectrum.
CN201911084413.1A 2019-11-06 2019-11-06 Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition Active CN110764147B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911084413.1A CN110764147B (en) 2019-11-06 2019-11-06 Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911084413.1A CN110764147B (en) 2019-11-06 2019-11-06 Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition

Publications (2)

Publication Number Publication Date
CN110764147A true CN110764147A (en) 2020-02-07
CN110764147B CN110764147B (en) 2020-11-03

Family

ID=69336776

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911084413.1A Active CN110764147B (en) 2019-11-06 2019-11-06 Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition

Country Status (1)

Country Link
CN (1) CN110764147B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112766044A (en) * 2020-12-28 2021-05-07 中海石油(中国)有限公司 Method and device for analyzing longitudinal and transverse wave speeds of loose sample and computer storage medium

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107589454A (en) * 2017-07-25 2018-01-16 西安交通大学 One kind is based on VMD TFPF compacting seismic prospecting random noise methods
US20180188399A1 (en) * 2017-01-05 2018-07-05 China Petroleum & Chemical Corporation Anisotropy matching filtering for attenuation of seismic migration artifacts
CN108919347A (en) * 2018-07-02 2018-11-30 东华理工大学 Seismic signal stochastic noise suppression method based on vmd
CN109061724A (en) * 2018-06-21 2018-12-21 北京化工大学 A kind of seismic data noise-reduction method based on adaptive variation mode decomposition
CN109977914A (en) * 2019-04-08 2019-07-05 哈尔滨工业大学 Self-adaptation noise reduction method based on VMD
CN110208856A (en) * 2019-06-05 2019-09-06 吉林大学 A kind of desert Complex Noise drawing method based on manifold subregion 2D-VMD

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180188399A1 (en) * 2017-01-05 2018-07-05 China Petroleum & Chemical Corporation Anisotropy matching filtering for attenuation of seismic migration artifacts
CN107589454A (en) * 2017-07-25 2018-01-16 西安交通大学 One kind is based on VMD TFPF compacting seismic prospecting random noise methods
CN109061724A (en) * 2018-06-21 2018-12-21 北京化工大学 A kind of seismic data noise-reduction method based on adaptive variation mode decomposition
CN108919347A (en) * 2018-07-02 2018-11-30 东华理工大学 Seismic signal stochastic noise suppression method based on vmd
CN109977914A (en) * 2019-04-08 2019-07-05 哈尔滨工业大学 Self-adaptation noise reduction method based on VMD
CN110208856A (en) * 2019-06-05 2019-09-06 吉林大学 A kind of desert Complex Noise drawing method based on manifold subregion 2D-VMD

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
方江雄 等: "基于变分模态分解的地震随机噪声压制方法", 《石油地球物理勘探》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112766044A (en) * 2020-12-28 2021-05-07 中海石油(中国)有限公司 Method and device for analyzing longitudinal and transverse wave speeds of loose sample and computer storage medium
CN112766044B (en) * 2020-12-28 2024-03-22 中海石油(中国)有限公司 Method and device for analyzing longitudinal and transverse wave speeds of loose sample and computer storage medium

Also Published As

Publication number Publication date
CN110764147B (en) 2020-11-03

Similar Documents

Publication Publication Date Title
CN103543469B (en) A kind of little yardstick Threshold Denoising Method based on wavelet transformation
Schimmel et al. The inverse S-transform in filters with time-frequency localization
CN107817527B (en) Seismic exploration in desert stochastic noise suppression method based on the sparse compressed sensing of block
CN102221708B (en) Fractional-Fourier-transform-based random noise suppression method
CN107102356B (en) Seismic signal high resolution data processing methods based on CEEMD
Tian et al. Variable-eccentricity hyperbolic-trace TFPF for seismic random noise attenuation
CN101598812B (en) Method for removing abnormal noise in single-point reception of seismic record by digital detector
CN109031422A (en) A kind of seismic signal noise suppressing method based on CEEMDAN and Savitzky-Golay filtering
WO2008112036A1 (en) Imaging of multishot seismic data
Leblanc et al. Denoising of aeromagnetic data via the wavelet transform
CN110208856B (en) Desert complex noise suppression method based on manifold partition 2D-VMD
CN104849757A (en) System and method for eliminating random noise in seismic signals
CN107436450A (en) A kind of seismic signal bandwidth broadning method based on continuous wavelet transform
CN113077386A (en) Seismic data high-resolution processing method based on dictionary learning and sparse representation
CN108961181A (en) A kind of ground penetrating radar image denoising method based on shearlet transformation
CN110764147B (en) Desert exploration weak signal recovery method based on VMD local F-X spectrum decomposition
Xi et al. Spurious signals attenuation using SVD-based Wiener filter for near-surface ambient noise surface wave imaging
CN107783191A (en) The method that hyperspace space-time time-frequency method cuts down seismic prospecting random noise
CN111308555B (en) Surface wave interference elimination method based on inter-frequency similarity low-frequency reconstruction
Zhang et al. Seismic random noise attenuation and signal-preserving by multiple directional time-frequency peak filtering
CN105487115A (en) Wavelet transform-based high frequency continuation method
CN110109179A (en) Bandwidth compensation processing method, device and equipment
CN112213785B (en) Seismic data desert noise suppression method based on feature-enhanced denoising network
CN113093282A (en) Desert data denoising method based on geometric modal characteristic parallel network
CN115236744A (en) Method for strong shielding stripping and weak signal energy compensation of igneous rock

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
GR01 Patent grant
GR01 Patent grant