郭天豪 解斐斐,2 霍志玲
(1. 山東科技大學(xué) 測繪科學(xué)與工程學(xué)院, 山東 青島 266590; 2. 山東省3S工程技術(shù)中心, 山東 青島 266590)
根據(jù)中國地震臺網(wǎng)數(shù)據(jù)顯示,世界上地震發(fā)生頻率非常頻繁,2019年4月至2020年4月僅強(qiáng)震(震級大于等于6.0)就發(fā)生了120余次。地震若發(fā)生在人口密度大、建筑物密集地區(qū),會嚴(yán)重危及當(dāng)?shù)厝嗣竦纳?cái)產(chǎn)安全[1]。在當(dāng)前技術(shù)手段下,預(yù)測地震難題依舊沒有解決。如何正確評估地震造成的地表破壞,特別是災(zāi)區(qū)建筑物破損統(tǒng)計(jì)分布情況,對政府相關(guān)部門制定救災(zāi)方案、最大限度地減少人員傷亡和財(cái)產(chǎn)損失等至關(guān)重要[2-3]。
合成孔徑雷達(dá)(Synthetic Aperture Radar,SAR)作為一種非接觸式測量手段,具有全天時、全天候地對地觀測優(yōu)勢,隨其不斷發(fā)展的SAR監(jiān)測手段,已成為監(jiān)測地震、火山等自然災(zāi)害的有效途徑[3-4]。為了提高地形形變精度,基于干涉原理的合成孔徑雷達(dá)干涉測量技術(shù)(Differential Interferometric Synthetic Aperture Radar, DInSAR)通過相位信息變化獲得形變信息,精度可達(dá)到厘米級甚至毫米級[5-7]。劉歡等基于成昆鐵路復(fù)線工程需要,利用DInSAR技術(shù)監(jiān)測垂直形變并驗(yàn)證其可行性與精度[8];焦昊分析了差分干涉方法,利用其監(jiān)測山西霍爾辛赫煤礦地形形變并提出了具體解決方案[9];Emanuela等基于DInSAR測量技術(shù)建立數(shù)學(xué)模型分析造成東阿塞拜疆地震的主要原因[10];Lorenzo等通過分析了250多篇論文綜述InSAR技術(shù)在意大利滑坡監(jiān)測研究方面的實(shí)際應(yīng)用[4],有力地論證了干涉技術(shù)在監(jiān)測地形形變的優(yōu)勢。
隨著高分辨率遙感技術(shù)的發(fā)展,利用高分辨率遙感影像為震害損毀建筑物提取提供了有利基礎(chǔ),成為當(dāng)前研究的熱點(diǎn)問題之一。然而,傳統(tǒng)的基于像元的方法由于“椒鹽噪聲”存在,解譯精度不高;面向?qū)ο蟮姆指罘椒ㄐ枰凑詹煌叨确指罱⒉煌某叨忍卣?自動化程度低,無法充分滿足抗震救災(zāi)、應(yīng)急指揮的要求[11]。近年來,基于機(jī)器學(xué)習(xí)的方法研究十分活躍,深度學(xué)習(xí)算法具有自學(xué)習(xí)和較強(qiáng)的容錯能力,被廣泛應(yīng)用于遙感圖像處理,為高效、自動化震區(qū)損毀建筑物提取提供可能性[12]。
本文利用覆蓋厄瓜多爾地區(qū)兩景Sentinel-1A數(shù)據(jù),采用差分干涉提取地震的形變場;基于高分辨率Worldview2光學(xué)衛(wèi)星數(shù)據(jù),利用ENVINet5深度卷積網(wǎng)絡(luò)對提取損壞建筑物,該研究為高效、自動化提取地震信息提供技術(shù)參考。
DInSAR是利用形變事件發(fā)生前后的兩期或多期合成孔徑雷達(dá)數(shù)據(jù),通過干涉圖差分處理來獲取地表微小形變的測量技術(shù)[13-14]。合成孔徑雷達(dá)生成的干涉條紋圖主要包含五項(xiàng)相位信息:地形相位、形變相位、大氣延遲相位、平地相位和相位噪聲,相位計(jì)算公式為[15-16]
φ=φdef+φtopo_res+φatm+φflat+φnoise
(1)
式中,φdef為視線向的形變相位;φtopo_res為地形相位;φatm為大氣延遲相位;φflat為平地相位;φnoise為噪聲相位。
(2)
式中,ΔR為形變信息;φp為兩幅SAR影像生成的干涉相位信息;φtp為地形相位信息。雙軌差分干涉原理如圖1所示。
圖1 雙軌法原理圖
卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Networks,CNN) 通過多層神經(jīng)元組合低層特征形成抽象的高層特征(屬性類別或特征),整合特征提取和分類器建模到一個學(xué)習(xí)框架下,減少人工/人為設(shè)計(jì)特征的工作。根據(jù)研究需要,本文選用ENVINet5卷積神經(jīng)網(wǎng)絡(luò)。
網(wǎng)絡(luò)的輸入是一張572×572的RGB影像,共有5個尺度,9個塊(block);編碼器結(jié)構(gòu)中的4個block均由3個有效卷積和一個最大值池化的下采樣組成,通過卷積和下采樣降低Feature map尺寸增加個數(shù),最終得到了32×32 尺寸的Feature map;解碼器結(jié)構(gòu)中的4個block(最后一個略有不同)由3個有效反卷積和1個上采樣組成,通過裁剪、合并相應(yīng)尺度的Feature map,最終得到了一個386×386 尺寸的類別激活柵格。類別激活柵格是一個灰度圖像,其像素大致代表給定特征的概率,圖像越明亮的地方匹配度越高[19]。
本文選取2016年4月厄瓜多爾地震作為研究對象,厄瓜多爾位于南美洲西北部,南接秘魯,西境瀕臨太平洋,國土面積約28萬km2。2016年4月16日18時53分37秒,厄瓜多爾地區(qū)發(fā)生7.8級地震,震源深度19.2 km,震源為0.371°N,79.940°W,造成668人死亡,超23萬人受傷。
本文的數(shù)據(jù)源有兩部分,一部分是用于地形變形監(jiān)測的Sentinel-1A雷達(dá)數(shù)據(jù)和數(shù)字高程模型(Digital elevation model,DEM)數(shù)據(jù),另一部分是進(jìn)行震害損毀建筑物提取的Worldview2光學(xué)衛(wèi)星數(shù)據(jù)(分辨率1.8 m)。其中雷達(dá)數(shù)據(jù)為2016年4月12日和24日兩幅厄瓜多爾區(qū)域Sentinel-1A 干涉寬幅模式斜距單視復(fù)數(shù)(IW SLC)衛(wèi)星數(shù)據(jù),極化方式為同向極化(VV),分辨率為5×20 m,具體參數(shù)如表1所示。DEM數(shù)據(jù)選用SRTM_V4 90 m DEM產(chǎn)品。
表1 SAR影像具體參數(shù)
本文采用雙軌差分干涉雷達(dá)技術(shù)對4月16日厄瓜多爾地震區(qū)域進(jìn)行震后形變場反演,具體流程如圖2所示。
圖2 形變反演流程圖
當(dāng)基線的垂直分量超出臨界線的極限時,相位信息不會被保留,相干性也將喪失,無法實(shí)現(xiàn)干涉測量,只有當(dāng)?shù)孛娣瓷渎示哂卸嘤趦蓚€天線重疊時才能產(chǎn)生干涉圖。
通過ENVI軟件進(jìn)行基線估計(jì)的運(yùn)算以獲得主從影像基礎(chǔ)信息及影像質(zhì)量,雷達(dá)波一次2π變化可以檢測到的形變變化為0.028 m,可以檢測到的地表形變敏感度為2.33 m。然后利用2期Sentinel-1A的干涉寬幅模式數(shù)據(jù),對距離相與方位相參數(shù)調(diào)整,結(jié)合DEM中的投影坐標(biāo)參數(shù),可以得到配準(zhǔn)后的干涉條紋圖。該干涉圖是由兩部分組成的:一是平地相位干涉條紋,由于兩次觀測的幾何關(guān)系不同,即便觀測的是平地,兩回波信號也會有一個相位差,形成干涉條紋;二是由實(shí)際觀測的地形產(chǎn)生的相位差。必須去掉前者,以得到純粹反映觀測地形的干涉圖,稱為去平地效應(yīng)。本文利用DEM數(shù)據(jù)對干涉相對去除平地效應(yīng)。
新形勢下計(jì)劃指標(biāo)使用管理的幾點(diǎn)思考(姜?dú)J杰) ........................................................................................1-42
由于三角函數(shù)的周期性,干涉相位的測量值都位于[-π,π],與真實(shí)相位相差2π的整數(shù)倍,相應(yīng)的相位圖呈條紋狀。為了將所測得的相位信息轉(zhuǎn)換成與傳播路徑直接相關(guān)的物理量,必須對干涉相位的測量值進(jìn)行二維相位展開處理。本文采用最小費(fèi)用流(Minimum Cost Flow)進(jìn)行相位解纏,該方法使用于當(dāng)有大面積的低相干或是其他限制增長的因素而使解纏困難時,解纏結(jié)果如圖3所示。為了實(shí)現(xiàn)控制測量,需要人工選擇控制點(diǎn),應(yīng)優(yōu)先在去平后的干涉圖上選擇控制點(diǎn),避免有地形相位沒有去除的區(qū)域和變化的區(qū)域(干涉條紋密集區(qū)域);盡量選擇相干性高的區(qū)域,遠(yuǎn)離形變區(qū),避免解纏錯誤的區(qū)域,不能位于解纏錯誤的相位躍變上(phase jump),如相位孤島等。本文利用三次多項(xiàng)式校正方法,選取26個點(diǎn),得到矯正后軌道均方根誤差為100.257。去除殘差后標(biāo)準(zhǔn)差為1.289,平均差為0.142。圖4為轉(zhuǎn)換到制圖坐標(biāo)系下的形變數(shù)據(jù),該形變數(shù)據(jù)表明厄瓜多爾地區(qū)整體地震形變在0~0.226 m之間。
圖3 相位解纏圖
圖4 形變結(jié)果圖
本文從DigitalGlobe提供的光學(xué)影像中選取震后特征明顯的建筑物區(qū)域,結(jié)合形變數(shù)據(jù)生成結(jié)果分析圖(圖5);從光學(xué)影像中可以看出,部分損壞建筑物特征明顯且大部分廢墟在形變變化處;形變大的區(qū)域地形損壞較為嚴(yán)重,建筑物坍塌也會導(dǎo)致救援困難,大形變區(qū)域和形變變化區(qū)域是損壞建筑物大概率存在處,也是災(zāi)后救援和財(cái)產(chǎn)評估的重點(diǎn)。
圖5 形變感興趣區(qū)結(jié)果圖
全損壞建筑物指研究區(qū)光學(xué)影像中所提取出來的全部損壞建筑物,廢墟是指已經(jīng)坍塌的建筑物。由于提取結(jié)果是概率柵格圖,因而需要選擇合適的概率值來區(qū)分不同損壞程度的建筑物,廢墟提取結(jié)果是在全損壞建筑物提取結(jié)果的基礎(chǔ)之上,通過目視解譯選擇合適的概率值篩選出廢墟,概率值68%以上為廢墟。實(shí)驗(yàn)數(shù)據(jù)選取DInSAR沉降形變大且房屋特征明顯的子影像區(qū)域,輸入原始影像和標(biāo)簽數(shù)據(jù)訓(xùn)練分類模型,調(diào)整網(wǎng)絡(luò)參數(shù),主要有損失權(quán)重、疊代次數(shù)、采樣密度等[19],提取流程如圖6所示。
圖6 提取損毀建筑物流程圖
基于ENVINet5的深度學(xué)習(xí)框架進(jìn)行研究區(qū)建筑物提取和損壞建筑物提取。結(jié)果表明:提取較大廢墟效果較好,對于損壞程度較小的建筑物只能提取小部分。因?yàn)槠淇蚣苁歉鶕?jù)感興趣區(qū)特征符合程度分類,特征越明顯,概率越大,而損壞程度較小的建筑物特征不明顯,概率就比較小,提取結(jié)果如圖7~8所示。
圖7 廢墟提取結(jié)果圖
圖8 損壞建筑物提取結(jié)果圖
研究區(qū)光學(xué)影像和形變結(jié)果的影像空間分辨率均為1.8 m,在經(jīng)過多次余震后,老房或舊房等建筑物在形變量較低的情況下,極易發(fā)生損壞情況。由圖9可知,對比光學(xué)影像和形變場結(jié)果影像的相同區(qū)域,損壞建筑物分布在形變的變化交界處附近,符合實(shí)際情況。
圖9 損壞建筑物和形變場對比圖
研究區(qū)損壞建筑物樣本數(shù)為27個,正確分類數(shù)為20個;建筑物樣本為79個,正確分類數(shù)為63個。由此可知,通過ENVINet5提取的損壞建筑物精度達(dá)到74%,建筑物精度達(dá)到79%。說明基于ENVINet5的深度學(xué)習(xí)網(wǎng)絡(luò)框架對建筑物進(jìn)行提取精度較高,為震后受災(zāi)情況分析提供了良好的科學(xué)依據(jù),具有較強(qiáng)的可行性。
本文以厄瓜多爾地區(qū)為研究區(qū)域,基于Sentinel-1A雷達(dá)數(shù)據(jù),采用雙軌干涉測量技術(shù),反演厄瓜多爾地震形變場;利用ENVI5提供的Deep Learning 模塊的ENVINet5深度卷積網(wǎng)絡(luò)對Worldview2光學(xué)衛(wèi)星數(shù)據(jù)提取震后損壞建筑物。DInSAR形變結(jié)果顯示,厄瓜多爾地區(qū)整體形變值在0~0.226 m,越靠近震中,地表凹陷程度越大;遠(yuǎn)離震中,地表形變趨于0。由損壞建筑物分布可以看出,在形變值變化區(qū)域有80%以上概率存在損壞度較大建筑物,地形形變明顯,道路阻塞嚴(yán)重,該區(qū)域應(yīng)該成為后續(xù)救援的重點(diǎn)。利用深度卷積神經(jīng)網(wǎng)絡(luò)ENVINet5提取損壞建筑物的過程為概率提取,相似程度越大,其概率越高。本文所選感興趣區(qū)提取分類精度達(dá)74%,這說明基于深度卷積網(wǎng)絡(luò)在高分辨率光學(xué)影像上提取災(zāi)后建筑物,精度高,效率高,節(jié)省人力、物力,有較好的可行性,對災(zāi)后房屋評估、災(zāi)區(qū)重建等提供了科學(xué)依據(jù)。