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

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

1 引言

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

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

2 基于邊界存儲(chǔ)的波場(chǎng)重建方法

2.1 CFS-NPML邊界條件

采用時(shí)間域一階應(yīng)力—速度方程進(jìn)行全波形反演。在二維情況下，彈性介質(zhì)一階應(yīng)力—速度波動(dò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ù)α主要影響對(duì)波的低頻成分的吸收。σ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)過(guò)程中引入的輔助變量，其對(duì)應(yīng)的控制方程為

(6)

(7)

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

(8)

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

(9)

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

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

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

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

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

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

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

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

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

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

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

3 全波形反演模型試算

為驗(yàn)證本文方法的正確性，采用Marmousi模型進(jìn)行試算?；谝浑A應(yīng)力—速度彈性波波動(dòng)方程，采用預(yù)條件共軛梯度法進(jìn)行FWI[27]。為了減少反演時(shí)間，將Marmousi模型等比例縮小，并在模型上面加入一海水層，最終模型網(wǎng)格數(shù)為300×116。反演初始模型(圖5)為真實(shí)模型的二維高斯函數(shù)光滑結(jié)果，光滑過(guò)程中縱橫向的相關(guān)長(zhǎng)度均取100m?？臻g網(wǎng)格間距為10m，時(shí)間步長(zhǎng)為0.5ms，記錄長(zhǎng)度為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。反演過(guò)程中每個(gè)頻帶設(shè)置的最低迭代次數(shù)為20次。 正演過(guò)程采用基于區(qū)域分解的并行計(jì)算[29]，本次數(shù)值測(cè)試一共采用10個(gè)進(jìn)程進(jìn)行計(jì)算。為了減弱多參數(shù)彈性波FWI中各個(gè)參數(shù)之間的相互耦合性，反演過(guò)程采用楊積忠等[30]提出的多參數(shù)反演策略。最終四個(gè)頻帶的縱、橫波速度和密度反演結(jié)果如圖6～圖8所示，可以看出隨著反演截止頻率的增加，模型的細(xì)節(jié)構(gòu)造越來(lái)越明顯。為了進(jìn)一步更加精確地對(duì)比反演結(jié)果，提取橫向750m、1250m和2000m處的速度進(jìn)行對(duì)比(圖9～圖11)，可以看出，反演得到的縱、橫波速度、密度與真實(shí)相值吻合。

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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