2018, Vol. 61
大地电磁测深法(Magnetotelluric, MT)通过地表电场和磁场观测，获得较宽频带天然电磁场响应数据，反演获得地球内部电性结构分布(Tikhonov, 1950; Cagniard, 1953; 陈乐寿和王光锷, 1990).过去几十年，MT反演一直以二维单点估计为主(De Groot-Hedlin and Constable, 1990; Siripunvaraporn and Egbert, 2000, 2007; Rodi and Mackie, 2001; 胡祖志等, 2005; 刘小军等, 2007; 蔡军涛等，2010; 韩波等, 2012; 叶涛等，2013; 陈小斌等，2014；陈小斌和郭春玲，2017；康敏等, 2017).但是，由于观测数据和模型化过程存在误差(Tarantola, 2005)，在可接受误差范围内的反演结果可能很多，单点估计较难评估反演模型解的不确定度，概率反演能弥补这种缺陷(Afonso et al., 2013).
大地电磁反演通常假设连续场分布，对介质分界面不太敏感.通行的基于Tikhonov正则化理论进行的反演加入光滑或平滑条件进行约束，使边界区分难度更大.因此常规的反演方式不太适用于尖锐边界的反演.常见的尖锐边界反演有两种，一种由Zhdanov等(Portniaguine and Zhdanov, 1999; Mehanee and Zhdanov, 2002; Auken and Christiansen, 2004; Zhang et al., 2010; 刘小军等, 2007; 张罗磊等, 2009)提出，基于最小梯度支撑(Minimum Gradient Support, MGS)方法进行反演.但是，通常情况下该种方式对异常块体的大小有一定压缩.另一种尖锐边界反演是针对二维层状结构较好的模型，通过加入横向平滑或光滑及层间厚度差的约束，提高层状性(Smith et al., 1999; De Groot-Hedlin and Constable, 2004; Chen et al., 2012；欧东新和王家林，2005; 叶益信等，2009；柴伦炜等, 2016)，该方式已成功用于层状结构较好的贝叶斯反演(Chen et al., 2012).但对于垂直结构较明显或者异常块体反演的结构，该方法作用有限.结合已知单点估计的手段，以区块的方式，直接对边界和各块体内的电阻率进行反演，能更好地分辨边界信息.
以贝叶斯理论为基础的尖锐边界反演，在反演出边界及电阻率的同时，直接对边界及电阻率的不确定度进行评价.前人通过贝叶斯理论对一维大地电磁反演问题进行了研究，说明了贝叶斯理论用于大地电磁反演的可行性(Tarits et al., 1994; Grandis et al., 1999; Cerv et al., 2007; Guo et al., 2011, 2014; Mandolesi et al., 2018; Xiang et al., 2018; 师学明和王家映, 1998; 杨迪琨和胡祥云, 2008; 郭荣文, 2011; 史学冬等, 2012；孙欢乐等，2016).受制于正演速度和采样量，高维的大地电磁贝叶斯反演研究相对较少，目前已知的只有Chen等(2012)进行了二维尖锐边界的大地电磁贝叶斯反演的尝试，针对层状结构较好的模型以电阻率和各层深度为模型参数实施反演，并将反演结果和三维单点估计结果对比，说明了方法的有效性.但是，该方法不适用于单独异常体或者垂直结构明显的模型.
本文基于贝叶斯理论开展二维MT数据的贝叶斯反演，计算尖锐边界的位置和各区域的电阻率，采用Markov Chain Monte Carlo(MCMC)方法进行采样.反演以Kelbert等(2014)的ModEM反演结果为基础，分析地下电阻率的大致分布，以边界位置及内部电阻率为模型参数对二维大地电磁数据进行反演.
1 方法介绍 1.1 贝叶斯理论贝叶斯理论能将所有先验及模型信息包含在后验概率密度分布(Posteriori Probability Density, PPD)中.关于公式的分析，可参考Tarantola(2005)、Ray和Key(2012)和史学冬等(2012).
如果数据误差服从高斯分布，我们可以得到后验概率密度分布σM(m)如下:
|
(1) |
上式中dobs为观测数据，g(m)为正演过程，C为误差协方差矩阵，ρM(m)为模型的先验概率分布，k为归一化常数，
分析σM(m)可获取模型相应信息，比如最大后验概率估计模型mbest，概率期望模型E(m)，模型协方差矩阵Cov(m)，边缘概率分布σM(mi)及相关度矩阵R(m)等:
|
(2) |
|
(3) |
|
(4) |
|
(5) |
|
(6) |
其中δ为Dirac delta函数，T为转置符号.若无先验信息或先验信息为高斯线性信息，且似然函数为高斯分布，对于线性反演问题，反演结果服从高斯分布，上式有解析解.大地电磁反演问题具有较强的非线性，上式无解析解，需要采用数值计算方法计算以上积分.本文选择MCMC方法对σM(m)分布进行采样，以数值方法计算上述积分.
1.2 正演正演是反演的核心，要求高效、可信、准确.本次正演部分基于目前国内外广泛采用的ModEM程序(Kelbert et al., 2014)，具体实现过程参阅Siripunvaraporn和Egbert(2000).
1.3 Metropolis-Hastings采样方法贝叶斯反演主要任务是计算1.1节中的积分，常用MCMC方法的Metropolis-Hastings(MH)采样对后验概率密度函数σM(m)进行采样，然后进行数值积分.MH方法由Metropolis等(1953)首次提出，经Hastings(1970)改进推广.根据建议分布函数p(xi|xj)选取候选模型，以一定概率接受或拒绝候选模型的方式完成采样.p(xi|xj)表示当前模型为xj时扰动生成候选模型xi的概率，π(x)为目标函数，按以下概率接受候选模型：
|
(7) |
xj扰动生成候选模型xi时，按照aj的概率接受候选模型：
|
(8) |
拒绝候选模型时，则当前模型将保存到采样样本(链)中.如果接受候选模型，则将候选模型保存到采样链中.
采样收敛时，建议分布函数不影响反演结果，但影响采样时的收敛效率.构造建议分布函数的方法很多(Tierney, 1994; Chib and Greenberg, 1995)，一种是Metropolis等(1953)提出的p(xi|xj)=p1(xj-xi).一种是Hastings(1970)提出的p(xi|xj)=p2(xi)，候选模型的选取与当前模型无关.一种是当π(x)∝φ(x)θ(x)时(φ(x)和θ(x)定义域需相同(Chib and Greenberg, 1994))，令p(xi|xj)=θ(xi)，此时对于ai其取值如下：
|
(9) |
以后验概率采样为例，若先验概率分布较易单独采样，则相对采样效率通常会高很多.还有自回归(Autoregressive)和拒绝采样(Rejection sampling)等构造方法，具体可参考Tierney(1994)、Chib和Greenberg(1995).接受率(Acceptance Rate)是评判建议分布函数的重要指标，接受率过高或过低都会导致采样效率低下，导致采样样本的混合性差.一般接受率建议为30%~50%(Chib and Greenberg, 1995; Tarantola, 2005; Bodin and Sambridge, 2009; Ray and Key, 2012).
1.4 反演在没有其他先验信息和其他约束时(为表示方便称为无先验信息)，采用的先验信息为均匀分布；存在其他先验信息时，将先验信息包含在后验概率密度函数中.本文目标函数(1)不加任何平滑或光滑类约束，在传统反演中约束本质是使结果更平稳、光滑，对于病态方程防止模型的小误差导致数据的巨大变化，即模型中可能存在一个较大的极值区域使得通常的NLCG等梯度方法无法越过，严重影响反演结果的有效性.概率采样求解过程，某种意义上是全局搜索，数据变化并不直接影响模型空间的跳动，而是影响跳动后是否接受新的跳动点.不加约束，可以更直接查看模型的敏感性，若某些参数分布集中，则该变量对数据的敏感性较大.
本文在ModEM二维反演的基础上，选取异常最明显的位置，设置相应模型情况.对于块体部分，设置上下位置深度为随机变量，实际中由于使用有限差分进行正演，设置块体数来反映更方便{h1, 1，h2, 1, …, hn-1, 1，hn, 1，h1, 2，h2, 2, …, hn-1, 2，hn, 2}，其中hi, 1表示i块体上边界所在层数，hi, 2表示i块体下边界所在层数，如图 1所示，同一块体中电阻率相同，不同块体设置不同值.若无先验信息，建议分布函数具有对称性(即p(xi|xj)=p(xj|xi))，则具体采样过程如下：
|
图 1 网格划分示意图 Fig. 1 Grid partition of anomalous body |
(1) 设初始模型为m0；
(2) 设总随机变量为N，当前模型为mi，生成1到N均匀分布的随机整数，对当前模型mi中该序号的模型参数按建议分布函数进行扰动，得到候选模型mnext；
(3) 计算mnext对应的数据dnext=g(mnext)，计算拟合残差Fnext=‖dobs-dnext‖2如下：
|
其中C为数据误差的关联矩阵，一般认为数据误差相互独立，且方差已知；
(4) 若Fnext < Fi，
则mi+1=mnext；
否则生成[0, 1]间均匀分布的随机数μ:
若μ < exp(-(Fnext-Fi)/2)
则mi+1=mnext
否则mi+1=mi；
(5) 按照(2), (3), (4)步重复计算直到达到采样数为止.
若有先验信息，步骤(3)将比较广义残差，对包含先验信息的σM(m)进行比较.当建议分布函数不对称时，按前述经Hastings(1970)修正的公式，将不对称的建议分布函数加入考虑.
本文用于模拟的模型如图 2所示.正演使用ModEM中内嵌的二维代码，计算过程是有限差分实现的.在计算区域内，横向为1500m每格，纵向每一层深度分别为[10, 30, 60, 100, 200, 300, 500, 800, 1000, 1200, 1500, 1800, 2100, 2500, 3000, 3600, 4300, 5000, 6500, 8000, 10000, 10000, 15000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 80000](m)，共31层.异常体所处网格共有40列，模型异常体底面深度在15层网格，深度为15.1km.电阻率共三组，分别记为ρi(i=1, 2, 3)，取自然对数分别为2.30259，6.90776，4.60517 Ωm.假设数据记录台站设置在-45.00 km到52.50 km之间的地表，每2.5 km取一个点，共40个采样点，频率为0.01 Hz~100 s按等比取12份，数据格式为阻抗，对所有数据点加5%标准差的高斯误差.根据趋肤深度计算公式
|
图 2 真实模型电阻率示意图，红色部分为低阻10 Ωm，蓝色部分为高阻1000 Ωm，背景为100 Ωm Fig. 2 True resistivity model. The resistivity of red part is 10 Ωm, blue part is high resistivity 1000 Ωm and background is 100 Ωm |
反演结果如图 3所示，从反演结果可以看出，该方法能对左右边界有较好的辨识度，但底边界的分辨能力明显不足，故将底边界和两部分异常体及背景电阻率作为模型参数.两部分异常块体大致分布在横向在-30~30 km的范围内，而网格划分中横向每格为1500 m，设置40个模型参数记录底边界位置，分别记为hi(i=1, 2, …, 40)，为了计算方便，hi均表示层数，取值范围为0~31的整数.将图 3中红色部分电阻率记为ρ1，蓝色为ρ2和背景为ρ3，由于正演代码中以电阻率的自然对数进行存储，直接以ln(ρ1)、ln(ρ2)和ln(ρ3)为模型参数，取值范围为(0, 10)(Ωm).故模型参数为hi(i=1, 2, 3, …, 40)和ln(ρi)(i=1, 2, 3)，共有43个模型参数，记为m=(h1, h2, …, h39, h40, ln(ρ1), ln(ρ2), ln(ρ3))，参数设置示意图如图 4所示.
|
图 3 尖锐边界反演与ModEM反演对比.背景为ModEM反演结果，红线为hi采样频次最高的值对应的深度，灰线为以此为基准扩展到60%样本点的位置的上下限位置，黑色虚线为实际块体的边界，白线为边界的概率期望模型解 Fig. 3 Sharp boundary inversion vs. ModEM inversion result. Background is the ModEM inversion result, red line is the maximum frequency in sampling of hi, which has been transferred to depth values, grey line is the upper and lower limit of boundary for 60% sample points, the dotted black line is true boundary, the white line is expectation of boundary |
|
图 4 模型参数设置示意图，其中两块体的电阻率为ρ1和ρ2，背景为ρ3，hi(i=1, 2, …, 40)表示边界所处层数 Fig. 4 Model variable, where ρ1 and ρ2 denote resistivity of two anomalous blocks, ρ3 denotes resistivity of background where hi(i=1, 2, …, 40) is the boundary, which counts on layer |
本文中，当前模型为mn时，根据建议分布函数选取候选模型mnext的方式如下.
生成1到43均匀分布的随机整数k：
a) 若k≤ 40，mn中的参数mn, k按如下建议分布函数选取mnext, k，其余参数不变，得到候选模型mnext：
|
(10) |
b) 若k>40，mn中的参数mn, k按如下建议分布函数选取mnext, k，其余参数不变，得到候选模型mnext：
|
(11) |
其中A~N(0.00, 0.01)，由于高斯分布的对称性，该点与所选点互选为候选模型的概率相同.采样过程按照前述采样过程所述步骤重复进行.采样的起始位置设在13层网格处，深度为9.6 km，(ln(ρ1), ln(ρ2), ln(ρ3))分别从3.3，3.5，5.5 Ωm开始，共进行了83万次采样.
2 反演结果归一化均方根(normalized Root Mean Square, nRMS)是跟踪采样过程、分析采样效率的重要手段，其计算方式为
|
(12) |
其中N为数据点数.此次采样的nRMS变化如图 5所示，nRMS在前三百次中快速下降，再迭代700次后达到稳定状态，可视为达到平稳分布，此时采样为有效采样，若无专门说明，本文模型分析以有效采样的样本完成.以真实模型带入算出nRMS为1.031，而迭代中存在最小的nRMS为1.0276，说明真实模型并不一定为最优模型.整体接受率为42.22%，满足30%~50%的接受率(Tarantola, 2005; Chib and Greenberg, 1995).结合nRMS分布和接受率判定此次采样是有效的.
|
图 5 拟合nRMS随迭代次数的变化 (a)表示所有迭代; (b)前300次迭代; (c) 301~1000次迭代; (d) 1000次迭代以后的情况. Fig. 5 nRMS changes with iteration times, (a) is for all, (b) is for the first 300 iterations, (c) is for 301th to 1000th, while (d) is the latter iterations |
根据式(5)算出hi的边缘概率分布，部分结果如图 6，电阻率的边缘概率分布如图 7所示.从图 6、图 7可以看出整体的分布是集中的，由于先验分布是均匀分布，不含任何约束，这种集中完全由数据拟合导致，采样中有部分层数达到31，说明采样是比较完全的.对于mbest的求解，ln(ρi)(i=1, 2, 3)是连续随机变量，直接从(2)式计算很困难，需要与传统的梯度方式相结合，如Gauss-Newton及其衍生算法，在集中分布附近设置初始位置即可求得全局最优解(Chen et al., 2008).
|
图 6 部分hi有效采样的边缘概率分布 hi表示采样下边界位置所处网格层数，真实模型的下边界均位于15层，图中横轴表示层数，纵轴表示频次，(a)—(d)分别表示h8，h16，h24，h32的边缘概率分布. Fig. 6 Marginal probability distribution of different sample layers hi is the layer of lower boundary, which is 15 for true, horizontal axis is the layer, while vertical axis is frequency.(a)—(d) is the samples of h8, h16, h24, h32. |
|
图 7 (a)、(b)、(c)分别表示ρ1、ρ2、ρ3有效采样的边缘概率分布情况， 其中灰线表示真值，横轴为电阻率取自然对数的坐标 Fig. 7 Marginal probability distribution of ρ1, ρ2 and ρ3. Grey line is the true value and horizontal axis is the logarithm of resistivity ρ |
hi出现最多频次及其占比情况见附录表 1，其中占比最少为29.75%，最多为64.14%，除h1, h2和h21出现频次最多的值为14，其余均为15，对边界刻画较好.按照(4)式计算模型协方差矩阵，再通过(6)式计算各模型参数之间的相关度，模型参数之间的相关度如图 8所示，从数据相关系数分布可以看出，模型参数之间的互相关很弱，故以边缘概率分布计算各自的置信区间是合理的.按照(3)式计算模型的概率期望模型，误差棒以均值为基准扩散至总采样点的60%，结果如图 9所示，左侧低阻部分边界的敏感性相对于右侧高阻部分明显会弱一些.电阻率的计算类似，对于ln(ρ1)，分布于区间[2.293, 2.313] (Ωm)的点占比为94.97%，对于ln(ρ3)，分布于[4.591, 4.611](Ωm)的点占比为96.15%，对于ln(ρ2)，分布于区间[6.898, 6.928] (Ωm)的点占比为96.12%.这种占比对应置信区间的范围，范围越小说明集中度越高，敏感度越高.尖锐边界反演结果与ModEM反演结果对比见图 3，选用边缘概率中采样频次最大的值会比均值的结果好，尖锐边界反演对边界约束较好.
|
图 8 模型参数互相关矩阵Ri, j，坐标轴中的1到40分别表示h1~h40，41~43表示ln(ρ1)、ln(ρ2)、ln(ρ3) Fig. 8 Model variables′ correlation coefficient Ri, j, in axis, 1~40 means h1~h40, 41~43 means ln (ρ1)、ln (ρ2)、ln (ρ3) |
|
图 9 下边界hi(i=1, 2, …, 40)所处层数及其误差，误差棒计算从中值处扩展至总采样点60% Fig. 9 Mean depth and error bar of hi, error bar counts to 60 percent |
部分数据拟合情况如图 10所示.层数以最大分布点、电阻率以均值进行计算，得到数据的拟合情况，考虑60%模型分布下对应的数据置信区间.数据拟合分为TE、TM和实部、虚部共四种情况展示，拟合结果较好.由于正演过程的非线性，极值点分布未知.重复在模型参数的置信区间内进行随机采样，计算每次采样对应的正演数据，记录各个数据点的最大值和最小值，以此为数据置信区间.
|
图 10 含置信区间的数据拟合，数据置信区间通过模型置信区间获得，以24号台站为例 Fig. 10 Curve fitting with confidence interval (e.g. station 24) |
设置不同起始点，若能在达到平稳过程后汇聚于同一区域，说明此处为全局最优，采样可信.本研究设置三组对照组，分别采样2万次，ln(ρ1)、ln(ρ2)、ln(ρ3)起始位置均设置为(3.3000, 5.5000, 3.5000)(Ωm).边界起始位置设置分别为：对照组1，设置左半侧hi从11层开始，右半侧hi从19层开始；对照组2，设置左半侧hi从19层开始，右半侧hi从11层开始；对照组3全部hi设置为从19层开始.图 11到图 13分别显示了对照组1到组3的nRMS变化过程.从图 11—13中可以看出，由于距离高概率区域选取较远，采样在1000次甚至更往后一段才能达到一个相对平稳的区域.在采样结束时采样点的位置见附录表 2，此时采样在真值附近进行.
|
图 11 对照组1的2万次采样nRMS变化细节 Fig. 11 nRMS change details of compared group 1 during 20000 sampling |
|
图 12 对照组2的2万次采样nRMS变化细节 Fig. 12 nRMS change details of compared group 2 during 20000 sampling |
|
图 13 对照组3的2万次采样nRMS变化细节 Fig. 13 nRMS change details of compared group 3 during 20, 000 sampling |
建议分布函数的选择对采样影响极大.建议分布函数从接受率和接受区域两个方面影响着Markov链的表现，其中接受率主要受到建议分布函数的步长设计影响，接受区域主要受到建议分布函数的整体表现影响.若选择建议分布函数不好，接受率可能极低，比如按Hastings(1970)提出的方式，选择全空间内均匀分布作为建议分布函数，候选模型直接在全空间随机选取.采样过程中的nRMS(图 14)初始下降速度较快，但是后段采样效率极低，整体接受概率只有0.84%.这是由于采样游走步调过大，导致真正高概率区域内的点较少被选中.实际中，高概率区域是最值得关注的部分，该区域对结果的解释更为重要.采样点有限时，更需要采样对高概率区域刻画充分.
|
图 14 建议分布函数为均匀分布且直接在全空间随机选取候选模型的nRMS.(a)所有nRMS变化，(b)1000次采样之后的变化 Fig. 14 nRMS changes of homogenous distribution as the proposal distribution function, while sampling randomly in whole model variable space. (a) is for all, while (b) is the nRMS change after 1000 times iteration |
上述采样不依赖于先验信息，从结果看已经能较好地描述真实模型，若有一定的先验信息，比如Chen等(2012)中提及的边界上下限或者边界平滑等约束，该过程或许能更好地对边界位置进行约束.作为对照组，加入hi上下限的约束，上限11层下限20层，设限后的方程可参考Chen等(2012)，此处将11与20层强行设置为相邻区域.新的定义域内目标函数同比增大，表现为归一化因子k变化.采样起始位置为左半边hi均为12层，右半边hi均为18层，电阻率取自然对数为(3.3000, 5.5000, 3.5000)(Ωm)，采样数为20万.此次采样过程中nRMS的表现与前述模型类似，故取1000次采样后的点为样本进行分析，采样后边界位置均值和60%的采样点分布如图 15所示.图 15与图 9相似度很高，主要因为加入的约束极弱，两相对比能合理推测h1、h2偏小及h5~h10偏大是由数据自身约束有限导致的.图 16为电阻率分布情况，从图中可以看出电阻率分布集中性较好.ln(ρ1)、ln(ρ2)、ln(ρ3)的计算均值分别为(2.3041，6.9101，4.5998)(Ωm)，其中ln(ρ1)取值在[2.204, 2.404] (Ωm)占比100.00%，ln(ρ2)取值在[6.8100, 7.0100] (Ωm)占比99.92%，ln(ρ3)取值在[4.591，4.611] (Ωm)占比98.60%.
|
图 15 限制边界反演结果，误差棒计算从中值处扩展至总采样点60% Fig. 15 Constrained boundary inversion, error bar counts to 60 percent |
|
图 16 (a)、(b)、(c)分别表示限制边界反演下ρ1、ρ2、ρ3的边缘概率分布情况，其中灰线表示真值，横轴为电阻率取自然对数的坐标，(a)为所有采样点计算结果，(b)为其中99.87%的采样点计算结果，(c)为其中98.60%的采样点计算结果 Fig. 16 Marginal distribution of ρ1、ρ2、ρ3 in constrained boundary inversion. Grey line is the true value and horizontal axis is the natural logarithm of resistivity. (a) has all sample points of ρ1, (b) has 99.87% sample points of ρ2, while (c) has 98.60% sample points of ρ3 |
模型结构的对比，以左侧为15层厚，右侧为14层厚构建新模型，ln(ρ1)、ln(ρ2)、ln(ρ3)为(2.30259，6.90776，4.60517)(Ωm)，合成数据及台站设置同前.采样的起始位置为左半边hi均为12层，右半边hi均为18层，电阻率取自然对数为(3.3000, 5.5000, 3.5000)(Ωm)，采样数为20万.以1000次采样后的点为样本进行分析，边界结果如图 17，当边界改变后合成数据用于反演时，反演的边界依然能做出类似的改变.图 18为电阻率分布情况，从图中可以看出电阻率分布集中性较好.ln(ρ1)、ln(ρ2)、ln(ρ3)的均值为(2.3040, 6.9156, 4.6000)(Ωm)，ln(ρ1)取值在[2.204, 2.404] (Ωm)占比100.00%，ln(ρ2)取值在[6.816, 7.016](Ωm)占比100.00%，ln(ρ3)取值在[4.500, 4.700](Ωm)占比99.98%.对不同模型该方法依然是有效的.
|
图 17 左侧15层，右侧14层边界反演结果，误差棒计算从中值处扩展至总采样点60% Fig. 17 Boundary inversion of left 15 layers and right 14 layers, error bar counts to 60 percent |
|
图 18 (a)、(b)、(c)分别表示左侧15、右侧14层的模型结构反演结果的ρ1、ρ2、ρ3的边缘概率分布情况，其中灰线表示真值，横轴为电阻率取自然对数的坐标，(a)为所有采样点计算结果，(b)为其中99.80%的采样点计算结果，(c)为其中99.21%的采样点计算结果 Fig. 18 Marginal distribution of ρ1、ρ2、ρ3 in left 15 layers and right 14 layers of new model. Grey line is true value and horizontal axis is the natural logarithm of resistivity. (a) has all sample points of ρ1; (b) has 99.80% sample points of ρ2; (c) has 99.21% sample points of ρ3 |
合成数据反演结果说明贝叶斯理论应用于二维大地电磁尖锐边界反演是可行的，即使边界对数据不敏感，依然能通过大量采样得出其分布.相较于传统反演方法，贝叶斯反演有三点优势：(一)反演结果不依赖于初始模型，初始模型仅仅影响模型收敛的速度.(二)给出了解的评价机制，即模型不确定度的估计.单点估计通常不会提供模型的不确定信息，常规的尖锐边界反演和贝叶斯反演均能提供模型的不确定信息，但提供方式和适用情况不同.常规尖锐边界反演在最优解处以其Jacobian矩阵构造协方差矩阵来度量不确定度(Auken and Christiansen, 2004)，适用于模型非线性程度较弱且解为全局最优.但是，对于非线性较强的MT反演上述方式不可行(Tarantola, 2005).贝叶斯反演通过分析采样结果便能获得模型参数的置信区间，对反演结果进行客观评价.(三)灵活性，由于模型化过程和观测数据存在误差，可能存在一定范围内的模型均可较好地满足数据拟合.贝叶斯反演提供的模型置信区间恰好是对这种范围的描述.当与传统结果的比较或结合时，可以使用最大分布点、中点、均值点等特殊值.
使用贝叶斯反演需要注意：首先，在多极值点存在的情况下，虽然贝叶斯反演能有效避免陷入局部极小值，但是，当局部极值点相对周边的概率较大，在采样不足的情况下，仍可能陷入局部极值(Vrugt, 2016).为了避免陷入局部极值，通常选用不同初始值产生Markov链，以稳定收敛后的散点图分布、比率诊断法、Teweke Test法等(朱新玲, 2009)对结果准确性进行评估.其次，关于采样效率，直接关联的是接受率，通常建议30%~50%的接受率(Tarantola, 2005; Chib and Greenberg, 1995)，图 14及多次模拟发现，对于概率相差过大的区域，在有限的计算次数内，不合适的采样方式不能有效刻画目标函数的分布，导致采样失败.通常解决方案为限制采样步长和选取合适的建议分布函数，有时直接限制边界上下限也会有效(Chen et al., 2012).但是，提供的上下限约束不强时，如文中所加约束，对采样结果几乎没有任何改善.此外，由于贝叶斯反演需要大量的计算资源和存储空间，计算成本是该方法在实际应用中需要面临的一个挑战.对于较为简单的基于MCMC采样方式进行计算的公式，通常至少需要数十万至数百万次的正演才能有足够的采样对结果进行解释.而对于复杂的模型结构，通常需要较多的随机变量，为了获得有效的模型信息，计算时间和存储空间需求将显著增加，因此，该方法的推广应用还有待高效并行算法的进一步研发以及计算资源的不断完善.
附录
|
附表 1 最多频次值、均值 TableS 1 The maximum frequency in sampling of hi and mean value of all variables |
|
|
附表 2 不同起始点采样2万次结束时所处位置 TableS 2 Last value after 20000 iterations for different starting value |
Afonso J C, Fullea J, Griffin W L, et al. 2013. 3-D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle. I:A priori petrological information and geophysical observables. Journal of Geophysical Research:Solid Earth, 118(5): 2586-2617. DOI:10.1002/jgrb.50124 |
Auken E, Christiansen A V. 2004. Layered and laterally constrained 2D inversion of resistivity data. Geophysics, 69(10): 752-761. DOI:10.1190/1.1759461 |
Bodin T, Sambridge M. 2009. Seismic tomography with the reversible jump algorithm. Geophysical Journal International, 178(3): 1411-1436. DOI:10.1111/j.1365-246x.2009.04226.x |
Cagniard L. 1953. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics, 18(3): 605-635. DOI:10.1190/1.1437915 |
Cai J T, Chen X B. 2010. Refined techniques for data processing and two-dimensional inversion in magnetotelluric Ⅱ:Which data polarization mode should be used in 2D inversion. Chinese Journal of Geophysics (in Chinese), 53(11): 2703-2714. DOI:10.3969/j.issn.0001-5733.2010.11.018 |
Cerv V, Menvielle M, Pek J. 2007. Stochastic interpretation of magnetotelluric data, comparison of methods. Annals of Geophysics, 50(1): 7-19. |
Chai L W, Li T L, Su X B, et al. 2016. Sharp boundary inversion used in Linjiang Liudaogou controlled source audio magnetotelluric sounding data processing. Progress in Geophysics, 31(3): 943-948. DOI:10.6038/pg20160302 |
Chen J S, Kemna A, Hubbard S S. 2008. A comparison between gauss-newton and Markov-chain Monte Carlo-based methods for inverting spectral induced-polarization data for Cole-Cole parameters. Geophysics, 73(6): F247-F259. DOI:10.1190/1.2976115 |
Chen J S, Hoversten G M, Key K, et al. 2012. Stochastic inversion of magnetotelluric data using a sharp boundary parameterization and application to a geothermal site. Geophysics, 77(4): E265-E279. DOI:10.1190/geo2011-0430.1 |
Chen L S, Wang G E. 1990. Magnetotelluric Sounding Method (in Chinese). Beijing: Geological Publishing House.
|
Chen X B, Cai J T, Wang L F, et al. 2014. Refined techniques for magnetotelluric data processing and two-dimensional inversion (Ⅳ):Statistical image method based on multi-site, multi-frequency tensor decomposition. Chinese Journal of Geophysics (in Chinese), 57(6): 1946-1957. DOI:10.6038/cjg20140625 |
Chen X B, Guo C L. 2017. Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅴ):Detecting the linear structures of the Earth by impedance tensor imaging. Chinese Journal of Geophysics (in Chinese), 60(2): 766-777. DOI:10.6038/cjg20170227 |
Chib S, Greenberg E. 1994. Bayes inference in regression models with ARMA (p, q) errors. Journal of Econometrics, 64(1-2): 183-206. DOI:10.1016/0304-4076(94)90063-9 |
Chib S, Greenberg E. 1995. Understanding the metropolis-Hastings algorithm. The American Statistician, 49(4): 327-335. DOI:10.2307/2684568 |
De Groot-Hedlin C, Constable S. 1990. Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics, 55(55): 1613-1624. DOI:10.1190/1.1442813 |
De Groot-Hedlin C, Constable S. 2004. Inversion of magnetotelluric data for 2D structure with sharp resistivity-contrast. Geophysics, 69(1): 78-86. DOI:10.1190/1.1649377 |
Grandis H, Menvielle M, Roussignol M. 1999. Bayesian inversion with Markov chains-I. The magnetotelluric one-dimensional case. Geophysical Journal of the Royal Astronomical Society, 138(3): 757-768. DOI:10.1046/j.1365-246x.1999.00904.x |
Guo R W. 2011. Nonlinearity and uncertainty analysis based on Bayeisan MT inversion[Ph. D. thesis] (in Chinese). Changsha: Central South University.
|
Guo R W, Dosso S E, Liu J X, et al. 2011. Non-linearity in Bayesian 1-D magnetotelluric inversion. Geophysical Journal International, 185(2): 663-675. DOI:10.1111/j.1365-246x.2011.04996.x |
Guo R W, Dosso S E, Liu J X, et al. 2014. Frequency-and spatial-correlated noise on layered magnetotelluric inversion. Geophysical Journal International, 199(2): 1205-1213. DOI:10.1093/gji/ggu329 |
Han B, Hu X Y, He Z X, et al. 2012. Mathematical classification of magnetotelluric inversion methods. Oil Geophysical Prospecting (in Chinese), 47(1): 177-188. |
Hastings W K. 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1): 97-109. DOI:10.1093/biomet/57.1.97 |
Hu Z Z, Hu X Y, Wu W L, et al. 2005. Compared study of two-dimensional magnetotelluric inversion methods. Coal Geology & Exploration (in Chinese), 33(1): 64-68. |
Kang M, Hu X Y, Kang J, et al. 2017. Compared of magnetotelluric 2D inversion methods. Progress in Geophysics (in Chinese), 32(2): 476-486. DOI:10.6038/pg20170205 |
Kelbert A, Meqbel N, Egbert G D, et al. 2014. ModEM:A modular system for inversion of electromagnetic geophysical data. Computers & Geosciences, 66: 40-53. DOI:10.1016/j.cageo.2014.01.010 |
Liu X J, Wang J L, Chen B, et al. 2007. Discussion on focus inversion algorithm of 2-D MT data. Oil Geophysical Prospecting (in Chinese), 42(3): 338-342. |
Mandolesi E, Ogaya X, Campanyà J, et al. 2018. A reversible-jump Markov chain Monte Carlo algorithm for 1D inversion of magnetotelluric data. Computers & Geosciences, 113: 94-105. |
Mehanee S, Zhdanov M. 2002. Two-dimensional magnetotelluric inversion of blocky geoelectrical structures. Journal of Geophysical Research:Solid Earth, 107(B4): EPM-1-EPM 2-11. DOI:10.1029/2001jb000191 |
Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. 1953. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6): 1087-1092. DOI:10.1063/1.1699114 |
Ou D X, Wang J L. 2005. The fast magnetotelluric inversion of 2-D block structure. Geophysical Prospecting for Petroleum (in Chinese), 44(5): 525-528. |
Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596 |
Ray A, Key K. 2012. Bayesian inversion of marine CSEM data with a trans-dimensional self parametrizing algorithm. Geophysical Journal International, 191(3): 1135-1151. DOI:10.1111/j.1365-246x.2012.05677.x |
Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics, 66(1): 174-187. DOI:10.1190/1.1444893 |
Shi X D, Liu J X, Guo R W, et al. 2012. Unbiased 1D magnetotelluric Bayesian inversion. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 34(4): 371-379. DOI:10.3969/j.issn.1001-1749.2012.04.01 |
Shi X M, Wang J Y. 1998. One dimensional magnetotelluric sounding inversion using simulated annealing. Earth Science-Journal of China University of Geosciences (in Chinese), 23(5): 542-546. |
Siripunvaraporn W, Egbert G. 2000. An efficient data-subspace inversion method for 2-D magnetotelluric data. Geophysics, 65(3): 791-803. DOI:10.1190/1.1444778 |
Siripunvaraporn W, Egbert G. 2007. Data space conjugate gradient inversion for 2-D magnetotelluric data. Geophysical Journal International, 170(3): 986-994. DOI:10.1111/j.1365-246x.2007.03478.x |
Smith T, Hoversten M, Gasperikova E, et al. 1999. Sharp boundary inversion of 2D magnetotelluric data. Geophysical Prospecting, 47(4): 469-486. DOI:10.1046/j.1365-2478.1999.00145.x |
Sun H L, Wang S B, Guo R W, et al. 2016. One dimensional magnetotelluric sounding apparent resistivity and phase inversion of adaptive simplex simulated annealing study. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 38(5): 584-592. |
Tarantola A. 2005. Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics, Philadelphia, ISBN 0-89871-572-5.
|
Tarits P, Jouanne V, Menvielle M, et al. 1994. Bayesian statistics of non-linear inverse problems:Example of the magnetotelluric 1-D inverse problem. Geophysical Journal International, 119(2): 353-368. DOI:10.1111/j.1365-246x.1994.tb00128.x |
Tierney L. 1994. Markov chains for exploring posterior distributions. The Annals of Statistics, 22(4): 1701-1762. DOI:10.1214/aos/1176325750 |
Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of the Earth's crust. Doklady Akademii Nauk SSSR, 73(2): 295-297. |
Vrugt J A. 2016. Markov chain Monte Carlo simulation using the DREAM software package:Theory, concepts, and MATLAB implementation. Environmental Modelling & Software, 75: 273-316. DOI:10.1016/j.envsoft.2015.08.013 |
Xiang E M, Guo R W, Dosso S E, et al. 2018. Efficient hierarchical trans-dimensional Bayesian inversion of magnetotelluric data. Geophysical Journal International, 213(3): 1751-1767. DOI:10.1093/gji/ggy071 |
Yang D K, Hu X Y. 2008. Inversion of noisy data by probabilistic methodology. Chinese Journal of Geophysics (in Chinese), 51(3): 901-907. |
Ye T, Chen X B, Yan L J. 2013. Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅲ):Using the impressing method to construct starting model of 2D magnetotelluric inversion. Chinese Journal of Geophysics (in Chinese), 56(10): 3596-3606. DOI:10.6038/cjg20131034 |
Ye Y X, Hu X Y, Jin G X, et al. 2009. Application analysis of sharp boundary inversion of magnetotelluric data for 2D structure. Progress in Geophysics (in Chinese), 24(2): 668-674. DOI:10.3969/j.issn.1004-2903.2009.02.041 |
Zhang L L, Yu P, Wang J L, et al. 2009. Smoothest model and sharp boundary based two-dimensional magnetotelluric inversion. Chinese Journal of Geophysics (in Chinese), 52(6): 1625-1632. DOI:10.3969/j.issn.0001-5733.2009.06.025 |
Zhang L L, Yu P, Wang J L, et al. 2010. A study on 2D magnetotelluric sharp boundary inversion. Chinese Journal of Geophysics, 53(3): 631-637. DOI:10.3969/j.issn.0001-5733.2010.03.017 |
Zhu X L. 2009. Research summary on MCMC methods. Statistics & Decision, (21): 151-153. |
蔡军涛, 陈小斌. 2010. 大地电磁资料精细处理和二维反演解释技术研究(二)——反演数据极化模式选择. 地球物理学报, 53(11): 2703-2714. DOI:10.3969/j.issn.0001-5733.2010.11.018 |
柴伦炜, 李桐林, 苏晓波, 等. 2016. 尖锐边界反演在临江市六道沟可控源音频大地电磁测深数据反演中的应用. 地球物理学进展, 31(3): 943-948. DOI:10.6038/pg20160302 |
陈乐寿, 王光锷. 1990. 大地电磁测深法. 北京: 地质出版社.
|
陈小斌, 蔡军涛, 王立凤, 等. 2014. 大地电磁资料精细处理和二维反演解释技术研究(四)——阻抗张量分解的多测点-多频点统计成像分析. 地球物理学报, 57(6): 1946-1957. DOI:10.6038/cjg20140625 |
陈小斌, 郭春玲. 2017. 大地电磁资料精细处理和二维反演解释技术研究(五)——利用阻抗张量成像识别大地线性构造. 地球物理学报, 60(2): 766-777. DOI:10.6038/cjg20170227 |
郭荣文. 2011. 贝叶斯MT反演的非线性和不确定度分析[博士论文]. 长沙: 中南大学.
|
韩波, 胡祥云, 何展翔, 等. 2012. 大地电磁反演方法的数学分类. 石油地球物理勘探, 47(1): 177-188. |
胡祖志, 胡祥云, 吴文鹂, 等. 2005. 大地电磁二维反演方法对比研究. 煤田地质与勘探, 33(1): 64-68. |
康敏, 胡祥云, 康健, 等. 2017. 大地电磁二维反演方法分析对比. 地球物理学进展, 32(2): 476-486. DOI:10.6038/pg20170205 |
刘小军, 王家林, 陈冰, 等. 2007. 二维大地电磁数据的聚焦反演算法探讨. 石油地球物理勘探, 42(3): 338-342. |
欧东新, 王家林. 2005. 二维块状结构大地电磁快速反演. 石油物探, 44(5): 525-528. |
史学东, 柳建新, 郭荣文, 等. 2012. 大地电磁法的1D无偏差贝叶斯反演. 物探化探计算技术, 34(4): 371-379. DOI:10.3969/j.issn.1001-1749.2012.04.01 |
师学明, 王家映. 1998. 一维层状介质大地电磁模拟退火反演法. 地球科学-中国地质大学学报, 23(5): 542-546. |
孙欢乐, 王世彪, 郭荣文, 等. 2016. 基于自适应纯形模拟退火法一维大地电磁测深视电阻率和相位反演研究. 物探化探计算技术, 38(5): 584-592. |
杨迪琨, 胡祥云. 2008. 含噪声数据反演的概率描述. 地球物理学报, 51(3): 901-907. |
叶涛, 陈小斌, 严良俊. 2013. 大地电磁资料精细处理和二维反演解释技术研究(三)——构建二维反演初始模型的印模法. 地球物理学报, 56(10): 3596-3606. DOI:10.6038/cjg20131034 |
叶益信, 胡祥云, 金钢燮, 等. 2009. 大地电磁二维陡边界反演应用效果分析. 地球物理学进展, 24(2): 668-674. DOI:10.3969/j.issn.1004-2903.2009.02.041 |
张罗磊, 于鹏, 王家林, 等. 2009. 光滑模型与尖锐边界结合的MT二维反演方法. 地球物理学报, 52(6): 1625-1632. DOI:10.3969/j.issn.0001-5733.2009.06.025 |
朱新玲. 2009. 马尔科夫链蒙特卡罗方法研究综述. 统计与决策, (21): 151-153. |