E-Ex广域电磁法横向约束反演

张继锋，周光裕，刘最亮，孙乃泉，张富翔，王中圣

(1.长安大学 地质工程与测绘学院，陕西 西安 710054；2.华阳新材料科技集团有限公司，山西 阳泉 045000)

传统的可控源音频大地电磁法(CSAMT)一般要求数据采集范围在远区，沿用大地电磁的卡尼亚视电阻率解释方法，在实际工作中低频信号常常达不到远区，这将导致数据处理和解释与地下介质的电性变化不一致。广域电磁法是在CSAMT基础上提出的，它采用大功率人工源以提高信噪比，摒弃了仅在远区进行观测的缺点，拓展了观测范围，加大了勘探深度，且可采用单分量进行数据采集和解释，提高了观测速度和野外效率。因此，广域电磁法自提出以后，就成为许多科研工作者研究的热点，各种不同装置形式和不同场分量的测量方法先后被研究，包括基于水平电流源的-,-，-和-广域电磁方法。考虑到野外实际情况，采用水平电流源测量平行于源方向的电场分量-方式应用最为广泛。

广域视电阻率是一种被广泛采用的、能够直观反映地下电性特征的参数，是目前广域电磁法的主要处理解释手段。地球物理反演是进行精细化勘探的主要手段，其作用是从有限的观测数据中尽可能地恢复原有的地下物性参数，这也是地球物理工作者的最终目标。因此，在电法勘探中，希望得到更接近于地下真实电阻率的成像结果。目前，2.5D反演和3D反演可以得到更为接近地下真实模型的结果，但计算量较大，效率较低，因此一维反演仍被广泛采用，特别是在近似层状地层构造中。

早期CSAMT一维反演多采用先将近区和过渡区数据校正到远区，再以MT的方法进行反演，或仅采用远区数据进行反演。BARTEL和JACOBSON、石昆法采用了前一种方案；而SASAKI等利用阻尼最小二乘法实现了后一种方案。之后，ZONGE等亦利用最小二乘法完成了不做近场校正的一维CSAMT全资料反演；ROUTH在ZONGE的基础上考虑了反演参考模型和反演结果的光滑度；谭捍东基于Occam反演思想完成了水平电偶源电场的反演；王若和王妙月经过分析后认为全资料反演为今后发展方向，并结合网格参数法和剥层法完成了一维CSAMT反演工作，该方法不依赖初始模型，且对噪声有一定的抑制作用；尚通晓和朱威等详细比较了Occam反演和阻尼最小二乘反演在双极源全区反演中的应用效果，2种方法均可较好地拟合真实模型，但阻尼最小二乘反演对初始模型的要求更高；王堃鹏等将MT中自适应正则化反演方法引入到CSAMT全区反演中，一定程度上加快了反演速度；周军等将Levenberg-Marquardt反演方法引入到CSAMT全区数据反演中，改善了矩阵求逆的条件，使反演收敛更快速稳定；苌云等直接采用Cagniard视电阻率作为反演参数，得到了较理想的反演结果。

为了得到更接近地下介质真实电性的结果，一些学者开始研究反演目标函数中添加约束项的反演方法。KERRY Key给出了带约束的反演迭代公式，认为将地震解释中的构造约束整合到CSEM反演中具有重大意义，同时指出2个分量联合反演不会提高反演精度；张凯飞在传统的改进的共轭梯度法(NLGG)基础上添加了电性约束项，得到了较一般反演方法更优的结果，但同时指出，加入约束项后计算速度降低，且在紧约束的条件下易陷入局部极值而反演失败；索光运等采用解析法计算雅克比矩阵及井震约束的方式实现了-一维并行反演，有效提高反演速度；唐传章等实现了CSAMT的一维正演的并行化后，采用边界约束的有限内存拟牛顿算法完成了CSAMT一维反演，反演分层能力较好。

频率电磁法中由于地表不均匀体会导致严重的静态位移现象，这会导致视电阻率断面图出现所谓的“挂面条”现象，即使进行校正后，也不能完全消除静态位移。另一方面，采用一维反演，由于地下构造的变化和每个点受到干扰不一致，一维反演常常出现电阻率断面图的不连续现象，这与实际地层不符。2004年，AUKEN等提出剖面横向约束反演方法，其主要思想是利用横向约束稀疏矩阵将单点反演组合成一个集合，以保证每个单点迭代拟合观测数据的过程中，可以保证剖面上横向电阻率和厚度变化尽可能的小，从而达到反演后的剖面横向上呈现连续性；其后2009年，VALLEE等对横向约束反演中的横向约束稀疏矩阵进行了改进，提出利用二阶差分代替一阶差分形成横向约束矩阵，以此表示让反演测点与左右两点间的差异尽可能的小；2014年，蔡晶等对频率域航空电磁数据进行横向约束反演，并提出加权横向约束的方法，以保证剖面根据先验信息可以实现不同程度的横向连续性；2016年，殷长春等把加权横向约束的思想应用于航空瞬变电磁反演中，并用该方法对实测资料进行处理，得到的反演结果比单点反演结果效果更好，证明了该方法的有效性。综上所述，横向约束反演方法相比于单点反演方法具有一定的优势，在电磁勘探中已经得到了验证。广域电磁法是一种新的电磁勘探方法，目前许多学者从该方法的理论、数值模拟以及应用效果等方面进行了分析，解释手段还主要以广域视电阻率和单点一维反演为主，没有考虑测点之间的连续性，在煤矿富水区勘探方面应用较少。

笔者基于横向约束的思想，采用相邻测点之间的约束或通过已知测井先验信息进行约束，最大程度的减少-广域电磁电阻率断面的横向的不连续性，实现一维-广域电磁法横向约束反演。对广域电磁法在山西某矿区的实际测量数据进行了反演，划分了富水区和采空区范围，把该方法拓展到煤炭电磁勘探领域，并取得了较好效果。

1 一维横向约束反演

1.1 目标函数的建立

地球物理反演是寻找一个最优模型，使得该模型计算出的理论值与实际观测值之差的平方和最小。根据加权最小二乘法的基本原理，定义实测值与模型参数对应的理论值的拟合差为

(1)

式中，为数据拟合差，越大，则表示实测值和模型的正演理论值相差越大，亦即愈不符合反演期望值;为观测数据个数；为第个频点对应的实测数据；[]为对于模型参数第个频点对应的正演数据；为第个频点对应的实测数据的均方根误差。

(2)

为了方便计算，把式(2)可改写为矩阵形式

(3)

1.2 横向约束

可控源电磁法中相邻点之间的反演结果应当较为接近，各层电阻率和厚度在横向上具有连续性，特别当工区位于沉积岩地层发育的地区，横向连续性更加明显。但在传统单点一维反演方法中，得到的电阻率和层厚度会出现不连续的现象，即使是相邻两点间反演结果都会出现电阻率突变的现象，这使得剖面成像结果，出现很多小的突变异常，而大的构造和层位显示不明显，给解释工作增加了难度。

因此，在反演过程中需要增加约束条件，在目标函数中横向约束项和模型粗糙度，以保证相邻点之间的反演结果差异不会过大，纵向上模型具有光滑性。模型粗糙度可采用模型参数相对深度的一阶导数平方的积分表示，写成矩阵形式为

=(∂)∂

(4)

其中，为模型粗糙度；∂为一阶差分算子。横向约束项可采用

(5)

将式(5)改写为

=[(-)][(-)]

(6)

式中，为带先验条件的模型参数；为加权矩阵。

加权矩阵的具体形式：

(7)

同样的，加权矩阵亦为一个对角阵，若某一层位不存在先验条件，则加权矩阵对应的对角线元素应直接给定为0。

=

(8)

那么，先验条件的矩阵形式改写为

(9)

采用拉格朗日乘子将粗糙度矩阵和初始目标函数连接起来,此时目标函数最终表达式为

(10)

式(10)即为确定好的横向约束目标函数表达式，式中将用作权衡数据拟合与模型光滑程度的相对权重参数。采用扰动法计算电磁场值对模型参数的偏导数矩阵，该方法通用性好，适用范围广。

1.3 拉格朗日乘子的选取

拉格朗日乘子作为权衡数据拟合与模型光滑程度的相对权重参数，直接影响着反演工作的进行。当→0时，目标函数可看做仅含数据拟合项，粗糙度矩阵被忽略；反之当→∞时，数据拟合项趋于0，此时目标函数中仅含反演模型的光滑程度及先验条件项。

Occam反演原则：在满足数据拟合的基础上，使得反演模型尽可能的光滑。因此拉格朗日乘子的取值十分重要，一般先采用进退法寻找到1个极值区间，然后再将拟合差看做是拉格朗日乘子的函数，采用一维线性搜索算法(黄金分割法)使得达到最小。一般情况下，拉格朗日乘子取值区间一般为[10,10]，但在该区间内满足拟合差取极小的可能不止1个。出于使得反演模型尽可能光滑的考虑，因此在实际计算中我们选取满足条件的中的最大值。为了进一步缩小搜索区间，下面给出进退法的步骤流程：

确定一个较小的搜索区间后，可使用黄金分割法进一步加速搜索拉格朗日乘子，下面给出黄金分割法线性搜索的步骤流程：

(1) 假定此时用进退搜索法搜索到的对数区间为(,)，给定黄金分割数=0618，并令=-(-)，=+(-)；

该方法虽然计算量较大，但由于每次都采用最佳的拉格朗日乘子计算新的模型，收敛稳定且其总的迭代次数不多。

1.4 分段反演策略

按照逐点反演的形式完成反演剖面工作的，因此当该条剖面某点位附近存在测井等先验资料时，可以以该点位为分界点分段反演，并使用该资料作为目标函数中的先验条件项约束反演结果，其余点位反演则使用前一个相邻的测点反演结果作为先验条件进行约束；若不存在先验资料，则可以在实测资料中寻找电(磁)场值或广域视电阻率较稳定的点位，再以该点位为分界点分段反演，因为其受到噪声干扰较小，从而使得起始点位反演结果误差较小，尽可能地避免了整条剖面反演结果出现误差积累的现象。综上分析，当某条剖面存在多个先验资料或存在多个电(磁)场值或广域视电阻率较稳定的点位时，可以把该剖面分为几段，如图1所示。每段反演可认为是相互独立的，且计算流程相同。

图1 分段反演示意Fig.1 Piecewise inversion diagram

针对每一小段剖面，则应着重提高单点反演的效率。此外，广域电磁法一维正演中各发射频率对应的电磁场响应计算之间同样没有逻辑关联的，因此也具备很高的并行性。

2 模型试验

2.1 一维模型试验

本节将完成一维地电模型的反演计算，并将计算结果与真实模型作对比。发射参数：发射源长度=1 000 m，电流强度=30 A，电流方向与轴正方向相同，发射频率取[10,10]Hz上对数等间隔的35个频点，接收点坐标为(4 km,8 km)。反演参数：层数35层，采用平移算法计算出的广域视电阻率作为反演初始模型，并在真实模型的基础上加上最大5%的随机扰动作，矩阵中的对角线元素均取2。

首先设计H型和K型地电模型进行反演测试(采用单分量)。型模型参数：自地表向下厚度依次为=600 m，=200 m，电阻率依次为=100 Ω·m，=10 Ω·m，=200 Ω·m；K型模型参数：自地表向下厚度依次为=100 m，=500 m，电阻率依次为=10 Ω·m，=100 Ω·m，=40 Ω·m。

图2和图3(b)中黑线表示给定的理论模型，蓝线是在模型基础加了误差扰动项作为参考模型，绿线为不加先验约束条件的反演结果，红线为加先验约束条件的反演结果，2者均能较好地反映出中间层的电阻率，反演曲线首端和尾端能分别趋近于第1层和最后一层真实电阻率，H型地电模型中带约束反演结果(红线)在电阻率界面处较无约束反演结果(绿线)波动较小，中间层反演厚度范围更接近于真实模型厚度，K型地电模型约束反演结果在电阻率界面处较无约束反演结果波动稍大，但中间层反演厚度和电阻率更接近于真实模型。图2(c)和图3(c)表示拉格朗日乘子随迭代次数的变化情况，综合来看，随迭代次数的增加而逐渐减小，当→0时，目标函数可看做仅含数据拟合项，粗糙度矩阵及先验条件项被忽略。图2和图3(d)，(e)分别给出了H型和K型模型在有无先验条件下，采用单分量反演方式各分量的相对误差，各模型反演结果中各分量相对误差均小于5%，这表明原模型正演响应与反演结果拟合较好；此外通过图2和图3(d)可以看出，含先验条件的收敛速度较快，在H型模型中尤为明显。

2.2 二维模型试验

采用基于模型降阶的三维矢量有限元正演程序对二维地质模型进行正演，该方法适合于多频电磁法的正演计算，由于采用模型降阶算法，正演的速度不再受频点制约，频率越多，加速比越明显，速度可提高至少一个数量级，本文二维地电模型是在三维模型的基础上令走向方向无限延伸得到的。在正演结果的基础上加入5%的高斯白噪声，然后利用上述反演方法进行计算，以分析其在二维模型中的应用效果。

图2 H型模型反演结果Fig.2 Inversion results of H-type model

图3 K型模型反演结果Fig.3 Inversion results of K-type model

图4 地堑模型反演结果对比Fig.4 Comparison of inversion results of graben model

对比图4(b)～(d)可以看出，分量计算出的广域视电阻率结果对低阻体的显示最弱，卡尼亚视电阻率结果稍强，2者对地堑的反映很很弱；不加横向约束的反演结果和加横向约束的反演结果均能较好的反映地下电性分布结构，且与地堑埋深对应较好，但是从图4(c)可以看出，不加横向约束的反演结果电阻率横向连续性较差，与实际地电模型界面差别比较大，如图4中红色矩形框标出的位置，采用横向约束反演后，地电分界面的连续性改善了很多，基本光滑连续，与实际模型相符合。

3 实测资料处理

笔者对山西某地区实测资料进行处理试验，以进一步研究一维约束反演的实用性和有效性。

该工区受井田构造影响，总体为一走向近EW、倾向S的单斜构造，倾角为2°～6°，并在此基础上发育次一级的褶曲构造。在该区域大部分被新近系、第四系松散层所覆盖，仅在井田部分山坡及山梁零星出露有三叠系、二叠系地层。根据地表出露情况及钻孔资料，该工区地层从老至新依次为：奥陶系、石炭系、二叠系、三叠系、新近系和第四系。其中，石炭系上统太原组和二叠系下统山西组为主要可采煤层地层。太原组中8号煤层和15号下煤层及山西组中6号煤层稳定程度较差，属于局部可采类型，而太原组中9号煤层和15号煤层及山西组3号煤层稳定程度较好，属于大部可采类型。

该区域自下而上可划分为6个含水层，分别为奥陶系中统石灰岩岩溶裂隙含水层、石炭系上统太原组碎屑岩类夹石灰岩岩溶裂隙含水层、二叠系下统山西组砂岩裂隙含水层、二叠系上下石盒子组砂岩裂隙含水层、基岩风化带裂隙含水层和第四系砂砾石层孔隙含水层组；隔水层自下而上主要可划分3个，分别为中石炭统本溪组隔水层、石炭系和二叠系含水层间的泥岩隔水层及第四系和新近系松散层隔水层。

此外，该区域地貌形态属中低山地貌，由于长期的冲蚀作用，沟谷纵横、地形复杂，地势总体呈南高北低、西高东低的趋势，平均海拔在约1 150 m，最大高差约180 m，因此该勘探区浅、表层地电条件一般；但通过对该区测井及邻区瞬变电磁勘探资料分析，本勘探区中深部地层电性差异明显，煤系地层的平均电阻率值比非煤系地层电阻率高，煤系地层基底奥灰系灰岩电阻率最高，这具有较好的电性标识作用，因此本区域具备了电法勘探的应用前提，地层结构和电阻率统计见表1。

表1 岩层电性参数Table 1 Electrical parameters of rock formation

本次勘探主要是对岩层的富水性和煤层采空区进行探究，野外数据采集工作由JSGY-2型广域电磁仪完成，发射频点采用2，3，4和5四个频组共计28个频点，发射源长度为1 170 m，方向收发距约9 km，测线共计151个测点，点距20 m。

从图5可以看出，广域视电阻率断面图在纵向上基本能够定性的反映地层的电性变化，中间的低阻层是二叠系泥岩的反映，但其受地表地形起伏的影响大，浅部横向出现了不连续性，局部静态效应严重，异常区域也非常多。值得指出的是视电阻率深度计算采用经验公式，其反映的是地下电性层的整体变化，仅用于定性解释，其深度有一定的偏差。

根据实测采集装置设置反演相关参数，反演层数设置为28层，并利用相邻点的实测频率-电场曲线的相似程度实时调整横向约束大小，最后根据已知测井信息，将整个剖面分为几段完成反演工作。

图6中的测线，已知测点(距离为1 410～1 870 m)为一岩层富水区，对应图6中B区域，可以看出存在较强的低阻显示；根据最新资料，出水点靠近辅助进风巷，最新的工作面出水点靠近辅助进风巷60 m，位于勘探测线1 620 m，出水量为25 m, 位置刚好位于反演电阻率低值区。测点(距离为3 050～3 450 m)为一采空区，采空区右边界有大量充水，对应图中D区域，但深度稍向上偏移，C区域的低阻显示推测是由于地表山脊地形和采空区富水共同引起的，因此推测该采空区已经部分充水。

图5 实测数据广域视电阻率断面Fig.5 Wide field apparent resistivity of survey data

此外，整个一维横向约束反演剖面的横向连续性较好，地层分层较为明显，低频段视电阻率最高，推测对应于奥陶系的灰岩；中频段频率从低至高视电阻率由高值变为低值，推测视电阻率高值区域对应石炭系地层，视电阻率低值对应二叠系地层，存在不同厚度的泥岩，整体电阻率比较低，形成一个低阻层。由于煤层被大量开采，因此富水性较强，在横向表现出局部相对低阻异常，高频段整体呈现低-中电阻反映，对应于新近系和第四系地层。

经过钻孔验证，钻孔1和钻孔2施工过程中无出水情况，其对应深度处的电阻率相对较高，与钻孔情况基本吻合。钻孔3施工至42.5 m中粒砂岩涌水量0.5 m/h；施工至93 m中粒砂岩涌水量增大至1.5 m/h，压力0.6 MPa，该钻孔对应深度处的电阻率相对较低，接近于K2砂岩富水边界。

图6 实测数据反演剖面Fig.6 Inversion profile of survey data

4 结 论

(1)把基于横向约束的反演算法引入到广域电磁法中，并应用于煤炭富水区和采空区探测，在反演结果剖面上表现出较好的连续性，与实际地层基本吻合，取得了一定效果。

(2)以广域视电阻率作为反演的初始模型的一维反演结果表明：带参考点约束反演收敛速度较无约束反演快，反演结果更趋近与真实电阻率，特别是在中间层电阻率部分。

(3)通过二维地堑模型反演，发现加横向约束的电阻率反演断面光滑连续性，与实际地电界面吻合的更好，而广域视电阻率只能反映出地电模型第一层顶界面，对下界面的分辨率不够。

(4)对二维长剖面可采用分段反演策略，分段的起点选在有已知测井先验信息附近或受干扰影响很小的测点附近，以保证起始测点的反演结果可靠，在反演过程中可采用并行算法提高计算效率。

(5)横向约束反演应用于-广域电磁法中，但对于其他分量的广域电磁法反演同样适应。