CN103617649A - 一种基于相机自标定技术的河工模型地形测量方法 - Google Patents

Abstract

本发明涉及一种基于相机自标定技术的河工模型地形测量方法，属于河工模型测量技术领域，该方法包括：制作和安装标定板，制作和撒布辅助标定片，对河工模型进行拍照，提取图片的特征点，将所有图片进行两两配对，对图片对中的特征点以特征距离进行匹配；剔除不合理的匹配同名点对，选择匹配同名点对最多的一对图片解算相机参数和匹配同名点对的世界坐标：使用最优化解法微调所有已解图片的相机参数，标定河工模型稀疏点云，生成河工模型稠密点云，根据得到的河工模型稠密点云的三维坐标数据生成河工模型地形。本发明具有精度较高、使用方便、成本低的优点。

Description

一种基于相机自标定技术的河工模型地形测量方法

技术领域

本发明属于河工模型测量技术领域，特别涉及一种基于相机自标定技术的河工模型地形测量方法。

背景技术

在重要的水利水电工程规划、设计和建设过程中，均需要进行河工模型试验(研究河床演变与河道整治)。对试验前后的模型地形进行测量，就可得到不同水流条件下河床的冲淤变化以及泥沙的冲淤量。早期的测量方法主要以使用水平仪人工测量为主，这类方法费时费力，而且测量点相对较少，得到的模型地形分辨率较低。随着自动测控技术的发展，以声光电原理为基础的河工模型地形测量方法大量出现，例如使用超声传感器或者激光测距仪测量地形。这类方法每次测量只能得到一个数据点，并且需要在各类测距仪上附加许多精密机械设备以控制其姿态及位置，测量精度依赖于这些机械结构的加工与安装调试精度，整个系统的成本较大，实际使用时的测量时间很长。最近二十年来，机器视觉技术发展迅速，出现了使用多相机同时拍摄河工模型图片，通过相机标定及图形图像技术从图像中重构河工模型地形的方法。这类方法需要多个相机，并且事先标定每个相机的参数，在测量过程中相机的位置、焦距等所有参数均不能变化，实际使用中设备架设好后便不能移动，而单套设备的测量范围较小，当需要测量规模较大的河工模型时，需要布置多套设备。因此，急需发明一种使用方便、精度较高且成本较低的河工模型地形测量方法。

近年来，不事先标定相机参数，而是直接根据所拍摄的景物根据一定原理来标定相机参数并得到物体空间坐标的自标定方法成为热点。自标定方法的理论基础是对极几何和摄像机投影模型，关键技术是特征点提取和匹配。对极几何和摄像机投影模型是机器视觉技术的理论基础，已经非常成熟。对于特征点提取和匹配，Lowe于2004年提出了SIFT算法（Lowe D.2004.Distinctive image features from scale-invariant key-points.InternationalJournal of Computer Vision,60(2),91–110），SIFT算法的实质是将计算分布在不同的尺度空间上，从而能在不同尺度上查找特征点。SIFT能在目标旋转、缩放、平移时保持不变，并能克服光照、遮挡以及图像噪声的影响。SIFT算法对其找到的每个特征点都附有一个描述符，描述其图像特征，便于在不同图像中对特征点进行匹配。当两张图片从不同位置和角度拍摄空间中的同一点时，不同图片中此点图像的SIFT描述符应该十分接近，因此可根据SIFT描述符对不同图片中的特征点进行匹配，匹配好的点对称为匹配同名点对，匹配同名点对是空间中同一点在不同图片中的图像。由于SIFT算法能够得到数量众多的特征点，所以需要快速高效的同名点匹配算法，最近邻结点kd树法（Arya S,et al.1998.Anoptimal algorithm for approximate nearest neighbor searching fixed dimensions.Journal of theACM,45(6),891–923）对特征点进行有效组织和排列，极大加快了匹配速度。根据对极几何，匹配同名点对的坐标满足极线约束，即匹配同名点对在各自图像中的坐标与基础矩阵F相乘之积为0。Fischler等提出了RANSAC算法，根据两张图片的匹配同名点对序列的图像坐标快速准确计算出基础矩阵F（Fischler MA,Bolles RC.1981.Random sampleconsensus:A paradigm for model fitting with applications to image analysis and automatedcartography.Communications of the ACM,24(6),381–395）。得到基础矩阵后，即可根据摄像机投影模型和相机参数计算出匹配同名点对代表的空间一点的坐标。自标定方法由于不需要精确设置相机位置、不需要事先标定相机参数等优势，在机器人视觉、导弹制导、水下地形测量、灾害现场重构、工程施工控制以及文物保护与研究等领域得到了广泛应用。

发明内容

本发明的目的是为克服已有河工模型地形测量方法的不足之处,提供一种基于相机自标定技术的河工模型地形测量方法。本发明利用相机自标定技术，操作者只需在河工模型上摆放标定板，之后手持普通照相机对准河工模型从不同角度拍照，将图片进行处理和坐标标定，即可获得准确的河工模型地形。具有精度较高、使用方便、成本低的优点。

为实现上述目的，本发明提出一种基于相机自标定技术的河工模型地形测量方法，其特征在于，该方法包括以下步骤：

步骤A：制作和安装标定板：制作由标定平板和固定支架组成的标定板，标定平板的表面设置有四个黑白相间的矩形格图案，固定支架与标定平板的底面固定，并能调节标定平板的高度及水平位置；记录矩形格图案的矩形格长、短边的数值；将标定板放置在待测河工模型中任意一不过水的位置上，将标定板的平板调整水平；调整好后标定板在整个测量过程中固定不动；

步骤B：制作和撒布辅助标定片：制作大量表面设置有黑白格的辅助标定片；将大量辅助标定片随机撒布在待测河工模型河床表面；

步骤C：对河工模型进行拍照：采用普通数码相机对河工模型进行拍照，使相机的镜头对准河工模型的中心，相机高度需使图片拍摄到河工模型整个范围；按照河工模型四周边及两对角线组成6条拍摄路径，依次移动相机进行拍摄，每条路径上拍摄至少10张图片/每米，每张图片中均包含至少3个辅助标定片；

步骤D：提取图片的特征点：将所得图片全部输入计算机，使用SIFT算法提取每张图片中的特征点，得到每个特征点在图片中的位置，并对每个特征点给出一个128个数的描述符；

步骤E：将所有图片进行两两配对，排列出所有图片对并记录每个图片对的序号；

步骤F：按顺序提取步骤E所得图片对的序号，对该图片对中的特征点以特征距离进行匹配；

步骤G：剔除不合理的匹配同名点对：根据步骤F中得到图片对的匹配同名点对，使用RANSAC算法计算图片对的基础矩阵；得到基础矩阵后，使用极线约束将与极线约束的误差大于设定值b的匹配同名点对剔除，当最终匹配同名点对的数目为设定值c以上时，接受此图片对中的所有匹配同名点对，将图片序号、匹配同名点对在各自图片中的位置以及各自的描述符记录下来。当匹配同名点对的数目少于c对时，认为该图片对无合理匹配同名点对；

步骤H：判断是否有未计算的图片对，若有，提取图片对重复步骤F到G，若无，继续进行步骤I；

步骤I：选择匹配同名点对最多的一对图片（I，J）解算相机参数和匹配同名点对的世界坐标：

步骤J：在未解图片中找到含有最多已解点的图片L；

步骤K：从已解点的世界坐标和已解点在图片L中的投影图像坐标（x_L，y_L），根据投影方程(5)使用最优化解法计算出图片L的相机参数，再采用已解图片的相机参数根据方程组（6）解算图片L中未解匹配同名点对的世界坐标，得到新的已解点；

步骤L：使用最优化解法微调所有已解图片的相机参数，使得已解点的总体投影误差最小；

步骤M：判断是否含有未解图片，若有，则转回步骤J，如没有，则进行步骤N；

步骤N：标定河工模型稀疏点云：所有图片均已解算出相机参数，所有匹配同名点对均已确定世界坐标，所有已解点组成未标定的河工模型稀疏点云；选取包含完整标定板图像的任意三张图片，将标定板中四个黑白相间的矩形格的9个角点的世界坐标代入未标定的河工模型稀疏点云，将该河工模型稀疏点云中的相机参数和已解点坐标均标定为真实世界坐标，得到已标定河工模型稀疏点云；

步骤O：生成河工模型稠密点云：以已解点为基础，在已解点附近寻找其它匹配同名点并确定世界坐标，最后得到河工模型稠密点云；

步骤P：根据河工模型稠密点云生成河工模型地形：利用河工模型稠密点云的三维坐标数据生成河工模型地形。

本发明具有如下优点:不需要购买价格昂贵的传感器，不需要对仪器进行精确位置固定和标定，手持普通相机拍摄图片即可得到精确的河工模型地形。相较于传统方法，本发明成本低廉、使用极其方便。增加图片数量时，能够得到毫米级精度的河工模型地形。

附图说明

图1为本发明方法流程框图；

图2为本发明中的标定板；

图3为本发明中的辅助标定片；

图4为本发明中相机拍摄路径；

图5为本发明的实施例河工模型地形；

图6为本发明与已有技术的测量精度对比。

具体实施方式

下面结合附图及实施例详细说明本发明技术方案中所涉及的各个细节问题。应指出的是,所描述的实施例仅旨在便于对本发明的理解,而对其不起任何限定作用。

本发明的一种基于相机自标定技术的河工模型地形测量方法，如图1步骤如下：

步骤A：制作和安装标定板：制作由标定平板和固定支架组成的标定板，标定平板的表面设置有四个黑白相间的矩形格图案，固定支架与标定平板的底面固定，并能调节标定平板的高度及水平位置；记录矩形格图案的矩形格长、短边的数值（作为步骤N中设置矩形格角点坐标的依据）；将标定板放置在待测河工模型中任意一不过水的位置上，将标定板的平板调整水平；调整好后标定板在整个测量过程中固定不动；

步骤B：制作和撒布辅助标定片：制作大量表面设置有黑白格的辅助标定片（辅助标定片边长在2～5厘米为宜，具体大小以在步骤C拍摄的每张图片中能清晰分辨出每个辅助标定片的黑白格为准，辅助标定片的作用为使后续SIFT算法能准确定位辅助标定片的中心点，增加测量精度）；将大量辅助标定片随机撒布在待测河工模型河床表面（撒布尽量均匀，各辅助标定片之间密度需保证所拍摄的每张图片中都包含至少3个辅助标定片）；

步骤C：对河工模型进行拍照：采用普通数码相机对河工模型进行拍照，使相机的镜头对准河工模型的中心，相机高度需使图片拍摄到河工模型整个范围（当河工模型规模很大时，可对河工模型进行分段拍摄，每段大小宜控制在5×5米以下）；按照河工模型四周边及两对角线组成6条拍摄路径，依次移动相机进行拍摄，每条路径上拍摄至少10张图片/每米，每张图片中均包含至少3个辅助标定片（当辅助标定片密度不满足要求时，可增加撒布辅助标定片）；

步骤D：提取图片的特征点：将所得图片全部输入计算机，使用SIFT算法提取每张图片中的特征点（提取出的特征点包含了辅助标定片的中心点），得到每个特征点在图片中的位置，并对每个特征点给出一个128个数的描述符m_i（i为1到128之间的整数，如果同一物体被拍进不同图片中，由于物体自身结构的一致性，不同图片中反映物体上同一点的那些特征点应具有相近的描述符m_i）；

步骤E：将所有图片进行两两配对，排列出所有图片对并记录每个图片对的序号；

步骤F：按顺序提取步骤E所得图片对的序号，对该图片对中的特征点以特征距离进行匹配；具体过程为：设匹配图片对序号为（I，J），定义图片I中特征点p和图片J中的特征点q之间的特征距离d_pq为：

$d_{\mathrm{pq}}=\sqrt{\underset{i=1}{\overset{128}{\mathrm{\Σ}}}{(m_{i,\mathrm{Ip}}-m_{i,\mathrm{Jq}})}^2}---\left(1\right)$

其中m_i,Ip和m_i,Jq分别为图片I中特征点p和图片J中的特征点q的描述符。对图片I中特征点p，遍历图片J中的所有特征点计算（1）式，得到最小和次小的特征距离值，若最小和次小特征距离值之比小于一设定值a（a取值在0.3到0.8之间，取值越小所得匹配点越少，匹配精度越高，相反则所得点越多精度越低），则认为该最小特征距离值对应的特征点和图片I中特征点p匹配，称为匹配同名点对（认为匹配的特征点是河工模型上的同一个点在不同图片中的反映），若大于a，则认为图片J中无特征点与图片I中特征点p匹配；

步骤G：剔除不合理的匹配同名点对：根据步骤F中得到图片对的匹配同名点对，使用RANSAC算法计算图片对的基础矩阵；得到基础矩阵后，使用极线约束将与极线约束的误差大于设定值b的匹配同名点对剔除（容许误差b的范围在0.06到0.005之间，值设置得越大，最终得到点越多，但误差越大，相反则点越少精度越高），当最终匹配同名点对的数目为设定值c以上时（c至少大于8，设置数值越大精度越高），接受此图片对中的所有匹配同名点对，将图片序号、匹配同名点对在各自图片中的位置以及各自的描述符mi记录下来。当匹配同名点对的数目少于c对时，认为该图片对无合理匹配同名点对；

步骤H：判断是否有未计算的图片对，若有，提取图片对重复步骤F到G，若无，继续进行步骤I；

步骤I：选择匹配同名点对最多的一对图片（I，J）解算相机参数和匹配同名点对的世界坐标：一对匹配同名点对中，在图片I中的图像坐标为（x_I，y_I），在图片J中的图像坐标为（x_J，y_J），该匹配同名点对对应的河工模型上同一点的世界坐标为（X，Y，Z），图像坐标与世界坐标之间满足投影关系，以图片I为例：

$\left\{\begin{array}{c}{\left[\begin{array}{ccc}{X}_I& Y_I& Z_I\end{array}\right]}^T=A_I\left[\begin{array}{cc}{R}_I& t_I\end{array}\right]\left[\begin{array}{ccc}X& Y& Z\end{array}\right]\\ \left[\begin{array}{c}{x}_c\\ y_c\end{array}\right]=\left[\begin{array}{c}{X}_I/Z_I\\ Y_I/Z_I\end{array}\right]\\ \left[\begin{array}{c}{x}_I\\ y_I\end{array}\right]=(1+k_1{r}^2+k_2{r}^4)\left[\begin{array}{c}{x}_c\\ y_c\end{array}\right]\end{array}\right.---\left(2\right)$

其中A_I为相机内参数，R_I为旋转矩阵，t_I为空间平移矢量，（X_I，Y_I，Z_I）为在以图片中心为坐标系原点、图片长短边为x和y轴的坐标系下，世界坐标为（X，Y，Z）的点的坐标；（x_c，y_c）为图像无畸变情况下点（X，Y，Z）在图片中所成像点的图像坐标，k₁、k₂为径向畸变参数，r²=x_c ²+y_c ²，（x_I，y_I）为点（X，Y，Z）在图片中所成像点的实际图像坐标；A_I、R_I、t_I、k₁、k₂均为待标定参数；A_I和R_I的具体形式为：

$A_I=\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]---\left(3\right)$

$\left\{\begin{array}{c}{R}_I=R_x\left({\mathrm{\θ}}_x\right)R_y\left({\mathrm{\θ}}_y\right)R_z\left({\mathrm{\θ}}_z\right)\\ R_x\left({\mathrm{\θ}}_x\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_x& {-\sin \mathrm{\θ}}_x\\ 0& {\sin \mathrm{\θ}}_x& {\cos \mathrm{\θ}}_x\end{array}\right]\\ R_y\left({\mathrm{\θ}}_y\right)=\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_y& 0& {\sin \mathrm{\θ}}_y\\ 0& 1& 0\\ {-\sin \mathrm{\θ}}_y& 0& {\cos \mathrm{\θ}}_y\end{array}\right]\\ R_z\left({\mathrm{\θ}}_z\right)=\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_z& {-\sin \mathrm{\θ}}_z& 0\\ {\sin \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_z& 0\\ 0& 0& 1\end{array}\right]\end{array}\right.---\left(4\right)$

其中f为相机焦距，θ_x、θ_y和θ_z为绕x、y、z轴旋转角度，t_I中包含三个未知数，加上k₁、k₂，投影关系(2)中一共有9个相机参数；将这9个相机参数组成一个9维矢量Θ_I，定义投影算子P

[x_I y_I]＝P(Θ_I,X,Y,Z)             （5）

算子P的计算过程为公式(2)。因此，对于一个匹配同名点对，可得到方程组

$\left\{\begin{array}{c}\left[\begin{array}{cc}{x}_I& y_I\end{array}\right]=P({\mathrm{\Θ}}_I,X,Y,Z)\\ \left[\begin{array}{cc}{x}_J& y_J\end{array}\right]=P({\mathrm{\Θ}}_J,X,Y,Z)\end{array}\right.---\left(6\right)$

从方程组（6）中根据世界坐标的X、Y、Z分别相等能得到3个独立方程，包含18个未知数。对于图片对（I，J）中的所有匹配同名点对，相机参数Θ_I和Θ_J均是相等的。使用最优化解法（例如Levenberg-Marquardt法，Nocedal J,Wright SJ.1999.Springer series inoperations research.Numerical optimization.New York:Springer），代入匹配同名对的图像坐标，解算出每张图片的相机参数ΘI和ΘJ（已解算出相机参数的图片称为已解图片）。得到两张已解图片的所有相机参数后，根据已知的Θ_I、Θ_J、（x_I，y_I）和（x_J，y_J）以及方程组(6)反算出此图片对中所有匹配同名点对所反映的河工模型上的点的世界坐标（X，Y，Z），这些已算出世界坐标的匹配同名点对称为已解点，其它尚未确定世界坐标的匹配同名点对称为未解点，未确定相机参数的图片称为未解图片。

步骤J：在未解图片中找到含有最多已解点的图片L；

步骤K：从已解点的世界坐标和已解点在图片L中的投影图像坐标（x_L，y_L），根据投影方程方程(5)使用最优化解法计算出图片L的相机参数Θ_L，再采用已解图片的相机参数根据方程组（6）解算图片L中未解匹配同名点对的世界坐标，得到新的已解点；

步骤L：使用最优化解法微调所有已解图片的相机参数，使得已解点的总体投影误差最小；

步骤M：判断是否含有未解图片，若有，则转回步骤J，如没有，则进行步骤N；

步骤N：标定河工模型稀疏点云：所有图片均已解算出相机参数，所有匹配同名点对均已确定世界坐标，所有已解点组成未标定的河工模型稀疏点云；选取包含完整标定板图像的任意三张图片，将标定板中四个黑白相间的矩形格的9个角点的世界坐标（根据标定板的实际尺寸给出，由于标定板水平放置，垂直坐标均为一已知值）代入未标定的河工模型稀疏点云，将该河工模型稀疏点云中的相机参数和已解点坐标均标定为真实世界坐标，得到已标定河工模型稀疏点云；

步骤O：生成河工模型稠密点云：以已解点为基础，在已解点附近寻找其它匹配同名点并确定世界坐标，最后得到河工模型稠密点云。

步骤P：根据河工模型稠密点云生成河工模型地形：利用河工模型稠密点云的三维坐标数据生成河工模型地形（可采用具有网格生成和插值功能的软件，例如Tecplot360）。

下面以黄河科学院黄河内蒙古段河工模型地形测量为实施例说明本发明的步骤，待测河工模型长约为20米，宽约为4米，在长方向上将河工模型分为5段区域，每段区域范围约为4×4米。本实施例的基于相机自标定技术的河工模型地形测量方法,如图1所示,具体包括以下步骤:

步骤A：制作和安装标定版：标定板由木质平板11、A4纸12、4根8毫米螺杆13和8个螺母14制成,如图2所示；在木质平板11四角打直径为10毫米的孔，用4根螺杆13穿过4个孔，并在孔的上下安装螺母14，将螺杆13和木质平板11固定为一体；在一张A4纸12上打印四个黑白相间的矩形格图案，每个矩形格的长边为130毫米，短边为90毫米，将打印好的A4纸12平贴在木质平板11朝上的平面；将标定板放置在待测河工模型中西柳沟汇入黄河主河道入汇口的不过水河滩上，将A4纸12的长边与黄河主河道大致摆放平行；使用水平尺，转动螺母14将标定板的木质平板11调整水平，调整好后标定板在整个测量过程中不能移动；

步骤B：制作和撒布辅助标定片2：在20张A4纸上双面打印正方形黑白格，每小正方形格边长为20毫米，裁剪为如图3所示辅助标定片2，其作用为使后续SIFT算法能准确定位辅助标定片2的中心点21，增加测量精度；将大量辅助标定片2略微湿润，随机撒布在待测河工模型河床表面，撒布尽量均匀；

步骤C：对河工模型进行拍照：手持松下DMC-LX5相机对每段区域31进行拍摄，拍摄高度和位置要对准本段区域的中心，同时兼顾尽量使图片拍摄到更大范围；按照如图4所示6条拍摄路径32、33、34、35、36、37移动相机进行拍摄，每条路径上均匀拍摄40张图片；每段区域31共拍摄240张图片（每张图片中都包含至少3个辅助标定片3，当辅助标定片3密度不满足要求时，按步骤B增加），整个河工模型总共拍摄1200张图片；

步骤D：使用SIFT算法对每张图片提取特征点：使用Dell Precision T7610型工作站处理数据(由于图片量大，工作站内存为8G)；将所得1200张图片输入工作站，使用SIFT算法提取每张图片中的特征点，每张图片中平均能提取到数千至一万多个特征点（提取出的特征点包含了辅助标定片2的中心点）；SIFT算法所得结果中包含每个特征点在图片中的位置，还对每个特征点给出一个128个数的描述符mi（i为1到128之间的整数），描述符mi描述特征点附近的灰度梯度特征；

步骤E：将所有图片进行两两配对，排列出所有图片对，并记录在数组Sequence中，数组Sequence每列记录一个图片对的序号，每列第1、2行分别为第1、2张图片的序号；数组Sequence一共有719400列；

步骤F：按顺序提取数组Sequence中每列记录的图片对，将该图片对中的特征点进行匹配：设匹配图片对序号为（I，J），对图片I中的特征点p，遍历图片J中的所有特征点计算（1）式，得到最小和次小的特征距离值，若最小和次小特征距离值之比小于a=0.6，则认为该最小特征距离值对应的特征点和图片I中的特征点p匹配，若大于0.6，则认为图片J中无特征点与图片I中的特征点p匹配（具体计算时使用最近邻结点kd树法）；

步骤G：使用极线约束删除不合理的匹配同名点对：根据步骤F中得到图片对（I，J）的匹配同名点对，使用RANSAC算法计算图片对的基础矩阵；进行每次RANSAC迭代时，容许误差设置为0.6%倍图片长边像素值（RANSAC迭代时的容许误差范围为0.3%到1%之间，值设置得越大，最终得到点越多，但误差越大，相反则点越少精度越高）。RANSAC算法给出基础矩阵后，将与极线约束的误差大于b=0.01的匹配同名点对剔除。当最终匹配同名点对多于c=20对时，接受此对图片中的所有匹配同名点对，将图片序列、匹配同名点对在各自图片中的位置以及各自的描述符mi记录在文件Candidate.txt中。当匹配同名点对的数目少于20时，认为图片对无合理匹配同名点，删除所有匹配同名点对的记录；

步骤H：判断数组Sequence中是否有未计算的图片对，若有，提取未计算的图片对重复步骤F到G，若无，继续进行步骤I；

步骤I：在文件Candidate.txt中选择匹配同名点对最多的一个图片对（I，J）解算相机参数和匹配同名点对的世界坐标：使用Levenberg-Marquardt法，代入图片对（I，J）中的所有匹配同名对的图像坐标，根据方程组（6）解算出每张图片的相机参数Θ_I和Θ_J；得到两张图片的相机参数后，根据已知的Θ_I、Θ_J、（x_I，y_I）和（x_J，y_J）以及方程组(6)反算出此对图片中所有匹配同名点对所反映的河工模型上的点的世界坐标（X，Y，Z），将其记录在文件Points.txt中；

步骤J：在未解图片中找到含有最多已解点的图片L；

步骤K：从已解点的坐标和已解点在图片L中的投影图像坐标（x_L，y_L），根据投影方程(5)使用Levenberg-Marquardt法计算出图片L的相机参数Θ_L，再根据已解图片的相机参数采用方程组（6）解算图片L中未解匹配同名点对的世界坐标，得到新的已解点；记录在文件Points.txt中；

步骤L：使用Levenberg-Marquardt法微调所有已解图片的相机参数，使得已解点的总体投影误差最小；

步骤M：判断文件Candidate.txt中是否含有未解图片，若有，则转回步骤J，如没有，则进行步骤N；

步骤N：标定河工模型稀疏点云：所有图片均已解算出相机参数，所有匹配同名点对均已确定世界坐标，此时得到354701个已解点，已解点组成未标定的河工模型稀疏点云；选取包含完整标定板图像的任意三张图片，将标定板中的9个角点121、122、123、124、125、126、127、128、129（包含在特征点中）的世界坐标分别为（X,Y,Z）=（0,0,0）、（130,0,0）、（260,0,0）、（0,90,0）、（0,180,0）、（130,90,0）、（130,180,0）、（260,90,0）和（260,180,0）代入未标定的河工模型稀疏点云，将该河工模型稀疏点云中的相机参数和已解点坐标均标定为以毫米为单位的世界坐标，得到已标定河工模型稀疏点云。

步骤O：生成河工模型稠密点云：将已标定河工模型稀疏点云与匹配同名点对所在图像序号、图像坐标以及描述符输入开源的PMVS程序（可在网上免费下载和使用），该程序以已解点为基础，在其附近寻找其它匹配同名点并确定世界坐标，最后得到河工模型稠密点云。河工模型稠密点云中共有12584351个已知世界坐标的点，PMVS将河工模型稠密点云中每个点的三维坐标数据记录在mode.pset文件中。

步骤P：根据河工模型稠密点云生成河工模型地形：将mode.pset文件中的点的三维坐标数据复制到文件mode.txt文件中，保存文件。使用软件Tecplot360直接读取mode.txt文件，利用Tecplot360的网格生成功能直接生成河工模型地形。

最终生成的河工模型地形如图5所示。图6为测量地形与实际地形对比，可见本方法所得最终河工模型地形精度在厘米量级，满足一般河工模型试验的要求。

Claims (3)

1.一种基于相机自标定技术的河工模型地形测量方法，其特征在于，该方法包括以下步骤：

步骤A：制作和安装标定板：制作由标定平板和固定支架组成的标定板，标定平板的表面设置有四个黑白相间的矩形格图案，固定支架与标定平板的底面固定，并能调节标定平板的高度及水平位置；记录矩形格图案的矩形格长、短边的数值；将标定板放置在待测河工模型中任意一不过水的位置上，将标定板的平板调整水平；调整好后标定板在整个测量过程中固定不动；

步骤B：制作和撒布辅助标定片：制作大量表面设置有黑白格的辅助标定片；将大量辅助标定片随机撒布在待测河工模型河床表面；

步骤C：对河工模型进行拍照：采用普通数码相机对河工模型进行拍照，使相机的镜头对准河工模型的中心，相机高度需使图片拍摄到河工模型整个范围；按照河工模型四周边及两对角线组成6条拍摄路径，依次移动相机进行拍摄，每条路径上拍摄至少10张图片/每米，每张图片中均包含至少3个辅助标定片；

步骤D：提取图片的特征点：将所得图片全部输入计算机，使用SIFT算法提取每张图片中的特征点，得到每个特征点在图片中的位置，并对每个特征点给出一个128个数的描述符；

步骤E：将所有图片进行两两配对，排列出所有图片对并记录每个图片对的序号；

步骤F：按顺序提取步骤E所得图片对的序号，对该图片对中的特征点以特征距离进行匹配；

步骤G：剔除不合理的匹配同名点对：根据步骤F中得到图片对的匹配同名点对，使用RANSAC算法计算图片对的基础矩阵；得到基础矩阵后，使用极线约束将与极线约束的误差大于设定值b的匹配同名点对剔除，当最终匹配同名点对的数目为设定值c以上时，接受此图片对中的所有匹配同名点对，将图片序号、匹配同名点对在各自图片中的位置以及各自的描述符记录下来。当匹配同名点对的数目少于c对时，认为该图片对无合理匹配同名点对；

步骤H：判断是否有未计算的图片对，若有，提取图片对重复步骤F到G，若无，继续进行步骤I；

步骤I：选择匹配同名点对最多的一对图片（I，J）解算相机参数和匹配同名点对的世界坐标：

步骤J：在未解图片中找到含有最多已解点的图片L；

步骤K：从已解点的世界坐标和已解点在图片L中的投影图像坐标（x_L，y_L），根据投影方程(5)使用最优化解法计算出图片L的相机参数，再采用已解图片的相机参数根据方程组（6）解算图片L中未解匹配同名点对的世界坐标，得到新的已解点；

步骤L：使用最优化解法微调所有已解图片的相机参数，使得已解点的总体投影误差最小；

步骤M：判断是否含有未解图片，若有，则转回步骤J，如没有，则进行步骤N；

步骤N：标定河工模型稀疏点云：所有图片均已解算出相机参数，所有匹配同名点对均已确定世界坐标，所有已解点组成未标定的河工模型稀疏点云；选取包含完整标定板图像的任意三张图片，将标定板中四个黑白相间的矩形格的9个角点的世界坐标代入未标定的河工模型稀疏点云，将该河工模型稀疏点云中的相机参数和已解点坐标均标定为真实世界坐标，得到已标定河工模型稀疏点云；

步骤O：生成河工模型稠密点云：以已解点为基础，在已解点附近寻找其它匹配同名点并确定世界坐标，最后得到河工模型稠密点云；

步骤P：根据河工模型稠密点云生成河工模型地形：利用河工模型稠密点云的三维坐标数据生成河工模型地形。

2.如权利要求1所述方法，其特征在于，所述步骤F：按顺序提取步骤E所得图片对的序号，对该图片对中的特征点以特征距离进行匹配；具体过程为：设匹配图片对序号为（I，J），定义图片I中特征点p和图片J中的特征点q之间的特征距离d_pq为：

$d_{\mathrm{pq}}=\sqrt{\underset{i=1}{\overset{128}{\mathrm{\Σ}}}{(m_{i,\mathrm{Ip}}-m_{i,\mathrm{Jq}})}^2}---\left(1\right)$

其中m_i,Ip和m_i,Jq分别为图片I中特征点p和图片J中的特征点q的描述符。对图片I中特征点p，遍历图片J中的所有特征点计算（1）式，得到最小和次小的特征距离值，若最小和次小特征距离值之比小于一设定值a（a取值在0.3到0.8之间，取值越小所得匹配点越少，匹配精度越高，相反则所得点越多精度越低），则认为该最小特征距离值对应的特征点和图片I中特征点p匹配，称为匹配同名点对（认为匹配的特征点是河工模型上的同一个点在不同图片中的反映），若大于a，则认为图片J中无特征点与图片I中特征点p匹配。

3.如权利要求1所述方法，其特征在于，所述步骤I：选择匹配同名点对最多的一对图片（I，J）解算相机参数和匹配同名点对的世界坐标：具体包括：一对匹配同名点对中，在图片I中的图像坐标为（x_I，y_I），在图片J中的图像坐标为（x_J，y_J），该匹配同名点对对应的河工模型上同一点的世界坐标为（X，Y，Z），图像坐标与世界坐标之间满足投影关系，以图片I为例：

$\left\{\begin{array}{c}{\left[\begin{array}{ccc}{X}_I& Y_I& Z_I\end{array}\right]}^T=A_I\left[\begin{array}{cc}{R}_I& t_I\end{array}\right]\left[\begin{array}{ccc}X& Y& Z\end{array}\right]\\ \left[\begin{array}{c}{x}_c\\ y_c\end{array}\right]=\left[\begin{array}{c}{X}_I/Z_I\\ Y_I/Z_I\end{array}\right]\\ \left[\begin{array}{c}{x}_I\\ y_I\end{array}\right]=(1+k_1{r}^2+k_2{r}^4)\left[\begin{array}{c}{x}_c\\ y_c\end{array}\right]\end{array}\right.---\left(2\right)$ 其中A_I为相机内参数，R_I为旋转矩阵，t_I为空间平移矢量，（X_I，Y_I，Z_I）为在以图片中心为坐标系原点、图片长短边为x和y轴的坐标系下，世界坐标为（X，Y，Z）的点的坐标；（x_c，y_c）为图像无畸变情况下点（X，Y，Z）在图片中所成像点的图像坐标，k₁、k₂为径向畸变参数，r²=x_c ²+y_c ²，（x_I，y_I）为点（X，Y，Z）在图片中所成像点的实际图像坐标；A_I、R_I、t_I、k₁、k₂均为待标定参数；A_I和R_I的具体形式为：

$A_I=\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]---\left(3\right)$

$\left\{\begin{array}{c}{R}_I=R_x\left({\mathrm{\θ}}_x\right)R_y\left({\mathrm{\θ}}_y\right)R_z\left({\mathrm{\θ}}_z\right)\\ R_x\left({\mathrm{\θ}}_x\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_x& {-\sin \mathrm{\θ}}_x\\ 0& {\sin \mathrm{\θ}}_x& {\cos \mathrm{\θ}}_x\end{array}\right]\\ R_y\left({\mathrm{\θ}}_y\right)=\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_y& 0& {\sin \mathrm{\θ}}_y\\ 0& 1& 0\\ {-\sin \mathrm{\θ}}_y& 0& {\cos \mathrm{\θ}}_y\end{array}\right]\\ R_z\left({\mathrm{\θ}}_z\right)=\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_z& {-\sin \mathrm{\θ}}_z& 0\\ {\sin \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_z& 0\\ 0& 0& 1\end{array}\right]\end{array}\right.---\left(4\right)$

其中f为相机焦距，θ_x、θ_y和θ_z为绕x、y、z轴旋转角度，t_I中包含三个未知数，加上k₁、k₂，投影关系(2)中一共有9个相机参数；将这9个相机参数组成一个9维矢量Θ_I，定义投影算子P

[x_I y_I]＝P(Θ_I,X,Y,Z)                 （5）

算子P的计算过程为公式(2)。因此，对于一个匹配同名点对，可得到方程组

$\left\{\begin{array}{c}\left[\begin{array}{cc}{x}_I& y_I\end{array}\right]=P({\mathrm{\Θ}}_I,X,Y,Z)\\ \left[\begin{array}{cc}{x}_J& y_J\end{array}\right]=P({\mathrm{\Θ}}_J,X,Y,Z)\end{array}\right.---\left(6\right)$

从方程组（6）中根据世界坐标的X、Y、Z分别相等能得到3个独立方程，包含18个未知数。对于图片对（I，J）中的所有匹配同名点对，相机参数Θ_I和Θ_J均是相等的。使用最优化解法，代入匹配同名对的图像坐标，解算出每张图片的相机参数Θ_I和Θ_J，得到两张已解图片的所有相机参数后，根据已知的Θ_I、Θ_J、（x_I，y_I）和（x_J，y_J）以及方程组(6)反算出此图片对中所有匹配同名点对所反映的河工模型上的点的世界坐标（X，Y，Z），这些已算出世界坐标的匹配同名点对称为已解点，其它尚未确定世界坐标的匹配同名点对称为未解点，未确定相机参数的图片称为未解图片。

Images