病态问题在控制网平差、卫星快速定位、航空重力向下延拓等领域的数据处理时经常碰到。当系数矩阵列向量为复共线性相关时，最小二乘估计(least squares, LS)虽然为无偏估计，但参数估值已严重失真，不是最优解^[1]。解决病态问题的常用方法有截断奇异值法^[2](truncated singular value decomposition, TSVD)、Tikhonov正则化法^[3-4]、岭估计法^[5]等。截断奇异值法是基于奇异值分解技术的病态方程的一种直接解法，其原理是将较小的奇异值删除，保留较大的奇异值进行解算^[6]。Tikhonov正则化法是在LS估计准则上加入稳定泛函约束条件，目的是将模型由严重病态变成病态性较弱或无病态，牺牲LS估计的无偏性，引入一定偏差使得方差降低，参数估值稳定可靠。正则化法和岭估计法都需要引入参数来调节法矩阵和正则化矩阵的相关关系，参数可由L-曲线^[7-8]、岭迹法^[9]、广义交叉核实法(generalized cross validation, GCV)^[10-11]等方法确定。岭估计法是正则化法的一种简化特例，即将正则化矩阵换成单位阵，通过合理地选择正则化参数和正则化矩阵，可以提高参数估计的可靠性和稳定性。

Tikhonov正则化法和岭估计法都是有偏估计，在降低方差的同时引入偏差。林东方等^[6]分析了岭估计的偏差和方差，提出一种确定正则化矩阵的方法，在降低方差的同时避免不必要的偏差，使得估值更为合理。shen等^[12]、徐天问等^[13]提出以均方误差作为衡量参数估值可靠性指标，定义均方误差是方差迹与偏差迹之和。林东方等^[14]提出引入偏差量小于降低方差量的思想，在降低方差的同时控制偏差的引入，更为合理地解决病态问题。根据均方误差最小原则，为合理分配方差与偏差，本文提出岭估计迭代法。首先对病态数据用岭估计法求解，然后将岭估计解代入观测模型，更新观测向量，再次用岭估计法计算参数估值，以此形成迭代。通过分析迭代估计结果的方差、估计偏差及其变化趋势将均方误差最小或收敛作为迭代终止条件，并用模拟算例验证该方法的有效性与可行性，最后分析估值的估计偏差和实际偏差的关系。

1 岭估计法

对于经典平差模型:

$ \mathit{\boldsymbol{L}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} = \mathit{\boldsymbol{AX}} $

(1)

式中，L为观测值，A为系数矩阵，X为待定参数，Δ为观测向量的噪声，Δ~N(0, σ₀² I)，σ₀²是单位权方差。其最小二乘估计及其协方差为：

$ \mathit{\boldsymbol{\hat X}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} $

(2)

$ {\mathop{\rm cov}} \left( {\mathit{\boldsymbol{\hat X}}} \right) = \sigma _0^2{\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}} \right)^{ - 1}} $

(3)

最小二乘估计是无偏估计，方差可由协方差矩阵之迹表示:

$ D\left( {\mathit{\boldsymbol{\hat X}}} \right) = {\rm{tr}}\left[ {{\mathop{\rm cov}} \left( {\mathit{\boldsymbol{\hat X}}} \right)} \right] = \sigma _0^2\sum\limits_{i = 1}^n {\frac{1}{{\lambda _i^2}}} $

(4)

式中, λ_i是系数矩阵A的奇异值。当法矩阵A^TA出现病态(条件数一般大于10³，条件数等于λ_max²/λ_min²)，即存在λ_i接近0，从式(4)可以看出，方差将变得异常大，并且法矩阵求逆将会表现得很不稳定，导致求解出的参数估值不可靠^[11]。

在解决系数矩阵病态问题时，Tikhonov正则化方法加入了稳定泛函Ω(X)，相应的准则为:

$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}}} \right\|_2^2 + \alpha \mathit{\Omega }\left( \mathit{\boldsymbol{X}} \right) = }\\ {\left\| {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}}} \right\|_2^2 + \alpha {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{RX}} = \min } \end{array} $

(5)

式中，R为正则化矩阵，‖·‖为欧氏2-范数，α为大于0的正则化参数。当R = I则为岭估计，岭估计是正则化方法的一种特殊简化形式，其平差结果为:

$ {{\mathit{\boldsymbol{\hat X}}}_0} = {\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} $

(6)

其中，N = A^TA，则岭估计的方差、偏差分别为:

$ \begin{array}{*{20}{c}} {D\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) = {\rm{tr}}\left[ {{\mathop{\rm cov}} \left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right)} \right] = {\rm{tr}}\left[ {{{\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}} \cdot } \right.}\\ {\left. {\mathit{\boldsymbol{N}}{{\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right] = \sigma _0^2\sum\limits_{i = 1}^n {\frac{{\lambda _i^2}}{{{{\left( {\lambda _i^2 + \alpha } \right)}^2}}}} } \end{array} $

(7)

$ {\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) = \mathit{\boldsymbol{\tilde X}} - E\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) = \alpha {\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)^{ - 1}}\mathit{\boldsymbol{\tilde X}} $

(8)

式中, $\mathit{\boldsymbol{\tilde X}}$是参数真值，根据均方误差定义^[13]:

$ {\rm{MSE}}\left( {\mathit{\boldsymbol{\hat X}}} \right) = E{\left( {\mathit{\boldsymbol{\tilde X}} - \mathit{\boldsymbol{\hat X}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{\tilde X}} - \mathit{\boldsymbol{\hat X}}} \right) $

(9)

则岭估计解均方误差为:

$ \begin{array}{*{20}{c}} {{\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) = D\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) + {\rm{tr}}\left( {{\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right){\rm{bias}}{{\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right)}^{\rm{T}}}} \right) = }\\ {\sigma _0^2{\rm{tr}}\left[ {{{\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}}\mathit{\boldsymbol{N}}{{\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right] + }\\ {\sigma _0^2{\rm{tr}}\left[ {{{\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}}\mathit{\boldsymbol{\tilde X}}{{\mathit{\boldsymbol{\tilde X}}}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right]} \end{array} $

(10)

对比分析式(4)、式(7)、式(10)可知，岭估计有效地降低了方差，却引入了偏差。因此，在处理病态数据时以均方误差作为衡量参数估值可靠性指标，降低方差的同时需要合理地控制偏差的引入量，让方差和偏差合理分配，使得参数估值稳定可靠。

2 岭估计迭代法

在有偏估计中，方差与偏差并存。林东方等^[14]提出引入偏差量小于降低方差量的思想，在降低方差的同时控制偏差的引入，更为合理地解决病态问题。参数估值的均方误差是方差之迹与偏差之迹的和，当参数估值的方差降低量大于偏差引入量时视为有效解。要得到最优有效解，则需均方误差最小。根据均方误差最小原则，合理分配方差与偏差，本文提出岭估计迭代解法。

观测模型如式(1)，则岭估计解为:

$ {{\mathit{\boldsymbol{\hat X}}}_0} = {\left( {\mathit{\boldsymbol{N}} + \alpha \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} $

(11)

将${\mathit{\boldsymbol{\hat X}}}$₀代入观测模型(1)中，更新观测向量L，得到${{\mathit{\boldsymbol{\hat L}}}_0} = \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_0}$，再次用岭估计法求解参数估值，得到迭代式:

$ {{\mathit{\boldsymbol{\hat X}}}_{k + 1}} = {\left( {\mathit{\boldsymbol{N}} + {\alpha _{k + 1}}\mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_k} $

(12)

迭代初值${{\mathit{\boldsymbol{\hat X}}}_0}$为岭估计解。

2.1 岭估计迭代法的收敛性及终止条件

若迭代矩阵的谱半径小于1，则迭代可以收敛^[15-16]。本文岭估计迭代法的迭代矩阵为(N+α_iI)^-1 A^TA，在对其特征值分解得到其特征值为Λ_f/(Λ_f+α_i)。谱半径是方阵的最大特征值，而迭代矩阵的特征值均小于1，因此迭代是可以收敛的。

与一般迭代公式终止条件不同，本文迭代式终止条件并不是‖${\mathit{\boldsymbol{\hat X}}}$_k-${\mathit{\boldsymbol{\hat X}}}$_k－1‖₂²。因为要在均方误差意义下求得最优解，本文方法的主要目的是求得迭代过程中的最小均方误差及此时的参数估值，而在迭代过程中，参数估值的方差不断减小而偏差却不断增加，因此为了求得最优有效解，该迭代式终止条件是均方误差最小或收敛。

3 岭估计迭代法方差、偏差、均方误差分析

岭估计迭代法相比于岭估计在迭代过程中逐渐降低方差，然而也逐渐提高了偏差。

由式(12)可以得到:

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_k} = {{\left( {\mathit{\boldsymbol{N}} + {\alpha _k}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{{\left( {\mathit{\boldsymbol{N}} + {\alpha _{k - 1}}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} \cdots }\\ {\left( {\mathit{\boldsymbol{N}} + {\alpha _0}\mathit{\boldsymbol{I}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}}} \end{array} $

(13)

根据协方差传播定律^[16]可得协方差矩阵:

$ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\mathop{\rm cov}} \left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) = {{\left( {\mathit{\boldsymbol{N}} + {\alpha _k}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{{\left( {\mathit{\boldsymbol{N}} + {\alpha _{k - 1}}\mathit{\boldsymbol{I}}} \right)}^{ - 1}} \cdot }\\ {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} \cdots {{\left( {\mathit{\boldsymbol{N}} + {\alpha _0}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{{\left( {\mathit{\boldsymbol{N}} + {\alpha _0}\mathit{\boldsymbol{I}}} \right)}^{ - 1}} \cdots } \end{array}}\\ {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{{\left( {\mathit{\boldsymbol{N}} + {\alpha _{k - 1}}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{{\left( {\mathit{\boldsymbol{N}} + {\alpha _k}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \end{array} $

(14)

方差可由协方差矩阵之迹表示为:

$ \begin{array}{*{20}{c}} {D\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) = {\rm{tr}}\left[ {{\mathop{\rm cov}} \left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right)} \right] = }\\ {\sigma _0^2\sum\limits_{i = 1}^n {\frac{{\mathit{\Lambda }_i^{2k + 1}}}{{{{\left[ {\left( {{\mathit{\Lambda }_i} + {\alpha _n}} \right)\left( {{\mathit{\Lambda }_i} + {\alpha _{n - 1}}} \right) \cdots \left( {{\mathit{\Lambda }_i} + {\alpha _0}} \right)} \right]}^2}}}} } \end{array} $

(15)

式中，Λ_i是法矩阵A^TA的特征值。与式(7)相比，可知：

$ D\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) < D\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) $

(16)

随迭代次数增加，岭估计迭代解将逐渐降低方差，而其偏差为:

$ \begin{array}{*{20}{c}} {{\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) = \mathit{\boldsymbol{\tilde X}} - E\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) = \mathit{\boldsymbol{\tilde X}} - {{\left( {\mathit{\boldsymbol{N}} + {\alpha _k}\mathit{\boldsymbol{I}}} \right)}^{ - 1}} \cdot }\\ {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{N}} + {\alpha _{k - 1}}\mathit{\boldsymbol{I}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} \cdots \left( {\mathit{\boldsymbol{N}} + {\alpha _0}\mathit{\boldsymbol{I}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A\tilde X}} = }\\ {\left[ {\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{G}}_i}\left( {1 - \frac{{\mathit{\Lambda }_i^{n + 1}}}{{\left( {{\mathit{\Lambda }_i} + {\alpha _n}} \right)\left( {{\mathit{\Lambda }_i} + {\alpha _{n - 1}}} \right) \cdots \left( {{\mathit{\Lambda }_i} + {\alpha _0}} \right)}}} \right)\mathit{\boldsymbol{G}}_i^{\rm{T}}} } \right]\mathit{\boldsymbol{\tilde X}}} \end{array} $

(17)

式中，G_i为法矩阵A^TA的特征向量，${\mathit{\boldsymbol{\tilde X}}}$为参数真值。值得注意的是，由于参数真值${\mathit{\boldsymbol{\tilde X}}}$通常是不可知的，参数估计偏差无法准确算出。Xu^[2]提出，可以采用参数的最小二乘估值或其他估值替代真值，而林方东等^[14]提到，${\mathit{\boldsymbol{\tilde X}}}$真值由某估值替代时, 均方误差实际上反映的是参数估值与真值替代值之间的离散程度, 因此真值替代值须具备较高的准确性、稳定性才能得到可靠的均方误差值。而最小二乘估值由于法矩阵病态导致严重偏离真值，故不适合采用。所以，在选择真值替代值时，可以选择岭估计估值代替真值，岭估计估值相比最小二乘估值更具准确性、稳定性。与式(8)相比可知，在迭代过程中，岭估计迭代解进一步降低方差的同时也再次引入了偏差，偏差将逐渐增大。

$ {\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) > {\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) $

(18)

因此，可以得到${\mathit{\boldsymbol{\hat X}}}$_k的均方误差为:

$ \begin{array}{*{20}{c}} {{\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) = D\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) + {\rm{tr}}\left[ {{\rm{bias}}{{\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right)}^{\rm{T}}}{\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right)} \right] = }\\ {\sigma _0^2\sum\limits_{i = 1}^n {\frac{{\mathit{\Lambda }_i^{2k + 1}}}{{{{\left[ {\left( {{\mathit{\Lambda }_i} + {\alpha _n}} \right)\left( {{\mathit{\Lambda }_i} + {\alpha _{n - 1}}} \right) \cdots \left( {{\mathit{\Lambda }_i} + {\alpha _0}} \right)} \right]}^2}}}} + }\\ {\sum\limits_{i = 1}^n {{{\mathit{\boldsymbol{\tilde X}}}^{\rm{T}}}{\mathit{\boldsymbol{G}}_i}} }\\ {\left( {1 - \frac{{\mathit{\Lambda }_i^{2k + 2}}}{{{{\left[ {\left( {{\mathit{\Lambda }_i} + {\alpha _n}} \right)\left( {{\mathit{\Lambda }_i} + {\alpha _{n - 1}}} \right) \cdots \left( {{\mathit{\Lambda }_i} + {\alpha _0}} \right)} \right]}^2}}}} \right)\mathit{\boldsymbol{G}}_i^{\rm{T}}\mathit{\boldsymbol{\tilde X}}} \end{array} $

(19)

由于各种参数均未知，本文方法与岭估计的均方误差无从比较。因此，为了体现本文方法的优势，将计算过程中均方误差达到最小值或收敛于某一值作为迭代终止条件：

$ {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_k}} \right) = D\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) - \sum\limits_{j = 1}^k {{d_j}} + B\left( {{{\mathit{\boldsymbol{\hat X}}}_0}} \right) + \sum\limits_{j = 1}^k {{b_j}} $

(20)

式中, D(${\mathit{\boldsymbol{\hat X}}}$₀)、B(${\mathit{\boldsymbol{\hat X}}}$₀)分别为岭估计方差和偏差，d_j、b_j(大于0)分别为本文方法迭代第j次过程中再次降低的方差量和再次引入的偏差量，可以通过一次差分求得。若d_j>b_j，则均方误差在逐渐减小，直到第k次迭代时，d_k=b_k, d_k+1 < b_k+1或d_k>b_k, d_k+1 < b_k+1，则此时均方误差最小，之后迭代d_j将小于b_j，均方误差将增大；或在多次迭代中，d_j总是大于b_j，则均方误差一直降低直至收敛。再以均方误差作为衡量参数估值标准，${\mathit{\boldsymbol{\hat X}}}$_k就是最优解。

下面给出本文方法的解算步骤：

1) 选择L-曲线法或其他方法确定岭参数，计算出岭估计解${\mathit{\boldsymbol{\hat X}}}$₀；

2) 将计算出的${\mathit{\boldsymbol{\hat X}}}$₀代入观测模型，更新观测向量L，${\mathit{\boldsymbol{\hat L}}}$₀= A${\mathit{\boldsymbol{\hat X}}}$₀；

3) 再次用L-曲线法或其他方法选择岭参数，再计算出岭估计解${\mathit{\boldsymbol{\hat X}}}$₁，并更新观测向量；

4) 重复步骤1)、2)、3)；

5) 因为要求得最小均方误差，而不知道第几次迭代得到最小值，所以计算时先任取一迭代次数，在迭代过程中计算出参数估值和求出其均方误差，并观察均方误差变化趋势，当均方误差最小或收敛时，此时的估值就为最优有效解。

4 算例及分析 4.1 算例1

采用文献[17]中的模拟病态问题算例，法矩阵条件数为2.083 8×10⁴，严重病态。未知参数有5个，X=[x₁  x₂  x₃  x₄  x₅]^T，真值为X=[1  1  1  1  1]^T，观测值L的噪声Δ~N(0, σ₀²I), σ₀=1(具体数据未列出，详情可见参考文献)。分别用最小二乘、岭估计、截断奇异值法和本文方法进行平差解算，本文方法按照上述方法的解算步骤计算，先确定迭代次数(10、100、1 000、3 000、5 000)，再根据均方误差变化趋势确定合理迭代次数。各方法结果、估计均方误差MSE(${\mathit{\boldsymbol{\hat X}}}$)、实际均方误差RMSE(${\mathit{\boldsymbol{\hat X}}}$)及差值范数‖Δ${\mathit{\boldsymbol{\hat X}}}$(‖Δ${\mathit{\boldsymbol{\hat X}}}$‖=‖ ${\mathit{\boldsymbol{\tilde X}}}$-${\mathit{\boldsymbol{\hat X}}}$‖)等见图 1~3及表 1。

4.2 算例2

采用文献[18]模拟空间测边网算例。P₁、P₂、…、P₁₀为10个已知点，其坐标具体数据略去，10个已知点到3个未知点P₁₁、P₁₂、P₁₃(假设模拟坐标真值分别为(0 m，0 m，0 m)、(68 m，-26 m，9 m)和(14 m，41 m，-11 m))的距离以及3个未知点间的距离假定已通过测量得到。设各距离为等精度观测，中误差为±0.01 m。根据33个距离观测值确定3个未知点坐标。计算中3个未知点坐标近似值分别取(0.03 m，-0.025 m，0.01 m)、(68.03 m，-25.97 m，8.98 m)和(14.04 m，40.97 m，-11.04 m)。该测边网所建立观测方程的系数阵A严重病态，法矩阵条件数为89 543。根据计算步骤，先确定迭代次数(10、100、1 000、3 000、5 000)，再根据均方误差变化趋势确定合理迭代次数。各种算法结果、估计均方误差、实际均方误差及差值范数等见图 4~6及表 2。

4.3 算例分析

由于本文算法是先任取迭代次数，经多次反复计算，算例1虽然在前几次迭代的结果已符合方法要求，但为了具体分析问题，最终选择30次迭代，记录结果并作出方差、偏差、均方误差及差分图表。算例2在迭代1 000次左右变化缓慢，从1 000次到3 000次迭代，结果变化非常小，从3 000次到5 000次迭代，结果基本不变，因此最终选择5 000次迭代。

从图 1(a)、图 4(a)可以看出，方差在降低，估计偏差与实际偏差在增加，变化情况与分析一致。算例1在前5次迭代计算中，方差下降较快，估计偏差与实际偏差增长也快；之后方差降低非常缓慢，而估计偏差与实际偏差依然增加。方差的降低量、估计偏差与实际偏差的增加量均大于方差的降低量，再继续计算均方误差将增大。算例2在前5次迭代计算中，方差下降较快，偏差增加较快；迭代几百次之后，方差降低较为缓慢，而偏差增加也非常缓慢，再之后的迭代中，方差与偏差变化均非常小，迭代结果仍有细小变化，一定的迭代次数之后就有可见变化。从图 4(b)可以看出，估计均方误差与实际均方误差均在减小且前几次迭代降低较快之后缓慢直至收敛。

图 2、图 3、图 5、图 6是方差与偏差的差分图，由差分图可进一步得出均方误差的变化情况。由图可见，方差、偏差及均方误差前几次迭代结果变化较快，之后则是均匀变化。在以均方误差作为衡量参数估值的标准时，算例1第5次迭代值为最优解；算例2迭代到5 000次时，结果基本不变，所以取第5 000次迭代值为最优解。从表 1可以看出，岭估计和本文方法的估计均方误差分别是0.991 7、0.048 0，实际均方误差分别为1.510 8、0.704 4；表 2的岭估计和本文方法的估计均方误差分别是22.704 7、0.496 5，实际均方误差分别为22.684 5、0.001 4，可见本文方法的均方误差比岭估计低；相比岭估计，算例1、算例2中本文方法差值范数也有所降低，最终结果比TSVD更优且均方误差更小。

本文在计算估计偏差时采用岭估计的估值代替真值，而从图 1(a)、图 2(b)、图 4(a)、图 5(b)可以看出，虽然估计偏差与实际偏差有一定的差值，但除前几次迭代偏差变化较大之外，之后迭代中两者的变化情况高度一致，规律相似。可以得出，估计均方误差与实际均方误差虽然有差值，但两者增加或减少的规律是一致的，在均方误差意义下，当估计均方误差达到最小时，其实际均方误差也将达到最小。因此，可以采用岭估计的估值代替真值求解均方误差。

在病态数据处理中，相较于岭估计，本文岭估计迭代法能有效地提高参数估值的稳定性、可靠性，进一步降低均方误差和差值范数，更具优势。

5 结语

本文岭估计迭代法在迭代计算中可使方差逐渐降低，偏差逐渐增加。根据林方东等^[14]提出的方差降低量大于偏差引入量的思想，若d_j>b_j，则结果有效，直到d_j≤b_j迭代终止，此时均方误差最小。通过合理次数的迭代计算得到最小或收敛的均方误差，取此时的参数估值为最优解，并用算例验证该方法，计算得到的均方误差、差值范数以及参数估值均优于LS估计、岭估计和TSVD，说明该方法有效、可行。

岭估计迭代法需要先选择迭代次数，可通过选择不同迭代次数，作出均方误差变化趋势图，一步步确定合理的迭代次数，根据图形选出最小或收敛的均方误差。迭代过程虽较为复杂，但可用MATLAB等计算软件编程实现，自动找出最小均方误差和相应的参数估值，故本文方法对计算效率无太大影响。在算例2中，本文为更好地说明问题而选择的迭代次数过多，在解决实际问题时，无需像算例2一样过多迭代，当迭代收敛时，设定一阈值终止迭代即可。b