一種標(biāo)定相機(jī)內(nèi)方位元素的迭代算法

王東杰 張繼友 馬麗娜

（北京空間機(jī)電研究所，北京 100094）

0 引言

相機(jī)標(biāo)定技術(shù)是從拍攝的圖像信息和標(biāo)定參照物中提取相機(jī)的內(nèi)外方位元素，從而建立空間三維坐標(biāo)系中位置、方向信息與圖像二維坐標(biāo)系間的精確映射關(guān)系。相機(jī)的內(nèi)方位元素指相機(jī)的主點(diǎn)、主距、畸變等參數(shù)，外方位元素指相機(jī)位置和方向等參數(shù)[1]。相機(jī)內(nèi)外方位元素的標(biāo)定在攝影測(cè)量、圖像校正、三維重建以及機(jī)器視覺(jué)等諸多領(lǐng)域具有重要作用。相機(jī)自標(biāo)定方法為近幾年興起的一種相機(jī)標(biāo)定技術(shù)，該方法由一臺(tái)或多臺(tái)相機(jī)從不同的位置和方向?qū)σ粋€(gè)尺寸已知的棋盤格模板拍攝多幅圖像，以模板和圖像的對(duì)應(yīng)關(guān)系求解出內(nèi)外方位元素，具有標(biāo)定操作過(guò)程簡(jiǎn)單易行、實(shí)驗(yàn)結(jié)果精度較高的優(yōu)點(diǎn)。

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

1 迭代算法

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

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

1.1 共線方程

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

式中

以上3個(gè)方程相乘后得：

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

在標(biāo)定過(guò)程中，模板中取棋盤角點(diǎn)的左上點(diǎn)為原點(diǎn)，坐標(biāo)值為實(shí)際距離值。對(duì)于理想相機(jī)模型，模板上一點(diǎn)（X,Y,Z）投影到圖像坐標(biāo)系后其坐標(biāo)為（x,y）。但由于圖像畸變，其實(shí)際坐標(biāo)為（x′,y′）。

模板坐標(biāo)系和圖像坐標(biāo)系之間的關(guān)系為

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

設(shè)棋盤模板上的棋盤角點(diǎn)數(shù)量為m，相機(jī)從不同位置拍攝的模板圖像數(shù)量為n。對(duì)于每個(gè)相機(jī)的放置位置，其對(duì)應(yīng)的圖像上能檢測(cè)到并能計(jì)算出坐標(biāo)的棋盤角點(diǎn)的數(shù)量為mi（i=1,2,3,…,n）個(gè)，有mi≤m，因此所有未知數(shù)數(shù)量為6n+9，所有共線方程的數(shù)量為，且不超過(guò)所拍攝所有圖像中的棋盤角點(diǎn)數(shù)量總數(shù)2mn，在檢驗(yàn)過(guò)程中的超定方程為

1.2 求解算法

1.2.1 迭代計(jì)算

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

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

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

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

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

1.2.3 收斂標(biāo)準(zhǔn)

文獻(xiàn)[6]中指出New ton-Raphson方法在初始估計(jì)的迭代值非常接近方程的解時(shí)，即迭代結(jié)果已滿足精度要求時(shí)收斂。本文結(jié)束迭代計(jì)算過(guò)程的收斂標(biāo)準(zhǔn)為

式中ε為一個(gè)視計(jì)算情況而定的一個(gè)無(wú)限小的數(shù)。

另一種收斂標(biāo)準(zhǔn)是檢驗(yàn)解向量 Δξk的L∞范數(shù)的變化[6]。對(duì)于上述兩種標(biāo)準(zhǔn)，都在ε= 10-8、迭代 10～20次時(shí)達(dá)到收斂。收斂結(jié)果不包含二次偏導(dǎo)項(xiàng)，與標(biāo)準(zhǔn)的New ton-Raphson算法不同，這是因?yàn)樵诓捎米钚《朔ń夥匠探M時(shí)，有些情況下雅克比矩陣可能有較大的條件數(shù)。

1.2.4 萬(wàn)向節(jié)死鎖

φ= π/2時(shí)，旋轉(zhuǎn)矩陣為

式（14）說(shuō)明，解向量?jī)H依靠ω和κ的變化，獨(dú)立未知量的數(shù)量減少了一個(gè)，這種現(xiàn)象稱為“萬(wàn)向節(jié)死鎖”。對(duì)于這種情況，檢驗(yàn)結(jié)果顯示收斂的解向量中ω和κ為隨機(jī)值，但二者的差值正確。對(duì)于φ= π/2的萬(wàn)向節(jié)死鎖情況無(wú)需特殊措施。由式（14）可以看出，相比于ω、κ各自的值，ω-κ的差值更重要，因?yàn)樽罱K要得到模板點(diǎn)坐標(biāo)和圖像點(diǎn)坐標(biāo)之間的映射關(guān)系[6-7]。

1.3 求解算法

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

1.3.1 焦距值的確定

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

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

則由雙空間幾何關(guān)系的特性可得到：

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

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

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

iiiiiii可得：

至此可得到2個(gè)主距的初始估計(jì)。

1.3.2 其他參數(shù)的確定

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

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

1.4 重投影誤差

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

從上式看出，目標(biāo)函數(shù)是由一組殘量平方和組成，若得到的目標(biāo)函數(shù)值越小，則數(shù)據(jù)擬合得越好，求得的參數(shù)值越準(zhǔn)確。本文算法不是直接求解非線性的畸變模型，而是將其轉(zhuǎn)化成求解非線性最小二乘問(wèn)題，通過(guò)非線性優(yōu)化算法多次迭代，最后得到使目標(biāo)函數(shù)值最小的參數(shù)值，降低了求解難度。

2 與張氏標(biāo)定法的比較

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

1）分別以圖像幅數(shù)范圍、網(wǎng)格板角度范圍等參數(shù)拍攝照片；

2）對(duì)所拍攝照片進(jìn)行圖像處理，用Harris角點(diǎn)提取算法提取角點(diǎn)；

3）利用棋盤格標(biāo)定板在世界坐標(biāo)系的坐標(biāo)和相應(yīng)的圖像角點(diǎn)坐標(biāo)計(jì)算每幅圖像的單應(yīng)性矩陣；

4）由所得單應(yīng)性矩陣，根據(jù)自標(biāo)定算法求解出內(nèi)方位元素矩陣初值；

5）設(shè)計(jì)畸變模型；

6）將由內(nèi)方位矩陣和畸變模型計(jì)算所得的角點(diǎn)坐標(biāo)和實(shí)際圖像角點(diǎn)坐標(biāo)代入數(shù)據(jù)擬合函數(shù)進(jìn)行擬合，并通過(guò)非線性優(yōu)化算法進(jìn)行迭代，得到內(nèi)方位矩陣和畸變系數(shù)的精確值。

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

1）確定內(nèi)方位元素矩陣初值的算法；

2）畸變模型；

3）非線性優(yōu)化迭代算法。

兩種算法的異同點(diǎn)見表1。

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

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

3 實(shí)驗(yàn)

對(duì)于本文提出的標(biāo)定方法，利用文獻(xiàn)[2]提供的測(cè)試數(shù)據(jù)進(jìn)行了實(shí)驗(yàn)，如圖3所示，并將測(cè)試結(jié)果進(jìn)行對(duì)比，見表2。張氏標(biāo)定法基于Levenberg-Marquardt算法對(duì)得到的內(nèi)方位元素初值進(jìn)行優(yōu)化迭代求解，并用相同的數(shù)據(jù)目標(biāo)擬合函數(shù)進(jìn)行重投影誤差計(jì)算，作為評(píng)價(jià)標(biāo)定結(jié)果精度的標(biāo)準(zhǔn)。

圖3 張氏標(biāo)定法使用的標(biāo)定圖像Fig.3 Calibration Images of Zhangmethod

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

由表2可以看出，雖然兩種算法用于優(yōu)化的初值相差較多，但經(jīng)過(guò)優(yōu)化后得到的終值十分相近。但張氏標(biāo)定法得到的重投影誤差為0.335個(gè)像元[14]，本文算法得到的重投影誤差為0.265個(gè)像元，這表示本文方法計(jì)算出的相機(jī)內(nèi)方位元素精度更高，更加接近實(shí)際相機(jī)成像模型。由于畸變模型不同，鏡頭畸變系數(shù)不具有可比性，但是本文的畸變模型有5項(xiàng)參數(shù)，能夠更加準(zhǔn)確反應(yīng)實(shí)際鏡頭的畸變情況。

4 結(jié)束語(yǔ)

本文介紹了一種利用相機(jī)在不同位置拍攝的多幅圖像中的特征點(diǎn)信息，計(jì)算相機(jī)內(nèi)外方位元素的迭代計(jì)算方法。該方法用B雙空間幾何的特性，利用消失點(diǎn)、消失線理論得到迭代方法的初值。在無(wú)需進(jìn)行干預(yù)初始估計(jì)值的情況下，通過(guò)迭代計(jì)算得到精度較高的內(nèi)方位元素值。通過(guò)與張氏標(biāo)定法的對(duì)比實(shí)驗(yàn)來(lái)比較兩標(biāo)定法的精度，實(shí)驗(yàn)結(jié)果表明，本文介紹的算法誤差更小，其標(biāo)定的相機(jī)參數(shù)也更加精確，能夠用于相機(jī)內(nèi)方位元素的標(biāo)定，且標(biāo)定精度較理想。

References)

[1] 譚啟蒙,胡成威,高升.空間機(jī)械臂視覺(jué)相機(jī)內(nèi)參標(biāo)定技術(shù)研究[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雙空間幾何中基于消隱點(diǎn)的攝像機(jī)標(biāo)定[J].計(jì)算機(jī)工程與應(yīng)用,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] 熊磊,劉國(guó)棟,劉小勇,等.數(shù)字圖像相關(guān)中基于非線性灰度改變模型的灰度梯度迭代算法[J].機(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] 吳聞.游戲程序設(shè)計(jì)中若干問(wèn)題的探討[J].電腦知識(shí)與技術(shù),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] 梅雪,夏良正,李久賢,等.一種三維場(chǎng)景的消失點(diǎn)檢測(cè)算法[J].信號(hào)處理,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.