蒲泓亦,李 鋒,徐 錚

(信息工程大學(xué),河南 鄭州 450001)

空中全景是利用旋翼無(wú)人機(jī)在空中懸停拍攝實(shí)景影像,對(duì)影像進(jìn)行融合處理制成的全方位全視角實(shí)景地圖,可將觀察視角從地面轉(zhuǎn)移到空中,實(shí)現(xiàn)從地平線到地平線的全景觀察,獲得空中俯視、遨游天空的體驗(yàn)感,對(duì)空間環(huán)境整體信息的掌握意義重大。國(guó)內(nèi)區(qū)域可根據(jù)任務(wù)需要實(shí)地拍攝獲取空中全景數(shù)據(jù),但受國(guó)際關(guān)系和法律法規(guī)制約,境外陌生地域尤其是敵對(duì)國(guó)家或地區(qū)無(wú)法實(shí)地拍攝獲取全景數(shù)據(jù)。隨著無(wú)人機(jī)技術(shù)的普及,國(guó)內(nèi)外無(wú)人機(jī)航拍愛(ài)好者利用消費(fèi)級(jí)無(wú)人機(jī)在高空拍攝全景影像,并上傳到網(wǎng)絡(luò)平臺(tái)供大眾瀏覽。這些眾源空中全景影像包含城市、海岸、港口、政府機(jī)構(gòu)等目標(biāo)信息,對(duì)境外戰(zhàn)場(chǎng)環(huán)境勘察有重要意義。但由于缺少準(zhǔn)確的地理位置和方向信息,眾源空中全景影像無(wú)法在地理空間中精確定位定向和融合匹配,限制了其應(yīng)用場(chǎng)景。

針對(duì)全景影像獲取與應(yīng)用問(wèn)題,國(guó)內(nèi)外已展開了大量研究。根據(jù)全景相機(jī)的成像原理,全景模型大致分為嚴(yán)格球面全景模型和理想球面全景模型[1-2]。嚴(yán)格球面全景模型強(qiáng)調(diào)全景成像過(guò)程中多鏡頭投影中心存在偏移和旋轉(zhuǎn),采用旋轉(zhuǎn)矩陣和平移向量構(gòu)建統(tǒng)一的球面全景坐標(biāo)系?；诖四Ｐ偷亩噻R頭全景相機(jī)被應(yīng)用到車載移動(dòng)測(cè)圖系統(tǒng)中,系統(tǒng)利用GPS+IMU提供的輔助數(shù)據(jù)平差計(jì)算球面全景影像位姿參數(shù)[3-4]。除了輔助地面載體進(jìn)行測(cè)量外,全景影像逐步應(yīng)用到無(wú)人機(jī)攝影測(cè)量工程實(shí)踐中。無(wú)人機(jī)全景傾斜攝影測(cè)量將空中全景影像投影至立方體表面,對(duì)立方體表面的影像實(shí)施空中三角測(cè)量,重建場(chǎng)景三維模型[5]。全景攝影測(cè)量則是以理想球面全景模型為基礎(chǔ)進(jìn)行測(cè)量,屬于近景攝影測(cè)量,主要用于建筑物三維重建[6-7]。算法方面,球面全景影像的快速位姿估計(jì)算法EPnP(efficient perspective-n-point)利用虛擬控制點(diǎn)表示像點(diǎn)和物方點(diǎn),求得虛擬控制點(diǎn)的像空間坐標(biāo)和物方空間坐標(biāo),根據(jù)歐式變換不變性解算影像位姿參數(shù)[8]。點(diǎn)-線特征聯(lián)合算法融合影像上點(diǎn)、線特征的解算數(shù)學(xué)模型進(jìn)行平差,求解位姿參數(shù)[9]?；陔p球面投影的相對(duì)定向-絕對(duì)定向影像位姿估計(jì)方法,利用雙目視覺(jué)原理計(jì)算視差,構(gòu)建可量測(cè)球形立體全景模型[10-11]。文獻(xiàn)[12]提出了適合球面全景成像特點(diǎn)的相對(duì)定向計(jì)算流程,結(jié)果符合攝影測(cè)量標(biāo)準(zhǔn)。當(dāng)前,國(guó)內(nèi)外全景影像位姿估計(jì)的大部分研究方法需要借助GPS+IMU輔助數(shù)據(jù),增加方法成本和計(jì)算復(fù)雜性,在眾源空中全景影像位姿解算方面尚無(wú)成熟的方法。

針對(duì)眾源空中全景影像缺少位姿參數(shù)的問(wèn)題,本文提出基于透視變換的眾源空中全景位姿解算方法。根據(jù)空間后方交會(huì)原理和球面模型特征,設(shè)計(jì)基于透視變換的空中全景位姿估計(jì)算法,提出解算流程;重點(diǎn)研究“平-球-平”變換中模型和坐標(biāo)變換,采用虛擬平面影像外方位元素描述全景姿態(tài);依據(jù)透視變換的變形特征,設(shè)計(jì)顧及誤差的迭代重加權(quán)最小二乘法,并通過(guò)案例驗(yàn)證方法的可行性和可靠性。

1 算法原理與步驟

1.1 算法原理

空間后方交會(huì)法根據(jù)成像特征建立共線條件方程,利用泰勒展開線性化方程,采用最小二乘法迭代求解6個(gè)外方位元素(XS,YS,ZS,φ,ω,κ)。其中,前3個(gè)是線元素,確定像片及其投影中心在物方空間坐標(biāo)系中的三維坐標(biāo);后3個(gè)是角元素,分別為偏角、傾角、旋角,確定像空間坐標(biāo)系三軸在物方空間坐標(biāo)系中的方向。普通影像的數(shù)學(xué)模型為垂直于像空間坐標(biāo)系z(mì)軸的平面,像點(diǎn)z軸坐標(biāo)為定值:z=-f。球面全景像點(diǎn)的z坐標(biāo)隨像點(diǎn)位置的變化而變化,直接代入公式會(huì)增加變量,加大解算的復(fù)雜程度。本文提出的方法使用透視變換將球面全景模型轉(zhuǎn)換為虛擬平面影像模型,通過(guò)解算虛擬平面影像外方位元素求得全景影像位姿參數(shù),簡(jiǎn)化模型復(fù)雜度。方法核心為兩個(gè)變換。

(1) “平-球-平”模型變換。平面全景模型是眾源空中全景影像的原始模型;球面全景模型由平面全景模型球面建模得到,其位姿參數(shù)是方法解算的根本對(duì)象;虛擬平面影像模型是球面模型經(jīng)球心透視變換產(chǎn)生的平面,作為連接球面像點(diǎn)和地面控制點(diǎn)的橋梁,其位姿參數(shù)是算法解算的直接對(duì)象。

(2)坐標(biāo)系變換。利用旋轉(zhuǎn)矩陣和平移向量將像點(diǎn)和控制點(diǎn)的坐標(biāo)系轉(zhuǎn)換為虛擬平面坐標(biāo)系。具體轉(zhuǎn)換過(guò)程如下:像素坐標(biāo)系向虛擬平面像空間坐標(biāo)系轉(zhuǎn)換;地理坐標(biāo)系向重心東北天坐標(biāo)系轉(zhuǎn)換;重心東北天坐標(biāo)系向虛擬平面像空間坐標(biāo)系變換。

考慮到透視變換的變形特征、控制點(diǎn)選擇誤差及測(cè)繪數(shù)據(jù)測(cè)量誤差,方法在約束平面全景像點(diǎn)選擇范圍的基礎(chǔ)上,設(shè)計(jì)顧及誤差的迭代重加權(quán)最小二乘法解算全景位姿,以提高算法準(zhǔn)確性和穩(wěn)定性。

1.2 解算流程

與實(shí)拍空中全景影像不同,眾源空中全景影像一般僅含有攝站所在縣市和景區(qū)的名稱屬性或攝站大致的地理位置信息。由于解算位姿的前提是已知一定數(shù)量的地面點(diǎn)物方空間坐標(biāo)和對(duì)應(yīng)像點(diǎn)像空間坐標(biāo),設(shè)計(jì)歷史測(cè)繪數(shù)據(jù)輔助的眾源空中全景位姿解算流程,如圖1所示。

圖1 解算流程

步驟1:基于遙感影像、DEM和傾斜攝影模型構(gòu)建三維地理環(huán)境,作為全景影像粗定位的地理空間。

步驟2:獲取眾源空中全景影像,篩選成像質(zhì)量較好的空中全景影像,預(yù)處理長(zhǎng)寬分辨率比達(dá)到2∶1。

步驟3:進(jìn)行球面全景建模。使用透視變換,將球面影像變換為虛擬平面影像,建立虛擬平面坐標(biāo)系。

步驟4:在三維地理環(huán)境中選擇地面控制點(diǎn),確定對(duì)應(yīng)像點(diǎn)。

步驟5:進(jìn)行像點(diǎn)和地面控制點(diǎn)坐標(biāo)系變換,將兩者統(tǒng)一到虛擬平面坐標(biāo)系。

步驟6:使用迭代重加權(quán)最小二乘法解算眾源空中全景影像位姿。

2 空中全景變換

2.1 球面全景建模

球面全景是將平面全景影像映射到球面上得到的球面影像[13],原理如下:以影像中心O′為原點(diǎn),U′軸平行于像素坐標(biāo)系U軸方向相同,V′軸平行于V軸方向相反,構(gòu)建輔助像素坐標(biāo)系。像素點(diǎn)P′像素坐標(biāo)(u,v)與輔助像素坐標(biāo)系下坐標(biāo)(u′,v′)的轉(zhuǎn)換關(guān)系為

(1)

式中,l和w分別為平面全景影像的長(zhǎng)和寬,比值為2∶1。

如圖2所示,采用等距圓柱投影將平面全景影像映射至球面時(shí),平面全景影像長(zhǎng)和寬分別覆蓋球面的經(jīng)度和緯度,即[0,l]對(duì)應(yīng)球面經(jīng)度[-π, π],[0,w]對(duì)應(yīng)球面緯度[-π/2, π/2]。以球心OS為原點(diǎn),X軸指向O′,Y軸垂直于X軸構(gòu)成球赤道平面XOSY,Z軸與X、Y軸構(gòu)成右手坐標(biāo)系,建立球空間坐標(biāo)系OS-XYZ。從Z軸正方向看,順時(shí)針球面經(jīng)度為負(fù),逆時(shí)針為正。球面投影點(diǎn)PS球面經(jīng)緯度(α,β)與點(diǎn)P′輔助像素坐標(biāo)(u′,v′)變換關(guān)系為

圖2 球面全景建模

(2)

2.2 全景模型變換

空間后方交會(huì)法根據(jù)成像特點(diǎn)構(gòu)造的共線條件方程為

(3)

式中,(x,y,z)為像點(diǎn)的像空間坐標(biāo);f為相機(jī)焦距;(X,Y,Z)為控制點(diǎn)物空間坐標(biāo);ai、bi、ci(i=1, 2, 3)為角元素的三角函數(shù)代數(shù)式。在球面全景影像中,像點(diǎn)PS的z值隨像點(diǎn)球面緯度變化,并不符合經(jīng)典空間后方交會(huì)的解算條件,通常的改進(jìn)方式是同時(shí)列出x、y、z的方程,大大增加計(jì)算成本。本文對(duì)球面全景模型進(jìn)行“球-平”模型變換,即利用透視變換將球面全景影像映射到虛擬平面,得到虛擬像平面像點(diǎn),采用平面空間后方交會(huì)法進(jìn)行計(jì)算,降低計(jì)算量。

如圖3所示,設(shè)球空間直角坐標(biāo)系Z軸負(fù)方向與球面交于點(diǎn)O1,過(guò)點(diǎn)O1并垂直于Z軸構(gòu)建切平面xO1y。球面像點(diǎn)PS(α,β)和虛擬平面映射點(diǎn)p(x,y,z)連線方向交地面于點(diǎn)P。點(diǎn)PS與點(diǎn)p映射關(guān)系為

圖3 全景模型變換

(4)

根據(jù)上述模型構(gòu)建原則可知:虛擬平面影像的外方位元素與球面全景影像的位置和姿態(tài)相同,即可以采用虛擬平面影像的外方位元素描述球面全景影像位姿。

3 空中全景位姿解算

3.1 建立物方空間坐標(biāo)系

(5)

式中,L和B分別為重心P0的經(jīng)度和緯度。

3.2 加權(quán)重迭代最小二乘法參數(shù)解算

圖4 “球-平”模型變換示意圖

(6)

(7)

(8)

根據(jù)式(8)可知,變換前后長(zhǎng)度產(chǎn)生變形。對(duì)Δd求導(dǎo)可知:當(dāng)β∈[-π/2,0]時(shí),β減小,Δd減小且變化速率降低;β增大,Δd增大且變化速率增大,即

(9)

(10)

(11)

式中,Δd′為Δd的導(dǎo)數(shù)。極點(diǎn)O1處長(zhǎng)度變形為0,赤道處單位弧長(zhǎng)與虛擬平面投影長(zhǎng)度的差值為無(wú)窮。式(9)說(shuō)明:對(duì)球面全景模型下半球進(jìn)行透視變換,實(shí)際得到一個(gè)垂直于Z軸且過(guò)極點(diǎn)O1的超平面。將遙感影像、DEM和傾斜攝影模型作為控制點(diǎn)選取空間,控制點(diǎn)的坐標(biāo)數(shù)據(jù)不可避免地包含了測(cè)量誤差,其數(shù)值大小通常不一致。同時(shí),人工選擇像點(diǎn)的過(guò)程會(huì)引入偶然誤差,尤其是選擇赤道附近像點(diǎn)時(shí),單位像素的偏差將引起虛擬平面像點(diǎn)極大的偏移,對(duì)位姿參數(shù)求解產(chǎn)生較大誤差。因此,“平-球-平”模型變換對(duì)控制點(diǎn)對(duì)(像點(diǎn)和對(duì)應(yīng)地面控制點(diǎn))的質(zhì)量提出更高要求,控制點(diǎn)對(duì)數(shù)目和分布都會(huì)對(duì)解算結(jié)果產(chǎn)生較大的影響。為盡量減少偶然誤差,方法設(shè)定選取像點(diǎn)的緯度區(qū)間為β∈[-π/2,δ](-π/2 <δ<0)。

盡管通過(guò)限制像點(diǎn)選擇范圍降低了偶然誤差,但仍然存在結(jié)果不準(zhǔn)確的現(xiàn)象。經(jīng)典的空間后方交會(huì)認(rèn)為每個(gè)控制點(diǎn)對(duì)解算結(jié)果的貢獻(xiàn)度相同,然而對(duì)于從歷史測(cè)繪數(shù)據(jù)中選取的控制點(diǎn)坐標(biāo),其測(cè)量誤差的分布與大小未知,如果不區(qū)別直接進(jìn)行相機(jī)位姿估計(jì),可能導(dǎo)致估計(jì)結(jié)果與真值相差甚遠(yuǎn)。因此采用顧及誤差的迭代重加權(quán)最小二乘法解算眾源空中全景影像位姿參數(shù),以達(dá)到提高高質(zhì)量控制點(diǎn)貢獻(xiàn),降低低質(zhì)量控制點(diǎn)貢獻(xiàn)的目的。在已知外方位元素近似值的情況下,按照式(3)一對(duì)控制點(diǎn)能夠組成2個(gè)方程式。理論上,利用3個(gè)不在同一直線上的平高控制點(diǎn)即可解算位姿參數(shù)近似值的改正數(shù)。實(shí)際計(jì)算中一般選取4對(duì)或4對(duì)以上控制點(diǎn),計(jì)算過(guò)程如下。

誤差方程式矩陣形式[15]為

AX-L=V

(12)

其中

(13)

(14)

(15)

式中,A中元素由3個(gè)角元素決定;X為外方位元素改正數(shù)向量;L中l(wèi)x=x-x計(jì),ly=y-y計(jì);V為誤差向量。

迭代重加權(quán)最小二乘法的目標(biāo)函數(shù)為

(16)

解算外方位元素改正數(shù)為

X=[ATWTWA]-1ATWTWL

(17)

式中,W為權(quán)值矩陣,每次迭代需不斷更新矩陣元素wi。

式(12)是線性化近似公式,因此必須采用逐次趨近的迭代計(jì)算。將前一次改正的位姿參數(shù)作為近似值,重復(fù)建立誤差方程求解,利用Xi(i=1,2,…,n)不斷糾正位姿參數(shù),直至改正數(shù)的絕對(duì)值小于閾值,或像點(diǎn)坐標(biāo)的測(cè)量值與計(jì)算值之間的差值小于規(guī)定限差。

3.3 權(quán)值選取

權(quán)值分配在控制點(diǎn)的測(cè)量精度這一問(wèn)題上應(yīng)該遵循控制點(diǎn)的可信度原則,通過(guò)權(quán)值的分配調(diào)整異常值對(duì)整體誤差的影響[16]。采用改進(jìn)的Huber函數(shù)[17]確定權(quán)重矩陣元素wi,既能夠快速收斂,又能夠平衡異常值對(duì)整體計(jì)算的影響。以每一次迭代中虛擬平面上的像方殘差為量度,根據(jù)殘差大的點(diǎn)所對(duì)應(yīng)的測(cè)量誤差和深度影響劇烈這一特點(diǎn)可得權(quán)重系數(shù)的表達(dá)式[18]為

(18)

4 應(yīng)用案例及分析

以河南省登封市紙坊水庫(kù)地區(qū)的空中全景影像為案例對(duì)本文方法的有效性進(jìn)行驗(yàn)證。試驗(yàn)采用全國(guó)底層遙感影像、DEM數(shù)據(jù)及紙坊水庫(kù)傾斜攝影模型建立虛擬地理環(huán)境。全局坐標(biāo)系采用WGS-84坐標(biāo)系,參考橢球長(zhǎng)半徑a=6378 137 m,短半徑b=6 356 752.314 2 m,扁率e=1/298.257 223 6。

從720云網(wǎng)址下載河南省登封市紙坊水庫(kù)地區(qū)的多張空中全景影像,分辨率為8192×4096像素。由于地面控制點(diǎn)極少情況出現(xiàn)在影像上半球,且需要限制控制點(diǎn)選擇范圍,設(shè)定緯度區(qū)間β∈[-π/2,-π/4],即像點(diǎn)坐標(biāo)分布在[0,8192]×[3072,4096]范圍。設(shè)攝站初始高差Z0=50 m,等效焦距f=28 mm,閾值為0.001°。圖5(a)以水庫(kù)全景影像為例,說(shuō)明了選取地面控制點(diǎn)并完成計(jì)算的過(guò)程,所選控制點(diǎn)見表1。

表1 控制點(diǎn)對(duì)坐標(biāo)

圖5 計(jì)算步驟和效果

計(jì)算得到全景影像外方位元素為

依據(jù)外方位元素,得到全景影像的地理坐標(biāo)為(*** .525 093°N,*** .105 866°E,*** .36 m)。如圖5(b)所示,根據(jù)所得到地理坐標(biāo)和角元素,將全景影像融合至虛擬地理環(huán)境內(nèi),發(fā)現(xiàn)影像定位在準(zhǔn)確位置,且進(jìn)入全景前后虛擬相機(jī)視線方向保持一致,說(shuō)明了本文算法的可行性和準(zhǔn)確性。

在驗(yàn)證算法準(zhǔn)確性的基礎(chǔ)上,本文對(duì)算法的復(fù)雜度進(jìn)行測(cè)試,使用迭代次數(shù)作為指標(biāo)。試驗(yàn)中,控制點(diǎn)數(shù)目從6增加至16,步長(zhǎng)為1,比較本文算法、球面后方交會(huì)算法[19]和球面迭代重加權(quán)后方交會(huì)算法收斂至閾值內(nèi)的迭代次數(shù)(如圖6所示)。其中,球面迭代重加權(quán)后方交會(huì)算法根據(jù)成像特點(diǎn)構(gòu)建x、y、z的共線條件方程,采用迭代重加權(quán)最小二乘法解算。

圖6 算法迭代次數(shù)對(duì)比

由圖6可知,與球面后方交會(huì)算法和球面迭代重加權(quán)后方交會(huì)算法相比,本文方法在有冗余控制點(diǎn)時(shí),明顯減少了算法的計(jì)算量。

5 結(jié) 語(yǔ)

解算影像位姿參數(shù)是眾源空中全景影像應(yīng)用于境外戰(zhàn)場(chǎng)環(huán)境勘察的基礎(chǔ)。本文提出了一種基于透視變換的眾源空中全景位姿解算方法,設(shè)計(jì)完整解算流程;采用透視變換簡(jiǎn)化數(shù)學(xué)模型,著重闡明模型變換和坐標(biāo)變換,減少了位姿解算計(jì)算量;使用顧及誤差的迭代重加權(quán)最小二乘法解算位姿參數(shù),增強(qiáng)了方法的穩(wěn)定性。最后,以登封市紙坊水庫(kù)眾源空中全景影像為例進(jìn)行驗(yàn)證。試驗(yàn)表明,本文方法能夠?qū)崿F(xiàn)眾源空中全景影像位姿解算;存在冗余控制點(diǎn)情況下,既能減少計(jì)算量又能保持較高的穩(wěn)定性。研究成果為眾源空中全景影像與地理信息空間融合等應(yīng)用提出一種新的實(shí)現(xiàn)方法,為境外戰(zhàn)場(chǎng)環(huán)境勘察提供重要參考。