测量控制网中的粗差是在数据获取、传送和加工过程中随机出现的大的误差，只有极个别的观测值可能含粗差。如果不能在平差中正确地发现、消除粗差，可能会较显著地歪曲网的平差结果而不能使用。粗差探测法就是在网平差中自动发现、定位和剔除粗差^{[1, 2]}。目前，粗差探测的方法主要有两种：一种是以观测值残差 (v_i) 为研究对象，包括假设检验法和抗差估计，其中假设检验的统计仍是最小二乘结果，因而其本身的抗差性能较差，而抗差估计的关键在于选择合适的极值函数，实际应用中，虽然抗差极值函数有几十种，但目前还没有被普遍适用的；另一种是以观测值真误差估值为研究对象，典型代表是欧吉坤的拟准检定法，事实上，拟准检定法第一步仍然是对所有观测值进行经典的最小二乘平差处理^[3]。刘根友提出了粗差探测的两步法：第1步是将观测值粗差作为未知参数，采用拟稳平差思想解算粗差；第2步是选取部分观测值为拟准观测值，采用部分最小二乘原理确定粗差^{[4, 5]}。本文在借鉴刘根友的两步法的基础上，进行了进一步研究，提出了在不等权情况下的带权平差因子概念和粗差定值定位的方法，并用一个水准网平差的例子进行了详细的说明。

1 带权平差因子

高斯-马尔柯夫线性模型^[6]为：

$ \mathit{\boldsymbol{ }}\!\!\Delta\!\!\rm{ }=\mathit{\boldsymbol{A\tilde{X}}}-\mathit{\boldsymbol{L}}, \mathit{\boldsymbol{ }}\!\!\Delta\!\!\rm{ }\sim N(0\ \ {{\sigma }^{2}}\mathit{\boldsymbol{Q}}) $

(1)

观测方程为：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{AX}}-\mathit{\boldsymbol{L}}\;\;\;或\;\;\;\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{\varepsilon }}, \mathit{\boldsymbol{P}} = {\mathit{\boldsymbol{Q}}^{-1}} $

(2)

其中，L为n维观测向量；Δ为观测值的真误差；V为观测值的改正数；ε为观测噪声；A为系数矩阵；X为m维未知参数向量；为$ \mathit{\boldsymbol{\tilde{X}}} $的真值；P为观测值的权阵。

由最小二乘原理V^TPV=min得：

$ {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{PV}} = 0 $

(3)

则

$ \mathit{\boldsymbol{X}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})^{-1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PL}} $

(4)

将式 (4) 代入式 (2) 得：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{A}}{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})^{- 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PL}}- \mathit{\boldsymbol{L}} = \left[{\mathit{\boldsymbol{A}}{{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{-1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}-\mathit{\boldsymbol{I}}} \right]\mathit{\boldsymbol{L}} $

式中，I为单位阵。定义：

$ \begin{array}{l} \mathit{\boldsymbol{J}} = \mathit{\boldsymbol{A}}{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})^{-1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}\\ \mathit{\boldsymbol{R}} = \mathit{\boldsymbol{I}}-\mathit{\boldsymbol{J}} \end{array} $

为带权平差因子，则有：

V=－RL

由于带权平差因子J、R具有最小二乘的数学意义^[7]，因此一旦使用，就表明对观测方程式 (2) 实施了最小二乘准则 (V^TPV=0)，由平差因子和观测值可以直接求得改正数V。

若对式 (1) 实施最小二乘准则 (Δ^TPΔ=min)，同样得到Δ=－RL，此时Δ不再是真误差。由于都实施了最小二乘准则，此时的Δ与V虽然符号不同，但都是代表真误差的估值。

2 基于拟稳平差思想进行粗差定位的原理

拟稳平差是周江文教授于1980年提出的一种新的平差方法，用于秩亏自由网平差，它的基本思想是：将整个网点分为两组，第1组是相对稳定的点，称之为拟稳点；第2组是相对不稳定的点，称之为动点，并对第1组的拟稳点施以最小范数条件X_F^TX_F=min，从而消除秩亏，解算出未知量^{[7, 9]}。本文利用拟稳平差思想，选取部分观测值为拟准观测值，并施以最小范数条件，使粗差出现明显分群现象，从而达到粗差定位的目的。

由于事先不知道哪些观测值含粗差、哪些不含粗差，这里假设所有观测值均含有粗差，设粗差为n维向量G(G=(g₁g₂…g_n)^T)，这种假设的合理性在于，当任意一观测值 (第i个) 不含粗差时，则令g_i=0。

若将粗差也作为待估参数，则同时含有粗差、偶然误差的观测方程表示为：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{AX}}-\mathit{\boldsymbol{G}}-\mathit{\boldsymbol{L}} $

(5)

式 (5) 为秩亏方程，秩亏数设为m。由最小二乘原理得:

$ \frac{\partial ({{\mathit{\boldsymbol{V}}}^{\rm{T}}}\mathit{\boldsymbol{PV}})}{\partial \mathit{\boldsymbol{X}}}=2{{\mathit{\boldsymbol{V}}}^{\rm{T}}}\mathit{\boldsymbol{P}}\frac{\partial \mathit{\boldsymbol{V}}}{\partial \mathit{\boldsymbol{X}}}=2{{\mathit{\boldsymbol{V}}}^{\rm{T}}}\mathit{\boldsymbol{PA}}=0 $

即

$ {{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PV}}=0 $

将式 (5) 代入上式得:

$ {{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PAX}}-{{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PG}}-{{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PL}}=0 $

(6)

由最小二乘原理得:

$ \frac{\partial ({{\mathit{\boldsymbol{V}}}^{\rm{T}}}\mathit{\boldsymbol{PV}})}{\partial \mathit{\boldsymbol{G}}}=2{{\mathit{\boldsymbol{V}}}^{\rm{T}}}\mathit{\boldsymbol{P}}\frac{\partial \mathit{\boldsymbol{V}}}{\partial \mathit{\boldsymbol{G}}}=-2{{\mathit{\boldsymbol{V}}}^{\rm{T}}}\mathit{\boldsymbol{P}}=0 $

即

$ \mathit{\boldsymbol{PV}}=0 $

(7)

将式 (5) 代入式 (7) 得：

PAX－PG－PL=0

上式两边左乘P^－1得：

G=AX－L

若A^TPA可逆，由式 (6) 得：

$ \mathit{\boldsymbol{X}}={{({{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{-1}}{{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{P}}(\mathit{\boldsymbol{G}}+\mathit{\boldsymbol{L}}) $

(8)

$ \mathit{\boldsymbol{RG}}=-\mathit{\boldsymbol{RL}} $

(9)

由G=AX－L可知，G与式 (2) 中的V相同。

由于R秩亏，无法直接求解粗差G，须附加其他约束条件。这里，根据逆稳平差的基本思想，对粗差进行约束，同时把不含粗差的观测值所对应的粗差未知数约束为0。约束方程个数不能小于m。考虑一般性，观测方程式 (5) 附加条件约束，可表示为：

$ \left\{ \begin{matrix} \mathit{\boldsymbol{V}}=\mathit{\boldsymbol{AX}}-\mathit{\boldsymbol{G}}-\mathit{\boldsymbol{L}}=(\begin{matrix} \mathit{\boldsymbol{A}} &-\mathit{\boldsymbol{E}} \\ \end{matrix})\left( \begin{matrix} \mathit{\boldsymbol{X}} \\ \mathit{\boldsymbol{G}} \\ \end{matrix} \right)-\mathit{\boldsymbol{L}} \\ \left( \begin{matrix} {{\mathit{\boldsymbol{C}}}_{1}} & {{\mathit{\boldsymbol{C}}}_{2}} \\ \end{matrix} \right)\left( \begin{matrix} \mathit{\boldsymbol{X}} \\ \mathit{\boldsymbol{G}} \\ \end{matrix} \right)=0 \\ \end{matrix} \right. $

(10)

由附加条件最小二乘原理得：

$ \left\{ \begin{align} & {{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PAX}}-{{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PG}}+{{\mathit{\boldsymbol{C}}}_{1}}^{\rm{T}}\mathit{\boldsymbol{K}}={{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PL}} \\ &-\mathit{\boldsymbol{AX}}+\mathit{\boldsymbol{G}}+{{\mathit{\boldsymbol{C}}}_{2}}^{\rm{T}}\mathit{\boldsymbol{K}}=-\mathit{\boldsymbol{L}} \\ & {{\mathit{\boldsymbol{C}}}_{1}}\mathit{\boldsymbol{X}}+{{\mathit{\boldsymbol{C}}}_{2}}\mathit{\boldsymbol{G}}=0 \\ \end{align} \right. $

(11)

如果K=0，条件方程系数阵一定是m行满秩，且 (C₁ C₂)(A －E)^T=0，则可导出：

X=(A^TPA)^－1A^TP(G+L)

G=AX－L或RG=－GL或V=0

当C₁=0，C₂=A^T，即$ \mathit{\boldsymbol{C}}=\left( \begin{matrix} \underset{m\times m}{\mathop{0}}\, & {{\underset{m\times n}{\mathop{\mathit{\boldsymbol{A}}}}\, }^{\rm{T}}} \\ \end{matrix} \right) $时，即A^TG=0，此时，C中部分列的元素全部为0，而秩仍为m，则对应最小范数解G^TG=min。则可得到：

$ \left\{ \begin{matrix} \mathit{\boldsymbol{X}}={{({{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{-1}}{{\mathit{\boldsymbol{A}}}^{\rm{T}}}\mathit{\boldsymbol{PL}} \\ \mathit{\boldsymbol{G}}=-\mathit{\boldsymbol{RL}} \\ \end{matrix} \right. $

X与最小二乘解相同，G即为最小二乘残差。

当$ $，且r>m时，即A_r^TG_r=0，这里的r称之为拟准观测数，相应观测值为拟准观测值，此时部分最小范数条件为：

G_r^TG_r=0

结合粗差检定的两种途径思路分析，确定本文粗差检定的步骤如下：

第1步，认为所有粗差都满足条件G^TG=min，即所有观测值均为拟准观测值，此时解算出的G就是经典的最小二乘解V。

第2步，计算粗差位值：

$ {{\sigma }_{s}}=b\frac{\sum{\left| {{g}_{i}} \right|}}{n} $

(12)

系数b取0.8, 可将所有观测值分为两类，即拟准观测值与非拟准观测值。取v_i小于σ_s的观测值为拟准观测值。

第3步，利用拟准观测值作部分约束，此时将矩阵C₂中q_i=0所对应的列元素都定义为零，重新解算法方程，得到验后单位权中误差为：

$ {{\sigma }_{r}}=\sqrt{\frac{{{\mathit{\boldsymbol{V}}}_{r}}^{\rm{T}}{{\mathit{\boldsymbol{V}}}_{r}}}{r-m}} $

(13)

取残差小于3σ_r的观测值为拟准观测值。

3 实例验证

某水准网如图 1所示，由两个已知水准点、3个水准点构成，观测了9段高差，其中h₆存在粗差^[10]。由于本方法中是将粗差作为未知参数，故有3+9=12个未知参数。设

$ \mathit{\boldsymbol{C}}=\left( \begin{matrix} {{\mathit{\boldsymbol{C}}}_{1}} & {{\mathit{\boldsymbol{C}}}_{2}} \\ \end{matrix} \right)=\left( \begin{matrix} \underset{3\times 3}{\mathop{0}}\, & \underset{3\times 9}{\mathop{{{\mathit{\boldsymbol{A}}}^{\rm{T}}}}}\, \\ \end{matrix} \right) $

引入辅助向量q_i(i=1, 2, …, 9), 当q_i=0，C₂中对应i列元素均为0；当q_i=1时，C₂中对应i列元素保持不变，表明相应观测值为拟准观测值。

具体计算过程如下：

1) 令q_i=1(i=1, 2, 3, …, 9)，即所有观测值为拟准观测值，计算结果与经典最小二乘平差一致，结果见表 1。

2) 计算粗差位值:

$ {{\sigma }_{s}}=0.8\frac{\sum{\left| {{g}_{i}} \right|}}{n}=2.71 $

取小于粗差位值σ_s的观测值为拟准观测值，q_i=1(i=1, 2, 7, 8, 9)。

3) 利用拟准观测值作部分约束，此时将矩阵C₂中q_i=0(i=3, 4, 5, 6) 所对应的列元素都定义为零，重新解算法方程，得验后单位权中误差σ_r=1.40，取小于3σ_r的观测值为拟准观测值，有q_i=1(i=1, 2, 3, 5, 7, 8, 9)。

4) 将矩阵C₂中q_i=0(i=4, 6) 所对应的列元素全定义为零，重新解算法方程，得到验后单位权中误差σ_r=0.82。大于3σ_r的观测值为非拟准观测值，即为含粗差观测值。

4 结束语

粗差是随机出现的大的误差，只有极个别的观测值可能含粗差。本文采用拟稳平差方法，假设所有观测值均含粗差，采用最小二乘原理，得到的是传统间接平差方法处理结果，由于观测值的粗差在平差后的观测值的改正数中会得到一定反映，且多小于其粗差，故可以根据观测值的改正数的大小来确定拟准观测值，对拟准观测值施以最小范数条件，直到拟准观测值和含粗差的观测值的改正数出现明显分群现象，最后相当于删除了含粗差的观测值，粗差根据拟准观测方程的参数平差值计算得到，因此粗差定位的可靠性较好。但前提条件是要有足够多的多余观测值，否则该方法也会失效。所以任何粗差探测和定值定位方法都不是万能的，对于控制网来说，其合理布设更为重要。