張凌遠(yuǎn),張宏兵,許可欣,祝新益,曾繁鑫,郭 強(qiáng)

(1.河海大學(xué)地球科學(xué)與工程學(xué)院,江蘇南京211100;2.中國礦業(yè)大學(xué)資源與地球科學(xué)學(xué)院,江蘇徐州221116)

近年來,油氣資源的勘探與開發(fā)逐漸向非常規(guī)油氣、海洋深水油氣等轉(zhuǎn)移。由于非常規(guī)油氣的埋藏與分布、品級與質(zhì)量相對于常規(guī)油氣藏均有所區(qū)別[1],因而油氣儲層的準(zhǔn)確識別與預(yù)測十分重要。疊前地震反演作為巖性識別、流體預(yù)測及油氣檢測最重要的技術(shù)之一,經(jīng)歷了長時間的發(fā)展。但疊前地震反演中地震數(shù)據(jù)與模型參數(shù)之間高度非線性特征,導(dǎo)致反演是一個病態(tài)的問題,反演結(jié)果不穩(wěn)定。并且,大部分疊前地震反演是基于各向同性介質(zhì)背景,而實(shí)際地下介質(zhì)情況十分復(fù)雜,它們由厚度不一的地層及各類構(gòu)造組成,具有明顯的各向異性特征,這對確保反演結(jié)果的精度帶來巨大的困難。

為解決反演的不適定問題,國內(nèi)外許多學(xué)者在反演過程中引入正則化加以約束,最廣泛應(yīng)用的是Tikhonov正則化方法[2-3],即在目標(biāo)函數(shù)中添加模型參數(shù)的零階或高階先驗(yàn)項(xiàng)。崔巖等[4]將帶先驗(yàn)知識的正則化引入到共軛梯度法,加快收斂速度,避免了傳統(tǒng)方法迭代不充分或迭代次數(shù)過多的問題。王振濤等[5]采用L2范數(shù)正則化方法并自適應(yīng)地獲取正則化參數(shù),借助Zoeppritz方程獲得雅可比矩陣,進(jìn)行線性化的疊前反演。但這類正則化方法使得反演結(jié)果過于平滑,難以體現(xiàn)地層與構(gòu)造邊界的特征,而Markov隨機(jī)域(MRF)作為一類聯(lián)合概率分布的圖像模型[6-7],通過控制模型梯度達(dá)到懲罰目標(biāo)函數(shù)結(jié)果,具有顯著的增強(qiáng)邊緣特性的能力,被廣泛應(yīng)用于圖像處理及修復(fù)[8]、非線性地球物理反演問題[9-11]等方面。因而,有學(xué)者提出了基于MRF的保護(hù)正則化方法[12-13]。

曹磊等[14]在火成巖預(yù)測中采用疊前地震反演并結(jié)合MRF,能夠直接預(yù)測火成巖的物理性質(zhì)。陳裕華[15]及TIAN等[16]采用Huber-Markov隨機(jī)域進(jìn)行疊前三參數(shù)反演,避免邊緣的模糊化。GUO等[17]基于地震數(shù)據(jù)驅(qū)動及MRF的正則化方法進(jìn)行疊前地震反演,能夠以較低的不確定性揭示構(gòu)造細(xì)節(jié),并獲得同步反演結(jié)果。WANG等[18]通過在MRF的不同方向改變疊前反演的勢函數(shù)以體現(xiàn)各向異性,對頁巖氣甜點(diǎn)識別進(jìn)行研究。但現(xiàn)有文獻(xiàn)大多采用各向同性的MRF,往往忽略了MRF的方向性特點(diǎn),尤其是對于構(gòu)造較為復(fù)雜的介質(zhì)等。在復(fù)雜介質(zhì)中,水平、垂直和對角方向上模型梯度值的先驗(yàn)作用差異很大。為此,有必要詳細(xì)分析在地震反演中MRF的各向異性特征。

本文針對各向異性介質(zhì),通過各向異性參數(shù)模型的先驗(yàn)概率分布分析,闡述了地震反演中MRF的各向異性特性,并基于各向異性擴(kuò)散法進(jìn)行擴(kuò)散因子和輔助模型的優(yōu)選并應(yīng)用于疊前地震反演。最終,通過合成數(shù)據(jù)與實(shí)際資料測試說明了方法的合理與有效。

1 疊前多參數(shù)同步非線性地震反演

1.1 疊前地震響應(yīng)數(shù)學(xué)模型

疊前角度道集的地震響應(yīng)數(shù)學(xué)模型采用褶積模型進(jìn)行表征,即：

Y(θ)=R(θ,e)*W(θ)+N

(1)

式中：Y(·)為疊前地震數(shù)據(jù);R(·)為地震反射系數(shù)序列;W(·)為疊前角度子波;*表示褶積運(yùn)算;e為表征地層的彈性參數(shù);N一般假設(shè)為高斯白噪,其數(shù)學(xué)期望為零,協(xié)方差為σ2;θ為入射角。

基于平面波AVO理論,平面縱波入射至地層界面可產(chǎn)生反射縱波、反射橫波、透射縱波及透射橫波這4類不同類型的平面波,其反射與透射能量分布可以用反射、透射系數(shù)進(jìn)行表征。對疊前地震反射系數(shù)R(θ)進(jìn)行精確求解通常采用描述平面波反射和透射系數(shù)隨入射角變化的Zoeppritz方程,即：

其中,

式中：θ1,θ2分別為縱波入射(反射)角與縱波透射角;φ1,φ2分別為橫波反射角與橫波透射角;R,T分別為反射系數(shù)和透射系數(shù),其第1個角標(biāo)指示入射波類型,第2個角標(biāo)指示反射波或透射波類型(P為縱波,S為橫波);vP,vS,ρ分別為縱波速度、橫波速度與密度,其角標(biāo)數(shù)字為地層代號,分界面上層為1,下層為2。

1.2 反演目標(biāo)函數(shù)

地震反演問題是尋找最匹配觀測數(shù)據(jù)的一組模型參數(shù)的過程,該過程可以通過對目標(biāo)函數(shù)的優(yōu)化來實(shí)現(xiàn)。但在實(shí)際觀測數(shù)據(jù)不完整或有缺失的情況下,(1)式關(guān)于模型參數(shù)的求取是欠定的。因此,地震反演是一個不適定問題,需要引入正則化進(jìn)行約束以保證其適定性。針對(1)式,在反演中引入基于貝葉斯框架的正則化,即：

式中：P(m|Y)為已知觀測數(shù)據(jù)Y和待反演參數(shù)m的后驗(yàn)概率;P(m)為待反演參數(shù)的先驗(yàn)概率;P(Y|m)為似然函數(shù),表示模擬數(shù)據(jù)與觀測數(shù)據(jù)的誤差;P(Y)為觀測數(shù)據(jù)Y的后驗(yàn)概率,為一常數(shù)。

由于(1)式中假設(shè)N為高斯白噪,其數(shù)學(xué)期望為零,協(xié)方差為σ2,則似然函數(shù)P(Y|m)可表示為：

(5)

而先驗(yàn)概率P(m)則是對反演過程進(jìn)行約束的正則化項(xiàng),需要對模型地層邊界特征進(jìn)行保護(hù)以避免過度平滑,同時約束反演結(jié)果以保證穩(wěn)定。因此,引入MRF對二維模型參數(shù)的特征進(jìn)行描述,通過MRF構(gòu)造待反演參數(shù)的先驗(yàn)概率P(m)為：

反演問題可看作最大后驗(yàn)概率問題,對待反演參數(shù)m的最大后驗(yàn)概率估算可等價于相應(yīng)對數(shù)函數(shù)的估算,那么反演的目標(biāo)函數(shù)可寫為：

式中：J1(m)為數(shù)據(jù)項(xiàng);J2(m)為先驗(yàn)項(xiàng)或正則項(xiàng);λ為權(quán)系數(shù)。對于多參數(shù)(縱波速度、橫波速度及密度)同步反演,J2(m)由3部分組成：

J2(m)=J2P(vP)+J2S(vS)+J2d(ρ)=

式中：J2P(vP),J2S(vS)和J2d(ρ)分別為模型參數(shù)vP,vS和ρ的先驗(yàn)項(xiàng);δP,δS和δd分別為模型參數(shù)vP,vS和ρ的刻度參數(shù),它們的取值不同。而對于勢函數(shù)的選擇,由于平方型勢函數(shù)ψQ、對數(shù)型勢函數(shù)ψHL等對大梯度的接受概率較低,因此該處勢函數(shù)采用邊界保護(hù)特性較好的φGM,即：

式中：x為模型參數(shù)的梯度大小。

1.3 反演優(yōu)化算法

疊前多參數(shù)同步反演獲取匹配觀測數(shù)據(jù)的模型參數(shù),可以看作是對目標(biāo)函數(shù)((8)式)進(jìn)行極小化優(yōu)化的過程。但觀測數(shù)據(jù)與模型參數(shù)之間的關(guān)系非常復(fù)雜,若采用線性優(yōu)化算法,會存在依賴初始模型或反演陷入局部極值等問題。因此,為獲取全局最優(yōu)解,需要采用完全非線性優(yōu)化算法,如模擬退火算法、粒子群算法、遺傳算法等,本文采用粒子群算法(PSO)進(jìn)行目標(biāo)函數(shù)極小值的求解。

設(shè)在D維搜索空間中,由p個粒子組成一個粒子群M,其中每個粒子所在的位置代表目標(biāo)函數(shù)的一個潛在解,第i個粒子在n+1次迭代時,其位置更新為：

Mi(n+1)=Mi(n)+c3·vi(n+1)

(11)

其速度更新為：

vi(n+1)=vi(n)+c1·rand1(·)·

(12)

其中,

各參數(shù)的擾動可表示為：

m(s+1)=m(s)+Δm(s)

(14)

式中：m為待反演參數(shù);Δm表示該參數(shù)擾動量,由優(yōu)化算法求出;s表示該輪的迭代次數(shù)。

2 地震反演中MRF的各向異性特征

2.1 MRF的各向異性特征

圖1 空間點(diǎn)S的MRF點(diǎn)集

傳統(tǒng)的MRF對真實(shí)地下介質(zhì)情況的刻畫有一定偏差,未考慮實(shí)際地層分布與構(gòu)造情況。因此,為更好地闡述各向異性MRF對反演效果的影響機(jī)理,構(gòu)建3類處于地層不同部位的一階MRF的參數(shù)模型(圖2),假設(shè)中心待反演點(diǎn)實(shí)際vP為3000m/s,其橫向vP梯度值(50～100m/s)明顯小于縱向vP梯度值(150～250m/s),符合水平層狀介質(zhì)特點(diǎn)。此外,就3類模型的vP梯度而言,模型1位于地層內(nèi)部,周圍參數(shù)梯度值較小,可稱為中心對稱模型;模型2位于地層界面處,水平方向vP梯度變化小,而縱向vP梯度變化較大,可稱為縱向非對稱模型;模型3位于地層尖滅處,其四周vP梯度變化均較大,可稱為完全非對稱模型。

圖2 各向異性MRF模型及其在地層內(nèi)的分布

針對上述3種模型,由(6)式分析采用各向同性的MRF待反演參數(shù)的先驗(yàn)概率分布特征。圖3展示了各向同性MRF計(jì)算模型點(diǎn)vP的先驗(yàn)項(xiàng)概率分布,可以看出,vP梯度值中心對稱的模型1出現(xiàn)最大先驗(yàn)概率的vP為3000m/s,與中心點(diǎn)的實(shí)際vP一致;模型2和模型3出現(xiàn)最大先驗(yàn)概率的vP分別為3030m/s和3090m/s,偏離了中心點(diǎn)實(shí)際vP。這表明,對于模型2和模型3,采用各向同性的MRF獲得的中心點(diǎn)vP有一定的誤差,并且隨著各方向vP梯度差異性增大誤差也會加大。為此,在實(shí)際資料反演時,如果采用各向同性MRF,會給反演結(jié)果帶來很大的誤差。

圖3 采用各向同性MRF計(jì)算模型中心點(diǎn)vP的先驗(yàn)概率分布

2.2 各向異性MRF的建立

針對上述各向同性MRF對反演結(jié)果的偏差問題,提出采用各向異性MRF(即改變各個方向MRF的權(quán)系數(shù))對先驗(yàn)項(xiàng)的概率進(jìn)行校正,以降低各向異性對反演精度的影響,則(9)式先驗(yàn)項(xiàng)可以重寫為：

式中：α為方位校正系數(shù),該系數(shù)是針對MRF中的vP梯度值進(jìn)行調(diào)整的權(quán)系數(shù),通過調(diào)整不同方向之間的vP梯度值,保證地層的水平連續(xù)性和垂向分辨率,從而體現(xiàn)出MRF的各向異性。一般地,水平方向較大的α對vP梯度的懲罰大,僅能接受較小的梯度值,保持地層的水平連續(xù)性;反之,垂直方向上較小的α對vP梯度的懲罰小,可以允許較大的垂向梯度,避免反演結(jié)果過度平滑或不穩(wěn)定?？紤]α的上述特性,我們引入各向異性擴(kuò)散法,通過各方向上的vP梯度值對不同方向的系數(shù)α進(jìn)行自適應(yīng)求?。?/p>

圖4 采用各向異性MRF計(jì)算模型中心點(diǎn)vP的先驗(yàn)概率分布

圖5 不同梯度下方位校正系數(shù)隨擴(kuò)散因子的變化特征

3 合成數(shù)據(jù)測試

為了檢驗(yàn)各向異性MRF改善反演效果的有效性和精度,測試了理論模型的合成地震數(shù)據(jù)。理論模型為一個二維斷層背斜模型(圖6a,圖6b,圖6c),其合成數(shù)據(jù)模型共有61個CDP道集,每道150個采樣點(diǎn),采樣率為2ms,每個CDP有15個角度道集,角度范圍為3°～45°,角度間隔為3°。近、中、遠(yuǎn)3個角度子波峰值頻率分別為35,33,31Hz,均為零相位Ricker子波。通過縱波速度與橫波速度及密度之間的函數(shù)擬合關(guān)系構(gòu)建縱波速度初始模型(圖6d)。

圖6 合成數(shù)據(jù)模型a 縱波速度; b 橫波速度; c 密度; d 縱波速度初始模型

在模型測試中,刻度參數(shù)(δp,δs,δd),權(quán)系數(shù)λ=0.326,采用小角度(入射角為0～10°)道集,采用粒子群算法進(jìn)行優(yōu)化,設(shè)置迭代次數(shù)n=500,粒子數(shù)目i=20。為說明各向異性MRF對反演效果的改善,采用一階(k=1)的各向同性MRF與各向異性MRF進(jìn)行反演,其結(jié)果如圖7所示。特別地,計(jì)算方位校正系數(shù)α選擇真實(shí)縱波速度模型(圖6a)作為輔助模型I,擴(kuò)散因子q固定為300。從圖7可以看出,采用各向異性MRF進(jìn)行反演的結(jié)果在地層界面處的水平連續(xù)性及斷層處的垂向分辨率均較好,與各向同性MRF的結(jié)果相比,各向異性MRF的地層分布與斷層特征更明顯。

圖7 基于各向同性MRF(a)、基于各向異性MRF(b)的疊前同步反演結(jié)果

此外,不同的擴(kuò)散因子q對反演結(jié)果也存在影響。采用上述方法計(jì)算擴(kuò)散因子q分別為300和850、輔助模型I仍為實(shí)際縱波速度模型時的多個方向的方位校正系數(shù),并對其大小與變化情況進(jìn)行分析。由于一階MRF共有8個方向(水平、垂直、對角),其中3個方位的方位校正系數(shù)結(jié)果如圖8所示。

圖8 基于各向異性擴(kuò)散法獲得的方位校正系數(shù)a 水平(q=300); b 垂直(q=300); c 對角(q=300); d 水平(q=850); e 垂直(q=850); f 對角(q=850)

從圖8可以看出,當(dāng)q=300時,在模型的水平層狀區(qū)域,水平方向的方位校正系數(shù)(圖8a)大于垂直方向的方位校正系數(shù)(圖8b),以保證層狀地層的橫向連續(xù)性和垂向分辨率;在復(fù)雜的斷背斜區(qū)域,水平方向的方位校正系數(shù)小于垂直方向方位校正的系數(shù),主要是由于該區(qū)域橫向連續(xù)性變差,橫向差異性是主流特征(斷層兩盤的儲層參數(shù)的差異性)。4個對角上的方位校正系數(shù)(圖8c)也反映了上述特征。當(dāng)q=850時,方位校正系數(shù)的分布特征與q=300時基本一致。但隨著擴(kuò)散因子增大,邊界特征在各個方向的方位校正系數(shù)上體現(xiàn)更為明顯,即大的擴(kuò)散因子刻畫大梯度的邊界特征。

4 實(shí)例分析

采用基于各向異性MRF的疊前地震反演處理某二維地震剖面(圖9),該剖面地震數(shù)據(jù)共計(jì)2542個CMP道集,每個CMP共有15個角度道集(如圖10),入射角為3°～45°,角度間隔為3°,時窗范圍為2200～2600ms,采樣率為2ms,并且經(jīng)過2口井,井2位于第918道,井1位于第2226道。從圖9中可以看出,剖面的橫向連續(xù)性較差。根據(jù)測井?dāng)?shù)據(jù)的縱橫波速度比與密度交會結(jié)果制作氣層識別模板(圖11),該模板能夠明顯區(qū)分氣層與水層,但是很難有效劃分氣層與氣水同層,為此,在后續(xù)的解釋中將氣層與氣水同層合并進(jìn)行處理。通過測井?dāng)?shù)據(jù)進(jìn)行插值獲得主要的地層參數(shù),建立了反演的初始模型(圖12)。

圖9 二維地震剖面

圖10 疊前地震數(shù)據(jù)的角度道集

圖11 氣層識別模板(密度與縱橫波速度比交會結(jié)果)

圖12 反演初始模型a 縱波速度; b 橫波速度; c 密度; d 縱橫波速度比

首先,根據(jù)反演的縱波阻抗剖面刻畫了砂體的分布特征(圖13)。結(jié)合測井層位解釋,可將縱波阻抗值低于8750g/cm3·m/s的區(qū)域圈定為目標(biāo)層中砂體分布范圍,如圖13中藍(lán)色的陰影部分所示,大致可分為上、下兩層,兩層之間波阻抗較高的區(qū)域考慮為泥巖分布。

由于該區(qū)塊疊前地震資料品質(zhì)較差,導(dǎo)致進(jìn)行疊前同步反演后vS結(jié)果不理想,并且小角度(3°～9°)與大角度(33°～39°)疊前同步反演vP,vS,ρ三參數(shù)結(jié)果均較差,同時,縱橫波速度比γ=vP/vS是區(qū)分氣層和水層的重要參數(shù),需要其結(jié)果以指導(dǎo)氣層識別。因此,對于該區(qū)塊我們采用中等角度道集(18°～24°)疊前地震數(shù)據(jù)進(jìn)行vP,γ,ρ三參數(shù)同步反演而非傳統(tǒng)的vP,vS,ρ三參數(shù),從而避免vS結(jié)果的不穩(wěn)定影響縱橫波速度比γ對氣層的識別效果。

我們采用3種不同的測試結(jié)果來說明各向異性MRF對反演結(jié)果的橫向連續(xù)性及縱向分辨率的改善：①各向同性MRF;②各向異性MRF(擴(kuò)散因子q為715,輔助模型I采用vP的初始模型);③各向異性MRF(擴(kuò)散因子q為715,輔助模型I采用vP的疊后反演結(jié)果)。反演中優(yōu)化算法及參數(shù)設(shè)置與測試的合成數(shù)據(jù)參數(shù)一致(其中縱橫波速度比的刻度參數(shù)δγ采用δP/δS的值),參與反射系數(shù)計(jì)算的vS則由該次迭代的vP/γ求得,疊前直接反演縱橫波速度比的結(jié)果如圖14所示。由于縱橫波速度比對于氣層識別敏感性較好,此處采用縱橫波速度比反演結(jié)果作為指示氣藏分布情況,并根據(jù)氣藏分布情況說明3種不同的MRF對復(fù)雜地質(zhì)邊界恢復(fù)情況的改善。需要說明的是,根據(jù)圖11縱橫波速度比介于1.55～1.70為氣層/氣水同層分布范圍,圖14中深藍(lán)色部分縱橫波速度比均介于1.55～1.70,反演結(jié)果指示氣層/氣水同層分布區(qū)域。

圖14 18°～24°道集疊前地震數(shù)據(jù)直接反演縱橫波速度比結(jié)果a 各向同性MRF; b 各向異性MRF(擴(kuò)散因子q為715,以vP初始模型為輔助模型); c 各向異性MRF(擴(kuò)散因子為715,以vP疊后反演結(jié)果為輔助模型)

由實(shí)際地質(zhì)資料與二維地震剖面可知,該區(qū)塊的主要特征是地層的橫向連續(xù)性差。對比圖14a與圖14b 可以看出,采用各向異性MRF后,能夠較為清晰地刻畫該區(qū)塊地質(zhì)邊界不連續(xù)的特征,而各向同性的MRF反演結(jié)果則較為平滑,與實(shí)際介質(zhì)差別較大。對比圖14b與圖14c可以看出,改變輔助模型為vP疊后反演結(jié)果后,CDP1～2000處時窗上、下的縱橫波速度比值明顯偏低,此外,反演結(jié)果的縱向分層能力也明顯下降,說明反演結(jié)果中存在一些不穩(wěn)定成分,會給基于各向異性擴(kuò)散法計(jì)算的方位校正系數(shù)帶來不確定性。

為進(jìn)一步說明基于各向異性擴(kuò)散法計(jì)算的各向異性MRF對反演結(jié)果的影響,將CMP第2226道反演結(jié)果與井1測井曲線進(jìn)行對比,結(jié)果如圖15所示,可以看出,反演結(jié)果與測井曲線的吻合較好,整體來說,采用各向異性MRF及以vP初始模型為輔助模型的反演結(jié)果吻合程度比其它兩者高,但仍然存在一些差異,主要是由于未進(jìn)行測井校正以及地震數(shù)據(jù)與聲波測井?dāng)?shù)據(jù)之間未進(jìn)行匹配處理。

圖15 疊前地震數(shù)據(jù)直接反演獲得的縱橫波速度比結(jié)果與井1測井曲線的對比a 各向同性MRF; b 各向異性MRF(擴(kuò)散因子q為715,以vP初始模型為輔助模型); c 各向異性MRF(擴(kuò)散因子q為715,以vP疊后反演結(jié)果為輔助模型)

5 結(jié)論

本文提出了基于各向異性MRF的邊界保護(hù)正則化疊前地震反演方法。該方法通過在MRF的各個方向引入方位校正系數(shù),有效地消除層狀地層與構(gòu)造復(fù)雜地區(qū)各向異性梯度對反演結(jié)果產(chǎn)生的影響,提高了反演結(jié)果的精度,盡可能地恢復(fù)地質(zhì)模型的構(gòu)造細(xì)節(jié)。此外,建立了擴(kuò)散因子及輔助模型的優(yōu)選方法,提升了各向異性MRF方位校正系數(shù)的準(zhǔn)確性,進(jìn)一步改善了復(fù)雜介質(zhì)條件下反演結(jié)果的精度。

1) 表征地層不同位置的3種模型基于各向異性MRF,可以準(zhǔn)確地搜索到待反演值,而各向同性MRF待反演值與真實(shí)值偏差較大,說明介質(zhì)的各向異性特征對反演結(jié)果有顯著影響。

2) 基于各向異性擴(kuò)散法對方位校正系數(shù)進(jìn)行定量確定,不同方位的校正系數(shù)差異較大,并且隨著擴(kuò)散因子q增大,邊界特征在各個方位的校正系數(shù)上的體現(xiàn)更為明顯。此外,擴(kuò)散因子q與輔助模型I對方位校正系數(shù)均有影響,即擴(kuò)散因子q應(yīng)選擇適當(dāng)大的值,q越大對模型梯度變化越敏感。而對于輔助模型I的選取,應(yīng)盡可能接近實(shí)際縱波速度模型或以層速度、均方根速度等地震屬性替代,采用初始模型會使反演結(jié)果較為平滑,而采用反演后的結(jié)果會造成不穩(wěn)定。

3) 合成數(shù)據(jù)測試與實(shí)例分析的結(jié)果說明,該方法相比各向同性MRF,對于縱波速度、密度、縱橫波速度比的同步反演,其結(jié)果的橫向連續(xù)性和縱向分辨率有明顯改善。此外,反演結(jié)果與測井曲線的匹配程度較好,說明了方法的有效性。