2. 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室, 江西 南昌 330013;
3. 武汉大学测绘学院, 湖北 武汉 430079
2. Key Laboratory of Mine Environmental Monitoring and Improving around Poyang Lake, Ministry of Natural Resources, Nanchang 330013, China;
3. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China
在大地测量数据处理领域中，测量平差理论是在基于加性误差的线性或非线性高斯-马尔可夫(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∈Rn×1表示观测值向量；f(·)表示未知参数的函数；β∈Rt×1表示未知参数向量；⊙表示Hadamard积符号；1∈Rn×1表示元素全为1的列向量；εm∈Rn×1表示服从正态分布的乘性误差列向量；εa∈Rn×1表示服从正态分布的加性误差列向量。
本文研究的函数f(·)是β的线性函数形式[4-7, 9]，则式(1)可表示为
(2)
式中，A∈Rn×t表示系数矩阵，且A=[a1, a2, …, at](ai∈Rn×1)。
假设乘性误差向量εm和加性误差向量εa彼此独立且各自具有相同的方差[4-6, 9]，即Dm=σm2In，Da=σa2In，cov(εm, εa)=0，根据期望和协方差的定义[21]，由式(2)可求得加乘性混合误差模型观测值的期望和协方差阵分别为[4-6, 9]
(3)
(4)
式中，
由式(4)可知，在加乘性混合误差模型中，观测值的协方差阵D(y)是参数估值β的非线性函数，这明显区别于传统的加性误差模型，即观测值的协方差阵D(y)与参数估值β无关，因此传统的加性误差模型平差理论不再适用于新型的加乘性混合误差模型。
1.2 加乘性混合误差模型的最小二乘解对式(2)应用最小二乘准则，可得目标函数为
(5)
根据求解目标函数自由极值的原理，将函数F1(β)对β求偏导，并令其为0，可得加乘性混合误差模型的最小二乘解为
(6)
由误差传播定律可得最小二乘估值
(7)
由式(6)和式(7)可得，最小二乘解
对式(2)应用加权最小二乘准则，可得目标函数为
(8)
式中，σ02表示先验单位权方差；
根据求解目标函数自由极值的原理，将函数F2(β)对β求偏导，并令其为0，可得
(9)
式中，
由上述分析可知，由于加乘性混合误差模型中观测值的权是参数估值的非线性函数，因此通过式(9)无法得到参数估值的解析解，只能使用数值逼近的方法进行迭代求得数值解，进而可以得到加乘性混合误差模型的加权最小二乘迭代解为
(10)
式中，P(y)i表示第i次迭代观测值的权；
由文献[5]可得加权最小二乘迭代解的单位权中误差和一阶协方差阵分别为
(11)
(12)
文献[5]通过公式推导证明了加权最小二乘解的偏差等于观测值的权对参数估值求偏导数项，因此，在加权最小二乘解的基础上删去观测值的权对参数估值求偏导数项(即式(10)右边第二项)得到了参数估值的二阶无偏估计，对应的偏差改正加权最小二乘解为
(13)
式(13)中观测值的权是参数估值的非线性函数，而参数估值又是观测值的权的非线性函数，因此式(13)也为一个迭代解，其迭代格式为
(14)
式中，
由文献[5]可得偏差改正加权最小二乘迭代解的单位权中误差和一阶协方差阵分别为
(15)
(16)
由式(16)可以看出，文献[5]是通过参数估值的迭代公式和协方差传播率来获取参数估值的精度信息，由于忽略了求偏导数的二阶及其高阶项，故只能达到一阶精度。
综上所述，本文将文献[5]中加乘性混合误差模型的偏差改正加权最小二乘迭代解法总结为算法1，具体步骤如下。
(1) 输入观测值y，系数矩阵A，计算最小二乘估值作为迭代初值
(17)
(2) 计算第i次迭代观测值的权P(y)i
(18)
(19)
(20)
(3) 由式(14)计算加乘性混合误差模型的偏差改正加权最小二乘迭代解。
(4) 重复步骤(2)至步骤(3)，当满足‖
(5) 由式(15)和式(16)分别计算偏差改正加权最小二乘迭代解的单位权中误差和一阶协方差阵。
2 加乘性混合误差模型精度评定的SUT法由算法1可以看出，偏差改正加权最小二乘迭代解法中每次迭代的参数估值与观测值的函数关系相同，因此在经过第i次迭代后，最终得到的参数估值的迭代过程可以表示为
(21)
式中，
式(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=[ΔyTb1ΔyΔyTb2Δy…ΔyTbtΔ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)
式中，
(2) 将采样点集{χi}(i=0, 1, …, 2n)代入函数Y=F(y)中，求得对应的解集
(29)
(3) 加权计算因变量Y的均值E(Y)SUT和协方差阵D(Y)SUT
(30)
(31)
式中，Wim和Wic表示的是样本点集vi(i=0, 1, …, 2n)的权值，具体可表示为[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阶海森矩阵；RmSUT和RcSUT分别代表四阶及以上的高阶项；vi=
由文献[24]可得tr(M1M2M3M4)=vec(M3T)(M2T⊗M4)vec(M1)，tr(M1M2)tr(M1M2)=tr(M1M2M1M2)。因此，式(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求得的偏差改正加权最小二乘迭代解
(2) 将2n+1组观测值{χpy}(p=0, 1, …, 2n)依次代入式(37)(即算法1)求得对应的解集
(37)
(3) 由式(30)和式(31)加权计算参数估值
(4) 由式(38)和式(39)计算参数估值的单位权中误差和标准差
(38)
(39)
由上述步骤可以看出，SUT法能够解决复杂或难以实现的求导计算，具有操作简单、易于实现等特点，且利用SUT法求解加乘性混合误差模型能够有效避免复杂的求导运算，所求得的参数估值及其协方差阵均能达到二阶精度。
3 算例与分析为验证本文提出方法在加乘性混合误差模型中精度评定方面的性能及探讨在大地测量领域中的应用，本文共设计了3个算例：通过算例1来展示本文方法在GPS高程拟合中的应用价值，算例2和算例3来展示本文方法在数字地面模型(DTM)中的应用价值，并验证本文方法在精度评定方面的优势。分别使用最小二乘解法、加权最小二乘迭代解法、文献[5]中的偏差改正加权最小二乘迭代解法(即算法1)及加乘性混合误差模型精度评定的SUT法(即算法2)进行对比解算以验证本文方法的可行性及优势，4种方案对应的方法分别列于表 1中。
| 方案 | 方法 |
| 方案1 | 最小二乘解法(LS) |
| 方案2 | 加权最小二乘迭代解法(WLS) |
| 方案3 | 算法1(bcWLS) |
| 方案4 | 算法2(SUT) |
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=0.3 m，随机误差由Matlab中的mvnrnd函数生成。表 2列出了地面点的高程模拟值，为了说明高程值受加乘性混合误差的影响程度，将不含随机误差时的高程真值及受加乘性混合误差干扰的高程值绘于图 1中。
| 点号 | 高程/m | 点号 | 高程/m | |
| 1 | 10.289 1 | 17 | 24.211 6 | |
| 2 | 10.791 2 | 18 | 28.222 4 | |
| 3 | 10.494 8 | 19 | 27.031 0 | |
| 4 | 10.663 6 | 20 | 33.351 9 | |
| 5 | 10.416 3 | 21 | 34.797 0 | |
| 6 | 10.202 8 | 22 | 40.181 8 | |
| 7 | 13.133 5 | 23 | 40.861 0 | |
| 8 | 11.937 4 | 24 | 46.468 6 | |
| 9 | 13.050 2 | 25 | 53.355 2 | |
| 10 | 15.946 3 | 26 | 54.275 3 | |
| 11 | 15.799 8 | 27 | 61.193 8 | |
| 12 | 16.668 9 | 28 | 66.493 2 | |
| 13 | 18.445 3 | 29 | 70.327 9 | |
| 14 | 17.754 6 | 30 | 73.176 6 | |
| 15 | 20.055 5 | 31 | 83.162 3 | |
| 16 | 20.196 6 |
|
| 图 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的结果可知，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表示未知待求参数；函数fi(·)(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=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 |
|
|
通过对比图 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表示观测值个数；htrue, i表示获取的DTM模型高程真值；cmod, 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法，进一步表明本文方法在实测数据处理中的优势。
|
4 结论
本文首先阐述了加乘性混合误差模型的统计性质及其误差结构，然后分析了已有加乘性混合误差模型的参数估计方法能达到二阶无偏，但精度评定方法只能达到一阶精度，且无法通过近似函数法来获取参数估值的二阶精度信息，最后通过结合一种无须求导、无须了解非线性函数构成的SUT法对加乘性混合误差模型进行精度评定，并设计了具体的算法流程。
通过对本文算例中的各方法对比分析得出，当乘性误差较小时，本文算法求得的参数估值的二阶精度信息略优于已有方法求得的参数估值的一阶精度信息；而当乘性误差较大时，此时通过SUT法求得的参数估值的二阶精度信息明显优于已有方法；此外，将SUT法应用于实测DTM模型中，可以得到比已有方法更小的均方根误差结果，且参数估值的精度更高，表明将SUT方法用于解决加乘性混合误差模型的精度评定问题是可行且有效的。
本文仅是在理论方面对加乘性混合误差模型精度评定问题进行探讨，并从理论上针对加乘性混合误差模型已有解法的不足给出了一种新的算法，为今后研究LiDAR数据、SAR图像等[44-45]实际观测数据的乘性误差特性，以及结合地空观测的实践与应用提供坚实的理论基础，而如何将本文方法与方差分量估计法[46]结合进而切实运用于实际中是接下来要研究的内容。
致谢: 感谢许光煜和邹传义提出的宝贵建议和帮助。
| [1] |
江勇. 最小二乘及其扩展方法在测绘中的应用[D]. 淮南: 安徽理工大学, 2016. JIANG Yong. Application of least square method and its extension method in surveying and mapping[D]. Huainan: Anhui University of Science and Technology, 2016. |
| [2] |
TEUNISSEN P J G, KNICKMEYER E H. Nonlinearity and least-squares[J]. CISM Journal ASCGC, 1988, 42(4): 321-330. DOI:10.1139/geomat-1988-0027 |
| [3] |
师芸, 徐培亮, 彭军还. 乘性误差模型平差理论研究进展概述[J]. 工程勘察, 2014, 42(6): 60-66. SHI Yun, XU Peiliang, PENG Junhuan. Multiplicative error models: an applications-oriented review of parameter estimation methods and statistical error analysis[J]. Geotechnical Investigation & Surveying, 2014, 42(6): 60-66. |
| [4] |
SHI Yun, XU Peiliang, PENG Junhuan. An overview of adjustment methods for mixed additive and multiplicative random error models[C]//Proceedings of 2015 International Association Ⅷ Hotine-Marussi Symposium on Mathematical Geodesy. Berlin, Germany: Springer, 2015, 142: 283-290.
|
| [5] |
XU Peiliang, SHI Yun, PENG Junhuan, et al. Adjustment of geodetic measurements with mixed multiplicative and additive random errors[J]. Journal of Geodesy, 2013, 87(7): 629-643. DOI:10.1007/s00190-013-0635-2 |
| [6] |
SHI Yun, XU Peiliang. Adjustment of measurements with multiplicative random errors and trends[J]. IEEE Geoscience and Remote Sensing Letters, 2021, 18(11): 1916-1920. DOI:10.1109/LGRS.2020.3010827 |
| [7] |
WANG Leyang, CHEN Tao. Virtual observation iteration solution and A-optimal design method for ill-posed mixed additive and multiplicative random error model in geodetic measurement[J]. Journal of Surveying Engineering, 2021, 147(4): 04021016. DOI:10.1061/(ASCE)SU.1943-5428.0000363 |
| [8] |
XU Peiliang, SHIMADA Seiichi. Least squares parameter estimation in multiplicative noise models[J]. Communications in Statistics - Simulation and Computation, 2000, 29(1): 83-96. DOI:10.1080/03610910008813603 |
| [9] |
师芸. 加乘性混合误差模型参数估计方法及其应用[J]. 武汉大学学报(信息科学版), 2014, 39(9): 1033-1037. SHI Yun. Least squares parameter estimation in additive/multiplicative error models for use in geodesy[J]. Geomatics and Information Science of Wuhan University, 2014, 39(9): 1033-1037. |
| [10] |
JULIER S J, UHLMANN J K, DURRANT-WHYTE H F. A new approach for filtering nonlinear systems[C]//Proceedings of 1995 American Control Conference. Seattle, WA, America: IEEE, 1995, 3: 1628-1632.
|
| [11] |
JULIER S J, UHLMANN J K. New extension of the Kalman filter to nonlinear systems[C]//Proceedings of 1997 International Society for Optics and Photonics Signal Processing, Sensor Fusion, and Target Recognition VI. Orlando, USA: Society of Photo Optical Instrumentation Engineers, 1997, 3068: 182-193.
|
| [12] |
WAN E A, MERWE R V D. The unscented Kalman filter[M]. New York, USA: Wiley, 2001.
|
| [13] |
MERWE R V D. Sigma-point Kalman filters for probabilistic inference in dynamic state-space models[D]. Oregon: Oregon Health & Science University, 2004.
|
| [14] |
GUSTAFSSON F, HENDEBY G. Some relations between extended and unscented Kalman filters[J]. IEEE Transactions on Signal Processing, 2012, 60(2): 545-555. DOI:10.1109/TSP.2011.2172431 |
| [15] |
MENEGAZ H M T, ISHIHARA J Y, BORGES G A, et al. A systematization of the unscented Kalman filter theory[J]. IEEE Transactions on Automatic Control, 2015, 60(10): 2583-2598. DOI:10.1109/TAC.2015.2404511 |
| [16] |
WANG Leyang, ZHAO Yingwen. Scaled unscented transformation for nonlinear error propagation: accuracy, sensitivity and applications[J]. Journal of Surveying Engineering, 2018, 144(1): 04017022. DOI:10.1061/(ASCE)SU.1943-5428.0000243 |
| [17] |
WANG Leyang, ZHAO Yingwen. Unscented transformation with scaled symmetric sampling strategy for precision estimation of total least squares[J]. Studia Geophysica et Geodaetica, 2017, 61(3): 385-411. DOI:10.1007/s11200-016-1113-0 |
| [18] |
王乐洋, 丁锐, 吴璐璐. SUT法偏差改正的Partial EIV模型方差分量估计及其精度评定[J]. 大地测量与地球动力学, 2019, 39(7): 711-716. WANG Leyang, DING Rui, WU Lulu. Partial EIV model variance component estimation and accuracy evaluation of SUT method by deviation correction[J]. Journal of Geodesy and Geodynamics, 2019, 39(7): 711-716. |
| [19] |
WANG Leyang, ZHAO Yingwen, ZOU Chuanyi, et al. Scaled unscented transformation method and adaptive Monte Carlo method used for effects of fault parameters estimation on green function matrix in slip distribution inversion[J]. Geodesy and Geodynamics, 2020, 11(5): 328-337. DOI:10.1016/j.geog.2020.05.002 |
| [20] |
WANG Leyang, DING Rui. Inversion and precision estimation of earthquake fault parameters based on scaled unscented transformation and hybrid PSO/simplex algorithm with GPS measurement data[J]. Measurement, 2020, 153: 107422. DOI:10.1016/j.measurement.2019.107422 |
| [21] |
武汉大学测绘学院测量平差学科组. 误差理论与测量平差基础[M]. 3版. 武汉: 武汉大学出版社, 2014. Surveying Adjustment Group of School of Geodesy and Geomatics, Wuhan University. Error theory and foundation of surveying adjustment[M]. 3rd ed. Wuhan: Wuhan University Press, 2014. |
| [22] |
MAGNUS J R, NEUDECKER H. Matrix differential calculus with applications in statistics and econometrics[M]. New York, USA: Wiley, 1999.
|
| [23] |
徐培亮. 非线性函数的协方差传播公式[J]. 武汉测绘科技大学学报, 1986, 11(2): 92-99. XU Peiliang. Variance-covariance propagation for a nonlinear function[J]. Geomatics and Information Science of Wuhan University, 1986, 11(2): 92-99. |
| [24] |
LUTKEPOHL H. Handbook of matrices[M]. New York, USA: Wiley, 1996.
|
| [25] |
JULIER S J, UHLMANN J K. A general method for approximating nonlinear transformations of probability distributions[M]. Oxford, UK: Oxford University Press, 1996.
|
| [26] |
YORK D, EVENSEN N M, MARTINEZ M L, et al. Unified equations for the slope, intercept, and standard errors of the best straight line[J]. American Journal of Physics, 2004, 72(3): 367-375. DOI:10.1119/1.1632486 |
| [27] |
SHEN Yunzhong, LI Bofeng, CHEN Yi. An iteration solution of weighted total least-squares adjustment[J]. Journal of Geodesy, 2011, 85(4): 229-238. DOI:10.1007/s00190-010-0431-1 |
| [28] |
王乐洋, 邹传义. 乘性误差模型参数估计及精度评定的Sterling插值方法[J]. 武汉大学学报(信息科学版), 2022, 47(2): 238-244. WANG Leyang, ZOU Chuanyi. Sterling interpolation method for parameter estimation and precision estimation in multiplicative error model[J]. Geomatics and Information Science of Wuhan University, 2022, 47(2): 238-244. |
| [29] |
王乐洋, 陈涛, 邹传义. 病态乘性误差模型的加权最小二乘正则化迭代解法及精度评定[J]. 测绘学报, 2021, 50(5): 589-599. WANG Leyang, CHEN Tao, ZOU Chuanyi. Weighted least squares regularization iteration solution and precision estimation for ill-posed multiplicative error model[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(5): 589-599. DOI:10.11947/j.AGCS.2021.20200126 |
| [30] |
陈杨. 乘性误差随机模型的粗差探测[D]. 西安: 西安科技大学, 2017. CHEN Yang. Outlier detection on random model of multiplicative error[D]. Xi'an: Xi'an University of Science and Technology, 2017. |
| [31] |
BOLKAS D, FOTOPOULOS G, BRAUN A, et al. Assessing digital elevation model uncertainty using GPS survey data[J]. Journal of Surveying Engineering, 2016, 142(3): 04016001. DOI:10.1061/(ASCE)SU.1943-5428.0000169 |
| [32] |
李志林, 朱庆, 谢潇. 数字高程模型[M]. 3版. 北京: 科学出版社, 2017. LI Zhilin, ZHU Qing, XIE Xiao. Digital elevation model[M]. 3rd ed. Beijing: Science Press, 2017. |
| [33] |
LIU Xiaoye. Airborne LiDAR for DEM generation: some critical issues[J]. Progress in Physical Geography, 2008, 32(1): 31-49. DOI:10.1177/0309133308089496 |
| [34] |
HLADIK C, ALBER M. Accuracy assessment and correction of a LiDAR-derived salt marsh digital elevation model[J]. Remote Sensing of Environment, 2012, 121(138): 224-235. |
| [35] |
AKIMA H. Algorithm 760: rectangular-grid-data surface fitting that has the accuracy of a bicubic polynomial[J]. Acm Transactions on Mathematical Software, 1996, 22(3): 357-361. DOI:10.1145/232826.232854 |
| [36] |
ZIMMERMAN D, PAVLIK C, RUGGLES A, et al. An experimental comparison of ordinary and universal Kriging and inverse distance weighting[J]. Mathematical Geology, 1999, 31(4): 375-390. DOI:10.1023/A:1007586507433 |
| [37] |
KIDNER D B. Higher-order interpolation of regular grid digital elevation models[J]. International Journal of Remote Sensing, 2003, 24(14): 2981-2987. DOI:10.1080/0143116031000086835 |
| [38] |
HILL C A, HARRIS M, RIDLEY K D, et al. LiDAR frequency modulation vibrometry in the presence of speckle[J]. Applied Optics, 2003, 42(6): 1091-1100. DOI:10.1364/AO.42.001091 |
| [39] |
KOBLER A, PFEIFER N, OGRINC P, et al. Repetitive interpolation: a robust algorithm for DTM generation from aerial laser scanner data in forested terrain[J]. Remote Sensing of Environment, 2007, 108(1): 9-23. DOI:10.1016/j.rse.2006.10.013 |
| [40] |
LEIGH C L, KIDNER D B, THOMAS M C. The use of LiDAR in digital surface modelling: issues and errors[J]. Transactions in GIS, 2010, 13(4): 345-361. |
| [41] |
SHI Yun, XU Peiliang, PENG Junhuan, et al. Adjustment of measurements with multiplicative errors: error analysis, estimates of the variance of unit weight, and effect on volume estimation from LiDAR-type digital elevation models[J]. Sensors, 2014, 14(1): 1249-1266. |
| [42] |
王乐洋, 许才军, 鲁铁定. 病态加权总体最小二乘平差的岭估计解法[J]. 武汉大学学报(信息科学版), 2010, 35(11): 1346-1350. WANG Leyang, XU Caijun, LU Tieding. Ridge estimation method in ill-posed weighted total least squares adjustment[J]. Geomatics and Information Science of Wuhan University, 2010, 35(11): 1346-1350. |
| [43] |
WANG Leyang, CHEN Tao. Ridge estimation iterative solution of ill-posed mixed additive and multiplicative random error model with equality constraints[J]. Geodesy and Geodynamics, 2021, 12(5): 336-346. |
| [44] |
HE Yufang, ZHU Wu, LEI Yang, et al. A comparative study of ionospheric correction on SAR interferometry: a case study of L'Aquila earthquake[J]. Journal of Geodesy and Geoinformation Science, 2022, 5(1): 5-13. |
| [45] |
ZHU Wu, LEI Yang, SUN Quan. Detection, estimation and compensation of ionospheric effect on SAR interferometry using azimuth shift[J]. Journal of Geodesy and Geoinformation Science, 2022, 5(1): 14-24. |
| [46] |
王乐洋, 温贵森. 偏差改正的Partial EIV模型方差分量估计[J]. 测绘学报, 2019, 48(4): 412-421. WANG Leyang, WEN Guisen. Bias-corrected variance components estimation of partial EIV model[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(4): 412-421. DOI:10.11947/j.AGCS.2019.20170693 |