在测量过程中，受设备、观测者及环境等影响，会产生不确定性，使观测值在真值左右波动^[1-2]。不确定性包含偶然误差、系统误差、粗差、数值和概念误差、可度量和不可度量误差^[3]。不确定度是不确定性的度量，可以用方差、均方差、误差区间、误差椭圆等表示^[4]。

文献[4-7]研究了不确定性和不确定度基础理论，文献[8-10]对不确定度的估计和应用作了研究。文献[11]对系数矩阵和观测向量同时存在不确定性的问题建立不确定性min-max平差准则，先利用奇异值分解(SVD)得到中间参数的方程，再解方程得到未知参数估值。特别地，对系数矩阵部分列存在不确定性的部分不确定性min-max平差准则，先用QR分解将其转换成一般形式，再用SVD和解方程得未知参数估值。文献[12]针对不确定性min-max平差准则，提出一种基于线性矩阵不等式(LMI)的参数估计方法。文献[5]使用不确定性min-max平差准则建立不确定性平差模型，而对于SVD后得到的中间参数方程，采用迭代算法得到未知参数估值。文献[13]给出部分不确定性平差模型，同样采用SVD后的迭代算法求解中间参数。整个求解过程过于复杂，且没有讨论迭代不收敛时的求解方法。

为提高部分不确定性平差模型的解算效率，本文在不使用QR分解和SVD的基础上，直接采用迭代算法进行求解，并讨论迭代不收敛时的求解方法。

1 平差模型及准则

设系数矩阵为且为满秩矩阵，观测向量为，未知参数向量为。顾及A中某列含有常数元素，文献[13]给出部分不确定性平差的函数模型：

(1)

式中，为A的常数元素部分；为A中含有不确定性元素的部分；为A₂的不确定性；为L的不确定性；η_L和η_A分别是ΔL和ΔA₂的不确定度^[5]，且均大于0；和为X的子块；‖·‖_F为F-范数。

针对式(1)，采用min-max准则^{[11, 13]}进行解算:

(2)

式中，符号“^”表示最优估值。该平差准则分为max(·)和min(·)内外两层，内层的目标函数max(·)表示在系数矩阵子块A₂和观测向量L分别受有界不确定性ΔA₂和ΔL干扰情况下，得到未知参数时，残差的最大不确定性；外层的目标函数min(·)表示使参数解中的不确定性最小化。

将ΔA₂和ΔL看成自变量，根据范数性质，其极大值点就是最大值点，设为(ΔA₂⁰，ΔL⁰)。当^[13]

(3)

时，max(·)取得极大值，亦即最大值。所以，平差准则(2)可等价转换成^{[11, 13]}：

(4)

针对平差准则(4)，文献[11]先对系数矩阵A用QR分解，将其转化为一般形式；再对目标函数求导，得到未知参数估计式；然后用SVD得到中间参数方程，并解出唯一正根；最后将其代入未知参数估计式，得到未知参数估值，简称QR分解-SVD-解方程法。而文献[13]直接对目标函数求导，对A进行SVD，采用迭代算法解出中间参数估值，简称SVD-迭代算法。以上方法解算过程都较为复杂，不易实现。本文提出一种更为简单的方法，直接利用迭代算法解算未知参数估值。

2 模型解算 2.1 直接迭代算法

当迭代算法收敛时，可直接采用迭代算法求解平差准则(4)。令目标函数为：

(5)

问题(4)转变为求解含有F-范数的二元函数在其定义域上的最小值。由范数性质易知，在全局定义域内连续可导，且二阶导数大于0，所以它的极小值点就是其最小值点，而函数的极小值点在其驻点或者不可导处取得。根据求函数自由极值的方法，分别对和求偏导数，并令各自偏导数等于0，得：

(6)

(7)

联合式(6)和式(7)，并参考文献[13]，先讨论的不可导点：1)假设，这意味着系数矩阵A和观测向量L的不确定性都不存在，与本文条件不符，所以假设不成立。2)假设X₂ =0，此时系数矩阵A中只剩下常数元素部分A₁参与解算。若观测向量L只含有随机误差，则可用最小二乘平差方法求解式(6)，得；若L还含有其他的不确定性，则选取相应的平差方法进行求解。

进一步讨论的驻点，令

由式(6)，得：

(8)

式中，A₁为列满秩矩阵，从而N₁₁为满秩矩阵，所以N₁₁的逆矩阵一定存在且唯一。由式(8)得：

(9)

由式(7)得：

(10)

式中，

(11)

将式(9)代入式(10)，整理得：

(12)

式中，I_k为k阶单位矩阵。式(12)两端都含有未知参数估值，可用迭代算法求解。

部分不确定性平差模型的迭代算法步骤为：

1) 给定的初始值，再联立式(8)、式(10)和式(11)，并将已知数值代入，可得未知参数初始估值；

2) 用获得的和已知值代入式(11)，得：

代入式(12)，得：

再将其代入式(9)，得：

3) 重复步骤2)，直到(ε为迭代算法的收敛限差)时，迭代终止。

2.2 QR分解-SVD-解方程法

当迭代不收敛时，本文对QR分解-SVD-解方程法进行讨论^[11]。对系数矩阵A进行QR分解：

(13)

式中，为正交矩阵，为上三角矩阵，其分块矩阵为，再对Q^TL进行分块：

(14)

式中，。将式(13)和式(14)代入式(4)，并由F-范数的正交不变性，式(4)等价于如下形式的平差准则：

(15)

在式(15)中，根据范数性质，若解算出未知参数估值，则求的最小值。令，得：

(16)

令，由式(15)得:

(17)

对A进行SVD：

(18)

式中，，均为正交矩阵；Σ=diag(σ₁, σ₂, …, σ_k)为对角矩阵，且对角元素σ₁≥σ₂…≥σ_k>0是A的全部奇异值。再对U^TL分块：

(19)

式中，。引入函数：

(20)

并令

(21)

设L₂≠0，若η_A≥τ，则；若η_A < τ，则，式中的是Φ(α)=0的唯一正实数根。

3 算例分析

在以下算例中，所用计算机CPU为Intel(R) Core(TM) i7-6500U，主频为2.50 GHz，内存为8 GB；采用软件Matlab R2012a(7.14.0.739)，64 bit。系统误差为常数s，随机误差由函数normrnd产生，期望为μ、标准差为σ。简记DI为本文提出的直接迭代算法，QSE为QR分解-SVD-解方程法，SI为SVD-迭代算法。取ε=10^－12。

3.1 直线拟合算例

用直线方程y=5.813 2－0.601 7x生成一组10个不含误差的点(x_i，y_i)，模拟真值见表 1。再向其中同时加入系统误差(s=－0.10)和随机误差(μ=0，σ=0.41)，并满足η_A=η_L=1.30，得到同时含有这两种误差的系数矩阵A和观测向量L：

共模拟1 000次。采用DI、QSE和SI分别进行参数估计，得到DI和SI的迭代次数以及3种算法的参数估值、计算时间，见表 2、图 1和图 2。

从表 2第2栏看出，DI得到的未知参数估值与另外两种算法一致，说明DI是正确可行的。

图 1统计了DI和SI的迭代次数，DI每次计算的迭代次数都要比SI少。表 2第3栏列出DI和SI计算1 000次的平均迭代次数，DI的迭代效率比SI提高约33.08%。图 2统计了3种算法每次计算的时间，DI每次都是最少的。表 2第4栏还列出它们计算1 000次的总时间，DI总的时间效率较QSE提高约94.41%，较SI提高约54.91%。分析可知，DI的解算效率比另外两种算法都有明显提高。

从图 1知，DI总体迭代次数的波动范围比SI小，说明DI的稳定性比SI好；从图 2知，DI计算时间的波动范围比QSE和SI都要窄，说明DI计算时间的稳定性也更好。

3.2 平面坐标转换算例

假设平面坐标转换的函数模型为：

给定转换参数真值[3.708 6 －5.216 7 0.988 2 －0.139 1]^T，及4个不含误差的旧坐标点(x_i，y_i)，由此生成4个新坐标点(x_i^′，y_i^′)，见表 3。给新旧坐标同时加入系统误差(s=－0.05)和随机误差(μ=0，σ逐次增大)，使得不确定度η_A=η_L的初始值为0.40，步长为0.40，终止值为3.60，共产生9组同时含有这两种误差的新旧坐标，得到系数矩阵A和观测向量L：

每组模拟1 000次，采用DI、QSE和SI分别对其进行参数估计，统计DI和SI迭代次数的均值以及3种算法的计算时间，见图 3和图 4。

从图 3可以看出，随着不确定度的增加，DI和SI的迭代次数也有不同程度的增加，且DI的迭代次数较少；它们的均值分别为约4.607 8和6.220 2，方差分别为约0.718 0和1.242 3，DI迭代次数的均值和方差都较小，说明DI较SI迭代效率较高，且更为稳定。

观察图 4，随着不确定度的逐次增加，DI、QSE和SI的计算时间也随之增加。它们计算时间的均值分别为约0.103 7、6.456 3和0.328 2，方差分别为约0.000 3、0.939 1和0.003 2。显然，DI计算时间的均值和方差都是最小的，较其他两种算法计算的时间效率最高，且最为稳定。

4 结语

本文推导了部分不确定性平差模型的直接迭代算法，并讨论了迭代不收敛时的算法。得到几点结论：1)直接迭代算法是正确可行的，且概念简单，易于实现。2)直接迭代算法有着更快的收敛速度、更高的解算效率、更好的稳定性。3)当迭代不收敛时，迭代算法会失效，此时可用QR分解-SVD-解方程法解算，但过程较为复杂，解算效率较低。

本文主要是从解算效率方面进行探讨，而系数矩阵和观测向量不确定度的获取以及部分不确定性平差模型的应用等，还需进一步研究。