截断时窗的低通滤波多尺度全波形反演方法

截断时窗的低通滤波多尺度全波形反演方法
【专利摘要】本发明涉及一种截断时窗的低通滤波多尺度全波形反演方法。该方法首先在低频段内，用截断时窗反演策略进行全波形反演，得到的反演结果作为初始模型用于中频段反演。然后将中频段的反演结果作为初始模型，用于常规时间域全波形反演并得到最终反演结果。截断时窗反演策略节约了大量正演时间。低通滤波多尺度反演策略，可将地震数据分为不同频段进行反演。将截断时窗与低通滤波组合使用，可在缺失低频和强高斯噪声干扰情况下，得到较好的反演结果。模型测试结果表明，该组合策略在提高计算效率的同时缓解了周跳现象，且截断时窗的低通滤波多尺度全波形反演方法在反演精度、计算效率、低频依赖、抗噪能力等方面优于常规时间域全波形反演。
【专利说明】
截断时窗的低通滤波多尺度全波形反演方法
技术领域
[0001] 本发明涉及一种地震勘探的数据处理方法，尤其是截断时窗和低通滤波多尺度反 演策略。整个计算过程利用超随机震源编码策略进行全波形反演，目的是加快整个全波形 反演的进程。
【背景技术】：
[0002] 随着石油工业的发展和勘探开发的不断深入，从构造勘探阶段逐渐走向岩性勘探 阶段，其勘探难度逐渐加大。为了顺应这种要求，地震数据处理技术不断完善，在这期间全 波形反演迅速发展起来，并成为当今地球物理界的研究热点。20世纪80年代，Tarantola提 出了时间域全波形反演理方法，并将目标函数对模型参数求导通过残差反传波场和正传波 场互相关运算得到，从而避免了求取Jacobi矩阵。20世纪90年代，Pratt将全波形反演推广 到了频率域，提出只需要几个离散的频率就可以得到高精度的反演结果，而且低频到高频 的反演策略可以解决陷入局部极小值的问题。在频率域全波形反演的基础上Shin等人发展 了 Laplace-Fourier域全波形反演，该方法能够在缺失低频的情况下可以获得一个高精度 初始模型。Wang等人(2011)提出在GPU上基于CUDA编程的加速全波形反演，该方法有效地加 快了全波形反演的计算效率。
[0003] Bunks(1995)提出基于低通滤波的多尺度全波形反演，该方法指出先对观测数据 进行低通滤波，再对估计的震源子波用以相同滤波器处理，可在一定程度上解决了全波形 反演的周跳问题。但是利用该多尺度方法进行全波形反演依然存在着一些问题，如:地层是 一个粘弹性介质，而且波动方程正演是一个复杂的过程，当地震波在地下传播的过程中会 出现频移效应或者其他复杂变化。利用低通滤波之后的震源子波进行正演，会出现正演记 录和观测记录不匹配的问题，关于这个问题Bunks等人并没有给出很好的解决方案。Chen和 Wu等人利用时间衰减滤波器同时处理观测数据和正演数据，即在地震数据上乘上一个随着 时间增加而衰减的指数函数，这样先演浅层信息，忽略深层构造。随着反演进行逐渐减小衰 减函数的影响，使得反演由浅层到深层逐层反演。该方法极大的减小了每一步的反演参数， 一定程度上解决了全波形反演的的周跳问题。但是该方法在正演过程中浪费了大量的时 间，因为正演得到的正演记录再乘以一个衰减函数，深层信息基本衰减至零。该方法没有利 用到深层信息，却浪费了大量时间来正演深部波场，很大程度上降低全波形反演的计算效 率。

【发明内容】
：
[0004] 本发明的目的是针对上述现有低通滤波之后数据不匹配问题，提供一种截断时窗 的低通滤波多尺度全波形反演方法。
[0005] 本发明的解决方案是:利用巴特沃斯最平滑低通滤波器同时处理观测记录和正演 记录，这样有效避免了因数据不匹配而产生错误的反演结果，该方法更适合实际地震数据 处理流程，可以得到精确的反演结果。针对时间域全波形反演计算效率低的问题，本发明提 出截断时窗全波形反演策略，即在刚开始反演的时候，只截取实际数据前部分时间进行反 演，所以也只需要正演对应的时窗数据进行数据匹配，这样通过正演截断时窗内的波场来 代替正演全部时间波场，极大地缩短了正演所需要的时间。随着反演的进行，逐渐增加反演 时窗长度，该截断时窗反演策略能够减少每一步反演参数的个数，同时增加了程序的稳定 性和反演结果的可靠性。
[0006] 本发明的目的是通过以下技术方案实现的：
[0007] 截断时窗的低通滤波多尺度全波形反演方法的核心是在MATLAB 2013a平台上完 成的。利用巴特沃斯滤波器对观测数据和正演数据进行低通滤波，通过设定巴特沃斯最平 滑滤波器不同的截断频率，来实现时间域多尺度反演策略。低通滤波器能够有效滤除地震 数据中的高频成分，保留地震数据中的低频成分，由于低频成分含有地下构造宏观信息，有 利于初始模型的建立，同时避免周跳现象的发生。截断时窗反演策略，即选择不同长度的时 间窗口来提取观测数据用于全波形反演，同时按照该时窗对应的时间长度进行正演，这样 可以避免每一步都正演所有时间数据，极大程度上节约了计算时间。
[0008] 截断时窗反演策略和低通滤波多尺度反演策略本身都具有压制周跳和提高反演 精度的作用。为了得到更稳定的反演能流程、更精确的反演结果和更高效的反演方法，本发 明是进一步把二者结合起来使用(如图1所示）。
[0009] 截断时窗的低通滤波多尺度全波形反演方法，包括以下步骤：
[0010] a、安装MATLAB2013a软件，并安装MATLAB语言编写的地震数据处理软件Crews工具 包；
[0011] b、对采集的地震数据进行预处理，并将观测记录和观测系统输入计算机进行时间 域全波形反演；
[0012] c、从观测记录中估计震源子波，利用估计得到的震源子波进行正演(本发明利用 雷克子波代替实际数据中的震源子波）；
[0013] d、建立线性递增初始模型，该速度模型作为反演的起始值，通过计算模型的更新 方向，实现由初始模型逐渐向真实模型靠近的过程；
[0014] e、根据最小二乘原理构造目标函数：
，其中V表示地下P波速 度，dobs为采集的观测记录，(1。31在速度模型通过正演得到的正演记录。在求目标函数的梯度 过程中对目标函数两端关于P波速度V求导，得到梯度表达式：
[0016] 其中Pb表示由伴随震源fs向地下传播得到的波场，Pf表示正演波场，s表示震源的 数目，r表示检波器的数目。
[0017] f、选取截断时窗的时间长度，截取观测记录前部分地震信号，并记录截断时窗长 度，该截断时窗的时间长度作为正演的时间长度；
[0018] g、按照声波方程的动力学规律，利用输入的观测系统、估计得到的震源子波和初 始速度模型进行地震波场正演。该步骤只需要正演截断时窗内的地震数据，不需要将所有 时间地震数据全部正演出来。存储正传波场数据和正演记录。
[0019] 根据要求设定截断时窗的低通滤波多尺度全波形反演方法相关参数，包括模型大 小nzXnx，网格距dx，dz，正演时窗长度It，时间采样间隔dt，雷克子波主频fo,巴特沃斯滤 波器的截断频率f?t，每个超随机震源的最大迭代次数itermax，更新梯度最小值gtol，目标 函数精度tol，反演结果的最大值vmax与最小值vmin;
[0020] h、对截断时窗内的观测记录和正演记录同时进行低通滤波，确保用相同的低通滤 波器处理数据。滤波之后，得到处理观测记录和处理正演记录；
[0021] i、将处理观测记录和处理正演记录做差，得到残差数据；
[0022] j、将残差数据作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的残 差反传波场；
[0023] k、正传波场与残差反传波场做互相关，得到模型更新梯度；
[0024] 1、用超记忆梯度优化算法计算模型更新方向，并通过Wolfe收敛准则寻找合适的 步长；
[0025] m、判断是够满足终止条件，若满足则输出截断时窗的低通滤波多尺度全波形反演 方法反演结果。若不满足终止条件，将反演结果作为下一个循环的初始模型，返回d步骤；
[0026] n、将上一步的输出结果作为常规时间域全波形反演的初始速度模型，然后进行常 规时间域全波形反演。迭代200次之后输出最终反演结果。
[0027] 有益效果:通过对原有理论方法的改进和相关地震全波形反演技术的整合，本发 明成功地将截断时窗和低通滤波等技术应用到了全波形反演中。
[0028] 提出的截断时窗的低通滤波多尺度全波形反演方法，该方法不仅能够缓解因低频 缺失带来的周跳问题，而且极大程度的节约了正演的时间。利用截断时窗反演策略显著提 高了全波形反演的计算效率和反演精度。利用巴特沃斯滤波器对观测数据和正演数据进行 低通滤波，滤除高频成分，保留低频成分，通过设定巴特沃斯最平滑滤波器不同的截断频 率，来实现频率由低到高的时间域多尺度反演策略，有效地压制了全波形反演的周跳现象。 该反演策略一定程度上提高了全波形反演的精度，且更适合实际地震数据处理，同时本发 明提出的低通滤波反演策略无需对地震子波进行滤波，有效避免因数据不匹配而造成的反 演错误。在整个计算过程中震源都采用超随机震源编码策略，目的是在加快反演速率的同 时压制超级炮内部的串扰噪声。与常规时间域全波形反演相比，本发明反演地下速度的过 程中用了截断时窗和低通滤波反演策略，解决了以下问题：
[0029] 1、首先将截断时窗和低通滤波策略组合起来，应用到全波形反演中，可快速得到 一个高精度的初始速度模型，该初始速度模型基本恢复出了Marmousi模型的宏观信息，利 用该初始速度模型进行常规全波形反演，可以缓解周跳问题。
[0030] 2、将截断时窗函数引入到全波形反演中，使得每次正演不需要计算所有时刻的波 场，只计算截断时窗内的波场即可，极大地缩短了正演所需要的计算时间。
[0031] 3、截断时窗低通滤波多尺度反演策略能有效减少每一次反演的未知数个数，根据 互相关求得的梯度和模型更新方向更加准确。由于在模型更新迭代的过程中，需要通过 Wolfe线性搜索寻找合适的迭代步长，当梯度方向或者下降方向不准确的情况下，Wolfe线 性搜索会浪费大量的时间，而本发明提出的截断时窗低通滤波反演策略能较准确的得到模 型更新的梯度方向，这样在线性搜索的过程中节约大量的计算时间。
[0032] 4、常规全波形反演需要存储正传波场和残差反传波场，存储这两个波场需要大量 的计算机内存。本发明提出的截断时窗反演策略，只需要存储时窗内的正传波场和残差反 传波场即可，一定程度上节约了所需计算机的内存空间，尤其在实际数据处理中。
[0033] 5、超随机震源编码策略能有效压制超级炮内部的串扰噪声，该方法在保证反演精 度的同时一定程度上减少了正演的次数，有效地缩短了全波形反演的计算时间。
[0034] 模型测试结果表明，即使在缺失低频数据和很强的高斯噪声干扰的情况下，利用 截断时窗低通滤波策略的全波形反演也可以得到较好的反演结果。该反演策略既能实现在 低频段由浅层到深层的逐层反演方法，又能在计算效率方面得到显著的改善，同时缓解了 因缺失低频数据而导致的周跳现象。可以看出本发明提出的截断时窗低通滤波反演策略可 为常规全波形反演提供一个高精度的初始速度模型。
【附图说明】
[0035] 图1截断时窗的低通滤波多尺度全波形反演方法流程图。
[0036]图2雷克子波
[0037](左)雷克子波波形图，（右)雷克子波频谱图。
[0038]图3巴特沃斯最平滑滤波器图。
[0039]图4速度模型
[0040] (左)线性递增初始速度模型图，（右)真实速度模型图。
[0041] 图5滤波波形图和频谱图
[0042] (a) 15Hz低通滤波波形图和真实地震记录波形图，
[0043] (b)25Hz低通滤波波形图和真实地震记录波形图，
[0044] (c)15Hz低通滤波，25Hz低通滤波和真实地震记录频谱图，
[0045] (d)lOHz高通滤波波形图，
[0046] (e)lOHz高通滤波频谱图。
[0047]图6超随机震源观测记录 [0048](左)无噪数据，（右)含噪数据。
[0049] 图7截断时窗示意图。
[0050] 图8观测记录低通滤波 [0051](左）15Hz低通滤波观测记录，
[0052](右)25Hz低通滤波观测记录。
[0053]图9常规时间域全波形反演结果图。
[0054]图10低通滤波多尺度全波形反演结果图 [0055] (a) 15Hz低通滤波全波形反演结果图，
[0056] (b)25Hz低通滤波全波形反演结果图，
[0057] (c)低通滤波最终反演结果图。
[0058]图11截断时窗全波形反演结果图。
[0059] 图12截断时窗低通滤波多尺度全波形反演结果图。
[0060] 图13常规时间域全波形反演结果图（缺失低频测试）。
[0061 ]图14低通滤波多尺度全波形反演结果图（缺失低频测试）。
[0062] (a) 15Hz低通滤波全波形反演结果图，
[0063] (b)25Hz低通滤波全波形反演结果图，
[0064] (c)低通滤波最终反演结果图。
[0065] 图15截断时窗全波形反演结果图（缺失低频测试）。
[0066] 图16截断时窗低通滤波多尺度全波形反演结果图（缺失低频测试）。
[0067] (a)截断时窗+15Hz低通滤波全波形反演结果图，
[0068] (b)25Hz低通滤波全波形反演结果图，
[0069] (c)截断时窗低通滤波最终反演结果图。
[0070]图17常规时间域全波形反演结果图（含高斯噪声测试）。
[0071 ]图18截断时窗全波形反演结果图（含高斯噪声测试）。
[0072] 图19低通滤波多尺度全波形反演结果图（含高斯噪声测试）。
[0073] 图20截断时窗低通滤波多尺度全波形反演结果图（含高斯噪声测试）。
【具体实施方式】
[0074]下面结合附图和实例对本发明进一步的详细说明
[0075]截断时窗的低通滤波多尺度全波形反演方法，包括以下步骤：
[0076] a、程序是在MATLAB2013a软件框架下编写完成，首先需要安装MATLAB2013a软件， 并安装MATLAB语言编写的地震数据处理软件Crews工具包；
[0077] b、采集的地震数据首先需进行预处理，并将观测记录和观测系统输入计算机进行 时间域全波形反演。其中预处理包括：
[0078] B1、多次波衰减、消除交混回响和压制虚反射等不能正演的地震波。
[0079] B2、低频保护去噪、缺失地震道补偿、保幅处理等。
[0080] B3、定义观测系统。
[0081] c、从观测记录中估计震源子波，该震源子波用于正演(本发明利用雷克子波代替 实际数据中的震源子波(图2));
[0082] d、建立线性递增初始模型（图4左），该速度模型作为反演的起始值，通过计算模型 的更新方向，来实现由初始模型逐渐向真实模型靠近的过程；
[0083] e、根据最小二乘原理构造目标函数：
，其中v表示地下P波速 度，cUs为采集的观测记录，(1。31在速度模型通过正演得到的正演记录。在求目标函数的梯度 过程中对目标函数两端关于P波速度v求导，得到梯度表达式：
[0085]其中Pb表示由伴随震源fs向地下传播得到的波场，pf表示正演波场，s表示震源的 数目，r表示检波器的数目。
[0086] f、选取截断时窗（图7)，截取观测记录前部分地震信号，并记录截断时窗长度，该 截断时窗的时间长度作为正演的时间长度；
[0087] g、按照声波方程的动力学规律，利用输入的观测系统、估计得到的震源子波和初 始速度模型进行地震波场正演。该步骤只需要正演截断时窗内的地震数据，不需要将所有 时间地震数据全部正演出来。存储正传波场数据和正演记录。
[0088] 根据要求设定截断时窗的低通滤波多尺度全波形反演方法相关参数，包括模型大 小nzXnx，网格距dx，dz，正演时窗长度It，时间采样间隔dt，雷克子波主频fo,巴特沃斯滤 波器的截断频率(f?t)，每个超随机震源的最大迭代次数itermax，更新梯度最小值gtol，目 标函数精度tol，反演结果的最大值vmax和最小值vmin;
[0089] 超随机震源编码:时间域全波形反演在单炮正演的过程中耗费大量的时间，有人 提出利用同时激发多个震源(超级炮)来代替单炮激发震源，这样能缩短正演所需要的时 间，但是利用同时激发震源进行全波形反演又给反演结果带来了串扰噪声。为了提高全波 形反演的计算效率同时压制串扰噪声本软件利用超随机震源编码策略。
[0090] 随机相位编码地表超级炮激发出的雷克子波，要求相位正负随机分配，同时每个 超级炮中的单炮激发延时不同，延时最大不超过400ms。设随机相位编码信息为：
[0092]其中s'US={1，2…Ns}，s中元素的个数固定，其表示一个超级炮中所含单炮数目 个数。超级炮中单炮数量可根据实际需求调整混合比例。s中元素大小随机分配。iGN+表示 随机的正整数。超随机震源编码可以表示为：
[0093] spf = Dy ? A ? r ? T ? L ? F ? f (4)
[0094]其中Dy表示动态震源编码函数(动态编码的意思就是在超记忆梯度法每迭代10次 之后重新更换随机相位编码信息，最终可实现地表全覆盖。该方法能够有效地增加地表炮 的覆盖密度，同时减弱串扰噪声对全波形反演的影响，最终实现计算效率上的巨大飞跃。）， A表示随机振幅编码函数，r表示随机相位编码函数，T表示随机延时编码函数，L表示随机 位置编码函数，F表示随机雷克子波主频编码函数，f震源函数。
[0095] h、对截断时窗内的观测记录和正演记录同时进行低通滤波，确保用相同的低通滤 波器(图3)处理数据。滤波之后，得到处理观测记录和处理正演记录；
[0096]巴特沃斯低通滤波器：以巴特沃斯函数来近似滤波器的系统函数。巴特沃斯滤波 器是根据幅频特性在通频带内具有最平滑特性定义的滤波器。
[0097]巴特沃思滤波器的低通模平方函数表示
[0099]巴特沃斯滤波器的主要特征如下：
[0102]把、|凡(^)|2是如勺单调下降函数；
[0103] H4、|Ha( j Q ) |2随着阶次N的增大而更接近于理想低通滤波器；
[0104] 滤波器的幅频特性随着滤波器阶次N的增加而变得越来越好，在截止频率n。处的 函数值始终为1/2的情况下，通带内有更多的频带区的值接近于1;在阻带内更迅速的趋近 于零。
[0105] i、将处理观测记录和处理正演记录做差，得到残差数据；
[0106] j、将残差数据作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的残 差反传波场；
[0107] k、正传波场与残差反传波场做互相关得到模型更新梯度；
[0108] 全波形反演的目标函数一般定义为正演数据与观测数据残差的L2泛函，形式如 下：
[0110] 其中0。31便是正演记录，Dcal表示实际野外的观测记录。将目标函数对速度v求导得 到：
[0112]本发明为了更清晰的阐述理论方法，假设密度为常数，则只需要对速度进行反演。 常密度声波方程为：
[0114]其中u为声波波场，f为震源函数。v表示地下介质的声速度。
[0115] 将声波方程简化写成算子的形式：
[0116] Lu = s (9)
[0118] 方程(9)两边同时对速度v求导得到(Pratt et al.，1998):
[0120]得到波场对地下速度v的导数为：
[0122]将方程(11)代入方程(7)中得到：
[0124] 其中I/1表示伴随算子或称为反传算子，表示残差反传波场。fs = Dcal-D〇bs称为伴随震源。
[0125] 1、用超记忆梯度优化算法计算模型更新方向，并通过Wolfe收敛准则寻找合适的 步长。
[0126] 超记忆梯度法：用当前梯度及之前多步梯度的信息来对当前梯度的信息进行校 正，以加快目标函数的收敛速度。其迭代形式为：
[0127] vk+i = Vk+akdk (13)
[0128] 其中Vk代表迭代第k部速度模型，ak表示步长，dk表示下降方向。
[0129] 超记忆梯度法下降方向可以通过如下形式来表示：
[0133] 其中，0<P<l，mem>0是一个正整数，mem我们称之为记忆度。不同记忆度目标函 数收敛速度及模型反演精度不同，记忆梯度数量过多则会导致梯度方向校正过度，目标函 数下降变慢，记忆梯度数量少则退化为共辄梯度法。经过实验验证得到，在地震全波形反演 中最好的记忆度大致在mem=4，5。超记忆梯度法有着一定的优势，但当梯度记忆过多，计算 速度变慢，内存需求增大。所以选择合适的记忆度才能发挥超记忆梯度法最大的优势。
[0134] m、判断是够满足终止条件，若满足则输出截断时窗的低通滤波多尺度全波形反演 方法反演结果。若不满足终止条件，将反演结果作为下一个循环的初始模型，返回d步骤。
[0135] n、将上一步的输出结果作为常规时间域全波形反演的初始速度模型，然后进行常 规时间域全波形反演。迭代200次之后输出最终反演结果。
[0136] 实施例1
[0137] 根据勘探要求，将MATLAB2013a在Windows 7旗舰版系统下进行安装，并安装 MATLAB语言编写的地震数据处理软件Crews工具包。
[0138] 本发明利用Marmousi P波速度模型进行测试，Marmousi速度模型具有复杂的构 造、细微的层理和深部储油藏构造，很适合进行理论方法测试。原始Marmousi速度模型较 大，考虑到硬件内存以及计算时间因素，本发明对原始模型进行抽稀处理，对抽稀之后的 Marmousi速度模型进行截断时窗的低通滤波多尺度全波形反演方法测试。真实模型（图4 右)和线性递增初始模型（图4左）。
[0139] 模型参数如下：
[0140] 表1截断时窗的低通滤波多尺度全波形反演方法测试参数(时间域）
[0142]截断时窗反演策略参数如下：
[0143] 对观测记录利用截断时窗进行处理，设定截断时窗长度It = 0.2:0.1: 2. Is，即选 择起始截断时窗〇-〇. 2s，然后每迭代10次截断时窗的时间长度增加0.1 s，且起始0时刻不变 (图7)，直到截断时窗的总长度为2. Is。该过程共20个截断时窗，全波形反演共迭代200次。 正演只需要正演对应时窗的时间长度，这样能很大程度上节约正演所需要的计算时间。截 断时窗反演策略由浅层到深层，逐层进行反演，极大地减小了每一次反演的参数，有效避免 了周跳现象的发生，同时增强了程序的稳定性。
[0144] 低通滤波多尺度反演策略参数如下：
[0145] 为了克服全波形反演过程中陷入局部极小和周跳等问题，同时也为了提高成像质 量，用低通滤波多尺度方法进行全波形反演，设定15Hz和25Hz两个截断频率。两个截断频率 将整个全波形反演过程分成三个频段：
[0146] I、0-15Hz:低频段。在低频段进行反演可以得到地下模型的宏观信息，在低频段的 反演结果作为下一个频率段反演的初始模型，这样可以逐渐把速度模型精确地反演出来。 但是实际数据中往往缺失低频信息，陆地地震数据缺失10Hz以下的数据，低通滤波只是在 存在低频数据的情况下才能发挥作用，当缺失低频数据时低通滤波多尺度全波形反演的效 果并不是很好(图13)。所以本发明在低频段内结合截断时窗反演策略来消除因缺失低频而 导致的周跳现象。
[0147] n、0-25Hz:中频段。该频段主要起到一个中间过度的作用，低频段的反演结果作 为中频段的初始速度模型。
[0148] m、〇-end:全频段。中频带的反演结果全频带的初始速度模型，该步骤即可以得到 时间域全波形的反演结果。
[0149] 基于以上参数，在表2所示的测试环境下进行全波形反演。
[0150] 表2计算机性能测试环境
[0152] 表3全波形反演方法对比
[0153]
[0154] 从表3可以得出，常规时间域全波形反演与本发明提出的截断时窗的低通滤波多 尺度全波形反演具有明显的差别(模型误差，计算时间方面）。常规时间域全波形反演的计 算时间为101.3分钟，而截断时窗的低通滤波多尺度全波形反演只需43.6分钟就可以完成 整个计算过程。从模型误差的角度结合反演结果图（图9，10，11，12)可以看出截断时窗的低 通滤波多尺度全波形反演反演结果精度更高，反演得到的速度结果更加准确，而常规时间 域全波形反演的结果甚至出现了错误，与真实模型相差很远。
[0155] 基于以上参数，在表2所示的测试环境下进行多种影响因素测试：
[0156] 表4不同影响因素下截断时窗的低通滤波多尺度全波形反演方法测试结果
[0158] 注:表3和表4中模型反演归一化误差计算公式为
，其中vtrue代 表真实速度模型，vk代表反演结果。
[0159] 从表4和图13,14,15中可以看出，缺失低频对全波形反演有着严重的影响。当地震 数据中缺失低频的时候，常规时间域全波形反演结果（图13)，低通滤波多尺度全波形反演 结果（图14)和截断时窗全波形反演结果（图15)都出现了明显的错误，模型反演误差明显偏 高。同时，由于缺失低频的因素导致模型更新梯度计算结果不正确，使得在进行Wolf e线性 搜索步长的过程中浪费了大量的时间，这就解释了为何相同的迭代次数，当缺失低频的时 候全波形反演的计算时间要明显长于完整频带全波形反演的计算时间。为了更好地解决因 缺失低频成分带来的反演错误问题，本发明将截断时窗反演策略和低通滤波多尺度策略结 合使用，从图16可以看出，该结合策略使得全波形反演的结果得到了明显的改善。模型误差 也说明了，该策略得到的反演结果与真实模型更加接近。
[0160]从表4和图17,18,19,20中可以看出，当地震数据中含有高斯噪声的时候，反演结 果变得不清晰了。常规时间域全波形反演（图17)和截断时窗全波形反演（图18)都出现了明 显的错误，说明截断时窗反演策略虽然在一定程度上加快了反演进程，但是其抗噪能力较 弱。低通滤波多尺度全波形反演（图19)较前两者具有明显的优势，证明了低通滤波多尺度 策略具有较强的抗噪能力。将截断时窗反演策略和低通滤波多尺度反演策略结合起来使用 (图20)，不仅加快了反演的进程，而且弥补了截断时窗抗噪能力弱的缺陷。
[0161]图1是整个反演过程的流程图，从流程图可以看出首先在低频段利用截断时窗反 演策略进行全波形反演，然后再把低频段截断时窗反演结果作为初始模型用于下一步中频 段反演，最后将中频段的反演结果最为常规时间域全波形反演的初始模型。该截断时窗低 通滤波多尺度全波形反演策略，在反演精度、计算效率、低频依赖、抗噪能力等方面优于常 规时间域全波形反演。
【主权项】
1. 一种截断时窗的低通滤波多尺度全波形反演方法，其特征在于，包括以下步骤： a、 安装MATLAB2013a软件，并安装MATLAB语言编写的地震数据处理软件Crews工具包； b、 对采集的地震数据进行预处理，并将观测记录和观测系统输入计算机，进行时间域 全波形反演； c、 从观测记录中估计震源子波，该估计得到的震源子波用于正演模拟； d、 建立线性递增初始模型，该速度模型作为反演的起始值，通过计算模型的更新方向， 实现由初始模型逐渐向真实模型逼近的过程； e、 根据最小二乘法原理构建目标函数其中v表示P波速度，(Us为采集的观测记录，dw在速度模型通过正演得到的正演记录， 在求目标函数的梯度过程中，对目标函数两端的P波速度v求导，得到梯度表达式：其中pb表示由伴随震源fs向地下传播得到的波场，pf表示正演模拟波场，s表示震源的 数目，r表示检波器的数目； f、 选取截断时窗的时间长度，截取观测记录前部分地震信号，并记录截断时窗长度，该 截断时窗的时间长度作为正演模拟的时间长度； g、 按照声波方程的动力学规律，利用输入的观测系统、估计得到的震源子波和初始速 度模型进行地震波场正演，存储正传波场数据和正演记录； 根据要求设定截断时窗的低通滤波多尺度全波形反演方法相关参数，包括模型大小nz Xnx，网格距dx，dz，正演模拟时窗长度It，时间采样间隔dt，雷克子波主频fo,巴特沃斯滤 波器的截断频率f?t，每个超随机震源的最大迭代次数iter max，更新梯度最小值gtol，目 标函数精度tol，反演结果的最大值v max和最小值v min; h、 对截断时窗内的观测记录和正演记录同时进行低通滤波，确保用相同的低通滤波器 处理数据，滤波之后，得到处理观测记录和处理正演记录； i、 将处理观测记录和处理正演记录做差，得到残差数据； j、 将残差数据作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的残差反 传波场； k、 正传波场与反传残差波场做互相关，得到模型更新梯度； l、 用超记忆梯度优化算法计算模型更新方向，并通过Wolfe收敛准则寻找合适的步长； m、 判断是否满足终止条件，若满足则输出截断时窗的低通滤波多尺度全波形反演方法 的反演结果;若不满足终止条件，将反演结果作为下一个循环的初始模型，返回d步骤； n、 将m步骤输出结果作为常规时间域全波形反演的初始速度模型，然后进行常规时间 域全波形反演，迭代200次之后输出最终反演结果。
【文档编号】G01V1/28GK106054244SQ201610429648
【公开日】2016年10月26日
【申请日】2016年6月16日
【发明人】胡勇, 韩立国, 张盼, 巩向博, 张凤蛟
【申请人】吉林大学