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

成景旺 毛寧波* 呂曉春 常鎖亮 嚴(yán) 皓 仲 華

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

1 引言

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

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

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

2.1 CFS-NPML邊界條件

采用時(shí)間域一階應(yīng)力—速度方程進(jìn)行全波形反演。在二維情況下，彈性介質(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ù)中的兩個(gè)參數(shù)，滿足κ≥1和α≥0。當(dāng)κ=1且α=0時(shí)就變?yōu)槌Ｒ?guī)的PML邊界條件。參數(shù)κ主要用于吸收廣角入射時(shí)內(nèi)邊界產(chǎn)生的瞬逝波，而參數(shù)α主要影響對波的低頻成分的吸收。σp為伸展坐標(biāo)系下的衰減因子。本文的參數(shù)計(jì)算公式[24]為

(4)

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

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

(5)

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

(6)

(7)

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

(8)

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

(9)

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

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

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

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

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

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

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

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

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

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

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

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

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

3 全波形反演模型試算

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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