陳 輝,張 卡,4,宿 東,王蓬勃
(1. 南京師范大學虛擬地理環(huán)境教育部重點實驗室,江蘇 南京 210023; 2. 江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心,江蘇 南京 210023; 3. 南京師范大學地理科學學院,江蘇 南京 210023; 4. 江蘇省地理環(huán)境演化國家重點實驗室培育建設點,江蘇 南京 210023)
城市建筑物高度對建筑物容積率、城市三維重建、城市環(huán)境等都具有顯著的影響[1]。近年來,隨著高分辨率遙感技術的日益發(fā)展,利用高分辨率遙感影像獲取建筑物的高度吸引了國內(nèi)外許多學者的研究。目前的研究方法總體上可分為兩類:一是利用高分辨率影像立體像對,根據(jù)空中三角測量原理來計算建筑物高度[2-4],該方法理論精度較高,但需要具有一定重疊度的兩幅或多幅影像;二是根據(jù)太陽、衛(wèi)星方位角度、陰影長度與建筑物的空間幾何關系來估計建筑物高度[5-8],該方法應用較多,計算量少且成本較低,但是受陰影質(zhì)量的影響較大,當遇到建筑物陰影相互遮擋等復雜情況時,建筑物高度估計的精度很難達到實際應用要求。
針對當前研究方法的不足,本文根據(jù)衛(wèi)星側(cè)視成像情況下因高差引起的建筑物像點位移與建筑物高度之間的幾何關系,提出基于建筑物側(cè)面輪廓線的高分辨遙感影像建筑物高度估算方法,利用實際高分辨率遙感影像進行試驗驗證,并與基于陰影建筑物高度提取法進行對比,試驗結(jié)果表明本文方法可以得到較高精度的建筑物高度估算結(jié)果。
本文根據(jù)衛(wèi)星側(cè)視成像時因高差引起的建筑物像點位移與建筑物高度的關系,通過提取影像中的建筑物側(cè)面輪廓線并計算其長度,從而反算建筑物的高度。本文方法的總體技術路線如圖1所示。
本文利用有理函數(shù)模型計算像點位移的角度技術流程為:先采用RPC模型的反解形式(見式(1)),取圖像上一點,坐標為(r0,c0),初始高程為該地區(qū)的平均高程Z0,經(jīng)過計算便得到一個虛擬的物點(X0,Y0,Z0),然后以10 m為間隔取Z0附近7個點;再利用RPC模型的正解形式(見式(2))求得對應的7個像點,對這些像點進行函數(shù)擬合,根據(jù)相關核線理論,航天遙感影像的核線在小范圍內(nèi)可看做直線[9-11],且WorldView系列影像核線具有很好的平行特性[12],因此這幾個像點可以用一條直線進行擬合,像點位移的方向便是該直線的方向α(α=tan-1k,k為擬合直線的斜率)。
(1)
(2)
式中,(r,c)和(X,Y,Z)分別為經(jīng)過歸一化處理的像點的像方坐標及其對應的物方點坐標;P1,P2,…,P8為有理多項式函數(shù),函數(shù)系數(shù)可從衛(wèi)星影像的RPC文件中獲得。
相關研究表明[13],線陣遙感衛(wèi)星側(cè)視成像時像移誤差公式為
(3)
式中,Dh為像點位移的大小;h為地面高度;H為衛(wèi)星軌道高度;f為傳感器焦距;Rn為中心像元到第n個成像像元的距離,若Rn以向外側(cè)CCD方向為正方向,取正,若以向內(nèi)側(cè)CCD方向為負方向,取負。一般而言,Rn相比ftanθ是一個極小量,因此在極小范圍內(nèi)Rn可看做常數(shù),即ΔDh與h正相關,由于影像中通過使用RPC模型計算得到的7個虛擬像點分別對應著不同的高程h,不同像點間的距離ΔDh同樣與其對應的物點高程之差Δh成正比,因此可求得衛(wèi)星側(cè)視成像的角度θ如下
(4)
利用前面求得的像點位移角度α將圖像進行旋轉(zhuǎn),使建筑物像點位移的方向沿水平方向。為了獲得圖像建筑物輪廓線,首先,采用Canny[14]邊緣檢測算法提取建筑物輪廓的二值圖像,然后對二值圖像進行形態(tài)學閉運算填充小型空洞;然后,構(gòu)建矩形形態(tài)學結(jié)構(gòu)元素SE(structure element)對二值圖像進行開運算提取圖像中的橫線,矩形形態(tài)學結(jié)構(gòu)元素的尺寸和方向直接決定檢出線狀地物的粗細和方向,由于建筑物輪廓線都是水平的,因此本文試驗主要采用寬度為1、長度一定的結(jié)構(gòu)元素,長度的選取與建筑物像點位移的距離有關,像點偏移越嚴重,長度越長;最后,利用Hough變換,設置角度閾值為90°檢測圖像中的橫線,并將橫線編碼,從而得到建筑物的輪廓線。
其中,t1、t2分別為直線坐標x、y距離閾值。橫線經(jīng)過合并之后,尋找距離端點處最近的角點,如果提取的橫線不能正確契合建筑物的側(cè)面輪廓線,無論是“超出”或是“過短”的情況,都可以利用與之一端最近的角點進行約束,從而對輪廓線的長度進行修正。同樣,如果同一建筑物有多條橫線,則取修正后最長的橫線為建筑物側(cè)面輪廓線的長度。
(5)
式中,f為傳感器焦距;H為衛(wèi)星軌道高度;Rn為中心像元到第n個成像像元的距離;θ為衛(wèi)星側(cè)視成像的角度。
為了驗證本文算法具有比陰影法更好的適用性及準確性,在Windows 10環(huán)境下,采用Visual C# .NET 2013程序設計語言結(jié)合Emgu CV函數(shù)庫對本文算法進行編程實現(xiàn),并采用兩組不同地區(qū)的遙感影像數(shù)據(jù)進行對比試驗。兩組數(shù)據(jù)都是WorldView-3衛(wèi)星4波段多光譜影像,分辨率為1.24 m,單波段全色影像分辨率為0.31 m及0.5 m的DSM數(shù)據(jù)(用作建筑物真值作為比對)。多光譜影像與全色影像融合之后試驗影像如圖3所示。其中,影像(a)位于巴西里約熱內(nèi)盧,尺寸為1302×1086像素,該影像中建筑物的陰影相互遮擋,無法有效提取出建筑物陰影,因此無法使用陰影法進行高度估算;影像(b)位于利比亞地區(qū),尺寸為1629×1606像素,該影像可以很好地提取陰影的長度。
在求出旋轉(zhuǎn)角度α將影像進行旋轉(zhuǎn)之后,為了更好地進行建筑物輪廓線的提取,本文將原圖像經(jīng)過雙邊濾波之后設置Canny算子高低閾值為200進行輪廓提取,隨后進行形態(tài)學閉運算,填充輪廓小空洞,使輪廓具有更好的連貫性。在進行建筑物的橫線提取時,形態(tài)學結(jié)構(gòu)元素的構(gòu)造對最終的試驗結(jié)果至關重要[16],因此試驗采用長度為15個像素的矩形結(jié)構(gòu)元素對經(jīng)過處理后的圖像進行開運算,使用合并精簡規(guī)則處理之后,最終的結(jié)果如圖4所示。最后設置FAST角點檢測器的閾值為60提取圖像的角點信息,以此修正橫線的距離,從而最終確定建筑物側(cè)面輪廓線的長度。
本次試驗每組數(shù)據(jù)各隨機選取了15座建筑物,其真實高度數(shù)值從DSM數(shù)據(jù)上獲取。在求得衛(wèi)星成像角度θ之后,利用建筑物側(cè)面輪廓線的長度反算建筑物高度,并與其真值進行絕對誤差和相對誤差的計算。里約地區(qū)試驗影像的建筑物高度提取結(jié)果見表1;對于利比亞地區(qū)的試驗影像,利用本文方法和文獻[17]中陰影法的建筑物高度提取對比試驗結(jié)果如圖5所示。
表1 里約地區(qū)建筑物高度測量結(jié)果與分析
根據(jù)試驗結(jié)果,里約地區(qū)影像建筑物計算結(jié)果絕對誤差的范圍為0.03~3.43 m,相對誤差范圍為0.12%~7.77%,平均誤差為0.74 m。而在利比亞影像的試驗結(jié)果中,本文方法的建筑物高度估算絕對誤差為0.29~2.06 m,相對誤差為0.52%~5.37%;陰影法的絕對誤差為2.06~5.26 m,相對誤差為2.63%~10.4%。因此,本文方法的建筑物高度估算精度明顯優(yōu)于陰影法。
另外,本文方法的誤差大小與建筑物底部被遮掩的情況有關。如果建筑物靠近地面的部分被植物遮擋或處于混合像元,則本文方法進行建筑物高度測算誤差相對較大。但是,在一定誤差允許的范圍內(nèi),利用本文方法計算的建筑物高度的精度仍然較高,有些甚至可以達到分米級,達到甚至超過了利用陰影法計算建筑物高度的精度。
本文提出了一種基于建筑物側(cè)面輪廓線的高分辨率影像建筑物高度計算的新方法,相對于已有的利用陰影求算建筑物高度的方法,本文方法只需要解算像點位移的方向,以及衛(wèi)星側(cè)視成像的角度(在此角度已知的情況下計算量會更少);而且,本文方法還可以應用于陰影混淆相互遮擋、衛(wèi)星夜間成像等陰影法無法應用的情況。通過使用實際遙感影像的試驗結(jié)果分析,本文提出的方法比基于陰影的建筑物高度估算方法的精度更高,且具有更廣泛的適應性。但是,面對建筑物過于密集、建筑物側(cè)面輪廓相互遮擋等復雜情況,本文方法還存在一定的不足。在未來的研究工作中,還需要考慮以上更加復雜的情況,以提高本文方法的適用性。