基于深度學習和高分遙感影像的青藏高原東部融凍泥流階地提取研究

王習敏， 黃榮剛， 徐志達， 焦志平， 江利明

（1.中國科學院精密測量科學與技術創(chuàng)新研究院大地測量與地球動力學國家重點實驗室，湖北 武漢 430077；2.中國科學院大學地球與行星科學學院，北京 100049）

0 引言

融凍泥流階地是一種冰緣區(qū)域斜坡受到反復凍融和重力作用形成的臺階狀堆積微地貌［1］。全新世以來青藏高原邊緣地帶（如東部地區(qū)）近地表經(jīng)歷了多次凍融循環(huán)的復雜歷史演化過程［2-3］，古融凍泥流階地分布廣泛，其空間分布對于重建第四紀古氣候環(huán)境和多年凍土分布范圍均具有重要意義［4-6］。

目前，針對青藏高原融凍泥流階地觀測較少，大范圍制圖研究基本空白。在先前的研究中，受到衛(wèi)星遙感影像分辨率過低的限制，實地調(diào)繪和航空遙感是融凍泥流階地制圖的兩種主要方式［7-11］。例如，1984—1990年，郭東信等［11］對青藏公路風火山埡口盆地融凍泥流階地先后進行了兩次地面調(diào)繪，現(xiàn)場測量還在很多泥流階地站點開展過實驗，如格陵蘭東北部［12］、斯瓦爾巴群島［13］、瑞典阿比斯庫［14］、科羅拉多前沿山脈［15］等。自20世紀90年代開始，航空遙感與攝影測量技術逐漸用于泥流階地提取和制圖。Walsh等［10］使用ADAR-5500航空影像繪制了美國冰川國家公園的泥流階地；Ridefelt等［7-8］利用航空攝影測量手段在瑞典阿比斯庫山進行了多期融凍泥流階地識別與分析。結合實地勘測調(diào)查和航空影像解譯，陳游東［16］獲取了中天山典型地區(qū)融凍泥流階地的分布信息，并分析了其空間分布規(guī)律及對公路工程的災害性影響。但是，以上手段只針對地塊和小范圍的泥流階地，難以滿足大范圍制圖的需求。

近年來，隨著遙感技術的發(fā)展，高分辨率衛(wèi)星遙感技術在融凍泥流階地大范圍識別中具有良好的應用潛力。然而，受融凍泥流階地幾何形態(tài)不一、表面覆蓋多樣（如裸地、草甸和稀疏植被）等因素影響，融凍泥流階地紋理復雜，基于遙感影像的自動提取依然面臨很大的挑戰(zhàn)。快速發(fā)展的深度學習方法，如U-Net［17］，SegNet［18］和DeepLab［19］，在遙感影像地物分類和目標識別等方面取得了突破性的進展，在地球科學和遙感領域得到了廣泛的應用［20-23］。其中，DeepLab V3+［24］深度學習模型，能獲取多尺度上下文信息和清晰的邊界，已用于青藏高原熱融滑塌制圖［25-26］，有望為復雜環(huán)境下基于遙感影像的大范圍泥流階地自動提取提供一種新途徑。

因此，本文提出了一種基于DeepLab V3+和高分辨率遙感影像的融凍泥流階地提取方法，并選擇甘孜州新都橋周邊為研究區(qū)，在青藏高原開展融凍泥流階地的遙感自動識別研究。首先，利用研究區(qū)Google Earth高分辨率遙感影像，通過人工勾勒方式繪制了泥流階地幾何邊界，建立深度學習的樣本數(shù)據(jù)集；然后，構建DeepLab V3+泥流階地深度學習模型，并將訓練好的模型用于高分辨率遙感影像進行泥流階地自動提取；最后，對提取結果進行精度評定，并分析了泥流階地空間分布特征及其與地形因子的關系。本研究為融凍泥流階地的大范圍遙感制圖提供了新思路。

1 研究區(qū)與數(shù)據(jù)

研究區(qū)位于四川省甘孜藏族自治州，如圖1，總面積為900 km2。該區(qū)域平均海拔3 945 m；年降水量約553 mm，降水量主要集中于6—7月，日降雨量最高達到50 mm；溫度范圍為-15～28℃，最高溫度分布于7月，最低溫度主要在1月和12月。該區(qū)域地表覆蓋主要是草地、少量灌木和農(nóng)田等。區(qū)域內(nèi)融凍泥流階地分布較多，主要集中在318國道、215省道沿線。

圖1 研究區(qū)位置Fig.1 Location of the study area

本文選用Google Earth遙感影像作為實驗數(shù)據(jù)，空間分辨率為0.59 m，數(shù)據(jù)采集時間為2016—2019年。此 外，選 用2008年 的12.5 m分 辨 率ALOS PALSAR DEM數(shù)據(jù)，用于泥流階地的地形因子分析。2020年30 m分辨率的地表覆蓋數(shù)據(jù)用于分析融凍泥流階地與地表下墊面的關系，數(shù)據(jù)來源于國家基礎地理信息中心全球地表覆蓋數(shù)據(jù)產(chǎn)品服務網(wǎng)站。

2 研究方法

2.1 DeepLab V3+模型

本文采用最近提出的DeepLab V3+網(wǎng)絡，開展融凍泥流階地遙感自動提取，該深度學習語義分割模型將深度卷積神經(jīng)網(wǎng)絡（DCNN）與概率圖模型結合，取得了比DCNN更好的檢測精度，并將在空間金字塔池化結構（ASPP）和解碼器部分（Encoder）引入深度可分離卷積降低網(wǎng)絡的計算復雜度［27］，提高了模型學習的效率，已在遙感圖像分割和目標識別中取得了廣泛應用［27-28］。

DeepLab V3+網(wǎng)絡結構如圖2所示，采用下采樣編碼和上采樣解碼結構。在編碼階段，通過骨干網(wǎng)絡，如Xception等［29］，結合多尺度空洞卷積，實現(xiàn)多尺度深度特征感知。在解碼階段，通過融合骨干網(wǎng)絡低層次特征和空洞卷積高層次特征，并采用簡單的上采樣操作，逐步提高語義分割的分辨率，實現(xiàn)與輸入圖像一致大小的分割結果。其流程主要包括：采用骨干網(wǎng)絡（如Xception65［29］）進行深度可分離卷積運算，然后通過ASPP模塊獲取多尺度特征（空洞率分別為1，6，12，18），并將ASPP空洞卷積層特征與全局池化層進行特征的融合，以此得到多尺度的高級語義特征輸出結果；將從主干網(wǎng)絡中提取的低級語義特征和編碼階段提取的高層次多尺度特征圖進行融合，并采用3×3的卷積運算和4倍上采樣，提高其空間分辨率。

圖2 DeepLab V3+網(wǎng)絡結構Fig.2 Structure of DeepLab V3+

2.2 基于DeepLab V3+的泥流階地自動提取

基于DeepLab V3+的泥流階地自動提取處理流程如圖3所示。首先，訓練樣本準備，制作泥流階地數(shù)據(jù)集。然后，設置參數(shù)，對深度學習網(wǎng)絡Deep-Lab V3+進行訓練。其次，利用訓練好的模型對研究區(qū)域Google Earth高分辨率影像進行預測、后處理。最后，對提取的泥流階地進行精度評估。

圖3 基于DeepLab V3+的泥流階地提取流程圖Fig.3 Flowchart of solifluction terraces extraction based on DeepLab V3+

2.2.1 多類型泥流階地正負樣本制作

綜合泥流階地明暗相間的紋理特征、長條狀幾何形態(tài)、集群分布等特點，結合Google Earth三維地形分析，采用ArcGIS矢量化工具勾勒泥流階地前沿陡坎邊緣線作為樣本。在區(qū)域A影像上勾勒泥流階地陡坎邊界，共計4 752條，主要包括高陡坎強紋理的泥流階地、低陡坎弱紋理的泥流階地、植被覆蓋的泥流階地和裸地覆蓋的泥流階地4種類型，如圖4（a）～4（d）所示。此外，本文還勾勒出易混淆的其他目標作為負樣本，如山脊線、山谷線、田地間陡坎、沖溝等線狀目標區(qū)域，如圖4（e）～4（h）所示，共計110個多邊形。最終勾勒的泥流階地和非泥流階地樣本分布如圖1所示。基于人工勾勒的正負樣本矢量多邊形，結合對應的Google Earth高分辨率影像，生成大小為513×513的正負樣本切片數(shù)據(jù)集，所得正負樣本切片數(shù)量分別為3 046和1 863。

圖4 泥流階地人工勾勒正負樣本示例Fig.4 Examples of manually delineated positive and negative samples：solifluction terraces with high risers and strong textures（a）；solifluction terraces with low risers and weak textures（b）；solifluction terraces covered by bare soil（c）；solifluction terraces covered by vegetation（d）；ridge lines（e）；valley lines（f）；raised paths through farmlands（g）；gullies（h）

2.2.2 DeepLab V3+模型訓練與泥流階地提取

基于生成的泥流階地數(shù)據(jù)集，對DeepLab V3+模型進行迭代訓練，獲取其模型參數(shù)。在訓練過程中，參數(shù)設置如下：批處理量為12；迭代次數(shù)為6×104次；初始學習率為1×10-3；權重衰減系數(shù)設置為4×10-5；采用的優(yōu)化器是Adam；骨干網(wǎng)絡為Xception65；ASPP結構中空洞率分別為［1，6，12，18］，ASPP模塊可有效提取不同尺度下的語義特征，一定程度上有助于提高不同尺度泥流階地的提取效果。

在模型訓練好后，將其用于整個實驗區(qū)的融凍泥流階地預測提取。首先，將實驗區(qū)所有Google Earth高分辨率影像切割成513×513大小的切片；然后，利用訓練好的模型對各切片進行泥流階地提取，將切片分割成泥流階地陡坎和其他背景像素。最后，使用GDAL進行鑲嵌處理，輸出完整的預測結果。

為了剔除錯誤提取的泥流階地，結合幾何特征、集群存在特征、與等高線相對平行的走勢關系等先驗知識，對預測結果進行后處理。后處理的規(guī)則主要包括：（1）剔除面積小于40 m2、長度小于10 m、寬度小于4 m的提取目標；（2）刪除條數(shù)少于10條的零散分布提取目標；（3）去除個別與等高線走勢相悖的目標。

2.3 精度評估

為評估泥流階地提取方法的性能，采用人工勾繪的結果作為參考，分別根據(jù)公式（1）、（2）和（3）計算以下三個統(tǒng)計指標對提取結果進行精度評估［30］。

式中：P（Precision）和R（Recall）分別為精確度、召回率，F(xiàn)1是綜合兩者的度量指標，TP、TP和FN分別正確提取、錯誤提取、漏提取的泥流階地陡坎面積。

3 結果與分析

3.1 泥流階地提取結果

圖5為基于深度學習的融凍泥流階地深度學習提取結果，共9 203條，總面積為2.12 km2。據(jù)提取結果發(fā)現(xiàn)，在本研究區(qū)內(nèi)，泥流階地沿山谷線兩側分布，主要位于國道等不同等級道路兩側，呈現(xiàn)了聚集分布。

圖5 基于DeepLab V3+模型的泥流階地自動提取結果Fig.5 Results of solifluction terraces based on DeepLab V3+

對于不同覆蓋類型具有明顯特征的泥流階地，本文方法在目標提取率和邊界吻合程度等方面均取得很好的效果，泥流階地的邊界提取效果較好，但類型間也存在差別。如圖6所示，植被覆蓋的泥流階地，條狀形態(tài)、紋理不大明顯的部分，提取結果有待提高；陡坎較高的泥流階地，邊界位置相對精準；陡坎紋理不清的泥流階地，深度學習模型僅能提取到部分結果。此外，相比較于形態(tài)尺度較小的泥流階地陡坎，形態(tài)尺度較大的泥流階地陡坎提取效果更好。

為了定量評估泥流階地深度學習提取質(zhì)量，本文分別統(tǒng)計了融凍泥流階地正確提取、錯誤提取以及漏提取的面積，并根據(jù)公式（1）、（2）和（3）計算深度學習自動提取結果的精度、召回率和F1（表1）。在訓練區(qū)域A范圍內(nèi)，精度F1達到0.738；在預測區(qū)域B范圍內(nèi)，精度F1達到0.68。綜合圖5圖6和表1看，本文方法總體提取精度較高，但也存在較多的錯檢和漏檢，主要原因有：（1）算法容易錯把少量梯田、植被、沖溝當成泥流階地陡坎；（2）由于樣本有限以及自動后處理面積閾值剔除等情況，在小面積泥流階地、弱紋理泥流階地（特別草甸覆蓋區(qū)域）等區(qū)域存在漏提取。

圖6 泥流階地提取結果與人工判讀比對Fig.6 Comparison between automated extraction solifluction terraces and the interpretation results：solifluction terraces covered by vegetation（a）；solifluction terraces with high risers（b）；solifluction terraces covered by meadow（c）；solifluction terraces covered by bare soil（d）

表1 融凍泥流階地提取結果精度評估Table 1 Evaluation of extracted result of solifluction terraces

綜上，泥流階地提取方法性能分析表明，Deep-Lab V3+深度學習模型能夠有效提取研究區(qū)泥流階地陡坎的邊界，且擁有良好的泛化能力，在形態(tài)、紋理相似的泥流階地區(qū)域（如G318國道沿線）具有推廣應用潛力。

3.2 實地考察驗證

為進一步驗證提取的結果，2020年夏季進行了實地考察［圖7（a）］，考察點位于G318國道沿線，如圖1中黑色定點標記所示。該處泥流階地陡坎特征不明顯，對其提取具有較大難度。但是通過自動提取結果和人工解譯結果對比［圖7（b）］表明，深度學習提取方法能夠較好地提取大部分泥流階地，但由于部分泥流階地陡坎特征紋理太弱，導致部分泥流階地無法有效提取［圖7（b）］。

圖7 提取結果的實地驗證Fig.7 Field validation of extraction results：on-site photos（a）；results comparison（b）

4 討論

4.1 地形因子統(tǒng)計分析

為分析泥流階地的空間分布特征，本文討論了研究區(qū)內(nèi)泥流階地與地形因子的統(tǒng)計關系。圖8為區(qū)域A、區(qū)域B的地形因子（坡度、坡向、高程）統(tǒng)計圖。結果顯示：

（1）在研究區(qū)A內(nèi)，泥流階地坡向主要集中在西北方向［圖8（a1）］。在高程統(tǒng)計方面，平均值為3 614.420 m，標準差為81.483 m，且主要集中在3 550～3 650 m區(qū)間，占研究區(qū)A內(nèi)泥流階地的43%

［圖8（b1）］。其中，最高的泥流階地，其高程為3 875 m，位于318國道沿線的水橋村北側920 m附近；最低的泥流階地高程值為3 438 m，位于營官寨村西北部2.5×103m附近。泥流階地坡度主要集中分布在10°～30°區(qū)間，占全部泥流階地數(shù)量的79%。

（2）在研究區(qū)B內(nèi)，泥流階地坡向主要集中在正西、西北和東北等陰坡方向，與研究區(qū)A內(nèi)的泥流階地具有類似的坡向分布特征［圖8（a2）］。在高程分布方面，泥流階地主要集中在3 650～3 750 m分布范圍內(nèi)，占研究區(qū)B內(nèi)泥流階地的43%；區(qū)域內(nèi)泥流階地平均高程為3 674.694 m，標準差為86.240 m。在坡度上，集中分布在10°～30°，并且在20°～25°區(qū)間的泥流階地數(shù)量最多。

圖8 泥流階地的地形因子統(tǒng)計分析Fig.8 Statistical analysis of topographic factors of solifluction terraces：aspect（a），elevation（b），slope（c）of Area A and B

總體而言，本研究區(qū)域內(nèi)泥流階地以西北、北等陰坡方向為主；在坡度上，集中分布在10°～30°之間；在高程上，3 650～3 750 m區(qū)間為主要分布范圍。

4.2 幾何因子統(tǒng)計分析

此外，本文還對泥流階地陡坎的幾何因子（面積、周長、長度、寬度）進行了統(tǒng)計分析和討論，如圖9所示。統(tǒng)計結果顯示：

圖9 泥流階地幾何因子統(tǒng)計Fig.9 Statistics of geometric factors of solifluction terraces：area（a）；perimeter（b）；length（c）；width（d）

（1）在研究區(qū)A中，泥流階地面積平均值為369 m2，標準差為330.79 m2。絕大部分泥流階地面積小于1 000 m2，占總數(shù)的95%。其中，小于500 m2的泥流階地占79%。泥流階地周長平均為128.46 m，標準差為79.59 m，最大值和最小值分別為890.65 m和20.50 m。長度、寬度大部分小于500 m，均值分別為54 m、26 m。

（2）在研究區(qū)B中，泥流階地面積均值463 m2，標準差為412.46 m2，小于1 000 m2的占91%；周長平均為146.84 m，標準差為86.36 m，100～200 m區(qū)間最為集中；長度、寬度均值分別為62 m、29 m。

總體而言，研究區(qū)A和研究區(qū)B內(nèi)的泥流階地在幾何形態(tài)上具有較高的相似性，但也存在一定差異。例如，研究區(qū)B的平均長度和寬度均大于研究區(qū)A。

4.3 地表覆蓋

研究區(qū)域的現(xiàn)今地表覆蓋空間分布如圖10所示，主要有草地、森林、耕地、水體、人造地表、積雪和冰川共6種。對各類型地表覆蓋的泥流階地進行統(tǒng)計分析，如圖11所示。結果表明：泥流階地主要分布在草地、森林、耕地三種地表覆蓋類型處。區(qū)域A中的草地類型占94%；區(qū)域B中的草地類型占95%。

圖10 泥流階地地表覆蓋空間分布圖Fig.10 Spatial distribution map of surface cover on solifluction terraces

圖11 泥流階地地表覆蓋類型統(tǒng)計Fig.11 Statistics on surface coverage types of solifluction terraces

5 結論

基于DeepLab V3+深度學習網(wǎng)絡和Google Earth高分辨率衛(wèi)星影像，本文提出了一種大范圍融凍泥流階地遙感自動提取的方法，并在四川省甘孜州新都橋地區(qū)開展了實驗研究，分析了該區(qū)泥流階地分布特征及其與地形因子、幾何因子、地表覆蓋等因素的關系。主要結論如下：

（1）與人工解譯結果對比，本方法提取結果的綜合精度達到0.68以上，并經(jīng)野外調(diào)查驗證了其有效性。陡坎較高的泥流階地提取精度最高，對于弱紋理的泥流階地精度較低但能實現(xiàn)其定位。

（2）共識別了9 203條融凍泥流階地，主要分布在新都橋鎮(zhèn)附近的山谷兩側，海拔高程主要分布在3 650～3 750 m區(qū)間，且主要位于20°～25°的陰坡上。

（3）泥流階地面積和周長分別在1 000 m2、150 m以內(nèi)，地表覆蓋類型主要是草地，占比高達95%。

致謝：本研究使用的數(shù)據(jù)來自于Google Earth平臺提供的遙感影像數(shù)據(jù)、日本JAXA提供的ALOS PALSAR DEM數(shù)據(jù)、國家基礎地理信息中心全球地表覆蓋數(shù)據(jù)產(chǎn)品中心提供的地表覆蓋數(shù)據(jù)、第二次青藏高原綜合科學考察項目提供的野外照片資料，本文訓練樣本制作中武漢工程大學王慧妮老師給予了大力支持，中國科學院精密測量科學與技術創(chuàng)新研究院陳圓圓同學提供了繪圖調(diào)色指導，在此表示衷心感謝。