-
在地球物理和大地测量领域的诸多场景中，数学模型呈现高度的病态性，导致参数的最小二乘估值变得很不稳定。解决病态问题的一般思路是对病态方程进行合适的正则化处理，以求得稳定且误差较小的有偏估计。最著名的正则化方法是Tikhonov正则化方法[1-2]，Tikhonov正则化方法的关键在于选取适当的正则化参数或正则化矩阵。通常，在测量之前针对待求参数缺乏先验信息，此时正则化矩阵一般被定为单位矩阵，而只调节正则化参数，这样的估计被称作岭估计[3]，如文献[4-5]提出的偏差修正的正则化算法，文献[6]中的广义岭估计直接解法也是岭估计的变形。岭估计中正则化参数一般采用L曲线算法确定[7]。针对不同应用领域的特定条件，学者们提出了多种其他Tikhonov正则化方案。如文献[8]结合卫星飞行轨道约束提出了加速度约束下的正则化算法。
GNSS快速精密定位领域中的模糊度解算问题是一个典型的病态问题，许多学者研究过正则化方法在此问题中的应用。比如文献[9-10]利用岭估计求解模糊度浮点解；文献[11]对基线和模糊度向量提出了双k岭估计；文献[12]提出了两步解算法；文献[13]利用基线先验值和先验精度信息选取合适的正则化参数。其中，文献[9]指出采用多个正则化参数和岭估计对比，效果没有明显差异。
本文结合贝叶斯(Bayes)估计理论，研究未知参数无先验信息条件下的Tikhonov正则化方法。Bayes估计理论的重要一环是先验分布的引入，文献[14]采用Bayes统计观点，引入未知参数的无信息先验分布研究卫星精密定位；文献[15]研究了正则化理论与Bayes估计之间的联系。本文引入未知参数的无信息先验分布，研究均方误差(mean square error, MSE)意义下最优的Tikhonov正则化矩阵；结合文献[16]中的参数估计方法提出了一种新的快速截断奇异值正则化方法，并探讨了该方法在GNSS模糊度解算中的应用。在实验中，通过一组GNSS观测数据比较了最小二乘(least squares，LS)方法、L曲线岭估计和本文所提方法的解算成功率和计算耗时。
-
测量系统观测方程经线性化和必要的矩阵变换后可转化为观测误差独立同分布的形式：
$$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{\varepsilon }} $$ (1) 式中，A为秩为n的m×n维设计矩阵；X为n×1待估参数向量；L为m×1观测向量；ε为m×1误差向量，$ E\left( {{\rm{ }}\mathit{\boldsymbol{\varepsilon }}{\rm{ }}} \right) = 0, E\left( {{\rm{ }}\mathit{\boldsymbol{\varepsilon \varepsilon }}{^{\rm{T}}}} \right) = \sigma _0^2\mathit{\boldsymbol{I}}$，其中E(·)为数学期望。Tikhonov正则化在最小二乘目标函数中加入合适的稳定性泛函改善参数估计的稳定性，其目标函数为：
$$ \min :{\left\| {\mathit{\boldsymbol{A\hat X}} - \mathit{\boldsymbol{L}}} \right\|^2} + \mathit{\Omega }\left( {\mathit{\boldsymbol{\hat X}}} \right) $$ (2) 式中，$\parallel \cdot {\parallel ^2} $是欧氏2范数；$\mathit{\Omega} \left( {\mathit{\boldsymbol{\hat X}}} \right) $称为稳定性泛函，$ \mathit{\Omega} \left( {\mathit{\boldsymbol{\hat X}}} \right) = {{\mathit{\boldsymbol{\hat X}}}^{\rm{T}}}\mathit{\boldsymbol{H\hat X}}{\rm{ }}, \mathit{\boldsymbol{H}}$是一个非负定对称矩阵，$ {\mathit{\boldsymbol{\hat X}}}$称为X正则化矩阵, 为X的估计值。一般情况下，在测量前待估计参数可以认为是完全未知的。因此，本文通过引入未知参数的无信息先验分布模型，求取均方误差最小意义下最优的Tikhonov正则化矩阵。
把待估参数看作随机变量是Bayes统计理论的一贯观点。由于X无先验信息，可认为待估计参数X的元素可以在其定义域中以等概率取任意值。对于元素为同一类型的未知参数向量，Bayes统计通常设定其先验分布模型为：X元素在其定义域中均匀分布，且不同元素之间相互独立[14-15]。数学模型为：
$$ \left. \begin{array}{l} f\left( {{x_i}} \right) \propto 1\\ f\left( {{x_i}{x_j}} \right) = f\left( {{x_i}} \right)f\left( {{x_j}} \right),i \ne j \end{array} \right\} $$ (3) 式中，xi和xj为X中的不同元素；f(·)表示概率密度函数，∝是正比于符号，f(xi)∝1表示xi在其取值范围内等概率分布。由式(3)可知X的元素独立同分布，假设元素方差为σx2，则其协方差矩阵为：
$$ E\left( {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}}} \right) = \sigma _x^2\mathit{\boldsymbol{I}} $$ (4) 借助该分布模型，本文可以证明如下命题。
命题1：对未知参数无先验信息的条件下，均方误差最小意义下最优的Tikhonov正则化矩阵R具有$ \mathit{\boldsymbol{R}}{\rm{ }} = \mathit{\boldsymbol{ V \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\mathit{\boldsymbol{R}}}}{\mathit{\boldsymbol{V}}^{\rm{T}}}$形式。其中$\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\mathit{\boldsymbol{R}}}} $为主对角矩阵，V为设计矩阵奇异值分解的右正交矩阵。
命题1等价于如下命题：对任意一个正则化矩阵H，若其对应Tikhonov正则解的均方误差为$\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) $，则必能找到另一正则化矩阵$\mathit{\boldsymbol{R}}{\rm{ }} = \mathit{\boldsymbol{ V \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\mathit{\boldsymbol{R}}}}{\mathit{\boldsymbol{V}}^{\rm{T}}} $，其对应Tikhonov正则解的均方误差为$ \ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}}} \right)$，使得$\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}}} \right) \le \ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) $。其中$\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\mathit{\boldsymbol{R}}}} $为主对角矩阵，V为设计矩阵奇异值分解的右正交矩阵。
证明：假设H为任一正则化矩阵，由式(1)、(2)可以得到其正则解：
$$ {{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} + \mathit{\boldsymbol{H}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\left( {\mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{\varepsilon }}} \right) $$ 对A进行奇异值分解$\mathit{\boldsymbol{A}} = \mathit{\boldsymbol{ U \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\sigma \mathit{\boldsymbol{}}}}{\mathit{\boldsymbol{V}}^{\rm{T}}} $，其中U和V均为酉矩阵。令$\mathit{\boldsymbol{K}}{\rm{ }} = {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{HV}} $ (因为H为对称矩阵，所以K也为对称矩阵；因为$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\sigma \mathit{\boldsymbol{}}}}$为对角矩阵，所以$\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma ^{\rm{T}}\mathit{\boldsymbol{ = \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\sigma \mathit{\boldsymbol{}}}} $。为了表达式的简洁，本文不对$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{}}_{\sigma \mathit{\boldsymbol{}}}}$与$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma ^{\rm{T}}$进行区分)。把${{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}} $用A的奇异值和矩阵K重新表达：
$$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}} = {{\left( {\mathit{\boldsymbol{V}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{V}}^{\rm{T}}} + \mathit{\boldsymbol{VK}}{\mathit{\boldsymbol{V}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{V}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma }{\mathit{\boldsymbol{U}}^{\rm{T}}}\left( {\mathit{\boldsymbol{U}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma }{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{X}} + } \right.}\\ {\left. \mathit{\boldsymbol{\varepsilon }} \right) = \mathit{\boldsymbol{V}}{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + } \right.}\\ {{{\left. \mathit{\boldsymbol{K}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma }{\mathit{\boldsymbol{U}}^{\rm{T}}}\mathit{\boldsymbol{\varepsilon }}} \end{array} $$ (5) 式中，Λσ2为主对角线上的值为σi2的n维对角矩阵。把式(5)代入均方误差表达式$ \ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = E[{\left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}} - {\rm{ }}\mathit{\boldsymbol{X}}{\rm{ }}} \right)^{\rm{T}}}({{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}} - {\rm{ }}\mathit{\boldsymbol{X}}{\rm{ }})]$，测量噪声ε服从均值为0的Gauss分布，所以交叉项为0，故：
$$ \begin{array}{*{20}{c}} {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = E\left\{ {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{V}}\left[ {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}} \right]\left[ {\mathit{\boldsymbol{I}} - } \right.} \right.}\\ {\left. {\left. {{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right]{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right\} + E\left[ {{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}\mathit{\boldsymbol{U}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma }\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + } \right.} \right.}\\ {\left. {{{\left. \mathit{\boldsymbol{K}} \right)}^{ - 1}}{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma }{\mathit{\boldsymbol{U}}^{\rm{T}}}\mathit{\boldsymbol{\varepsilon }}} \right]}\\ {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = {\rm{tr}}\left\{ {\left[ {\mathit{\boldsymbol{I}} - {{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right]{\mathit{\boldsymbol{V}}^{\rm{T}}}E\left( {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}}} \right) \cdot } \right.}\\ {\left. {\mathit{\boldsymbol{V}}\left[ {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}} \right]} \right\} + {\rm{tr}}\left[ {{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma } \cdot } \right.}\\ {\left. {{\mathit{\boldsymbol{U}}^{\rm{T}}}E\left( {\mathit{\boldsymbol{\varepsilon }}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}} \right)\mathit{\boldsymbol{U}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\sigma }{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}} \right]} \end{array} $$ (6) 式中，tr(·)表示矩阵的迹。由上文可知$E\left( {{\rm{ }}\mathit{\boldsymbol{\varepsilon }}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}} \right) = \sigma _0^2\mathit{\boldsymbol{I}}, E\left( {{\rm{ }}\mathit{\boldsymbol{XX}}{^{\rm{T}}}} \right) = \sigma _x^2\mathit{\boldsymbol{I}} $，代入式(6)得：
$$ \begin{array}{*{20}{c}} {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = \sigma _x^2 \cdot {\rm{tr}}\left\{ {\left[ {\mathit{\boldsymbol{I}} - {{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right]\left[ {\mathit{\boldsymbol{I}} - } \right.} \right.}\\ {\left. {\left. {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}} \right]} \right\} + \sigma _0^2 \cdot {\rm{tr}}\left[ {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + } \right.} \right.}\\ {\left. {{{\left. \mathit{\boldsymbol{K}} \right)}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} + \mathit{\boldsymbol{K}}} \right)}^{ - 1}}} \right]} \end{array} $$ (7) 设$ {\left( {{\rm{ }}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}} + {\rm{ }}\mathit{\boldsymbol{K}}{\rm{ }}} \right)^{ - 1}} = {\rm{ }}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + {\rm{ }}\mathit{\boldsymbol{P}}$，其中$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}$是由$ {\left( {{\rm{ }}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}} + {\rm{ }}\mathit{\boldsymbol{K}}{\rm{ }}} \right)^{ - 1}}$主对角线上的元素组成的主对角矩阵；P是由$ {\left( {{\rm{ }}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}} + {\rm{ }}\mathit{\boldsymbol{K}}{\rm{ }}} \right)^{ - 1}}$非主对角线上的元素组成的矩阵。则$ \mathit{\boldsymbol{P}}{\rm{ }} = {\left( {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}} + {\rm{ }}\mathit{\boldsymbol{K}}{\rm{ }}} \right)^{ - 1}} - {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}$，且P的主对角线上的元素皆为0。式(7)可变为：
$$ \begin{array}{*{20}{c}} {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = \sigma _x^2 \cdot {\rm{tr}}\left\{ {\left[ {\mathit{\boldsymbol{I}} - \left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + \mathit{\boldsymbol{P}}} \right){\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right]\left[ {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} \cdot } \right.} \right.}\\ {\left. {\left. {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + \mathit{\boldsymbol{P}}} \right)} \right]} \right\} + \sigma _0^2 \cdot {\rm{tr}}\left[ {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + \mathit{\boldsymbol{P}}} \right){\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + \mathit{\boldsymbol{P}}} \right)} \right]} \end{array} $$ (8) 因为${\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}}} $和${\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} $都是对角矩阵，所以它们也都是对称矩阵，考虑到K为对称矩阵，所以P也是对称矩阵。故将式(8)展开可以得到：
$$ \begin{array}{l} \ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = \sigma _x^2 \cdot {\rm{tr}}\left[ {\mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^4}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + 2\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^4}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + } \right.\\ \;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^4}}}{\mathit{\boldsymbol{P}}^{\rm{T}}} - 2{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}} - 2\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right] + \sigma _0^2 \cdot \\ \;\;\;\;\;\;\;{\rm{tr}}\left[ {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + \mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{P}}^{\rm{T}}} + 2\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}} \right] \end{array} $$ (9) 由于矩阵P主对角线上的元素为0，则$\mathit{\boldsymbol{P \boldsymbol{\varLambda} }}{_{{\sigma ^4}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}、\mathit{\boldsymbol{P \boldsymbol{\varLambda} }}{_{{\sigma ^2}}}、{\rm{ }}\mathit{\boldsymbol{P \boldsymbol{\varLambda} }}{_{{\sigma ^2}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} $的主对角线元素均为0。式(9)可以化简为：
$$ \begin{array}{*{20}{c}} {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) = \sigma _x^2 \cdot {\rm{tr}}\left[ {\mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^4}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + 2\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^4}}}{\mathit{\boldsymbol{P}}^{\rm{T}}} - } \right.}\\ {\left. {2{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right] + \sigma _0^2 \cdot {\rm{tr}}\left[ {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} + \mathit{\boldsymbol{P}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{P}}^{\rm{T}}}} \right]} \end{array} $$ (10) 式中，${\rm{tr}}\left( {{\rm{ }}\mathit{\boldsymbol{P \boldsymbol{\varLambda} }}{_{{\sigma ^4}}}\mathit{\boldsymbol{P}}{^{\rm{T}}}} \right) = \sum\limits_{i = 1}^n {} \sum\limits_{j = 1}^n {} {\rm{ }}p_{ij}^2\sigma _j^4 \ge 0, {\rm{tr}}({\rm{ }}\mathit{\boldsymbol{P \boldsymbol{\varLambda} }}{_{{\sigma ^2}}}\mathit{\boldsymbol{P}}{^{\rm{T}}}) = \sum\limits_{i = 1}^n {} \sum\limits_{j = 1}^n {} {\rm{ }}p_{ij}^2\sigma _j^2 \ge 0 $；其中pij是矩阵P中的第i行第j列元素。所以有如下不等式成立：
$$ \begin{array}{*{20}{c}} {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) \ge \sigma _x^2 \cdot {\rm{tr}}\left[ {\mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^4}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}} - 2{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}} \right] + }\\ {\sigma _0^2 \cdot {\rm{tr}}\left[ {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{\sigma ^2}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}} \right]} \end{array} $$ (11) 当且仅当P为零矩阵时式(11)等号成立。由矩阵P的定义可知，当P为零矩阵时，即$ {\left( {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}} + {\rm{ }}\mathit{\boldsymbol{K}}{\rm{ }}} \right)^{ - 1}} = {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}$，因为$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}}$为主对角矩阵，则$\mathit{\boldsymbol{K}}{\rm{ }} = \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{P}}^{ - 1} - \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{\sigma ^2}}} $也是一个主对角矩阵。若记K为$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{R}}}$，则由矩阵K的定义知，此时正则化矩阵为$ \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{R}}}{\mathit{\boldsymbol{V}}^{\rm{T}}}$，记为R。式(11)不等号右侧即为当正则化矩阵取$\mathit{\boldsymbol{R}}= \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{R}}}{\mathit{\boldsymbol{V}}^{\rm{T}}}$时的均方误差，即有：
$$ \ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{H}}}} \right) \ge \ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}}} \right) $$ (12) 命题1得证。下面根据命题1，进一步求解最优正则化矩阵，即求解其主对角矩阵$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{R}}}$。设$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_\mathit{\boldsymbol{R}}}$的元素为ri(i=1, 2…n)，Tikhonov正则化解为$ {{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}}\mathit{\boldsymbol{}} = {\left( {{\rm{ }}\mathit{\boldsymbol{A}}{^{\rm{T}}}\mathit{\boldsymbol{A}}{\rm{ }} + \mathit{\boldsymbol{R}}{\rm{ }}} \right)^{ - 1}}\mathit{\boldsymbol{A}}{^{\rm{T}}}\left( {{\rm{ }}\mathit{\boldsymbol{AX}}{\rm{ }} + {\rm{ }}\mathit{\boldsymbol{\varepsilon }}{\rm{ }}} \right)$，对A进行奇异值分解$\mathit{\boldsymbol{A}} = {\rm{ }}\mathit{\boldsymbol{U \boldsymbol{\varLambda} }}{_\sigma }{\mathit{\boldsymbol{V}}^{\rm{T}}}$，设ui和vi分别为酉矩阵U和V的第i个列向量，由该奇异值分解把${{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}} $的表达式展开可得：
$$ {{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}} = \sum\limits_{i = 1}^n {\left( {\frac{{\sigma _i^2}}{{\sigma _i^2 + {r_i}}}{\mathit{\boldsymbol{v}}_i}\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}} + \frac{{{\sigma _i}}}{{\sigma _i^2 + {r_i}}}{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}\varepsilon } \right)} $$ (13) 由式(13)和均方误差的定义可得：
$$ \begin{array}{*{20}{c}} {\ell \left( {{{\mathit{\boldsymbol{\hat X}}}_\mathit{\boldsymbol{R}}}} \right) = E\left[ {\sum\limits_{i = 1}^n {{{\left( {\frac{{{r_i}}}{{\sigma _i^2 + {r_i}}}} \right)}^2}{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{v}}_i}\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}} + } \right.}\\ {\left. {{{\left( {\frac{{{\sigma _i}}}{{\sigma _i^2 + {r_i}}}} \right)}^2}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}\mathit{\boldsymbol{\varepsilon }}} \right]} \end{array} $$ (14) 式(14)右边每一项对ri求导，令导数等于0，可解得其唯一驻点：
$$ {r_i} = \frac{{E\left( {\mathit{\boldsymbol{u}}_i^{\rm{T}}\mathit{\boldsymbol{\varepsilon }}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}{\mathit{\boldsymbol{u}}_i}} \right)}}{{E\left( {\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{v}}_i}} \right)}} $$ (15) 式(15)函数的导数在其驻点的左邻域小于零，右邻域大于零，故该驻点就是均方误差取最小值时的岭参数。
-
§1虽然求得了最优Tikhonov正则化矩阵表达式，但在式(15)中，ε和X都是未知量，必须对其进行合理变形和估计。再次使用X的无信息分布模型，把$ E\left( {{\rm{ }}\mathit{\boldsymbol{\varepsilon }}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}} \right) = \sigma _0^2\mathit{\boldsymbol{I}}, E\left( {{\rm{ }}\mathit{\boldsymbol{XX}}{^{\rm{T}}}} \right) = \sigma _x^2\mathit{\boldsymbol{I}}$代入式(15)得:
$$ {r_i} = \frac{{\sigma _0^2}}{{\sigma _x^2}} $$ (16) 此时正则化矩阵简化为正则化参数与单位阵的乘积，Tikhonov正则化简化为岭估计。σ02是测量噪声方差，在病态问题中，最小二乘法对测量误差的估计受设计矩阵的条件数影响不大，仍然可以用最小二乘残差来估计测量误差：
$$ \hat \sigma _{0{\rm{LS}}}^2 = \frac{{{{\left( {\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - \mathit{\boldsymbol{L}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - \mathit{\boldsymbol{L}}} \right)}}{{m - n}} $$ (17) 同时，在病态问题中，最小二乘法对X估计的准确性太低，不适合用${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} $估计σx2。文献[16]中提出了一种TSVD正则化截断参数选取方法，可以求得比较准确的正则解。文献[16]通过选取局部最优截断参数和全局最优截断参数两个步骤确定TSVD正则化截断参数。满足式(18)的参数k被认为是局部最优的截断参数：
$$ \left. \begin{array}{l} \frac{{2{{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}\mathit{\boldsymbol{\hat \varepsilon }}}}{{\chi _{1 - \frac{\alpha }{2}}^2\left( {m - n} \right)}} < {\left( {\mathit{\boldsymbol{u}}_k^{\rm{T}}\mathit{\boldsymbol{L}}} \right)^2}\\ \frac{{2{{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}\mathit{\boldsymbol{\hat \varepsilon }}}}{{\chi _{\frac{\alpha }{2}}^2\left( {m - n} \right)}} > {\left( {\mathit{\boldsymbol{u}}_{k + 1}^{\rm{T}}\mathit{\boldsymbol{L}}} \right)^2} \end{array} \right\} $$ (18) 式中，$ {\mathit{\boldsymbol{\hat \varepsilon }}}$是误差向量的估计值，可选用误差的最小二乘估计；α为$ {\chi ^2}(m - n)$分布函数的分位数，选为0.5。在局部最优截断参数集合里选取任意两个，通过Δ值比较选取全局最优截断参数:
$$ \Delta = \frac{{\mathit{\boldsymbol{\hat X}}_{{k_p}}^{\rm{T}}{{\mathit{\boldsymbol{\hat X}}}_{{k_p}}}\left( {{k_q} - {k_p}} \right)}}{n} - \sum\limits_{i = {k_p} + 1}^{{k_q}} {\frac{{\hat \sigma _{{\rm{0LS}}}^2}}{{\sigma _i^2}}} $$ (19) 式中，kq和kp为任意两个局部最优的截断参数；$ {{\mathit{\boldsymbol{\hat X}}}_{{k_p}}}$是以kp为截断参数的正则化解。若Δ≤0，则认为kp是均方误差意义下比kq更优的截断参数，反之亦然。选取完全局最优的截断参数后，即可求得最终的正则解$ {{\mathit{\boldsymbol{\hat X}}}_t}$ (限于篇幅，该算法推导细节本文不展开介绍，有兴趣的读者可以参考原文献)。本文利用${{\mathit{\boldsymbol{\hat X}}}_t} $估计σx2:
$$ \hat \sigma _{xt}^2 = \frac{{\mathit{\boldsymbol{\hat X}}_t^{\rm{T}}{{\mathit{\boldsymbol{\hat X}}}_t}}}{n} $$ (20) 再把$\hat \sigma _{0{\rm{LS}}}^2$与$\hat \sigma _{xt}^2{{\mathit{\boldsymbol{\hat X}}}_t} $代入式(16)中求得正则化参数的估计值。
-
GNSS精密定位观测模型是一个混合整数模型：
$$ \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{A}}_a}\mathit{\boldsymbol{a}} + \mathit{\boldsymbol{Bb}} + \mathit{\boldsymbol{e}} $$ (21) 式中，y 是载波相位与伪码观测向量；a是模糊度向量；b是实参数向量；Aa和B是联系未知参数和观测量之间的设计矩阵；e是误差向量，其协方差矩阵为Qyy。由于向量a和b是不同类型的参数，它们的先验分布并不一致，本文通过投影变换去除b，只对模糊度向量a进行正则化估计。令${\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}} = {\rm{ }}\mathit{\boldsymbol{B}}{\rm{ }}{\left( {{\rm{ }}\mathit{\boldsymbol{B}}{{\rm{}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{yy}^{ - 1}\mathit{\boldsymbol{B}}{\rm{ }}} \right)^{ - 1}}\mathit{\boldsymbol{B}}{{\rm{}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{yy}^{ - 1}, {\rm{ }}\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot = {\rm{ }}\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}} $，式(21)两边同乘$\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot $得：
$$ \mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot {\mathit{\boldsymbol{A}}_a}\mathit{\boldsymbol{a}} + \mathit{\boldsymbol{e}} $$ (22) 对Qyy进行满秩分解：$ \mathit{\boldsymbol{Q}}{_{yy}} = {\rm{ }}\mathit{\boldsymbol{D}}{^{\rm{T}}}\mathit{\boldsymbol{D}}$，然后对式(22)继续做变换：
$$ {\mathit{\boldsymbol{D}}^{ - {\rm{T}}}}\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{D}}^{ - {\rm{T}}}}\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot {\mathit{\boldsymbol{A}}_a}\mathit{\boldsymbol{a}} + {\mathit{\boldsymbol{D}}^{ - {\rm{T}}}}\mathit{\boldsymbol{e}} $$ (23) 令$\mathit{\boldsymbol{D}}{^{{\rm{ - T}}}}\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot \mathit{\boldsymbol{y}}{\rm{ }} = {\rm{ }}\mathit{\boldsymbol{L}}{\rm{ }}{\rm{ }}\mathit{\boldsymbol{D}}{^{{\rm{ - T}}}}\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}^ \bot \mathit{\boldsymbol{A}}{_a} = \mathit{\boldsymbol{A}}, {\rm{ }}\mathit{\boldsymbol{D}}{^{{\rm{ - T}}}}\mathit{\boldsymbol{e}}{\rm{ }} = {\rm{ }}\mathit{\boldsymbol{\varepsilon }} $，则上式成为如式(1)的标准观测模型。可直接对式(23)使用本文所提出的Tikhonov正则化方法求得模糊度浮点解${{\mathit{\boldsymbol{\hat a}}}_r} $和协方差矩阵$ {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{a}}_r}}}$。进而确定求取模糊度固定解的目标函数：
$$ \min :{\left( {\mathit{\boldsymbol{a}} - {{\mathit{\boldsymbol{\hat a}}}_r}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{a}}_r}}^{ - 1}\left( {\mathit{\boldsymbol{a}} - {{\mathit{\boldsymbol{\hat a}}}_r}} \right) $$ (24) 式(24)可以用最小二乘模糊度降相关(least-squares ambiguity decorrelation adjustment，LAMBDA)方法求最优解，也可以用Bootstraping等次优方法求取次优解[16]。
-
实验基于美国全球定位系统(Global Positioning System, GPS)连续运行参考站网的一组实测数据，用MC01和MC02两个相距9 km的站点于世界协调时2015年5月1日24 h观测数据，观测量为L1载波相位与P1伪码，采样间隔30 s，卫星的截止高度角设为15°，随机模型采用等权模型，伪码和相位观测噪声分别设为0.3 m和0.003 m。为了提高可靠性，本文实验设定共视卫星数大于或等于8的历元进行模糊度解算，共1 267个有效历元。
需要指出的是，在日常GNSS精密定位解算中，所用随机模型都是简化的，并不能100%正确反映观测噪声的真实情况。正则化算法对随机模型的准确性要求较高，为了排除随机模型不准确所造成的干扰，本文实验采用Monte Carlo仿真观测噪声。这种做法保留了真实的卫星结构和设计矩阵，同时避免了不准确的随机模型造成的干扰，被很多文献采用[9, 17-18]。分别用最小二乘估计、L曲线岭估计和本文提出的新正则化估计，逐个历元解算模糊度浮点解，并进一步求取固定解，比较不同方法解算的成功率和计算耗时。
实验中两站点的共视卫星数和位置精度因子(position dilution of precision，PDOP)值的变换情况如图 1所示。分别用3种方法求得模糊度浮点解，图 2画出了各历元浮点解与真值之间的均方误差，可见LS方法所得解均方误差明显大于其他两种方法；而表 1则对所有历元的浮点解求取了平均均方误差，其显示本文方法平均均方差略小于L曲线岭估计。用这些浮点解的协方差矩阵，采用Monte Carlo仿真方法计算出每个历元模糊度解算的理论成功率[9]。图 3画出了所有解算历元成功率值的概率累积分布统计图。由图 3可知，用LS方法计算的成功率较小的历元所占比例远大于其他两种方法，而L曲线法又略大于本文方法。图 4是与图 3对应的概率密度分布统计条形图，从图 4中也可以看出两种正则化方法在解算成功率上差异不大，但本文方法中成功率较高(如大于90%)的历元数比例稍高。图 5画出了实际解得的模糊度固定解均方差情况，其中模糊度正确固定时均方差为0，可以看出两种正则化方法模糊度正确固定的数目远高于LS方法。表 1列出了成功固定数目和平均解算成功率情况，亦可见本文方法成功固定历元数稍高于L曲线岭估计。图 6画出了3种方法的计算耗时，表 1也列出了平均耗时，虽然3种方法实时性都比较高，但是L曲线岭估计的计算耗时高于其余两种方法一个数量级以上。这在需要处理大量历元数据的情况下其劣势将显得比较明显。
图 1 各历元的可见卫星数和PDOP值
Figure 1. PDOP and Visible Satellites for Each Epoch
图 2 3种方法所求得的模糊度浮点解均方误差
Figure 2. MSE of Float Solutions Solved by 3 Methods
表 1 3种方法解算性能平均值统计
Table 1. Average Computational Properties of 3 Methods
算法 平均浮点解均方差/周 成功固定历元数/个 平均解算成功率/% 平均耗时/ms LS估计 9.89 497 39.2 1.06 L曲线岭估计 1.70 971 76.6 41.08 本文方法 1.67 1 003 79.2 4.15 
图 3 3种方法求得成功率的概率累积分布统计
Figure 3. Cumulative Distribution of SR Solved by 3 Methods
图 4 3种方法求得成功率的概率密度分布统计
Figure 4. Density Distribution of SR Solved by 3 Methods
图 5 3种方法所求得的模糊度固定解均方误差
Figure 5. MSE of the Fixed Solutions Solved by 3 Methods
图 6 3种方法计算耗时
Figure 6. Time Consumption of 3 Methods
Optimal Tikhonov Regularization Matrix and Its Application in GNSS Ambiguity Resolution
-
摘要: 首先，用贝叶斯（Bayes）统计理论的观点，把未知参数看作随机变量，引入未知参数的无信息先验分布函数，从数学上推导了均方误差最小意义下的正则化矩阵；然后，结合最优正则化矩阵和快速截断奇异值算法，提出了一种新的正则化方法；最后，探讨了新方法在全球卫星导航系统（Global Navigation Satellite System，GNSS）模糊度解算中的应用。通过一组GNSS模糊度解算实验，比较了最小二乘（least squares，LS）方法、L曲线岭估计和新方法的性能。结果表明，新方法解算成功率略高于L曲线岭估计，远高于LS方法；计算耗时略大于LS方法，远小于L曲线岭估计。
-
关键词:
- 病态问题 /
- Tikhonov正则化 /
- Bayes统计理论 /
- GNSS /
- 模糊度解算
Abstract: This contribution can be mainly divided into 3 aspects:(1) Based on Bayesian theory, unknown parameters are treated as random varies and their non-informative prior distribution function is introduced. Mathematical analysis is carried out to drive the optimal Tikhonov regularization matrix in the sense of minimizing the mean square error (MSE) of the solutions. (2) Combining the efficient truncated singular value decomposition (eTSVD), a new regularization method is proposed. (3) Global Navigation Satellite System(GNSS) ambiguity resolution application of the new method is discussed. Least squares (LS) estimation, ridge estimation of L curve and the new algorithm are compared by a group of GNSS ambiguity resolution experiments. The results show that the MSE of the new algorithm is slightly smaller than ridge estimation of L curve and much smaller than LS, however, the computational cost of the new algorithm is slightly more than LS but much less than ridge estimation of L curve.-
Key words:
- ill-posed problem /
- Tikhonov regularization /
- Bayesian statistical theory /
- GNSS /
- ambiguity resolution
-
表 1 3种方法解算性能平均值统计
Table 1. Average Computational Properties of 3 Methods
算法 平均浮点解均方差/周 成功固定历元数/个 平均解算成功率/% 平均耗时/ms LS估计 9.89 497 39.2 1.06 L曲线岭估计 1.70 971 76.6 41.08 本文方法 1.67 1 003 79.2 4.15 
下载: 导出CSV
-
[1] Tikhonov A N. Regularization of Ill-Posed Problems[J]. Dokl Akad Nauk SSSR, 1963, 151(1):49-52 http://d.old.wanfangdata.com.cn/OAPaper/oai_doaj-articles_eae560fae48c249b50ae0ce63aa285d9 [2] Tikhonov A N. Solution of Incorrectly Formulated Problems and the Regularization Method[J]. Dokl Akad Nauk SSSR, 1963, 151(3):501-504 [3] Hoerl A E, Kennard R W. Ridge Regression:Bia-sed Estimation for Non-orthogonal Problems[J]. Technometrics, 1970, 12:55-67 doi: 10.1080/00401706.1970.10488634 [4] Xu P L, Rummel R. A Simulation Study of Smoothness Methods in Recovery of Regional Gra-vity Fields[J]. Geophys J Int, 1994, 117:472-486 doi: 10.1111/gji.1994.117.issue-2 [5] Shen Y Z, Xu P L, Li B F. Bias-Corrected Regulari-zed Solution to Inverse Ill-Posed Models[J]. J Geod, 2012, 86:597-608 doi: 10.1007/s00190-012-0542-y [6] 游扬声, 王新洲, 刘星.广义岭估计的直接解法[J].武汉大学学报·信息科学版, 2002, 27(2):175-178 http://ch.whu.edu.cn/CN/abstract/abstract4940.shtml You Yangsheng, Wang Xinzhou, Liu Xing. Direct Solution to Generalized Ridge Estimate[J]. Geomatics and Information Science of Wuhan University, 2002, 27(2):175-178 http://ch.whu.edu.cn/CN/abstract/abstract4940.shtml [7] Hansen P C. Analysis of Discrete Ill-posed Problems by Means of the L-curve[J]. SIAM Review, 1992(1):561-580 http://www.emeraldinsight.com/servlet/linkout?suffix=b4&dbid=16&doi=10.1108%2F03321640110383933&key=10.1137%2F1034115 [8] 柳响林, Pavel Ditmar.基于B-spline和正则化算法的低轨卫星轨道平滑[J].地球物理学报, 2006, 49(1):99-105 doi: 10.3321/j.issn:0001-5733.2006.01.014 Liu Xianglin, Pavel D. Smoothing a Satellite Orbit on the Basis of B-spline and Regularization[J]. Chinese Journal of Geophysics, 2006, 49(1):99-105 doi: 10.3321/j.issn:0001-5733.2006.01.014 [9] 郭海涛, 张保明, 归庆明.广义岭估计在解算单线阵CCD卫星影像外方位元素中的应用[J].武汉大学学报·信息科学版, 2003, 28(4):444-447 http://ch.whu.edu.cn/CN/abstract/abstract4868.shtml Guo Haitao, Zhang Baoming, Gui Qingming. Application of Generalized Ridge Estimate to Computing the Exterior Orientation Elements of Satellite Linear Array Scanner Imagery[J]. Geomatics and Information Science of Wuhan University, 2003, 28(4):444-447 http://ch.whu.edu.cn/CN/abstract/abstract4868.shtml [10] Li B F, Shen Y Z, Feng Y M. Fast GNSS Ambi-guity Resolution as an Ill-Posed Problem[J]. J Geod, 2010, 84:683-698 doi: 10.1007/s00190-010-0403-5 [11] Shen Y Z, Li B F. Regularized Solution to Fast GPS Ambiguity Resolution[J]. J Surv Eng, 2007, 133:168-172 doi: 10.1061/(ASCE)0733-9453(2007)133:4(168) [12] Gui Q M, Han S H. New Algorithm of GPS Rapid Positioning Based on Double-k-Type Ridge Estimation[J]. J Surv Eng, 2007, 133:67-72 [13] 王振杰, 欧吉坤, 柳林涛.一种解算病态问题的方法:两步解法[J].武汉大学学报·信息科学版, 2005, 30(9):821-824 http://ch.whu.edu.cn/CN/abstract/abstract2283.shtml Wang Zhenjie, Ou Jikun, Liu Lintao. A Method for Resolving Ill-Conditioned Problems:Two-Step Solution[J]. Geomatics and Information Science of Wuhan University, 2005, 30(9):821-824 http://ch.whu.edu.cn/CN/abstract/abstract2283.shtml [14] Zhu J, Ding X, Chen Y. Maximum-Likelihood Ambiguity Resolution Based on Bayesian Principle[J]. J Geod, 2001, 75:175-187 doi: 10.1007/s001900100167 [15] Wu Z M, Bian S F. GNSS Integer Ambiguity Validation Based on Posterior Probability[J]. J Geod, 2015, 89:961-977 doi: 10.1007/s00190-015-0826-0 [16] 陈伟, 王子亭.传统正则化方法和Bayes统计理论之间的联系[J].河北理工大学学报(自然科学版), 2011, 33(3):104-106 doi: 10.3969/j.issn.1674-0262.2011.03.024 Chen Wei, Wang Ziting. The Relation Between the Traditional Regularization and the Bayesian Method[J]. Journal of Hebei Polytechnic University (Natural Science Edition), 2011, 33(3):104-106 doi: 10.3969/j.issn.1674-0262.2011.03.024 [17] Teunisen P J G. The Success Rate and Precision of GPS Ambiguities[J]. J Geod, 2000, 74(3/4):321-326 doi: 10.1007-s001900050289/ [18] Verhagen S, Teunissen P J G. The Ratio Test for Future GNSS Ambiguity Resolution[J]. GPS Solut, 2013, 17(4):535-548 doi: 10.1007/s10291-012-0299-z -
点击查看大图
图(6) / 表(1)
计量
- 文章访问数: 1308
- HTML全文浏览量: 124
- PDF下载量: 278
- 被引次数: 0
