周 煒,關(guān)洪軍
(陸軍工程大學(xué) 野戰(zhàn)工程學(xué)院, 南京 210007)
微地貌單元是重要的戰(zhàn)場地理要素之一,也是影響景觀格局的重要因素[1]。微地貌通常指規(guī)模相對較小的地貌形態(tài),包含豐富的覆被屬性信息、水源信息、地質(zhì)信息等,對分析軍事地質(zhì)和特殊目標(biāo)具有重要意義。近年來,微地貌的提取技術(shù)主要基于數(shù)字高程模型數(shù)據(jù)(Digital Elevation Model)實現(xiàn)。麻莉[2]基于DEM數(shù)據(jù)提出了一種針對洪泛平原濕地微地貌單元的數(shù)字分析方法。Baade等[3]應(yīng)用TerraSAR-X雷達(dá)數(shù)據(jù)監(jiān)測秘魯南部海岸沙漠地區(qū)微地貌變化,比較了利用DEM提取微地貌的局限性。沖洪積扇微地貌廣義上講包括沖積扇與洪積扇兩種。沖積扇沉積體具有構(gòu)造“ 敏感性”,其復(fù)雜的沉積機制是解釋盆地邊緣構(gòu)造活動的鑰匙,也是“源-聚-匯” 系統(tǒng)中的“源” 位置[4]。其前沿多為潛水溢出帶和淺埋地下水分布區(qū),對野戰(zhàn)給水保障具有重要意義。目前,沖洪積扇微地貌提取技術(shù)多基于遙感影像的光學(xué)特征,DEM多用于輔助篩選泥石流分布區(qū)的溝谷和縱坡分析[5],這些微地貌識別方法缺少地學(xué)因子的參考,不適應(yīng)復(fù)雜多變的微地貌邊界的識別。本文以青藏高原西部,北緯32°26′,東經(jīng)79°36′為研究區(qū)域,提出了一種基于CART分類器的沖洪積扇提取模型。該區(qū)域全年降雨量少,季節(jié)性強,氣候寒冷干燥[6]。寒冷、晝夜溫差大導(dǎo)致冰川發(fā)育廣泛,形成了高山、溝谷、冰蝕、冰磧、沖積扇等地貌類型,為沖洪積扇提取提供了大量的實驗對象。該模型以GF-2影像、Landsat 8影像、DEM影像為數(shù)據(jù)源,結(jié)合沖洪積扇微地貌的光學(xué)特征、幾何特征、紋理特征、地學(xué)特征位置特征,能夠準(zhǔn)確界定沖洪積扇扇緣,實現(xiàn)沖洪積扇高精度提取。
數(shù)據(jù)源包括Landsat8影像、GF-2衛(wèi)星影像、數(shù)字高程模型。其中,數(shù)字高程模型根據(jù)1∶5萬地形圖生成;Landsat8影像數(shù)據(jù)來源于中國科學(xué)院計算機網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn);GF-2影像為1 T級產(chǎn)品,參數(shù)如表1所示。
表1 GF-2衛(wèi)星數(shù)據(jù)波段參數(shù)
數(shù)據(jù)處理包括幾何校正、濾波降噪、圖像融合3部分。幾何校正主要針對GF-2遙感影像、DEM影像、Landsat 8影像進(jìn)行圖像配準(zhǔn),依靠重采樣對圖像亮度值進(jìn)行插值計算,建立新的圖像矩陣,結(jié)果誤差控制在0.5個像元以內(nèi)。濾波降噪通過中值濾波器實現(xiàn),濾波窗口大小為5×5。中值濾波能夠有效濾除椒鹽噪聲,盡量保護(hù)邊緣信息的完整[7]。圖像融合的主要目的在于借助GF-2全色段影像的高分辨率優(yōu)勢,提高洪積扇邊緣的識別精度。融合采用主成分分析法,融合效果如圖1所示。
沖洪積扇具有獨特的地學(xué)特征和影像特征,深入研究這些特征因子是構(gòu)建沖洪積扇提取模型的必要前提。
1) 地學(xué)特征分析
① 沖洪積扇頂位于河流出山口,地學(xué)特征表現(xiàn)為高程的快速下降,位置關(guān)系表現(xiàn)為鄰接山體與平原。
② 沖洪積扇體在遙感影像上的幾何形態(tài)通常呈現(xiàn)為近似扇形。其發(fā)育范圍取決于山區(qū)河流的匯水面積、巖體風(fēng)化特征、水文特征等。
③ 沖洪積扇體的物質(zhì)組成取決于其河流上游匯水范圍的巖體類型;不同的物質(zhì)成分在遙感影像上呈現(xiàn)不同的光譜特征。
④ 沖洪積扇體的土壤顆粒較粗,主要為砂、礫,通常植被稀疏,不利于植物生長。扇體邊緣和沖洪積扇前緣的土壤顆粒較細(xì),一般為沙、粉沙及亞黏土,通常植被發(fā)育,植被覆蓋率高。
⑤ 沖洪積扇在地形上呈現(xiàn)為向盆地傾斜的山前傾斜平原或傾斜臺地,其坡度為5°~10°,扇體前緣的坡度為1°~2°。
2) 影像特征分析
沖洪積扇體的反射率隨波長變化的趨勢與沖洪積扇土體顆粒變化基本一致。由于扇體前緣土體顆粒變細(xì),以沙、粉沙以及亞黏土為主,土體粒徑約在0.45~0.25 μm之間,土壤孔隙率相對增高,土壤含水率也明顯增高,其在各波段的反射率略低于沖洪積扇主體的反射率。沖洪積扇前緣地下水溢流帶產(chǎn)生的“微陰影”以及相對高的土壤含水率,導(dǎo)致扇體前緣部分的反射率、植被指數(shù)和地表溫度的突變,形成沖洪積扇前緣分界帶。
1) 植被特征指數(shù)
歸一化差分植被指數(shù)[8](NDVI)常用于監(jiān)測植被生長的狀態(tài)、植被覆蓋度,能夠有效排除植被密集區(qū),計算公式如下:
(1)
式(1)中:R表示影像的紅波段,NIR表示影像的近紅外波段。
2) 溫度植被干旱指數(shù)
溫度植被干旱指數(shù)[9](TVDI)可用于干旱檢測以及研究干旱程度的空間變化特征,計算公式如下:
(2)
式中:Ts表示給定像元的地表溫度,Tsmin、Tsmax表示給定像元相同NDVI值下的最低溫度和最高溫度。根據(jù)Landsat 8熱紅外波段影像反演研究區(qū)溫度植被干旱指數(shù),指數(shù)值越接近0,土壤濕度越大。反演結(jié)果如圖2所示。
3) 紋理特征指數(shù)
灰度共生矩陣是一種描述紋理特征的常用方法,是建立在估計圖像的二階組合條件概率密度函數(shù)基礎(chǔ)上的統(tǒng)計方法,能夠有效提取紋理特征[10]。通過灰度共生矩陣一定程度上能夠識別沖洪積扇頂?shù)暮拥兰y理。選用7×7窗口提取5種統(tǒng)計量,分別為一致性(Homogeneity)、對比度(Contrast)、相異性(Dissimilarity)、熵(Entropy)、均值(Mean)。其計算公式如下:
(3)
式(3)中:P(i,j)表示i行j列像元的亮度值,u為7×7窗口內(nèi)的像元亮度均值。
4) 坡度、坡向
坡度(Slope)、坡向(Aspect)能充分呈現(xiàn)高程的變化率以及變化方向,有效識別沖洪積扇頂?shù)牡貙W(xué)特征。
(4)
式(4)中:m表示點在x方向上高程變化率,n表示點在y方向上高程變化率。坡度、坡向計算一般是在最小像元的基礎(chǔ)上,以5×5窗口為單位進(jìn)行局部微分或曲面回歸擬合得出。
5) 光譜特征參數(shù)
沖洪積扇的物質(zhì)組成主要以河流攜帶的沙礫石和粘土為主,分選性較好,植被極少。其在多光譜影像上色調(diào)比較均勻、亮度值變化差異較小。對沖洪積扇及該區(qū)域其他地物進(jìn)行樣本采集分析,結(jié)果如表2所示。
表2 沖洪積扇及其他地物光譜值
決策樹是一種經(jīng)典的分類器,其擴(kuò)展通過節(jié)點的類別情況,每個類別導(dǎo)致一個結(jié)果[11]。分類器使用改進(jìn)于傳統(tǒng)決策樹的CART(Classification And Regression Tree)分類方法,即分類與回歸樹[12]。
CART為在給定輸入x的條件下,輸出隨機變量y的條件概率分布。CART算法二分每個特征(包括標(biāo)簽特征以及連續(xù)特征),通過最優(yōu)二分特征及其最優(yōu)二分特征值的選擇、切分,生成二叉樹并進(jìn)行剪枝[13]。分類過程中,回歸CART樹用誤差平方和準(zhǔn)則選擇特征,分類CART樹用基尼系數(shù)準(zhǔn)則選擇特征,并遞歸調(diào)用構(gòu)建二叉樹過程生成CART樹?;貧w樹實現(xiàn)流程如圖3所示。
1) 計算DEM影像的坡度值。以3°和10°為閾值,將坡度低于3°的區(qū)域視為沖洪積扇緣發(fā)育區(qū);將坡度為3°~10°的區(qū)域視為沖洪積扇體發(fā)育區(qū);將坡度高于10°的區(qū)域視為沖洪積扇頂發(fā)育區(qū)。
2) 計算目標(biāo)區(qū)域多光譜影像的紋理特征指數(shù)。沖洪積扇頂?shù)挠跋窦y理特征與山體存在明顯區(qū)別,條帶狀河流沖刷痕跡明顯。為縮減灰度共生矩陣計算量,基于Matlab軟件平臺對影像的灰度級數(shù)進(jìn)行壓縮,即從影像最小像元值起,每5級壓縮為1級,如1~5視為1,6~10視為2。
3) 計算NDVI、TDVI指數(shù)。研究區(qū)地處高原,多鹽堿地,植被指數(shù)主要與土壤類型相關(guān),參考表2的平原臺地地光譜特征指數(shù),以0作為區(qū)分閾值。根據(jù)沖洪積扇體表面植被覆蓋度極低的原則,將NDVI指數(shù)大于0的區(qū)域視為植被存在區(qū)域,排除在沖洪積扇可能存在的區(qū)域之外。從沖洪積扇根部體向扇體邊緣的土壤含水率逐漸增高,到扇前緣的地下水水溢出帶土壤含水率有明顯躍升,多泉水發(fā)育。根據(jù)TDVI指數(shù),以0.3作為閾值提取影像中的高表層土壤含水率的區(qū)域。
4) 根據(jù)影像光譜特征參數(shù),以藍(lán)波段像元值720為閾值,區(qū)分鹽堿地與其他地物,排除在沖洪積扇可能存在的區(qū)域之外。
5) 基于決策樹法排除有明顯特征的非沖洪積扇地物,界定山體與山間盆地的分界線?;贑ART樹獲取最優(yōu)分割特征值,實現(xiàn)不同源頭沖洪積扇的提取與區(qū)分。
6) 基于沖洪積扇與山體之間的鄰接關(guān)系判定沖洪積扇,實現(xiàn)沖洪積扇提取。
為檢驗?zāi)P托Ч?,分別使用本文提出的模型、最大似然法、SVM法對沖洪積扇進(jìn)行提取,并通過混淆矩陣對提取結(jié)果進(jìn)行精度評價,結(jié)果如圖4、表3所示。
注:圖中白色表示沖積扇A,灰色表示沖積扇B,黑色表示洪積扇C
表3 沖積扇分類結(jié)果
從表3看,本文提出的模型提取效果明顯優(yōu)于最大似然法和SVM法,對影像中3種不同源頭的扇體提取結(jié)果的Kappa系數(shù)分別為0.95、0.86、0.90,總體精度達(dá)到92.37%。結(jié)合圖4所示,模型通過沖洪積扇與山體之間的拓?fù)溧徑雨P(guān)系,有效濾除了 “同譜異物”的其他背景地物,極大地提高了提取精度。
1) 以Landsat8影像、GF-2衛(wèi)星影像、數(shù)字高程模型為數(shù)據(jù)源,基于CART分類器的沖洪積扇提取模型,實現(xiàn)了沖洪積扇微地貌單元的精確提取,為區(qū)域地貌環(huán)境的監(jiān)測提供技術(shù)支持。該模型提取沖積扇的總精度為92.37%。
2) DEM數(shù)據(jù)是沖洪積扇提取的重要參數(shù),結(jié)合光學(xué)、地學(xué)特征的模型提取精度遠(yuǎn)優(yōu)于僅基于光學(xué)特征的最大似然法和支持向量機方法。
3) 溫度植被干旱指數(shù)能有效表達(dá)土壤濕度,對識別沖洪積扇頂部、前緣具有重要作用,能夠準(zhǔn)確劃分沖洪積扇前緣地下水溢出帶。