基于量子退火Metropolis-Hastings算法的疊前隨機反演

張廣智 趙 晨 涂奇催 劉 江 張佳佳 裴忠林

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071； ③中海石油(中國)有限公司上海分公司，上海 200030)

1 引言

地震反演可以分為確定性反演與隨機反演兩大類。與傳統的確定性反演方法相比，隨機反演對于具有薄層特征的油氣藏識別具有一定優勢[1]。隨機地震反演方法的技術關鍵是分析并擬合儲層地球物理特性的分布規律并對不同地球物理參數進行研究，以此獲得這些參數與地層巖性的關系[2]。隨機反演方法可以同時擬合實際地震觀測數據和測井數據，且可以有效地利用測井數據中所包含的高頻信息提高反演的分辨率[3,4]。

作為常用的反演方法，馬爾科夫鏈蒙特卡洛方法(MCMC)可以獲得后驗概率密度的一系列樣本，通過對這些樣本進行統計分析，可以獲得滿足要求的反演結果?；贛etropolis算法[5]，Hasting[6]提出了Metropolis-Hastings (MH)算法，是最常用的MCMC方法； Smith等[7]首次提出利用MCMC方法獲取參數反演的后驗分布； Alberto等[8]利用直流電阻率測深數據反演一維地球模型，首次將MCMC方法應用于完全非線性反演問題； Chen等[9]發現MCMC反演方法能提供大量與未知參數有關的全局信息，與確定性反演方法相比，能得到更好的估計值； 朱嵩等[10]提出了動態多鏈搜索策略，在多鏈并行的過程中，通過逐步減少鏈的數目提高算法的計算效率；張廣智等[11]研究了基于MCMC方法的疊后及疊前地震反演方法； 李遠等[12]將Haario等[13,14]提出的AM-MCMC方法引入地震反演，主要通過修改候選值的產生方式來提高算法的收斂效率； 張廣智等[15]、Pan等[16]將AM (Adaptive Metropolis)策略與DR(Delayed Rejection)策略相結合，提出了AMDR-MCMC算法，將DR策略融入AM策略，克服AM策略過度依賴初始協方差的劣勢。綜合前人研究成果，可以發現傳統的MCMC算法較為依賴初始模型及搜索策略,當計算時間有限或搜索策略設置不當時,對于一個較為復雜的參數空間,MH算法往往不能對反演參數空間進行充分的搜索。

針對該問題，前人借鑒量子退火最優化思想對傳統的MH方法進行改進。魏超等[17]依據模擬退火算法與量子蒙特卡洛理論，提出量子退火最優化算法(QA)，并利用簡單的單道模型數據進行了驗證;Alulaiw等[18]將基于QA算法的反演方法應用于實際工區。量子退火算法的核心在于通過引入Hamilton量的概念修改模擬退火的接受概率。因此，同樣對MH算法的接受概率形式進行修改，引入一個逐漸減小的正變量調整狀態接受的概率，提高算法的穩定性和收斂性。本文借鑒量子退火算法的改進措施，提出了量子退火MH算法，并將其應用于疊前隨機反演。利用測井數據，獲取反演參數的先驗信息，利用反演參數的正演關系構建似然函數，利用量子退火MH算法對后驗概率密度進行抽樣，得到最終的反演結果。

2 理論與方法

在貝葉斯理論框架下，通過先驗信息及似然函數構建與后驗概率密度相關的目標函數，利用量子退火MH算法進行反演，得到反演結果。

2.1 MH算法

MCMC方法的核心是構造一個平穩分布且與所求后驗分布相同的馬爾科夫鏈，反復迭代至平穩狀態，從而得到后驗分布的樣本，再基于這些樣本做各種統計、推斷[19]。

MH方法的主要步驟為：首先由建議分布q(xt,x*)產生一個潛在的轉移xt→x*，然后根據概率α(xt,x*)來決定是否接受。從[0,1]的均勻分布上抽取隨機數u，則馬爾科夫鏈下一時刻的狀態為[20]

(1)

式中：t表示馬爾科夫鏈的當前時刻；xt表示在時刻t的值；x*表示由建議分布產生的候選抽樣值；α(xt,x*)代表轉移核函數，常用的形式為

(2)

式中，π指馬爾科夫鏈平穩分布表達式。

當目標分布取為似然函數L(x)，且建議分布為對稱分布時，則式(2)可改寫為

(3)

因此，對于反演問題，MH算法的接受概率可改寫為

α(xt,x*)=exp{min[0,g(x*)-g(xt)]}

(4)

式中g表示所構建的目標函數。

經過多次迭代，可以獲得后驗概率分布的一系列樣本，對這些樣本進行篩選，可獲得滿足要求的反演結果。

如搜索策略設置不當，傳統的MH算法收斂速度較慢。由于搜索策略的設置需要針對不同問題具體分析，因此很難找到最佳的搜索策略。本文利用量子退火MH算法應對搜索策略的問題。

2.2 量子退火MH算法

為提高反演算法的效率與精度，參照量子退火最優化方法，改進MH算法，提出量子退火MH算法，該方法的基礎在于Hamilton量的引入。

當有外力作用于體系，此時系統的Hamilton量H為

(5)

將目標函數之差ΔE=E(m(l+1))-E(m(l))看作體系的動能H0，則

H=ΔE+CΓ(t)

(6)

式中C為常數。

量子最優化算法利用系統的Hamilton量表達式替換傳統的模擬退火算法中的目標函數之差，CΓ(t)的引入使迭代的反演結果在接近模型參數時依舊存在一定的接受概率，快速接近最優結果[21]。

量子退火算法屬于全局隨機搜索的最優化算法，能有效避免線性化反演的缺陷，具有較大的發展潛力[22]，然而其接受概率最終趨近于0，因此僅能獲得唯一的反演結果，無法對反演結果進行不確定性分析；MH算法可獲得收斂于后驗概率分布的一系列反演結果的樣本，因此可以對反演結果進行不確定性分析。因此，依據量子退火算法的改進思想將式(4)改寫為

α(xt,x*)

(7)

式中φ是一個隨著迭代次數增加而逐漸趨近于0的較小的正變量。

對于量子退火MH算法來說，φ的引入可以適當減少當前狀態不必要的轉移，使馬爾科夫鏈能夠更加快速、穩定地收斂于后驗概率密度，避免局部極值的出現。

參照傳統的MH算法，量子退火MH算法的實現流程見圖1，其中“rand”為0～1之間的隨機數。其與MH算法的主要區別在于新的量子退火MH算法采用了新的接受概率表達形式。

中國礦業大學（北京）是煤炭特色高等教育的全國重點院校，擁有我國首家以能源與安全為特色的科技園——“中關村能源與安全科技園”和“中國礦業大學留學人員創業園”，并與北京市共建能源安全產業技術研究院，組成了學校產學研用及科技成果轉化體系；另外科技園就位于校內，學生不需要高頻校外甚至異地往返，避免了學校頻繁組織交通車、組織成本高且可能存在安全等方面的問題，為創建校企合作的煤礦特色機械虛擬仿真實驗平臺提供了地理和資源的雙重優勢。

圖1 量子退火MH算法流程圖

2.3 基于量子退火MH算法的疊前反演方法

疊后地震反演只能獲得地下阻抗的相關信息，對于較為復雜的地層，傳統的疊后反演無法有效識別巖性及流體，且疊后反演是基于疊后地震數據的，忽略了部分地震信息。相比疊后反演，疊前反演具有較高的保真度，且能得到泊松比、拉梅參數以及孔隙度、泥質含量、含流體飽和度等多種參數，進而為儲層預測提供更可靠的信息[23]。

疊前地震反演的理論基礎是Zoeppritz方程，然而精確的Zoeppritz方程極其復雜，不利于反演[24]。因此，前人對其進行了簡化，提出了Zoeppritz近似方程。Aki等[25]首次提出了Aki-Richards近似公式，本文基于該近似公式進行基于隨機反演的疊前反演研究。

Aki-Richards近似公式具體形式為

(8)

要利用量子退火MH算法進行疊前反演，需要構建后驗概率表達形式，即目標函數。貝葉斯理論是其基礎，它將先驗信息通過似然函數轉化為后驗信息[26]，利用量子退火MH算法對后驗概率密度進行抽樣，便可獲得后驗概率的一系列樣本。

本文所研究的反演問題可以表述為

p(vP,vS,ρ|d)=f(vP,vS,ρ)+e

(9)

式中：d為觀測地震數據；f代表正演算子；e代表觀測噪聲。待反演參數vP、vS以及ρ的后驗概率密度可寫作

p(vP,vS,ρ|d)=p(vP,vS,ρ)·p(d|vP,vS,ρ)

(10)

對于似然函數，假設地震噪聲滿足均值為0、方差為σn的正態分布，則似然函數可表示為

(11)

式中N表示待反演參數的樣本數量。

假設待反演參數vP、vS以及ρ均服從高斯分布，且相互獨立，那么待反演參數的先驗信息可表示為

p(vP,vS,ρ)=p(vP)·p(vS)·p(ρ)

(12)

(13)

最后我們利用量子退火MH算法進行反演，最終得到收斂于后驗概率密度分布的馬爾科夫鏈，再對馬爾科夫鏈進行統計分析，即可獲得反演參數vP、vS及ρ。

3 模型測試

3.1 一維模型測試

首先利用一維模型進行反演測試，該一維模型來自實際井資料。子波選取25Hz的零相位雷克子波，設置炮檢距為100～1000m。利用炮檢距、速度、時間深度等參數通過計算得到角度信息，然后利用式(8)計算反射系數，與子波進行褶積，并加入一定噪聲(SNR=3)，最終得到觀測地震記錄。

圖2為傳統MH算法反演結果與模型數據的對比，圖3為在相同搜索策略下量子退火MH算法反演結果與模型數據的對比。可以發現當搜索策略相同時，量子退火MH算法反演結果要優于傳統的MH算法反演。

圖4為觀測地震記錄、未加噪聲的地震記錄及傳統MH算法和量子退火MH算法反演結果合成的地震記錄，由圖可見，量子退火MH算法反演結果所合成的地震記錄與實際地震記錄更加接近。這是由于Γ的引入能夠適當調整原有的接受概率，減少當前狀態不必要的轉移，使算法更加穩定地收斂于后驗概率分布。

圖2 傳統MH算法反演結果與模型數據的對比 綠線為初始模型，紅線為模型數據，藍線為反演結果 (a)vP; (b)vS; (c)ρ

圖3 量子退火MH算法反演結果與模型數據的對比 綠線為初始模型，紅線為模型數據，藍線為反演結果 (a)vP; (b)vS; (c)ρ

圖4 觀測地震記錄(a)、未加噪聲的地震記錄(b)及傳統MH算法(c)和量子退火MH算法(d)反演結果合成的地震記錄 圖c、圖d中紅線為反演結果合成的記錄，黑線為未加噪聲的合成記錄

從一維模型中選取某一位置處的采樣點，比較兩種方法在該采樣點處密度值的迭代結果的變化。圖5為分別利用傳統MH算法和量子退火MH算法在某采樣點的縱波速度值隨迭代次數的變化曲線。

圖5 分別利用傳統MH算法(a)和量子退火MH 算法(b)的縱波速度值隨迭代次數的變化

由圖5可見，傳統MH算法在迭代到20000次左右才開始在真值附近區域內采樣，共耗時1230.1s，而量子退火MH算法在4000次左右就開始在真值

附近區域內采樣，共耗時228.1s，且后者的穩定性要強于前者。這說明量子退火MH算法相比傳統的MH算法具有更高的計算效率，且反演結果更加穩定。

利用量子退火MH算法進行疊前反演，可以同時獲得同一個采樣點的縱、橫波速度以及密度的多個反演結果(圖6)，進行概率統計與不確定性估算。由圖可見，縱、橫波速度與密度反演結果的概率統計均表現為高斯分布，這與給出的先驗假設一致。

3.2 二維模型測試

為了進一步驗證基于量子退火MH算法的隨機反演方法的可行性，利用部分Marmous2二維模型進行反演測試，選取模型的時窗為1000～1600ms，地震子波選取頻率為30Hz的雷克子波，設置炮檢距為100～1000m。利用炮檢距、速度、時間深度等參數計算入射角等信息，之后利用量子MH算法進行反演。

圖7為反演結果(SNR=3)與模型數據的對比。由圖可見，即使存在一定的噪聲，反演結果與模型數據基本吻合,模型中的薄層也能很好地識別。

為了進一步說明反演結果的有效性，分別比較信噪比為1和3時，模型第61道的反演結果與模型數據(圖8)。從圖8可以發現，無論信噪比為1還是3，反演結果與模型數據均吻合較好，說明了該算法具有一定的抗噪能力。

圖9為該二維模型某采樣點縱波速度、橫波速度和密度多個反演結果的概率直方圖。由圖可見，縱、橫波速度與密度反演結果的概率統計均表現為高斯分布，這與先驗假設一致。因此，可同樣選擇縱、橫波速度和密度在某采樣點處反演的均值作為最大后驗概率估計。

圖6 某采樣點的縱波速度vP(a)、橫波速度vS(b)和密度ρ(c)的概率直方圖

圖7 反演結果與模型數據對比 (a)vP模型； (b)vP反演結果； (c)vS模型； (d)vS反演結果； (e)ρ模型； (f)ρ反演結果

圖8 vP、vS、ρ反演結果與模型數據對比 (a)SNR=1； (b)SNR=3

圖9 某采樣點縱波速度vP(a)、橫波速度vS(b)和密度ρ(c)的概率直方圖

圖10 縱波速度vP(a)、橫波速度vS(b)和密度ρ(c)反演結果

4 實際資料分析

為了測試疊前隨機反演方法在實際資料應用中的可行性，應用實際資料進行反演測試。該資料來自中國北部，最大入射角為35°左右，主要目標為多層系含油的典型復式油氣聚集區，根據測井資料解釋結果，密度可以很好地反映砂體分布。

為了提高反演的計算效率，將地震炮檢距道集分為3°～13°、14°～24°、25°～35°三個角度部分疊加的地震剖面。為分析反演方法的效果，應用基于量子退火MH算法的疊前反演方法進行儲層預測。

圖10為反演結果。由圖可見，反演結果與測井數據吻合較好，且密度反演結果能大體反映出砂體的分布特征，驗證了該反演方法對實際數據進行反演的可行性。

5 結論

(1)借鑒量子退火的思路對MH算法進行改進，模型測試結果和實際數據的分析均表明，該反演方法能夠獲得縱、橫波速度和密度參數，且收斂速度及穩定性相較于傳統的MH算法有一定的提升；

(2)量子退火MH算法基于貝葉斯理論，將先驗信息與似然函數相結合，提高了反演結果的穩定性，同時可獲得收斂于后驗概率分布的一系列樣本，便于進行解的不確定分析；

(3)量子退火MH方法的核心在于φ的引入，它能夠適當調整原有的接受概率，減少當前狀態不必要的轉移，使馬爾科夫鏈能夠更加快速穩定地收斂于后驗概率密度，但是φ的選取原則及衰減函數還需進一步探索。

[1] 張繁昌,肖張波,印興耀.地震數據約束下的貝葉斯隨機反演.石油地球物理勘探,2014,49(1):176-182. Zhang Fanchang,Xiao Zhangbo,Yin Xingyao.Baye-sian stochastic inversion constrained by seismic data.OGP,2014,49(1):176-182.

[2] 李寧.基于模擬退火的地質統計學反演方法研究[學位論文].山東青島：中國石油大學(華東),2013.

[3] 孫瑞瑩.先驗信息構建與地震隨機反演方法研究[學位論文].山東青島：中國石油大學(華東),2015.

[4] 王保麗,孫瑞瑩,印興耀等.基于Metropolis抽樣的非線性反演方法.石油地球物理勘探,2015,50(1):111-117. Wang Baoli,Sun Ruiying,Yin Xingyao et al.Nonli-near inversion based on Metropolis sampling algorithm.OGP,2015,50(1):111-117.

[5] Metropolis N,Rosenbluth A,Rosenbluth M et a1.Equation of state calculations by fast computing machines.The Journal of Chemical Physics,1953,21(6):1087-1092．

[6] Hastings W K.Monte Carlo sampling methods using Markov chains and their applications. Biometrika,1970,57(1):97-109.

[7] Smith A F and Roberts G O.Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods.Journal of the Royal Statistical Society,1993,55(1):3-23.

[8] Alberto M，Carlos T V.Bayesian inversion of DC electrical measurements with uncertainties for reservoir monitoring.Inverse Problems,2000,16(5):1343-1356.

[9] Chen J S,Kemna A,Hubbard S S.A comparison between Gauss-Newton and Markov-chain Monte Carlo-based methods for inverting spectral induced-polarization data for Cole-Cole parameters.Geophysics,2008,73(6):F247-F259.

[10] 朱嵩,毛根海,劉國華等.改進的MCMC方法及其應用.水力學報,2009,40(8):1019-1023. Zhu Song,Mao Genhai,Liu Guohua et al.Improved MCMC method and its application.Journal of Hydraulic Engineering,2009,40(8):1019-1023.

[11] 張廣智,王丹陽,印興耀.利用MCMC方法估算地震參數.石油地球物理勘探,2011,46(4):605-609. Zhang Guangzhi,Wang Danyang and Yin Xingyao．Seismic parameter estimated using Markov Chain Monte Carlo method.OGP,2011,46(4):605-609.

[12] 李遠,張廣智.基于改進MCMC的地震參數反演方法.中國地球物理學會第二十九屆年會,2013.

[13] Harrio H,Saksman E and Tamminen J.Adaptive proposal distributions for random walk Metropolis algorithm.Computational Statistics,1999,14(3):375-396.

[14] Harrio H,Saksman E and Tamminen J.An adaptive Metropolis algorithm.Bernoulli,2001,7(2):223-242.

[15] 張廣智,潘新朋,孫昌路等.縱橫波聯合疊前自適應MCMC反演方法.石油地球物理勘探,2016,51(5):938-946. Zhang Guangzhi,Pan Xinpeng,Sun Changlu et al.PP- & PS-wave prestack nonlinear inversion based on adaptive MCMC algorithm.OGP,2016,51(5):938-946.

[16] Pan Xinpeng,Zhang Guangzhi.Zeroppritz-based non-linear AVO inversion using AMDR-MCMC method.SEG Technical Program Expanded Abstracts,2016,35:572-576.

[17] 魏超,朱培民,王家映.量子退火反演的原理和實現.地球物理學報,2006,49(2):577-583. Wei Chao,Zhu Peimin,Wang Jiaying.Quantum annealing inversion and its implementation.Chinese Journal of Geophysics,2006,49(2):577-583.

[18] Alulaiw B,Sen M K.Prestack seismic inversion by quantum annealing: application to Cana Field.SEG Technical Program Expanded Abstracts,2015,34:3507-3511.

[19] 李遠.基于AM-MCMC的地震反演方法研究[學位論文].山東青島：中國石油大學(華東),2014.

[20] 王丹陽.基于MCMC方法的疊前反演方法研究[學位論文].山東青島：中國石油大學(華東),2012.

[21] 魏超,李小凡,張美根.量子退火最優化與地球物理反演方法.地球物理學進展,2007,22(3):785-789. Wei Chao,Li Xiaofan,Zhang Meigen.Quantum annealing optimization and geophysical inverse method.Progress in Geophysics,2007,22(3):785-789.

[22] 方中于,王麗萍,杜家元等.基于混合智能優化算法的非線性AVO反演.石油地球物理勘探,2017,52(4):797-804. Fang Zhongyu,Wang Liping,Du Jiayuan et al.Nonlinear AVO inversion based on hybrid intelligent optimization algorithm.OGP,2017,52(4):797-804.

[23] 李建華,劉百紅,張延慶等.疊前AVO反演在儲層含油氣性預測中的應用.石油地球物理勘探,2016,51(6):1180-1186. Li Jianhua,Liu Baihong,Zhang Yanqing et al.Oil-bearing reservoir prediction with prestack AVO inversion.OGP,2016,51(6):1180-1186.

[24] 張璐.基于巖石物理的地震儲層預測方法應用研究[學位論文].山東青島：中國石油大學(華東),2009.

[25] Aki K,Richards P G.Quantitative Seismology: Theory and Methods.WH Freeman & Co,San Francisco，1980.

[26] V.Vapnik著;張學工譯.統計學習理論的本質.北京:清華大學華夏出版社,2000.

[27] 潘新朋.優化MCMC方法在地震反演中的應用研究[學位論文].山東青島：中國石油大學(華東),2016.