0 引言

图像分析中的边缘检测技术被广泛用于地震数据裂缝检测。目前地震数据边缘检测方法主要有基于图像空间域微分算子(Prewitt算子、Robert算子、Canny算子、Sobel算子、高斯偏导滤波器(LOG)等)的边缘检测法及以变换域小波分析为基础的多尺度边缘检测法^[1]。实际勘探经验表明，油气藏中普遍分布孔隙和裂缝，因此可利用地震资料检测边缘特征。Wu等^[2-4]、丁燕等^[5]利用深度学习算法提取三维地震图像中的断层及其产状信息。

希尔伯特变换广泛应用于检测地震资料不连续边缘信息。党志敏等^[6]、熊晓军等^[7]利用广义希尔伯特变换检测含噪地震资料的边缘信息。陈学华等^[9-10]将高阶伪希尔伯特变换用于地震资料边缘检测。谢静等^[11]提出了时间域加窗希尔伯特变换边缘检测方法。基于三维地震资料在深度方向的多尺度特征，李斌等^[12]提出了多尺度加窗希尔伯特变换的地震资料体边缘检测方法。在前人的研究^[13-17]基础上，Lyu等^[18-19]实现了基于二维希尔伯特变换的三维地震资料体边缘检测方法。利用上述方法检测裂缝时，往往没有考虑地震资料不连续信息在方向上的多样化与复杂性。为此，周元茂等^[20]联合应用各向异性高斯迭代平滑滤波方法与地震体曲率分析方法，形成了组合型方向体曲率分析方法。

基于希尔伯特变换的边缘检测算法用于实际地震资料时，提取的往往是裂缝综合信息，信息繁杂且会弱化部分有效信息，难以分辨某一特定方向的数据特点，不利于分析裂缝走向。因此，考虑到断层、裂缝等不连续信息的复杂性，为了分析地震资料中各个方向的不连续特征，本文构建了一种任意角度旋转的加窗希尔伯特变换(arbitrarily rotated windowed Hilbert transform，ARWHT)边缘检测算子，将其用于三维地震资料体边缘检测，不仅能提取地震资料在任意方向的不连续特征，还可突出局部异常信息，减少噪声的影响，从而获得更精确的储层裂缝走向及分布信息。

1 方法原理 1.1 ARWHT边缘检测算子

为了充分提取、分析三维地震资料中任意方向的不连续信息并突出局部化特征，本文构建了ARWHT边缘检测算子。

若将希尔伯特算子所有点表示为一个几何信息矩阵Q，将其与方向检测模板T相乘﹐就得到旋转变换后的几何信息矩阵Q′，即

$ {\mathit{\boldsymbol{Q}}}\prime = {\mathit{\boldsymbol{QT}}} $

(1)

其中

$ {\mathit{\boldsymbol{T}}} = \;\left( \begin{array}{l} \;\;\;\;{\rm{cos}}\theta \;{\rm{sin}}\theta \\ \;\; - {\rm{sin}}\theta \;{\rm{cos}}\theta \end{array} \right) $

设加窗后的希尔伯特算子序列的样点坐标为P[n，h′(n)]，将其逆时针旋转角度θ后，新算子序列的样点坐标为P′[n_R，h′_R(n)]，则有

$ \left\{ \begin{array}{l} {n_{\rm{R}}} = an + ch\prime \left( n \right)\\ {{h'}_{\rm{R}}}\left( n \right) = bn + dh\prime \left( n \right) \end{array} \right. $

(2)

其中

$ h\prime \left( n \right) = h\left( n \right) \cdot w\left( n \right) $

(3)

式中：h(n)为原始希尔伯特变换算子，n为采样点号；w(n)为对算子作局部化处理的窗函数；a、b、c、d为设置的方向检测模板初始形式的组成参数。将式(2)改为矩阵形式

$ \left[ {{n_{\rm{R}}}\;h{\prime _{\rm{R}}}\left( n \right)} \right] = \left[ {n\;h\prime \left( n \right)} \right]\left( \begin{array}{l} a\;b\\ c\;d \end{array} \right) $

(4)

设点P与坐标原点O的连线OP与x轴的夹角为α，则有

$ \left\{ \begin{array}{l} h\prime \left( n \right) = \overline {OP} \cdot{\rm{sin}}\alpha \\ n = \overline {OP} \cdot{\rm{cos}}\alpha \end{array} \right. $

(5)

$ \left\{ \begin{array}{l} {n_{\rm{R}}} = \overline {OP\prime } {\rm{cos}}\left( {\alpha + \theta } \right) = n{\rm{cos}}\theta - h\prime \left( n \right){\rm{sin}}\theta \\ h{\prime _{\rm{R}}}\left( n \right) = \overline {OP\prime } {\rm{sin}}\left( {\alpha + \theta } \right) = n{\rm{sin}}\theta + h\prime \left( n \right){\rm{cos}}\theta \end{array} \right. $

(6)

式中 OP、OP′分别为OP、OP′的长度，则式(4)变为

$ \left[ {{n_{\rm{R}}}\;h{\prime _{\rm{R}}}\left( n \right)} \right] = \left[ {n\;h\prime \left( n \right)} \right]\left( \begin{array}{l} \;\;{\rm{cos}}\theta \;{\rm{sin}}\theta \\ - {\rm{sin}}\theta \;{\rm{cos}}\theta \end{array} \right) $

(7)

为强调地震资料的局部化特征、突出不连续信息，本文选择高斯窗对希尔伯特变换算子进行局部化处理，使检测对象的有效边缘信息更突出、非边缘信息的差异更明显。高斯窗函数的表达式为

$ w\left( n \right) = \frac{1}{{\sigma \sqrt {2\pi } }}{{\rm{e}}^{ - \frac{{{n^2}}}{{2{\sigma ^2}}}}} $

(8)

式中σ为决定窗函数时间宽度的尺度因子，即σ越大，窗函数越宽。图 1为ARWHT算子。

将ARWHT算子用于边缘检测，得到一个新的基于ARWHT的边缘提取公式，边缘提取结果为

$ \begin{array}{l} \tilde x\left( {{n_{\rm{R}}}} \right) = x\left( {{n_R}} \right) * h{\prime _{\rm{R}}}\left( n \right)\\ \;\;\;\;\;\;\;\;\; = \sum\limits_{k = 0}^{N - 1} x \left[ {{n_{\rm{R}}} - \left( {2k - 1} \right)} \right] \cdot h{\prime _{\rm{R}}}\left( {2k - 1} \right) \end{array} $

(9)

其中

$ {n_{\rm{R}}} = n{\rm{cos}}\theta - h\prime \left( n \right){\rm{sin}}\theta $

(10)

$ h{\prime _{\rm{R}}}\left( n \right) = n{\rm{sin}}\theta + h\prime \left( n \right){\rm{cos}}\theta $

(11)

式中N为x(n_R)的长度。

1.2 基于ARWHT的地震资料体边缘检测算法的实现

由式(9)~式(11)计算任意特定方向的边缘提取结果，进而由

$ D\left( {{m_0}, k, \theta } \right) = \sum\limits_{\lambda = - k}^k {\frac{{\left[ {x\left( {{m_0} + \lambda \cdot\Delta m, {n_\theta }} \right) * {{h'}_\theta }\left( n \right)} \right] + \left| {{{\tilde x}_{\theta , \lambda - \min }}} \right|}}{{{{\tilde x}_{\theta , \lambda - {\rm{max}}}} - {{\tilde x}_{\theta , \lambda - {\rm{min}}}}}}} $

(12)

计算基于ARWHT的地震资料体边缘检测结果。式中：D(m₀，k，θ)表示算子旋转角度为θ时的地震资料体边缘检测结果，m₀为目的层的时间位置，k为三维地震资料在m₀上、下的采样点数，则计算的地震资料层段厚度为2k·Δm(Δm为在深度时窗上相邻样点之间的采样间隔)；x(m₀+λ·Δm，n_θ)为沿层地震数据；$ {{\tilde x}_{\theta , \lambda - {\rm{min}}}}、{{\tilde x}_{\theta , \lambda - {\rm{max}}}}$分别为第λ个切片裂缝提取结果的最小值、最大值。由希尔伯特变换边缘检测原理可知，当算子旋转角度为θ时，可提取θ±90°方向的边缘信息。

对于所有的算子旋转角度θ∈[0°，180°]，取每个空间位置的D(m₀，k，θ)最大值，由

$ D'\left( {{m_0}, k} \right) = \mathop {{\rm{max}}}\limits_\theta \left[ {D\left( {{m_0}, k, \theta } \right)} \right] $

(13)

得到融合了所有方向的基于ARWHT的三维地震资料体边缘检测结果。

2 模型试算与分析

为了验证本文算法的效果，设置一个正八边形模型(图 2a)，利用ARWHT算子提取各方向的边缘信息(图 2b、图 2c)。结果表明：逆时针旋转45°的ARWHT算子检测结果(图 2b)着重突出了0°、90°、135°的边缘信息(135°边缘信息绝对值最大)，有效过滤了45°的边缘信息；逆时针旋转135°的ARWHT算子检测结果(图 2c)着重突出了0°、45°、90°的边缘信息(45°边缘信息绝对值最大)，有效过滤了135°边缘信息。

为了验证本文方法提取复杂边缘特征信息的有效性，对经典的Lena图像提取边缘信息(图 3)。结果表明：①原始图像(图 3a)背景中含有较多不同角度的线条，有利于验证算法的有效性。②本文算法弱化了与算子旋转角度一致的边缘信息，同时突显了与其垂直的边缘信息(图 3b、图 3c、图 3e、图 3f)，且边缘信息提取效果优于常规希尔伯特变换(图 3d)。如0°算子对帽子的阴影部分表现更好(图 3b红色长方形框处)，90°算子提取的帽檐上边缘明显更连续(图 3c黄色箭头处)，且唇部线条也更突出。③由45°(图 3e)和135°(图 3f)算子检测结果可见，对于背景中不同方向的边缘信息，不同旋转角度算子侧重提取不同方向的边缘信息(绿色和蓝色箭头处)。④由于眉毛角度的不同，不同旋转角度的算子检测结果对眉毛的表现效果也有差异(红色椭圆框处)。

3 实际地震资料处理

为了说明本文方法在实际地震资料处理中的应用效果，以中国KL地区的三维地震资料为例，分别使用旋转0°、45°、90°、135°的算子提取体边缘信息。由KL地区沿目的层振幅切片(图 4)可见，目的层内河道和断层分布极为复杂。

基于ARWHT的实际地震资料体边缘检测结果(图 5)表明：①不旋转(图 5a)与逆时针旋转90°(图 5b)算子的方向相互垂直，检测结果差别明显，有利于刻画细小裂缝的方向。表现为：前者着重突出了纵向构造特征，断层及裂隙连续性明显加强，显著减弱了横向构造信息；后者的横向不连续性信息以及构造异常特征明显加强，明显压制了纵向边缘信息，显示出原始地震资料中难以用肉眼分辨的大型横向连续裂缝(红色长方形虚线框处)。另外，还可以看到一条隐蔽的横向河道，在其他提取结果中(图 5a、图 5c、图 5d)难以分辨(蓝色箭头处)。②逆时针旋转45°算子检测结果(图 5c)突显了135°的体边缘信息(黄色箭头处)；逆时针旋转135°算子检测结果(图 5d)着重表现了45°的不连续信息(绿色箭头处)，还可以追踪到一条走向大致呈45°方向的蜿蜒河道(绿色点状线所示)，在其他提取结果(图 5a、图 5b、图 5c)中难以分辨。③综合四个方向的提取结果来看，边缘检测结果的不同区域(黄色方形虚线框和蓝色椭圆虚线框)数据的特征不同，包含的不连续信息方向特征与ARWHT算子提取边缘的规律一致，更好地刻画了地震异常信息的空间分布和走向特征。

在利用本文方法处理实际地震资料时，应考虑数据的具体情况选择算子的旋转角度。对于已知走向信息的大断裂或者某一区域内分布方向大致相同的裂缝系统，应选择与其走向或分布方向垂直的算子旋转角度提取边缘信息，从而加强走向和分布方向的不连续信息特征；对于复杂多变的裂缝系统，其包含各个方向的信息，此时应选择多个旋转角度提取边缘信息进行融合，以显示更丰富的有效信息，利于地震资料精细解释。

为了进一步分析本文方法的应用效果及优势，根据式(13)融合图 5中各方向的边缘信息提取结果(图 6c)，并与常规希尔伯特变换边缘检测结果(图 6a)、RMS振幅体切片(图 6b)对比、分析。结果表明：图 6c综合了图 5各结果的优点，包含了各个方向的体边缘特征信息(图中红色虚线框与箭头处)；与图 6a和图 6b相比，图 6c的有效信息更全面，不仅清楚地展示了河道和断裂系统的分布特征，还刻画了裂缝不连续地质异常的边界信息，有利于地震资料的精细解释。

为了更清晰地展示河道、断层、断裂带等信息，将本文方法的裂缝检测结果(图 6c)与图 6b融合，融合结果(图 6d)不仅展示了多方向的裂缝信息，还刻画了河道走向、裂缝边界信息，展示了丰富的地质信息。

综上所述，本文方法能突显地震资料中任意特定方向的地质异常信息，更全面地展示了地震不连续信息的分布特征，检测结果体现了算法的可行性和优越性。

4 结束语

本文在加窗希尔伯特变换理论的基础上，构建了ARWHT边缘检测算子，将其应用于实际地震资料边缘检测，可以得到任意特定方向的不连续异常特征，实现了提取各方向不连续特征的三维地震资料体边缘检测。

地震资料的ARWHT体边缘检测方法在表征储层和地质结构边缘信息时保留了有效的方向特征。因此，可以在实际地震数据中充分挖掘并突显任意特定方向的地质特征。对于具有不同特征的数据，应选择不同旋转角度的算子进行边缘检测。实际地震数据测试表明，ARWHT体边缘检测方法可描述河道、断层、断裂带、构造、储层裂缝分布，强化地震数据在任意特定方向的构造信息，更好地突出地质构造特征，为解释构造、了解不连续地质异常信息及其展布规律、预测储层等提供依据。