分享

论文推荐 | 胡川, 方兴, 赵立都:融合正交几何信息的非线性等式约束整体最小二乘平差及迭代算法

 沐沐阅览室 2020-08-08

《测绘学报》

构建与学术的桥梁        拉近与权威的距离

复制链接,关注《测绘学报》抖音!

【测绘学报的个人主页】长按复制此条消息,长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]



融合正交几何信息的非线性等式约束整体最小二乘平差及迭代算法
引文格式:胡川, 方兴, 赵立都. 融合正交几何信息的非线性等式约束整体最小二乘平差及迭代算法. 测绘学报,2020,49(7):816-823. DOI: 10.11947/j.AGCS.2020.20190112.
基金项目:国家自然科学基金(41774009);重庆市基础科学与前沿技术研究(一般)项目(cstc2017jcyjAX0102);重庆市教委科学技术研究(KJ1705132);中央高校自主科研学科交叉项目(2042018kf0230);2016、2017年重庆交通大学高层次人才科研启动项目(16JDKJC-A025;17JDKJC-A027)
阅读全文:http://xb./article/2020/1001-1595/2020-7-816.htm
全文概述
与只考虑因变量误差影响的加权最小二乘(weighted least squares,WLS)平面线形拟合模型不同,变量误差模型(error-in-variables,EIV)兼顾观测向量(因变量)和系数矩阵(自变量)的误差,平差模型更合理,理论更严密,是当前测绘领域的研究热点之一。

解算EIV模型的流行方法是整体最小二乘法[1](total least squares,TLS),其主要的计算方法有奇异值分解法(singular value decomposition,SVD)和迭代法[2-4]两种。根据是否需要对函数式进行线性化,迭代法可进一步划分为线性化法和非线性化法[5]。目前对TLS的研究主要集中在计算算法开发[5-8]、应用范围扩展[9-15]和解算方法拓展[16-21]等方面。在数学领域,将误差矩阵附有与系数矩阵结构一致性约束的TLS称作约束整体最小二乘(constrained total least squares,CTLS)[22],也常被称作结构TLS(structured total least squares,STLS),也有学者将其称作偏整体最小二乘(partial total least squares,PTLS)[7]。根据已知的先验信息,给TLS附加上参数恒等式约束,可将其扩展为附有等式约束的TLS(equality constrained total least squares,ECTLS)[23]。当约束等式只有固定量时,其被扩展为附有固定等式约束的TLS[24];当约束等式含有随机量时,其被扩展为附有随机等式约束的TLS[25]。它们都属于线性等式约束问题,有关此问题的综合讨论可参阅文献[26]。文献[25, 27]讨论了附有二阶等式约束的TLS平差问题及其解算方法,它属于非线性固定等式约束问题,不含有误差改正数。

正交距离最小二乘和加权整体最小二乘(weighted total least squares,WTLS)在拟合平面线形时都考虑了自变量的误差干扰。但是它们是两种相互独立的方法,前者从几何角度描述拟合准则,后者从代数角度描述拟合准则。几何准则的优点是能保证测量点到拟合对象的正交距离的平方和最小,代数准则的优点是能灵活处理自变量和因变量的随机信息。考虑到基于代数准则的WTLS并不能完全确保其正交距离的平方和为极小值,笔者将几何准则中的正交信息加入其中,建立一种约束方程含有误差改正数的非线性等式约束整体最小二乘平差模型(nonlinear equality constrained total least squares,NECTLS),并研究其解算算法以及在直线拟合中的应用。

1 NECTLS平差法  

常规WTLS的函数模型和随机模型可分别表示为[1-2]

 (1)

 (2)

式中,y和ey分别是n×1的观测向量和误差向量;A和EA分别是n×m的系数矩阵和误差矩阵;x和ex=vec(EA)分别是nm×1的观测值向量和误差向量;ξ是m×1的参数向量;符号E和D分别表示期望和方差;σ02表示未知的单位权方差;Qx和Qy分别是nm×nm和n×n的协因数矩阵。在常规模型中,系数矩阵没有常数元素或(和)重复元素,不具有结构特性,同时也不对参数附加额外的等式约束。

当系数矩阵Α含有重复元素(包括随机量和常数)时,ex≠vec(EA),表现出明显的结构特征。为消除其影响,对函数式(1)稍作变换,整理可得

 (3)

式中,A(x)表示A是关于x的矩阵函数;⊗表示克罗内克积;N表示nm×n的转换矩阵[7],根据向量求导法则可推得其计算式[5],即

 (4)

为能将几何准则中的正交信息融入WTLS,首先将该正交信息转换为可执行的数学表达式,然后将其作为约束方程加入WTLS,最后综合组成NECTLS平差的函数模型。根据正交几何信息,可以得到如下两个具有等价性的数学表达式,即

 (5)

 (6)

式中,带波浪线和带小冒的量分别表示预测值和估计值,误差预测值即是误差改正数;i=1, 2, …, n。

综上所述,NECTLS平差的误差方程、约束方程和平差准则可分别描述为

 (7)

 (8)

 (9)

 (10)

式中,Py=Qy-1,Px=Qx-1。

2 NECTLS的解算  

考虑到约束方程式(8)和式(9)通常具有等价性,笔者在此仅对带有约束方程式(8)的情况进行讨论。

将整体最小二乘平差看成是一个非线性参数估计问题[69],令

 (11)

 (12)

则非线性方程式(7)可被线性化为[6]

式中,Y=y-Aξ(0),Z=NT(ξ(0)⊗In)。将式(11)和式(12)代入非线性误差方程式(8),对其进行线性化,可得

 (14)

式中,

 (15)

 (16)

它们分别是n×1和n×m的矩阵。值得注意的是,只有当非线性模型接近线性时,即非线性程度较弱时,这种处理才比较合理[28]。

根据线性化后的数学方程和拉格朗日乘数法,可将求解参数最优估值的目标函数定义为

 (17)

式中,λ和μ是n×1的拉格朗日乘积因子。

将目标函数式(17)对各量求偏导,并令它们等于零,则有

 (18)

 (19)

 (20)

 (21)

 (22)

由式(18)和式(19)可分别推得预测值的计算公式,即

 (23)

 (24)

它们是关于拉格朗日乘积因子的函数。将式(23)和式(24)代入式(20)和式(21),整理后可得

 (25)

 (26)

将式(25)和式(26)代回式(23)和式(24),可得误差预测值关于ξ(0)和的函数表达式,即

 (27)

 (28)

将式(25)和式(26)再代回式(22)中,整理后可得的计算公式,即

 (29)

式中,

根据式(27)和式(28)计算的观测误差预测值,可求得误差的加权平方和,考虑到误差预测值是关于拉格朗日乘积因子的函数,故也可用式(30)来计算其值,即

 (30)

将其与多余观测数(n-m)相除,得到验后单位权方差的估计值,其计算公式为

 (31)

根据协因数传播定律,由式(12)、式(29)和式(31)可推得评定参数估计精度的计算式(32)

 (32)

3 NECTLS的计算算法  

虽然约束条件式(8)和式(9)的求解过程类似,但它们求解参数估值的计算公式略有不同。因此,尽管它们可以采用相同的计算流程,但是需要采用与各自相对应的计算公式。本节以式(8)为约束条件推求的计算公式为基础,设计NECTLS的计算算法,具体计算过程如下:

(1)根据已知观测量,给定y,x,Py,Px,N。

(2) 用x组成系数矩阵A,采用WLS估计参数初值:ξ(0)=(ATPyA)-1ATPyy,并令EA(0)=0。

(3) 按线性化公式依次计算下列各值

(4) 按下列公式计算参数修正量和误差预测值

(5) 如果,则,EA(i+1)=vec-1(Nex(i+1)),将值代回到第(3)—(4)步迭代计算;反之,则令,并分别用式(31)和式(32)计算验后单位权方差和参数估值的方差-协方差,停止迭代,输出平差结果。

4 算例及分析  

本文以直线参数估计为例,验证笔者提出的方法和算法。设待拟合直线的方程为

 (33)

式中,a是斜率;b是截距;(xi, yi)是一组观测数据,即测量点。观测值yi组成观测向量y,观测值xi组成系数矩阵A=[x  1];协因数矩阵Qx=[diag(wx1, wx2, …, wxn)-1,Qy=[diag(wy1, yx2, …, wyn)]-1。

要采用NECTLS估计直线参数,需要将测量点与拟合点的连线垂直于拟合直线这一正交几何信息转换成非线性等式约束方程,即

 (34)

 (35)

对于直线拟合而言,这两种约束具有等价性,参数估计结果完全相同。

对误差约束方程式(34)和式(35)进行线性化,可以分别得到

 (36)

 (37)

 (38)

 (39)

将它们分别代入式(15)、式(16)和式(41)、式(42),可求得C、H和K、J的值。在迭代计算时,它们都需要利用前一步的计算结果重新计算。为能从统计上和个例上对本文提出的方法进行检验,设置如下两个试验。

试验1:假设直线的斜率a=0.35,截距b=5.6,给定8个无观测误差的x坐标值(5.88,5.44,4.47,4.62,3.51,3.75,2.86,2.81),以此计算对应的y坐标值,单位为m。分别给x和y坐标,附加上期望为零,方差分别为0.01 m2和0.006 4 m2的随机误差,获得含有随机误差的坐标观测值,模拟1000次, 每次都采用WLS、WTLS(文献[2]的算法)和NECTLS进行解算。

用3种方法估计的参数值减去给定值,分析后发现:①以式(34)或式(35)为约束条件计算的参数估值相同;②就单次模拟计算而言,没有明确规律表明哪一种方法估计的参数值最接近给定值。将这些差值作直方图统计,并拟合出对应的正态分布曲线,如图 1所示。从中可以发现:①WTLS和NECTLS的标准差比WLS的略大,WTLS的比NECTLS略大;②WTLS和NECTLS的平均值比WLS的小,NECTLS的最小。

图 1 估计参数与真值差值的直方分布Fig. 1 Histogram of difference between the estimated parameters and the given value

图选项 

将3种方法每次计算的验后单位权方差相互作差,用WLS减去WTLS(WLS-WTLS),NECTLS减去WLS(NECTLS-WLS),NECTLS减去WTLS(NECTLS-WTLS),结果如图 2所示。

图 2 验后单位方差之差Fig. 2 Difference of variance of unit weight

图选项 

WLS-WTLS和NECTLS-WTLS的值始终为正值,说明WLS和NECTLS的验后单位权方差始终比WTLS的大;NECTLS-WLS的值始终为负值,说明NECTLS的验后单位权方差始终比WLS的小。因此,3种方法计算的验后单位权方差的大小关系为:WLS>NECTLS>WTLS,3种方法模拟1000次的平均值分别为1.184 07、1.013 75、0.994 19。

3种方法得到的参数估值的方差也存在明显差异,将它们作差,用NECTLS减去WLS、WTLS的结果,分别用NECTLS-WLS和NECTLS-WTLS表示,WTLS减去WLS的结果用WTLS-WLS表示,结果如图 3所示。

图 3 参数估值方差之差Fig. 3 Difference of variance of estimated value

图选项 

比较分析可以发现:①尽管NECTLS和WLS的值非常接近,但前者还是略小于后者;②NECTLS-WTLS的值始终为负,说明NECTLS计算的值比WTLS的小;③WTLS-WLS的值始终为正,说明WTLS计算的值比WLS的大。综合而言,NECTLS的值最小。

每次模拟计算完成后,计算各测量点到拟合点的距离dp和测量点到拟合直线的距离dl,并计算它们差值的平均值:∑(dp-dl)/8。如图 4所示,NECTLS的差值等于零,表明这两种距离相等;WLS和WTLS的差值大于零,表明这两种距离存在明显差异,且dp大于dl,同时WLS的差异最大。

图 4 dp减去dl的平均值Fig. 4 Mean of difference between dp and dl

图选项 

在每次模拟计算过程中,根据计算的点到直线的距离dl计算距离的平方和(正交距离平方和),并将不同方法计算的结果做差,用NECTLS减去WLS和WTLS的结果,分别用NECTLS-WLS和NECTLS-WTLS表示,如图 5所示。从图上不难发现,它们的差值始终为负值,说明NECTLS的正交距离平方和最小,与WLS和WTLS相比,NECTLS实现了正交距离平方和最小的参数估计。

图 5 正交距离平方和差值Fig. 5 Difference of ∑dl2 between NECTLS and WLS, NECTLS and WTLS

图选项 

试验2:采用整体最小二乘研究中常用的直线观测数据进行验证分析,观测数据和相应权值如文献[2]的表 1所示,它们分别取自文献[29]和文献[30]。WLS、WTLS和NECTLS(以式(34)或(35)为约束条件,解算的结果相同)的参数估计结果见表 1。对结果进行比较分析,可以发现:①3种方法估计的模型参数值明显不同;②WTLS计算的验后单位权方差最小;③NECTLS计算的参数估值的方差最小。它们与试验1的统计结果一致。

表 1 参数估计结果Tab. 1 Estimated parameters

参数方法
WLSWTLSNECTLS
â-0.610 8-0.480 5-0.576 5
6.100 15.479 95.867 9
4.293 21.483 37.133 1
0.003 90.005 00.000 8
0.179 80.129 10.015 8
-0.026 0-0.024 4-0.002 7

表选项 

表 2列出了3种方法计算的dp和dl,分析表中数据可以发现:①NECTLS计算的dp和dl相等,WLS和WTLS的不相等;②WLS、WTLS和NECTLS计算的正交距离平方和分别为:0.824 0、0.835 3、0.667 4,NECTLS的值最小。它们与试验1的统计结果一致。

表 2 测量点到拟合点、拟合直线的距离Tab. 2 Distance from measured point to the fitted point and the fitted line respectively

idp
dl
WLSWTLSNECTLS
WLSWTLSNECTLS
10.200 10.420 00.027 8
0.170 80.378 60.027 8
20.150 40.352 40.044 1
0.128 30.317 80.044 1
30.600 60.214 60.372 7
0.512 60.193 70.372 7
40.088 00.368 60.200 1
0.075 10.333 00.200 1
50.584 40.385 70.403 2
0.498 70.355 30.403 2








    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多