基于地質統計先驗信息的隨機地震反演

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

（中國石油大學（華東）地球科學與技術學院，青島 266580）

基于地質統計先驗信息的隨機地震反演

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

（中國石油大學（華東）地球科學與技術學院，青島 266580）

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

地質統計先驗信息；貝葉斯理論；高分辨率；計算效率

0 引言

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

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

1 方法原理

利用貝葉斯理論尋找反演問題的解，貝葉斯理論的一個優勢就是在反演彈性參數的同時，可以得到不確定性的估計。基本思路就是在貝葉斯理論框架下，通過SGS算法得到高斯概率密度，再利用GDM擾動更新先驗信息，然后根據地質統計先驗信息和原始地震資料確定似然函數，最后利用Metropolis算法對后驗概率密度進行抽樣得到反演問題的解。圖1為基于地質統計先驗信息的隨機地震反演流程圖。

圖1 基于地質統計先驗信息的隨機地震反演流程圖Fig．1 Flowchart of stochastic seismic inversion based on the geostatistical priori information

1．1 貝葉斯理論

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

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

假設地震褶積模型為：

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

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

1．2 地質統計先驗信息

序貫高斯模擬（SGS）是常用的先驗信息構建方法，但是SGS存在計算量較大的缺點，所以采用一種新的逐點模擬方式對測井數據進行序貫高斯模擬，為獲得更多的模擬實現，結合GDM更新算法控制先驗概率密度的微擾，而不是利用SGS重復模擬得到多個實現，這在一定程度上進一步提高了計算效率。因此可以結合SGS和GDM獲得有效的先驗信息，這就是所謂的基于地質統計學的先驗信息。

1．2．1 SGS算法

序貫模擬的基本思想是沿隨機模擬路徑，將某一鄰域內的所有已知數據（原始數據和先前已模擬的數據）作為條件數據進行模擬，獲取每個網格節點的條件累積概率密度（CCDF），進而從CCDF得到模擬值。序貫高斯模擬適用于高斯模型，因此需確保原始數據服從高斯分布，若不服從，必須進行高斯變換，模擬結束后，進行高斯反變換。

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

1．2．2 GDM更新算法

傳統的隨機反演方法在進行模擬時，即使達到了Metropolis抽樣的接受條件，也不能保證模擬結果一定會使目標函數收斂，這里在利用新的SGS模擬方式提高計算效率的同時，引入了逐步變形更新算法，擾動更新SGS的實現，以保證模擬搜索時目標函數快速地收斂，從而達到提高反演精度的目的。

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

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

1．3 似然函數構建

對后驗信息進行抽樣時，接受概率的建立需要使用似然函數，似然函數采用如下形式：

其中：si代表模擬地震波形的振幅，是由彈性參數計算得到的合成記錄；siobs是觀測數據的采樣點；σ是期望數據不確定性的標準差；Ri代表反演的彈性參數；Ri0是確定性反演中井位置處的低頻約束（平滑約束和點約束）［15］；α是可調參數，當地震資料信噪比較高時，可選用較小的α值，當地震資料中含有多種噪音時，適當地增大α值。

1．4 Metropolis抽樣

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

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

在貝葉斯理論框架下，通過新的SGS方式得到高斯概率密度，然后利用GDM算法擾動更新模擬實現，就可以得到基于地質統計學的先驗信息。根據確定性反演中井位置處的低頻約束構建似然函數，最后利用Metropolis算法對后驗概率密度進行抽樣得到了反演問題的解。

2 模型試算

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

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

圖2 單道反演結果對比Fig．2 Comparison of inversion results and real data

圖3 原始數據Fig．3 The original data

3 實際資料測試

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

反演速度較快，當反演較厚的儲層時，可以應用稀疏脈沖反演；而當實際工區需要進行薄層的識別時需要應用隨機地震反演方法。

圖4 隨機地震反演結果Fig．4 Stochastic seismic inversion results

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

4 結論與認識

序貫高斯模擬是經典的地質統計學模擬方法，而GDM更新算法能夠連續更新儲層參數實現使目標函數快速收斂，并且隨機地震反演能夠有效融合測井的高頻信息，提高反演結果的分辨率，但計算效率是其面臨的一個問題。結合序貫高斯模擬和GDM更新算法進行隨機反演，新的模擬方式結合GDM有效提高了計算效率，并且通過模型試算和實際資料分析結果表明，該反演方法能夠穩定地收斂，反演結果比較可靠。

需要強調的是，基于地質統計先驗信息的隨機地震反演方法中，變差函數的影響非常重要，特別是二維和三維的反演過程中，各向異性變差函數對反演結果的影響會很明顯，可考慮通過各向異性變差函數改善隨機反演的橫向連續性。

圖6 稀疏脈沖反演與隨機反演的對比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］．石油大學學報：自然科學版，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］印興耀，劉永社．儲層建模申地質統計學整地震數據的方法及砰究進展［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］．東營：中國石油大學（華東），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項目（2013CB228604）；國家科技重大專項（2011ZX05009）；山東省自然科學基金（ZR2011DQ013）；國家自然科學基金（41204085）；中國石化地球物理重點實驗室（WTYJY－WX2013－04－07）

葉端南（1984－），男，博士，主要研究方向為地震儲層預測，E-mail：yedn＠upc．edu．cn。