大角度立体像对相对定向的混合共轭梯度算法

2019-04-11李佳田王聪聪贾成林牛一如张文靖吴华静

测绘学报 2019年3期

李佳田，王聪聪，贾成林，牛一如，王瑜，张文靖，吴华静，李键

1. 昆明理工大学国土资源工程学院，云南昆明 650093； 2. 昆明理工大学云南省高校高原山区空间信息测绘技术应用工程研究中心，云南昆明 650093

立体像对相对定向是恢复摄影时刻相邻两像片摄影光束的相互关系，从而使同名光线对对相交[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)

Xk+1=Xk+sign·random·dx

(4)

式中，random为[0，1]均匀分布的随机数；当x≥0时，sign=1，否则sign=-1；dx为给定的扰动步长，控制待估参数的修改幅度。

为进一步加快算法收敛速度，在局部搜索过程中，考虑采用共轭梯度法进行计算，即混合共轭梯度算法，计算过程如图2所示。

(1) 按式(4)扰动当前待估参数值来产生新值。

(2) 判断产生的新值是否符合要求。

(3) 若符合要求，则接受新值，并以新值为初值进行共轭梯度法迭代，否则，返回步骤(1)继续扰动产生新值。

(4) 判断是否满足全局收敛条件，满足则结束计算，否则，则返回(1)。

结合混合共轭梯度算法的特点及1.1节大角度相对定向存在的3个问题，本文具体改进是：

(1) 采用具有全局收敛性的SHC算法进行全局搜索，为局部搜索提供迭代初值。

(2) 不采用小角度μ、ν近似，BX已知的简化模型，而直接对BX、BY、BZ、φ、ω、κ这6个相对定向元素求导。

(3) 采用共轭梯度法代替最速下降法，加快局部搜索的收敛速度。

图2 混合共轭梯度法Fig.2 Hybrid conjugate gradient

1.3 共轭梯度法

由上文可知，求解式(1)中6个相对定向元素BX、BY、BZ、φ、ω、κ，至少需要6对同名像点。为提高解算精度，常有多余观测。设有n对同名像点(n>6)，则有由n个式(1)构成的方程组

H(p)=[h1(p),h2(p),…,hn(p)]T=0

(5)

式中，p=[BX,BY,BZ,φ,ω,κ]T为相对定向元素。式(5)是典型的非线性优化过程[23]。共轭梯度法的求解思路是：给定一个相对定向元素初值p0，假定初值p0在真值p*邻域内，则相对定向元素可通过式(6)迭代计算

pk+1=pk+αkPk

(6)

式中

P0=-g0

(7)

Pk+1=-gk+1+βkPk

(8)

式中，αk为1维搜索步长，通过线性搜索获得；Pk为搜索方向；βk为参数。不同的βk对应不同算法，文中采用Polak-Ribiere-Polyak算法[24]

(9)

(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)。

2 试验与分析

为验证本文算法的有效性，分别进行模拟试验和实测试验。

2.1 模拟试验

已知摄像机焦距f及立体像对左像片像点坐标，以左像片为基准，采用交向摄影模拟产生3组右像片分别与左像片构成虚拟立体像对，相对定向元素模拟值见表1，其中仅考虑了大角度相对定向。然后分别计算不同相对定向元素下右像片的像点坐标，并在其中加入均值为0，方差为0.05的高斯噪声。为模拟大角度立体像对，假设成像不受图幅限制。以左右像片的像点坐标、摄像机焦距f为已知条件，采用本文算法、最速下降法及文献[21]中基于最速下降法的混合遗传算法(简称为算法3)分别对相对定向元素求解，并将相对定向元素计算值与模拟值进行比较。试验中，相对定向线元素和角元素的模型扰动步长dx分别为0.05 mm和0.02 rad。

由表1可知，由于相对定向元素模拟值倾角较大，为保证最速下降法正确收敛，试验中采用相对定向线性变换法[4]为其提供初值。而本文算法和算法3中相对定向元素初值均设为0，即BX=BY=BZ=φ=ω=κ=0。3种算法计算结果见表2，本文算法和算法3的迭代次数见表3。

表1 相对定向元素模拟值

表2 相对定向元素计算结果

表3 迭代次数

由表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

立体像对算法BX/mmBY/mmBZ/mmφ/radω/radκ/rad精度/μm像片2/像片1最速下降法-1.3675-0.01620.0524-1.06470.0698-1.57082.1算法3-1.3674-0.01620.0521-1.06450.0698-1.57062.4本文算法-1.3675-0.01630.0525-1.06480.0699-1.57072.2像片3/像片1最速下降法-1.39110.02450.0411-1.02970.02540.04361.8算法3-1.39100.02450.0413-1.02960.02540.04382.1本文算法-1.39110.02440.0410-1.02970.02550.04371.8