吳浩霖 聶高眾 范熙偉 魏本勇 安基文
摘要:選取新疆阿圖什市下轄的瓊哈拉峻村為研究區(qū),將小型旋翼無人機拍攝的圖像作為數(shù)據(jù)源,分別采用面向?qū)ο笠约懊嫦蛳裨?種影像分析方法對研究區(qū)的房屋建筑進行提取,并對2種算法的提取結(jié)果進行對比,分析了各自的優(yōu)勢。結(jié)果表明:面向?qū)ο蠓椒梢杂行У厝コ符}噪聲對分類帶來的影響,保證房屋形態(tài)的完整性,但影響內(nèi)部相似的光譜、紋理信息若對應(yīng)多種物體則會導(dǎo)致影像對象的錯分。在面向像元的提取方法中加入了改進的數(shù)學(xué)形態(tài)學(xué)算法,可有效的抑制椒鹽噪聲,保持建筑物邊緣的連續(xù)性與完整性,較好地解決了面向?qū)ο蠓椒ㄖ胁糠洲r(nóng)田與房屋出現(xiàn)錯分的問題。
關(guān)鍵詞:地震災(zāi)害損失預(yù)評估;無人機;地表建筑提取;形態(tài)學(xué)圖像處理
中圖分類號:P315.943 ??文獻標識碼:A ??文章編號:1000-0666(2019)02-0236-09
0 引言
地震災(zāi)害損失預(yù)評估是震后及時、科學(xué)地進行地震應(yīng)急指揮的關(guān)鍵之一(安基文等,2015;聶高眾等,2002)。地震導(dǎo)致人員傷亡的因素中有80%是房屋破壞造成的。在進行地震災(zāi)害損失預(yù)評估的過程中,關(guān)鍵要做好評估區(qū)域房屋建筑狀況,以及人口規(guī)模的調(diào)查(聶高眾等,2012)。而對于研究區(qū)域的房屋建筑提取工作則是上述調(diào)查的有力保證。
由于無人機遙感技術(shù)相較于星載遙感技術(shù)擁有諸多優(yōu)勢(Berra et al,2016;Dandois,Ellis,2013;Messinger et al,2016),在地震災(zāi)害損失預(yù)評估的實際工作中無人機逐漸開始代替?zhèn)鹘y(tǒng)的星載遙感技術(shù)對地表建筑物信息進行獲?。ǜ恫?,2018;臧克等,2010;李祖?zhèn)鞯龋?010;彭大雷等,2017;周洋等,2017;陳晉等,2018)。在調(diào)查過程中每一個區(qū)域的調(diào)查點繁多,且每一個調(diào)查點的調(diào)查時間有限,所攜帶的大多數(shù)為飛行準備時間短、操作靈活、易于攜帶的小型旋翼無人機,而大多數(shù)小型無人機所搭載的傳感器為數(shù)字相機,拍攝得到的影像雖然圖像分辨率高,但是光譜分辨率很低,即所返回的影像數(shù)據(jù)波段數(shù)較少。本文利用面向?qū)ο笠约懊嫦蛳裨?種影像分析方法對研究區(qū)的房屋建筑進行提取,并對2種算法的提取結(jié)果進行對比,分析了在該數(shù)據(jù)條件下2種方法各自的優(yōu)勢。
1 數(shù)據(jù)來源及影像拼接
本文選取的研究區(qū)為新疆阿圖什市下轄的哈拉峻鄉(xiāng)瓊哈拉峻村(40.178°N,76.832°E)。哈拉峻鄉(xiāng)90 m空間分辨率數(shù)字高程模型(Digital Elevation Model,DEM)疊加相應(yīng)的山地陰影Hillshade圖,如圖1a所示;無人機飛行路線規(guī)劃如圖1b所示。
本文使用小型四旋翼無人機,為大江精靈3專業(yè)版,相機有效像素1 240萬,在新疆瓊哈拉峻村由軟件控制自動飛行獲取影像數(shù)據(jù)。為了便于對拍攝得到的影片進行拼接,保證拍攝時一定的航向與旁向重疊度(程效軍等,2005),設(shè)置旁向重疊度為75%,航向重疊度為50%,通過飛行拍攝得到69幅高分辨率航拍影片。
實際飛行時將無人機飛行高度(H)設(shè)置為相對地面100 m,由于無人機相機焦距(h=4 mm)和感光元件尺寸(d=0.001 58 mm)已知,參考公式:
式中:R為地面空間分辨率。無人機拍攝影像的地面空間分辨率約為3.9 cm,優(yōu)于航天平臺可見光遙感圖像的空間分辨率,達到高分辨率遙感影像的標準。
采用北京天地智繪科技有限公司開發(fā)的EasyUAV軟件對拍攝得到的影像進行拼接,該軟件采用的圖像拼接技術(shù)主要分為3個步驟:
(1)通過圖像預(yù)處理對圖像中幾何畸變進行矯正,使參考圖像和待拼接圖像不存在明顯的幾何畸變。本次飛行拍攝得到的影像為正射影像,其最常見的幾何畸變是由像主點或中心點向外擴展時光軸傾斜度增大、視野變寬、比例尺縮小所帶來的。由于本次飛行的飛行高度較低(100 m),單幅圖像的像主點與邊緣點相差的實際距離不大,在處理該幾何畸變時利用畸變的圖像提取陰影信息,沿著圖像失真的逆過程恢復(fù)圖像的原貌。實際的復(fù)原過程是設(shè)計一個反向濾波器,使其能從失真圖像中計算得到失真圖像的估值,根據(jù)規(guī)定的誤差準則,最大程度地接近真實圖像。預(yù)處理過程中還可以對圖像噪聲進行抑制,去除數(shù)據(jù)傳輸中帶來的不可預(yù)測的隨機信號,但這種隨機信號極易形成圖像上的噪聲點。噪聲點的去除通常由軟件在圖像輸出時完成,本次輸出的圖像分辨率較高,人為觀察沒有發(fā)現(xiàn)明顯的噪聲點,如果此時再用傳統(tǒng)的高斯濾波、均值濾波、中值濾波等低通濾波器對圖像進行二次處理,反而會使圖像的分辨率降低,降低提取精度。
(2)圖像配準主要對參考圖像和待拼接圖像中的匹配信息進行提取,在提取出的信息中尋找最佳的匹配,完成圖像之間的對齊。由于無人機拍攝影像分辨率較高,拍攝地區(qū)房屋數(shù)量較多,所以房屋的角點可以很好地作為配準點。使用逐一配準法對圖像進行配準,其原理為:在搜索圖S中以某點作為基點(搜索圖S的長寬分別為M,N),截取一個與模板T(長寬分別為U,V)大小一樣的分塊圖像,這樣的分塊圖像有(M-U+1)×(N-V+1)個,配準的目標就是在(M-U+1)×(N-V+1)個分塊圖像中找一個與待配準圖像最相似的圖像,這樣的基準點就是最佳配準點。
(3)在完成圖像配準后對圖像進行縫合,對縫合邊界進行平滑處理,讓縫合自然過渡,最終生成研究區(qū)DOM影像(圖2)。從圖2中可以看出, ???由于拍攝時間在冬季,新疆阿圖什市下轄瓊哈拉峻村內(nèi)農(nóng)田、草地等其他地物的光譜、紋理特征與房屋較為相似,光譜分辨率低極易出現(xiàn)同譜異物現(xiàn)象(田新光,2007;陳濟才等,2014),而這種情況可能會對接下來的分類產(chǎn)生一定影響。
2 研究方法
2.1 面向?qū)ο蟮挠跋穹治龇椒?/p>
面向?qū)ο笥跋穹治龇椒▽⑼|(zhì)性相近的臨近像元進行分割合并成影像對象,以影像對象作為分析處理的基本單元,對影像進行地物信息的提?。˙aatz,2000)。本文利用圖像中的光譜因子與形態(tài)因子計算該圖像中影像對象之間的異質(zhì)性(其中形態(tài)因子由光滑度和緊致度決定),計算公式為:
式中:F為異質(zhì)性;w為光譜權(quán)重;(1-w)為形狀權(quán)重。
光譜異質(zhì)性hcolor的計算公式為:
式中:c為圖層標識;wc為參與圖像分割的圖層權(quán)重;SMerge為合并后對象面積;σMergec為新生成的對象在圖層c中的光譜標準差;Sobj1,Sobj2分別為待合并的2個相鄰對象的面積;σobj1c,σobj2c分別為這2個對象在圖層c中的光譜標準差。
形狀異質(zhì)性hshape主要由緊致度hcompact和光滑度hsmooth2個方面組成,hcompact和hsmooth取決于構(gòu)成對象的像素個數(shù):
式中:wcompact為緊致度權(quán)重;(1-wcompact)為光滑度權(quán)重;nMerge,nobj1,nobj2分別表示合并對象的最小外包正方形、待合并的2個相鄰對象的最小外包正方形;lMerge,lobj1,lobj2分別表示這些對象的邊長。
對利用以上公式計算得到的異質(zhì)性與所給閾值進行對比,作為判斷某區(qū)域是否合并的標識。
2.2 面向像元的形態(tài)學(xué)處理技術(shù)
2.2.1 基于重建的膨脹與腐蝕方法
膨脹是求局部最大值的操作,對于前景物體中一些細小的斷裂處,如果小于結(jié)構(gòu)元素的大小,這些斷裂的地方就會被連接起來。而腐蝕與膨脹相反,是求局部最小值的操作。對于前景物體中一些細小的連接處,如果小于結(jié)構(gòu)元素的大小,這些連接處就會被斷開。膨脹與腐蝕的原理如圖3所示。A被B膨脹的數(shù)學(xué)表達式可記為AB,作為集合操作,定義為:
式中: 為空集;B為結(jié)構(gòu)元;B^表示B的反射,z表示把B的反射用z平移,可以表示為z(zx,zy),x和y表示平移的橫縱坐標,等式右邊可以解釋為B的反射被所有z平移后與A至少有一個非零公共元素。
A被B腐蝕的數(shù)學(xué)表達式可記為AB,作為集合操作,定義為:
等式右邊可以解釋為用z平移的B包含在A中的所有的點z的集合。
本文首先定義一種基于重建的膨脹方法,為下述的開閉操作進行準備。定義F為標記圖像,記錄膨脹的起點,G為模板圖像用來對膨脹進行約束且FG,則一次迭代的膨脹過程可以定義為:
則F關(guān)于G的n次迭代膨脹可定義為:
在該遞推公式中,集合求交在每一步都執(zhí)行,由于每一次迭代模板G會限制標記F的生長,所以經(jīng)過有限數(shù)量的迭代步驟后n總會收斂:
基于重建的膨脹算法原理如圖4所示?;谥亟ǖ母g原理與上述過程類似。定義F為標記圖像,記錄腐蝕的起點,G為模板圖像用來對膨脹進行約束且FG,則一次迭代的膨脹過程可以定義為:
式中:并集操作在每一次迭代中必須被執(zhí)行,保證每一次腐蝕過后的圖像仍然大于等于模板圖像,由于每一次迭代模板G會限制標記F的生長,所以經(jīng)過有限數(shù)量的迭代步驟后n總會收斂:
2.2.2 基于迭代的開操作與閉操作
在傳統(tǒng)的形態(tài)學(xué)開操作中,腐蝕典型地去除小的物體,且隨后的膨脹趨向于恢復(fù)保留的物體形狀。然而,這種恢復(fù)的精確度取決于形狀和結(jié)構(gòu)元之間的相似性。本文在基于重建的膨脹與腐蝕的基礎(chǔ)上,利用迭代的方式進行開操作能夠準確地恢復(fù)腐蝕之后的物體形狀。
A被B形態(tài)學(xué)開操作表示為A B,定義為A被B腐蝕,然后再用B膨脹腐蝕結(jié)果:
幾何學(xué)上,A B是B在A中完全匹配平移的并集。圖5a顯示了集合A和圓形結(jié)構(gòu)元B,圖5b顯示了在A內(nèi)完全匹配B的平移,所有這些平移的并集即圖5c為形態(tài)學(xué)開操作的結(jié)果,剩下被割去的區(qū)域就是結(jié)構(gòu)元不能完全在A中匹配的區(qū)域。形態(tài)學(xué)開操作除去了所有不能包含結(jié)構(gòu)元的部分,平滑目標的輪廓,斷開了細的連接部分,去除細的突出。A被B形態(tài)學(xué)閉操作表示為A·B,定義為A被B膨脹,然后再用B腐蝕膨脹結(jié)果:
幾何學(xué)上,A·B執(zhí)行所有不與A重疊的B平移的補集。圖5d顯示了一些與A不重疊的B的平移,通過完成以上平移的并集操作,可以得到圖5e所示的區(qū)域,該區(qū)域代表著完成閉操作后生成的區(qū)域。像開操作一樣,形態(tài)學(xué)閉操作趨向于平滑物體的輪廓。然而,不同于開操作的是,閉操作一般連接窄的斷裂并填滿細長的縫隙,填滿比結(jié)構(gòu)元小的洞。
本文所使用的迭代開操作可以定義為:記F為標記圖像,記錄腐蝕的起點,G為模板圖像用來對膨脹進行約束且FG,在每一次迭代過程中結(jié)構(gòu)元B先對模板進行腐蝕,之后再對腐蝕后的圖像進行膨脹重建??梢员硎緸椋?/p>
而對于閉操作,對上述開操作的結(jié)果求補即為閉操作結(jié)果。
3 實驗過程
3.1 面向?qū)ο蟮挠跋穹治龇椒?/p>
本文采用的圖像分割方法如下:基于前文所列公式計算得到每一次分割時每個影像對象的異質(zhì)性為F,利用基于異質(zhì)性最小原則的區(qū)域合并算法(即多尺度分割的思想)進行圖像的分割。首先將單個像元合并為一個個較小的對象,然后小的對象可以經(jīng)過若干步驟逐漸合并成大的對象,合并過程必須保證影像對象的異質(zhì)性小于給定的閾值。
在實際的分割過程中,由于房屋建筑的光譜特征差異性較小,設(shè)置光譜因子的權(quán)重為0.2,形狀因子的權(quán)重為0.8,由于目視觀察發(fā)現(xiàn)房屋屋頂紋理較為平滑,設(shè)置光滑度權(quán)重為0.7,緊密度權(quán)重0.3,經(jīng)多次嘗試閾值s設(shè)置在350 em為最佳。以任意一個像元為中心開始分割,第一次分割結(jié)束后由式(2)計算相鄰像元的異質(zhì)性F是否小于給定的閾值s,如果不滿足閾值條件,用上一次分割出的影像對象為基礎(chǔ)進行第二次分割,重復(fù)上述步驟直到相鄰2個影像對象的異質(zhì)性小于給定閾值s,完成最后的分割。利用多尺度分割的分割結(jié)果如圖6a所示。
由圖6a可以看出,在分割過程加入了緊致度以及表示屋頂?shù)墓饣纫蜃涌梢韵鄬τ行У貙ρ芯繀^(qū)不同地物進行分割,在計算異質(zhì)性的過程中弱化光譜信息在分割中所占權(quán)重,提高形狀因子所占權(quán)重,使得每一個影像對象的邊緣較為完整。對分割后的圖像提取每個影像對象3個波段的灰度均值以及標準差作為光譜特征,利用每個對象內(nèi)部的灰度共生矩陣計算ASM與IDM值,以1∶ 1加權(quán)求和作為紋理特征。ASM值反映了每一個對象內(nèi)灰度分布均勻程度和紋理粗細度,計算公式如下:
IDM反映圖像紋理的同質(zhì)性,度量圖像紋理局部變化的多少,計算公式如下:
利用提取出的特征向量對分割結(jié)果進行分類,如圖6b所示,其中綠色表示農(nóng)田,藍色表示房屋建筑,黃色表示道路,黑色表示陰影。
3.2 面向像元的形態(tài)學(xué)處理的算法實現(xiàn)
首先對原始圖像進行圖像增強處理,提取原始圖像每一點的灰度值做頻數(shù)直方圖,如圖7a所示,在圖中可以看出該圖像的灰度分布集中分布在50~255,0~50除了0灰度存在多數(shù)像元(圖像處理之后處于邊緣的無效值),其他灰度值幾乎無像元分布。針對這種情況采用分段線性拉伸的方法對圖像進行增強,將圖像中灰度50~255的區(qū)域進行拉伸。原始圖像中房屋的灰度值主要集中在250~255區(qū)域且有一個明顯的波峰。利用灰度拉伸后的原始影像作為底圖,提取250~255灰度區(qū)域的圖像作為前景圖像,對前景圖像二值化,結(jié)果如圖7b所示。
在前景圖像的提取過程中,容易提取出較多非房屋類型的像元,而用簡單的去噪方法無法有效地去除這種大面積噪音,還可能誤將有效的房屋像元一同去掉。鑒于這一問題,本文利用數(shù)學(xué)形態(tài)學(xué)圖像處理的方法去除圖像中復(fù)雜的背景,具體操作如下:①考慮到無關(guān)像元均為細碎斑點,較少數(shù)連續(xù)的無關(guān)像元面積也均遠小于房屋面積這一事實,采用長水平線作為結(jié)構(gòu)元B進行重建的開操作。在重建過程中,圖像首先被腐蝕,如同標準的形態(tài)學(xué)開操作一樣,然而這個被腐蝕的圖像在圖像重建過程中被用作標記圖像F,原始圖像作為模板G重新執(zhí)行形態(tài)學(xué)重建過程:將標記圖像F初始化為h1(標記F必須是G的子集,也就是FG)。建立結(jié)構(gòu)元B=ones(1,11)定義連通性,本文使用8連通。重復(fù)hk+1=(hkB)I G,直到hk+1=hk,完成重建的開操作。②對重建后的圖像執(zhí)行頂帽重建,即從原始圖像中減去重建的開操作,執(zhí)行結(jié)果記為f-thr。③將結(jié)構(gòu)元換成長豎直線B=ones(11,1),對原圖像進行頂帽重建,重建結(jié)果記為f-thi。④以f-thr作為模板,以min(f-thi,f-thr)作為標記,最后執(zhí)行形態(tài)學(xué)重建,作為最終的去噪結(jié)果。對重建之后的圖像用canny算子進行邊緣檢測,其中閾值向量T設(shè)定為0.04~0.1,平滑濾波器的標準差sigma設(shè)定為0.1。
對圖像中房屋的內(nèi)部孔洞、裂隙進行處理:本文基于重建的膨脹方法,利用迭代的方式進行開操作,準確地恢復(fù)腐蝕之后的物體形狀。在實際操作中,由于研究區(qū)內(nèi)的房屋多為矩形,且房屋面積均大于5 m2,因此在執(zhí)行開操作時選擇矩形結(jié)構(gòu)元且閾值設(shè)為5。對圖像進行迭代的開操作,將房屋中孔洞及裂縫進行填補,對房屋內(nèi)部的形態(tài)進行重建,之后同樣利用長寬分別為5的矩形結(jié)構(gòu)元對圖像進行迭代的形態(tài)學(xué)閉操作,對房屋外部輪廓進行重建,去除細的突出,使房屋邊界平整。最后對形態(tài)學(xué)重建后的圖像進行填充。整體算法的操作流程如圖8所示。加入形態(tài)學(xué)處理后的邊緣檢測結(jié)果與最終的提取結(jié)果如圖9所示。
3.3 精度評價
地表真實數(shù)據(jù)的獲取方法如下:以數(shù)字正射影像圖為底圖對研究區(qū)進行矢量化構(gòu)圖操作,利用目視判讀的方式提取出房屋所在位置以及面積信息,生成Shapefile文件,并以此文件作為模板在研究區(qū)域截取相關(guān)房屋數(shù)據(jù),對截取出的圖像進行二值化。為了對研究中提取的分類結(jié)果進行精度評價,將截取出的圖像進行重分類,將分類后的圖像(圖10)作為真實的房屋位置與面積數(shù)據(jù)。利用混淆矩陣對分類結(jié)果進行評價:分別將面向?qū)ο蟮姆课萏崛〗Y(jié)果與數(shù)學(xué)形態(tài)學(xué)方法進行房屋提取的結(jié)果作為分類數(shù)據(jù),將矢量化后的房屋位置及面積信息(圖10)作為地表真實數(shù)據(jù)建立混淆矩陣,計算制圖精度、用戶精度、漏分誤差以及錯分誤差,計算結(jié)果如表1所示。
4 討論與結(jié)論
在新疆瓊哈拉峻村內(nèi)通過小型無人機獲取影像數(shù)據(jù),分別利用面向?qū)ο蠛兔嫦蛳裨姆椒ㄟM行分類,從面向?qū)ο蟮挠跋穹诸惤Y(jié)果可以看出:面向?qū)ο蠓椒梢杂行У厝コ符}噪聲對分類帶來的影響;在分割過程中合理地加入地物的形態(tài)、紋理信息可以使分割后每一個對象的邊緣保持連續(xù)性,這也保證了影像分類時房屋形態(tài)的完整性。但是,在結(jié)果中出現(xiàn)了部分區(qū)域錯分的情況。觀察錯分區(qū)域的正射影像可以發(fā)現(xiàn):錯分區(qū)域多為農(nóng)田,該區(qū)域不僅光譜特征與房屋相似,還擁有與房屋屋頂相似的紋理以及矩形形狀。由于面向?qū)ο蟮姆诸愂腔诿恳粋€影像對象的分類,而在影像的多尺度分割時只能考慮相鄰影像對象之間的異質(zhì)性差異,不能考慮非相鄰對象之間可能會出現(xiàn)特征相似的現(xiàn)象。研究區(qū)農(nóng)田與房屋存在多數(shù)特征相似的情況,若使用面向?qū)ο蟮姆绞教崛∮跋駥ο髢?nèi)部特征進行分類,則容易導(dǎo)致影像對象的錯分。
利用面向像元的方法分離房屋則會存在提取的地物完整性不高、易出現(xiàn)椒鹽噪聲的問題。在實驗過程中首先需要考慮如何在不影響真實房屋像元的情況下剔除噪聲,且在房屋提取過程中應(yīng)盡量保持房屋的形態(tài)完整。本文在前人提出的數(shù)學(xué)形態(tài)學(xué)方法上做出改進,較好地實現(xiàn)了像元級房屋的提取,解決了面向?qū)ο蠓椒ㄖ谐霈F(xiàn)的由于影像對象內(nèi)部特征相似導(dǎo)致的部分農(nóng)田與房屋出現(xiàn)錯分現(xiàn)象這一問題。
參考文獻:
安基文,徐敬海,聶高眾,等.2015.高精度承災(zāi)體數(shù)據(jù)支撐的地震災(zāi)情快速評估[J].地震地質(zhì),37(4):1225-1241.
陳濟才,文學(xué)虎,李國明.2014.基于面向?qū)ο蟮母叻钟跋竦乇砀采w典型要素快速提取對比研究[J].遙感信息,(4):37-40.
陳晉,習(xí)聰望,陳文凱,等.2018.基于無人機,高分衛(wèi)星遙感影像的甘肅省隴南市建筑物空間化研究[J].地震研究,41(2):192-200.
程效軍,朱鯉,劉俊領(lǐng).2005.基于數(shù)字攝影測量技術(shù)的三維建模[J],同濟大學(xué)學(xué)報(自然科學(xué)版),33(1):37-41.
付博.2018.基于無人機正射影像的建筑物震害識別研究[D].北京:中國地震局地質(zhì)研究所.
李祖?zhèn)?,馬建文,張睿,等.2010.利用融合紋理與形態(tài)特征進行地震倒塌房屋信息自動提取[J].武漢大學(xué)學(xué)報(信息科學(xué)版),35(4):446-450.
聶高眾,安基文,鄧硯.2012.地震應(yīng)急災(zāi)情服務(wù)進展[J].地震地質(zhì),34(4):782-791.
聶高眾,高建國,馬宗晉,等.2002.中國未來10~15年地震災(zāi)害的風(fēng)險評估[J].自然災(zāi)害學(xué)報,11(1):68-73.
彭大雷,許強,董秀軍,等.2017.無人機低空攝影測量在黃土滑坡調(diào)查評估中的應(yīng)用[J].地球科學(xué)進展,32(3):319-330.
田新光.2007.面向?qū)ο蟾叻直媛蔬b感影像信息提取[D].北京:中國測繪科學(xué)研究院.
臧克,孫永華,李京,等.2010.微型無人機遙感系統(tǒng)在汶川地震中的應(yīng)用[J].自然災(zāi)害學(xué)報,19(3):162-166.
周洋,明小娜,楊艷珠,等.2017.災(zāi)評新技術(shù)在云龍5.0級地震烈度調(diào)查中的應(yīng)用[J].地震研究,40(1):161-166.