基于全区视电阻率的瞬变电磁一维Occam反演中雅克比矩阵的解析算法

郭嵩巍,刘小畔,郑凯,张磊

(1.中国地质大学(北京) 地球科学与资源学院，北京 100083； 2.浙江省地矿勘探院，浙江 杭州 310000； 3.长江大学 地球物理与石油资源学院，湖北 武汉 430100； 4.内蒙古国土资源勘查开发有限责任公司,内蒙古 呼和浩特 010020)

0 引言

瞬变电磁法作为一种常见的电磁测深物探方法，广泛应用于水文勘查、工程勘查以及矿产勘查[1-2]。在目前的工程应用中，除“烟圈”反演快速成像法外，一维反演结果仍然是地质解释的主要依据。对于中心回线装置，可以解析推导出一维正演响应的计算公式[3]，其数值计算的核心是，内层积分的汉克尔变换[4]、外层积分的余弦变换[5]均可采用数值滤波算法，随着滤波系数的不断推陈出新，计算量明显减少，大大提高了正演的计算效率和精度。在反演方面，目前常见的反演方法，如阻尼最小二乘法(马奎特法)[6]、Occam法[7-10]、正则化法[11-13]等，都需要计算正演响应对模型的导数，即雅克比矩阵，因此雅克比矩阵的计算效率和效果就成为反演计算的关键因素之一。计算雅克比矩阵最简便的方法是差分法，该算法程序实现较为简单，缺点是不仅调用正演次数多，计算效率低，而且有可能增加累积误差。因此，本文详细推导了中心回线装置瞬变电磁一维正演全区视电阻率的雅克比矩阵的解析计算公式，并通过Occam反演，对比分析了理论模型和实测数据的反演效果，证明推导正确，值得采用。

1 瞬变电磁正演响应之全区视电阻率

目前，瞬变电磁仪器大多采集的数据为感应电动势，对其归一化后可转换为晚期视电阻率，然而晚期视电阻率是假设在时间趋于无穷大时近似推导所得，因此在响应早期，视电阻率曲线变形严重，呈明显高阻假象，成为了地质解释重要的干扰因素之一，因此很多学者提出了全区视电阻率(也称全期视电阻率)的概念[14]。由于通过磁场强度计算的全区视电阻率更为接近真正的视电阻率定义[15]，可通过感应电动势转换为磁场强度计算得出。因此，本文将该全区视电阻率作为Occam反演的拟合参数。

对于中心回线瞬变电磁的一维正演，可直接求解磁场强度。假设n层地电模型，第j层电阻率为ρj、层厚度为hj，频率域的磁场强度响应为：

(1)

其中：a为发射线框等效半径；u的递推公式[3]为

(2)

対于发射波形的选择，大多瞬变电磁仪普遍采用占空比1∶1的方波，如V8电法工作站的TD50方波，在正演模拟中只需考虑该波形的半个周期，可视为阶跃波：

(3)

式中I0为发射电流。对上式进行傅里叶变换，得到频率域中发射电流表达式：

(4)

于是，瞬变电磁响应可化简为：

(5)

(6)

对于均匀半空间地电模型，可以推导出Hz的解析表达式[3]：

(7)

其中erf(x)为误差函数。归一化Hz得到函数：

(8)

最后求取Z(x)的反函数x，得出基于磁场强度的全区视电阻率：

(9)

由于函数Z(x)为隐函数，无法推导出其反函数，因此无法通过解析算法求解，但可利用函数在0～1之间单调递增的特点，通过二分查找的办法得到数值解。笔者曾采用该方法实现了全区视电阻率的计算以及对应的“烟圈”反演，取得了不错的效果[15]。

可以看到，瞬变电磁正演响应计算可简化为：首先计算核函数，然后再对内层积分(计算频率域的电磁响应)做汉克尔变换，对外层积分(频率域到时间域的转换)做余弦变换(式(6))，得到Hz和Z(x)，最后采用二分查找法通过Z(x)求出x，从而得到正演响应的全区视电阻率。

2 雅克比矩阵的解析推导

式(9)对模型参数求导，根据复合函数的链式求导法则，层状介质模型任意层j的模型电阻率ρj及厚度hj的导数为：

(10)

(11)

利用反函数导数与原函数导数的倒数关系，得到：

(12)

根据Z(x)与Hz的关系式(8)，得到

(13)

将式(11)～式(13)带入式(10)，化简得到：

(14)

其中：

(15)

对核函数导数的计算，根据链式求导法则得：

(16)

(17)

令：

(18)

(19)

其中，

(20)

于是模型导数式(14)可表示为：

(21)

(22)

(23)

(24)

根据正演公式(式(1)、式(2))，进而推出第n层的模型导数关系：

(25)

(26)

3 Occam反演

Constable等提出Occam反演法[7]，其基本原理是寻找在目标函数相对极小情况下的光滑模型，它具有迭代过程稳定、收敛较快且对初始模型依赖度较小等优点。该方法在反演过程中每迭代一次需要对模型求一次导数，因此由模型导数构成的雅克比矩阵的计算精度和计算效率将决定Occam反演的计算效果和计算效率。Occam反演原理简述如下。

m为正演模型向量，d为正演响应向量，表示为：

dj=Fj[m],j=1,2,…,M

(27)

数据拟合差：

χ2=‖Wd-WF[m]‖2。

(28)

根据约束极小化理论，利用拉格朗日乘子形成一个无约束的函数U:

(29)

式中∂为粗糙度矩阵，可表示为：

(30)

式中右端第一项是粗糙度，第二项是由拉格朗日乘子加权的拟合差；取梯度，引起U不变的向量m遵循

μ-1(WJ)TWJm-μ-1(WJ)TWd+∂T∂m=0 ，

(31)

式中M×N的矩阵J，即是雅可比矩阵：

(32)

可表示为：

(33)

上式可通过式(21)、(22)、(26)计算得出。

最终，假设第k次反演迭代已完成，k+1次迭代的模型向量为：

mk+1(μ)=[μ∂T∂+(WJk)TWJk]-1(WJk)TWdk，

(34)

用一系列μ值计算模型mk+1(μ)的拟合差：

χk+1(μ)=‖Wd-WF[mk+1(μ)]‖，

(35)

其中：μ为拉格朗日乘子，χ为拟合差。

最终找到一个μ值相对最小，拟合差相对最小，光滑度相对高的模型，即是反演结果。

4 理论模型算例

为验证解析算法的正确性，本文采用与中间差分算法

(36)

做对比。从理论上说，差分步长越短，结果越接近真实值，但步长过短，会由于正演的精度不够而偏差增大，一般采用的经验差分步长(Δm/m)取0.001。

对典型三层地电模型做一维Occam反演：不失一般性，Occam反演的初始模型选择均匀半空间，电阻率取响应的平均值，模型层采用对数等间距递增，起始厚度5 m，对数递增间距为0.1，共分10层，最大反演深度687.2 m。

H型地电模型：层电阻率分别为100、10、100 Ω·m，层厚度分别为60、30 m，计算结果见图1、表1。K型地电模型：层电阻率分别为10、100、10 Ω·m，层厚度分别为30、60 m，计算结果见图2、表2。

从上述算例可以看出，差分算法不仅计算效率低，而且曲线拟合不佳，反演模型与真实模型差距大，而解析算法则表现优秀。差分算法效果不佳的原因很可能与全区视电阻率的计算有关，因为全区视电阻率在计算隐函数Z(x)的反函数时，是通过二分查找后插值所得，本身误差较大，再对其进行差分求导，放大了误差。因为若差分步长较大，差分本身的误差较大；而步长较小，正演响应Hz变化不大，二分查找的结果非常接近而趋于相等，差分结果趋近于零，总之会放大误差。随着反演不断迭代，误差会进一步累积放大，最终导致曲线拟合困难、反演模型震荡的结果(图1、图2)。因此，对于全区视电阻率的反演不建议使用差分算法。

表1 H型模型反演结果Table 1 Inversion result parameters of H-type layered model

图1 H型反演结果对比Fig.1 Comparison of H-type layered Model inversion results

图2 K型反演结果对比Fig.2 Comparison of K-type layered Model inversion results

表2 K型模型反演结果Table 2 Inversion result parameters of K-type layered model

5 实测数据实例

解析算法对理论模型的反演表现良好。为进一步验证解析算法Occam反演对实测数据的有效性，采用对实测剖面数据进行逐点一维反演的方法，即拟二维反演，绘制成电阻率反演断面图。与全区视电阻率拟断面、“烟圈”反演的结果进行效果对比，并通过地质钻探来验证方法的有效性。

5.1 地质概况

原始数据采集自内蒙古自治区凉城县西水塘镇附近，如图3所示，为内蒙古凉城县某地热勘查项目中的一条剖面(400线)，测线位于山前倾斜平原，布置该测线的目的是为划分工作区地层，推断基底埋深，并进一步验证地质推断的构造，为地热找水提供物探依据。工作区地形平坦，空旷且开阔，无工业电流干扰，非常有利于开展地面瞬变电磁工作。

根据区域地质资料显示[16]，工作区地层主要有太古宇界桑干群，中生界侏罗系及新生界古近系、新近系、第四系(图3)。区域地质构造属阴山东西向构造带的南缘，主要有NE和NW两个方向，其中NE向构造属高角度压性断层，断层连续性强，是地下热水和温泉形成的主要通道。侵入岩主要分布在蛮汉山区及平原区深部，主要为太古宙代早期侵入岩，岩性以苏长岩和似斑状花岗岩为主。

1—中低山; 2—山前倾斜平原; 3—第四系全新统冲洪积物; 4—侏罗系火山岩; 5—太古宇桑干群片麻岩; 6—太古宙似斑状花岗岩; 7—断裂破碎带; 8—测试瞬变电磁测深测线; 9—地质界线; 10—地貌界线; 11—已有钻孔；12—验证地热钻孔(DR2); 13—瞬变电磁发射线框位置; 14—瞬变电磁接收线圈位置1—middle low mountain; 2—piedmont inclined plain; 3—Quaternary Holocene alluvial proluvial; 4—Jurassic volcanic rock; 5—Archaean Sanggan group gneiss; 6—Archaean porphyry like granite; 7—fracture fracture Fracture Zone; 8—test TEM sounding line; 9—geological boundary; 10—geomorphic boundary; 11—existing drilling；12—verification of geothermal drilling (DR2); 13—TEM emission line frame position; 14—TEM receiving coil position图3 工作区地质图Fig.3 Geological map of work area

测线所在工作区全部为第四系覆盖。据区内已有钻孔显示(图3，图4a)，岩性从下至上分别为：太古宙似斑状花岗岩(γ)，多为灰色中粗粒似斑状花岗岩；古近系(E)为深棕红色泥岩、泥质砂岩、砂质泥岩夹砂砾岩；新近系(N)上部为浅棕红色泥岩、泥质砂岩、砂质泥岩夹砂砾岩，下部为土黄色、灰白色砂砾岩，与古近系地层呈假整合接触；第四系(Q)岩性为浅黄棕色、黄褐色粉土、粉质黏土与砂砾卵石、含卵砾中粗砂、中粗砂等。

图4 工作区钻孔柱状图Fig.4 Borehole histogram of work area

5.2 Occam反演及地质解释

工作装置的设置：根据工作区内已知钻孔显示(图4a)，工作区地层分层明显，侵入岩、沉积岩具有明显的纵向电性差异，具备了瞬变电磁测深施工的物性条件。测线所在位置基底埋深推断为在200～300 m之间，且地层中包括较厚的一套泥岩地层，方法试验也证实了这套低阻泥岩的存在。 为穿透这套低阻地层探测到基底花岗岩，在充分利用场地宽度的情况下，设置发射装置中心回线为400 m×600 m。

本着从已知到未知的推断原则，设定Occam反演初始模型。初始模型采用深度对数等间隔模型，让反演模型与瞬变电磁响应的采样间隔方式一致，这样一方面有利于曲线拟合，另一方面可以尽可能多分层，提高反演的纵向分辨能力。经过多次实验，以兼顾反演计算效率和突出局部异常为原则，确定初始模型起始深度为20 m，避开瞬变电磁响应早期的盲区，对数间隔1.15，分25层，最小层厚3 m，最大层厚65 m，最大反演深度498 m；初始模型电阻率参考全区视电阻率的均值，定为60 Ω·m。

瞬变电磁测深视电阻率曲线大多表现为KH、KQH型，反演结果总体可分为5个大电性层(图5)。第一层为中低阻，电阻率在 50 Ω·m左右，反演深度在0～30 m左右，与钻孔揭露的第四系粉土、黏土对应。第二层为中高阻层，电阻率在50～80 Ω·m之间，反演深度在30～80 m之间，对应以砂岩、砂砾岩、粉砂岩等为主的第四系地层，透水性好，为含水层。第三层为低阻层，电阻率缓慢下降，在10～50 Ω·m之间，反演深度80～150 m左右，与第四系地层粉质黏土对应，隔水性较好，为隔水层。第四层为低层，电阻率在10 Ω·m左右，反演深度150～300 m之间，与新近系，古近系的红色泥岩对应，为隔水层。第五层为高阻，电阻率大于100 Ω·m，对应基底太古代似斑状花岗岩。

图5 钻孔投影点820点反演结果Fig.5 Inversion result of site 820 where is projectionposition of borehole

上述地质推断，经钻孔DR2验证(图4b)，投影到剖面位置为820点，与实际情况基本吻合。该井定位于构造中，钻探深度 308.6 m，成井后自流，自流量为4 483.64 m3/d，水温 36 ℃，矿化度 1.01 g/L，是内蒙古凉城县地热勘查中取得的突出成果[17]。

5.3 反演结果对比分析

利用本实测数据地层分层明确、纵向电性差异明显的特点，可测试反演算法在纵向上的分辨能力，因此绘制了全区视电阻率拟断面和“烟圈”反演断面图，与Occam反演结果对比分析(图6)：拟断面和“烟圈”反演总体表现为“高—低—高”三层模型，中间层表现为巨厚的低阻层，对地层的分辨能力有限，而Occam反演则反演出第四系次高阻层，对应含砂岩、粉砂岩较多的地层，该类地层渗透性好，一般为含水层，对找水工作意义重大；对于高阻基底，上述两种的方法的反演结果基本一致，与钻探结果基本吻合。

a—拟断面;b—“烟圈”反演结果；c—Occam反演结果a—pseudo section; b—smoke ring inversion section; c—Occam inversion section图6 400线瞬变电磁反演电阻率断面对比Fig.6 Comparison of TEM inversion sections of line 400

6 结论

采用基于瞬变电磁全区视电阻率的Occam反演，有效避免了晚期视电阻率早期变形严重而造成的假高阻。在理论正演模型算例中显示：雅克比矩阵的解析算法在计算效果和计算效率两方面都优于差分算法；在实测数据算例中可以看到：Occam反演结果与钻孔更为吻合，对地层的纵向分辨能力高于烟圈反演，特别是在第四系低阻地层中次高阻砂岩层的反映，对瞬变电磁找水工作具有较高的实用价值。上述算例证明，解析算法公式推导正确，基于全区视电阻率的Occam反演不仅反演效率高，而且效果好。

此外，本文不仅推导了Occam反演中雅克比矩阵所需的对模型层电阻率的求导公式，还推导出了对模型层厚度的求导公式，可应用于对模型层厚度反演的方法，如马奎特反演等。

致谢：感谢内蒙古国土资源勘查开发院以及内蒙古自治区第一水文地质工程地质勘查院提供的宝贵实测数据及相关地质资料;感谢审稿老师提出的宝贵意见。