2. 中国地质大学地球物理与信息技术学院, 北京 100083
2. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
IMPULSE系统是中国国土资源航空物探遥感中心在2002年由加拿大引进的频率域吊舱式直升机航空电磁测量系统.该设备拥有水平共面与垂直同轴两组线圈,每组线圈可以同时发射三个不同频率的电磁信号,具有工作频率多、信息量大和分辨率高等技术特点[1].目前,该系统已经完成了多个地区矿产勘查与科学研究任务,提高了国内航空电磁应用水平,取得了初步的地质效果.目前,国际上有关航空电磁法一维正反演以及快速成像方面的研究工作已比较成熟,普遍使用一维两层以上和近似计算的解释软件.由于国内整体研究水平不高,还未对航空电磁数据处理解释方法开展过系统的研究工作,一定程度上限制了航空电磁资料的解释效果.伴随IMPULSE系统引进的解释软件功能较为单一,仅能开展航空视电阻率计算及一维近似反演工作,计算精度低,应用效果有限.相对于地面电磁测量而言,频率域航空电磁测量在工作方式、收录信息(地面电磁测量通常收录总磁场值,航空电磁测量仅收录归一化的二次磁场值)、有效工作频率数(航空电磁测量每种装置仅有三个工作频率)等方面具有明显的差异性,需要在理论计算公式、初始模型赋值、反演控制参数选取等方面进行进一步的分析与研究,形成实用化的算法,实现频率域吊舱式直升机航空电磁资料的马奎特反演计算,为地质解释人员提供一种有效的解释手段.
2 基本原理与方法 2.1 基本原理马奎特反演也叫阻尼最小二乘反演,是一种基于最小二乘准则的反演方法,它是通过最小化仅包含数据残差的一个目标函数来求得模型“最优解”的反演方法.马奎特反演过程只追求模拟数据与原始测量数据的最大拟合,因而具有算法简单、计算速度快的特点.为了方便对马奎特反演原理进行说明,首先需要对一维正演方程进行简化表示.
对于水平层状介质上空,频率域吊舱式直升机航空电磁测量系统(图 1)一维正演计算方程如下:
(1) |
垂直同轴装置
(2) |
(3) |
(4) |
(5) |
(6) |
(7) |
其中,H2/H1为接收线圈处二次磁场强度与一次磁场强度的比值,r为收发线圈距,h为飞行高度,RTE表示TE激化模式下的反射系数,J0()、J1()分别为零阶和一阶贝塞尔函数,μ0为自由空间磁导率,即μ0=4π×10-7H/m,ω为发射线圈工作圆频率,tl是第l层介质厚度,σl为第l层介质电导率.方程(1)、(2)可以简化为
(8) |
d=(d1,d2,…,dN)为数据向量,m=(m1,m2,…,mM)为模型参数向量(电导率、层厚度),F为一维正演方程.
于是,马奎特反演目标函数定义为[5]
(9) |
其中,σ1,σ2,…,σN为各观测数据对应的噪声均方差.
由于频率域航空电磁法一维正演方程F是高度非线性的方程,需要对其在初始模型m0处展开为泰勒展开式,并且忽略二次以上导数项,得到差分等式:
(10) |
其中Δm为模型增量,Δd为观测数据与模型响应差值,即Δd=d-F(m0),J为N×M阶雅可比矩阵,其元素如下:
(11) |
对于方程(10),当J为非奇异矩阵时,有唯一解:
(12) |
当J为奇异或接近奇异矩阵时,用(12)式求解误差会很大,可能导致反演过程不收敛.据此,需要对对称矩阵JTJ正交分解,把JTJ分解为RΛRT,这里的Λ是JTJ的特征值构成的对角线矩阵,而R是JTJ的特征向量矩阵且满足
(13) |
为了改善JTJ的求逆条件,对其变换如下:
(14) |
(15) |
λi是对称矩阵JTJ的第i个特征值,α为一个大于零的实数.此时矩阵的条件数为[6]
(16) |
这就是著名的马奎特法又称阻尼最小二乘法,此时,模型改变量为
(17) |
其中,α为阻尼系数或叫加权因子,I为N×N阶单位矩阵.对大多数地球物理反演而言,α为零或足够大都难以取得模型的最优解.这里不存在一种简单的计算阻尼系数的方法,具体算法编程中,阻尼系数是通过“尝试法”在迭代过程中不断加以修改[5],从而得到一个合适的α值.事实上,马奎特反演也正是在对目标函数Φ全面最优解的搜索过程中,通过调整α的大小达到调整了搜索方向和步长的目的.
通过以上方法,经过k次迭代反演后得到的新模型可以表示为
(18) |
通过将新模型代入目标函数(9)求取新的拟合差值Φk,并与期望拟合差值Φ*相比较,当拟合差达到期望拟合差值(Φk≤Φ*)或者迭代次数达到设定最大迭代次数(通常为5~10)时迭代过程停止,这时得到的模型mk(电阻率、层厚度)即为马奎特反演结果模型.
从理论上讲,马奎特反演底层模型为均匀半空间,没有深度上的限制,该方法的实际反演深度取决于系统探测深度和测量数据的精度.
2.2 雅可比矩阵计算雅可比矩阵是地球物理反演中极为重要的一个关键环节,其计算精度高低将直接影响到反演的精度和反演结果的正确性.雅可比矩阵的计算一方面要保持尽可能高的精度,另一方面要满足较快的速度,雅可比矩阵的求取通常有以下两种方法:
2.2.1 差分法差分法是求一阶导数最常用的方法.该方法求一个偏导数需要两次正演计算,对于计算参数相对较少,正演速度快的一维算法来说,该方法具有编程简单、计算速度快的特点.缺点是计算精度较低,在精度要求不太高时选用差分法求取雅可比矩阵会是一种很好的选择.从一维层状反演应用效果看,模型改变量通常选取1%基本上就可以满足反演计算精度的要求.
2.2.2 解析法由于一维正演计算公式相对简单,可以通过解析法直接求取雅可比矩阵.由公式(1)~(3)可以得到不同装置类型电磁响应对模型变量的偏导数如下:
水平共面装置
(19) |
垂直同轴装置
(20) |
这里pl是第l层介质的模型参数,它可以是电导率σ也可以是层厚度t.事实上,频率域航空电磁法一维正演公式中地下介质模型参数(电导率、层厚度)仅在反射系数RTE中出现,我们可以把对各种装置类型电磁响应求偏导数简单归结为对反射系数RTE求偏导,这样雅可比矩阵的求取实际上转化为反射系数RTE对模型参数求偏导数,后者可以通过下面通用的解析公式简单求出:
(21) |
(22) |
(23) |
公式(22)中ul对模型参数pl的偏导数可以由公式(5)、(6)直接求得.
公式(19)~(23)就是解析法求频率域航空电磁法雅可比矩阵J的解析公式.相对而言,解析法求解雅可比矩阵具有更高的计算精度,本次软件编程就采用该算法.
2.3 初始模型赋值解的非唯一性是地球物理资料反演中最棘手的难题之一.由于地球物理多解性的存在,马奎特反演对初始模型要求很高.如果初始模型选择不合适,反演过程很可能出现“发散”.因此,只有给定了合适的初始模型,马奎特反演才有可能收敛于实际地电模型附近.可见,初始模型赋值是马奎特反演最重要的技术前提之一.
理论上,马奎特方法可反演水平多层地电模型参数.事实上,由于直升机航空电磁系统每种装置仅有三个工作频率(收录6个实、虚分量),对于三层以上的地电模型得到的是“欠定方程”,会造成反演过程的不稳定,无法形成实用化的算法.另一方面,考虑到该类装置的理论探测深度通常在120~180 m之间[7~9],三层模型已基本可以满足大多数地质问题的需要,本次反演研究仍以三层模型为准(二层模型可视为三层模型的特例).初始模型的赋值方法在传统航空电磁资料反演中缺少系统的研究,因而会造成反演结果的不确定性,本次研究中将采用质心深度近似反演结果对初始模型进行赋值.这样做的好处在于:首先,质心深度算法简单,耗费计算资源少;其次,这样的赋值结果在一定程度上反映了实际电阻率分布特征,模型更为合理,计算结果更为准确.经过反复的对比测试,本次研究中初始模型采用以下方法赋值:
(24) |
(25) |
(26) |
(27) |
(28) |
其中,ρ1、ρ2、ρ3分别为第一、二、三层的电阻率,d1、d2分别为第一、二层的厚度;ρah、ρam、ρal分别为实虚分量法计算得到的高、中、低频视电阻率,zah、zal分别为高、低频计算的质心深度[1, 10].
上面给出的赋值方法是在大量理论模型对比分析的基础上提出来的,对于常见的地电结构是适用的.地球物理反演中初始模型的选择常带有一定的偶然性和经验性,没有一种成熟方法.一般在实测数据的反演中,当前一测点数据能稳定收敛时,用其反演结果作为下一测点的初始模型.
2.4 反演模型选择和其他地球物理反演一样,航空电磁数据反演首先也需要确定选择什么样的反演解释模型.反演模型通常分为两类:形态模型和物性模型[11].形态模型是指模型的形状、物性参数均为变量,反演过程要通过改变模型形态、物性等来达到模型数据与观测数据的最大拟合;物性模型是将场源区域划分成小的单元组合(一维反演中为层状模型),在反演过程中单元的形态不变,通过物性变化勾画出场源范围.
对于一维层状反演而言,形态模型反演参数多(电阻率、层厚度),模型变化灵活,反演结果为物性差异明显的层状介质;缺点是对初始模型依赖程度高,反演过程不好控制.物性模型反演则相反,具有模型参数少(仅包含电阻率参数)、初始模型赋值简单(选用均匀半空间模型即可)、反演方法技术受限制条件少以及不涉及复杂的形态变化等优点;缺点是物性连续分布,反演结果只是实际电阻率分布特征的一种近似,电性层间界限模糊.
本次研究对以上两种算法均进行了编程,并在用户界面(窗体)中提供开关选项,供解释人员自由选择.
3 程序流程图频率域吊舱式直升机航空电磁系统马奎特反演计算流程如图 2.
为了检验马奎特反演应用效果和反演算法的正确性,以IMPULSE吊舱式航空电磁系统为例,运用新开发算法对理论模型和实测数据进行测试计算.
4.1 模型计算模型一:由一系列首层厚度逐渐增大的一维三层模型来模拟高阻介质中的倾斜低阻层(图 3a).其中第一、三层为高阻层,即ρ1=ρ3=500Ωm;第二层为低阻层,电阻率ρ2=30Ωm,层厚度t2=20m;模型第一层厚度逐渐增大,t2=10~100 m;飞行高度(吊舱离地高度)h=30m;点距s=2m.
模型二:由一系列首层厚度逐渐增大的一维三层模型来模拟一个良导介质中的倾斜高阻层(图 4a).其中第一、三层为低阻层,即ρ1=ρ3=30Ωm;第二层为高阻层,电阻率ρ2=500Ωm,层厚度t2=20m;第一层厚度逐渐增大,t2=10~100m;飞行高度(吊舱离地高度)h=30m;点距s=2m.
首先,对上述理论模型进行正演计算[12, 13],以其计算结果作为观测数据,运用质心深度近似反演结果对初始模型进行赋值;然后,运用新开发算法进行马奎特反演计算.反演过程中,当前一测点数据稳定收敛后,即将该测点计算结果作为下一测点初始模型进行新的迭代计算.同时,为了对比不同类型反演应用效果,分别对两个模型进行了形态模型与物性模型的马奎特反演计算.模型及反演结果分别如图 3、4所示.
从两个理论模型的反演结果看,在初始模型给定合适的情况下,形态模型反演能够很好地反映出电阻率的分布情况,并且能够很快收敛到“最优解”附近.反演结果电性层之间的界限清楚,与理论模型对应较好.物性模型反演也基本反映出了模型电性分布特征,并且反演结果受初始模型影响小,反演过程更为稳定;缺点是反演结果与实际电阻率值有一定差异,且电性层界限模糊,不易区分.
4.2 实例分析实测数据测试选用IMPULSE系统在我国西北某地区获得的部分实测数据.首先,对实测剖面数据进行质心深度近似反演计算,将计算结果作为输入参数,运用本次研究给定方法(公式(24)~(28))计算初始模型电阻率与层厚度,作为形态模型反演初始模型;以均匀半空间(ρ=100Ωm)作为物性模型反演初始模型.运用新开发算法分别进行两种反演模型(形态模型与物性模型)马奎特反演计算,并对计算结果进行对比分析.从实测数据反演结果(图 5)看:在3400~3700m处存在明显的高阻异常,视电阻率300~600Ωm,地面踏勘证实为地下煤火区(X)的反映.通过电磁定量解释地下煤火宽400m左右,从地表向下延深约60m,电阻率大于200Ωm,推断由地下煤火燃烧产生的裂隙、空洞(群)引起;在1600~1800m处存在相对高阻异常,为另一处煤火区(Ⅷ)在本剖面的反映.通过电磁定量解释地下煤火宽度约200m左右,从地表向下延深约30 m,电阻率约90Ωm以上;在6800~7600 m处的高电阻率异常,在地质剖面上有奥陶系(O1s)出露,认为是高阻的奥陶系地层的反映[14, 15].
从实测数据的反演过程及不同反演类型(形态模型、物性模型)计算结果对比,还可以看出以下几点:
(1)在通常情况下,采用本次研究中给定的初始模型赋值方法,选用合适的反演参数,马奎特反演能够稳定收敛,并且计算速度较快.
(2)两种不同类型反演结果大体一致,局部差别源于算法上的差异和显示方式上的不同.
(3)形态模型反演对初始模型依赖性强,反映出的异常信息更为丰富;物性模型反演初始模型赋值简单,但反演结果较为粗略,体现了真实模型的某种近似.
(4)通过对反演剖面叠加上地形信息,使图形显示更为直观,也更符合解释人员的使用习惯.
5 结论本文详细介绍了马奎特反演基本原理及实现方法,重点分析了雅可比矩阵的不同求解方法,并对初始模型赋值方法进行探讨,给出了行之有效的计算办法;随后,详细对比了形态模型反演和物性模型反演技术特点,通过软件编程形成实用化的算法;最后,通过理论模型和实测数据的对比测试,分析算法选择的合理性和计算结果的正确性.
事实上,地下电阻率的分布是非常复杂的,仅仅用简单的水平层状模型(通常为三层模型)来模拟实际三维分布的地质体是有局限性的.如果初始模型选择不当,还容易造成反演过程出现“发散”,也就是说,马奎特反演对测区地质条件的要求是比较高的.另外,马奎特反演效果还会受到地质情况复杂程度和收录数据的精度等因素影响[16~19].通常在某一测点两侧一定范围内的介质为近似水平层状时,即测点两侧的水平地层延伸足够远,非水平层的影响可以忽略不计时,马奎特反演总能得到比较满意的结果.当然,马奎特反演还有其最显著的优点:如算法简单,计算速度快、精度高,反演结果(形态反演)电性层界面清晰等,这是其他反演方法(如质心深度近似反演、Occam反演等)所不具备的,在探测冻土层厚度、寻找地下淡水、探测海水入侵范围等工程及地质勘查中拥有广阔的应用前景.
[1] | 周道卿. 频率域航空电磁数据的质心深度近似反演. 物探与化探 , 2007, 31(3): 242–244. Zhou D Q. Inversion of FAEM data with the centriod depth theory. Geophysical and Geochemical Exploration (in Chinese) , 2007, 31(3): 242-244. |
[2] | 米萨克N纳比吉安著.勘查地球物理--电磁法.赵经祥译.北京:地质出版社, 1992 Misac N.Nabighian.Electromagnetic Methods in Applied Geophysics (in Chinese).Zhao J X tran.Beijing:Geological Publishing House, 1992 |
[3] | Zhang Zhiyi, Partha S Routh, Douglas W Oldenburg. Reconstruction of 1-D conductivity from dual-loop EM data. Geophysics , 2000, 65(2): 492-501. DOI:10.1190/1.1444743 |
[4] | Huang Haoping, Fraser Douglas C. Inversion of helicopter electromagnetic data to a magnetic conductive layered earth. Geophysics , 2003, 68(4): 1211-1223. DOI:10.1190/1.1598113 |
[5] | 高惠璇. 统计计算. 北京: 北京大学出版社, 2003 . Gao H X. Statistical Computing (in Chinese). Beijing: Peking University Press, 2003 . |
[6] | 王家映. 地球物理反演理论. 北京: 高等教育出版社, 2002 . Wang J Y. Inverse Theory in Geophysics (in Chinese). Beijing: Higher Education Press, 2002 . |
[7] | 王卫平, 王守坦. 直升机频率域航空电磁系统在均匀半空间上方的电磁响应特征与探测深度. 地球学报 , 2003, 24(3): 285–288. Wang W P, Wang S T. Electromagnetic response character of helicopter frequency domain EM system above uniform half space and its prospecting depth. Acts Geoscientica Sinica (in Chinese) , 2003, 24(3): 285-288. |
[8] | Huang Haoping, Fraser Douglas C. The differential parameter method for multifrequency airborne resistivity mapping. Geophysics , 1996, 61(1): 100-109. DOI:10.1190/1.1574674 |
[9] | Mundry E. On the interpretation of airborne electromagnetic data for the two-layer case. Geophysical Prospecting , 1984, 32(2): 144-172. |
[10] | Sengpiel K P. Approximate inversion of airborne EM data from a multilayered ground. Geophysics Prospecting , 1988, 36(4): 446-459. DOI:10.1111/gpr.1988.36.issue-4 |
[11] | 姚长利, 郝天珧, 管志宁. 重磁反演约束条件及三维物性反演技术策略. 物探与化探 , 2002, 26(4): 253–257. Yao C L, Hao T Y, Guan Z N. Restrictions in gravity and magnetic inversion and technical strategy of 3D properties inversion. Geophysical and Geochemical Exploration (in Chinese) , 2002, 26(4): 253-257. |
[12] | 周道卿, 谭捍东, 王卫平. 频率域航空电磁资料Occam反演研究. 物探与化探 , 2006, 30(2): 162–165. Zhou D Q, Tan H D, Wang W. P. Occam's inversion method in FAEM datum processing.Geophysical and Geochemical Exploration (in Chinese) , 2006, 30(2): 162-165. |
[13] | Verma Saurabh Kumar, Sharma S P. Focused resolution of thin conducting layers by various dipole EM systems. Geophysics , 1995, 60(2): 381-389. DOI:10.1190/1.1443774 |
[14] | 中国国土资源航空物探遥感中心编.内蒙古乌达地区直升机航空电磁测量试生产成果报告.2007 China Aero Geophysical Survey and Remote Sensing Center for Land and Resources.Report on trial production with the helicopter airborne electromagnetic system in Wuda area of Inner Mongolia (in Chinese).2007 |
[15] | 熊盛青, 陈斌, 于长春, 等. 地下煤层自燃遥感与地球物理探测技术. 北京: 地质出版社, 2006 . Xiong S Q, Chen B, Yu C C, et al. Remote Sensing and Geophysical Exploration Technology on the Spontaneous Combustion of Underground Coal (in Chinese). Beijing: Geological Publishing House, 2006 . |
[16] | Yin Changchun. Fraser Douglas C, Attitude corrections of helicopter EM data using a superposed dipole model.. Geophysics , 2004, 69(2): 431-439. DOI:10.1190/1.1707063 |
[17] | Farquharson Colin G, Oldenburg Douglas W, et al. Simultaneous 1D inversion of loop-loop electromagnetic data for magnetic susceptibility and electrical conductivity. Geophysics , 2003, 68(6): 1857-1869. DOI:10.1190/1.1635038 |
[18] | Sattel Daniel. Inverting airborne electromagnetic (AEM) data with Zohdy's method. Geophysics , 2006, 70(4): G77-G85. |
[19] | Beard Les P, Nyquist Jonathan E. Simultaneous inversion of airborne electromagnetic data for resistivity and magnetic permeability. Geophysics , 1998, 63(5): 1556-1564. DOI:10.1190/1.1444452 |