立体视觉入门指南（3）：相机标定之张式标定法【超详细值得收藏】

亲爱的同学们，我们的世界是3D世界，我们的双眼能够观测三维信息，帮助我们感知距离，导航避障，从而翱翔于天地之间。而当今世界是智能化的世界，我们的科学家们探索各种机器智能技术，让机器能够拥有人类的三维感知能力，并希望在速度和精度上超越人类，比如自动驾驶导航中的定位导航，无人机的自动避障，测量仪中的三维扫描等，都是高智机器智能技术在3D视觉上的具体实现。

立体视觉是三维重建领域的重要方向，它模拟人眼结构用双相机模拟双目，以透视投影、三角测量为基础，通过逻辑复杂的同名点搜索算法，恢复场景中的三维信息。它的应用十分之广泛，自动驾驶、导航避障、文物重建、人脸识别等诸多高科技应用都有它关键的身影。

本课程将带大家由浅入深的了解立体视觉的理论与实践知识。我们会从坐标系讲到相机标定，从被动式立体讲到主动式立体，甚至可能从深度恢复讲到网格构建与处理，感兴趣的同学们，来和我一起探索立体视觉的魅力吧！

本课程是电子资源，所以行文并不会有太多条条框框的约束，但会以逻辑清晰、浅显易懂为目标，水平有限，若有不足之处，还请不吝赐教！
个人微信：EthanYs6，加我申请进技术交流群 StereoV3D，一起技术畅聊。
CSDN搜索：Ethan Li 李迎松，查看网页版课程。
随课代码，将上传至github上，地址：StereoV3DCode：https://github.com/ethan-li-coding/StereoV3DCode

同学们好，在上一篇立体视觉入门指南（2）：关键矩阵中，我们对立体视觉的3个关键矩阵：本质矩阵 E E E、基础矩阵 F F F、单应性矩阵 H H H作了较为详细的描述，同时给出了本质矩阵、单应性矩阵的求解方法以及本质矩阵分解外参 R ， t R，t R，t的具体公式。更加难能可贵的是，我们在博文最后提供了几个作业题并在Github开源了参考答案代码【我知道很多心理都在默念李博666】【当然肯定也有一些人在默念这太easy了李博能不能上点难度】。无论如何，博主觉得这是一件有意义的事情，只希望没有误人子弟。

而本篇的内容，则是立体视觉的绝对核心模块：相机标定。虽然它因为技术相对成熟，如今研究的人不多，也容易被人忽略，往往用一个开源算法库如Opencv或者Matlab标定工具箱就直接搞定，但实际在立体视觉工程化、产品化时，开源工具由于其精度不高、灵活度低而不建议直接使用，企业往往是自己开发相机标定算法。相机标定作为立体视觉的核心模块，掌握其理论是相当必要的，对我们深入理解立体视觉技术大有帮助，有开山辟路之效。

相机标定，目的是确定相机的内参矩阵 K K K、外参矩阵 R , t R,t R,t 和 畸变系数 d ( k 1 , k 2 , k 3 , p 1 , p 2 ) d(k_1,k_2,k_3,p_1,p_2) d(k1,k2,k3,p1,p2)。我们暂且先不考虑畸变系数（后面再讨论它），从前篇可知，图像坐标 p p p和世界坐标 P w P_w Pw之间，通过内外参建立投影关系，具体公式如下：
λ p = K [ R t ] [ P w 1 ] (1) \lambda p=K\left[\begin{matrix}R&t\end{matrix}\right] \left[\begin{matrix}P_w\\1\end{matrix}\right] \tag 1 λp=K[Rt][Pw1](1)
所谓标定，即是由大量观测值拟合参数模型的过程，且在此拟合的参数模型是已知的，所以应尽可能探索能便捷获取大量观测值的方案，如果观测值之间还满足一些其他的几何约束就更有助于求解具体单个参数值。

今天所述的Zhang式标定法¹即提供了一种便捷获取大量观测值的的方案，同时观测值之间还满足一类明显的几何约束（即平面约束），可直接求解出内外参。其操作方式非常简单，只需要拍摄带有标定板图案的平面，即可完成相机标定，使标定难度极大降低，如果不追求高精度，打印一张棋盘格标定板图案粘贴到近似平的硬纸板上即可完成标定，加快了立体视觉的入门和普及，影响深远，是相机标定领域绝对的经典。

本篇即带大家深入了解Zhang式相机标定法，掌握本篇对立体视觉的理论掌握及工程实战来说都是非常必要的。

实施方法

Zhang式标定法能够被广泛应用，其中一个重要原因是其实施方法十分简单，不需要专业的工艺制作即可完成。

第一步，设计一张具有明显角点特征，且已知每个角点二维坐标的图案作为标定图案，常见的图案有三种：

棋盘格

ChArUco

圆点

规则的图案设计可以方便的计算出角点在图案内的二维坐标，拿棋盘格来说，角点之间的间隔像素数是固定的，假设左上角角点的坐标为 ( 0 , 0 ) (0,0) (0,0)，则其他角点的像素坐标都可以通过格子的偏移量计算出来，而一张已知DPI的标定板图像，在打印后每个角点的二维空间坐标也是完全已知的（通过像素换算成空间尺寸）。

第二步，将标定板图案以某种方式置于一个平面上。比如最简单的方式是将标定图原尺寸打印出来，然后找一块近似平的平板，将打印后的标定图案粘贴至平板上；更专业高精度的方式是找专业厂家制作高精平板（如陶瓷板）并将标定图案以某种工艺刻印到平板上。这一步的目的是让标定图案的角点都位于一个平面上。

图片来源：https://calib.io/blogs/knowledge-base/calibration-patterns-explained

如此，第一步所描述的二维坐标可以转换成第三维 Z Z Z坐标等于0的三维坐标（将世界坐标系的原点放在标定板的某个角点，Z轴垂直于标定板）。

第三步，移动相机到N（N>=3）个不同的位姿拍摄标定板图案。

第四步，对上一步拍摄的标定板图案进行角点提取，解算标定参数。

matlab软件标定示意图

以上便是相机标定的实施步骤，总结来说，在一个平面标定板上有一组已知空间坐标的角点，相机在多个不同位姿下拍摄角点图案并提取角点的像素坐标，即可完成相机内外参数的解算。这句总结蕴含2个重要的知识点：

标定图案中角点的空间坐标是已知的，且它们都位于一个平面上， Z Z Z坐标等于0
相机需要在多个不同的位姿拍摄角点图案并提取像素坐标

可以看到Zhang式标定法确实是易于实施的方法，其中蕴含的标定理论，是非常有价值的，对于立体视觉的初学者来说，掌握其理论很有必要，还请务必阅读本篇接下来的内容。

理论基础

定义

为了保持和论文的公式一致，我们改变下前两篇博客中的命名习惯，假设图像上某点像素坐标为 m = [ u , v ] T \text m=[u,v]^T m=[u,v]T，空间中某点三维坐标为 M = [ X , Y , Z ] T \text M=[X,Y,Z]^T M=[X,Y,Z]T，两者的齐次表达式分别为 m ~ = [ u , v , 1 ] T , M ~ = [ X , Y , Z , 1 ] T \widetilde{\text m}=[u,v,1]^T,\widetilde{\text M}=[X,Y,Z,1]^T m =[u,v,1]T,M =[X,Y,Z,1]T。则文章片头的投影公式可表达为：
s m ~ = A [ R t ] M ~ (2) s\widetilde{\text m}=\text A\left[\begin{matrix}\text R&\text t\end{matrix}\right]\widetilde{\text M} \tag 2 sm =A[Rt]M (2)
其中， s s s为尺度因子， R , t \text R,\text t R,t为外参矩阵，用于将世界坐标系坐标转换成相机坐标系坐标，前者为旋转矩阵，后者为平移矢量。 A \text A A为相机的内参矩阵，具体的，
A = [ α γ u 0 0 β v 0 0 0 1 ] A=\left[\begin{matrix}\alpha&\gamma&u_0\\0&\beta&v_0\\0&0&1\end{matrix}\right] A=⎣⎡α00γβ0u0v01⎦⎤
其中， ( u 0 , v 0 ) (u_0,v_0) (u0,v0)为像主点坐标， α , β \alpha,\beta α,β分别为图像 u u u轴和 v v v轴方向的焦距（像素单位）值， γ \gamma γ为像元轴的倾斜因子。

单应性矩阵

在实施方法中，我们一定能够关注到一个非常关键的信息：标定图案被置于一个平面上。它的目的是为了让标定图案中的角点都位于一个空间平面上，从而当相机拍摄角点成像后，空间平面和像平面之间存在一个单应性变换关系，即可通过一个单应性矩阵 H \text H H将角点的空间坐标转换成图像坐标。

一般情况下，我们将世界坐标系置于标定板某个角点上，并让 Z Z Z轴垂直于标定板，此时角点的 Z Z Z坐标将全部等于0，则三维坐标可表示为 [ X , Y , 0 , 1 ] T [X,Y,0,1]^T [X,Y,0,1]T。同时定义旋转矩阵 R \text R R的第 i i i个列向量为 r i \text r_i ri，则公式 s m ~ = A [ R t ] M ~ s\widetilde{\text m}=\text A\left[\begin{matrix}\text R&\text t\end{matrix}\right]\widetilde{\text M} sm =A[Rt]M 可表示为
s [ u v 1 ] = A [ r 1 r 2 r 3 t ] [ X Y 0 1 ] = A [ r 1 r 2 t ] [ X Y 1 ] (3) s\left[\begin{matrix}u\\v\\1\end{matrix}\right]=\text A\left[\begin{matrix}\text r_1&\text r_2&\text r_3&\text t\end{matrix}\right]\left[\begin{matrix}X\\Y\\0\\1\end{matrix}\right]=\text A\left[\begin{matrix}\text r_1&\text r_2&\text t\end{matrix}\right]\left[\begin{matrix}X\\Y\\1\end{matrix}\right] \tag 3 s⎣⎡uv1⎦⎤=A[r1r2r3t]⎣⎢⎢⎡XY01⎦⎥⎥⎤=A[r1r2t]⎣⎡XY1⎦⎤(3)
如前所述，标定平面和像平面之间可用单应性矩阵 H \text H H进行变换，即有
s m ~ = H M ~ (4) s\widetilde{\text m}=\text H\widetilde{\text M}\tag 4 sm =HM (4)
由于空间点 Z Z Z坐标为0，所以上式中 M ~ = [ X , Y , 1 ] T \widetilde{\text M}=[X,Y,1]^T M =[X,Y,1]T

结合式（3），可得
H = λ A [ r 1 r 2 t ] (5) \text H=\lambda \text A\left[\begin{matrix}\text r_1&\text r_2&\text t\end{matrix}\right] \tag 5 H=λA[r1r2t](5)
H \text H H是一个 3 × 3 3\times 3 3×3的矩阵， λ \lambda λ 是一个任意标量（ λ \lambda λ的存在是因为齐次坐标的尺度不变性，你也可以认为 λ \lambda λ就等于1）。

对内参的约束

承上，我们定义 H = [ h 1 h 2 h 3 ] \text H=\left[\begin{matrix}\text h_1&\text h_2&\text h_3\end{matrix}\right] H=[h1h2h3]，则有
[ h 1 h 2 h 3 ] = λ A [ r 1 r 2 t ] (6) \left[\begin{matrix}\text h_1&\text h_2&\text h_3\end{matrix}\right]=\lambda \text A\left[\begin{matrix}\text r_1&\text r_2&\text t\end{matrix}\right] \tag 6 [h1h2h3]=λA[r1r2t](6)

由于旋转 R \text R R矩阵是正交阵， r 1 \text r_1 r1和 r 2 \text r_2 r2是标准正交基，也就是它们两不仅相互垂直，而且模长均为1。即满足 r 1 T r 2 = 0 \text r_1^T\text r_2=0 r1Tr2=0， r 1 T r 1 = r 2 T r 2 = 1 \text r_1^T\text r_1=\text r_2^T\text r_2=1 r1Tr1=r2Tr2=1。

而上式可推得 r 1 = 1 λ A − 1 h 1 \text r_1=\frac 1 \lambda \text A^{-1}\text h_1 r1=λ1A−1h1， r 2 = 1 λ A − 1 h 2 \text r_2=\frac 1 \lambda \text A^{-1}\text h_2 r2=λ1A−1h2。

很容易得到：
h 1 T A − T A − 1 h 2 = 0 h 1 T A − T A − 1 h 1 = h 2 T A − T A − 1 h 2 (7) \begin{aligned} \text h_1^T\text A^{-T}\text A^{-1}\text h_2&=0\\ ~\\ \text h_1^T\text A^{-T}\text A^{-1}\text h_1&=\text h_2^T\text A^{-T}\text A^{-1}\text h_2 \end{aligned} \tag 7 h1TA−TA−1h2 h1TA−TA−1h1=0=h2TA−TA−1h2(7)
由上式可知，单应性矩阵 H \text H H和内参矩阵 A \text A A的元素之间满足两个线性方程约束。如果能够计算出单应性矩阵，直观上感觉应该可通过解线性方程解出内参矩阵。具体是否可以呢？我们继续往下看。

相机参数求解

本小节将带大家了解相机参数的具体解法，回答上节遗留的问题，即是否能够通过公式（7）来求解相机参数。

本小节将分为3个部分，第一部分介绍直接由式（7）求解相机参数的闭合解公式；第二部分介绍相机参数的最大似然解法；第三部分介绍考虑相机畸变后的畸变参数解法。

闭合解

首先定义 B = A − T A − 1 = [ B 11 B 12 B 13 B 12 B 22 B 23 B 13 B 23 B 33 ] \text B=\text A^{-T}\text A^{-1}=\left[\begin{matrix}B_{11}&B_{12}&B_{13}\\B_{12}&B_{22}&B_{23}\\B_{13}&B_{23}&B_{33}\end{matrix}\right] B=A−TA−1=⎣⎡B11B12B13B12B22B23B13B23B33⎦⎤。

具体的，可以计算出 B \text B B的每个元素：
B = [ 1 α 2 − γ α 2 β v 0 γ − u 0 β α 2 β − γ α 2 β γ 2 α 2 β 2 + 1 β 2 − γ ( v 0 γ − u 0 β ) α 2 β 2 − v 0 β 2 v 0 γ − u 0 β α 2 β − γ ( v 0 γ − u 0 β ) α 2 β 2 − v 0 β 2 ( v 0 γ − u 0 β ) 2 α 2 β 2 + v 0 2 β 2 + 1 ] (8) \text B=\left[\begin{matrix}\frac 1 {\alpha^2}&-\frac \gamma {\alpha^2\beta}&\frac {v_0\gamma-u_0\beta}{\alpha^2\beta}\\-\frac \gamma {\alpha^2\beta}&\frac {\gamma^2}{\alpha^2\beta^2}+\frac 1 {\beta^2}&-\frac {\gamma(v_0\gamma-u_0\beta)}{\alpha^2\beta^2}-\frac {v_0} {\beta^2}\\\frac {v_0\gamma-u_0\beta}{\alpha^2\beta}&-\frac {\gamma(v_0\gamma-u_0\beta)}{\alpha^2\beta^2}-\frac {v_0} {\beta^2}&\frac {(v_0\gamma-u_0\beta)^2}{\alpha^2\beta^2}+\frac {v_0^2}{\beta^2}+1\end{matrix}\right] \tag 8 B=⎣⎢⎡α21−α2βγα2βv0γ−u0β−α2βγα2β2γ2+β21−α2β2γ(v0γ−u0β)−β2v0α2βv0γ−u0β−α2β2γ(v0γ−u0β)−β2v0α2β2(v0γ−u0β)2+β2v02+1⎦⎥⎤(8)
（建议同学们还是动手推导下公式（8），虽说看上去比较复杂，但推导其实很简单）

显然， B \text B B是一个对称矩阵，我们可以用一个6维的矢量来定义：
b = [ B 11 B 12 B 22 B 13 B 23 B 33 ] T (9) \text b=\left[\begin{matrix}B_{11}&B_{12}&B_{22}&B_{13}&B_{23}&B_{33}\end{matrix}\right]^T \tag 9 b=[B11B12B22B13B23B33]T(9)
定义单应性矩阵 H \text H H的第 i i i列向量为 h i = [ h i 1 h i 2 h i 3 ] T \text h_i=\left[\begin{matrix}h_{i1}&h_{i2}&h_{i3}\end{matrix}\right]^T hi=[hi1hi2hi3]T，则有
h i T B h j = v i j T b (10) \text h_i^T\text B\text h_j=\text v_{ij}^T\text b \tag {10} hiTBhj=vijTb(10)
其中， v i j = [ h i 1 h j 1 h i 1 h j 2 + h i 2 h j 1 h i 2 h j 2 h i 3 h j 1 + h i 1 h j 3 h i 3 h j 2 + h i 2 h j 3 h i 3 h j 3 ] T \text v_{ij}=\left[\begin{matrix}h_{i1}h_{j1}&h_{i1}h_{j2}+h_{i2}h_{j1}&h_{i2}h_{j2}&h_{i3}h_{j1}+h_{i1}h_{j3}&h_{i3}h_{j2}+h_{i2}h_{j3}&h_{i3}h_{j3}\end{matrix}\right]^T vij=[hi1hj1hi1hj2+hi2hj1hi2hj2hi3hj1+hi1hj3hi3hj2+hi2hj3hi3hj3]T。

依次，公式(7)可以转换为：
[ v 12 T ( v 11 − v 22 ) T ] b = 0 (11) \left[\begin{matrix}\text v_{12}^T\\(\text v_{11}-\text v_{22})^T\end{matrix}\right]\text b=0 \tag {11} [v12T(v11−v22)T]b=0(11)

这是一个典型的线性方程组，系数矩阵的行数为2，求解6维向量 b \text b b。由公式（4），当相机在1个位姿下拍摄标定板图案后，经过角点的像素坐标提取，可得所有角点的世界坐标系和像素坐标系的对应关系，进而通过线性方程组的最小二乘解法求解当前位姿下的单应性变换矩阵 H \text H H，可得公式（11）的具体表达式。

但公式（11）的系数矩阵只有2行，要求解6维向量 b \text b b是不够的。所以我们需要相机在 n n n个位姿下拍摄标定图案，得到 n n n个单应性矩阵，以及行数为 2 n 2n 2n的系数矩阵，当 n > = 3 n>=3 n>=3时，便可求解6维向量 b \text b b。也就是说至少3张图片才能完成相机标定。最后得到的总方程组可表达为：
Vb = 0 (12) \text V\text b=0 \tag {12} Vb=0(12)

V \text V V是一个 2 n × 6 2n\times6 2n×6的系数矩阵。

求解公式（12）常用的方法，其一是矩阵 V T V \text V^T\text V VTV的最小特征值对应的特征向量即为方程解，其二是对系数矩阵 V \text V V进行奇异值分解 V = U S V \text V=USV V=USV，分解矩阵 V V V的第三列即为方程解。

需要注意的是，公式（12）是尺度等价的，求解出的 b \text b b乘以任何倍数依旧是正确解，所以由 b \text b b组成的矩阵 B \text B B并不是严格满足 B = A − T A − 1 \text B=\text A^{-T}\text A^{-1} B=A−TA−1，而是存在一个任意的尺度因子 λ \lambda λ，满足 B = λ A − T A − 1 \text B=\lambda\text A^{-T}\text A^{-1} B=λA−TA−1。通过推导，可通过公式（13）来计算所有的相机内参数：
v 0 = ( B 12 B 13 − B 11 B 23 ) / ( B 11 B 22 − B 12 2 ) λ = B 33 − [ B 13 2 + v 0 ( B 12 B 13 − B 11 B 23 ) ] / B 11 α = λ / B 11 β = λ B 11 / ( B 11 B 22 − B 12 2 ) γ = − B 12 α 2 β / λ u 0 = γ v 0 / β − B 13 α 2 / λ (13) \begin{aligned} v_0&=(B_{12}B_{13}-B_{11}B_{23})/(B_{11}B_{22}-B_{12}^2)\\ \lambda&=B_{33}-[B_{13}^2+v_0(B_{12}B_{13}-B_{11}B_{23})]/B_{11}\\ \alpha&=\sqrt {\lambda/B_{11}}\\ \beta&=\sqrt{\lambda B_{11}/(B_{11}B_{22}-B_{12}^2)}\\ \gamma&=-B_{12}\alpha^2\beta/\lambda\\ u_0&=\gamma v_0/\beta -B_{13}\alpha^2/\lambda \end{aligned} \tag {13} v0λαβγu0=(B12B13−B11B23)/(B11B22−B122)=B33−[B132+v0(B12B13−B11B23)]/B11=λ/B11 =λB11/(B11B22−B122) =−B12α2β/λ=γv0/β−B13α2/λ(13)
当 A \text A A求解出后，每个位姿下的外参数 R , t \text R,\text t R,t可以通过公式（14）计算：
r 1 = λ A − 1 h 1 r 2 = λ A − 1 h 2 r 3 = r 1 × r 2 t = λ A − 1 h 3 (14) \begin{aligned} \text r_1&=\lambda \text A^{-1} \text h_1\\ \text r_2&=\lambda \text A^{-1} \text h_2\\ \text r_3&=\text r_1\times \text r_2\\ \text t&=\lambda \text A^{-1}\text h_3 \end{aligned}\tag {14} r1r2r3t=λA−1h1=λA−1h2=r1×r2=λA−1h3(14)
其中， λ = 1 / ∣ ∣ A − 1 h 1 ∣ ∣ = 1 / ∣ ∣ A − 1 h 2 ∣ ∣ \lambda=1/||\text A^{-1}\text h_1||=1/||\text A^{-1}\text h_2|| λ=1/∣∣A−1h1∣∣=1/∣∣A−1h2∣∣。

需要注意的是，由于噪声的存在，计算出来的旋转矩阵 Q = [ r 1 , r 2 , r 3 ] \text Q=[\text r_1,\text r_2,\text r_3] Q=[r1,r2,r3]并不满足正交性，需要进一步计算和 Q \text Q Q最接近的正交矩阵 R \text R R，可通过奇异值分解来计算。最佳的正交矩阵 R \text R R应该满足矩阵 R − Q \text R-\text Q R−Q的Frobenius范数最小，即求解如下问题：
min ⁡ R ∣ ∣ R − Q ∣ ∣ F 2 subject to R T R = I (15) \min_\text R||\text R-\text Q||_F^2 \ \ \ \ \text{subject to} \ \text R^T\text R=\text I \tag {15} Rmin∣∣R−Q∣∣F2 subject to RTR=I(15)
因为
∣ ∣ R − Q ∣ ∣ F 2 = t r a c e ( ( R − Q ) T ( R − Q ) ) = 3 + t r a c e ( Q T Q ) − 2 t r a c e ( R T Q ) \begin{aligned} ||\text R-\text Q||_F^2&=trace((\text R-\text Q)^T(\text R-\text Q))\\ &=3+trace(\text Q^T\text Q)-2trace(\text R^T\text Q) \end{aligned} ∣∣R−Q∣∣F2=trace((R−Q)T(R−Q))=3+trace(QTQ)−2trace(RTQ)
所以，问题转换为求解让 t r a c e ( R T Q ) trace(\text R^T\text Q) trace(RTQ)最大的 R \text R R。
我们定义矩阵 Q \text Q Q的的奇异值分解为 US V T \text U\text S\text V^T USVT，其中 S = d i a g ( σ 1 , σ 2 , σ 3 ) \text S=diag(\sigma_1,\sigma_2,\sigma_3) S=diag(σ1,σ2,σ3)。如果我们定义正交矩阵 Z = V T R T U \text Z=\text V^T\text R^T\text U Z=VTRTU，则有
t r a c e ( R T Q ) = t r a c e ( R T US V T ) = t r a c e ( V T R T US ) = t r a c e ( ZS ) = ∑ i = 1 3 z i i σ i ≤ ∑ i = 1 3 σ i \begin{aligned} trace(\text R^T\text Q)&=trace(\text R^T\text U\text S\text V^T)=trace(\text V^T\text R^T\text U\text S)\\ &=trace(\text Z\text S)=\sum_{i=1}^3z_{ii}\sigma_{i}\le\sum_{i=1}^3\sigma_i \end{aligned} trace(RTQ)=trace(RTUSVT)=trace(VTRTUS)=trace(ZS)=i=1∑3ziiσi≤i=1∑3σi
显然，当 Z = I \text Z=\text I Z=I时， t r a c e ( R T Q ) trace(\text R^T\text Q) trace(RTQ)达到最大值，可得 R = U V T \text R=\text U\text V^T R=UVT。

以上，内外参数全部解出，整个闭合解法完成。

最大似然估计

闭合解法可得到代数距离最小的解，并没有考虑到各个参数实际的几何含义，由于噪声的存在，解并不会非常精确。我们可以通过最大似然估计法来获取更精确的解。

如果我们观测到标定板 n n n张拍摄图像中的 m m m个点对。假设所有点存在独立的等尺度的噪声，则最大似然估计即最小化如下表达式：
∑ i = 1 n ∑ j = 1 m ∣ ∣ m i j − m ^ ( A , R i , t i , M j ) ∣ ∣ 2 (16) \sum_{i=1}^n \sum_{j=1}^m ||\text m_{ij} - \hat {\text m} (\text A,\text R_i,\text t_i,\text M_j)||^2 \tag {16} i=1∑nj=1∑m∣∣mij−m^(A,Ri,ti,Mj)∣∣2(16)
其中， m ^ ( A , R i , t i , M j ) \hat {\text m}(\text A,\text R_i,\text t_i,\text M_j) m^(A,Ri,ti,Mj)是空间点 M j \text M_j Mj在图像 i i i上的投影点，可由公式（2）计算。旋转矩阵 R \text R R可通过三维向量 r \text r r来表达， r \text r r和旋转轴平行且模长为旋转角， R \text R R和 r \text r r可通过罗德里格式公式来相互转换。

式（16）是一个非线性表达式，所以这是一个非线性最小化求解问题，可以用Levenberg-Marquardt算法来求解。想必部分同学知道求解该类问题需要一个相对准确的初始化值，这个值便可以使用前一节得到的闭合解来获取。

在工程实践中，我们通常用ceres-solver库来搞定这部分。

相机畸变

到目前为止，我们还没有考虑相机畸变，所有推导都是基于无畸变的理想情况，而实际情况是相机必然会有多多少少的畸变，主要包括两种：径向畸变和切向畸变。一般情况下，会考虑3项径向畸变 k 1 , k 2 , k 3 k_1,k_2,k_3 k1,k2,k3和2项切向畸变 p 1 , p 2 p_1,p_2 p1,p2。在Zhang式标定法中，只考虑了2项径向畸变 k 1 , k 2 k_1,k_2 k1,k2，而实际应用时，我们会考虑更多项，原理相同，依次类推即可。

定义 ( u , v ) (u,v) (u,v)为理想的像素坐标， ( u ˘ , v ˘ ) (\breve u,\breve v) (u˘,v˘)为实际观测到的像素坐标，则考虑2项径向畸变后，两者满足如下公式：
u ˘ = u + ( u − u 0 ) [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] v ˘ = v + ( v − v 0 ) [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] (17) \begin{aligned} \breve u&=u+(u-u_0)[k_1(x^2+y^2)+k_2(x^2+y^2)^2]\\ \breve v&=v+(v-v_0)[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \tag {17} \end{aligned} u˘v˘=u+(u−u0)[k1(x2+y2)+k2(x2+y2)2]=v+(v−v0)[k1(x2+y2)+k2(x2+y2)2](17)
其中， x = u − u 0 x=u-u_0 x=u−u0， y = v − v 0 y=v-v_0 y=v−v0。

式（17）是一个线性方程组，求解 k 1 , k 2 k_1,k_2 k1,k2看上去是一件简单的事儿，但实际上却面临一个问题，即无法获取无畸变的理想坐标 ( u , v ) (u,v) (u,v)，它需要已知准确的内外参通过投影公式（2）来求解，同学们可能会问：前两节不是解出了内外参数吗！？不要忘了，前面是没有考虑相机畸变的，内外参的求解实际是不精确的（把观测值当做无畸变坐标来求解的内外参，是有偏差的）。所以这成了一个鸡生蛋蛋生鸡的矛盾问题。

那就干脆不奢求通过线性方程求出精确畸变，将闭合解得到的内内外参代入公式（2）求出近似的理想坐标 ( u , v ) (u,v) (u,v)，从而由公式（17）建立线性方程组求解近似的 k 1 , k 2 k_1,k_2 k1,k2。

然后，我们建立新的最大似然估计表达式：
∑ i = 1 n ∑ j = 1 m ∣ ∣ m i j − m ˘ ( A , k 1 , k 2 , R i , t i , M j ) ∣ ∣ 2 (18) \sum_{i=1}^n \sum_{j=1}^m ||\text m_{ij} - \breve {\text m}(\text A,k_1,k_2,\text R_i,\text t_i,\text M_j)||^2 \tag {18} i=1∑nj=1∑m∣∣mij−m˘(A,k1,k2,Ri,ti,Mj)∣∣2(18)
同样通过非线性求解方法来求解所有内外参数和畸变系数。而畸变系数的初值由上面所描述的方法来求解。

实际上，因为畸变系数是非常小值，所以直接将畸变系数的初值全部设置为0也是可以的，这就不需要解线性方程组了。

总结

下面对整个Zhang式标定的流程做一个总结：

打印标定图案并粘贴至一个平面上，称之为标定板。
通过移动相机或移动标定板在不同的位姿拍摄多张标定板图像（图像数>=3）。
在所有图像上检测特征点（角点或者圆心点）。
使用闭合解法求解所有内参数和外参数。
通过线性方程组求解近似的畸变系数（或者直接赋值为0）。
通过非线性优化计算精确的内外参数和畸变系数。

以上便是Zhang式标定法的所有理论讲解，希望同学们读完此篇能有醍醐灌顶之效，并能融会贯通，在实际工程应用中得心应手。

在Zhang式标定法的原文中，还有一些值得学习的高阶知识点并未在本篇中体现，包括关键公式（7）的几何解释，以及退化情况下的讨论，感兴趣的同学可以继续阅读原文来了解。

下一篇我们将带来相机标定的另一种方法：DLT直接线性变换法。

作业

这里为大家准备了一些练习题，可以通过实践加深理解：

1 通过opencv开源库提供的接口完成相机标定，可使用opencv提供的图像，也可使用自己的图像。
2 更高阶的是，你能够自己不依赖opencv库写一套相机标定算法吗？或者只使用opencv来检测角点坐标，其他步骤自己来实现。

参考答案地址：https://github.com/ethan-li-coding/StereoV3DCode [代码持续更新，大家感兴趣可以star和watch，本篇代码暂缺]

参考文献

- [1] Zhang Z . A Flexible New Technique for Camera Calibration[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000, 22(11):1330-1334.