留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

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

嵇昆浦 沈云中

嵇昆浦, 沈云中. TSVD正则化解法的单位权方差无偏估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
引用本文: 嵇昆浦, 沈云中. TSVD正则化解法的单位权方差无偏估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
JI Kunpu, SHEN Yunzhong. Unbiased Estimation of Unit Weight Variance by TSVD Regularization[J]. Geomatics and Information Science of Wuhan University, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
Citation: JI Kunpu, SHEN Yunzhong. Unbiased Estimation of Unit Weight Variance by TSVD Regularization[J]. Geomatics and Information Science of Wuhan University, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270

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

doi: 10.13203/j.whugis20180270
基金项目: 

国家自然科学基金 41731069

详细信息
    作者简介:

    嵇昆浦, 硕士生, 主要从事大地测量数据处理研究。1575540259@qq.com

    通讯作者: 沈云中, 博士, 教授。yzshen@tongji.edu.cn
  • 中图分类号: P207

Unbiased Estimation of Unit Weight Variance by TSVD Regularization

Funds: 

The National Natural Science Foundation of China 41731069

More Information
    Author Bio:

    JI Kunpu, postgraduate, specializes in geodetic data processing.1575540259@qq.com

    Corresponding author: SHEN Yunzhong, PhD, professor.yzshen@tongji.edu.cn
图(7) / 表(2)
计量
  • 文章访问数:  1094
  • HTML全文浏览量:  333
  • PDF下载量:  36
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-01-13
  • 刊出日期:  2020-04-05

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

doi: 10.13203/j.whugis20180270
    基金项目:

    国家自然科学基金 41731069

    作者简介:

    嵇昆浦, 硕士生, 主要从事大地测量数据处理研究。1575540259@qq.com

    通讯作者: 沈云中, 博士, 教授。yzshen@tongji.edu.cn
  • 中图分类号: P207

摘要: 截断奇异值(truncated singular value decomposition,TSVD)法通过截掉病态观测方程系数矩阵的小奇异值来改善模型的病态性,提高参数估值的稳定性和精度。然而,截除小奇异值后,改变了观测方程的结构,不仅参数估值有偏,残差估值也是有偏的;因此,其单位权方差不能用传统的估计公式计算。针对此,导出了TSVD正则化解的单位权方差无偏公式,并以第一类Fredholm积分方程和病态测边网为算例验证了公式的正确性。

English Abstract

嵇昆浦, 沈云中. TSVD正则化解法的单位权方差无偏估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
引用本文: 嵇昆浦, 沈云中. TSVD正则化解法的单位权方差无偏估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
JI Kunpu, SHEN Yunzhong. Unbiased Estimation of Unit Weight Variance by TSVD Regularization[J]. Geomatics and Information Science of Wuhan University, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
Citation: JI Kunpu, SHEN Yunzhong. Unbiased Estimation of Unit Weight Variance by TSVD Regularization[J]. Geomatics and Information Science of Wuhan University, 2020, 45(4): 626-632. doi: 10.13203/j.whugis20180270
  • 病态观测方程广泛存在于大地测量解算模型中,如卫星定位[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)

      式中,LRm为观测向量;eRm为误差向量;ARm×n为系数矩阵;XRn为待估参数。其统计性质如下:

      $$ \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)

      式中,URm×m, VRn×n且满足UTU=UUT=Im, VTV=VVT=InImIn分别表示m×mn×n阶单位阵; ΣRm×n, $\Sigma = \left[ {\begin{array}{*{20}{c}} D\\ 0 \end{array}} \right]$, Dn×n阶对角阵,其对角元λi(1≤in)为A的奇异值,且按降序排列。若将UV用列向量表示成U=(u1 u2um), V=(v1 v2vn),则A可展开成:

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

      将式(7)代入式(5),可得XLSD(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项,截掉后nk项,得到新的系数矩阵:

      $$ {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);ΣkRm×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)

      式中,ΣTRn×n为对角阵;其对角元为diag($\frac{1}{{\lambda _1^2}}\frac{1}{{\lambda _2^2}} \cdots \frac{1}{{\lambda _k^2}}0 \cdots 0$);V2V的第k+1~n列向量构成的矩阵,即V2=(vk+1 vk+2vn)。

    • U1=(u1 u2uk), U2=(uk+1 uk+2um), 顾及式(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]。

      取采样间隔为Δxy=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个已知点P1P2P9,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%,显然优于传统公式所估的结果。

参考文献 (24)

目录

    /

    返回文章
    返回