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

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

0 引言

地震反演是估計(jì)地下儲(chǔ)層參數(shù)最重要的技術(shù)之一，在儲(chǔ)層預(yù)測(cè)和巖相分類中發(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)；另一種是基于后驗(yàn)概率的非確定性反演(Tarantola, 2005).由于數(shù)據(jù)測(cè)量的不確定性和目標(biāo)函數(shù)病態(tài)問題等原因，使得非確定性反演對(duì)地震反演問題求解更具優(yōu)勢(shì)(Grana et al., 2017).

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

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

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

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

1 方法原理

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

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

(1)

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

(2)

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

1.2 布谷鳥搜索算法

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

布谷鳥算法流程如下：

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

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

(8)

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

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

(9)

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

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

③通過Levy飛行產(chǎn)生新解，對(duì)其他鳥巢位置與狀態(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代所有個(gè)體中任意兩個(gè)個(gè)體.

⑤計(jì)算種群適應(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)計(jì)學(xué)隨機(jī)反演

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

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

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

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

(16)

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

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

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

(17)

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

(18)

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

2 實(shí)驗(yàn)分析

2.1 模型測(cè)試

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

在圖2中，綠色線條代表實(shí)際井?dāng)?shù)據(jù)，藍(lán)色、淺藍(lán)色線和紅色線條分別代表GMM-c、GMM-v和GMM-CSMCMC反演結(jié)果.為了定量對(duì)比，我們還計(jì)算了3種反演方法結(jié)果與井?dāng)?shù)據(jù)相關(guān)系數(shù)和均方根誤差以驗(yàn)證反演效果有效性，如表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é)果對(duì)比Fig.2 Comparison of single trace results of GMM-c (a), GMM-v (b), GMM-CSMCMC (c) inversion using synthetic seismic records

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

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

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

進(jìn)一步為了測(cè)試算法抗噪性，我們?cè)诤铣捎涗浀幕A(chǔ)上添加信噪比等于10和信噪比為4的高斯噪聲.單道測(cè)試結(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é)果對(duì)于高斯噪聲有很好的抗噪性，而GMM-c反演結(jié)果受噪聲影響較大.從圖6和圖7巖相結(jié)果對(duì)比可以發(fā)現(xiàn)，添加SNR=10的高斯噪聲后，不同反演結(jié)果的獲得巖相分類正確個(gè)數(shù)分別為：43，45，50；在信噪比為4時(shí)不同反演結(jié)果的巖相分類正確個(gè)數(shù)分別為：46，47，50，從參數(shù)數(shù)值結(jié)果精確度與巖相分類正確個(gè)數(shù)來看，GMM-CSMCMC反演結(jié)果都是最優(yōu)的，因此上述實(shí)驗(yàn)說明了在不同信噪比條件下，GMM-CSMCMC反演方法的抗噪性要優(yōu)于其他兩種方法.

圖4 添加SNR=10的高斯噪聲的合成地震記錄種GMM-c (a)、GMM-v (b)、GMM-CSMCMC (c) 反演結(jié)果對(duì)比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é)果與真實(shí)井曲線誤差及相關(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é)果與真實(shí)井曲線誤差及相關(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é)果對(duì)比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é)果對(duì)比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é)果對(duì)比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 實(shí)際資料測(cè)試

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

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

為了考慮參數(shù)選擇對(duì)反演結(jié)果的影響，我們進(jìn)行算法參數(shù)測(cè)試.根據(jù)圖9分析可知，當(dāng)馬爾科夫鏈長(zhǎng)度過短時(shí)會(huì)導(dǎo)致誤差下降不充分，相關(guān)系數(shù)變化不穩(wěn)定，隨著馬爾科夫鏈長(zhǎng)度逐步增加，反演結(jié)果的誤差在逐漸下降，相關(guān)系數(shù)整體呈上升趨勢(shì)，長(zhǎng)度為8時(shí)誤差與相關(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ù)變化對(duì)比Fig.8 Error vs. iteration number comparison chart of GMM-c inversion (a), GMM-v inversion (b), and GMM-CSMCMC inversion (c)

圖9 不同長(zhǎng)度馬爾科夫鏈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

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

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

宿主發(fā)現(xiàn)異常卵的概率不同，結(jié)果也會(huì)有所不同.因此我們?cè)谶x擇最佳種群個(gè)數(shù)的基礎(chǔ)上，選擇不同的發(fā)現(xiàn)異常卵概率進(jìn)行結(jié)果測(cè)試.與初始種群個(gè)數(shù)相似，對(duì)于每個(gè)發(fā)現(xiàn)異常卵概率，進(jìn)行10次獨(dú)立測(cè)試，通過比較每個(gè)發(fā)現(xiàn)異常卵概率得到的平均誤差，來判斷最優(yōu)發(fā)現(xiàn)異常卵概率.在本次實(shí)驗(yàn)中，經(jīng)過測(cè)試可以設(shè)置該參數(shù)值分別為0.25,0.5,0.75(表5)，可以發(fā)現(xiàn)，概率為0.25時(shí)，結(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í)際運(yùn)用中，可以將迭代次數(shù)設(shè)置為200.因此，GMM-CSMCMC反演方法建議最佳參數(shù)設(shè)置如表6所示.

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

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

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

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

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

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

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