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

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

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