随着油气田的持续开发, 寻找复杂油气藏和隐蔽油气藏成为当下油气勘探的重点。在三维地震数据处理中, 充分利用数据中的体信息, 能够有效提高地震资料的信噪比, 改善地震数据中地质体边界的检测精度^[1-2]。传统的压制随机噪声的方法可分为空间域和变换域方法, 前者主要包括均值滤波、中值滤波^[3-4]及各向异性扩散滤波^[5]等方法, 后者主要包括傅里叶变换、Radon变换、小波变换、曲波变换阈值去噪方法^[6-8]等。地震资料噪声压制的效果会直接影响后续的处理与解释, 而传统去噪方法应用于三维地震数据处理往往存在去噪效果差、弱有效信号被压制、伪影干扰等缺点^[9], 如曲波变换能够压制地震数据中的随机噪声, 但它有一般变换域方法所固有的一个缺点, 即对复杂地震数据的适应性较差, 处理后的地震数据中往往会引入假信号^[10]。为了改善传统去噪方法应用于三维地震数据时的处理效果, 张华等^[11]研究了多尺度、多方向性的曲波变换三维地震数据去噪方法, 实现了二维曲波变换对三维地震数据的去噪处理, 相较于三维曲波变换, 提高了运算效率^[12]。随着计算机存储能力及计算水平的不断提高, 针对三维地震数据的体处理方式正逐步取代传统按剖面处理的方式, 可有效避免二维处理方法对三维数据体结构特征的破坏^[12]。三维地震数据中存在较强的相似性信息, 张岩等^[13]研究了利用多道相似组稀疏表示方法压制随机噪声, 能更好地保持复杂区域地震数据同相轴的局部特征。

充分利用地震数据中的地质体边界信息, 能够帮助处理和解释人员更准确地开展地震资料的处理和解释工作。边缘检测最早于1959年提出, 之后为了追求更好的检测效果, 各种检测方法也陆续出现^[14]。常规的边缘检测方法主要包括微分算子法、最优算子法、多尺度方法、数学形态学方法及神经网络方法等^[15]。Canny边缘检测方法于1986年首次提出^[16], CANNY对一些经典边缘检测方法和应用进行了总结, 在此基础上提出了边缘检测准则, 采用最优化数值方法, 得到一个实用的边缘检测算法。

为了充分利用地震数据中的自相似信息和冗余信息, 在实现去噪的同时能保持地震数据的细节特征, 突出地质体边界信息, 本文研究了针对三维地震数据的自相似块匹配去噪和边缘检测方法, 该去噪方法是一种非局部变换域去噪方法, 充分利用了三维地震数据的自相似性和冗余性信息, 理论上去噪结果可保留更多的有效信号, 去除多余的噪声。采用自相似块匹配Canny边缘检测算法, 对去噪后数据的时间切片进行边缘检测, 理论上使用自相似块匹配去噪后的数据可以保留更多边界信息, 使得边缘检测的效果更好。将该方法应用于实际三维地震资料, 并与传统的基于三维曲波变换的去噪方法相比较, 验证了方法的有效性。

1 自相似块匹配基本理论

自相似块匹配去噪算法是一种利用地震数据自身相似性的去噪方法, 首先通过块匹配将三维地震数据中的相似块聚集, 继而在变换域中衰减噪声, 最后通过逆变换得到去噪后的剖面^[17-18]。该方法能够保留原始数据中断层边缘、裂缝等区域的不连续性信息, 同时又能保持有效的去噪效果^[19-20]。

1.1 自相似块匹配信噪分离

自相似块匹配算法针对含噪数据z, 实现对无噪数据y和噪声数据η的分离^[21]。算法主要包括两个阶段, 即阈值滤波阶段和维纳滤波阶段, 每个阶段又包括分组、联合滤波和聚集3个步骤。

1.1.1 阈值滤波阶段

1) 分组。处理三维地震数据时, 首先按照指定的大小将三维地震数据分成大小相等的数据块, 从中选出一块作为参考块, 其它块为候选块, 从候选块中找出与参考块相似的块, 将它们进行组合。

令C_xR^z表示大小为L×L×L的数据块, 在阈值阶段, 从含噪数据z中获得三维坐标x_R∈X。两个数据块之间的相似性表示为:

$ \begin{array}{c} d\left(C_{x_{i}}^{z}, C_{x_{j}}^{z}\right)=\frac{\left\|C_{x_{i}}^{z}-C_{x_{j}}^{z}\right\|_{2}^{2}}{L^{3}} \\ i, j \in R \end{array} $

(1)

式中:‖·‖₂²表示两个数据块之间的欧式距离; 分母L³为标准化因子; 数据块的大小L根据实际数据进行选定, 为充分保留数据的细节特征, 其取值不宜过大, 一般在10以内即可^[22]。如果两个数据块的相似性小于或等于预设阈值τ, 则认为它们是相似的。

2) 联合滤波。对相似块组成的四维数组进行联合滤波, 去除数据中的随机噪声^[23]。

首先利用离散小波变换对每个数组ASxRzz进行四维T4Dt变换, 然后采用阈值算子Υ对其谱系数进行收缩去噪, 去噪后的数组Υ[T4Dt(ASxRzz)]表示为$\hat A^y_{x_R}$, 再将其进行逆T4Dt变换:

$ T_{4 \rm D}^{t-1}\left\{\varUpsilon \left[T_{4 \mathrm{D}}^{t}\left(A_{s_{x_{R}}^{z}}^{z}\right)\right]\right\}= \hat A_{s_{x_{R}}^{z}}^{y}=\prod\limits_{x_{i} \in S_{x_{R}}^{z}} \hat{C}_{x_{i}}^{y} $

(2)

公式(2)中的每个数据组$\hat C^y_{x_i}$是对无噪数据y中的数据块C_xi^y的估计; S_xR^z为从含噪数据z中提取的与C_xR^z相似的数据块集合。

3) 聚集。将逆变换处理后的数据块还原到组合前的位置, 实现对数据的重构。得到的去噪数据组中, 单个区域可能分布在不同的数据组中, 导致同一区域的数据块可能会重叠, 因此, 在重叠的区域内不同的数据块为同一位置提供了多个互不相同的数值估计^[24]。在聚集步骤中, 对重叠部分进行加权平均, 利用这种重叠后的冗余信息得到去噪后的硬阈值估计结果:

$ \hat{y}^{t}=\frac{\sum_{x_{R} \in X}\left(\sum_{x_{i} \in S_{x_{R}}^{z}} \omega_{x_{R}}^{t} \hat{C}_{x_{i}}^{y}\right)}{\sum_{x_{R} \in X}\left(\sum_{x_{i} \in S_{x_{R}}^{z}} \omega_{x_{R}}^{t} \chi_{x_{i}}\right)} $

(3)

式中:ω_xR^t是与数据组相关的权系数; χ_xi是$\hat C^y_{x_i}$域的指标函数, 在$\hat C^y_{x_i}$区域内χ_xi=1, 在区域外χ_xi=0。(3)式中的聚集权系数定义为:

$ \omega_{x_{R}}^{t}=\frac{1}{\sigma^{2} N_{x_{R}}^{t}} $

(4)

式中:σ是观测数据z中噪声的标准差; NxRt表示(2)式Υ[T4Dt(ASxRzz)]中非零系数的个数。表现出高度相关性的数据块具有较大的权系数, 而其它具有噪声的数据块用小权系数进行压制。

1.1.2 维纳滤波阶段

维纳滤波是对阈值估计结果$\hat y^t$进行处理, 其噪声含量较低, 可以对其进行更有效的稀疏化表示, 从而获得更好的去噪结果^[25], 因此对x_R∈X内的每个参考数据块$C^{\hat y^t}_{x_R}$建立相似数据块集。

联合滤波采用经验维纳滤波器完成, 首先建立第一个相似的四维数据块集$A^{\hat y^t}_{S^{\hat y^t}_{x_R}}$, 然后根据其频谱的能量, 定义经验维纳滤波器收缩系数为:

$ W_{S^{\hat y^t}_{x_R}}= \frac{{|T^{\rm wie}_{\rm 4D}(A^{\hat y^t}_{S^{\hat y^t}_{x_R}})|^2}}{{|T^{\rm wie}_{\rm 4D}(A^{\hat y^t}_{S^{\hat y^t}_{x_R}})|^2+σ^2}} $

(5)

式中:T_4D^wie表示四维线性变换。对含噪数据z中相同条件建立的第二个相似块集$A^{z}_{S^{\hat y^t}_{x_R}}$的频谱与收缩系数进行逐元素相乘实现协同滤波, 采用逆T_4D^wie变换得到四维数据体估计:

$ \hat{A}_{S_{{{x}_{R}}}^{{{{\hat{y}}}^{t}}}}^{y}=T_{4\text{D}}^{\text{wie}{{\text{ }}^{-1}}}\left[ {{W}_{S_{{{x}_{R}}}^{{{{\hat{y}}}^{t}}}}}\cdot T_{4\text{D}}^{\text{wie }}\left( A_{S_{{{x}_{R}}}^{{{{\hat{y}}}^{t}}}}^{z} \right) \right] $

(6)

在四维变换之后, 对于每个x_R∈X都能计算出他们的权重值:

$ \omega_{x_{R}}^{\text {wie }}=\sigma^{-2}\left\|W_{S_{x_{R}}^{\hat{y} t}}\right\|_{2}^{-2} $

(7)

再把所有多维相似块集的基本估计结果聚集就能得到最终的维纳滤波估计:

$ \hat{y}^{\text {wie }}=\frac{\sum_{x_{R} \in X}\left(\sum_{x_{i} \in {S}_{x_{R}}^{z}} \omega_{x_{R}}^{\text {wie }} \hat{A}_{S_{x_{R}}^{\hat y^t}}^{y}\right)}{\sum_{x_{R} \in X}\left(\sum_{x_{i} \in S_{x_{R}}^{z}} \omega_{x_{R}}^{\text {wie }} \chi_{x_{i}}\right)} $

(8)

该方法可以利用每个多维数据集中数据块的局部相关性, 以及不同数据集中同一数据块的非局部相关性, 分离有效信号和噪声^[26]。基于自相似块匹配去噪的技术流程如图 1所示。

1.2 自相似块匹配Canny边缘检测

边缘检测算法通常具有良好的边缘检测性能和定位性能以及单一边界响应特性^[27]。为了满足上述特性, 本文采用自相似块匹配去噪处理代替常规Canny边界算法中高斯平滑滤波处理, 从而增强对有效信号的保护和检测效果。

自相似块匹配Canny边缘检测算法包括以下4个步骤。

1) 对数据进行去噪处理。使用前述自相似块匹配去噪方法对数据进行去噪处理, 尽可能减少噪声对边缘检测结果的影响。

2) 计算梯度的幅值和方向。采用2×2邻域一阶偏导的有限差分方法计算去噪后数据的梯度方向和梯度幅值, 再由此确定数据点的梯度强度G和方向θ:

$ \begin{aligned} &G=\sqrt{G_{x}^{2}+G_{y}^{2}}\\ &\theta=\arctan \left(G_{y} / G_{x}\right) \end{aligned} $

(9)

式中:G_x和G_y分别为去噪后数据水平和垂直方向的一阶导数值。

3) 压制梯度幅值。非极大值压制是一种边界稀疏技术, 将数据点的梯度幅值与沿正负梯度方向上两个数据点的幅值进行比较, 当该数据点的梯度幅值大于另外两个数据点的梯度幅值时, 则认为该点为边界点。通常为了更加精确的计算, 在跨越梯度方向的两个相邻数据点之间使用线性插值来得到要比较的数据点的梯度幅值。

4) 双阈值检测和边界连接。设置两个阈值τ_h和τ_l, 其中, τ_h为高阈值, τ_l为低阈值, 若数据点处的梯度幅值大于τ_h, 则认为该点为边界点; 若数据点处的梯度幅值小于τ_l, 则认为该点不是边界点; 若梯度幅值介于两者之间, 则认为该点为弱边界点, 可能是真实边界或者是由噪声干扰造成的。通常, 由真实边界引起的弱边界数据点连接到强边界数据点, 而噪声干扰引起的弱边界点不连接。为了跟踪边界连接, 通过查看弱边界数据点及其邻域数据点, 只要邻域数据点其中一个为强边界数据点, 则该弱边界数据点就可以保留为真实边界。

图 2a和图 2b分别给出了对同一含噪Marmousi合成地震记录采用常规Canny边缘检测方法和自相似块匹配Canny边缘检测方法得到的检测结果。图 2中黑色线条为边缘检测结果。本模型数据浅层和部分深层的有效信号能量较弱, 且受到噪声干扰。由图 2a和图 2b可见, 常规Canny边缘检测方法采用高斯平滑滤波, 该方法只能检测到部分弱信号; 而自相似块匹配Canny边缘检测方法能对这些弱信号进行有效检测, 边缘检测结果准确细致, 能够清晰地分辨边缘位置。

2 理论模型试算 2.1 二维切片去噪与自相似块匹配去噪结果对比

为了测试自相似块匹配去噪算法的去噪效果以及自相似块匹配Canny边缘检测算法对时间切片的边缘检测效果, 对含噪合成地震数据分别进行了传统二维切片去噪和自相似块匹配去噪处理, 并对其去噪后的自相似块匹配Canny边缘检测结果进行了比较(图 3)。图 3a为去噪前的单炮地震数据, 由图 3a可见, 随机噪声存在于三维数据空间; 图 3b为采用二维切片去噪后的单炮地震数据, 由图 3b可见, 垂直切片方向依然存在随机噪声; 图 3c为采用自相似块匹配方法得到的去噪结果, 由图 3c可见, 存在于三维数据空间的噪声得到了消除, 各个方向上均无噪声残留; 图 3d为无噪数据时间切片的自相似块匹配Canny边缘检测结果; 图 3e为对应的含噪数据时间切片自相似块匹配Canny边缘检测结果, 图 3d和图 3e中红色线条是检测的边界, 对比图 3d和图 3e可以看出, 在存在噪声的情况下, 本文的自相似块匹配Canny边缘检测算法也能准确识别出信号边界。

图 4给出了低信噪比条件下采用自相似块匹配方法得到的去噪结果。图 4a为含强随机噪声的地震记录, 该地震记录中加入了高斯白噪声, 信噪比为4.51 dB; 图 4b为采用自相似块匹配方法去噪后的地震记录; 图 4c为去除的噪声。由图 4b和图 4c可以看出, 噪声均得到有效压制, 残差剖面中无有效信号残留。由此可见, 对于含高比例噪声的情况, 本文提出的自相似块匹配方法也能够得到有效的去噪结果。

为了对比本文方法与曲波变换阈值法的去噪效果, 我们对图 5a所示含噪地震记录分别采用两种方法进行去噪处理, 得到的结果如图 5b和图 5c所示。由图 5b和图 5c可见, 采用曲波变换阈值法去噪后的剖面中出现了伪影现象, 而且在去噪的同时损害了部分有效信号; 而采用本文方法去噪后的剖面上没有残留随机噪声, 有效信号也得到了很好的保护, 没有产生多余的虚假信息, 说明本文方法具有更好的保真度。

边缘保持指数(EPI)^[28]用来定量评价边缘保持程度, 其计算公式为:

$ \mathrm{EPI}=\frac{\sum\left[\left|I_{\rm r}(i, j)-I_{\rm r}(i+1, j)\right|+\left|I_{\rm r}(i, j)-I_{\rm r}(i, j+1)\right|\right]}{\sum\left[\left|I_{\rm o}(i, j)-I_{\rm o}(i+1, j)\right|+\left|I_{\rm o}(i, j)-I_{\rm o}(i, j+1)\right|\right]} $

(10)

式中:I_r表示去噪处理后的数据; I_o代表原始数据。边缘保持指数越大, 说明数据中边缘和细节信息保持得越好。不同信噪比数据的EPI如图 6所示。由图 6可以看出, 自相似块匹配方法对数据的边缘保持效果更好。综上所述, 自相似块匹配方法无论是在提高信噪比还是保真度方面, 其效果都优于传统的曲波变换方法。

2.2 自相似块匹配边界保持效果分析

建立如图 7a所示的含断层的三维地震数据模型, 在合成地震记录中加入随机噪声, 如图 7b所示, 分别采用FXdecon方法和自相似块匹配去噪方法对含噪地震记录进行去噪处理, 结果如图 7c和图 7d所示。由图 7c和图 7d可以看出, 两种方法均可以压制原始数据中的随机噪声, 但FXdecon方法去噪结果中, 断层边缘处同相轴发生畸变(如图中箭头所示), 而自相似块匹配方法对地震数据边缘处的不连续性信息保持效果较好, 能够为边缘检测提供稳定的地震数据。

3 实际资料处理分析

对某地区三维叠后地震数据进行了自相似块匹配去噪处理。图 8a为去噪前的地震数据, 图 8b为采用自相似块匹配方法得到的去噪后的地震数据, 图 8c为去除的随机噪声。由图 8可见, 去噪后的剖面信噪比得到有效提高, 尤其是复杂构造处, 剖面品质得到改善, 地质体边界的清晰度得到增强, 残差剖面上没有有效信号。说明本文提出的自相似块匹配去噪方法可以在去噪的同时保护有效信号不被破坏。

图 9a为图 8a所示地震数据的一个等时切片; 图 9b为采用常规Canny边缘检测方法对时间切片的边缘检测结果; 图 9c为采用自相似块匹配Canny边缘检测方法对时间切片的边缘检测结果。图 9b和图 9c中红色线条为边缘检测结果。对比图 9b和图 9c可以看出, 本文方法能够有效保持地震信号的边界, 即使是弱边界也保持得很好, 检测结果与实际切片中的边界位置吻合。

4 结论

本文研究了针对三维地震数据的自相似块匹配去噪和边缘检测方法。自相似块匹配方法将每个三维数据块分配到多个组中, 进行多次协同估计, 可以充分利用原始数据的自相似性和冗余性信息。对自相似块匹配去噪处理后的数据进行Canny边缘检测, 能够有效检测断裂、裂缝等边缘特征。理论模型及实际资料的应用结果表明:

1) 自相似块匹配去噪方法可以有效去除地震资料中的随机噪声, 与传统曲波变换阈值法相比, 本文方法去噪后的地震数据信噪比及保真度得到提升, 去噪结果中无伪影现象。去噪后得到的高信噪比、高保真度的地震数据能够为地震资料的后期处理和解释提供高质量和准确的地下信息。

2) 在边界保持方面, 自相似块匹配去噪方法的边界保持能力较好, 断层等同相轴不连续处在去噪处理后不会产生畸变, 断面清晰, 能够为边缘检测提供稳定、准确的去噪数据, 为断层等不连续处的边缘检测等提供了技术支撑。