基于构造导向滤波的高斯束层析速度建模方法及其应用

郑浩，蔡杰雄，王静波

(1.中国石化石油物探技术研究院,江苏 南京 211103; 2.中国石化勘探分公司,四川 成都 610041)

0 引言

近年来随着“两宽一高”地震采集和层析反演方法取得突破进展，地震成像分辨率得到有效拓宽[1-2], 高分辨率层析反演速度建模方法取得了长足的进步。传统的深度域射线层析反演速度建模方法的分辨率存在中波数段信息缺失问题[3]，其可以较好地恢复地下速度的背景低波数分量，但难以准确反演中波数段的高精度速度模型[4-5]。学者们提出波动方程层析反演方法[6]，将速度反演分辨率向中高波数端扩展。该类方法基于波动理论，避免了高频近似假设引入的误差，反演结果精度高，但该类方法存在计算效率问题，且反演精度非常依赖于数据品质，当数据较差时存在反演稳定性问题[7-8]，在实际应用中推广较少。为了降低射线类层析矩阵的稀疏性，同时克服波动方程层析的计算量与稳定性问题，诸多学者提出介于常规射线理论和波动方程理论之间的方法。早在1995年，Vasco[9]便提出一种介于射线和波动方程之间的胖射线层析方法，该方法通过将射线加宽，提高射线路径的投影范围，从而提高层析精度，但该方法缺乏严格的理论推导，仅是对波动理论感性认识，应用较为局限。菲涅尔体层析[10]及高斯束层析[11-12]的出现实现了对波动方程线性化近似的“束”层析。菲涅尔体层析利用射线周围的波动性质来拟合走时和波场振幅，相比于射线层析的高频近似假设，该方法引入了菲涅尔体的有限频效应，更接近真实的波场传播，这降低常规射线层析方法中矩阵的稀疏性，提高计算稳定性[13]。高斯束层析[14]实际上是将临近区域空间内的多条射线合并为一条高斯束射线从而建立层析方程，对于观测系统中的一个炮检对，通过提高射线覆盖范围来提高反演精度及分辨率，降低层析方程的病态性，且该方法相对于波动方程层析，计算速度快、计算结果稳定。但早期的高斯束层析理论推导不够严格，蔡杰雄[15]从高斯束偏移角度道集出发，在波动方程一阶Born近似和Rytov近似下，推导出了基于高斯束算子的成像域走时层析方程及其敏感度核函数，完善了高斯束层析反演的理论基础，推动了高斯束层析反演的实用化。

层析反演速度建模的另一个重要研究内容是各种约束条件的利用，即在反演过程中加入包括层位、构造产状、井速度等多种已知信息约束以实现正则化，提高反演稳定性及分辨率。目前，层析反演中的构造约束方式主要有两类，一类是在速度层析过程中引入地质层位信息，通过层约束提高反演精度，但这类方法需要解释人员进行地质构造的精细化解释，对解释结果有较高的精度要求，若解释结果出现错误，会使得层析结果误差较大。该类方法大都基于主观认知强行加入反演过程中，比如层位约束层析反演就是在本轮的反演中不修改层位而仅更新层位内速度值，待完成本轮反演后重新偏移成像，在新剖面上重新进行层位拾取，之后再加入到新一轮的层析速度更新中，通常我们称这种约束为“强”约束[16-17]。强约束的反演具有双刃性，它一方面可能快速促进收敛，但另一方面可能引入解释人员的主观错误信息，导致反演结果出现误差；另一类方法通过引入含有构造信息的正则化算子，用于导引速度更新方向，避免了解释人员的主观错误信息，但构造约束精度取决于正则化算子的精度。我们称这种约束为“软”约束。目前有多种预条件正则化算子，包括Sobel算子、Reborts算子、基于斜率的平滑滤波算子、基于扩散方程解析解的构造导向滤波算子等。这些方法通过提取地震偏移图像的构造特征，从而构建预条件算子，实现构造约束层析反演。

本文从高斯束层析核函数出发，通过引入含有构造信息的预条件模型正则化算子，推导得到构造导向滤波约束下的高斯束层析速度建模方法。在预条件模型正则化算子的设计上，本文引入结构张量算法，用于提取偏移图像中的平坦、边缘及角点区域，通过扩散方程解析解实现构造导向平滑滤波，建立层析反演正则化算子，实现反演过程的自动“软”约束，提高高斯束速度层析反演的精度与稳定性。

1 成像域高斯束层析方法

频率域高斯束成像条件可以表示为：

IGBM(x,θ,φ,w)=S(x,pS,w;xS)R*(x,pR,w;xS)

(1)

式中:IGBM(x,θ,φ,w)表示共成像点道集,S(x,pS,w;xS)表示来自激发点的高斯束下行波场,R(x,pR,w;xS)表示检波点接收的高斯束上行波场,x=(x,y,z)表示共成像点位置的坐标，pS和pR分别代表激发点和接收点位置的慢度矢量，θ表示反射张角，φ表示成像点的方位角，w=2πf表示角频率，*表示共轭。

从角度域高斯束偏移成像条件(1)出发，在波动方程的一阶Born近似下，可得到剩余时差与速度修正量的表达式，再引入Rytov近似可进一步简化得到剩余时差与速度修正量的线性关系式[18]：

(2)

其中，Δt(x,θ,φ,ω)为高斯束偏移的剩余时差，KT表示高斯束层析核函数，其可以表示为两部分：

基于式(2)表示的层析核函数，就可以建立类似于射线层析的方程:

AΔm=Δt。

(4)

图1 深度域高斯束层析反演核函数Fig.1 Gaussian beam tomography kernel function in imaging domain

A是利用式(2)层析核函数计算得到的敏感度核函数矩阵，Δm是待反演的速度更新量，Δt是共成像点道集(common image gathers,CIG)的走时残差，在反射点局部，可以假定炮点和检波点的慢度pS和pR是常数，那么该点的慢度可以表示为：

Δpz=Average(pS+pR) 。

(5)

通过在CIG道集上拾取剩余深度差Δz(RMO)，则可利用式(5)得到走时残差Δt：

Δt=ΔpzΔz。

(6)

通常，式(4)中的矩阵A是极度稀疏的，直接求解稳定性较差，且多解性较强，往往需要引入正则化项，加入一些先验约束，使解在约束条件下收敛。这里采用预条件模型正则化的方式加入先验信息，待反演的模型扰动量可表示为：

Δm=Du。

(7)

其中算子D表示预条件算子，即引入的先验信息，这样可以将Δm的求解问题转化为算子u的求解问题，将式(7)代入到式(4)中，可得：

ADu=Δt。

(8)

式(8)即为构建的预条件层析方程，其阻尼最小二乘方程可表示为：

DTATADu+εu=DTATΔt,

(9)

求解上式得到u，代入到式(7)中即可得到速度更新量Δm。

通过上述推导，实现了波动方程线性化的深度域走时层析反演方法，再与叠前深度偏移相结合进行迭代：速度误差带来成像不聚焦误差，表现为成像道集走时差；再利用成像误差反过来修正速度，实现两者的迭代收敛，获得更精确的速度场和成像结果。

2 基于结构张量构造导向滤波正则化技术

从式(9)中可以看出，预条件模型正则化的高斯束层析反演关键是算子D的构建，传统方法常选择空间平滑算子，主要是通过减小层析系数矩阵的条件数，缓解反问题病态性，在矩阵求解的稳定性方面有一定作用，但该方法容易模糊构造边界，其反演结果对于边界的刻画较差。本文采用各向异性扩散方程解析解[20]构建预条件算子D，其隐式表达式为如下的偏微分方程：

g(x)-α·T(x)g(x)=f(x) ,

(10)

其中:f(x)是待滤波的地震数据，g(x)是滤波后的地震数据，α是表示平滑力度强弱的常数；T(x)是从地震剖面中计算得到的结构张量。结构张量实际上相当于沿构造梯度的一个平滑解，它可以用于构造的倾角扫描，且平滑过程沿构造进行，达到保边界的效果。对于一个三维数据体而言，其某点的结构张量T(x)可以通过一个3×3的对称正定矩阵表示：

(11)

其中:dx、dy、dz分别为数据沿x，y，z三个方向的梯度。通过矩阵的特征值分解，上式可以表示为:

T(x)=λuuuT+λvvvT+λwwwT,

(12)

其中:λu、λv和λw为分布在0～1的正实数，表示矩阵T(x)的特征值，u、v和w是对应的特征向量。通常定义λu≥λv≥λw≥0，这种情况下，特征向量u表示梯度最大的方向，即表示垂直于构造倾角的方向。特征向量v和w分别表示inline和crossline平行于构造的方向，λ越大意味沿对应特征方向的平滑尺度越大。

这里以二维数据为例，对特征向量进行说明。二维情况下，式(12)只有前两项，图2a、b中黄色的椭圆分别是地震切片与剖面上结构张量计算结果，实际上在无构造(各向同性)区域，该点的结构张量形态是圆形，如图3a所示，此时u和v的长度相等；当存在构造形态(各向异性)时，就会出现椭圆形态的结构张量，如图3b所示，其中v表示平行于构造方向的特征向量，u表示垂直于构造方向的特征向量，这样我们就可以计算出所有点的结构张量。

图2 地震切片(a)及剖面(b)计算得到的结构张量Fig.2 Structural tensor calculated on seismic slice (a) and section (b)

得到每点的结构张量后，通过合理调整λu、λv、λw，从而实现沿层位和断层的构造导向滤波处理，平滑效果如图4所示。

图3 各向同性(a)和各向异性(b)的结构张量算子Fig.3 Isotropic(a) and anisotropic(b) structural tensor operators

图4 构造导向滤波前(a)后(b)对比，可同时加强层位与断层特征Fig.4 The comparison of (a) before and (b) after structural-guided filtering,which can enhance the structure and fault characteristics

引入有限差分近似，预条件算子D可写作：

D=(I+TT)-1,

(13)

其中T即为式(12)所求的结构张量，单位矩阵I的作用是在λu、λv、λw均为零时保持图像不改变，将式(13)所得结果代入式(9)中即可实现基于构造导向滤波的高斯束层析反演。此外，相比于其他的构造约束算法，该方法无需指定平滑角度，完全自动提取，计算灵活；且方法的抗噪性较好，即便是在信噪比极低的情况下，构造导向滤波算法也可以得到相对精确的结果。

3 模型及实际资料测试

3.1 二维层状模型数据测试

设计二维层状地质模型验证本文方法的精度，模型中填充速度属性值。所设计模型为速度沿纵向递增的层状模型， 其中模型中部有一高速隆起构造(图5a)，纵向采样点为601，采样间隔为10 m，横向共有1201道，道间距为10 m。利用声波波动方程正演得到单炮记录(图6)，观测系统采用中心放炮，两边接收，每炮801道。初始模型采用常梯度模型，见图5，其中第一层速度给成准确值。图5c、5d展示了常规高斯束层析与构造约束预条件高斯束层析的反演结果，显然，常规的高斯束层析对于浅层的水平层状构造反演结果尚可，但对于中部隆起构造反演结果较差，分辨率明显不足，且中深部反演结果误差较大；相比较而言，基于构造导向滤波的高斯束层析算法所得结果精度较高，与真实模型基本吻合，更有利于精确成像。

为了更加直观地看出造构造导向滤波高斯束层析的优势，在层析结果竖直抽取某一道进行对比，如图7a所示。图7b表示了理论速度值、平滑速度值和反演的速度值的对比曲线。，显然，基于构造导向滤波的高斯束层析反演方法不仅可以反演得到中波数分量的结果，提供更加丰富的速度信息，且能够达到保边界的效果，构造边界刻画更为清晰，这说明基于构造导向滤波的高斯束层析对于构造的刻画明显优于常规高斯束层析方法，且反演精度更高，减少了由于人为干预带来的误差，结果更加合理。

a—理论模型;b—初始模型;c—常规高斯束层析反演结果;d—构造导向滤波高斯束层析反演结果a—velocity model;b—initial model;c—updated velocity model by conventional Gaussian beam tomography;d—updated velocity model by Gaussian beam tomography based on structural-guided filtering图5 理论地质模型(速度值)及反演结果Fig.5 Velocity model and the inverted results

图6 模型正演的单炮记录Fig.6 The modelling common shot gathers

a—抽取模型中心道;b—真实值、常规平滑正则化约束层析反演结果及构造约束正则化层析反演结果单道对比a—the location of the extracted trace in the model;b—comparison of theoretical result,conventional smoothing regularization tomography inversion result and structural-guided regularization tomography result图7 反演结果单道对比分析Fig.7 Single trace comparison

深度域高斯束层析反演技术可以缓解层析矩阵的病态问题，反演得到速度模型的中波束分量，提高地震勘探的成像精度。基于此，通过引入基于结构张量的构造导向滤波算法用于约束速度更新量的计算，能反演出高精度的地下速度模型，增强断层突变位置的反演效果，弥补地震成像的中波数段缺失，提高地震成像精度。

3.2 实际地震资料应用

实际资料选自中国西南部复杂山地页岩气某工区。该地区地表起伏大，悬崖峭壁密布，山高谷深，数据信噪比非常低，且地下结构复杂，速度横向变化剧烈，速度建模难度大，成像困难，具有典型的复杂地表和复杂构造的“双复杂”特性，其对地震资料成像精度要求远高于常规勘探。

首先对该区域的资料进行预处理、叠前时间偏移、深度域初始建模及偏移成像等基础工作，得到高品质的CIG道集，在此基础上，通过基于构造导向滤波的高斯束层析反演得到深度域速度模型，结合克希霍夫叠前深度偏移得到最终的成像剖面。图8是某实际数据深度域常规高斯束层析建模和构造导向滤波高斯束层析速度建模结果对比，从两张对比图上可以看出，构造导向滤波高斯束层析速度建模结果表现了速度模型更多细节，清晰地刻画了速度纵横向空间变化，提高了速度建模的精度。

图8 某实际数据深度域常规高斯束层析建模(a)及构造导向滤波高斯束层析建模(b)结果对比Fig.8 Comparison of conventional Gaussian beam tomography (a) and structure-guided filter Gaussian beam tomography (b) in depth domain

速度模型的更新和剩余速度分析采用速度更新量和剩余速度谱联合质量控制方法，图9是质控图,图10是更新前后的CIG道集。从更新前的道集上看，存在同相轴未校平现象，更新后这些同相轴变平，图上前头所指。由于存在各向异性，还有一些同相轴未校平。另外，该区域存在多次波的影响(如图中红框区域所示)，因此仍有部分道集下拉现象。

图9 速度更新量(a)与剩余谱(b)联合质控Fig.9 Joint quality control of velocity update (a) and residual spectrum (b)

图10 速度更新前(a)后(b)CIG道集Fig.10 The CIGs before (a) and after (b) velocity update

图11是对应于图8速度模型范围的偏移结果，显然，采用相同的成像方法，构造导向滤波下的高斯束层析反演对应的偏移结果中断点更加连续，复杂构造区信噪比更高，整体成像质量明显优于常规方法所得结果，这对后续的构造落实有较好的指导意义。通过该实际资料处理过程可以认识到，基于构造导向滤波的高斯束层析反演速度建模能够反演得到速度的中高波数分量，所得成像结果细节更加丰富，信噪比高，成像质量好。

a—常规高斯束层析反演对应的偏移结果;b—构造导向滤波高斯束层析反演对应的偏移结果a—the imaging results corresponding to conventional Gausssian beam tomography;b—the imaging results corresponding to structural-guided filter Gaussian beam tomography图11 某实际数据PSDM成像结果Fig.11 PSDM imaging results

4 结论

基于波动方程线性化近似的深度域高斯束层析速度建模方法具有较高的反演精度与实用性。在此基础上，本文引入结构张量算法，推导了构造导向滤波正则化技术。该技术通过计算地震数据的结构张量，并利用各向异性光滑算子将其应用于高斯束层析反演中，作为正则化约束项，实现了构造导向滤波约束下的高斯束层析速度建模方法。相比于常规的高斯束层析速度建模技术，该方法完全基于数据驱动，实现了对层析反演的“软约束”，既缓解了层析反演求解的多解性，同时提高反演分辨率与稳定性，得到更加符合地质认识的高精度速度模型，有效改善成像效果，提高了高斯束层析反演对复杂模型的反演分辨率，为高斯束实用化推广铺平了道路。

进一步可考虑将基于结构张量的构造导向滤波算法应用于各向异性层析建模，通过更高精度且自动化地计算构造的方位角及倾角推进各向异性高斯束层析方法的实用化。