王美霞,張曉慧,張 釙,唐 冰,徐 昇

(1.Statoil Gulf Services,Houston 77042;2.前斯塔托爾北京技術(shù)服務(wù)有限公司,北京100022)

復(fù)雜構(gòu)造區(qū)的地震勘探,如大角度地層或鹽丘等需要高質(zhì)量的速度模型和高精度的成像方法。逆時(shí)偏移[1-4]利用雙程波傳播方程,能夠準(zhǔn)確模擬反射波、折射波、多次波等各種地震波,是目前最為精確的成像技術(shù)之一。

逆時(shí)偏移算法在20世紀(jì)80年代就已經(jīng)提出,但由于其巨大的計(jì)算量和超量的計(jì)算資源需求,直到21世紀(jì)初才有實(shí)際應(yīng)用。逆時(shí)偏移應(yīng)用面臨著理論和計(jì)算兩方面的挑戰(zhàn)。理論方面主要有兩點(diǎn):①逆時(shí)偏移低頻噪聲問題[4]。逆時(shí)偏移成像是通過對(duì)震源波場(chǎng)和檢測(cè)器波場(chǎng)應(yīng)用傳統(tǒng)的零延遲互相關(guān)成像條件[5]得到的,如果存在強(qiáng)反射,這個(gè)過程不僅在實(shí)際反射界面處產(chǎn)生圖像,也會(huì)沿著波路徑產(chǎn)生嚴(yán)重的低頻噪聲。這些噪聲可能掩蓋真實(shí)的構(gòu)造,給解釋工作造成困難。將這個(gè)傳統(tǒng)的成像條件和真振幅偏移方法[6-11]做比較,我們發(fā)現(xiàn)主要區(qū)別在于依賴于震源和檢波器波場(chǎng)夾角的權(quán)重函數(shù)。實(shí)際上,真振幅偏移中關(guān)于反射角的函數(shù)可以壓制很多反向散射噪聲,關(guān)鍵在于如何計(jì)算逆時(shí)偏移中的反射角。②快速生成偏移共成像點(diǎn)道集問題。共成像點(diǎn)道集是質(zhì)量控制、速度建模、屬性分析等的重要輸入數(shù)據(jù),被廣泛應(yīng)用于復(fù)雜地區(qū)數(shù)據(jù)處理。常用的共成像點(diǎn)道集是共偏移距道集,也就是使用地面偏移距作為索引的道集。然而,在復(fù)雜情況下,共偏移距道集具有很多由于多路徑傳播而引入的偏移假像,而且很難直接由逆時(shí)偏移得到高質(zhì)量的共偏移距道集,需要進(jìn)行特殊處理[12]。相比之下,由逆時(shí)偏移產(chǎn)生角道集[13-15]就比較自然。此時(shí)角道集由地下反射角作為索引,偏移噪聲也相對(duì)較弱。如果使用真振幅逆時(shí)偏移,角道集還可以為AVA分析提供輸入數(shù)據(jù)。角道集可以使用拓展成像條件的方法得到[16-19],也可以直接由波場(chǎng)產(chǎn)生[11-20]。這些方法在二維情況下計(jì)算效率都比較高,但是在三維情況下計(jì)算量都很大,所以在實(shí)現(xiàn)逆時(shí)偏移算法時(shí)首先要考慮算法效率以及算法的并行計(jì)算和優(yōu)化等。計(jì)算效率方面,引入Poynting矢量方法提高了生成角道集的計(jì)算效率,但是由于檢波器波場(chǎng)比較復(fù)雜,一般Poynting矢量?jī)H被用于震源波場(chǎng),使用逆時(shí)偏移成像或反射界面傾角而不是檢波器端波傳播方向信息[21-22]。在三維道集逆時(shí)偏移應(yīng)用中,面對(duì)存儲(chǔ)量大、計(jì)算量高等挑戰(zhàn),Poynting矢量方法更具吸引力。YOON等[23]利用Poynting矢量方法估計(jì)波傳播方向,該方法的缺點(diǎn)是在一個(gè)時(shí)空點(diǎn)處,Poynting矢量只能給出一個(gè)方向,因此,當(dāng)波場(chǎng)復(fù)雜時(shí),不同傳播方向的波混疊在一起,Poynting矢量就不能準(zhǔn)確計(jì)算方向。若用于角道集計(jì)算,就會(huì)使得有些能量被成像在錯(cuò)誤的角度上,進(jìn)而降低角道集質(zhì)量或?qū)е陆堑兰嬖诓糠帜芰咳笔?這嚴(yán)重影響Poynting矢量方法在檢波器端波場(chǎng)的應(yīng)用。

本文進(jìn)一步分析了Poynting矢量方法。我們認(rèn)為只需要在使用Poynting矢量計(jì)算波傳播方向之前對(duì)波場(chǎng)進(jìn)行分解,便可以將其同時(shí)用于震源和檢波器波場(chǎng)。波場(chǎng)分解越精細(xì),Poynting矢量就越精確,但計(jì)算量也會(huì)隨之增加。對(duì)于一般的地面觀測(cè)系統(tǒng),初至波主要是震源的下行波,反射波主要是檢波器的下行波。如果模型中存在強(qiáng)速度差界面,波場(chǎng)中初至波和反射波就會(huì)疊加,因此進(jìn)行簡(jiǎn)單的上、下行波場(chǎng)分解[2-4]就可以顯著提高Poynting矢量的精度。

傳統(tǒng)的波場(chǎng)分解算法[2-4]一般是通過傅里葉變換[25]在頻率-波數(shù)域進(jìn)行的。分解計(jì)算效率取決于快速傅里葉變換的實(shí)現(xiàn)效率。然而,傅里葉變換是全局運(yùn)算,計(jì)算耗時(shí)。另一種方法是使用希爾伯特變換(HT)[24,26-28]進(jìn)行波場(chǎng)分解。下行波和上行波可以通過原始波場(chǎng)和在時(shí)間-深度方向進(jìn)行希爾伯特變換后的波場(chǎng)組合得到[2-7]。這和利用傅里葉變換得到的結(jié)果一致。優(yōu)點(diǎn)是希爾伯特變換的實(shí)現(xiàn)比較靈活;可以通過快速傅里葉變換實(shí)現(xiàn);也可以在時(shí)間-空間域使用有限長度卷積算子實(shí)現(xiàn),以便進(jìn)一步進(jìn)行算法優(yōu)化以提高計(jì)算效率。本文使用SSE和多線程優(yōu)化,提高利用希爾伯特變換進(jìn)行波場(chǎng)分解的計(jì)算效率。

利用希爾伯特變換進(jìn)行波場(chǎng)分解時(shí),需要考慮時(shí)間和空間方向的希爾伯特變換。有兩種方法可以實(shí)現(xiàn)波場(chǎng)在時(shí)間方向的希爾伯特變換。一種是利用希爾伯特變換后的震源函數(shù)進(jìn)行波傳播得到一個(gè)新的波場(chǎng)[27],該方法會(huì)加倍波傳播計(jì)算量和存儲(chǔ)量,因此讀寫時(shí)間也會(huì)增加。另一種辦法是只使用原始波場(chǎng),在計(jì)算道集的過程中邊計(jì)算邊對(duì)原始波場(chǎng)進(jìn)行希爾伯特變換。若效率和算法不是瓶頸,后者需要的數(shù)據(jù)讀寫較少。因此,本文采用后者。

本文首先簡(jiǎn)單回顧真振幅逆時(shí)偏移成像和角道集理論以及Poynting矢量方法。其次,我們證明通過傅里葉變換或希爾伯特變換都能夠進(jìn)行波場(chǎng)分解,且理論上等價(jià),進(jìn)一步使用單指令多數(shù)據(jù)流(SIMD)優(yōu)化希爾伯特變換算法,并使用數(shù)值例子驗(yàn)證了優(yōu)化后算法的精度和效率。然后,將優(yōu)化的希爾伯特算法用于波場(chǎng)分解,進(jìn)一步使用多線程優(yōu)化并行算法。利用數(shù)值計(jì)算例子說明這種方法在得到和傳統(tǒng)傅里葉方法一致結(jié)果同時(shí),能夠顯著提高數(shù)值計(jì)算效率。最后,將這種加速的波場(chǎng)分解算法應(yīng)用于逆時(shí)偏移產(chǎn)生角道集,利用二維Sigsbee2A和三維SEAM TTI計(jì)算實(shí)例說明本文方法的有效性和穩(wěn)定性。

1 方法理論

利用雙程波方程進(jìn)行波傳播的數(shù)值計(jì)算[29-35],可以得到震源波場(chǎng)us(x,t)和檢波器波場(chǎng)ur(x,t)。這里t表示時(shí)間,向量x表示空間位置。在二維介質(zhì)中,x=(z,x)。在三維介質(zhì)中,x=(z,x,y)。下面我們首先回顧成像條件和角道集理論,之后引出一種快速的波場(chǎng)分解算法。

1.1 成像條件

偏移成像的傳統(tǒng)互相關(guān)成像條件[5]可以表示為:

(1)

式中:xs和ys表示震源位置;ω為角頻率。此成像條件在反射點(diǎn)處產(chǎn)生像,但同時(shí)也在滿足成像條件的路徑上產(chǎn)生假像。假像和真像的區(qū)別在于在假像處震源波場(chǎng)和檢波器波場(chǎng)的傳播方向相異。YOON等[23]提出只有當(dāng)震源和檢波器波傳播方向所成角度在一定范圍內(nèi)時(shí)才使用互相關(guān)條件成像,且波傳播方向由Poynting矢量計(jì)算得到,該方法的一個(gè)缺點(diǎn)是很難選取用于成像的最大角度。COSTA等[36]提出了依賴于角度的光滑函數(shù)作為偏移權(quán)重,在一定程度上克服了最大成像角度選取的困難。與真振幅偏移[6-7,10-11]相比,其和傳統(tǒng)成像條件的區(qū)別在于依賴于震源和檢波器波場(chǎng)夾角的權(quán)重函數(shù)。因此,一個(gè)更加適當(dāng)?shù)某上駰l件可以表達(dá)為[6,11,36-38]:

(2)

式中:θ表示反射角。由于偏移噪聲是當(dāng)震源和檢波

器波傳播方向成較大角度時(shí)產(chǎn)生的,所以權(quán)重函數(shù)cos2θ可以有效壓制偏移中的低頻噪聲。

1.2 使用Poynting矢量產(chǎn)生真振幅逆時(shí)偏移角道集

共炮逆時(shí)偏移的真振幅三維角度域共成像點(diǎn)道集[11,39,40]可以按下式計(jì)算:

(3)

式中:R(x,θ0,φ0)表示在空間位置x處方位角φ0和反射角θ0方向上的反射系數(shù);δ(θ-θ0)和δ(φ-φ0)指狄拉克函數(shù),它們將積分限制在要求的角度θ0和φ0,v(x)為速度。

反射角和方位角計(jì)算如下:

其中,p=1/2(ps+pr),x,z為坐標(biāo)軸的單位向量;ps和pr分別表示震源和檢波器的波傳播方向。震源波傳播方向可以按下式利用Poynting矢量進(jìn)行計(jì)算:

(6)

類似地,檢波器波傳播方向pr也可以使用Poynting矢量計(jì)算得到。

Poynting矢量計(jì)算是在時(shí)間-空間域進(jìn)行的,當(dāng)波場(chǎng)比較復(fù)雜時(shí),Poynting矢量并不精確。圖1a展示了從Sigsbee2A模型得到的一個(gè)波場(chǎng)快照,可以看到波場(chǎng)很復(fù)雜,幾種不同類型的波交織在一起,這使得Poynting矢量很難給出一個(gè)正確的方向。然而,如果我們僅關(guān)注下行波(圖1b),波場(chǎng)就變得簡(jiǎn)單很多,Poynting矢量可以較準(zhǔn)確地計(jì)算波傳播方向。

圖1 由Sigsbee2A得到的波場(chǎng)快照(a)及其下行波場(chǎng)分量(b)

因此,波場(chǎng)分解,特別是上、下行波場(chǎng)分解(對(duì)于地表觀測(cè)系統(tǒng)),能顯著提高Poynting矢量計(jì)算方向的精度,這也可以進(jìn)一步提高角道集的質(zhì)量。

1.3 波場(chǎng)分解

1.3.1 利用傅里葉變換進(jìn)行波場(chǎng)分解

傳統(tǒng)方法利用傅里葉變換在頻率-波數(shù)域進(jìn)行波場(chǎng)分解。約定時(shí)間-深度方向的二維傅里葉變換如下:

(7)

于是下行波在頻率-波數(shù)域可表示為:

(8)

上行波在頻率-波數(shù)域可表示為:

UU(kz,ω)=U(kz,ω)-UD(kz,ω)

(9)

在公式(8)和公式(9)中,對(duì)于滿足ω·kz=0的U(kz,ω),我們使用了系數(shù)1/2以便使得UD(kz,ω)+UU(kz,ω)=U(kz,ω)。這是合理的,因?yàn)闈M足ω·kz=0的波既非上行波也非下行波。

1.3.2 利用希爾伯特變換進(jìn)行波場(chǎng)分解

SHEN等[27]指出上、下行波可以使用希爾伯特變換進(jìn)行波場(chǎng)分離。他們使用一對(duì)波場(chǎng)實(shí)現(xiàn)分解:震源函數(shù)傳播得到的原始波場(chǎng)和進(jìn)行希爾伯特變換后的震源函數(shù)傳播得到的波場(chǎng)。這種方法很有效,但計(jì)算效率不高,它需要兩倍的波傳播計(jì)算量。本文提出另一種實(shí)現(xiàn)方法,僅使用原始波場(chǎng),在計(jì)算道集過程中進(jìn)行希爾伯特變換。重點(diǎn)以上、下行波分解來描述算法,并將該方法拓展到其它方向的波場(chǎng)分解。

通過傅里葉變換和希爾伯特變換之間的關(guān)系,最終下行波可以表示為:

(10)

上行波可以表示為:

(11)

式中:Ht[u(x,t)]表示波場(chǎng)u(x,t)在時(shí)間方向的希爾伯特變換;Hz[u(x,t)]表示深度方向的希爾伯特變換。公式(10)和公式(11)與SHEN等[27]提出的公式一致,具體推導(dǎo)過程見附錄A。

公式(10)和公式(11)提供了另一種波場(chǎng)分解的辦法。希爾伯特變換可以簡(jiǎn)單通過傅里葉變換實(shí)現(xiàn),但是這樣計(jì)算效率與傳統(tǒng)頻率-波數(shù)域的波場(chǎng)分解一樣。為了提高計(jì)算效率,我們提出直接在時(shí)間-空間域通過卷積進(jìn)行希爾伯特變換,并進(jìn)一步使用SSE優(yōu)化該算法。

對(duì)于一維輸入數(shù)組aj,j=0,1,…,n-1,使用(2m+1)長度卷積算子進(jìn)行的離散希爾伯特變換如下:

(12)

我們注意到進(jìn)行希爾伯特變換的卷積算子是反對(duì)稱的,即:hk=-h2m-k,k=0,1,…,m。而且大約一半的系數(shù)hk都是0?；谶@兩條性質(zhì),公式(12) 可以進(jìn)一步改寫為如下形式以減少計(jì)算量:

(13)

式中:k+=2表示k的步長為2。

以下的實(shí)現(xiàn)和數(shù)值計(jì)算實(shí)例都是基于公式(13)。

1.3.3 希爾伯特變換波場(chǎng)分解算法的優(yōu)化

若將波場(chǎng)分解用于逆時(shí)偏移中進(jìn)行上、下行波分解,我們需要時(shí)間和深度方向的希爾伯特變換(如公式(10)和公式(11))。約定波場(chǎng)在內(nèi)存中的存儲(chǔ)規(guī)律為最快維為深度,最慢維為時(shí)間。因此,提高波場(chǎng)分解算法速度的關(guān)鍵是分別優(yōu)化快維和慢維方向的希爾伯特變換。下面先說明如何優(yōu)化二維數(shù)組的快維和慢維的希爾伯特變換,然后,將其用于波場(chǎng)分解。

給定二維數(shù)組ai,j,i=0,1,…,n2-1,且j=0,1,…,n1-1(下角標(biāo)i表示慢維索引,j表示快維索引),記數(shù)組{ai,j}的希爾伯特變換為數(shù)組{bi,j}。傳統(tǒng)方法在快維方向的希爾伯特變換可直接表示為:

(14)

若要計(jì)算慢維方向的希爾伯特變換,傳統(tǒng)方法是先通過轉(zhuǎn)置將慢維變?yōu)榭炀S,再由公式(14)計(jì)算,最后再做一次數(shù)組轉(zhuǎn)置。對(duì)數(shù)據(jù)做兩次轉(zhuǎn)置增加了計(jì)算量。我們提出一種不用轉(zhuǎn)置的算法,即按照公式(15)直接對(duì)慢維進(jìn)行數(shù)值運(yùn)算:

(15)

由于我們使用相對(duì)較短的卷積算子,k相對(duì)于慢維長度n2較小,就可以使用SIMD沿快維j方向加速算法,即將每一列ai-m+k,*-ai+m-k,*看作一個(gè)元素,同一個(gè)系數(shù)h2m-k與該元素相乘,這樣就可以使用SSE同時(shí)處理多個(gè)數(shù)據(jù)。我們將這種優(yōu)化的算法記為慢維方向的向量化希爾伯特變換(vectorized Hilbert transform,VHT)。

通過使用SSE優(yōu)化,慢維方向的希爾伯特變換不再需要對(duì)數(shù)組進(jìn)行轉(zhuǎn)置。理論上,數(shù)值計(jì)算效率可以提高到4倍(SSE)或8倍(AVX,Advanced Vetor Extensions)。但是考慮到從緩存中獲取數(shù)據(jù)也需要時(shí)間,實(shí)際提高水平可能略低于理論值。

對(duì)于快維方向的希爾伯特變換,我們發(fā)現(xiàn)如下使用SSE優(yōu)化的算法比直接利用公式(14)計(jì)算更加高效:

(16)

在以上過程中,二維轉(zhuǎn)置可以使用優(yōu)化的SSE轉(zhuǎn)置函數(shù)(_MM_TRANSPOSE4_PS)實(shí)現(xiàn)。盡管使用了轉(zhuǎn)置,SSE向量化仍然可以大幅提高計(jì)算效率。我們稱此為快維方向的VHT。

由于重點(diǎn)是上、下行波分解,所以進(jìn)行波場(chǎng)分解算法的核心是二維波場(chǎng)的分解。對(duì)于三維情況,y方向是最慢維,只需對(duì)每一個(gè)y切面波場(chǎng)進(jìn)行二維波場(chǎng)分解。在使用SSE向量化的同時(shí),還可以使用OpenMP進(jìn)行多線程并行計(jì)算以提高效率。對(duì)于地震成像應(yīng)用,目前的緩存一般足夠容納二維子波場(chǎng)。因此,基于公式(10)和公式(11),我們提出以下流程進(jìn)行高效波場(chǎng)分解:

1) 對(duì)每一個(gè)ix,將二維子波場(chǎng)u(z,ix,t)從內(nèi)存加載到緩存中,并使用慢維方向的向量化希爾伯特變換計(jì)算Ht[u(z,ix,t)];

2) 對(duì)每一個(gè)it,將二維子波場(chǎng)Ht[u(z,x,it)]從內(nèi)存加載到緩存中,并使用快維方向的VHT計(jì)算HzHt[u(z,x,it)];

3) 使用公式(10)(或公式(11))得到下行(或上行)波場(chǎng)。

對(duì)于三維波場(chǎng)的分解,我們只需要對(duì)每一個(gè)iy切片進(jìn)行以上二維波場(chǎng)分解。

1.3.4 數(shù)值算例

先以一個(gè)二維波場(chǎng)快照u(z,x)為例說明VHT的計(jì)算效率及其在波場(chǎng)分解中的應(yīng)用。輸入波場(chǎng)是由常速度v=1790m/s模型通過偽譜法[41]模擬產(chǎn)生,震源函數(shù)是主頻為14Hz的Ricker子波,我們對(duì)子波進(jìn)行了3Hz的高通濾波。波傳播網(wǎng)格為dx=dz=15m,nz=nx=701,且深度z為快維。震源在區(qū)域中心。我們使用一個(gè)相對(duì)較長的卷積算子,m=40,即81個(gè)點(diǎn),進(jìn)行希爾伯特變換以便得到較好的精度。

首先,我們驗(yàn)證與傳統(tǒng)方法相比,慢維方向的VHT能夠提高計(jì)算效率。圖2是t=2s時(shí)刻的波場(chǎng)快照。圖3a是使用向量化希爾伯特變換計(jì)算得到的二維波場(chǎng)在慢維x方向的希爾伯特變換波場(chǎng)。圖3b比較了慢維方向的VHT和傳統(tǒng)希爾伯特變換的計(jì)算效率。為了得到可靠的計(jì)算時(shí)間統(tǒng)計(jì),我們重復(fù)了5000次計(jì)算。從圖3b可以看到,對(duì)于慢維方向的希爾伯特變換,向量化算法比傳統(tǒng)算法效率高三倍多。為檢驗(yàn)VHT的精度,我們?nèi)D3a 中一條深度方向的記錄與傅里葉方法得到的結(jié)果進(jìn)行比較,見圖3c。由于傅里葉方法對(duì)于所有波數(shù)都是精確的,可以認(rèn)為其得到的結(jié)果為理想結(jié)果。圖3c中VHT結(jié)果與傅里葉方法的結(jié)果十分吻合,說明本文的VHT方法準(zhǔn)確有效。在進(jìn)行VHT時(shí),卷積算子越長,精度越高。在實(shí)際應(yīng)用中,如果不需要特別高的精度,可以使用一個(gè)較短的卷積算子,這樣可以進(jìn)一步提高計(jì)算效率。

圖2 t=2s時(shí)刻的二維波場(chǎng)快照

其次,我們?cè)衮?yàn)證快維方向的VHT。圖4a是圖3a 所示波場(chǎng)快照在深度方向的希爾伯特變換。圖4b 比較了VHT方法和傳統(tǒng)HT方法用于計(jì)算快維方向希爾伯特變換的計(jì)算效率。我們同樣進(jìn)行了5000次重復(fù)計(jì)算。可以看到向量化方法可以減少大概一半的計(jì)算時(shí)間。

圖3 使用VHT得到的波場(chǎng)快照(a)、VHT和傳統(tǒng)希爾伯特變換(HT)計(jì)算效率對(duì)比(b)以及深度z=5250m處VHT結(jié)果和傅里葉變換結(jié)果(FT)比較(c)

圖4 快維(z方向)的希爾伯特變換a 使用VHT得到的波場(chǎng)快照在z方向的希爾伯特變換; b z 方向的VHT和傳統(tǒng)希爾伯特變換(HT)的效率比較(計(jì)算重復(fù)了5000次)

從以上兩個(gè)實(shí)例(圖3b和圖4b),我們可以看到:傳統(tǒng)HT方法在計(jì)算慢維方向希爾伯特變換時(shí)比計(jì)算快維方向希爾伯特變換時(shí)需要時(shí)間更多;相反地,向量化方法在計(jì)算快維方向希爾伯特變換時(shí)比計(jì)算慢維方向希爾伯特變換時(shí)需要更多時(shí)間。這均是由于是否需要數(shù)據(jù)轉(zhuǎn)置所致。

最后,我們進(jìn)一步檢驗(yàn)VHT用于波場(chǎng)分解的綜合計(jì)算效率,并將其和傳統(tǒng)的利用傅里葉變換的波場(chǎng)分解進(jìn)行比較。為了合理地進(jìn)行比較,在進(jìn)行傅里葉變換時(shí)使用了FFTW(the Faster Fourier Transform in the West)[42]。波場(chǎng)時(shí)間方向采樣為:nt=660,dt=4ms。圖5a和圖5b分別是由本文算法和傅里葉方法得到的下行波場(chǎng)在t=2s時(shí)的波場(chǎng)快照。圖5c 比較了本文方法和傅里葉變換方法的計(jì)算效率。兩種方法都是用16核計(jì)算機(jī)進(jìn)行計(jì)算。為得到可靠的計(jì)算時(shí)間統(tǒng)計(jì),運(yùn)算重復(fù)了100次。可以看到本文提出的使用VHT的分解方法能夠很大地提高計(jì)算效率,比傳統(tǒng)使用傅里葉變換進(jìn)行分解的方法快了6倍多。為進(jìn)行更詳細(xì)的精度比較,圖6比較了從圖5a和圖5b中取出的x=4950m處一條深度方向的記錄,可以看到兩種方法得到的結(jié)果非常一致。

圖5 上、下行波分解示例(采用t=2s 時(shí)刻的下行波分量)a VHT結(jié)果; b FFTW計(jì)算結(jié)果; c 兩者效率比較(計(jì)算重復(fù)100次)

圖6 下行波場(chǎng)中x= 4950m 處深度方向記錄的比較

1.3.5 波場(chǎng)分解拓展

以上主要闡述了如何進(jìn)行上、下行波分解,我們也可以將該方法拓展到其他方向的波場(chǎng)分解。為了理論的完整性,附錄B列出了其他方向波分量的公式,包括左行波(ul)、右行波(ur)、左上行波(uul)、右上行波(uur)、左下行波(udl)、右下行波(udr)、左上后行波(uulb)、左上前行波(uulf)、右上后行波(uurb)、右上前行波(uurf)、左下后行波(udlb)、左下前行波(udlf)、右下后行波(udrb)、右下前行波(udrf)。

同理,用于上、下行波分解的優(yōu)化方法也可以用于其它方向波場(chǎng)的分解。這里給出一個(gè)左、右行波場(chǎng)分解的例子來說明其有效性。我們同樣使用圖2所使用的模型。圖7a和圖7b分別是t=2s時(shí)刻的左行波和右行波。圖7c比較了兩種方法的計(jì)算效率。可以看到,本文提出的VHT方法能夠顯著提高左、右行波分解的計(jì)算效率。值得注意的是本文方法用于左、右行波分解時(shí)不需要任何數(shù)據(jù)轉(zhuǎn)置,這一點(diǎn)也優(yōu)于傳統(tǒng)方法。

圖7 向量化希爾伯特變換得到的t=2s 時(shí)刻的左、右行波分解結(jié)果a 左行波分量; b 右行波分量; c 兩者效率比較(計(jì)算重復(fù)100次)

2 應(yīng)用實(shí)例

將波場(chǎng)分解算法用于逆時(shí)偏移成像來檢驗(yàn)其在逆時(shí)偏移和角道集生成中的作用。偏移像可以通過疊加角道集得到,這使得成像更加靈活,例如可以只疊加某一角度范圍內(nèi)的道集或使用合適的關(guān)于角度的權(quán)重函數(shù)。下面介紹了二維Sigsbee2A模型和三維SEAM TTI 模型例子。

2.1 二維Sigsbee2A例子

首先使用Sigsbee2A模型和數(shù)據(jù)。偏移中共使用500炮的合成地震數(shù)據(jù)。最大頻率為36Hz。為了準(zhǔn)確模擬波場(chǎng),我們使用偽譜法[41]進(jìn)行波傳播,并對(duì)時(shí)間頻散進(jìn)行了校正[34]。圖8a是不使用波場(chǎng)分解時(shí)得到的角道集,可以觀察到鹽丘下邊界之上道集并不連續(xù)也不清晰。圖8b為使用VHT進(jìn)行波場(chǎng)分解得到的結(jié)果,道集連續(xù)性變好,噪聲也得到壓制。這說明波場(chǎng)分解能夠提高Poynting矢量方向計(jì)算的精度。圖9是使用互相關(guān)成像條件(方程(1))和未進(jìn)行分解的波場(chǎng)得到的疊加圖像,從圖上可以看到嚴(yán)重的低頻噪聲,以致掩蓋了真實(shí)的反射界面。圖10a是使用包含角度權(quán)重的成像條件(方程(2))和未進(jìn)行分解的波場(chǎng)得到的疊加圖像,圖10b是使用包含角度權(quán)重的成像條件(方程(2))和下行波得到的疊加圖像。比較圖9和圖10a,可以看到依賴于角度的權(quán)重函數(shù)可以壓制很多噪聲,鹽丘以下的反射層也變得清楚。但是圖10a中還存在一些噪聲,使鹽丘下邊界之上的圖像比較模糊,如圖中紅色箭頭所示區(qū)域。若進(jìn)一步使用波場(chǎng)分解,疊加圖像就變得比較清晰,質(zhì)量得到提高,如圖10b 所示。這也說明波場(chǎng)分解之后使用包含角度權(quán)重的成像條件(方程(2))可以有效降低偏移噪聲。

圖8 Sigsbee2A模型的角道集(道集抽取位置6004.56～25816.56m，間隔1524m)a 不使用波場(chǎng)分解得到的角道集; b 使用VHT進(jìn)行波場(chǎng)分解得到的角道集

圖9 利用互相關(guān)成像條件得到的二維Sigsbee2A模型疊加圖像(沒有使用波場(chǎng)分解)

圖10 利用0～60°角道集疊加得到二維Sigsbee2A模型疊加圖像(成像條件為依賴于角度的權(quán)重函數(shù))a 未使用波場(chǎng)分解; b 采用VHT對(duì)下行波場(chǎng)進(jìn)行波場(chǎng)分解

2.2 三維SEAM TTI例子

實(shí)際上,我們更加重視三維應(yīng)用的計(jì)算效率。為此,我們使用三維SEAM (SEG Advanced Modeling) TTI 模型來說明本文算法在三維逆時(shí)偏移應(yīng)用中的計(jì)算效率。SEAM 模型區(qū)域?yàn)?0km×35km×15km。包含一個(gè)復(fù)雜的鹽丘[43],這在墨西哥灣地區(qū)非常典型。不僅模型大小符合實(shí)際情況,合成地震數(shù)據(jù)也比較符合真實(shí)情形。由于計(jì)算資源有限,我們選取了215炮數(shù)據(jù)進(jìn)行測(cè)試。圖11展示了炮的位置。在算法中,我們使用了XU 等[33]提出的擬P波方程,這樣不僅效率高而且沒有S波產(chǎn)生的噪聲。測(cè)試中,偏移最高頻率為18Hz。以下我們選取位置y=32250m

處觀測(cè)偏移結(jié)果。圖12a和圖12b分別是不使用和使用波場(chǎng)分解得到的結(jié)果?？梢钥吹?圖12a中,特別是紅色橢圓標(biāo)記區(qū)域,丟失了很多成像細(xì)節(jié);圖12b相對(duì)清晰很多,也有更多細(xì)節(jié)信息,特別是在鹽丘之上。這個(gè)比較說明即便對(duì)于非常復(fù)雜的三維模型,波場(chǎng)分解也能夠提高Poynting矢量的精度,從而提高角道集的質(zhì)量。圖13是不使用波場(chǎng)分解得到的疊加偏移圖像(使用成像條件(1)),可以看到嚴(yán)重的噪聲,真實(shí)的構(gòu)造被掩蓋。圖14是使用成像條件(2)得到的圖像,圖14a和圖14b分別使用波場(chǎng)分解前后得到的結(jié)果。比較圖14a 和圖13,可見關(guān)于角度的權(quán)重函數(shù)能夠壓制一部分噪聲,但是圖像質(zhì)量還不是很高。進(jìn)一步使用波場(chǎng)分解,可以得到更加清晰、細(xì)節(jié)更豐富的圖像(圖14b)。因此依賴于角度的權(quán)重函數(shù)可以壓制一些噪聲,波場(chǎng)分解可以進(jìn)一步提高成像質(zhì)量,去除偏移噪聲。

圖11 三維 SEAM TTI模型實(shí)例(選取215 炮的位置)

圖12 三維 SEAM TTI模型角道集(角道集位置11～25km，間隔為1km)a 未進(jìn)行波場(chǎng)分解得到的角道集; b 利用VHT波場(chǎng)分解后得到的角道集

圖13 三維 SEAM TTI模型疊加偏移圖像(采用互相關(guān)成像條件(方程(1)),未進(jìn)行波場(chǎng)分解)

圖14 利用0～60°的角道集疊加得到三維 SEAM TTI模型疊加偏移圖像(采用依賴于角度的權(quán)重函數(shù) (公式(2)成像)a 未使用波場(chǎng)分解; b 采用VHT進(jìn)行波場(chǎng)分解

3 結(jié)論

角度域共成像點(diǎn)道集對(duì)復(fù)雜構(gòu)造的地震成像十分重要,可以使用Poynting矢量高效計(jì)算波傳播方向,但是當(dāng)波場(chǎng)復(fù)雜時(shí),傳統(tǒng)的Poynting方法精度較低。為提高Poynting矢量的精度,可以在計(jì)算Poynting矢量之前將波場(chǎng)進(jìn)行分解,再進(jìn)行Poynting矢量計(jì)算,得到高精度的傳播方向,從而得到高質(zhì)量的角道集和偏移圖像。我們提出了VHT和多核并行的波場(chǎng)分解算法。經(jīng)過此優(yōu)化,波場(chǎng)分解算法加快大約5倍。二維Sigsbee2A和三維SEAM TTI的數(shù)值計(jì)算實(shí)例均說明本文的波場(chǎng)分解算法能夠有效去除逆時(shí)偏移圖像和角道集中的噪聲,提高精度和計(jì)算效率。該方法并不局限于上、下行波場(chǎng)分解,可以很容易拓展到其它方向的波場(chǎng)分解。本文算法也可以進(jìn)一步應(yīng)用于最小二乘偏移和全波形反演以消除成像噪聲和提高成像質(zhì)量。

致謝:感謝Hongbo Zhou 和Jun Mu 的幫助和討論。感謝MIT的FFTW,SMAART JV的Sigsbee2A模型和數(shù)據(jù),以及SEG的SEAM模型和數(shù)據(jù)。

附錄A 公式(10)和公式(11)的推導(dǎo)

我們首先定義如下兩個(gè)輔助函數(shù):

通過使用輔助函數(shù)(A1)式和 (A2)式,可以將頻率-波數(shù)域下行波的定義公式(8)改寫成如下形式:

(A3)

同理,使用輔助函數(shù)(A1)式和 (A2)式,可以將頻率-波數(shù)域上行波的定義公式(9)改寫成如下形式:

(A4)

對(duì)于任意信號(hào),其希爾伯特變換和傅里葉變換之間存在如下關(guān)系:

(A5)

式中:ξ表示頻率或波數(shù),F(ξ)表示信號(hào)的傅里葉變換,Fh(ξ) 表示信號(hào)的希爾伯特變換。根據(jù)我們關(guān)于希爾伯特變換和傅里葉變換的定義,σ(ξ) 是如下函數(shù):

(A6)

根據(jù)函數(shù)定義(A1)式、(A2)式和 (A6)式,我們可得到以下g1(ξ),g2(ξ)關(guān)于σ(ξ)的表示:

將(A7)式和(A8)式代入(A3)式可以得到下行波的表達(dá)式:

(A9)

同理,將(A7)式和 (A8)式代入(A4)式可以得到上行波的表達(dá)式:

(A10)

(A9)式和 (A10)式是下行波和上行波在頻率-波數(shù)域的表達(dá)式。再使用傅里葉變換和希爾伯特變換的關(guān)系(A5)式,可以觀察到下行波和上行波都可表示為原始波場(chǎng)及其在時(shí)間和深度方向進(jìn)行希爾伯特變換后波場(chǎng)的組合。因此,在時(shí)間-空間域中,下行波場(chǎng)和上行波場(chǎng)可以分別表示為公式(10)和公式(11)。

附錄B波場(chǎng)分解拓展

如公式(10)和公式(11)所示,上、下行波場(chǎng)都可以通過希爾伯特變換得到。在此附錄中,我們將此理論加以延伸。首先推導(dǎo)左行波、右行波和左上行波的表達(dá)式。其它方向波分量可以同理得到,我們僅給出最終表達(dá)式。

在頻率-波數(shù)域中,右行波可以定義為:

(B1)

左行波可以定義為:

(B2)

類似附錄A中上、下行波的推導(dǎo),可以得到以下由希爾伯特變換表示的左、右行波表達(dá)式:

左上行波可以通過進(jìn)一步將上行波uu(x,t)進(jìn)行左、右行波分解得到:

(B5)

將上行波表達(dá)式(公式(11))代入此方程可以得到:

(B6)

同樣,右上行波可以表達(dá)為:

(B7)

左下行波可以表達(dá)為:

(B8)

右下行波可以表達(dá)為:

(B9)

左上后行波可以表達(dá)為:

(B10)

左上前行波可以表達(dá)為:

(B11)

右上后行波可以表達(dá)為:

(B12)

右上前行波可以表達(dá)為:

(B13)

左下后行波可以表達(dá)為:

(B14)

左下前行波可以表達(dá)為:

(B15)

右下后行波可以表達(dá)為:

(B16)

右下前行波可以表達(dá)為:

(B17)