解决病态问题的常用方法有截断奇异值法^[1-2](truncated singular value decomposition, TSVD)、Tikhonov正则化法^[3-4]、岭估计法^[5-6]等。截断奇异值法是基于奇异值分解技术的直接解算方法, 仅利用较大奇异值包含的观测信息进行参数估计，将严重扩大方差的小奇异值截去^[7]。Tikhonov正则化法是在LS估计准则上加入稳定泛函，将模型由严重病态变成病态性较弱或无病态。正则化法需要合理选择并引入正则化参数和正则化矩阵得到稳定的参数估值。在正则化参数的选择上，可使用L-曲线法^[8-9]、岭迹法^[10]及广义交叉核实法(generalized cross validation, GCV)^[11-12]等。而对于正则化矩阵的选取，目前研究较少。

通常将单位阵作为正则化矩阵，称作岭估计。这种方法通过修正所有奇异值来达到改善法矩阵病态性的目的，得到稳定解。岭估计是有偏估计^[13]。林方东等^[14]对岭估计的方差和偏差进行分析，提出用法矩阵小奇异值对应的特征向量构造正则化矩阵，在降低方差的同时避免不必要的偏差，使估值更为合理。王振杰等^[15]提出将参数的后验方差的逆矩阵主对角线元素作为正则化矩阵，得到参数估值的均方误差小于最小二乘(least squares, LS)估计，这个正则化矩阵可以看作是待定参数的权重，可是待定参数不是观测值，也没有实际意义上的权重^[14]。

在法矩阵病态的情况下，参数估计的稳定性常采用均方误差作为衡量指标，是方差与偏差之和^[15-16]。在均方误差意义下，针对正则化矩阵的选取问题，本文首先根据法矩阵奇异值分解得到的小奇异值所对应的特征向量构造一个对称矩阵，再选取该矩阵的主对角线元素构造一个对角矩阵，将这个对角矩阵与单位阵组合出一种新的正则化矩阵，最后利用模拟数据进行验证。

1 岭估计原理

对于经典平差模型：

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

(1)

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

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

(2)

$ D(\mathit{\boldsymbol{\hat X}}) = \sigma _0^2{\rm{tr}}{({\mathit{\boldsymbol{A}}^{\rm T}}\mathit{\boldsymbol{A}})^{ - 1}} = \sigma _0^2\sum\limits_{i = 1}^n {\frac{1}{{\lambda _i^2}}} $

(3)

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

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

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

(4)

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

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

(5)

将式(5)变形得到：

$ {{\mathit{\boldsymbol{\hat X}}}_1} = {[\mathit{\boldsymbol{G}}(\mathit{\boldsymbol{ \boldsymbol{\varLambda} }} + \alpha \mathit{\boldsymbol{I}}){\mathit{\boldsymbol{G}}^{\rm T}}]^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm T}}\mathit{\boldsymbol{L}} $

(6)

式中，G为法矩阵 A^TA 所有特征值所对应的特征矩阵，Λ是法矩阵 A^TA 的特征值，岭估计则相当于在法矩阵每个特征值上加上α，修正所有奇异值，以此改善法矩阵的病态性。

2 正则化矩阵构造新方法

在文献[15]中，正则化矩阵构造是由岭估计得到的均方差矩阵MSEM(${{\mathit{\boldsymbol{\hat X}}}_1}$)求逆，然后取其主对角线元素，即R =diag(MSEM(${{\mathit{\boldsymbol{\hat X}}}_1}$)^－1)。

在文献[14]中，正则化矩阵是基于小奇异值特征向量构造而成：

$ {\mathit{\boldsymbol{R}}_1} = \sum\limits_{i = k}^n {{\mathit{\boldsymbol{G}}_i}\mathit{\boldsymbol{G}}_i^{\rm{T}}} $

(7)

式中，G_i是小奇异值对应的特征向量。

文献[15]得到的矩阵R是对角矩阵，且矩阵的对角元素非常大，容易使偏差异常大。因此，要得到可靠的参数估计，就要求正则化参数非常准确，而正则化参数往往难以准确确定^[13]。在文献[14]中，矩阵R₁只修正小奇异值，参数则与岭参数相同，参数估计无明显优势，其偏差相应降低。因此，综合上述2种方法，本文提出构造正则化矩阵的新方法，不需要再求正则化参数，直接代入公式计算即可。

本文构造正则化矩阵的方法是由式(7)生成的矩阵R₁主对角线元素乘以正则化参数α，得到矩阵M，则矩阵M可以变换成：

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{R}}_1} + \left( {\alpha - 1} \right){\mathit{\boldsymbol{R}}_2} $

(8)

式中，R₂为R₁主对角线元素构成的对角矩阵, 将式(8)中的R₁替换成I，则构成本文所提正则化矩阵R₃：

$ {\mathit{\boldsymbol{R}}_3} = \mathit{\boldsymbol{I}} + \left( {\alpha - 1} \right){\mathit{\boldsymbol{R}}_2} $

(9)

平差结果为：

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

(10)

3 新正则化矩阵适用性分析

加入新正则化矩阵的参数估计为：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_2} = {{\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{R}}_3}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} = }\\ {{{\left[ {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} + \alpha \mathit{\boldsymbol{I}} + \left( {1 - \alpha } \right)\left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{R}}_2}} \right)} \right]}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}}} \end{array} $

(11)

令A^TA +α I = A₁^TA₁，R₄=(1－α)(I － R₂)，根据矩阵反演公式^[7]可得:

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

(12)

整理得到：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_2} = {{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} - }\\ {\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}{{\left( {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + \mathit{\boldsymbol{R}}_4^{ - 1}} \right)}^{ - 1}}} \right]{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}}} \end{array} $

(13)

从式(13)可看出，${{\mathit{\boldsymbol{\hat X}}}_2}$是${{\mathit{\boldsymbol{\hat X}}}_1}$的一个线性变换，R₂是主对角元素均小于1大于0的正定矩阵：

$ {\mathit{\boldsymbol{R}}_2} = \left[ {\begin{array}{*{20}{c}} {{H_1}}&{}&{}&{}\\ {}&{{H_2}}&{}&{}\\ {}&{}& \cdots &{}\\ {}&{}&{}&{{H_n}} \end{array}} \right] $

(14)

其中，${H_f} = \sum\limits_{j = k}^n {g_{fj}^2} $，而$\sum\limits_{j = k}^n {g_{fj}^2} = 1$，g_fj是特征向量组成的特征矩阵G上的元素，所以H_f < 1，也得到(I － R₂)为对角元素小于1的正定矩阵。

在文献[15, 18]中，均方误差采用了简化形式，虽然这不是严格意义上的均方误差，但该简化形式便于比较均方误差大小而得到采用，而在算例需要计算均方误差时，则采用其实际均方误差计算公式^[16]。因此，采用文献[15, 18]的均方误差简化形式来比较均方误差大小，可得${{\mathit{\boldsymbol{\hat X}}}_2}$的均方误差阵为：

$ \begin{array}{*{20}{c}} {{\rm{MSEM}}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right) = \hat \sigma _1^2{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1} + {\mathit{\boldsymbol{R}}_4}} \right)}^{ - 1}} = }\\ {\hat \sigma _1^2\left\{ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} - {{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + } \right.} \right.}\\ {\left. {{{\left. {\mathit{\boldsymbol{R}}_4^{ - 1}} \right]}^{ - 1}}{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}} \right\} = }\\ {{\rm{MSEM}}\left( {{{\mathit{\boldsymbol{\hat X}}}_1}} \right) - \hat \sigma _1^2\left\{ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + } \right.} \right.}\\ {\left. {{{\left. {\mathit{\boldsymbol{R}}_4^{ - 1}} \right]}^{ - 1}}{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}} \right\}} \end{array} $

(15)

所以，${{\mathit{\boldsymbol{\hat X}}}_2}$的均方差为：

$ \begin{array}{*{20}{c}} {{\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right) = {\rm{tr}}\left( {{\rm{MSEM}}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right)} \right) = }\\ {{\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_1}} \right) - \hat \sigma _1^2{\rm{tr}}\left\{ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + } \right.} \right.}\\ {{{\left. {\mathit{\boldsymbol{R}}_4^{ - 1}} \right]}^{ - 1}}{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}} \end{array} $

(16)

式中，$\hat \sigma _1^2$为单位权方差。

1) 当0 < α < 1，R₄=(1－α)(I － R₂)和A^TA均是正定矩阵，所以，

$ \begin{array}{*{20}{c}} {\hat \sigma _1^2{\rm{tr}}\left\{ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + } \right.} \right.}\\ {\left. {{{\left. {\mathit{\boldsymbol{R}}_4^{ - 1}} \right]}^{ - 1}}{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}} \right\} > 0} \end{array} $

(17)

由式(17)可知：

$ {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right) < {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_1}} \right) $

(18)

所以，在0 < α < 1的条件下，新方法正则化解在均方差意义下优于岭估计。

2) 当α=1，得到R₄=0, 则A₁^TA₁+ R₄= A^TA + I，新正则化矩阵退化成岭估计的单位阵，新方法与岭估计相同，2种方法求得的结果也相同。因此在正则化参数等于1时，岭估计是本文方法的一种特例。

3) 当α>1，R₄=(1－α)(I － R₂)是负定矩阵，证明得到的[(A₁^TA₁)^－1+ R₄^-1]也为负定矩阵。证明如下：

${(\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1})^{ - 1}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{G}}_i}(\frac{1}{{{\mathit{\Lambda }_i} + \alpha }})\mathit{\boldsymbol{G}}_i^{\rm{T}}} $, A^TA是正定矩阵，所以$\frac{1}{{{\mathit{\Lambda }_i} + \alpha }} < 1$。由式(14)知，R₂是对角元素小于1的对角矩阵，任取R₂的一个主对角线元素H_i，0 < H_i < 1，与R₄的主对角线元素L_i对应，根据R₄=(1－α)(I － R₂)得到L_i < 0，再求出：

$ {L_i} = 1 - {H_i} - \alpha + \alpha {H_i} $

(19)

两边加上负号并添加限制条件得到：

$ - {L_i} = - 1 + {H_i} + \alpha - \alpha {H_i} < \alpha < {\mathit{\Lambda }_j} + \alpha $

(20)

两边再求倒数并移项：

$ \frac{1}{{{L_i}}} < - \frac{1}{{{\mathit{\Lambda }_j} + \alpha }} $

(21)

可知，R₄^－1是主对角线元素$\frac{1}{{{L_i}}} < - \frac{1}{{{\mathit{\Lambda }_j} + \alpha }}$的对角矩阵，令$\mathit{\boldsymbol{R}}_4^{ - 1} = - \frac{1}{{{\mathit{\Lambda }_j} + \alpha }}\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{F}}$, 其中，令Λ_j是法矩阵 A^TA 的最小特征值，F是一假设矩阵，得到$\frac{1}{{{\mathit{\Lambda }_j} + \alpha }}$最大。由式(21)可知，F是正定矩阵，可从式(6)得：

$ \begin{array}{*{20}{c}} {\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + R_4^{ - 1}} \right] = \left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} - \frac{1}{{{\mathit{\Lambda }_j} + \alpha }}\mathit{\boldsymbol{I}} - } \right.}\\ {\left. \mathit{\boldsymbol{F}} \right] = \left[ {\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{G}}_i}\left( {\frac{1}{{{\mathit{\Lambda }_i} + \alpha }} - \frac{1}{{{\mathit{\Lambda }_j} + \alpha }}} \right)\mathit{\boldsymbol{G}}_i^{\rm{T}} - \mathit{\boldsymbol{F}}} } \right]} \end{array} $

(22)

从式(22)可以看出，[(A₁^TA₁)^－1+ R₄^-1]是负定矩阵，${{\mathit{\boldsymbol{\hat X}}}_2}$的均方误差为：

$ \begin{array}{*{20}{c}} {{\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right) = {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_1}} \right) - }\\ {\hat \sigma _1^2{\rm{tr}}\left\{ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}{{\left[ {\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)} \right]}^{ - 1}} + {{\left. {\mathit{\boldsymbol{R}}_4^{ - 1}} \right]}^{ - 1}}{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}} \right\}} \end{array} $

(23)

从而得到：

$ \begin{array}{*{20}{c}} {\hat \sigma _1^2{\rm{tr}}\left\{ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}\left[ {{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}} + } \right.} \right.}\\ {\left. {{{\left. {\mathit{\boldsymbol{R}}_4^{ - 1}} \right]}^{ - 1}}{{\left( {\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}} \right)}^{ - 1}}} \right\} < 0} \end{array} $

(24)

由式(24)可得：

$ {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right) > {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_1}} \right) $

(25)

所以，在α>1的条件下，新方法正则化解在均方误差意义下劣于岭估计。

在大地测量数据处理中，病态主要是系数矩阵 A 存在接近于0的奇异值，而导致法矩阵求逆表现不稳定。由于法矩阵病态，最大奇异值与最小奇异值相差几个数量级，由式(6)可以看出，较小的正则化参数就可以在较大程度上修正奇异值，改善矩阵病态性和降低条件数。

综上所述，本文正则化矩阵是根据正则化参数来确定，即先确定参数再确定矩阵，用L-曲线法或其他方法求出参数，再根据式(9)确定新正则化矩阵，最后将此矩阵代入式(10)直接参与计算而不需要再次求取正则化参数，计算出参数估值。同时，比较了新方法与岭估计在不同参数区间下的优劣。下文算例中将求出本文方法下的实际均方误差，并与岭估计均方误差进行比较，以验证本文理论。

4 算例及分析 4.1 算例1

采用文献[19]中的模拟病态问题算例，法矩阵条件数为2.083 8×10⁴，严重病态。未知参数有5个，X =[x₁x₂x₃x₄x₅]^T，真值为X =[1 1 1 1 1]^T，中误差为±0.1(具体数据未列出，详情可见文献[19])。为比较新方法与岭估计解算结果的优劣，分别用最小二乘法、截断奇异值法、岭估计和新方法进行平差解算，并比较参数估值的均方误差MSE(${\mathit{\boldsymbol{\hat X}}}$)及估值与真值的差值范数$\left\| {\Delta \mathit{\boldsymbol{\hat X}}} \right\|$(‖$\left\| {\Delta \mathit{\boldsymbol{\hat X}}} \right\| = \left\| {\mathit{\boldsymbol{\tilde X}} - \mathit{\boldsymbol{\hat X}}} \right\|$，${\mathit{\boldsymbol{\tilde X}}}$是参数真值)，结果见表 1。

4.2 算例2

采用文献[1]中的空间测边网算例，P₁、P₂、…、P₉为9个已知点坐标，9个已知点到2个未知点P₁₀、P₁₁(模拟真值分别为(0 m，0 m，0 m)和(7 m，10 m，-5 m))的观测距离也已知。2个未知点之间的观测距离为d_{10, 11}=13.107 85 m，各距离为等精度观测，中误差为±0.001 m。根据19个等精度观测距离确定2个未知点坐标(具体数据未列出，详情可见文献[1])。在该算例中，法矩阵条件数为4.601 6×10³，属于病态问题。为验证新方法的可行性，分别用LS估计、TSVD、岭估计和新方法进行计算并比较运算结果均方误差和差值范数的优劣，详细结果见表 2。

4.3 算例分析

根据表 1、2，在法矩阵A严重病态下，将各种方法求得的参数估值与参数真值对比可以发现，LS估计与真值相差最大，新方法的参数估值与真值更为接近，相比其他3种方法更为可靠。差值范数$\left\| {\Delta \mathit{\boldsymbol{\hat X}}} \right\|$可以更好地体现估值与真值的接近程度。算例1中4种方法的差值范数依次为1.308 8、0.813 0、0.854 7、0.494 5；算例2中4种方法的差值范数依次为10.353 3、0.287 5、0.166 7、0.128 1。可以看出，LS估计差值范数最大，新方法的差值范数最小，说明新方法估值结果更接近真值。然而，测量中参数真值往往是不可知的，差值范数很难求得，所以需要更为科学的评判标准评价参数估计。

许多学者^{[15-16, 19]}都将均方误差作为衡量病态方程求解参数估计的指标，即均方误差越小，说明参数估计越优。算例1中4种方法的均方误差MSE(${\mathit{\boldsymbol{\hat X}}}$)依次为42.877 6、4.309 3、0.991 7、0.395 2；算例2中4种方法的均方误差依次为18.282 3、0.293 0、0.079 6、0.053 0。可以看出，由于系数矩阵病态，LS估计的均方误差远大于另外3种方法。为克服病态带来的影响，TSVD、岭估计、新方法都在一定程度上抑制法矩阵A^TA求逆时的不稳定，主要表现在对法矩阵的所有特征值增加一个数，使得法矩阵条件数降低，病态性削弱。对比3种方法求解出的估值的均方误差可以看出，新方法的均方误差比TSVD、岭估计低，所以可以认为，本文新方法的参数估值更加可靠。并且在构成新正则化矩阵之后，新方法不需要再求正则化参数而直接计算，计算方便、高效。2个算例均表明，新方法参数估值可靠、稳定，方法有效、可行，在处理病态问题时具有一定优势。

5 结语

本文提出一种构造正则化矩阵的新方法，利用单位阵和一个特定对角矩阵进行组合，在计算时直接解算不需要求正则化参数，因此不会对解算效率产生影响。

经过推理，并利用2组模拟数据进行验证，在正则化参数小于1的情况下求得本文的正则化矩阵，参数估值的均方误差小于LS估计、TSVD和岭估计，得到的运算结果的差值范数也小于LS估计、TSVD和岭估计，并且得到的参数估值也更优，处理病态问题采用本文所提正则化矩阵进行解算具有一定的优势。

当正则化参数等于1时，求得的正则化矩阵就变成单位阵，新方法就与岭估计相同，且求出的参数估值、均方误差和差值范数也相同。因此在这种条件下，岭估计是本文新方法的特例。

在正则化参数大于1时，式(19)~(25)证明了新方法正则化解的均方误差将大于岭估计，表明新方法劣于岭估计。因此，本文新方法适用于参数小于1的情况，在参数大于1的情况下如何选择合理的正则化矩阵有待进一步研究。