VTI介質(zhì)彈性體波模式解耦高效算法

阮 倫 程玖兵

(①同濟(jì)大學(xué)海洋地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092； ②同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院，上海 200092)

0 引言

地震體波在地球內(nèi)部非均勻介質(zhì)中的傳播過程中存在縱波(P波)和橫波(S波)兩種模式(即波型)之間的相互轉(zhuǎn)換與耦合。為了獲得具有明確物理含義的彈性波偏移圖像，縱橫波模式解耦被認(rèn)為是至關(guān)重要的步驟[1-4]。模式解耦預(yù)條件對(duì)于彈性全波形反演也很有意義[5-6]。

在各向同性介質(zhì)中，根據(jù)亥姆赫茲分解理論，彈性矢量波場(chǎng)可以分解成相應(yīng)的縱波和橫波成分[7-8]，馬德堂等[9]、李振春等[10]提出了P波和S波近似解耦的彈性波方程；吳瀟等[11]比較了一系列各向同性P波和S波解耦方法。然而，地球內(nèi)部介質(zhì)普遍存在不同程度的彈性各向異性，表現(xiàn)為傳播速度與偏振特征的方向相關(guān)性。此時(shí)，縱波的偏振方向與傳播方向不平行，橫波偏振方向與傳播方向不垂直，而橫波還會(huì)分裂成偏振方向相互正交的快橫波和慢橫波[12]，它們依次被稱為qP波、qS1波和qS2波。qS1和qS2波的偏振方向隨著傳播方向的變化通常是不連續(xù)的，而且當(dāng)兩種qS波相速度相等(即發(fā)生橫波奇異性)時(shí)，無法確定各自的偏振方向，故而難以區(qū)分兩種qS波模式。在對(duì)稱性較高的橫向各向同性(TI)介質(zhì)中，qS1和qS2的偏振方向要么在對(duì)稱軸平面以內(nèi)(qSV波)，要么與對(duì)稱軸平面垂直(qSH波)[13-16]。

各向異性彈性波模式解耦理論和算法研究有著近三十年的發(fā)展歷程。一方面，人們一直試圖從方程解耦的角度實(shí)現(xiàn)qP和qS波場(chǎng)的解耦，代表性的工作包括擬聲學(xué)近似[17]、純模式近似[18]和偽純模式波動(dòng)方程[19-20]等，這些近似方程都存在動(dòng)力學(xué)誤差。另一方面，受亥姆赫茲分解理論的啟發(fā)，人們也對(duì)各向同性介質(zhì)中縱、橫波分離的算法進(jìn)行了擴(kuò)展。最初，Dellinger等[21]提出在波數(shù)域求取偽散度和偽旋度算子，通過偏振投影分離qP與qS波。為了適應(yīng)非均勻各向異性介質(zhì)，Yan等[14]將偏振投影算子構(gòu)建成非穩(wěn)相空間濾波形式，魏石磊等[22]在此基礎(chǔ)上研究了空間域?yàn)V波算子的特征。然而非均勻各向異性介質(zhì)中基于偏振投影的模式解耦計(jì)算成本很高[16]，為了突破計(jì)算瓶頸，Yan等[23]設(shè)計(jì)了一種類似于單程波偏移成像中的相移加插值(PSPI)的混合域解耦算法。而后，Cheng等[16]將模式分解表達(dá)成空間—波數(shù)域的廣義傅里葉積分算子，并借用矩陣低秩近似技術(shù)[24]，提出了一種分離qP、qSV和qSH波標(biāo)量和矢量場(chǎng)的快速算法；Wang 等[25]建議對(duì)介質(zhì)模型進(jìn)行分塊，再利用低秩近似解耦算法。此外，也有學(xué)者將Yan等[14]提出的非穩(wěn)相空間濾波算子簡(jiǎn)化成更經(jīng)濟(jì)的偽散度/偽旋度算子，實(shí)現(xiàn)qP和qS波近似解耦[26]?？傊瑥木?、成本和可操作性等角度發(fā)展各向異性彈性波模式解耦新算法仍然非常必要。

本文在回顧各向異性彈性波模式解耦的基本原理基礎(chǔ)上，首先改進(jìn)了Zhou等[27]提出的基于主能量傳播與偏振方向偏差角構(gòu)造的qP/qS波解耦偽散度、偽旋度算法，一定程度上考慮偏差角隨傳播方向的變化，進(jìn)而建立精度更高的偽散度、偽旋度算子。然后，探討低秩近似解耦算法及其模型分區(qū)、CPU-GPU異構(gòu)平臺(tái)實(shí)現(xiàn)策略。最后，應(yīng)用SEG Hess VTI模型檢驗(yàn)和評(píng)價(jià)本文的模式解耦高效算法。

1 彈性波模式解耦基本原理

在各向同性介質(zhì)中，彈性波場(chǎng)P/S波分離通常采用散度和旋度運(yùn)算[8]

(1)

(2)

式中：K=(Kx,Ky,Kz)T表示波矢量；U或U為u或u的傅里葉變換。在各向同性介質(zhì)中，波的傳播方向與P波的偏振方向一致，因此上面的模式解耦也可視為偏振投影的范疇。

在各向異性介質(zhì)中，采用散度和旋度運(yùn)算會(huì)造成模式分離不徹底，出現(xiàn)模式泄漏或串?dāng)_。為此，Dellinger等[21]提出以qP波偏振方向?yàn)榛A(chǔ)的偽散度和偽旋度算子，從理論上解決了各向異性qP與qS波分離問題，即

(3)

式中aP(K)為歸一化的qP波偏振矢量。

為了求取偏振矢量，需要求解平面波解對(duì)應(yīng)的Christoffel方程

(4)

式中：m為qP、qS1、qS2；I為單位對(duì)角矩陣；G為Christoffel張量

Gik=Cijklnjnli,j,k,l=1,2,3

(5)

在VTI介質(zhì)中，qS波一般分解成偏振方向連續(xù)變化的qSV和qSH波，前者在對(duì)稱軸平面內(nèi)垂直qP波偏振，后者在對(duì)稱軸平面法向方向上偏振[13]。因此，TI介質(zhì)彈性波模式解耦對(duì)應(yīng)波數(shù)域的偏振投影公式為

U(m)(K)=iam(K)·U(K)

(6)

式中qSV、qSH波歸一化偏振矢量滿足

(7)

式中s為TI介質(zhì)的對(duì)稱軸方向。

不難發(fā)現(xiàn)，如同散度和旋度運(yùn)算一樣，Dellinger等[21]構(gòu)造的偽散度、偽旋度運(yùn)算分離得到標(biāo)量縱波和矢量橫波場(chǎng)(其中橫波場(chǎng)在二維情況下僅剩一個(gè)分量)。而式(6)表示偏振投影解耦得到了三種波模式的標(biāo)量場(chǎng)。為了從總的彈性波場(chǎng)分離出各個(gè)波模式量綱一致、振幅和相位相對(duì)保真的矢量場(chǎng)，可采用Zhang 等[15]提出的矢量分解算子

U(m)(K)=am(K)[am(K)·U(K)]

(8)

對(duì)于非均勻VTI介質(zhì)，偏振投影算子比較復(fù)雜，因?yàn)槠穹较蛟诳臻g上隨介質(zhì)參數(shù)變化而改變。例如，Yan等[14]把式(6)改寫為空間域的非穩(wěn)相濾波(投影)算子

(9)

2 偏振投影空間域?yàn)V波近似算法

在各向異性介質(zhì)中，qP波傳播方向與偏振方向通常是不一致的，它們之間存在一個(gè)偏差角。一旦給定傳播方向，可以借助Christoffel方程計(jì)算精確的偏振方向。為了避免采用計(jì)算復(fù)雜度非常高的精確的偏振投影算子，有學(xué)者利用相對(duì)容易估算的相(傳播)方向或群方向去估計(jì)傳播與偏振方向之間的夾角，進(jìn)而獲得偏振方向的近似估計(jì)，從而構(gòu)造出比較容易計(jì)算的、近似的qP/qS波模式分離算子[27]。

三維情況下，qP波偏振矢量aP(K)的方向可用方位角α0(K)和傾角β0(K)表示，波矢量K的方向可用方位角α1(K)和傾角β1(K)表示，因此，可以用φ(K)和θ(K)表示aP(K)與K的方向偏差角，即φ(K)=α0(K)-α1(K)，θ(K)=β0(K)-β1(K)。理論上，除了極個(gè)別特殊方向(如對(duì)稱軸方向)，這兩個(gè)偏差角隨傳播方向變化，且都不為零。利用偏差角可從傳播方向估計(jì)qP波的偏振方向，即

(10)

將式(10)施加到偏振投影算子并變換到空間域，就與式(9)完全等價(jià)。對(duì)非均勻TI介質(zhì)而言，空間域非穩(wěn)相濾波算子遠(yuǎn)比散度/旋度算子復(fù)雜，計(jì)算成本極高。

2.1 基于主能量傳播與偏振方向夾角的偽散度與偽旋度算子

為了在精確的偏振投影非穩(wěn)相濾波算子和散度/旋度算子之間找到比較恰當(dāng)?shù)恼壑校琙hou等[27]提出了qP波偏振方向的近似估算方法，進(jìn)而構(gòu)造出了一種非常經(jīng)濟(jì)的偽散度與偽旋度算子。其思路是，首先用散度運(yùn)算獲得近似qP波標(biāo)量場(chǎng)u(P)，然后估計(jì)其能流密度矢量(坡印廷矢量)e=(e1，e2，e3)T

(11)

表征空間各點(diǎn)主能量傳播方向K。針對(duì)一階速度—應(yīng)力方程計(jì)算的矢量波場(chǎng)，也可采用

(12)

計(jì)算的群方向去估計(jì)相方向，從而避免估計(jì)qP波標(biāo)量場(chǎng)引起的誤差。式中：τjk表示應(yīng)力分量；vk表示振動(dòng)速度分量；j，k=1，2，3。再由Christoffel方程計(jì)算主能量方向qP波的偏振方向以及它與傳播方向的夾角。為了簡(jiǎn)化投影濾波算子，相同空間點(diǎn)上可采用不隨傳播方向變化的偏差角去逼近qP波的偏振方向[27]，即

(13)

上述近似雖可以由傳播方向估計(jì)qP波的偏振方向，但無法估計(jì)qSV和qSH波的偏振方向，至多只能分離qP與qS波。另外，當(dāng)qP、qS波從不同方向同時(shí)達(dá)到該點(diǎn)時(shí)，忽視偏差角隨傳播方向變化也會(huì)引入明顯誤差。

根據(jù)式(13)，偏振投影(濾波)算子可表示成旋轉(zhuǎn)后的散度/旋度算子(即偽散度/偽旋度算子)，即

(14)

式中

(15)

與旋度算子類似，偽旋度算子僅從矢量總波場(chǎng)中分離出矢量形式的qS波，無法區(qū)分qSV和qSH波。

在Zhou等[27]僅推導(dǎo)了二維情況下的偽散度/偽旋度算子

(16)

以及相應(yīng)的模式解耦非穩(wěn)相濾波算子

(17)

本質(zhì)上，這些近似的偏振投影運(yùn)算用傳播與偏振方向夾角相關(guān)的系數(shù)去顯式地對(duì)波場(chǎng)的一階偏導(dǎo)數(shù)進(jìn)行加權(quán)求和。由于這些近似的權(quán)系數(shù)可以提前計(jì)算，偏導(dǎo)數(shù)波場(chǎng)在傳播模擬過程中也預(yù)先計(jì)算好，這種空間非穩(wěn)相濾波的計(jì)算復(fù)雜度大幅度降低。值得注意的是，由于傳播方向是從每個(gè)時(shí)刻波場(chǎng)的坡印廷矢量估計(jì)的，因此每個(gè)延拓時(shí)間步上都需要構(gòu)造并施加相應(yīng)的偽散度/偽旋度算子。

如圖1所示，在一個(gè)兩層VTI介質(zhì)模型中，震源位于(2000m，0)，通過彈性波方程數(shù)值模擬計(jì)算得到900ms時(shí)刻的質(zhì)點(diǎn)速度波場(chǎng)。然后采用主能量方向偏振投影近似投影算法，基本實(shí)現(xiàn)大部分空間qP、qSV波的分離。只是在一些波前面交叉區(qū)域(圖1d方框所示)，存在明顯的模式泄漏。其原因在于，傳播與偏振矢量夾角隨傳播方向變化，僅依據(jù)主能量方向估計(jì)的一個(gè)偏差角構(gòu)建偏振投影算子，誤差較大，不能處理不同方向同時(shí)達(dá)到多種波模式的分離問題。其中采用主能量方向偏振投影近似算法得到的是標(biāo)量qP和qSV波。同時(shí)，結(jié)合式(1)和式(2)，波場(chǎng)空間偏導(dǎo)數(shù)運(yùn)算相當(dāng)于在波數(shù)域乘以iK，虛數(shù)單位i會(huì)使得縱、橫波相位改變90°，與非歸一化的波數(shù)K相乘會(huì)使二者的振幅發(fā)生改變(對(duì)比圖1b和圖1d)。為了保證分離得到的縱、橫波在相位和振幅上盡可能保真，可以采用以式(8)為基礎(chǔ)的矢量分解模式解耦算法，圖1f～圖1k為精確偏振投影算法分離的矢量qP和qSV波，可以看出矢量qP和qSV波的水平分量和垂直分量相位上分別與質(zhì)點(diǎn)速度的水平分量和垂直分量保持一致。

圖1 兩層VTI介質(zhì)模型偏振投影近似算法分離結(jié)果

2.2 基于多個(gè)主能量傳播與偏振方向夾角的偽散度與偽旋度算子

一般來講，波傳播到某一點(diǎn)時(shí)能量隨球面擴(kuò)散，傳播過程中qP波與qS波在同一時(shí)刻都經(jīng)過空間某點(diǎn)也是非常普遍的。為了提升Zhou等[27]方法的精度，本文適當(dāng)考慮傳播與偏振方向偏差角隨傳播方向的變化，通過多方向子波場(chǎng)的偽散度與偽旋度運(yùn)算實(shí)現(xiàn)更精確的模式解耦。

首先，將原始波場(chǎng)各個(gè)分量變換到波數(shù)域，按角度面元拆分成若干方向扇區(qū)(Ωi)的子波場(chǎng)[28]

(18)

式中u(x)和U(K)分別表示空間域、波數(shù)域的子波場(chǎng)。N一般取4～8，具體取值視波場(chǎng)復(fù)雜性和精度要求而定。然后，分別在每個(gè)方向扇區(qū)計(jì)算相應(yīng)的坡印廷矢量ei，指示多個(gè)主能量傳播方向。再利用Christoffel方程估計(jì)不同帶區(qū)的偏差角φi和θi，進(jìn)而建立各個(gè)方向扇區(qū)內(nèi)偏振方向的估算公式

(19)

據(jù)此構(gòu)造出每個(gè)方向扇區(qū)內(nèi)的偏振投影算子，將其作用于相應(yīng)帶區(qū)的子波場(chǎng)就獲得解耦后的qP與qS波子波場(chǎng)。最后所有方向帶區(qū)解耦子波場(chǎng)疊加就得到最終分離的qP和qS波場(chǎng)。

如圖2所示，對(duì)8個(gè)主能量方向子波場(chǎng)施加相應(yīng)的偏振投影再求和，明顯提升了qP與qS波重疊區(qū)域的模式分離精度。與精確偏振投影算法處理結(jié)果(圖1f和圖1g)相比，整體上效果是可接受的。不足的是，基于多個(gè)主能量傳播與偏振方向夾角的偽散度與偽旋度算子分離的標(biāo)量qP和qSV波在相位和振幅上仍存在一定的誤差(對(duì)比圖1f和圖1g與圖2a和圖2b)。

圖2 兩層VTI介質(zhì)模型改進(jìn)的偏振投影近似算法分離結(jié)果

3 偏振投影低秩近似算法及其模型分區(qū)與GPU加速策略

在非均勻VTI介質(zhì)中，彈性波模式解耦的偏振投影對(duì)應(yīng)空間域非穩(wěn)相濾波算子。無論是Yan等[23]的PSPI算法，還是基于主能量傳播與偏振方向夾角估算的非穩(wěn)相濾波算法，都會(huì)損失偏振投影算子的精度，不足以應(yīng)對(duì)非均勻性和各向異性都很強(qiáng)的地質(zhì)情況。

3.1 偏振投影低秩近似算法

為了既能夠保證模式解耦的可靠性，又能夠大幅度降低計(jì)算復(fù)雜度，Cheng等[16]提出了精確表征偏振投影的空間—波數(shù)域積分算子及其低秩近似快速算法。對(duì)非均勻VTI介質(zhì)，偏振投影被表達(dá)成廣義傅里葉積分公式

(20)

式中am=(amx，amy，amz)T表示某種波模式對(duì)應(yīng)的偏振向量。假設(shè)空間或波數(shù)采樣點(diǎn)數(shù)為M，則這種混合域積分的計(jì)算復(fù)雜度為O(M2×lgM)，現(xiàn)有計(jì)算條件下很難承擔(dān)這樣的任務(wù)。根據(jù)式(8)，還可以構(gòu)建縱、橫波矢量分解的精確偏振投影算子，表示為廣義傅里葉積分形式(以x方向?yàn)槔?

(21)

根據(jù)式(20)和式(21)可知，式中的積分核同時(shí)依賴于空間域模型參數(shù)和傅里葉域波數(shù)成分。廣義傅里葉積分算子可以通過核函數(shù)近似來降低積分計(jì)算的成本。通常來說，地球介質(zhì)具有分層或分區(qū)變化的彈性特征，物理上也認(rèn)為該混合域核函數(shù)矩陣是低秩的。因此根據(jù)低秩近似理論[24]，該混合域核函數(shù)矩陣可以表示為三個(gè)矩陣的乘積

(22)

式中：W表示混合域核函數(shù)矩陣；B和C分別表示降維之后的波數(shù)和空間維度的混合域矩陣；Dmn為L(zhǎng)×N矩陣，其中L和N分別代表該算子矩陣在低秩近似意義下的列秩與行秩。如果給定矩陣分解的誤差門檻值為1.0×10-6，常見各向異性介質(zhì)模型偏振投影算子矩陣的秩均較小，通常為幾或十幾這樣的數(shù)量級(jí)。對(duì)于圖1a中兩層VTI模型，L和N均等于2。偏振投影模式解耦算子矩陣之所以存在上述低秩近似形式，主要是因?yàn)槭艿刭|(zhì)作用的控制，大多數(shù)地層都呈現(xiàn)一定的層狀或分塊特征，從而使該算子被那些最能代表模型非均勻性和各向異性特征的空間位置和傳播方向所主導(dǎo)[16]。

依據(jù)上述低秩分解(式(22))，式(20)中的積分式可退化成

(23)

式中j={x，y，z}。它們主要的計(jì)算量體現(xiàn)在N個(gè)矢量波場(chǎng)的正、反傅里葉變換，以及空間波場(chǎng)的N和L次加權(quán)求和。因?yàn)長(zhǎng)和N遠(yuǎn)遠(yuǎn)小于空間或波數(shù)采樣點(diǎn)數(shù)M，所以計(jì)算復(fù)雜度降低至O(L×M×lgM)或O(N×M×lgM)。這樣就以較小的計(jì)算成本合理逼近了式(20)中的廣義傅里葉積分運(yùn)算[16]。更重要的是，偏振投影算子矩陣的低秩分解完全是根據(jù)介質(zhì)模型本身特點(diǎn)自動(dòng)進(jìn)行的，不存在Yan等[23]的PSPI算法中設(shè)置參考模型類似的人為因素，這種模型自適應(yīng)特點(diǎn)對(duì)實(shí)際應(yīng)用非常重要。

3.2 模型分區(qū)策略

Wang[25]等認(rèn)為，對(duì)于一些大規(guī)模應(yīng)用問題，可把模型空間劃分成若干區(qū)域，分別在每個(gè)區(qū)域構(gòu)造低秩近似偏振投影算子，進(jìn)一步提高算法效率和并行擴(kuò)展?jié)摿ΑＨ鐖D3a所示，首先需要將模型進(jìn)行分區(qū)(黃色部分)，區(qū)域之間保留一定的重疊(綠色部分)，模型外圍有PML吸收邊界(藍(lán)色部分)，然后分別對(duì)每塊區(qū)域求取偏振矢量并構(gòu)建低秩近似投影算子，對(duì)彈性波場(chǎng)也進(jìn)行同樣地分區(qū)并施加窗函數(shù)。如圖3b所示，重疊區(qū)域的邊界衰減函數(shù)采用的是正弦和余弦函數(shù)，且滿足二者權(quán)重系數(shù)求和為1； 最后，對(duì)不同區(qū)域的波場(chǎng)分別施加低秩近似偏振投影，并將分離后的波場(chǎng)按區(qū)域拼接。

圖3 模型分區(qū)示意圖(a)及重疊區(qū)域邊界衰減系數(shù)設(shè)置(b)

盡管模型分區(qū)一定程度上可減小低秩近似解耦算法中FFT的總體成本，且一定程度降低了每個(gè)分區(qū)內(nèi)偏振投影算子矩陣的秩，有助于節(jié)省單個(gè)時(shí)間片波場(chǎng)的模式解耦計(jì)算成本。但是，模型和波場(chǎng)分區(qū)需要每個(gè)分區(qū)周圍的鑲邊處理，既帶來額外計(jì)算開銷，又在邊界附近引入一定的誤差，這些問題對(duì)彈性波逆時(shí)偏移或全波形反演等實(shí)際應(yīng)用中大量延拓時(shí)間步上的模式解耦來說是不可忽略的。

3.3 GPU加速策略

GPU并行加速在大規(guī)模波場(chǎng)數(shù)值模擬、逆時(shí)偏移成像領(lǐng)域已經(jīng)得到廣泛應(yīng)用，常采用的有OpenCL和CUDA兩種并行編程框架，現(xiàn)已開始使用OpenACC[29]。

從前文分析可知，正、反傅里葉變換在低秩近似模式解耦算法中占據(jù)大部分計(jì)算量，因此提高傅里葉變換計(jì)算效率非常關(guān)鍵。利用基于CUDA架構(gòu)的快速傅里葉變換并行算法庫cuFFT很好地解決了這個(gè)問題。其次，彈性波波模式解耦算法本身具有很好的數(shù)據(jù)獨(dú)立性，因此可利用GPU進(jìn)行細(xì)粒度的并行加速。

對(duì)于模型空間太大、單卡內(nèi)存不足的情況，需要通過多卡并行方式來解決。多GPU并行算法的基礎(chǔ)是沿著慢維方向?qū)δＰ涂臻g進(jìn)行分區(qū)，而核心問題是GPU間的通信及其效率。單節(jié)點(diǎn)內(nèi)GPU間的實(shí)時(shí)通信可以通過P2P技術(shù)實(shí)現(xiàn)，本文以單節(jié)點(diǎn)雙卡并行作為示例。如圖4所示，整個(gè)模型空間分成兩部分，分別對(duì)應(yīng)在GPU1和GPU2上計(jì)算。每個(gè)計(jì)算分區(qū)都包含計(jì)算區(qū)域、數(shù)據(jù)發(fā)送區(qū)域(S1、S2)和數(shù)據(jù)接收區(qū)域(R1、R2)，其中S1(S2)、R1(R2)分別對(duì)應(yīng)為GPU1(GPU2)上的數(shù)據(jù)交換(發(fā)送、接收)內(nèi)存區(qū)域。數(shù)據(jù)傳輸過程如圖中箭頭所示。P2P技術(shù)避免以CPU作為數(shù)據(jù)傳輸中轉(zhuǎn)站，減小了數(shù)據(jù)傳輸延遲。

圖4 雙卡并行模型分區(qū)示意圖

4 數(shù)值試驗(yàn)

本文以SEG Hess VTI介質(zhì)二維模型為例，測(cè)試本文發(fā)展的模式解耦高效算法。圖5為模型右半部分，尺寸為11.25km×9.38km(對(duì)應(yīng)的空間網(wǎng)格數(shù)為1800×1500)，震源位于(4.75km，3.81km)處。圖6為1400ms時(shí)的質(zhì)點(diǎn)速度水平和垂直分量的波前快照。圖7分別展示了本文依據(jù)4個(gè)主能量方向傳播與偏振矢量夾角構(gòu)建的偽散度/偽旋度算子、偏振投影低秩近似解耦算法以及基于模型分區(qū)(以分成四個(gè)區(qū)域?yàn)槔?的偏振投影低秩近似解耦算法得到的qP和qSV波場(chǎng)，可見幾種算法都取得了比較好的模式解耦效果。

圖5 SEG Hess VTI模型

圖6 SEG Hess VTI模型彈性波質(zhì)點(diǎn)速度場(chǎng)1400ms時(shí)刻的水平(a)和垂直(b)分量

圖7 SEG Hess VTI模型改進(jìn)的偏振投影近似算法和低秩近似算法波場(chǎng)分離結(jié)果

如果采用偽散度/偽旋度算子進(jìn)行模式解耦，需要在每個(gè)時(shí)間步計(jì)算坡印廷矢量，進(jìn)而構(gòu)建每個(gè)時(shí)間步的空間非穩(wěn)相濾波器。相反，對(duì)于偏振投影低秩近似解耦算法，僅需要提前一次性構(gòu)建出解耦濾波器，然后在各個(gè)時(shí)間步直接應(yīng)用。在單核CPU(Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz)平臺(tái)上，以2000個(gè)時(shí)間步彈性qP/qSV波模式解耦為例，采用依據(jù)4個(gè)主能量方向傳播與偏振矢量夾角構(gòu)建的偽散度/偽旋度算子，平均每個(gè)時(shí)間片波場(chǎng)模式解耦耗用時(shí)間為6.8s，而偏振投影低秩近似算法的平均耗時(shí)為5.8s，兩種算法在計(jì)算效率方面比較接近。

對(duì)于SEG Hess VTI模型低秩近似模式解耦，基于模型分區(qū)的低秩近似解耦GPU算法相比于常規(guī)的低秩近似解耦算法，單GPU卡(Nvidia Tesla K80)上改進(jìn)后的低秩近似算法實(shí)現(xiàn)的加速比達(dá)到了4以上，同時(shí)從分離結(jié)果上可以看出前者精度也基本沒有損失(適當(dāng)選取區(qū)域劃分的個(gè)數(shù))。對(duì)于彈性波逆時(shí)偏移和全波形反演等處理而言，需要每個(gè)時(shí)間的波場(chǎng)都施加模式解耦，基于模型分區(qū)和GPU加速的低秩近似解耦算法在保證精度的前提下明顯縮短了處理周期。此外，對(duì)于圖1的兩層VTI模型而言，基于模型分區(qū)和GPU加速的低秩近似解耦算法反而比常規(guī)的耗時(shí)多，原因在于CPU-GPU異構(gòu)平臺(tái)涉及它們之間的數(shù)據(jù)通信開銷，只有當(dāng)模型規(guī)模達(dá)到一定程度時(shí)，基于模型分區(qū)和GPU加速算法優(yōu)勢(shì)才能凸顯。

對(duì)比圖7a、圖7b與圖7g～圖7j可知，偏振投影低秩近似矢量解耦算法可以使得分離后的縱橫波場(chǎng)的相位和振幅保真，同時(shí)結(jié)合模型分區(qū)和GPU加速策略，使算法的精度和計(jì)算效率達(dá)到最佳。

5 結(jié)論

針對(duì)VTI介質(zhì)qP與qS波模式解耦高昂的計(jì)算成本，本文一方面基于現(xiàn)有的由主能量傳播與偏振方向偏差角構(gòu)造的偽散度/偽旋度算法，通過波場(chǎng)方向分解估計(jì)的多個(gè)主能量方向逼近偏振投影算子，實(shí)現(xiàn)了更精確的空間域非穩(wěn)相濾波算法，較好地解決了原算法在qP與qS波同時(shí)經(jīng)過相同空間區(qū)域時(shí)的模式泄漏問題。由于該方法利用坡印廷矢量估計(jì)的傳播與偏振矢量夾角逼近qP波的偏振方向(無法分別估計(jì)qSV和qSH波的偏振方向)，它構(gòu)建的模式解耦濾波只能實(shí)現(xiàn)qP與qS波分離，無法區(qū)分VTI介質(zhì)中的qSV與qSH波。另一方面，基于現(xiàn)有的偏振投影低秩近似算法，通過采用模型分區(qū)和GPU加速策略進(jìn)一步提高了算法效率。將模型空間分區(qū)和線程塊分配有機(jī)結(jié)合，利用CUDA架構(gòu)下的cuFFT快速傅里葉變換算法以及P2P數(shù)據(jù)通信技術(shù)，可以縮短這種高精度解耦算法的計(jì)算周期。對(duì)于二維VTI介質(zhì)彈性波模式解耦問題，本文兩種高效算法性價(jià)比相當(dāng)。第二種算法因?yàn)槟軌驅(qū)S波進(jìn)一步解耦，因此三維擴(kuò)展性更強(qiáng)。今后將結(jié)合三維情況下的彈性波逆時(shí)偏移和全波形反演實(shí)際應(yīng)用發(fā)揮該算法的優(yōu)勢(shì)。