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

包乾宗 戴 雪 梁 雪

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

0 引言

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

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

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

1 原理方法

1.1 交錯網格有限差分震源波場重建方法

交錯網格有限差分離散格式為[30]

(1)

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

圖1 交錯網絡差分模板示意圖

(2)

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

(3)

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

基于平面波理論推導頻散關系，將波場表示為

pi=p0eikxh

(4)

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

(5)

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

(6)

(7)

1.2 重建系數優化

(8)

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

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

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

(9)

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

ψ1(kxh)]-1

(10)

建立基于無窮范數的極小化目標函數

(11)

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

(1)求解線性方程組

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

(12)

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

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

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

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

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

(6)基于約束方程式(8)得到優化的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為待求常數。初始采樣點(N+j+1個)可通過對波數范圍[0,βj]均勻采樣得到。

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

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

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

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

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

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

(14)

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

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

表2 本文方法的重建系數(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直接決定重建方法的精度。當M增大時，x=0和x=jh(j=2,3,…,M-1)層的精度增加，但x=h層的精度幾乎不變；當N增大時，x=h層的精度大幅改善。為了更好地比較，表3給出不同方法的有效帶寬βf。常規方法中采用優化差分系數[29]。隨著N的增大，有效波數范圍變大；當N=1時，本文方法的帶寬與常規方法(存儲M層)非常接近(約為2)。當一個波長采兩個樣點時，β=kxh達到尼奎斯特極限(等于π)。通過G=2π/βf將帶寬βf轉換為一個波長需要的采樣點數G，表4給出不同方法的G值。由表可知：本文方法的G值略大于常規方法。隨著N增大，G逐漸減小，當N=1時，一個波長只需要3個采樣點就可使重建誤差小于0.001。

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

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

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

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

3 穩定性分析

基于馮·諾依曼分析[27]可推導本文方法的穩定性條件。二維穩定性因子為

(15)

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

表5 不同方法的穩定性因子s

4 模型算例

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

4.1 聲波逆時偏移

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

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

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

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

圖4 聲波Marmousi模型

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

4.2 彈性波全波形反演

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

圖7 彈性Sigsbee2A模型

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

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

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

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

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

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

5 結論

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

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