（地球探测与信息技术专业论文）中梯装置双频激电法电阻率数据地形校正研究.pdf

中南大学硕士学位论文 摘要 摘要 双频激电法属于频率域激电法，是一种非常有效的地球物理勘探 方法，它具有抗干扰能力强、仪器轻便、稳定、观测精度高等优点。 该方法能够在地形复杂、环境恶劣的自然条件下开展工作，并在我国 的矿产资源勘查中得到了越来越广泛的应用。 长期以来，起伏地形下的地电异常是电法勘探的一个难题，因为 这个问题不仅本身复杂，而且目前的科技发展水平也达不到完全解决 的地步。近年来经过地球物理专家们的努力，起伏地形下的二维电阻 率正演及反演成像问题已经基本解决。本文主要针地形起伏地形下， 双频激电中梯装置电阻率数据的地形改正进行研究，希望通过对双频 激电采集到的电阻率数据进行地形改正，使改正后的电阻率数据与幅 频率f 。一样，也成为判断地下异常的依据。 本文研究了双频激电法电阻率数据点电源边界单元解法以及二、 三维起伏地形的地形校正，采用边界单元法在二、三维地形条件下对 边界进行剖分，并在前人的工作基础上进行了更加完善的公式推导。 编制了二维地形条件下视电阻率地形校正的程序，建立了纯地形模型 并对各种不同二维均匀起伏地形进行了研究，总结了规律，建立了非 均匀条件下二维起伏地形模型，对低阻直立板状体与高阻水平板状体 模型在各种不同地形下的影响进行了分析，验证了本文算法的正确性 以及程序的可靠性，为二、三维地形条件下的反演成像研究奠定了基 础。 关键字：双频激电法，边界单元法，地形校正，异常电位 中南大学硕士学位论文 摘要 a b s t r a c t d u a l 仔e q u e n c y i n m i c e d p o l a r i z a t i o n m e t h o d b e l o n g sf 诧q u e n c y d o m a i ni pm e t h o d ，a n d “i sav e qe f i e c t i v eg e o p h y s i c a lm e t h o d ，“h a s m o s ta d v a n t a g e s ，s u c ha sa n t i - i n t e r f e r e n c ea b i l i t ) ，l i 曲t ，s t a b l ea n d h i g l l p r e c i s i o no b s e r v a t i o n ，s oi tc a nb eu s e di nt h ec o n d i t i o no fc o m p l e x t e n ! a i na n db a dn a c u r a lc o n d i t i o n a n dt h i sm e t l l o dh a sb e e nu s e dm o r e a n dm o r ew i d e l yi nm i n e r a le ) 【p l o r a t i o n t h er e s e a r c ho fg e o e l e c t r i ca n o m a l yo nt h eu n d u l a t et o p o 铲印h yh a s b e e nad i m c u l tp r o b l e mi nt h ef i e l do fe l e c t r i c a lp r o s p e c t i n gf o ral o n g t i m e ，i ti sn o to n l yc o m p l i c a t e di t s e l b u ta l s ot h el e v e lo ft e c h n o l o g y d e v e l o p m e n tc a n n o tb ea c h i e v e d i nr e c e n ty e a r s ，t b r o u 曲t h ee 筋r t so f g e o p h y s i c a le x p e r t s ， t h e p r o b l e m o f2 - do nu n d u l a t et o p o g r a p h y r e s i s t i v 时f o n v a r da n di n v e r s ei m a g i n gh a sb e e ns o l v e d i nt h j sp a p e r ，w e m a i n l vr e s e a r c ht h et e 啪i nc o l l r e c t i o no fr e s i s t i v i t yd a t aa b o u td u a l 舭q u e n c yi n d u c e dp o l a r i z a t i o nm e t h o d o nu n d u l a t et o p o g r a p h y ，a n d m a d ei tb e c o m e sa n o t h e re v i d e n c et oi u d g eu n d e r 星，o u n dm e d i a t h i sp a p e rr e s e a r c h e dt h eb o u n d a 珂e l e m e n tm e t h o df 研r e s i s t i v i 谚 m o d e l i n gw i map o i n t c u r r e n ts o u r c eo fd u a l f r e q u e n c y i n d u c e d p o l 砌i o n m e t h o da n dt e r r a mc o r r e c t i o ni 芏l2 d 3 du n d u l a t et o p o 孕a p h y ， w eu s e db e mt o p a r t i t e t l l e b o u n d 撒y ， a n dc a r 巧 o u tam o r e c o m p r e h e n s i v ef o 彻u l ad e r i v a t i o no nt h eb a s i so fp r e v i o u s ，as u i to f2 一d f - o 刑a r da l g o r i t h j nc o d e si sa l s oa c c o m p l i s h e d ，w ee s t a b l i s ha p u r em o d e l a n dr e s e a r c hd i f f e r e n t2 du n d u l a t et o p o 铲a p h ym o d e l s ，a n de s t a b i i s h e d n o n u n i f o 咖2 du n d u l a t i n gt e r r a i nm o d e l s ，觚a l y s i sl o wv e n i c a lt a b u l a r b o d ya n dh i g hl e v e lp l a t eb o d ym o d e l s ，i tv e r i f i e dt l l e c o r r e c to fm e a l g o r i t h l l l ，a sw e na sm er e l i a b i l i 锣o ft l l ep r o c e d u r e t h e s er e s e a r c h a c h i e v e m e n t so fd cr e s i s t i v i t vh a se s t a b l i s h e df o 明d a t i o nf o r 向r t h e r s t u d y i n 9 2 d0 r3 dr e s i s t i v i 够f o n a r dm o d e l i n ga n di n v e r s i o ni m a g i n go n u n d u l a t et o p o g r a p h y 中南大学硕士学位论文 摘要 k e y w o r d sd u a l 能q u e n c ym d u c e d p o l a r i z a t i o nm e t h o d ，b o u n d 哪 e l e m e n tm e m o d ，t e m i nc o r r e c t i o n ，a n o m a l o u s p o t e n t i a l 中南火学硕士学位论文第一章绪论 第一章绪论 1 1双频激电法电阻率数据地形改正研究的意义 双频激电法是由中南大学何继善院士提出的一种频率域激电观测方法，它是 第一个由中国人自己提出方法原理，又由中国人自己发明仪器的电法勘探方法。 双频激电法在我国除台湾以外的地区得到了推广和应用，找到了大量的铜、铅、 锌等有色金属矿产，为我国的有色资源金属勘探做出了很大的贡献u 。3 1 ，目前中 南大学基于双频激电法的基本原理，开发了一系列双频激电观测系统，并且在该 基础上研制出了伪随机多功能电法仪。 仪器采集地下数据以后，我们还要对数据进行各种处理，才能准确地得到探 测区地下的地质与地球物理特征。双频激电法与其它电法勘探方法相比具有许多 本质优点，但目前在数值模拟和野外资料处理方面还存在很多不足，如对许多不 规则地形还无法得到它的解析解或数值解，资料处理仅仅局限于一些定性解释方 法和一些经验方法上，使得解释结果与实际情况差别较大h 1 。而且目前在双频激 电法找矿中，地形起伏一般较大，以至于激电中梯装置采集到的电阻率数据基本 不可用，幅频率f 。只能作为定性或半定量解释，因此必须对测得的电阻率数据 进行二维地形改正，使采集的电阻率数据也能成为判断地下介质属性的依据。电 阻率数据的地形校正问题是目前双频激电法进一步发展的瓶颈，深入研究该问题 对于双频激电法的理论发展和增强实用很有必要。 起伏地形对二、三维电阻率数据影响很大，导致真实异常“畸变”、“偏移” 甚至是“面目全非”，给实测资料的解释带来了相当大的困难。只有从理论上正 确认识地形对电阻率数据的影响，才能在野外资料的处理中给出正确的解释。数 值模拟手段是进行地球物理方法研究的重要武器，用其研究起伏地形条件下电阻 率异常规律有重要意义。 本论文的目的就是通过对双频激电法所采集到的电阻率数据进行二维地形 改正研究，使电阻率数据也能成为判断地下介质属性的重要依据，从而使双频激 电法在今后的野外勘探中能充分地发挥其优势。如今我国许多金属矿山已知储量 已经非常有限，继续寻找新的接替资源，尤其是我国西南、西北等地区，地形条 件极其复杂，地形起伏一般较大，通过我们对双频激电所采集到的电阻率数据进 行地形改正研究，使电阻率数据成为可用数据，并与幅频率数据相参照来做后期 的数据处理，将对异常的解释将会更加可靠。因此我认为深入研究该问题将对双 频激电法的发展起到积极的推动作用，该课题具有极大的研究价值。 l 中南大学硕士学位论文第一章绪论 1 2 研究现状与水平 1 2 1 电阻率法数值模拟研究现状 地形影响问题是2 0 世纪5 0 年代伴随着直流电法在矿产资源勘查等方面的大 规模应用开始受到重视。目前主流的数值方法有：积分方程法、有限差分法、有 限单元法和边界单元法四种，其中有限差分法和有限单元法属于“网格法”范畴， 积分方程法和边界单元法属于“非网格法”范畴，这四种方法各有优缺点。 ( 1 ) 积分方程法( i n t e g r a le q u a t i o n ) 。积分方程法是从场参数( 场强或场 位) 所满足的微分方程边值问题出发，通过某些变换导出有关参数( 场强、位、积 累电荷密度或散射电流密度) 所满足的积分方程。然后用近似计算方法求此积分 方程的数值解。积分方程法是一种经典的无网格数值模拟方法，其基础是关于物 理场的表现定理m 1 。在异常体( 不均匀或散射体) 的尺度不大时，与网格法相 比积分方程法具有物理意义明确、计算量小和计算速度快的优点，因此特别适用 于解决无界半空间中有限地质体对入射场的散射问题。 l a l f o n o ，k v o z o f f 哺1 及g v k e l l e r 和f c f r i s c h k n e c h t 等阻1 ，早在5 0 年代 末到6 0 年代初就提出了用积分方程法计算任意形状构造上地电体的视电阻率异 常；1 9 6 9 年k d i e t e r 等刚提出用积分方程法计算二维体的激电和视电阻率异常曲 线的具体方法；1 9 7 3 年a d e y 等1 和1 9 7 5 年g m h o h m a n n n 2 1 从麦克斯韦方程出发， 建立了更普遍的积分方程，把积分方程的应用领域推广到电磁法中；1 9 7 6 年 j j l a j o i e 等n 引，1 9 8 l 、1 9 8 2 年m o k a b e n 盯和u c d a s n 副等进一步发展了积分方程 法。 在国内电磁法领域最早应用积分方程法解决电法数值模拟问题的学者有周 熙襄n 朝n 、罗延钟n 踟等，随着积分方程的发展，出现了边界积分方程，在求解过 程中借用有限单元法的离散技术，进而发展成边界单元法。 ( 2 ) 有限差分法( f i n i t ed i f f e r e n c em e t h o d ，或简称f d m ) 。有限差分法是 一种经典的数值模拟计算方法，其基本原理就是用差商代替微商，将待求解的连 续微分方程变换为离散的差分方程，并通过求解差分方程得到原微分方程的近似 解。 有限差分法在地球物理领域有着广泛的应用，具有简单和直接的优点，而且 便于编程实现。最重要的是，有限差分法对任何复杂介质模型都适用。当然，有 限差分法也有很多不足，如果处理不恰当，可能会导致错误的计算结果。现有研 究表明，在有些条件下有限差分格式除了累计误差外还会引起数值频散、数值各 向异性、网格散射( 绕射) 以及会导致非指数不稳定性的漂移解。另外在计算方 2 中南大学硕士学位论文 第一章绪论 案中加入自由边界条件时，会由于数值微分算子在自由边界上的非对称性和低阶 性而导致整个计算过程的不稳定n 引。 6 0 年代末有限差分法开始应用于直流电法领域。j e p s e n ( 1 9 6 9 ) 讨论了二维地 电条件下点源电阻率法和激发极化法问题；l a m o n t a g n e 幢叫等( 1 9 7 1 ) 对矩形薄板的 电磁法响应进行了数值模拟研究；s w i f t 乜等对二维非均匀介质的大地电磁理论 响应进行了数值模拟；j e p s o n 乜2 1 ( 1 9 7 5 ) 利用有限差分法对电阻率法正演问题进行 数值计算；m u f t i ( 1 9 7 6 ) 改进了有限差分网格剖分方法，采用非均匀网格计算 任意形状二维地电体响应；7 0 年代末期，d e y 和m o r r i s o n 幢盯( 1 9 7 9 ) 将混合边界条 件引入有限差分法中，对二维任意形状地质体的电阻率进行正演模拟，提高了计 算精度，至今混合边界条件处理方案还在广泛使用；s c r i b e 盼射( 1 9 8 1 ) 对三维断面 的电场响应进行了计算；z h d a n o w 心引等( 1 9 8 2 ) 对二维地质体的电磁响应进行数值 模拟；g 1 d m a n 等( 1 9 8 3 ) 对三维断面的瞬变电磁响应进行数值模拟； l e e p i n 旧1 ( 1 9 9 2 ) 用有限差分法解决了2 5 维矩形电流回线的时间域电磁响应问题； s p i t z e r 睁驯等在有限差分的离散化和收敛速度方面进行了研究，对三维正演模 拟的精度和速度进行了分析。 我国的地球物理作者从8 0 年代开始研究有限差分法正演问题，如周熙襄n 等，罗延钟等2 1 为有限差分法的发展做出了相当大的贡献：刘树才b 3 1 等、刘正栋 b 引等在差分格式、网格剖分方法等方面做了很多工作。刘树才等对电法二维正演 中的若干问题作了详细论述，他指出有限差分法和有限单元法是解决二维问题的 理想方法：吴小平等汹基于有限差分原理利用共轭梯度法求解了点电源二维正演 模拟问题。 有限单元法和边界单元法是最常使用的两种地形改正方法，它们各有优缺点。 有限单元法汹1 是近二十年来随有限差分法之后迅速发展起来的又一种有效的数 值模拟方法，它属于区域型的数值计算方法，具有自动满足内部边界条件的功能， 无需对内部边界单独处理，所以在解决物性复杂分布的地球物理问题时有显著的 优点，但当涉及无限区域或三维问题时，有限单元法所需的节点大大增加，需要 很大的计算机内存，计算速度非常慢。 ( 3 ) 有限单元法( f i n i t ee l e m e n tm e t h o d 或简称f e m ) 。有限单元法是一种 求解场的变分问题的数值方法。用这种方法求解稳定电流场电位时，首先要建立 给定边界条件下电位的微分方程，并通过变分原理，等价地变换成变分方程，然 后将求解区域离散化，即按一定的规则将求解区域剖分成相互连接的网格单元， 导出单元节点电位的高阶线性方程，集成整个区域单元的节点方程，通过求解方 程组得到所有节点的电位。与有限差分法相比，有限单元法不仅能求出网格节点 的未知函数值，而且还可以求解出计算区域内任意点上的未知函数值。 3 中南大学硕士学位论文第一章绪论 有限单元法在电阻率地形改正上有很大的发展。c o g g o n ( 1 9 7 1 ) 针对直立导电 极化体，用有限元法计算得到平底、斜边的山谷和山脊2 - d 地形下电阻率及极化 率响应曲线；c o g g o n ( 1 9 7 3 ) 用有限元技术模拟比较极化率双偶极和三级装置的 地形影响；r i c h a r dcf o x 等( 1 9 8 0 ) 利用2 一d 有限元算法做了源于地形的偶极一偶 极视电阻率异常和地形( 山谷、山脊) 对极化率影响的系统研究。周熙襄等( 1 9 8 6 ) 应用有限元法计算偶极装置下台阶、山脊、山谷2 一d 地形对视电阻率异常曲线的 影响及在地表下存在低阻体时山脊、山谷2 一d 地形的组合视电阻率异常曲线。 p a n a g i o t i s 、t s o u r l o s 等( 1 9 9 9 ) 用2 5 维有限元模拟方法考察几种常用的电阻率 装置的地形影响并采用比值法进行校正，杨华和李金铭( 1 9 9 9 ) m 1 以山脊、山谷 地形为例，通过点源2 一d 有限元数值模拟给出了起伏地形条件下地形对充电法电 场的影响，并用比值法对地形影响进行校正，ht r u i n a nh 0 1 c o m b e 和g e o r g er j i r a c e k ( 1 9 8 4 ) 给出3 一d 有限元算法，对该算法进行多种测试并展示实际模型精确 结果，采用比值法校正。刘海飞( 2 0 0 8 ) 研究了起伏地形电导率连续变化的三 维激电数据有限元数值模拟，并取得了较好效果。 ( 4 ) 边界单元法( b o u n d a r ye l e m e n tm e t h o d ，或简称为b e m ) 。边界单元法 的前身就是边界积分方程法，它属于边界型的数值方法啪1 。随着边界单元法的兴 起，其单元划分和插值函数的概念引入到了边界积分方程法，发展成为边界单元 法。边界单元法最早在英国的s o u t h a m p t o n 大学土木工程系开始使用，并逐渐应 用到各个领域。边界单元法的优点是只在研究区域的边界上剖分单元，这样可以 降低求解问题的维数，使最后形成的代数方程组的元数减少，特别适合地球物理 中遇到的无限区域问题和三维问题；边界单元法的缺点是当区域内存在复杂的介 质分布时，用边界单元法求解各个介质的边界积分方程是非常麻烦的。 边界单元法在地形改正上的发展：在8 0 年代，f o x m l ，1 9 8 0 ：o k a b e h ，1 9 8 1 ： n a r d i n ia n db r e b b i a 池1 ，1 9 8 3 ：o p p l i g e r h 驯，1 9 8 4 ，为边界单元法的发展做出了 卓越的贡献。黄兰珍等( 1 9 8 6 ) ，介绍点源场电阻率法二维地形改正的边界元法， 并计算了山脊、山谷、角域地形上联合剖面法视电阻率异常曲线，展示用比值法 对实测曲线进行地改的应用效果。田宪谟等( 1 9 8 6 ) 介绍点源电阻率法3 维地形改 正的边界元法。给出用比值法电阻率三维地形校正的例子。徐世浙1 等利用边 界单元法成功地解决了二维、三维地电断而的正演模拟问题。徐世浙等：马钦忠 等，在资料解释方面也取得了一些成果。钱家栋和马钦忠( 1 9 9 2 ) 探讨了用边界元 法解决层状介质中存在二维横向不均匀情况下点源场视电阻率计算问题。黄兰珍 和田宪谟h 耵( 1 9 9 7 ) 介绍边界单元法应用于工程地质勘查电阻率法的地形改正实 例。刘继东( 1 9 9 8 ) 介绍求解异常电位的二维边界单元法对电测深资料进行比值法 地形改正的应用效果。 4 中南大学硕士学位论文第一章绪论 总体来说，以上四种数值模拟方法在电阻率法正演计算方面都有一定的优势 和不足，边界单元法的优势是正演速度快，内存需求少，主要用于地形改正和地 下少量地质体的正演模拟。有限差分法的优点是原理简单，编程容易实现，计算 速度相对比较快，其缺点是当物性参数复杂分布或场域的几何特征不规则时，适 应性比较差。有限单元法在处理复杂几何形状时，其灵活性和适用性比其他方法 要好，缺点是计算量比较大，计算时间长。 1 2 2 地形改正国内外研究水平与现状 自2 0 世纪7 0 年代以来，复杂地表条件下的地球物理场数值模拟一直受到研究 者的重视。迄今为止，取得了一些重要的成果，但是并未彻底得到解决。就复杂 地表条件下的地球物理场数值模拟来讲，目前的状态是发展与问题并存，机遇与 挑战并存h 引。 根据定义，勘探地球物理的目的是通过对观测数据的数学物理解释和地质分 析来得到关于地下地质构造特征以及岩石物理性质方面的精确信息，并以此为基 础来确定埋藏在地下的矿产资源分布以及储量情况，从原理上来说，无论是数学 物理解释还是地质解释，都是一种从已知到未知的过程，具有反问题的性质和特 点。因此，对各种条件下地球物理场的时空分布特征和强度的研究就构成了对地 球物理观测资料进行数学物理和地质解释的基础矧。 在直流电法领域，国内的研究者早在2 0 世纪8 0 年代就研究了2 5 维地电模型 中的起伏地形问题，先后利用过有限单元法和边界单元法。从理论上将，2 5 维 复杂地表条件下稳定电流场的数值模拟问题已经基本解决，剩下的问题是如何使 得理论问题实用化。但对于三维问题研究较少，这与国际上的研究状态基本类似 【柏】 o 处理起伏地形问题常用的方法有网格映射法一1 、插值法协1 、曲网格法娜刀 等，解决起伏地表能力也各有千秋。网格映射法要求计算地形函数相关偏导数， 对相对光滑或曲率连续的起伏地形有很好的适应性；插值法对起伏地形进行一定 的近似，可以适用于较陡峭的地形特点，但在处理过程中边界条件的确定要划分 多种类型，应用起来比较复杂；曲网格法在物理意义上较为明确，比较贴近地形 本身的特点，但同样处理起来比较复杂。 f o x 旧1 ( 1 9 8 0 ) 阐述了地形对电阻率和激发极化的影响，c h o u t e a u m l ( 1 9 8 8 ) 研究了二维大地电磁法的地形改正问题；徐世浙h 等人实现了用边界单元法进 行二、三维地形改正的问题，因为边界单元法可以降低空间维数，且具有数据量 少、计算速度快、精度高等优点，所以该法拟合地形可以达到较高精度，是二、 5 中南大学硕士学位论文第一章绪论 三维地形改正中一种有效地方法；l u n t a ot o n ga n dc h i e h h o uy a n g m 3 用有限 元作正演计算，用二阶m a r q u a r d t 最小二乘法实现了带入不平地形的反演计算， 最终实现二维起伏地形反演的算法，合成数据和实测数据显示效果较好；阮百尧 等n 3 1 改进了日本应用广泛的电阻率激发极化法二维有限元反演程序，可起伏地形 反演电阻率激发极化法模型剖面；吴小平h 钉于2 0 0 5 年实现了非平坦地形条件下 电阻率三维反演：吕玉增等 5 1 也通过四面体剖分实现了复杂地形条件下电阻率三 维有限元的数值模拟；刘海飞 引，强健科等人实现了起伏地形下电导率连续变 化的三维激电数据有限元数值模拟。 1 9 9 2 年，t e s s m e r 等嘲1 提出了一种可以被称为整体间接法或是曲化平方法的 处理二维起伏地表的方法。在这种方法中，起伏地表首先通过一定的映射规则转 化为水平地表，然后再对经过曲化平映射后的波动方程运用伪谱法和形式时间积 分法进行数值分析；1 9 9 4 年t e s s m e r 和k o s l o f f 旧1 又进一步将这种方法的基本思想 用于解决三维起伏地表问题；此后，h e s t h o l m 等旧确1 对这种曲化平方法进行了改 造并将其与有限差分法相结合，研究了起伏地形条件下粘弹性介质模型中地震波 数数值模拟问题。在对自由边界条件的处理上，t e s s m e r 等人采用的是特征变量 法，而h e s t h 0 1 m 及其合作者采用从局部坐标到整体坐标的坐标变换法。 1 3 研究内容和创新点 1 3 1 研究内容 本文首先从边界单元法的基本原理出发，基于格林公式、6 函数和基本解的 问题，研究了双频激电法电阻率数据点电源的边界单元解法。通过确定点电源场 基本方程和边界条件，采用边界单元法进行带地形的二维电阻率法数值模拟。 采用编制出的二维起伏地形下视电阻率地形改正的程序，对一些典型地电模 型进行模拟，通过对比地改前后的视电阻率曲线，验证了程序的效果，说明该套 程序的算法是正确的，并通过纯地形的异常模拟，总结了山脊、山谷、陡坎以及 各种组合地形条件下的视电阻率影响曲线以及规律。 研究了复杂介质起伏地形二维电阻率数值模拟。采用“比值法对地下介质 中低阻直立板状体与高阻水平板状体在各种典型地形条件下的地形改正，分析了 改正前与改正后的电阻率曲线，并验证了程序处理非均匀半空间起伏地形的地形 改正也是行之有效的，为今后双频激电法电阻率数据三维复杂条件下的应用提供 了依据。 6 中南大学硕士学位论文 第一章绪论 1 3 2 创新点 ( 1 ) 首次提出基于双频激电法电阻率数据的地形改正研究，利用边界单元法， 实现了起伏地形条件下二维电阻率法正演问题，可以对复杂地电模型进行模拟。 ( 2 ) 该正演程序基于c + + 开发环境编译，在程序编译过程中，采用插值法对 数据进行插值拟合，使正演所得到的曲线更加准确。 7 中南大学硕士学位论文第二章边界单元法的基本原理 第二章边界单元法的基本原理 边界单元法( b o u n d a r ye l e m e n tm e t h o d ) 是新发展起来的一种数值计算方 法，它比“区域”型方法，如有限单元法，有重大的优越性。有限单元法在区域 内划分单元，单元和节点的数量多，因而方程组的规模浩大，边界单元法在边界 上划分单元，单元和节点的数量少得多，因而方程组的规模也小得多。有限单元 法不适于无界区域，在遇到无界区域时，有限单元法要用一个足够大的有界区 域代替无界区域边界单元法适于无界区域。地球物理勘探中的大部分数值计算是 在无界区域内进行，因此边界单元法特别合适。 2 1 边界单元法的数学基础 2 1 1 格林公式 格林公式是边界单元法用于地球物理勘探的最基本的公式啪1 ，在将边值问题 转化为积分方程的过程中，格林( g r e e n ) 公式是非常有用的工具。设区域q 中 有两个函数u 和巾，格林公式为： i 胛2 伊一妒2 p q = 彳( 警一缈筹卜 c 2 ) 其中是f 区域q 的边界，对于三维区域，r 是边界面，对于二维区域，r 是 边界线，刀是r 的外法向。格林公式简证如下： 根据场论中哈密顿算子v 的运算规则有： v ( v 9 ) = v 伊2 + v v 伊 ( 2 2 ) 根据区域积分与边界积分的关系有： ，。( 胛缈) 船彳驷 汜3 ) 对于v ( 咖) 也是同样结果，两者相减，则可以得到格林公式( 2 1 ) 。 8 中南大学硕士学位论文第二章边界单元法的基本原理 2 1 26 函数 狄拉克6 函数在边界单元法中是个十分方便的工具，以p 点为中心的6 函 数定义如下： 万( p ) = o ，当q 中的点不在p 时 万( p ) = ，当q 中的点在时 p ( p 矽q = 1 ，当p 在q 内部时 n 对于三维6 函数，上式中的q 是三维区域，而对于二维6 函数，q 是二维区 域。 根据上述定义，有以下几个重要的性质： ( 1 ) 当p 点在q 的边界r 上，且p 处的边界光滑时，有： l 万( p ) d q = 丢 ( 2 ) 设函数甜在p 点连续，当p 点在q 内部时， 弘万( p ) d q = ( p ) o 其中( p ) 是点p 的函数值。 ( 3 ) 当p 点在q 的边界上时，且p 处边界光滑，有 2 1 3 基本解 弘艿( p 矽q = 寻( p ) q _ 微分算子l 对函数进行某种微分运算，构成微分方程： 三( ) = o ( 2 4 ) ( 2 5 ) ( 2 6 ) ( 2 7 ) 若某函数缈经过l 的微分运算，得到一个一万( p ) 函数，则称妒为三( ) = o 的 基本解。 我们分别给出二维和三维拉普拉斯方程的基本解： 9 中南大学硕士学位论文 第二章边界单元法的基本原理 ( 1 ) 二维拉普拉斯方程，v 2 = 0 ，其中v 是二维哈密顿算子，它的基本接是 缈：士l n ! ( 2 8 ) 2 万， 、。7 其中，是二维区域q 中某点p 至q 中任意点的距离。 ( 2 ) 二维拉普拉斯方程，v 2 = 0 ，其中v 是三维哈密顿算子， l 矽= 一 4 万， 其中，是三维区域q 和某点p 至q 中任意点的距离。 它的基本解是 ( 2 9 ) 中南大学硕士学位论文第三章起伏地形下二维边界元数值模拟 第三章起伏地形下二维边界元数值模拟 双频激电法( d f i p ) 是由中国人自己提出方法原理，又由中国人自己发明 仪器的电法勘探方法，它是一种频率域激电观测方法。双频激电法在我国除台湾 以外的地区得到了推广和应用，找到了大量的铜、铅、锌等有色金属矿产，为我 国的有色资源金属勘探做出了很大的贡献。 应用双频激电法在寻找地下矿产资源过程中，中梯装置观测方式囚其效率高、 对异常反映明显，在实际工作中常常被采用以确定异常体的平面分布。然而由于 受起伏地形的影响，采集的原始视电阻率数据不能直接使用，通常仅利用单一的 幅频率数据进行解释。鉴于此，为更充分的发挥双频激电中梯装置在找矿中的应 用效果，我们针对双频激电中梯装置采集的视电阻率数据二维地形改正做了一些 研究，试图去消除地形起伏的影响，双频激电电阻率数据与一般的直流电法电阻 率数据处理是类似的。 3 1边值问题 为了研究问题方便，来建立一个二维地电断面模型。设地下介质均匀，电阻 率为p = l 。选取z 轴平行地面，x 轴垂直z 轴并保持水平，j ，轴垂直向上。选取 地表的4 点为坐标原点，在彳点供入i = 1 的点电源，则电位u 的基本方程是： v 3 u 三一2 万( 彳) ( 3 1 ) 其中万( 么) 是以彳点为中心的三维万函数。 以彳为中心，以厂= 为半径，平行地形走向在地下画假想的半圆柱面r 。， 与地面r 形成闭合曲面，则有下面电位的边界条件成立： 丝i ：o ( 3 2 ) 圳l 其中力为外法线向量，半圆柱面r 。上的电位是点电源的电位， 甜l ：，兰 ( 3 3 ) 甜i r j2 丽 其中c 是比例系数。 通过傅氏变换，将三维微分方程转化为二维微分方程，因此我们定义下列傅 1 1 中南大学硕士学位论文 第三章起伏地形下二维边界元数值模拟 氏变换m 1 u ( 五弘七) = j c o 甜( x ，只七) c o s 舷击 其中u = 研u 】，因此有 f c 窘，= 害 f r 争害 f 埘u 另外有 j c o 删) c o s 乜出= 掣 其中万( 彳) 是以a 为中心的二维狄拉克函数代入( 3 1 ) 有， v 2 u 一后2 u = 一及a ) 其中v 是二维哈密顿算子，若对足够多的七求解方程( 3 7 ) ， 变换， 材( x ，少，z ) = 昙f u ( x ，y ，七) c 。s 拓出 即可就得电位函数甜。 为求解( 3 7 ) ，我们还需要确定u 的边界条件： 在地表，法向，2 与z 无关，因而有 u i l2 蜊= c j c o 游 = ( 扫 ( 3 4 ) ( 3 5 ) ( 3 6 ) ( 3 7 ) 再进行反付氏 ( 3 8 ) ( 3 9 ) 根据边界条件( 3 2 ) 有 丝i - o 锄j r ( 3 1 0 ) 此处l 为地形剖面线。 对r 。上的电位z f 进行反付氏变换，根据( 3 3 ) ，有 中南大学硕士学位论文 第三章起伏地形下二维边界元数值模拟 u l = j c o 如删= c j c o 寿 = 呱( b ) 其中k 是第二类零阶修正贝塞尔函数，r 。是砂平面上的半圆周。 我们采用二维格林公式 ( 3 1 1 ) ( 叩2 矿一妒2 d 伫唾( 喀一伊警) 订 ( 3 1 2 ) 其中q 是砂平面上的二维区域，如图2 1 ，r 是q 的边界线，r = r ，+ r 。， 刀是r 的外法向，令 口：墨盟 ( 3 1 3 ) 2 万 图孓1 二维地形和坐标 其中，是q 中某点尸到q 中任意一点的距离。矿有如下特性： v 2 伊一七2 妒= 一伊( 力 其中万( p ) 是以尸点为中心的迪拉克函数。 将( 3 7 ) ，( 3 1 4 ) 代入( 3 1 2 ) 的左侧，得 ( u v 2 缈一妒2 u ) d q = 一u 万( p ) d q + 妒( 么) d q 根据艿函数的积分性质，有 ( 3 1 4 ) ( 3 1 5 ) 中南大学硕士学位论文 第三章起伏地形下二维边界元数值模拟 筘( 么) 纶竿= 其中r 是彳到尸点的距离。 1 k ( 欣) 22 万 在边界单元法解题过程中，若边界不光滑，则尸点的积分为 州彬q = 乏) ( 3 1 6 ) ( 3 1 7 ) 其中u ( p ) 是尸点的u ，哆是尸点对区域q 的张角。若p 点边界光滑，则有 国f a 。 将( 3 1 6 ) ( 3 1 7 ) 代入( 3 1 5 ) ，得 ( 孵9 一妒2 州似象) + 警 c 嗉一伊筹灯= l 曝r 掣坍 2 l 喙c 乌笋，亳订= - - u 警c 。s 胁 其中k 是第二类一阶修正贝塞尔函数。 对于边界r 。来说，由于厂t - o o ，所以，于是积分 c 喀一妒等班= c u 警一伊券沙 兰丝! 堑2 苎! 堑2 + 兰丝! 堑2 堡! 堑2 2 万2 7 = o ( 3 1 8 ) ( 3 1 9 ) ( 3 2 0 ) 将( 3 1 8 ) ，( 3 1 9 ) ，( 3 2 0 ) 代入( 3 1 2 ) ，得 鲁m ，= 警+ l u 掣c o s 触 组2 。) 当地形水平时，有胪：9 0 0 ，c o s 卸，q 2 石，上式被积函数为o ，所以 砜( p ) = 掣堕 ( 3 2 2 ) 这是地形水平时u ( p ) 的分布规律。 将( 3 2 1 ) 写成下式， 生u ( p ) ：( p ) + 【u 丝丛堕c 。s 肛玎 1 4 中南大学硕士学位论文 第三章起伏地形下二维边界元数值模拟 通过求解3 2 3 式的积分方程则可得到地表各点的电位u 。 3 2 边界的网格剖分 如图3 2 所示，用刀个节点将r ，剖分成疗一1 个单元，每个单元足够小，近似 为直线。电源点在某单元l 的内部。将边界积分分解为各个单元积分之和，由 ( 3 2 3 ) 式可得， 丢u = 扛+ 薯u 警c o s 脚 慨2 4 ， 其中是地形水平时的电位。 图3 2 ：地表的边界剖分 单元两端节点的编号为_ ，七，其坐标为( ，乃) ，( 以，儿) ，为了进行单元积分， 故作如下的坐标变换： 孝，= 型丘业 x i y t x k y i 考t = 1 一考； 毒和善，是x ，y 的线性函数，而且满足： t ，点：= l ，彘= o k 点： 孝，= 0 ，磊= l 除了电源所在单元l 外，假定u 在各单元上是线性变化的，用q ，表示 _ ，j j 的u ，则单元上的u 可表示为 u = 缶，( ，= ，七) 中南大学硕士学位论文 第三章起伏地形下二维边界元数值模拟 单元r ，上的积分为 l u 挚c o s 朋 2 戳专掣c o s 肛 = 厶u ( 3 2 5 ) 其中 石= 旬掣c o s r ，( h 问 ( 3 2 6 ) 上式积分计算解析式是很困难的，我们求助于高斯求积公式计算。 3 3 积分方程求解 在电源所在单元l 上u 为线性变化的假设和实际情况有重大出入，因此l 的积 分按照如下方法计算比较好。 单元l 上u 近似为 u ：鱼盟+ c 2 石 ( 3 2 7 ) 其中c 是地形影响的附加项，为待定常数。由于_ ，点的r = 和u = u ，得 到4 段的c ， c ：u 。一螋 92文 脚o k r危 k 土幼 + q l l 为 u u的段 4 以 所 中南大学硕士学位论文 第三章起伏地形下二维边界元数值模拟 e u 半c o s 胛 = 鹏+ 型喾型】 堕盟c 。s 脚 = u 丘掣c o s 触+ l 型喾型半c o s 肛 公式( 3 3 0 ) 为4 段的积分。同理，4 段积分为 i l l u 半c o s 砸 = 以l 挈c o s 腑+ l 型雩型掣c o s 卿 将4 和4 段积分相加，得 p 半吣8 弧= f p l 七眠如 其中 石= f ，u 半c o s 胁 厶= lu 半c o s 触 e = l 型喾型掣c o s 胁 + l 型雩型半c o s 肌 各单元积分之和： 喜lu 挚c o s 触= 【f 。弓吒】 u 珥虬 r + e = f u 斗c j ( 3 3 0 ) ( 3 3 1 ) ( 3 3 2 ) ( 3 3 3 ) ( 3 3 4 ) 其中巧是歹节点两侧单元的乃之和；只= 【f - 弓瓦】；u = 配q ) 。 1 7 中南大学硕士学位论文第三章起伏地形下二维边界元数值模拟 对节点f ，( 3 2 7 ) 式可以写成 詈q = + 善l u 挚c o s 触电岬+ c ( 3 3 5 ) 对每个节点都得到如上式的一个方程，由全部节点得到方程组，记作 其中 ( 三国一f ) u ：+ c 万 f q i 缈2 i 。 io l 巧= 【互巧e 】t = ，) r c = c 。c ，c 。) t ( 3 3 6 ) ( 3 3 7 ) 方程( 3 3 7 ) 含有行个方程和刀个未知数u 虬。用高斯消去法解此线性代 数方程组，可得到各个节点的u ，进而得到各点地形改正以后的视电阻率数据。 我们在双频激电法中梯装置电阻率数据地形改正工作中发现当剖面上地形 起伏较剧烈时使用边界元法得到的数值解，在地形变化最大的山脊部位会出现不 稳定，有时还会出现负电位。有专家们认为，这一现象的出现并非方法本身的缺 陷，而是地形起伏剧烈，节点数较大，线性代数方程组的阶数增加，系数矩阵条 件变坏，使一般的解方程法不能取得稳定合理的解所造成。 中南大学硕士学位论文第四章起伏地形下三维边界元数值模拟 第四章起伏地形下三维边界元数值模拟 起伏地形三维边界元数值模拟是解决实际问题中地形校正问题的关键，所以 我们所讨论的三维问题。三维问题可以通过傅里叶变换降为二维问题来处理，较 之各种域法，边界元法具有解决问题所需的方程组的阶数低、数据量少、计算速 度快、精度高等优点m 1 。实践证明，边界元法是求解电源三维地电模型位场问题 的一种实际可行的新的数值计算方法。 用电阻率法在起伏不平三维地形地区进行工作，消除地形对观测位场的畸 变影响，提高测量资料的释效果，是国内外地球物理工作者所关注的问题。 4 1边界单元法地形改正原理 假定地形线下的介质均匀，其电阻率p = 1 ， 电源。若彳处地表光滑，电位甜的基本方程为 v 2 z f = _ 2 万( 彳) 其中万( 彳) 是以a 为中心的狄拉克函数。 电位的边界条件是： l ，电流沿着地表1 ，流过，则有 丝i ：o 锄i r 成立，其中n 是r s 的外法向； 2 ，r 。上的电位是点电源的电位 l c “l l 2 7 在地表放置电流强度为l 的点 ( 4 1 ) ( 4 2 ) ( 4 3 ) 其中c 为比例系数。 确定了边界条件和电位的初值以后，我们再次利用格林公式： l ( 伊一妒2 “) 施= 4 鼍一伊豢灯 ( 4 4 ) 格林公式可以将区域问题转化为边值问题，其中q 为地下整个区域，r 为 q 的边界面，r = r 。+ r 。，刀是r 的外法向。 由前面可以知道，三维拉普拉斯方程v 3 = 0 的基本解为 中南大学硕士学位论文第四章起伏地形下三维边界元数值模拟 口：二 ( 4 5 ) 口= 一 l 4 aj 4 ，r ， 其中r 为q + r 中某点j p 至任意点的距离。 且不难证明， v 2 缈= 一万( p ) ( 4 6 ) 成立 将( 4 1 ) ，( 4 5 ) ( 4 6 ) 代入( 4 4 ) 式的左边，当点尸在边界上时， 得到： ( 彬缈一妒2 州q = 一扣( p ) d q + 挚