基于混合高斯先驗(yàn)分布的地質(zhì)統(tǒng)計(jì)學(xué)反演

賀東陽(yáng)，李海山，何 潤(rùn)，王 偉

（中國(guó)石油勘探開(kāi)發(fā)研究院西北分院，蘭州 730020）

0 引言

利用地表觀測(cè)的地球物理資料來(lái)推斷地下介質(zhì)物理參數(shù)是一個(gè)地球物理反演問(wèn)題。地震反演通過(guò)結(jié)合地震、測(cè)井、鉆井資料及已知的地質(zhì)規(guī)律來(lái)估計(jì)儲(chǔ)層彈性參數(shù)，對(duì)巖性預(yù)測(cè)和流體識(shí)別具有重要的指導(dǎo)意義。常規(guī)的確定性反演只能提供一個(gè)最優(yōu)解［1-3］，無(wú)法對(duì)模型參數(shù)進(jìn)行不確定性分析。而概率統(tǒng)計(jì)反演方法將反演結(jié)果以概率的形式表示，可以獲得模型參數(shù)的多次實(shí)現(xiàn)，能夠?qū)Ψ囱萁Y(jié)果進(jìn)行不確定性分析。因此，概率統(tǒng)計(jì)反演在一定程度上可以避免反演結(jié)果不確定性導(dǎo)致的預(yù)測(cè)風(fēng)險(xiǎn)。概率統(tǒng)計(jì)反演分為貝葉斯反演和地質(zhì)統(tǒng)計(jì)學(xué)反演兩類(lèi)方法。

在貝葉斯反演中，將模型參數(shù)的解表示為后驗(yàn)概論密度函數(shù)［4-6］。后驗(yàn)概率密度函數(shù)的解析解通常無(wú)法求取，可利用馬爾科夫蒙特卡洛（Markov chain Monte Carlo，McMC）方法進(jìn)行抽樣求解［7-8］。張廣智等［9］將該算法應(yīng)用于疊前地震反演，獲得了縱橫波阻抗和密度反射系數(shù)。王輝等［10］利用該算法在疊前域反演了縱橫波速度和密度。但對(duì)于高斯線(xiàn)性假設(shè)，即模型參數(shù)的先驗(yàn)概率密度函數(shù)是高斯分布并且正演算子可以線(xiàn)性化，則可以獲得后驗(yàn)概率密度函數(shù)的解析解［5］。Buland 等［11］通過(guò)線(xiàn)性高斯假設(shè)并利用疊前地震資料進(jìn)行縱橫波速度和密度參數(shù)反演。滕龍等［12］通過(guò)線(xiàn)性高斯假設(shè)進(jìn)行彈性參數(shù)與物性參數(shù)的聯(lián)合反演。為表征離散巖性和流體對(duì)模型參數(shù)的分布影響，Grana 等［13-14］將模型參數(shù)的先驗(yàn)分布表示為混合高斯分布，通過(guò)線(xiàn)性假設(shè)得到混合高斯線(xiàn)性反演解。

在地質(zhì)統(tǒng)計(jì)學(xué)反演中，通過(guò)序貫?zāi)M或非序貫?zāi)M算法來(lái)構(gòu)建模型參數(shù)的先驗(yàn)信息，然后不斷擾動(dòng)模型參數(shù)來(lái)匹配觀測(cè)數(shù)據(jù)，從而得到模型參數(shù)反演解，這種方法通常計(jì)算量非常大，且不能保證反演解收斂。Bortoli 等［15］和Haas 等［16］最早將序貫高斯模擬算法應(yīng)用于地質(zhì)統(tǒng)計(jì)學(xué)反演。余振等［17］利用序貫高斯模擬和序貫指示模擬算法進(jìn)行地質(zhì)統(tǒng)計(jì)學(xué)反演，并實(shí)現(xiàn)了高分辨率的流體預(yù)測(cè)。由于序貫高斯模擬通常須要對(duì)模型參數(shù)進(jìn)行正態(tài)變換，Soares 等［18］、Caetano 等［19］、Azevedo 等［20］利用直接序貫?zāi)M算法進(jìn)行地震隨機(jī)反演。為快速構(gòu)建模型參數(shù)的先驗(yàn)信息，孫瑞瑩等［21］、王保麗等［22］、Yang等［23］利用非序貫?zāi)M的FFT-MA 譜模擬算法進(jìn)行隨機(jī)反演。針對(duì)兩點(diǎn)地質(zhì)統(tǒng)計(jì)學(xué)的不足，González等［24］和劉興業(yè)等［25］利用多點(diǎn)地質(zhì)統(tǒng)計(jì)學(xué)進(jìn)行地震隨機(jī)反演。由于傳統(tǒng)的序貫?zāi)M算法通過(guò)求解克里金方程組來(lái)獲得待模擬點(diǎn)的均值和方差［26-29］，再利用蒙特卡洛的方法逐點(diǎn)隨機(jī)模擬該點(diǎn)的值。當(dāng)克里金方程組較大時(shí)，求解過(guò)程非常耗時(shí)。因此，Hansen 等［30］將線(xiàn)性高斯理論與序貫?zāi)M相結(jié)合來(lái)求取待模擬點(diǎn)的后驗(yàn)均值和方差，避免求解大型協(xié)方差矩陣和擾動(dòng)優(yōu)化過(guò)程，直接獲得模型參數(shù)的多次實(shí)現(xiàn)，但是該方法沒(méi)有考慮離散巖性和流體對(duì)模型參數(shù)的先驗(yàn)分布影響，只能單一地反演連續(xù)變量的模型參數(shù)。

在Hansen 等［30］的研究基礎(chǔ)上，將模型參數(shù)的先驗(yàn)分布表示為受離散巖性影響的混合高斯分布，并將線(xiàn)性混合高斯理論與序貫?zāi)M相結(jié)合，同時(shí)獲得聲波阻抗與巖性的多次實(shí)現(xiàn)，以期反演結(jié)果具有較好的空間連續(xù)性和穩(wěn)定性。

1 序貫?zāi)M的線(xiàn)性混合高斯反演

1.1 混合高斯先驗(yàn)分布

由于離散巖性的影響，模型參數(shù)如彈性阻抗、孔隙度、含水飽和度等連續(xù)變量呈現(xiàn)多峰的混合高斯分布。模型參數(shù)服從多峰的混合高斯概率密度函數(shù)，可以表示為

式中：L為高斯分量的個(gè)數(shù)，這里表示離散巖性的個(gè)數(shù)；pl為第l個(gè)高斯分量的先驗(yàn)權(quán)值；N為混合高斯的第l個(gè)高斯分量，表示為

1.2 混合高斯線(xiàn)性反演解

假設(shè)觀測(cè)數(shù)據(jù)d和模型參數(shù)m之間的正演關(guān)系可以用線(xiàn)性算子G表示，即d=G m+ε，地震噪聲ε服從高斯分布。如果模型參數(shù)m的先驗(yàn)信息服從混合高斯分布，那么模型參數(shù)m的后驗(yàn)條件分布也服從混合高斯分布［13-14］，如下：

1.3 線(xiàn)性混合高斯序貫采樣

為了提高反演解的精度和穩(wěn)定性，將低頻模型df作為待采樣點(diǎn)的條件數(shù)據(jù)，低頻模型可以通過(guò)井插值或者其他方法獲取。低頻模型可以表示為［31］：

式中：εf是一個(gè)協(xié)方差為∑f的噪聲，描述了反演結(jié)果與模型參數(shù)的相似度。

將d=G m+ε和式（7）聯(lián)合可以寫(xiě)成：

式中：E表示單位矩陣。

如果d′=，則式（8）可以寫(xiě)成：

根據(jù)Hansen 等［30］提出的線(xiàn)性高斯序貫采樣理論，在線(xiàn)性混合高斯序貫采樣過(guò)程中，須要估算待采樣點(diǎn)mx在三類(lèi)條件數(shù)據(jù)（mk，d及df）下的后驗(yàn)條件分布，這里的mk表示已采樣點(diǎn)的模型參數(shù)自身。該后驗(yàn)條件分布mx|(mkd,df), 也服從混合高斯分布，即mx|(mk,d′) 服從混合高斯分布：

其均值和協(xié)方差分別變?yōu)椋?/p>

后驗(yàn)條件權(quán)值變?yōu)椋?/p>

線(xiàn)性混合高斯序貫采樣的步驟如下：

（1）定義模擬網(wǎng)格，利用隨機(jī)種子隨機(jī)訪問(wèn)模擬網(wǎng)格中的待采樣點(diǎn)mx。

（2）利用式（11）—（13）分別計(jì)算待采樣點(diǎn)的后驗(yàn)條件均值、協(xié)方差以及后驗(yàn)條件權(quán)值

（3）比較后驗(yàn)條件權(quán)值系數(shù)確定待采樣點(diǎn)的離散巖性，即高斯分量，然后從確定的高斯分量中隨機(jī)抽取一個(gè)值，作為當(dāng)前待采樣點(diǎn)的連續(xù)變量的模型參數(shù)值。

（4）隨機(jī)訪問(wèn)下一個(gè)待采樣點(diǎn)，并將上一次采樣點(diǎn)的值作為該點(diǎn)的條件數(shù)據(jù)。

（5）重復(fù)步驟（1）—（4），直到訪問(wèn)完所有采樣點(diǎn)。

這種逐點(diǎn)序貫采樣的方法，在反演每一個(gè)點(diǎn)的時(shí)候，由于考慮了模型參數(shù)、觀測(cè)數(shù)據(jù)及低頻模型等多類(lèi)條件數(shù)據(jù)，實(shí)際上是對(duì)空間相關(guān)的后驗(yàn)分布進(jìn)行采樣，使得反演結(jié)果具有較好的空間連續(xù)性。通過(guò)選取合理的反演時(shí)窗和模擬鄰域，可以避免大型協(xié)方差矩陣的計(jì)算，從而提高計(jì)算效率。

2 模型測(cè)試

為驗(yàn)證本文提出的方法的有效性，這里利用二維模型來(lái)驗(yàn)證本文提出方法的有效性。

圖1（a）是聲波阻抗模型，圖1（b）是巖相模型，圖1（c）是真實(shí)記錄，一共206 道，每道302 個(gè)采樣點(diǎn)，并假設(shè)每個(gè)網(wǎng)格大小為1 m×1 m。模型中的黑線(xiàn)表示抽取的偽井，偽井可以用來(lái)求取協(xié)方差矩陣，并在反演過(guò)程中充當(dāng)條件數(shù)據(jù)。圖2 是從抽取的偽井中統(tǒng)計(jì)的變差函數(shù)，利用球狀理論變差函數(shù)來(lái)擬合實(shí)驗(yàn)變差函數(shù)，聲波阻抗在橫向上的變程為110 m，縱向上的變程為75 m。利用變差函數(shù)與協(xié)方差函數(shù)的關(guān)系：γ（h）=1-C（h），可以求取聲波阻抗的空間協(xié)方差。首先進(jìn)行無(wú)噪測(cè)試，反演結(jié)果如圖3 所示?？梢钥吹剑囱莸膬纱温暡ㄗ杩?、巖性及合成記錄與模型數(shù)據(jù)吻合較好。其中2 次巖性反演結(jié)果的歸類(lèi)正確率分別為93.63%和93.53%。圖3（g）和圖3（h）是2 次反演的砂巖概率，可以用于巖性的不確定性分析，然后進(jìn)行信噪比為4 的加噪測(cè)試，反演結(jié)果如圖4 所示?？梢钥吹?，2 次的反演結(jié)果與合成記錄的精度較圖3 有所下降，但仍然與模型數(shù)據(jù)能夠較好地吻合。其中，2 次巖性反演結(jié)果的歸類(lèi)正確率分別為92.18%和92.24%。圖4中砂巖概率的2 次反演結(jié)果較圖3 的結(jié)果出現(xiàn)了一定的波動(dòng)性，這是加入的噪聲所引起的，但這并沒(méi)有對(duì)巖性的歸類(lèi)結(jié)果產(chǎn)生大的影響。

圖1 模型數(shù)據(jù)Fig.1 Model data

圖2 聲波阻抗的變差函數(shù)統(tǒng)計(jì)Fig.2 Variogram statistics of acoustic impedance

圖3 無(wú)噪聲時(shí)的兩次反演結(jié)果Fig.3 Two inversion results without noise

圖4 信噪比為4 時(shí)的兩次反演結(jié)果Fig.4 Two inversion results when SNR is 4

3 實(shí)際應(yīng)用

為了驗(yàn)證該方法的實(shí)際效果，以國(guó)內(nèi)東部某油田的工區(qū)資料進(jìn)行了應(yīng)用。該工區(qū)以砂泥巖為主，且砂巖的阻抗較高，泥巖的阻抗較低。圖5（a）為原始的二維疊后地震記錄，共有139 道，反演時(shí)窗為2 450～2 680 ms，采樣率為2 ms，其中A 井用作驗(yàn)證井。圖5（b）和圖5（c）分別是聲波阻抗和巖性的反演結(jié)果，可看到反演結(jié)果具有較好的橫向展布，將反演的聲波阻抗與子波合成的地震記錄如圖5（d）所示，與真實(shí)的地震記錄吻合較好。圖5（e）是砂巖概率的反演結(jié)果，可以預(yù)測(cè)有利儲(chǔ)層的位置。為了驗(yàn)證本文提出的方法的可靠性，將A 井位置處抽取的反演結(jié)果和A 井?dāng)?shù)據(jù)對(duì)比，如圖6 所示，可以看到聲波阻抗吻合較好，巖性結(jié)果除了少數(shù)點(diǎn)被錯(cuò)誤歸類(lèi)外，大部分與真實(shí)巖性吻合較好，歸類(lèi)正確率達(dá)82.76%。

圖5 地震記錄和反演結(jié)果Fig.5 Real seismogram and inversion results

4 結(jié)論

（1）混合高斯模型能夠表征模型參數(shù)的多峰分布，將其與地質(zhì)統(tǒng)計(jì)學(xué)的序貫?zāi)M相結(jié)合可以同時(shí)反演連續(xù)變量和離散巖性。相比傳統(tǒng)的地質(zhì)統(tǒng)計(jì)學(xué)反演，這種方法在序貫采樣的過(guò)程中，同時(shí)考慮了模型參數(shù)自身、觀測(cè)數(shù)據(jù)及低頻模型，屬于邊模擬邊反演的過(guò)程，而不是先模擬再反演的過(guò)程。

（2）低頻模型的加入能讓反演結(jié)果更加穩(wěn)定，不易受噪聲的影響，空間相關(guān)的后驗(yàn)序貫采樣使得反演結(jié)果具有較好地空間連續(xù)性。這種方法還可以用于彈性阻抗反演、疊前AVO 反演。

（3）需要指出的是，這種邊模擬邊反演的方法由于直接以地震數(shù)據(jù)作為條件數(shù)據(jù)，缺少了模型參數(shù)的小尺度變化性，使得這種方法的垂向分辨率受到一定的影響。