蒲泓亦,李 鋒,徐 錚

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

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

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

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

1 算法原理與步驟

1.1 算法原理

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

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

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

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

1.2 解算流程

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

圖1 解算流程

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

步驟2:獲取眾源空中全景影像,篩選成像質量較好的空中全景影像,預處理長寬分辨率比達到2∶1。

步驟3:進行球面全景建模。使用透視變換,將球面影像變換為虛擬平面影像,建立虛擬平面坐標系。

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

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

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

2 空中全景變換

2.1 球面全景建模

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

(1)

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

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

圖2 球面全景建模

(2)

2.2 全景模型變換

空間后方交會法根據成像特點構造的共線條件方程為

(3)

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

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

圖3 全景模型變換

(4)

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

3 空中全景位姿解算

3.1 建立物方空間坐標系

(5)

式中,L和B分別為重心P0的經度和緯度。

3.2 加權重迭代最小二乘法參數解算

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

(6)

(7)

(8)

根據式(8)可知,變換前后長度產生變形。對Δd求導可知:當β∈[-π/2,0]時,β減小,Δd減小且變化速率降低;β增大,Δd增大且變化速率增大,即

(9)

(10)

(11)

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

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

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

AX-L=V

(12)

其中

(13)

(14)

(15)

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

迭代重加權最小二乘法的目標函數為

(16)

解算外方位元素改正數為

X=[ATWTWA]-1ATWTWL

(17)

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

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

3.3 權值選取

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

(18)

4 應用案例及分析

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

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

表1 控制點對坐標

圖5 計算步驟和效果

計算得到全景影像外方位元素為

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

在驗證算法準確性的基礎上,本文對算法的復雜度進行測試,使用迭代次數作為指標。試驗中,控制點數目從6增加至16,步長為1,比較本文算法、球面后方交會算法[19]和球面迭代重加權后方交會算法收斂至閾值內的迭代次數(如圖6所示)。其中,球面迭代重加權后方交會算法根據成像特點構建x、y、z的共線條件方程,采用迭代重加權最小二乘法解算。

圖6 算法迭代次數對比

由圖6可知,與球面后方交會算法和球面迭代重加權后方交會算法相比,本文方法在有冗余控制點時,明顯減少了算法的計算量。

5 結 語

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