病态问题广泛存在于大地测量的数据处理中，如GPS快速定位、大地反演中常存在不适定问题，其原因是观测信息不充分或者对参数物理信息和先验信息不了解^[1-2]。为解决病态问题，学者们提出一系列有偏估计方法，如岭估计法、Tikhonov正则化法、截断奇异值法，这些方法主要是利用数学手段缩小病态观测方程的系数阵(或法矩阵)的条件数，但很少进行病态原因的分析和从机制上采取有效措施，没有重视算法在实际问题中的物理意义是否合理^[3-4]。解决病态问题的有效途径是充分利用参数的附加信息或先验信息对部分参数进行约束，从而保证参数解的唯一性和稳定性，如等式约束、不等式约束、区间约束以及椭球约束等^[5-8]。等式约束和不等式约束在大地测量实际问题中是一种简单的约束信息，相关文献较多，如谢建等^[6]建立附等式约束的病态模型，推导出参数估值的计算公式；冯光财等^[7]给出不等式约束和广义岭参数之间的转换方法。有学者提出利用椭球集合描述不确定先验信息，通过集员估计方法进行参数估计，如Fogel等^[8]提出的椭球集员估计算法。利用椭球约束描述复杂的先验信息，在测绘数据处理中是一种新的尝试，由于带有椭球约束的线性模型与岭估计有较强的相关性^{[7, 9]}，因此用集员估计方法处理病态问题比较合适，此方法已经在其他领域得到广泛应用，但在测量平差数据处理中还需要进一步研究^[10]。

本文将利用椭球约束描述观测向量和参数的有界不确定先验信息，把它们融入到测量平差模型中，建立平差准则对不适定性进行抑制，研究新的有界椭球不确定性平差方法在病态问题中的应用。最后，结合测边网病态数据，比较多种方法的效率。

1 区间约束和有界椭球不确定性平差模型

参数的有界不确定性可以用范数形式来描述，如‖X‖≤c；也可用区间来描述，如X∈[l, u]。这些描述方法都可以统一用一个椭球集合来表示。对于参数带区间约束的平差模型：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}},\mathit{\boldsymbol{e}} \sim N\left( {0,\sigma _0^2\mathit{\boldsymbol{P}}_0^{ - 1}} \right)}\\ {{l_i} \le {x_i} \le {u_i},i = 1,2, \cdots ,n} \end{array}} \right. $

(1)

式中，l={l₁, …, l_n}是参数X的下界，u={u₁, …, u_n}是参数X的上界。观测噪声e可以转化为区间形式，进而转化为椭球集合形式。对于X∈[l, u]，当区间不关于零点对称时，令$\mathit{\boldsymbol{c}} = \frac{1}{2}\left( {\mathit{\boldsymbol{l}} + \mathit{\boldsymbol{u}}} \right), \mathit{\boldsymbol{r}} = \frac{1}{2}\left( {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{l}}} \right)$，可以将区间约束转化为椭球约束^[7]：

$ {\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{c}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{ - 1}}\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{c}}} \right) \le 1 $

(2)

式中，$ \mathit{\boldsymbol{Q}} = {\rm{diag}}\left( {q_1^2, \cdot \cdot \cdot , q_n^2} \right)$是椭球的特征矩阵，用来刻画椭球的形状特征，则参数的有界不确定性的椭球集合表示为：

$ \begin{array}{*{20}{c}} {E\left( {\mathit{\boldsymbol{c}},\mathit{\boldsymbol{Q}}} \right) = }\\ {\left\{ {\mathit{\boldsymbol{X}}:{{\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{c}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{ - 1}}\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{c}}} \right) \le 1} \right\}} \end{array} $

(3)

观测噪声e可通过取极限误差作为上下边界转化为区间形式，同理，观测噪声e的椭球集合为：

$ E\left( \mathit{\boldsymbol{e}} \right) = \left\{ {\mathit{\boldsymbol{e}}:{\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{e}} \le 1} \right\} $

(4)

式中，P为m阶正定矩阵，是椭球的特征矩阵，用来刻画椭球的形状特征，类似于e的协方差阵描述e的特征。在测量平差中可取X₀=X_LS=(A^TA)^-1A^T L，使得L≈ AX₀，结合L= AX + e可得：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{e}} = {{\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right) \approx }\\ {{{\left( {\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_0}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_0}} \right)} \end{array} $

(5)

因此，e的有界不确定性可以近似地表示为：

$ E\left( \mathit{\boldsymbol{e}} \right) = \left\{ {\mathit{\boldsymbol{X}}:{{\left( {\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_0}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_0}} \right) \le 1} \right\} $

(6)

对于下面的有界椭球不确定性平差模型：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}}\\ {{\rm{s}}{\rm{. t}}{\rm{. }}\mathit{\boldsymbol{e}} \in E\left( \mathit{\boldsymbol{e}} \right),\mathit{\boldsymbol{X}} \in E\left( {\mathit{\boldsymbol{c}},\mathit{\boldsymbol{Q}}} \right)} \end{array}} \right. $

(7)

利用式(3)和式(6)的约束条件可以写成：

$ \mathit{\boldsymbol{X}} \in E\left( \mathit{\boldsymbol{e}} \right) \cap E\left( {\mathit{\boldsymbol{c}},\mathit{\boldsymbol{Q}}} \right) $

(8)

故(7)也可以表示成：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}}\\ {{\rm{s}}{\rm{. t}}{\rm{. }}\mathit{\boldsymbol{X}} \in E\left( e \right) \cap E\left( {\mathit{\boldsymbol{c}},\mathit{\boldsymbol{Q}}} \right)} \end{array}} \right. $

(9)

假设模型(9)的可行解集E=E(e)∩E(c, Q)的外包椭球为：

$ \begin{array}{*{20}{c}} {E\left( {{{\mathit{\boldsymbol{\hat x}}}_U},{\mathit{\boldsymbol{P}}_U}} \right) = }\\ {\left\{ {\mathit{\boldsymbol{X}}:{{\left( {\mathit{\boldsymbol{X}} - {{\mathit{\boldsymbol{\hat x}}}_U}} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}_U^{ - 1}\left( {\mathit{\boldsymbol{X}} - {{\mathit{\boldsymbol{\hat x}}}_U}} \right) \le 1} \right\}} \end{array} $

(10)

模型(9)的参数估计$\mathit{\boldsymbol{\hat{x}}}_U$和P_U^-1可以表示为^[10]：

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat x}}}_U} = {\left[ {a{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}} + (1 - a){\mathit{\boldsymbol{Q}}^{ - 1}}} \right]^{ - 1}} \cdot \\ \;\;\;\;\;\left[ {a{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_0} + (1 - a){\mathit{\boldsymbol{Q}}^{ - 1}}\mathit{\boldsymbol{c}}} \right] \end{array} $

(11)

$ \mathit{\boldsymbol{P}}_U^{ - 1} = \frac{1}{{1 - \rho }}\left[ {a{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}} + (1 - a){\mathit{\boldsymbol{Q}}^{ - 1}}} \right] $

(12)

式中，0 < a < 1，0≤ρ < 1。$\mathit{\boldsymbol{\hat{x}}}_U$是式(9)的一个解，从式(11)可以看出，$\mathit{\boldsymbol{\hat{x}}}_U$是岭估计算法的一种推广，形式上就是Schaffrin等^[11]提出的加权混合估计。令

$ \mathit{\boldsymbol{K}} = \mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{AQ}}{\mathit{\boldsymbol{A}}^{\rm{T}}} + \frac{{1 - a}}{a}\mathit{\boldsymbol{P}}} \right)^{ - 1}} $

(13)

有：

$ {{\mathit{\boldsymbol{\hat x}}}_U} = \mathit{\boldsymbol{c}} + \mathit{\boldsymbol{K}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right) $

(14)

$ {\mathit{\boldsymbol{P}}_U} = \beta \left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{KA}}} \right)\mathit{\boldsymbol{Q}} $

(15)

式中，$\beta = \frac{{\left( {1 - \rho } \right)}}{{1 - a}}$，ρ的表达式为^[10]：

$ \begin{array}{*{20}{c}} {\rho = \left( {1 - a} \right){{\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{AQ}}{\mathit{\boldsymbol{A}}^{\rm{T}}} + } \right.}\\ {{{\left. {\frac{{1 - a}}{a}\mathit{\boldsymbol{P}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right)} \end{array} $

(16)

因此，

$ \begin{array}{*{20}{c}} {\beta = \frac{1}{{1 - a}} - {{\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{AQ}}{\mathit{\boldsymbol{A}}^{\rm{T}}} + \frac{{1 - a}}{a}\mathit{\boldsymbol{P}}} \right)}^{ - 1}} \cdot }\\ {\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right) = 1 + \frac{a}{{1 - a}} - {{\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right)}^{\rm{T}}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{AQ}}{\mathit{\boldsymbol{A}}^{\rm{T}}} + \frac{{1 - a}}{a}\mathit{\boldsymbol{P}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{Ac}}} \right)} \end{array} $

(17)

最后的测量更新椭球与系数a有关，为了保证椭球交的外包椭球的最小性，即求得最小迹外包椭球需要使系数a满足：

$ a = \arg\min\;{\rm tr}\left( {{\mathit{\boldsymbol{P}}_U}} \right) $

(18)

因为0 < a < 1，0≤ρ < 1，可知a必须满足：

$ \beta = \frac{{(1 - \rho )}}{{1 - a}} > 0 $

(19)

由式(11)可知，P_U是一个正定矩阵，因此a还必须满足：

$ {\rm tr}\left[ {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{KA}}} \right)\mathit{\boldsymbol{Q}}} \right] > 0 $

(20)

系数a在满足(19)和(20)的条件下，利用直接搜索算法(a从0开始到1止，通过增量Δa逐步搜索，得到使tr(P_U)达到最小的a)求出a的值，从而求出参数估计$\mathit{\boldsymbol{\hat{x}}}_U$和P_U。

2 算例分析 2.1 算例1

采用文献[2]中的病态矩阵A与观测向量L，L中加入了e~N(0, σ²I)、σ=0.1的随机误差。如式(21)所示，参数的真$ \mathit{\boldsymbol{\tilde X}} = {\left\{ {1, 1, 1, 1, 1} \right\}^{\rm T}}$，目标函数为$F = {\left( { \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A \hat X}}} \right)^{\rm T}} \mathit{\boldsymbol{P}}\left( { \mathit{\boldsymbol{L - A\hat X}}} \right)$，参数估值与真实值间的残差平方和为${{m}}\left( {\mathit{\boldsymbol{\hat X}}} \right) = {\left\| {\mathit{\boldsymbol{\hat X}} - \mathit{\boldsymbol{\tilde X}}} \right\|^2}$。

$ \mathit{\boldsymbol{A = }} \\ \left[ {\begin{array}{*{20}{c}} {2.000\;0}&{ - 5.000\;0}&{1.000\;0}&{1.000\;0}&{ - 9.000\;0}\\ { - 2.000\;0}&{4.000\;0}&{1.000\;0}&{ - 1.050\;0}&{8.500\;{\rm{0}}}\\ { - 2.000\;0}&{1.000\;0}&{1.000\;0}&{ - 1.000\;0}&{{\rm{2}}{\rm{.400}}\;{\rm{0}}}\\ { - 1.000\;0}&{2.500\;0}&{4.000\;0}&{ - 0.500\;0}&{{\rm{7}}{\rm{.000}}\;{\rm{0}}}\\ { - 1.000\;0}&{3.200\;0}&{4.000\;0}&{ - 0.500\;0}&{{\rm{8}}{\rm{.400}}\;{\rm{0}}}\\ {1.000\;0}&{1.000\;0}&{ - 3.000\;0}&{0.400\;0}&{{\rm{0}}{\rm{.490}}\;{\rm{0}}}\\ {3.000\;0}&{7.000\;0}&{ - 3.000\;0}&{1.500\;0}&{{\rm{12}}{\rm{.700}}\;{\rm{0}}}\\ {5.000\;0}&{ - 1.000\;0}&{ - 2.000\;0}&{2.500\;0}&{ - 3.000\;0}\\ {4.000\;0}&{2.000\;0}&{ - 2.000\;0}&{2.010\;0}&{3.000\;0}\\ {4.000\;0}&{3.000\;0}&{ - 2.000\;0}&{2.000\;0}&{5.000\;0} \end{array}} \right], \\ \mathit{\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} { - 10.653}\\ {10.7535}\\ {1.400\;0}\\ {12.000\;0}\\ {14.100\;0}\\ { - 0.110\;0}\\ {21.200\;0}\\ {1.500\;0}\\ {9.010\;0}\\ {12.00\;0} \end{array}} \right] $

(21)

已知e~N(0、σ²I)、σ=0.1，取极限误差-3σ和3σ作为边界值，即e_l={-0.3, -0.3, -0.3, -0.3, -0.3}^T，e_u={0.3, 0.3, 0.3, 0.3, 0.3}^T，e∈[e_l, e_u]，则P=diag(0.9, 0.9, 0.9, 0.9, 0.9)，e的椭球集合为E(e) ={ e:e^TP^-1e≤1}；对于参数X，直接取X∈[l, u]，l ={-2, -2, -2, -2, -2}^T，u={4, 4, 4, 4, 4}^T，则Q=diag(45, 45, 45, 45, 45)。

经典最小二乘法的参数估值$\mathit{\boldsymbol{\hat X}}_{\rm LS}$=(1.094 8, 3.663 4, 1.658 3, 0.845 7, -0.315 0)^T，与真值相比偏离很大，残差平方和为m=9.289 2，严重失真；有界椭球不确定性的结果$\mathit{\boldsymbol{\hat x}}_U$=(1.012 0, 1.026 1, 1.007 1, 0.995 0, -0.997 4)^T，残差平方和为m=0.000 9，与经典最小二乘法相比，结果改善较多，且尽管先验区间的上下界与真值相差较大，但$\mathit{\boldsymbol{\hat x}}_U$与真值的偏差非常小，这也说明了有界椭球不确定性算法在处理病态问题上的有效性。

2.2 算例2

测边网数据引自文献[12]，该测边网得到的系数矩阵是病态矩阵，P1~P10点是已知点，其坐标以及与未知点之间的观测边长如表 1所示，K1、K2、K3为未知点(真值分别为(68, -26, 9)、(14, 41, -11)以及(0, 0, 0))，各观测量中误差为e~N(0, σ²I)、σ=0.01。3个未知点之间的观测边长分别为d_K1K2=88.340 2，d_K1K3=73.355 1，未知点改正数参数估值为$\mathit{\boldsymbol{\hat X}} = {\left\{ {{{\hat x}_1}, {{\hat y}_1}, {{\hat z}_1}, {{\hat x}_2}, {{\hat y}_2}, {{\hat z}_2}, {{\hat x}_3}, {{\hat y}_3}, {{\hat z}_3}} \right\}^{\rm{T}}}$，单位为m。测边网如图 1所示。

已知e~N(0, σ²I)、σ=0.01，同算例1一样取上下边界值，并进一步求得椭球特征矩阵P=diag(0.019 8, 0.019 8, 0.019 8, 0.019 8, 0.019 8)，则e的椭球集合为E(e) ={ e:e ^TP^-1 e≤1}。对于3个未知点的三维坐标，将前期观测得到的三维坐标作为近似坐标，分别为(68.010, -25.990, 8.980)、(14.020, 40.990, -11.010)以及(0.010, -0.010, 0.020)，近似坐标对应的点位精度分别为(0.01, 0.01, 0.02)、(0.02, 0.01, 0.01)以及(0.01, 0.01, 0.020)，取3倍点位精度作为边界值，求得的椭球特征矩阵Q=diag(0.008 1, 0.008 1, 0.032 4, 0.032 4, 0.008 1, 0.008 1, 0.008 1, 0.008 1, 0.032 4)。

表 2是测边网算例的比较结果，可以看出：

1) LS的结果与真值偏离较大，$m(\mathit{\boldsymbol{\hat X}})$=32.367 5，严重失真。

2) 有界椭球不确定性算法的残差平方和最小，而且远小于利用L曲线确定岭参数的岭估计。已知有界椭球不确定性算法是岭估计的一种推广，其对应的参数是利用先验信息来确定的，在先验信息可靠的前提下，可以说明有界椭球不确定性算法充分利用了参数的有界不确定性，比只使用数学原理的岭估计能更好地解决病态问题。

3) 椭球约束是将椭球特征矩阵的逆作为广义岭估计的岭参数，同样利用了参数的有界不确定性，其残差平方和与有界椭球不确定性相差不大，验证了充分利用参数先验信息处理病态问题的有效性。

4) 虽然TSVD的残差平方和也很小，但其截断参数取k=1，而系数矩阵的奇异值为2.694 3>2.494 4>2.061 7>1.898 9>1.224 1>1.074 7>0.093 5>0.006 1，明显截断了过多信息，不符合实际。

3 结语

对于病态问题，充分利用参数的先验信息可以有效抑制系数矩阵的复共线性，降低解的不稳定性。基于这种思想，本文将有界椭球不确定性应用于解决病态问题。不确定性先验信息一般是以范数、区间、集合或模糊数形式表达的，而不是一个具体的数值，有界椭球不确定性将参数和观测向量的先验信息用2个椭球集合来表示，以这2个椭球外包椭球的特征矩阵的迹最小为准则进行不确定性平差，可以充分利用复杂先验信息，从而能够保持解的稳定性。