CN117607966B - Two-dimensional mode decomposition based on curvelet domain method for enhancing weak signal of seismic data - Google Patents
Two-dimensional mode decomposition based on curvelet domain method for enhancing weak signal of seismic data Download PDFInfo
- Publication number
- CN117607966B CN117607966B CN202410097182.2A CN202410097182A CN117607966B CN 117607966 B CN117607966 B CN 117607966B CN 202410097182 A CN202410097182 A CN 202410097182A CN 117607966 B CN117607966 B CN 117607966B
- Authority
- CN
- China
- Prior art keywords
- curvelet
- seismic data
- signal
- coefficient
- scale
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 162
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 38
- 230000002708 enhancing effect Effects 0.000 title claims abstract description 7
- 230000009466 transformation Effects 0.000 claims abstract description 18
- 238000004458 analytical method Methods 0.000 claims abstract description 13
- 238000006243 chemical reaction Methods 0.000 claims abstract description 10
- 230000000694 effects Effects 0.000 claims abstract description 10
- 230000008569 process Effects 0.000 claims description 50
- 239000011159 matrix material Substances 0.000 claims description 34
- 238000006073 displacement reaction Methods 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000010219 correlation analysis Methods 0.000 claims description 3
- 230000002146 bilateral effect Effects 0.000 description 9
- 238000001914 filtration Methods 0.000 description 9
- 238000005452 bending Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
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. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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. analysis, for interpretation, for correction
- G01V1/32—Transforming one recording into another or one representation into another
-
- 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. analysis, for interpretation, for correction
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Abstract
The invention discloses a two-dimensional method based on curvelet domain A method for enhancing weak signals of seismic data by mode decomposition, in particular to the technical field of geophysics. The invention obtains the original curvelet coefficient and carries out scale analysis by carrying out curvelet forward conversion on the original seismic data, after the coarse scale signal curvelet coefficient, the middle scale signal curvelet coefficient and the fine scale signal curvelet coefficient are determined, hard threshold processing and weak effective signal compensation processing are carried out on the middle scale signal curvelet coefficient, and performing hard threshold processing on the small-scale signal curvelet coefficients, recombining the processed curvelet coefficients of each scale signal to obtain new curvelet coefficients, performing curvelet inverse transformation, obtaining the seismic data after denoising and weak signal enhancement, and evaluating the denoising effect. The method solves the problem that the curved wave transformation hard threshold method is discontinuous at the threshold, denoising the seismic data, guaranteeing the continuity of weak and same-phase axes of the seismic data, compensating the weak signal energy and effectively improving the quality of seismic data processing.
Description
Technical Field
The invention relates to the technical field of geophysics, in particular to a seismic data weak signal enhancement method based on curvelet domain two-dimensional mode decomposition.
Background
High precision exploration requires processing of seismic data with higher fidelity, better identification of weak effective signals and maximum utilization of low frequency information. Most seismic data processing methods cannot extract signals according to time sequence features of the signals, particularly deep data, and the complexity of the deep stratum data is higher due to interference of noise in the data acquisition process. When the existing denoising method is adopted to process the seismic data, weak signals in the seismic data are more blurred and even disappear, original characteristics of the seismic data signals are destroyed, and accordingly quality of the seismic data is reduced. Therefore, improving the resolution of seismic data and compensating for weak effective signal energy has become one of the difficulties and important points in seismic data processing.
At present, as the curvelet hard threshold method is discontinuous at the threshold, partial weak effective signals in the seismic data imaging are blurred, the continuity of the signals is reduced, and the resolution of the seismic data is affected. The method for combining the curvelet domain and the bilateral filtering considers the space domain and the value domain simultaneously in the denoising process, and effectively reserves the detailed information such as weak effective signals and the like while realizing image denoising. However, the method of combining the curvelet domain and the bilateral filtering lacks of self-adaptability, so that artifacts appear after denoising.
Disclosure of Invention
Aiming at the problem that the curvelet domain combined bilateral filtering method lacks self-adaptability, the invention provides a seismic data weak signal enhancement method based on curvelet domain two-dimensional mode decomposition, and the problem that a curvelet transformation hard threshold method is discontinuous at a threshold value is solved by introducing the two-dimensional mode decomposition into the curvelet domain, so that the denoising of the seismic data is realized, meanwhile, the continuity of a weak phase axis of the seismic data is improved, the weak signal energy is compensated, and the quality of seismic data processing is effectively improved.
In order to achieve the above purpose, the invention adopts the following technical scheme:
a seismic data weak signal enhancement method based on curvelet domain two-dimensional mode decomposition comprises the following steps:
step 1, acquiring seismic data, and performing curved wave forward conversion processing on the seismic data to obtain an original curved wave coefficient;
step 2, performing scale analysis on the original curvelet coefficient, and dividing the original curvelet coefficient into a coarse scale signal curvelet coefficient, a middle scale signal curvelet coefficient and a fine scale signal curvelet coefficient;
step 3, performing hard threshold processing and weak effective signal compensation processing on the curvelet coefficient of the intermediate scale signal to obtain a processed curvelet coefficient and using the processed curvelet coefficient as a new intermediate scale signal;
step 4, performing hard threshold processing on the curve wave coefficient of the fine-scale signal to obtain a processed fine-scale signal serving as a new fine-scale signal;
step 5, recombining the curvelet coefficient of the coarse-scale signal, the new intermediate-scale signal and the new fine-scale signal to obtain a new curvelet coefficient;
step 6, performing inverse curvelet transformation on the new curvelet coefficient to obtain the seismic data after denoising and weak signal enhancement;
and 7, evaluating the denoising effect of the processed seismic data based on the signal-to-noise ratio and the peak signal-to-noise ratio of the seismic data.
Preferably, the seismic data includes effective signals and noise, as shown in formula (1):
（1）
in the method, in the process of the invention,in order to be able to take time,in the case of noisy seismic data,is the original seismic data that is free of noise,is noise;
according to the number of seismic channels and the number of sampling points of the acquired noise-containing seismic data, determining the total scale number as follows:
（2）
in the method, in the process of the invention,as the number of degrees of the total scale,is the number of seismic traces of noisy seismic data,the number of samples for noisy seismic data,in order to perform the rounding operation upwards,as a function of the minimum value;
fourier transform results due to mother curved waveEqual to the fourier window function valueAnd combining the displacement factor sequence, the rotation angle and the position, and translating and rotating the mother curvelet to obtain a curvelet function, wherein the curvelet function is shown in a formula (3):
（3）
wherein,
rotation matrixThe method comprises the following steps:
（4）
in the method, in the process of the invention,is a spatial domain factor;as a scale factor of the dimensions of the device,，is an integer of the number of the times,for the reason of scaleMaximum value of the seed;as a factor of the angle,64，is an integer;for the sequence of displacement factors,，for the first sequence of displacement factors,for the second sequence of displacement factors,is an integer set;is a curvelet function;is a mother curvelet;is the position;in order for the angle of rotation to be a function of,，，is a downward rounding operation;
based on the curvelet function, in real number domainPerforming curved wave forward transformation on the acquired noise-containing seismic data to obtain original curved wave coefficientsThe method comprises the following steps:
（5）
in the method, in the process of the invention,is the fourier transform of noisy seismic data,as a function of the fourier window,as a factor of the frequency domain,is imaginary.
Preferably, in the step 2, the curved wave coefficients corresponding to each scale of the seismic data are respectively subjected to curved wave inverse transformation to obtain a spectrogram of each scale of information;
analyzing a spectrogram of each scale information, and obtaining a coarse scale signal, a middle scale signal and a fine scale signal through scale analysis, wherein the coarse scale signal only contains low-frequency information, the middle scale signal contains a large number of effective signals, and the fine scale signal contains high-frequency noise signals;
dividing the original curvelet coefficient into coarse-scale signal curvelet coefficients based on the scale analysis resultMesoscale signal curvelet coefficientAnd a fine scale signal curvelet coefficientAs shown in formula (6):
（6）。
preferably, the hard thresholding preserves significant strong signals during mid-scale signal wavelet coefficientsPerforming hard threshold denoising, wherein the curvelet function reserved by each scale information is not more than 20% during denoising, retaining effective strong signals and removing noise signals to obtain denoised curvelet coefficientsThe method comprises the following steps:
（7）
in the method, in the process of the invention,a threshold value set when denoising the curvelet coefficient of the intermediate scale signal;
in the weak effective signal compensation processing process, the number of the curvelet coefficients reserved by each scale information is controlled to be 30% -50% according to the decomposition condition of a two-dimensional empirical mode decomposition method and the factors for reducing noise, and a hard threshold coefficient matrix is constructedAs shown in formula (8):
（8）
in the method, in the process of the invention,is a coefficient matrix threshold.
Preferably, the number of the two-dimensional empirical mode decomposition is set for the hard threshold coefficient matrixAnd (3) performing two-dimensional empirical mode decomposition to obtain:
（9）
in the method, in the process of the invention,to decompose intoThe corresponding modal function matrix is then used,as a result of the residual signal,in order to decompose the number of the components,number of decompositionCorresponding to different modes;
respectively carrying out matrix and denoised curvelet coefficient on each mode functionPerforming cross-correlation analysis and calculating to obtain cross-correlation coefficientsAs shown in formula (10):
（10）
in the method, in the process of the invention,is the mean value of each modal function matrix,as a matrix of the modal functions,for the denoised curvelet coefficient matrix,the mean value of the denoised curvelet coefficient matrix;
acquiring an eigen mode function corresponding to the maximum cross correlation coefficient in each angle to obtain:
（11）
in the method, in the process of the invention,for the number of sequences of eigenmode functions,the function is used for solving the eigenvalue function corresponding to the maximum cross-correlation coefficient;
for denoised curvelet coefficientCompensating weak effective signal to obtain compensated curvelet coefficientAs shown in formula (12):
（12）
in the method, in the process of the invention,is the eigenmode function corresponding to the maximum cross-correlation coefficient.
Preferably, the fine-scale signal bending coefficient is hard thresholded as shown in equation (13):
（13）
in the method, in the process of the invention,is the wavelet coefficient of the fine-scale signal after hard thresholding.
Preferably, the new curvelet coefficientThe method comprises the following steps:
（14）
performing inverse wavelet transform on the new wavelet coefficients to obtain denoised and weak signal enhanced seismic data, as shown in formula (15):
（15）
in the method, in the process of the invention,is the seismic data after denoising and weak signal enhancement.
Preferably, in the step 7, a signal-to-noise ratio calculation formula of the seismic data is:
（16）
in the method, in the process of the invention,is seismic dataIs set in the number of rows of (a),and is also provided withIs an integer of the number of the times,is the total number of line numbers;is seismic dataIs used for the sequence number of (c),and is also provided withIs an integer of the number of the times,is the total number of column ordinals;is seismic dataIs used for the signal-to-noise ratio of (c),is seismic dataIs a denoising result;
the calculation formula of the peak signal-to-noise ratio is as follows:
（17）
wherein,
（18）
in the method, in the process of the invention,is seismic dataIs used for the peak signal-to-noise ratio of (c) in the (c),as a function of the maximum value,is the variance of the seismic data.
The beneficial technical effects brought by the invention are as follows:
the invention provides a seismic data weak signal enhancement method based on a curved wave domain two-dimensional mode decomposition, which is characterized in that an original curved wave coefficient is obtained by performing curved wave forward conversion processing on an acquired original seismic, after the original curved wave coefficient is subjected to scale analysis, a coarse scale signal curved wave coefficient, an intermediate scale signal curved wave coefficient and a fine scale signal curved wave coefficient are determined, the two-dimensional mode decomposition is introduced into a curved wave domain, the intermediate scale signal curved wave coefficient is processed by the two-dimensional mode decomposition method, and meanwhile, the fine scale signal curved wave coefficient is subjected to hard threshold processing, the processed curved wave coefficients after the weak signal compensation are obtained by recombining the processed scale signal coefficients, and the seismic data after the weak signal compensation is obtained by inverse curved wave conversion.
The method applies the two-dimensional empirical mode decomposition method to the curvelet transformation domain, has the characteristics of strong self-adaptability and high effective signal purity, and can better nonlinear and nonstationary signals. The method effectively solves the problem that the curvelet transformation hard threshold method is discontinuous at the threshold value, improves the continuity of weak in-phase axes of the seismic data and compensates the weak signal energy while denoising the seismic data, is beneficial to improving the quality of seismic data processing, reduces the artifact, enriches the deep weak signal information of the stratum, and provides high-quality seismic data for model construction and interpretation analysis of the stratum deep seismic data.
Drawings
FIG. 1 is a flow chart of a method of seismic data weak signal enhancement based on curvelet domain two-dimensional pattern decomposition.
Fig. 2 is a flow chart of signal processing for each scale in the original curvelet coefficient scale.
Detailed Description
The present invention will be described in further detail with reference to the accompanying drawings and examples.
Example 1
The embodiment discloses a seismic data weak signal enhancement method based on curvelet domain two-dimensional mode decomposition, which specifically comprises the following steps as shown in fig. 1:
step 1, acquiring seismic data, and performing curved wave forward conversion processing on the seismic data to obtain an original curved wave coefficient.
The seismic data includes a valid signal and noise as shown in equation (1):
（1）
in the method, in the process of the invention,in order to be able to take time,in the case of noisy seismic data,is the original seismic data that is free of noise,is noise.
According to the number of seismic channels and the number of sampling points of the acquired noise-containing seismic data, determining the total scale number as follows:
（2）
in the method, in the process of the invention,as the number of degrees of the total scale,is the number of seismic traces of noisy seismic data,the number of samples for noisy seismic data,in order to perform the rounding operation upwards,as a function of the minimum value.
Fourier transform results due to mother curved waveEqual to the Fourier window functionAnd combining the displacement factor sequence, the rotation angle and the position, and translating and rotating the mother curvelet to obtain a curvelet function, wherein the curvelet function is shown in a formula (3):
（3）
wherein,
rotation matrixThe method comprises the following steps:
（4）
in the method, in the process of the invention,is a spatial domain factor;as a scale factor of the dimensions of the device,，is an integer of the number of the times,is the maximum value of the scale factor;as a factor of the angle,64，is an integer;for the sequence of displacement factors,，for the first sequence of displacement factors,for the second sequence of displacement factors,is an integer set;is a curvelet function;is a mother curvelet;is the position;in order for the angle of rotation to be a function of,，，is a round-down operation.
Based on the curvelet function, the collected noise-containing seismic data is subjected to curvelet forward conversion in a real number domain to obtain an original curvelet coefficientThe method comprises the following steps:
（5）
in the method, in the process of the invention,is the fourier transform of noisy seismic data,as a function of the fourier window,as a factor of the frequency domain,is imaginary.
And 2, performing scale analysis on the original curvelet coefficient, and dividing the original curvelet coefficient into a coarse-scale signal curvelet coefficient, an intermediate-scale signal curvelet coefficient and a fine-scale signal curvelet coefficient, as shown in fig. 2.
In the step 2, the inverse curvelet transformation is respectively carried out on curvelet coefficients corresponding to each scale of the seismic data, so as to obtain a spectrogram of each scale of information; and analyzing a spectrogram of each scale information, and obtaining a coarse scale signal, a middle scale signal and a fine scale signal by scale analysis, wherein the coarse scale signal only contains low-frequency information, the middle scale signal contains a large number of effective signals, and the fine scale signal contains high-frequency noise signals.
Since the coarse-scale signal has only one scale and one angle and contains fewer curved wave coefficients, the coarse-scale signal is not processed.
Dividing the original curvelet coefficient into coarse-scale signal curvelet coefficients based on the scale analysis resultMesoscale signal curvelet coefficientAnd a fine scale signal curvelet coefficientAs shown in formula (6):
（6）。
and step 3, performing hard threshold processing and weak effective signal compensation processing for retaining effective strong signals on the curvelet coefficients of the intermediate scale signals to obtain processed curvelet coefficients which are used as new intermediate scale signals.
In the process of preserving effective strong signals in the hard thresholding, the medium-scale signal curvelet coefficientsPerforming hard threshold denoising, wherein the curvelet function reserved by each scale information is not more than 20% during denoising, retaining effective strong signals and removing noise signals to obtain denoised curvelet coefficientsThe method comprises the following steps:
（7）
in the method, in the process of the invention,and (5) a threshold value set for denoising the curved wave coefficient of the intermediate scale signal.
In the process of the weak effective signal compensation processing, the quantity of the curvelet coefficients reserved by each scale information is controlled to be 30% -50% according to the decomposition condition of the two-dimensional empirical mode decomposition method and the factors for reducing noise, and if the small curvelet coefficients are not removed, the magnitude of the curvelet coefficients of the effective signal and the magnitude of the curvelet coefficients of the noise signal are increased simultaneously after the two-dimensional empirical mode decomposition processing, so that the display effect of the weak effective signal is disturbed.
Construction of a hard threshold coefficient matrixAs shown in formula (8):
（8）
in the method, in the process of the invention,is a coefficient matrix threshold.
Setting the number of the two-dimensional empirical mode decomposition and for a hard threshold coefficient matrixAnd (3) performing two-dimensional empirical mode decomposition to obtain:
（9）
in the method, in the process of the invention,to decompose intoThe corresponding modal function matrix is then used,as a result of the residual signal,in order to decompose the number of the components,number of decompositionCorresponding to different modalities.
Respectively carrying out matrix and denoised curvelet coefficient on each mode functionPerforming cross-correlation analysis and calculating to obtain cross-correlation coefficientsAs shown in formula (10):
（10）
in the method, in the process of the invention,is the mean value of each modal function matrix,as a matrix of the modal functions,for the denoised curvelet coefficient matrix,is the mean value of the denoised curvelet coefficient matrix.
Acquiring an eigen mode function corresponding to the maximum cross correlation coefficient in each angle to obtain:
（11）
in the method, in the process of the invention,for the number of sequences of eigenmode functions,is a function of the eigenmode function corresponding to the maximum cross correlation coefficient.
For denoised curvelet coefficientCompensating weak effective signal to obtain compensated curvelet coefficientAs shown in formula (12):
（12）
in the method, in the process of the invention,is the eigenmode function corresponding to the maximum cross-correlation coefficient.
Step 4, hard threshold processing is carried out on the bending wave coefficient of the fine-scale signal, the processed fine-scale signal is obtained and used as a new fine-scale signal, and the method is as shown in a formula (13):
（13）
in the method, in the process of the invention,is the wavelet coefficient of the fine-scale signal after hard thresholding.
Step 5, recombining the curvelet coefficient of the coarse-scale signal, the new intermediate-scale signal and the new fine-scale signal to obtain a new curvelet coefficient, wherein the new curvelet coefficientThe method comprises the following steps:
（14）。
step 6, performing inverse curvelet transformation on the new curvelet coefficient to obtain the seismic data after denoising and weak signal enhancement, as shown in a formula (15):
（15）
in the method, in the process of the invention,is the seismic data after denoising and weak signal enhancement.
And 7, evaluating the denoising effect of the processed seismic data based on the signal-to-noise ratio and the peak signal-to-noise ratio of the seismic data, wherein the denoising effect of the seismic data is better when the signal-to-noise ratio is higher, and the denoising effect of the seismic data is better when the peak signal-to-noise ratio is higher.
In the step 7, the signal-to-noise ratio calculation formula of the seismic data is:
（16）
in the method, in the process of the invention,is seismic dataIs set in the number of rows of (a),and is also provided withIs an integer of the number of the times,is the total number of line numbers;is seismic dataIs used for the sequence number of (c),and is also provided withIs an integer of the number of the times,is the total number of column ordinals;is seismic dataIs used for the signal-to-noise ratio of (c),is seismic dataIs a denoising result;
the calculation formula of the peak signal-to-noise ratio is as follows:
（17）
wherein,
（18）
in the method, in the process of the invention,is seismic dataIs used for the peak signal-to-noise ratio of (c) in the (c),as a function of the maximum value,is the variance of the seismic data.
Example 2
The method for enhancing the weak seismic data signals based on the curvelet-domain two-dimensional mode decomposition, which is provided in the embodiment 1, is applied to actual weak seismic data signal enhancement processing, and specifically comprises the following steps:
step 1, collecting seismic data, wherein the number of sampling points is 3501, sampling intervals are 2ms, and performing curved wave forward conversion processing on the collected seismic data to obtain an original curved wave coefficient, wherein the number of the collected original seismic data is 2001 in total.
And 2, performing scale analysis on the original curvelet coefficient, and dividing the original curvelet coefficient into a coarse scale signal curvelet coefficient, an intermediate scale signal curvelet coefficient and a fine scale signal curvelet coefficient.
And 3, performing hard threshold processing and weak effective signal compensation processing for retaining effective strong signals on the medium-scale signal curvelet coefficients, setting different thresholds on the medium-scale signal curvelet coefficients of each scale, and retaining only the curvelet coefficients larger than the thresholds under each scale to ensure the two-dimensional empirical mode decomposition to obtain a hard threshold coefficient matrix.
And carrying out two-dimensional empirical mode decomposition on the hard threshold coefficient matrix to obtain an intrinsic mode function, carrying out curved wave forward conversion on the selected actual seismic data, and then carrying out two-dimensional empirical mode decomposition method processing on the curved wave coefficients of the 2 nd scale and the 3 rd angle to obtain a decomposition result of the original curved wave coefficients.
And calculating cross-correlation coefficients between the eigenmode functions and the denoised curvelet coefficients respectively, selecting the eigenmode function corresponding to the maximum cross-correlation coefficient, and adding the eigenmode function and the denoised intermediate scale curvelet coefficient to obtain a new intermediate scale signal curvelet coefficient.
And step 4, performing hard threshold processing on the wavelet coefficients of the fine-scale signals to obtain new fine-scale signals.
And 5, integrating the curvelet coefficient of the coarse-scale signal, the new intermediate-scale signal and the new fine-scale signal to obtain a new curvelet coefficient.
And 6, performing inverse curvelet transformation on the new curvelet coefficient to obtain the seismic data after denoising and weak signal enhancement.
Comparing the denoising and weak signal enhancement effects of the method and the curvelet transformation hard threshold method and the curvelet transformation-combined bilateral filtering method, the method finds that weak effective signals in seismic data are removed while denoising after the curvelet transformation hard threshold method is adopted, and partial effective signals and resolution ratio reduction can be seen according to a denoising difference value; after being processed by a curvelet transformation-combined bilateral filtering method, the seismic data is denoised while weak effective signals are reserved, but some artifacts are reserved; when the method is used for processing the seismic data, the noise of the seismic data is removed, meanwhile, weak effective signals are reserved, compared with a curvelet transformation-combined bilateral filtering method, the method reduces artifacts, particularly enriches deep weak signal information, and provides high-quality seismic data for subsequent modeling and interpretation of deep data.
Comparing the original seismic data with the same channel of seismic data by adopting a curvelet transformation hard threshold method, a curvelet transformation-combined bilateral filtering method and the method, the method has the advantages that compared with the curvelet transformation hard threshold method and the curvelet transformation-combined bilateral filtering method, when the method is adopted to enhance the weak signals of the seismic data based on curvelet domain two-dimensional mode decomposition, the continuity of the weak effective signals is further improved by properly reducing the energy of the effective strong signals while removing the noise of the seismic data.
And 7, evaluating the denoising effect of the processed seismic data based on the signal-to-noise ratio and the peak signal-to-noise ratio of the seismic data.
The signal-to-noise ratio value and the peak signal-to-noise ratio value of the seismic data processed by the method are respectively calculated and compared by the method and the curvelet transformation hard threshold method and the curvelet transformation-combined bilateral filtering method, and the signal-to-noise ratio value and the peak signal-to-noise ratio value of the seismic data processed by the method are improved after comparison, and the signal-to-noise ratio value and the peak signal-to-noise ratio value processed by the method are the largest in the three methods.
Therefore, the method is obviously superior to the traditional curvelet domain processing method for processing actual seismic data, solves the problem that a hard threshold method is discontinuous at a threshold, improves the continuity of weak phase axes while removing the noise of the seismic data, and obviously improves the signal to noise ratio of the seismic data.
In the description of the present invention, it should be noted that, unless explicitly specified and limited otherwise, the terms "disposed," "mounted," "connected," and "fixed" are to be construed broadly, and may be, for example, fixedly connected, detachably connected, or integrally connected; can be mechanically or electrically connected; can be directly connected or indirectly connected through an intermediate medium, and can be communication between two elements. The specific meaning of the above terms in the present invention can be understood by those of ordinary skill in the art in a specific case.
It should be understood that the above description is not intended to limit the invention to the particular embodiments disclosed, but to limit the invention to the particular embodiments disclosed, and that the invention is not limited to the particular embodiments disclosed, but is intended to cover modifications, adaptations, additions and alternatives falling within the spirit and scope of the invention.
Claims (4)
1. The seismic data weak signal enhancement method based on curvelet domain two-dimensional mode decomposition is characterized by comprising the following steps of:
step 1, acquiring seismic data, and performing curved wave forward conversion processing on the seismic data to obtain an original curved wave coefficient;
step 2, performing scale analysis on the original curvelet coefficient, and dividing the original curvelet coefficient into a coarse scale signal curvelet coefficient, a middle scale signal curvelet coefficient and a fine scale signal curvelet coefficient;
step 3, performing hard threshold processing and weak effective signal compensation processing on the curvelet coefficient of the intermediate scale signal to obtain a processed curvelet coefficient and using the processed curvelet coefficient as a new intermediate scale signal;
step 4, performing hard threshold processing on the curve wave coefficient of the fine-scale signal to obtain a processed fine-scale signal serving as a new fine-scale signal;
step 5, recombining the curvelet coefficient of the coarse-scale signal, the new intermediate-scale signal and the new fine-scale signal to obtain a new curvelet coefficient;
step 6, performing inverse curvelet transformation on the new curvelet coefficient to obtain the seismic data after denoising and weak signal enhancement;
step 7, evaluating the denoising effect of the processed seismic data based on the signal-to-noise ratio and the peak signal-to-noise ratio of the seismic data;
the seismic data includes a valid signal and noise as shown in equation (1):
（1）
in the method, in the process of the invention,for time (I)>Is noise-containing seismic data->Is noise-free raw seismic data, +.>Is noise;
according to the number of seismic channels and the number of sampling points of the acquired noise-containing seismic data, determining the total scale number as follows:
（2）
in the method, in the process of the invention,for total degree of scale, +.>Number of seismic traces, which is noisy seismic data, +.>Sample points for noisy seismic data, +.>For the rounding up operation, ++>As a function of the minimum value;
fourier transform results due to mother curved waveEqual to the fourier window function value->And combining the displacement factor sequence, the rotation angle and the position, and translating and rotating the mother curvelet to obtain a curvelet function, wherein the curvelet function is shown in a formula (3):
（3）
wherein,
rotation matrixThe method comprises the following steps:
（4）
in the method, in the process of the invention,is a spatial domain factor; />As a scale factor, < >>，/>Is an integer>Is the total degree of scale; />Is an angle factor->64，/>Is an integer; />For the sequence of shift factors, +.>，/>For the first displacement factor sequence,/a.sub.f.>For the second shifting factor sequence,/a. About.>Is an integer set; />Is a curvelet function; />Is a mother curvelet;is the position; />For the rotation angle +.>，/>，/>Is a downward rounding operation;
based on the curvelet function, the collected noise-containing seismic data is subjected to curvelet forward conversion in a real number domain to obtain an original curvelet coefficientThe method comprises the following steps:
（5）
in the method, in the process of the invention,is the Fourier transform of noisy seismic data, +.>For Fourier window function, +.>Is a frequency domain factor, +.>Is an imaginary number;
in the step 2, the inverse curvelet transformation is respectively carried out on curvelet coefficients corresponding to each scale of the seismic data, so as to obtain a spectrogram of each scale of information;
analyzing a spectrogram of each scale information, and obtaining a coarse scale signal, a middle scale signal and a fine scale signal through scale analysis, wherein the coarse scale signal only contains low-frequency information, the middle scale signal contains a large number of effective signals, and the fine scale signal contains high-frequency noise signals;
dividing the original curvelet coefficient into coarse-scale signal curvelet coefficients based on the scale analysis resultMidrange signal curvelet coefficient->And the wavelet coefficient of the fine-scale signal +.>As shown in formula (6):
（6）；
in the process of preserving effective strong signals in the hard thresholding, the medium-scale signal curvelet coefficientsPerforming hard threshold denoising, wherein the curvelet function reserved by each scale information is not more than 20% during denoising, retaining effective strong signals and removing noise signals to obtain a denoised curvelet coefficient matrix->The method comprises the following steps:
（7）
in the method, in the process of the invention,a threshold value set when denoising the curvelet coefficient of the intermediate scale signal;
in the weak effective signal compensation processing process, the number of the curvelet coefficients reserved by each scale information is controlled to be 30% -50% according to the decomposition condition of a two-dimensional empirical mode decomposition method and the factors for reducing noise, and a hard threshold coefficient matrix is constructedAs shown in formula (8)：
（8）
In the method, in the process of the invention,is a coefficient matrix threshold;
setting the number of the two-dimensional empirical mode decomposition and for a hard threshold coefficient matrixAnd (3) performing two-dimensional empirical mode decomposition to obtain:
（9）
in the method, in the process of the invention,for decomposing number +.>Modal function matrix corresponding to time, +.>For the residual signal>For decomposing number of parts>Decomposing number->Corresponding to different modes;
respectively carrying out the mode function matrix and the denoised curvelet coefficient matrixPerforming cross-correlation analysis, and calculating to obtain cross-correlation coefficient +.>As shown in formula (10):
（10）
in the method, in the process of the invention,for the mean value of the mode function matrix +.>Is a modal function matrix>Is a denoised curvelet coefficient matrix, < ->The mean value of the denoised curvelet coefficient matrix;
acquiring an eigen mode function corresponding to the maximum cross correlation coefficient in each angle to obtain:
（11）
in the method, in the process of the invention,for the number of eigenmode functions, +.>The function is used for solving the eigenvalue function corresponding to the maximum cross-correlation coefficient;
to the denoised curvelet coefficient matrixCompensating weak effective signal to obtain compensated curvelet coefficient +.>As shown in formula (12):
（12）
in the method, in the process of the invention,is the eigenmode function corresponding to the maximum cross-correlation coefficient.
2. The method for enhancing a weak signal of seismic data based on a curvelet domain two-dimensional pattern decomposition according to claim 1, wherein said fine-scale signal curvelet coefficient hard thresholding is as shown in formula (13):
（13）
in the method, in the process of the invention,is the wavelet coefficient of the fine-scale signal after hard thresholding.
3. The method for enhancing a weak seismic data signal based on a curvelet-domain two-dimensional pattern decomposition according to claim 2, wherein said new curvelet coefficientsThe method comprises the following steps:
（14）
performing inverse wavelet transform on the new wavelet coefficients to obtain denoised and weak signal enhanced seismic data, as shown in formula (15):
（15）
in the method, in the process of the invention,is the seismic data after denoising and weak signal enhancement.
4. The method for enhancing a weak signal of seismic data based on a curvelet-domain two-dimensional pattern decomposition according to claim 1, wherein in the step 7, a signal-to-noise ratio calculation formula of the seismic data is:
（16）
in the method, in the process of the invention,for seismic data->Is (are) the number of lines of->And->Is an integer>Is the total number of line numbers; />For seismic data->Ordinal number of->And->Is an integer>Is the total number of column ordinals; />For seismic data->Signal to noise ratio of>For seismic data->Is a denoising result;
the calculation formula of the peak signal-to-noise ratio is as follows:
（17）
wherein,
（18）
in the method, in the process of the invention,for seismic data->Peak signal to noise ratio,/->As a maximum function>Is the variance of the seismic data.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410097182.2A CN117607966B (en) | 2024-01-24 | 2024-01-24 | Two-dimensional mode decomposition based on curvelet domain method for enhancing weak signal of seismic data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410097182.2A CN117607966B (en) | 2024-01-24 | 2024-01-24 | Two-dimensional mode decomposition based on curvelet domain method for enhancing weak signal of seismic data |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117607966A CN117607966A (en) | 2024-02-27 |
| CN117607966B true CN117607966B (en) | 2024-04-09 |
Family
ID=89948381
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410097182.2A Active CN117607966B (en) | 2024-01-24 | 2024-01-24 | Two-dimensional mode decomposition based on curvelet domain method for enhancing weak signal of seismic data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117607966B (en) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009051899A1 (en) * | 2007-10-17 | 2009-04-23 | Exxonmobil Upstream Research Company | Method to adapt a template dataset to a target dataset by using curvelet representations |
| CN104007469A (en) * | 2014-05-24 | 2014-08-27 | 长江大学 | Weak seismic signal reconstruction method based on curvelet transform |
| CN110031899A (en) * | 2018-01-11 | 2019-07-19 | 中石化石油工程技术服务有限公司 | Compressed sensing based weak signal extraction algorithm |
| CN114779343A (en) * | 2022-05-05 | 2022-07-22 | 中国石油大学(华东) | Seismic data denoising method based on curvelet transform-joint bilateral filtering |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150120199A1 (en) * | 2013-10-31 | 2015-04-30 | Matthew S. Casey | Multiple Domain Processing For Combining Reservoir Models and Seismic Data |
| US20170212274A1 (en) * | 2015-08-21 | 2017-07-27 | Halliburton Energy Services, Inc. | Borehole Acoustic Logging Receiver Quality Control and Calibration |
-
2024
- 2024-01-24 CN CN202410097182.2A patent/CN117607966B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009051899A1 (en) * | 2007-10-17 | 2009-04-23 | Exxonmobil Upstream Research Company | Method to adapt a template dataset to a target dataset by using curvelet representations |
| CN104007469A (en) * | 2014-05-24 | 2014-08-27 | 长江大学 | Weak seismic signal reconstruction method based on curvelet transform |
| CN110031899A (en) * | 2018-01-11 | 2019-07-19 | 中石化石油工程技术服务有限公司 | Compressed sensing based weak signal extraction algorithm |
| CN114779343A (en) * | 2022-05-05 | 2022-07-22 | 中国石油大学(华东) | Seismic data denoising method based on curvelet transform-joint bilateral filtering |
Non-Patent Citations (2)
| Title |
|---|
| Curvelet threshold denoising joint with empirical mode decomposition;Lieqian Dong et al.;SEG Houston 2013 Annual Meeting;20131231;全文 * |
| 基于CEEMD的曲波阈值去噪方法;董烈乾 等;地球物理学进展;20161231;第31卷(第6期);第2521页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117607966A (en) | 2024-02-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Mohan et al. | A survey on the magnetic resonance image denoising methods | |
| Neelamani et al. | Wavelet-based deconvolution for ill-conditioned systems | |
| Mohan et al. | MRI denoising using nonlocal neutrosophic set approach of Wiener filtering | |
| CN110174702B (en) | Method and system for recovering low-frequency weak signals of marine seismic data | |
| Zhao et al. | 2-D TFPF based on Contourlet transform for seismic random noise attenuation | |
| CN111158051B (en) | Joint constraint random noise suppression method based on sparse regularization | |
| CN108961181A (en) | A kind of ground penetrating radar image denoising method based on shearlet transformation | |
| CN111242854A (en) | Image denoising method | |
| CN109724693B (en) | Fusion spectrum denoising method based on stationary wavelet | |
| CN112578471A (en) | Method for removing clutter noise of ground penetrating radar | |
| Hu et al. | The complex data denoising in MR images based on the directional extension for the undecimated wavelet transform | |
| Tang et al. | Adaptive threshold shearlet transform for surface microseismic data denoising | |
| CN114779343A (en) | Seismic data denoising method based on curvelet transform-joint bilateral filtering | |
| Zhang et al. | Seismic signal enhancement and noise suppression using structure-adaptive nonlinear complex diffusion | |
| CN117607966B (en) | Two-dimensional mode decomposition based on curvelet domain method for enhancing weak signal of seismic data | |
| Michailovich et al. | Iterative reconstruction of medical ultrasound images using spectrally constrained phase updates | |
| Routray et al. | MRI Denoising using sparse based Curvelet transform with variance stabilizing transformation framework | |
| CN116612032A (en) | Sonar image denoising method and device based on self-adaptive wiener filtering and 2D-VMD | |
| CN110032968A (en) | Denoising method based on dual-tree complex wavelet and adaptive semi-soft threshold model | |
| CN113160063A (en) | Method for denoising tunnel image based on wavelet and median filtering | |
| CN113484913B (en) | Seismic data denoising method for multi-granularity feature fusion convolutional neural network | |
| CN112363215B (en) | Seismic exploration signal enhancement method based on wave atom decomposition adaptive filtering | |
| CN114076986B (en) | Multi-scale dictionary learning sparse denoising method | |
| Bhuiyan et al. | A bidimensional empirical mode decomposition method for color image processing | |
| Wang et al. | An adaptive PolSAR trilateral filter based on the mechanism of scattering consistency |
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 |