2. 中国石油新疆油田分公司勘探开发研究院, 新疆 克拉玛依 834000
2. Research Institute of Exploration and Development, Xinjiang Oilfield Company, PetroChina, Karamay, Xinjiang 834000, China
全波形反演(Full Waveform Inversion, FWI)是利用模拟数据来匹配观测数据中所有波形信息,并获得获取地下模型参数的高分辨率的反演结果.20世纪80年代,Tarantola(1984)提出了时间域FWI方法,并通过伴随状态法求取目标函数的梯度,大大缩短了FWI梯度计算所需要的时间.20世纪90年代,Pratt将FWI推广到了频率域,提出只需要几个离散的频率就可以得到高精度的反演结果.2009年Shin等提出了Laplace-Fourier域FWI,该反演方法主要是利用Laplace域的衰减系数来实现从浅层到深层的频率域FWI,最终可以得到一个相对较好的初始速度模型,一定程度上缓解了FWI对低频成分的依赖.但是常规FWI是一个强非线性优化问题,且实际数据中常常缺失低频成分,利用局部优化算法在反演过程中很容易陷入局部极小值,这给FWI在实际应中带来了巨大的挑战(Virieux and Operto, 2009;胡光辉等,2013;胡光辉等,2014).由于地震数据采集到低频分量造价昂贵,近年来FWI方法研究主要集中于如何构建一个弱非线性目标函数来增大观测数据与模拟数据波形匹配的容错范围,使得整个反演流程具有更强的鲁棒性,同时能够尽量缓解波形匹配过程中出现错位现象(周波跳跃).
通过波形匹配的方式来缓解FWI周波跳跃主要分为三类:第一类是通过重建地震数据中的低频成分,来增大两个波形的可匹配范围,例如Bozda等(2011),Wu等(2014),Chi等(2014),Luo和Wu(2015),敖瑞德等(2015),罗静蕊等(2016),董良国等(2017)利用地震波形的宏观包络信息来获得两个波形更大的可匹配范围,同时在此基础上Hu等(2017)和胡勇等(2017)提出了解调包络和波形模态分解的方法来重建地震数据中缺失的低频分量,并直接利用解调包络的残差数据进行波动方程反向传波,避免利用链式法则求取伴随震源,打破了基于弱散射原理求取伴随震源的缺陷(Wu and Chen, 2017).Bharadwaj等(2016)通过类似于包络凸函数方式来恢复地震数据中的长波长成分,可以得到模型的宏观结构.Zhang等(2017)和Li和Demanet(2016)在原始地震数据上通过提取有效信号的脉冲信息来重建地震数据,重建以后的地震数据中含有丰富的低频成分,将其应用FWI中可以得到一个较好的初始速度模型.但是该类低频重建的方法都是通过对数据的数学变换来提取地震数据中的低频信息进而反演速度模型中的长波长成分.这些方法虽然在FWI中取得了一定的效果,但是它们所用到的挖掘低频分量的方法多来自于信号处理,与描述地震波场传播的数学物理方程没有直接的关系,因此所得到的所谓低频信息的数学物理含义不明确(陈生昌和陈国新,2016);第二类是通过波形校正技术,使得观测数据与模拟数据足够的相似,以达到缓解周波跳跃的问题.例如:Warner和Guasch(2014)提出自适应FWI方法,Zhu和Fomel(2016)提出的自适应匹配滤波目标函数和白璐等(2016)提出的最小二乘匹配滤波FWI.该类反演方法主要是通过改造地震数据的波形信息来实现观测数据与模拟数据更好地匹配.Métivier等(2017)将数学领域优化运输问题应用到FWI中,构建了一个类似于波形搬运的线性目标函数,成功避免FWI过程陷入局部极小点.该类波形校正技术使得观测数据与模拟数据的波形相差在一定的范围内,能够有效缓解FWI的周波跳跃现象;第三类是利用多尺度分解策略来获得不同尺度的地震数据信息(Bunks et al., 1995;Shin and Cha, 2009;Chen et al., 2015;Hu et al., 2015;Song et al., 2011;Qu et al., 2016),该类方法大多是通过尺度分解得到从低频到高频的数据,或者是从浅层到深层的地震数据,逐步逐层的反演思想很好地减弱了FWI的非线性.Xie (2015)利用角度域波数滤波的方法将地震波场数据按照波数的尺度进行分类,先用低波数成分恢复模型的宏观信息,再用高波速成分更新模型的细节信息,一定程度地缓解了FWI的周波跳跃现象.地震数据是一个复频率信号,是多频率相互叠加形成变化复杂的波形.当利用尺度分解策略提取低频段信号时,可以得到变化平缓的地震波形,其对应着地下模型的宏观构造信息.先反演速度模型的宏观构造信息,再反演细节信息,最终能够实现高精度的反演结果.但是该类多尺度方法在地震数据中缺失低频分量的时候仍然会出现周波跳跃现象.
常规相位反演方法仅考虑了地震波形的相位信息(Fichtner et al., 2008),忽略了地震波形的振幅信息,这样一定程度上缓解了FWI的非线性问题.但是地震波形的相位信息会随着时间的变化存在相位缠绕问题,即相位值会随着时间轴存在突变点,这给FWI的稳定性带来了巨大挑战.关于时间域瞬时相位缠绕问题Luo等(2016)提出新的解缠相位的策略,很好地缓解了地震波形相位缠绕对FWI的影响.但是基于相位FWI方法仍然存在着一些问题,例如:地震波形的相位变化与振幅变化无关,这一点使得相位FWI目标函数在强反射界面和弱反射界面具有相同的值,这将可能使得强反射界面对应的梯度与弱反射界面对应的梯度具有相同的数值(Fichtner et al., 2008),导致速度模型更新量不对称,并在后续的FWI迭代中整个反演可能会向错误的方向更新.当地震数据中存在微弱的噪声干扰,这些微弱的噪声在地震波的瞬时相位图上具有与有效信息相等的权重,这可能会使得模拟数据的相位不断向微弱噪声的相位匹配,导致整个反演错误,甚至得不到地下模型的任何信息.因此地震信号的振幅信息在反演的过程中也是决定FWI反演精度的一个重要的因素.
本文将自适应非稳态相位校正技术应用到了FWI中,并引入自适应相位校正因子,可以根据数据的需求来调整相位校正量的大小,同时在反演的过程还结合时频域目标函数以及低通滤波多尺度反演策略,成功缓解了FWI在低频缺失的情况下出现的周波跳跃现象.全文所涉及的基础理论主要有时频域非稳态相位校正和Gabor正反变换.在这些理论的基础上构建一个非稳态时频域相位校正目标函数,再通过链式法则推导FWI的伴随震源表达式,用所得到的伴随震源求得反传波场,并与正传波场进行零延迟互相关得到自适应非稳态相位校正时频域FWI的更新梯度.数值算例对比结果表明,当地震数据中含有低频信息时,低通滤波多尺度FWI与非稳态相位校正时频域多尺度FWI都可以得到很好的反演结果.但是当地震数据缺失低频时,低通滤波多尺度FWI出现了明显的周波跳跃现象,而非稳态相位校正时频域多尺度FWI仍然能够得到高精度的反演结果.本文还对震源子波不准确的情况做了测试,当震源子波存在135°相位差时,整个波形接近翻转,波形差异程度较为明显,这给FWI带来了巨大的挑战.在震源子波不准确的情况下,并与常规FWI对比测试,证明了基于自适应非稳态相位校正的时频域多尺度FWI方法对震源子波有一定的抗干扰能力.
1 方法原理当初始速度模型与真实速度模型相差较大的情况下,观测数据与模拟数据可能存在明显的相位差异,常规FWI方法无法将相位相差大于半个周期的波形正确归位,这便导致了周波跳跃问题.本文提出基于自适应非稳态相位校正的时频域多尺度FWI目标函数,该反演方法在利用振幅信息的同时还能缓解观测数据与模拟数据相位匹配过程中因大于半个周期而导致的周波跳跃问题,可为常规FWI构建一个较好的初始速度模型.
1.1 非稳态相位校正在数字信号处理领域,通常利用复地震数据来得到观测数据与模拟数据的瞬时相位差,但是这些方法直接得到的瞬时相位变化范围都是在[-π, π]之间,且都是相互缠绕在一起,即瞬时相位随着时间的变化存在突变点.为了避免因瞬时相位缠绕问题影响FWI的稳定性以及反演精度,本文利用高斯窗函数截取非稳态的地震信号,并进行分时窗分析地震数据对应频率的相位特性,同时利用时频域互相关函数来间接求得观测数据与模拟数据在时频域的相位差(Fichtner et al., 2008),并将其应用到非稳态相位校正时频域多尺度FWI中.
为了方便表示时频域地震数据的互相关,首先做出如下定义:
(1) |
(2) |
其中u(τ), d(τ)为时间域模拟数据与观测数据;g(τ)为高斯窗函数,定义为
(3) |
其中σ为调节高斯窗大小的参数,它决定高斯窗函数的宽窄,即傅里叶变换对应的时间分辨率.τ表示高斯窗函数的时间长度.因此互相关函数可以表示为
(4) |
对互相关结果进行傅里叶变换:
(5) |
其中F[·]表示傅里叶变换;*表示取复共轭;ϕu, ϕd分别表示模拟数据和观测数据的相位.两个数据的相位差可以用下式间接表示:
(6) |
其中,
从图 1a中可以看出,当两个波形只存在相位上差别的时候,利用非稳态相位校正技术可以完全把偏差的相位校正回来,让观测数据的波形变化与模拟数据保持一致.如图 1b给出的相位校正测试结果,当观测数据与模拟数据相位存在一定差异,同时波形具有一定延时的情况下,通过非稳态相位校正技术,可以让校正后观测数据的波形变化趋势更加接近模拟数据.在FWI过程中如果直接利用存在较大相位差异的两个波形匹配,且初始速度模型较差的情况下,在反演过程中很容易出现波形匹配错误的问题(周波跳跃),影响最终的反演精度.从图 1的测试结果中可以看出,我们可以在观测数据与模拟数据之间加一个相位校正的桥梁,让两个数据在FWI过程中波形相位上的差异小一些,这样使得当波形相差在半个周期内的时候,FWI就能够有效地缓解周波跳跃现象的发生.
图 2a为了便于对比,在整个地震数据上加了振幅增益,经过增益以后的地震数据在深部的波形差异也可以明显的看出来,利用非稳态相位校正以后的观测数据与模拟数据在波形上相似度更高.局部放大图如图 3中所示,尤其在低频段(图 3a)非稳态相位校正效果更加明显,且校正后的观测数据与模拟数据基本重合.图 2b为原始观测数据与模拟数据的时频域振幅残差,随着深度的增加地震波能量逐渐衰减,但是地震数据波形相位差异会随着深度的增加而逐渐加大.因此,图 2b能量和相位两者综合因素导致在0.6 s附近时频域振幅差异达到了最大值.图 2c为原始观测数据与模拟数据的相位残差,可以看出浅层直达波附近因波形相似且不受地下介质的影响,相位差异较小.同时可以看出在深部区域的相位差异与浅部相位差异基本在一个数量及上,没有明显的大小之分.这就会产生上文说的强反射界面与弱反射界面出现数量级相同的梯度,使得速度模型更新量不对称.并且在图 2c中可以明显看出相位突变的现象,这是波形相位之间相互缠绕所导致的.该相位图是对公式(6)取对数得到的,如果直接利用Δϕ=ϕu-ϕd来求取相位差,那么相位缠绕现象会更加明显.
FWI是一个不断用模拟数据去拟合实际数据的过程,最终找到一个速度模型使得两个数据残差的平方和最小.它具有成像精度高、反演模型参数准确、复杂构造成像效果好等优点.但是FWI具有较强的非线性,当地震数据中缺失低频分量时,很容易出现周波跳跃现象.为了缓解FWI过程出现的周波跳跃问题,本文利用相位校正策略来构建时频域FWI目标函数.当观测数据经过非稳态相位校正以后更接近模拟数据,减弱了两者波形上的差异,可以有效地缓解了FWI过程中波形不匹配的问题.
本文采用Gabor变换(见附录A)来得到时频域的地震数据.Gabor变换的基本思想是将信号加上一个滑动窗函数.当利用高斯窗函数来截取非稳态信号时,在这个时窗内的信号可看成是稳态的,再将加窗后的信号进行傅里叶变换,并不断地将窗函数沿着时间轴移动,依次进行傅里叶变换,该过程即为对时间域地震信号的时频分析.随着分析的进行,可得到信号的频谱在时间轴上的变化规律称其为该信号的时频谱,通过观察时频谱的特征,更能直观的反映信号频率随时间变化的规律.非稳态相位校正时频域FWI目标函数可以定义为用时频域模拟数据不断去拟合相位校正以后的观测数据的过程,则目标函数表示为
(7) |
其中
图 4是利用图 1a所示的波形计算得到的目标函数值,其中观测数据和模拟数据存在着135°的相位差,当我们直接利用观测数据和模拟数据来进行波形匹配时,常规FWI对应的目标函数如图 4中的虚线所示,其中的全局极小值并不是观测数据与模拟数据最佳匹配的位置,称之为“伪全局极小点”,当FWI陷入伪全局极小点时,目标函数不再下降,无法对模型进一步更新,这会使得FWI无法得到最终高精度的反演结果.当我们利用非稳态相位校正以后的观测数据与模拟数据进行匹配时,如图 4中实线所示,可以看出,目标函数的全局极小点就是模拟数据与观测数据完全匹配的点.通过非稳态相位校正技术,将波形中存在的相位偏转和延时等差异降到FWI方法允许的差异范围之内,这样能够更好实现波形匹配,缓解目标函数中存在的伪全局极小点的问题.
FWI过程中目标函数梯度的计算是整个反演流程的关键问题,早先基于扰动理论来计算FWI梯度的方法由于计算量巨大而无法适用于大模型.基于伴随状态方法可以将每个检波点的伴随震源同时反传到模型空间,减小了正演模拟的次数,节约了计算时间.本文利用伴随状态法来求取非稳态相位校正时频域FWI目标函数的梯度,其中核心问题是利用链式法则将目标函数对速度参数求偏导数.
(1) 常规FWI目标函数的梯度及伴随震源
FWI通过模拟数据拟合观测数据来得到地下模型高精度的速度信息,关于常规时间域FWI目标函数可以定义为
(8) |
其中u表示时间域模拟数据;d表示时间域观测数据.目标函数对模型速度参数v的偏导数有(见附录B)
(9) |
其中T表示对Jacobian矩阵的转置;常规时间域FWI的伴随震源可以表示为:fs=u-d.本文采用伴随状态法(Pratt et al., 1998)计算目标函数对速度参数的梯度,即正传波场与反传残差波场的零延迟互相关.则常规FWI目标函数的梯度最终可以表示为
(10) |
其中Pb表示由伴随震源向地下传播得到的反传波场,Pf表示由震源发出的正演模拟波场.
(2) 自适应非稳态相位校正时频域FWI梯度及伴随震源
本文利用非稳态相位校正方法,目的是减小观测数据与模拟数据在波形上的差异,缓解FWI过程中出现的周波跳跃现象,同时在地震数据缺失低频分量时能够有效地提高FWI的反演精度.将公式(6)代入到公式(7)中得到自适应非稳态相位校正时频域FWI的目标函数:
(11) |
其中
(12) |
其中
(13) |
关于梯度及伴随震源的公式详细推到过程见附录C,最终目标函数的梯度可以按照常规FWI梯度的形式简化为
(14) |
其中Pb表示由自适应非稳态相位校正时频域FWI伴随震源向地下传播得到的反传波场,Pf表示由震源发出的正演模拟波场.
1.4 基于自适应非稳态相位校正时频域多尺度FWI流程关于自适应非稳态相位校正时频域FWI的核心思想是通过相位校正后的观测数据来引导模拟数据与观测数据匹配的过程,以达到缓解FWI周波跳跃的目的.首先,根据模拟数据来对观测数据的相位进行校正,让观测数据在波形上更加接近模拟数据,当两个数据的波形相差较小的时候,通过最小二乘反演让模拟数据进一步地向观测数据靠近,并不断对速度模型进行迭代更新.由于自适应相位校正因子的作用,当模拟数据每向观测数据靠近一步,相位校正后的观测数据便会向原始观测数据靠近一步,因此在整个反演过程中模拟数据不断地小幅度地向观测数据靠近,最终实现自适应引导模拟数据与观测数据匹配的过程,并能够逐步引导FWI跳出局部极小值.通过调整自适应相位校正因子的大小,可以控制相位校正后的观测数据与模拟数据的相似程度,当给一个较大的自适应相位校正因子时,对观测数据的相位校正量较大,这时相位校正以后的观测数据在波形上更加接近模拟数据;当给一个较小的自适应相位校正因子时,这时对观测数据的相位校正量较小.因此可以通过观测数据与模拟数据的差异程度以及初始速度模型的好坏程度来调整相位校正量的大小.
本文为了更好地模拟缺失低频的观测数据,利用Butterworth高通滤波器滤除观测数据中7 Hz以下的低频信息,同时确保观测数据与模拟数据经历相同的滤波处理过程,这样能有效避免因人为处理而导致数据不匹配问题.当地震数据中缺失低频分量时,基于低通滤波类多尺度反演策略并不能从根本上解决周波跳跃问题.针对地震数据低频缺失问题,本文利用自适应非稳态相位校正时频域多尺度FWI来得到一个较好的初始速度模型(宏观模型).然后将得到较好的初始速度模型应用到常规FWI中,这样就能够满足了FWI的弱散射条件假设,最终得到高精度的反演结果.
如图 5所示,基于自适应非态相位校正时频域多尺度FWI流程中的核心步骤是利用Gabor正反变换来求取伴随震源,并将伴随震源反传到整个模型空间得到反传波场.进而反传波场与正传波场进行零延迟互相关得到速度模型的更新梯度,然后结合优化算法便可以实现自适应非稳态相位校正时频域多尺度FWI.将得到的高精度初始速度模型应用到常规时间域FWI模块中,得到最终的反演结果.
为了验证基于自适应非稳态相位校正的时频域多尺度FWI方法的正确性和有效性,本文采用抽稀的Marmousi速度模型进行测试.经过抽稀后的速度模型网格纵横大小为69×192,模型网格间距为12.5 m,横向距离为2.4 km,纵向距离为0.8625 km.震源沉放深度为50 m,共20炮数据,炮点均匀分布.检波器沉放深度为12.5 m,每个网格对应着一个检波器,共192个检波器.震源选用22 Hz主频的雷克子波,采样间隔为0.001 s,时间采样总长度为2.0 s.初始速度模型选用如图 6b所示的水平线性速度模型,其中速度变化范围从1.5 km·s-1到4 km·s-1.
Bunks等在1995年提出低通滤波多尺度策略,很好地将地震数据按照频率尺度特性拆分开来逐步反演,该反演策略先利用低频信息重建速度模型的宏观部分,再逐渐增加地震数据中的高频段信息,实现高精度反演.如图 7a所示,利用Butterworth低通滤波器(缓解滤波导致的吉普斯效应)滤除15 Hz以上的高频信息,只利用低频段的数据可以很好地构建出Marmousi速度模型的宏观信息.与图 7a相比,图 7b是利用低频段的地震数据结合了自适应非稳态相位校正技术得到的反演结果,可以看出所得到的反演结果更能体现速度模型宏观结构信息.当低频段的地震数据经过相位校正以后使得观测数据与模拟数据波形变化趋势更加接近,原先在地震数据上反映的一些细节信息随着相位校正而被去除,留下的仅仅是对应地下速度模型的宏观信息.这样就使得经过非稳态相位校正以后的FWI目标函数得到的反演结果更加能够体现出模型的宏观结构.图 8是利用图 7的反演结果作为初始速度模型进行的常规FWI结果,对比图 8a和8b可以发现两者的反演结果并没有明显的差距,这一点可以证明图 7提供的初始速度模型都能很好缓解FWI过程中的周波跳跃问题.为了进一步验证基于自适应非稳态相位校正的时频域多尺度FWI方法在重建初始速度模型和缓解周波跳跃方面的优越性,下文将对其进行缺失低频测试,并与常规多尺度FWI方法进行对比分析.
由于常规方法采集到的野外地震数据常常缺失低频成分,导致地震数据波形变化剧烈,使得在反演的过程中观测数据与模拟数据波形匹配困难,这给FWI在实际数据应用中带来了巨大的挑战.为了证明本文提出的基于自适应非稳态相位校正的时频域多尺度FWI方法在重建初始速度模型和缓解周波跳跃方面的优势,首先对地震数据利用Butterworth高通滤波器滤除7 Hz以下的低频信息(Butterworth滤波器的截断频率为13 Hz),利用缺失低频信息的地震数据分别进行常规多尺度FWI和非稳态相位校正时频域多尺度FWI测试.
图 9给出了不同自适应相位校正因子在低频段的反演结果,其中自适应相位校正因子为0时表示常规多尺度FWI结果,若自适应相位校正因子越大,代表波形中相位校正量越大,校正后的观测数据的波形变化趋势更加接近模拟数据.从图 9中可以看出当地震数据缺失低频成分,常规多尺度FWI方法在低频段没办法得到速度模型正确的宏观信息.在模型的左侧出现了明显的错误,随着自适应相位校正因子的逐渐增大,反演结果左侧错误的地方逐渐有所改进.当自适应相位校正因子为0.8时,原先左侧错误区域得到了很好地改善,最终的反演结果与含有低频成分的常规多尺度FWI低频段结果基本一致(图 7a).证明了本文提出的基于自适应非稳态相位校正的时频域多尺度FWI方法在构建初始速度模型方面有着一定的优势.
图 10是利用多尺度策略最终得到的反演结果.为了进一步证明基于自适应非稳态相位校正的时频域多尺度FWI反演结果的准确性,分别抽取了图 10中的反演速度剖面,如图红线所标注的地方.距离为0.75 km处对应的速度剖面图,如图 11所示.可以看出,当地震数据缺失低频分量时,常规多尺度FWI反演结果在浅层的速度值明显偏离真实速度值;当自适应相位校正量小于0.6的时候,对于0.75 km处的速度值都没有明显的改善;当自适应相位校正量等于0.8时,FWI反演得到的速度值与真实速度值更加接近.但当自适应相位校正量等于1时,观测数据与模拟数据的相位差被完全校正,这样的话不利于目标函数下降,很难达到反演要求的精度.综合考虑波形匹配和目标函数差异两个方面因素,同时经过数值实验测试表明,当ε=0.8时最终的反演结果相对较好.
提取速度剖面1.75 km处的速度值做对比,如图 13所示.FWI对于Marmousi速度模型的右侧反演相对容易,可以看出即使是常规多尺度FWI结果相比于Marmousi模型的右侧反演结果也相对左侧偏好.这主要是因为右侧为倾斜构造,地震波在地下传播过程中存在大角度散射,由于观测系统的设置,这些大角度散射能量恰好能被检波器接收到,所以观测数据含有右侧模型丰富的低波数成分.模型左侧大部分为水平层理构造,主要为薄层的弱反射信息,我们知道常规FWI对于反射波数据很容易发生周波跳跃的问题.由于观测数据中含有的信息成分不同,随着不断地迭代,模型右侧的反演结果就会相对于模型左侧精度高一些.但是深部区域常规多尺度FWI效果不是很好,这主要是因为深部区域模型误差较大,致使观测数据与模拟数据波形上存在较大的相位差,当对应的波形存在半个周期以上的相位差时,常规FWI方法将无法通过局部优化算法来匹配对应的波形.这使得深部区域在反演的过程中无法得到正确的模型更新量,最终在深部的反演结果不理想.随着自适应相位校正量的增加,深部反演的速度值逐渐接近真实速度值,因为经过相位校正以后观测数据与模拟数据波形相似程度增强,波形错位问题得到了缓解,直接体现在模型更新量上.这样经过不断地小幅度地模型更新以后,最终在深部区域速度值较常规FWI方法精度高一些.
为了证明反演结果与真实速度模型的相似程度,我们在反演结果速度模型中的0.75 km处设置一个震源(此处反演的速度模型差异较为明显)来正演模拟数据,并与观测数据做差进行对比.图 13是模拟数据与观测数据的残差,其中模拟数据是在最终反演得到的速度模型(图 10)上进行正演模拟得到,观测数据是在真实速度模型上正演模拟得到的.从图 13中可以看出,当相位校正因子ε=0, 0.3, 0.6时,地震数据残差的左侧存在明显的绕射现象,这主要是在图 10反演结果的左侧存在的明显错误,导致了模拟数据无法与观测数据相匹配的问题.当反演陷入局部极小值以后,没办法再利用常规FWI来收敛到全局最优解.当ε=0.8时,残差数据中能量较弱,无明显的绕射现象,证明反演结果与真实速度模型在构造和速度值方面都得到了很好地匹配.再一次证明了基于自适应非稳态相位校正的时频域多尺度FWI能够引导反演过程逐步跳出局部极小点,达到重建正确的初始速度模型来缓解周波跳跃的目的.
2.3 震源子波不准测试为了进一步证明本文提出的基于自适应非稳态相位校正的时频域多尺度FWI方法相比于常规多尺度FWI的优势,选取了不同相位的震源子波来进行测试(如图 14所示).图 14a中的真实子波是由雷克子波经过135°相位旋转后得到,两者波形相位存在较大差异.如果直接用这两个震源子波在FWI过程中进行数据匹配,即使在相同的速度模型上正演模拟,得到的地震数据也会存在明显差异.这将直接影响FWI迭代的更新方向,以及最终的反演精度.本节测试的真实震源子波和模拟震源子波只存在一个135°的相位差,所以二者具有相同的频谱图,都是22 Hz主频的震源子波,如图 14b所示.
为了更直观的看出不同震源子波在伴随震源中的表现形式,以及非稳态相位校正方法在震源不准确情况下的优势,我们假设初始速度型与真实速度模型完全相同,只保留震源这一个变化因素.图 15所示的伴随震源,是由图 14所示的真实震源子波和模拟震源子波经过有限差分正演模拟得到的(伴随震源详细推导过程见附录C).从图 15a中可以看出,即使当速度模型完全一样的情况下,真实震源子波对应的观测数据与模拟震源子波对应的模拟数据仍然有着较大的差异.如果直接利用图 15a所示的伴随震源来求取模型更新梯度,可能导致整个反演向错误的方向更新.图 15b,15c和15d分别是利用相位校正因子ε=0.3, 0.6, 0.8得到的伴随震源,其实图 15中伴随震源中的能量来自于震源子波的差异.当利用较大的相位校正因子时,校正后的观测数据与模拟数据差异变小,伴随震源的能量相对较弱.从图 15d中可以看出,虽然相位校正因子ε=0.8没能把伴随震源的能量完全校正,但相比于图 15a的波形差异已经有了明显的改善.将经过相位校正后的伴随震源利用到FWI中,能够一定程度上缓解震源子波错误引起的周波跳跃问题.
利用2.2节缺失低频测试参数并加上错误的震源子波来对本文提出的方法进行测试.测试结果如图 16所示.由于本节主要测试自适应相位校正时频域FWI方法相比于常规FWI方法的优越性,所以只对比了低频段的反演结果.图 16a结合图 15a可以看出,常规多尺度FWI在震源子波波形相差过大的时候无法得到模型的任何信息.图 16b,16c和16d分别是利用相位校正因子ε=0.3, 0.6, 0.8得到的反演结果,并结合图 15b,15c和15d我们可以看出,随着相位校正因子的增加,由不准确震源子波带来的影响逐渐减小,模型的深部构造逐渐被反演出来,浅层的干扰也逐渐得到了缓解.当相位校正因子ε=0.8时,Marmousi速度模型的宏观信息基本得到了恢复,虽然相比图 9d反演结果有一定的差距,但是与常规多尺度FWI结果(图 16a)相比,已经有了很大程度的改善.证明了基于自适应非稳态相位校正方法能够有效缓解缓解FWI过程中由震源子波不准确带来的周波跳跃问题.
(1) 本文利用Gabor变换得到时频域地震数据,进而在时频域对观测数据进行相位校正,并建立对应的时频域目标函数.从单道地震数据的非稳态相位校正结果可以看出,该方法有效避免了校正后的波形出现突变点的问题,同时能够让观测数据波形与模拟数据波形相似程度提高,增大了FWI方法在波形上的可匹配范围.数值测试结果证明了该方法能够通过非稳态相位校正后的观测数据来引导模拟数据,并逐渐帮助FWI跳出局部极小点,最终得到高精度的反演结果.
(2) 自适应非稳态相位校正时频域FWI结合了低通滤波多尺度反演策略,可以从低频到高频逐渐反演模型中的细节信息.在低频段,该方法可以得到一个较好的初始速度模型,然后再利用常规FWI方法,最终得到高精度的反演结果.数值实验测试结果表明,即使地震数据中缺失7 Hz以下的低频数据,该多尺度反演方法也能得到高精度的FWI反演结果.
(3) 自适应非稳态相位校正技术缓解了因速度模型差异引起波形相位偏差问题的同时,还能在一定程度上缓解由于震源子波不准确带来的周波跳跃现象.震源子波不准确测试结果表明,非稳态相位校正时频域FWI能够在震源子波波形相位差异较大的情况下得到速度模型的宏观信息,缓解震源不准确带来的影响.
当然该方法还可以与其它FWI技术相结合,例如:Laplace域FWI、不依赖震源子波反演策略和解调包络重建低频等方法.通过与其它新方法相互结合,能够进一步缓解FWI对初始速度模型、震源子波和低频成分的依赖性.
附录A地震观测数据和模拟数据傅里叶变换可以表示为
(A1) |
(A2) |
其中F表示傅里叶变换.地震观测数据和模拟数据傅里叶逆变换可以表示为
(A3) |
(A4) |
其中F-1表示傅里叶逆变换.地震观测数据和模拟数据Gabor变换可以表示为
(A5) |
(A6) |
其中G表示Gabor变换.地震观测数据和模拟数据Gabor逆变换可以表示为
(A7) |
(A8) |
其中G-1表示Gabor逆变换.
附录B关于常规时间域FWI目标函数可以定义为
(B1) |
其中v表示模型的速度参数,u表示模拟数据,d表示观测数据,ns表示震源的个数,nr表示检波器的个数.本文只讨论常密度声波方程:
(B2) |
其中P代表正传波场,s代表震源函数.将声波方程写成离散的矩阵形式:
(B3) |
(B4) |
根据Pratt (1998)在频率域给出的伴随状态方法,将其推广到时间域FWI伴随震源的表达式中,同时将时间域波动方程用矩阵的形式来表示出来.公式(B3)两边同时对速度参数求偏导数,则模拟数据对速度模型的导数可以表示为
(B5) |
其中P(x, t)=u表示模拟数据,将公式代入得到模拟数据对速度模型的偏导数:
(B6) |
其中(L-1)T表示伴随算子,用于将伴随震源反传到模型空间,对于二阶声波方程如公式(B2)存在(L-1)T=(L-1).T表示矩阵的转置.本文利用常密度声波方程进行数值模拟,则目标函数对模型速度参数v的偏导数有:
(B7) |
其中Jacobian矩阵定义为:
(B8) |
其中Pb表示由伴随震源向地下传播得到的反传波场,Pf表示震源发出的正演模拟波场.
附录C自适应非稳态相位校正时频域FWI的目标函数:
(C1) |
其中
(C2) |
其中为了方便用∇v来表示对速度参数求梯度.对于时频域数据绝对值的求导问题(Fichtner et al., 2008),可以有如下表示:
(C3) |
公式(C3)中,复数的偏导数与其复共轭的偏导数其虚部正好相反,因此在求和的过程中只剩下实部成分.其中Re[·]表示取计算结果的实部.按照公式(C3)的计算步骤,公式(C2)可以表示为
(C4) |
通过链式法则进一步求导得到:
(C5) |
对公式(C5)提取公因式得到
(C6) |
公式(C6)中涉及到时频域的模拟数据的复共轭对速度模型求偏导数,关于这个问题可以将其转换到时间域来进行计算,如下所示:
(C7) |
其中u(τ), g(τ-t)为实数域模拟数据和高斯窗函数,因此其复共轭与原始值保持一致,存在下式成立:
(C8) |
根据公式(C8),时频域数据的复共轭可以表示成
(C9) |
将(C9)式带入到(C6)式中,得到
(C10) |
根据公式(A7),进一步发现(C10)可以通过调整积分的先后顺序,可以进一步利用Gabor逆变换来对公式进行简化.将(A7)式代入到(C10)式中,可以将时频域模拟数据对速度参数的导数转换成时间域模拟数据对速度参数的导数:
(C11) |
进而自适应非稳态相位校正时频域FWI的伴随震源可以表示成
(C12) |
按照常规FWI的表达习惯,可将自适应非稳态相位校正时频域FWI梯度表示为
(C13) |
(L-1)Tfs表示自适应非稳态相位校正时频域FWI伴随震源的反传波场.因此目标函数的梯度可以简化为
(C14) |
其中Pb表示由自适应非稳态相位校正时频域FWI伴随震源向地下传播得到的反传波场,Pf表示震源发出的正演模拟波场.
Ao R D, Dong L G, Chi B X. 2015. Source-independent envelope-based FWI to build an initial model. Chinese Journal of Geophysics (in Chinese), 58(6): 1998-2010. DOI:10.6038/cjg20150615 |
Bai L, Han L G, Zhang P, et al. 2016. Time-domain full waveform inversion based on least square filter. Oil Geophysical Prospecting (in Chinese), 51(4): 721-729. |
Bharadwaj P, Mulder W, Drijkoningen G. 2016. Full waveform inversion with an auxiliary bump functional. Geophysical Journal International, 206(2): 1076-1092. DOI:10.1093/gji/ggw129 |
Bozda E, Trampert J, Tromp J. 2011. Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements. Geophysical Journal International, 185(2): 845-870. DOI:10.1111/gji.2011.185.issue-2 |
Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics, 60(5): 1457-1473. DOI:10.1190/1.1443880 |
Chen G X, Wu R S, Chen S C. 2015. Full waveform inversion in time domain using time-damping filters. //2015 SEG Annual Meeting. New Orleans, Louisiana: SEG, 1451-1455.
|
Chen S C, Chen G X. 2016. Full waveform inversion of the second-order time integral wavefield. Chinese Journal of Geophysics, 59(10): 3765-3776. DOI:10.6038/cjg20161021 |
Chi B X, Dong L G, Liu Y Z. 2014. Full waveform inversion method using envelope objective function without low frequency data. Journal of Applied Geophysics, 109: 36-46. DOI:10.1016/j.jappgeo.2014.07.010 |
Dong L G, Huang C, Chi B X, et al. 2015. Strategy and application of waveform inversion based on seismic data subset. Chinese Journal of Geophysics (in Chinese), 58(10): 3735-3745. DOI:10.6038/cjg20151024 |
Fichtner A, Kennett B L N, Igel H, et al. 2008. Theoretical background for continental-and global-scale full-waveform inversion in the time-frequency domain. Geophysical Journal International, 175(2): 665-685. DOI:10.1111/gji.2008.175.issue-2 |
Hu G H, Jia C M, Xia H R, et al. 2013. Implementation and validation of 3D acoustic full waveform inversion. Geophysical Prospecting for Petroleum (in Chinese), 52(4): 417-425. |
Hu G H, Wang L X, Fang W B, et al. 2014. Full Waveform Inversion Method and Application (in Chinese). Beijing: Petroleum Industry Press.
|
Hu Y, Zhang D, Yuan J Z, et al. 2015. An efficient multi-scale waveform inversion method in Laplace-Fourier domain. Petroleum Exploration & Development, 42(3): 369-378. |
Hu Y, Han L G, Xu Z, et al. 2017. Adaptive multi-step full waveform inversion based on waveform mode decomposition. Journal of Applied Geophysics, 139: 195-210. DOI:10.1016/j.jappgeo.2017.02.017 |
Hu Y, Han L G, Xu Z, et al. 2017. Demodulation envelope multi-scale full waveform inversion based on precise seismic source function. Chinese Journal of Geophysics (in Chinese), 60(3): 1088-1105. DOI:10.6038/cjg20170321 |
Li Y E, Demanet L. 2016. Full-waveform inversion with extrapolated low-frequency data. Geophysics, 81(6): R339-R348. DOI:10.1190/geo2016-0038.1 |
Luo J R, Wu R S. 2015. Seismic envelope inversion:reduction of local minima and noise resistance. Geophysical Prospecting, 63(3): 597-614. DOI:10.1111/1365-2478.12208 |
Luo J R, Wu R S, Gao F C. 2016. Time-domain full-waveform inversion using instantaneous phase with damping. //2016 SEG International Exposition and Annual Meeting. Dallas, Texas: SEG, 1472-1476. http://www.researchgate.net/publication/307880413_Time-domain_full-waveform_inversion_using_instantaneous_phase_with_damping
|
Luo J R, Wu R S, Gao J H. 2016. Local minima reduction of seismic envelope inversion. Chinese Journal of Geophysics (in Chinese), 59(7): 2510-2518. DOI:10.6038/cjg20160716 |
Métivier L, Brossier R, Mérigot Q, et al. 2017. Measuring the misfit between seismograms using an optimal transport distance:application to full waveform inversion. Geophysical Journal International, 205(1): 345-377. |
Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x |
Qu Y, Li Z, Huang J, Li J. 2018. Multi-scale full waveform inversion for areas with irregular surface topography in an auxiliary coordinate system. Exploration Geophysics: 49. |
Shin C, Cha Y H. 2009. Waveform inversion in the Laplace-Fourier domain. Geophysical Journal International, 177(3): 1067-1079. DOI:10.1111/gji.2009.177.issue-3 |
Song J Y, Zheng X D, Qin Z, et al. 2011. Multi-scale seismic full waveform inversion in the frequency-domain with a multi-grid method. Applied Geophysics, 8(4): 303-310. DOI:10.1007/s11770-011-0304-2 |
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367 |
Warner M, Guasch L. 2014. Adaptive waveform inversion: Theory. //2014 SEG Annual Meeting. Denver, Colorado, USA: SEG.
|
Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model. Geophysics, 79(3): WA13-WA24. DOI:10.1190/geo2013-0294.1 |
Wu R S, Chen G X. 2017. New Frechet derivative for envelope data and multi-scale envelope inversion. //79th EAGE Conference and Exhibition 2017. EAGE. http://www.researchgate.net/publication/317550992_New_Frechet_Derivative_for_Envelope_Data_and_Multi-scale_Envelope_Inversion
|
Xie X B. 2015. An angle-domain wavenumber filter for multi-scale full-waveform inversion. //2015 SEG Annual Meeting. New Orleans, Louisiana: SEG. http://www.researchgate.net/publication/299905983_An_angle-domain_wavenumber_filter_for_multi-scale_full-waveform_inversion
|
Zhang P, Han L G, Xu Z, et al. 2017. Sparse blind deconvolution based low-frequency seismic data reconstruction for multiscale full waveform inversion. Journal of Applied Geophysics, 139: 91-108. DOI:10.1016/j.jappgeo.2017.02.021 |
Zhu H J, Fomel S. 2016. Building good starting models for full-waveform inversion using adaptive matching filtering misfit. Geophysics, 81(5): U61-U72. DOI:10.1190/geo2015-0596.1 |
敖瑞德, 董良国, 迟本鑫. 2015. 不依赖子波、基于包络的FWI初始模型建立方法研究. 地球物理学报, 58(6): 1998-2010. DOI:10.6038/cjg20150615 |
白璐, 韩立国, 张盼, 等. 2016. 最小平方滤波法时间域全波形反演. 石油地球物理勘探, 51(4): 721-729. |
陈生昌, 陈国新. 2016. 时间二阶积分波场的全波形反演. 地球物理学报, 59(10): 3765-3776. DOI:10.6038/cjg20161021 |
董良国, 黄超, 迟本鑫, 等. 2015. 基于地震数据子集的波形反演思路、方法与应用. 地球物理学报, 58(10): 3735-3745. DOI:10.6038/cjg20151024 |
胡勇, 韩立国, 许卓, 等. 2017. 基于精确震源函数的解调包络多尺度全波形反演. 地球物理学报, 60(3): 1088-1105. DOI:10.6038/cjg20170321 |
胡光辉, 贾春梅, 夏洪瑞, 等. 2013. 三维声波全波形反演的实现与验证. 石油物探, 52(4): 417-425. |
胡光辉, 王立歆, 方伍宝, 等. 2014. 全波形反演方法及应用. 北京: 石油工业出版社.
|
罗静蕊, 吴如山, 高静怀. 2016. 地震包络反演对局部极小值的抑制特性. 地球物理学报, 59(7): 2510-2518. DOI:10.6038/cjg20160716 |