應用球型差分模板的低秩有限差分法純qP 波逆時偏移

楊禮勝，黃金強*，，高國超，，夏鵬，，何云川，吳浩

（1. 貴州大學資源與環境工程學院，貴州貴陽 550025；2. 喀斯特地質資源與環境教育部重點實驗室（貴州大學），貴州貴陽 550025）

0 引言

地下介質普遍具有各向異性特征，其主要表現為地震波傳播速度隨傳播方向變化，因此在實際資料處理中采用傳統各向同性偏移方法會產生較大的成像誤差［1］。近年來，為提高深部儲層及復雜構造區的成像質量，進一步挖潛剩余油氣儲量，各向異性偏移已成為地震資料常規處理的重要環節。鑒于各向異性彈性波逆時偏移存在參數較多且精度有限、波場分離難度大、計算成本高等難題，在實際生產中，通常采用更具可行性的qP（quasi-P）波逆時偏移［2-3］。

對于各向異性介質，用于實現其逆時偏移算法的qP 波方程大多源于擬聲波近似。Alkhalifah［2］假設沿對稱軸方向的橫波速度為零，簡化了VTI（Vertical Transversely Isotropic）介質相速度公式，導出了VTI介質四階擬聲波波動方程。在此基礎上，人們利用不同方法發展了多種形式的擬聲波方程［4-5］。然而，基于擬聲波近似思想的qP 波方程一直存在qP 波與qSV（quasi-SV）波耦合的現象，即使在均勻介質的波場快照中，也包含有殘余qSV 波能量，利用該類方程進行逆時偏移可能出現嚴重的偏移假象和數值不穩定，成像效果不佳。為了徹底消除偽橫波干擾、避免不穩定現象，有學者考慮從qP 波和qSV 波耦合頻散關系出發，推導出不含qSV 波的純qP波方程［6-7］；與耦合波動方程相比，純qP 波能夠實現更精確的波場模擬及逆時偏移，但其控制方程中包含復雜的擬微分算子，不易直接采用數值方法進行求解。因此，針對純qP 波方程的高效求解方法仍是當前的研究熱點。

迄今為止，已發展了多種形式的純qP 波方程及相應的求解方法。Chu等［8］利用泰勒系數展開簡化了純qP 波方程中的擬微分算子，推導了TTI（Titled Transversely Isotropic）介質純qP 波波動方程，并給出了偽譜法求解方案，然而該方法在波場延拓的每一時間步需進行多次Fourier 變換（二維7 次，三維22次），計算效率極低［9］。為此，Zhan 等［10］聯合偽譜法與有限差分法求解擬微分算子，減少了計算過程中所需的Fourier 變換次數，在一定程度上提升了偽譜法的計算效率，但是在面對大規?；蛉S復雜模型的數值模擬時，計算成本仍會限制該類方法在實際生產中的應用［11］。此外，Xu 等［12］另辟蹊徑，將偽微分算子分解為一個橢圓微分算子和一個標量算子，并采用有限差分法進行求解，該策略計算效率較高，易于實現，已推廣至各類復雜各向異性介質波場模擬和偏移成像［13］。用有限差分法進行純qP 波波場延拓需求解波場的空間偏導數，較大空間步長時易出現數值頻散現象，與偽譜法相比，計算精度較低。近年來，為平衡純qP 波逆時偏移的計算精度與計算效率，還發展了快速泊松算法［14］、低秩分解算法［15-16］等。

為克服上述純qP 波求解方法的諸多不足，本文將兼具計算效率與計算精度的低秩有限差分法用于純qP 波波場延拓。Fomel等［17］利用低秩分解算法來近似空間—波數混合域地震波傳播算子，實現了各向同性和各向異性介質的波場模擬，并指出低秩分解算法可以準確地描述地震波的運動學特征。隨后，Song等［18］將該方法用于正交各向異性介質的正演模擬；Sun等［19］聯合低秩分解與一步法延拓公式實現了黏彈介質的正演模擬及逆時偏移。在此基礎上，顧漢明等［16］實現了基于低秩一步法延拓的黏聲各向異性純qP波正演模擬。由上述研究成果可知：低秩分解算法是利用頻散關系構建不同介質的波場傳播矩陣，進而實現精確的波場模擬，且無需推導顯式波動方程。因此，該方法在各向異性介質正演模擬及逆時偏移中具有計算穩定、無數值頻散以及無偽橫波干擾等優勢；然而，該方法在波場延拓的每一時間步也需要進行多次Fourier變換（一次正變換和多次反變換，反變換次數與模型復雜度有關）［15］，其計算效率較低。針對這一不足，Song 等［20］提出低秩有限差分算法，即利用低秩分解求取差分系數，有效地降低了波場模擬的計算成本，同時該方法還繼承了低秩分解和有限差分法的優勢，且具有較高的并行性。隨后，方剛等［21］將低秩分解與交錯網格有限差分相結合，提出了交錯網格低秩有限差分法，實現了較大時間步長和空間步長下的波場延拓；袁雨欣等［22］將該方法用于彈性波正演模擬，有效地提高了常規交錯網格有限差分的數值模擬精度；黃金強等［23］通過構建緊致差分模板求解低秩有限差分系數，實現了復雜TI介質純qP 波正演模擬及逆時偏移。然而，目前對低秩有限差分算法的研究還主要集中于二維情形，尤其是二維正演模擬的探討，未見有應用于三維各向異性純qP 波逆時偏移的相關文獻。

針對上述問題，本文借鑒前人利用低秩分解求解二維差分系數的思路［20，23］，構建了一種適用于三維各向異性介質的球型差分模板，隨后結合低秩分解求解差分系數，將低秩有限差分法由二維拓展至三維，并實現了三維低秩有限差分法純qP 波正演模擬和逆時偏移。首先結合波動方程偽解析解和各向異性介質純qP 波頻散關系，推導了各向異性介質純qP 波偽解析解波場延拓公式；隨后，構建了三維球型差分模板，結合低秩分解求取與模型相適應的差分系數，由此實現三維純qP 波波場延拓；針對計算效率問題，利用GPU（Graphics Processing Unit）并行對波場延拓及成像過程進行加速，實現了基于GPU 加速的低秩有限差分法純qP 波逆時偏移；最后，采用典型二維模型和三維模型進行逆時偏移成像測試，驗證了本方法在二維各向異性介質成像方面的優勢以及對于三維復雜模型成像的適用性。

1 方法原理

1.1 純qP 波波場延拓公式

三維常密度聲波波動方程可以表示為

式中：p(x，t)為壓力波場，x=(x，y，z)表示空間坐標；v=v(x)為速度波場；?2為拉普拉斯算子；t為時間。對于均勻介質，利用Fourier 變換將式（1）變換至波數域

式中：k=(k1，k2，k3)為波數矢量，k1、k2、k3分別為x、y、z方向的波數，表示波數域波場，與空間域波場p(x，t)滿足

引入一個復數波場

式中Q為P的Hilbert 變換。在頻率域，Q(k，ω)=為圓頻率。利用頻散關系將其反變換到時間域

結合式（5），式（2）可等價為

對于變速介質，?=v(x)|k|=?(x，k)，變為空間—波數域混合形式，此時復數波場P?仍滿足式（6），其時間—空間域形式為

式中Δt為時間延拓步長。通過Fourier變換，將式（8）轉換至波數域

式（9）即為波數域地震波場的時間正、反向一步法延拓公式。式（9）涉及復數運算，且?(x，k)為空間—波數域混合形式，實現過程較為復雜，因此本文引入兩步法延拓公式。將式（9）正、反向算子對應相加，可消除復數波場的虛部Q，再對實部P進行反Fourier變換，得到空間域實波場p的兩步法延拓公式

若基于上述公式直接求解，則需進行多次傅里葉變換。為此，本文利用非均勻介質頻散關系，忽略高階項，可對?作如下近似

在各向同性介質中，相速度只與空間坐標x有關，即v=v(x)；而在各向異性介質中，地震波傳播速度隨傳播方向而變化，因此相速度v不僅與x有關，還與波數矢量k有關，即v=v(x，k)。

由上述推導可知，地震波傳播算子與圓頻率函數ω有關，因此，可選用不同介質的頻散關系構建對應的波場延拓算子。為構建VTI 介質純qP 波延拓算子，引入VTI介質精確頻散關系

式中：vP和vS分別為沿對稱軸方向的qP 波、qSV 波相速度；ε和δ為Thomsen 各向異性參數與各向同性面內的波矢量有關；vH為各向同性面內的qP 波相速度，為qP 波動校正速度。

擬聲波近似下的耦合型qP 波方程在正演模擬中之所以會出現偽橫波，是因為：①在推導顯式qP 波波動方程過程中，為了消去頻散關系中的根號，防止分數階算子的出現，對qP 波頻散關系式進行了平方處理，從而使方程產生了增根［25］；②在各向異性介質中，令沿對稱軸方向的橫波速度為零，這一做法并不能使其他方向的橫波相速度和群速度為零，因此在波場模擬時會出現偽橫波干擾。

低秩方法利用頻散關系構建波場傳播算子，無需推導顯式的qP 波波動方程，也無需對頻散關系進行平方處理，所以從根本上避免了偽橫波的產生。低秩方法對精確頻散關系和近似頻散關系均有很好的適用性，在橫波速度未知的情況下，可以假設橫波速度為零，采用擬聲波近似下的頻散關系模擬qP 波波場，這也是已有近似方法中最精確的。利用上述兩種頻散關系模擬的qP 波波場特征十分接近［15］，也即假設橫波速度為零對qP 波相速度沒有太大的影響，因此本文利用qP波近似頻散關系來構建純qP波波場傳播算子。假設vS=0，式（12）簡化為式中η=(ε-δ)/(1+2δ)。

再結合式（10）與式（11）即可得出VTI 介質純qP 波波場延拓算子，該算子能夠有效地消除偽橫波干擾，同時可以準確地描述VTI 介質純qP 波的運動學特征。對該算子的波數分量作如下坐標旋轉變換

最后，將式（14）代入式（13），再結合式（10）與式（11）即可得到TTI 介質純qP 波兩步法延拓公式。由式（13）可見：在時間方向延拓一次需進行一次Fourier 變換和N（N=Nx×Ny×Nz，Nx、Ny、Nz分別為三個方向的模型網格點數）次Fourier反變換，計算復雜度為O(N2)，用于復雜模型波場正演時會造成極高的計算成本。對此，本文引入低秩有限差分法，并結合兩步法波場延拓公式，推導出基于低秩有限差分的波場延拓公式，以降低計算成本。

1.2 低秩有限差分算法

首先，依據式（10）與式（13）構建VTI 介質純qP波波場傳播矩陣

式中W為N×M階矩陣，方向的波數個數（本文取三個方向的波數個數相等）。該矩陣的任一元素2cos(ωi，jΔt)表示當(x，k)取(xi，kj)時，采用式（13）計算的qP 波角頻率。可以看出：該矩陣由相速度v、各向異性參數、波矢量k及時間步長Δt共同決定。

利用低秩分解將該空間—波數域傳播矩陣W近似分解為三個子矩陣的乘積［17］

式中：W1是由全矩陣W中U個線性無關的列向量組成，其中U個列向量分別由{km}m=1，2，…，U計算得到；W2是由全矩陣W中V個線性無關的行向量組成，其中V個行向量分別由{xn}n=1，2，…，V計算得到；{km}、{xn}分別為k、x的子集；A={amn}為連接兩個子矩陣W1、W2的系數矩陣。

對子矩陣W2作進一步分解［23］，有

式中：L表示模板長度，即模板中參與計算的網格點總數；C(xn，ξl)=W2(xn，k)·B?(ξl，k)表示系數矩陣，B?等于差分模板矩陣B的廣義逆矩陣，矩陣B的具體形式為

最后，將式（17）代入式（10），根據Fourier 變換的時移特性，式（10）可改寫為［20］

式中：gl為G的第l列元素；

式（20）即為基于低秩有限差分的VTI介質純qP波兩步法波場延拓公式，對上式進行移項，可以得到其對應的波場正傳與反傳公式

對比式（10）與式（22）可知：①低秩有限差分法在波場延拓時能夠有效減少Fourier 逆變換的次數，在時間方向延拓一次的計算復雜度為O(LN)，計算成本顯著降低；②利用低秩分解法即可求取差分系數，實現波場的穩定延拓，無需計算各方向混合偏導數，有效提升了計算效率；③低秩有限差分延拓公式可以自由設定波矢量的取值范圍，保證地震波場的模擬精度。

本文借鑒二維低秩有限差分模板思想［20，23］，構建了適用于三維各向異性介質的球型差分模板，將前人利用低秩分解求解二維差分系數的思路拓展至球型差分模板的差分系數求解中，最終將低秩有限差分法由二維拓展至三維。本文構建的球型差分模板除了采用坐標軸方向的網格點參與波場遞推計算外，在坐標軸以外的相應網格點也參與計算，與傳統三維交錯網格差分模板相比更具一般性和靈活性。假設模板中心網格點的坐標為(i，j，h)，其余參與計算的網格點坐標(i′，j′，h′)滿足如下關系即所有參與計算的網格點都分布在以網格中心點為球心、以R為半徑的球內，因此，一旦選定R，便可確定差分模板類型、模板長度L及對應參與計算的網格點索引數組

圖1展示了不同階數的二維及三維差分模板示意圖，圖1a中黑色圓圈所標識的字母分別表示該位置處的差分系數G(x，ξl)。以圖1a 左為例，圖中字母a 處的網格點索引為依次類推，字母b 處同一字母（區分大小寫）表示的差分系數相同，觀察可見：差分系數關于原點a 位置呈中心對稱分布；而大寫字母與小寫字母符合軸對稱分布（如：D 與d），在各向同性介質中，兩者相等，在各向異性介質中，兩者不相等。由此可知：二維差分模板表現出中心對稱性，在計算過程中無需存儲差分模板中所有網格點處的差分系數。同理，三維差分模板也表現出中心對稱規律，在計算過程中無需存儲差分模板中所有網格點處的差分系數，與基于傳統Taylor 展開的差分系數不同，Taylor差分系數與空間位置無關，而低秩有限差分需要存儲所有空間位置處的差分模板所對應的差分系數，因此利用該特性能夠成倍降低差分系數的存儲量。低秩有限差分系數是基于低秩分解計算的一種優化差分系數，對于不同介質的不同空間位置，需分別求取對應的差分系數，因此該方法不受模型參數限制。

圖1 二維（a）和三維（b）不同階次低秩有限差分模板示意圖

圖2 基于低秩有限差分的逆時偏移成像流程

1.3 逆時偏移流程及GPU 算法

逆時偏移的實現過程主要包括震源波場正傳、檢波點波場反傳以及互相關成像三個步驟：首先將震源波場沿時間正方向延拓，并將所有時刻波場值存入硬盤；隨后將檢波點波場沿時間反方向延拓，同時讀取硬盤中相應時刻的震源波場值；最后應用互相關成像條件

即可提取逆時偏移成像剖面。式中：I(x，y，z)為最終成像結果；pf(x，y，z，t) 表示正傳震源波場；pb(x，y，z，t)表示反傳檢波點波場。

基于GPU 并行加速的低秩有限差分逆時偏移實現流程如圖2 所示，對于波場延拓與互相關成像部分利用GPU 進行加速計算（圖2 中紅色虛線框所示），即分別利用GPU 加速計算式（22）和式（24）；并且在波場正傳過程中利用自相關條件獲取照明能量，以實現對成像結果的照明補償。本文還采用了基于GPU 并行的震源波場重構策略［26］，即在震源波場正向延拓后，僅存儲每一時刻的邊界波場值，在檢波點波場反向延拓的同時，通過邊界波場值實現震源波場重構，利用同一時刻的兩個波場值即可實現互相關成像。這一策略能夠有效地降低逆時偏移過程中的存儲要求和I/O 傳輸次數，并且增加的計算量也是可以接受的。

基于GPU 并行的低秩有限差分逆時偏移具體算法過程如下：

（1）輸入帶道頭的炮數據文件；

（2）搜索觀測系統及計算參數，確定觀測系統類型，激發總炮數、炮間距及炮點坐標位置，每炮的接收點個數、道間距及接收點坐標位置，時間采樣點數及采樣間隔；

（3）讀入速度文件及各向異性參數文件；

（4）進入炮循環，根據觀測系統加載當前炮的地震數據，并將模型參數與各向異性參數讀入內存，同時生成地震子波；

（5）定義差分模板類型（二維或三維）及模板長度L，并計算差分模板對應的網格索引數組再通過內循環，計算空間所有位置處的差分系數G(x，ξl)；

（6）在計算區域內分配GPU 線程，首先計算震源波場并保存震源波場的邊界值。具體步驟：根據差分模板的網格點索引數組計算當前時刻的兩波場之和p（xL，t）+p（xR，t）；再將差分系數與兩波場之和相乘；最后加上震源時間函數，即可實現震源波場正向時間延拓。該過程采用GPU 加速，由零時刻逐步延拓至最大時刻，還需保存最后時刻的震源波場；

（7）其次，反傳檢波點波場，與步驟（6）相似，根據式（22）便可沿逆時間更新檢波點波場，實現檢波點波場的反向時間延拓，該過程仍然采用GPU 加速，所不同的是，震源加載變為地震數據加載，由最大時刻反向延拓至零時刻；同步進行的還有震源波場重構，通過加載震源波場的邊界值即可重構震源波場，該過程同樣也利用GPU 加速計算；

（8）將同一時刻的檢波點波場與重構的震源波場進行互相關計算便可提取該時刻的成像值，將所有時刻的成像值進行累加，可得到當前炮的單炮成像結果；

（9）逐炮偏移，并將單炮成像結果進行疊加，直到完成所有炮；

（10）對所有炮的成像結果進行拉普拉斯濾波等處理后，便可輸出最終成像剖面。

2 模型試算

2.1 均勻模型波場快照

本文首先開展了三個均勻VTI 介質正演試驗來對比說明不同正演方法的特點，以驗證本文方法的正確性。模型參數如表1所示，模型縱、橫向網格點數均為301，網格間距為10 m；震源采用主頻為20 Hz 的Ricker子波，在模型中心（1.5 km，1.5 km）處激發；時間采樣點數為1501，時間步長為1 ms。

表1 均勻VTI 介質模型參數

圖3a 為常規擬聲波方程交錯網格高階有限差分（SGFD，Staggered Grid Finite Difference）法的正演結果，圖3b 為純qP 波低秩有限差分（LFD，Low-rank Finite Difference）法的正演結果。對比可見：常規擬聲波方程只有在ε≥δ時才能實現波場的穩定傳播（圖3a左、中所示）。需要說明的是，當ε＞δ時，波場中會出現菱形的偽橫波假象，在成像中會形成偏移噪聲和假象，影響成像結果的同相軸連續性和降低剖面信噪比；當ε＝δ時，代表橢圓各向異性介質，此時波場中無菱形干擾；當ε<δ時會出現嚴重的數值不穩定現象（圖3a右所示），在此情形下該方法失效，而純qP 波法則突破了這一參數限制，在任意類型的各向異性介質中均能實現波場的穩定傳播，且正演結果都不存在偽橫波干擾及數值不穩定性。

圖3 三個均勻模型不同方法模擬的0.3 s 時刻波場快照

圖4 展示了圖3a 左和圖3b 左在x=1500 m 處的波場記錄，對比可見：除偽橫波干擾外，擬聲波SGFD法、純qP 波LFD 法與qP 波的解析解高度吻合，驗證了本文方法的正確性。

2.2 Sigsbee2a 模型

為了驗證基于LFD 的逆時偏移成像算法在計算效率上的優勢，選用如圖5所示的Sigsbee2a各向同性模型進行測試。該模型主要由一個高速鹽丘和幾組正、逆斷層構成，模型參數如下：模型網格數為3201×1201，網格間距均為7.62 m，模型尺寸為24.38 km×9.14 km。采用主頻為25 Hz 的Ricker 子波作為激發震源，共激發500 炮，第一炮位于x=281.94 m，炮間距為45.72 m，炮點深度均為7.62 m；采用348道檢波器等間距移動接收，接收點的道間距為22.86 m，接收點深度為7.62 m，炮檢距范圍為0～7932.42 m。從第355 炮開始采用變觀采集，接收道數變為347 道，隨后逐炮按2 道等差遞減，最后一炮的接收道數變為57 道，最后一炮的最大炮檢距為1280.16 m，因此接收總道數為152684；時間采樣點數均為1500，時間步長為8 ms。

圖5 Sigsbee2a 模型

圖6為國際標準地震數據的不同單炮記錄，可以看出：地震記錄中淺部反射波同相軸連續，且淺部能量強、深部能量弱，與實際采集的地震記錄特征相同。分別應用基于SGFD 和基于LFD 的逆時偏移方法進行成像試處理，成像結果如圖7所示?？梢姡涸谀Ｐ妥髠燃s0～6 km 范圍內及鹽丘右上部，因構造簡單，兩種方法都清晰地刻畫出地層的起伏形態和薄層厚度，成像效果較好；相比于圖7a的SGFD法成像結果，在圖7b所示的LFD法成像結果中，鹽丘內部的成像噪聲更小，反射界面的同相軸連續性更好，鹽丘輪廓的成像分辨率更高，淺中深層的能量更加均衡，因此總體成像質量更好。

圖6 Sigsbee2a 模型單炮正演記錄

圖7 Sigsbee2a 模型兩種方法逆時偏移結果對比

在偏移成像測試過程中，兩種方法均采用了GPU 并行加速計算和震源波場重構策略，對其單炮偏移的耗時和內存占用進行了對比，結果如表2所示?？梢姡噍^于SGFD 法，LFD 方法的計算耗時更少、內存占用更低。測試結果表明，基于LFD 的逆時偏移成像方法能夠同時兼顧成像精度和計算效率。

表2 Sigsbee2a 模型兩種方法單炮逆時偏移性能對比

2.3 VTI-Hess 模型

在驗證了本方法正確性和有效性的基礎上，采用VTI-Hess 國際標準模型進一步驗證LFD 逆時偏移成像方法對復雜各向異性介質的適應性。VTI-Hess模型如圖8所示，模型主要由三部分組成，即左側的高速巖體、中部的強各向異性體和右側的陡傾角斷層。模型網格數為3617×1500，網格間距為6.1 m，模型尺寸為22.04 km×9.14 km?？紤]到實際地震數據中常伴有強烈的多次波，因此本文采用該模型開源的兩套地震數據分別進行成像試算，其中第一套數據不含多次波，而第二套數據含有多次波。

圖8 VTI-Hess 模型

第一套數據的觀測系統如圖9a所示，該數據的激發震源是主頻為25 Hz 的Ricker 子波，激發炮數是721，第一炮位于x=0 處，炮間距為30.48 m，等間距激發，采用656道檢波器等間距移動接收，接收點的道間距為12.19 m，炮檢距范圍介于0～7984.45 m；從第463、464 炮開始變觀采集，接收道數分別變為653、651，隨后依次每兩炮分別按3、2 道等差遞減，最后兩炮道數分別為13、11 道，最后一炮的最大炮檢距為121.92 m，因此總接收道數為388422；時間采樣點數為1332，時間步長為6 ms。

圖9 兩套觀測系統示意圖

第二套數據的觀測系統如圖9b 所示，該數據的激發震源是18 Hz的Ricker子波，激發炮數是361 炮，第一炮位于x=0處，炮間距為60.96 m，等間距激發，采用500 道檢波器等間距移動接收，接收點的道間距為12.19 m，炮檢距范圍介于0～6082.81 m；從第263炮開始做了變觀采集，接收道數變為498，隨后逐炮按5 道等差遞減，最后一炮道數變為8，其最大炮檢距為85.34 m，因此總接收道數為156047；時間采樣點數為2000，時間步長為4 ms。

圖10 展示了不含多次波和含有多次波的單炮地震記錄。圖中可見：當不含多次波時，深部反射信號的能量較強；而當含有多次波時，深部反射信號被多次波干擾，很難觀察到明顯的反射波同相軸，此外，由于在頂部未使用吸收邊界條件，表面多次波信息較豐富。

圖10 VTI-Hess 模型單炮地震記錄

當不含多次波時，不同方法的偏移成像結果如圖11 所示，由圖可知：①在VTI 介質中，與基于SGFD逆時偏移成像結果相比，在LFD 成像結果中，左側高速巖體輪廓更加清晰，巖體下部地層成像能量更強，對上部地層及中部各向異性體的分辨率更高；②對比圖11b 與圖11c 可知，用各向同性LFD 進行逆時偏移成像時，地層界面無法聚焦，分辨率降低，右側陡傾斷層成像效果變差（圖11c 中黑色虛線框所示）；③理論情況下采用四周全接收的LFD 逆時偏移結果，在真正意義上實現了波場的完全重構，此時成像剖面的分辨率顯著提高，同相軸明顯變細，能量也十分均衡，但在濾波處理時仍會出現部分偏移噪聲（如圖11d 黑色箭頭所指）。對比上述成像結果可見：LFD 逆時偏移對于高速巖體及各向異性體的成像效果與采用四周接收的LFD 逆時偏移成像效果最接近，能夠獲得高質量的VTI 介質成像結果，這表明該方法能夠應用于復雜各向異性介質的逆時偏移成像處理。

圖11 不含表層多次波時不同方法VTI-Hess 模型逆時偏移結果

當地震數據含有多次波時，VTI-Hess 模型的SGFD、LFD 逆時偏移成像結果分別如圖12所示。圖中可見：由于存在多次波干擾，成像結果中均含有較多的成像噪聲及偏移假象（圖12中黑色虛線框所示）；與圖12a 中基于SGFD 的逆時偏移成像結果相比，基于LFD 的成像結果信噪比更高，對上部地層及中部強各向異性體具有更高的分辨率，總體成像質量更好。由此進一步驗證了本文方法對于各向異性介質逆時偏移成像的適應性。此外，SGFD 和LFD 逆時偏移單炮成像耗時分別為330.17、196.75 s，可見基于LFD 的逆時偏移成像方法對于各向異性介質具有明顯的效率優勢。

圖12 含表層多次波時兩種方法VTI-Hess 模型逆時偏移結果對比

2.4 三維Overthrust 模型

為進一步驗證基于LFD 的三維各向異性介質逆時偏移方法的適用性，本文選用三維Overthrust 模型進行測試。圖13 為速度模型，該模型結構復雜，橫向速度變化劇烈，含有一個較大的逆沖推覆構造，介質非均質性強，地層厚度小，且界面起伏大，是檢驗成像算法優劣的國際標準模型。但該模型缺少各向異性參數，因此本文設置各向異性參數均勻不變，分別為ε=0.24、δ=0.10、傾角θ=30°、方位角φ=45°。模型網格點數為501×301×187，三個方向的網格間距均為25 m；采用主頻為9 Hz 的Ricker 子波作為激發震源，在x和y方向分別激發51、21 炮，共激發1071 炮，炮間距分別為150、250 m，第一炮位于（2500 m，1250 m）處，激發點深度均為125 m；單炮在x方向和y方向分別有201、101 道接收，道間距為25 m，接收點深度為100 m；時間采樣點數1601，時間步長為1.5 ms。

圖13 三維Overthrust 速度模型

圖14 為三維Overthrust模型LFD 逆時偏移成像結果，由圖可見：本方法對水平地層及傾斜地層均可以實現準確成像，并且能夠清晰刻畫逆掩斷層的形態和斷面的實際位置，對于復雜地質構造具有較高的成像分辨率。

圖14 三維Overthrust 模型（左）與LFD 逆時偏移結果（右）的對比

本文采用基于GPU 加速和單CPU 兩套程序進行三維逆時偏移成像測試。測試過程中GPU 型號為NVIDIA Quadro K6000，顯存為12 GB；CPU 型號為Intel Xeon E5-2680 v4 @ 2.40 GHz，內存為58 GB。三維Overthrust模型逆時偏移過程中，GPU 單炮偏移耗時59.69 s，而CPU 單炮偏移耗時2578.28 s；完成1071炮偏移，GPU 總耗時約為1065 min，CPU 總耗時約為46022 min?？梢?，與基于單CPU 串行的LFD 逆時偏移方法相比，基于GPU 并行的LFD 逆時偏移方法能夠有效減少成像的耗時，顯著提高成像效率，兩種方法的用時之比約為1∶43。

3 結論與認識

本文推導了各向異性介質純qP 波波場延拓公式，隨后引入三維球型差分模板求取與模型相適應差分系數，并將其應用于三維純qP 波波場正、反向延拓，最后利用GPU 并行進行加速計算，實現了基于GPU 加速的三維LFD 純qP 波逆時偏移成像。由模型測試及理論分析，得出以下結論與認識。

（1）本文構建的純qP 波延拓公式克服了模型參數的限制，在任意模型參數下均能實現波場穩定傳播，并且消除了偽橫波干擾及數值不穩定現象。

（2）與SGFD 逆時偏移方法相比，LFD 法對復雜各向異性介質具有更高的成像分辨率，當模型含有多次波時，仍然能夠取得較好的成像效果；并且在保證成像精度的前提下，該方法能夠有效降低逆時偏移計算耗時及內存占用，表明LFD 法在各向異性介質逆時偏移成像中具有一定的優勢。

（3）基于LFD 的逆時偏移方法對于三維各向異性介質的逆時偏移成像同樣具有較好的適用性，且本文采用的GPU 并行加速策略能夠有效提高三維逆時偏移計算效率，其加速比可達43 倍左右，有利于推廣和應用于實際地震數據。

將LFD 應用于三維黏聲各向異性介質逆時偏移成像是下一步的研究重點。