（地球探测与信息技术专业论文）瞬变电磁勘察煤田采空区的可行性分析.pdf

摘要 采空区是重要的煤田自然灾害，它源于煤田开采并将威胁到煤田 的后续开采及煤田范围的地质安全。瞬变电磁方法是一种建立在电磁 感应原理基础上的，通过分析提取并研究人工源激励下地质体的响应 信息，达到探测地下地质体目的的时间域人工源电磁探测法。应用瞬 变电磁勘察煤田采空区已有几十年的历史，在仪器不断更新的今天， 讨论瞬变电磁勘察煤田采空区的可行性很有必要。 本文从瞬变电磁2 5 维正演原理出发，首先用f o r t r a n 语言编写 程序，实现瞬变电磁2 5 维有限元正演，通过计算已有解析解的一维 模型证明算法的正确性，然后对不同充水情况的采空区进行模型抽象 并用2 5 维有限元正演算法进行数值模拟，得到各种情况下的瞬变电 磁响应并换算成实际工作中的测量参数a v ，进而讨论在现有仪器精 度及地质噪声情况下，电磁响应a v 的可分辨程度，实现对煤田采空 区不同赋水情况下的瞬变电磁响应的讨论，分析了瞬变电磁勘察煤田 采空区的可行性。 本文的创新点在于，改变过去用磁矩来讨论瞬变电磁可行性分析 的模式，直接用实测参数来做可行性分析，对采空区的讨论从一维上 升n - - 维。 关键字瞬变电磁，有限单元法，2 5 维正演，煤田采空区 a bs t r a c t g o a fi sa ni m p o r t a n tn a t u r a ld i s a s t e ri nc o a l f i l e d i ti sat h r e a t c o m i n gf r o m t h ec o a l f i e l d m i n i n ga n dt h r e a t i n gf o l l o w u pc o a l f i e l d m i n i n ga n dt h eg e o l o g ys a f e t yo fc o a l f i e l d t r a n s i e n te l e c t r o m a g n e t i ci sa t i m e d o m a i nd e t e c t i o nm e t h o dw h i c hi sb u i l to nt h eb a s i so ft h ep r i n c i p l e o fe l e c t r o m a g n e t i ci n d u c t i o n i td e t e c t st h eu n d e r g r o u n dg e o l o g i c a lb o d y b ya n a l y z i n g ，e x t r a c t i n ga n ds t u d y i n gi n s p i r e dg e o l o g i c a li n f o r m a t i o ni n r e s p o n s et oa r t i f i c i a le l e c t r o m a g n e t i cs o u r c e a p p l i c a t i o no ft r a n s i e n t e l e c t r o m a g n e t i ci ns u r v e y i n gt h eg o a fo fc o a l f i e l dh a sd e c a d e so fh i s t o r y w i mt h e u p d a t i n go ft h et r a n s i e n te l e c t r o m a g n e t i ca p p a r a t u s ，i ti s n e c e s s a r yf o ru st od i s c u s sf e a s i b i l i t yo ft r a n s i e n te l e c t r o m a g n e t i ci n s u r v e y i n gt h eg o a fo fc o a l f i e l d i nt h i s a r t i c l e ，b a s i n g o nt h e p r i n c i p l e o f2 5 dt r a n s i e n t e l e c t r o m a g n e t i cf o r w a r d ，f i r s to fa l l ，t h ew r i t e ra c h i e v e s2 5 df i n i t e e l e m e n tf o r w a r db yf o r t ra np r o g r a m m i n g s e c o n d l y , t h r o u g ht h e c a l c u l a t i o no fo n e d i m e n s i o n a lm o d e lw h i c hh a sa n a l y t i c a ls o l u t i o n t h e c o r r e c t n e s so f2 5 一df i n i t ee l e m e n tn u m e r i c a ls i m u l a t i o na l g o r i t h m si s p r o v e d a f t e rt h a t ，t h ew r i t e rg i v e sg o a fa b s t r a c tm o d e l si nd i f f e r e n t w a t e r r i n gs i t u a t i o na n dc a l c u l a t e st h e i rt r a n s i e n te l e c t r o m a g n e t i cr e s p o n s e u s i n g2 5 - df i n i t ee l e m e n tn u m e r i c a ls i m u l a t i o na l g o r i t h m sa n dc o n v e r t s t h e mi n t ot h ea c t u a lm e a s u r e dp a r a m e t e ra g ，t h e nd i s c u s s e sw h e t h e rt h e g o a f se l e c t r o m a g n e t i cr e s p o n s eo fd i f f e r e n tw a t e r i n gs i t u a t i o nc a nh e d i s t i n g u i s h e di nt h ec o n d i t i o no ft h ea c c u r a c yo ft h ee x i s t i n ge q u i p m e n t a n dg e o l o g i c a ln o i s e a tt h es a m et i m e ，t h ed i s c u s s i o no ff e a s i b i l i t yo f t r a n s i e n te l e c t r o m a g n e t i cs u r v e y i n gc o a l f i e l dg o a fi sr e a l i z e d t h ei n n o v a t i o no ft h i sa r t i c l ei sc h a n g i n gt h ef e a s i b i l i t ya n a l y s i s m o d eo ft r a n s i e n te l e c t r o m a g n e t i cf r o mm a g n e t i cm o m e n tt ot h ep r a c t i c a l p a r a m e t e r sa n di m p r o v i n gt h ed i s c u s s i o no fg o a ff r o mo n e d i m e n s i o n a lt o t w o d i m e n s i o n a l k e yw o r dt r a n s i e n t e l e c t r o m a g n e t i c ，f i n i t e e l e m e n t m e t h o d ， 2 5 - d i m e n s i o n a lf o r w a r d ，c o a l f i e l dg o a f 一i i 原创性声明 本人声明，所呈交的学位论文是本人在导师指导下进行的研究 工作及取得的研究成果。尽我所知，除了论文中特, 另j d n 以标注和致谢 的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不 包含为获得中南大学或其他单位的学位或证书而使用过的材料。与我 共同工作的同志对本研究所作的贡献均已在论文中作了明确的说明。 作者签名：粱之壁】 日期：丝竖年月堑日 学位论文版权使用授权书 本人了解中南大学有关保留、使用学位论文的规定，即：学校 有权保留学位论文并根据国家或湖南省有关部门规定送交学位论文， 允许学位论文被查阅和借阅；学校可以公布学位论文的全部或部分内 容，可以采用复印、缩印或其它手段保存学位论文。同时授权中国科 学技术信息研究所将本学位论文收录到中国学位论文全文数据库， 并通过网络向社会公众提供信息服务。 作者签名：雪幽导师签 中南大学硕士学位论文 第一章绪论 1 1 引言 第一章绪论 瞬变电磁方法是一种建立在电磁感应原理基础上的时问域人工源电磁探测 法。它利用不接地回线( 磁源) 或接地线源( 电偶源) 向地下发送一次脉冲磁场 ( 通常为一次场) ，在其激发下，地下地质体中激励起的感应涡流将产生随时问 变化的感应电磁场( 通常称二次场) 。由于二次场包含有地下地质体的丰富信息， 在一次脉冲磁场间歇期间，利用线圈或接地电极观测二次场( 或称响应场) ，通 过对这些响应信息的提取与分析，达到探测地下地质体的目的。 煤田采空区是在煤田开采过程中形成的。煤层被挖空使得上覆岩层的应力平 衡被破坏，加之不同程度的充水，形成采空区相对围岩表现出高阻或低阻异常口1 ， 这是瞬变电磁勘察煤田采空区的地球物理前提。 1 2 瞬变电磁正演问题的研究现状 时间域电磁法应用于探测地质构造问题起始于前苏联。我国于2 0 世纪7 0 年 代初期开始研究瞬变电磁，先后投入研究的单位有长春地质学院( 现吉林大学) 、 中国地质科学院物化探勘查研究所、桂林有色矿产地质研究院、中南工业大学( 现 中南大学) 、西安地质学院( 现长安大学) 、西安交通大学、北京有色矿产地质 研究院、中国地质大学等。国内学者主编的主要书籍有朴化荣电磁法测深原理、 牛之琏时间域电磁法原理、方文藻瞬变电磁测深原理、蒋邦远实用近区 磁源瞬变电磁法勘探、李貅瞬变电磁测深的理论与应用等。 瞬变电磁数值模拟方法主要分为三种：有限差分法、积分方程法和有限单元 法，这三种方法各有优缺点。 ( 1 ) 、有限差分法( f i n i td if f e r e n c em e t h o d 或简称为f d m ) 是一种经典的数 值模拟计算方法，其基本原理就是用差商代替微商，将待求解的连续微分方程变 换为离散的差分方程，并通过求解差分方程得到原微分方程的近似解。 ( 2 ) 、积分方程法( i n t e g r a le q u a t i o n ) 是从场参数( 场强或场位) 所满足的微 分方程边值问题出发，通过某些变换导出有关参数( 场强、位、积累电荷密度或 散射电流密度) 所满足的积分方程，然后用近似计算方法求此积分方程的数值解。 中南大学硕士学位论文第一章绪论 ( 3 ) 、有限单元法( f i n i t ee l e m e n tm e t h o d 或简称f e m ) 是将要分析的连续场 分割为很多较小的区域，它们的集合代表原来的场，然后建立每个单元上待求场 量的近似式，再结合起来进而求得连续场的解。从数学角度上来讲，它是从变分 原理出发，通过区域剖分和分片插值，把二次泛函的极值问题化为多元二次函数 的极值问题，后者等价为求解一组多元线性代数方程组，是一种从部分到整体的 方法，可使分析过程大为简化。 总的来说，三种数值模拟方法在瞬变电磁正演计算方面都有一定的优势和不 足，有限差分法的优点是方法简便易算，其缺点是当物性参数复杂分布或场域的 几何特征不规则时，适应性比较差。而积分方程法的优势是正演速度快、内存需 求少，主要用于地形改正和地下少量地质体的正演模拟。有限单元法与前述方法 相比，在瞬变电磁正演方面有独到的优势： ( 1 ) 在变分问题中，自然边界条件已经隐含地得到满足，只需考虑加强边 界条件( 第一类边界条件) ，推导过程简单； ( 2 ) 在处理复杂的几何形状时，其灵活性和适应性比其他方法要好； ( 3 ) 适用于多种介质和非均匀连续介质问题。由于多种介质和非均匀介质 是地球物理场域的基本特征，这是其他数值模拟难以胜任之处； ( 4 ) 有限单元法方程的系数矩阵是正定的，保证了解的存在唯一性； ( 5 ) 对于二阶偏微分方程，其变分问题只含有一阶导数，大大降低了偏微 分的处理难度。当然，有限单元法缺点是计算量比较大，计算时间长，但目前计 算机的内存容量和c p u 的迅速发展，使有限单元法在解决许多工程领域的数学物 理问题中，成为一种高效、通用的计算方法。将有限单元法引入地球物理计算方 法中，已经解决了许多以前无法计算的地球物理问题。 2 5 维人工源时间域电磁场数值模拟( 正演) 问题是目前国际上尚未妥善解决 的计算地球物理疑难问题之一。对它的研究始于2 0 世纪7 0 年代中期，但直到9 0 年代才为人们所重视。l e e p i nm 5 1 ( 1 9 9 2 年) 用有限差分法( f d m ) 求解了2 5 维 矩形电流回线的时间域电磁响应问题；t o r r e s - v e r d i n 等1 ( 1 9 9 4 年) 基于非均 匀介质中电磁散射问题的局部非线性近似，用积分方程法研究了频率域中垂直磁 偶极源的2 5 维正演，但积分方程法涉及到求张量格林函数和散射电流的褶积问 题，只适用于模拟为数不多的局部异常体，否则计算量会成倍增加；u n s w o r t hm j 等盯1 ( 1 9 9 3 年) 采用x 方向布源的“挖源 总场( t m 模式) 算法，实现了频率域电 流偶极源电磁场的2 5 维有限单元法模拟，当异常较小时，计算精度欠佳； e v e r e t tme 等喁1 ( 1 9 9 3 年) 和y ul i m i n g 呻1 ( 1 9 9 4 年) 则在计算洋中脊模型的 电偶极源瞬变电磁响应时，采用二次场算法，避免了上述问题。另外还用g s 变 换作拉氏逆变换，将复数计算化为纯实数计算，提高了计算速度。这些算法有一 一2 一 中南大学硕士学位论文第一章绪论 定的计算精度，但在微机上实现，需要的时间仍然较长( 用8 0 2 5 的网格 e l e r o n 5 5 0 ，1 2 8 m 内存的微机上计算1 个场源1 6 个时间需4 小时) 。m i t s u h a t ay j 阻1 ( 2 0 0 0 年) 提出了一种描述电流偶极源的拟d e l t a 函数法来模拟2 5 维问题的 总场等参单元( i s o p a r a m e t r i c ) 有限单元算法，实现了地形不平的可控源电磁场 模拟，并经数值试验证明了其代码的有效性。由于其对源的近似性，这种算法在 场源附近的计算结果将变得不精确，无法实现近区场源的瞬变电磁响应计算。 我国从2 0 世纪8 0 年代开始着手研究2 5 维电磁场的数值模拟。目前已用二 次场算法实现了二维地电构造上谐变电偶极子电磁场的有限单元算法1 、电偶源 c s a m t 法二维正演的有限单元算法n2 1 、主剖面上时间域和频率域电偶源瞬变电磁 场的2 5 维有限单元正演模拟n 3 儿1 ，时间谱电阻率法的二维正演算法n5 1 ，并对大 粒度一时间级的电磁场数值模拟的微机并行计算进行了探讨们。 针对2 5 维电磁法的正演模拟计算量偏大，精度低，对生产实践中最常用的 重叠回线和中心回线的研究较少的状况，国内王华军( 2 0 0 1 ) 、熊彬( 2 0 0 4 ) 实 现了2 5 维瞬变电磁有限单元法中心回线装置的正演，并为反演做了一些准备。 熊彬( 2 0 0 6 ) 进一步作了电导率线性变化情况下的2 5 维瞬变电磁有限单元法， 其正演都是在普通配置的电脑上实现的，单个测深点运行时间约3 0 0 0 s 。肖明顺 ( 2 0 0 8 ) 对2 5 维磁偶源及电偶源瞬变电磁的理论公式作了详细的推导，并且实 现了带地形的2 5 维瞬变电磁有限单元法j 下演。 本论文的工作首先用f o r t r a n 编程实现中心回线瞬变电磁2 5 维有限元正演， 然后用它对不同赋水情况下的采空区进行数值模拟，并结合现有仪器的精度及一 般情况的地质噪声来讨论瞬变电磁勘察煤田采空区的可行性。 1 3 采空区勘察的现状 采空区问题的普遍性和对人类生命财产安全的威胁，已引起各国的重视。有 关采空区问题的研究，国内外主要是在煤炭、冶金、铁道和军事等部门进行，积 累了大量宝贵的经验；由于采空区具有隐蔽性、复杂性、突然性和长期性等显著 特点，所以至目的为止，针对采空区问题的研究往往是经验介绍，在采空区的勘 察、稳定性评价、治理及其质量监控等方面尚未形成成熟的理论体系。采空区隐 蔽性主要指地下地质条件和受力状况复杂多变，其特征难以弄清；复杂性是指采 空区受多种自然和人为因素影响，其沉降的机理和过程以及对地表、地表上的构 筑物和环境的影响规律非常复杂；突然性是指采空区的失稳和塌陷的时间难以准 确估计预测；长期性是指采空区的活化失稳是一个长期的过程，可能持续几年或 几十年，甚至上百年，它是一个长期缓慢变形过程n 引。针对采空区的上述特点， 3 一 中南大学硕士学位论文第一章绪论 采空区勘察一般遵循如下原则： l 、总原则是以采矿调查为主，有条件时进行井下测量，辅以物探和必要的 钻探；2 、采空区各种探测方法互有优缺点，需互相配合和补充；3 、电法、电磁 法、地震、微重力、放射性等物探技术具有多解性和偏重性，应采用综合物探技 术；4 、单一的勘测手段，往往得出的只是片面的信息，只有把地质、物探、钻 探三者有机结合起来，才能起到事半功倍的效果。 采空区物探方法选用原则：应用于地下空洞探测的地球物理勘探技术种类繁 多，主要有电法、电磁法、地震波法、声波法、重力法等几类。如何从中选取信 息量最大的、最可靠的方法和确定其应用顺序，如何分配各种方法的经费以获取 最大效果就成为首要的问题。各种方法都有其特点、一定的适用条件和应用范围。 因此，下伏采空区的勘察必须根据该类工程的特殊技术要求、场地地质条件及物 探方法的特点与适用条件，选择相应的物探方法及其组合。 采空区综合物探技术探测原理综合物探技术就是应用几种物探方法对同一 工程进行探测，相互配合，发挥各种物探方法的优越性，以获得满意的效果及精 度。目前，采空区勘察常用物探方法的探测原理如下： 1 、电法勘探 电法勘探是以岩土体的导电性差异为基础，研究在外加电场的作用下，地层 中传导电流的分布规律。岩层与岩层之间、岩层与煤层之间、采空区与其围岩之 间的电阻率差异是采用电阻率法探测采空区的地球物理前提。近年来发展较快的 是高密度电法，这种方法的优点是：( 1 ) 探测深度大，可达1 5 0 m ，是地质雷达、瑞 雷波等方法所不能比拟的；( 2 ) 较普通电测深法拾取的信息量大。 2 、地球物理测井 常规测井包括密度测井、声波测井、电阻率测井和放射性测井，以获取地层 参数为主。利用孔内常规测井和钻孔超声成像测井的组合，可以探测采空区及其 上部地层的性质，获耿孔壁周围一定范围内地层的各种物性、电性参数，确定采空 区“三带”在钻孔中的分布发育情况。 特别是超声成像测井是一种较为先进和精密的地球物理测井系统，具有测量 速度快、精度高、方便灵活、随测即显，当场就能看到直观感较强的孔壁图像的 特点。通过对图像及数据资料的分析，便可获得地层产状、钻孔几何形态和其他 物理参量。特别是对无岩芯钻或少芯钻时钻孔质量的检查，地层层理、裂隙的分 析、统计及评价，岩性分辨尤为方便。 3 、地震勘探 地震勘探主要是研究人工激发的地震波在地质体中的传播规律，以探测浅部 地层和构造的分布，或测定岩、土的力学参数特征等。采空区探测中常用的包括 4 中南人学硕十学位论文第一章绪论 s a s w 法、弹性波c t 法和高分辨率浅层地震法等。地震折射波法对探测浅部采空区 效果较好，但因地表干扰影响较大，一般探测深度有限；高分辨率反射波法也是 探测空洞的有效方法，特别是采用横波勘探，可提高分辨率；钻孔弹性波是近年 来新兴的一项技术，可以直观地以剖面形式给出两钻孔间采空区等地质异常体赋 存的状态，从而确定异常范围；s a s w 法探测空洞主要是利用测试曲线的异常( 如 间断、错乱) 来判断采空区的存在，而且其可代替测定k 值，大大减小了波速勘察 的成本。 4 、电磁勘探 电磁法就是通过在地表、井中或空中观测电磁场源( 人工的或天然的) 所产生 的电磁场各分量的变化，达到解决岩土工程问题的目的。电磁法勘探种类较多， 包括瞬变电磁法、地质雷达、频率电磁测深、甚低频电磁法、声频大地电磁法等， 在采空区等空洞探测中常用的方法主要是瞬变电磁法和地质雷达。电磁法对于采 空区大面积普查有很好的效果，可在一定程度上确定空洞的特征，并可探测深部 采空区。 5 、放射性勘探 放射性探测是基于测量放射核素放出的a 、b 、c 射线及中子的一类方法。放 射性探测可分为两大类：天然放射性方法和人工放射性测量方法。目前在采空区 探测中主要使用氡射气勘探，通过a 粒子强度的异常区( 线状异常或环状异常) ，来 判定采空区平面位置和范围，但应注意由于该法定向性不高，而且地表浅层干扰 较大，常会出现较多的异常，故需与临近地段进行比较，主要作为一种大面积普 查的方法。 6 、微重力勘探 微重力勘探的地球物理基础是采空区与围岩之间的密度差和剩余质量。该法 在探测埋深浅、范围小的地质体方面发挥了较好的作用，但由于异常量级和范围 一般都非常有限，因此对其使用要求比较严格。 为了探索采空区勘察的有效物探组合方法，在煤矿采空区的勘察过程中，往 往应用多种物探手段的综合手段。 1 4 研究的目的与意义 瞬变电磁方法对低阻体的响应比较灵敏，且探测深度大，是未来勘察煤田采 空区最有生命力的一种地球物理方法。用瞬变电磁法勘察煤田采空区已有几十年 的历史，取得了一定的成果。但是瞬变电磁跟其他电( 磁) 法一样，在分辨率尤 其是纵向分辨率方面还存在很多问题，野外的好多采空勘察实践多多少少得益于 一5 中南大学硕士学位论文 第一章绪论 采空区的塌陷，相当于采空区规模变大，埋深变浅。在这种情况下我们有必要做 瞬变电磁勘察煤田采空区的可行性分析，以备指导深层采空区的勘察及在地面实 现采空区下覆煤田开采时的探水防水工作。 本文旨在经过数值模拟得到地下采空区的瞬变电磁响应的实测参数a v ，结 合现有仪器的精度及一般情况的地质噪声来讨论瞬变电磁勘察煤田采空区的可 行性。从而避免使用瞬变电磁的盲目性，一定程度上也是在促进并保护瞬变电磁 的发展。 1 5 本文的主要工作及结构 本文的主要工作由以下组成： l 、首先推导、整理了瞬变电磁边值问题和变分问题满足的微分方程，然后 在拉氏傅氏( s ，m ) 域，进行有限单元法计算，由变分把求解二维大地瞬变电磁变 分问题转化为求解多元线性代数方程组，再根据电磁场分量的微分关系求得我们 关心的垂直磁场分量，最后进行拉氏逆变换和傅氏逆变换，最终实现时域瞬变电 磁的数值模拟。这样就完成了论文的程序编写部分。 2 、对所关心的煤田采空区进行2 5 维瞬变电磁数值模拟，分析讨论煤田采 空区在不同赋水情况下的瞬变电磁响应。 论文共分五章，各章内容如下： 第一章绪论：介绍了瞬变电磁法有限元数值模拟的国内外发展现状，介绍 了论文题目的意义及研究目的，最后给出了本论文的主要工作及结构。 第二章瞬变电磁2 5 维正演的基本原理：叙述从麦克斯韦方程组，到有限 元正演所依据的泛函极值问题。推导矩形剖分双线性插值情况下的刚度矩阵，求 变分，把变分方程转化成线性方程组。利用分量之间的微分关系求得我们关心的 垂直磁场分量，通过拉氏逆变换和傅氏逆变换，最终实现时域瞬变电磁数值模拟。 第三章程序实现：介绍了瞬变电磁2 5 维有限单元法正演实现过程中的几 个关键点。 第四章程序验证及采空区模型抽象与数值模拟：用编写的2 5 维程序计算 一维模型，检验程序的正确性和精度。对煤田采空区不同赋水情况进行模型抽象、 数值模拟，进而讨论瞬变电磁勘察煤田采空区的可行性。 第五章结论与建议：通过对几个典型模型的计算结果进行分析，得出一些 结论。 6 中南人学硕十学位论文第二章瞬变电磁2 5 维正演的基本原理 第二章瞬变电磁2 5 维正演的基本原理 王华军，熊彬口0 1 ，肖明顺对瞬变电磁2 5 维有限元正演基本原理，已经 做过相应的分析研究，本文进行了重新的梳理和规整。 2 1 电磁场的偏微分方程定解 瞬变电磁2 5 维正演算法实现的过程为：从时域麦克斯韦方程组及边界条件 出发，经拉普拉斯变换和傅立叶变换转换到拉氏傅氏域，在拉氏傅氏域求解完后， 再经拉氏逆变换和傅氏逆变换转回到时间域。 2 1 1 时域麦克斯韦方程组及边界条件 在三维分区均匀和各向同性介质中，磁 域麦克斯韦方程表述为【1 7 l 【1 9 1 1 2 3 】： v x h ( 垆s 掣+ o r e ( f ) + j 。( ，) a 矿e ( ，) ：一掣+ j ( r ) ( 电) 偶极子形成的电磁场满足的时 ( 2 - 1 ) ( 2 - 2 ) 在介质分界面上满足的边界条件为： n le “( ，) 一e “o ) =(1-0 2 - 3 ) n ih “( r ) 一h 心( r ) l ，= 0 ( 2 - 4 ) 其中，为时间，占、盯、t 分别为介电常数、电导率和磁导率；i i 为分界面f 的法线方向矢量；e 和h 为电场和磁场强度，其上标l 和2 分别表示分界面的两 侧；j 、j 。分别为场源磁流密度矢量和场源电流密度矢量。若场源位于笛卡尔 坐标原点，场源表示为： j 。9 ) = e 万( x ) 艿( y ) 万( z ) o ) ( 2 5 ) j 。( ，) = - a o p m s ( x ) 6 ( y ) 6 ( z ) d f ( t ) 衍 ( 2 6 ) 其中，硒为真空中的磁导率，艿( x ) 、万( ) ，) 、占( z ) 为狄拉克源函数；e = l xi 为源的电流偶极矩；p 埘= s xi 为源的磁偶极矩；l 为供电偶极距矢量；s 为供电 回线面积矢量，指向回线面法线方向；i 为供电电流；z ( ，) 和厶( ，) 分别为电性 源和磁性源随时间嫂化的形函数，在时间域( 瞬变) 电磁场研究中，通常取为上 7 一 中南人学硕士学位论文第二章瞬变电磁2 5 维正演的基本原理 卜舻s 掣卅瓦 7 ， ke s ( 垆一掣“毋) l o t 矿 e ) 一e 乳) f = o ( 2 8 ) 矿 h o ( f ) 一h 乳f ) r = o ( 2 9 ) e ，( f ) = e ( f ) 一e ，( f ) w ：昆掣+ 衍叫：衍e p ( f ) 讲 删= 一酗掣= 。 ( 2 - 1 1 ) ( 2 - 1 2 ) ( 2 - 1 3 ) 其中下标p 和s 分别代表一次场和二次场，昆 - e c p ，衍= 仃一o p ， 印= 一作，计算时忽略了和s p ，和作之间的差异，即认为昆= 0 ，酗= 0 ， 仅考虑电导率差异衍。从( 2 7 ) 式可以看出，纯二次电磁场的场源来自异常体 ( 8 0 - 0 ) 处，且仅与一次场有关，因此，对于中心回线的瞬变响应的数值模拟， 必须计算矩形回线源在均匀半空间中产生的一次电磁场。 8 中南人学硕士学位论文 第二章瞬变电磁2 5 维止演的基本原理 2 1 2 拉氏傅氏域偏微分方程边值问题 为降维( 消除时间变量f ) ，对电磁场作拉氏变换，将时间域i 司趑变成拉氏域 ( 或称频率s 域) 中的问题，然后考虑到二维地电条件下s 、和仃与坐标轴少无 关，沿y 轴对电磁场作傅氏变换，用e ，办表示( s ，所) 域中的电场和磁场，并引入： 锐p 2 e p 。y ，v p2 一i h p yu = e s ，y v = h s 。y 砖= 作s ( o s + ) ，k 2 = 声( 岱+ 盯) c p s + 仃口 m 一, u p s 2 耐2 而i 。j 再可 n 4 p = kp s + 0 1 n s p = 一泌， q = 鬻一南一寿川，a s = - 芦 q2 磊再 吒2 荔砰 口32 磊砰嘞2 够j + d ， 芦 魏= a 2 a 2 p 衍，6 2 一- - - a 2 a 3 j 口衍= a 3 a 2 p 衍，6 3 = a 3 a 3 j 口衍，钆= 衍。 其中肌为y 方向的波数，k 为传播常数。利用矢量运算规则，可以导出仅与 ( s ，m ) 域中电场和磁场沿地电构造走向( y ) 分量有关的l l 和、，的偏微分方程组： 孙o u 拿岛誓+ 也针一 昙( q 瓦o u 一口：象+ 岛鲁一如警) = 吼甜+ 钆 * 知秘鲁一如针 治 鲁f 吩謇+ 呸瓦o u + 玩鲁一如誓 = 呜v 和在与构造走向垂直的( x o z ) 坐标平面上，无一次场场源的介质分界线1 - 两侧，“ 和，满足的边值条件： ( z ，1 - - ( 2 ) ) r = o ( 2 1 6 ) ( v n 一，2 ) r = o ( 2 1 7 ) 9 中南大学硕士学位论文 第二章瞬变电磁2 5 维正演的基本原理 掣鲨o n 智，警+ f 6 f 1 ) 警唧，等 1 锄 l 1 锄 1 锄j = ( 掣世，) 警+ ( 掣世，) 等 掣警智等+ 等攒等) 制，坤) 等+ ( 掣一噬2 ，) 警 ( 2 - 1 8 ) ( 2 - 1 9 ) 式中，上标( 1 ) 和( 2 ) 分别标识分界线r 两侧的参数；r 表不分界线r 的切线 方向矢量。于是，电磁场的计算问题归结为求( s ，m ) 域电磁场j ，分量定解问题 ( 2 1 4 ) 和( 2 - 1 5 ) 式的数值解。在计算出材和，后，可以借助于( 2 - 7 ) 式导出 的关系式： 铲：罢+ ) + ，( 6 2 誓+ 岛鲁 2 。， ( 一q 妄) + z ( b 2a u p 一喀 2 - ， ”( 一吼罢) + 一誓一嗉丁 2 2 ， 札一( q 罢+ ) 一( 6 l 誓+ 喙) 2 3 ， 计算( s ，所) 域中的其他电磁场分量，通常它们均不满足形如( 2 1 6 ) 和( 2 1 7 ) 式的连续性边值条件，故不能直接用有限元法求解。 对于给定的s 和朋值，( 2 - 1 4 ) 和( 2 1 5 ) 式定解问题是二维标量问题，且 为纯实数运算。它的数值解法相对于求解原先的四维矢量问题要简单些。解题思 路如下：分别对若干给定的拉氏和傅氏变换变量s 和朋值，在( s ，m ) 域中计算出 矩形框回线源在空间中的一次电磁场甜p 和，并用有限元方法计算甜和1 ，及巳， 巳，：，忽，和玩，：的数值解，然后经拉氏逆变换和傅氏逆变换便可计算出三维空间 域中申磁场分量的瞬蛮电磁响惠。 一l o 中南大学硕士学位论文第二章瞬变电磁2 5 维正演的基本原理 2 2 拉氏傅氏域二次电磁场的泛函极值 用有限早兀法求( 2 一1 4 ) 和( 2 1 5 ) 瓦疋膦i 叫题阴数值膦。田边但条仟l z 一1 6 ) 和( 2 1 7 ) 式，可导出在有限元法算域( 网格) 外边界g 满足第三类边值条件： 卜q 知乱= 睇+ 如筑( 2 - 2 4 ) 卜吗知鲁l - - 睇+ 包鲁 g ( 2 - 2 5 ) 其中，系数7 7 “和叩”可近似认为是由在均匀半空间条件下，由二次场场源所产生 的电磁场分量来计算。 仁一雌q 新鲁一如剀卢) g 浯2 6 ) 仁一肚知昙+ ( 6 3 鲁气o - 刈l l j g 4 协2 7 ， 由于二次场场源来自地下异常体，计算相当复杂。但考虑到二次场场源本身 比较小( 和一次场相i ：k ) ，在外边界g 上，其数值就更加微小( 甚至可与计算机的 截断误差相当) ，计算时可将其取为零。可用逐次迭代修正办法进行计算，即第 一次计算时设为零，待用有限元求出近似的u 和 ，以后，再来计算7 7 “和7 7 ”，第二 次则在考虑边界条件的情况下再做一次有限元计算，即可求得较为精确的和，如 此反复，直到满意为止。 至此推导出与偏微分方程边值问题( 2 1 4 ) 和( 2 1 5 ) 等价的泛函极值问题 为： ，c 甜，v ，= g q ( 罢) 2 + 口t ( 老) 2 + 口，( 塞) 2 + 口，( 塞) 2 + 以“2 + 呜v + 2 。：罢笔一2 呸老罢“hi 戗u pa 出u - + 2 岛誓罢+ 2 玩誓塞 。2 2 8 ， 啦誓老+ 2 以u p u + 2 如誓罢一2 包誓老+ 2 如誓老 - 2 6 2 誓罢卜砌c ，( r v v 2 + r u u 2 ) 讲 中南大学硕士学位论文第二章瞬变电磁2 5 维正演的基本原理 2 3 变分问题的有限元求解算法 2 3 1 有限元分析 2 3 1 1 区域剖分 对给定的水平地形二维断面，由程序自动剖分有限元网格，剖分的原则为： l 、在介质分界面和异常体范围加密网格，向外逐渐加大间距。 2 、编号以z 方向优先采用蛇状顺序进行，以减少等带宽存贮时的矩阵空间。 用矩形单元对区域进行剖分( 见图2 - 1 ) 瞳4 | ，每个单元的四个顶点为节点， 并给每个单元的电导率盯赋值，矩形单元的宽度为口，高度为b ，单元的四个节 点的编号如图2 - 1 所示： 2 3 + 一n x 个间隔 6 1 0 7ll l ，2 为节点编号 ，为单元编号 4 81 2 图2 - 1节点和单元的排列次序( 双线性插值) 2 3 1 2 双线性插值 取图2 2 ( 8 ) 所示的正方形单元( 母单元) ，四个顶点的编号及其坐标如图上 所示，构造如下的形函数： m = 丢( 1 一孝) ( 1 + 7 ) ， 3 = 三( 1 + 孝) ( 1 一，7 ) ， 或统一写成： j = 丢( 1 + 舌孝) ( 1 + 仇7 7 ) ( 2 - 3 0 ) 其中毒，r i 是点i ( i = l ，4 ) 的坐标。单元中“的插值函数是 1 2 )，： 们 唧 一， + 卜 0n ：，汀从八八v 约 苫 弋 + 卜 亿、n：， 知一4 一 产 中南大学硕士学位论文第二章瞬变电磁2 5 维正演的基本原理 “= ，坼 形函数( 2 - 3 0 ) 式满足 帅) = 骺二 的要求，其中j 代表点号。 将形函数展开： n i = a , 孝r + a e # r + n 蠢q + a 4 其中a i ，6 4 是常数，m 含有勃，孝，7 7 和常数项。 y ( b ) 图2 2 双线性插值的母单元与子单元 ( a ) 母单元( b ) 子单元 图2 - 2 ( b ) 上任意四边形单元( 子单元) 的四个顶点的坐标分别是： ( x i ，y 。) ，( x 4 ，儿) ，函数值是甜l ，一，u 4 ，子单元上的“，x m y 可表示成： 甜= i v , + 2 “2 + 3 坞+ 4 l x = i 毛+ 2 x 2 + 3 x 3 + 4 x 4 y = n l n + n 2 y 2 + n 3 y 3 + n 4 y 4 将形函数( 2 - 3 3 ) 式代入上式，整理后，得 “= 4 孝7 7 + 4 孝+ 4 刁+ 4l x = 骂善7 7 + 岛孝+ b 叩+ 蜀 y = c l 勃+ c 2 孝+ c 3 7 7 + c 4j 其中4 ，c 4 是常数。 令7 7 = c ，这在勃平面上是一条水平线，代入( 2 3 2 ) 式，得： z ，= a m # + a 2 ，x = 6 l 孝+ 6 2 ，y = q 孝+ 乞 1 3 一 ( 2 - 3 1 ) ( 2 - 3 2 ) ( 2 - 3 3 ) ( 2 - 3 4 ) ( 2 - 3 5 ) ( 2 - 3 6 ) 中南人学硕士学位论文第二章瞬变电磁2 5 维正演的基本原理 其中q ，c 2 是常数，所以z f ，x 和y 均是孝的线形函数，当善从一1 变到1 时，“是线性变化的，砂平面上的点也将沿着一条直线变化。同样，令孝= c 时， “，x 和j ，均是r 的线形函数，所以( 2 2 9 ) 式为双线性函数，( 2 3 5 ) 式为双线性 插值，与母单元的四条边对应的子单元的四条边都是直线。 图2 - 3 ( a ) 是母单元，( b ) 是子单元，两个单元间的坐标变换关系为 x ：+ 一a ，y ：y o + 竺(一36)xoy or 23 6x2 + i 亏，少2+ i ( 一 二二 其中，儿是子单元中点的坐标，口，b 是子单元的两个边长，微分关系为 d x = 昙d 孝，咖= 罢咖 ( 2 _ 3 7 ) 2 蛐= 去谢q ( 2 - 3 8 ) ( b ) x 图2 3 矩形单元剖分、双线性插值的母单元与子单元 ( a ) 母单元( b ) 子单元 2 3 1 3 单元分析 当同一单元中的口，0 = l 5 ) 和6 。( 6 = l 4 ) 是常数时，“对x 的偏导： 罢李，= ( 封u e = u e t ( 罢) 其帆哳“，= ( 挚譬i a a , 卜所以 “ ( 2 = u e t ( 等) ( 篆) r u 。 将( 2 - 2 8 ) 式中的区域积分分解为各单元积分之和： ( 1 ) 单元积分 1 4 ( 2 - 3 9 ) ( 2 - 4 0 ) 中南大学硕士学位论文 第二章瞬变电磁2 5 维正演的基本原理 舡( 罢) 2 + ( 老) 2p q = u ；k ：，。u 。 其中k ：，。= ( 吒) ，k = 屯 ( 2 - 4 1 ) 二； e 鐾圣；裂然鲁) 号 d f d 力c 2 4 2 ，= 卅( 等警) ( 考警) + ( 筹鲁) ( 等鲁) 等d f 却 对( 2 3 7 ) 求善和7 7 的微商，代入上式积分后，即可得屯，具体计算公式如 毛l = 2 a + 毛2 = 屯l 恕l = a 一2 p丘：= 岛。 岛l = 一口一毛3 = 墨。 也l = - 2 a + 心3 = 如l 乞：= 向。k = 墨。 其中口= 石b ，= 盖 ( 2 ) 单元积分 舡( 削册疆。v 。 其中畦，= k ：，。 ( 3 ) 单元积分 肌2 d q = v t2 。 其中。= ( 乃) ，k = k 下： 铲。，n , n ja ，b d n ( 2 - 4 3 ) ( 2 - 4 4 ) ( 2 - 4 5 ) ( 2 - 4 6 ) 对( 2 3 7 ) 求f 或7 7 的微商，代入上式积分后，即可得毛，具体计算公式如 1 5 中南大学硕士学位论文第二章瞬变电磁2 5 维正演的基本原理 墨l = 4 a 乞l = 2 口 毛l = 口 缸l = 2 a 岛：= 墨。 其中口= 等 ( 4 ) 单元积分 小2 d q = v t2 ，v v 。 屯：= 乞。 屯：= 毛。 毛，= 毛。 缸，= 毛。 = 毛。 其中略= k 删2 ( 5 ) 单元积分 舡罢塞一老塞p q = t 属3 以几， 其中k 叫3 = ( 毛) ，k = 一乃 吒= ( 百a n , ) a n , j - 剐i 2j 似l - - d - = ( 盟d e 堂、( 生d r i 塑d 2 ) ( 型d , 7 却d z ) ( q d e 蟛 愕t d f 却 x f f ( 2 3 7 ) 求孝和r 的微商，代入上式积分后，可得气 毛l = 0 k 2 l = - 0 5 岛l = 0 颤t = 0 5 屯：= 毛 岛：= 岛。 屯：= 毛。 岛，= 毛， 心，= 也， k = 毛 ( 2 - 4 7 ) ( 2 - 4 8 ) ( 2 - 4 9 ) ( 2 5 0 ) ( 2 5 1 ) 为了使总的刚度矩阵对称，在形成总的刚度矩阵时需要对其做适当变换。 同理可得： ( 6 ) 单元积分 舡誓罢毒罢卜啦t 。，u 其中k i 脚= k ：，。 ( 7 ) 单元积分 一1 6 ( 2 - 5 2 ) 中南大学硕十学位论文 第二章瞬变电磁2 5 维正演的基本原理 虹誓塞+ 誓珈吨t 咄儿，p 其中k k ，= k ：，。 ( 8 ) 单元积分 m p 谢q = u 。t k 删2 p 其中k 乙，。= k ：，。 ( 9 ) 单元积分 盟誓罢一誓老卜咄：以p 其中k ：，唧= k o ( 1 0 ) 单元积分 瓤誓老一誓妻卜啪t 。， ( 2 5 3 ) 其中k ：，础= 一儿， ( 1 1 ) 单元积分 4 ( 7 7 。，2 + 7 7 “材2 ) d ，对于二次场算法而言，二次场场源是地下的异常体，对 _ 、 a q 边界的积分是相当复杂的，好在是异常体相对我们研究的区域较小，在外边界上 二次场已经衰减的很微小。正规的做法是先设为零，经有限元计算后加以修正再 进行有限元计算，反复多次至到满意为止。本文在计算自己的模型时将其取为零， 而在