地球物理学进展  2014, Vol. 29 Issue (6): 2723-2729   PDF    
引用本文
罗勇, 陆从德, 王宇航. 2014. 时间域航空电磁一维阻尼特征参数反演方法[J]. 地球物理学进展, 29(6): 2723-2729, doi: 10.6038/pg20140638  
LUO Yong, LU Cong-de, WANG Yu-hang. 2014. 1D inverse method based on damped eigenparameter for airborne time domain electromagnetic data. Progress in Geophysics, 29(6): 2723-2729, doi: 10.6038/pg20140638   
时间域航空电磁一维阻尼特征参数反演方法
罗勇1, 陆从德1 , 王宇航3    
1. 成都理工大学地球探测与信息技术教育部重点实验室, 成都 610059;
2. 中国地质大学地球物理与空间信息学院, 武汉 430074;
3. 四川省地震局, 成都 610000
收稿日期: 2014-03-18 修回日期: 2014-05-28 .
基金项目：国家863计划主题项目课题(2013AA063904)和四川省教育厅重点项目(11ZA045)联合资助.
作者简介：罗勇, 男, 1987年生, 硕士生, 研究方向为航空电磁正反演方法及应用.(E-mail:316242438@qq.com)
通讯作者：陆从德, 男, 1972年生, 博士, 副教授, 研究方向为地球物理信号处理.(E-mail:cdlu@163.com)
摘要：时间域航空电磁数据反演是当前航空电磁领域的研究热点和难点,甚至一维反演问题仍然没有得到很好的解决.本文应用阻尼特征参数算法对时间域航空电磁数据一维反演进行研究,即将阻尼因子与相对奇异值结合,通过改变相对奇异值修正阻尼因子,达到约束反演过程的目的.由于阻尼特征参数算法综合了奇异值截断和马奎特方法的优点,所以具有较高的计算效率和精度.为了验证该反演方法的有效性,使用MAXWELL软件中的2.5D正演方法获得反演所需的模拟数据.实验主要对三层和四层地电模型的模拟数据进行反演,结果表明该方法收敛稳定,对低阻体的反演结果较好.
关键词： 阻尼特征参数算法     奇异值分解     阻尼因子     余弦变换    
1D inverse method based on damped eigenparameter for airborne time domain electromagnetic data
LUO Yong1, LU Cong-de1 , WANG Yu-hang3    
1. Key Lab of Earth Exploration &Information Techniques of Ministry of Education, Chengdu University of Technology, Chengdu 610059, China;
2. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
3. Sichuan Earthquake Administration, Chengdu 610000, China
Abstract: The inverse problem of time-domain airborne electromagnetic data is still one of the hot and difficult spots in the field of airborne electromagnetic prospecting,even one-dimensional inverse problems cannot be solved very well.In this paper,the inverse problems of one-dimensional electromagnetic data may be solved by the damped eigenparameter method, which constrains the inverse process by combining the damping factors with the relative singular values and correcting iteratively damping factor. The method proposed is with low computational cost and high accuracy because it integrates the advantages of the Marquardt and Truncation estimates. Synthetic data generated by the 2.5D modeling program of the MAXWELL software are applied to demonstrate the validity of the inverse method shown. The experiments upon the typical geoelectric models with the three-layer and the four-layer media respectively show that the proposed method is with high accuracy and stable convergence speed, especially for low resistivity object.
Key words: damped eigenparameter inversion     SVD     damping factor     cosine transform    
0 引 言

由于时间域航空电磁法具有可大面积勘探、勘探深度较大以及相对成本较低等优点，因此该方法已广泛应用于金属矿产勘探、水资源勘查和地质填图等应用领域(Anderson et al., 1993；Worrall et al., 1998；Sattel and Kgotlhang， 2004).但是时间域航空电磁法的解释水平较低，其大规模推广应用受到了一定的限制，为此，时间域航空电磁法的反演解释研究成为了当前的研究热点和难点.由于时间域航空电磁法的数据量很大，在反演解释研究中，二、三维反演方法计算量又较大，因此目前主要集中在一维反演方法.

一维时间域航空电磁的解释方法主要包括如下几种方法：电导率深度成像(CDI)(Macnae et al., 1991；Liu and Asten， 1993；Wolfgram and Karlik， 1995)，薄板反演(Liu and Asten， 1993; Christensen，2002)，τ域成像(Stolz and Macnae， 1997)，迭代反演(Farquharson and Oldenburg， 1993；Christiansen and Christensen， 2003)等.前三种方法属于快速成像方法，不需要设置初始模型也不用进行迭代计算，计算速度快但结果精度较低，特别是在地形复杂地区，往往不能得到理想的效果；最后一种方法通常是对超定问题的求解，通过不断迭代拟合得到最符合的模型，一般都能得到较好的效果，但是计算结果受初始模型的选择影响较大，如阻尼最小二乘方法(Raiche et al., 1985；Huang and Palacky， 1991)，模型交替调整反演(李永兴等，2010)，神经网络方法(朱凯光等，2010)，自适应正则化反演(毛立峰等，2011)等.鉴于地下介质是三维分布，为了获得更加符合实际的反演结果，Wilson等(2006)和Ellis(1995)尝试用高斯-牛顿法分别对2.5维和三维时间域航空电磁数据反演做了研究，求解小模型的反演问题.尽管2.5维和三维时间域航空电磁反演具有较高的计算精度，但是其计算效率太低的问题一直没有得到较好的解决，还处于理论研究阶段，因此一维反演仍然是时间域航空电磁法解释的主要手段.阻尼特征参数算法结合了马奎特方法和奇异值截断的优点(Jupp and Vozoff， 1975)，将经奇异值分解后得到的特征值分为重要参数和不重要参数，截断小的奇异值，通过改变奇异值来修正阻尼因子.相对于阻尼最小二乘方法从两次计算结果来改变阻尼因子，该方法综合所有的数据特性，从整体上指导反演过程的进行，具有快速稳定的特点，在直流电阻率法和大地电磁的联合反演中取得了较好的效果(Vozoff and Jupp， 1975).本文首先计算了水平层状大地上方垂直磁偶极源的正演响应，通过余弦变换和Guptasarma线性数值滤波转换得到时间域的结果，然后用阻尼特征参数算法实现了时间域航空电磁一维反演. 1 正演理论

时间域航空电磁正演模型如图 1所示，在水平层状大地上水平放置发射线圈Tx，接收线圈Rx，水平收发距为r，发射线圈和接收线圈距地面高度分别为h和z.第n层的电阻率为ρn，厚度为dn.当向发射线圈供电流I时，可以得到磁场的垂直分量Hz的表达式(罗延钟等，2003)

其中，m表示磁矩，J0是零阶贝塞尔函数，λ是积分变量，RTE是反射系数.
图 1 时间域航空电磁正演模型 Fig. 1 Sketch of ATEM forward model

在正演计算中有两个关键的问题：频时转换和汉克尔积分变换.这两个问题解决的好坏直接影响着正演计算结果，下面分别对这两方面进行简要分析. 1.1 频时转换

为了得到时间域响应必须将频率域正演计算结果进行转换.频时转换的方法很多，早期的计算采用快速傅里叶变换，但是该方法计算速度慢.Knight和Raiche(1982)用G-S变换计算了层状大地频率域到时间域的转换，后来又发展了余弦变换和数值线性滤波等方法(Guptasarma，1982).这些方法针对不同的地电条件和适用范围都有各自的优缺点，本文选取余弦变换计算频率域到时间域的转换，余弦变换公式为

其中ω表示频率，lm[E(ω)]表示电场的虚部.令x=ωt，根据Abramowitz和Stegun(1964)，，并带入(2)式，化简后得

上式中Wi表示滤波系数，合理的滤波系数选择不仅跟待变换函数有关，同时也要考虑计算精度和效率.针对时间域航空电磁数据特点，本文中n的取值范围为-200~99之间的300个系数，保证在满足精度要求前提下提高计算速度. 1.2 汉克尔积分变换

对于汉克尔积分变换

其中，K(λ)是待变换的函数，Ji是i阶贝塞尔函数，r是给定值，λ表示坐标轴中的横坐标，即λi=(1/r)×10[a+(i－1)s]，i=1，2，…，n，且a表示初始位移，s表示采样间隔.由此可以得到汉克尔变换的计算公式为

上式中Wi是滤波系数.Guptasarma和Singh(1997)给出了零阶和一阶的滤波系数，并与Anderson给出的系数进行了对比分析，结果表明该系数稳定，计算效率高. 2 反演方法

设时间域航空电磁观测数据向量和模型数据向量分别为 y i和xi，xi表示电阻率或厚度，且 y i和x i可分别表示为 y i=[y1，y 2，… y M]T，x i=[ x 1，x 2，… x N]T.

地球物理数据反演的目的是寻求观测数据和理论模型响应在某种度量准则下的最小值，使理论模型响应最大程度地吻合实测数据.从公式(1)可以看出，时间域航空电磁正演是一个非线性方程.为了简化计算，将Hz在x0处按泰勒级数展开并线性化，可得

y0i是初始模型x0j的正演响应，为书写方便，将上式写成矩阵形式

其中ε是残差向量，且ε=yi－y0i，J是雅克比矩阵，J=时间域航空电磁中的雅克比矩阵通常是一个奇异矩阵，对于任意的矩阵用奇异值分解可以写成 J=USVT，这里 U是数据 M×N的特征向量矩阵，V是模型N×N的参数特征向量矩阵，P是分解后的奇异值，且奇异值的大小为Sj(S1>S2>S3…).UUT=IM，VVT= IN，IM和IN是秩分别为M和N的酉阵.当采用(7)式计算时，如果奇异值Sj较小，计算会出现不稳定甚至不收敛的情况，为解决该问题，可以通过阻尼因子T改善反演过程，于是获得模型参数改正量的表达式为

其中P是雅克比矩阵J的秩，J+是J的广义逆，

在求取δP 的过程中，如何选择阻尼因子以及给出合适的迭代准则和反演结束标志，是稳定收敛、快速反演需要解决的关键问题.2.1 阻尼因子的选择

阻尼因子是用来控制修正方向和步长，选择合适的阻尼因子能减小条件数，起到阻止发散的作用，因此阻尼因子的选择非常重要.Jupp和Vozoff(1975)给出阻尼因子表达式为

其中，ki=Si/S1，μ是相对奇异值的阈值.通常先设定μ为较大阈值，如μ=0.1，避免小奇异值在反演过程中带来的振荡，与较大奇异值对应的模型参数称为重要参数.随着反演过程的进行，均方相对误差逐渐减小，相对奇异值阈值也减小，更多数据进入反演，此时的模型参数称为不重要参数.当相对奇异值阈值减小到事先给定的信噪比NSR=0.01时就不再改变.地球物理反演问题中，重要参数在反演过程中起着决定性的作用，但是不重要参数会导致迭代过程中的不稳定，引起解的振荡使其偏离真值.为了说明阻尼因子的变化情况，图 2给出一个三层模型的阻尼因子的变化曲线示意图，其中，模型参数设置为：ρ1=50 Ωm，ρ2=200 Ωm，ρ3=50 Ωm，d1=100 m，d2=50 m.从图中可以看出，在前14次迭代过程中，阻尼因子基本保持不变；但是，从第15次迭代开始，由于相对奇异值阈值的减小，因此阻尼因子也迅速减小，使更多的数据进入反演计算，直到阈值等于信噪比就不再变化.
图 2 K型地电模型的阻尼因子变化曲线 Fig. 2 Damping factor of the K-type
2.2 反演迭代准则和结束标志

对于反演计算中的迭代过程，总是希望在保证收敛的前提下，取较小的阻尼因子，以获得较大的步长，加快收敛速度.迭代成功与否的判断通常是采用某种度量准则比较两次迭代，若值减小则说明该次迭代成功，减小阻尼因子继续迭代，否则增大阻尼因子.本文采用均方根相对误差作为度量准则，若两次的均方相对误差之差小于两次残差的平方，则该次迭代成功.对单点反演，均方根相对误差综合了所有时 间道的数据，相对于阻尼最小二乘能更好反映数据的总体特性，从整体上指导反演过程.所采用的均方根相对误差公式为

VD和VM分别表示观测数据和理论数据，Wi是权重.对于反演结束通常有很多方法来进行判断，本文采用如下三种反演结束标准，即反演达到最大迭代次数，均方根相对误差小于事先给定的值和(RMSi+1－RMSi)<0.001RMSi+1.为了避免反演过程陷入死循环，只要满足其中条件之一就停止迭代，反演结束. 3 反演算例

用阻尼特征参数算法对三层和四层地电模型的2.5维的正演计算结果进行反演，其中地电模型可参见图 1.在所有的反演模型中，飞行高度为35 m，发射线圈面积531 m2，线圈匝数4匝，发射和接收线圈的水平收发距120 m.发射波形采用峰值电流为273 A的双极性梯形波，采样时间道为26道，OFFTIME时间是9.167 ms.为了说明反演方法的有效性，用Maxwell软件中2.5维时间域航空电磁正演计算结果作为一维反演的模拟数据. 3.1 三层模型

三层模型的反演计算包括H型和K型地电模型，其中H型地电模型中第一层至第三层的电阻率分别是100 Ωm，20 Ωm和200 Ωm，前两层的厚度分别为100 m和50 m；K型 地电模型中第一层至第三层的电阻率分别是50 Ωm，200 Ωm和50 Ωm，前两层的厚度分别为100 m和50 m.使用阻尼特征参数反演方法对上述的H型和K型地电模型进行反演，所得的反演结果如图 3所示，左右两列分别表示H型和K型地电模型反演结果，从上到下依次是：模型反演结果，反演模型响应，均方根相对误差，反演结果的残差变化.其中反演模型响应是对反演得到的模型参数进行正演计算得到的响应.由于基底的厚度是无限大，为了和真实模型对比，a，b中给出了50 m厚度的基底.在图 1a、1b中，蓝色表示真实模型，红色表示反演结果.从这两个图中可以发现，第一层中的反演曲线尽管还有些跳动但基本都能反映真实模型；H型曲线在低阻层出现过拟合的情况，而K型曲线则在高阻层表现出欠拟合.此外，由于时间域航空法对低阻层比较敏感，响应较大，而对高阻层的响应较小，并且受上下两层的影响，这可能会使得在反演过程中低阻层电阻率得到更小的值，高阻层的电阻率低于真值.为了说明反演结果的有效性，图 1c和1d分别给出了H型和K型的反演模拟响应和模拟数据，其中虚线表示模拟数据，实线为反演模拟响应.这两个图所给出的结果表明：反演结果逼近了真实模型，且相对误差小于1%.图 1e、1f和1g、1h分别给出了H型曲线和K型曲线的均方根相对误差和反演结果的残差变化.这些图所给出的结果表明：反演结果较好，且H型和K型的均方根相对误差分别为0.76%和1.51%，残差分别为3.78×10-3和1.0×10－3.此外，均方根相对误差和残差向量基本是一条衰减曲线，均方根相对误差的初始误差较大，这和初始模型的选择有关，随着反演过程的进行，均方根相对误差和残差都迅速减小直到反演停止.

图 3 H型和K型地电模型的反演结果
(a)和(b)表示模型反演结果；(c)和(d)表示反演模型响应；(e)和(f)表示均方根相对误差；(g)和(h)表示残差变化. Fig. 3 Inverse results of H-type and K-type The left and right column show the inverse results of H-type and K-type models respectively
(a) and (b)inverse results;(c) and (d)response of inverse model;(e) and (f)relative RMS error;(g) and (h)residual.
3.2 四层模型

四层模型的反演计算包括HK型和KH型地电模型，其中HK型中的第一层至第四层的电阻率分别是200 Ωm，10 Ωm，100 Ωm和10 Ωm，前三层的厚度均为50 m；KH型第一层至第四层的电阻率分别是10 Ωm，200 Ωm，5 Ωm和200 Ωm，前三层的厚度均为50 m.使用阻尼特征参数反演方法对上述的HK型和KH型地电模型进行反演，所得的反演结果如图 4所示，其中左右两列分别表示HK型和KH型的反演结果，从上到下依次是：模型反演结果(参见a和b)，反演模型响应(参见c和d)，均方根相对误差(参见e和f)，反演结果的残差变化(参见g和h).其中反演模型响应是对反演得到的模型参数进行正演计算得到的响应.对于HK型和KH型的反演结果，曲线形态大致符合真实模型，但是在两种模型的高阻层的反演结果要低于真值.两种模型基本在反演进行到11次的时候就停止迭代，反演模型响应，均方根相对误差和残差都呈快速衰减形态.反演结束时，HK和KH型曲线反演模型响应与正演模型响应的相对误差都小于1%，残差分别为3.18×10－3和7.21×10－3，均方根相对误差最后为0.34%和0.87%.从图 4可以看出，在反演结束时，无论是正反演模型响应的相对误差、残差值还是均方根相对误差都已经很小，但是从反演结果上来看，和真实模型相比有一定的误差.此外，对于阻尼特征参数反演方法，不重要参数将会使解产生震荡现象，有时候甚至会出现不收敛的情况，图 4中HK型和KH型曲线的高阻层解的震荡比较严重，受两边低阻层的影响，反演结果很难达到真值.后续的研 究工作需要对不重要参数进行深入分析和研究，以降低不重要参数的对反演计算的影响，提高反演的稳定性和收敛速度.

图 4 HK型和KH型地电模型的反演结果
(a)和(b)表示模型反演结果；(c)和(d)表示反演模型响应；(e)和(f)表示均方根相对误差；(g)和(h)表示残差变化. Fig. 4 Inverse results of HK-type and KH-type The left and right column show the inverse results of HK-type and KH-type models respectively
(a) and (b)inverse results;(c) and (d)response of inverse model;(e) and (f)relative RMS error;(g) and (h)residual.
4 结 论

本文用阻尼特征参数算法对典型的三层和四层地电模型进行了反演计算，从反演结果看基本能和真实模型相吻合，所有模型的均方根相对误差都小于2%，反演模型响应和模拟数据的相对误差小于1%，这说明将阻尼因子和奇异值相结合，能更加有效地约束反演过程.此外，单点反演的耗时一般在10s左右，这也说明阻尼特征参数反演是一种快速的反演算法.当然，如果在反演时能给出合适的初始模型和参数约束条件，那么将极大地提高反演精度.

参考文献
[1] Abramowitz M, Stegun I A. 1964. Handbook of mathematical functions with Formulas, Graphs, and Mathematical tables [M]. New York: Dover Publications.
[2] Anderson H F, Duncan A C, Lynch S M. 1993. Geological mapping capabilities of the QUESTEM airborne electromagnetic system for mineral exploration-Mt. Isa Inlier, Queensland [J].Exploration Geophysics, 24(3-4): 333-340, doi: 10.1071/EG993333.
[3] Christensen N B. 2002. A generic 1-D imaging method for transient electromagnetic data[J]. Geophysics, 67(2): 438-447, doi: 10.1190/1.1468603.
[4] Christiansen A V, Christensen N B. 2003. A quantitative appraisal of airborne and ground-based transient electromagnetic (TEM) measurements in Denmark[J]. Geophysics, 68(2): 523-534, doi: 10.1190/1.1567220.
[5] Ellis R G. 1995. Airborne electromagnetic 3D modelling and inversion [J].Exploration Geophysics, 26(3): 138-143, doi: 10.1071/EG995138.
[6] Farquharson C G, Oldenburg D W. 1993. Inversion of time-domain electromagnetic data for a horizontally layered earth [J]. Geophys. J. Int., 114(3): 433-442, doi: 10.1111/j.1365-246X.1993.tb06977.x.
[7] Guptasarma D, Singh B. 1997. New digital linear filters for Hankel J0 and J1 transforms [J]. Geophysical Prospecting, 45(5): 745-762, doi: 10.1046/j.1365-2478.1997.500292.x.
[8] Guptasarma G. 1982. Optimization of short digital linear filters for increased accuracy[J]. Geophysical Prospecting, 1982, 30(4): 501-514, doi: 10.1111/j.1365-2478.1982.tb01320.x.
[9] Huang H, Palacky G J. 1991. Damped least-squares inversion of time-domain airborne EM data based on singular value decomposition [J]. Geophysical Prospecting, 39(6): 827-844, doi: 10.1111/j.1365-2478.1991.tb00346.x.
[10] Jupp D L B, Vozoff K. 1975. Stable iterative methods for the inversion of geophysical data [J].Geophys. J. R. Astr. Soc., 42(3): 957-976, doi: 10.1111/j.1365-246X.1975.tb06461.x.
[11] Knight J H, Raiche A P. 1982. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method [J]. Geophysics, 47(1): 47-50, doi: 10.1190/1.1441280.
[12] Li Y X, Qiang J K, Tang J T. 2010.A research on 1-D forward and inverse airborne transient electromagnetic method [J].Chinese Journal of Geophysics(in Chinese), 53(3): 751-759, doi: 10.3969/j.issn.0001-5733.2010.03.031.
[13] Liu G M, Asten M W. 1993. Fast approximate solutions of transient EM response to a target buried beneath a conductive overburden [J]. Geophysics, 58(6): 810-817, doi: 10.1190/1.1443466.
[14] Liu G, Asten M. 1993. Conductance-depth imaging of airborne TEM data [J]. Exploration Geophysics, 24(4): 655-662, doi: 10.1071/EG993655.
[15] Luo Y Z, Zhang S Y, Wang W P. 2003.A research on onedimensional forward for aerial electromagneticmethod in time domain [J].Chinese Journal of Geophysics(in Chinese), 46(5): 719-724.
[16] Macnae J C, Smith R, Polzer B D, et al. 1991. Conductivity-depth imaging of airborne electromagnetic Step-response data [J]. Geophysics, 56(1): 102-114, doi: 10.1190/1.1442945.
[17] Mao L F, Wang X B, Chen B. 2011.Study on an adaptive regularized 1D inversion method of helicopter TEM data[J]. Progress in Geophys. (in Chinese), 26(1): 300-305, doi: 10.3969/j.issn.1004-2903.2011.01.035.
[18] Raiche A P, Jupp D L B, Rutter H, et al. 1985. The joint use of coincident loop transient electromagnetic and Schlumberger sounding to resolve layered structures[J]. Geophysics, 50(10): 1618-1627, doi: 10.1190/1.1441851.
[19] Sattel D, Kgotlhang L. 2004.Groundwater exploration with AEM in the Boteti area, Botswana [J].Exploration Geophysics, 35(2): 147-156, doi: 10.1071/EG04147.
[20] Stolz E M, Macnae J C. 1997. Fast approximate inversion of TEM data [J]. Exploration Geophysics, 28(3): 317-322, doi: 10.1071/EG997317.
[21] Vozoff K, Jupp D L B. 1975.Joint inversion of geophysical data [J].Geophys. J. Int., 42(3): 977-991, doi: 10.1111/j.1365-246X.1975.tb06462.x.
[22] Wilson G A, Raiche A P, Sugeng F. 2006. 2.5D inversion of airborne electromagnetic data [J]. Exploration Geophysics, 37(4): 363-371, doi: 10.1071/EG06363.
[23] Wolfgram P, Karlik G. 1995. Conductivity-depth transform of GEOTEM data [J]. Exploration Geophysics, 26(3): 179-185, doi: 10.1071/EG995179.
[24] Worrall L, Munday T, Green A. 1998. Beyond bump finding-airborne electromagnetics for mineral exploration in regolith dominated terrains [J]. Exploration Geophysics, 29(2): 199-203, doi: 10.1071/EG998199.
[25] Zhu K G, linJ, Han Y H, etal. 2010.Research on conductivity depth imaging of time domain helicopter-borne electromagnetic data based on neural network [J].Chinese Journal of Geophysics(in Chinese), 53(3): 743-750, doi: 10.3969/j.issn.0001-5733.2010.03.030.
[26] 李永兴, 强建科, 汤井田. 2010. 航空瞬变电磁法一维正反演研究[J]. 地球物理学报, 53(3): 751-759, doi: 10.3969/j.issn.0001-5733.2010.03.031.
[27] 罗延钟, 张胜业, 王卫平. 2003. 时间域航空电磁法一维正演研究[J]. 地球物理学报, 46(5): 719-724.
[28] 毛立峰, 王绪本, 陈斌. 2011. 直升机航空瞬变电磁自适应正则化一维反演方法研究[J]. 地球物理学进展, 26(1): 300-305, doi: 10.3969/j.issn.10042903.2011.01.035.
[29] 朱凯光, 林君, 韩悦慧,等. 2010. 基于神经网络的时间域直升机电磁数据电导率深度成像[J]. 地球物理学报, 53(3): 743-750, doi: 10.3969/j.issn.0001-5733.2010.03.030.