瞬变电磁一维正反演算法研究

学位是每个学校教学科研的重要学术成果，是反映学校学术特点和学术水平的文献资源 分类号 密 级 U D C 编 号 硕 士 学 位 论 文 名和副题名 中心回线瞬变电磁一维正反演算法研究 作 者 姓 名 郭嵩巍 指导教师姓名及职称 王绪本 教授 申请学位级别 硕士 专业名称 地球探测与信息技术 论文提交日期 2010年6月 论文答辩日期 2010年6月 学位授予单位和日期 成 都 理 工 大 学（ 年 月） 答辩委员会主席 评阅人 二○一○年六月 内 封 格 式 分类号 学校代码：10616 U D C 密级 学号：2007020443 成都理工大学硕士学位论文 中心回线瞬变电磁一维正反演算法研究 指导教师姓名及职称 王绪本教授 申请学位级别 硕士 专业名称 地球探测与信息技术 论文提交日期 2010年6月 论文答辩日期 2010年6月 学位授予单位和日期 成 都 理 工 大 学（ 年 月） 答辩委员会主席 评阅人 二○一○年六月 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得 成都理工大学 或其他教育机构的学位或证而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意。 学位论文作者签名： 年 月 日 学位论文版权使用授权书 本学位论文作者完全了解 成都理工大学 有关保留、使用学位论文的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁盘，允许论文被查阅和借阅。本人授权 成都理工大学 可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 （保密的学位论文在解密后适用本授权书） 学位论文作者签名： 学位论文作者导师签名： 年 月 日 TOC \o "1-3" \h \z \u 摘 要 I Abstract III 第1章 引 言 2 1.1瞬变电磁法的研究现状及发展趋势 2 1.1.1 方法简介 2 1.1.2 研究现状 3 1.1.3 瞬变电磁的发展趋势 6 1.2本文研究的主要内容 7 1.2.1层状均匀介质瞬变电磁响应正演理论 7 1.2.2瞬变响应曲线特征 8 1.2.3对几种全区视电阻率的定义方法的比较和分析 8 1.2.4反演成像 8 第2章 中心回线瞬变电磁一维正演理论 9 2.1电磁学基本理论 9 2.1.1 时间域和频率域麦克斯韦方程 9 2.2中心回线层状模型的正演理论 9 2.3 快速汉克尔变换的数值计算 15 2.3.1 Anderson数值滤波汉克尔变换 15 2.3.2 Guptasarma和Singh线性数值滤波 16 2.3.3 两种汉克尔变化的比较 16 2.4 余弦变换 18 2.4.1 折线逼近法求余弦变换 18 2.4.2 对折线逼近法求余弦变换的研究 19 2.4.3 数值滤波余弦变换 21 2.4.4 小结 21 2.5本章小结 22 第3章 瞬变电磁响应的正演模拟 23 3.1不同观测参数条件下的瞬变响应特征分析 23 3.1.1不同大地电阻率的瞬变响应特征的比较 23 3.1.2不同发射线框边长的瞬变响应特征的比较 23 3.1.3不同发射电流的瞬变响应特征的比较 25 3.2 不同地电参数模型的瞬变响应特征 25 3.2.1 D型地电模型 25 3.2.2 G型地电模型 26 3.2.3 A型地电模型 27 3.2.4 K型地电模型 27 3.2.5 H型地电模型 28 3.2.6 Q型地电模型 29 3.2.7 小结 29 3.3拟二维正演模拟 29 3.4本章小结 30 第4章 视电阻率的定义 31 4.1 早晚期视电阻率的定义 32 4.1.1 早期视电阻率 32 4.1.2 晚期视电阻率 32 4.2 全区视电阻率的定义 33 4.2.1 通过磁场强度法计算全区视电阻率 33 4.2.2 感应电动势法计算全区视电阻率 35 4.2.3 平移算法计算全区视电阻率 35 4.3 三种全区视电阻率的比较 37 4.3.1感应电动势法和平移算法全区视电阻率曲线重合原因的分析 37 4.3.2全区视电阻率曲线特征的分析 38 4.4 本章小结 41 第5章 成像与反演 43 5.1 烟圈反演 43 5.1.1 H型地电模型的烟圈反演 45 5.1.2 K型地电模型的烟圈反演 46 5.1.3 多层地电模型的烟圈反演 47 5.2 阻尼最小二乘反演 48 5.3 算例 50 5.3.1 H型地电模型的反演 50 5.3.2 K型地电模型的反演 51 5.1.1 多层地电模型的反演 51 5.4 小结 52 结 论 53 致 谢 55 参考文献 56 攻读学位期间取得学术成果 59 中心回线瞬变电磁一维正反演算法研究 作者简介：郭嵩巍，性别：男，1984年6月生，师从成都理工大学王绪本教授，2010年06月毕业于成都理工大学地球探测与信息技术专业，获得工学硕士学位。 摘 要 瞬变电磁法(Transient Electromagnetic Method，简称TEM)，有时也称时间域电磁法（TDEM）。它是利用不接地回线或接地线源向地下发送一次脉冲磁场，通过观测及研究二次涡流场随时间的变化规律来探测介质的电性。该方法是关断一次场后观测纯二次场，具有分辨率高、勘探的深度大等优点，是目前在能源、矿产、水文、工程、环境等领域广泛应用的物探方法。 本文首先回顾了瞬变电磁法的发展历程，简要的介绍了瞬变电磁法的特点、装置类型，并且从方法理论及仪器等方面介绍了国内外发展现状及发展趋势。 正演计算方面。采用将麦克斯韦方程变换到频率域中进行推导，得到中心回线装置的瞬变响应公式。鉴于瞬变电磁探测通常采用阶跃波，以及瞬变响应实部、虚部的奇偶性，化简瞬变响应公式，得到正演计算公式，其中包括两个关键步骤：汉克尔变换、余弦变换。关于汉克尔变换，本文对比分析了Anderson线性数字滤波法汉克尔变换和Guptasarma和Singh线性数值滤波汉克尔变换的计算精度和计算速度，证明在中心回线瞬变电磁一维正演中Guptasarma和Singh线性数值滤波汉克尔变换计算精度高、计算速度快，值得采用；余弦变换，本文对比分析了折线逼近法余弦变换和数字滤波法余弦变换，对折线逼近法选取不同的频点数和不同的频宽的情况下，计算精度、速度做了细致的分析；相比之下，相同的计算精度，数值滤波余弦变换计算速度明显优于折线逼近法余弦变换，建议采用数值滤波余弦变换。 利用上述正演数值算法，编写中心回线瞬变电磁一维正演程序，并进行了大量的模型模拟，具体包括D、G、 A、 H、 K、 Q 型地电模型，此外还利用一维正演做拟二维正演模拟。分析了正演的曲线形态，做了大量的对比试验，对瞬变电磁的处理解释有指导意义。 视电阻率方面。常规的中心回线瞬变电磁视电阻率定义一般为早期、晚期。就全区视电阻率的定义，本文给出了磁场强度法、感应电动势法以及平移算法三种全区视电阻率的定义方式。发现感应电动势法和平移算法曲线完全重合，并给出了曲线重合的原因。对比分析了几种全区视电阻率对低阻薄层、高阻薄层的反映情况、对低阻、高阻的灵敏程度。证明瞬变电磁响应对低阻层反映灵敏、高阻层反映一般。此外，还对感应电动势法全区视电阻率存在有假异常现象进行了分析，给出了假异常现象的特点，这将有助于更正确的利用全区视电阻率对瞬变电磁资料进行解释。 反演方面。烟圈反演是针对瞬变电磁视电阻率一种时深转换方法，本文烟圈反演基于全区视电阻率，通过对比分析了几种全区视电阻率的烟圈反演，证明烟圈反演对异常的反映效果要优于其全区视电阻率。最优化反演，本文选用了阻尼最小二乘法反演对几种典型地电模型进行反演。经过对比分析，建议综合运用以上两种方法，以烟圈反演作为快速反演，做定性解释；再以烟圈反演的结果作为阻尼最小二乘反演的初始模型，做精细反演，试验证明两种方法相结合不仅可以提高反演的计算效率，计算结果也更接近模型。 关键词：中心回线；快速汉克尔变换；余弦变换；全区视电阻率；烟圈反演 Study about 1D forward and inversion algorithm of centrl loop transient electromagnetic method Introduction of the author: Guo Songwei, male, was born in June, 1984 whose tutor was Professor Wang Xuben. He graduated from Chengdu University of Technology with Earth exploration and Information major and was granted the Master Degree in June, 2010. Abstract Transient Electromagnetic Method (TEM), also known as time domain electromagnetic (TDEM), is based on the process of transmitting the primary electromagnetic impulse to underground with ungrounded loop or ground source and analyzing changes of secondary field versus time to get the electrical characters of the medium. Because this method is only receiving the secondary field after turning down the first field power, it has the advantages of resolution, exploration depth and so on. TEM is now widely used in mineral exploration, water exploration, engineering and environmental geophysical exploration. Above all we review the development of transient electromagnetic method and briefly introduce the advantage and device type of TEM. We also introduce the current development and the future trend at the aspects of the theories and devices. On the aspect of forward calculation, we used the central loop transient response calculation, which is derived from the Maxwell's equations. The on-off step wave is usually used in transient electromagnetic prospecting. The real part and imaginary part of transient response is also used to simplify the formulas and acquire forward calculation formulas. There are two key steps --the fast Hankel transform and cosine transform. In this paper we analyze and compare the accuracy and computing speed of the fast Hankel transform between the linear digital filter of Anderson and the linear digital filter of Guptasarma and Singh. It's proved that the filter of Guptasarma and Singh, which we'll adopt finally, is better than that of Anderson on the aspects of accuracy and computing speed. We also analyze and compare the Cosine transform of the broken line approximation method and digital filter method. We analyze the accuracy and computing speed of the broken line approximation cosine transform briefly in different frequency points and bandwidth. On the comparison, numerical filtering cosine transform method has the higher speed than the broken line approximation method, so numerical filtering cosine transform method is advised to use. With the numerical forward calculation algorithm above, we make a computer program for central loop TEM forward calculations. Using this program, we simulate numbers of geo-electric models as the model types of D, G, A, H, K, Q. In addition we also use the one-dimensional forward to simulate the two-dimensional forward modeling which is called pseudo two-dimensional forward. With analyzing the shapes of these forward curves and comparing different experiments, it'll be a help for processing and interpretation of TEM data. On the aspect of apparent resistivity, the conventional central loop transient electromagnetic resistivity definition is classified as early and late apparent resistivity. In this paper, three definitions of all-time apparent resistivity are introduced as the magnetic field strength method, electromotive force method and translation algorithm method. We found that curves of the magnetic field strength method and electromotive force method are complete overlap, and then give the conclusion of the reason. Comparing the TEM showed it sensitivity between low resistivity layer and high resistivity layer and analyzed the sensitivity to low resistivity thin layer and high resistivity thin layer of different models. In addition, we also analyze the pseudo-anomies of electromotive force method and conclude the features of it, which will be benefit for explanation of TEM by using all-time apparent resistivity. On the aspect of inversion, the smoke loop inversion we used is based on smoke theory, which is a simple time-depth conversion method for TEM. In this paper, smoke loop inversion is based on all-time apparent resistivity. It is proved that the result of resistivity reflected from smoke loop inversion is better than all-time apparent resistivity. On the aspect of optimization inversion, we choose the LMF method and take inversions of several typical models with this method and conclude that it's better to use smoking loop inversion at first for qualitative explanation, to use the result of it as the initialized model of LMF inversion, and to use the result of LMF inversion for qualitative explanation. After experiments based the combination of the two methods it is proved that it not only improves the efficiency of the inversion calculation, but also makes the results more closer to real models. Keywords: Central loop; Fast Hankel transform; Cosine transform;All-time apparent resisitivity;Smoke loop inversion 第1章 引 言 1.1瞬变电磁法的研究现状及发展趋势 1.1.1 方法简介 瞬变电磁法[1][2](Transient Electromagnetic Method，简称TEM)，有时也称时间域电磁法（TDEM）。它是利用不接地回线或接地线源向地下发送一次脉冲磁场，地下的导电地质体在一次脉冲磁场的激发下，被感应出涡流，其大小取决于地质体的导电程度，并在空间形成二次瞬变磁场。一次脉冲磁场随着脉冲电流的关断而崩溃，但二次瞬变磁场并不立即消失，它将有一个过渡(衰减) 过程，随时间按指数规律衰减，衰减的快慢与地质体的电性参数、结构构造、体积形态等有关。一般导体越大、导电性越好，涡流的热损耗越小，二次瞬变磁场衰减就越慢。在一次脉冲磁场的间歇期间，利用线圈或接地电极观测二次涡流场，由于二次涡流场随时间变化，因而在其周围又产生新的磁场，称二次磁场。由于导电矿体内感应电流的热损耗，二次磁场大致按指数规律随时间衰减，形成瞬变磁场。二次磁场主要来源于良导电矿体内的感应电流，因此包含着与矿体有关的地质信息（纯异常）。二次磁场通过接收回线观测，并对所观测的数据进行分析和处理，据此，解释地下矿体及相关物理参数。 瞬变电磁探测[17]既包含剖面测量技术，又包含了测深技术，可以用于解决各种不同的地质问题。根据地质任务选择不同的技术，TEM 常用的测量技术按场源不同可分为电性源和磁性源两类，按接收、发射线圈的排列不同又可以分为同点偶极，偶极—偶极装置，大定源回线装置、中心回线装置，重叠回线装置等等（如图1-1 所示）。 图1-1 TEM剖面装置：(a)重叠或中心回线 (b)偶极偶极 (c)大回线源 TEM测深装置：(1)电偶源(2)线源(3)磁偶源(4)中心回线 相对于大地电磁测深（MT）和可控源音频大地电磁测深（CSAMT）等其他频率域电磁法相比，瞬变电磁方法的具有以下几个特点： （1）TEM分辨率高 由于观测的是纯二次场，可在近区（晚期）观测，甚至可采用同点组合(重叠回线、中心回线法)法观测，与探测目标藕合好，取得的异常幅度强、形态简单、分层能力强。 （2）TEM勘探的深度范围大 在目前的技术条件下，勘探范围浅可至几米、深可达几千米；随着仪器的进步，通过加大发射功率增大二次场的幅度，提高信躁比，勘探深度范围还能进一步的扩大。 （3）TEM地形影响较小 在高阻围岩地区地形起伏引起的假异常较小;在低阻围岩地区，由于是多道观测，早期道地形影响也较易识别。 （4）TEM信息丰富 一个脉冲激发可以观测整条衰减曲线，频率成分丰富；剖面测量与测深工作同时进行，提供的地电信息丰富，便于资料的解释。 （5）TEM资料信噪比高 人工源方法，且可很容易地实现多次叠加(可达几百次至几千次)，资料信噪比高；且通过选择不同的时间窗口进行观测，可以有效地压制噪声。 （6）等值现象较弱 对局部异常的反映突出，对低阻层具有较高的分辨能力，具有比其他频率域电磁方法（MT，CSAMT）更高的分辨率。 （7）TEM适应环境能力强 由于可采用不接地回线发射与接收，不存在接地电阻问题，在基岩出露区、冻土带、沙漠、水泥路面、江河湖海水面上均可进行测量。具有适用范围广、施工方便、工作效率高的特点。 （8）TEM应用领域广 瞬变电磁法可解决的地质问题有：能源、矿产勘查、水文、工程、环境地质调查、考古探测等。 1.1.2 研究现状 （1）正反演方面 前苏联的时间域电磁法始于20世纪50年代，发展了建场测深法，并建立了瞬变电磁法解释理论和野外施工方法技术，理论研究方面也一直走在世界前列，并成功地完成了瞬变电磁法的一维正、反演。该法使用电偶源激发，在距源r处用接收线圈测量垂直磁场，主要用于地震探测油气田效果不理想的地区[3]。而利用瞬变电磁法找导电矿体的概念最早是由加拿大地球物理学家J.R.Wait于1951年提出的，他在1951年发表了关于导电介质中瞬变电磁场传播的论文[9]。随着理论研究的深入和电子技术的进步，七十年代以来该方法又得到了进一步的发展。A.A.Kaufmann等人后来把瞬变电磁测深法引入美国，推动了一维瞬变电磁理论进行的系统研究[5]。1971年J. R. Wait给出了层状均匀介质上的测深理论[9]。关于层状半空间瞬变场定量解释方法也由H. F. Morrison等人提出[18]，但是直到1981年Raiche和Spies给出了适合于延时、电导率和层深改变的二层均匀大地的理论曲线[19]。 由于实际TEM的观测数据均为三维响应，研究二维、三维模拟很有必要。随着计算机技术的发展，自二十世纪80年代以来，随着计算技术的发展，欧美各国在瞬变电磁法的二，三维正演摸拟技术方面(有限元，有限差分，积分方程及混合方法直接解时间域热传导方程或者先解频域亥姆霍兹方程，再进行域的转换)亦做了大量的计算和研究工作。1980年，Kuo J. T.和Cho D. H.首次用有限元法解时域中的变分方程，求任意二维地电断面的瞬变响应并用中心差分来代替热传导方程中对时间的导数[10]。1986年，Goldman等人用有限单元法求解任意二维电阻率分布的TEM响应，求解的是总场而非二次场。1992年，Leppin解决了长久以来未曾解决的二维非均质体、三维源的时间域电磁模拟问题。他间接求解时间域散射磁场的扩散方程，通过傅氏变换转到空间一波数域中，化为二维问题:用显式格式方法计算了微分方程后，用逐次超松驰选代求解大型线性方程组，之后再通过数字滤波返回到时间域[10]。当采用60×60个网格剖分、300次时间步长迭代时，求解6612个线性方程(取9-15个波数重复计算)，在COMPAQ386机上(带NDP FORTRAN编辑器)计算需60小时，可见计算量相当大[10]。1985年，San Filpo和Homann用积分方程法求解了导电半空间三维体的TEM响应。1993年，Wang和Honmann解决了“真正意义”上的三维TEM正演问题[10]。采用了三维空间网格的有限差分法直接在时间域中求解任意源的三维瞬变场响应。用交错网格法(staggeredgrid)求解电磁方程的时间分步，避开求解二阶Maxwell方程而直接求一阶的方程，无需求各参数的空间导数，故而解决了许多有关数值计算上的难点，计算的磁场晚期结果准确。当区域剖分100×100×50个网格节点，选取最大电阻率达100Ωm，最小网格间距为lΩm，计算l0m的场，在IBM3090 / 600S机上约需3. 5小时[10]。 此外，基于电磁场本身的叠加原理，从麦克斯韦方程组出发，导出了中心回线瞬变电磁2.5维二次场(纯异常)的有限单元计算公式。该算法采用三角形有限元网格，在尽可能拟合地下电性断面的情况下减少有限元网格的节点数和单元数，用选主元的LU分解法求解线性方程组。使用1057个节点，2070个单元，用P3-733,128 M内存，175 min(其中一次场计算占50以上)完成一条9个测点、15个时间点的2 .5维剖面的正演计算。 我国的瞬变电磁法研究起于70年代初，由长春地院、地矿部物化探所、中南工业大学、中国地质大学等单位分别在方法理论、仪器及野外试验方面做了大量工作，在一维正、反演及方法技术理论上取得了大量的成果，长期从事这方法研究工作的有：方文藻、朴化荣、蒋邦远、罗延钟、牛之链、王庆乙等。近十余年来，在多维正演模拟中也取得了很大的进展，如:殷长春、刘斌(1994)利用并矢格林函数理论和积分方程方法计算两层大地中三维异常体的频率域电磁响应，并利用余弦变换方法将其转换为时间域电磁响应[12]。从“七五”和“八五”至今，我国在一维瞬变电磁正反演理论方面所做的工作比较突出。朴化荣系统总结和研究了电磁测深理论，并实现了通过快速汉克尔变换在频率域进行层状大地的正演计算，再利用G-S逆拉氏变换法[32][34]或正（余）弦变换法转到时间域响应的算法[35]，其计算的精度和速度都与国外当时的水平相近 从瞬变电磁法的发展与研究现状来看，多维时域电磁法的正演计算速度是一个关键问题，特别是2.5维、完全三维问题的正演模拟，计算量很大，关于多维时域电磁法反演问题的研究进展缓慢。 （2）视电阻率定义方面 在视电阻率定义方面，由感应电动势按通常方法定义的全区视电阻率,难以解析的给出全时间段视电阻率，广泛采用的多是早、晚期近似的，不能连续直观地给出电性层变化的信息，对反演及定性解释带来了困难。 于是全区视电阻率的定义和计算成为了研究人员急需解决的难题，为此利用长谷川健[45]的思想实现了全区视电阻率的定义，许多学者做了大量的研究。白登海（2003年）[37]给出了一种瞬变电磁法视电阻率的数值计算方法,该方法根据中心方式磁场垂直分量时间变化率的核函数的表现特征,把整个瞬变过程分为早期阶段、早期到晚期的转折点和晚期阶段。分别得到早期视电阻率和晚期视电阻率,然后通过转折点构成一条完整的全程视电阻率曲线。杨生（2003年）[39]给出了中心回线装置发射电流为斜阶跃波形条件下全区视电阻率迭代反演计算方法。李建平（2007年）[40]把回线分解为水平电偶极子,然后给出利用电偶极子求取全区视电阻率的方法。王华军（2008年）[42]依据均匀半空间瞬变响应曲线随地下电导率、发射回线边长与观测时间具有平移伸缩特性,提出了一种直接计算全区视电阻率的方法。付志红等人（2008年）[41]研究了斜阶跃场源瞬变电磁法的全区视电阻率数值计算。虽然研究的人员很多,但直到现在,还没有一种比较完善的计算全区视电阻率的方法。本文通过正演模拟，重点对比分析了其中颇具代表性的三种数值计算方法，总结了它们的优缺点。应该说，对全区视电阻率的定义和计算方法的改进仍有很大的空间。 （3）反演成像方面 一维层状瞬变响应资料的反演解释主要沿着两个方向发展：一是基于严格的正演公式，用最优化方法进行反演。即给定某个初始模型，采用多次迭代的方法不断修正模型，使其正演结果与观测数据达到最小二乘意义下的最佳拟合，从而得到反演结果。1985年Raiche还提出用共框TEM测深数据和对称四极电阻率测深数据进行联合反演[9]。由于是基于严格的正演公式，层数较多时正演计算时间较长，反演计算速度较慢；此外与直流电测法一样，也存在等值性的问题，反演存在多解性。另一种用作解释与反演的理论是“烟圈”理论，M. N. Nabighian(1979)[3][47]提出：瞬变响应可近似地用向地下扩散的电流环来等效，这些电流环好象是发射回线吹出的“烟圈”，其形状与发射回线相同，随着延时的增加而向外、向下扩一散。根据这个原理，可把地表测得的随时间变化的瞬变响应转化为电阻率随深度变化的函数曲线，从而实现瞬变响应的反演，这是一种快速、有效的近似反演方法。 PerryA. Eaton及Gerald W.Hohmann发展了一种类似的快速反演方法用于解释三维瞬变响应数据，效果也较好[53]。但近似反演的结果往往不能拟合观测数据，精度还需要进一步提高。S. Lee根据电磁波与地震波的相似性，利用有限差分法实现了电磁波场的二维偏移和电磁数据的成像。 （4）仪器方面 近三十年来，欧美各国许多大学和科研机构也竞相发展了瞬变电磁系统，1958年加拿大Barringer公司研制的航空INPUT系统在世界各地进行了矿产普查，而作为真正的地面用的瞬变电磁勘探用的商业性仪器直到上个世纪七十年代才出现。电法仪器过去多为单一功能的仪器，如澳大利亚的SITOTEM系统和EM-37系统，主要用于浅层勘探。加拿大的UTEM系统是功能比较先进的仪器，但是由于使用感应源，只能探测良道题，同时探测深度也不大。到了九十年代，由于计算机技术的应用，使得在数据采集、资料处理与解释方面达到了较高的水平，世界上已开发多功能、集成化的电法仪器，如加拿大凤凰公司的V8系统。现代仪器的发展趋势是向大功率（100千瓦），GPS同步，多道叠加、24位模数转换的高精度多道测量系统发展，这些先进技术为瞬变电磁法的发展和应用提供了广阔的天地[13][14]。 总之，我国的瞬变电磁法的理论、方法及仪器研制等方面都与国外有相当大的差距。瞬变电磁一维正、反演计算理论已趋于成熟，虽然学者做了大量的二、三维的正、反演模拟，但仍然无法应用于生产，在生产中仍然以一维为主。此外，到目前为止，瞬变电磁视电阻率的定义依然没有很好的解决，对二三维正演模拟的结果的视电阻率的定义依然沿用一维的定义。 1.1.3 瞬变电磁的发展趋势 由于瞬变电磁法是寻找铜多金属矿的有效方法之一,尤其是在深部矿产勘查方面,地-井TEM法已成为常规的勘探方法,并已取得显著的效果。对瞬变电磁法来说,随着矿业的发展,将被更加广泛的应用,要求瞬变电磁法具备较大勘探深度和较高的纵向分辨能力。对瞬变电磁仪器而言,大功率、大动态范围、高密度时序序列数据采样、三分量同步观测、低噪声仪器性能将是国外先进TEM仪器发展的主要趋势。国内TEM仪器目前应下大力气解决稳定性问题,完善仪器的功能,如石英钟同步、发射电流下降沿线性且无振荡、发射电流下降沿时间自动调整或可自测、配套探头、降低仪器内噪声等技术功能。 在解释技术方面,急需解决的重要内容之一是数据处理技术,其中关键为消除干扰噪声和装置过渡过程的方法技术,这二方面是影响解释成果质量的重要因素。三维反演解释技术将是长期的研究目标,在近期达到实用的程度可能较小。一维反演和二维电性分布成像仍将是近期主要的解释手段。三分量电阻率成像是研究的主要目标之一,目前尽管有多种仪器已实现了三分量数据采集,但对于水平分量信息的应用仍处于曲线形态特征比较的水平,水平分量能较准确的圈定导电地质体的边界和产状。电磁偏移成像和电磁波场转换解释技术需要发展和完善,它可实现复杂地电条件下的快速确定导电地质体的空间位置和基本形态。地-井TEM资料的解释方法目前主要是依据曲线形态推断目标体的空间深度和形态,以及用三分量交汇法显示导体距钻孔的距离,因此,三分量地-井TEM资料的电阻率成像和三维反演也将是研究的主要方向。 高温超导磁强计用于TEM法的传感器,已显示出明显的优势,可增大勘探深度1. 5倍左右,但目前由于其抗干扰能力较差,应用范围受环境因素限制,解决好高温超导磁强计的抗干扰问题,使之能在多数干扰环境下应用,将大大推动TEM法在勘查应用领域的进步[13][14][15]。 1.2本文研究的主要内容 1.2.1层状均匀介质瞬变电磁响应正演理论 正演不仅可以对特定的地电模型进行模拟，同时也是资料处理解释的基础。本文介绍了中心回线瞬变电磁法阶跃响应的正演数值计算方法：首先推导得出中心回线瞬变电磁法阶跃响应在频率域中的响应表达式，通过汉克尔变换得到频率域瞬变响应，再做余弦变换得到时间域瞬变电磁响应。 汉克尔变换，本文就Anderson线性数字滤波汉克尔变换和Guptasarma和Singh线性数值滤波汉克尔变换两种汉克尔变换，对它们在中心回线瞬变电磁一维正演中的计算速度、计算精度做了对比分析。 余弦变换，本文介绍了折线逼近法余弦变换和数字滤波法余弦变换，对折线逼近法在选取不同的频点数和不同的频宽的情况下，计算精度、速度做了细致的分析。 通过对几种汉克尔变换、余弦变换的对比分析，最终给出了一种在保证一定精度情况下，计算效率最高的正演数值计算方案。 1.2.2瞬变响应曲线特征分析 不同观测参数条件下瞬变电磁响应曲线特征分析中心回线瞬变电磁响应不仅与地层参数有关，还与观测参数有关，因此有必要分析其在不同观测参数条件下的响应特征。本文分别对不同大地电阻率、发射线框、发射电流、进行了模拟；此外还对几种常见的地电模型进行了模拟。具体包括二层 D、G，三层 A、K、H、Q；四层模型包括KH，五层模型包括HKH。 此外，本文还对直立板状体进行了拟二维模拟。 1.2.3对几种全区视电阻率的定义方法的比较和分析 本文对三种全区视电阻率的定义方式进行了对比分析，包括利用磁场强度法，感应电动势法以及平移算法，分析它们的曲线形态特征。 1.2.4反演成像 本文在成像方面首先采用了烟圈反演，它是TEM最常用、计算速度最快的一种快算成像方法。本文的烟圈反演是基于全区视电阻率做的成像，并对比分析了不同全区视电阻率的烟圈反演的成像效果。 本文在反演方面还采用基于一维正演的阻尼最小二乘法反演，对几种地电模型进行了反演，通过对正演模拟数据的反演结果分析了它们的特点及效果。 结合两种反演方法的特点，采用烟圈反演结果作为阻尼最小二乘法的初始模型，对几种模型进行反演。 第2章 中心回线瞬变电磁一维正演理论 2.1电磁学基本理论 2.1.1 时间域和频率域麦克斯韦方程 一个电磁场可用定义域中的4个矢量函数E--电场强度、B--磁感应强度、D--电位移、H--磁场强度来确定，所有的电磁现象都服从麦克斯韦方程[7]，时间域中麦克斯韦方程[7]有如下形式： （2-1） （2-2） （2-3） （2-4） 式中j为电流密度，ρ为电荷密度。方程(2-1)至(2-4)是一组关于5个矢量 、 、 、 和 的互相独立的微分方程，可通过描述场矢量之间关系的状态方程(2-5),(2-6)和(2-7)把它们联系起来。 对时间域的麦克斯韦方程(2-1)至(2-4)做一维傅氏变换，并利用状态方程 (2-5),(2-6)，(2-7)。 （2-5） （2-6） （2-7） 其中 为介电常数， 为磁化率， 为电导率，得到频率域麦克斯韦方程组： （2-8） （2-9） （2-10） （2-11） 其中 为频率，以弧度/秒为单位。 2.2中心回线层状模型的正演理论 图2-1 中心回线地电模型示意图 对于中心回线装置而言[18]，如图2-1所示，发射线框半径为a，置于距地面高h., 假定圆柱坐标系下，在距离发射回线r产生的感应电流密度为： （2-12） 其中： 方向为电流线的任意方向； 代表发射线框电流的频谱表达式； r接收线框到发射线框的距离（偏移距）； 沿 方向积分，得到场源引起电流密度为： （2-13） 对于磁偶极子场可由 、 、 来表示。线圈水平放置，所以 ，于是通过傅立叶变换将麦克斯韦方程变换到频率域仅有一下三个方程： （2-14） （2-15） （2-16） 可以看到，消去 ，解上述偏微分方程组（2-14）、（2-15）、（2-16）式，便可以得到 及 ；同时，将（2-14）、（2-15）式带入（2-16）式，消去 、 ，得到： （2-17） 其中，一般来说 （2-18） 定义汉克尔变换对： （2-19） （2-20） 对（2-17）式变换,得： （2-21） 其中： （2-22） 方便起见: （2-23） 对于任意层来说，都是均匀介质，于是电场可表示成： （2-24） 对于一个n层地电模型，引入波阻抗的概念，第i层为： （2-25） 根据（2-14）式，波阻抗还可表示才成： （2-26） 我们便得到了从第n层介质到第一层介质的波阻抗递推表达式： （2-27） 其中：hi为第i层的厚度，数值计算从最底层（n层）开始: （2-28） 逐层递推，最终算到第一层。 因此，任意层介质的波阻抗我们都可以通过λ和ω求出。利用（2-23）、（2-24）、（2-25）式，我们可以得到如下关系式： （2-29） （2-30） 联立（2-29），（2-30）式，解出： （2-31） 由（2-23）式，得到： （2-32） 替换（2-23）（2-21）式，并写成汉克尔反变换形式，得到半空间响应的方程： (2-33) (2-34) 利用（2-15）式，得到垂直方向的磁场响应： (2-35) 式中： （2-36） 由于上半空间为空气层，导电率σ=0，因此： 波阻抗 最后，通过傅里叶反变换将频率域瞬变响应变换回时间域： （2-37） 事实上，（2-33）、（2-34）、（2-35）式中含有零阶和一阶贝塞尔函数，依然很难进行数值计算，因此需要对其进行进一步化简： 对于低频而言，只考虑位移电流： 对于中心回线而言（偏移距为零）： 对于地面接收而言： 于是： （2-38） 将（2-38）式代入（2-35）式，（2-35）式化简为： （2-39） 因此地面中心回线在频率域的响应为：， （2-40） 其中： a为圆形发射回线的半径(米), 若发射回线是边长为L(米)的方形回线时,可通过两者面积相等来计算其等效半径[2]: （2-41） 从目前的仪器来看，瞬变电磁多采用占空比为1：1阶跃方波（如图2-2），又称 TD50波，本质上还是阶跃波[21][22]（如图2-2）。使用TD50波目的是为实现多次叠加，正反向供电为降低电磁耦合，提高信噪比。 阶跃波： （2-42） 对其进行傅里叶变换： （2-43） 图2-2阶跃波示意图 图2-3 TD50波形示意图 将（2-43）式带入（2-40）式，再带入（2-37）式做傅里叶反变换，即可得到实现数值计算。不过，针对表达式的特点，利用 虚部、实部的对称性（实部偶函数、虚部奇函数），我们可将傅里叶变换简化，便得到可以利用正余弦变换的积分表达式[20]（本文在数值模拟中选择了余弦变换）： （1）利用场强虚部做计算 （2-44） （2-45） （2）利用场强实部做计算 （2-46） （2-47） 最后，根据法拉第电磁感应定律[7]，得到二次场的感应电动势为： （2-48） 式中：S为接收线框面积，N为接收线框扎数。 于是我们就得到了中心回线瞬变电磁法的数值计算表达式（2-44、2-45），可以看出：数值计算的关键就是解决汉克尔变换及余弦变换的数值计算问题，这将在本章后面的小结中详细的介绍，正演数值计算的流程图如下： 图2-4 中心回线瞬变电磁正演数值计算流程图 2.3 快速汉克尔变换的数值计算 前文推导中已经提到了汉克尔变换对，且最终的表达式中含有一阶贝塞尔函数。贝塞尔函数不仅是一个高频振荡函数，而且是一个慢衰减函数。这种积分除了在均匀半空间条件下可以求出解析解外，其它情况下一般求不出解析解，故只能采用数值积分方法或数字滤波方法[33]。一般来说，数值积分方法的计算时间相对于数字滤波方法一般都比较长，因此本文采用数字滤波方法计算。汉克尔变换数值滤波算法众多，本文选择了其中两个算法，其原理如下： 2.3.1 Anderson数值滤波汉克尔变换 Anderson[24][25][27][28]采用了线性数字滤波法，取得了良好的效果，其原理如下： 考虑如下一阶贝塞尔函数的汉克尔变换: （2-49） 令 、 ,即 ， ,带入(2-49)式，并等式两边同乘ex得： （2-50） 即： （2-51） 该式为一个卷积式k(e-x)为输入函数，exf(ex)为输出函数，[ex J1(ex)]为滤波项，于是求出滤波系数，包含贝塞尔积分的汉克尔变换可采用数字滤波法计算。 由傅里叶变换理论，上述卷积式等价于变换域中的乘积，即： （2-52） F(s)、K (S)、H(S)分别为 exf(ex)、k(e-x)、[ex J1(ex)]的傅氏变换式，则： （2-53） 对H(S)做反傅里叶变换即可得到汉克尔变换滤波系数，Walter L.Anderson选取的汉克尔变化对为： （2-54） 求得汉克尔变换滤波系数之后，（2-49）式得积分变为一个离散求和式近似 的表示为： （2-55） 式中：Hi为汉克尔滤波系数，Ai-x为移动的横坐标。 2.3.2 Guptasarma和Singh线性数值滤波 Guptasarma和Singh线性数值滤波[31][33]计算公式: （2-56） （2-57） 数值计算的精度取决于积分区间的长度n、抽样点的位置 和加权系数 ，通常n越大计算精度越高。Guptasarma和Singh (1997年) 给出的61点汉克尔 变换线性滤波器和47点汉克尔 变换线性滤波器系数。此外，Guptasarma和Singh还给出了120点汉克尔 变换线性滤波器和140点汉克尔 变换线性滤波器系数。本文的数值计算存在 汉克尔变换，采用47点线性滤波系数。 2.3.3 两种汉克尔变化的比较 在第四章第一节中我们将具体介绍早晚期视电阻率的定义，这里我利用早期视电阻率的早期结果可靠、晚期视电阻率的晚期结果可靠的特点，将早晚期视电阻率曲线作为判断正演计算的有效范围的评价依据。 采用均匀半空间模型，大地电阻率取1.0欧姆米；发射线框200米*200米，匝数为1，发射电流强度为10安培；接收线框面积为1平方米，匝数为1，接收时间范围从1.0e-6秒到1.0e+2秒，指数等间隔0.1，共81个时刻点。计算结果如下表： 在余弦变换相同的情况下： 表2-1 两种汉克尔变换计算结果 汉克尔滤波系数 计算时间（秒） 时间序列有效范围（秒） Guptasarma（47点） 0.578 1.0e-6~1.0e+3 Anderson（283点） 3.11 1.0e-6~1.0e+1 说明：计算机双核cpu，主频2.0GHz，后文中模拟均使用该计算机。 图2-5 早期视电阻率对比图 图2-6 晚期视电阻对比图 从以上的模拟结果可以得出： 一般来说，加长滤波系数能明显提高计算精度,使正演结果更接近真实模型;但加长滤波系数却意味着增加计算量，使计算效率降低。 在计算精度方面，针对不同的函数应该说计算精度不尽相同，因此本文就中心回线瞬变电磁一维正演的模拟利用以上两种汉克尔变换方法进行计算，两种算法都得到了较高的精度（如表2-1），感应电动势曲线基本重合，但Guptasarma的计算速度明显快的多（0.578秒）；从早、晚期视电阻率来看，早期的准确范围相差不大，晚期Guptasarma汉克尔变换明显好于Anderson汉克尔变换（如图2-6、表2-1）。 因此，Guptasarma和Singh线性数值滤波算法较Anderson的汉克尔变换不仅计算速度要快，而且计算精度高。建议采用Guptasarma和Singh汉克尔变换滤波系数（下文的数值模拟中的汉克尔变换全部采用Guptasarma和Singh线性数值滤波法）。 2.4 余弦变换 在地球物理数值模拟的过程中，经常遇到频率域到时间域转换，一般来说采用傅里叶变换。这时候往往会根据函数本身的特殊性改用其它变换，已达到提高计算精度、计算速度的目的。在TEM一维正演中，可以根据响应函数实部偶函数、虚部奇函数的特点，将傅里叶变换化简为正余弦变换[20][35]： （2-58） （2-59） 当ω非常大时，该计算变得非常困难，主要原因是由于正余弦函数的震荡的很强烈，通常的近似积分公式很难保证其数值计算的精度。为此人们提出了各种数值计算方法，主要包括数字滤波法[35]、折线逼近法[8]等。本文在数值模拟中采用了余弦变换，因此对就不对正弦变换具体阐述了。 2.4.1 折线逼近法求余弦变换 针对（2-58）式的余弦变换式 根据傅氏变换的性质： （2-60） 如果 K (ω)是足够光滑的函数，则可将上式积分区间分段，并在每个段内用折线来逼近 K(ω)，这时二阶导数K ''(ω)变成一系列δ函数： （2-61） 带入上式进行积分，并重新组合成： （2-62） 上式中取了有限项，N 的大小取决于函数 K (ω)的分布形态及计算精度。当： ，及 可以用于计算磁场强度及磁场强度变化率（关于时间的导数），利用（2-44）式、（2-45）式： 得到： （2-63） （2-64） 以上两式适用于所有频率域电磁法正演结果到时间域瞬变电磁法的转换。此外，余弦变换还可以用滤波方法来求解。 2.4.2 对折线逼近法求余弦变换的研究 与汉克尔变换类似，我们依然选用将早晚期视电阻率曲线作为判断正演计算的有效范围的评价依据。 （1）对不同频点数的研究 模型采用均匀半空间模型，大地电阻率取10欧姆米；发射线框200米*200米，匝数为1，发射电流强度为10安培；接收线框面积为1平方米，匝数为1，接收时间范围从1.0e-10秒到1.0e+2秒，指数等间隔0.2，共61个时刻点。 正演中采用Guptasarma和Singh汉克尔变换，滤波系数47点；余弦变换频率范围1.0e-5到1.0e+7，频点数分别取120，1200,12000，计算结果如下表： 表2-2 不同频点余弦变换计算结果 选择频点数 有效时间范围（指数） 计算时间（秒） 12000 -9 ~1 1.797 1200 -7 ~-2 0.1870 120 -7 ~-3 0.016