劉定進,張永升,段心標,陶永慧

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

地震成像的目標是估計出地下介質(zhì)的幾何形態(tài)和彈性參數(shù)值,為精確描述油氣藏提供可靠的資料。當前生產(chǎn)中常用的偏移成像技術(shù),包括Kirchhoff射線類方法[1-4]、單程波波動方程類方法[5-8]和雙程波波動方程類方法[9-11],其偏移算子都是波場正演傳播算子的共軛轉(zhuǎn)置,而非逆算子。由于地震數(shù)據(jù)觀測孔徑有限、地下模型復(fù)雜以及地震波帶寬有限,采用常規(guī)偏移技術(shù)僅僅能夠得到比較模糊的構(gòu)造信息,無法對復(fù)雜油氣藏進行精細成像,并且成像振幅不均衡,影響后續(xù)巖性參數(shù)反演的可靠性和精度。最小二乘偏移是一種線性化的地震反演成像技術(shù),在最小二乘反演框架下估計一個最優(yōu)的偏移成像結(jié)果,相當于對常規(guī)偏移結(jié)果作用一個Hessian逆矩陣算子,可以降低或者去除采集照明不均或者采集腳印所導(dǎo)致的偏移假象,校正偏移結(jié)果中存在的振幅誤差,提高地震成像的分辨率。因此,近年來最小二乘偏移成為地震成像領(lǐng)域的研究熱點,相繼發(fā)展了最小二乘Kirchhoff偏移[12]、最小二乘單程波偏移[13]和最小二乘逆時偏移[14-16]等技術(shù)。但是,最小二乘偏移基于數(shù)據(jù)擬合方式逐次迭代尋求最優(yōu)成像值,波場傳播算子不能準確預(yù)測出實際數(shù)據(jù),迭代反演計算量巨大,收斂比較困難,三維最小二乘逆時偏移尤為嚴重。國內(nèi)外學(xué)者為此開展了大量的研究工作,如TANG[17]采用隨機相位編碼方式求取Hessian矩陣,并提出了面向目標的最小二乘偏移方法;GUITTON[18]利用初始偏移結(jié)果和反偏移-再偏移結(jié)果估計非穩(wěn)相濾波器來近似Hessian逆矩陣,提高了最小二乘偏移反演的效果;王彥飛[19]討論了最小二乘偏移的預(yù)條件與正則化方法,提高了反演的穩(wěn)定性;ZENG等[20]針對預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)的不一致性,引入了一個可信度參數(shù)作為反演的預(yù)條件;ZHANG等[21]、DUTTA[22]、DUAN等[23]將互相關(guān)誤差泛函引入到最小二乘偏移成像框架中,降低了最小二乘偏移方法對震源子波和振幅因素的敏感性。為了提高最小二乘偏移的計算效率,DAI等[24]提出了平面波編碼的最小二乘逆時偏移方法;黃金強等[25]將平面波最小二乘偏移思路應(yīng)用于各向異性介質(zhì);田文彬等[26]采用導(dǎo)向濾波算子壓制平面波最小二乘偏移假象,同時提高了收斂效率和成像精度;柯璇等[27]則將GPU加速器應(yīng)用于二維最小二乘逆時偏移算法中,加速效果顯著。

本文針對三維最小二乘逆時偏移的數(shù)據(jù)匹配和計算效率問題,采用基于雙程波傅里葉有限差分的反偏移方法、共軛梯度反演步長三點拋物搜索尋優(yōu)方法和GPU加速偏移/反偏移核函數(shù)計算方法,形成了一套具有實際處理能力的三維數(shù)據(jù)域最小二乘逆時偏移技術(shù)方案,利用實際資料驗證了方法的有效性。

1 方法原理

最小二乘偏移成像是指由已知的地震觀測數(shù)據(jù)、假設(shè)的波傳播方程和背景速度模型對地下介質(zhì)的高波數(shù)擾動進行估計。由于觀測數(shù)據(jù)信息不完整、波傳播方程無法精確描述真實波場傳播、震源子波未知以及觀測數(shù)據(jù)存在噪聲等因素影響,地震反演無法獲得精確的參數(shù)模型?，F(xiàn)代地震反演理論采用貝葉斯估計思想對地震反演過程及反演解進行評價。貝葉斯公式為:

(1)

式中:m為地下參數(shù),d為地震數(shù)據(jù),兩者關(guān)系由波場傳播算子L建立,d=Lm;P(m/d)是后驗概率密度函數(shù),P(d)是數(shù)據(jù)空間的概率密度函數(shù),P(m)是模型空間的先驗概率密度函數(shù),P(d/m)是條件概率密度函數(shù)。假設(shè)P(m)和P(d/m)符合高斯分布,P(d)符合標準的正態(tài)分布(觀測系統(tǒng)是規(guī)則的),那么后驗概率密度可表示為:

(2)

式中,CM表示模型空間的協(xié)方差矩陣,CD表示數(shù)據(jù)空間的協(xié)方差矩陣,mprior表示先驗?zāi)Ｐ托畔?dobs表示觀測數(shù)據(jù)。反演問題的求解可以認為是在P(m)、P(d)、P(d/m)的約束下,使P(m/d)最大化。記代價函數(shù)為:

(3)

那么,可將后驗概率密度最大化問題轉(zhuǎn)化為代價函數(shù)最小化問題,于是反演問題變成了一個非線性的、無約束的最優(yōu)化問題。假設(shè)模型空間、數(shù)據(jù)空間的元素相互獨立,每個元素的方差相等,且不考慮先驗?zāi)Ｐ托畔?則得到經(jīng)典的最小二乘偏移目標函數(shù):

(4)

其物理意義是使得觀測數(shù)據(jù)與預(yù)測數(shù)據(jù)的L2范數(shù)誤差最小。

最小二乘反演問題通常采用數(shù)據(jù)域迭代算法(如最速下降法、共軛梯度法等)求解,在一定假設(shè)條件下由初始值逐步逼近真實解。共軛梯度法最小二乘逆時偏移成像迭代公式為:

mk+1=mk-μk+1dk

(5)

式中,mk、mk+1分別是第k、k+1次迭代后的成像值。步長μk+1、共軛梯度方向dk可分別表示為:

(6)

(7)

此處,gk是泛函梯度,由數(shù)據(jù)殘差偏移成像得到。αk是共軛梯度方向修正因子,可采用下式構(gòu)造:

(8)

數(shù)據(jù)域共軛梯度法最小二乘逆時偏移具體實現(xiàn)流程如圖1所示。

圖1 最小二乘逆時偏移流程

2 傅里葉有限差分反偏移方法

基于雙程波方程的逆時偏移和反偏移是最小二乘逆時偏移算法的核心,通常采用高階差分波場傳播算子實現(xiàn)。然而,由于實際資料成像計算網(wǎng)格比較大,采用有限差分方法存在較為嚴重的數(shù)值頻散,計算得到的反偏移數(shù)據(jù)與實際地震采集數(shù)據(jù)存在比較明顯的頻率差異。最小二乘偏移期望正演模擬算子能夠完全模擬觀測數(shù)據(jù),避免由于算子精度不夠?qū)е碌臄?shù)據(jù)不一致問題。如果無法避免這種問題出現(xiàn),應(yīng)該盡量減少頻散效應(yīng)。為了降低反偏移計算中的數(shù)值頻散,我們可以通過提高有限差分的階數(shù)或者減小計算網(wǎng)格的大小,也可以采用更高精度的波場傳播算法,如偽譜法[28-33]、擬解析法[34]和傅里葉有限差分法[35]等來實現(xiàn)。本文利用傅里葉有限差分法,在時間-空間域與時間-波數(shù)域交替進行波場傳播,從而實現(xiàn)雙程波反偏移正演模擬。該方法的主要思想是利用速度分解,求解如下波場傳播方程[36]:

u(x,t+Δt)+u(x,t-Δt)-2u(x,t)=

(9)

圖2是某工區(qū)實際資料反偏移記錄的局部放大結(jié)果,3個方向的計算網(wǎng)格大小分別是25、25、10m。由圖2可見,采用常規(guī)高階有限差分法反偏移的結(jié)果存在較為嚴重的數(shù)值頻散,在記錄上表現(xiàn)為頻散噪聲嚴重(圖2a),而基于傅里葉有限差分法的反偏移結(jié)果則較好地壓制了頻散噪聲(圖2b)。

圖2 某工區(qū)實際資料反偏移記錄a 常規(guī)有限差分法反偏移; b 傅里葉有限差分法反偏移

3 共軛梯度反演步長尋優(yōu)方法

三維數(shù)據(jù)域最小二乘逆時偏移運算成本非常高,除了因為雙程波方程波場傳播算子計算量大以外,另一個更重要的原因是迭代收斂速度慢。共軛梯度法相對于最速下降法而言,是一種收斂速度更快的反演迭代算法。然而,經(jīng)典的共軛梯度法反演步長計算公式是基于目標泛函梯度為0推導(dǎo)得到的,與實際情況不符,迭代更新效果往往不是很好。本文采用三點拋物曲線擬合步長尋優(yōu)法,假設(shè)誤差泛函在解的鄰域為拋物線型,利用3個試探步長所對應(yīng)的3個點擬合一條拋物線,將拋物線極值點對應(yīng)的步長作為迭代更新的最優(yōu)步長。

假設(shè)3個步長α0、α1、α2位于一個單峰區(qū)間內(nèi),對應(yīng)的目標函數(shù)值依次是Eα0、Eα1、Eα2,于是可以構(gòu)造一個拋物函數(shù)。計算其拋物線最小值對應(yīng)的步長:

(10)

構(gòu)造拋物曲線時,3個試探步長點選擇方法如下:α0=0,Eα0為上一次成像更新后的誤差,α1是經(jīng)典共軛梯度法公式計算出的步長,α2=2α1。如果不滿足Eα2>Eα1,則α2=4α1,逐步增大α2,直到Eα2>Eα1為止。

從上述過程可以看出,步長計算中最少需要兩次試探步長后的反偏移模擬來擬合拋物線,如果所有炮都進行反偏移,則計算量巨大,對于三維最小二乘偏移而言不可接受。為了減少計算量,可以隨機選擇部分炮進行正演模擬和曲線擬合。圖3對比了二維模型數(shù)據(jù)最小二乘偏移迭代中全部炮的誤差曲線和部分炮的誤差曲線,可見兩者總體趨勢一致,說明通過部分炮的誤差來擬合拋物曲線是可行的。

圖3 全部炮與部分隨機炮的殘差曲線

4 GPU加速計算方法

在最小二乘逆時偏移中,無論是高階有限差分逆時偏移算子還是傅里葉有限差分反偏移算子,都存在計算量巨大的問題。GPU加速是應(yīng)對這一問題的有效途徑。

高階有限差分算法的GPU并行方案比較簡單,在Y方向做循環(huán),在X、Z方向做多線程并行計算。逆時偏移中采取震源波場存儲方案[37],僅存儲邊界波場,在記錄波場反傳的同時利用邊界數(shù)據(jù)重構(gòu)震源波場,這大大降低了震源波場的直接存儲量。然而,由于GPU內(nèi)存較小,在每一步波場延拓中需要將邊界波場數(shù)據(jù)傳送到CPU進行存儲。

對于傅里葉有限差分算法的GPU加速,由于每個時刻波場延拓包含正、反兩次三維傅里葉變換,計算效率很低,因此采用GPU計算平臺提供的三維傅里葉變換——cuFFT函數(shù)實現(xiàn)FFT,加速20倍以上。圖4是基于GPU的三維傅里葉有限差分法反偏移流程,整體反偏移計算效率提高12倍以上。

圖4 基于GPU的三維傅里葉有限差分法反偏移流程

5 實際資料應(yīng)用

5.1 西部探區(qū)碎屑巖成像

將本文三維數(shù)據(jù)域最小二乘逆時偏移技術(shù)應(yīng)用于西部某工區(qū)實際地震資料處理。該區(qū)中淺層存在一套碎屑巖儲層,對成像分辨率要求較高。實際數(shù)據(jù)共4.8×104炮,計算網(wǎng)格為30m×30m×10m,成像深度為1×104m,Inline、Xline偏移孔徑分別為6000m、5000m,共完成了3次最小二乘逆時偏移成像迭代。圖5對比了逆時偏移和最小二乘逆時偏移的效果,圖6 對比了逆時偏移和最小二乘逆時偏移成像剖面的頻譜。

圖5 西部某工區(qū)成像剖面對比a 逆時偏移; b 最小二乘逆時偏移

圖6 西部某工區(qū)成像剖面頻譜對比

從圖5、圖6可以看出,最小二乘逆時偏移的分辨率更高,儲層內(nèi)部刻畫更清楚。圖7是基于逆時偏移和最小二乘逆時偏移成像剖面的沿層相干屬性分析結(jié)果,可以看出最小二乘逆時偏移結(jié)果中雁形斷裂特征更加清晰。

5.2 南方山地地震成像

南方某山地地震資料整體信噪比偏低,常規(guī)成像效果不佳。選擇面積為60km2的數(shù)據(jù)進行最小二乘偏移成像試處理,數(shù)據(jù)量為7190炮,計算網(wǎng)格為25m×25m×10m,成像深度為1.2×104m,Inline、Xline偏移孔徑分別為5000m、4000m,共完成5次迭代。圖8對比了逆時偏移和最小二乘逆時偏移效果,可見經(jīng)過最小二乘逆時偏移迭代反演,成像剖面信噪比有一定提升,同相軸連續(xù)性變好,斷層成像更加清晰。

圖7 西部某工區(qū)沿層相干屬性分析對比a 基于逆時偏移結(jié)果; b 基于最小二乘逆時偏移結(jié)果

圖8 南方某山地地震資料成像剖面對比a 逆時偏移; b 最小二乘逆時偏移

6 結(jié)束語

本文從反演方法理論、反偏移方法、共軛梯度步長尋優(yōu)算法、GPU加速方法等幾個方面闡述了一套具有實際處理能力的三維數(shù)據(jù)域最小二乘逆時偏移技術(shù)實施方案,并通過兩個探區(qū)實際地震資料處理驗證了方法的效果和優(yōu)勢,獲得以下結(jié)論和認識:

1) 基于雙程波傅里葉有限差分的反偏移方法能夠有效壓制數(shù)值頻散,提高反偏移模擬數(shù)據(jù)與實際觀測數(shù)據(jù)的匹配程度,有利于改善最小二乘逆時偏移成像效果。

2) 導(dǎo)致數(shù)據(jù)域最小二乘偏移效率低下的一個重要原因是反演迭代次數(shù)多,必須提高共軛梯度步長尋優(yōu)的效果和效率,部分隨機炮誤差曲線擬合和三點拋物搜索是達到這一目的的有效技術(shù)手段。

3) 雙程波波場傳播算子計算量大是制約三維數(shù)據(jù)域最小二乘逆時偏移效率的另一個重要因素,借助GPU并行計算能夠顯著提升三維數(shù)據(jù)域最小二乘逆時偏移的計算效率,提升該技術(shù)實際應(yīng)用的可行性。

4) 最小二乘逆時偏移能夠提高地震成像分辨率,改善同相軸的連續(xù)性,整體成像效果優(yōu)于逆時偏移。

由于三維數(shù)據(jù)域最小二乘逆時偏移計算效率低、實際數(shù)據(jù)匹配困難,目前其規(guī)?；茝V應(yīng)用仍然較為困難。下一步的研究方向是先驗信息的利用、Hessian矩陣的近似等,以逐步打造一個規(guī)?；嵱没淖钚《似品囱莩上窦夹g(shù)。