CN109782346B - Acquisition footprint pressing method based on morphological component analysis - Google Patents

Acquisition footprint pressing method based on morphological component analysis Download PDF

Info

Publication number
CN109782346B
CN109782346B CN201910033429.3A CN201910033429A CN109782346B CN 109782346 B CN109782346 B CN 109782346B CN 201910033429 A CN201910033429 A CN 201910033429A CN 109782346 B CN109782346 B CN 109782346B
Authority
CN
China
Prior art keywords
dimensional
signal
discrete cosine
transform
phi
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.)
Active
Application number
CN201910033429.3A
Other languages
Chinese (zh)
Other versions
CN109782346A (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.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong 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 Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201910033429.3A priority Critical patent/CN109782346B/en
Publication of CN109782346A publication Critical patent/CN109782346A/en
Application granted granted Critical
Publication of CN109782346B publication Critical patent/CN109782346B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention discloses a method for pressing a collected footprint based on morphological component analysis, which comprises the following steps: step 101: constructing two-dimensional local discrete cosine transform according to the acquired footprint waveform morphological characteristics in the seismic record of the stratosphere slice or the isochronous slice of the three-dimensional seismic data volume, and combining the two-dimensional local discrete cosine transform with the two-dimensional stationary wavelet transform to form a super-complete dictionary; step 102: performing preliminary pressing of the collected footprints in each layer slice or isochronal slice by using a morphological component analysis-based method in a layer-by-layer execution mode on the original data of the seismic recording signals; step 103: determining a low cut-off frequency of the two-dimensional local discrete cosine transform; step 104: and repeating the step 102 and the step 103 until all the slice data are processed, performing signal-noise separation by using a method based on morphological component analysis, and suppressing the acquired footprint noise, thereby finally realizing the suppression of the acquired footprint noise in the three-dimensional seismic data volume. The method is used for carrying out acquisition footprint pressing on the three-dimensional seismic data volume, and the aim of improving the signal-to-noise ratio of seismic data is fulfilled.

Description

Acquisition footprint pressing method based on morphological component analysis
Technical Field
The invention belongs to the technical field of data processing in seismic exploration, and particularly relates to a method for suppressing a collected footprint in seismic collected data.
Background
Acquisition footprints are typically a periodic amplitude artifact generated in seismic imaging due to the rolling arrangement of the seismic survey system and incomplete sampling caused by the source and receiver side-line spacing, and are typically observed in time or depth slices of three-dimensional seismic data. In the exploration process of stratum-lithology and other concealed oil and gas reservoirs, in order to meet the demand of complex reservoir prediction, fine structural interpretation and lithology trap interpretation need to be carried out on seismic data. However, the collected footprint noise can seriously interfere with the seismic interpretation process, and can even be misjudged as geological structure or lithology abnormity, which seriously affects the accuracy and reliability of seismic data interpretation. Therefore, the method for collecting the footprints by pressing and attenuating is very important to research, and has great theoretical significance and market value.
With the development of the signal sparsity theory, Starck et al propose a mixed signal decomposition method of morphological component analysis. The morphological component analysis means that two transformation dictionaries with different atomic features form an ultra-complete dictionary according to the morphological features of the composition components of the complex signal, so that a more sparse representation mode and more effective information recognition capability of the complex signal are realized, and the separation of the two components is realized. The dictionary is usually selected or constructed from known mathematical transformations empirically, and the assumption whether the selected dictionary sufficiently satisfies morphological analysis is the key to success of the morphological analysis method. The purpose of suppressing harmonic noise is achieved by separating signals and harmonic noise by using a morphological component analysis method, and proper sparse transformation needs to be constructed according to waveform morphological characteristics of effective signals and harmonic noise. It is appropriate to have a literature study to choose to sparsely represent the effective signal using continuous wavelet transform. However, harmonic noise is complex, and a transform which is suitable for matching needs to be constructed to be sparsely represented, and specific parameters of the transform are selected and determined, so that related researches are few at present.
The prior art is as follows:
and optimizing the acquisition method. The wide azimuth observation system is adopted, so that the nonuniformity of azimuth distribution of data recording can be reduced, and omnibearing observation is carried out on the ground, so that wave field characteristics of underground reflection points in all directions are recorded by seismic data, and complete seismic wave field information is obtained; the wide azimuth observation system also helps to reduce shot backscatter noise and improve static correction effects, which all help to reduce acquisition footprints.
The prior art has the following disadvantages:
the data collected before can not be processed only aiming at the newly collected data; different designs need to be made for different acquisition systems, and the transportability is poor.
Disclosure of Invention
The invention aims to provide a method for pressing a collected footprint based on morphological component analysis, so as to solve the technical problem. Aiming at the difference of the waveform morphological characteristics of the acquisition footprint and the effective wave signal in the layered slice or the isochronous slice of the three-dimensional seismic data body, proper sparse representation transformation is constructed, namely two-dimensional local discrete cosine transformation is selected to represent the acquisition footprint, two-dimensional stationary wavelet transformation is selected to represent the effective wave signal, the two waveform morphological dictionaries are combined to form a super-complete dictionary, two-dimensional stationary wavelet transformation parameters are selected and determined according to the data characteristics, and a coordinate blocking relaxation algorithm method is used for realizing signal-noise separation so as to achieve the purpose of effectively suppressing the acquisition footprint in the seismic data of the layered slice or the isochronous slice of the three-dimensional seismic data body.
In order to achieve the purpose, the invention adopts the following technical scheme:
the acquisition footprint pressing method based on morphological component analysis comprises the following steps:
step 101: constructing two-dimensional local discrete cosine transform according to the acquired footprint waveform morphological characteristics in the seismic record of the edge slice or the isochronous slice of the three-dimensional seismic data volume, combining the two-dimensional local discrete cosine transform with the two-dimensional stationary wavelet transform to form a super-complete dictionary, and simultaneously determining the decomposition series J of the two-dimensional stationary wavelet transform and the analysis window size W of the two-dimensional local discrete cosine transform;
step 102: performing preliminary pressing of the collected footprints in each layer slice or isochronal slice by using a morphological component analysis-based method in a layer-by-layer execution mode on the original data of the seismic recording signals;
step 103: determining a low cut-off frequency of the two-dimensional local discrete cosine transform;
step 104: and repeating the steps 02-03 until all the slice data are processed, performing signal-noise separation by using a method based on morphological component analysis, suppressing the acquired footprint noise, and finally realizing the suppression of the acquired footprint noise in the three-dimensional seismic data volume.
Furthermore, in step 101, two transformation dictionaries for morphological component analysis are determined according to the waveform morphological characteristics of the effective signals and the waveform morphological characteristics of the collected footprints, two-dimensional stationary wavelet transformation is selected for the effective signals, local discrete cosine transformation is selected for the collected footprints, and an ultra-complete dictionary is formed.
The object of applying morphological component analysis is to contain two components with different morphological characteristics:
s＝s1+s2，
in the formula: s represents a signal to be analyzed; s1、s2Representing two components of the signal having different morphological characteristics; respectively extract s1、s2These two components are the targets of morphological component analysis; suppose s1And s2Can be respectively composed of a dictionary phi1And phi2Efficient sparse representation, but with phi2Sparse representation s1And by phi1Sparse representation s2Poor temporal sparsity;
selecting a two-dimensional stationary wavelet transform as a transformation dictionary for sparsely representing effective signal components, wherein the two-dimensional stationary wavelet forward transform is:
Figure BDA0001945046010000031
Figure BDA0001945046010000032
Figure BDA0001945046010000033
Figure BDA0001945046010000034
in the formula HjAnd GjRespectively, representing the filter banks of the j-th layer decomposition.
The inverse transformation of the two-dimensional stationary wavelet transform is:
Figure BDA0001945046010000041
in the formula
Figure BDA0001945046010000042
And
Figure BDA0001945046010000043
respectively generation by generationWatch HjAnd GjThe dual filter bank of (1).
Selecting a two-dimensional local discrete cosine transform (type IV) as a transform dictionary for sparsely representing collected footprints, wherein a forward transform of the two-dimensional local discrete cosine transform:
Figure BDA0001945046010000044
wherein f (i, j) represents the signal to be analyzed,
Figure BDA0001945046010000045
two-dimensional local discrete cosine transform coefficients, k, representing a signal to be analyzed1,k2＝0,1,...N-1。
The inverse of the two-dimensional local discrete cosine transform is:
Figure BDA0001945046010000046
according to morphological component analysis theory, using the above selected dictionary phi1I.e. two-dimensional stationary wavelet transform and phi2Namely, two-dimensional local discrete cosine transform, forms an ultra-complete dictionary, sparsely represents a signal s, and calculates a sparse representation coefficient:
Figure BDA0001945046010000047
in the formula: x is the number of1To reconstruct the sum phi of the coefficients1A corresponding portion; x is the number of2To reconstruct the sum phi of the coefficients2A corresponding portion.
Figure BDA0001945046010000048
Is a lagrange multiplier.
Further, in step 102, the block coordinate relaxation algorithm is used for processing the raw seismic record data along-layer slicing or isochronous slicing to realize the suppression of the collected footprints, and the optimization problem can be solved through the coordinate block relaxation algorithm. The basic idea of the coordinate block relaxation algorithm is to compute x in alternating iterations1And x2. The method mainly comprises the following steps:
initialization: initial iteration step number k is 0, initial solution
Figure BDA0001945046010000051
Figure BDA0001945046010000052
For representing the signal component s1The initial solution of the coefficients of (a) to (b),
Figure BDA0001945046010000053
for representing the signal component s2Initial solution of coefficients of (a);
iteration: and increasing the iteration step number k by 1 in each step of iteration, and calculating:
Figure BDA0001945046010000054
Figure BDA0001945046010000055
in the formula, TλIs a hard threshold function;
Figure BDA0001945046010000056
and phi1A pair of positive and negative conversion is formed,
Figure BDA0001945046010000057
and phi2A pair of positive and negative conversion is formed;
termination conditions were as follows: when in use
Figure BDA0001945046010000058
When the value is smaller than the preset value, the influence of continuous iteration on the result is small enough, and the iteration is terminated;
and (3) outputting:
Figure BDA0001945046010000059
Figure BDA00019450460100000510
for the separated signal component s1The final transform coefficients of (a) are,
Figure BDA00019450460100000511
for the separated signal component s2The final transform coefficients of (a);
in the block coordinate relaxation algorithm, the hard threshold function formula is as follows:
Figure BDA00019450460100000512
in the formula:
Figure BDA00019450460100000513
as a function of the hard threshold, λ is the hard threshold,
Figure BDA00019450460100000514
is a coefficient matrix
Figure BDA00019450460100000515
N, N is the size of the coefficient matrix.
Further, step 103, determining a low-frequency-cutoff implementation of the two-dimensional local discrete cosine transform specifically includes:
step 301: taking M isochronous slices of the three-dimensional seismic data volume according to the interval time delta T and using the slices to test the low-frequency-cutoff parameters of the reconstruction transformation of the local cosine transformation, wherein the positions of the slices are T0,tΔT,t2ΔT,…,t(M-1)ΔT；
Step 302: for position tkΔTThe k is 0,1, M-1, a series of low-frequency-cut parameters of reconstruction transformation of two-dimensional local discrete cosine transformation are given, and the separation of obtaining effective wave signals and collecting footprint noise is solved for each low-frequency-cut parameter; determining a low-cut-frequency parameter f for obtaining the best signal-noise separation effect at the current time position by comparing the signal-noise separation effects under the parameterskΔT；
Step 303: benefit toBy linear interpolation, according to the time position tkΔTK-0, 1.., M-1 and its corresponding low-cutoff parameter fkΔTAnd k is 0, 1., M-1, and obtaining low-frequency-cut parameters of reconstruction transformation of two-dimensional local discrete cosine transformation of other time positions of the whole three-dimensional seismic data volume;
further, each time slice of the three-dimensional seismic data volume is processed layer by layer from the shallow layer to the deep layer in step 04 through a solving formula
Figure BDA0001945046010000061
Obtaining effective wave signal and representing vector x for collecting footprint noiseeAnd xfTo obtain the collected footprint noise s removed from the current time slicef＝ΦfxfAnd the corresponding effective wave signal is se＝s-sf。
Further, the preset value is 10-6。
Further, hard threshold λ is taken
Figure BDA0001945046010000062
In a large to small arrangement
Figure BDA0001945046010000063
A value.
Compared with the prior art, the invention has the following beneficial effects:
the invention utilizes the morphological component analysis theory, regards the collected footprint noise as a signal component in the post-stack seismic wave field, constructs proper sparse representation transformation, namely two-dimensional local discrete cosine transformation, forms a super-complete dictionary with two-dimensional stationary wavelet transformation, selects and determines two-dimensional stationary wavelet transformation parameters according to data characteristics, and realizes signal-noise separation by using a coordinate block relaxation algorithm method to achieve the purpose of suppressing harmonic noise. The invention not only can effectively suppress and collect the footprints and remove part of random noise and offset processing false images, but also has higher fidelity of effective signals.
Drawings
FIG. 1 is a diagram of a partial cosine based time domain waveform;
fig. 2 is a bell-shaped window function (k ═ 0);
FIG. 3 is a two-dimensional 8 × 8 discrete cosine transform atom;
FIG. 4 is a two-dimensional stationary wavelet transform atom;
FIG. 5 is a flow chart of the method of the present invention for processing data;
FIG. 6 is a flow chart for determining low-cutoff parameters for a partial cosine transform reconstruction transform;
FIG. 7A is a cross-section of a main line of actual data; FIG. 7B is the active wave signal isolated from FIG. 7A by the method of the present invention; FIG. 7C illustrates the collected footprint noise removed from FIG. 7A by the method of the present invention;
FIG. 8A is a time slice of 1.7s in an actual data volume; FIG. 8B is a schematic representation of the isolated active wave signal of FIG. 8A; fig. 8C is a graph of the collected footprint noise removed from fig. 8A by the method of the present invention.
Detailed Description
In order to make the purpose and technical solution of the present invention more apparent, the present invention will be described in further detail with reference to the accompanying drawings and the detailed description. The exemplary embodiments and descriptions of the present invention are provided to explain the present invention, but not to limit the present invention.
In an embodiment of the present invention, a method for collecting footprint pressing based on morphological component analysis is provided, which includes the following steps:
step 101: constructing two-dimensional local discrete cosine transform according to the acquired footprint waveform morphological characteristics in the seismic record of the edge slice or the isochronous slice of the three-dimensional seismic data volume, combining the two-dimensional local discrete cosine transform with the two-dimensional stationary wavelet transform to form a super-complete dictionary, and simultaneously determining the decomposition series J of the two-dimensional stationary wavelet transform and the analysis window size W of the two-dimensional local discrete cosine transform;
step 102: performing preliminary pressing of the collected footprints in each layer slice or isochronal slice by using a morphological component analysis-based method in a layer-by-layer execution mode on the original data of the three-dimensional seismic data volume;
step 103: determining a low cut-off frequency of the two-dimensional local discrete cosine transform;
step 104: and repeating the step 102 and the step 103 until all the slice data are processed, performing signal-noise separation by using a method based on morphological component analysis, and suppressing the acquired footprint noise, thereby finally realizing the suppression of the acquired footprint noise in the three-dimensional seismic data volume.
In step 01, two-dimensional local discrete cosine transform is constructed according to the morphological characteristics of acquired footprint waveforms in the three-dimensional seismic record, and the two-dimensional local discrete cosine transform and the two-dimensional stationary wavelet transform form a super-complete dictionary, and the method specifically comprises the following steps:
two-dimensional local discrete cosine transform is selected as a transformation dictionary for sparsely representing collected footprints, the commonly used local cosine bases in the discrete cosine transform include an I type and an IV type, the invention mainly adopts an IV type cosine base, and the definition of a local IV type cosine base function (L CB-IV) is as follows:
Figure BDA0001945046010000081
in the above formula, bmn(x) Representing a local cosine basis having a wavenumber index m and a position index n, as shown in FIG. 1, where
Figure BDA0001945046010000082
In order to be the starting position of the interval,
Figure BDA0001945046010000083
in order to be the end position of the interval,
Figure BDA0001945046010000084
is the interval length. Bell window function Bn(x) Is defined in a closed interval
Figure BDA0001945046010000085
Is a smooth function of, and
Figure BDA0001945046010000086
and' are the overlapping radii of the left and right sides of the interval respectively,
Figure BDA0001945046010000087
wherein β (x) is a contour function (incremental or decremental transformation) whose expression is defined recursively as
Figure BDA0001945046010000088
Where k is ≧ 0, and
Figure BDA0001945046010000089
contour function βk+1(x) Will increase with increasing k, e.g. a profile function of 0
Figure BDA00019450460100000810
The corresponding bell-shaped window function at this time is shown in fig. 2.
For function f ∈ CNThe discrete type IV cosine transform (DCT-IV) is defined as follows:
Figure BDA0001945046010000091
and its inverse transform (IDCT-IV) is defined as:
Figure BDA0001945046010000092
selecting a two-dimensional stationary wavelet transform as a transformation dictionary for sparsely representing effective signals, wherein for a given one-dimensional scale function phi (t) and wavelet function phi (t), the scale function and wavelet function of the two-dimensional stationary wavelet transform are obtained by phi (t) and phi (t) in a tensor product mode:
two-dimensional scale function φ (x, y): phi (x, y) phi (x) phi (y);
two-dimensional horizontal wavelet function ΨH(x,y)：ψH(x,y)＝φ(x)ψ(y)；
Two-dimensionalVertical direction wavelet function psiV(x,y)：ψV(x,y)＝ψ(x)φ(y)；
Two-dimensional diagonal direction wavelet function psiD(x,y)：ψD(x,y)＝ψ(x)ψ(y)。
The two-dimensional stationary wavelet transform decomposes the low-frequency part of the signal of the j layer into a low-frequency part of the j +1 layer and high-frequency parts in the vertical, horizontal and diagonal directions, wherein the low-frequency part of the signal is correspondingly arranged as a low-frequency signal and is listed as a low-frequency signal; the horizontal high-frequency part of the signal corresponds to a signal with low frequency and high frequency; the vertical high-frequency part of the signal corresponds to a signal with high frequency and low frequency; the diagonal high frequency portion of the signal corresponds to a row high frequency, column high frequency signal. Using a multi-aperture algorithm to implement a stationary wavelet transform, filter banks H and G are defined, then HjAnd GjFilter banks each representing a j-th layer decomposition by inserting 2 between the respective coefficients of H and Gj-1 zero, for any j ≧ 0, the positive transform of the stationary wavelet transform of level j:
Figure BDA0001945046010000101
Figure BDA0001945046010000102
Figure BDA0001945046010000103
Figure BDA0001945046010000104
if H is presentjAnd GjAre respectively a dual filter bank
Figure BDA0001945046010000105
And
Figure BDA0001945046010000106
then the inverse transform of the stationary wavelet transform of the j-th layer can be obtained as：
Figure BDA0001945046010000107
As shown in fig. 3 and 4, the two-dimensional discrete cosine transform atom and the two-dimensional stationary wavelet transform atom are two-dimensional 8 × 8, respectively, wherein the two-dimensional discrete cosine transform atom in fig. 3 has gradually increasing frequency in the horizontal direction from the left to the right and gradually increasing frequency in the vertical direction from the top to the bottom, and fig. 4 is a two-dimensional stationary wavelet transform atom with three different decomposition scales, each of which has three directions of horizontal, vertical and diagonal, completely meeting the requirements of morphological component analysis theory, and achieving the expected sparse representation and the expected separation effect for different signal components.
According to morphological component analysis theory, using the above selected dictionary phi1I.e. two-dimensional stationary wavelet transform and phi2Namely, two-dimensional local discrete cosine transform constitutes a sparse representation signal s of the overcomplete dictionary, and a sparse representation coefficient is calculated:
Figure BDA0001945046010000108
in the formula: x is the number of1To reconstruct the sum phi of the coefficients1A corresponding portion; x is the number of2To reconstruct the sum phi of the coefficients2A corresponding portion.
Figure BDA0001945046010000109
Is a lagrange multiplier. The optimization problem is solved by a coordinate block relaxation algorithm.
In step 02, the block coordinate relaxation algorithm is used for processing the original seismic record data by slicing along the layer or slicing at equal time, so that the compression of the acquired footprint is realized, and the optimization problem is solved by the coordinate block relaxation algorithm. The basic idea of the coordinate block relaxation algorithm is to compute x in alternating iterations1And x2. The method mainly comprises the following steps:
initialization: initial iteration step number k is 0, initial solution
Figure BDA0001945046010000111
Figure BDA0001945046010000112
The initial solution of coefficients used to represent signal component 1,
Figure BDA0001945046010000113
the initial solution of coefficients used to represent signal component 2;
iteration: and increasing the iteration step number k by 1 in each step of iteration, and calculating:
Figure BDA0001945046010000114
Figure BDA0001945046010000115
in the formula, TλIs a hard threshold function;
Figure BDA0001945046010000116
and phi1A pair of positive and negative conversion is formed,
Figure BDA0001945046010000117
and phi2A pair of positive and negative conversion is formed;
termination conditions were as follows: when in use
Figure BDA0001945046010000118
When the value is smaller than the preset value, the influence of continuous iteration on the result is small enough, and the iteration is terminated;
and (3) outputting:
Figure BDA0001945046010000119
Figure BDA00019450460100001110
for the final transform coefficients of the separated signal component 1,
Figure BDA00019450460100001111
the final transform coefficients for the separated signal component 2;
in the block coordinate relaxation algorithm, the hard threshold function formula is as follows:
Figure BDA00019450460100001112
in the formula:
Figure BDA00019450460100001113
is a hard threshold function, and λ is a hard threshold, usually taken
Figure BDA00019450460100001114
In a large to small arrangement
Figure BDA00019450460100001115
The value of the one or more of,
Figure BDA00019450460100001116
is a coefficient matrix
Figure BDA00019450460100001117
N, N is the size of the coefficient matrix, Φ is the transform dictionary, and s is the time-domain signal.
As shown in fig. 6, the determining the low cutoff frequency of the reconstructed transform of the two-dimensional local discrete cosine transform in step 103 mainly includes the following steps:
step 301: taking M isochronous slices of the three-dimensional seismic data volume according to the interval time delta T and using the slices to test the low-frequency-cutoff parameters of the reconstruction transformation of the local cosine transformation, wherein the positions of the slices are T0,tΔT,t2ΔT,…,t(M-1)ΔT；
Step 302: for position tkΔTThe k is 0,1, M-1, a series of low-frequency-cut parameters of reconstruction transformation of two-dimensional local discrete cosine transformation are given, and the optimization problem is solved for each low-frequency-cut parameter to obtain separation of effective wave signals and collected footprint noise; by comparing the signal-to-noise separation effect under various parametersObtaining a low-frequency-cutting parameter with the best signal-noise separation effect at the current time position;
step 303: by linear interpolation, according to the time position tkΔTK-0, 1.., M-1 and its corresponding low-cutoff parameter fkΔTAnd k is 0,1,.. times, M-1, and low-frequency-cutoff parameters of reconstruction transformation of the two-dimensional local discrete cosine transformation of other time positions of the whole three-dimensional seismic data volume are obtained.
In the embodiment, the acquisition footprint pressing method based on morphological component analysis is provided, so that the purpose of effectively pressing the acquisition footprints in the three-dimensional seismic data is achieved. The invention has the following beneficial effects:
the method provided by the invention is used for carrying out acquisition footprint pressing on the three-dimensional seismic data volume, so that the acquisition footprints can be effectively pressed, part of random noise and offset processing false images are removed, and effective signals have higher fidelity.
The above transformation structure and parameter determination method are specifically described below with reference to a specific embodiment, and the embodiment is only for better illustrating the present invention and does not limit the present invention.
FIG. 7A is a main survey line section of marine three-dimensional seismic data of an oil field, where the data is collected at an earlier time, and the data collection quality is greatly different between different survey lines and the far cable collection quality is not good due to the limitation of the current collection technology and collection instrument, so that the offset data is affected by the collection footprint very seriously; in addition, in the marine seismic acquisition process, the temperature and salinity of seawater are greatly changed due to seasonal alternation and ocean current influence, so that the propagation speed of seismic waves is easily changed, and obvious acquisition footprints are generated in marine seismic data.
As shown in fig. 7B and 7C, the seismic record (fig. 7A) containing the acquisition footprint is processed by the method of the present invention to obtain a corresponding effective signal and removed acquisition footprint noise. The method can effectively suppress the collected footprint noise, also removes partial random noise and offset processing false images, and ensures that the processed data becomes clearer and the continuity of the same phase axis is stronger, thereby proving that the method has better effective wave signal fidelity.
Fig. 8A shows a 1.7s time slice of the three-dimensional seismic data volume in which the acquisition footprints are predominantly along the inline direction and are so pronounced that it is difficult to determine from the original time slice whether such amplitude anomalies are caused by lateral variations in the subsurface medium or the effects of the acquisition footprints.
As shown in fig. 8B-8C, as a result of processing the three-dimensional data volume by applying the method of the present invention, the method of the present invention suppresses the footprint noise well and maintains various fine structural features of the effective wave signal in the cross section well.
In the model and actual data calculation example, the sparse representation transformation of the acquired footprint is constructed by using the method of the invention, and the acquired footprint suppression is carried out on the seismic data, so that the acquired footprint noise can be effectively suppressed, and effective signals have higher fidelity, thereby laying a foundation for the analysis of subsequent data.
Finally, it should be noted that the above models and practical data examples provide further verification for the purpose, technical solution and beneficial effects of the present invention, which are only specific examples of the present invention and are not intended to limit the scope of the present invention. Any modification, improvement or equivalent replacement made within the spirit and principle of the present invention shall fall within the protection scope of the present invention.

Claims (3)

1. A collection footprint pressing method based on morphological component analysis is characterized by comprising the following steps:
step 101: constructing two-dimensional local discrete cosine transform according to the acquired footprint waveform morphological characteristics in the seismic record of the isochronous slice of the three-dimensional seismic data volume, combining the two-dimensional local discrete cosine transform with the two-dimensional stationary wavelet transform to form a super-complete dictionary, and simultaneously determining the decomposition series J of the two-dimensional stationary wavelet transform and the analysis window size W of the two-dimensional local discrete cosine transform;
step 102: performing preliminary pressing of the collected footprints in each isochronous slice by using a morphological component analysis-based method in a layer-by-layer execution mode on the original data of the three-dimensional seismic data volume;
step 103: determining a low cut-off frequency of the two-dimensional local discrete cosine transform;
step 104: repeating the step 102 and the step 103 until all the slice data are processed, performing signal-noise separation by using a method based on morphological component analysis, and suppressing the acquired footprint noise to finally realize the suppression of the acquired footprint noise in the three-dimensional seismic data volume;
step 101, comprising:
the object of applying morphological component analysis is to contain two components with different morphological characteristics:
s＝s1+s2，
in the formula: s represents a signal to be analyzed; s1、s2Representing two components of the signal having different morphological characteristics; respectively extract s1、s2These two components are the targets of morphological component analysis; suppose s1And s2Can be respectively composed of a dictionary phi1And phi2Efficient sparse representation, but with phi2Sparse representation s1And by phi1Sparse representation s2Poor temporal sparsity;
selecting a two-dimensional stationary wavelet transform as a transformation dictionary for sparsely representing effective signal components, wherein the two-dimensional stationary wavelet forward transform is:
Figure FDA0002502811530000011
Figure FDA0002502811530000012
Figure FDA0002502811530000013
Figure FDA0002502811530000014
in the formula HjAnd GjFilter banks respectively representing the j-th layer decomposition;
the inverse transformation of the two-dimensional stationary wavelet transform is:
Figure FDA0002502811530000021
in the formula
Figure FDA0002502811530000022
And
Figure FDA0002502811530000023
each represents HjAnd GjThe dual filter bank of (1);
selecting a two-dimensional local discrete cosine transform as a transformation dictionary for sparsely representing the collected footprints, wherein a forward transform of the two-dimensional local discrete cosine transform is:
Figure FDA0002502811530000024
wherein f (i, j) represents the signal to be analyzed,
Figure FDA0002502811530000025
two-dimensional local discrete cosine transform coefficients, l, representing a signal to be analyzed1,l2N-1, N being the total number of rows or columns of the signal f to be analyzed;
the inverse of the two-dimensional local discrete cosine transform is:
Figure FDA0002502811530000026
according to morphological component analysis theory, using the above selected dictionary phi1And phi2Forming an ultra-complete dictionary, sparsely representing a signal s, and calculating a sparse representation coefficient:
Figure FDA0002502811530000027
in the formula: x is the number of1To reconstruct the sum phi of the coefficients1A corresponding portion; x is the number of2To reconstruct the sum phi of the coefficients2A corresponding portion;
Figure FDA0002502811530000028
is a lagrange multiplier; the optimization problem is solved through a block coordinate relaxation algorithm;
step 102 specifically includes:
the block coordinate relaxation algorithm comprises the following steps:
initialization: initial iteration step number k is 0, initial solution
Figure FDA0002502811530000031
Figure FDA0002502811530000032
The initial solution of coefficients used to represent signal component 1,
Figure FDA0002502811530000033
the initial solution of coefficients used to represent signal component 2;
iteration: and increasing the iteration step number k by 1 in each step of iteration, and calculating:
Figure FDA0002502811530000034
Figure FDA0002502811530000035
in the formula, TλIs a hard threshold function;
Figure FDA0002502811530000036
and phi1A pair of positive and negative conversion is formed,
Figure FDA0002502811530000037
and phi2A pair of positive and negative conversion is formed;
termination conditions were as follows: when in use
Figure FDA0002502811530000038
When the value is smaller than the preset value, the influence of continuous iteration on the result is small enough, and the iteration is terminated;
and (3) outputting:
Figure FDA0002502811530000039
Figure FDA00025028115300000310
for the final transform coefficients of the separated signal component 1,
Figure FDA00025028115300000311
the final transform coefficients for the separated signal component 2;
in the block coordinate relaxation algorithm, the hard threshold function formula is as follows:
Figure FDA00025028115300000312
in the formula:
Figure FDA00025028115300000313
is composed of
Figure FDA00025028115300000314
Is a hard threshold, λ is a hard threshold,
Figure FDA00025028115300000315
is a coefficient matrix
Figure FDA00025028115300000316
N, N is the total number of rows or columns of the signal f to be analyzed, Φ is the transformation dictionary, and s is the time domain signal;
step 103 specifically comprises:
step 301: taking M isochronous slices of the three-dimensional seismic data volume according to the interval time delta T and using the slices to test the low-frequency-cutoff parameters of the reconstruction transformation of the local discrete cosine transformation, wherein the positions of the slices are T0,tΔT,t2ΔT,…,t(M-1)ΔT；
Step 302: for position thΔTH-0, 1., M-1, a series of low-frequency-cut parameters of reconstruction transformation of two-dimensional local discrete cosine transformation are given, and the separation of obtaining effective wave signals and collecting footprint noises is solved for each low-frequency-cut parameter; determining a low-cut-frequency parameter f for obtaining the best signal-noise separation effect at the current time position by comparing the signal-noise separation effects under the parametershΔT；
Step 303: by linear interpolation, according to the time position thΔTH-0, 1.., M-1 and its corresponding low-cutoff parameter fhΔTAnd h is 0,1, M-1, and obtaining low-frequency-cutoff parameters of reconstruction transformation of the two-dimensional local discrete cosine transformation of other time positions of the whole three-dimensional seismic data body.
2. The method according to claim 1, wherein the preset value is 10-6。
3. The method of claim 1, wherein the hard threshold λ is selected from the group consisting of
Figure FDA0002502811530000041
In a large to small arrangement
Figure FDA0002502811530000042
A value.
CN201910033429.3A 2019-01-14 2019-01-14 Acquisition footprint pressing method based on morphological component analysis Active CN109782346B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910033429.3A CN109782346B (en) 2019-01-14 2019-01-14 Acquisition footprint pressing method based on morphological component analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910033429.3A CN109782346B (en) 2019-01-14 2019-01-14 Acquisition footprint pressing method based on morphological component analysis

Publications (2)

Publication Number Publication Date
CN109782346A CN109782346A (en) 2019-05-21
CN109782346B true CN109782346B (en) 2020-07-28

Family

ID=66500652

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910033429.3A Active CN109782346B (en) 2019-01-14 2019-01-14 Acquisition footprint pressing method based on morphological component analysis

Country Status (1)

Country Link
CN (1) CN109782346B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112329595B (en) * 2020-11-02 2022-11-15 中南大学 Rock joint surface geometric morphology frequency spectrum analysis and reconstruction method

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8942467B2 (en) * 2012-03-23 2015-01-27 Mitsubishi Electric Research Laboratories, Inc. Method for reducing blocking artifacts in images
CN105182417B (en) * 2015-09-11 2018-06-22 合肥工业大学 A kind of surface wave separation method and system based on anatomic element analysis
CN107356967B (en) * 2017-07-26 2019-04-12 西安交通大学 A kind of compacting seismic data shields by force the sparse optimization method of interference

Also Published As

Publication number Publication date
CN109782346A (en) 2019-05-21

Similar Documents

Publication Publication Date Title
Zhu et al. Seismic data denoising through multiscale and sparsity-promoting dictionary learning
Gan et al. Compressive sensing for seismic data reconstruction via fast projection onto convex sets based on seislet transform
Chen et al. Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization
Xue et al. Amplitude-preserving iterative deblending of simultaneous source seismic data using high-order radon transform
Herrmann et al. Sparsity-and continuity-promoting seismic image recovery with curvelet frames
CN107817527B (en) Seismic exploration in desert stochastic noise suppression method based on the sparse compressed sensing of block
Yuan et al. Inversion-based 3-D seismic denoising for exploring spatial edges and spatio-temporal signal redundancy
Zhou et al. Spike-like blending noise attenuation using structural low-rank decomposition
CN110780348B (en) Primary wave and multiple combined imaging method and system based on stereo imaging conditions
Li et al. Wavelet-based higher order correlative stacking for seismic data denoising in the curvelet domain
CN109738950B (en) The noisy-type data primary wave inversion method of domain inverting is focused based on sparse 3 D
CN107179550B (en) A kind of seismic signal zero phase deconvolution method of data-driven
CN110031899A (en) Compressed sensing based weak signal extraction algorithm
WO2011103297A2 (en) Estimating internal multiples in seismic data
Zhang et al. 3D simultaneous seismic data reconstruction and noise suppression based on the curvelet transform
Zu et al. 3D deblending of simultaneous source data based on 3D multi-scale shaping operator
Zhou et al. A hybrid method for noise suppression using variational mode decomposition and singular spectrum analysis
Geng et al. Relative time seislet transform
Bai et al. Least-squares Gaussian beam transform for seismic noise attenuation
CN109782346B (en) Acquisition footprint pressing method based on morphological component analysis
Cheng et al. Deblending of simultaneous-source seismic data using Bregman iterative shaping
Guo et al. Seismic data denoising under the morphological component analysis framework by dictionary learning
CN102509268B (en) Immune-clonal-selection-based nonsubsampled contourlet domain image denoising method
WANG et al. Inversion of ground‐penetrating radar data for 2D electric parameters
Lv et al. Learning dictionary in the approximately flattened structure domain

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