全聲波方程波場深度延拓倏逝波壓制研究

魏俊廷, 尤加春, 趙金泉, 陳軍兆, 代中奎

成都理工大學地球物理學院, 成都 610059

0 引言

地震深度偏移是地震勘探中的一個重要研究領域(劉喜武和劉洪, 2002;劉旭明等, 2021).理論上,深度偏移包括了震源波場正向延拓和檢波器波場反向延拓,然后再使用互相關成像條件來得到成像剖面(馬在田, 2002).在石油勘探中,深度偏移的主要目的是為了確定油氣資源分布并提供一個高質量的地下構造剖面(李振春, 2014;張選朋等, 2022).一般而言,深度偏移方法大致可以分為三類:單程波方程深度偏移方法(Gazdag and Sguazzero, 1984; Wu, 1994; Ristow and Rühl, 1994; Le Rousseau and de Hoop, 2001; 張宇, 2006)、逆時偏移方法(RTM)(Baysal et al., 1983; 劉紅偉等, 2010;楊仁虎等, 2010;楊仁虎,2021)和全聲波方程波場深度延拓及成像方法(Kosloff and Baysal, 1983; You et al., 2016, 2018; You and Cao, 2020).

全聲波方程波場深度延拓及成像方法和RTM均是基于全聲波方程,而傳統單程波方程深度偏移方法是基于單程波方程,其中單程波方程是對全聲波方程的一種近似.因此,與其他兩種深度偏移方法相比,單程波偏移方法難以有效地描述波場傳播中的動力學特征(張宇等, 2007).全聲波方程波場深度延拓及成像方法和RTM在算法方面都有著各自的特點:(1)雖然RTM和全聲波方程波場深度延拓及成像方法都是基于全聲波方程,均是對聲波方程進行求解,但兩者的內涵完全不一樣.RTM是在時間域求解聲波方程,波場延拓是沿著時間軸展開;而全聲波方程波場深度延拓及成像方法是在深度域求解聲波方程,波場延拓是沿著深度軸展開,如圖1所示.因此在深度域開展波場深度延拓與在時間域開展波場時間延拓在理論上有明顯的優勢.RTM需要至少存儲一個傳播時間方向上的全部波場,對二維情況而言,即要存儲(nx,nz,tn)維度的數據(其中nx和nz分別表示速度模型在x和z兩個方向的采樣點,tn表示時間采樣點),而深度延拓只需要存儲一個深度上的波場(nx,ωn)(其中ωn為頻率采樣點,ωn

圖1 深度偏移與RTM波場延拓示意圖

因此在分析前人關于倏逝波壓制理論和方法的基礎之上,本文嘗試提出新的倏逝波壓制策略,期望能夠在計算效率和計算精度上達到較好的平衡.因此本論文的組織結構如下:首先,在You和Cao(2020)提出的波場深度延拓矩陣計算的基礎上,本文通過引入一個新的壓力波場來代替傳播算子中壓力對深度的偏導數,推導建立了一種新的波場深度延拓策略;其次,基于本文所提新的波場深度延拓策略,在頻率-波數域開發兩種新的倏逝波壓制策略;最后,建立一系列理論模型并開展一系列數值實驗,進一步證明與其他倏逝波壓制方法相比,本文所提出的全聲波方程波場深度延拓及成像方法在對倏逝波壓制和偏移成像方面具有更好的穩定性和準確性.

1 基本原理

1.1 全聲波方程波場深度延拓策略

二維情況下常密度頻率域聲波方程可以表示為

(1)

(2)

(3)

(4)

其中y=kzΔz,i為虛數,即i2=-1.

1.2 倏逝波的產生原理

為了引出倏逝波,本文以二維情況下聲波方程為例(Sandberg and Beylkin, 2009):

(5)

(6)

因為當k>2πω/v時kz為虛數,可以觀察到式(6)中的一部分函數呈指數型增長,導致方程解也呈指數型增長,呈指數型增長的這一部分函數所對應的波場就是倏逝波.對于特征值λk,其由λk=-4π2ω2/v2+k2(k=0,1,…)所給出,當k>2πω/v時,特征值λk為正,對應著倏逝波的產生.倏逝波會使得雙程波深度偏移算法出現不穩定,甚至會導致雙程波深度偏移方法難以有效成像.從式(6)可知,倏逝波的產生是與特征值密切相關,在實際應用中特征值的計算是一個比較費時費力的過程,特別是對于大型矩陣,因此,開發一種快速、準確的倏逝波壓制技術對于雙程波方程波場深度延拓及成像是一個亟待解決的關鍵問題.

1.3 全聲波方程波場深度延拓中的倏逝波壓制

對基于全聲波方程波場深度延拓而言,其要面臨一個特殊的挑戰:如何準確、高效地壓制倏逝波.目前針對全聲波方程波場深度延拓中的倏逝波壓制問題,有兩種解決方案:第一種解決方案是利用最大速度值來構建低通濾波器,該方法計算效率高,但也會去除一些對陡傾角構造成像有價值的高波數分量(Kosloff and Baysal, 1983),詳見附錄B;第二種解決方案是譜投影法(Sandberg and Beylkin, 2009);它可以去除Helmholtz算子中與倏逝波相關的負特征值成分,但該方法計算效率低,不利于實際應用,詳見附錄C.

基于本文所建立的新波場深度延拓策略即式(4),本文提出了兩種新的倏逝波壓制策略.策略I是廣義低通濾波器,策略II是能量守恒濾波器.

1.4 策略I:廣義低通濾波器

為了更好地介紹廣義低通濾波器,接下來,本文首先證明單程波深度偏移算子和上述全聲波方程深度延拓算子之間的相似性.然后在此基礎上,本文開發了一種新策略來壓制倏逝波,該策略是由經典低通濾波器的擴展形成(Kosloff and Baysal, 1983),更適用于強橫向速度變化介質模型的準確波場延拓計算.

根據歐拉公式,式(4)中的cosy和isiny函數可以寫為

(7)

基于式(7),全聲波方程波場深度延拓即式(4)可以修改為

(8)

在式(8)中可以發現,雙程波方程中的波場深度延拓可表示為兩個指數函數的線性組合.在式(8)中,當使用歐拉公式將三角函數轉化為指數函數形式時,可以發現波場深度延拓能由指數函數線性表示.通過這種轉變,由式(4)表示的矩陣算子可等效地由指數函數表示,其中指數函數與倏逝波具有直接關系.與單程波方程深度偏移方法類似,本文提出利用指數函數在頻率-波數域中開展倏逝波的識別與壓制工作,其算法可以進一步描述為

(9)

在式(9)中,當速度在橫向(x方向)上發生變化,僅用一個速度值很難精確解決倏逝波壓制問題,特別是在橫向速度變化劇烈的情況.因此,本文在延拓深度的每一個橫向速度上利用式(9)開展倏逝波壓制.為了應對在速度橫向變化情況下倏逝波的壓制問題,本文提出了一種新的倏逝波壓制算法,該算法具有以下關鍵操作——在每個頻率ω和每個延拓深度z,對x(橫向)方向上每個速度進行獨立的倏逝波識別壓制,具體算法策略如下:

從算法的執行策略上看,本文所提新算法是對經典低通頻率-波速域濾波器方法的擴展,特別適用于并行計算.

1.5 策略II:能量守恒濾波器

能量守恒濾波器是以能量守恒原理為基礎,式(4)可以等價地表示為如下形式:

(10)

從式(10)可以看出,本文的傳播算子A是一個相移矩陣,在深度延拓中,它僅改變波場的相移量,而不改變波場的能量(波場矢量范數).因此,對于正常的波場傳播過程,能量在傳播過程中會根據能量守恒定律而保持不變;但對于倏逝波,能量因為不遵守能量守恒定律而出現振幅指數增長情況,導致波場深度延拓算法的不穩定.

(11)

其中α是比例因子,在本文的數值實驗中,比例因子設置為1.02～1.03.

基于上述描述,本文提出了一種“能量守恒濾波器”方法,為了描述該算法的實現,在每個頻率和每個延拓深度都使用以下步驟來壓制倏逝波:

(3)應用式(11)中表示的能量范數方法在頻率-波數域中識別并壓制倏逝波.

對于策略I,在實際應用中,本文可逐點對x方向上每一速度點都使用該策略進行倏逝波的識別和壓制,但是采用“逐點”策略的情況下,如果速度模型的橫向速度維度很大(x方向的速度點很多),會影響算法的執行效率.為了解決該問題,本文采用“群組”或者“逐點-群組”的方式來提高算法的執行效率.其中“群組”算法可以理解為:將每個深度延拓的速度曲線分為k組,數據點總數為N,每一組具有Nk=N/k個數據點.在每一組速度中,用局部最大速度來表示這一組數據,再采用插值法實現對群組內每一個速度點的波場計算.而對“逐點-群組”策略算法可以理解為是一種混合算法,即在橫向速度變化較弱的區域采用“群組”算法策略,在橫向速度變化較強的區域采用“逐點”算法策略.因為在實際應用中,偏移過程中使用的速度模型的分辨率受速度分析精度的影響,通過速度分析方法得到的速度模型相當于是對精確模型進行平滑處理后所得到的結果,因此在本文中使用的“群組”策略是合理的,后面的數值實驗也證明了該策略執行的有效性.

此外,因為本文所提倏逝波壓制策略都是對每一個速度點進行波場計算,因此可以進一步采用線程并行化的方式來提高計算效率.其中對指數函數的并行化計算示意圖如圖2所示.在并行化實現過程中,本文給每一個內核分配一個并行化池來處理橫向速度變化中的一個速度值及該速度值對應的指數函數.在后面的數值實驗中,本文還將深入討論并行化計算對偏移成像處理計算效率的影響.

圖2 并行計算模式圖

2 數值實驗

在本節中,本文建立若干個理論速度模型來研究所提策略在壓制倏逝波上的計算效果.為了更好地對比不同算法的計算性能,本文將所提算法與若干個常規偏移方法進行性能對比,這些常規偏移方法包括:單程波廣義屏傳播算子(GSP)、基于低通濾器的全聲波方程深度偏移方法和基于譜投影方法的全聲波方程深度偏移方法.

2.1 梯度速度模型中的脈沖響應

為了測試各種偏移方法對橫向速度變化模型的適應性,本文構建了一個梯度速度模型,如圖3所示.該速度模型的維度大小尺寸為1500 m×3000 m,其水平和垂直網格間距均為5.0 m.用于計算脈沖響應的震源位于x=750 m和z=0 m處,波場模擬使用的震源為20 Hz的雷克子波.

圖3 梯度速度模型圖

為了對比不同偏移方法計算脈沖響應的差異,本文使用了經典有限差分方法(FD)、單程波廣義屏方法(GSP)和本文所提策略來計算波場傳播,結果如圖4所示.在計算的波場中,由于有限差分方法能適應任意橫向速度變化的介質,因此本文將有限差分方法計算的結果作為參考標準.本文將由有限差分方法計算的波前面提取出來,與單程波廣義屏方法和本文所提策略計算出的波前面繪制在一起,便于進行對比研究.

圖4 不同方法計算的脈沖響應圖

從圖4中可以清晰地觀察到,除了在非常大的傳播角度區域外,本文所提策略計算出的波前面與有限差分方法計算出的波前面吻合良好,而單程波廣義屏方法所計算出的波前面僅在一定角度內與有限差分方法所計算出的波前面擬合較好.由于該速度模型具有較強烈的速度變化,對于單程波偏移方法而言,在傳播角度非常大的情況下難以準確描述波場傳播,這也是單程波偏移方法面臨的一個難題.至于在非常大的傳播角度(接近90°)處出現本文所提策略計算的波前面與有限差分方法計算出的波前面不一致的情況,本文認為在此處kz接近于0,在這時將部分有用波場當成倏逝波進行了去除,導致出現兩個波前面吻合不好的情況.

在該梯度速度模型中,本文也對所提策略開展了并行化性能測試.在該并行化測試中,本文分別設置了4、8、16和28個CPU內核來并行化計算波場傳播并統計其運行時間.并行化計算時間見表1和表2.顯然,隨著并行計算核心數的逐步增加,波場計算的運行時間也相應減少.

表1 基于策略I的全聲波方程深度偏移方法在不同CPU核心數的運行時間

表2 基于策略II的全聲波方程深度偏移方法在不同CPU核心數的運行時間

從表1和表2可以充分地展示本文設計的算法能充分地發揮多核CPU的并行化計算能力.利用并行計算能力,甚至可以將本文所提策略移植到Cuda GPU上并行運行.從統計的運行時間上不難發現本文所提策略II比策略I所需更少的運算時間,也體現了本文所提能量守恒濾波器方法在壓制倏逝波中的優越性.

2.2 鹽丘模型成像

在數值實驗中,本文使用鹽丘模型來比較不同策略在壓制倏逝波方面的性能,因為鹽丘模型中有一個高速陡傾角鹽丘體,其速度大約是圍巖的2倍,所以其是評價成像性能的經典模型.本文使用兩種鹽丘模型,第一種是準確的鹽丘模型,用于產生偏移成像的炮道集;第二種是使用3×3高斯濾波窗口產生的平滑速度模型,用于偏移成像.兩種速度模型如圖5所示.速度模型小大為4400 m(z)×14000 m(x),水平和垂直方向網格間距均為20.0 m.在波場模擬中,采用左邊放炮右邊接收的觀測系統,共計算224個炮記錄,炮間距60.0 m,每炮設置240個檢波器,檢波器間距20.0 m,最小偏移距0 m.波場模擬使用的震源為10 Hz的雷克子波,時間采樣間隔為0.001 s,最長記錄時間為5.0 s.

圖5 鹽丘速度模型圖

在本偏移成像處理中,本文對比了5種偏移方法的成像性能,其中包括:經典的單程波廣義屏方法、基于低通濾波器方法的全聲波方程深度偏移方法、基于譜投影方法的全聲波方程深度偏移方法、基于策略I和II的全聲波方程深度偏移方法.使用不同偏移方法計算的成像剖面如圖6所示.

圖6 利用不同偏移方法計算得到的成像剖面

由于鹽丘模型的橫向速度變化較大,因此單程波廣義屏方法在波場計算中的能力有限,尤其是對于大角度傳播的波場,這一現象可以在圖6a中清晰地觀察到.因此,單程波廣義屏方法無法對鹽丘模型下的一些斷層進行有效成像,如圖6c、d中紅色箭頭所示;并且單程波廣義屏方法對于水平構造的成像不如其他偏移成像方法那樣準確,會出現彎曲的情況,如圖6a中黑色箭頭所示.此外,如圖6a中,也能觀察到鹽丘下方出現了一些虛假成像軸.當使用經典的低通濾波方法壓制倏逝波時,一些有用的高波數域成分也被去除,這樣會導致難以對鹽丘模型中的一些陡傾角構造結構進行有效地成像,如圖6b中的藍色箭頭所示.

相反,譜投影方法可以有效地壓制倏逝波,從而可以很好地保留有效波場.然而,由于譜投影方法涉及了密集矩陣的乘法運算,因此譜投影方法的計算效率限制了其廣泛運用.本文使用的策略I和策略II在每一個水平方向的速度值上對倏逝波進行準確地識別和壓制.因此,本文所提策略比經典的低通濾波方法保留了更多有用的高波數域信息.如圖6d、e所示,本文所提策略可以成功地壓制倏逝波并對鹽丘模型進行準確成像.與低通濾波器成像的結果相比(圖6b),不難發現本文所提策略計算得出的剖面質量遠高于經典低通濾波方法.

然而,在對比本文所提策略I和策略II計算得出的成像剖面時,在速度急劇變化的構造位置,策略II計算的偏移結果比策略I計算的偏移結果會產生更多的虛假成像,如圖6e中的綠色箭頭所示.對于速度橫向變化不劇烈的介質,例如圖5a中1500 m以上的構造,策略II表現出良好的成像性能.

除了關注成像質量之外,本文還關注基于不同倏逝波壓制策略的全聲波方程深度偏移方法的計算效率.在本數值實驗中,本文使用了16個CPU計算核心來完成單炮記錄的深度偏移,并統計不同偏移方法的計算時間,如表3所示.

表3 不同偏移方法的計算時間

從表3中可以明顯地發現,基于低通濾波的聲波方程深度偏移方法比基于譜投影方法的聲波方程深度偏移方法的計算時間要少得多.顯然,基于密集的遞歸矩陣乘法運算是相當費時的,因此基于譜投影方法的全聲波方程深度偏移方法需要更多的計算時間來識別和壓制倏逝波.

與這兩種經典的倏逝波壓制方法相比,本文所提策略I和II比譜投影方法需要更少的計算時間,同時實現了與譜投影方法相似的計算精度;與策略I相比,本文所提策略II由于操作簡單,因此具有更少的計算時間,有更高的計算效率.

為了進一步證明本文所提策略在壓制倏逝波方面的可行性,本文提取了z=2000 m深度處用隨機函數所生成的隨機波場,然后使用低通濾波器、譜投影和本文所提策略來延拓該隨機波場,圖7a為深度z=2000 m的鹽丘速度曲線.為了便于對比分析,本文使用經過譜投影方法處理得到的波場作為參考標準.為了進行定量地對比分析,本文計算由譜投影方法生成的波場與由低通濾波器和本文所提策略生成的波場之間的均方根誤差,波形誤差如圖7b所示.通過觀察這些曲線,本文可以很清晰地發現,與低通濾波器和策略II相比,由本文所提策略I計算的曲線具有最低的誤差.對于圖中的高速部分(圖7a中的鹽丘位置),策略I和低通濾波器方法的曲線非常接近.這是因為低通濾波器方法使用的是速度模型中最大的速度來進行濾波.策略II在高速部分(鹽丘體)顯示出相對較高的誤差,這也解釋了策略II會在圖6e中鹽丘部分產生噪聲的原因.

圖7 均方根誤差圖

本文除了關注在鹽丘模型的成像精度和質量,還關注所提策略在實際應用中的計算效率.在數值實驗中,本文發現 “群組”策略不僅在計算精度上與“逐點”策略類似,同時還具有更高的計算效率.在對鹽丘模型的偏移成像實驗中,速度模型的水平方向共有14000個速度樣本點,本文分別將速度樣本點分成40、60和80個“群組”進行測試;然后對同一炮的偏移結果,本文統計了其各自的計算時間,如表4所示.

表4 基于策略I的不同“群組”大小的計算時間

圖8a、b、c分別是采用不同“群組”大小的成像結果.將圖8的這些成像結果與“逐點”圖(6d)策略計算得出的成像結果進行對比,可以清楚地發現采用“群組”策略算法顯著提高了“逐點”策略算法的計算效率,如表4所示;但是在采用“群組”策略算法時當“群組”的大小太大,會導致在橫向速度變化劇烈位置的成像不準確,如圖8中的紅色箭頭所示.同時,對于橫向速度變化較弱的區域,“群組”策略可以在保持較高成像精度的情況下,有效提高計算效率.為了進行更深入地對比研究,本文對混合“逐點-群組”策略算法開展了成像性能和效率測試,成像得到的結果如圖8d所示,計算時間見表4.通過對比成像結果和計算精度,本文所提混合“逐點-群組”策略算法也對鹽丘模型進行了準確地成像,同時也具有較高的計算效率,從而達到計算精度與計算效率的平衡.

圖8 基于策略I,通過使用不通的“群組”大小所得到的成像剖面

3 結論

雙程波方程波場深度延拓及成像是地震成像中重要研究領域,相較于其他偏移方法,倏逝波問題是雙程波方程波場深度延拓中面臨的特殊挑戰之一,嚴重制約該方法在實際中的生產應用,為實現高效、準確的倏逝波壓制,本文提出了廣義低通濾波器(策略I)和能量守恒濾波器(策略II).由于上述兩種所提出策略均可單獨地處理橫向變化的每一個速度值,因此,本文建立了上述兩種策略的并行程序,為實際應用提供基礎工具.通過對鹽丘模型的成像,本文所提出策略克服了基于經典低通濾波器方法難以對陡傾角構造成像的缺點,可實現譜投影方法相當的倏逝波壓制精度而且計算效率更高.進一步深入開展對鹽丘模型的成像質量和成像效率的對比,本文發現策略I的成像質量要優于策略II的成像質量,但是策略II計算效率要高于策略I.為了進一步提高策略I的計算效率,本文提出了混合“逐點-群組”策略算法,通過對該算法的數值試驗,本文發現混合“逐點-群組”策略算法能對計算精度和計算效率做到更好的平衡,具有更加廣闊的應用前景.

附錄A 雙程波方程波場深度延拓中的快速波場分解

二維情況下常密度頻率域聲波方程(式(1))的通解形式可以寫為

×e-ikz(z-z0),

(A1)

對式(A1)中深度變量求偏導數,則有:

(A2)

×e-ikz(z-z0).

(A3)

對式(A1)、(A3)進行加法和減法運算,則有:

2pd(x,z,ω)eikzΔz,

(A4)

2pu(x,z,ω)e-ikzΔz,

(A5)

附錄B 最大速度低通濾波器方法

(B1)

其中kx為水平波數,ω為角頻率,cmax為延拓深度處的最大速度值.

附錄C 譜投影方法

以二維情況下常密度頻率域聲波方程為例,常規的雙程波方程波場深度延拓方程可以寫為

(C1)

Helmholtz算子L包含有正特征值和負特征值,而負特征值對應的是對倏逝波的傳播.為了對倏逝波進行壓制,Sandberg和Beylkin(2009)提出了一種矩陣迭代計算的譜投影方法,其中譜投影算子表示為

P=[I-sign(L)]/2,

(C2)

其中I為單位矩陣.

在譜投影方法壓制倏逝波的計算中,為了計算矩陣的符號函數sign(L),Sandberg和Beylkin(2009)提出采用矩陣迭代計算的方式,該方法具體執行流程可寫為:

(1)初始化:S0=L/‖L‖2;

(3)滿足迭代條件‖Sk+1-Sk‖≤ε時,Sk→sign(L).其中ε為收斂函數,通常設置為10-6.