时间域航空电磁横向分段组合约束反演

林兴元 陆从德* 钟炳辉 范崇霄

(①成都理工大学地球物理学院，四川成都 610059； ②成都理工大学信息科学与技术学院，四川成都 610059)

0 引言

时间域航空电磁法(Airborne Time-domain Electromagnetic Method,ATEM)利用飞机上搭载的发射线圈向地下发送一次脉冲磁场(一次场)，地下的地质体由一次场激发出感应涡流，产生随时间变化的感应场(二次场)，同时使用飞机上搭载的接收线圈采集从地下返回的电磁信号，通过分析采集到的电磁信号，探测地下地质体[1]。ATEM具有良好的灵活性和机动性，勘查速度快、成本相对较低，受地形、地貌条件的影响较小等优势，能够应用于油气[2]和矿产资源勘探[3]、水文地质调查[4]以及环境监测[5-6]等领域。由于ATEM数据量巨大，且二维[7]和三维反演[8-10]需要大量的计算时间，因此很少应用于实际工程。一维反演因其计算量相对较小，在实际的ATEM数据处理和解释中有着广泛的应用[4-5]。

ATEM数据是全时连续采集，数据密度大，测点间距小，在大部分地质环境中，尤其是沉积环境下，相邻测点间的反演结果应具有较强的横向连续性。传统的一维反演方法，如Occam反演[5,11]、自适应正则化反演[12-13]、阻尼特征参数反演[14]等只考虑单个测点的情况，既使通过选取合理的初始值对反演进行优化[15]，其反演结果的横向连续性也较差。为此，Auken等[16]在反演直流电法数据时利用相邻测点间的横向粗糙度对模型进行约束(Laterally constrained inversion,LCI)，获得了较好的反演结果；此外，Auken等[17]还应用LCI对古河道模型的瞬变电磁数据进行反演，并与一维单点反演结果进行对比，结果表明LCI可以在一定程度上提高反演分辨率，同时还能抑制数据噪声的影响。Vallee等[18]将横向约束应用于ATEM数据，并详细地讨论了横向约束中一阶导数约束和二阶导数约束对反演结果的影响。为了突出层参数横向连续的重要性，蔡晶等[19]引入加权横向约束的概念。然而，由于LCI同时对多个测点进行计算,随着测点数的增多，模型参数的数量也随之增加，导致了计算时间急剧增加。Siemon等[20]对频率域航空电磁数据进行反演，并讨论了LCI模型参数的数量与计算时间之间的关系：当模型参数数量超过1000时，反演时间会急剧增加。Yu等[21]用前一个测点的反演结果作为先验信息约束当前测点，以此保证反演结果在横向上的连续性。但是这种横向先验信息约束很大程度上依赖于前一个测点的反演结果，因此反演结果具有明显的方向性。

ATEM勘探中，单条测线的测点数量从几千到几万不等，数据量巨大，因此对整条测线上所有测点同时进行计算，会耗费大量的时间。为了减少计算时间和提高反演精度，本文对整条测线上的测点进行分段计算，对每一段同时使用横向和垂向粗糙度约束以确保反演结果的连续性，段与段之间使用Yu等[21]提出的横向先验信息进行约束，确保反演结果在段间的连续性。

本文首先详细阐述这种横向分段组合约束反演的实现方案，然后使用模拟数据和实测数据对本文反演方案进行测试，最后分析上述三种约束对反演结果的影响，验证了本文反演方案的实用性和有效性。

1 反演方法

1.1 数据拟合

ATEM数据反演通过不断调整迭代模型或迭代模型参数m，使正演响应g(m)与观测数据dobs实现最佳拟合。将观测数据通过一阶泰勒级数展开，建立模型参数与观测数据之间的线性逼近关系[16,22]

dobs≅g(m)+G(mtrue-m)

(1)

式中：mtrue表示真实模型数据；G是雅可比矩阵，其第u行、第v列元素为

(2)

式中：du表示第u个正演响应数据；mv表示第v个模型参数。令

δm=mtrue-m

(3)

δdobs=dobs-g(m)

(4)

则式(1)变为

Gδm-δdobs=eobs

(5)

式中eobs为截断误差。本文采用阻尼最小二乘法对式(5)进行求解。由于只考虑了数据拟合误差函数，因此反演结果具有多解性。为此，除了拟合误差函数，还需对模型进行约束，以降低反演的多解性。

1.2 横向粗糙度约束

假设地质模型在横向上的电性特征变化较平缓，为了保证反演结果在横向上具有较好的连续性，需要施加横向粗糙度约束[16]

Rlsmtrue=els

(6)

式中：els表示模型参数在横向上的误差，其值越小，代表模型在横向方向上越光滑；Rls是模型的横向约束矩阵。将式(6)的左、右两边同时减去Rlsm，得到

Rlsmtrue-Rlsm=-Rlsm+els

(7)

则有

Rlsδm+Rlsm=els

(8)

横向约束矩阵Rls采用模型在横向上的一阶差分矩阵，即在横向相邻的两模型参数间使用一阶差分，其他位置均为0

Rls=

(9)

式中：S1=(M-1)N;T1=MN。其中M为测点个数，N为模型层数。

1.3 横向先验信息约束

横向先验信息约束是将前一个点的反演结果作为先验信息约束当前点的反演[21]。先验信息约束假设模型与先验信息的差距尽可能小，即

mtrue=mprior+elp

(10)

式中：mprior为模型的先验信息，本文指上一段最后一个测点的反演结果；elp是模型与先验信息的误差，elp越小，代表模型与先验信息越接近。将式(10)左、右两边同时减去m，有

mtrue-m=mprior-m+elp

(11)

得到

δm=mprior-m+elp

(12)

1.4 垂向粗糙度约束

垂向粗糙度约束[23]能使模型在垂直方向上具有较好的连续性。其推导过程与横向约束一样，因此与式(8)类似，即

Rvsδm+Rvsm=evs

(13)

式中下标“VS”代表垂直方向，各变量的含义与式(8)相同。与式(9)类似，Rvs的取值为

(14)

式中S2=MN,T2=M(N-1)。

1.5 横向分段组合约束反演

图1 横向分段组合约束反演示意图

根据图1中的横向分段组合约束反演的计算策略，将数据拟合方程(式(5))、横向粗糙度约束方程(式(8))、横向先验信息约束方程(式(12))以及垂向粗糙度约束方程(式(13))联合，可以得到总的反演方程

(15)

令

(16a)

(16b)

则式(15)可以写为

G′δm-δd′=e′

(17)

根据文献[24],式(17)的解为

δmest=[(G′)TWG′]-1(G′)TWδd′

(18)

式中W为模型参数和数据项的权值矩阵，可以写成

(19)

其中Wobs、Wls、Wvs和Wlp均为对角加权矩阵，Wobs为数据项的加权矩阵，Wls为模型项的横向约束加权矩阵，Wvs为模型项的垂向约束加权矩阵，Wlp为先验模型的加权矩阵。为了使式(18)解更稳定，引入阻尼因子，将其改为迭代形式

(20)

式中：λn为第n次迭代的阻尼因子，求解方法见文献[25]；mn为第n次迭代的模型参数。将式(16)和式(19)代入式(20)，有

(21)

2 实验验证

2.1 理论数据反演

为了验证该方法的有效性，本文使用2.5维正演程序[26]模拟ATEM响应数据。ATEM数据采集系统见图2，发射和接收系统主要参数为：发射线圈面积为123.9m2，匝数为5；发射基频为25Hz，发射波形为三角波，上升沿时长为1.5ms，下降沿时长为1.50ms，峰值电流为300A；接收和发射线圈离地高度均为30m。理论模型包括3层：第1层的厚度为40m，电阻率为100Ω·m；第2层电阻率为5Ω·m，两边厚度为40m，中间厚度为100m，在中间位置有一宽度大约为200m、厚度逐渐变化的区域；第3层是基底，电阻率为100Ω·m。测线长度为1600m，测点间距为25m，测点个数为65。在测量发射电流关断后0.04～4.00ms时间段内按照对数等间距抽取21道的感应电动势(图3)。为了尽量接近实测数据，对模拟的电磁响应数据添加5%的高斯白噪声。

图2 航空电磁采集系统示意图

图3 理论模型的模拟电磁响应曲线 图中曲线从上到下分别对应衰减早期到晚期的感应电动势

反演初始模型是电阻率为100Ω·m的均匀半空间，共26层，每层厚度均为10m。图4b～图4g为不同方法的反演结果，其中图4b～图4d和图4g采取分段反演方法，即测线上的所有测点分为7段，前6段的测点数均为10，最后一段的测点个数为5，而图4e采用单测点依次反演，图4f则是所有测点同时反演。表1为不同反演方法的计算时间及反演结果与理论模型的均方根误差。本文所用的计算平台CPU型号为E3-1505M(8核)，内存为16GB。

从图4b～图4d、图4g可以看出：无约束的分段反演结果(图4b)整体效果较差，边界模糊，中间部位电阻率值波动较大。进行了垂向粗糙度约束(图4c)的反演结果得到了一定的改善，中间部位基本连成一块，但是横向连续性依然不理想。图4d在图4c的基础上增加了横向粗糙度约束，可见反演结果在横向上的连续性得到了改善，但是由于段与段之间没有进行横向约束，反演结果出现了分块现象。图4g是在图4d的基础上在段与段之间增加了横向先验信息约束，其结果与图4a中的理论模型比较接近。对比图4d与图4g可以发现，横向先验信息约束可以消除由分段引起的反演结果的分块现象，提高反演的精度。

图4 模型数据反演结果

(a)理论模型; (b)无约束分段反演; (c)垂向粗糙度约束分段反演; (d)垂向和横向粗糙度约束分段反演; (e)垂向粗糙度约束和横向先验信息约束单点反演; (f)垂向和横向粗糙度约束反演; (g)横向分段组合约束反演

图中黑线为理论模型界面

结合表1可以发现，在ATEM数据反演中，增加适当的约束可在提高反演结果准确性的同时减少计算时间。

对比图4e与图4g可以看出： 图4e整体连续性较好，但由于采用前一个测点的反演结果作为当前测点的先验信息来保持模型在横向上的连续性，这只是单向约束，导致反演结果具有明显的方向性，因而与图4a中理论模型相差较大； 图4g的横向分段组合约束反演是在段与段之间采用单向的横向先验信息约束，而在段内采用双向约束的横向粗糙度约束，因此，反演结果受单向横向先验信息约束的影响较小，没有出现类似于图4e中明显的方向性。因此，图4e虽计算时间较少(表1)，但是反演结果的异常体的边界不准确，误差较大； 图4g的计算时间稍长，但是反演精度较高，能真实地反映地下介质的电性分布特征。

对比图4f与图4g可以看出，二者与理论模型(图4a)较为接近，边界轮廓比较清晰，整体的连续性也较好。图4f的计算时间是图4g的近80倍，且误差也稍大。因此，横向分段组合约束反演可以在确保反演精度的同时大幅度减少计算时间。

表1 理论模型不同反演方法计算时间及均方根误差(RMSE)

2.2 实测数据反演

为了验证该横向分段组合约束反演方案的实用性，对美国地质调查局(USGS)网站上公开的野外ATEM勘探数据(图5)[27]进行反演。

每10个测点分为一段，反演结果见图6。对比图6a～图6d可以发现，无约束的反演结果(图6a)整体连续性较差，仅能够反应地质体的大致位置，边界较模糊；增加了垂向约束的反演结果(图6b)在垂直方向上的连续性增强，但是地质体的形状和轮廓仍然较为模糊，且第三层的反演效果不佳；增加了横向粗糙度约束的反演结果(图6c)的整体连续性增加，但是分段导致相邻段之间的反演结果差别较大，不能很好地突出地质体的边界和轮廓；图6d整体较光滑，但是地质体的位置和形状不能很好地得到反映；横向分段组合约束的反演结果(图6e)所反映的地质体位置、形状和边界较清晰。

图5 实测电磁响应曲线

图6 实测数据不同方式的反演结果

(a)无约束的分段反演结果; (b)垂向粗糙度约束的分段反演结果; (c)垂向和横向粗糙度约束的分段反演结果; (d)垂向粗糙度约束和横向先验信息约束的单点反演结果; (e)横向分段组合约束反演的反演结果

对比文献[26]提供的地质资料，对图6分析如下：

(1)图中①处存在着一条山脉，分布着大量的基性火成岩，其电阻率大于1000Ω·m，与反演结果基本一致；

(2)图6e中①、②之间的A处存在断层,而图6e中的电阻率在①、②之间的断层A处的左右对比明显，能够很好地展现出断层的垂向变化；

(3)图中的②、③、④位于盆地中，其中②处的表层是由松散砂砾石组成的冲积扇沉积物，电阻率为80～150Ω·m,反演结果能够很好地体现，③位于露头花岗岩附近，图中这个位置也呈现较高电阻率，④处富含粘土，具有较低的电阻率；

(4)图中的⑤处位于盆地边缘，主要由长英质火山岩组成，其电性与花岗岩相近，即具有较高的电阻率，与反演剖面结果一致；

(5)图中B处存在一个横切的断裂带，这可能是导致⑤处内部电阻率等值线不连续的原因；

(6)图中⑥处表层黄色位置为长英质火山岩，具有较高的电阻率。蓝色位置由较新的中等冲积扇沉积物组成，主要包含分选不良、致密的胶结砂和砾石，具有较低的电阻率(<10Ω·m)。

由此可见，实测ATEM数据反演所得到的构造、断层及地质体的位置、形状和边界较清晰，反演结果与实际的地质资料吻合较好，说明本文提出的反演策略可有效地应用于实测ATEM数据的反演。

3 结论

本文将垂向粗糙度约束、横向粗糙度约束以及横向先验信息约束三者相结合，对ATEM数据进行分段反演。模拟数据和实际数据的反演结果说明：

(1)适当的约束条件可提高ATEM数据反演结果的准确性，同时也能减少计算时间；

(2)横向先验信息约束可消除由分段引起的段间不连续现象；

(3)横向粗糙度约束能有效地降低由横向先验信息约束引起的方向性。