频率域航空电磁法已被广泛应用于地质填图、资源勘探及环境工程等领域(雷栋等,2006).电阻率成像技术以其快速高效的特点在航空电磁数据处理与解释中发挥着重要作用(Sengpiel,1988;Huang and Fraser, 1996).然而,为了获得更加准确、定量化的地下电性分布信息,需要对观测数据进行反演处理.传统的梯度反演方法首先将非线性反演问题线性化,然后沿着目标函数最速下降方向搜索极小点,以期通过迭代技术获得最优解(Egbert and Kelbert, 2012).此类方法具有较高的工作效率,但是解的收敛性受初始模型影响很大,尤其是针对目标函数为多峰值的复杂情况,反演解容易陷入局部极小,严重影响了反演的精确性和稳定性,给进行正确的地质解释造成困难.
近年来,完全非线性反演方法在地球物理数据解释中越来越受到重视.师学明等(2000)提出大地电磁数据的多尺度逐次逼近遗传反演算法.通过模拟生物进化过程,采用随机优化技术,基于选择、变换和变异三种基本操作来处理模型参数,并以较大的概率搜索全局最优解.采用多尺度逐次逼近遗传算法,解决了有效基因丢失和早熟收敛问题.徐海浪和吴小平(2006)实现了电阻率二维神经网络反演.该方法是一种模拟人脑处理信息功能的新型人工智能技术,由神经元、网络结构和学习规则三个基本要素构成,利用大量简单的、高度互连的处理元素组成复杂的网络计算系统.利用神经网络可以解决许多传统人工智能方法无法解决的难题(姚姚,2002).Yin和Hodges(2007)将模拟退火反演方法应用到航空电磁数据处理中,通过模拟流体冷却结晶的热力学过程,按照Boltzmann分布沿上行和下行双向搜索全局最小值,获得准确的反演结果.
然而,上述反演方法都只能提供单一的“最佳”反演结果,并不能获得反演模型参数的不确定度信息.考虑到反演问题解的非唯一性,传统方法无法对反演结果的可靠性进行评价.为了解决这一问题,贝叶斯反演(Bayesian Inversion)方法受到越来越多的关注.经验贝叶斯方法和分级贝叶斯方法已经被应用于一维大地电磁数据(Guo et al,2011)和2.5维接地长导线源频率域电磁数据(Mitsuhata,2002)的反演当中,通过加入模型约束,采用高效的优化方法,加快收敛速度,改善了反演效果.变维数贝叶斯反演方法是一种新发展起来的概率反演方法,利用可逆跳跃马尔科夫链蒙特卡洛方法(RJMCMC),按照建议分布在模型空间进行采样,根据接受概率接受或拒绝采样模型,最终获得反演参数的概率分布.此方法受初始模型影响小,收敛稳定,可以提供更可靠的反演结果.相对前两种贝叶斯反演,在变维数贝叶斯反演方法中,模型层数作为反演参数也参与反演计算,其固有的吝啬性使得在满足数据拟合条件下使用最少的层数,获得最简单的反演模型.
针对变维数贝叶斯反演方法,国外的地球物理学家已经进行了卓有成效的研究.Green(1995)对传统的MCMC方法进行改进,提出了模型参数个数可变的RJMCMC方法,并与贝叶斯理论相结合,为变维数贝叶斯反演的发展奠定了理论基础.Malinverno(2002)首次将RJMCMC方法应用于直流电阻率一维层状介质反演中,为RJMCMC方法在电磁数据反演领域中的应用开辟了先河.Sambridge等(2006)从贝叶斯理论角度对变维数反演问题进行了深入研究和细致讨论.Agostinetti 等(2010)采用变维数蒙特卡洛采样方法对远震接收函数数据进行处理,获得地下介质的弹性性质.Minsley(2011)利用变维数的贝叶斯方法进行频率域航空电磁反演,通过施加纵向光滑约束有效地改善了反演效果.Ray和Key(2012)将变维数贝叶斯反演方法应用于海洋电磁各向异性数据反演解释中.
相比之下,国内对于变维数贝叶斯方法在电磁数据反演领域的研究相对较少.本文对变维数贝叶斯方法进行改进并将其应用于航空电磁数据反演解释中.针对传统的变维数贝叶斯反演方法对深部良导层反演效果不佳的问题,通过合理引入先验信息加权系数,调整反演模型的约束强度,在很大程度上改善了反演效果.除此之外,通过改进模型统计方法,只将满足数据拟合要求的模型纳入统计范围,削弱了不合理模型对统计结果的干扰.我们首先在文中简要介绍一维频率域航空电磁正演理论,然后对变维数贝叶斯算法的基本原理做详细阐述,最后通过对包含高斯噪声的理论数据和实测数据进行反演,并与Occam反演结果进行对比分析,验证通过调整先验信息权重改善反演结果的有效性.
2 航空电磁1D正演理论目前,频率域航空电磁系统主要采用水平共面(HCP)装置.根据Yin和Hodges(2007),垂直磁偶极发射源对应的磁场垂直分量为

,h0为发射源高度,z为接收点纵坐标,r为发射源到接收机的水平距离,m是发射线圈的偶极矩,T(z-)是汉克尔积分,由下式给出,

反演处理中所使用的数据为去除一次场并做归一化的数据:

根据贝叶斯公式,后验概率可以由先验概率和似然函数给出:



如果数据噪声服从高斯分布,似然函数定义为多维正态分布,

根据Minsley(2011),层数和发射源高度的先验概率分别定义为(kmin,kmax)和(htxmin,htxmax)上的均匀分布,一般kmin取为1,kmax的取值要足够大,以满足数据拟合的要求.相应的界面深度先验概率分布可表示为

电阻率的先验概率分布定义为多维正态分布,

在贝叶斯反演中,模型约束信息包含在先验分布中,然而在以往的变维数贝叶斯反演中没有充分考虑改变先验分布权重对反演结果的影响.针对埋深较大的低阻层反演效果较差,不能准确反映真实电阻率和界面位置这一问题,本文在先验概率分布中引入权系数,通过合理选择权系数,改变先验分布的权重,达到调整模型约束强度,改善反演效果的目的.方程(9)中的p(ρ |k)可以重新定义为

在进行贝叶斯反演之前首先要对模型参数进行初始化,给定反演模型参数的变化范围:最小和最大 层数(kmin,kmax)、层界面最小和最大深度(zmin,zmax)、 最小和最大发射源高度(htxmin,htxmax)以及电阻率变化范围(ρmin,ρmax).初始模型一般选为均匀半空间,层数为两层,层界面位于(zmax-zmin)/2+zmin处,初始发射源高度为实际测量高度.为了确保反演 参数为正,同时各反演变量可以自由变化,我们对各模型参数取对数,在对数域内进行反演计算(Yin,2000).
3.3 建议分布在变维数贝叶斯反演中,候选模型根据建议分布函数定义,且只与当前模型有关,而与之前的反演结果无关,所以贝叶斯反演受初始模型的影响小.建议分布可以表示为

理论上,建议分布的选取对反演结果没有影响,然而事实证明合理选取建议分布对提高工作效率有很大帮助.理想的建议分布应与后验分布相同,这样 可以以较大的概率接受候选模型(Malinverno,2002).
根据Green(1995)和Minsley(2011),候选模型使用RJMCMC方法采样获得.层数和层界面位置由新层生成—旧层灭亡过程定义,此过程包括四种状态:(1)新层生成(p=1/6):在确保大于最小层厚的条件下,一个新的层界面在zmin和zmax之间随机产生,层数加1;(2)旧层灭亡(p=1/6):随机选择一个层界面并将其删除,层数减1;(3)扰动更新(p=1/6):层数保持不变,随机选择一个层界面,其位置在(-hmin,hmin)范围内扰动;(4)停滞不变(p=1/2):层数和界面位置都保持不变.
电阻率的建议分布定义为以当前模型电阻率为均值的多维正态分布,

为后验电阻率方差,

为权系数.后验电阻率方差的选取十分重要,它决定采样步长的大小,影响采样效率.如果采样步长太大,会造成接受概率过低,大量样本被拒绝;如果采样步长太小,会造成接受概率过高,收敛速度下降(Ray and Key, 2012).合理施加放 缩系数fac,可以调节搜索步长,进而提高反演效率,

候选发射源高度服从以当前发射源高度为均值的一维正态分布(Minsley,2011):

根据建议分布产生候选模型后,按照接受概率判断新模型是否被接受,接受概率由下式给出(Green,1995):

,一般取为1(Agostinetti and Malinverno, 2010).
根据(4)—(15)式,考虑到均匀分布和对称分布的特点,可以得到简化后的接受概率表达式

在频率域航空电磁法中,模型参数与电磁响应之间存在十分复杂的非线性关系,目标函数可能出现多个峰值.单独使用一条马尔科夫链往往不能对整个模型空间进行充分搜索.为解决这一问题,本文使用3条马尔科夫链,采用不同的初始模型和随机种子同时对反演模型进行搜索.使用OpenMP并行技术,计算效率得到大幅提升.
在以往的贝叶斯反演中,当数据拟合差小于事先设定的阀值时,“预热”阶段(burn-in period)结束,将之后所有被接受的采样模型输出并计入统计.然而,为了能对模型空间进行充分搜索,贝叶斯反演在“预热”阶段结束后,仍接受一些不满足数据拟合要求的采样模型以跳出局部极小.如果将这些样本纳入统计范围,将会在一定程度上影响统计结果.本文在遵循原有的模型采样方法和接受标准的基础上,只将满足数据拟合要求的模型纳入统计范围,削弱不合理模型对统计结果的干扰.
本文采用Guitton和Hoversten(2011)所使用的势尺度衰减因子(potential scale reduction factor-PSRF)作为判断收敛的条件.PSRF综合考查每条链内部以及链与链之间采样样本的统计分布,为判断变维数贝叶斯反演收敛提供可靠依据.对于简单模型情况,当PSRF<1.1时,判断反演收敛.然而对于复杂模型,应适当放宽收敛条件,根据Chen(2008),定义PSRF<1.2为判断收敛条件.
4 理论数据反演为了验证通过调整先验信息权重改善反演结果的有效性,本文设计了两个水平层状模型,针对HCP装置计算5个频率(386、1538、6257、25790 Hz和100264 Hz)的垂直磁场,并在理论数据中加入3%的高斯随机噪声,利用变维数贝叶斯反演方法对其进行处理.发射源到接收机的距离为10 m,发射源距地面高度为30 m.
4.1 三层模型反演本文首先设计了一个三层介质模型(图 1),第一层是厚度为30 m、电阻率为300 Ωm的高阻层;第二层是厚度为5 m、电阻率为10 Ωm的低阻薄层;第三层是电阻率为100 Ωm的高阻基底.在垂直磁场的理论模拟数据中加入3%的随机高斯噪声(如图 2所示).最大层数设为30层;最小界面深度 为1 m,最大界面深度为120 m;最小电阻率为0.1 Ωm,最大电阻率为10000 Ωm.下面分两种情况对变维数贝叶斯反演效果进行讨论:第一种情况假设发射源高度保持不变(图 3);第二种情况假设发射源高度参与反演,可在20 m到40 m之间变化(图 4).
![]() | 图 1 理论模型Fig. 1 Synthetic model |
![]() | 图 2 反演数据Fig. 2 Inversion data |
![]() | 图 3 发射源高度固定不同权系数反演结果对比图:(A1)—(A3)为传统变维数贝叶斯反演结果,(B1)—(B3)为λ=30的反演结果,(C1)—(C3)为λ=100的反演结果Fig. 3 Comparison of inversion results for different weights with fixed transmitter altitude:(A1)—(A3)for traditional trans-dimensional Bayesian inversion,(B1)—(B3)for inversions with λ=30,while(C1)—(C3)for inversions with λ=100 |
![]() | 图 4 发射源高度参与反演时不同权系数反演结果对比图:(A1)—(A4)为λ=1的反演结果,而(B1)—(B4)为λ=100的反演结果Fig. 4 Comparison of inversion results for different weights with variable transmitter altitude:(A1)—(A4)for inversions with λ=1,while(B1)—(B4)for inversions with λ=100 |
图 3为发射源高度固定,权系数不同的变维数贝叶斯反演结果.其中(A1)—(A3)为利用传统变维数贝叶斯方法获得的反演结果,(B1)—(B3)和(C1)—(C3)分别为λ等于30和100时的反演结果.(A1)、(B1)和(C1)是RJMCMC方法采样模型的复合分布,灰黑色阴影表示反演模型的概率分布,颜色越深概率越大,阴影部分的面积在一定程度上反映了反演结果的不确定性,不确定性与阴影面积成正比.图中黑色实线代表真实模型,黑色虚线代表反演结果中概率最大模型,灰色虚线为Occam反演结果,灰色实线是针对不含噪声数据的Occam反演结果.(A2)、(B2)、(C2)和(A3)、(B3)、(C3)分别是界面深度和层数的概率统计结果,图中黑色实线指示了真实界面位置和真实层数.三组反演结果均显示地下模型分为3层,且第二层低阻层,但是传统变维数贝叶斯反演方法对于低阻薄层的反演效果不佳,电阻率和界面位置都不准确,而且灰色阴影面积较大,表明反演结果的不确定性较强.相比之下,B组和C组在引入权系数λ后,反演效果得到很大改善.界面位置与实际模型吻合较好,电阻率更接近真实值.对比B组和C组的结果,可发现随着λ增大,模型约束强度也随着增强,阴影面积也相应减小,采样模型更加集中,反演结果对低阻薄层的反映更加准确,可靠性提高.对比Occam方法对不含噪声和含有高斯噪声两组数据的反演结果,可以看出噪声对Occam反演结果影响较大,相比之下,贝叶斯反演结果更加稳定.除此之外,与Occam反演结果相比,变维数贝叶斯反演结果可以划分出清晰的层界面,不产生强烈震荡引起的假异常.
4.1.2 发射源高度可变的情况在实际工作中,由于树木和植被等原因,航空电磁系统难以获得精确的发射源高度.为解决这一问题,本文将发射源高度作为一个变量参与反演计算.图 4为λ=1和λ=100的两组反演结果的对比.
图 4中,(A4)和(B4)是发射源高度的统计分布,其余图示与图 3中意义相同.从图中可以看出,由于发射源高度不确定,使得变维数贝叶斯反演变得复杂,反演效果受到很大影响.当λ=1时,模型层数、界面位置和发射源高度均不准确,与真实模型差距较大.对比A、B两组反演结果,可以发现通过增大权系数加强模型约束,反演效果得到很大改善.模型层数和发射源高度分布较为集中,与真实的模型参数吻合很好,界面位置基本准确.与(A1)相比,图(B1)中采样模型分布更加集中,阴影面积较小,没有出现大面积低概率渲染区,反演结果的可靠性得到很大提高.
4.2 四层模型反演为了进一步验证通过调整先验信息权重改善反演结果的有效性,本文设计了更复杂的四层介质模型.第一层是厚度为10 m、电阻率为100 Ωm的高阻层;第二层是厚度为5 m、电阻率为10 Ωm的低阻薄层;第三层是厚度为50 m、电阻率为300 Ωm的高阻层;第四层是电阻率为50 Ωm的均匀半空间.我们在理论垂直磁场数据中加入了3%的随机高斯噪声.假设最大层数为30层;最小界面深度为 1 m,最大界面深度为120 m;最小电阻率为0.1 Ωm,最大电阻率为10000 Ωm; 发射源高度参与反演,在20 m到40 m之间变化.
从图 5中可以看出,对于复杂的四层模型,λ=1时反演效果并不理想,发射源高度和界面位置反演都不准确,最大概率模型与真实模型之间差距很大.大面积低阻渲染区的出现说明反演结果的不确定性较强.相比之下,当λ调整为70时,反演效果得到明显改善.发射源高度与实际高度一致,最大概率模型和界面位置与真实模型基本吻合.采样模型以较大概率集中分布在真实模型附近,表明反演模型的可靠性得到很大的提高.
![]() | 图 5 发射源高度参与反演时不同权系数反演结果对比图:(A1)—(A4)为λ=1的反演结果,而(B1)—(B4)为λ=70的反演结果Fig. 5 Comparison of inversion results for different weights with variable transmitter altitude:(A1)—(A4)for inversions with λ=1,while(B1)—(B4)for inversions with λ=70 |
本文利用改进变维数贝叶斯反演方法对某地近地表勘查的实测数据进行反演,验证此方法的有效性.实测数据采用HCP装置测量垂直磁场分量,观测频率为300 Hz~100 kHz.
图 6为Occam反演结果,初始模型为10 Ωm的均匀半空间.由图可以看出,整条剖面的层连续性较好,反映地下的电性结构大体分为3层:第一层电阻率为8~18 Ωm,层厚约为10 m的高阻层,经实际勘查,确定为泥沙层;第二层电阻为3~8 Ωm,层厚约为10 m的低阻层,确定为粘土层;第三层电阻率为8~18 Ωm的高阻基底,确定为含砂质地层.
![]() | 图 6 航空电磁实测数据Occam反演结果Fig. 6 Inversion results of survey data with Occam′s method |
本文在上述剖面中选取了具有代表性的A、B、C、D四个测点进行变维数贝叶斯反演.根据实际资料,假设最大层数为10层;最小界面深度为1 m,最大界面深度为40 m;最小电阻率为1 Ωm,最大电阻 率为100 Ωm;发射源高度参与反演,在20 m到40 m 之间变化.初始模型为30 Ωm的均匀半空间,针对四个测点取权系数λ分别为30,20,15和30.
图 7是A、B、C和D四个测点的贝叶斯反演结果,其中(A1—A4)—(D1—D4)分别是测点A—D的电阻率—深度、界面位置、层数和发射源高度的反演结果.在(A1—D1)中黑色实线为Occam的反演结果,而黑色虚线为最大概率模型.可以看出,两种方法反演结果的整体趋势基本吻合,相比之下,贝叶斯反演结果能划分出清晰的层界面,不产生强烈的震荡.通过施加合适的加权因子,采样模型的分布较为集中,反演模型的可靠性较高.在(A4—D4)中黑色实线表示发射源高度的实测值,可以发现反演的发射源高度与实测高度基本相同,只是略有差别.
![]() | 图 7 航空电磁实测数据变维数贝叶斯反演结果:(A1)—(A4),(B1)—(B4),(C1)—(C4),(D1)—(D4)分别为测点A、B、C和D的反演结果Fig. 7 Inversion results of survey data with trans-dimensional Bayesian method:(A1)—(A4),(B1)—(B4),(C1)—(C4),(D1)—(D4)are results for station A,B,C and D respectively |
鉴于传统梯度反演方法受初始模型影响大,容易陷入局部极小的缺点,本文采用一种全新的完全非线性反演方法—变维数贝叶斯反演.通过对模型空间进行全局搜索,有效克服了初始模型对反演解的影响.根据建议分布采样获得候选模型,并根据接受概率筛选合理的采样模型并计入统计,最终得到反演模型的概率分布和模型参数的不确定度信息,为判断反演结果的可靠性提供依据.针对传统变维数贝叶斯方法不能准确反映深部低阻层的电性分布信息,本文在先验信息中引入权系数,通过改变权系数调整对反演模型的约束强度,有效地改善了反演效果.在此基础上,对模型统计方法进行改进,避免了不合理模型的干扰.最后对含有高斯噪声的理论数据和实测航空电磁数据进行反演处理,验证了通过调整权系数来改善反演结果的有效性,同时有力地证实本文提出的加权约束变维数贝叶斯反演是一种航空电磁数据处理的有效反演手段.
致谢 作者向贲放和黄威博士在文章准备过程中提供的帮助表示感谢.我们要特别感谢各位审稿人对本文提出的宝贵意见.| [1] | Agostinetti N P, Malinverno A. 2010. Receiver function inversion by trans-dimensional Monte Carlo sampling. Geophys. J. Int., 181(2): 858-872, doi:10.1111/j.1365-246X.2010. 04530.x. |
| [2] | Chen J S, Kemna A, Hubbard 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. |
| [3] | Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems. Geophys. J. Int., 189(1): 251-267, doi:10.1111/j.1365-246X.2011.05347.x. |
| [4] | Green P J. 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82(4): 711-732. |
| [5] | Guo R W, Dosso S E, Liu J X, et al. 2011. Non-linearity in Bayesian 1-D magnetotelluric inversion. Geophys. J. Int., 185(2): 663-675, doi:10.1111/j.1365-246X.2011.04996.x. |
| [6] | Huang H P, Fraser D C. 1996. The differential parameter method for multifrequency airborne resistivity mapping. Geophysics, 61(1): 100-109. |
| [7] | Lei D, Hu X Y, Zhang S F. 2006. Development status of Airborne Electromagnetic. Contributions to Geology and Mineral Resources Research, 21(1):40-53. |
| [8] | Malinverno A. 2002. Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem. Geophys. J. Int., 151(3): 675-688. |
| [9] | Minsley B J. 2011. A trans-dimensional Bayesian Markov chain Monte Carlo algorithm for model assessment using frequency-domain electromagnetic data. Geophys. J. Int., 187(1): 252-272. |
| [10] | Mitsuhata Y, Uchida T, Amano H. 2002. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source. Geophysics, 67(6): 1753-1768. |
| [11] | Ray A, Key K. 2012. Bayesian inversion of marine CSEM data with a trans-dimensional self parametrizing algorithm. Geophys. J. Int., 191(3): 1135-1151, doi:10.1111/j.1365-246X. 2012.05677.x. |
| [12] | Sambridge M, Gallagher K, Jackson A, et al. 2006. Trans-dimensional inverse problems, model comparison and the evidence. Geophys. J. Int., 167(2): 528-542, doi:10.1111/j.1365-246X.2006.03155.x. |
| [13] | Sengpiel K P. 1988. Approximate inversion of airborne EM data from a multilayered ground. Geophysical Prospecting, 36(4): 446-459. |
| [14] | Shi X M, Wang J Y, Zhang S Y. 2000. Multiscale genetic algorithm and its application in magnetotelluric sounding data inversion. Chinese J. Geophys. (in Chinese), 43(1): 122-130. |
| [15] | Trainor-Guitton W, Hoversten G M. 2011. Stochastic inversion for electromagnetic geophysics: Practical challenges and improving convergence efficiency. Geophysics, 76(6): F373-F386, doi:10.1190/GEO2010-0223.1. |
| [16] | Xu H L, Wu X P. 2006. 2D resistivity inversion using the neural network method. Chinese J. Geophys. (in Chinese), 49(2): 584-589. |
| [17] | Yao Y. 2002. The Basic Theory and Application of Geophysical Inversion Methods. Beijing: China University of Geosciences Press. |
| [18] | Yin C C. 2000. Geoelectrical inversion for a one-dimensional anisotropic model and inherent non-uniqueness. Geophys. J. Int., 140: 11-23. |
| [19] | Yin C C, Hodges G. 2007. Simulated annealing for airborne EM inversion. Geophysics, 72(4): F189-F196. |
| [20] | 雷栋, 胡祥云, 张素芳. 2006. 航空电磁法的发展现状. 地质找矿论丛, 21(1): 40-53. |
| [21] | 师学明, 王家映, 张胜业等. 2000. 多尺度逐次逼近遗传算法反演大地电磁资料. 地球物理学报, 43(1): 122-130. |
| [22] | 徐海浪, 吴小平. 2006. 电阻率二维神经网络反演. 地球物理学报, 49(2): 584-589. |
| [23] | 姚姚. 2002. 地球物理反演基本理论与应用方法. 北京: 中国地质大学出版社. |
2014, Vol. 57









