顧及物方坐標(biāo)誤差的張正友標(biāo)定算法

瞿 鍇,郭中佑,唐育林,陳抒錄,陳朋明

(1.中機(jī)國際工程設(shè)計研究院有限責(zé)任公司,長沙 410021;2.湘潭大學(xué) 土木工程學(xué)院，湖南 湘潭 411100;3.湘潭市勘測設(shè)計院，湖南 湘潭 411100)

0 前 言

近年來，無人機(jī)攝影測量技術(shù)得到了迅猛發(fā)展，在生產(chǎn)實(shí)踐中廣泛應(yīng)用于農(nóng)業(yè)、林業(yè)、海洋等領(lǐng)域[1]。在無人機(jī)攝影測量中普遍采用非量測型數(shù)碼相機(jī)，但其存在著相機(jī)光學(xué)畸變大的缺點(diǎn)，故需要對相機(jī)成像參數(shù)予以檢校。相機(jī)檢校是通過對已知坐標(biāo)或相對位置關(guān)系明確的物體進(jìn)行拍照，利用物體相片量測坐標(biāo)與已知坐標(biāo)之間的對應(yīng)關(guān)系來求取相機(jī)成像參數(shù)。為保證參數(shù)解算的可靠性，通常需要進(jìn)行大量的多余觀測，這些數(shù)據(jù)被認(rèn)為服從高斯-馬爾可夫模型(Gauss-Markov Theory)，可構(gòu)建眾多的觀測數(shù)據(jù)方程，在最小二乘準(zhǔn)則約束下求解方程獲得參數(shù)的最優(yōu)估值。相機(jī)檢校常用的Tsai算法、張正友標(biāo)定算法、滅點(diǎn)法和直接線性變換(Direct Line transform，DLT)算法在計算相機(jī)參數(shù)時都采用了G-M模型[2]。但G-M模型只考慮了觀測向量包含誤差，則將導(dǎo)致應(yīng)用傳統(tǒng)最小二乘準(zhǔn)則估計的參數(shù)不是最優(yōu)解。

考慮到G-M模型的缺點(diǎn)，整體最小二乘準(zhǔn)則被提出并逐漸應(yīng)用到上所述相機(jī)檢校工作中所列方程求解的約束條件[3]。如陶葉青[4]、楊會軍[5]、孔建[6]、魯鐵定[7]和孫培芪等[8]將整體最小二乘應(yīng)用到解決方案中去，實(shí)驗(yàn)結(jié)果表明，整體最小二乘在擬合函數(shù)、轉(zhuǎn)換參數(shù)的求解問題中更加具備優(yōu)勢。其次是檢校場類型對相機(jī)檢校起著很重要的作用，倪靜、曹良忠等[9-10]研究了相對條件較為苛刻的野外檢校場、田雷等[11]應(yīng)用了室內(nèi)檢校場檢校、李平[12]簡化了三維空間檢校場，應(yīng)用二維檢校場進(jìn)行相機(jī)檢校中比較出名的是張正友檢校法[13]，這種方法只需要一個規(guī)則的格網(wǎng)狀平板，通過移動相機(jī)從不同位置對其拍照，即可計算出內(nèi)外方位元素。

考慮到野外檢校場和室內(nèi)三維檢校場建立與維護(hù)的不易，而無人機(jī)攝影測量相機(jī)又需要頻繁檢校，本文采用自制的格網(wǎng)型模板，先以尼康D800為例，利用控制變量法分析相片數(shù)量和拍攝角度對參數(shù)解算精度的影響，確定相片數(shù)量和拍攝方式后，使用Nikon D800相機(jī)與索尼ILCE-7R相機(jī)拍攝模板獲取數(shù)據(jù)。接著，考慮到整體最小二乘算法較最小二乘算法的模型更為完善，顧及到了模型中系數(shù)矩陣可能存在的誤差，使用整體最小二乘算法對張正友檢校算法予以改進(jìn)，并通過試驗(yàn)進(jìn)行驗(yàn)證，以提高非量測相機(jī)的檢校精度，保證測繪生產(chǎn)獲得較高質(zhì)量的成果。

1 研究方法

自檢校方法是Maybank[14]提出的一種與相機(jī)移動和空間幾何信息無關(guān)的相機(jī)檢校方法。這方面國內(nèi)的工作有孟曉橋等[15]從多個方向拍攝模板圖像，應(yīng)用圓檢校算法，利用射影的不變性，解算出了相機(jī)參數(shù)。從計算原理講，張正友檢校法也可以歸于自檢校方法。相機(jī)的自檢校方法是指不需要依據(jù)一個已經(jīng)建立好的檢校場，是通過影像與影像之間的對應(yīng)關(guān)系，計算出相機(jī)的參數(shù)。此種方法的優(yōu)點(diǎn)是價格低廉，可以隨時隨地對相機(jī)進(jìn)行檢校；缺點(diǎn)則是精度會比較低。張正友相機(jī)標(biāo)定算法是將打印的棋盤格貼在平面上，從不同角度進(jìn)行拍攝，利用像點(diǎn)和物點(diǎn)的幾何關(guān)系確定相機(jī)內(nèi)外方位元素的初值，然后將其作為非線性相機(jī)幾何模型的初值，求出單應(yīng)性矩陣，以此解算出攝像機(jī)內(nèi)外參數(shù)和畸變系數(shù)，模型如圖1所示。

圖1 張正友算法模型

1.1 單應(yīng)性矩陣

設(shè)世界坐標(biāo)某點(diǎn)的物理坐標(biāo)為M=[U,V,W,1]T，其對應(yīng)的像素坐標(biāo)系下的像素坐標(biāo)點(diǎn)為m=[u,v,1]T，故標(biāo)定用的棋盤格平面到圖像平面的單應(yīng)性關(guān)系為：

(1)

其中:s為尺度因子;A為相機(jī)內(nèi)參矩陣;包括仿射變換和透視投影;(u0,v0)為像主點(diǎn)坐標(biāo);α,β為焦距與像素橫縱比的融合;γ為徑向畸變參數(shù);R為旋轉(zhuǎn)矩陣;t為平移向量;[Rt]也被成為相機(jī)的外參矩陣。

設(shè)棋盤格平面W=0，可得：

(2)

A[r1,r2,t]為單應(yīng)性矩陣，即內(nèi)參矩陣和外參矩陣的積，記作H=[r1,r2,t]，得：

(3)

所以每一格標(biāo)定板的角點(diǎn)可以提供兩個坐標(biāo)約束方程：

(4)

上式對于同一張圖片上所有的角點(diǎn)均成立，其中[u,v]是像素坐標(biāo)系下標(biāo)定板的角點(diǎn)坐標(biāo)，[U,V]是世界坐標(biāo)系下標(biāo)定板的角點(diǎn)坐標(biāo)。

1.2 內(nèi)參矩陣A和外參矩陣的計算

(5)

化簡得：

(6)

公式(6)中：v是僅包含H元素的矩陣，b=[B11,B12,B22,B23,B33]T。由于矩陣H已知，矩陣v又全部由矩陣H的元素構(gòu)成，因此矩陣v已知。每張標(biāo)定板圖片可以提供一個約束關(guān)系，采用最小二乘擬合最佳的向量b，得到矩陣B，所以有：

(7)

相機(jī)的內(nèi)參矩陣A取決于相機(jī)的內(nèi)部參數(shù)，無論標(biāo)定板和相機(jī)的位置關(guān)系怎么變換，相機(jī)的內(nèi)參矩陣是不變的。而外參矩陣反映的是標(biāo)定板和相機(jī)的位置關(guān)系，對于被拍攝的不同圖片，標(biāo)定板和相機(jī)的位置關(guān)系已經(jīng)改變，此時每一張圖片對應(yīng)的外參矩陣都是不相等的。

因?yàn)镠=A[r1,r2,t]，求解得到H和A后，可通過式(8)得到外參矩陣。

[r1,r2,t]=A-1H

(8)

1.3 畸變系數(shù)k的計算

對于畸變系數(shù)k，張氏標(biāo)定法只關(guān)注了影響最大的徑向畸變，數(shù)學(xué)表達(dá)式為[18]：

(9)

k=[k1,k2]T=(DTD)-1DTd

(10)

1.4 整體最小二乘原理

在實(shí)際工作中，系數(shù)矩陣中也會包含誤差[7-8]，故真誤差的數(shù)學(xué)期望往往不為0。求解Ax=b時，應(yīng)該同時考慮C和c二者的誤差和擾動，令B矩陣的誤差擾動為E，向量b的誤差向量為e，矩陣方程為：

(C+E)x=c+e

(11)

可得：

c+e-(xT×I)E-Cx=0

(12)

式中:xT×I表示兩個矩陣做克羅內(nèi)克積，令F=[i,-(xT×I)],v=[e,E]T，可得：

Fv+c-Cx=0

(13)

整體最小二乘準(zhǔn)則為：

minφ=vTPv-2λT(Fv+c-Cx)

(14)

式中:P為權(quán)矩陣;λ為Largrange乘數(shù)，分別進(jìn)行偏導(dǎo)：

(15)

可得：

x=[(C+I)T(FP-1F)-1C]-1(C+I)T(FP-1F)-1c

(16)

對公式(6)和公式(10)采用整體最小二乘計算，得到最優(yōu)估計值。

2 檢校場與計算軟件的設(shè)計

2.1 室內(nèi)二維檢校場的實(shí)現(xiàn)

這里選用尼康D800和索尼ILCE-7R兩種相機(jī)獲取影像數(shù)據(jù)，拍攝時設(shè)置為手動曝光方式，物鏡對焦于無窮遠(yuǎn)，相機(jī)傳感器參數(shù)見表1。

表1 尼康D800和索尼ILCE-7R相機(jī)傳感器參數(shù)

計算機(jī)型號為聯(lián)想Lenovo G50-70 20351筆記本電腦，CPU酷睿i5-4258U，內(nèi)存DDR3L 1600 MHz。實(shí)驗(yàn)選取棋盤格作為相機(jī)檢校的模板，二維控制場采用基于純平液晶顯示器的檢校方法[19]，使所有控制點(diǎn)在一個平面上。

為研究拍攝角度對標(biāo)定精度的影響，從不同角度進(jìn)行拍攝，得到圖像，棋盤格大小為9×5個方格，每個方格尺寸為45 mm × 45 mm(長×寬)，角點(diǎn)數(shù)量為32個，部分影像如圖2所示。

圖2 不同方向拍攝的棋盤格控制場

2.2 相機(jī)檢校軟件的設(shè)計

這里根據(jù)計算原理利用Visual Studio平臺和C#設(shè)計相機(jī)檢校軟件，主要分為如下幾項功能：文件讀取數(shù)據(jù)、最小二乘約束、整體最小二乘約束、畸變改正、查看報告和退出等功能。程序操作非常簡便，首先點(diǎn)擊可執(zhí)行文件(.exe)，啟動相機(jī)檢校程序，在相同目錄下添加不同角度拍攝的棋盤格影像，可以為tif、jpg、bmp等格式，使用對應(yīng)算法自動提取棋盤格角點(diǎn)的坐標(biāo)作為后續(xù)計算的數(shù)據(jù)，程序解算會自動將最小二乘約束與整體最小二乘約束解算所得到的參數(shù)值作為初始值使用，直接點(diǎn)擊計算輸出按鈕可得到最終的相機(jī)參數(shù)報告，包括最小二乘約束下和整體最小二乘約束下的不同相機(jī)檢校算法的相機(jī)參數(shù)。參數(shù)解算部分結(jié)果如圖3所示。前面部分為相機(jī)的檢校參數(shù)值，后面部分為各參數(shù)的解算誤差。

圖3 相機(jī)檢校軟件

3 實(shí)驗(yàn)結(jié)果與分析

3.1 相片數(shù)量對檢校的影響

為探究相片數(shù)量對相機(jī)檢校精度的影響，依次使用3、5、7、9、11不同數(shù)量的相片參與計算，結(jié)果如表2所示，參與解算相片數(shù)量的不同會導(dǎo)致檢校參數(shù)的精度有所差異。相機(jī)參數(shù)包括：f(焦距)，x0、y0(主點(diǎn)坐標(biāo))和k1、k2(畸變參數(shù))。y值在3張相片解算時最優(yōu)，后面呈波動分布。f值在9張相片解算時達(dá)到最優(yōu)，11張相片為35.85445，相比而言，值下降了0.066。x0、k1和k2的值隨著參加解算相片的數(shù)量越多，精度越高，分別為3661.5、-0.0114和0.0403。3張相片進(jìn)行解算時，f、k1和k2與參考值相差最大，這是因?yàn)橄嗥臄?shù)量太少導(dǎo)致不能夠充分反應(yīng)出鏡頭的畸變情況。隨著后面參與解算的相片數(shù)量越來越多，值變化的越來越趨于平穩(wěn)狀態(tài)，說明相片并不是越多越好，但也不能太少，過少會使得結(jié)果精度不高，過多則會導(dǎo)致計算量的增加，從而導(dǎo)致效率降低。

表2 尼康D800相片數(shù)量與檢校參數(shù)計算結(jié)果

3.2 拍攝角度對檢校的影響

為探究拍攝相片的角度對檢校結(jié)果的影響，采用兩種方式拍攝相片，先拍攝一組近似平行的照片計算檢校參數(shù)，再以不同的角度拍攝相片計算檢校參數(shù)，分析得到的檢校結(jié)果，兩組相片的檢校結(jié)果如表3所示。當(dāng)相機(jī)位置處于平行狀態(tài)時，f和x0的精度較高，分別為35.96734和3654.2；相片處于不同的角度時，y0、k1和k2精度較高，依次為2404.7、-0.006、0.019。一般認(rèn)為，近似平行的相片計算的檢校參數(shù)精度會低于不同角度拍攝的照片計算的檢校參數(shù)精度。因?yàn)閺堈褭z校方法是通過不同相片之間的關(guān)系計算的，所以拍攝的相片近似平行的時候，有的參數(shù)之間會存在很強(qiáng)的相關(guān)性，就會影響最后的結(jié)果。在使用張正友標(biāo)定法對相機(jī)進(jìn)行檢校時，應(yīng)該以不同的角度對棋盤格進(jìn)行拍攝，從而保證高精度的校驗(yàn)結(jié)果。y0、k1和k2的精度與預(yù)期表現(xiàn)相符，而f和x0的結(jié)果精度恰恰相反，可能與相片拍攝的質(zhì)量有關(guān)。

表3 尼康D800相機(jī)拍攝角度與檢校參數(shù)計算結(jié)果

3.3 改進(jìn)前后方法的對比

由以上可知，相機(jī)拍攝相片時應(yīng)盡量從不同方向進(jìn)行拍攝，數(shù)量根據(jù)區(qū)域大小確定，本次使用11張不同方向拍攝的相片參與解算。為檢驗(yàn)改進(jìn)前后張正友方法在尼康D800和索尼ILCE-7R兩種相機(jī)中的情況，分別使用兩款相機(jī)隨機(jī)在5個攝站對棋盤格進(jìn)行拍攝，每張相片不重復(fù)，角點(diǎn)數(shù)量為45個，部分圖片如圖2所示。分別采用經(jīng)典最小二乘算法與整體最小二乘算法求解單應(yīng)性矩陣參數(shù)，進(jìn)而解算相機(jī)參數(shù)，并將結(jié)果進(jìn)行對比分析。表4為尼康D800和索尼ILCE-7R相機(jī)內(nèi)定向參數(shù)計算結(jié)果，在尼康D800數(shù)據(jù)集上，LS和TLS法得到的焦距f與參考值分別相差0.421 mm、0.259 mm；主點(diǎn)坐標(biāo)(x0，y0)分別相差176.117像素、44.159像素；畸變參數(shù)k1/k2分別相差0.037/0.012、0.466/0.150；計算時間相差0.345。在索尼ILCE-7R數(shù)據(jù)集上，LS和TLS法得到的焦距f與參考值分別相差0.435 mm、0.120 mm；主點(diǎn)坐標(biāo)(x0，y0)分別相差155.731像素、93.362像素；k1/k2分別相差0.1225/0.110、0.150/0.128；計算時間相差0.185。相比LS，TLS計算得到的參數(shù)精度均得到了不同程度的提升，說明TLS確實(shí)能夠較好的考慮觀測值和系數(shù)存在的誤差，解算出更為精確的結(jié)果。雖然相機(jī)參數(shù)精度均有提升，但增加了計算時間，這是必然的結(jié)果，因?yàn)橛嬎懔康脑黾?，勢必會?dǎo)致計算效率的降低?？偟膩碚f，整體最小二乘方法相比傳統(tǒng)的最小二乘方法在張正友標(biāo)定算法相機(jī)參數(shù)解算中表現(xiàn)更優(yōu)。

4 結(jié) 論

本文在傳統(tǒng)張正友的算法上進(jìn)行改進(jìn)，考慮顧及物方控制點(diǎn)坐標(biāo)誤差，用整體最小二乘約束求解單應(yīng)性矩陣，將其精度得到提高的結(jié)果參與后面相機(jī)參數(shù)的運(yùn)算。為確定最佳解算相片數(shù)據(jù)，先以尼康D800為例，使用控制變量法探究拍攝角度和相片數(shù)量對解算結(jié)果的影像。然后采用改進(jìn)前后的張正友標(biāo)定算法在尼康D800和索尼ILCE-7R相機(jī)開展實(shí)驗(yàn)，結(jié)果表明：

(1)解算參數(shù)x0、k1和k2的值隨著參加解算相片的數(shù)量越多，精度越高，值逐漸趨于平穩(wěn)，說明盲目增加相片數(shù)量不會一直提高解算精度，反而會增大計算量。

(2)以不同角度拍攝的照片計算的檢校參數(shù)y0、k1和k2精度高于近似平行的相片計算的檢校參數(shù)精度，不同角度的相片能夠更好提供約束條件，得到更優(yōu)解。

(3)在尼康D800和索尼ILCE-7R數(shù)據(jù)上，相比傳統(tǒng)張正友算法，f(焦距)、x0、y0(主點(diǎn)坐標(biāo))和k1、k2(畸變參數(shù))的解算精度均得到了提升，但計算量提升了一倍，這在預(yù)期之內(nèi)。

本次實(shí)驗(yàn)的某些相機(jī)參數(shù)精度與預(yù)期不相符，可能是因?yàn)橄嗥|(zhì)量沒有做抗差估計，在進(jìn)行解算時無法保證每張相片的質(zhì)量統(tǒng)一，從而導(dǎo)致結(jié)果出現(xiàn)差異性。后面可以針對不同相機(jī)檢校算法驗(yàn)證本文所提方法的可行性，并為了保證數(shù)據(jù)質(zhì)量統(tǒng)一標(biāo)準(zhǔn)，對拍攝相片進(jìn)行抗差處理。