当观测值不含粗差、观测误差服从零均值的分布时，经典最小二乘估计具有最优无偏性。然而，在实际的数据采集中经常会出现粗差，由于最小二乘不具有抗差性，对粗差观测值十分敏感，少量的粗差就可能对最小二乘估计结果产生破坏性影响^[1-6]。因此，在粗差干扰出现的情况下，选择适当的估计方法，尽量降低或者去除粗差对参数估值的影响，得出正常模式下的最优或接近最优的参数估值，即为测量数据处理领域稳健估计的基本思想。对于现代空间测量技术如GNSS、摄影测量与遥感等，利用稳健估计处理自动化采集、海量的测量数据中混入的粗差具有重要的理论和实际意义。

1953年，文献[7]首次提出了稳健估计(又称为抗差估计)这一专业术语。文献[8]提出了污染分布模式。文献[9]阐述了定位参数的抗差估计。文献[10]提出了影响函数和崩溃点的概念。这些都为稳健估计理论奠定了基础。1980年，文献[11]将稳健估计理论引入测量界，提出了著名的“丹麦法”。此外，国内学者在稳健估计方面也取得重大突破，如文献[12]提出等价权的概念；文献[13]提出自适应的稳健估计算法等。目前，稳健估计可以分为3大类：M估计、L估计、R估计。其中，M估计是经典极大似然估计的推广，是测量数据处理中稳健估计最重要的估计方法。

稳健估计的选权迭代算法属于M估计，在测量中应用广泛。基本原理为：首先采用最小二乘法得到观测值的估计残差；根据残差确定各观测值新的权因子，进行下一轮平差计算；反复迭代，直至前后两次估计值的变化值小于限差为止。算法每次迭代求解法方程均需进行矩阵的求逆运算，矩阵求逆的复杂度是其维数的三次方。因此，平差模型的待估参数多，即法方程正交矩阵维数高，算法的计算量和计算时间急剧上升，算法效率难以满足现代测量中大规模测量方案、海量数据及平差模型实时解算等要求。

本文基于矩阵逆的运算法则，提出了一种有效改进现有选权迭代算法的方法。迭代计算中，仅需计算由前一次残差确定的新权阵下的模型未知参数估计值的变化量，避免了法方程正交矩阵的求逆，显著提高了算法的效率，并能保持原稳健权函数的抗差性。本文以水准网为例，通过增加估计参数的数量和粗差数量，验证了当模型的估计量增加时，改进的选权迭代算法的计算效率要远高于经典选权迭代算法。

1 经典选权迭代算法

设平差的线性函数模型和随机模型为

式中，l为观测值向量；B为系数矩阵并且列满秩；x为未知参数向量；Δ为误差向量；D为对角方差-协方差矩阵；σ₀为单位权中误差；Q为对角协因数矩阵；P为对角权阵。

独立不等权观测情况下的M估计准则为

式中，b_i=[b_i1   b_i2  …  b_it]；l_i表示第i个观测值，相应的误差方程为

对式(2)中x求导，记，可得

令p_i=p_iw_i，，则

式中，P为等价权阵；p_i为等价权元素；w_i为权因子^[12]。

将式(3)代入式(5)，得到M估计下的法方程

解式(6)得

式中，为的协因数矩阵；；w=B^TPl，将式(7)代入式(3)可得观测值的改正数v。

根据残差v确定各观测值新的权因子，再进行下一轮平差，反复迭代直至前后两次估计值的变化小于限差，计算结束。算法的核心是确定权函数ρ(或φ)，即由每次迭代的残差估计值确定新的权阵，以保证估值的抗差性和效率^[13]。以下介绍几种常见的权函数。

(1) Huber法^[14-17]

(2) Danish法^[18-20]

(3) IGGⅢ方案^[21-23]

c₁=1.5, c₂=3.0

以上权函数中，w(u)为权函数，c₁和c₂为相应的调和系数。

2 改进选权迭代算法

经典选权迭代算法中，每次由残差v的估计值确定新的权阵计算时均需求逆运算，矩阵求逆的复杂度约达到矩阵维数的三次方的量级，远高于矩阵的线性运算。在平差模型的待估参数的个数多、粗差数多的情况下，算法的计算效率低、计算时间大大增加。改进算法采用矩阵逆的运算法则，将分解为两部分，一部分为迭代前的待估参数的协因数阵，另一部分为根据权函数确定新的权阵后协因数阵的增量矩阵，从而避免了经典选权迭代计算中的协因数阵的求逆运算。

选权迭代算法中，每次迭代运算确定了观测值新的权值后(每次仅改变一个权值)，将其增量记为对角矩阵dP，则待估参数新的协因数阵为

式中

将式(9)代入式(8)得

根据矩阵求逆引理^[24-25]

将式(10)展开

式中，表示权阵更新前的参数的协因数阵；表示权阵更新后的协因数阵的增量矩阵，可表示为

式(7)中，令w_old、w_new分别表示权阵更新前后的权阵，则

式中，dw为

令

式中，表示迭代前参数估计值；表示新权阵下的参数增量估计值；表示权阵更新后的参数估计值。

由式(12)—式(16)得

式中，；α=d_i^Tb_i；β=d_i^Tw_old；γ=dp_i^－1+b_i^Td_i，λ=dp_il_i。

则式(13)可表示为

改进选权迭代算法步骤总结如下：

(1) 估计最小二乘初值：

根据平差模型的输入量l、B、P(Q)，采用最小二乘估计初值。

(2) 选权迭代计算：

由v和式(15)计算权阵更新后的dw，式(18)计算；

由式(16)、式(17)计算和；

直到小于设定的限值，迭代终止，得到最终估计结果。

改进选权迭代算法每次迭代只需计算增量和dw，避免了矩阵的求逆，在平差模型的待估参数或者粗差数量较多的情况下，大大降低了计算量，提升了计算效率。改进算法不改变权函数，保持了权函数的抗差性，对独立不等精度观测值的选权迭代算法具有普适性。

3 实例分析

为了从数值上说明改进选权迭代算法在计算效率上的优势，笔者设计了一个条状、结构规整的水准网实例，如图 1所示。点A、B、C、D为已知控制点，其余为未知点，未知点个数为t(t为偶数)，观测值个数为n, 且n=1.5t+2。

水准网设计试验数据数据说明如下：

控制点高程：H_A=1 m, H_B=0 m, H_C=1.5 m, H_D=0.5 m

高差观测值：

h₁=h₂=h₇=h₈=…=(0.5+Δ)m

h₄=h₅=h₁₀=h₁₁=…=(－0.5+Δ)m

h₃=h₆=h₉=h₁₂=…=(1+Δ)m

式中，Δ为符合正态分布(μ=0，σ=0.005)的随机误差，设计粗差的量级约为1 m。

当水准网包含不同数量级的未知点t=100、500、1000、5000、10 000时，且分别包含1到5个粗差的情况下，笔者采用Matlab编程对改进选权迭代算法与经典算法的计算效率进行了比较。程序采用Matlab自带的tic/toc函数，记录了算法改进前后程序运行所需时间，试验结果见表 1。

试验结果表明：

(1) 改进算法计算效率要高于经典选权迭代算法(图 2、图 3)，尤其是在平差模型中待求参数个数t较大或者多个粗差的情况下，计算效率显著提升。当未知数个数t=10 000且含5个粗差时，改进算法所需计算时间仅为经典算法的1/30。

(2) 改进算法的计算效率取决于采用最小二乘估计初值的计算时间，迭代计算无需求逆，当粗差个数增加时，计算时间仅有微小增长。例如，当未知数个数t=10 000，粗差从1个增长到5个时，经典算法计算时间急剧增长，由829 s增加4倍，达到3285 s；而改进算法从100 s仅增加到109 s，改进算法对于处理多个粗差情况非常高效。

4 结论

经典选权迭代算法属于M估计，是当前运用最为广泛的稳健估计方法之一。本文采用矩阵逆运算法则，提出了一种改进的选权迭代算法。改进算法在迭代过程中，仅需计算新权阵下的最小二乘估计改正值，避免了经典算法中每步对平差模型法方程的求逆运算，显著提高了算法的计算效率。改进选权迭代算法适用于现代空间测量技术下对自动观测的海量数据的实时或者准实时的处理需求，如GNSS、遥感等观测数据，当这些数据含有一定数量的粗差时，改进算法可以大大降低经典选权迭代算法的计算时间。本文算法基于独立不等精度观测值的情况，笔者下一步的研究目标是针对海量数据，讨论如何提高相关观测值^{[5, 17, 26]}的选权迭代算法的计算效率。