基于大数据量的初至层析成像算法优化

1.北京大学地球与空间科学学院, 北京 1000871; 2.中国石油青海油田公司, 敦煌 736200; †通信作者, E-mail: zxb@pku.edu.cn (张献兵), tianyue@pku.edu.cn (胡天跃)

摘要三维陆地地震资料采集技术向着高密度和宽方位的方向发展, 导致产生海量地震数据。初至旅行时层析成像的经典算法需要过多内存和计算时间, 不适合处理大的地震数据量。为解决这一问题, 提出一种基于优化数学公式的初至旅行时层析成像算法, 既减少 Frechet 矩阵和 Hessian 矩阵的存储对内存的需求, 也减少Hessian 矩阵的求逆时长, 能快速和高效地解决大数据量的层析反演问题, 适合高密度和宽方位采集的陆地地震资料和精细参数化的网格模型, 并且不影响精度, 易于实现并行计算。模型试验和实际数据都验证了该方法的有效性, 当初至旅行时达到一定的数量时, 该方法能提供可靠的层析成像结果, 进行静校正。

近年来, 中国西部地区逐渐成为油气勘探的重点区域。但是, 这一区域存在复杂的近地表结构, 包括剧烈起伏的地表、不稳定的高速层顶结构以及浅层的横向速度变化, 静校正问题突出。目前的折射静校正方法要求存在稳定的高速层顶界面, 在近地表结构复杂的中国西部地区应用效果欠佳。采用基于初至旅行时的层析成像, 可有效地获得浅层速度结构, 并计算静校正量[1‒2], 在中国西部复杂地表工区的应用效果证明, 与折射静校正相比, 层析静校正可以得到更准确的静校正量[3‒4]。

近地表速度建模是地震资料初至层析静校处理的重要基础, 建模质量直接决定静校正量的准确性。速度建模包含通过求解程函方程实现射线追踪计算最短路径的正演过程, 以及获取 Frechet 导数矩阵和 Hessian 矩阵, 并通过层析反演得到速度模型更新的反演过程。其中, 正演过程是通过给定的初始速度模型计算理论走时。早期的试射法和弯曲法等射线追踪方法[5]存在多路径问题, 可能导致反演结果不稳定, 随后流行的有限差分法[6]和快速扫描法[7‒8]仅适用于速度横向变化的模型。近来发展的后验最短路径射线追踪法, 计算精度高且计算结果稳定, 广泛应用于层析成像正演计算中[9‒10]。现在的地震数据采集系统部署数以千计的震源和检波器, 产生百万级甚至更多的采集道数, 同时勘探工区面积的扩大导致速度模型包含千万个网格数量[11]。此外, 复杂近地表结构往往包含强烈的垂向速度变化, 要求非常密集的垂向网格划分。对于二维观测系统, 经典的初至旅行时层析成像算法一般不受计算能力的限制, 但对于三维观测系统, 则要考虑计算的时间消耗和内存需求量。

解决上述限制的方法包括减少初至拾取的数量或离散化速度模型, 但都会造成有效信息的缺失, 进而导致速度模型的精度和分辨率降低。如何减小或解除 Frechet 导数矩阵和 Hessian 矩阵对内存的需求, 是亟待解决的关键问题之一。与此同时, 反演问题中 Hessian 矩阵的大小是网格参数的平方, 对于三维大数据量的采集, Hessian 矩阵十分庞大, 难以直接解析求逆, 是非常棘手的问题。一般的折射旅行时层析静校正方法是结合基于后验射线追踪的程函方程来计算 Frechet 导数矩阵, 然后用拟牛顿优化方案(如 LSQR 算法[12])来解决线性优化的反演问题。经典的拟牛顿法不计算真实的 Hessian 矩阵, 而是通过迭代方法, 不断更新对 Hessian 矩阵的估计。最常用且有效的拟牛顿法是 BFGS 方法。为减小记录 Hessian 矩阵的内存消耗, Nocedal 等[13]提出Limited BFGS (L-BFGS)法, 只记录模型变化量和目标函数梯度这两组向量, 用于计算每个步骤的 Hes-sian 矩阵。另一种减小层析反演内存消耗的方法是Taillandier 等[14]提出的伴随状态法, 通过引入伴随矩阵来计算目标函数的梯度, 不需要计算 Frechet导数矩阵。谢春等[15]随后发展了伴随状态走时方法, 验证了其有效性。

本文提出一种新的基于优化数学公式的初至旅行时层析成像算法, 可以减少内存消耗和计算时长, 获得满意的层析成像结果。该方法基于后验的射线追踪, 直接求取速度模型, 不拟合方程梯度的局部最优值和它对应的 Hessian 矩阵, 估算的元素被耦合到拟牛顿法中, 用以确保目标函数的最小化。这样直接求取 Frechet 矩阵和 Hessian 矩阵的主要优势是能够快速和高效地处理大数据量的反演问题。该方法适合于高密度采集和精细参数化的网格模型, 并易于实现并行计算。经合成数据实例验证, 该方法可以得到相对准确的速度建模结果。将本文提出的方法应用于处理柴达木盆地的三维实际地震资料, 工区初至数量达到百万数量级, 同样能高效地实现层析建模, 且与高程静校正和折射静校正相比, 炮集和叠加剖面中的同相轴连续性明显加强。

1 方法原理

1.1 简化的求逆方法

初至波层析建模是一个反演问题, 即利用观测的初至波到时反演得到近地表速度模型。如果可以找到使计算得到的初至旅行时与观测值最接近的速度模型, 就认为已得到正确的速度模型。因此, 构建如下具有最小二乘法意义的目标函数:

式中, s是慢度模型, tobs(s)是观测的初至旅行时, t(s)是根据慢度模型正演计算得到的初至旅行时。当t(s)与tobs(s)相等时, J(s)达到最小值, 因此, 初至旅行时层析反演就转化为最小化目标函数J(s)的计算问题。

为求解反演问题, 首先要解决由慢度模型得到计算旅行时的正演问题。采用一维有限差分程函方程正演计算旅行时t(s):

式(2)和(3)中, x是空间坐标, xshot 是炮点位置。式(3)规定炮点处的旅行时为零, 该算法具有快速并支持强速度差异的特点。通过局部牛顿迭代方法求解使目标函数 J 最小化的模型:

式中, ∇J 是目标函数对慢度的梯度变化值, H–1是相关 Hessian 矩阵的逆, n 表示迭代次数。牛顿方法的特点是, 当目标函数是二次方程形式, 即需要求解的问题是线性问题时, 仅进行一次迭代即可收敛。因此, 求解每次迭代的 H–1 和∇J 就可以得到速度模型。

考虑近地表的一个笛卡尔网格的慢度模型 s=sj, 第 i 个初至旅行时 ti 与慢度模型 sj 的关系可以由每个离散网格内的射线长度lij相联系:

其矩阵形式为 t=Ls。根据这种线性化近似的正演问题, 将之前定义的目标函数(式(1))改写为

基于上述近似的正演方法, 用有限差分法求解目标函数的一阶偏导数∇J:

由以上计算可得, 目标函数梯度为∇J=LTδt, 其中 δt = t – tobs 是旅行时剩余量。用有限差分法求解目标函数的二阶偏导数 H:

由以上计算可得, 目标函数的步长 Hessian 矩阵为 H=LTL, 所以这里的 Hessian 矩阵是对称矩阵。射线长度 L 组成的矩阵对应于 Frechet 导数矩阵, L的大小等于观测旅行时数量乘上模型参数的数量, H 的大小只与模型参数数量有关, 因此式(4)的牛顿方程迭代过程改写为

其中, δsn = sn+1 – sn 为慢度模型的更新量。

每一次迭代过程都要求计算目标函数梯度∇J和 Hessian 矩阵 H 的逆, 为避免存储整个 Frechet 和Hessian 矩阵, 减小内存消耗, 本文采取如下局部化估计方法。

1)估计单炮的∇J: 对每组炮检点组合, ∇J是网格内的射线长度与剩余旅行时的乘积, 全局的∇J是对所有单炮估计的∇J求和。

2)在牛顿迭代方法中, 目标函数梯度对应目标函数减小的方向, 目标函数梯度对应的Hessian矩阵的逆给出步长, 简化的Hessian矩阵的逆的计算方法如下:

其中, 对角元素为Hjj=(LTL)jj=∑ilij2, 且大于非对角元素。因此, Hessian矩阵的逆可以被近似为一个对角矩阵, 其对角项可写成

1.2 简化求逆方法正确性验证

为证明用只含对角线元素的近似Hessian矩阵代替Hessian矩阵逆的正确性, 本文设计数值模型来测试和对比。

图1显示一组实验结果, 图1(a)中射线长度矩阵L由随机稀疏矩阵给出, 满足初至波射线仅穿过模型中少量网格的条件; 图1(b)是所构成的Hessian矩阵H=LTL, 它的对角元素大于非对角元素, 该实例中对角元素是非对角元素的22倍; 图1(c)是Hessian矩阵的逆Hinv_accu, 由摩尔‒彭若斯广义逆算出; 图1(d)是利用式(8)的优化方法估计的Hes-sian矩阵的逆Hinv_sim, 由 width=54.8,height=17.25

算得出。对Hinv_accu和Hinv_sim, 用下式求相似度:

式中,

和

分别为A和B的矩阵元素平均值。当输入初至旅行时为2000, 网格模型大小为90×90时, r=0.9780。也就是说, Hinv_accu和Hinv_sim的对角元素几乎相同。两者对角元素差值的二范数占Hessian矩阵逆对角元素二范数的比率为0.079。进一步得到的模型慢度修正量中, 剩余旅行时差由随机数生成, 两条曲线基本上重合, 由此验证优化的求逆方法在该实例中可以得到相对准确的慢度修正量。

表 1 列出改变输入初至数量和网格模型大小的多组试验数据结果。实验数据显示, 对于同样大小的网格模型, 输入的初至数量越大, 近似Hessian矩阵的逆越接近准确的 Hessian 矩阵的逆。对于inline 和 crossline 方向上大小为 90×90 的网格, 当初至输入量达到 30 万及以上时, 误差仅有1/10000。与直接求逆算法相比, 优化的逆矩阵算法大大地缩短了计算时间(表 1 第 7 列所示)。在本文的局部化估计过程中, 计算得出每组初至对应的 Frechet 矩阵后, 可以直接计算局部∇J 和 Hessian 矩阵的逆, 并直接累加到全局结果中, 不需要在内存中存储 Fre-chet 矩阵。因此, 本文计算方法要求的内存消耗只取决于离散化模型的大小, 与输入的初至数量无关, 大大地节约了内存(表 1 第 8 列所示), 适用于海量地震初至的层析反演。

1.3 算法流程

本文层析反演算法的总流程如下。

1)拾取全部输入道集的初至到时。静校正问题由起伏地表和近地表速度变化引起, 这类复杂近地表结构往往带来采集施工难度和各类噪声, 导致原始数据的初至信噪比很低。为保证初至拾取的准确性, 往往需要在计算机自动拾取初至的基础上手动拾取。针对海量地震数据, 可能包含百万道地震道, 手动拾取将耗费大量的时间和人力成本。因此, 对输入数据做适当的预处理(如滤波法或地震干涉法)后再应用自动初至拾取算法, 可以高效且准确地拾取初至[16‒17]。本文采用相似度加权的超虚干涉法来加强初至波信号[18]。

2)建立初始速度(或慢度)模型。初始模型越精确, 反演收敛所需的迭代次数越少。本文假设近地表模型具有层状结构, 采用折射静校正的层状模型作为初始模型, 利用拾取的初至到时, 线性地反演每一层的层厚和速度, 得到初始速度模型。

3)利用当前速度模型进行射线追踪。通过求解一维有限差分程函方程(式(5)和(7)), 得到单炮的Frechet 矩阵和旅行时残差。计算单炮的目标函数梯度∇J, 并累加到总∇J中, 利用优化的Hessian矩阵求逆方法求解该炮对应的Hessian矩阵逆(式(11)),并累加到总的Hessian矩阵逆中。

4)对所有炮的输入数据重复步骤3, 得到总的∇J和总的Hessian矩阵逆, 利用式(9)求解模型更新量 δsn, 并更新速度模型。

5)重复步骤3和4, 直至δsn小于设定的阈值或达到最大迭代次数, 得到最终反演的速度模型。

2 合成数据实例

在模型实验中, 由于小范围数据采集时, 可将大部分地下介质视为层状介质, 所以本文首先选择水平层模型进行实验, 用于验证本文方法的正确性。我们设置一个大小为0.5×4km2的起伏层模型, 如图2所示, 模型介质的声波传播速度为一个四层介质, 其中有两个凹陷, 每层的速度随着深度依次递增, 自上而下的速度分别为2300, 3000, 3600和4200m/s, 密度取值为2000kg/m3。震源位于地表, 震源间距为200m, 检波器同样依次排列于地表, 检波器间距为10m。

计算过程中, 首先根据观测系统的设置, 对地震数据的初至走时进行拾取。以地震数据的初至走时为基础, 手动生成初始模型(图3, 其中第一层是一个地表低速层, 速度为1000m/s, 真实地层为其下的3层)。根据初始模型, 用多层模型的线性化反演法(layered model update by linearized inversion, LMI)计算折射静校正, 经过迭代后得到图3所示的结果。图3(b)的折射静校正结果说明, 本文建立的凹陷模型的反演结果在速度方面有所改进。

将折射静校正的结果作为本文算法的初始模型, 经过平滑处理后, 再采用本文算法进行反演计算, 最终得到图4(b)所示的层析成像结果。可以看出, 本文算法的计算结果更加稳定, 且平滑模型结果更加趋于真实模型, 证明了本文算法的稳定性和有效性。

我们截取图2所示起伏层模型的凹陷部分为对象(图5(a)), 研究初至层析反演算法对初始模型的依赖性。本研究进行以下4组试验。

1)根据该区块的多次采集数据等先验信息, 给出比较准确的初始速度模型(图5(b)), 其浅层速度为2300m/s, 梯度为4m/s², 速度随深度增加而增大。8次迭代后得到的反演模型(图5(c))中, 第二层凹陷刻画明显, 浅层刻画不准确。

2)在区域内选取多个控制点, 每个控制点给定合适的半径, 根据统计的初至时间, 建立初始速度模型(图5(d)), 8次迭代后得到的反演模型(图5(e))与真实模型吻合程度最好, 初至时间的均方根误差为3.2 ms。

3)对工区内已有的速度模型进行平滑处理, 作为初始速度模型(图5(f)), 8次迭代后得到的反演模型(图5(g))快速收敛, 迭代稳定, 印证了第2组试验结果的准确性。

4)假设在对介质模型没有调查的情况下给出浅层速度为3000m/s, 梯度为2m/s², 速度随深度增加而增大的初始速度模型(图5(h))经过8次迭代后得到错误的反演模型(图5(i)), 反演结果不稳定(如图6中灰色虚线所示), 证明初至层析反演要求有比较准确的初始模型。

4种情况下层析反演的收敛速度以及最终初至拟合的均方根值见图6, 可以得到以下认识: 由试验4得出稳定和准确的层析反演结果需要可靠的初始速度模型, 由试验1～3的最终收敛情况可知, 由于层析反演的是背景波场, 属于低波数成分, 所以根据先验信息和初至时间建立的初始速度场, 一般需要经过平滑处理, 才能作为最终反演的初始速度模型。

3 实际数据实例

为了衡量本文提出的初至旅行时层析成像算法的效果, 我们将它应用于实际数据。资料源于柴达木盆地西部的英东三维工区, 该工区的地面海拔高度在 2800～3500m 之间, 地形起伏大, 沟壑纵横, 地表呈东、西高, 中间低的马鞍型。工区南部比较平缓, 地表被巨厚砾石层覆盖, 主要障碍物为 315 国道、花格输油管线、大乌斯输油泵站、通信光缆和高压输电线路等。工区中部和北部为山地, 地表出露灰黄色与浅灰色砂泥岩互层的上第三系地层, 地表风化和剥蚀严重, 第四系风化黄土层的厚度为10～50cm, 干燥疏松, 冲沟发育, 冲沟走向的展布十分杂乱。该工区的施工难度创柴达木盆地之最。

布设于该工区的三维陆地地震资料采集观测系统有数以万计的炮点和检波点, 完全自动拾取的初至旅行时数据有几百万条。我们首先用LMI方法建立初始层状速度模型, 然后将层状模型转化为网格模型, 然后用本文算法的初始速度模型进行迭代更新。初始速度模型与用两种方法更新后的速度模型的对比如图7所示。

接下来, 基于上述速度模型进行层析反演和地表一致性折射静校正。这些静校正结果最终会应用到地震道, 然后通过求取叠加来对比图7中两种方法得到的速度模型的准确性。单炮和叠加剖面的对比如图8～11所示, 可以看出, 在应用高程+层析静校正后, 单炮初值变得圆滑, 扭曲变形的同相轴得到改善, 双曲线更明显; 叠加剖面同相轴的连续性增加, 高陡山地的复杂构造成像质量明显改善, 清晰度和连续性有很好的改进。用本文初至旅行时层析成像算法推断出的速度模型, 提供了比较好的折射静校正量。

4 结论

本文提出一种基于优化的Hessian矩阵求逆算法的初至旅行时层析成像方法, 该算法具有以下优势: 1)减少了Frechet矩阵和Hessian矩阵对内存的需求; 2)实现并行运算, 提高了运算效率; 3)当输入的初至旅行时达到一定数量时, 能保证层析成像的精度。

合成数据和实际资料都验证了算法的运算效率, 将层析成像的结果应用于折射静校正, 得到较好的结果。同时, 本文论证了层析反演需要可靠的初始模型。

基于优化数学公式的初至旅行时层析成像算法能有效地解决海量地震数据层析反演问题, 针对高密度、宽方位的三维陆地地震资料采集资料以及精细化的网格模型, 能快速地得到较为可靠的层析成像结果, 解决长波长静校正问题。

[1] Docherty P. Solving for the thickness and velocity of the weathering layer using 2-D refraction tomography. Geophysics, 1992, 57(10): 1307‒1318

[2] 韩晓丽, 杨长春, 麻三怀, 等. 复杂山区初至波层析反演静校正. 地球物理学进展, 2008, 23(2): 475‒483

[3] 张林, 杨勤勇, 张兵, 等. 复杂近地表初至波层析反演静校正技术研究. 地球物理学进展, 2017, 32 (2): 0816‒0821

[4] 周衍, 饶莹. 黄土塬覆盖区的层析反演静校正方法研究. 地球物理学报, 2019, 62(11): 4393‒4400

[5] Zelt C A, Barton P J. Three-dimensional seismic refraction tomography: a comparison of two methods applied to data from the Faeroe Basin. J Geophysics Res, 1998, 103(B4): 7187–7210

[6] Vidale J E. Finite-difference calculation of traveltime. Bulletin of the Seismological Society of America, 1990, 78(6): 2062‒2076

[7] Zhao Y S, Gu H M, Shi X M. SIRT tomography inversion of 2D sonic data based on Moser’s curved ray tracing method. Chinese Journal of Engineering Geophysics, 2005, 2(3): 167‒176

[8] Qian J L, Zhang Y T. Fast weeping methods for eikonal equitions on triangular meshes. SIAM Journal on Numerical Analysis, 2007, 45(1): 81‒107

[9] Zhang J, Toksoz M N. Nonlinear refraction traveltime tomography. Geophysics, 1998, 63(5): 1496‒1823

[10] 黄国娇, 孙江兵, 白超英, 等. 三维 TI 介质中多波走时层析成像. 石油地球物理勘探, 2018, 53(1): 63‒72

[11] 王喜双, 赵邦六, 董世泰, 等. 油气工业地震勘探大数据面临的挑战及对策. 中国石油勘探, 2014, 19(4): 43‒47

[12] Paige C, Saunders M A. LSQR: sparse linear equa-tions and least squares problems, Part I and Part II. ACM Trans Math Soft, 1982, 8(1): 43–71

[13] Nocedal J, Wright S J. Numerical optimization. New York: Springer, 2000

[14] Taillandier C, Noble M, Chauris H, et al. First arrival traveltime tomography based on the adjoint-state me-thod. Geophysics, 2009, 74(6): WCB57–WCB66

[15] 谢春, 刘玉柱, 董良国, 等. 伴随状态法初至波走时层析. 石油地球物理勘探, 2014, 49(5): 877‒883

[16] Alshuhail A, Aldawood A, Hanafy S. Application of super-virtual seismic refraction interferometry to en-hance first arrivals: a case study from Saudi Arabia. The Leading Edge, 2012, 31(1): 34‒39

[17] An Shengpei, Hu Tianyue, Liu Yimou, et al. Automa-tic first-arrival picking based on extended super-virtual interferometry with quality control procedure. Exploration Geophysics, 2017, 48: 124–130

[18] 吕雪梅, 安圣培, 胡天跃, 等. 相似度加权的超虚干涉法加强初至波信号. 北京大学学报(自然科学版), 2018, 54(1): 87‒93

Algorithm Optimization of First-Break Tomography Statics Based on Large Datasets

1. School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Qinghai Oilfield, PetroChina, Dunhuang 736200; † Corresponding authors, E-mail: zxb@pku.edu.cn (ZHANG Xianbing), tianyue@pku.edu.cn (HU Tianyue)

Abstract The development of 3D land seismic data acquisition in the direction of wide azimuth and high density will lead to huge dataset. The classical first arrival time tomography algorithms are not suit for processing huge seismic dataset due to very high memory request and computing time cost. In order to solve this problem, the authors develop an optimal mathematical formula from the classical first break travel time tomography method to avoid the memory occupation that required by Frechet derivative matrix and Hessian matrix, and reduce the time cost of computing Hessian matrix inversion. This method can efficiently solve the tomography inversion problem for huge datasets. It is suitable for huge and high-density land seismic exploration, and not affect the dataset and model accuracy. It is easy for parallel processing. Both the model and real data examples confirm the effectiveness of this method. It can provide reliable tomographic results for static correction when the first breaks reaches a certain amount.

国家重点研发计划(2018YFA0702503)、国家自然科学基金(41674122)和国家科技重大专项(2016ZX05004003)资助