三维初至波旅行时层析速度反演算法优化

0 引言 中国西北部塔里木盆地蕴含着丰富的油气资源，勘探潜力巨大，但地表结构复杂，地貌以山前带、戈壁、沙漠等为主。盆地内地表起伏剧烈、高速砾岩分布不均匀、岩层风化剥蚀严重，如此复杂的近地表结构导致地震波传播路径非常复杂，难以得到准确的近地表速度结构，对后续静校正、中深层反射波层析反演及偏移成像等产生极大影响。为克服这些问题，高精度静校正[1-3]、起伏地表射线追踪、基于波动方程反演成像[4-6]等处理技术应运而生，而通过复杂地表条件下初至波走时层析反演获得高精度的浅表层速度则是这些技术应用的基本保障。初至波层析反演对准确描述浅表层速度结构具有重要意义，同时也是地震数据处理和成像反演中至关重要的一步[7-8]。初至波旅行时是初至波走时层析反演中关键的信息，初至拾取的精度是获得高精度走时层析反演结果的前提[9-10]。目前，国内外已提出很多初至拾取方法，其中采用图像分割方法进行初至拾取是深度学习在地球物理学应用方面的热点[11-12]。相较以往仅将地震剖面分割为初至前与初至后的分类方法，将初至拾取转化为图像分割问题来解决，新增了包含初至窄带的一类图像，并对这3类分割图像赋予了相应权重以突出包含初至窄带的像素点。在正演走时计算方面，Vidale[13-14]提出了有限差分算法计算地震波旅行时，但存在计算不稳定问题；Sethian等[15-17]提出了快速推进算法（Fast Marching Method，FMM），该方法基于迎风差分方案和快速排序技术数值求解非线性Eikonal方程，可应用于地震初至波走时计算，解决了Vidale方法的计算不稳定问题；Zhao[18]提出了快速扫描法（Fast Sweeping Method，FSM），该方法较FMM适用范围更广且计算效率更高；Hassouna等[19]提出了改进版本的FMM算法，又称为多模板快速推进算法（Multi-stencils Fast Marching Method，MSFM），通过旋转坐标系建立多个差分模板，并沿多个模板求解Eikonal方程，计算每个网格节点的旅行时，最后选择满足物理意义的解；王飞等[20]使用MSFM计算波前传播时间，采用最速下降法反向计算射线路径，以此提高旅行时计算的精度和效率；丁鹏程等[21-22]将MSFM成功应用于三维起伏地表条件下的走时计算。 现阶段“两宽一高”采集技术已常规化，获得海量的实际地震数据，使得初至拾取的精度与效率、射线类层析反演需存储的大型Hessian矩阵成为射线类层析反演发展面临的亟待解决的难题，因此，初至走时层析速度反演在效率与精度方面不得不采取折中的处理手段。为此，部分专家致力于快速、高效地解决大数据量的层析反演问题[23]；有的学者为解决西部浅表层地质情况复杂的问题，着力于提升层析反演的精度[24-26]；还有的专家针对初至波层析反演的某一核心环节进行改进[27-30]。在提升计算效率方面，采用分布式存储器编程MPI（基于信息传递的编程模型），结合MSFM技术与窄带延拓策略计算正演模型走时，并采用Runge-kutta算法完成反向射线追踪。在提升层析反演精度方面，保证初至拾取效率与精度的同时，针对初至层析反演流程中的核心环节，提出一套多信息约束的初至波层析反演优化策略，即增加视慢度与正则化对层析反演目标函数进行约束，降低目标方程的病态性、多解性，使用微测井信息对层析反演初始速度场加以约束。 1 机器学习初至拾取 此次研究采用的神经网络是基于U-Net（图 1）。正向传播过程中，卷积核的大小为（9，3），零填充大小为（4，1），步长为1；降采样过程采用二维平均池进行传播，内核大小为2，步幅为2；反向传播过程中，通过转置卷积将通道数量减半进行升采样，卷积核的大小为（7，3），零填充大小为（3，1），步长为2。由于特征图的大小没有变化（1 000×24），因此通过复制，将特征从收缩转移到扩展路径。从训练结束到初至拾取分为3类映射，通过二维卷积进行，内核大小为1，步长为1。 下载原图 图 1 神经网络架构 Fig. 1 Neural network architecture diagram 在开始训练之前，将包含45 000个模型的数据集按所必需比例（8∶1∶1）拆分为训练、验证和测试子集，每隔100个训练批次进行验证。在训练神经网络时，采用Adam作优化器，设置恒定的学习率为0.000 1。为及时停止学习过程，在验证数据集上计算IoU（交并比）指数，如果IoU指数在一定数量的验证中没有增大，则神经网络学习将停止。将上述方法应用于预测包含坏道的地震数据初至拾取（图 2），其中蓝线为处理人员手动拾取的真实初至，绿线为通过神经网络预测拾取的初至，对比2种初至拾取结果可以看出，在随机噪声以及地震坏道出现的情况下，通过神经网络预测拾取初至的精度可满足工业生产需求，精度与人工拾取无异。在拾取效率方面，对包含80% 随机噪声的576 000道的合成地震数据训练神经网络需要约3 h，通过神经网络预测拾取初至耗时仅约5 min，而处理人员人工手动逐道拾取则耗时约10 d，效率提升近80倍。另外，存在人工无法拾取的连续坏道的情况下（图 2a的第10~12道；图 2b的第18~20道；图 2c的第13~ 15道；图 2d的第15~17道），神经网络预测拾取初至时，通过相邻地震道拾取结果在坏道上插值拾取，仍可满足精度需求。 下载原图 图 2 机器学习预测地震初至拾取与人工拾取效果对比 Fig. 2 Comparison between machine learning prediction and artificial first break picking 2 三维初至波正演走时计算 2.1 多模板快速推进算法（MSFM） 多模板快速推进算法通过求解沿多个模板的Eikonal方程计算网格单元节点的走时，从中选择满足具有真实物理意义的解。在二维情况下，2个模板可以覆盖8个相邻节点；在三维空间中，6个迎风差分模板可以覆盖26个相邻网格节点。 以二阶差分格式下的模板F2为例（图 3），U1，U2，U3可以近似为 $ U_v=\max \left[\frac{3}{2\left\|x-x_v\right\|}\left(T_{i, j}^k-T_v\right), 0\right], v=1, 2, 3 $ (1) 下载原图 图 3 6个三维多模板快速推进算法（MSFM）笛卡尔域模板（据文献[19]修改） 注：T1，T2，T3分别为每个模板不同模板臂方向的最小走时。 Fig. 3 Six three-dimensional MSFM Cartesian domain templates 式中：xv为对角线方向较小走时的网格点；U1，U2，U3分别为沿着$\overrightarrow{r_1}, \overrightarrow{r_2} \text { 和 } \overrightarrow{r_3}$的方向导数；$\overrightarrow{r_1}=\left[r_{11} r_{12} r_{13}\right]^{\mathrm{T}}$, $\overrightarrow{r_2}=\left[r_{21} r_{22} r_{23}\right]^{\mathrm{T}}, \overrightarrow{r_3}=\left[r_{31} r_{32} r_{33}\right]^{\mathrm{T}}$分别为沿着l1l2，p1p2，q1q2方向的单位向量；l1l2，p1 p2，q1q2为网格相交点；Tv为对应网格点位置的走时，ms；T1，T2和T3分别为每个模板不同模板臂方向的最小走时，ms，可表示为 $ \begin{aligned} & T_1=\min \left(\frac{4 T_{i-1, j}^k-T_{i-2, j}^k}{3}, \frac{4 T_{i+1, j}^k-T_{i+2, j}^k}{3}\right) \\ & T_2=\min \left(\frac{4 T_{i, j+1}^{k-1}-T_{i, j+2}^{k-2}}{3}, \frac{4 T_{i, j-1}^{k+1}-T_{i, j-2}^{k+2}}{3}\right) \\ & T_3=\min \left(\frac{4 T_{i, j+1}^{k+1}-T_{i, j+2}^{k+2}}{3}, \frac{4 T_{i, j-1}^{k-1}-T_{i, j-2}^{k-2}}{3}\right) \end{aligned} $ (2) 由于多模板快速推进算法通过多个笛卡尔域模板计算网格节点的走时，考虑了网格节点周围26个相邻节点的走时信息，相较于单模板在对角线方向可更有效提高走时计算的精度。 2.2 起伏地表窄带延拓算法 通过多模板迎风差分算法仅可实现局部网格节点走时计算，采用三维窄带延拓技术模拟地震波波前扩展过程可实现全局走时计算（图 4）。根据波前到达时节点是否完成走时计算可将网格节点分为4类：①完成点，走时已计算且无法修改；②远离时还未计算；③窄带点，走时已计算，仍可更新；④地表以上点，位于地表之上，不参与波前扩展。 下载原图 图 4 三维网格节点旅行时计算状态 Fig. 4 Three-dimensional calculation status of grid nodes 最初，所有网格边界点都被标记为已知。通过求解计算它们的旅行时，将它们最近的网格节点标记为窄带。具体实现流程如下： ①在所有窄带点中，找到旅行时最小的点，将其状态标记变更为已知点。 ②找出距离最小旅行时节点最近的窄带点或远离点。 ③求解更新所有节点的旅行时间，即远离点标记为窄带点或窄带点被赋予更小的旅行时。 ④回到第一步循环迭代。 2.3 三维射线追踪 走时层析反演过程中需要进行射线追踪，即采用走时对模型慢度的微商来构建层析反演的敏感核函数，通过MSFM计算全局旅行时之后，再将全局走时计算结果反投影即可完成射线追踪。在各向同性介质条件下，地震波传播方向为波前面法线方向，即走时梯度方向。因此，可从MSFM得到的正演全局走时入手，以检波点的位置作为初始点，沿着走时梯度方向，按一定的步长反向进行射线追踪直至到达激发点，求取射线路径。 $ x_{n+1}=x_n+\alpha \nabla t $ (3) 式中：α为射线追踪的步长；$\nabla t$为走时梯度，ms，通过MSFM算法可获得；xn为射线追踪的当前位置；xn + 1为下一个射线追踪点的位置。这里以x方向为例（y，z方向计算方式同理），采用差分公式计算旅行时梯度$\nabla t$的一阶分量 $ \nabla t_x=\frac{\left(t_{i, j, \mathrm{k}}-t_{i-1, j, \mathrm{k}}\right)}{\Delta x} $ (4) 式中：$\nabla t$为坐标轴x方向上的走时梯度，ms；Δx为x方向网格间距大小，m；ti，j，k表示网格节点（i，j，k）处的旅行时，ms。 考虑到一阶分量的计算精度较低，为提升射线追踪的精度，使其适用于较复杂介质情况，推导给出了旅行时梯度的二阶与三阶分量作为参考： 二阶分量， $ \nabla_{t_x}=\frac{\left(3 t_{i, j, k}-4 t_{i-1, j, k}+t_{i-2, j, k}\right)}{2 \Delta x} $ (5) 三阶分量， $ \nabla_{t_x}=\frac{\left(11 t_{i, j, k}-18 t_{i-1, j, k}+9 t_{i-2, j, k}-2_{t-3, j, k}\right)}{6 \Delta x} $ (6) 通过上式可获得所有网格节点的走时梯度。为进一步提升全局走时梯度的计算精度，引入Rungekutta算法，通过近似射线方程的Taylor展开，使用多个网格节点的走时梯度加权获得当前射线追踪位置处xn的梯度值 $ x_{n+1}=x_n+\frac{h}{6}\left(K_1+2 K_2+2 K_3+K_4\right) $ (7) 式中：h为计算步长；K1，K2，K3，K4可分别表示为 $ \begin{aligned} & K_1=\nabla_t\left(x_n\right) \\ & K_2=\nabla_t\left(x_n+\frac{h}{2} K_1\right) \\ & K_3=\nabla_t\left(x_n+\frac{h}{2} K_2\right) \\ & K_4=\nabla_t\left(x_n+h K_3\right) \end{aligned} $ (8) 通过联立式（7）与式（8）便可完成走时梯度计算。设定射线追踪的终止阈值为射线追踪路径长度远小于射线追踪步长，即满足 $ \left\|x_{n+1}-x_n\right\|<\varepsilon \cdot h $ (9) 则射线追踪结束，式中ε为一常数，根据不同精度要求可取不同值。 结合不同阶数旅行时梯度计算公式、RungeKutta算法及步长h可灵活调节的优势，该方法可以满足复杂地表三维实际资料层析反演的要求，并且兼具计算效率与稳定性高的优势。 3 基于多信息约束的三维初至波层析反演 3.1 视慢度约束 走时层析反演通过拟合真实走时d与正演走时G（m），逐步更新慢度模型。其反演过程可表示为 $ \varPhi(m)=\|\bar{d}-\bar{G}(m)\|^2 $ (10) 式中：m为慢度，s/m；Φ（m）为反演目标函数；d = d/l，表示对真实初至走时的归一化，s；G (m) = G(m) /l表示对模型正演走时的归一化，s；l为射线路径长度，m。 使用归一化处理的目的是防止在走时较大的情况下反演过程中走时残差收敛太快，而此时浅表层速度尚未准确反演，从而影响浅层速度模型的反演效果。对于包含海量走时数据的三维走时层析反演而言，模型网格较多，射线数目的不足会导致层析反演方程组欠定，加剧反演问题的病态性。因此，需加入先验信息来约束反演解空间，以加速迭代收敛，方可获得满足实际应用需求的近似解（最大似然解）。针对这一问题，Zhang等[31]提出将视慢度信息加入目标函数以约束慢度模型反演平均慢度d，有助于提高浅层速度的精度，反演视慢度则有助于恢复深层速度结构。因此，通过最小化平均慢度与视慢度的不适配来降低反演多解性。在引入视慢度信息后，目标函数Φ（m）可表示为 $ \varPhi(m)=(1-\omega)\|\bar{d}-\bar{G}(m)\|^2+\omega\left\|D_x[d-G(m)]\right\|^2 $ (11) $ \Phi(m)=(1-\omega)\|\bar{d}-\bar{G}(m)\|^2+\omega\|\hat{d}-\hat{G}(m)\|^2 $ (12) 式中：Dx是旅行时相对于距离的微分算子；$\hat{d}=D_x d$为走时与偏移距曲线的梯度（也称为视慢度）；ω为平均慢度不适配与视慢度不适配之间的权重因子；$\bar{d} \stackrel{\text { def }}{\approx} \frac{d}{l}$表示平均慢度数据；$\hat{d} \stackrel{\text { def }}{\approx} \frac{\partial d}{\partial x}$表示视慢度数据。引入视慢度信息约束后，层析反演目标函数是对最小化平均慢度和视慢度的数据不适配，这就意味着层析反演目标函数不仅是对初至走时的反演，而且结合了走时曲线进行反演。借助视慢度信息可反演中、深层速度结构，式（12）中权重因子ω可控制不同慢度信息所占比重，其值越大，中深层反演效果会更佳，反之亦然。当ω = 0时，表示没有视慢度约束。 3.2 正则化信息约束 应用Tikhonov正则化[32]需要使用模型导数算子，并可能产生较为“平滑”的解，这便产生了求解层析反演问题的稳定性与反演结果准确性的折中问题。使用较低阶数的正则化算子对减弱层析反演方程病态性的效果微乎其微，而使用较高阶数的正则化算子会使得层析反演结果变得特别平滑，会造成反演结果损失大量有效信息，因而降低层析结果的可信度。Zhang等[31]经过研究，给出了不同阶数的Tikhonov正则化约束公式，其中一阶导数算子在数值上相当于应用线性插值方法，二阶导数算子则对应于三次样条插值，等价于非线性插值方法的应用。在模型慢度变化为线性时，一阶光滑性准则对于正则化是足够的。然而，实际应用中我们经常会遇到慢度非线性变化的情况，因此，必须应用非线性插值方法。 正如所有反问题一样，初至波走时层析反演的解不是唯一的，存在很多与真实解逼近的模型解，但模型解并非都有实际物理意义。为获得最优解，常见策略是在符合多数物理概念及先验条件时，为有物理意义的数据和模型拟合建立相关定量标准。在视慢度信息约束的基础上，应用Tikhonov正则化来约束模型的平滑度，最小化目标函数Φ(m)可改写为 $ \begin{aligned} \varPhi(m)= & (1-\omega)\left\|C_l[d-G(m)]\right\|^2+ \\ & \omega\left\|D_x[d-G(m)]\right\|^2+\tau\|R m\|^2= \\ & (1-\omega)\|\bar{d}-\bar{G}(m)\|^2+\omega\|\hat{d}-\hat{G}(m)\|^2+ \\ & \tau\|R m\|^2=(1-\omega) S_1+\omega S_2+\tau S_3 \end{aligned} $ (13) 式中：τ为平滑权重因子；Cl为对应射线长度l对旅行时进行缩放的算子，R为正则化算子（例如导数算子、平滑算子）。 采用Scales[33]提出的将高斯-牛顿（GN）法最小化目标函数相关的稳态方程，应用共轭梯度（CG）算法完成迭代更新，得到 $ \begin{gathered} {\left[(1-\omega) A_k^{\mathrm{T}} A_k+\omega B_k^{\mathrm{T}} B_k+\tau R^{\mathrm{T}} R+\varepsilon_k I\right] \Delta m_k=} \\ (1-\omega) A_k^{\mathrm{T}}[\bar{d}-\bar{G}(m)]+\omega B_k^{\mathrm{T}}\left[\hat{d}-\hat{G}\left(m_k\right)\right]- \\ \tau R^{\mathrm{T}} R m_k \end{gathered} $ (14) $ A \stackrel{\text { def }}{=} \frac{\partial \bar{G}}{\partial m}=\frac{1}{l} \frac{\partial G}{\partial m} $ (15) $ B \stackrel{\text { def }}{=} \frac{\partial \hat{G}}{\partial m}=\frac{\partial^2 G}{\partial m \partial x} $ (16) $ m_{k+1}=m_k+\Delta m_k \quad(k=1, 2, 3, \cdots, N) $ (17) 式中：Ak为平均慢度核函数；Bk为视慢度核函数；εkI为变量阻尼项。 采用迭代反演方法线性化一个非线性逆问题在数值计算领域广泛应用。在Morgan[34]提出的GN方法中加入一个可变阻尼参数来解释非线性，可将模型平滑度显式表达为最小化层析反演目标函数。一般由于初始反演模型较差，只有在反演初期才会出现较大的非线性，当目标函数中S1和S2的数据不适配较小，S3（模型的平滑度）占比较大时，层析反演目标函数趋于二次型。 4 实际应用示例 4.1 塔里木盆地甫沙4线束应用示例 应用示例1位于塔里木西南区块（图 5）。为落实地质任务、评价甫沙4油气藏，同时也为塔西南后续地震勘探寻求新思路、新方法奠定基础，在过甫沙4号构造圈部署了甫沙4线束，其目的层为古近系与白垩系。激发参数为：960次覆盖；处理面元大小为7.5 m×15.0 m；选取的最大炮检距为3 000 m；激发线距、激发点距与接收线距均为30 m；接收点距为15 m。结合老资料和野外施工的表层调查资料分析可知，工区表层结构与地形地貌关系密切，工区低速层厚度变化较大，大多为2~80 m。工区南部黄土山体区低降速带速度为300~1 600 m/s；工区西北部和东北部主要为三层结构，高速层速度一般大于1 800 m/s，降速带速度变化较大，为1 300~ 1 600 m/s，低速带速度一般为600~1 300 m/s，高速层顶随地表起伏变化。 下载原图 图 5 塔里木盆地甫沙4线束观测系统 Fig. 5 Observation system of Fusha 4 wiring bundle in Tarim Basin 初至旅行时层析反演中偏移距选取过大时会引起射线路径随反演深度的增大而存在阴影区，造成层析结果失真，降低结果的可信度；另一方面偏移距选取过大，地震数据包含初至量增加，初至错误拾取率相对增加，会从源头影响层析反演结果。因此，处理时仅保留偏移距在3 000 m以内的地震数据进行初至拾取。采用机器学习方法拾取近偏移距地震数据（偏移距小于3 000 m）所有道的初至起跳点，先进行质控，再去除掉初至拾取异常值（异常大或0等），本次拾取结果的精度可满足工业生产需求（图 6）。 下载原图 图 6 塔里木盆地甫沙4线束机器学习方法拾取地震初至起跳 Fig. 6 First break picking using machine learning of Fusha 4 wiring bundle in Tarim Basin 因工区地震数据量巨大，为验证算法的有效性，仅截取工区西部共100炮地震记录进行初至层析反演算法测试，分别采用常规线性梯度方法（图 7a）与微测井信息约束的方法建立反演初始速度场（图 7b）。采用多模板快速推进算法（MSFM）正演计算初始速度场的全局网格走时，再采用起伏地表条件下的波前扩展策略计算全局走时，采用RuggeKutta算法进行射线追踪构建层析核函数，完成反演方程组的构建。 下载原图 图 7 塔里木盆地区甫沙4线束层析反演初始速度场 Fig. 7 Initial velocity used in tomography inversion of Fusha 4 wiring bundle in Tarim Basin 常规初至走时层析反演算法与优化后的反演算法所采用的网格相同，大小为15 m×15 m×10 m。该区块的常规初至层析反演结果存在严重的“牛眼”现象，包含较多高速异常体，低降速带及高速层顶界面均不清晰（图 8a）。优化后初至波层析反演结果基本去除“牛眼”现象，对起伏地表条件下的低降速带刻画效果较好，低速带速度约为800 m/s，降速带速度约为1 300 m/s，且高速层顶界面清晰，速度约为2 000 m/s，均与表层调查结果吻合（图 8b）。 下载原图 图 8 塔里木盆地甫沙4线束优化前、后层析反演速度场对比 Fig. 8 Comparison of tomography inversion velocity field before and after optimization of Fusha 4 wiring bundle in Tarim Basin 常规层析反演结果射线分布杂乱且不合理（图 9a），优化反演结果的射线路径分布更为合理且更均匀（图 9b）。从射线密度对比可看出优化后射线在深度500 m附近明显集中，且在深度1 000 m处存在明显的折射界面，深度大于1 000 m时基本没有射线分布，这与只采用偏移距小于3 000 m的地震数据有关。 下载原图 图 9 塔里木盆地甫沙4线束优化前（a）、后（b）射线密度对比 Fig. 9 Comparison of ray distribution before(a)and after optimization(b)of Fusha 4 wiring bundle in Tarim Basin 4.2 塔里木盆地英买区块应用示例 示例2为塔里木盆地英买区块复杂三维实际地震数据。优化前、后初至走时层析反演结果对比可以看出，常规反演速度场在纵向切片中存在速度异常，并且由于射线分布不均匀且深度较小，造成边界成层阴影现象较重；多信息约束优化后层析反演速度场中，由于正则化的因素，速度更为平滑，消除了“牛眼”异常速度块，边界阴影现象得到明显改善（图 10）。从优化前、后的射线密度场可以看出，优化后反演的射线分布较常规算法的反演射线分布更加均匀、相同位置处射线密度更大，降低了层析反演方程的稀疏性，因而速度反演结果也更为可信（图 11）。 下载原图 图 10 塔里木盆地英买区块优化前（a）、后（b）层析反演速度场对比 Fig. 10 Comparison of tomography inversion velocity before(a)and after optimization(b)inYingmai block Tarim Basin 下载原图 图 11 塔里木盆地英买区块优化前（a）、后（b）射线密度对比 Fig. 11 Comparison of ray distribution before(a)and after(b)optimization in Yingmai block, Tarim Basin 为了量化对比优化前、后反演的地层速度，分别在2种层析反演结果的相同位置抽取速度曲线并叠合绘制，其中绿色曲线为测井速度曲线（井点未在工区内，该井为距工区最近的井位），橙色为多信息约束优化后反演的速度曲线，褐色为常规初至波层析反演的速度曲线（图 12）。从图上可以看出优化前、后的反演速度在深度为200~800 m处差异较小，深度大于800 m时优化后速度曲线与测井速度曲线吻合度更高，即可说明优化后的反演结果更准确，更接近实际的地层速度。 下载原图 图 12 塔里木盆地英买区块优化前、后速度与测井速度曲线对比 Fig. 12 Comparison of velocity curves before and after optimization in Yingmai block, Tarim Basin 5 结论 （1）使用机器学习方法拾取初至可在保证初至拾取精度大于99% 的同时，效率提升80余倍，但使用模拟数据训练的网络用于实际资料时，迁移应用需多次调试才可满足精度需求。 （2）使用视慢度信息与三阶Tikhonov正则化对目标方程进行约束，可增强层析反演方程组求解的稳定性，约束模型解空间，提高了反演结果的准确性。 （3）实际应用显示采用微测井约束的层析反演速度场中浅表层速度较常规反演结果更准确，与测井数据的吻合度更高，可用于后续处理中的深度域速度建模。