《测绘学报》
构建与学术的桥梁 拉近与权威的距离
大角度立体像对相对定向的混合共轭梯度算法
李佳田1,2, 王聪聪1,2
1. 昆明理工大学国土资源工程学院, 云南 昆明 650093; 2. 昆明理工大学云南省高校高原山区空间信息测绘技术应用工程研究中心, 云南 昆明 650093
收稿日期：2017-11-27；修回日期：2018-10-20
基金项目：国家自然科学基金（41561082；41161061）
第一作者简介：李佳田(1975-), 男, 博士, 教授, 研究方向为数值最优化方法与机器场景理解。E-mail:ljtwcx@163.com
通信作者：王聪聪, E-mail：1083719493@qq.com
摘要：无初值依赖的快速收敛是大角度相对定向解算的关键所在。为此，本文提出一种混合共轭梯度算法，具体过程是：①采用随机爬山算法对给定的相对定向元素初值进行随机扰动，产生保证优化方向的初值；②局部优化中以超线性收敛的共轭梯度法取代相对定向中的最速下降法，以提高其收敛速度；③全局收敛条件为计算误差小于规定的限差。对比试验表明，混合共轭梯度算法无初值依赖性，具有较高的解算精度和较少的迭代次数。
关键词：相对定向 大角度 全局收敛 随机爬山算法 共轭梯度法
A hybrid conjugate gradient algorithm for solving relative orientation of big rotation angle stereo pair
LI Jiatian1,2, WANG Congcong1,2, JIA Chenglin1,2, NIU Yiru1,2, WANG Yu1,2, ZHANG Wenjing1,2, WU Huajing1,2, LI Jian1,2
1. Faculty of Land Resource Engineering, Kunming University of Science and Technology, Kunming 650093, China;2. Surveying and Mapping Geo-Informatics Technology Research Center on Plateau Mountains of Yunnan Higher Education of Kunming University of Science and Technology, Kunming 650093, China
Foundation support: The National Natural Science Foundation of China (Nos. 41561082; 41161061)
First author: LI Jiatian(1975—), male, PhD, professor, majors in numerical optimization method and machine scene understanding.E-mail:ljtwcx@163.com
Corresponding author: WANG Congcong, E-mail: 1083719493@qq.com
Abstract: The fast convergence without initial value dependence is the key of large angle relative directional solution. Therefore, a hybrid conjugate gradient algorithm is proposed in this paper. The concrete process is:① stochastic hill climbing(SHC) algorithm is used to make random disturbance to the given initial value of the relative directional element, and the new value to guarantee the optimization direction is generated; ② In local optimization, super-linear convergent conjugate gradient method is used to replace the steepest descent method in relative orientation to improve its convergence rate; ③ The global convergence condition is that the calculation error is less than the prescribed limit error. The comparison experiment shows that the method proposed in this paper is independent of initial value, has higher accuracy and fewer iterations.
Key words: relative orientation big rotation angle global convergence stochastic hill climbing algorithm conjugate gradient algorithm
立体像对相对定向是恢复摄影时刻相邻两像片摄影光束的相互关系，从而使同名光线对对相交[1]。相对定向算法研究具有典型的实际意义。
在航空摄影测量中，立体像对相对定向时像片间的旋偏角一般小于6°，可对相对定向元素中的基线分量近似，以简化相对定向模型，在相对定向元素初值为0的前提下，通常采用最速下降法迭代求解[1-2]。随着倾斜摄影测量和低空摄影测量的发展，立体像对相对定向时像片间的旋偏角可能达到甚至超过10°，最速下降法常因相对定向元素初值无法获取而失败[3]，一般采用相对定向线性变换法求解[4]。上述两种算法的不足之处是：最速下降法具有线性收敛速度、求解精度高，但过分依赖相对定向元素初值[5-6]，使得大角度相对定向时收敛速度减慢、收敛至局部极值甚至不收敛[7]；相对定向线性变换是一种直接解法，虽然无须预知相对定向元素的初值，但求解精度不高，一般作为光束法平差的初值[8]。因此，两种算法均较难兼顾收敛性与求解精度两个方面的要求。
关于初值依赖问题，国内外学者已取得丰硕的成果，主要分为代数法和矩阵分解法两类[9]。代数法是根据共面方程的约束关系，将相对定向元素引入约束方程，构造含有未知数的约束方程组[7, 10-16]。其中，文献[7]将计算机视觉中的8点算法引入相对定向中，无须初值，仅需8个同名点即可解算相对定向元素；较之8点算法，5点算法因有较高精度而被关注，因其采用多项式求解，存在多个解及误解排除的问题[13-16]。矩阵分解法是通过分解基础矩阵或本质矩阵解算相对定向元素[17-19]。文献[18]在顾及基础矩阵元素非独立条件下实现无须初值的相对定向元素求解。文献[19]从基础矩阵中直接导出立体像对像片间的平移和旋转量，较难得到精确的相对定向元素。近年来，针对最速下降法具有较高计算精度，依赖于初值的特点，有学者提出混合算法以克服初值依赖并提高相对定向的解算精度。文献[20]通过分解本质矩阵得到相对定向元素初值，然后采用最速下降法迭代求解。文献[21]结合遗传算法的全局收敛性与最速下降法的局部收敛性，利用两者的混合算法求解大角度立体像对相对定向。
初值依赖性、数值收敛性及求解精度是解算相对定向需考虑的3个问题。文献[22]对常规模拟退火法产生新的待估参数及接受概率进行改进，提出了随机爬山(stochastic hill climbing，SHC)算法。它具有全局收敛性，且较之遗传算法与模拟退火法，其具有收敛速度快的优点。因此，为解决大角度相对定向时，最速下降法对相对定向元素初值依赖的问题，并加快算法的收敛速度，本文提出一种基于SHC算法与共轭梯度的混合共轭梯度算法，采用超线性收敛的共轭梯度法替代最速下降法，加快局部收敛速度。当相对定向元素初值未知时，为避免共轭梯度法陷入局部极值并加快全局收敛速度，利用SHC算法为其提供迭代初值。
1 混合共轭梯度法解算大角度相对定向1.1 相对定向模型
连续像对相对定向是以左像片为基准，确定右像片相对于左像片的相对定向元素(BX, BY, BZ,φ,ω,κ)，如图 1所示。
图 1 相对定向示意图Fig. 1 Schematic diagram of relative orientation
图 1中以左像片的像空间坐标系为像对的像空间辅助坐标系，记为S1-X1Y1Z1，以右摄影中心S2为原点建立另一个像空间辅助坐标系，记为S2-X2Y2Z2。两坐标系互相平行，且左右像点a1,a2在各自像空间辅助坐标系中的坐标分别为(X1,Y1,Z1)和(X2,Y2,Z2)，像点坐标分别为(x1,y1)和(x2,y2，S2)在(S1-X1Y1Z1)中的坐标为(BX, BY, BZ)，则每对同名像点的共面条件方程为[1]
(1)
式中
(2)
式中，f为摄像机焦距；R为相对定向φ、ω、κ角元素方向余弦构成的矩阵。
相对定向中，由于立体像对像片间的旋偏角较小，一般采用简化模型求解，即基线分量BX已知，BY和BZ用小角度μ、ν近似表示
(3)
对式(3)中5个相对定向元素BY、BZ、φ、ω、κ分别求导，采用最速下降法迭代求解。在大角度相对定向情况之下，存在如下问题：
(1) 立体像对像片间的旋偏角较大，式(3)的简化模型不再适用。
(2) 最速下降法具有线性收敛速度，因此，收敛速率方面仍有改进的空间。
(3) 最速下降法是局部收敛算法，在相对定向元素初值难以获取的情况下，计算会陷入局部极值，甚至不收敛。
1.2 混合共轭梯度法
SHC算法是一种改进后的模拟退火法，同样具有全局收敛的性质，它通过对模拟退火算法中新的待估参数产生过程和概率重新修改，获得比模拟退火法更快的收敛速度。设待估参数当前值和新值分别为Xk和Xk 1，f为目标函数，则其只接受f(Xk 1) <f(Xk)的参数，并由下式产生新的待估参数[22]
式中，random为[0, 1]均匀分布的随机数；当x≥0时，sign=1，否则sign=-1；dx为给定的扰动步长，控制待估参数的修改幅度。
为进一步加快算法收敛速度，在局部搜索过程中，考虑采用共轭梯度法进行计算，即混合共轭梯度算法，计算过程如图 2所示。
图 2 混合共轭梯度法Fig. 2 Hybrid conjugate gradient
(1) 按式(4)扰动当前待估参数值来产生新值。
(2) 判断产生的新值是否符合要求。
(3) 若符合要求，则接受新值，并以新值为初值进行共轭梯度法迭代，否则，返回步骤(1)继续扰动产生新值。
(4) 判断是否满足全局收敛条件，满足则结束计算，否则，则返回(1)。
结合混合共轭梯度算法的特点及1.1节大角度相对定向存在的3个问题，本文具体改进是：
(1) 采用具有全局收敛性的SHC算法进行全局搜索，为局部搜索提供迭代初值。
(2) 不采用小角度μ、ν近似，BX已知的简化模型，而直接对BX、BY、BZ、φ、ω、κ这6个相对定向元素求导。
(3) 采用共轭梯度法代替最速下降法，加快局部搜索的收敛速度。
1.3 共轭梯度法
由上文可知，求解式(1)中6个相对定向元素BX、BY、BZ、φ、ω、κ，至少需要6对同名像点。为提高解算精度，常有多余观测。设有n对同名像点(n>6)，则有由n个式(1)构成的方程组
式中，p=[BX, BY, BZ,φ,ω,κ]T为相对定向元素。式(5)是典型的非线性优化过程[23]。共轭梯度法的求解思路是：给定一个相对定向元素初值p0，假定初值p0在真值p*邻域内，则相对定向元素可通过式(6)迭代计算
式中
式中，αk为1维搜索步长，通过线性搜索获得；Pk为搜索方向；βk为参数。不同的βk对应不同算法，文中采用Polak-Ribiere-Polyak算法[24]
(9)
式中，gk=∇H(pk)是函数H(p)在当前相对定向元素估计值pk处的梯度矩阵，即
(10)
进一步由式(1)可得
(11)
因此，在给定1维搜索步长αk的情况下，若相对定向元素初值充分接近真值，则相对定向问题可利用共轭梯度法多次迭代式(6)解算。
1.4 SHC算法[25]
在大角度相对定向中，若相对定向元素初值充分接近真值，则大角度相对定向问题可由超线性收敛的共轭梯度法解算。但其相对定向元素初值一般较难获取，若将相对定向将初值设为0，则会使共轭梯度法收敛于局部极小值，甚至不收敛。为此，在大角度相对定向元素初值未知的情况下，采用SHC算法计算相对定向元素初值。具体过程如下：
(1) 已知左右像片n对同名像点坐标和摄影机焦距f，给定相对定向元素初值pk。
(2) 由式(2)计算左右像对同名像点在各自像空间辅助坐标系的坐标(X1,Y1,Z1)和(X2,Y2,Z2)。
(3) 由式(1)、式(5)计算得到H(pk)。
(4) 按式(4)修改当前相对定向元素pk以产生新的相对定向元素pnew，计算得到ΔH=H(pnew)-H(pk)。
(5) 若ΔH>0，则拒绝pnew，返回步骤(4)；否则接受pnew，即pk=pnew，并以更新后的pk作为初值按共轭梯度法进行相对定向元素计算，直至收敛于局部极值，设收敛的局部极值为pk*，对应的目标函数值为H(pk*)。
(6) 用收敛的局部极值pk*更新当前相对定向元素，即pk=pk*。若未达到全局收敛条件(若真值p*已知，则全局收敛条件为H(pk)=H(p*)；若真值未知，则为算法是否达到指定的迭代次数或规定的限差)，则返回步骤(3)继续执行；否则计算结束，此时的pk即为所求的相对定向元素。
2 试验与分析
为验证本文算法的有效性，分别进行模拟试验和实测试验。
2.1 模拟试验
已知摄像机焦距f及立体像对左像片像点坐标，以左像片为基准，采用交向摄影模拟产生3组右像片分别与左像片构成虚拟立体像对，相对定向元素模拟值见表 1，其中仅考虑了大角度相对定向。然后分别计算不同相对定向元素下右像片的像点坐标，并在其中加入均值为0，方差为0.05的高斯噪声。为模拟大角度立体像对，假设成像不受图幅限制。以左右像片的像点坐标、摄像机焦距f为已知条件，采用本文算法、最速下降法及文献[21]中基于最速下降法的混合遗传算法(简称为算法3)分别对相对定向元素求解，并将相对定向元素计算值与模拟值进行比较。试验中，相对定向线元素和角元素的模型扰动步长dx分别为0.05 mm和0.02 rad。
表 1 相对定向元素模拟值Tab. 1 Simulative value of relative orientation element
| 基线向量/mm | 旋转角度/rad | ||||
| BX | BY | BZ | φ | ω | κ |
| 0.989 7 | -0.086 5 | 0.195 2 | 0.808 5 | -0.483 3 | 0.675 1 |
| -0.956 2 | -0.097 5 | -0.186 2 | 0.390 6 | -0.204 8 | -0.903 6 |
| 0.975 3 | 0.070 3 | -0.245 9 | -0.949 1 | -0.501 3 | -0.661 8 |
由表 1可知，由于相对定向元素模拟值倾角较大，为保证最速下降法正确收敛，试验中采用相对定向线性变换法[4]为其提供初值。而本文算法和算法3中相对定向元素初值均设为0，即BX=BY=BZ=φ=ω=κ=0。3种算法计算结果见表 2，本文算法和算法3的迭代次数见表 3。
表 2 相对定向元素计算结果Tab. 2 Calculating results of relative orientation element
| 像对编号 | 算法 | 基线向量/mm | 旋转角度/rad | ||||
| BX | BY | BZ | φ | ω | κ | ||
| 左-右1 | 本文算法 | 0.989 6 | -0.086 4 | 0.195 1 | 0.808 4 | -0.483 3 | 0.675 2 |
| 最速下降法 | 0.989 6 | -0.086 3 | 0.195 0 | 0.808 4 | -0.483 2 | 0.675 2 | |
| 算法3 | 0.989 5 | -0.086 1 | 0.194 8 | 0.808 3 | -0.483 1 | 0.675 4 | |
| 左-右2 | 本文算法 | -0.956 3 | -0.097 4 | -0.186 2 | 0.390 6 | -0.204 8 | -0.903 6 |
| 最速下降法 | -0.956 3 | -0.097 5 | -0.186 1 | 0.390 6 | -0.204 8 | -0.903 6 | |
| 算法3 | -0.956 0 | -0.097 6 | -0.186 0 | 0.390 7 | -0.204 6 | -0.903 7 | |
| 左-右3 | 本文算法 | 0.975 2 | 0.070 3 | -0.245 8 | -0.949 1 | -0.501 2 | -0.661 7 |
| 最速下降法 | 0.975 2 | 0.070 3 | -0.245 8 | -0.949 2 | -0.501 3 | -0.661 8 | |
| 算法3 | 0.975 5 | 0.070 4 | -0.246 1 | -0.948 9 | -0.501 2 | -0.661 6 | |
表 3 迭代次数Tab. 3 Number of iteration
| 像对编号 | 算法 | 迭代次数 |
| 左-右1 | 本文算法 | 45 |
| 算法3 | 56 | |
| 左-右2 | 本文算法 | 40 |
| 算法3 | 62 | |
| 左-右3 | 本文算法 | 43 |
| 算法3 | 60 |
由表 2、表 3可知，本文算法的计算精度高于算法3，与最速下降法相当，但最速下降法依赖于相对定向元素初值，为保证其正确收敛，需预知相对定向元素初值；而本文算法无须提供初值，且较之算法3，收敛速度较快。原因在于：混合共轭梯度算法中采用了共轭梯度法，具有超线性收敛速度，虽然可能陷入局部极值，但由于采用SHC算法通过扰动产生新的随机相对定向元素，并按目标函数值对产生的相对定向元素进行取舍，保证收敛方向和速度。因此，在多次迭代中，能够跳出共轭梯度法陷入的局部极值。试验中，当相对定向元素初值距离模拟值较远时，混合共轭梯度算法在不同极值处进行多次局部优化，并在每次局部搜索过程中，因SHC算法对所产生的相对定向元素进行判断，越来越接近模拟值，经多次迭代后，最终收敛。
2.2 实测试验
如图 3所示，为云南省某地居民楼房两条基线下交向摄影获得的3幅像片，构成2对连续立体像对。航拍相机为Canon 5D Mark Ⅱ，镜头焦距为35 mm，相对航高445 m，图像大小为449×300像素，像素大小为6.4 μm。为减少同名点定位带来的外部误差，采用人工量测的8个同名像点来进行相对定向元素的解算。采用最速下降法、算法3和本文算法分别对立体像对进行相对定向元素解算。由于立体像对像片间的旋偏角较大，试验中采用相对定向线性变换法[4]为最速下降法提供初值，而本文算法和算法3中，相对定向元素初值设为0，即BX=BY=BZ=φ=ω=κ=0，3种算法的计算结果见表 4，本文算法和算法3的迭代次数见表 5。试验中相对定向线元素和角元素的模型扰动步长dx分别为0.05 mm和0.02 rad。
图 3 像片序列Fig. 3 Sequence of images
表 4 相对定向元素计算结果Tab. 4 Calculating results of relative orientation element
| 立体像对 | 算法 | BX/mm | BY/mm | BZ/mm | φ/rad | ω/rad | κ/rad | 精度/μm |
| 像片2/像片1 | 最速下降法 | -1.367 5 | -0.016 2 | 0.052 4 | -1.064 7 | 0.069 8 | -1.570 8 | 2.1 |
| 算法3 | -1.367 4 | -0.016 2 | 0.052 1 | -1.064 5 | 0.069 8 | -1.570 6 | 2.4 | |
| 本文算法 | -1.367 5 | -0.016 3 | 0.052 5 | -1.064 8 | 0.069 9 | -1.570 7 | 2.2 | |
| 像片3/像片1 | 最速下降法 | -1.391 1 | 0.024 5 | 0.041 1 | -1.029 7 | 0.025 4 | 0.043 6 | 1.8 |
| 算法3 | -1.391 0 | 0.024 5 | 0.041 3 | -1.029 6 | 0.025 4 | 0.043 8 | 2.1 | |
| 本文算法 | -1.391 1 | 0.024 4 | 0.041 0 | -1.029 7 | 0.025 5 | 0.043 7 | 1.8 |
表 5 迭代次数Tab. 5 Number of iteration
| 立体像对 | 算法 | 迭代次数 |
| 像片2/像片1 | 本文算法 | 58 |
| 算法3 | 67 | |
| 像片3/像片1 | 本文算法 | 54 |
| 算法3 | 61 |
从表 4、表 5可知，本文算法计算精度与最速下降法相当，相对定向达到约±0.5像素的精度，而最速下降法因依赖于相对定向元素初值，需要相对定向线性变换法提供初值；较之算法3，本文算法具有计算精度高且收敛快的优点。
3 结论
提出了一种全局收敛的混合共轭梯度算法。试验表明，在无须预知相对定向元素初值的前提下，算法能够正确解算大角度立体像对相对定向，并具有解算精度高和收敛速度快的特点。此外，合理选取SHC算法中的扰动步长将会进一步加快收敛速度。针对大角度相对定向解算，扰动步长的定量分析与选取有待进一步研究。
【引文格式】李佳田, 王聪聪, 贾成林, 等. 大角度立体像对相对定向的混合共轭梯度算法. 测绘学报，2019，48(3)：322-329. DOI: 10.11947/j.AGCS.2019.20170672
精
彩
回
顾
《测绘学报》2019年第3期网刊发布
关于人类历史上首张黑洞照片，你想知道的都在这里
开放报名！首届DataEarth《测绘学报》杯开发者大赛万元奖金等你来拿
今年的硕士、博士确定能毕业么？教育部今年拟抽检6000篇学位论文
2019年QS世界大学学科排名出炉
杨必胜、张小红、赵齐乐等测绘信息领域专家入选第四批国家“万人计划”入选人员
关于召开“测绘前沿科技大讲堂与科技期刊论文写作理论与方法高级研修班”的（一号）通知
权威 | 专业 | 学术 | 前沿
微信投稿邮箱 | song_qi_fan@163.com
微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。
欢加入《测绘学报》作者QQ群： 297834524
进群请备注：姓名 单位 稿件编号
,
