0 引言

通常断块型油藏的地质特征较复杂，其地震响应经过处理后仍存在噪声，特别是断层附近的地震反射数据的噪声能量更强。当前的绝大多数断层检测方法主要利用断层处的地震反射信息急速变化识别与描述断层，对噪声较敏感，因此需要预先压制地震噪声，同时保护断层等的边缘信息。

采用传统的滤波方法(中值滤波、F-X滤波、高斯滤波等)^[1-2]去除地震噪声，不能有效保护断层等不连续的边界信息，因而人们广泛研究了保边滤波技术。目前，各向异性扩散滤波是保边去噪处理中应用较广泛的一种滤波技术，其扩散方程的形式为偏微分方程，具有较高的计算效率和稳定性。Perona等^[3]在各向同性扩散方程的基础上提出了各向异性扩散方程，使扩散滤波在不同方向上具有不同的扩散速度，被称为Perona-Malik(简称P-M)扩散模型；Weickert^[4]系统研究了各向异性扩散理论，并结合相干体技术提出了相干增强扩散算法，模型试算及实际数据应用都取得了较好的效果；Fehmers等^[5]首次将各向异性扩散滤波技术用于地震勘探领域，在保护断层等构造信息的同时，压制了噪声；Hale^[6-7]提出了基于结构导向的各向异性扩散滤波，可以有效提高资料信噪比；问雪等^[8]将构造导向平滑方法应用于断层解释，提高了断层解释精度；Lavialle等^[9]改进了三维扩散张量特征值的设计方法，提出了保护断层信息的扩散滤波方法，使断层信息更突出；Chopra等^[10]利用各向异性扩散滤波对地震数据预处理，取得了很好的效果；杨培杰等^[11]提出了一种方向性边界信息保护的滤波方法，在保护断层信息的同时，提高了资料信噪比；孙夕平等^[12]提出了一种具有最优旋转不变性的各向异性扩散滤波方法，可有效提高地震资料信噪比，突出地层接触关系，增强地震数据对层序体内部结构的成像能力；张尔华等^[13]将各向异性扩散滤波方法拓展到三维；严哲等^[14]在各向异性扩散滤波算法中将相干值作为一个参数引入具体算法，提高了对横向不连续点的保护性能；杨千里等^[15]、姚振岸等^[16-17]应用各向异性扩散滤波对地震数据进行保边去噪，可以清晰、准确地识别断层；李福强等^[18]提出了基于断层算子的各向异性扩散滤波方法，可有效去除噪声，自适应增强断点处的能量，提高了断层识别精度；李军等^[19]提出了图像熵各向异性扩散滤波方法，可定量给出构建结构张量矩阵需要添加的二阶导数占比，更合理地判断扩散方向。

本文在前人的研究基础上，分析了结构张量特征值与三维地震图像的局部结构特征之间的关系，利用断层置信度参数重新构建扩散张量特征值，进而重构扩散张量，实现基于不连续性的三维各向异性扩散滤波方法，可有效提高地震资料信噪比，保护断层等不连续的边界信息。

1 方法原理 1.1 各向异性扩散滤波方程

各向异性扩散滤波的扩散过程等价于热传导过程，表示为^[20]

$ \frac{{\partial U}}{{\partial t}} = {\rm{div}}\left( {g\nabla U} \right) $

(1)

式中：U为输入的地震数据振幅信息；t为扩散时间；g为扩散系数。g的取值有三种情况：①当g取常数时，式(1)为线性各向同性方程，对数据进行滤波后，能够提高信噪比，但不能保护断层等地质体的边界及方向信息；②当g为非负单调递减函数^[3]时，式(1)为非线性各向同性方程，滤波后容易产生阶梯效应和针孔效应，而在强噪声区域无去噪效果；③使用扩散张量D作为扩散系数^[4]，可以根据各向异性信息对数据进行自适应滤波，其三维各向异性扩散滤波方程为

$ \left\{ \begin{array}{l} \frac{{\partial U\left( {x, y, z, t} \right)}}{{\partial t}} = {\rm{div}}\left[ {\mathit{\boldsymbol{D}}\cdot\nabla U\left( {x, y, z, t} \right)} \right]{\rm{ }}\\ U\left( {x, y, z, t} \right) = {U_0}\left( {x, y, z} \right) \end{array} \right. $

(2)

式中U₀(x, y, z)为原始三维地震数据，(x，y，z)为某点的空间坐标。

1.2 结构张量的构建

在构建D的过程中，一般采用具有对称半正定的结构张量S_ρ计算局部方向信息，其形式为

$ {{\bf{S}}_\rho } = {G_\rho }*(\nabla U_\sigma ^{\rm{T}}\nabla {U_\sigma }) $

(3)

式中：$ \nabla {U_\sigma }$为数据的梯度信息，U_σ=G_σ*U，“*”为褶积运算符，σ为噪声尺度；ρ为整合尺度，一般大于σ；G_ρ(或G_σ)为高斯核函数，即

$ \begin{array}{l} {G_\rho }\left( {x, y, z} \right) = \frac{1}{{{{(2{\rm{ \mathit{ π} }}{\rho ^2})}^{\frac{3}{2}}}}} \times \\ \;\;\;{\rm{exp}}\left[ { - \frac{{{{(x - {x_0})}^2} + {{(y - {y_0})}^2} + {{(z - {z_0})}^2}}}{{2{\rho ^2}}}} \right] \end{array} $

(4)

式中(x₀, y₀, z₀)为参考点空间坐标。在断层位置处经过高斯滤波后地震数据对空间的二阶导数一般具有局部极大值，将该信息添加到S_ρ的计算公式中，可将式(3)改写为

$ \begin{array}{*{20}{l}} {{{\bf{S}}_\rho } = {G_\rho }*[\nabla U_\sigma ^{\rm{T}}\nabla {U_\sigma } + \alpha {{(\frac{{{\partial ^2}{U_\sigma }}}{{\partial {x^2}}}\frac{{{\partial ^2}{U_\sigma }}}{{\partial {y^2}}}\frac{{{\partial ^2}{U_\sigma }}}{{\partial {z^2}}})}^{\rm{T}}} \times }\\ {\;\;\;\;\;\;\;\;(\frac{{{\partial ^2}{U_\sigma }}}{{\partial {x^2}}}\frac{{{\partial ^2}{U_\sigma }}}{{\partial {y^2}}}\frac{{{\partial ^2}{U_\sigma }}}{{\partial {z^2}}})]} \end{array} $

(5)

式中α∈[0, 1]为稳定系数，用于控制高斯滤波后地震数据对空间的二阶导数的贡献率。

对S_ρ进行特征分解，得

$ {\mathit{\boldsymbol{S}}_\rho } = \varepsilon [{\mathit{\boldsymbol{\zeta }}_1}\;\;{\mathit{\boldsymbol{\zeta }}_2}\;\;{\mathit{\boldsymbol{\zeta }}_3}]\left[ {\begin{array}{*{20}{c}} {{\lambda _1}}&0&0\\ 0&{{\lambda _2}}&0\\ 0&0&{{\lambda _3}} \end{array}} \right]{[{\mathit{\boldsymbol{\zeta }}_1}\;\;{\mathit{\boldsymbol{\zeta }}_2}\;\;{\mathit{\boldsymbol{\zeta }}_3}]^{\rm{T}}} $

(6)

式中：λ₁、λ₂和λ₃为S_ρ的三个特征值，且λ₁>λ₂>λ₃；ζ₁、ζ₂和ζ₃分别为λ₁、λ₂和λ₃对应的特征向量，ζ₁指示梯度变化最大的方向，ζ₃指示梯度变化最小的方向；ε为横向不连续因子，在地层连续处ε趋近于1，在断层位置处ε趋近于0。

Bakker^[21]指出，对于三维地震数据来说，当结构张量的特征值λ_i(i=1，2，3)为0时，表示图像沿该特征向量的方向移不变。当沿一个方向移不变时，为线状线性结构；当沿两个方向移不变时，为面状线性结构。根据特征值的大小及其对应的不同的特征向量可以将地震图像结构表述为以下四类：

(1) λ₁≈λ₂≈λ₃≈0，三个特征向量方向的变化很小或者基本没有变化，指示平滑区域，表示地质体呈各向同性；

(2) λ₁＞λ₂≈λ₃≈0，在ζ₁方向梯度变化较大，另外两个特征方向的梯度变化很小且几乎相等，指示平面线性地质体，对应模型为近似面状线性结构；

(3) λ₁＞λ₂＞λ₃≈0，其中一个方向梯度变化近似为零，指示在断层等不连续的边界处出现不连续的间断信息，对应模型为线状结构模型；

(4) λ₁＞λ₂＞λ₃＞0，三个特征向量方向的梯度变化都不同，并且均不为零，说明区域各向异性特征明显，表示地质体的不整合结构。

van Kempen等^[22]根据结构张量特征值与地震图像结构的关系，用线状置信度量M_line和面状置信度量M_plane表征图像结构^[18]，即

$ \left\{ \begin{array}{l} {M_{{\rm{line}}}} = \frac{{{\lambda _2} - {\lambda _3}}}{{{\lambda _2} + {\lambda _3}}}\\ {M_{{\rm{plane}}}} = \frac{{{\lambda _1} - {\lambda _2}}}{{{\lambda _1} + {\lambda _2}}} \end{array} \right. $

(7)

M_line∈[0, 1]、M_plane∈[0, 1]。M_line描述某点邻域与线状结构的相似程度，当M_line=1时，指示线状结构；M_plane描述某点邻域与面状结构的相似程度，当M_plane=1时，指示面状结构。由M_line、M_plane得到检测断层等不连续结构边界信息的断层置信度量M_fault^{[21, 23]}

$ {M_{{\rm{fault}}}} = {M_{{\rm{line}}}}(1 - {M_{{\rm{plane}}}}) = \frac{{2{\lambda _2}\left( {{\lambda _2} - {\lambda _3}} \right)}}{{\left( {{\lambda _1} + {\lambda _2}} \right)({\lambda _2} + {\lambda _3})}} $

(8)

表 1为置信度参数与地层结构以及特征值之间的关系。

在地震剖面上，如果地震同相轴的连续性较好，则M_line趋于0、M_plane趋于1、M_fault趋于0；同理，在存在断层等构造的同相轴不连续区域，M_line趋于1、M_plane趋于0、M_fault趋于1。

1.3 利用M_fault设计扩散张量特征值

构建扩散张量D时，为了保证扩散方向沿着构造方向，D的特征向量应与S_ρ的一致，因此

$ \mathit{\boldsymbol{D}} = \varepsilon [{\mathit{\boldsymbol{\zeta }}_1}\;\;{\mathit{\boldsymbol{\zeta }}_2}\;\;{\mathit{\boldsymbol{\zeta }}_3}]\left[ {\begin{array}{*{20}{c}} {{\gamma _1}}&0&0\\ 0&{{\gamma _2}}&0\\ 0&0&{{\gamma _3}} \end{array}} \right]{[{\mathit{\boldsymbol{\zeta }}_1}\;\;{\mathit{\boldsymbol{\zeta }}_2}\;\;{\mathit{\boldsymbol{\zeta }}_3}]^{\rm{T}}} $

(9)

式中γ₁、γ₂和γ₃为D的三个特征值，对应的特征向量也为ζ₁、ζ₂和ζ₃。

在三维各向异性扩散滤波中，扩散滤波的强度取决于γ₁、γ₂和γ₃。特征值越接近0，扩散滤波的强度越小；特征值趋于1，扩散滤波强度最大

为了克服常规扩散滤波方法不能保护数据边缘信息的缺点，在D的特征值设计中引入M_line、M_plane以及M_fault，即

$ \left\{ \begin{array}{l} {\gamma _1} = a{\rm{ }}\\ {\gamma _2} = \left\{ \begin{array}{l} a\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\lambda _1} = {\lambda _2}\\ a + \left( {1 - a} \right){h_\tau }({M_{{\rm{plane}}}})\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right.\\ {\gamma _3} = \left\{ \begin{array}{l} {\gamma _2}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\lambda _2} = {\lambda _3}\\ {\gamma _2} + ({M_{{\rm{fault}}}} - {\gamma _2}){h_\tau }({M_{{\rm{line}}}})\;\;\;\;\;\;\;其他 \end{array} \right. \end{array} \right. $

(10)

式中：a为一个接近于0的正数；$ {h_\tau }\left( x \right) = \frac{{{\rm{tanh}}\left[ {s\left( {x - \tau } \right)} \right] + {\rm{1}}}}{{{\rm{tanh}}\left[ {s\left( {1 - \tau } \right)} \right] + 1}}$为S型函数^[24-25]，s和τ分别为门槛值和斜率，h_τ(x)可以有效控制扩散滤波在两个特征向量之间的强度大小。

式(10)表明：①沿梯度变化最大的方向，即沿ζ₁方向，γ₁=a，该方向的扩散系数趋于0，即扩散滤波的强度较小，可以保护断层等不连续结构的边界信息；②沿ζ₂和ζ₃的方向，即垂直于最大梯度的方向，可以根据M_line、M_plane以及M_fault调整特征值。如果M_fault趋于0，地层为面状结构，γ₂和γ₃的值趋于1，沿ζ₂和ζ₃组成的平面扩散模型滤波强度大；如果M_fault趋于1，γ₁≈γ₂≈a，沿ζ₁、ζ₂方向扩散滤波强度小，可以有效保护断层等不连续的边界信息。

综上所述，利用M_fault设计D特征值的方法可在增强扩散滤波作用的同时，较好地保护断层等特殊地质结构的边界信息。

1.4 方法实现

计算出D后，通过迭代法求解三维各向异性扩散方程，最终得到扩散滤波后的地震数据

$ \begin{array}{*{20}{l}} {{U^{(k + 1)}} = {U^{(k)}} + \Delta t \cdot {\rm{div}}[\mathit{\boldsymbol{D}} \cdot \nabla {U^{(k)}}]}\\ {\;\;\;\;\;\;\;\;\; = {U^{(k)}} + \Delta t\left[ {\frac{{\partial {U^{(k)}}}}{{\partial x}}\frac{{\partial {U^{(k)}}}}{{\partial x}}\frac{{\partial {U^{(k)}}}}{{\partial z}}} \right] \times }\\ {\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{D}}{{\left[ {\frac{{\partial {U^{(k)}}}}{{\partial x}}\frac{{\partial {U^{(k)}}}}{{\partial y}}\frac{{\partial {U^{(k)}}}}{{\partial z}}} \right]}^{\rm{T}}}} \end{array} $

(11)

式中：U^(k)和U^(k+1)分别为第k和k+1次迭代计算得到的数据；Δt为迭代步长。在实际应用中为了保证迭代过程的收敛以及稳定性，Δt要足够小，但是若Δt过小，会增加迭代次数，从而增加计算时间。通常选择0.05≤Δt≤1即可满足地震数据处理要求。

本文方法的实现流程为：

(1) 定义初始迭代次数k为0，用原始三维数据U₀(x, y, z)作为迭代初始值并设置最大迭代次数为K；

(2) 使用噪声尺度为σ的高斯核函数对数据进行高斯滤波处理；

(3) 重构S_ρ并求取特征值λ₁、λ₂和λ₃以及对应的特征向量ζ₁、ζ₂和ζ₃；

(4) 利用M_fault设计D的特征值，并确定横向连续性因子ε；

(5) 构建D后，根据式(11)计算第k+1次迭代结果；

(6) 如果k＜K，转到步骤(2)继续迭代，否则结束迭代、输出结果。

2 理论模型测试

为验证基于M_fault的各向异性扩散滤波方法的应用效果，选用三维Qdome模型(图 1)^[26]进行验证。根据Qdome原始地震数据(图 2a)的频带对随机噪声进行低通滤波，然后将其加入到理论模型中得到加噪模型(图 2b)。

通过前文分析可知，各向异性扩散滤波涉及的参数较多，主要包括高斯核函数的噪声尺度σ、整合尺度ρ、迭代步长Δt以及最大迭代次数K。为了研究各个参数对滤波效果的影响，利用图 2b对参数选择进行试验。

2.1 滤波参数的选择 2.1.1 σ的影响

图 3展示了不同噪声尺度的滤波结果。由图可见：①在其他滤波参数固定的情况下，随σ取值增大，信噪比随之变大(图 3a、图 3b)，但当信噪比达到最大值(SNR=10.8947)后，随σ取值增大，信噪比开始下降(图 3c)；②σ取值过大，断层边界模糊(图 3d箭头处)。因此当σ∈[1.0，2.0]时滤波效果较好。

2.1.2 ρ的影响

图 4为不同整合尺度的滤波结果。由图可见：在其他滤波参数固定的情况下，随ρ取值增大，信噪比随之变大(图 4a、图 4b)，但当信噪比达到最大值(SNR=10.8947)后，开始出现不稳定变化(图 4c、图 4d)。因此，ρ取值不能太大。

2.1.3 Δt的影响

图 5为不同迭代步长的滤波结果。由图可见：Δt过小，导致计算时间增加；Δt过大，资料信噪比并未明显提高，而且造成地震数据的能量损失，断层边缘处出现模糊(图 5d红色箭头处)。因此，当Δt∈[0.05，1.00]时滤波效果较好，且计算效率较高。

2.1.4 K的影响

图 6为不同最大迭代次数的滤波结果。由图可见，K越大，计算时间越长，且会破坏断层的边界特征(图 6d红色箭头处)。因此，当K∈[20, 50]时，可明显提高信噪比，且断层边缘特征保持较好(图 6b、图 6c)。若地震资料信噪比太低，可以增大K，以取得较好的滤波效果。

2.2 模型试算

通过参数实验确定了滤波参数的选择范围，分别应用本文方法和各向异性扩散滤波方法^[15]对图 2b进行去噪处理，结果表明：本文方法滤波后数据(图 7a左)的信噪比高于各向异性扩散滤波方法(图 7b左)；本文方法滤波后地震数据基本没有能量损失(图 7a右)，各向异性扩散滤波方法滤波后地震数据在断层边缘处有轻微能量损失(图 7b右红色箭头处)。图 8为本文方法滤波前、后数据频谱。由图可见，滤波处理基本不影响地震数据频带。模型试算说明本文方法在去噪和保护地质体边缘信息方面的有效性。

3 实际地震资料应用

选取断裂众多、结构复杂、勘探难度较大的胜利油田X断块地震资料作为研究数据，应用文中方法进行滤波处理并计算相干属性。图 9为本文方法滤波前、后的Inline3570剖面。由图可见，滤波后地震剖面的同相轴更连续、断点更清晰(图 9b黄色和蓝色箭头处)，断层更易识别，说明本文方法在去噪的同时能有效保护断层等不连续性边界的信息。图 10为滤波前、后的相干时间切片。由图可见，利用滤波后的相干时间切片(图 10b)可准确地描述断层形态及展布，如在滤波前的相干时间切片中一些不易识别的小断层(图 10a中黄色和红色箭头处)，在滤波后的相干时间切片上得到准确、清晰地识别，这也证明了该滤波方法在提高资料信噪比和保持断层等边缘信息方面的有效性和和准确性。

4 结论

针对三维地震资料保边滤波处理的需要，本文研究了三维各向异性扩散滤波方法，并采用结构张量分析了三维地质结构的局部特征。基于断层置信度参数重新设计了扩散张量的特征值，可以控制地震数据沿不同方向的扩散强度，在各向异性扩散滤波处理中，能够更好地保护断层等边界信息。理论模型测试和实际地震资料应用证明了该方法在提高资料信噪比的同时，可有效保护断层等的边界信息，并在地层连续区域增强同相轴的连续性，为断层解释提供良好的基础数据。