CN104020492B - 一种三维地震资料的保边滤波方法 - Google Patents

Abstract

本发明公开了一种三维地震资料的保边滤波方法，包括：1)计算梯度结构张量；2)计算地层横向不连续性度量；3)构造各向异性拉普拉斯滤波器；4)各向异性拉普拉斯滤波处理。该方法从常规的非线性各向异性扩散方程出发，将非线性各向异性扩散方程的右端空间导数项分裂为关于待处理数据的二阶导数项和关于扩散张量一阶导数项，对两项进行取舍推导出各向异性拉普拉斯滤波器模型，能够有效衰减三维地震资料的随机噪声，增强地震同相轴的空间一致性，并且能够保护有效信号和断层、裂缝等地层边缘和细节信号结构，该技术方案易于实现，可操作性强。<pb pnum="1" />

Description

一种三维地震资料的保边滤波方法

技术领域

本发明属于地震勘探技术领域，涉及一种叠后三维地震资料的保边缘随机噪声衰减方法，尤其是一种利用各向异性拉普拉斯滤波器进行三维地震资料保边滤波。

背景技术

地震记录中的随机噪声通常是多种因素综合产生的无规则干扰波，在地震记录的时间、空间域内随机分布，其频带相对较宽，完全覆盖了有效波的频段。随着油气勘探程度的不断提高，挖掘油气储量主要来源于那些隐蔽性强、空间延伸小且非常薄的储层，相应增加了震资料解释的复杂程度，由简单的地质构造解释逐步转向精细的结构解释与储层预测。随着地震勘探目标日益复杂化，要求叠后地震资料的构造解释和储层预测工作更为精细，而地震资料的低信噪比是影响地震解释与储层预测可靠性的主要因素。常规地震资料处理中使用的去噪技术有小波变换去噪^[1]、f-x域预测滤波^[2]、KL变换^[3]、SVD分解^[4]等，这些方法大多利用了信号的空间相干特性，以牺牲横向分辨率达到提高信噪比的目的，容易使倾斜和弯曲同相轴受到衰减，并且会模糊和压制一些细微的信号结构，如小断裂、细河道等，甚至还能引起较大断距的断层两侧同相轴的错误连接，给精细的地震结构解释带来困难。因而，能够较好地保护地震数据的细微结构与断层、裂缝等边缘结构的叠后噪声衰减技术受到地球物理工作者的重视^[5-7]。

基于各种偏微分方程的滤波方法是上世纪九十年代初兴起的一种数字图像增强技术，尤其基于非线性各向异性扩散方程的滤波方法，由于其特有的边缘保持性能和结构倾向去噪，很快在多个领域得到广泛关注和应用。Fhemers和Hocker^[8]首先将扩散滤波技术用于地震资料的解释性处理中，提出称为VAN-GOGH的各向异性扩散滤波技术。Lavialle等^[9]结合面型和线型扩散滤波方程给出专门针对保护和增强断层构造的断层保持扩散滤波器。Kadlec等^[10]探索了能够保持层序特征的扩散滤波方法，用于增强地震数据中河道砂体特征的连续性，提高后续的河道自动识别和分割能力。由于非线性各向异性扩散滤波器的滤波性能受到多个参数选择的影响，Mueller^[11]给出扩散滤波器中保边缘项的新构造途径，降低了扩散滤波器对参数选择的依赖性。孙习平等^[12]通过考察有限差分法求解演变方程显式方法的优劣势，给出具有最优旋转不变性的各向异性扩散滤波步骤，可明显改善低信噪比资料的品质，突出地层接触关系，增强地震数据对层序体内部结构的成像能力。张尔华等^[13]利用梯度结构张量直接计算地震数据的横向不连续性，结合一种指数函数，灵活地控制各向异性扩散滤波器的边缘保护性能。Gomez^[14]认为基于各向异性扩散方程的滤波方法通过扩散时间的迭代推移形成对原始数据的滤波作用，由于不好估计出获得最佳滤波效果的扩散时间，导致无法准确预测扩散滤波器的滤波性态。

基于偏微分方程的滤波理论最初由Perona和Malik^[16]在数字图像处理领域提出，他们构造了如下扩散滤波器模型：

$\left\{\begin{array}{c}\frac{\&part;U}{\&part;t}=div\left(g\right(\left|\right|\&dtri;U\left|\right|)\&CenterDot;\&dtri;U)\\ U{|}_{t=0}=U_0\end{array}\right.,---\left(1\right)$

其中，t为扩散滤波方程的扩散时间，为扩散滤波器的滤波尺度参数；U为t时刻的扩散滤波结果；div为向量场的散度算子；g(||▽U||)为扩散滤波器的扩散系数，非负单调递减，且满足g(0)＝1，g(∞)＝0；U₀表示t＝0时刻的原始数据，作为扩散方程的初始条件。

在此基础上发展出很多非线性扩散滤波器模型，根据各自的特点应用于不同领域的图像、视频数据的保边缘去噪和增强处理。但是这些扩散模型所选取的扩散系数均为标量型，没有考虑到信号纹理的方向特性。为了能够保持和增强图像中的线型纹理结构，Weickert^[17,18]通过对图像信号的进行结构分析，为非线性扩散滤波器模型构造了张量型扩散系数，得到非线性各向异性扩散滤波方程：

$\left\{\begin{array}{c}\frac{\&part;U}{\&part;t}=div(D\&CenterDot;\&dtri;U)\\ U{|}_{t=0}=U_0\end{array}\right..---\left(2\right)$

其中，D为各向异性扩散滤波器的张量型扩散系数，即扩散张量。对于处理三维数据的扩散滤波方程，扩散张量D为3×3的对称半正定矩阵：

$D=\left(\begin{array}{ccc}{D}_{11}& D_{12}& D_{13}\\ D_{21}& D_{22}& D_{23}\\ D_{31}& D_{32}& D_{33}\end{array}\right)=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}{\mathrm{\&mu;}}_1& 0& 0\\ 0& {\mathrm{\&mu;}}_2& 0\\ 0& 0& {\mathrm{\&mu;}}_3\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T.---\left(3\right)$

其中，μ₁、μ₂和μ₃为扩散张量的三个非负特征值，在区间[0,1]内取值，分别定义了扩散滤波器沿特征方向v₁、v₂和v₃的滤波强度，μ_i＝0表示沿特征方向v_i无滤波操作，μ_i＝1表示沿特征方向v_i进行完全滤波操作。

以上现有技术具有以下缺点：

(1)常规各向异性扩散方程通过构造扩散张量约束扩散滤波操作，会导致扩散滤波操作的空间结构不明确，在一些特殊情况下会引起二义性；

(2)常规各向异性扩散方程将原始地震数据作为初始条件，通过扩散时间的迭代推移进行滤波，无法估计出获得最佳滤波效果的扩散时间。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种各向异性拉普拉斯滤波器的三维地震资料保边滤波方法，其从常规的非线性各向异性扩散方程出发，将非线性各向异性扩散方程的右端空间导数项分裂为关于待处理数据的二阶导数项和关于扩散张量一阶导数项，对两项进行取舍推导出各向异性拉普拉斯滤波器模型，能够有效衰减三维地震资料的随机噪声，增强地震同相轴的空间一致性，并且能够保护有效信号和断层、裂缝等地层边缘和细节信号结构，该技术方案易于实现，可操作性强。

本发明的目的是通过以下技术方案来解决的：

本发明三维地震资料的保边滤波方法，包括以下步骤：

1)计算梯度结构张量

首先利用有限差分方法计算三维地震数据U的梯度向量▽U，将梯度向量▽U与其转置向量(▽U)^T相乘得到初始张量矩阵，然后对初始张量矩阵的各个分量进行低通滤波处理，得到三维地震数据U对应的梯度结构张量J(▽U)；

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

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

$J(\&dtri;U)=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}{\mathrm{\&lambda;}}_1& 0& 0\\ 0& {\mathrm{\&lambda;}}_2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(10\right)$

式中，v₁、v₂和v₃为梯度结构张量的三个特征向量，构成一个局部正交坐标系，v₁指向信号的梯度方向，而由向量v₂和v₃张成的平面平行于局部构造平面，并且v₃对应于局部最一致性方向，即能量变化最小的方向；λ₁≥λ₂≥λ₃≥0对应于梯度结构张量的三个非负特征值，它们反应了信号沿特征方向的平均能量变化程度；

然后基于地下介质的层状结构分布假设，利用梯度结构张量J(▽U)的三个特征值λ₁、λ₂和λ₃计算三维地震数据的线型结构和面型结构的置信度量：

$C_{\mathrm{line}}=\frac{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}$ 和 $C_{\mathrm{plane}}=\frac{{\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_2}.---\left(11\right)$

式中，C_line为线型结构置信度量，C_plane为面型结构置信度量，两者均在区间[0,1]之间取值；

最后结合线型结构置信度量C_line与面型结构置信度量C_plane计算能够反应断层、裂缝等地质体边缘的断层置信度量：

C_fault＝C_line(1-C_plane)     (12)

式中，C_fault为断层置信度量，(1-C_plane)项为信号局部结构特征相对面型结构的背离程度；

3)构造各向异性拉普拉斯滤波器

根据梯度结构张量J(▽U)在高斯邻域G_ρ内得到的断层结构置信度量C_fault，构造自适应于信号结构的各向异性拉普拉斯滤波核函数为

$G\left(x_0,x\right)=\exp \left(-\frac{1}{2}\underset{i=1}{\overset{3}{\&Sigma;}}\frac{{\left(\left(x-x_0\right)\&CenterDot;v_i\left(x_0\right)\right)}^2}{{\mathrm{\&sigma;}}_i^2\left(x_0\right)}\right)=\exp \left(-\frac{{\left(x-x_0\right)}^T{D}^{-1}\left(x-x_0\right)}{2{\mathrm{\&sigma;}}_{\max }^2}\right).---\left(13\right)$

式中，x₀＝(x₀,y₀,t₀)为当前滤波点空间、时间位置坐标，x＝(x,y,t)为滤波邻域内采样点相对当前滤波点的空间、时间位置坐标，σ_i(x₀)，i＝1,2,3分别表示当前滤波点的各向异性拉普拉斯滤波核函数在由梯度结构张量J(▽U)所确定的局部正交空间的滤波尺度；D为各向异性拉普拉斯滤波方程的扩散张量，由下式描述：

$D=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& {\mathrm{\&sigma;}}_2/{\mathrm{\&sigma;}}_{\max }& 0\\ 0& 0& {\mathrm{\&sigma;}}_3/{\mathrm{\&sigma;}}_{\max }\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(15\right)$

式中，σ_max为尺度参数归一化因子；

根据信号结构自适应滤波的目的，利用三维地震数据的断层置信度量C_fault，构造各向异性拉普拉斯值滤波核函数的三个尺度参数σ₁(x₀)、σ₂(x₀)和σ₃(x₀)分别为：

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x_0\right)=0\\ {\mathrm{\&sigma;}}_2\left(x_0\right)={\mathrm{\&sigma;}}_{\min }+\left(1-C_{fault}\left(x_0\right)\right)\left({\mathrm{\&sigma;}}_{iso}-{\mathrm{\&sigma;}}_{\min }\right)\\ {\mathrm{\&sigma;}}_3\left(x_0\right)={\mathrm{\&sigma;}}_{iso}^2/{\mathrm{\&sigma;}}_2\left(x_0\right)\end{array}\right.,---\left(14\right)$

式中，σ_min为断层、裂缝等各向异性结构区域拉普拉斯滤波核函数横跨各向异性结构的最小尺寸，用以保证这些地层边缘区域的滤波性能；σ_iso为平坦地层区域各向同性的拉普拉斯滤波核的最大尺寸，用以避免过大的各向同性滤波核压制局部地层的起伏变化；

4)各向异性拉普拉斯滤波处理

将原始地震数据U作为初始条件，采用Newman边界条件，通过如下各向异性拉普拉斯滤波方程的扩散时间t的迭代推移，实现对原始地震数据U的滤波处理：

$\frac{\&part;U}{\&part;t}=trace\left(DH\right)={\mathrm{\&mu;}}_1{U}_{v_1{v}_1}+{\mathrm{\&mu;}}_2{U}_{v_2{v}_2}+{\mathrm{\&mu;}}_3{U}_{v_3{v}_3}---\left(7\right)$

式中，trace(·)为矩阵求迹算子，H为信号U的Hessian矩阵，为U沿特征方向v_i的二阶方向偏导数，μ₁＝0、μ₂＝σ₂/σ_max和μ₃＝σ₃/σ_max为扩散张量D的三个非负特征值，在区间[0,1]内取值，分别定义了各向异性拉普拉斯滤波器沿特征方向v₁、v₂和v₃的滤波强度，μ_i＝0表示沿特征方向v_i无滤波操作，μ_i＝1表示沿特征方向v_i进行完全滤波操作；

另外，还可以通过将原始地震数据U与各向异性拉普拉斯滤波核函数G(x₀,x)逐点卷积实现对原始地震数据U的滤波处理：

$U(x,t)=U(x,t=0)\&CircleTimes;G\left(x\right)---\left(8\right)$

式中，为卷积算符；

进一步，上述步骤1)中，梯度结构张量J(▽U)按照以下计算：

$J(\&dtri;U)=LP\&CircleTimes;J_0(\&dtri;U)=LP\&CircleTimes;(\&dtri;U{\left(\&dtri;U\right)}^T)---\left(9\right)$

式中， $\&dtri;U=\left[\begin{array}{ccc}\frac{\&PartialD;U}{\&PartialD;x}& \frac{\&PartialD;U}{\&PartialD;y}& \frac{\&PartialD;U}{\&PartialD;t}\end{array}\right]$ 为梯度向量，LP为三维空间的低通滤波函数，取为尺度为ρ的高斯低通滤波函数G_ρ(r)＝(2πρ²)^-3/2exp(-|r|²/(2ρ²))，为卷积算符，T为矩阵转置算子。

本发明具有以下有益效果：

本发明能够有效衰减三维地震资料的随机噪声和采集脚印等部分相干噪声，增强地震同相轴的空间一致性，并且能够保护有效信号和断层、裂缝等地层边缘和细节信号结构特征；该技术方案的算法避免了常规扩散滤波器的扩散时间选择问题，易于实现，可操作性好；利用梯度结构张量对地层结构的度量判定，对平面型地层区域设计大尺寸各向同性滤波核函数，提高平面地层区域的随机噪声衰减性能，对断裂区域设计狭长的滤波核函数，增强断层和裂缝等地层边缘结构的显示。

附图说明

图1为含有断层的三维模型5dB加噪数据滤波处理对比显示图；

图2为某油田三维地震资料的ANISO LAP滤波处理对比显示图；

图3为某油田纯波资料经ANISO LAP滤波处理前后的高分辨处理对比分析图；

图4为滤波结果的高分辨率处理剖面。

具体实施方式

该三维地震资料的保边滤波方法，包括以下步骤：

1)计算梯度结构张量

首先利用有限差分方法计算三维地震数据U的梯度向量▽U，将梯度向量▽U与其转置向量(▽U)^T相乘得到初始张量矩阵，然后对初始张量矩阵的各个分量进行低通滤波处理，得到三维地震数据U对应的梯度结构张量J(▽U)：

$J\left(\&dtri;U\right)=LP\&CircleTimes;J_0(\&dtri;U)=LP\&CircleTimes;(\&dtri;U{\left(\&dtri;U\right)}^T)---\left(9\right)$

式中， $\&dtri;U=\left[\begin{array}{ccc}\frac{\&PartialD;U}{\&PartialD;x}& \frac{\&PartialD;U}{\&PartialD;y}& \frac{\&PartialD;U}{\&PartialD;t}\end{array}\right]$ 为梯度向量，LP为三维空间的低通滤波函数，取为尺度为ρ的高斯低通滤波函数G_ρ(r)＝(2πρ²)^-3/2exp(-|r|²/(2ρ²))，为卷积算符，T为矩阵转置算子；

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

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

$J(\&dtri;U)=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}{\mathrm{\&lambda;}}_1& 0& 0\\ 0& {\mathrm{\&lambda;}}_2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(10\right)$

式中，v₁、v₂和v₃为梯度结构张量的三个特征向量，构成一个局部正交坐标系，v₁指向信号的梯度方向，而由向量v₂和v₃张成的平面平行于局部构造平面，并且v₃对应于局部最一致性方向，即能量变化最小的方向；λ₁≥λ₂≥λ₃≥0对应于梯度结构张量的三个非负特征值，它们反应了信号沿特征方向的平均能量变化程度；

然后基于地下介质的层状结构分布假设，利用梯度结构张量J(▽U)的三个特征值λ₁、λ₂和λ₃计算三维地震数据的线型结构和面型结构的置信度量：

$C_{\mathrm{line}}=\frac{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}$ 和 $C_{\mathrm{plane}}=\frac{{\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_2}.---\left(11\right)$

式中，C_line为线型结构置信度量，C_plane为面型结构置信度量，两者均在区间[0,1]之间取值；

最后结合线型结构置信度量C_line与面型结构置信度量C_plane计算能够反应断层、裂缝等地质体边缘的断层置信度量：

C_fault＝C_line(1-C_plane)     (12)

式中，C_fault为断层置信度量，(1-C_plane)项为信号局部结构特征相对面型结构的背离程度；

3)构造各向异性拉普拉斯滤波器

根据梯度结构张量J(▽U)在高斯邻域G_ρ内得到的断层结构置信度量C_fault，构造自适应于信号结构的各向异性拉普拉斯滤波核函数为

$G\left(x_0,x\right)=\exp \left(-\frac{1}{2}\underset{i=1}{\overset{3}{\&Sigma;}}\frac{{\left(\left(x-x_0\right)\&CenterDot;v_i\left(x_0\right)\right)}^2}{{\mathrm{\&sigma;}}_i^2\left(x_0\right)}\right)=\exp \left(-\frac{{\left(x-x_0\right)}^T{D}^{-1}\left(x-x_0\right)}{2{\mathrm{\&sigma;}}_{\max }^2}\right).---\left(13\right)$

式中，x₀＝(x₀,y₀,t₀)为当前滤波点空间、时间位置坐标，x＝(x,y,t)为滤波邻域内采样点相对当前滤波点的空间、时间位置坐标，σ_i(x₀)，i＝1,2,3分别表示当前滤波点的各向异性拉普拉斯滤波核函数在由梯度结构张量J(▽U)所确定的局部正交空间的滤波尺度；D为各向异性拉普拉斯滤波方程的扩散张量，由下式描述：

$D=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& {\mathrm{\&sigma;}}_2/{\mathrm{\&sigma;}}_{\max }& 0\\ 0& 0& {\mathrm{\&sigma;}}_3/{\mathrm{\&sigma;}}_{\max }\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(15\right)$

式中，σ_max为尺度参数归一化因子；

根据信号结构自适应滤波的目的，利用三维地震数据的断层置信度量C_fault，构造各向异性拉普拉斯值滤波核函数的三个尺度参数σ₁(x₀)、σ₂(x₀)和σ₃(x₀)分别为：

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x_0\right)=0\\ {\mathrm{\&sigma;}}_2\left(x_0\right)={\mathrm{\&sigma;}}_{\min }+\left(1-C_{fault}\left(x_0\right)\right)\left({\mathrm{\&sigma;}}_{iso}-{\mathrm{\&sigma;}}_{\min }\right)\\ {\mathrm{\&sigma;}}_3\left(x_0\right)={\mathrm{\&sigma;}}_{iso}^2/{\mathrm{\&sigma;}}_2\left(x_0\right)\end{array}\right.,---\left(14\right)$

式中，σ_min为断层、裂缝等各向异性结构区域拉普拉斯滤波核函数横跨各向异性结构的最小尺寸，用以保证这些地层边缘区域的滤波性能；σ_iso为平坦地层区域各向同性的拉普拉斯滤波核的最大尺寸，用以避免过大的各向同性滤波核压制局部地层的起伏变化；

4)各向异性拉普拉斯滤波处理

将原始地震数据U作为初始条件，采用Newman边界条件，通过如下各向异性拉普拉斯滤波方程的扩散时间t的迭代推移，实现对原始地震数据U的滤波处理：

$\frac{\&part;U}{\&part;t}=trace\left(DH\right)={\mathrm{\&mu;}}_1{U}_{v_1{v}_1}+{\mathrm{\&mu;}}_2{U}_{v_2{v}_2}+{\mathrm{\&mu;}}_3{U}_{v_3{v}_3}---\left(7\right)$

式中，trace(·)为矩阵求迹算子，H为信号U的Hessian矩阵，为U沿特征方向v_i的二阶方向偏导数，μ₁＝0、μ₂＝σ₂/σ_max和μ₃＝σ₃/σ_max为扩散张量D的三个非负特征值，在区间[0,1]内取值，分别定义了各向异性拉普拉斯滤波器沿特征方向v₁、v₂和v₃的滤波强度，μ_i＝0表示沿特征方向v_i无滤波操作，μ_i＝1表示沿特征方向v_i进行完全滤波操作；

另外，还可以通过将原始地震数据U与各向异性拉普拉斯滤波核函数G(x₀,x)逐点卷积实现对原始地震数据U的滤波处理：

$U(x,t)=U(x,t=0)\&CircleTimes;G\left(x\right)---\left(8\right)$

式中，为卷积算符。

下面结合附图对以上本发明的方法做进一步详细解释：

导向拉普拉斯滤波器

考虑式(2)中常规的张量型各向异性扩散方程，其右端的空间导数项经过简单的数学推导可以分裂为如下两项：

$\begin{array}{c}div\left(D\&dtri;U\right)={\mathrm{\&mu;}}_1{U}_{v_1{v}_1}+{\mathrm{\&mu;}}_2{U}_{v_2{v}_2}+{\mathrm{\&mu;}}_3{U}_{v_3{v}_3}+{\left(\&dtri;U\right)}^T\stackrel{\&RightArrow;}{div}\left(D\right)\\ =trace\left(DH\right)+{\left(\&dtri;U\right)}^T\stackrel{\&RightArrow;}{div}\left(D\right)\end{array}---\left(4\right)$

在式(4)中，为U沿特征方向v_i的二阶偏导数，H为信号U的Hessian矩阵，算子定义为

$\stackrel{\&RightArrow;}{div}\left(D\right)=\left(\begin{array}{c}div\left({\left[\begin{array}{ccc}{D}_{11}& D_{12}& D_{13}\end{array}\right]}^T\right)\\ div\left({\left[\begin{array}{ccc}{D}_{21}& D_{22}& D_{23}\end{array}\right]}^T\right)\\ div\left({\left[\begin{array}{ccc}{D}_{31}& D_{32}& D_{33}\end{array}\right]}^T\right)\end{array}\right).---\left(5\right)$

分析式(4)得到分解项，第一项trace(DH)是由扩散张量D确定的正交空间对信号U的加权拉普拉斯分解，这与经典的热传导方程右端的空间导数项相吻合^[19]，对应于由扩散张量D的各个特征分量所确定的局部滤波操作；第二项为信号U的梯度向量加权扩散张量D由式(5)所定义的矩阵散度，这项与扩散张量的空间变化相关。在非线性各向异性扩散方程中，采用扩散张量D确定的空间滤波结构，约束滤波过程以扩散张量的特征向量所确定的方向和其特征值所确定的滤波强度进行。因而，式(4)中第二项和各向异性扩散滤波器的设计初衷不相吻合，该项的存在会造成扩散滤波器的空间滤波结构不明确。这里我们给出一个典型的例子，分别选取两种完全不相同的扩散张量

$D_1=\frac{Id}{\left|\right|\&dtri;U\left|\right|}$

和

$D_2=\frac{(\&dtri;U\&dtri;U^T)}{\left|\right|\&dtri;U|{|}^3},$

此时，式(2)中扩散方程的右端项变为

$div(D_1\&CenterDot;\&dtri;U)=div(D_2\&CenterDot;\&dtri;U)=div\left(\frac{\&dtri;U}{\left|\right|\&dtri;U\left|\right|}\right)---\left(6\right)$

从式(6)式可见，各向同性扩散张量D₁和沿梯度方向的扩散张量D₂最终导出同一个扩散滤波器模型，即全变差最小的正则化模型。

为了使扩散滤波器具有更为清晰地空间滤波结构，使得滤波操作和扩散张量D所确定的空间滤波结构一一对应，本文去除扩散方程右端项中与扩散张量D的空间变化相关联项，得到如下的扩散滤波模型

$\frac{\&part;U}{\&part;t}=trace\left(DH\right)={\mathrm{\&mu;}}_1{U}_{v_1{v}_1}+{\mathrm{\&mu;}}_2{U}_{v_2{v}_2}+{\mathrm{\&mu;}}_3{U}_{v_3{v}_3}---\left(7\right)$

方程(7)就是导向拉普拉斯发展方程，其滤波性能由扩散张量D唯一确定。通过求解导向拉普拉斯发展方程在t时刻的解，即获得对应尺度的滤波结果。可从不同途径求解方程(7)，一种是采用经典的数值计算方法，采用中心差分法计算U的空间导数项得到梯度向量和Hessian矩阵，而一阶时间导数可采用前向差分法或者后向差分法计算，从而通过时间步长的迭代得到滤波结果。此外，方程(7)还存在如下解析解：

$U(x,t)=U(x,t=0)*G_t^D\left(x\right)=U(x,t=0)*\left\{\frac{1}{4\mathrm{\&pi;}t}\exp (-\frac{x^T{D}^{-1}x}{4t})\right\}---\left(8\right)$

此时，通过由扩散张量D确定的定向高斯核函数与数据U逐点卷积就获得滤波结果。

结构自适应拉普拉斯滤波器

导向拉普拉斯滤波器的滤波性能主要取决于扩散张量的构造，设计合理的扩散张量是对地震数据作保边滤波处理的关键。本文采用梯度结构张量法估计地震反射层的方向，度量地层结构的规则性，在此基础上给出能够压制干扰噪声和增强地质结构特征的扩散张量。对于三维地震数据U，构造其对应的梯度结构张量为：

J(▽U)＝LP*J₀(▽U)＝LP*(▽U(▽U)^T)     (9)

其中，LP为三维空间的低通滤波函数，取为尺度为ρ的高斯函数；J₀(▽U)表示初始梯度张量，由▽U和其转置向量构成。

高斯低通滤波函数G_ρ(r)＝(2πρ²)^-3/2exp(-|r|²/(2ρ²))与初始张量J₀(▽U)的每个分量通过卷积作用，确定了梯度结构张量J(▽U)可度量的信号特征的最大尺度。

梯度结构张量J(▽U)为3×3的对称半正定矩阵，对其作矩阵的特征分解：

$J(\&dtri;U)=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}{\mathrm{\&lambda;}}_1& 0& 0\\ 0& {\mathrm{\&lambda;}}_2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(10\right)$

其中，v₁、v₂和v₃为梯度结构张量的三个特征向量，构成一个局部正交坐标系，v₁指向信号的梯度方向，而由向量v₂和v₃张成的平面平行于局部构造平面，并且v₃对应于局部最一致性方向，即能量变化最小的方向；λ₁≥λ₂≥λ₃≥0对应于梯度结构张量的三个非负特征值，它们反应了信号沿特征方向的平均能量变化程度。基于地下介质的层状结构分布假设，梯度结构张量方法可有效识别地震数据中的线型结构和面型结构，Bakker^[21]结合梯度结构张量的三个特征值给出的信息定义了地震数据中线型结构和面型结构的置信度量：

$C_{\mathrm{line}}=\frac{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}$ 和 $C_{\mathrm{plane}}=\frac{{\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_2}.---\left(11\right)$

联合线型结构置信度量C_line和面型结构置信度量C_plane，Bakker给出能够反应断层、裂缝等地质体边缘的断层置信度量：

C_fault＝C_line(1-C_plane)     (12)

在同相轴延展良好的区域，由于C_line→0且C_plane→1，使得C_fault→0；而在断层、裂缝等地层边缘区域，由于C_line→1且C_plane→0，使得C_fault→1。

根据在梯度结构张量的高斯邻域G_ρ内获得的不同信号结构的置信度量值，构造自适应于信号结构的各向异性拉普拉斯滤波核函数为

$G\left(x_0,x\right)=\exp \left(-\frac{1}{2}\underset{i=1}{\overset{3}{\&Sigma;}}\frac{{\left(\left(x-x_0\right)\&CenterDot;v_i\left(x_0\right)\right)}^2}{{\mathrm{\&sigma;}}_i^2\left(x_0\right)}\right)=\exp \left(-\frac{{\left(x-x_0\right)}^T{D}^{-1}\left(x-x_0\right)}{2{\mathrm{\&sigma;}}_{\max }^2}\right).---\left(13\right)$

式(13)中，x₀＝(x₀,y₀,t₀)为当前的分析点坐标，σ_i(x₀)，i＝1,2,3分别表示当前分析点的各向异性拉普拉斯滤波函数在由梯度结构张量J(▽U)所确定的局部正交空间的滤波尺度，根据信号结构自适应滤波的目的，构造如下形式的尺度参数：

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x_0\right)=0\\ {\mathrm{\&sigma;}}_2\left(x_0\right)={\mathrm{\&sigma;}}_{\min }+\left(1-C_{fault}\left(x_0\right)\right)\left({\mathrm{\&sigma;}}_{iso}-{\mathrm{\&sigma;}}_{\min }\right)\\ {\mathrm{\&sigma;}}_3\left(x_0\right)={\mathrm{\&sigma;}}_{iso}^2/{\mathrm{\&sigma;}}_2\left(x_0\right)\end{array}\right.,---\left(14\right)$

式(14)中，σ_min是在断层、裂缝等各向异性结构区域拉普拉斯滤波核横跨各向异性结构的最小尺寸，用以保证这些地层边缘区域的滤波性能；σ_iso为在平坦地层区域各向同性的拉普拉斯滤波核的最大尺寸，用以避免过大的各向同性滤波核压制局部地层的起伏变化。此时，对应的扩散张量D由下式描述：

$D=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& {\mathrm{\&sigma;}}_2/{\mathrm{\&sigma;}}_{\max }& 0\\ 0& 0& {\mathrm{\&sigma;}}_3/{\mathrm{\&sigma;}}_{\max }\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(15\right)$

联合分析式(13)和(14)，在断层、裂缝等地层边缘结构区域，由于C_fault→1，故有σ₂(x₀)→σ_min和即拉普拉斯滤波核沿局部地层边缘方向拉伸为强各向异性的线条状形态，使得滤波操作仅仅沿着地层边缘进行；在反射同相轴延展良好的区域，由于C_fault→0，故有σ₂(x₀)→σ_iso和σ₃(x₀)→σ_iso，即拉普拉斯滤波核沿局部地层平面均匀拉伸为各向同性的圆饼状形态，使得该区域的滤波程度最大化；对于同相轴内存在振幅平缓变化或者同相轴有一定程度起伏的区域，由于0＜C_fault＜1，拉普拉斯滤波核沿局部地层结构呈现为弱各向异性形态，从而既保证这些区域的滤波性能，又不造成对有效能量变化或者地层几何形态的压制。在实际地震资料的处理中，可将σ_min取为地震资料中的随机噪声干扰的尺度，而σ_iso选取为与计算梯度结构张量的低通高斯函数G_ρ的尺度ρ相当。

数值仿真结果

合成模型资料

图1a和b分别为加噪三维合成模型的垂直和水平切面。该合成模型沿主测线、联络测线以及垂直方向都有150个采样点，由两种同频、沿不同传播方向的平面波组成，在两种平面波的交界处形成弯曲的倾斜断层；此外，该模型还有一处与此倾斜断层相交的垂直断层。向原始模型数据加不同信噪比的高斯白噪声，分别采用经典的VAN GOGH扩散滤波器^[8]和本文的各向异性拉普拉斯滤波器(ANISO LAP)处理加噪的合成模型。两种方法中梯度结构张量的计算取相同的参数ρ＝3.0，并当滤波结果的整体信噪比达到最大时终止VAN GOGH扩散滤波器的迭代过程，而ANISOLAP滤波器的参数σ_iso＝3.0，σ_min＝0.5。为了精细对比两种滤波器的滤波性能和保边缘性能，本文采用梯度结构张量计算断层置信度量C_fault，通过对C_fault设定阈值将合成数据分为断层和非断层两个区域进行分析。表1中给出两种滤波器的处理结果在断层和非断层区域以及整体信噪比的对比结果，图1给出5dB信噪比情况下两种滤波器处理结果。从两种滤波结果的信噪比来看，本文的ANISO LAP滤波器对断层的保持效果明显优于VAN GOGH扩散滤波器，且在非断层区域取得与VAN GOGH扩散滤波器相当的滤波性能，这说明ANISO LAP滤波器不会因为保边缘约束而过度降低了在非边缘结构处的滤波性能；对比滤波结果的纵向和水平切片，两种滤波器均在整体上取得优异的滤波效果，但是，比较两种滤波器得到的纵向切片可以看出，VANGOGH扩散滤波器已经引起断层两侧同相轴的错误连接，而ANISO LAP滤波器相对较好地保持了两种平面波边界的尖锐性；另外，比较两者得到的水平切片可以看出，VAN GOGH扩散滤波器对圆圈选定的边缘突起结构造成一定压制，而ANISO LAP滤波器较好地保持了该结构的本来面目。因而，相比经典的VAN GOGH扩散滤波器，本文的ANISO LAP滤波器更好地兼顾了不连续性保持和优异的滤波性能。

表1 VAN GOGH和各向异性LAP滤波器滤波结果的信噪比

实际地震资料

将本文提出的ANISO LAP滤波器用于某油田三维地震资料的噪声衰减与结构特征增强。在该三维区块中，一个复杂的断层网络横截几条河道砂体沉积特征，因而要求滤波器具有优异的边缘和细微结构保持能力，否则很容易引起细小断层和河道的损伤和压制。图2a是该三维纯波地震资料的一条主测线，图2b是本文的ANISOLAP滤波器的处理结果。对比处理前后的剖面可以看出，原始地震剖面内的随机噪声和采集脚印等相干干扰噪声得到很好地衰减，处理后的地震剖面的信噪比得到明显的提高，地层的连续性得到增强，断层和河道变得更为清晰；对比图2c和d中的滤波前后的地震相干体剖面可以看出，经过ANISO LAP滤波处理后，地震相干剖面对于断层的刻画更加清晰准确，相干剖面的信噪比得到增强，对于一些由采集脚印噪声引起的虚假断层信息也得到很好地抑制。进一步对比的三维纯波地震资料处理前后的某油层组的沿层振幅切片来分析ANISO LAP滤波器的滤波效果。图2e是原始三维纯波地震资料的沿层振幅切片，从图中可以看到两条北东向的河流贯穿整个工区，在工区的西北和西南还有小型的网状河道；此外，从沿层切片的许多地方存在东西向的灰白色条带，这不应该是实际的沉积地质现象，而是采集脚印等高频噪声的反映。由于受到随机噪声和采集脚印的严重干扰，图2e中的沉积现象难以被准确识别。图2f是ANISO LAP滤波器处理结果的沿层振幅切片，与图2e相比较容易发现，随机噪声和采集脚印等高频干扰得到有效压制，切片中的沉积地质特征变得非常清晰，尤其工区西侧的网状河流可以很容易被识别。

将本文提出的ANISO LAP滤波器用于某油田海上工区的三维叠后地震资料的噪声衰减与结构特征增强，并分析其对叠后地震资料高分辨率处理效果产生的影响。图3a是该三维工区的一条测线，图3b是ANISO LAP的滤波结果。对比滤波前后的地震剖面可以看出，干扰噪声得到很大程度的衰减，几处明显的采集脚印噪声得到很好地压制，地层的连续性得到明显增强，地质体的边缘变得更为清晰。对滤波前后的数据分别进行高分辨率处理，得到的处理结果分别如图3c所示，可以看出两种结果中地震剖面的纵向分辨率均得到明显提高，但是对滤波前的地震资料直接进行高分辨率处理，由于受到随机噪声和采集脚印等高频噪声的干扰，高分辨率处理会进一步放大这些高频干扰噪声，使得高分辨率结果的信噪比很低，地层的连续性变得更差。对比图3c和图4中经过高分辨率处理的两个剖面的视分辨率可见，滤波前的高分辨率处理结果与滤波后的高分辨处理结果的纵向分辨几乎完全一致，在图3c中能够观察到反射同相轴这表明本文的ANISO LAP滤波方法在滤波过程能够很好地保持信号的频谱结构，不会降低地震资料的纵向分辨率。进一步比较图3c中黑框选定的目标地层区域，虽然滤波前的剖面经过高分辨率处理后两组地层沉积情况变得更加细致，但是地层间的接触关系由于受到高频噪声的干扰变得模糊不清，甚至会引起误判；而从滤波后地震资料的高分辨率处理结果可以看出，由于干扰噪声得到预处理压制，高分辨率处理结果仍然保持了较高的信噪比，黑框选定的两组沉积地层的层间接触关系清晰可辨，有利于对该目标区域的储层情况进行深入分析。

Claims (2)

1.一种三维地震资料的保边滤波方法，其特征在于，包括以下步骤：

1)计算梯度结构张量

首先利用有限差分方法计算三维地震数据U的梯度向量将梯度向量与其转置向量相乘得到初始张量矩阵，然后对初始张量矩阵的各个分量进行低通滤波处理，得到三维地震数据U对应的梯度结构张量

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

首先对梯度结构张量进行矩阵特征分解：

$J(\&dtri;U)=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}{\mathrm{\&lambda;}}_1& 0& 0\\ 0& {\mathrm{\&lambda;}}_2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(10\right)$

式中，v₁、v₂和v₃为梯度结构张量的三个特征向量，构成一个局部正交坐标系，v₁指向信号的梯度方向，而由向量v₂和v₃张成的平面平行于局部构造平面，并且v₃对应于局部最一致性方向，即能量变化最小的方向；λ₁≥λ₂≥λ₃≥0对应于梯度结构张量的三个非负特征值，它们反应了信号沿特征方向的平均能量变化程度；

然后基于地下介质的层状结构分布假设，利用梯度结构张量的三个特征值λ₁、λ₂和λ₃计算三维地震数据的线型结构和面型结构的置信度量：

$C_{line}=\frac{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}$ 和 $C_{plane}=\frac{{\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_2}.---\left(11\right)$

式中，C_line为线型结构置信度量，C_plane为面型结构置信度量，两者均在区间[0,1]之间取值；

最后结合线型结构置信度量C_line与面型结构置信度量C_plane计算能够反应断层、裂缝地质体边缘的断层置信度量：

C_fault＝C_line(1-C_plane)          (12)

式中，C_fault为断层置信度量，(1-C_plane)项为信号局部结构特征相对面型结构的背离程度；

3)构造各向异性拉普拉斯滤波器

根据梯度结构张量在高斯邻域G_ρ内得到的断层结构置信度量C_fault，构造自适应于信号结构的各向异性拉普拉斯滤波核函数为

$G(x_0,x)=\exp (-\frac{1}{2}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}\frac{{\left(\right(x-x_0)\&CenterDot;v_i(x_0\left)\right)}^2}{{\mathrm{\&sigma;}}_i^2\left(x_0\right)})=\exp (-\frac{{(x-x_0)}^T{D}^{-1}(x-x_0)}{2{\mathrm{\&sigma;}}_{\max }^2}).---\left(13\right)$

式中，x₀＝(x₀,y₀,t₀)为当前滤波点空间、时间位置坐标，x＝(x,y,t)为滤波邻域内采样点相对当前滤波点的空间、时间位置坐标，σ_i(x₀)，i＝1,2,3分别表示当前滤波点的各向异性拉普拉斯滤波核函数在由梯度结构张量所确定的局部正交空间的滤波尺度；D为各向异性拉普拉斯滤波方程的扩散张量，由下式描述：

$D=\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& {\mathrm{\&sigma;}}_2/{\mathrm{\&sigma;}}_{\max }& 0\\ 0& 0& {\mathrm{\&sigma;}}_3/{\mathrm{\&sigma;}}_{\max }\end{array}\right){\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)}^T---\left(15\right)$

式中，σ_max为尺度参数归一化因子；

根据信号结构自适应滤波的目的，利用三维地震数据的断层置信度量C_fault，构造各向异性拉普拉斯值滤波核函数的三个尺度参数σ₁(x₀)、σ₂(x₀)和σ₃(x₀)分别为：

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1\left(x_0\right)=0\\ {\mathrm{\&sigma;}}_2\left(x_0\right)={\mathrm{\&sigma;}}_{min}+(1-C_{fault}(x_0))({\mathrm{\&sigma;}}_{iso}-{\mathrm{\&sigma;}}_{min}),\\ {\mathrm{\&sigma;}}_3\left(x_0\right)={\mathrm{\&sigma;}}_{iso}^2/{\mathrm{\&sigma;}}_2\left(x_0\right)\end{array}\right.---\left(14\right)$

式中，σ_min为断层、裂缝各向异性结构区域拉普拉斯滤波核函数横跨各向异性结构的最小尺寸，用以保证这些地层边缘区域的滤波性能；σ_iso为平坦地层区域各向同性的拉普拉斯滤波核的最大尺寸，用以避免过大的各向同性滤波核压制局部地层的起伏变化；

4)各向异性拉普拉斯滤波处理

将原始地震数据U作为初始条件，采用Newman边界条件，通过如下各向异性拉普拉斯滤波方程的扩散时间t的迭代推移，实现对原始地震数据U的滤波处理：

$\frac{\&part;U}{\&part;t}=trace\left(DH\right)={\mathrm{\&mu;}}_1{U}_{v_1{v}_1}+{\mathrm{\&mu;}}_2{U}_{v_2{v}_2}+{\mathrm{\&mu;}}_3{U}_{v_3{v}_3}---\left(7\right)$

式中，trace(·)为矩阵求迹算子，H为信号U的Hessian矩阵，i＝1,2,3为U沿特征方向v_i的二阶方向偏导数，μ₁＝0、μ₂＝σ₂/σ_max和μ₃＝σ₃/σ_max为扩散张量D的三个非负特征值，在区间[0,1]内取值，分别定义了各向异性拉普拉斯滤波器沿特征方向v₁、v₂和v₃的滤波强度，μ_i＝0表示沿特征方向v_i无滤波操作，μ_i＝1表示沿特征方向v_i进行完全滤波操作；

另外，通过将原始地震数据U与各向异性拉普拉斯滤波核函数G(x₀,x)逐点卷积实现对原始地震数据U的滤波处理：

$U(x,t)=U(x,t=0)\&CircleTimes;G\left(x\right)---\left(8\right)$

式中，为卷积算符。

2.根据权利要求1所述的三维地震资料的保边滤波方法，其特征在于，步骤1)中，梯度结构张量按照以下计算：

$J(\&dtri;U)=LP\&CircleTimes;J_0(\&dtri;U)=LP\&CircleTimes;(\&dtri;U{(\&dtri;U)}^T)---\left(9\right)$

式中，为梯度向量，LP为三维空间的低通滤波函数，取尺度为ρ的高斯低通滤波函数G_ρ(r)＝(2πρ²)^-3/2exp(-|r|²/(2ρ²))，为卷积算符，T为矩阵转置算子。