顧及物方坐標誤差的張正友標定算法

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

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

0 前 言

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

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

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

1 研究方法

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

圖1 張正友算法模型

1.1 單應性矩陣

設世界坐標某點的物理坐標為M=[U,V,W,1]T，其對應的像素坐標系下的像素坐標點為m=[u,v,1]T，故標定用的棋盤格平面到圖像平面的單應性關系為：

(1)

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

設棋盤格平面W=0，可得：

(2)

A[r1,r2,t]為單應性矩陣，即內參矩陣和外參矩陣的積，記作H=[r1,r2,t]，得：

(3)

所以每一格標定板的角點可以提供兩個坐標約束方程：

(4)

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

1.2 內參矩陣A和外參矩陣的計算

(5)

化簡得：

(6)

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

(7)

相機的內參矩陣A取決于相機的內部參數，無論標定板和相機的位置關系怎么變換，相機的內參矩陣是不變的。而外參矩陣反映的是標定板和相機的位置關系，對于被拍攝的不同圖片，標定板和相機的位置關系已經改變，此時每一張圖片對應的外參矩陣都是不相等的。

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

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

(8)

1.3 畸變系數k的計算

對于畸變系數k，張氏標定法只關注了影響最大的徑向畸變，數學表達式為[18]：

(9)

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

(10)

1.4 整體最小二乘原理

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

(C+E)x=c+e

(11)

可得：

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

(12)

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

Fv+c-Cx=0

(13)

整體最小二乘準則為：

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

(14)

式中:P為權矩陣;λ為Largrange乘數，分別進行偏導：

(15)

可得：

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

(16)

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

2 檢校場與計算軟件的設計

2.1 室內二維檢校場的實現

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

表1 尼康D800和索尼ILCE-7R相機傳感器參數

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

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

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

2.2 相機檢校軟件的設計

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

圖3 相機檢校軟件

3 實驗結果與分析

3.1 相片數量對檢校的影響

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

表2 尼康D800相片數量與檢校參數計算結果

3.2 拍攝角度對檢校的影響

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

表3 尼康D800相機拍攝角度與檢校參數計算結果

3.3 改進前后方法的對比

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

4 結 論

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

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

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

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

本次實驗的某些相機參數精度與預期不相符，可能是因為相片質量沒有做抗差估計，在進行解算時無法保證每張相片的質量統一，從而導致結果出現差異性。后面可以針對不同相機檢校算法驗證本文所提方法的可行性，并為了保證數據質量統一標準，對拍攝相片進行抗差處理。