張廣智,姜嵐杰,孫昌路,黃義雙

(1.中國石油大學(華東),山東青島266580;2.海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島266235;3.中國天辰工程有限公司,天津300400;4.中國石油化工股份有限公司勘探分公司,四川成都610041)

基于照明預處理的分步多參數(shù)時間域聲波全波形反演方法研究

張廣智1,2,姜嵐杰1,2,孫昌路3,黃義雙4

(1.中國石油大學(華東),山東青島266580;2.海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島266235;3.中國天辰工程有限公司,天津300400;4.中國石油化工股份有限公司勘探分公司,四川成都610041)

密度是地震勘探中最重要的信息之一,在巖性解釋、儲層流體預測等多個方面起到不可替代的作用。但在全波形反演(FWI)中,由于密度和速度串擾的影響,很難反演出理想的密度信息。針對該問題,采取分步多參數(shù)全波形反演的策略,將反演的高精度速度結果作為初始速度模型,聯(lián)合初始密度模型進行下一步多參數(shù)同時反演,通過提高初始速度模型的精度,得到更為準確的密度結果。同時,為了進一步平衡梯度能量,減小地震波傳播過程中幾何擴散的影響,采用照明預處理L-BFGS法提高反演精度。模型測試結果表明,分步多參數(shù)全波形反演及照明預處理L-BFGS法能很好地提高反演精度。

多參數(shù)全波形反演;分步反演;L-BFGS法;照明分析;聲波介質

隨著油氣勘探開發(fā)復雜程度的不斷加深,如何提高儲層預測的精度逐漸成為專家學者研究的重點,而尋找一種適用于復雜介質的反演方法尤其重要。全波形反演(FWI)利用疊前地震波場中的運動學及動力學信息反演得到地下巖石物理特性,能夠解釋復雜地質背景下的構造細節(jié)信息[1-2],自1984年[3]問世后得到了飛速發(fā)展。PRATT等將全波形反演的思想引入頻率域,推導了頻率域內最速下降法全波形反演梯度公式[4-5]。SHIN等提出Laplace-Fourier域全波形反演方法,通過引入復頻率,在缺少低頻的情況下反演出較好的結果[6]。此外,SIRGUE等將全波形反演應用于混合域,得到了較好的反演結果[7]。在這些方法中,時間域全波形反演可對數(shù)據(jù)進行靈活預處理,適用于多節(jié)點并行運算,已被成功地應用于實際地震數(shù)據(jù)處理[8]。

在單參數(shù)全波形反演中,常假設密度隨空間位置變化緩慢,從而將密度設為一已知恒定值,反演得到速度結果。但在實際地層中,密度是指示儲層流體的重要信息,在油氣解釋中起著重要作用。密度的變化會引起反射率等變化,影響地震波傳播及反演結果[9],如何從地震數(shù)據(jù)中得到有效的密度信息越來越成為精細化地震勘探的研究重點。多參數(shù)全波形反演是一種強非線性問題,且速度與密度之間存在串擾現(xiàn)象[10],嚴重影響反演結果的精度,因此,選擇合適的反演方法和策略是重建高精度參數(shù)模型的關鍵。FORGUES等通過輻射模式分析指出,很難從全波形反演中得到密度[11]。JEONG等提出頻率域密度反演策略,先將密度設置為常數(shù),反演得到彈性參數(shù),再進行速度、密度聯(lián)合反演,提高了反演精度[12]。PRIEUX等首先采用長偏移距數(shù)據(jù)反演速度,然后利用全偏移數(shù)據(jù)反演速度和密度,提高了密度反演的穩(wěn)定性,但是密度的反演結果仍然不是很理想[13]。董良國等通過分析目標函數(shù)隨密度和速度擾動的變化規(guī)律,指出全波形反演對密度的較小尺度攝動具有較好敏感度,但是無法反演中長波長的密度攝動[14]。RAKNES等討論了彈性波全波形反演過程中反演參數(shù)的組合方式,并提出多參數(shù)全波形反演策略——逐一反演各個參數(shù)[15]。楊積忠等提出分步變密度聲波方程多參數(shù)反演策略,得到了比較好的速度和密度反演結果[16]。本文將該方法運用于時間域聲波全波形反演中，利用照明預處理L-BFGS方法使不同深度的梯度能量分布更為均衡，進一步提高模型參數(shù)收斂速度和反演精度。

1 照明預處理L-BFGS

全波形反演是Bayes估計理論在勘探地球物理領域的一個應用范例[17],它將波動方程正演模擬得到的地震記錄與實際觀測地震記錄進行匹配,通過使數(shù)據(jù)之差最小建立目標函數(shù),尋找最佳模型參數(shù)。時間域全波形反演算法在實現(xiàn)過程中主要包括兩個過程:以地震子波為震源的正演過程和以地震記錄殘差為震源的波場反傳過程[18-19]。其目標函數(shù)可以寫為:

(1)

式中:E為誤差泛函;δd為波場殘差。模型參數(shù)m可以通過下式進行迭代更新:

(2)

式中:αk是迭代步長;k是迭代次數(shù);hk是迭代方向。全波形反演的目標函數(shù)高度非線性,在迭代過程中存在強烈的非線性問題,因此要尋找合適的優(yōu)化方法提高迭代的穩(wěn)定性。

牛頓類優(yōu)化算法中通過對二階偏導海森算子的計算,均衡梯度算子修正量,可以加速迭代收斂效率。但是海森矩陣占用大量的存儲空間和計算時間,對計算機的存儲能力和計算能力都提出了很高的要求,導致大規(guī)模全波形反演無法開展。1989年,LIU等提出的L-BFGS方法[20]是對牛頓類優(yōu)化方法的改進。該方法對海森矩陣進行近似代替,不直接保存海森矩陣,而是在迭代過程中保存前n次迭代的參數(shù)修正值及梯度信息,大大減小了存儲量及計算量[18,20-22]:

(3)

利用近似海森矩陣對梯度的預處理可以校正波場幾何擴散對不同深度梯度的影響,增加參數(shù)在模型深部的校正量,同時降低淺層不合理的低頻修正量,提高建模、成像的分辨率。但是,在L-BFGS法前n次迭代過程中,不能計算得到近似海森矩陣,梯度的能量不平衡,影響反演結果。全波形反演通過正傳波場和反傳波場在時間上的零延遲互相關來計算模型參數(shù)的更新量,因此梯度能量的不均衡很大程度上是由地震波傳播過程中幾何擴散的影響引起的。而照明分析可以研究地震波在地下介質中傳播的能量分布情況,進行能量補償[23-24],同時,全波形反演需要大量的計算量和存儲空間,為了在得到較好反演效果的同時盡量減少全波形反演的計算負擔,我們采用單向照明強度近似代替雙向照明強度的方法,通過正演模擬得到模型各個點的波場,將其作為該點單向照明強度。但檢波器接收到的地震波是雙程旅行的過程[23],為了便于計算,我們將模型中各個點波場值的平方作為雙向照明強度,對前n次L-BFGS法計算的梯度值進行預處理。則N炮照明總強度IS可用下式表示:

(4)

式中:ui(x,z)為t時刻在模型坐標(x,z)處的波場值。因此,參數(shù)m的前n次迭代可以用下式來表示:

(5)

2 多參數(shù)全波形反演

時間域二階變密度聲波方程如下:

(6)

式中:κ表示體積模量;ρ表示密度;p表示壓力場;xs為震源坐標;s為震源函數(shù)。則速度和密度的梯度可分別用下式表示:

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

多參數(shù)全波形反演因速度和密度之間的串擾影響非常不穩(wěn)定,難以獲得精確的反演結果[10]。速度的變化會引起地震波傳播的旅行時和振幅發(fā)生擾動,而密度的變化會影響地震波的振幅大小,且隨著角度的變大,波場變化逐漸減小[9]。同時,地震反射波對密度的較小尺度攝動非常敏感,而中長波長的密度攝動對反射波基本沒有影響[16]。因此,多參數(shù)全波形反演雖然可以重建較好的速度模型,但只能重建密度

的高頻成分[9,14,16]。與此同時,全波形反演,尤其是多參數(shù)全波形反演,對初始模型的依賴性較強,初始速度模型越準確,密度反演結果越精確。

基于以上分析,我們將楊積忠等提出的分步反演策略[16]用于多參數(shù)時間域全波形反演,得到速度和密度反演結果。第一步是利用多參數(shù)全波形反演得到較為精確的速度模型,此時密度反演結果較差;第二步是將第一步反演的速度結果作為初始速度模型,結合初始密度模型進行速度、密度聯(lián)合反演,此時由于初始速度模型精度的提高,密度反演的精度也有很大提高。為了進一步提高反演精度,可將第二步得到的速度反演結果作為初始模型,利用初始密度模型再進行一輪速度、密度聯(lián)合反演。

3 模型測試

3.1 照明預處理L-BFGS方法測試

為了測試照明預處理L-BFGS方法的有效性,我們將其應用于Marmousi模型，進行速度參數(shù)的全波形反演,并與最速下降法和L-BFGS法進行對比。實際速度模型如圖1a所示,網格點數(shù)為260×130,網格間距10m。地表放置38個震源,震源間距70m,并放置260個檢波器,檢波器間距10m。所用震源為主頻30Hz的雷克子波。圖1b為初始速度模型。

圖2是最速下降法、L-BFGS法、照明預處理L-BFGS法速度反演結果,可以看到:最速下降法反演效果不理想,雖然淺層速度反演得很好,但是深層速度值偏差較大;L-BFGS方法在反演出較好的淺層速度的基礎上,由于自身可以校正波場幾何擴散的影響,深部速度的校正量增加,速度反演更為準確;而經過照明預處理后的L-BFGS方法進一步彌補了波場傳播時深層能量的缺失,得到了更好的深部速度反演結果。圖3是三種方法在不同水平位置處的反演速度隨深度變化的曲線,可以看到照明預處理L-BFGS法得到的速度結果與真實值擬合最好,反演精度最高。圖4是三種反演方法的目標函數(shù)收斂曲線,可以看到照明預處理L-BFGS法收斂最快,失配函數(shù)最小,說明照明預處理L-BFGS方法可以提高收斂速度及精度,減小運算量及運算時間。

圖1 實際速度模型(a)和初始速度模型(b)

圖2 三種方法速度反演結果比較a 最速下降法; b L-BFGS法; c 照明預處理L-BFGS法

圖3 水平方向900m(a)、1200m(b)、1700m(c)、2250m(d)處不同方法反演的速度隨深度變化曲線

3.2 分步多參數(shù)全波形反演測試

比較最速下降法、傳統(tǒng)L-BFGS法及照明預處理L-BFGS法的效果,我們選擇反演精度最高的照明預處理L-BFGS方法,采用與單參數(shù)反演模型相同的參數(shù),進行速度、密度的聯(lián)合反演測試。圖5a是實際速度模型,圖5b是用Gardner公式和實際速度模型計算得到的實際密度模型;圖6a和圖6b是初始速度模型和初始密度模型,圖7a和圖7b是速度和密度第一步反演結果。從圖7可以看出,速度和密度的反演精度都不高,淺層速度反演較好,構造也得到很好的重建,但深部速度和實際值仍有很大偏差,密度的反演結果精度不高,不能滿足反演精度的要求。

為了得到更好的反演結果,我們舍棄密度反演結果,將第一步的速度反演結果(圖7a)作為初始速度模型,聯(lián)合初始密度模型再進行一次多參數(shù)同時反演(第二步反演),結果如圖8所示?？梢钥吹?速度反演結果精度有了明顯提高,中深層速度值與真實值吻合較好,密度反演精度較之前也有了提高,證明分步多參數(shù)全波形反演在提高反演精度方面是有效的。但是,受到初始速度模型精度的制約,密度精度的提高有限(圖8b)。要進一步重建更為精確的密度模型,需要更為準確的初始速度模型,因此我們將第二步的速度反演結果作為初始速度模型,結合初始密度模型,又進行了一次多參數(shù)聯(lián)合反演(最終反演),結果如圖9所示。由圖9可以看到，速度和密度反演的精度都有了明顯改善。比較分步多參數(shù)全波形反演過程中不同水平位置處速度、密度值縱向剖面圖(圖10), 可以看到速度、密度聯(lián)合反演的精度每一次都在提高,最終的模型反演結果與實際模型值擬合較好。

圖4 不同反演方法的目標函數(shù)收斂曲線對比

圖5 實際速度模型(a)和實際密度模型(b)

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

圖7 照明預處理L-BFGS分步多參數(shù)全波形反演第一步反演結果a 速度; b 密度

圖8 照明預處理L-BFGS分步多參數(shù)全波形反演第二步反演結果a 速度; b 密度

圖9 多參數(shù)全波形反演最終反演結果a 速度; b 密度

圖10 分步多參數(shù)全波形反演速度(a,b)、密度(c,d)隨深度變化曲線

4 結論

在多參數(shù)全波形反演中,由于速度和密度的耦合關系影響了反演的精度和效率,密度反演結果很不理想,選擇合適的反演算法和反演策略對全波形反演至關重要。將分步多參數(shù)全波形反演策略應用于時間域多參數(shù)聲波全波形反演中,通過建立更為準確的速度初始模型,減小速度和密度間的串擾影響,得到了更為準確的密度反演結果。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] 王西文,楊午陽,呂彬,等.把握世界物探技術發(fā)展現(xiàn)狀,促進地震資料處理、解釋技術進步[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] 董良國,遲本鑫,陶紀霞,等.聲波全波形反演目標函數(shù)性態(tài)[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] 楊積忠,劉玉柱,董良國.變密度聲波方程多參數(shù)全波形反演策略[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] 楊勤勇,胡光輝,王立歆.全波形反演研究現(xiàn)狀及發(fā)展趨勢[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] 陳永芮,李振春,秦寧,等.波動方程雙向照明優(yōu)化的全波形反演[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] 李萬萬.基于波動方程正演的地震觀測系統(tǒng)設計[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)、國家重點基礎研究發(fā)展計劃(973計劃)項目(2014CB239201-7HZ)、中央高?；究蒲袠I(yè)務費專項資金(15CX08002A)聯(lián)合資助。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