二维直流电阻率与音频大地电磁自适应渐进联合反演

李 曼，于 鹏，张志勇，张宝松

（1.同济大学海洋与地球科学学院，上海200092；2.东华理工大学地球物理与测控技术学院，江西南昌330013；3.中国地质调查局南京地质调查中心，江苏南京210016）

联合反演可有效降低反演的多解性，提高反演结果的精度和可靠性［1-2］。直流电阻率（DC）与大地电磁（MT）联合反演研究始于20世纪70年代，一维研究表明，联合反演可避免单一方法反演的不确定性［3］；Constable等［4］采用平滑约束的直流测深与MT一维欧克姆（OCCAM）反演研究表明，MT、DC数据误差不一致会造成两个数据集拟合程度的不同，同时实际模型的非一维特征对反演存在影响。在二维方面，基于矩形剖分网格、平滑模型约束、互换定理计算偏导数矩阵的MT与DC联合反演实践表明，联合反演可以减少单独反演的不确定性［5-6］。类似MT方法，通过观测正交的电场与磁场分量，计算阻抗并定义视电阻率的射频电磁法即RMT方法（radiomagnetotelluric，频段10 kHz～1 MHz），由于其勘探深度与直流电阻率有很大的重叠，推动了联合反演工作。研究表明，联合反演RMT与DC数据可解决静态效应对RMT的影响［7］。单独反演与联合反演的对比研究表明，DC对高阻与低阻地质体均有较好响应，而RMT对低阻体更灵敏，通过合理选择两个数据集在反演中的权重可以得到更佳的反演结果［8］。实践工作表明，二维联合反演可以得到比三维张量反演更精细的地下结构［9］。采用相同模型约束条件下的概率反演与确定性反演的效果相当，联合反演改善了RMT反演效果［10］。

鉴于DC与MT、DC与RMT联合反演可有效改善反演效果，本文进一步开展非结构网格剖分条件下DC与音频大地电磁（AMT）二维联合反演的研究工作。开发了以模型灵敏度信息为依据的反演网格优化技术，以生成高质量反演网格；构建由粗网格到细网格的渐进反演策略，减少反演对稳定因子的依赖，从而简化正则化因子计算搜索的算法；通过最小二乘方法求解模型空间变化趋势方程系数，进而求解非结构三角网格的二阶梯度算子，构建非结构化网格条件下最小结构稳定因子；采用高斯-牛顿法优化求解正则化目标函数，利用双共轭梯度稳定算法（BICGSTAB）求解高斯-牛顿方程确保反演稳定，减少内存需求；最后，通过合成数据与实测数据的反演计算，验证了上述方法的有效性。

1 DC与AMT正则化联合反演

为更好地拟合地表与地下地质单元，实现野外工作点距差别大、勘探深度与分辨率不同的AMT与DC两种方法的高效、高精度正演，本文采用非结构化三角网格进行区域剖分；应用有限单元法进行DC与AMT正演［11-14］；采用基于最小结构模型的正则化方法构造非结构化三角网格条件下的模型二阶平滑计算方法，采用高斯-牛顿最优化方法求解反演目标函数，并通过双共轭梯度稳定算法（BICGSTAB）不完全求解高斯-牛顿方程，以确保反演过程稳定，同时减少反演对内存的需求。

1.1 正则化目标函数

地球物理反问题通常为病态问题，正则化是稳定求解反问题的有效手段［15-16］。正则化反演的目标函数为

式中：m为反演参数向量；dobs为反演数据向量；μ为正则化因子；φd(m，dobs)为数据误差函数；φm(m)为模型误差函数。

DC与AMT联合反演的变量为介质电阻率，采用以下变换函数［17-18］以确保反演电阻率在实际岩、矿石的取值范围内：

式中：mi为反演空间参数；ρi为反演单元的电阻率值；ρmin、ρmax分别为研究区电阻率的下限和上限。

反演数据由AMT与DC数据组成，dobs={ρTM，a，ρTE，a，φTM，a，φTE，a，ρDC，s}，其中 ，ρTM，a、φTM，a分别为音频大地电磁TM模式视电阻率与相位；ρTE，a、φTE，a分别为TE模式视电阻率与相位；ρDC，s为直流电阻率法观测的视电阻率。数据误差函数可表示为如下统一形式［16］：

式中：G(m)为正演函数，分别由AMT与DC两种方法观测数据组成；Wd为数据方差矩阵，

其中，ε为确保分母不为零的一个正实数；χi为数据方差；αi为数据权重系数，其作用是平衡AMT与DC数据在反演中的权重，以避免由于两种数据间误差范围的不同造成某一数据集的权重过大，而另一数据不起作用。在数据误差为高斯噪声的假设下，反演后数据误差函数期望值为NMT+NDC，NMT为AMT参与反演的数据个数、NDC为DC参与反演的数据个数。如果希望AMT与DC两个数据集在反演中取相同比重，数据权重系数可表示为

先验信息通常根据模型空间分布特征进行构造，根据其反演后模型的特点可大致分为平滑与陡变两种约束。其中平滑模型一般利用模型参数空间梯度进行计算。支持陡变模型的约束方法有多种，如最大变化量支持模型、最小支持模型、最小梯度支持模型等；另外，通过模型误差的L1与L2范数形式可以实现模型分块与平滑形式的反演约束。为确保反演的稳定性，本文采用最小结构模型约束。

式中：βs，βx，βy为模型误差三部分之间的比例系数；mapr为先验模型；模型误差可统一表示为

式中：Wm为统一定义形式的模型误差矩阵。

式（1）中μ为平衡数据误差与模型误差的系数（正则化因子），可采用广义交叉检验（GCV）、“L－Curve”方法、贝叶斯方法、经验选择方法等进行计算。正则化因子的选择关系反演的成败，若采用GCV、LCurve、贝叶斯方法等选取正则化因子，计算量较大，严重影响反演计算效率。本文研究了渐进反演技术，在反演过程中尽可能减少对稳定因子的依赖，因此采用简单的经验选择方法即可实现μ的选择。研究工作所采用的经验选择方法以数据误差、数据误差期望值等信息为依据，进行正则化因子的动态调整，基本策略如下：首次反演迭代设置正则化因子初值，设定数据误差的期望值；当迭代次数大于1时，如数据误差大于其期望值，按下式进行调整：

式中：γ＞1.0为调整系数；Du、Dd为数据误差下降速度的限制参数，当数据误差下降过快，增加正则化因子，若下降过慢则减少正则化因子，本次研究工作取Du=2.0，Dd=1.3。如数据误差小于期望值，则按比例增大正则化因子；另外，考虑到正演计算误差，在渐进反演初期采用较大的数据误差期望，随反演网格的细化减小数据误差期望。

1.2 非结构化三角网格模型梯度计算

矩形剖分时，模型粗糙度一般通过水平与垂直方向差分计算；而非结构三角网格，由于单元边的法向方向一般不在水平与垂直方向且各单元之间也不一致，所以计算更为复杂。根据单元的空间梯度与方向梯度存在线性关系，Lelièvre等［19］给出了非结构网格一阶梯度的计算方法；Key［18］采用了与计算单元存在共用顶点的三角单元进行粗糙度计算的方法，并在单元中心距离的计算中引入方向惩罚因子，实现水平与垂直方向不同强度的约束；Özyildirim等［20］采用了与计算单元具有共用边的三角单元，并通过边长与周长比值进行差分定义的粗糙度计算方法。由于Lelièvre算法意义明确，可将其扩展以实现二阶梯度计算。取第i个三角单元中心坐标为(xi，yi)，介质物性值为ωi，如图1所示。

图1 非结构三角网格梯度Fig.1 Gradient calculation of unstructured triangular mesh

假设物性分布为空间坐标的二次函数，有

通过与计算单元具有共用顶点的三角形组成的集合进行二阶梯度计算，用最小二乘求解式（9）的系数为

式中：q=(a，b，c，d，e，f)T；

进而求得二阶梯度为

式中：x0、y0为单位方向矢量。

由于采用了与计算三角单元具有共用顶点的所有单元进行计算，实际计算过程中式（10）一般为超定问题；如果式（10）欠定，可进一步向周边扩展计算集合，以满足最小二乘计算要求。

1.3 目标函数最优化求解

反演目标函数的求解为最优化问题，对于反演参数数量不大的一维问题可以直接采用奇异值分解；而对于二维或三维反问题常采用优化解法，如非线性共轭梯度法、高斯-牛顿法、拟牛顿法等。由于高斯-牛顿求解效率高，本次研究采用高斯-牛顿法求解。将式（3）、（7）代入式（1），可得第k次高斯-牛顿迭代公式为

式中：

其中，JMT为AMT数据对模型参数的偏导数矩阵；JDC为DC数据对模型参数的偏导数矩阵。偏导数矩阵是满阵，需要很大的存储量，而且其计算往往也是反问题求解中最耗时的部分。为有效地减少灵敏度矩阵计算的工作量，通常采用互换定理计算偏导数矩阵；基于伴随方程计算灵敏度矩阵的计算量与互换定理相当，采用均匀半空间解析解作为伴随方程解的近似方法可大大减少计算量。研究工作采用迭代法解式（12），应用了一种隐式的偏导数矩阵求解方案，将偏导数矩阵及其转置与向量的乘积转换为正演计算，避免显式存储偏导数矩阵与海森矩阵，进而有效地减少内存需求。式（12）中，数据变化量

为第k次正演得到的计算数据与观测数据的差。为确保每一次迭代目标函数可以充分下降，对模型更新步长进行线性搜索，有

χ为线性搜索步长，满足

则

求解式（12），由于右端矩阵条件数大，直接解法很难求得稳定解。采用正交分解方法进行求解可提高反演稳定性，但其计算量及内存需求较大；采用共轭梯度法（CG）、广义最小余量法（GMRES）等迭代解法进行不精确求解，稳定性较好；由CG发展而来的双共轭梯度稳定算法（BICGSTAB），计算效率高、稳定性好，本文采用该算法求解高斯-牛顿方程。

1.4 渐进反演算法

正则化技术依赖于先验信息，当先验信息不足或者不正确时将会产生错误的反演结果［21］，减少反演网格单元的数量与提高反演单元质量是实现稳定反演的重要途径。生成高质量网格的难度在于如何选择单元细化或组合的标准，最直接的方法是通过分析以下模型分辨率矩阵［15］来实现，然而对于二维和三维反问题，由于其计算量过大，难于实现。

研究工作利用模型灵敏度信息，进行自适应网格生成，定义模型灵敏度为

模型灵敏度Sj是模型mj改变时对整个观测数据集的影响。不考虑数据方差与稳定因子的条件下，Sj实质上反映了高斯-牛顿优化过程中海森矩阵的主对角元素。通过Sj对网格进行优化实质是调整海森矩阵的主对角元素，从而改善海森矩阵性质，进而得到高质量的反演网格。

基于研发的网格优化算法，进一步研究了一种变网格反演技术。变网格算法一般采用两种思路，一是先产生粗网格然后通过物性梯度［22］、多尺度算法［23］逐步细化；另一种则是先产生细网格，然后组合具有相似物性的单元［24］。此外，还有学者研究了基于角点与边检测技术的智能网格反演方法［25］。本文采用的渐进网格反演以模型灵敏度为依据进行网格细化，反演由粗网格到细网格逐步进行，反演流程如图2所示。计算过程如下：

Step 1读入数据，设置渐进网格反演迭代最大次数、网格优化比例、最小单元面积等信息。

Step 2从数据信息生成正演模型的几何描述，设定反演区域。

Step 3进行非结构化三角单元剖分。

Step 4当前网格开始最优化反演。设定当前网格高斯-牛顿反演迭代次数、BICGSTAB迭代次数、收敛条件、正则化因子初值，设置初始模型和参考模型（初始网格采用均匀半空间值，后面采用前一重网格反演结果作为初始与参考模型）。

Step 5高斯-牛顿反演迭代，求取模型改变量。

Step 6线性搜索，并更新模型，调整正则化因子。

Step 7当前网格进行高斯-牛顿反演终止条件判断，满足进入下一步，不满足返回Step 5。

Step 8渐进网格反演终止条件判断，满足退出反演，不满足进入Step 9。

Step 9互换定理计算雅克比矩阵，并计算模型灵敏度Sj。

Step 10分析模型灵敏度的平均值、方差，选择待优化单元。

Step 11返回Step 3，生成优化网格。

图2 渐进网格反演流程Fig.2 Flowchart of step by step inversion

渐进反演从粗网格开始，通过减少单元数量来减少反演对正则化因子的依赖；网格细化后，以前一重网格反演得到的模型作为下次网格的初始模型与参考模型，让反演逐步进行，进一步保障反演的稳定性，同时避免反演过程陷入局部极小或其他伪解。网格细化以模型灵敏度为依据进行局部细化，具体算法：对模型灵敏度进行统计分析，按设定比例选择模型灵敏度大的单元，本文研究取反演单元数量的10%进行优化；同时，为了避免过度细化网格，通过设定最小单元面积进行约束，即只对三角单元面积大于预设值的单元进行细化。

2 理论模型算例

设计如图3所示的二维模型，背景电阻率为100 Ω·m，共设置38个异常体，异常体均为截面为矩形的四棱柱，低阻体电阻率为10 Ω·m，高阻体电阻率为1 000 Ω·m。其中，接近地表交替有14个低阻和14个高阻异常体，异常体截面宽10 m、高30 m、上顶埋深为20 m，相邻两个异常体间隔为100 m；稍深交替有4个低阻与3个高阻体，异常体截面宽100 m、高100 m、上顶埋深为100 m，相邻两个异常体间隔为400 m；再向下，位于断面中部有1个低阻体，异常体截面宽400 m、高200 m、上顶埋深为300 m；最深部有1个低阻和1个高阻异常体，异常体截面宽500 m、高500 m，上顶埋深为500 m，两个异常体间隔为600 m。

近地表异常体用于模拟AMT的静态效应，静态效应严重干扰视电阻率，通常采用校正方法进行压制［26-27］；为避免校正方法引起的数据误差，可通过其他浅层方法进行补充［28-31］，进而开展综合解释和联合反演以提高解释精度。

图3 理论模型Fig.3 Synthetic model

直流电阻率法采用二极装置，共布设电极320根，相邻电极距10 m，最大极距200 m，数据集点数6 000个；AMT测点距40 m，共80个测点，测量频段1～4 096 Hz，按2的指数共13个频点，将正演数据加入3%随机噪声作为反演数据集。

图4为正演得到的TM模式卡尼亚电阻率与相位。二维条件下TM模式卡尼亚电阻率受静态效应影响严重，断面图出现条带状异常，无法识别深部异常体；静态效应对相位影响小，相较视电阻率断面图可以分辨更多的异常体。

采用本文介绍的渐进反演算法，分别对DC、AMT、DC与AMT数据进行渐进反演。反演过程对网格进行5次细化，反演结果如图5所示。

图4 理论模型TM模式正演结果Fig.4 TM-mode synthetic data

图5 理论模型反演断面Fig.5 Sections of synthetic data inversion

图5 a为直接反演DC数据所得断面，反演结果没有得到深部异常体；近地表的28个异常体反演效果较好，但受稍深的7个异常体影响存在畸变。图5b为AMT反演结果，从反演断面可以判断出几乎所有异常体，但是很多异常相互联结无法分离；近地表异常反演效果不理想，异常体位置有高、低阻体出现，但反演所得电阻率值与真实值相差较大，说明模拟设定的点距与频率条件下AMT对浅层不均匀体有一定的反演能力。图5c为联合反演两个数据集的反演结果，反演有效地给出了浅层与深部所有异常的位置；相比于单独反演，近地表28个异常体的电阻率与空间位置均较单独反演更接近真实值；稍深的7个异常体形态与设定模型有一定的出入，但异常体的位置与反演电阻率值均较单独反演精确；只有联合反演可分辨断面中部独立的低阻异常体；深部低阻异常的反演效果优于高阻，分析原因为在本文所采用的模型设置下，低阻二度体相较高阻二度体更易形成沿走向方向的电流场，有利于反演。

图6 反演区网格变化Fig.6 Mesh of inversion domain

图7 正则化因子、数据均方根误差、模型误差变化曲线Fig.7 Curves of optimum multiplier,RMS,and stabilizer

表1 渐进反演反演次数、RMS与模型误差统计表Tab.1 Iteration times,RMS,and model error of AMT inversion

为进一步分析渐进反演过程，现对AMT反演进行分析，反演过程网格变化如图6所示。表1给出了各次反演网格单元数量、总单元数量，反演迭代开始和结束时的数据均方根误差、模型误差值。图7为正则化因子、数据均方根误差（RMS）、模型误差变化曲线。反演区单元数从836到1 932，逐渐增加，第1、2次网格反演，RMS值快速下降，后面优化网格RMS值下降量逐渐减少，趋于稳定；反演过程中，初始网格正则化因子出现增大现象，随后稳定下降，表明渐进网格方法对正则化因子的初始依赖不大，可以通过自动调整确保反演稳定进行；最后两重网格数据误差下降不多，更多的是对模型进行优化；由于地表的灵敏度较大，网格总体由上到下逐渐加密，从图6c～6f看，由电阻率分布引起的网格变化过程清晰，体现了基于模型灵敏度的网格优化特点。

3 实测数据反演

采用本文的反演算法，对某沉积岩区的DC与AMT勘探数据开展联合反演研究，并与单独反演DC与AMT的结果进行对比分析。勘探剖面长度3 500 m，其中直流电阻率法采用10 m电极距、温纳α装置测量，最大供电极距900 m，共采集有效直流电阻率测点7 474个；AMT采用50 m点距，频率范围1～9 600 Hz，共41个频点，共采集有效测深点70个。

反演结果如图8所示，反演断面电阻率变化范围1～1 000 Ω·m，主体电阻率不高。DC反演断面，测线前600 m表现为中间高的三层结构，且近地表存在不连续低阻；而600 m以后，出现连续层状低阻体，其上部存在连续性较差的高阻，近地表出现不连续低阻；考虑到所采用的工作装置最大供电极距为900 m，勘探深度有限，300 m以下与初始模型电阻率相差不大的区域应当达到了直流反演的最大深度。AMT反演断面，在测线起点端左下方出现了一个中高阻体，并向浅部发展，在接近水平坐标500 m的位置接近地表；近地表500～680 m之间出现近直立低阻与高阻相伴异常；680 m后，电阻率性质与DC反演断面分布形态相似，但AMT反演断面，低阻的分布区明显比直流更大，这与直流的勘探深度不足有关。相较DC与AMT单独反演断面，联合反演得到了更为合理的底部高阻，高阻顶界面相比AMT变化缓，与实质地质条件相符。100～600 m深度出现了更精细的结构，分析原因为：①当极距增大时直流温纳α装置本身更适合反演层状结构。②工区浅部电阻率低，电磁法的趋肤深度小、垂向分辨率高。③AMT本身的横向分辨率较高，所以联合反演取得了理想的浅层识别能力；水平坐标500～680 m近直立低阻与高阻相伴异常形态更清晰；测线起点浅部高阻体边界更为明显。

图8 实测数据反演断面Fig.8 Sections of field data inversion

4 结论

为提高反演精度，本文开展了DC和AMT自适应渐进联合反演研究。通过对反演方法的研究和模型试算，取得以下几点认识：

（1）非结构化三角网格剖分适用于直流电阻率与音频大地电磁两种空间采样与分辨率差别较大的方法进行联合反演时对剖分的需求。

（2）基于模型灵敏度的渐进网格反演算法，减小了正则化方法对稳定因子的依赖，提高了反演的稳定性，同时减少了对正则化因子进行搜索的计算量。

（3）联合反演相较单一方法提高了反演模型的分辨率，可实现AMT产生静态效应异常体的直接反演。