一種標定相機內方位元素的迭代算法

王東杰 張繼友 馬麗娜

（北京空間機電研究所，北京 100094）

0 引言

相機標定技術是從拍攝的圖像信息和標定參照物中提取相機的內外方位元素，從而建立空間三維坐標系中位置、方向信息與圖像二維坐標系間的精確映射關系。相機的內方位元素指相機的主點、主距、畸變等參數，外方位元素指相機位置和方向等參數[1]。相機內外方位元素的標定在攝影測量、圖像校正、三維重建以及機器視覺等諸多領域具有重要作用。相機自標定方法為近幾年興起的一種相機標定技術，該方法由一臺或多臺相機從不同的位置和方向對一個尺寸已知的棋盤格模板拍攝多幅圖像，以模板和圖像的對應關系求解出內外方位元素，具有標定操作過程簡單易行、實驗結果精度較高的優點。

本文介紹一種利用消失點理論得到迭代算法初值以求解相機內外方位元素的B雙空間幾何自標定方法。不同于張氏標定法[2]和Tsai的兩步標定法僅利用標定圖像的點信息，B雙空間幾何法基于消失點理論，利用二維棋盤格標定圖像中點、線、面間固有的幾何關系獲得相機參數的初始估計值[3]，采用牛頓—拉夫遜（New ton-Raphson）算法對非線性的共線方程進行迭代求解。New ton-Raphson迭代算法能夠在保證計算效率的同時獲得較高的結果精度[4-5]，圖像中的棋盤點坐標由Harris角點提取算法得到，作為已知條件，通過使用張氏標定法中的實驗數據進行對比性實驗，實驗結果表明，兩種方法得到的結果標定精度十分接近，證實本文算法能夠用于相機內外方位元素的標定，標定精度較理想。

1 迭代算法

本文利用New ton-Raphson算法對由相機參數的初始估計值得到的非線性共線方程進行迭代求解，以得到精確的相機參數標定結果。相機參數的初始估計值由B雙空間幾何理論得到，而本文所采用的畸變模型包含3個徑向和2個切向畸變系數，比張氏標定法只包括2個徑向畸變系數的畸變模型能夠更精準地給出畸變值。整個迭代算法的流程如圖1所示。

圖1 迭代算法流程Fig.1 The flow chartof the iterativemethod

1.1 共線方程

設模板坐標系和圖像坐標系中點的坐標分別為（X,Y,Z）和（x,y,z）。圖像坐標系的原點在模板坐標系中的坐標為（XC,YC,ZC）。圖像坐標系相對于模板坐標系的方向由一個3×3的旋轉矩陣R表示。R的9個元素rij（i=1,2,3；j=1,2,3）非獨立，通過歐拉角變換，R可以視為分別繞3個坐標軸旋轉的3個旋轉矩陣的乘積[6]。設繞x,y,z軸旋轉的角度分別為ω，φ，κ，則

式中

以上3個方程相乘后得：

由式（2）可知，R的各行（列）向量相互正交。

在標定過程中，模板中取棋盤角點的左上點為原點，坐標值為實際距離值。對于理想相機模型，模板上一點（X,Y,Z）投影到圖像坐標系后其坐標為（x,y）。但由于圖像畸變，其實際坐標為（x′,y′）。

模板坐標系和圖像坐標系之間的關系為

式中 Δx,Δy為包括對稱和非對稱畸變的畸變模型；xp,yp為圖像坐標系中的主點；fx,fy為主距，fx=fdx，fy=fdy，f為相機焦距，dx,dy分別為像面上X,Y方向上的像元尺寸。Δx,Δy的計算式為

設棋盤模板上的棋盤角點數量為m，相機從不同位置拍攝的模板圖像數量為n。對于每個相機的放置位置，其對應的圖像上能檢測到并能計算出坐標的棋盤角點的數量為mi（i=1,2,3,…,n）個，有mi≤m，因此所有未知數數量為6n+9，所有共線方程的數量為，且不超過所拍攝所有圖像中的棋盤角點數量總數2mn，在檢驗過程中的超定方程為

1.2 求解算法

1.2.1 迭代計算

與內外方位元素有關的超定方程可使用最小二乘法求解。首先將方程改寫成一個簡易的形式，定義一個包括所有未知量的向量為

式中l=mi-1+j，m0=0，j=1,2,3,…,mi，i=1,2,3,…,n；χj=(Xj,Yj,Zj)為模板坐標系中第j個棋盤點的坐標向量；χC,i=(XC,i,YC,i,ZC,i)為相機第i個位置所拍攝圖像的坐標原點在模板坐標系中的坐標向量。單位向量rk,i為第i組歐拉角(ωi,φi,κi)對應的旋轉矩陣R中的第k列列向量[6]。

對式（9）用迭代算法進行求解，假設解向量ξ在第k次迭代得到，則式（9）左邊可以用第k+1次迭代時的有限泰勒級數表示：

式中 Δξk為一個包含6n+10個元素的列向量；δ為一個包含個元素的列向量；行、6n+10列的雅克比矩陣。根據New ton-Raphson迭代算法，令δk+1=0，得到：

式中wj為由矩陣A的奇異值組成對角矩陣中的第j個元素。V■■di ag(1/wj)■■UT為矩陣A的偽逆矩陣。在本文中，雅克比矩陣的偽逆陣通過奇異值分解法得到，在ξk+1處的解向量通過最小二乘法得到。使用奇異值分解法的優點是當雅克比矩陣為奇異矩陣或病態矩陣時，上述過程依然有效。計算雅克比矩陣的條件數，即最大與最小奇異值的比值，若條件數的倒數小于10–p（p為條件數評價標準的指數），則將所有小于最大奇異值10–p倍的奇異值設為0，可以排除舍入誤差對解造成的干擾，本文中的計算為雙精度計算，令p=12。

1.2.3 收斂標準

文獻[6]中指出New ton-Raphson方法在初始估計的迭代值非常接近方程的解時，即迭代結果已滿足精度要求時收斂。本文結束迭代計算過程的收斂標準為

式中ε為一個視計算情況而定的一個無限小的數。

另一種收斂標準是檢驗解向量 Δξk的L∞范數的變化[6]。對于上述兩種標準，都在ε= 10-8、迭代 10～20次時達到收斂。收斂結果不包含二次偏導項，與標準的New ton-Raphson算法不同，這是因為在采用最小二乘法解方程組時，有些情況下雅克比矩陣可能有較大的條件數。

1.2.4 萬向節死鎖

φ= π/2時，旋轉矩陣為

式（14）說明，解向量僅依靠ω和κ的變化，獨立未知量的數量減少了一個，這種現象稱為“萬向節死鎖”。對于這種情況，檢驗結果顯示收斂的解向量中ω和κ為隨機值，但二者的差值正確。對于φ= π/2的萬向節死鎖情況無需特殊措施。由式（14）可以看出，相比于ω、κ各自的值，ω-κ的差值更重要，因為最終要得到模板點坐標和圖像點坐標之間的映射關系[6-7]。

1.3 求解算法

在透視投影下，三維空間中的平行線映射到圖像平面上相交于一點，稱為消失點，不同消失點的連線稱為消失線。不同于成像平面上的其他特征點，消失點和消失線蘊含了直線的方向信息，比邊緣、拐角等特征更具魯棒性。對消失點和消失線的分析可以提供大量的場景三維結構和方向信息[7]。本文利用B雙空間幾何方法對相機內外方位元素進行求解，作為迭代算法的初始估計值。B雙空間幾何法通過B雙空間的轉換將消失點和消失線映射成容易計算的有限向量[8-9]，利用平行和正交的向量關系可以計算出消失點的值，進而得到相機內外方位元素，作為迭代算法的初值[3,9-11]。

1.3.1 焦距值的確定

B雙空間幾何中點線面關系如圖2所示，假設正方形ABCD經透視投影后，正方形的4個點在以相?機主點為原點的坐標系中的坐標分別為，其坐標值（x1,y1,1），（x2,y2,1），（x3,y3,1），（x4,y4,1）由角點提取算法給出，坐標值中的“1”為在齊次坐標系中的有限遠平面，點V1為向量的消失點，點V2為向量的消失點，為與面ΠH相關的消失線。

圖2 雙空間幾何中點線面關系Fig.2 The relationship among points,lines and planes in B-dual-space geometry

則由雙空間幾何關系的特性可得到：

由于ABCD為正方形，其對角線相互垂直，利用這條約束可得到兩對角線與消失線間的關系為

由于V1⊥V2、V3⊥V4，可得：

式中V?［abc］T為各消失點的像元坐標，a、b、c為像元坐標值，i=1,2,3,4。由式（17）

iiiiiii可得：

至此可得到2個主距的初始估計。

1.3.2 其他參數的確定

由文獻[8,11-12]可知，在利用棋盤模板進行標定的過程中，不能確定主點位置（xp,yp）。因此在初始估計階段，將主點位置假設為像面的中心，即為像面X、Y方向上的像元個數。通過迭代算法求解出主點的實際值。

根據文獻[8]中雙空間幾何理論中總結的消失線特性可知，消失線向量正比于平面ΠH的法向量，因此可由其得到棋盤模板平面相對于圖像平面的方向，即可通過向量得到外方位元素中表示相機方向的角度值(ω，φ，κ)。此外，文獻[8]中指出模板坐標系的原點與主點之間的距離也可通過消失線向量λH求出，于是可得到外方位元素中表示相機位置向量（XC，YC，ZC），相機的外方位元素的初始估計值為(XC,YC,ZC,ω,φ,κ)。

1.4 重投影誤差

作為標定結果精度的衡量標準，重投影誤差定義為圖像點檢測值（xi，yi）和根據內外方位元素(k1,k2,k3,p1,p1,xp,yp,fx,fy)、模板點坐標（Xi,Yi）帶入畸變模型后進行計算得到的圖像坐標點（xi′,yi′）之間的偏差。重投影誤差的計算通過建立數據擬合目標函數得到：

從上式看出，目標函數是由一組殘量平方和組成，若得到的目標函數值越小，則數據擬合得越好，求得的參數值越準確。本文算法不是直接求解非線性的畸變模型，而是將其轉化成求解非線性最小二乘問題，通過非線性優化算法多次迭代，最后得到使目標函數值最小的參數值，降低了求解難度。

2 與張氏標定法的比較

張正友提出了基于二維平面模板的標定方法（張氏標定法）[13-14]，張氏標定法只要從不同角度對同一標定平面(標定板)拍攝2幅以上的圖像，就可以求出相機的內外方位元素。由于該方法不需要知道標定板移動的具體位置信息，而且標定板的制作簡單，因此這種方法簡單、靈活。其整個自標定過程可分解為以下幾個步驟：

1）分別以圖像幅數范圍、網格板角度范圍等參數拍攝照片；

2）對所拍攝照片進行圖像處理，用Harris角點提取算法提取角點；

3）利用棋盤格標定板在世界坐標系的坐標和相應的圖像角點坐標計算每幅圖像的單應性矩陣；

4）由所得單應性矩陣，根據自標定算法求解出內方位元素矩陣初值；

5）設計畸變模型；

6）將由內方位矩陣和畸變模型計算所得的角點坐標和實際圖像角點坐標代入數據擬合函數進行擬合，并通過非線性優化算法進行迭代，得到內方位矩陣和畸變系數的精確值。

本文提出的基于B雙空間幾何的相機自標定方法與張氏標定法類似，所需靶標同樣為尺寸精準的棋盤格標定板，標定過程也同樣需要上述6個步驟。但不同之處在于：

1）確定內方位元素矩陣初值的算法；

2）畸變模型；

3）非線性優化迭代算法。

兩種算法的異同點見表1。

表1 兩種算法的比較Tab.1 Comparison between Zhang’s and our calibration algorithms

由表1可知，本文提出的算法通過基于消失點理論的B雙空間幾何法，充分利用二維棋盤格標定圖像中點、線、面間固有的幾何關系獲得相機參數的初始估計值；通過New ton-Raphson迭代算法在保證非線性優化迭代計算效率的同時獲得較高的結果精度；畸變模型能更加詳細地給出畸變值。

3 實驗

對于本文提出的標定方法，利用文獻[2]提供的測試數據進行了實驗，如圖3所示，并將測試結果進行對比，見表2。張氏標定法基于Levenberg-Marquardt算法對得到的內方位元素初值進行優化迭代求解，并用相同的數據目標擬合函數進行重投影誤差計算，作為評價標定結果精度的標準。

圖3 張氏標定法使用的標定圖像Fig.3 Calibration Images of Zhangmethod

表2 張氏標定法與本文算法的結果及誤差比較Tab.2 Comparison between Zhang’sand our calibration algorithm s 像元

由表2可以看出，雖然兩種算法用于優化的初值相差較多，但經過優化后得到的終值十分相近。但張氏標定法得到的重投影誤差為0.335個像元[14]，本文算法得到的重投影誤差為0.265個像元，這表示本文方法計算出的相機內方位元素精度更高，更加接近實際相機成像模型。由于畸變模型不同，鏡頭畸變系數不具有可比性，但是本文的畸變模型有5項參數，能夠更加準確反應實際鏡頭的畸變情況。

4 結束語

本文介紹了一種利用相機在不同位置拍攝的多幅圖像中的特征點信息，計算相機內外方位元素的迭代計算方法。該方法用B雙空間幾何的特性，利用消失點、消失線理論得到迭代方法的初值。在無需進行干預初始估計值的情況下，通過迭代計算得到精度較高的內方位元素值。通過與張氏標定法的對比實驗來比較兩標定法的精度，實驗結果表明，本文介紹的算法誤差更小，其標定的相機參數也更加精確，能夠用于相機內方位元素的標定，且標定精度較理想。

References)

[1] 譚啟蒙,胡成威,高升.空間機械臂視覺相機內參標定技術研究[J].航天返回與遙感,2013,34(6):74-81.TAN Qimeng,HU Chengwei,GAO Sheng.Research on Calibration of Intrinsic Parameters for Space Manipulator Camera Based on 2D Planar Pattern[J].Spacecraft Recovery&Remote Sensing,2013,34(6):74-81.(in Chinese)

[2] ZHANG Zhengyou.A Flexible New Technique for Camera Calibration[EB/OL].[2014-03-20].http://research.microsoft.com/users/zhang/Calib/Calibration/Calib.txt.

[3] 吳剛,唐振民.B雙空間幾何中基于消隱點的攝像機標定[J].計算機工程與應用,2009,45(24):4-8.WU Gang,TANG Zhenmin.Camera Calibration Based on Vanishing Points in B-dual-space Geometry[J].Computer Engineering and Applications,2009,45(24):4-8.(in Chinese)

[4] 熊磊,劉國棟,劉小勇,等.數字圖像相關中基于非線性灰度改變模型的灰度梯度迭代算法[J].機床與液壓,2012,40(13):70-72.XIONG Lei,LIU Guodong,LIU Xiaoyong,et al.Gray-gradient Iterative Algorithm Based on Nonlinear Intensity Change Model in Digital Image Correlation[J].Machine Tool&Hydraulics,2012,40(13):70-72.(in Chinese)

[5] Marzan G,Karara H.A Computer Program for Direct Linear Transformation Solution of the Collinearity Condition and Some Applications of it[C].Proc.APSSymp.on Close-Range Photogrammetric Systems,Urbana,1975:420-476.

[6] Ravi S.A Method to Solve Interior and Exterior Camera Calibration Parameters for Image Resection[J/OL].[2014-03-20].http://www.nas.nasa.gov/Research/Reports/Techreports/1999/nas-99-003-abstract.htm l.

[7] 吳聞.游戲程序設計中若干問題的探討[J].電腦知識與技術,2009,5(16):4331-4336.WU Wen.Discussion of Some Problems in Game Programm ing[J].Computer Know ledge and Technology.2009,5(16):4331-4336.(in Chinese)

[8] 梅雪,夏良正,李久賢,等.一種三維場景的消失點檢測算法[J].信號處理,2007,23(6):924-926.MEIXu,XIA Liangzheng,LI Jiuxian,et al.A Vanishing Point Detection Algorithm for 3D Scene[J].Signal Processing,2007,23(6):924-926.(in Chinese)

[9] Bouguetj JY.VisualMethods for Three-dementional Modeling[D].California:California Institute of Technology,1999.

[10] Jean Y B,Pietro P.Camera Calibration from Points and Lines in Dual-Space Geometry[EB/OL].http://www.vision.caltech.edu/bouguetj/calib.doc.

[11] Jean Y B,Pietro P.Closed-form Camera Calibration in Dual-space Geometry[EB/OL].http://www.vision.caltech.edu/bouguetj/calib.doc.

[12] Tsai R Y.A Versatile Camera Calibration Technique for High Accuracy 3D Machine Vision Metrology Using off-the-shelf TV Cameras and Lenses[J].IEEE J.Rovotics Automat,1987,3(4):323-344.

[13] WANG L L,TsaiW H.Camera Calibration by Vanishing Lines for 3D Computer Vision[J].PAM I,1991,13(4):370-376.

[14] ZHANG Zhengyou.Flexible Camera Calibration by View ing a Plane From Unknown Orientations[J].IEEE International Conference on Computer Vision,1999(1):666-673.

[15] Abdel-Aziz Y,Karara H.Direct Linear Transformation from Comparator Coordinates into Object Space Coordinates[C]//Proc.APSSymp.on Close-Range Photogrammetric Systems,Urbana,IL.1971:1-48.

[16] Bell JH,M cLachlan BG.Image Registration for Pressure-Sensitive PaintApplications[J].Experiments in Fluids,1996(22):78-86.