含地形起伏的MT模型非规则四边形FEM正演与反演

2018-04-11冯德山刘金宝王珣

中南大学学报（自然科学版） 2018年3期

冯德山，刘金宝，王珣

冯德山1, 2，刘金宝1, 2，王珣1, 2

(1. 中南大学地球科学与信息物理学院，湖南长沙，410083； 2．有色资源与地质灾害探查湖南省重点实验室，湖南长沙，410083)

从大地电磁(MT)边值问题满足的变分方程出发，采用非规则四边形网格、双线性插值有限单元法(FEM)开展复杂起伏地形MT模型正演，探讨二维Jacobian变换行列式的计算方法，推导任意非规则四边形单元的插值方式及单元系数矩阵表达式，实现起伏地形MT模型的高精度正演。然后，介绍光滑约束的Tikhonov正则化反演算法，针对反演中正则化参数选取困难的问题，将广泛应用的L-curve法引入反演的正则化参数选取中。研究结果表明：L-curve法的曲线中曲率最大的拐点准确地指示了最优正则化参数；L-curve法选取的最优正则化参数对应的反演结果与原模型所示结果吻合度最高，反演效果最好。

大地电磁；非规则四边形网格；L-curve法；正则化约束反演

大地电磁是以天然交变电磁场为场源，通过在地表观测相互正交的电场、磁场分量来获取地电构造信息，并以岩矿石的电性差异为基础的一种重要的地球物理勘探方法[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所示不规则四边形单元对模型进行剖分。不规则四边形单元能够灵活地拟合复杂起伏地形，并且具备普通四边形网格的优势，在计算速度上接近普通矩形剖分。

由二维大地电磁正演的边值问题可以推导出如下变分问题[2]：

图1 起伏地形不规则四边形网格剖分示意图

图2 四边形剖分双线性插值母和子单元节点编号示意图

用u表示四边形单元上的电场或磁场强度，则电磁场值与子单元上的坐标可以表示为

由式(2)可以得出：

将式(6)代入式(5)可以得出

其中：

(12)

TE模式：

TM模式：

2 光滑度约束Tikhonov正则化反演

2.1 Tikhonov正则化反演

MT反演大多是高度不适定的病态问题，数据微小误差或较小扰动就会导致所求近似解严重偏离真解。为了提高反演结果的稳定性，使反演结果更加真实可靠，通常会在反演过程中引入Tikhonov正则化方法，其目标函数可以表示为

式中r为包含先验信息的参考模型；为加权矩阵；不同的对应不同类型的模型约束。本文采用差分矩阵来定义模型的粗糙度，并以此来约束模型的光滑程度。基于先验模型的光滑约束最小二乘正则化反演的迭代公式可以写为

2.2 正则化参数选取L-curve准则

通过求式(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 二维起伏地形模型图

正演基于不规则四边形的FEM方法，地下介质水平网格数为68个，垂直网格数为40个，单元网格的长×宽为100 m×200 m，频率选取10−1~103Hz范围内的21个频点。TM模式在横向上有更高的分辨率，在纵向上的分辨率较低。TE模式在纵向上有更高的分辨率，在横向的分辨率较低。因此，本文采用TM和TE模式联合反演，使反演结果优于单一模式的反演，具有更高的分辨率[21]。

图4 L-curve法求得的最优正则化参数

选用图5所示的正则化参数分别为0.05，0.25和4.00的反演结果。从图5可见：因不同导致反演结果存在显著差异，其中图5(c)所示的=0.25时反演剖面中异常体的位置和大小均清晰地得到反映，低阻异常体电阻率准确得到体现，但高阻异常体电阻率偏低。这是由于高阻异常体的正演分辨率较低，导致反演对高阻体不敏感。图5(a)所示为正则化参数=0.05时的反演剖面，由于正则化参数选取过小，背景平均视电阻率出现波动，且异常体的特征也有所形变。图5(b)所示为正则化参数=4.00时的反演剖面，当正则化参数选取过大时，电阻率异常与背景过于平缓，从而导致异常面积过大，模型约束过度导致反演结果不准确。

图5 二维起伏地形模型在不同的正则化参数条件下反演剖面

图6所示为不同正则化参数下的反演剖面与原模型剖面之间的误差曲线。按照反演结果的网格将原模型绘制出，计算出所有网格的均方误差，计算公式如下：

式中：n为模型网格总数；为原模型每个网格的电阻率；为反演网格的电阻率。以正则化参数递减的顺序计算每次结果与原模型的误差，可以从图6中看出L-curve法选取的最佳正则化参数的反演误差为曲线最低点，这说明通过L-curve法选取最优正则化参数能有效提高反演质量。

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.

(编辑陈灿华)

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)

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.

magnetotelluric(MT); irregular quadrilateral mesh; L-curve method; smoothness constrained inversion

P631

10.11817/j.issn.1672-7207.2018.03.016

1672−7207(2018)03−0626−07

2017−03−10；

2017−05−16

国家自然科学基金资助项目(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