测绘学报 | 王乐洋：加乘性混合误差模型精度评定的SUT法

本文内容来源于《测绘学报》2022年第11期（审图号GS京（2022）1229号）

加乘性混合误差模型精度评定的SUT法王乐洋^1,2, 陈涛^1,3     1. 东华理工大学测绘工程学院, 江西 南昌 330013;
2. 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室, 江西 南昌 330013;
3. 武汉大学测绘学院, 湖北 武汉 430079基金项目：国家自然科学基金(42174011;41874001); 东华理工大学研究生创新专项资金(DHYC-202020)摘要：已有加乘性混合误差模型参数估计方法能达到二阶精度, 但精度评定方法只能达到一阶精度, 若通过传统泰勒级数展开近似函数法来获取参数估值的二阶精度信息, 由于加乘性混合误差模型中参数估值与观测值为一个复杂的非线性关系, 必然需要复杂的求导运算。针对该问题, 本文使用一种无须求导、无须了解非线性函数构成的比例无迹变换(scaled unscented transformation, SUT)法来计算参数估值的二阶精度信息。通过算例分析表明, 利用SUT法求解加乘性混合误差模型能够有效避免复杂的求导运算, 所求得的参数估值及其协方差阵均能达到二阶精度, 从而验证了本文方法的可行性和优势。关键词：加乘性混合误差模型    加权最小二乘    非线性函数    SUT法    精度评定   

引文格式：王乐洋, 陈涛. 加乘性混合误差模型精度评定的SUT法[J]. 测绘学报，2022，51(11)：2303-2316. DOI: 10.11947/j.AGCS.2022.20200514 WANG Leyang, CHEN Tao. The SUT method for precision estimation of mixed additive and multiplicative random error model[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(11): 2303-2316. DOI: 10.11947/j.AGCS.2022.20200514 

阅读全文：http://xb.chinasmp.com/article/2022/1001-1595/20221107.htm

引 言

在大地测量数据处理领域中，测量平差理论是在基于加性误差的线性或非线性高斯-马尔可夫(Gauss-Markov, GM)^[1-2]模型基础上发展起来的。随着现代测绘技术的高速发展，以电磁波为代表的观测值的随机误差表现为乘性误差或加乘性混合误差^[3-4]，例如LiDAR观测值的随机误差表现为加乘性混合误差，电磁波测距(electronic distance measurement, EDM)、全球导航卫星系统(GNSS)和甚长基线干涉测量(very long baseline interferometry, VLBI)基线精度计算中，乘常数和加常数本质上也是加乘性混合误差^[3-7]。传统的加性误差模型平差理论不再适用于新型加乘性混合误差模型，若直接采用基于加性误差的平差理论进行解算，模型不能正确反映观测值随机误差的特征，此时求得的结果将不可靠，因此需要针对加乘性混合误差模型进行更深入的研究来丰富现代大地测量数据处理理论。

目前，加乘性混合误差模型的解法包括拟极大似然法和最小二乘法^[5]。虽然拟极大似然法是广义线性模型或乘性误差模型参数估计的主要方法，但它存在以下问题^[5]：①没有明确定义的目标函数；②拟似然解是由一个非线性方程组决定的，可能存在多组解，无法确定最似然解。乘性误差模型的最小二乘法最早是由文献[8]提出，文献[5]在乘性误差模型的基础上最早给出了加乘性混合误差模型的函数模型和随机模型，并基于最小二乘原理推导了最小二乘解法和加权最小二乘迭代解法，此外还通过公式推导将文献[8]中乘性误差模型的偏差改正加权最小二乘迭代解法推广到加乘性混合误差模型，但该方法是通过参数估值的迭代公式根据误差传播定律求得参数估值的精度信息，由于在进行泰勒函数展开时忽略了求偏导数的二阶及其高阶项，故只能达到一阶精度。文献[9]给出了乘性误差和加性误差彼此独立且各自具有相同方差情况下的偏差改正加权最小二乘迭代解法，但该方法也只能求得参数估值的一阶精度。因此需要针对加乘性混合误差模型参数估值的精度评定问题作进一步的研究。

无迹变换(unscented transformation, UT)方法是一种确定性采样策略的近似非线性函数概率密度分布法，最早由文献[10—11]在非线性Kalman滤波领域中首次提出，其通过选择特定的采样策略来近似非线性函数的概率分布，无须求导、无须了解非线性函数的构成，通过加权的方法便可计算非线性函数的均值和协方差阵。经过20多年的发展，不同的学者提出了多种采样策略，扩展了其应用范围。文献[12—15]提出了一种适应性强的比例对称采样策略，相应的UT即为SUT。文献[13]证明了通过SUT法计算的非线性函数的均值和协方差阵均可以达到二阶精度。文献[16]证明了将SUT法应用在测量非线性函数误差传播中是可行和有效的，并通过与二阶近似函数法进行对比，证明了SUT法求得的参数估值及其精度信息均可以达到二阶精度。文献[17]使用SUT法研究了变量误差(errors-in-variables, EIV)模型的参数估计和精度评定问题。文献[18]研究了部分变量误差(partial errors-in-variables, Partial EIV)模型方差分量估计精度评定问题。文献[19]研究了Okada模型中震源参数估值对滑动分布反演中格林函数矩阵的影响分析并进行精度评定。文献[20]研究了Okada模型的震源参数的精度评定问题。EIV模型、Partial EIV模型及Okada模型中的参数估值和观测值之间为一个复杂的非线性迭代过程，加乘性混合误差模型中参数估值和观测值之间也是一个非常复杂的非线性迭代过程，因此，SUT法可以较好地用于加乘性混合误差模型的精度评定。

综上所述，针对已有加乘性混合误差模型精度评定方法只能达到一阶精度，且加乘性混合误差模型中参数估值与观测值为一个复杂的非线性关系，无法通过求导运算和公式推导进行精度评定的问题，本文使用一种无须求导、无须了解非线性函数构造的SUT法对其进行精度评定，旨在求得参数估值的二阶精度信息；最后，通过算例验证本文将SUT法应用于加乘性混合误差模型精度评定问题中的可行性及优势。

1 加乘性混合误差模型及其解法

1.1 加乘性混合误差模型

加乘性混合误差模型可以表示为^[5]

 (1)

式中，y∈R^n×1表示观测值向量；f(·)表示未知参数的函数；β∈R^t×1表示未知参数向量；⊙表示Hadamard积符号；1∈R^n×1表示元素全为1的列向量；ε_m∈R^n×1表示服从正态分布的乘性误差列向量；ε_a∈R^n×1表示服从正态分布的加性误差列向量。

本文研究的函数f(·)是β的线性函数形式^[4-7, 9]，则式(1)可表示为

 (2)

式中，A∈R^n×t表示系数矩阵，且A=[a₁, a₂, …, a_t](a_i∈R^n×1)。

假设乘性误差向量ε_m和加性误差向量ε_a彼此独立且各自具有相同的方差^[4-6, 9]，即D_m=σ_m²I_n，D_a=σ_a²I_n，cov(ε_m, ε_a)=0，根据期望和协方差的定义^[21]，由式(2)可求得加乘性混合误差模型观测值的期望和协方差阵分别为^[4-6, 9]

 (3)

 (4)

式中，。

由式(4)可知，在加乘性混合误差模型中，观测值的协方差阵D(y)是参数估值β的非线性函数，这明显区别于传统的加性误差模型，即观测值的协方差阵D(y)与参数估值β无关，因此传统的加性误差模型平差理论不再适用于新型的加乘性混合误差模型。

1.2 加乘性混合误差模型的最小二乘解

对式(2)应用最小二乘准则，可得目标函数为

 (5)

根据求解目标函数自由极值的原理，将函数F₁(β)对β求偏导，并令其为0，可得加乘性混合误差模型的最小二乘解为

 (6)

由误差传播定律可得最小二乘估值_LS的协方差阵为

 (7)

由式(6)和式(7)可得，最小二乘解_LS为无偏估计，但其方差已不是最优的，不具有方差最小性^[9]。

1.3 加乘性混合误差模型的加权最小二乘迭代解

对式(2)应用加权最小二乘准则，可得目标函数为

 (8)

式中，σ₀²表示先验单位权方差；表示观测值的权。

根据求解目标函数自由极值的原理，将函数F₂(β)对β求偏导，并令其为0，可得

 (9)

式中，_WLS表示参数估值的加权最小二乘估值。根据矩阵微分定理^[22]可得，(i=1,2...,t)。

由上述分析可知，由于加乘性混合误差模型中观测值的权是参数估值的非线性函数，因此通过式(9)无法得到参数估值的解析解，只能使用数值逼近的方法进行迭代求得数值解，进而可以得到加乘性混合误差模型的加权最小二乘迭代解为

 (10)

式中，P(y)ⁱ表示第i次迭代观测值的权；(r=1,2,...,t)。

由文献[5]可得加权最小二乘迭代解的单位权中误差和一阶协方差阵分别为

 (11)

 (12)

1.4 加乘性混合误差模型的偏差改正加权最小二乘迭代解

文献[5]通过公式推导证明了加权最小二乘解的偏差等于观测值的权对参数估值求偏导数项，因此，在加权最小二乘解的基础上删去观测值的权对参数估值求偏导数项(即式(10)右边第二项)得到了参数估值的二阶无偏估计，对应的偏差改正加权最小二乘解为

 (13)

式(13)中观测值的权是参数估值的非线性函数，而参数估值又是观测值的权的非线性函数，因此式(13)也为一个迭代解，其迭代格式为

 (14)

式中，_bcWLSⁱ⁺¹表示参数估值的偏差改正加权最小二乘迭代解。

由文献[5]可得偏差改正加权最小二乘迭代解的单位权中误差和一阶协方差阵分别为

 (15)

 (16)

由式(16)可以看出，文献[5]是通过参数估值的迭代公式和协方差传播率来获取参数估值的精度信息，由于忽略了求偏导数的二阶及其高阶项，故只能达到一阶精度。

综上所述，本文将文献[5]中加乘性混合误差模型的偏差改正加权最小二乘迭代解法总结为算法1，具体步骤如下。

(1) 输入观测值y，系数矩阵A，计算最小二乘估值作为迭代初值

 (17)

(2) 计算第i次迭代观测值的权P(y)ⁱ

 (18)

 (19)

 (20)

(3) 由式(14)计算加乘性混合误差模型的偏差改正加权最小二乘迭代解。

(4) 重复步骤(2)至步骤(3)，当满足‖_bcWLSⁱ⁺¹－_bcWLSⁱ‖≤10^－6时停止循环，输出参数估值_bcWLSⁱ⁺¹。

(5) 由式(15)和式(16)分别计算偏差改正加权最小二乘迭代解的单位权中误差和一阶协方差阵。

2 加乘性混合误差模型精度评定的SUT法

由算法1可以看出，偏差改正加权最小二乘迭代解法中每次迭代的参数估值与观测值的函数关系相同，因此在经过第i次迭代后，最终得到的参数估值的迭代过程可以表示为

 (21)

式中，_LS表示由式(6)求得的最小二乘估值作为迭代初值；φ(·)表示偏差改正加权最小二乘迭代解法中参数估值与观测值进行一次迭代的函数关系，即式(14)；F(·)表示最终得到的参数估值与观测值之间的函数关系；y∈R^n×1表示观测值即自变量，其均值和协方差阵分别为E(y)和D(y)；Y∈R^t×1表示最终求得的参数估值，即因变量，其均值和协方差阵分别为E(Y)和D(Y)。

式(21)表明最终得到的参数估值与观测值之间为一个复杂的非线性函数关系，此时若通过近似函数法来进行精度评定，必然需要复杂的求导运算。为了避开复杂的求导运算和公式推导，同时得到参数估值的二阶精度信息，本文在偏差改正加权最小二乘迭代解法中引入无须求导、无须了解非线性函数构成的SUT采样法进行精度评定，旨在获得参数估值的二阶精度信息。由于式(21)是一个复杂的非线性函数，无法将其具体表达式表示出来，但其本质上也是一个非线性函数，因此，不失一般性，现从理论上证明通过SUT法求得一般的函数Y=F(y)的因变量的均值和协方差阵均可以达到二阶精度，具体如下。

在式(21)中，因变量Y对自变量y的Jacobian矩阵和Hessian矩阵可表示为^[23]

 (22)

 (23)

 (24)

进一步对函数Y=F(y)进行二阶泰勒展开可得

 (25)

式中，Δy=y－E(y)；G=[Δy^Tb₁ΔyΔy^Tb₂Δy…Δy^Tb_tΔy]^T。

根据期望和协方差阵的定义^[21]，对式(25)进行求期望和协方差阵运算可得

 (26)

 (27)

式中，⊗表示Kronecker积^[24]。

由式(26)和式(27)可以看出，若通过传统泰勒级数展开的近似函数方法进行运算，需要求得Jacobi矩阵a和Hessian矩阵c，但当非线性函数复杂且非线性强时，很难进行求取。因此，为了确保因变量的均值和协方差阵在误差传播中能达到二阶精度，并且避免求导，利用SUT法求取因变量的均值和协方差阵，对于函数Y=F(y)，SUT法精度评定的具体步骤为^[12-13]。

(1) 基于自变量y的先验均值E(y)和协方差阵D(y)依据式(28)生成采样点集{χ_i}(i=0, 1, …, 2n)

 (28)

式中，表示对D(y)进行Cholesky分解生成的下三角矩阵；为矩阵的第i列；参数λ=α²(n+κ)－n，其中α取值范围为[10^－4, 1]^[12]，本文中α的取值为10^－3[14]；κ为另一个比例常数，在本文中κ的取值为0^[14, 25]。

(2) 将采样点集{χ_i}(i=0, 1, …, 2n)代入函数Y=F(y)中，求得对应的解集

 (29)

(3) 加权计算因变量Y的均值E(Y)^SUT和协方差阵D(Y)^SUT

 (30)

 (31)

式中，W_i^m和W_i^c表示的是样本点集v_i(i=0, 1, …, 2n)的权值，具体可表示为^[13] ，，(i=1, 2, …, 2n)；对于正态分布，γ=2为最优取值^[13]。

将通过SUT法求得的因变量Y的均值E(Y)^SUT和协方差阵D(Y)^SUT，即式(30)和式(31)具体展开可得^[16]

 (32)

 (33)

式中，J为Y在E(y)处展开的n维一阶偏导向量；H为Y在E(y)处展开的n×n阶海森矩阵；R_m^SUT和R_c^SUT分别代表四阶及以上的高阶项；v_i=(i=1, 2, …, n)；I^j_δiF为χ_i在E(y)处展开的j阶偏导。

由文献[24]可得tr(M₁M₂M₃M₄)=vec(M₃^T)(M₂^T⊗M₄)vec(M₁)，tr(M₁M₂)tr(M₁M₂)=tr(M₁M₂M₁M₂)。因此，式(32)和式(33)进一步可表示为

 (34)

 (35)

由于在SUT法中，α=10^－3，γ=2，κ=0，因此式(35)进一步可表示为

 (36)

对比式(26)与式(34)可知，两者表达式的形式完全相同，表明利用SUT法求得函数Y=F(y)的均值可以精确地达到二阶精度；同时，对比式(27)与式(36)可知，两者表达式的形式也完全相同，表明利用SUT法求得函数Y=F(y)的协方差阵也可以精确地达到二阶精度。因此，表明本文利用SUT法求解加乘性混合误差模型，所求得的参数估值及其协方差阵均能达到二阶精度。

综上所述，本文将加乘性混合误差模型精度评定的SUT法总结为算法2，具体步骤如下。

(1) 依据文献[20, 26—29]的做法，将加乘性混合误差模型的观测值y作为采样的先验统计均值E(y)；将通过算法1求得的偏差改正加权最小二乘迭代解_bcWLSⁱ⁺¹代入式(4)求得的观测值的后验协方差阵D(y)作为采样的先验协方差阵。然后，将均值E(y)和协方差阵D(y)代入式(28)生成采样点集{χ_p^y}(p=0, 1, …, 2n)(即2n+1组观测值)。

(2) 将2n+1组观测值{χ_p^y}(p=0, 1, …, 2n)依次代入式(37)(即算法1)求得对应的解集(p=0, 1, …, 2n)(即2n+1组参数估值的解)

 (37)

(3) 由式(30)和式(31)加权计算参数估值_SUT及其协方差阵D(_SUT)。

(4) 由式(38)和式(39)计算参数估值的单位权中误差和标准差

 (38)

 (39)

由上述步骤可以看出，SUT法能够解决复杂或难以实现的求导计算，具有操作简单、易于实现等特点，且利用SUT法求解加乘性混合误差模型能够有效避免复杂的求导运算，所求得的参数估值及其协方差阵均能达到二阶精度。

3 算例与分析

为验证本文提出方法在加乘性混合误差模型中精度评定方面的性能及探讨在大地测量领域中的应用，本文共设计了3个算例：通过算例1来展示本文方法在GPS高程拟合中的应用价值，算例2和算例3来展示本文方法在数字地面模型(DTM)中的应用价值，并验证本文方法在精度评定方面的优势。分别使用最小二乘解法、加权最小二乘迭代解法、文献[5]中的偏差改正加权最小二乘迭代解法(即算法1)及加乘性混合误差模型精度评定的SUT法(即算法2)进行对比解算以验证本文方法的可行性及优势，4种方案对应的方法分别列于表 1中。

表 1 4种方案及其对应的方法

Tab. 1 Four schemes and their corresponding methods

表选项 

3.1 算例1

为初步验证本文方法在加乘性混合误差模型中的可行性和正确性，算例1采用文献[30]中的GPS高程拟合算例，该算例中高程真值与离某一坐标原点的距离除以100满足以下函数模型

 (40)

式中，y表示高程真值；x表示离坐标原点的距离，研究区域为300 m之内，间隔为10 m；未知待求参数为式(40)的4个系数，即β=[10 1 2 2]^T。

假设高程真值受加乘性混合误差的影响，具体可表示为^[7]

 (41)

式中，Y表示31维高程观测值向量；y表示31维高程真值向量；1表示元素全为1的31维列向量；ε_m表示独立且服从正态分布的31维乘性误差列向量；ε_a表示独立且服从正态分布的31维加性误差列向量。

根据文献[5—7]，本算例将乘性误差的标准差设置为0.05，加性误差的标准差设置为0.15 m，单位权中误差设置为σ₀=0.3 m，随机误差由Matlab中的mvnrnd函数生成。表 2列出了地面点的高程模拟值，为了说明高程值受加乘性混合误差的影响程度，将不含随机误差时的高程真值及受加乘性混合误差干扰的高程值绘于图 1中。

表 2 地面点的高程值

Tab. 2 The elevations of ground points

表选项 

图 1 不含随机误差时的高程真值及受加乘性混合误差干扰的高程值

Fig. 1 The true elevations without random error and the elevations disturbed by mixed additive and multiplicative random error

图选项 

由图 1可知，当高程受到加乘性混合误差的干扰时，观测值会受到一定的影响，进而与原始的高程真值产生了偏离。

分别使用表 1中的4种算法计算，图 2给出了WLS法和bcWLS法的流程，图 3给出了SUT法的流程。将参数真值、参数估值及其与真值的二范数和单位权中误差列于表 3，参数估值的标准差列于表 4(表中表示参数估值，表示参数估值与真值之间的二范数，表示单位权中误差，表示参数估值的标准差)。

图 2 WLS法和bcWLS法的流程

Fig. 2 Flowchart of WLS method and bcWLS method

图选项 

图 3 SUT法的流程

Fig. 3 Flowchart of SUT method

图选项 

表 3 算例1中LS法、WLS法、bcWLS法及SUT法求得的参数估值、估值与真值之间的二范数和单位权中误差

Tab. 3 The parameter estimations, two norms and unit weight errors obtained by the LS method, WLS method, bcWLS method and SUT method of the first example

表选项 

表 4 算例1中LS法、WLS法、bcWLS法及SUT法求得的参数估值的标准差

Tab. 4 The standard deviations of parameter estimations obtained by the LS method, WLS method, bcWLS method and SUT method of the first example

表选项 

由表 3的结果可知，LS法由于在平差时未考虑观测值的权，因此求得的二范数结果最大，且单位权中误差也严重偏离真值；WLS法由于在LS法的基础上考虑了观测值的权，因此求得的二范数及单位权中误差结果均优于LS法；bcWLS法由于在WLS法的基础上减去了其偏差，因此求得的二范数的结果优于WLS法，表明文献[5]中方法的有效性。对比bcWLS法和SUT法可以看出，两者求得的二范数结果一样，而文献[5]已通过严格的数学公式推导证明了bcWLS法求得的参数估值能达到二阶精度，说明本文通过SUT法求得的参数估值也能达到二阶精度，表明本文方法的可行性。

由表 4求得的标准差结果可知，LS法求得的结果最大，WLS法的结果小于LS法，而bcWLS法的结果略小于WLS法，精度最高，这与文献[5]得到的结论一致，进一步表明在平差时考虑观测值权的必要性。然而，由于bcWLS法是基于参数估值的最后一步迭代表达式，通过协方差传播率求得参数估值的精度信息，因此只能达到一阶精度。对比bcWLS法和SUT法求得的结果可以看出，SUT法求得的标准差最小，精度最高，而上文已证明利用SUT法求得加乘性混合误差模型参数估值的协方差阵可以达到二阶精度，说明通过SUT法求得的参数估值的二阶精度信息的情形下中误差是减小的，这与文献[28—29]得到的结论一致。此外，若参数估值的二阶精度信息情形下中误差是增大的，此时参数估值的二阶标准差将大于bcWLS法求得的标准差，而WLS法求得的标准差也大于bcWLS法求得的标准差，说明通过WLS法求得的标准差为二阶精度，但这显然是不合理的，因此进一步说明通过SUT法求得参数估值的二阶精度信息的情形下中误差是减小的，从而验证了本文方法的可行性和优势。

3.2 算例2

为进一步验证本文方法在加乘性混合误差模型中的可行性和正确性，算例2是一个DTM模型算例，该模型通过有限的地形高程数据实现对地面地形的数字化模拟，被广泛用于工程、军事、遥感与制图、地理分析等方面，例如在工程项目中的挖填方计算、线路勘测设计、航天遥感数字图像定量解释中的应用等^[31-32]。目前，机载LiDAR被广泛用于DTM数据的获取^[33-34]，但是以LiDAR技术为数据源构建DTM模型都是基于加性误差来进行处理的^[35-37]，然而LiDAR数据的噪声已被证明具有乘性噪声性质^[38-40]，且文献[5—7, 9]已将LiDAR数据的误差当作加乘性混合误差进行处理，因此本算例将仿真设计受加乘性混合误差干扰的LiDAR DTM模型试验，以更好地验证本文方法的优势。本算例模拟的含4个未知参数的加乘性混合误差的DTM模型观测方程如下

 (42)

 (43)

式中，(x, y)表示地面点坐标，横坐标x和纵坐标y都在0~100 m之间，以步长为2 m取值；H(x, y)表示DTM模型高程观测值；h(x, y)表示DTM模型高程真值；1表示元素全为1的列向量；ε_m表示独立且服从正态分布的乘性误差列向量；ε_a表示独立且服从正态分布的加性误差列向量；β_i表示未知待求参数；函数f_i(·)(i=1, 2, 3, 4)的具体形式为

 (44)

根据文献[38—40]中对LiDAR数据的理论和实际测试及已有关于加乘性混合误差模型的文献[5—7]，本算例将加性误差的标准差固定为σ_a=0.15 m，将乘性误差的标准差分别设置为σ_m=0.005、σ_m=0.05和σ_m=0.1这3组值进行对比解算。为了说明DTM模型受加乘性混合误差的影响程度，将不含随机误差时的DTM模型绘于图 4中，将受到加乘性混合误差干扰的DTM模型分别绘于图 5、图 6和图 7中。令单位权中误差σ₀=0.3 m^[5-7]，分别使用表 1中的4种算法计算，将参数真值、参数估值及其与真值的二范数和单位权中误差列于表 5，参数估值的标准差列于表 6。

图 4 不含误差时的数字地形模型

Fig. 4 The digital terrain model without random errors

图选项 

图 5 乘性误差σ_m=0.005时的数字地形模型

Fig. 5 The digital terrain model in the case σ_m=0.005

图选项 

图 6 乘性误差σ_m=0.05时的数字地面模型

Fig. 6 The digital terrain model in the case σ_m=0.05

图选项 

图 7 乘性误差σ_m=0.1时的数字地面模型

Fig. 7 The digital terra in model in the case σ_m=0.1

图选项 

表 5 算例2中LS法、WLS法、bcWLS法及SUT法求得的参数估值、估值与真值之间的二范数和单位权中误差

Tab. 5 The parameter estimations, two norms and unit weight errors obtained by the LS method, WLS method, bcWLS method and SUT method of the second example

表选项 

表 6 算例2中LS法、WLS法、bcWLS法及SUT法求得的参数估值的标准差

Tab. 6 The standard deviations of parameter estimations obtained by the LS method, WLS method, bcWLS method and SUT method of the second example

表选项 

通过对比图 4、图 5、图 6和图 7可以发现，随机误差对高程产生了较大的影响，并且随着乘性误差变大，DTM模型的高程也严重偏离真值，因此针对加乘性混合误差的DTM模型需要进行更深入的研究来丰富现代大地测量数据处理理论。

通过比较表 5中4种算法求得的参数估值与真值之间的二范数和单位权中误差结果可以看出，相比于LS法，WLS法、bcWLS法及SUT法求得的参数估值更接近于参数真值，且单位权中误差也更接近真值，表明考虑观测值的权能够得到较优的结果。同时，相比于WLS法，bcWLS法能够在一定程度上减小参数估值的偏差得到较优的参数估值，表明文献[5]中方法的有效性。通过对比bcWLS法和SUT法的结果可以看出，SUT法能够得到与文献[5]中方法相同的参数估值结果，这与算例1得到的结论一致，进一步表明本文通过SUT法求得的参数估值能达到二阶精度。此外，对比不同乘性误差求得的结果可以发现，当乘性误差σ_m=0.005时，4种算法所得参数估值与真值较为接近；当乘性误差σ_m=0.05和σ_m=0.1时，4种算法所得参数估值逐渐偏离参数真值，表明乘性误差的大小会影响所求得的参数估值结果。

通过比较表 6中4种算法求得的参数估值的标准差可以看出，LS法求得的结果最大，WLS法求得的结果远小于LS法，bcWLS法求得的结果略小于WLS法的结果，而SUT法求得的结果最小，精度最高，这与算例1得到的结论一致，进一步表明文献[5]中方法只能求得参数估值的一阶精度，而本文方法能求得参数估值的二阶精度。同时，对比不同乘性误差求得的结果可以发现，随着乘性误差变大，4种算法所得的参数估值的精度也逐渐变差，但SUT法总能求得最优的参数估值的精度信息，进一步表明将SUT法与加乘性混合误差模型结合是可行且有效的。

3.3 算例3

为进一步验证本文方法在实测数据中的适用性及优势，算例3通过在网站http://srtm.csi.cgiar.org/srtmdata/提供的LiDAR型DTM模型数据下载了塔里木盆地864个点的三维坐标数据进行解算^[7, 29]。同算例2一样，受加乘性混合误差干扰的DTM模型观测方程见式(42)，不同的是DTM模型函数，算例3采用文献[7, 9, 29, 41]中的DTM模型函数拟合为

 (45)

式中，(x, y)表示地面点坐标；h(x, y)表示DTM模型高程值；β_i表示未知待求参数。

目前，由于已有关于加乘性混合误差模型的研究成果较少，表现为加乘性混合误差特征的观测值未得到足够重视，且实际观测数据中加性误差和乘性误差的分析是很复杂的，需要进行专门的论文研究，因此本算例暂不讨论分离实际观测数据中的加性误差和乘性误差，而是将获取的实测数据视为真值通过模拟误差进行解算。为了模拟出实际地形，本算例将乘性误差的标准差设置为0.002，加性误差的标准差设置为σ_a=0.15 m，单位权中误差同算例1和2一样。将不含随机误差时的DTM模型绘于图 8中，受到加乘性混合误差干扰的DTM模型绘于图 9中。由算例1和2的结果可知，LS法和WLS法的结果不足以反映参数估值及其真实精度，因此算例3不再使用LS法和WLS法求解参数估值及其精度信息，仅通过bcWLS法与SUT法进行对比来验证本文方法的可行性及优势。针对求解过程中出现的病态奇异情况，采用岭估计法^[42-43]来解决法矩阵求逆不稳定的问题，所得的参数估值的精度信息被列于表 7，拟合得到的DTM模型如图 10和图 11所示。在本算例中，将获取的DTM模型高程真值与模型反演值之间的均方根误差(root mean square error, RMSE)的大小来比较两种方法的优势，具体可表示为

 (46)

图 8 不含误差时的数字地形模型

Fig. 8 The digital terrain model without random errors

图选项 

图 9 受加乘性混合误差干扰后的数字地形模型

Fig. 9 The digital terrain model disturbed by mixed additive and multiplicative random error

图选项 

图 10 不含误差时的数字地形模型(点)及bcWLS法求得的数字地面模型(圈)

Fig. 10 The digital terrain model without random errors (point) and the digital terrain model obtained by the bcWLS method (circle)

图选项 

图 11 不含误差时的数字地形模型(点)及SUT法求得的数字地面模型(圈)

Fig. 11 The digital terrain model without random errors (point) and the digital terrainmodel obtained by the SUT method (circle)

图选项 

式中，n表示观测值个数；h_true, i表示获取的DTM模型高程真值；c_mod, i表示模型反演的高程值。

由图 8可以看出，本文获取的DTM模型为一个陡峭的曲面，高程最高点和最低点的偏差为18.892 8 m；同时，通过对比图 8和图 9中的DTM模型可以看出，尽管本算例加入的乘性误差的标准差仅为0.002，但其对DTM模型的高程产生了明显的偏差。此外，由图 10和图 11可以看出，两种方法均能较好地拟合原始的DTM模型，但通过仔细对比也可以发现，SUT法较bcWLS法能够更好地拟合原始的DTM模型，表明本文方法在实测数据处理中的优势。

对比表 7中的均方根误差可以看出，SUT法求得的结果较bcWLS法在数值上更小，表明本文算法能够更好地拟合原始的DTM模型。同时，对比参数估值的标准差可以看出，两种方法均能求得较优的参数估值的精度信息，但由于SUT法为二阶精度评定方法，因此求得的参数估值的精度明显优于bcWLS法，进一步表明本文方法在实测数据处理中的优势。

表 7 算例3中bcWLS法和SUT法求得的参数估值的标准差及均方根误差

Tab. 7 The standard deviations and mean square matrix error of parameter estimations obtained by the bcWLS method and SUT method of the third example

表选项 

4 结论

本文首先阐述了加乘性混合误差模型的统计性质及其误差结构，然后分析了已有加乘性混合误差模型的参数估计方法能达到二阶无偏，但精度评定方法只能达到一阶精度，且无法通过近似函数法来获取参数估值的二阶精度信息，最后通过结合一种无须求导、无须了解非线性函数构成的SUT法对加乘性混合误差模型进行精度评定，并设计了具体的算法流程。

通过对本文算例中的各方法对比分析得出，当乘性误差较小时，本文算法求得的参数估值的二阶精度信息略优于已有方法求得的参数估值的一阶精度信息；而当乘性误差较大时，此时通过SUT法求得的参数估值的二阶精度信息明显优于已有方法；此外，将SUT法应用于实测DTM模型中，可以得到比已有方法更小的均方根误差结果，且参数估值的精度更高，表明将SUT方法用于解决加乘性混合误差模型的精度评定问题是可行且有效的。

本文仅是在理论方面对加乘性混合误差模型精度评定问题进行探讨，并从理论上针对加乘性混合误差模型已有解法的不足给出了一种新的算法，为今后研究LiDAR数据、SAR图像等^[44-45]实际观测数据的乘性误差特性，以及结合地空观测的实践与应用提供坚实的理论基础，而如何将本文方法与方差分量估计法^[46]结合进而切实运用于实际中是接下来要研究的内容。

致 谢

感谢许光煜和邹传义提出的宝贵建议和帮助。

作者简介

第一作者简介：王乐洋(1983—), 男, 博士, 教授, 研究方向为大地测量反演及大地测量数据处理。E-mail: [e-mail]

通信作者：陈涛, E-mail：[e-mail]

初审：张艳玲

复审：宋启凡
终审：金   君

往期推荐

资讯