DOI: 10.11817/j.issn.1672-7207.2018.03.016
含地形起伏的MT模型非规则四边形FEM正演与反演
冯德山1, 2，刘金宝1, 2，王珣1, 2
(1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；
2．有色资源与地质灾害探查湖南省重点实验室，湖南 长沙，410083)
摘要：从大地电磁(MT)边值问题满足的变分方程出发，采用非规则四边形网格、双线性插值有限单元法(FEM)开展复杂起伏地形MT模型正演，探讨二维Jacobian变换行列式的计算方法，推导任意非规则四边形单元的插值方式及单元系数矩阵表达式，实现起伏地形MT模型的高精度正演。然后，介绍光滑约束的Tikhonov正则化反演算法，针对反演中正则化参数选取困难的问题，将广泛应用的L-curve法引入反演的正则化参数选取中。研究结果表明：L-curve法的曲线中曲率最大的拐点准确地指示了最优正则化参数；L-curve法选取的最优正则化参数对应的反演结果与原模型所示结果吻合度最高，反演效果最好。
关键词：大地电磁；非规则四边形网格；L-curve法；正则化约束反演
中图分类号：P631 文献标志码：A 文章编号：1672-7207(2018)03-0626-07
MT finite element method forward modeling and inversion using irregular quadrilateral mesh of steep topography model
FENG Deshan1, 2, LIU Jinbao1, 2, WANG Xun1, 2
(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
2. Key Laboratory of Non-ferrous Resources and Geological Detection, Ministry of Hunan Province, Changsha 410083, China)
Abstract: Based on the variational problem derived from magnetotelluric(MT) boundary value problems, the finite element method using irregular quadrilateral mesh and bilinear interpolation was used to solve MT forward problem of steep topography model. The detailed calculation of Jacobi transformation matrix was discussed and interpolation method of arbitrary irregular quadrilateral unit and unit coefficient matrix expression was derived which achieved high precision forward of steep topography MT model. Then, the smooth constrained Tikhonov regularization inversion algorithm was introduced. To overcome the difficulty in determination of regularization parameters for MT, L-curve method was used to calculate regularization parameters in MT inversion. The results show that inflection of L-curve can indicate the optimal regularization parameter accurately. Inversion results using optimal regularization parameter calculated by L-curve method is most similar to those of the original model, and better inversion results can be obtained.
Key words: magnetotelluric(MT); irregular quadrilateral mesh; L-curve method; smoothness constrained inversion
大地电磁是以天然交变电磁场为场源，通过在地表观测相互正交的电场、磁场分量来获取地电构造信息，并以岩矿石的电性差异为基础的一种重要的地球物理勘探方法[1]。由于探测对象电性结构的复杂性以及地形起伏等影响因素，导致MT探测环境复杂。在地形起伏较大的复杂山区开展MT勘探时，起伏地形对MT探测数据影响很大，若不考虑起伏地形引起的大地电磁场畸变，则将使MT探测数据产生显著误差，其解释结果必然偏离实际结果[1-7]。为了正确认识起伏地形对MT探测数据的影响，提高复杂地形条件下MT探测数据的解释精度，需要开展带地形条件下快速、高精度MT正反演。徐世浙[2]将有限单元法(FEM)引入MT的正演计算中；CHOUTEAU 等[3]采用FEM模拟了起伏地形对大地电磁响应的畸变效应，并应用畸变张量法进行了二维起伏地形校正；赵广茂等[4]实现了二次场起伏地形MT响应的FEM数值模拟，使计算结果更符合野外真实情况；此外，刘树才等[5-6]指出在地形起伏的突变处，将会导致局部电磁场突变，并开展了带地形的MT正演。为了进一步提高起伏地形的MT正演精度，柳建新等[7]应用自适应FEM开展地堑、地垒、山谷、山脊及组合地形的MT正演，总结了起伏地形对MT电磁响应曲线的影响规律；赵慧等[8]将完全非结构化三角形自适应FEM算法应用于海洋MT正演；韩骑等[9]采用非结构三角网格自适应FEM开展任意二维地质结构的正演，并将自适应FEM与Occam反演相结合，有效指导复杂地形下MT数据处理与解释。以上研究均表明地形的多变性对MT正演工作有重要影响，在MT正演工作中应全面考虑地形因素。尽管带地形的MT正演能够对MT特征起到一定的指导作用，但为了使MT资料解释能更直接与直观，需要开展带起伏地形的复杂MT模型反演，然而，大多数MT反演是不适定的病态问题，常用的做法是在反演计算中引入正则化方法[10]。在正则化方法应用中，正则化参数的选取比较困难。正则化参数控制着正则化误差项和观测数据误差项在目标函数中所占的比例，它的选取对解的形态起着关键的作用[11-12]。目前，已有多种方法求取最适宜的正则化参数：陈小斌等[13]选取数据拟合泛函与模型稳定泛函的比值作为每次迭代的正则化因子；刘小军等[14]利用共轭梯度法求解MT反问题，在构建目标函数时引入正则化思想，且在每次迭代过程中能根据目标的收敛情况更新正则化因子；韩波[15]开展了拟Newton法MT二维反演，讨论了正则化参数的选取、灵敏度矩阵的计算；向阳等[16]对比了多种正则化因子选取方法在给定共轭梯度法求解MT反问题中的表现，提出了改进的自适应正则化方案；相鹏[17]通过采用正则化参数的自适应更新算法，构建性状更接近Gauss-Newton法Hessian矩阵的预条件算子，提出一种改进的预条件非线性共轭梯度法，兼顾了反演的稳定性和计算效率；尹银轩等[18]利用广义交叉验证方法确定了最优的正则化参数，保持数据目标函数和模型约束目标函数之间的平衡，使模型方差和分辨率之间达到最佳折中，将MT反演的不适定问题转化为适定问题进行求解。在MT的正则化反演中，正则化因子的选取方法是多样的，计算和比较各种正则化因子的选取方法以选取最适当算法将会极大改善反演结果，其中L曲线法结果直观，计算相对简单，并且是一种具有很强适应性的算法。为了兼顾起伏地形条件下MT模型正演计算的效率与精度，本文作者采用任意非规则四边形网格双线性插值的FEM算法进行MT正演。该非规则四边形剖分方式对复杂模型的拟合既具有较高的精度，又保持结构化矩形网格剖分的计算速度。在反演计算中，首先通过L-curve法选取最优正则化参数，然后应用该最优正则化因子开展光滑约束的Tikhonov正则化MT反演，不会出现因正则化参数过小，导致背景正电阻率及异常体都产生形变，或者因正则化参数过大，导致电阻率异常与背景过度平缓、异常面积显得过大等问题，有效地提高了反演结果的准确性与稳定性。
1 非规则四边形网格FEM大地电磁正演
应用FEM开展MT正演时，网格剖分是FEM前处理技术的一项重要研究内容，直接影响着FEM的求解效率与计算精度。常用的网格剖分方式有结构化网格及非结构化网格。非结构化网格易于离散外形复杂的区域，网格单元的大小、形状和网格节点的位置更灵活，但生成的工作量大，离散求解速度也慢，程序编写复杂，实现较困难。普通四边形结构化网格插值相对于三角网格的面积插值可以在2个方向上双线性拉伸，同时具有较快的计算速度，计算简单，但对起伏地形、复杂模型拟合不好，计算精度不高。为此，采用图1所示不规则四边形单元对模型进行剖分。不规则四边形单元能够灵活地拟合复杂起伏地形，并且具备普通四边形网格的优势，在计算速度上接近普通矩形剖分。
图1 起伏地形不规则四边形网格剖分示意图
Fig. 1 Sketch map of irregular quadrilateral grids subdivision on rugged terrain
由二维大地电磁正演的边值问题可以推导出如下变分问题[2]：
；
；
(1)
CD和AB分别表示区域
的下边界和上边界，
为区域外边界。在TE模式中，
，
；在TM模式中，
，
。其中u在2种模式中分别为电场和磁场在走向方向的分量；
为电磁波的角频率；
为介质的介电常数；
为磁导率；
为电导率；
。
使用四边形单元对整个区域进行剖分，采用双线性插值，编号和坐标见图2[2]。
图2 四边形剖分双线性插值母和子单元节点编号示意图
Fig. 2 Sketch maps of node number of parent element and sub-element based on irregular quadrilateral
用
和
表示四边形单元节点i (i=1, 2, 3, 4)的坐标，则双线性插值的形函数可表示为如下形式：
(2)
用ui表示四边形单元上的电场或磁场强度，则电磁场值与子单元上的坐标可以表示为
(3)
单元刚度矩阵的计算中将会用到Jacobian变换行列式，为此，首先讨论二维Jacobian变换行列式的计算，则面元
可变换为
(4)
其中：Jacobian变换行列式
可表示为
(5)
由式(2)可以得出：
(6)
将式(6)代入式(5)可以得出
(7)
其中：
(8)
(9)
，
，
，因此，Jacobian行列式可表示为
和
的线性函数，它只与单元4个顶点的坐标有关。利用上述形函数以及Jacobian行列式的计算进行单元分析，将式(1)中的区域积分分解为各单元积分之和：
(10)
将式(10)中的
，
和
求出即可得各单元积分，将式(2)和式(9)代入可得：
(11)
其中：
(12)
单元矩阵均为4×4的矩阵，大小与每个单元的顶点坐标有关。将所得单元矩阵扩展成全体单元的系数矩阵，再代入上边界条件，可以解出所有单元节点的电场值或磁场值，利用数值方向求出场值沿垂向的偏导数
，在TE和TM模式中分别对应
和
。将其代入式(13)和(14)，即可求得视电阻率
和相位
。
TE模式：
；
(13)
TM模式：
；
(14)
2 光滑度约束Tikhonov正则化反演
2.1 Tikhonov正则化反演
MT反演大多是高度不适定的病态问题，数据微小误差或较小扰动就会导致所求近似解严重偏离真解。为了提高反演结果的稳定性，使反演结果更加真实可靠，通常会在反演过程中引入Tikhonov正则化方法，其目标函数可以表示为
(15)
式中：
为总目标函数；
为数据目标函数，即观测数据与预测数据拟合差的平方和，
(16)
为模型约束目标函数，根据不同的约束要求可以选择不同的构建方法：
(17)
式中：mr为包含先验信息的参考模型；C为加权矩阵；不同的C对应不同类型的模型约束。本文采用差分矩阵来定义模型的粗糙度，并以此来约束模型的光滑程度。基于先验模型的光滑约束最小二乘正则化反演的迭代公式可以写为
(18)
其中：
为模型的第k次的迭代结果；
为待求模型的修正矢量。
(19)
为观测数据与计算数据间的误差矢量；为待求模型的修正矢量；J为正演过程中求得的雅克比矩阵。为实现TM和TE模式的联合反演，可将模型参数m表示为如下形式：
(20)
式中：
为第i层第j个网格单元的电阻率；M为网格层数；N为水平网格数目。对TM和TE模式的数据目标函数，可以对视电阻率和阻抗相位进行同时拟合，其表达式为
(21)
式中：NT为周期个数；Nx为测点数；
为实测视电阻率；
为模型响应视电阻率；
为实测阻抗相位；
为模型响应阻抗相位。
2.2 正则化参数选取L-curve准则
在反演的目标函数中，正则化参数
的选择十分关键，对反演结果有很大影响。当取值过小时，强调数据的拟合而忽略模型的光滑度，会出现反演的模型十分复杂且与实际情况不相符，甚至可能会过度拟合数据中的误差部分；当过大时，强调模型的光滑度而忽略数据的拟合程度，加大了原问题与新问题之间的偏离程度，得到的模型可能会过于简单，数据得不到充分拟合。而L-curve法是由HANSEN等[19-20]提出的正则化参数选取方法，因绘制出的曲线形状呈“L”形而得名。该方法通过绘制正则解的范数
和正则化解残差的范数
的对数曲线，并通过选取曲线曲率最大的拐点来确定最佳正则化参数。曲线中的垂直部分表示数据的拟合状态良好但模型的契合程度不足，此时，对正则化参数的变化很敏感。曲线的水平部分表示过正则化状态，解的范数趋于稳定，而正则化解残差范数对的变化非常敏感。L-curve法选择垂直部分和水平部分的交点来确定最佳正则化参数，在数据拟合度和解的稳定性之间作出最佳平衡。
L曲线法通过求解(
，
)曲线的最大曲率，即L曲线的拐角点，以确定最佳正则化参数。因此，L曲线法的关键是定位L曲线上曲率最大的那个点。以对数形式推导公式，令
，
，两边取对数，得
，
(22)
则L曲线由许多点(
，
)拟合而成。
用
和
分别表示
和
的一阶导数，
和
分别表示和的二阶导数，这里，，，，和均是正则化参数的函数，则L曲线的曲率
的计算公式为
(23)
由于
，
，
，
，
，代入式(23)可得
(24)
通过求式(24)的最大曲率，即可得到所求的正则化参数。L曲线的优点是不需要预先知道误差估计，易于找出极小值点，具有很强的适应性，不受扰动误差的影响。
3 反演计算实例
设计如图3所示的二维起伏地形模型，地表共有2个山脊地形，其最高点高度分别为200 m和400 m，大地介质的视电阻率为100Ω·m。地下存在3个矩形异常体，其宽度与高度均为600 m，左侧2个是电阻率为10Ω·m的低阻体，其顶端埋深分别为1.6 km和2.1 km；最右侧是电阻率为1.0 kΩ·m的高祖异常体，其顶端埋深为1.5 km。模型总水平跨度为16.4 km。地下深度小于4.0 km。
图3 二维起伏地形模型图
Fig. 3 Geoelectric model of 2D undulating terrain
正演基于不规则四边形的FEM方法，地下介质水平网格数为68个，垂直网格数为40个，单元网格的长×宽为100 m×200 m，频率选取10-1~103 Hz范围内的21个频点。TM模式在横向上有更高的分辨率，在纵向上的分辨率较低。TE模式在纵向上有更高的分辨率，在横向的分辨率较低。因此，本文采用TM和TE模式联合反演，使反演结果优于单一模式的反演，具有更高的分辨率[21]。
正则化参数选用4.00，2.00，1.00，0.50，0.25，0.15，0.10和0.05，将这8组参数分别代入Tikhonov正则化目标函数，可以得到8组不同的反演结果。将每次的数据拟合差和模型稳定项代入L-curve法中进行计算，得到图4中所示的L型曲线。图4中曲线最大拐点处为0.25。
图4 L-curve法求得的最优正则化参数
Fig. 4 Comparison of L-curve criterion for computing optimal regularization parameter
选用图5所示的正则化参数λ分别为0.05，0.25和4.00的反演结果。从图5可见：因λ不同导致反演结果存在显著差异，其中图5(c)所示的λ=0.25时反演剖面中异常体的位置和大小均清晰地得到反映，低阻异常体电阻率准确得到体现，但高阻异常体电阻率偏低。这是由于高阻异常体的正演分辨率较低，导致反演对高阻体不敏感。图5(a)所示为正则化参数λ=0.05时的反演剖面，由于正则化参数选取过小，背景平均视电阻率出现波动，且异常体的特征也有所形变。图5(b)所示为正则化参数λ=4.00时的反演剖面，当正则化参数选取过大时，电阻率异常与背景过于平缓，从而导致异常面积过大，模型约束过度导致反演结果不准确。
图5 二维起伏地形模型在不同的正则化参数条件下反演剖面
Fig. 5 Inversion results of 2D rugged terrain model using different regularization parameters
图6所示为不同正则化参数下的反演剖面与原模型剖面之间的误差曲线。按照反演结果的网格将原模型绘制出，计算出所有网格的均方误差，计算公式如下：
(25)
式中：n为模型网格总数；
为原模型每个网格的电阻率；
为反演网格的电阻率。以正则化参数递减的顺序计算每次结果与原模型的误差，可以从图6中看出L-curve法选取的最佳正则化参数的反演误差为曲线最低点，这说明通过L-curve法选取最优正则化参数能有效提高反演质量。
图6 不同正则化参数的模型各网格拟合均方误差对比
Fig. 6 Comparison of model fitting MSE with different regularization parameters
4 结论
1) 不规则四边形网格能够灵活地模拟复杂起伏地形，且计算过程简单，编程难度较非结构网格的低，是一种较合适的FEM网格剖分方式。
2) 开展了复杂起伏地形MT模型的不规则四边形双线性插值FEM正演及光滑度约束的Tikhonov反演。本文采用的反演算法能准确地反映异常体的特征，对复杂地形具有良好的适用性，说明了正反演算法的有效性和可行性。
3) L-curve法通过绘制正则解的范数和正则化解残差范数的双对数曲线，然后选取曲线中曲率最大的拐点来确定最优正则化参数。该最优正则化参数能够有效地平衡反演目标函数中模型约束项和数据拟合项，使正则化MT反演结果与原模型所得结果更相符，显著提高反演质量。
参考文献：
[1] 柳建新, 童孝忠, 郭荣文, 等. 大地电磁测深勘探: 资料处理反演与解释[M]. 北京: 科学出版社, 2012: 1-12.
LIU Jianxin, TONG Xiaozhong, GUO Rongwen, et al. Magnetotelluric sounding exploration: data processing inversion and interpretation[M]. Beijing: Science Press, 2012: 1-12.
[2] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 229-241.
XU Shizhe. Finite element in geophysics[M]. Beijing: Science Press, 1994: 229-241.
[3] CHOUTEAU M, BOUCHARD K. Two-dimensional terrain correction in magnetotelluric surveys[J]. Geophysics, 1988, 53(6): 854-862.
[4] 赵广茂, 李桐林, 王大勇, 等. 基于二次场二维起伏地形MT有限元数值模拟[J]. 吉林大学学报(地球科学版), 2008, 38(6): 1055-1059.
ZHAO Guangmao, LI Tonglin, WANG Dayong, et al. Secondary fied-based two-dimension topographic numerical simulation in magnetotelluric by finite element method[J]. Journal of Jilin University(Earth Science Edition), 2008, 38(6): 1055-1059.
[5] 刘树才, 何昭友, 刘志新. 适合地形起伏的二维有限元数值模拟[J]. 物探化探计算技术, 2005, 27(2): 131-134.
LIU Shucai, HE Zhaoyou, LIU Zhixin. 2D FEM numerical simulation adapted to the fluctuant topography[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2005, 27(2): 131-134.
[6] 刘云, 王绪本. 大地电磁二维自适应地形有限元正演模拟[J]. 地震地质, 2010, 32(3): 382-391.
LIU Yun, WANG Xuben. FEM using adaptive topography in 2D MT forward modeling[J]. Seismology and Geology, 2010, 32(3): 382-391.
[7] 柳建新, 孙丽影, 童孝忠, 等. 基于自适应有限元的起伏地形MT二维正演模拟[J]. 地球物理学进展, 2012, 27(5): 2016-2023.
LIU Jianxin, SUN Liying, TONG Xiaozhong, et al. Undulating terrain 2D MT forward modelling with adaptive finite element[J]. Progress in Geophysics, 2012, 27(5): 2016-2023.
[8] 赵慧, 刘颖, 李予国. 自适应有限元海洋大地电磁场二维正演模拟[J]. 石油地球物理勘探, 2014, 49(3): 578-585.
ZHAO Hui, LIU Ying, LI Yuguo. Adaptive finite element forward modeling for two-dimensional marine magnetotelluric fields[J]. Oil Geophysical Prospecting, 2014, 49(3): 578-585.
[9] 韩骑, 胡祥云, 程正璞, 等. 自适应非结构有限元MT二维起伏地形正反演研究[J]. 地球物理学报, 2015, 58(12): 4675-4684.
HAN Qi, HU Xiangyun, CHENG Zhengpu, et al. A study of two-dimension MT inversion with steep topography using the adaptive unstructured finite element method[J]. Chinese Journal of Geophysics, 2015, 58(12): 4675-4684.
[10] TIKHONOV A N, ARSENINV Y. Solutions of ill-posed problems[M]. New York: Wiley, 1977: 1-30.
[11] MARTIN H. Limitations of the L-curve method in ill-posed problems[J]. Bit Numerical Mathematics, 1996, 36(2): 287-301.
[12] HANSEN P C, JENSEN T K, RODRIGUEZ G. An adaptive pruning algorithm for the discrete L-curve criterion[J]. Journal of Computational and Applied Mathematics, 2007, 198(2): 483-492.
[13] 陈小斌, 赵国泽, 汤吉, 等. 大地电磁自适应正则化反演算法[J]. 地球物理学报, 2005, 48(4): 937-946.
CHEN Xiaobin, ZHAO Guoze, TANG Ji, et al. An adaptive regularized inversion algorithm for magnetotelluric data[J]. Chinese Journal of Geophysics, 2007, 48(8): 937-946.
[14] 刘小军, 王家林, 吴健生. 二维大地电磁正则化共轭梯度法反演算法[J]. 上海地质, 2007, 101(1): 71-74.
LIU Xiaojun, WANG Jialin, WU Jiansheng. Inversion algorithm of 2-D magnetotelluric data using regularized conjugate gradient method[J]. Shanghai Geology, 2007, 101(1): 71-74.
[15] 韩波. 大地电磁二维拟牛顿反演研究[D]. 武汉: 中国地质大学地球物理与空间信息学院, 2012: 46-47.
HAN Bo. Quasi-Newton Method for two-dimension magnetotelluric inversion[D]. Wuhan: China University of Geosciences. Institute of Geophysics & Geomatics, 2012: 46-47.
[16] 向阳, 于鹏, 陈晓, 等. 大地电磁反演中改进的自适应正则化因子选取[J]. 同济大学学报(自然科学版), 2013, 41(9): 1429-1434.
XIANG Yang, YU Peng, CHEN Xiao, et al. An improved adaptive regularized parameter selection in magnetotelluric inversion[J]. Journal of Tongji University(Natural Science), 2013, 41(9): 1429-1434.
[17] 相鹏. 一种改进的二维MT预条件非线性共轭梯度反演方法[J]. 中国石油大学学报(自然科学版), 2014, 38(4): 42-49.
XIANG Peng. Two-dimensional MT inversion method based on an improved preconditioned nonlinear gradient algorithm[J]. Journal of China University of Petroleum(Edition of Natural Science), 2014, 38(4): 42-49.
[18] 尹银轩, 张连伟, 童孝忠. 大地电磁二维正则化反演研究[J]. 铁道工程学报, 2014, 187(4): 19-23.
YIN Yinxuan, ZHANG Lianwei, TONG Xiaozhong. Research on two-dimensional magnetotelluric inversion[J]. Journal of Railway Engineering Society, 2014, 187(4): 19-23.
[19] HANSEN P C. Analysis of discrete ill-posed problems by means of the L-curve[J]. Society for Industrial and Applied Mathematics Review, 1992, 34(4): 561-580.
[20] HANSEN P C, O’LEARY D P. The use of the L-curve in the regularization of discrete ill-posed problems[J]. Society for Industrial and Applied Mathematics, 1993, 14(6): 1487-1503.
[21] 冯德山, 王珣. 大地电磁双二次插值FEM正演及最小二乘正则化联合反演[J]. 中国有色金属学报, 2013, 23(9): 2524-2531.
FENG Deshan, WANG Xun. Magnetotelluric finite element method forward based on biquadratic interpolation and least squares regularization joint inversion[J]. The Chinese Journal of Nonferrous Metals, 2013, 23(9): 2524-2531.
(编辑 陈灿华)
收稿日期：2017-03-10；修回日期：2017-05-16
基金项目(Foundation item)：国家自然科学基金资助项目(41574116)；中南大学创新驱动项目(2015CX008)；中南大学教师研究基金资助项目2014JSJJ001)，中南大学升华育英人才计划项目(2012)；湖湘青年创新创业平台培养对象共同资助项目(2013) (Project(41574116) supported by the National Natural Science Foundation of China; Project(2015CX008) supported by Innovation Driven Item of Central South University; Project(2014JSJJ001) supported by the Research Foundation of Central South University; Project(2012) supported by Shenghua Yuying Project of Central South University; Project(2013) supported by Youth Innovation and Entrepreneurship Training Platform of Hunan Province)
通信作者：冯德山，博士，教授，博士生导师，从事地球物理正反演与数据处理研究；E-mail: fengdeshan@126.com
