CN104019799A - 一种利用局部参数优化计算基础矩阵的相对定向方法 - Google Patents

Abstract

本发明提供了一种利用局部参数优化计算基础矩阵的相对定向方法，包括如下步骤：设置多个空间特征点；拍摄并获取所述多个空间特征点在多个位置点处的像面坐标；用线性方法求基础矩阵初值；用本质矩阵和相机内参数表示基础矩阵，根据本质矩阵的具有两个相同的特征值且秩为2的特性，将表征两相机相对位置关系的外极线几何方程参数化表示，将其表示为像点坐标与一系列旋转矩阵相乘的形式，通过高斯牛顿方法迭代更新包含外参数的矩阵，精确解算外参数。本发明的相对定向方法简化了外极线几何误差数学模型，而且在求解误差函数的雅可比矩阵时让各旋转角的初值都取0，使得建模和计算都变得简单高效率。

Description

一种利用局部参数优化计算基础矩阵的相对定向方法

技术领域

本发明涉及摄影测量技术领域，特别涉及一种基础矩阵精确解算的非线性迭代方法。

背景技术

近年来，随着以机电重工、汽车、船舶、通信电线，航空航天等为代表的现代大型装备制造业不断发展，大范围现场高精度测量技术在装备、对接和质量控制等任务环节中的需求越来越多，因此便携、快速、高精度、大尺寸测量成为数字化摄影精密测量技术的一个重要发展方向。摄影测量系统的外方位参数解算是直接影响测量系统的关键因素。然而传统定向方法的定向精度单纯依赖于定向靶的加工和校准精度，且在测量过程中定向靶位置固定，测量范围受限。当测量范围较大时，就需要在测量构件不同位置多次摆放靶标才能完成对大尺寸或超大尺寸构件的整体测量，测量效率低且测量精度不高。为了克服传统定向方式的缺陷，扩展测量范围和测量灵活性，从计算机视觉中引入了基于外极线几何约束的两图像相对定向方法。计算机视觉中最经典的双图像相对定向方法是八点算法，利用外极线几何约束关系最少用两幅图像中公共的8个对应像点就可以计算出包含两幅图像间外方位关系的基础矩阵。基础矩阵是分析同一景物两幅图像外极线几何约束关系的有力工具，它独立于场景结构，包含了相机的内参数和外方位参数全部信息，被广泛应用于图像匹配、运动估计、三维重建、相机标定及虚拟现实等计算机视觉实际应用问题中，其较高精度的求解是后续工作的基础，因此，精确、鲁棒的估计基础矩阵是非常重要的。

由于实际图像观测值存在误差，基础矩阵的求解问题转化成最小化误差优化问题，基础矩阵估计的关键问题是：选择最小化优化误差函数，基础矩阵参数化表示方法及估计方法。

目前计算基础矩阵的常用方法主要分为三类：线性方法，非线性方法和鲁棒方法。线性方法运算速度快，但忽略基础矩阵各参数间的约束关系且最小化模型缺乏物理意义，对噪声非常敏感，所以精度和鲁棒性较差。而非线性方法和鲁棒方法可以得到更精确的结果，但是相当耗时，不满足在线测量要求。

线性算法中最常用的是Hartley提出的归一化8点算法，经过对图像坐标归一化和基础矩阵的秩为2约束的加入，能够有效消除由于噪声干扰产生的图像观测值误差，与其他线性迭代算法相比，计算简单，运算速度较快。然而，由于基础矩阵是外极线几何的数学表达形式，包含了相机的内参数和外参数全部信息，理论上所有极线交于一点，基础矩阵秩为2，但由于计算误差及像点坐标定位误差，实际中的极点并不交于一点，而是近似交于一点，所以实际中基础矩阵秩不为2，这就需要通过归一化八点算法将求出的基础矩阵强加入秩为2的约束，然后采用非线性迭代方法对求解出的初始基础矩阵进行迭代更新，使其无限逼近最优解。

常用的非线性方法和鲁棒方法主要有梯度法、M估计法、随机一致性采样法(RANSAC)，但是这些算法由于数学模型比较复杂，运算速度慢，不能很好的满足在线测量要求。因此，需要一种高精度并且运算简便快捷的解算基础矩阵的方法。

发明内容

本发明的目的在于提供一种高效的解算基础矩阵的相对定向方法，即一种利用局部参数优化计算基础矩阵的相对定向方法，所述方法包括如下步骤：

(1).在测量场中布置多个空间特征点；

(2).在大于或等于两个的不同站位由相机对所述多个空间特征点进行拍摄测量，获取在每个站位处拍摄的每一幅图像中的n个空间特征点各自的像面坐标，其中n≥8；

(3).获得相机内参数，从上述图像中选择具有至少8个公共空间特征点的两幅图像，则两图像外极线几何约束用本质矩阵和像点坐标描述为如下式所示：

m^'TEm＝0      (1)

其中m(m,n,c)^T和m'(m',n',c')^T代表同一空间特征点在两幅图像中对应像点在相机坐标系下的坐标，E是表征两图像中相机外方位参数的本质矩阵；

(4).用局部参数优化的方法对拍摄获得的每一组对应像点坐标构成的外极线几何约束式(1)中的本质矩阵进行如下参数优化过程：

1)首先将本质矩阵E进行奇异值分解，即其中对角矩阵U、V都是酉矩阵；再将本质矩阵E参数化表示，即，将U、V视为旋转矩阵，利用正交坐标系的三坐标轴上的旋转角，分别将U、V参数化表示成欧拉角分解的形式，即其中，R_u,R_x是绕x轴的旋转矩阵，R_v,R_y是绕y轴的旋转矩阵，R_w,R_z是绕z轴的旋转矩阵，R_z与具有互换性，将R_w和R_z两个旋转矩阵进行合并，将误差函数的自由参数降至五个，这时本质矩阵用旋转矩阵相乘表示为：

$E=R_u{R}_v{R}_w\hat{I}{R}_y^T{R}_x^T---\left(2\right)$

其中

$R_u=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos u& -\sin u\\ 0& \sin u& \cos u\end{array}\right],R_v=\left[\begin{array}{ccc}\cos v& 0& \sin v\\ 0& 1& 0\\ -\sin v& 0& \cos v\end{array}\right],R_w=\left[\begin{array}{ccc}\cos w& -\sin w& 0\\ \sin w& \cos w& 0\\ 0& 0& 1\end{array}\right]$

2)用高斯牛顿迭代法，迭代更新酉矩阵U、V，定义为误差函数，m_i＝(m_i,n_i,c_i)^T与m'_i＝(m′_i,n′_i,c′_i)^T是第i个空间特征点在两幅图像中的对应像点在相机坐标系下的坐标，定义角度参数变量θ＝(θ_u,θ_v,θ_w,θ_x,θ_y)^T，计算误差函数ε_i在角度初值θ＝(0,0,0,0,0)^T时的雅可比矩阵：

$J_i{|}_0=\frac{{\∂\mathrm{\ϵ}}_i}{\∂\mathrm{\θ}}{|}_{\mathrm{\θ}=0}=(n_i{c}_i^{\′},-m_i{c}_i^{\′},m_i{n}_i^{\′}-n_i{m}_i^{\′},c_i{n}_i^{\′},-c_i{m}_i^{\′})---\left(3\right)$

a)利用n个空间点的像面坐标和相机内参数，用归一化8点算法解算本质矩阵初值，并对该初值进行奇异值分解利用分解出的酉矩阵U和V进行像点坐标转换：即，将归一化后的对应像点坐标m_i(m_i,n_i,1)^T，m′_i(m′_i,n′_i,1)^T(i＝1,...,n,n≥8)左乘U^T、V^T表示为X'＝U^Tm′_i，X＝V^Tm_i并代入方程(1)，得到误差函数

b)将转换后的像点坐标X＝(X_i,Y_i,Z_i)^T与X'＝(X′_i,Y′_i,Z′_i)^T代入上面的误差函数，则方程简化为

c)将转化后的像点坐标X、X’代入(3)式得到如下雅可比矩阵J，由Jδ_θ＝-ε_i推得参数变量δ_θ＝-(J'J)^-1J'ε_i；

$J=\frac{{\∂\mathrm{\ϵ}}_i}{\∂\mathrm{\θ}}{|}_{\mathrm{\θ}=0}=(Y_i{Z}_i^{\′},-X_i{Z}_i^{\′},X_i{Y}_i^{\′},-Y_i{X}_i^{\′},Z_i{Y}_i^{\′},-Z_i{X}_i^{\′})---\left(4\right)$

d)用求解出来的δ_θ＝(δ_u,δ_v,δ_w,δ_x,δ_y)^T迭代更新原始的矩阵U、V得到新的U_new和V_new如下：

U_new＝UR_u(δ_u)R_v(δ_v)R_w(δ_w)           (5)

V_new＝VR_x(δ_x)R_y(δ_y)

重复上述步骤，当||ε||满足预设的精度要求或收敛条件时，停止迭代，输出最终获得的本质矩阵

3)利用Hartley提出的Cheirality约束条件从本质矩阵E中解算出正确的两相机外方位参数中的平移向量T和旋转矩阵R，完成基于两图像的相对定向方法的流程，然后对其他具有至少8个公共特征点的两幅图像重复上述步骤1)和2)，完成所有图像间的相对定向。

优选的，所述空间特征点为具有特定的空间约束关系的编码点，直径为3mm。

优选的，所述空间特征点均匀分布在被测场中，且所有特征点不能共面，保证两相邻站位间的公共空间特征点数不能少于8个。

优选的，在对空间特征点进行拍摄时采用单相机系统或双相机系统，在所述单相机系统中，物体不动，相机运动，以从不同站位拍摄，在双相机系统中，物体运动而两相机不动。

优选的，所述首先对本质矩阵进行奇异值分解，然后对本质矩阵的两个相同特征值归一化处理，将本质矩阵分解出的特征值矩阵表示成对角矩阵最后再将本质矩阵参数化表示，简化误差函数的雅可比矩阵。

优选的，所述迭代停止条件为图像中的所有像点到对应外极线距离的均方根误差不再收敛或δ_θ≤10^-7。

优选的，所述误差函数的雅可比矩阵在每次迭代中的角度初值都取0，完成误差函数模型在0角度的局部迭代优化。

优选的，在所述双相机系统测量动态物体情况下，通过对贴有编码点的基准尺在测量场不同位置进行拍摄测量，完成两相机间的初始相对定向。

优选的，在所述单相机系统的情况下，完成所有图像间的相对定向后，选取包含测量场中编码点最多的图像作为世界坐标系，建立所有图像的相对方位转换链，将所有图像统一到一个世界坐标系下。

优选的，在所述需要用经过标定的贴有编码点的基准尺作为尺度基准，来唯一确定外参数中平移向量T的比例因子。

应当理解，前述大体的描述和后续详尽的描述均为示例性说明和解释，并不应当用作对本发明所要求保护内容的限制。

附图说明

参考随附的附图，本发明更多的目的、功能和优点将通过本发明实施方式的如下描述得以阐明，其中：

图1示出了双目立体视觉系统示意图。

图2示出了根据本发明具体实施方式的相对定向方法的流程图。

图3示出了各成像站位相对方位转换图。

图4a-4e分别示出了实验中使用的相机、基准尺、编码特征点、回光反射标志点及定向靶。

图5示出了特征像点到对应外极线距离误差。

具体实施方式

通过参考示范性实施例，本发明的目的和功能以及用于实现这些目的和功能的方法将得以阐明。然而，本发明并不受限于以下所公开的示范性实施例；可以通过不同形式来对其加以实现。说明书的实质仅仅是帮助相关领域技术人员综合理解本发明的具体细节。

在下文中，用参考附图描述本发明的实施例。在附图中，相同的附图标记代表相同或类似的部件，或者相同或类似的步骤。

在理想成像条件下，每个像点都是由空间特征点反射光经透镜中心到达像面，反过来，同一个特征点的所有成像，从各个像面沿光线反方向投射到空间，将汇聚到同一个空间特征点。

如图1所示，由于至少需要用两个像平面上的点的汇聚才可唯一确定一个空间点，因此根据本发明的方法用于由多幅图像获取物体三维几何信息。

在图1所示的双目立体视觉系统中，两个相机坐标系分别为O_c1X_c1Y_c1Z_c1和O_c2X_c2Y_c2Z_c2。空间点P通过光线O_c1P成像于左像面上一点p₁，通过光线O_c2P成像于右像面上一点p₂，p₁和p₂是一对同名点。平面PO_c1O_c2与左右两个像平面分别交于直线l₁和l₂。由于p₁的同名点p₂既位于右像平面上，又位于平面PO_c1O_c2上，因此p₂必位于PO_c1O_c2于右像平面的交线l₂上；同理，p₂的同名点p₁必位于交线l₁上。l₂称为右图上对应于p₁点的极线，l₂称为左图上对应于p₂点的极线。随着空间点P位置的变化，像点和对应的极线在图像上的位置和角度也发生变化，但是，由于所有的PO_c1O_c2平面都相交于直线O_c1O_c2，而O_c1O_c2交两个像平面于固定两点e₁和e₂，故左像平面上所有的极线相交于e₁，右像平面上所有极线相交于e₂。e₁是右相机光心O_c1在左像面的像点，叫做左极点；e₂是左相机光心O_c2在右像面的像点，叫做右极点。这就是外极线几何约束条件。上述外极线几何理论是根据本发明的相对定向方法的数学模型基础。

图1中，p′₁和p′₂是由于系统畸变等原因造成的实际成像点，可以看出，由于像点误差的存在，实际成像点偏离理想成像位置，使得同名点的空间光线在空间中无法交于一点，而是构成异面直线，间隔距离D。

下面举例说明外极线约束与基础矩阵及外方位参数的关系。

基础矩阵F是外极线几何约束的数学表达形式，成3×3的矩阵形式，它包含相机所有的内参数和两相机的外方位参数。外极线约束可以用下面方程表示：

x^'TFx＝0     …(1)

式中，x、x'代表空间特征点P对应的两个像面的像点经过归一化后的齐次像面坐标(这里可采用Hartley的归一化方法或者其他方法)。

把齐次像面坐标x＝(x,y,1)，x'＝(x',y',1)代入方程(1)，则方程(1)改写为：

Uf＝0         …(2)

式中f＝(F₁₁,F₁₂,F₁₃,F₂₁,F₂₂,F₂₃,F₃₁,F₃₂,F₃₃)^T,U＝(xx',yx',x',xy',yy',y',x,y,1)^T。

由于实际采集图像中存在噪声(即像点坐标定位有误差)，使得空间对应点P的特征像点p、p'不能恰好位于对应极线l、l'上，实际上的极点和极线是有一定距离误差的，所以实际极线方程(1)并不等于0，需要考虑误差，这里用残差ε＝|x^'TFx|表示误差的大小。

因此，对于基础矩阵F估计问题就转化为用最小二乘法求解残差ε最小的无约束最优化问题，即：

$\underset{F}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(x^{\′T}\mathrm{Fx}\right)}^2=\underset{F}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i^2...\left(3\right)$

由于(2)式中平凡解f＝0没有实际意义，在最小化求解过程中需要加入某些约束避免平凡解的出现(如加入||f||²＝f^Tf＝1约束)，最常用的方法是奇异值分解法，但是也可以采用别的最小二乘法来对基础矩阵进行修正。

奇异值分解法具体过程如下：

1)图像坐标归一化处理

首先用像点坐标提取算法提取采集到的两图像的像面坐标，然后用预先标定的内参数对两幅图像中编码点像面坐标进行畸变校正，最后用归一化变换矩阵T，T’对两幅图像同时进行归一化处理。参照Hartley的像面坐标归一化方法，第一步先对每幅图像中的像点坐标进行平移使点集的形心移至像素坐标系原点，然后对坐标系进行缩放使得像点到原点的平均距离等于即“平均”特征像点的齐次坐标向量为(1,1,1)^T，称为归一化，这时像面坐标归一化为：x＝Tx,x'＝T'x'，其中T，T’就是归一化变换矩阵，由平移和缩放矩阵组成。

2)初始化矩阵U

式(3)也可写成对应||f||²＝f^Tf＝1，由已知n个对应像点齐次像面坐标由式(2)可以得到n个线性方程组，n≥8，

其中 $A=\left[\begin{array}{ccccccccc}{x}_1{x}_1^{\′}& y_1{x}_1^{\′}& x_1^{\′}& x_1{y}_1^{\′}& y_1{y}_1^{\′}& y_1^{\′}& x_1& y_1& 1\\ .& .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .& .\\ x_n{x}_n^{\′}& y_n{x}_n^{\′}& x_n^{\′}& x_n{y}_n^{\′}& y_n{y}_n^{\′}& y_n^{\′}& x_n& y_n& 1\end{array}\right.$

所以A奇异值分解后对应的最小奇异值的奇异矢量即为基础矩阵F的最优解。

3)强迫约束：加入秩为2的约束，采用奇异值分解修正求出的基础矩阵F，设矩阵F奇异值分解为F＝UDV^T，用代替奇异值分解出的对角矩阵D＝diag(σ₁,σ₂,σ₃)，基础矩阵的降秩近似为

4)解除归一化：恢复基础矩阵此时的矩阵F是对应原始像点坐标x、x’的基础矩阵。

由于加入强迫约束会带来计算误差，需要进一步用非线性的方法对基础矩阵进行修正。

当相机内参数已知时，基础矩阵只有5个自由度，可以用本质矩阵E描述，设两相机的内参数矩阵分别为K₁、K₂，其中C₁，C₂为两相机的主距

$K_1=\left[\begin{array}{ccc}{-C}_1& 0& 0\\ 0& {-C}_1& 0\\ 0& 0& 1\end{array}\right],K_2=\left[\begin{array}{ccc}-C_2& 0& 0\\ 0& {-C}_2& 0\\ 0& 0& 1\end{array}\right];$

将代入(1)式的外极线约束方程可以表示为：

$x^{,T}\mathrm{Fx}=x^{,T}{K}_2^{-T}{\mathrm{EK}}_1^{-1}x=0...\left(4\right)$

把归一化的像点坐标都用内参数矩阵左乘，即像点坐标表示成 则方程外极线约束方程(4)可以表示为：

m^’TEm＝0          …(5)

根据本发明的用局部参数优化的方法对式(5)中的本质矩阵进行迭代优化，具体方法步骤如下：

1)将本质矩阵进行奇异值分解其中因为U和V都是酉矩阵，可以用欧拉角把U、V参数化表示成旋转矩阵相乘的形式即其中R_u和R_x是绕x轴的旋转矩阵，R_v和R_y是绕y轴的旋转矩阵，R_w和R_z是绕z轴的旋转矩阵。R_z的形式决定了它与的可交换性，将R_z移到左边，将R_w和R_z两个旋转矩阵进行合并，这时误差函数有五个自由参数，分别是θ＝(θ_u,θ_v,θ_w,θ_x,θ_y)^T，与本质矩阵的自由度相匹配。这时本质矩阵可以表示为：

$E=R_u{R}_v{R}_w\hat{I}{R}_y^T{R}_x^T...\left(6\right)$

其中

$R_u=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos u& -\sin u\\ 0& \sin u& \cos u\end{array}\right],R_v=\left[\begin{array}{ccc}\cos v& 0& \sin v\\ 0& 1& 0\\ -\sin v& 0& \cos v\end{array}\right],R_w=\left[\begin{array}{ccc}\cos w& -\sin w& 0\\ \sin w& \cos w& 0\\ 0& 0& 1\end{array}\right]$

2)用如下高斯牛顿迭代法，迭代更新酉矩阵U、V，定义误差函数为其中m_i＝(m_i,n_i,c_i)^T，m'_i＝(m′_i,n′_i,c′_i)^T为第i个空间特征点在两幅图像中对应像点在相机坐标系下的坐标，参数变量θ＝(θ_u,θ_v,θ_w,θ_x,θ_y)^T，误差函数ε_i在初值θ＝(0,0,0,0,0)^T时的雅可比矩阵：

$J_i{|}_0=\frac{{\∂\mathrm{\ϵ}}_i}{\∂\mathrm{\θ}}{|}_{\mathrm{\θ}=0}=(n_i{c}_i^{\′},-m_i{c}_i^{\′},m_i{n}_i^{\′}-n_i{m}_i^{\′},c_i{n}_i^{\′}.-c_i{m}_i^{\′})...\left(7\right)$

其中高斯牛顿迭代法的具体步骤如下：

a)相机坐标系下归一化后的像点坐标m_i＝(m_i,n_i,1)^T，m′_i＝(m′_i,n′_i,1)^T，左乘U^T,V^T表示为：X'＝U^Tm',X＝V^Tm，其中X＝(X_i,Y_i,Z_i)^T，X'＝(X′_i,Y′_i,Z′_i)^T；

b)奇异值分解本质矩阵初值

c)计算误差函数

d)用X’，X坐标表示误差函数的雅可比矩阵J，由Jδ_θ＝-ε_i推导出δ_θ＝-(J'J)^-1J'ε_i，其中，J＝(Y_iZ′_i,-X_iZ′_i，X_i,X_iY′_i-Y_iX′_i,Z_iY′_i.-Z_iX′_i)；

e)用d)求解出来的参数变量δ_θ＝(δ_u,δ_v,δ_w,δ_x,δ_y)^T将E₀奇异值分解出的酉矩阵U、V更新为U’、V’，用新的U’、V’代替U、V继续进行迭代运算；

U'＝UR_u(δ_u)R_v(δ_v)R_w(δ_w)         …(8)

V'＝VR_x(δ_x)R_y(δ_y)

f)当||ε||满足设定的精度要求，停止迭代，输出本质矩阵

需要注意的是理论上这种迭代方法可以一直迭代到匹配像点无限接近极线，但是实际中由于像点坐标定位精度的限制，过多迭代次数对解算外参数R、T是没有意义的，假设像点坐标的定位精度为0.02pixel，当像点到其对应极线距离计算到小于0.002pixel时，外参数的求解精度就不会再有提高，此时限制因素是像点坐标的定位精度，因此可以根据实际精度要求设置迭代次数。由于本文用的是最速下降法，计算效率非常高，每次迭代过程仅耗时不超过2μs，迭代2-3次就可以达到外参数求解的极限精度，非常接近线性算法的迭代速度。另外由于舍入误差存在，经过几次迭代后||ε||值在收敛域附近可能会增加，这时保留前一次迭代中U、V的值，然后跳出循环。

3)根据以上步骤获得的U、V解算外方位参数中的平移向量T和旋转矩阵R，参照Hartley提出的Cheirality约束条件，利用直接线性变换重建空间特征点的三维坐标，通过空间特征点与两相机的位置关系确定E分解出四组解中的正确解，用迭代更新后的U、V解算出的四组外方位参数如下：

$\left[\begin{array}{cc}R& T\end{array}\right]=\left\{\begin{array}{c}\left[\begin{array}{cc}{\mathrm{UQV}}^T& T\end{array}\right]\\ \left[\begin{array}{cc}{\mathrm{UQV}}^T& -T\end{array}\right]\\ \left[\begin{array}{cc}{\mathrm{UQ}}^T{V}^T& T\end{array}\right]\\ \left[\begin{array}{cc}{\mathrm{UQ}}^T{V}^T& -T\end{array}\right]\end{array}\right.,T=U\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],Q=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]...\left(9\right)$

如上过程即为基于两图像的相对定向方法的流程，如图2所示，以此类推，将所有图像定向完毕后，以所有图像中包含编码点最多的成像站位作为世界坐标系，建立所有成像站位间的相对方位转换链，如图3确定站位1相对于站位2的方位(T_C1C2,R_C1C2)和站位3相对于站位2的方位(T_C3C2,R_C3C2)，将相对定向得到的所有图像的外参数转换到同一个世界坐标系下,完成整个测量场的多站位全局外参数标定。

上述过程可以采用单相机系统或双相机系统。

在单相机系统中，物体不动，相机运动，并且从不同角度拍摄以得到足够的信息还原被测物体。在双相机系统中，物体是运动的，相机不动，两台相机固定好位置后就不会再改变，然后由两个相机同时拍摄图像。在上述过程中，两者的不同之处在于：计算中，单相机系统中，内参数K1＝K2，双相机系统中，K1不等于K2。双相机系统的约束比单相机系统的约束少，精度低，但只需测量一次即可，而且不用建立所有图像相对方位转换链以完成所有图像间的相对定向，只需完成两相机图像的相对定向即可。本领域技术人员可以根据测量物体的运动状态和相机数量或者测量环境等具体应用情况来选择适当的相机测量模式系统。

下面通过具体实验来说明根据本发明的方法的技术效果。

应用上述根据本发明的方法，通过具体实验来验证本发明评价方法的有效性。实验条件如下，使用尼康D2Xs相机(图4a)，其内参数已标定，见表I。像元尺寸为6μm×6μm，分辨率为4288×2848。实验中分别选择2m×1m、4m×3m、6m×4m3个不同测量范围的测量场对本文的定向方法进行验证，标志点为反射级别为2、直径分别为3mm及6mm的GSI编码点(Geodetic Systems,Inc(GSI)，大地测量系统公司)及圆形RRT(回光反射标志点)，成像圆尺寸分别约为5×5pixels及11×11pixels，如图4b，4c。选取图4d基准尺为尺度基准，其标定长度为1050.681mm，基准尺上两特征点之间的距离为892.5572mm，精度为0.01mm。GSI自动定向靶如图4e所示。

表1相机内参数

其中，C代表：主距，x0、y0代表：主点；K1、K2、K3代表：径向畸变系数；P1、P2代表：偏心畸变系数；AP1、AP2代表：非正交畸变系数。

实验步骤如下进行：

步骤1：按测量场的大小均匀布置编码点和RRT标志点，所有编码点不能共面，在场中心放置定向靶，为了唯一确定相机站位的相对方位，需要在测量空间中设置尺度基准。

步骤2；调整相机参数，使成像的背景噪声小于0.5个灰度级别，即背景灰度基本为零。根据测量场的大小选取合适的站位对测量场拍摄照片，每次实验拍摄约60张左右的图像。

步骤3：用高精度摄影测量系统V-STARS软件(测量精度为12万分之一)对尼康D2Xs相机拍照得到的图像进行计算，得到测量场中相机的外参数和编码点的空间坐标精确值并以此为评价标准衡量不同定向方法的定向精度。

步骤4：利用灰度质心法提取尼康D2Xs相机拍照得到的像点坐标，分别采用定向靶后方交会的绝对定向方法，基于归一化8点算法的相对定向方法及本文提到的基于局部参数优化基础矩阵的相对定向方法三种定向方法解算相机的外方位参数和编码点的空间坐标。

步骤5：将V-STARS及三种定向方法得到的外方位参数，空间坐标都转化到同一个相机坐标系下，以V-STARS计算结果为约定真值，计算三种定向方法的解算误差。

步骤6：重复性测量，重复10次，得到10组数据，求平均值。

基于上述测量和计算结果，进行数据分析。用像点到其对应外极线距离的误差来评价归一化8点算法和本文提出的算法解算基础矩阵的精度，该误差定义为：

$S=\sqrt{\frac{1}{N-1}\underset{i}{\overset{N}{\mathrm{\Σ}}}(d^2(m_i^{\′},{\mathrm{Fm}}_i)+d^2(m_i,F^T{m}_i^{\′}))}=\sqrt{\frac{1}{N-1}\underset{i}{\overset{N}{\mathrm{\Σ}}}\left((\frac{1}{l_1^2+l_2^2}+\frac{1}{l_1^{\′2}+l_2^{\′2}}){\left(u_i^{\′T}{\mathrm{Fu}}_i\right)}^2\right)}...\left(10\right)$

设l＝F^Tu′_i＝(l₁,l₂,l₃)^T，l’＝Fu_i＝(l’₁,l’₂,l’₃)^T分别表示图像I和II的外极线描述，比较结果如图5，以V-STARS计算得到的相机外参数和编码点空间坐标为约定真值，求得三种算法解算旋转矩阵R和平移矢量T的绝对误差精度及编码点空间坐标的相对误差精度，比较结果如表2。

表2外方位参数及空间坐标的计算值与v-stars值比较

其中：算法标号ARCT代表用定向靶空间坐标后方交会的定向结果；算法标号N8P代表用归一化8点线性算法解算基础矩阵的定向结果；算法标号RN8P代表用本发明优化的8点算法解算基础矩阵的定向结果。

表2中，Ω，κ，Xt，Yt，Zt为相机的6个外方位参数。其中，Xt，Yt，Zt是平移向量，用来描述投影中心C相对于空间坐标系的位置；Ω，κ是旋转角度，用来描述影像平面在摄影瞬间的空中旋转姿态，即对于像空间坐标系X、Y、Z轴，以Y为轴旋转角，然后绕X轴旋转Ω角，最后绕Z轴旋转κ角。每个测量范围的第一行数据为用定向靶空间坐标后方交会的定向结果，第二行为用归一化8点线性算法解算基础矩阵的定向结果，第三行为用本发明优化的8点算法解算基础矩阵的定向结果。外参数的绝对误差是指与V-STARS计算得到的外参数结果进行比对后计算均方根RMS误差：

$\mathrm{RMS}=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(e_i-e)}^2}$

其中N表示拍摄的照片数量，e_i表示由定向测量模型计算得到的第i张图像的外参数，e表示V-STARS计算得到的外参数)，编码点空间坐标的相对误差是指与V-STARS计算得到的编码点空间坐标结果进行比对后计算得到的相对于测量空间的坐标误差(编码点空间坐标的均方根误差RMS与测量场范围的比值)，每个测量场重复测量10次，将得到的10组数据取平均值，得到表2的数据。

由表2中显示的实验结果可知，在三种测量范围内，本发明提出的基于优化的8点算法的定向精度是三者中最高的，归一化8点算法的定向精度次之，而用Autobar进行空间坐标交会得到的定向精度最低，需要进一步需要用光束平差方法修正。

需要特别指出的是，外方位参数的定向精度是由环境条件、被摄目标大小、成像站位、成像系统性能及定向算法共同造成的，因此同一算法在不同条件下可得到不同的精度。在实际应用中，可用本发明的定向方法确定某场景下的相机外方位参数。

需要特别指出的是本发明的相对定向方法虽然优化效率很高，但不能处理误匹配的情况，需要具有特定约束关系的编码点预先建立匹配信息才可以进行定向，当有误匹配点存在时，需要用鲁棒性算法(如随机抽样一致性RANSAC)来去除误匹配点，然后再用本文的优化算法精确的解算相机外参数。

根据本发明的方法利用外极线约束关系提出了一种高效率优化基础矩阵的非线性迭代方法，以归一化8点算法为基础，求得初始基础矩阵，当内参数已知时，用本质矩阵和像点坐标表示外极线几何约束方程并参数化表示，通过迭代的方式更新参数，进而更新精度更高的基础矩阵，实验数据表明，相比归一化8点算法，该算法具有更高的精度，相比其他的非线性迭代方法，该算法简单，运算速度更快。

本发明先将本质矩阵进行奇异值分解再参数化表示，将其视作一系列刚体旋转运动，本质矩阵有一个明显的特点，它有两个相同的特征值，而且矩阵秩为2，因此本发明将本质矩阵进行参数局部化，简化了数学模型，而且在求雅可比矩阵时让旋转角的初值都取0，使得建模和计算都变得简单高效率。

根据本发明的方法优化效率很高而且很稳定，可以大幅度提高基础矩阵的精度，使求得的基础矩阵无限逼近最优解。根据本发明的方法精度高，速度快，但局限性在于需要应用中像点定位误差保持在预定范围。一般摄影测量中的像点坐标定位误差是在0.02～0.05像素(实验时像素单位为6um)，换算成微米0.12um～0.3um，误差非常小，因此非常适合这种快速迭代方法。

结合这里披露的本发明的说明和实践，本发明的其他实施例对于本领域技术人员都是易于想到和理解的。说明和实施例仅被认为是示例性的，本发明的真正范围和主旨均由权利要求所限定。

Claims (10)

1.一种利用局部参数优化计算基础矩阵的相对定向方法，所述方法包括如下步骤：

(1).在测量场中布置多个空间特征点；

(2).在大于或等于两个的不同站位由相机对所述多个空间特征点进行拍摄测量，获取在每个站位处拍摄的每一幅图像中的n个空间特征点各自的像面坐标，其中n≥8；

(3).获得相机内参数，从上述图像中选择具有至少8个公共空间特征点的两幅图像，则两图像外极线几何约束用本质矩阵和像点坐标描述为如下式所示：

m^'TEm＝0         …(1)

其中m(m,n,c)^T和m'(m',n',c')^T代表同一空间特征点在两幅图像中对应像点在相机坐标系下的坐标，E是表征两图像中相机外方位参数的本质矩阵；

(4).用局部参数优化的方法对拍摄获得的每一组对应像点坐标构成的外极线几何约束式(1)中的本质矩阵进行如下参数优化过程：

1)首先将本质矩阵E进行奇异值分解，即其中对角矩阵U、V都是酉矩阵；再将本质矩阵E参数化表示，即，将U、V视为旋转矩阵，利用正交坐标系的三坐标轴上的旋转角，分别将U、V参数化表示成欧拉角分解的形式，即其中，R_u,R_x是绕x轴的旋转矩阵，R_v,R_y是绕y轴的旋转矩阵，R_w,R_z是绕z轴的旋转矩阵，R_z与具有互换性，将R_w和R_z两个旋转矩阵进行合并，将误差函数的自由参数降至五个，这时本质矩阵用旋转矩阵相乘表示为：

$E=R_u{R}_v{R}_w\hat{I}{R}_y^T{R}_x^T...\left(2\right)$

其中

$R_u=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos u& -\sin u\\ 0& \sin u& \cos u\end{array}\right],R_v=\left[\begin{array}{ccc}\cos v& 0& \sin v\\ 0& 1& 0\\ -\sin v& 0& \cos v\end{array}\right],R_w=\left[\begin{array}{ccc}\cos w& -\sin w& 0\\ \sin w& \cos w& 0\\ 0& 0& 1\end{array}\right]$

2)用高斯牛顿迭代法，迭代更新酉矩阵U、V，定义为误差函数，m_i＝(m_i,n_i,c_i)^T与m'_i＝(m′_i,n′_i,c′_i)^T是第i个空间特征点在两幅图像中的对应像点在相机坐标系下的坐标，定义角度参数变量θ＝(θ_u,θ_v,θ_w,θ_x,θ_y)^T，计算误差函数ε_i在角度初值θ＝(0,0,0,0,0)^T时的雅可比矩阵：

$J_i{|}_0=\frac{{\∂\mathrm{\ϵ}}_i}{\∂\mathrm{\θ}}{|}_{\mathrm{\θ}=0}=(n_i{c}_i^{\′},-m_i{c}_i^{\′},m_i{n}_i^{\′}-n_i{m}_i^{\′},c_i{n}_i^{\′},-c_i{m}_i^{\′})...\left(3\right)$

a)利用n个空间点的像面坐标和相机内参数，用归一化8点算法解算本质矩阵初值，并对该初值进行奇异值分解利用分解出的酉矩阵U和V进行像点坐标转换：即，将归一化后的对应像点坐标m_i(m_i,n_i,1)^T，m′_i(m′_i,n′_i,1)^T(i＝1,...,n,n≥8)左乘U^T、V^T表示为X'＝U^Tm′_i，X＝V^Tm_i并代入方程(1)，得到误差函数

b)将转换后的像点坐标X＝(X_i,Y_i,Z_i)^T与X'＝(X′_i,Y′_i,Z′_i)^T代入上面的误差函数，则方程简化为

c)将转化后的像点坐标X、X’代入(3)式得到如下雅可比矩阵J，由Jδ_θ＝-ε_i推得参数变量δ_θ＝-(J'J)^-1J'ε_i；

$J=\frac{{\∂\mathrm{\ϵ}}_i}{\∂\mathrm{\θ}}{|}_{\mathrm{\θ}=0}=(Y_i{Z}_i^{\′},-X_i{Z}_i^{\′},X_i{Y}_i^{\′},-Y_i{X}_i^{\′},Z_i{Y}_i^{\′},-Z_i{X}_i^{\′})...\left(4\right)$

d)用求解出来的δ_θ＝(δ_u,δ_v,δ_w,δ_x,δ_y)^T迭代更新原始的矩阵U、V得到新的U_new和V_new如下：

U_new＝UR_u(δ_u)R_v(δ_v)R_w(δ_w)         …(5)

V_new＝VR_x(δ_x)R_y(δ_y)

重复上述步骤，当||ε||满足预设的精度要求或收敛条件时，停止迭代，输出最终获得的本质矩阵

3)利用Hartley提出的Cheirality约束条件从本质矩阵E中解算出正确的两相机外方位参数中的平移向量T和旋转矩阵R，完成基于两图像的相对定向方法的流程，然后对其他具有至少8个公共特征点的两幅图像重复上述步骤1)和2)，完成所有图像间的相对定向。

2.如权利要求1所述的相对定向方法，其中所述空间特征点为具有特定的空间约束关系的编码点，直径为3mm。

3.如权利要求1所述的相对定向方法，其中所述空间特征点均匀分布在被测场中，且所有特征点不能共面，保证两相邻站位间的公共空间特征点数不能少于8个。

4.如权利要求1所述的相对定向方法，其中在对空间特征点进行拍摄时采用单相机系统或双相机系统，在所述单相机系统中，物体不动，相机运动，以从不同站位拍摄，在双相机系统中，物体运动而两相机不动。

5.如权利要求1所述的相对定向方法，其中所述首先对本质矩阵进行奇异值分解，然后对本质矩阵的两个相同特征值归一化处理，将本质矩阵分解出的特征值矩阵表示成对角矩阵最后再将本质矩阵参数化表示，简化误差函数的雅可比矩阵。

6.如权利要求1所述的相对定向方法，其中所述迭代停止条件为图像中的所有像点到对应外极线距离的均方根误差不再收敛或δ_θ≤10^-7。

7.如权利要求1所述的相对定向方法，其中所述误差函数的雅可比矩阵在每次迭代中的角度初值都取0，完成误差函数模型在0角度的局部迭代优化。

8.如权利要求4所述的相对定向方法，其中在所述双相机系统测量动态物体情况下，通过对贴有编码点的基准尺在测量场不同位置进行拍摄测量，完成两相机间的初始相对定向。

9.如权利要求4所述的相对定向方法，其中在所述单相机系统的情况下，完成所有图像间的相对定向后，选取包含测量场中编码点最多的图像作为世界坐标系，建立所有图像的相对方位转换链，将所有图像统一到一个世界坐标系下。

10.如权利要求1所述的相对定向方法，其中所述需要用经过标定的贴有编码点的基准尺作为尺度基准，来唯一确定外参数中平移向量T的比例因子。