牛麗萍,胡華鋒,周單,鄭曉東,耿建華

(1.中石化石油物探技術(shù)研究院有限公司,江蘇 南京 211103;2. 中國石油勘探開發(fā)研究院,北京 100083;3. 同濟大學(xué)海洋地質(zhì)國家重點實驗室,上海 200092;4.同濟大學(xué) 海洋與地球科學(xué)學(xué)院,上海 200092;5.同濟大學(xué) 海洋資源研究中心,上海 200092)

0 引言

油氣地震勘探的目的是探明地下介質(zhì)的幾何形態(tài)、巖性和物性甚至流體分布,進而為鉆井井位部署提供依據(jù)。目前,石油工業(yè)界廣泛采用疊前AVO反演獲得油氣儲層彈性參數(shù),在巖石物理分析的基礎(chǔ)上,對油氣儲層進行定量刻畫。疊前AVO反演的理論基礎(chǔ)是描述平面波在地下介質(zhì)分界面處發(fā)生反射和透射的Zoeppritz方程[1]。20世紀(jì),許多學(xué)者提出了大量的AVO線性近似公式[2-8],以此來解決Zoeppritz方程求解的非線性問題,并提高反演的計算效率。然而,這些線性AVO近似基于弱反差介質(zhì)與小入射角的假設(shè),在強反差介質(zhì)及大角度入射情形下會產(chǎn)生較大的計算誤差[9-11]。隨著油氣勘探的深入,勘探目標(biāo)轉(zhuǎn)向非均質(zhì)性強的復(fù)雜油氣藏,而長排列地震采集以及角度域成像技術(shù)的應(yīng)用也使得有效地震數(shù)據(jù)的入射角越來越大,因此迫切需要發(fā)展適應(yīng)強反差介質(zhì)以及大入射角情形下的AVO反演技術(shù)。

近些年,很多學(xué)者開始研究基于精確Zoeppritz方程的擬線性及非線性反演方法。劉福平等[12]首次給出了縱波反射系數(shù)對縱橫波速度一階偏導(dǎo)數(shù)的精確表達式;Zhi等[13]利用精確Zoeppritz方程和PP-PS聯(lián)合反演得到了界面兩側(cè)的密度比和3個速度比;Zong等[14]推導(dǎo)了Zoeppritz方程所描述的縱波反射系數(shù)的精確表達式,并結(jié)合非線性反演方法估計出了反射界面上下層介質(zhì)的6個彈性參數(shù);Zhi等[15]提出了基于Zoeppritz方程的迭代正則化Levenberg-Marquardt反演方法;Yin等[16]利用縱波速度反射率、橫波速度反射率和密度反射率重新表示Zoeppritz方程,然后基于逆算子估計實現(xiàn)三參數(shù)的非線性反演;張凌遠等[17]利用快速模擬退火算法實現(xiàn)了基于Zoeppritz方程的橫波速度、縱橫波速度比反演。盡管基于精確Zoeppritz方程的疊前反演方法獲得了較大的發(fā)展,但是這類反演方法的計算復(fù)雜度高,在實際應(yīng)用中,特別是針對低信噪比數(shù)據(jù),反演結(jié)果的準(zhǔn)確性以及穩(wěn)定性等依然存在較大挑戰(zhàn)[18-20]。

MCMC模擬是一種啟發(fā)式的全局尋優(yōu)算法,可以直接模擬任意復(fù)雜形式的概率分布。利用MCMC算法進行地震反演,不僅可以避免假設(shè)儲層彈性參數(shù)的分布形式服從高斯分布,而且也可以避免對彈性參數(shù)與地震數(shù)據(jù)之間的物理關(guān)系做線性假設(shè)。此外,基于MCMC模擬所獲得的后驗概率分布能夠完整地刻畫參數(shù)估計的不確定性,從而為更科學(xué)地利用反演結(jié)果對儲層進行刻畫提供了依據(jù)。因此,MCMC算法是一種在復(fù)雜儲層與大入射角地震數(shù)據(jù)條件下實現(xiàn)AVO彈性參數(shù)反演的理想選擇[21-22]。

盡管MCMC算法具有諸多優(yōu)點,但是在將其應(yīng)用于疊前彈性參數(shù)反演時仍然存在以下兩方面的問題:①對于實際復(fù)雜地下介質(zhì),彈性參數(shù)的真實概率分布往往十分復(fù)雜,很難尋找一個合適的概率密度函數(shù)來描述,常規(guī)利用高斯分布刻畫彈性參數(shù)的統(tǒng)計特征存在較大的局限性;②由于地下模型參數(shù)空間巨大以及地震數(shù)據(jù)中噪聲等因素的影響,MCMC對彈性參數(shù)后驗概率分布的搜索過程極易受到局部極值的影響,這使得MCMC采樣過程難以穩(wěn)定收斂,并且增加了反演結(jié)果的不確定性,特別是在低信噪比資料情況下,這種問題更加明顯。為了避免這些問題,大多數(shù)學(xué)者利用Zoeppritz方程的線性近似公式來描述地震正演過程[23-27],并未將MCMC算法應(yīng)用于求解基于精確Zoeppritz方程的疊前地震反演問題。Pan 等[28]提出了基于精確Zoeppritz方程和MCMC-DRAM(delayed rejection adaptive Metropolis)算法的疊前反演方法,但是DRAM算法無法有效地遍歷具有很多局部極值的復(fù)雜高維參數(shù)空間,而且也難以判斷馬爾科夫鏈何時達到收斂[29]。為了降低后驗概率分布的復(fù)雜度以提高MCMC采樣的穩(wěn)定性與收斂性,Pan等假設(shè)彈性參數(shù)相互獨立且服從高斯分布,然而,該假設(shè)條件在實際復(fù)雜地下儲層中是難以滿足的。

本文針對實際復(fù)雜儲層及低信噪比資料條件下基于精確Zoeppritz方程的疊前地震反演問題,從以下3個方面來提高MCMC疊前AVO反演的穩(wěn)定性與收斂性:①利用低頻模型約束,將待反演參數(shù)轉(zhuǎn)換為模型參數(shù)的擾動量,從而使得先驗概率分布由實際復(fù)雜概率分布轉(zhuǎn)換為高斯分布,以降低后驗概率分布的復(fù)雜度;②通過對似然函數(shù)取對數(shù),并利用低頻模型來約束地震正演過程;③利用基于隨機子空間采樣的多鏈模擬算法對疊前非線性反演問題進行全局尋優(yōu),以避免采樣過程過早地收斂到局部極值。低信噪比的模擬數(shù)據(jù)和實際數(shù)據(jù)測試表明,本文所提方法能夠有效提高反演結(jié)果的準(zhǔn)確性以及穩(wěn)定性。

1 方法原理

1.1 基于精確Zoeppritz方程的貝葉斯疊前地震反演

貝葉斯理論通過將模型參數(shù)視為隨機變量,把模型參數(shù)的先驗分布與似然函數(shù)結(jié)合起來,從而獲得同時反映模型合理性與數(shù)據(jù)擬合程度的后驗概率分布。當(dāng)給定地震觀測數(shù)據(jù)d時,模型參數(shù)x的后驗概率分布可以寫為[30]:

P(x|d)∝P(x)P(d|x),

(1)

式中:先驗分布P(x)表示獨立于實際地震觀測數(shù)據(jù)d的模型參數(shù)的已知信息;似然函數(shù)P(d|x)描述了實際地震觀測數(shù)據(jù)與由模型參數(shù)正演所得模擬地震數(shù)據(jù)之間的擬合程度。

對于實際地下介質(zhì)而言,彈性參數(shù)的先驗概率分布形式通常十分復(fù)雜,特別是對于復(fù)雜巖性儲層[31]。為了充分發(fā)揮先驗低頻模型在疊前隨機反演過程中的約束作用,我們將待反演參數(shù)由彈性參數(shù)本身轉(zhuǎn)換為彈性參數(shù)的擾動量。由于彈性參數(shù)的擾動量圍繞彈性參數(shù)的低頻趨勢m0而變化,因此,不管真實的彈性參數(shù)m服從簡單單峰分布還是復(fù)雜多峰分布,彈性參數(shù)擾動量Δm一般表現(xiàn)為單峰分布,其先驗概率分布可以用高斯分布來表示:

P(Δm)=Ν(Δm;μΔm,ΣΔm),

(2)

式中:μΔm和ΣΔm分別為擾動量的先驗期望和先驗協(xié)方差。根據(jù)貝葉斯理論,彈性參數(shù)擾動量的后驗概率分布可以表示為:

P(Δm|d)∝P(Δm)P(d|Δm),

(3)

假設(shè)地震數(shù)據(jù)中的噪聲服從高斯分布,同時考慮到數(shù)值計算的穩(wěn)定性,我們采用對數(shù)似然函數(shù):

(4)

dsyn=Sc(m)+e,

(5)

式中:S為地震子波矩陣;e為服從零均值高斯分布的噪聲;c(m)為利用精確Zoeppritz方程(見附錄A)計算得到的縱波反射系數(shù),且m=m0+Δm。

1.2 基于隨機子空間采樣的多鏈模擬

假設(shè)待反演地震道的時間采樣點數(shù)為T,則待反演的彈性參數(shù)擾動量個數(shù)d=3×T,d即為模型參數(shù)空間的維數(shù)。為了避免采樣過程過早地收斂到局部極值,本文通過同時運行多條初始狀態(tài)不同的馬爾科夫鏈來搜索模型空間。此外,在對高維模型參數(shù)空間的隨機搜索過程中,如果對所有模型參數(shù)同時進行更新,通常情況下生成的大量候選點會被拒絕,導(dǎo)致馬爾科夫鏈?zhǔn)諗糠浅＞徛?。為了提高候選點的接受率以加速馬爾科夫鏈的收斂過程,在每次迭代時,隨機選擇某些模型參數(shù)組成子空間,只對子空間內(nèi)的模型參數(shù)進行更新,而對子空間之外的模型參數(shù)保留原采樣點。

1)隨機生成彈性參數(shù)擾動量的初始樣本xi,結(jié)合該道的低頻模型得到相應(yīng)的3個彈性參數(shù)(縱波速度、橫波速度和密度),然后基于精確Zoeppritz方程和褶積模型正演得到該道的模擬地震記錄,根據(jù)式4計算該樣本所對應(yīng)的似然函數(shù),并將其與彈性參數(shù)擾動量的先驗分布(式2)相結(jié)合得到該樣本xi所對應(yīng)的后驗概率密度p(xi);

首先，農(nóng)村的“空心化”引發(fā)的人才短缺，在城鎮(zhèn)化、工業(yè)化不斷發(fā)展的當(dāng)下，大量農(nóng)村青壯年勞力不斷涌入城市，農(nóng)村常住人口多是一些“老弱婦幼”人員，并且隨著農(nóng)村人口逐漸空心化的同時，逐漸演變發(fā)展為農(nóng)村土地、產(chǎn)業(yè)以及各種基礎(chǔ)設(shè)施的整體“空心化”，當(dāng)下農(nóng)村嚴(yán)重缺乏能夠滿足現(xiàn)代服務(wù)業(yè)要求的高素質(zhì)勞動力人口、以及相應(yīng)基礎(chǔ)設(shè)施。

2)根據(jù)原采樣點xi計算第i條鏈的候選點zi:

zi=xi+dxi。

(6)

式中,dxi為模型參數(shù)的更新量,可以表示為:

(7)

式中:B為由隨機選擇的d′個模型參數(shù)組成的子空間;δ為馬爾科夫鏈的組數(shù);μ和ν均為從{1,…,i-1,i+1,…,N}中隨機取出δ個整數(shù)所組成的向量,且數(shù)值無重復(fù);C為尺度因子;ε為從高斯分布Nd′(0,c)中生成的隨機數(shù),通常可令c=10-6;

3)計算候選點所對應(yīng)的后驗概率密度p(zi),然后根據(jù)Metropolis準(zhǔn)則計算接受概率α(xi,zi)[32]:

(8)

4)判斷是否接受轉(zhuǎn)移。如果候選點被接受,則xi=zi;否則,xi=xi;

5)對每一個時刻k重復(fù)執(zhí)行步驟2)～4)直至完成N條馬爾科夫鏈的采樣,然后進入時刻k+1進行采樣直至達到最大迭代次數(shù)K。

圖1 對彈性參數(shù)擾動量的隨機采樣流程Fig.1 Flow chart of stochastic sampling for the elastic parameter perturbations

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

圖2 彈性參數(shù)的測井曲線及其低頻模型Fig.2 Real well logs of elastic parameters and the corresponding LFMs

圖3 彈性參數(shù)的交會及利用高斯混合分布擬合得到的概率密度等值線Fig.3 Crossplots of the elastic parameters and the corresponding probability density contours obtained by fitting with Gaussian mixture distributions

圖4 彈性參數(shù)擾動量的交會及利用高斯分布擬合得到的概率密度等值線Fig.4 Crossplots of elastic parameter perturbations and the corresponding probability density contours obtained by fitting with Gaussian distributions

利用精確Zoeppritz方程計算該測井?dāng)?shù)據(jù)的縱波反射系數(shù),并將其與主頻為45 Hz的Ricker子波相褶積得到無噪聲的地震角道集(圖5a),然后加入高斯隨機噪聲得到最終的模擬地震角道集(圖5b),模擬地震數(shù)據(jù)的信噪比為3,入射角范圍為3°～48°,角度采樣間隔為3°。圖6對比了本文所提改進后的反演方法與改進前的基于Zoeppritz方程反演的效果,可以看出,本文所提方法反演得到的彈性參數(shù)與實際測井曲線的擬合度更好。圖7為馬爾科夫鏈的收斂性指標(biāo)隨迭代次數(shù)的變化,在迭代約32 000次后,馬爾科夫鏈已經(jīng)達到收斂。圖8和圖9分別為1 630 ms和1 672 ms處3個彈性參數(shù)擾動量的后驗統(tǒng)計直方圖,不同位置處3個彈性參數(shù)擾動量均穩(wěn)定收斂到后驗高斯分布。

a—無噪聲;b—信噪比為3a—no noise;b—S/N is 3

圖6 含噪聲地震數(shù)據(jù)的彈性參數(shù)反演結(jié)果Fig.6 The inversion results of the elastic parameters with noisy seismic data

圖7 收斂性指標(biāo)隨迭代次數(shù)的變化Fig.7 The variation of convergence diagnostic with the number of iterations

圖8 時間為1 630 ms處彈性參數(shù)擾動量的后驗統(tǒng)計直方圖Fig.8 Statistical histogram of the sampling results for the elastic parameter perturbations at time 1630 ms

圖9 時間為1 672 ms處彈性參數(shù)擾動量的后驗統(tǒng)計直方圖Fig.9 Statistical histogram of the sampling results for the elastic parameter perturbations at time 1672 ms

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

實際二維地震數(shù)據(jù)的疊加剖面及測井位置如圖10所示,該數(shù)據(jù)的時間范圍為2 100～2 350 ms,采樣率是2 ms,地震子波主頻為20 Hz。圖11展示了不同CDP位置處的地震角道集,角度范圍為1°～31°,角度間隔為2°。從圖10和圖11可以看出,該剖面左側(cè)地震數(shù)據(jù)的信噪比較低,而在剖面的中間及右側(cè)地震數(shù)據(jù)的信噪比較高。

圖10 疊加地震剖面Fig.10 Poststack seismic section

圖11 不同CDP位置處的實際地震角道集Fig.11 Actual seismic angle gathers at different CDP

在地震層位的約束下,我們首先根據(jù)測井?dāng)?shù)據(jù)構(gòu)建彈性參數(shù)的低頻模型(圖12),然后通過對測井位置處彈性參數(shù)擾動量的統(tǒng)計分析得到其先驗期望和先驗協(xié)方差,圖13展示了彈性參數(shù)擾動量的交會及其高斯分布擬合結(jié)果,可以看出,高斯分布能夠很好地擬合彈性參數(shù)擾動量的統(tǒng)計分布特征。

a—縱波速度;b—橫波速度;c—密度a—P-velocity;b—S-velocity;c—density

圖14為利用本文所提方法反演得到的彈性參數(shù)剖面,圖15和圖16分別展示了W1井和W2井處的反演結(jié)果。從兩口井的測井曲線可以看出,該工區(qū)的彈性參數(shù)在2 200 ms附近存在一個非常明顯的突變界面,在該界面以下,兩口井的測井曲線存在明顯的差異,這反映了該工區(qū)地層橫向變化較大。二維剖面的反演結(jié)果與該先驗認(rèn)識一致。從井旁道的反演結(jié)果來看,本文所提方法的反演結(jié)果能夠較好地擬合實際測井曲線,而且基于反演結(jié)果的模擬地震角道集也與實際地震角道集相吻合。

a—縱波速度;b—橫波速度;c—密度a—P-velocity;b—S-velocity;c—density

a—縱波速度;b—橫波速度;c—密度;d—井旁地震角道集a—P-velocity;b—S-velocity;c—density;d—the seismic angle gather at well

a—縱波速度;b—橫波速度;c—密度;d—井旁地震角道集a—P-velocity;b—S-velocity;c—density;d—the seismic angle gather at well

a—縱波速度; b—橫波速度; c—密度

圖18 不同CDP位置處收斂性指標(biāo)隨迭代次數(shù)的變化Fig.18 The variation of convergence diagnostic with the number of iterations at different CDP

4 結(jié)論

本文針對實際復(fù)雜儲層及低信噪比資料條件下基于精確Zoeppritz方程的疊前地震反演問題,通過低頻模型約束先驗概率分布以及地震正演模擬過程,并結(jié)合對數(shù)似然函數(shù)和基于隨機子空間采樣的多鏈模擬算法,提出了一種改進的MCMC疊前地震反演方法。低信噪比的模擬數(shù)據(jù)測試表明,相較于改進前的反演方法,本文所提方法能夠獲得更加符合實際測井曲線的彈性參數(shù)反演結(jié)果。二維實際資料的應(yīng)用效果表明,盡管本文所提方法為逐道反演,但是剖面上相鄰兩道的反演結(jié)果之間具有良好的一致性,這很好地說明了本方法的穩(wěn)定性;此外,當(dāng)?shù)卣鹳Y料的信噪比越低時,彈性參數(shù)反演結(jié)果的不確定性越大,本文所提方法能夠給出可信的、定量的不確定性估計。

附錄A Zoeppritz方程

在各向同性的彈性介質(zhì)中,當(dāng)平面縱波入射到兩層半無限介質(zhì)的分界面處時,反射系數(shù)、透射系數(shù)與介質(zhì)彈性參數(shù)和入射角之間的關(guān)系可以表示為[1]:

式中:VP,VS和ρ分別表示縱波速度、橫波速度和密度,其下標(biāo)1和2分別表示上層和下層介質(zhì);θ1和φ1分別為P波和SV波的反射角;θ2和φ2分別為P波和SV波的透射角;RPP和RPS分別為P波和SV波的反射系數(shù);TPP和TPS分別為P波和SV波的透射系數(shù)。

附錄B 收斂性檢驗

對于N條采樣長度為T的馬爾科夫鏈,每條鏈內(nèi)的方差可以表示為:

其中,?·」表示取整運算;而所有鏈之間的方差可以表示為:

那么,馬爾科夫鏈的收斂性檢驗指標(biāo)可以表示為: