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

史學(xué)東，柳建新，郭榮文，童孝忠，劉文劼，曹創(chuàng)華

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

史學(xué)東1，柳建新2，郭榮文2，童孝忠2，劉文劼2，曹創(chuàng)華2

（1．中原油田對(duì)外經(jīng)濟(jì)貿(mào)易總公司，濮陽(yáng) 457001；2．中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083）

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

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

0 前言

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

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

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

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

1 貝葉斯理論

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

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

其中P（m｜d）表示后驗(yàn)概率密度函數(shù)。

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

其中 積分區(qū)域v表示多維模型空間；φ（m）表示包括了數(shù)據(jù)擬合和先驗(yàn)信息的廣義數(shù)據(jù)擬合函數(shù)，其表達(dá)式見(jiàn)式（4）。

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

其中δ表示Dirac delta函數(shù)；T表示轉(zhuǎn)置符號(hào)。

要提取以上信息，需要用優(yōu)化方法求解好的初始模型，并對(duì)整個(gè)模型區(qū)域進(jìn)行積分。對(duì)于線性反演問(wèn)題，如果已知數(shù)據(jù)誤差服從Gaussian分布，且先驗(yàn)信息服從均勻概率分布或Gaussian分布，那么PPD本身也將服從Gaussian分布，以上屬性問(wèn)題存在解析解。對(duì)于非線性反演問(wèn)題，可以采用線性化過(guò)程和線性反演方法來(lái)求解，其近似程度取決于反演問(wèn)題的非線性程度。另一種替代方法就是通過(guò)數(shù)值方法求解，它沒(méi)有線性化近似誤差，具有普遍性，但是需要更高強(qiáng)度的計(jì)算量為代價(jià)。

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

2 Metroplis－Hastings采樣數(shù)值積分

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

其中f（m）表示一般函數(shù)。

式（11）定義的積分可以通過(guò)解析方法或數(shù)值方法求解，作者主要考慮采用基于Monte Carlo的積分求解法。在標(biāo)準(zhǔn)Monte Carlo積分中［20］，模型是隨機(jī)從定義在積分區(qū)域上的均勻分布中產(chǎn)生。給定一個(gè)含Q個(gè)模型的樣本，每個(gè)樣本的概率可通過(guò)式（12）計(jì)算。

其中i＝1，…，Q

積分式（11）通過(guò)式（13）進(jìn)行估計(jì)。

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

如果積分式集中在模型空間的某些局部區(qū)域，那么從均勻分布上進(jìn)行采樣，其效率不高?？梢圆捎弥攸c(diǎn)采樣解決均勻采樣效率低的問(wèn)題，該方法對(duì)積分貢獻(xiàn)較大的區(qū)域進(jìn)行重點(diǎn)采樣。假設(shè)g（m）代表重點(diǎn)采樣的采樣分布函數(shù)，Q個(gè)模型從該分布上產(chǎn)生，其歸一化條件見(jiàn)式（14）。

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

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

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

因此，通過(guò)MHS采樣，式（11）的積分簡(jiǎn)化成在樣本集上計(jì)算f（m）的平均值。

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

3 主分量空間的擾動(dòng)

如果模型空間的相關(guān)性比較強(qiáng)，那么基于算法（Markov Chain Monte Carlo）的優(yōu)化方法和采樣方法的效率將受到嚴(yán)重的影響。在這種情況下，可以將模型空間轉(zhuǎn)換到主分量空間中以減輕這種影響。物理空間與主分量參數(shù)空間的轉(zhuǎn)換，可通過(guò)式（17）實(shí)現(xiàn)。

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

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

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

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

在“burn－in”采樣過(guò)程中，協(xié)方差矩陣Cm的起始線性化估計(jì)自適應(yīng)地被基于MHS采樣的非線性估計(jì)代替，它可以更好地代表參數(shù)空間的總體協(xié)方差（這個(gè)階段的采樣沒(méi)有保存用于積分估計(jì)）。

4 采樣溫度

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

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

假設(shè)一個(gè)由Q個(gè)模型組成的樣本，服從exp［－φ（m）／T］定義的分布（等價(jià)于在溫度“T”時(shí)的MHS采樣），那么式（21）可變?yōu)槭剑?2）。

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

5 反演結(jié)果分析

5．1 理論模型反演

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

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

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

5．2 實(shí)測(cè)數(shù)據(jù)反演

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

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

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

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

6 結(jié)論

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

作者在文中分別對(duì)未知數(shù)據(jù)誤差條件下的五層模型及實(shí)測(cè)數(shù)據(jù)進(jìn)行了無(wú)偏差不確定度分析。從無(wú)偏差貝葉斯反演的邊緣概率分布，得出各參數(shù)不確定度和參數(shù)間的相關(guān)關(guān)系，以及不同深度模型的不確定度分布。不同溫度下的COPROD1數(shù)據(jù)反演結(jié)果表明，在T＝1時(shí)采樣徹底，反演結(jié)果得到雙峰結(jié)構(gòu)。非地電參數(shù)當(dāng)成隨機(jī)變量在反演中得到了約束，說(shuō)明數(shù)據(jù)本身不僅包含地球物理模型信息（電導(dǎo)率），還包含了這些非地電參數(shù)的信息。

［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

國(guó)家科技支撐計(jì)劃項(xiàng)目（2011BAB04B08）；中國(guó)地質(zhì)調(diào)查局科研項(xiàng)目（資［2011］03－01－64）；有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室項(xiàng)目（2010TP4012－6）

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

1001—1749（2012）04—0371—09

史學(xué)東（1962－），男，碩士，高級(jí)工程師，現(xiàn)主要從事物探工作。