李 博,劉志成,李小愛,許 璐

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.中國石油化工股份有限公司西北油田分公司,新疆烏魯木齊830011)

基于雙程波動方程的逆時偏移方法在實際工業(yè)生產(chǎn)中能夠應(yīng)對各種復雜的地震波成像問題[1-2],因此一直都是石油勘探行業(yè)的研究熱點。該方法采用有限差分方法進行震源波場正傳與檢波點波場反傳,并將兩種波場沿時間方向做零延遲互相關(guān)成像[3-4]。利用雙程波波場模擬算法能夠正確描述反射波、折射波、回轉(zhuǎn)波、多次波等復雜波場現(xiàn)象,處理復雜構(gòu)造地區(qū)的地震成像問題。然而,由于波場傳播算子沒有明確區(qū)分波的傳播方向,引入了偏移噪聲[5],因此當入射波場與反射波場傳播方向相同時,將產(chǎn)生物理上不存在的成像噪聲,表現(xiàn)為低頻強振幅掩蓋了真實的反射界面[6]。針對RTM成像低頻噪聲問題的處理方法主要可以分為兩類:①減少背向反射波,降低噪聲能量。如,BAYSAL等[7]采用無反射方程來減少背向反射波,但該方法一般只對垂直入射波應(yīng)用效果明顯;LOEWENTHAL等[8]采用大于波長長度的窗函數(shù)對模型慢度做平滑,減少了背向反射波;FLETCHER等[9]通過在波動方程中引入背向反射波衰減項壓制偏移噪聲;還有學者通過上下行波分解限制逆時偏移成像條件,對入射波場和反射波場分別成像以有效減少低頻噪聲的產(chǎn)生[10-14],但是該方法需要進行高維平面波分解,計算量巨大,進行多次高維傅里葉變換難以滿足實際生產(chǎn)對工期的要求。②壓制大反射角成像能量。如,YOON等[15]采用Poynting矢量成像條件在成像前進行震源波場和檢波點波場的傳播角度估計,根據(jù)入射線與反射線的夾角進行加權(quán)疊加,達到壓制偏移噪聲的目的;ZHANG等[1]、XU等[16]提出疊后拉普拉斯濾波等價于對逆時偏移角度道集進行大角度切除,在偏移完畢后進行拉普拉斯濾波壓制低頻噪聲;杜啟振等[17]從計算效率和處理效果兩個方面對壓制低頻噪聲的各種策略進行了深入分析,指出拉普拉斯濾波是目前工業(yè)應(yīng)用中最經(jīng)濟的方法,波場分解成像是效果最好的方法。

但是,上述兩類方法都有一定的適用條件。如減少背向反射波方法必須進行波動方程改造,增加衰減項或輔助變量,計算效率很低,或者需要光滑的背景速度模型,難以符合復雜的地質(zhì)情況;疊后拉普拉斯濾波方法實現(xiàn)簡單、計算效率高,在實際應(yīng)用中得到了廣泛應(yīng)用,但該方法保幅性較差,同時會損失成像剖面的低頻信息。本文提出一種先利用復數(shù)波場重構(gòu)地震數(shù)據(jù),再通過復數(shù)波場運算實現(xiàn)方向波場分解的成像方法,并利用GPU平臺在幾乎不增加計算量的同時完成波場保幅成像。

1 復數(shù)波場構(gòu)建及波場分解

通常認為采用聲波方程直接模擬得到的全波場進行偏移成像具有成像傾角無限制的優(yōu)勢[1-4],但由于目前地震成像理論的基礎(chǔ)仍然是CLAERBOUT提出的褶積成像條件或反褶積成像條件[18],不是地震反射波成像的充分必要條件,在逆時偏移中除了產(chǎn)生低頻成像噪聲外,還有可能產(chǎn)生虛假成像結(jié)果,因此必須考慮地震波的瞬時傳播方向,選擇真實的地下一次反射波進行成像。

首先,將炮點與檢波點波場進行方向波分解。以上下行波為例,可以將全波場表達為:

(1)

式中:S(t,x)、R(t,x)分別表示震源、檢波點全波場;SD(t,x)、RD(t,x)分別為震源、檢波點的下行波場;SU(t,x)、RU(t,x)分別為震源、檢波點的上行波場。若采用(1)式的波場表達方式,則逆時偏移成像條件可以寫為:

(2)

式中:IDU(x)與IUD(x)為傳播方向相反波場的成像結(jié)果,反映了地下界面成像信息;IDD(x)與IUU(x)為傳播方向相同波場的成像結(jié)果,對最終成像結(jié)果不會產(chǎn)生貢獻,只產(chǎn)生偏移噪聲和假象。因此可以直接舍棄IDD(x)和IUU(x),公式(2)寫為:

(3)

從公式(3)可見,實現(xiàn)方向波分解成像必須同時獲得公式(1)中的所有行波分量。雙程波模擬通常采用有限差分方法,只能獲得全耦合的波場,但是可以通過頻率-波數(shù)域分解,分別獲得行波分量。以公式(1)中的震源波場為例,通過傅里葉變換可以得到震源波場的頻率-波數(shù)域表達式:

(4)

其中,S(ω,k)、SD(ω,k)、SU(ω,k)分別是S(t,x)、SD(t,x)、SU(t,x)的頻率-波數(shù)譜。根據(jù)方向波傳播特性,可以獲得上下行波分量的表達式:

(5)

其中,kz表示深度z方向的波數(shù)。可見,已知S(ω,k),就可以很容易獲得上下行波分量。類似地,可以得到檢波點反傳波場的上下行波場表達式。由(5)式可知,波場的方向由頻率-波數(shù)的符號決定,若直接進行波場分解成像,則需要進行沿時間和深度方向的正、反傅里葉變換,對于三維逆時偏移算法而言,其計算量非常大[19-20]。LIU等[11]、YOON等[15]采用Poynting矢量方法在每個成像時刻估計波場的瞬時傳播方向,對于簡單模型取得了較好的應(yīng)用效果,但是該方法在復雜波場情況下無法正確識別同一位置的多個傳播角度和波峰波谷的瞬時傳播方向,實際應(yīng)用中仍存在較多的噪聲影響。本文提出通過希爾伯特變換得到波場的虛部,將原始波場作為實部,構(gòu)建復數(shù)波場優(yōu)化波場分解的高效計算方法。復數(shù)波場的表達形式如下:

(6)

(7)

分析公式(5)可見,如果能確定頻率ω的符號,則只需要判斷kz的符號即可獲得分解后的波場頻率-波數(shù)譜。將復數(shù)波場的頻率-波數(shù)譜公式(7)代入公式(5),即可得到上下行波的復數(shù)波場表達形式:

(8)

由(8)式可知,只要構(gòu)建了復數(shù)波場,在每次成像條件應(yīng)用前進行一次深度方向的傅里葉變換和反變換,即可獲得上下行波分解后的波場。將分解后的上下行波場分別代入公式(2)成像條件,即可獲得成像結(jié)果。

2 復數(shù)波場逆時偏移實施方案

根據(jù)公式(6)進行復數(shù)波場的構(gòu)建需要進行全波場在時間方向的希爾伯特變換,但是常規(guī)逆時偏移僅僅保存幾個時間步的波場信息,無法進行時間方向的希爾伯特變換。LIU等[11]根據(jù)波動方程的有限差分算子具有線性時不變性,將地震數(shù)據(jù)先進行希爾伯特變換再實施波場延拓,獲得了復數(shù)波場的虛部,但是該方法需要進行兩次波動方程的延拓,計算量和存儲量各增加一倍[13]。

依據(jù)公式(8)對復數(shù)波場上下行波的分解性質(zhì)(kz≥0表示上行波,kz≤0表示下行波),從全波場S(ω,k)出發(fā),先考慮kz的符號再區(qū)分頻率的符號位,同樣可以獲得上下行波分解方程:

(9)

其中,

(10)

將公式(9)所示的上下行波分解方程同時應(yīng)用于震源波場和檢波點波場,并代入公式(3)所示的成像條件,可以得到如下復數(shù)波場成像條件:

(11)

其中,SHz(t,x)和RHz(t,x)分別表示震源波場S(t,x)和檢波點波場R(t,x)在深度z方向的希爾伯特變換。公式(11)所示的復數(shù)波場成像條件可以推廣到左右行波分解成像,只需將希爾伯特變換的方向替換為水平方向即可:

(12)

其中,SHx(t,x)和RHx(t,x)分別表示震源波場S(t,x)和檢波點波場R(t,x)在水平方向的希爾伯特變換。

從公式(11)和公式(12)可以看出,在常規(guī)逆時偏移實施過程中,只要對成像前的震源和檢波點波場進行空間的希爾伯特變換,就可以在不增加存儲量僅僅增加很少計算量的條件下完成不同方向的行波分解成像,具體實施方案如圖1所示。目前逆時偏移技術(shù)在實際應(yīng)用中可借助GPU的強大計算能力來實現(xiàn)[21],希爾伯特變換可以通過空間域褶積[20]在GPU上實現(xiàn)并行計算,因此本文復數(shù)波場逆時偏移方法的計算量與常規(guī)逆時偏移方法基本相當。

圖1 復數(shù)波場逆時偏移實施方案

3 模型數(shù)據(jù)測試

采用三維鹽丘模型數(shù)據(jù)對方法進行了測試。該模型縱橫方向13500m ×13500m,垂直深度方向5000m,數(shù)據(jù)采樣點數(shù)為451×451×501,采樣間隔為30m×30m×10m,等間隔激發(fā)3333炮,檢波器在深度z=0全網(wǎng)格共16×104道接收,子波為主頻20Hz的雷克子波,速度模型和單炮記錄見圖2。從單炮記錄看,直達波、反射波、繞射波豐富,寬方位采集方式對于鹽下隱蔽構(gòu)造的照明更好。

圖2 三維鹽丘模型和單炮記錄a 三維鹽丘速度模型; b 采集的單炮數(shù)據(jù); c 單炮數(shù)據(jù)放大顯示的3個檢波器排列

實際計算使用10塊NVIDIA-GPU K10卡,進行了常規(guī)逆時偏移成像、Laplace濾波、復數(shù)波場逆時偏移成像,結(jié)果如圖3、圖4所示,計算效率如表1所示。由圖3可見,在復雜鹽丘模型中利用本文方法進行波場分解可以分別獲得上下、左右行波的波場。從分解效果來看,分解后的波場(圖3c至圖3f)相比雙程波波場(圖3b)幾乎沒有殘余干擾信息,能夠滿足地震成像的需求。由圖4可見,常規(guī)逆時偏移成像結(jié)果被強低頻噪聲覆蓋(圖4a),難以準確識別地層位置;通過Laplace濾波對常規(guī)逆時偏移剖面進行低頻噪聲壓制后(圖4c),可以清晰地看到地震反射層位置。工業(yè)界逆時偏移軟件通常采用RTM+Laplace濾波這種組合方式進行成像[22],由于Laplace濾波是空間二階導數(shù),與RTM直接輸出的振幅存在數(shù)量級的差異,因此不能進行噪聲提取得到噪聲剖面,只能獲得成像剖面。本文復數(shù)波場逆時偏移方法根據(jù)公式(10)、公式(11)所示的聯(lián)合成像條件實現(xiàn)波場分解成像,不僅可以輸出成像剖面(圖4d),而且能夠輸出噪聲剖面(圖4b)。

圖3 鹽丘模型的波場快照a 鹽丘局部模型; b 雙程波波場; c 上行波波場; d 下行波波場; e 左行波波場; f 右行波波場

圖4 成像剖面及噪聲 a 常規(guī)全波場RTM剖面; b 復數(shù)波場逆時偏移所壓制的噪聲; c Laplace濾波后的RTM剖面; d 復數(shù)波場RTM成像剖面

表1 三維鹽丘模型計算參數(shù)與效率統(tǒng)計

從低頻噪聲的壓制角度來分析,RTM+Laplace濾波的方式在壓制低頻噪聲的同時損失了一些低頻有效信息,尤其是高陡反射界面(鹽丘右側(cè))處損失嚴重。本文復數(shù)波場逆時偏移方法則能夠在有效壓制低頻噪聲的同時,保留低頻有效信號,使鹽丘邊界的高陡構(gòu)造成像更準確,水平同相軸的成像更聚焦,波阻特征更明顯,信噪比更高。圖4b所示的噪聲剖面上沒有反射界面的能量損失,表明本文復數(shù)波場逆時偏移成像方法保幅性更高,并且解決了逆時偏移的低頻噪聲壓制所帶來的信號損失問題。

從表1可見,偏移孔徑設(shè)定為13500m、深度5000m,完全覆蓋了整個計算模型。在保證穩(wěn)定性條件的前提下取1ms延拓步長進行波場模擬,時間長度為8s,延拓步數(shù)為7999步。計算過程中,單個模型參數(shù)的三維體大小為388MB,常規(guī)的RTM方法需要5個模型大小的GPU存儲量,復數(shù)波場RTM方法需要增加一個模型空間用于存儲波場希爾伯特變換的結(jié)果,因此是6個模型大小的存儲量,與常規(guī)RTM相比增加了20%的GPU存儲需求。在計算效率方面,復數(shù)波場RTM單炮計算時間增加了約7.5%,整個工區(qū)3333炮數(shù)據(jù)在10塊GPU上完成,常規(guī)RTM需要6.4d,復數(shù)波場RTM需要6.9d,增加的時間很少,基本可以忽略不計。

4 實際數(shù)據(jù)應(yīng)用

為了驗證本文方法在實際資料處理中的適用性,選擇我國東部地區(qū)某實際地震資料進行了測試(圖5、圖6)。工區(qū)內(nèi)地下構(gòu)造復雜,斷裂發(fā)育,尤其是深大斷層邊界刻畫是該工區(qū)需要重點落實的目標。針對高陡斷面波的成像,逆時偏移方法理論上具有90°傾角的成像能力,但是采用常規(guī)RTM+Laplace濾波方法進行成像時,由于Laplace濾波在壓制低頻噪聲的同時會損失低頻高陡構(gòu)造斷面波能量,低頻高陡構(gòu)造成像效果不佳(圖5a,圖6a)。圖5b、圖6b為本文復數(shù)波場逆時偏移方法直接輸出的成像結(jié)果,相比常規(guī)方法的成像結(jié)果在斷面波成像方面有明顯的優(yōu)勢。實際資料測試結(jié)果表明,本文方法在低頻信號的保護方面具有一定的優(yōu)勢,成像保幅性更好。

圖5 實際資料逆時偏移結(jié)果對比(主測線)a RTM+Laplace濾波后的成像剖面; b 復數(shù)波場逆時偏移成像剖面

圖6 實際資料逆時偏移結(jié)果對比(聯(lián)絡(luò)測線)a 常規(guī)RTM+Laplace濾波后的成像剖面; b 復數(shù)波場逆時偏移成像剖面

5 結(jié)束語

本文研究了復數(shù)波場逆時偏移方法,解決了逆時偏移過程中同時進行波場分解成像的計算及存儲問題,在計算量微幅增加并保證計算精度的條件下,實現(xiàn)了上下、左右行波分別成像,給出了復數(shù)波場逆時偏移在GPU平臺上的實現(xiàn)流程。模型數(shù)據(jù)測試結(jié)果表明,本文方法能在計算量和存儲量少量增加的條件下,實現(xiàn)地震波場的分解成像,在低頻保護方面優(yōu)于傳統(tǒng)的RTM方法,成像保幅性更好。實際資料測試結(jié)果表明,本文方法針對復雜高陡斷裂系統(tǒng)的成像效果優(yōu)于傳統(tǒng)RTM方法,尤其在低頻斷面波成像方面具有較明顯的優(yōu)勢。實際地震資料成像中可以兼顧成像效果與計算效率,采用本文方法進行低頻保護的行波分解成像,提高深層斷裂系統(tǒng)的識別能力,落實構(gòu)造圈閉。