紀(jì)永禎,朱立華,林正良,胡華鋒

(中國(guó)石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

疊前振幅隨角度變化(AVA)反演技術(shù)是利用疊前地震數(shù)據(jù)的振幅隨角度變化特征,采用反演的手段獲取巖石彈性參數(shù)的技術(shù)[1-4]。相較于疊后反演,由于利用了不同入射角數(shù)據(jù)的信息,可以獲得更多樣的彈性參數(shù)估計(jì),因而在儲(chǔ)層描述中起著至關(guān)重要的作用[5]。但由于利用的疊前地震數(shù)據(jù)與彈性參數(shù)的映射關(guān)系更加復(fù)雜,信噪比更低,因而疊前反演方法需要應(yīng)對(duì)更多難題[6-8]。其中一個(gè)難題就是降低解的非唯一性,通常的解決辦法是在目標(biāo)函數(shù)中加入正則化項(xiàng)或在貝葉斯框架下加入待估計(jì)參數(shù)的先驗(yàn)信息約束,在許多可能的解中選擇一個(gè)最優(yōu)或最符合假設(shè)的解,且在貝葉斯框架下加入待求解參數(shù)的先驗(yàn)信息約束,在一定條件下與正則化約束具有相似的內(nèi)涵,可以推導(dǎo)成相同形式[9]。

基于貝葉斯理論框架的疊前反演研究將似然函數(shù)與先驗(yàn)地質(zhì)信息結(jié)合,利用求解最大后驗(yàn)概率來建立反演的目標(biāo)函數(shù),進(jìn)而反演彈性參數(shù)。其中先驗(yàn)信息的正則約束作用可以較好地解決疊前反演中所存在的非唯一性和病態(tài)問題[9]。在構(gòu)建反演目標(biāo)函數(shù)時(shí),似然函數(shù)的構(gòu)建一般體現(xiàn)反演參數(shù)的正演結(jié)果與地震數(shù)據(jù)的匹配關(guān)系,但先驗(yàn)信息的構(gòu)建并沒有統(tǒng)一的假設(shè)及論述。

BULAND等[10]通過概率統(tǒng)計(jì)分析認(rèn)為彈性參數(shù)(特別是聲波速度)大體符合高斯分布,基于此構(gòu)建出貝葉斯反演方法并被廣泛應(yīng)用。DOWNTON等[11]將長(zhǎng)尾分布作為先驗(yàn)分布,獲得了稀疏的反射系數(shù)估計(jì),并認(rèn)為常規(guī)高斯分布先驗(yàn)(這類先驗(yàn)約束與最小二乘本質(zhì)相同)獲得的結(jié)果比較光滑,不利于高分辨率分析和解釋,發(fā)現(xiàn)柯西分布先驗(yàn)比高斯分布先驗(yàn)分辨率更高。此后,利用“長(zhǎng)尾分布”作為先驗(yàn)受到關(guān)注,其中較為常用的是柯西分布及其變體(如修正柯西分布等)。THEUNE等[12]分析了彈性參數(shù)反射系數(shù)稀疏假設(shè)下的兩個(gè)先驗(yàn)?zāi)Ｐ?發(fā)現(xiàn)可微拉普拉斯分布先驗(yàn)可獲得一個(gè)凸的目標(biāo)函數(shù),而柯西分布先驗(yàn)獲得的是非凸函數(shù),因此可微拉普拉斯先驗(yàn)更具優(yōu)勢(shì)。ALEMIE等[13]將多元柯西分布先驗(yàn)引入貝葉斯反演,這種多元先驗(yàn)分布可以合并模型參數(shù)之間的相關(guān)性,并獲得高分辨率的結(jié)果。隨后,Huber分布、t分布等先驗(yàn)先后被引入到貝葉斯反演中,均獲得了較好的效果[9,14-16]。以上先驗(yàn)均假設(shè)反射系數(shù)是隨機(jī)變量,其分布滿足某一特定分布,然而實(shí)際彈性參數(shù)反射系數(shù)的分布受獨(dú)特的地質(zhì)條件影響,常表現(xiàn)出獨(dú)特的特征,無法用某一特定分布準(zhǔn)確描述;不同彈性參數(shù)的反射系數(shù)往往具有不同的分布,對(duì)擬反演的不同彈性參數(shù)的反射系數(shù)進(jìn)行同時(shí)反演時(shí),采用同一分布進(jìn)行約束也不夠合理;并且,特定分布一般由測(cè)井資料統(tǒng)計(jì)獲取,一經(jīng)確定則不隨地震數(shù)據(jù)變化,而地質(zhì)特征的橫向變化引起的反射系數(shù)分布變化難以得到反映。因此,彈性參數(shù)反射系數(shù)的反演結(jié)果可能出現(xiàn)誤差。自動(dòng)相關(guān)性判別(automatic relevance determination,ARD)先驗(yàn)可以解決上述問題,這種先驗(yàn)信息不再假設(shè)反射系數(shù)是符合特定分布的隨機(jī)變量,而是給每一個(gè)反射系數(shù)分配一個(gè)參數(shù),通過估計(jì)該參數(shù)獲得對(duì)反射系數(shù)的最終估計(jì)。該參數(shù)作為一個(gè)待估計(jì)的特定值,由迭代學(xué)習(xí)自適應(yīng)確定,并隨著不同地震道的特征和不同彈性參數(shù)反射系數(shù)的特征變化。這種先驗(yàn)會(huì)保持反射系數(shù)的稀疏性[17-18],從而使反演結(jié)果具有較高的分辨率,其在文字處理、圖像分類和動(dòng)態(tài)光散射分析的應(yīng)用中顯現(xiàn)了較好的潛力[19-21]。

本文通過引入ARD先驗(yàn),將彈性參數(shù)反射系數(shù)的先驗(yàn)信息融合到反演中,降低反演的不適定性和非唯一性。為了降低噪聲對(duì)彈性參數(shù)反射系數(shù)的振幅保持的影響,在反演過程中,引入低頻趨勢(shì)約束。并且,考慮到時(shí)間域疊前反演(利用地震數(shù)據(jù)的所有頻率成分)常受數(shù)據(jù)中5～100Hz以外低信噪比信息的影響,采用了頻域疊前反演策略。另外,考慮到常規(guī)基于彈性阻抗的疊前反演方法盡管具備較為穩(wěn)定的優(yōu)勢(shì)[22-23],但反演過程中需要先獲取分角度彈性阻抗估計(jì),再轉(zhuǎn)換為彈性參數(shù)的估計(jì),會(huì)存在一定的轉(zhuǎn)換誤差,因此,采用同時(shí)反演策略,在反演過程中利用各角度信息直接同時(shí)獲取不同彈性參數(shù)的估計(jì),可以減少?gòu)椥宰杩罐D(zhuǎn)換為彈性參數(shù)時(shí)的轉(zhuǎn)換誤差。

1 方法原理

具有不同入射角度的地震疊前數(shù)據(jù)可以看成是地層界面處的反射系數(shù)與不同角度子波的褶積,地層界面處的反射系數(shù)隨入射角度的變化規(guī)律可以用Zoeppritz方程及其線性近似式描述[24],其中Fatti近似式如下[25]:

(1)

式中:RPP(θ)是隨角度變化的地層界面反射系數(shù),由入射角和彈性參數(shù)反射系數(shù)決定;RIP,RIS,Rρ分別是縱波反射系數(shù)、橫波反射系數(shù)和密度反射系數(shù);γ是橫波速度與縱波速度的比值;θ是入射角。考慮到大入射角地震數(shù)據(jù)較難獲取,導(dǎo)致密度反演不準(zhǔn)確,進(jìn)而會(huì)影響到縱橫波阻抗的反演穩(wěn)定性,且Fatti公式可以在省略密度項(xiàng)的同時(shí)保持對(duì)地層反射系數(shù)的正演精度[25],因此本文將采用兩項(xiàng)Fatti公式作為正反演中采用的近似式。

將兩項(xiàng)Fatti公式寫成矩陣形式:

RPP(θ)=[A(θ)B(θ)][RIPRIS]T

(2)

式中:A(θ)=(1+tan2θ)/2,B(θ)=-4γ2sin2θ。在等式兩邊同時(shí)進(jìn)行傅里葉變換,可以獲得單一入射角疊前數(shù)據(jù)在頻率域的正演矩陣形式:

(3)

式中:

如果采用N個(gè)入射角的疊前地震數(shù)據(jù),即θ∈[θ1,θN],那么正演矩陣形式如下:

(4)

式中:

其中,diag[·]代表對(duì)角矩陣。

為了提高反演過程的穩(wěn)定性,在以上正演矩陣中加入趨勢(shì)約束,形式如下:

(5)

式中:Ψ=F-1ΛFC∈RK×K,Λ=diag[H],Λ為以漢寧窗函數(shù)向量(H)作為對(duì)角元素的對(duì)角矩陣,C是積分矩陣,F和F-1分別為正、反傅里葉變換矩陣。趨勢(shì)約束矩陣Ψ的物理意義是,首先利用積分矩陣C將彈性參數(shù)反射系數(shù)轉(zhuǎn)換成對(duì)應(yīng)的阻抗對(duì)數(shù),然后利用正反傅里葉變換矩陣F和F-1以及漢寧窗矩陣Λ將低頻趨勢(shì)提取出來,與模型趨勢(shì)進(jìn)行匹配約束[3]。LowIP=F-1ΛF[log(IP)]∈RK×1,LowIS=F-1ΛF[log(IS)]∈RK×1,分別是從模型中獲取的縱波阻抗和橫波阻抗的趨勢(shì),IP和IS分別為縱波阻抗和橫波阻抗;λ1,λ2分別是趨勢(shì)約束的權(quán)重系數(shù)。

為便于公式推導(dǎo)和理解,將公式(5)寫成如下簡(jiǎn)單形式:

dMN+2K=G(MN+2K)×2Km2K+n

(6)

式中:d代表公式(5)中等號(hào)左邊的項(xiàng);G包含傅里葉變換矩陣、角度相關(guān)矩陣和低頻趨勢(shì)提取矩陣;m代表待估計(jì)的縱橫波阻抗參數(shù)反射系數(shù)向量;n代表噪聲向量。

在貝葉斯框架下,假設(shè)噪聲向量符合均值為0,方差為σ2的高斯分布,那么地震數(shù)據(jù)的似然函數(shù)有如下形式[26]:

p(d|m,σ2,θ)=(2πσ2)-MNexp·

(7)

引入ARD先驗(yàn)信息,形式如下[18,27]:

(8)

式中:h=[h1,h2,…,h2K]T,其包含2K個(gè)獨(dú)立參數(shù),這些參數(shù)決定了每一個(gè)時(shí)間采樣點(diǎn)處的反射系數(shù)的幅值,物理意義是每一個(gè)反射系數(shù)mk均服從零均值高斯分布,但每一個(gè)反射系數(shù)的方差均不相同,由獨(dú)立參數(shù)hk表示。TIPPING[27]和WIPF等[18]證實(shí)了該先驗(yàn)可以有效保持待估計(jì)彈性參數(shù)反射系數(shù)的稀疏特征。根據(jù)貝葉斯理論,待估計(jì)參數(shù),即彈性參數(shù)反射系數(shù)的后驗(yàn)概率有如下形式:

(9)

Σ=(H+σ-2GTG)-1

(10)

μ=σ-2ΣGTd

(11)

其中,C為常數(shù),H=diag(h1,h2,…,h2K)是一個(gè)對(duì)角矩陣,μ是待估計(jì)彈性參數(shù)反射系數(shù)的后驗(yàn)概率均值,即由彈性參數(shù)反射系數(shù)組成的反演結(jié)果向量。由公式(10)和公式(11)可見,為了獲取μ的估計(jì),需要估計(jì)H和σ2。本文利用輸入的疊前地震數(shù)據(jù)估計(jì)方差σ2;采用第二型最大似然估計(jì)方法,將(12)式的邊緣似然函數(shù)最大化獲得H的估計(jì)。

p(d|h,σ2,θ)=(2π)MN|Q|exp(dTQ-1d)

(12)

式中:Q=σ2I+GH-1GT。最大化邊緣似然函數(shù)的求解方法可以分為3類,具體可以參考TIPPING[27],MACKAY[28]和YUAN等[29]文獻(xiàn)。最后,將獲得的H和σ2的估計(jì)代入公式(11)獲得μ,即彈性參數(shù)反射系數(shù)[RIP,RIS]T的反演結(jié)果。

獲得了縱波阻抗和橫波阻抗反射系數(shù)的估計(jì)后,采用公式(13),可以獲得縱波阻抗和橫波阻抗的估計(jì)[3]。

(13)

式中:IP0和IS0分別是縱波阻抗和橫波阻抗的初值。

從方法的原理介紹中可以得出:引入的ARD先驗(yàn),將反射系數(shù)先驗(yàn)信息融合到反演中,不再假設(shè)反射系數(shù)是符合特定分布的隨機(jī)變量,而是給每一個(gè)反射系數(shù)分配一個(gè)參數(shù),通過估計(jì)該參數(shù)獲得反射系數(shù)的最終估計(jì)。該方法對(duì)反射系數(shù)的估計(jì)不依賴初始模型,完全由地震數(shù)據(jù)驅(qū)動(dòng)獲得。在反演過程中加入的趨勢(shì)約束雖然在一定程度上依賴于初始模型,但也可以從偏移速度等數(shù)據(jù)中獲取,使得反演結(jié)果不會(huì)過多地受到井?dāng)?shù)據(jù)影響。采用在頻率域進(jìn)行反演的策略,可以通過選取優(yōu)勢(shì)頻帶避開疊前數(shù)據(jù)中的低信噪比成分,提高反演的準(zhǔn)確度。由疊前地震數(shù)據(jù)直接獲得彈性參數(shù)的估計(jì),避免了傳統(tǒng)基于彈性阻抗的反演方法中,從疊前數(shù)據(jù)到彈性阻抗再到彈性參數(shù)的轉(zhuǎn)換誤差。

2 模擬數(shù)據(jù)測(cè)試

設(shè)計(jì)的彈性參數(shù)模型和采用的趨勢(shì)約束(參數(shù)模型的0～6Hz成分)如圖1a所示,基于該模型,利用公式(1)獲得了不同入射角的反射系數(shù),并與一個(gè)主頻為30Hz,采樣間隔為1ms的雷克子波進(jìn)行褶積得到模擬的不同入射角的疊前數(shù)據(jù),如圖1b所示。加入20%的隨機(jī)噪聲得到含噪聲的模擬數(shù)據(jù),如圖1c所示。圖1d和圖1e分別給出了未加趨勢(shì)約束的無噪聲和含噪聲數(shù)據(jù)反演結(jié)果;圖1f和圖1g分別給出了加趨勢(shì)約束的無噪聲和含噪聲數(shù)據(jù)反演結(jié)果。反演時(shí),分別選取5～120Hz和10～70Hz作為無噪聲和含噪聲數(shù)據(jù)的頻帶輸入。由圖1d和圖1f可見,在無噪聲情況下,是否加入趨勢(shì)約束并未對(duì)反演結(jié)果產(chǎn)生過多影響,反演結(jié)果均非常好地重現(xiàn)了模型特征。圖1g的反演結(jié)果優(yōu)于圖1e,體現(xiàn)了趨勢(shì)約束具有一定的抗噪性。加趨勢(shì)約束的反演結(jié)果(圖1f、圖1g)較好地重現(xiàn)了模型特征,展示了趨勢(shì)約束的有效性。由模擬數(shù)據(jù)的反演實(shí)驗(yàn)可見,本文提出的方法具有可行性和抗噪性;ARD先驗(yàn)使得反演結(jié)果(無論加入趨勢(shì)約束與否)分界面清晰,“塊化”效果明顯,具有保持邊界的優(yōu)勢(shì),可提供分辨率更高,更有利于后續(xù)儲(chǔ)層和地質(zhì)解釋的反演結(jié)果;趨勢(shì)約束的加入可以有效減小由于隨機(jī)噪聲引起的彈性參數(shù)反射系數(shù)的幅值畸變,有利于反演結(jié)果保持穩(wěn)定并且使反演結(jié)果的準(zhǔn)確性更高。

圖1 無噪和含噪模擬數(shù)據(jù)對(duì)本文方法的測(cè)試結(jié)果a 彈性參數(shù)模型和趨勢(shì)約束; b 無噪聲模擬數(shù)據(jù); c 含噪聲模擬數(shù)據(jù); d 未加趨勢(shì)約束的無噪聲數(shù)據(jù)反演結(jié)果; e 未加趨勢(shì)約束的含噪聲數(shù)據(jù)反演結(jié)果; f 加趨勢(shì)約束的無噪聲數(shù)據(jù)反演結(jié)果; g 加趨勢(shì)約束的含噪聲數(shù)據(jù)反演結(jié)果

采用經(jīng)典的推覆體2D模型進(jìn)一步測(cè)試本文方法,利用30Hz雷克子波并加入10%的隨機(jī)噪聲獲得了含噪聲的、采樣間隔為1ms的不同角度模擬數(shù)據(jù)(入射角度分別為10°、20°和30°),模擬數(shù)據(jù)的縱橫波阻抗參數(shù)模型如圖2a和圖2b所示;利用ARD先驗(yàn),但未加入趨勢(shì)約束(參數(shù)模型的0～6Hz成分)的反演結(jié)果如圖2c和圖2d所示;加入趨勢(shì)約束的反演結(jié)果如圖2e和圖2f所示?？傮w上看,反演獲得的縱波阻抗結(jié)果優(yōu)于橫波阻抗結(jié)果;趨勢(shì)約束的加入提高了反演結(jié)果的橫向連續(xù)性,并且減輕了隨機(jī)噪聲引起的“掛面條”現(xiàn)象,這一提高對(duì)于穩(wěn)定性相對(duì)較差的橫波阻抗反演結(jié)果更為明顯。2D模型的測(cè)試結(jié)果進(jìn)一步證實(shí)了本文方法在獲得高分辨率反演結(jié)果和保持反演結(jié)果橫向穩(wěn)定性方面的有效性和優(yōu)勢(shì)。

圖2 2D推覆體模型模擬數(shù)據(jù)對(duì)本文方法的測(cè)試結(jié)果a 縱波阻抗參數(shù)模型; b 橫波阻抗參數(shù)模型; c 未加趨勢(shì)約束的縱波阻抗反演結(jié)果; d 未加趨勢(shì)約束的橫波阻抗反演結(jié)果; e 加趨勢(shì)約束的縱波阻抗反演結(jié)果; f 加趨勢(shì)約束的橫波阻抗反演結(jié)果

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

實(shí)際數(shù)據(jù)來自川西地區(qū),該地區(qū)目的層發(fā)育有多期河道儲(chǔ)層。首先,選取了一個(gè)井旁道地震數(shù)據(jù)作為方法的實(shí)驗(yàn)數(shù)據(jù)。不同入射角的疊前井旁道地震數(shù)據(jù)如圖3a所示,采用的入射角度分別為5°、15°和25°。反演時(shí)選取10～70Hz作為輸入數(shù)據(jù)的頻帶范圍。對(duì)應(yīng)的趨勢(shì)約束曲線(彈性參數(shù)測(cè)井曲線的低頻成分)、彈性參數(shù)測(cè)井曲線與反演結(jié)果如圖3b和圖3c 所示。由圖3b和圖3c可見,反演結(jié)果展示了與測(cè)井曲線較好的匹配關(guān)系,并且反演結(jié)果具有較好的邊界特征,較好地反映了地層彈性參數(shù)的變化特征,其中藍(lán)色虛線框中較薄的儲(chǔ)層也被較好地展現(xiàn)出來,有利于儲(chǔ)層識(shí)別和后續(xù)地質(zhì)解釋。盡管在180～220ms薄層的連續(xù)變化特征沒有被完美地重現(xiàn),反演結(jié)果仍然對(duì)后續(xù)的儲(chǔ)層和地質(zhì)解釋具有重要意義。

圖3 實(shí)際井旁地震道數(shù)據(jù)反演結(jié)果a 不同入射角井旁地震道數(shù)據(jù); b 疊合了測(cè)井曲線和趨勢(shì)約束的縱波阻抗反演結(jié)果; c 疊合了測(cè)井曲線和趨勢(shì)約束的橫波阻抗反演結(jié)果

選取實(shí)際數(shù)據(jù)中某二維剖面作為實(shí)驗(yàn)數(shù)據(jù)。不同入射角的疊前部分疊加地震數(shù)據(jù)如圖4所示,采用的入射角分別為7°、14°和21°。采用的低頻趨勢(shì)約束由提取測(cè)井彈性參數(shù)曲線插值模型的低頻成分(0～6Hz)獲得。通過井震標(biāo)定獲得不同角度子波進(jìn)行反演。反演時(shí)選取10～70Hz作為輸入數(shù)據(jù)的頻帶范圍,反演后將結(jié)果中高于150Hz的成分濾除。圖5a和圖5b分別給出了采用本文方法反演得到的縱波阻抗和橫波阻抗。由圖5a和圖5b可見,本文方法得到的反演結(jié)果分辨率較高,地層界面清晰,與地質(zhì)認(rèn)識(shí)較為吻合,為后續(xù)儲(chǔ)層和地質(zhì)解釋提供了較好的基礎(chǔ)資料。圖5c和圖5d分別給出了采用某商業(yè)軟件得到的縱波阻抗和橫波阻抗反演結(jié)果。由圖5可見,本文方法具有更高的分辨率,尤其是橫波阻抗的反演結(jié)果優(yōu)勢(shì)明顯,具有更好的可解釋性,有利于后續(xù)儲(chǔ)層與地質(zhì)解釋。不足的是,盡管橫波阻抗反演結(jié)果的改善很大,但剖面局部仍存在一些由于噪聲和振幅畸變?cè)斐傻牟环€(wěn)定現(xiàn)象,后續(xù)研究可以重點(diǎn)嘗試?yán)枚嗟啦呗约訌?qiáng)反演結(jié)果的橫向穩(wěn)定性。本文方法的計(jì)算效率在可接受的范圍內(nèi),但略遜于商業(yè)軟件,后續(xù)研究可以進(jìn)行針對(duì)性的算法優(yōu)化以提高效率。

圖4 具有不同入射角的疊前部分疊加地震數(shù)據(jù)a 入射角為7°; b 入射角為14°; c 入射角為21°

圖5 實(shí)際地震數(shù)據(jù)采用本文方法和某商業(yè)軟件得到的橫波阻抗和縱波阻抗反演結(jié)果a 本文方法反演得到的縱波阻抗; b 本文方法反演得到的橫波阻抗; c 某商業(yè)軟件反演得到的縱波阻抗; d 某商業(yè)軟件反演得到的橫波阻抗

4 結(jié)論

通過引入自動(dòng)相關(guān)判別先驗(yàn),在頻率域提出了一種新的疊前同時(shí)反演方法。自動(dòng)相關(guān)判別先驗(yàn)假設(shè)每個(gè)反射系數(shù)都由一個(gè)待估計(jì)的參數(shù)決定,避免了對(duì)反射系數(shù)的總體分布進(jìn)行假設(shè),并保證了反演結(jié)果的較高分辨率和對(duì)巖性分界面的識(shí)別能力。同時(shí),該參數(shù)可以隨著不同地震道和彈性參數(shù)反射系數(shù)的特征自適應(yīng)變化,提高了反演的準(zhǔn)確性。反演中加入的趨勢(shì)約束有效提高了反演的穩(wěn)定性和對(duì)噪聲的適應(yīng)能力。本文方法利用優(yōu)勢(shì)頻帶獲取反演結(jié)果,避免了利用信噪比低的頻帶信息,有利于提高反演結(jié)果的準(zhǔn)確性。模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)實(shí)驗(yàn)結(jié)果證實(shí)了本文方法的可行性、優(yōu)勢(shì)和應(yīng)用潛力。需要注意的是,基于離散傅里葉變換的低頻約束會(huì)在一定程度上降低反演的效率,在實(shí)際應(yīng)用中,可以考慮采用頻譜分析縮減反演利用的頻率區(qū)間、精細(xì)劃分層位控制參與反演的采樣點(diǎn)數(shù)等方面進(jìn)一步優(yōu)化算法、提高效率。