基于GF-2衛星影像的沖積扇地貌識別方法

周 煒，關洪軍

(陸軍工程大學 野戰工程學院， 南京 210007)

微地貌單元是重要的戰場地理要素之一，也是影響景觀格局的重要因素[1]。微地貌通常指規模相對較小的地貌形態，包含豐富的覆被屬性信息、水源信息、地質信息等，對分析軍事地質和特殊目標具有重要意義。近年來，微地貌的提取技術主要基于數字高程模型數據(Digital Elevation Model)實現。麻莉[2]基于DEM數據提出了一種針對洪泛平原濕地微地貌單元的數字分析方法。Baade等[3]應用TerraSAR-X雷達數據監測秘魯南部海岸沙漠地區微地貌變化，比較了利用DEM提取微地貌的局限性。沖洪積扇微地貌廣義上講包括沖積扇與洪積扇兩種。沖積扇沉積體具有構造“ 敏感性”，其復雜的沉積機制是解釋盆地邊緣構造活動的鑰匙，也是“源-聚-匯” 系統中的“源” 位置[4]。其前沿多為潛水溢出帶和淺埋地下水分布區，對野戰給水保障具有重要意義。目前，沖洪積扇微地貌提取技術多基于遙感影像的光學特征，DEM多用于輔助篩選泥石流分布區的溝谷和縱坡分析[5]，這些微地貌識別方法缺少地學因子的參考，不適應復雜多變的微地貌邊界的識別。本文以青藏高原西部，北緯32°26′，東經79°36′為研究區域，提出了一種基于CART分類器的沖洪積扇提取模型。該區域全年降雨量少，季節性強，氣候寒冷干燥[6]。寒冷、晝夜溫差大導致冰川發育廣泛，形成了高山、溝谷、冰蝕、冰磧、沖積扇等地貌類型，為沖洪積扇提取提供了大量的實驗對象。該模型以GF-2影像、Landsat 8影像、DEM影像為數據源，結合沖洪積扇微地貌的光學特征、幾何特征、紋理特征、地學特征位置特征，能夠準確界定沖洪積扇扇緣，實現沖洪積扇高精度提取。

1 數據處理

數據源包括Landsat8影像、GF-2衛星影像、數字高程模型。其中，數字高程模型根據1∶5萬地形圖生成；Landsat8影像數據來源于中國科學院計算機網絡信息中心地理空間數據云平臺(http://www.gscloud.cn)；GF-2影像為1 T級產品，參數如表1所示。

表1 GF-2衛星數據波段參數

數據處理包括幾何校正、濾波降噪、圖像融合3部分。幾何校正主要針對GF-2遙感影像、DEM影像、Landsat 8影像進行圖像配準，依靠重采樣對圖像亮度值進行插值計算，建立新的圖像矩陣，結果誤差控制在0.5個像元以內。濾波降噪通過中值濾波器實現，濾波窗口大小為5×5。中值濾波能夠有效濾除椒鹽噪聲，盡量保護邊緣信息的完整[7]。圖像融合的主要目的在于借助GF-2全色段影像的高分辨率優勢，提高洪積扇邊緣的識別精度。融合采用主成分分析法，融合效果如圖1所示。

2 沖洪積扇特征分析

沖洪積扇具有獨特的地學特征和影像特征，深入研究這些特征因子是構建沖洪積扇提取模型的必要前提。

1) 地學特征分析

① 沖洪積扇頂位于河流出山口，地學特征表現為高程的快速下降，位置關系表現為鄰接山體與平原。

② 沖洪積扇體在遙感影像上的幾何形態通常呈現為近似扇形。其發育范圍取決于山區河流的匯水面積、巖體風化特征、水文特征等。

③ 沖洪積扇體的物質組成取決于其河流上游匯水范圍的巖體類型；不同的物質成分在遙感影像上呈現不同的光譜特征。

④ 沖洪積扇體的土壤顆粒較粗，主要為砂、礫，通常植被稀疏，不利于植物生長。扇體邊緣和沖洪積扇前緣的土壤顆粒較細，一般為沙、粉沙及亞黏土，通常植被發育，植被覆蓋率高。

⑤ 沖洪積扇在地形上呈現為向盆地傾斜的山前傾斜平原或傾斜臺地，其坡度為5°～10°，扇體前緣的坡度為1°～2°。

2) 影像特征分析

沖洪積扇體的反射率隨波長變化的趨勢與沖洪積扇土體顆粒變化基本一致。由于扇體前緣土體顆粒變細，以沙、粉沙以及亞黏土為主，土體粒徑約在0.45～0.25 μm之間，土壤孔隙率相對增高，土壤含水率也明顯增高，其在各波段的反射率略低于沖洪積扇主體的反射率。沖洪積扇前緣地下水溢流帶產生的“微陰影”以及相對高的土壤含水率，導致扇體前緣部分的反射率、植被指數和地表溫度的突變，形成沖洪積扇前緣分界帶。

3 沖洪積扇提取模型構建

3.1 特征指數構建

1) 植被特征指數

歸一化差分植被指數[8](NDVI)常用于監測植被生長的狀態、植被覆蓋度，能夠有效排除植被密集區，計算公式如下：

(1)

式(1)中:R表示影像的紅波段，NIR表示影像的近紅外波段。

2) 溫度植被干旱指數

溫度植被干旱指數[9](TVDI)可用于干旱檢測以及研究干旱程度的空間變化特征，計算公式如下：

(2)

式中:Ts表示給定像元的地表溫度，Tsmin、Tsmax表示給定像元相同NDVI值下的最低溫度和最高溫度。根據Landsat 8熱紅外波段影像反演研究區溫度植被干旱指數，指數值越接近0，土壤濕度越大。反演結果如圖2所示。

3) 紋理特征指數

灰度共生矩陣是一種描述紋理特征的常用方法，是建立在估計圖像的二階組合條件概率密度函數基礎上的統計方法，能夠有效提取紋理特征[10]。通過灰度共生矩陣一定程度上能夠識別沖洪積扇頂的河道紋理。選用7×7窗口提取5種統計量，分別為一致性(Homogeneity)、對比度(Contrast)、相異性(Dissimilarity)、熵(Entropy)、均值(Mean)。其計算公式如下：

(3)

式(3)中:P(i,j)表示i行j列像元的亮度值，u為7×7窗口內的像元亮度均值。

4) 坡度、坡向

坡度(Slope)、坡向(Aspect)能充分呈現高程的變化率以及變化方向，有效識別沖洪積扇頂的地學特征。

(4)

式(4)中:m表示點在x方向上高程變化率，n表示點在y方向上高程變化率。坡度、坡向計算一般是在最小像元的基礎上，以5×5窗口為單位進行局部微分或曲面回歸擬合得出。

5) 光譜特征參數

沖洪積扇的物質組成主要以河流攜帶的沙礫石和粘土為主，分選性較好，植被極少。其在多光譜影像上色調比較均勻、亮度值變化差異較小。對沖洪積扇及該區域其他地物進行樣本采集分析，結果如表2所示。

表2 沖洪積扇及其他地物光譜值

3.2 分類器設計

決策樹是一種經典的分類器，其擴展通過節點的類別情況，每個類別導致一個結果[11]。分類器使用改進于傳統決策樹的CART(Classification And Regression Tree)分類方法，即分類與回歸樹[12]。

CART為在給定輸入x的條件下，輸出隨機變量y的條件概率分布。CART算法二分每個特征(包括標簽特征以及連續特征)，通過最優二分特征及其最優二分特征值的選擇、切分，生成二叉樹并進行剪枝[13]。分類過程中，回歸CART樹用誤差平方和準則選擇特征，分類CART樹用基尼系數準則選擇特征，并遞歸調用構建二叉樹過程生成CART樹。回歸樹實現流程如圖3所示。

3.3 沖洪積扇提取模型構建

1) 計算DEM影像的坡度值。以3°和10°為閾值，將坡度低于3°的區域視為沖洪積扇緣發育區；將坡度為3°～10°的區域視為沖洪積扇體發育區；將坡度高于10°的區域視為沖洪積扇頂發育區。

2) 計算目標區域多光譜影像的紋理特征指數。沖洪積扇頂的影像紋理特征與山體存在明顯區別，條帶狀河流沖刷痕跡明顯。為縮減灰度共生矩陣計算量，基于Matlab軟件平臺對影像的灰度級數進行壓縮，即從影像最小像元值起，每5級壓縮為1級，如1～5視為1，6～10視為2。

3) 計算NDVI、TDVI指數。研究區地處高原，多鹽堿地，植被指數主要與土壤類型相關，參考表2的平原臺地地光譜特征指數，以0作為區分閾值。根據沖洪積扇體表面植被覆蓋度極低的原則，將NDVI指數大于0的區域視為植被存在區域，排除在沖洪積扇可能存在的區域之外。從沖洪積扇根部體向扇體邊緣的土壤含水率逐漸增高，到扇前緣的地下水水溢出帶土壤含水率有明顯躍升，多泉水發育。根據TDVI指數，以0.3作為閾值提取影像中的高表層土壤含水率的區域。

4) 根據影像光譜特征參數，以藍波段像元值720為閾值，區分鹽堿地與其他地物，排除在沖洪積扇可能存在的區域之外。

5) 基于決策樹法排除有明顯特征的非沖洪積扇地物，界定山體與山間盆地的分界線。基于CART樹獲取最優分割特征值，實現不同源頭沖洪積扇的提取與區分。

6) 基于沖洪積扇與山體之間的鄰接關系判定沖洪積扇，實現沖洪積扇提取。

4 結果驗證與分析

為檢驗模型效果，分別使用本文提出的模型、最大似然法、SVM法對沖洪積扇進行提取，并通過混淆矩陣對提取結果進行精度評價，結果如圖4、表3所示。

注：圖中白色表示沖積扇A，灰色表示沖積扇B，黑色表示洪積扇C

表3 沖積扇分類結果

從表3看，本文提出的模型提取效果明顯優于最大似然法和SVM法，對影像中3種不同源頭的扇體提取結果的Kappa系數分別為0.95、0.86、0.90，總體精度達到92.37%。結合圖4所示，模型通過沖洪積扇與山體之間的拓撲鄰接關系，有效濾除了 “同譜異物”的其他背景地物，極大地提高了提取精度。

5 結論

1) 以Landsat8影像、GF-2衛星影像、數字高程模型為數據源，基于CART分類器的沖洪積扇提取模型，實現了沖洪積扇微地貌單元的精確提取，為區域地貌環境的監測提供技術支持。該模型提取沖積扇的總精度為92.37%。

2) DEM數據是沖洪積扇提取的重要參數，結合光學、地學特征的模型提取精度遠優于僅基于光學特征的最大似然法和支持向量機方法。

3) 溫度植被干旱指數能有效表達土壤濕度，對識別沖洪積扇頂部、前緣具有重要作用，能夠準確劃分沖洪積扇前緣地下水溢出帶。