CN103489159B - Based on the 3D seismic data image denoising method of three limit structure directing filtering - Google Patents

Abstract

The invention discloses a kind of 3D seismic data image denoising method based on three limit structure directing filtering, based on the framework of traditional structure Steerable filter, improve the structure of structure tensor, and in conjunction with protecting the good three limit filtering of limit effect, significantly promote the performance of original Gaussian kernel anisotropic diffusion filtering, make denoising effect good, structural information loss is few, in the noise reduction of the guarantor limit of seismic image, bring into play very big effect, can make seismic image in denoising, retain more structural informations, be convenient to explanation and the subsequent treatment of seismic image.

Description

Three-dimensional seismic data image noise reduction method based on trilateral structure guided filtering

Technical Field

The invention relates to the field of seismic image processing, in particular to a three-dimensional seismic data image noise reduction method based on trilateral structure guided filtering.

Background

The seismic data contain a lot of geological structure information, such as geological units of rock strata, lenticles, salt domes, pincushion bodies and the like, and structural units of anticline, syncline, fault, lithologic invader, unconformity and the like, which are important indicators for searching oil and gas reservoirs and are target areas searched by oil and gas seismic exploration. However, due to the complexity of the ground seismic wave signal acquisition environment, the seismic data have strong noise, and in addition, the complex geological structures reflect the seismic wave signals weakly, so that a target area in an imaging graph of the seismic data is fuzzy and unclear, and even submerged in the noise, the target area cannot be identified.

The denoising method of the seismic data is divided into two categories of linear filtering and nonlinear filtering. The linear filtering includes Gaussian filtering, F-x prediction filtering, F-K domain filtering, multi-channel coherent filtering and the like, and the nonlinear filtering includes simple median filtering, classical adaptive filtering, filtering based on mathematical morphology, filtering based on wavelet theory, filtering based on partial differential equation and the like. When the linear filtering method is used for removing background noise interference, the fuzzy of discontinuous structural characteristics such as fracture, geologic body edge and the like can be caused, the edge protection capability is not provided, and the edge protection effect of the nonlinear filtering method is good but is respectively long. The most typical nonlinear filtering method is filtering based on partial differential equation-structure oriented filtering, which has excellent edge-preserving effect because of fine filtering depending on the structural information of the original image, and is widely applied to edge-preserving and noise-reducing of seismic data.

The structure-oriented filtering method is characterized in that a structure tensor is obtained by carrying out Gaussian low-pass filtering on a gradient tensor of a seismic image, a diffusion tensor is built by depending on the structure tensor, an iterative equation of the structure-oriented filtering is built according to the diffusion tensor for filtering, the filtering method is also named, anisotropic diffusion filtering reflects that the filtering is a self-diffusion process of the image, and self-diffusion filtering is carried out under the structural constraint of the image.

A filtering method for researching anisotropic diffusion by geological experts at home and abroad provides a series of methods:

(1) foreign Fhemers and Hocker firstly introduce a diffusion filtering technology into rapid processing and interpretation of seismic data in 2003, and provide an anisotropic diffusion filtering technology of structural constraint and edge protection, which is used for enhancing the structural characteristics of seismic data to facilitate interpretation; laviale et al propose a diffusion filtering method for seismic fault retention in 2007, so that the fault in the processed seismic data is more prominent; kadlec et al propose an anisotropic diffusion method for maintaining sequence characteristics in 2008 for enhancing continuity of river channel sand body characteristics in seismic data and improving subsequent automatic river channel identification and segmentation capabilities.

(2) The application of the consistency enhancement diffusion filtering technology in nonlinear anisotropic diffusion in two-dimensional seismic section edge-preserving filtering is explored respectively in 2004 by Sunshinping et al, Wangzhou pine and Yangchun in 2006 in China. The nonlinear anisotropic diffusion filtering method for the three-dimensional seismic data is researched in 2010 by Zhang et al, and the discontinuous structure confidence measurement factor is introduced to regulate and control diffusion filtering, so that a better processing result is obtained.

Although the traditional Gaussian kernel structure-oriented filtering method can perform structure protection on the noise reduction of the three-dimensional seismic image, the construction of the structure tensor depends on Gaussian filtering, and the Gaussian filtering is a smoothing method without edge preserving capability, so that the generated structure tensor lacks the detailed structure information of a plurality of images, the filtered image structure information is still more, and the method is very unfavorable for the weak structure characteristics in the seismic image.

Disclosure of Invention

The technical problem to be solved by the invention is to provide a three-dimensional seismic data image denoising method based on trilateral structure guided filtering, so that more structural information is kept while denoising is carried out on a three-dimensional seismic image, convenience is brought to searching of a target area of the three-dimensional seismic image, and interpretation and subsequent processing of the three-dimensional seismic image are facilitated.

In order to solve the technical problems, the technical scheme adopted by the invention is as follows: a three-dimensional seismic data image noise reduction method based on trilateral structure guided filtering is disclosed, and the method comprises the following steps:

1) setting trilateral structure guiding filtering parameters of the three-dimensional seismic image data: filtering to compute structure tensorParameter sigma₁The time step of diffusion filtering, namely, the time step of diffusion filtering, and the iteration times of trilateral structure-oriented filtering, namely, iterTimeTimes;

2) importing three-dimensional seismic image data dataIn, and calculating the gradient of each point in the three-dimensional seismic image to obtain three gradient component bodies g_x、g_yAnd g_z: the gradient at the three-dimensional seismic image point (m, n, k) is calculated as follows:

$\left\{\begin{array}{c}{g}_x(m,n,k)=\frac{dataIn(m+1,n,k)-dataIn(m-1,n,k)}{2};\\ g_y(m,n,k)=\frac{dataIn(m,n+1,k)-dataIn(m,n-1,k)}{2};\\ g_z(m,n,k)=\frac{dataIn(m,n,k+1)-dataIn(m,n,k-1)}{2};\end{array}\right.$

wherein dataIn (m, n, k) represents a pixel value of a point (m, n, k) in the three-dimensional seismic image; g_x(m, n, k) represents the first component of the gradient at point (m, n, k), i.e. the x-direction component in the three-dimensional coordinate system, g_y(m, n, k) represents the second component of the gradient at point (m, n, k), i.e. the y-direction component in the three-dimensional coordinate system, g_z(m, n, k) represents a third component body of the gradient at (m, n, k), i.e., a z-direction component in the three-dimensional coordinate system; dataIn (m +1, n, k), dataIn (m-1, n, k), dataIn (m, n +1, k), dataIn (m, n-1, k), dataIn (m, n, k +1), dataIn (m, n, k-1) are pixel values corresponding to the neighborhood of point (m, n, k);

3) using the gradient component g_x、g_yAnd g_zCalculating a gradient tensor GT of the three-dimensional seismic image data dataIn:

$GT=g\·g^T=\left[\begin{array}{ccc}{g}_x^2& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y^2& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z^2\end{array}\right]=\left[\begin{array}{ccc}{\mathrm{GT}}_{11}& {\mathrm{GT}}_{12}& {\mathrm{GT}}_{13}\\ {\mathrm{GT}}_{21}& {\mathrm{GT}}_{22}& {\mathrm{GT}}_{23}\\ {\mathrm{GT}}_{31}& {\mathrm{GT}}_{32}& {\mathrm{GT}}_{33}\end{array}\right],$

g＝[g_x,g_y,g_z]^T,

wherein, GT₁₁,GT₁₂,GT₁₃,GT₂₂,GT₂₃,GT₃₃Six component bodies of the gradient tensor GT respectively;

4) for the six component bodies GT of the calculated gradient tensor GT₁₁,GT₁₂,GT₁₃,GT₂₂,GT₂₃,GT₃₃Respectively carrying out trilateral filtering to obtain a structure tensor ST of the three-dimensional seismic image;

5) constructing a diffusion tensor DT of the three-dimensional seismic image according to the structure tensor ST, and calculating the flow of each point of the three-dimensional seismic image; wherein: the expression of the diffusion tensor DT is:

$\left\{\begin{array}{c}DT=\mathrm{\ϵ}\·(v_2\·v_2^T+v_3\·v_3^T)\\ \mathrm{\ϵ}=\frac{Tr\left(J_{0}J\right)}{Tr\left(J_0\right)Tr\left(J\right)}\\ \[v_1,v_2,v_3\]=eigv\left(ST\right)\end{array}\right.$

the three component volumes flux _ x, flux _ y, flux _ z of the flux (m, n, t) at the three-dimensional seismic image point (m, n, t) are expressed as:

$flux=\mathrm{\ϵ}\·DT\·\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]$

wherein, take valuesThe range is 0-1; tr (-) represents the trace of the matrix; j. the design is a square₀For the initial three-dimensional seismic image structure tensor, i.e. ST₀(ii) a J is the structure tensor of the current iterative three-dimensional seismic image, namely ST; eigV (-) represents the solution of the eigenvector, and eigenvector v₁,v₂,v₃The corresponding characteristic value satisfies lambda₁>＝λ₂>＝λ₃；flux＝[flux_x，flux_y，flux_z]^T；

6) Calculating the divergence of the flow at the three-dimensional seismic image point to obtain the increment I of each point in the diffusion process in unit time_Δ(ii) a Wherein the divergence is the sum of the differences of the flow components in the corresponding directions;

7) calculating the result of trilateral structure guided filtering of the current iteration according to the time step timestep set in the step 1), namely calculating the output three-dimensional seismic image data dataOut:

dataOut＝dataIn+timeSetp·I_Δ；

8) judging whether the iteration times iterTimes of the trilateral structure guided filtering set in the step 1) are reached, and if so, outputting dataOut; otherwise, using dataOut as the three-dimensional seismic image data to be imported, and returning to the step 2).

In step 1), the filter parameter σ of the structure tensor is calculated₁The value range of (A) is 0.5-5; the value range of the diffusion filtering time step timestep is 0.01-1; the iteration number iterTimes of the trilateral structure guided filtering ranges from 2 to 10.

In step 4), any component ST of the structure tensor ST of the three-dimensional seismic image_ijThe calculation method comprises the following steps:

1) finding a component GT of a gradient tensor GT_ijGradient of (D) to obtain GT_ijThree gradient component bodies GT_ij_g_x，GT_ij_g_y，GT_ij_g_z；GT_ijGradient GT at point (m, n, k)_ijG (m, n, k) is calculated as:

$\left\{\begin{array}{c}{\mathrm{GT}}_{ij}\_g_x(m,n,k)=\frac{{\mathrm{GT}}_{ij}(m+1,n,k)-{\mathrm{GT}}_{ij}(m-1,n,k)}{2};\\ {\mathrm{GT}}_{ij}\_g_y(m,n,k)=\frac{{\mathrm{GT}}_{ij}(m,n+1,k)-{\mathrm{GT}}_{ij}(m,n-1,k)}{2};\\ {\mathrm{GT}}_{ij}\_g_z(m,n,k)=\frac{{\mathrm{GT}}_{ij}(m,n,k+1)-{\mathrm{GT}}_{ij}(m,n,k-1)}{2};\end{array}\right.$

wherein, GT_ij_g(m,n,k)＝[GT_ij_g_x(m,n,k),GT_ij_g_y(m,n,k),GT_ij_g_z(m,n,k)]^T；(i,j)＝{(1,1),(1,2),(1,3),(2,2),(2,3),(3,3)}；

2) According to GT_ij_g_x，GT_ij_g_y，GT_ij_g_zDetermining trilateration filter internal parameters sigma₂：

σ₂＝β·||g_max-g_min||，

Wherein, g_max,g_minAre respectively GT_ijβ is 0.05-0.5, preferably 0.15;

3) according to the parameter σ₁And σ₂Using bilateral filtering to obtain GT_ijLocal mean gradient G of_θ(x)；

G_θ(x) Three component bodies G_θ_x，G_θ_y，G_θThe formula for z is calculated as:

$\left\{\begin{array}{c}{G}_{\mathrm{\θ}}\left(x\right)=\frac{1}{k_{\▿}\left(x\right)}{\∫}_{\mathrm{\Ω}}w\left(p\right)\·c\left(y\right)\·\▿I_{in}(x+p)dp\\ k_{\▿}\left(x\right)={\∫}_{\mathrm{\Ω}}w\left(p\right)\·c\left(y\right)dp\\ y=\left|\right|\▿I_{in}(x+p)-\▿I_{in}\left(x\right)\left|\right|\end{array}\right.,$

wherein, $w\left(p\right)=\frac{1}{\sqrt{2{\mathrm{\π\σ}}_1}}\·\exp (-\frac{p^2}{2{\mathrm{\σ}}_1^2}),c\left(y\right)=\frac{1}{\sqrt{2{\mathrm{\π\σ}}_2}}\·\exp (-\frac{y^2}{2{\mathrm{\σ}}_2^2}),$ x represents the position of a gradient tensor point of the currently computed three-dimensional seismic image; p represents the offset of the pixel point participating in bilateral filtering to the point x; w (p) represents the spatial similarity weight coefficient of the pixel point participating in bilateral filtering to the x point; y represents GT_ijGradient at (x + p) point with GT_ijThe magnitude of the delta of the gradient at the x-point; is a local mean gradient G_θ(x) Normalizing the denominator value in calculation;

4) calculating GT_ijTrilateral filter delta at each point

$\left\{\begin{array}{c}{I}_{\stackrel{\‾}{\mathrm{\Δ}}}\left(x\right)=\frac{{\∫}_{\mathrm{\Ω}}w(x+p)c\left(I_{\mathrm{\Δ}}\right(x+p\left)\right)I_{\mathrm{\Δ}}(x+p)N(x+p)dp}{k_{\▿}^{\′}\left(x\right)}\\ k_{\▿}^{\′}\left(x\right)={\∫}_{\mathrm{\Ω}}w(x+p)c\left(I_{\mathrm{\Δ}}\right(x+p\left)\right)N(x+p)dp\end{array}\right.,$

Wherein, $\left\{\begin{array}{c}{I}_{\mathrm{\&Delta;}}(x+p)=I_{in}(x+p)-I_p(x+p)\\ I_p(x+p)=I_{in}\left(x\right)+G_{\mathrm{\&theta;}}\left(x\right)\&CenterDot;p\end{array}\right.,$ I_p(x + p) is I_pValue at (x + p), I_pIs GT_ijA local tangent plane at point x; i is_Δ(x + p) is GT_ijAt point (x + p) for I_pAn increment of (x + p); n (x + p) is a gradient approximation screening function, $N(x+p)=\{\begin{array}{cc}1,& \left|\right|G_{\mathrm{\&theta;}}(x+p)-G_{\mathrm{\&theta;}}\left(x\right)\left|\right|<R\\ 0,& otherwise\end{array},$ r is a selection threshold, which may be set to R ═ σ₂,For trilateral filtering incrementsNormalizing the denominator value in calculation;

5) combining GT with a tube_ijAnd trilateral filter deltaSuperposing to obtain ST_ij。

Compared with the prior art, the invention has the beneficial effects that: the method combines the trilateral filtering with good edge-preserving effect, greatly improves the performance of the original Gaussian kernel anisotropic diffusion filtering, has good denoising effect and less structural information loss, plays a great role in edge-preserving and denoising of the seismic image, has very outstanding structure-preserving performance, brings great convenience to searching of a target area of the seismic image, and further facilitates interpretation and subsequent processing of the three-dimensional seismic image.

Drawings

FIG. 1 is a block diagram of the process flow of the present invention;

FIG. 2 is a graph comparing the time slicing effects of Gaussian kernel structure guided filtering and trilateral structure guided filtering; the data is selected from data of a Dutch F3 work area disclosed on a network, and FIG. 2(a) is F3 time slice before filtering; FIG. 2(b) shows the result of Gaussian structure-oriented filtering; FIG. 2(c) is the trilateral structure-guided filtering result; FIG. 2(d) shows guided filtering of noise by Gaussian structure; FIG. 2(e) illustrates trilateral guided filtering to remove noise; as can be seen from fig. 2(b) and 2(c), the image structure obtained by the trilateral structure guided filtering is clearer, and as can be seen from fig. 2(d) and 2(e), the loss of structural information of the trilateral structure guided filtering during denoising is very small, and the structure keeping effect is obviously improved compared with that of the traditional gaussian kernel structure guided filtering.

Detailed Description

Referring to fig. 1, the process of the method of the present invention is as follows:

firstly, setting input parameters, namely setting trilateral structure guiding filtering parameters of three-dimensional seismic image data. The parameters being three, including the filter parameter σ for calculating the structure tensor₁Time step of diffusion filtering, timeset, iteration number iterTimes of trilateral structure-oriented filtering. Wherein the structure tensor is calculated by trilateral filtering, wherein the filtering parameter sigma is₁The three-edge filtering parameter has a value range of 0.5-5 generally, and influences the calculation of the structure tensor of the three-dimensional seismic image. The trilateral structure guided filtering of the three-dimensional seismic image is realized through multiple iterations, the increment of the three-dimensional seismic data input during the iteration is calculated during each iteration, the time step length timestep is the increment coefficient and is used for adjusting the increment amplitude and generally takes a value within 0.01-1, when the trilateral structure guided filtering is applied, a small point is taken to be better, such as 0.02,0.05 and 0.1, and the parameter controls the filtering strength of one iteration; the iteration number iterTimes is the filtering degree of the trilateral structure-oriented filtering, and it is sufficient that the iteration number iterTimes is generally 4 or 5.

And secondly, importing three-dimensional seismic image data dataIn and calculating the gradient of each point of the three-dimensional seismic image.

The structural information of the image is transmitted by gradient, and the gradient reflects the direction and intensity of the maximum change between adjacent pixel points on the image. In structured regions of the image, the pixel values of the image vary significantly, and this variation is reflected clearly in the gradient. The trilateral structure-oriented filtering utilizes the structural information to denoise and protect the structural information.

The gradient is specifically calculated by applying forward difference to calculate the gradient of each point to obtain three gradient components, g_x、g_yAnd g_zThe gradient at point (m, n, k) is calculated as follows:

$\left\{\begin{array}{c}{g}_x(m,n,k)=\frac{dataIn(m+1,n,k)-dataIn(m-1,n,k)}{2};\\ g_y(m,n,k)=\frac{dataIn(m,n+1,k)-dataIn(m,n-1,k)}{2};\\ g_z(m,n,k)=\frac{dataIn(m,n,k+1)-dataIn(m,n,k-1)}{2};\end{array}\right.---(2-1)$

in equation 2-1, dataIn (m, n, k) represents a value of a three-dimensional seismic image at a position (m, n, k), i.e., a pixel value, g_x(m, n, k) represents a first component of the gradient at (m, n, k), which can be considered as the x-direction component in a three-dimensional coordinate system, g_y(m, n, k) represents a second component of the gradient at (m, n, k), which may be considered as a y-direction component in a three-dimensional coordinate system, g_z(m, n, k) represents a third component of the gradient at (m, n, k), which can be considered as a z-direction component in the three-dimensional coordinate system.

Their calculation is calculated using the difference of the pixel values of the neighboring points, dataIn (m +1, n, k), dataIn (m-1, n, k), dataIn (m, n +1, k), dataIn (m, n-1, k), dataIn (m, n, k +1), dataIn (m, n, k-1) being the pixel value of the neighboring point of the point (m, n, k).

And thirdly, calculating the gradient tensor according to the calculation result of the second step.

The gradient tensor and the gradient are equivalent in reflecting structure information, but the more accurate local structure information of the image can be obtained by smoothing the gradient tensor (obtaining local average), and the situation that positive and negative directions are cancelled in smoothing can be avoided. Meanwhile, when the gradient tensor is smoothed, a plurality of change direction information of the image local area can be introduced, and the structure information of the image local area is reflected.

The specific calculation formula is as follows:

$GT=g\&CenterDot;g^T=\left[\begin{array}{ccc}{g}_x^2& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y^2& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z^2\end{array}\right]=\left[\begin{array}{ccc}{\mathrm{GT}}_{11}& {\mathrm{GT}}_{12}& {\mathrm{GT}}_{13}\\ {\mathrm{GT}}_{21}& {\mathrm{GT}}_{22}& {\mathrm{GT}}_{23}\\ {\mathrm{GT}}_{31}& {\mathrm{GT}}_{32}& {\mathrm{GT}}_{33}\end{array}\right]---(2-2)$

point-by-point computation, due to the symmetry of the gradient tensor matrix, only six component volumes are needed for storage, GT respectively₁₁,GT₁₂,GT₁₃,GT₂₂,GT₂₃,GT₃₃。

And fourthly, calculating the structure tensor.

As described in the previous step, the structure tensor is generated by gradient tensor smoothing. The method for selecting the smoothness is very critical, and the quality of the image structure tensor directly influences the quality of structure information retention in the trilateral structure guided filtering. The structure tensor is generated by using the trilateral filtering with good edge-preserving effect, and the structure tensor with high quality can be generated.

The specific implementation is that the trilateral filtering is applied to respectively filter the six component bodies of the gradient tensor obtained in the second step to obtain the structure tensor which is also the six component bodies, ST₁₁,ST₁₂,ST₁₃,ST₂₂,ST₂₃,ST₃₃。

With three-edge filteringThe steps are complicated with ST₁₁The calculation of (2) is described in detail as an example, and the generation of other component bodies is the same as this.

Step 4-1, for component GT of gradient tensor₁₁Obtaining a gradient in the same manner as the second step (2-1), obtaining GT₁₁_g_x，GT₁₁_g_y，GT₁₁_g_zThree gradient component volumes. The calculation formula is as follows:

$\{\begin{array}{c}{\mathrm{GT}}_{11}\_g_x(m,n,k)=\frac{{\mathrm{GT}}_{11}(m+1,n,k)-{\mathrm{GT}}_{11}(m-1,n,k)}{2};\\ {\mathrm{GT}}_{11}\_g_y(m,n,k)=\frac{{\mathrm{GT}}_{11}(m,n+1,k)-{\mathrm{GT}}_{11}(m,n-1,k)}{2};\\ {\mathrm{GT}}_{11}\_g_z(m,n,k)=\frac{{\mathrm{GT}}_{11}(m,n,k+1)-{\mathrm{GT}}_{11}(m,n,k-1)}{2};\end{array}---(2-3)$

the equation (2-3) is the same as the equation (2-1), and the input data is different, where the input data is the component GT of the gradient tensor₁₁Calculated is GT₁₁Of the gradient of (c).

Trilateration also protects the structural information of the input data by means of gradients.

Step 4-2, solving parameter sigma according to the component obtained in step 4-1₂。

Parameter sigma₂An internal parameter is used for trilateral filtering, and the parameter is used for calculating the pixel point proximity weight coefficient participating in trilateral filtering. Sigma₂The smaller the value, the higher the required similarity of the pixel points participating in the trilateral filtering.

Specifically, the calculation is firstly searched in the step 4-1Discharged GT₁₁The gradient g of the maximum and minimum modulus is searched out_max,g_minThen calculates σ₂：σ₂＝β·||g_max-g_minAnd l, β is an internal fixed parameter, the selectable value range is 0.05-0.5, and generally β is preferably set as 0.15.

And 4-3, obtaining a local average gradient.

The local gradient of the image is obtained by bilateral filtering.

The specific calculation method is based on the parameter sigma₁And sigma₂To GT₁₁Using bilateral filtering to obtain GT₁₁Local mean gradient G of_θ(x) It consists of three components, G_θ_x，G_θ_y，G_θAnd (iv) is (z). The calculation formula is as follows:

$\left\{\begin{array}{c}{G}_{\mathrm{\&theta;}}\left(x\right)=\frac{1}{k_{\&dtri;}\left(x\right)}{\&Integral;}_{\mathrm{\&Omega;}}w\left(p\right)\&CenterDot;c\left(y\right)\&CenterDot;\&dtri;I_{in}(x+p)dp\\ k_{\&dtri;}\left(x\right)={\&Integral;}_{\mathrm{\&Omega;}}w\left(p\right)\&CenterDot;c\left(y\right)dp\\ y=\left|\right|\&dtri;I_{in}(x+p)-\&dtri;I_{in}\left(x\right)\left|\right|\end{array}\right.---(2-4)$

wherein, $w\left(p\right)=\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}_1}}\&CenterDot;\exp (-\frac{p^2}{2{\mathrm{\&sigma;}}_1^2}),c\left(y\right)=\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}_2}}\&CenterDot;\exp (-\frac{y^2}{2{\mathrm{\&sigma;}}_2^2}),$ x represents the position of the current calculation point, p represents the offset of the pixel point participating in filtering (namely, (x + p) point) to the point x, w (p) represents the spatial similarity weight coefficient of the pixel point participating in filtering (namely, (x + p) point) to the point x, and the smaller p is, the closer p is to the point x, the larger w (p) is; y represents GT₁₁Gradient at (x + p) point for GT₁₁The modulus of the gradient difference at point x, w (p) is the gradient similarity weight coefficient, the smaller y, the closer the gradient at these two points, the larger c (y),is a local mean gradient G_θ(x) The denominator values in the calculation are normalized.

Representing a gradient, i.e. GT₁₁_g_x，GT₁₁_g_y，GT₁₁_g_z. Bilateral filtering may be applied to the three component bodies, respectively, or may be applied to the three component bodies simultaneously.

And 4, calculating the trilateral filtering increment of the input data at each point.

The input data here is trilaterally filtered input data, i.e. the gradient tensor component GT₁₁. Trilateral filtering is accomplished by finding an increment at each point of the input data and superimposing the input data with the increment. The increment calculation is the average of tangent plane increments of each point in the neighborhood, and the average calculation strategy of bilateral filtering is adopted. Such calculation has the effect of protecting the gradient information in the structure information.

Specifically, GT is calculated₁₁Increment per point obtained in filteringTo preserve the gradient information, the increment at each point is the calculated increment of the point in the neighborhood relative to the center tangent plane, calculated as follows:

$\left\{\begin{array}{c}{I}_{\stackrel{\&OverBar;}{\mathrm{\&Delta;}}}\left(x\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;}}w(x+p)c\left(I_{\mathrm{\&Delta;}}\right(x+p\left)\right)I_{\mathrm{\&Delta;}}(x+p)N(x+p)dp}{k_{\&dtri;}^{\&prime;}\left(x\right)}\\ k_{\&dtri;}^{\&prime;}\left(x\right)={\&Integral;}_{\mathrm{\&Omega;}}w(x+p)c\left(I_{\mathrm{\&Delta;}}\right(x+p\left)\right)N(x+p)dp\end{array}\right.---(2-5)$

wherein, $\left\{\begin{array}{c}{I}_{\mathrm{\&Delta;}}(x+p)=I_{in}(x+p)-I_p(x+p)\\ I_p(x+p)=I_{in}\left(x\right)+G_{\mathrm{\&theta;}}\left(x\right)\&CenterDot;p\end{array}\right.,$ I_Δ(x + p) is GT₁₁Value I at (x + p)_in(x + p) for GT at point x₁₁Inclined plane (partial tangent plane) of (1)_pValue I at (x + p)_pIncrement of (x + p), N (x + p) is the gradient approximation screening function, which plays the role of further screening the gradient approximation point, and has $N(x+p)=\{\begin{array}{cc}1,& \left|\right|G_{\mathrm{\&theta;}}(x+p)-G_{\mathrm{\&theta;}}\left(x\right)\left|\right|<R\\ 0,& otherwise\end{array},$ R is a selection threshold, which may be set to R ═ σ₂，For trilateral filtering incrementsThe denominator values in the calculation are normalized.

And 4, step 4-5, overlapping the increment and calculating the trilateral filtering value.

Specifically, for example, ST is calculated₁₁Using data GT before trilateral filtering₁₁Adding the increments obtained in steps 4-4Can obtain ST₁₁。

And fifthly, constructing a diffusion tensor according to the structure tensor calculated in the third step, and calculating the flow at each point.

The structure tensor is used for representing local structure information of an image, so that a diffusion tensor is constructed, trilateral structure-guided filtering can be carried out under the constraint of structure guidance, and structural features are reserved. The diffusion tensor controls the direction of diffusion and has fine control ability.

Specifically, the diffusion tensor and the flow rate at each point are calculated in sequence. At each point, the structure tensor is eigenvalued decomposed:

$ST=\left(\begin{array}{ccc}v1& v2& v3\end{array}\right)\left(\begin{array}{ccc}{\mathrm{\&lambda;}}_1& 0& 0\\ 0& {\mathrm{\&lambda;}}_2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3\end{array}\right)\left(\begin{array}{c}v{1}^T\\ v{2}^T\\ v{3}^T\end{array}\right)---(2-6)$

in the formula (2-6), lambda₁>＝λ₂>＝λ₃。

The diffusion tensor DT is built up based on the structure tensor as follows:

$\left\{\begin{array}{c}DT=\mathrm{\&epsiv;}\&CenterDot;(v_2\&CenterDot;v_2^T+v_3\&CenterDot;v_3^T)\\ \mathrm{\&epsiv;}=\frac{Tr\left(J_{0}J\right)}{Tr\left(J_0\right)Tr\left(J\right)}\end{array}\right.---(2-7)$

in the formula (2-7), the eigenvalue λ of the structure tensor is selected₂,λ₃The corresponding eigenvector is used as the eigenvector of the diffusion tensor, namely the diffusion direction, the layered structure of the seismic image is protected, the adjustment and control are performed by reflecting the continuity of the image, the value range is 0-1, each point of the image has a value, the value of the point is close to 1, the image is relatively continuous at the point, the value of the point is close to 0, the discontinuous structure such as a river channel, a fracture and the like is shown at the point, the diffusion tensor DT can be further adjusted and controlled by the value, the filtering discontinuity at the discontinuous structure is avoided, the fracture, the river channel and other structures are stored, and J₀Is the initial image structure tensor, i.e., GT, J is the structure tensor, i.e., ST, of the current iteration image. Tr (-) denotes the trace of the matrix, i.e. the sum of the diagonal elements.

The flux (m, n, t) at each point is, three component bodies of flux, flux _ x, flux _ y, flux _ z are obtained. The calculation formula is as follows:

$flux=\mathrm{\&epsiv;}\&CenterDot;DT\&CenterDot;\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]---(2-8)$

the gradient g calculated in the second step is used here_x，g_y，g_yAnd the continuity factor and diffusion tensor DT just calculated.

Sixthly, calculating the increment I of each point in unit time_Δ(m,n,t)。

And (4) obtaining the increment of each point in the diffusion process per unit time by calculating the divergence of the flow.

The calculation formula is as follows:

$\left\{\begin{array}{c}{I}_{\mathrm{\&Delta;}}(m,n,t)=\&dtri;flux(m,n,t)=flux\_\mathrm{\&Delta;}x+flux\_\mathrm{\&Delta;}y+flux\_\mathrm{\&Delta;}z;\\ flux\_\mathrm{\&Delta;}x=\frac{flux\_x(m+1,n,t)-flux\_x(m-1,n,t)}{2};\\ flux\_\mathrm{\&Delta;}y=\frac{flux\_y(m,n+1,t)-flux\_y(m,n-1,t)}{2};\\ flux\_\mathrm{\&Delta;}z=\frac{flux\_z(m,n,t+1)-flux\_z(m,n,t-1)}{2};\end{array}\right.---(2-9)$

the calculation formula shows that the divergence is the sum of the differences of the flow components in the corresponding directions.

And seventhly, calculating the result of the current iterative trilateral structure guided filtering according to the set time step length.

The calculation formula is as follows:

dataOut＝dataIn+timeSetp·I_Δ(2-10)

wherein dataOut is the output result, dataIn is the input seismic image, timesep is the set time step, I_ΔThe increment calculated in the fifth step.

And step eight, judging whether the iteration times are reached, if the iteration times are completed, outputting a result, and if not, sending the result generated in the step seven to the step two.

FIG. 2 is a graph comparing the time slicing effect of Gaussian kernel structure guided filtering and trilateral structure guided filtering according to the present invention, and FIG. 2(a) is selected from data of a Netherlands F3 work area disclosed on a network; as can be seen from fig. 2(b) and 2(c), the image structure obtained by the trilateral structure guided filtering is clearer, and as can be seen from fig. 2(d) and 2(e), the loss of structural information of the trilateral structure guided filtering during denoising is very small, and the structure keeping effect is obviously improved compared with that of the traditional gaussian kernel structure guided filtering.

Claims (9)

1. A three-dimensional seismic data image noise reduction method based on trilateral structure guided filtering is characterized by comprising the following steps:

1) setting trilateral structure guiding filtering parameters of the three-dimensional seismic image data: calculating the filter parameter σ of the structure tensor₁The time step of diffusion filtering, namely, the time step of diffusion filtering, and the iteration times of trilateral structure-oriented filtering, namely, iterTimeTimes;

2) importing three-dimensional seismic image data dataIn, and calculating the gradient of each point in the three-dimensional seismic image to obtain three gradient component bodies g_x、g_yAnd g_z: the gradient at the three-dimensional seismic image point (m, n, k) is calculated as follows:

$\left\{\begin{array}{c}{g}_x(m,n,k)=\frac{dataIn(m+1,n,k)-dataIn(m-1,n,k)}{2};\\ g_y(m,n,k)=\frac{dataIn(m,n+1,k)-dataIn(m,n-1,k)}{2};\\ g_z(m,n,k)=\frac{dataIn(m,n,k+1)-dataIn(m,n,k-1)}{2};\end{array}\right.$

wherein dataIn (m, n, k) represents a pixel value of a point (m, n, k) in the three-dimensional seismic image; g_x(m, n, k) represents the first component of the gradient at point (m, n, k), i.e. the x-direction component in the three-dimensional coordinate system, g_y(m, n, k) represents the second component of the gradient at point (m, n, k), i.e. the y-direction component in the three-dimensional coordinate system, g_z(m, n, k) represents a third component body of the gradient at (m, n, k), i.e., a z-direction component in the three-dimensional coordinate system; dataIn (m +1, n, k), dataIn (m-1, n, k), dataIn (m, n +1, k), dataIn (m, n-1, k), dataIn (m, n, k +1), dataIn (m, n, k-1) are pixel values corresponding to the neighborhood of point (m, n, k);

3) using the gradient component g_x、g_yAnd g_zCalculating a gradient tensor GT of the three-dimensional seismic image data dataIn:

$GT=g\&CenterDot;g^T=\left[\begin{array}{ccc}{g}_x^2& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y^2& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z^2\end{array}\right]=\left[\begin{array}{ccc}{\mathrm{GT}}_{11}& {\mathrm{GT}}_{12}& {\mathrm{GT}}_{13}\\ {\mathrm{GT}}_{21}& {\mathrm{GT}}_{22}& {\mathrm{GT}}_{23}\\ {\mathrm{GT}}_{31}& {\mathrm{GT}}_{32}& {\mathrm{GT}}_{33}\end{array}\right],$

g＝[g_x,g_y,g_z]^T,

wherein, GT₁₁,GT₁₂,GT₁₃,GT₂₂,GT₂₃,GT₃₃Six component bodies of the gradient tensor GT respectively;

4) for the six component bodies GT of the calculated gradient tensor GT₁₁,GT₁₂,GT₁₃,GT₂₂,GT₂₃,GT₃₃Respectively carrying out trilateral filtering to obtain a structure tensor ST of the three-dimensional seismic image data dataIn;

5) constructing a diffusion tensor DT of the three-dimensional seismic image according to the structure tensor ST, and calculating the flow of each point of the three-dimensional seismic image; wherein: the expression of the diffusion tensor DT is:

$\left\{\begin{array}{c}DT=\mathrm{\&epsiv;}\&CenterDot;(v_2\&CenterDot;v_2^T+v_3\&CenterDot;v_3^T)\\ \mathrm{\&epsiv;}=\frac{Tr\left(J_{0}J\right)}{Tr\left(J_0\right)Tr\left(J\right)}\\ \&lsqb;v_1,v_2,v_3\&rsqb;=eigV\left(ST\right)\end{array}\right.$

the three component volumes flux _ x, flux _ y, flux _ z of the flux (m, n, t) at the three-dimensional seismic image point (m, n, t) are expressed as:

$flux=\mathrm{\&epsiv;}\&CenterDot;DT\&CenterDot;\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]$

wherein the value range is 0-1; tr (-) represents the trace of the matrix; j. the design is a square₀For the initial three-dimensional seismic image structure tensor, i.e. ST₀(ii) a J is the structure tensor of the current iterative three-dimensional seismic image, namely ST; eigV (-) represents the solution of the eigenvector, and eigenvector v₁,v₂,v₃The corresponding characteristic value satisfies lambda₁>＝λ₂>＝λ₃；flux＝[flux_x，flux_y，flux_z]^T；

6) Calculating the divergence of the flow at the three-dimensional seismic image point to obtain the increment I of each point in the diffusion process in unit time_Δ(ii) a Wherein the divergence is the sum of the differences of the flow components in the corresponding directions;

7) calculating the result of trilateral structure guided filtering of the current iteration according to the time step timestep set in the step 1), namely calculating the output three-dimensional seismic image data dataOut:

dataOut＝dataIn+timeSetp·I_Δ；

8) judging whether the iteration times iterTimes of the trilateral structure guided filtering set in the step 1) are reached, and outputting the iterTimes if the iteration times iterTimes of the trilateral structure guided filtering set in the step 1) are reached; otherwise, using dataOut as the three-dimensional seismic image data to be imported, and returning to the step 2).

2. The method for denoising three-dimensional seismic data images based on trilateral structure-oriented filtering according to claim 1, wherein in step 1), the filtering parameter σ of the structure tensor is calculated₁The value range of (A) is 0.5-5.

3. The three-dimensional seismic data image noise reduction method based on trilateral structure-oriented filtering according to claim 2, wherein in the step 1), the time step timestep of diffusion filtering ranges from 0.01 to 1.

4. The three-dimensional seismic data image noise reduction method based on trilateral structure-oriented filtering according to claim 3, wherein in the step 1), the range of iterTimes of iteration times of trilateral structure-oriented filtering is 2-10.

5. The method for denoising three-dimensional seismic data image based on trilateral structure-oriented filtering according to one of claims 1 to 4, wherein in step 4), any component ST of the structure tensor ST of the three-dimensional seismic image data dataIn_ijThe calculation method comprises the following steps:

1) finding a component GT of a gradient tensor GT_ijGradient of (D) to obtain GT_ijThree gradient component bodies GT_ij_g_x，GT_ij_g_y，GT_ij_g_z；GT_ijGradient GT at point (m, n, k)_ijG (m, n, k) is calculated as:

$\left\{\begin{array}{c}{\mathrm{GT}}_{ij}\_g_x(m,n,k)=\frac{{\mathrm{GT}}_{ij}(m+1,n,k)-{\mathrm{GT}}_{ij}(m-1,n,k)}{2};\\ {\mathrm{GT}}_{ij}\_g_y(m,n,k)=\frac{{\mathrm{GT}}_{ij}(m,n+1,k)-{\mathrm{GT}}_{ij}(m,n-1,k)}{2};\\ {\mathrm{GT}}_{ij}\_g_z(m,n,k)=\frac{{\mathrm{GT}}_{ij}(m,n+1,k)-{\mathrm{GT}}_{ij}(m,n-1,k)}{2};\end{array}\right.$

wherein, GT_ij_g(m,n,k)＝[GT_ij_g_x(m,n,k),GT_ij_g_y(m,n,k),GT_ij_g_z(m,n,k)]^T；(i,j)＝{(1,1),(1,2),(1,3),(2,2),(2,3),(3,3)}；

2) According to GT_ij_g_x，GT_ij_g_y，GT_ij_g_zDetermining trilateration filter internal parameters sigma₂：

σ₂＝β·||g_max-g_min||，

Wherein, g_max,g_minAre respectively GT_ijβ the range of value is 0.05-0.5;

3) according to the parameter σ₁And σ₂Using bilateral filtering to obtain GT_ijLocal mean gradient G of_θ(x)；

4) Calculating GT_ijTrilateral filter delta at each point

5) Combining GT with a tube_ijAnd trilateral filter deltaSuperposing to obtain ST_ij。

6. The method of claim 5, wherein β is 0.15.

7. The method of claim 6, wherein GT is a method of denoising three-dimensional seismic data images based on trilateral-guided filtering_ijLocal mean gradient G of_θ(x) Three component bodies G_θ_x，G_θ_y，G_θThe formula for z is calculated as:

$\left\{\begin{array}{c}{G}_{\mathrm{\&theta;}}\left(x\right)=\frac{1}{k_{\&dtri;}\left(x\right)}{\&Integral;}_{\mathrm{\&Omega;}}w\left(p\right)\&CenterDot;c\left(y\right)\&CenterDot;\&dtri;I_{in}(x+p)dp\\ k_{\&dtri;}\left(x\right)={\&Integral;}_{\mathrm{\&Omega;}}w\left(p\right)\&CenterDot;c\left(y\right)dp\\ y=\left|\right|\&dtri;I_{in}(x+p)-\&dtri;I_{in}\left(x\right)\left|\right|\end{array}\right.,$

wherein, $w\left(p\right)=\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}_1}}\&CenterDot;\exp (-\frac{p^2}{2{\mathrm{\&sigma;}}_1^2}),c\left(y\right)=\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}_2}}\&CenterDot;\exp (-\frac{y^2}{2{\mathrm{\&sigma;}}_2^2}),$ x represents the position of a gradient tensor point of the currently computed three-dimensional seismic image; p represents the offset of the pixel point participating in bilateral filtering to the point x; w (p) represents the spatial similarity weight coefficient of the pixel point participating in bilateral filtering to the x point; y represents GT_ijGradient at (x + p) point with GT_ijThe magnitude of the delta of the gradient at the x-point; is a local mean gradient G_θ(x) The denominator values in the calculation are normalized.

8. The method of claim 7 in which GT is a three-dimensional seismic data image denoising method based on trilateral-guided filtering_ijTrilateral filter delta at each pointThe calculation formula of (A) is as follows:

$\left\{\begin{array}{c}{I}_{\stackrel{\&OverBar;}{\mathrm{\&Delta;}}}\left(x\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;}}w(x+p)c\left(I_{\mathrm{\&Delta;}}\right(x+p\left)\right)I_{\mathrm{\&Delta;}}(x+p)N(x+p)dp}{k_{\&dtri;}\left(x\right)}\\ k_{\&dtri;}^{\&prime;}\left(x\right)={\&Integral;}_{\mathrm{\&Omega;}}w(x+p)c\left(I_{\mathrm{\&Delta;}}\right(x+p\left)\right)N(x+p)dp\end{array}\right.,$

wherein, $\left\{\begin{array}{c}{I}_{\mathrm{\&Delta;}}(x+p)=I_{in}(x+p)-I_p(x+p)\\ I_p(x+p)=I_{in}\left(x\right)+G_{\mathrm{\&theta;}}\left(x\right)\&CenterDot;p\end{array}\right.,$ I_p(x + p) is I_pValue at (x + p), I_pIs GT_ijA local tangent plane at point x; i is_Δ(x + p) is GT_ijAt point (x + p) for I_pAn increment of (x + p); n (x + p) is a gradient approximation screening function, $N(x+p)=\left\{\begin{array}{cc}1,& \left|\right|G_{\mathrm{\&theta;}}(x+p)-G_{\mathrm{\&theta;}}\left(x\right)\left|\right|<R\\ 0,& otherwise\end{array}\right.,$ r is the selection threshold, and R is the selection threshold,for trilateral filtering incrementsThe denominator values in the calculation are normalized.

9. The method of claim 8, wherein the selection threshold R- σ is selected based on a trilateral guided filtering based three-dimensional seismic data image denoising method₂。