大地电磁法的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－），男，硕士，高级工程师，现主要从事物探工作。