陳樹民,劉禮農,張劍峰,裴江云,陳志德,王 成

(1.中國石油天然氣股份有限公司大慶油田有限責任公司勘探開發(fā)研究院,黑龍江大慶163712;2.中國科學院地質與地球物理研究所中國科學院油氣資源研究重點實驗室,北京100029)

地震波傳播過程中,一部分能量會轉化為熱能,相應的地震波的幅值產生衰減效應,該過程稱為吸收。地震波不同頻率成分的幅值衰減是不同的,頻率越高,衰減越快,這是因為高頻成分的波長相對低頻成分的波長較短,對于一個固定的傳播距離,相當于低頻成分少個波長,高頻成分則多個波長,而地震波每傳播一個波長的距離,能量損失的程度恒定。這也導致接收到的反射地震資料的有效頻帶隨反射深度增大逐漸變窄;而不同頻率成分以不同的速度傳播,也導致了地震子波的頻散,且反射層位越深,頻散越嚴重。由于常規(guī)偏移方法沒有補償黏性吸收導致的幅值衰減,也沒有校正頻散,因而使得偏移成像結果的分辨率較低,反射層位越深,分辨率越低。油氣勘探中,對薄砂體、小斷層的識別是認識油氣疏導體系、識別有利儲層的重要環(huán)節(jié)。因此,在油氣勘探中,提高地震成像分辨率一直是反射地震資料處理過程中的一個重要研究內容。

提高地震成像分辨率的方法很多,包括針對疊后剖面的譜白化反褶積[1]、非穩(wěn)態(tài)反褶積[2-3]、基于統(tǒng)計假設或測井資料的各類拓寬頻帶技術以及反Q濾波[4-5]和黏彈性疊前深度偏移[6-8]等。譜白化反褶積可在有效頻帶內將振幅譜拉平,但其提高分辨率的效果較大程度上依賴于是否獲取了合適的參數[1]。各類拓頻技術通過引入地震記錄以外的信息來提高地震方法的分辨率,在獲得更高視分辨率的同時,也需對地震數據基本頻帶進行提高信噪比的處理;此外,各類拓頻技術使用的前提是地震記錄的子波是時不變的,因此即使應用此類技術,也必須在預處理中首先補償地震波的吸收衰減,實現地震子波的一致性[9]。

非穩(wěn)態(tài)反褶積是針對黏性吸收導致的地震分辨率降低而發(fā)展的提高分辨率方法。該方法基于嚴謹的物理基礎,但在估計空變的非穩(wěn)態(tài)子波方面存在較大困難,且通常不能同時實現頻散校正。目前,非穩(wěn)態(tài)反褶積僅在疊加剖面上能得到較好應用效果。穩(wěn)定性和噪聲放大是該方法在實際應用中面臨的問題。

反Q濾波可同時應用于疊前和疊后地震資料。該方法從補償地震波幅值的黏性吸收出發(fā),具有堅實的物理基礎。反Q濾波應用于疊前地震資料時忽略了地震波傳播路徑對幅值衰減的影響,用于非均勻Q值模型時存在較大誤差;應用于疊后地震資料時可處理層狀Q值模型情況,但由于疊加過程已將不同偏移距或者入射角的地震道進行了疊加,此時,應用反Q濾波是將存在不同程度吸收衰減的幅值作同一處理,因而并不能完全消除吸收衰減的影響。

為了補償地震波吸收衰減,提高地震成像分辨率,在偏移成像過程中恢復地震波被衰減的高頻成分是提高地震成像分辨率的關鍵。恢復地震波被衰減的高頻成分可真正地提高地震勘探方法對小斷層和薄砂體的識別能力。

黏彈性疊前深度偏移準確考慮了地震波的傳播和黏性衰減過程,理論上可較好地補償地震波的黏性吸收,但黏彈性深度偏移方法需要深度域地層Q值模型,難度更大,因此黏彈性深度偏移的實際應用還存在較大的困難。

疊前時間偏移是疊前偏移成像中的另一種重要方法。該方法可對斷層較為復雜但速度橫向變化不是很劇烈的地質構造進行較好的成像。與疊前深度偏移方法相比,除具有較高的計算效率外,其主要優(yōu)點是只需使用疊加(均方根)速度,這樣可簡單地通過速度掃描等方式得到合適的偏移速度模型,回避了使用疊前深度偏移方法面臨的層速度建模困難的問題。因此,疊前時間偏移方法已成為地震勘探領域廣泛應用的關鍵技術之一。但現在常用的疊前時間偏移方法不具有補償地震波吸收衰減的能力,所以如何結合疊前時間偏移過程恢復地震波被衰減的高頻成分,成為一個具有重要應用價值的問題。

近年來,通過結合疊前深度偏移和疊前時間偏移方法,人們提出了一系列基于等效參數的廣義疊前時間偏移方法[10-13],在疊加速度的基礎上引入新的等效參數,拓寬了現行疊前時間偏移方法解決實際復雜問題的能力和范圍。通過引入第二個描述黏性吸收的等效Q值參數,提出了可直接補償吸收衰減的疊前時間偏移方法[11-12],用較簡單的技術方案實現了補償地震波傳播過程中吸收衰減的目標。本文詳細闡述了這一結合疊前時間偏移過程恢復地震波被衰減的高頻成分的方法,以及在保持偏移計算穩(wěn)定性、介質吸收參數建模、補償噪聲壓制以及計算效率提高等方面提出的一系列技術方案。

1 補償介質吸收疊前時間偏移

在提高地震數據分辨率的研究中,通常假定品質因子Q與頻率無關或者隨頻率弱變化[14],這便是恒Q模型,在地震學觀測頻率范圍內(0.001～100.000Hz),Q值可以近似看作常數。基于此,引入復速度c(ω)與實速度vr(ω)表達式[15]:

式中:ω是角頻率;ωc是高截頻;j為虛數單位。當頻率趨向于ωc時,實速度將趨近于常值。在此,用主頻ω0來代替高截頻ωc,原因在于主頻對應的實速度恰是我們通過速度估計方法可以得到的。作如下代換和近似:

(2)

將(2)式代入(1)式得:

(3)

式中:v為常規(guī)的偏移速度;Q代表與頻率無關的恒Q值。

若將單個地震道看作是僅有一個接收道的單炮記錄,則其炮域偏移的炮點至成像點的正傳波場為[11]:

(4)

檢波點至成像點的反傳波場為:

(5)

實際上,我們無法知道準確的S(ω),而現行的處理流程中的反褶積處理可以認為是已剔除了S(ω)·[ω/(2π)]exp(jπ/2)這些有關震源子波的影響。因此,忽略震源子波的影響,應用反褶積成像條件,可得單道數據的成像結果,也即脈沖響應,即:

(6)

式中:(τs/τg)2是成像權系數,它補償了地震波的球面擴散影響。公式(6)表明,成像點處的等效Q值唯一決定該點的成像幅值。因此,我們可以采用類似于確定均方根偏移速度的方法,用掃描方式來確定該點的等效Q值。

(7)

式中:f′(t)為該地震道的一階導數。

當Qeff不隨時間深度和水平位置變化時,若假設:

(8)

則(6)式可寫為:

(9)

由(8)式可知,G(t)恰好是該地震道的反Q濾波結果。(9)式表明,“反Q濾波+疊前時間偏移”是補償吸收衰減的偏移方法在均勻Q值情況下的特例。這為我們應用掃描方法快速建立等效Q值模型提供了理論基礎。

φ(η)=exp(η),η≤lnG

(10a)

φ(η)=G[1-lnG-2.5(lnG)2]+G(1+5lnG)η-

2.5Gη2, lnG<η≤lnG+0.2

(10b)

φ(η)=1.1G,η>ln(G)+0.2

(10c)

式中:η=ω(τs+τg)/(2Qeff);G為給定的閾值。

2 等效Q值估計

Q值估計涉及到隨頻率變化的幅值,與偏移速度估計相比困難得多,目前主要是利用透射波的信息求取,即利用上行波的VSP測井資料或者井間資料,根據主頻移動、頻譜形狀等信息確定Q值。在實際應用中,其局限性很明顯:通常情況下,地震勘探目標工區(qū)的井資料有限,難以建立非均勻的Q值模型,此外,Q值是一個與地震信號主頻相關的變量,將依據透射資料得到的Q值模型應用到反射地震資料時,因為地震反射波與透射波主頻不同,還需做進一步校正。為此,我們發(fā)展了利用反射資料根據補償效果和頻帶寬度確定等效Q值的方法。

由(9)式可知,“反Q濾波+疊前時間偏移”是補償吸收衰減的疊前時間偏移方法在均勻Q值情況下的特例,而其計算量僅是后者的幾十分之一,我們只需對預先給定的等效Q值序列分別進行“反Q濾波+常規(guī)疊前時間偏移”計算,就可得到進行等效Q值掃描的基礎數據集。而在實際地震處理中確定合適的等效Q值時,往往是在信噪比可接受的前提下追求頻譜的最寬,因此,在每一分析時窗內,我們以頻譜的最寬化作為等效Q值選取的準則。

具體來講,在每個等效Q值拾取窗口,對系列Q值中的Qi,從其成像剖面最終結果中計算頻譜包絡,拾取-20dB對應的最大和最小頻率f1a和f1b(工業(yè)中通常認為-20dB以下的信號不再有效,這一分貝值在軟件實現時可作為參數由用戶設置)。為了防止-20dB遇到頻譜凹陷,采用分貝區(qū)間[d1,d2]代替-20dB點。區(qū)間[d1,d2]包含-20dB,d1和d2圍繞-20dB上下分布,如果頻譜上下震蕩比較劇烈則d1和d2波動也較大,將區(qū)間等間距的分為N份,那么第i個點值是d1i=d1+(i-1)×(d2-d1)/(N-1),此時-20dB對應的最大、最小頻率即為di對應的頻率平均值。

(11)

式中:f1ai,f1bi分別表示d1i對應的最大、最小頻率。同樣,設f2a和f2b分別表示-10dB對應的最大和最小頻率,用區(qū)間[d3,d4]代替-10dB點,區(qū)間等間距分為N份,第i個點值為d2i=d3+(i-1)×(d4-d3)/(N-1)。那么-10dB對應的最大、小頻率即為di對應的頻率平均值。

(12)

式中:f2ai,f2bi分別表示d2i對應的最大、最小頻率。

從小偏移距的部分偏移疊加剖面的頻譜包絡中,拾取-10dB對應的最大和最小頻率f2a1和f2b1,從大偏移距偏移疊加剖面的頻譜包絡中,拾取-10dB對應的最大和最小頻率f2a2和f2b2,計算頻率偏差;基于f1ai,f1bi和f2ai,f2bi得到4條隨1/Q變化的曲線,首先厘定較小的區(qū)間,在該區(qū)間中,主要觀察-20dB對應的高、低頻曲線,尋找將高頻提高得最多同時低頻損失較少的那個Q值,這樣才能確定Q值達到了提高有效頻帶的目標,最終確定該Q值窗中心處的Q值。如圖1所示,這是某一窗口的高、低頻曲線,兩條紅線(或藍線)之間的寬度就是-20dB(或-10dB)對應的頻譜寬度。圖1中,綠線對應的橫坐標表示所確定的Q值,此時高頻較高,同時低頻損失較少。

圖1 某一分析時窗不同分貝下的頻帶曲線 圖中,紅線為-20dB對應的最大、最小頻率隨1/Q變化的曲線，藍線為-10dB對應的最大、最小頻率隨1/Q變化的曲線。

3 結合穩(wěn)相偏移方法壓制偏移噪聲

針對偏移噪聲壓制,基于菲涅耳帶疊加實現偏移成像就是一個最佳選擇,這也是各類偏移方法一直努力的目標[16-18]。由于實際地質構造中速度和反射界面的復雜性,直接從地質模型出發(fā)估計準確的菲涅耳帶幾乎不可能實現。針對這一問題,我們在偏移過程中構建傾角域偏移道集,直接將菲涅耳帶形象地展示在這一偏移道集中,可容易地從中確定菲涅耳帶。

3.1 傾角道集

如圖2所示,I為成像點,s和g分別代表炮點和接收點,v和vrms分別為成像點I處的層速度與均方根速度,n為反射界面法向量,p和q分別為成像點I處的入射向量和反射向量,dx和dy分別為反射界面與xOz平面和yOz平面的交線。定義反射界面在xOz平面的視傾角為θx,其旅行時相關的視傾角為φx;反射界面在yOz平面的視傾角為θy,其旅行時相關的視傾角為φy,則有:

圖2 傾角道集計算示意

若在每一成像點,不考慮偏移距,僅將偏移結果按φx和φy的大小進行分選和疊加,就可在每個水平位置(即CRP點)形成一對分別針對沿測線方向傾角和與測線垂直方向傾角的一維傾角域偏移道集I(x,y,T0,φx),I(x,y,T0,φy):

圖3 某CDP點的inline方向(a)與crossline方向(b)的傾角道集(圖中白線為拾取的菲涅耳帶邊界)

3.2 結合穩(wěn)相偏移壓制噪聲

在三維補償介質吸收疊前時間偏移計算中,通過計算擬成像反射界面的傾角來判斷傾角是否在菲涅耳帶內,若在菲涅耳帶內,則參與穩(wěn)相疊加,從而實現穩(wěn)相疊前時間偏移。這樣,可保證僅是菲涅耳帶內的結果參與偏移成像。如此,既可避免常規(guī)偏移中采用較小的偏移孔徑導致的陡傾角構造成像缺失,又可避免采用過大的偏移孔徑帶來的偏移成像低信噪比問題,實現基于菲涅耳帶的最優(yōu)孔徑成像[12]。

與常規(guī)疊前時間偏移相比,穩(wěn)相疊前時間偏移在對每一地震道的偏移計算中,先利用傾角公式((13)式和(14)式)計算在每一成像點處的與旅行時相關的傾角,同時讀取該成像點的菲涅耳帶在兩個傾角方向上的上、下界角度。當該地震道在這一成像點處的與旅行時相關的傾角在菲涅耳帶內時,將計算結果累加到存放成像點I(x,y,T0)上,否則不參與成像,這樣可得到僅對傾角域偏移道集中的菲涅耳帶部分疊加成像的高信噪比偏移結果。圖4給出了未結合和結合穩(wěn)相算法的補償介質吸收疊前時間偏移成像的水平切片,可以看出,穩(wěn)相偏移水平切片在保持構造形態(tài)不變的情況下信噪比得到很大提高。此外,需要指出的是,穩(wěn)相偏移不但可以實現合理的偏移孔徑空變以得到復雜構造高精度精細成像結果,而且還可以提高偏移計算的效率。圖4所示的補償介質吸收疊前時間偏移計算實例,結合穩(wěn)相偏移方案可節(jié)約30%左右的計算時間。

圖4 疊前時間偏移的水平切片a 結合穩(wěn)相算法的補償介質吸收; b 未結合穩(wěn)相算法的補償介質吸收

4 基于GPU加速和異步雙緩存的并行計算

補償介質吸收疊前時間偏移是在頻率域內進行,較常規(guī)時間域偏移而言,計算量明顯增加。增加的計算量主要體現在兩個方面:一是頻率域的補償算法包括相位校正、振幅補償,且這些計算大多都是對數和指數運算,對應多個計算機指令,耗時較長;二是與常規(guī)疊前時間偏移相比,這一偏移算法在頻率域進行,幅值信息必須將所有頻率成分分別計算得到的結果累加。因此,我們基于GPU技術提高其計算效率[19-21]。

由于實際三維地震資料數據量龐大,地震道數多,采用前文所述的單道成像方案會使數據在CPU和GPU之間頻繁拷貝,而且數據處理的前置準備工作也不可避免,其中包括數據從硬盤讀入內存、導數求取、能量調節(jié)等,針對單個地震道而言,這些環(huán)節(jié)和核心成像過程必須串行,并行效率低。為此,我們提出采用多流技術的異步雙緩存方案,該方案是在原始單道成像基礎上,采用單次批量進入多道數據,多條成像線結果輸出。這樣在GPU進行成像計算時,由于異步執(zhí)行,所以CPU是就緒狀態(tài)的,可以準備第二道數據的前置計算,這樣就使得CPU的前置工作和實際GPU成像(頻率域累加)并行執(zhí)行。對用戶而言,隱藏了CPU的工作時間,使得GPU可以充分發(fā)揮其計算性能,將必要的CPU工作掩蓋在GPU的工作過程中。

由于拷貝和核函數的異步執(zhí)行可能會導致程序中不同流之間成像孔徑的計算拷貝及成像計算產生錯誤,因此,我們提出了雙緩存方案。通過增加一個變號器,在每一次異步拷貝時用變號器控制成像孔徑拷貝,保證了隨后成像計算的正確進行。

5 大慶油田實際應用

為推動補償介質吸收疊前時間偏移技術的工業(yè)應用,我們集成上述理論研究成果,開發(fā)了一套大規(guī)模實用化的地震疊前高分辨率成像GPU軟件系統(tǒng),解決了應用中面臨的參數選擇、數據多樣性、資料噪聲、海量數據存取與計算以及魯棒性等問題。該偏移成像方法已在大慶油田廣泛應用,及時有效地指導了大慶油田扶余油層致密油的勘探部署。

大慶油田扶余油層以陸相河流—三角洲沉積體系為主,埋深超過1500m,砂體厚度薄(單層厚度2～5m),橫向速度變化快(河道寬度小于600m),扶余油層致密油水平井勘探主要針對2～5m單砂體。上覆地層為砂泥巖薄互層,地層壓實程度低、非均質性強。上覆地層對地震波的吸收是影響地震資料分辨率和保真度的關鍵因素。

圖5為永樂工區(qū)1085線常規(guī)疊前時間偏移剖面和補償介質吸收疊前時間偏移剖面。圖5中方框區(qū)域放大結果如圖6所示。可以看出,本文方法有效消除了地層吸收衰減的影響,提高了地震資料的縱向分辨率,從而提高了該區(qū)地震資料刻畫薄儲層的能力。圖7 為放大結果對應的頻譜,可以看出,補償介質吸收疊前時間偏移在振幅-20dB處頻譜展寬20Hz。

圖5 常規(guī)疊前時間偏移(a)和補償介質吸收疊前時間偏移(b)剖面對比

圖6 常規(guī)疊前時間偏移(a)與補償介質吸收疊前時間偏移剖面(b)局部放大顯示

圖7 常規(guī)疊前時間偏移與補償介質吸收疊前時間偏移剖面局部放大結果對應的頻譜分析曲線

大慶油田三肇地區(qū)扶余油層垂向上在T2反射層以下100ms之內,是目前致密油勘探水平井部署的重點區(qū)域。圖8是針對扶余油層3～5m甜點目標部署的水平井ZP15井的常規(guī)疊前時間偏移剖面與補償介質吸收疊前時間偏移剖面。部署的主要參考直井為Z9-38井和Z57-F1井,甜點層在剖面1.48s的位置,地震反射特征對應中強反射軸。在水平井軌跡的位置(A,B,C,D分別代表水平井軌跡4個控制點,A點為入靶點)。常規(guī)疊前時間偏移剖面上,扶余油層的頻帶為6～70Hz,垂向分辨率低,不能識別位于T2反射下波谷內的F11-3油層組目標砂巖。補償介質吸收疊前時間偏移剖面上扶余油層的頻帶為6～90Hz,垂向分辨率提高,F11-3油層組目標砂巖對應于T2反射下的第一個波峰位置,可以連續(xù)追蹤。振幅橫向變化能夠反映砂體厚度的橫向變化。依據補償介質吸收疊前時間偏移成像資料,指導了ZP15井水平井鉆探,ZP15井在A靶點準確入靶,水平段長度1050m,其中,砂巖1010m,泥巖40m,砂巖鉆遇率達93.16%。

補償介質吸收疊前時間偏移成像方法已在大慶油田累計完成16個三維井區(qū)合計590km2的目標處理,提交16口水平井部署設計,已實施的9口水平井均準確入靶,水平段平均砂巖鉆遇率達到86.0%,比前期應用常規(guī)偏移結果提高了15.1%,油層鉆遇率提高了14.4%,有效地支撐了大慶油田扶余油層致密油勘探。

圖8 常規(guī)疊前時間偏移(a)和補償介質吸收疊前時間偏移(b)過ZP15井剖面

6 結論

本文介紹了補償介質吸收疊前時間偏移技術及其在大慶油田的應用實踐。等效Q值參數的引入和對應的基于地表觀測數據的等效Q值建模方法,解決了黏性偏移應用中的Q值模型建模的困難;利用黏性補償因子的閾值穩(wěn)定性控制方法,解決了黏性吸收補償的穩(wěn)定性問題。將基于傾角道集的穩(wěn)相(菲涅耳帶疊加)偏移方案結合到黏彈性疊前時間偏移方法中,有效地解決了黏性補償導致的偏移噪聲放大問題;而偏移實現GPU算法的優(yōu)化也較大程度地提高了這一補償算法的計算效率。大慶油田多個三維實際地震資料的處理結果證實,本文方法在保持反射振幅關系的前提下展寬地震頻帶20Hz,有效提高了地震資料的垂向分辨率,是松遼盆地中淺層薄互層地質條件下高分辨率成像的有效技術。