葉端南，印興耀，孫瑞瑩，王保麗

（中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，青島 266580）

基于地質(zhì)統(tǒng)計(jì)先驗(yàn)信息的隨機(jī)地震反演

葉端南，印興耀，孫瑞瑩，王保麗

（中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，青島 266580）

基于地質(zhì)統(tǒng)計(jì)先驗(yàn)信息的隨機(jī)地震反演方法是一種基于蒙特卡洛的非線性反演方法。在貝葉斯理論框架下，通過序貫高斯模擬方法（sequential Gaussian simulation，SGS）和逐漸變形算法（Gradual Deformation Method，GDM）得到基于地質(zhì)統(tǒng)計(jì)學(xué)的先驗(yàn)信息，然后構(gòu)建似然函數(shù)，最終利用Metropolis算法實(shí)現(xiàn)后驗(yàn)概率密度的抽樣，得到反演問題的解。與確定性反演結(jié)果相比，該方法能夠有效地融合測井資料中的高頻信息，提高反演結(jié)果的分辨率。數(shù)值模擬試驗(yàn)表明：本方法的反演結(jié)果與理論模型吻合較好，具有較高的分辨率；序貫高斯模擬采用一種新的逐點(diǎn)模擬方式，并結(jié)合GDM，有效提高了隨機(jī)反演的計(jì)算效率。

地質(zhì)統(tǒng)計(jì)先驗(yàn)信息；貝葉斯理論；高分辨率；計(jì)算效率

0 引言

石油勘探和開發(fā)已經(jīng)從簡單的構(gòu)造識(shí)別轉(zhuǎn)向復(fù)雜構(gòu)造、薄儲(chǔ)層以及老油區(qū)剩余油的研究。目前工業(yè)目標(biāo)地質(zhì)體的厚度一般在10m甚至10m以下，而實(shí)際的反演方法較難區(qū)分10m以下的薄層頂?shù)祝?］，所以提高分辨率以反演薄層信息在儲(chǔ)層識(shí)別和油藏描述中起到了重要的作用。隨著地球物理研究的深入，非線性反演方法更適用于復(fù)雜隱蔽性儲(chǔ)層的研究。地球物理反演中普遍存在非線性和不適定性的問題，若采用線性方法來反演非線性反演問題通常較難得到其真值［2］，而隨機(jī)反演能夠在空間相關(guān)性和井的約束下，模擬出地震頻帶以外的信息，分辨率高于常規(guī)的確定性反演方法。因此，基于統(tǒng)計(jì)反演理論的非線性隨機(jī)地震反演方法得到了進(jìn)一步的發(fā)展［］。

Hass等［4］最早將地質(zhì)統(tǒng)計(jì)學(xué)中的序貫高斯模擬技術(shù)與地震反演方法相結(jié)合，并以與實(shí)際地震數(shù)據(jù)的相似程度作為判別準(zhǔn)則外推整個(gè)波阻抗數(shù)據(jù)體，這就是隨機(jī)地震反演的雛形；Dubrule等人［5－6］進(jìn)一步研究和發(fā)展了隨機(jī)地震反演方法。序貫高斯模擬為經(jīng)典的條件模擬方法，它是基于蒙特卡洛模擬思想發(fā)展起來的［7－8］，但是序貫高斯模擬方法的計(jì)算速度較慢。這里考慮采用Zou等人［9］提出的新的序貫高斯模擬實(shí)現(xiàn)方式，并結(jié)合Hu［10］提出的逐步變形更新算法（GDM）構(gòu)建地質(zhì)統(tǒng)計(jì)先驗(yàn)信息，以提高計(jì)算效率，從而實(shí)現(xiàn)隨機(jī)地震反演。Kj?nsberg等人［11］建議采用貝葉斯理論進(jìn)行反演以獲取后驗(yàn)概率密度的樣本，這里假定基于地質(zhì)統(tǒng)計(jì)學(xué)的先驗(yàn)信息和似然函數(shù)都服從高斯分布，那么得到的后驗(yàn)信息就服從非高斯分布，即不能得到后驗(yàn)概率密度的解析解。對(duì)于后驗(yàn)概率密度無法用公式表達(dá)的情況，需要對(duì)后驗(yàn)概率密度進(jìn)行抽樣來求解反問題。作者采用的抽樣方法是Metropolis算法，基于蒙特卡洛的地震反演方法結(jié)合了Metropolis抽樣算法、SGS算法和GDM擾動(dòng)更新算法，可以為先驗(yàn)信息提供多參數(shù)擾動(dòng)，從而加快收斂速度，提高反演精度。

1 方法原理

利用貝葉斯理論尋找反演問題的解，貝葉斯理論的一個(gè)優(yōu)勢就是在反演彈性參數(shù)的同時(shí)，可以得到不確定性的估計(jì)?；舅悸肪褪窃谪惾~斯理論框架下，通過SGS算法得到高斯概率密度，再利用GDM擾動(dòng)更新先驗(yàn)信息，然后根據(jù)地質(zhì)統(tǒng)計(jì)先驗(yàn)信息和原始地震資料確定似然函數(shù)，最后利用Metropolis算法對(duì)后驗(yàn)概率密度進(jìn)行抽樣得到反演問題的解。圖1為基于地質(zhì)統(tǒng)計(jì)先驗(yàn)信息的隨機(jī)地震反演流程圖。

圖1 基于地質(zhì)統(tǒng)計(jì)先驗(yàn)信息的隨機(jī)地震反演流程圖Fig．1 Flowchart of stochastic seismic inversion based on the geostatistical priori information

1．1 貝葉斯理論

AVO反演的理論基礎(chǔ)是Zoeppritz方程及其近似式。這里提出的隨機(jī)地震反演方法在計(jì)算反射系數(shù)時(shí)，直接采用精確的Zoeppritz方程計(jì)算縱波反射系數(shù)，表示形式為：

其中：Rpp表示P波反射系數(shù)；Rps表示S波反射系數(shù)；Tpp表示P波透射系數(shù)；Tps表示S波透射系數(shù)；θ1和θ2分別為P波反射角和透射角；φ1和φ2分別為S波反射角和透射角；Δα、Δβ和Δρ分別表示界面兩側(cè)的縱、橫波速度差及密度差，Δα＝（α2－α1），Δβ＝（β2－β1），Δρ＝（ρ2－ρ1）；α、β和ρ分別表示縱波速度、橫波速度及密度；下標(biāo)1和2分別表示界面上、下兩側(cè)的參數(shù)信息。

假設(shè)地震褶積模型為：

其中：d為觀測地震數(shù)據(jù)；r為反射系數(shù)序列；G為子波褶積矩陣；e為噪音。根據(jù)貝葉斯理論可知［12］：

其中：ρM（R）表示儲(chǔ)層彈性參數(shù)的先驗(yàn)信息；L（d／R）代表似然函數(shù)，表示模型與數(shù)據(jù)的匹配程度；σM（R）就是貝葉斯后驗(yàn)分布；k是概率歸一化常數(shù)。這里直接采用精確的Knott-Zoeppritz方程計(jì)算正演合成地震記錄，在一定程度上，可以減少Zoeppritz方程近似引起的誤差。

1．2 地質(zhì)統(tǒng)計(jì)先驗(yàn)信息

序貫高斯模擬（SGS）是常用的先驗(yàn)信息構(gòu)建方法，但是SGS存在計(jì)算量較大的缺點(diǎn)，所以采用一種新的逐點(diǎn)模擬方式對(duì)測井?dāng)?shù)據(jù)進(jìn)行序貫高斯模擬，為獲得更多的模擬實(shí)現(xiàn)，結(jié)合GDM更新算法控制先驗(yàn)概率密度的微擾，而不是利用SGS重復(fù)模擬得到多個(gè)實(shí)現(xiàn)，這在一定程度上進(jìn)一步提高了計(jì)算效率。因此可以結(jié)合SGS和GDM獲得有效的先驗(yàn)信息，這就是所謂的基于地質(zhì)統(tǒng)計(jì)學(xué)的先驗(yàn)信息。

1．2．1 SGS算法

序貫?zāi)M的基本思想是沿隨機(jī)模擬路徑，將某一鄰域內(nèi)的所有已知數(shù)據(jù)（原始數(shù)據(jù)和先前已模擬的數(shù)據(jù)）作為條件數(shù)據(jù)進(jìn)行模擬，獲取每個(gè)網(wǎng)格節(jié)點(diǎn)的條件累積概率密度（CCDF），進(jìn)而從CCDF得到模擬值。序貫高斯模擬適用于高斯模型，因此需確保原始數(shù)據(jù)服從高斯分布，若不服從，必須進(jìn)行高斯變換，模擬結(jié)束后，進(jìn)行高斯反變換。

以前的序貫高斯模擬是采用每次模擬整道數(shù)據(jù)或者逐點(diǎn)多次模擬的方式來匹配地震數(shù)據(jù)，這里使用Zou等人［9］提出的新的序貫高斯模擬實(shí)現(xiàn)方式，每次僅模擬一個(gè)點(diǎn)，并且每個(gè)網(wǎng)格只模擬一次，直至模擬完所有的網(wǎng)格節(jié)點(diǎn)，得到新的模型，該模擬方式在一定程度上可以提高反演結(jié)果的計(jì)算效率。

1．2．2 GDM更新算法

傳統(tǒng)的隨機(jī)反演方法在進(jìn)行模擬時(shí)，即使達(dá)到了Metropolis抽樣的接受條件，也不能保證模擬結(jié)果一定會(huì)使目標(biāo)函數(shù)收斂，這里在利用新的SGS模擬方式提高計(jì)算效率的同時(shí)，引入了逐步變形更新算法，擾動(dòng)更新SGS的實(shí)現(xiàn)，以保證模擬搜索時(shí)目標(biāo)函數(shù)快速地收斂，從而達(dá)到提高反演精度的目的。

逐步變形算法最早由Hu等［13］提出，用來逐步修改高斯分布的儲(chǔ)層模型，其后被擴(kuò)展到非高斯分布儲(chǔ)層的序貫指示模擬［14］。假設(shè)存在兩個(gè)獨(dú)立的高斯隨機(jī)函數(shù)Z1和Z2，GDM算法可用于得到一個(gè)新的高斯場，即表示為兩個(gè)獨(dú)立的高斯隨機(jī)函數(shù)的線性組合：

其中：p的取值范圍為0～1／2；Z、Z1和Z2分別為更新的高斯白噪、當(dāng)前待更新的高斯白噪和新加入的高斯白噪。由于高斯白噪Z的存在使得SGS算法具有隨機(jī)性，不管模擬網(wǎng)格的大小，只需要對(duì)變形參數(shù)進(jìn)行擾動(dòng)就可以對(duì)整個(gè)模型進(jìn)行修改，但由于沒有改變協(xié)方差結(jié)構(gòu)，所以這種擾動(dòng)不會(huì)影響數(shù)據(jù)的空間變差結(jié)構(gòu)。采用新的SGS模擬方式和GDM更新算法構(gòu)建地質(zhì)統(tǒng)計(jì)先驗(yàn)信息，可以有效地提高計(jì)算效率。

1．3 似然函數(shù)構(gòu)建

對(duì)后驗(yàn)信息進(jìn)行抽樣時(shí)，接受概率的建立需要使用似然函數(shù)，似然函數(shù)采用如下形式：

其中：si代表模擬地震波形的振幅，是由彈性參數(shù)計(jì)算得到的合成記錄；siobs是觀測數(shù)據(jù)的采樣點(diǎn)；σ是期望數(shù)據(jù)不確定性的標(biāo)準(zhǔn)差；Ri代表反演的彈性參數(shù)；Ri0是確定性反演中井位置處的低頻約束（平滑約束和點(diǎn)約束）［15］；α是可調(diào)參數(shù)，當(dāng)?shù)卣鹳Y料信噪比較高時(shí)，可選用較小的α值，當(dāng)?shù)卣鹳Y料中含有多種噪音時(shí)，適當(dāng)?shù)卦龃螃林怠?/p>

1．4 Metropolis抽樣

Metropolis抽樣是一種基于蒙特卡洛的抽樣方法［16］。需要對(duì)后驗(yàn)概率密度進(jìn)行Metropolis抽樣來得到反演問題的解。Metropolis準(zhǔn)則表示為：

其中：mpropose是先驗(yàn)概率密度的信息；maccept是先前接受模型的擾動(dòng)信息；Paccept是接受概率；其中的L就是似然函數(shù)。Gelman等人［17］發(fā)現(xiàn)，對(duì)高維分布來說接受概率應(yīng)該約為23%，對(duì)于較大的接受概率，算法探索后驗(yàn)概率密度的速度過于緩慢。另一方面，對(duì)于較小的接受概率，計(jì)算試驗(yàn)量較大，所以應(yīng)該調(diào)整各參數(shù)，例如變差函數(shù)的各種參數(shù)［18］、GDM擾動(dòng)區(qū)域和GDM算法中p的取值等，使接受概率保持在23%左右。

在貝葉斯理論框架下，通過新的SGS方式得到高斯概率密度，然后利用GDM算法擾動(dòng)更新模擬實(shí)現(xiàn)，就可以得到基于地質(zhì)統(tǒng)計(jì)學(xué)的先驗(yàn)信息。根據(jù)確定性反演中井位置處的低頻約束構(gòu)建似然函數(shù)，最后利用Metropolis算法對(duì)后驗(yàn)概率密度進(jìn)行抽樣得到了反演問題的解。

2 模型試算

首先對(duì)如圖2所示的一維數(shù)據(jù)進(jìn)行試算分析，從反演結(jié)果可以看出，無論是反射系數(shù)、波阻抗還是合成地震記錄，都和真實(shí)值比較匹配。

然后采用作者的反演方法對(duì)圖3（a）所示的二維模型進(jìn)行反演，該模型存在三個(gè)薄層，從中抽取10道數(shù)據(jù)作為偽井，如圖3（b）所示，采用的子波是主頻為30Hz、2ms采樣的雷克子波，將模型的合成地震記錄作為實(shí)際地震記錄進(jìn)行反演。圖4為二維反演結(jié)果，其中（a）為隨機(jī)道反演（反演的地震道數(shù)是隨機(jī)的）的結(jié)果，（b）為逐道反演（按順序，逐道進(jìn)行反演）的結(jié)果，可以看出，利用這兩種反演方式所得的結(jié)果和真實(shí)模型比較相似（似然函數(shù)中加入了較多的低頻約束信息），并且能較好地識(shí)別薄層。利用基于地質(zhì)統(tǒng)計(jì)先驗(yàn)信息的隨機(jī)地震反演方法能夠有效融合測井資料的高頻信息，提高反演結(jié)果的分辨率，識(shí)別較薄的儲(chǔ)層。

圖2 單道反演結(jié)果對(duì)比Fig．2 Comparison of inversion results and real data

圖3 原始數(shù)據(jù)Fig．3 The original data

3 實(shí)際資料測試

以國內(nèi)某油田某區(qū)塊的實(shí)際資料驗(yàn)證該方法的應(yīng)用效果。圖5是疊后地震數(shù)據(jù)，圖6是稀疏脈沖反演與隨機(jī)反演的對(duì)比，其中圖6（a）是稀疏脈沖反演結(jié)果，圖6（b）是隨機(jī)地震反演結(jié)果，黃色為含油砂巖（對(duì)應(yīng)較高阻抗），藍(lán)色為泥巖和含水砂巖（對(duì)應(yīng)較低阻抗）。從圖6中可以明顯地看出，隨機(jī)反演結(jié)果的分辨率高于稀疏脈沖反演。稀疏脈沖反演最多可識(shí)別7m左右的儲(chǔ)層，僅比地震數(shù)據(jù)的分辨率高一些，隨機(jī)反演可以識(shí)別3m左右的儲(chǔ)層，反演精度比較高。如圖6中黑色橢圓中所示，井A上有一個(gè)6．6 m的油層，是一個(gè)薄互層（由兩個(gè)薄油層構(gòu)成），利用稀疏脈沖反演方法不能分辨，而隨機(jī)反演可以將兩者區(qū)分開來。井B黑色圓圈內(nèi)上方有一個(gè)3m左右的薄油層，稀疏脈沖反演沒有反演出來，而隨機(jī)反演卻分辨出了該薄層。因此通過實(shí)際資料分析可以得到，隨機(jī)反演的分辨率高于稀疏脈沖反演，這里并不是說隨機(jī)反演方法優(yōu)于確定性反演，而是兩者各有其適用性，需要根據(jù)具體問題具體分析，因?yàn)榇_定性

反演速度較快，當(dāng)反演較厚的儲(chǔ)層時(shí)，可以應(yīng)用稀疏脈沖反演；而當(dāng)實(shí)際工區(qū)需要進(jìn)行薄層的識(shí)別時(shí)需要應(yīng)用隨機(jī)地震反演方法。

圖4 隨機(jī)地震反演結(jié)果Fig．4 Stochastic seismic inversion results

圖5 疊后地震剖面Fig．5 The post－stack seismic profile

4 結(jié)論與認(rèn)識(shí)

序貫高斯模擬是經(jīng)典的地質(zhì)統(tǒng)計(jì)學(xué)模擬方法，而GDM更新算法能夠連續(xù)更新儲(chǔ)層參數(shù)實(shí)現(xiàn)使目標(biāo)函數(shù)快速收斂，并且隨機(jī)地震反演能夠有效融合測井的高頻信息，提高反演結(jié)果的分辨率，但計(jì)算效率是其面臨的一個(gè)問題。結(jié)合序貫高斯模擬和GDM更新算法進(jìn)行隨機(jī)反演，新的模擬方式結(jié)合GDM有效提高了計(jì)算效率，并且通過模型試算和實(shí)際資料分析結(jié)果表明，該反演方法能夠穩(wěn)定地收斂，反演結(jié)果比較可靠。

需要強(qiáng)調(diào)的是，基于地質(zhì)統(tǒng)計(jì)先驗(yàn)信息的隨機(jī)地震反演方法中，變差函數(shù)的影響非常重要，特別是二維和三維的反演過程中，各向異性變差函數(shù)對(duì)反演結(jié)果的影響會(huì)很明顯，可考慮通過各向異性變差函數(shù)改善隨機(jī)反演的橫向連續(xù)性。

圖6 稀疏脈沖反演與隨機(jī)反演的對(duì)比Fig．6 Comparison of CSSI and SI

［1］CHOPRA S，CASTAGNA J，PORTNIAGUINE O．Seismic resolution and thin－bed reflectivityinversion［J］．CSEG recorder，2006，31（1）：19－25．

［2］ASTER R C，BORCHERS B，THURBER C H．Parameter estimation and inverse problems［M］．Academic Press，2013．

［3］FRANCIS A．Limitations of deterministic and advantages of stochastic seismic inversion［J］．CSEG Recorder，2005，30：5－11．

［4］HAAS A，DUBRULE O．Geostatical inversion－a sequential method of stochastic reservoir modelling constrined by seismic data［J］．First break，1994，13（12）：561－569．

［5］DUBRULE O，THIBAUT M，LAMY P，et al．Geostatistical reservoir characterization constrained by 3Dseismic data［J］．Petroleum Geoscience，1998，4（2）：121－128．

［6］ROTHMAN D H．Geostatistical inversion of 3D seismic data for thin－sand delineation［J］．Geophysics，1998，51（2）：332－346．

［7］印興耀，賀維勝，黃旭日．貝葉斯－序貫高斯模擬方法［J］．石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2006，29（5）：28 －32．

YIN X Y，HE W S，HUANG X R．Bayesian sequential Gaussian simulation method［J］．Journal of the University of Petroleum，2006，29（5）：28－32．（In Chinese）

［8］印興耀，劉永社．儲(chǔ)層建模申地質(zhì)統(tǒng)計(jì)學(xué)整地震數(shù)據(jù)的方法及砰究進(jìn)展［J］．石油地球物理勘探，2002，37（4）：423－430．

YIN X Y，LIU Y S．Methods and development of integrating seismic data in reservoir model－building［J］．OGP，2002，37（4）：423－430．（In Chinese）

［9］ZOU Y M，LIU W L，ZHOU H，et al．A new implementation procedure of sequential Gaussian simulation in stochastic seismic inversion［C］／／2013 SEG Annual Meeting，2013．

［10］HU L Y．Gradual deformation and iterative calibration of Gaussian－related stochastic models［J］．Mathematical Geology，2000，32（1）：87－108．

［11］KJ?NSBERG H，HAUGE R，KOLBJ?RNSEN O，et al．Bayesian Monte Carlo method for seismic predrill prospect assessment［J］．Geophysics，2010，75（2）：09－019．

［12］JOHN A，SCALES，MARTIN L．SMITH，SVEN TREITEL．Introductory Geophysical Inverse Theory ［M］．Golden：Samizdat Press，2001．

［13］HU L Y，BLANC G．Constraining a reservoir facies model to dynamic data using agradual deformation method［C］．6th European Conference on the Mathematics of Oil Recovery．1998．

［14］HU L Y，LE RAVALEC M，BLANC G，et al． Reducing uncertainties in production forecasts by constraining geological modeling to dynamic data［C］．SPE Annual Technical Conference and Exhibition．Society of Petroleum Engineers，1999．

［15］楊培杰．地震子波盲提取與非線性反演［D］．東營：中國石油大學(xué)（華東），2008．

YANG P J．SeismicWavelet Blind Extraction and Non －Linear Inversion［D］．Dongying：China University of Petroleum（Huadong），2008．（In Chinese）

［16］MOSEGAARD K，TARANTOLA A．Monte Carlo sampling of solutions to inverse problems［J］．Journal of Geophysical Research：Solid Earth（1978–2012），1995，100（B7）：12431－12447．

［17］GELMAN A，ROBERTS G，GILKS W．Efficient metropolis jumping hules［J］．Bayesian statistics，1996，5：599－608．

［18］ZHANG G Z，LIU Y S，YIN X Y．The Effects of Factors of Spatial Correlation Analysis On Kriging Estimate In Geostatistical Reservoir Characterization［C］．2005SEG Annual Meeting．Society of Exploration Geophysicists，2005，1457－1460．

Stochastic seismic inversion based on the geostatistical priori information

YE Duan-nan，YIN Xing－yao，SUN Rui-ying，WANG Bao-li
（School of Geosciences，China University of Petroleum，Qingdao 266580，China）

Stochastic seismic inversion based on the geostatistical priori information is a Monte Carlo based strategy for non －linear inversion．It is formulated in a Bayesian framework．Firstly，the priori information can be obtained through sequential Gaussian simulation（SGS）and gradual deformation method（GDM）．Then we can construct the likelihood function．Finally，we apply Metropolis sampling algorithm in order to obtain an exhaustive description of the posteriori probability density and get the inversion results．Compared with the deterministic inversion，the inversion method we proposed can effectively integrate the high－frequency information of well－logging data and have a higher resolution．According to the numerical calculations，the final results match the model well and have a high resolution．In addition，we use the sequential Gaussian simulation（SGS）in a new implementation way and combine with GDM，which can improve the calculation efficiency of inversion method effectively．

geostatistical priori information；Bayesian theory；high-resolution；calculation efficiency

P 631．4

A

10．3969／j．issn．1001-1749．2015．03．12

1001-1749（2015）03-0341-07

2014-02-25 改回日期：2014-07-01

國家973項(xiàng)目（2013CB228604）；國家科技重大專項(xiàng)（2011ZX05009）；山東省自然科學(xué)基金（ZR2011DQ013）；國家自然科學(xué)基金（41204085）；中國石化地球物理重點(diǎn)實(shí)驗(yàn)室（WTYJY－WX2013－04－07）

葉端南（1984－），男，博士，主要研究方向?yàn)榈卣饍?chǔ)層預(yù)測，E-mail：yedn＠upc．edu．cn。