射线追踪方程与程函方程的初至旅行时层析对比

射线追踪方程与程函方程的初至旅行时层析对比

郭振波;孙鹏远;李培明;任晓乔;钱忠平;唐博文

【摘要】初至旅行时层析反演是目前应用最广泛的近地表建模方法.按照反演中正演算子的不同,可将基于射线理论的层析反演方法分为基于射线追踪方程与程函方程两类.本文分别从理论及数值测试上对两种方法在反演精度、计算效率等方面进行了系统的对比分析.结果表明:①两种方法可在统一的反演框架下推导得到,两者主要的差别均是由反演中正演算子的不同引起的;②两者具有相似的反演精度,由于后者的核函数是带限的,在复杂地区的层析反演中更加稳定;③前者的计算效率、内存占用等依赖于检波点个数,后者则依赖于模型大小,检波点稀疏时优选前者,否则优选后者;④前者具有射线密度等质控手段,后者缺少类似的质控手段.

【期刊名称】《石油地球物理勘探》

【年(卷),期】2019(054)003

【总页数】8页(P558-564,576)

【关键词】层析;射线追踪;程函方程;近地表建模;地震反演

【作者】郭振波;孙鹏远;李培明;任晓乔;钱忠平;唐博文

【作者单位】东方地球物理公司物探技术研究中心,河北涿州072751;东方地球物理公司物探技术研究中心,河北涿州072751;东方地球物理公司物探技术研究中心,河北涿州072751;东方地球物理公司物探技术研究中心,河北涿州072751;东方地球物理公司物探技术研究中心,河北涿州072751;东方地球物理公司物探技术研究中心,河北涿州072751

【正文语种】中文

【中图分类】P631

0 引言

对于油气地震勘探，需建立较精确的近地表模型以消除地表起伏及近地表速度异常对反射信号的影响，特别是对于山地等复杂探区，近地表模型的精度直接决定了最终的成像效果[1-2]。由于初至旅行时的稳健性及旅行时层析的高效性，目前利用初至旅行时层析反演进行近地表建模是应用最广泛的方法。近地表模型通常用于静校正[3]、深度域建模[4]等处理流程。

经典的初至层析反演方法基于射线追踪方程，在反演中需要显式计算炮点到检波点的旅行时及射线路径[5]，本文将这类经典的初至层析反演方法称为基于射线追踪方程的初至层析反演方法。反演中通常采用两种方式求取旅行时及射线路径：①采用两点射线追踪[6]或类似方法同时求取旅行时及射线路径；②采用波前追踪类方法，如基于程函方程的有限差分类方法[7-8]、基于线性旅行时假设的线性旅行时插值方法[9]、基于图理论的最短路径方法[10]等，该类方法需在求解得到初至时间场的基础上从检波点开始反向追踪得到射线路径[9，11-12]。反演通常采用反向投影[13]、最小二乘奇异值分解(LSQR)[11]、同时迭代重构(SIRT)[14]等方法。基于射线追踪方程的层析反演方法已在实际生产中广泛应用，目前主要研究方向是如何与波形反演等精度较高的反演方法相结合[15-16]。

除了上述经典算法外，还存在一类层析反演算法，本文将其称为基于程函方程的层析反演方法。该类方法仅需要计算初至时间场而不需要计算射线路径，借助伴随状态方法[17]求取模型更新量。该类方法最早由Sei等[18]于1994年提出，近年来得益于波形反演的兴起，伴随状态方法为地球物理学家所熟知，进而开展了广泛的

研究。Leung等[19]、Taillandier等[20]、谢春等[21]采用快速扫描方法进行初至旅行时层析； Huang等[22]利用该方法进行透射及反射旅行时层析；Benaichouche等[23]、李勇德等[24]讨论了该类层析方法中预处理方法对层析的影响； Waheed等[25]将该类方法用于各向异性介质的初至旅行时层析。

虽然两类方法已得到了广泛研究，但目前还没有从理论及数值测试上进行系统的对比分析。本文首先基于统一的反演框架从方法原理上对两种方法进行了详细的对比分析，并做了定性的评价；然后通过一个数值测试实例，对两种方法的反演精度、计算效率、内存占用等情况进行了定量的对比分析。

1 方法原理及对比分析

1.1 层析反演基本理论

旅行时层析通过迭代反演寻找正演旅行时与观测旅行时最优拟合的模型，若不考虑正则化项，L2模的目标函数可写为

(1)

式中： m为模型介质参数向量； ns为总炮数；为在当前介质模型m下第s炮对

应所有检波点处的正演初至时间；为在第s炮观测记录上通过拾取得到的初至时间。将目标函数在m0处展开且只保留至二阶项，可整理为

O(m)≈O(m0)+[|m=m0]Tδm+

|m=m0δm

(2)

通常将目标函数相对于介质参数的一阶导数称为梯度，即

(3)

将目标函数相对于介质参数的二阶导数称为Hessian矩阵，即

(4)

由于求取目标函数最小，需要求取目标函数的驻点，即g=0的点。对式(2)求导并将其设为零，可得

H|m0δm=-g|m0

(5)

δm=-(H|m0)-1g|m0

(6)

式(6)为在当前模型m0下对应模型更新量的显式表示。利用δm更新当前模型参数m0并将其作为下次迭代的初始模型，可实现非线性反演。当目标函数小于给定的阈值时，停止迭代，得到最终的反演模型。

由于Hessian矩阵的逆很难直接求取，通常通过迭代反演的方式求取式(5)所示的线性方程，此时反演方法通常称为高斯—牛顿法。对于尺度较大的问题，Hessian 矩阵H的求取非常耗时且存储困难，通常需要对其进行近似处理。若将Hessian 矩阵近似为一单位矩阵，则退化为最速梯度方法；若将Hessian矩阵近似为对角或带限矩阵，则退化为预处理的最速梯度方法。

1.2 基于射线追踪方程的初至层析反演

对基于射线追踪方程的初至层析反演，需要求取每个炮检对的射线路径，在得到射线路径之后，可将初至时间表示为

(7)

式中： nm为模型离散网格总数； mi为第i个网格内的慢度；为第s炮、第r个接收点对应射线在第i个网格内射线段的长度。可将其写成矩阵的形式