（地球探测与信息技术专业论文）频率域可控源二维有限元正演数值模拟.pdf

摘要 频率域可控源大地电磁法由于人工源的加入使其克服了大地电 磁法中天然场源信号微弱的缺点，因此在矿产普查、油气勘探等方面 得到了广泛的应用。本文主要对频率域可控源电磁法二维正演问题进 行研究，并通过全区视电阻率对过渡区以及近区数据进行校正。 本文采用有限单元法进行数值模拟。用m a t l a b 为程序编译工具， 采用双线性插值和双二次插值的有限单元法，推导出了相对应的边值 问题和变分问题满足的微分方程。为了模拟无穷远边界并满足计算机 的内存需求，在保证计算精度的情况下设计了非均匀网格剖分。在程 序编制过程中，只存储有限元系数矩阵的非零元素，大大节省了存储 空间。通过计算几个典型模型的卡尼亚视电阻率拟断面图和相位拟断 面图，得出了不同模型下电磁场响应特征的分布规律，验证了算法的 正确性。 针对频率域可控源电磁法中卡尼亚电阻率在过渡区和近区畸变 的问题，给出了全区视电阻率的迭代公式，并对典型的一维层状模型 以及简单二维模型进行了计算。过渡区和近区数据经过校正后，可以 正确反映出模型的地电特征，证明了线源下近区勘探的可能性。 关键词：频率域可控源电磁法，有限单元法，正演，全区视电阻率 a b s t r a c t b e c a u s eo ft h ea r t i f i c i a ls o u r c e ，f c s e mh a so v e r c o m et h ew e a k n e s s o fn a t u r a ls o u r c e ，w h i c hi st h es h o r t c o m i n go fm t , s o i th a sb e e nw i d e l y u s e di nt h em i n e r a l ，o i la n dg a se x p l o r a t i o na n ds oo n t h i sp a p e rh a s d o n es o m es t u d i e so nt h e2 df o r w a r dp r o b l e mo ff c s e m ，a n du s e st h e w h o l er a n g ea p p a r e n tr e s i s t i v i t y n e a rt h es o u r c e t or e c t if yt h ed a t ew h i c ha r ec o ll e c t e d t h ef i n i t ee l e m e n tm e t h o dw a su s e df o rn u m e r i c a ls i m u l a t i o n b a s e d o nm a t l a ba sp r o g r a m m i n gt o o l ，w ea d o p t e dt h ef i n i t ee l e m e n tm e t h o do f l i n e a ri n t e r p o l a t i o na n ds e c o n di n t e r p o l a t i o ni nar e c t a n g u l a re l e m e n tf o r s o l v i n gt h ef c s e mf o r w a r dp r o b l e m a n dd e r i v e dt h ec o r r e s p o n d i n g f o r m u l a sf o rc a l c u l a t i n g i no r d e rt os i m u l a t et h ei n f i n i t yb o r d e ra n dm e e t t h ed e m a n do fc o m p u t e rm e m o r y , t h en o n - u n i f o r mg i r di sd e s i g n e dt o e n s u r i n gt h ea c c u r a c y i np r o g r a m m i n g ，w eo n l y s t o r e dt h en o n z e r o e l e m e n t so ft h ef i n i t ee l e m e n tm a t r i xf o rs a v i n gt h ec a l c u l a t i o nm e m o r y b vn u m e r i c a ls i m u l a t i o no fs e v e r a lt y p i c a lm o d e l s ，w ev e r i f i e d t h e c o r r e c t n e s so ft h ef o r w a r dm o d e l i n ga l g o r i t h m t h ei t e r a t i v ef o r m u l ao ft h ew h o l er a n g ea p p a r e n tr e s i s t i v i t yw a s d e v e l p p e dt os o l v i n gt h ep r o b l e mo fd i s t o r t i o n o fc a g n i a r da p p a r e n t r e s i s t i v i t yw h e n t h ed a t ea r ec o l l e c t e dn e a rt h es o u r c e t y p i c a l1dl a y e r e d m o d e la n ds i m p l e2 dm o d e la r ec a l c u l a t e d a f t e rr e c t i f yt h ed a t e ，i tc a n a c c u r a t e l yr e f l e c tt h ec h a r a c t e r i s t i c so f t h em o d e lw h i c hv e r i f i e dt h ef a c t t h a tt h ee x p l o r a t i o nn e a rt h es o u r c ei sp o s s i b l e k e y w o r d ：f r e q u e n c yd o m a i nc o n t r o l l e d - s o u r c ee l e c t r o m a g n e t i c ，f i n i t e e l e m e n tm e t h o d ，f o r w a r dp r o b l e m ，t h ew h o l er a n g ea p p a r e n tr e s i s t i v i t i l 中南人学硕十学位论文 第一章绪论 1 1 引言 第一章绪论 频率域町控源大地电磁( 测深) 法“f r e q u e n c yd o m a i nc o n t r o l l e d s o u r c e e l e c t r o m a g n e t i c ”最早是加拿大多伦多大学的d w s t r a n g w a y 教授和他的研究 生m y r o ng o l d t e i n 于1 9 7 1 年提出的。针对大地电磁法场源的随机性和信号微弱， 以致观测十分困难这一状况，他们提出了一种改变方案采用可以控制的人工 场源。他们的论文于1 9 7 5 年公开发表，从理论和实践两方面奠定了可控源大地 电磁的基础。可控源电磁法得到实际应用，一些公司相继生产用于可控源电磁法 测量的仪器和解释用的软件。特别自8 0 年代以来，方法理论和仪器都得到了很 大的发展，应用领域也扩展到普查、勘探石油、天然气、地热、金属矿产、水文、 环境等各个方面，从而成为受人重视的一种地球物理方法i 。 可控源电磁法和大地电磁法相比，最大的区别就是加入了人工场源，这可以 使得电磁波的信号大大的增强，弥补了大地电磁法电磁波信号微弱的缺点。在实 际工作中，可控源电磁法具有工作效率高，勘探范围大，水平以及垂直分辨率较 高，不受高阻层屏蔽等优点。 1 2 频率域可控源正演研究现状 随着计算机技术的迅猛发展，数值模拟技术在科研中变得越来越重要，在大 型工程建设、资源勘探等关系国计民生的重要领域，数值模拟技术已成为最为关 键环节之一。地球物理勘探领域的数值模拟技术主要包括两个方面的内容：正演 和反演。在实际情况下，反演是核心。在今天，资料处理可以说已经等同于反演。 然正演是反演的基础，反演所取得的重大进展往往来自于j 下演方面的突破。 从j 下演方面来讲，主要的数值模拟方法主要分为两类：微分方程法和积分方 程法。微分方程法的一个主要优点在于其适应性强。给定一个求解区域，无论模 型复杂与否，都以统一的方法给出求解结果，这一点对于实际问题非常重要。而 相应的积分方程法对于复杂问题无能为力。因此实践中通用的还是微分方程法。 本文将主要对微分方程法中的有限元法作讨论。 有限单元法最早是在力学中得到应用。最早将有限元法应用于电磁法的文章 中南大学硕士学位论文 第一章绪论 是c o g g o n l 2 1 于1 9 71 年发表的( ( e l e c t r o m a g n e t i ca n de l e c t r i c a lm o d e l i n gb yt h e f i n i t e e l e m e n t m e t h o d ) ) ，在这篇文章中他用有限单元法对二维模型进行了数值模 拟。在此之后，有限单元法在电磁法数值模拟中逐渐普及丌来，s i l v e s t e r 也将有 限元应用于大地电磁测深中。w i l l i a m l 3 1 ( 1 9 7 6 ) 在用有限单元法进行数值模拟的时 候，将双线性插值和双二次插值进行了对比，并得出，高次插值数值模拟精度会 有显著提高的结论。硒i o r 4 1 ( 1 9 7 7 ) 在其相关论文中提出了通用性网格的理论，这 使得在正演计算中，计算的速度和精度都得到了很大的提高。w a n n a m a k e r i s l ( 1 9 8 6 ) 通过三角剖分，模拟大地电磁中的遇到的地形问题，并且提出了使用的地形校正 方法。d el u g a o 等【6 1 ( 1 9 9 6 ) ；币i f r a n k 等1 7 1 ( 2 0 0 4 ) 也通过不同形式的有限元研究y - - 维大地电磁的正演问题。m o g i l 8 1 ( 1 9 9 6 ) 、z y s e m a n 等1 9 1 ( 2 0 0 0 ) 、s o n 等l i o j ( 2 0 0 3 ) 和 m i t s u h a t a 等1 ! ( 2 0 0 4 ) 在三维研究中取得成功。l i ( 2 0 0 6 ) 对海底可控源二维模型 进行了相关的研究。 在国内，有限单元法最早是由冯康提出来。随着8 0 年代初期，李大潜的专 著有限元素法在电法测井中的应用的发表，有限单元法被引入到电磁法中。 自8 0 年代中期以来，在国内有限元法的应用得到迅速推广，范围涉及电磁勘探 的所有领域。一大批优秀的论文和专著涌现出来。在m t 方面，有限元二维j 下演 问题一再作为硕士论文被提出。在可控源勘探领域，有限元正演问题屡屡被列为 国家自然教育基金项目。胡建德等1 1 2 】( 1 9 8 2 ) 将有限元法和大地电磁二维正演 相结合。此后，徐世浙等f 1 3 - 1 6 1 ( 1 9 8 3 ，1 9 9 4 ，1 9 9 5 ，1 9 9 6 ) 对有限元法进行了深 入的研究，发表了一系列有限单元法在电磁法中的应用的论文。陈小斌l l7 j ( 2 0 0 0 ) 提出了有限元直接迭代法，可大大提高运算效率。阮百尧等f 1 8 1 【2 4 1 ( 2 0 0 0 、2 0 0 1 ) 对三维情况下的地电响应进行了数值模拟，并对带地形的复杂三维模型做了相关 的研究。王若等【2 5 1 ( 2 0 0 7 ) 在c s a m t 三维正演中采用有限单元法进行数值模拟， 取得较好结果。其他的一些学者【2 6 l - 1 3 2 1 也做了相关的研究，推动了有限单元法在 电磁法中的应用。 有限单元法和其他的数值方法相比较，主要具有以下几个优点： ( 1 ) 有限单元法剖分的网格具有很大的灵活性，可以根据不同的地电情况 采取不同的剖分手段，以适应不同的地质构造。 ( 2 ) 从有限单元法得出的方程组是稀疏对称的矩阵，可以使得求解大大简 化。 ( 3 ) 在有限单元法中，边界条件容易处理，便于在程序中实现。一 所以，现在有限单元法已成为电磁勘探中一种主要而通用的数值模拟方法。 大多数反演程序是选择以有限单元法作正演。这充分体现了有限元法的强大的生 2 中南大学硕士学位论文第一章绪论 命力。 1 3 全区视电阻率校正研究现状 目前，在f c s e m 中，一般计算视电阻率采用卡尼亚视电阻率的计算公式， 该公式的优点是计算简单、方便，在远区能够j 下确客观的反映出地电断面的垂线 变化，但是卡尼亚视电阻率在过渡区以及近区计算视电阻率是会出现明显的偏 差，影响我们对资料的解释。近区视电阻率也存在类似的问题。在实际野外勘探 中，我们一般尽量把测区放在远区，但是由于地形条件或发射功率等原因，我们 在实际勘测中有时候不得不在过渡区甚至于近区进行工作。如果对于近区或过渡 区的数据全部放弃，会造成不必要的浪费。因此，更广泛地研究对近区以及过渡 区数据的校正方法是必要的，并具对实际的生产有很大的现实意义1 3 3 l 。 曹昌祺幽l 曾经在水平层状大地上的交流视电阻率一文中对波区视电阻 率的概念提出了自己的改进意见，s p i e s l 3 5 1 和w i l t l 3 6 l 等也分别提出了求取视电阻 率的方法，殷长春和朴化荣f 37 】提出关于电偶源的波区视电阻率的求取，汤井田 和何继善【3 8 1 通过分析水平电偶源提出了全区视电阻率的概念，并得到了很好的 效果。 这些等效的电阻率统一了远区、过渡区和近区的视电阻率的定义，可以较好 的反映出地下电性特征。 1 4 研究的目的与意义 虽然可控源大地电磁法是一种成熟的地球物理方法，但该种方法自身的特点 决定了，只有在波区，其电、磁场水平分量才随距离r 的变化具有相同的变化规 律，各种视电阻率才能客观地反应地电断面垂向变化，而在实际工作条件以及地 质条件的制约下，我们又不得不进入中间区甚至近区观测，因此定义全区视电阻 率可改善可控源大地电磁测深的应用效果。 正演计算是我们认识各种条件下可控源大地电磁响应特征的良好途径，通过 对不同模型的正演研究，可总结出典型条件下电磁场的分布规律。可通过正演， 在不同条件下，对比卡尼亚视电阻率和全区视电阻率的效果，在非波区全区视电 阻率是否可取得比较好的结果。而这些由正演得到的结论，对于资料分析处理、 反演、解释具有很重要的作用。 在目前几种主要的数值模拟方法中，本文选择了有限单元法，从变分原理出 3 中南人学硕十学位论文第一章绪论 发，将泛函的极值问题转化为多元的线性方程组，这使得分析的过程大大的简化； 同时，在求解微分方程的变分形式的时候，内边界条件可以相互抵消，这样就可 以不用去考虑复杂的内边界条件，简化了有限单元法的运算。 目前在实际地球物理工作中所遇到的问题可近似看为二维问题的占很大的 比例。对于频率域可控源大地电磁二维问题的正演，本人欲采用有限元法进行数 值模拟，并对结果进行分析。 1 5 本文的主要工作及结构 本文的主要工作由以下方面组成：首先对二维可控源大地电磁的变分问题满 足的微分方程进行推导，用矩形单元对研究区域进行剖分，每个剖分单元分别采 用双线性插值和双二次插值进行计算。把求解变分问题转化为求解一组多元线性 代数方程。并同时计算全区视电阻率，编出实用的正演程序，通过计算电磁场的 分布，研究不同模型电磁场的响应特征以及分布规律。 论文共分五章，各章内容如下： 第一章绪论：介绍频率域可控源大地电磁法有限元数值模拟以及非波区视 电阻率校正国内外发展现状，介绍论文研究的目的以及意义，并给出本文的主要 工作及结构。 第二章二维频率域可控源大地电磁正演模拟基本原理：从基本的电磁场公 式出发，推导二维频率域可控源的边值问题以及变分问题的方程。采用矩形单元 对所研究的区域进行剖分，在剖分的时候采用双线性插值和双二次插值两种方 法，然后求解出数值解。并对卡尼亚视电阻率及其相位和全区视电阻率的公式进 行了推导。 第三章i f 演模型的计算：对均匀半空间，典型一维层状模型以及典型二维 模型进行了计算，并对计算出来的结果进行了分析。 第四章全区视电阻率在近区数据中的运用：建立典型的一维和二维模型， 对近区数据用全区视电阻率公式进行计算，并和m t 的结果相比较，并对对比结 果进行分析。 第五章结论与建议：通过第二章的理论分析以及第三章和第四章典型模型 的计算，得出结论并给出建议。 4 中南大学硕士学位论文第二章二维频率域可控源大地电磁正演模拟基本原理 第二章二维频率域可控源大地电磁正演模拟基本原理 对于线源频率域可控源电磁响应正演问题早在2 0 世纪7 0 年代就有研究， h o h m a n n 将积分方程退化为矩阵方程得到了线源作用下异常体的电磁场数值解妇钔 下面讨论有限元方程的推导过程。 2 1 边值问题 设电磁场的时问因子为e 却，则其麦克斯韦方程组可写为1 3 9 4 2 1 ： v x e = f 掣日 ( 2 - 1 ) v x h = ( 仃一泐占) + 以 ( 2 2 ) 在式中，e 是电场矢量；h 为磁场矢量；缈为角频率；为磁导率；仃为电 导率；g 为介电常数；以为外加的源矢量。 假设地下的电性结构是二维的，将原点取在地面，其走向方向假定为沿x 轴，垂直向下的方向为z 轴，位于地面同时和x 轴相垂直的是y 轴。外加源的方 向是沿着x 轴的，此时存在不随走向改变而改变的x 方向的电场。 将2 2 式展开，则可得到： 冬一孕却一泐) t + 厶 ( 2 - 3 ) 口y 宓 将2 1 式展开，并只考虑x 方向的电场，可以得到： 誓= 卸以 ( 2 _ 4 ) 一孕：f 掣h ： ( 2 - 5 ) 将2 4 和2 5 代入2 3 中，可以得到： 晏( 士誓) + 导( 土粤) + ( 盯一泐) 巨：也 ( 2 6 ) c vz 缈口c v d zz 反堆lo z 令p = ；q = 口一f 国s ；材= 巨；g 三一厶。则式2 - 6 可以简化为： 5 中南人学硕十学位论文 第二章二维频率域可控源人地电磁正演模拟基本原理 晏( p 豢) + 昙( p 害) + 掣：g 瓦( p 瓦) + 瓦( p 瓦) + 掣2 9 2 1 1 外边界条件 取外边界足够大，如图2 1 所示， 处为零。设图中“ 处为源的位置， 质和不均匀体两种不一样的介质。 ( 2 7 ) 可以使得局部不均匀体的异常场在外边界 圆圈为局部不均匀体，1 和2 表示周围介 ； y a d i o 图2 1 研究区域示意图 r i 为不均匀体的边界，在空中边界d a b c 上，电磁波传播的方程为： 甜= u o e 一妒 ( 2 8 ) 式中，k o 为空气中的波数，k o = 缈2 鳓；勿为圆频率；气和鳓分别为空 气中的介电常数和磁导率；r 为源到任意一点的距离，对2 8 式求导可得： 譬= 一p 却 ( 2 9 ) u r 对2 8 式两边分别乘以，并与2 - 9 式相加可得： _ 0 u + k 0 = 0 ( 2 1 0 ) 6 中南大学硕士学位论文 第二章二维频率域可控源大地电磁正演模拟基本原理 因为宇：c o s 0 ，其中，n 为外边界的外法线方向，要表示沿n 方向的方向 a ” u n 导数，秒为源到边界上任意一点的矢径与边界法线方向的夹角，代入上式，可以 得到空气中外边界上的边界条件： 坐+ 厶c o s 锄：0 ( 2 1 1 ) 锄 。 在地层中的边界e d 和f c 、e f 上，电磁波的传播也应满足传播方程，同样 可以推导出类似于2 1 l 式的边界条件，可以写成： 坐+ k c o s 0 ”：0 ( 2 一1 2 ) 锄 式中，k 为地层中的波数，后= 缈2 c f l + i m p o - ；彩为圆频率；占和分别为 地层中的介电常数和磁导率，对于电磁法来说，由于频率足够低，可以忽略位移 电流的影响，地层中的波数可以写成k = 咖仃。 2 1 2 内边界条件 在介质的分界面上电场的切向分量连续，即 i 2 u 2 磁场的切向分量连续 q = 以c o s 口+ 皿s i n a = 面1t i a ec o s 口- - s m ) = 壶鲁 ( 2 - 1 3 ) ( 2 1 4 ) 式中，口为内边界的切向分量与y 轴的夹角，所以有： 玟盥：p 堕 ( 2 一1 5 ) a 22i 2 小1 ) ， u n u n 综上所述，可以将二维电磁场的边值问题联合写为： v ( p v u ) + q u = g q 学+ c o s o u = 0 d a b c u r l 坐 4 - k c o s o u ：0d e f c 。吃 a 抛， a 2 仍苫 o n o n f l f l 7 中南大学硕士学位论文第二二章二维频率域可控源大地电磁正演模拟基本原理 2 2 变分问题 对于二维电磁场的边值问题，我们构造如下泛函1 3 9 4 2 1 ： 成： l ( u ) = 【主p ( v 甜) 2 一芝1 2 + g u 】d q q 。 。 其变分为 万，( 甜) = p v z ，v 6 u d q 一如甜一g ) 6 u d q qq 应用恒等式 v ( 彳缈) = v 彳伊+ a v 缈 ( 2 1 6 ) ( 2 1 7 ) ( 2 1 8 ) 有 甜( z f ) = f v ( p v z ，万材) d q f v ( p v 甜) + 鼋l l - - g d u d f 2 ( 2 一1 9 ) qq 由前文易知，上式第二项为零。 再由散度定珲，v ( 彳缈矽q = 叮4 彬1 1 ，上式可以写成： 删2 曼局挚印r + 叮np z o u i 2 掣r ( 2 2 0 ) 式中，r 是区域的外边界；f 。是内边界，根据内边界条件有 a _ o n i 她：仍磐6 u 2 ( 2 2 1 ) o n , o n 、 但是外法线向量和n 2 方向相反，所以在内边界上相互抵消，2 2 0 可以写 甜( 炉即- - 鲁s u l d fr 将外边界条件代入上式，有： 万，( 甜) = 叮一p c 。s o u 6 u d f + 叮一p k c 。s o u 6 u d f d a b cd e f c ( 2 2 2 ) = 万叮一争屯c 。s 锄2 订曼争后c 。s 砌2 d r c 2 彩) d a b c。defc。 变形为： 8 中南大学硕士学位论文第二章二维频率域可控源人地电磁止演模拟基本原理 o q c 甜，+ 童，争毛c 。s 砌2 订+ 万曩扣c o s 砌2 订，= 。( 2 - 2 4 ) 所以边值问题和下列变分问题等价 f ) = 且吾p ( v u ) 2 一了1g 2 + g u a q + q 尹1 k c 。s 砌2 d r + 万叮吾肚c 。s p 甜2 j r q 。 。dabfj-debr 万f ( ”) = 0 2 3 有限单元法 2 3 1 有限单元法插值理论 ( 2 2 5 ) 2 3 1 1 双线性插值 取图2 2 所示的正方形单元( 母单元) ，四个顶点的编号及其坐标，如图 上所示，构造如下的形函数1 1 4 1 【删： l = 三( 1 一似l + 矾2 = 丢( 1 一o ( 1 - u ) 3 = 丢( 1 + 孝) ( 1 - r ) ，4 = 丢( 1 + 纵1 + 7 7 ) 或统一写成： m 2 言( 1 + 参踟+ 玩巧) 其中毒，r , 是点i ( i = 1 ，e eee ee p4 ) 的坐标。 形函数2 2 6 满足 删，= 髋二 的要求，其中i 代表点号。 将形函数展开： n l = n 蠢q + a 玎+ a 3 9 r l + a i 其中q ，a 4 是常数，m 含有勿，孝，7 和常数项。 9 ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) 中南大学硕1 ：学位论文第二章二维频率域可控源人地电磁正演模拟基本原理 ( 一l ，1 ) ( 1 ，1 ) 图2 2 双线性插值的母单元 ( 1 ，1 ) 图2 3 上任意四边形单元( 子单元) 的四个顶点的坐标分别是( 而，m ) ， ( _ ，欺) 函数值是，子单元上的u ，x 和y 可表示成： 将形函数代入上式，整理后，可得： f 甜= 4 勃+ 4 孝+ 4 刁+ 4 x = 蜀勿+ 岛f + 最7 + 反 【y = c l 勃+ c 2 孝+ c 3 刁+ c 4 其中4 1 ，c 4 是常数。 1 0 ( 2 3 0 ) ( 2 3 1 ) 拱 m m + + + 吩黾憾 + + + 屹恐影 + + + 薯m m m m i i l i 一一 甜 x y 中南大学硕士学位论文第二章二维频率域可控源人地电磁正演模拟基本原理 y x 图2 3 双线性插值的子单元 令r l = c ，这在勃平面上是一条水平线，代入2 3 l 式，得： u = a + a 2 ，x = 6 l f + 也，y = q # + c 2 ( 2 3 2 ) 其中q ，c z 是常数，所以u ，x 和y 均是善的线性函数，当f 从一l 变 到l 时，u 是线性变化的，x y 平面上的点也将沿着一条直线变化。同样，令孝= c 时，u ，x 和y 均是r 的线性函数，所以2 2 6 是双线性函数，2 3 0 是双线性插值， 与母单元的四条边对应的子单元的四条边都是直线。 2 3 1 2 双二次插值 如图2 4 所示的正方形的母单元，八个节点的编号以及坐标如图所示，构造 如下的形函数： 5 = 6 = n 7 = 8 = 易得，形函数2 - 3 3 满足 删，= 骺二 的要求，将形函数展开可得 m = q f 2 r l + a 2 孝r 1 2 + 吩孝2 + a 4 r 1 2 + 吩勿+ 吼f + 伤矽+ 口8 上式含八个系数。 图2 4 双二次插值的母单元 1 2 ( 2 3 3 ) ( 2 3 4 ) ( 2 3 5 ) d d 一 一 穆 矽 + 一 手 善 1 一 ，l，l 棚 功 + 一i ，王 乒 孝 一 + o o l一4一4 = = 心 m 枞 枞 十 一 l l x f 一 + o o ，一4一4 = = m m 中南大学硕十学位论文第二章二维频率域可控源人地电磁止滨模拟基本原理 图2 9 矩形单元剖分、双线性插值的子单元 2 3 2 3 单兀分析 在实际做单元分析的时候，对于外边界条件，我们采用无穷远边界条件，也 就可认为电场在外边界已经衰减到非常小，接近于零，对结果的影响不大，可以 忽略。 将变分问题重写如下： f ( 甜) = 踟( v 甜) 2 d q 一踟z ，2 饱+ i g u d q ( 2 - 4 4 ) q-n-q 将2 4 3 代入上式，整理后，对于单元e ，第一项积分可以写成单元积分的形 式： 多( v ) 2 d q = 圭嘭【彬】u = 1 2 l r ，e 7 k k 虬 ( 2 - 4 5 ) 其中， 膨= 肛( 等等+ 警警舭= “p t 虿a n 虿a n 7 4 + 型a 7 叭a r b 4 2 ) 舢考r 却 1 6 ( 2 - 4 6 ) 中南人学硕士学位论文第二章二维频率域可控源人地电磁正演模拟基本原理 = l 巍2 a + 2 f l 2 a + 2 f l a - 2 f l2 a j k k ：l 口一2 2 口+ 2 l l 一口一 一2 口+ i l 2 口+ 一口+ 2 ， 其中， pb a = 1 - 6 口 侈：旦竺 6b 第二项的单元积分 岿砌q = 兰畦咛】配= 互1 嘭饥 其中， e f e q n , n a y 出戈沁n j 百a b d 考如 p = 协4 其中， a b 口。磊g j o 第三项的单元积分 p d q 2 广坝y m 耐q = 一7 尸 其中， p = 【三，o oo 】 ( 2 4 7 ) ( 2 4 8 ) ( 2 - 4 9 ) ( 2 5 0 ) ( 2 5 1 ) ( 2 5 2 ) 将k k 、k 殄、p 扩展成全体节点的矩阵( 扩展矩阵上面加“一表示) ，然 后将其代入2 _ 4 4 中可以得到： f ( 甜) = 丢u7 夏k u 一言u r i 知u 一三u 7 尸 对上式取变分可得： 5 f ( 材) ：_ 1 8 u i h u p ：0 1 7 ( 2 5 3 ) ( 2 5 4 ) ! 塑叁竺堡：兰堂堡笙苎 。第二章二二维频率域可控源人地电磁正演模拟基本原理 二二二二= ：= ：= ：三：! = ：= 可以得到线性方程组： k u = p ( 2 5 5 ) 其中，k = k - k 知。解这个方程组，可以得到各个节点的电场值。 2 3 3 矩形单兀的双二次插值 2 3 3 1 区域剖分 用矩形单元对区域进行剖分，如图2 1 0 所示，每个单元的四个角点和四边 形的中点，共八个点作为节点。并给每个单元的电阻率赋值，矩形单元的宽度为 a ，高度为b ，单元的八个节点的编号如图2 1 1 所示，在单元内进行双二次插值。 l ，2 ，代表节点 ，代表单元 图2 1 0 双二次插值的节点和单元的排列次序 图2 11 双二次插值单元中的节点号次序 2 3 3 2 双= 次插值 图2 1 2 是母单元，图2 1 3 是子单元。两个单元之间的坐标变换关系为： 1 8 中南大学硕士学位论文第二章二维频率域可控源人地电磁上 演模拟基本原理 为： x ：h + 竺f ，y ：+ 一b ，7 ( 2 5 6 ) x 。h + 互亏，y 2 + 芝矽 ( 2 。 其中，是子单元中点的坐标，a ，b 是子单元的两个边长，微分关系 出= 詈嘶，咖= 兰却，蚴= 警蜊，7 双二次插值的形函数见2 - 3 3 ，单元中n 的插值函数为： = m 其中m ( i = l ，8 ) 是单元中8 个节点的待定u 值。 5 2 图2 1 2 矩形单元剖分双二次插值的母单元 1 9 ( 2 5 7 ) ( 2 5 8 ) 中南大学硕士学位论文第二章二维频率域可控源大地电磁正演模拟基本原理 x 图2 1 3 矩形单元剖分双二次插值的子单元 2 3 3 3 早兀分析 对2 4 4 第一项积分： 咖( v 妒d q = 1 u t k l e u e , = j 1 蟛眈 其中， 磅= p t c 等芸，c 考参+ c 筹多，c 等挚等d 跏 具体计笪公式如下： ( 2 5 9 ) ( 2 6 0 ) 中南人学硕士学位论文第二章二维频率域可控源大地电磁j i ：演模拟基本原理 f 毛l = 5 2 a + 5 2 , 3 ，k 2 l = 1 7 a + 2 8 f l l k 3 l = 2 3 a + 2 3 f l ，缸i = 2 8 a + 1 7 f l l 屯l = 6 a s o p ，丸i = _ 4 呶一6 i k 7 i = - - 6 a s o p ，毛l = 一8 0 a + 6 f l ik = 毛。，k = 缸 i 屯：= 七3 ，k = k 5 。 l k 6 2 = 缸，= 岛。 lk = 晚。，k = 毛。 i k 4 3 = k 2 l ，k s 3 = k 7 i l k 6 ，= 毛。，岛，= 毛。 l k 8 3 = 屯。，k = & 1 k = 岛。，k = k 6 。 l 岛。= 屯。，k = 。 l 屯5 = 4 8 a + 1 6 0 5 ，k = 0 i k 7 5 = - - 4 8 a + s o p ，k = 0 lk = 1 6 0 a + 4 8 3 ，k 7 6 = 0 i k = 8 0 a r 一4 8 1 ，k 7 7 = k 【，= o 毛。= k 口：旦竺，：旦旦 对于第二项积分： 多砌q = 扣巧。 w e = 三2 卿眈 其中砖。的具体计算公式如下： k 知：a 6 26 对 32 6 232 6 6 6 8 83 2 称 一8 6 6 82 03 2 8 8 6 61 62 03 2 6 8 8 62 01 62 0 3 2 其中口= 蔫g 2 l ( 2 6 1 ) ( 2 6 2 ) ( 2 6 3 ) 中南人学硕士学位论文 第二章二维频率域可控源大地电磁正演模拟基本原理 了。 总体系数矩阵的集成，解线性方程组等和双线性插值一样，这里就不在重复 2 4 计算视电阻率和相位 2 4 1 卡尼亚视电阻率和相位的计算 在得到各个节点的电场之后，可以求出地面上与电场垂直的磁场： 日 ，_ 上堡 ( 2 6 4 ) j t po z 其中电场的偏导数计算可以根据文献1 4 1 4 3 计算 警l 瑚= 去( _ 1 1 e x ，+ 1 8 e ，2 - 9 耻2 ( 2 6 5 ) 其中，e x 。为地面上节点的电场值，巨：、巨，、巨。分别为地面下第一、二、 三个节点的电场值。 当地面上相互垂直的电场与磁场已知以后，可以通过下面的式子计算视电阻 率： p = 万i 瓦g 1 4 掣，。 令z l 专 相位可以通过下面的式子计算： 疙= a r c t a n 器 2 4 2 全区视电阻率的计算 ( 2 6 6 ) ( 2 6 7 ) 对于均匀半空间，供电电源为线源，在地面，地面的观测场为【4 5 1 ： e = 丽i 巧# a i 驴【f 蟛t ( i k o y ) 一f 锥( 匆) 】( 2 - 6 8 ) 其中k ( z ) 为第二类一阶变形的贝塞尔函数，k 为波数，y 为偏移距，当y 较大的时候，定义远区视电阻率为： 中南人学硕十学位论文第一二章二维频率域可控源人地电磁止演模拟基本原理 一万- 1 ， 拯i x l l p 2 硒丽面 ( 2 6 9 ) 对于远区视电阻率，只有当y 足够大的时候，计算出来的远区视电阻率才和 实际的电阻率相近。为了解决这种情况，将2 6 8 和2 6 9 稍作变形定义全区视电 阻率如下： 砰：碍一掣 i k o k , ( 一曦工) 一i k i k t ( 一i k , x ) 】 , r t l ，x ( 2 7 0 ) p ：f 牟 ( 2 7 1 ) 胪陶 吃川 先预先设定一个初始的电阻率值( 任意给出) ，计算出毛，将这个毛代入到 2 - 7 0 中求出新的毛，再将新的毛代入到2 7 1 中求出新的电阻率，当误差满足一 个值的时候，迭代停止，所得到的电阻率值就是要求的全区视电阻率。 中南人学硕十学位论文 第二章正演模拟的计算 第三章正演模型的计算 首先为了验证数值解法的正确性，对比均匀半空间解析解和数值解的电场 值，如图3 1 所示。 图3 1 均匀半空间解析解和数值解的电场值比较 ( 背景电阻率为1q m ，频率为i o h z ) 从图3 1 可以看出，数值解和解析解的结果十分接近，证明了数值解的正确 性。下面对具体模型进行计算。 模型1 图3 2 模型l 网格剖分单元数目：x 轴( 水平) 方向为5 5 个，y 轴( 垂直) 方向为4 5 个( 其中地下3 5 个，空气中为1 0 个) ，源位于水平方向第1 4 个点处。 x 轴网格宽度( 单位：m ) 3 0 0 0 15 0 0 5 0 0 2 0 0 10 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 3 0 2 0 3 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 ，4 中南人学硕+ 学位论文第二章正演模拟的计算 l0 0 1o o 1o o 2 0 0 3 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 10 0 0 10 0 0 2 0 0 0 2 0 0 0 2 0 0 0 2 0 0 0 3 0 0 0 3 0 0 0 3 0 0 0 3 0 0 0 3 0 0 0 3 0 0 0 3 0 0 0 y 轴网格宽度( 单位：m ) 1 0 1 0 1 0 1 0 1 0 2 0 3 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 10 0 2 0 0 4 0 0 8 0 0 16 0 0 3 2 0 0 6 4 0 0 10 0 0 0 2 0 0 0 0 3 0 0 0 0 10 3 0 10 0 3 0 0 10 0 0 3 0 0 0 10 0 0 0 3 0 0 0 0 10 0 0 0 0 3 0 0 0 0 0 、 表3 1 模型l 双线性双= 坎电阻率与相位值对比 ( 偏移距j6 0 0 0 r a ) 频率 双线性电阻率双二次电阻牢双线性相位双二次相位 ( t t z ) ( q m )( q m ) ( 。)( 。) 0 1 2 54 7 1 1 0 2 2 6 6 24 8 8 5 6 1 9 6 4 93 6 4 4 9 1 6 7 4 13 6 3 8 4 2 7 3 6 0 2 5 3 8 7 8 4 9 2 7 1 3 4 0 2 6 6 9 3 2 2 33 3 4 9 0 7 4 4 2 43 3 5 6 9 3 4 5 8 o 531 2 8 3 0 9 3 2 43 2 6 5 6 0 3 3 5 83 2 1 8 4 7 3 3 3 63 2 0 7 7 3 8 3 12 5 3 9 0 5 9 5 4 22 6 3 7 8 5 6 1 3 43 1 1 9 2 0 3 0 3 73 0 8 7 9 7 3 9 6 22 0 1 1 4 8 4 7 92 0 6 91 4 9 8 6 33 0 0 7 3 3 7 8 8 52 9 6 3 4 3 8 4 4 415 0 9 4 7 8 9 7 81 5 4 6 9 5 2 1 7 12 9 8 4 9 2 2 3 8 92 9 4 4 7 8 5 8 3 811 3 2 2 7 1 8 6 61 1 4 6 4 1 0 4 73 2 8 3 2 3 7 8 73 2 2 7 3 0 5 1 2 1 69 5 0 5 4 0 8 2 8 79 5 4 9 2 9 2 8 2 93 9 0 1 2 7 0 5 3 83 8 5 5 1 8 0 2 9 3 29 8 4 8 1 7 9 0 4 39 7 7 9 7 6 0 6 3 84 3 4 6 l1 6 6 8 94 3 2 7 1 8 9 0 6 6 41 0 0 3 2 1 0 1 0 510 0 0 2 5 2 8 7 54 3 9 7 1 5 6 6 1 94 4 0 0 9 8 1 6 5 1 2 89 9 9 7 0 9 8 7 2 49 9 9 5 2 0 8 4 3 64 4 4 4 7 5 4 5 7 94 4 5 0 0 7 3 1 3 2 5 69 9 8 7 7 6 9 0 7 39 9 9 8 8 0 4 7 94 4 6 8 1 6 2 1 5 34 4 7 5 3 3 2 6 3 5 1 29 9 7 4 2 0 2 7 9 71 0 0 0 0 1 5 2 5 44 4 7 8 7 6 9 5 8 44 4 8 8 1 0 1 7 2 1 0 2 4 9 9 51 8 0 4 6 5 71 0 0 0 1 2 2 2 2 84 4 8 4 9 2 2 6 7 4 4 4 9 4 6 5 3 8 2 0 4 89 9 2 8 5 6 1 5 5 41 0 0 0 2 1 7 8 4 14 4 9 2 9 8 4 6 2 34 4 9 8 1 7 1 5 1 4 0 9 69 9 2 3 8 7 9 4 511 0 0 0 2 3 9 7 8 54 5 0 4 4 4 2 5 3 24 5 0 0 4 7 0 1 2 中南大学硕士学位论文第三章正演模拟的计算 图3 3 模型1 双线性、双二次电阻率、相位对比 ( 偏移距6 0 0 0 m ) 从图3 3 可以看出，在高频段的数值解和真实的地电参数吻合的非常好，随 着频率的减小，误差增大，特别是在过渡区以及近区，误差很大，很好的反映了 频率域可控源在近区、过渡区以及