1、一种利用特征值性质的 MT 阻尼最小二乘反演 唐荣江 王绪本 甘露 成都理工大学地球物理学院 地球勘探与信息技术教育部重点实验室 摘 要： 阻尼最小二乘反演的方法总会得到一个最终迭代方程, 该方程系数矩阵的特征值对反演结果有重要影响:小特征值相对大特征值对模型修正量的贡献更大, 能够得到更高的分辨率, 但同时更加不稳定。为此, 以一维层状介质的大地电磁阻尼最小二乘反演为基础, 采用特定单位的电阻率参数和层厚参数进行运算, 同时改进阻尼因子的自适应调节方法, 可以巧妙地利用小特征值能提高模型空间分辨率的特点, 得到一种改进的阻尼最小二乘反演方法, 并与传统阻尼最小二乘法进行了对比。合成数据测试结
2、果表明, 该方法能够得到电性陡变的反演模型, 弥补了平滑模型反演不能有效识别电性突变界面的不足, 对电性变化较为剧烈的地电结构的识别具有较好效果, 分辨率有一定提高。关键词： 阻尼最小二乘反演; 大地电磁; 特征值; 阻尼因子; 分辨率; 突变界面; 作者简介：唐荣江 (1992) , 男, 硕士在读, 主要从事电磁地球物理正、反演方法研究。收稿日期：2016-11-08基金：国家自然科学基金项目 (41674078) 资助A damped least square inversion for MT utilizing eigenvalue propertyTANG Rongjiang WAN
3、G Xuben GAN Lu College of Geophysics, Chengdu University of Technology; Abstract： 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 e
4、igenvalues 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 eige
5、nvalues 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 elect
6、rical 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, enabl
7、ing higher resolution.Keyword： damped least square inversion; magnetotelluric (MT) ; eigenvalue; damping factor; resolution; electrical abrupt interface; Received： 2016-11-08阻尼最小二乘法最早由 LEVENBERG 提出, 1963 年 MARQUARDT1对其做了完善, 1971 年 FLETCHER2做了进一步改进, 改进的阻尼最小二乘法是目前求解无约束最优化问题最优秀的算法之一。1988 年陈钟琦3讨论了高阻尼因子的
8、影响并提出了改进方案。传统的大地电磁反演方法大都是基于最大平滑理论, 反演出的地电图像都是渐变的, 不能准确反映出电性突变界面的位置。例如 CON-STABLE 等4提出的 Occam 反演, RODI 等5提出的非线性共轭梯度反演 (nonlinear conjugate gradient inversion, NLCG) 都是光滑模型反演。最小二乘反演将极小值问题转化成求解线性方程组, MARQUARDT1提出阻尼最小二乘法并改进了方程组, 增大了矩阵的特征值, 降低了条件数, 改善了收敛性。奇异值分解法仍是目前求病态解线性方程组最有力的工具之一, 并在近几年发展了基于奇异值分解法的地震数
9、据随机噪声压制方法和子波提取方法6-7。根据奇异值分解理论, WIGGINS8曾用截断一些小奇异值的方法来减少解, 以取得与阻尼最小二乘法类似的效果, 但是应当截除什么样的奇异值是奇异值截断法中的主要困难, 且截断小奇异值后分辨率也随之降低。李平等9探讨了奇异值分解中广义逆解的可靠性估计, 给出了更加合理的奇异值修改方案。阻尼最小二乘法最终迭代矩阵为满秩对称方阵, 相对于雅可比矩阵来说形式更简易, 且雅可比矩阵的奇异值和迭代方阵的特征值对反演结果具有相同影响。阻尼因子的大小直接影响矩阵特征值的大小, 小特征值与高分辨率对应, 大特征值与较小的方差和较好的收敛性对应。如果能有效利用和改进小特征值
10、, 同时注意大特征值以保证方程的收敛性, 是否能在一定程度上提高阻尼最小二乘法的分辨率, 或针对不同的地质模型采取不同的特征值变化方案以达到约束反演的目的?通常情况下, 大地电磁一维反演的分辨率要高于二维、三维反演的分辨率, 马奎特一维反演中层厚和电阻率要同时参与运算, 这为提高一维反演的分辨率和优化反演结果提供了更大的改进空间。据此, 本文以一维层状模型的大地电磁阻尼最小二乘反演为基础, 进行特征值分析, 探讨其如何影响模型空间, 据此提出改进方法:不直接修改特征值, 而是使用特定单位的反演参数和改变阻尼因子的自适应方案以间接达到改变特征值的目的, 从而优化反演结果。理论数据实验表明, 改进
11、后的反演方法能够有效识别电性突变界面, 弥补了光滑模型反演的缺陷, 一定程度上提高了大地电磁的分辨率。1 反演基本原理大地电磁阻尼最小二乘法的目标函数为:式中:d 为观测数据; 为阻尼因子;m= 1 j N, h1hjhN-1为模型参数, 其中, j为第 j 层电阻率, h j为第 j 层厚度, N 为反演层数;A (m) 为正演算子;W d和 Wm分别为数据空间和模型空间加权矩阵。为了将非线性问题线性化, 考虑:用泰勒公式对 A (m1) 进行一阶展开, F 为 A 的弗雷歇导数, 若 Wd和 Wm为单位矩阵, A (m 1) 表示为 , 目标函数变为:求目标函数对 m 的导数:根据最小化条
12、件要求导数为 0, 可得:通过解方程 (5) , 可得到 m, 多次迭代之后可得到满足精度的反演解。在大地电磁反演中, 为了增强方程的稳定性, 陈乐寿等10对目标函数做了改进:对数据空间取对数, 并对模型空间取比例。据此, 方程 (1) 变为:式中:lnA (m) -lnd为数据拟合泛函或拟合差, 后文用 1 表示;为稳定泛函, 用 2 表示。对 (6) 式求导, 取其导数为零, 得最终的迭代方程:式中: , i 为频点个数循环量, j 为模型参数循环量;m 在这里为当前迭代的模型参数。在常规的阻尼最小二乘法中, 为了避免电阻率本身作为模型参数数量级差异太大造成反演不稳定, 反演对电阻率和层厚
13、作为模型参数时取其对数:那么雅可比矩阵中:根据上述理论, 常规阻尼最小二乘反演具体实现步骤10如下。1) 设置初始模型 m, 初始阻尼因子 (这里设置为 1) , 加权系数 k (通常设置为 10) , 给定数据拟合误差的阈值 E。4) 判断条件 1 (m+m) 1 (m) 时, 才被迫取较大的阻尼因子, 使迭代稳定收敛。当 给定后, 逐步尝试增加或减少。若采用本文方法的反演参数, 层厚数值较小, 在 较小时方程稳定性降低, 层厚参数极易变为负值, 使反演失败。此外, 根据2.2 节, 为了使 与 h 尽快接近以达到层厚改变的目的, 必须保证阻尼因子的快速下降。基于以上目的, 我们改进阻尼因子
14、的自适应变化方案, 反演的步骤 4) 变为:4) m=m+m, 判断各层的 mi是否大于 0, 如果是, 则 =k, 返回步骤 2) , 否则 =k, 返回步骤 3) 。两种方法的不同之处在于判断 减少或增大的条件:传统方法是根据 1和 1的大小来确定, 而本文方法采用判断反演的参数是否为正来确定。由于传统方法每次反演迭代得到的解需要取指数, 指数值则始终大于 0, 所以不用考虑参数变为负的问题。为了说明本文方法能够保证方程的稳定性和收敛性以及阻尼因子的快速降低, 我们设计一个简单 3 层模型:电阻率为 100m 的低阻层埋深 1.0km, 厚度为 0.5km, 背景电阻率为 1 000m。正
15、演数据在 0.0016 000.000Hz 采用 30 个计算频点, 无随机噪声。反演深度 3km, 建立 15 层模型, 迭代 20 次, 阻尼因子初始值设置 10, 采用不同的阻尼因子变化方案和模型参数得到拟合差结果如表 1 所示。由表 1 可知: (1) 本文方法的 在第 7 次迭代之后就远远低于传统方法的 , 说明新的收敛条件可以保证 更快的下降速度; (2) 两种方法的拟合差均逐步减少, 本文方法的拟合差没有传统方法的小, 但阻尼因子更低, 导致特征值更低, 观测数据的微小变化将导致模型参数的更大变化; (3) 本文方法中第 2 次迭代到第 3 次迭代拟合差不降反升, 且后期震荡变化
16、, 整体为降低的趋势, 说明阻尼因子起着调节特征值和稳定性的作用, 当层厚修正量过大使层厚出现负值时, 则提高阻尼因子, 降低层厚修正量, 保证方程稳定性;若层厚为正值, 则降低阻尼因子, 保证分辨率。以上结论说明通过判断模型参数是否为负值来保证反演的稳定性和收敛性是合理的。表 1 不同迭代次数下阻尼因子和拟合差变化规律 下载原表 2.4 分辨率提高的原因前文的理论说明:采用特定的反演参数单位保证了反演迭代方程具有更小的特征值, 使得层厚参数相对于电阻率有更高的修正量;本文方法保证了反演的稳定性和收敛性, 要说明的是, 新的改进能够达到更高的分辨率, 迭代后期层厚参数修正量必须以负值为主, 使
17、平滑的反演模型得到压缩, 突出异常。但使用定量分析证明模型修正量为负值较为困难, 可作定性分析:在保证数据拟合差极小的情况下, 相对于增加层厚, 减少层厚来保证数据的高度拟合明显要容易得多。为了说明其正确性, 这里同样采用 2.3 节中的模型和参数, 使用本文改进方法对不同迭代次数的模型修正量成图, 并对比两种方法的反演结果。在图 1 中, 图 1a 和图 1b 的纵坐标分别为模型修正量 m 和模型修正量 m 与模型参数 m (包括电阻率 和层厚 h) 之比, 横坐标为反演参数序号, 图 1a 和图 1b 中的反演参数序号 115 对应第 1 层到 15 层的电阻率, 1629 对应第 1 层
18、到 14 层的层厚。由图可见, 图 1a 中 h 远小于 , 且两者在迭代后期都趋于 0;图 1b 中 h/h 在早期总体小于 /, 随后逐渐增大甚至大于 /, 两者在后期都减小, 但 h/h 总体仍大于 /, 且主要为负值。应当注意, 引起层厚变化率大于电阻率变化率的并不是模型修正量本身, 而是修正量与模型参数的比例。图 2a 和图 2b 分别为对单个低阻层模型采用传统方法和本文方法进行 20 次迭代时的反演结果, 可以看出, 相对于本文方法, 传统方法的反演结果较为平缓, 不能准确反映电性突变界面的地下结构。从图 2b 可以看出, 电性过度界面节点较为密集, 说明该处的层厚受到压缩, 在一
19、定程度上提高了分辨率, 同时表明, 利用层厚的小特征值可以提高模型空间的分辨率, 而传统方法拟合差更低。该结论与小特征值对应高分辨率, 大特征值意味着以牺牲一些分辨率为代价来换取较小的方差的结论一致。图 1 不同迭代次数下采用本文反演方法的模型修正量 下载原图a 模型修正量;b 模型修正量与模型参数比值图 2 单个低阻层模型传统方法 (a) 与本文方法 (b) 进行 20 次迭代的反演结果 下载原图本文方法通过压缩电性界面突变区层厚来突出异常, 所以对于平滑模型的反演, 本文方法得出的结果往往与实际情况不吻合, 此时选用传统阻尼最小二乘反演能取得更为理想的结果。所以, 在实际情况中, 针对不同
20、的反演目标, 应当选择合适的反演方法来解决地质问题。2.5 初始模型和解方程组从本质上来说, 本文方法仍然是阻尼最小二乘法的延续, 所以本文方法和传统阻尼最小二乘反演法的初始模型、方程组解法、收敛条件等的选择原则一致。本文解方程组采用共轭梯度法12, 该方法可以快速、稳定求解对称矩阵的线性代数方程组;收敛条件上可以采用设定阈值或者固定迭代次数的方式, 反演深度和加权系数在各理论模型的反演中均有说明。初始模型建立可以采用均匀半空间和 Bostic 变换, 本文采用类似于均匀半空间的改进方法:以大地电磁最高频视电阻率当作初始模型第 1 层电阻率, 将最低频视电阻率当作初始模型最后一层电阻率, 其它
21、层的电阻率位于两者之间渐变分布, 各层厚度采用等比增长。3 合成数据测试及分析设计了 2 个理论模型来验证本文算法的有效性。图 3 的 6 层理论模型与文献13中的模型一致, 图 4 的 8 层理论模型与文献14中的模型一致, 无随机噪声, 反演模型设定为 35 层, 6 层模型反演深度 15km, 8 层模型反演深度为 40km, 阻尼因子初始值都为 10, 加权系数 k 为 10, 迭代次数为 20 次。图 3 采用本文方法和传统阻尼最小二乘法的 6 层模型反演结果对比 下载原图a 传统方法;b 本文方法;c 拟合曲线图 4 采用本文方法和传统阻尼最小二乘法的 8 层模型反演结果对比 下载
22、原图a 传统方法;b 本文方法;c 拟合曲线图 3 和图 4 分别显示了采用本文方法、传统阻尼最小二乘法的 6 层模型和 8 层模型的反演结果。由图 3 和图 4 可见, 采用传统方法能基本反演出各层地电模型, 吻合程度相对较低, 电性突变界面的反演层过渡相对平缓, 不能准确识别突变界面的具体位置 (图 3a, 图 4a) ;采用本文方法不但能有效反演出各层地电模型, 而且在很大程度上与理论模型吻合, 能有效识别电性突变界面, 相对传统方法来说, 分辨率有一定提高 (图 3b, 图 4b) ;两种算法都能较好地拟合原始数据, 虽然传统方法具有更高的拟合精度, 但从图像上无法直接看出, 这说明可
23、以接受本文算法的拟合精度, 且一定程度上牺牲了拟合精度来提高分辨率是值得的 (图 3c, 图 4c) 。4 结束语本文以迭代方程特征值的性质为理论基础, 进行反演参数和阻尼因子自适应方法的相关改进, 并将改进后的大地电磁一维阻尼最小二乘算法有效地应用于解决电性突变地质结构的反演问题。研究表明, 小特征值对应大的模型修正量和更高的分辨率, 大特征值对应小的模型修正量和更小的拟合差;较大的特征值差异导致较大模型修正量差异, 较小的特征值差异导致较小模型修正量差异。使用了数量级悬殊的电阻率 (m) 和层厚 (km) 直接反演, 可以导致对应的特征值差异。阻尼因子较大时, 较大, h 较小, 随着 的
24、降低, 和h 逐渐接近, 数量级较小的层厚更易发生改变, 且修正量主要为负值。同时, 采取判断模型参数是否全为正来修正阻尼因子 , 力求在保证反演收敛情况下, 最快速度下降阻尼因子, 达到使 和 h 逐渐接近的目的。一维层状理论模型的反演结果表明:改进后的反演方法能够有效识别电性突变界面, 针对地下电性突变的结构, 在一定程度上提高了大地电磁测深的分辨率, 且该方法可以拓展应用到可控源电磁法的反演上, 但是对于地下电性平滑的界面, 采用传统方法反演则更有效。附录 A以下所有矩阵的特征值满足降序排列 1 2 n。定理 1:若 A, B 为 n 阶 Hermite 矩阵, A, B 的特征值分别为
25、 i (A) , i (B) , 则对于 i=1, 2, , n, 有:详细证明过程请参见文献11, 根据定理 1, 这里得出如下推论。推论:若 A, B 为 2 N-1 阶实对称矩阵, B 满足前 N 个特征值远大于后 N-1 个特征值, 且远大于 A 的特征值, 则有: i (A+B) 的前 N 个特征值 s1 (B) 远大于后 N-1 个特征值 s2 (B) ;且随着 B 矩阵所有特征值成比例逐渐减小, i (A+B) 的前 N 个特征值与后 N-1 个特征值的差逐渐变小。证明:当 i1, N时, 根据定理 1, 不等式左端最大值条件为 2 Nr+s3 N-1, r (A) + s (B
26、) 要取得最大值, B 的特征值必在前 N 个中取值;不等式右端最小值条件为 1r+sN+1, r (A) + s (B) 要取得最小值, 由于 r1, 必有 sN, 则 B 的特征值只能在前 N 个中取值。当 iN+1, 2 N-1时, 根据定理 1, 不等式左端最大值条件为 3 Nr+s4 N-2, r (A) + s (B) 要取得最大值, 由于 r2 N-1, 必有 sN+1, 则 B 的特征值只能在后 N-1 个中取值;不等式右端最小值条件为 N+2r+s2 N, r (A) + s (B) 要取得最小值, 由于 r1, 必有 sN+1, 则 B 的特征值只能在后 N-1 个中取值。
27、综上可知, i (A+B) 的前 N 个特征值的上、下界均包含 s1 (B) , 且1s 1N; i (A+B) 的后 N-1 个特征值的上、下界均包含 s2 (B) , 不含 s1 (B) , 且 N+1s 22 N-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 个特征值的差逐渐变小。参考文献
