病态观测方程广泛存在于大地测量解算模型中，如卫星定位^[1-2]、大地测量反演^[3-4]、重力场向下延拓^[5-7]等。观测方程的病态性使最小二乘解对观测误差极为敏感，导致其参数估值极不稳定。为了求得稳定的参数估值，需要根据病态观测方程的特点，构建正则化解。常用的正则化方法有Tikhonov正则化、岭估计和截断奇异值（truncated singular value decomposition，TSVD）法。Tikhonov正则化方法通过构造稳定泛函准则并引入正则化参数来求解稳定的参数估值^[8-9]；岭估计由Hoerl和Kennard^[10]于1970年提出，它是Tikhonov正则化方法的特例；TSVD法是基于奇异值分解的一种病态观测方程的解算方法，通过截掉系数矩阵的小奇异值避免放大高频观测误差对参数估值的影响, 仅利用大的奇异值及相应的特征向量构建参数估值^[11-12]。由于TSVD法只需要确定截断的奇异值数k, 因此广泛应用于病态观测方程的解算。

然而，TSVD法截掉了系数矩阵的小奇异值，尽管改善了观测方程的病态性，但改变了观测方程系数矩阵的结构，导致TSVD正则化解的参数估值和残差均是有偏的，若还用传统单位权方差估计公式计算单位权方差，其结果必定也是有偏的。目前，关于TSVD法的研究主要集中于分析参数估值的有偏性和解算更好的参数估值^[13-15], 应用其有偏的残差估值计算无偏单位权方差则未见相关的研究文献。文献[16-17]根据有偏残差的二次型期望公式分别导出了病态方程Tikhonov正则化解的无偏单位权方差估计式和半参数模型补偿最小二乘解的无偏单位权方差估计式。本文根据类似的思路，利用TSVD法的有偏残差导出了其正则化解无偏单位权方差估计式，并以第一类Fredholm积分方程和病态测边网为算例验证公式的正确性。

1 病态模型的最小二乘解及其谱分解式

常用的G-M模型为：

$ L = AX + e $

(1)

式中，L∈R^m为观测向量；e∈R^m为误差向量；A∈R^m×n为系数矩阵；X∈Rⁿ为待估参数。其统计性质如下：

$ \left\{ \begin{array}{l} E(L) = AX\\ E(e) = 0\\ D(e) = \sigma _0^2{P^{ - {\rm{ }}1}} \end{array} \right. $

(2)

式中，$\sigma _0^2$为先验单位权方差；P为权阵。

对式(1)作如下变换：

$ \sqrt P L = \sqrt P AX + \sqrt P {\rm{ }}e $

(3)

式中，$\sqrt P $为P阵的平方根分解矩阵。令:

则式(3)变为：

$ \bar L = \bar A{\rm{ }}X + {\rm{ }}\bar e $

(4)

这说明不等权观测方程(1)可变换为等权观测方程(4)。因此，为了推导方便，可设式(1)的权阵为单位阵。根据最小二乘（least square，LS）准则可导出式(1)的最小二乘参数估值X_LS及其协方差阵D(X_LS)：

$ \left\{ \begin{array}{l} {X_{{\rm{LS}}}} = {\left( {{A^{\rm{T}}}A} \right)^{ - {\rm{ }}1}}{A^{\rm{T}}}, \\ D\left( {{X_{{\rm{LS}}}}} \right) = \sigma _0^2{\left( {{A^{\rm{T}}}A} \right)^{ - {\rm{ }}1}} \end{array} \right. $

(5)

对A进行奇异值分解：

$ A = U\Sigma {V^{\rm{T}}} $

(6)

式中，U∈R^m×m, V∈R^n×n且满足U^TU=UU^T=I_m, V^TV=VV^T=I_n，I_m和I_n分别表示m×m和n×n阶单位阵; Σ∈R^m×n, $\Sigma = \left[ {\begin{array}{*{20}{c}} D\\ 0 \end{array}} \right]$, D为n×n阶对角阵，其对角元λ_i(1≤i≤n)为A的奇异值，且按降序排列。若将U、V用列向量表示成U=(u₁ u₂⋯u_m), V=(v₁ v₂⋯v_n)，则A可展开成:

$ A = \sum\limits_{i = {\rm{ }}1}^n {{\lambda _i}{u_i}{v_i}^{\rm{T}}} $

(7)

将式（7）代入式(5)，可得X_LS和D(X_LS)的谱分解式^[18]：

$ \left\{ \begin{array}{l} {X_{{\rm{LS}}}} = \sum\limits_{i = {\rm{ }}1}^n {\frac{{{{\begin{array}{*{20}{l}} u \end{array}}_i}^{\rm{T}}L}}{{{\lambda _i}}}} {v_i} = {\rm{ }}(\sum\limits_{i = {\rm{ }}1}^n {\frac{{\begin{array}{*{20}{l}} {{v_i}u} \end{array}_i^{\rm{T}}}}{{{\lambda _i}}}} )L = {A^ + }L\\ D\left( {{X_{{\rm{LS}}}}} \right) = \sigma _0^2\sum\limits_{i = {\rm{ }}1}^n {\frac{{\begin{array}{*{20}{l}} {{v_i}v_i^{\rm{T}}} \end{array}}}{{\lambda _i^2}}} \end{array} \right. $

(8)

2 TSVD正则化解及其单位权方差无偏估计 2.1 TSVD正则化解

当系数矩阵A病态时，其奇异值单调趋向于零。由式(8)可知，由于小奇异值出现在分母上，显著放大了观测误差对参数估值及其方差的影响，导致最小二乘解极不稳定。TSVD法的基本思想是将引起解不稳定的小的奇异值截掉，保留较大的奇异值，重新构造系数矩阵，以避免较小的奇异值对参数估值及其方差的影响。根据A的奇异值分解式(7)，保留前k项，截掉后n−k项，得到新的系数矩阵：

$ {A_k} = \sum\limits_{i = {\rm{ }}1}^k {{\lambda _i}{u_i}{v_i}^{\rm{T}}} = {\rm{ }}U{\Sigma _k}{V^{\rm{T}}} $

(9)

式中，U、V意义同式(6)；Σ_k∈R^m×n, ${\Sigma _k} = \left[ {\begin{array}{*{20}{c}} {{D_k}}\\ 0 \end{array}} \right]$，其中D_k为对角阵diag(λ₁ λ₂⋯λ_k 0⋯0)。根据式(8)可得TSVD正则化解及其协方差阵：

$ \left\{ \begin{array}{l} {X_T} = \sum\limits_{i = {\rm{ }}1}^k {\frac{{\begin{array}{*{20}{l}} u \end{array}_i^{\rm{T}}L}}{{{\lambda _i}}}} {v_i} = \left( {\sum\limits_{i = {\rm{ }}1}^k {\frac{{\begin{array}{*{20}{l}} {{v_i}u_i^{\rm{T}}} \end{array}}}{{{\lambda _i}}}} } \right)L = A_k^ + L\\ D\left( {{X_T}} \right) = \sigma _0^2\sum\limits_{i = {\rm{ }}1}^K {\frac{{\begin{array}{*{20}{l}} {{v_i}v_i^{\rm{T}}} \end{array}}}{{\lambda _i^2}}} \end{array} \right. $

(10)

式中，$A_k^ + = \sum\limits_{i = {\rm{ }}1}^k {\frac{{\begin{array}{*{20}{l}} {{v_i}u_i^{\rm{T}}} \end{array}}}{{{\lambda _i}}}} = V\sum {_k^ + } {U^{\rm{T}}}$; $\sum {_k^ + } = \left[ {\begin{array}{*{20}{c}} {D_k^{ - 1}}\\ 0 \end{array}} \right]$，其中${D_k^{ - 1}}$为对角阵diag($\frac{1}{{{\lambda _1}}}\frac{1}{{{\lambda _2}}} \cdots \frac{1}{{{\lambda _k}}}0 \cdots 0$)。由式(10)可知，TSVD法通过截掉小的奇异值来消除小奇异值对方差的扩大，有效减少估值的方差，因此其参数估值更为稳定。但TSVD法由于改变了观测方程系数矩阵的结构，所求的参数解是有偏的，其偏差和用于评定精度的均方误差为：

$ \left\{ \begin{array}{l} {\rm{Bias}}({X_T}) = E({X_T}) - X = A_k^ + E(L) - X = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (A_k^ + A - {I_n})X = {\rm{ }} - {V_2}V_2^{\rm{T}}X\\ {\rm{MSE }}({X_T}) = E[{({X_T} - X)^{\rm{T}}}({X_T} - X){\rm{ }}] = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D({X_T}) + {\rm{ Bias}}({X_T}){\rm{ Bias}}{({X_T})^{\rm{T}}} = \\ {\rm{ }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\sigma _0^2V{\Sigma _T}{V^{\rm{T}}}} \end{array} + {V_2}V_2^{\rm{T}}X{X^{\rm{T}}}{V_2}V_2^{\rm{T}} \end{array} \right. $

(11)

式中，Σ_T∈R^n×n为对角阵；其对角元为diag($\frac{1}{{\lambda _1^2}}\frac{1}{{\lambda _2^2}} \cdots \frac{1}{{\lambda _k^2}}0 \cdots 0$)；V₂为V的第k+1~n列向量构成的矩阵，即V₂=(v_k+1 v_k+2…v_n)。

2.2 正则化解的单位权方差无偏估计

记U₁=(u₁ u₂⋯u_k), U₂=(u_k+1 u_k+2⋯u_m), 顾及式(6)、(7)和式(9)，可得：

$ \begin{array}{l} AA_k^ + - {I_m} = U\Sigma {V^{\rm{T}}}V\Sigma _k^ + {U^{\rm{T}}} - {I_m} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} U\Sigma \Sigma _k^ + {U^{\rm{T}}} - {I_m} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {{U_1}}&{{U_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{I_k}}&{}\\ {}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {U_1^{\rm{T}}}\\ {U_2^{\rm{T}}} \end{array}} \right] - {I_m} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {U_1}U_1^{\rm{T}} - {I_m} = - {U_2}U_2^{\rm{T}} \end{array} $

(12)

将式(10)代入式(1)，并顾及式(12)得TSVD法的残差为：

$ {{\hat e}_k} = A{X_T} - L = \left( {AA_k^ + - {I_m}} \right)L = - {U_2}U_2^{\rm{T}}L $

(13)

其数学期望为：

$ E\left( {{{\hat e}_k}} \right) = - {U_2}U_2^{\rm{T}}AX $

(14)

根据式（13）和协方差传播定律，可得残差的协方差阵为：

$ D\left( {{{\hat e}_k}} \right) = {U_2}U_2^{\rm{T}}{U_2}U_2^{\rm{T}}\sigma _0^2 = \sigma _0^2{U_2}U_2^{\rm{T}} $

(15)

根据式（14）和式（15）并利用矩阵迹的性质得：

$ \left\{ \begin{array}{l} {\rm{tr}}\left[ {E\left( {{{\hat e}_k}} \right)E{{\left( {{{\hat e}_k}} \right)}^{\rm{T}}}} \right] = E{\left( {{{\hat e}_k}} \right)^{\rm{T}}}E\left( {{{\hat e}_k}} \right) = \\ {\rm{ }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {X^{\rm{T}}}{A^{\rm{T}}}{U_2}U_2^{\rm{T}}{U_2}U_2^{\rm{T}}AX = \\ {\rm{ }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{ }}{X^{\rm{T}}}{A^{\rm{T}}}{U_2}U_2^{\rm{T}}AX\\ {\rm{tr}}D\left( {{{\hat e}_k}} \right) = \sigma _0^2{\rm{tr}}\left( {{U_2}U_2^{\rm{T}}} \right) = \sigma _0^2{\rm{tr}}\left( {U_2^{\rm{T}}{U_2}} \right) = \\ {\rm{ }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma _0^2{\rm{tr}}\left( {{I_{m - k}}} \right) = \sigma _0^2\left( {m - k} \right) \end{array} \right. $

(16)

已知二次型的数学期望公式为：

$ \begin{array}{l} E\left( {\hat e_k^{\rm{T}}{{\hat e}_k}} \right) = {\rm{tr}}\left[ {E\left( {{{\hat e}_k}\hat e_k^{\rm{T}}} \right)} \right] = \\ {\rm{ }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{tr}}\left[ {D\left( {{{\hat e}_k}} \right)} \right] + {\rm{tr}}\left[ {E\left( {{{\hat e}_k}} \right)E{{\left( {{{\hat e}_k}} \right)}^{\rm{T}}}} \right] \end{array} $

(17)

将式(16)代入式(17)右端，并将期望值$E\left( {\hat e_k^{\rm{T}}{{\hat e}_k}} \right)$用估值${\hat e_k^{\rm{T}}{{\hat e}_k}}$代替，$\sigma _0^2$也相应地用估值$\hat \sigma _0^2$代替得：

$ \hat e_k^{\rm{T}}{{\hat e}_k} = \hat \sigma _0^2\left( {m - k} \right) + {X^{\rm{T}}}{A^{\rm{T}}}{U_2}U_2^{\rm{T}}AX $

(18)

解得单位权方差的无偏估计公式为：

$ \hat \sigma _0^2 = \frac{{\hat e_k^{\rm{T}}{{\hat e}_k} - {X^{\rm{T}}}{A^{\rm{T}}}{U_2}U_2^{\rm{T}}AX}}{{m - k}} $

(19)

由于无偏公式(19)包含有真值X，而真值本身是未知的，因此若用估值代替真值进行计算，求得的单位权方差并非真正无偏的。在实际应用中，真值可用TSVD正则化解X_T代替。

2.3 截断参数k的确定

影响TSVD解算结果的关键因素是截断参数k的选择。本文采用L曲线法^[19]确定TSVD正则化法的截断参数，该法的基本思想是：分别计算取不同截断参数下的正则化解范数$\left\| {{{\hat X}_k}} \right\|$及其残差范数$\left\| {{{\hat e}_k}} \right\|$, 得到若干散点($\left\| {{{\hat e}_k}} \right\|$, $\left\| {{{\hat X}_k}} \right\|$), 经拟合后得到一条曲线，因曲线形状呈L形状，故称为L曲线，离曲线上曲率最大点最近的散点所对应的截断值即为所求的截断参数。

由于奇异值的跨度通常为多个数量级，故对横纵坐标取对数可更容易观察L曲线。记$\eta = {\left\| {{{\hat X}_R}} \right\|^2}$, $\rho = {\left\| {{{\hat e}_r}} \right\|^2}$, 分别对其取对数得：

$ \left\{ \begin{array}{l} \hat \eta = {\rm{lg}}\eta = 2{\rm{lg}}\left\| {{{\hat X}_R}} \right\|\\ \hat \rho = {\rm{lg}}\rho = 2{\rm{lg}}\left\| {{{\hat e}_r}} \right\| \end{array} \right. $

(20)

则L曲线由若干散点($\hat \rho /2$, $\hat \eta /2$)拟合而成，曲线上点的曲率计算公式为：

$ \kappa = 2\frac{{\hat \rho '\hat \eta '' - \hat \rho ''\hat \eta '}}{{{{\left[ {{{\left( {\hat \rho '} \right)}^2} + {{\left( {\hat \eta '} \right)}^2}} \right]}^{3/2}}}} $

(21)

式中，ρ'、η'、ρ''、η''分别表示一阶和二阶导数。计算式(21)的最大值即可得到曲线上曲率的最大点，离最大点最近的散点所对应的截断值即为所求的截断参数。文献[20]提供了一个正则化工具箱，其中包含了使用L曲线法确定截断参数的子程序l_curve。

3 算例验证 3.1 数值算例

第一类Fredholm积分方程是典型的病态问题，广泛应用于卫星重力反演^[21]、地面重力反演^[22]和形变反演^[23]，其通用表达形式如下：

$ z(y) = \smallint _a^bK(x, y)f(x){\rm{ }}dx $

(22)

本文采用文献[16]中的算例，核函数取为：

$ K(x, y) = \frac{{1.0}}{{\begin{array}{*{20}{l}} {1.0{\rm{ }} + {\rm{ }}100{{(y - x)}^2}} \end{array}}} $

(23)

精确解函数f(x)为：

$ f(x) = \frac{{{\rm{exp}}({\beta _1}) + {\rm{exp}}({\beta _2})}}{{10.955{\kern 1pt} {\kern 1pt} {\kern 1pt} 040{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}} - 0.052{\kern 1pt} {\kern 1pt} {\kern 1pt} 130{\kern 1pt} {\kern 1pt} {\kern 1pt} 913 $

(24)

式中, ${\beta _1} = - \frac{{{{(x - {\rm{ }}0.3{\rm{ }})}^2}}}{{0.03}}$; ${\beta _2} = - \frac{{{{(x - {\rm{ }}0.7{\rm{ }})}^2}}}{{0.03}}$; x∈[0, 1], y∈[−2, 2]。

取采样间隔为Δx=Δy=0.02，对式（22）右端用复合梯形公式离散化得：

$ \begin{array}{l} z({y_j}) = {\rm{ }}\frac{1}{{100}}{\rm{ }}[K({x_1},{y_j})f({x_1}) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{2}}\sum\limits_{i = 2}^{50} K ({x_1},{y_j})f({x_i}) + K({x_{51}},{y_j})f({x_{51}}) \end{array} $

(25)

式中, j=1, 2⋯201。将式(25)改写成矩阵形式得：

$ AX = L $

(26)

式中,

对L模拟零均值，单位权方差为0.01，满足正态分布的偶然误差，法矩阵A^TA的条件数为6.463 8×10¹², 严重病态。分别采用最小二乘法和TSVD正则化法求解参数并与真实值相比，图 1为使用L曲线法确定截断参数的示意图，由图 1可知，距离曲线曲率最大的点最近的散点对应的截断值为16，故截断参数k取16。

两种算法解算的参数值分别见图 2和图 3。可知，由于法矩阵病态性的影响，最小二乘解极不可靠，严重偏离真值；而TSVD正则化法截掉了小的奇异值，有效改善了法矩阵的病态性，其解与真实值非常接近。

为验证本文提出的TSVD正则化解的单位权方差无偏估计公式的正确性, 模拟1 000组数据，每组数据均按ε~N(0, 0.01)模拟偶然误差，采用传统不考虑偏差的公式${{\hat \sigma }^2}_{{\rm{LS}}} = \frac{{\hat e_k^{\rm{T}}{{\hat e}_k}}}{{m - n}}$(图 4中LS)、本文式（19）（图 4中TVSD）以及文献[16]提出的Tikhonov正则化解的单位权方差无偏估计式(图 4中Tikhonov)分别计算1 000个单位权方差。由图 4可知，传统公式算得的1 000个单位权方差均值为0.011 124 61，与模拟值相差11.25%；而本文公式和文献[16]中公式算得的1 000个单位权方差均值分别为0.010 035 11和0.010 041 27，与模拟值只相差0.35%和0.41%, 二者的计算结果精度相当，且均优于传统公式计算的结果。

3.2 病态测边网算例

为进一步验证本文公式的正确性，采用文献[24]中的病态测边网算例加以验证。该算例中共有11个点位, 其中9个已知点P₁、P₂⋯P₉，2个未知点P₁₀, P₁₁, 图 5为测边网的点位平面分布图，表 1给出了已知点坐标和距离观测值。两个未知点的真实坐标分别为(0, 0, 0)和(7, 10, -5)，其观测距离d_{10, 11}=13.107 8 m, 各观测距离均为等精度观测, 两个未知点坐标可根据19个距离观测值建立误差方程平差得到。

该算例中，法矩阵A^TA的条件数为4.585 1×10³, 存在病态。利用该数据，分别采用最小二乘法和TSVD正则化法求解参数，图 6为使用L曲线法确定截断参数的示意图。由图 6可知，截断参数取4, 两种算法求得的参数估值与真值的差值范数$\left\| {\Delta \hat X} \right\| = \left\| {\hat X - X} \right\|$和相对偏差范数$\kappa = \left\| {\hat X - X} \right\|{\rm{/}}\left\| {\hat X} \right\|$见表 2。由表 2可知，由于法矩阵病态性影响，最小二乘估值极不可靠，严重偏离真值；TSVD正则化法将较小的奇异值截掉，改善了法矩阵的病态性，其解算结果与真值之差明显小于最小二乘解的结果。

利用本算例数据，模拟1 000组观测误差，每组均按e~N(0, 0.05)模拟偶然误差，分别采用传统公式、本文公式和文献[16]中公式计算单位权方差，结果如图 7所示。由图 7可知，传统公式计算得到的1 000个单位权方差均值为0.057 665 89，与模拟值相差15.33%，明显偏大；而本文公式和文献[16]中公式计算得到的单位权方差均值分别为0.051 045 99和0.050 907 73, 与模拟值只相差2.09%和1.82%，二者的计算结果精度相当，且均优于传统公式计算的结果。

4 结语

截断奇异值法通过截掉病态观测方程系数矩阵的小奇异值来改善模型的病态性，以提高参数估值的稳定性和精度。然而，截掉小奇异值后，由于改变了观测方程系数矩阵的结构，导致参数估值和残差估值均有偏，若再用传统的单位权方差估计式估计单位权方差，其结果必定有偏。本文基于残差二次型的数学期望公式，利用有偏的残差导出了TSVD正则化解的无偏单位权方差估计式。由于该公式中包含参数真值，在实际计算时可用TSVD正则化解代替。

算例分析结果表明，传统公式计算的单位权方差明显偏离真值，1 000次模拟实验结果的均值与真值分别相差11.25%（数值算例）和15.33%（病态测边网算例），而本文公式和文献[16]提出的Tikhonov正则化解的单位权方差无偏估计式计算的单位权方差精度相当，数值算例中1 000次模拟实验结果的均值与真值分别相差0.35%和0.41%，病态测边网算例中1 000次模拟实验结果的均值与真值分别相差2.09%和1.82%，显然优于传统公式所估的结果。