無物方信息的未檢校近景影像的方位元素估計

周擁軍，寇新建，朱建軍

（1.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海200240；2.中南大學(xué) 測繪與國土信息工程系，湖南 長沙410083）

內(nèi)外方位元素的估計問題是近景攝影解析處理的核心問題，近景攝影測量的短攝影基線、構(gòu)形不規(guī)則及大角度攝影、多重疊等特點增加了定向的難度.對于量測相機(jī)或經(jīng)預(yù)檢校的普通相機(jī)，相機(jī)的內(nèi)方位元素已知，若有一定數(shù)量的控制點，可以通過直接線性變換、空間后方交會等傳統(tǒng)解析方法確定外方位元素.對于無物方控制點的立體影像序列則只能實現(xiàn)度量重建［1］，還需要通過絕對定向才能得到物方坐標(biāo).未檢校相機(jī)的內(nèi)方位元素未知，若采用傳統(tǒng)方法，需要通過自檢校［2］或廣義相對定向［3］恢復(fù)內(nèi)方位元素，然后經(jīng)相對定向得到獨(dú)立模型，通過模型連接得到統(tǒng)一模型，最后經(jīng)絕對定向獲取物方坐標(biāo).由于自檢校得到的內(nèi)方位元素精度低，加之多個獨(dú)立模型的連接誤差累積等影響，導(dǎo)致定向結(jié)果誤差較大.

本文采用改進(jìn)的層次重建方法，攝影過程中保持內(nèi)方位元素不變，選用其中的1張影像的像空間坐標(biāo)作為參考坐標(biāo)系，通過3張以上影像的同名像點恢復(fù)投影矩陣和結(jié)構(gòu)，然后采用基于絕對二次曲面的自檢校方法恢復(fù)相機(jī)內(nèi)方位元素，從而在無物方法得到物方坐標(biāo)，該方法既可應(yīng)用于低精度的三維重建，也可以作為高精度光束法平差的初值，是一種區(qū)別于傳統(tǒng)定向方法的新方法.但由于無物方控制和不考慮畸變，因此定向的精度相對較低，只能用于低精度的測量或作為光束法平差的初值。

1 層次重建原理

現(xiàn)實世界的空間關(guān)系用歐氏空間來描述，若不考慮鏡頭畸變和其他誤差的影響，相機(jī)成像模型滿足射影變換關(guān)系，攝影測量中將成像模型表達(dá)為共線方程，若寫成齊次坐標(biāo)的形式，與視覺領(lǐng)域的針孔（Pinhole）成像模型是一致的，為表達(dá)方便，本文均用齊次坐標(biāo)的形式表示.

設(shè)任一空間點的物方坐標(biāo)為＝表示其投影得到的像點坐標(biāo)，不考慮相機(jī)的畸變，則投影模型可表示為［1］

式中：A為相機(jī)的內(nèi)部矩陣，若不考慮傾斜因子，共由4個獨(dú)立參數(shù)構(gòu)成，包括x、y方向上的比例因子αx，αy和像主點坐標(biāo)u0，v0；Ri表示旋轉(zhuǎn)矩陣；ti＝（tx tytz）T表示投影中心的物方坐標(biāo)，Ri、ti統(tǒng)稱為外方位元素，由6個獨(dú)立參數(shù)組成；P為3×4的投影矩陣.

計算機(jī)視覺領(lǐng)域?qū)缀巫儞Q分4個層次，即射影（Projective）變換、仿射 （Affine）變換、度量（Metric）變換和歐氏（Euclidean）變換，分別對應(yīng)15，12，7，6個獨(dú)立變換參數(shù)［4－5］.設(shè)某空間點變換前后的齊次坐標(biāo)分別表示為變換關(guān)系表示為變換矩陣T為4×4的非奇異矩陣，對于不同層次的變換，T應(yīng)滿足不同的約束條件.層次重建的思路就是先利用立體影像的同名點進(jìn)行射影重建，然后利用立體影像間的約束關(guān)系逐層過渡到歐氏重建.即先確定一個投影矩陣P和射影結(jié)構(gòu)，使其滿足很顯然射影重建的結(jié)果不唯一，因為對于任一非奇異4×4的矩陣Q＝λPQQ－1也成立，層次重建的思想就是找到變換矩陣T，使其滿足PT＝A［R｜－Rt］，其中，~XE為歐氏結(jié)構(gòu)，R為旋轉(zhuǎn)矩陣.近年來國內(nèi)外學(xué)者采用改進(jìn)的層次重建方法實現(xiàn)歐氏重建，如以平行透視或準(zhǔn)透視代替通用的透視變換［6－7］，或考慮內(nèi)方位元素改變的情況等［8］.

2 基于分解算法的射影重建

由投影關(guān)系式可知影像坐標(biāo)受兩種因素的影響，其一是投影矩陣，也稱為相機(jī)的運(yùn)動（motion）參數(shù)，其二是坐標(biāo)，即結(jié)構(gòu)（structure），從序列圖像中恢復(fù)結(jié)構(gòu)和運(yùn)動是計算機(jī)視覺的一個基本問題.圖像坐標(biāo)是相機(jī)姿態(tài)和結(jié)構(gòu)的雙線性組合，如果其中一個因素保持不變，則像點坐標(biāo)是另一因素的線性形式.分解方法就是基于這種考慮而得到的一種簡潔高效的算法，它通過矩陣的奇異值分解（SVD）方法將圖像測量矩陣分解成相機(jī)運(yùn)動和三維結(jié)構(gòu)分量.分解算法最早是由 Tomasi等［9］提出，Sturm，Heyden等改進(jìn)了該方法［10－11］，并用于從序列圖像中恢復(fù)運(yùn)動和結(jié)構(gòu).

設(shè)有n個空間點，每個點在m幅圖像中均有成像表示點i在第j幅影像的齊次坐標(biāo)表示點i的射影結(jié)構(gòu)，λij表示投影深度，則存在投影關(guān)系：

將所有像點組成一個聯(lián)合矩陣，令P，X分別表示所有的投影矩陣和結(jié)構(gòu)，則有：

由式（3）可以看出，矩陣Ws可以表示成3m×4的矩陣P和4×n的矩陣X的線性組合，因此Ws的秩應(yīng)滿足rank（Ws）≤4，當(dāng)空間點的分布不存在奇異構(gòu)形（共面或共線）的情況下，投影矩陣的秩應(yīng)等于4.由于測量誤差的影響使得Ws的秩大于4，需要找到一個秩為4的矩陣＾Ws逼近Ws，這可通過Ws的SVD分解實現(xiàn)，這就是基于分解算法的射影重建的思想，將Ws進(jìn)行奇異值分解：

其中，k＝min（3m，n），σ1≥σ2≥…≥σk，U，V均為經(jīng)svd分解后的矩陣，設(shè)U′和V′T分別表示U、V的前4列向量，則：

基于分解算法的射影重建過程歸納為以下步驟：① 將坐標(biāo)歸一化，即對所有像點uij（i＝1，…，m，j＝1，…，n）作相似變換或按文獻(xiàn)［6］的規(guī)則選取歸一化矩陣N；② 設(shè)所有尺度參數(shù)初值組成測量矩陣W；③ 將W進(jìn)行奇異值分ss解，得到＾Ws；④ 根據(jù)計算尺度參數(shù)，并組成新的測量矩陣判斷‖Ws－‖是否小于給定的限差，若大于限差則轉(zhuǎn)入③繼續(xù)計算，否則計算結(jié)束；⑤ 將投影矩陣進(jìn)行逆變換得到N－1P.

3 方位元素估計

3.1 內(nèi)方位元素估計

對于無任何物方控制的立體影像，需要通過自檢?；謴?fù)內(nèi)方位元素，視覺領(lǐng)域的自檢校方法常采用基于絕對二次曲線和絕對二次曲面的自檢校方法［5］，這里采用基于絕對二次曲面的自檢校方法.射影重建恢復(fù)了射影結(jié)構(gòu)，還需要進(jìn)一步變換得到歐氏結(jié)構(gòu)，若用表示第i張影像的歐氏投影矩陣，則其必須滿足：

若存在一個非奇異的4×4的變換矩陣TPE可以將式（2）的射影重建轉(zhuǎn)變?yōu)闅W氏重建，則有：

式中，Ri為第i張影像的旋轉(zhuǎn)矩陣.

取物方坐標(biāo)為第一張影像的像空間坐標(biāo)系，則有R1＝I，t1＝0，代入式（9），則TPE必須滿足：

d＝（abc）T是一未知參數(shù)向量，其物理意義表示無窮遠(yuǎn)平面的法向量［2］.從射影變換到歐氏變換的工作就是利用已有的射影矩陣Pi，確定參數(shù)矩陣A、d，使所有影像均滿足式（9）.

用T3表示TPE的前三列向量，pij為投影矩陣Pi的元素，代入式（9），有：

考慮旋轉(zhuǎn)矩陣的正交性，將式（11）兩邊取轉(zhuǎn)置后相乘：

令，則式（12）可表示為

其中：

由于矩陣的對稱性，每張影像按式（13）可組成6個方程，對于m個影像，未知參數(shù)為αx，αy，u0，v0，a，b，c，λi（i＝1，…，m），共有m＋7個，由于P1已經(jīng)滿足了式（13），因此至少需要已知兩張像片的投影矩陣才能解出所有的未知參數(shù)，即共需要3張以上影像.

解算時先用線性求解參數(shù)初值，然后用最小二乘法求精，對于普通數(shù)碼相機(jī)，像主點的坐標(biāo)可近似選擇影像中心，并將影像坐標(biāo)平移至像主點位置，此時u0＝v0＝0，而絕大多數(shù)數(shù)碼相機(jī)水平和垂直方向的比例系數(shù)近似相等，令αx＝αy＝α，，代入式（13），得：

不考慮各參數(shù)間的相關(guān)性，令Θ＝（α2aαbα c‖d‖21），將式（14）寫成齊次方程的形式：

若共選用m張影像，則B為4（m－1）×6的矩陣，當(dāng)m≥3時，可以在‖BΘ‖＝min的最小二乘準(zhǔn)則下求解，此時Θ的值對應(yīng)于矩陣B的最小特征值對應(yīng)的特征向量.

以上計算沒有考慮參數(shù)間的相關(guān)性，若要提高估計精度，需要用最小二乘法求精，可令目標(biāo)函數(shù)F最?。?/p>

其中‖‖F(xiàn)表示矩陣的Frobenius范數(shù)，式（16）可以在已知參數(shù)初值的前提下采用最小二乘迭代方法，解算出內(nèi)方位元素以及附加參數(shù)向量d，并進(jìn)一步得到變換矩陣TPE.

3.2 外方位元素估計

通過上面的計算，恢復(fù)了內(nèi)方位元素A和變換矩陣TPE，第1張影像的外方位元素已知，還需要恢復(fù)其他影像的外方位元素，對于第i（i＞1）張影像，令Q取PiTPE前三列向量，由式（9）得到：

但射影重建過程中未考慮旋轉(zhuǎn)矩陣的特性，并不能保證旋轉(zhuǎn)矩陣滿足det（Ri）＝1以及＝I的條件，可以先通過det（Ri）＝1的約束條件求出λi.

式（17）得到的Ri還未滿足正交條件，需找正交矩陣＾Ri使其最接近Ri，為此將Ri進(jìn)行SVD分解，設(shè)svd（Ri）＝UΣVT，則滿足最小二乘‖＾Ri－Ri‖2F為最小的旋轉(zhuǎn)矩陣為［12］：

4 算例

算例采用某近景攝影測量工程數(shù)據(jù)，該影像采用分辨率為1536×1024像素的Kodak DCS 420相機(jī)拍攝，拍攝過程中保持內(nèi)方位元素不變，共拍攝18張照片，經(jīng)光束法平差得到相機(jī)的內(nèi)方位元素為：f＝1703.489，u0＝764.821，v0＝509.368，k1＝－3.79×10－8，k2＝2.60×10－14.控制點坐標(biāo)用全站儀精確測量，測量精度為±1mm，選擇其中的3張照片，經(jīng)量測得到同名像點坐標(biāo)（表1），量測精度為±0.5pixel.計算時假設(shè)主點位于影像中心，取η＝1／1000將影像坐標(biāo)歸一化，然后按分解算法進(jìn)行射影重建，得到第2、3張影像的投影矩陣為

根據(jù)上述投影矩陣采用本文方法得到的相機(jī)內(nèi)方位元素見表2，從表中可以看出內(nèi)方位元素和采用傳統(tǒng)的光束法平差結(jié)果極為接近，表明基于絕對二次曲面的自檢校法能達(dá)到較高的精度.由于無物方控制，攝影基線相差一個尺度參數(shù)，將基線的X分量設(shè)為1，將得到的外方位元素和影像1－2和1－3相對定向參數(shù)比較，得到表3，可以看出差值不大，可用于確定定向參數(shù)的初值或用于低精度的近景攝影測量之中.需要指出的是，本文只選用3張影像，參與計算的像點數(shù)量和覆蓋度也不太理想，若采用多張點位構(gòu)形好、重疊度大的影像將有助于提高定向精度.

表1 像點的量測坐標(biāo)及物方坐標(biāo)Tab.1 The image coordinates and the 3Dcoordinates

表2 內(nèi)方位元素及附加參數(shù)估計值Tab.2 Interior parameters and the additional unknowns estimation

表3 外方位元素值Tab.3 Exterior parameters estimation

5 結(jié)論和建議

本文采用了一種基于層次重建的近景影像定向方法，可以在無任何物方控制的條件下利用3張以上的未檢校影像的同名點恢復(fù)相機(jī)的內(nèi)外方位元素，是一種有別于傳統(tǒng)攝影測量解析處理方法的新方法，但由于無物方控制和未考慮畸變，因此定向的精度相對較低，只能用于低精度的測量或作為光束法平差的初值.該方法的核心是射影重建和自檢校，而自檢校的精度取決于射影重建的精度，射影重建的精度除了相機(jī)本身的因素外，還與同名點的數(shù)量、分布、重疊度、量測精度等有關(guān)，因此在作業(yè)過程中要避免粗差，盡可能提高量測精度并保持較好的成像構(gòu)形，否則可能導(dǎo)致計算失敗.

［1］Heyden A，?tr?sm K.Euclidean reconstruction from constant intrinsic parameters［C］／／Proceeding International Conference of Pattern Recognition.Vienna：［s.n］，1996：339－343.

［2］Hartley R L.Self－calibration of stationary cameras［J］.Inter Journal of Computer Vision，1997，22（1）：5.

［3］Pan H P.A direct closed－form solution to general relative orientation of two stereo views［J］.Digital Signal Processing，1999，9（27）：195.

［4］Faugeras O.Stratification of 3D vision：projective，affine，and metric representations［J］.Journal of the Optical Society of America，1995，12（3）：465.

［5］Pollefeys M.Self－calibration and metric 3D reconstruction from uncalibrated image sequences ［D］.Leuven： Katholieke Universiteit.Department of Electrical Engineering，1999.

［6］邱志強(qiáng)，于起峰.從建筑物序列圖像恢復(fù)三維結(jié)構(gòu)［J］.武漢大學(xué)學(xué)報：工學(xué)版，2003，36（3）：13.QIU Zhiqiang，YU Qifeng.3D building structure from image sequences［J］.Engineering Journal of Wuhan University：Engineering Edition，2003，36（3）：13.

［7］Wang G H，Wu J.Quasi－perspective projection model：theory and application to structure and motion factorization from uncalibrated image sequences［J］.Inter Journal of Computer Vision，2010，87（3）：213.

［8］Wang G H，Wu J.Perspective 3－D euclidean reconstruction with varying camera parameters［J］.IEEE Trans on Circuits and Systems for Video Technology，2009，19（12）：1793.

［9］Tomasi C，Kanade T.Shape and motion from image stream under orthography：a factorization method［J］.International Journal of Computer Vision，1992，9（2）：137.

［10］Sturm P，Triggs B.A factorization based algorithm for multiimage projective structure and motion［C］／／Proceeding 4th European Conference on Computer Vision.［S.l.］：Springer，1996：709－720.

［11］Heyden A，Berthilsson R，Sparr G.An iterative factorization method for projective structure and motion from image sequences［J］.Image and Vision Computing，1999，17（13）：981.

［12］Zhang Z.A flexible new technique for camera calibration［J］.IEEE Trans Pattern Analysis and Machine Intelligence，2002，22（11）：1330.