-
病态观测方程广泛存在于大地测量解算模型中，如卫星定位[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积分方程和病态测边网为算例验证公式的正确性。
-
常用的G-M模型为：
$$ L = AX + e $$ (1) 式中，L∈Rm为观测向量；e∈Rm为误差向量；A∈Rm×n为系数矩阵；X∈Rn为待估参数。其统计性质如下：
$$ \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阵的平方根分解矩阵。令:
$$ \bar L = \sqrt P L,\bar A = \sqrt P A,\bar e = \sqrt P e $$ 则式(3)变为：
$$ \bar L = \bar A{\rm{ }}X + {\rm{ }}\bar e $$ (4) 这说明不等权观测方程(1)可变换为等权观测方程(4)。因此，为了推导方便，可设式(1)的权阵为单位阵。根据最小二乘（least square，LS）准则可导出式(1)的最小二乘参数估值XLS及其协方差阵D(XLS)：
$$ \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∈Rm×m, V∈Rn×n且满足UTU=UUT=Im, VTV=VVT=In，Im和In分别表示m×m和n×n阶单位阵; Σ∈Rm×n, $\Sigma = \left[ {\begin{array}{*{20}{c}} D\\ 0 \end{array}} \right]$, D为n×n阶对角阵，其对角元λi(1≤i≤n)为A的奇异值，且按降序排列。若将U、V用列向量表示成U=(u1 u2⋯um), V=(v1 v2⋯vn)，则A可展开成:
$$ A = \sum\limits_{i = {\rm{ }}1}^n {{\lambda _i}{u_i}{v_i}^{\rm{T}}} $$ (7) 将式（7）代入式(5)，可得XLS和D(XLS)的谱分解式[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) -
当系数矩阵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∈Rm×n, ${\Sigma _k} = \left[ {\begin{array}{*{20}{c}} {{D_k}}\\ 0 \end{array}} \right]$，其中Dk为对角阵diag(λ1 λ2⋯λ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∈Rn×n为对角阵；其对角元为diag($\frac{1}{{\lambda _1^2}}\frac{1}{{\lambda _2^2}} \cdots \frac{1}{{\lambda _k^2}}0 \cdots 0$)；V2为V的第k+1~n列向量构成的矩阵，即V2=(vk+1 vk+2…vn)。
-
记U1=(u1 u2⋯uk), U2=(uk+1 uk+2⋯um), 顾及式(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正则化解XT代替。
-
影响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。
-
第一类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) 式中,
$$ X = {[\begin{array}{*{20}{c}} {f({x_1})}&{f({x_2})}& \cdots &{f({x_{51}})} \end{array}]^{\rm{T}}} $$ $$ L = {[\begin{array}{*{20}{c}} {z({y_1})}&{z({y_2})}& \cdots &{z({y_{51}})} \end{array}]^{\rm{T}}} $$ $$ A = \frac{1}{{100}}{\left[ {\begin{array}{*{20}{c}} {K({x_1}, {y_1})}& \cdots &{K({x_1}, {y_{201}})}\\ {2K({x_2}, {y_1})}& \cdots &{2K({x_2}, {y_{201}})}\\ \vdots &{}& \vdots \\ {2K({x_{50}}, {y_1})}& \cdots &{2K({x_{50}}, {y_{201}})}\\ {K({x_{51}}, {y_1})}& \cdots &{K({x_{51}}, {y_{201}})} \end{array}} \right]^2} $$ 对L模拟零均值，单位权方差为0.01，满足正态分布的偶然误差，法矩阵ATA的条件数为6.463 8×1012, 严重病态。分别采用最小二乘法和TSVD正则化法求解参数并与真实值相比，图 1为使用L曲线法确定截断参数的示意图，由图 1可知，距离曲线曲率最大的点最近的散点对应的截断值为16，故截断参数k取16。
图 1 L曲线法确定截断参数的示意图（数值算例）
Figure 1. Illustration of Determining the Truncation Parameter by Using L-curve Method(Numerical Example)
两种算法解算的参数值分别见图 2和图 3。可知，由于法矩阵病态性的影响，最小二乘解极不可靠，严重偏离真值；而TSVD正则化法截掉了小的奇异值，有效改善了法矩阵的病态性，其解与真实值非常接近。
图 2 最小二乘解结果
Figure 2. Result of Least Squares Solution
图 3 TSVD正则化解结果与真值
Figure 3. Result of TSVD Regularization and Actual Values
为验证本文提出的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%, 二者的计算结果精度相当，且均优于传统公式计算的结果。
图 4 3种算法算得的1 000个单位权方差对比（数值算例）
Figure 4. Comparison of 1 000 Unit Weight Variance Estimated by Three Algorithms(Numerical Example)
-
为进一步验证本文公式的正确性，采用文献[24]中的病态测边网算例加以验证。该算例中共有11个点位, 其中9个已知点P1、P2⋯P9，2个未知点P10, P11, 图 5为测边网的点位平面分布图，表 1给出了已知点坐标和距离观测值。两个未知点的真实坐标分别为(0, 0, 0)和(7, 10, -5)，其观测距离d10, 11=13.107 8 m, 各观测距离均为等精度观测, 两个未知点坐标可根据19个距离观测值建立误差方程平差得到。
图 5 空间测边网的点位平面分布图
Figure 5. Point Position Distribution Map of the Space Net in XY Plane
表 1 控制点的坐标及距离观测值
Table 1. Coordinates and Distance Observations of Control Points
点号 坐标 距离观测值 X Y Z di, 10 di, 11 P1 23.000 10.000 0.010 25.078 69 16.765 17 P2 −10.000 9.990 0.000 14.134 51 17.719 65 P3 35.000 10.010 −0.010 36.415 88 28.442 94 P4 100.000 19.990 0.005 101.479 43 93.168 39 P5 −36.000 10.005 0.000 37.364 22 43.299 05 P6 0.000 10.010 −0.005 10.010 04 8.600 60 P7 56.000 9.995 0.010 56.996 06 49.256 18 P8 −15.000 10.015 −0.010 18.035 90 22.559 66 P9 −1.700 10.008 0.015 10.150 63 10.043 82 该算例中，法矩阵ATA的条件数为4.585 1×103, 存在病态。利用该数据，分别采用最小二乘法和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正则化法将较小的奇异值截掉，改善了法矩阵的病态性，其解算结果与真值之差明显小于最小二乘解的结果。
图 6 用L曲线法确定截断参数的示意图(病态测边网算例)
Figure 6. Illustration of Determining the Truncation Parameter by Using L-curve Method(Example of Ill-Posed Trilateration Network)
表 2 病态测边网计算结果
Table 2. Results of Ill-Posed Trilateration Network
方法 ${{\hat X}_{10}}$ ${{\hat Y}_{10}}$ ${{\hat Z}_{10}}$ ${{\hat X}_{11}}$ ${{\hat Y}_{11}}$ ${{\hat Z}_{11}}$ $\left\| {\hat X - X} \right\|$ κ 真值 0.000 0 0.000 0 0.000 0 7.000 0 10.000 0 −5.000 0 0.000 0 0.000 0 最小二乘解 −0.041 8 −0.054 1 11.004 5 9.031 7 11.027 7 −4.302 7 11.259 4 0.608 1 TSVD正则化解 0.036 7 0.048 1 0.411 2 6.912 5 10.376 8 −4.872 4 0.581 9 0.043 5 利用本算例数据，模拟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%，二者的计算结果精度相当，且均优于传统公式计算的结果。
图 7 种算法算得的1 000个单位权方差对比(病态测边网算例)
Figure 7. Comparison of 1 000 Unit Weight Variance Estimated by Three Algorithms(Example of Ill-Posed Trilateration Network)
-
截断奇异值法通过截掉病态观测方程系数矩阵的小奇异值来改善模型的病态性，以提高参数估值的稳定性和精度。然而，截掉小奇异值后，由于改变了观测方程系数矩阵的结构，导致参数估值和残差估值均有偏，若再用传统的单位权方差估计式估计单位权方差，其结果必定有偏。本文基于残差二次型的数学期望公式，利用有偏的残差导出了TSVD正则化解的无偏单位权方差估计式。由于该公式中包含参数真值，在实际计算时可用TSVD正则化解代替。
算例分析结果表明，传统公式计算的单位权方差明显偏离真值，1 000次模拟实验结果的均值与真值分别相差11.25%（数值算例）和15.33%（病态测边网算例），而本文公式和文献[16]提出的Tikhonov正则化解的单位权方差无偏估计式计算的单位权方差精度相当，数值算例中1 000次模拟实验结果的均值与真值分别相差0.35%和0.41%，病态测边网算例中1 000次模拟实验结果的均值与真值分别相差2.09%和1.82%，显然优于传统公式所估的结果。
-
摘要: 截断奇异值（truncated singular value decomposition，TSVD）法通过截掉病态观测方程系数矩阵的小奇异值来改善模型的病态性，提高参数估值的稳定性和精度。然而，截除小奇异值后，改变了观测方程的结构，不仅参数估值有偏，残差估值也是有偏的；因此，其单位权方差不能用传统的估计公式计算。针对此，导出了TSVD正则化解的单位权方差无偏公式，并以第一类Fredholm积分方程和病态测边网为算例验证了公式的正确性。Abstract: The truncated singular value method (TSVD) improves the morbidity of model by truncating small singular values of the ill-posed observational equation matrix and increases the stability and accuracy of the parameter estimation. However, the structure of observation equation has been changed after truncating small singular value, which makes parameter valuation and residual biased, the unit weight variance cannot be calculated by using traditional estimation formula. This paper derives the unbiased formula of unit weight variance for TSVD regularization, and uses the first Fredholm integral equation and ill-posed trilateration network as examples to verify the correctness of the formula.
-
图 1 L曲线法确定截断参数的示意图（数值算例）
Figure 1. Illustration of Determining the Truncation Parameter by Using L-curve Method(Numerical Example)
图 4 3种算法算得的1 000个单位权方差对比（数值算例）
Figure 4. Comparison of 1 000 Unit Weight Variance Estimated by Three Algorithms(Numerical Example)
图 5 空间测边网的点位平面分布图
Figure 5. Point Position Distribution Map of the Space Net in XY Plane
图 6 用L曲线法确定截断参数的示意图(病态测边网算例)
Figure 6. Illustration of Determining the Truncation Parameter by Using L-curve Method(Example of Ill-Posed Trilateration Network)
图 7 种算法算得的1 000个单位权方差对比(病态测边网算例)
Figure 7. Comparison of 1 000 Unit Weight Variance Estimated by Three Algorithms(Example of Ill-Posed Trilateration Network)
表 1 控制点的坐标及距离观测值
Table 1. Coordinates and Distance Observations of Control Points
点号 坐标 距离观测值 X Y Z di, 10 di, 11 P1 23.000 10.000 0.010 25.078 69 16.765 17 P2 −10.000 9.990 0.000 14.134 51 17.719 65 P3 35.000 10.010 −0.010 36.415 88 28.442 94 P4 100.000 19.990 0.005 101.479 43 93.168 39 P5 −36.000 10.005 0.000 37.364 22 43.299 05 P6 0.000 10.010 −0.005 10.010 04 8.600 60 P7 56.000 9.995 0.010 56.996 06 49.256 18 P8 −15.000 10.015 −0.010 18.035 90 22.559 66 P9 −1.700 10.008 0.015 10.150 63 10.043 82 
下载: 导出CSV
表 2 病态测边网计算结果
Table 2. Results of Ill-Posed Trilateration Network
方法 ${{\hat X}_{10}}$ ${{\hat Y}_{10}}$ ${{\hat Z}_{10}}$ ${{\hat X}_{11}}$ ${{\hat Y}_{11}}$ ${{\hat Z}_{11}}$ $\left\| {\hat X - X} \right\|$ κ 真值 0.000 0 0.000 0 0.000 0 7.000 0 10.000 0 −5.000 0 0.000 0 0.000 0 最小二乘解 −0.041 8 −0.054 1 11.004 5 9.031 7 11.027 7 −4.302 7 11.259 4 0.608 1 TSVD正则化解 0.036 7 0.048 1 0.411 2 6.912 5 10.376 8 −4.872 4 0.581 9 0.043 5
下载: 导出CSV
-
[1] 周忠谟, 易杰军. GPS卫星测量原理及其应用[M].北京:测绘出版社, 1992 Zhou Zhongmo, Yi Jiejun. GPS Satellite Measurement Principle and Its Application[M]. Beijing:Surveying and Mapping Press, 1992 [2] Li B, Shen Y, Feng Y. Fast GNSS Ambiguity Resolution as an Ill-posed Problem[J]. Journal of Geodesy, 2010, 84(11):683-698 doi: 10.1007/s00190-010-0403-5 [3] 党亚民, 陈俊勇.大地测量反演模型优化问题的研究[J].地壳形变与地震, 1998, 18(2):35-40 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199800045082 Dang Yaming, Chen Junyong. Study on Model Optimization for Geodetic Inversion[J]. Crustal Deformation and Earthquake, 1998, 18(2):35-40 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199800045082 [4] 顾勇为, 归庆明, 张璇, 等.大地测量与地球物理中病态性问题的正则化迭代解法[J].测绘学报, 2014, 43(4):331-336 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=chxb201404002 Gu Yongwei, Gui Qingming, Zhang Xuan, et al. Iterative Solution of Regularization to Ill-conditioned Problems in Geodesy and Geophysics[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(4):331-336 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=chxb201404002 [5] 徐新禹, 李建成, 王正涛, 等. Tikhonov正则化方法在GOCE重力场求解中的模拟研究[J].测绘学报, 2010, 39(5):29-34 http://d.old.wanfangdata.com.cn/Periodical/chxb201005005 Xu Xingyu, Li Jiancheng, Wang Zhengtao, et al. The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Mission[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(5):29-34 http://d.old.wanfangdata.com.cn/Periodical/chxb201005005 [6] 邓凯亮, 黄谟涛, 暴景阳, 等.向下延拓航空重力数据的Tikhonov双参数正则化法[J].测绘学报, 2011, 40(6):690-696 http://www.cqvip.com/QK/90069X/201106/40330127.html Deng Kailiang, Huang Motao, Bao Jingyang, et al. Tikhonov Two-parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6):690-696 http://www.cqvip.com/QK/90069X/201106/40330127.html [7] Zhou W. Normalized Full Gradient of Full Tensor Gravity Gradient Based on Adaptive Iterative Tikhonov Regularization Downward Continuation[J]. Journal of Applied Geophysics, 2015, 118:75-83 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=cc060d9a4883c11e1dbfdbb1b6a82925 [8] Tikhonov A N. Solution of Incorrectly Formulated Problems and the Regularization Method[J]. Soviet Math Dokl, 1963, 5(4):501-504 https://www.scirp.org/(S(351jmbntvnsjt1aadkposzje))/reference/ReferencesPapers.aspx?ReferenceID=1088489 [9] Tikhonov A N. On the Solution of Ill-Posed Problems and the Method of Regularization[J]. Doklady Akademii Nauk, 1963, 151(3):501-504 http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_math%2f0511354 [10] Hoerl A E, Kennard R W. Ridge Regression:Applications to Nonorthogonal Problems[J]. Technometrics, 1970, 12(1):69-82 doi: 10.1080/00401706.1970.10488635 [11] Hansen P C. The Truncated SVD as a Method for Regularization[J]. BIT Numerical Mathematics, 1987, 27(4):534-553 doi: 10.1007-BF01937276/ [12] Xu P. Truncated SVD Methods for Discrete Linear Ill-Posed Problems[J]. Geophysical Journal International, 1998, 135(2):505-514 doi: 10.1046-j.1365-246X.1998.00652.x/ [13] Sourbron S, Luypaert R, van Schuerbeek P, et al. Choice of the Regularization Parameter for Perfusion Quantification with MRI[J]. Physics in Medicine & Biology, 2004, 49(14):3307-3031 https://ui.adsabs.harvard.edu/abs/2004PMB....49.3307S/abstract [14] 林东方, 朱建军, 宋迎春.顾及截断偏差影响的TSVD截断参数确定方法[J].测绘学报, 2017, 46(6):679-688 http://d.old.wanfangdata.com.cn/Periodical/chxb201706003 Lin Dongfang, Zhu Jianjun, Song Yingchun. Truncation Method for TSVD with Account of Truncation Bias[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(6):679-688 http://d.old.wanfangdata.com.cn/Periodical/chxb201706003 [15] Wu Y, Zhang Y, Zhang Y, et al. The Regularization Method Based on TSVD for Forward-Looking Radar Angular Superresolution[C]. IEEE Geoscience and Remote Sensing Symposium (IGARSS), Fort Worth, Texas, USA, 2017 [16] 沈云中, 刘大杰.正则化解的单位权方差无偏估计公式[J].武汉大学学报·信息科学版, 2002, 27(6):604-606 http://ch.whu.edu.cn/CN/abstract/abstract5018.shtml Shen Yunzhong, Liu Dajie. Unbiased Estimation Formula of Unit Weight Mean Square Error in Regularization Solution[J]. Geomatics and Information Science of Wuhan University, 2002, 27(6):604-606 http://ch.whu.edu.cn/CN/abstract/abstract5018.shtml [17] 张俊, 独知行, 杜宁, 等.补偿最小二乘解的单位权方差的无偏估计[J].大地测量与地球动力学, 2014, 34(1):153-156 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dkxbydz201401035 Zhang Jun, Du Zhixing, Du Ning, et al. Unbiased Estimation Formula of Unit Weight Variance in Solution by Penalized Least Square[J]. Journal of Geodesy and Geodynamics, 2014, 34(1):153-156 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dkxbydz201401035 [18] 沈云中, 许厚泽.不适定方程正则化算法的谱分解式[J].大地测量与地球动力学, 2002, 22(3):10-14 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dkxbydz200203003 Shen Yunzhong, Xu Houze. Spectral Decomposition Formula of Regularization Solution for Ill-Posed Equation[J]. Journal of Geodesy and Geodynamics, 2002, 22(3):10-14 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dkxbydz200203003 [19] Hansen P C. Analysis of Discrete Ill-Posed Problems by Means of the L-curve[J]. SIAM Review, 1992, 34(4):561-580 doi: 10.1137/1034115 [20] Hansen P C. Regularization Tools:A MATLAB Package for Analysis and Solution of Discrete IllPosed Problems[J]. Numerical Algorithms, 1994, 6(1):1-35 doi: 10.1007/BF02149761 [21] 游为, 范东明, 黄强.卫星重力反演的短弧长积分法研究[J].地球物理学报, 2011, 54(11):2745-2752 http://d.old.wanfangdata.com.cn/Periodical/dqwlxb201111004 You Wei, Fan Dongming, Huang Qiang. Analysis of Short-arc Integral Approach to Recover the Earth's Gravitational Field[J]. Chinese Journal of Geophysics, 2011, 54(11):2745-2752 http://d.old.wanfangdata.com.cn/Periodical/dqwlxb201111004 [22] 宋淑云.基于重力异常的密度界面反演研究[D].合肥: 中国科学技术大学, 2008 Song Shuyun. Density Interface Inversion Based on Gravity Anomaly[D]. Hefei: University of Science and Technology of China, 2008 [23] 姜刚.龙门山断裂带地壳形变特征及地震断层参数反演研究[D].西安: 长安大学, 2016 Jiang Gang. Research in Crustal Deformation in Longmenshan Fault Zone and Earthquake Fault Parameters Inversion[D]. Xi'an: Chang'an University, 2016 [24] 郭建锋.测量平差系统病态性的诊断与处理[D].郑州: 信息工程大学, 2002 Guo Jianfeng. The Diagnosis and Process of Ill-Conditioned Adjustment System[D]. Zhengzhou: Information Engineering University, 2002 -
点击查看大图
图(7) / 表(2)
计量
- 文章访问数: 1094
- HTML全文浏览量: 333
- PDF下载量: 36
- 被引次数: 0
