CN103364835A - Stratum structure self-adaption median filtering method - Google Patents

Abstract

The invention discloses a stratum structure self-adaption median filtering method. The method comprises the steps of: 1, calculating a gradient structure tensor; 2, calculating a stratum lateral incontinuity measurement; 3, structuring a structure self-adaption median filtering device; and 4, performing structure self-adaption median filtering and processing. By utilizing the method, earthquake random noise and part coherent noise can be effectively attenuated, and stratum edge and detailed structure characteristics such as effective signals, geologic faults and amplitude abnormalities can be maintained in a filtering process; the algorithm of the technical scheme is easy to realize and has good operability; and moreover, a selection problem of a median filtering processing window size is avoided, the method meets the requirements of the median filtering processing of complicated signals, for data areas with good event extension, the filtering performance of the median filtering device is improved through adaptively controlling the direction characteristics of a filtering window, and for areas with geologic faults and angular unconformity, the maintaining performance of the median filtering device to stratum edge structures such as the geologic faults is improved through adaptively adjusting the size of the filtering window.

Description

A kind of stratal configuration adaptive median filter method

Technical field

The invention belongs to the seismic exploration technique field, relate to a kind of earthquake random noise attenuation method, especially a kind of stratal configuration adaptive median filter decay earthquake random noise of utilizing.

Background technology

In petroleum natural gas exploration, the deposition characteristicses such as the architectural feature such as tomography, passage and crack and river channel sand are the bases of finding and describe oil reservoir, differentiate and analyze that these architectural features have great importance in the seismic data.Yet, because the complicacy of these areal geology forms so that these zones are subject to the interference of noise in the geological data, adopt manual interpretation still the automatic detection algorithm such as relevant and rim detection describe and portray these zones comparatively difficulty that all becomes.Traditional seismic data denoising method, such as F-K territory tendency filtering, f-x territory predictive filtering and multichannel coherence filtering etc., mostly be the horizontal consistance of utilizing the seismic reflection lineups, easily cause the compacting to inclination lineups amplitude, also the discontinuous constructions such as little tomography and crack can be blured, even the incorrect link of large tomography both sides lineups can be caused.How can be preferably evenly the relation of seismic noise decay and useful signal protection received widely concern ^[1-6], and the poststack noise attenuation technique that can protect and strengthen the edge structures such as fracture especially is subject to seismic interpretation worker's attention ^[7-9]

Luo etc. ^[10]Thought based on order statistics provides a kind of geological data guarantor limit filtering (EPS, Edge-Preserving Smoothing) algorithm, the method adopts many windows of Kuwahara analytical technology, calculate present analysis point on every side average and the variance of the interior data of each sub-window, the average of the sub-window that minimum variance is corresponding is as the filtering output of current point.AlBinHassan etc. ^[11]Developed this guarantor limit filter thought, provided the sub-window subdivision method under the three-dimensional situation, the guarantor limit filtering that the EPS algorithm is used for 3D seismic data is processed.Lu etc. ^[12]Promoted the EPS algorithm of one dimension based on the thought of fitting of a polynomial, the method comprises in the sub-window of a series of translations of present analysis point signal the fitting of a polynomial error of given order is selected the Signal estimation value of sub-window corresponding to minimum error of fitting as the filtering output of current point by calculating.Yang Peijie etc. ^[13]Also analyze thought based on this many windows and provide a kind of directivity Edge keeping tomography enhancing technology, so that the faultage image in the coherent body obtains more clear demonstration.Liu etc. ^[14,15]Two-dimentional multistage median filtering device during application image is processed is eliminated the earthquake random noise, and has analyzed the impact of wave filter scale parameter selection on noise attentuation and signal detail protection.

Fehmers etc. ^[16]Introduce the guarantor limit filtering processing that differential equation filter method is used for seismic data, the method adopts the gradient-structure tensor to estimate the stratum direction, and the continuity on tolerance stratum retrains filtering direction and the filtering degree of diffusing filter equation with this.Lavialle etc. ^[17]Keep diffusion filter for the tomography of protection and enhancing fault tectonic specially in conjunction with line style and face type diffusing filter model.Sun Xi equality ^[18]With Wang Xusong etc. ^[19]By in diffusion equation, introducing Derivative Terms, Zhang Erhua etc. ^[20]Further improve the edge-protected performance of nonlinear anisotropic diffusion wave filter by the guarantor limit bound term in the modified diffusion equation.Yet, based on the filtering method of anisotropic diffusion equation with original earthquake data as starting condition, pass the formation filter action by the iteration of diffusion time, because the bad diffusion time that estimates acquisition optimum filtering effect, causing can not its filtering condition of Accurate Prediction ^[21]

Classical two dimension median filter method has excellent noise attentuation performance, is particularly useful for eliminating the peak noise in the non-stationary signal ^[22]Yet because its filter window keeps constant in filtering, and filtering operation do not have directivity characteristics, so that adopt median filtering method to process when having the seismic signal of layer structure, will inevitably cause the fuzzy of the damage of useful signal and discontinuous construction feature.People make great efforts to seek to have concurrently the wave filter of details protection and noise suppression feature always; wherein; the multistage median filtering device is to protect one of median filter from a kind of the most representative details that image processing field proposes, and the basic structure of matching image details---the sub-window of filtering is distinguished signal structure and noise effectively preferably by one group.

If u (x, t) is the discrete two-dimensional signal, be positioned at the box filter window of (2N+1) * (2N+1) of (x, t) for a center, can define a kind of basic subwindow and be

$\left\{\begin{array}{cc}{W}_1(x,t)=u(x+k,t)& -N\≤k\≤N\\ W_2(x,t)=u(x,t+k)& -N\≤k\≤N\\ W_3(x,t)=u(x+k,t+k)& -N\≤k\≤N\\ W_4(x,t)=u(x+k,t-k)& -N\≤k\≤N\end{array}\right.---\left(1\right)$

Then the output of multistage median filtering is defined as

y(x,t)=median{y _max(x,t),y _min(x,t),u(x,t)} （2）

Wherein,

$y_{\max }(x,t)=\underset{1\≤i\≤4}{\max }\left\{z_i(x,t)\right\},$

$y_{\min }(x,t)=\underset{1\&le;i\&le;4}{\mathrm{<figure-callout\; id="1"\; label="min"\; filenames="HDA00003443981600012.png,HDA00003443981600071.png"\; state="\{\{state\}\}">min</figure-callout>}}\left\{z_i(x,t)\right\},$

z _i(x,t)=median{u(x,t)∈W _i(x,t)},i=1,2,3,4，

Here, median{} represents conventional medium filtering operation.

Above prior art has following shortcoming:

(1) too much when the sub-window size Selection of two-dimentional multistage median filtering device, can't mate the signal structure feature in the pending data neighborhood, can cause obvious compacting to useful signal;

(2) get too smallly when the sub-window size Selection of two-dimentional multistage median filtering device, can't mate the signal structure feature in the pending data neighborhood, can cause the noise attentuation ability significantly to reduce.

Summary of the invention

The object of the invention is to overcome the shortcoming of above-mentioned prior art, provide a kind of stratal configuration adaptive median filter method, it can effective attenuation earthquake random noise and partial coherence noise, and utilize filtering size and the directivity characteristics of the confidence measure index control filter window of two kinds of stratal configurations, kept the filtering performance of median filter to lineups extension good data zone, improved the hold facility of filtering to the edges, stratum such as useful signal and tomography, amplitude anomaly and detailed structure feature, and be easy to realize, workable.

The objective of the invention is to solve by the following technical programs:

Stratal configuration adaptive median filter method of the present invention may further comprise the steps:

1) compute gradient structure tensor

At first utilize the gradient vector ▽ u of finite difference method Two-dimensional geological data u (x, t), with gradient vector ▽ u and its transposed vector (▽ u) ^TMultiplying each other obtains initial tensor matrix, then each component of initial tensor matrix is done the Gassian low-pass filter processing of yardstick ρ, thereby is obtained gradient-structure tensor S corresponding to two-dimension earthquake data u (x, t) _ρ(▽ u):

2) calculate the horizontal uncontinuity tolerance in stratum

At first to gradient-structure tensor S _ρ(▽ u) carries out matrix character and decomposes:

$S_{\mathrm{\&rho;}}(\&dtri;u)=\left(\begin{array}{cc}{v}_1& v_2\end{array}\right)\left(\begin{array}{cc}{\mathrm{\&mu;}}_1& 0\\ 0& {\mathrm{\&mu;}}_2\end{array}\right){\left(\begin{array}{cc}{v}_1& v_2\end{array}\right)}^T.---\left(4\right)$

In the formula, μ ₁And μ ₂Be gradient-structure tensor S _ρTwo eigenwerts of (▽ u), and μ ₁〉=μ ₂〉=0; v ₁And v ₂Be two proper vectors of gradient-structure tensor, and v ₁Change maximum direction corresponding to local amplitude, i.e. Signal gradient direction, v ₂Change minimum direction, the i.e. orientation of reflection line-ups corresponding to local amplitude;

Then utilize gradient-structure tensor S _ρTwo eigenwert μ of (▽ u) ₁And μ ₂Calculate the confidence measure of line style signal structure:

$\mathrm{CL}=\frac{{\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_2}{{\mathrm{\&mu;}}_1+{\mathrm{\&mu;}}_2}.---\left(5\right)$

In the formula, CL is the linear structure confidence measure, value between interval [0,1];

At last in conjunction with linear structure confidence measure CL and gradient-structure tensor S _ρThe eigenwert μ of (▽ u) ₂Calculate the horizontal uncontinuity confidence measure of earthquake:

CI=(1-CL)μ ₂. （6）

In the formula, CI is horizontal uncontinuity confidence measure, (1-CL) deviates from degree, μ for the relative linear structure of signal partial structurtes feature ₂Be under the square error least meaning signal along the energy variation intensity of locally consistent direction;

3) structural texture adaptive median filter

Structure has the structure adaptive medium filtering window of oval filter window:

$M(x,y)=\{y\&Element;(\frac{{((y-x)\&CenterDot;n\left(x\right))}^2}{{\mathrm{\&sigma;}}_1^2\left(x\right)}+\frac{{((y-x)\&CenterDot;n_{\&perp;}\left(x\right))}^2}{{\mathrm{\&sigma;}}_2^2\left(x\right)}\&le;1)\},---\left(8\right)$

In the formula, M (x, y) is current Filtering position x place structure adaptive medium filtering window, and y is space, the time location coordinate of the sampled point that comprises of current filter window, is the inner product of vectors operator, and n (x) is the direction vector of dip direction, n _⊥(x) be the direction vector of stratum gradient direction, σ ₁(x) be the major axis of oval filter window, σ ₂(x) be the minor axis of oval filter window;

Utilize earthquake linear structure confidence measure CL and horizontal uncontinuity confidence measure CI, determine the scale parameter σ relevant with the directional selectivity of the filter scale of structure adaptive median filter and filtering operation ₁(x) and σ ₂(x):

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x\right)=R_{\max }\&CenterDot;g\left(\mathrm{CI}\left(x\right)\right)\\ {\mathrm{\&sigma;}}_2\left(x\right)=(1-\mathrm{CL}\left(x\right)){\mathrm{\&sigma;}}_1\left(x\right)\end{array}\right.---\left(9\right)$

In the formula, R _MaxBe the full-size parameter of oval spectral window, g () is the monotone decreasing function about CI (x), and limit its span for (0,1], it is exponential function that the present invention gets g ()

$g\left(\mathrm{CI}\left(x\right)\right)=\exp (-\frac{\mathrm{CI}\left(x\right)}{\mathrm{\&beta;}})---\left(10\right)$

In the formula, β is the threshold parameter for CI (x), the rate of decay of control characteristic function;

4) the structure adaptive medium filtering is processed

At first according to the distribution situation of transverse energy variation on time and space in the pending seismic section, seismic section Ω is divided into the some overlapping subregion Ω in certain border=∪ Ω that have _i, at every sub regions Ω _iIn get horizontal uncontinuity confidence measure CI (x) maximum value, namely

${\mathrm{CI}}_{\max ,i}=\max \left\{\mathrm{CI}\left(x\right){|}_{{x\&Element;\mathrm{\&Omega;}}_i}\right\},---\left(11\right)$

Then the threshold parameter β that obtains in the regional by following formula is

In the formula, α is that number percent is adjusted the factor, and thr is the ground noise threshold value;

Determine at last every sub regions Ω _iThe scale parameter σ of the medium filtering window that interior each space, time location x=(x, t) locate ₁(x) and σ ₂(x), thus carry out the operation of following structure adaptive medium filtering:

$\hat{u}\left(x\right)=\mathrm{median}\{u\left(y\right),y\&Element;M(x,y)\}.---\left(7\right)$

In the formula,

Be the Output rusults of the structure adaptive medium filtering at position x place, u (y) is each sampled point value in the medium filtering window M (x, y) at position x place.

Further, in the above step 1), gradient-structure tensor S _ρThe computing formula of (▽ u) is:

$S_{\mathrm{\&rho;}}(\&dtri;u)=G_{\mathrm{\&rho;}}\&CircleTimes;(\&dtri;u{(\&dtri;u)}^T)=\left(\begin{array}{cc}{G}_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;x}\right)}^2& G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;x}\frac{\&PartialD;u}{\&PartialD;t}\right)\\ G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;t}\frac{\&PartialD;u}{\&PartialD;x}\right)& G_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}^2\end{array}\right),---\left(3\right)$

In the formula, $\&dtri;u=\left[\begin{array}{cc}\frac{\&PartialD;u}{\&PartialD;x}& \frac{\&PartialD;u}{\&PartialD;t}\end{array}\right]$ Be gradient vector, G _ρFor having the dimensional Gaussian low-pass filter function G of scale parameter ρ _ρ(x, t)=exp ((x ²+ t ²)/2 ρ ²), Be the convolution operator, T is the matrix transpose operator.

The present invention has following beneficial effect:

The present invention can effective attenuation earthquake random noise and partial coherence noise, and in filtering can remain valid signal and edge, stratum and the detailed structure features such as tomography, amplitude anomaly; The algorithm content of this technical scheme is easy to realization, and operability is good; Avoided simultaneously the selection problem of medium filtering processing window size, be applicable to the demand that the sophisticated signal medium filtering is processed, to the good data area of lineups extension, the directivity characteristics of adaptive control filter window improves the filtering performance of median filter, to tomography and unconformability of dip zone, the size that self-adaptation is adjusted filter window improves median filter to the retention of the stratum marginal textures such as tomography.

Description of drawings

Fig. 1 is the schematic diagram of the present invention's structure adaptive median filter window of constructing;

Fig. 2 is the major axis control function figure of structure adaptive medium filtering window;

Fig. 3 is the synthetic seismic data demonstration figure that contains the bad ground structure of tomography;

Fig. 4 is the signal to noise ratio (S/N ratio) change curve of structure adaptive median-filtered result corresponding to different threshold parameters;

Fig. 5 is the window parameter demonstration figure for the structure adaptive median filter that adds the generated data of making an uproar;

Fig. 6 is that the filtering result who adds the generated data of making an uproar contrasts demonstration figure;

Fig. 7 is that actual seismic data filtering result contrasts demonstration figure;

Fig. 8 is as a result mean amplitude spectrum contrast of actual seismic data filtering demonstration figure;

Embodiment

Stratal configuration adaptive median filter method of the present invention may further comprise the steps:

1) compute gradient structure tensor

At first utilize the gradient vector ▽ u of finite difference method Two-dimensional geological data u (x, t), with gradient vector ▽ u and its transposed vector (▽ u) ^TMultiplying each other obtains initial tensor matrix, then each component of initial tensor matrix is done the Gassian low-pass filter processing of yardstick ρ, thereby is obtained gradient-structure tensor S corresponding to two-dimension earthquake data u (x, t) _ρ(▽ u): gradient-structure tensor S _ρThe computing formula of (▽ u) is:

$S_{\mathrm{\&rho;}}(\&dtri;u)=G_{\mathrm{\&rho;}}\&CircleTimes;(\&dtri;u{(\&dtri;u)}^T)=\left(\begin{array}{cc}{G}_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;x}\right)}^2& G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;x}\frac{\&PartialD;u}{\&PartialD;t}\right)\\ G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;t}\frac{\&PartialD;u}{\&PartialD;x}\right)& G_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}^2\end{array}\right),---\left(3\right)$

In the formula, $\&dtri;u=\left[\begin{array}{cc}\frac{\&PartialD;u}{\&PartialD;x}& \frac{\&PartialD;u}{\&PartialD;t}\end{array}\right]$ Be gradient vector, G _ρFor having the dimensional Gaussian low-pass filter function G of scale parameter ρ _ρ(x, t)=exp ((x ²+ t ²)/2 ρ ²),

Be the convolution operator, T is the matrix transpose operator.

2) calculate the horizontal uncontinuity tolerance in stratum

At first to gradient-structure tensor S _ρ(▽ u) carries out matrix character and decomposes:

$S_{\mathrm{\&rho;}}(\&dtri;u)=\left(\begin{array}{cc}{v}_1& {\mathrm{<figure-callout\; id="2"\; label="v"\; filenames="HDA00003443981600071.png"\; state="\{\{state\}\}">v</figure-callout>}}_2\end{array}\right)\left(\begin{array}{cc}{\mathrm{\&mu;}}_1& 0\\ 0& {\mathrm{\&mu;}}_2\end{array}\right){\left(\begin{array}{cc}{v}_1& v_2\end{array}\right)}^T.---\left(4\right)$

In the formula, μ ₁And μ ₂Be gradient-structure tensor S _ρTwo eigenwerts of (▽ u), and μ ₁〉=μ ₂〉=0; v ₁And v ₂Be two proper vectors of gradient-structure tensor, and v ₁Change maximum direction corresponding to local amplitude, i.e. Signal gradient direction, v ₂Change minimum direction, the i.e. orientation of reflection line-ups corresponding to local amplitude;

Then utilize gradient-structure tensor S _ρTwo eigenwert μ of (▽ u) ₁And μ ₂Calculate the confidence measure of line style signal structure:

$\mathrm{CL}=\frac{{\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_2}{{\mathrm{\&mu;}}_1+{\mathrm{\&mu;}}_2}.---\left(5\right)$

In the formula, CL is the linear structure confidence measure, value between interval [0,1];

At last in conjunction with linear structure confidence measure CL and gradient-structure tensor S _ρThe eigenwert μ of (▽ u) ₂Calculate the horizontal uncontinuity confidence measure of earthquake:

CI=(1-CL)μ ₂. （6）

In the formula, CI is horizontal uncontinuity confidence measure, (1-CL) deviates from degree, μ for the relative linear structure of signal partial structurtes feature ₂Be under the square error least meaning signal along the energy variation intensity of locally consistent direction;

3) structural texture adaptive median filter

Structure has the structure adaptive medium filtering window of oval filter window:

$M(x,y)=\{y\&Element;(\frac{{((y-x)\&CenterDot;n\left(x\right))}^2}{{\mathrm{\&sigma;}}_1^2\left(x\right)}+\frac{{((y-x)\&CenterDot;n_{\&perp;}\left(x\right))}^2}{{\mathrm{\&sigma;}}_2^2\left(x\right)}\&le;1)\},---\left(8\right)$

In the formula, M (x, y) is current Filtering position x place structure adaptive medium filtering window, and y is space, the time location coordinate of the sampled point that comprises of current filter window, is the inner product of vectors operator, and n (x) is the direction vector of dip direction, n _⊥(x) be the direction vector of stratum gradient direction, σ ₁(x) be the major axis of oval filter window, σ ₂(x) be the minor axis of oval filter window;

Utilize earthquake linear structure confidence measure CL and horizontal uncontinuity confidence measure CI, determine the scale parameter σ relevant with the directional selectivity of the filter scale of structure adaptive median filter and filtering operation ₁(x) and σ ₂(x):

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x\right)=R_{\max }\&CenterDot;g\left(\mathrm{CI}\left(x\right)\right)\\ {\mathrm{\&sigma;}}_2\left(x\right)=(1-\mathrm{CL}\left(x\right)){\mathrm{\&sigma;}}_1\left(x\right)\end{array}\right.---\left(9\right)$

In the formula, R _MaxBe the full-size parameter of oval spectral window, g () is the monotone decreasing function about CI (x), and limit its span for (0,1], it is exponential function that the present invention gets g ()

$g\left(\mathrm{CI}\left(x\right)\right)=\exp (-\frac{\mathrm{CI}\left(x\right)}{\mathrm{\&beta;}})---\left(10\right)$

In the formula, β is the threshold parameter for CI (x), the rate of decay of control characteristic function;

4) the structure adaptive medium filtering is processed

At first according to the distribution situation of transverse energy variation on time and space in the pending seismic section, seismic section Ω is divided into the some overlapping subregion Ω in certain border=∪ Ω that have _i, at every sub regions Ω _iIn get horizontal uncontinuity confidence measure CI (x) maximum value, namely

${\mathrm{CI}}_{\max ,i}=\max \left\{\mathrm{CI}\left(x\right){|}_{{x\&Element;\mathrm{\&Omega;}}_i}\right\},---\left(11\right)$

Then the threshold parameter β that obtains in the regional by following formula is

In the formula, α is that number percent is adjusted the factor, and thr is the ground noise threshold value;

Determine at last every sub regions Ω _iThe scale parameter σ of the medium filtering window that interior each space, time location x=(x, t) locate ₁(x) and σ ₂(x), thus carry out the operation of following structure adaptive medium filtering:

$\hat{u}\left(x\right)=\mathrm{median}\{u\left(y\right),y\&Element;M(x,y)\}.---\left(7\right)$

In the formula,

Be the Output rusults of the structure adaptive medium filtering at position x place, u (y) is each sampled point value in the medium filtering window M (x, y) at position x place.

Below in conjunction with accompanying drawing above method of the present invention is done and to be further explained in detail explanation:

Dip direction is estimated and horizontal uncontinuity tolerance

Have several different methods to can be used for the systematicness analysis of stratal configuration, such as waveform coherence analysis, gradient-structure tensor analysis and gray level co-occurrence matrixes analysis] etc., wherein the gradient-structure Tensor Method is a kind of simple and effective approach.For two-dimension earthquake section u (x, t), corresponding gradient-structure tensor is described by following formula

$S_{\mathrm{\&rho;}}(\&dtri;u)=G_{\mathrm{\&rho;}}\&CircleTimes;(\&dtri;u{(\&dtri;u)}^T)=\left(\begin{array}{cc}{G}_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;x}\right)}^2& G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;x}\frac{\&PartialD;u}{\&PartialD;t}\right)\\ G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;t}\frac{\&PartialD;u}{\&PartialD;x}\right)& G_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}^2\end{array}\right),---\left(3\right)$

In the formula, ▽ is gradient operator; G _ρExpression two-dimensional Gaussian function G _ρ(x, t)=exp ((x ²+ t ²)/2 ρ ²), ρ is scale parameter;

Be the convolution operator; T is the matrix transpose operator.

In the definition of gradient-structure tensor, the out to out of the measurable signal characteristic of gradient-structure tensor has been determined in the gauss low frequency filter effect that the tensor matrix that is formed by gradient vector ▽ u and transposed vector thereof and yardstick are ρ.Should symmetry positive semidefinite matrix S _ρ(▽ u) makes matrix character and decomposes

$S_{\mathrm{\&rho;}}(\&dtri;u)=\left(\begin{array}{cc}{v}_1& {\mathrm{<figure-callout\; id="2"\; label="v"\; filenames="HDA00003443981600071.png"\; state="\{\{state\}\}">v</figure-callout>}}_2\end{array}\right)\left(\begin{array}{cc}{\mathrm{\&mu;}}_1& 0\\ 0& {\mathrm{\&mu;}}_2\end{array}\right){\left(\begin{array}{cc}{v}_1& v_2\end{array}\right)}^T.---\left(4\right)$

In the formula (4), eigenwert μ ₁〉=μ ₂〉=0 has reacted the Gauss neighborhood G of signal at locus (x, t) _ρ(x, t) interior amplitude variations intensity along characteristic direction; Proper vector v ₁And v ₂Provide local orthogonal orientation, v ₁Change maximum direction corresponding to local amplitude, i.e. Signal gradient direction, and v ₂Change minimum direction, the i.e. orientation of reflection line-ups corresponding to local amplitude.According to the value condition of two eigenwerts of gradient-structure tensor, distinguishable different signal structure feature, Bakker(2002) the middle confidence measure of introducing the line style signal structure:

$\mathrm{CL}=\frac{{\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_2}{{\mathrm{\&mu;}}_1+{\mathrm{\&mu;}}_2}.---\left(5\right)$

This confidence measure is value between interval [0,1], weighs the similarity degree of local signal feature and linear structure, also characterizes the order of accuarcy that local orientation is estimated simultaneously.If CL → 1 represents that then the local signal feature is tending towards linear structure, the accuracy of the stratigraphic dip of estimation is higher; Otherwise if CL → 0 then characterizes the local signal structure and deviated from the line style hypothesis, this may be caused by horizontal uncontinuity such as interference noise and tomography, crack etc.Can identify preferably the reflection line-ups that the part is linear structure according to the linear structure confidence measure; But, because what this degree of confidence value was judged is the relative intensity of variation of signal energy, easily be subject to the interference of noise, and value is nonsensical in without the zone of clear signal structure, thereby is difficult to from jamming pattern, identify clearly the laterally discontinuous construction such as tomography, crack.

Based on above-mentioned analysis, the present invention provides a kind of degree of confidence value that contains signal transverse energy change intensity information on the basis of linear structure confidence measure, is used for the laterally discontinuous construction such as identification tomography, crack.

CI=(1-CL)μ ₂. （6）

The relative linear structure of (1-CL) expression signal partial structurtes feature deviates from degree, μ in the formula (6) ₂Be characterized in the given tolerance neighborhood under the square error least meaning signal along the energy variation intensity of locally consistent direction ^[30]The zone that is linear structure at reflection line-ups, because CL → 1, i.e. (1-CL) → 0, this moment is along the signal intensity less of dip direction, i.e. μ ₂→ 0, thus so that CI → 0.Containing the laterally zone of uncontinuity such as tomography, crack, because CL → 0, i.e. (1-CL) → 1; And at horizontal discontinuity zone, because signal structure has deviated from local linear structure hypothesis, the consistance of signal orientation is indefinite herein, thereby also relatively large along the signal intensity intensity of locally coherence direction, i.e. μ ₂＞＞0, thus so that CI＞＞0.Therefore, can identify well the laterally discontinuous construction such as tomography and crack according to horizontal uncontinuity confidence measure CI.Because the signal transverse energy changes a μ ₂Existence, be weaker than the situation of signal energy for noise energy, CI has good anti-noise ability; In addition, CI also can differentiate zone without the clear signal structure (this moment μ ₂→ 0, i.e. CI → 0), this can be avoided the generation of some filtering illusions to adopting CI to instruct filtering particularly important.

The structure adaptive median filter

For two-dimension earthquake section u (x), note x=(x, t) is space, the time location of sampled point, and the structure adaptive medium filtering at definition position x place is operating as:

$\hat{u}\left(x\right)=\mathrm{median}\{u\left(y\right),y\&Element;M(x,y)\}.---\left(7\right)$

In the formula (7),

The Output rusults of the structure adaptive medium filtering at expression x place, position; (see Fig. 1 a), y is space, the time location coordinate of the sampled point that comprises of current filter window to the structure adaptive filter window at M (x, y) expression x place, position.Structure adaptive spectral window M (x, y) is the function of Filtering position x, and its size and shape depends on that the x place analyzes the signal structure feature in the neighborhood:

$M(x,y)=\{y\&Element;(\frac{{((y-x)\&CenterDot;n\left(x\right))}^2}{{\mathrm{\&sigma;}}_1^2\left(x\right)}+\frac{{((y-x)\&CenterDot;n_{\&perp;}\left(x\right))}^2}{{\mathrm{\&sigma;}}_2^2\left(x\right)}\&le;1)\},---\left(8\right)$

Wherein, expression inner product of vectors operator, n (x) is the direction vector of dip direction, σ ₁(x) be the major axis of oval filter window, n _⊥(x) be the direction vector of stratum gradient direction, σ ₂(x) be the minor axis of oval filter window.

In the formula (8), the major axis σ of oval filter window ₁(x) and minor axis σ ₂(x) filter scale of structure adaptive median filter and the directional selectivity of filtering operation have jointly been determined.Therefore, the design of wave filter scale parameter is the key that affects structure adaptive median filter condition, requirement according to the structure adaptive filtering operation, this paper provides the following adaptively selected strategy of wave filter scale parameter in conjunction with linear structure confidence measure CL and horizontal uncontinuity confidence measure CI:

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x\right)=R_{\max }\&CenterDot;g\left(\mathrm{CI}\left(x\right)\right)\\ {\mathrm{\&sigma;}}_2\left(x\right)=(1-\mathrm{CL}\left(x\right)){\mathrm{\&sigma;}}_1\left(x\right)\end{array}\right.---\left(9\right)$

Wherein, R _MaxStipulated the full-size of oval filter window, g () is the monotone decreasing function about CI (x), its span be defined as (0,1], getting g () in the literary composition is exponential function

$g\left(\mathrm{CI}\left(x\right)\right)=\exp (-\frac{\mathrm{CI}\left(x\right)}{\mathrm{\&beta;}})---\left(10\right)$

Wherein, β is the threshold parameter for CI (x), the rate of decay of control characteristic function.The β value is less, and the decline of exponential function is faster; Otherwise the decay of exponential function slows down.

The structure adaptive median filter has flexibly filter scale regulatory function and set direction characteristic according to the partial structurtes feature of signal, controls the σ of maximum filter scale ₁(x) determined by horizontal uncontinuity confidence measure CI (x), and the σ of adjustment set direction characteristic ₂(x) determined by linear structure confidence measure CL (x).Horizontal discontinuity zone is owing to CI (x)＞＞0, so that the major axis σ of ellipse window function in tomography, crack etc. ₁(x) become short; And in other zone with CI (x) → 0, σ is arranged ₁(x) → R _MaxBe the linear structure zone at reflection line-ups, because CL (x) → 1, so that the minor axis σ of ellipse window function ₂(x) atrophy; And in other zones with CL (x) → 0, σ is arranged ₂(x) → σ ₁(x).Therefore, shown in Fig. 1 b, be the zone that line style is extended in the reflection line-ups part, filter window is long and narrow anisotropy window (position x along the lineups orientation stretching ₁The zone); At horizontal discontinuity zones such as tomographies, adaptively atrophy of filter window is that fenestella is with coupling discontinuous construction feature (position x herein ₂The zone); And in the noise zone without the clear signal structure, filter window extends and is large near-isotropic window (position x ₃The zone).

Association type (9) and (10) as seen, the filtering performance of structure adaptive median filter is subject to parameters R _MaxJointly retrain the full-size parameters R with β _MaxThe Lubricity of control wave filter, threshold parameter β then adjusts wave filter to the maintenance degree of the edge structures such as tomography.Consider that the present invention adopts the gradient-structure tensor to estimate inclination angle and the regular degree on stratum, the scale parameter ρ of gradient-structure tensor has determined its tolerance aperture to local signal structure.Therefore, in order to ensure credible based on the filtering of two kinds of confidence measure constraints, require the full-size R of filter window _MaxBe no more than 2 ρ, namely get R _Max≤ 2 ρ.For threshold parameter β, as seen from Figure 2, when the β value is less than normal, slightly large horizontal uncontinuity confidence measure CI (x) namely causes the size reduction of filter window, thereby the edge structures such as some weak tomographies, crack are effectively protected, but also can weaken wave filter simultaneously at the filter effect in horizontal uncontinuity zone; When the β value is bigger than normal, only have very large CI (x) just can cause the filter window atrophy to enough little, thereby the noise attentuation ability of wave filter is enhanced in this case, but can causes to a certain extent the fuzzy and compacting of weak marginal texture.

Protect limit condition self-adaptation control method

The marginal texture retention of structure adaptive median filter is controlled by threshold parameter β, and selecting suitable β is the key that guarantees can effectively keep in the filtering edges, stratum such as tomography, crack.Yet seismic data from the shallow-layer to the deep layer, even the dynamic range that the signal transverse energy changes in same interval fluctuation also very greatly, is difficult to select uniformly a global threshold parameter beta usually.The present invention provides a kind of adaptively selected strategy of threshold parameter of Region Segmentation, namely according to the distribution situation of transverse energy variation on time and space in the pending section, section Ω is divided into the some overlapping subregion Ω in certain border=∪ Ω that have _i, at every sub regions Ω _iIn get horizontal uncontinuity confidence measure CI (x) maximum value, namely

${\mathrm{CI}}_{\max ,i}=\max \left\{\mathrm{CI}\left(x\right){|}_{{x\&Element;\mathrm{\&Omega;}}_i}\right\},---\left(11\right)$

Thereby the adaptive threshold parameter beta that obtains in this zone by following formula is

Wherein, α is that number percent is adjusted the factor, determines according to the requirement overall situation of protecting edge filter; Thr is the ground noise threshold value, is determined by the overall noise interference level of section.

By pending section is divided into some subregions, to in the formula (10) the absolute value of threshold parameter β be converted into relative value to the number percent factor-alpha by formula (11) and (12), thereby the edge that can unify to adjust the structure adaptive median filter keeps condition.In real data is processed, do not need accurately to determine the size of subregion, a kind of typical selection is to comprise 100 roads in every sub regions, per pass comprises 150 sampled points, can satisfy the requirement that filtering keeps the edge features such as tomography, crack.

Numerical Simulation Results

The synthetic model data

Validity for the protection edge denoising of verifying the structure adaptive median filter; this paper chooses the Ricker wavelet that dominant frequency is 20Hz; adopt ray-tracing scheme to obtain a synthetic seismic data that contains the bad ground structural model of tomography and (see Fig. 3 a); this random noise of band limit for height that adds 5dB in generated data obtains the noisy data of Fig. 3 b.

Adopt the structure adaptive median filter to process and add the generated data of making an uproar, the signal to noise ratio (S/N ratio) of considering generated data is very low, in the calculating of gradient-structure tensor, the dimensional Gaussian low-pass filter function is got larger scale parameter ρ=4.0, and the maximum filtering size of structure adaptive median filter is taken as R _Max=ρ=4.0, namely corresponding out to out is 9 filter window.For the selected representative region of black rectangle frame among Fig. 3 b, analyze the threshold percentage factor-alpha to filtering result's impact.Adjust α since 10% with 10% interval variation to 200%, record each filtering result's signal to noise ratio (S/N ratio), obtain this regional filtering result's signal to noise ratio (S/N ratio) with the change curve of threshold percentage factor-alpha, as shown in Figure 4.As can be seen from the figure, the signal to noise ratio (S/N ratio) curve raises with the increase of α, and curve slows down gradually in the rear rising of α=90%.This is because along with the increase of α, guarantor's edge performance of wave filter weakens, as α〉after 90%, filtering has caused the damage of section part waveform, and offset to a certain extent filtering performance and strengthened the signal to noise ratio (S/N ratio) of bringing and increase.Thereby when processing this and add the generated data of making an uproar, select threshold percentage factor-alpha=90%, the major axis σ of the anisotropy window function of the structure adaptive median filter that this moment is corresponding ₁(x) and minor axis σ ₂(x) value distributes shown in Fig. 5 a and b.Comparison diagram 5a and Fig. 5 b are as seen, the scale parameter of structure adaptive median filter can mate the signal structure feature of generated data well, the major axis that is the oval filter window in lineups zone of linear structure in the part stretches, the minor axis atrophy, is formed with the anisotropy window of strong set direction characteristic; In fault region major axis and simultaneously atrophy of minor axis, form the small scale analysis window; And in the zone of no signal feature, major axis and minor axis all are stretched, and form the isotropy window of large scale.

The noise sections that adds the structure adaptive median-filtered result of the generated data of making an uproar and obtain is respectively shown in Fig. 6 a and b.In order to contrast the treatment effect of this paper method, use classical two dimension median filter device to process with the two-dimentional multistage median filtering device with signal detail protective capability and add the generated data of making an uproar.The window size of getting two kinds of median filters is identical with the maximum window size of structure adaptive median filter, be that the two dimension median filter device is selected 9 * 9 filter window, two dimension multistage median filtering device is got 9 filter length, and the filtering result of both correspondences and the noise sections that obtains are respectively shown in Fig. 6 c and e and d and f.

From Fig. 6 a as seen, behind the structure adaptive medium filtering, random noise is decayed largely, and the consistance of reflection line-ups is enhanced, and it is very clean that section integral body becomes; From Fig. 6 b as seen, only have the interference noise of stochastic distribution in the noise sections that the structure adaptive median filter obtains, without obvious useful signal structure, fault structure has obtained effective maintenance in filtering.The filtering result of two kinds of median filter method among comparison diagram 6c and the e, 9 * 9 two dimension median filter device has been eliminated most random noise, but also cause the damage of useful signal, from its noise sections Fig. 6 d that obtains as seen, filtering causes the fuzzy of fault information, also cause the damage of reflection line-ups, especially more obvious to the damage of inclined reflection structure; Though the two-dimentional multistage median filtering device of 9 filter lengths has stronger signal detail protective capability, but from its noise sections Fig. 6 f that obtains as seen, two-dimentional multistage median filtering device is limited to the random noise attenuation degree.

From above-mentioned analysis, relatively can find out: under the condition of getting identical filter window yardstick, compare two dimension median filter device and two-dimentional multistage median filtering device, the structure adaptive median filter can be according to the signal structure feature of seismic section, structure and the filter window that signal structure matches have been obtained good compromise between the edge structures such as Removing Random No and protection tomography adaptively; The noise attentuation performance of structure adaptive median filter is near 9 * 9 two dimension median filter device, and its protective capability to tomography and useful signal is better than the two-dimentional multistage median filtering device of 9 filter lengths.

The actual seismic data

Based on theoretical analysis result, the filtering that further the structure adaptive median filter is used for stacked seismic data is processed.Provide the poststack pure wave seismic data in certain oil field among Fig. 7 a, contain abundant discontinuous construction feature in the section, except two main tomographies, at central part the mixed and disorderly reflector space in one place is arranged, there is amplitude anomaly zone, a place in section on the lower side on the right side.Because be subject to random noise and gather the interference of footprint, the continuity of lineups is under some influence, be unfavorable for layer position a horizontal tracking, discontinuous construction automatic detection and follow-up processing and the explanation such as pick up.Use the structure adaptive median filter of this paper to process this section, select the subregion segmentation strategy of 100 roads, 150 sampled points of per pass; In the calculating of gradient-structure tensor, get the scale parameter ρ of Gaussian filter=3.0; The maximum filtering size of structure adaptive median filter is taken as R _Max=4.0, namely correspondence is 9 filter window size to the maximum; In order to reach preferably edge retention, get threshold percentage factor-alpha=50%, corresponding filtering result and the noise sections that obtains are respectively shown in Fig. 7 b and c.As a comparison, adopt equally two kinds of median filter method to process this section.Because it is too much that the filter window of two dimension median filter device is selected, serious to the useful signal damage, and select too smallly, can cause noise attentuation inadequate, so locating to get the window size of two dimension median filter device is 5 * 5, and two-dimentional multistage median filtering device adopts 9 filter lengths, and namely the maximum filtering size with the structure adaptive median filter is suitable, and the noise sections that corresponding filtering result and they obtain is respectively shown in Fig. 7 d and f, e and g.

Sectional view 7a and b before and after the comparison structure adaptive median filter can find out, structure adaptive medium filtering process has been eliminated most interference noise, has improved the continuity of reflection line-ups, and fault tectonic is enhanced, and the amplitude anomaly zone becomes clear.The noise sections that the analytical structure adaptive median filter obtains (such as Fig. 7 c) as seen, the interference noise of shallow-layer is very strong, the interference noise of mid-deep strata relatively a little less than, except random noise, also have the coherent noise of some wide-angles in the section; Original seismic section among the comparison diagram 7a be can't see obvious useful signal structure in the noise sections.The filtering result of two dimension median filter device among the analysis chart 7d, from seeing that visually noise is pressed substantially, signal to noise ratio (S/N ratio) is greatly improved, but the marginal texture of signal especially tomography, mixed and disorderly reflection etc. seriously blured, this point can be confirmed from its noise sections that obtains (such as Fig. 7 e).Can find out from Fig. 7 f and g, two-dimentional multistage median filtering device only has faint noise attentuation effect along the lineups trend, and filter capacity is very limited.Above result shows, the structure adaptive median filter all can effectively be suppressed the very noisy of shallow-layer and the small noise of deep layer, not only can random noise attenuation, also the partial coherence noise be can eliminate, and useful signal and edge, stratum and the detailed structure features such as tomography, amplitude anomaly kept well; The structure adaptive median filter not only has the noise attentuation ability of two dimension median filter device excellence, and can not cause because of the demand of protection edge, stratum and detailed structure again the obvious reduction of filter capacity as two-dimentional multistage median filtering device.

Result to original seismic section and three kinds of filtering methods is carried out spectrum analysis, and Fig. 8 provides the mean amplitude spectrum of corresponding section.Frequency spectrum discovery before and after contrast is processed, before the high band random noise amplitude after three kinds of methods are processed all will be lower than processing, and the main structure of frequency spectrum does not have destroyed; Yet structure adaptive median filter and two-dimentional multistage median filtering device have all kept the spectral amplitude of dominant frequency and useful information frequency range preferably, and the two dimension median filter device is comparatively serious to the damage of the spectral amplitude of dominant frequency and useful information frequency range.Above presentation of results structure adaptive median filter can effectively be suppressed the earthquake random noise, and damages hardly the frequency spectrum of useful information frequency range.

Claims (2)

1. a stratal configuration adaptive median filter method is characterized in that, may further comprise the steps:

1) compute gradient structure tensor

At first utilize the gradient vector ▽ u of finite difference method Two-dimensional geological data u (x, t), with gradient vector ▽ u and its transposed vector (▽ u) ^TMultiplying each other obtains initial tensor matrix, then each component of initial tensor matrix is done the Gassian low-pass filter processing of yardstick ρ, thereby is obtained gradient-structure tensor S corresponding to two-dimension earthquake data u (x, t) _ρ(▽ u);

2) calculate the horizontal uncontinuity tolerance in stratum

At first to gradient-structure tensor S _ρ(▽ u) carries out matrix character and decomposes:

$S_{\mathrm{\&rho;}}(\&dtri;u)=\left(\begin{array}{cc}{v}_1& v_2\end{array}\right)\left(\begin{array}{cc}{\mathrm{\&mu;}}_1& 0\\ 0& {\mathrm{\&mu;}}_2\end{array}\right){\left(\begin{array}{cc}{v}_1& v_2\end{array}\right)}^T.---\left(4\right)$

In the formula, μ ₁And μ ₂Be gradient-structure tensor S _ρTwo eigenwerts of (▽ u), and μ ₁〉=μ ₂〉=0; v ₁And v ₂Be two proper vectors of gradient-structure tensor, and v ₁Change maximum direction corresponding to local amplitude, i.e. Signal gradient direction, v ₂Change minimum direction, the i.e. orientation of reflection line-ups corresponding to local amplitude;

Then utilize gradient-structure tensor S _ρTwo eigenwert μ of (▽ u) ₁And μ ₂Calculate the confidence measure of line style signal structure:

$\mathrm{CL}=\frac{{\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_2}{{\mathrm{\&mu;}}_1+{\mathrm{\&mu;}}_2}.---\left(5\right)$

In the formula, CL is the linear structure confidence measure, value between interval [0,1];

At last in conjunction with linear structure confidence measure CL and gradient-structure tensor S _ρThe eigenwert μ of (▽ u) ₂Calculate the horizontal uncontinuity confidence measure of earthquake:

CI＝(1-CL)μ ₂. （6）

In the formula, CI is horizontal uncontinuity confidence measure, (1-CL) deviates from degree, μ for the relative linear structure of signal partial structurtes feature ₂Be under the square error least meaning signal along the energy variation intensity of locally consistent direction;

3) structural texture adaptive median filter

Structure has the structure adaptive medium filtering window of oval filter window:

$M(x,y)=\{y\&Element;(\frac{{((y-x)\&CenterDot;n\left(x\right))}^2}{{\mathrm{\&sigma;}}_1^2\left(x\right)}+\frac{{((y-x)\&CenterDot;n_{\&perp;}\left(x\right))}^2}{{\mathrm{\&sigma;}}_2^2\left(x\right)}\&le;1)\},---\left(8\right)$

In the formula, M (x, y) is current Filtering position x place structure adaptive medium filtering window, and y is space, the time location coordinate of the sampled point that comprises of current filter window, is the inner product of vectors operator, and n (x) is the direction vector of dip direction, n _⊥(x) be the direction vector of stratum gradient direction, σ ₁(x) be the major axis of oval filter window, σ ₂(x) be the minor axis of oval filter window;

Utilize earthquake linear structure confidence measure CL and horizontal uncontinuity confidence measure CI, determine the scale parameter σ relevant with the directional selectivity of the filter scale of structure adaptive median filter and filtering operation ₁(x) and σ ₂(x):

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x\right)=R_{\max }\&CenterDot;g\left(\mathrm{CI}\left(x\right)\right)\\ {\mathrm{\&sigma;}}_2\left(x\right)=(1-\mathrm{CL}\left(x\right)){\mathrm{\&sigma;}}_1\left(x\right)\end{array}\right.---\left(9\right)$

In the formula, R _MaxBe the full-size parameter of oval spectral window, g () is the monotone decreasing function about CI (x), and limit its span for (0,1], it is exponential function that the present invention gets g ()

$g\left(\mathrm{CI}\left(x\right)\right)=\exp (-\frac{\mathrm{CI}\left(x\right)}{\mathrm{\&beta;}})---\left(10\right)$

In the formula, β is the threshold parameter for CI (x), the rate of decay of control characteristic function;

4) the structure adaptive medium filtering is processed

At first according to the distribution situation of transverse energy variation on time and space in the pending seismic section, seismic section Ω is divided into the some overlapping subregion Ω in certain border=∪ Ω that have _i, at every sub regions Ω _iIn get horizontal uncontinuity confidence measure CI (x) maximum value, namely

${\mathrm{CI}}_{\max ,i}=\max \left\{\mathrm{CI}\left(x\right){|}_{{x\&Element;\mathrm{\&Omega;}}_i}\right\},---\left(11\right)$

Then the threshold parameter β that obtains in the regional by following formula is

In the formula, α is that number percent is adjusted the factor, and thr is the ground noise threshold value;

Determine at last every sub regions Ω _iThe scale parameter σ of the medium filtering window that interior each space, time location x=(x, t) locate ₁(x) and σ ₂(x), thus carry out the operation of following structure adaptive medium filtering:

$\hat{u}\left(x\right)=\mathrm{median}\{u\left(y\right),y\&Element;M(x,y)\}.---\left(7\right)$

In the formula,

Be the Output rusults of the structure adaptive medium filtering at position x place, u (y) is each sampled point value in the medium filtering window M (x, y) at position x place.

2. stratal configuration adaptive median filter method according to claim 1 is characterized in that, in the step 1), and gradient-structure tensor S _ρThe computing formula of (▽ u) is:

$S_{\mathrm{\&rho;}}(\&dtri;u)=G_{\mathrm{\&rho;}}\&CircleTimes;(\&dtri;u{(\&dtri;u)}^T)=\left(\begin{array}{cc}{G}_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;x}\right)}^2& G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;x}\frac{\&PartialD;u}{\&PartialD;t}\right)\\ G_{\mathrm{\&rho;}}\&CircleTimes;\left(\frac{\&PartialD;u}{\&PartialD;t}\frac{\&PartialD;u}{\&PartialD;x}\right)& G_{\mathrm{\&rho;}}\&CircleTimes;{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}^2\end{array}\right),---\left(3\right)$

In the formula, $\&dtri;u=\left[\begin{array}{cc}\frac{\&PartialD;u}{\&PartialD;x}& \frac{\&PartialD;u}{\&PartialD;t}\end{array}\right]$ Be gradient vector, G _ρFor having the dimensional Gaussian low-pass filter function G of scale parameter ρ _ρ(x, t)=exp ((x ²+ t ²)/2 ρ ²),

Be the convolution operator, T is the matrix transpose operator.

Images