CN104463884A - 最小化空间距离的线性交会测量方法 - Google Patents

Abstract

本发明属于图像处理、摄像测量、计算机视觉和三维重建等技术领域，特别涉及光学非接触测量中的交会测量领域。其技术方案是：最小化空间距离的线性交会测量方法，包括以下步骤：A.安装摄像机；B.摄像机标定；C.图像采集；D.图像处理；E.交会。本发明的交会测量方法能够同时保证测量的高精度和实时性，同样适用于光学测量、摄影测量、视觉测量或计算机视觉等测量领域，具有广泛的应用价值。

Description

最小化空间距离的线性交会测量方法

技术领域

本发明属于图像处理、摄像测量、计算机视觉和三维重建等技术领域，特别涉及光学非接触测量中的交会测量领域。

背景技术

摄像测量技术是近年发展起来的具有重大工程实用价值的新兴技术之一，具有非接触、高精度、可测点多、实时动态测量等显著特点。交会是摄像测量中的关键步骤之一。通过交会测量可以由图像二维信息恢复待测目标三维信息或三维结构性质，有些地方称之为结构计算。

交会可分为基于点的线线交会和基于线的面面交会。其中，线线交会方法研究和使用最多，下文中简称之交会。交会按是否线性分为线性与非线性两种；按是否双目分为双目交会和多目交会两种，其中双目交会是最基本最简单的交会方法。

已知的主要几种交会方法有如下五种：直接线性交会方法、平差物点的迭代方法、双目中垂线段中点法、双目Sampson近似方法、双目最优三角形方法。其中第一种属于线性方法，第二种属于迭代的非线性方法，后三种属于双目交会测量方法。

1.直接线性交会方法

直接线性交会方法是对共线方程进行直接线性变换后，用常用数值方法求解线性最小二乘问题。若为齐次方程形式，一般采用奇异值分解方法；若为非齐次方程形式，一般用正规方程求逆方法。

对第一幅图像，将共线方程x＝PX展开(x，P，X分别为像点齐次坐标向量、物点齐次坐标向量和3×4投影矩阵，斜体表示标量，加粗表示的为向量，^T符号意义为转置)，然后进行直接线性变换得：

$\left[\begin{array}{c}w{P}^{2T}-y{P}^{3T}\\ -w{P}^{1T}+x{P}^{3T}\\ {\mathrm{yP}}^{1T}-x{P}^{2T}\end{array}\right]x=0---\left(1\right)$

一幅图像中物点与像点的对应提供三个方程，其中只有两个是独立的。一个物点的齐次坐标有4个变量，非齐次坐标有3个变量，因此至少需要两幅图像点的对应才有可能交会。在像机参数与像点坐标已知条件下，在||X||＝1的约束下联立求解齐次线性方程组，称之为直接线性交会方法。

2.平差物点的迭代方法

平差物点的迭代方法是以线性交会结果作为初始值，以最小化图像上的重投影残差为目标函数，以物点(像机参数不准确时也优化像机参数)为待优化参数的方法。它具有明确的几何意义，且能提高交会的精度，被认为是目前为止最好的交会方法。该方法在射影几何空间中仍然适用。

3.双目中垂线段中点法

双目中垂线段中点法是比较纯粹的欧氏几何方法，目的是计算到两条射线的平均距离最小的物点坐标。其方法是将两条反投影射线的表达转换成在世界系中的表达，然后求解两射线之间的中垂线段中点，作为交会点。

左视图的反投影射线用点向式表示，它过左像机光心

$\tilde{C}=-M^{-1}{p}_4---\left(2\right)$

其中M和p₄分别为投影矩阵P的前三列和第四列，为光心非齐次坐标。射线方向在世界系中的归一化表示为：

d＝M^-1x/||M^-1x||       (3)

其中x为像点齐次坐标。对应的，右视图的反投影射线过光心且方向为两射线之间的距离为：

$e=\±(\tilde{C}-{\tilde{C}}^{\′})\·(d\×d^{\′})/\left|\right|d\×d^{\′}\left|\right|---\left(4\right)$

记d×d′＝[A B C]^T，重排列后记为g＝[C A B]^T。中垂线段的两端点为：

$X=\tilde{C}+\mathrm{\λd},X^{\′}={\tilde{C}}^{\′}+{\mathrm{\λ}}^{\′}{d}^{\′}---\left(5\right)$

其中 $\mathrm{\λ}=-[(\tilde{C}-{\tilde{C}}^{\′})\·(g\×d^{\′})]/[g\·(d\×d^{\′})]$ 和 ${\mathrm{\λ}}^{\′}=-[(\tilde{C}-{\tilde{C}}^{\′})\·(g\×d)]/[g\·(d\×d^{\′})].$ 交会点取为中垂线段中点 $(\tilde{X}+{\tilde{X}}^{\′})/2.$

双目中垂线段中点法最大的特点是，不需要求解方程，交会点能够从显示表达式直接得到。在大型光电经纬仪交会测量中常用此方法，因为两条反投影射线已事先统一在大地坐标系中进行表达。

4.双目Sampson近似方法

双目Sampson近似方法在测量误差相对于测量值很小时是有效的，它通过对两幅图像上测量像点的Sampson矫正，得到矫正点的一阶近似：

$\left[\begin{array}{c}\hat{x}\\ \hat{y}\\ {\hat{x}}^{\′}\\ {\hat{y}}^{\′}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ x^{\′}\\ y^{\′}\end{array}\right]-\frac{x^{\′T}\mathrm{Fx}}{\left|\right|\left[\begin{array}{cccc}{x}^{\′T}{F}_1& x^{\′T}{F}_2& F^{1T}x& F^{2T}x\end{array}\right]{\left|\right|}^2}\left[\begin{array}{c}{x}^{\′T}{F}_1\\ x^{\′T}{F}_2\\ F^{1T}x\\ F^{2T}x\end{array}\right]---\left(6\right)$

式中x＝[x y 1]^T和x′＝[x′ y′ 1]^T为两幅图像上测量像点的规范化齐次坐标，F为3×3两视图之间的基础矩阵。用矫正完的像点坐标进行线性交会得到交会点坐标。如果在每一幅图像上的矫正量很小(小于1pixel)，则该近似是准确的，且计算量小，但是矫正的点不会准确地满足对极几何关系，或者说两条反投影射线在空间中不一定相交。

5.双目最优三角形方法

双目最优三角形方法是在双目情形下，求解使重投影误差最小并且满足对极几何约束的物点。此方法的目标函数等价于最小化图像上的测量点到对应极线的距离平方和，通过将极平面用单参数表示，然后利用一元6次多项式求极值，最后得到的最优像点满足三角形相交。它是用非迭代算法实现全局最小化重投影残差的方法，如果高斯噪声模型的假定是正确的，可以证明此方法是最优的。

传统交会测量模式为先直接线性交会然后迭代平差物点，其中平差物点的方法为非线性迭代方法，以最小化图像上的重投影残差为目标函数，其初始值来自直接线性交会的线性解。该测量模式虽然能够保证精度，但是运算量大，不能保证实时，常用于事后处理。所以目前在实时测量时，通常只采用直接线性交会方法，但该方法却难以保证测量精度。

因此，发明一种具有上述两种方法的共同优点的方法，即既能保证测量的高精度，又能保证实时性的交会测量方法，将具有重大的意义和广泛的应用价值。。

发明内容

本发明的目的是：提供一种同时具有高精度和实时性的交会测量方法。

本发明的技术方案是：最小化空间距离的线性交会测量方法，包括以下步骤：

A.安装摄像机

安装用于交会测量的m台摄像机，m≥2；

B.摄像机标定

对摄像机内外参数进行标定，得到每台摄像机的投影矩阵P_j，j＝1，2，...，m；

P_j为3×4的投影矩阵，由公式P_j＝K[R T]求得，其中： $K=\left[\begin{array}{ccc}{f}_x& s& c_x\\ & f_y& c_y\\ & & 1\end{array}\right],$ 为3×3的摄像机内参数矩阵；f＝[f_x f_y]^T，为行列两个方向的等效焦距；c＝[c_x c_y]^T为主点；s为扭曲系数； $R=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],$ 为摄像机外参数旋转矩阵；A＝[α β γ]^T，为与所述摄像机外参数旋转矩阵对应的欧拉角向量；T＝[t_x t_y t_z]^T，为摄像机外参数平移向量，其物理意义为世界系原点在摄像机坐标系中的坐标；

C.图像采集

D.图像处理

D1.图像预处理；

D2.图像提取：

在每台摄像机拍摄的图像上提取测量像点坐标，得到m个测量像点的齐次坐标；分别记为x_j，x_j为3×1的列向量，x_j＝[x_j y_j 1]^T，T符号意义为转置；

D3.图像匹配：对于每一个物点，匹配得到其在所有图像上对应的图像点；

E.交会

以上述步骤处理得到的结果作为输入：

E1.构造齐次方程组：

$\left(\frac{1}{{\left|\right|M_j^{-1}{x}_j\left|\right|}^2}\right)M_j^T{\left[x_j\right]}_x{M}_j{M}_j^T{\left[x_j\right]}_x{P}_{j}X=0,j=\mathrm{1,2},...,m$

其中：M_j为投影矩阵P_j的前三列， $M=\left[\begin{array}{ccc}{P}_{11}& P_{12}& P_{13}\\ P_{21}& P_{22}& P_{23}\\ P_{31}& P_{32}& P_{33}\end{array}\right];$ T符号意义为转置，-1符号意义为矩阵求逆；[x_j]_x表示列向量x_j的反对称矩阵，大小3×3；表示向量的模，为标量；

E2.求解齐次方程组

记所述齐次方程组形式为：AX＝0，其中：A为3m×4的矩阵，X为4×1的列向量，方程右侧为3m×1的零向量；

对A用通用数值计算方法作奇异值分解，得到：A＝UDV^T，其中V为4×4的正交矩阵；

记v为V的最后一列，则：X＝v/v(4)，其中v(4)为v的第四个元素；

输出参数为：X＝[X Y Z 1]^T，为物点的齐次坐标。

本发明能够同时保证测量的高精度和实时性，精度与平差物点的迭代方法相当，计算量与直接线性交会方法相当。同样适用于光学测量、摄影测量、视觉测量或计算机视觉等测量领域，具有广泛的应用价值。

具体实施方式

最小化空间距离的线性交会测量方法，包括以下步骤：

A.安装摄像机

安装用于交会测量的m台摄像机，m≥2；

B.摄像机标定

对摄像机内外参数进行标定，得到每台摄像机的投影矩阵P_j，j＝1，2，...，m；

P_j为3×4的投影矩阵，由公式P_j＝K[R T]求得，其中： $K=\left[\begin{array}{ccc}{f}_x& s& c_x\\ & f_y& c_y\\ & & 1\end{array}\right],$ 为3×3的摄像机内参数矩阵；f＝[f_x f_y]^T，为行列两个方向的等效焦距；c＝[c_x c_y]^T为主点；s为扭曲系数； $R=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],$ 为摄像机外参数旋转矩阵；A＝[α β γ]^T，为与所述摄像机外参数旋转矩阵对应的欧拉角向量；T＝[t_x t_y t_z]^T，为摄像机外参数平移向量，其物理意义为世界系原点在摄像机坐标系中的坐标；

C.图像采集

D.图像处理

D1.图像预处理；

D2.图像提取：

在每台摄像机拍摄的图像上提取测量像点坐标，得到m个测量像点的齐次坐标；分别记为x_j，x_j为3×1的列向量，x_j＝[x_j y_j 1]^T，T符号意义为转置；

D3.图像匹配：对于每一个物点，匹配得到其在所有图像上对应的图像点；

E.交会

以上述步骤处理得到的结果作为输入：

E].构造齐次方程组：

$\left(\frac{1}{{\left|\right|M_j^{-1}{x}_j\left|\right|}^2}\right)M_j^T{\left[x_j\right]}_x{M}_j{M}_j^T{\left[x_j\right]}_x{P}_{j}X=0,j=\mathrm{1,2},...,m$

其中：M_j为投影矩阵P_j的前三列， $M=\left[\begin{array}{ccc}{P}_{11}& P_{12}& P_{13}\\ P_{21}& P_{22}& P_{23}\\ P_{31}& P_{32}& P_{33}\end{array}\right];$ T符号意义为转置，-1符号意义为矩阵求逆；[x_j]_x表示列向量x_j的反对称矩阵，大小3×3；表示向量的模，为标量；

E2.求解齐次方程组

记所述齐次方程组形式为：AX＝0，其中：A为3m×4的矩阵，X为4×1的列向量，方程右侧为3m×1的零向量；

对A用通用数值计算方法作奇异值分解，得到：A＝UDV^T，其中V为4×4的正交矩阵；

记v为V的最后一列，则：X＝v/v(4)，其中v(4)为v的第四个元素；

输出参数为：X＝[X Y Z 1]^T，为物点的齐次坐标。

关于所构造齐次方程合理性的说明如下：

仅考虑一个摄像机所联系的物点与像点的对应。设摄像机的投影矩阵为P，像点坐标为x，光心到像点的反投影射线方向在世界系中的归一化表示为：

d＝M^-1x/||M^-1x||         (7)

其中M为投影矩阵P的前三列；

摄像机光心到物点的有向矢量为：

$v=\tilde{X}-\tilde{C}---\left(8\right)$

其中和分别为物点和光心的非齐次坐标。则物点到反投影射线的距离矢量为：

$\mathrm{\Δd}=(\tilde{X}-\tilde{C})-\left[d^T(\tilde{X}-\tilde{C})\right]d=(I-{\mathrm{dd}}^T)(\tilde{X}-\tilde{C})---\left(9\right)$

其中I为3×3的单位矩阵；

将式(2)的光心表达式p₄为P的第四列；将光心表达式代入式(9)中，得到式(9)的第二项为：

$\tilde{X}-\tilde{C}=\tilde{X}+M^{-1}{p}_4=M^{-1}\left[\begin{array}{cc}M& p_4\end{array}\right]\left[\begin{array}{c}\tilde{X}\\ 1\end{array}\right]=M^{-1}\mathrm{PX}---\left(10\right)$

将式(7)代入式(9)中，得到式(9)的第一项为

I-dd^T＝-[d]_x[d]_x＝-[M^-1x]_x[M^-1x]_x/||M^-1x||²     (11)

利用[t]_xM＝M^-T[M^-1t]_x得

[M^-1x]_x＝M^T[x]_xM      (12)

将式(12)代入式(11)，然后和式(10)一起代入式(9)得(省略负号)

$\mathrm{\Δd}=\left(\frac{1}{{\left|\right|M^{-1}x\left|\right|}^2}\right)M^T{\left[x\right]}_x{\mathrm{MM}}^T{\left[x\right]}_x\mathrm{PX}---\left(13\right)$

当没有误差时，光心到物点的有向矢量与反投影射线重合，物点到反投影射线的距离Δd为零；当存在误差时，该距离Δd不为零；

给定物点和像点之间的一个对应关系，利用式(13)建立的一个方程组由三个方程组成，它们分别对应空间中三个正交方向上的距离残差。给定m组(m＞＝2)物点和像点之间的对应关系，利用式(13)可建立m个方程组，共3m个方程。使用最小二乘方法求解该联合线性方程组的物理意义明确，即最小化物点到所有3m个反投影射线的空间距离平方和。因此，只需一步解算即可得到使空间距离残差最小的物点三维坐标解。

Claims (1)

1.最小化空间距离的线性交会测量方法，包括以下步骤：

A.安装摄像机

安装用于交会测量的m台摄像机，m≥2；

B.摄像机标定

对摄像机内外参数进行标定，得到每台摄像机的投影矩阵P_j，j＝1，2，...，m；

P_j为3×4的投影矩阵，由公式P_j＝K[R T]求得，其中： $K=\left[\begin{array}{ccc}{f}_x& s& c_x\\ & f_y& c_y\\ & & 1\end{array}\right],$ 为3×3的摄像机内参数矩阵；f＝[f_x f_y]^T，为行列两个方向的等效焦距；c＝[c_x c_y]^T为主点；s为扭曲系数； $R=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],$ 为摄像机外参数旋转矩阵；A＝[α β γ]^T，为与所述摄像机外参数旋转矩阵对应的欧拉角向量；T＝[t_x t_y t_z]^T，为摄像机外参数平移向量，其物理意义为世界系原点在摄像机坐标系中的坐标；

C.图像采集

D.图像处理

D1.图像预处理；

D2.图像提取：

在每台摄像机拍摄的图像上提取测量像点坐标，得到m个测量像点的齐次坐标；分别记为x_j，x_j为3×1的列向量，x_j＝[x_j y_j 1]^T，T符号意义为转置；

D3.图像匹配：对于每一个物点，匹配得到其在所有图像上对应的图像点；

E.交会

以上述步骤处理得到的结果作为输入：

E1.构造齐次方程组：

$\left(\frac{1}{{\left|\right|M_j^{-1}{x}_j\left|\right|}^2}\right)M_j^T{\left[x_j\right]}_x{M}_j{M}_j^T{\left[x_j\right]}_x{P}_{j}X=0,j=\mathrm{1,2},...,m$

其中：M_j为投影矩阵P_j的前三列， $M=\left[\begin{array}{ccc}{P}_{11}& P_{12}& P_{13}\\ P_{21}& P_{22}& P_{23}\\ P_{31}& P_{32}& P_{33}\end{array}\right];$ T符号意义为转置，-1符号意义为矩阵求逆；[x_j]_x表示列向量x_j的反对称矩阵，大小3×3；表示向量的模，为标量；

E2.求解齐次方程组

记所述齐次方程组形式为：AX＝0，其中：A为3m×4的矩阵，X为4×1的列向量，方程右侧为3m×1的零向量；

对A用通用数值计算方法作奇异值分解，得到：A＝UDV^T，其中V为4×4的正交矩阵；

记v为V的最后一列，则：X＝v/v(4)，其中v(4)为v的第四个元素；

输出参数为：X＝[X Y Z 1]^T，为物点的齐次坐标。