CN103926616B - A kind of multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection - Google Patents

Abstract

The present invention relates to a kind of multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection, it includes step: 1) the prestack CRP road collection of input is carried out Regularization；2) adopt two dimension Mallat algorithm to carry out multi-resolution decomposition, and initialize every sub-section after decomposition；3) try to achieve diffusion threshold value through diffusion coefficient method, enter step 6)；After carrying out trying to achieve diffusion tensor parameter based on the anisotropy parameter of diffusion tensor method, substitute into Anisotropic Nonlinear equation and carry out anisotropic diffusion filtering, and enter step 4)；4) the sub-section after each iteration is calculated SNR, MSE and PSNR, it is preferable that go out the sub-section of earthquake of the best；5) each sub-section after anisotropic diffusion filtering being processed carries out step 4) process, enters step 6)；6) adopt two dimension Mallat restructing algorithm that each sub-profile information after iteration carries out lossless reconstruction prestack CRP road collection, after reconstruct, export the prestack CRP road collection of the best.Present invention can be widely used in the processing procedure of various oil exploration seismic data.

Description

A kind of multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection

Technical field

The present invention relates to the seismic data processing technique in a kind of oil exploration process, especially with regard to the one multiple dimensioned anisotropic diffusion filtering method based on prestack CRP (common reflection point) road collection.

Background technology

In actual production, along with the raising of accuracy of seismic exploration, people are more and more higher to the requirement of seismic data quality." three high " (high-resolution, high s/n ratio, the high fidelity) treatment technology of seismic data is always up focus and the difficult point of field of geophysical exploration research, wherein, signal to noise ratio governs the raising of resolution, and fidelity directly affects the reliability of data interpretation achievement.Along with exploration targets changes to lithology, hidden-type oil-gas reservoir gradually from conventional oil gas reservoir, prestack inversion becomes important means, the major influence factors of this means is signal to noise ratio, therefore improves the success rate of Subtle reservoir exploration, and one of key is to improve the signal to noise ratio of prestack CRP road set information.Effective noise-removed filtering method while compacting noise, should retain effective lineups information better.

Anisotropic diffusion filtering is a kind of Noise Elimination method, it is possible to be equivalent to the solution procedure of a physical thermal conduction.Partial differential equation application in image procossing is derived by classical physical problem, and PDE model is initially applied to image denoising.Assuming that piece image is the different two-dimensional flat plate of one piece of heat, on plate, the temperature of zones of different is different.The pixel that in the frequency domains such as minutia in image, noise, value is bigger can be regarded as high temperature dot, and the low frequency pixel of homogeneity, smooth region is considered as low warm spot.Image denoising process is considered as the process making entire plate temperature tend to balance: high temperature dot to low warm spot carry energy, high temperature dot cooling and low warm spot accept energy heat up.When the process of energy transmission terminates, the temperature in entire plate realizes balance.In seismic profile, the process of noise compacting is just as temperature equilibrating process, is processed the signal to noise ratio of raising data by diffusing filter.

Perona et al. proposes P-M model (nonlinear anisotropic diffusion equation), and it remains the boundary information of image, but its well-posedness and stability have problems, and do not consider the structural orientation of topography.Weickert introduces diffusion tensor and replaces diffusion coefficient to extract the directivity information of picture structure, improve P-M model in edge details, lose the defect of effective information, the method proposing again later there is the structure diffusion tensor of rotational invariance, overcome the restriction of time step, but structure masterplate is complex.

Summary of the invention

For the problems referred to above, it is an object of the invention to provide a kind of prestack CRP road of can effectively suppressing and concentrate random noise and improve the multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection of data signal to noise ratio quality.

For achieving the above object, the present invention takes techniques below scheme: a kind of multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection, it comprises the following steps: 1) the prestack CRP road collection of input is carried out Regularization, utilizes the gaussian kernel function G of different variances sigma value_σWith seismic profile convolution；2) adopt two dimension Mallat algorithm that the seismic profile after Regularization is done multi-resolution decomposition, using symmetrical wavelet sym6 as mother wavelet, original image is carried out N shell wavelet decomposition；And every sub-section after decomposing is initialized；3) each sub-section of multi-resolution decomposition after initialization is tried to achieve diffusion threshold value through diffusion coefficient method, enter step 5)；Meanwhile, after being undertaken each sub-section of multi-resolution decomposition after initializing trying to achieve diffusion tensor parameter based on the anisotropy parameter of diffusion tensor method, substitution Anisotropic Nonlinear equation carries out anisotropic diffusion filtering, and enters step 4)；4) the sub-section that iteration each time is obtained, calculates signal to noise ratio, mean square deviation and Y-PSNR, finally determines best iterations, obtains the sub-section of earthquake of the best according to this iterations；5) adopt two dimension Mallat restructing algorithm that each sub-profile information after iteration carries out lossless reconstruction prestack CRP road collection, after reconstruct, export the prestack CRP road collection of the best；Wherein best prestack CRP road collection is:

$C_j=H_r^{*}{H}_c^{*}{C}_{j+1}+H_r^{*}{H}_c^{*}{D}_{j+1}^{*}+H_r^{*}{H}_c^{*}{D}_{j+1}^2+H_r^{*}{H}_c^{*}{D}_{j+1}^3,$

In formula, C_j+1For low frequency component,For horizontal high frequency component,For vertical high frequency component,For diagonal angle high fdrequency component,WithRespectively representative with scaling filter coefficient to arrayThe conjugation of operator of row and column effect.

Described step 2) in, described multi-resolution decomposition comprises the following steps: (1) setsFor picture signal to be analyzed, its two dimension is approached image and is: $A_{j}f=A_{j+1}f+D_{j+1}^{1}f+D_{j+1}^{2}f+D_{j+1}^{3}f,---\left(1\right)$ Wherein,

Wherein, A_jFor the result of the jth layer after decomposing, C_j+1(m, n) for low frequency component,For horizontal high frequency component,For vertical high frequency component,For diagonal angle high fdrequency component,For wavelet basis function, m is the translational movement n of wave filter in column direction is the translational movement by line direction wave filter；Utilize the orthogonality of scaling function and wavelet function, formula (1) and (2) obtain:

$C_{j+1}(m,n)=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)h(l-2n)C_j(k,l)$

$\left\{\begin{array}{c}{D}_{j+1}^1=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)g(l-2n)C_j(k,l)\\ D_{j+1}^2=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)h(l-2n)C_j(k,l)\\ D_{j+1}^3=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)g(l-2n)C_j(k,l)\end{array}\right.$

To low frequency component C_j+1, horizontal high frequency componentVertical high frequency componentWith diagonal angle high fdrequency componentCarry out two dimension Mallat to decompose:

$\left\{\begin{array}{c}{C}_{j+1}=H_r{H}_c{C}_j\\ D_{j+1}^1=H_r{G}_c{C}_j\\ D_{j+1}^2=G_r{H}_c{C}_j\\ D_{j+1}^3=G_r{G}_c{G}_j\end{array}\right.,j=\mathrm{0,1},...,J$

In formula, H_rAnd H_cRespectively expression with scaling filter coefficient to arrayThe operator of row and column effect, G_rAnd G_cThe expression wavelet filter coefficient operator to row and column effect respectively；(2) effective information of the seismic section image after decomposing is concentrated mainly on low frequency component C_j+1In, directly to low frequency component C_j+1Carry out anisotropic diffusion filtering process；For the high fdrequency component of seismic section image, it is carried out threshold value shrink process；Threshold value chooses formula:

$\mathrm{\η}=\mathrm{\δ}\·2\sqrt{\log N}/\sqrt{\log (l+1)}\mathrm{\δ}=\mathrm{median}\left(u\right)/0.6745$

Wherein, η is the wavelet threshold set, and the standard deviation that δ is high fdrequency component is estimated, u is image wavelet component, and median represents and takes mediant, and l is Decomposition order, and N is pixel number in component.

In described step 3), described diffusion tensor method comprises the following steps: (1) set diffusing filter equation as:

$\frac{\∂u(x,y,t)}{\∂t}=\▿\·(D\left(\right|{\▿u}_{\mathrm{\σ}}\left|\right)\▿u)$

Wherein,For diffusion tensor, σ is scaling function,For rightCarry out Regularization, wherein

${\▿u}_{\mathrm{\σ}}=\▿(K_{\mathrm{\σ}}*u(x,y,t)),K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\σ}}^2})$

In formula, K_σRepresent that yardstick is the Gaussian function of σ；(2) to each sub-section of multi-resolution decomposition after initializing along the derivation of room and time direction, it is determined that fabric tensor, and solve the eigenvalue of fabric tensor and unit character is vectorial, to determine seismic profile local image structure: 1. fabric tensor T for:Assume fabric tensor

$T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right),$

T is the symmetric positive definite matrix of 2 × 2；2. adopting the Gauss-Seidel method with accelerated factor that fabric tensor T is solved, the eigenvalue obtaining fabric tensor T matrix is:

$u_1=\frac{1}{2}(T_{11}+T_{12}+\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

$u_2=\frac{1}{2}(T_{11}+T_{12}-\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

Unit character vector v₁, v₂It is parallel to gradient direction and unit character vector v₁Meet

$\left(\begin{array}{c}\cos \mathrm{\α}\\ \sin \mathrm{\α}\end{array}\right)=v_1\left|\right|\left(\begin{array}{c}2T_{12}\\ T_{22}-T_{11}+\sqrt{{(T_{11}-T_{22})}^2+{{4T}_{12}}^2}\end{array}\right);$

(3) according to the characteristic vector of the fabric tensor T matrix tried to achieve in step (2) build diffusion tensor D, diffusion tensor D be also symmetric positive definite matrix and its characteristic vector identical with the characteristic vector of fabric tensor T, the eigenvalue of diffusion tensor D is λ₁And λ₂, solve the parameter in diffusion tensor D matrix: diffusion tensor D is:

$D=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right);$

(4), after obtaining diffusion tensor parameter a, b, c, d, it is achieved solving of Anisotropic Nonlinear equation, diffusion tensor equation is as follows:

$\frac{\∂u}{\∂t}=\mathrm{div}(D\·{\▿u}^k)---\left(3\right)$

The left side for diffusion tensor equation (3) solves, and replaces differential to obtain the discrete form of diffusion tensor equation by difference:

$\frac{\∂u}{\∂t}=\frac{u^{k+1}-u^k}{\mathrm{\Δt}}$

U in formula^kAnd u^k+1Respectively section result after k time and k+1 iteration, Δ t is step-length；Diffusion tensor equation (3) the right has been solved:

$\begin{array}{c}\mathrm{div}(D\·{\▿u}^k)=\left(\begin{array}{cc}{\∂}_x& {\∂}_y\end{array}\right)\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)\left(\begin{array}{c}{\∂}_x{u}^k\\ {\∂}_y{u}^k\end{array}\right)\\ ={\∂}_x\left({a\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left({b\∂}_x{u}^k\right)+{\∂}_y\left({c\∂}_y{u}^k\right)\end{array}$

Diffusion tensor equation (1) the right and left simultaneous is got up for:

$u^{k+1}=u^k+\mathrm{\Δt}({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left(b{\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right))---\left(4\right)$

Wherein, and Δ t ∈ (0,0.3], realize each layer section after decomposing is diffused Filtering Processing by formula (4), the random noise that compacting prestack CRP road is concentrated strengthens the seriality of lineups；(5) root-mean-square mutation operator Θ is built according to fabric tensor T and diffusion tensor D:

$\mathrm{\Θ}=\sqrt{\frac{\mathrm{\Σ}\left(\mathrm{diag}\left(T^{0}T\right)\right)}{\mathrm{\Σ}\left(\mathrm{diag}\left(T^0\right)\right)\·\mathrm{\Σ}\left(\mathrm{diag}\left(T\right)\right)}};$

Wherein, diag () is matrix T⁰The leading diagonal of T, T⁰Initial construction tensor, T is current iteration fabric tensor, Θ ∈ [0,1], and near sudden change, Θ levels off to 0；During away from sudden change, Θ levels off to 1；Being added in diffusion tensor equation (3) by root-mean-square mutation operator Θ, obtaining Anisotropic Nonlinear equation is:

$u^{k+1}=u^k+\mathrm{\Δt\Θ}\·({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right)).$

In described step (3), four parameter a in described diffusion tensor D matrix, b, c and d solve as follows: 1. solve the eigenvalue λ of diffusion tensor D₁And λ₂, require to adopt sudden change enhancing coefficient method and the relevant two kinds of methods of coefficient method that strengthen to estimate eigenvalue λ according to actual seismic profile₁And λ₂: sudden change strengthens coefficient method: take eigenvalue λ₂=1；Characteristic vector v₁Characteristic of correspondence value λ₁For fabric tensor T eigenvalue u₁Monotonous descending function, and

λ₂=1

Relevant enhancing coefficient method: make eigenvalue λ₁And λ₂Value is as follows:

λ₁=c₁

Wherein,For the class natural impedance factor；2. eigenvalue λ₁And λ₂Solve the parameter a of diffusion tensor D, b, c and d:

$D={(v_1,v_2)}^T\·\left(\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right)\·(v_1,v_2)=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)$

$\left\{\begin{array}{c}a={\mathrm{\λ}}_1{\cos }^{2}a+{\mathrm{\λ}}_2{\sin }^{2}a\\ b=({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)\cos a\sin a\\ c={\mathrm{\λ}}_1{\sin }^{2}a+{\mathrm{\λ}}_2{\cos }^{2}a\end{array}\right..$

In described step 3), described diffusion coefficient method comprises the following steps: (1) sets diffusion coefficient: using Biweight norm as diffusion coefficient

In formula, k is diffusion threshold value；(2) adopting neighborhood inside gradient absolute deviation method of estimation to solve diffusion threshold value according to diffusion coefficient, diffusion thresholding k (l) is relevant to Decomposition order l:

$k\left(l\right)=\frac{1}{0.6745}\·\mathrm{median}\left(\right||\▿u\left(l\right)-\mathrm{median}\left(\right||\▿u\left(l\right)|\left|\right)|\left|\right)\sqrt{2\log (l+1)/l}.$

Due to the fact that and take above technical scheme, it has the advantage that 1, due to the fact that employing soft-threshold small echo processes high fdrequency component, therefore the damage to section effective information can be reduced, improve the accuracy of abrupt information on diffusion process midship section, namely can distinguish the gray-value variation caused by noise and the gray-value variation caused of suddenling change in strong noise environment.2, the present invention proposes the diffusion equation with specific response characteristic according to the feature of pending section, and seek the balance between noise reduction and details maintenance, therefore each layer section after decomposing can be diffused Filtering Processing, the random noise that compacting prestack road is concentrated, strengthen the seriality of lineups simultaneously, the internal memory taken is few, and computational efficiency is high.3, the present invention is by improving diffusion parameter and the condition of convergence adaptability to different statistical property sections, namely estimates the parameter such as diffusion thresholding, diffusion iterations adaptively, is therefore effectively improved the robustness of diffusion result.4, seismic profile is processed by the application multiple dimensioned anisotropic diffusion filtering method of the present invention, not only makes noise effectively be suppressed, and the signal to noise ratio of data improves；And effectively enhance the seriality of lineups.5, due to the fact that introducing root-mean-square mutation operator, make abrupt information be retained, be therefore really achieved the effect being effectively retained original effective information and compacting noise.6, the present invention is in seismic data processing practically, can strengthen lineups seriality, can keep again the abrupt information that prestack CRP road is concentrated, and strengthens the identification ability to abrupt information；Section seismic wave characteristic after adopting the present invention to process is natural, and signal to noise ratio is improved, and lateral continuity improves, and tectonic information is abundanter.Present invention can be widely used in the processing procedure of various oil exploration seismic data.

Accompanying drawing explanation

Fig. 1 is anisotropic diffusion filtering schematic flow sheet of the present invention

Fig. 2 is the result schematic diagram adopting prior art Focus denoising module to filter random noise

Fig. 3 is the result schematic diagram adopting the present invention to filter random noise

Fig. 4 is that the section lineups information after adopting filtering noise of the present invention is with other method comparing result schematic diagram

Fig. 5 is the curve synoptic diagram that interface, top, Fig. 4 Zhong Ge CRP Dao Ji coal seam amplitude is followed the trail of

Fig. 6 is the stacked section contrast schematic diagram of the CRP road collection different angles before and after denoising

Fig. 7 is that somewhere A well stream journey adopts prior art to optimize the VP/VS inverting generalized section of forward and backward pre-stack data

Fig. 8 is that somewhere A well stream journey adopts the present invention to optimize the VP/VS inverting generalized section of forward and backward pre-stack data

Fig. 9 is based on old process and the present invention optimizes the work area VP/VS inverse time section contrast schematic diagram that forward and backward data obtains

Detailed description of the invention

Below in conjunction with drawings and Examples, the present invention is described in detail.

As it is shown in figure 1, the multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection of the present invention comprises the following steps:

1) the prestack CRP road collection of input is carried out Regularization, namely utilize the gaussian kernel function G of different variances sigma value_σ(σ controls smoothness, and the more big image of σ value is more smooth, chooses according to demand in reality) and seismic profile convolution, so can eliminate part singular value, obtain flattened section so that subsequent treatment is more stable.

2) two dimension Mallat(horse traction is adopted) seismic profile after Regularization does multi-resolution decomposition by algorithm, using symmetrical wavelet sym6 as mother wavelet, original image is carried out N shell wavelet decomposition；And every sub-section after decomposing is initialized, such as given iterations and iteration step length.

Wherein, multi-resolution decomposition comprises the following steps:

(1) setFor picture signal to be analyzed, its two dimension is approached image and is:

$A_{j}f=A_{j+1}f+D_{j+1}^{1}f+D_{j+1}^{2}f+D_{j+1}^{3}f,---\left(1\right)$

Wherein,

Wherein, A_jFor the result of the jth layer after decomposing, C_j+1(m, n) for low frequency component,For horizontal high frequency component,For vertical high frequency component,For diagonal angle high fdrequency component,For wavelet basis function, m is the translational movement n of wave filter in column direction is the translational movement by line direction wave filter.

Utilize the orthogonality of scaling function and wavelet function, formula (1) and (2) can obtain:

$C_{j+1}(m,n)=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)h(l-2n)C_j(k,l)---\left(3\right)$

$\left\{\begin{array}{c}{D}_{j+1}^1=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)g(l-2n)C_j(k,l)\\ D_{j+1}^2=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)h(l-2n)C_j(k,l)\\ D_{j+1}^3=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)g(l-2n)C_j(k,l)\end{array}\right.---\left(4\right)$

To low frequency component C_j+1, horizontal high frequency componentVertical high frequency componentWith diagonal angle high fdrequency componentCarry out two dimension Mallat to decompose:

$\left\{\begin{array}{c}{C}_{j+1}=H_r{H}_c{C}_j\\ D_{j+1}^1=H_r{G}_c{C}_j\\ D_{j+1}^2=G_r{H}_c{C}_j\\ D_{j+1}^3=G_r{G}_c{G}_j\end{array}\right.,j=\mathrm{0,1},...,J---\left(5\right)$

In formula, H_rAnd H_cRespectively expression with scaling filter coefficient to arrayThe operator of row and column effect, G_rAnd G_cThe expression wavelet filter coefficient operator to row and column effect respectively.

(2) effective information of the seismic section image after decomposing is concentrated mainly on low frequency component C_j+1In, can directly to low frequency component C_j+1Carry out anisotropic diffusion filtering process (i.e. diffusion tensor method)；For the high fdrequency component of seismic section image, mainly some trickle catastrophe point or random noises, therefore to carry out threshold value shrink process (i.e. diffusion coefficient method) to it.Wherein threshold value is chosen formula and is:

$\mathrm{\η}=\mathrm{\δ}\·2\sqrt{\log N}/\sqrt{\log (l+1)},\mathrm{\δ}=\mathrm{median}\left(u\right)/0.6745---\left(6\right)$

Wherein, η is the wavelet threshold set, and the standard deviation that δ is high fdrequency component is estimated, u is image wavelet component, and median represents and takes mediant, and l is Decomposition order, and N is pixel number in component.

The threshold value contraction method chosen is that soft-threshold small echo processes, and so can reduce the damage to section effective information, then soft-threshold sizeFor:

$\hat{u}=\left\{\begin{array}{cc}\mathrm{sgn}\left(u_j\right)\left(\right|u_j-\mathrm{\η}\left|\right)& \left|u_j\right|\≥\mathrm{\η}\\ 0& \left|u_j\right|\≤\mathrm{\η}\end{array}\right..---\left(7\right)$

3) each sub-section of multi-resolution decomposition after initialization is tried to achieve diffusion threshold value through diffusion coefficient method, enter step 5)；Meanwhile, after being undertaken each sub-section of multi-resolution decomposition after initializing trying to achieve diffusion tensor parameter based on the anisotropy parameter of diffusion tensor method, substitution Anisotropic Nonlinear equation carries out anisotropic diffusion filtering, and enters step 4)；Wherein Anisotropic Nonlinear equation is:

$\frac{\∂u}{\∂t}=\mathrm{div}[c\left(\right||\▿(G_{\mathrm{\σ}}*u_k)|\left|\right)\·\▿u_k],---\left(8\right)$

In formula, u_kFor sub-section,For image gradient, G_σRepresenting that variance is the gaussian kernel function of σ, section has smoothing effect, div is divergence.

4) the sub-section that iteration each time is obtained, calculates signal to noise ratio snr, mean square deviation MSE and Y-PSNR PSNR, finally determines best iterations, obtains the sub-section of earthquake of the best according to this iterations.

5) adopt two dimension Mallat restructing algorithm that each sub-profile information after iteration carries out lossless reconstruction prestack CRP road collection, the prestack CRP road collection of the best after reconstruct, can be exported；Wherein best prestack CRP road collection is:

$C_j=H_r^{*}{H}_c^{*}{C}_{j+1}+H_r^{*}{H}_c^{*}{D}_{j+1}^{*}+H_r^{*}{H}_c^{*}{D}_{j+1}^2+H_r^{*}{H}_c^{*}{D}_{j+1}^3---\left(9\right)$

In formula, C_j+1For low frequency component,For horizontal high frequency component,For vertical high frequency component,For diagonal angle high fdrequency component,WithRespectively representative with scaling filter coefficient to arrayThe conjugation of operator of row and column effect.

Above-mentioned steps 3) in, the filtering method based on diffusion tensor is according to part section information architecture fabric tensor, then is obtained diffusion tensor by fabric tensor, and diffusion tensor method comprises the following steps:

(1) set diffusing filter equation as:

$\frac{\∂u(x,y,t)}{\∂t}=\▿\·(D\left(\right|{\▿u}_{\mathrm{\σ}}\left|\right)\▿u)---\left(10\right)$

Wherein,For diffusion tensor, σ is scaling function, it is possible to the mean square deviation Automatic adjusument according to noise,For rightCarry out Regularization, be equivalent to effect gaussian filtering, be converted into a well-posed problem, as follows:

${\▿u}_{\mathrm{\σ}}=\▿(K_{\mathrm{\σ}}*u(x,y,t)),K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\σ}}^2})---\left(11\right)$

In formula, K_σRepresent that yardstick is the Gaussian function of σ.

(2) to each sub-section of multi-resolution decomposition after initializing along the derivation of room and time direction, and then determine fabric tensor, and solve eigenvalue and the unit character vector of fabric tensor, to determine seismic profile local image structure:

1. fabric tensor T is:

$T(\▿u)=\▿u_{\mathrm{\σ}}\⊗{\▿u}_{\mathrm{\σ}}---\left(12\right)$

Assume fabric tensor $T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right),$ T is the symmetric positive definite matrix of 2 × 2；

2. owing to fabric tensor T is symmetric positive definite matrix, therefore adopt Gauss-Seidel(Gauss with accelerated factor-Saden you) fabric tensor T solves by method, the advantage adopting Gauss-Seidel method is in that: as long as fabric tensor T is symmetric positive definite, must restrain；SolvingProcess is constantly usedSubstitute in timeNumber before, the internal memory so taken is few, and computational efficiency is high；Its solution procedure is as follows:

$x_i^{k+1}=(b_i-\underset{j=1}{\overset{i-1}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^{k+1}-\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^k)/a_{\mathrm{ii}},i=\mathrm{1,2},...n---\left(13\right)$

In order to improve computational efficiency, introduce best accelerating factor α in an iterative process:

$x_i^{k+1}=x_i^k+\mathrm{\α}(x_i^{k+1}-x_i^k),\mathrm{\α}=\frac{2}{1+\sqrt{1-{\mathrm{\ρ}}^2\left(T\right)}}---\left(14\right)$

Wherein, ρ is the spectral radius of Iterative Matrix.

By solving it can be seen that the eigenvalue of fabric tensor T matrix is:

$\begin{array}{c}{u}_1=\frac{1}{2}(T_{11}+T_{12}+\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})\\ u_2=\frac{1}{2}(T_{11}+T_{12}-\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})\end{array}---\left(15\right)$

Unit character vector v₁, v₂It is parallel to gradient direction and unit character vector v₁Meet

$\left(\begin{array}{c}\cos \mathrm{\α}\\ \sin \mathrm{\α}\end{array}\right)=v_1\left|\right|\left(\begin{array}{c}2T_{12}\\ T_{22}-T_{11}+\sqrt{{(T_{11}-T_{22})}^2+{{4T}_{12}}^2}\end{array}\right)---\left(16\right)$

u₁、u₂For two eigenvalues of fabric tensor T matrix, and there is eigenvalue u₁Be more than or equal to eigenvalue u₂, therefore eigenvalue u₁Direction will point to the direction that data fluctuations is bigger, i.e. the direction of vertical prestack CRP road collection lineups, the abrupt information that instruction road is concentrated；Eigenvalue u₂Point to the direction that data fluctuations is less, be namely parallel to the direction of lineups.v₁、v₂For two characteristic vectors of fabric tensor T matrix, mutually orthogonal between the two, characteristic vector v₁With eigenvalue u₁Corresponding, its direction is also vertical with lineups direction, characteristic vector v₂With eigenvalue u₂Corresponding, it is parallel to lineups direction.

(3) according to the characteristic vector of the fabric tensor T matrix tried to achieve in step (2) build diffusion tensor D, diffusion tensor D be also symmetric positive definite matrix and its characteristic vector identical with the characteristic vector of fabric tensor T, the eigenvalue of diffusion tensor D is λ₁And λ₂, and solve the parameter in diffusion tensor D matrix.Wherein, diffusion tensor D is:

$D=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right).$

Four parameter a in diffusion tensor D matrix, b, c and d solve as follows:

1 zero eigenvalue λ solving diffusion tensor D₁And λ₂, require to adopt sudden change enhancing coefficient method and the relevant two kinds of methods of coefficient method that strengthen to estimate eigenvalue λ according to actual seismic profile₁And λ₂, sudden change Enhancement Method can retain more details abrupt information, and relevant Enhancement Method can make the seriality of prestack CRP road collection lineups strengthen.

Sudden change strengthens coefficient method: in order to reach to strengthen the purpose of sudden change detailed information, take eigenvalue λ₂=1；Characteristic vector v₁Characteristic of correspondence value λ₁For fabric tensor T eigenvalue u₁Monotonous descending function, and

λ₂=1

Relevant enhancing coefficient method: in order to reach to strengthen the successional purpose of lineups, make eigenvalue λ₁And λ₂Value is as follows:

λ₁=c₁

Wherein,For the class natural impedance factor, expression 2-D data has the character of invariable rotary in a certain direction, avoids the process finding optimum invariable rotary matrix, makes diffusion process be greatly simplified；

2. eigenvalue λ₁And λ₂Solve the parameter a of diffusion tensor D, b, c and d:

$D={(v_1,v_2)}^T\·\left(\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right)\·(v_1,v_2)=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)---\left(19\right)$

$\left\{\begin{array}{c}a={\mathrm{\λ}}_1{\cos }^{2}a+{\mathrm{\λ}}_2{\sin }^{2}a\\ b=({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)\cos a\sin a\\ c={\mathrm{\λ}}_1{\sin }^{2}a+{\mathrm{\λ}}_2{\cos }^{2}a\end{array}\right.---\left(22\right)$

Eigenvalue λ₁、λ₂It is used to control diffusion, λ₁It is a less arithmetic number, represents that diffusion is only small, λ₂Numerical values recited and Iterative Matrix eigenvalue u₁、u₂In close relations, as Iterative Matrix eigenvalue u₁=u₂Time, represent that diffusion is only small equally, but work as u₁≠u₂, it is meant that there is Characteristics of Mutation in prestack CRP road collection lineups, diffusion need to carry out along sudden change direction, eigenvalue λ₂More big, it is meant that the enhancing of diffusion；Characteristic vector v₁、v₂It is used for controlling dispersal direction, by the diffusion strengthening diffusion along collection lineups direction, prestack CRP road, reduce on vertical lineups direction, reaches the purpose strengthening prestack CRP road collection lineups.

(4), after obtaining diffusion tensor parameter, it is achieved solving of Anisotropic Nonlinear equation, diffusion tensor equation is as follows:

$\frac{\∂u}{\∂t}=\mathrm{div}(D\·\▿u^k)---\left(20\right)$

The left side for diffusion tensor equation (20) solves, it is possible to replace differential by difference, obtains the discrete form of diffusion tensor equation:

$\frac{\∂u}{\∂t}=\frac{u^{k+1}-u^k}{\mathrm{\Δt}}---\left(21\right)$

U in formula^kAnd u^k+1Respectively section result after k time and k+1 iteration, Δ t is step-length.

Diffusion tensor equation (20) the right has been solved:

$\begin{array}{c}\mathrm{div}(D\·{\▿u}^k)=\left(\begin{array}{cc}{\∂}_x& {\∂}_y\end{array}\right)\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)\left(\begin{array}{c}{\∂}_x{u}^k\\ {\∂}_y{u}^k\end{array}\right)\\ ={\∂}_x\left({a\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left({b\∂}_x{u}^k\right)+{\∂}_y\left({c\∂}_y{u}^k\right)\end{array}---\left(22\right)$

Diffusion tensor equation (20) the right and left simultaneous is got up for:

$u^{k+1}=u^k+\mathrm{\Δt}({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left(b{\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right))---\left(23\right)$

Wherein, and Δ t ∈ (0,0.3], achieved that by formula (23) and each layer section after decomposing is diffused Filtering Processing, the random noise that compacting prestack CRP road is concentrated strengthens the seriality of lineups.

(5) build root-mean-square mutation operator Θ, root-mean-square mutation operator Θ according to fabric tensor T and diffusion tensor D and can strengthen the abrupt information that lineups seriality can keep again prestack CRP road to concentrate, strengthen the identification ability to abrupt information:

$\mathrm{\Θ}=\sqrt{\frac{\mathrm{\Σ}\left(\mathrm{diag}\left(T^{0}T\right)\right)}{\mathrm{\Σ}\left(\mathrm{diag}\left(T^0\right)\right)\·\mathrm{\Σ}\left(\mathrm{diag}\left(T\right)\right)}}---\left(24\right)$

Wherein, diag () is matrix T⁰The leading diagonal of T, T⁰Initial construction tensor, T is current iteration fabric tensor, Θ ∈ [0,1], and near sudden change, Θ levels off to 0；During away from sudden change, Θ levels off to 1.

Being added in diffusion tensor equation (20) by root-mean-square mutation operator Θ, obtaining Anisotropic Nonlinear equation is:

$u^{k+1}=u^k+\mathrm{\Δt\Θ}\·({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right)).---\left(25\right)$

Above-mentioned steps 3) in, diffusion coefficient method comprises the following steps:

(1) diffusion coefficient is set: for making diffusion process that parameter to arrange more robustness, using Biweight(than watt) norm is as diffusion coefficient

In formula, k is diffusion threshold value, and when section neighborhood gradient is less than diffusion thresholding k, current pixel and its neighborhood territory pixel belong to same smooth domain, then neighborhood territory pixel is bigger to the contribution that current pixel is smooth；When neighborhood gradient is more than diffusion thresholding k and less than thresholdingTime, the probability that current pixel and its neighborhood territory pixel belong to same smooth domain is less, then the contribution that current pixel is smooth is reduced by neighborhood territory pixel；When neighborhood gradient is more than diffusion thresholdingTime, current pixel and its neighborhood territory pixel are not belonging to same smooth domain, then current pixel is smooth not by the impact of its neighborhood territory pixel.

(2) (MAD) method solves diffusion threshold value to adopt neighborhood inside gradient absolute deviation to estimate according to diffusion coefficient: in order to control to carry out, on different resolution level image, the filtering that yardstick is different, adopting MAD method, diffusion thresholding k (l) is relevant to Decomposition order l:

$k\left(l\right)=\frac{1}{0.6745}\·\mathrm{median}\left(\right||\▿u\left(l\right)-\mathrm{median}\left(\right||\▿u\left(l\right)|\left|\right)|\left|\right)\sqrt{2\log (l+1)/l}---\left(27\right)$

In above-mentioned steps 4) in, iterations quantitative analysis includes following three kinds of methods:

(1) mean square deviation (MSE) method:

$\mathrm{MSE}=\frac{1}{M*N}\underset{(i,j)=1}{\overset{M,N}{\mathrm{\Σ}}}{(u*(i,j)-u(i,j))}^2---\left(28\right)$

In formula, MSE is mean square deviation, u^*Filtered image and original image is represented respectively with u.MSE value is more little, and noise reduction is more good.

(2) signal to noise ratio (snr) method:

$\mathrm{SNR}\left(\mathrm{dB}\right)=10{\log }_{10}\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}u{(i,j)}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(u^{*}(i,j)-u(i,j))}^2}---\left(29\right)$

u^*Represent filtered image and original image respectively with u, after process, the more big explanation treatment effect of SNR is more good.

(3) Y-PSNR (PSNR) method:

$\mathrm{PSNR}=10{\log }_{10}\left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(30\right)$

Wherein, MSE is mean square deviation, and Y-PSNR PSNR is more big represents that the effect of filtering is more obvious.

According to above-mentioned three kinds of quantitative analysis standards, it is preferable that go out three and reach the iterations that statistics is optimum, what make to reach smooth lineups is maintained with edge detail information, keeps good fidelity.

As shown in Figure 2, wherein Fig. 2 (a) is original road collection, Fig. 2 (b) is the prestack CRP road collection after adopting Focus denoising module to filter random noise, the Fig. 2 (c) noise result for leaching, it follows that Focus denoising module is not satisfactory to the filtration result of random noise from figure, from the noise section filtered, the method compromises effective reflected signal of target zone and the lineups information of superficial part significantly.

As it is shown on figure 3, wherein Fig. 3 (a) is original road collection, Fig. 3 (b) is the prestack CRP road collection after adopting the present invention to filter random noise, and Fig. 3 (c) is the noise result leached；It follows that the present invention has not only filtered random noise preferably, protected effective lineups information from figure; and the upper and lower lineups feature of target zone is clear; lateral continuity is better improved; lineups have have by force weak; wave group feature becomes apparent from; the energy coincidence of far and near offset distance is strengthened, and the fidelity of amplitude is better.It will be seen that the present invention is less to the damage of effective information from the noise section filtered, overcome the prior art Focus module defect to the bigger damage of lineups.

As shown in Figure 4, wherein, a is well-log information, and b is the Zheng Yan road collection utilizing ray model to generate according to aboveground reflection coefficient, and c adopts Ya Zaohou road of the present invention collection, and d is road collection after well control amplitude recovery, and e is original road collection；As can be seen from the figure: the signal to noise ratio of original road collection e is relatively low, to the tracking of lineups difficulty, after well control amplitude recovery, collection d in road has recovered energy preferably, the method but the seriality of lineups, signal to noise ratio substantially diffusing filter pressure not as improving is made an uproar.Zheng Yan road collection b and employing Ya Zaohou road of the present invention collection not only have the good goodness of fit near reference lamina, and the lineups information of all the other degree of depth also has the good goodness of fit, and true subsurface reflector is effectively identified；But after well control amplitude recovery almost these effective lineups information of None-identified on road collection d and original road collection e.Therefore present invention improves section signal to noise ratio, and its resolution is true and reliable.

As it is shown in figure 5, wherein, A line is original road collection amplitude curve, and B line is ray model curve, and C line is road collection amplitude curve after well control amplitude recovery, and D line is that pressure of the present invention is made an uproar curve；By contrasting each bar curve it can be seen that ray model curve B is characterized by the increase along with offset distance, amplitude is gradually increased, and trend is comparatively mild；The trend of original road collection amplitude curve A and model are not inconsistent substantially, and near migration range place amplitude is relatively stable, along with the increase of offset distance has a decline by a relatively large margin, has a degree of rise afterwards, and remote offset distance also has fluctuation；Well control amplitude recovery curve C is close at variation tendency and the original amplitude curve of near migration range, in, the variation tendency of remote offset distance consistent with model curve, but locally deviation oscillation vibrates bigger；The pressure of the present invention curve D that makes an uproar is higher near migration range and the ray model curve B goodness of fit, it is slightly less than ray model curve B at middle offset distance, at remote offset distance higher than ray model curve B, but overall trend trend is substantially identical with ray model curve B, similarity coefficient is higher, and slickness is better, vibrates less.Each amplitude curve after being processed by contrast, it can be seen that make an uproar curve and the ray model curve B goodness of fit of pressure of the present invention is significantly high, has and protects width characteristic preferably, is conducive to extracting AVO attribute and the inverting that prestack CRP road is concentrated.

As shown in Figure 6, wherein, A line is chosen respectively: 0-15 °, B line: 15 °-30 °, C line: 30 ° of-48 ° CRP road collection is overlapped, and Fig. 6 (a) is stacked section schematic diagram before denoising, Fig. 6 (b) is stacked section schematic diagram after employing denoising of the present invention；Find by comparing stack result before and after denoising, before denoising, remote, in, the energy of near migration range, phase place concordance in the horizontal poor, there is bigger difference in waveform；After adopting denoising of the present invention, remote, in, the seriality of nearly deflection energy, phase place be improved, the useful signal energy of different signal to noise ratios is effectively recovered, and after denoising, the fidelity of amplitude is better.

As shown in Figure 7, Figure 8, comparison diagram 7 and Fig. 8 are known, based on the VP/VS inverting section that the integrated process of the present invention obtains, the sand shale lateral continuity of its reservoir, definition are higher than adopting prior art denoising inverting section, in conjunction with Well logs correlation, the lineups information that the present invention is finally inversed by after processing is authentic and valid, and sand shale distribution characteristics also improves simultaneously.

As shown in Figure 9, wherein Fig. 9 (a) adopts the inverting section that old process obtains, show that 6 wells are probably mud stone in certain time depth, 5 wells are probably sandstone, and the upper result shown of inverting section is just the opposite after Fig. 9 (b) adopts denoising of the present invention, the probability that 6 wells (solid circles) destination layer is sandstone is bigger, and the probability that 5 wells (dotted line circle) are sandstone is less.Real well data shows, 6 well target zones have oil/gas show, and 5 well target zones are then without oil/gas show.It can be said that the resolution of pre-stack data, signal to noise ratio and fidelity are all achieved good effect by the bright employing present invention.

The various embodiments described above are merely to illustrate the present invention; wherein the structure of each parts, connected mode and processing technology etc. all can be varied from; every equivalents carried out on the basis of technical solution of the present invention and improvement, all should not get rid of outside protection scope of the present invention.

Claims (5)

1., based on a multiple dimensioned anisotropic diffusion filtering method for prestack CRP road collection, it comprises the following steps:

1) the prestack CRP road collection of input is carried out Regularization, utilize the gaussian kernel function G of different variances sigma value_σWith seismic profile convolution；

2) adopt two dimension Mallat algorithm that the seismic profile after Regularization is done multi-resolution decomposition, using symmetrical wavelet sym6 as mother wavelet, original image is carried out N shell wavelet decomposition；And every sub-section after decomposing is initialized；

3) each sub-section of multi-resolution decomposition after initialization is tried to achieve diffusion threshold value through diffusion coefficient method, enter step 5)；Meanwhile, after being undertaken each sub-section of multi-resolution decomposition after initializing trying to achieve diffusion tensor parameter based on the anisotropy parameter of diffusion tensor method, substitution Anisotropic Nonlinear equation carries out anisotropic diffusion filtering, and enters step 4)；

4) the sub-section that iteration each time is obtained, calculates signal to noise ratio, mean square deviation and Y-PSNR, finally determines best iterations, obtains the sub-section of earthquake of the best according to this iterations；

5) adopt two dimension Mallat restructing algorithm that each sub-profile information after iteration carries out lossless reconstruction prestack CRP road collection, after reconstruct, export the prestack CRP road collection of the best；Wherein best prestack CRP road collection is:

$C_j=H_r^{*}{H}_c^{*}{C}_{j+1}+H_r^{*}{H}_c^{*}{D}_{j+1}^{*}+H_r^{*}{H}_c^{*}{D}_{j+1}^2+H_r^{*}{H}_c^{*}{D}_{j+1}^3,$

In formula, C_j+1For low frequency component,For horizontal high frequency component,For vertical high frequency component,For diagonal angle high fdrequency component,WithRespectively representative with scaling filter coefficient to arrayThe conjugation of operator of row and column effect.

2. the multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection as claimed in claim 1, it is characterised in that: described step 2) in, described multi-resolution decomposition comprises the following steps:

(1) setFor picture signal to be analyzed, its two dimension is approached image and is:

$A_{j}f=A_{j+1}f+D_{j+1}^{1}f+D_{j+1}^{2}f+D_{j+1}^{3}f,---\left(1\right)$

Wherein,

Wherein, A_jFor the result of the jth layer after decomposing, C_j+1(m, n) for low frequency component,For horizontal high frequency component,For vertical high frequency component,For diagonal angle high fdrequency component,For wavelet basis function, m is the translational movement n of wave filter in column direction is the translational movement by line direction wave filter；

Utilize the orthogonality of scaling function and wavelet function, formula (1) and (2) obtain:

$C_{j+1}(m,n)=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)h(l-2n)C_j(k,l)$

$\left\{\begin{array}{c}{D}_{j+1}^1=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)g(l-2n)C_j(k,l)\\ D_{j+1}^2=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)h(l-2n)C_j(k,l)\\ D_{j+1}^3=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)g(l-2n)C_j(k,l)\end{array}\right.$

To low frequency component C_j+1, horizontal high frequency componentVertical high frequency componentWith diagonal angle high fdrequency componentCarry out two dimension Mallat to decompose:

$\left\{\begin{array}{c}{C}_{j+1}=H_r{H}_c{C}_j\\ D_{j+1}^1=H_r{G}_c{C}_j\\ D_{j+1}^2=G_r{H}_c{C}_j\\ D_{j+1}^3=G_r{G}_c{G}_j\end{array}\right.,j=\mathrm{0,1},...,J$

In formula, H_rAnd H_cRespectively expression with scaling filter coefficient to arrayThe operator of row and column effect, G_rAnd G_cThe expression wavelet filter coefficient operator to row and column effect respectively；

(2) effective information of the seismic section image after decomposing is concentrated mainly on low frequency component C_j+1In, directly to low frequency component C_j+1Carry out anisotropic diffusion filtering process；For the high fdrequency component of seismic section image, it is carried out threshold value shrink process；Threshold value chooses formula:

$\mathrm{\η}=\mathrm{\δ}\·2\sqrt{\log N}/\sqrt{\log (l+1)},\mathrm{\δ}=\mathrm{median}\left(u\right)/0.6745$

Wherein, η is the wavelet threshold set, and the standard deviation that δ is high fdrequency component is estimated, u is image wavelet component, and median represents and takes mediant, and l is Decomposition order, and N is pixel number in component.

3. the multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection as claimed in claim 1 or 2, it is characterised in that: in described step 3), described diffusion tensor method comprises the following steps:

(1) set diffusing filter equation as:

$\frac{\∂u(x,y,t)}{\∂t}=\▿\·(D\left(\right|{\▿u}_{\mathrm{\σ}}\left|\right)\▿u)$

Wherein,For diffusion tensor, σ is scaling function,For rightCarry out Regularization, wherein

${\▿u}_{\mathrm{\σ}}=\▿(K_{\mathrm{\σ}}*u(x,y,t)),K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\σ}}^2})$

In formula, K_σRepresent that yardstick is the Gaussian function of σ；

(2) to each sub-section of multi-resolution decomposition after initializing along the derivation of room and time direction, it is determined that fabric tensor, and the eigenvalue of fabric tensor is solved and unit character is vectorial, to determine seismic profile local image structure:

1. fabric tensor T is:

$T(\▿u)=\▿u_{\mathrm{\σ}}\⊗{\▿u}_{\mathrm{\σ}}$

Assume fabric tensor $T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right),$ T is the symmetric positive definite matrix of 2 × 2；

2. adopting the Gauss-Seidel method with accelerated factor that fabric tensor T is solved, the eigenvalue obtaining fabric tensor T matrix is:

$u_1=\frac{1}{2}(T_{11}+T_{12}+\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

$u_2=\frac{1}{2}(T_{11}+T_{12}-\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

Unit character vector v₁, v₂It is parallel to gradient direction, and unit character vector v 1 meets

$\left(\begin{array}{c}\cos \mathrm{\α}\\ \sin \mathrm{\α}\end{array}\right)=v_1\left|\right|\left(\begin{array}{c}2T_{12}\\ T_{22}-T_{11}+\sqrt{{(T_{11}-T_{22})}^2+{{4T}_{12}}^2}\end{array}\right);$

(3) according to the characteristic vector of the fabric tensor T matrix tried to achieve in step (2) build diffusion tensor D, diffusion tensor D be also symmetric positive definite matrix and its characteristic vector identical with the characteristic vector of fabric tensor T, the eigenvalue of diffusion tensor D is λ₁And λ₂, solve the parameter in diffusion tensor D matrix: diffusion tensor D is:

$D=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right);$

(4), after obtaining diffusion tensor parameter a, b, c, d, it is achieved solving of Anisotropic Nonlinear equation, diffusion tensor equation is as follows:

$\frac{\∂u}{\∂t}=\mathrm{div}(D\·{\▿u}^k)---\left(3\right)$

The left side for diffusion tensor equation (3) solves, and replaces differential to obtain the discrete form of diffusion tensor equation by difference:

$\frac{\∂u}{\∂t}=\frac{u^{k+1}-u^k}{\mathrm{\Δt}}$

U in formula^kAnd u^k+1Respectively section result after k time and k+1 iteration, Δ t is step-length；

Diffusion tensor equation (3) the right has been solved:

$\begin{array}{c}\mathrm{div}(D\·{\▿u}^k)=\left(\begin{array}{cc}{\∂}_x& {\∂}_y\end{array}\right)\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)\left(\begin{array}{c}{\∂}_x{u}^k\\ {\∂}_y{u}^k\end{array}\right)\\ ={\∂}_x\left({a\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left({b\∂}_x{u}^k\right)+{\∂}_y\left({c\∂}_y{u}^k\right)\end{array}$

Diffusion tensor equation (1) the right and left simultaneous is got up for:

$u^{k+1}=u^k+\mathrm{\Δt}({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left(b{\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right))---\left(4\right)$

Wherein, and Δ t ∈ (0,0.3], realize each layer section after decomposing is diffused Filtering Processing by formula (4), the random noise that compacting prestack CRP road is concentrated strengthens the seriality of lineups；

(5) root-mean-square mutation operator Θ is built according to fabric tensor T and diffusion tensor D:

$\mathrm{\Θ}=\sqrt{\frac{\mathrm{\Σ}\left(\mathrm{diag}\left(T^{0}T\right)\right)}{\mathrm{\Σ}\left(\mathrm{diag}\left(T^0\right)\right)\·\mathrm{\Σ}\left(\mathrm{diag}\left(T\right)\right)}};$

Wherein, diag () is matrix T⁰The leading diagonal of T, T⁰Initial construction tensor, T is current iteration fabric tensor, Θ ∈ [0,1], and near sudden change, Θ levels off to 0；During away from sudden change, Θ levels off to 1；

Being added in diffusion tensor equation (3) by root-mean-square mutation operator Θ, obtaining Anisotropic Nonlinear equation is:

$u^{k+1}=u^k+\mathrm{\Δt\Θ}\·({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right)).$

4. the multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection as claimed in claim 3, it is characterised in that: in described step (3), four parameter a in described diffusion tensor D matrix, b, c and d solve as follows:

1. the eigenvalue λ of diffusion tensor D is solved₁And λ₂, require to adopt sudden change enhancing coefficient method and the relevant two kinds of methods of coefficient method that strengthen to estimate eigenvalue λ according to actual seismic profile₁And λ₂:

Sudden change strengthens coefficient method: take eigenvalue λ₂=1；Characteristic vector v₁Characteristic of correspondence value λ₁For fabric tensor T eigenvalue u₁Monotonous descending function, and

λ₂=1

Relevant enhancing coefficient method: make eigenvalue λ₁And λ₂Value is as follows:

λ₁=c₁

Wherein,For the class natural impedance factor；

2. eigenvalue λ₁And λ₂Solve the parameter a of diffusion tensor D, b, c and d:

$D={(v_1,v_2)}^T\·\left(\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right)\·(v_1,v_2)=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)$

$\left\{\begin{array}{c}a={\mathrm{\λ}}_1{\cos }^{2}a+{\mathrm{\λ}}_2{\sin }^{2}a\\ b=({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)\cos a\sin a\\ c={\mathrm{\λ}}_1{\sin }^{2}a+{\mathrm{\λ}}_2{\cos }^{2}a\end{array}\right..$

5. the multiple dimensioned anisotropic diffusion filtering method based on prestack CRP road collection as claimed in claim 1 or 2, it is characterised in that: in described step 3), described diffusion coefficient method comprises the following steps:

(1) diffusion coefficient is set: using Biweight norm as diffusion coefficient:

In formula, k is diffusion threshold value；

(2) adopting neighborhood inside gradient absolute deviation method of estimation to solve diffusion threshold value according to diffusion coefficient, diffusion thresholding k (l) is relevant to Decomposition order l:

$k\left(l\right)=\frac{1}{0.6745}\·\mathrm{median}\left(\right||\▿u\left(l\right)-\mathrm{median}\left(\right||\▿u\left(l\right)|\left|\right)|\left|\right)\sqrt{2\log (l+1)/l}.$