全波形反演中波場(chǎng)存儲(chǔ)量的優(yōu)化分析

全波形反演中波場(chǎng)存儲(chǔ)量的優(yōu)化分析

0波場(chǎng)重建策略地球物理勘探的主要目的是獲得地球上的精確圖像。近年來(lái),隨著計(jì)算硬件性能的提升以及GPU(GraphicProcessingUnit)等高性能計(jì)算設(shè)備在地震數(shù)據(jù)成像中的成功應(yīng)用,具有高分辨能力的全波形反演(FullWaveformInversion,FWI)方法受到了地球物理學(xué)家們青睞時(shí)間域的FWI方法最先由Tarantola實(shí)現(xiàn),隨后,Pratt將FWI發(fā)展到頻率域除了初始模型的問題會(huì)干擾FWI的應(yīng)用效果,巨大的計(jì)算量和存儲(chǔ)量同樣也成為影響FWI實(shí)用化推廣的棘手問題。在FWI迭代計(jì)算過(guò)程中,反演梯度的計(jì)算與逆時(shí)偏移(ReverseTimeMigration,RTM)方法流程類似,因此,采用互相關(guān)成像方法時(shí),震源波場(chǎng)的存儲(chǔ)量決定著RTM和FWI方法對(duì)計(jì)算硬件存儲(chǔ)性能的要求。同理,針對(duì)RTM方法的存儲(chǔ)策略同樣可借鑒應(yīng)用至FWI方法中。Clapp提出在速度外圍構(gòu)建隨機(jī)速度層作為邊界條件在前人研究工作的基礎(chǔ)上,本文提出了一種基于波場(chǎng)重建的時(shí)域FWI方法。借助拉格朗日乘子法,首先推導(dǎo)了FWI的梯度表達(dá)式,并引入解析步長(zhǎng)方法和混合共軛梯度優(yōu)化方法求解該反問題。針對(duì)FWI成像梯度計(jì)算所需的巨大存儲(chǔ)量問題,我們引入了一種能大幅降低存儲(chǔ)成本的波場(chǎng)重建策略,為時(shí)域FWI的實(shí)用化推廣提供一套切實(shí)可行的理論方案。最后通過(guò)勻速模型和Marmousi2模型的波場(chǎng)重建和FWI測(cè)試驗(yàn)證了本文方法的有效性。1方法原理1.1震源波場(chǎng)模擬時(shí)空域二階聲波常密度波動(dòng)方程為:式中,x是空間坐標(biāo);x式中,x全波形反演(FWI)通過(guò)最小化模擬數(shù)據(jù)和野外觀測(cè)數(shù)據(jù)之間的誤差,推導(dǎo)地下的速度參數(shù),其目標(biāo)函數(shù)通常定義為模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù)殘差的L2范數(shù)式中,ns和ng分別為炮點(diǎn)和檢波點(diǎn)數(shù)目;p該目標(biāo)函數(shù)受到模擬數(shù)據(jù)滿足的波動(dòng)方程的限制,我們采用拉格朗日乘子法法求解該約束性優(yōu)化問題,建立拉格朗日乘子符為:式中,λ(x,t;x該伴隨方程與聲波方程(1)具有相同的表達(dá)形式,因此只需要將震源項(xiàng)替換為模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù)的殘差,即可利用相同的數(shù)值求解方法模擬震源波場(chǎng)和伴隨波場(chǎng)。因此,傳統(tǒng)的聲波全波形反演的流程可以簡(jiǎn)述為:1)沿時(shí)間正向求解方程(1),得到正向傳播的震源波場(chǎng)p2)沿時(shí)間反方向求解方程(6),得到反向傳播的檢波點(diǎn)波場(chǎng)λ(x,t;x3)從計(jì)算機(jī)內(nèi)存中提取保存的波場(chǎng),應(yīng)用公式(5),構(gòu)建FWI的梯度。4)采用合適的反演算法,計(jì)算優(yōu)化的步長(zhǎng),更新速度參數(shù)v。5)重復(fù)1~4步驟,直到目標(biāo)函數(shù)值小于設(shè)定的值或滿足最大迭代次數(shù)。1.2快速搜索步長(zhǎng)方法上一節(jié),我們已經(jīng)推導(dǎo)了FWI的梯度,給出了傳統(tǒng)FWI的實(shí)現(xiàn)流程。然而,為了得到合理的反演結(jié)果,最優(yōu)化的搜索步長(zhǎng)也尤其重要,常規(guī)的步長(zhǎng)求取方法有很多,如:固定步長(zhǎng)將第k+1次迭代的目標(biāo)函數(shù)進(jìn)行二階泰勒展開,可得:式中,α式中,<,>表示兩個(gè)矢量之間的內(nèi)積;我們近似了Hessian矩陣H式中,f為式(1)所示的正演算子,ε為試探步長(zhǎng),它滿足條件1.3嘴唇交叉試驗(yàn)常規(guī)最速下降方法在迭代尋優(yōu)的過(guò)程,收斂速度非常緩慢,會(huì)帶來(lái)鋸齒折疊效應(yīng)。為了在迭代過(guò)程中,求解一個(gè)合理的下降方向d其中,β在得到共軛梯度方向和優(yōu)化步長(zhǎng)后,速度參數(shù)則能夠迭代地更新,更新公式可以表示為:2逆時(shí)波場(chǎng)重建在傳統(tǒng)FWI的實(shí)現(xiàn)過(guò)程中,需要對(duì)背景震源波場(chǎng)和伴隨波場(chǎng)進(jìn)行零延遲互相關(guān)來(lái)構(gòu)建梯度。但二者在時(shí)間方向上互逆,因此可預(yù)先保存某一方向所有時(shí)刻的波場(chǎng)值,當(dāng)另一個(gè)波場(chǎng)傳播時(shí),提取相應(yīng)時(shí)刻保存的波場(chǎng)進(jìn)行互相關(guān)運(yùn)算,但由此產(chǎn)生巨大的存儲(chǔ)需求和I/O(Input/Output)耗時(shí)通常是難以接受的。因此,本文提出一種適用于FWI的逆時(shí)重建震源波場(chǎng)的方法,旨在降低了FWI方法的存儲(chǔ)需求和I/O耗時(shí)。本文采用高階有限差分法求解波動(dòng)方程(1):波動(dòng)方程差分公式如圖1所示,假設(shè)真實(shí)波場(chǎng)范圍為圖1中的A區(qū)域,如果采用空間12階差分計(jì)算格式求解波動(dòng)方程,僅需保存A區(qū)域外側(cè)6層波場(chǎng)數(shù)據(jù),如圖1中B區(qū)域,即可滿足差分條件,C區(qū)域?yàn)槲者吔鐓^(qū)域。為了實(shí)現(xiàn)震源波場(chǎng)的逆時(shí)重建,震源波場(chǎng)沿正向傳播時(shí),我們需要記錄每時(shí)刻的B區(qū)域的波場(chǎng)信息以及最后時(shí)刻的A區(qū)域的波場(chǎng)信息。利用這些保存的信息,我們可以實(shí)現(xiàn)逆時(shí)波場(chǎng)重建:首先通過(guò)波動(dòng)方程的逆向運(yùn)算實(shí)現(xiàn)震源波場(chǎng)的逆時(shí)傳播,同時(shí),將對(duì)應(yīng)時(shí)刻的B區(qū)域內(nèi)的波場(chǎng)使用正傳過(guò)程中保存的波場(chǎng)更新,以滿足差分運(yùn)算條件。在震源波場(chǎng)逆時(shí)重建過(guò)程中,結(jié)合公式5即可進(jìn)行梯度的構(gòu)建,直到所有時(shí)刻和所有炮完成疊加。其余反演步驟與傳統(tǒng)方法一致,綜上所述,本文提出的基于波場(chǎng)重建的時(shí)域FWI流程如圖2所示。文中首先采用簡(jiǎn)單的均勻各向同性的介質(zhì)進(jìn)行震源波場(chǎng)逆時(shí)重建方法測(cè)試,該模型的水平和垂直方向網(wǎng)格點(diǎn)均為301,沿水平和垂直方向的網(wǎng)格尺寸均為10m。該模型的縱波傳播速度為2800m/s,為了吸收來(lái)自由于邊界截?cái)喈a(chǎn)生的反射波,我們應(yīng)用了完美匹配層吸收邊界(PML)圖3為均勻各向同性介質(zhì)不同時(shí)刻的波場(chǎng)快照,其中圖3a和圖3b分別顯示了正向傳播和本文方法逆時(shí)重建的震源波場(chǎng),從左自右的時(shí)刻分別為0.15,0.30,0.45,0.60s。我們不難發(fā)現(xiàn),模型區(qū)域內(nèi)的震源波場(chǎng)已經(jīng)被很好地重建出來(lái),由于PML邊界內(nèi)的波場(chǎng)不參與梯度的計(jì)算,因此并沒有重建。我們進(jìn)一步對(duì)波場(chǎng)快照抽取單道進(jìn)行對(duì)比,正傳波場(chǎng)和重建波場(chǎng)的單道數(shù)據(jù)對(duì)比如圖4所示,抽取位置在t=0.45s,x=2km處。觀察可知,正傳的單道圖4a和逆時(shí)重建的單道圖4b幾乎沒有差別。兩者的誤差單道圖4c所示,可知單道誤差10為了證明文章中所述方法能夠節(jié)約存儲(chǔ),我們對(duì)比了傳統(tǒng)全波場(chǎng)存儲(chǔ)的方法和文章所提出方法。假設(shè),該模型的傳播總時(shí)間間隔為5000,波場(chǎng)數(shù)據(jù)在GPU端的數(shù)據(jù)類型為4字節(jié)浮點(diǎn)型,全波場(chǎng)存儲(chǔ)的計(jì)算方法為:nx·hz·nt·4byte,本文的方法的存儲(chǔ)為z·(nx+nz)·nt·6·4byte,其中6表示為有限差分方法的空間差分階數(shù)的1/2,nx和nz為模型沿水平和垂直方向的空間網(wǎng)格點(diǎn)數(shù)。存儲(chǔ)消耗對(duì)比如表1所示,這驗(yàn)證了本文所述方法能夠大幅降低存儲(chǔ)。特別地,當(dāng)模型為3D時(shí),文中所述的方法將大幅降低存儲(chǔ)消耗,加速反演的計(jì)算速度。3模型測(cè)試為了驗(yàn)證文中所述波場(chǎng)重建方法和全波形反演方法的有效性,分別對(duì)Marmousi2模型進(jìn)行的波場(chǎng)重建測(cè)試和FWI方法測(cè)試。3.1震源波場(chǎng)重建首先,進(jìn)行了基于Marmousi2模型的波場(chǎng)重建測(cè)試,速度模型如圖5所示。模型沿水平和垂直方向的網(wǎng)格點(diǎn)數(shù)分別為500和200,沿水平和垂直方向的網(wǎng)格間距均為10m,為了壓制了來(lái)自邊界的虛假反射,我們?cè)谀Ｐ退闹芗虞d了厚度為50層PML吸收邊界。震源位置位于水平方向2.5km,垂直方向0.65km處,震源函數(shù)為主頻20Hz的雷克子波,時(shí)間采樣間隔為1ms。分別截取正傳至t=0.20,0.40,0.60,0.80s時(shí)刻的震源波場(chǎng),如圖6a所示,當(dāng)?shù)卣鸩▊鞑ブ磷畲髸r(shí)刻5s時(shí),進(jìn)行波場(chǎng)逆時(shí)重建,重建所得的t=0.20,0.40,0.60,0.80s時(shí)刻的波場(chǎng)如圖6b所示。根據(jù)圖6對(duì)比可知,本文所述的逆時(shí)重建方法能夠很好地恢復(fù)出在模型區(qū)域內(nèi)的震源波場(chǎng)。進(jìn)一步對(duì)波場(chǎng)快照抽取單道數(shù)據(jù)進(jìn)行對(duì)比,所選數(shù)據(jù)為t=0.40s,x=3km處的波場(chǎng)數(shù)據(jù),對(duì)比結(jié)果如圖7所示。正傳波場(chǎng)和逆時(shí)重建波場(chǎng)的的單道數(shù)據(jù)分別如圖7a和7b所示,二者并無(wú)肉眼可見差別,對(duì)兩者做相減運(yùn)算,分析所得誤差如圖7c所示,誤差為103.2迭代次數(shù)對(duì)fwi方程力學(xué)性能的影響為了驗(yàn)證本文提出的FWI方法對(duì)復(fù)雜模型的有效性,對(duì)Marmousi2模型進(jìn)行了FWI測(cè)試,精確速度模型如圖5所示。模型的基本信息與波場(chǎng)重建測(cè)試相同。震源和檢波器均置于地面,震源的范圍在0~5km,炮間距為250m,共20炮,采用主頻為10Hz的雷克子波,共500個(gè)檢波器,檢波點(diǎn)間距為10m。時(shí)間采樣間隔為1ms,總時(shí)間長(zhǎng)度為5s。由于實(shí)際地震勘探無(wú)法獲得精確的速度模型,因此在進(jìn)行FWI測(cè)試時(shí),首先將圖5所示的精確速度進(jìn)行平滑,如圖8所示,以此作為初始速度模型。不同迭代次數(shù)的反演結(jié)果如圖9所示。隨著迭代次數(shù)的增加,反演結(jié)果逐漸接近真實(shí)速度,對(duì)比不同迭代次數(shù)的反演結(jié)果可進(jìn)一步發(fā)現(xiàn),初期的反演結(jié)果對(duì)中淺部的速度構(gòu)造更新效果較好,當(dāng)?shù)螖?shù)超過(guò)50次后,深層的速度構(gòu)造也得以更新,趨近真實(shí)速度。分別提取了水平方向1,2,3,4km處的反演結(jié)果的單道數(shù)據(jù),如圖10所示,黑色實(shí)線表示初始速度,藍(lán)色實(shí)線表示真實(shí)速度,綠色和紅色實(shí)線分別表示迭代20次和100次的反演結(jié)果。曲線對(duì)比的結(jié)果與前文所述結(jié)論一致,執(zhí)行20次迭代計(jì)算后,文中所述的方法可以很好地反演出中淺層的速度,隨著迭代次數(shù)的增加,深層構(gòu)造也得到了有效更新。圖11顯示了目標(biāo)函數(shù)隨迭代次數(shù)的下降情況,隨著迭代計(jì)算,理論模擬數(shù)據(jù)與觀測(cè)數(shù)據(jù)誤差逐漸減少,進(jìn)一步驗(yàn)證了本文方法的有效性。4重建策略與模型試算波場(chǎng)數(shù)據(jù)的存儲(chǔ)策略對(duì)全波形反演方法的應(yīng)