三维地形直流电阻率有限元法模拟

点电源电位满足的基本方程 ,详细讨论见文献[ 18 ] .
Δ
Δ
Δ
·(σ· u) = - ·j ,
(1)
其中 ,σ为电导率 , u 为电位 , j 为电流密度. 对 (1)
图 1 网格剖分示意图 Fig. 1 Schematic diagram for finite element mesh
1608
Z 方向的步长)
,设为
Δ z
.
三个
棱面都是平行四边形 ,其中 ,平面1254在坐标平面
本文研究以三棱柱为单元的有限单元正演算 法 ,可以实现任意起伏地形的电阻率模拟计算. 三 棱柱具有以下优点 :只要给定地面上任何三点的高 程 ,就可以惟一确定一块小面积 ,而且对高程的变化 没有特别限制 ,这就意味着三棱柱单元可以模拟任 意坡度的地形. 求解变分方程时 ,采用有限单元法 把目标区域离散成一系列三棱柱单元. 含有地形特 征的不规则三棱柱 ,有五个面六个节点 ,单元刚度矩 阵含有 6 ×6 个元素 ,其中有 21 个非零元素 ,每个元 素由 30 多项含有 5 个变量的三重积分组成 ,若手工 计算非常复杂且容易出错 ,本文通过编制积分算法 程序 ,只要计算一次就可得出 21 项非零元素的数学 表达式 ,可节约推导时间.
(3)
∫ F( u)
=
Ω
1 2
σ(
ቤተ መጻሕፍቲ ባይዱ
Δ
u) 2 - 2 Iδ( A) u dΩ
∫ +
1 2
Γ∞
σcos
(r r
,
n)
u2
dΓ
,
(4)
δF ( u) = 01
3 有限单元法
311 网格剖分 为较好地模拟不平地形 ,同时又尽可能使网格
简单化和规则化 ,并考虑以后反演方便 ,正演模拟中 采用了三棱柱体网格单元. 三维有限单元模拟的三 棱柱 网 格 由 中 心 区 网 格 和 扩 展 区 网 格 两 个 部 分 组成.
长 ,布置 10～15 个节点 ,形成扩展区网格.
312 任意三棱柱单元分析
按前一节所述 ,三维有限单元模拟的网格被剖
分为一系列直立三棱柱体. 图 2 示出了其中的一个
有代表性的三棱柱体单元. 单元的三个棱14 ∥25 ∥
36 ∥OZ ,即都沿铅垂方向 ,平行于 Z 轴 ; 它们的长
度相同 (即为单元沿
5 期
强建科等 :三维地形上直流电阻率有限元法模拟
1607
1 引 言
在进行电法勘探时 ,起伏地形严重影响直流电 阻率法或电磁法的应用 ,有时甚至使异常面目全非. 因此 ,用数值方法计算地形对视电阻率的影响并设 法消除它 ,是提高资料地质解释效果的关键之一.
国内外专家已经作过一些有关起伏地形三维电 阻率正演模拟问题的研究工作 ,Fox 等[1] . 阐述了地 形对电阻率和极化率的影响问题 ; Truman H T 和 George R J [2] 研究了三维地形对电阻率的影响并进 行校正 ; Yuguo L 和 Klaus Spitzer[3] 比较了三维直流 电阻率有限单元和有限差分算法研究; 黄兰珍 等[4～7] 、徐世浙等[8～12] 用边界元法较好地解决了地 形影响电阻率的正演问题 ,但当地下物性分布复杂 时 ,边界元法并不具优势. Sasaki[13] 、熊彬等[14] 采用 四面体单元剖分目标区的有限单元法实现了三维起 伏地形电阻率的正演模拟计算 ,缺点是单元数增大 6 倍 ,不利于进一步做反演研究. 吴小平等[15～17] 利 用有限差分法 、有限单元法开展了起伏地形下电阻 率Π激发极化 32D 正 、反演技术研究 ,但该算法的正 演计算建立在以长方体为单元的基础上 ,对实际地 形的模拟存在明显不足 ,因为地表四个点不能惟一 地确定一个面.
三维地形直流电阻率有限元法模拟
强建科1 ,2 ,罗延钟2
1 中南大学信息物理工程学院 ,长沙 410083 2 中国地质大学地球物理与空间信息学院 ,武汉 430074
摘 要 基于稳定电流场的基本方程 、三维区域满足的边值问题以及相应的变分问题 ,研究了三维起伏地形条件
下电阻率的有限单元数值模拟算法. 离散积分区域时 ,以三棱柱为最小研究单元 ,推导了含有地形特征信息的三线
垂线上按测量电极距等距分布的点 ,作为中心区网
格的地下节点. 当地面水平时 ,这样形成的中心区
网格 ,为节点均匀分布的正方体网格 ; 当地面起伏
不平时 ,水平方向呈正方形分布的四个相邻地面节
点可能不在一个平面上 ,因而不能组成正方体网格.
为避免出现上述情况 ,连接上述四个相邻地面测点
的两个对角点 (约定连接左前方到右后方的对角
第 50 卷 第 5 期 2007 年 9 月
地 球 物 理 学 报
CHINESE JOURNAL OF GEOPHYSICS
Vol. 50 , No. 5 Sep. , 2007
强建科 ,罗延钟. 三维地形直流电阻率有限元法模拟. 地球物理学报 ,2007 ,50 (5) :1606～1613 Qiang J K, Luo Y Z. The resistivity FEM numerical modeling on 32D undulating topography. Chinese J . Geophys. (in Chinese) , 2007 ,50 (5) :1606～1613
点) ,将地面节点在水平方向上划分成一系列直角三
角形 ,从而将中心区网格剖分为一系列直立三棱柱
体单元. 这些三棱柱体的三个棱沿铅垂方向 ,上底
和下底为相互平行的 、随地面起伏变化的倾斜直角
三角形.
为减小网格边界的影响 ,从中心区网格分别向
左 、右 、前 、后和下方 ,按公比值 113 的等比数列取步
topography
基金项目 国家自然科学基金项目 (60672042) ,中国地质大学 (北京) 地下信息探测技术与仪器教育部重点实验室开放基金 ( GDL0502) ,中南大 学博士后基金联合资助.
作者简介 强建科 ,男 ,1967 年生 ,中南大学在读博士后 ,主要从事地球物理电磁法正反演算法研究. E2mail :qiangjianke @1631com
地 球 物 理 学 报 (Chinese J . Geophys. )
50 卷
格节点的个数为 ( M + 1) ×( N + 1) ×( L + 1) .
在三维电阻率法勘查时 (如高密度) ,通常采用
正方形测网 ,即测线距和测点距相同. 以水平方向
呈正方形分布的地面测点 ,作为三维有限单元模拟
中心区网格的地面节点 ;由其向下引铅垂线 ,并取铅
算精度满足误差要求. 在二维山脊上的二极剖面或三维山谷上的中间梯度剖面上 ,其三维计算结果与相应模型的
土槽实验结果或边界元法计算结果也非常接近.
关键词 电阻率法 ,有限单元法 ,三维数值模拟 ,起伏地形
文章编号 0001 - 5733(2007) 05 - 1606 - 08
中图分类号 P631
Abstract Based on stable electricity flow field equation , three2dimension boundary question and its corresponding variation question , the algorithm of resistivity FEM numerical modeling on a 32D topography is studied. The integral area is divided into many tri2prism units with anomalous shapes and landform characteristics. The function of tri2linear interpolation and the matrix of tri2prism unit are then deduced. Using varied band matrix and one2dimension array to collect the large sparse matrix , the integral coefficient of all the nodes , we could save memory space. Before solving the system of linear equations , the Cholesky algorithm was adopted to decompose the large sparse matrix only once , and then we obtained the potentials of all nodes of hundreds of different current electrodes by backward substitutions. This is an important step to improve the speed of forward modeling when there are more sounding2points. The results of this modeling show that the 32D numerical solution is in agreement with the analytic solution on level stratified medium , on 22D mountain ridge with pole to pole array , or on 32D mountain valley with central gradient array. The 32D numerical solution is also in agreement with the results from physical experiment or boundary element method. Keywords Resistivity method , Finite element method ( FEM) , 32D numerical modeling , Undulating
式做数学变换 ,再依据半空间δ函数的积分性质得
到三维微分方程 :
9 σ9u 9x 9x
+
9 9y
σ9 u 9y
+
9 9z
σ9 u 9z
= - 2 Iδ( xA )δ( yA )δ( zA ) .
(2)
三维电场的边值问题是
9u 9n
=
0,
9u 9n
+
cos
(
r r
,
n)
u
=
0,
等价的变分问题是
∈ΓS (地表) ∈Γ∞ (地下边界) .
收稿日期 2006 - 12 - 20 ,2007 - 05 - 17 收修定稿
The resistivity FEM numerical modeling on 32D undulating topography
QIANG Jian2Ke1 ,2 , LUO Yan2Zhong2
1 School of Info2physics and Geomatics Engineering , Central South University , Changsha 410083 , China 2 Institute of Geophysics & Geomatics , China University of Geosciences , Wuhan 430074 , China
性插值型函数以及单元刚度矩阵. 采用变带宽 、一维数组方式只存储稀疏刚度矩阵中非零元素 ,能够节约内存. 利
用 Cholesky 分解法只分解一次大型稀疏矩阵 ,通过回代可以求出方程组的全部解 ,当求解有多个供电点的测深问题
时可以缩短计算时间. 模型计算表明 ,在水平层状介质模型上 ,三维计算结果与解析解或二维数值解十分吻合 ,计
考虑到电法的野外工作方式 、计算精度 、计算量 以及使程序简化 ,采用图 1 形式的单元剖分方法 :在 每一个长方体单元内 ,进一步剖分出交叉对称的三 棱柱单元. 其中 ,每个三棱柱单元中的电导率是均 匀的. 三棱柱单元个数为 M ×N ×L ×2 ,三棱柱网
2 稳定电流场电位满足的基本方程
为了后面讨论方便 , (1) 式列出有关稳定电流场