賀東陽，李海山，何 潤，王 偉

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

0 引言

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

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

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

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

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

1.1 混合高斯先驗分布

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

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

1.2 混合高斯線性反演解

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

1.3 線性混合高斯序貫采樣

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

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

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

式中：E表示單位矩陣。

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

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

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

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

線性混合高斯序貫采樣的步驟如下：

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

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

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

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

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

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

2 模型測試

為驗證本文提出的方法的有效性，這里利用二維模型來驗證本文提出方法的有效性。

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

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

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

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

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

3 實際應(yīng)用

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

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

4 結(jié)論

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

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

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