2014, Vol. 57
航空电磁法(Airborne Electromagnetic，简称AEM )起始于20世纪40年代末.50年代，频率域航空电磁法在加拿大被首次应用于矿产勘查，并发挥重要作用(王卫平，2011).经过60多年的发展，频率域航空电磁已取得很大进步，被广泛应用于地质填图、金属矿普查、油气、地热、水文及环境工程勘查等领域(雷栋等，2006).
由于航空电磁法数据量巨大，通常采用成像技术快速获得地下电性分布信息.Sengpiel(1988)提出了质心深度成像法，这种方法利用传播函数与视电阻率之间的关系得到一个近似电阻率剖面.质心深度成像虽然在良导介质的航空电磁数据处理中取得了较好的效果，然而该方法不能准确地确定地质体的真实埋深.Huang and Fraser(1996)提出了微分电阻率成像方法，该方法在Sengpiel成像法的基础上，将视电阻率和视深度转化为微分电阻率和深度，从而得到电阻率随深度变化光滑的地电成像模型.与Sengpiel法相比，微分电阻率法能够更加准确的确定地下异常体的埋深.
为了得到更加准确的地电信息，需要对航空电磁数据进行反演.Marquardt方法是一种应用十分广泛的非线性最小二乘反演方法.通过加入阻尼因子提高反演的稳定性，加快了收敛速度.这种反演方 法反演速度快，电性界面清晰，但是对初始模型的要 求较高，容易陷入局部极小(周道卿等，2010). Constable等(1987)首次将Occam技术应用于电磁数据反演当中.Occam反演在满足数据拟合的同时，对反演模型的光滑度设限.该方法反演稳定，对初始模型的要求较低，收敛速度快，但地层分界面不够清晰.Yin and Hodges(2007)提出航空电磁模拟退火反演，这种方法模拟流体冷却结晶的热力学过程，从Boltzmann分布沿上行和下行双向搜索，从而可以跳出局部极小值搜索到全局最小值，得到真实的反演解.这种方法对反演初始模型要求不高，且无需计算雅克比矩阵，大幅减少反演的计算量.然而，该方法在模拟退火过程中初始温度及下降步长的选择取决于模型参数，增加了反演难度.
以上的成像技术和反演方法都是基于横向电阻率均匀条件下的单点数据处理.实际地质情况并不是电阻率横向均匀的，利用简单的水平层状模型模拟实际二、三维地质体分布是不合理的，应采用二维或三维电磁反演方法.然而，二维和三维反演十分复杂，计算量大，特别是对于移动场源和接收的航空电磁法，实现二、三维数据反演目前仍非常困难.为了能较准确、快速灵活地反演地下电性分布信息，提出了横向约束反演方法(LCI).LCI方法是一种拟二维反演方法.它是通过将剖面数据集成反演，并施加横向约束，同时反演出同一测线上多测点的地电参数.
横向约束反演方法首先由Auken等(2000) 引入到直流电法数据反演中； Auken and Christiansen(2004) 利用LCI方法对2D电磁数据进行了反演处理；Santos (2004)进一步将LCI方法应用于EM34 剖面数据的反演之中； Christiansen and Auken(2004)利用Broyden 更新近似计算雅克比矩阵，对二维LCI反演进行优化，并利用LCI方法对电阻率数据进行 分段反演处理(Auken et al., 2005)； Siemon等(2008) 利用LCI反演方法对直升机吊舱数据进行处理；2008年，Auken再次利用LCI方法处理瞬变 电磁数据进行古河道勘查，而Triantafilis and Santos (2010)将LCI方法应用于石油勘探.
虽然国外对横向约束反演方法已经进行了卓有成效的研究，国内该方法还没有引起足够重视.本文在国内首次提出改进的加权横向约束反演方法(WLCI)，将其应用于频率域航空电磁数据反演中.考虑到航空电磁法数据采样密集，具有较高的横向分辨率，这为最大限度的发挥WLCI方法的优点提供了前提.本文中首先阐述频率域航空电磁法正演理论，详细介绍WLCI方法的原理及其在航空电磁数据反演中的应用技术，重点介绍如何通过层参数加权约束改善反演的地下电性界面的连续性.最后，通过对理论模型和实测数据的反演处理，同时和传统的反演方法进行对比，以检验加权横向约束反演方法的有效性.
目前，用于环境和工程勘查的航空电磁系统主要采用水平共面(HCP)装置.如图1所示，假设坐标原点位于发射线圈正下方的地面上，Z轴向下为正.磁偶极子高度为h0，地下水平地层的电阻率和 厚度分别为ρi(i=1,2,…，N)和hi(i=1,2,…，N-1), N为模型层数.
 | 图1 航空电磁模型示意图 Fig.1 Airborne EM system over a layered earth |
根据Yin和Hodges(2007)，对于水平共面(HCP)装置， 在测点处沿x方向的水平磁场分量为
而磁场垂直分量为
其中，z±=h0±z，R+= r2+z2+ ，r为收发距，m是 发射线圈的偶极矩，T(z－)是汉克尔积分，由下式给出：
其中,η是核函数，与空间波数和地电参数有关(Yin et al.,2007).
航空电磁通常采用HCP装置，观测若干个频 率的磁场数据.观测结果经一次场归一化后，以PPM (Parts Per Million)为单位表示.本文以5个频率航空电磁观测系统为例.将磁场观测结果作为反演数据，以矢量形式表示为
其中, d obsi表示第i个测点的实测数据，Hzij表示第i个测点第j个频率的垂直磁场响应值.由此，整个剖面的反演数据可以表示为 d obs=( d obs1, d obs2,…, d obsM)T， 其中M表示测点个数.
反演的模型参数由地下各层的电阻率和层厚度构成.为了确保反演参数为正， 同时各反演变量可在(－∞,+∞)之间变化，我们对各模型参数取对数 (Yin，2000)，将第i个测点的模型参数写成矢量形式：
其中, p i表示第i个测点的模型参数向量，ρij表示第i个测点第j层的电阻率，hij表示第i个测点第j层的层厚度.由此，全部反演模型参数可表示为 m =(m1,m2,…,mM*(2N－1))T=( p 1, p 2,…, p M)T.
雅克比矩阵或灵敏度矩阵，表示模型参数对正演结果的影响程度.它由正演响应函数对模型参数的偏导数组成.对航空电磁系统，沿飞行测量剖面第k个测点的雅克比矩阵可以写成:
式中，对于HCP装置
其中,Hzi为第i个频率的垂直磁场，pj是第j个模型参数.
反演过程是通过调整模型参数实现理论模型正演结果与实测数据拟合的过程.如上所述，电磁响应与地电参数之间存在复杂的非线性关系，为简化计算，我们对其做一阶泰勒展开：
其中,F是电磁响应函数， d obs是观测数据， m 是反演迭代模型， m 0是初始或参考模型， e obs是截断误差， G 是雅克比矩阵.根据Auken(2004)，上式可以简写为
其中,
Gi表示第i个测点的雅克比矩阵，参见(7)式.
LCI反演的中心思想是使相邻测点地电参数的 差别尽可能小，以达到光滑目的.根据Auken(2004)， 我们有
两边减去 R p m 0得到
其中
e rp 表示相邻测点模型参数间的差异，横向约束矩阵 R p 是由1和-1构成的稀疏矩阵，可表示为
其中 S =(M-1)×(2N-1)， T =M×(2N-1)，N是模型层数，M 是测点个数.
在实际反演中，为突出层参数横向连续的重要性，可通过引入加权因子改变对各层光滑度的横向约束强度，实现加权横向约束反演WLCI.为此，将(16)式两边乘以权系数矩阵，得到
其中 W p是横向约束加权矩阵，(21)式中的加权等价于(19)式中的 R p矩阵相应的行乘以加权系数.参数加权的大小根据实际需要设定.
上面讨论的约束条件只是对单个层参数进行.为了保证层界面的光滑和连续性，有时需要对层界 面深度进行约束.根据Auken(2004)，我们可以得到
两边减去 R t m 0得
其中 e rt表示相邻测点层界面深度差异， R t是一个与深度和层厚有关的稀疏矩阵，可表示为
其中hi,j表示第i个测点第j层的厚度，ti,j表示第i 个测点，第j层下界面的深度， S =(M-1)×(N-1), T =M×(2N-1).
对(24)式中各层深度横向约束的加权可采取前面参数约束相同的模式，通过选取适当的加权系数来调整各层层界面连续的重要性.为此，我们引入加权矩阵，重组(24)式得到
其中
将数据拟合方程(11)、电阻率和厚度约束方程(20)以及深度约束方程(27)集成, 得到在加权横向约束条件下的总体反演方程:
进一步简化为
实际反演中可利用最小二乘法求解(31).通过对误差函数的平方求极小，我们得到如下法方程:
其解为
(32)和(33)式中的系数矩阵通常存在奇异性，导致反演解不稳定.为此，在系数矩阵中加入阻尼因子λ2，将法方程变为
其中， I 是单位矩阵.为求解(34)式，将系数矩阵 A 进行奇异值(SVD)分解：
其中 U和V 分别为数据和参数矩阵， Λ 为奇异值矩阵.将(35)式带入(34)式并化简得
为了选取最优的阻尼因子，将λ的初始值设为奇异值矩阵 Λ 对角线元素的最大值，然后按照每次减小 2 的规律进行搜索，当拟合差不再减小时，此时的λ即为此次迭代所选取的阻尼因子.考虑到实际反演问题的非线性，从一个选定的初始模型开始，利用方程(36)进行迭代，直到误差函数小于一个事先设定的阀值.实际WLCI反演流程如图2.
 | 图2 WLCI反演流程图 Fig.2 Flowchart of WLCI inversion |
为了验证WLCI方法的有效性，本文首先对理论模型数据进行反演，并与SVD单点反演结果进行 对比.理论模型数据是针对HCP装置计算的垂直 磁场 响应.频率设为 380、1600、6300、25 kHz和120 kHz， 收发距为8 m，飞行高度为30 m.如图3a所示，设计的理论模型是四层水平层状大地，第一层是电阻率为200 Ωm，层厚为20 m，在横向30—50号点之间存在着50 Ωm的低阻异常区域；第二层是电阻率为100 Ωm的高阻层，层厚为30 m；第三层是电阻率为 5 Ωm低阻层，层厚为10 m;第四层是电阻率1000 Ωm 的高阻基底.初始模型选择电阻率为200 Ωm的均匀半空间，单点SVD反演和LCI反演结果如图3(b—d).实际反演解释中针对不同电性参数(ρ1,h1,t1,ρ2,h2,t2,ρ3,h3,t3,ρ4)采用加权以调整不同地电参数横向连续 性：图3c中的加权因子 取为(0.9, 2, 2, 0.9, 0.9, 0.9, 0.9, 10, 0.9, 0.9)， 而图3d中的加权因子取为(0.9, 2, 2, 0.9, 0.9, 0.9, 0.9, 2, 0.9, 0.9).
从图3可以看出: (1) 利用SVD单点反演方法和WLCI方法均能反映出地下电性分布情况，显示为四层地电结构，第一层中间部分和第三层的低阻异常均得到不同程度的反映.值得指出的是，为了改善SVD单点反演效果，我们将前一测点的反演结果 作为后一测点的初始模型，这在一定程度上改善了SVD反演结果的横向连续性.然而，从图可以看出，SVD单点反演结果中各层电阻率的连续性仍然较差，各层层界面比较粗糙，电阻率条带明显；相比之下， WLCI反演结果中各层的连续性得到很大改善，层界面光滑，反演的各层电阻率横向分布比较均匀; (2) 虽然图3c和3d的WLCI反演均得到光滑的反演界面，然而，由于采用的加权因子不同，取得的效果不同.图3c中，由于加权因子过大，各层的电性沿水平方向被过度平均，导致第一、二层的电性分布模糊，第三层厚度明显大于理论值；图3d由于采用的加权因子适中，反演的电性分界面连续，电性分布和层厚度非常接近真实模型.这说明选择合适加权因子的重要性.实际反演中，加权因子的选取基于地质条件，对于连续性较好的地层，可适当地加大加权因子；而对于成层性不明显或连续性差的地层，应减小加权因子.根据理论模型和实际资料反演得到的经验，加权因子的合适选取范围为0.1~10.
 | 图3 层状模型及反演结果 (a)理论模型；(b)SVD单点反演；(c) WLCI反演(大加权因子)；(d) WLCI反演(加权因子适中). Fig.3 Layered earth model and inversion results (a) Synthetic model; (b) SVD inversion; (c) WLCI inversion with strong weighting；(d) WLCI inversion with approriate weighting. |
通过对航空电磁用于路基稳定性调查的观测数据反演, 验证WLCI反演方法的有效性. 实际观测 采用HCP装置测量垂直磁场分量，频率为 380 Hz～ 100 kHz.图4给出了利用SVD、WLCI反演方法和CDT成像技术(Macnae et al,1991)对实测数据的反演结果.各测点的反演初始模型均假设相同的三层模型：第一层电阻率为8 Ωm，层厚20 m，第二层 电阻率为4 Ωm，层厚10 m，第三层电阻率为15 Ωm. 针对不同电性参数(ρ1,h1,t1,ρ2,h2,t2,ρ3)采用加权因子(1, 1，10, 1, 1, 1, 1).由图可以看出，利用各种反演方法得到的结果均不同程度显示地下电性分 布：第一层是高阻层，电阻率为5~9 Ωm，厚度为 20 m左右，通过对路基的实际勘查，确认为路基顶 部的泥沙层；第二层是低阻薄层，电阻率为2~5 Ωm， 厚度为10 m左右，确认为路基中部的粘土层.该层是本次勘探的目标，其存在直接影响路基的稳定性.第三层是高阻基底，电阻率为11~17 Ωm，确认为路基底部的含沙地层.对比各种反演方法的结果不难发现，虽然SVD单点反演的结果基本反映出了地下电性分布，但各层的连续性较差，地层界面不光滑，下半空间的电阻率反演结果中不均匀条带尤其明显；WLCI反演结果中各层的连续性得到很大改善，地层界面比较光滑；下半空间电阻率不均匀条带明显减少，很好地反映了路基的稳定地质条件.另外，从施加的约束条件和地电断面反演结果可以进一步看出，对第一层底界面深度设定强约束条件(加权因子较大)可有效地改善整个反演结果，不仅各层界面的连续性得到改善，基岩中横向电性不均匀性也得到较大改善，反演结果与CDT成像结果吻合较好.
 | 图4 航空电磁实测数据反演结果 (a) SVD反演；(b) WLCI反演; (c) CDT成像. Fig.4 The inversion results of survey data (a) SVD inversion; (b) WLCI inversion; (c) CDT imaging. |
通过系统阐述频率域航空电磁的正演和加权横向约束反演理论，本文建立了一套有效的航空电磁数据反演的方法技术.与传统的单点非线性最小二乘反演相比，加权横向约束反演WLCI可有效改善反演稳定性，相邻测点的电性反演结果比较连续，界面清晰.特别是通过引入层参数加权，实际反演地电断面的层界面更加光滑，横向电性分布均匀，电阻率假异常得到有效压制.横向约束加权因子的选择影响反演结果，加权因子太小对反演结果的约束作用很小；反之加权因子选的太大，电性沿横向方向被过度平均.实际反演中，加权因子的选择应遵循从已知到未知的原则：通过对钻孔附近剖面性观测资料进行反演，确定合适的加权因子，推广应用到整个测区.对于没有相关地质资料的盲区，可根据航空电磁数据的分辨能力和一般地质规律，地层埋藏越深，地电参数横向相对越光滑，加权因子越大.和其他电磁成像(比如CDT成像)结果进行对比反演，也有益于选择最优的横向约束加权因子.
对实测航空电磁数据的反演不仅证实了WLCI反演技术的有效性，同时也证实了航空电磁法在工程勘查方面的有效性.然而，WLCI反演的前提是相邻测点地下电性结构变化不大，因而可以通过施加横向约束来保证反演参数的连续性和反演界面的光滑性.由于航空电磁系统采用密集型采样，具有很高的横向分辨率，为WLCI方法的成功应用提供保证.WLCI方法对于其他地球物理数据反演的适应性，取决于其数据采样的密集程度.
WLCI反演虽然不属于严格的二维反演，然而，由其反演所获得的地电模型较之于传统方法更接近真实模型，可作为复杂的二维电磁反演的初始模型，改善反演解的稳定性.未来的研究方向将侧重于二维面积性电性参数的横向约束，以期为未来的三维电磁反演提供好的初始模型.
致 谢 向刘云鹤博士对于本文准备和撰写过程中所提供的帮助表示感谢.同时对各位审稿人对本文提出的修改意见表示衷心感谢.| [1] | Auken E, Christiansen A V. 2004. Layered and laterally constrained 2D inversion of resistivity data. Geophysics, 69(3): 752-761. |
| [2] | Auken E, Christiansen A V, Jacobsen Bo H, et al. 2005. Piecewise 1D laterally constrained inversion of resistivity data. Geophysics, 53(4): 497-506. |
| [3] | Auken E, Christiansen A V, Jacobsen L H, et al. 2008. A resolution study of buried valleys using Laterally constrained inversion of TEM data. Journal of Applied Geophysics, 65(1): 10-20, doi: 10.1016/j.jappgeo.2008.03.003. |
| [4] | Auken E, Srensen K I, Thomsen P. 2000. Lateral constrained inversion (LCI) of profile oriented data—the resistivity case. Proceedings of the EEGS-ES, Bochum, Germany, EL06. |
| [5] | Christiansen A V, Auken E. 2004. Optimizing a layered and laterally constrained 2D inversion of resistivity data using Broyden's update and 1D derivatives. Journal of Applied Geophysics, 56(4): 247-261, doi: 10.1016/j.jappgeo.2004.07.004. |
| [6] | Constable S C, Parker R L, Constable C G. 1987. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3): 289-300. |
| [7] | Huang H P, Fraser D C. 1996. The differential parameter method for multifrequency airborne resistivity mapping. Geophysics, 61(1): 100-109. |
| [8] | Lei D, Hu X Y, Zhang S F. 2006. Development status of Airborne Electromagnetic. Contributions to Geology and Mineral Resources Research, 21(1): 40-44, 53. |
| [9] | Macnae, J, Smith R, Polzer B, et al. 1991. Conductivity-depth imaging of airborne electromagnetic step-response data. Geophysics, 56(1): 102-114. |
| [10] | Santos F A M. 2004. 1-D laterally constrained inversion of EM34 profiling data. Journal of Applied Geophysics, 56(2): 123-134, doi:10.1016/j.jappgeo.2004.04005. |
| [11] | Sengpiel K P. 1988. Approximate inversion of airborne EM data from a multilayered ground. Geophysical Prospecting, 36(4): 446-459. |
| [12] | Siemon B, Auken E, Christiansen A V. 2009. Laterally constrained inversion of helicopter-borne frequency-domain electromagnetic data. Journal of Applied Geophysics, 67(3): 259-268, doi: 10.1016/j.jappgeo.2007.11.003. |
| [13] | Triantafilis J, Santos F A M. 2010. Resolving the spatial distribution of the true electrical conductivity with depth using EM38 and EM31 signal data and a laterally constrained inversion model. Australian Journal of Soil Research, 48(5): 434-446, doi:10.1071/SR09149. |
| [14] | Wang W P. 2011. Application and Prospect of the Research about Frequency Domain Airborne Electromagnetic Method. Twenty-seventh Annual Meeting Proceedings of the Chinese Geophysical Society, 675-676. |
| [15] | Yin C C. 2000. Geoelectrical inversion for a one-dimensional anisotropic model and inherent non-uniqueness. Geophys. J. Int., 140(1): 11-23. |
| [16] | Yin C C, Hodges G. 2007. Simulated annealing for airborne EMinversion. Geophysics, 72(4): F189-F196. |
| [17] | Zhou D Q, Tan L, Tan H D, et al. 2010. Inversion of frequency domain helicopter-borne electromagnetic data with Marquardt's method. Chinese J. Geophys. (in Chinese), 53(2): 421-427, doi: 10.3969/j.issn.0001-5733.2010.02.020. |
| [18] | 雷栋,胡祥云, 张素芳. 2006. 航空电磁法的发展现状. 地质找矿丛论, 21(1): 40-44, 53. |
| [19] | 王卫平. 2011. 频率域航空电磁法发展研究现状与应用研究. 中国地球物理学会第二十七届年会论文集, 675-676. |
| [20] | 周道卿, 谭林, 谭捍东, 等. 2010. 频率域吊舱式直升机航空电磁资料的马奎特反演. 地球物理学报, 53(2): 421-427, doi: 10.3969/j.issn.0001-5733.2010.02.020. |