基于照明预处理的分步多参数时间域声波全波形反演方法研究

张广智,姜岚杰,孙昌路,黄义双

(1.中国石油大学(华东),山东青岛266580;2.海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛266235;3.中国天辰工程有限公司,天津300400;4.中国石油化工股份有限公司勘探分公司,四川成都610041)

基于照明预处理的分步多参数时间域声波全波形反演方法研究

张广智1,2,姜岚杰1,2,孙昌路3,黄义双4

(1.中国石油大学(华东),山东青岛266580;2.海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛266235;3.中国天辰工程有限公司,天津300400;4.中国石油化工股份有限公司勘探分公司,四川成都610041)

密度是地震勘探中最重要的信息之一,在岩性解释、储层流体预测等多个方面起到不可替代的作用。但在全波形反演(FWI)中,由于密度和速度串扰的影响,很难反演出理想的密度信息。针对该问题,采取分步多参数全波形反演的策略,将反演的高精度速度结果作为初始速度模型,联合初始密度模型进行下一步多参数同时反演,通过提高初始速度模型的精度,得到更为准确的密度结果。同时,为了进一步平衡梯度能量,减小地震波传播过程中几何扩散的影响,采用照明预处理L-BFGS法提高反演精度。模型测试结果表明,分步多参数全波形反演及照明预处理L-BFGS法能很好地提高反演精度。

多参数全波形反演;分步反演;L-BFGS法;照明分析;声波介质

随着油气勘探开发复杂程度的不断加深,如何提高储层预测的精度逐渐成为专家学者研究的重点,而寻找一种适用于复杂介质的反演方法尤其重要。全波形反演(FWI)利用叠前地震波场中的运动学及动力学信息反演得到地下岩石物理特性,能够解释复杂地质背景下的构造细节信息[1-2],自1984年[3]问世后得到了飞速发展。PRATT等将全波形反演的思想引入频率域,推导了频率域内最速下降法全波形反演梯度公式[4-5]。SHIN等提出Laplace-Fourier域全波形反演方法,通过引入复频率,在缺少低频的情况下反演出较好的结果[6]。此外,SIRGUE等将全波形反演应用于混合域,得到了较好的反演结果[7]。在这些方法中,时间域全波形反演可对数据进行灵活预处理,适用于多节点并行运算,已被成功地应用于实际地震数据处理[8]。

在单参数全波形反演中,常假设密度随空间位置变化缓慢,从而将密度设为一已知恒定值,反演得到速度结果。但在实际地层中,密度是指示储层流体的重要信息,在油气解释中起着重要作用。密度的变化会引起反射率等变化,影响地震波传播及反演结果[9],如何从地震数据中得到有效的密度信息越来越成为精细化地震勘探的研究重点。多参数全波形反演是一种强非线性问题,且速度与密度之间存在串扰现象[10],严重影响反演结果的精度,因此,选择合适的反演方法和策略是重建高精度参数模型的关键。FORGUES等通过辐射模式分析指出,很难从全波形反演中得到密度[11]。JEONG等提出频率域密度反演策略,先将密度设置为常数,反演得到弹性参数,再进行速度、密度联合反演,提高了反演精度[12]。PRIEUX等首先采用长偏移距数据反演速度,然后利用全偏移数据反演速度和密度,提高了密度反演的稳定性,但是密度的反演结果仍然不是很理想[13]。董良国等通过分析目标函数随密度和速度扰动的变化规律,指出全波形反演对密度的较小尺度摄动具有较好敏感度,但是无法反演中长波长的密度摄动[14]。RAKNES等讨论了弹性波全波形反演过程中反演参数的组合方式,并提出多参数全波形反演策略——逐一反演各个参数[15]。杨积忠等提出分步变密度声波方程多参数反演策略,得到了比较好的速度和密度反演结果[16]。本文将该方法运用于时间域声波全波形反演中，利用照明预处理L-BFGS方法使不同深度的梯度能量分布更为均衡，进一步提高模型参数收敛速度和反演精度。

1 照明预处理L-BFGS

全波形反演是Bayes估计理论在勘探地球物理领域的一个应用范例[17],它将波动方程正演模拟得到的地震记录与实际观测地震记录进行匹配,通过使数据之差最小建立目标函数,寻找最佳模型参数。时间域全波形反演算法在实现过程中主要包括两个过程:以地震子波为震源的正演过程和以地震记录残差为震源的波场反传过程[18-19]。其目标函数可以写为:

(1)

式中:E为误差泛函;δd为波场残差。模型参数m可以通过下式进行迭代更新:

(2)

式中:αk是迭代步长;k是迭代次数;hk是迭代方向。全波形反演的目标函数高度非线性,在迭代过程中存在强烈的非线性问题,因此要寻找合适的优化方法提高迭代的稳定性。

牛顿类优化算法中通过对二阶偏导海森算子的计算,均衡梯度算子修正量,可以加速迭代收敛效率。但是海森矩阵占用大量的存储空间和计算时间,对计算机的存储能力和计算能力都提出了很高的要求,导致大规模全波形反演无法开展。1989年,LIU等提出的L-BFGS方法[20]是对牛顿类优化方法的改进。该方法对海森矩阵进行近似代替,不直接保存海森矩阵,而是在迭代过程中保存前n次迭代的参数修正值及梯度信息,大大减小了存储量及计算量[18,20-22]:

(3)

利用近似海森矩阵对梯度的预处理可以校正波场几何扩散对不同深度梯度的影响,增加参数在模型深部的校正量,同时降低浅层不合理的低频修正量,提高建模、成像的分辨率。但是,在L-BFGS法前n次迭代过程中,不能计算得到近似海森矩阵,梯度的能量不平衡,影响反演结果。全波形反演通过正传波场和反传波场在时间上的零延迟互相关来计算模型参数的更新量,因此梯度能量的不均衡很大程度上是由地震波传播过程中几何扩散的影响引起的。而照明分析可以研究地震波在地下介质中传播的能量分布情况,进行能量补偿[23-24],同时,全波形反演需要大量的计算量和存储空间,为了在得到较好反演效果的同时尽量减少全波形反演的计算负担,我们采用单向照明强度近似代替双向照明强度的方法,通过正演模拟得到模型各个点的波场,将其作为该点单向照明强度。但检波器接收到的地震波是双程旅行的过程[23],为了便于计算,我们将模型中各个点波场值的平方作为双向照明强度,对前n次L-BFGS法计算的梯度值进行预处理。则N炮照明总强度IS可用下式表示:

(4)

式中:ui(x,z)为t时刻在模型坐标(x,z)处的波场值。因此,参数m的前n次迭代可以用下式来表示:

(5)

2 多参数全波形反演

时间域二阶变密度声波方程如下:

(6)

式中:κ表示体积模量;ρ表示密度;p表示压力场;xs为震源坐标;s为震源函数。则速度和密度的梯度可分别用下式表示:

式中:pf是正传波场;pb是反传波场;v表示速度。

多参数全波形反演因速度和密度之间的串扰影响非常不稳定,难以获得精确的反演结果[10]。速度的变化会引起地震波传播的旅行时和振幅发生扰动,而密度的变化会影响地震波的振幅大小,且随着角度的变大,波场变化逐渐减小[9]。同时,地震反射波对密度的较小尺度摄动非常敏感,而中长波长的密度摄动对反射波基本没有影响[16]。因此,多参数全波形反演虽然可以重建较好的速度模型,但只能重建密度

的高频成分[9,14,16]。与此同时,全波形反演,尤其是多参数全波形反演,对初始模型的依赖性较强,初始速度模型越准确,密度反演结果越精确。

基于以上分析,我们将杨积忠等提出的分步反演策略[16]用于多参数时间域全波形反演,得到速度和密度反演结果。第一步是利用多参数全波形反演得到较为精确的速度模型,此时密度反演结果较差;第二步是将第一步反演的速度结果作为初始速度模型,结合初始密度模型进行速度、密度联合反演,此时由于初始速度模型精度的提高,密度反演的精度也有很大提高。为了进一步提高反演精度,可将第二步得到的速度反演结果作为初始模型,利用初始密度模型再进行一轮速度、密度联合反演。

3 模型测试

3.1 照明预处理L-BFGS方法测试

为了测试照明预处理L-BFGS方法的有效性,我们将其应用于Marmousi模型，进行速度参数的全波形反演,并与最速下降法和L-BFGS法进行对比。实际速度模型如图1a所示,网格点数为260×130,网格间距10m。地表放置38个震源,震源间距70m,并放置260个检波器,检波器间距10m。所用震源为主频30Hz的雷克子波。图1b为初始速度模型。

图2是最速下降法、L-BFGS法、照明预处理L-BFGS法速度反演结果,可以看到:最速下降法反演效果不理想,虽然浅层速度反演得很好,但是深层速度值偏差较大;L-BFGS方法在反演出较好的浅层速度的基础上,由于自身可以校正波场几何扩散的影响,深部速度的校正量增加,速度反演更为准确;而经过照明预处理后的L-BFGS方法进一步弥补了波场传播时深层能量的缺失,得到了更好的深部速度反演结果。图3是三种方法在不同水平位置处的反演速度随深度变化的曲线,可以看到照明预处理L-BFGS法得到的速度结果与真实值拟合最好,反演精度最高。图4是三种反演方法的目标函数收敛曲线,可以看到照明预处理L-BFGS法收敛最快,失配函数最小,说明照明预处理L-BFGS方法可以提高收敛速度及精度,减小运算量及运算时间。

图1 实际速度模型(a)和初始速度模型(b)

图2 三种方法速度反演结果比较a 最速下降法; b L-BFGS法; c 照明预处理L-BFGS法

图3 水平方向900m(a)、1200m(b)、1700m(c)、2250m(d)处不同方法反演的速度随深度变化曲线

3.2 分步多参数全波形反演测试

比较最速下降法、传统L-BFGS法及照明预处理L-BFGS法的效果,我们选择反演精度最高的照明预处理L-BFGS方法,采用与单参数反演模型相同的参数,进行速度、密度的联合反演测试。图5a是实际速度模型,图5b是用Gardner公式和实际速度模型计算得到的实际密度模型;图6a和图6b是初始速度模型和初始密度模型,图7a和图7b是速度和密度第一步反演结果。从图7可以看出,速度和密度的反演精度都不高,浅层速度反演较好,构造也得到很好的重建,但深部速度和实际值仍有很大偏差,密度的反演结果精度不高,不能满足反演精度的要求。

为了得到更好的反演结果,我们舍弃密度反演结果,将第一步的速度反演结果(图7a)作为初始速度模型,联合初始密度模型再进行一次多参数同时反演(第二步反演),结果如图8所示。可以看到,速度反演结果精度有了明显提高,中深层速度值与真实值吻合较好,密度反演精度较之前也有了提高,证明分步多参数全波形反演在提高反演精度方面是有效的。但是,受到初始速度模型精度的制约,密度精度的提高有限(图8b)。要进一步重建更为精确的密度模型,需要更为准确的初始速度模型,因此我们将第二步的速度反演结果作为初始速度模型,结合初始密度模型,又进行了一次多参数联合反演(最终反演),结果如图9所示。由图9可以看到，速度和密度反演的精度都有了明显改善。比较分步多参数全波形反演过程中不同水平位置处速度、密度值纵向剖面图(图10), 可以看到速度、密度联合反演的精度每一次都在提高,最终的模型反演结果与实际模型值拟合较好。

图4 不同反演方法的目标函数收敛曲线对比

图5 实际速度模型(a)和实际密度模型(b)

图6 初始速度模型(a)和初始密度模型(b)

图7 照明预处理L-BFGS分步多参数全波形反演第一步反演结果a 速度; b 密度

图8 照明预处理L-BFGS分步多参数全波形反演第二步反演结果a 速度; b 密度

图9 多参数全波形反演最终反演结果a 速度; b 密度

图10 分步多参数全波形反演速度(a,b)、密度(c,d)随深度变化曲线

4 结论

在多参数全波形反演中,由于速度和密度的耦合关系影响了反演的精度和效率,密度反演结果很不理想,选择合适的反演算法和反演策略对全波形反演至关重要。将分步多参数全波形反演策略应用于时间域多参数声波全波形反演中,通过建立更为准确的速度初始模型,减小速度和密度间的串扰影响,得到了更为准确的密度反演结果。L-BFGS是一种平衡梯度能量十分有效的方法,但是得到的反演结果分辨率仍有上升空间,通过对L-BFGS法前n次的计算结果进行照明预处理,可以在不增加过多计算量的基础上更有效地平衡模型修正量的能量,提高反演精度和收敛效率,减少反演计算量和计算时间。

[1] 杨午阳,王西文,雍学善,等.地震全波形反演方法研究综述[J].地球物理学进展,2013,28(2):766-776 YANG W Y,WANG X W,YONG X S,et al.The review of seismic full waveform inversion method[J].Progress in Geophysics,2013,28(2):766-776

[2] 王西文,杨午阳,吕彬,等.把握世界物探技术发展现状,促进地震资料处理、解释技术进步[J].地球物理学进展,2013,28(1):224-239 WANG X W,YANG W Y,LU B,et al.Grasp of the world geophysical technology development status,to promote seismic data processing and interpretation of technological progress[J].Progress in Geophysics,2013,28(1):224-239

[3] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[4] PRATT R G.Inverse theory applied to multi-source crosshole tomography,part 2:elastic wave-equation method[J].Geophysical Prospecting,1992,38(3):311-329

[5] PRATT R G,SHIN C,HICK G J.Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

[6] SHIN C,CHA Y H.Waveform inversion in the Laplace-Fourier domain[J].Geophysical Journal International,2009,177(3):1067-1079

[7] SIRGUE L,PRATT R G.Efficient waveform inversion and imaging:a strategy for selecting temporal frequencies[J].Geophysics,2004,69(24):231-248

[8] YOON K,SANG S,CAI J,et al.Improvements in time domain FWI and its applications[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[9] 张广智,孙昌路,潘新朋,等.变密度声波全波形反演中密度影响因素及反演策略[J].吉林大学学报(地球科学版),2016,46(5):1550-1560 ZHANG G Z,SUN C L,PAN X P,et al.Influence factors and strategy of inversion for density of acoustic full waveform inversion with variable density[J].Journal of Jilin University ( Earth Science Edition),2016,46(5):1550-1560

[10] PRZEBINDOWSKA A,KURZMANN A,KÖHN D,et al.The role of density in acoustic full waveform inversion of marine reflection seismics[R].Copenhagen:74thEAGE Annual Conference & Exhibition,2012

[11] FORGUES E,LAMBARÉ G.Parameterization study for acoustic and elastic ray plus born inversion[J].Journal of Seismic Exploration,1997,6(2/3):253-277

[12] JEONG W,LEE H Y,MIN D J.Full waveform inversion strategy for density in the frequency domain[J].Geophysical Journal International,2012,188(3):1221-1242

[13] PRIEUX V,BROSSIER R,OPERTO S,et al.Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the valhallfield,part 1:imaging compressional wave speed,density and attenuation[J].Geophysical Journal International,2013,194(3):1640-1664

[14] 董良国,迟本鑫,陶纪霞,等.声波全波形反演目标函数性态[J].地球物理学报,2013,56(10):3445-3460 DONG L G,CHI B X,TAO J X,et al.Objective function behavior in acoustic full-waveform inversion[J].Chinese Journal of Geophysics,2013,56(10):3445-3460

[15] RAKNES E B,ARNTSEN B.Strategies for elastic full waveform inversion[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:1222-1226

[16] 杨积忠,刘玉柱,董良国.变密度声波方程多参数全波形反演策略[J].地球物理学报,2014,57(2):628-643 YANG J Z,LIU Y Z,DONG L G.A multi-parameter full waveform inversion strategy for acoustic media with variable density[J].Chinese Journal of Geophysics,2014,57(2):628-643

[17] 杨勤勇,胡光辉,王立歆.全波形反演研究现状及发展趋势[J].石油物探,2014,53(1):77-83 YANG Q Y,HU G H,WANG L X.Research status and development trend of full waveform inversion[J].Geophysical Prospecting for Petroleum,2014,53(1):77-83

[18] 张文祥,吴国忱,王玉梅.基于特征能量梯度算子的时间域全波形反演[J].地球物理学进展,2015,30(1):363-371 ZHANG W X,WU G C,WANG Y M.Time domain full waveform inversion based on eigen energy gradient equation[J].Progress in Geophysics,2015,30(1):363-371

[19] 梁展源,吴国忱.基于复频率衰减反传算子的全波形反演梯度构建[J].地震学报,2016,38(2):232-243,328 LIANG Z Y,WU G C.Full waveform inversion gradient construction based on the complex frequency attenuation back propagation operator[J].Acta Seismologica Sinica,2016,38(2):232-243,328

[20] LIU D C,NOCEDAL J.On the limited memory BFGS method for large scale optimization[J].Mathematical Programming,1989,45(1/3):503-528

[21] MARTIN G S,WILEY R,MARFURT K J.Marmousi2:an elastic upgrade for Marmousi[J].The Leading Edge,2006,25(2):156-166

[22] 任浩然,黄光辉,王华忠,等.地震反演成像中的Hessian算子研究[J].地球物理学报,2013,56(7):2429-2436 REN H R,HUANG G H,WANG H Z,et al.A research on the Hessian operator in seismic inversion imaging[J].Chinese Journal of Geophysics,2013,56(7):2429-2436

[23] 陈永芮,李振春,秦宁,等.波动方程双向照明优化的全波形反演[J].地球物理学进展,2013,28(6):3015-3021 CHEN Y R,LI Z C,QIN N,et al.Full waveform inversion with wave equation bi-directional illumination optimization[J].Progress in Geophysics,2013,28(6):3015-3021

[24] 李万万.基于波动方程正演的地震观测系统设计[J].石油地球物理勘探,2008,43(2):134-141 LI W W.Design of seismic geometry based on wave equation forward simulation[J].Oil Geophysical Prospecting,2008,43(2):134-141

(编辑：戴春秋)

The stepped multi-parameter FWI of acoustic media in time-domain by L-BFGS method with illumination analysis

ZHANG Guangzhi1,2,JIANG Lanjie1,2,SUN Changlu3,HUANG Yishuang4

(1.ChinaUniversityofPetroleum,Qingdao266580,China;2.FunctionLaboratoryofMarineGeo-ResourceEvaluationandExplorationTechnology,Qingdao266235,China;3.ChinaTianchenEngineeringCorporation,Tianjin300400,China;4.SinopecExplorationCompany,Chengdu610041,China)

Density can predict fluid saturation of reservoir and plays an important role in reservoir interpretation and hydrocarbon prediction.Due to the cross-talk effects of velocity and density,density is difficult to reconstruct in multi-parameter full waveform inversion.To solve this problem,the strategy of stepped multi-parameter FWI is chosen which takes more accurate velocity inversion result as the initial velocity model and combines with the initial density model for next step simultaneous multi-parameter inversion.The more accurate the initial velocity model is,the more accuracy the density result is.In the meantime,to further balance the energy of gradient and reduce the influence of geometric diffusion in seismic wave propagation,the preprocessing L-BFGS method is carried out based on the illumination analysis.The numerical examples testify the feasibility of the stepped multi-parameter FWI and the pretreatment L-BFGS method that they can improve accuracy of result.

multi-parameter FWI,stepped inversion,L-BFGS,illumination analysis,acoustic media

2016-10-17;改回日期：2016-11-25。

张广智(1971—),男,教授,主要从事岩石物理、储层预测、流体识别方面的研究工作。

国家自然科学基金(41674130)、国家科技重大专项(2016ZX05027004-001,2016ZX05002-005-009)、国家重点基础研究发展计划(973计划)项目(2014CB239201-7HZ)、中央高校基本科研业务费专项资金(15CX08002A)联合资助。This research is financially supported by the National Natural Science Foundation of China (Grant No.41674130),the National Science and Technology Major Project of China (Grant Nos.2016ZX05027004-001,2016ZX05002005-009) (973 Program),the National Key Basic Research and Development Program of China (Grant No.2014CB239201-7HZ) and the Fundamental Research Funds for the Central Universities (Grant No.15CX08002A).

P631

A

1000-1441(2017)01-0031-07

10.3969/j.issn.1000-1441.2017.01.004