CN104062683A - Combined attenuation random noise processing method based on curvelet transform and total variation - Google Patents
Combined attenuation random noise processing method based on curvelet transform and total variation Download PDFInfo
- Publication number
- CN104062683A CN104062683A CN201410107764.0A CN201410107764A CN104062683A CN 104062683 A CN104062683 A CN 104062683A CN 201410107764 A CN201410107764 A CN 201410107764A CN 104062683 A CN104062683 A CN 104062683A
- Authority
- CN
- China
- Prior art keywords
- data
- curvelet
- scale
- total variation
- threshold
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 8
- 238000000034 method Methods 0.000 claims abstract description 26
- 238000005516 engineering process Methods 0.000 claims abstract description 11
- 238000012545 processing Methods 0.000 claims description 20
- 230000009466 transformation Effects 0.000 claims description 20
- 230000005284 excitation Effects 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 8
- 238000012937 correction Methods 0.000 claims description 6
- 241000287196 Asthenes Species 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000004880 explosion Methods 0.000 claims description 3
- 239000002360 explosive Substances 0.000 claims description 3
- 230000010355 oscillation Effects 0.000 claims description 3
- 230000003068 static effect Effects 0.000 claims description 3
- 238000011478 gradient descent method Methods 0.000 claims description 2
- 239000013598 vector Substances 0.000 claims description 2
- 230000006866 deterioration Effects 0.000 claims 1
- 238000004458 analytical method Methods 0.000 abstract description 4
- 230000007547 defect Effects 0.000 abstract description 3
- 230000000694 effects Effects 0.000 description 15
- 238000013508 migration Methods 0.000 description 11
- 230000005012 migration Effects 0.000 description 11
- 238000004088 simulation Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 235000009508 confectionery Nutrition 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000011426 transformation method Methods 0.000 description 1
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention relates to the technical field of seismic exploration, in particular to a combined attenuation random noise processing method based on curvelet transform and total variation. The method comprises the following steps: acquiring single-shot earthquake data; performing curvelet transform on the acquired single-shot data or stacking of single-shot data; performing multi-scale curvelet threshold denoising; denoising by adopting a total-variation denoising technology; outputting earthquake data through graphic display software. According to the method, an optimal threshold is selected according to the distribution rule of random noise in a curvelet domain to maximize the signal to noise ratio of data, and the aim of optimally denoising is fulfilled through curvelet transform; according to the total variation minimization technology; a curvelet coefficient is adjusted through a total variation minimization technology, thereby overcoming the defect of pseudo-curve caused by separate use of curvelet transform, and making displayed stratum data more real and reliable in order to further perform stratum analysis and obtain more accurate analysis results about the oil content and ore content.
Description
Technical Field
The invention relates to the technical field of seismic exploration, in particular to a method for suppressing random noise in seismic data processing.
Background
Curvelet transform is a relatively new multi-scale geometric transformation algorithm. In 1999, candies and Donoho proposed Curvelet transform based on Ridgelet transform, the first generation Curvelet transform; in 2002, Cand's et al proposed a second generation Curvelet transform; in 2005, canddes et al proposed two fast dispersion methods based on the second generation Curvelet transformation theory: 1) a two-dimensional FFT algorithm (USFFT) of non-uniform spatial sampling; 2) the Wrapping algorithm (Wrapping-Based Transform). Curvelet transforms were first applied to the seismic data processing field in 2004 as Curvelet transforms have been widely used in the seismic data processing field due to the superior locality of Curvelet to frequency and direction.
In seismic exploration, a conventional random noise attenuation method has a large damage to effective signals, such as RPF (radial predictive denoising), AMCOD (dip coherent addition), polynomial fitting denoising, median filtering, SVD methods, and RNA denoising methods commonly used in the industry, and although RNA is a better random noise attenuation method in many methods, the method also has a certain damage to effective signals, and in order to better remove random noise, Neelamani et al introduced a multi-scale transformation method, namely, curvelet transformation, to attenuate random noise in 2008, and a better effect is obtained. Although Curvelet transformation has many advantages in denoising, Curvelet transformation also has inherent defects, namely strong energy clusters are easy to generate in a seismic section during denoising, and meanwhile, the unsmooth phenomenon is generated at the edges of the same phase axis.
Disclosure of Invention
The invention aims to provide a combined attenuation random noise processing method based on curvelet transformation and total variation, which can effectively suppress random noise influencing the quality of pre-stack single shot data and post-stack data, furthest reduce the damage to effective signals and is an effective and efficient random noise attenuation technology.
In order to achieve the technical purpose, the invention can adopt the following technical scheme: a combined attenuation random noise processing method based on curvelet transform and total variation sequentially comprises the following steps:
1) acquiring seismic data:
firstly, acquiring single shot data: arranging excitation point positions and a plurality of receiving point positions, then generating seismic waves through explosive explosion excitation embedded in the excitation points, receiving the seismic waves reflected by an underground reflection interface by detectors arranged at all the receiving points, and forming single-shot data by receiving the seismic wave data excited by the same shot by all the detectors;
2) performing curvelet transformation on the obtained single shot data or the superposition of the single shot data
For superimposed data or single shot datafThe curvelet transform is carried out to obtain curvelet transform coefficientsc（j，l，k）Whereinc（j，l，k）Can be expressed as:
in the formula,CT（f）is composed offThe transformation of the curved wave of (a) is carried out,the function of the curvelet wave is expressed,j，l，krespectively show the scale and directionAnd a location parameter.k=(k 1 ，k 2 ) ，k 1 =1,2,...M，k 2 =1,2,…NWhereinMIn order to be the number of seismic traces,Nsampling the number and scale of points in one timej=ceil(log 2 min(M,N)-3)Direction of rotation；
3）Multi-scale curvelet threshold denoising
After the data is subjected to curvelet transformation, dividing the data into j scale layers, wherein each scale layer comprises a plurality of directional data, and the innermost layer is called a low-frequency coefficient layer; the outermost layer is called a high-frequency coefficient layer; the middle scale layer is called a medium-high frequency coefficient layer; in different scale layers, the coefficient distribution of the effective signal and the random noise is different, and in the low-frequency coefficient layer, the coefficients of the effective signal and the random noise are not obviously demarcated, so that all the coefficients are reserved in the scale layer, and each of the high-frequency coefficient layer and the medium-high frequency coefficient layer respectively selects the optimal threshold with the highest signal-to-noise ratio to achieve the optimal denoising effect;
the threshold value calculation method is as follows:
wherein,here, theIs the length of the seismic data;noise in one direction for one dimensionAcoustic variance whileαIs not more than 2.
Then, all curvelet coefficients and threshold values of a certain direction under a certain scaleT α Comparing the sizes of the two components, and reserving the size not less than the threshold valueT α And then reconstructing (also called inverse curvelet transform) the remaining curvelet coefficients as shown in the following equation:
the inverse of the curvelet transform is represented,indicating the passage of a thresholdThe coefficient of the curvelet after the treatment,representing a curvelet function;
transformation ofValue, get differentFor all ofAnalyzing the signal-to-noise ratio to obtain different signal-to-noise ratios，Is composed ofCorresponding signal to noise ratio value such thatHighest value ofThreshold value corresponding to valueThe value is the optimal threshold value under the angle of the scale;
4) total variation denoising
SignalFrom above and below the optimum thresholdConsists of two parts:
for the signalPerforming multi-scale multidirectional threshold processing in a curvelet transform domain, wherein the optimal multi-scale multidirectional threshold isAnd reconstructing to obtain:
thus obtainedResult of noiseThe phenomenon of pseudo Gibbs oscillation is easy to occur, and the phenomenon can be restrained to a certain extent by a total variation minimization technology; here, the data after denoising the curvelet threshold valueIs shown asThen the total variation of the data is as follows:
whereinRepresenting dataThe gradient of (a) of (b) is,for the purpose of the gradient representation,representing dataThe total variation of (a) to (b),as dataIs supported by the support frame, and the support frame,as coordinate vectors of data；
The gradient field of (a) is:
wherein ， As datafIn that(i+1,j)The numerical value of the point, and the like;
the total variation based denoising method can be implemented by minimizing the following function:
the first term is an approximation term, so that the denoised image can still better approximate to the original image and has certain fidelity; the second term is a total variation regularization term, and lambda is a Lagrangian constant, so that a balance is played between an approximation term and the regularization term;
the above objective functionIs thatIs a convex function of which a sufficient requirement for the existence of an extremum isFrom which correspondences can be derivedThe Euler-Largrange equation is:
the equation is non-linear, in whichAssuming that the equation satisfies the Neumann boundary condition for divergence, the data is iterated repeatedly by a gradient descent method until a stable solution is obtained, so as to obtain the denoised data, wherein the iterative formula is as follows:
wherein,to representAs a result of the sub-iteration(s),f n+1 is thatn+1The result of the sub-iteration is initialized，The step size of the iteration is indicated,represents a total variation function inA sub-gradient of (d);
is usually usedTo replaceGet itAnd is a small positive value, which isThe difference is at least two orders of magnitude to solve the irreconcilability of the total variation function at certain points;
using the seismic data and the curvelet coefficient after the best threshold denoising as input, adjusting the curvelet coefficient by utilizing a total variation minimization technology, and setting the maximum iteration times and the initial value，Calculating the sub-gradientStep length is takenCalculating(ii) a GetJudgment ofIf the number of iterations is equal to the maximum number of iterations, ending the iterations, otherwise, continuing the iterations until the iterations are ended;
5) and outputting seismic data: the seismic data is output by graphical display software.
In the step 2), the method for acquiring the overlay data comprises the following steps: and respectively carrying out surface wave and abnormal amplitude attenuation, surface consistency energy compensation, surface consistency deconvolution and surface consistency residual static correction on the obtained single shot data, then sorting to a common central point data domain, carrying out speed explanation, carrying out dynamic correction by using the explained speed, and superposing the dynamically corrected data in the common central point domain to obtain superposed data.
In step 5), the graphic display software may include one of a graphic display tool in Omega processing software, a graphic display tool in Cgg processing software, a graphic display tool in Promax processing software, a graphic display tool in Focus processing software, a graphic display tool in Geoeast processing software, a graphic display tool in Landmark interpretation software, a graphic display tool in geotrace interpretation software, and a wigb program in Matlab.
The invention has the beneficial effects that: the invention denoises through the multi-scale curvelet threshold, for the attenuation of random noise, the optimal threshold is selected according to the distribution rule of the random noise in the curvelet domain, for the problems of different scales and different directions, the curvelet coefficient distribution rule of the random noise is different, the optimal threshold is selected according to the characteristics of each scale and angle data, the signal-to-noise ratio of the data is the highest, and thus the purpose of obtaining the optimal denoising effect by using curvelet transformation is achieved; after denoising of the curvelet threshold, a total variation minimization technology is introduced, the curvelet coefficient is adjusted, and the defects of a 'curvelet-shaped' pseudo curve and the like caused by singly using curvelet transformation are overcome. On the whole, the combined denoising technology is superior to the effect of using Curvelet denoising alone in the prior art. The displayed formation data is more real and reliable, so that formation analysis can be further carried out, and more accurate analysis results of oil and ore are obtained.
Drawings
FIG. 1 is a flow chart of the present invention.
Fig. 2-7 show the results of the present invention in performing random noise attenuation on a single shot seismic data generated by a model forward simulation.
Fig. 2 is simulated single shot seismic data without random noise, fig. 3 is simulated single shot seismic data with random noise added, fig. 4 is a result of suppressing random noise by a conventional single curvelet transform, fig. 5 is a comparison graph of a result of performing multi-scale decomposition on the simulated single shot seismic data with noise in fig. 3 by using curvelet transform, fig. 6 is a result of suppressing random noise by using the present invention, and fig. 7 is error data obtained by subtracting the result of suppressing random noise by using the present invention in fig. 6 from the simulated single shot seismic data without random noise in fig. 2.
Fig. 8-12 show the results of migration stacking seismic data processing in a work area in the northwest exploration area of the oil field in Jiangsu province.
Fig. 8 is original migration stacking data, fig. 9 is a result of a mainstream software Omega in the seismic data processing field after random noise is suppressed by using an RNA random noise attenuation module, fig. 10 is error data obtained by subtracting the original migration stacking data in fig. 8 from the result data after Omega processing in fig. 9, fig. 11 is a processing result obtained by applying the present invention to attenuate random noise, and fig. 12 is error data obtained by subtracting the original migration stacking data in fig. 8 from the result data obtained by applying the present invention in fig. 11.
In each figure, the abscissa is the seismic trace number and the ordinate is the recording time (in ms).
Detailed Description
Example 1
A combined attenuation random noise processing method based on curvelet transform and total variation comprises the following steps:
1) acquiring seismic data:
acquisition of simulation data: single-shot seismic data are simulated through forward modeling of a wave equation model, and technicians provide migration stacking data.
Acquiring actual measurement data: the method comprises the steps of arranging excitation point positions and a plurality of receiving point positions, then exciting to generate seismic waves through explosion of explosives embedded in the excitation points, receiving the seismic waves reflected by an underground reflection interface by receivers arranged on all the receiving points, and forming single-shot data by receiving the seismic wave data excited by the same shot and received by all the receivers.
2) Seismic data of single shot to be acquiredfPerforming curvelet transform
For seismic datafPerforming curvelet transform to obtain curvelet transform coefficientc（j，l，k）Whereinc（j，l，k）Can be expressed as:
in the formula,is composed offThe transformation of the curved wave of (a) is carried out,the function of the curvelet wave is expressed,respectively representing scale, orientation and position parameters.WhereinMIn order to be the number of seismic traces,Nsampling the number and scale of points in one timeDirection of rotation；
The seismic data of the single simulated cannon is 512 tracks and 2048 sampling points, the actual migration stacking data is 200 tracks and 451 sampling points, and the single simulated cannon obtained by utilizing the calculation formula of the scale and the direction can be decomposed intoA unit of dimension andfor each direction, the actual offset stack data can be decomposed intoA unit of dimension andand (4) direction.
3) Multi-scale curvelet optimal threshold denoising
Fixing a certain scale layerjAll the coefficients of the curvelet of other scale layers are set to 0, andjperforming curvelet reconstruction on all curvelet coefficients of the scale layer to obtainjThe scale layer separately reconstructs the data. And performing the same operation on all the scale layers to obtain the data which is independently reconstructed by all the scale layers.
After being subjected to curvelet transformation, the simulated single shot seismic data are divided into 6 scales and 8 directions, and the migration and stacking data are divided into 5 scales and 6 directions; in different scale layers, the distribution of the curvelet coefficients of the effective signals and the random noise is different, and in the low-frequency coefficient layer, the curvelet coefficients of the effective signals and the random noise are not obviously demarcated, so that all the coefficients are reserved in the scale layer, and in each of the high-frequency coefficient layer and the medium-high frequency coefficient layer, an optimal threshold value with the highest noise-to-noise ratio after denoising is selected respectively to achieve the optimal denoising effect.
The method of calculating the threshold is as follows:
wherein,here, theIs the length of the seismic data;the variance of the noise in a certain direction for a certain scale, andαis not more than 2.
Regarding the simulation single shot seismic data, through testing, for the scale layer 1, all the curvelet coefficients are reserved, namely the threshold processing is not carried out on the scale curvelet coefficients; with respect to the scale layer 2, it is,αwhen 2=1.0, the obtained threshold denoising effect is optimal; with respect to the scale layer 3,αwhen 3=1.0, the obtained threshold denoising effect is optimal; with respect to the scale layer 4, it is,αwhen 4=1.2, the obtained threshold denoising effect is optimal; with respect to the scale layer 5, it is,αwhen 5=1.2, the obtained threshold denoising effect is optimal; with respect to the scale layer 6, it is,αwhen 6=0.8, the obtained threshold denoising effect is optimal; with respect to the scale layer 7, it is,αwhen 7=0.5, the obtained threshold denoising effect is optimal; for the scale layer 8, since there is no effective signal and almost random noise, the coefficient of the curved wave of the scale layer is set to 0 in its entirety.
With respect to actual migration stacked seismic data, it was tested that for scale layer 1, all the curvelet coefficients were retained, i.e.(ii) a With respect to the scale layer 2, it is,αwhen 2=1.2, the obtained threshold denoising effect is optimal; with respect to the scale layer 3,αwhen 3=1.2, the obtained threshold denoising effect is optimal; with respect to the scale layer 4, it is,αwhen 4=0.9, the obtained threshold denoising effect is optimal; with respect to the scale layer 5, it is,αwhen 5=0.6745, the obtained threshold denoising effect is optimal.
Then, performing curvelet reconstruction on the curvelet coefficients after performing optimal threshold processing on each scale layer, wherein the simulation single shot data is represented by the following formula:
performing curvelet reconstruction after actual migration superposition data optimal threshold processing and the like, thereby obtaining seismic data after multi-scale curvelet optimal threshold denoising。
4) Total variation denoising
In the process of denoising the optimal threshold of curvelet transformation, the seismic data can be considered as the seismic data after the optimal threshold is givenFrom above and below the optimum thresholdConsists of two parts:
for the signalPerforming multi-scale multidirectional optimal threshold processing in a curvelet transform domain, wherein the multi-scale multidirectional optimal threshold isAnd reconstructing to obtain:
the denoising result obtained therebyThe phenomenon of pseudo Gibbs oscillation is easy to occur, and the linearity can be restrained to a certain degree by a total variation minimization technology; here, the data after denoising processing is performed by making the optimal threshold of curveletIs shown asThen the total variation of the data is as follows:
and carrying out total variation processing by using an iterative formula for solving total variation minimization, wherein the formula is as follows:
wherein,to representAs a result of the sub-iteration(s),is thatThe result of the sub-iteration is initialized，The step size of the iteration is indicated,represents a total variation function inA sub-gradient of (d);
is usually usedTo replace，And is a small positive value, which isAt least two orders of magnitude apart to account for the non-differentiability of the total variation function at certain points, taken here ；
By passing throughAdjusting best threshold denoised seismic dataCoefficient of sum-waveAs input, the total variation minimization technique is used to adjust the curvelet coefficient, and the maximum iteration number is set asInitial value of，Calculating the sub-gradientStep length is takenBy passingThis formula calculates(ii) a Then performing loop iteration to obtainJudgment ofWhether or not to approachIf so, ending iteration, otherwise, continuing iteration until ending iteration, and performing 30 iterations on the simulated single-shot seismic data and the actual migration stack seismic data to obtain a final processing result。
4) And outputting seismic data: the seismic data is output by graphical display software.
The processing results are shown in figures 2-7,
the simulated single shot seismic data is a hyperbolic model, the data size is 512 tracks, 2048 sampling points are sampled for 1ms, the original data without random noise is shown in figure 2, the data after random noise is added is shown in figure 3, the obtained signal-to-noise ratio is 0.28, and the signal-to-noise ratio is expressed by the ratio of the square sum of the signal amplitude to the square sum of the noise amplitude. Fig. 4 shows the result of the test using the conventional curvelet transform method, and the signal-to-noise ratio is 1.2. FIG. 5 is seismic data at six scales after flexural wave transformation of noisy seismic data. Fig. 6 shows the de-noising result of the combined de-noising method based on the curvelet transform and the total variation, the specific threshold value is selected as described above during the curvelet transform, the iteration number of the total variation is 30, the signal-to-noise ratio of the de-noised profile of the invention reaches 4.2, fig. 7 shows the de-noised result of the invention and the error profile of the original noise-free data in fig. 2, and it can be seen from the figure that the combined de-noising technology based on the curvelet transform and the total variation of the invention is far superior to the traditional method for suppressing random noise based on the curvelet transform, and the invention not only suppresses the random noise well, but also protects the effective signal well.
Example 2
As shown in fig. 8-12, the results of the seismic data processing are superimposed for the migration of one work area in the northwest exploration area of the oil field in Jiangsu province.
The difference from embodiment 1 is that in step 2), the curvelet transform is performed by using the superimposed data, and the method for acquiring the superimposed data is as follows: and respectively carrying out surface wave and abnormal amplitude attenuation, surface consistency energy compensation, surface consistency deconvolution and surface consistency residual static correction on the obtained single shot data, then sorting to a common central point data domain, carrying out speed explanation, carrying out dynamic correction by using the explained speed, and superposing the dynamically corrected data in the common central point domain to obtain superposed data. The rest steps and the processing method are the same, and only the specific parameter values are different.
The results after treatment were as follows:
FIG. 8 shows the original offset stack data, with a data size of 200 tracks, 451 sample points, and a sampling interval of 1 ms; the curvelet transform and total variation combined denoising method provided by the invention is compared with an RNA random noise attenuation module of a mainstream software Omega in the field of seismic data processing, wherein fig. 9 is a denoised result of the RNA random noise attenuation module of the Omega, fig. 10 is error data of the denoised data of the RNA in fig. 9 and original offset superposition data in fig. 8, local area data fault development can be seen from the graph, the quality of the data is influenced due to the existence of random noise, the random noise in an original section is well suppressed after the RNA denoising, the signal-to-noise ratio of the section is obviously improved, and meanwhile, a breakpoint of the section is well reserved. However, in a complex-structured area, the effective signal is damaged to a certain extent, the existence of the effective signal can be obviously seen in an error profile, and meanwhile, in a deep weak signal area, the fidelity of the effective signal is relatively poor. FIG. 11 shows the result of denoising by the present invention, and FIG. 12 shows the error data of the denoised data of the present invention in FIG. 11 and the original offset superimposed data in FIG. 8, which can be seen from the denoised section and the error section.
The present invention is not limited to the above-mentioned embodiments, and based on the technical solutions disclosed in the present invention, those skilled in the art can make some substitutions and modifications to some technical features without creative efforts according to the disclosed technical contents, and these substitutions and modifications are all within the protection scope of the present invention. For example, the graphic display software may be one of a graphic display tool in Omega processing software, a graphic display tool in Cgg processing software, a graphic display tool in Promax processing software, a graphic display tool in Focus processing software, a graphic display tool in Geoeast processing software, a graphic display tool in Landmark interpretation software, a graphic display tool in Geoframe interpretation software, and a wigb program in Matlab.
Claims (3)
1. A combined attenuation random noise processing method based on curvelet transform and total variation is characterized by sequentially comprising the following steps:
1) acquiring seismic data:
firstly, acquiring single shot data: arranging excitation point positions and a plurality of receiving point positions, then generating seismic waves through explosive explosion excitation embedded in the excitation points, receiving the seismic waves reflected by an underground reflection interface by detectors arranged at all the receiving points, and forming single-shot data by receiving the seismic wave data excited by the same shot by all the detectors;
2) performing curvelet transformation on the obtained single shot data or the superposition of the single shot data
For superimposed data or single shot dataPerforming curvelet transform to obtain curvelet transform coefficientWhereinCan be expressed as:
in the formula,is composed ofThe transformation of the curved wave of (a) is carried out,the function of the curvelet wave is expressed,respectively representing scale, direction and position parameters;
whereinFor seismic tracesThe number of the first and second groups is,sampling the number and scale of points in one timeDirection of rotation；
3）Multi-scale curvelet threshold denoising
After the data is subjected to curvelet transformation, dividing the data into j scale layers, wherein each scale layer comprises a plurality of directional data, and the innermost layer is called a low-frequency coefficient layer; the outermost layer is called a high-frequency coefficient layer; the middle scale layer is called a medium-high frequency coefficient layer; the low-frequency coefficient layer keeps all coefficients, and each layer in the high-frequency coefficient layer and the middle-high frequency coefficient layer selects the optimal threshold with the highest signal-to-noise ratio;
the threshold value calculation method is as follows:
wherein,here, theIs the length of the seismic data;a noise variance in a direction for a certain scale;
then, all curvelet coefficients and threshold values of a certain direction under a certain scaleT α Comparing the sizes of the two components, and reserving the size not less than the threshold valueT α Then on the remaining coefficient of the curveletReconstruction (also known as inverse curvelet transform) is shown as follows:
the inverse of the curvelet transform is represented,indicating the passage of a thresholdThe coefficient of the curvelet after the treatment,representing a curvelet function;
transformation ofValue, get differentFor all ofAnalyzing the signal-to-noise ratio to obtain different signal-to-noise ratios，Is composed ofCorresponding signal to noise ratio value such thatHighest value ofThreshold value corresponding to valueThe value is the optimal threshold value in the direction of the scale;
4) total variation denoising
SignalFrom above and below the optimum thresholdConsists of two parts:
for the signalPerforming multi-scale multidirectional threshold processing in a curvelet transform domain, wherein the optimal multi-scale multidirectional threshold isAnd reconstructing to obtain:
the denoising result obtained therebyIs susceptible to pseudo-Gibbs oscillation phenomenaThe deterioration minimizing technique can suppress such a phenomenon to some extent; here, the data after denoising the curvelet threshold valueIs shown asThen the total variation of the data is as follows:
whereinRepresenting dataThe gradient of (a) of (b) is,for the purpose of the gradient representation,representing dataThe total variation of (a) to (b),as dataIs supported by the support frame, and the support frame,is a coordinate vector of the data;
the gradient field of (a) is:
wherein As dataIn thatThe numerical value of the point, and the like;
the total variation based denoising method can be implemented by minimizing the following function:
the first term is an approximation term, so that the denoised image can still better approximate to the original image and has certain fidelity; the second term is the total variation regularization term, λ is the Lagrangian constant;
the above objective functionIs thatIs a convex function of which a sufficient requirement for the existence of an extremum isFrom this, it can be derived that the corresponding Euler-largrage equation is:
the equation is non-linear, in whichAssuming that the equation satisfies the Neumann boundary condition for divergence, the data is iterated repeatedly by a gradient descent method until a stable solution is obtained, so as to obtain the denoised data, wherein the iterative formula is as follows:
wherein,to representAs a result of the sub-iteration(s),is thatThe result of the sub-iteration is initialized，The step size of the iteration is indicated,represents a total variation function inA sub-gradient of (d);
by usingTo replaceGet itAnd is a small positive value, which isBy at least two orders of magnitude;
using the seismic data and the curvelet coefficient after the best threshold denoising as input, adjusting the curvelet coefficient by utilizing a total variation minimization technology, and setting the maximum iteration times and the initial value，Calculating the sub-gradientStep length is takenCalculating(ii) a GetJudgment ofIf the number of iterations is equal to the maximum number of iterations, ending the iterations, otherwise, continuing the iterations until the iterations are ended;
5) and outputting seismic data: the seismic data is output by graphical display software.
2. The method for processing the joint attenuation random noise based on the curvelet transform and the total variation as claimed in claim 1, wherein the method for obtaining the superposition data in the step 2) comprises: and respectively carrying out surface wave and abnormal amplitude attenuation, surface consistency energy compensation, surface consistency deconvolution and surface consistency residual static correction on the obtained single shot data, then sorting to a common central point data domain, carrying out speed explanation, carrying out dynamic correction by using the explained speed, and superposing the dynamically corrected data in the common central point domain to obtain superposed data.
3. The method according to claim 1 or 2, wherein in step 5), the graphic display software comprises one of a graphic display tool in Omega processing software, a graphic display tool in Cgg processing software, a graphic display tool in Promax processing software, a graphic display tool in Focus processing software, a graphic display tool in Geoeast processing software, a graphic display tool in Landmark interpretation software, a graphic display tool in geofield interpretation software, and a wigb program in Matlab.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410107764.0A CN104062683A (en) | 2014-03-21 | 2014-03-21 | Combined attenuation random noise processing method based on curvelet transform and total variation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410107764.0A CN104062683A (en) | 2014-03-21 | 2014-03-21 | Combined attenuation random noise processing method based on curvelet transform and total variation |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN104062683A true CN104062683A (en) | 2014-09-24 |
Family
ID=51550481
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410107764.0A Pending CN104062683A (en) | 2014-03-21 | 2014-03-21 | Combined attenuation random noise processing method based on curvelet transform and total variation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104062683A (en) |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104932017A (en) * | 2015-05-26 | 2015-09-23 | 中国石油天然气集团公司 | Anomalous amplitude suppression method and apparatus |
| CN106597539A (en) * | 2016-12-28 | 2017-04-26 | 中国石油化工股份有限公司 | Curvelet domain Radon transform noise suppression method for loess tableland region |
| CN106680876A (en) * | 2017-01-22 | 2017-05-17 | 中国石油大学(华东) | Joint denoising method for seismic data |
| CN107589450A (en) * | 2017-09-01 | 2018-01-16 | 中国科学院地质与地球物理研究所 | Seismic data noise attenuation method and apparatus based on warp wavelet and cluster |
| CN109917459A (en) * | 2017-12-12 | 2019-06-21 | 中国石油天然气股份有限公司 | Method, device and system for suppressing earthquake noise |
| CN110007346A (en) * | 2019-03-20 | 2019-07-12 | 中国石油化工股份有限公司 | Adaptive intermediate value very noisy based on wavelet decomposition interferes seismic data noise attenuation method |
| CN110244361A (en) * | 2019-06-21 | 2019-09-17 | 同济大学 | Seismic data separation method and system based on simultaneous processing of interpolation and anti-aliasing |
| CN110244360A (en) * | 2019-06-21 | 2019-09-17 | 同济大学 | Seismic data separation method and system based on anti-aliasing in effective frequency wavenumber domain |
| CN111158051A (en) * | 2020-01-07 | 2020-05-15 | 自然资源部第一海洋研究所 | Joint constraint random noise suppression method based on sparse regularization |
| CN111273351A (en) * | 2019-11-21 | 2020-06-12 | 西安工业大学 | Structural guide direction generalized total variation regularization method for seismic data denoising |
| CN112147692A (en) * | 2019-06-28 | 2020-12-29 | 中国石油化工股份有限公司 | Secondary seismic source attenuation method for marine seismic data |
| CN113740908A (en) * | 2020-05-29 | 2021-12-03 | 中国石油化工股份有限公司 | Two-dimensional variation analysis method for seismic slice, electronic device, and medium |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009051900A2 (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 |
| US20110069582A1 (en) * | 2009-09-22 | 2011-03-24 | David Nichols | Slant-stack for time gather to angle gather conversion |
| CN102141634A (en) * | 2010-12-17 | 2011-08-03 | 西南交通大学 | Method for suppressing interference of neutral line of prestack seismic signal based on curvelet transform |
-
2014
- 2014-03-21 CN CN201410107764.0A patent/CN104062683A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009051900A2 (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 |
| US20110069582A1 (en) * | 2009-09-22 | 2011-03-24 | David Nichols | Slant-stack for time gather to angle gather conversion |
| CN102141634A (en) * | 2010-12-17 | 2011-08-03 | 西南交通大学 | Method for suppressing interference of neutral line of prestack seismic signal based on curvelet transform |
Non-Patent Citations (2)
| Title |
|---|
| 俞华 等: "曲波变换及全变差最小化技术联合去噪", 《石油物探》 * |
| 薛永安 等: "改进的曲波变换及全变差联合去噪技术", 《物探与化探》 * |
Cited By (20)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104932017A (en) * | 2015-05-26 | 2015-09-23 | 中国石油天然气集团公司 | Anomalous amplitude suppression method and apparatus |
| CN106597539A (en) * | 2016-12-28 | 2017-04-26 | 中国石油化工股份有限公司 | Curvelet domain Radon transform noise suppression method for loess tableland region |
| CN106597539B (en) * | 2016-12-28 | 2019-07-12 | 中国石油化工股份有限公司 | For the bent wave zone Radon converter noise drawing method of Huangtuyuan area |
| CN106680876A (en) * | 2017-01-22 | 2017-05-17 | 中国石油大学(华东) | Joint denoising method for seismic data |
| CN106680876B (en) * | 2017-01-22 | 2019-04-12 | 中国石油大学(华东) | A kind of seismic data joint denoising method |
| CN107589450A (en) * | 2017-09-01 | 2018-01-16 | 中国科学院地质与地球物理研究所 | Seismic data noise attenuation method and apparatus based on warp wavelet and cluster |
| CN107589450B (en) * | 2017-09-01 | 2019-01-04 | 中国科学院地质与地球物理研究所 | Seismic data noise attenuation method and apparatus based on warp wavelet and cluster |
| CN109917459A (en) * | 2017-12-12 | 2019-06-21 | 中国石油天然气股份有限公司 | Method, device and system for suppressing earthquake noise |
| CN110007346A (en) * | 2019-03-20 | 2019-07-12 | 中国石油化工股份有限公司 | Adaptive intermediate value very noisy based on wavelet decomposition interferes seismic data noise attenuation method |
| CN110244360A (en) * | 2019-06-21 | 2019-09-17 | 同济大学 | Seismic data separation method and system based on anti-aliasing in effective frequency wavenumber domain |
| CN110244361A (en) * | 2019-06-21 | 2019-09-17 | 同济大学 | Seismic data separation method and system based on simultaneous processing of interpolation and anti-aliasing |
| CN110244361B (en) * | 2019-06-21 | 2020-07-07 | 同济大学 | Seismic data separation method and system based on simultaneous processing of interpolation and anti-aliasing |
| CN110244360B (en) * | 2019-06-21 | 2020-07-07 | 同济大学 | Seismic data separation method and system based on effective frequency wave number domain de-aliasing |
| CN112147692A (en) * | 2019-06-28 | 2020-12-29 | 中国石油化工股份有限公司 | Secondary seismic source attenuation method for marine seismic data |
| CN112147692B (en) * | 2019-06-28 | 2022-08-05 | 中国石油化工股份有限公司 | Secondary seismic source attenuation method for marine seismic data |
| CN111273351A (en) * | 2019-11-21 | 2020-06-12 | 西安工业大学 | Structural guide direction generalized total variation regularization method for seismic data denoising |
| CN111273351B (en) * | 2019-11-21 | 2022-04-08 | 西安工业大学 | Structural guide direction generalized total variation regularization method for seismic data denoising |
| CN111158051A (en) * | 2020-01-07 | 2020-05-15 | 自然资源部第一海洋研究所 | Joint constraint random noise suppression method based on sparse regularization |
| CN113740908A (en) * | 2020-05-29 | 2021-12-03 | 中国石油化工股份有限公司 | Two-dimensional variation analysis method for seismic slice, electronic device, and medium |
| CN113740908B (en) * | 2020-05-29 | 2024-05-07 | 中国石油化工股份有限公司 | Two-dimensional variogram analysis method, electronic equipment and medium for seismic slice |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104062683A (en) | Combined attenuation random noise processing method based on curvelet transform and total variation | |
| Kaur et al. | Seismic ground‐roll noise attenuation using deep learning | |
| US20200225377A1 (en) | Methods for simultaneous source separation | |
| Innocent Oboué et al. | Robust damped rank-reduction method for simultaneous denoising and reconstruction of 5D seismic data | |
| CN110879416B (en) | Method for suppressing multiple between layers | |
| AU2015383134B2 (en) | Multistage full wavefield inversion process that generates a multiple free data set | |
| CN102831588B (en) | De-noising processing method for three-dimensional seismic images | |
| CN110954945B (en) | A full waveform inversion method based on dynamic stochastic source coding | |
| CN114488302B (en) | In-situ anisotropic ground stress field prediction method and system | |
| US20160320507A1 (en) | Time lapse seismic data processing | |
| US10330809B2 (en) | Device and method for optimization of 4D and 3D seismic data | |
| Liu et al. | Source-independent time-lapse full-waveform inversion for anisotropic media | |
| US10634804B2 (en) | System and method for dip-guided seismic image stacking | |
| CN116719086B (en) | Sparse seabed four-component data high-resolution imaging method based on point spread function | |
| Zhu et al. | Stretch-free migration imaging condition | |
| US20160202376A1 (en) | Stack Ghost Suppression | |
| Tian et al. | Joint sparsity recovery for noise attenuation | |
| CN109655916B (en) | Method and system for separating effective waves and multiples in seismic data | |
| US20150142831A1 (en) | Computerized method and a computer program rpoduct for determining a resulting data set representative of a geological region of interest | |
| Sanchis et al. | Multicomponent streamer noise characteristics and denoising | |
| Wang et al. | Passive seismic data primary estimation and noise removal via focal‐denoising closed‐loop surface‐related multiple elimination based on 3D L1‐norm sparse inversion | |
| Robertsson et al. | Efficient snapshot-free reverse time migration and computation of multiparameter gradients in full-waveform inversion | |
| Wang et al. | Practical strategies for data-driven land multiple attenuation | |
| CN115951402B (en) | Offset imaging method, system and medium based on vector reflectivity forward operator | |
| Sun et al. | Least-squares reverse time migration in pseudo-depth domain and its application |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20140924 |
|
| RJ01 | Rejection of invention patent application after publication |