基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演-冯晅.pdf

《基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演-冯晅.pdf》由会员分享，可在线阅读，更多相关《基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演-冯晅.pdf（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 第卷 第期年月 吉林大学学报(地球科学版)Journal of Jilin University( Earth Science Edition) Vol No Nov 冯晅,鲁晓满,刘财,等基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演吉林大学学报(地球科学版) , ( ) : doi: / j cnki jjuese Feng Xuan, Lu Xiaoman, Liu Cai, et al Frequency-Domain Full W aveform Inversion of D Viscous Acoustic W ave Equa-tion Using Decreasin

2、g Random Shot Subsampling M ethod Journal of Jilin University ( Earth Science Edition) , , ( ) : doi: / j cnki jjuese 基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演冯 晅,鲁晓满,刘 财,周 超,金泽龙,张明贺吉林大学地球探测科学与技术学院,长春 摘要:全波形反演方法利用叠前地震波场的运动学和动力学信息重建地下速度结构,具有揭示复杂地质背景下构造与岩性细节信息的潜力.然而,巨大的计算量是阻碍其发展的一个瓶颈问题.为此,研究者们提出了震源编码技术来减少计算量,但是此方法

3、在模型更新过程中会引进随机串扰噪声,降低反演结果准确性.所以,在保证计算精度的情况下,本文提出了采用逐减随机震源采样的方法来高效计算全波形反演问题.笔者将此方法应用于频率域二维黏滞声波波动方程全波形反演,开始了在频率域进行随机震源采样类方法的研究,计算过程中共使用了依次增大的个频率段;并应用Overthrust模型来验证此类随机震源采样法的正确性.实验结果表明:基于逐减随机震源采样法的反演结果与实际Overthrust模型的拟合误差为 ,而应用基于全部震源的全波形反演方法得到的反演结果与实际Overthrust模型的拟合误差为 ,两者差别不大;但计算用时由 min减少到 min,即计算效率提高

4、了 倍.为了更好地确定方法的有效性,将其应用于M armousi模型进行试算.模型试算结果表明:基于逐减随机震源和基于全部震源得到的反演结果与实际M armousi模型的拟合误差分别为 和 ,相差不大;但计算用时由 min减少到 min,计算效率提高了 倍.综上,在保证反演精度的情况下,基于逐减随机震源采样法的频率域全波形反演方法大大减少了计算量,具有不可替代的计算优势,并且没有引进随机串扰噪声.关键词:逐减随机震源采样法;频率域;二维黏滞声波方程;全波形反演;无记忆拟牛顿算法; L BFGSdoi: / j cnki jjuese 中图分类号:P 文献标志码:A收稿日期: 作者简介:冯晅(

5、) ,男,教授,博士生导师,主要从事探地雷达、地震数据处理和解释方面的研究, E-mail: fengxu-an jlu edu cn基金项目:国家自然科学基金重点项目( ) ;国家重点基础研究发展计划( “ ”计划)项目( CB)Supported by the State Key Program of National Natural Science of China ( ) and the State Key DevelopmentProgram for Basic Research of China ( CB)Frequency-Domain Full W aveform Invers

6、ion of D Viscous Acoustic W aveEquation Using Decreasing Random Shot Subsampling M ethodFeng Xuan, Lu Xiaoman, Liu Cai, Zhou Chao, Jin Zelong, Zhang M ingheCollege of GeoExploration Science and Technology , Jilin University , Changchun , ChinaAbstract: The full waveform inversion ( FW I) method util

7、izes kinematic and dynamic information ofpre-stack seismic data to rebuild underground velocity structure However, the large amount of万方数据calculation is a bottleneck problem for its development Therefore, Researchers proposed the techniquesof the source-encoding in order to reduce calculation amount

8、, but this method would introduce randomcrosstalk in model updates The paper uses decreasing random shot subsampling method to invert fullwaveform effectively in the case of guaranteeing the calculation accuracy The paper applied the methodinto frequency-domain full waveform inversion of D viscous a

9、coustic wave equation, and research thetype of random shot sampling method in the frequency domain W e totally use eight gradually increasedfrequencies in FW I process, and apply overthrust model to demonstrate the validity of such random shotsampling method Through the experiment, we can see that t

10、he fitting error between inversion resultsbased on decreasing random shot subsampling method and the actual overthrust model is , andthe fitting error between inversion results based on all the shots and the actual overthrust model is Namely the inversion results of the two methods have little diffe

11、rence However, thecalculated time reduces from min to min Namely the computational efficiency increased times In order to better determine the validity of the method, we apply it to M armousi model Thefitting error between inversion results based on decreasing random shot subsampling method and thea

12、ctual M armousi model is , and the fitting error between inversion results based on all the shotsand the actual M armousi model is Namely the inversion results of the two methods have littledifference However, the calculated time reduces from min to min Namely thecomputational efficiency increased t

13、imes In conclusion, under the condition of the inversionaccuracy, the full waveform inversion based on the decreasing random shot sampling method infrequency domain greatly reduces the amount of calculation Namely it has irreplaceable computingadvantage Therefore there is no introduction of random c

14、rosstalk noiseKey words: decreasing random shot subsampling method; frequency-domain; D viscous acousticwave equation; full waveform inversion; memory-less quasi-newton algorithm; L BFGS 引言在地震数据采集过程中,为了得到更多的地下构造和岩性信息,通常采用多炮激发多检波器接收的观测系统.但是,实际数据处理的时候,有时不需要如此稠密的震源采样 .所以,应该寻求合适的方法,在保证计算精度的情况下减少采用的震源数,即

15、减少每次迭代的正演次数,减少计算量,提高计算效率.众所周知,全波形反演计算量巨大.因为一次模型迭代更新至少需要三次正演,而正演计算量不容小觑 .所以,全波形反演方法一度因为计算量问题而停滞不前.计算机技术的发展缓解了计算量问题,给全波形反演带来了曙光,使其进入了快速发展时期.至今多年,全波形反演取得了许多突破性成果 .但是,计算量大的问题仍然没有得到彻底解决,至今仍是研究人员需要突破的关键性问题.国内研究人员主要从算法的收敛性和GPU(众核处理器)技术来加快全波形反演的速度.例如,邓武斌等 提出将L BFGS( limited-memory Broyden-Fletcher-Goldfarb-

16、Shanno)算法和SPGL( spectralprojected gradient for L minimization)算法联合应用于全波形反演中;苏超等 使用GPU并行计算技术对数据进行处理;龙桂华等 在反演过程中用Hessian矩阵的对角线元素来做梯度类方法的预条件算子,能够吸收高斯牛顿法的二次收敛优势,使得算法具有较快的收敛速度;刘璐等 提出利用一种新的拟牛顿公式对DFP ( Davidon-Fletcher-Powell)和BFGS算法进行修正,改进后的BFGS算法在近似Hessian矩阵逆矩阵时,不仅考虑了梯度和模型信息,还加入了目标函数本身的信息,在保证反演精度的同时,明显提高

17、了反演效率;成景旺等 提出伪Hessian矩阵的预处理方法,该方法吸收了高斯牛顿法的二次收敛优势,在不增加计算量的前提下,加快收敛速度.国外学者为减少计算量也提出了很多有效的方法. Li等 提出使用压缩感知方法来加快全波形反演的速度.震源编码技术通过同时激发震源成为典型的减少计算量的方法之 吉林大学学报(地球科学版) 第卷 万方数据一 ;但是,其在加快运算速度的同时引进了震源之间的串扰噪声,从而降低了反演的精度 .为此,研究人员为了发挥编码技术的优点,发现编码技术的运用效果取决于超级炮中编码的震源数 .所以,在应用此技术的时候要注意选择合适的震源数参与编码.为了克服串扰噪声,对合适的震源编码步

18、骤,如随机编码方案 等进行研究是必然的.Schuster等 提出相编码技术来压制串扰噪声.全部震源参与编码较部分震源参与编码有效率,但它也是对串扰噪声最敏感的方法 .值得注意的是,应用L范数时,震源编码技术不适于拖拽检波器数据.Choi和Alkhalifah 引进了卷积波场并发展了应用互相关目标函数的震源编码方法来反演海洋检波器数据.本文是关于应用逐减随机震源采样法减少全波形反演计算量问题的研究.减少每次迭代中参与运算的震源数量就会增加全波形反演的效率.随机震源采样方法,即每次迭代时从全部震源中随机抽取一部分震源参与运算,每次迭代所用震源不同,获得的信息几乎包含全部震源参与运算时获得的信息,这

19、样就在保证了计算精度的同时提高了计算效率.许多学者已成功将随机震源法应用于全波形反演中并取得了良好效果. 年, M oldovean 在D模型中实现了时间域规则和随机震源采样方法;年, Ha等 提出循环震源采样方法并应用于Laplace域声波方程; 年, W ang等 提出不规则震源采样方法并应用于基于VTI ( verticaltransverse isotropy)介质的声波波动方程.相较于震源编码技术,逐减随机震源采样法不会产生串扰噪声.但是,若是采用周期采样方法会产生信号混叠问题.所以,在选取采样方法的时候一般偏向于选择逐减随机震源采样法.因为它可以将混叠噪声转化为低级随机噪声 .本文

20、提出了逐减随机震源采样法,并将其应用于频率域二维黏滞声波波动方程中;还将其与无记忆拟牛顿算法、 L BFGS算法结合,以证明无记忆拟牛顿算法和L BFGS进行全波形反演过程中应用逐减随机震源采样法的可行性. 方法原理 逐减随机震源采样方法本文介绍的逐减随机震源采样方法如图所示.即每次迭代时,只从所有震源中选取一定数量的震源,每次选取时,每个震源被选中的概率是一样的;但是,这是概率事件,或许会发生使用某些震源次数比另一些多的情况.随机震源采样方法加入了随机因素,使得到的地下构造和岩性信息具有随机性,即每次迭代会得到一些不同的信息,直至迭代完毕使用了几乎所有得到的地下信息,使得在降低计算量的同时亦

21、没有失去太多地下信息,保证了反演精度几乎不变.随着迭代次数的增加,逐减随机震源采样方法参与计算的震源数逐渐减少.此方法在模型迭代开始时参与计算的震源数目较多,这样可以获取更多信息,得到好的反演模型,为下一次迭代奠定好的初始模型基础.黑点表示被选中的震源.图 逐减随机震源采样策略Fig Decreasing random shot subsampling methodschemes 频率域全波形反演之所以选择在频率域进行基于逐减随机震源采样法的全波形反演,是因为频率域反演具有计算高效、频率选择灵活、易引入衰减因子、适于多炮计算等优点 . 频率域黏滞声波方程正演虽然对二维声波方程模拟波在介质中传播

22、的研究已经很多,但是实践证明,二维黏滞声波方程能更好地模拟波的传播;因为其加入了品质因子,考虑了地层的黏滞吸收作用,更好地贴合了实际介质.所以本文研究的波动方程为频率域二维黏滞声波波动方程: Kx, z( )P x, z, ( ) x x, z( ) P x, z, ( ) x z x, z( ) P x, z, ( ) z x, z, ( ) . ( )式中: K x , z( ) x , z( )v x , z , ( ) ,是体积模量; v x , z( ) v x , z( ) iQ x , z( ) ,是复速 第期 冯 晅,等:基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演

23、万方数据度, v x , z( )是声波方程中的实速度; Q x , z( )是品质因子; x , z( )是密度; P x , z , ( )是波场值;是圆频率; x , z , ( )是震源项.所用边界条件为PM L(完全匹配层)吸收边界条件,正演模拟方法为交错网格有限差分法.对方程( )进行差分离散后可得到简化后的形式:AP F . ( )式中: A是阻抗矩阵,取决于模拟频率和介质性质;P是波场向量; F是震源向量.求解矩阵方程( ) ,可得到对应于圆频率的波场值. 基于逐减随机震源采样法的频率域全波形反演流程本文在频率域进行了随机震源采样类方法的实现.在此给出频率域逐减随机震源采样法的

24、全波形反演流程(图) .其中,逐减随机震源采样法应用于“计算目标函数对模型参数的梯度和海森矩阵”这一步.每个震源计算出来的梯度的加和即为最终得到的梯度值.每次迭代过程中,通过随机震源采样方法选取一部分震源参与运算,求得梯度值的一个近似,然后用此近似梯度值参与到后面的运算. 全波形反演原理本文将全波形反演中的目标函数设置为C m( ) i P Ti S d, i P i i P obsi . ( )式中: C m( )是反演的目标函数; P i是第i炮模拟波场 P cali和理论波场 P obsi的差; S d, i是第i炮与偏移距相关的权系数; T代表矩阵转置; 代表矩阵的复共轭; m为待反演

25、参数.目标函数关于反演参数的一阶导数如下:C m( ) i Re P cali( )T A m TA S d, i P cali( ) i P obsi .( )本文使用拟牛顿类局部优化算法求取搜索方向: m H C m( ) . ( )式中: H是海森矩阵; m是模型扰动量.由于使用的是拟牛顿类优化算法,那么H 是通过最近几次迭代得到的模型参数值和梯度值构造而得,满足拟牛顿条件即可.本文所用拟牛顿方法为无记忆拟牛顿算法, Liu等 首次将其应用于全波形反演中,其原理及公式可以参考文献 ,在此不再赘述.图 基于逐减随机震源采样法的频率域全波形反演流程Fig Frequency-domain f

26、ull waveform inversion usingdecreasing random shot subsampling method得到模型扰动量后,更新模型可通过下式得到:m k ( ) m k( ) k m k( ) . ( )式中, k是迭代步长,可通过多种方式得到,本文是采用抛物线拟合方式求取得到的. 模型试算本文分别截取了Overthrust模型和M armousi模型具有典型构造特征的部分作为试算模型,如图所示.截取的模型大小均为 .其中:Overthrust模型的纵横网格间距均 m;M armousi模型的纵横网格间距均为 m.观测系统为:在z m深度处布设炮震源,Over

27、thrust模型和M armousi模型炮间隔分别为m和 m;检波器布设在地表,共个,两模型各检波器之间的距离分别为 m和 m;震源位置 吉林大学学报(地球科学版) 第卷 万方数据与检波器位置一一对应,分布在网格点上;两模型最大偏移距分别为 km和 km.本文所用随机方法为:在Overthrust模型中,参与计算的震源最多个,最少个,中间随着迭代次数增加而从递减;在M armousi模型中,参与计算的震源最多个,最少个,中间随着迭代次数增加而从递减.经过大量实验发现,两模型各自选择震源的方式达到的效果最好.本文在频率域进行反演.震源主频为 Hz.反演频率从小到大依次为 , , , , , ,

28、, Hz,共个频率;每个频率最大迭代次数为次.迭代终止条件为反演速度参数值变化率小于 .初始模型由理论模型通过高斯平滑得到,如图所示. Overthrust模型图为利用Overthrust模型、基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法反演得到的速度模型.比较图a、 b可以看出,针对Overthust模型,两种方法反演得到的速度模型几乎一致,没有差别.比较模型中具有典型构造特征的部分(圈出来的部分)可以看出,两方法均反演出了断裂带和推覆构造这样的构造分布信息.并且,由图b可以看出,逐减随机震源采样法反演得到的速度模型没有引进随机串扰噪声.图 截取的部分Overthust

29、(a)和M atmousi(b)理论速度模型Fig Capture parts of the Overthust(a)and the M atmousi(b)theoretical velocity modela Overthust; b M atmousi.图 用于反演的初始速度模型Fig Starting velocity models for the inversions 第期 冯 晅,等:基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演万方数据利用Overthrust模型对两方法进行试验比较.图 基于全部震源(a)和逐减随机震源采样(b)的全波形反演方法反演得到的速度模型Fig

30、 Inversion velocity models obtained using all shots(a)and decreasing random shot subsampling method(b)为了更直观地比较两种方法的反演精度,从两种方法反演得到的Overthrust模型和理论Overthrust模型中抽取地表第个点和第个点处、即位于水平位置x m处和x m处的纵向反演曲线进行对比,结果如图所示.由图中圈出来的部分可以看出,两种方法对于模型上部的反演都很精确,对于底部速度突变点处,利用全部震源反演得出的速度值较精确些,但是相差不大,可以接受.现在从数值上精确确定基于全部震源的全波形

31、反演方法和基于逐减随机震源采样的全波形反演方法的计算精度.本文利用反演得到的模型参数值与理论模型参数值的相对误差M err来判定反演模型的精确度,其计算公式为M err N m cal m truemtrue. ( )式中: m cal是反演模型参数; m true是理论模型参数;N是模型的总网格点数.本文所用计算机运行内存为 G, CPU为运行频率为 GHz的Intel( R) Xeon( R) ,所有计算均使用同一台计算机.通过对Overthrust模型进行试算,并结合式( ) ,得到了基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法的总计算用时和相对误差(表) .从表可

32、以看出,逐减随机震源采样法频率域全波形反演在保证反演精度的情况下,能快速减少计算量,较之前的计算效率提高了 倍.理论速度; 逐减随机震源采样法反演速度; 利用全部震源反演速度.图 针对Overthrust模型的地表第个点(a)和第个点(b)处反演速度曲线对比Fig Inversion speed curve contrast of surface th point(a)and th point(b)for Overthrust model表 两方法的总计算用时和相对误差结果Table Total time and relative error results of the twomethodO

33、verthrust模型总计算用时/ min相对误差基于全部震源的全波形反演方法 基于逐减随机震源采样的全波形反演方法 吉林大学学报(地球科学版) 第卷 万方数据利用M armousi模型对两方法进行试验对比.图 基于全部震源(a)和逐减随机震源采样法(b)的全波形反演方法反演得到的速度模型Fig Inversion velocity models obtained using all shots(a)and decreasing random shot subsampling method(b)本文将无记忆拟牛顿算法与随机震源采样法结合并将其应用于Overthrust模型进行试算,取得了不错的

34、反演结果. M armousi模型图为利用M armousi模型、基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法反演得到的速度模型.比较图a、 b可以看出,针对M armousi模型,两种方法反演得到的速度模型几乎一致.比较模型中具有典型构造特征的部分(圈出来的部分)可以看出,两方法均反演出了构造断裂带和速度突变点等信息.并且,由图b可以看出,逐减随机震源采样法反演得到的速度模型没有引进随机串扰噪声.为了更直观地比较两方法的反演精度,从两种方法反演得到的M armousi模型和理论M armousi模型中抽取地表第个点和第个点处、即位于水平位置x m处和x m处的纵向反演曲

35、线进行对比,结果如图所示.在图中,z m处反演速度拟合的误差很大.这是因为,观测系统布设在地表,能量随深度逐减,使得处于地表的检波器记录到的深部信息量减少,反演时信息不足导致不能准确得出此处精确速度值;并且,此处速度强烈突变,不利于反演.除此之外,两种方法对于模型上部和下部的反演都很精确.通过M armousi模型试算,并结合式( ) ,得到了基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法的总计算用时和相对误理论速度; 逐减随机震源采样法反演速度; 利用全部震源反演速度.图 针对M armousi模型的地表第个点(a)和第个点(b)处反演速度曲线对比Fig Inversio

36、n speed curve contrast of surface th point(a)and th point(b)for M armousi model差(表) .从表可以看出,逐减随机震源采样法频率域全波形反演在保证反演精度的情况下,能快速减少计算量,较之前的计算效率提高了 倍.在M armousi模型试算中,本文使用的优化算法是L BFGS方法,也取得了满意的反演结果.综上,从两模型的实验结果可以看出,随机震源采样方法具有不可替代的计算优势. 第期 冯 晅,等:基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演万方数据表 两方法的总计算用时和相对误差结果Table Total

37、time and relative error results of the twomethodM armousi模型总计算用时/ min相对误差基于全部震源的全波形反演方法 基于逐减随机震源采样的全波形反演方法 结论和展望)模型试算的结果表明:基于逐减随机震源采样法的频率域全波形反演可以通过有限的频带反演出精度很高的速度值;此方法合理地减少了参与计算的震源数,即减少了正演次数,提高了计算效率,缓解了全波形反演计算量大的问题.)通过分析反演得到的速度图,可以看出该方法没有引进随机串扰噪声,即该方法具有一定的抗噪能力,表现出较好的可靠性、有效性和稳定性.)逐减随机震源采样法与无记忆拟牛顿算法、L

38、 BFGS算法结合完成了全波形反演工作,证明该方法可以与各种反演优化算法结合,说明其稳定可靠;利用Overthrust模型和M armousi模型进行试算,说明该方法适用模型不单一,再次证明了该方法的稳定可靠.基于随机震源采样法的全波形反演用于D模型所能提高的效率有限,预计若将其应用于D模型计算效率会大幅度提高.本文所做研究只是为D模型提供借鉴及经验.在将来的发展中,研究如何利用最少的震源获取最多的地下地质信息是今后努力的方向.参考文献(References): Herrmann F J Randomized Sampling and Sparsity:Getting M ore Inform

39、ation from Fewer Samples J Geophysics, , ( ) : W B W B M ilton A, Trickett S, Burroughs L Reducing Ac-quisition Costs with Random Sampling andM ultidimensional Interpolation J / / SEG TechnicalProgram Expanded Abstracts, , ( ) : 杨庆节,刘财,耿美霞,等交错网格任意阶导数有限差分格式及差分系数推导 J 吉林大学学报(地球科学版) , , ( ) : Yang Qingj

40、ie, Liu Cai, Geng M eixia, et al StaggeredGrid Finite Difference Scheme and CoefficientsDeduction of Any Number of Derivatives J Journalof Jilin University ( Earth Science Edition) , , ( ) : Tarantola A Inversion of Seismic-Reflection Data inthe Acoustic Approximation J Geophysics, ,( ) : Sirgue L,

41、Pratt R G Efficient W aveform Inversion andImaging: A Strategy for Selecting TemporalFrequencies J Geophysics, , ( ) : 邓武兵,韩立国,李翔,等 L-BFGS和SPGL在全波形反演中的联合应用 C / /中国地球物理学会第二十七届年会长沙:中国地球物理学会, : Deng W ubing, Han Liguo, Li Xiang, et al Appply theCombination of L-BFGS and SPGL to Full-W aveformInversion

42、C / / The th Annual M eeting of ChineseGeophysical Society Changsha: Chinese GeophysicalSociety, : 苏超,周辉,林鹤时间域声波全波形反演及GPU加速 C / /中国地球物理学会第二十八届年会北京:中国地球物理学会, : Su Chao, Zhou Hui, Lin He Acoustic Full W aveformInversion in Time-Domain and Its Acceleration byGPU C / / The th Annual M eeting of ChineseG

43、eophysical Society Beijing: Chinese GeophysicalSociety, : 龙桂华,李小凡,张美根,等频率域黏弹性声波透射波形速度反演 J 地震学报, , ( ) : Long Guihua, Li Xiaofan, Zhang M eigen, et al Vis-coacoustic Transmission W aveform Inversion forVelocity Structure in Space-Frequency Domain J ActaSeismologica Sinica, , ( ) : 刘璐,刘洪,张衡,等基于修正拟牛顿公式

44、的全波形反演 J 地球物理学报, , ( ) : Liu Lu, Liu Hong, Zhang Heng, et al Full W aveformInversion Based on M odified Quasi-Newton Equation J Chinese Journal of Geophysics, , ( ) : 成景旺,顾汉明,刘春成,等频率域反射波全波形速度反演 J 地球科学:中国地质大学学报, , ( ) : Cheng Jingwang, Gu Hanming, Liu Chuncheng, etal Full W aveform Inversion for Velocity Structurefrom Reflected W ave Seismic Data in the Frequency-Domain J Earth Science: Journal of ChinaUniversity of Geoscience, , ( ) : Li X, Aravkin Y A, van Leeuwen T Fast Rando- 吉林大学学报(地球科学版) 第卷 万方数据mized Full-W aveform Inversion with Compres