复杂地形电导率线性变化二维大地电磁有限单元法正演模拟_欧东新

[
1 2
]，
（ 14 ）
uT um ） 。 l = （ uj ， （ 6）
T 其中： l 为三角形单元位于边界 CD 上的边长； u l 为
三角形单元位于边界 CD 上的节点场值向量。 3. 2 TM 模式系数矩阵 TM 模式离散的积分形式为
F（ u） =
∑ {
Ω
△
1 · （ u） [1 2 σ － iωε
成，但算法复杂， 并且为了变带宽系数矩阵的存 储还要对网格节点编号进行优化，以便半带宽最 小。刘云等采用矩形网格变形技术
［10 ］
2变分问题根据来自 1a 坐标系，徐世浙［13 ］
模拟起伏地
给出二维大地电
形，这种方法相对于 Delaunary 算法等自动剖分方 法简单有效，并且容易对节点单元编号， 形成较 小的半带宽，但没有给出模拟垂直地形的方法 。 下面给出垂直间断面地形的三角单元网格剖 分具体过程： ①测量 地 形 高 程 和 各 个 垂 直 间 断 面 的 断 距， 形成地形拉平断距不变的矩形网格 （ 图 2a） 。 ②按照地形起伏改变地形线上节点的 y 坐标， 同一条竖线上的所有节点同时在纵向移动相同距 离，形成平行四边形单元 （ 图 2b） 。 ③将平行四边形剖分成三角形， 把上下边界
其中： Ω 为 ABCD 内所有计算区域； E z 为电场强度； 变平，最终形成图 1b 的三角网格。对于 TM 模式， H z 为磁场强度； ω 为圆频率； μ 为磁导率； σ 为电导 空气层是不必要的， 但是为了编号的方便可以采 率； ε 为介电常数； k = 槡－ iωμσ ， 为传播系数。 。 用较薄的空气层 ④将空气层的单元的各个节点电导率赋为 0 ； 地面节点电导率值有 2 个， 分属不同单元。 如果 不作这步处理， 在垂直间断面的转角四边形单元 中，只有一个节点电导率为 0 ，单元内电导率插值 将形成连续地形而不是垂直地形 。
ai ai aa ［ A i， j］ = j i am ai bi bi bb B i， ［ j］ = j i bm bi

△
1 ［ B j］ dxdy； σ i N i + σ j N j + σ m N m － iωε i， （ 18 ） 对称 aj aj am aj 对称 bj bj bm bj ， am am 。 bm bm （ 19 ）
（ 10） 其中： σ i ～ σ j 为三角形顶点的电导率值； △ 为三角 形面积。 式 （ 4 ） 中第 3 项积分为



[( ) ( ) ]
(
)

k e2 =
△ · 60
k2 =
对称 6σi + 2σj + 2σm 2σi + 2σj + σm 2σi + 6σj + 2σm 。 2σi + σj + 2σm σi + 2σj + 2σm 2σi + 2σj + 6σm
1 4 △2
2
－
1 iωμu2 dxdy + 2
]
}
y i ～ y m 为三角形顶点坐标； △ 为三 其中： x i ～ x m ，
∑∫
CD
l
1 1 · · － iωμσ ·u2 dl， 2 σ － iωε 槡
（ 15）
在单元内利用面积坐标对电导率进行线性插
值，由于电导率位于分母，无法直接求出 TM 模式 很容易推导出 TE 模式下线性插值系数矩阵， 下系数矩阵解析解， 采用高斯数值积分方法获得 这里直接给出 TE 模式单元积分系数矩阵公式： 单元系数矩阵。 式（ 4 ） 中第 1 项积分为 式 （ 15 ） 中第 1 项单元积分为 1 1 1 · （ u ） 2 d x d y = u T （ 7） e k e1 u e ； 1 1 iωμ 2 △ 2 · （ u ） 2 d x d y σ － iωε △ 2 1 1 2 2 k e1 = ［ k t］ ，k s， （ a a + bs bt ） ， t = 1 1 u + u d x d y iωμ s， 2△ s t = · σ － iωε x △ 2 y s， t = i， j， m。 （ 8） 1 1 T = uT （ k + k2 ） u e = u T k u； （ 16 ） 其中： u e = （ u i ，u j ，u m ） 为三角形顶点电场值向 2 e 1 2 e e1 e k e1 为 3 × 3 的系数矩阵。 量。 对电导率和场值在单元内线性插值得 式 （ 4 ） 中第 2 项积分为 1 1 k1 = ［ A i， dxdy； j］ 2 1 T 1 2 4 △ σ i N i + σ j N j + σ m N m － i ωε △ （ 9） － σu dxdy = － u e k e2 u e ， 2 2 △ （ 17 ）
收稿日期： 基金项目： 作者简介： 引文格式：
自动生
2013 － 04 － 09 广西自然科学基金项目 （ 0848021 ） ； 广西地质工程中心重点实验室开放基金项目 （ 桂科能 07109011 － K003 ） 欧东新 （ 1974 —） ，男，博士，副教授，研究方向： 地球物理正反演，odx@ glut. edu. cn。 ． 桂林理工大学学报，2014 ， 欧东新，韦者良． 复杂地形电导率线性变化二维大地电磁有限单元法正演模拟 ［J］ 34 （ 1 ） ： 30 － 38．
第 34 卷 第 1 期 2014 年 2 月
桂 林 理 工 大 学 学 报
Journal of Guilin University of Technology
Vol. 34 No. 1 Feb. 2014

文章编号： 1674 － 9057 （ 2014 ） 01 － 0030 － 09 doi： 10. 3969 / j． issn. 1674 － 9057. 2014. 01. 005
磁的变分问题的公式 1 1 2 2 F（ u） = τ （ u ） － λ u d Ω + 2 2 Ω 1 2 τku dΓ， 2 CD = 1， u AB δF（ u） = 0 。
[
]
∫
（ 1）
对于 E 型波 u = E z ，τ = 1 / （ iωμ） ，λ = σ － iωε。 对于 H 型波 u = Hz ， τ = 1 / （ σ － iωε） ， λ = iωμ。 （ 3 ） （ 2）
具体三角单元高斯积分公式如下14个求积节点三角单元高斯积分系数14tablesevengaussintegralnodestriangleelement求积节点个数n求积节点坐标代数精度12对称23采用底层电阻率均匀的边界条件e3的求法与te模式的式14类似视电阻率的计算对上面的单元积分矩阵进行扩展叠加变分形成线性方程组加入边界条件解线性方程组等常规处理后可以得到各个节点的场值13的方法利用主场值的导数求辅助场从而求得复数视电阻率
式 （ 4 ） 中第 4 项积分为 1 · · － iωμσ ·u dl = ∫1 2 iωμ 槡
2
其中： σ i ， σj ， σ m 为三角形单元顶点电导率值； N i ，
1 T u k u 。 （ 13） 2 l e4 l
积分在计算区域的 CD 边进行， 采用 CD 边以下电导 率均匀的边界条件， 则有限单元法系数矩阵为 k e4 = 2 l 1 · · 槡－ iωμσ 6 iωμ 1
目前大地电磁有限单元法正演基本上都是单 元内电导率不变
［ 1 － 2］
但是在 TM 模式下电导率不是线性变化，而是电导 率和介电常数的组合参数线性变化 ； 此外， 文献 本文采用三角单元剖分， 单元内电导率和场 值线性变化，推导出大地电磁 TE 模式有限单元法 系数矩阵，利用高斯积分计算 TM 模式下有限单元 法系数矩阵。 三角单元利用矩形单元的对角线进 行剖分，给出了利用矩形网格变形技术模拟具有 垂直间断 面 的 二 维 地 形 的 方 法， 实 现 了 TE （ 横 二维复杂地形有限单元法正演模拟 。
［ 6］ ［ 5］
、 电） 、TM （ 横磁 ） 模式电导率线性变化大地电磁
利用三角单元进行了大地电磁地
［7 － 8 ］
形模拟，徐世浙等
利用边界单元法和有限单
［ 9］
元法进行了大地电磁地形校正研究， 王绪本等
1
二维模型网格剖分
如图 1a 所示，计算区域采用矩形， 只要在计
利用有限单元法进行了大地电磁二维地形校正研 究，刘云进行了平行四边形二维有限单元法正演 模拟
复杂地形电导率线性变化二维大地 电磁有限单元法正演模拟
欧东新，韦者良
（ 桂林理工大学 a. 广西矿冶与环境科学实验中心； b. 广西隐伏金属矿产勘查重点实验室，广西 桂林 541004 ）
摘
要： 为模拟复杂地形及地质结构，利用三角单元进行剖分，给出了具有垂直断面地形的简单网格剖分
方法； 电导率和电磁场值在单元内线性变化，推导出大地电磁 TE （ 横电） 模式有限单元法系数矩阵，利 用高斯积分计算 TM （ 横磁） 模式有限单元法系数矩阵； 计算辅助场时，在垂向采用二次拉格朗日插值计 算主场值的垂向偏导数； 用有限单元法对一维及二维模型进行视电阻率及相位计算，与前人计算结果进行 了对比检验，计算结果与解析解的相对均方差小于 0. 17% ； 计算了垂直间断面地形的视电阻率断面图。 关键词： 大地电磁； 二维； 垂直断面地形； 电导率线性变化； 有限单元法； 正演 中图分类号： P631 文献标志码： A
［10 ］
算区域内形成网格单元就可以用于有限单元法的 正演计算。图 1b 为用本文网格剖分方法给出的具 有垂直间断 面 地 形 的 三 角 网 格， 粗 线 为 地 形 线。
［12］
，可以模拟连续起伏地形， 他们的模拟方
［ 11 ］
法中电导率都是是分块均匀的； 随后刘云等
给
出了电性参数线性变化的 有 限 单 元 法 正 演 方 法， 三角网格的形成可以采用 Delaunary 算法
+
3
有限单元法系数矩阵
TE 模式系数矩阵 TE 模式离散的积分形式为
在空气和地面的节点上电导率是有突变的， 因此 3. 1
F（ u） =
∑ {
Ω
△
1 · （ u） [1 2 iωμ
2
－
1 （ σ － iωε） u2 dxdy 2
]
}
∑∫
CD
l
1 1 · · － iωμσ ·u2 dl， 2 iωμ 槡
（ 4）
32
桂
林
理
工
大
学
学
报
2014 年
其中： l 为三角形单元位于边界 CD 上的边。 在三角单元内利用面积坐标对电导率及场值 进行线性插值
［ 13］
k e3
， 以电导率为例的插值公式 （ 5）
l
1 / 12 = 2 △·iωε 1 / 24 1 / 24
对称 1 / 12 1 / 24
， 而场值根据不同的插值函
10 － 11］ 中没有给出垂直地形的模拟方法 。 数变化。但是野外的地电分布一般是连续变化的， ［ 因此有必要研究电性参数 连 续 变 化 的 正 演 方 法。 徐世浙、 李予国等
［3 － 4 ］
采用矩形网格剖分， 进行
电性参数分块连续变化、 双线性插值、 双二次插 值等有限元数值模拟， 但是采用矩形单元， 无法 准确模拟起伏地形。 由于地形起伏会引起严重的 静态效应，给资料解释造成极大困难， 因此需要 进行带 地 形 的 大 地 电 磁 正 演 计 算。 胡 建 德 等 Wannamaker 等
。 1 / 12
（ 12 ）
σ = σi Ni + σj Nj + σm Nm ， Nj ， N m 为插值系数， 具体形式如下 Ni = 1 （ ai x + b i y + ci ） ， 2△ 1 （ a x + b j y + cj ） ， Nj = 2 △ j 1 N m = 2 △（ a m x + b m y + cm ） ； aj = ym － yi ， am = yi － yj ； ai = yj － ym ， bj = xi － xm ， bm = xj － xi ； bi = xm － xj ， c = x y － x y ， j m m j cj = x m y i － x i y m ， i cm = x i y j － x j y i ； 1 △ = 2 （ ai bj － aj bi ） 。 角形面积。
第1 期
欧东新等： 复杂地形电导率线性变化二维大地电磁有限单元法正演模拟
31
图1
Fig. 1
计算区域 （ a） 及三角网格剖分 （ b）
Computation region （ a） and triangle mesh （ b）
图2
Fig. 2
垂直间断面地形的网格剖分示意图
Grids division of discontinuity topography