CN105096328A - 一种基于ios-pso的摄像机自标定方法 - Google Patents

Abstract

本发明公开了一种基于IOS-PSO的摄像机自标定方法，该方法不仅算法比较简单、效率较高，而且寻优结果精度较高。主要包括三个步骤：S1、建立摄像机基于内参数的代价函数；S2、IOS-PSO优化代价函数；S3、IOS-PSO搜索结果限定规则。该方法的具体步骤为：采集同一场景的不同角度的四幅图片；进行图像特征点提取、匹配；利用匹配点对，采用RANSAC方法求取三个基础矩阵F，然后通过基础矩阵的SVD分解，构造代价函数作为粒子群算法的适应值函数；对设定的初始搜索区间划分，并按照最优值的取值情况设定搜索规则；将要标定的摄像机五个未知内参数f_u,f_v,u₀,v₀,s作为粒子群粒子的元素；采用IOS-PSO优化搜索；若搜索结果符合搜索规则，则取最终搜索结果，否则重新搜索。通过四种算法的对比实验，验证了本发明标定结果精确、稳定、高效。

Description

一种基于IOS-PSO的摄像机自标定方法

技术领域

本发明涉及粒子群优化算法，属于计算机视觉领域，具体涉及一种基于IOS-PSO的摄像机自标定方法。

背景技术

摄像机标定，是从摄像机所拍摄的平面二维图像重建物体在空间中的三维信息所必不可少的步骤。假设空间中某物体上的一点A在摄像机所拍摄的平面图像的对应点为a,那么二者存在如下的线性关系：[像点a]＝K[Rt][空间点A]。其中R是3*3矩阵，t是3*1矩阵，分别是世界坐标系到摄像机坐标系的旋转矩阵与平移矩阵；K是摄像机内参数矩阵，确定K的过程即为摄像机标定。

现有摄像机自标定方法利用同一场景的若干幅图像之间的代数几何约束来解算摄像机内参数，可以实现在线标定，比较灵活。本发明采用的自标定方法是基于Kruppa方程的自标定方法，主要有两种，一种是采用数学规划法求解Kruppa方程，这类方法主要通过建立对Kruppa方程中未知参数的约束关系来进行求解，这类方法由于最小点附近呈现大范围的平坦区，故很难找到最小点且目标函数的最小点易受图像噪声的干扰，如SezginTM,DavisR的《HandlingOvertracedStrokesinHand-DrawnSketches》。另一类是利用SVD分解直接求解Kruppa方程的方法，这类方法由于Kruppa方程未知比例因子的存在，因此是一个非线性优化问题。胡华亮，江泽涛的《基于遗传算法的摄像机自标定方法》中，假定了摄像机内参数矩阵中的倾斜因子与比例因子为理想值，这样建立的代价函数自变量为三维，搜索复杂性降低；在非理想情况下，实际搜索变量为五维，经过实验验证发现，遗传算法与粒子群算法虽然在高维优化搜索时均会出现搜索结果的震荡性，但遗传算法中的变异、交叉会增加种群的多样性，这使得遗传算法的搜索结果更加不稳定，不精确。耿利川等的《ANewCameraSelf-calibrationMethodBasedonCSA》算法比较耗时，不利于目标实时跟踪。

粒子群算法(PSO)是一种群智能算法，它模仿鸟类的捕食行为进行全局寻优。不仅规则简单容易实现而且能比较灵活的应用于函数寻优，能很迅速的收敛到全局最优。它从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质。但是，粒子群搜索后期易陷入局部极小，精度差。

PSO算法如下：

1)在d维空间中初始化Par个粒子，第i个粒子初始化为：

${p_i}^g={\left(\begin{array}{cccccc}{p}_{i1}& p_{i2}& ...& p_{ik}& ...& p_{id}\end{array}\right)}^T$

第i个粒子的个体最优：

${{\mathrm{iBest}}_i}^g={\left(\begin{array}{cccccc}{p}_{I1}& p_{I2}& ...& p_{Ik}& ...& p_{Id}\end{array}\right)}^T$

第i个粒子的速度：

${v_i}^g={\left(\begin{array}{cccccc}{v}_{i1}& v_{i2}& ...& v_{ik}& ...& v_{id}\end{array}\right)}^T$

其中g是迭代次数。

2)计算适应值，比较得出全局最优：

gBest^g＝(P₁P₂...P_k...P_d)^T

3)更新粒子:

$\begin{array}{c}{v_{id}}^{g+1}={{\mathrm{wv}}_{id}}^g+c_1{r}_1\left({{\mathrm{iBest}}_{id}}^g-{p_{id}}^g\right)+c_2{r}_2\left({{\mathrm{gBest}}_d}^g-{p_{id}}^g\right)\\ {p_{id}}^{g+1}={p_{id}}^g+{v_{id}}^{g+1}\end{array}$

其中w是非负值，决定前一速度对当前速度的影响；c₁,c₂用于调整粒子在运动中自身经验与群体经验的权重；r₁,r₂是[0,1]上的任意值。

4)判断全局最gBest是否满足规定的全局最优值或是达到规定的迭代次数。否则继续寻优。

但是将PSO应用于本发明建立的代价函数，优化结果会出现解的大幅度震荡性。以代价函数其中一个独立自变量f_u的搜索结果为例说明：

设f_u的定义域为[0,f_u]，当用PSO算法在上述定义域进行搜索时，出现波动，如图4。即，在区间[0,f_u]上进行搜索寻优时，每次运行结果变化幅度较大，可能会收敛到误差较大的P₁,P₂,P₃。

发明内容

有鉴于此，本发明提供了一种基于IOS-PSO的摄像机自标定方法。该方法不仅算法比较简单、效率较高，而且寻优结果精度较高、稳定。

为了解决上述技术问题，本发明是这样实现的：

一种基于IOS-PSO的摄像机自标定方法，包括如下步骤：

步骤1、建立摄像机基于内参数的代价函数f_cost＝f(f_u,f_v,u₀,v₀,s)，其中f_u和f_v分别表示以像素为单位的图像平面坐标系中x方向和y方向的等效焦距；u₀和v₀是图像平面的主点坐标；s是倾斜因子；

步骤2、将f_u和f_v各自的搜索区间分别均分为N个搜索单元，f_u的搜索单元与f_v的搜索单元一一对应；每一个f_u的搜索单元与对应的f_v搜索单元以及u₀、v₀和s的搜索区间形成一个子区间block，则共形成N个子区间block；N为大于或等于4的整数；

步骤3、子区间并行叠进式搜索；

从第一个子区间block开始，数M个子区间block，形成第一个搜索子区间x-search；从第二子区间block开始，数M个子区间block，形成第二个搜索子区间x-search；以此类推，形成多个搜索子区间x-search；M为大于或等于3的整数；

步骤4、采用粒子表征摄像机的5个内参数f_u、f_v、u₀、v₀和s，在每个搜索子区间x-search的各个子区间block中分别初始化粒子，然后同时采用粒子群优化算法对各子区间block进行搜索，得到各子区间block的局部最优粒子以及对应的适应值；

步骤5、步骤4中获得的适应值最小的局部最优粒子则为搜索子区间x-search的局部最优粒子；按照步骤4的方法获得各个搜索子区间x-search的局部最优粒子，再确定最小适应值及对应的局部最优粒子，即为整个搜索区间对应的全局最优粒子；

步骤6、将该全局最优粒子对应的内参数代入到所述摄像机的代价函数中，则完成对摄像机的标定。

较佳的，在所述步骤4得到各子区间block的局部最优粒子以及对应的适应值后，判断所述各个搜索子区间x-search的各子区间block的适应值是否同时单调递减或单调递增：如果是，返回所述步骤4，对各个子区间block中粒子重新进行初始化，搜索最优子区间，然后继续依次执行步骤5和6；否则，直接依次执行步骤5和步骤6。

较佳的，在所述步骤4得到各子区间block的局部最优粒子以及对应的适应值后，当所述各个搜索子区间x-search的局部最优粒子的适应值不满足均为单调递减或单调递增关系时，所述步骤5中，获得全局最优粒子后，先判断全局最优粒子是否满足如下条件：

条件1、当全局最优粒子落在一些搜索子区间x-search的重叠区域时，并且，其它没有重叠区域的搜索子区间x-search对应的局部最优粒子所在子区间block应为所在搜索子区间x-search中最靠近所述全局最优粒子的子区间block；

条件2、当全局最优粒子只落在其中一个搜索子区间x-search中时，其它不包含全局最优粒子的搜索子区间x-search对应的局部最优粒子所在子区间block应为所在搜索子区间x-search中最靠近所述全局最优粒子的子区间block；

如果满足条件1和条件2中的任意一个，则采用所述全局最优粒子的内参数进行摄像机标定；否则，返回所述步骤4，对各个个子区间block中粒子重新进行初始化，搜索局部最优粒子，然后继续执行步骤5和步骤6。

本发明具有如下有益效果：

(1)本发明所提供的方法，先建立代价函数，再根据相机参数确定5个未知量的搜索区间，对搜索区间进行子区间划分，在各子搜索区间上u,v,s参数的区间保持一致；采用改进的粒子群寻优方法IOS-PSO采用子区间并行叠进式搜索，确定搜索区间局部最优粒子，最终确定全局最优粒子，方法原理简单，易于实现，搜索结果精度高、稳定，同时又能弥补粒子群算法高维寻优时的不足。

(2)每次搜索确定全局最优粒子后，根据搜索规则判定最优粒子是否满足精度要求，如果满足，可根据最优粒子对相机进行标定，如果不满足，继续采用粒子群优化算法搜索最优粒子，如此，可避免搜索粒子的波动性，使得摄像机标定结果较PSO精度较大幅度提高。

附图说明

图1为本发明所提供的方法的流程图；

图2为粒子群算法流程图；

图3为粒子群分区间并行叠进搜索流程图；

图4为原始粒子群算法搜索情况；

图5为粒子群分区间并行叠进搜索情况；

图6为粒子群分区间并行叠进搜索方式；

图7为粒子群分区间并行叠进搜索方式具体实施方式；

图8为二视图几何；

图9为本实验所用相机图；

图10为拍摄方式；

图11为本发明实验拍摄的场景图片；

图12为本发明提供的四种算法标定结果得到的距离相对误差分布。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明的标定方法包括如下步骤：

步骤1、建立摄像机基于内参数的代价函数；同一个场景的两幅图像对应的像点满足：

x'Fx＝0(0.1)

其中F为3×3的基本矩阵，秩为2。

经过推导求得如下的Kruppa方程：

${\mathrm{FCF}}^T=c{\[e^{\′}\]}_{\×}C{\[e^{\′}\]}_{\×}^T---\left(0.2\right)$

其中 $C={\mathrm{KK}}^T,K=\left(\begin{array}{ccc}{f}_u& k& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right),$ e'为两幅图像中x'所在图像的极点；[.]_×代表反对称矩阵；c为未知的正比例因子。

将矩阵F进行SVD分解的结果，

$F={\mathrm{UDV}}^T=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\σ}}_i^F{u}_i{v}_i^T---\left(0.3\right)$

其中是F的第i个奇异值；u_i和v_i分别是左右奇异向量。

由基础矩阵与极点的关系得：

F^Te'＝0(0.4)

将(1.3)(1.4)带入式(1.2)，得到：

$\begin{array}{c}{f}_1=\frac{{\left({{\mathrm{\σ}}_1}^F\right)}^2\[k_1{v_{11}}^2+k_r{v_{21}}^2+{v_{31}}^2+2\left(k_2{v}_{11}{v}_{21}+k_3{v}_{11}{v}_{31}+k_5{v}_{21}{v}_{31}\right)\]}{\left(k_1{v_{12}}^2+k_4{u_{22}}^2+{u_{31}}^2\right)+2\left(k_2{u}_{12}{u}_{22}+k_3{u}_{12}{u}_{31}+k_5{u}_{22}{u}_{32}\right)}\\ f_3=\frac{{\left({{\mathrm{\σ}}_2}^F\right)}^2\[k_1{v_{12}}^2+k_4{v_{22}}^2+{v_{32}}^2+2\left(k_2{v}_{12}{v}_{22}+k_3{v}_{12}{v}_{32}+k_5{v}_{22}{v}_{32}\right)\]}{\left(k_1{v_{11}}^2+k_4{u_{21}}^2+{u_{31}}^2\right)+2\left(k_2{u}_{11}{u}_{21}+k_3{u}_{11}{u}_{31}+k_5{u}_{21}{u}_{31}\right)}\\ f_2=\frac{{{\mathrm{\σ}}_1}^F{{\mathrm{\σ}}_2}^F\[k_1{v}_{11}{v}_{12}+k_4{v}_{21}{v}_{22}+v_{31}{v}_{32}+k_2\left(v_{11}{v}_{22}+v_{21}{v}_{12}\right)+k_3\left(v_{11}{v}_{32}+v_{31}{v}_{12}\right)+k_5\left(v_{21}{v}_{32}+v_{31}{v}_{22}\right)\]}{-\[k_1{u}_{11}{u}_{12}+k_4{u}_{21}{u}_{22}+u_{31}{u}_{32}+k_2\left(u_{11}{u}_{22}+u_{21}{u}_{12}\right)+k_3\left(u_{11}{u}_{32}+u_{31}{u}_{12}\right)+k_5\left(u_{21}{u}_{32}+u_{31}{u}_{22}\right)\]}\end{array}---\left(0.5\right)$

其中k1＝f_u ²+u₀ ²,k2＝u₀v₀,k3＝u₀,k4＝f_v ²+v₀ ²,k5＝v₀；

Kru个Kruppa方程可以建立如下代价函数：

$\begin{array}{c}{f}_{\cos t}={\left(f_{11}-f_{12}\right)}^2+{\left(f_{11}-f_{13}\right)}^2+{\left(f_{12}-f_{13}\right)}^2+\\ {\left(f_{21}-f_{22}\right)}^2+{\left(f_{21}-f_{23}\right)}^2+{\left(f_{22}-f_{23}\right)}^2+\\ \mathrm{......}\\ {\left(f_{Kru1}-f_{Kru2}\right)}^2+{\left(f_{Kru1}-f_{Kru3}\right)}^2+{\left(f_{Kru2}-f_{Kru3}\right)}^2\end{array}---\left(0.6\right)$

本发明中，进行摄像机标定获得5个未知参数，而一个Kruppa方程仅提供关于未知参数的两个独立约束，所以取Kru值为3。

步骤2、IOS-PSO优化代价函数；

此方法是应用于照相机自标定，将代价函数(1.6)表示为f_cost＝f(f_u,f_v,u₀,v₀,s)，其中f_u,f_v分别表示以像素为单位的图像平面坐标系中x方向和y方向的等效焦距；u₀,v₀是图像平面的主点坐标，以像素为单位；s是倾斜因子，该倾斜因子值近似为0。其中f_u,f_v的搜索区间分别设为因为现在相机的设计制造工艺比较高，所以主点坐标近似为摄像机的像素中心点，它们的搜索区间可以定义的更精确[u₁,u₂],[v₁,v₂]；s的搜索区间设为[0,s]。本发明中所用相机满足所以设定中的

1)IOS-PSO区间划分；

将f_u,f_v的搜索区间分别平均分为N个搜索单元，而N个搜索单元各自对应的u,v,s的搜索区间[u₁,u₂],[v₁,v₂],[0,s]保持一致，分后如下：

区间1：[u₁,u₂],[v₁,v₂],[0,s]

区间2：[u₁,u₂],[v₁,v₂],[0,s]

......

区间N：[u₁,u₂],[v₁,v₂],[0,s]

2)子区间并行叠进式搜索；

以一维变量f_u的搜索为例，搜索方式如图6。

整个搜索过程进行了m次X-search1,X-search2,…,X-searchm，任意一次搜索X-search所包含的子区间个数为n(n>＝3)任意相邻的两次X-searchx,X-search(x+1),后者相对于前者前进了一个子区间。

X-search1:同时进行搜索，在X-search1的n个子区间分别初始化par个粒子，然后同时进行PSO优化搜索，直到满足收敛条件，得到n个子区间的局部最优以及对应的适应值，分别为对应的最优区间1_Block_x；

X-search2:同时进行搜索，在X-search2的n个子区间分别初始化par个粒子，然后同时进行PSO优化搜索，直到满足收敛停止条件，得到n个子区间的局部最优以及对应的适应值，分别为最优区间2-Block_y。

……

X-searchm:同时进行搜索，在X-searchm的n个子区间分别初始化par个粒子，然后同时进行PSO优化搜索，直到满足收敛停止条件，得到n个子区间的局部最优以及对应的适应值，分别为最优区间m-Block_y。

其它的X-search按照上述方式逐步叠进搜索,得到m个最优子区间。可以看到，X-search1与X-search2搜索时有两个区间是重叠的。叠进搜索的原理在于若X-search1的重叠区域搜索出现误差；那么X-search2进一步搜索的同时，对重叠区域重复搜索，可以对之前的误差进行校正，同理，X-search1也可以对X-search2进行校正。并且，在重叠区域两次搜索的搜索结果都正确的情况下，若X-search2的搜索结果优于X-search1，可以替换X-search1的结果，使结果更优化。所以n至少是3，若小于3重叠区域起不到相互更正误差的作用，当然随着n以及m的增大(区间分的越精细)，耗时是随着增加的。

步骤3、IOS-PSO搜索结果限定规则；

首先，判断此次搜索是否为Abnor-1或Abnor-2；若是则跳至S2，否则继续判断是否为Situation1，Situation2中的任何一种情况，若满足，则此次搜索是正确的，然后取m次搜索得到的m个最优子区间中的最优者作为全局最优；否则跳至S2重新进行搜索。

1)IOS-PSO搜索结果限定；

根据最优值所在区间假设，本发明设定了优化结果的两种情况。

Situation1:当其中一些搜索子区间search具有重叠区域时，该全局最优粒子应落在重叠区域中，并且，其它没有重叠区域的搜索子区间X-search对应的局部最优粒子所在子区间block应为所在搜索子区间X-search中最靠近所述全局最优粒子的子区间block；这样的搜索结果是正确的。

Situation2:当全局最优粒子只落在其中一个搜索子区间X-search中时，其它不包含全局最优粒子的搜索子区间X-search对应的局部最优粒子所在子区间block应为所在搜索子区间X-search中最靠近所述全局最优粒子的子区间block；这样的搜索结果是正确的。

2)IOS-PSO异常结果；

X-search1:经过第一次搜索，个子区间的全局最优对应的适应值，依次设为

X-search2:经过第二次搜索，个子区间的全局最优对应的适应值，依次设为

......

X-searchm:经过第m次搜索，个子区间的全局最优对应的适应值，依次设为

Abnor-1:且且

Abnor-2:且且

实施例：

一种基于IOS-PSO的摄像机自标定方法，该方法源于对独立五参数代价函数的优化搜索。主要包括三个步骤：S1、建立摄像机基于内参数的代价函数；S2、IOS-PSO优化代价函数；S3、IOS-PSO搜索结果限定规则。该方法的具体介绍为：

S1.建立摄像机基于内参数的代价函数；

如图8所示，左图中任意一点x在右图中都有相应的一条对极线I'，二者的映射关系表示为x→I'，引入基础矩阵F，将此映射表示为：

I'＝Fx(0.7)

因为左图中点x在右图中的对应点x'位于对极线I'上，所以同一个场景的两幅图像对应的像点满足：

x'Fx＝0(0.8)

其中F为3×3的基本矩阵，秩为2。

经过推导求得如下的Kruppa方程：

${\mathrm{FCF}}^T=c{\[e^{\′}\]}_{\×}C{\[e^{\′}\]}_{\×}^T---\left(0.9\right)$

其中 $C={\mathrm{KK}}^T,K=\left(\begin{array}{ccc}{f}_u& k& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right);$ e'为两幅图像中x'所在图像的极点；[.]_×代表反对称矩阵；c为未知的正比例因子。

消去比例因子c,可以得到如下等比方程：

$\begin{array}{c}\frac{{\left({\mathrm{FCF}}^T\right)}_{11}}{{\left({\[e^{\′}\]}_{\×}C{{\[e^{\′}\]}_{\×}}^T\right)}_{11}}=\frac{{\left({\mathrm{FCF}}^T\right)}_{12}}{{\left({\[e^{\′}\]}_{\×}C{{\[e^{\′}\]}_{\×}}^T\right)}_{12}}=\\ \frac{{\left({\mathrm{FCF}}^T\right)}_{13}}{{\left({\[e^{\′}\]}_{\×}C{{\[e^{\′}\]}_{\×}}^T\right)}_{13}}=\frac{{\left({\mathrm{FCF}}^T\right)}_{22}}{{\left({\[e^{\′}\]}_{\×}C{{\[e^{\′}\]}_{\×}}^T\right)}_{22}}=\\ \frac{{\left({\mathrm{FCF}}^T\right)}_{23}}{{\left({\[e^{\′}\]}_{\×}C{{\[e^{\′}\]}_{\×}}^T\right)}_{23}}=\frac{{\left({\mathrm{FCF}}^T\right)}_{33}}{{\left({\[e^{\′}\]}_{\×}C{{\[e^{\′}\]}_{\×}}^T\right)}_{33}}\end{array}---\left(0.10\right)$

将矩阵F进行SVD分解的结果：

$F={\mathrm{UDV}}^T=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\σ}}_i^F{u}_i{v}_i^T---\left(0.11\right)$

其中是F的第i个奇异值；u_i和v_i分别是左右奇异向量。

如图8所示e'，是右图中的极点，为基础矩阵的右零极点，所以它满足：

F^Te'＝0(0.12)

由(1.12)可得：

$e^{\′}\≅u_3\≅U\·{\left(\begin{array}{ccc}0& 0& 1\end{array}\right)}^T---\left(0.13\right)$ 将(1.11)(1.13)带入(1.9)得：

${\mathrm{DV}}^T{\mathrm{CVD}}^T\≅{\mathrm{MU}}^T{\mathrm{CUM}}^T---\left(0.14\right)$

按照式(1.14)的形式展开得：

$\begin{array}{c}\left(\begin{array}{ccc}{\left({{\mathrm{\σ}}_1}^F\right)}^2{v_1}^T{\mathrm{Cv}}_1& {{\mathrm{\σ}}_1}^F{{\mathrm{\σ}}_2}^F{v_1}^T{\mathrm{Cv}}_2& 0\\ {{\mathrm{\σ}}_1}^F{{\mathrm{\σ}}_2}^F{v_2}^T{\mathrm{Cv}}_1& {\left({{\mathrm{\σ}}_2}^F\right)}^2{v_2}^T{\mathrm{Cv}}_2& 0\\ 0& 0& 0\end{array}\right)\≅\left(\begin{array}{ccc}{u_2}^T{\mathrm{Cu}}_2& -{u_2}^T{\mathrm{Cu}}_1& 0\\ -{u_1}^T{\mathrm{Cu}}_2& {u_1}^T{\mathrm{Cu}}_1& 0\\ 0& 0& 0\end{array}\right)\\ \frac{{\left({{\mathrm{\σ}}_1}^F\right)}^2{v_1}^T{\mathrm{Cv}}_1}{{u_2}^T{\mathrm{Cu}}_2}=\frac{{{\mathrm{\σ}}_1}^F{{\mathrm{\σ}}_2}^F{v_1}^T{\mathrm{Cv}}_2}{-{u_2}^T{\mathrm{Cu}}_1}=\frac{{\left({{\mathrm{\σ}}_2}^F\right)}^2{v_2}^T{\mathrm{Cv}}_2}{{u_1}^T{\mathrm{Cu}}_1}\end{array}---\left(0.15\right)$

令

$\begin{array}{c}{f}_1=\frac{{\left({{\mathrm{\σ}}_1}^F\right)}^2{v_1}^T{\mathrm{Cv}}_1}{{u_2}^T{\mathrm{Cu}}_2},f_2=\frac{{{\mathrm{\σ}}_1}^F{{\mathrm{\σ}}_2}^F{v_1}^T{\mathrm{Cv}}_2}{-{u_2}^T{\mathrm{Cu}}_1}\\ f_3=\frac{{\left({{\mathrm{\σ}}_2}^F\right)}^2{v_2}^T{\mathrm{Cv}}_2}{{u_1}^T{\mathrm{Cu}}_1}\end{array}---\left(0.16\right)$

$\begin{array}{ccc}{v}_1=\left(\begin{array}{c}{v}_{11}\\ v_{12}\\ v_{13}\end{array}\right),& v_2=\left(\begin{array}{c}{v}_{21}\\ v_{22}\\ v_{23}\end{array}\right),& v_3=\left(\begin{array}{c}{v}_{31}\\ v_{32}\\ v_{33}\end{array}\right),\\ u_1=\left(\begin{array}{c}{u}_{11}\\ u_{12}\\ u_{13}\end{array}\right),& u_2=\left(\begin{array}{c}{u}_{21}\\ u_{22}\\ u_{23}\end{array}\right),& u_3=\left(\begin{array}{c}{u}_{31}\\ u_{32}\\ u_{33}\end{array}\right)\end{array}---\left(0.17\right)$

将(1.17)带入式(1.16)，得到：

$\begin{array}{c}{f}_1=\frac{{\left({{\mathrm{\σ}}_1}^F\right)}^2\[k_1{v_{11}}^2+k_r{v_{21}}^2+{v_{31}}^2+2\left(k_2{v}_{11}{v}_{21}+k_3{v}_{11}{v}_{31}+k_5{v}_{21}{v}_{31}\right)\]}{\left(k_1{v_{12}}^2+k_4{u_{22}}^2+{u_{31}}^2\right)+2\left(k_2{u}_{12}{u}_{22}+k_3{u}_{12}{u}_{31}+k_5{u}_{22}{u}_{32}\right)}\\ f_3=\frac{{\left({{\mathrm{\σ}}_2}^F\right)}^2\[k_1{v_{12}}^2+k_4{v_{22}}^2+{v_{32}}^2+2\left(k_2{v}_{12}{v}_{22}+k_3{v}_{12}{v}_{32}+k_5{v}_{22}{v}_{32}\right)\]}{\left(k_1{v_{11}}^2+k_4{u_{21}}^2+{u_{31}}^2\right)+2\left(k_2{u}_{11}{u}_{21}+k_3{u}_{11}{u}_{31}+k_5{u}_{21}{u}_{31}\right)}\\ f_2=\frac{{{\mathrm{\σ}}_1}^F{{\mathrm{\σ}}_2}^F\[k_1{v}_{11}{v}_{12}+k_4{v}_{21}{v}_{22}+v_{31}{v}_{32}+k_2\left(v_{11}{v}_{22}+v_{21}{v}_{12}\right)+k_3\left(v_{11}{v}_{32}+v_{31}{v}_{12}\right)+k_5\left(v_{21}{v}_{32}+v_{31}{v}_{22}\right)\]}{-\[k_1{u}_{11}{u}_{12}+k_4{u}_{21}{u}_{22}+u_{31}{u}_{32}+k_2\left(u_{11}{u}_{22}+u_{21}{u}_{12}\right)+k_3\left(u_{11}{u}_{32}+u_{31}{u}_{12}\right)+k_5\left(u_{21}{u}_{32}+u_{31}{u}_{22}\right)\]}\end{array}---\left(0.18\right)$

其中k1＝f_u ²+u₀ ²,k2＝u₀v₀,k3＝u₀,k4＝f_v ²+v₀ ²,k5＝v₀

3个上述Kruppa方程可以建立如下代价函数：

$\begin{array}{c}{f}_{\cos t}={\left(f_{11}-f_{12}\right)}^2+{\left(f_{11}-f_{13}\right)}^2+{\left(f_{12}-f_{13}\right)}^2+\\ {\left(f_{21}-f_{22}\right)}^2+{\left(f_{21}-f_{23}\right)}^2+{\left(f_{22}-f_{23}\right)}^2+\\ {\left(f_{31}-f_{32}\right)}^2+{\left(f_{31}-f_{33}\right)}^2+{\left(f_{32}-f_{33}\right)}^2\end{array}---\left(0.19\right)$

S2.IOS-PSO优化代价函数；

此方法是应用于照相机的内参数自标定，将式(1.19)表示成f_cost＝f(f_u,f_v,u₀,v₀,s)，其中f_u,f_v分别表示以像素为单位的图像平面坐标系中x方向和y方向的等效焦距；u₀,v₀是图像平面的主点坐标，以像素为单位；s是倾斜因子，该倾斜因子值近似为0。其中本发明中所用相机满足f_u≈f_v，f_u,f_v的搜索区间分别设为[0,2000],[0,2000]；因为现在相机的设计制造工艺比较高，主点近似为摄像机像素中心点，本发明所用相机的分辨率为640*480，它们的搜索区间可以定义的更精确为[310,330],[230,250]；s的搜索区间设为[0,0.01]。

1)IOS-PSO算法区间划分

我们将f_u,f_v的定义区间[0,2000],[0,2000]分别平均分为四个子区间，而四个子区间中的u,v,s的定义区间[310,330],[230,250],[0,0.01]保持一致，分后如下：

区间1：[0,500],[0,500],[310,330],[230,250],[0,0.01]

区间2：[500,1000],[500,1000],[310,330],[230,250],[0,0.01]

区间3：[1000,1500],[1000,1500],[310,330],[230,250],[0,0.01]

区间4：[1500,2000],[1500,2000],[310,330],[230,250],[0,0.01]

2)IOS-PSO算法子区间并行叠进搜索；

按照图7的优化方式进行优化，在X-search1的

区间1：[0,500],[0,500],[310,330],[230,250],[0,0.01]

区间2：[500,1000],[500,1000],[310,330],[230,250],[0,0.01]

区间3：[1000,1500],[1000,1500],[310,330],[230,250],[0,0.01]

并行搜索，每个区间30个粒子；

Step1.第i个粒子初始化为：

${p_i}^g={\left(\begin{array}{ccccc}{f_i}^u& {f_i}^v& u_i& v_i& s_i\end{array}\right)}^T$

第i个粒子的个体最优：

${{\mathrm{iBest}}_i}^g={\left(\begin{array}{ccccc}{f}_{pi}^u& f_{pi}^v& u_{pi}& v_{pi}& s_{pi}\end{array}\right)}^T$

第i个粒子的速度：

${v_i}^g={\left(\begin{array}{ccccc}{v}_i^{fu}& v_i^{fv}& v_i^u& v_i^v& v_i^s\end{array}\right)}^T$

其中g是迭代次数。

Step2.将初始化的粒子带入式(1.19)计算适应值，两个粒子适应值小的为更优个体，比较得出第k个子区间的全局最优：

${\mathrm{gBest}}_{zk}^g={\left(\begin{array}{ccccc}{f}_{zk}^u& f_{zk}^v& u_{zk}& v_{zk}& s_{zk}\end{array}\right)}^T$

Step3.若迭代次数满足本发明，则停止，否则，更新粒子：

$\begin{array}{c}{v_{id}}^{g+1}={{\mathrm{wv}}_{id}}^g+c_1{r}_1\left({{\mathrm{iBest}}_{id}}^g-{p_{id}}^g\right)+c_2{r}_2\left({{\mathrm{gBest}}_d}^g-{p_{id}}^g\right)\\ {p_{id}}^{g+1}={p_{id}}^g+{v_{id}}^{g+1}\end{array}$

其中w是非负值其值取为0.3，决定前一速度对当前速度的影响；c₁,c₂用于调整粒子在运动中自身经验与群体经验的权重，c₁,c₂分别取值0.5,0.1；r₁,r₂是[0,1]上的任意值。

Step4.计算粒子适应值，比较更新个体最优与全局最优。

如果各子区间第g+1代的个体适应值比g代的个体最优适应值大，即则iBest_iF ^g值保持不变iBest_iF ^g+1＝iBest_iF ^g，否则如果第g+1代全局最优比第g代的全局最优适应值大，即则gBest^g保持不变,gBest^g+1＝gBest^g，否则gBest^g+1＝gBest^g+1。跳至Step3继续执行。

Step5.满足收敛条件，得到X-search1三个区间的局部最优以及三个区间的最优搜索结果的适应值，分别为最优区间为1_Block_x；

Step6.在X-search2的三个区间

区间2：[500,1000],[500,1000],[310,330],[230,250],[0,0.01]

区间3：[1000,1500],[1000,1500],[310,330],[230,250],[0,0.01]

区间4：[1500,2000],[1500,2000],[310,330],[230,250],[0,0.01]

分别初始化30个粒子，然后按照步骤Step1、Step2、Step3、Step4同时搜索，直到满足收敛停止条件。设三个区间的最优搜索结果的适应值，分别为最优区间为2_Block_y。

S3.IOS-PSO搜索结果限定规则；

首先，判断此次搜索是否为Abnor-1或Abnor-2；若是则跳至S2，否则继

继续判断(1_Block_x,2-Block_y)是否为表1中的任何一个，若满足，则此次搜索是正确的，然后取六个子区间中的最优者作为全局最优；否则跳至S2重新进行搜索。

1)IOS-PSO搜索结果限定；

根据最优值所在区间假设，本发明设定了优化结果的四种情况，IOS-PSO搜索结果须满足其中的任意一种，如表1；

以其中一种情况为例进行说明：假设最优值位于区间

X-search1:经过第一次搜索，三个区间的全局最优对应的适应值，依次设为其中最优适应值应该是第二个区间的全局最优对应的适应值。记本次搜索结果为1-2(其中第一个数字1代表第一次搜索X-search1；第二个数字2代表本次最优值所在区间为第2个区间)。

X-search2:经过第二次搜索，三个区间的全局最优对应的适应值，依次设为其中最优适应值应该是第一个区间的全局最优对应的适应值。记本次搜索结果为2-2(其中第一个数字2代表第二次搜索X-search2；第二个数字2代表本次最优值所在区间为第2个区间)。

2)IOS-PSO异常结果；

X-search1:经过第一次搜索，三个区间的全局最优对应的适应值，依次设为

X-search2:经过第二次搜索，三个区间的全局最优对应的适应值，依次设为

Abnor-1:且

Abnor-2:且

实验用相机为SonyDev-50，如图9，选取校园内三辆停放在一起的汽车做为本次实验的场景目标。拍摄时，从东向西依次拍摄pic(a)、pic(b)、pic(c)、pic(d)，如图10。拍摄得到的四幅图片，如图11。经过测量,目标物体距离照相机的距离约为43.4m；pic(a)与pic(b)的距离约为2.72m，用于计算目标到相机的距离。相机拍摄时分辨率为640*480。f_u,f_v,u,v,s的搜索区间选取上述[0,2000],[0,2000],[310,330],[230,250],[0,0.01]值。注意拍摄时四幅图像的共同区域尽可能多，以提取尽可能多的匹配点；为了后期计算目标距离拍摄相机的距离，拍摄过程中尽可能保持相机的姿态一致。

首先，从第一幅图像，pic(a)中选定要实际进行定位的目标物体——三辆小汽车部分，并将选定部分记为pic(a)’；然后，对pic(a)’、pic(b)、pic(c)、pic(d)进行特征点提取，并将pic(a)’分别与pic(b)、pic(c)、pic(d)匹配，利用匹配特征点建立约束来分别求取pic(a)’与pic(b)、pic(c)、pic(d)之间的基础矩阵共三个；最后利用具体实施措施S1的方法建立代价函数。

将内参数矩阵的五个未知参数f_u,f_v,u₀,v₀,s作为IOS-PSO的输入参数，将式(1.19)代价函数的取值作为输出值，亦适应值，训练次数设置为6000次。实验结果如表2。因为f_u≈f_v，为了结果更精确，将f_u,f_v的搜索结果取平均f作为最后两者的取值结果。为了说明本算法的高精度，本发明同时提供了粒子群算法(PSO)、遗传算法(GA)、共轭梯度法(CG)的自标定结果，三个算法的标定结果分别为表3、表4、表5。四种算法中，每种算法均运行10次，求解目标与摄像机距离时均带入f_u,f_v的平均值f。图12则提供了四种算法10次标定结果求得的距离相对误差的分布对照

$relativeerror=\frac{calculateddis\tan ce-actualdis\tan ce}{actualdis\tan ce}.$

表2IOS-PSO标定结果

表3PSO标定结果

表4GA标定结果

表5CG标定结果

1.由计算的距离分析精度；

由表2、表3、表4、表5可以看到IOS-PSO、PSO、GA、CG最终标定得到的结果求得的10次距离平均值分别为42.500m、75.9716m、54.927m、63.7969m。而本发明目标的距离约为43.4m。可以得出，IOS-PSO算法的标定精度最高。由图12可以看到四种算法中，IOS-PSO的十次标定误差波动最小，标定结果最稳定。说明本发明公开的算法不仅标定结果稳定，而且误差小，精度高。

2.u₀,v₀,s分析；

由表2、表3、表4、表5可以看到，四种算法10次运算的u₀,v₀,s值均比较平稳精确。

3.f_u,f_v分析；

由表2、表3、表4、表5可以看到，IOS-PSO标定的f_u,f_v值，在10次运算中比较平稳，而后三者震荡比较大。

由图12可以看到，GA算法的波动性最大。这是因为，搜索算法本身在优化高维函数时具有比较的波动性，GA算法中交叉、变异的作用是增加种群的多样性，这就导致GA的搜索更不易收敛。CG算法属于数学迭代，波动比较小，但是它容易收敛到距离最优值稍微远一点的区域；这是因为，CG算法是按照梯度值进行迭代寻优的，而本发明的代价函数的值量级非常小，这就导致了CG算法最终因为寻优步长太小，无法进一步寻优就停止搜索。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于IOS-PSO的摄像机自标定方法，其特征在于，包括如下步骤：

步骤1、建立摄像机基于内参数的代价函数f_cost＝f(f_u,f_v,u₀,v₀,s)，其中f_u和f_v分别表示以像素为单位的图像平面坐标系中x方向和y方向的等效焦距；u₀和v₀是图像平面的主点坐标；s是倾斜因子；

步骤2、将f_u和f_v各自的搜索区间分别均分为N个搜索单元，f_u的搜索单元与f_v的搜索单元一一对应；每一个f_u的搜索单元与对应的f_v搜索单元以及u₀、v₀和s的搜索区间形成一个子区间block，则共形成N个子区间block；N为大于或等于4的整数；

步骤3、子区间并行叠进式搜索；

从第一个子区间block开始，数M个子区间block，形成第一个搜索子区间x-search；从第二子区间block开始，数M个子区间block，形成第二个搜索子区间x-search；以此类推，形成多个搜索子区间x-search；M为大于或等于3的整数；

步骤4、采用粒子表征摄像机的5个内参数f_u、f_v、u₀、v₀和s，在每个搜索子区间x-search的各个子区间block中分别初始化粒子，然后同时采用粒子群优化算法对各子区间block进行搜索，得到各子区间block的局部最优粒子以及对应的适应值；

步骤5、步骤4中获得的适应值最小的局部最优粒子则为搜索子区间x-search的局部最优粒子；按照步骤4的方法获得各个搜索子区间x-search的局部最优粒子，再确定最小适应值及对应的局部最优粒子，即为整个搜索区间对应的全局最优粒子；

步骤6、将该全局最优粒子对应的内参数代入到所述摄像机的代价函数中，则完成对摄像机的标定。

2.如权利要求1所述的一种基于IOS-PSO的摄像机自标定方法，其特征在于，

在所述步骤4得到各子区间block的局部最优粒子以及对应的适应值后，判断所述各个搜索子区间x-search的各子区间block的适应值是否同时单调递减或单调递增：如果是，返回所述步骤4，对各个子区间block中粒子重新进行初始化，搜索最优子区间，然后继续依次执行步骤5和6；否则，直接依次执行步骤5和步骤6。

3.如权利要求2所述的一种基于IOS-PSO的摄像机自标定方法，其特征在于，在所述步骤4得到各子区间block的局部最优粒子以及对应的适应值后，当所述各个搜索子区间x-search的局部最优粒子的适应值不满足均为单调递减或单调递增关系时，所述步骤5中，获得全局最优粒子后，先判断全局最优粒子是否满足如下条件：

条件1、当全局最优粒子落在一些搜索子区间x-search的重叠区域时，并且，其它没有重叠区域的搜索子区间x-search对应的局部最优粒子所在子区间block应为所在搜索子区间x-search中最靠近所述全局最优粒子的子区间block；

条件2、当全局最优粒子只落在其中一个搜索子区间x-search中时，其它不包含全局最优粒子的搜索子区间x-search对应的局部最优粒子所在子区间block应为所在搜索子区间x-search中最靠近所述全局最优粒子的子区间block；

如果满足条件1和条件2中的任意一个，则采用所述全局最优粒子的内参数进行摄像机标定；否则，返回所述步骤4，对各个个子区间block中粒子重新进行初始化，搜索局部最优粒子，然后继续执行步骤5和步骤6。