非分裂完全匹配層邊界存儲時間域全波形反演

成景旺 毛寧波* 呂曉春 常鎖亮 嚴 皓 仲 華

(①油氣資源與勘探技術(shù)教育部重點實驗室(長江大學(xué))，湖北武漢 430100； ②華北水利水電大學(xué)資源與環(huán)境學(xué)院，河南鄭州 450011； ③太原理工大學(xué)礦業(yè)工程學(xué)院，山西太原 030024； ④中海石油(中國)天津分公司，天津 300452)

1 引言

全波形反演集地震子波估計、初始模型建立、正演模擬、反演于一體，是一套完整的理論體系，已被證明是一種建立高精度速度模型的有效方法[1]，可在時間域、頻率域、Laplace域?qū)崿F(xiàn)[2]。Tarantola[3]利用伴隨狀態(tài)法高效率地求取梯度，實現(xiàn)了二維時間域全波形反演；Pratt[4]將全波形反演理論推廣到頻率域；Shin等[5]針對地震數(shù)據(jù)帶寬有限、初始模型獲取困難等問題，提出利用阻尼波場零頻分量反演低頻模型作為頻率域波形反演的初始模型，即Laplace域全波形反演。頻率域全波形反演由于其固有的多尺度特征使其理論研究和實際應(yīng)用得到快速發(fā)展[6-10]。然而頻域率反演最大的問題就是大型稀疏方程組的存儲以及求解過程中巨大的內(nèi)存需求，大多借用MUMPS線性方程組求解軟件包進行求解。但對于大規(guī)模尤其是三維情況下，頻率域正演對內(nèi)存的超大要求，一般的計算機集群仍無法滿足要求，因此近年來時間域全波形反演研究成為熱點[11，12]，尤其是時間域正演聯(lián)合頻率域反演的混合算法[13-15]。而制約時間域全波形反演的關(guān)鍵問題之一就是正演波場的存儲或重建。由于全波形反演的一次梯度求取過程與逆時偏移過程相同，故可將逆時偏移的邊界存儲策略應(yīng)用于全波形反演。國內(nèi)外學(xué)者在解決逆時偏移巨大內(nèi)存需求上已做了大量的研究，Symes[16]提出了采用設(shè)置檢查點的方法以降低逆時偏移的存儲量。Dussaud等[17]指出檢查點技術(shù)雖然對存儲量的需求最小，但是卻明顯地增加了計算量。Clapp[18]提出了邊界存儲策略，該策略要求存儲邊界網(wǎng)格層內(nèi)所有時刻的波場和整個空間最后時刻的波場，反傳時作為邊界條件和初值條件，只需額外正演一次即可重建正演波場；Clapp[19]又提出了隨機邊界方法，只保存最后一個時刻所有空間點的波場值作為反傳時的初值，存儲量進一步減少但是在偏移剖面中會引入噪聲。在所有這些策略中，Clapp的邊界儲存法具有存儲量小、計算量相對較少且適用于任何邊界條件等優(yōu)點[20]，因此該邊界存儲法被廣泛研究并被用于逆時偏移。

不同的邊界條件的吸收效果以及存儲量都會不同。完全匹配層(PML)邊界條件被認為是最好的吸收邊界條件，實現(xiàn)方法主要有全局分裂完全匹配層(SPML)邊界條件、局部SPML邊界和非分裂完全匹配層(NPML)邊界條件。其中全局SPML對邊界區(qū)域和計算區(qū)域使用同樣形式的PML波動方程，編程簡單，但是由于波場分成垂直邊界和水平邊界方向兩部分(二維為2個方向，三維為3個方向)，因此在整個計算過程中所需存儲量是常規(guī)波動方程的2倍甚至3倍；局部SPML邊界條件只有在PML吸收層內(nèi)采用PML波動方程，而在計算區(qū)域采用常規(guī)波動方程，因此只有在吸收層內(nèi)需要額外增加存儲量。但局部SPML在計算區(qū)域和邊界區(qū)域使用不同形式的波動方程，因此實現(xiàn)起來較繁瑣，除了計算區(qū)域外，需要考慮模型邊界、邊角的問題，編程較困難。目前局部SPML邊界條件已被廣泛用于基于邊界儲存的波場重建[20，21]。但是傳統(tǒng)的SPML邊界條件存在一定的缺陷，對大角度入射產(chǎn)生的掠射波吸收效果不佳，為此引入了復(fù)頻移(CFS)伸展函數(shù)改進坐標變換，該邊界條件稱為復(fù)頻移(CFS-PML)吸收邊界條件[22-24]。采用復(fù)頻移伸展函數(shù)后，PML不易采用傳統(tǒng)的分裂形式實現(xiàn)，而采用不分裂卷積算法時需要進行大量卷積計算。為此，Komatitsch等[25]采用了不分裂遞推卷積方法實現(xiàn)了CFS-NPML吸收邊界，并進行了彈性波的數(shù)值模擬；Drossaret等[26]提出了基于遞歸積分的非分裂CFS-NPML邊界條件，通過引入輔助變量實現(xiàn)遞歸積分，同樣避免了卷積計算。這兩種NPML邊界條件雖實現(xiàn)方法不同，但最終計算使用的離散公式形式相同。綜合考慮，CFS-NPML不需要對波場進行分裂，實現(xiàn)簡單，邊界處所需變量個數(shù)少，因此更加適合于基于邊界存儲的波場重建。本文分析了CFS-NPML條件的實現(xiàn)原理和吸收效果，提出基于CFS-NPML邊界條件的有效邊界存儲策略進行正演波場重建，實現(xiàn)時間域全波形反演，并通過理論模型論證其可行性。

2 基于邊界存儲的波場重建方法

2.1 CFS-NPML邊界條件

采用時間域一階應(yīng)力—速度方程進行全波形反演。在二維情況下，彈性介質(zhì)一階應(yīng)力—速度波動方程可表示為

(1)

(2)

式中sp為伸展函數(shù)。CFS-NPML與SPML邊界相比，主要的區(qū)別就是引入了復(fù)頻移伸展(CFS)函數(shù)。CFS函數(shù)定義為

(3)

其中κp(p=x，z)和αp(p=x，z)為復(fù)頻移伸展函數(shù)中的兩個參數(shù)，滿足κ≥1和α≥0。當κ=1且α=0時就變?yōu)槌Ｒ?guī)的PML邊界條件。參數(shù)κ主要用于吸收廣角入射時內(nèi)邊界產(chǎn)生的瞬逝波，而參數(shù)α主要影響對波的低頻成分的吸收。σp為伸展坐標系下的衰減因子。本文的參數(shù)計算公式[24]為

(4)

式中：L為PML邊界網(wǎng)格層數(shù)；l表示與內(nèi)邊界的網(wǎng)格點距離；Pd、Pκ、Pα為多項式衰減因子的系數(shù)，取值范圍一般為[1，4]，通常情況下取2；κmax通常介于1～20；αmax=πf0，f0是震源的主頻；cp是PML內(nèi)縱波的傳播速度；Rc為邊界的吸收系數(shù)，一般取0.0001。

引入兩個輔助變量，經(jīng)詳細推導(dǎo)可得到時間域CFS-NMPL邊界條件的一階應(yīng)力速度方程[22]為

(5)

式中Ω和Ψ為推導(dǎo)過程中引入的輔助變量，其對應(yīng)的控制方程為

(6)

(7)

將式(6)和式(7)寫成統(tǒng)一形式的一階微分方程

(8)

則輔助變量可通過下式迭代求解[24]

(9)

式中： Δt為時間步長；n為時間采樣序號。比較式(5)與式(1)可以看出，通過引入輔助變量，將波動方程右端分成了正常項和衰減項兩部分。在計算區(qū)域內(nèi)，式(5)演變?yōu)槭?1)，可以按照正常波動方程進行計算。而在邊界吸收層內(nèi)可先計算正常項，然后減去由輔助變量表示的衰減項即可。因此CFS-NPML邊界條件不需要對波場進行分裂，只需在邊界吸收層內(nèi)額外增加幾個輔助變量的存儲量即可。對于SPML邊界條件，需要對波場分量分為x方向和z方向兩部分(三維情況下為三部分)，如vx=vxx+vxz。對全局SMPL、局部SPML以及CFS-NPML在邊界吸收層內(nèi)需要申請的變量個數(shù)進行對比，如表1所示。其中全局SPML在二維情況下需要申請10個變量，三維情況下需要申請27個。但是由于全局SPML是整個模型保存變量，因此所需內(nèi)存最大。局部SPML在二維情況下邊界內(nèi)需申請15個變量，三維情況下需申請36個變量；CFS-NMPL在二維情況下邊界內(nèi)需保存13個變量，三維情況下需申請27個變量。所以CFS-NPML邊界條件不僅實現(xiàn)方便，而且在吸收層內(nèi)所需的變量個數(shù)最少。

表1 不同PML吸收邊界條件邊界層內(nèi)所需變量個數(shù)統(tǒng)計

注：i=1代表x方向；i=2代表z方向；i=3代表y方向

為了對比CFS-CPML邊界與局部SPML邊界吸收效果，建立一個網(wǎng)格數(shù)為300×300的均勻介質(zhì)理論模型，縱波速度為3300m/s，橫波速度為1905m/s，空間步長為5m，時間步長為0.5ms。震源采用主頻為20Hz的Ricker子波，位于模型中心(750m，750m)，選擇較少的邊界網(wǎng)格點數(shù)進行邊界吸收效果對比。當邊界吸收網(wǎng)格點數(shù)為10時(圖1)，局部SPML邊界條件出現(xiàn)了較弱的邊界反射，而CFS-NPML的邊界吸收依然較好，即使能量放大到10倍也看不見邊界反射，能量放大到100倍時可以看到微弱的反射，明顯好于局部SPML邊界條件。當邊界吸收網(wǎng)格點數(shù)為5時(圖2)，局部SPML和CFS-NPML均有邊界反射，但CFS-NPML的邊界反射能量較弱，即使能量放大10倍也要優(yōu)于邊界網(wǎng)格點數(shù)為10時的局部SPML條件。在(750m，300m)處取單道波形進行對比(圖3)，未做邊界吸收處理的反射波與直達波能量相當，其他四種方式的邊界條件均對邊界反射(0.35s后)有吸收作用。將該道記錄中的邊界反射局部放大，可以清楚地看到CFS-NMPL邊界條件的反射能量均小于SPML邊界條件。設(shè)E為未作邊界處理的邊界反射波最大振幅，E′為吸收邊界處理后的邊界反射波最大振幅，將比值(E-E′)/E作為分析邊界反射好壞的標準。采用局部SPML(邊界網(wǎng)格數(shù)為5)的邊界反射吸收率為88%，采用局部SPML(邊界網(wǎng)格為10)的邊界反射吸收率為92.9%； 采用CFS-PML(邊界網(wǎng)格數(shù)為5)的邊界反射吸收率為98.1%，采用CFS-PML(邊界網(wǎng)格數(shù)為10)的邊界反射吸收率為99.9%。因此為了進一步減少邊界存儲內(nèi)存，在盡量減少邊界網(wǎng)格點數(shù)的前提下，CFS-NMPL比常規(guī)PML邊界條件具有更好吸收效果。除此之外CFS-NPML邊界條件不需要對波場進行分裂，計算效率高，編程難度小，是時間域全波形反演波場重建最有效的正演邊界處理方法。

圖1 邊界網(wǎng)格點數(shù)為10時SPML(上)和CFS-NPML(下)邊界條件吸收效果對比

2.2 基于邊界存儲的波場重建實現(xiàn)原理

基于邊界儲存的波場重建就是要利用最后一個時刻(二階時間離散)的波場值通過反傳計算出前面任意時刻的波場值。無論是哪種方式的PML邊界條件，其在時間方向上都是不可逆的，因此需要存儲邊界區(qū)域每個時刻的全部波場值。反向傳播作為邊界條件替換邊界處的波場值，此時只需計算內(nèi)部區(qū)域，不需再做邊界處理就可以完全計算出前面任意時刻的波場值。

基于邊界存儲的波場重建可按如下步驟進行：

圖2 邊界網(wǎng)格點數(shù)為5時SPML和CFS-NPML邊界條件吸收效果對比

圖3 不同邊界條件(750m，300m)處的地震記錄

(7)若n>1，則返回步驟(2)；n=1結(jié)束，完成波場重建。

在上述波場重建過程中，應(yīng)特別注意應(yīng)力分量和速度分量的更新迭代次序一定要與正向傳播時的相反，否則隨著時間的迭代，在計算區(qū)域與邊界區(qū)域交界處會產(chǎn)生新的繞射源，導(dǎo)致重建波場中存在干擾波。為了驗證上述CFS-NPML邊界存儲波場重建的有效性，應(yīng)用理論模型進行測試。將Marmousi模型上面加入一海水層，并等比例縮小。最終模型網(wǎng)格數(shù)為300×116，網(wǎng)格間距為10m，時間步長為0.5ms，震源采用位于(1500m，120m)處20Hz的Ricker子波，為了盡量減少邊界保存所需內(nèi)存，取PML吸收層網(wǎng)格點數(shù)為10。圖4為正向傳播和逆向重建的200ms、450ms和800ms波場快照。圖4左邊一列為正向傳播波場，中間一列為速度分量和應(yīng)力分量的更新次序與正向傳播一樣時得到的重建波場，最右邊一列速度分量與應(yīng)力分量更新次序與正向傳播相反的重建波場?？梢钥闯?，利用上述重建步驟得到的重建波場(右列)與正向傳播波場完全一致，說明利用該邊界存儲策略完全能夠重建正演波場。而當速度分量和應(yīng)力分量更新次序與正向傳播一致時，可以看出得到的重建波場大致與正向傳播波場相同。800ms時刻的重建波場與正向傳播波場一致，但隨著波場由邊界向里邊傳播，在計算區(qū)域與邊界交界處出現(xiàn)了新的干擾波場，并隨著正常波場一起向模型內(nèi)部傳播。從200ms和450ms的波場快照圖中可明顯看出存在干擾波，這說明在波場重建過程中速度分量和應(yīng)力分量的更新次序?qū)χ亟ú▓龅闹匾浴?/p>

圖4 200(上)、450(中)、800ms(下)正向傳播波場與重建波場快照對比

3 全波形反演模型試算

為驗證本文方法的正確性，采用Marmousi模型進行試算?；谝浑A應(yīng)力—速度彈性波波動方程，采用預(yù)條件共軛梯度法進行FWI[27]。為了減少反演時間，將Marmousi模型等比例縮小，并在模型上面加入一海水層，最終模型網(wǎng)格數(shù)為300×116。反演初始模型(圖5)為真實模型的二維高斯函數(shù)光滑結(jié)果，光滑過程中縱橫向的相關(guān)長度均取100m?？臻g網(wǎng)格間距為10m，時間步長為0.5ms，記錄長度為2s。震源采用主頻為20Hz的Ricker子波，CFS-NMPL邊界網(wǎng)格數(shù)取為10。

震源點和接收點深度均設(shè)置為140m。一共激發(fā)19炮，炮間距為150m，每一炮都有固定的275道接收，道間距為10m。由于時間域FWI是全頻帶數(shù)據(jù)反演，反演容易陷入局部極小值，因此采用濾波器進行多尺度反演[28]，從低到高依次給定截止頻率實現(xiàn)不同頻帶的反演。本文采用巴特沃斯低通濾波器進行濾波，給定的低通截止頻率分別為5Hz、9Hz、16.36Hz和29.2Hz。反演過程中每個頻帶設(shè)置的最低迭代次數(shù)為20次。 正演過程采用基于區(qū)域分解的并行計算[29]，本次數(shù)值測試一共采用10個進程進行計算。為了減弱多參數(shù)彈性波FWI中各個參數(shù)之間的相互耦合性，反演過程采用楊積忠等[30]提出的多參數(shù)反演策略。最終四個頻帶的縱、橫波速度和密度反演結(jié)果如圖6～圖8所示，可以看出隨著反演截止頻率的增加，模型的細節(jié)構(gòu)造越來越明顯。為了進一步更加精確地對比反演結(jié)果，提取橫向750m、1250m和2000m處的速度進行對比(圖9～圖11)，可以看出，反演得到的縱、橫波速度、密度與真實相值吻合。

圖5 Marmousi模型(左)及反演初始模型(右)

表2給出了常規(guī)保存波場FWI和本文方法所需的內(nèi)存存儲量和梯度計算所用時間。根據(jù)縱、橫波速度和密度的反演梯度，需要保存vx、vz、τxz、τxx和τzz共5個變量，其中常規(guī)FWI需保存所有時刻所有網(wǎng)格點的波場值，而本文方法只需保存邊界處(邊界網(wǎng)格點數(shù)僅為10)的波場值即可(計算中波場采用單精度浮點數(shù))。

表2 兩種反演方法內(nèi)存和梯度計算時間對比

從表2可以看出，基于非分裂PML邊界存儲的FWI在內(nèi)存需求上明顯減少，僅僅是常規(guī)FWI內(nèi)存需求的23.9%。由于常規(guī)保存波場FWI在一次梯度求解中相當于進行兩次正演模擬，而基于邊界存儲法在一次梯度求取中相當于要進行三次正演模擬，因此基于邊界存儲的全波形反演必然會增加計算時間。但從表2中可以看出，由于采用并行計算技術(shù)，所有的梯度計算時間總和并沒有顯著增加，其計算時間比常規(guī)FWI增加了16%。本文提出的基于邊界存儲的FWI既可有效減少內(nèi)存需求，又能保證計算效率。

圖6 不同截止頻率反演的縱波速度

圖7 不同截止頻率反演的橫波速度

圖8 不同截止頻率反演的密度

圖9 不同橫向位置處的縱波速度真實模型(黑色)、初始模型(紅色)和反演結(jié)果(藍色)對比

圖11 不同橫向位置處的密度真實模型(黑色)、初始模型(紅色)和反演結(jié)果(藍色)對比

4 結(jié)論

本文將逆時偏移中的邊界存儲重建波場技術(shù)引入時間域全波形反演，并采用CFS-NPML邊界條件代替?zhèn)鹘y(tǒng)PML邊界條件解決了時間域全波形反演的巨大存儲量問題。CFS-NPML邊界條件不需要對波場進行分裂，在邊界吸收層內(nèi)所需的變量個數(shù)最少，且在邊界網(wǎng)格數(shù)較少的情況下，CFS-NPML邊界條件的吸收效果要優(yōu)于傳統(tǒng)PML邊界條件。因此基于CFS-NPML邊界條件實現(xiàn)波場重建，可最大化地減少邊界存儲內(nèi)存需求，是時間域全波形反演波場重建的有效方法。同時，應(yīng)用并行計算技術(shù)，可保證迭代過程中的梯度計算時間沒有顯著增加。

本文提出的基于邊界存儲的全波形反演實現(xiàn)方便簡單，既有效減少了內(nèi)存需求，又保證了計算效率。該邊界條件還可用于黏彈性或各向異性等復(fù)雜介質(zhì)的正演波場重建，也可直接推廣到三維。基于CFS-NPML邊界存儲的復(fù)雜介質(zhì)三維時間域全波形反演是進一步的研究方向。