（地球探测与信息技术专业论文）航空瞬变电磁法一维正演模拟与反演解释.pdf

摘要 航空瞬变电磁法( a t e m ) 具有效率高，成本低等优点，可以广 泛应用于地质填图、直接或间接找矿、水工环勘查等各个领域。近些 年随着国家的重视，航空瞬变电磁法得到了迅速的发展。 本文主要研究了航空瞬变电磁法的一维正演和一维反演算法。首 先从层状半空间上的磁偶极子电磁场公式出发，推导出了偶极一偶极 装置的一维正演公式，并应用数学计算领域中一些最新的研究成果对 算法进行了改进，编写了正演程序。通过对不同装置类型和不同参数 的层状大地模型的理论响应进行对比计算，分析了航空瞬变电磁法探 测的特点。 论文系统研究了现有的a t e m 数据反演方法，将现有方法分为近 似反演方法和定量反演方法两大类。对近似反演方法，本文实现了两 种算法：用迭代法求取全区视电阻率法和其中一种电导率快速转换方 法浮动薄板法。对定量反演方法，实现了具有光滑模型约束能力 的奥康姆反演法和具有较高反演精度的模型交替调整反演法。 通过上述研究工作，得出以下几点结论： 1 ) 由正演算例表明，a t e m 水平共面装置的探测信号最强，因 此在实际探测中应用得最多；a t e m 对低阻异常体的探测能力要远远 强于对高阻异常体的探测能力。 2 ) 近似反演方法因其速度快，是解决现场a t e m 海量数据处理 最好的方法，而定量反演方法效果好，但其速度较慢，只适用于室内 后期研究之用。 3 ) 对近似反演方法，浮动薄板法计算速度快，反映低阻异常的 效果明显，电阻率下降最快的深度可以认为为低阻层的上顶界面。 4 ) 在定量反演方法中，能获得光滑模型的奥康姆反演方法适合 于电磁法的特点，其稳定性较高；模型交替调整反演算法是一种技巧 性方法，模型算例证明了其效果较好，精度较高。 关键词航空瞬变电磁法，一维正演，全区视电阻率，浮动薄板法， 奥康姆反演法，模型交替调整反演法 a b s t r a c t a l r b o m et r a n s i e n te l e c t r o m a 印e l i c ( a t e m ) ，w 汕i t s h i 曲e 伍c i e n c ya n dl o wc o s t ， c o u l db ew i d e l y 印p l i e di nf i e l d ss u c ha s g e o l o 西c a lm a p p i n g ，d i r e c to ri n d i r e c t m l n e r a lp r o s p e c t i n g ，h y d r 0 1 0 百c a li n v e s t i g a t i o n ，e n v i r o l l m e l l t a le x p l o r a t i o na n ds oo n w i t ht h en a t i o n a la t t e n t i o ni nr e c e n ty e a r s ，a t e mh a sb e e i ld e v e l o p e dr a p i d l y i nt h l sp 印e r t h el da t e m 南n a r dm o d e l i n ga i l d i m e 单r e t a t i o nm e t h o da r e s t l l d i e d f i r s t l y ，t h el df o n a r dm o d e l i n gf o 珊u l ao fd i p o l e d i p o l ea n a yi sd 硎v e d 舶mt h ee l e c t r o m a 印e t i cf i e l d e x p r e s s i o no fam a 印e t i cd i p o l es o u r c ea b o v em e l a y 即e de a r c h ，a n dt h ea l g o t h mi so p t i m i z e dw i t hn e wa c h i e v e m e n t si nm a t h e m a t i c s i h ep r 0 黟锄i sm a d e u p ，a n dt h ef e a t l l r e so fa t e me x p l o r a t i o nm e t h o da r es t l l d i 酣b v a i l a l y s i n gt h er e s p o n s e so fl a y e r e de a 吡 p a r 姗e t e r sa n dm o d e lp a r 锄e t e r s m o d e l sw i t l ld i f r 醯e n ta m yt y p e s ，a r r a y a l l t h ea t e md a t ai n t 唧r e t a t i o nm e t h o d sa r e 咖d i e d s y s t e n l a t i c a l l v 锄d c l a s s i f i e di n t o a p p r o x i m a t i o ni n t e 叩r e t a t i o nt e c h l l i q u e sa i l d q u a i l t i t a t i v ei n v e r s i o n m e m o d s f o rt h ef o m t h i sp a p e rp r e s e n t s 铆ot e c h n i q u e s ： a l l t i m e a p p a r e l l t r e s i s t i v i t ya _ b t a i n e db yi t e r a t i o nm e t h o da 1 1 do n eo ft h e c o n d u c t i v i t y - d e p t h 仃a n s f o n n a t i o n ( c d t ) m e t h o d s ：s - i n v e r s i o n f o rt h el a t t t h e r ea r ea l s ot 、) l ，om e t h o d s p r e s e n t e d ：m eo c c 锄i n v e r s i o nm e t h o db yw h i c has m o o t h e s tm o d e lc o u l db e g e n e r a t e d ，a 1 1 dt h ea l t e m a t em o d e la d j u s t m e i l tm e t h o dw h o s er e s u l t sc o u l di n t e r p r e t t h ed a t aw i t hh i 曲p r e c i s i o n b yt h es t u d i e sa b o v e ，t h ef o l l o w i n gc o n c l u s i o n sc a i lb ea c q u i r e d ： 1 ) t h ef o n v a r dm o d e l i n gs h o w st h a tt h es i 印a lo ft h eh o r i z o n t a lp a r a l l e la r r a yi s s t r o n g e rt h a l la n yo t h e ra 1 1 r a y ，a n dt h a ti sw h yt h i sk i n di sm o r ea d o p t i v ei na p p l i c a t i o n ； t h ea n o m a l yi n d u c e db yc o n d u c t a l l c ei sm u c hm o r ep o w e r 如lt h a n t i l eo n eb y r e s i s t a n c e 2 )a p p r o x i m a t i o ni n t e 叩r e t a t i o n t e d m i q u e s a r cs u i t a b l e f o ri m m e d i a t e i n t e 叩r e t a t i o n f o rm a s sf i e l dd a t am l et oi t s f 弧t p r o c e s s i n gs p e e d ，w h i l e ，m e q u a l l t i t a t i v ei n v e r s i o nm e t h o d sw h i c h 百v eb 甜e rr e s u l t sw i ml o w e re m c i e n c yc a l lb e u s e di nl a t e rd a t ap r o c e s s i n ga n d a n a l y s i s 3 ) t h er e s u l t so ft h es - i n v e r s i o nm e t h o dr e n e c tt h ec o n d u c t a n c ea i l o m a l v o b v i o u s l y ，a n di nw h o s ep r o f i l e ，t h ed 印t hw h e r er e s i s t i v i t yd r o p sm o s ts h 唧l yc 锄b e c o n s i d e r e da st h et o po fc o n d u c t a n c e 1 a y e r 4 ) o c c 锄i n v e r s i o nm e t h o d ，w h i c hg e n e r a t eas m o o t h e s tm o d e l ，i ss u i t a b l e e l e c t r o m a g n e t i ce x p l o r a t i o n ，锄di t ss t a b i l i t yi sd 锄o n s t r a t e db ym o d e le x p e r i m e n t s ； t h ea l t 锄a t em o d e la 由u s t n l e n tm e t h o di sd e v e l o p e db ye x p e r i e n t i a ls k i n ，a l l dt h e e x 锄p l e sp r o v ei t sh i 曲e rp r e c i s i o na 1 1 db e t t e rr e s u l t s k e yw o r d s ：a i r b o m e 仃习m s i e i l te l e c t r o m a 莎e t i c ( a t e m ) ，1 - df o n v a r dm o d e l i n g ， a l lt i m e 叩p a r e n tr e s i s t i v i t y ，n o a t i n gs h e e tm e t h o d ，o c c 锄i n v e r s i o n ，a l t e m a t em o d e l a d j u s t l l l e mm e t l l o d l i l 中南大学硕+ 学位论文第一章绪论 第一章绪论 航空瞬变电磁法( a t e m ，a i r b o m et r a l l s i e i l te l e c t r o m a 妒e t i c ) 或称时问域航 空电磁法( a i r b o m et i m e d o m a i ne l e c t i o m a 鲫e t i c ) ，具有速度快，成本低，效果 好的特点，能在地面电磁法难以进入的山区、沙漠、丛林、湖泊等地区开展工作。 这种方法在地面上空间发射脉冲电磁波，观测由脉冲电磁场激发的二次电磁场， 一般为感应电动势随时间的衰减曲线，通过分析观测响应，可获知地下介质的电 性分布情况。航空瞬变电磁法广泛应用于地质填图、直接或间接找矿、水文地质、 工程地质勘查、环境监测等各个领域。 1 1a e r ot e m 正反演计算的发展与研究现状 国际上对航空瞬变电磁数值模拟的研究开展较早，2 0 世纪6 0 年代就有学者 研究了时间域的理论响应，并进行了相应的仿真研究，到9 0 年代已经实现了二、 三维的正演模拟【l 3 1 。如1 9 9 2 年x i o n 毋l 】利用积分方程法实现了航空瞬变电磁法 的三维j 下演模拟，同年l e p p i n 【2 】也研究了a t e m 的三维正演算法，他所用的方法 是有限差分法。二三维的正演模拟是二三维反演的基础，也为当今较常用的拟断 面成像方法研究提供了试验数据，他们的研究成果有很重要的意义。 国内航空瞬变电磁法的发展较早，但由于各种条件的限制发展缓慢。近些年， 随着国家对航空物探的重视，航空电磁法的理论研究和仪器开发研制等多项项目 己被列入国家重大科研项目中，例如国家8 6 3 计划“航空地球物理勘查技术系 统”、“航空瞬变电磁法快速正反演研究专题”等，航空瞬变电磁法的理论研究得 到了较快的发展。 2 0 0 3 年，罗延钟等【1 6 】研究了时间域航空电磁法的维正演算法，探讨了电 磁响应衰变特性，吊舱高度，大地导电性和收发距对响应的影响，为我国a t e m 的理论研究迈出了重要的一步。此f j ，国内地面瞬变地磁法的理论研究发展较快， 王华军等【4 】从矩形中心回线入手导出了2 5 维异常场有限单元算法，并给出时频 域转换关键技术。熊彬等【5 ，6 】对傅里叶逆变换和拉普拉斯逆变换做了详细分析， 进一步开发出基于独立电磁场的2 5 维瞬变电磁响应数值模拟算法，大幅削减了 求解稀疏线性方程组的规模。地面瞬变电磁法与航空瞬变电磁法有不同之处，但 他们的研究成果对a t e m 的正演研究具有重要的借鉴意义。 国际上，航空瞬变电磁法已经是一种普遍使用的勘查方法，实际应用的需要 对数据解释方法提出的要求不断提高，促进了a t e m 数据解释方法的发展。最 中南人学硕十学位论文第一章绪论 早使用的a t e m 数据处理方法是，直接将多个测点的感应电动势对时间的关系 曲线拼凑成拟断面图，直接依靠经验观察感应电动势相对相邻测点的关系判断高 低异常，这种“数据解释”方法，其实是迫于理论与技术的限制而采用的方法， 早已被淡忘。2 0 世纪六、七十年代，计算机应用还没有普及的时候，学者们研 究出了一些简单的近似反演方法【2 2 ，2 5 3 6 】由观测曲线导出的参数进行解释，如 视电阻率解释法，视纵向电导及视探测深度，视时间常数及视综合参数等，这些 参数一般可以通过对观测数据进行简单的运算得到，由于简单方便，至今很多场 合仍然使用( 详见第三章第一、二节) 。 随着计算机应用技术的发展，后来出现了很多数据解释的方法，目前国际上 普遍应用的a t e m 一维反演方法可以分为两类：一类是电导率深度转换 ( c d t ，c o n d u c t i v i t y d 印t ht r a i l s f o m a t i o n s ) 的近似方法，另一类是层状大地反 演( l e i ，l a y e r e d e a r c hi n v e r s i o n ) 。c d t 法包括f u l l a g a r 和r e i d 【6 1 】提出的e m a x a i r 法、m a c n a e 掣6 2 3 提出的e m f l o w 法、z h d a l l o v 等【4 6 4 7 1 提出的s 反演法，还有 c “s t e n s e n 【6 4 1 ，e a t o n 【6 3 1 ，l i u 和a s t e l l 【州，w b l f 西锄和k a r l i k 【删也提出了各自的 反演方法。以上的c d t 法还可以分为两种体系，一种源予n a b i 曲i a n 的“烟圈 理论，另一源于电导薄板模型理论( 详见第三章第三节) 。l e i 法即用传统的反 演理论拟合观测数据，有h u a n g 和p a l a c k y 6 9 1 ，s a t t e l 【7 0 1 ，c o n s t a b l e 等【7 2 1 ，和 f a r q u h a r s o n 等【7 3 】提出的算法( 详见第四章第一节) 。另外，s a t t e l 【8 0 】引入了z o h d y 法，该方法能取得比较好的效果，而且比一般的l e i 法有速度的优势。 近年，国际上的研究更是向高精度，智能化和二、三维反演的方向发展。2 0 0 4 年，s a s a l ( i 等1 7 j 提出了一种考虑地形，高度的反演方法，这种方法解决勘查过程 中高度测量不准确的问题以及地形的影响。2 0 0 6 年，b r o d i e 等【8 】研究出了一种三 维反演方法，可以在反演同时进行校正处理等，这种方法命名为“整体反演法 ( h o l i s t i c a p p r o a c h ) ，刚开始仅应用于频率域航空电磁法，其后有学者引入到了 时间域电磁法中。w | o l 琦锄等1 9 j 研究了拟二维反演方法，该方法基于此前 c s t e n s e n 【1 0 l 在波恩近似方法方面的研究工作，在现阶段，拟二维反演方法有可 能发展成为可以推广应用的方法，这项研究是非常有意义的尝试。另外，还有不 少学者也提出了各自的三维反演算法叫引。 国内，航空瞬变电磁法数据解释方法的研究，也随着国家重大项目的启动快 速发展起来。早期国内的研究为我们提供宝贵的经验，当时国内已经掌握简单的 参数解释方法，现在的研究将沿着前人的轨迹继续往f i 发展。而地面t e m 数据 解释方法的研究经过很多学者的努力得到很大的发展，如严良俊等【5 0 】从薄板理 论出发，通过利用感应电动势及其一阶导数表达式作简单数学推导，给出了全区 视纵向电导与视深度的反演算法研究了中心回线瞬变电磁测深全区视纵向电导 2 中南人学硕七学1 = 7 ：论文第一章绪论 解释方法，对应用实例表明该方法可以有效划分电性层；薛国强等【5 ，考虑了 隧道的特定环境，引入浮动薄板法，应用二次电导微分参数成像作隧道超前预报， 应用实例表明，该方法对隧道前方水体病害效果明显。我们将从这些地面t e m 数据解释方法的研究中得到启发。 总体来说，我国航空瞬变电磁法正反演的研究水平跟国际上的先进水平还有 很大的差距。我国a t e m 正反演研究的发展路线应该是先学习国外的先进技术， 然后吸收消化，再根据应用实际寻求改善。 1 2a e r ot e m 一维正反演研究的目的与意义 首先，航空瞬变电磁法具有速度快，成本低，效果好，能复杂克服地形情况 的优点，跟地面瞬变电磁法相比具有明显的优势和发展前景，在国际上已经是一 种普遍使用的勘查方法，而我国的发展水平较为落后。另一方面，我国经济的飞 速发展，资源相对紧缺，对资源勘查提出了更高的要求，而且我国要加快西部大 开发的进程，西部的条件恶劣，地形条件复杂，这些问题也要求我们发展更加可 行而高效的勘查方法。因此，现实情况需要我们开展航空瞬变电磁法的应用研究。 一维正演模拟是航空瞬变电磁法理论研究的重要手段，a t e m 一维正演模拟 可以为a t e m 的应用研究，如相关仪器和装置的设计，观测信号特征等研究提 供理论依据。一维正演模拟也是a t e m 一维反演解释的基础，它既是一维反演 迭代计算的必要过程，也为一维反演解释算法的研究提供试验数据。因此，要开 展a t e m 应用和解释方法的研究，就必须先实现a t e m 的一维正演模拟，尽可 能提高一维j 下演模拟的精度和速度。 国际上已经开展了a t e m 二、三维反演的研究，但由于a t e m 的海量观测 数据和计算条件的限制，目前应用最多的还是一维反演解释方法，包括c d t 法， l e l 法和z o h d v 法等。尽管这些方法现在已经相当成熟，应用普遍，但是深入研 究后，发现这些方法也不是十分理想，还有很多需要改进的地方。我国a t e m 数据反演解释方法的研究更是处于起步阶段，应该对现有的数据解释方法进行系 统研究，对有应用前景的方法进行改进，以满足a t e m 观测数据解释的需要。 1 3 本文的研究内容与方法 本文第二章将研究a t e m 一维正演的理论计算，首先根据层状半空间上的磁 偶极子公式推导出一维正演计算公式，并应用当前数学计算领域中一些最新的研 究成果对一维算法进行改进。编写了程序，通过对不同装置类型、装置参数和层 状大地模型的理论响应作比较分析，研究a t e m 的探测特点。 中南大学硕十学位论文第一章绪论 系统研究现有的a t e m 数据解释方法，将所有数据解释方法分为近似反演方 法和定量反演方法，近似反演方法的特点是速度快，定量解释方法的特点是给出 观测响应和模型理论响应之间的拟合差。本文把c d t 法和上述的参数解释方法 划入近似反演方法，把l e i 法和新提出的模型交替调整反演法归类为定量反演方 法。 第三章将研究两种近似反演方法，全区视电阻率法和浮动薄板解释法。其中 视电阻率法是早期的参数解释方法之一，文中将阐述视电阻率的定义，该解释方 法的优缺点，和已有的求取视电阻率方法，由i j f 人提出的迭代求取视电阻率的思 想，应用到a t e m 视电阻率计算中，就是本文的全区视电阻率法。 再系统地分析现有的c d t 法，根据其理论来源，分为两大体系，一种源于 n a b i 出a n 提出的“烟圈”理论，另一种源于电导薄板模型理论。后者就是浮动 薄板法及其变种。本文将从电导薄板模型的理论计算出发，推导出浮动薄板法的 计算公式，对模型理论响应进行计算验证其效果并分析其特点。 第四章将研究两种定量反演方法，奥康姆反演法和模型交替调整反演法。其 中奥康姆反演法属于l e i 法，本文首先探讨这种方法的思想，分析它的原理，论 证光滑模型的优点，把这种常用于地球物理反演的方法应用到a t e m 数据解释 中，通过对模型算例进行分析验证其应用效果。 本文把其后出现的z o h d v 法，和新提出的模型交替调整反演法根据其特点， 定义为一类技巧性方法，归类为定量反演方法。本文将通过实例，详细说明模型 交替调整反演法的计算步骤，分析其获得高精度的原理，最后通过模型实验算例 验证其效果。 4 中南人学硕十学位论文 第二章a e r o t e m 一维正演模拟 第二章a e r o t e m 一维正演模拟 现存的a e r o t e m 装置可以分为两种，一种是发射线圈与接收线圈分开有一 定距离，可以看作是两个磁偶极子，这种观测系统可以简化为偶极偶极装置模 型；另一种是采用大回线源发射，中心线圈或者重合回线接收，而这种观测系统 要作为回线源装置进行处理。罗延钟等【1 6 】推导出了层状大地条件下偶极偶极装 置的正演计算公式和算法，并做了相关实验，对a e r o t e m 的探测能力和探测条 件作了研究。刘桂芬【1 。7 】在其硕士期间实现了a e r o t e m 回线源装置的一维正演计 算，并探讨了发射边长和发射电流与电磁响应之间的关系。 本文研究一维正演，选用的装置是偶极偶极装置。当发射线圈面积不是很大 时，可以看作是一个磁偶极子，罗延钟等所做的研究得出平行装置( 当发射线圈 和接收线圈均为水平或者垂直时) 适合用小的收发距观测( k 1 5 m ) ，当收发距很小 时，可以看成重合回线装置，因此本课题选择偶极一偶极装置进行研究。 图2 一l 中，图2 1 a 为a e r o t e m 观测系统的示意图，图2 1 b 为层状大地上空 该观测系统的简化磁偶极子模型，t x 和r x 分别是发射和接收线圈，由发射线 圈提供一个瞬变的信号，在接收线圈观测由大地产生的响应值；五和z 分别是发 射和接收线圈的高度；r 为收发距；设大地为甩层水平层状模型，各层的电阻率 和厚度为：力，d ，；力，出；岛，磊，磊一。设圆柱坐标系原点o 位于发送 线圈的正下方的地面上，z 轴垂直地面向上。 t 隐n s m i 髓e r r e c e l 玎 i l i e a r l t ls u r 眨t c e ( a ) i 乙 l ( t x 9 r x h i i z o r l p l ，d l p 2 d ： p n ，d n 寸 ( b ) 图2 1 a t e m 磁偶板子场的简化模型 为计算时间域电磁场，我们采用先在拉普拉斯变换域中作计算，再借助 g a v e r s t e h f e s t 概率变换算法做逆拉氏变换计算电磁场的瞬变过程。本文首先实 5 中南大学硕十学位论文 第二章a e r o t e m 一维正演模拟 现已有的算法，然后研究最新的数值计算方法成果，对已有算法进行改进。以下 第一节推导出一维正演的计算公式；第二，三节给出汉克尔变换与逆拉氏变换的 方法和新的研究成果；第四节给出一维正演的算例，并讨论a e r o t e m 的探测特 点与探测能力。 2 1 一维正演的理论计算 由w a r d 和h o h m a n l l m “，推导出的公式司有，在空中垂直磁偶檄源( 水半线圈) 所产生磁场的水平分量h ；和垂直分量h ；分别为， 研= 罢j c o p 1 “外舢一h 沙。山谚州办) 以 ( 2 - 1 ) 彰= 昙f f 扛“) + 镌水“) 】丢厶( 办) 以 ( 2 - 2 ) 水平磁偶极源( 垂直线圈) 所产生的磁场的水平分量日；哿垂直分量日；分别为， 形= 一昙f 旷。一镌p “】【名山( 一手以( 刎m ( 2 - 3 ) 日；2 昙f k l “z + 舯一镌p o 。_ 妒以( 办) 掀 ( 2 4 ) 式中聊= ，为发射线圈的磁偶极矩( 单位为a m 2 ) ，为发射线圈的有效面积， ，为发射线圈中的供电电流强度( 单位为a ) ，厶和分别是零阶和一阶贝塞尔函 数，a 为积分变量，为反射系数 九一”口 任= i 1 九+ “q ( 2 5 ) “q 可以由以下递推公式求出， ，2瑞其中，：疗一。，以一2，；甜。：“。26) l f f = 历彳，砰叫掣。，圳1 ， 文献中z 轴方向向下，对以上公式在z 前加负号，用拉氏变换量s 替换f 国， 并在右端除以s ，水平线圈发射，在接收线圈处磁场的径向分量和垂直分量为， 日；( s ) 2 景f 【扩( 2 i 舯一e ”。扛“谚以( 办) 烈 ( 2 7 ) 州= 杀p 矿1 署舭 ( 2 8 ) 6 中南人学硕士学位论文第二章a e r o t e m 一维正演模拟 垂直线圈发射时，径向分量和垂直分量为， 啪) = 嚣肛舯一咿州圳私山( ”手似m ( 2 9 ) 其中， 彰( s ) 2 杀r 【们一饪p 1 以妒( 办) 以 ( 2 1 。) ， 一五一“。 2 蔬 l 中“一 ( 2 1 1 ) 甜。由以下递推公式求出， ，2咋描其中，：胛一。，疗一2，。；“：“。2。2， 吩= 肛彳，砰= 一眠肛，圳，1 扩1 蹦g h t 和r a i c h e 2 0 1 认为以p 也，d 一代替( 2 1 2 ) 式中的t h ( “，z ) ， 数值的稳定性，由此，可用下列递推公式代替， 2 篇 能更好地获得 ( 2 1 3 ) 试中e 由下列递推公式算出， 降也一鬻，向- 2 ，f p 也川“k ( 2 - 。4 ) 1 1 + r f + 1 f = 力一2 ，l 一3 ，1 只一l = p 叫一4 川见一l ( ，1 4 、 尺二：啦 “，+ ，f = 0 ，l ，z 实测的感应电动势与磁场分量有以下关系， y 一警一盹警d tn t ( 2 1 5 ) ( 2 一1 6 ) 其中s 。为接收线圈的面积，b 为磁感应强度分量。 由( 2 7 ) ( 2 1 0 ) 式可知，为计算感应电动势y 需要计算包含零阶和一阶贝塞 尔函数的无限积分( 汉克尔变换) ，为方便计算，引入归一化的无量纲参数， x = 打z r = z | r hr = h | r q 2 p 庐k p 、戈、h 戳。d i | r 。 墨5 岛房，s2 吩允= ( 1 + k q ) 2 ，( f = l ，2 “一，z )( 2 1 7 ) 由此可以写出感应电动势在拉氏变换域中的表达式为 7 中南大学硕十学位论文第二章a e r o t e m 一维正演模拟 哪) - - 警胪咐) 一谚川内( x 膨 ( 2 - 1 8 ) 鬈( s ) = 丛笔严n 州) + e - ( w 引2 厶( x ) 捌 ( 2 1 9 ) 髟( s ) = 一学f 乙弘朋+ 镌e _ ( n 砌x 2 山( x ) 一弛( x ) 】蚜 ( 2 2 。) 形( s ) = 一警f p ( 卜扪一镌p - ( z | n 乩引2 以( x ) 批 ( 2 2 1 ) ( 2 - 1 3 ) 式中的r 和e 按下列公式计算， x 怒，向。棚，1 一。玎磷k ( 2 _ 2 2 ) l + r 曩+ 1 江万一2 ，厅一3 ，1 c l = p 1 一“胁一r if ，、 p 一墨二显! 。s + & l ，f = o ，l ，刀( 2 2 3 ) 以上( 2 1 8 ) ( 2 2 1 ) 式是在拉氏变换域下面的公式，对该公式作逆拉氏变换， 可算出时间域感应电动势的上阶跃响应，即充电过程；取负值就可得到感应电动 势的下阶跃响应，即放电过程。( 2 1 8 ) ( 2 2 1 ) 式中积分号内的第一项，是与逆拉 氏变换变量s 无关的常量，而常数的逆拉氏变换在零时刻以外皆为零，故在实际 计算时可删除该项；积分号内的第二项只有反射系数镌与s 有关，所以，实际 上只需对做逆拉氏变换。 ( 2 - 1 8 ) ( 2 2 1 ) 式中存在的另一个问题是，以( 办) 是贝塞尔函数，其中胛= o 时为零阶，l = l 时为一阶，含有贝塞尔函数的积分式称为汉克尔型积分，可通 过汉克尔变换求得该积分的值。 以下两节将分别讲拉氏变换与汉克尔变换，和一些相关文献中提出的新理论 与新方法，并据此作了相关的对比计算。 2 2 逆拉氏变换 目前有多种从频率域转换到时间域的方法，k n i 出和r a i c h e l 2 0 】提出了用 g a v * s t e h f e s t ( 简写为g s ) 逆拉普拉斯变换的方法，跟其它方法相比，这种方 法较简单，是纯实数计算，而且对s 值计算的个数少，计算速度较快。w i o o d 掣堪1 对逆拉氏变换算法进行了研究，提供了g s 变换程序，对常用的逆拉氏变 换算法作了对比计算，和p i e s s e n s ，t a l b o t 等算法比较，在大多数情况下都有更 高的精度。 g s 变换算法是应用于地球物理计算中，常用的逆拉氏变换算法。其原理可 中南人学硕十学位论文 第二章a e m t e m 一维止演模拟 参考文献 2 5 ，2 6 】，其计算方法如下：对给定的时间t ，感应电动势的值可由( 2 2 4 ) 式计算， m ，= 等薹以矿( 半历) ， 协2 4 ， 其中，z 是决定于计 k = ( 一1 ) 的j 下偶整数；是变换系数， ： f 动) ! 可按( 2 2 5 ) 式算出， f f 磊 ( 2 2 5 ) 式中，求和下限f l 是( 聊+ 1 ) 2 的整数部分。当然，在做正演计算时，为减少计 算量，应先算出一组的值在程序中直接引用。对甩的取值，现存不同的意见， 罗延钟等认为，z 取1 6 ，罗延钟和昌彦君【2 6 】，昌彦君和张桂割2 3 】认为甩一般取 1 2 ；w 6 0 d e l l 等f 1 8 】却认为以取8 或1 0 较合适，当咒的取值增大时，变换会变得 不稳定，结果抖动将会明显。针对此问题，现对一些验算公式，用不同以值进行 对比计算，研究其相对误差的变化规律。 t ist ls 图2 2 g s 变换中胛不同取值的误差比较 图2 2 a ，为一个比较算例：把尺s ) = a r c t a n ( 1 居) 代入式( 2 - 2 4 ) ，对玎= 8 ，1 0 ， 1 2 ，1 6 四种情况进行计算，取f 的范围为1 0 击到1 0 ，( a t e m 常用的时问探测 范围为1 0 5 到l o 2 s ) ，计算的结果与解析解只力= s i n ( f ) f 之间的进行比较，通过 式( 2 2 6 ) 计算出数值解和解析解之间的相对误差： i ( f ) = l l 1 0 0 ( 2 2 6 ) 式中，c 。一，为数值解，匕“为解析解。由图中可以看出，随行的值增大，相 对误差变大，但是误差变得稳定；虽然以的值小，误差小，但随f 增大( 在靠近 0 1 s 的区段) ，误差增大了几个数量级，比行= 1 6 的误差都要大。 9 数卜嘲 位椭 虮 ! 奴 算 肿 中南大学硕十学位论文第二章a e r o t e m 一维正演模拟 图2 2 b ，为另一个比较算例：取耶) = l “l 忉) ，解析解为尺力= 一，利用上 述方法，计算出对不同以取值的相对误差。从曲线可以得出，刀的取值小的时候， 相对误差不稳定，l = 8 ，1 0 时，想对误差随f 增大，迅速上升；但由变化趋势， 可以从图中推测在f 更小的区段( f l o 6 ) ，相对误差较小；而刀= 1 6 时，相对误 差跳动较大，此与w b o d e n 所得到的结论也是吻合的。 由以上两个比较算例，可以得出，对，z 的值越大，应用的计算区间越广，但 是与解析解之间的相对误差也较大，还对另外一些验算公式进行的试验，可以得 到相同的结论；通过不同刀取值的对比计算，笔者认为，综合考虑误差与适用区 间，刀= 1 0 ，1 2 为应用于a t e m 正演计算较为合适的取值。 2 3 汉克尔变换 汉克尔变换是在地球物理领域，常用于电阻率计算的数字滤波算法。其计算 方法如下：以( 办) 是贝塞尔函数，其中刀= 0 时为零阶，胛= l 时为一阶，含有贝 塞尔函数的积分式称为汉克尔型积分， 似) 亍j ：k ( ，) 以( 办) 掀 ( 2 2 7 ) 可以通过线性滤波的方式求出， m ) = k ( 以) 彬 ，= l ， 见：三lo 【n + ( h ) ，l 7 ， ，f = l ，2 ，万( 2 2 8 ) 式中，彬是滤波系数。 目前常用的滤波系数是由a n d e r s o n 【2 4 】引入，2 8 3 个系数和8 0 1 个系数。 g u p t a s a 彻a 和s i n 曲【1 9 】提出一组新的汉克尔变换滤波系数，其中零阶6 1 或1 2 0 个系数，一阶4 7 或1 4 0 个系数( 系数较多的一组，精度要高，但是计算速度慢) ， 其系数比a n d e r s o n 的数值滤波算法少，计算速度快，但g u p t a s a 锄a 作了对比计 算，其精度却比a n d e r s o n 系数的高。因此，这组汉克尔线性滤波系数具有非常 高效的特点。 上文中，由于式( 2 1 8 ) ( 2 2 1 ) 进行了参数归一化以后，汉克尔变换的变量， 的取值总为l ，因此，对g u p t a s a 珊a 所做的对比计算结果，只需要考查1 0 9 l o ( ，) = o 点位置的计算精度，从文献 1 9 的结果中，大部分验算公式可以满足上述要求。 所以，采用g u p t a s a h n a 的汉克尔变换系数，在理论上可以提高a t e m 正演的精 度。 l o 中南大学硕十学位论文 第- 二章a e r o t e m 一维正演模拟 2 4 发送波形的影响 通过上述两种数值滤波方法解决了式( 2 一1 8 ) ( 2 2 1 ) 中存在的两个问题，可以 编写程序算出a t e m 的阶跃响应，但还必须考虑发射电流波形的问题。现考虑 发送波形为占空比为1 的正负方波，设上述公式算出的阶跃响应为1 ，( f ) ，则j 下负 方波的响应y 即( f ) 为： 4 _ 1 ，护( f ) = ( 一旷删“m 5 5 1 h ( f 一1 ) 乙】 式中，岛为正负方波的脉宽；i i l t 为取整函数；为某一足够大的正整数，程序 编写中取= 2 。 2 5 一维正演的算例与分析 根据上述方法及所考虑到的问题，可以用f o r t m 语言编写出正演程序。为 研究a t e m 的装置对观测曲线的影响，和a t e m 的探测能力，现作以下实验， 并对实验结果展开分析与讨论。实验中，只计算接收线圈的有效灵敏度d 剧d f 的 时间响应，即发射磁矩为l a m 2 ，接收线圈等效面积为1 m 2 时的感应电动势。有 以下三个实验：( 1 ) 首先更换装置类型，探讨不同装置类型对同一模型的响应的 特点；( 2 ) 然后更换装置参数，找出不同参数对响应影响的规律；( 3 ) 最后计算 不同大地模型的响应，和均匀大地模型响应之间的相对异常。 以下进行的计算，不失一般性，这里选取某些参数： 发射线圈的有效面积：品v = l 聊2 ； 发射电流：j = 1 彳； 接受线圈的有效面积：& v = 1 聊2 ： 供电波形为占空比为1 的正负方波；供电方波脉宽为1 0m s ，采样点数为1 4 ； 发射线圈高度和接收线圈高度取3 0 5 0 m ；收发距r 取o 1 5 0 m 。 装置类型编号用n a 表示： a = l ，水平线圈发射，水平线圈接收； a = 2 ，垂直线圈发射，垂直线圈接收； a = 3 ，水平线圈发射，垂直线圈接收； a = 4 ，垂直线圈发射，水平线圈接收。 2 5 1 不同装置类型的响应比较 先给出一个三层地电断面模型：力= 3 ，助= 2 0 ，d = 3 q m ；刃= l o o ，玉= 3 0 0 m ， 以一。装置参数为，矗= z = 3 0 m ，= 2 5 m ，利用4 种探测装置对该模型理论响 中南人学硕十学位论文第二章a e r o t e m 一维正演模拟 应曲线如图2 3 a 所示，可以看出，虽然不同装置测得的感应电动势( 或磁感应强 度对时间的变化率) 在数值上各不相同；但它们随时间变化的性态( 瞬变响应曲 线的形态) 基本上是一致的，即随时间增大而衰减；在小收发距的情况下，感应 电动势的强度按第一种装置( a = 1 ) ，第二种装置( a = 2 ) ，第三和第四种( a = 3 ， 4 ) 装置依次减弱；前两种装置的衰减速度相同，后两种装置的衰减速度较快。应 该指出，第三和第四种装置的响应完全相同，实际上它们的公式( 2 1 8 ) 和( 2 2 1 ) 完全相同，因而两种计算结果自然一样，即两种装置是完全等效的。 2 4 2 不同装置参数的响应比较 装置参数共有， ，z ，j ，踟，地六个，从公式( 2 1 8 ) ( 2 2 1 ) 可以看出j ， s 叔，锨三个参数与瞬变响应的幅值成正比关系，因此这里讨论飞机高度与收发 距对瞬变响应的影响。 图2 3 a 是对发射线圈高度办= 3 0 m 和接收线圈高度z = 3 0 m 算得的，对| l = 4 0 m ，z = 2 0 m 和j i = 2 0 m ，z = 4 0 m 的情况作了计算，其结果和图2 3 完全相同。 这可由瞬变响应的公式( 2 1 8 卜( 2 2 1 ) 得到解释，上文中已提到只需计算上述四个 公式中积分号内的第二项，而该项与发送高度和接收高度之和( 日一孙) 有关，即 a t e m 的响应只决定于发射线圈高度与接收线圈高度之和，而与单个取值无关。 这一结论是基于层状大地模型上的理论公式的，脱离层状大地模型，这一结论不 再成立。 图2 3 不同装置类型( a ) 和不同收发距( b ) 的响应比较 图2 3 b 为变换收发距，r = 0 1 ，1 ，5 ，8 ，1 5 ，3 0 ，5 0 m 等7 个值，对上述 模型，利用第二种类型装置( a = 2 ) ，在 = 3 0 m 和z = 3 0 m 条件下的响应曲线。 由图可知，当收发距很大时( 例如图中曲线l ，= 5 0 m ) ，瞬变响应的早期出现负 值( 虚线) ，瞬变曲线变得比较复杂，先负后正，最初由小变大，取得最大值后 1 2 中南人学硕+ 学位论文第二章a e r o t e m 一维正演模拟 又由大变小，最后衰减至零；随着收发距减小，感应电动势和相对于均匀大地时 的异常值都变大；而当收发距r 1 5 m 时，已接近于，_ o 时的渐近值。对四 种装置类型的计算可以得出，第一种装置瞬变响应曲线随r 的变换规律与第二种 装置的相似，第一和第二种装置在此称为平行装置。第三种装置正交装置的 瞬变响应随收发距，的变换形态，与第一和第二种装置的完全不同，特别是在， 1 5 m 时，正交装置的感应电动势和绝对异常值与厂成正比减小，相对异常值 保持不变，且当，专0 时，j 下交装置的响应值和异常值都趋于零。 由以上实验可以得出，平行装置适用于小收发距，甚至同点装置；而正交装 置则适宜采用较大的收发距不宜采用同点装置。 1 0 0 术 = 一1 0 q 2 0 0 246 81 0 f ，瞧 图2 4 不同覆盖厚度的k 型模型( a ) 和h 型模型( b ) 的相对异常比较 ( a ) k 型模型的参数：朋= 3 ，见= 2 0 ，胁= 3 q m ，如= 3 0 0 m ；( b ) h 型模型的参数：n