[发明专利]一种地层结构自适应中值滤波方法

权利要求书

1.一种地层结构自适应中值滤波方法，其特征在于，包括以下步骤：

1）计算梯度结构张量

首先利用有限差分方法计算二维地震数据u(x,t)的梯度向量▽u，将梯度向量▽u与其转置向量(▽u)^T相乘得到初始张量矩阵，然后对初始张量矩阵的各个分量作尺度ρ的高斯低通滤波处理，从而得到二维地震数据u(x,t)对应的梯度结构张量S_ρ(▽u)；

2）计算地层横向不连续性度量

首先对梯度结构张量S_ρ(▽u)进行矩阵特征分解：

Sρ(▿u)=v1v2μ100μ2v1v2T.---(4)]]>

式中，μ₁和μ₂为梯度结构张量S_ρ(▽u)的两个特征值，且μ₁≥μ₂≥0；v₁和v₂为梯度结构张量的两个特征向量，且v₁对应于局部振幅变化最大的方向，即信号梯度方向，v₂对应于局部振幅变化最小的方向，即反射同相轴的取向；

然后利用梯度结构张量S_ρ(▽u)的两个特征值μ₁和μ₂计算线型信号结构的置信度量：

CL=μ1-μ2μ1+μ2.---(5)]]>

式中，CL为线型结构置信度量，在区间[0,1]之间取值；

最后结合线型结构置信度量CL与梯度结构张量S_ρ(▽u)的特征值μ₂计算地震横向不连续性置信度量：

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

式中，CI为横向不连续性置信度量，(1-CL)项为信号局部结构特征相对线型结构的背离程度，μ₂项为在均方误差最小意义下信号沿局部一致方向的能量变化强度；

3）构造结构自适应中值滤波器

构造具有椭圆形滤波窗口的结构自适应中值滤波窗：

M(x,y)={y∈(((y-x)·n(x))2σ12(x)+((y-x)·n⊥(x))2σ22(x)≤1)},---(8)]]>

式中，M(x,y)为当前滤波位置x处结构自适应中值滤波窗，y为当前滤波窗口包含的采样点的空间、时间位置坐标，·为向量内积算符，n(x)为地层倾向的方向向量，n_⊥(x)为地层梯度方向的方向向量，σ₁(x)为椭圆形滤波窗口的长轴，σ₂(x)为椭圆形滤波窗口的短轴；

利用地震线型结构置信度量CL和横向不连续性置信度量CI，确定与结构自适应中值滤波器的滤波尺度及滤波操作的方向选择性相关的尺度参数σ₁(x)和σ₂(x)：

σ1(x)=Rmax·g(CI(x))σ2(x)=(1-CL(x))σ1(x)---(9)]]>

式中，R_max为椭圆形滤波窗的最大尺寸参数，g(·)为关于CI(x)的单调减函数，且限定其取值范围为(0,1]，本发明取g(·)为指数函数

g(CI(x))=exp(-CI(x)β)---(10)]]>

式中，β为针对CI(x)的阈值参数，控制指数函数的衰减速度；

4）结构自适应中值滤波处理

首先根据待处理地震剖面内横向能量变化在时间和空间上的分布情况，将地震剖面Ω分割为若干具有一定边界重叠的子区域Ω＝∪Ω_i，在每个子区域Ω_i内取横向不连续性置信度量CI(x)的极大值，即

CImax,i=max{CI(x)|x∈Ωi},---(11)]]>

然后通过下式得到各个区域内的阈值参数β为

式中，α为百分比调整因子，thr为基底噪声阈值；

最后确定每个子区域Ω_i内各空间、时间位置x＝(x,t)处的中值滤波窗口的尺度参数σ₁(x)和σ₂(x)，从而进行如下的结构自适应中值滤波操作：

u^(x)=median{u(y),y∈M(x,y)}.---(7)]]>

式中，为位置x处的结构自适应中值滤波的输出结果，u(y)为位置x处的中值滤波窗M(x,y)内的各个采样点取值。

2.根据权利要求1所述的地层结构自适应中值滤波方法，其特征在于，步骤1）中，梯度结构张量S_ρ(▽u)的计算公式为：

Sρ(▿u)=Gρ⊗(▿u(▿u)T)=Gρ⊗(∂u∂x)2Gρ⊗(∂u∂x∂u∂t)Gρ⊗(∂u∂t∂u∂x)Gρ⊗(∂u∂t)2,---(3)]]>

式中，▿u=∂u∂x∂u∂t]]>为梯度向量，G_ρ为具有尺度参数ρ的二维高斯低通滤波函数G_ρ(x,t)＝exp(-(x²+t²)/2ρ²)，为卷积算符，T为矩阵转置算子。