张正友标定数学原理

张正友标定算法原理详解
张正友标定算法是一种多摄像机标定的有效方法，它使用了一个具有内部参考点的3D物体来检测至少6个相互独立的2D-3D配准，从而确定每个摄像机之间的关系。

该算法的优点在于它能够利用棋盘格子上的多个点，从而使得标定更加可靠。

张正友标定算法主要分为三步：
(1)首先获取所有摄像机的内参和外参，这些参数表示每个摄像机的视角和位置，包括焦距、畸变参数和旋转参数等信息。

(2)计算每个摄像机的平面坐标系，这些平面坐标系将求解的不同摄像机之间的关系。

(3)把实际的物体点映射到每个摄像机的平面坐标系中，并计算所有摄像机之间的关系。

最后，张正友标定算法可以得出每个摄像机的位置和姿态，从而实现多摄像机视角的标定。

张正友标定算法解读一直以来想写篇相机标定方面的东西，最近组会上也要讲标定方面东西，所以顺便写了。

无论是OpenCV还是matlab标定箱，都是以张正友棋盘标定算法为核心实现的,这篇PAMI的文章<<A Flexible New Technique for Camera Calibration>>影响力极大，张正友是zju的机械系出身，貌似现在是微软的终身教授了，有点牛的。

我就简单的介绍下算法的核心原理，公式的推理可能有点多。

一基本问题描述：空间平面的三维点与相机平面二维点的映射假设空间平面中三维点：（齐次坐标，世界坐标系）相机平面二维点：（齐次坐标，相机坐标系）那么空间中的点是如何映射到相机平面上去呢？我们用一个等式来表示两者之间关系：（1）注：A为相机内参矩阵，R，t分别为旋转和平移矩阵，s为一个放缩因子标量。

我们把等式（1）再简化下：（2）因为张正友算法选取的是平面标定，所以令z=0,所以平移向量只有r1,r2即可。

H就是我们常说的单应性矩阵，在这里描述的是空间中平面三维点和相机平面二维点之间的关系。

因为相机平面中点的坐标可以通过图像处理的方式（哈里斯角点，再基于梯度搜索的方式精确控制点位置）获取，而空间平面中三维点可以通过事先做好的棋盘获取。

所以也就是说每张图片都可以计算出一个H矩阵。

二内参限制我们把H矩阵(3*3)写成3个列向量形式，那么我们把H矩阵又可写成：（3)注:lamda是个放缩因子标量，也是s的倒数。

那么现在我们要用一个关键性的条件：r1和r2标准正交。

正交：（4）单位向量（模相等）：（5）这个两个等式是非常优美的，因为它完美的与绝对二次曲线理论联系起来了，这里就不展开了。

三相机内参的求解我们令：（6）我们可知B矩阵是个对称矩阵，所以可以写成一个6维向量形式：（7）我们把H矩阵的列向量形式为：（8）那么根据等式（8）我们把等式（4）改写成：（9）（10）最后根据内参数限制条件(等式（4）(5)):(11)即,(12)V矩阵是2*6矩阵，也就是说每张照片可建立起两个方程组，6个未知数。

相机标定一、实验原理相机标定就是求解相机的内参数以及畸变参数的过程。

相机的标定主要有两种：传统的摄像头标定方法和摄像头自标定方法，典型的有：（1）Tsai（传统的标定方法）；（2）张正友（介于传统和自标定之间）。

1999年，微软研究院的张正友提出了基于移动平面模板的相机标定方法。

此方法是介于传统标定方法和自标定方法之间的一种方法，传统标定方法虽然精度高设备有较高的要求，其操作过程也比较繁琐，自标定方法的精度不高，张正友标定算法克服了这两者的缺点同时又兼备二者的优点，因此对办公、家庭的场合使用的桌面视觉系统(DVS)很适合。

张正友标定方法由于简单、效果好而得到广泛使用。

张正友标定法的标定步骤：1、打印一张模板并贴在一个平面上；2、从不同角度拍摄若干张模板图像；3、检测出图像中的特征点；4、求出摄像机的外参数（单应性矩阵）和内参数（最大似然估计）；5、求出畸变系数；6、优化求精。

张正友标定方法的主要思想是：1、相机内参矩阵其中，q 的坐标系是默认的OpenCV 的像素坐标系，Q 的坐标系是标定板坐标系，Z 轴为0，原点在标定板的某个内角点上（标定板上角点的坐标均为[*,*,0]的形式），在OpenCV 3.0中使用的是（[i ∗Squres_Size ，j ∗Square_Size ，0]的形式）。

其中fx 和fy 表示相机x 轴和y 轴的焦距，s 表示成像平面x 轴和y 轴的不正交性。

2、基础公式对于不同位置的棋盘格到相机的成像，可以使用下面的公式进行表示：其中，[R|t]表示棋盘格坐标系相对于相机坐标系的位姿。

把矩阵R 和M ~写开，如下式所示：进行化简得：其中[u v 1]是已知量，[X Y 1]也是已知量，A 和[r1 r2 t]是未知量。

其中H=A[r1 r2 t]又叫做单应性矩阵，可以使用下面的3中所述的方法求解。

3、单应矩阵求解这里使用的方法基于最大似然准则：假设提取的m 存在均值为0，噪声协方差矩阵为的高斯白噪声。

文章编号:1006 2610(2022)06 0079 06顾及物方坐标误差的张正友标定算法瞿　锴1,郭中佑2,唐育林2,陈抒录3,陈朋明2(1.中机国际工程设计研究院有限责任公司,长沙　410021;2.湘潭大学土木工程学院,湖南湘潭　411100;3.湘潭市勘测设计院,湖南湘潭　411100)摘　要:针对传统张正友相机标定算法未顾及系数矩阵偶然误差的问题,引入以变量误差模型为基础的整体最小二乘平差方法以提高其解算结果精度,同时开发相机参数计算软件以提高解算效率㊂先以尼康D800为例,使用控制变量法分析相片数量和拍摄角度对参数解算精度的影响,确定相片数量和拍摄方式后,将改进前后的张正友标定算法在尼康D800和索尼ILCE-7R 两种常用非量测型相机上开展实验㊂结果表明:随着相片数量的增加,解算精度不会一直提高,而是逐渐趋于平稳;不同角度拍摄的相片参与解算比单一角度的精度更高;与传统算法相比,采用整体最小二乘得到的参数精度更好㊂可为相关生产单位的相机检校工作提供有用的价值参考㊂关键词:相机检校;张正友标定算法;整体最小二乘;单应性矩阵中图分类号:P237 文献标识码:A DOI :10.3969/j.issn.1006-2610.2022.06.013Zhang Zhengyou 's Calibration Algorithm with consideration of Object -squared Coordinate ErrorsQU Kai 1,GUO Zhongyou 2,TANG Yulin 2,CHEN Shulu 3,CHEN Pengming 2(1.China Machinery International Engineering Design &Research Institute Co.,Ltd ,Changsha　410021,China ;2.School of Civil Engineering and Mechanics ,Xiangtan University ,Xiangtan　411100,China ;3.Xiangtan Institute of Surveying Designing ,Xiangtan 411100,China )Abstract :Aiming at the problem that the traditional Zhang Zhengyou camera calibration algorithm does not take into account the accidental error of thecoefficient matrix ,the overall least squares leveling method based on the variable error model is introduced to improve the accuracy of the solution re⁃sults ,and the camera parameter calculation software is developed to improve the solution efficiency.Taking the Nikon D800as an example ,the influ⁃ence of the number of photos and shooting angle on the parameter solving accuracy is analyzed by using the control variable method ,and after determi⁃ning the number of photos and shooting method ,the Zhang Zhengyou calibration algorithm before and after the improvement is applied on the Nikon D800and Sony ILCE-7R ,two commonly used non-measurement cameras.The results show that with the increase of the number of photos ,the solution accuracy does not continuously increase ,but gradually level off.Photos taken from different angles participate in the solution with higher accuracy than those taken from a single angle ;Compared with the traditional algorithm ,the parameter accuracy obtained by using the overall least squares is better.The study results can provide a useful value reference for the camera calibration work of the relevant production units.Key words :camera correction ;Zhang Zhengyou calibration algorithm ;overall least squares ;homography matrix 收稿日期:2022-07-14 作者简介:瞿锴(1996-),男,湖南省常德市人,助理工程师,主要从事工程测绘工作. 基金项目:湖南省自然科学基金项目(14JJ7039).0　前　言近年来,无人机摄影测量技术得到了迅猛发展,在生产实践中广泛应用于农业㊁林业㊁海洋等领域[1]㊂在无人机摄影测量中普遍采用非量测型数码相机,但其存在着相机光学畸变大的缺点,故需要对相机成像参数予以检校㊂相机检校是通过对已知坐标或相对位置关系明确的物体进行拍照,利用物体相片量测坐标与已知坐标之间的对应关系来求取相机成像参数㊂为保证参数解算的可靠性,通常需要进行大量的多余观测,这些数据被认为服从高斯-马尔可夫模型(Gauss -Markov Theory),可构建众多的观测数据方程,在最小二乘准则约束下求解方程获得参数的最优估值㊂相机检校常用的Tsai 算法㊁张正友标定算法㊁灭点法和直接线性变换(Direct Line transform,DLT)算法在计算相机参数时都采用了G-M 模型[2]㊂但G-M 模型只考虑了观测向量包含误差,则将导致应用传统最小二乘准则估计的参数不是最优解㊂97西北水电㊃2022年㊃第6期===============================================考虑到G-M模型的缺点,整体最小二乘准则被提出并逐渐应用到上所述相机检校工作中所列方程求解的约束条件[3]㊂如陶叶青[4]㊁杨会军[5]㊁孔建[6]㊁鲁铁定[7]和孙培芪等[8]将整体最小二乘应用到解决方案中去,实验结果表明,整体最小二乘在拟合函数㊁转换参数的求解问题中更加具备优势㊂其次是检校场类型对相机检校起着很重要的作用,倪静㊁曹良忠等[9-10]研究了相对条件较为苛刻的野外检校场㊁田雷等[11]应用了室内检校场检校㊁李平[12]简化了三维空间检校场,应用二维检校场进行相机检校中比较出名的是张正友检校法[13],这种方法只需要一个规则的格网状平板,通过移动相机从不同位置对其拍照,即可计算出内外方位元素㊂考虑到野外检校场和室内三维检校场建立与维护的不易,而无人机摄影测量相机又需要频繁检校,本文采用自制的格网型模板,先以尼康D800为例,利用控制变量法分析相片数量和拍摄角度对参数解算精度的影响,确定相片数量和拍摄方式后,使用Nikon D800相机与索尼ILCE-7R相机拍摄模板获取数据㊂接着,考虑到整体最小二乘算法较最小二乘算法的模型更为完善,顾及到了模型中系数矩阵可能存在的误差,使用整体最小二乘算法对张正友检校算法予以改进,并通过试验进行验证,以提高非量测相机的检校精度,保证测绘生产获得较高质量的成果㊂1　研究方法自检校方法是Maybank[14]提出的一种与相机移动和空间几何信息无关的相机检校方法㊂这方面国内的工作有孟晓桥等[15]从多个方向拍摄模板图像,应用圆检校算法,利用射影的不变性,解算出了相机参数㊂从计算原理讲,张正友检校法也可以归于自检校方法㊂相机的自检校方法是指不需要依据一个已经建立好的检校场,是通过影像与影像之间的对应关系,计算出相机的参数㊂此种方法的优点是价格低廉,可以随时随地对相机进行检校;缺点则是精度会比较低㊂张正友相机标定算法是将打印的棋盘格贴在平面上,从不同角度进行拍摄,利用像点和物点的几何关系确定相机内外方位元素的初值,然后将其作为非线性相机几何模型的初值,求出单应性矩阵,以此解算出摄像机内外参数和畸变系数,模型如图1所示㊂图1　张正友算法模型1.1　单应性矩阵设世界坐标某点的物理坐标为M=[U,V,W, 1]T,其对应的像素坐标系下的像素坐标点为m= [u,v,1]T,故标定用的棋盘格平面到图像平面的单应性关系为:sm=A[R　t]M,A=α　γ　u00　β　v0éëêêêêùûúúúú0　0　1(1)其中:s为尺度因子;A为相机内参矩阵;包括仿射变换和透视投影;(u0,v0)为像主点坐标;α,β为焦距与像素横纵比的融合;γ为径向畸变参数;R为旋转矩阵;t为平移向量;[R t]也被成为相机的外参矩阵㊂设棋盘格平面W=0,可得:suvéëêêêùûúúú1=A[r1,r2,r3,t]UVéëêêêêêùûúúúúú1=A[r1,r2,t]UVéëêêêùûúúú1(2) A[r1,r2,t]为单应性矩阵,即内参矩阵和外参矩阵的积,记作H=[r1,r2,t],得:uvéëêêêùûúúú1=1sH11　H12　H13H21　H22　H23H31　H32　Héëêêêêùûúúúú33UVéëêêêùûúúú1(3) 所以每一格标定板的角点可以提供两个坐标约束方程:u=H11U+H12V+H13H31U+H32V+H33v=H21U+H22V+H23H31U+H32V+Hìîíïïïï33(4)08瞿锴,郭中佑,唐育林,陈抒录,陈朋明.顾及物方坐标误差的张正友标定算法=============================================== 上式对于同一张图片上所有的角点均成立,其中[u ,v ]是像素坐标系下标定板的角点坐标,[U ,V ]是世界坐标系下标定板的角点坐标㊂1.2　内参矩阵A 和外参矩阵的计算旋转矩阵R 存在单位正交的关系,即r T1r 2=0㊁r T 1r 2=r T 2r 2=1,令A -T A -1=B ,B 为对称阵,可得:H T1B H 2=0H T 1B H 1=H T 2B H 2{=1(5) 化简得:v T 11v T 11-v T éëêêùûúú22b =vb =0(6)公式(6)中:v 是仅包含H 元素的矩阵,b =[B 11,B 12,B 22,B 23,B 33]T ㊂由于矩阵H 已知,矩阵v 又全部由矩阵H 的元素构成,因此矩阵v 已知㊂每张标定板图片可以提供一个约束关系,采用最小二乘拟合最佳的向量b ,得到矩阵B ,所以有:v 0=B 12B 13-B 11B 23B 11B 22-B 212α=1B 11,β=B 11B 11B 22-B 212γ=-B 12α2β,u 0=γv 0β-B 13α2(7) 相机的内参矩阵A 取决于相机的内部参数,无论标定板和相机的位置关系怎么变换,相机的内参矩阵是不变的㊂而外参矩阵反映的是标定板和相机的位置关系,对于被拍摄的不同图片,标定板和相机的位置关系已经改变,此时每一张图片对应的外参矩阵都是不相等的㊂因为H =A [r 1,r 2,t ],求解得到H 和A 后,可通过式(8)得到外参矩阵㊂[r 1,r 2,t ]=A -1H(8)1.3　畸变系数k 的计算对于畸变系数k ,张氏标定法只关注了影响最大的径向畸变,数学表达式为[18]:(u -u 0)(x 2+y 2)(u -u 0)(x 2+y 2)2(v -v 0)(x 2+y 2)(v -v 0)(x 2+y 2)éëêêùûúú2k 1k éëêêùûúú2^u -u ^éëêêùûúúv -v (9)公式(9)中:(u ,v )是理想无畸变的像素坐标;(^u ,^v )是实际畸变后的像素坐标;(u 0,v 0)代表主点;(x ,y )是理想无畸变的连续图像坐标;(^x,^y )是实际畸变后的连续图像坐标;k 1和k 2为前两阶的畸变参数,记作Dk =d ,可得:k =[k 1,k 2]T =(D T D )-1D T d(10)1.4　整体最小二乘原理在实际工作中,系数矩阵中也会包含误差[7-8],故真误差的数学期望往往不为0㊂求解Ax =b 时,应该同时考虑C 和c 二者的误差和扰动,令B 矩阵的误差扰动为E ,向量b 的误差向量为e ,矩阵方程为:(C +E )x =c +e (11) 可得:c +e -(x T ×I )E -Cx =0(12)式中:x T ×I 表示两个矩阵做克罗内克积,令F =[i ,-(x T ×I )],v =[e ,E ]T ,可得:Fv +c -Cx =0(13) 整体最小二乘准则为:min φ=v T Pv -2λT (Fv +c -Cx )(14)式中:P 为权矩阵;λ为Largrange 乘数,分别进行偏导:əφəv =v T P -λT F =0əφəx =(C +I )Tλ=0əφəλ=Fv +c -Cx =ìîíïïïïïï0(15) 可得:x =[(C +I )T (FP -1F )-1C ]-1(C +I )T (FP -1F )-1c(16) 对公式(6)和公式(10)采用整体最小二乘计算,得到最优估计值㊂2　检校场与计算软件的设计2.1　室内二维检校场的实现这里选用尼康D800和索尼ILCE-7R 两种相机获取影像数据,拍摄时设置为手动曝光方式,物镜对焦于无穷远,相机传感器参数见表1㊂表1　尼康D800和索尼ILCE-7R 相机传感器参数相机参数类型传感器尺寸/mm 分辨率/像素像素大小/μm 焦距/mm 尼康D80035.9×247360×4912 4.8835索尼ILCE-7R35.9×244800×32004.8835计算机型号为联想Lenovo G50-7020351笔记本电脑,CPU 酷睿i5-4258U,内存DDR3L 1600MHz㊂实验选取棋盘格作为相机检校的模板,二维控制场采用基于纯平液晶显示器的检校方法[19],使18西北水电㊃2022年㊃第6期===============================================所有控制点在一个平面上㊂为研究拍摄角度对标定精度的影响,从不同角度进行拍摄,得到图像,棋盘格大小为9×5个方格,每个方格尺寸为45mm ×45mm(长×宽),角点数量为32个,部分影像如图2所示㊂图2　不同方向拍摄的棋盘格控制场2.2　相机检校软件的设计这里根据计算原理利用Visual Studio 平台和C #设计相机检校软件,主要分为如下几项功能:文件读取数据㊁最小二乘约束㊁整体最小二乘约束㊁畸变改正㊁查看报告和退出等功能㊂程序操作非常简便,首先点击可执行文件(.exe),启动相机检校程序,在相同目录下添加不同角度拍摄的棋盘格影像,可以为tif㊁jpg㊁bmp 等格式,使用对应算法自动提取棋盘格角点的坐标作为后续计算的数据,程序解算会自动将最小二乘约束与整体最小二乘约束解算所得到的参数值作为初始值使用,直接点击计算输出按钮可得到最终的相机参数报告,包括最小二乘约束下和整体最小二乘约束下的不同相机检校算法的相机参数㊂参数解算部分结果如图3所示㊂前面部分为相机的检校参数值,后面部分为各参数的解算误差㊂图3　相机检校软件3　实验结果与分析3.1　相片数量对检校的影响为探究相片数量对相机检校精度的影响,依次使用3㊁5㊁7㊁9㊁11不同数量的相片参与计算,结果如表2所示,参与解算相片数量的不同会导致检校参数的精度有所差异㊂相机参数包括:f (焦距),x 0㊁y 0(主点坐标)和k 1㊁k 2(畸变参数)㊂y 值在3张相片解算时最优,后面呈波动分布㊂f 值在9张相片解算时达到最优,11张相片为35.85445,相比而言,值下降了0.066㊂x 0㊁k 1和k 2的值随着参加解算相片的数量越多,精度越高,分别为3661.5㊁-0.0114和0.0403㊂3张相片进行解算时,f ㊁k 1和k 2与参考值相差最大,这是因为相片的数量太少导致不能够充分反应出镜头的畸变情况㊂随着后面参与解算的相片数量越来越多,值变化的越来越趋于平稳状态,说明相片并不是越多越好,但也不能太少,过少会使得结果精度不高,过多则会导致计算量的增加,从而导致效率降低㊂表2　尼康D800相片数量与检校参数计算结果　f /mmx 0/像素y 0/像素k 1k 23张相片5张相片7张相片9张相片11张相片参考值36.425743236.230751936.269847936.0533535.8544535.986923631.93631.03679.53642.93661.53654.0312474.42416.02455.52412.32403.22486.250-0.0440-0.0315-0.0247-0.0142-0.01146.72960E-050.14330.14290.10800.05060.0403-5.6425E-0528瞿锴,郭中佑,唐育林,陈抒录,陈朋明.顾及物方坐标误差的张正友标定算法===============================================3.2　拍摄角度对检校的影响为探究拍摄相片的角度对检校结果的影响,采用两种方式拍摄相片,先拍摄一组近似平行的照片计算检校参数,再以不同的角度拍摄相片计算检校参数,分析得到的检校结果,两组相片的检校结果如表3所示㊂当相机位置处于平行状态时,f和x0的精度较高,分别为35.96734和3654.2;相片处于不同的角度时,y0㊁k1和k2精度较高,依次为2404.7㊁-0.006㊁0.019㊂一般认为,近似平行的相片计算的检校参数精度会低于不同角度拍摄的照片计算的检校参数精度㊂因为张正友检校方法是通过不同相片之间的关系计算的,所以拍摄的相片近似平行的时候,有的参数之间会存在很强的相关性,就会影响最后的结果㊂在使用张正友标定法对相机进行检校时,应该以不同的角度对棋盘格进行拍摄,从而保证高精度的校验结果㊂y0㊁k1和k2的精度与预期表现相符,而f和x0的结果精度恰恰相反,可能与相片拍摄的质量有关㊂表3　尼康D800相机拍摄角度与检校参数计算结果　f/mm x0/像素y0/像素k1k2近似于平行处于不同角度参考值35.9673435.4957535.986923654.23624.73654.0312404.62404.72486.250-0.0282-0.00606.72960E-050.15710.0190-5.6425E-053.3　改进前后方法的对比由以上可知,相机拍摄相片时应尽量从不同方向进行拍摄,数量根据区域大小确定,本次使用11张不同方向拍摄的相片参与解算㊂为检验改进前后张正友方法在尼康D800和索尼ILCE-7R两种相机中的情况,分别使用两款相机随机在5个摄站对棋盘格进行拍摄,每张相片不重复,角点数量为45个,部分图片如图2所示㊂分别采用经典最小二乘算法与整体最小二乘算法求解单应性矩阵参数,进而解算相机参数,并将结果进行对比分析㊂表4为尼康D800和索尼ILCE-7R相机内定向参数计算结果,在尼康D800数据集上,LS和TLS法得到的焦距f 与参考值分别相差0.421mm㊁0.259mm;主点坐标(x0,y0)分别相差176.117像素㊁44.159像素;畸变参数k1/k2分别相差0.037/0.012㊁0.466/0.150;计算时间相差0.345㊂在索尼ILCE-7R数据集上,LS 和TLS法得到的焦距f与参考值分别相差0.435 mm㊁0.120mm;主点坐标(x0,y0)分别相差155.731像素㊁93.362像素;k1/k2分别相差0.1225/0.110㊁0.150/0.128;计算时间相差0.185㊂相比LS,TLS 计算得到的参数精度均得到了不同程度的提升,说明TLS确实能够较好的考虑观测值和系数存在的误差,解算出更为精确的结果㊂虽然相机参数精度均有提升,但增加了计算时间,这是必然的结果,因为计算量的增加,势必会导致计算效率的降低㊂总的来说,整体最小二乘方法相比传统的最小二乘方法在张正友标定算法相机参数解算中表现更优㊂表4尼康D800和索尼ILCE-7R相机内定向参数计算结果参数尼康D800LS法TLS法参考值索尼ILCE-7RLS法TLS法参考值f/mm36.4082035.7278835.9869235.7278836.0430936.1626 x0/像素3742.33612.23654.0313556.13612.13702.249 y0/像素2319.72472.12486.2502472.12442.62418.318 k1-0.0372-0.01236.72960E-050.12250.1099-3.99663E-05 k20.4657-0.1497-5.6425E-05-0.1497-0.12781.54691E-07计算时间/s0.3727460.717842　0.2951500.480264　4　结　论本文在传统张正友的算法上进行改进,考虑顾及物方控制点坐标误差,用整体最小二乘约束求解单应性矩阵,将其精度得到提高的结果参与后面相机参数的运算㊂为确定最佳解算相片数据,先以尼康D800为例,使用控制变量法探究拍摄角度和相片数量对解算结果的影像㊂然后采用改进前后的张正友标定算法在尼康D800和索尼ILCE-7R相机开展实验,结果表明:(1)解算参数x0㊁k1和k2的值随着参加解算相片的数量越多,精度越高,值逐渐趋于平稳,说明盲目增加相片数量不会一直提高解算精度,反而会增大计算量㊂(2)以不同角度拍摄的照片计算的检校参数y0㊁k1和k2精度高于近似平行的相片计算的检校参数精度,不同角度的相片能够更好提供约束条件,得到更优解㊂(3)在尼康D800和索尼ILCE-7R数据上,相比传统张正友算法,f(焦距)㊁x0㊁y0(主点坐标)和k1㊁k2(畸变参数)的解算精度均得到了提升,但计算量提升了一倍,这在预期之内㊂本次实验的某些相机参数精度与预期不相符,可能是因为相片质量没有做抗差估计,在进行解算时无法保证每张相片的质量统一,从而导致结果出现差异性㊂后面可以针对不同相机检校算法验证本38西北水电㊃2022年㊃第6期===============================================文所提方法的可行性,并为了保证数据质量统一标准,对拍摄相片进行抗差处理㊂参考文献:[1]　王柯,付怡然,彭向阳,等.无人机低空遥感技术进展及典型行业应用综述[J].测绘通报,2017(S1):79-83.[2]　Cheng C L,SchneeweißH.On the Polynomial Measurement Error Model[J].Springer Netherlands,2002.[3]　汪奇生.线性回归模型的总体最小二乘平差算法及其应用研究[D].昆明:昆明理工大学,2014.[4]　陶叶青.总体最小二乘模型及其在矿区测量数据处理中的应用研究[D].徐州:中国矿业大学,2015.[5]　Huijun Y.GNSS Height Fitting Based on Total Least Squares[J].GNSS World of China,2015,40(06):99-101.[6]　Ian K J,Yibin Y,Han W.Iterative Method for Total Least -Squares[J].Geomatics and Informa tion Science of Wuhan University,2010,35(06):711-714.[7]　Tieding L,Ian Z S.An Iterative Algorithm for Total Least SquaresEstimation[J].Geomatics and Informa tion Science of Wuhan Uni⁃ver sity,2010,35(11):1351-1354.[8]　Sun P Q,Bu J Z,Tao T Y,et al.Improvement of camera calibra⁃tion DLT algorithm[J].Journal of Hefei University of Technology(Natural Science),2020,43(09):1229-1233,1280.[9]　倪静.非量测相机检校方法与系统实现[D].东营:中国石油大学(华东),2016.[10]　Liang-Zhong C,Liao Y,Pei-Tao K,et al.Calibration of non-metric digital camera based on calibration field [J].Science of Surveying and Mapping,2015,40(02):132-137.[11]　田雷,马然.无人机非量测相机检校方法研究[J].测绘通报,2016(07):81-83,89.[12]　李平.非量测型数码相机的检校[D].西安:西安科技大学,2009.[13]　Zhang,Zhengyou.A Flexible New Technique for Camera Cali⁃bration[J].IEEE Transactions on Pattern Analysis &MachineIntelligence,2000,22(11):1330-1334.[14]　Maybank S J,Faugeras O D.A theory of self -calibration of amoving camera [J].International Journal of Computer Vision,1992,8(02):123-151.[15]　孟晓桥,胡占义.摄像机自标定方法的研究与进展[J].自动化学报,2003,29(01):110-124.[16]　徐刚.基于灭点几何和径向准直约束的数码相机标定研究[D].桂林:桂林电子科技大学,2011.[17]　张伟航,杨姗姗,李建刚.近景摄影测量相机检校模型综述[J].山西建筑,2013,39(20):207-208.[18]　吴飞.总体最小二乘拟合推估方法及其应用研究[D].南昌:东华理工大学,2015.[19]　Jingnan L,Wenxian Z,Peiliang X.Overview of total leastsquares methods[J].Geomatics and information science of wu⁃han university,2013,38(05): 505-512. (上接第70页)筛选部分物理指标作为主控因子,将力学参数作为目标参数,采用数理统计手段可拟合出合适的回归方程,用以推测选取滑带土的力学参数㊂(2)以白龙江水泊峡水电站Ⅶ号滑坡为研究对象,在获取大量滑带土试验数据的基础上,将物理指标ω㊁ρs 及e 作为自变量,采用多元线性回归拟合内摩擦角φ的预测方程;分别采用多元线性回归和非线性回归建立黏聚力c 的预测模型,经对比分析,非线性拟合的计算结果更贴合实测数据㊂(3)沿Ⅶ号滑坡主滑方向剖切计算剖面,根据滑带土的物理指标与预测模型计算力学参数,计算得Ⅰ区滑带土c 值为6.9kPa,φ值为21.3°;Ⅱ区滑带土c 值为28.6kPa,φ值为28.9°㊂稳定性计算结果表明:天然工况下K =1.434,处于稳定状态;地震工况下,K =1.063,处于整体变形~滑动状态㊂计算结果基本与实际情形吻合㊂参考文献:[1]　吴博,赵法锁,王雷,等.黄土-基岩滑坡滑带土特性及其演化过程[J].煤田地质与勘探,2018,46(04):129-134,141.[2]　游敏,聂德新.利用大型岩质滑坡形成的环境条件及重力圧密原理分析滑面形态与强度参数[J].工程地质学报,2010,18(04):470-476.[3]　孙淼军,唐辉明,王潇弘,等.蠕动型滑坡滑带土蠕变特性研究[J].岩土力学,2017,38(02):385-391,399.[4]　闫玉平,肖世国.考虑滑带强度参数分区取值的堆积层滑坡稳定性分析方法[J].中国地质灾害与防治学报,2020,31(02):44-49.[5]　郑莉,王卫.联合水库右岸滑坡滑带土强度参数分析[J].三峡大学学报(自然科学版),2019,41(S1):134-138.[6]　赵宝平,罗先启,毕强,等.库水作用下滑带土力学参数降低的模拟[J].西北水电,2009(01):66-69.[7]　田桂莉.金沙江上游滑带土强度特性的环剪试验研究[D].吉林:吉林大学,2018.[8]　李鹏.Hoek-Brown 准则在软岩~中硬岩坝基岩体力学指标参数选取中的应用 以陕西省延安市龙安水利枢纽工程为例[J].西北水电,2019(04):16-20.[9]　杨期祥.大型深层滑坡滑带土力学强度参数综合选取研究[J].铁道勘察,2022,48(04):61-66.[10]　颜清,彭小平.工程实验数据的非线性拟合方法[J].计算机与应用化学,2015,32(03):365-368.48瞿锴,郭中佑,唐育林,陈抒录,陈朋明.顾及物方坐标误差的张正友标定算法===============================================。

张氏标定法张正友针对径向畸变问题提出了一个新的求解摄像机内外参数的方法，即张氏标定法，该方法是介于传统标定和自标定之间的一种方法，它只需要摄像机对某个标定板从不同方向拍摄多幅图片，通过标定板上每个特征点和其像平面的像点间的对应关系，即每一幅图像的单应矩阵来进行摄像机的标定，由于该方法模板制作容易，使用方便，成本低，鲁棒性好，准确率高，因此得到了较为广泛的应用【49】。

该算法也属于两步法，摄像机和模板可以自由的移动，不需要知道运动参数。

本文即采用张氏摄像机标定的方法。

（1） 单应矩阵H 的计算由前面所述内容可知，根据摄像机针孔成像模型，可以得到下面的等式：[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡11][1321W W W W W W Z Y X t r r r A Z Y X t R A v u s 现在世界坐标系平面置于标定模板所在的平面，即0=W Z 。

则上式可变为如下形式：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡110][121W W W W Y X t r r A Y X t R A v u s 其中，i r 表示旋转矩阵R 的第i 列向量。

令[]T Y X M 1~=，[]T v u m1~=，则上式可简写为：M H m s ~~=其中：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡===1][323123222113121132121h h h h h h h h h h h t r r A H H 即为单应性矩阵。

有M H ms ~~=式子可以推出： ⎪⎩⎪⎨⎧++=++=++=13231232221131211h X h s h Y h X h sv h Y h X h su 从而推得：⎩⎨⎧++=++++=++23222132311312113231h Y h X h v vYh vXh h Y h X h u uYh uXh 令：[]3231232221131211'h h h h h h h h h =则：0'10000001=⎥⎦⎤⎢⎣⎡------h v vYvXYXu uY uX Y X 上式可以看作0'=Sh ，那么矩阵S S T最小特征值所对应的特征向量就是该方程的最小二乘解。

相机标定一、实验原理相机标定就是求解相机的内参数以及畸变参数的过程。

相机的标定主要有两种：传统的摄像头标定方法和摄像头自标定方法，典型的有：（1）Tsai（传统的标定方法）；（2）张正友（介于传统和自标定之间）。

1999年，微软研究院的张正友提出了基于移动平面模板的相机标定方法。

此方法是介于传统标定方法和自标定方法之间的一种方法，传统标定方法虽然精度高设备有较高的要求，其操作过程也比较繁琐，自标定方法的精度不高，张正友标定算法克服了这两者的缺点同时又兼备二者的优点，因此对办公、家庭的场合使用的桌面视觉系统(DVS)很适合。

张正友标定方法由于简单、效果好而得到广泛使用。

张正友标定法的标定步骤：1、打印一张模板并贴在一个平面上；2、从不同角度拍摄若干张模板图像；3、检测出图像中的特征点；4、求出摄像机的外参数（单应性矩阵）和内参数（最大似然估计）；5、求出畸变系数；6、优化求精。

张正友标定方法的主要思想是：1、相机内参矩阵其中，q 的坐标系是默认的OpenCV 的像素坐标系，Q 的坐标系是标定板坐标系，Z 轴为0，原点在标定板的某个内角点上（标定板上角点的坐标均为[*,*,0]的形式），在OpenCV 3.0中使用的是（[i ∗Squres_Size ，j ∗Square_Size ，0]的形式）。

其中fx 和fy 表示相机x 轴和y 轴的焦距，s 表示成像平面x 轴和y 轴的不正交性。

2、基础公式对于不同位置的棋盘格到相机的成像，可以使用下面的公式进行表示：其中，[R|t]表示棋盘格坐标系相对于相机坐标系的位姿。

把矩阵R 和M ~写开，如下式所示：进行化简得：其中[u v 1]是已知量，[X Y 1]也是已知量，A 和[r1 r2 t]是未知量。

其中H=A[r1 r2 t]又叫做单应性矩阵，可以使用下面的3中所述的方法求解。

3、单应矩阵求解这里使用的方法基于最大似然准则：假设提取的m 存在均值为0，噪声协方差矩阵为的高斯白噪声。

张正友算法原理及其改进由于世界坐标系的位置可以任意选取，我们可以假定世界坐标系和摄像机坐标系重合，故定义模板平面落在世界坐标系的0W Z =平面上。

用i r 表示R 的每一列向量，那么对平面上的每一点，有：[][]123120111W W W W X u X Y s v A r r r t A r r t Y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦(3.1) 这样，在模板平面上的点和它的像点之间建立了一个单应性映射H ，又称单应性矩阵或投影矩阵。

如果已知模板点的空间坐标和图像坐标，那么就已知m 和M ，可以求解单应性矩阵H 。

)1,,(w w Y X )1,,(v u因为11W W u X s v H Y ⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦,其中11121321222331321h h h H h h h h h ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,可推出： 11121321222331321W W W W W W su h X h Y h sv h X h Y h s h X h Y =++⎧⎪=++⎨⎪=++⎩(3.2) 故，1112133132212223313211W W W W WW W W h X h Y h u h X h Y h X h Y h v h X h Y ++⎧=⎪++⎪⎨++⎪=⎪++⎩(3.3) 将分母乘到等式左边，即有31321112133132212223W W W W W W W W uX h uY h u h X h Y h vX h vY h v h X h Y h ++=++⎧⎨++=++⎩ (3.4)又令[]Th h h h h h h h h 3231232221131211='，则1000001W W W W WWWW X Y uX uY u h X Y vX vY v --⎡⎤⎡⎤'=⎢⎥⎢⎥--⎣⎦⎣⎦(3.5) 多个对应点的方程叠加起来可以看成Sh d '=。

张正友算法原理及其改进由于世界坐标系的位置可以任意选取，我们可以假定世界坐标系和摄像机坐标系重合，故定义模板平面落在世界坐标系的0W Z =平面上。

用i r 表示R 的每一列向量，那么对平面上的每一点，有：[][]123120111W W W W X u X Y s v A r r r t A r r t Y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦(3.1) 这样，在模板平面上的点和它的像点之间建立了一个单应性映射H ，又称单应性矩阵或投影矩阵。

如果已知模板点的空间坐标和图像坐标，那么就已知m 和M ，可以求解单应性矩阵H 。

)1,,(w w Y X )1,,(v u因为11W W u X s v H Y ⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦,其中11121321222331321h h h H h h h h h ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,可推出： 11121321222331321W W W W W W su h X h Y h sv h X h Y h s h X h Y =++⎧⎪=++⎨⎪=++⎩(3.2) 故，1112133132212223313211W W W W W W W W h X h Y h u h X h Y h X h Y h v h X h Y ++⎧=⎪++⎪⎨++⎪=⎪++⎩(3.3) 将分母乘到等式左边，即有31321112133132212223W W W W W W W W uX h uY h u h X h Y h vX h vY h v h X h Y h ++=++⎧⎨++=++⎩ (3.4)又令[]Th h h h h h h h h 3231232221131211='，则1000001W W W W WWWW X Y uX uY u h X Y vX vY v --⎡⎤⎡⎤'=⎢⎥⎢⎥--⎣⎦⎣⎦(3.5) 多个对应点的方程叠加起来可以看成Sh d '=。