CN102855620B - 基于球形投影模型的纯旋转摄像机自标定方法 - Google Patents

Abstract

一种基于球形投影模型的纯旋转摄像机自标定方法。针对针孔摄像机内参数的标定任务，本发明提出了一种新颖的基于球形投影模型的纯旋转自标定方法。首先构造了针孔摄像机的球形投影模型,之后分析了空间静止点对应的球面投影点之间的距离在摄像机纯旋转时不变；然后根据该性质构造了内参数的约束方程组；进而以非线性最小二乘算法求解该方程组。相比现有方法，本发明利用两幅图像上对应的点特征即可得到内参数，因而无需复杂的矩阵数值运算，并且仅需两幅图像上的4个匹配点即可完成对摄像机4个内参数的标定，且均适用于在线与离线标定。仿真与实验结果表明，本发明简单实用并且标定精度高，且对图像噪声与平移噪声具有很好的鲁棒性。

Description

基于球形投影模型的纯旋转摄像机自标定方法

技术领域

本发明属于计算机视觉与摄像机标定的技术领域，特别是涉及一种基于球形投影模型的纯旋转摄像机自标定方法。

背景技术

摄像机内参数的标定是指将其获得的过程，它是计算机视觉方面的经典问题之一，同时也是基于视觉反馈的控制系统的一项重要技术。摄像机的标定方法大致可分为传统标定法与自标定法。传统的标定法以Tsai的三维标定块标定法[1]以及张正友的平面标定法[2]为代表，传统方法标定精度高，然而均为离线方法且需要精确度高的标定块等装置。自标定方法由Faugeras提出[3]，是指不需标定块等装置，仅通过图像点间的对应关系来进行标定的过程[4-5]。Faugeras在文献[3]中提出了利用Kruppa约束方程对内参数进行求解，考虑到Kruppa方程难以求解且不稳定情况，诸如Pollefeys模约束[6]等分层逐步方法被提出，但这些方法仍存在计算复杂等问题。

对于做特殊运动的摄像机，自标定算法的复杂度会降低且往往能获得线性解[4]。对于其特殊运动的研究，集中在了纯旋转等方面[7-10]。Hartley在文献[7]中提出了经典的基于纯旋转的自标定法，然而该方法需要至少3幅图像间的点对应；Wang等针对基于单应矩阵(Homography)的纯旋转自标定方法，从理论上分析了微小平移所造成的标定误差[11]；Zhang与Wong利用转台图像序列对应的消影点不变等性质，得到了3参数的自标定结果[12]。然而以上的纯旋转自标定方法均以绝对二乘曲线与对极几何为基础，因而不可避免地需要利用复杂的矩阵数值计算。最近，方勇纯等人从控制理论角度出发，设计了一种基于非线性观测器的纯旋转自标定方法，该方法无需矩阵数值计算，并且得到了内参数的全局指数收敛性能[13]，然而该方法需要多幅图像才能使内参数收敛。因而如何设计出一种避免矩阵数值计算且需要较少图像即可得到内参数的方法，是亟待解决的问题。

近年来，全景摄像机(Omnidirectional Cameras)的球形投影模型受到了研究者的关注。首先，Geyer等人将各式的全景摄像机以及平面投影摄像机(即针孔摄像机)以统一的球形投影模型进行了表达[14]，因而在很大程度上方便了对各式摄像机的分析。进而Mariottini等人利用该球形投影模型的自极点(auto-epipolar)性质以全景相机实现了移动机器人的视觉镇定控制[15]；Becerra等人由球形投影模型的1维三焦点张量(1D trifocal tensor)，并结合滑模控制实现了移动机器人的镇定任务[16]；Fomena等人利用该模型以针孔摄像机完成了操作臂的视觉伺服(Visual Servoing)任务[17]，并且得到了相比于经典IBVS方法(Image-Based Visual Servoing)更好的性能。本发明根据该全景摄像机球形投影模型受到启发，完成了对针孔摄像机的内参数自标定任务。

发明内容

本发明的目的是解决现有自标定技术存在的上述不足，提供一种基于球形投影模型的纯旋转摄像机自标定方法。

本发明提出了一种新颖的基于球形投影模型的纯旋转摄像机自标定方法。该方法最大的特点是直接利用两幅图像上对应的点特征即可得到摄像机的内参数。因而避免了现有自标定方法需要复杂的矩阵数值运算的问题，并且仅需两幅图像上的4个匹配点即可完成对4内参数摄像机的标定，而且均适用于在线标定与离线标定。具体而言，首先描述了本文定义的针孔摄像机的球形投影模型(Spherical Projection Model)。然后分析了对于纯旋转摄像机，空间静止点对应的球面投影点之间的距离不变；之后，根据该性质构造了关于内参数的约束方程组；进而利用非线性最小二乘算法对该方程组进行求解。仿真与实验结果表明，本文方法不仅简单实用，并且标定精度高，且对图像噪声与平移噪声具有较好的鲁棒性，因而具有很好的实际应用意义。

本发明提供的基于球形投影模型的纯旋转摄像机自标定方法包括：

第1，构造针孔摄像机的球形投影模型

定义P_i，P_j分别表示第i，j个空间点。像素坐标系的横坐标轴与纵坐标轴分别以u，v表示。以表示摄像机坐标系，其中的原点在摄像机光心位置，的z轴与摄像机光轴重合，x轴方向与u轴方向相同，y轴方向与v轴方向相同。f为摄像机焦距，f的单位为米；cp_i，cp_j表示点P_i，P_j对应的图像像素点在下的位置。表示以的原点为球心的单位虚拟球面；s_i，s_j分别为cp_i，cp_j对应在上的投影点，称其为球面投影点；

对于作纯旋转运动的摄像机，与分别表示摄像机在参考位姿处与经过纯旋转运动后的坐标系；s_i，s_j与s_i′，s_j′分别表示点P_i，P_j在与下的球面投影点；关于s_i,s_j与s_i′，s_j′有定理1所述性质：

定理1:球面投影点之间向量的模长在摄像机作纯旋转运动时不变，如式(1)所示：

||s_i-s_j||₂＝||s_i′-s_j′||₂    (1)

定理1呈现了球面投影点在纯旋转摄像机下的性质，根据摄像机的针孔投影模型，其内参数包括f_x，f_y，u₀，v₀；其中f_x，f_y分别为焦距对应于u，v方向的像素块个数，即f_x＝f/d_x，f_y＝f/d_y；其中d_x，d_y分别为单个像素块在u，v方向的长度，单位为米；(u₀，v₀)为图像主点坐标；因此，本发明的目的为根据空间特征点并利用定理1，对纯旋转运动下的摄像机作内参数自标定；

第2，纯旋转运动下的摄像机内参数自标定

第2.1，构造约束方程组

首先，推导并构造了含有摄像机内参数的约束方程为：

$\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}=\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}}---\left(9\right)$

其中

$\left\{\begin{array}{c}{a}_{\mathrm{ij}}=(u_i-u_0)(u_j-u_0)+(v_i-v_0)(v_j-v_0){\mathrm{\γ}}^2\\ b_{\mathrm{ij}}={(u_i-u_0)}^2+{(v_i-v_0)}^2{\mathrm{\γ}}^2;\\ c_{\mathrm{ij}}={(u_j-u_0)}^2+{(v_j+v_0)}^2{\mathrm{\γ}}^2\\ l_{\mathrm{ij}}=({u_i}^{\′}-u_0)({u_j}^{\′}-u_0)+({v_i}^{\′}-v_0)({v_j}^{\′}-v_0){\mathrm{\γ}}^2\\ m_{\mathrm{ij}}={({u_i}^{\′}-u_0)}^2+{({v_i}^{\′}-v_0)}^2{\mathrm{\γ}}^2;\\ {{n_{\mathrm{ij}}=({u_j}^{\′}-u_0)}^2+{({v_j}^{\′}-v_0)}^2\mathrm{\γ}}^2\end{array}\right.$

其中定义了像素块在v方向与u方向的长度比(u_i，v_i)，(u_j，v_j)分别为点P_i与P_j对应的图像像素坐标，(u_i′，v_i′)，(u_j′，v_j′)分别为点P_i与P_j经过摄像机纯旋转之后对应的图像像素坐标；

然后可以利用4个空间点分别得到6个公式(9)形式的约束方程，构成约束方程组，其中设定任意两空间点与摄像机光心不共线；

第2.2，利用非线性最小二乘算法对约束方程组进行求解

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少4个空间点最小化如下的目标函数J_whole(·)，得到u₀，v₀，f_x，γ的解：

$J_{\mathrm{whole}}(u_0,v_0,f_x,\mathrm{\γ})=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{(\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}-\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}})}^2---\left(20\right)$

最后利用f_y＝f_x/γ得到f_y，进而完成纯旋转运动下的摄像机内参数自标定。

本发明方法的理论依据及推导过程

第1，定义了针孔摄像机的球形投影模型

附图1给出了本发明定义的针孔摄像机球形投影模型，其中点P_i，P_j分别表示第i，j个空间点，其中所需空间点的个数大于4即可；像素坐标系的横坐标轴与纵坐标轴分别以u，v表示。以表示摄像机坐标系，其中的原点在摄像机光心位置，的z轴与摄像机光轴重合，x轴方向与u轴方向相同，y轴方向与v轴方向相同。f为摄像机焦距(单位为米)。cp_i，cp_j为点P_i，P_j对应的图像像素点在下的坐标。表示以的原点为球心的单位虚拟球面。s_i，s_j分别为cp_i，cp_j对应在上的投影点，本发明称其为球面投影点。

对于作纯旋转运动的摄像机(以绕y轴旋转为例)，如附图2所示。其中(u₀，v₀)为图像主点坐标，与分别表示摄像机在参考位姿处与经过纯旋转运动后的坐标系。s_i，s_j与s_i′，s_j′分别表示点P_i，P_j在与下的球面投影点。关于s_i，s_j与s_i′，s_j′有定理1所述性质。

定理1：球面投影点之间向量的模长在摄像机作纯旋转运动时不变，如式(1)所示(定理证明见第3小节的附录A)：

||s_i-s_j||₂＝||s_i′-s_j′||₂    (1)

定理1呈现了球面投影点在纯旋转摄像机下的性质。根据摄像机的针孔投影模型，其内参数矩阵为：

$A=\left[\begin{array}{ccc}{f}_x& 0& u_0\\ 0& f_y& v_0\\ 0& 0& 1\end{array}\right]---\left(2\right)$

其中f_x，f_y分别为焦距对应于u，v方向的像素块个数：

f_x＝f/d_x，f_y＝f/d_y    (3)

其中d_x，d_y分别为单个像素块在u，v方向的长度，单位为米。因此，本发明的目的为根据空间特征点并利用定理1，对纯旋转运动下的摄像机作内参数自标定。

第2，纯旋转运动下的摄像机内参数自标定

本节根据球面投影点在摄像机纯旋转运动下的性质，得到关于内参数的约束方程组；然而鉴于该方程组难以得到解析解，本发明将采用非线性最小二乘方法进行数值最优化求解。

第2.1，构造约束方程组

根据摄像机的投影模型，点P_i，P_j对应的图像像素点在下的坐标为：

${\mathrm{cp}}_i=\left[\begin{array}{c}(u_i-u_0)d_x\\ (v_i-v_0)d_y\\ f\end{array}\right],$ ${\mathrm{cp}}_j=\left[\begin{array}{c}(u_j-u_0)d_x\\ (v_j-v_0)d_y\\ f\end{array}\right]---\left(4\right)$

其中(u_i，v_i)，(u_j，v_j)分别为P_i与P_j对应的图像像素坐标，将cp_i，cp_j投影到上，可得

$s_i=\frac{{\mathrm{cp}}_i}{{\left|\right|{\mathrm{cp}}_i\left|\right|}_2},$ $s_j=\frac{{\mathrm{cp}}_j}{{\left|\right|{\mathrm{cp}}_j\left|\right|}_2}---\left(5\right)$

其中||cp_i||₂根据下式计算得到：

${\left|\right|{\mathrm{cp}}_i\left|\right|}_2=\sqrt{{(u_i-u_0)}^2{d_x}^2+{(v_i-v_0)}^2{d_y}^2+{\left(f\right)}^2}---\left(6\right)$

同样地，点P_i，P_i对应的图像像素点在下的坐标为：

${{\mathrm{cp}}_i}^{\′}=\left[\begin{array}{c}({u_i}^{\′}-u_0)d_x\\ ({v_i}^{\′}-v_0)d_y\\ f\end{array}\right],$ ${{\mathrm{cp}}_j}^{\′}=\left[\begin{array}{c}({u_j}^{\′}-u_0)d_x\\ ({v_j}^{\′}-v_0)d_y\\ f\end{array}\right]---\left(7\right)$

将cp_i′，cp_j′投影到上，可得：

${s_i}^{\′}=\frac{{{\mathrm{cp}}_i}^{\′}}{\left|\right|{{\mathrm{cp}}_i}^{\′}\left|\right|},$ ${s_j}^{\′}=\frac{{{\mathrm{cp}}_j}^{\′}}{\left|\right|{{\mathrm{cp}}_j}^{\′}\left|\right|}---\left(8\right)$

根据定理1中描述的||s_i-s_j||₂＝||s_i′-s_j′||₂性质，将式(4)-(8)代入并整理，可得：

$\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}=\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}}---\left(9\right)$

其中

$\left\{\begin{array}{c}{a}_{\mathrm{ij}}=(u_i-u_0)(u_j-u_0)+(v_i-v_0)(v_j-v_0){\mathrm{\γ}}^2\\ b_{\mathrm{ij}}={(u_i-u_0)}^2+{(v_i-v_0)}^2{\mathrm{\γ}}^2\\ c_{\mathrm{ij}}={(u_j-u_0)}^2+{(v_j-v_0)}^2{\mathrm{\γ}}^2\end{array}\right.---\left(10\right)$

且：

$\left\{\begin{array}{c}{l}_{\mathrm{ij}}=({u_i}^{\′}-u_0)({u_j}^{\′}-u_0)+({v_i}^{\′}-v_0)({v_j}^{\′}-v_0){\mathrm{\γ}}^2\\ m_{\mathrm{ij}}={({u_i}^{\′}-u_0)}^2+{({v_i}^{\′}-v_0)}^2{\mathrm{\γ}}^2\\ n_{\mathrm{ij}}={({u_j}^{\′}-u_0)}^2+{({v_j}^{\′}-v_0)}^2{\mathrm{\γ}}^2\end{array}\right.---\left(10\right)$

其中定义了像素块在v方向与u方向的长度比：

$\mathrm{\γ}\stackrel{\mathrm{\Δ}}{=}\frac{d_y}{d_x}=\frac{f_x}{f_y}---\left(12\right)$

由以上分析可知，需要至少4个式(9)形式的等式以构成对内参数的约束方程组才能得到u₀，v₀，f_x，f_y的唯一解。由于该方程组为高次多元非线性方程组，因而难以得到解析解，因此本文采用Levenberg-Marquardt(LM)^[18]非线性最小二乘方法进行数值优化求解。

注1：利用4个空间点可以得到6个式(9)形式方程，构成约束方程组，当其存在4个线性无关方程时，即可求解出4个内参数。但是当存在两个空间点与摄像机光心共线时，至多存在3个线性无关方程，即出现了退化情况。因此为避免方程组的退化，本文假设任意两空间点与摄像机光心不共线。易知该假设是合理的，并且其它基于纯旋转的自标定方法亦需要该假设。

第2.2，利用非线性最小二乘算法对约束方程组进行求解

在f_x＝f_y情况下，可知γ＝1，以此求出u₀，v₀，f_x作为LM优化算法的初值。首先，对于u₀，v₀的初值u_0init,v_0init采用图像像素的中心坐标：

$u_{0\mathrm{init}}=\frac{u_{\max }}{2},$ $v_{0\mathrm{init}}=\frac{v_{\max }}{2}---\left(13\right)$

易知仅利用1对对应点并根据式(9)即可得到f_x的初值f_xinit，，其计算方法为：利用空间点P₁，P₂，并将γ＝1代入(9)式，经整理可得：

A₆f_x ⁶+A₄f_x ⁴+A₂f_x ²+A₀＝0    (14)

其中

A₆＝2a₁₂-2l₁₂+(m₁₂+n₁₂)-(b₁₂+c₁₂)    (15)

A₄＝m₁₂n₁₂-b₁₂c₁₂+2a₁₂(m₁₂+n₁₂)

                                     (16)

    -2l₁₂(b₁₂+c₁₂)+a₁₂ ²-l₁₂ ²

A₂＝a₁₂ ²(m₁₂+n₁₂)-l₁₂ ²(b₁₂+c₁₂)

                                     (17)

    +2a₁₂m₁₂n₁₂-2l₁₂b₁₂c₁₂

A₀＝a₁₂ ²m₁₂n₁₂-l₁₂ ²b₁₂c₁₂            (18)

高次方程(14)有绝对值相等的1个正实根与1个负实根，4个实部为0的复数根，因此取该正实根作为解即可。

接下来分两步进行优化：首先在γ≡1的情况下，利用LM方法与至少3个点对u₀，v₀，f_x进行优化求解，即该LM方法在假设f_x＝f_y的前提下进行。方法是最小化如下的目标函数J_local(·)：

$J_{\mathrm{local}}(u_0,v_0,f_x,\mathrm{\γ}\≡1)=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{(\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}-\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}})}^2---\left(19\right)$

因而得到用于第二步骤LM优化的初值u_0init′，v_0init′，f_xinit′。第二步骤为利用LM算法与至少4个点最小化如下的目标函数J_whole(·)，因此得到了u₀，v₀，f_x，γ的解。

$J_{\mathrm{whole}}(u_0,v_0,f_x,\mathrm{\γ})=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{(\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}-\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}})}^2---\left(20\right)$

最后利用(12)式中f_y＝f_x/γ得到f_y，进而完成纯旋转运动下的摄像机内参数自标定。

注2：文献[4]指出，对于纯旋转摄像机，其自标定方法在如下情况下存在多义性：(1)若摄像机绕y轴旋转，则仅能标定出f_x，无法标定出f_y；(2)若绕x轴旋转，仅能标定出f_y，无法标定出f_x。(3)若绕z轴旋转，则仅能标定出γ。针对这些特殊旋转情况，第3小节的附录B给出了相应的球形投影下的自标定方法。

注3：本文分析的是4内参数摄像机模型，易知3参数的摄像机(即f_x＝f_y)是本文模型的特例，因而最小化目标函数J_local(·)即可得到标定结果，且除纯绕z轴旋转外，无需考虑注2中提到的限制条件。

第3，附录

第3.1，附录A，定理1证明

本发明在此给出定理1的证明。

证明：向量s_i-s_j的模长d_ij由下式计算：

$d_{\mathrm{ij}}\stackrel{\mathrm{\Δ}}{=}{\left|\right|s_i-s_j\left|\right|}_2=\sqrt{2-{{2s}_i}^T{s}_j}---\left(21\right)$

对其关于时间求导得：

${\stackrel{\·}{d}}_{\mathrm{ij}}=-\frac{1}{d_{\mathrm{ij}}}({s_j}^T{\stackrel{\·}{s}}_i+{s_i}^T{\stackrel{\·}{s}}_j)---\left(22\right)$

对于空间点P_i，其关于的运动学关系为^[19]：

${\stackrel{\·}{P}}_i=-v-w\×P_i---\left(23\right)$

其中v，w分别表示摄像机在其自身坐标系下的角速度与线速度。由于P_i＝s_i||P_i||，且||P_i||的变化率为v在P_i方向上的投影，即：

$\frac{d\left|\right|P_i\left|\right|}{\mathrm{dt}}=-v^T{s}_i---\left(24\right)$

因此由上两式并加以整理可得：

${\stackrel{\·}{s}}_i=-\frac{I_3-s_i{s_i}^T}{\left|\right|P_i\left|\right|}v+{\left[s_i\right]}_{\×}w---\left(25\right)$

同样地可以得到将它们带入并化简可得：

${\stackrel{\·}{d}}_{\mathrm{ij}}=\frac{1}{d_{\mathrm{ij}}}(\frac{{s_j}^T-{s_j}^T{s}_i{s_i}^T}{\left|\right|P_i\left|\right|}+\frac{{s_i}^T-{s_i}^T{s}_j{s_j}^T}{\left|\right|P_j\left|\right|})v---\left(26\right)$

因此可知仅与v相关，与w无关。即当摄像机作纯旋转运动时，d_ij不变。因此得到

||s_i-s_j||₂＝||s_i′-s_j′||₂    (27)

                                                        ■

第3.2，附录B，特殊纯旋转下的自标定

这种特殊旋转情况是指摄像机绕单一坐标轴旋转。由注2可知，在该情况下，无法将全部4个内参数标定出，因而需要摄像机绕其它轴做第二次旋转。但这种特殊旋转的优势在于它会给出更强的约束，从而使算法更加简单。

B.1：绕y轴旋转

对于绕y轴旋转情况，至少需要再绕x轴旋转才能标定出4个内参数。下面给出该种情况的自标定方法(绕x轴旋转的方法类似)：

当摄像机绕y轴旋转时，球面投影点的y坐标不变，即：

s_iy＝s_iy′    (28)

因此针对点P_i，根据式(5)(8)，整理并化简可得

(v_i-v₀)²(u_i′-u₀)²+(v_i′-v₀)²(u_i-u₀)²＝[(v_i′-v₀)²-(v_i-v₀)²]f_x ²    (29)

因此根据3个空间点，利用LM算法求解如下目标函数J_y axis(·)的最小值可解出u₀，v₀，f_x。

$J_{\mathrm{yaxis}}(u_0,v_0,f_x)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({(v_i-v_0)}^2{({u_i}^{\′}-u_0)}^2+{({v_i}^{\′}-v_0)}^2{(u_i-u_0)}^2----\left(30\right)$

$-{({({v_i}^{\′}-v_0)}^2-{(v_i-v_0)}^2{f_x}^2)}^2$

再将摄像机在参考坐标系下绕x轴旋转，则球面投影点的x坐标不变，即

s_ix＝s_ix′    (31)

因此根据(5)(8)有：

(u_i-u₀)²(v_i′-v₀)²+(u_i′-u₀)²(v_i-v₀)²＝((u_i′-u₀)²-(u_i-u₀)²)f_y ²    (32)

因此可以将f_y以解析的方式求出。

B.2：绕z轴旋转

对于绕z轴旋转情况，球面投影点的z坐标不变：

s_iz＝s_iz′    (33)

因此针对点P_i，根据式(5)(8)，整理后可得：

(u_i-u₀)²-(u_i′-u₀)²＝((v_i′-v₀)²-(v_i-v₀)²)γ²    (34)

因此根据3个空间点，利用LM算法求解如下目标函数J_z axis(·)的最小值可解出u₀，v₀，γ。

$J_{\mathrm{zaxis}}(u_0,v_0,\mathrm{\γ})=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({(u_i-u_0)}^2-{({u_i}^{\′}-u_0)}^2-({({v_i}^{\′}-v_0)}^2-{(v_i-v_0)}^2){\mathrm{\γ}}^2)}^2---\left(35\right)$

之后在参考坐标系下绕x轴(或y轴)，以(29)式(或(32)式)求出f_x(或f_y)以完成对4个内参数的求取。

本发明的优点和有益效果

本发明提出了一种基于球形投影模型的纯旋转摄像机自标定方法。本发明直接利用两幅图像上对应的点特征即可得到摄像机的内参数。因而无需复杂的矩阵数值运算，并且仅需两幅图像上的4个匹配点即可完成对摄像机4个内参数的标定，且均适用于在线标定与离线标定。仿真与实验结果表明，本发明简单实用并且标定精度高，且对图像噪声与平移噪声具有很好的鲁棒性。

附图说明：

图1为针孔摄像机的球形投影模型；

图2为摄像机纯旋转运动对应的球面投影点；其中(a)图表示摄像机在纯旋转运动之前的自标定场景，(b)图表示在左图基础上，沿y轴进行纯旋转所得自标定场景；

图3为仿真场景图；其中(a)图表示经过纯旋转运动的自标定场景；(b)图表示空间特征点的像素坐标（圆形点表示摄像机在原位置的图像像素点，菱形点表示纯旋转后的图像像素点），(c)图表示摄像机在原位置的球面投影点，(d)图表示纯旋转后的球面投影点。

图4为带有图像噪声的标定误差仿真结果；其中(a)图表示不同图像噪声水平下的u₀误差，(b)图表示不同图像噪声水平下的v₀误差，(c)图表示不同图像噪声水平下的f_x误差，(d)图表示不同噪声水平下的f_y误差；

图5为带有平移噪声的标定误差仿真结果；其中(a)图表示不同平移噪声水平下的u₀误差，(b)图表示不同平移噪声水平下的v₀误差，(c)图表示不同平移噪声水平下的f_x误差，(d)图表示不同噪声水平下的f_y误差；

图6表示本发明自标定方法的实验场景；

图7表示本发明自标定方法的实验结果；其中(a)图表示实验中的u₀结果，(b)图表示实验中的v₀结果，(c)图表示实验中的f_x结果，(d)图表示实验中的f_y结果；

图8表示第1,6,11次实验中的图像,其中(a)为实验1中第1幅图像,(b)为实验1中第2幅图像,(c)为实验6中第1幅图像,(d)为实验6中第2幅图像,(e)为实验11中第1幅图像,(f)为实验111中第2幅图像.

具体实施方式：

实施例1：

第1，构造针孔摄像机的球形投影模型

定义点P_i，P_j分别表示第i，j个空间点。像素坐标系的横坐标轴与纵坐标轴分别以u，v表示。以表示摄像机坐标系，其中的原点在摄像机光心位置，的z轴与摄像机光轴重合，x轴方向与u轴方向相同，y轴方向与v轴方向相同。f为摄像机焦距，f的单位为米；cp_i，cp_j表示点P_i，P_j对应的图像像素点在下的位置。表示以的原点为球心的单位虚拟球面；s_i，s_j分别为cp_i，cp_j对应在上的投影点，称其为球面投影点；

对于作纯旋转运动的摄像机，与分别表示摄像机在参考位姿处与经过纯旋转运动后的坐标系；s_i，s_j与s_i′，s_j′分别表示点P_i，P_j在与下的球面投影点；关于s_i，s_j与s_i′，s_j′有定理1所述性质：

定理1:球面投影点之间向量的模长在摄像机作纯旋转运动时不变，如式(1)所示：

||s_i-s_j||₂＝||s_i′-s_j′||₂    (1)

定理1呈现了球面投影点在纯旋转摄像机下的性质，根据摄像机的针孔投影模型，其内参数包括f_x，f_y，u₀，v₀；其中f_x，f_y分别为焦距对应于u，v方向的像素块个数，即f_x＝f/d_x，f_y＝f/d_y；其中d_x，d_y分别为单个像素块在u，v方向的长度，单位为米；(u₀，v₀)为图像主点坐标；因此，本发明的目的为根据空间特征点并利用定理1，对纯旋转运动下的摄像机作内参数自标定；

第2，纯旋转运动下的摄像机内参数自标定

第2.1，构造约束方程组

首先，推导并构造了含有摄像机内参数的约束方程为：

$\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}=\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}}---\left(9\right)$

其中

$\left\{\begin{array}{c}{a}_{\mathrm{ij}}=(u_i-u_0)(u_j-u_0)+(v_i-v_0)(v_j-v_0){\mathrm{\γ}}^2\\ b_{\mathrm{ij}}={(u_i-u_0)}^2+{(v_i-v_0)}^2{\mathrm{\γ}}^2;\\ c_{\mathrm{ij}}={(u_j-u_0)}^2+{(v_j+v_0)}^2{\mathrm{\γ}}^2\\ l_{\mathrm{ij}}=({u_i}^{\′}-u_0)({u_j}^{\′}-u_0)+({v_i}^{\′}-v_0)({v_j}^{\′}-v_0){\mathrm{\γ}}^2\\ m_{\mathrm{ij}}={({u_i}^{\′}-u_0)}^2+{({v_i}^{\′}-v_0)}^2{\mathrm{\γ}}^2;\\ {{n_{\mathrm{ij}}=({u_j}^{\′}-u_0)}^2+{({v_j}^{\′}-v_0)}^2\mathrm{\γ}}^2\end{array}\right.$

其中定义了像素块在v方向与u方向的长度比(u_i，v_i)，(u_j，v_j)分别为点P_i与P_j对应的图像像素坐标，(u_i′，v_i′)，(u_j′，v_j′)分别为点P_i与P_j经过摄像机纯旋转之后对应的图像像素坐标；

然后可以利用4个空间点分别得到6个公式(9)形式的约束方程，构成约束方程组，其中设定任意两空间点与摄像机光心不共线；

第2.2，利用非线性最小二乘算法对约束方程组进行求解

采用Levenberg-Marquardt(LM)非线性最小二乘方法进行数值最优化求解；利用至少4个空间点最小化如下的目标函数J_whole(·)，得到u₀，v₀，f_x，γ的解：

$J_{\mathrm{whole}}(u_0,v_0,f_x,\mathrm{\γ})=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{(\frac{a_{\mathrm{ij}}+{f_x}^2}{\sqrt{b_{\mathrm{ij}}+{f_x}^2}\sqrt{c_{\mathrm{ij}}+{f_x}^2}}-\frac{l_{\mathrm{ij}}+{f_x}^2}{\sqrt{m_{\mathrm{ij}}+{f_x}^2}\sqrt{n_{\mathrm{ij}}+{f_x}^2}})}^2---\left(20\right)$

最后利用f_y＝f_x/γ得到f_y，进而完成纯旋转运动下的摄像机内参数自标定。

第3，仿真与实验效果描述

第3.1，仿真结果

本节在存在图像噪声与平移噪声情况下，对本文算法进行仿真验证。对于图像分辨率为740×582像素的待标定摄像机，设定其真实内参数为：

$A=\left[\begin{array}{ccc}1003.1& 0& 369.8\\ 0& 995.4& 306.3\\ 0& 0& 1\end{array}\right]---\left(36\right)$

4个空间点在下的坐标分别为：

(0.4，0.3，1.5)，(0.1，0.2，1.5)，(0.4，-0.3,1.3)，(0.2，-0.15,1.4)    (37)

图像像素点获取过程如下：首先使摄像机在姿态下拍摄第一幅图像，然后使摄像机先绕对应的y轴旋转π/49rad，再绕旋转后的x轴转π/19rad，此时的摄像机坐标系即为继而在姿态下拍摄第二幅图像，如附图3的(a)图所示。

空间点对应的图像像素点如附图3的(b)图所示。其中圆形点与菱形点分别为第一幅与第二幅图像的图像点。附图3的(c)图与(d)图分别对应了摄像机在原始姿态与旋转后姿态下的空间点对应的球面投影点(该两图中的坐标轴为比例值，因而无单位)。易知该4幅图均在内参数假定已知的情况下作出。

在仿真中，4个特征点的图像坐标分别加入了均值为0，幅值为0~1个像素的均匀分布白噪声。为了更准确地反映噪声对实际自标定的影响，对于每个噪声水平，都进行了1000次独立的实验，然后用误差的平均值来对算法的性能进行评估。所得结果如附图4所示。可知在1个像素的噪声之内，4个内参数的相对误差最高为6%，因此得到的结果精度较高。由于在实际场景中，像素坐标的噪声水平一般不会超过0.5个像素，且所用特征点一般大于4个，因而可以得到精度更高的标定结果。

在实际应用中，当云台等机构控制摄像机作纯旋转运动时，经常会伴有非常小的平移量(本发明称其为平移噪声)。附图5给出了在均值为0.0，幅值为0.00-0.20cm的均匀分布平移噪声下的标定误差。由于仿真验证中的空间点的深度在1.4m左右。易知若采用深度更大的空间点，则对平移噪声鲁棒性更强，而这样的空间点在实际场景中是易于找到的。

第3.2，实验结果

在实验中，待标定摄像机采用大恒公司的SV400FC型摄像头，通过IEEE1394接口线与上位机相连，其图像分辨率与仿真中的相同。之后在Visual Studio 2005与OpenCV2.0.0环境下进行了编程实现。

附图6为本文自标定方法的实验场景，空间特征点采用图上所示的4个平面特征点(非平面亦可)。在实验中利用手动使摄像机做纯旋转运动，实验中的三脚架提供了摄像机位置的参考，减少因手动旋转而给摄像机带来的平移噪声.

为了测试图像噪声与平移噪声对标定结果的影响，共进行了3组对比实验，每组实验均进行了5次，实验结果如附图7所示。其中第1~5次实验表示大深度且大旋转时的标定结果(图中星形)，第6~10次实验对应小深度且小旋转情况(图中三角形)，第11~15次实验表示小深度且大旋转情况(图中方块)。附图8对于这三组对比实验，分别给出了第1、6、11次实验中对应的图像。

为检验本文方法的精度，以张正友平面标定法的结果为参考值并作比较，该平面标定法的结果如附图7中粗虚线所示。由附图7可以看出，在大深度且大旋转情况下，算法对噪声鲁棒性更强，易知这是由于在该种情况下，图像与平移噪声分别对于特征点像素值与深度的相对值减小了，因而呈现了更好的算法性能。表1给出了该情况下的标定误差，其中对应于u₀,v₀,f_x,f_y的平均误差分别为0.83%，1.49%，2.47%，0.28%，因而呈现出很好的标定精度。

值得指出的是，本文使摄像机在手持情况下进行自标定，而这会给摄像机带来相当大的平移噪声，然而本文算法仍能获得很好的实验结果，可见本文方法简单实用且具有标定精度高的性能。

表1：第一组实验的标定误差

参考文献

1.Tsai R.Y.,A versatile camera calibration technique for high-accuracy 3D machine visionmetrology using off-the-shelf TV cameras and lenses,IEEE J.Robot.Autom.,1987,3(4):323-344.

2.Zhang Z.,A flexible new technique for camera calibration,IEEE Trans.Pat.Ana.Mach.Intell.,2000,22(11):1330-1334.

3.Faugeras O.,Luong Q.T.,Maybank S.,Camera self-calibration:theory and experiments,inProc.2nd Euro.Conf.Comput.Vis.,1992:321-334.

4.Hartley R.,Zisserman A.,Multiple view geometry in computer vision(2nd Edition),UK:Camb.Univ.Press,2004.

5.孟晓桥,胡占义,摄像机自标定方法的研究进展,自动化学报,2003,29(1):110-124.

6.Pollefeys M.,Gool L.V.,Stratified self-calibration with the modulus constraint,IEEE Trans.Pat.Ana.Mach.Intell.,199921(8):707-724.

7.Hartley R.I.,Self-calibration of stationary cameras,Int.J.Computer Vis.,1997,22(1):5-23.

8.Ma S.,A self-calibration technique for active vision systems,IEEE Trans.Robot.Autom.,1996,12(1):114-120.

9.吴福朝,李华,胡占义,基于主动视觉系统的摄像机自定标方法研究,自动化学报,2001,27(6):752-762.

10.杨长江,孙凤梅,胡占义,“基于二次曲线的纯旋转摄像机自标定”,自动化学报,2001,27(3):310-317.

11.Wang L.,Kang S.B.,Shum H.–Y.,Xu G.,Error analysis of pure rotation-basedself-calibration,IEEE Trans.Pat.Ana.Mach.Intell.,2004,26(2):275-280.

12.Zhang H.,Wong K.–Y.K.,Self-calibration of turntable sequences from silhouettes,IEEETrans.Pat.Ana.Mach.Intell.,2009,31(1):5-14.

13.方勇纯,刘玺,李宝全,张雪波,一种全局指数收敛的摄像机内参数观测器,控制理论与应用,2011,28(9):1082-1090.

14.Geyer C.,Daniilidis K.,Catadioptric Projective Geometry，Int.J.Comput.Vis.,2001,45(3):223–243.

15.Mariottini G.L.and Prattichizzo D.,Image-based visual servoing with central catadioptriccamera,Int.J.Robot.Res.,2008,27(1):41–56.

16.Becerra H.M.,Lopez-Nicolas G.,Sagues C.,Omnidirectional visual control of mobilerobots based on the 1D trifocal tensor,Robot.Auton.Sys.,2010,58(6):796-808.

17.Fomena R.T.,Tahri O.,Chaumette F.,Distance-based and orientation-based visual servoingfrom three points,IEEE Trans.Robot.,2011,27(2):256-267.

18.NocedalJ.,Wright S.J.,Numerical Optimization,New York:Springer,1999.

19.Luca A.D.,Oriolo G.,Giordano P.R..Feature depth observation for image-based visualservoing:theory and experiments,Int.J.Robot.Res.,2008,27(10):1093-1116。

Claims (1)

1.一种基于球形投影模型的纯旋转摄像机自标定方法，其特征在于该方法包括：

第1，构造针孔摄像机的球形投影模型

定义P _i ，P _j 分别表示第i，j个空间点；像素坐标系的横坐标轴与纵坐标轴分别以u，v表示；以表示摄像机坐标系，其中的原点在摄像机光心位置，的z轴与摄像机光轴重合，x轴方向与u轴方向相同，y轴方向与v轴方向相同；f为摄像机焦距，f的单位为米；cp _i ，cp _j 表示点P _i ，P _j 对应的图像像素点在 下的位置； 表示以的原点为球心的单位虚拟球面；s _i ，s _j 分别为cp _i ，cp _j 对应在 上的投影点，称其为球面投影点；

对于作纯旋转运动的摄像机， 与 分别表示摄像机在参考位姿处与经过纯旋转运动后的坐标系；s _i ，s _j 与s _i '，s _j '分别表示点P _i ，P _j 在 与 下的球面投影点；关于s _i ，s _j 与s _i '，s _j '有定理1所述性质：

定理1：球面投影点之间向量的模长在摄像机作纯旋转运动时不变，如式(1)所示：

||s _i -s _j ||₂＝||s _i '-s _j '||₂    (1)

定理1呈现了球面投影点在纯旋转摄像机下的性质，根据通用的摄像机的针孔成像模型，其内参数包括f _x ，f _y ，u ₀，v ₀；其中f _x ，f _y 分别为焦距对应于u，v方向的像素块个数，即f _x ＝f/d _x ，f _y ＝f/d _y ；其中d _x ，d _y 分别为单个像素块在u，v方向的长度，单位为米；(u₀，v₀)为图像主点坐标； 

第2，纯旋转运动下的摄像机内参数自标定

第2.1，构造约束方程组

首先，推导并构造了含有摄像机内参数的约束方程为：

其中

其中定义了像素块在v方向与u方向的长度比 (u _i ,v _i )，(u _j ,v _j )分别为点P _i 与P _j 对应的在像素坐标系下的坐标，(u _i ',v _i ')，(u _j ',v _j ')分别为点P_i与P_j经过摄像机纯旋转之后对应的在像素坐标系下的坐标；

然后利用4个空间点能够得到6个公式(9)形式的约束方程，构成约束方程组，其中设定任意两空间点与摄像机光心不共线；

第2.2，利用非线性最小二乘算法对约束方程组进行求解

采用Levenberg-Marquardt非线性最小二乘方法进行数值最优化求解；利用至少4个空间点最小化如下的目标函数J_whole(·)，得到u ₀，v ₀，f _x ，γ的解：

最后利用f _y ＝f _x /γ得到f _y ，进而完成纯旋转运动下的摄像机内参数自标定。