2010, Vol. 53
航空瞬变电磁法（ATEM，Airborne Transient Electromagnetic）或称时间域航空电磁法（Airborne Time-domain Electromagnetic），具有速度快，成本低，能在地面电磁法难以进入的山区、沙漠、丛林、湖泊等地区开展工作.因此，为加速我国西部大开发，解决地下水资源和矿产勘查的需要，非常有必要开展对ATEM一维正反演的方法研究.
目前国际上研究和应用的ATEM反演方法可以分为两类; 一类是电导率-深度转换（CDT，Conductivity-Depth Transformations）的近似方法，另一类是层状大地反演（LEI，Layered-Earth Inversion）.CDT法包括Fullagar和Reid[1]提出的EmaxAir法、Macnae等[2]提出的EMFlow法、Zhdanov等[3]提出的S反演法，还有Christensen[4]，Eaton[5]，Liu和Asten[6]，Wolfgram和Karlik[7]也提出了各自的反演方法.CDT法反演速度快，但是有两个缺点，一是容易受到噪声的影响，二是没有给出响应比较和误差估计，因此不好对反演效果进行评价.LEI法直接用模型的响应来拟合观测数据，用最小二乘法建立一个目标函数，使这一函数的值达到最小.Huang和Palacky[8]，Sattel[9]对LEI方法做过研究，为了简化算法，他们提出的方法中，建立目标函数时只考虑了数据的误差项，这种方法需要解决一个超定问题，得到的反演模型一般比较凹凸，效果并不理想；Constable等[10]，和Farquharson等[11]研究了Occam反演法，目标函数引入了模型项，增加了平滑模型的限制条件，但是需要非常耗时的参数正则化计算，因此速度比较慢.另外，Sattel[12]把常用于直流电法的Zohdy法引入到ATEM反演计算中，该方法首先把响应曲线转换为电导率和深度的关系，再调整模型参数来拟合观测数据，使均方差达到最小，能够取得较好的效果，但是这种方法的反演深度过小，忽略了深部异常体的响应信号，当要求反演较深处的地质体时产生较大的误差.
基于以上方法存在的问题，本文提出一种模型交替调整反演算法，解决了普通LEI反演方法速度较慢的问题，又避免了模型深度的过分压缩，从而保留了大地较深部的响应信息，编写了相关的程序，对若干典型模型做了试算，取得了较理想的效果.
2 ATEM一维正演罗延钟等[13]对ATEM一维正演算法作了研究，采用的是先在拉普拉斯变换域（下称拉氏域）中作计算，然后借助Gaver-Stefest概率变换算法（G-S变换）做逆拉氏变换转换到时间域的方法.
Wooden等[14]对逆拉氏变换做过研究，提供了G-S变换程序，对目前常用的逆拉氏变换算法做了对比计算，和Piessens，talbot等算法相比，对大多验算公式都可以取得更高的精度.G-S变换的具体算法可以参考文献[15]和[16].对正演公式中的贝塞尔函数Jn（λr），可以通过汉克尔变换的方法求出含有贝塞尔函数的汉克尔型积分式.Guptasarma和Singh[17]提出了两组新的汉克尔变换滤波系数，其中零阶61或120个系数，一阶47或140个系数，其系数个数比Anderson的数值滤波算法的少，因此计算速度快，但是Guptasarma做了验证和对比计算，其精度比Anderson提出的系数要高.这两组新的汉克尔线性滤波系数具有非常高效的特点，因此，编写正演程序时采用新的汉克尔变换可以得到更高的精度.汉克尔变换的具体算法可参阅文献[16]和[18].
本文正演理论计算基于文献[13]，采用新的汉克尔变换系数后，提高了计算精度.以下给出两组计算算例.第一组计算算例; 取发送电流I=1A，发送线圈等效面积STX=接收线圈等效面积SRX=1 m2，发送线圈高度h=接收线圈高度z=30m，收发距r=2.5m等装置参数，对参数为ρ1=3，ρ2=20，ρ3=3Ωm；d1=100，d2=300m，d3→∞的层状大地模型，用4种装置做正演计算.NA为装置类型号（NA=1，2，3，4，分别代表水平（线圈）发送，水平（线圈）接收；垂直发送，垂直接收；水平发送，垂直接收；垂直发送，水平接收；），在10-2~10ms之间等对数间隔取14个点采样，得到4条响应曲线，如图 1a所示.由图 1a可以看出，4种装置响应的衰减形态都是一致的，NA=1类型装置的响应值最强，NA=2稍低，NA=3，4最弱；NA=3，4两种装置的曲线完全一样，由等效原理可以推知它们的结果是一样的，而实际上它们的公式是完全相同的；从仪器测量上讲，1型装置的信号最强，测得的信号也就最可靠，实际上也用得最多，因此以下另外一个实验和反演研究中都用1型装置.
|
图 1 一维正演程序计算算例 （a）不同装置类型的响应曲线；（b）不同埋深低阻异常体的响应相对异常，模型参数，ρ1=20，ρ2=3，ρ3=20Ωm；d2=300m. Fig. 1 Examples of forward program.(a) Response curves with different array; (b) Relative anomalies of the response from low resistivity anomality with differen covered thicknesses |
第二组算例为计算相对异常的实验，如图 1b所示; 装置参数不变，对H型模型做正演计算，模型参数在图说明中给出，变换第一层的厚度，计算出它们的响应值，并和电阻率为20Ωm的均匀大地正演响应做对比计算，利用下式算出各个模型响应的相对异常;
|
(1) |
式中，Δ表示相对异常值，B′mod表示H模型的响应值，B′ave表示电阻率为20 Ωm的均匀大地的响应值，下标表示第i道.
由图 1b可以看出，厚度为300m的低阻层埋藏在300m以下仍可以产生超过23%的相对异常，而且还有闭合的趋势，甚至到600m以下相对异常仍超过10%.以上计算的结果和文献[13]所得到的结果相符，在计算汉克尔变换时采用新的一组系数，理论上提高了计算精度，相对异常值要比文献中的结果大，因此我们认为ATEM的实际探测能力原来被低估了，ATEM可以探测得更深.
3 ATEM一维反演通过ATEM正演程序计算一些典型模型的正演响应，作为实测数据开展ATEM反演算法的研究.
ATEM模型交替调整反演算法的过程可以分为以下三个步骤; （1）根据装置类型和装置参数选择或者建立一个相应的插值函数表；（2）根据每一道的观测响应在插值函数表中通过插值法取得相应的视电阻率值，并计算深度值，形成一个初始模型；（3）最后对模型参数进行调整，包括深度调整和电阻率值的调整两个过程交替地进行.
3.1 建立插值函数表反演方法的效果很多情况下依赖于初始模型的选取，本文提出的模型交替调整反演方法能够根据ATEM正演响应曲线每一测道的幅值，自动选取初始模型，这种方法选取的初始模型可以大致反映地下电阻率的变化趋势，因此在反演计算过程中更有利于收敛，得到更小的误差.为此，必须先建立一组选取视电阻率的函数，称为插值函数表.为计算方便，反演中所采用的装置类型和装置参数与上文正演算例图 1b中所用到的一致.
这个函数表的建立方法如图 2所示，首先计算给定装置类型和装置参数条件下，不同电阻率值的均匀大地模型如10Ωm的正演响应曲线，如图 2a所示，为了程序的通用性，可以在电阻率为10-1.5~103.5Ωm之间等对数间隔（这个范围包含了常见的电阻率值），取若干个电阻率值，计算出它们的响应曲线；然后，整理出每一道中响应值和电阻率的关系，作为一个函数保存下来，如图 2b所示，channal-1曲线为第一道响应值和电阻率的关系曲线；正演算法中设计有14个测道，由此可作出14条曲线.这些函数以表格形式保存下来，在后面计算中以插值方式取得视电阻率值，因此称为插值函数表.
|
图 2 建立取视电阻率函数表方法 （a）电阻率为10Ωm的均匀大地模型的响应曲线；（b）各道响应和电阻率的关系曲线，图中所示从上往下14条曲线分别代表 1到14道的函数曲线顺序分布. Fig. 2 Construct the fucntions between response amplitudes and apparent conductivity (b) by responses of homogenous models (a).the 14 lines from the top to the bottom displayed in figure2b represent the functions of the 14 channals in sequence |
为方便解释获得初始模型的过程，本文给出了一个算例，如图 3所示，图 3a为一个H型模型的正演响应曲线，模型的参数为ρ1=100，ρ2=10，ρ3=100Ωm；d1=50，d2=100m，d3→∞.判断最后反演的效果，在于最终反演的模型与原模型的吻合程度，和最终模型的计算响应与原始响应之间的误差.
|
图 3 求初始模型的步骤（a）模型响应，模型参数ρ1=100，ρ2=10，ρ3=100Ωm；d1=50，d2=100m，d3→∞；（b）根据响应数据插值得到视电阻率的值；（c）根据趋肤深度公式（2）算出初始深度值；（d）由公式（5）算出初始模型的电阻率. Fig. 3 The process to access the first model (a) Response from the true model; (b) Apparent resistivity derived from the function-table; (c) Depth value derived with equation (2);(d) First model resistivity derived with equation (5). |
观察图 2b中早期测道的响应值和电阻率的关系，响应值随电阻率的增大有一个先升后降的过程；而晚期道则单调递减，这是由于发送波形的周期性所致.因此早期道的取值有不确定性，可能取到双值.为了解决这个问题，应当从晚期道开始，倒序进行取值（晚期道取值惟一），到早期道时，即使取到双值，可以由视电阻率变化的连续性进行判断，取一个与稍晚一道接近的值，以此可以消除早期道的不确定性.如图 3b就是用这种方法推导出的视电阻率对时间的关系曲线.
用取得的视电阻率值和该点时间可以由下式求出相应的深度值;
|
(2) |
如图 3c所示，就是视电阻率对深度的曲线.那么初始层厚度;
|
(3) |
假设上层视电阻率可以反映真实模型的电阻率，则可以通过下式推导出初始模型电阻率;
|
(4) |
也可以从下层往上推导;
|
(5) |
图 3d所示曲线就是利用（5）式计算出来的初始模型电阻率和深度的关系曲线.由于早期道插值函数曲线的形态不为单调，出现极值，而其反函数在极值附近的导数很大，因此插值所得的视电阻率值误差比较大，用（4）式从上层开始推导，容易使误差也逐层传递，以至整个初始模型误差都较大，因此建议使用（5）式.对比试验结果表明，应用（5）式推导的效果比（4）式的好.
3.3 模型参数调整模型参数调整可分为两部分内容; （1）电阻率不变，深度进行调整；（2）深度不变，电阻率进行调整.两个步骤交替地进行可以得到更好的反演结果.
3.3.1 模型参数调整的方法第一步，先固定电阻率不变，所有层厚度值均减小一定比例，比如说5%，10%等，然后计算压缩前后模型的正演响应，再计算响应值和观测响应之间的拟合均方根误差;
|
(6) |
式中，m为测道数，也是模型的层数，Yid是第i道的观测响应，Yic是由模型正演计算出的第i道响应；如果压缩处理后的δ值减小了，就重复这一步处理，直到均方根误差δ取得最小值；如果刚开始压缩模型就使δ增大，那么这一步就应该反向，尝试伸展模型，直到δ取得最小值.厚度压缩或者伸展的比例5%或10%称为模型伸缩因子.用上述方法对图 3d中的初始模型进行处理，可得到调整后的模型，如图 4a所示，虚线代表初始模型，实线代表经过深度调整后得到的模型.
|
图 4 模型调整 （a）电阻率不变，深度调整；（b）深度不变，电阻率调整 Fig. 4 Model after depth adjustment (a) and after resistivity adjustment (b) |
第二步，深度调整完了就进行电阻率的调整.调整过程可以用以下的方法迭代地进行处理，模型电阻率ρiold乘以观测响应的视电阻率ρiapp_observed和模型计算响应的视电阻率ρiapp_calculated的比值，得到新的模型电阻率ρinew，可用以下公式表示;
|
(7) |
各层电阻率值同时通过上式修正，然后计算新模型的正演响应Yic，再利用（6）式算出均方根误差δ，判断若δ减小了，就重复计算，直到δ取得最小值.如图 4b所示，粗实线是真实模型，细实线是最终反演模型，由图可以看出最终模型基本上可以反映真实模型的形态，能够识别低阻层的埋深和电阻率值，最后算出的均方根误差δ的值是4.81%.
3.3.2 交替调整方法以上结果是对深度和电阻率进行了一轮调整，要想得到更理想的效果，应该对模型进行连续交替的调整，其流程如图 5所示.首先给定深度调整的次数nd，电阻率调整的次数nr，和一个足够小的值Δ，如（2%，5%或者0）；求出初始模型后，进入上述的第一步厚度调整，调整了nd次以后，即使δ还没有达到最小值也退出.当然，如果此过程中nd已经达到最小，直接退出；进行上述第二步电阻率调整，同样的方法连续进行nr次调整，直到求得δ的值不能再减小或者nr次调整都完成后退出.这样一轮调整结束后，如果δ没有减小或者已经小于给定的Δ值，如果给出的Δ=0，则重复以上计算过程，直至使得δ不能再减小.在nd，nr和Δ都给得非常合理的情况下，只要一轮计算就可以使δ收敛至很好的值.
|
图 5 交替调整算法流程图 Fig. 5 Flow chart of alternate adjustment |
其原理是，初始模型由视电阻率推导得出，而视电阻率的值是由均匀大地模型的响应值取得，当存在覆盖层的影响时，视电阻率的变化就变得平缓而不明显；利用趋肤深度公式计算模型的深度，由于公式中就含有视电阻率的值，模型异常层的界面跟真实模型之间有差别，因此要对模型的厚度进行伸缩处理，使模型的电阻率开始变化的层与真实模型相对应，此时再进行电阻率的迭代处理就可以得到很好的反演效果.观察图 4a，当δ达到最小时，由于模型尾部的响应信号不明显，模型电阻率的极小值都压缩到了真实的低阻界面深度处，所以，使δ值最小这一判断标准并不太合理.然而，不可能知道真实模型的形态，因此，nd，nr和Δ等参数也不可能非常合理地给出，应该用交替调整的方法试探性地寻找收敛点.
如图 6所示就是对上述例子的最终反演结果，取伸缩因子为10%，nd=3，nr=30，Δ=0，经过了三轮调整，由图可以看出，反演曲线的形态更加贴近真实模型，其边界、极小值的位置和电阻率的值都更加精确，最后均方根误差δ达到1.22%.
|
图 6 用交替调整方法反演出的最终模型 Fig. 6 Inversed model (broken line) and true model (solidline) |
为了更好地评价反演的效果，下面对图 6中的最终反演模型求出正演响应和观测响应，即反演结果与真实模型的响应做比较，求出各测道幅值的相对误差.两条曲线如图 7所示，各测道在图中标记出，每个测道对应时间点上标有“+”的曲线表示正演响应，标有“°”的曲线表示最终模型的观测响应.由图可以看出，两条曲线几乎重合在一起，反演模型响应曲线跟观测响应曲线拟合得非常好.
|
图 7 观测响应与反演模型响应曲线对比 Fig. 7 Observed response (with crosses) and calculated response (with circles) |
表 1给出了模型的正演响应和观测响应各测道的幅值和它们之间的相对误差，观察表 1中数据可以看出，最后一道的误差值比较大，达到3.34%，其次是较靠后的4个点的误差在1~3%之间，其余各道的相对误差都在1%以下，甚至有的在0.1%以下.可以得出，模型交替调整反演方法可以得到比较高的精度.
|
表 1 模型响应与观测响应的相对误差 Table 1 Relative error of calculated response to observed response |
为了进一步验证模型交替调整反演方法的效果，下面给出一些典型的算例，包括H型，K型，四层HK型和KH型模型.
3.4.1 H型模型算例以上算例是一个H型模型，取得了比较好的效果，下面给出另一个H型模型，加大异常体和覆盖层的厚度，再尝试其反演效果.真实模型的参数为; ρ1=100，ρ2=10，ρ3=100Ωm；d1=100，d2=150m，d3→∞，如图 8a中粗线所示.
|
图 8
三层模型反演算例
（a）H型模型参数; ρ1=100，ρ2=10，ρ3=100Ωm；d1=100，d2=150m，d3→∞； （b）K型模型参数ρ1=10，ρ2=100，ρ3=10Ωm；d1=50，d2=100m，d3→∞； Fig. 8 New method and Zohdy's method compare with Hmodel (a) and Kmodel (b) |
取伸缩因子为10%，nd=3，nr=30，Δ=0，经过了三轮调整，可得到如图 8a中实线所示的反演模型.可以看出，曲线基本可以反映真实模型的形态，电阻率的峰值和数值都和真实模型很接近，均方根误差达到1.07%，本文还用Zohdy法作对比，Zohdy法的反演结果如图中虚线所示.由图 8a可以看出，当深度增大时，Zodhy法的反演效果并不理想.对模型交替调整反演方法反演的模型做了响应对比和误差分析，可以得到与上面算例同样的结论，两响应曲线基本重合，各道的相对误差很小.
3.4.2 K型模型算例K型模型的算例中真实模型的参数为; ρ1=10，ρ2=100，ρ3=10Ωm；d1=50，d2=100m，d3→∞，如图 8b中粗线所示.取伸缩因子为10%，nd=1，nr=30，Δ=0，可得到图 8b中实线所示的模型，虚线为Zohdy法计算出来的结果，可以看出实线更能反映真实模型的形态，效果明显优于后者，虽然它们最后的正演响应与原始响应之间均方根误差差不多.还对nd的取值做过一些实验，分别对nd=2，3，4，5都作了计算，发现对该模型，nd越小反演的效果越好.
本文还对覆盖层更厚一点的K型模型进行过计算，发现埋藏深度增大到100m时，反演的结果可以显示出高阻异常，但是模型尾部已经不闭合，跟H型模型比较起来，其反演深度不深.这也验证了正演研究得出的结论，ATEM对低阻的探测能力比高阻要强.
3.4.3 四层KH型和HK型模型算例本文还利用模型交替调整反演方法计算了四层KH型和HK型模型响应，模型参数和计算结果如图 9所示.
|
图 9 四层KH模型和HK模型反演算例 （a）KH型模型参数; ρ1=20，ρ2=100，ρ3=10，ρ4=100Ωm；d1=40，d2=80，d3=100m，d4→∞.（b）HK型模型参数; ρ1=100，ρ2=10，ρ3=100，ρ4=20Ωm；d1=40，d2=80，d3=100m，d4→∞. Fig. 9 New method applied for KH model (a) and HK model (b) |
由图 9可以看出，反演的结果在形态和数值上都和真实模型相吻合，也和Zohdy法的反演结果作了比较，得到和之前三层模型相同的结论，模型交替调整反演法更能反映原始模型的形态，可以得到更小的均方根误差，效果更好.
图 9a中为KH型模型，细实线表示反演模型，粗实线表示真实模型，最后的均方根误差δ达到0.16%，模型匹配的效果非常好.
图 9b中为HK型模型最后均方根误差δ为0.7%，从图 9b可看出，反演模型的电阻率比真实模型参数中第三层高阻体电阻率值低，跟真实模型不能很好地吻合，一般情况下反演高阻层的效果没有低阻的效果那么好，如图 8b和9b，但能够反映电阻率的变化趋势，再次验证了ATEM探测高阻没有探测低阻的效果好的结论.
4 结论（1）优化了ATEM一维正演算法，取得了更高的计算精度.
（2）提出了一种模型交替调整反演方法.一些典型模型算例表明，该方法能够取得更好的反演效果.
（3）模型交替调整算法属于一种人机交互反演方法，要求使用者有一定经验，对不同的模型给出合适的伸缩因子、深度调整次数nd、电阻率调整次数nr等.考虑到时间域航空电磁法的数据量非常大，研究表明，伸缩因子为10%，nr给定一个较大的值（如30），nd取2，一般情况下可以得到较为理想的反演结果.
（4）模型交替反演方法跟CDT类方法相比较，因为包含迭代拟合和判断的过程，对每一测点反演所需的时间肯定比CDT类方法多得多，但具有精度高的优势，适合应用于后期深入的数据处理.
| [1] | Fullagar P K, Reid J E.Emax conductivity-depth transformation of airborne TEM data.15th Conference and Exhibition, Australian Society of Exploration Geophysicists.Extended Abstracts, 2001 |
| [2] | Macnae J, King A, Stolz N, Oskamkoff A, Blaha A. Fast AEM processing and inversion. Exploration Geophysics , 1998, 29: 163-169. DOI:10.1071/EG998163 |
| [3] | Zhdanov M S, Pavlov D A, Ellis R G. Localized S-inversion of TDEM data. Geophysics , 2002, 67: 1115-1125. DOI:10.1190/1.1500372 |
| [4] | Cristensen N B. A generic 1-D imaging method for transient electromagnetic data. Geophysics , 2002, 67: 438-447. DOI:10.1190/1.1468603 |
| [5] | Eaton A P. Application of an improved technique for interpreting transient electromagnetic data. Exploration Geophysics , 1998, 29: 1175-183. |
| [6] | Liu G, Asten M. Conductance-depth imaging of airborne TEM data. Exploration Geophysics , 1993, 24: 655-662. DOI:10.1071/EG993655 |
| [7] | Wolfgram P, Karlik G. Conductivity-depth transform of GEOTEM data. Exploration Geophysics , 1995, 26: 179-185. DOI:10.1071/EG995179 |
| [8] | Huang H, Palacky G J. Damped least-squares inversion of time-domain airborne EM data based on singular value decomposition. Geophysical Prospecting , 1991, 39: 827-844. DOI:10.1111/gpr.1991.39.issue-6 |
| [9] | Sattel D. Condustivity information in three dimensions. Exploration Geophysics , 1998, 29: 157-162. DOI:10.1071/EG998157 |
| [10] | Constable S C, Parker R L, Contrable C F. Occam's inversion:A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics , 1987, 52: 289-300. DOI:10.1190/1.1442303 |
| [11] | Farquharson C G, Oldenburg D W, Routh P S. Simultaneous 1D inversion of loop-loop electromagnetic data for both magnetic susceptibility and electrical conductivity. Geophysics , 2003, 68: 1857-1869. DOI:10.1190/1.1635038 |
| [12] | Sattel D. Inverting airborne electromagnetic (AEM) data with Zohdy's Method. Geophysics , 2005, 80: 77-85. |
| [13] | 罗延钟, 张胜业, 王卫平. 时间域航空电磁法一维正演研究. 地球物理学报 , 2003, 46(5): 719–724. Luo Y Z, Zhang S Y, Wang W P. A research on one-dimention forward for aerial electromagnetic method in time domain. Chinese J.Geophys. (in Chinese) , 2003, 46(5): 719-724. |
| [14] | Wooden B, Azari M, Soliman M. Well test analysis benefits from new method of Laplace space inversion. Oil & Gas Journal , 1992, 90(29): 108-110. |
| [15] | Knight J H, Raiche A P. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method. Geophysics , 1982, 47(1): 47-50. DOI:10.1190/1.1441280 |
| [16] | 朴化荣. 电磁测深法原理. 北京: 地质出版社, 1990 : 139 -161. Piao H R. Principle of Electromagnetic Sounding Method (in Chinese). Beijing: Geological Publishing House, 1990 : 139 -161. |
| [17] | Guptasarma D, Singh B. New digital linear filters for Hankel J_0 and J_1 transforms. Geophysical Prospecting , 1997, 45: 745-762. DOI:10.1046/j.1365-2478.1997.500292.x |
| [18] | 方文藻, 李予国. 瞬变电磁测深法原理. 西安: 西北工业大学出版社, 1993 : 54 -58. Fang W Z, Li Y G. Principle of Transient Electromagnetic Sounding Method (in Chinese). Xi'an: Northeast College of Industry Press, 1993 : 54 -58. |