王乐洋1 ，陈涛1，邹传义1，2 　　

　　

　　1.东华理工学院测绘工程学院，江西南昌330013； 　　

　　

　　2.湖北武汉测绘学院430079 　　

　　

　　基金项目：国家自然科学基金(41874001；41664001);东华理工大学研究生创新项目(DHYC-202020) 　　

　　

　　摘要:针对乘性误差模型的病态问题，引入Tikhonov正则化方法，推导出病态乘性误差模型的加权最小二乘正则化解。考虑到加权最小二乘正则化方法求解病态乘性误差模型时，参数估计与观测值之间存在复杂的非线性关系，本文采用无需求导数来评估病态乘性误差模型的精度，利用尺度无迹变换(SUT)方法，通过加权计算非线性函数的均值和均方误差矩阵的比例对称采样。仿真和实例结果表明，本文提出的加权最小二乘正则化迭代解法能有效削弱病态模型，基于SUT方法的精度评估方法能获得比现有方法更合理的精度信息，具有较强的适用性。 　　

　　

　　关键词：病态乘性误差模型的精度评估吉洪诺夫正则化L曲线法SUT法 　　

　　

　　王乐洋、陈涛、邹传义。不适定乘性误差模型的加权最小二乘正则化迭代求解及精度估计。大地测量学与制图学学报，2021，50(5): 589-599。1947年11月10日，澳大利亚地球物理学会。63636363667 　　

　　

　　在大地测量数据处理领域，如果观测误差随观测值的大小或物理信号的强弱而变化，这种误差就是乘性误差1。如现代观测手段中合成孔径雷达(SAR)观测值的随机误差为乘性误差2-3，电子测距(EDM)和全球定位系统(GPS)的观测误差为乘性误差或加性混合误差4-5。目前在大地测量学领域，关于乘性误差模型的参数估计和精度评价的研究成果相对较少，6-9，没有文献研究乘性误差模型中的病态问题，因此如何估计病态乘性误差模型的参数和评价其精度是一个值得研究的问题。 　　

　　

　　测绘数据处理中普遍存在弊病。模型有病会引起参数解的不稳定，甚至严重偏离真值10。现有的总体最小二乘11、大地反演12和控制网平差13的不适定问题都是基于加性误差模型的，而对于乘性误差模型的不适定问题还没有发现研究成果。文献6首先提出了乘性误差模型的最小二乘法。文献7在文献6的基础上，总结了乘性误差模型的最小二乘法、加权最小二乘法和偏差修正加权最小二乘法，并推导了参数估计的精度评定公式。如果在不考虑病态系数矩阵的情况下，直接用文献7中的三种方法处理病态乘性误差模型，参数估计会有偏差且不稳定。4为了避免系数矩阵的病态，模拟点集中在数字地形模型的顶点(DTM)，其他点均匀分布，只能用于处理模拟数据，不能应用于真实观测数据。系数矩阵病态会导致系数矩阵列向量之间的多重共线性，从而导致法方程的条件数过多，解不稳定。目前处理病态问题的方法主要有Tikhonov正则化方法14-16和虚拟观测解方法17。确定正则化参数的方法主要有广义交叉检验法18、岭迹法19和L曲线法20-22。广义交叉验证法难以获得最优解，适用性较弱；岭迹法虽然计算简单，但具有一定的主观性；l曲线法适用性强，可以获得合理的正则化参数。它广泛应用于选择 　　

　　

　　通过分析发现，用加权最小二乘正则化方法求解病态乘性误差模型时，参数估计是观测值的非线性函数，而观测值的协变量矩阵是参数估计的非线性函数。而且在逐步迭代的过程中，参数估计的每一步都是随机的，使得最终的参数估计与观测值之间存在复杂的非线性关系，只能通过数值迭代求解，无法得到解析解。为此，引入了不受需求引导、不需要获得非线性函数具体表达式的SUT方法23-24来评估病态乘性误差模型的精度。 　　

　　

　　综上所述，鉴于现有文献尚未考虑病态乘性误差模型，本文首先利用Tikhonov正则化方法推导出病态乘性误差模型的加权最小二乘正则化解，然后利用SUT方法评估病态乘性误差模型的精度。最后，通过两个仿真数据实例和一个真实数据实例验证了该方法的适用性和优越性。 　　

　　

　　给定一组具有乘性误差的观测值，相应的乘性误差模型可以表示为7。 　　

ic/img.php?k=srt项目的可行性分析怎么写,srt项目怎么做2.jpg">(1)

　　

随机模型为

　　

进一步可将式(1)改写为如下形式

　　

本文研究的f()是β的线性函数形式，即fi(β)=aiTβ，其中aiT为n×t维向量，则式(3)可以写为

　　

由式(4)可知，在乘性误差模型中，观测值精度与系数矩阵和参数估值的乘积有关：信号Aβ越强，观测值y的误差越大；反之，信号Aβ越弱，误差越小。这明显区别于传统的加性误差模型，即观测值y的误差与系数矩阵和参数估值无关。目前，已有关于乘性误差模型的研究均不考虑系数矩阵A病态<7-8>，然而在大地测量实际数据处理中，系数矩阵A很容易出现病态，此时会引起法方程条件数过大，观测量微小的扰动将会引起参数较大的变化，导致解不稳定。例如，文献<25>针对SAR型噪声的特性采用乘性误差模型进行影像去噪，去噪效果达到70%~93%。此时，若系数矩阵A病态，会引起法方程条件数过大，导致解不稳定，进而影响影像去噪。因此需要针对病态乘性误差模型进行更深入的研究来丰富现代大地测量数据处理理论。

　　

为了更方便地进行公式推导，将式(4)改写为如下形式

　　

根据协因数传播率，得到式(6)中e和y的协因数阵为

　　

根据最小二乘准则求解式(6)，得

　　

β的最小二乘估值为

　　

根据前文中的分析可知，当法方程N=ATA病态时，由式(9)求得的最小二乘解不可靠。因此，引入正则化因子α，构建如下病态乘性误差模型的正则化准则

　　

对式(10)中的β求偏导，并令其为0，可得

　　

将式(6)代入式(11)，可得病态乘性误差模型的加权最小二乘正则化解

　　

式中，较是一个极小量，可在计算过程中舍弃。

　　

进而可以得到

　　

式中，为加权最小二乘正则化解；α为正则化参数；In为单位阵。

　　

与文献<7>中的偏差改正加权最小二乘解相比，式(13)右端法方程求逆部分增加了αIn项，有效抑制了法方程的病态性，进而可以得到可靠、稳定的参数解。

　　

由式(10)可知，和都是关于正则化参数α的函数，L曲线法确定病态乘性误差模型参数估计正则化参数的主要步骤是<20-21>：以为横坐标、为纵坐标作图，得到多组坐标点，通过曲线拟合将这些坐标点绘成一条类似“L”形状的曲线，曲率最大的那个点所对应的α值作为所求的正则化参数值α。

　　

在确定本文病态乘性误差模型中的正则化参数时，由于观测值的协因数阵为参数估值的非线性函数，因此观测值的权阵也为参数估值的非线性函数，这使得正则化参数值产生动态变化。

　　

根据文献<8>，病态乘性误差模型的单位权方差估值可表示为

　　

由式(13)可知，加权最小二乘正则化解与观测值是线性的关系，由协因数传播率可得的协因数阵为

　　

由偏差的定义<26>，对式(13)两边取期望，可得参数估值的偏差为

　　

利用矩阵反演公式<27>，将正则化解作为参数真值<28-29>，进一步可将式(16)表示为

　　

病态问题中，由于产生了偏差项，此时的均方误差不等于协方差不能作为一个良好的精度评判标准。虽然最小二乘估值是无偏，但已不是最优值<7>。根据均方误差的定义<30-31>，可知均方误差的表达式为

　　

在求解式(13)时，协因数阵Qe是参数估值的非线性函数，而求得的参数估值又是协因数阵Qe的非线性函数，只能使用数值逼近的方法进行迭代求得数值解，而无法得到参数估值的解析解。因此，在由式(13)求解参数估值时，需利用上一次计算得到的参数估值更新协因数阵Qe和正则化参数α，进而求解参数估值，迭代循环直至前后两次参数估值之差小于某一特定的阈值ε，则迭代终止；本文取ε为10-10。

　　

式(13)的迭代公式如下所示

　　

式中，计算得到；αi为第i次迭代的正则化参数。

　　

综上所述，本文将病态乘性误差模型的加权最小二乘正则化法迭代解法命名为算法1，具体步骤如下：

　　

由算法1可以看出，病态乘性误差模型的加权最小二乘正则化迭代解法中，每一步迭代参数估值都是观测值的函数，而观测值的协因数阵又是参数估值的函数，使得最终的参数估值与观测值为一个复杂的非线性函数。本文将这种非线性关系表示为<32-33>

　　

第1步参数迭代估值可以表示为<32-33>

　　

最终的参数估值表达式为<32-33>

　　

式(23)表明，病态乘性误差模型的加权最小二乘正则化法的迭代过程使得参数估值和观测值之间的非线性关系非常复杂，且迭代的过程使参数的每一步的估值都具有随机性。为了避免复杂的公式运算，本文在病态乘性误差模型加权最小二乘正则化法中引入SUT采样法，旨在获得更精确的参数估值和精度信息。

　　

本文将病态乘性误差模型精度评定的SUT法命名为算法2，步骤如下：

　　

本文首先对病态乘性误差模型加权最小二乘正则化迭代解法进行改化，将病态乘性误差模型中参数估值与观测值表示为嵌套函数的形式，然后引入无需求导、无需了解非线性函数构造的SUT法求解病态乘性误差模型中参数估值的均值、均方误差矩阵和单位权方差。

　　

为了比较不同方法下参数估值和精度信息的差异，分别使用文献<7>中的未考虑病态性的三种方法即最小二乘法、加权最小二乘法和偏差改正加权最小二乘法，以及本文提出算法1和算法2进行解算。5种方案及对应的方法见表 1。表 1 5种方案及对应的方法

　　

Tab. 1 Five schemes and corresponding methods

　　

　　

　　

为了验证本文公式推导以及算法的有效性，算例1根据文献<34>改编得到，通过利用GPS测量某一区域道路中心线地面点的高程，假定该区域道路中心线地面点高程真值与某一坐标原点的距离除以100符合以下函数模型

　　

本算例中x的取值范围为0~300 m，在x的取值范围内等间距选取31个点，根据式(29)计算出高程点的真值，参数真值为式(29)右侧的5个系数，分别为10，2，1，1和0.5。

　　

假设地面点高程真值被乘性误差干扰，其中乘性误差向量εm相互独立，且服从均值为0、标准差为0.1的正态分布，即εm~0, 0.12×I31，相应的乘性误差模型的观测方程为

　　

由式(30)得到一组地面高程点的数据模拟值见表 2。为了说明高程点受乘性误差的影响程度，本文将地面高程点受乘性误差干扰前后的数据绘于图 1中。

　　

表 2 数据模拟值

　　

Tab. 2 The simulated data m

　　

　　

　　

图 1 地面高程点受乘性误差干扰前后分布

　　

Fig. 1 GPS elevation points before and after disturbed by multiplicative error

　　

　　

本算例选取的单位权中误差σ0为0.3。使用表 1中的5种方案计算，将参数真值、5种方案计算的参数估值、参数估值与真值的2范数及单位权中误差估值列于表 3；5种方案计算的参数估值的均方根误差列于表 4。

　　

表 3 参数估值、参数估值与真值之间的2范数和单位权中误差

　　

Tab. 3 Parameter estimation, 2-norm between the parameter estimation and true value and the unit weight error

　　

　　

　　

　　

表 4 参数估值的均方根误差

　　

Tab. 4 Root mean square error of parameter estimation

　　

　　

　　

由上文分析，本文算法中的正则化参数值会产生动态变化，正则化参数α随迭代次数的变化如图 2所示。经计算，本算例法矩阵N=的条件数为9.896 2×104，具有病态性，使用本文方法求得的法矩阵条件数随迭代次数的变化如图 3所示。

　　

图 2 正则化参数α随迭代次数的变化

　　

Fig. 2 The regularization parameter with the number of iterations

　　

　　

图 3 法矩阵条件数随迭代次数的变化

　　

Fig. 3 The condition number of the normal matrix with the number of iterations

　　

由图 1可知，由于地面高程点受到乘性误差干扰，导致点位产生了严重偏离。由表 3中参数估值与真值的2范数可知，LS方法求得的2范数结果最大；WLS方法由于考虑了观测值的权得到的结果优于LS法，其2范数为6.058 4；bcWLS方法对WLS方法进行了偏差改正，2范数较WLS方法减小了0.119 4，为5.939 0，但由于未考虑模型的病态性，仍与真值偏离较大。本文算法1和算法2求得的2范数分别为0.978 4和0.833 7，参数估值更接近真值，说明本文算法对于降低病态性有一定的效果。其中，由于算法2考虑了正则化解法逐步迭代过程中每一步的随机性和非线性迭代的过程对参数估值的影响，求得的2范数较算法1减小了0.144 7，表明了本文使用SUT法的可行性。由表 3中的单位权中误差可知，由于文献<7>中的方法未考虑模型的病态性，导致求得的单位权中误差偏离真值，而本文算法考虑了这些影响，求得的单位权中误差更接近真值，进一步表明了本文算法的优势。

　　

由表 4可得，文献<7>中3种方法求得的参数估值的均方根误差较大，这是因为法矩阵中存在接近于0的特征值，导致求逆不稳定，进而使得均方根误差结果较大，而本文算法1和算法2求得的均方根误差均比文献<7>中的3种方法小，说明了结果的可靠性。算法2考虑到迭代过程对参数估值的影响，使用SUT法进行精度评定，求得的参数估值的均方根误差最小，其中的均方根误差比算法1小了0.298 0，这也验证了本文算法的适用性和优势。

　　

算例2为一个DTM模型的算例，本文主要考虑模型为病态的情况，因此模拟一个病态DTM模型。文献<4>在模拟DTM模型时，为防止法矩阵病态情况出现，在模拟的DTM模型的每个峰值周围采样更多点，并使峰值点之间远离，虽然可以在一定程度上降低病态性，但仅可用来处理模拟数据，无法适用于真实的观测数据。本文依据文献<4>的思路，使用插值法来模拟DTM模型，插值法是模拟DTM模型的主要方法之一<35-36>，本文模拟的DTM模型由下面的插值函数生成

　　

(31)

　　

函数fi(x, y)(i=1, 2, 3, 4)分别为如下

　　

(32)

　　

(33)

　　

(34)

　　

(35)

　　

模拟的DTM模型如图 4所示，相应的乘性误差模型的观测方程为

　　

(36)

　　

图 4 模拟的DTM模型

　　

Fig. 4 Simulated DTM model

　　

在本文中，假设误差向量εm相互独立、且服从均值为0，标准差为0.1的正态分布，即εm~0, 0.12×I1681。为了说明DTM模型受乘性误差的影响程度，本文将不含误差的观测值绘于图 4中，受乘性误差干扰的观测值绘于图 5中。

　　

图 5 受乘性误差干扰的DTM模型

　　

Fig. 5 DTM model disturbed by multiplicative error

　　

本算例选取的单位权中误差σ0为0.3。使用表 1中的5种方案计算，将参数真值、5种方案计算的参数估值、参数估值与真值的2范数及单位权中误差估值列于表 5；5种方案计算的参数估值的均方根误差列于表 6。正则化参数α随迭代次数的变化如图 6所示。经计算，法矩阵N=的条件数为5.190 5×104，具有病态性，使用本文方法求得的法矩阵条件数随迭代次数的变化如图 7所示。

　　

表 5 参数估值、参数估值与真值之间的2范数和单位权中误差

　　

Tab. 5 Parameter estimation, 2-norm between parameter estimation and true value and the unit weight error

　　

　　

　　

表 6 参数估值的均方根误差

　　

Tab. 6 Root mean square error of parameter estimation

　　

　　

　　

图 6 正则化参数α随迭代次数的变化

　　

Fig. 6 The regularization parameter with the number of iterations

　　

　　

图 7 法矩阵条件数随迭代次数的变化

　　

Fig. 7 The condition number of the normal matrix with the number of iterations

　　

由图 4和图 5可以发现，尽管本文模拟的DTM模型加入的乘性误差标准差仅为0.1，但对高程产生了较大的影响，此时法矩阵的条件数为5.190 5×104，严重病态。因此，针对病态乘性误差模型数据处理理论需要进行更深入的研究。

　　

由表 5可得，由于未考虑模型的病态性，文献<7>中的方法求得的参数估值严重偏离真值，而本文算法求得的结果更接近真值，参数估值的2范数分别为0.817 7和0.770 5，比bcWLS方法小了10.031 8和10.079 0，进一步说明本文算法对于降低病态性有一定的效果。其中，本文算法2考虑到非线性迭代的过程对参数估值的影响，求得的2范数最小，这与算例1得到的结果一致。由表 5中的单位权中误差可知，LS方法求得的单位权中误差0.700 6，结果严重偏离真值；WLS方法和bcWLS方法由于考虑了观测值的权，求得的单位权中误差比LS方法更接近真值，为0.293 0和0.294 4；本文算法考虑了病态性对单位权中误差的影响，求得的单位权中误差最接近真值。

　　

由表 6可得，本文使用正则化方法得到的均方根误差远小于文献<7>中的方法，表明了结果的可靠性。在本文算法中，算法1求得的参数估值、、的精度信息有了较大改善，其均方根误差远远小于bcWLS方法，但是的均方根误差只比bcWLS方法减少了0.148 9，这可能是由于正则化迭代过程对参数估值带来的影响。本文算法2将算法1正则化迭代过程中参数估值与观测值的关系视为一个非线性嵌套函数，考虑到非线性迭代过程对参数估值的影响，通过SUT法进行精度评定，得到参数估值的均方根误差最小，且的均方根误差比算法1减少了0.355 3，这也验证了本文方法是可行且有效的。

　　

为进一步说明本文方法的可行性和适用性，算例3为一个病态数字高程模型真实数据算例。数字高程模型是一种承载地面高程信息的空间数据模型，是DTM模型的一个分支，本文主要利用SRTM(shuttle radar topography mission)提供的免费DEM数据集，通过在网站http://srtm.csi.cgiar.org/srtmdata/下载了2304个点的三维坐标数据来说明本文算法在实际应用中的优势。SRTM数据是由航天飞机雷达测量生成的LiDAR影像数据<37>，且LiDAR数据的噪声已被证明具有乘性噪声性质<38>，因此本文含有6个未知参数的数字高程DEM乘性误差模型如下

　　

由于SRTM数据均是在等精度测量中获取，因此在定权时假设未知的单位权中误差和乘性误差相等，即乘性误差的权阵为单位阵。本文获取的受乘性误差干扰的数字高程模型如图 8所示，分别使用表 1中的5种方案计算，WLS方法和bcWLS方法由于受病态性的影响，系数矩阵奇异，参数估值不收敛；LS方法、本文算法1和2方法拟合的数字高程模型如图 9、图 10和图 11所示。将3种收敛方案计算的参数估值和单位权中误差估值列于表 7，参数估值的均方根误差列于表 8。

　　

图 8 受乘性误差干扰的数字高程模型

　　

Fig. 8 Digital elevation model disturbed by multiplicative error

　　

　　

图 9 LS方法求得的数字高程模型

　　

Fig. 9 Digital elevation model obtained by the LS method

　　

　　

图 10 RWLS方法求得的数字高程模型

　　

Fig. 10 Digital elevation model obtained by the RWLS method

　　

　　

图 11 SUTRWLS方法求得的数字高程模型

　　

Fig. 11 Digital elevation model obtained by the SUTRWLS method

　　

　　

表 7 参数估值和单位权中误差

　　

Tab. 7 Parameter estimation and the unit weight error

　　

　　

　　

表 8 参数估值的均方根误差

　　

Tab. 8 Root mean square error of parameter estimation

　　

　　

　　

由图 8和图 9可以发现，LS方法求得拟合高程的结果严重偏离原始的数字高程模型，这是由于LS方法未考虑观测值的权且模型为病态性所导致的。本文算法考虑了病态性的影响，在处理实际数据时更具优势，从图 10和图 11来看，求得拟合高程的结果能够较好地拟合图 8的高程。说明现有文献未考虑病态性的方法无法处理实际数据，因此顾及模型病态性带来的影响是十分必要的。

　　

由表 7可以看出，LS方法求得的参数估值和单位权中误差过大，不符合实际情况，进一步说明了LS方法无法处理实际数据，本文算法在处理实际数据的优势。由表 8也可以看出，LS方法求得参数估值的根均方误差最大，本文算法远小于LS方法。其中，在本文算法中，由于算法2考虑了非线性迭代过程对参数估计带来的影响，使用SUT法求得的参数估值的均方根误差小于算法1，这也验证了本文算法的适用性和优势。对比算法1和算法2在参数估值上并没有很大差别，这可能是因为参数估值较小；而参数估值的均方根误差却差别很大，这是因为算法1根据误差传播定律求得的均方根误差只能反映参数估值的一阶精度信息，而算法2采用的SUT法为二阶精度评定方法，求得参数估值的精度信息为二阶精度。

　　

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

　　

　　

初审：张艳玲

　　

《测绘学报》专刊征稿函 | 滑坡灾害监测预警

　　

　　

《测绘学报（英文版）》（JGGS）专刊征稿 | Call for Papers：空间人文与社会地理计算 (SHGSS)

　　

　　

重磅 | 第八届“测绘科学前沿技术论坛”会议通知（一号）

　　

　　

自然资源部办公厅关于印发《自然资源部2021年政务公开工作要点》的通知

　　

　　

　　

<