高密度电法-专著

2 高密度电阻率法 高密度电阻率法仍然是以岩 土导电性的差异为基础 研究人工施加稳定 电流场的作用下地中传导电流分布规律的一种电探方法 因此 它的理论基础 与常规电阻率法相同 所不同的是方法技术 高密度电阻率法野外测量时只需将 全部电极 几十至上百根 置于观测剖面的各测点上 然后利用程控电极转换 装置和微机工程电测仪便可实现数据的快速和自动采集 当将测量结果送入微 机后 还可对数据进行处理并给出关于地电断面分布的各种图示结果 显然 高密度电阻率勘探技术的运用与发展 使电法勘探的智能化程度大大向前迈进 了一步 由于高密度电阻率法的上述特点 相对于常规电阻率法而言 它具有 以下特点 1 电极布设是一次完成的 这不仅减少了因电极设置而引起的故障和干扰 而且为野外数据的快速和自动测量奠定了基础 2 能有效地进行多种电极排列方式的扫描测量 因而可以获得较丰富的关 于地电断面结构特征的地质信息 3 野外数据采集实现了自动化或半自动化 不仅采集速度快 大约每一测点 需2 5s 而且避免了由于手工操作所出现的错误 4 可以对资料进行预处理并显示剖面曲线形态 脱机处理后还可自动绘制和打 印各种成果图件 5 与传统的电阻率法相比 成本低 效率高 信息丰富 解释方便 一 高密度电阻率法采集系统 早先的高密度电阻率法采集系统采用集中式电极转换方式 如图 4 1 所示 进行现场测量时 用多芯电缆将各个电极连接到程控式电极转换箱上 电极转 换箱是一种由微片机控制的电极自动转换装置 它可以根据需要自动进行电极 装置形式 极距及测点的转换 电极转换箱开关由电测仪控制 电信号由电极 转换箱送入电测仪 并将测量结果依次存入存储器 图 4 1 高密度电阻率法测量系统结构示意图 集中式 随着技术的发展 高密度电法仪日趋成熟 表现在 采用嵌入式工控机 大大提高系统的稳定性与可靠性 采用笔记本硬盘存储数据 可以满足野外长 时间施工的工作需求 系统采用视窗化 嵌入式实时控制与处理软件 便于野 外操作 可实现多种工作模式的转换 计算机与电测仪一体化 携带方便 新 一代高密度电法仪多采用分布式设计 所谓分布式是相对于集中式而言的 是 指将电极转换功能放在电极上 分布式智能电极器串联在多芯电缆上 地址随 机分配 在任何位置都可以测量 实现滚动测量和多道 长剖面的连续测量 图 4 1 高密度电阻率法测量系统结构示意图 分布式 系统可以做高密度电阻率测量 又可以同时做高密度极化率测量 应用范围宽 图 1 常用装置 高密度电阻率法在一条剖面上布置一系列电极时可组合出十多种装置 高 密度电阻率法的电极排列原则上可采用二极方式 即当依次对某一电极供电时 同时利用其余全部电极依次进行电位测量 然后将测量结果按需要转换成相应 的电极方式 但对于目前单通道电测仪来讲 这样测量所费时间较长 其次 当测量电极逐渐远离供电电极时 电位测量幅值变化较大 需要不断改变电源 不利于自动测量方式的实现 高密度电阻率法常用的装置见图 包括 温纳装置 Wenner Wenner Wenner 偶极 偶极装置 Dipole Dipole 三极装置 Pole Dipole Dipole Pole 温纳 斯伦贝谢装置 Wenner Schlumberger 等 图 a 温纳 装置 b 温纳 装置 c 温纳 装置 d 偶极 偶极装置 e 三极装置 f 温纳 斯伦贝谢装置 2 装置特点及视参数的计算 1 温纳装置 在高密度电阻率法中 由于温纳装置与异常对应关系好 是常用的装置之 一 最早的高密度电阻率法一般使用三电位电极系 所谓三电位电极系就是将 温纳装置 偶极装置和微分装置按一定方式组合后构成的一种测量系统 这是 由于电极转换需要时间 因此当连接好等距的 AMNB 四个电极后 可以作三次组 合 依次构成温纳装置 偶极装置和微分装置 或称为温纳 装置 温纳 装置和温纳 装置 这样在某一测点就可以获得三个电极排列的测量参数 温纳装置对电阻率的垂向变化比较敏感 一般用来探测水平目标体 温纳 装置的装置系数是 相比于其它装置而言是最小的 因而同样情况下 可a 2 观测到较强的信号 可以在地质噪声较大的地方使用 另一方面 由于它的装 置系数小 因此在同样电极布置情况下 它的探测深度也小 另外 温纳装置 的边界损失较大 温纳 装置 温纳 装置和温纳 装置三种排列形式 见图 视电阻率参数及计算公式为 I U k s ak 2 I U k s ak 6 I U k s ak 3 根据三种电极排列的电场分布 三者之间的视电阻率关系 sss 3 2 3 1 对高密度电阻率法而言 由于一条剖面地表电极总数是固定的 因此 当 极距扩大时 反映不同勘探深度的测点数将依次减少 图 显示了温纳 装置测点分布 图 温纳 装置测点分布示意图 最小电极距 n 间隔系数x 由图可见 剖面上的测点数随剖面号增加而减小 断面上测点呈倒梯形分 布 任意剖面上测点数可由下式确定 Dn Psum Pa 1 n 式中 n 间隔系数 Dn 剖面上测点数 Psum 实接电极数 Pa 装置电极数 对三电位电极系而言 Pa 4 对三极装置 Pa 3 如对温纳装置而言 设有 30 路电极 则 当 n 1 时 第一条3n30Dn 剖面上的测点数 令 可求出最大间隔系数为 总测点57D1 1 n D9 max n 数剖面数而言 总测点数为 N 9 1 330 n nN 2 偶极 偶极装置 偶极 偶极装置高灵敏度区域出现在发射偶极和接收偶极下方 这意味着本 装置对每对偶极下方电阻率变化的分辨能力是比较好的 同时 灵敏度等值线 几乎垂直的 因此偶极 偶极装置水平分辨率比较好 一般用来探测向下有一定 延伸的目标体 相对于温纳装置 偶极 偶极装置观测的信号要小一些 I U k MN s annnk 2 1 2 3 三极装置 三极装置有更高的灵敏度和分辨率 Sasaki 1992 同时 三极装置的两 个电位电极在网格内 因此受电噪声干扰也相对小一些 与偶极 偶极装置相比 三极装置所测信号要强一些 另外 三极装置可以进行 正向 单极 偶极 和 反向 偶极 单极 测量 因此边界损失小 I U k MN s annk 1 2 4 温纳 斯伦贝谢装置 温纳 施伦贝谢是一个变种 Pazdirek and Blaha 1996 其高灵敏度值 出现在测量电极之间的正下方 有适当的水平和纵向分辨率 但探测深度小 在三维电法难以单一使用 I U k MN s annk 1 可以联合使用这些装置 有的程序可联合反演 二 资料处理与反演解释 1 统计处理 统计处理包括以下内容 1 利用滑动平均计算视电阻率的有效值 例如三点平均 3 1 1 iiii sssx 式中 1 2 3 为 点的视电阻率有效值 i n D i x i 2 计算整个测区或某一断面的统计参数 平均值 为某一测区或某一断面上的测点数 N i xx i N 1 1 N 标准差 相等马 nni N i sA 2 x 2 1 ni N i xsA 2 1 3 计算电极调整系数 LiLK sx 其中为电极距为 L 时全部视电阻率观测数据平均值 L s 4 计算相对电阻率 LiiLKi sxxxy 通过计算相对电阻率 可以在一定程度上消除地点断面由上到下水平地层 的对变化 因此 相对电阻率断面图照顾要反映地电体沿剖面的横向变化 5 对视参数分级 为了对视参数进行分级 首先必须按平均值和标准差关系视参数的分级间 隔 间隔太小 等级过密 间隔太大 等级过稀 都不利于反映地电体的分布 一般情况下 以采用五级制为宜 即根据平均值和标准差的关系划分四个界限 Ax D 1 3 2Ax D 3 3Ax D Ax D 4 利用上述视参数的分级间隔 可将断面上各点的或划分成不同 i s i y 的等级用不同的符号或灰阶 灰度 表示时 便得到视参数异常灰度图 如 低阻 1 Di s 较低阻 21 DDi s 中等 32s D D i 较高阻 43 DDi s 高阻 4 Di s 视参数的等级断面图在一定条件下能比较直观和形象反映地点面的分布特 征 统计处理原则上适应于三电位电极系中各种电极排列的测量结果 只是在 考虑视电阻率参数图示时 由于偶极和微分两种排列的异常和地电体之间具有 复杂得对应关系 因此一般只对温纳装置的测量结果进行统计处理 当然 随着现代高密度电法仪装置的增加 温纳 斯伦贝谢装置的测量结果也可进行统 计处理 2 比值参数 高密度电阻率法的野外观测结果出了可以绘制相应装置的视参数断面图外 根据需要还可绘制两种比值参数图 考虑到三电位电极系中三种视参数异常的 分布规律 选择了温纳 装置和温纳 装置两种装置的测量结果为基础的一 类比值参数 该比值参数的计算公式为 iiiT ss 由于温纳 和温纳 这两种装置在同一地电体上锁获得的视参数总是具 有相反的变化规律 因此用该参数绘制的比值断面图 在反映地电结构的分布 形态方面 远比相应装置的视电阻率断面图清晰和明确的多 图 是对所谓地下石林模型的正演模拟结果 模型的电性分布已如图所 示 其中温纳 装置的拟断面图几乎没有反映 而比值断面图则清楚地 sT 反映了上述模型的电性分布 另一类比值参数是利用联合三极装置的测量结果为基础组合而成的 其表 达式为 1 1 1 ii ii ii B s A s B s A s 式中和分别表示剖面上相邻两点视电阻率值 计算结果示于 和 i s 1 i s i 点之间 比值参数反映了联合三极装置歧离带曲线沿剖面水二乘向的变1 i 化率 图 表征比值参数在反映地电结构能力方面所作的模拟实验 视电阻率断面图只反映了基底的起伏变化 而比值断面图却同时反映了基 s 底起伏中的低阻构造 2 高密度电阻率法二维地形边界元数值解法 高密度电阻率法是常规电法的一个变种 就其原理而言 与常规电法完全 相同 它仍然以岩 矿石的电性差异为基础 通过观测和研究人工建立的地中 稳定电流分布规律 解决水文 环境与工程地质问题 高密度电阻率法的正演 问题就是传导类电法的正演问题 也就是求解稳恒点电源电流场的边值问题 对二维地形 设起伏地面下均匀各向同性介质的电阻率为 具有电流强 1 度为 I 的稳恒点电流源位于地面任一点 A x 0 z 域 的边界由和组成 1 2 见图 1 根据位场理论可知 在有源域内及其边界任意一点 M x y z 处的电位 U x y z 满足 控制微分方程 M 1 FzyxU 2 自然边界条件 M 2 zyxQ n U zyxQ 1 本质边界条件 M 3 zyxUzyxU 2 式中 和分别是边界和上已知边值函数 这里 2 1 AMIF QU 1 2 为源点 A 到场点 i 间的距离 0 Q Ai r I U 1 2 1 Ai r 可以看到 地形是二维的 即沿 y 轴无变化 而电位 U 是 y 偶函数 所以 我们也把上边值问题称为 2 5 维问题 对 1 2 3 式进行余弦傅立叶变 换可得 控制微分方程 M 4 fzxuzxu 22 自然边界条件 M 5 zxqzxq 1 本质边界条件 M 6 zxuzxu 2 式中 为第二类零 1AA zzxxIf 0 q 2 0 1 A rK I u 0A rK 阶修正贝塞尔函数 是余弦傅立叶变换量或波数 这样 便将三维偏微分方程 1 变成了而二维偏微分方程 4 即将三维 空间的电位 U x y z 变换为二维空间的变换电位 u x z 为求得 u 可采 用边界元法求解 u 所满足的亥姆霍兹方程 4 式 借助格林公式及二维介质亥 姆霍兹方程的基本解 即可把 u 满足的亥姆霍兹方程及边界条件等价地归化为 如下的边界积分方程 7 dquduqrK I u Ai i 11 0 1 42 式中 为边界点 i 对区域的张角 为亥姆霍兹方程的基本解 i u k 为波数 为点到点 2 1 0 rKu cos 2 1 rnkrK k n u q r ii zx x z 的矢径 为边界的外法线方法 为第二类一阶修正贝塞尔函数 n 1 krK 采用边界元离散技术 将域 的边界进行剖分 分成个单元 根据 1 1 N 积分的可加性 7 式中对边界的积分可化为对每个单元上的积分之和 1 j 8 dquduqBuc jj N j N j iii 1 1 11 式中 2 i i c i Ai rK I B 0 1 4 方程组 8 式仅是含有未知量的线性方程组 解此方程组即可求得变 1 N 换电位值 u x z 然后按下式 0 cos 2 dyzxuzyxU 进行傅氏逆变换 即可求得电位值 U x y z 根据所采用的高密度电阻率法装置类型 逐点计算出某记录点处的纯地形 异常视电阻率值 然后用 比较法 进行地形改正 地形改正公式如下 D s 1 D ss G s 式中地形改正后的视电阻率值 该记录点实测的视电阻率值 纯 G s s D s 地形影响值 它是一个无量纲的标量 一般取 1 1 m 利用边界单元法计算高密度电阻率法地形边界位场问题是很有效的 但是 在算法引入时 必须针对高密度电阻率法的特点 作一些技术处理 高密度电 阻率法电极排列密集 并且采用了差分装置 所有这些特点都要求计算精度高 运算速度快 另外 所形成的矩阵也因测量电极到供电电极的距离变化很大呈 带状分布 并且当波数较大时 矩阵中的系数几乎都接近于零 造成解的不稳 定 为了解决这一问题 采用了增广矩阵法求解方程组的效果较为满BHU 意 为了保证精度 同时又减少运算次数 除采用九波数傅氏反变换外 还采 用了不等分单元剖分方案 具体做法是 在测线外 越远则单元剖分长度越大 且为最小电极距的整数陪 在测线段 则以最小电极距长度划分边界单元 为 了避免 r 等于零时贝塞尔函数无穷大的问题 剖分结点应不与电极点位置重合 最好选取相邻电极的中点为结点 2 12 1 二维地电构造中点电流源场的正演有限元算法的基本原理二维地电构造中点电流源场的正演有限元算法的基本原理 4 4 有限单元的思想最早是由 Courant 于 1943 年提出 倪光正 20 世纪 50 年代初期 由于工程分析的需要 有限元法在复杂的航空结构分析中最先得到 应用 而有限元法 Finite Element Method 这个名称则由 Clough 于 1960 年在 其著作中首先提出 半个世纪以来 以变分原理为基础建立起来的有限元法 因其理论依据的普遍性 不仅广泛地应用于各类结构工程 而且作为一种声誉 很高的数值计算方法已被普遍推广并成功地用来解决其他工程领域中的问题 如热传导 渗流 流体力学 空气动力学 岩土力学 机械零件强度分析 电 气工程问题等 20 世纪 70 年代初 J H Coggon 首先将有限元法用于电法勘探 后来 L R Rijo 完善了有限元法数值计算方法 使之成为正演模拟计算的有效方法 罗延钟 80 年代初 我国的周熙襄等在引进 Coggon Rijo 的算法时 将 Dey 和 Morrison 用于有限差分模拟的混合边值条件引入了有限单元法 从而发展了 Coggon Rijo 的有限单元算法 此后 罗延钟等在选用边值条件和反傅氏变换 的算法及波数取值等方面 又有了一些新的发展 使整个算法更臻完善 能在 不做任何校正的情况下 对相当大范围内变化的电极距 取得较高精度的计算 结果 90 年代中期 杨进提出了迭代有限元算法 它不仅模拟复杂地球物理模 型的能力强 模拟精度高 且占用计算机内存小 是对有限元方法的又一大发 展 以前的二维地电构造中点电流源场的有限元算法都是建立在常规直流电法 基础之上的 而高密度电阻率法因其装置的特殊设计 所以用有限元进行正演 计算时必须兼顾速度和精度 本文详细地研究了二维地电构造中点电流源场的 正演有限元法的理论基础 并结合高密度电阻率法的特点 在有限元算法本身 及其相关技术上作了许多改进 使之成为高密度电阻率法正演计算的有效算法 2 1 12 1 1 稳定电流场微分方程边值问题及其相应的变分问题稳定电流场微分方程边值问题及其相应的变分问题 高密度电阻率法二维正演计算问题即点 源场二维地电断面的边值问题 取三维笛卡 儿坐标系的 X 轴为垂直于地质体走向的方向 Y 为轴为地质体走向的方向 Z 轴垂直指向下 我们知道 二维情况下一般设地下导电性沿 Y 轴方向无变化 即 点电流源 zx 位于 XOZ 平面上 这时 虽然地下 AA ZxA 介质的电导率是二维的 而场源是三维的 电位函数 u x y z 所满足的微分方程及边界 条件为 f z zyxu zx zy zyxu zx yx zyxu zx x 2 1 式中 I 为 A 点的电流强度 AA zzyxxIf 在地面边界上应满足第二类边界条件 即诺依曼条件 在电法勘探中也可 1 称为绝缘边值条件 0 1 n u 2 2 在其余边界上 例如在离电流源很远的地方 可设 2 2 3 0 2 u 这种给定已知电位的边界条件称为第一类边界条件 即狄里希莱条件 或称极 限边界条件 也可以称作强加边界条件 第三类边值条件 实际计算时 计算网格毕竟是有限的 式 2 3 中的无穷远边界条件很难 满足 因此提出了第三类边值条件 对均匀半空间点源二维问题 地中任一点 电位的变化总有下面的一般形式 图 2 1 点源二维区域示意图 Figure2 1 Sketch map of point element for 2 D area 2 4 r c zyx c zyxu 222 式中 c 为常数 并设电源点放在坐标原点 r 为点源到计算点的径向距离 因 此有 2 5 cos 2 r u nr r c n u 式中为和边界外法线之间的夹角 上式也可以写为 r n 2 6a 0 r u n u 式中 上式便是第三类边值条件 或曰混合边界条件 cos 但是对于非均匀地电断面 便导不出 2 6a 式 为此 需采用更一般形 式的混合边界条件 即 0 r u n u 2 6b 式中 为与地电断面不均匀性等有关的一个系数 大量的试算发现取值在 0 1 之间 均值约为 1 2 关于系数的取值问题暂不作讨论 在下面的公式 推导中 为记述的方便 不妨暂时取 1 在所论由边界和包围的区域内 两种具有导电率为和的介质的交 1 2 1 2 界面处 电位和电流密度法向分量应满足以下衔接条件 或曰普遍边界条件 a 由于电位的连续性 在交界面处 有 2 7 21 zyxuzyxu b 由于电流法向分量的连续性 在交界面处 有 2 8 n u n u 2 2 1 1 其中为交界面的法线方向的单位矢量 n 前已谈到 对于点源二维问题 介质的电阻率是二维分布 而场源是三维 的 为使问题简化 我们常把三维问题通过傅氏变换化为二维问题 由于电位 是关于 y 轴的偶函数 故取余弦变换对 2 9 0 0 cos 2 cos 2 dxFxf dxxxfF 对 2 1 式按 2 9 式取傅氏变换 得 2 10 2 fzxvzx z zxv zx zx zxv zx x 式中 称为傅氏电位 变换电位或傅氏变 2 1 AA zzxxIf zxv 换电位 这样 便将三维微分方程 2 1 式变成了二维微分方程 2 10 式 当分 别对若干个给定的波数值求解方程 2 10 式 计算出傅氏电位 v x y 后 再按 2 9 式作傅氏反变换 即可计算出所求的电位 2 11 dyyzxvzyxu cos 2 0 在求解偏微分方程 2 10 时 于求解边界上 宜采用如下边值条件 在地面边界上 应用绝缘边值条件 即第二类边值条件 1 2 12 0 1 n v 在其余边界上 采用前述第三类边界条件 即混合边界条件 其导出过程 2 如下 由于求解区的左 右和底边界远离场源和电性异常体 故在那里的电场分 布可近似看成与均匀大地情况相同 设场源是位于 电流强度为 kkk zoxA 的个点电流源 1 2 3 则对于上述边界上的某点 P x y z 其 k INkN 电位可写成如下形式 2 13 N k k k yr IQ zyxu 1 22 式中 均匀大地的电阻率 2 Q 22 kkk zzxxr 对 2 13 式作傅氏变换 得 0 cos dyyzyxuzxv dyy yr QI N k k k cos 0 1 22 2 14 N k kk rKIQ 1 0 式中的 是零阶修正贝塞尔函数 对 2 14 式沿界外法线求方向导 0k rK n 数 并考虑到修正贝塞尔函数的求导性质 可得 10 xKdxxdK 2 15 N k kkk rKIQ n zxv 1 1 cos 式中 为一阶修正贝塞尔函数 为二维断面 x z 坐标面 上 场 1k rK k 源到所论点 P 的矢径与 P 点边界外法线方向之间的夹角 k A k r n 从 2 14 和 2 15 式中消去常数 Q 可得上述边界上的混合边值条件 2 16 0 2 n v v 式中 n k kk N k kkk rKI rKI 1 0 1 1 cos 于是 所提出的问题可归结为对若干个给定的波数 值 求解电位的傅氏 变换电位 v 所满足的下列二维偏微分方程的边值问题 v 满足二维偏微分方程 2 17 fvzx z v zx zx v zx x 2 其中 N k kkk zzxxIf 1 2 1 2 18 在地面边界上 1 2 19 0 1 n v 在其余边界上 2 2 20a 0 2 n v v 其中 1 20b n k kk N k kkk rKI rKI 1 0 1 1 cos 与上述二维偏微分方程边值问题 2 17 式 2 19 式 2 20a 式 等价的变分问 题为 极值 2 21 2 22222 2 dlvdsfvv z v x v vJ 从二维条件下的欧拉方程边值问题出发 我们可以证明 2 17 2 19 2 20a 与 2 21 式的等价性 并且地面边值条件 2 2 式在 泛函极值问题 2 21 式中没有出现 这是因为这一条件在泛函极值问题中能 自然得到满足 故称此边界条件为变分问题的自然边界条件 此外 在求解区 由边界围成的区域 内存在有限个电导率突变面时的内边界条件 21 2 7 2 8 式 在变分问题中也是自然满足的 不需另作处理 至此 用 有限单元法计算二维地电构造中的点电流源场的核心 就是用数值方法解变分 问题 2 21 式 求解变分方程 2 21 式 就是要找出一个傅氏电位的空间坐标函数 以使泛函最小 有限单元法就是用来求解这一问题的数值解的 zxv vJ 计算技术 它依据泛函欧拉方程边值问题与泛函极值问题的等价性 将微分方 程的边值问题转化为相应的泛函极值问题 对于稳定电流场 根据电场总能量 最小化原理 泛函的极值问题 2 21 式 应满足 02 2 22222 dlvdsvfv z v x v vJ 2 22 即泛函取得极值的必要条件是它的变分为零 如果将整个求解区间 划分为若干个单元 参见图 2 2 足够小 以致可以认为傅ee 氏电位在各单元内呈简单函数形式 zxv 比如 在二维条件下 可假设在 zxv 各三角单元内呈线性变化 于是 泛函的极 图 2 2 区域剖分示意图 Figure2 2 Sketch map of area dissecting 值问题 或变分问题 又可简化为多元函数的极值问题 而多元函数的极值问 题是大家所熟知的 这就是有限单元法的基本思想 求出傅氏电位后 利用傅氏反变换 2 11 式 即可求出电位函 zxv 数 对于高密度电阻率法来讲 一般 即要求求出主剖面上的电 zyxu0 y 位 再利用公式 zxu 2 23 I U K MN s 即可计算出各种装置的视电阻率值 有限单元法的优点是对形状不规则的地质体和地形起伏模拟能力强 但它 的计算程序编制相当复杂 要编制出一个效度好 适应能力强的高密度电阻率 法正演有限元算法程序 需要做很多理论研究和实际工作才行 2 1 22 1 2 变分问题的离散化变分问题的离散化 为了用有限单元法求解偏微分方程的边值问题 首先需要将其等价的变分 问题离散化 也就是要将 2 21 式离散化 这包括求解区的离散化 网格 剖分 和泛函的离散化 导出线性方程 两个组成部分 vJ 网格剖分 由于受计算机的计算量 储存量及计算速度的限制 所有的数值计算方法 都需要将无限大地中的电场分布限定在有限的求解区内进行讨论 并将连续的 求解区域离散化 即网格剖分 原则上讲 网格剖分可以根据所研究的地电条件灵活地将求解区限定为任 何形状 并灵活采用任意适当的网格来剖分求解区 在二维情况下 最常用的 办法是将求解区剖分为一系列互不重迭的三角单元 每个单元的顶点称为节点 剖分的规则是 1 整个剖分范围应该是越大越好 剖分范围越大 计算精度也 就越高 但剖分得越大 计算量增大 计算速度降低 因此 在作网格剖分时 既要考虑精度 又要照顾速度 2 各单元的节点只能与相邻单元的节点相重 而不能成为相邻单元的边内点 3 要使每个单元内介质的电导率为常数 在不 同电导率介质的内分界面上或整个求解区的外边界上 使三角形单元的一个边 相互衔接 以这样的折线去逼近边界线 4 三角形单元最好是接近正三角形 不要使其中一个角很小或出现很大的钝角 5 在电场变化剧烈 电位参数的二 阶或高阶导数大的地方 如 在异常体边沿 场源附近等 单元宜剖分的细一 些 而在电场变化平缓的地方 如 远离异常体和场源 单元面积可取得大一 些 当然 单元面积的变化应是渐变的 为使程序简化 L R Rijo 采用了图 2 3 所示比较规则的矩形求解区和三角 形剖分网格 试算结果表明 这种在矩形网格中布置交叉对称三角形剖分网格 可以足够近似地模拟一般常见的不平地形和电性异常体 同时又能节省计算量 是经典的一种剖分方法 在完成上述单元剖分的同时 于求解区内形成了若干个节点 包括内节点 和边节点 待定的傅氏电位 函数在这些节点上的 zxv 值 称为节点函数值 自然它 们是待定的未知值 这样 便 用求解区有限个待定节点函数 值来近似表征待定函数在连续 空间中的分布 这些待定的节 点函数值 1 v 2 v M v 总结点数 组成一组独M 立变量 所谓求待定函数的数 值解 就是要确定这些节点的函数值 泛函的离散化 1 线性插值 为了计算二维变分问题 2 21 式中积分形式的泛函 需要知道整个 vJ 求解区 内的 函数值 前面说到 求解区经过网格剖分后 可以用求解区内v 有限个待定节点函数值来近似表征待定函数在连续空间中的分布 这样 对于 二维情况 就可以利用各节点的函数值在各单元内作线性插值来逐个单元计算 也就是说 将 表示成该单元三个节点之函数值的函数 具体算法如下vv 设第 e 个单元的三个节点按逆时针方向的编号分别为 i j m 其坐标记 图 2 3 L R Rijo 网格剖分示意图 Figure2 3 Sketch map of Mesh dissecting by L R Rijo 为 对应的节点函数值为 参见图 2 2 ii zx jj zx mm zx i v j v m v 假设各单元内 函数是线性变化的 即 zxv 因为积分变量 是一个常数 不妨把 记为 czbxazxv zxv zxv 2 24 式中的的系数可由单元上的三个节点的函数值和坐标算出 对于单元 cba e 三个节点上都满足 2 24 式 故有 2 25 mmm jjj iii czbxav czbxav czbxav 按克莱姆 Cramer 法则求解上述线性方程组 得 mm jj ii mmm jjj iii zx zx zx zxv zxv zxv a 1 1 1 2 26a 2 mmjjii vavava 同样 有 2 26b 2 mmjjii vbvbvbb 2 26c 2 mmjjii vcvcvcc 式中 2 27 ijmjimijjim mijimjmiimj jmimjijmmji xxczzbzxzxa xxczzbzxzxa xxczzbzxzxa 单元面积 2 28 2 1 1 1 1 2 1 ijji mm jj ii cbcb zx zx zx 将 2 26a 2 26b 2 26c 式代入 2 24 式 便得到单元 e 内函数 线v 性插值的近似表示式 mmmmjjjjiiii vzcxbavzcxbavzcxbazxv 2 1 2 29 或 2 30 mmjjii vzxNvzxNvzxNzxv 式中的 称为基函数 分别为 i N j N m N 2 31 2 2 2 zcxbazxN zcxbazxN zcxbazxN mmmm jjjj iiii ezx 如果在每个单元上都作出了的这种近似 就能得到整个求解区内 zxv 的总体近似函数 这个函数在各个单元内是线性的 即它是一个分片为 zxv 线性的函数 对任意两个相邻单元来说 近似函数在公共边上的值 被两个端 节点的函数值唯一决定 所以 总体近似函数在整个求解区自然是连续的 2 单元分析 至此 二维变分问题 2 21 式中的泛函 整个求解区 上的积分 vJ 可以分解为各个单元 e 上的积分之和 即 2 32 vJvJ e e 而不难根据各单元 内函数 的线性插值近似表示式 2 29 或 2 30 vJeev 算出 因此 首先应从分析一个三角单元 e 出发 对于当前问题的变分提法 2 21 式 三个顶点分别为 i j m 的三角形单元 e 参见图 2 2 上的泛函 e e dlvdsfvv z v x v J 2 22222 2 ee ee dsvds z v x v 2222 2 33 e N k kkk dlvdszzxxIv 1 2 2 式中 表示单元 内的电导率值 由于它在每一个单元是常数 因此可以 e e e 提到积分外 同样 积分变量也可以提到积分外 根据 2 24 式及 2 26b 2 26c 有 2 34a 2 mmjjii vbvbvbb x v 2 34b 2 mmjjii vcvcvcc z v 可见 它们只与节点坐标及节点函数值有关 在一个单元内是常数 故可提到 积分号外 于是 可将 2 33 式中第一项积分表示为 ds z v x v e 22 e dxdz z v x v 22 2 35 22 4 1 mmjjiimmjjii vcvcvcvbvbvb 同时 根据狄拉克函数的性质 有 2 36 eA eA dxdzzzxx A e A 0 1 于是 2 33 式中积分的第三项为 2 37 e N k m mm j jj i ii kkk n vI n vI n vI dszzxxIv 1 式中 分别表示在 点的供电电流强度 若场源 供电电 i I j I m Iijm 极 不在某个节点上 则该点上的电流值为零 分别表示以节点 i n j n m ni 为顶点的单元数 2 33 式的第二项积分jm e mmjjii e dxdzvyxNvyxNvyxNdsv 2 2 e jj e ii dxdzzxNvdxdzzxNv 22 2 2 2 e imim e mjmj e jiji e mm dxdzzxNzxNvv dxdzzxNzxNvv dxdzzxNzxNvvdxdzzxNv 2 2 2 22 38 为了计算上式中的积分 我们首先来看基函数 zxNi 图 2 4 基函数几何意义示图 Figure2 4 Sketch map of base function geometric meaning 的几何意义 对于三角单元 e i j m 中的任一点 P x y zxN j zxNm 见图 2 4 由 2 27 式可算得 pjm zx zx zx z x x x z z zx zx zcxba mm jj m j m j mm jj iii 2 1 1 1 1 1 1 1 式中 为三角形的面积 故由 2 31 式 得pjm pjm 2 39 pjmzcxba zxN iii i 2 即表示以 P x z 为一顶点 为对边的三角形面积 与单元面 zxNijmpjm 积之比值 和 的几何意义可以类推 由此可得 zxN j zxNm zxNi 和 的下述性质 zxN j zxNm 2 40 1 zxNzxNzxN mji 1 41 0 0 1 mmi jji iii zxN zxN zxN 以及 1 42 0 1 0 mmj jjj iij zxN zxN zxN 另一方面 由 2 31 式 可得 1 43 2 1 2 1 zcxbaN zcxbaN jjjj iiii 可将 2 43 式看作是从 x z 到 的坐标变换 由 i N j N 2 41 和 2 42 式可知 这种变换将 x z 平面上的 点分别变换为 平面上的 mmjjii zxmzxjzxi i N j N 点 见图 2 4 其面积单元之比等于雅可比 Jacobian 0 0 1 0 0 1 mji 式 4 1 22 22 2 ijji jj ii jj ii ji cbcb cb cb z N x N z N x N dxdz dNdN 考虑到 2 28 式 2 ijji cbcb 故有 2 44 jidN dNdxdz 2 由此可计算 e e ijm mji jiii dNdNNdxdzzxN 22 2 1 0 1 0 2 2 i N j idNdNN i 1 0 2 6 1 2 iii dNNN 2 45 同样 e e ijm mji jijiji dNdNNNdxdzzxNzxN 2 1 0 1 0 2 1 0 2 1 2 22 ii i i N jj idNN N dNdNNN i 1 0 32 12 2 iiii dNNNN 2 46 可以将以上两式类推到 2 38 式中的其余积分 并写成如下一般形式 2 47 1 12 pqq e p dxdzzxNzxN 式中 p q 分别表示 i j 或 m 0 1 qp qp pq 将 2 47 式代入 2 38 式得 2 48 6 2222 immjjimji e vvvvvvvvvdsv 现在 只剩下 2 33 式中最后一项 线 积分了 而此项积分可仿照第二 项 面 积分的方法求得 应该指出 此 线 积分是由求解区边界上的边 2 值条件引入的 故只需对属于 的单元边界计算该积分 设单元边界 在 2 jm 求解区边界 上 可以近似地将 上的系数视为常数 取为 j 或 m 点的 2 jm 值 或 同时 电导率亦为常数 皆可提到积分号外 于是 问题 j m e 归结为计算积分 jmjm mmjjii dlvzxNvzxNvzxNdlv 2 2 根据前述的几何意义 可以写出它们在 边上的表示式 参见 mji NNN jm 图 2 5 0 pjmzxNi jmj llpmizxN 1 jmm llpijzxN 故有 jm l m jm j jm jm dlv l l v l l dlv 0 22 1 jmjmjm ll jmjm mj jm m l jmjm j dl l l l l vvdl l l vdl l l l l v 00 2 2 2 2 2 0 2 2 2 2 21 图 2 5 含边界的三角单元上线积分计算用图 Figure2 5 Map for Liner intergration of trigon unit on boundary 3 22 mjmj jm vvvv l 2 jm 2 49 如果或 位于边界上 则可类推得ijmi 2 2 50 jiji ij ij vvvv l dlv 222 3 2 ij 或 2 51 3 222 imim mi mi vvvv l dvv 2 mi 将 2 35 2 37 2 48 2 49 2 50 和 2 51 代入 2 33 式 得单元 e 上的泛函 22 4 mmjjiimmjjii e e vcvcvcvbvbvbJ 6 222 2 immjjimj i e vvvvvvvvv m mm j jj i ii n vI n vI n vI 3 3 3 2 22 2 22 2 22 mivvvvl jmvvvvl ijvvvvl imimmimie jimjjmjme jijiijije 2 52 这样 便将单元泛函表示成该单元各节点函数值 和的多元函数 i v j v m v 将所有单元的相加 便得到整个求解区的 vJe vJ 2 53 e e vJvJ 于是 前述变分问题 泛函的极值问题 就变为多元函数的极值问题 我们知 道 多元函数取得极值的必要条件是其偏导数为零 即对每一变量 vJ vJ 1 v 的偏导数为零 2 v M v 2 54 0 vJ v vJ vv vJ e e r e e rr 2 1 Mr 我们仍就从一个单元分析入手 由 2 52 式可计算包含三个节点mji 之单元 e 的泛函 Je 对各节点的函数值的一阶偏导数 0 1 v Je i mii iji iiiie i e v mil jm ijl ccbb v J 3 2 0 3 2 3 2 1 2 2 2 2 j ijj jijie v mi jm ijl ccbb 0 0 3 6 2 1 2 2 2 2 i i m mimi mimie n I v mil jm ij ccbb 3 0 0 6 2 1 2 2 2 2 iim e imj e iji e ii nIvKvKvK jjm e jmj e jji e ji j e nIvKvKvK v J mmm e mmj e mji e mi m e nIvKvKvK v J 0 M e v J 可以将以上各式写成矩阵形式 i 列 j 列 m 列 行 行 行 m j i nI nI nI v v v v v KKK KKK KKK v J v J v J v J v J mm jj ii M m j i e mm e mj e mi e jm e jj e ji e im e ij e ii M e m e j e i e e 1 1 矩阵 及列矢量 Ie 中的虚点均为零元素或记为 e K 2 55 eee IVKJ 式中 维列矢量 的元素为M e J r 1 2 r ee r v J J M 2 56 维列矢量的元素为各节点的函数值 r 1 2 M V r vM 维列矢量 Ie 的元素为M 2 57 0 mjir mjirnI I rre r 阶单元矩阵 的非零元素为MM e K 2 58 36 2 1 2 sr s srsre e rs lG G ccbbK 1 2 sr sr Grsmjirrs 表示式中最后一项仅当单元 e 的边在求解区边界上才出现 其中的 e rs K sr 2 是按 2 20b 式对节点 s 算出的系数值 s 3 总体合成 2 55 式是整个求解区中某一个单元 e 的泛函对各节点的函数值的一阶 e J 偏导数所形成的线性方程组的矩阵形式 将所有单元的相加 便得到整个 e J 求解区的泛函的变分 0 e e JJ 2 59 写成维列矢量M 2 60 e e JJ0 2 1 Mr v J JJ r r 的元素为 将 2 55 代入 2 60 式得 0 e e e e IVK 或写成 2 61 IVK 式中 刚度矩阵 其元素为所有单元矩阵的相应元素之和 e e KK e K 2 62 e e rsrs KK r s 1 2 M 场源列矢量 由 2 55 可知 其元素为节点 r 上的供电 e e II e e rr II 电流强度 r 1 2 如果供电电极不在节点 r 上 则 M0 r I 2 61 式便是由变分问题 2 22 式离散化后导出的线性方程组 解此 方程组便可对给定的波数确定各节点电位的傅氏变换电位 并可通过反傅 v 氏变换计算各节点的电位值 然后 根据公式 2 23 即可计算出视电阻率u 值 2 22 2 二维地电构造中点电流源场的正演有限元算法计算的相关二维地电构造中点电流源场的正演有限元算法计算的相关 技术技术 2 2 12 2 1 方程组的简化方程组的简化 采用图 2 3 所示的三角形网格能对求解区作比较精细的剖分 因而模拟地 形或电性异常体及单元内作线性插值的近似 程度都较好 但其缺点是节点数目较多 所 形成的线性方程组阶数较高 因而计算量较 大 下面介绍一种算法 把图 2 3 中所示网 格中位于各矩形中心之节点函数值从方程 2 61 中消除 以降低待求解的线性方程 组的阶数 从而大大减少计算量 我们考察图 2 3 所示网格中的一个小矩 形单元 图 2 6 其四个角的节点编号设为 i j m k 显然有E 2 63 jMm iMk 式中 为矩形网格的纵向节点数 参见图 2 3 该矩形单元的中心节点设为M a 其与四个角节点相连 将此矩形单元分为四个三角形单元 依次表示为 1 e 图 2 6 网格中的一个小矩形单元 Figure2 6 Little rectangle unit in mesh 和 它们的面积都相同 且为 2 e 3 e 4 e 2 64 zxh h 4 1 其中 和分别为矩形单元沿和方向的步长 x h z hxz 设各三角形单元内电导率是均匀的 但各三角形单元之间 电导率可以不 相同 分别以 表示 内介质的电导率 按 1 2 3 4 1 e 2 e 3 e 4 e 2 58 式可写出各三角形单元 的矩阵 和 1 e 2 e 3 e 4 e 1 e K 2 e K 3 e K 它们之和便给出矩形单元的矩阵 4 e KE E K 4321 eeeeE KKKKK 2 65 在计算这些矩阵的元素时应注意 在当前情况下 和 r i j m k 或 a 都 r b r c 可表示为十分简单的形式 比如 在中 按 2 26 式有 1 e 2 zaji hzzb 2 xiai hxxc 2 ziaj hzzb 2 xaiaj hxxc xjia hzzb 0 ija xxc 于是 按 2 58 式可算出的各非零元素为 1 e K 3 2 2 1 22 1 1 ii e ii cbK 43 44 2 1 2 1 22 1 zxxz zx hhhh hh zx z x x z hh h h h h 12 2 2 11 3 2 2 1221 1 jj e jj cbK 43 44 2 1 2 1 22 1 zxxz zx hhhh hh zx z x x z hh h h h h 12 2 2 11 3 2 2 1221 1 aa e aa cbK 43 0 2 1 2 1 2 1 zx z zx hh h hh zx x z hh h h 12 2 2 1 1 6 2 2 11 11 jiji e ji e ij ccbbKK 46 44 2 1 2 1 22 1 zxxz zx hhhh hh zx z x x z hh h h h h 24 2 2 11 6 2 2 11 11 aiai e ai e ia ccbbKK 46 0 2 2 1 2 1 2 1 zxz zx hhh hh zx x z hh h h 24 2 1 1 6 2 2 11 11 ajaj e aj e ja ccbbKK 46 0 2 2 1 2 1 2 1 zxz zx hhh hh zx x z hh h h 24 2 1 1 同样可计算出 的各非零元素 2 e K zx z x x ze jj hh