阻尼最小二乘法最早由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 反演基本原理

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

$ P\left( {\mathit{\boldsymbol{m}},\mathit{\boldsymbol{d}}} \right) = {\left\| {{\mathit{\boldsymbol{W}}_d}\mathit{\boldsymbol{A}}\left( \mathit{\boldsymbol{m}} \right) - {\mathit{\boldsymbol{W}}_d}\mathit{\boldsymbol{d}}} \right\|^2} + \alpha {\left\| {{\mathit{\boldsymbol{W}}_m}\Delta \mathit{\boldsymbol{m}}} \right\|^2} $

(1)

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

$ {\mathit{\boldsymbol{m}}_1} = {\mathit{\boldsymbol{m}}_0} + \Delta \mathit{\boldsymbol{m}} $

(2)

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

$ \begin{array}{l} P\left( {{\mathit{\boldsymbol{m}}_1},\mathit{\boldsymbol{d}}} \right) = {\left[ {\mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{m}}_0}} \right) + {\mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}}\Delta \mathit{\boldsymbol{m}} - \mathit{\boldsymbol{d}}} \right]^{\rm{T}}}\left[ {\mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{m}}_0}} \right) + } \right.\\ \;\;\;\;\;\;\left. {{\mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}}\Delta \mathit{\boldsymbol{m}} - \mathit{\boldsymbol{d}}} \right] + \alpha \Delta {\mathit{\boldsymbol{m}}^{\rm{T}}}\Delta \mathit{\boldsymbol{m}} \end{array} $

(3)

求目标函数对Δm^T的导数:

$ \frac{{\partial P}}{{\partial \Delta {\mathit{\boldsymbol{m}}^{\rm{T}}}}} = \mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}^{\rm{T}}\left[ {\mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{m}}_0}} \right) + {\mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}}\Delta \mathit{\boldsymbol{m}} - \mathit{\boldsymbol{d}}} \right] + \alpha \Delta \mathit{\boldsymbol{m}} $

(4)

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

$ \left( {\mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}^{\rm{T}}{\mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}} + \alpha } \right)\Delta \mathit{\boldsymbol{m}} = \mathit{\boldsymbol{F}}_{{\mathit{\boldsymbol{m}}_0}}^{\rm{T}}\left[ {\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{m}}_0}} \right)} \right.} \right] $

(5)

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

$ P\left( {\mathit{\boldsymbol{m}},\mathit{\boldsymbol{d}}} \right) = {\left\| {\ln \mathit{\boldsymbol{A}}\left( \mathit{\boldsymbol{m}} \right) - \ln \mathit{\boldsymbol{d}}} \right\|^2} + \alpha \sum\limits_{j = 0}^{2N - 1} {{{\left( {\frac{{\Delta {m_j}}}{{{m_j}}}} \right)}^2}} $

(6)

式中:‖ln A(m)-lnd‖²为数据拟合泛函或拟合差, 后文用Φ₁表示; $\alpha \sum\limits_{j = 0}^{2N - 1} {{{\left( {\frac{{\Delta {m_j}}}{{{m_j}}}} \right)}^2}} $为稳定泛函, 用Φ₂表示。对(6)式求导, 取其导数为零, 得最终的迭代方程:

$ \left( {\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}^{\rm{T}}{\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}} + \alpha \mathit{\boldsymbol{D}}} \right)\Delta \mathit{\boldsymbol{m}} = \mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}^{\rm{T}}\left[ {\ln \mathit{\boldsymbol{d}} - \ln \mathit{\boldsymbol{A}}\left( \mathit{\boldsymbol{m}} \right)} \right] $

(7)

式中:${\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}} = {\mathit{\boldsymbol{U}}_{ij}} = \frac{1}{{{A_i}}} \cdot \frac{{\partial {A_i}}}{{\partial {m_j}}},\mathit{\boldsymbol{D}} = {\rm{diag}}(1/m_j^2)$, i为频点个数循环量, j为模型参数循环量; m在这里为当前迭代的模型参数。

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

$ \begin{array}{l} \mathit{\boldsymbol{m}} = \left[ {\ln {\rho _1}\;\;\; \cdots \;\;\;\ln {\rho _j}\;\;\; \cdots \;\;\;\ln {\rho _N},\ln {h_1}\;\;\; \cdots } \right.\\ \;\;\;\;\;\;\left. {\ln {h_j}\;\;\; \cdots \;\;\;\ln {h_{N - 1}}} \right] \end{array} $

(8)

那么雅可比矩阵中:

$ \frac{{\partial {A_i}}}{{\partial \left( {\ln {\rho _j}} \right)}} = {\rho _j}\frac{{\partial {A_i}}}{{\partial {\rho _j}}} $

(9)

$ \frac{{\partial {A_i}}}{{\partial \left( {\ln {h_j}} \right)}} = {\rho _j}\frac{{\partial {A_i}}}{{\partial {h_j}}} $

(10)

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

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

2) 求变化后的雅可比矩阵${\mathit{\boldsymbol{U}}_{\mathit{\boldsymbol{m}}ij}} = \frac{1}{{{A_i}}} \cdot \frac{{\partial {A_i}}}{{\partial {m_j}}}$和矩阵D=diag(1/m_j²)。

3) 通过解方程$(\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}^T{\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}} + \alpha \mathit{\boldsymbol{D}})\Delta \mathit{\boldsymbol{m}} = \mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}^T[{\rm{ln}}\mathit{\boldsymbol{d}} - {\rm{ln}}\mathit{\boldsymbol{A}}(\mathit{\boldsymbol{m}})]$, 求出Δm。

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

5) 判断条件‖lnA(m+Δm)-lnd‖＜E是否满足, 如果是, 反演终止, 否则返回步骤2)。

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

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

7) 式中U_m^TU_m为非奇异实对称方阵, 作奇异值分解, 存在正交矩阵V和对角矩阵Q, 使得:

$ \mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}^{\rm{T}}{\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}} = \mathit{\boldsymbol{VQ}}{\mathit{\boldsymbol{V}}^{\rm{T}}} $

(11)

定义矩阵H=U_m^TU_m, 同样存在正交矩阵R和对角矩阵Λ, 使得:

$ \begin{array}{l} {\left( {\mathit{\boldsymbol{H}} + \alpha \mathit{\boldsymbol{D}}} \right)^{ - 1}} = {\left( {\mathit{\boldsymbol{VQ}}{\mathit{\boldsymbol{V}}^{\rm{T}}} + \alpha \mathit{\boldsymbol{D}}} \right)^{ - 1}} = {\left( {\mathit{\boldsymbol{R \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)^{ - 1}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^{ - 1}}{\mathit{\boldsymbol{R}}^{\rm{T}}} \end{array} $

(12)

Λ的具体表达式为:

$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = {\rm{diag}}\left[ {\begin{array}{*{20}{c}} {{\lambda _1}}& \cdots &{{\lambda _N}}& \cdots &{{\lambda _{N + 1}}}& \cdots &{{\lambda _{2N - 1}}} \end{array}} \right] $

(13)

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

$ \begin{array}{l} \Delta {m_j} = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^{ - 1}}{\mathit{\boldsymbol{R}}^{\rm{T}}}\mathit{\boldsymbol{U}}_\mathit{\boldsymbol{m}}^{\rm{T}}\left[ {\left( {\ln \mathit{\boldsymbol{d}} - \ln \mathit{\boldsymbol{A}}\left( \mathit{\boldsymbol{m}} \right)} \right.} \right]\\ \;\;\;\;\;\;\; = \sum\limits_{i = 1}^{2N - 1} {{R_{ji}}\frac{1}{{{\lambda _i}}}{K_i}} \end{array} $

(14)

式中:R为U_m^TU_m+αD作奇异值分解后的正交方阵, 阶数等于模型参数的个数。其中:

$ {K_i} = \sum\limits_{k = 1}^{2N - 1} {{R_{ik}}\left[ {\sum\limits_{l = 1}^M {{U_{lk}}\left( {\ln {d_l} - \ln {A_l}} \right)} } \right]} $

(15)

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

2.2 反演参数的改进

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

αD的特征值为α/m_j², 对于同一层的电阻率和层厚, 数值上当ρ_j≫h_j时, αD的层厚对应特征值λ_h(αD)远大于电阻率对应的特征值λ_ρ(αD), 且远大于H的特征值。若将所有的特征值按降序排列, 那么根据附录A中的推论(定理证明见文献[11]), 有:λ_h(H+αD)≫λ_ρ(H+αD), 记为λ_h′≫λ_ρ′, 此时Δρ≫Δh; 当α减小时, α/ρ²和α/h²的数量级等比例减小, 则λ_h′和λ_ρ′的差也逐渐变小, Δρ与Δh相对接近。

综上所述, λ_ρ′与λ_h′数值上逐渐接近使得Δρ与Δh也相对接近, 此时, 相对于电阻率, 数量级较小的层厚更易发生改变。避免取对数的原因在于: Δln ρ≫Δlnh的条件不成立。

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

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

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

两种方法的不同之处在于判断α减少或增大的条件:传统方法是根据Φ₁^(k+1)和Φ₁^(k)的大小来确定, 而本文方法采用判断反演的参数是否为正来确定。由于传统方法每次反演迭代得到的解需要取指数, 指数值则始终大于0, 所以不用考虑参数变为负的问题。为了说明本文方法能够保证方程的稳定性和收敛性以及阻尼因子的快速降低, 我们设计一个简单3层模型:电阻率为100Ω·m的低阻层埋深1.0km, 厚度为0.5km, 背景电阻率为1000Ω·m。正演数据在0.001~6000.000Hz采用30个计算频点, 无随机噪声。反演深度3km, 建立15层模型, 迭代20次, 阻尼因子初始值设置10, 采用不同的阻尼因子变化方案和模型参数得到拟合差结果如表 1所示。

由表 1可知:①本文方法的α在第7次迭代之后就远远低于传统方法的α, 说明新的收敛条件可以保证α更快的下降速度; ②两种方法的拟合差均逐步减少, 本文方法的拟合差没有传统方法的小, 但阻尼因子更低, 导致特征值更低, 观测数据的微小变化将导致模型参数的更大变化; ③本文方法中第2次迭代到第3次迭代拟合差不降反升, 且后期震荡变化, 整体为降低的趋势, 说明阻尼因子起着调节特征值和稳定性的作用, 当层厚修正量过大使层厚出现负值时, 则提高阻尼因子, 降低层厚修正量, 保证方程稳定性; 若层厚为正值, 则降低阻尼因子, 保证分辨率。以上结论说明通过判断模型参数是否为负值来保证反演的稳定性和收敛性是合理的。

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可以看出, 电性过度界面节点较为密集, 说明该处的层厚受到压缩, 在一定程度上提高了分辨率, 同时表明, 利用层厚的小特征值可以提高模型空间的分辨率, 而传统方法拟合差更低。该结论与小特征值对应高分辨率, 大特征值意味着以牺牲一些分辨率为代价来换取较小的方差的结论一致。

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

2.5 初始模型和解方程组

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

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

3 合成数据测试及分析

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

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

4 结束语

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

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

附录A

以下所有矩阵的特征值满足降序排列λ₁＞λ₂＞…＞λ_n。

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

${\rm{ma}}{{\rm{x}}_{r + s = n + i}}\left[ {{\lambda _r}\left( \mathit{\boldsymbol{A}} \right) + {\lambda _s}\left( \mathit{\boldsymbol{B}} \right)} \right] ≪ {\lambda _i}\left( {\mathit{\boldsymbol{A}} + \mathit{\boldsymbol{B}}} \right)≪{\rm{mi}}{{\rm{n}}_{r + s = i + 1}}[{\lambda _r}\left( \mathit{\boldsymbol{A}} \right) + {\lambda _s}\left( \mathit{\boldsymbol{B}} \right)]$

(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≤s₁≤N; λ_i(A+B)的后N-1个特征值的上、下界均包含λ_s2(B), 不含λ_s1(B), 且N+1≤s₂≤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个特征值的差逐渐变小。