均匀半空间视电阻率MATLAB代码,基于Matlab的直流电阻率三维数值模拟及可视化

电阻率MATLAB代码,基于Matlab的直流电阻率三维数值模拟及可视化"/>

均匀半空间视电阻率MATLAB代码,基于Matlab的直流电阻率三维数值模拟及可视化

基于Matlab的直流电阻率三维数值模拟及可视化

第2 9 卷第3 期 成 都大 学 学 报( 自然科 学 版) 2 0 1 0 年 9 月 J o u r n a l o f C h e n g d u U n i v e r s i t y ( N a t u r a l S c i e n c e E d i t i o n ) V_01 ． 29 No． 3 S e p ． 2 01 0 文章编号： 1 0 0 4—5 4 2 2 ( 2 0 1 0 ) 0 3 —0 2 6 1 — 0 4 基于 Ma t l a b的直流 电阻率三维数值模拟及可视化 彭 淼 ，熊 杰 ( 中国地质大学( 北京) 地球物理与信息技术学院，北京 1 0 0 0 8 3 ) 摘要： 针对直流电阻率法的三雏正演问题， 采用积分方程法的数值解法， 特别地， 对极状体三雏模型做 了具 体的网格剖分， 在 Y m t h b下编程实现了求解其电位分布的算法． 根据电位分布分别采用中闻梯度法， 三极剖面 法和联合剖面法计算 了地面的视电阻率分布． 结果表明， 视 电阻率异常图充分显示了积分方程求解的正确性 ． 对于该结果， 在 M a t l a b中采用 G U I 编程生成了更为方便的人机交互式正演界面， 该界面具有直观、 高效的特点． 关键词： Ma t l a b ； 积分方程； 直流电阻率法； 三维正演 ； G U I 中图分类号 ： P 6 3 1 ． 3 ’ 2 5 文献标识码 ： A 0 引 言 直流电阻率法是地球 物理勘探中一种常用 的方 法． 它是以岩矿石导电性差异为基础 ， 通过观测 、 研 究人工建立的地下稳定电流场 的分布规律以达到找 矿和解决其他地质问题为目的的一种电法勘探方 法u ] ． 研究直流电阻率三维正演问题是电法勘探中 的一项重要任务， 近 3 0 年来国内外已有不少学者对 该问题采用有限单元法、 有限差分法、 积分方程法等 数值模拟方法进行了算法研究， 成果显著 一 】 ． 本文采用积分方程法进行直流电阻率的三维正 演 ． 对于该正演问题 ， 周熙襄等⋯详细阐述 了积分方 程的近视求解方法 ， 并对具体模型作了正演计算 ． 他 们实现该算法是用 F O R T R A N语言编写的， 而本文 采用可视化功能更为强大的 M a tl a b 来实现． 笔者不 仅最终成功实现了直流电阻率的积分方程法三维正 演算法的高效计算 ， 而且还采用 G U I 编程开发了正 演用户界面， 这使得方便、 快捷的人机交互式正演成 为可能 ． 1 积分方程及其求解方法 假设在无限大半空间的均匀介质中， 存在一个 三维地质体， 它的电性与周围介质的电性不同． 位于 地面的一个电流源供给稳定电流 ， ， 均匀介质的电 阻率为 ID ， 三维地质体的电阻率为 P ， 表面电荷密 度为 g ． 根据基本的电场和电磁场定律及边界条件， 可 以导 出如 F 积分 方程 ： 孥 q ( P ) = [杀 (去 ) ] + 杀 (击 ) d s + q杀 c ㈤ ．； L： ( 2 ) P2+ P P ) = ／P l + 景 d s + Ⅱ ( 3 ) 式( 1 ) 、 ( 2 ) 为电阻率突变界 面积累电荷密度 q 所满足的积分方程式 ． 当外 电源 和电阻率分 布给定 时， 求解式( 1 ) 、 ( 2 ) 即可得到 g的分布． 通过已知的 g分布求解式 ( 3 ) 便可以计算出空间任意点的电位． 对于上面的积分方程， 要给出其解析解是很困 难的， 本文采用数值计算的近似解法， 其具体步骤如 下 ： 首先 ， 将积分域离散化为有限个小面积的组合 ． 然后 ， 将每一个小面积用一个简单的平面近似代替 ． 设小平面内的口 为常数， 这样， 积分方程就变成了一 组线性方程组： 2 ； 7／“ 口 = I p 丌 , [ - a ；( 1 ) ] + 砉 杀 ( 去 ) d s + 骞 杀 c ㈤ 式中， g i ( 或 ) 为第 个( 或第．『 个) 小平面内的积累 收稿 日期 ： 2 0 1 0—0 7—2 5 ． 作者简介：彭 淼( 1 9 6 3一 ) ， 男，博士研究生 ， 从事电法勘探技术研究 ． · 2 6 2 · 成都 大学学报( 自然科 学版 ) 第 2 9 卷 电荷密度 ， 它在该小平面内为一常数， ， ( 或 △一S j ) 为小面元的面积 ． 式( 4 )中， 当 i=1 ， 2 ， ⋯， ／7 ． 时， 可以得到 几个线 性方程 ， 组成一个 n阶线性方程组， 其矩阵方程为 ， A = ， 其中， 系数矩阵 A的对角线元素为 ， 一 Ⅱ O ( 1 ( 5 ) △ 』 “ 非对角线元素为， 矗 c 一 杀 c ㈤ 右端项为 ， = [杀(击 ) ，杀( 击 一 ，杀(击 (7 ) 然后 ， 对系数矩阵各元素中的面积分项进行 计 算 ， 有以下近似表达式 ： Ⅱ 若 ( 枷s = O 1 (8 ) 杀 c = 毒 将式( 3 ) ( 9 ) 联立起来求解 ， 就可以得到 q的 分布， 此即为积分方程( 1 ) 的数值近似解 ． 2 板状体三维模 型 上述积分方程的数值近似求解方法适用于任意 形体的三维模型计算． 本文采用的是板状体三维模 型( 见图 1 ) ． 图 1中 1 - _ 部的长方体代表虚源． ， ( 。 ． 一： J ) B 图 1 板状体三维模型示意图及表面网格剖分 对于一个长、 宽、 高分别为 A 、 B 、 C的板状体， 本文采取以下方式划分表面的小面元 ： 把长方体 的 每条棱划分成 Ⅳ段 ， 其中 口=A ／ N， b= B ／ N， c= C ／ N， 那么共 出现3 种不同的小面元 ， 0×b ， b×c 和 c×n ， 一共有 6 X N z 个小面元 ． 这样划分的优点是 避免 了奇异面元的讨论 ， 从而能很方便地实现算法 ． 由于每一个划分得足够小 的面元 内的积累电荷 密度 q 都可以视为常数 ， 故对 于式( 4 ) 、 ( 5 ) 、 ( 6 )中 的积分项 ， 可以采用下列近似公式求解 ： 旦 ‘ 1 on , )： c o s a+ c o s 』9+ ～ L ， ∞鲫 + ∞。 + c 。 s y ( 1 0 ) ( )： c o s + c o s J8+ d R R3 R3 。 。 。 ’ ， ( 1 1 ) 3 矗 (击 ) = co s¨ c o s．19 + __ 一 一 ∞ 。 + 可 ∞ 。 + c 。 s y ( 1 2 ) 通过对板状体上顶面、 底面、 前侧面、 后侧面、 左 侧面以及右侧面的分类讨论 ， 计算系数矩阵和右端 项 ， 并求解线性方程组可得到 q的分布， 然后利用下 式计算地表任意点的电位 ： P)= +2 s ) 最后 ， 本文在 Ma t l a b的编程环境 中实现了计算 地表电位分布的算法 ， 并将程序模块保存于相应文 件 中． 3 视电阻率的计算及算例 为了检验程序的有效性和正确性 ， 本文分别采 用中间梯度法 、 三极剖面法以及联合剖面法计算得 到