CN105046704A - 基于圆形目标阵列与观测方程的视觉测量标定方法 - Google Patents

Abstract

本发明公开了一种基于圆形目标阵列与观测方程的视觉测量标定方法,主要包括实验室标定和现场标定,实验室标定中基于拍摄的圆形目标阵列图像求取摄像机的内部参数：水平尺度因子，垂直尺度因子，摄像机拍摄图像中心的横纵坐标，以及径向畸变系数；现场标定中基于修正误差方程对获得的摄像机外部参数及焦距的初始值进行修正，求取精确的焦距和外部参数：摄像机方位角，倾斜角，旋转角和光心坐标。与现有技术相比，本发明中摄像机内部参数只需在实验室标定一次，如果摄像机的位姿发生变化，只需对其外部参数进行重新测量，并带入误差方程进行修正，得到更精确的参数值。本发明适用于多次标定，需标定的参数较少，计算量较小，运行速度较快。

Description

基于圆形目标阵列与观测方程的视觉测量标定方法

技术领域

本发明涉及一种基于圆形目标阵列与观测方程的视觉测量标定方法，属于摄影测量的技术领域。

背景技术

本发明涉及一种基于圆形目标阵列与观测方程的视觉测量标定方法，属于摄影测量的技术领域。

海洋水文测量的对象包括：海水温度、盐度、海流、潮汐、海浪、透明度、水色、海发光、海冰以及海洋—大气相互作用等。近年来，海浪的测量已成为海洋工程领域最为重要的动力荷载之一，对工程结构的安全和造价有极大的影响。海浪的观测方法有多种：如航空测波法、立体摄影法、平面单应法等。而航空法所需的观测成本较高，立体摄影法计算量较大，精度较低。因此，平面单应法成为较好的选择。平面单应法中又涉及单应模型的选取、摄像机标定方法的选取。摄像机标定方法的好坏直接影响到最终观测结果的精度以及观测的效率。目前，面向海洋水文测量的摄像机标定方法已成为新的研究热点。

摄像机标定过程中，由于环境比较恶劣，摄像机不适合架设在船上，而用GPS系统来定位GCP控制点代价又相对较高。即便经过适当的设计，还是很难在如此恶劣的环境中维持稳定。摄像机位置一旦改变需要重新标定，且标定参数较多，计算量较大。

为克服上述限制，本发明设计了一种新的标定方法，用以求解两步法中摄像机的各种外部及内部参数。

发明内容

发明目的：本发明的目的在于解决现有技术中存在的缺陷，提供一种基于圆形目标阵列与观测方程的视觉测量标定方法。

技术方案：本发明提供一种基于圆形目标阵列与观测方程的视觉测量标定方法，包括实验室标定和现场标定两大步骤；

所述实验室标定步骤中基于拍摄的圆形目标阵列图像求取摄像机的内部参数，所述内部参数包括摄像机的水平尺度因子λ_u，垂直尺度因子λ_v，摄像机拍摄图像中心的横坐标u₀和纵坐标v₀，以及径向畸变系数；

所述现场标定步骤中基于修正误差方程对获得的摄像机外部参数及焦距f的初始值进行修正，求取精确的外部参数和焦距；所述外部参数包括摄像机方位角φ，倾斜角τ，旋转角σ和光心坐标(x_c,y_c,z_c)。

所述实验室标定步骤中,具体包括：

(11)用摄像机正直拍摄圆形目标阵列，得到圆形目标阵列n×m的图像；

(12)将圆形目标阵列图像阈值化为二值图像，计算各圆形目标的中心，求出相应的畸变图像坐标；

(13)根据圆形目标阵列的图像坐标和对象空间坐标的关系：和构造包含单应矩阵未知元素的方程组，基于任意四组圆形目标的图像和空间坐标，求得标定系数L1，L2，L4，L5，L6，L8，L9，L10的最小二乘解；其中x，y表示圆形目标阵列中每个目标在空间中的横纵坐标,u，v表示圆形目标阵列中每个目标在图像中的横纵坐标；

(14)根据步骤(13)求得的标定系数，计算圆形目标的无畸变图像的图像坐标估计值，并基于径向畸变模型，求得径向畸变系数；

(15)对畸变图像进行校正，根据校正后的目标点间像素距离与米制距离的比值，求得水平尺度因子和垂直尺度因子。

所述现场标定步骤中，具体包括：

(21)获取GCP的世界坐标(x,y)和图像坐标(u,v)，为每对GCP坐标建立如下函数关系：

$\left\{\begin{array}{c}F(u^{*},\mathrm{\φ},\mathrm{\τ},\mathrm{\σ},f,x_c,y_c,z_c)={\mathrm{qu}}^{*}+of\\ G(v^{*},\mathrm{\φ},\mathrm{\τ},\mathrm{\σ},f,x_c,y_c,z_c)={\mathrm{qv}}^{*}+pf\end{array}\right.---\left(8\right)$

其中，q＝m₃₁Δx+m₃₂Δy+m₃₃Δz，o＝m₁₁Δx+m₁₂Δy+m₁₃Δz，Δx＝x-x_c，Δy＝y-y_c，Δz＝z-z_c，p＝m₂₁Δx+m₂₂Δy+m₂₃Δz，u^*＝(u-u₀)λ_u，v^*＝(v-v₀)λ_v，m_ij为方向余弦，由φ、τ和σ连续旋转推导得到；

(22)将函数F和G按仅剩一阶项的泰勒序列展开得到：

$\left\{\begin{array}{c}0=F_0+\frac{\∂F}{\∂u^{*}}{\mathrm{du}}^{*}+\frac{\∂F}{\∂\mathrm{\φ}}d\mathrm{\φ}+\frac{\∂F}{\∂\mathrm{\τ}}d\mathrm{\τ}+\frac{\∂F}{\∂\mathrm{\σ}}d\mathrm{\σ}+\frac{\∂F}{\∂f}df+\frac{\∂F}{\∂x_c}{\mathrm{dx}}_c+\frac{\∂F}{\∂y_c}{\mathrm{dy}}_c+\frac{\∂F}{\∂z_c}{\mathrm{dz}}_c\\ 0=G_0+\frac{\∂G}{\∂v^{*}}{\mathrm{dv}}^{*}+\frac{\∂G}{\∂\mathrm{\φ}}d\mathrm{\φ}+\frac{\∂G}{\∂\mathrm{\τ}}d\mathrm{\τ}+\frac{\∂G}{\∂\mathrm{\σ}}d\mathrm{\σ}+\frac{\∂G}{\∂f}df+\frac{\∂G}{\∂x_c}{\mathrm{dx}}_c+\frac{\∂G}{\∂y_c}{\mathrm{dy}}_c+\frac{\∂G}{\∂z_c}{\mathrm{dz}}_c\end{array}\right.---\left(9\right)$

其中，F₀和G₀为将实验室标定获得的内部参数及现场标定初始观测获得的外部参数和焦距带入式(9)得到的初始近似值，dφ、dτ、dσ、df、dx_c、dy_c、dz_c分别为φ、τ、σ、f、x_c、y_c、z_c的修正值；设偏导数和均等于q,du^*＝Δu及dv^*＝Δv，得到参数的误差修正方程：

$\left\{\begin{array}{c}\mathrm{\Δ}u-\frac{F_0}{q}=b_{11}d\mathrm{\φ}+b_{12}d\mathrm{\τ}+b_{13}d\mathrm{\σ}+b_{14}df-b_{15}{\mathrm{dx}}_c-b_{16}{\mathrm{dy}}_c-b_{17}{\mathrm{dz}}_c\\ \mathrm{\Δ}v-\frac{G_0}{q}=b_{21}d\mathrm{\φ}+b_{22}d\mathrm{\τ}+b_{23}d\mathrm{\σ}+b_{24}df-b_{25}{\mathrm{dx}}_c-b_{26}{\mathrm{dy}}_c-b_{27}{\mathrm{dz}}_c\end{array}\right.---\left(10\right)$

其中，b_ij系数表示函数F和G以q为尺度相对于未知数的偏导数；

(23)基于式(10)求得所有参数的初始值的修正量，将求得的修正量加入各参数的初始值得到修正后的近似值，将修正后的值带入式(8)并根据计算结果的绝对值进行判断，若小于设定的阈值，则修正结束，求得各未知参数精确的估计值，否则再次计算修正量进行迭代计算。

所述步骤(13)中圆形目标阵列的对象空间坐标的定义为：(x^ij,y^ij)＝(i,j)，其中i和j分别为圆形目标所在的行号和列号，x^ij，y^ij表示圆形目标阵列中每个目标在空间中的相对位置横纵坐标，其值与行列号对应相等。

所述步骤(15)中在圆形目标阵列中分别均匀地取若干对水平和垂直的相邻的圆形目标，将两目标间的米制距离除以像素距离，并将这些比值做加权平均，分别得到横向尺度因子λ_u和纵向尺度因子λ_v。

所述步骤(14)中径向畸变模型为：

Δr＝k₁r³+k₂(7)

其中，是畸变图像上的点(u_d,v_d)到图像中心(u₀,v₀)的距离， $\mathrm{\Δ}r=\sqrt{{(u_d-u_0)}^2+{(v_d-v_0)}^2}-\sqrt{{(u_p-u_0)}^2+{(v_p-v_0)}^2},$ 表示畸变引起的像素位移，(u_p,v_p)为(u_d,v_d)对应的无畸变图像的坐标估计值，k₁和k₂为两个径向畸变系数。

有益效果：与现有技术相比，本发明具有以下优点：

1、内部参数只需在实验室标定一次，如果摄像机的位姿发生变化，只需对其外部参数进行重新测量，并带入参数修正方程进行求精，最终转换成DLT系数，得到精确的观测结果。

2、进行多次标定时，需标定的参数较少，计算量较少，运行速度较快。

综上所述，本发明能够满足面向海洋水文参数测量的要求。

附图说明

为了使本发明的内容更容易被清楚的理解，下面根据的具体实施例并结合附图，对本发明作进一步详细的说明，其中：

图1为本发明的方法总体流程图；

图2为本发明的实验室标定流程图；

图3为本发明的现场标定流程图；

图4为摄像机拍摄圆形目标阵列示意图。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

摄影测量领域中常用的摄像机共线方程：

$\left\{\begin{array}{c}u-u_0=-C_u\[\frac{m_{11}(x-x_c)+m_{12}(y-y_c)+m_{13}(z-z_c)}{m_{31}(x-x_c)+m_{32}(y-y_c)+m_{33}(z-z_c)}\]\\ v-v_0=-C_v\[\frac{m_{21}(x-x_c)+m_{22}(y-y_c)+m_{23}(z-z_c)}{m_{31}(x-x_c)+m_{32}(y-y_c)+m_{33}(z-z_c)}\]\end{array}\right.---\left(1\right)$

该方程中，C_u＝f/λ_u和C_v＝f/λ_v是水平及垂直尺度因子λ_u和λ_v相关于有效焦距的系数。各元素m_ij被称为方向余弦，可由角度φ(方位角)、τ(倾斜角)及σ(旋转角)这三种连续旋转推导得到，(x_c,y_c,z_c)表示摄像机光心坐标，(u,v)表示图像坐标，(x,y,z)表示世界坐标。

$\left\{\begin{array}{c}{m}_{11}=cos\mathrm{\φ}cos\mathrm{\σ}+sin\mathrm{\φ}cos\mathrm{\τ}sin\mathrm{\σ}\\ m_{12}=-sin\mathrm{\φ}\cos \mathrm{\σ}+cos\mathrm{\φ}cos\mathrm{\τ}\sin \mathrm{\σ}\\ m_{13}=sin\mathrm{\τ}\sin \mathrm{\σ}\\ m_{21}=-cos\mathrm{\φ}\sin \mathrm{\σ}+sin\mathrm{\φ}cos\mathrm{\τ}cos\mathrm{\σ}\\ m_{22}=sin\mathrm{\φ}\sin \mathrm{\σ}+\cos \mathrm{\φ}cos\mathrm{\τ}\cos \mathrm{\σ}\\ m_{23}=sin\mathrm{\τ}\cos \mathrm{\σ}\\ m_{31}=\sin \mathrm{\φ}\sin \mathrm{\τ}\\ m_{32}=\cos \mathrm{\φ}\sin \mathrm{\τ}\\ m_{33}=-cos\mathrm{\τ}\end{array}\right.---\left(2\right)$

将式(1)简化为：

$\left\{\begin{array}{c}u=\frac{L_{1}x+L_{2}y+L_{3}z+L_4}{L_{9}x+L_{10}y+L_{11}z+1}\\ v=\frac{L_{5}x+L_{6}y+L_{7}z+L_8}{L_{9}x+L_{10}y+L_{11}z+1}\end{array}\right.---\left(3\right)$

如图1所示，本发明基于圆形目标阵列与观测方程的视觉测量标定方法主要包括实验室标定和现场标定，其中，实验室标定基于拍摄的圆形目标阵列图像求取摄像机的内部参数λ_u、λ_v、u₀、v₀、k₁、k₂。λ_u表示摄像机的水平尺度因子，λ_v表示摄像机的垂直尺度因子，u₀、v₀分别表示摄像机拍摄图像中心的横坐标和纵坐标，k₁、k₂表示径向畸变系数。现场标定基于修正误差方程对获得的摄像机外部参数及焦距f的初始值进行修正，精确求取摄像机的τ、φ、σ、x_c、y_c、z_c6个外部参数及f一个内部参数。其中φ表示摄像机方位角，τ表示摄像机倾斜角，σ表示摄像机旋转角，初始时刻，由摄像机云台上的角偏器读出，x_c、y_c、z_c表示摄像机的光心坐标，并且现场标定前初始值由全站仪粗略测量出，f表示摄像机的焦距，初始时刻由摄像机选用镜头的参数给出。

实验室标定时，对式(3)的共线方程施加一个约束条件：所有控制点都位于一个与摄像机焦平面平行的公共平面，此时平面上所有目标的高度都为零。因此可得到如下式子：

$\left\{\begin{array}{c}{L}_{1}x+L_{2}y+L_4-{\mathrm{uL}}_9-{\mathrm{uL}}_{10}y=u\\ L_{5}x+L_{6}y+L_8-{\mathrm{vL}}_{9}x-{\mathrm{uL}}_{10}y=v\end{array}\right.---\left(4\right)$

在实验室中，选用横纵向分别21个和15个控制点组成的圆形目标阵列，阵列中水平方向和垂直方向的圆形目标之间的距离都为2cm。给出圆形目标阵列的定义空间相对位置坐标的定义：(x^ij,y^ij)＝(i,j)，其中i和j分别为圆形目标所在的行号和列号。x^ij，y^ij表示圆形目标阵列中每个目标在空间中的相对位置横纵坐标，其值与行列号对应相等。

如图2所示，实验室标定步骤具体包括：

步骤(1)摄像机正直拍摄圆形目标阵列，调节摄像机与圆形目标阵列间的距离，使得拍摄的图像正好包含圆形目标阵列的边界点(如图4)。

步骤(2)将拍摄的图像阈值化为二值图像，从低强度的背景中分离出接近白色的目标。计算各目标的中心，求出相应的畸变图像坐标

步骤(3)获取圆形目标阵列n×m的图像和对象空间坐标，按照和构造包含单应矩阵未知元素的方程组：

$\left[\begin{array}{cccccccc}1& 1& 1& -1u_d^{11}& -1u_d^{11}& 0& 0& 0\\ 2& 1& 1& -2u_d^{12}& -2u_d^{12}& 0& 0& 0\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ m& n& 1& -{\mathrm{mu}}_d^{nm}& -{\mathrm{nu}}_d^{nm}& 0& 0& 0\\ 0& 0& 0& -1v_d^{11}& -1v_d^{11}& 1& 1& 1\\ 0& 0& 0& -2v_d^{11}& -2v_d^{11}& 2& 1& 1\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ 0& 0& 0& -{\mathrm{mv}}_d^{nm}& -{\mathrm{nv}}_d^{nm}& m& n& 1\end{array}\right]\left[\begin{array}{c}{L}_1\\ L_2\\ L_3\\ L_4\\ L_5\\ L_6\\ L_7\\ L_8\\ L_9\\ L_{10}\end{array}\right]=\left[\begin{array}{c}{u}_d^{11}\\ u_d^{12}\\ .\\ .\\ .\\ u_d^{nm}\\ v_d^{11}\\ v_d^{12}\\ .\\ .\\ .\\ v_d^{nm}\end{array}\right]---\left(5\right)$

方程中表示圆形目标的畸变图像坐标，任意获取四组圆形目标的图像和空间坐标，代入式(5)，得到其系数L1，L2，L4，L5，L6，L8，L9，L10的最小二乘解。

步骤(4)求得标定系数后，无畸变图像的坐标(u_p,v_p)可估计为：

$\left\{\begin{array}{c}{u}_p^{ij}=\frac{L_{1}j+L_{2}i+L_4}{L_{9}j+L_{10}i+1}\\ v_p^{ij}=\frac{L_{5}j+L_{6}i+L_8}{L_{9}j+L_{10}i+1}\end{array}\right.---\left(6\right)$

为求解畸变系数k₁和k₂，根据受约束的DLT解表示的无畸变模型，我们要计算对应于目标阵列的观测与预测图像坐标间的偏差。

我们选择采用两个系数的奇数阶多项式对径向畸变建模：

Δr＝k₁r³+k₂(7)

其中，是点(u_d,v_d)到图像中心(u₀,v₀)的距离，若获取的图像长宽的像素分别为M、N，则图像中心坐标为(M/2,N/2)； $\mathrm{\Δ}r=\sqrt{{(u_d-u_0)}^2+{(v_d-v_0)}^2}-\sqrt{{(u_p-u_0)}^2+{(v_p-v_0)}^2},$ 表示畸变引起的像素位移。

求得的两个畸变系数k₁和k₂被作为多项式(7)对于观测值r和Δr的最佳拟合解。

步骤(5)根据步骤(4)得到的畸变系数k₁和k₂获取圆形目标畸变校正后的图像坐标，在畸变校正后的圆形目标阵列中分别均匀地取9对水平和垂直的相邻圆形目标，将两目标间的米制距离除以像素距离，并将这些比值做加权平均，得到横向和纵向的尺度因子。

至此，完成了摄像机的实验室标定。

如图3所示，摄像机的现场标定步骤具体包括：

步骤(1)获取GCP(Groundcontrolpoints)的世界坐标(x,y)和图像坐标(u,v)，为每对GCP坐标建立如下两个函数：

$\left\{\begin{array}{c}F(u^{*},\mathrm{\φ},\mathrm{\τ},\mathrm{\σ},f,x_c,y_c,z_c)={\mathrm{qu}}^{*}+of\\ G(v^{*},\mathrm{\φ},\mathrm{\τ},\mathrm{\σ},f,x_c,y_c,z_c)={\mathrm{qv}}^{*}+pf\end{array}\right.---\left(8\right)$

其中，q＝m₃₁Δx+m₃₂Δy+m₃₃Δz，o＝m₁₁Δx+m₁₂Δy+m₁₃Δz，Δx＝x-x_c，Δy＝y-y_c，Δz＝z-z_c，p＝m₂₁Δx+m₂₂Δy+m₂₃Δz，u^*＝(u-u₀)λ_u，v^*＝(v-v₀)λ_v。

步骤(2)将函数F和G按仅剩一阶项的泰勒序列展开得：

$\left\{\begin{array}{c}0=F_0+\frac{\∂F}{\∂u^{*}}{\mathrm{du}}^{*}+\frac{\∂F}{\∂\mathrm{\φ}}d\mathrm{\φ}+\frac{\∂F}{\∂\mathrm{\τ}}d\mathrm{\τ}+\frac{\∂F}{\∂\mathrm{\σ}}d\mathrm{\σ}+\frac{\∂F}{\∂f}df+\frac{\∂F}{\∂x_c}{\mathrm{dx}}_c+\frac{\∂F}{\∂y_c}{\mathrm{dy}}_c+\frac{\∂F}{\∂z_c}{\mathrm{dz}}_c\\ 0=G_0+\frac{\∂G}{\∂v^{*}}{\mathrm{dv}}^{*}+\frac{\∂G}{\∂\mathrm{\φ}}d\mathrm{\φ}+\frac{\∂G}{\∂\mathrm{\τ}}d\mathrm{\τ}+\frac{\∂G}{\∂\mathrm{\σ}}d\mathrm{\σ}+\frac{\∂G}{\∂f}df+\frac{\∂G}{\∂x_c}{\mathrm{dx}}_c+\frac{\∂G}{\∂y_c}{\mathrm{dy}}_c+\frac{\∂G}{\∂z_c}{\mathrm{dz}}_c\end{array}\right.---\left(9\right)$

其中F₀和G₀为初始近似值，初始时刻，将实验室标定得到的6个内参及现场标定测量得到的7个参数带入式子(9)得到的值，dφ、dτ、dσ、df、dx_c、dy_c、dz_c分别为φ、τ、σ、f、x_c、y_c、z_c的修正值。

假设偏导数和均等于q,du^*＝Δu及dv^*＝Δv，则参数的误差修正方程如下：

$\left\{\begin{array}{c}\mathrm{\Δ}u-\frac{F_0}{q}=b_{11}d\mathrm{\φ}+b_{12}d\mathrm{\τ}+b_{13}d\mathrm{\σ}+b_{14}df-b_{15}{\mathrm{dx}}_c-b_{16}{\mathrm{dy}}_c-b_{17}{\mathrm{dz}}_c\\ \mathrm{\Δ}v-\frac{G_0}{q}=b_{21}d\mathrm{\φ}+b_{22}d\mathrm{\τ}+b_{23}d\mathrm{\σ}+b_{24}df-b_{25}{\mathrm{dx}}_c-b_{26}{\mathrm{dy}}_c-b_{27}{\mathrm{dz}}_c\end{array}\right.---\left(10\right)$

其中b_ij系数表示函数F和G(以q为尺度)相对于未知数的偏导数。

步骤(3)由于所有参数对函数F和G的初始近似是基于观测值的，初始值的修正通过式(10)的最小二乘解获得。在初次求解后，将求得的修正量加入初始近似值来获得修正后的近似值。重复该过程(迭代)直到修正量的大小可忽略，当把修正过后的7个参数τ、φ、σ、x_c、y_c、z_c、f,带入式子(8)，计算得到的结果的绝对值小于设定的阈值(如0.0001)时，修正量就可忽略，于是便求得未知数更精确的估计值。

上述现场标定过程中由于每组已知的控制点对可得到两个方程，所有七个余下的摄像机标定参数可通过最少四个可见的GCP来求解。通过采用预先求解的值来约束摄像机坐标(x_c,y_c,z_c)，未知数的数量可降至四个，且方程组可唯一地由两个GCP求解。更多的控制点将产生具有冗余信息的超定方程组，可在最小二乘的意义下求解。可用各种其它已知和未知现场参数的组合来减少所需GCP的数量或增加冗余度。

例如，通过水平线的测量(由此可求得倾斜和转动角)及摄像机位置的知识，可仅用一个勘测的GCP，采用模型解算得到焦距和方位角。

当然，不管输入个数是多少，任意现场参数都可简单地用预先求解的值来约束。

一旦求得了所有必须的共线参数值，我们通常将这些估计值转换为等价的

DLT(DirectLinearTransformation)标定系数，转换关系如下：

$\left\{\begin{array}{c}L=-(x_c{m}_{31}+y_c{m}_{32}+z_c{m}_{33})\\ L_1=(u_0{m}_{31}+f{m}_{11})/\left({\mathrm{\λ}}_{u}L\right)\\ L_2=(u_0{m}_{32}+f{m}_{12})/\left({\mathrm{\λ}}_{u}L\right)\\ L_3=(u_0{m}_{33}+f{m}_{13})/\left({\mathrm{\λ}}_{u}L\right)\\ L_4=-(L_1{x}_c+L_2{y}_c+L_3{z}_c)\\ L_5=m_{31}/L\\ L_6=m_{32}/L\\ L_7=m_{33}/L\\ L_8=(v_0{m}_{31}+f{m}_{21})/\left({\mathrm{\λ}}_{v}L\right)\\ L_9=(v_0{m}_{32}+f{m}_{22})/\left({\mathrm{\λ}}_{v}L\right)\\ L_{10}=(v_0{m}_{33}+f{m}_{23})/\left({\mathrm{\λ}}_{v}L\right)\\ L_{11}=-(L_8{x}_c+L_9{y}_c+L_{10}{z}_c)\end{array}\right.$

通过上述实施例可以看出，本发明的中摄像机只需在实验室标定1次，获得的内部参数可以重复使用，进行2次标定时，求取参数少，计算量小计算过程简单。

Claims (8)

1.基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:包括实验室标定和现场标定两大步骤；

所述实验室标定步骤中基于拍摄的圆形目标阵列图像求取摄像机的内部参数，所述内部参数包括摄像机的水平尺度因子λ_u，垂直尺度因子λ_v，摄像机拍摄图像中心的横坐标u₀和纵坐标v₀，以及径向畸变系数；

所述现场标定步骤中基于修正误差方程对获得的摄像机外部参数及焦距f的初始值进行修正，求取精确的外部参数和焦距；所述外部参数包括摄像机方位角φ，倾斜角τ，旋转角σ和光心坐标(x_c,y_c,z_c)。

2.根据权利要求1所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:所述实验室标定步骤中,具体包括：

(11)用摄像机正直拍摄圆形目标阵列，得到圆形目标阵列的图像；

(12)将圆形目标阵列图像阈值化为二值图像，计算各圆形目标的中心，求出相应的畸变图像坐标；

(13)根据圆形目标阵列的图像坐标和对象空间坐标的关系：和构造包含单应矩阵未知元素的方程组，基于任意四组圆形目标的图像和空间坐标，求得标定系数L1，L2，L4，L5，L6，L8，L9，L10的最小二乘解；其中x，y表示圆形目标阵列中每个目标在空间中的横纵坐标,u，v表示圆形目标阵列中每个目标在图像中的横纵坐标；

(14)根据步骤(13)求得的标定系数，计算圆形目标的无畸变图像的图像坐标估计值，并基于径向畸变模型，求得径向畸变系数；

(15)对畸变图像进行校正，根据校正后的目标点间像素距离与米制距离的比值，求得水平尺度因子和垂直尺度因子。

3.根据权利要求1所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:所述现场标定步骤中，具体包括：

(21)获取GCP的世界坐标(x,y)和图像坐标(u,v)，为每对GCP坐标建立如下函数关系：

$\left\{\begin{array}{c}F(u^{*},\mathrm{\φ},\mathrm{\τ},\mathrm{\σ},f,x_c,y_c,z_c)={\mathrm{qu}}^{*}+of\\ G(v^{*},\mathrm{\φ},\mathrm{\τ},\mathrm{\σ},f,x_c,y_c,z_c)={\mathrm{qv}}^{*}+of\end{array}\right.---\left(8\right)$

其中，q＝m₃₁Δx+m₃₂Δy+m₃₃Δz，o＝m₁₁Δx+m₁₂Δy+m₁₃Δz，Δx＝x-x_c，Δy＝y-y_c，Δz＝z-z_c，p＝m₂₁Δx+m₂₂Δy+m₂₃Δz，u^*＝(u-u₀)λ_u，v^*＝(v-v₀)λ_v，m_ij为方向余弦，由φ、τ和σ连续旋转推导得到；

(22)将函数F和G按仅剩一阶项的泰勒序列展开得到：

$\left\{\begin{array}{c}0=F_0+\frac{\∂F}{\∂u^{*}}{\mathrm{du}}^{*}+\frac{\∂F}{\∂\mathrm{\φ}}d\mathrm{\φ}+\frac{\∂F}{\∂\mathrm{\τ}}d\mathrm{\τ}+\frac{\∂F}{\∂\mathrm{\σ}}d\mathrm{\σ}+\frac{\∂F}{\∂f}df+\frac{\∂F}{\∂x_c}{\mathrm{dx}}_c+\frac{\∂F}{\∂y_c}{\mathrm{dy}}_c+\frac{\∂F}{\∂z_c}{\mathrm{dz}}_c\\ 0=G_0+\frac{\∂G}{\∂v^{*}}{\mathrm{dv}}^{*}+\frac{\∂G}{\∂\mathrm{\φ}}d\mathrm{\φ}+\frac{\∂G}{\∂\mathrm{\τ}}d\mathrm{\τ}+\frac{\∂G}{\∂\mathrm{\σ}}d\mathrm{\σ}+\frac{\∂G}{\∂f}df+\frac{\∂G}{\∂x_c}{\mathrm{dx}}_c+\frac{\∂G}{\∂y_c}{\mathrm{dy}}_c+\frac{\∂G}{\∂z_c}{\mathrm{dz}}_c\end{array}\right.---\left(9\right)$

其中，F₀和G₀为将实验室标定获得的内部参数及现场标定初始观测获得的外部参数和焦距带入式(9)得到的初始近似值，dφ、dτ、dσ、df、dx_c、dy_c、dz_c分别为φ、τ、σ、f、x_c、y_c、z_c的修正值；设偏导数和均等于q,du^*＝Δu及dv^*＝Δv，得到参数的误差修正方程：

$\left\{\begin{array}{c}\mathrm{\Δ}u-\frac{F_0}{q}=b_{11}d\mathrm{\φ}+b_{12}d\mathrm{\τ}+b_{13}d\mathrm{\σ}+b_{14}df-b_{15}{\mathrm{dx}}_c-b_{16}{\mathrm{dy}}_c-b_{17}{\mathrm{dz}}_c\\ \mathrm{\Δ}v-\frac{G_0}{q}=b_{21}d\mathrm{\φ}+b_{22}d\mathrm{\τ}+b_{23}d\mathrm{\σ}+b_{24}df-b_{25}{\mathrm{dx}}_c-b_{26}{\mathrm{dy}}_c-b_{27}{\mathrm{dz}}_c\end{array}\right.---\left(10\right)$

其中，b_ij系数表示函数F和G以q为尺度相对于未知数的偏导数；

(23)基于式(10)求得所有参数的初始值的修正量，将求得的修正量加入各参数的初始值得到修正后的近似值，将修正后的值带入式(8)并根据计算结果的绝对值进行判断，若小于设定的阈值，则修正结束，求得各未知参数精确的估计值，否则再次计算修正量进行迭代计算。

4.根据权利要求2所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:所述步骤(13)中圆形目标阵列的对象空间坐标的定义为：(x^ij,y^ij)＝(i,j)，其中i和j分别为圆形目标所在的行号和列号，x^ij，y^ij表示圆形目标阵列中每个目标在空间中的相对位置横纵坐标，其值与行列号对应相等。

5.根据权利要求2所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:所述步骤(15)中在圆形目标阵列中分别均匀地取若干对水平和垂直的相邻的圆形目标，将两目标间的米制距离除以像素距离，并将这些比值做加权平均，分别得到横向尺度因子λ_u和纵向尺度因子λ_v。

6.根据权利要求2所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:所述步骤(14)中径向畸变模型为：

Δr＝k₁r³+k₂(7)

其中，是畸变图像上的点(u_d,v_d)到图像中心(u₀,v₀)的距离， $\mathrm{\Δ}r=\sqrt{{(u_d-u_0)}^2+{(v_d-v_0)}^2}-\sqrt{{(u_p-u_0)}^2+{(v_p-v_0)}^2},$ 表示畸变引起的像素位移，(u_p,v_p)为(u_d,v_d)对应的无畸变图像的坐标估计值，k₁和k₂为两个径向畸变系数。

7.根据权利要求3所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于:通过最少四个可见的GCP来求解摄像机焦距f，方位角φ，倾斜角τ，旋转角σ和光心坐标(x_c,y_c,z_c)。

8.根据权利要求3所述的基于圆形目标阵列与观测方程的视觉测量标定方法，其特征在于：采用预先求解的值来约束摄像机光心坐标(x_c,y_c,z_c)，通过最少两个可见的GCP来求解摄像机焦距f，方位角φ，倾斜角τ和旋转角σ。