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

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

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

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

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

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

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

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

1 方法理論

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

1.1 成像條件

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

(1)

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

(2)

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

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

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

共炮逆時偏移的真振幅三維角度域共成像點道集[11,39,40]可以按下式計算:

(3)

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

反射角和方位角計算如下:

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

(6)

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

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

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

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

1.3 波場分解

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

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

(7)

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

(8)

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

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

(9)

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

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

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

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

(10)

上行波可以表示為:

(11)

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

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

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

(12)

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

(13)

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

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

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

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

給定二維數(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)

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

(15)

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

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

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

(16)

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

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

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

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

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

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

1.3.4 數(shù)值算例

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

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

圖2 t=2s時刻的二維波場快照

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

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

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

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

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

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

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

1.3.5 波場分解拓展

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

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

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

2 應(yīng)用實例

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

2.1 二維Sigsbee2A例子

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

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

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

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

2.2 三維SEAM TTI例子

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

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

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

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

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

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

3 結(jié)論

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

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

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

(A3)

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

(A4)

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

(A5)

式中:ξ表示頻率或波數(shù),F(ξ)表示信號的傅里葉變換,Fh(ξ) 表示信號的希爾伯特變換。根據(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)式,可以觀察到下行波和上行波都可表示為原始波場及其在時間和深度方向進(jìn)行希爾伯特變換后波場的組合。因此,在時間-空間域中,下行波場和上行波場可以分別表示為公式(10)和公式(11)。

附錄B波場分解拓展

如公式(10)和公式(11)所示,上、下行波場都可以通過希爾伯特變換得到。在此附錄中,我們將此理論加以延伸。首先推導(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)