基于常分数阶拉普拉斯算子的...与衰减参数的全波形反演方法_胡博涛.pdf

第 卷 第期 年月地球物理学报 ，胡博涛，黄超，董良国等 基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法地球物理学报，（）：，：，.（），（）：，：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法胡博涛，黄超，董良国，张建明同济大学海洋地质国家重点实验室，上海 摘要地震波在地下介质传播过程中由于非弹性衰减的存在将导致能量损失和相位变化，精确的速度与衰减参数建模对油气识别、提高强衰减介质中地震波成像的质量都起着至关重要的作用常分数阶拉普拉斯算子黏声方程由于完全分离的速度频散项与振幅衰减项的优势，以及在强非均质衰减介质中可以高精度求解的特点，已被应用于速度与衰减参数的建模中本文将二阶常分数阶拉普拉斯算子黏声方程拆分为等价的一阶方程组，并在此一阶方程组的基础上推导出新的梯度公式与伴随方程，建立了一种新的速度与衰减参数同时重建的全波形反演方法相较于原二阶常分数阶拉普拉斯算子黏声方程建立的全波形反演流程，数值实验表明，新建立的反演流程可以有效避免原梯度数值计算中的噪声，尤其是可以有效提高衰减参数梯度的反演精度，从而显著提高反演的收敛速度与反演精度关键词品质因子；分数阶拉普拉斯算子；全波形反演；黏声方程 ：中图分类号 收稿日期 ，收修定稿基金项目国家自然科学基金（，），同济大学海洋地质国家重点实验室自主课题（）以及中央高校基础研究基金（）资助第一作者简介胡博涛，男，年生，硕士研究生，主要从事地震波衰减与反演等方面的研究 ：通讯作者黄超，男，研究员，主要从事全波形反演、微地震监测与二氧化碳封存监测方面的研究 ：，地 球 物 理 学 报（）卷 ；引言地下介质通常是黏弹性的，这种非弹性的衰减将引起地震波传播过程中振幅的衰减与相位的延迟（，）地震勘探中通常采用品质因子来表征地震波在地下介质中的这种衰减效应 建立高精度的速度与模型将对准确补偿地震波振幅损失和相位延迟起到关键作用，可以显著改善带补偿的偏移或者最小二乘偏移的成像精度（，；，；刘伟和张剑锋，；，）此外，品质因子还可用于岩性识别和储层预测，对油气检测也具有重要意义 常规的衰减参数建模方法利用地震波衰减随频率的提高而增强的特点，通过地震数据的频移或者能量谱比等方法来建立粗略的模型（邹鹏和程玖兵，）但此类方法精度不高，且严重依赖速度建模的精度，速度建模的误差将严重影响模型的反演精度（，）基于数据匹配的全波形反演方法（）由于充分挖掘了地震信号中的波形信息，理论上可以获得目前地球物理反演方法中最高分辨率的反演结果（，），因而成为目前高精度地质参数建模的重要手段，并在近十几年得到了迅速的发展（，）利用全波形反演技术重建地下模型的建模方法，前人也有了不少的研究（，；，；，；，；，；，；，）波形反演参数建模方法借助数值求解波动方程来预测地震数据，因此利用全波形反演技术实现衰减参数的建模，必然也需要数值求解能够描述地震波衰减机制的黏声或者黏弹波动方程求解黏声或黏弹方程可以在频率域或者时间域进行，在频率域求解用复速度（，）表达的黏声或黏弹方程，因为需要进行大型矩阵的 分解，内存消耗巨大，特别在三维情况下，目前的机器性能往往难以满足，因而限制了该类方法在实际中的应用在时间域求解，则涉及到应力与应变之间的卷积运算，这要求在求解当前时刻的应力和应变时存取以往所有时刻的应力和应变量，对内存的需求同样巨大 为了提高求解的精度并降低内存的消耗，不同的学者利用不同的近似模型得到了不同的近似方程，主要可分为两大类：一类是通过不同物理模型等效近似地球介质的黏弹效应，具有代表性的有 模型（，）、模型（，）以及标准线性体模型（）（，）；另一类则是通过数学模型来近似等效介质在地震波频带范围内所呈现的常特征（，）在这两类模型中，第一类模型对应的黏弹或黏声方程在数值求解过程中依然增加了方程和表征参数的个数，且表征参数与品质因子往往没有显式的数学表达关系，增加了参数反演的难度（，）第二类模型中比较有代表性的方程是 （）提出的分数阶方程以及由此衍生出的一系列修正方程（，；，；，）但在这些方程中，普遍存在空间导数的阶数随衰减参数空间变化的问题，使得在求解该类方程时要求介质的衰减参数在空间平缓变化，因而无法准确模拟强非均质性介质情况下的地震波传播 非均质介质条件下求解空变阶数的拉普拉斯算子可以采用加权近似（，）或者低秩分解方法（，；，），但需要进行额外的矩阵低秩分解 为解决此问题，等（）首次尝试利用常分数阶拉普拉斯算子近似空变阶数的拉普拉斯算子，数值结果表明该方程可以较好地模拟 介质中传播距离小于 时的地震波传播，但对于 的介质其模拟精度有所降低 随后，等（）基于该常分数阶拉普拉斯算子黏声方程讨论了速度双参数全波形反演方法 和 （）则利用级数展开以及系数优化技术提出了一种基于常分数阶拉普拉斯算子的黏声方程该方程中分数阶拉普拉斯算子导数为常数，可以在强非均质性介质中方便求解 此外，上述常分数阶拉普拉斯算子黏声方程均保留了主控振幅衰减和速度频散项分离的特点，为带补偿的偏移成像和波形反演提供了便利为了将常分数阶拉普拉斯黏声方程应用到速度与参数的建模中，和 （）基于他们提出的黏声方程进一步推导了相应的梯度和伴随方期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法程，建立了黏声全波形反演流程 但直接基于该二阶常分数阶黏声方程建立的全波形反演流程在梯度的数值计算中容易出现强烈的数值噪声，需要进行额外的滤波处理（，）处理的效果取决于参数的经验选择，且附加的处理不仅带来计算量的增加，同时也将影响梯度的精度，进而降低反演的精度，特别是对弱参数建模的影响尤为明显（，）此外，由于参数与速度参数的耦合特性，参数建模精度的降低最终也将反过来影响速度参数的反演精度，从而影响后续偏移成像的质量 针对这个问题，本文通过数值测试，分析了梯度计算中数值噪声的产生机制，基于降阶的思想，提出一套新的黏声方程全波形反演流程 利用新的反演流程进行全波形反演，可以有效避免反演过程中的噪声问题，实现了在不额外增加计算量的前提下，提高速度与衰减参数同时建模的精度与反演的收敛速率本文将从 和 （）的常分数阶拉普拉斯算子黏声方程出发，首先分析基于二阶黏声方程导出的梯度公式在计算中产生数值噪声的原因然后针对该噪声的形成机制，基于降阶的思路，提出基于降阶后的一阶方程组推导对应的双参数反演的梯度与伴随方程，建立新的速度与衰减参数反演流程 最后，通过数值实验验证方法的正确性与有效性方法原理 考虑介质黏滞效应的全波形反演方法考虑介质黏滞效应的全波形反演方法，要求所采用的黏声控制方程能够准确描述地震波在传播过程中的吸收衰减效应 和 （）提出的常分数阶拉普拉斯算子黏声方程，由于其与空间参数无关的分数阶拉普拉斯算子和分离的速度频散项与振幅衰减项，具有数值求解精度高、振幅补偿方便的特点，可以方便地应用到全波形反演中 该常分数阶拉普拉斯算子黏声方程可记为：（，）（，）（）（）（，）（）（，），（）其中，为地震波传播速度，（，）为压力场，为参考 角 频 率，为 与 品 质 因 子相 关 的 参 数 （），.（）.鉴于该方程的优势，目前已有了基于该方程建立全波形反演流程来重建速度和衰减参数的研究（，；，）基于方程（）建立全波形反演流程，首先建立如下的基于数据匹配的目标函数：（，）（）（），（）其中，为检波点坐标，为观测数据，为理论合成数据 然后基于黏声方程（）推导出介质参数对应的梯度公式（）和（）来求解方程（）这样一个最优化问题：（，）（）（，）（）（，）（）（，）（，），（）（）（，）（）（，）（）（，）（，），（）其中，（，）为伴随波场，该伴随波场的计算利用了与方程（）相同形式的伴随方程梯度（）和（）基于方程（）推导得到，本文记为基于二阶方程的梯度 但直接通过方程（）和（）数值计算的梯度往往存在较强的噪声干扰（，），影响了反演的效果 下面采用一个简单数值实验来分析此类噪声的产生机制本实验使用均匀介质模型，模型尺寸为 ，离散网格的大小为 实验采用单炮检对的观测方式，炮点位于（，）（，），检波点位于（，）真实模型为均匀介质，其中地震波传播速度为 ，品质地 球 物 理 学 报（）卷因子为 初始模型采用速度为 、品质因子为 的均匀模型 图 和 为采用公式（）和（）计算的第一代梯度，从图 中可以看到明显的数值噪声 为进一步分析此噪声的产生机制，将梯度公式（）和（）分别分解为三部分：与（）有关的项、与（）有关的项和剩余项（如表所示）图展示了梯度 和 分解后对应的三个组成部分 从图可以清楚地看出，数值噪声主要来自于与（）有关的部分（图 和），其产生原因在于高阶导数的数值求解带来的较大数值误差 而图 中的速度梯度之所以看不到明显的噪声污染，主要在于图 的数值量级远高于图，在三者相加后压制了图 中的数值噪声 但是，衰减参数梯度中含（）项（图）与含（）项（图）在数值上为同一数量级，二者比不含分数阶项（图）的数值高出两个数量级，含（）项中的较高的数值噪声严重污染了最终的梯度（图）可见，图 梯度中的数值噪声主要来源于含（）项，减少噪声的关键在于如何避免对高阶导数的直接数值求解图（）利用公式（）计算的速度的梯度；（）利用公式（）计算的参数的梯度 （）（）；（）（）表二阶梯度各组分表达式 二阶梯度成分分解速度衰减参数含（）项（）（，；）（）（，；）（）（，；）（）（，；）不含分数阶的项（，；）（，；）含（）的项（）（，；）（）（，；）基于一阶黏声方程组的全波形反演方法为避免在梯度公式中出现高分数阶导数，可以考虑将方程（）拆分为三个一阶方程构成的方程组：（，）（，），（）（，）（，），（）（，）（）（）（，）（）（，）（，）.（）其中，（，），（），（，）为质点速度，（，）为声学密度，（，）为震源，为参考密度 将（）式代入（）式可得：（，）（）（）（，）（）（，）（，）.（）拆分的方程组（）、（）、（）与方程（）等价，但相比于方程（），拆分后的方程组利用变量替换的方式将原方程中的高分数阶导数项（如（）降阶为低分数阶导数项（如（）和（），避免了对波场直接求取（）的高分数阶导数，期望以此来实现减少梯度中数值噪声的目标下面，将基于此一阶方程组导出相应的全波形反期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法图基于二阶方程的梯度成分分解：梯度中含（）项的成分（）速度，（）参数；梯度中不含分数阶项的成分（）速度，（）参数；梯度中含（）项的成分（）速度，（）参数 ：（）（）（）；（）（）；（）（）（）演方法基于常分数阶拉普拉斯算子的一阶黏声方程组（）、（）、（），利用拉格朗日乘子法来推导新的速度和衰减参数对应的梯度公式和相应的伴随方程 首先将方程（）、（）和（）作为约束条件引入目标函数（），构造新的目标函数：（，）（，）（）（）（）（）（），（）其中，、为伴随波场 由于各变 量 相 互 独立，当 目 标 函 数取 得 最 小 值 时，满 足 以 下关系：，（），（）地 球 物 理 学 报（）卷 ，（），（）.（）将（）式代入（）（）式，化简后可得目标函数对速度和与品质因子相关的参数的梯度：（）（）（），（）（）（）（）.（）将（）、（）式转换到频率域，即可得到频率域反演中的梯度公式：?（）（）?（）?，（）?（）（）?（）?，（）其中?和?以及?分别代表和以及变换到频率域的波场，?代表?在频率域的共轭波场 方程（）和（）为基于一阶方程组（）、（）、（）导出的速度和衰减参数对应的梯度，本文记为基于一阶方程的梯度 观察公式（）和（）可以看出，该梯度公式中分数阶拉普拉斯算子也均为低阶，有效避免了梯度公式中高阶导数项的出现而利用公式（）、（）计算梯度，首先需要计算公式中的伴随波场（，），而伴随波场的计算则涉及到求解对应的伴随方程 接下来，将推导伴随波场（，）对应的伴随方程 将（）式代入（）（）式，利用分部积分化简可得到新的伴随方程（详细的推导见附录）：（，）（，），（）（，）（，）（，），（）（，）（）（）（，）（）（，）.（）其中，为伴随震源 从形式上看，新的伴随方程（）、（）、（）与原方程组（）、（）、（）还存在一些差别，这表明一阶黏声方程组（）、（）、（）并非自伴随方程，计算伴随波场时必须采用新导出的伴随方程（）、（）、（）来求解 若还采用原方程计算伴随波场，将得不到正确的伴随波场，必然影响后续的参数更新 至此，我们已经基于一阶黏声方程组，重新建立了一套新的速度与衰减参数的全波形反演流程接下来，以图中的参数设置为例，采用同样的观测数据和初始模型，利用新推导得到的梯度公式（）和（）计算的第一代梯度如图所示 对比图和图中展示的梯度可以看到，基于一阶方程组导出的新梯度在数值计算中未产生强烈的数值噪声，极大改善了衰减参数对应梯度的精度，实现了我们预期的目标图（）利用公式（）计算的速度的梯度；（）利用公式（）计算的参数的梯度 （）（）；（）（）期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法数值实验在上述方法原理部分，我们基于降阶后的一阶方程组表示的黏声方程，推导出了相应的梯度和伴随方程，建立了新的速度与衰减参数全波形反演流程 为验证该反演流程在重建高精度速度与衰减参数的优势，下面将通过两个理论模型实验来证明该方法的正确性与有效性 在数值实验中，均采用交错网格有限差分法与伪谱法相结合的方式在时间域数值求解黏声方程，而反演过程则在频率域进行，且采用 的梯度优化算法对梯度进行预条件 球状异构模型为验证新反演流程在重建速度和衰减参数上的优势，本节使用球状异构侵入体模型来对比基于一阶方程与二阶方程梯度的波形反演方法的优劣介质模型中速度与衰减的异常体直径均为 其中中心位于（，）处的球状异常体在参考频率为 情况下的真实相速度为 ，背景速度为 中心位于（，）球状异常体的真实品质因子为，背景为 图所示的介质模型大小为 ，离散网格的大小为 实验采用井间观测系统来采集地震数据，在左侧的炮井中激发 炮，震源深度从 到 均匀分布，炮间距 每炮采用 个检波器接收，均匀分布在介质右侧的观测井中，深度从 到 均匀分布，检波器间距为 震源选用主频为 的 子波 反演中选择频率从 到 共 个频率成分的地震数据进行反演，初始模型均采用均匀背景介质模型图展示了反演中利用基于一阶和二阶方程的梯度计算的第一代梯度 比较速度的梯度，基于一阶与二阶方程的梯度接近，但图 中对应速度异常区域的梯度的边界分辨率稍高 而对于衰减参数的梯度，采用基于二阶方程的梯度公式（）计算的梯度在炮点附近出现了强烈的数值噪声（如图），从而影响衰减参数的准确更新 而利用公式（）计算的梯度则未受到数值噪声的影响，可以更好地实现对衰减参数的反演更新 图和图分别展示了两种方法的最终反演结果 从图中可以看出，与梯度的分析一致，两种方法对速度异常体的重建结果（图、）大体近似 而基于二阶方程的反演，由于梯度中数值噪声的存在，降低了品质因子的反演精度但采用基于一阶方程的梯度方法进行反演，品质因子的反演结果精度得到了较大改善（图）图所示的反演的速度和衰减参数的抽线结果，也进一步验证了我们的结论 图展示了两种反演方法的目标函数收敛曲线 从收敛曲线可以明显看到，采用基于一阶方程的梯度的反演方法相较于二阶方程的梯度的反演方法收敛速度更快，目标函数下降程度更低，从另一个方面证明了基于二阶方程的梯度中的数值噪声对反演收敛速度和反演精度的影响图真实模型（）速度模型；（）品质因子的倒数模型 （）（）地 球 物 理 学 报（）卷图迭代反演过程中第一代梯度对比对应于（）二阶方程与（）一阶方程的速度的梯度；对应于（）二阶方程与（）一阶方程的参数的梯度 （）（），（）（）气烟囱模型上述异构侵入体的数值例子验证了基于一阶方程的梯度反演方法较基于二阶方程的梯度反演方法的优势 下面将采用气烟囱模型（，）来验证本文提出的基于一阶方程的梯度反演方法在复杂模型反演中的有效性模型中海底火山喷出的气体在近海 左右形成了一个低速强衰减的气云区，这个气云区的存在，将严重影响气云下方的结构成像，为此，需要准确的速度与模型来帮助提高气云区下方的成像质量和精度 测试中采用的模型（图）尺寸为 ，离散网格的大小为 实验中共激发 炮，震源位置从 到 ，均匀分布在海水表面，炮间距 每炮采用均匀布设在地表的 个检波器接收，检波器间距为 震源选用主频为 的 子波 反演中选取频率从 到 的地震数据，频率间隔为 ，共 个频率成分同时反演 速度初始模型为高斯平滑后得到的光滑模型（图 ），品质因子初始模型是值为 的均匀模型（图 ）图 展示的是采用基于一阶方程的梯度方法重建的速度与品质因子模型从图 中可以看出，通过基于一阶方程的梯度方法的反演重建，速度参数获得了高精度的恢复，特期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法图基于二阶方程的梯度方法的反演结果（）速度；（）品质因子的倒数 （）；（）图基于一阶方程的梯度方法的反演结果（）速度；（）品质因子的倒数 （）；（）别是浅层的一些细微结构都得到了较好地重建（如图 ）而从品质因子的反演结果看（图 ），气云区对应的强衰减边界都被很好地刻画出来了，且值的大小也与真实值接近 这样的反演结果表明，本文提出的基于一阶方程的梯度的 反演方法可以对复杂模型的速度和品质因子的精细结构实现高精度的同时重建 图 展示了某单炮观测数据分别与利用初始模型合成的地震数据和利用最终反演模型合成的地震数据的残差从数据残差的减小程度，也进一步说明利用反演结果模拟的合成数据与观测数据中的主要波现象都得到了充分的匹配 而图 目标函数收敛曲线（反演频带为 ，频率间隔为 ，共 个频率成分）也从数值上进一步验证了数据的匹配情况 鉴于实际地震数据的有效频带一般都在 以上，为了贴近实际情况，我们还进行了利用 地震数据重建速度地 球 物 理 学 报（）卷与衰减参数的反演测试，图 展示了最终的反演结果 与图 的结果比较发现，两者比较接近，说明图反演结果在 处的抽线对比（）速度；（）品质因子的倒数 （）；（）图目标函数下降曲线 在当前的初始模型（图）情况下，本方法不管从 开始还是 开始反演，均未发生频率跳周现象，可以较好地重建地下介质的速度和衰减参数模型本实验中重建高精度的速度与品质因子参数的目标，是为了实现对气云区下方区域的结构进行高质量的成像（如图 虚线框所示）图 展示了利用初始模型（图）对数据进行时间域带补偿的逆时偏移成像结果 从图上可以看出，由于速度与衰减参数误差的存在，气云区下方的成像同相轴连续性差，且能量微弱 而通过基于一阶方程的梯度的反演方法反演以后，获得了高精度的速度与衰减参数反演结果，利用该结果重新进行逆时偏移成像，其结果如图 所示，同相轴的连续性清晰可见，且气云区下方的成像能量得到了很好地补偿，进一步验证了 反演结果的可靠性实际数据往往容易受到外界噪声的干扰，为了验证方法对噪声的敏感程度，下面将进行含噪数据的反演测试 本实验仍然利用气烟囱模型实验中的参数（反演频带也为 到 ，频率间隔为 ），只是按照不同的信噪比在观测数据中加入噪声，加入噪声的单炮记录如图 所示 利用加入噪声后的地震数据进行反演，得到图 的反演结果 从图中可以看到，随着数据信噪比的降低，反演结果受到的影响越来越强 但即使数据的信噪比低到 时，反演依然对初始模型进行了一定程度的更新，表明了该反演方法具有一定的抗噪能力，但随着信噪比的降低，反演结果的精度也呈现降低的趋势 从图 的目标函数下降曲线也可以看出，随着噪声在数据中占比的提高，反演中目标函数的下降比例越来越低，自然对反演的影响越来越严重，这与图 的反演结果一致图 真实模型（）速度；（）品质因子的倒数 （）（）期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法图 初始模型（）速度；（）品质因子的倒数 （）（）图 基于一阶方程的梯度方法的反演结果（）速度；（）品质因子的倒数 （）；（）图 单炮记录数据残差比较（）观测炮记录；（）观测记录与利用初始模型合成记录的残差；（）观测记录与利用反演结果合成记录的残差 （）；（）；（）图 目标函数下降曲线 结论本文基于降阶的一阶黏声方程组，推导得到了全新的梯度公式以及与其对应的伴随方程，建立了同时重建速度与衰减参数的全新的黏声全波形反演方法 该方法在梯度计算中可以有效避免由于求解高阶导数项带来的数值噪声对梯度的干扰，对衰减参数梯度的改善尤其明显提出的这种基于一阶方程梯度的全波形反演方法，可以实现对速度和品质地 球 物 理 学 报（）卷图 数据频带为 情况下的一阶方程的梯度方法的反演结果（）速度；（）品质因子的倒数 ：（）；（）图 带衰减补偿的逆时偏移成像结果对比（）基于初始模型（图）；（）基于反演结果（图）（）（）（）图 第 炮地震记录（直达波已切除）（）无噪声；（）添加信噪比为 的高斯噪声；（）添加信噪比为 的高斯噪声 （）；（）（）；（）因子双参数的高精度同时重建，可为后续成像精度的提高提供准确的速度模型和可靠的模型附录 基于分数阶拉普拉斯算子的黏声方程伴随方程的推导从正文中的方程（）出发，当方程（）取得最小值时，其中的，和应满足以下三个方程：，（），（）.（）将方程（）分别代入方程（）、（）和（），可以得到：期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法图 添加高斯噪声生成不同信噪比地震数据情况下的反演结果（）信噪比为 的速度反演结果；（）信噪比为 的品质因子的倒数反演结果；（）信噪比为 的速度反演结果；（）信噪比为 的品质因子的倒数反演结果 （）；（）；（）；（）（，）（，）（，）（，）（，）（，），（）（，）（，）（，）（，）（，）（）（）（，）（，）（）（，）（，），（）（，）（，）（，）（，）（，）（）（）（，）（，）（）（，）（，）.（）利用分部积分可得 （，）（，），（）地 球 物 理 学 报（）卷 （，）（，）（，），（）（，）（）（）（，）（）（，）.（）化简可得对应一阶方程组（）、（）、（）的伴随方程组：（，）（，），（）（，）（，）（，），（）（，）（）（）（，）（）（，）.（）图 不同信噪比地震数据进行反演的目标函数下降曲线 ，：，（）：，：，：，（）：，：，（）：，：，（）：，：，（）：，（）：，（）.，（）：，：，：，（）：，：，：，（）：，：，（）：，：（），（）：，：，：期胡博涛等：基于常分数阶拉普拉斯算子的黏声方程重建速度与衰减参数的全波形反演方法 ，（）：，：，（）：，（）：，：，（）：，：，（）：，：，（）：，（）：，：，（）：，：，（）：，：，（），（）：，：附中文参考文献刘伟，张剑锋 黏弹性叠前时间偏移：陡倾角构造成像与实际应用地球物理学报，（）：，：邹鹏，程玖兵 黏声方程值反射波反演地球物理学报，（）：，：（本文编辑何燕）