Edge-preserving diffusion filtering method based on fault confidence parameter controlTechnical Field
The invention relates to a fault and fault fracture zone identification and description method, belongs to the field of seismic data interpretation, and particularly relates to an edge-preserving diffusion filtering method based on fault confidence parameter control.
Background
With the continuous improvement of the exploration and development degree of oil fields, the position of a fault block type oil and gas reservoir is more and more important, the geological features of the fault block type oil reservoir are more complex in general, the acquired data also have a certain degree of noise, and particularly, the near fault seismic reflection data have stronger noise compared with other area data. Most fault detection methods currently used mainly use local rapid changes at fault positions to identify and describe, so that the methods are sensitive to noise. This requires that we have noise suppressed the seismic data in advance while preserving the edge information.
The applicant previously applied for CN108415077a new low-order fault identification method for edge detection (publication No. CN108415077 a). The new edge detection low-order fault identification method comprises the following steps: step 1, suppressing random noise according to a denoising technology combining three-dimensional multi-stage blind source separation and structure-preserving filtering in a time space domain of a post-stack seismic data body, and improving the signal to noise ratio of the seismic data; step 2, performing dip angle search to estimate a stratum dip angle value in a time domain by using a multi-channel coherent algorithm, and estimating the stratum dip angle by using time delay represented by the stratum dip angle on the seismic record in a frequency domain and performing Fourier transform on the time delay characteristic of the stratum dip angle; step 3, calculating and analyzing the intra-window structure guide Canni attribute on the basis of local layer leveling along the stratum dip angle based on time-frequency domain mixed dip angle scanning; and 4, carrying out seismic identification of the low-order faults in the forms of a continuous well section, a surface slice and the like according to the structure-oriented Canni attribute data body. According to the method, random noise is suppressed by using a denoising technology combining three-dimensional multi-stage blind source separation and structure-preserving filtering in a time space domain according to the post-stack seismic data body, and the signal to noise ratio of the seismic data is improved; performing dip angle search by adopting a multi-channel coherent algorithm to estimate a stratum dip angle value; and calculating and analyzing the structure guide Canni attribute in the window on the basis of local layer leveling along the stratum dip angle to form an effective low-order fault identification method.
However, the invention only solves the problem of low-order fault identification, has poor applicability, and can cause the loss of information data of discontinuous boundaries such as original faults in the data processing process.
Disclosure of Invention
The invention aims to solve the problems in the prior art and provides a novel edge-preserving diffusion filtering method based on fault confidence parameter control, which is suitable for denoising seismic data and protecting fault and other geological boundary information.
The method is mainly used for improving the signal-to-noise ratio of seismic data and protecting the stratum boundary and breakpoint position information of geologic bodies such as faults, the method adopts Gaussian kernel functions to preprocess the seismic data, obtains characteristic values and corresponding characteristic vectors through reconstruction structure tensor matrix calculation, redesigns the characteristic values of diffusion tensors by combining fault confidence coefficient parameters, solves a three-dimensional anisotropic diffusion equation through an iteration method, and finally obtains filtered data.
The technical proposal is as follows:
an edge preserving diffusion filtering method based on fault confidence parameter control, comprising:
step 1: firstly, preprocessing seismic data by using a Gaussian kernel function with a specific noise scale;
step 2: reconstructing a structure tensor matrix;
step 3: obtaining eigenvalues and corresponding eigenvectors of the structure tensor matrix;
step 4: obtaining fault confidence parameters;
step 5: redesigning a diffusion tensor eigenvalue according to the fault confidence coefficient parameters obtained in the step 4;
step 6: after the diffusion tensor is calculated, the three-dimensional anisotropic diffusion equation is solved through an iteration method, and finally the seismic data after diffusion filtering is obtained.
The above scheme further includes:
the step 2 of reconstructing the structure tensor matrix is as follows:
wherein S isρ Tensor of structure, Uσ ＝Kσ *U,Kσ Is a Gaussian kernel function, and the expression is shown in a formula (2); * Is a convolution operator, and sigma is a noise scale;representing the transpose of the matrix multiplied by the matrix; k (K)ρ As Gaussian kernel function, rho integration scale should be generally larger than noise scale sigma, [ V ] Uσ For gradient information of data, alpha is a stability factor, the size is between 0 and 1, and is used for controlling the contribution rate of the second derivative,>is second derivative information;
the step 3 is to calculate the eigenvalue of the structure tensor matrix and the corresponding eigenvector to the structure tensor Sρ Performing characteristic decomposition to obtain;
for the structure tensor Sρ Beta, as1 、β2 、β3 For its three characteristic values and beta1 ≥β2 ≥β3 The corresponding feature vector is ζ1 、ζ2 And zeta3 ，ζ1 Indicating the direction of maximum gradient change ζ3 Representing the direction in which the gradient change is minimal, epsilon represents the transverse discontinuity factor.
According to the magnitudes of the three eigenvalues and the corresponding different eigenvectors, the structure of the seismic image is expressed into four types:
beta (one)1 ≈0、β2 ≈0、β3 The three feature vectors have little or no change in direction and are expressed as smooth areas, which indicates that the geologic body is isotropic;
(II) beta1 >0、β2 ≈0、β3 Approximately 0, indicated in ζ1 The change rate of the direction is larger, and the change rates of the other two characteristic directions are very small and almost equal, which represents a plane linear geologic body;
(III) beta1 >0、β2 >0、β3 About 0, indicating that one of the direction change rates is approximately zero, indicating discontinuous information appearing at a discontinuous boundary such as a fault, and the corresponding model is a linear structure model;
(IV) beta1 >0、β2 >0、β3 >And 0, the change rates of the three eigenvector directions are different, and the change rates of all directions are not zero, which indicates that the area has no obvious seismic structure information and the seismic reflection is disordered.
The step 4 of obtaining fault confidence coefficient parameters comprises linear confidence coefficient Mline Sum-of-planar confidence metric Mplane (formula (4)) for detecting the relationship between the image structure of the seismic data and the structure tensor eigenvalue:
Mplane and Mline In the range of 0 to 1, the two represent the similarity degree of the neighborhood homoplanar structure and the linear structure of a certain point respectively, and M is calculatedline And Mplane Combining together to obtain confidence metrics M for detecting discontinuity structure boundary information including faultsfault ：
Step 5 is to redesign the diffusion tensor eigenvalue according to the fault confidence coefficient parameters obtained in step 4, and in order to ensure that the diffusion direction is along the construction direction, the eigenvector of the diffusion tensor D should be matched with the structure tensor Sρ Is obtained by:
wherein gamma is1 、γ2 And gamma3 Three eigenvalues of the diffusion tensor D correspond to the eigenvector ζ1 、ζ2 And zeta3 。
In the three-dimensional anisotropic diffusion filtering, the magnitude of the diffusion filtering intensity depends on the eigenvalue gamma of the diffusion tensor D1 、γ2 And gamma3 The closer the eigenvalue is to 0, the smaller the intensity of the diffusion filtering; the characteristic value approaches 1, and the diffusion filtering strength is larger.
Will Mline 、Mplane Fault confidence measure Mfault The following form is obtained by introducing the characteristic value of the diffusion tensor D into the design:
where a is a positive number close to 0,s and τ represent threshold and slope, respectively, as Sigmoid function, function hτ (x) The intensity of diffusion filtering between two eigenvectors can be effectively ensured.
In summary, the design method of the characteristic value can enhance the diffusion filtering effect and simultaneously maintain boundary information of special geological structures such as faults and the like.
After the diffusion tensor D is calculated in the step 6, solving a three-dimensional anisotropic diffusion equation by an iteration method to finally obtain seismic data after diffusion filtering:
Uk and Uk+1 Respectively representing the data obtained by the k and k+1 iterative computation, wherein Deltat is the iterative step length, and the data is stable in practical applicationAnd the value of the property should be small.
Compared with the prior art, the fault confidence coefficient parameter is introduced into the anisotropic filtering, and the diffusion coefficient in the direction of maximum gradient change, namely in the direction of maximum feature vector, approaches to 0, so that the strength of the diffusion filtering is smaller, and boundary information of discontinuous structures such as faults can be protected; along the directions of the other two feature vectors, the sizes of the feature values can be adjusted according to three confidence parameters, namely planar confidence, linear confidence and fault confidence. If the fault confidence degree tends to 0, the stratum structure is characterized by a planar structure, and the diffusion model is along gamma2 And gamma3 The formed plane filtering strength is high; if the fault confidence tends to be 1, gamma2 And gamma3 The value of (2) is an integer close to 0, and the diffusion filtering strength is weak along the directions of the two eigenvectors; the information of discontinuous boundaries such as faults can be effectively protected.
Drawings
FIG. 1 is a flow chart of an edge preserving diffusion filtering method based on fault confidence parameter control.
Fig. 2 is an effect diagram of theoretical verification on a three-dimensional Qdome noise adding model (signal to noise ratio is 5, 3 and 1 respectively) by adopting the method, wherein: (a) an original seismic profile (snr=1), (b) a post-filter seismic profile (snr=1), (c) a pre-filter similarity attribute (snr=1), (d) a post-filter similarity attribute (snr=1), (e) an original seismic profile (snr=3), (f) a post-filter seismic profile (snr=3), (g) a pre-filter similarity attribute (snr=3), (h) a post-filter similarity attribute (snr=3), (i) an original seismic profile (snr=5), (j) a post-filter seismic profile (snr=5), (k) a pre-filter similarity attribute (snr=5), (l) a post-filter similarity attribute (snr=5).
FIG. 3 is a view of an actual seismic section before and after filtering using the present method, wherein: (a) A pre-filter inline seismic profile, (b) a post-filter inline seismic profile.
FIG. 4 is a cross-sectional view of the similarity attributes of actual seismic data before and after filtering using the present method, wherein: (a) A pre-filter similarity attribute, (b) a post-filter similarity attribute.
FIG. 5 is a coherent slice of real seismic data before and after filtering using the present method, wherein: (a) Pre-filter coherence time slicing, (b) post-filter coherence time slicing.
Detailed Description
The present invention will be described in further detail with reference to the following examples in order to make the objects, technical solutions and advantages of the present invention more apparent. The specific embodiments described herein are for purposes of illustration only and are not intended to limit the invention.
Example 1
A three-dimensional anisotropic diffusion filtering method based on confidence parameter control comprises the following steps.
Step 1: seismic data is first preprocessed using a gaussian kernel function of a particular noise scale.
Step 2: reconstructing the structure tensor matrix.
In Uσ ＝Kσ *U,Kσ Is a Gaussian kernel function, and the expression is shown in a formula (2); * Is a convolution operator, and sigma is a noise scale;representing the transpose of the matrix multiplied by the matrix; k (K)ρ As Gaussian kernel function, rho integration scale should be generally larger than noise scale sigma, [ V ] Uσ Alpha is a stable coefficient with the size between 0 and 1 for controlling the contribution rate of the second derivative,is second derivative information.
Step 3: and obtaining the eigenvalue and the corresponding eigenvector of the structure tensor matrix.
For the structure tensor Sρ Performing characteristic decomposition to obtain;
for the structure tensor Sρ Beta, as1 、β2 、β3 For its three characteristic values and beta1 ≥β2 ≥β3 The corresponding feature vector is ζ1 、ζ2 And zeta3 ，ζ1 Indicating the direction of maximum gradient change ζ3 Indicating the direction in which the gradient change is minimal. Epsilon represents the transverse discontinuity factor.
The seismic image structure may be expressed in several types as described in table 1 according to the magnitude of the eigenvalues and their corresponding different eigenvectors.
TABLE 1 relationship of Structure tensor eigenvalues to seismic image Structure
Step 4: and (5) obtaining fault confidence parameters.
Two confidence parameters-Linear confidence metrics Mline Sum-of-planar confidence metric Mplane (formula (4)) for detecting the relationship between the image structure of the seismic data and the structure tensor eigenvalue:
Mplane and Mline In the value range of 0 to 1, the two represent the similarity degree of the neighborhood homoplanar structure and the linear structure of a certain point respectively, Mplane And Mline The magnitude relation between the same characteristic values is shown in table 2. Will Mline And Mplane Confidence measure M for detecting discontinuous structure boundary information such as faultfault ：
On the seismic section, if the continuity of the seismic event is good, the linear confidence measure Mline Tend to 0, planar confidence measure Mplane Trend 1, fault confidence measure Cfault Trend towards 0; similarly, in a region where the phase axis is discontinuous, such as a fault, the linear confidence Mline Tend to 1, planar confidence Mplane Trend to 0, fault confidence parameter Mfault Tend to be 1.
Step 5: redesigning the diffusion tensor eigenvalue according to the fault confidence coefficient parameters obtained in the step 4
In order to ensure that the diffusion direction is along the construction direction, the eigenvectors of the diffusion tensor D should be equal to the structure tensor Sρ Can be derived from the agreement of (1):
wherein gamma is1 、γ2 And gamma3 Three eigenvalues of the diffusion tensor D correspond to the eigenvector ζ1 、ζ2 And zeta3 。
In the three-dimensional anisotropic diffusion filtering, the magnitude of the diffusion filtering intensity depends on the eigenvalue gamma of the diffusion tensor D1 、γ2 And gamma3 . The closer the eigenvalue is to 0, the smaller the intensity of the diffusion filter; the characteristic value approaches 1, and the diffusion filtering strength is larger.
To overcome the defect that the conventional diffusion filtering method can not keep the data edge information, M is addedline 、Mplane Fault confidence measure Mfault The characteristic value of the diffusion tensor D is introduced into the design, and the following form can be obtained：
Where a is a positive number close to 0,s and τ represent threshold and slope, respectively, as Sigmoid function, function hτ (x) The intensity of diffusion filtering between two eigenvectors can be effectively ensured.
The principle of the characteristic value design scheme of the diffusion tensor shown in the formula (7) is as follows: along the direction of greatest gradient change, i.e. along the eigenvector ζ1 Direction, gamma1 The value of (a) is a, and the diffusion coefficient in the direction approaches to 0, so that the strength of diffusion filtering is smaller, and boundary information of discontinuous structures such as faults can be protected; along the eigenvector ζ2 And zeta3 The magnitude of the characteristic value can be based on Mplane 、Mline M is as followsfault These three confidence parameters are adjusted. If Mfault 0, the stratum structure is characterized by a planar structure, gamma2 And gamma3 Has a value close to 1, and the diffusion model is along ζ2 And zeta3 The formed plane filtering strength is high; if Mfault →1，γ1 And gamma2 The value of (a) is an integer a close to 0, and the diffusion filtering strength is weak along the directions of the two eigenvectors; the information of discontinuous boundaries such as faults can be effectively protected.
In summary, the design method of the characteristic value can enhance the diffusion filtering effect and simultaneously maintain boundary information of special geological structures such as faults and the like.
Example 2
Based on the above embodiment 1, step 6: after the diffusion tensor D is calculated, the three-dimensional anisotropic diffusion equation can be solved through an iteration method, and finally the seismic data after diffusion filtering is obtained:
Uk and Uk+1 The data obtained by the k and k+1 iterative calculation are respectively represented, Δt is the iterative step length, and in practical application, the value of Δt should be smaller for stability.
Example 3
With reference to the accompanying drawings, the edge preserving diffusion filtering method based on fault confidence parameter control mainly comprises the following steps:
s1: seismic data is first preprocessed using a gaussian kernel function of a particular noise scale.
S2: reconstructing the structure tensor matrix.
S3: and obtaining the eigenvalue and the corresponding eigenvector of the structure tensor matrix.
The seismic image structure can be expressed into four types according to the magnitude of the eigenvalues and the corresponding different eigenvectors.
S4: and (5) obtaining fault confidence parameters.
Mplane And Mline And in the range of 0 to 1, the two values respectively represent the similarity degree of the neighborhood homoplanar structure and the linear structure of a certain point.
On the seismic section, if the continuity of the seismic event is good, the linear confidence measure Mline Tend to 0, planar confidence measure Mplane Trend 1, fault confidence measure Cfault Trend towards 0; similarly, in a region where the phase axis is discontinuous, such as a fault, the linear confidence Mline Tend to 1, planar confidence Mplane Trend to 0, fault confidence parameter Mfault Tend to be 1.
S5: and (4) redesigning the diffusion tensor eigenvalue according to the fault confidence coefficient parameters obtained in the step (4).
Where a is a positive number close to 0,sigmoid function proposed for Terebes et al (2005), s and τ represent the threshold and slope, respectively, function hτ (x) The intensity of diffusion filtering between two eigenvectors can be effectively ensured.
The principle of the characteristic value design scheme of the diffusion tensor shown in the formula (3) is as follows: along the direction of greatest gradient change, i.e. along the eigenvector ζ1 Direction, gamma1 The value of (a) is a, and the diffusion coefficient in the direction approaches to 0, so that the strength of diffusion filtering is smaller, and boundary information of discontinuous structures such as faults can be protected; along the eigenvector ζ2 And zeta3 The magnitude of the characteristic value can be based on Mplane 、Mline M is as followsfault These three confidence parameters are adjusted. If Mfault 0, the stratum structure is characterized by a planar structure, gamma2 And gamma3 Has a value close to 1, and the diffusion model is along ζ2 And zeta3 The formed plane filtering strength is high; if Mfault →1，γ1 And gamma2 The value of (a) is an integer a close to 0, and the diffusion filtering strength is weak along the directions of the two eigenvectors; the information of discontinuous boundaries such as faults can be effectively protected.
S6: after the diffusion tensor is calculated, the three-dimensional anisotropic diffusion equation can be solved through an iteration method, and finally the seismic data after diffusion filtering is obtained.
The invention obtains the final diffusion-filtered seismic data. According to the invention, the signal-to-noise ratio of the seismic data can be effectively improved, the form and the spread of faults on the filtered seismic section can be more easily and accurately described, and small faults which are difficult to identify in the original data can be accurately and clearly identified after the filtering.
Example 4
The following are specific examples of applications of the present invention:
the method is applied to the three-dimensional post-stack seismic data, and finally the inline seismic profile shown in fig. 2, the similarity attribute profile shown in fig. 3 and the coherence slice shown in fig. 4 are obtained.
As can be seen from FIG. 2, the filtering method of the invention can better control the diffusion intensity, can better protect the fault edge information while denoising, has clearer fault boundary after filtering, has obvious improvement of the similarity attribute after filtering, and is easier to identify the fault edge on the attribute section.
As can be seen from FIG. 3, on the actual seismic section after the filtering by the method, the same phase axis is more continuous, and faults are easier to identify. As can be seen from fig. 4 and fig. 5, the form and the spread of the fault are more easily and accurately described on the filtered attribute map, and some small faults which are not easy to identify in the original data are accurately and clearly identified on the form spread after filtering, which also proves the effectiveness of the filtering method in maintaining the fault edge information.
The related technical content which is not described in the mode can be realized by adopting or referring to the prior art.
It is noted that one skilled in the art, with the benefit of this disclosure, may make modifications, such as equivalents, or obvious variations thereof. The above-mentioned variants are all within the scope of the invention.