韓如冰 郎 超

(①中國地質(zhì)科學(xué)院地質(zhì)研究所自然資源部深地動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100037； ②中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074)

0 引言

全波形反演(Full-waveform Inversion)以波動(dòng)方程作為數(shù)學(xué)模型模擬地震波的傳播規(guī)律，能夠更充分地利用地震數(shù)據(jù)信息，具有成像更準(zhǔn)確、分辨率更高的特點(diǎn)[1-2]。Tarantola[1]首先提出了基于廣義最小二乘的時(shí)間域全波形反演理論和方法，隨后Pratt[3-4]將其推廣到頻率域內(nèi)。相比于時(shí)間域，頻率域全波形反演具有計(jì)算過程穩(wěn)定、占用內(nèi)存小、不存在累計(jì)誤差、緩和局部極小、降低反演的非線性、易于處理衰減頻散效應(yīng)、易于并行計(jì)算等特點(diǎn)[5-10]，近些年得到廣泛發(fā)展與應(yīng)用[11-12]。

在頻率域計(jì)算中，選擇合適的數(shù)值方法離散波動(dòng)方程，處理區(qū)域邊界、求解大型線性方程組是必要的。常見的數(shù)值離散方法包括有限差分法[13-14]、有限元法[15-16]、譜元法[17-18]等。當(dāng)網(wǎng)格不夠密集時(shí)，常規(guī)的有限差分(Ordinary Finite-difference，OFD)方法容易出現(xiàn)明顯的數(shù)值頻散現(xiàn)象，盡管通過加密離散網(wǎng)格可以在一定程度上壓制頻散，但是相應(yīng)的計(jì)算量增大，所需時(shí)間和內(nèi)存也會(huì)增加[19]。

針對以上問題，Yang等[20-21]引入了近似解析離散化(Nearly Analytic Discrete，NAD)方法，并應(yīng)用于時(shí)間域地震波模擬，隨后Lang等[8]將其引入到頻率域正演模擬中，其基本思想是利用波場位移與其梯度值的聯(lián)合近似得到位移的高階偏導(dǎo)數(shù)，進(jìn)而離散波動(dòng)方程。相比于其他差分方法，NAD算子攜帶了更多的波場信息，特別是含有刻畫波場變化特征的梯度信息。前人研究表明，針對較粗的離散網(wǎng)格，NAD方法在壓制數(shù)值頻散方面具有很好的效果[8]，此外，與同階差分格式相比，NAD數(shù)值格式計(jì)算量小，易于并行計(jì)算。然而現(xiàn)有的頻率域四階NAD(NAD4)方法對于計(jì)算效率的提升還很有限，需要考慮在此基礎(chǔ)上構(gòu)造更加有效的數(shù)值方法。

為提高全波形反演的計(jì)算效率，本文構(gòu)造了頻率域八階NAD(NAD8)格式離散二維聲波方程，采用完全匹配層吸收邊界(Perfectly Matched Layer，PML)[22]條件。首先詳細(xì)介紹了NAD8格式的具體離散過程；然后使用一類不精確旋轉(zhuǎn)分塊三角預(yù)處理算子(IRBTP)加速廣義極小殘量方法(GMRES)[23]求解離散后的大型稀疏線性方程組；通過數(shù)值頻散分析以及波場模擬與六階NAD(NAD6)、八階OFD(OFD8)方法進(jìn)行對比，結(jié)果表明NAD8方法在壓制頻散和提高計(jì)算效率等方面具有優(yōu)勢；最后，基于共軛梯度法[24]，首次構(gòu)造了基于NAD8方法的頻率域全波形反演算法，對兩種典型分層介質(zhì)模型和經(jīng)典Marmousi模型進(jìn)行了反演試算，獲得了高保真、高分辨率的成像結(jié)果，進(jìn)而驗(yàn)證了所提方法的有效性和適用性。

1 方法原理

歸結(jié)起來，頻率域NAD方法的離散過程可概括為以下幾步：對頻率域波動(dòng)方程關(guān)于各個(gè)方向求偏導(dǎo)，并在計(jì)算區(qū)域的邊界處施加吸收邊界條件，進(jìn)而得到偏微分方程組，采用NAD網(wǎng)格差分模板離散其中的高階偏導(dǎo)數(shù)項(xiàng)，再按一定節(jié)點(diǎn)排序規(guī)則形成大型稀疏線性方程組。

假定地球?yàn)橥耆珡椥越橘|(zhì)時(shí)，可以用聲波方程研究地震波的傳播規(guī)律。常密度介質(zhì)下的二維頻率域聲波方程為

(1)

式中：u(x，z，ω)表示位移或者壓力場；ω=2πf表示角頻率；c為地震波波速；s表示震源項(xiàng)。

對式(1)兩端分別關(guān)于x方向與z方向求偏導(dǎo)，有

(2)

與時(shí)間域不同，在頻率域施加PML邊界條件的過程中，不需要進(jìn)行反褶積計(jì)算，實(shí)現(xiàn)更方便、準(zhǔn)確。首先引入復(fù)坐標(biāo)控制衰減，以x方向?yàn)槔?/p>

(3)

將衰減函數(shù)取為

(4)

(5)

式中

(6)

與NAD4方法[19]不同，對于NAD8格式，任一點(diǎn)(i,j)的高階偏導(dǎo)數(shù)需要由(i-2，j)、(i-1，j)、(i，j)、(i+1，j)、(i+2，j)五個(gè)點(diǎn)處的波場及其空間梯度的加權(quán)組合近似。求解微分方程中高階偏導(dǎo)數(shù)的離散格式為

(7)

(8)

(9)

(10)

128ui-1,j+1-128ui+1,j-1+128ui-1,j-1-

7ui+2,j-2+7ui-2,j-2)+

(11)

1408ui+1,j+1-1408ui-1,j+1-62ui+2,j-

2816ui+1,j+2816ui-1,j+62ui-2,j+

31ui+2,j-2+1408ui+1,j-1-1408ui-1,j-1-

(12)

1408ui+1,j+1-2816ui,j+1+1408ui-1,j+1-

1408ui+1,j-1+2816ui,j-1-1408ui-1,j-1-

31ui+2,j-2+62ui,j-2-31ui-2,j-2)+

(13)

式中h表示空間離散步長。將離散結(jié)果依次代入式(5)，并寫成統(tǒng)一線性方程形式

(14)

(15a)

(15b)

(16a)

(16b)

(17a)

(17b)

由式(14)的離散結(jié)果發(fā)現(xiàn)，NAD8差分格式共涉及周圍17個(gè)節(jié)點(diǎn)，顯然算子長度小于同階精度的其他格式[25]，具體差分模板如圖1所示。

為快速有效求解離散后的大型稀疏線性方程組，本文采用一類不精確旋轉(zhuǎn)分塊三角(IRBT)預(yù)處理算子加速廣義極小殘量方法(GMRES)。該方法是一種專門用于非對稱線性方程組求解的Krylov子空間方法，充分利用了系數(shù)矩陣T的稀疏結(jié)構(gòu)，結(jié)合不精確預(yù)處理的思想，在迭代次數(shù)并未顯著增加的前提下，可節(jié)省較多的運(yùn)算時(shí)間與存儲(chǔ)量，從而能明顯提高計(jì)算效率[23]。

圖1 NAD8差分網(wǎng)格模板

2 正演計(jì)算

從數(shù)值頻散分析和波場模擬兩方面考查NAD8方法相對于NAD6和常規(guī)八階有限差分方法的模擬效率。震源選取雷克(Ricker)子波，根據(jù)Fourier變換的性質(zhì)，其頻率域表達(dá)式[9]為

(18)

式中：A表示振幅，波場模擬時(shí)統(tǒng)一設(shè)置為1；f0是主頻。根據(jù)頻譜分析[8]，主要能量分布在0～3f0。定義網(wǎng)格頻率(Gf)為每個(gè)波長包含的離散網(wǎng)格點(diǎn)數(shù)，以度量網(wǎng)格剖分的大小[26]，其計(jì)算公式為

(19)

式中w表示波長。從式(19)可以看出，對于同一計(jì)算區(qū)域，Gf越大，每個(gè)波長內(nèi)網(wǎng)格點(diǎn)數(shù)目會(huì)越多，從而計(jì)算精度越高，但是相應(yīng)地增大了計(jì)算量。

2.1 理論頻散曲線分析

首先從理論上論證各種數(shù)值方法在壓制數(shù)值頻散方面的能力?；舅枷霝椋和ㄟ^計(jì)算數(shù)值速度與實(shí)際波速的偏差考察各種離散方法的準(zhǔn)確性，偏差越小說明頻散越小。

為計(jì)算NAD方法的數(shù)值速度，首先引入均勻介質(zhì)中的波動(dòng)方程及其關(guān)于x、z方向的偏導(dǎo)

(20)

與式(5)的離散過程類似，通過NAD方法可將式(20)離散為矩陣形式

(21)

(22a)

(22b)

式中：a=exp(ikxh)；b=exp(ikzh)。

圖2 θ=0(a)和θ=π/6(b)時(shí)三種數(shù)值方法的頻散曲線

2.2 簡單模型的波場模擬

以均勻介質(zhì)模型(M1)和雙層介質(zhì)模型(M2)，進(jìn)一步討論NAD8方法的計(jì)算精度與數(shù)值效率。

M1模型的速度v=4km/s； M2模型的速度v1=4km/s、v2=5km/s，界面位于z=4km處，其他參數(shù)設(shè)置如表1所示。

表1 M1和M2模型的基本參數(shù)

圖3為頻率f=10和30Hz時(shí)， 本文構(gòu)造的NAD8方法針對以上兩種模型得到的頻率域單頻波波場快照。圖4為NAD8、NAD6和OFD8三種方法模擬的t=0.5s時(shí)刻時(shí)間域波場快照。從圖中可以看出，在相同網(wǎng)格規(guī)模的情況下， NAD8方法有效地壓制了數(shù)值頻散，而后兩種方法則出現(xiàn)了明顯的頻散現(xiàn)象，OFD8方法頻散最為嚴(yán)重，表明在相同網(wǎng)格頻率下， NAD8方法較其他方法精度較高。根據(jù)Shannon采樣定理[27]，為了準(zhǔn)確地刻畫波場，網(wǎng)格頻率至少大于2.0， 而NAD8方法突破了這個(gè)極限，Gf最小值可到達(dá)1.98，這是因?yàn)镹AD方法除了利用位移信息，還利用了波場梯度信息，大大提高了波場模擬精度，降低了對網(wǎng)格頻率的要求。

圖3 由八階NAD方法計(jì)算的M1(左)和 M2(右)模型頻率域單頻波波場快照(a)10Hz； (b)30Hz

圖4 三種方法計(jì)算的M1(左)和M2(右) 模型t=0.5s時(shí)刻的時(shí)間域波場快照(a)NAD8； (b)NAD6； (c)OFD8

最后，針對以上兩個(gè)模型，在不產(chǎn)生數(shù)字頻散的情況下，統(tǒng)計(jì)三種方法運(yùn)行時(shí)間(表2，10次計(jì)算結(jié)果的平均)，其中NAD8(Gf=1.98)用時(shí)最少，比NAD6方法(Gf=2.40)約省時(shí)12%，比OFD8方法(Gf=3.50)約省時(shí)25%。

表2 不同方法簡單模型波場模擬運(yùn)行時(shí)間對比

2.3 復(fù)雜模型的波場模擬

選取經(jīng)典的Marmousi模型[28](圖5)檢驗(yàn)NAD8方法的適用性。模型尺寸為17.25km×5.63km；網(wǎng)格間距h=18.75m，網(wǎng)格頻率Gf=3.0；四周均為10層PML邊界，震源位于(8.63km，0.38km)，主頻為20Hz。

圖6為NAD8方法模擬的頻率f=10Hz和15Hz時(shí)Marmousi模型的頻率域單頻波波場快照，圖7為t=1.00、1.67s時(shí)刻的時(shí)間域波場快照?？梢钥闯觯瑢τ趶?fù)雜的Marmousi模型，沒有明顯的頻散，可見NAD8方法適用于復(fù)雜介質(zhì)模型的正演模擬。

圖5 Marmousi速度模型

圖6 NAD8方法模擬的Marmousi模型單頻波波場快照(a)f=10Hz； (b)f=15Hz

圖7 NAD8方法模擬的Marmousi模型時(shí)間域波場快照(a)t=1.00s； (b)t=1.67s

3 頻率域全波形反演

3.1 基于NAD8方法的頻率域全波形反演

所謂反演即從初始模型出發(fā)，進(jìn)行正演模擬，匹配計(jì)算波場與真實(shí)波場，然后不斷修正模型，最終達(dá)到特定的精度。本文基于NAD8正演算法對經(jīng)典模型進(jìn)行反演，反演采用非線性共軛梯度法，其優(yōu)勢在于：與最速下降法相比，不用引入顯著的正演計(jì)算量便可以得到更加準(zhǔn)確、穩(wěn)健的計(jì)算結(jié)果，與牛頓類方法相比，不用計(jì)算反演目標(biāo)函數(shù)的二階偏導(dǎo)數(shù)矩陣(Hessian矩陣)，收斂性較強(qiáng)。將離散后的頻率域聲波方程(式(14))改寫為

T(ω)u(ω)=s(ω)

(23)

式中：T表示系數(shù)矩陣；u表示波場項(xiàng)；s表示震源項(xiàng)，且三者都與角頻率ω有關(guān)。合成波場與觀測波場數(shù)據(jù)殘差的L2范數(shù)可表示為

(24)

式中：m為模型參數(shù)； “*”表示復(fù)共軛；ns為炮數(shù)；nr為每炮接收點(diǎn)數(shù)。δdij表示第i炮的第j個(gè)檢波點(diǎn)對應(yīng)的模擬波場值u與觀測波場d的殘差。選擇合適的方法求解E(m)的最小值是問題的關(guān)鍵，基于梯度算法，通過推導(dǎo)計(jì)算可得相鄰步數(shù)之間的迭代基本格式為

(25)

通過對波動(dòng)方程關(guān)于模型參數(shù)求偏導(dǎo)，進(jìn)而可導(dǎo)出

(26)

式中v表示殘差反傳波場?？梢娪?jì)算反演目標(biāo)函數(shù)梯度時(shí)候，不需要計(jì)算矩陣J的具體表達(dá)式，整個(gè)問題只需計(jì)算TTv=δd*。復(fù)雜的矩陣計(jì)算可以變成簡單的求解一次線性方程組的問題，大大降低了計(jì)算量。

通過以上過程計(jì)算出梯度之后，進(jìn)而可以確定搜索方向。本文選用計(jì)算效率較高的拋物線擬合方法選擇步長，即近似認(rèn)為目標(biāo)函數(shù)是關(guān)于步長的二次函數(shù)，確定一個(gè)試探步長r1，使其滿足在目標(biāo)函數(shù)的下降方向上，不斷重復(fù)搜索步驟，直到目標(biāo)函數(shù)不降反增，確定另一步長r2； 結(jié)合初始步長r0，三點(diǎn)確定拋物線，求取極值點(diǎn)便為搜索步長。最后，反演頻率的選擇參照非等距頻率選擇策略[29]，這樣頻率間距隨著頻率的增加變得越來越大，從而整體上越來越松散，但不損失反演結(jié)果的分辨率，使用更少的頻率參與反演，可節(jié)省計(jì)算時(shí)間。

此外，定義標(biāo)準(zhǔn)化殘差的概念來衡量反演誤差的大小，即

(27)

式中：mt表示真實(shí)模型； ‖·‖2表示L2范數(shù)。顯然給定某一頻率，κ是一個(gè)隨迭代步數(shù)變化的曲線，曲線越趨近于0表示反演結(jié)果越好。

3.2 簡單模型反演

基于NAD8方法的正演程序，依據(jù)上文所述的頻率域反演原理，首先對兩個(gè)簡單的速度模型進(jìn)行反演與分析。

模型一(圖8)參數(shù)設(shè)置如下：模型尺寸為2.5km×2.5km，背景速度為4.0km/s，中央有一速度為4.5km/s、邊長為0.5km的異常體，網(wǎng)格間距h=25m； 將四周設(shè)置為10層PML邊界，主頻設(shè)置為10Hz，震源振幅為1.0×105(為了避免波場值太小而引起的機(jī)器出錯(cuò))，初始速度設(shè)置為4.0km/s。炮點(diǎn)與接收點(diǎn)布設(shè)兩組，位于計(jì)算區(qū)域內(nèi)。在頻率域接收點(diǎn)的多少并不會(huì)增大反演的計(jì)算量，而震源數(shù)目的增多會(huì)使計(jì)算量成倍增大。為保證計(jì)算效率，需要選擇適量的震源數(shù)目，而盡量多布設(shè)接收點(diǎn)，橫向間距為一個(gè)網(wǎng)格距離。為了得到深部更高分辨率的成像結(jié)果，分別設(shè)置兩排震源和接收器，設(shè)置震源數(shù)目為上、下各11個(gè)，共計(jì)22個(gè)，縱向上分別位于0.275km和2.225km處，橫向上為0.25～2.25km，間距為0.2km。上、下接收點(diǎn)同時(shí)接收傳播信號(hào)。

圖8 速度模型一及震源與接收點(diǎn)布設(shè)示意圖

首先選取5個(gè)頻率(1、7、13、19、25Hz)，低頻勾勒大致輪廓，高頻細(xì)化內(nèi)部結(jié)構(gòu)，每個(gè)頻率的最大迭代步數(shù)為30步，將反演過程記為S1，最終f=25Hz的反演結(jié)果如圖9a 所示?？梢钥闯龇囱莸乃俣饶Ｐ团c理論模型整體輪廓基本一致，但是分辨率相對較低，特別在靠近速度異常體周邊的數(shù)值與理論值差異較大。進(jìn)一步加密震源，每行放置21個(gè)，即2×21個(gè)，增加反演頻數(shù)到30個(gè)(1～30Hz，間隔1Hz)，每個(gè)頻率最大的迭代步數(shù)為50，將反演過程記為S2，最終f=30Hz的反演結(jié)果如圖9b所示。與圖9a相比，圖9b的分辨率明顯提高，細(xì)節(jié)刻畫更清晰，對異常體的刻畫更準(zhǔn)確。

圖9 模型一反演結(jié)果(a)S1； (b)S2

圖10中藍(lán)色曲線是頻率取25Hz時(shí)S1反演的誤差曲線，可以看出，當(dāng)步數(shù)到達(dá)最終30時(shí)，誤差曲線還有下降趨勢；紅色曲線為頻率取25Hz時(shí)S2反演的誤差曲線，表現(xiàn)為前幾步急劇下降，后期逐步平緩，停滯在小于0.1范圍內(nèi)且基本不再下降，1～30Hz中其他頻率的誤差曲線類似。

圖10 模型一反演誤差曲線

模型二(圖11)為雙層模型，尺寸為2.5km×2.5km，上層速度為4.00km/s，下層速度為4.50km/s在分界面處有一速度為4.25km/s、邊長為0.5km的異常體，網(wǎng)格間距為h=25m。將四周設(shè)置10層的PML邊界，震源主頻設(shè)為10Hz，震源振幅設(shè)置為1.0×105。反演的初始速度設(shè)置為4.00km/s。接收器盡可能多的放置，間距為一個(gè)網(wǎng)格距離，震源數(shù)目為上、下各布設(shè)41個(gè)，共計(jì)82個(gè)，縱向上分別位于0.275km和2.225km處，橫向上從0.25km到2.25km，間距為0.05km。上、下接收點(diǎn)同時(shí)接收傳播信號(hào)。

設(shè)置反演頻數(shù)為30個(gè)(1～30Hz，間隔為1Hz)，每個(gè)頻點(diǎn)的最大迭代次數(shù)為100。圖12a為f=30Hz的反演結(jié)果，可以看出，整體分辨率較高，細(xì)節(jié)刻畫清晰，特別是在速度分界面與異常體的形態(tài)反映較為真實(shí)，整體效果較為理想；圖12b為f=30Hz的反演誤差曲線，開始時(shí)下降明顯，40次后穩(wěn)定在0.1以內(nèi)，符合高精度反演的誤差曲線形態(tài)?？梢姡疚牡幕贜AD8的頻率域反演方法對于簡單的模型一和模型二取得了較為理想的結(jié)果。

圖11 速度模型二及震源與接收器布設(shè)示意圖

3.3 Marmousi模型反演

針對Marmousi模型(圖5)的反演，參數(shù)設(shè)置如下：網(wǎng)格間距h=15m，反演區(qū)域?yàn)?.47km×1.14km，PML層數(shù)為10，主頻設(shè)置為10Hz，震源振幅設(shè)為1.0×105。初始模型(圖13)是對真實(shí)模型進(jìn)行Gauss平滑后的結(jié)果，只具有大致的分層結(jié)構(gòu)。震源布設(shè)44個(gè)，縱向上位于0.17km處，橫向上從0.15km到3.35km， 間距為0.075km。僅在地表放置一排。接收器排列與震源的間隔相同，縱向位置上位于0.15km處，橫向上往外多擴(kuò)了兩個(gè)網(wǎng)格距離。

選擇30個(gè)頻點(diǎn)(1～30Hz，間隔為1Hz)、每個(gè)頻點(diǎn)迭代200次進(jìn)行反演，圖14為10和30Hz的反演結(jié)果。由圖可以看出，10Hz低頻反演結(jié)果基本給出了模型的大體輪廓，但是整體分辨率較低，出現(xiàn)了多處異常亮點(diǎn)，特別是深部高速異常體的精度較低。30Hz反演結(jié)果分辨率明顯提高、細(xì)節(jié)更加清晰，特別是深部高速異常體的刻畫更加精確。

圖12 模型二反演結(jié)果(a)及反演誤差曲線(b)

圖13 初始Marmousi速度模型及震源與接收點(diǎn)布設(shè)示意圖

圖14 Marmousi模型10(a)和30Hz(b)的速度反演結(jié)果

圖15是x=2.0km處、頻率為30Hz的速度反演曲線與真實(shí)、初始速度曲線的對比，可以看出，整體的反演結(jié)果較為理想，特別是淺部低速體保真度較高，但是深層精度略有下降，導(dǎo)致反演速度與真實(shí)模型速度之差的L2范數(shù)略高(約為0.8)，這是因?yàn)檎鹪次挥诘乇硖?，波傳播到深層時(shí)振幅降低、誤差累積，從而分辨率下降。

圖16為頻率為30Hz的反演誤差曲線，可以看出，整體上一開始下降明顯，并逐漸穩(wěn)定在0.1以內(nèi)，反演效果較好，其他頻率反演誤差曲線基本相似?？梢?，基于NAD8的頻率域全波形反演方法針對復(fù)雜模型依然適用，且分辨率、保真度都較高，效果良好。

圖15 x=2.0km處反演曲線與真實(shí)模型、初始模型的對比

圖16 Marmousi模型反演誤差曲線

4 結(jié)論

為進(jìn)一步提高頻率域全波形反演的計(jì)算效率，本文首次構(gòu)造了頻率域NAD8差分格式離散波動(dòng)方程，給出了PML邊界條件與頻率域NAD方法相結(jié)合的正演算法，并詳細(xì)推導(dǎo)總結(jié)了整個(gè)數(shù)值離散過程，得到了大型線性代數(shù)方程組，根據(jù)方程組系數(shù)矩陣的數(shù)學(xué)結(jié)構(gòu)，采用一類不精確旋轉(zhuǎn)分塊三角預(yù)處理算子加速Krylov迭代方法，提高了正演計(jì)算效率。數(shù)值頻散分析以及均勻介質(zhì)和雙層介質(zhì)模型的波場模擬結(jié)果表明：①與OFD8方法、NAD6方法相比，NAD8方法計(jì)算精度最高，在壓制頻散方面效果最明顯； ②在有效壓制頻散的情況下，NAD8方法所需網(wǎng)格頻率最小(Gf=1.98)，突破了Shannon采樣定理可達(dá)到的理論最小值(Gf=2.0)，所用的計(jì)算時(shí)間最短，大約比NAD6方法(Gf=2.4)省時(shí)12%，比OFD8方法(Gf=3.5)省時(shí)25%。

對于兩種典型的分層介質(zhì)模型及經(jīng)典Mar-mousi模型，本文基于NAD8格式的頻率域全波形反演算法獲得了較高分辨率和保真度的結(jié)果，結(jié)合反演誤差曲線驗(yàn)證了方法的正確性和適用性。