Modeling the 3D Terrain Effect on MT by the Boundary Element Method
-
摘要: 提出了一种用边界元法计算大地电磁场三维地形影响的数值模拟方法.首先用矢量积分理论和电磁场边界条件, 将上半空间(空气)和下半空间(地下介质)两个区域电磁场边值问题变为仅对地形界面的两个矢量面积分方程, 其中一个计算磁场, 称磁场方程; 另一个计算电场, 称电场方程.然后将对地形界面的积分剖分为一系列的三角单元积分.在三角单元积分中, 假设单元中电磁场为水平均匀大地空间电磁场与地形影响的迭加, 并假设地形影响为常项, 这样既保证了计算精度又使得计算方法简便.通过分解和计算, 每一个矢量面积分方程分解为对应3个坐标方向的3个常量线性方程, 这些线性方程组成了对角占优的线性方程组, 可用SSOR方法求解.文中给出了2个三维地形上大地电磁视电阻率曲线的计算结果.Abstract: A numerical method is put forward in this paper, using the boundary element method (BEM) to model 3D terrain effects on magnetotelluric (MT) surveys. Using vector integral theory and electromagnetic field boundary conditions, the boundary problem of two electromagnetic fields in the upper half space (air) and lower half space (earth medium) was transformed into two vector integral equations just related to the topography: one magnetic equation for computing the magnetic field and the other electrical equation for computing the electrical field. The topography integral is decomposed into a series of integrals in a triangle element. For the integral in a triangle element, we suppose that the electromagnetic field in it is the stack of the electromagnetic field in the homogeneous earth and the topography response which is a constant; so the computation becomes simple, convenient and highly accurate. By decomposition and computation, each vector integral equation can be calculated by solving three linear equations that are related to the three Cartesian directions. The matrix of these linear equations is diagonally dominant and can be solved using the Symmetric Successive Over-Relaxation (SSOR) method. The apparent resistivity curve of MT on two 3D terrains calculated by BEM is shown in this paper.
-
Key words:
- 3D terrain /
- MT /
- boundary element method /
- numerical modeling
-
大地电磁场地形影响数值模拟国内外已有不少研究, 如Chouteau and Bouchard (1988)和Wannamaker et al. (1986)用有限元法模拟了二维地形对大地电磁测深信号的影响; Xu et al. (1997b)用边界元法模拟了二维地形影响; 后来, 徐世浙(1995)、徐世浙等(1997a)、阮百尧和王有学(2005)又用边界元法模拟了大地电磁测深的三维地形影响.
本文在徐世浙、阮百尧和周辉等人研究基础上, 对原来需要在地上空气介质和地下导电介质2个区域的边界面积分通过迭加方法变为只对地形界面的面积分, 对原来只考虑磁场方程变为同时考虑磁场方程和电场方程, 对原来三角单元电磁场使用零次插值改为仅地形影响项采用零次插值.由于这些改进, 不仅使方法的计算精度得到了提高, 而且使最后形成的线性方程组具有明显的对角占优特征, 可用SSOR方法求解, 因此与原方法相比大大提高了计算速度.
1. 电磁场边界积分方程
MT中的初始波是垂直入射的平面波, 设初始磁场H的偏振方向沿x轴(图 1), 时间因子为e-iwt.选取足够大的六面体区域, 期间被三维地形Γ所分割, 地形上部记为Ω1, 是空气; 地形下部记为Ω2, 是充满电导率为σ的均匀岩石.设地形起伏局限于区域的中部, 由于区域足够大, 地形在六面体边界Γ1和Γ2上产生的异常电磁场为零.若取AB处的磁场Hx为单位磁场, 可得到磁场和电场的边界条件如下:
图 1 区域与边界Ω1.空气; Ω2.地下均匀岩石; Γ.三维地形; Γ1和Γ2.分别表示上下部无穷远边界; n.表示法向Fig. 1. Domain and boundary地形以上边界Γ1
地形以下边界Γ2
其中
设地形以上Ω1区域中电磁场场强为E1、H1, 从Maxwell方程经过矢量分析, 可推算得如下积分方程(徐世浙等, 1997a) :
(3) 式中p1为区域Ω1中一点; ωp1是p1点对区域Ω1所张的立体角, 当p1点在区域Ω1内时ωp1=4π, Γ为地形边界; n1垂直地面向下;
(4) 上式是边界影响项, 由于边界Γ1上场值已知式(1), 因此E1b(p1) 和H1b (p1)是已知项, 可以通过高斯数值积分求得.
同理, 若地形以下Ω2区域场强为E2、H2, 有
(5) 式中p2为区域Ω2中一点, ωp2是p2点对区域Ω2所张的立体角, Γ为地形边界, n2垂直地面向上,
(6) 同样, 上式是已知的边界影响项, 可以通过高斯数值积分求得.
根据界面上电场和磁场切向分量连续、电流密度法向分量连续, 当p1、p2点位于地面上p点时, 有
(7) 
(8) 
(9) 式中:
用边界条件(7)、(8)、(9) 将(5) 式中的电磁场E2, H2变换成E1, H1, 然后与(3) 式相加, 就可得到下面的三维地形问题电磁场边界积分方程:
(10) 式中, 由于电场与磁场下标都为1, 即统一到上半空间区域的电场和磁场, 所以为了简化, 去掉电场和磁场下标1.其中En和Hn为电场和磁场的法向分量, np为地形界面上p点的法向矢量.通常(10) 式第一个方程称为电场方程, 第二个方程称为磁场方程.解上述边界积分方程, 即可求得三维地形的大地电磁响应.
2. 边界元方法
用三角形将三维地形剖分成M个单元, 如图 2所示, 并取单元的中心为节点, 这时, 节点上的ωp=2π, 将(10) 式中对地面Γ的积分分解为各三角单元Γj积分之和, 这样对p点(假设为第i个节点) (10) 式可写成:
(11) 式中基本解φ1和φ2的下标i表示其中的r是从节点i至积分点的距离,
三角单元Γj中任意一点p的电磁场E(p)和H(p) 为水平大地空间上的电磁场Ep(p)和Hp(p)以及地形影响之和.为计算简便, 假设单元内地下空间影响项为常量, 即单元内电场和磁场满足以下关系:
(12) 将(12) 式代入(11) 式, 有
(13) 其中
(14) 上式中基本函数、水平大地上的电场和磁场都是已知的, 所以可用高斯数值积分法求得Eci和Hci.
边界积分方程(13) 是矢量方程, 可分解成x, y, z方向上的6个常量线性方程, 其中3个由电场方程所得, 3个由磁场方程所得.由于每个节点的积分方程(13) 都可得到6个常量线性方程, 则M个单元(节点) 可得到6M个线性方程.这些方程组成明显对角占优的线性方程组, 可用超松弛迭代法(SSOR) 求解.
求得地表各单元节点上的电磁场后, 便可用积分方程(3)、(4) 计算出地表上方空气介质中或地表下方岩石介质中任意一点的电磁场.
3. 视电阻率的计算
Hx极化时MT视电阻率公式为:
(15) 4. 算例
考虑到地形起伏处电磁场变化较大, 设计采用以地形起伏中心为中心点, 先进行x方向和y方向的直线剖分, 再用斜线以围绕中心的方向将矩形单元剖分为三角单元, 如图 2所示, 然后读入节点上的高程.地下介质电阻率取100 Ω·m.MT中的初始磁场H的偏振方向沿x轴, T取0.1 s.
算例1, 均匀水平半空间大地.由于本文所设边界条件的限制, 本文所给算法不能计算二维地形上的视电阻率曲线.因此, 为了观测本文方法的计算精度, 我们用对均匀水平空间大地的视电阻率计算作为检验标准.取x方向和y方向节点数均为19, 取Δx和Δy为等距: 5 m.计算结果见表 1.
表 1 水平均匀大地上MT边界元法计算结果(T=0.1 s)Table Supplementary Table Result obtained by the MT boundary element method in the horizontal homogeneous half-space
算例2, 三维正地形, 平面图见图 3.纵剖面(y=0.33 m) 视电阻率曲线见图 4a, 横剖面(x=0.33 m) 见图 4b.从图中可见, 三维正地形影响, 在地形隆起处, 视电阻率值变小.
图 4 三维正地形上MT视电阻率曲线(T=0.1s)Fig. 4. MT apparent resistivit y plot in t he uplifted t opography算例3, 三维负地形, 平面图与算例2相同, 见图 3.从计算结果(图 5) 可见, 由于三维负地形的影响, 在地形凹陷处, 视电阻率变大.
图 5 三维负地形上MT视电阻率曲线(T=0.1 s)Fig. 5. MT apparent resistivity plot in the 3D depressed topography5. 结论
本文介绍的三维地形大地电磁场的边界元计算方法, 通过(1) 利用电磁场边界条件, 将上半空间(空气) 和下半空间(地下介质) 两个区域电磁场边值问题变为仅对地形界面的两个矢量面积分方程; (2) 在边界元积分中, 假设单元中电磁场为水平均匀大地电磁场与地形影响的迭加, 并假设地形影响为常项等手段, 使计算方法显得简便快速, 并保证了计算精度.
-
图 1 区域与边界
Ω1.空气; Ω2.地下均匀岩石; Γ.三维地形; Γ1和Γ2.分别表示上下部无穷远边界; n.表示法向
Fig. 1. Domain and boundary
图 4 三维正地形上MT视电阻率曲线(T=0.1s)
Fig. 4. MT apparent resistivit y plot in t he uplifted t opography
图 5 三维负地形上MT视电阻率曲线(T=0.1 s)
Fig. 5. MT apparent resistivity plot in the 3D depressed topography
表 1 水平均匀大地上MT边界元法计算结果(T=0.1 s)
Table 1. Result obtained by the MT boundary element method in the horizontal homogeneous half-space
-
[1] Chouteau, M., Bouchard, K., 1988. Two-di mensional terrain correction in magnetotelluric surveys. Geophysics, 53: 854-862. doi: 10.1190/1.1442520 [2] Ruan, B. Y., Wang, Y. X., 2005. A boundary element modeling method for the electromagnetic field by artificial source in frequency domain with 3-D topography. Chinese J. Geophy. , 48 (5): 1197-1204 (in Chinese with English abstract). [3] Wannamaker, P. E., Stodt, J. A., Rijo, L., 1986. Two-di mensional topographic responses in magnetotelluric model using finite elements. Geophysics, 51: 2131-2144. doi: 10.1190/1.1442065 [4] Xu, S. Z., 1995. The boundary element method in geophysics. Science Press, Beijing (in Chinese). [5] Xu, S. Z., Ruan, B. Y., Zhou, H., et al., 1997a. Numerical modeling of 3-D terrain effect on MT field. Science in China (Series D), 43 (2): 269-275 (in Chinese). [6] Xu, S. Z., Ruan, B. R., Zhou, H., et al., 1997b. Modeling the 2D terrain effect on MT by the boundary element method. Geophysical Prospecting, 45: 931-943. doi: 10.1046/j.1365-2478.1997.610301.x [7] 阮百尧, 王有学, 2005. 三维地形频率域人工源电磁场的边界元模拟方法. 地球物理学报, 48 (5): 1197-1204. doi: 10.3321/j.issn:0001-5733.2005.05.031 [8] 徐世浙, 1995. 地球物理中的边界元法. 北京: 科学出版社. [9] 徐世浙, 阮百尧, 周辉, 等, 1997a. 大地电磁场三维地形影响的数值模拟. 中国科学(D辑), 43 (2): 269-275. https://www.cnki.com.cn/Article/CJFDTOTAL-JDXK199701002.htm -
下载:
下载:
点击查看大图
计量
- 文章访问数: 3387
- HTML全文浏览量: 89
- PDF下载量: 4
- 被引次数: 0
