（地球探测与信息技术专业论文）基于ansys直流电测深法的边坡勘探正演.pdf

中南大学硕士学位论文 摘要 摘要 直流电测深法边坡勘探时，电测深曲线的异常并不仅仅是边坡的 隐患所引起的，这就对电测深资料的解释带来了很大的困难，为了更 好的了解边坡隐患在电测深曲线上的异常特征，本文用a n s y s 软件对 这一问题作了探讨和研究，实现了电测深法对边坡地电模型的有限元 数值模拟，并与相应的物理模拟结果进行对比分析。 本文首先回顾了电法勘探正演数值模拟的研究现状，以及前人对 边坡勘探所做的研究；总结了电法勘探正演问题现有的数值模拟方 法，简单介绍了a n s y s 有限元软件的主要功能，在进行直流电测深法 边坡勘探数值模拟过程中的主要步骤及其注意的关键问题。然后介绍 了点源场电位边值、变分问题；通过使用a n s y s 本身提供的a p d l 语 言编制相应的计算程序，以几种典型的理论边坡模型为例，进行了一 系列的直流电测深地电模型的数值模拟计算，得到了一些边坡隐患的 电位分布规律，结合物理模拟进行了对比和验证，结果对今后进行边 坡电测深勘探具有很好的参考价值。 关键词：电测深；数值模拟；边坡；有限元； 中南人学硕士学位论文a b s t r a c t a b s t r a c t w h e nw ee x p l o r et h es l o p eb yr e s i s t i v i t ys o u n d i n g ，t h ea b n o r m i t yo f t h e s o u n d i n gc u r v i l i n e a rb yt h es l o p eh i d d e nt r o u b l e i ti st a k em u c h d i f f i c u l tt oe x p l a i nt h e r e s i s t i v i t ys o u n d i n gt oe x p l a i nt h er e s i s t i v i t y s o u n d i n gd a t ei nt h ei n t e r e s to fr e a l i z et h ea b n o r m i t yo ft h es l o p eh i d d e n t r o u b l ei nt h es o u n d i n gc u r v e t h ep a p e rt od i s c u s sa n dr e s e a r c ht h e q u e s t i o nb yn u m e r i c a ls i m u l a t i o n ，a n dc a r r yo u tt h ef e mt oi m p l e m e n t n u m e r i c a lm o d e l i n go f r e s i s t i v i t ys o u n d i n gb ya n s y s s o f t w a r e t h e nw e c o m p a r e di tw i t ht h er e s u l to fp h y s i c sm o d e l i n g f i r s t l y , t h ep a p e rr e v i e w s t h er e s e a r c ha c t u a l i t yt h e t o p o g r a p h y i n f l u e n c ei nt h ep r o s p e c t i n go fs l o p e ，a n dt h ef o r w a r dr e s e a r c ht h a tp e o p l e h a v ed o n e i tc o n c l u d e dt h es e v e r a lk i n d so fm e t h o d si nt h ee l e c t r i c s o u n d i n gp r o s p e c t i n g ，a n ds i m p l yi n t r o d u c e dt h em o s t l yf u n c t i o no f a n s y sf e mp r o g r a ma n dt h ec e n t r a l s t e p s i nt h e p r o c e s s i n g o f n u m e r i c a lm o d e l i n g t h eb o u n d a r yv a l u ea n dv a r i a t i o ne q u a t i o no ft o t a l p o t e n t i a lo np o i n ts o u r c ee l e c t r i c a lf i e l do fg e o e l e c t r i cs e c t i o nw e r ea l s o i n t r o d u c e d t ou s et h ea p d ll a n g u a g e ，w h i c ho f f e r e db ya n s y si t s e l f , w ew o r k e do u t c o r r e s p o n d i n gc o m p u t i n gp r o g r a m i nt h et y p i c a lt h e o r e t i c s l o p em o d e l ，w ec a r r i e do u tas e r i e so fn u m e r i c a lm o d e l i n ga n df o u n do u t t h er u l eo ft o p o g r a p h yi n f l u e n c ei nt h ep r o s p e c t i n go fs l o p e i ti s c o m p a r e db yp h y s i c sm o d e l i n g s e c o n d l y , w em a k et h ee r r o ra n a l y s i sb y t h er e s u l to fp h y s i c sm o d e l i n ga n dn u m e r i c a lm o d e l i n g f i n a l l y , w em a k e ac o n c l u s i o no ft h ep a p e rt h a th a sp o s i t i v em e a n i n gi nt h ep r a c t i c ew o r k k e y w o r d s ：r e s i s t i v i t ys o u n d i n g ； n u m e r i c a ls i m u l a t i o n ；c o a s t ； i n f i n i t ee l e m e n t ； l l 原创性声明 本人声明，所呈交的学位论文是本人在导师指导下进行研究工 作及所取得的研究成果。尽我所知，除了论文中特别加以标注和致谢 的地方外，论文中不包含其他人已经发表或是撰写过的研究成果， 也不包含为获得中南大学或其它单位的材料。与我共同工作的同志对 本研究所做的贡献均已在论文中作了明确的说明。 作者签名：塑魍日期：z 丑年上月j 竺且 关于学位论文使用授权说明 本人了解中南大学有关保留、使用学位论文的规定，学校有权保 留学位论文，允许学位论文被查阅和借阅；学校可以公布学位论文的 全部和部分内容，可以采用复印、缩印或其它手段保存学位论文；学 r 日 一j 中南大学硕士学位论文第1 章绪论 1 1 引言 第1 章绪论 边坡是自然或人工形成的斜坡，它的形成与演变，与人类生产、活动有着极 为密切的关系。作为全球性三大地质灾害之一的边坡失稳破坏已严重危机到人们 的生命安全和财产安全。目前对边坡地质灾害的勘探手段有很多，地球物理方法 是探测边坡隐患的一种既经济又实用的方法。直流电测深法是电法勘探中应用最 为广泛的一种方法，在野外勘探中具有独特的优势。但是在直流电测深法进行边 坡勘探时，由于地形、跑极方向等因素的影响，使得电测深曲线的形态发生了变 化，也就是产生了所谓的假异常，这给资料解释带来了极大的困难。因而我们有 必要针对边坡勘探时边坡隐患的电测深曲线异常形态做些研究和讨论。 1 2 边坡隐患勘探的研究现状 依据滑坡体组成物质的不同可将滑坡分为土层滑坡和基岩滑坡两大类，其 中，土层滑坡包括粘土滑坡、黄土滑坡和坡积土滑坡；基岩滑坡包括破碎岩体 滑坡和完整基岩滑坡。不同种类的滑坡具有不同的地球物理条件，只要滑坡与 滑动面下部的岩( 土) 之间存在明显的波速或电阻率差异，便可采用地震勘探 和电法勘探进行探测。用于边坡失稳调查常用的地球物理方法有 4 - 1 7 体2 】：电 法、浅层地震，探地雷达、地热法、放射性法、重力法、声波法和土壤气体测 量法。特别是雷达法、声波法和土壤气体测量法是近几年来使用的最新方法。 对于边坡失稳的物探方法选择，要针对具体问题进行选择。赵安宁 s 5 1 等通过某 水电工程边坡岩体卸荷厚度探测的应用实例，说明瞬态瑞雷波勘探技术在查明边 坡岩体松动、卸荷程度和深度分布是一种便捷和行之有效的方法；吴向明等【5 7 】通 过实际的滑坡勘察工作，介绍了采用综合测井方法在确定滑坡滑面位置和滑坡体 厚度上的应用；李大兴口卅介绍探地雷达利用脉冲电磁反射波的波形特征了解了 地下介质的分布；高志亮【5 4 1 等通过研究得出环境地球物理在滑坡勘探中也是有 效果的；金为民等陋引利用地震反射波法，地震折射波法等浅层地震勘探方法也 能找到滑坡的滑动面。 如果差异较小则应用物探方法探测就比较困难。如滑移式边坡失稳，用重 力测量圈定滑坡范围，测量滑坡体下滑推力以及滑动面的状况；用电阻率法或 地质雷达圈出滑坡范围和深度以及确定滑坡带的位置、自然电场监测滑坡动态、 电磁波测量监测滑坡动态；用地震法确定滑坡范围、深度以及滑动界面状态， 中南大学硕士学位论文 第1 章绪论 利用微动观测监测滑坡体的震动频谱，确定滑坡滑动方向和滑动面运动等；用 地温测量与滑坡有关的地下水动态；用放射性法测量断裂及滑动面状态，采用 氡气测量判断坡度的稳定性，圈定滑坡体并监测滑坡发展的过程。如唐深海运 用对称四极电测深法对川西高原康定白土坎滑坡进行勘探，查明了滑体结构， 滑面埋深及地下水分布状况；杨宝健采用g d t 高分辩率地质探测仪和g d s 高分 辩测深法在高阶地前缘斜坡地区探测黄土滑坡，取得较好效果；刘崧运用电法、 地震法、重力法等对山峡库区滑坡进行研究，确定了滑体外边界、划分了不同 的物质结构区，以及确定滑床形态、埋深；夏望林，王国利运用浅层地震反射 波多次覆盖技术对万州区和平广场滑坡群进行勘查，查清了滑体堆积层厚度及 其展布情况，揭示了覆盖层厚度及所测基岩面埋深；叶树林，黄国夫用放射性 勘查方法对广元市境内一老滑坡进行勘查，大致圈定了滑坡边界，控制了滑坡 的滑动方向。王谦身等运用微重力法对紫阳城滑坡体进行探测，根据微重力的 探测结果，推算出了滑坡体的边界和浅层应力场；吴庆曾在深圳特区环境地质 调查中，对岩体软弱层和破碎带进行声波测试，取得良好效果；原地矿部水文 地质工程地质技术方法所先后在重庆万洲关塘口滑坡、巴东县新址滑坡体进行 声波测井勘探中，划分了岩性、评价了其完整性，确定了滑带、破碎带位置。 航测和卫星照片是目前世界上用于判断滑坡规模的先进的、可靠的方法。我国 在不少的滑坡观测上，用到了该项技术从照片上反映的地貌形态、水文的分布、 色调及植被等即可判断出滑坡的范围和大小；若辅以红外线的扫描，还可以判 断出滑坡体的细节。在日本，利用遥感测量技术，成功地对奈川渡水库进行调 查，圈出滑坡或崩塌达4 1 9 处。因此，对边坡失稳地质灾害隐患探测方法的研 究也是一个不断探索、发展的过程。 1 3 电阻率法正演模拟的研究现状 与其它地球物理方法的发展一样，电阻率法的研究也经历了三个过程的发 展；求解析解、物理模拟试验和数值模拟计算。( 1 ) 求解析析解是指对地下比 较规则的地质体，如球体和板状体、水平层，可以通过解析法推导出电场的解析 表达式，该方法得到的结果很精确，但公式的推导过程非常复杂，应用范围不广， 对于一些比较复杂的场源分布或不规则地质体，解析解很难得到，不能满足实 际生产的需要，一般用来验证其他正演方法的j 下确性。( 2 ) 物理模拟法，也叫 模型试验法，主要使用各种物理模拟设备如水槽、土槽、电阻网络、导电纸和薄 水层等方法，它们均以物理现象的相似性味基础。对于稳定电流场，物理模拟的 2 中南人学硕十学位论文第1 章绪论 相似性准则比较简单：在模拟地电断面时，模型的几何尺寸和电阻率都可分别按 一定的比例尺进行模拟、缩放。通过测量各种模型的场值得分布与大小，便 可以得到正演模拟结果。不过，由于模型制作的困难，有时难以模拟比较复杂的 模型，因此该方法应用不广。( 3 ) 数值模拟计算，根据地球物理中的偏微分方 程和边界条件，用数值计算的方法求解场值的近似解。虽然这是一种近似的方法， 但它适用于复杂物性分布和复杂边界形状的地球物理模型求解问题，所以适用很 广，而且计算速度很快，因此已经成为现代地球物理研究的有力工具。 解析法和物理模拟法经过大半个世纪的发展已经比较成熟和完善，但对新出 现的复杂地球物理问题却束手无策；数值模拟计算是一种有效途径，近数十年来 在正演模拟方面的发展很快。下面就电阻率法的数值模拟发展的情况作以下介 绍。 电阻率法正演问题的数值计算方法一般可归纳为四种：有限差分发法、有限 单元法、积分方程法，这四种方法各有优缺点： ( 1 ) 有限差分法( f i n i td i f f e r e n c em e t h o d 或f d m ) 。这是一种经典的数 值模拟计算方法，其基本原理就是用差商代替微商，将待求解的连续微分方程变 换为离散的差分方程，并通过求解差分方程得到原微分方程的近似解。 第一次将有限差分法利用到电磁法领域是在6 0 年代末期，j e p s e n ( 1 9 6 9 ) 讨论了二维地电条件下点电源电阻率法和激发极化法问题，l a m o n t a g n e 【1 1 等 ( 1 9 7 1 ) 对矩形薄板的电磁法响应进行了数值模拟，s w i f tm 1 对二维非均匀介质 的大地电磁的理论响应进行了分析计算；j e p s e n 3 1 在博士论文中( 1 9 7 5 ) 对有限 差分在电阻率正演模拟中的应用做了详细论述，m u f t i 3 1 ( 1 9 7 6 ) 改进了有限差 分网格剖分方法，采用非均匀网格计算任意形状二维体取得了很好的效果。7 0 年代木期，d e y 和m o r r i s o n 5 1 ( 1 9 7 9 ) 将混合边界提条件引入有限差分法中，对任 意三维形状的地质体的电阻率进行了正演模拟，s c r i b a 【6 1 在1 9 8 1 年对三维断面 的电场响应也进行了论述。z h d a n o wc 7 1 等( 1 9 8 2 ) 对三维地质体的电磁响应进行 了数值模拟。l e e p i n ，m 8 1 ( 1 9 9 2 ) 用有限差分法求解了2 5 维矩形电流回线的 时间域电磁响应问题；s p i t z e r 郫等( 2 0 0 1 ) 在有限差分的离散化和收敛速度方 面进行了深入研究，对三维正演模拟的精度和速度进行了分析。 我国的地球物理工作者从8 0 年代开始研究有限差分法j 下演计算问题，如周 熙襄等 1 0 1 ( 1 9 8 3 ) 、罗延钟等i j ( 1 9 8 4 ，1 9 8 6 ) 为有限差分法的发展做出了相 当大的贡献；刘树才等2 1 ( 1 9 9 5 ) 、刘正栋等3 1 ( 2 0 0 0 ) 在差分格式、网格剖 分方法等方面做了很多工作。刘树才等 1 4 1 ( 2 0 0 4 ) 对电法三维正演中的若干问 题作了详细论述，他指出有限差分法和有限单元法是解决三维问题的理想方法； 吴小平，徐果明等卅( 1 9 9 8 ，1 9 9 9 ) 基于有限差分原理利用共轭梯度法求解了 3 中南大学硕士学位论文第1 章绪论 点电源三维正演模拟问题。刘志新 1 6 1 等( 2 0 0 2 ) 、刘树才【1 7 j ( 2 0 0 4 ) 对矿井中 三维直流电阻率法的正演算法作了改进，使原有的算法能适应条件。 ( 2 ) 有限单元法( f i n i te l e m e n tm 或简称为f e m ) 是一种求解场的变分问 题的数值方法。用这种方法求解稳定电流场电位时，首先要建立给定边值条件下 电位的微分方程，并通过变分原理，等价地变换成扁分方程；然后将求解区域离 散化，即按一定的规则将求解剖分成相互连接的网格单元，导出单元节点点电位 的高阶线形方程，集成整个区域单元的节点方程，解方程组可求出所有节点电位。 该方法5 0 年代首次在力学领域发展起来的一种模拟方法，最初用于结构力 学和应用力分析( 0 c 齐基威茨，y k 邱，1 9 6 7 ) u s ；1 9 7 1 年，c o g g o n b g 发 表了著名论文拉丌了有限元法在电法勘探领域正演模拟的序幕。r o d i 2 0 1 ( 1 9 7 6 ) 发展了有限单元法的剖分方法，采用矩形网格剖分，以解决二维大地电磁测深正 演问题，使有限元向前发展了一步；r i j o 2 h ( 1 9 9 7 ) 引入了一个通用性网格剖 分方法，使有限元正演的精度和速度得到了大幅度提高，成为计算二维地电条件 下电法j 下演模拟的有效工具，使有限单元法正式进入实用阶段。7 0 年代末，在 甚低频法正演( p k a i k k o n e n 2 2 1 ，1 9 9 7 ) 、时间域电磁法的正演模拟( j o n et k u o e ta 1 ，1 9 8 0 ) 【2 3 】中得到应用；d e ya n dm o r r i s o n 【2 4 1 ( 1 9 7 9 ) 将混合边界条件引 入到有限单元法的正演模拟中使其得到进一步发展；1 9 8 1 年，p r i d m o r e 【2 5 1 等对 有限单元法的三维地电断面的电法正演问题进行论述，w a n n a m a k e r 等 2 6 1 ( 1 9 8 6 ) 用有限元法模拟大地电磁中的二维地形响应。u n s w o r t h m j 等1 2 7 j 发表了频率域 电流偶极源电磁场的2 5 维有限单元模拟； 7 0 年代木有限单元法引入我国( 朱伯芳，1 9 8 6 ) 2 s 1 ，8 0 年代初，我国数学 家李大潜发表专著有限元素法在电法测井中的应用【2 9 1 ，有限单元法正式应 用到电法领域，从此地球物理学家周熙襄m 】1 3 t 1 ( 1 9 8 0 ，1 9 8 6 ) 、罗延钟【3 2 1 ( 1 9 8 7 ) 和徐世浙 3 3 1 【3 4 1 ( 1 9 8 2 ，1 9 9 4 ) 等对有限单元法进行了深入研究，发表了一系列 有限单元法在电阻率法中的应用的论文或专著，研究领域涉及直流电法，电磁法 领域，网格剖分也由原来的简单剖分发展到三角形单元和三角矩形综合剖分，周 熙襄等 在二维电阻率有限元法正演中采用混合边界条件等优化措施，提高了计算精 度和速度，对有限单元法的发展起到了推动作用；从八十年代开始着手研究2 5 维电磁场的数值模拟d s 1 3 6 ，到目前，已用二次场算法实现了二维地电构造上谐 变电偶极子电磁场的有限单元算法【3 7 】、电偶源c s a m t 法二维正演的有限单元算 法【3 8 】、电偶源瞬变电磁场的2 5 维有限元正演模拟 3 9 1 【4 0 】 4 h ，时问谱电阻率法 的二维正演算法【4 2 】；阮百尧等【4 3 】【“1 ( 2 0 0 1 ，2 0 0 2 ) 利用有限单元法在个人微机 上实现了三维地电断面的j 下演模拟，开启研究三维地质构造的新起点。 4 中南人学硕十学位论文第1 章绪论 ( 3 ) 积分方程法( i n t e g a le q u a t i o n ) 是从场参数( 场强或场位) 所满足 的微分方程边值问题出发，通过某些变换导出有关参数( 场强、位、积累电荷密 度或散射电流密度) 所满足的积分方程。然后用近似计算方法求此积分方程的数 值解。 l a i f o n o 4 5 1 ，k v o z o f f 4 6 1 及g v k e l l e r 和f c f r i s c h k n e c h t 【4 7 1 等，早 在五十年代末到六十年代初就提出了用积分方程法计算任意形状构造上的视电 阻率。后来，k d i e t e r 【4 8 l 等于1 9 6 9 年提出具体用积分方程法计算了三维体的激 电和视电阻率异常曲线。1 9 7 3 年a d e y 【4 9 i 等从麦克斯韦方程出发，建立了更普 遍的积分方程，把积分方程的应用范围推广到计算电磁法异常。 在我国电磁法领域最早应用积分方程法解决电法数值模拟问题的学者有周 熙襄”0 1 5 q 、罗延钟【5 2 1 等。随着积分方程的发展，出现了边界积分方程，在求解 过程中借用有限单元法的离散技术，进而发展为边界单元法。 ( 4 ) 边界单元法( b o u n d a r ye l e m e n tm e t h o d 或简称b e m ) 。边界单元法的 前身就是边界积分方程法，随着有限单元法的兴起，其单元划分和插值函数的概 念引入到了边界积分方程法中，发展成为边界单元法。边界单元法最早是由英国 的s o u t h a m p t o n 大学土木工程系开始使用，并逐渐应用到各个领域。边界单元法 的优点是只需要在电性异常体的边界上剖分单元，使求解问题的维数降低，最后 形成的代数方程组得元数大为减少，尤其对于三维问题和无限区域问题优点很突 出。 这四种数值模拟方法各有其优缺点，有限差分法方法简便易算，但是当物性 参数复杂分布或是场源的几何特征不规则时，适应性比较差。边界单元法正演速 度快，内存需求少，主要用于地形改正和地下少量地质体的正演模拟。有限单元 法适用范围比较广，但是其缺点是计算量大，计算时间长。 1 4 研究的主要内容 一绪论：主要介绍了论文的研究意义以及边坡勘探和电阻率正演模拟的研 究现状。 二数值模拟方法的基本原理：首先介绍几种常用的数值模拟方法，然后着 重讲述了点电源2 5 维，点源三维有限单元法的基本原理。 三有限元模拟在a n s y s 中的实现：详细讲述了如何应用a n s y s 软件实现电 有限单元法模拟的步骤和注意事项。 四边坡的地球物理背景及模型计算；从边坡的类型和边坡具有的独特地球 物理性质阐述了直流电测深法勘探边坡隐患的可能性。依据边坡的类型和地球物 理性质构建边坡模型加以计算，并与物理模拟的结果对比，验证数值模拟的可靠 5 中南大学硕士学位论文 第1 章绪论 性。 五结论和建议：总结了通过数值模拟归纳的直流电测深法边坡勘探时，电 测深曲线的异常形态。 6 中南人学硕七学位论文第2 章正演数值模拟方法的基本原理 第2 章正演数值模拟方法的基本原理 2 1 电阻率法正演模拟的方法 物理模拟方法，虽然比数学解析法更易于解决比较复杂地电条件的问题， 但仍有其缺点。比如，对比较复杂的地电条件构建模型也十分复杂，此外，实 验结果受电位差观测误差的影响，往往精度不高；最后，物理模拟的效率较低， 并且实验结果不便于用于自动化数据处理。现代数值模拟方法可以避免这些缺 点。 电场的数值模拟，就其原理来说，是早已解决了的问题；但由于其计算工 作量大，所以实现比较困难。近年来，由于电子计算技术的飞速发展，为对电 场做快速和精确的数值模拟提供了技术条件，使之得以实现。目前，用来对稳 定电流场作数值模拟的方法主要有以下几种，即有限差分法，有限单元法，边 界单元法，积分方程法。 2 1 1 有限差分法 有限差分法是从电场( 电位) 或电磁场所满足的偏微分方程和边界条件出 发，将微分方程转变为差分方程，其研究步骤是：首先将研究区域按一定方式 离散化，然后在每个单元内设位、场呈线性变化，电性为均匀的，因而微分方 程的微分就可用差分来代替，于是就可以建立一组线性差分方程，最后求解此 线性方程组即得相应的位场分布。 通常解决稳定电流场空问分布的方法是根据电磁场理论导出一定的微分方 程，并用解析法求解；有限差分法则是将导出的微分方程离散化，得到相应的 差分方程以求解。具体来说，在稳定电流场中，电流密度，、电场强度e 、电 位v 和电导率盯之间有如下关系： ，：仃e ，e ：一v y ( 2 1 ) 即： ；：一仃v y ( 2 2 ) 在电流源以外的点上，根据电荷守恒定律有： v ( 仃v 矿) = 0( 2 - 3 ) 或者经过矢量运算，将上式改写为： 7 中南人学硕士学位论文 第2 章正演数值模拟方法的基本原理 v 2 ( ) + o - v 2 v 一网2 盯= 0 ( 2 4 ) 上式便是在非均匀介质中，稳定电流场电位v 满足的微分方程。在均匀或分区 域均匀介质中，它简化为众所周知的拉普拉斯方程：v 2 = o 。 在用有限差分法求微分方程( 2 - 4 ) 的数值解时，首先要在求解区( 或模拟 区) 建立网格，以便将微分方程离散化。对于模拟区( 即网格) 边缘上的节点 ( 外边界点) ，其电位值可根据具体的地电条件给出，或者表示为该点坐标的已 知函数，或者表示为相邻节点电位的线性函数。这样，网格上的每一个节点都 有一个表示该点电位与相邻点电位之间关系得线性方程，总合起来，便得到一 个未知数( 即各点的电位值) 和方程个数等于网格总节点数的高阶线性方程组。 现代计算技术有多种方法可以用来解这个方程组。由此，便可算出网格中各节 点的电位值，即求得模拟区内空问各点的电位分布。 2 1 2 有限单元法 有限元法是从位、场所满足的偏微分方程出发，根据微分方程解与泛函极 小问题的等价性，将微分方程和其边界条件转化为相应泛函的变分问题，其研 究的步骤是将研究区域按一定方式离散化，设单元内位、场呈线性变化，电性 参数均匀，这时泛函是各节点位、场的二次函数，利用求极小的必要条件，即 泛函对各节点位、场的变分为零，二次函数的变分为一次函数，由此得到一个 线性方程组，解此线性方程组便可得各节点的位、场值。 有限单元法是近二十多年来随有限差分法之后迅速发展起来的又一种有效 的数值模拟方法。七十年代初，h j 科冈首先将其用于电法模拟。以后其他一些 学者进一步完善了它，使之成为求解复杂条件下电阻率法问题的有力工具。 有限单元法求解稳定电流场电位v ，是根据变分法原理或最小总电能原理 将给定边界条件下求解( v ) 微分方程的问题，等价地转换成求v 的泛函的极 值问题；再经过离散化，得到由空间各点未知电位值组成的高阶线性方程组； 最后在数字电子计算机上求解该方程组，以确定空间各点的电位。 在包括电流源在内的不均匀三度空间中，稳定电流场的电位v 满足泊松方 程： - 9 , v ( c r v y ) = 一v 丘 ( 2 5 ) - 0 式中口为介质电导率；v 为稳定电流场电位；v t 一电流源体密度。它们 都是空间坐标的函数。在直角坐标系中，公式( 2 5 ) 可写成： 中南大学硕士学位论文第2 章正演数值模拟方法的基本原理 兰( 仃尝) + 晏p 豢) + 兰( 盯娑) ：一v z ( 2 _ 6 ) o ) co x钟铆d zo z 根据变分原理，若在区域q 的一部分边界s 上，v 已知；在其余边界s ，上， a t 厂i v 满足“自然边界条件”：l，其中n 是s ：的外法线方向。则求解微分方程 咖i j 2 ( 2 6 ) 可等价转换为泛函极值问题，即找出一个函数v ，使其泛函j 取得极小值： ，= 皿主 仃 c 尝，2 + c 芳，2 + c 尝，2 一y v z ，蛐出= 极小 。2 7 ) v l s 。= 以知：s 。+ 最为q 的全部边界 有限单元法和有限差分法是根据不同的原理作数值模拟，实践证明，两种 方法可以得到彼此相近的模拟结果。 有限差分法的程序编制比较简单，采用这种方法对平面规则形体( 如板状 体、层状体等) 进行数值模拟比较方便；但它对内节点、内边界点及两种类型 的外边界点要分别用不同的电位计算格式。所以，当不同电导率的介质较多或 不均匀体表面积为曲面时，使用这种方法不甚有利。 有限单元法易于模拟具有不规则表面的地电模型；此外，这种方法在整个 求解区可采用同一种计算格式，因而可避免采用多种计算格式时给模拟结果带 来的误差。但有限单元法程序编制较复杂，占用计算机的内存量也较大。 本文所用的地电模型主要是有地形的边坡模型，故采用有限元法进行实现。 2 1 3 边界单元法 一以加权剩余格式建立边界积分方程 假设q 是一个域，它的边界r 被分成两部分1 1 。和r ：，总边界r = r 。+ 1 1 ：位函 v 2 u o = f 数所满足的控制方程为本质边界条件： 1 1 0 ：s 自然边界条件：g ：誓：g 。 o n ( 2 8 ) 式中为场位的精确解，s 和g 为定义在f 。和1 1 ：上的已知函数。现以加权剩 余格式由控制方程及边值条件出发，建立边界积分方程。我们使用以下符号 9 中南大学硕十学位论文第2 章正演数值模拟方法的基本原理 ( 2 - 9 ) 上二式分别表示函数u 和v 的乘积再域q 和边界r 的积分。称( “，v ) q 和 “，d r 为u 和v 的内积。 首先定义l ( u ) 与权函数的内积 q + ( s ( u ) - - s ，w r 。+ ( g ( “) 一g ，w ) r ：= 0 ( 2 1 7 ) 如果在域q 及其各边界段分别设权函数w 为 r l ( z ，y ，z ) q w = g ( w 1 ) ( x ，y ，z ) e f l ( 2 一1 8 ) i s ( 2 ) ( z ，y , z ) r 2 利用式( 2 1 3 和( 2 - 1 8 ) ，可将( 2 - 1 7 ) 式写成： ( 厂，) n = ( “，( w ) ) q + ( s ( ) ，g ( “) ) r 一( s ( “) ，g ( ) ) r ( 2 1 9 ) 在边界元法中，采用算子的基本解作为权函数w ，即采用满足如下方程的奇异 1 0 、rj 【 ， 别 讲 ，如p = i l c = r y y m “ 中南人学硕十学位论文第2 章正演数值模拟方法的基本原理 函数 三( 矿) + 4 = o ( 2 2 0 ) 式中点是狄拉克函数，它在所考虑的“i ”点的积分等于1 ，而在其他地方的 积分为零，且有 = uf ( 2 2 1 ) 将( 2 - 2 0 ) ，( 2 - 2 1 ) 代入方程( 2 - 1 9 ) 中，即得： c i “，= - n 一 r ( 2 2 2 ) 式中q 是与q 的边界所确定的张角成正比的系数。 若记e = ( f ，) n 则( 2 2 2 ) 式可写成 c i “，+ ( s ( “) ，g ( w ) ) ，= ( s ( w ) ，g ( “) r b ， ( 2 2 3 ) 公式( 2 - 2 3 ) 就是待求得边界积分方程。 二确定离散方程与基本方程 为了利用边界积分方程( 2 2 3 ) 计算待求电位u ；需要将其离散化为线性代 数方程组。为此将边界r 划分成许多的小单元，称之为边界元。 因为线性元的u 和q 在单元上呈线性变化，故在两个直线元的交点( 端结 点) 上赋予己知值。假定在r 。的端结点上，q 是己知的，在r ：的端结点上u 是 已知的。根据积分的可加性，对于结点“i ”，将式( 2 - 2 3 ) 中的积分化为每个 单元1 1 ，的积分之和，即得离散方程： 掣，+ 二nf fg ：塞正，g 钉一b ，(224)j u i j u 对于平面域的边界元r ，因为u 和q 在单元上是线性变化的，所以u 和q 在单 元上任一点的值，可借助予端结点的值和两个线性插值基函数来确定： 甜( f ) 2 仍“- + 仍“z2 【仍仍】 ：， 2 2 5 ) g ( 孝) 2 仍g - + 仍g z2 【仍仍】t ；：， 2 2 6 ) 式中f 是刚量为一坐标，f ：2 z j ，而仍：旦碧，仍：罢，一1 f 1 将 ( 2 2 5 ) 和( 2 2 6 ) 式代入( 2 2 4 ) 中南人学硕士学位论文第2 章正演数值模拟方法的基本原理 后，具甲的早兀积分口j 写成： 工r ，u j q * d f = f ， 仍仍k 。d r 笼) = 嘭1 磅2 ) 芝) i lq j u * d f = f ， 仍仍k + d r 妻 = 巧1 谬 袭) 这里 乃o = f ，仍g + d f ，巧2 = f i 仍g d f ， 9 9 ) _ 工，伊。u * d f ，g 扩) _ 工，仍甜+ 订， 要写出对应于节点“i ”的离散方程式，需要把两个相邻单元1 1 和r 卜，的 贡献加起来，合成一项，以决定j4 1 丁- 4 - 点上“，和 q ，的系数： 万扩= 磁。+ 磅，g 玎= g 黔。+ g ：， 于是，式( 2 2 4 ) 对麻干点“j ”的装配方程可写成： n n c f “f + h 玎“= g 玎q 一b ( 2 2 7 ) i = lf = l 对于空间域q 的边界面元f _ ，假设三角形单元f _ ，的三个顶点坐标 ( x ，y ，z ) 为 | = 涨 “( f ) = 点嘶+ 岛“：+ 彘“，= 依心 3 g ( 善) = 氧吼+ 彘g ：+ 磊吼- - z 纯吼 k = l 式中仇= 氕。 工，吩g d r = 工，t 仍仍仍，g + d r t。u“，t=喜彬“。 c2 2 8 ) 1 2 中南人学硕十学位论文第2 章正演数值模拟方法的基本原理 f ，q ：u d f = f ，c 仍鲠仍，“d r 圣 = 喜爵g 。 c2 2 9 ) 式中：嘭= 工依g + d f ，酵= s i - 败“d r ( k = l ，2 ，3 ) ) m m c f “f + h 盯“，= g 玎g ，一b f ( 2 3 0 ) 对于方程( 2 2 7 ) 和( 2 3 0 ) 若记= ( f ) 和= + q ( f = ) ， 嘞甜，= g f ，9 j 哆 ( 2 3 1 ) “= g i f ，劬- b , ( 2 3 2 ) 【日】 u = 【g 】 q 一 b ( 2 3 3 ) 积分方程法是从场参数，所满足的微分方程边值问题出发，通过某些变换 导出有关参数所满足的积分方程，然后用近似计算方法求此积分方程的数值解， 并由此得出或进一步计算场参数的分布。以面积分法为例： 在电阻率为岛的半无限介质中，埋有电阻率为见的三维地质体，由场论中 静电场与稳定电流场的类比关系可知，为求解稳定的电流场，可利用电阻率分 界面上积累电荷的面分布。为此，设地下某一点m 的电位u ( m 1 为 u ( m ) = g o ( m ) + 眈( m ) + 虬( m ) ( 2 3 4 ) 式中砜( m ) 为供电电极a 在m 点产生的正常电位，即 玑f m l ：坦二 ( 2 3 5 ) ”、7 2 zr uf m ) 为地质体表面的积累电荷对m 点产生的异常电位： 中南大学硕士学位论文第2 章正演数值模拟方法的基本原理 v o ( m ) = l 掣d s 其中仃( q ) 为地下矿体表面s _ l ：4 ：r - - - 点q 的积累电荷面密度。 体表面i 积累电荷对m 点产生的异常电位，有 方。( m ) ：j 蛾i r 因而，地下任一点m 的点位表达式( 2 3 4 ) 可写成为： 吧) = 尝去j 警+ j j 掣d i ( 2 3 6 ) 眈( m ) 为镜像矿 根据边界条件求出积累电荷的面分布，便可由上式计算出电位值。 2 2 点源二维有限元法的基本原理 2 2 1 基本方程 ( 2 3 7 ) ( 2 3 8 ) ( 1 ) 位函数所满足的微分方程 在三维笛卡尔坐标系( x ，y ，z ) 中，假设电源供电点a ( x ay 。，z 。) 位于地表 时，电流强度为i ，产生稳定电流场的位函数u ( x ，y ，z ) 满足三维微分方程： v ( a v u ) = 一2 舾( x x a ) 6 ( j ，一y 一) 艿( z z ) ( 2 - 3 9 ) 对均匀介质，即电导率仃为常数。上式便为泊松方程： v 2 u = 一2 1 p 6 ( x x a ) 6 ( y y 一) 万( z z a ) ( 2 - 4 0 ) 若在无源空间，上式变为拉普拉斯方程： v 2 u = 0 ( 2 - 4 1 ) 对于三维地电条件，电导率仃= a ( x ，y ，z ) ，电位u = u ( x ，y ，z ) 。设地下介 质的导电性沿走向( y 轴方向) 无变化，即仃( x ，y ，z ) = a ( x ，z ) ；点电流源a ( i ) 位于坐标平面x o z 上，这时电场仍然是三维的，电位函数u = 饥五弘z ) 所满足的 微分方程变为： 珈琊，掣k 卜z ，学 + 昙l 盯( z ，z ) 掣l = 一2 硒( x x 一) 万( j ，) 万( z z 4 ) ( 2 4 2 ) 为了消除沿走向坐标y ( 使变为二维问题) ，故对上式两端分别沿y 轴方向 1 4 中南大学硕士学位论文 第2 章正演数值模拟方法的基本原理 作余弦付氏变换，上式变为： 昙( 盯娑) + 昙( 盯娑) 一a s o v ：f ( 2 - 4 3 ) ocoz仍 式中：盯= a ( x ，y ，z ) ；f = 一l ( 罗( x - - x 一) 万( z z 一) ；v ( x , 2 ，z ) 为电位u ( x ，y ，z ) 的付 氏电位，这样便将三维偏微分方程变成了二维偏微分方程。 ( 2 ) 位函数所满足的边界条件 由于电法勘探研究的稳定电流场分布于整个地下半空间，为了减少计算 量，在解正演问题时，通常把计算范围限定在一个有限的求解区内，这样便需 要在求解区的边界上对付氏电位v ( x ，a ，z ) 给定定解问题的边界条件，地表面r f 上( 如图2 - 1 ) ，电场强度的法向分量为零，故其相应的边界条件为： 地面 r x r | oo r i z 图2 - 1 所研究地电断面示意l 奎i 娑l ：0 ( 2 4 4 ) 锄lr 、 对于其它边界上，根据不同的假设，分别有三类边界条件： 第一类：= 0 ( 2 4 5 ) 第二类：_ a v l ：0 ( 2 - 4 6 ) o n r 第三类：( 娑+ 7 矿) l ：0 ( 2 - 4 7 ) o n r 1 其中r = 名k 。( 万) c o s ( 尹元) k o ( 办) ，k l ，k o 分别为一阶、零阶修正贝塞 尔函数。根据前人的研究，采用第一类边界条件时，计算的v 和u 值比理论值 偏小：采用第二类边界条件时，计算的v 和u 值比理论值偏大。于是周熙襄等 人采用第三类边界条件计算v 值，获得了比采用前两类边界条件更高的精度。 但当电极距较大时，计算的v 和u 值与理论值仍然有一定的偏差，故罗延钟教 授又对混合边佰条件式( 2 - 4 7 ) 中系数r l 的取佰进行了改进，其表示式变为： 中南人学硕士学位论文第2 章正演数值模拟方法的基本原理 刁= 兄 ，七k l ( a r k ) c o s ( v k 元) ，七k o ( 机) (2-48) 盯打 k=l|k = l 式中：，七一第k 个供电电极的电流强度； 五一第k 个供电电极的到边界点的矢径； 2 2 1 稳定电流场的边值问题 在二维地电条件下，由上述微分方程式( 2 4 1 ) 和边界条件式( 2 4 2 ) 所组成的，。 对若干个给定波数兄求解点电流源场的傅氏变换电位v ( x ，五，z ) 所满足的如下二 维偏微分方程的边值问题： 昙( 盯豢) + 昙p 娑) 卅仃y ：，q o xo xo zo z 娑：0 ，1 1 ，( 2 - 4 9 ) o n o - v + c o s ( r , n ) v ：o ，r 。 o nr 式中f = 一i 。6 ( x - x k ，z z k ) ，i k 为第k 个点电流源强度，f ，和f 。分别为 研究区域的地蓑边界和地下无穷远边