應(yīng)用球型差分模板的低秩有限差分法純qP 波逆時(shí)偏移

楊禮勝，黃金強(qiáng)*，，高國(guó)超，，夏鵬，，何云川，吳浩

（1. 貴州大學(xué)資源與環(huán)境工程學(xué)院，貴州貴陽(yáng) 550025；2. 喀斯特地質(zhì)資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室（貴州大學(xué)），貴州貴陽(yáng) 550025）

0 引言

地下介質(zhì)普遍具有各向異性特征，其主要表現(xiàn)為地震波傳播速度隨傳播方向變化，因此在實(shí)際資料處理中采用傳統(tǒng)各向同性偏移方法會(huì)產(chǎn)生較大的成像誤差［1］。近年來(lái)，為提高深部?jī)?chǔ)層及復(fù)雜構(gòu)造區(qū)的成像質(zhì)量，進(jìn)一步挖潛剩余油氣儲(chǔ)量，各向異性偏移已成為地震資料常規(guī)處理的重要環(huán)節(jié)。鑒于各向異性彈性波逆時(shí)偏移存在參數(shù)較多且精度有限、波場(chǎng)分離難度大、計(jì)算成本高等難題，在實(shí)際生產(chǎn)中，通常采用更具可行性的qP（quasi-P）波逆時(shí)偏移［2-3］。

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

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

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

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

1 方法原理

1.1 純qP 波波場(chǎng)延拓公式

三維常密度聲波波動(dòng)方程可以表示為

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

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

引入一個(gè)復(fù)數(shù)波場(chǎng)

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

結(jié)合式（5），式（2）可等價(jià)為

對(duì)于變速介質(zhì)，?=v(x)|k|=?(x，k)，變?yōu)榭臻g—波數(shù)域混合形式，此時(shí)復(fù)數(shù)波場(chǎng)P?仍滿足式（6），其時(shí)間—空間域形式為

式中Δt為時(shí)間延拓步長(zhǎng)。通過(guò)Fourier變換，將式（8）轉(zhuǎn)換至波數(shù)域

式（9）即為波數(shù)域地震波場(chǎng)的時(shí)間正、反向一步法延拓公式。式（9）涉及復(fù)數(shù)運(yùn)算，且?(x，k)為空間—波數(shù)域混合形式，實(shí)現(xiàn)過(guò)程較為復(fù)雜，因此本文引入兩步法延拓公式。將式（9）正、反向算子對(duì)應(yīng)相加，可消除復(fù)數(shù)波場(chǎng)的虛部Q，再對(duì)實(shí)部P進(jìn)行反Fourier變換，得到空間域?qū)嵅▓?chǎng)p的兩步法延拓公式

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

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

由上述推導(dǎo)可知，地震波傳播算子與圓頻率函數(shù)ω有關(guān)，因此，可選用不同介質(zhì)的頻散關(guān)系構(gòu)建對(duì)應(yīng)的波場(chǎng)延拓算子。為構(gòu)建VTI 介質(zhì)純qP 波延拓算子，引入VTI介質(zhì)精確頻散關(guān)系

式中：vP和vS分別為沿對(duì)稱軸方向的qP 波、qSV 波相速度；ε和δ為T(mén)homsen 各向異性參數(shù)與各向同性面內(nèi)的波矢量有關(guān)；vH為各向同性面內(nèi)的qP 波相速度，為qP 波動(dòng)校正速度。

擬聲波近似下的耦合型qP 波方程在正演模擬中之所以會(huì)出現(xiàn)偽橫波，是因?yàn)椋孩僭谕茖?dǎo)顯式qP 波波動(dòng)方程過(guò)程中，為了消去頻散關(guān)系中的根號(hào)，防止分?jǐn)?shù)階算子的出現(xiàn)，對(duì)qP 波頻散關(guān)系式進(jìn)行了平方處理，從而使方程產(chǎn)生了增根［25］；②在各向異性介質(zhì)中，令沿對(duì)稱軸方向的橫波速度為零，這一做法并不能使其他方向的橫波相速度和群速度為零，因此在波場(chǎng)模擬時(shí)會(huì)出現(xiàn)偽橫波干擾。

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

再結(jié)合式（10）與式（11）即可得出VTI 介質(zhì)純qP 波波場(chǎng)延拓算子，該算子能夠有效地消除偽橫波干擾，同時(shí)可以準(zhǔn)確地描述VTI 介質(zhì)純qP 波的運(yùn)動(dòng)學(xué)特征。對(duì)該算子的波數(shù)分量作如下坐標(biāo)旋轉(zhuǎn)變換

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

1.2 低秩有限差分算法

首先，依據(jù)式（10）與式（13）構(gòu)建VTI 介質(zhì)純qP波波場(chǎng)傳播矩陣

式中W為N×M階矩陣，方向的波數(shù)個(gè)數(shù)（本文取三個(gè)方向的波數(shù)個(gè)數(shù)相等）。該矩陣的任一元素2cos(ωi，jΔt)表示當(dāng)(x，k)取(xi，kj)時(shí)，采用式（13）計(jì)算的qP 波角頻率?？梢钥闯觯涸摼仃囉上嗨俣葀、各向異性參數(shù)、波矢量k及時(shí)間步長(zhǎng)Δt共同決定。

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

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

對(duì)子矩陣W2作進(jìn)一步分解［23］，有

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

最后，將式（17）代入式（10），根據(jù)Fourier 變換的時(shí)移特性，式（10）可改寫(xiě)為［20］

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

式（20）即為基于低秩有限差分的VTI介質(zhì)純qP波兩步法波場(chǎng)延拓公式，對(duì)上式進(jìn)行移項(xiàng)，可以得到其對(duì)應(yīng)的波場(chǎng)正傳與反傳公式

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

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

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

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

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

1.3 逆時(shí)偏移流程及GPU 算法

逆時(shí)偏移的實(shí)現(xiàn)過(guò)程主要包括震源波場(chǎng)正傳、檢波點(diǎn)波場(chǎng)反傳以及互相關(guān)成像三個(gè)步驟：首先將震源波場(chǎng)沿時(shí)間正方向延拓，并將所有時(shí)刻波場(chǎng)值存入硬盤(pán)；隨后將檢波點(diǎn)波場(chǎng)沿時(shí)間反方向延拓，同時(shí)讀取硬盤(pán)中相應(yīng)時(shí)刻的震源波場(chǎng)值；最后應(yīng)用互相關(guān)成像條件

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

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

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

（1）輸入帶道頭的炮數(shù)據(jù)文件；

（2）搜索觀測(cè)系統(tǒng)及計(jì)算參數(shù)，確定觀測(cè)系統(tǒng)類型，激發(fā)總炮數(shù)、炮間距及炮點(diǎn)坐標(biāo)位置，每炮的接收點(diǎn)個(gè)數(shù)、道間距及接收點(diǎn)坐標(biāo)位置，時(shí)間采樣點(diǎn)數(shù)及采樣間隔；

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

（4）進(jìn)入炮循環(huán)，根據(jù)觀測(cè)系統(tǒng)加載當(dāng)前炮的地震數(shù)據(jù)，并將模型參數(shù)與各向異性參數(shù)讀入內(nèi)存，同時(shí)生成地震子波；

（5）定義差分模板類型（二維或三維）及模板長(zhǎng)度L，并計(jì)算差分模板對(duì)應(yīng)的網(wǎng)格索引數(shù)組再通過(guò)內(nèi)循環(huán)，計(jì)算空間所有位置處的差分系數(shù)G(x，ξl)；

（6）在計(jì)算區(qū)域內(nèi)分配GPU 線程，首先計(jì)算震源波場(chǎng)并保存震源波場(chǎng)的邊界值。具體步驟：根據(jù)差分模板的網(wǎng)格點(diǎn)索引數(shù)組計(jì)算當(dāng)前時(shí)刻的兩波場(chǎng)之和p（xL，t）+p（xR，t）；再將差分系數(shù)與兩波場(chǎng)之和相乘；最后加上震源時(shí)間函數(shù)，即可實(shí)現(xiàn)震源波場(chǎng)正向時(shí)間延拓。該過(guò)程采用GPU 加速，由零時(shí)刻逐步延拓至最大時(shí)刻，還需保存最后時(shí)刻的震源波場(chǎng)；

（7）其次，反傳檢波點(diǎn)波場(chǎng)，與步驟（6）相似，根據(jù)式（22）便可沿逆時(shí)間更新檢波點(diǎn)波場(chǎng)，實(shí)現(xiàn)檢波點(diǎn)波場(chǎng)的反向時(shí)間延拓，該過(guò)程仍然采用GPU 加速，所不同的是，震源加載變?yōu)榈卣饠?shù)據(jù)加載，由最大時(shí)刻反向延拓至零時(shí)刻；同步進(jìn)行的還有震源波場(chǎng)重構(gòu)，通過(guò)加載震源波場(chǎng)的邊界值即可重構(gòu)震源波場(chǎng)，該過(guò)程同樣也利用GPU 加速計(jì)算；

（8）將同一時(shí)刻的檢波點(diǎn)波場(chǎng)與重構(gòu)的震源波場(chǎng)進(jìn)行互相關(guān)計(jì)算便可提取該時(shí)刻的成像值，將所有時(shí)刻的成像值進(jìn)行累加，可得到當(dāng)前炮的單炮成像結(jié)果；

（9）逐炮偏移，并將單炮成像結(jié)果進(jìn)行疊加，直到完成所有炮；

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

2 模型試算

2.1 均勻模型波場(chǎng)快照

本文首先開(kāi)展了三個(gè)均勻VTI 介質(zhì)正演試驗(yàn)來(lái)對(duì)比說(shuō)明不同正演方法的特點(diǎn)，以驗(yàn)證本文方法的正確性。模型參數(shù)如表1所示，模型縱、橫向網(wǎng)格點(diǎn)數(shù)均為301，網(wǎng)格間距為10 m；震源采用主頻為20 Hz 的Ricker子波，在模型中心（1.5 km，1.5 km）處激發(fā)；時(shí)間采樣點(diǎn)數(shù)為1501，時(shí)間步長(zhǎng)為1 ms。

表1 均勻VTI 介質(zhì)模型參數(shù)

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

圖3 三個(gè)均勻模型不同方法模擬的0.3 s 時(shí)刻波場(chǎng)快照

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

2.2 Sigsbee2a 模型

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

圖5 Sigsbee2a 模型

圖6為國(guó)際標(biāo)準(zhǔn)地震數(shù)據(jù)的不同單炮記錄，可以看出：地震記錄中淺部反射波同相軸連續(xù)，且淺部能量強(qiáng)、深部能量弱，與實(shí)際采集的地震記錄特征相同。分別應(yīng)用基于SGFD 和基于LFD 的逆時(shí)偏移方法進(jìn)行成像試處理，成像結(jié)果如圖7所示?？梢?jiàn)：在模型左側(cè)約0～6 km 范圍內(nèi)及鹽丘右上部，因構(gòu)造簡(jiǎn)單，兩種方法都清晰地刻畫(huà)出地層的起伏形態(tài)和薄層厚度，成像效果較好；相比于圖7a的SGFD法成像結(jié)果，在圖7b所示的LFD法成像結(jié)果中，鹽丘內(nèi)部的成像噪聲更小，反射界面的同相軸連續(xù)性更好，鹽丘輪廓的成像分辨率更高，淺中深層的能量更加均衡，因此總體成像質(zhì)量更好。

圖6 Sigsbee2a 模型單炮正演記錄

圖7 Sigsbee2a 模型兩種方法逆時(shí)偏移結(jié)果對(duì)比

在偏移成像測(cè)試過(guò)程中，兩種方法均采用了GPU 并行加速計(jì)算和震源波場(chǎng)重構(gòu)策略，對(duì)其單炮偏移的耗時(shí)和內(nèi)存占用進(jìn)行了對(duì)比，結(jié)果如表2所示。可見(jiàn)，相較于SGFD 法，LFD 方法的計(jì)算耗時(shí)更少、內(nèi)存占用更低。測(cè)試結(jié)果表明，基于LFD 的逆時(shí)偏移成像方法能夠同時(shí)兼顧成像精度和計(jì)算效率。

表2 Sigsbee2a 模型兩種方法單炮逆時(shí)偏移性能對(duì)比

2.3 VTI-Hess 模型

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

圖8 VTI-Hess 模型

第一套數(shù)據(jù)的觀測(cè)系統(tǒng)如圖9a所示，該數(shù)據(jù)的激發(fā)震源是主頻為25 Hz 的Ricker 子波，激發(fā)炮數(shù)是721，第一炮位于x=0 處，炮間距為30.48 m，等間距激發(fā)，采用656道檢波器等間距移動(dòng)接收，接收點(diǎn)的道間距為12.19 m，炮檢距范圍介于0～7984.45 m；從第463、464 炮開(kāi)始變觀采集，接收道數(shù)分別變?yōu)?53、651，隨后依次每?jī)膳诜謩e按3、2 道等差遞減，最后兩炮道數(shù)分別為13、11 道，最后一炮的最大炮檢距為121.92 m，因此總接收道數(shù)為388422；時(shí)間采樣點(diǎn)數(shù)為1332，時(shí)間步長(zhǎng)為6 ms。

圖9 兩套觀測(cè)系統(tǒng)示意圖

第二套數(shù)據(jù)的觀測(cè)系統(tǒng)如圖9b 所示，該數(shù)據(jù)的激發(fā)震源是18 Hz的Ricker子波，激發(fā)炮數(shù)是361 炮，第一炮位于x=0處，炮間距為60.96 m，等間距激發(fā)，采用500 道檢波器等間距移動(dòng)接收，接收點(diǎn)的道間距為12.19 m，炮檢距范圍介于0～6082.81 m；從第263炮開(kāi)始做了變觀采集，接收道數(shù)變?yōu)?98，隨后逐炮按5 道等差遞減，最后一炮道數(shù)變?yōu)?，其最大炮檢距為85.34 m，因此總接收道數(shù)為156047；時(shí)間采樣點(diǎn)數(shù)為2000，時(shí)間步長(zhǎng)為4 ms。

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

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

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

圖11 不含表層多次波時(shí)不同方法VTI-Hess 模型逆時(shí)偏移結(jié)果

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

圖12 含表層多次波時(shí)兩種方法VTI-Hess 模型逆時(shí)偏移結(jié)果對(duì)比

2.4 三維Overthrust 模型

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

圖13 三維Overthrust 速度模型

圖14 為三維Overthrust模型LFD 逆時(shí)偏移成像結(jié)果，由圖可見(jiàn)：本方法對(duì)水平地層及傾斜地層均可以實(shí)現(xiàn)準(zhǔn)確成像，并且能夠清晰刻畫(huà)逆掩斷層的形態(tài)和斷面的實(shí)際位置，對(duì)于復(fù)雜地質(zhì)構(gòu)造具有較高的成像分辨率。

圖14 三維Overthrust 模型（左）與LFD 逆時(shí)偏移結(jié)果（右）的對(duì)比

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

3 結(jié)論與認(rèn)識(shí)

本文推導(dǎo)了各向異性介質(zhì)純qP 波波場(chǎng)延拓公式，隨后引入三維球型差分模板求取與模型相適應(yīng)差分系數(shù)，并將其應(yīng)用于三維純qP 波波場(chǎng)正、反向延拓，最后利用GPU 并行進(jìn)行加速計(jì)算，實(shí)現(xiàn)了基于GPU 加速的三維LFD 純qP 波逆時(shí)偏移成像。由模型測(cè)試及理論分析，得出以下結(jié)論與認(rèn)識(shí)。

（1）本文構(gòu)建的純qP 波延拓公式克服了模型參數(shù)的限制，在任意模型參數(shù)下均能實(shí)現(xiàn)波場(chǎng)穩(wěn)定傳播，并且消除了偽橫波干擾及數(shù)值不穩(wěn)定現(xiàn)象。

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

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

將LFD 應(yīng)用于三維黏聲各向異性介質(zhì)逆時(shí)偏移成像是下一步的研究重點(diǎn)。