布谷鳥馬爾科夫鏈蒙特卡洛混合高斯地質統計學隨機反演

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

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

0 引言

地震反演是估計地下儲層參數最重要的技術之一，在儲層預測和巖相分類中發揮著至關重要的作用 (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).地震反演問題依據求解方式不同，主要分為兩種方式：一種是基于最優化思想的確定性反演(Aster et al., 2011; Sen and Stoffa, 2013)；另一種是基于后驗概率的非確定性反演(Tarantola, 2005).由于數據測量的不確定性和目標函數病態問題等原因，使得非確定性反演對地震反演問題求解更具優勢(Grana et al., 2017).

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

現有混合高斯地質統計學隨機反演方法將不同巖相數據高斯分布以加權形式進行整合，權系數按照測井數據中統計不同巖相比例確定，在反演過程中一般固定不變.這種固定初始巖相分配比例的方案可能會因測井數據分布不均和測井數據樣本選擇不全面等原因造成誤差，進而影響最終反演結果.如果將混合高斯模型權重和儲層參數、儲層巖相都作為待反演參數進行同步反演，必然能夠避免初始設定誤差，提升反演精度.如果巖相權重變為待反演參數，則反演目標函數變為非線性函數，且反演參數變多導致多解性變強，需要提供一種穩定性更強、能較好抑制多解性的方法對目標函數進行求解.

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

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

1 方法原理

1.1 混合高斯地質統計學隨機反演

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

(1)

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

(2)

(3)

儲層參數m與地震數據d之間的似然函數反映了地震數據與儲層參數對應函數關系，如果可以考慮為地震波正演形式，可以描述為

(4)

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

(5)

結合上述公式可以構建混合高斯儲層參數的后驗概率密度函數為

(6)

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

(7)

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

1.2 布谷鳥搜索算法

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

布谷鳥算法流程如下：

①根據解空間維度隨機生成Ω個鳥巢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維數據的最大值與最小值.此外，我們需要設置迭代次數、發現概率等參數；

②計算鳥巢的適應度，得到最優鳥巢Mbest，適應度函數可以反映種群中個體之間性能優劣，一般情況下選擇均方根誤差作為評判標準；

③通過Levy飛行產生新解，對其他鳥巢位置與狀態進行更新.Levy飛行代表布谷鳥位置更新，公式如下：

(10)

(11)

(12)

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

(13)

(14)

(15)

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

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

⑥如果沒有達到最大迭代次數或者沒有達到誤差降低標準，則返回③，直至輸出最優目標解.

1.3 布谷鳥MCMC混合高斯地質統計學隨機反演

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

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

①根據地質統計學建模方法，我們可以得到目標地層的初始分布，該模型粗略地反映地層分布，使用該結果作為反演的初始模型；

②根據地質統計學的結果進行種群初始化，即馬爾科夫鏈初始化：

(16)

其中，i代表個體，j代表個體i的維度，m_ini是由地質統計學建模求得的地層參數分布，var_ini是地質統計學隨機建模求得的對應位置的方差值，eLb是Lb長度之內的隨機數，K1和K2是常數項，Lb和Ub代表下界與上界.地質統計學建模可以得到地層模型的初始分布，因此該初始化策略可以使得初始種群更加接近于實際地層參數分布；

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

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

(17)

⑥求出后驗概率分布函數為

(18)

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

2 實驗分析

2.1 模型測試

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

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

表1 合成記錄誤差與相關系數Table 1 Error and correlation coefficient of synthetic records

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

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

圖3 井數據(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的分相結果對比Fig.3 Comparison lithology results of well data (a), GMM-c inversion (b), GMM-v inversion (c), GMM-CSMCMC inversion (d)

類正確個數分別為：43，48，50，表明本文所提方法可得出更為準確的巖相分類結果.

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

圖4 添加SNR=10的高斯噪聲的合成地震記錄種GMM-c (a)、GMM-v (b)、GMM-CSMCMC (c) 反演結果對比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)反演結果與真實井曲線誤差及相關系數Table 2 Error and correlation coefficient of inversion results and real well curves with Gaussian noise (SNR=10)

表3 添加高斯噪聲(SNR=4)反演結果與真實井曲線誤差及相關系數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) 反演結果對比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高斯噪聲下井數據(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的分相結果對比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高斯噪聲情況下井數據(a)、GMM-c反演(b)、GMM-v反演(c)、GMM-CSMCMC反演(d)的分相結果對比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反演之所以能夠得到高穩定性的結果，是因為布谷鳥算法特殊的搜索方式：在前期搜索步長較大，從而可以增加種群的多樣性，此時可以有效地跳出局部最優解；后期搜索步長較小，可以得到某個區域內的較高精度解，因此布谷鳥搜索是一種更方便得出全局最優解的搜索方式.為了驗證該理論，我們選擇東部某過井地震數據，測試GMM-CSMCMC反演方法與GMM-CSMCMC反演方法的迭代誤差下降趨勢，結果如圖8所示.

根據圖8可知，GMM-c反演與GMM-v反演誤差下降趨勢呈階梯狀，而GMM-CSMCMC反演誤差一直下降直到穩定，能跳出由于參數設置不合理等引起的局部最優現象.可見，GMM-CSMCMC能夠更為快速的下降到更低的穩態誤差，因此我們可以證明GMM-CSMCMC反演方法可以得到更準確的反演結果.

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

圖8 GMM-c反演(a)、GMM-v反演(b)，與GMM-CSMCMC反演(c)的誤差隨迭代次數變化對比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單道反演結果的誤差(a)與相關系數(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反演不同溫度單道反演結果的誤差(a)與相關系數(b)曲線Fig.10 Errors (a) and correlation coefficient curves (b) of inversion results with different temperature of GMM-CSMCMC inversion

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

表4 不同種群個數的結果誤差表Table 4 Results error table for different population numbers

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

表5 宿主發現異常卵不同概率的結果誤差表Table 5 Error table of results for hosts with different probability of finding abnormal eggs

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

表6 反演最佳參數設置表Table 6 Optimal parameter setting table

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

圖11 實際地震數據剖面(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反演結果的合成地震剖面(a)、GMM-v反演結果的合成地震剖面(b)與GMM-CSMCMC反演結果合成地震剖面(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

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

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

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

表9 單井反演結果與真實井誤差及相關系數Table 9 Single trace inversion results with real well errors and correlation coefficients using actual seismic data

3 結論

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

致謝感謝東方地球物理公司提供實際數據和指導.