立體視覺入門指南（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[R?t?][Pw?1?](1)
所謂標定，即是由大量觀測值擬合引數模型的程序，且在此擬合的引數模型是已知的，所以應盡可能探索能便捷獲取大量觀測值的方案，如果觀測值之間還滿足一些其他的幾何約束就更有助于求解具體單個引數值，

今天所述的Zhang式標定法¹即提供了一種便捷獲取大量觀測值的的方案，同時觀測值之間還滿足一類明顯的幾何約束（即平面約束），可直接求解出內外參，其操作方式非常簡單，只需要拍攝帶有標定板圖案的平面，即可完成相機標定，使標定難度極大降低，如果不追求高精度，列印一張棋盤格標定板圖案粘貼到近似平的硬紙板上即可完成標定，加快了立體視覺的入門和普及，影響深遠，是相機標定領域絕對的經典，

本篇即帶大家深入了解Zhang式相機標定法，掌握本篇對立體視覺的理論掌握及工程實戰來說都是非常必要的，

實施方法

Zhang式標定法能夠被廣泛應用，其中一個重要原因是其實施方法十分簡單，不需要專業的工藝制作即可完成，

第一步，設計一張具有明顯角點特征，且已知每個角點二維坐標的圖案作為標定圖案，常見的圖案有三種：

棋盤格

ChArUco

圓點

規則的圖案設計可以方便的計算出角點在圖案內的二維坐標，拿棋盤格來說，角點之間的間隔像素數是固定的，假設左上角角點的坐標為 ( 0 , 0 ) (0,0) (0,0)，則其他角點的像素坐標都可以通過格子的偏移量計算出來，而一張已知DPI的標定板影像，在列印后每個角點的二維空間坐標也是完全已知的（通過像素換算成空間尺寸），

第二步，將標定板圖案以某種方式置于一個平面上，比如最簡單的方式是將標定圖原尺寸列印出來，然后找一塊近似平的平板，將列印后的標定圖案粘貼至平板上；更專業高精度的方式是找專業廠家制作高精平板（如陶瓷板）并將標定圖案以某種工藝刻印到平板上，這一步的目的是讓標定圖案的角點都位于一個平面上，

如此，第一步所描述的二維坐標可以轉換成第三維 Z Z Z坐標等于0的三維坐標（將世界坐標系的原點放在標定板的某個角點，Z軸垂直于標定板），

第三步，移動相機到N（N>=3）個不同的位姿拍攝標定板圖案，

第四步，對上一步拍攝的標定板圖案進行角點提取，解算標定引數，

以上便是相機標定的實施步驟，總結來說，在一個平面標定板上有一組已知空間坐標的角點，相機在多個不同位姿下拍攝角點圖案并提取角點的像素坐標，即可完成相機內外引數的解算，這句總結蘊含2個重要的知識點：

1. 標定圖案中角點的空間坐標是已知的，且它們都位于一個平面上， Z Z Z坐標等于0
2. 相機需要在多個不同的位姿拍攝角點圖案并提取像素坐標

可以看到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[R?t?]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?γβ0?u0?v0?1????
其中， ( 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[R?t?]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[r1??r2??r3??t?]?????XY01??????=A[r1??r2??t?]???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[r1??r2??t?](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=[h1??h2??h3??]，則有
[ 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 [h1??h2??h3??]=λA[r1??r2??t?](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 r1T?r2?=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 r1T?r1?=r2T?r2?=1，

而上式可推得 r 1 = 1 λ A ? 1 h 1 \text r_1=\frac 1 \lambda \text A^{-1}\text h_1 r1?=λ1?A?1h1?， r 2 = 1 λ A ? 1 h 2 \text r_2=\frac 1 \lambda \text A^{-1}\text h_2 r2?=λ1?A?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 \tag 7 \end{aligned} h1T?A?TA?1h2? h1T?A?TA?1h1??=0=h2T?A?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=???B11?B12?B13??B12?B22?B23??B13?B23?B33?????，

具體的，可以計算出 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=[B11??B12??B22??B13??B23??B33??]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?=[hi1??hi2??hi3??]T，則有
h i T B h j = v i j T b (9) \text h_i^T\text B\text h_j=\text v_{ij}^T\text b \tag 9 hiT?Bhj?=vijT?b(9)
其中， 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?=[hi1?hj1??hi1?hj2?+hi2?hj1??hi2?hj2??hi3?hj1?+hi1?hj3??hi3?hj2?+hi2?hj3??hi3?hj3??]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 \text B=\text A^{-T}\text A B=A?TA，而是存在一個任意的尺度因子 λ \lambda λ，滿足 B = λ A ? T A \text B=\lambda\text A^{-T}\text A B=λA?TA，通過推導，可通過公式（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 \tag {13} \end{aligned} v0?λαβγu0??=(B12?B13??B11?B23?)/(B11?B22??B122?)=B33??[B132?+v0?(B12?B13??B11?B23?)]/B11?=λ/B11? ?=λB11?/(B11?B22??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 \tag {14} \end{aligned} r1?r2?r3?t?=λ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∑3?zii?σ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∑n?j=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∑n?j=1∑m?∣∣mij??m?(A,k1?,k2?,Ri?,ti?,Mj?)∣∣2(18)
同樣通過非線性求解方法來求解所有內外引數和畸變系數，而畸變系數的初值由上面所描述的方法來求解，

實際上，因為畸變系數是非常小值，所以直接將畸變系數的初值全部設定為0也是可以的，這就不需要解線性方程組了，

總結

下面對整個Zhang式標定的流程做一個總結：

1. 列印標定圖案并粘貼至一個平面上，稱之為標定板，
2. 通過移動相機或移動標定板在不同的位姿拍攝多張標定板影像（影像數>=3），
3. 在所有影像上檢測特征點（角點或者圓心點），
4. 使用閉合解法求解所有內引數和外引數，
5. 通過線性方程組求解近似的畸變系數（或者直接賦值為0），
6. 通過非線性優化計算精確的內外引數和畸變系數，

以上便是Zhang式標定法的所有理論講解，希望同學們讀完此篇能有醍醐灌頂之效，并能融會貫通，在實際工程應用中得心應手，

在Zhang式標定法的原文中，還有一些值得學習的高階知識點并未在本篇中體現，包括關鍵公式（7）的幾何解釋，以及退化情況下的討論，感興趣的同學可以繼續閱讀原文來了解，

下一篇我們將帶來相機標定的另一種方法：DLT直接線性變換法，

作業

這里為大家準備了一些練習題，可以通過實踐加深理解：

1 通過opencv開源庫提供的介面完成相機標定，可使用opencv提供的影像，也可使用自己的影像，
2 更高階的是，你能夠自己不依賴opencv庫寫一套相機標定演算法嗎？或者只使用opencv來檢測角點坐標，其他步驟自己來實作，

參考答案地址：https://github.com/ethan-li-coding/StereoV3DCode [代碼持續更新，大家感興趣可以star和watch，本篇代碼暫缺]

參考文獻