董天成,楊肖,李卉,張志,齊睿
(1.中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,武漢 430074; 2.中國地質(zhì)大學(xué)(武漢)地球科學(xué)學(xué)院,武漢 430074; 3.32023部隊,大連 116023)
青藏高原分布有地球上海拔最高、數(shù)量最多、面積最大,以鹽湖和咸水湖集中為特色的高原內(nèi)陸湖群,是我國湖泊分布密集的地區(qū)之一,其湖泊總面積約占全國湖泊總面積的一半[1]。2019年7月亞洲水塔國際研討會上,中國科學(xué)家提出的以青藏高原冰川和湖泊為核心的第三極水塔計劃[2]受到廣泛關(guān)注,青藏高原冰川及湖泊變化情況已成為國際水資源研究與實踐探索的熱點之一。實現(xiàn)冰川和湖泊的準(zhǔn)確監(jiān)測是青藏高原水資源科學(xué)化管理的必要前提[3]。由于青藏高原區(qū)域自然環(huán)境惡劣,實地調(diào)查測繪等工作受制約情況明顯,利用遙感影像進行高原湖泊邊界及面積監(jiān)測成為目前最高效的手段之一。合成孔徑雷達(Synthetic Aperture Radar,SAR)主要反映地物的后向散射信息[4],使得以鏡面反射為主的湖泊呈現(xiàn)出完全不同于冰川和陸地等其他地物的特征。由于不受天氣、晝夜和季節(jié)的影響,SAR圖像已成為現(xiàn)階段湖泊勘測的重要數(shù)據(jù)源。
基于閾值分割的SAR圖像湖泊信息提取方法是國內(nèi)外較早一批發(fā)展完善的高效提取方法,2010年,安成錦等[5]利用典型Otsu算法閾值分割處理RadarSat-1影像,證明了多閾值SAR圖像水域分割的使用價值,最終分割結(jié)果精度達93.47%; 2013年,李智慧等[6]利用混合閾值法處理Envisat ASAR影像實現(xiàn)了陽澄湖輪廓提取,精度達99.41%,證明了其精度高于傳統(tǒng)單純一維或二維算法。值得注意的是,閾值分割方法簡單、運算速度快,能進行快速水體提取,但在區(qū)分與水體散射特征相似的地物時誤分情況較為明顯且缺乏空間特征考量。
為進一步提高湖泊提取精度,部分學(xué)者開始將SAR圖像中的水體紋理和形態(tài)特征引入用于辨識水體。Hahmann T等[7]利用Snake主動輪廓模型實現(xiàn)了TerraSAR-X數(shù)據(jù)中的湖泊分割,精度相較于閾值分割法提升明顯; 王慶等[8]將紋理特征與第一主成分閾值相結(jié)合,實現(xiàn)基于SAR數(shù)據(jù)的鄱陽湖水體變化監(jiān)測; Jin H等[9]和鄧瀅等[10]均提出結(jié)合紋理特征與極化特征進行極化SAR圖像水體識別,精度達98.89%,同一實驗區(qū)相較于基于像素的分類算法提高約6%; 冷英等[11]利用模糊C均值聚類法(Fuzzy C-Means,FCM)提取哨兵1A數(shù)據(jù)中的水體區(qū)域,實現(xiàn)速度提升明顯,精度可達97.60%; Sghaier M O等[12]提出局部紋理描述子與形態(tài)學(xué)相結(jié)合提取湖泊輪廓曲線的算法,利用RADARSAT-2圖像實現(xiàn)了加拿大多個地區(qū)湖泊提取,精度達到98.50%; Li N等[13]提出了噪聲過濾器與改進的幾何輪廓模型組合算法,基于高分三號和Sentinel-1A數(shù)據(jù)實現(xiàn)了丹江口水庫的水線提取與動態(tài)監(jiān)測。值得注意的是,上述多種方法還缺乏高級語義特征的加入,且實驗區(qū)差異性較大,各算法推廣性存在較大障礙。
由于對深層高維特征的挖掘能力和高應(yīng)用泛化能力,深度學(xué)習(xí)算法逐步成為時下基于SAR圖像湖泊提取的熱點算法。Zhang Y等[14]利用神經(jīng)網(wǎng)絡(luò)算法實現(xiàn)了基于高分三號影像的中國山西省岱海湖泊邊界提取,精度相較于傳統(tǒng)算法提升明顯; Shang F等[15]利用四元數(shù)神經(jīng)網(wǎng)絡(luò)實現(xiàn)了PolSAR圖像的復(fù)雜地表下湖區(qū)提??; Yang F Y等[16]利用Mask R-CNN實例分割算法[17]對包括SAR數(shù)據(jù)在內(nèi)的多種數(shù)據(jù)集進行訓(xùn)練,生成檢測模型在提取規(guī)則形狀的水體時準(zhǔn)確度達90%,提取不規(guī)則水體準(zhǔn)確度達76%,但在結(jié)論中發(fā)現(xiàn)Mask R-CNN在提取多水體樣本時檢測效果不佳。
上述各類算法均為低維或高維特征單一分析水體提取,已實現(xiàn)較高準(zhǔn)確度的SAR圖像湖泊提取,但當(dāng)進行大范圍多湖泊提取時,部分算法會受相似特征地物干擾,影響最終提取準(zhǔn)確度。實現(xiàn)高維特征與低維特征的結(jié)合會使分類精度繼續(xù)有提升的可能?;诖?,本文提出結(jié)合Faster R-CNN和MorphACWE的湖泊提取算法,使用歐洲航天局哥白尼計劃(GMES)中的Sentinel-1A衛(wèi)星的C波段合成孔徑雷達數(shù)據(jù)作為卷積神經(jīng)網(wǎng)絡(luò)的訓(xùn)練樣本和測試樣本數(shù)據(jù)源,通過Faster R-CNN的目標(biāo)檢測算法分析高維特征,實現(xiàn)大范圍區(qū)域內(nèi)的湖泊中心點確定,再通過MorphACWE輪廓模型結(jié)合低維湖泊紋理特征完成湖泊邊界劃分,綜合實現(xiàn)大范圍區(qū)湖泊端對端語義分割,在綜合干擾下多湖泊提取結(jié)果精度達99.71%。
研究區(qū)位于中國西藏自治區(qū)那曲市南部至日喀則市北部,如圖1所示,范圍為30°40′~ 33°01′N、86°37′~ 92°19′E。覆蓋面積約為17萬km2,平均海拔約4 600 m,擁有面積50 km2以上的湖泊38個,1 km2以上湖泊70多個。湖泊多為咸水湖,分布較為密集,同時又由于藏西地區(qū)植被稀少,較小時間跨度下湖泊變化情況不明顯,環(huán)境干擾因素較少。重點檢測區(qū)共包含1 km2以上湖泊6個,其中面積最大的為面積244.7 km2的達則措。
圖1 研究區(qū)概況示意圖Fig.1 Research area overview
實驗數(shù)據(jù)選擇歐空局Sentinel-1A干涉寬幅模式(TOPS Mode)的斜距單視復(fù)數(shù)產(chǎn)品(SLC),空間分辨率為5 m×20 m。共選擇研究區(qū)范圍內(nèi)VV和VH極化影像13幅,成像時間為2019年7月1日至2019年8月18日,數(shù)據(jù)列表如表1所示。上述數(shù)據(jù)均進行輻射定標(biāo)、地形校正、地理編碼和噪聲去除等預(yù)處理。研究區(qū)內(nèi)影像為目標(biāo)檢測模塊訓(xùn)練樣本集數(shù)據(jù)來源,重點檢測區(qū)內(nèi)影像為輸出模型檢測數(shù)據(jù)。
表1 研究使用數(shù)據(jù)Tab.1 Data in the study area
參考數(shù)據(jù)為歐空局Sentinel-2A的多光譜影像經(jīng)正射校正和亞像元級幾何精校正后的大氣表觀反射率產(chǎn)品(Level-1C)數(shù)據(jù),該影像共13個波段,空間分辨率為10 m。由于云霧影響,僅選取4幅影像,成像時間為2019年7月22日和2019年9月22日。完成輻射定標(biāo)和大氣校正等預(yù)處理后,對其進行人機交互解譯,分類解譯成果是后期進行多方法分類結(jié)果精度評價的重要參考數(shù)據(jù)。
2.1.1 基于VGG16的Faster R-CNN
在R-CNN和Fast R-CNN目標(biāo)檢測算法基礎(chǔ)上,Ren S Q等提出了Faster R-CNN算法[18]。該算法由兩個重要部分組成,分別是基于全卷積的區(qū)域生成網(wǎng)絡(luò) (region proposal networks,RPN)和基于區(qū)域的目標(biāo)識別網(wǎng)絡(luò)Fast R-CNN。兩個網(wǎng)絡(luò)結(jié)構(gòu)共享卷積層參數(shù),使得Faster R-CNN算法的檢測速度和精度都獲得大幅提高。
Faster R-CNN的工作流程分4步,如圖2所示,第一步將圖像輸入VGG16卷積神經(jīng)網(wǎng)絡(luò),提取出圖像的各項特征,輸出特征圖像; 第二步將特征圖像輸入RPN網(wǎng)絡(luò)中,通過Anchors(候選框種類k=9)輸出大量候選區(qū)域框; 第三步將第一步的特征圖像和第二步RPN輸出的候選區(qū)域框共同輸入Fast R-CNN框架; 第四步結(jié)合Softmax層的候選框中各類別概率分析和全連接層的候選區(qū)域精確校正,最終輸出目標(biāo)檢測結(jié)果。
圖2 Faster R-CNN結(jié)構(gòu)圖Fig.2 Structure of Faster R-CNN
2.1.2 改進型特征提取器
VGG16網(wǎng)絡(luò)的池化層位于卷積層后,主要目的為降低圖像維度,改善過擬合情況。原版特征提取器VGG16中所使用的多個最大池化層可以較好地保留圖像區(qū)域特征和紋理信息。但值得注意的是由于本文所需處理的影像數(shù)據(jù)為SAR數(shù)據(jù),該數(shù)據(jù)中極易出現(xiàn)大量噪聲信息,雖然對原始數(shù)據(jù)已進行LEE濾波等去噪處理,但離完全去除噪聲仍有部分距離,所以此時利用最大池化層進行池化操作可能會導(dǎo)致部分高噪聲區(qū)域的信息誤判,由于前向傳播的連續(xù)性,這一信息誤判會一直影響到最終結(jié)果的輸出; 此外池化層在增大感受野的同時會導(dǎo)致部分信息缺失,這對基于像元任務(wù)的結(jié)果預(yù)測有一定的影響[19]。基于上述問題,本文將原版特征提取器VGG16中的4個最大池化層更換為步長為2的卷積層,卷積核大小為3×3,如圖3所示。更換后的特征提取器可以更好地保留圖像的背景信息,也可以在一定程度上實現(xiàn)去噪功能,3×3的卷積核可以在保持原始數(shù)據(jù)的紋理信息基礎(chǔ)上減少特征損失。
圖3 改進型VGG16結(jié)構(gòu)圖Fig.3 Structure of Improved VGG16
無邊界主動輪廓模型(Active Contours without Edges,ACWE)是一種全局最優(yōu)的輪廓提取模型,于2001年由Chan和Vese提出[20],故又稱CV模型。CV模型假設(shè)圖像可以分為均勻同質(zhì)區(qū)域,利用目標(biāo)和背景的灰度值與各自對應(yīng)的平均灰度值差的平方來構(gòu)建驅(qū)動函數(shù)。
假設(shè)一幅圖像為I(x),則其無邊界主動輪廓能量泛函可表示為:
(1)
式中:C為輪廓曲線內(nèi)區(qū)域;c1和c2分別為輪廓線內(nèi)外部區(qū)域的灰度均值;μ為約束邊界長度的參數(shù); ν為約束閉合輪廓面積的常數(shù);λ1和λ2為控制輪廓線內(nèi)外能量權(quán)重的參數(shù); inside(C)為目標(biāo)區(qū)域,outside(C)為背景區(qū)域; 長度約束L和面積約束A用于控制輪廓線的光滑度和規(guī)則度。
公式右側(cè)第一和第二項均為輪廓曲線內(nèi)力,在迭代過程中可以控制過分割和誤分割產(chǎn)生的多余輪廓; 第三項和第四項為外力,用來引導(dǎo)曲線收斂到目標(biāo)邊界處。
當(dāng)F(C)取最小值時,輪廓線即為實際區(qū)域邊界。引入水平集可簡化極小化能量函數(shù),演化曲線用水平集函數(shù)u表示,同時加入一維Dirac測度函數(shù)及Heaviside函數(shù),可得:
(2)
值得注意的是原始CV模型需要使用數(shù)值積分方法來求解偏微分方程和水平集,計算成本相對較高,且存在一定穩(wěn)定性問題?;诖耍?014年Marquez-Neila P等[21]提出一種基于曲率形態(tài)算子的無邊界主動輪廓模型(MorphACWE) ,通過在二元水平集上定義一組具有等價無窮小行為的形態(tài)算子來實現(xiàn)偏微分方程的數(shù)值解。
假設(shè)第n次迭代后輪廓曲線變?yōu)閡n,則第n+1次函數(shù)變?yōu)椋?/p>
(3)
式中:Dd為擴張運算;Ed為收縮運算;SId°ISd為平滑操作中的曲率形態(tài)算子。
MorphACWE模型相較于原版CV模型運算速度更快,計算成本更低,穩(wěn)定性更強,相較于Snake模型和GAC模型則具有更好的抗噪性能[22],并且通過獲取全局最優(yōu)值可以更好地處理有模糊邊緣甚至離散狀邊界的目標(biāo),故成為現(xiàn)階段最為高效的輪廓提取算法之一。但CV系列模型的缺陷在于在對區(qū)域分割時分割結(jié)果過于依賴初始值[23],MorphACWE模型的初始值確定對分割精度和速度影響較大。
現(xiàn)有算法通過低維或高維特征單一分析水體提取,已實現(xiàn)較高準(zhǔn)確度的SAR圖像湖泊提取,但當(dāng)處理大范圍多湖泊提取時,部分算法會受相似特征地物干擾,影響最終提取準(zhǔn)確度; 此外嘗試將高維深層特征與低維特征結(jié)合可以使分類精度進一步提升。
基于上述問題,本文提出的FR-MorphACWE(Faster Region-based Convolution Neural Network-MorphACWE)算法以改進的Faster R-CNN目標(biāo)檢測算法和MorphACWE輪廓模型為基礎(chǔ)進行構(gòu)建。首先利用改進型Faster R-CNN對訓(xùn)練集數(shù)據(jù)進行訓(xùn)練,利用訓(xùn)練結(jié)果模型對檢測數(shù)據(jù)進行分類檢測,刪去原有算法中的識別輸出層,對接MorphACWE模型,直接將識別區(qū)域的中心點坐標(biāo)位置作為輪廓曲線膨脹初始點,再通過輪廓曲線的內(nèi)外力引導(dǎo)膨脹,最終迭代出最佳輪廓線,精準(zhǔn)劃分湖泊邊界,實現(xiàn)湖泊與非湖泊的語義分割,架構(gòu)如圖4所示。上述模型設(shè)計既可以保留原有Faster R-CNN目標(biāo)識別的準(zhǔn)確度,同時由于直接輸出識別區(qū)域中心點,避免了候選框與初始回歸位置的偏差; 此外利用MorphACWE模型提取輪廓可以較好地保證目標(biāo)的內(nèi)部同質(zhì)性和空間一致性。
圖4 FR-MorphACWE算法架構(gòu)圖Fig.4 Structure of FR-MorphACWE
為驗證本文方法有效性,具體實驗步驟(圖5)共分以下3步。第一,對研究區(qū)內(nèi)Sentinel-1A影像進行預(yù)處理,裁剪為小塊后建立L-SAR數(shù)據(jù)集進行檢測模型訓(xùn)練; 第二,利用引言中所提到SAR圖像湖泊提取發(fā)展各階段典型算法: OTSU閾值分割算法、模糊C均值分類、Mask R-CNN算法和本文提出FR-MorphACWE算法從綜合干擾下的多湖泊提取角度進行分類實驗; 第三通過計算準(zhǔn)確率、精準(zhǔn)率、召回率、F1分?jǐn)?shù)和Kappa系數(shù)實現(xiàn)湖泊提取效果定量評價。
圖5 實驗流程Fig.5 Experimental procedures
3.2.1 L-SAR數(shù)據(jù)集
樣本集是深度學(xué)習(xí)模型的血液。本文以PASCAL VOC數(shù)據(jù)集為藍本構(gòu)建高原湖泊SAR數(shù)據(jù)集(L-SAR數(shù)據(jù)集),數(shù)據(jù)源為13景完成預(yù)處理后的Sentinel-1A影像,選取大小湖泊共73個,利用圖像裁剪和數(shù)據(jù)增強處理共生成樣本數(shù)據(jù)共1 200景,訓(xùn)練湖泊累計1 872個,單景樣本數(shù)據(jù)尺寸在147像元×285像元到357像元×500像元范圍內(nèi)且保持每景中存在1~5個湖泊。為使最終輸出模型擁有更好魯棒性,L-SAR數(shù)據(jù)集包含不同極化方式、尺寸和區(qū)域位置等成像條件的樣本數(shù)據(jù)。將數(shù)據(jù)集中的960景作為Faster R-CNN的訓(xùn)練數(shù)據(jù)集(占比70%),240景作為測試數(shù)據(jù)集(占比20%),120景作為驗證數(shù)據(jù)集(占比10%)。
3.2.2 模型超參數(shù)及訓(xùn)練
本文模型基于Tensorflow框架,進行CUDA8.0加速后在GPU(NVIDIA GTX 1050ti 4 GB顯存)中實現(xiàn)檢測模型訓(xùn)練,具體參數(shù)如表2,訓(xùn)練時長為12 h25 min。模型損失率折線圖如圖6,損失率最終收斂至0.047 3~0.051 4。
表2 FR-MorphACWE超參數(shù)設(shè)置表Tab.2 FR-MorphACWE super parameter settings
圖6 損失率折線圖Fig.6 Loss rate line chart
多湖泊提取的主要問題在于多干擾因素的綜合影響,主要干擾因素包括高原冰川和湖心島,高原冰川多呈雪白色且整體表面光滑[24],在接收到微波傳感器發(fā)射的信號時后向散射極低; 湖心島是湖泊中露出水面的陸地,其上土壤水分含量相對高,后向散射系數(shù)也與水體近似[25]。此外研究區(qū)內(nèi)湖泊數(shù)量、大小差異性較大,形態(tài)復(fù)雜度和水質(zhì)變異度較高。受影像分辨率限制,研究范圍內(nèi)僅選取6個1 km2以上湖泊作為分類目標(biāo)進行綜合精度評價,其中最大的湖泊為達則錯,面積244.7 km2,面積最小的湖泊位于研究區(qū)中部偏東南方向,面積為4.95 km2。
(a) 強度值對比 (b) 湖泊影像1 (c) 湖泊影像2
對比各湖泊的強度值不難發(fā)現(xiàn),在選取各湖泊深水區(qū)的情況下,其強度值存在部分差異,其中湖泊4的強度峰值最小,整體湖泊平均值為0.004 806; 湖泊3的強度峰值最大,整體湖泊平均值為0.015 124; 同時存在冰川、湖心島等其他多種干擾地物影響下,各分類算法均出現(xiàn)小范圍誤分情況,各分類方法結(jié)果如表3所示; 各分類方法結(jié)果準(zhǔn)確率、精準(zhǔn)率、召回率、F1指標(biāo)和Kappa系數(shù)如表4所示; 各分類方法結(jié)果與真實值差異圖如圖8所示。
表3 各分類方法提取結(jié)果對比Tab.3 Comparison of extraction results of different classification methods
表4 各分類方法精度對比Tab.4 Accuracy comparison ofdifferent classification methods
總結(jié)表3、表4和圖8各項指標(biāo)可得,如表4中 OTSU閾值分割算法在處理大范圍且存在冰川、湖心島等干擾地物的區(qū)域時召回率僅為69.244 7%,大量干擾地物無法與湖泊區(qū)分出來。如表3研究區(qū)1號范圍內(nèi)的大量冰川和山體陰影等非湖泊區(qū)域被劃分為湖泊,但研究區(qū)2號范圍內(nèi)的湖心島整體識別情況較為良好。由于OTSU閾值分割算法是基于像素和閾值的分類算法,極易受到影像椒鹽噪聲的影響,分類結(jié)果湖泊內(nèi)部的同質(zhì)性受影響降低明顯,整體準(zhǔn)確率、精準(zhǔn)率和召回率相較于后續(xù)深度學(xué)習(xí)算法還存在明顯差距。
模糊C均值算法分類結(jié)果準(zhǔn)確率低于其他3種分類算法,召回率低于50%,研究區(qū)北部幾乎無法實現(xiàn)干擾排除,大量高原冰川被劃分為湖泊,且多湖泊內(nèi)部同質(zhì)性受噪聲影響明顯。主要原因為模糊C均值算法是一種無監(jiān)督聚類算法,通過迭代求解目標(biāo)函數(shù)最小化的過程來確定每一個樣本分類隸屬度,其對初始值和噪聲值非常敏感,迭代極易陷入局部極值, 導(dǎo)致偏離全局最優(yōu)。
Mask R-CNN實例分割算法是基于Faster R-CNN的重要改進算法,整體框架與Faster R-CNN保持一致,其創(chuàng)新性地加入了全連接的分割子網(wǎng),實現(xiàn)了目標(biāo)的實例分割。由于挖掘地物高維特征信息,結(jié)果準(zhǔn)確度相較于OTSU閾值分割算法和模糊C均值算法有大幅提升,湖泊內(nèi)部同質(zhì)性非常強。其準(zhǔn)確度高達98.79%,精準(zhǔn)率和召回率也均高于94%,證明了提升高維特征挖掘能力對提取精度的顯著影響。但從圖8中不難發(fā)現(xiàn),其分類結(jié)果中的湖泊邊界較為圓滑,與真實情況差距較大,未實現(xiàn)精準(zhǔn)湖泊邊界勾畫,湖泊內(nèi)湖心島也被誤分為湖泊。其原因主要由于分割子網(wǎng)利用反卷積操作來實現(xiàn)湖泊面提取,無法提供精度更高的邊界信息。
結(jié)合深度學(xué)習(xí)高維特征分析和邊緣輪廓模型目標(biāo)邊緣分割算法的FR-MorphACWE算法最大特點在于自適應(yīng)能力較強,在處理綜合干擾下的多湖泊提取情況下分類精度達到99%以上,且精準(zhǔn)率和召回率均保持在98%以上,均明顯高于其他3種分類算法。由于結(jié)合了MorphACWE輪廓模型的邊緣提取算法,湖泊外部輪廓清晰、內(nèi)部同質(zhì)性較強,不會受到影像椒鹽噪聲的影響。此外本文提出算法運行步驟較少,人為后續(xù)干預(yù)度低,優(yōu)勢更為明顯。
本文提出的FR-MorphACWE算法綜合改進型Faster R-CNN目標(biāo)檢測算法的高維特征分析和MorphACWE輪廓模型的邊界提取算法,實現(xiàn)了SAR圖像的高原湖泊端對端語義分割。本文在西藏自治區(qū)那曲市南部至日喀則市北部展開實驗,實現(xiàn)了多干擾、大尺度影響下的高精度湖泊提取,通過與OTSU閾值分割、模糊C均值和Mask R-CNN算法的對比可明顯發(fā)現(xiàn),本文提出的FR-MorphACWE算法整體參數(shù)設(shè)置少、需求數(shù)據(jù)量小、運行周期短、自適應(yīng)能力強、泛化性能好,可以作為SAR圖像湖泊提取的新算法加以推廣和應(yīng)用。
本文FR-MorphACWE算法的普適性還有待在中國地形二級階梯和三級階梯湖泊密集區(qū)做進一步檢驗,同時在下一階段的研究中可選擇以多源遙感數(shù)據(jù)作為數(shù)據(jù)源,在高準(zhǔn)確度提取的基礎(chǔ)上,實現(xiàn)多時相變化檢測。