大地電磁法的1D無偏差貝葉斯反演

史學東，柳建新，郭榮文，童孝忠，劉文劼，曹創華

大地電磁法的1D無偏差貝葉斯反演

史學東1，柳建新2，郭榮文2，童孝忠2，劉文劼2，曹創華2

（1．中原油田對外經濟貿易總公司，濮陽 457001；2．中南大學 有色資源與地質災害探查湖南省重點實驗室，長沙 410083）

應用貝葉斯理論對一維（1D）大地電磁反演問題進行無偏差不確定度分析。在貝葉斯理論中，測量數據和先驗信息包含在后驗概率密度函數（PPD）中，它可以解釋成模型的單點估計和不確定度等貝葉斯推斷，這些信息的獲取需要對反演問題進行優化求最優模型和在高維模型空間中對PPD進行采樣積分。采樣的完全、徹底和效率，對反演結果有著重要的影響。為了使采樣更有效、更完全，數值積分采用主分量參數空間的Metropolis Hastings采樣，并采用了不同的采樣溫度。在反演中，同時采用了欠參數化和超參數化方法，數據誤差和正則化因子被當成隨機變量。反演結果得到各參數的不確定度、參數間的相關關系和不同深度模型的不確定度分布。COPROD1數據的反演結果表明模型空間中存在雙峰結構。非地電參數在反演中得到了約束，說明數據本身不僅包含地球物理模型信息（電導率等），還包含了這些非地電參數的信息。

貝葉斯反演；采樣；大地電磁法

0 前言

一維（1D）大地電磁法反演問題是通過地表測量的電磁場，估計層狀大地的地電模型，在地球物理界已經獲得了廣泛的關注［1～5］。大部份大地電磁反演方法是基于線性化參數估計［6～8］，而很少有人關注非線性分析和不確定度分析。貝葉斯推斷理論為這些研究提供了天然的框架，可以有效地解決反演結果的不確定估計。

非線性貝葉斯反演需要優化求最佳起始模型和對整個模型空間進行積分。Monte Carlo積分是在積分域上隨機采樣，它可以提供無偏差的積分估計，但收斂非常慢。多MAP采樣通過多次運行優化算法來收集樣本［9］，它采樣效率高，但得到的積分估計結果有偏差。Markov－chain Monte Carlo方法（MCMC），比如Metropolis Hastings采樣和Gibbs采樣，由于它是漸近地對PPD本身進行采樣，提供了無偏差積分估計［10、11］。采樣效率和完整性是MCMC方法的兩個重要方面，其目的是獲得“充分混合”的Markov鏈，以有效地對參數空間進行采樣，避免小的、無效擾動和大的、高淘汰率的擾動。特別是對于相關性強的參數空間，由于它的高概率區傾斜于參數坐標，會導致樣本的“不充分混合”。

參數非線性不確定度的量化估計，需要知道數據的誤差信息和合理的模型參數化信息。特別是對于擬合值的計算，需要已知準確的方差信息，但在實際中卻無法獲得。在貝葉斯框架下，未知的方差及正則化因子可以通過優化近似得到（“經驗”貝葉斯方法），或把它看成隨機變量包含在采樣過程中（“分級”貝葉斯方法）。后者的優點是可以把方差的影響也包含在模型的不確定度中［12］。

Tarits et al．［13］采用Monte Carlo積分方法進行1D大地電磁法的貝葉斯反演，以獲得每層的厚度與電導率的邊緣概率分布（通過對限定先驗區間的剖分）。得到的邊緣概率分布是單峰的Gaussian分布，但也有些強烈的非Gaussian分布。Grandis et al．［14］把反演模型看成是一系列的固定厚度的層狀介質組成，并加入了平滑先驗信息，采用Gibbs采樣方法對PPD進行采樣。在反演過程中，數據誤差看成是隨機變量，通過若干不同的正則化因子的選擇，研究了正則化因子對反演的影響。Cerv et al．［15］用三種不同方法對一維大地電磁反演問題的不確定度進行估計，包括多MAP估計［9］，Gibbs采樣和Neighbourhood算法［16］。他們指出雖然MAP估計和Neighbourhood采樣效率高，但是產生的不確定度估計是有偏差的。Gibbs采樣提供了無偏差采樣，但是Cerv et al．［15］指出對于強相關參數空間，他們的算法效率低，甚至有可能導致不收斂。

1 貝葉斯理論

這里簡單地描述貝葉斯反演基本理論和方法，介紹貝葉斯反演推斷信息。m表示未知的M維模型參數空間，d表示N維的數據空間，在貝葉斯反演中都看成是隨機變量。為了更具一般性，假設模型參數化θ也是未知隨機變量，它們之間的關系通過貝葉斯框架進行關聯。

其中P（m｜d，θ）表示m對d和θ的條件概率密度函數；P（d｜m，θ）表示d對m和θ的條件概率密度函數；P（d，θ）和P（m，θ）分別表示數據和模型與θ的聯合先驗概率密度函數。在實際應用中，如果數據d是預先給定，參數化θ也是預先確定的，為了表示方便把θ省略，這時P（d）是常量。P（d｜m）可以理解為已知數據d條件下隨m變化的函數，也稱為似然函數，表示成L（m）。此時，式（1）可寫成式（2）。

其中P（m｜d）表示后驗概率密度函數。

通常似然函數可以寫成L（m）∝exp［-E（m）］，其中E（m）表示數據擬合函數（將在后面討論）。因此，后驗概率密度函數可以寫成式（3）。

其中 積分區域v表示多維模型空間；φ（m）表示包括了數據擬合和先驗信息的廣義數據擬合函數，其表達式見式（4）。

在貝葉斯反演中，式（3）中表示的多維空間后驗概率密度函數，可以用來表述一個反演問題所有解的信息。從模型的后驗概率密度函數中提取模型信息，需要計算它的模型估計，參數不確定度以及參數間的相關度等屬性。例如，最大后驗概率估計模型，概率期望模型，模型協方差Cm，1D和2D邊緣概率分布P（mi｜d）和P（mi，mj｜d），以及相關度矩陣Rm，它們的定義如下

其中δ表示Dirac delta函數；T表示轉置符號。

要提取以上信息，需要用優化方法求解好的初始模型，并對整個模型區域進行積分。對于線性反演問題，如果已知數據誤差服從Gaussian分布，且先驗信息服從均勻概率分布或Gaussian分布，那么PPD本身也將服從Gaussian分布，以上屬性問題存在解析解。對于非線性反演問題，可以采用線性化過程和線性反演方法來求解，其近似程度取決于反演問題的非線性程度。另一種替代方法就是通過數值方法求解，它沒有線性化近似誤差，具有普遍性，但是需要更高強度的計算量為代價。

對于大多數反演問題，一個模型的合理參數化事先是不知道的。在這種情況下，兩種常用的方法可以用于參數化反演模型：①欠參數化方法。該方法是用數據所能解決的最少參數（比如層數）來代表一個模型，越簡單越好；②超參數化方法。該方法用超出數據解決能力參數量（比如層數）來參數化一個模型，但在求解過程中強加了某種期望的結構來獲得穩定的解，比如平整結構、平滑結構等。也可以采用Green用于1995年提出的把參數化看成是一個隨機變量包含在后驗概率密度函數中［17］，讓數據本身去決定參數化問題—可逆跳躍Markov Chain Monte Carlo采樣。作者在本文中主要考慮兩種常用的參數化反演方法，并且在反演中，數據誤差（Cd＝sc2diag（s12，s22，…，sN2）中si表示數據誤差的粗略估計，sc表示不相關數據誤差的矯正或縮放因子［18、19］）及正則化因子看成是未知量。詳細的理論和方法可參考文獻［17～19］，由于篇幅限制，本文不再累述。

2 Metroplis－Hastings采樣數值積分

模型的后驗概率分布包含了所有的反演信息，貝葉斯反演就是從概率后分布中提取有用的模型信息。這些模型信息包括式（6）～式（8）所定義的參數不確定度和相關度信息等屬性，提取這些信息需要對后概率分布進行數值積分。將以上積分寫成一般形式表示為式（11）。

其中f（m）表示一般函數。

式（11）定義的積分可以通過解析方法或數值方法求解，作者主要考慮采用基于Monte Carlo的積分求解法。在標準Monte Carlo積分中［20］，模型是隨機從定義在積分區域上的均勻分布中產生。給定一個含Q個模型的樣本，每個樣本的概率可通過式（12）計算。

其中i＝1，…，Q

積分式（11）通過式（13）進行估計。

其中V代表M－D積分體積。

如果積分式集中在模型空間的某些局部區域，那么從均勻分布上進行采樣，其效率不高。可以采用重點采樣解決均勻采樣效率低的問題，該方法對積分貢獻較大的區域進行重點采樣。假設g（m）代表重點采樣的采樣分布函數，Q個模型從該分布上產生，其歸一化條件見式（14）。

式（11）可以寫成式（15）。

式（13）的Monte Carlo積分代表重點采樣的一個特例，即在均勻分布g（m）＝1／V上的重點采樣。需要指出的是，多重MAP積分法［9］的模型樣本產生于優化過程，其采樣函數是未知的（如模擬退火法和遺傳算法），因此沒有、也無法對式（15）中的除數g（m）進行校正。采用這種采樣方法即使對模型空間進行徹底采樣，計算出的積分估計值與實際積分值也會產生潛在的、嚴重的偏差。

基于Markov chain Monte Carlo方法（MCMC）的Metropolis－Hastings采樣（MHS［10，21，22］）的采樣函數為g（m）＝P（m｜d），把g（m）代入式（15）中，積分估計可以簡化成式（16）。

因此，通過MHS采樣，式（11）的積分簡化成在樣本集上計算f（m）的平均值。

雖然基于MCMC采樣的最終積分結果與采樣的起始模型沒關系，但一個好的起始模型（最大概率模型）可以加速采樣的收斂，例如采用線性化或非線性反演方法獲得的MAP估計作為起始模型。采樣的收斂性判斷有很多方法［23］，例如通過對比兩組或更多組平行的樣本的積分估計（或最大邊緣概率分布）：當幾個樣本的積分估計差小于某個預設的閥值時，認為采樣過程已經收斂，最后把所有的樣本合并成一組樣本進行積分估計。

3 主分量空間的擾動

如果模型空間的相關性比較強，那么基于算法（Markov Chain Monte Carlo）的優化方法和采樣方法的效率將受到嚴重的影響。在這種情況下，可以將模型空間轉換到主分量空間中以減輕這種影響。物理空間與主分量參數空間的轉換，可通過式（17）實現。

其中U是協方差矩陣Cm的列特征向量，其分解如式（18）所示。

其中 Λ＝diag［λ1，…，λM］是特征值矩陣；λi表示投影到向量ui上的方差。

由于反演開始無法獲得全局協方差矩陣Cm的信息，作者在本文采用了局部近似方法。該方法基于非線性反演問題的線性化估計，可以為高效的非線性采樣算法提供有用的局部信息［23］。Cm的線性化近似為式（19）。

其中J是在MAP解（初始采樣模型）處的Jacobian矩陣。

在“burn－in”采樣過程中，協方差矩陣Cm的起始線性化估計自適應地被基于MHS采樣的非線性估計代替，它可以更好地代表參數空間的總體協方差（這個階段的采樣沒有保存用于積分估計）。

4 采樣溫度

以上討論的數值積分是基于T＝1時的MHS采樣，然而對于多峰PPD反演問題，存在一些被低概率區隔開的高概率區，對這些區域的采樣可能不能成功地在模型空間中來回跳躍，也就是說部份模型空間被“凍結”。為了更好、更徹底地對模型空間進行采樣，有必要采用T＞1，然后將樣本矯正到T＝1。參考文獻［23］式（11）可以寫成式（20）。

其中Z1和ZT分別表示溫度為“1”和“T”時的歸一化量（也叫剖分函數）。

假設一個由Q個模型組成的樣本，服從exp［－φ（m）／T］定義的分布（等價于在溫度“T”時的MHS采樣），那么式（21）可變為式（22）。

利用式（22），任意溫度的MHS采樣都可以通過矯正得到T＝1的積分估計，使采樣更完全。當T＝1時，式（22）可簡化為標準的MHS采樣，也就是式（16）。當T→∞時，式（22）變成均勻分布的標準Monte Carlo采樣，見式（13），它是最一般的采樣方式。當T增大時，采樣更完全，更多的時間在低密度區域跳躍。對于單峰模型，高溫時的采樣相對低溫時收斂更慢，采樣時間更長；對于多峰模型，如果低溫時采樣不能徹底在整個空間中跳躍，高溫采樣將比低溫采樣效率更高。如果溫度的增加，使MHS采樣的邊緣概率分布或聯合概率密度發生了明顯的變化，這說明在更低溫度的采樣不完全，有些模型區域被“凍結”；反之，如果溫度的變化，使邊緣概率分布變化不明顯，就說明采樣完全，積分估計可靠。作者采用了一系列溫度T＝1、2、3進行MHS采樣，以確保采樣的完全性、徹底性。

5 反演結果分析

5．1 理論模型反演

作者首先考慮五層模型的貝葉斯反演，在反演中假設數據誤差（sc）是未知的，并考慮欠參數化和超參數化兩種方法。五層模型結構如表1所示，其深部存在高電導率層。在0．002 5 s～25 s間產生二十五個等對數間隔周期的合成數據，并加入百分之二的Gaussian噪聲。所有參數的先驗信息都是這樣給出的：每層的電導率和厚度分別服從定義在［0，1］S／m和［100，10000］m上的均勻分布。

圖1（見下頁）顯示在欠參數化下，五層模型參數的邊緣概率分布。圖1中顯示第一層的電導率和厚度、第二層和第五層的電導率以及sc的不確定度窄，說明在反演中能很好地得到解決。而其它參數不確定度稍微更寬，對數據響應的靈敏度相對更差。部份參數（包括誤差信息）的2D邊緣概率密度分布顯示在圖2（見下頁）中，其中h2和σ2、σ3、σ4、h1、h3之間存在著較強的相關性。圖2中還顯示出大部份參數間存在非線性關系，但是這種非線性比較微弱，可以近似為線性。

見后面，圖3和圖4分別顯示了在不同參數化下的邊緣概率剖面。在圖3中，欠參數化反演結果顯示除了第四層低阻層外，其它幾層的參數在反演中都得到了較好的約束，而第四層不確定度大。同時在經過反演，其結果還是可以區分第四層與背景地電結構的差異。圖4中的超參數化反演結果表明，淺部參數的不確定度相對深部參數的小。由于在反演中加入了平滑信息，邊緣概率剖面也反映現強加的結構信息。

5．2 實測數據反演

這里將對大地電磁法COPROD1實測數據［24、25］進行反演，這些數據已廣泛被其他研究者用作反演解釋［14、15、26］，測量周期為28．5 s～1 960．7 s間的十五個周期。為了使反演更具一般性，用Constable等［26］提供的數據誤差作為起始標準方差，實際誤差為標準方差乘上方差比例因子sc，電導率用對數坐標表示。所有參數的先驗信息是這樣給出：電導率和厚度分別服從定義在［10－6，1］S／m和［1，1000］km上的均勻分布。

反演開始，采用貝葉斯信息準則（BIC）確定反演模型的層數，作者選用四層模型，其具體過程可見文獻［18］。后面圖5顯示在不同溫度下四層模型的邊緣概率分布，每個小圖從上到下，依次表示溫度T＝1、T＝2和T＝3時的反演結果。它們的邊緣概率分布沒有明顯的差異，說明對模型空間的采樣已經完全、徹底。從參數的邊緣概率分布上，可以觀察到h3和σ4的非線性反演結果是雙峰模型，σ1和h2的邊緣概率分布較寬。h3和σ4的兩個截然不同的高概率分布區為：一個對應于峰值在更小的h3和更低的σ4處，為了表示方便簡稱“1區”；另一個對應于峰值在更大的h3和更高的σ4處，簡稱“2區”，模型參數間的相關度信息顯示在圖6（見后面）中。圖6中顯示有些參數間存在強烈的相關性，比如強負相關的有h2、σ2及h1、h2等，強正相關的有h1、σ2及h3、σ4等。

圖7（見下頁）是在欠參數化下COPROD1數據反演的邊緣概率剖面。與圖5對應，也存在彼此不相連的高概率區。有趣的是先前的COPROD1數據反演結果［15，25］與圖5所示的“2區”基模型相似（本文第一層電導率偏小，這是由先驗區間過寬造成的），它是從“2區”產生的，而“1區”的基模型具有更小的廣義數據擬合值，代表了MAP解。以往的反演結果（四層模型）都是單模結構，這可能是因采樣的不完全而引起的。

對于超參數化反演方法，在0 km～103km上，反演模型被剖分成固定厚度，成對數增加的五十層。圖8（見下頁）顯示在超參數化下反演的邊緣概率剖面，淺部不確定度大，說明數據對該深度范圍的結構不敏感，中間段模型結構在反演中很好地得到約束（約小于1 000 km），再往深部誤差越大。下頁圖9顯示μ和sc2的邊緣概率分布。sc2的1D邊緣概率分布服從期望的χ2分布，它的最大概率處被約束在“1”附近。μ在反演中得到了較好的約束（最大概率處窄），μ和sc2的2D邊緣概率分布顯示在下頁圖9（右圖部份），很明顯它們之間存在較強的相關性。

6 結論

我們應用貝葉斯理論對1D大地電磁法問題進行無偏差貝葉斯反演，計算了模型（電導率和層厚）的邊緣概率分布，并綜合邊緣概率分布計算出電導率～深度不確定度剖面，不確定度估計是通過對PPD的MCMC采樣得到。為了克服采樣過程中由于參數間的強相關性而導致的收斂慢及采樣不完全等問題，作者在本文采用了主軸空間的MHS采樣。在采樣過程中，采用了不同的采樣溫度以保證采樣的全局性、完整性，不同溫度下的采樣通過矯正得到單位溫度的無偏差積分估計。在反演中，考慮了欠參數化和超參數化反演，并把數據誤差和或正則化因子看作是未知隨機變量，最終得到該模型的無偏差不確定度分布。

作者在文中分別對未知數據誤差條件下的五層模型及實測數據進行了無偏差不確定度分析。從無偏差貝葉斯反演的邊緣概率分布，得出各參數不確定度和參數間的相關關系，以及不同深度模型的不確定度分布。不同溫度下的COPROD1數據反演結果表明，在T＝1時采樣徹底，反演結果得到雙峰結構。非地電參數當成隨機變量在反演中得到了約束，說明數據本身不僅包含地球物理模型信息（電導率），還包含了這些非地電參數的信息。

［1］TIKHONOV A N．On the determination of the electric characteristics of deep layers of the Earth’s crust［J］．Doklady Akademii Nauk SSSR，1950，151（2）：295．

［2］PARKER R L．The inverse problem of electrical conductivity in the mantle［J］．Geophysical Journal of the Royal Astronomical Society，1970，22（1），38．

［3］WEIDELT P．The inverse problem of geomagnetic induction［J］．Z．Geophys．，1972（38）：257．

［4］ROKITYANSKY I I．Geoelectromagnetic Investigation of the Earth’s Crust and Mantle［M］．New York：Springer，1982．

［5］WHITTALL K P，Oldenburg D W．Inversion of magnetotelluric data for a one－dimensional conductivity Geophys Monogr．Ser．5 ed．Fitterman D V，Society of Exploration［A］．OK：Geophysicists Tulsa，1992．

［6］BECHER W D，Sharpe C B．Asynthetic approach to magnetotelluric exploration［J］．Radio Science，1969（3）：1089．

［7］OLDENBURG D W．One－dimensional inversion of natural－source magnetotelluric observations［J］．Geophysics，1979（44）：1218．

［8］SMITH J T，BOOKER J R．Magnetotelluric inver－sion for minimum structure．Geophysics［J］．1988，53（12）：1565．

［9］SEN M K，STOFFA P L．Global Optimization Methods in Geophysical Inversion［M］．Amsterdam：Elsevier，1995．

［10］GILKS W R，RICHARDSON S，SPIEGELHALTER G J．Markov Chain Monte Carlo in Practice［M］．London：Chapman and Hall，1996．

［11］SAMBRIDGE M，MOSEGARRD K．Monte Carlo methods in geophysical inverse problems［J］．Review of Geophysics，2002，40，doi：10．1029／2000RG000089．

［12］MALINVERNO A，BRIGGS V A．Expanded uncertainty quantification in inverse problems：hierarchical Bayes and empirical Bayes［J］．Geophysics，2004，69（4）：1005．

［13］TARITS P，MENVIELLE J V，ROUSSIGNOL M．Bayesian statistics of non－linear inverse problems：example of the magnetotelluric 1－D inverse problem［J］．Geophysical Journal International，1994，119（2）：353．

［14］GRANDIS H，MENVIELLE M，ROUSSIGNOL M．Bayesian inversion with Markov chains－I．The magnetotelluricone－dimensional case［J］．Geophysical Journal International，1999，138（3）：757．

［15］CERV V，MENVIELLE M，PEK J，Stochastic interpretation of magnetotelluric data，comparison of methods［J］．Annals of Geophysics，2007，50（1）：122．

［16］SAMBRIDGE M．Geophysical inversion with a neighbourhood algorithm—i．Searching a parameter space［J］．Geophysical Journal International，1999，138（2）：479．

［17］GREEN P J．Reversible jump Markov chain Monte Carlo computation and Bayesian model determination［J］．Biometrika，1995，82（4）：711．

［18］GUO RONGWEN，DOSSO STAN E，LIU JIANXIU，et al．Non－linearity in Bayesian 1－D magnetotelluric inversion［J］．Geophysical Journal International，2011，185（2）：663．

［19］MITSUHATA Y，UCHIDA T．2－D Inversion of frequency－domain EM data caused by a 3－D source［M］．in Three－Dimensional Electromagnetics，153～172，eds Zhdanov M S，Wannamaker P E，New York：Elsevier，2002．

［20］PRESS W H，TEUKOLSKY S A，VETTERLING W T，et al．Numerical Recipes in FORTRAN［M］．2nd edn．Cambridge：Cambridge Univ．Press，1992．

［21］METROPOLIS N，ROSENBLUTH A W，ROSENBLUTH M N，et al．Equation of State Calculations by Fast Computing Machines［J］．Journal of Chemical Physics，1953，21（6）：1087．

［22］HASTINGS W K．Monte Carlo sampling methods using Markov chains and their applications［J］．Biometrika，1970（57）：97．

［23］DOSSO S E，WILMUT M J．Uncertainty estimation in simultaneous Bayesian tracking and environmental inversion［J］．Journal of the Acoustical Society of A－merica，2008（124）：82．

［24］JONES A G，HUTTON R A，multi－station magnetotelluric study in southern Scotland–I．Fieldwork，data analysis and results［J］．Geophysical Journal of the Royal Astronomical Society，1979，56（2）：329．

［25］JONES A G，HUTTON R A．multi－station magnetotelluric study in southern Scotland–II．Monte－Carloinversion of the data and its geophysical and tectonic implications［J］．Geophysical Journal of the Royal Astronomical Society，1979，56（2）：351．

［26］CONSTABLE S C，PARKER R L，CONSTABLE C G．Occam′s inversion：A practical algorithm for generating smooth models from electromagnetic sounding data［J］．Geophysics，1987，52（3）：289．

P 631．3＋25

A

10．3969／j．issn．1001－1749．2012．04．01

國家科技支撐計劃項目（2011BAB04B08）；中國地質調查局科研項目（資［2011］03－01－64）；有色資源與地質災害探查湖南省重點實驗室項目（2010TP4012－6）

2011－09－30改回日期：2011－12－25

1001—1749（2012）04—0371—09

史學東（1962－），男，碩士，高級工程師，現主要從事物探工作。