布谷鳥馬爾科夫鏈蒙特卡洛混合高斯地質(zhì)統(tǒng)計學(xué)隨機(jī)反演

王峣鈞， 邢凱， 厙斌， 劉宇， 陳挺， 胡光岷， 吳秋波

1 電子科技大學(xué)資源與環(huán)境學(xué)院， 成都 6117312 太原理工大學(xué)求實學(xué)院， 太原 0300243 東方地球物理公司物探技術(shù)研究中心， 河北 涿州 072751

0 引言

地震反演是估計地下儲層參數(shù)最重要的技術(shù)之一，在儲層預(yù)測和巖相分類中發(fā)揮著至關(guān)重要的作用 (Castagna, 2001; Buland and Omre, 2003; Russell et al., 2003, 2011; Downton, 2005; González et al., 2008; Zhang, 2017; Yin and Zhang, 2014; Yin et al., 2015; Connolly and Hughes, 2016; Aleardi, 2018).地震反演問題依據(jù)求解方式不同，主要分為兩種方式：一種是基于最優(yōu)化思想的確定性反演(Aster et al., 2011; Sen and Stoffa, 2013)；另一種是基于后驗概率的非確定性反演(Tarantola, 2005).由于數(shù)據(jù)測量的不確定性和目標(biāo)函數(shù)病態(tài)問題等原因，使得非確定性反演對地震反演問題求解更具優(yōu)勢(Grana et al., 2017).

非確定性反演中應(yīng)用最為廣泛的是基于貝葉斯框架的地質(zhì)統(tǒng)計學(xué)隨機(jī)反演.該方法利用測井?dāng)?shù)據(jù)和地震數(shù)據(jù)構(gòu)建先驗分布和似然函數(shù)，通過抽樣獲取后驗概率估計未知儲層參數(shù)模型(Scales and Tenorio, 2001; Ulrych et al., 2001; Tarantola, 2005).在常規(guī)地質(zhì)統(tǒng)計學(xué)反演問題中，通常假設(shè)儲層參數(shù)先驗信息服從高斯分布.但是，實際地下儲層由多種巖性沉積物構(gòu)成，儲層參數(shù)具有多峰統(tǒng)計特征，如果對地下介質(zhì)按照相同的分布進(jìn)行預(yù)測顯然存在問題.研究表明，不同巖相下參數(shù)之間存在差異，但是相同巖相條件下參數(shù)近似符合高斯分布，可以采用混合高斯分布(Grana and Rossa, 2010; Dubreuil-Boisclair et al., 2012; Sauvageau et al., 2014; Amaliksen, 2014)近似刻畫地下儲層參數(shù)分布特征(Hastie and Tibshirani, 1996).基于這一實際情況，相關(guān)學(xué)者提出了基于混合高斯的地質(zhì)統(tǒng)計學(xué)隨機(jī)反演方法(Grana et al., 2017; Li et al., 2019)，該方法可以有效提高反演精度并輸出儲層巖相分類結(jié)果.

現(xiàn)有混合高斯地質(zhì)統(tǒng)計學(xué)隨機(jī)反演方法將不同巖相數(shù)據(jù)高斯分布以加權(quán)形式進(jìn)行整合，權(quán)系數(shù)按照測井?dāng)?shù)據(jù)中統(tǒng)計不同巖相比例確定，在反演過程中一般固定不變.這種固定初始巖相分配比例的方案可能會因測井?dāng)?shù)據(jù)分布不均和測井?dāng)?shù)據(jù)樣本選擇不全面等原因造成誤差，進(jìn)而影響最終反演結(jié)果.如果將混合高斯模型權(quán)重和儲層參數(shù)、儲層巖相都作為待反演參數(shù)進(jìn)行同步反演，必然能夠避免初始設(shè)定誤差，提升反演精度.如果巖相權(quán)重變?yōu)榇囱輩?shù)，則反演目標(biāo)函數(shù)變?yōu)榉蔷€性函數(shù)，且反演參數(shù)變多導(dǎo)致多解性變強(qiáng)，需要提供一種穩(wěn)定性更強(qiáng)、能較好抑制多解性的方法對目標(biāo)函數(shù)進(jìn)行求解.

目前隨機(jī)反演主要采用馬爾科夫鏈蒙特卡洛(MCMC)方法進(jìn)行求解，其基本流程是首先在目標(biāo)解空間內(nèi)進(jìn)行Metropolis-Hastings隨機(jī)采樣，根據(jù)采樣前后狀態(tài)求出后驗概率分布，然后利用變量的馬爾可夫性求出狀態(tài)轉(zhuǎn)移概率，當(dāng)狀態(tài)轉(zhuǎn)移概率趨于穩(wěn)定時得到最終結(jié)果.該算法主要受迭代步長、馬爾科夫鏈長度、迭代終止條件等因素影響，一旦迭代次數(shù)或馬爾科夫鏈長度選擇不合理，會導(dǎo)致采樣結(jié)果陷入局部最優(yōu)(Geyer, 2005).如果采用目前的MCMC方法進(jìn)行前文所述混合高斯模型權(quán)重和儲層參數(shù)、儲層巖相非線性同步反演，由于反演維度和參數(shù)進(jìn)一步增加，會進(jìn)一步加劇MCMC方法不穩(wěn)定性.為此，本論文引入布谷鳥搜索方法和MCMC算法進(jìn)行整合，實現(xiàn)可控步長和搜索方向的布谷鳥蒙特卡洛優(yōu)化算法(CSMCMC)，跳出局部最優(yōu)解，獲得穩(wěn)定的多參數(shù)反演結(jié)果.

布谷鳥算法是2009年提出的一種新型啟發(fā)式算法(Rajabioun, 2011; Yang and Deb, 2009,2010,2013)算法是根據(jù)布谷鳥幼雛不斷與宿主進(jìn)行適應(yīng)性斗爭，最終使得宿主相信的進(jìn)化過程所提出的一種全局優(yōu)化算法.引入布谷鳥算法，可以利用該算法中Levy飛行變步長搜索特點提升MCMC中馬爾科夫鏈的多樣性，使得抽樣具有更好的動態(tài)調(diào)節(jié)性，從而提升對非線性多參數(shù)問題求解全局優(yōu)化性.在本文工作中，我們基于混合高斯模型地質(zhì)統(tǒng)計學(xué)隨機(jī)反演框架，采用布谷鳥MCMC混合優(yōu)化算法實現(xiàn)了混合高斯分布后驗概率的求解，同步獲得混合高斯模型不同巖相比例系數(shù)、巖相和儲層參數(shù)，實現(xiàn)了儲層巖相的精確劃分，避免初始巖相比例固定導(dǎo)致的反演誤差.通過模型和實際資料驗證表明了所提出算法的有效性.

1 方法原理

1.1 混合高斯地質(zhì)統(tǒng)計學(xué)隨機(jī)反演

地震反演可以通過后驗概率分布來表示(Grana et al., 2017)：

(1)

其中p(m|d)是在地震數(shù)據(jù)d約束下儲層參數(shù)的后驗概率分布，p(d|m)是儲層參數(shù)m與地震數(shù)據(jù)之間的似然函數(shù)，p(m)是儲層參數(shù)的先驗概率分布.在混合高斯地質(zhì)統(tǒng)計學(xué)反演中先驗信息可以分不同巖相進(jìn)行描述.假設(shè)巖相可以分為K類，儲層參數(shù)先驗分布函數(shù)可以描述為不同巖相內(nèi)儲層參數(shù)先驗分布概率密度之和(Li et al., 2019)，即：

(2)

(3)

儲層參數(shù)m與地震數(shù)據(jù)d之間的似然函數(shù)反映了地震數(shù)據(jù)與儲層參數(shù)對應(yīng)函數(shù)關(guān)系，如果可以考慮為地震波正演形式，可以描述為

(4)

其中λt是與Cm相關(guān)的常量，Cm是當(dāng)前儲層參數(shù)的方差，Wt是正演算子，此處代表地震子波矩陣，D是差分矩陣，分別表示為

(5)

結(jié)合上述公式可以構(gòu)建混合高斯儲層參數(shù)的后驗概率密度函數(shù)為

(6)

將(3)式代入并取對數(shù)并化簡，可以得到后驗概率密度函數(shù)的對數(shù)形式：

(7)

上述目標(biāo)函數(shù)中，通常假設(shè)λk為固定值以便降低變量數(shù)量以通過MCMC方法進(jìn)行求解.如果λk為未知參數(shù)，該目標(biāo)函數(shù)參數(shù)量和非線性程度進(jìn)一步增加，采用MCMC算法難以實現(xiàn)穩(wěn)定求解，因此本論文提出引入布谷鳥算法與MCMC結(jié)合實現(xiàn)該反問題的求解.

1.2 布谷鳥搜索算法

布谷鳥搜索算法(Cuckoo search, CS)主要包括種群生成、Levy飛行搜索更新和適應(yīng)性選擇等幾個主要步驟，其中Levy飛行有重要作用.Levy飛行與普通隨機(jī)搜索的區(qū)別在于該方法是一種變步長的搜索方法，可以避免結(jié)果陷入局部最優(yōu).

布谷鳥算法流程如下：

①根據(jù)解空間維度隨機(jī)生成Ω個鳥巢Mi，每個鳥巢有n個個體

Mi=(m1,m2,m3,…,mn),i=1,2,…,Ω

(8)

其中，第i個鳥巢第j個個體的初始值如下：

i=1,2,…,Ω;j=1,2,…,n

(9)

其中Lj_max,Lj_min是第j維數(shù)據(jù)的最大值與最小值.此外，我們需要設(shè)置迭代次數(shù)、發(fā)現(xiàn)概率等參數(shù)；

②計算鳥巢的適應(yīng)度，得到最優(yōu)鳥巢Mbest，適應(yīng)度函數(shù)可以反映種群中個體之間性能優(yōu)劣，一般情況下選擇均方根誤差作為評判標(biāo)準(zhǔn)；

③通過Levy飛行產(chǎn)生新解，對其他鳥巢位置與狀態(tài)進(jìn)行更新.Levy飛行代表布谷鳥位置更新，公式如下：

(10)

(11)

(12)

其中u和v服從正態(tài)分布，u=randnφ,v=randn，β通常取1.5.而φ可以由下式得出：

(13)

(14)

(15)

其中σ是縮放因子，σ～U(0,1)，γt,t代表第t代所有個體中任意兩個個體.

⑤計算種群適應(yīng)度，選出新的最佳巢穴Mnewbest，將Mnewbest與上一最優(yōu)值Mbest進(jìn)行比較，如果滿足適應(yīng)度下降的接受條件，則改變當(dāng)前最優(yōu)值，并且更新種群適應(yīng)度值.

⑥如果沒有達(dá)到最大迭代次數(shù)或者沒有達(dá)到誤差降低標(biāo)準(zhǔn)，則返回③，直至輸出最優(yōu)目標(biāo)解.

1.3 布谷鳥MCMC混合高斯地質(zhì)統(tǒng)計學(xué)隨機(jī)反演

本文所提算法核心是在布谷鳥的尋找當(dāng)前最優(yōu)巢穴的過程中采用MCMC算法中的Gibbs采樣獲得多條馬爾科夫鏈，然后通過Levy飛行方式對多條馬爾科夫鏈進(jìn)行優(yōu)化，最后根據(jù)多條馬爾科夫鏈的最優(yōu)適應(yīng)度選擇輸出結(jié)果.這種多條馬爾科夫鏈同步優(yōu)化方法可以擴(kuò)大樣本范圍，提升方法全局優(yōu)化能力.此外，對于種群初始化，由于傳統(tǒng)初始化策略使得每個初始種群都與實際參數(shù)模型分布有很大的差異，導(dǎo)致反演過程中收斂速度慢，考慮不采取隨機(jī)初始化方式，而是利用地質(zhì)統(tǒng)計學(xué)建模獲得的儲層參數(shù)初始分布作為算法種群初始解將傳統(tǒng)布谷鳥搜索中種群初始化進(jìn)行優(yōu)化，使得初始化的種群仍然保留儲層參數(shù)分布的大致特點.基于布谷鳥搜索(CS)和MCMC結(jié)合的混合高斯地質(zhì)統(tǒng)計學(xué)隨機(jī)反演方法(GMM-CSMCMC)流程如圖1所示.具體流程為

圖1 基于布谷鳥搜索馬爾科夫鏈蒙特卡洛(MCMC)混合高斯地質(zhì)統(tǒng)計學(xué)隨機(jī)反演流程Fig.1 The workflow of GMM-CSMCMC algorithm

①根據(jù)地質(zhì)統(tǒng)計學(xué)建模方法，我們可以得到目標(biāo)地層的初始分布，該模型粗略地反映地層分布，使用該結(jié)果作為反演的初始模型；

②根據(jù)地質(zhì)統(tǒng)計學(xué)的結(jié)果進(jìn)行種群初始化，即馬爾科夫鏈初始化：

(16)

其中，i代表個體，j代表個體i的維度，m_ini是由地質(zhì)統(tǒng)計學(xué)建模求得的地層參數(shù)分布，var_ini是地質(zhì)統(tǒng)計學(xué)隨機(jī)建模求得的對應(yīng)位置的方差值，eLb是Lb長度之內(nèi)的隨機(jī)數(shù)，K1和K2是常數(shù)項，Lb和Ub代表下界與上界.地質(zhì)統(tǒng)計學(xué)建?？梢缘玫降貙幽Ｐ偷某跏挤植?，因此該初始化策略可以使得初始種群更加接近于實際地層參數(shù)分布；

④針對所有馬爾科夫鏈，根據(jù)適應(yīng)度函數(shù)計算最優(yōu)值.記第i個種群馬爾科夫鏈的適應(yīng)度為fitMi，最小的適應(yīng)度為bestfitM，對應(yīng)標(biāo)號為Lbestfit，則該種群當(dāng)前最優(yōu)解Mbest=M(Lbestfit).適應(yīng)度函數(shù)為

Fit(Mi)=(GMi-d)T(GMi-d)

(17)

⑥求出后驗概率分布函數(shù)為

(18)

⑦如果誤差沒有下降到目標(biāo)范圍或者迭代次數(shù)不足，則迭代次數(shù)t=t+1，重復(fù)上述步驟⑤—⑥，直至迭代次數(shù)達(dá)到最大為止，輸出彈性參數(shù)、巖相比例和巖相分類反演結(jié)果.

2 實驗分析

2.1 模型測試

首先使用40 Hz雷克子波和真實測井?dāng)?shù)據(jù)合成記錄作為實際地震數(shù)據(jù)進(jìn)行測試.根據(jù)算法參數(shù)設(shè)置經(jīng)驗，設(shè)置初始種群個數(shù)(鳥巢數(shù)量)為25，宿主發(fā)現(xiàn)異常卵的概率設(shè)置為0.25.選擇我國東北某工區(qū)3口井?dāng)?shù)據(jù)進(jìn)行反演分析.為方便比較，此處我們將權(quán)系數(shù)固定不變的常規(guī)混合高斯地質(zhì)統(tǒng)計學(xué)反演簡記為GMM-c，權(quán)系數(shù)變化但是采用MCMC算法進(jìn)行求解的方法稱為GMM-v，采用新方法求解的稱為GMM-CSMCMC，圖2為結(jié)果對比圖.

在圖2中，綠色線條代表實際井?dāng)?shù)據(jù)，藍(lán)色、淺藍(lán)色線和紅色線條分別代表GMM-c、GMM-v和GMM-CSMCMC反演結(jié)果.為了定量對比，我們還計算了3種反演方法結(jié)果與井?dāng)?shù)據(jù)相關(guān)系數(shù)和均方根誤差以驗證反演效果有效性，如表1所示.從反演結(jié)果與井?dāng)?shù)據(jù)的吻合程度來看，GMM-v反演結(jié)果要比GMM-c反演結(jié)果更優(yōu)，而GMM-CSMCMC反演結(jié)果與井?dāng)?shù)據(jù)的吻合程度要明顯優(yōu)于其他兩類反演方法.

表1 合成記錄誤差與相關(guān)系數(shù)Table 1 Error and correlation coefficient of synthetic records

圖2 GMM-c (a)、GMM-v (b)、GMM-CSMCMC (c) 反演的單道結(jié)果對比Fig.2 Comparison of single trace results of GMM-c (a), GMM-v (b), GMM-CSMCMC (c) inversion using synthetic seismic records

選取第一口井巖相反演結(jié)果做對比，如圖3所示.GMM-c、GMM-v和GMM-CSMCMC的巖相分

圖3 井?dāng)?shù)據(jù)(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的分相結(jié)果對比Fig.3 Comparison lithology results of well data (a), GMM-c inversion (b), GMM-v inversion (c), GMM-CSMCMC inversion (d)

類正確個數(shù)分別為：43，48，50，表明本文所提方法可得出更為準(zhǔn)確的巖相分類結(jié)果.

進(jìn)一步為了測試算法抗噪性，我們在合成記錄的基礎(chǔ)上添加信噪比等于10和信噪比為4的高斯噪聲.單道測試結(jié)果如圖4和5所示，三種反演結(jié)果與井?dāng)?shù)據(jù)的均方根誤差與相關(guān)系數(shù)如表2和表3所示，巖相分類結(jié)果如圖6和7所示.可以發(fā)現(xiàn)，在含噪聲數(shù)據(jù)中，采用變權(quán)值的方法(GMM-v)反演結(jié)果要更加穩(wěn)定，而GMM-CSMCMC反演結(jié)果與井?dāng)?shù)據(jù)的吻合程度要高于GMM-v，可見GMM-CSMCMC反演結(jié)果對于高斯噪聲有很好的抗噪性，而GMM-c反演結(jié)果受噪聲影響較大.從圖6和圖7巖相結(jié)果對比可以發(fā)現(xiàn)，添加SNR=10的高斯噪聲后，不同反演結(jié)果的獲得巖相分類正確個數(shù)分別為：43，45，50；在信噪比為4時不同反演結(jié)果的巖相分類正確個數(shù)分別為：46，47，50，從參數(shù)數(shù)值結(jié)果精確度與巖相分類正確個數(shù)來看，GMM-CSMCMC反演結(jié)果都是最優(yōu)的，因此上述實驗說明了在不同信噪比條件下，GMM-CSMCMC反演方法的抗噪性要優(yōu)于其他兩種方法.

圖4 添加SNR=10的高斯噪聲的合成地震記錄種GMM-c (a)、GMM-v (b)、GMM-CSMCMC (c) 反演結(jié)果對比Fig.4 Comparison of GMM-c (a), GMM-v (b) and GMM-CSMCMC (c) inversion results using synthetic seismic records with Gaussian noise (SNR=10)

表2 添加高斯噪聲(SNR=10)反演結(jié)果與真實井曲線誤差及相關(guān)系數(shù)Table 2 Error and correlation coefficient of inversion results and real well curves with Gaussian noise (SNR=10)

表3 添加高斯噪聲(SNR=4)反演結(jié)果與真實井曲線誤差及相關(guān)系數(shù)Table 3 Error and correlation coefficient of inversion results and real well curves with Gaussian noise (SNR=4)

圖5 添加SNR=4的高斯噪聲的合成地震記錄種GMM-c (a)、GMM-v (b)、GMM-CSMCMC (c) 反演結(jié)果對比Fig.5 Comparison of GMM-c (a), GMM-v (b), and GMM-CSMCMC (c) inversion results using synthetic seismic records with Gaussian noise (SNR=4)

圖6 添加SNR=10高斯噪聲下井?dāng)?shù)據(jù)(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的分相結(jié)果對比Fig.6 The comparison of the lithology results of well data (a), GMM-c inversion (b), GMM-v inversion (c), GMM-CSMCMC inversion (d) with Gaussian noise (SNR=10)

圖7 添加SNR=4高斯噪聲情況下井?dāng)?shù)據(jù)(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的分相結(jié)果對比Fig.7 The comparison of the lithology results of well data (a), GMM-c inversion (b), GMM-v inversion (c), GMM-CSMCMC inversion (d) with Gaussian noise (SNR=4)

2.2 實際資料測試

GMM-CSMCMC反演之所以能夠得到高穩(wěn)定性的結(jié)果，是因為布谷鳥算法特殊的搜索方式：在前期搜索步長較大，從而可以增加種群的多樣性，此時可以有效地跳出局部最優(yōu)解；后期搜索步長較小，可以得到某個區(qū)域內(nèi)的較高精度解，因此布谷鳥搜索是一種更方便得出全局最優(yōu)解的搜索方式.為了驗證該理論，我們選擇東部某過井地震數(shù)據(jù)，測試GMM-CSMCMC反演方法與GMM-CSMCMC反演方法的迭代誤差下降趨勢，結(jié)果如圖8所示.

根據(jù)圖8可知，GMM-c反演與GMM-v反演誤差下降趨勢呈階梯狀，而GMM-CSMCMC反演誤差一直下降直到穩(wěn)定，能跳出由于參數(shù)設(shè)置不合理等引起的局部最優(yōu)現(xiàn)象.可見，GMM-CSMCMC能夠更為快速的下降到更低的穩(wěn)態(tài)誤差，因此我們可以證明GMM-CSMCMC反演方法可以得到更準(zhǔn)確的反演結(jié)果.

為了考慮參數(shù)選擇對反演結(jié)果的影響，我們進(jìn)行算法參數(shù)測試.根據(jù)圖9分析可知，當(dāng)馬爾科夫鏈長度過短時會導(dǎo)致誤差下降不充分，相關(guān)系數(shù)變化不穩(wěn)定，隨著馬爾科夫鏈長度逐步增加，反演結(jié)果的誤差在逐漸下降，相關(guān)系數(shù)整體呈上升趨勢，長度為8時誤差與相關(guān)系數(shù)變化趨于穩(wěn)定.根據(jù)圖10分析可知，當(dāng)溫度下降至260 ℃附近，反演結(jié)果已經(jīng)趨于穩(wěn)定，說明反演過程更快達(dá)到穩(wěn)定狀態(tài).綜合分析可知，GMM-v反演融合布谷鳥算法后，減小了參數(shù)設(shè)置的約束，更容易得到合適的算法參數(shù).

圖8 GMM-c反演(a)、GMM-v反演(b)，與GMM-CSMCMC反演(c)的誤差隨迭代次數(shù)變化對比Fig.8 Error vs. iteration number comparison chart of GMM-c inversion (a), GMM-v inversion (b), and GMM-CSMCMC inversion (c)

圖9 不同長度馬爾科夫鏈GMM-CSMCMC單道反演結(jié)果的誤差(a)與相關(guān)系數(shù)(b)曲線Fig.9 Errors (a) and correlation coefficient curves (b) of inversion results with Markov chains of different lengths of GMM-CSMCMC inversion

圖10 GMM-CSMCMC反演不同溫度單道反演結(jié)果的誤差(a)與相關(guān)系數(shù)(b)曲線Fig.10 Errors (a) and correlation coefficient curves (b) of inversion results with different temperature of GMM-CSMCMC inversion

在算法中初始種群個數(shù)(鳥巢個數(shù))的不同也會導(dǎo)致結(jié)果出現(xiàn)偏差.為了保證算法效率，設(shè)置參數(shù)時不宜過大.我們選擇其中一口井的數(shù)據(jù)，分別對初始種群個數(shù)為5,10,15,20,25,30的參數(shù)進(jìn)行測試.根據(jù)圖8c所示，當(dāng)?shù)螖?shù)在200左右時誤差已經(jīng)趨于穩(wěn)態(tài)，因此此處迭代次數(shù)設(shè)置為200.對于每個初始種群個數(shù)，都記錄每次反演結(jié)果與實際測井曲線的均方根誤差，然后進(jìn)行10次平均，用來比較判斷種群個數(shù)對結(jié)果的影響，如表4所示.可以發(fā)現(xiàn)，初始種群個數(shù)過大(30)或者過小(5)，結(jié)果誤差都比較大，當(dāng)鳥巢個數(shù)為25時，結(jié)果的誤差最小.通常實際計算時初始種群個數(shù)可在25附近調(diào)節(jié).

表4 不同種群個數(shù)的結(jié)果誤差表Table 4 Results error table for different population numbers

宿主發(fā)現(xiàn)異常卵的概率不同，結(jié)果也會有所不同.因此我們在選擇最佳種群個數(shù)的基礎(chǔ)上，選擇不同的發(fā)現(xiàn)異常卵概率進(jìn)行結(jié)果測試.與初始種群個數(shù)相似，對于每個發(fā)現(xiàn)異常卵概率，進(jìn)行10次獨立測試，通過比較每個發(fā)現(xiàn)異常卵概率得到的平均誤差，來判斷最優(yōu)發(fā)現(xiàn)異常卵概率.在本次實驗中，經(jīng)過測試可以設(shè)置該參數(shù)值分別為0.25,0.5,0.75(表5)，可以發(fā)現(xiàn)，概率為0.25時，結(jié)果誤差最小.

表5 宿主發(fā)現(xiàn)異常卵不同概率的結(jié)果誤差表Table 5 Error table of results for hosts with different probability of finding abnormal eggs

根據(jù)圖8c所示，迭代次數(shù)在200左右誤差下降速度已經(jīng)很緩慢，所以在實際運用中，可以將迭代次數(shù)設(shè)置為200.因此，GMM-CSMCMC反演方法建議最佳參數(shù)設(shè)置如表6所示.

表6 反演最佳參數(shù)設(shè)置表Table 6 Optimal parameter setting table

進(jìn)一步對我國東部某實際工區(qū)進(jìn)行測試，該工區(qū)數(shù)據(jù)范圍為iline∈[1,142]，xline∈[1,110]，屬于碎屑巖儲層，地震采樣間隔1 ms，巖相主要為砂泥巖.根據(jù)測井?dāng)?shù)據(jù)統(tǒng)計可以得出平均砂巖、泥巖兩種巖相分別所占比例為λ1=0.38,λ2=0.62，并將該比例作為初始值代入算法進(jìn)行計算.由于MCMC方法屬于全局優(yōu)化算法，對初始巖相比例依賴性不大，此處初始巖相比例也可自由設(shè)置，對結(jié)果影響不大.對于三種不同的反演方法，我們設(shè)置相同的共有參數(shù)，得到的結(jié)果過井剖面如圖11所示.圖11a是過井地震剖面圖，圖11(b,c,d)分別是GMM-c、GMM-v和GMM-CSMCMC反演的阻抗結(jié)果.從剖面分析可以發(fā)現(xiàn)，三種方法都能得到相比地震數(shù)據(jù)更高分辨率的波阻抗反演結(jié)果.相較于GMM-c和GMM-v反演結(jié)果而言，GMM-CSMCMC反演剖面參數(shù)橫向展布更加均勻，信噪比更高.為了進(jìn)一步驗證本實驗結(jié)果正確性，我們選擇該剖面對應(yīng)的實際地震數(shù)據(jù)，與三種反演結(jié)果的合成地震數(shù)據(jù)進(jìn)行對比.圖12(a—c)分別為 GMM-c、GMM-v和GMM-CSMCMC反演結(jié)果的合成地震剖面.與圖11a實際地震剖面對比可看出，GMM-v和GMM-CSMCMC方法得到結(jié)果合成地震數(shù)據(jù)剖面與實際地震剖面更接近.為了定量化描述這一結(jié)論，我們計算三種結(jié)果合成數(shù)據(jù)剖面與實際地震數(shù)據(jù)剖面之間的均方根誤差，結(jié)果見表7.根據(jù)定量分析可以看出，GMM-CSMCMC反演結(jié)果的合成地震數(shù)據(jù)誤差最小，最接近真實地震剖面，表明考慮巖相比例變化并采用新方法求解得到結(jié)果更符合實際地質(zhì)情況.此外，針對三種反演方法得到巖相分類結(jié)果進(jìn)行了對比，結(jié)果如圖13(a—c)所示，分別表示GMM-c、GMM-v和GMM-CSMCMC反演的巖相分類結(jié)果，巖相比例見表8.從巖相剖面對比可以看出，GMM-CSMCMC方法得到結(jié)果巖相信噪比更高.

圖11 實際地震數(shù)據(jù)剖面(a)、GMM-c反演(b)、GMM-v反演(c)與GMM-CSMCMC反演過井剖面圖(d)Fig.11 Profile of seismic data (a), GMM-c inversion (b), GMM-v inversion (c) and GMM-CSMCMC inversion profile (d)

圖12 GMM-c反演結(jié)果的合成地震剖面(a)、GMM-v反演結(jié)果的合成地震剖面(b)與GMM-CSMCMC反演結(jié)果合成地震剖面(c)Fig.12 Synthetic seismic profile of GMM-c inversion results (a), GMM-v inversion results (b) and GMM-CSMCMC inversion results (c)

表7 三種合成地震剖面與實際地震剖面的均方根誤差表Table 7 RMSE table of synthetic seismic profiles and actual seismic profiles

圖13 GMM-c反演(a)、GMM-v反演(b)與GMM-CSMCMC反演(c)剖面巖相分類圖Fig.13 Lithofacies classification map of GMM-c inversion (a), GMM-v inversion (b) and GMM-CSMCMC inversion (c)

表8 最終巖相比例表Table 8 Final lithofacies proportion table

為進(jìn)一步定量說明問題，本文抽取上述剖面中三口井反演結(jié)果進(jìn)行分析，反演波阻抗結(jié)果如圖14所示，定量分析如表9所示.由三口井波阻抗反演對比圖可見，GMM-CSMCMC反演與井?dāng)?shù)據(jù)的吻合程度高于GMM-c反演、GMM-v反演與井?dāng)?shù)據(jù)的吻合程度.從定量對比可以發(fā)現(xiàn)GMM-CSMCMC反演結(jié)果誤差要小于GMM-c和GMM-v反演結(jié)果的誤差，相關(guān)系數(shù)也更高.為了說明巖相分類精度，選擇第一口井?dāng)?shù)據(jù)巖相反演結(jié)果進(jìn)行對比，如圖15所示，可以發(fā)現(xiàn)三種反演方法得到巖相分類正確數(shù)分別為：47，49，50，GMM-CSMCMC反演準(zhǔn)確度最高.綜上所述，從各種評價指標(biāo)分析，GMM-CSMCMC反演要優(yōu)于其他兩種反演效果，驗證了本文所提方法的有效性.

圖14 三口井GMM-c、GMM-v、GMM-CSMCMC過井反演結(jié)果對比Fig.14 Comparison of GMM-c, GMM-v, GMM-CSMCMC inversion results for three well logs

圖15 真實井?dāng)?shù)據(jù)巖相分類(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的巖相分類結(jié)果對比Fig.15 Comparison of lithology of real well data(a), GMM-c inversion(b), GMM-v inversion(c), GMM-CSMCMC inversion(d)

表9 單井反演結(jié)果與真實井誤差及相關(guān)系數(shù)Table 9 Single trace inversion results with real well errors and correlation coefficients using actual seismic data

3 結(jié)論

論文提出了基于可變巖相參數(shù)的混合高斯模型地質(zhì)統(tǒng)計學(xué)反演框架，同時輸出巖相和波阻抗反演結(jié)果.針對該反演目標(biāo)函數(shù)的高維多參數(shù)非線性求解問題，提出了布谷鳥搜索和MCMC混合的全局優(yōu)化求解算法.通過實驗對比發(fā)現(xiàn)基于變巖相的混合高斯模型地質(zhì)統(tǒng)計學(xué)反演可以得到精度更高的巖相和反演結(jié)果.論文所提出的布谷鳥蒙特卡洛馬爾可夫鏈優(yōu)化算法利用Levy飛行可實現(xiàn)可變步長搜索過程，在算法前期搜索過程中，較大步長搜索增加了種群多樣性，算法后期搜索過程中步長變小，可提高搜索精度，便于在某個小區(qū)域內(nèi)得到全局最優(yōu)解，因此可以避免結(jié)果陷入局部最優(yōu)解.綜合所提出算法和變巖相混合高斯統(tǒng)計學(xué)反演目標(biāo)函數(shù)，可得到更精確的巖相分類和儲層參數(shù)反演結(jié)果.

致謝感謝東方地球物理公司提供實際數(shù)據(jù)和指導(dǎo).