（地球探测与信息技术专业论文）25维电阻率断面的最小构造反演.pdf

硕t ：学位论文摘要 摘要 电法勘探是传统的地球物理勘探手段之一，已广泛应用于地质勘察、地下水 开采、金属矿产勘探等领域，兼顾电剖面法横向分辨率高和电测深法纵向分辨率 高的优点，电阻率层析成像成为继地震波层析成像和电磁波层析成像以后地学c t 的重要手段。但实际资料处理中，如水文、考古学、建筑等领域，由于解释速度 及有效性的要求，2 维、2 5 维电阻率剖面数值计算和反演解释仍然是我们的主 要手段，这有赖于对反演方法进一步研究。本文共分为四章，分别对电阻率断面 正演的基本原理和有限元数值模拟方法，以及最小构造反演目标函数和偏导数矩 阵求取等多个方面进行讨论。 本文首先在点源二维的正演模拟一章中，从地电断面电位满足的微分方程出 发，导出二维地电断面点源场总电位的边值问题、变分问题，程序实现中的各种 细节说明，包括矩阵的存储、大型稀疏方程组的求解方法选取、电阻率断面的正 演加速方法等。通过对地下一些简单地电模型的模拟试算，得出各种装置下的计 算断面图及相关误差分析，正演模拟的可靠性；为后期多方法反演奠定了基础。 其次在电阻率反演的一章中，一方面阐述了传统有限单元电阻率反演方法中的偏 导数求取，包括有限元标准方法、差分法、补偿电路法、格林函数法和摄动法等； 反演中的目标函数的各种可选形式，最小二乘、广义逆法，及阻尼因子和粗糙度 矩阵的引入；非线性问题下最优化过程中灵敏度矩阵的推导及深入分析；反演形 成的大型病态线性方程组的求解，共轭梯度法和奇异值分解算法；理论数据和野 外实践数据的反演效果图。另一方面，提出了与有限单元法相结合的偏导数的解 析法求取，即用可以求得解析解的球体单元来替代有限元网格剖分单元来求 取偏导数矩阵。最后在结论中得出，球体替代求取偏导数的方法基本达到预期的 效果。将地电断面的全域反演改为目标区反演，能够大大提高反演的效率，最小 构造二维反演迭代中基本不会引入数据分辨不出的多余构造干扰，利用共轭梯度 的求取法方程计算速度快、精度高、内存需求少，可适应于实际复杂结构的反演。 关键词数值模拟，有限元法，2 5 维，电阻率反演 硕士学位论文 摘要 a b s t r a ct e l e c t r i c a lp r o s p e c t i n gi so n eo ft r a d i t i o n a lg e o p h y s i c a le x p l o r a t i o n m e t h o d s i th a sb e e nw i d e l ya p p l i e di nt h ef i e l d so fg e o l o g i c a l e x a m i n a t i o n ，g r o u n dw a t e rd e v e l o p m e n t ，m i n i n ge x p l o r a t i o n ，e ta 1 f r o m 19 8 0 s h a v i n gh i g h e rh o r i z o n t a l r e s o l u t i o no fe l e c t r i c a lp r o f i l i n ga n d h i g h e rv e r t i c a lr e s o l u t i o no fe l e c t r i c a ls o u n d i n g r e s i s t i v i t yt o m o g r a p h y h a sb e c o m ea ni m p o r t a n ta p p r o a c ho fg e o s c i e n c ec tt h a th a v eh a d s e i s m i cw a v et o m o g r a p h ya n de l e c t r o m a g n e t i cw a v et o m o g r a p h y i nt h e d i s p o s a l o f p r a c t i c a l m a t e r i a l ss u c ha s h y d r o l o g y ，a r c h a e o l o g y ， a r c h i t e c t o n i c s ，e ta l ，w ea r es t i l lc a r r y i n go u t2 da n d2 5 dn u m e r i c a l c a l c u l a t i o na n di n v e r s i o nf o rt h es p e e da n dv a l i d i t y s ow em u s tp e n e t r a t e i n t ot h ei n v e r s i o nr e s e a r c h t h ed i s s e r t a t i o nw a sd i v i d e di n t o5c h a p t e r s t o t a l l y , i n c l u d i n gt h es e c t o r so ft h eb a s i cp r i n c i p l eo fr e s i s t i v i t ym o d e l i n g f i i l i t ee l e m e n tm e t h o d t h eo b j e c tf u n c t i o no ft h el e a s ts t r u c t u r er e s i s t i v i t y i n v e r s i o na n dt h ec a l c u l a t i o no fp a r t i a ld e r i v a t i v e sm a t r i x w h i c ha r e d i s c u s s e dr e s p e c t i v e l y f i r s t l y ，i nt h ec h a p t e ro f2 5d i m e n s i o n a lr e s i s t i v i t ym o d e l i n g ，t h e b o u n d a r yv a l u e ，v a r i a t i o ne q u a t i o no ft o t a lp o t e n t i a lo np o i n ts o u r c e e l e c t r i c a lf i e l do ft w od i m e n s i o n a lg e o e l e c t r i cs e c t i o na n ds o m ep r o b l e m o fp r o g r a m ( i n c l u d i n gt h es t o r a g eo fm a t r i x ，t h es o l u t i o nm e t h o do fb i g s p a r s ee q u a t i o na n dt h es p e e d u pm e t h o do fm o d e l i n g ) ，w e r ei n t r o d u c e d b e g i n n i n gw i 也t h ed i f f e r e n t i a le q u a t i o no fp o t e n t i a l w ec a l c u l a t e ds o m e s i m p l em o d e lb yu s i n gf e mm e t h o d ，g e tt h ep r o f i l ef i go fs o m ek i n d c o n f i g u r a t i o na n di m p r o v et h er e l i a b i l i t yo ft h ep r o g r a m s e c o n d l y i nt h e c h a p t e ro ft h ep r o c e s s i n go fr e s i s t i v i t yi n v e r s i o n ，o n eh a n d ，t h et r a d i t i o n a l c a l c u l a t i o no fp a r t i a ld e r i v a t i v e sm a t r i x ( i n c l u d i n gt h ef e ms t a n d a r d m e t h o d ，d i f f e r e n c em e t h o d , c o m p e n s a t i n gc i r c u i tm e t h o da n dg r e e n f u n c t i o nm e t h o di n2 da n d3 di n v e r s i o n ) ，t h eo b j e c tf u n c t i o no f r e s i s t i v i t yi n v e r s i o n ( i n c l u d i n gt h el e a s ts t r u c t u r ei n v e r s i o n ，t h eb r o a d s e n s ec o n v e r s em e t h o d ，a n dt h ei n t r o d u c eo fd a m pa n dr o u g h n e s sf a c t o r ) ， t h ec a l c u l a t ：i o no fp a r t i a ld e r i v a t i v e sm a t r i xi nt h eo p t i m i z a t i o n p r o g r e s s i n go fn o n l i n e a rp r o b l e m , t h es o l u t i o no fb i gm o r b i d i t ye q u a t i o n ( i n c l u d i n g t h e c o n j u g a t eg r a d m e t h o da n d s i n g u l a r i t y v a l u e d e c o m p o s i t i o n ) w h i c ha r ed i s c u s s e dr e s p e c t i v e l y a i m i n ga tt h er e s t r i c t i o n o fi n v e r s i o nb a s e do nf e m g e t t i n gt h el o g a r i t h mo fa c t u a lr e s i s t i v i t ya n d t h e o r e t i c r e s i s t i v i t y , c a l li m p r o v et h es t a b i l i t yo fi n v e r s i o n u s i n g e x p o n e n t i a lf o r mt oc o r r e c tv a l u eo f t h em o d e lc a nf e t c hu p 也eb u gw h i c h 2 e tm i n u sr e s i s t i v i t ya n dr e d u c ei t e r a t i o n si ni n v e r s i o n h a n g i n gf u l l i n v e r s i o ni n t ot a r g e ti n v e r s i o nc a l lr e d u c et h et i m eo fc a l c u l a t i o n s oi n t h i sp a p e r w ea d v i s et oc a l c u l a t et h ep a r t i a ld e r i v a t i v e sb yu s i n gt h e s p h e r i f o r me l e m e n t ，w h i c hh a st h e o r e t i c a lp a r t i a ld e r i v a t i v e s ，t or e p l a c e 硕上学位论文摘要 t h eg r i de l e m e n t t h ew a yu s i n gi nt h i sp a p e ri s g e tf r o mc o n j u g a t e g r a d i e n tm e t h o d ，w h i c hs p e e da n dp r e c i s i o ni sa p p r o p r i a t e ，f i t st o c o m p l e xs t r u c t u r a li n v e r s i o n k e yw o r d s n u m e r i c a lm o d e l i n g ，f i n i t ee l e m e n tm e t h o d ，2 5 一d r e s i s t i v i t yi n v e r s i o n ， i 原创性声明 本人声明，所呈交的学位论文是本人在导师的指导下进行的研究 工作及取得的研究成果。尽我所知，除了论文中特别加以标注和致谢 的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不 包含已获得中南大学或其他单位的学位或证书所使用过的材料。与我 共同工作的同志对本研究所作的贡献均已在论文中作了明确的说明。 作者签名：么罐日期：丛年上月上日 关于学位论文使用授权说明 本人了解中南大学有关保留、使用学位论文的规定，即：学校 有权保留学位论文，允许学位论文被查阅和借阅；学校可以公布学位 论文的全部内容，可以采用复印、缩印或其它手段保存学位论文；学 校可根据国家或湖南省有关部门规定送交学位论文。 作者签名：微导师签名：辽弛日期：越年上月上日 颈i ：学位论文第一章绪论 第一章绪论 1 1 前言 电法勘探是研究地球内部电性结构和电场分布的一种探测方法，在金属矿、 煤f f l 、石油天然气、地下水等探测中得到广泛且有效的应用，近年来在浅层工程、 环境、考古等与国民经济建设、人类社会生活密切相关的探测领域也得到很大的 拓展。在最近几年，尤其以高密度电阻率法的迅猛发展，已经被广泛的应用于隧 道洞身病害探测、建筑物基础检测、堤坝管涌探测、高速公路岩溶勘查、煤阳采 空区及地下古墓探测等方面，并且取得了良好的应用效果。随着计算机的迅速发 展及用数值方法模拟地电响应取得实质性突破，更加证实了该方法的可行性和有 效性，它的应用f j i 景也将因此变得更加光明。 反问题是地球物理中最为核心的问题和地球物理工作者的最终目的。其原理 就是根据地面上的观测信号推测地球内部与信号有关部位的物理状态。因此，反 问题的求解方法成为地球物理反演的主要研究对象。地球物理反问题主要是根据 最小二乘原理，利用正演模型和实测数据构造目标函数，并使其达到极小，然后 通过反复叠代，修改初始模型，最终得到满足误差条件的解。 电阻率法反演发展到今天，主要的手段是将非线性地球物理问题化为线性问 题。使非线性反演简化为线性反演问题，主要采用置换参数法、泰勒级数展开法 等；地球物理的非线性问题线性化后，使用最优化方法求解最常用的求解方法有 最d - - 乘法、广义逆法等。在地球物理反演中，传统最小二乘迭代方法得到非常 广泛、有效应用川( j u p pa n dv o z o f f , 1 9 7 6 ：m a c k i ee ta l ，1 9 8 8 ；s i cp a r k , 1 9 9 1 ； o l a y i n k a , 2 0 0 0 ) 。由于反演参数太多，传统的阻尼最小二乘反演往往导致过于复 杂的模型，即产生所谓多余构造，它是数据本身所不要求的或是不可分辨的构造 信息，给解释带来困难，在最小二乘准则中加入光滑约束，形成所谓的最小构造 反演方法拉1 ( s t e v e na n dc o n s t a b l e ，1 9 8 7 ：杨长福等，2 0 0 2 ) 反演求得光滑模型， 提高了解的稳定性。层析成像技术在地球物理反演中也得到广泛的关注口1 ( d a i l y ， 1 9 9 1 ；s a s a k i ，1 9 9 2 ：白登海，1 9 9 5 ) ，电阻率层析成像的主要原理是应用反投 影技术由测得的测区周围在各个不同方向产生电位差重构探测区内电阻率分布， 因此多应用于跨孔电阻率层析成像；近年来，共扼梯度( c g ) 方法是解大型最优化 问题最有效方法之一【4 】【5 】( 吴小平，2 0 0 0 5s a s a k i ，2 0 0 1 ) ，它为反演中需要存储大 型稀疏矩阵和耗用大量内存提出了很好的解决办法，直接从反演目标函数出发， 无须直接解正则化方程，可避免偏导数矩阵的计算和存储，提高计算速度、节省 机器内存。直流电阻率3 d 非线性反演方法，首先利用有限元法计算出正演结果， 硕i ：学位论文 第一章绪论 反演问题采用了快速的预条件共扼梯度算法，大大节约了计算时间和存储内存等 问题；其它诸如高斯一牛顿法( t a m b u f r i n 0 2 0 0 0 、张美玲等2 0 0 2 ) ，采用一种新的 雅克比算法建立了变分问题的离散牛顿方程以克服反演算子的高病态性问题以 及广义逆法( 陈乐寿等，1 9 8 9 ：杨长福等，2 0 0 2 ) 等在电祖率反演中都有相关的研 究。另外，用单一的勘探装置准确估算地下构造参数一直是非常困难的事，由此 联合反演的提出有效的克服了单一手段难以解决的多解性问题( m a t i a sd el a v e g a ，2 0 0 3 ) ，利用w e n n g , l - , 祠i 单极一单极装置联合反演区域内汽油污染的土壤电 阻率问题取得了很好的效果。尽管这些方法取得一定的成效，但局部线性化的迭 代方法不可避免的存在多解性和易陷入局部极小的问题，加之偏导数矩阵( j a c o b i 矩阵) 的求取及其存储的困难，使得电阻率二维反演难于较好实现。 而在进行地球物理反问题的过程中，最为关键的问题就是偏导数矩阵 ( j a c o b i n am a t r i x ) 的求取，大部分的反演时间都是耗在此，一些算法力图避开 灵敏度矩阵在多次的反演迭代进行计算，只是在反演的丌始才计算一次偏导数矩 阵，这就违背了我们反演的初衷；或是以某些可选择的系数或因子代之，或者对 野外数据提出特殊要求等，这些就会给反演工作增加许多变数，其反演结果将因 人而异。因此本文想在偏导数矩阵的求取方法上做一定的研究正演过程仍然 采用通用的有限元数值模拟方法，反演的时则不按照传统的“视电阻率对模型参 数的偏导数”的求取方法，而是尝试用可以求得解析解的球体单元来替代有 限元网格剖分单元，通过分析这两种方法的误差估计，评估该方法的有效性和同 用性。 1 2 本文研究的主要内容和成果 1 、有限单元法的正演模拟过程，程序实现中的各种细节说明，包括矩阵的 存储、大型稀疏方程组的求解方法选取、电阻率断面的正演加速方法等。各种装 置下的计算断面图及相关误差分析，为后期多方法反演奠定了基础。 2 、传统有限单元法中的电阻率反演方法，即：视电阻率对模型参数的偏导 数求取；介绍反演中的目标函数的各种可选形式，着重介绍最小二乘、广义逆法， 及阻尼因子和粗糙度矩阵的引入；非线性问题下最优化过程中灵敏度矩阵的推导 及深入分析；反演形成的大型病态线性方程组的求解，重点介绍共轭梯度法和奇 异值分解算法；理论数据和野外实践数据的反演效果图。 3 、与有限单元法相结合的偏导数的解析法求取，该方法的可行性论证、理 论推导、编程实现及相关误差分析，与相关商用软件或前人工作结果的对比分析。 4 、结论。 第二章点源二维有限元法的正演模拟 电阻率正演计算需要解决的问题是己知研究区域的二维、或三维电性( 电阻 率) 分布，求解地面观测的电位( 差) 或视电阻率数据。由于地下介质常常是不均匀 的，因此该问题不存在解析解，只能采用数值解法求解，有限元法和有限差分法 是目前常用的主要数值方法。最早将有限元法用于电法勘探的是j h c o g g o n 7 1 ， 1 9 7 1 年他从电磁场总能量最小原理出发，建立了用有限单元法进行电和电磁模 拟的算法；1 9 7 7 年l r r i j o t 8 】进一步完善了它的算法，使之成为计算二维地电条 件下电阻率法和激发极化法异常的有效方法；1 9 8 1 年d e p r i d m o r e t 9 】等发表了用 有限单元法作三维电法和电磁模拟的研究成果；周熙襄等f l o 1 9 8 3 引进 c o g g o n - r i j o 的算法及d e y 和m o r r i s o n 于1 9 7 9 用于有限差分模拟的混合边值条 件，罗延钟等】选用边值条件和反傅氏变换的算法及波数取值做了改进，使整 个算法更臻完善。 有限单元法( f i n i t ee l e m e n tm e t h o d 或简称为f e m ) 是将要分析的连续场分割 为很多较小的区域，它们的集合代表原来的场，然后建立每个单元上待求场量的 近似式，再结合起来进而求得连续场的解。从数学角度上来讲，它是从变分原理 出发，通过区域剖分和分片插值，把二次泛函的极值问题化为多元二次函数的极 值问题，后者等价为求解一组多元线性代数方程组，是一种从部分到整体的方法。 在目前常用的这两种数值方法中，有限元法可以模拟任意地形和复杂的电性 分布，缺点是计算量大、内存要求多；有限差分法求解简单但精度与差分步长有 关。因此，将采用有限元法求解点源二维电阻率正演问题。 2 1有限单元法的基本原理 2 1 1 基本方程 ( 1 ) 位函数所满足的微分方程【1 2 1 3 】 在三维笛卡尔坐标系( x ，y ，z ) 中，假设电源供电点a ( ，y ，z 。) 位于地表 时，电流强度为i ，产生稳定电流场的位函数u ( x ，y ，z ) 满足三维微分方程： v ( c r v u ) = - 2 舾( x 一工彳) 艿( y y ) 万( z z _ ) ( 2 1 ) 对均匀介质，即电导率叮为常数。上式便为泊松方程： v 2 u = 2 ，p 万o 一) 艿( j ，一儿矽0 一乙) ( 2 2 ) 若在无源空间，上式变为拉普拉斯方程： v 2 u = 0 ( 2 3 ) 对于三维地电条件，电导率o r = o r ( x ，y ，z ) ，电位u = u ( x ，y ，z ) 。设地下介质 的导电性沿走向( y 轴方向) 无变化，即o r ( x ，y ，z ) = o r ( x ，z ) ；点电流源a ( i ) 位于 坐标平面x o z 上，这时电场仍然是三维的，电位函数u = 陬n z ) 所满足的微分 方程变为： 珈琊，掣愕卜z ，掣 + 昙k 力掣i ：一2 1 8 ( z 巩) m ) m 一) ( 2 - 4 ) 为了消除沿走向坐标y ( 使变为二维问题) ，故对上式两端分别沿y 轴方向 作余弦付氏变换，上式变为： 兰修娑) + 昙p 芒) 一2 ：o r v ：( 2 - 5 ) 2 z o r vfi ( 盯i ) + i ( 盯i ) 一 = 式中：仃= o r ( x , y ，z ) ；f = 一协( x 一) 万( z z ) ；z ( x ，五，z ) 为电位u ( x , y ，z ) 的付 氏电位，这样便将三维偏微分方程变成了二维偏微分方程。 ( 2 ) 位函数所满足的边界条件【1 0 1 由于电法勘探研究的稳定电流场分布于整个地下半空间，为了减少计算量， 在解正演问题时，通常把计算范围限定在一个有限的求解区内，这样便需要在求 解区的边界上对付氏电位v ( x , 2 ，z ) 给定定解问题的边界条件，地表面r l 上( 如 图2 1 ) ，电场强度的法向分量为零，故其相应的边界条件为： 绝曩 r x n oo r z avj：o(2-6) o n i t , 对于其它边界r 2 上，根据不同的假设，分别有三类边界条件： 第类：v l r 2 = 0 ( 2 7 ) 第二类：掣l ：o ( 2 - 8 ) o n i r , 第三类- o v + q v l r = 。 ( 2 - 9 ) r 其中玎：名k ( 万) e o s ( f 元) k ( 符) ，k ，。k n 分别为一阶、零阶修正贝塞尔 雨数。根据前人的研究，采用第一类边界条件时，计算的v 和t 落己e 霉谵值偏小； 采用第二类边界条件时，计算的v 和u 值比理论值偏大。于是周熙襄等人采用第 三类边界条件计算v 值，获得了比采用前两类边界条件更高的精度。但当电极距 较大时，计算的v 和u 值与理论值仍然有一定的偏差，故罗延钟教授又对混合边 值条件式( 2 - 9 ) 中系数1 1 的取值进行了改进，其表示式变为： r|n r = 五 厶k i ( 饥) c o s ( 五酬，t k o ( 饥) ( 2 1 0 ) k = l女= l 式中：j 。一第k 个供电电极的电流强度； 五一第k 个供电电极的到边界点的矢径； 但对于电阻率成像而言，上述几类边界条件的v 值都有一定的偏差，为了提 高模拟精度和电阻率成像的反演效果，文献作者提出，将式( 2 - 9 ) 作进一步的修 正，即采用更为一般的边界条件，即修正为： a v ( k 三+ 玎y ) l = 0 ( 2 - 1 1 ) o n l r 2 k 为边界修正参数。它是一个与装置形式、网格剖分和地电断面等有关的一个参 数。根据文献作者试算结论，对于双异性点电流源施仑贝尔装置，k 一般在0 8 1 2 之间，均值为1 0 。但由于野外地质断面的特征千变力化，如何根据装置类 型和网格剖分形式的不同来选择k 值，还有待于进一步研究。 2 1 2 稳定电流场的边值问题 在二维地电条件下，由上述微分方程式( 2 3 ) 和边界条件式( 2 _ 4 ) 所组成的，对 若干个给定波数五求解点电流源场的傅氏变换电位y ( 五五，z ) 所满足的如下二维 偏微分方程的边值问题： 兰( 盯尝) + 兰p 娑) 卅盯y ：厂 q 班次o z院 娑：o l o n 娑+ c o s ( r , n ) v = 0 r 。 d 押 r 【2 1 2 ) 式中厂= 一厶艿“一，7 - - z k ) ，厶为第k 个点电流源强度，r j 和r 。分别为 研究区域的地蓑边界和地下无穷远边界，n 为地下无穷远边界的外法线方向，仃 为电导率，r 为点电流源a 到边界的距离，i 是供电电流。由于研究区域内部的 电性差异界面满足自然边界条件，故在上述边界条件中不需考虑。 亟：! ：堂丝途塞箜三重 盛塑= 三丝鱼匣五这曲正这搓塑 2 1 3 稳定电流场的变分问题 与二维偏微分方程边值问题等价的变分问题为： ，c 矿，= 蚁盯 ( 芸) 2 + ( 警) 2 + 旯2 矿2 + 2 y ) 凼 + ，r ：盯彳y 2 d = m i n ( 2 1 3 ) 求出变换电位矿( x ，五，z ) 之后，便可按下式作傅罩叶逆变换计算电位 u ( 圳，z ) = 吾j c o 哪，彬) c o s ( 纠掀( 2 - 1 4 ) 2 1 4 区域离散及网格剖分 在用有限单元法求解电法勘探正演问题时，由于受计算量和存储量的限制， 需要将无限地下半空间中的电场分布限定在有限的求解区内，并将连续的求解区 和变分问题进行离散化，即网格剖分。 ，- ， ， x ， ， ， ，7 ， ， ， j - ， j - ，一， ， j - ， - - 。，。 ， ，。 - - ， ，j ， ， ， ， ， ? r p - j ，_ ， r p _ ， ，7 j ， j ， ， ， ，。， ， ， ， ， j ， ，二外区 ，j， ， ，7 ， ， ，。 ， ， ， ， ， ， ， jj， ， ， ， ， z 图2 - 2 网格剖分示意图 1 、网格划分基本原则： 对于二维电法问题，最常用的是三角单元，这种在矩形网格中布置交叉对称 三角形的剖分网格可以足够近似地模拟一般常见的不平地形和电性异常体，同时 又能节省计算量。根据有限元法及二维电场的特点，网格剖分应注意以下几个基 本原则： 各三角单元不能重叠、不能有公共内点。 网格剖分要遍及整个待解场域，当边界为直线时可直接取为线元，当边 界为曲线时，应先将曲线边界分段线性化后再取为线元。 网格剖分愈细，计算精度愈高，但计算量也愈大，所以在满足一定精度 要求下，应尽量减少网格单元的数目。 对非均匀介质模型，每个三角单元中只能含有一种介质，导电率为常数。 为节省计算机内存，网格节点编号应使有限元方程的总系数矩阵的带宽 尽可能小。 在电场变化剧烈的地方，网格应划分细一些，而在电场变化平缓的地方， 网格应划分的粗一些，这样可以减少单元总量。但网格的疏密应是逐渐过渡的。 ( 2 ) 网格编排： 节点编号按从左向右，从上到下依次编排。网格剖分方式不同，数值模拟的 精度也不同。由于泛函积分区域是网格区域，所以一般的网格区域越大，并且网 格区域边界远离不均匀体，都能得到可靠的模拟精度，但增加网格区域势必导致 节点数增加，相应的也将占用大量的计算机内存和增加了计算量。对于网格单元 来说，单元越小，计算精度也越高，然而网格单元越小，网格剖分节点数就会成 倍增加，将使微机的计算效率大大降低。为克服网格剖分区域、剖分单元的大小 和模拟精度及计算量之间的矛盾，本文选择了图2 - 2 的网格剖分形式。在成像区 ( 粗线区域内) ，采用较小的均匀网格，在外区采用较大的不均匀网格，并且由内 到外逐渐增大网格单元。采用这种剖分形式，即克服了上述矛盾，又兼顾了了稳 定电流场的分布规律，完全可以满足电阻率成像反演的需要。 2 1 5 线性插值 为了计算二次函数j ( v ) 需要知道求解区域内的v 函数值，通常利用各节点 函数值在各个单元内作线性内插来求v 值，如图2 - 3 所示三角形单元。 ， 朋 图2 - 3 有限元网格剖分中的三角形单元 设第e 个单元三个节点按逆时针方向编号，依次y , ji ，j ，l q l ，其坐标为 “，毛) ，( 而，乃) ，( ，乙) 对应的节点函数值为k 、巧、圪，单元内函数y 是线性 变化的： 即： y ( z ) = a + 缸+ 口( 2 - 1 5 ) 单元r 4 线性插值函数v ( x ，z ) 就可以近似表示为三个顶点的傅氏电位k 、匕、 圪的线性函数： k = 口+ 魄+ i 巧= 口+ + q ( 2 _ 1 6 ) 圪= 口+ b x + c z mj 用克莱姆法则求解上述方程组并整理可得单元e 的函数值： v ( x ，z ) = m “z ) k + ，( 了，z ) 巧+ 虬( x ，z ) 圪 ( 2 1 7 ) 式中，m ，帆为形函数，其表达式为： 1 1 m 2 五( q + + q 弓) m = 去( 乃+ t _ + c j z j ) 虬= 瓦1 ( 口辨+ 吃+ 气乙) 其中，p 为单元e 的面积： 口，= x s z 一x , z s ，岛= z s 一，c ： = 一t a s2x m z i x | z m ，b j2z m z , , c j 2 x i x m a ：= x t z j x j z , ，b ：= z l z j ，c ：= x j x i ( 2 1 8 ) ( 2 - 1 9 ) 将每个单元都做出这种近似表示式，就得到整个求解区域内v 的总体近似函 数，此函数在各单元内是线性的，对任意的两个单元来讲，近似函数在公共边上 的值被两端点的节点函数值唯一确定，故总体近似函数在整个求解区内是连续 的。 2 1 6 变分问题的离散 二维变分问题中的泛函j ( v ) 是对整个求解区域积分的，它可表示为各个单 元的积分以( 矿) 之和： ，( 矿) = 以( 功 ( 2 - 2 0 ) 式中j ：( v ) 可l l i l t t i l 各单元e 内函数v 的线性插值近似式求得。 ( 1 ) 单元分析 现讨论正( y ) 对v 的一阶导数， 哿= 2 口 警壶( 篆) + 警参( 警) + 2 簧一2 孑 砒 + 吖譬矿挚 ( 2 - 2 1 ) 由娑，i a v 的表达式以及( 4 1 7 ) 式可求得： o 眩 m ( j ，z ) = 心y t 七ny j n 3 n 。i ( 2 - 2 2 ) 将上式代入( 2 - 2 0 ) , 并考虑到血j ：f d x d z = 血和厂2 秒i 黼对( 2 - 2 0 ) 中 括号中的被积函数积分整理得下式： 恚 ( 留+ 0 心+ ( + 巳c ) “，+ ( 包+ q 气) ”， ( 2 - 2 3 ) 被积函数第二项为： ( 2 - 2 4 ) 驴厂簧砒2 旷袱卜，卜m 砒 ( 2 - 2 5 ) 对于线积分项，当以j ，m 节点代表边界节点，因此j ，m 节点线积分可整 理为： 2 吖护挚= 2 叭c 扣圳 2 吨辛划c “圪) q 乏6 式中，= f ( _ 一) 2 + ( z j - - z m ) 2 1 2 为j ，m 的边长。 总结以上公式，讨论券的各项计算公式，对j ，m 可类似推导得出公式。加上 对边界单元的讨论，可以归纳如下：对与任何一个小三角形单元警的计算有 9 = y k 矿一k a a a a 矿 ，-il-_i_j-_lil-l + 哆 + k仁 丝6 = 砒 y k a a 口 舻煮项 三 第数函积被 = 眨 式中单元系数矩阵的个元素为： 心= 乏( + 俐+ 三腽 屯= 恚( + c f j ) + a e + 詈吖 k = 毫( + 气) + l , a a e + 詈吖 毛= b = 乏( 岛岛+ q c j ) + l f l a p k = 吒，= 乏( 岛k + q ) + 否1 肚 = = 盖( 乞k + 巳勺) + a e + 三吖 ( 2 ) 总体合成 ( 2 - 2 7 ) ( 2 - 2 8 ) 将d 域中各单元的系数矩阵元素计算出来，再进行对应元素的叠加，可得到 总体系数矩阵并建立方程，将与电流i 有关的项移到右端，得到将所有单元的 j , ( r ) 相j i ，得整个求解区域泛函j ( v ) 。 笔 曼 = 三 。2 2 9 ， 2 1 7 刚度矩阵系数方程组的解法 在点源二维地电条件下，对于给定若干个波数入，求解电位u 的付氏电位 v 所满足的二维偏微分方程的边值问题为： + 昙( 娑) 一) 1 , 2 c r v o r c r v ： + l 1 一 2 彪、玉7 t 筹州卜 伽 l r ( 2 3 0 ) 一、，一 ，0 0 ，l = i、，一k _ ，。、，跏枷 ， 舭一职觇一嵋舭一叱 k ： 鼬助：伽，。-一 警卸 口 忆 ，l 枣竺锄 式中= 一j 。8 ( x - x k ) 8 ( z - z 。) ，l k 为第k 个电流源强度。与其等价的变 分问题为： - ，c 矿，= 三仃 ( 篆) 2 + ( 警) 2 + y 2 + 厂y 如 + j 1f 2 叼y 2 刃= m i n ( 2 - 3 1 ) 在用有限元法作正演模拟时，要将整个求解区的泛涵j ( v ) 问题( 2 - 3 1 ) 转化 为各个剖分单元的积分之和。即 j ( y ) = j 。( y ) ( 2 - 3 2 ) e 然后对各个单元泛涵以( 矿) 求极值，得到单元刚度矩阵，再将各个单元刚度 矩阵合成总体刚度矩阵，进而形成以节点函数值为未知数的方程组： k 】【明= 【i 】( 2 - 3 3 ) 其中【k 】为总体刚度矩阵；【v 】为待求的付氏电位；【，】为与供电点有关的列 矢量。由变分问题离散化所得到的线性方程组的系数矩阵( 即刚度矩阵【k 】) 具有对称、正定等特点：还有刚度矩阵【k 】的阶数等于求解区的总点数，通常求 解电法勘探的正演问题时，方程组的阶数可达到n * 1 0 2 - - - n * 1 0 3 ( 甚至更大) ，所 以瞵】是大型的；由于一般情况下求解区沿水平方向较长，节点数目( n ) 较多， 而沿垂直( z ) 方向较短，节点数目( m ) 较少，故将节点编号从左边开始，逐列 从上至下进行编号，使中心节点与相邻节点的编号的差值尽量小，则矩阵【k 】的 非零元素将比较集中的分布在对角线附近的带状区域内，并且矩阵【k 】的半带宽 w 最小为m + i 或m + 2 ( 这主要与网格的剖分形式有关) 。 由于用有限元法求解电法勘探正演问题时，大部分计算量都是花在解线性 方程组上，所以恰当的选用求解线性方程组的方法是十分重要的。在有限元法中， 常用到的求解线性方程组的方法有直接法和迭代法，但是在解决电剖面法或电测 深法的正演问题时，需要计算同一断面不同供电位置的电场分布，供电位置改变 时，只改变方程( 2 3 3 ) 的右端场源项【，】，而其系数矩阵【k 】保持不变。这样在计 算不同供电位置的节点的函数值时，只需形成和分解一次刚度矩阵【足】，通过回 代就可计算不同供电位置下的所有节点的函数值，这样可以节约很多计算量，所 以选用直接法进行求解是比较经济的。 a = 仃。一厂一 - j _ 一 l 图2 - 4 原矩阵( a ) 及紧凑压缩存储格式( b ) 示意图 ，s x n 实际计算时，矩阵【以】是按以下顺序存放在内存中，即： ( 口i i ，a 2 2 ，口肌，id 2 1 ，a 3 2 ，口 ，l i ，0i ，ia s , i ，口j + i 2 ，口疗 j l 一5 一i ，o ，o ) 这时，矩阵 彳】的非零元素口j ，与 以】的元素口5 j ，_ ，之间有如下对应关系： a i 嘞：卜川沦， ( 2 - 3 4 ) , j2 口_ ，2 1 口，产件l f f _ ， 2 。) 4 对于对称j 下定方程组： a x = b ( 2 3 5 ) 采用c h o l e s k y 分解法进行求解1 。以压缩存储为依据( 为使算式直观，暂 不考虑二维存储a ) ，设置二维数组a 和一维数组b 与t ，得到如下比较实用的递 推算法： 对于i = 2 , 3 ，刀计算到第步 f _ ，卜一i - “ i k = l n u 卜t j | q 珏 卜一口诸 i = i i - i 岛卜6 j 一b l ( 歹= 1 , 2 ，f 1 ) 求“ ( _ ，= 1 ，2 ，i 1 ) i v 求4 ( f = l ，2 i9 订) i 荥y i 一 岛卜6 f a 豇一口矗( f = 刀，n - 1 , - - - , 1 ) 求毛 k = i + l 最后在数组b 中得到a x = b 的解工 对于上述二维压缩存储格式，虽然存储空间减少了一半多，但是半带宽矩 阵中的零元素仍然占据了相当可观的存储空间( 约r * ( 3 - - 4 ) 个) 。可以看出，随 着剖面长度的增加，即水平方向网格剖分数越多，零元素所占的额外空间越多，所 以考虑采用一维压缩存储格式，仅存储半带宽矩阵中的非零元素，对上述存储方 法作进一步的改进是很有必要的。 2 1 8 反傅氏变换的计算 在对若干个不同的波数名值司n 刚n 7 1 分别求出各节点的傅氏电位v 后，对其进 行反傅氏变换，可得各节点电位u 值。通常在求解地面电法的正演问题时，并不 需要确定求解区所有节点的电位，只要计算供电点周围一定范围内地面节点的电 位值。实际上只能对有限个离散的波数五值求出节点函数值v ，计算v 的波数越多， 分布越密，则由v 计算u 的近似性越好( 精度越高) 。因为改变波数五时，刚度矩阵 也随之改变，故对于每一个波数五值，都需要重新形成和分解刚度矩阵及重新作 回代，以计算相应的节点函数值v ，因此计算量近于随所用波数的个数成正比地增 大。所以反傅氏变换问题归结于用尽量少的波数及相应的节点函数v 值，计算出 符合精度要求的节点电位u 值。通常只计算和供电点在同一断面内之节点电位值 y = 0 9m u ( x ，0 ，z ) = 二【v ( x , a ，z ) d 2 ( 2 3 6 ) 石” 因函数v 随波数五迅速衰减，我们可用数值积分方法对有限个离散波数a 求出节点函数值v ，从而得到上式的近似值。 2 2 有限元法的特点和优点 有限差分、有限元和边界积分三种数值模拟方法在电阻率正演计算方面都有 一定的优势和不足；有限差分法的优点是方法简便易算，其缺点是，当物性参数 复杂分布或场域的几何特征不规则时，适应性比较差。而边界单元法的优势是正 演速存需求少，主要用于地形改正和地下少量地质体的正演模拟。有限单元法与 前述方法相比，在电阻率法j 下演方面有独到的优势： ( 1 ) 把二次泛函的