（地球探测与信息技术专业论文）三维电阻率测深数据近似反演研究.pdf

桂林工学院硕士学位论文 摘要 详述了三维电阻率测深数据有限元正演计算方法和步骤，并简单介绍维 z o h d v 近似反演方法的原理和基本思想。针对三维电阻率测深数据采用z o h d y 方 法进行近似反演，编制了三维电阻率测深数据近似反演程序。初始模型由测量视 电阻率数据给出，通过比较实测视电阻率和预测模型计算的视电阻率对数差来修 改模型网格电阻率 为了解决任意电极距测深数据的反演，采用大、小双网格剖 分。大网格反映地下电性分布情况，小网格用于实际有限元正演计算。讨论了标 准的z o h d y 方法在反演过程中存在的缺陷并对其进行了改进，用不同深度转换因 子反演一维h 型断面和三维低阻凹陷模型理论数据，验证了0 5 作为深度转换因 子反演模型能快速收敛；在电阻率修改公式中加一个迭代系数，反演的收敛速度 明显加快；在反演前对视电阻率数据进行水平滤波，反演含不同随机噪声的三维 电阻率测深数据，结果显示该方法可以压制部分噪声，对离信噪比的数据反演能 够稳定收敛。r 对几个不同类型加5 随机噪声的三维地电模型理论视电阻事测深断面数据 和实测数据进行反演，反演得到的结果表明改进的z o h d y 方法反演三维电阻率测 深数据是可行的，所编制的程序是正确的。该反演方法需要的内存小、计算量少， 而且在迭代初期误差下降较快，可以用较少的迭代次数就可以得到与实际模型十 分相近的电阻率分布。但是近似反演是一种直接拟合方法，反演得到的电阻率分 布与实际地电模型仍然存在差距，反演复杂的和信噪比低的视电阻率断面数据时 效果较差，也不能用于地形起伏三维电阻率测深数据的反演，应用于野外实测数 据的反演还有一定的差距。 关键词：三维电阻率测深数据；z o h d yj j 去 数 r p 一 深度转换因子；水平滤波；迭代系 u 一 桂林工学院硕士学位论文 a b s t r a c t 1 1 1 ea p p r o a c ha n dp r o c e d u r ef o rt h r e e - d i m e n s i o n a lr e s i s t i v i t yf o r w a r ds o l u t i o n b a s e do nf i n i t ee l e m e n tm e t h o da r es p e c i f i e d a n dt h ep r i n c i p l ea n df u n d a m e n t a li d e a o fz o h d yi n v e r s i o na r ed e s c r i b e & a na p p r o x i m a t em e t h o dw h i c hc a l lb cu s e dt o i n v e r s et h r e e d i m e n s i o n a lr e s r i v i t yd a t ai sp r e s e n t e d a n dt h ec o r r e s p o n d i n g p r o g r a m i sw o r k e do u t 1 1 1 em e t h o di sb a s e de l lf i n i t ee l e m e n tf o r w a r ds o l u t i o no ft h ee l e c t r i c a l f i e l do ft h r e e - - - d i m e n s i o n a le a r t h t h ei n f t i a le l e c t r i c a l p a r a m e t e r so ft h e e a r t ha r e o b t m n e df r o mo b s e r v e dd a t ao rf o r w a r dd a t ao fe a r t hm o d e la n dm o d i f i e db yt h e l o g a r i t h i cd i f f e r e n c eb e t w e e no b s e r v e dd a t aa n dt h et h e o r e t i c a ld a t ac a l c u l a t e db yt h e f i n i t ee l e m e n tf o r w a r ds o l u t i o n t w ok i n d so f 鲥d s y s t e mi se m p l o y e dt os o l v et h e p r o b l e ma r o s ef r o mt h ea r b i t r a r ye l e c t r o d e - s p a c i n g 伯eb i gg r i dp i c t u r e sr e s i s t i v i t y d i s t r i b u t i o no fe a r t hm o d e la n dt h es m a l lo n ei su s e dt od e f i n ee l e m e n t si nf o r w a r d s o l u t i o n 佻es h o r t c o m i n g so fs t a n d a r dz o h d y m e t h o d , s u c ha sl o wc o n v e r g e n c ea n d u n s t a b l es o l u t i o n , a r ed i s c u s s e da n di m p r o v e d 1 r h ee l e c t r o d e - s p a c i n gt o d e p t hs h i f t f a c t o r s r a n g i n g f r o m 0 3t o0 6a r e a p p l i e d t oi n v e r s ea h - t y p ep r o f i l e a n d t h r e e d i m e n s i o n a lv a l l e ym o d e lt h er e s u l t ss h o w0 5i sa p p r o p r i a t ev a l u eo ft h es h i f t f a c t o r i na d d i t i o n , t h ei t e r a t i o nc o e f f i c i e n ti sa d d e dt or e s i s t i v i t ya d j u s t m e n tf o r m u l a s ot h a tt h ec o n v e r g e n c e v e l o c i t yh a so b v i o u si m p r o v e m e n ti nt h ec o u l s eo fi n v e r s i o n n e p l a n ea v e r a g eo fr e s i s t i v i t y d a t ab e a r i n gs o m er a n d o mn o i s es h o u l db et a k e n b e f o r ei n v e r s i o n t h es e v e r a lt h r e e d i m e n s i o n a le a r t hm o d e l sa r e d e s i g n e dt o t e s tt h ei n v e r s i o n p r o g r a m b e f o r ei n v e r s i o n , t h e5 r a n d o mn o i s e i sa d d e dt ot h et h e o r e t i c a ld a t a o b t a i n e df r o mt h em o d e l s t h er e s u l t sc o n f i r mt h a tt h em e t h o di se f f e c t i v ea n dt h e p r o g r a mi s c o r r e c t a n dl e s ss t o r a g ea n dc o m p u t a t i o nt i m ea r er e q u i r e dt h a no t h e r i n v e r s i o nm e t h o d s i nt h ei n i t i a l s t a g eo fi n v e r s i o ni t e r a t i o n s ，t h er o o tm e a ns q u a r e d e l t o rd e c r e a s e s r a p i d l y t h e r e f o r e ，a f t e rs e v e r a li t e r a t i o n s ，t h ei n v e r s i o nr e s u l t s a p p r o a c h t h er e a l r e s i s t i v i t y d i s t r i b u t i o n v e r yw e l l b u tt h ez o h d ya p p r o x i m a t e m 桂林工学院硕士学位论文 i n v e r s i o ni sad i r e c ts o l u t i o ns ot h a tt h e r ei sd i f i e f e n c eb e t w e e nr e s u l t e d1 1 1 0 d e i sa n d r e a lo n e s m o r e o v e r , t h em e t h o di s n tf i tf o rd a t ao f c o m p l e x e a r t hm o d e l so rl o ws n ， a n dc a n ti n v e r s et h e s eo b s e r v e dd a t ao f g r o u n du n d u l a t i o n i naw o r dt h e r ei ss o m e w a y t op u tt h i sm e t h o di n t op r a c t i c a lu s e k e yw o r d s ：t h r e e - d i m e n s i o n a lr e s i s t i v i t ys o u n d i n gd a t a ；z o h d ym e t h o d ；p l a n e a v e r a g e ；t h ed e p t h s h i f tf a c t o r ；i t e r a t i o nc o e f f i c i e n t i v 桂林工学院硕士学位论文 独创性声明 本入声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研 究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他 入已经发表或撰写过的已经成果，也不包含为获得桂林工学院或其它教育机构的 学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已 在论文中作了明确韵说明并表示了谢意。 签名：强强磕日期：丝里。! ：兰 关于论文使用授权的说明 本人完全了解桂林工学院有关保留、使用学位论文的规定，即：学校有权保 留送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分 内容，可以采用影印、缩印或其他复制手段保存论文。 ( 保密论文在解密后应遵守此规定) 桂林工学院硕士学位论文 第1 章引言 电阻率法是以地下不同岩矿石的导电性( 电阻率) 差异为基础，通过观测与 研究人工建立的地下稳定电流场的分布规律，以达到找矿和解决其它地质问题目 的的一种勘探方法【4 1 】1 4 2 1 1 4 3 。在普查金属、非金属矿产和研究地质构造方面都取 得了很好的效果，随着观测方式的不断改进、测量精度的提高，特别是高密度电 测系统的出现，使得电阻率法的效率更高、获得的数据点更多，开始在浅层的水 文、工程、环境、考古等与国民经济建设、人类社会生活密切相关的领域得到越 来越广泛的应用1 4 4 1 1 5 5 1 。 电阻率法是通过在目标区测量视电阻率值来获得其内部结构的电阻率( 导电 率) 分布情况，其中如何从外边界测量的视电阻率值来得到内部电阻率值( 电阻 率反演) 是该方法的核心问题。电阻率反演也是地球物理勘探中最早应用计算技 术的方法之一。1 9 3 3 年，l b s l i c h t e r 提出了电测深血线的“直接解释法”，但由 于受计算技术的限制，未能在生产实践中推广使用。计算机的出现给电阻率法提 供了计算条件，开始用量板法对视电阻率曲线进行一维解释p 7 】。随着计算技术的 不断提高，电阻率反演方法经历了一维、二维到三维的发展过程，在早期基于核 函数递推视电阻率正演计算方法基础上，研究出一维拟合反演方法，并实现了电 阻率的自动反演。随着勘探实践的深入和生产需求，要求的地电模型更加复杂， 一维电阻率反演难以满足要求。计算技术的飞速发展，能够利用很多方法进行电 阻率二维正演模拟，这样电阻率反演从一维过度到二维。p e l t o n 6 1 ( 1 9 7 8 ) 等提出了 一种确定地下二维模型网格电阻率、极化率大小的最小二乘最优化反演方法，之 后s a s a l 【i 【l l 】、s m i t h 9 1 、t r i p p t 7 1 等) k 发展了该方法，使得电阻率法更加实用。但 是地下探测目标多表现为三维电性结构，s p i t z e r 1 3 1 ( 1 9 9 8 ) 通过电阻率法灵敏度矩 阵揭示出：对三度体用电剖面数据做二维解释不可避免会有偏差。进入2 0 世纪 8 0 年代，三维直流电阻率法正演计算的实现，特别是有限元计算程序为三维电 阻率反演打下了基础。 地球物理反问题是不适定的，二、三维电阻率反演更是如此1 4 9 脚m 】。反演 过程基本是依赖于正演计算，正演是基于有限差分或有限元一类计算方法的数值 桂林工学院硕士学位论文 解，反演是求数量巨大的网格单元电阻率值。反演过程中不可避免地要遇到如何 求偏导数矩阵、如何解决计算速度慢及庞大计算内存需求的问题，另外网格单元 太多时反演常表现出不稳定性和高度非唯一性。随着计算机和数值模拟技术的发 展，近年来，已有不少三维电阻率反演成果发表。 电阻率三维反演国内外起步不久，主要有【1 6 】【l7 】：基于b o m 近似的三维反演、 电阻率反投影技术( 电阻率层析成像) 、t o r a n t o l a 反演以及传统的最d , - 乘反演 等方法。b o r n 近似反演认为二次场仅是一次场在电性不均匀处产生的积累电荷 激发，面忽略积累电荷间的相互影响。该方法需要电阻率背景值的估计，b o r n 近似反演只有当地下介质导电性均匀时，近似的精确程度较高，反之，介质电阻 率变化越大，精确程度越低。可见b o r n 近似反演法分辨率低，应用于复杂电性 结构的反演受到限制。电阻率反投影技术是根据测区周围的各个不同方向视电阻 率差来重构探测区内的电阻率分布，反演的介质电阻率仅是实际电阻率的一级近 似，分辨率较低，一般作为其它反演方法的初始模型。虽然经过许多人( 如y o r k e r 、 d a i l y ) 的改进，但该方法另一限制是要求在目标地质体的周围都有观测数据点， 才能保证反投影过程的稳定可靠。t o r a n t o l a 反演方法是一种基于概率分布的方 法，模型修改过程中需要给出协方差矩阵和先验模型。直流电阻率法测量一般不 能提供协方差矩阵，即使给出一个先验模型，由于三维网格单元太多，以及其电 阻率值变化可能达到几个数量级，又不能给出合理的协方差矩阵来限定三维反演 迭代过程中模型的修改，使反演收敛到可靠的解也是比较困难的。传统的最小二 乘反演方法是地球物理反演中的高精度方法，该方法大大提高电阻率反演的精 度。但是最小二乘反演不可避免雅戈比矩阵的求取和大型矩阵的计算，无论从计 算速度还是机器内存需求而言，均只能适宜网格单元数量较少的模型反演。近几 年，国内不少研究者对上述方法进行了改进，如毛先进1 3 8 】应用积分方程法进行电 阻率三维反演，吴小平hs 】1 3 9 】利用共轭梯度法的电阻率三维反演，这些方法虽然降 低了计算量，但反演的效果难以令人满意。 1 9 8 9 年，z o h d y l l 】提出了一种直流电阻率测深数据的拟合反演方法，这种方 法不需给出初始模型，初始模型由测深曲线直接给出，通过模型计算的视电阻率 与实测视电阻率的差异比较来修改模型：她不需要计算雅戈比矩阵，计算量少， 需要的内存少，稳定性好。该方法已经在一维、二维电阻率法反演中取得了良好 的效果。本论文将用z o h d y 法对三维电阻率测深数据进行近似反演。 2 桂林工学院硕士学位论文 本文主要进行如下几个方面的研究和工作： 基于三维有限元正演计算程序，编写了三维电阻率测深数据z o h d y 近似反演 程序，采用大、小双网格系统剖分单元，分别用于反映目标区电阻率的分布和实 际的视电阻率测深数据的正演计算，来解决任意极距情况下多测深点数据的正演 计算。 针对地球物理反演过程中的收敛难和稳定性差的缺陷，对z o h d y 近似反演方 法进行了讨论和改进： 1 ) 用不同深度转换因子来进行反演计算和收敛速度的对比，说明选取0 5 作为转换因子，反演基本稳定收敛，而且得到的结果较好。 2 ) 为了加速反演的收敛速度，在标准的z o h d y 反演方法调整模型网格电阻 率的公式中加入一个迭代系数，使模型计算的视电阻率值快速趋于观测值， 用较少的迭代次数就可以得到很好的反演效果。 3 ) 针对含噪声的观测数据，引入水平平均( 滤波) 方法，可以在一定程度 上压制噪声，保持反演的稳定。 用该程序对几个基本模型的数值模拟数据和实测数据进行了反演，说明 z o h d y 近似反演方法可以用于三维电阻率测深数据的反演，计算量少，稳定性好。 但是由于z o h d y 反演方法是直接对视电阻率数据进行拟合，难以用于复杂的地电 模型反演，要求观测的数据必须具有高信噪比，且模型的视电阻率断面不能够太 复杂。 桂林工学院硕士学位论文 第2 章三维电阻率测深数据的正反演 地球物理数据的计算存在两类基本问题：正演计算和反演解释。正演是给定 场源及地球物理模型参数分布，求解场值的大小，基本途径有 4 5 】：1 ) 解析法， 只适合简单数学物理方程的模型计算，2 ) 模型试验法，各种因素影响较大且复 杂的异常很难模拟。3 ) 数值模拟法，根据地球物理场的偏微分方程和边界条件， 用数值方法求场值的近似解，适合各种物性参数分布和复杂边界形状的地球物理 模型计算。对于三维电阻率测深数据正演而言，要求能够计算复杂的异常分布， 特别是用于反演过程中，近年来发展起来的有限元数值方法，具有很强的适应性， 随着计算机的迅速发展，解决了计算工作量大、存储量大的困难，是三维电阻率 测深数据正演计算的主要方法。 三维电阻率测深数据反演问题是根据目标区测量的视电阻率值来确定地下 介质的电性参数分布，其基本方法是拟合法，即对假设模型计算的视电阻率值与 实测的视电阻率值进行比较来调整模型参数，直至两者差别不大为止，在多次对 比迭代过程中一定要用到正演计算。由此可见，正演计算的好坏是反演的关键。 2 1 三维电阻率测深数据有限元法数值模拟“5 m 有限元方法不仅适合复杂物性参数分布的正演问题，其解题过程也比较规 范，是一种很好的三维电阻率测深数据的数值模拟手段。其基本过程如下： 1 、边值问题 对三维地电断面双点电源场的边值问题的求法有两种情况，一是计算总电 位，二是计算异常电位。 总电位v 满足： v ( o w ) - - - 2 1 0 ( a ) + 2 t 6 ( b ) o v ：0 ， a ，l 7 o v + c o s ( r , n ) v ：o ， o n r 4 1 2 ； r ：( 2 - - 1 ) t 桂林工学院硕士学位论文 异常电位“满足 v ( t l v u ) 一一v o “。) ， - = o u 。0 ， 0 n 竺+ 竺! 匕生“，0 ， d nr 口： f s ； ( 2 2 ) f * 式中： 口计算区域，c 为区域口的地面边界，c 为区域q 的地下边界，。为 边界的外法线方向，仃为介质的电导率，吼为电源点处电导率，仃；盯一a o 蔓t j 介质的异常电导率， u o 为正常电位，且。上f 三一三1 ，r a 、r b 分别是测点 。2 行0 0i ，= i j 到电源点a 、b 的距离。 由于总电位法正演计算时，在电源点处奇异，所以本文采用计算异常电位法 来进行三维电阻率测深数据的正演计算。 2 、变分问题 与三维双点电源场的边值问题( 2 2 ) 式等价的变分问题为： 荆- c f o 幽l o ( w 罗+ o “ e 睦掣趔“ + “。“d f ； ( 2 3 ) j 3 、有限元计算方法 1 ) 单元剖分，采用正六面体对目标区f 2 进行剖分，如图2 - - 1 所示。 5 寰 。 b 一 呵 、 桂林工学院硕士学位论文 k 图2 1 区域d 剖分厦单元示意 2 ) 插僵，采用三线性插值，即 一了n , u 4 式中，n i ( i = l ，2 8 ) 为形函数，且 i 2 砉( 1 + 笔x 1 + 卯t x 1 + 鹭；) 其中皇、叩f 、茧是点i 的坐标，亭、叩、善与戈、弘z 的关系为 x = + 兰亭，y = y 。+ 互b ，7 ， z ；z j + 兰； x o 、y o 、z o 是于单兀e 的三个中点，口、b 、c 是子单元e 的三个边长。 3 ) 单元分析 将全区域口的积分分解为单元e 的积分之和 f ( “) a 副主仃陬罗叫+ 善咖w v 枷 + 。 ( 2 - - 4 ) 瓤( 鼍掣“2 卢+ 引号掣郴卜 分别计算得： 仁仃陬) z 蛔= - 萼u r k 。u 。 6 桂林工学院硕士学位论文 伊v u 。v u d s 2 ；o u ：1 e u 扛掣2 d f = 詈u k 。u 。 f r , c r c o s ( r , n ) ) h 。耐厂；盯u ：k ：。u 。 式中u 产( “0 7 ，u o e = ( 7 ，k l 。= ( k l l j ) ，岛i 尸，k 2 产化0 ，七2 i 尸锄，, j f f i l ，2 8 ，k i 。、 k 2 e 的详细计算见参考文献【2 0 】，则单元p 内 e ) 一号h f ( k ，+ k 2 , ) u + 盯u ：( k ，+ k ：。) u 。 一去u ：k u + u r k u 0 ( 2 5 ) 一去u ：蟊u + u ：可。 其中k ；o ( k 。+ k ：。l k ：= 盯( k 。+ k ：) ，u 、嘞分别是全体节点“、 苎璧的列向量，k 。、霹分别是n ( n 为剖分的网格节点数) 阶k 。、k ：的扩 矩阵。 4 ) 总体分析 将全部单元的e 相加，得 f o ) 4 辜只o ) - 与27 氍7 辜霰：u 。一i u r k u 川k ，u 。 ( 2 _ 6 ) k 、k 为总体( 或系数) 矩阵。 5 ) 求变分 掰0 ) - d u r k u + u r k u 0 1 = o 求解得方程组 k u = 一k u 。( 2 7 ) 6 ) 解线性方稗组 7 桂林工学院硕士学位论文 上式的系数矩阵很大，为了节约存储空间，采用变带宽矩阵存储，求解方 程组( 2 7 ) 式，得到目标区域各个单元节点上的异常电位值“。 7 ) 计算视电阻率 由各个单元节点上的电位值计算任一测深点的视电阻率 舻k 型粤螋+ p 0 ( 2 - 8 ) 式中k 为装置系数，l ( 旧、叼分别是测量电极m 、上的异常电位值，p 。为 2 2 z o h d y 近似反演方法 z o h d y 方法是基于对s c h l u m b e r g e r 和w e n n e r 对称四极装置测深曲线的解释 面提出来的，该方法是一种简单的拟合反演方法。具有初始模型简单，计算方法 简单，反演后的结果接近真实电阻率的分布情况等优点，是一种非常实用的对称 四极电阻率测深数据的解释方法。 墨 器 图2 2 z o h d y 反演方法的主要步骤 a 一观测敷据和初始模型b 一转换分层和计算的测深曲线 c 一最终分层和计算的电阻率疆口深曲线结果 8 桂林工学院硕士学位论文 其基本思路：假定初始模型的电性层数等于实测视电阻率的电极距数，并始 终保持不变；初始模型的深度等于对应电极距乘以某一常数( 深度转换因子) ， 电阻率值等于相应电极距的视电阻率值；计算预测模型的视电阻率；将计算的视 电阻率与实测视电阻率进行对比，并进行电阻率调整，重复迭代直至误差达到最 小。在这一过程中需要确定深度转换因子和电阻率的调整关系。其主要步骤如图 2 2 所示 1 、深度转换因子的确定： 1 ) 设深度转换因子意产1 0 ，初始模型的电阻率值等于相应电极距的视电阻 率值； 2 ) 计算该模型的理论视电阻率测深数据； 3 ) 计算理论视电阻率测深数据与实测视电阻率测深数据的均方根误差r ： r j ； ( 2 9 ) i 一是迭代次数； 广第j 层或第j 个电极距； p 。( j ) 一第j 个电极距上的实测视电阻率值； p 。( ，) 一第f 次迭代第j 个电极距上的计算视电阻率值。 4 ) 缩小深度转换因子岛( 如置产0 9 也) ，重复2 ) 、3 ) 步，并比较均方根误 差，直至r i r 。，这时深度转换因子也为最合适的值。 2 、电阻率的调整 确定了深度转换因子之后，初始模型的深度等于对应电极距乘以深度转换 因子，电阻率值等于相应电极距的视电阻率值。 1 ) 计算该模型的理论视电阻率测深数据； 2 ) 计算理论视电阻率测深数据与实测视电阻率测深数据的均方根误差r i ； 3 ) 如果风满足要求( 如r z 2 ) ，这时模型参数( 电阻率、深度) 即为反演 的结果，否则按式2 - - 1 0 调整电阻率，重复1 ) 、2 ) 步直至均方根误差满 足要求。 n + 。( ，) = n ( _ ) + p 。( j ) p 。( ，) ( 2 - 1 0 ) 式中： 桂林工学院硕士学位论文 _ - _ _ - - - _ _ - _ _ - - _ _ _ - _ _ _ _ _ - _ _ - _ _ _ _ _ _ _ _ _ _ _ _ _ _ - - _ _ _ - _ - _ _ - _ _ _ _ 一l l l l l - - _ _ _ _ _ _ - _ _ o ， i 一迭代次数 - 第f 层或第j 个电极距 p 。( j ) 一第j 个电极距上的实测视电阻率值 p 。( j ) 一第i 次迭代第j 个电极距上的计算视电阻率值 n u ，一第i 次迭代第f 层的真电阻率值 z o h d y 对大量地电模型进行反演，发现对于s c h l u m b e r g e r 对称四极测深装置 的数据反演，深度转换因子为0 3 0 6 ，而且用某一深度转换因子( 如0 5 ) 去反 演不同的模型，大部分模型都能够收敛，只不过收敛的速度受到一定影响。因此， 二维、三维电阻率反演不需要重复迭代求深度转换因子，可以节约大部分的计算 量。 z o h d y 反演方法出现之后，由于其独特优点雨引起了广大地球物理学者的广 泛关注，并将其引入到二维电阻率测深数据的反演中，如b a r k e r 2 1 3 1 、王兴泰嘲 用于w e a n e r 对称四极测深装置的二维电阻率图像重建，阮百尧i 矧用于 s c h l u m b c r g e r 对称四极测深装置的二维电阻率测深数据近似反演，都取得了很好 的效果。 l o 桂林工学院硕士学位论文 第3 章三维电阻率测深数据近似反演方法 3 1 三维电阻率测深数据近似反演的基本步骤 近似反演的基本目的是将三维视电阻率断面数据转变成比较接近实际的真 电阻率分布断面数据。 3 1 i 初始模型的建立与正演计算 三维视电阻率测深数据的反演，一般需要正演计算多个测深点上的电阻率测 深数据，这就要求所啬畸分网格既能够反映地下电性分布情况，又能够方便有限元 正演计算。通常的做法是使用一个网格，进行细密的网格剖分，这样需要众多的 网格节点，要求计算机有较大的内存，计算量大。为了减小计算量，本次论文采 用双网格系统，用大网格反映地下电性变化和测深点的位置，小网格进行实际视 电阻率测深数据的有限元计算，来解决测线剖面方向的网格蒯分。垂直于剖面的 水平方向上，大小网格的剖分是相同的；在剖面的测线方向上，大网格均匀剖分， 小网格根据极距的变化而变化；在深度方向上大小网格可以任意的剖分，而且是 相同的。在视电阻率数据有限元正演计算时，将小网格的中心点放在大网格的节 点上来求该点的测深曲线，如图3 1 所示。 当正演计算时，其中小网格单元的电阻率数据从大网格单元中取值，由下式 给出： 2 ：dx z ( 3 - i ) p 。越 i i 鲁- i p 。 p - 4 , 网格单元中的电阻率； p j7 一大网格单元中的电阻率； n - - 剖面方向上小网格单元经过大网格单元数量； 嚣一剖面方向上小网格单元在电阻率为oi 大网格单元中所占据的长度。 初始模型电阻率值是由测量视电阻率数据给出，在一维z o h d y 反演中，深度 1 1 桂林工学院硕士学位论文 转换因子是经过多次迭代求得，但在三维视电阻率数据反演中，多次迭代花费的 时间太长，根据z o h d y 得出的结论，一般取固定的深度转换因子( 如0 5 ) ，这时 h i = o 5 + 嚣，在深度方向上，真电阻率等于相应极距的视电阻率值，采用水平平均 ( 滤波) 视电阻率测深数据，如图3 1 ( d ) 所示，即某一块区域各相连测深点 视电阻率数据为c o , k , o 、c ( + j ，毛d 、c u “j ，f ) 、c ( + j ，k + l ，f ) ，那么实际用于反演 的视电阻率测深数据是e q , g o m 式( 3 - 2 ) 给出，所以赋给网格初始真电阻率值为 ，。( 五岛，) = 工焉z ) 。 。愀，n ，幽! 地坦! ! 圣塑塑丛塑丛! 兰墨坳 ( 3 2 ) * _ _ - _ t - 一 _ _ _ - _ t - 龃) 图3 1 网格剖分示意图 a 一大网格示意图b 一小网格示意图 c 一目标区视电阻率断面和电阻率断面位置d 一测深点和实际的反演电阻率点的相对住置 在深度方向上，网格的剖分一般根据极距的变化来剖分的，( 1 ) 如果测量极 距是均匀的，那么网格的剖分按极距乘以一个常数( 深度转换因子) 。( 2 ) 测量 极距是不等距的，小极距密，大极距稀疏，这时如果按方案( 1 ) 来剖分，就会 出现浅部剖分太密，同时在深度方向上为了保证计算的精度，需要增加网格数， 从而增加计算量。因此网格仍然按测量极距剖分，由于反演深度等于极距乘以常 数，而网格是固定的，在深度方向上，反演的等效深度与网格剖分相对位置发生 痧扩 桂林工学院硕士学位论文 变化，这时相应网格中的电阻率赋值需按式( 3 3 ) 计算给出。 4h，(3-3) p 。趣 乍，p 。 p 一大网格单元中的电阻率； j d i7 一反演后对于等效深度的等效电阻率； 一大网格单元中的包含等效深度点的数量： 4 鬼一等效深度在大网格单元中所占据的长度( 深度) 。 实际视电阻率正演计算是采用小网格来进行的，小网格的电阻率从大网格中 取出，采用异常电位法三维有限元正演程序进行数值模拟，改变电源点位置，从 而获得各剖面的理论视电阻率拟断面数据。 3 1 2 反演模型修改 初始模型经过三维有限元法正演计算后，得到了理论视电阻率断面，将其与 实测的数据进行比较，由( 3 - - 4 ) 式调整该模型的每个网格单元的电阻率，根据 姐鞘对数关系崦( 甓学筹卜s ( 爱嬲) 来进行修改刺按下式来迭代修 改网格电阻率值： “舭，f ) 吲舭，f ) 爱嬲( 3 - - 4 )p d l j ，k ，1 ) f 一迭代次数： u 毛f ) 一网格的剖面方向、垂直于剖面的水平方向、深度方向的单元标号； 肛( ，k ，z ) 一第i 次迭代g 毛2 ) 单元的真电阻率； p 。( ，k ，f ) 一第g 的点上第f 极距的测量视电阻率： p d ( ，k ，z ) 一第o 七) 点上第f 极距的正演计算视电阻率。 迭代修改得到网格真电阻率值，再经过三维有限元法正演计算后，按( 3 - - 2 ) 式处理得到网格上的电阻率测深数据，如果按极距均匀剖分，则上述电阻率直接 赋给对应大网格单元，否则按式( 3 - - 3 ) 赋电阻率值给大网格中。重复迭代，直至 满足误差要求。 桂林工学院硕士学位论文 3 1 3 反演拟合判断标准 在实际工作中，用测量值和计算视电阻率值的均方根误差作为判断反演结果 是否达到要求的标志： 6l 彘萋铡咝辫) 2 ( 3 5 ) 一般情况下，6 小于给定的误差m ( 如m 取值2 ) ，即可满足要求，但由 于三维正演计算的非常耗时，迭代的次数越多，需要的时间越长，给一个预定的 误差往往需要迭代很多次，般迭代4 - 8 次基本能够大致反映电阻率的分布情 况。 实现z o h d y 近似反演的程序流程如图3 2 开始 输入实测视电阻串值p 。 对视电阻率拟断面p 。进行处理 网橹剖分及赋值 至塑计算求理论税电阻率拟断面p 。 反演求单元真电阻辜值p 桫岁 输出真电阻奉数值p i 是 图3 2 实现三维近似反演程序基本流程 1 4 桂林工学院硕士学位论文 3 2 三维电阻率测深数据近似反演方法的讨论与改进 地球物理反演普遍存在收敛速度慢、稳定性差的问题，三维视电阻率反演也 存在同样的问题，下面对从几方面对反演方法进行讨论：深度转换因子的选择、 收敛速度的改进、稳定性情况。首先给出几个理论模型，对上述几个方面进行讨 论，其中模型的参数基本如下： 大网格目标区： 测线5 条，线距l o r e ，测线长度l o o m ，点距l o r e ，1 1 个测深点测线。 小网格： 小网格上最大极距l o o m ，共1 6 个极距：1 、3 、5 、7 、1 0 、1 5 、2 0 、2 5 、3 0 、 4 0 、5 0 、6 0 、7 0 、8 0 、9 0 、1 0 0 ：沿测线方向网格节点数3 3 个，垂直测线的水平 方向网格节点数1 6 个；其深度方向按极距剖分，见图3 1 。 3 2 1 深度转换因子的选择 为了讨论的的方便，首先给出一个一维h 型地电断面模型，来说明用三维近 似反演的效果及不同深度转换因子对反演的影响。 模型1 ：一维三层h 性断面，其电阻率分布如下： h = 2 0 mp = 1 0 0 q m p = 5 0 0 q m 图3 3 模型1 示意图 考虑到这是一维地电断面，它的5 条视电阻率拟断面基本相同，所以只列出 其中一条视电阻率拟断面图，如图3 - - 4 所示。 深度转换因子的选择是影响整个反演过程的重要因素，必须尽可能地选择合 适的深度转换因子，如何把电极距与深度很好地联系起来，以便使计算得出的视 电阻率值与观测数据保持“同相”是z o h d y 反演的关键问题，z o h d y 在一维反演 桂林工学院硕士学位论文 中，在不考虑电阻率值的情况下，先对深度转换因子在一定范围内进行选择，即 寻找一个常数( 深度与极距的比例关系) ，经过一系列迭代，使迭代误差达到最 小时的常数值作为深度转换因子。对于三维反演来说，如果直接选择一个合理的 深度转换因子，可以省去一半的迭代时间。 位置( m ) 圈3 4 模型l 理论视电阻率拟断面 其实探测深度和电极距的关系很复杂，z o h d y 通过对一维模型的反演发现， 深度转换因子随模型的变化而变化，但基本分布在o 3 加6 范围内，对于大部分 模型来说使用不同的深度转换因子反演基本都是收敛的，只影响收敛速度。 e d w a r d s 通过研究发现，均匀介质中反演等效深度与对应极距的一半有关，对于 模型1 使用不同的深度转换因子迭代反演的误差随迭代次数减小的情况，如图3 - - 5 所示。 图3 5 用不同转换因子z o h d y 方法反演模型1 的误差曲线 从图3 - - 5 可以看出，z o h d y 方法反演模型得到的误差曲线都在迭代开始时 桂林工学院硕士学位论文 下降很快，随着迭代次数的增加，其收敛速度逐渐放慢。对于用不同深度转换因 子去反演，其误差曲线是不同的，深度转换因子为o 3 、0 4 时，随着迭代次数的 增加，其误差有可能增加，可能由于此时电阻率值的“相位”超前于视电阻率曲 线造成的，深度转换因子超过0 5 时，反演过程能够稳定收敛，当深度转换因子 继续增大，收敛相对减缓。o 5 是一个比较合适的值，这与e d w a r d s 所在均匀介 质中定义的深度转换因子是一致的。 模型2 ：低阻凹陷 其视电阻率断面如图3 6 ，由于给出的地电模型是对称的，所以只需给出1 、 2 、3 测线的视电阻率拟断面。 翻域 l i 。 图3 6 低阻凹陷模型示意图 1 7 桂林工学院硕士学位论文 图3 7 理论视电阻率拟断面图 a 一1 测量线的视电阻率断面，b 一2 测量线的视电阻率断面，b 一3 测量线的视电阻率断面 低阻凹陷是一个三维模型，其误差曲线如图3 8 ，其特征与模型1 的结果基 本相似，说明对于不同转换因子进行z o h d y 反演，都能收敛，其收敛速度是不同 的。对于大的转换因子( 如0 5 、0 6 ) ，其收敛速度下降慢，但能够使电阻率值 与视电阻率曲线保持“同相”，迭代多次反演仍然收敛，这对于任意模型反演时， 取一个固定的深度转换因子来保持反演收敛是至关重要的。 桂林工学院硕士学位论文 3 0 2 4 均 亲 1 8 误 差 ( ) 1 2 8 o 02 468 选代 i ：数 图3 8 用不同转换因子z o h d y 方法反演模型1 的误差曲线 b a r k e r 等人通过对二维电阻率测深数据反演，认识到z o h d y 方法使用0 5 作 为转换因子，对大量模型反演，都稳定收敛并取到了很好的效果，通过对上面两 个模型的反演也可发现，0 5 作为深度转换因子来反演三维模型也能够得到快速 收敛。 用0 5 作为转换因子，模型1 的反演结果如图3 - - 9 ，由于电测深方法的测量 1 口 深度一般只能达到等，所以只显示深度3 0 米以上的范围：反演后几条电阻 j二 率断面基本相同，只需显示任意一个断面。 图3 9 模型1 反演后的电阻率断面 1 9 桂林工学院硕士学位论文 从图上可以看出：反演的电阻率分布明显优于视电阻率拟断面，其电阻率分 布与地电断面基本一致，大致在5 2 0 米范围出现了低阻层( 5 0q r f l 1 5 0q m ) ： 但这是一个近似反演，其反演后电阻率分布与实际电阻率分布存在一定的差别， 在浅层，除了局部的毛刺外，其电阻率基本是均匀的，与实际的电阻率分布几乎 一致；到了深部，由于受上部模型参数的影响，反演是将整体来考虑，其电阻率 呈渐变的趋势，而且到了3 0 米还未达到5 0 0o m 值，这是z o h d y 反演的缺陷。 z o h d y 反演的最大优点是稳定性好，从上述两个模型可以看出，反演都能稳 定地收敛；另外反演不需要计算偏导数矩阵，其计算量明显少于其它方法，在 p e n t i u m 41 5 g 计算机上反演极距1 6 个，测点1 1 ，测线5 条的视电阻率断面，迭代 一次的z o h d y 反演时间只有几十秒，主要计算时间用在正演计算上。由于三维正 演的计算量大，花费的时间长，几次迭代一般很难达到一定的均方根误差f 如2 ) 要求，但是几次迭代( 如5 次迭代) 就能够大致反映断面电阻率的分布，模型2 经过5 次的结果如图3 1 0 ，由于模型是对称的( 1 、4 电阻率断面相同，2 、3 电阻率断砸相同) ，所以只列出1 、2 电阻率断面。 从反演的结果可以看出。表层的低阻体( 2 0 时，取工( ，k ，z ) = 2 o ；当 ，( ，k ，f ) 一一卸一 一柏一帅一柏柏一 - - 、_ _ _ _ _ _ _ - _ - _ 、一，_ _ _ _ h 一 摄 距 ( m ) 槎 距 ( 呻 d b 图3 14 模型2 加5 随机噪声的视电阻率断面图 a 、b 、c 、d 、e 分别为第l 、2 、3 、4 、5 奈测残的视电阻率断面 桂林工学院硕士学位论文 图3 1 5 模型2 加5 随机噪声反演后的电阻率断面图 3 、b 、c 、d 分别为t 、2 、3 、4 、5 条刹线问对应曲电阻率断面e h - 1 5 米住置的水平切片 桂林工学院硕士学位论文 援 距1 0 ( m ) 位置( m ) j 幻曲一一 柏一柏一柏一 二一：孑二一 d 图3 16 模型2 加5 随机噪声的视电阻率断面图 a 、b 、c 、d 、e 分别为第1 、2 3 ，4 、5 条测线的视电阻率断面 桂林工学院硕士学位论文 图3