地震几何属性的快速算法实现.pdf

书书书第 卷 第期 年月地球物理学报 ，张文，伍新明，漆杰 地震几何属性的快速算法实现地球物理学报，（）：，：，犆 犺 犻 狀 犲 狊 犲 犑犌 犲 狅 狆 犺 狔 狊（），（）：，：地震几何属性的快速算法实现张文，伍新明，漆杰中国科学技术大学地球和空间科学学院，合肥 ，休斯顿 摘要地震几何属性被广泛用于提取地震数据中的几何结构特征，从而辅助解释相关的地质沉积和构造过程本文提出基于多种递归滤波和构造导向滤波的地震几何属性快速算法，能显著提高地震相干和体曲率的计算效率和分辨率递归滤波的计算效率远高于传统的加权求和计算，且其计算成本与平均窗口的大小无关；同时，利用高维滤波的可分离特性可以将其分解为多个一维滤波来进一步提升计算效率并有利于多线程并行运算此外，使用构造导向滤波实现相干计算，可以有效地消除倾斜构造对结果的影响而无需传统方法中耗时的插值运算三维实际地震数据的应用结果表明，本文提出的快速算法能将传统的相干和曲率属性计算速度提高 倍且对断层的刻画更加完整和连续关键词地震属性；相干属性；曲率；递归滤波器 ：中图分类号 收稿日期 ，收修定稿基金项目国家自然科学基金“地震资料中断层的自动化智能解释”（）资助第一作者简介张文，女，年生，在读博士研究生，主要从事特征工程等方面的研究 ：通讯作者伍新明，男，年生，教授，主要从事人工智能地球物理等方面的研究 ：犉 犪 狊 狋 犮 狅 犿 狆 狌 狋 犪 狋 犻 狅 狀狅 犳 狊 犲 犻 狊 犿 犻 犮犵 犲 狅 犿 犲 狋 狉 犻 犮犪 狋 狋 狉 犻 犫 狌 狋 犲 狊 ，犛 犮 犺 狅 狅 犾 狅 犳犈 犪 狉 狋 犺犪 狀 犱犛 狆 犪 犮 犲犛 犮 犻 犲 狀 犮 犲 狊，犝 狀 犻 狏 犲 狉 狊 犻 狋 狔狅 犳犛 犮 犻 犲 狀 犮 犲犪 狀 犱犜 犲 犮 犺 狀 狅 犾 狅 犵 狔狅 犳犆 犺 犻 狀 犪，犎 犲 犳 犲 犻 ，犆 犺 犻 狀 犪犌 犲 狅 狆 犺 狔 狊 犻 犮 犪 犾 犐 狀 狊 犻 犵 犺 狋 狊，犎 狅 狌 狊 狋 狅 狀 ，犜犡，犝 犛 犃犃 犫 狊 狋 狉 犪 犮 狋 ，期张文等：地震几何属性的快速算法实现 ，犓 犲 狔 狑 狅 狉 犱 狊 ；引言地震属性是通过数值计算来量化地震数据的振幅和 形 态 特 征，进 而 提 取 地 质 特 征 的 有 效 手 段（，）地震属性中的几何属性主要用来解释构造和地层特征，借助几何属性，地质解释工作者能够全面而有效地解释不连续 特 征、层 序 边 界、沉 积 环 境 和 层 理 模 式 等（，）倾角和方位角、相干属性、曲率属性以及其他地震振幅或能量横向变化的体积估计统称为几何属性（，）其中，相干性和曲率可用于刻画地震数据中的地质构 造 和 地 层 特 征（，）自 年提出以来，地震相干属性已被广泛用于提取河道、裂缝、岩溶塌陷、断层等地质构造特征（杨金政，；，；，；，；，；，；白博等，；，；，）为了提高结果的分辨率和稳定性，目前已经提出了三代相干算法：经典的归一化互相关（）、多道相似性（）和基于特征结构的相干（）（，；，；，）虽然 算法的计算效率高并且公式简单易于实现，但是该方法对数据中的噪声、计算窗口的选择敏感，且结果只局限于目标地震道与其测线方向上的相邻地震道之间的相似性评估，没有消除地震同相轴倾角的影响 算法的计算成本高于，但是 算法的抗噪能力强，而且可以通过同相轴倾角校正来得到更合理的结果相比前两种算法，算法得到的结果分辨率最高，并且算法本身已经消除了构造倾角的影响，但是 的计算成本非常高，在处理大体量地震数据时耗时长因此，近些年很多研究关注于提高相干算法的计算效率（，；，；，）另一个重要的几何属性是曲率，它能帮助解释人员识别构造和地层特征（，）构造曲率属性的计算最初是基于手动选择的层位，这就不可避免地会受到解释人员层位识别能力的影响，尤其是面对被噪声污染的数据集时，识别层位的工作就会变得更加困难（，；，）直到 和 （）、和 （）提出了体曲率的计算方法，这个问题才得以解决体曲率是对计算窗口内同相轴倾角的一阶导数进行一系列计算得到的，并不需要层位信息因此，体曲率的提出为曲率属性的计算提供了极大的便利，并且提高了计算结果的客观性和准确性除此之外，对于同一个目标地质体，通过调整曲率计算中的波长，还可以描绘不同尺度的地质特征（，；，）一般来说，短波长曲率能够刻画局部裂缝发育带的细节特征，而长波长曲率可以凸显大范围的弯曲特征，这些特征通常与低于地震分辨率的断裂带、崩塌和成岩作用相关（，；，）在众多体曲率属性中，最正曲率和最负曲率是最有用的，因为它们公式简单、计算简便，并且可以描绘出常见的地质结构（，；张宏伟等，）目前对曲率的体积计算方法主要是基于对同相轴倾角的估计，因此许多研究着重关注于通过改进倾角估计方法来提高体曲率的分辨率（，；，；，；，；，；李海山和杨午阳，）然而这些方法的计算成本通常很高，更高效的体 曲 率 实 现 方法 有待 进 一 步研 究（，）针对目前相干和曲率属性计算效率较低的问题，本文提出使用递归滤波器和构造导向滤波器提高相干和曲率属性的计算效率、改善计算结果本文主要探究将递归平滑滤波器和构造导向的平滑滤波器应用于 相干算法时的计算效率和结果，以及提出应用递归导数滤波器计算体曲率来提高曲率属地 球 物 理 学 报（）卷性的计算效率和结果首先，使用简单的一维测试模型介绍了几种平滑滤波器的原理和特性其次，介绍了高斯滤波器和指数滤波器的递归算法，通过理论推导和简单二维模型测试，证明了高斯滤波器和指数滤波器的可分离性，并通过计算证明了利用可分离性可以实现多维数组的滤波计算提速除此之外，可分离性使得递归滤波器的并行计算在多维数据中的应用成为可能带有丰富方向信息的地震图像适合使用构造导向的平滑滤波，利用简单的二维合成地震数据的方向信息进行构造导向平滑滤波，验证了将构造导向平滑滤波器应用于地震数据的可行性以梵高的著名作品 星夜 为例，对比了递归平滑滤波器和构造导向滤波器的计算效率和平滑效果基于以上对于递归平滑滤波器和构造导向滤波器的理解，分别在二维和三维地震数据集上应用基于递归平滑滤波器和构造导向滤波器的 相干算法以及传统的 相干算法，并比较它们的计算效率和结果最后，应用递归导数滤波器计算了最常用的最正和最负体曲率，并与传统的曲率实现方式对比计算效率和结果 平滑滤波器本文主要研究应用递归滤波器和构造导向滤波器求取地震几何属性时的计算效率和结果首先，在本节简要介绍本文中涉及到的几种滤波器的实现原理，并评估它们的运算效率和平滑效果涉及到矩形滤波器、高斯滤波器、指数滤波器和构造导向滤波器，它们都可以实现低通滤波或进行平均计算，但它们的计算效率和平滑效果各不相同滤波器的递归实现能够显著提升计算效率（，）；递归高斯滤波器和递归指数滤波器的可分离性可以极大的提高多维数据滤波的计算速度此外，本节证明了构造导向 滤波 器适用 于平滑 富 含 方 向 信 息 的 地 震图像 平滑滤波器的一维冲激响应这一部分将简要介绍不同滤波器的定义，并分析它们各自的特点其中，矩形滤波器的一维冲激响应序列表达式如下：狔犻犕犕犼犕（狓犻犼），（）其中狓和狔分别表示滤波器的原始输入和最终的平滑结果，犕表示计算窗口的半宽度矩形滤波器以相等的非零权重犕计算分析窗口中的所有样本而对于计算窗口之外的样本，计算权重都为零公式（）表明，矩形滤波器的功能是对分析窗口内所有样本的值进行平均计算，并将平均值作为最终输出当犕 时，矩形滤波器的加权序列（或冲激响应）如图中用黑色表示但是在计算均值时，像这样从非零权重到零权重的突变是没有意义的，因此引入权重分配更加合理的高斯滤波器和指数滤波器高斯滤波器的一维冲激响应序列是将一维冲激函数与高斯函数卷积得到的如下数学表达式：狔犻犮犼狓犻犼犼，（）其中犕（犕）槡表示与矩形滤波器在平滑半径为犕时的滤波效应相当的高斯滤波窗口半径，犮为常数算子犼表示高斯滤波器的权重公式（）说明了高斯滤波器的计算权重取决于样本与中心的距离：样本点距离滤波窗口的中心越近，则权重越大，距离越远则权重越小，且权重始终大于零地震相干属性刻画的是与目标地震道相邻的几个地震道之间的相似性大小，与高斯滤波器的权重分配一致，因此高斯滤波器适合于计算相干等地震属性图中红色序列表示与半窗宽犕 的矩形滤波序列滤波效应相当的高斯滤波的冲激响应犕（犕）槡（）类似于高斯滤波器，指数滤波器的一维冲激响应序列是将一维冲激函数与指数函数卷积得到的，具有以下形式：狔犻犪犪犼狓犻犼犪犼，（）其中犪的范围是，），且由指定的半窗宽确定：犪槡，其中，犕（犕）槡算子犪犼表示指数滤波器的权重公式（）表明，除了输入和输出向量边缘附近的样本外，指数平滑滤波序列的权重分布是对称的，从其峰值到两侧呈指数衰减，且指数滤波器的脉冲响应始终大于零与高斯滤波序列相同，与半窗宽为犕 的矩形滤波器滤波效应相当的指数滤波器的一维冲激响应在图中用蓝色表示（）从以上三种常用滤波器的冲激响应（图）来看，矩形滤波器在平滑窗口内具有相同的权重系数，高斯和指数滤波器的权重系数均随着距离增加而衰减，且指数滤波器的衰减趋势快于高斯滤波器一个期张文等：地震几何属性的快速算法实现合理的滤波器应该在其平滑窗口中心位置处的样本分配相对多的权重，而远离中心位置处的样本分配相对少的权重（，）显然，相对于矩形滤波器，高斯和指数滤波器具有更合理的加权平均效应因此，我们将利用高斯和指数滤波器来实现相干属性计算中基于窗口的加权平均计算然而，高斯和指数滤波器中的空变加权系数的计算会增加滤波或加权平均的计算复杂度，为此我们考虑利用高效的递归实现算法来对它们进行提速图滤波器的一维冲激响应矩形滤波器（黑色）、高斯滤波器（红色）和指数滤波器（蓝色）的一维冲激响应（引自 （）；对于这些滤波器，每个输出样本是其附近输入样本的加权平均值；在此显示的是每个滤波器的标准化的权重分布，权重的和均为 （），（），（），（）；，；，滤波器的递归实现相对于矩形滤波器，高斯滤波器和指数滤波器的运算过程较为繁琐，因此引入递归算法来实现高效的计算对占用大量内存的输入数组，递归滤波器的优势尤为明显矩形滤 波 器 的 递 归 算 法 可 以 通 过 以 下 方 式实现：狔狀狔狀（狓狀犕狓狀犕）（犕）（）对于高斯滤波器的递归实现，可以通过结合因果系统和反因果系统来构造具有近似高斯函数对称脉冲响应的高阶递归系统 （）提到，无限脉冲响应（）近似高斯滤波的计算代价与计算半窗口宽度无关，类似于矩形滤波器的递归算法（），高斯滤波器可以递归地解决如下的差分方程：狔狀犫狓狀犫狓狀犫狓狀犪狔狀犪狔狀，（）其中，滤波器系数犫、犫、犫、犪和犪是高斯滤波半窗宽的函数然而，由于单一差分方程的递归解只有个参数，无法使系统的脉冲响应很好地拟合高斯函数此外，公式（）系统的脉冲响应是单侧的，而高斯函数是对称的因此，对称过滤需要因果系统和反因果系统的一个解序列公式（）中 滤波器的狕变换为：犎（狕）犫犫狕 犫狕 犪狕 犪狕（）假设犎（狕）的所有极点都位于复狕平面的单位圆内，那么该二阶系统是因果稳定的相应的稳定反因果二阶系统如下：犎（狕）犫犫狕犫狕犪狕犪狕，（）犎（狕）的所有极点都位于复狕平面的单位圆外这两个二阶递归系统都不足以拟合高斯滤波，但通过把像这样的因果系统和反因果系统结合起来，就可以构造具有近似高斯函数的对称脉冲响应的高阶递归系统 （）以及 等（）提出了以下高斯滤波器的递归实现方式，已广泛应用于计算机视觉领域（，）（）将递归高斯滤波器构造为因果系统和反因果系统的和得到以下复合系统：犎（狕）犎（狕）犎（狕）（）对于 的四阶滤波器，因果系统为：犎（狕）犫犫狕 犫狕 犫狕 犪狕 犪狕 犪狕 犪狕，（）反因果系统为：犎（狕）犫狕犫狕犫狕犫狕犪狕犪狕犪狕犪狕（）为了实现复合系统犎（狕），需要对输入序列狓狀 应用因果和反因果过滤器，并将结果按照时间顺序逐个累积到输出序列狔狀 中其中，因果过滤器和反因果过滤器可以并行实现公式（）和（）中的系数犪，犪，犪和犪只与滤波器极点的位置有关，并且对于因果系统犎（狕）和反因果系统犎（狕），这些系数是相同的公式（）里因果系统犎（狕）中的系数犫、犫、犫和犫依赖于滤波器零值的位置因为复合滤波器犎（狕）必须是对称的，即犎（狕）犎（狕）犎（狕）犎（狕）犎（狕）犎（狕），那么公地 球 物 理 学 报（）卷式（）中反因果系统犎（狕）的系数犫、犫、犫和犫可以通过系数犪，犪，犪、犪、犫、犫、犫和犫的简单的加减乘法运算得到也就是说，因果系统犎（狕）的个系数决定了复合系统犎（狕）的对称脉冲响应，而这个系数取决于因果系统犎（狕）的极点和零值位置 通过计算脉冲响应犺狀 与半窗宽为的采样高斯函数犵（狀；）的平方差总和犈的最小值来获取滤波器极点和零值位置其中，犈的表达式如下：犈犖狀犺狀犵（狀；）（）等（）将递归高斯滤波器构造为因果系统和反因果系统的乘积得到以下复合系统：犎（狕）犎（狕）犎（狕）（）对于 等人的四阶滤波器，因果系统为：犎（狕）犫犪狕 犪狕 犪狕 犪狕，（）反因果系统为：犎（狕）犫犪狕犪狕犪狕犪狕（）为了实现复合系统犎（狕），首先需要对输入序列狓狀 应用因果滤波器，得到中间输出序列狔狀，然后对序列狔狀 应用反因果滤波器，得到最终的输出序列狔狀这个实现过程是将因果滤波器和反因果滤波器进行串联应用同样的，公式（）和（）中的系数犪，犪，犪和犪只取决于滤波器极点的位置，且对于因果系统犎（狕）和反因果系统犎（狕）而言，这些系数是相同的复合滤波器犎（狕）的对称性可以通过犎（狕）犎（狕）犎（狕）犎（狕）得到验证 等人通过最小化滤波器频率响应的均方根误差或最大误差来计算滤波器的极点对于高斯滤波器而言，通过最小化均方根和最大误差计算得到的极点没有显著差异（，）指数递归滤波器的实现比高斯递归滤波器更加高效 （）提出，可以通过 行语言代码实现与指数平滑滤波效果近似的、标准化权重（权重之和为）的递归指数平滑滤波（附录）其中，程序的输入是一个浮点参数犪（即公式（）中的犪）和一个由狀个浮点值组成的序列狓犻（犻，狀），输出是一个由狀个长度相同的浮点数狔犻 组成的序列该算法优势除了体现在程序简洁之外，程序的实现还可以使输出值狔犻 与输入值狓犻 存储在同一个指针，这样在对大型数组进行指数平滑滤波时，可以极大的节约内存此外，对于任意大小的计算窗口，递归指数滤波器的浮点型操作为，因此，相对于递归矩形滤波器和递归高斯滤波器，递归指数滤波器的浮点型操作是最少的，且与计算窗口大小或平滑效果无关通过以下实例计算直观的描述递归高斯和递归指数滤波器的计算效率：对于长度为狀的一维数组输入，计算半窗口宽度为犕，设狀 ，犕，根据公式（），如果使用最常用的矩形滤波器，将要进行狀（犕）次浮点型操作，即 次操作；如果使用递归矩形滤波器，将要进行（犕）（狀）次浮点型操作，即 次操作；如果使用递归高斯滤波器，不论计算窗口大小是多少，每个点的平滑均需要 次浮点型操作，即对于输入数组总共将进行 次浮点型操作；如果使用递归指数滤波器，每个点的平滑均需要次浮点型操作，与计算窗口大小无关，即对于该输入数组总共将进行 次浮点型操作并且，在增大计算半窗口宽度犕时，矩形滤波器的浮点型操作次数将随之明显增多，而计算窗口的增大并不影响递归滤波器的浮点型操作次数或计算效率然而，即使递归矩形滤波器所需要的浮点型操作相对递归高斯和递归指数滤波器是最少的，根据 节中对于三种滤波器滤波效果的图件描述结合地震属性的计算需求，通常在计算地震属性的时候优先考虑使用递归高斯和递归指数滤波器高斯滤波器和指数滤波器的递归实现可以极大的节约计算成本，提高计算效率对于二维或三维输入数据的平滑操作，在递归的基础上，还可以利用滤波器的可分离性进行并行操作，从而实现更高效的计算 多维高斯和指数滤波器的可分离性为了验证递归高斯滤波器和递归指数滤波器的可分离性，以及探索二维情况下不同平滑滤波器的平滑效果，使用维度为 的二维零矩阵作为输入，并将中心点的值赋为 ，如图 所示滤波器的可分离性指的是：先后沿着不同坐标轴方向进行平滑滤波等效于对计算窗口内所有点进行整体的平滑滤波处理以二维高斯平滑滤波为例：犵狓 狔犮犼犻犘狓犻 狔犼犻狓犼（）狔犮犼（犻犘狓犻 狔犼犻狓）犼狔，（）期张文等：地震几何属性的快速算法实现图滤波器可分离性及平滑效果对比（）为输入图像；（）（）和（）分别为经过矩形滤波器、高斯滤波器、指数滤波器平滑处理后的输出图像；（）为对（）沿垂直方向进行高斯平滑滤波的结果；（）为在（）的基础上沿水平方向进行高斯平滑滤波的结果；（）为对（）沿垂直方向进行指数平滑滤波的结果；（）为在（）的基础上沿水平方向进行指数平滑滤波的结果 （）；（）（）（），；（）（），（）（）；（）（），（）（）地 球 物 理 学 报（）卷（）式说明，对于计算窗口内进行的整体高斯平滑处理（图）等价于先沿竖直方向进行高斯平滑（图），然后对于该平滑结果再进行沿水平方向的高斯平滑滤波（图）同样的，多维情况下指数平滑滤波器也具有像这样的可分离性（图，和）分别使用矩形滤波器、递归高斯滤波器和递归指数滤波器对输入图像（图）进行平滑处理，结果如图、和 所示，其中，矩形滤波器的半窗宽设置为犕，递归高斯滤波器和递归指数滤波器的半窗宽是将犕 代入前文所述滤波器半窗宽公式犕（犕）槡中计算得到的为了更直观地对比矩形滤波器和与之对应的其他滤波器的实际平滑范围，图 中不为零值的红色区域边界在图 和 中以红色方框表示图 和 中均有非零值部分在红色方框外，这表明矩形滤波器的实际平滑范围小于与之滤波效应相当的递归高斯滤波器和递归指数滤波器的平滑范围此外，递归指数滤波器从中心向四周衰减的速度明显快于递归高斯滤波器这些特点均与各滤波器的一维冲激响应特点一致利用递归高斯滤波器和递归指数滤波器的可分离性可以加速二维或三维数据的计算进程在二维情况下，设计算窗口大小为犖犖，如果对计算窗口进行整体的高斯或指数平滑滤波处理，那么需要对窗口内所有样本点进行与对应权重的相乘计算，以及带权重的样本点之间的求和计算，那么对于每一个计算窗口，浮点型操作次数为犖；同样的，如果对计算窗口进行分离的高斯或指数平滑滤波处理，那么对于计算窗口中的每一个方向进行平滑时所需要的浮点型操作次数为犖，即二维窗口中共需要（犖）次浮点型操作类似的，对于三维数据，设计算窗口大小为犖犖犖，对整个计算窗口进行高斯或指数平滑滤波处理需要进行的浮点型操作次数为犖，而对计算窗口进行分离的高斯或指数平滑滤波处理共需要（犖）次浮点型操作显然，当犖时，对于每个计算窗口，先后沿不同方向进行平滑处理所需要的浮点型操作远少于对整个计算窗口进行平滑滤波处理，并且犖越大，两种方式的浮点型操作次数相差越大（表）基于滤波器的可分离实现，可以将高维的高斯和指数滤波器分解为对输入高维数组依次沿各个维度方向的一维滤波由于一维滤波操作只与其滤波方向数值相关，而与其他维度方向数值无关，因而更表不同半窗宽下对于每个样本点，使用整体的递归高斯指数滤波器和分离的递归高斯递归指数滤波器所需的浮点型操作次数对比犜 犪 犫 犾 犲 犆 狅 犿 狆 犪 狉 犻 狊 狅 狀狅 犳 犳 犾 狅 犪 狋 犻 狀 犵 狆 狅 犻 狀 狋 狅 狆 犲 狉 犪 狋 犻 狅 狀犮 狅 狌 狀 狋 狊狆 犲 狉 狊 犪 犿 狆 犾 犲狆 狅 犻 狀 狋（犉 犔 犗 犘 犛）犳 狅 狉狑 犺 狅 犾 犲 狉 犲 犮 狌 狉 狊 犻 狏 犲犌 犪 狌 狊 狊 犻 犪 狀犲 狓 狆 狅 狀 犲 狀 狋 犻 犪 犾 犳 犻 犾 狋 犲 狉 狊犪 狀 犱 狊 犲 狆 犪 狉 犪 狋 犲 犱狉 犲 犮 狌 狉 狊 犻 狏 犲犌 犪 狌 狊 狊 犻 犪 狀犲 狓 狆 狅 狀 犲 狀 狋 犻 犪 犾 犳 犻 犾 狋 犲 狉 狊狌 狀 犱 犲 狉犱 犻 犳 犳 犲 狉 犲 狀 狋犺 犪 犾 犳 狑 犻 狀 犱 狅 狑狑 犻 犱 狋 犺 狊半窗宽犕（采样点）整体的递归高斯指数滤波器分离的递归高斯指数滤波器二维情况 三维情况 加方便多线程并行实现，进而可以实现对滤波过程的进一步加速 构造导向的平滑滤波器 （）提出，由张量场参数化的构造导向的平滑滤波器通常比一般的平滑滤波器更通用构造导向的平滑滤波器可以用来实现沿倾角的 相干算法以及构造导向的 相干算法这一部分将以二维合成地震图像（图）为例，介绍构造导向的平滑滤波器的实现效果在构造导向的滤波中，方向和维度是两个关键因素，可以由结构张量提供，即：犜犵 犵，（）其中犵表示输入图像的梯度，表示高斯平滑算子在二维情况下，特征分解可以使用归一化特征向量狌和狏以及相应的特征值狌和狏来表示：犜狌狌 狌狏狏 狏，（）其中，特征向量狌的方向垂直于原始图像的切线方向（图 黄色短线），而特征向量狏的方向与原始图像的切线方向一致（图 黄色短线）结合结构张量场犜，可以通过求解以下偏微分方程实现高效的构造导向的平滑滤波（，）：犵（狓）犇（狓）犵（狓）犳（狓），（）其中犳（狓）表示输入图像，犵（狓）表示平滑后的图像，为常数参数，犇（狓）为平滑张量场，犇（狓）与结构张量犜共享特征向量狌和狏，通过赋予特征值狌期张文等：地震几何属性的快速算法实现图构造导向平滑滤波器在二维合成地震数据中的应用实例（）为二维合成地震数据；（）和（）分别为垂直于构造方向的特征向量狌以及平行于构造方向的特征向量狏；（）和（）分别是对于测试数据沿特征向量狌以及特征向量狏进行各向异性平滑的结果；（）为对地震数据沿构造方向的平滑结果 （）；（）（）狌 狏 ，；（）（）狌 狏，；（）和狏不同值来实现沿不同方向的平滑以一个人为设定的行列 个元素值为、其余元素值为的二维数组为例，取图 中的结构张量场，设置狌，狏 ，对该二维数组进行构造导向的平滑滤波，可以得到效果与图 黄色短线方向一致的图；当设置狌 ，狏时，得到结果与图 黄色短线方向一致的图 像这样通过给特征值狌和狏赋予不同值，可以实现沿不同方向的平滑滤波同样的，对于二维合成地震数据（图），设置狌 ，狏，可以对该地震数据实现沿构造方向的平滑（图）此外，参数可以被犇（狓）吸收，但是为了方便控制平滑的程度，一般将它们分开表示具体而言，当 时，根据公式（）可以得到犵（狓）犳（狓），这就意味着没有进行平滑处理与半窗宽为犕的矩形滤波器滤波效应相当的构造导向平滑滤波器的半窗宽为犕（犕）滤波器的计算效率及平滑效果对比根据 和 节所述，矩形滤波器的实际计算窗口比与之相当的其他平滑滤波器的实际计算窗口小因此，为了更好的比较几种滤波器的计算效率和平滑效果，在这一部分的实验中，递归高斯滤波器、递归指数滤波器和构造导向的平滑滤波器的半窗宽使用稍小的窗口犿犕 计算，其中犕 是本次实验中矩形滤波器的半窗宽大小梵高著名作品 星夜（图）中平滑弯曲的线条包含丰富的方向信息，适合用来评估平滑滤波器尤其是构造导向的平滑滤波器的性能在维度为 的原始图像（图）中，均匀的抽取 个子样本并将它们分配到维度为 的二维零矩阵的对应位置处作为输入图像（图）使用不同的滤波器对输入图像进行平滑处理，得到的结果分别如图、和 所示对于矩形滤波器的处理结果（图），除了以 个子样本点为中心的矩形范围内的值等于对应分析窗口内所有样本的平均值以外，其余区域的值都为零经过递归高斯滤波器处理后得到的图像（图）包含 个以子样本点为中心的白色非零值区域，说明子样本点及其邻近的样点分配到较大的权重经过递归指数滤波器处理后的图像（图）中白色区域的分布相对于经过递归高斯滤波器处理后的图像（图）而言更集中，这意味地 球 物 理 学 报（）卷图不同滤波器对于测试模型的平滑处理结果（）为原始图像；（）（）分别为从原始图像中抽取部分样点得到的输入图像、；（）（）（）和（）分别为使用矩形滤波器、递归高斯滤波器、递归指数滤波器和构造导向的平滑滤波器对输入图像进行平滑处理所得结果；（）为使用构造导向的平滑滤波器对输入图像进行平滑处理所得结果 （）；（）（）；（）（）（）（），；（）着指数滤波器会给子样本点周围较小范围内的元素分配相对大的权重根据 节所述，在地震数据处理过程中通常需要考虑构造信息，可以通过使用构造导向平滑滤波器，调节特征值的大小来实现对地震数据的各向异性平滑滤波因此，结合原始图像（图）的结构张量场，对（图）进行构造导向平滑滤波处理，设狌、狏 （或其他极小的值），可以实现沿原始图像中线条切线方向的对于图 的各向异性平滑，如图 所示为了更好地展示构造导向平滑滤波器的特性，以每 个像素点为间隔，对原始图像（图）进行采样，然后将抽取到的子样本点的值分别赋值到维度为 的二维零矩阵的对应位置处，作为新的输入图像（图）对进行与上述相同的构造导向的各向异性平滑处理，输期张文等：地震几何属性的快速算法实现出图像（图）中的线条方向与原始图像中曲线的局部方向一致为了测试各平滑滤波器的计算效率，通过不断增加矩形滤波器半窗宽犕的大小来比较相应的各个滤波器的计算耗时情况，递归高斯滤波器、递归指数滤波器和构造导向的平滑滤波器的半窗宽随着犕的改变而改变，其中，根据 节中对于各滤波器的介绍，与半窗宽为犕的矩形滤波器滤波效应相当的递归高斯滤波器和递归指数滤波器的半窗宽相同，皆为犿（犿）槡，与半窗宽为犕的矩形滤波器滤波效应相当的构造导向平滑滤波器的半窗宽犿（犿）如表所示，当矩形滤波器半窗宽犕不断增大时，在本次对比实验的半窗宽的变化范围内，矩形滤波器的耗时明显增加，递归高斯滤波器和构造导向的平滑滤波器的耗时增加较小，递归指数滤波器的耗时基本不变矩形滤波器的耗时始终大于计算窗口与之相当的递归平滑滤波器和构造导向滤波器，并且半窗宽越大，耗时差距越明显由此可见，递归滤波器和构造导向滤波器在计算效率方面表现出很大优势，可以将其应用到 相干算法以及体曲率的计算来提高计算效率表不同计算窗口下平滑滤波器的耗时（单位：狊）犜 犪 犫 犾 犲犆 狅 犿 狆 狌 狋 犪 狋 犻 狅 狀犮 狅 狊 狋 狊狅 犳 狊 犿 狅 狅 狋 犺 犻 狀 犵犳 犻 犾 狋 犲 狉 狊 犳 狅 狉犱 犻 犳 犳 犲 狉 犲 狀 狋 狊 犻 狕 犲 狊狅 犳犪 狀 犪 犾 狔 狊 犻 狊狑 犻 狀 犱 狅 狑 狊（狌 狀 犻 狋：狊）半窗宽犕（采样点）矩形滤波递归高斯滤波递归指数滤波构造导向的平滑滤波 第二代相干属性 的快速算法实现相干属性能够增强地震数据中由沉积过程、储层边界或构造运动等造成的不连续特征（，）传统方法中第二代相干属性 被定义为分析窗口中平均地震道的能量与所有地震道的平均能量的比值，测量的是相邻地震道之间的波形形状和振幅的相似性（，；，）（）提出，传统定义方法中的“平均”计算可以被加权，并将第二代相干属性重新定义为输入的地震数据图像的平滑序列的平方与输入数据平方的平滑序列的比值 二维场景将前面提到的滤波器应用于图 所示的二维合成地震剖面，使用传统计算方法以及 （）提出的方法，在不同的分析窗口下进行 的计算二维情况下 的传统算法可以表示为（，；，）：犛犻，犼犼犖狕犼犖犻犕狓犻犕犳狓，狕犕熿燀燄燅犖犼犖狕犼犖犻犕狓犻犕（犳狓，狕）犕烄烆烌烎犖犼犖狕犼犖犻犕狓犻犕犳狓，狕（犕）犼犖狕犼犖犻犕狓犻犕（犳狓，狕），（）其中，犳狓，狕 表示振幅，犛犻，犼 表示输出的计算窗口中心点（图 中黄色点）的 结果，该公式对应的计算窗口形状如图 中红色网格所示使用更简练的矩阵表示法，公式（）可以写成外积 的形式：犛狌犆 狌（犆），（）其中，犆表示协方差矩阵，且有犆犼犖狕犼犖犳 犳，（犆）表示协方差矩阵犆的主对角线上的元素之和（从左上角到右下角），狌是归一化向量（狌）：狌（犕）熿燀燄燅（）为了更简单的表达以及实现更高效的计算，公式（）中求均值的部分可以用不同权重分布的平滑算子 来代替（，），例如：犻犕狓犻犕犳狓，狕犕犳狓，（）那么，地 球 物 理 学 报（）卷图使用不同算法求取二维合成地震剖面的 相干属性（）中红色窗格表示不考虑倾角的分析窗口；（）中红色窗格表示倾角导向的分析窗口；（）中红色窗格表示构造导向的分析窗口；（）为使用一般方法计算得到的 属性；（）和（）分别为使用递归高斯滤波器和递归指数滤波器得到的 属性；（）和（）分别为使用递归高斯滤波器和递归指数滤波器计算得到的倾角导向的 属性；（）为使用构造导向的平滑滤波器计算得到的构造导向的 属性（）（）和（）中黄色虚线椭圆框标识出陡倾构造造成的伪影；（）和（）中蓝色箭头指示断层的片段式显示；（）中洋红色箭头指示断层假象 （）；（）；（）；（）；（）（），；（）（），；（）（）（）（），；（）（）；（）犼犖狕犼犖犻犕狓犻犕犳狓，狕犕熿燀燄燅犖 犳狓狕（）由此，公式（）可以改写成：犛狓，狕 犳狓狕 犳狓狕，（）其中，狓和 狕分别表示沿着狓方向和狕方向的平滑算子，这些算子可以用第节介绍的矩形滤波器、递归高斯滤波器以及递归指数滤波器来实现，对应的结果分别如图、和 所示但是图、和 中，在构造陡倾的位置通常伴随有明显的伪影（黄色虚线圈内）出现，这就会影响对断层的刻画效果为了消除这些伪影，需要使用构造导向的平滑期张文等：地震几何属性的快速算法实现滤波器调整计算窗口的形状从而实现倾角矫正将公式（）分子和分母内层的平滑算子由沿狓轴方向改为沿倾角方向（图 黄色短线）并用构造导向的平滑滤波器来实现，这样窗口的形状变为如图 所示的空变红色网格，即窗口横向上弯曲并与地震构造方向保持平行，公式（）也相应的变为：犛 犳狏狕 犳狏狕，（）由此得到倾角导向的 属性犛 当沿狕轴方向计算的算子 狕分别由递归高斯滤波器和递归指数滤波器实现时，相干结果分别如图 和 所示，图、和 中黄色虚线圈里标出的由倾斜构造造成的伪影在图 和 中基本消除，但是倾角导向的 属性中断层的显示是片段式的（图 和 中蓝色箭头），这可能是由于分子和分母中的外部平滑算子仅负责垂直方向上的平滑因此，倾角导向的 仍不足以处理具有水平和垂直方向特征的图像为此，将公式（）分子和分母外层的平滑算子由沿狕轴方向改为垂直于倾角方向（图 黄色短线），并且也用构造导向的平滑滤波器来实现，这样窗口的形状变为如图 所示的空变红色网格，即窗口横向上与构造平行，竖直方向上与构造保持垂直，公式（）也相应的变为：犛 犳狏狌 犳狏狌，（）由此得到构造导向的 属性犛 ，相干结果如图 所示相对于图 和 而言，图 中提取的断层变得更加连续，但是只有当断层附近的构造倾角方向垂直于断层面时，断层的连续性才会被加强，否则就会产生图 中洋红色箭头所示意的假象因为在地震解释工作中，研究对象往往是三维地震数据，并且在三维窗口中分析目标地震道及其相邻地震道的相似性将会更加全面和准确，故基于以上二维情况下的理论和结果的分析，需要将其引申到三维场景，进行进一步的对比实验从而获取更有实际意义的信息 三维场景在这一部分介绍三维情况下的 相干算法的计算方法，并将传统的 相干算法以及基于递归滤波器的快速算法应用到三维实际地震数据，比较运算效率和实验结果与二维场景类似，三维情况下 的传统算法如下：犛犻，犼，犽犽犔狕犽犔犼犖狔犼犖犻犕狓犻犕犳狓，狕（犕）（犖）犽犔狕犽犔犼犖狔犼犖犻犕狓犻犕（犳狓，狕），（）同样的，公式（）也可以用协方差矩阵表示，得到外积的，并通过 软件实现为了提高计算效率，公式（）中计算均值的部分可以进行加权计算并用平滑算子 代替，从而改写成便于用递归滤波器实现的简洁形式：犛狓，狔，狕 犳狓，狔狕 犳狓，狔狕，（）其中，平滑算子 可以使用矩形滤波器、递归高斯滤波器和递归指数滤波器实现根据地震数据及其倾角信息，可以计算倾角导向的，从而消除倾斜构造带来的影响：犛 犳狏，狑狕 犳狏，狑狕，（）其中，狏，狑通过在使用构造导向的平滑滤波器时，设置狏，狑 和狌 实现沿构造方向的平滑，而 狕可以选择使用递归高斯滤波器或递归指数滤波器实现为了使得显示的断层更加连续，通过将 狕换为 狌，得到构造导向的：犛 犳狏，狑狌 犳狏，狑狌（）三维定量基准实验用到的输入数据是截取的荷兰北海近海采集的高质量三维归档地震数据中富含断层构造的部分（图），数据维度为 将前文提到的三维 传统算法以及基于滤波器的快速算法应用到图 所示的数据体，以比较这几个 方 法 的 运 算效 率 和 实现 效果本 文 实验 中，和 方向半窗宽犕采用个样本点，时间轴方向半窗宽犕采用个样本点应用不同方法计算 的耗时情况如表所示，在不考虑倾角矫正的情况下，使用递归滤波器可以极大的加速公式（）的运算过程相对而言，递归指数滤波器要比递归高斯滤波器更加高效，但是高斯滤波器的权重分布相对更合理图 所示的地震图像犳（狓，狓，狓）在 和 方向的倾角信息分别为：狆（狓，狓，狓）（犳 狓）（犳 狓）以及狆（狓，狓，狓）地 球 物 理 学 报（）卷图三维地震数据体和使用不同算法求取的 相干属性（）（）和（）分别为三维地震数据体及其 方向的倾角和 方向的倾角；（）（）和（）分别为使用不考虑倾角的矩形滤波器、递归高斯滤波器和递归指数滤波器计算的 属性；（）为使用递归高斯滤波器和构造导向的平滑滤波器求得的倾角导向的 属性；（）是使用构造导向的平滑滤波器求得的构造导向的 属性；（）是使用传统方法求得的倾角导向的 属性（）、（）和（）中洋红色和蓝色箭头分别指示两处伪影，而（）、（）和（）中洋红色和蓝色箭头分别指示的这两处位置无伪影 （）（）（），；（）（）（），；（）；（）；（）（），（）（），（），（）（）（犳 狓）（犳 狓），相关程序可参考 ：同时输入地震数据（图）以及地震数据的倾角信息（图 和），根据公式（）可以使用构造导向的平滑滤波器和递归滤波器来实现倾角导向的 算法，这比使用协方差矩阵计算并通过插值进行倾角补偿的传统算法的计算效率高 倍以上根据公式（）使用构造导向的平滑滤波器计算构造导向的 属性，其计算效率比传统的倾角导向的 算法高倍左右对于由不同算法求得的相干属性，如图所示，在不考虑倾角矫正的情况下，使用矩形滤波器、递归高斯滤波器和递归指数滤波器得到的相干结果相近，在倾斜构造附近伴随有伪影，影响断层的识别；而使用递归滤波器和构造导向的平滑滤波器计算得到的倾角导向的 属性以及使用构造导向的平滑滤波器得到的构造导向的 属性都与传统方法计算得到的 属性结果相似因此可以借助递归滤波器和构造导向的平滑滤波器简单快速地实现 相干算法，从而提高地震构造解释的效率期张文等：地震几何属性的快速算法实现表应用不同方法计算犆 相干属性时的耗时比较犜 犪 犫 犾 犲 犆 狅 犿 狆 犪 狉 犻 狊 狅 狀狅 犳 犮 狅 犿 狆 狌 狋 犪 狋 犻 狅 狀犮 狅 狊 狋 狊狑 犺 犻 犾 犲犪 狆 狆 犾 狔 犻 狀 犵犱 犻 犲 狉 犲 狀 狋犪 狆 狆 狉 狅 犪 犮 犺 犲 狊 犻 狀狋 犺 犲犮 犪 犾 犮 狌 犾 犪 狋 犻 狅 狀狅 犳犆 犮 狅 犺 犲 狉 犲 狀 犮 犲犪 狋 狋 狉 犻 犫 狌 狋 犲 类型实现途径耗时（）犛狓，狔，狕基于矩形窗的常规实现 犛狓，狔，狕多线程并行实现的递归高斯滤波器 犛狓，狔，狕多线程并行实现的递归指数滤波器 犛 多线程并行实现的递归高斯滤波器、构造导向的平滑滤波器 犛 多线程并行实现的递归指数滤波器、构造导向的平滑滤波器 犛 协方差矩阵、插值（实现）（，；，）犛 多线程并行实现的构造导向的平滑滤波器 工作站配置：（）（），核数；三维实际数据维度：曲率属性的快速实现体曲率能够增强地震数据中的裂缝、断层等线性构造特征，可以与相干属性互补，帮助地震解释工作者更加全面和准确地认识工区中的地质构造情况传统方法中通常使用倾角的一阶导数的体计算来表示曲率，结合“分数导数”调节波长大小，最终实现不同波长的曲率的计算，为解释工作者提供多样的视角，从而获取局部或者区域性的多尺度地质构造信息但是在实际工作中通常面临大批量的工区地震数据，使用传统的曲率算法将会带来较大的计算成本本文提出了使用递归导数滤波器来取代“分数导数”实现多波长体曲率的快速计算方法，旨在提高计算效率、节约实际生产工作中的计算成本 （）提出，对于一个提取出来的层位曲面，可以使用最小二乘拟合或其他近似方法确定如下形式的二次曲面狕（狓，狔）：狕（狓，狔）犪 狓犮 狓 狔犫 狔犱 狓犲 狔犳，（）通过对公式（）求一阶导数和二阶导数，参数犪、犫、犮、犱、犲和犳可以用同相轴倾角狆 狕 狓和狇 狕 狔的一阶导数表示：犪 狆 狓，（）犫 狇 狔，（）犮 狇 狓