一種基于新差分模板和無窮范數(shù)的震源波場重建方法

包乾宗 戴 雪 梁 雪

(①長安大學(xué)地質(zhì)工程與測繪學(xué)院地球物理系，陜西西安 710054；②海洋油氣勘探國家工程研究中心，北京 100028；③東方地球物理公司，河北涿州 072750)

0 引言

伴隨狀態(tài)法[1]通過兩次正演運算獲取模型參數(shù)的梯度，可避免對所有震源和接收點格林函數(shù)的直接計算和存儲，已廣泛應(yīng)用于地震波動方程偏移和反演，如逆時偏移(RTM)[2-5]、全波形反演(FWI)[6-10]、最小二乘逆時偏移(LSRTM)[11-13]、波動方程旅行時反演(WETI)[14-16]等。

伴隨狀態(tài)法通過震源波場和伴隨波場間的相互運算(如相關(guān))計算參數(shù)梯度/敏感核。但震源波場和伴隨波場分別沿時間正向和逆向傳播。為同時訪問震源波場和伴隨波場，震源波場需先保存在計算機(jī)內(nèi)存中或在波場反向傳播中進(jìn)行重建。存儲所有時刻、所有空間網(wǎng)格點的震源波場是不現(xiàn)實的，特別是三維情況。最優(yōu)檢測點方法(Optimal CheckPointing)在波場正向傳播中僅需保存幾個時刻(檢測點)的震源波場，在波場反向傳播中將存儲的檢測點波場作為初值重新計算相鄰檢測點之間的震源波場[17-19]。最優(yōu)檢測點方法可大幅減少計算機(jī)內(nèi)存消耗，但需額外的正演運算。邊界值方法采用位于邊界區(qū)域幾層(與差分精度有關(guān))網(wǎng)格點的波場和最后兩個時刻的快照重建震源波場[20-22]。常規(guī)邊界值方法需要存儲每個時刻M層網(wǎng)格點的邊界波場(M為有限差分算子長度)，存儲量仍然較大。Feng等[23]采用變階數(shù)有限差分算子重建震源波場，只需存儲一層空間網(wǎng)格點的邊界波場。但由于降階處理，該方法的重建精度不高。Tan等[24]提出基于拉格朗日多項式或高階波動方程外推法，該方法需保存一層或兩層網(wǎng)格點的邊界波場，內(nèi)存需求約為常規(guī)方法的 37.5%。Liu等[25]采用邊界波場的線性組合重建震源波場，所需存儲量僅為常規(guī)方法的1/M。段沛然等[26]將該方法推廣到交錯網(wǎng)格有限差分波場重建。Liu等[25]的方法存在如下問題：基于L2范數(shù)準(zhǔn)則求取差分系數(shù)，有效波數(shù)范圍/帶寬較窄；重建系數(shù)的自由度較小(最小為2)，很難精確逼近重建頻散關(guān)系。為保證重建誤差小于 0.0025，每個波長需要的采樣點數(shù)至少為 5。Ren等[27]設(shè)計了新的差分模板，采用邊界區(qū)域的N層波場和M-N層波場的線性組合來重建內(nèi)部區(qū)域的震源波場(0 ≤N≤M)。該方法通過調(diào)整N的大小實現(xiàn)重建精度和內(nèi)存需求的折衷。當(dāng)N=0時，退化為Liu等[25]的方法。Ren等[27]的方法每個波長只需3個采樣點就可保證重建誤差小于0.001。

Ren等[27]的方法是基于規(guī)則網(wǎng)格設(shè)計的，無法直接用于變密度聲波或彈性波速度—應(yīng)力方程RTM或FWI。本文發(fā)展了基于交錯網(wǎng)格有限差分的震源波場重建方法，設(shè)計了新的交錯網(wǎng)格差分模板，仍通過存儲邊界區(qū)域N層波場和M-N層波場的線性組合來重建震源波場；構(gòu)建了無窮范數(shù)極小化目標(biāo)函數(shù)，采用Remez交換算法[28-29]優(yōu)化了重建系數(shù)?；诟倪M(jìn)的差分模板和優(yōu)化的重建系數(shù)，新方法可大幅提高震源波場的重建精度，且內(nèi)存需求僅為常規(guī)方法的(N+1)/M。本文通過精度和穩(wěn)定性分析、聲波RTM和彈性波FWI算例驗證了方法的有效性。

1 原理方法

1.1 交錯網(wǎng)格有限差分震源波場重建方法

交錯網(wǎng)格有限差分離散格式為[30]

(1)

式中：ai為差分系數(shù)；h為空間網(wǎng)格間隔；p為地震波場(壓強(qiáng)或速度)。常規(guī)邊界值方法采用邊界區(qū)域的邊界波場p-i+0.5(i=1,2,…,M)重建內(nèi)部區(qū)域的震源波場[20-22]。為減小內(nèi)存需求，設(shè)計如圖1所示的差分模板。內(nèi)部區(qū)域各層(x=jh，j=0,1,…,M-1)的空間導(dǎo)數(shù)通過下式計算

圖1 交錯網(wǎng)絡(luò)差分模板示意圖

(2)

式中：ai(i=1,2,…,M)和bj,i(j=1,2,…,M-1,i=1,2,…,N+j+1)統(tǒng)稱為重建系數(shù)。式(2)可改寫為

(3)

波場正向延拓中存儲邊界波場p-i+0.5(i=1,2,…,N)和波場的線性組合A，波場反向延拓中基于式(3)重建內(nèi)部區(qū)域不同層(x=jh，j=0,1,…,M-1)的空間偏導(dǎo)數(shù)。所提方法只需保存N層邊界波場和波場線性組合A，內(nèi)存需求為常規(guī)邊界值方法的(N+1)/M。該方法可通過調(diào)節(jié)N控制重建精度和存儲需求。隨著N增大，重建系數(shù)的個數(shù)增加(自由度增大)，精度提高，但存儲量隨之增加。當(dāng)N=M或0時，本文方法退化為常規(guī)邊界值方法[20-22]或Liu等[25]的方法。

基于平面波理論推導(dǎo)頻散關(guān)系，將波場表示為

pi=p0eikxh

(4)

式中kx為x方向的波數(shù)。將式(4)代入式(3)，有

(5)

式中j=1,2,…,M-1。式(5)為新重建方法的頻散關(guān)系，定義如下的相對誤差用于后續(xù)重建精度分析

(6)

(7)

1.2 重建系數(shù)優(yōu)化

(8)

式(8)為約束方程，可保證低波數(shù)成分的重建精度。

將式(8)代入式(6)，可得

δ0(kxh)=φ1(kxh)+

(9)

δj(kxh)=φ1(kxh)+

ψ1(kxh)]-1

(10)

建立基于無窮范數(shù)的極小化目標(biāo)函數(shù)

(11)

式中βj為最大波數(shù)(j=0,1,…,M-1)。利用式(11)無法直接求取導(dǎo)數(shù)，需采用Remez交換算法求解，詳細(xì)步驟如下。

(1)求解線性方程組

δ0(kxlh)=(-1)lλ0l=1,2,…,M

(12)

得到ai(i=2,3,…,M)和λ0。式中：λ0為待求常數(shù)；kxl為kx的第l個采樣點。M個初始采樣點通過對波數(shù)范圍[0,β0]均勻采樣得到。

(2)采用式(11)計算E0。當(dāng)E0>η(η為最大允許誤差)，執(zhí)行步驟(3)；當(dāng)E0≤η，執(zhí)行步驟(6)。

(3)式(12)中重建誤差在采樣點處正負(fù)交替，M個采樣點之間存在M-1個根，M-1個根將[0,β0]分成M個區(qū)間，搜索得到M個區(qū)間的極值點。

(4)將M個極值點作為新的采樣點，重復(fù)步驟(1)～步驟(3)，直到滿足收斂條件E0≤η，輸出優(yōu)化的ai(i=2,3,…,M)和β0。

(5)若達(dá)到最大迭代次數(shù)，減小β0，重復(fù)步驟(1)～步驟(4)。

(6)基于約束方程式(8)得到優(yōu)化的a1。

(7)求解下式得到bj,i(i=2,3,…,N+j+1)和λj(j=1,2,…,M-1)。

δj(kxlh)=(-1)lλjj=1,2,…,M-1；

l=1,2,…,N+j+1

(13)

式中λj為待求常數(shù)。初始采樣點(N+j+1個)可通過對波數(shù)范圍[0,βj]均勻采樣得到。

(8)采用式(11)計算Ej(j=1,2,…,M-1)。當(dāng)Ej>η，執(zhí)行步驟(9)；當(dāng)Ej≤η，執(zhí)行步驟(12)。

(9)式(13)中重建誤差在采樣點處正負(fù)交替，N+j+1個采樣點之間存在N+j個根，N+j個根將[0,βj]分成N+j+1個區(qū)間，搜索得到N+j+1個區(qū)間的極值點。

(10)將N+j+1個極值點作為新的采樣點重復(fù)步驟(7)～步驟(9)，直到滿足收斂條件Ej≤η，輸出優(yōu)化的bj,i和βj。

(11)若達(dá)到最大迭代次數(shù)，減小βj，重復(fù)步驟(7)～步驟(10)。

(12)基于約束方程式(8)得到優(yōu)化的bj,1。

對于不同的層x=jh(j=0,1,…,M-1)，優(yōu)化的最大波數(shù)βj是不同的。用于評價所提方法的重建精度的有效波數(shù)范圍/帶寬為

(14)

上述優(yōu)化步驟可保證所提方法的重建誤差在[0,βf]范圍內(nèi)均小于最大允許誤差η。表1和表2分別給出N=0和1時的重建系數(shù)。

表1 本文方法的重建系數(shù)(N=0，M=6和η=0.001)

表2 本文方法的重建系數(shù)(N=1，M=6和η=0.001)

2 精度分析

采用式(6)分析重建方法的精度。圖2和圖3分別為N=0和1時的頻散曲線，可見本文方法的頻散曲線呈波紋狀分布。與δ0和δj(j=2,3,…,M-1)相比，δ1(x=h層的精度)的幅值最先超出最大允許誤差(η= 0.001)。因此，δ1直接決定重建方法的精度。當(dāng)M增大時，x=0和x=jh(j=2,3,…,M-1)層的精度增加，但x=h層的精度幾乎不變；當(dāng)N增大時，x=h層的精度大幅改善。為了更好地比較，表3給出不同方法的有效帶寬βf。常規(guī)方法中采用優(yōu)化差分系數(shù)[29]。隨著N的增大，有效波數(shù)范圍變大；當(dāng)N=1時，本文方法的帶寬與常規(guī)方法(存儲M層)非常接近(約為2)。當(dāng)一個波長采兩個樣點時，β=kxh達(dá)到尼奎斯特極限(等于π)。通過G=2π/βf將帶寬βf轉(zhuǎn)換為一個波長需要的采樣點數(shù)G，表4給出不同方法的G值。由表可知：本文方法的G值略大于常規(guī)方法。隨著N增大，G逐漸減小，當(dāng)N=1時，一個波長只需要3個采樣點就可使重建誤差小于0.001。

圖2 本文方法不同有限差分算子的頻散曲線(N=0，η=0.001)

圖3 本文方法不同有限差分算子的頻散曲線(N=1和η=0.001)

表3 不同方法的有效帶寬βf(η=0.001)

表4 不同方法一個波長需要的采樣點數(shù)G(η=0.001)

3 穩(wěn)定性分析

基于馮·諾依曼分析[27]可推導(dǎo)本文方法的穩(wěn)定性條件。二維穩(wěn)定性因子為

(15)

式中：j=1,2,…,M—1，s為允許的最大庫朗數(shù)，穩(wěn)定性隨著s增大逐漸改善。表3給出不同方法的穩(wěn)定性因子。與常規(guī)方法相比，本文方法的穩(wěn)定性條件更加嚴(yán)格。隨著N的增加，穩(wěn)定性略微變差。

表5 不同方法的穩(wěn)定性因子s

4 模型算例

將本文方法應(yīng)用于聲波逆時偏移和彈性波全波形反演進(jìn)行檢驗。時間二階、空間十二階優(yōu)化交錯網(wǎng)格有限差分[29]用于波場正、反向延拓。常規(guī)邊界值重建方法采用優(yōu)化差分系數(shù)[29]。正向震源波場、常規(guī)方法得到的像和反演結(jié)果作為參考解，計算不同方法結(jié)果的最大絕對誤差εMAV和均方根誤差εRMS。

4.1 聲波逆時偏移

聲波Marmousi模型如圖4,時間步長為1.0ms，記錄時間為4.0s，空間網(wǎng)格大小為10m×10m，網(wǎng)格點數(shù)為501×353，震源為20Hz主頻的Ricker子波。51炮和501個檢波點均勻分布于地表。圖5為不同方法重建的波場，炮點位于(9.5km,0)。表6給出不同方法重建波場的εMAV和εRMS。由表可知，常規(guī)方法的重建誤差非常小，本文方法的εMAV和εRMS分別在10-2和 10-4數(shù)量級。重建波場的誤差隨著N的增大有所減小。圖6為逆時偏移中采用不同重建方法得到的像，表7為本文方法的εMAV和εRMS，可見，N=0或1都能保證成像結(jié)果的均方根誤差小于10-3。當(dāng)N=0或1時，本文方法的存儲量為常規(guī)方法的16.6%或33.3%。為兼顧重建精度和內(nèi)存需求，聲波RTM建議采用N=0。

圖6 聲波Marmousi模型不同方法的逆時偏移結(jié)果

表6 聲波Marmousi模型不同方法重建波場的ε MAV和ε RMS

表7 聲波Marmousi模型本文方法不同N的ε MAV和ε RMS

圖4 聲波Marmousi模型

圖5 聲波Marmousi模型不同方法重建的波場

4.2 彈性波全波形反演

彈性波Sigsbee2A模型如圖7所示?？臻g網(wǎng)格尺寸為16m×16m，網(wǎng)格點數(shù)為351×186，時間步長為1.5ms，記錄長度為4.5s。爆炸震源激發(fā)主頻為15Hz的Ricker子波。36個震源和351個檢波點均勻分布于地表。初始縱波和橫波速度模型為隨深度線性變化的一維模型(圖8)。采用多尺度反演(0～5Hz、0～10Hz和0～15Hz)策略，每個尺度上迭代10次。圖9給出不同方法得到的波場，炮點位于(2720m，0)，表8給出不同方法重建波場的εMAV和εRMS。常規(guī)方法(存儲M層)的誤差在10-6數(shù)量級，本文方法的誤差隨N的增大而減小。圖10為采用不同重建方法進(jìn)行全波形反演得到的縱、橫波速度模型。表9為本文方法反演結(jié)果的εMAV和εRMS。當(dāng)N=0時，反演結(jié)果的精度較差；當(dāng)N=1時，反演精度得到大幅提高(εRMS<0.5 %)。本文方法的存儲需求僅為常規(guī)方法的(N+1)/M。波場重建誤差對模型參數(shù)梯度的精度有影響，為保證反演精度，全波形反演建議采用N=1。

圖7 彈性Sigsbee2A模型

圖8 彈性Sigsbee2A模型的初始模型

圖9 彈性Sigsbee2A模型不同方法重建的1.5s波場

圖10 彈性Sigsbee2A模型不同方法重建波場的全波形反演結(jié)果

表 8 彈性Sigsbee2A模型不同方法重建波場的ε MAV和ε RMS

表9 彈性Sigsbee2A模型本文方法不同N時反演結(jié)果的ε MAV和ε RMS

與常規(guī)方法相比，本文方法可降低存儲量，但重建波場的誤差有所增加(詳見精度分析和模型算例)。實際應(yīng)用中，可通過調(diào)整N的大小進(jìn)行重建精度和內(nèi)存需求的折衷。由式(1)和式(3)可知，常規(guī)方法的計算復(fù)雜度為3M2，本文方法計算復(fù)雜度為(7M2+2MN+3M-4)/2。當(dāng)N為0或1時，本文方法的計算量約為常規(guī)方法的7/6倍。與波場正、反向延拓相比，本文震源波場重建方法增加的計算量微不足道。

5 結(jié)論

本文提出了一種交錯網(wǎng)格有限差分震源波場重建方法。該方法通過存儲N層邊界波場和M—N層波場的線性組合重建震源波場。為提高重建精度，基于無窮范數(shù)準(zhǔn)則和Remez交換算法優(yōu)化了重建系數(shù)。詳細(xì)分析了提出方法的精度和穩(wěn)定性，并將其應(yīng)用于聲波逆時偏移和彈性波全波形反演。新方法可通過調(diào)整N的大小實現(xiàn)重建精度和存儲需求的平衡。新方法能夠獲得較好的重建波場、像和反演結(jié)果，且內(nèi)存需求僅為常規(guī)方法的(N+1)/M。

本文方法可直接應(yīng)用于三維逆時偏移和全波形反演或其他基于伴隨狀態(tài)法的地球物理問題中，如最小二乘逆時偏移、波動方程旅行時反演或全球尺度波形層析等。但新方法無法用于基于黏聲或黏彈波動方程的成像和反演中。如何避免衰減介質(zhì)波場重建的不穩(wěn)定性需要進(jìn)一步考慮。