一种利用特征值性质的MT阻尼最小二乘反演

唐荣江,王绪本,2,甘 露

(1.成都理工大学地球物理学院,四川成都610059;2.地球勘探与信息技术教育部重点实验室,四川成都610059)

唐荣江,王绪本,甘露.一种利用特征值性质的MT阻尼最小二乘反演[J].石油物探,2017,56(6):-904

TANG Rongjiang,WANG Xuben,GAN Lu.A damped least square inversion for MT utilizing eigenvalue property[J].Geophysical Prospecting for Petroleum,2017,56(6):-904

一种利用特征值性质的MT阻尼最小二乘反演

唐荣江1,王绪本1,2,甘 露1

(1.成都理工大学地球物理学院,四川成都610059;2.地球勘探与信息技术教育部重点实验室,四川成都610059)

阻尼最小二乘反演的方法总会得到一个最终迭代方程,该方程系数矩阵的特征值对反演结果有重要影响:小特征值相对大特征值对模型修正量的贡献更大,能够得到更高的分辨率,但同时更加不稳定。为此,以一维层状介质的大地电磁阻尼最小二乘反演为基础,采用特定单位的电阻率参数和层厚参数进行运算,同时改进阻尼因子的自适应调节方法,可以巧妙地利用小特征值能提高模型空间分辨率的特点,得到一种改进的阻尼最小二乘反演方法,并与传统阻尼最小二乘法进行了对比。合成数据测试结果表明,该方法能够得到电性陡变的反演模型,弥补了平滑模型反演不能有效识别电性突变界面的不足,对电性变化较为剧烈的地电结构的识别具有较好效果,分辨率有一定提高。

阻尼最小二乘反演;大地电磁;特征值;阻尼因子;分辨率;突变界面

阻尼最小二乘法最早由LEVENBERG提出,1963年MARQUARDT[1]对其做了完善,1971年FLETCHER[2]做了进一步改进,改进的阻尼最小二乘法是目前求解无约束最优化问题最优秀的算法之一。1988年陈钟琦[3]讨论了高阻尼因子的影响并提出了改进方案。传统的大地电磁反演方法大都是基于最大平滑理论,反演出的地电图像都是渐变的,不能准确反映出电性突变界面的位置。例如CONSTABLE等[4]提出的Occam反演,RODI等[5]提出的非线性共轭梯度反演(nonlinear conjugate gradient inversion,NLCG)都是光滑模型反演。

最小二乘反演将极小值问题转化成求解线性方程组,MARQUARDT[1]提出阻尼最小二乘法并改进了方程组,增大了矩阵的特征值,降低了条件数,改善了收敛性。奇异值分解法仍是目前求病态解线性方程组最有力的工具之一,并在近几年发展了基于奇异值分解法的地震数据随机噪声压制方法和子波提取方法[6-7]。根据奇异值分解理论,WIGGINS[8]曾用截断一些小奇异值的方法来减少解,以取得与阻尼最小二乘法类似的效果,但是应当截除什么样的奇异值是奇异值截断法中的主要困难,且截断小奇异值后分辨率也随之降低。李平等[9]探讨了奇异值分解中广义逆解的可靠性估计,给出了更加合理的奇异值修改方案。

阻尼最小二乘法最终迭代矩阵为满秩对称方阵,相对于雅可比矩阵来说形式更简易,且雅可比矩阵的奇异值和迭代方阵的特征值对反演结果具有相同影响。阻尼因子的大小直接影响矩阵特征值的大小,小特征值与高分辨率对应,大特征值与较小的方差和较好的收敛性对应。如果能有效利用和改进小特征值,同时注意大特征值以保证方程的收敛性,是否能在一定程度上提高阻尼最小二乘法的分辨率,或针对不同的地质模型采取不同的特征值变化方案以达到约束反演的目的？

通常情况下,大地电磁一维反演的分辨率要高于二维、三维反演的分辨率,马奎特一维反演中层厚和电阻率要同时参与运算,这为提高一维反演的分辨率和优化反演结果提供了更大的改进空间。据此,本文以一维层状模型的大地电磁阻尼最小二乘反演为基础,进行特征值分析,探讨其如何影响模型空间,据此提出改进方法:不直接修改特征值,而是使用特定单位的反演参数和改变阻尼因子的自适应方案以间接达到改变特征值的目的,从而优化反演结果。理论数据实验表明,改进后的反演方法能够有效识别电性突变界面,弥补了光滑模型反演的缺陷,一定程度上提高了大地电磁的分辨率。

1 反演基本原理

大地电磁阻尼最小二乘法的目标函数为:

(1)

式中:d为观测数据;α为阻尼因子;m=[ρ1…ρj…ρN,h1…hj…hN-1]为模型参数,其中,ρj为第j层电阻率,hj为第j层厚度,N为反演层数;A(m)为正演算子;Wd和Wm分别为数据空间和模型空间加权矩阵。为了将非线性问题线性化,考虑:

(2)

用泰勒公式对A(m1)进行一阶展开,F为A的弗雷歇导数,若Wd和Wm为单位矩阵,A(m1)表示为A(m0)+Fm0Δm,目标函数变为:

(3)

求目标函数对ΔmT的导数:

(4)

根据最小化条件要求导数为0,可得:

(5)

通过解方程(5),可得到Δm,多次迭代之后可得到满足精度的反演解。在大地电磁反演中,为了增强方程的稳定性,陈乐寿等[10]对目标函数做了改进:对数据空间取对数,并对模型空间取比例。据此,方程(1)变为:

(6)

(7)

在常规的阻尼最小二乘法中,为了避免电阻率本身作为模型参数数量级差异太大造成反演不稳定,反演对电阻率和层厚作为模型参数时取其对数:

(8)

那么雅可比矩阵中:

根据上述理论,常规阻尼最小二乘反演具体实现步骤[10]如下。

1) 设置初始模型m,初始阻尼因子α(这里设置为1),加权系数k(通常设置为10),给定数据拟合误差的阈值E。

4) 判断条件Φ1(m+Δm)<Φ1(m)是否满足,如果是,m=m+Δm,α=αk-1,执行步骤5),否则α=αk,返回步骤3)。

5) 判断条件‖lnA(m+Δm)-lnd‖

2 特征值的利用与反演方法的改进

2.1 特征值分析

反演迭代方程的特征值对反演结果具有重要影响,下面分析特征值的大小如何影响模型参数。

(11)

(12)

Λ的具体表达式为:

(13)

式中:λi为H+αD的特征值。根据(7)式,模型修正量Δm的求和形式为:

(14)

(15)

式中:N为层厚;M为测点数。由(14)式可知,特征值λi越小,对模型参数Δmj的影响越大,同时也使方程更加不稳定。

2.2 反演参数的改进

反演参数的改进为:避免取对数,直接使用以千米为单位的层厚和欧姆米为单位的电阻率进行反演。这种改进巧妙地利用了小特征值对模型修正量贡献更大的特点,使得层厚修正量比例大于电阻率修正量比例,以达到持续改变反演模型的目的。

2.3 阻尼因子的改进与反演的收敛性

实际计算中,阻尼因子α不能取太大,否则步长太小,收敛太慢。仅当Φ1(m+Δm)>Φ1(m)时,才被迫取较大的阻尼因子,使迭代稳定收敛。当α(0)给定后,逐步尝试增加或减少。若采用本文方法的反演参数,层厚数值较小,在α较小时方程稳定性降低,层厚参数极易变为负值,使反演失败。此外,根据2.2节,为了使Δρ与Δh尽快接近以达到层厚改变的目的,必须保证阻尼因子的快速下降。基于以上目的,我们改进阻尼因子的自适应变化方案,反演的步骤4)变为:

4)m=m+Δm,判断各层的mi是否大于0,如果是,则α=αk-1,返回步骤2),否则α=αk,返回步骤3)。

由表1可知:①本文方法的α在第7次迭代之后就远远低于传统方法的α,说明新的收敛条件可以保证α更快的下降速度;②两种方法的拟合差均逐步减少,本文方法的拟合差没有传统方法的小,但阻尼因子更低,导致特征值更低,观测数据的微小变化将导致模型参数的更大变化;③本文方法中第2次迭代到第3次迭代拟合差不降反升,且后期震荡变化,整体为降低的趋势,说明阻尼因子起着调节特征值和稳定

表1 不同迭代次数下阻尼因子和拟合差变化规律

性的作用,当层厚修正量过大使层厚出现负值时,则提高阻尼因子,降低层厚修正量,保证方程稳定性;若层厚为正值,则降低阻尼因子,保证分辨率。以上结论说明通过判断模型参数是否为负值来保证反演的稳定性和收敛性是合理的。

2.4 分辨率提高的原因

前文的理论说明:采用特定的反演参数单位保证了反演迭代方程具有更小的特征值,使得层厚参数相对于电阻率有更高的修正量;本文方法保证了反演的稳定性和收敛性,要说明的是,新的改进能够达到更高的分辨率,迭代后期层厚参数修正量必须以负值为主,使平滑的反演模型得到压缩,突出异常。但使用定量分析证明模型修正量为负值较为困难,可作定性分析:在保证数据拟合差极小的情况下,相对于增加层厚,减少层厚来保证数据的高度拟合明显要容易得多。为了说明其正确性,这里同样采用2.3节中的模

型和参数,使用本文改进方法对不同迭代次数的模型修正量成图,并对比两种方法的反演结果。

在图1中,图1a和图1b的纵坐标分别为模型修正量Δm和模型修正量Δm与模型参数m(包括电阻率ρ和层厚h)之比,横坐标为反演参数序号,图1a和图1b 中的反演参数序号1～15对应第1层到15层的电阻率,16～29对应第1层到14层的层厚。由图可见,图1a中Δh远小于Δρ,且两者在迭代后期都趋于0;图1b中Δh/h在早期总体小于Δρ/ρ,随后逐渐增大甚至大于Δρ/ρ,两者在后期都减小,但Δh/h总体仍大于Δρ/ρ,且主要为负值。应当注意,引起层厚变化率大于电阻率变化率的并不是模型修正量本身,而是修正量与模型参数的比例。

图2a和图2b分别为对单个低阻层模型采用传统方法和本文方法进行20次迭代时的反演结果,可以看出,相对于本文方法,传统方法的反演结果较为平缓,不能准确反映电性突变界面的地下结构。从图2b 可以看出,电性过度界面节点较为密集,说明该处的层厚受到压缩,在一定程度上提高了分辨率,同时表明,利用层厚的小特征值可以提高模型空间的分辨率,而传统方法拟合差更低。该结论与小特征值对应高分辨率,大特征值意味着以牺牲一些分辨率为代价来换取较小的方差的结论一致。

图1 不同迭代次数下采用本文反演方法的模型修正量 a 模型修正量; b 模型修正量与模型参数比值

图2 单个低阻层模型传统方法(a)与本文方法(b)进行20次迭代的反演结果

本文方法通过压缩电性界面突变区层厚来突出异常,所以对于平滑模型的反演,本文方法得出的结果往往与实际情况不吻合,此时选用传统阻尼最小二乘反演能取得更为理想的结果。所以,在实际情况中,针对不同的反演目标,应当选择合适的反演方法来解决地质问题。

2.5 初始模型和解方程组

从本质上来说,本文方法仍然是阻尼最小二乘法的延续,所以本文方法和传统阻尼最小二乘反演法的初始模型、方程组解法、收敛条件等的选择原则一致。本文解方程组采用共轭梯度法[12],该方法可以快速、稳定求解对称矩阵的线性代数方程组;收敛条件上可以采用设定阈值或者固定迭代次数的方式,反演深度和加权系数在各理论模型的反演中均有说明。

初始模型建立可以采用均匀半空间和Bostic变换,本文采用类似于均匀半空间的改进方法:以大地电磁最高频视电阻率当作初始模型第1层电阻率,将最低频视电阻率当作初始模型最后一层电阻率,其它层的电阻率位于两者之间渐变分布,各层厚度采用等比增长。

3 合成数据测试及分析

设计了2个理论模型来验证本文算法的有效性。图3的6层理论模型与文献[13]中的模型一致,图4的8层理论模型与文献[14]中的模型一致,无随机噪声,反演模型设定为35层,6层模型反演深度15km,8层模型反演深度为40km,阻尼因子初始值都为10,加权系数k为10,迭代次数为20次。

图3 采用本文方法和传统阻尼最小二乘法的6层模型反演结果对比 a 传统方法; b 本文方法; c 拟合曲线

图4 采用本文方法和传统阻尼最小二乘法的8层模型反演结果对比 a 传统方法; b 本文方法; c 拟合曲线

图3和图4分别显示了采用本文方法、传统阻尼最小二乘法的6层模型和8层模型的反演结果。由图3和图4可见,采用传统方法能基本反演出各层地电模型,吻合程度相对较低,电性突变界面的反演层过渡相对平缓,不能准确识别突变界面的具体位置(图3a,图4a);采用本文方法不但能有效反演出各层地电模型,而且在很大程度上与理论模型吻合,能有效识别电性突变界面,相对传统方法来说,分辨率有一定提高(图3b,图4b);两种算法都能较好地拟合原始数据,虽然传统方法具有更高的拟合精度,但从图像上无法直接看出,这说明可以接受本文算法的拟合精度,且一定程度上牺牲了拟合精度来提高分辨率是值得的(图3c,图4c)。

4 结束语

本文以迭代方程特征值的性质为理论基础,进行反演参数和阻尼因子自适应方法的相关改进,并将改进后的大地电磁一维阻尼最小二乘算法有效地应用于解决电性突变地质结构的反演问题。研究表明,小特征值对应大的模型修正量和更高的分辨率,大特征值对应小的模型修正量和更小的拟合差;较大的特征值差异导致较大模型修正量差异,较小的特征值差异导致较小模型修正量差异。使用了数量级悬殊的电阻率(Ω·m)和层厚(km)直接反演,可以导致对应的特征值差异。阻尼因子较大时,Δρ较大,Δh较小,随着α的降低,Δρ和Δh逐渐接近,数量级较小的层厚更易发生改变,且修正量主要为负值。同时,采取判断模型参数是否全为正来修正阻尼因子α,力求在保证反演收敛情况下,最快速度下降阻尼因子,达到使Δρ和Δh逐渐接近的目的。

一维层状理论模型的反演结果表明:改进后的反演方法能够有效识别电性突变界面,针对地下电性突变的结构,在一定程度上提高了大地电磁测深的分辨率,且该方法可以拓展应用到可控源电磁法的反演上,但是对于地下电性平滑的界面,采用传统方法反演则更有效。

[1] MARQUARDT D W.An algorithm for least-squares estimation of nonlinear parameters[J].Journal of the Society for Industrial and Applied Mathematics,1963,11(2):431-441

[2] FLETCHER R.A modified Marquardt subroutine for non-linear least squares:UKAEA Research Group Report AERF R-6799[R].UK:Atomic Energy Research Establishment,1971

[3] 陈钟琦.高阻尼因子对阻尼最小二乘法效果的影响和克服[J].现代地质,1988,2(2):255-264

CHEN Z Q.The influence of the high damping factor to the result of damped least square method and it’s overcome[J].Geoscience,1988,2(2):255-264

[4] CONSTABLE S C,PARKER R L,CONSTABLE C G.Occam’s inversion:a practical algorithm for generating smooth models from electromagnetic sounding data[J].Geophysics,1987,52(3):289-300

[5] RODI W,MACKIE R L.Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[J].Geophysics,2001,66(1):174-187

[6] 马继涛,王建花,刘国昌.基于频率域奇异值分解的地震数据插值去噪方法研究[J].石油物探,2016,55(2):205-213

MA J T,WANG J H,LIU G C.Seismic data noise attenuation and interpolation using singular value decomposition in frequency domain[J].Geophysical Prospecting for Petroleum,2016,55(2):205-213

[7] 王蓉蓉,戴永寿,李闯,等.基于奇异值分解的时变子波提取准确性评价方法[J].石油物探,2015,54(5):531-540

WANG R R,DAI Y S,LI C,et al.An evaluation criterion on the accuracy of time-varying wavelet extraction based on singular value decomposition[J].Geophysical Prospecting for Petroleum,2015,54(5):531-540

[8] WIGGINS R A.The general linear inverse problem:implication of surface waves and free oscillations for earth structure[J].Reviews of Geophysics,1972,10(1):251-285

[9] 李平,王椿镛,徐厚泽,等.地球物理反演中奇异值分解应用的若干问题探讨[J].自然科学进展,2001,11(8):891-896

LI P,WANG C Y,XU H Z,et al.Several problems discussion in the application of singular value decomposition in geophysical inversion[J].Progress in Natural Science,2001,11(8):891-896

[10] 陈乐寿,刘任,壬天生.大地电磁测深资料处理与解释[M].北京:石油工业出版社,1989:1-236

CHEN L S,LIU R,REN T S,et al.Magnetotelluric sounding data processing and interpretation[M].Beijing:Petroleum Industry Publishing House,1989:1-236

[11] 廖安平,刘建州.矩阵论[M].长沙:湖南大学出版社,2005:97-101

LIAO A P,LIU J Z.Theory of matrices[M].Changsha:Hunan University Publishing House,2005:97-101

[12] NOCEDAL J,WRIGHT S.Numerical optimization 2nd [M].The United States:Springer Science & Business Media,2006:101-133

[13] 管贻亮,李予国,胡祖志,等.大地电磁非线性共轭梯度一维反演[J].石油物探,2014,53(6):752-759

GUAN Y L,LI Y G,HU Z Z,et al.Nonlinera conjugate gradients algorithm for 1D magnetotelluric inversion[J].Geophysical Prospecting for Petroleum,2014,53(6):752-759

[14] 陈小斌,赵国泽,汤吉,等.大地电磁自适应正则化反演算法[J].地球物理学报,2005,48(4):937-945

CHEN X B,ZHAO G Z,TANG J,et al.An adaptive regularized inversion algorithm for magnetotelluric data[J].Chinese Journal of Geophysics,2005,48(4):937-945

附录A

以下所有矩阵的特征值满足降序排列λ1>λ2>…>λn。

定理1:若A,B为n阶Hermite矩阵,A,B的特征值分别为λi(A),λi(B),则对于i=1,2,…,n,有:

(A1)

详细证明过程请参见文献[11],根据定理1,这里得出如下推论。

推论:若A,B为2N-1阶实对称矩阵,B满足前N个特征值远大于后N-1个特征值,且远大于A的特征值,则有:λi(A+B)的前N个特征值λs1(B)远大于后N-1个特征值λs2(B);且随着B矩阵所有特征值成比例逐渐减小,λi(A+B)的前N个特征值与后N-1个特征值的差逐渐变小。

证明:

当i∈[1,N]时,根据定理1,不等式左端最大值条件为2N≤r+s≤3N-1,λr(A)+λs(B)要取得最大值,B的特征值必在前N个中取值;不等式右端最小值条件为1≤r+s≤N+1,λr(A)+λs(B)要取得最小值,由于r≥1,必有s≤N,则B的特征值只能在前N个中取值。

当i∈[N+1,2N-1]时,根据定理1,不等式左端最大值条件为3N≤r+s≤4N-2,λr(A)+λs(B)要取得最大值,由于r≤2N-1,必有s≥N+1,则B的特征值只能在后N-1个中取值;不等式右端最小值条件为N+2≤r+s≤2N,λr(A)+λs(B)要取得最小值,由于r≥1,必有s≥N+1,则B的特征值只能在后N-1个中取值。

综上可知,λi(A+B)的前N个特征值的上、下界均包含λs1(B),且1≤s1≤N;λi(A+B)的后N-1个特征值的上、下界均包含λs2(B),不含λs1(B),且N+1≤s2≤2N-1;所以λi(A+B)的前N个特征值远大于后N-1个特征值。随着B矩阵所有特征值成比例逐渐减小,λi(A+B)的前N个特征值的上、下界大幅减小,因为λs(B)在λr(A)+λs(B)中占有主要比重;后N-1个特征值的上、下界小幅度降低,因为减小幅度受λr(A)和λs(B)的比例共同控制,所以λi(A+B)的前N个特征值与后N-1个特征值的差逐渐变小。

(编辑:顾石庆)

AdampedleastsquareinversionforMTutilizingeigenvalueproperty

TANG Rongjiang1,WANG Xuben1,2,GAN Lu1

(1.CollegeofGeophysics,ChengduUniversityofTechnology,Chengdu610059,China;2.KeyLaboratoryofEarthexplorationandInformationTechniquesofMinistryofEducation,Chengdu610059,China)

Damping least squares inversion always involves a final iterative equation,and the eigenvalues of its coefficient matrix have an important effect on the inversion results.Namely,compared to large eigenvalues,small eigenvalues contribute more to model corrections,resulting in higher resolution,but leading to lower stability.An improved damped least square inversion method is proposed that is based on a one-dimensional magnetotelluric inversion scheme.The method cleverly uses the characteristic that small eigenvalues can improve the spatial resolution,which uses the resistivity and layer thickness parameters of a specific unit for inversion,while improving the adaptive adjustment method of the damping factor.Synthetic data test results showed that,the proposed method could invert a model with an electrical abrupt change,which makes up for the shortcoming of smooth model inversion,which cannot effectively identify an electrical abrupt interface.Compared with traditional methods,this method is more effective for the identification of electrical structures with geoelectrical abrupt change,enabling higher resolution.

damped least square inversion,magnetotelluric (MT),eigenvalue,damping factor,resolution,electrical abrupt interface

2016-11-08;改回日期2017-06-12。

唐荣江(1992—),男,硕士在读,主要从事电磁地球物理正、反演方法研究。

国家自然科学基金项目(41674078)资助。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41674078).

P631

A

1000-1441(2017)06-0898-07

10.3969/j.issn.1000-1441.2017.06.016