（地球探测与信息技术专业论文）椭球异常体中梯激电正反演研究.pdf

s c h o olo fe a r t hs cin c e g u iinu niv e r sit yo ft e c h n oio g y s e p t e m b e r ，2 0 0 8t oa p r ii 2 0 1 0 研究生学位论文独创性申明和版权使用授权说明 独创性声明 本人声明：所呈交的论文是我个人在阮百尧教授指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包 含其他人已经发表或撰写过的研究成果，也不包含为获得桂林理工大学或其它教 育机构的学位或证书而使用过的材料。对论文的完成提供过帮助的有关人员已在 论文中作了明确的说明并致以了谢意。 学位论文作者( 签字，描瓜垩 签字日期：丞! ! ：王：笸 关于论文使用授权的说明 本人完全了解桂林理工大学关于收集、保存、使用学位论文的规定，即：按 照学校要求提交学位论文的印刷本和电子版本；学校有权保留送交论文的印刷本 和电子版，并提供目录检索与阅览服务：学校可以采用影印、缩印或其他复制手 段保存论文；在不以赢利为目的的前提下，学校可以公开论文的部分或全部内容。 ( 保密论文在解密后应遵守此规定) 学位论文作者( 签字) ：避 指导教师签名：卫乞五至至一 签字日期2 0 l o 5 厶j 一 中文摘要 以及区域内 的大致位置 及类型，但是对异常体的埋深、规模和物性参数等信息仍难以确定。本文以中梯 激电的理论为基础，选择椭球体作为研究的地质异常体，采用解析和数值计算相 结合的方法对不同场源的椭球异常体进行了正演拟合反演，以确定异常体的各模 型参数。 在正演方面，针对野外地质体形态多样的特点，本文推导了均匀场中水平大 地任意形态椭球体激电异常表达式，计算并分析了典型椭球体的主剖面及旁侧剖 面的视电阻率和视极化率异常特征。在点源场正演计算时，采用四面体单元剖分 总电位有限单元法来计算起伏地形的激电异常，使用交错剖分技术使得系数矩阵 的非零元更少，提高了程序计算速度，研究了利用边界单元法计算地形校j 下因子， 提出了任意地形下利用数值计算中已有网格坐标来计算各个节点立体角的计算 方法，通过比值法对地形影响进行改正，以典型地形验证了所编写的地形改正程 序正确性。用f o r t r a n 9 0 开发工具编制了相应的计算程序，并用算例验证了程序 的可靠性和运行效率。 在反演部分，首先参考层状大地电阻率、极化率最优化反演的理论和方法， 依据均匀场中任意形态椭球体的激电异常解析计算式，建立了以椭球体的轴长、 埋深、倾斜角度、电阻率和极化率为模型参数的最优化反演目标函数，通过奇异 值分解法求得模型参数修正量。给出了合理的修正步长。然后分别把解析解和数 值解作为观测数据讨论了各种情况下的适用条件、反演效果及效率。 反演算例表明，本方法精度高、速度快、对初始模型有一定适应性、可用于 复杂地形中梯激电数据的反演分析和解释。 关键词：中间梯度法；激发极化法；椭球异常体；有限单元法；边界单元法；最 优化反演 桂林理工大学硕士学位论文 a b s t r a c t a b s tr a c t a l t h o u g ht h ei n f o r m a t i o no ft h ea b n o r m a lb o d y sl o c a t i o na n dt y p ec a nb e a c c e s s e dd i r e c t l yf r o mt h em a po fp o f i l ea n dc o n t o u r ，t h ed e p t h ，s i z ea n dp h y s i c a l p a r a m e t e r so ft h ea n o m a l o u sb o d ya r es t i l ld i 伍c u l tt od e t e r m i n e i nt h i sp a p e r , t h e e l l i p s o i db o d y i ss e l e c t e da sas t u d yt a r g e to fg e o l o g i c a la n o m a l yt oi d e n t i f yt h em o d e l p a r a m e t e r sb yf o r w a r df i t t i n gi n v e r s i o nc o m b i n i n ga n a l y t i c a la n dn u m e r i c a lc a l c u l a t i o n s f o rd i f f e r e n tf i e l ds o u r c ei nt e r m so ft h et h e o r yo fi p i nf o r w a r d ，i pa n o m a l yf o r m u l ai sd e r i v e do fe l l i p s o i da n o m a l o u sb o ( 1 yw i n l a r b i t r a r ys h a p ea c c o r d i n gt ot h ev a r i o u sf o r m so fg e o l o g i c a lb o d y t h ef e a m r e so ft h e m a i np r o f i l ea n dl a t e r a lp r o f i l e s sa p p a r e n tr e s i s t i v i t ya n di pa r ec a l c u l a t e da n d a n a l y z e da tt h es a m et i m e i np o i n ts o u r c ef i e l df o r w a r dc a l c u l a t i o n , i tm a k e st h e u n d u l a t i n gt e r r a i no ft h ea n o m a l o u si pd a t ae a s i e rt oo b t a i nb yt h et o t a lp o t e n t i a l s t a g g e r e dt e t r a h e d r a lm e s hf e m ，a n dr e q u i r e sl e s ss t o r a g eo fn o n - z e r oe l e m e n t ， i m p r o v e st h ec a l c u l a t i o ns p e e do fp r o g r a ma n dl a y st h es o l i df o u n d a t i o nf o rr a p i d i n v e r s i o na sar e s u l t f o ri pd a t ai nc o m p l e xt e r r a i no b t a i n e db ym e a n so fp o i n t s o u r c ef i e l dt h e o r ya n dn u m e r i c a lc a l c u l a t i o n ，a r es t u d i e db yu s i n gt h eb o u n d a r y e l e m e n tm e t h o dt oc a l c u l a t et o p o g r a p h i cc o r r e c t i o nf a c t o ri na n yg r o u n d 鲥dw i t h n o d ec o o r d i n a t e si nt h ec a l c u l a t i o no fs o l i da n g l e ，c o r r e c t e dt h ee f f e c t so ft e r r a i nb y t h er a t i om e t h o d t h ea c c u r a c yo ft h et e r r a i nc o r r e c t i o np r o g r a mi sd e m o n s t r a t e db y a n a l y z i n gt h et y p i c a lt o p o g r a p h y ，d i s c u s s i n gt h ei n v e r s i o ne f f e c to ft h em a i np r o f i l e a n dl a t e r a lp r o f i l e s t h ec o r r e s p o n d i n gc o m p u t e rp r o g r a md e v e l o p e db yt h et o o l so f f 9 0i no r d e rt ot e s ta n dv e r i f yt h er e l i a b i l i t ya n de f f i c i e n c yo ft h ep r o g r a mi nl i n ew i m e x a m p l e sr e s u l t s i ni n v e r s i o n t h ei n v e r s i o no b j e c t i v ef u n c t i o na r ee s t a b l i s h e di nw h i c ht h ea x i s d e p t h ，t i l ta n g l e ，r e s i s t i v i t ya n di n d u c e dp o l a r i z a t i o nr a t ea r eu s e d 器t h ei n v e r s i o n m o d e lp a r a m e t e r si nr e f e r e n c et ot h eo p t i m i z a t i o ni n v e r s i o nt h e o r ya n dm e t h o d so f 1 a y e r e de a r t hr e s i s t i v i t ya n dp o l a r i z a t i o nt h e c o r r e c t i o nv a l u eo fm o d e lp a r a m e t e r s a r eo b t a i n e db yu s i n gt h em e t h o do fs i n g u l a rv a l u ed e c o m p o s i t i o n t h ea m e n d m e n t s t e pi sg i v e nr e a s o n a b l e t h e n ，t h ep r o c e d u r e sa r ed i s c u s s e dt h ea p p l i c a b l ec o n d i t i o n s ， i n v e r s i o ne f f e c ta n de f f i c i e n c yt a k i n gt h ea n a l y t i c a la n dn u m e r i c a ls o l u t i o n s r e s p e c t i v e l ya st h eo b s e r v a t i o nd a t ai nav a r i e t yo fs i t u a t i o n s i n v e r s i o ne x a m p l e ss h o wt h a tt h em e t h o di sa c c u r a t e ，f a s t ，h a ss o m ef l e x i b i l i t yo n t h ei n i t i a lm o d e la n dc a nb eu s e df o ri pd a t aa n a l y s i sa n di n t e r p r e t a t i o ni nc o m p l e x t e r r a i n k e yw o r d s ：c e n t r a lg r a d i e n ta r r a y , i n d u c e dp o l a r i z a t i o n ，e l l i p s o i da b n o r m a lb o d y , f i n i t ee l e m e n tm e t h o d ，b o u n d a r ye l e m e n tm e t h o d ，o p t i m i z a t i o n a li n v e r s i o n i i 桂林理工大学硕士学位论文 目录 目录 摘要i a b s t r a c t i i 第1 章前言1 1 1 弓i 言1 1 2 直流电法中梯激电椭球异常体正反演的研究现状、实际应用效果以及存在的问题1 1 3 论文主要研究内容和成果3 第2 章任意形态椭球体中梯激电异常计算5 2 1 水平大地椭球体激电异常解析计算5 2 2 算例9 2 3 本章小结1 0 第3 章复杂地形三维电阻率、极化率快速正演研究l l 3 1 有限单元法的边值问题和变分问题l l 3 2 不同剖分形式下有限元系数矩阵非零元数量比较1 l 3 3 三维点源场有限单元法正演计算1 3 3 3 1 插值计算1 3 3 3 2 单元积分1 4 3 3 3 边界积分1 4 3 3 4 总体合成1 5 3 3 5 系数矩阵按m s r 非零元素存储代码1 5 3 3 6 方程组求解1 9 3 4 本章小结2 0 第4 章三维电场直流电阻率边界单元法地形改正2 1 4 1 边界元模拟法2 l 4 1 1 边值问题2l 4 1 2 边界积分方程2 l 4 1 3 边界元法。2 1 4 1 4 立体角的计算方法2 3 4 2 算例2 4 4 2 1 算例l 2 4 4 2 2 算例2 2 4 4 2 3 算例3 2 5 4 3 本章小结2 6 第5 章椭球异常体最优化反演方法2 7 5 1 椭球异常体最优化反演方法。2 7 5 2 椭球异常体最优化反演流程2 8 5 3 本章小结2 8 第6 章水平大地均匀场椭球异常体激电数据反演效果2 9 6 1 模型反演效果分析2 9 6 1 1 球体模型2 9 6 1 2 椭球体模型3 0 6 1 3 多椭球模型31 6 2 反演误差分析3 3 6 2 1 将旁侧剖面数据当作主剖面时的反演结果。3 3 桂林理工大学硕士学位论文 目录 6 2 2 误差数据的反演情况3 4 6 2 3 异常体电阻率极化率估计误差对反演结果的影响3 5 6 2 4 只用视极化率数据或视电阻率数据对反演结果的影响3 6 6 2 5 初始模型对反演结果的影响3 6 6 3 本章小结3 7 第7 章点源场椭球异常体中梯激电数据最优化反演分析3 8 7 1 起伏地形椭球异常体中梯激电数据反演。3 8 7 1 1 模型一3 8 7 1 2 模型二4 0 7 2 本章小结4 l 第8 章结论和建议4 3 8 1 主要研究成果4 3 8 2 今后需要继续研究的问题4 4 j s 【谢4 5 参考文献4 6 个人简历4 9 桂林理j 【大学硕士学位论文 第1 章前言 1 1 引言 第1 章前言 直流电法中的激发极化法【1 - 2 1 ( 简称激电法) 是以岩石、矿石激电效应的差 异为物理基础，通过观测和研究大地激电效应，来探查地下地质情况的一种物探 方法。中间梯度装置( 简称中梯法) 是将正负供电电极置于相距较远处，而观测 是在其中间的三分之一地段进行，在岩石为均匀、各向同性的情况下观测区可以 近似为均匀电场。中梯法异常直观，工作效率高，因此在野外数据采集中使用极 为广泛。 近年来，随着观测技术的不断发展，观测水平的不断提高，激电法的应用领 域不断拓宽，已在寻找金属、非金属矿产和地质构造中发挥着重要的作用，在水 文工程、地质工程、环境、考古、地质灾害评价等与国民经济建设密切相关的领 域也得到广泛的应用。与此同时直流电法的迅速发展也对相应的理论与技术提出 了更高的要求，因此发展更加可靠快捷的反演方法是众多的地球物理学者孜孜追 求的目标。 对于实测中梯激电数据，通常根据曲线异常特征就可以大致判断异常的类型 和位置，但是对于异常体的形态特征和几何特征等重要信息却难以确定。如何提 高对实测资料解释的精度和可靠性，这是新的反演方法有待解决的问题。 椭球体可以更好地近似成球体、板状体、柱状体等野外常见异常形体，因此 把椭球异常体作为研究目标，进行正反演研究具有理论和现实意义。 1 2 直流电法中梯激电椭球异常体正反演的研究现状、实际应用效果 以及存在的问题 椭球异常体中梯激电正演计算的方法有三种：解析法、模型实验法和数值模 拟法。解析法主要适用于水平大地均匀场计算，数值模拟法适用于复杂电性分布 和复杂边界状况的点源场的数值计算。在解析法正演计算方面，a h 3 a 6 0 p o b c k l 4 f i 和傅良魁的著作中研究了基于直流电法理论的一定形态椭球体中 间梯度电位和极化率计算式【l 】【3 1 。肖宏跃( 1 9 9 1 年) 计算了板状体及直立椭球体 的视极化率并给出了基于b a s i c 语言的计算程序h 1 。刘崧( 1 9 9 4 年) 、强建科( 2 0 0 7 桂林理工大学硕士学位论文 第1 章前言 年) 等研究了椭球体频域激电异常特征呻1 。但是这些椭球体的形状都有一定限 制，并不能真正的表示空间任意形状的椭球体。在数值计算方面，2 0 世纪六、 七十年代，数值积分法、有限差分法、边界单元法和有限单元法逐渐成为电法的 正演计算和解释的重要方法【7 1 引，( a l f a n o1 9 5 9 ，d i e t e r1 9 6 9 ， g e o s i e n c e1 9 6 5 ， p a t e r s o na n dg r a n t19 6 9 ，j h c o g g o n19 71 ，h o h m a n n19 7 5 ，s n y g e r19 7 6 ) ，首先 用数值方法解决电法地球物理问题。随后，边界单元法和有限单元法迅速发展 1 3 - 2 3 】( w l r o d i1 9 7 6 ，p r i d m o r e1 9 7 8 ，f o x ，h o h m a n na n dk i l l p a c k1 9 8 0 ，o k a b e 1 9 8 1 ，n a r d i n ia n db r e b b i a1 9 8 3 ，o p p l i g e r1 9 8 4 等) 。9 0 年代以来，我国学者徐世 浙、马钦忠、李金铭、张献民、阮百尧、吴小平等人在前人工作的基础上，发展 了适合复杂地形三维地质模型正演模拟的数值计算方法。2 0 0 6 年吕玉增【2 牝5 】等 人又将四面体交叉剖分技术、m s r 存贮技术和s s o r p c g 技术成功运用到地一井、 井一地直流电法的正演研究中。 2 0 世纪5 0 年代随着直流电法在矿产资源勘查等方面的大规模应用，地形影 响问题研究开始受到重视。因为起伏地形不但可以引起假异常，而且会掩盖地下 由矿体或目标物引起的真异常，如果不能正确的认识和消除，就会导致完全错误 的解释结果。对起伏地形资料解释，人们一般采用经验消除、带地形直接反演和 比值法消除地形，经验有时并不可靠，带地形直接反演不仅难度增加，而且可靠 性也降低。因而用比值法来消除地形影响更为常见，所谓比值法即是用实测视电 阻率除以“地形校正因子( 地形模型视电阻率地形模型电阻率) ，从而得到改正 后或水平地形下的视电阻率。比值法中对地形校正因子计算的方法很多。早期的 有基于物理实验的定性解释法，基于数学公式的解析计算方法，比如保角变换法。 但由于定性解释的局限性以及数学公式的复杂性和只能计算简单地形的缺点，这 些方法的应用受到一定的限制。近些年来，随着计算机技术的发展和物探解释方 法的进步，人们研究了适合于起伏地形的2 d 、3 d 数值模拟和解释技术，从而使 地形影响问题的模拟、识别和校正成为可能【2 】。这些方法主要有：面积分方程 法，边界元方法，有限差分法和有限元方法等。其中的边界单元法利用格林函数 将3 d 边值问题转换成二维计算，因而计算量大为减少、计算速度大为提高，更 有利于大型数值计算以及快速反演的实现，从而得到了更为广泛的应用。徐世浙 等【4 h 5 】使用立体角方法导出位方程，提出模拟3 d 地形影响的新方法，用边界元 方法来解积分方程。该法将计算区域剖分成多个三角小单元，这样所有的网格节 点均是孤立点，因而如何计算这些点的立体角是需要解决的问题。 直流电法的反演经历了一维、二维到三维的发展过程。1 9 7 8 年，p e l t o n 4 6 】 提出了一种确定地下模型网格电阻率、极化率的最d - - 乘最优化反演方法，利用 参考数据库中数据插值来计算偏导数，该方法不能用于解释复杂地电断面。 2 桂林理工大学硕士学位论文 第l 章前言 s a s a k i t 4 。7 】随后通过电位函数和模型参数间的关系来计算偏导数，减少了计算工作 量。之后，s a s a k i t 4 8 】将解的光滑性加入目标函数中，使得解的稳定性更高。阮百 尧掣4 9 】提出了一种针对层状大地的电阻率、极化率数据最优化反演方法，取得 了良好的效果，在生产实际中得到了广泛的应用。1 9 9 9 年，阮百尧【5 0 j 在s a s a k i 等人的基础上，给出了一种电阻率和极化率分块双线性变化的有限元最小二乘二 维反演方法。该方法在正演计算中，采用了三角单元剖分，使实测数据在反演前 不需要进行地形改正；并在目标函数中加入了最简单模型和背景场等先验信息， 反演结果更接近实际情况，并在找矿、找水和工程勘察等方面，取得了很好的应 用【5 1 , 5 2 】。近年来，国内外开展电法三维反演研究【5 3 5 5 1 ，有基于b o r n 近似的三维 反演、层析成像反演、t a r a n t o l a 反演以及最d - - 乘等反演。b o r n 近似反演忽略 积累电菏间的相互影响，因此适应性比较差，特别对电性变化较大的复杂电性结 构的反演的分辨率较低。电阻率层析成像是日本学者s h i m a 5 6 】( 1 9 8 7 ) 首先提出 的，此方法要求目标体的周围都要有观测数据点，才能保证反演过程的收敛，因 而实用性受到了很大的限制。t a r a n t o l a 反演是基于概率分布的反演方法，但由于 三维模型的数据量大，反演收敛到可靠的解比较困难。最小二乘法是地球物理反 演中的精确反演方法，但是反演过程不可避免要求取雅可比矩阵和求解大型方程 组，对计算机的运算速度和内存要求很高，三维反演的计算量大的惊人。近年来， 国内外不少研究者对上述方法进行了改进，如吴小平利用共扼梯度法进行电阻率 三维反演【57 1 ，毛先进用积分方程法进行电阻率三维反演1 5 8 1 ，这些方法尽管减少 了计算量，但反演效果难以令人满意。黄俊革等人用异常电位法求取偏导数矩阵， 在目标函数中，加入了体积因子等先验信息，提高了反演精度和解的稳定性【5 4 , 5 9 】， 但是仅使用于水平地形的情况。 三维反演由于方程欠定会引起的多解性问题，虽然加入了先验信息可以解决 这个问题，但是得到的解也只是在约束条件下的最优解，在某些时候可能未必是 真解。 1 3 论文主要研究内容和成果 综合上述问题，本文将椭球体的几何尺寸和物性参数作为反演目标函数中的 待定量，避免了三维反演因为模型参数量过多引起的多解性问题。研究基于椭球 异常体的最优化反演方法、适用条件和反演效果，以对异常体的大小、规模、埋 深等形态特征和物性参数进行确定。本文主要进行如下相应的研究和工作： 1 、椭球体正演解析计算 3 桂林理工大学硕士学位论文第1 章前言 在经典激电理论的基础上，通过三次坐标旋转，推导了水平大地均匀场中任 意形态椭球体中梯激电异常的解析计算式，讨论了不同形态椭球体的主剖面及旁 侧剖面异常特征，编写了解析计算程序。 2 、椭球体正演数值计算 在吕玉增、阮百尧等人编写的三维有限元程序基础上。用交错四面体单元来 剖分计算区域，使有限单元法生成的系数矩阵所需存贮单元更少，计算时间仅为 十几秒，在p c 机上实现了椭球体模型的快速正演计算，为接下来的快速反演解 释奠定基础。 3 、起伏地形影响改正 利用边界单元法消除数值计算中地形影响问题，对不满足线性关系的电源点 附近的节点电位关系做了修正，推导边界单元法地形改正公式中的立体角的计算 方法，并把该方法应用到地形影响改正中，编写相应的地改程序，用算例验证程 序的可靠性和效率问题，并对典型地形的地形影响问题进行了讨论。 4 、最优化反演的方法及效果 结合最优化反演的理论和方法，将椭球体的几何尺寸和物性参数作为反演目 标函数中的待定量，研究基于椭球异常体的最优化反演方法，编制相应的反演程 序。针对均匀场水平地形解析数据和起伏地形数值计算数据的反演进行了详尽的 讨论和分析，对反演的效果和效率及适应条件进行了讨论。 算例表明：本文提出的反演方法能够有效的对中梯激电椭球异常体的资料进 行反演解释，但反演对数据质量有较高要求，误差超过5 的数据会导致反演不 收敛或不可信。反演对初始模型有一定适应性。先验信息的加入对于避免反演陷 入局部收敛十分重要，特别是异常体的电阻率或极化率，甚至是反演成功的必要 条件。 4 桂林理工大学硕士学位论文 第2 章任意形态椭球体中梯激电异常计算 第2 章任意形态椭球体中梯激电异常计算 为使计算形体更符合野外实际异常体，本章在参考前人成果的基础上推导了 经三次坐标旋转后具有不同角度和三轴尺寸的多个椭球异常体的电位、视极化率 和视电阻率解析表达式，以卡尔松完全椭圆积分实现了解析式中标准勒让德椭圆 积分，用f o r t r a n 9 0 编制了相应程序，利用该程序计算了两类典型椭球体的主剖 面及旁侧剖面的视电阻率和视极化率异常。 2 1 水平大地椭球体激电异常解析计算 如图2 1 所示，水平大地均匀场中一水平椭球异常体，地电模型参数：背景 电阻率为, o l ，极化率为r h ；椭球异常体的电阻率为：岛，极化率分别为：7 7 2 ， 中心点坐标分别为：而，y o ，z o ，三轴尺寸分别为：a ，b ，c 。均匀外电场岛 沿x 轴方向。 一 图2 1 模型示意图图2 2 视电阻率微分计算示意图 图2 1 所示的模型当沿x 方向作梯度测量时，其一次场电位k 、总场电位y 、 视极化率瑶分别为【l 】： y x ：- j 。p - - - - - 坚x ( 1 2 ，r ) l r l l 。 ( 2 1 ) 儿竹2 ( 1 ) 呓昙( 斌) 咖 甫2 v x l x1 一。二l l 其中w 、屹、口、t o 是与模型电阻率、极化率和椭球体几何尺寸有关的量， 电异常计算 电阻率的微 ( 2 2 ) 图2 3 坐标变换示意图 若将图2 1 所示的水平椭球体分别在坐标面x o y 、x o z 和o z ”按逆时针逐 次旋转角度口、和y 后，则地面某测点p 在x y z m 中的坐标为： 工m = 一z s i n + ( x c o s 口一j ，s i n 口) c o s y 卅= ( 工s i n 口+ y c o s 口) c o s 7 一 z c o s + ( j c o s 口一) ，s i n 口) s i n s i n y z m = ( x s i n o t + y c 。s 口) s i n y + z c o s f l + ( x c o s 口一y s i n 口) s i n c 。s y ( 2 3 ) 相应地，一次场电位、总场电位分别为： k 。= - j o p li x - 2 x ”c o s a c o s p v ；。l x 。 - 2 y ”( s i n 口c o s 厂一c o s 口s i i l s i n y ) 吖。l y 。 一2 z ( s i n 口s i n 厂+ c o s 口s i n p c o s y ) 巧。r 。 6 桂林理工大学硕士学位论文 第2 章任意形态椭球体中梯激电异常计算 = - j o p l ( 1 - r l t ) e x - 2 z ”c o s 口c o s v 。r 。 - 2 y ”( s i n 口c o s 7 一c 0 s 口s i n s i i l 7 ) 1 ，y 。 ( 2 4 ) 一2 z ”( s i n 口s i n y + c o s 口s i n c o s ，) 矿。r 。 由叠加原理，若水平大地均匀场中存在行个任意形态椭球异常体，在忽略这刀 个异常体之间相互影响时，电位计算只需要将异常电位部分相加则可，则相应的 视电阻率和视极化率计算式为： 伸喜p 和2 甄 ( 1 _ 州一， 挈 “ 主小妒2 篁 , t f v vo ( 矿f f f ) = 尚卜窆i = 1b 妒魏y ，挈l l j ( 2 6 ) 其中， - = 以一c o s 口c o s 0 = 乃( s i n a c o s y c o s a s i n f l s i n y ) 乞= 五( s i n 口s i n 厂+ c o s 口s i n c o s y ) 其中，f 是x ”，y ，z 之一，是坐标旋转变换相应值。 4 7 r a b c s 2 舯： 一啜 ( 9 2 + 岛) ( g 2 + t o ) v s = s s j 其中，g 、g 是石m s 、s t 方向椭球体半轴长， 甲= 一焉+ 南+ 南 以下主要针对三类典型地质特征体分析式( 2 4 ) 和( 2 5 ) 中的r 表达式。 i ) 球体 = 竿t 量：f ( 2 7 ) 3 ( # + 订+ 彳2 ) 3 “ 7 桂林理：r = 大学硕士学位论文第2 章任意形态椭球体中梯激电异常计算 i i ) 等轴椭球体 有两轴相等，对于呸= 6 c ， 当a i q 时 l x , = 2 7 r a j b i c i = b i 当a i 岛 c j ，对于任意形态椭球体可以通 过坐标旋转来满足该关系式。将r 中标准勒让德椭圆积分改为卡尔松完全椭圆积 分可得【l 】， 三2川 ，f、 桂林理工大学硕士学位论文第2 章任意形态椭球体中梯激电异常计算 f = 丽4 7 r a 丽, b l c i ( b 一q 2 ) 七2s i n 3 3 ( 口f 2 6 f 2 ) 如一 痂神饥 其中，母为第一类、第二类卡尔松( c a r l s o n ，1 9 7 9 年) 完全椭圆积 分，函数的定义如下 2 2 算例 如( w ，z ) = 互3f 瓜而而 砰( 训，z ) = 圭 ( 2 1 2 ) 2 2 1 模型一 为了分析水平椭球体各剖面的异常变化情况，本例设计了如下的模型，如图 2 4 所示，水平大地均匀场中一低阻高极化水平椭球体参数为：a = 2 m ，b = l m ， c = 0 5 m ，中心埋深, o = 3 m ，局= l o o f 2 m ，局= 5 q m ，r l = 1 ，仍= 5 0 ，第 一到第四旁侧剖面距离主剖面分别为l m ，2 m ，3 m ，4 m 。应用式( 2 1 ) 和式( 2 2 ) 计算的中梯装置各剖面视极化率和视电阻率曲线如图2 4 所示。对于低阻高极化 椭球体，视极化率曲线表明主剖面异常最大，旁侧剖面离主剖面越远异常越小， 视电阻率曲线主剖面异常最低，旁侧剖面离主剖面越远异常越大。 图2 4 模型一及视极化率和视电阻率曲线图 9 桂林理1 = 大学硕士学位论文第2 章任意形态椭球体中梯激电异常计算 2 2 2 模型二 为了分析水平椭球体经旋转后异常变化情况，如图2 5 所示，本例设计了模 型参数与模型一相同的一水平椭球体，分别旋转0 0 ，3 0 0 ，6 0 0 ，9 0 0 后，水平大 地均匀场中各参数，应用式( 2 5 ) 和式( 2 6 ) 计算的中梯装置主剖面视极化率 和视电阻率曲线如图2 5 所示。对于低阻高极化椭球体，视极化率曲线表明不旋 转时的异常最小，随着旋转角度的增大椭球体最高点距离地面越近，异常随之逐 渐变大，而视电阻率曲线表明当椭球体不旋转时的异常最大，随着旋转角度的增 大异常逐渐减小。 r l , j ，_ 一 l 。 - - 0 - - o 一兰车一 廿一 z 三趱i；?：蓬蓦4 | - - 昂= ：甏鹱 、蕊矿 。- ：l _ _ “* - 工j m _ 2 3 本章小结 图2 5 模型二及视极化率和视电阻率曲线图 本章在激电理论的基础上，经三次坐标旋转后推导了具有不同角度和三轴尺 寸的多个椭球异常体的地面电位和视极化率解析表达式，并通过等效电阻率法的 微分计算式推导了视电阻率计算公式，在此基础上编写了相应的程序，算例给出 了典型水平椭球体及旋转椭球体的异常特征，为下一步中梯激电数据的快速反演 工作奠定了基础。 1 0 桂林理工大学硕士学位论文第3 章复杂地形三维电阻率、极化率快速正演研究 第3 章复杂地形三维电阻率、极化率快速正演研究 本章采用有限单元法来实现对电阻率、极化率的数值模拟。首先比较了不同 剖分方式下有限单元法生成的系数矩阵中非零元数量的多少，从而选择一种最节 省内存资源的剖分方法，以达到提高程序计算效率的目的。然后从电场满足的边 值问题出发，给出了电场满足的变分方程，用m s r ( m o d i f i e ds p a r s er o w ) 存 贮系数矩阵中的非零元，最后用s s o r p c g ( s y m m e t r i cs u c c e s s i v eo v e rr e l a x a t i o n p r e c o n d i t i o n e dc o n j u g a t eg r a d i e n t ) 法来求解方程。 3 1 有限单元法的边值问题和变分问题 双点电源中电位满足的边值问题可归结如下【删： v ( o v u ) = - ( 4 z c o a ) 1 6 ( a ) + ( 4 7 r c o s ) 1 6 ( b ) q o u o n = 0 r ， ( 3 1 ) 蛳+ 兰r 2 - r 1 融元矿一r 2 ，- ) r 。l 吒 眨 j 其中，在彳、b 分别供单位正、负点电源，c 为区域q 的地面边界，l 为 区域q 的地下边界，n 为边界的外法向方向，仃为介质的电导率，u 为总电位， ( 9 a 、分别是么、b 点对地下区域q 张的立体角。 与( 3 1 ) 式等价的变分问题为： f ( 甜) = 1 2 盯( v “) 2 - ( 4 ，r c o a ) 1 6 ( a ) “+ ( 4 z r c ) 1 6 ( b ) “】d q + 1 2f 盟i 垒c o s ( 相一- - rc o s ( r 一2 ，元) k ( 3 2 ) ( 吃一r ll 眨 j d f ( u ) = 0 3 2 不同剖分形式下有限元系数矩阵非零元数量比较 m s r ，即改进的行压缩存储，对称矩阵时仅存贮行下标小于等于列下标的 桂林理- t 大学硕士学位论文第3 章复杂地形三维电阻率、极化率快速正演研究 非零元。下面以图3 1 为例，依次按z 、x 、y 方向，2 3 x 4 的网格来分析不同 剖分形式下有限单元法形成的系数矩阵的非零元数量。 以网格中的第1 7 号点为例，如果采用六面体来剖分计算区域，则比1 7 号点 小又与该点共单元的点有1 ，2 ，3 ，4 ，5 ，6 ，7 ，8 ，9 ，1 3 ，1 4 ，1 5 和1 6 ，共 计1 3 个点。如果采用四面体单元且按照图3 2 ( a ) 或3 2 ( b ) 的剖分方式来剖 分计算区域，则比1 7 号点小又与该点共单元的点有2 ，4 ，5 ，6 ，8 ， 1 3 ，1 4 ， 1 5 和1 6 ，共计9 个点；如果采用四面体单元且按照图3 2 ( a ) 和3 2 ( b ) 的剖 分方式交替剖分计算区域，则比1 7 号点小又与该点共单元的点有5 ，1 4 和1 6 ， 共计3 个点。 从上面的分析可见，对于m s r 存贮方法，用六面体剖分计算区域时需要存 贮的非零元最多，四面体不交错剖分需要存贮的非零元居其次，而四面体交错剖 分需要存贮的非零元最少。以9 0 6 0 4 5 的网格为例，六面体剖分需要存贮 3 4 6 3 1 4 1 个非零元，四面体非交错剖分需要存贮2 4 6 6 8 4 1 个非零元，四面体交错 剖分需要存贮1 7 4 9 9 9 1 个非零元。由此可见选择四面体交错剖分需要的计算时间 会更少。 图3 1 规则网2 3 x 4 网格剖分示意图 1 2 正演研究 体单元e ( 3 3 ) n i = a i x - f 岛少+ c f z + d ，= 等 ( 3 4 ) 其中，v 是四面体单元体积，k 是插值点p ( x ，y ，z ) 与四面体其他3 个顶点 ( _ ，= 1 ，2 ，3 ，4 ，i ) 所组成的四面体体积，a f 、岛、c f ( f = 1 , 2 ，3 ，4 ) 是与四面体单 元顶点坐标有关的常数。以：，。( k ) 为例，各系数的表达式为： 同理，可得： ， 1 y = 一 6 五月 艺耽 x 3y 3 x 4y 4 i 奶 q = i 乃 l 儿 z l 1 乞1 z 3 1 z 4 1 z 2 1 z 3 1 z 4 1 ，：，。：k ：= 1。 o ，岛= 一 x y z1 x 2y 2z 2 1 x ，y z j 1 x 4y 4z 4 1 斗= 巨 y 2 1 y 3 1 y 4 1 快速正演研究 图3 3 四面体单元 3 3 2 单元积分 式( 3 2 ) 第一项单元积分 珐盯( v u 心q 2 睦砥尝) 2 + ( 詈) 2 + ( 尝) 2 d x d y d z e 。e 。 i ， 。 = 要u j k l e