现代测绘数据采集的特点包括数据量大、种类繁多等，因此在观测数据中难免会出现粗差。目前粗差探测与识别方法主要有两类：一类将粗差纳入随机模型，认为粗差与正常观测值等期望，但方差变大，由此发展了抗差估计理论^[1]。还有一类将粗差纳入函数模型，认为粗差导致该观测值期望变化，由此发展了数据探测法^[2]。施闯等^[3]提出基于相关分析进行粗差定位，后来不断扩充到由偏相关系数和复相关系数进行多维粗差定位。关于粗差定位与定值，欧吉坤^[4]提出拟准检定法(QUAD法)。於宗俦等^[5]提出LEGE方法进行多维粗差的探测和估计；林国庆等^[6]提出利用部分最小二乘(PLS)用于多维粗差的定位与定值。笔者提出一种新的确定粗差大小的方法，即先验方差待定参数法PVUP(prior-variance undecided parameter method)，其在先验方差准确时可确定粗差的大小。

1 基本原理与数学模型 1.1 基本原理

假设有一组不包含粗差的观测数据{L₁, L₂, L₃, L₄…L_n－1, L_n}，其对应的误差方程为：

$ \underset{n\times 1}{\mathop{\boldsymbol{v}}}\,=\underset{n\times t}{\mathop{\boldsymbol{A}}}\,\underset{t\times 1}{\mathop{\boldsymbol{x}}}\,-{{\left[ {{L}_{1}}\ \ \ldots \ \ {{L}_{n}} \right]}^{\text{T}}} $

(1)

这时假定L₁观测值含有粗差。去掉L₁观测值得到新的误差方程：

$ \underset{(n\text{-1})\times 1}{\mathop{{\boldsymbol{{v}'}}}}\,=\underset{(n-1)\times t}{\mathop{{\boldsymbol{{A}'}}}}\,\underset{t\times 1}{\mathop{\boldsymbol{x}}}\,-{{\left[ \begin{matrix} {{L}_{2}}\ \ \ldots \ \ {{L}_{n}} \\ \end{matrix} \right]}^{\text{T}}} $

(2)

根据式(2)可计算单位权方差$\widehat \sigma $_0old。为了估计L₁的粗差值，令y=L₁为待定参数代入式(1)：

$ \underset{n\times 1}{\mathop{\boldsymbol{v}}}\,=\underset{n\times t}{\mathop{\boldsymbol{A}}}\,\underset{t\times 1}{\mathop{\boldsymbol{x}}}\,-{{[\begin{matrix} y\ \ \ldots \ \ {{L}_{n}} \\ \end{matrix}]}^{\text{T}}} $

(3)

计算式(3)的单位权方差$\widehat \sigma $_0new= $\widehat \sigma $_0new(y)为y的函数。此时，若$\widehat \sigma $_0old与$\widehat \sigma $_0new相差不大甚至近似相等(这种假设的合理性在数学模型部分给出论述)，便可确定y的值，即L₁观测值的应有值，从而估计L₁观测值上的粗差值。这便是先验方差待定参数法的基本思想。接下来推导其数学模型。

1.2 数学模型

先验方差待定参数法需要确定先验方差。可以从3种方法获取先验方差信息：1)根据拟准观测值^[4]进行计算(如何选取拟准观测值参考文献[7])；2)根据仪器的观测精度和观测值的先验噪声大小确定；3)根据已知稳固的数据(历史数据)确定。

假设已经将观测数据正确地分成两类：1)不含粗差的数据；2)可能含有粗差待确定的数据。假设根据剔除所有怀疑含有粗差的数据计算出的单位权方差$\widehat \sigma $_0new与不含粗差观测值所计算出的单位权方差$\widehat \sigma $_0old一致(事实上，这种假设与实际也是相符的，因为粗差观测值相对于正确的观测值所占有的比例很小，不足以影响单位权方差的计算。而且多数观测值精度通常都是已知的，也可以用来作为平差的先验信息)。

假设第1类观测值向量为L, 对应的系数矩阵为A′，根据最小二乘计算出的单位权方差为$\widehat \sigma $₀。接下来以探测第n个观测值为例，说明先验方差待定参数法的具体算法。

根据AX = L + Δ可以得到误差方程:

$ \mathit{\boldsymbol{V'}} = \mathit{\boldsymbol{A'X}} - \mathit{\boldsymbol{L}} $

(4)

并根据$\sqrt {{{\mathit{\boldsymbol{V'}}}^{\rm{T}}}\mathit{\boldsymbol{V'}}/(n - m)} $计算出单位权方差估值。接下来为了探测出第n个观测值的近似准确观测值L_n，将其设为待定参数加入到式(4)中一起平差，即在A′第1行中加入L_n对应的系数向量，A′系数矩阵变成A。设y=L_n为待定参数，并将其加入到L的第1行中。令

$ \mathit{\boldsymbol{L'}} = {[\begin{array}{*{20}{c}} y&{{\mathit{\boldsymbol{L}}^{\rm{T}}}} \end{array}]^{\rm{T}}} $

(5)

则式(4)变成：

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

(6)

根据最小二乘原理可得参数X的估值：

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

(7)

令$\mathit{\boldsymbol{J}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}$，则：

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{JL'}} $

(8)

将式(8)代入式(6)得：

$ {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}} = {{\rm{(}}\mathit{\boldsymbol{L'}}{\rm{)}}^{\rm{T}}}{\rm{(}}\mathit{\boldsymbol{E}} - \mathit{\boldsymbol{AJ}}{\rm{)}}\mathit{\boldsymbol{L'}} $

(9)

其中E为单位矩阵。再令M = E － AJ得：

$ {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}} = {(\mathit{\boldsymbol{L'}})^{\rm{T}}}\mathit{\boldsymbol{ML'}} $

(10)

将M进行分块：

$ \mathit{\boldsymbol{M}} = \left( \begin{array}{l} \mathop {{\mathit{\boldsymbol{M}}_{11}}}\limits_{1 \times 1} \\ \mathop {{\mathit{\boldsymbol{M}}_{21}}}\limits_{(n - 1) \times 1} \end{array} \right.\left. \begin{array}{l} \mathop {{\mathit{\boldsymbol{M}}_{12}}}\limits_{1 \times (n - 1)} \\ \mathop {{\mathit{\boldsymbol{M}}_{22}}}\limits_{(n - 1) \times (n - 1)} \end{array} \right) $

(11)

其中, n为A的行数, 也即式(7)中的观测值数。将式(5)、(11)代入式(10)得：

$ \begin{array}{l} {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}} = {\mathit{\boldsymbol{M}}_{11}}{y^2} + ({\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_{21}} + {\mathit{\boldsymbol{M}}_{12}}\mathit{\boldsymbol{L}})y + \\ {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_{22}}\mathit{\boldsymbol{L}} = (n - m)\hat \sigma _{\rm{0}}^{\rm{2}} \end{array} $

(12)

式中最右端已经计算出来，为已知值。令

$ \left\{ \begin{array}{l} a = {\mathit{\boldsymbol{M}}_{11}}\\ b = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_{21}} + {\mathit{\boldsymbol{M}}_{12}}\mathit{\boldsymbol{L}}\\ c = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_{22}}\mathit{\boldsymbol{L}} - (n - m)\hat \sigma _0^2 \end{array} \right. $

可将(12)式化为一元二次方程。解此一元二次方程，便得到y的值，即L_n。L_n减去实际观测值便得到精确的粗差大小。这里需说明3个问题：1)该方法是直接确定观测值应有值。不难发现，如果计算得到应有值与实际相差很大，可以肯定该观测值确有粗差；相反，若相差不大，则之前分类可能出现问题或者该观测值并不含粗差。2)该一元二次方程有2个值，具体选取哪一个值可以根据观测值的合理性选取或根据其他粗差估计方式获得的值当作初值进行选取。3)这里待定参数的估值y其实质是根据不含粗差观测值的精度求得的，即根据应有的观测精度去推求含有粗差的观测值应有值，其精度自然与不含粗差的观测值精度一致，从而达到准确探测粗差大小的目的。

2 算例分析

算例1  取文献[8]的水准网，网形及相应观测值如图 1(单位m)所示。分别在l₁观测值模拟3 m、4 m、5 m、10 m粗差，且假设观测值等权独立。由文献[8]知，当观测值等权独立时，数据探测法、PLS法、QUAD法以及LEGE法是等价的，故利用本文提出的先验方差待定参数法估计l₁上的粗差，并与上述4种方法进行对比(表 1)。从表 1可以看出，当观测值独立等权的时候，先验方差待定参数法能够和其他4种方法达到同样的效果。但其先验方差待定参数并不与上述方法等价(参见算例2)，其具体的证明方法还需要进一步的研究。但是可以确定的是，先验方差精度越高，PVUP法越适用。

算例2 采用文献[4]中的测角网，用先验方差待定参数法求得的粗差估值与拟准检定进行对比(因为有学者已经证明观测值等权独立时，数据探测法、PLS法、QUAD法以及LEGE法是等价的，下文仅与QUAD法比较)。从表 2可以看出，本文的先验待定参数法与QUAD法并不等价，且优于QUAD法估计的粗差精度。

算例3  给出先验待定参数法在周跳探测里的应用。含有周跳的GPS观测值显然是一种粗差。现有原始相位观测值和含模拟周跳的相位观测值如表 3，并在第35号观测值加入10周周跳。接下来说明先验方差待定参数方法在周跳探测中的运用。根据前5个相位观测数据计算单位权方差(可利用多项式拟合法计算该方差)，该算例的单位权方差为 $\widehat{\sigma } $₀=0.093 2；然后利用第2部分数学模型的方法，计算得第35号观测值524 686.126 7，可见得到了第35号相位观测值的应有值。需要注意的是，这里与35号观测值整周数相等为巧合(因为该组相位观测值较平稳，实际可能较为复杂)。但是该算例说明了先验方差待定参数法的有效性和合理性，也说明了该方法不一定需要先确定可疑的粗差位置，只要确定先验方差即可。

3 结语

本文提出一种新的粗差估计方式——先验方差待定参数法，其在先验方差准确时可达到确定粗差大小的目的。待定参数法也可以与其他粗差探测方法结合，不局限于拟准观测，且不一定要准确探测出粗差位置，也可估计粗差估值(如算例3)。准确确定粗差值，可以提高数据利用率，同时其准确性也能用于数据质量评估等。但该方法仅考虑同精度独立观测值的情况，不同精度或者相关观测值的适用性还需进一步研究。