岭估计在二次曲面高程拟合模型数据处理中的应用

The Application of Ridge Estimation in Data Processing of Quadric Surface Elevation Fitting Model

摘要: 在测量工作中，由于参数过多、观测结构不合理等各种原因，会导致各参数之间存在近似线性关系，此时经典的最小二乘模型解算结果会失真。本文通过二次曲面高程拟合模型着重讨论岭估计在最小二乘模型失效的情况下解算结果的可靠性，总结岭参数的求解方法，并验证其解算的可靠性。

Abstract: In the measurement work, due to too many parameters, unreasonable observation structure and other reasons, there will be an approximately linear relationship between the parameters. At this time, the solution results of the classical least squares model will be distorted. In this paper, through the elevation fitting quadric surface model, the reliability of the solution result of ridge estimation in the case of the failure of the least square model is discussed, the solution method of ridge parameters is summarized, and the reliability of its solution is verified.

1. 引言

通过GPS测量可以得到待测点的大地高，如果测区地势平坦、范围较小，可以采用多项式曲面拟合模型来拟合待测点正常高，其中较为常见的为二次曲面模型。但在实际生产中，由于参数的选择、观测条件的限制、结构设计的不合理等会导致参数之间有近似共线性，这使得法方程矩阵 $N=B^{\text{T}}B$ 接近奇异，此时我们称法方程为病态方程，为解决法方程病态时的最小二乘估计问题，有许多学者进行过探索，如基于Neumann级数的有偏估计方法、靶向改变法方程矩阵和基于LIU估计的迭代新方法等 [1] [2] [3]。

但是在众多有偏估计之中，岭估计、广义岭估计、主成分估计的影响最大。在岭估计的探索当中，学者对岭参数的选择方法不断创新，使得岭估计在实际应用中硕果累累，尤其是在周跳探测、高频GNSS钟差、GPS数据处理、高程拟合等领域获得了较高的可靠性 [4] [5] [6] [7] [8]。

2. 岭估计

2.1. 岭估计的原理

岭估计(Ridge Estimation)是Hoerl和Kennard与1970年提出来的，是目前最有影响力的有偏估计方法。其模型为：

$L=BX+\Delta ,\Delta ~N\left(0~{\sigma }^{2}I\right)$ (1)

其对应的解为：

$\hat{X}\left(k\right)={\left(N+kI\right)}^{-1}{B}^{\text{T}}PL$ (2)

式中k为大于零的任意常数，称为岭参数，而 $N=B^{\text{T}}PB$。可以看出，当k的选择不同，就会得到不同的岭估计，当 $k=0$ 时，为最小二乘估计。

2.2. 岭估计的性质

为了讨论岭估计的性质，将(1)式改写为以下形式：

$L=BX+\Delta =AY+\Delta$ (3)

其中： $A=BG$， $Y=G^{\text{T}}X$， $G=\left(G_1,G_2,\cdots ,G_t\right)$。

$G_i\left(i=1,2,\cdots ,t\right)$ 为 $B^{\text{T}}B$ 对应于特征值 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_t$ 的标准正交化特征向量，且满足

$A^{\text{T}}A={\left(BG\right)}^{\text{T}}\left(BG\right)=G^{\text{T}}{B}^{\text{T}}BG=\Lambda =\text{diag}\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_t\right)$ (4)

式(4)称之为Gauss-Markov模型的典则形式，Y称为典则参数。

则典则参数Y的最小二乘估计和岭估计分别为：

$\begin{array}{l}{\hat{Y}}_{\text{LS}}={\left(A^{\text{T}}A\right)}^{-1}{A}^{\text{T}}L={\Lambda }^{-1}{G}^{\text{T}}{B}^{\text{T}}L\\ \hat{Y}\left(k\right)={\left(A^{\text{T}}A+kI\right)}^{-1}{A}^{\text{T}}L={\left(\Lambda +kI\right)}^{-1}{G}^{\text{T}}{B}^{\text{T}}L\end{array}$ (5)

2.3. 岭参数的选择

在岭参数的选择上，国内学者进行了相关研究，较为常见的有U曲线法、行列式法和正则化法，还有部分学者将以上方法进行比较，确定最终岭参数 [9] [10] [11] [12]。

引进岭估计的目的是减少均方误差，统计学家己经证明了当 $k>0$ 时，

$\text{MSE}\left(\hat{X}\left(k\right)\right)<\text{MSE}\left({\hat{X}}_{\text{LS}}\right)$ (6)

因此在均方误差意义下，岭估计可改进最小二乘估计。对于岭参数k的选择，一直困扰着学者们，由于岭参数的选择在于使均方误差更小，即

$\text{MSE}\left(\hat{X}\left(k\right)\right)={\sigma }_0^2{\displaystyle \underset{i=1}{\overset{t}{\sum }}\frac{{\lambda }_i}{{\left({\lambda }_i+k\right)}^2}}+k^2{\displaystyle \underset{i=1}{\overset{t}{\sum }}\frac{y_i^2}{{\left({\lambda }_i+k\right)}^2}}$ (7)

最小，而Y，X未知，故不能用求极值的方法求取岭参数的值。国内外学者们经过不断探索，提出了许多岭参数的选择和求解的方法，如岭迹法、双h公式法。

1) 岭迹法：

就是以岭估计 ${\hat{X}}_i$ 的分量 ${\hat{X}}_i\left(k\right)\left(i=1,2,\cdots ,t\right)$ 作为岭参数k的函数，将t条岭迹画函数图像，选择使t条岭迹都处于较为稳定状态下的那个k值作为岭参数。

2) 双h公式法：

$k=\frac{h_1{\hat{\sigma }}_0^2}{{\hat{X}}_{\text{LS}}^{\text{T}}G{\hat{X}}_{\text{LS}}+h_2{\hat{\sigma }}_0^2}$ (8)

当取 $G=I$， $h_1=t$， $h_2=0$ 时，上式就变为了Hoerl-Kennard-Baldwin公式：

$k=\frac{t{\hat{\sigma }}_0^2}{{\hat{X}}_{\text{LS}}^{\text{T}}{\hat{X}}_{\text{LS}}}$ (9)

2.4. 二次曲面模型

通过GPS测量获得的高程属于大地高，已知大地高和正常高的关系为：

$H=H_r+\xi$ (10)

其中：H代表大地高， $H_r$ 代表正常高， $\xi$ 代表高程异常。

二次曲面模型可以表示为：

$\xi \left(x,y\right)=A_1+A_{2}x+A_{3}y+A_{4}xy+A_5{x}^2+A_6{y}^2$ (11)

由式(10)可知：

$\xi \left(i\right)=H_i-H_{ri},i=1,2,3,4,5,\cdots ,n$ (12)

其中： $H_i$ 代表第i个点的大地高， $H_{ri}$ 代表第i个点的正常高。

根据式(11)和式(12)，可列出误差方程：

$V=BA-\xi$ (13)

式(13)中各项矩阵的表达式为：

$\xi =\left[\begin{array}{c}{\xi }_1\\ {\xi }_2\\ ⋮\\ {\xi }_n\end{array}\right]$， $V=\left[\begin{array}{c}{V}_1\\ V_2\\ ⋮\\ V_n\end{array}\right]$， $A=\left[\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6\end{array}\right]$， $B=\left[\begin{array}{cccccc}1& x_1& y_1& x_1{y}_1& x_1^2& y_1^2\\ 1& x_2& y_2& x_2{y}_2& x_2^2& y_2^2\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& x_n& y_n& x_n{y}_n& x_n^2& y_n^2\end{array}\right]$ (14)

3. 岭估计的算例分析

本算例取自某工程实测数据，已知12个控制点，经GNSS静态观测获得这12个点的大地高，通过四等三角高程测量获得这12个点的正常高。这12个控制点的坐标及高程如表1所示。

取H2、H4、H5、H16、H17、H45、H46和H49共8个点进行模型参数计算点，H15、H18、H19和H21共4个点作为检核点，已知模型参数为6个，则多余观测量为2个。

$B^{\text{T}}B$ 的6个特征值为：

${\lambda }_1=4.8177\times {10}^{-9},{\lambda }_2=4.2623\times {10}^3,{\lambda }_3=1.8878\times {10}^9,{\lambda }_4=3.1617\times {10}^{17}$

${\lambda }_5=1.3791\times {10}^{22},{\lambda }_6=6.1291\times {10}^{26}$

Table 1. Coordinates of known point

表1. 已知点坐标

按照条件数的计算公式可知，条件数为：

$p=\frac{{\lambda }_{\max }}{{\lambda }_{\min }}=\frac{6.1291\times {10}^{26}}{4.8177\times {10}^{-9}}=1.2722\times {10}^{35}$

由于该算例中条件数远大于1000，可以认为该方程组为病态方程组 [13]，具有严重的复共线性。

3.1. 岭估计的解算

本文采用岭迹法和双h公式法对本算例进行解算：

1) 岭迹法

利用Matlab编程计算最小二乘拟合结果和岭估计拟合结果，将岭参数的范围控制在0~2内，步长设置为0.01，得到的岭迹图如图1(a)所示。由岭迹图可以看出，当k仅有极小扰动时，各参数解算结果就基本稳定。为便于观察，将k = 0去除后重新绘制得到岭迹图如图1(b)所示，根据图1(b)可知，在k = 0.2时，各数值已基本稳定，可取k = 0.2。

(a) (b)

Figure 1. Ridge trace

图1. 岭迹图

Table 2. Solution results

表2. 解算结果

2) 双h公式法

取 $G=I$， $h_1=t$， $h_2=0$，由式(8)可计算出k = 3.0。

现将按照经典最小二乘模型及二种岭估计的k值求解方程组，结果见表2。

3.2. 结果分析

根据解算结果可知：

1) 该模型复共线性强，最小二乘拟合结果失真，均方误差很大，计算结果不准确。

2) 该模型岭估计求解时，就均方误差而言，相对最小二乘模型均有很大优化，但岭迹法小于双h公式法解算的结果。

3) 该模型岭估计求解时，就解算精度而言，相对最小二乘模型均有很大优化，但岭迹法优于双h公式法解算的结果。

4. 结论

本文通过岭估计原理和二次曲面模型的阐述、岭估计参数的选择和实例计算，在岭估计的选择和可靠性方面得出如下结论：

第一，在病态方程组中，如果最小二乘拟合结果失真，那么均方误差将很大。采用岭估计方法对病态方程组进行处理后，可减小均方误差，并提高解算结果的可靠性。

第二，由双h公式法的原理可知，其求取的k值比较单一，通过该方法得到的高程残差约为0.06 m，精度较低，难以满足实际生产工作需要。而岭迹法在确定k值时，通过调整参数变化的步长，绘制岭迹图，所对应的解能较好地呈现其总体解的情况，本文通过选取k = 0.2，得到的高程残差绝对值小于0.035 m，可以满足实际生产工作需要，效果比双h公式法好。所以岭迹法相对于双h公式法，可直接选择均方误差较小的值作为岭参数，能够获得较高精度的解。

第三，在采用岭迹法的时候，判断岭迹图何时趋于平稳受人为因素影响，没有统一的标准，所以导致同种方法所产生的结果也可能因人而异。