随着油气勘探开发程度的不断深入, 地质目标的尺度越来越小、埋藏越来越深、所面临的地震地质条件越来越复杂, 对地震资料处理(尤其是速度建模)的要求也越来越高, 迫切需要理论更先进、精度更高的地震速度建模及成像技术提供支撑。目前, 深度域速度建模起着至关重要的作用, 直接影响地震成像精度与质量, 进而影响后续整个探区地质构造认识、储层预测及工程钻井施工。因此, 深度域速度建模研究是一项更加综合、系统、也更加迫切的工作。无论是在学术界还是工业界, 许多国内外专家学者在此方面开展了大量的研究工作。其中, 射线类层析方法由于其成熟的理论基础及较高的计算效率^[1-2], 已被广泛应用于工业界, 但单纯的射线类层析方法存在明显的缺陷, 其只能反演速度的中低波数分量。随着波动理论的逐步成熟, 基于波动方程的层析反演算法出现^[3], 该类算法避免了射线类层析方法的高频假设问题, 理论上反演精度更高, 但在实际应用中, 由于计算量问题及对资料品质要求较为苛刻, 至今未能推广应用^[4-5]。为了减少波动类算法的计算量, 同时提高速度建模精度, 一种介于波动类与射线类的射线束算法应运而生, 该类方法由早期缺乏理论依据的胖射线层析发展到目前理论成熟的波动方程线性化层析反演算法(菲涅尔体、高斯束等), 发展较为迅速, 部分技术已实现工业化应用^[6-9]。

实用化的深度域层析反演是一套完整的技术流程, 除了层析反演核心算法之外, 往往还需要层位自动拾取、正则化约束等相应的配套技术才能应用于实际资料的处理, 并大规模生产应用。其中, 正则化约束对于层析反演的结果影响较大, 尽管正则化方式较为单一, 但具体实现方式较多, 效果也不尽相同。目前较为普遍的做法是在反演过程中引入构造、层位等先验信息, 使得反演结果更具有地质意义, 从而达到提高反演精度的目的^[10-11]。这类方法俗称"硬约束", 效果较为明显, 但会引入解释人员的主观认知, 一旦出现认识上的错误, 反演结果可能更差。与之相对应的"软约束"正则化算法通过完全的数据驱动引入地质信息, 实现构造约束。但不管是"软"约束还是"硬"约束, 其约束效果均依赖于正则化算子的精度及约束信息的可靠性。如何将断层信息作为正则化约束项引入层析反演过程已成为新的研究热点^[12-15]。针对复杂断裂系统, 目前已经有不少学者提出针对性的建模成像技术^[16-17], 对于深度域建模而言, 将断裂属性作为先验信息应用到正则化约束中, 常规做法是将人工解释的断裂信息应用到速度模型中, 该策略可以在一定程度上提升断裂周边速度建模精度, 但断裂信息依赖于解释人员的主观认识, 容易引入误差, 断裂建模成像精度提升有限。

基于此, 本文采用数据驱动的方式实现基于断裂属性约束的深度域层析速度建模技术, 以高斯束层析反演技术为基础, 通过最优路径寻优算法计算数据体的断裂属性, 以"软约束"的方式将其作为先验信息实现预条件正则化, 搭建基于断裂属性体约束的深度域层析速度建模框架, 实现针对断裂的高精度速度建模, 提高断裂发育区域的速度建模精度。

1 深度域层析反问题的建立

深度域层析反演通常可以表示为如下矩阵形式:

$ \Delta \mathit{\boldsymbol{t}} = \mathit{\boldsymbol{K}}{}\Delta \mathit{\boldsymbol{s}}{} $

(1)

式中: Δt表示剩余时差; Δs表示待求的慢度更新量; K是敏感度核函数矩阵, 不同的求解算法, 其表现形式不同。采用高斯束算子计算层析核函数, 传播方向为p的核函数可以表示为^[18]:

$ \begin{array}{l} \boldsymbol{K}\left( {x, {x_{\rm{r}}}, \boldsymbol{p};{x_{\rm{s}}}} \right) = \int W (\omega ){\mathop{\rm Im}\nolimits} \left[ {\frac{{2\omega }}{{v_0^3(x)}} \cdot } \right.\\ \;\;\;\;\;\;\;\;\;\;\left. {\frac{{{G_0}\left( {{x_{\rm{r}}}, \boldsymbol{p}, \omega ;x} \right){G_0}\left( {x, \boldsymbol{p}, \omega ;{x_{\rm{s}}}} \right)}}{{{G_0}\left( {{x_{\rm{r}}}, \boldsymbol{p}, \omega ;{x_{\rm{s}}}} \right)}}} \right]{\rm{d}}\omega \end{array} $

(2)

其中, K(x, x_r, p; x_s)表示波场中某点x到反射点x_r的层析核函数, G₀(x_r, p, ω; x)表示波场中某点x到反射点x_r的高斯束传播算子, G₀(x, p, ω; x_s)表示波场中炮点x_s到某点x的高斯束传播算子, G₀(x_r, p, ω; x_s)表示波场中炮点x_s到反射点x_r的高斯束传播算子, Im表示取虚部, v₀(x)代表x点的初始速度, W(ω)表示归一化加权函数, ω表示归一化加权函数的自变量。本文采用高斯函数$g\left( \omega \right) = \exp \left\{ { - \left[ {{{\left( {\omega - {\omega _0}} \right)}^2} - \left( {2{\delta ^2}} \right)} \right]} \right\}/\left( {\sqrt {2{\rm{ \mathit{ π} }}} \delta } \right)$来求解W(ω), δ表示标准差, 以ω₀为高斯函数的期望值(对称中心), Δω为高斯函数的标准差(半宽度), 那么$W\left( \omega \right) = g\left( \omega \right)/\mathop \smallint \limits_{{\omega _0} - \Delta \omega }^{{\omega _0} + \Delta \omega } g\left( \omega \right){\rm{d}}\omega $。图 1展示了反射波高斯束算子传播路径。相比于射线, 高斯束展布范围更宽, 反演范围更大。实际上, 高斯束层析作为一种波动方程线性化方法, 本质上属于波动类算法, 该方法更接近波的实际传播情况, 精度更高。

通常, 公式(1)中的核函数是一个典型的巨型稀疏矩阵, 目前有较多的求解算法, 例如联合迭代(SIRT)、最小平方QR分解(LSQR)及共轭梯度(CGNR)算法, 这些算法均可以稳定求解。但在实际应用过程中, 为了提高算法的收敛速度, 往往需要引入一些构造形态等的先验信息^[19-20], 进行适当的正则化约束, 使解朝着期望方向逼近。目前正则化方法中应用较多的主要有两种, 一种为贝叶斯框架下经典的Tikhonov正则化, 若引入此正则化, 那么求解方程(1)可得:

$ \Delta {\boldsymbol{s}} = {\left( {{{\boldsymbol{K}}^{\rm{T}}}{\boldsymbol{K}} + {{\boldsymbol{\Gamma }}^{\rm{T}}}{\boldsymbol{\Gamma }}} \right)^{ - 1}}{{\boldsymbol{K}}^{\rm{T}}}\Delta {\boldsymbol{t}} $

(3)

其中, Γ表示Tikhonov正则化项。对于公式(3)而言, 每次迭代需要求解的矩阵变为:

$ \left[ {\begin{array}{*{20}{l}} {\boldsymbol{K}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \end{array}} \right]\Delta {\boldsymbol{s}} = \left[ {\begin{array}{*{20}{c}} {\Delta \boldsymbol{t}}\\ 0 \end{array}} \right] $

(4)

显然, 这种正则化方式会增加矩阵维度, 内存占用量及计算量均随之增加。

相比较而言, 另一种预条件正则化将更加灵活, 增加正则化项可保持矩阵维度不变。假设Δs=FΔu, 其中, F表示预条件正则化算子。那么预条件正则化下方程(1)的解则表示为:

$ \Delta {\boldsymbol{u}} = {\left( {{{\boldsymbol{F}}^{\rm{T}}}{{\boldsymbol{K}}^{\rm{T}}}{\boldsymbol{KF}}} \right)^{ - 1}}{{\boldsymbol{F}}^{\rm{T}}}{{\boldsymbol{K}}^{\rm{T}}}\Delta {\boldsymbol{t}} $

(5)

其中, Δu表示预条件的解。

对于公式(5), 每次迭代需要求解的矩阵为:

$ [{\boldsymbol{KF}}]\Delta {\boldsymbol{u}} = [\Delta \boldsymbol{t}] $

(6)

显然, 预条件正则化未改变矩阵的维度, 仅增加了矩阵相乘的计算量。实际上两种算法在一定条件下是等价的, 采用相同的先验信息, 两者结果基本一致, 但预条件正则化方法计算效率更高。对于公式(6), 本文采用CGNR算法迭代求解^[21]。

2 基于最优路径寻优(OSV)的断控正则化技术

目前针对断裂属性的检测算法很多, 其中应用较为广泛的是相干属性及似然性属性, 这些算法能够较好地刻画断裂形态, 对地震资料解释有较大的指导意义。但这些属性往往对噪声较为敏感, 提取的属性信噪比较低, 无法直接应用于速度建模的正则化约束。基于此, 本文采用了基于最优路径寻优(OSV)的断层提取技术^[22], 该技术通过对常规断裂属性进行优化, 得到信噪比更高的断裂属性, 断点更加清晰, 路径寻优公式如下:

$ D[i, j] = \mathop {{\mathop{\rm argmax}\nolimits} }\limits_j \sum\limits_{i = 1}^N s [i, j] $

(7)

式中: D[i, j]是利用最优路径寻优算法得到的第i行第j列的断裂属性; s[i, j]表示第i行第j列的常规断裂属性。采用相干属性计算s[i, j]:

$ s[i, j] = \frac{{\sum\limits_{i - {M_1}}^{i + {M_1}} {{{\left( {\sum\limits_{j - {M_2}}^{j + {M_2}} f \left[ {i - {M_1}, j - {M_2}} \right]} \right)}^2}} }}{{\left( {2{M_2} + 1} \right)\sum\limits_{i - {M_1}j - {M_2}}^{i + {M_1}j + {M_2}} {{{\left( {f\left[ {i - {M_1}, j - {M_2}} \right]} \right)}^2}} }} $

(8)

其中, f[i－M₁, j－M₂]表示时窗(2M₁+1)×(2M₂+1)内的成像点, M₁和M₂分别是两个方向的半时窗长度。采用公式(7)和公式(8)即可确定沿断裂方向的最优路径, 从而得到更加干脆清晰的断裂属性。图 2a和图 2b分别显示了地震剖面和利用最优路径寻优算法提取的断裂属性。可以看出, 利用最优路径寻优算法提取的断裂属性断点清晰, 断面连续, 信噪比较高, 更加适用于速度建模的正则化约束。

利用上述方法得到断裂属性后, 即可构建预条件正则化算子F。采用基于构造导引的平滑算法^[23]引入断裂属性约束, 具体公式如下:

$ {{\boldsymbol{V}}_{{\rm{out }}}} - \alpha {\nabla ^{\rm{T}}}{\boldsymbol{D}}\nabla {{\boldsymbol{V}}_{{\rm{out}}}} = {{\boldsymbol{V}}_{{\rm{in}}}} $

(9)

式中: D为断裂属性; ▽为梯度算子。该公式相当于一个断裂预条件约束滤波器, 其中, V_in和V_out分别为滤波器的输入和输出。在本文中, V_in是CGNR算法求解过程中未进行断裂预条件约束的中间变量, V_out是CGNR算法求解过程中进行了断裂预条件约束后的中间变量, α是控制约束力度的参数, 当α=0, 即V_out=V_in时, 没有进行断裂属性约束。可以看出, 方程(9)实际上是对输入数据沿着断裂方向进行构造平滑, 从而达到断裂约束的目的, 这样可以高效稳定地引入断裂属性进行正则化约束。公式(9)可改写为:

$ {{\boldsymbol{V}}_{{\rm{out }}}}\left( {{\boldsymbol{I}} - \alpha {\nabla ^{\rm{T}}}{\boldsymbol{D}}\nabla } \right) = {{\boldsymbol{V}}_{{\rm{in}}}} $

(10)

即断控预条件算子F为:

$ {\boldsymbol{F}} = {\left( {{\boldsymbol{I}} - \alpha {\nabla ^{\rm{T}}}{\boldsymbol{D}}\nabla } \right)^{ - 1}} $

(11)

其中, I表示单位矩阵。对于公式(10)而言, 小规模的数值计算可以直接采用公式(11)的求逆计算, 但大规模生产应用时, 公式(10)往往是大型稀疏矩阵, 常规求解算法稳定性较差, 计算效率较低。可以采用Z变换来近似表示公式(10)中的偏微分算子, 进而采用有限差分算法进行稳定求解, 具体求解过程在文献[23]中已经有详细探讨, 这里不再赘述。将预条件算子代入公式(5)即可实现基于断裂属性约束的预条件正则化深度域速度建模。

3 基于断裂属性约束的深度域层析速度建模流程

基于断裂属性约束的深度域层析速度建模技术流程(图 3)与常规射线层析反演的流程基本一致, 只是在正则化项上引入了断裂属性约束。

从流程图可以看出, 高斯束层析的输入要求从叠前深度偏移剖面上自动拾取控制点, 并在这些控制点对应的成像道集上自动拾取剩余时差Δt, 再利用高斯束层析算子计算层析核函数K; 同时从偏移剖面中提取断裂属性D, 进而得到预条件算子F, 最终构建形成层析反演矩阵, 利用CGNR算法进行高效求解。

4 模型及实际资料测试 4.1 二维层状模型数据测试

为验证本文方法的精度, 设计了如图 4所示的二维正断层模型。模型横纵向采样点为302×641, 采样间隔10 m×10 m, 模型的速度范围为2 000~4 300 m/s, 断层倾向为左倾, 倾角为70°, 垂向断距为0~40 m, 采用一个更加接近实际地质情况的断裂带来刻画断层, 该断裂带横向宽度为40 m, 地震速度为3 200 m/s。采用中间激发、双边接收的观测系统进行声波方程正演, 共激发50炮, 炮间距为20 m, 每炮201道, 道间距为10 m, 正演采样率为1 ms, 最小偏移距为0。正演结果如图 5所示。从图 5中明显可见地层反射波及断面绕射波, 波场较为复杂。利用该数据进行断裂属性约束深度域速度建模, 建模结果如图 6所示。

图 6a至图 6d分别为理论地质模型、初始模型、常规网格层析反演结果和本文方法反演结果。由图 6可以看出, 常规网格层析反演结果整体构造趋势基本准确, 但对于高陡断裂的刻画精度较差; 断裂属性约束层析反演方法得到的速度模型精度更高, 能够清晰地看到断裂的影子, 与真实模型更加接近。

图 7显示了图 6a红线位置的单道速度对比结果。显然, 与初始模型相比, 网格层析反演结果在一定程度上提升了速度模型精度, 引入断裂属性约束的层析反演结果精度更高, 更接近真实模型, 断层刻画更为清晰。说明断裂属性约束的层析速度建模技术在断层刻画方面明显优于常规网格层析反演方法, 反演结果更加准确。

图 8a和图 8b分别为常规网格层析和断裂属性约束层析反演的速度模型对应的成像道集。图 8中抽取的6个道集对应图 4中红色三角所指位置。从图 8可以看出, 常规网格层析反演方法可以将成像道集两侧简单构造区域的道集拉平, 但中部断层附近的道集并没有完全拉平, 证明断裂周边速度精度还有待提升; 断裂属性约束层析反演速度模型能够将道集整体拉平, 证明速度场更加准确。另外, 图 8a中红色箭头位置是断面附近的成像道集, 可以明显看出断层对成像道集的影响。图 9a和图 9b分别显示了采用常规网格层析和断裂属性约束层析方法反演的速度模型对应的RTM成像结果。从图 9a中可见, 由于断层两侧速度不准确, 使得断层下盘明显存在构造假象, 成像质量较差。而图 9b中的成像结果更优, 断层下盘构造假象得到消除, 成像更加准确。

4.2 实际资料应用

实际资料来自中国东北某工区, 该工区地形比较平缓, 相对高差较小, 但该工区地下波场复杂, 存在多组小尺度断裂系统, 速度变化剧烈, 速度场求取难度大, 深度域精确速度建模困难。

图 10a和图 10b分别显示了常规网格层析和断裂属性约束层析反演方法速度建模结果。可以看出, 常规网格层析反演方法建模结果能够准确刻画大尺度构造, 但所得速度场较为平滑, 建模结果中多为速度的低波数分量, 缺乏细节信息; 断裂属性约束层析反演方法除了能够反演得到大尺度速度背景场外, 还能够反演得到速度的中高波数分量, 速度中多组断裂形态清晰可见, 更多细节信息得到展示, 速度模型精度明显提升。

为了更清晰地展示图 10a和图 10b的速度精度差异, 抽取CDP=100位置的单道速度进行对比, 如图 11所示。从图 11可以看出, 常规网格层析反演的速度趋势与测井速度大致吻合, 但精度较差。断裂属性约束层析反演速度更加接近测井曲线, 细节更加丰富, 模型精度更高。

图 12a和图 12b分别显示了常规网格层析、断裂属性约束层析反演的速度模型对应的RTM成像结果。可以看出, 对于大尺度构造成像, 两者结果相当, 但图 12a对细节的刻画能力较差, 特别是对于小断裂成像; 图 12b中的地层和断面归位更准确, 断面更加干脆, 断裂形成期次及交切关系更加清楚, 基底成像更加清晰。

图 13是图 12蓝框内的局部放大图。对比图 13a和图 13b可以看出, 断裂属性约束层析反演方法反演的速度模型的RTM成像结果对断层刻画更加合理, 断点清晰, 断面连续, 成像更加精确, 更有利于后续对小尺度断裂的精细解释。实际资料测试结果还表明, 断裂属性约束层析反演方法可以反演得到更加精细的速度模型, 对应的成像结果也更加精确合理。

5 结论

深度域速度建模是一项系统工程, 面向复杂断裂的速度建模在精度与合理性方面都面临着巨大挑战。本文提出了一种基于断裂属性约束的层析速度建模方法, 通过添加含有断层信息的正则化约束项, 实现了断裂属性约束速度建模。理论模型和实际资料测试结果表明, 本文方法建立的速度模型对复杂断裂系统刻画更加清晰, 反演的速度模型细节更加丰富, 成像品质整体显著提升, 这也证明了基于断裂属性约束的层析速度建模技术具有较大实际应用前景。