大地电磁交错采样有限差分二维正反演研究

图1 MT二维交错网格有限差分离散示意

Fig.1 MT 2-D staggered grid finite difference separation diagram

1.1 TE模式正演

对于TE模式,以E_y(i,k)为中心,式(1)Maxwell第1个方程展开为:

(3)

\begin{array}{l} H_{z} (i, k) Δ z (k) + H_{x} (i, k + 1) Δ x (i) - \\ H_{z} (i + 1, k) Δ z (k) - H_{x} (i, k) Δ x (i) = \\ σ (i, k) E_{y} (i, k) Δ z (k) Δ x (i) 。 \end{array}

对式(3)中的4个磁场,由式(2)又可分别用周围的电场求得:

(4)

E_{y} (i, k) - E_{y} (i - 1, k) = i ω μ_{0} H_{z} (i, k) \frac{Δ x (i) + Δ x (i - 1)}{2},

(5)

E_{y} (i, k - 1) - E_{y} (i, k) = i ω μ_{0} H_{x} (i, k) \frac{Δ z (k) + Δ z (k - 1)}{2},

(6)

E_{y} (i + 1, k) - E_{y} (i, k) = i ω μ_{0} H_{z} (i + 1, k) \frac{Δ x (i) + Δ x (i + 1)}{2},

(7)

\begin{array}{l} E_{y} (i, k) - E_{y} (i, k + 1) = \\ i ω μ_{0} H_{x} (i, k + 1) \frac{Δ z (k) + Δ z (k + 1)}{2} 。 \end{array}

将这4个关系式代入式(3)中即可消去磁场分量,最终得到E_y(i,k)与其周边4个电场的关系式。用C₁~C₅的系数表示为

(8)

\begin{array}{l} C_{1} E_{y} (i - 1, k) + C_{2} E_{y} (i, k) + C_{3} E_{y} (i + 1, k) + \\ C_{4} E_{y} (i, k - 1) + C_{5} E_{y} (i, k + 1) = 0, \end{array}

其中:

$\begin{array}{l} C_{1} = - \frac{2 Δ z (k)}{Δ x (i) + Δ x (i - 1)}, \\ C_{2} = - C_{1} - C_{5} - C_{4} - C_{3} - σ (i, k) i ω μ_{0} Δ z (k) Δ x (i), \\ C_{3} = - \frac{2 Δ z (k)}{Δ x (i) + Δ x (i + 1)}, \\ C_{4} = - \frac{2 Δ x (i)}{Δ z (k) + Δ z (k - 1)}, \\ C_{5} = - \frac{2 Δ x (i)}{Δ z (k) + Δ z (k + 1)} 。 \end{array}$

1.2 TM模式正演

对于TM模式,以H_y(i,k)为中心,式(2)Maxwell第2个方程展开为

(9)

\begin{array}{l} E_{z} (i, k) Δ z (k) + E_{x} (i, k + 1) Δ x (i) - \\ E_{z} (i + 1, k) Δ z (k) - E_{x} (i, k) Δ x (i) = \\ i ω μ_{0} H_{y} (i, k) Δ z (k) Δ x (i) 。 \end{array}

对式(9)中的4个电场,由式(1)又可分别用周围的磁场求得:

(10)

H_{y} (i, k) - H_{y} (i - 1, k) = \frac{σ (i, k) Δ x (i) + σ (i - 1, k) Δ x (i - 1)}{2} E_{z} (i, k),

(11)

H_{y} (i, k - 1) - H_{y} (i, k) = \frac{σ (i, k) Δ z (k) + σ (i, k - 1) Δ z (k - 1)}{2} E_{x} (i, k),

(12)

H_{y} (i + 1, k) - H_{y} (i, k) = \frac{σ (i, k) Δ x (i) + σ (i + 1, k) Δ x (i + 1)}{2} E_{z} (i + 1, k),

(13)

H_{y} (i, k) - H_{y} (i, k + 1) = \frac{σ (i, k) Δ z (k) + σ (i, k + 1) Δ z (k + 1)}{2} E_{x} (i, k + 1) 。

将这4个关系式代入式(9)中即可消去电场分量,最终得到H_y(i,k)与其周边4个磁场的关系式。用C₁~C₅的系数表示为:

(14)

\begin{array}{l} C_{1} H_{y} (i - 1, k) + C_{2} H_{y} (i, k) + C_{3} H_{y} (i + 1, k) + \\ C_{4} H_{y} (i, k - 1) + C_{5} H_{y} (i, k + 1) = 0, \end{array}

其中:

$\begin{array}{l} C_{1} = - \frac{2 Δ z (k)}{σ (i, k) Δ x (i) + σ (i - 1, k) Δ x (i - 1)}, \\ C_{2} = - C_{1} - C_{5} - C_{3} - C_{4} - i ω μ_{0} Δ z (k) Δ x (i), \\ C_{3} = - \frac{2 Δ z (k)}{σ (i, k) Δ x (i) + σ (i + 1, k) Δ x (i + 1)}, \\ C_{4} = - \frac{2 Δ x (i)}{σ (i, k) Δ z (k) + σ (i, k - 1) Δ z (k - 1)}, \\ C_{5} = - \frac{2 Δ x (i)}{σ (i, k) Δ z (k) + σ (i, k + 1) Δ z (k + 1)} 。 \end{array}$

1.3 方程组

根据前述TE和TM离散方案,最终可以获得关于E_y和H_y的离散方程:

(15)

\begin{array}{l} A_{e} E_{y} = b_{e}, \\ A_{h} H_{y} = b_{h} 。 \end{array}

式中:E_y和H_y为求解向量,b_e、b_h对应的右端项向量,A_e、A_h分别为式(8)、式(14)中不同单元相关的C系数组成。

要想求解研究区域内所有采样点上的场值,还需给出研究区边界上场值,交错采样网格二维边界条件可简要概括为表1所示形式^[2],其中顶边界为第一类边界条件(其中TE模式顶面位于空气顶层,TM模式顶面位于地表),侧边界为第二类边界条件,底边界对应第三类边界条件(其中σ为底层介质电导率,Z_bottom为底层顶面阻抗)。对构建的方程组,由于直接解法可获得精确解,因此这里采用基于LU分解的直接法求解。

表1 边界条件形式

Table 1 Forms of boundary conditions

模式位置	顶边界 (z=z_min)	底边界 (z=z_max)	侧边界 (y=y_min\|y_max)
TE模式	E_y=1	$\frac{\partial E_{y}}{\partial z}$ =iωμ $\frac{E_{y}}{Z_{bottom}}$	$\frac{\partial E_{y}}{\partial x}$ =0
TM模式	H_y=1	$\frac{\partial H_{y}}{\partial z}$ =σ $\frac{H_{y}}{Z_{bottom}}$	$\frac{\partial H_{y}}{\partial x}$ =0

新窗口打开| 下载CSV

2 LBFGS反演

根据Tikhonov的正则化反演理论^[14],反演问题可表述为

(16)

\begin{array}{l} ψ (m) = ψ_{d} (m) + λ ψ_{m} (m) = \\ (W_{d} {Δ d)}^{T} (W_{d} Δ d) + λ (W_{m} {m)}^{T} (W_{m} m) 。 \end{array}

式中:m为模型;ψ_d(m)为数据目标函数;ψ_m(m)为模型约束目标函数;λ为正则化因子;Δd为观测数据与模型响应的差向量;W_d为数据加权矩阵;W_m为模型约束矩阵,目的是使反演得到的模型具备所设定的相应结构特征,一般采用二阶拉普拉斯算子。

在确立反演目标函数后,接下来的工作就是选择合适的方法进行最优化反演。大地电磁反演方法较多,维性、反演效率以及对计算资源需求等因素影响着反演方法的选择。经典的最优化方法基本都是基于梯度下降法和牛顿法的改进或重构,在大地电磁中几乎所有的最优化反演方法都是这两种方法的实用化或变种,其中有限内存拟牛顿法(LBFGS)因其特殊的优势,近年来逐渐应用广泛并受到认可^{[12,13,14,15]}。LBFGS算法流程如下:

1)设定k=0,初始模型m₀和精度阈值ε;

2)设置初始的正定矩阵H₀=I;

3)计算搜索方向d_k=- $H_{k}^{- 1}$ g_k;

4)令初始步长γ_k为1,采用公式γ_k=arg $\min_{γ \in R}$ f(m_k-γ $H_{k}^{- 1}$ g_k)线搜索得到最佳的步长γ_k,再计算m_k₊₁=m_k-γ_k $H_{k}^{- 1}$ g_k;

5)若‖g_k₊₁‖<ε,则算法结束;

6)利用m个已保存的s、y向量对计算H_k₊₁:

$\begin{array}{l} H_{k + 1} = ({V^{T}}_{k} \dots {V^{T}}_{k - m}) H_{0} (V_{k - 1} \dots V_{k}) + \\ ρ_{k - m} ({V^{T}}_{k} \dots {V^{T}}_{k - m + 1}) s_{k - m} {s^{T}}_{k - m} (V_{k - m + 1} \dots V_{k}) + \\ ρ_{k - m + 1} ({V^{T}}_{k} \dots {V^{T}}_{k - m + 2}) s_{k - m + 1} {s^{T}}_{k - m + 1} (V_{k - m + 2} \dots V_{k}) + \\ \dots + ρ_{k} s_{k} {s^{T}}_{k} 。 \end{array}$

7)迭代次数k=k+1,跳转到步骤3继续。

3 算例分析

3.1 层状模型正演

由于层状介质模型可通过一维正演求得其解析解,因此二维正演算法的数值解可采用一维解析解进行对比验证。层状模型第一层电阻率100 Ω·m,厚1 km,第二层电阻率1 000 Ω·m,厚9 km,第三层电阻率10 Ω·m。

计算过程中非扩展区域网格数为横向50列,每列宽度10 km,纵向80层,层深度对数递增,网格剖分单元数4 000个,测点均匀分布在地表每列单元中心共100个,计算了10⁴~10^-4 Hz共61个频率的视电阻率和阻抗相位响应。图2给出了解析解与剖面中心位置测点的视电阻率和阻抗相位对比曲线,可以看出数值结果和解析解曲线基本重合,仅在约1 000 Hz以上高频部分存在细小差异,但这是有限差分数值计算普遍均存在的现象。二维数值解与解析解视电阻率平均相对误差为1.19%,相位平均相对误差为0.72%,表明相对误差较小,二维正演结果精度较高。

图2

图2 层状模型二维正演中心测点视电阻率和阻抗相位

Fig.2 Apparent resistivity and impedance phase of 2-D forward modeling center of layered model

3.2 异常体模型正反演

为进一步测试正演算法的正确性以及反演算法的稳定性和效率,设计了较为复杂的异常体模型。模型的设计参考借鉴Egbert对Modem算法的测试模型^[12],设计由浅至深为2层,横向上3列,共6个交错分布的高低阻异常体,模型共80层×50列=4 000个单元,50个地表测点按100 等距分布,图3为模型的具体形态。模型背景电阻率100 Ω·m,纵向上在150~300 m、800~1 500 m深度存在异常体,位置交错分布,电阻率值分别为1 000 Ω·m、10 Ω·m。计算了10⁴~1 Hz共30个频点的正演响应,图4为正演响应拟断面,可以看到,无论视电阻率还是阻抗相位,都对高阻和低阻异常有较明显的反映。

图3

图3 异常体模型

Fig.3 Anomalous body model

图4

图4 异常体模型正演拟断面

Fig.4 Forward section of abnormal body model

对异常体模型进行反演计算。反演前在正演数据中加入3%的高斯随机噪音,并将其作为数据协方差,反演正则化因子设为 0.01,初始模型为100 Ω·m的均匀半空间。反演终止条件为拟合差RMS不大于1,这里的拟合差计算公式为:

(17)

RMS = \sqrt[]{\frac{\overset{N_{数据总数}}{\sum_{i = 1}} [(d_{i}^{观测值} - d_{i}^{理论值}) / ε_{i}^{观测值误差}]^{2}}{N_{数据总数}}} 。

分别对合成数据开展TE、TM和TE+TM反演,最终反演的步长曲线和拟合差曲线如图5所示。其中,TE模式反演步长经过8次后稳定为1,TM模式9次后稳定,TE+TM模式29次迭代后稳定为1。3种模式反演的拟合差RMS从初始的约11开始,经过几次迭代后快速衰减到3左右,随后缓慢降低,其中TE、TM和TE+TM分别在第17、23、40次迭代后满足收敛条件。从反演步长曲线和拟合差曲线看出反演较为稳定,符合LBFGS收敛特征。

图5

图5 异常体模型二维反演步长和拟合差曲线

Fig.5 2-D inversion step length and fitting difference curve of abnormal body model

TE、TM及TE+TM3种模式LBFGS最终反演结果如图6所示。与真实模型对比看出,3种模式的反演结果总体上与真实模型都较为吻合:TE模式对异常体厚度反映较好,但横向边界相对一般;TM模式对异常横向范围反映较好,但纵向分辨相对一般;TE+TM模式综合了TE和TM的优点,对异常体反映较好。上述特征符合MT二维反演规律,表明算法可靠。

图6

图6 异常体模型二维反演

Fig.6 2-D inversion of anomalous body model

另外,为分析LBFGS反演的效率,将该模型同时开展了NLCG二维反演,表2给出了2种算法的时间对比。对于NLCG反演,TE、TM、TE+TM的10次迭代耗时分别为75、57、133 s,而LBFGS的10次迭代耗时分别为59、45、88 s,相对有约27%的效率提升。NLCG拟合差达到1的计算时间分别为175、92、506 s,而LBFGS只用了125、73、371 s,相对有约25%的效率提升。这表明LBFGS反演效率总体上比NLCG高。

表2 LBFGS和NLCG反演效率对比

Table 2 Comparison of inversion efficiency between LBFGS and NLCG

模式	$\frac{迭代 10 次耗时 / s}{NLCGLBFGS}$		$\frac{拟合差 < 1 的耗时 / s}{NLCGLBFGS}$
TE	75	59	175	125
TM	57	45	92	73
TE+TM	133	88	506	371

新窗口打开| 下载CSV

3.3 实测资料反演

将本文的算法应用到了甘肃省陇南市宕昌县官鹅沟地热勘查项目。勘查区布设了1条5 km长、150 m点距的MT测线(图7),数据的野外观测采用张量观测方式,单测点观测时间不低于6 h。勘查区主要出露地层为下三叠统(T₁)、中三叠统(T₂)和第四系(Q),其中T₁岩性以灰岩为主,夹粉砂质板岩、泥岩;T₂中的 $T_{2}^{a}$ 岩性以灰岩为主, $T_{2}^{b}$ 岩性主要为粉砂岩、粉砂质泥岩;第四系沉积物主要为更新统黄土、冰碛砾石层和全新统砂砾石层。

图7

图7 MT测线分布

Fig.7 MT line distribution

图8为采用本文算法对实测MT资料的TM模式反演结果,分析看出,剖面范围内发育多条断层,断层位置和形态与地表地质调查以及氡气测量结果较为一致。此外,地表的诸多断裂在约1 000 m深度汇集,形成F1和F2两个主要断层,发育深度均大于2 000 m,其中F1倾向南西、 F2倾向北东, 2个断裂在深部交汇,断层深部影响带范围较大。F1和F2断层带中间低阻现象较为明显,表明断层带上盘局部岩体破碎带较发育,为地热流体的赋存提供了一定的储集空间,具有一定的导水、导热作用。地热流体沿断层带上升后储存在热储含水层内,在上覆隔水保温盖层的作用下,形成地热资源。通过实测MT资料的反演和解释,可以看出本文算法的反演结果揭露的特征,与地质和相关资料较为一致。

图8

图8 MT反演结果

Fig.8 MT inversion results

4 结论及讨论

本文借鉴三维大地电磁正演模拟中使用的交错采样网格方法,采用交错采样网格开展了大地电磁二维有限差分正演研究,模型计算表明该正演算法具有较高的计算精度。基于交错采样网格大地电磁二维有限差分正演,利用有限内存拟牛顿算法,实现了大地电磁二维快速反演。理论模型计算表明有限内存拟牛顿反演算法具有较好的稳定性,与非线性共轭梯度算法对比发现,LBFGS计算效率至少大于NLCG效率的1.2倍。将反演算法在宕昌县官鹅沟地热勘查项目中开展了大地电磁资料二维反演解释,结果表明本文算法与现有资料和认识较为一致,具有一定的实用价值。

由于算法暂未考虑地形以及各向异性介质的影响,建议下一步开展基于交错采样网格的起伏地形、各向异性正反演研究。此外,尽管LBFGS具有较高的计算效率,但随着数据量和模型单元的增大,并行计算研究仍然值得开展。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

胡祥云, 霍光谱, 高锐, 等.

大地电磁各向异性二维模拟及实例分析

[J]. 地球物理学报, 2013,56(12):4268-4277.

Hu X

, Huo G

, Gao

, et al.

Two dimensional simulation and case analysis of magnetotelluric anisotropy

[J]. Journal of Geophysics, 2013,56(12):4268-4277.

[2]

童孝忠, 吴思洋, 谢维.

基于有限差分法的大地电磁二维倾子响应计算

[J]. 工程地球物理学报, 2018,15(4):93-99.

Tong X

, Wu S

, Xie

Calculation of magnetotelluric two-dimensional tilter response based on finite difference method

[J]. Journal of Engineering Geophysics, 2018,15(4):93-99.

[3]

张弘宇, 温建亮, 张富明.

有限差分法二维地电构造模型大地电磁法正演模拟

[J]. 矿产勘查, 2019,10(8):1988-1992.

Zhang H

, Wen J

, Zhang F

Forward modeling of two-dimensional geoelectric structural model by finite difference method

[J]. Mineral Exploration, 2019,10(8):1988-1992.

[4]

Mitsuhata

, Uchida

3D magnetotelluric modeling using the T-Ω finite-element method

[J]. Geophysics, 2004,69(1):108-119.

[5]

谭捍东, 余钦范, Booker

, 等.

大地电磁法三维交错采样有限差分数值模拟

[J]. 地球物理学报, 2003,46(5):705-711.

URL [本文引用: 1]

系统地论述了大地电磁三维交错采样有限差分数值模拟算法实现过程中交错网格剖分、积分公式离散化、边界条件、方程组求解、三维张量阻抗的计算等内容. 由于提出了简洁的边界条件，采用了解大型系数矩阵方程组的双共轭梯度稳定解法，所实现的三维交错采样有限差分数值模拟算法具有迭代收敛稳定、计算精度高、速度快等特点. 通过两个理论模型的计算结果检验了算法的正确性和计算精度. 所实现的三维交错采样有限差分数值模拟算法为研究三维反演问题奠定了基础.

Tan H

, Yu Q

, Booker

, et al.

Three dimensional staggered sampling finite difference numerical simulation of magnetotelluric method

[J]. Chinese Journal of Geophysics, 2003,46(5):705-711.

[6]

李焱, 胡祥云, 杨文采, 等.

大地电磁三维交错网格有限差分数值模拟的并行计算研究

[J]. 地球物理学报, 2012,55(12):160-167.

, Hu X

, Yang W

, et al.

Research on parallel computing of three-dimensional staggered grid finite difference numerical simulation of magnetotelluric

[J]. Chinese Journal of Geophysics, 2012,55(12):160-167.

[7]

韩波

大地电磁二维拟牛顿反演研究

[D]. 武汉:中国地质大学(武汉), 2012.

Han

Quasi-Newton method for two-dimensional magnetotelluric inversion

[D]. Wuhan:China University of Geosciences(Wuhan), 2012.

[8]

Constable S

, Parker R

, Constable C

Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data

[J]. Geophysics, 1987,52(3):289-300.

[9]

Rodi

, Mackie R

Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion

[J]. Geophysics, 2001,66(1):174-187.

[10]

Dennis J

, Schnabel R

Numerical methods for unconstrained optimization and nonlinear equations

[M]. Upper Saddle River: Prentice-Hall, 1983.

[11]

阮帅

三维大地电磁有限内存拟牛顿反演

[D]. 成都:成都理工大学, 2015.

Ruan

3D magnetotelluric finite memory quasi Newton inversion

[D]. Chengdu: Chengdu University of Technology, 2015.

[12]

王堃鹏

张量CSAMT三维主轴各向异性正反演研究

[D]. 北京:中国地质大学(北京), 2017.

Wang K

Research on the forward and inverse anisotropy of three-dimensional principal axis of tensor CSAMT

[D]. Beijing: China University of Geosciences(Beijing), 2017.

[13]

凡友华, 刘雪峰, 陈晓非.

面波频散反演地下层状结构的拟牛顿法

[J]. 物探与化探, 2006,10(5):83-86.

[14]

Tikhonov A

, Arsenin V

Solutions of Ill-posed problems

[M]. Washington DC: Winston and Sons, 1977.

Fan Y

, Liu X

, Chen X

Quasi Newton method for surface wave dispersion inversion of underground layered structure

[J]. Geophysical and Geochemical Exploration, 2006,20(5):83-86.