近几年航空电磁发展迅速,不仅可应用于矿产和油气资源勘查,也可用于环境工程、地下水、地热资源和海洋地形调查等(殷长春等,2015).特别是近年发展起来的时间域航空电磁法,由于其同时具备发射功率、探测深度大和横向分辨率高等优势,逐渐成为航空电磁勘查的主流技术.
时间域航空电磁数据处理和反演已经从早期的近似成像技术发展到基于层状介质模型的一维反演及逐步走向实用化的二、三维反演.然而,航空电磁数据量巨大,同时考虑到横向采样率高导致相邻测点地质条件变化较小的特点,一维反演仍是航空电磁数据解释的主要技术手段.常用的一维反演技术大体分为最速下降搜索和全局搜索方法.基于最速下降的阻尼最小二乘法要求最大程度的拟合数据,同时沿负梯度方向搜索收敛速度快,但对初始模型的依赖性较强、容易陷入局部极小,对剖面数据反演会出现纵向和横向上电阻率不连续的现象.Constable等(1987)提出Occam反演,通过在目标函数中引入粗糙度来约束模型单元之间垂向电性变化,从而得到纵向上光滑的反演结果.Occam反演稳定可靠、不受初始模型限制,但该方法耗时较大,而且由于方法本身的局限不能清晰刻画地下地层界面.LCI(Laterally-Constrained Inversion)反演可以通过对相邻测点的电阻率和厚度施加横向约束,以保证相邻测点电性连续性,同时可以加入剖面上任何测点处的先验信息并有助于控制整个测线的反演结果(Auken and Christiansen, 2004;Auken et al., 2008).蔡晶等(2014)提出加权横向约束反演方法,并针对频率域航空电磁数据与阻尼最小二乘反演结果的对比,指出加权横向约束反演地电断面的层界面更加光滑.以上这些梯度类反演方法的有效性严重依赖于初始模型,不合理的初始模型容易导致反演收敛于局部极小.模拟金属淬火反演(Simulated Annealing)和贝叶斯(Bayesian Inversion)反演是求解全局最小的反演方法.模拟金属淬火是通过模拟熔融金属冷却结晶的热力学过程,按照玻尔兹曼分布进行上行和下行搜索全局最小值的方法(Yin and Hodges, 2007).变维数贝叶斯反演方法由殷长春等(2014)引入到频率域航空电磁数据反演.该方法根据建议分布,利用蒙特卡洛方法充分搜索模型空间,只接受同时满足接受概率和数据拟合要求的候选模型,最终获得反演模型的最大概率分布和不确定信息,但该方法在随机搜索过程中需要多次计算雅克比矩阵会产生巨大的计算量.
在地下水和地热资源勘查、海侵、垃圾填埋场渗滤液及废弃矿山等地质调查时,对于层界面的刻画显得至关重要,但一般要求地下介质分层明显且电阻率差异较大.此时,Occam反演或其他基于L2范数的光滑反演难以对地下陡变的层界面进行有效刻画.为了提高界面的分辨率,学者们提出了L1范数反演和聚焦反演技术.Farquharson和Oldenburg(1998)提出在电磁数据反演中最小化模型垂直梯度的L1范数;Farquharson和Oldenburg(2003)、Loke等(2003)在二维直流电阻率反演中通过对模型的水平和垂直梯度的L1范数进行最小化,反演结果清晰刻画出目标体的水平和垂直界面.Portniaguine和Zhdanov(1999)提出的聚焦反演是基于最小梯度支撑稳定泛函,实现模型参数变化剧烈和不连续区域的最小化,得到清晰和聚焦的地质结构,并成功应用于重、磁、电磁三维反演中(Portniaguine and Zhdanov, 2002;Zhdanov,2009).除了传统的Lp范数反演,徐义贤和王家映(1998)利用小波变换的多尺度特性来处理大地电磁数据,改善了传统反演方法容易陷入局部极小的弊端.Chiao和Kuo(2001)提出了基于多分辨率小波表示的模型参数化方法,在地震数据处理中有效地提高了反演的分辨率,然而目前尚未在电磁反演中得到应用.
本文针对航空电磁用于环境工程和地下水勘查中对垂向分层的特殊需求,并考虑到传统阻尼最小二乘反演无法实现准确刻画地下电性界面的技术难题,而全局搜索反演方法因计算量巨大难以实用化的困难,我们提出基于广义模型约束的电磁反演方法.该方法通过改变模型约束项或模型的求解域构造反演方法,获得陡变层界面的反演结果,主要包含Lp范数反演、聚焦反演和基于小波变换的稀疏约束反演.我们以时间域航空电磁一维反演为例,通过设定不同的理论模型进行测试,以展示它们对于陡变层界面刻画的优势和缺点.进而,通过对时间域航空电磁实测数据反演,进一步验证方法的有效性.
1 广义模型约束反演原理Tikhonov1963年提出的正则化理论已被广泛应用于解决电磁反演问题.本文从正则化理论出发提出基于广义模型约束的反演方法,其目标函数定义为
(1) |
其中φd(u)为数据拟合项,φm(k)为模型约束项,β为正则化因子,u和k的具体形式为
(2) |
(3) |
式中,dobs为观测数据向量,dprd是由模型m计算的电磁响应,mref为参考模型,Wd为一个对角阵,其元素是观测响应中噪声的倒数.
φd和φm的一般形式可表示为
(4) |
其中,x表示(2)和(3)式中的u或k.第n次迭代的目标函数可表示为
(5) |
而有:
(6) |
(7) |
式中,dn-1是由上一次迭代得到的模型mn-1计算得到的响应向量,δm=mn-mn-1,J是灵敏度矩阵.对目标函数求极小等价于使其导数为零,此时得到的模型参数向量m即为反演问题的最优解.
将(4)式对模型参数求导(Farquharson and Oldenburg, 1998)可得:
(8) |
式(8)也可表示为
(9) |
式中, ∂φ/∂δm=(∂φ/∂δm1, …, ∂φ/∂δmN)T,N为介质层数,Bij=∂xi/∂δmj,q=(ψ′(x1), …, ψ′(xN))T.引入对角矩阵R将公式(9)重写为
(10) |
其中,广义约束矩阵定义为R=diag{ψ′(x1)/x1, …, ψ′(xN)/xN}.对于目标函数中的数据拟合项,(10)式中B=-WdJ;而对于模型粗糙度项,B=Wm.因此,对(5)式求导并令其为零,可得每次反演迭代中求解的线性方程组为
(11) |
其中,Rd为元素为2的对角阵(对应于数据拟合项通常采用L2范数),Rm将在后续讨论中给出细节.广义模型约束反演中约束项形式灵活多变,即公式(4)中的ψ(x)可选择不同形式.通过改变模型约束项可以构造出不同反演方法,比如Lp范数反演和聚焦反演;通过改变模型的求解域可以构造出基于小波变换的稀疏约束反演.下面我们讨论这些反演方法的原理及各自特点,并通过算例验证其有效性.
1.1 Lp范数反演不同于常规的Lp范数,Ekblom(1987)提出的Lp范数形式为
(12) |
式中ε是一个小值.上述的Lp范数相对于常规Lp范数形式更具有适用性.其原因是当p=1时,在x=0处它的导数存在.相应地,公式(11)中Rm矩阵的对角元素可以表示为
(13) |
当p=2时,Rii=2,可构造出L2范数反演;当p=1时,Rii=(xi2+ε2)-1/2,可构造出L1范数反演;而当p=0.8时,
众所周知,L2范数反演具有光滑连续的特点,而L1范数反演则可以实现陡变界面的有效反演.这是由于这两种范数的差别造成的.事实上,L2范数是指向量中各元素的平方和求平方根,L2范数中每个元素自乘(求平方)意味着值大的元素对整体贡献较大,因此最小化一个模型的L2范数会产生相对缓慢变化的模型.Occam反演本质上是L2范数反演,它的模型粗糙度项表示的是相邻单元电阻率的差,最小化模型粗糙度导致相邻单元电阻率差值最小化,因此反演结果光滑.L1范数也被称为“稀疏规则算子”,是指向量中各个元素的绝对值求和,选择L1范数是因为L1范数可以实现稀疏规则化.稀疏化算子的引入可以在反演中去掉没有信息的特征,即把这些特征对应的权重设为零.有信息的特征对应地下介质的异常区域,而没有信息的特征对应地下介质的背景区域.因此,稀疏化突出了地下异常区域的分布特征.为了便于说明,假设Lp范数的常规形式为‖w‖p,矢量w包含两个数w1和w2.在(w1, w2)平面上可以画出数据拟合和模型约束的等值线,参见图 1.对于L2范数,满足关系w12+w22=C(C为一常数),对应的模型约束等值线是一个圆形;对于L1范数,满足关系|w1|+|w2|=C,对应的模型约束等值线是一个正方形;我们同样可画出L0.8范数对应的模型约束等值线图.目标函数是由数据拟合和模型约束两部分构成的,因此目标函数的最小解对应的是这两项的共同解,即这两项等值线相切的位置.由图可以看出,L1和L2范数的不同之处在于L1范数在和每个坐标轴相交的地方都出现“角”,在角的位置上容易产生稀疏性;而对于L2范数,目标函数等值线难以与L2范数等值线相交在角的位置,因此不容易产生稀疏性;而L1范数与L0.8范数相比,目标函数等值线更容易与L0.8范数在角的位置相交,这就说明L0.8比L1范数更容易实现稀疏化.因此本文提出L0.8范数反演可以得到更加陡变的一维反演结果.
Last和Kubik(1983)及Portniaguine和Zhdanov(1999)采用如下形式的泛函稳定器ψ(x)实现聚焦反演,即:
(14) |
式中,当ε很小时,该约束项基本上与矢量x中非零元素的数量成正比,其对应的Rm矩阵的元素为
(15) |
图 2说明了上述几种不同形式ψ(x)的变化曲线,其中a、b、c曲线分别代表L2范数,L1范数和聚焦反演泛函稳定器,d曲线代表Ekblom提出的Lp范数形式.对模型约束项进行最小化实质上等价于ψ(x)的最小化(趋于零).由于在零附近聚焦反演泛函稳定器曲线的斜率比L2范数和L1范数大,因此其趋近于零的速率更大.这就意味着最小化模型约束项时,会导致模型参数梯度快速趋于零,由此模型变得连续;当x值很大时,ψ(x)的值趋近于1,此时模型的梯度最大,可以很好反演模型参数的陡变特征.对于聚焦反演,模型粗糙度项只存在趋于1和0两种情况,分别代表异常区域和背景区域.因此,与L1范数反演相比,聚焦反演结果在背景区域更连续,而在异常区域会产生陡变特征.
小波变换可以将空间域模型转换成两组小波域系数.其中,一组为近似系数,而另一组为细节系数,它们分别包含了模型的全局信息和局部细节信息.小波变换可以表示为
(16) |
式中,
(17) |
式中,Ww-1是逆小波变换算子.由于小波系数通常是稀疏的,因此本文使用小波系数的稀疏约束直接在小波域进行反演求解.根据Donoho(2006),L1范数可以产生稀疏解,本文使用小波系数的L1范数实现反演,迭代过程中求解的线性方程组变为
(18) |
由于本文对数据拟合项使用L2范数,对模型约束项使用L1范数,因此上式中Rd仍然是元素为2的对角阵.对于Rm矩阵,本文使用Ekblom(1987)提出的具体形式,由公式(12)给出.由于小波变换本身具有一定的光滑约束,因此这里不再需要光滑算子Wm.
小波类型的选择对于模型特征提取至关重要.Haar小波由Alfréd Haar于1909年提出,基于Haar小波的小波变换也是最简单的一种变换.Haar小波的基函数是一组由分段常数组成的函数集合,其形状是一个定义在区间[0, 1)上的矩形.Haar小波基函数可表示为
(19) |
而图 3给出Haar小波基函数的波形图.
对于地下水平层状介质,如何很好地反演出层厚度和层界面位置至关重要.Haar小波基函数是阶跃波,因此基于Haar小波变换的小波域反演更容易得到陡变的层界面的反演结果.Haar小波变换的原理如下:假设模型参数为一维数组,首先利用一次Haar小波变换将模型参数转换成小波域系数,即依次从数组中取两个数,计算它们的均值(两数之和除以2)以及差值(两数之差除以2),并将均值和差值依次保存在数组的前半部分和后半部分.前半部分的均值称为近似系数,其包含全局平滑信息;后半部分的差值称为细节系数,其包含局部信息.通过定义一个对分界线来区分这两组系数.进而,我们再对前半部分近似系数进行一次Haar小波变换,进一步提取近似系数中的近似系数和细节系数,如此不断地进行变换,直到不可再分为止.最后再对所有系数进行反演,并将反演后的小波域系数恢复为空间域的模型参数.利用小波变换的这种多尺度特性,将空间域的模型转换到小波域进行反演,这样不仅可以得到高分辨率的反演结果,同时也会提高反演速度.
2 理论模型反演算例本文研究的时间域航空电磁系统采用CGG公司的HELITEM MULTIPULSE系统,发射波形是半正弦波加上一个梯形波,半正弦波形宽度4 ms,峰值电流435 A,梯形波在半正弦波断电后10.5 ms开始发射,上升沿宽度0.05 ms,平稳供电0.9 ms,下降沿0.05 ms,发射线圈半径17.5 m,接收机在发射线圈上方0.5 m.为了验证Lp范数反演、聚焦反演和基于小波变换的稀疏约束反演在时间域航空电磁数据处理中对于层界面的刻画能力,我们利用添加3%的高斯随机噪声的理论模拟数据进行反演,并与Occam反演结果进行对比.我们设计如下初始模型:电阻率值为50 Ωm的均匀半空间,位于地表的初始层厚为2 m,往下各层厚度依次递增.
2.1 三层模型反演本文首先讨论三层H型模型的反演结果.第一层的电阻率为200 Ωm,厚度为50 m,第二层电阻率为10 Ωm,厚度为30 m,第三层基底层电阻率为100 Ωm.由图 4给出的反演结果可以看出:与Occam反演相比,本文提出的反演算法均能较好地刻画层界面,其中尤以聚焦反演和基于小波变换的稀疏约束反演结果与设计的模型最接近,对于层界面的刻画非常准确.与L1范数相比,L0.8范数反演在层界面处得到更为陡变的反演结果.计算速度方面:Occam反演用时80.5 s;L1范数反演用时21.4 s;L0.8范数反演用时26.5 s;聚焦反演用时36.3 s;基于小波变换的稀疏约束反演用时25.2 s.
为进一步验证本文提出的反演算法的有效性,我们设计一个四层模型.第一层的电阻率为50 Ωm,厚度为30 m,第二层电阻率为300 Ωm,厚度为70 m,第三层电阻率为10 Ωm,厚度为30 m,第四层基底电阻率为100 Ωm.从图 5给出的反演结果可以得出与图 4相似的结论,本文提出的所有反演算法对地下电性界面均能很好地刻画,界面位置与真实模型吻合很好.计算速度方面:Occam反演用时94.1 s;L1范数反演用时17.5 s;L0.8范数反演用时24.6 s;聚焦反演用时23.4 s;基于小波变换的稀疏约束反演用时22.2 s.
我们最后设计一个复杂的五层模型.第一层电阻率为200 Ωm,厚度为50 m,第二层是电阻率为10 Ωm的低阻层,厚度为30 m,第三层电阻率为100 Ωm,厚度为60 m,第四层是电阻率为5 Ωm的低阻层,厚度为50 m,基底电阻率为20 Ωm.图 6给出的反演结果再次证实本文提出的算法在反演地下介质陡变层界面时的有效性.计算速度方面:Occam反演用时75.3 s;L1范数反演用时19.6 s;L0.8范数反演用时25.8 s;聚焦反演用时35.2 s;基于小波变换的稀疏约束反演用时17.2 s.
综合上面各理论模型反演结果,我们可以得出如下结论:Occam反演稳定,但得到的是光滑反演结果,不能准确刻画出层界面的厚度和位置.L1范数反演与Occam反演相比可以很好地刻画层界面位置,但是L1范数反演不能刻画陡变的拐角,同时L1范数反演不容易收敛.相比之下,L0.8范数反演可以得到陡变的拐角,同时反演结果界面清晰.聚焦反演可以得到更为清晰的界面,但层界面位置没有L0.8范数反演刻画的准确.基于小波变换的稀疏约束反演具有反演速度快且易于收敛的特点,同时也可以获得陡变的反演结果,层界面反映准确,但反演结果中会出现局部冗余结构,这是因为Haar小波波形的特点,在对分边界处的小波系数较少,导致对分边界处约束较弱,因此在此处容易产生冗余结构.在计算速度方面,由于Occam反演每次迭代需要多次正演计算来搜索拉格朗日乘子,因此耗时较大;同时对于贝叶斯反演等求解全局最小的反演方法,在搜索过程中需要多次计算雅克比矩阵会产生巨大的计算量,所以非常耗时,因此本文提出的几种反演方法具有计算速度快的优点.
3 实测数据反演本节我们对加拿大Alberta某测区实测航空电磁数据进行反演.数据由CGG公司使用HELITEM MULTIPULSE系统采集.根据Chen等(2015),测区位于Fort McMurray矿区北东方向60 km, 测区表层为35 m中等导电的冰川沉积,中间夹杂良导的粘土层.冰川沉积下伏厚度为35~40 m相对高阻的Grand Rapids砂岩层,其下部Clearwater页岩地层为该测区的主要良导层,厚度60~80 m.Clearwater页岩下部为Fort McMurray地层,为油砂储集层.Fort McMurray地层下部为高阻石灰岩基底,底部存在薄盐层.该测区具有电导率变化明显的分层地质背景,因此可以很好地对本文提出的各种算法进行测试.我们选取10320号测线上248个测点,反演测线长度为9960 m.每个测点初始模型是电阻率为50 Ωm的均匀半空间模型.图 7给出上述五种反演方法的反演结果.由图 7a给出的Occam反演结果可以看到完整和光滑连续的低阻层,但层界面位置和层厚度不能清楚地刻画.图 7b给出L1范数反演结果,在地层的边界处反演结果较为尖锐,可以大致确定层界面位置.图 7c给出L0.8范数反演结果,分层现象明显、低阻层完整且连续,低阻层和高阻层界面非常清晰.图 7d给出聚焦反演结果.虽然聚焦反演可以得到非常陡变的层界面,但是由于聚焦反演难以准确地确定层界面位置,反演的低阻层在横向上不连续,地层厚度也不能准确地刻画.图 7e给出基于小波变换的稀疏约束反演结果.受小波波形影响,反演结果为一些连续的块状结构,在测区左侧较好地恢复了层界面,但右侧区域的地层结构没有得到很好地刻画.
本文在对广义模型约束反演进行系统分析的基础上,成功地将Lp范数反演,聚焦反演和基于小波变换的稀疏约束反演应用于时间域航空电磁数据处理中.通过对理论和实测数据反演,并与Occam反演结果对比得出如下结论:
(1) 提出的反演算法,可以准确地刻画地下层界面,同时反演速度快,是时间域航空电磁用于环境工程和地下水勘查时有效的反演手段.
(2) 对比L0.8、L1范数和聚焦反演结果发现,L1范数可以准确反演层界面的深度信息,而聚焦反演可以获得清晰的层界面(电阻率反演准确),但界面位置不够准确,而L0.8范数反演方法既能得到清晰的层界面,同时也可以准确反演层界面的位置.
(3) 基于小波变换的稀疏约束反演具有高效和高分辨率的特点,可以准确地确定地下电性分界面,但受小波波形特点的影响,在对分边界处有时会产生冗余结构.如何消除对分边界处的冗余结构将是我们未来重要研究方向.
致谢 感谢吉林大学电磁研究团队成员在本文算法研究、模型试验及论文准备过程中给予的帮助,同时特别感谢各位审稿人对本文提出的宝贵意见.
Auken E, Christiansen A V. 2004. Layered and laterally constrained 2D inversion of resistivity data. Geophysics, 69(3): 752-761. DOI:10.1190/1.1759461 |
Auken E, Christiansen A V, Jacobsen L H, et al. 2008. A resolution study of buried valleys using laterally constrained inversion of TEM data. Journal of Applied Geophysics, 65(1): 10-20. DOI:10.1016/j.jappgeo.2008.03.003 |
Cai J, Qi Y F, Yin C C. 2014. Weighted Laterally-constrained inversion of frequency-domain airborne EM data. Chinese J. Geophys. (in Chinese), 57(3): 953-960. DOI:10.6038/cjg20140324 |
Chen T Y, Hodges G, Miles P. 2015. MULTIPULSE-high resolution and high power in one TDEM system. Exploration Geophysics, 46(1): 49-57. DOI:10.1071/EG14027 |
Chiao L Y, Kuo B Y. 2001. Multiscale seismic tomography. Geophysical Journal International, 145(2): 517-527. DOI:10.1046/j.0956-540x.2001.01403.x |
Constable S C, Parker R L, Constable C G. 1987. Occam's inversion:A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3): 289-300. DOI:10.1190/1.1442303 |
Donoho D L. 2006. For most large underdetermined systems oflinear equations the minimal L1-norm solution is also the sparsest solution. Communications on Pure and Applied Mathematics, 59(6): 797-829. DOI:10.1002/(ISSN)1097-0312 |
Ekblom H. 1987. The L1-estimate as limiting case of an Lp-or Huber-estimate.//Doge Y. Statistical Data Analysis based on the L1-Norm and Related Methods. Amsterdam: Elsevier, 109-116.
|
Farquharson C G. 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversion. Geophysics, 73(1): K1-K9. |
Farquharson C G, Oldenburg D W. 1998. Non-linear inversion using general measures of data misfit and model structure. Geophysical Journal International, 134(1): 213-227. DOI:10.1046/j.1365-246x.1998.00555.x |
Farquharson C G, Oldenburg D W. 2003. Constructing piece-wise-constant models using general measures in non-linear, minimum-structure inversion algorithms.//6th International Symposium, Society of Exploration Geophysicists of Japan, Expanded Abstracts, 240-243.
|
Last B J, Kubik K. 1983. Compact gravity inversion. Geophysics, 48(6): 713-721. DOI:10.1190/1.1441501 |
Loke M H, Acworth I, Dahlin T. 2003. A comparison of smooth and blocky inversion methods in 2D electrical imaging surveys. Exploration Geophysics, 34(3): 182-187. DOI:10.1071/EG03182 |
Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596 |
Portniaguine O, Zhdanov M S. 2002. 3-D magnetic inversion with data compression and image focusing. Geophysics, 67(5): 1532-1541. DOI:10.1190/1.1512749 |
Press W H, Teukolsky S A, Vetterling W T, et al. 1996. Numerical Recipes in Fortran 90:The Art of Parallel Scientific Computing. Cambridge: Cambridge University Press.
|
Xu Y X, Wang J Y. 1998. A multiresolution inversion of one-dimensional magnetotelluric data. Acta Geophysica Sinica (in Chinese), 41(5): 704-711. |
Yin C C, Hodges G. 2007. Simulated annealing for airborne EM inversion. Geophysics, 72(4): F189-F195. DOI:10.1190/1.2736195 |
Yin C C, Qi Y F, Liu Y H, et al. 2014. Trans-dimensional Bayesian inversion of frequency-domain airborne EM data. Chinese J.Geophys. (in Chinese), 57(9): 2971-2980. DOI:10.6038/cjg20140922 |
Yin C C, Zhang B, Liu Y H, et al. 2015. Review on airborne EM technology and developments. Chinese J. Geophys. (in Chinese), 58(8): 2637-2653. DOI:10.6038/cjg20150804 |
Zhdanov M S. 2009. New advances in regularized inversion of gravity and electromagnetic data. Geophysical Prospecting, 57(4): 463-478. DOI:10.1111/gpr.2009.57.issue-4 |
蔡晶, 齐彦福, 殷长春. 2014. 频率域航空电磁数据的加权横向约束反演. 地球物理学报, 57(3): 953-960. DOI:10.6038/cjg20140324 |
徐义贤, 王家映. 1998. 大地电磁的多尺度反演. 地球物理学报, 41(5): 704-711. DOI:10.3321/j.issn:0001-5733.1998.05.014 |
殷长春, 齐彦福, 刘云鹤, 等. 2014. 频率域航空电磁数据变维数贝叶斯反演研究. 地球物理学报, 57(9): 2971-2980. DOI:10.6038/cjg20140922 |
殷长春, 张博, 刘云鹤, 等. 2015. 航空电磁勘查技术发展现状及展望. 地球物理学报, 58(8): 2637-2653. DOI:10.6038/cjg20150804 |