張 艷 可，王 金 亮*，蘇 懷，程 峰

(1.云南師范大學(xué)地理學(xué)部，云南 昆明 650500；2.云南省高校資源與環(huán)境遙感重點實驗室，云南 昆明 650500；3.云南省地理空間信息技術(shù)工程技術(shù)研究中心，云南 昆明 650500)

0 引言

地貌是地球表面各個圈層系統(tǒng)間相互作用的重要一環(huán)[1]，對地貌進行分類研究不僅可以客觀了解當?shù)氐匦蔚膶嶋H情況，而且可為土地利用、水土保持研究等多個領(lǐng)域提供依據(jù)[2，3]。早期研究主要采用人工方法輔助野外實地勘察獲取地貌信息[4]。隨著GIS和RS技術(shù)的發(fā)展，國內(nèi)外學(xué)者開始研究地貌類型的自動或半自動劃分方法。目前地貌類型自動劃分方法主要包括：1)基于地貌單元的劃分方法。主要有基于柵格單元[5-8]和面向?qū)ο骩9-12]兩種統(tǒng)計方法，前者會破壞地理實體的完整性，后者雖可在一定程度上保證地表單元的完整性，但基于圖像分割技術(shù)劃分的圖像單元的地理意義不明確[13]。小流域作為自然形成的地理單元，在地表形態(tài)和地貌演變方面具有明確的地理意義，能夠很好地契合地貌和水文之間的關(guān)系[13]。但目前研究多應(yīng)用單一流域尺度[13-16]，當流域單元過小時，雖能保證不同地貌類型分界處界線清晰，但易造成相同地貌單元被過度分割，導(dǎo)致分類結(jié)果破碎；若流域單元過大，不同類型的地貌單元可能被歸并為同一地貌類型，即出現(xiàn)漏分現(xiàn)象，在地勢起伏較大的山地會更明顯。2)基于地形因子的劃分方法。有學(xué)者應(yīng)用絕對高程[17]、相對高程[18]或二者的組合[4,19]進行地貌劃分，基本可以篩選出某地區(qū)的地貌類型，但未考慮微觀地形單元；也有學(xué)者考慮到更多因子對地貌劃分的影響，如利用海拔、坡度、平面曲率、剖面曲率等多個地形因子劃分地貌類型[20],或者從地形因子統(tǒng)計量、特征點線統(tǒng)計量等方面選取地形因子[13，14]，使地貌類型包含的信息量更豐富，但參與分類的因子太多會造成信息冗余，影響分類結(jié)果。

綜上，在地貌劃分單元選取中，為避免劃分結(jié)果破碎或漏分現(xiàn)象，本研究擬采用雙尺度流域單元對地貌類型進行劃分：即先利用小流域單元提取出劃分效果明顯的不同地貌單元分界處，再采用較大流域單元對破碎地貌類型進行二次整合劃分，保證不同地貌類型的異質(zhì)性和相同地貌類型的同質(zhì)性。另外，以往地貌劃分方法通常只能處理單一或少量因子[1,8，21]，而CART決策樹作為典型的二叉樹算法，可以處理非數(shù)值型數(shù)據(jù)，具有很強的穩(wěn)健性和對多源數(shù)據(jù)的深度挖掘能力[22]，在基于遙感數(shù)據(jù)的分類研究中取得了較好的分類結(jié)果[23-25]，但目前鮮見利用該方法進行地貌自動分類研究。本研究為適應(yīng)多個地形因子參與分類，選用CART決策樹算法進行地貌劃分，選取北回歸線(云南段)地區(qū)為研究區(qū)域，基于ASTER GDEM 提取的多個地形因子進行最佳地形因子組合篩選，得出研究區(qū)主要地貌類型并進行檢驗，以期為研究區(qū)旅游開發(fā)、植被生長、氣象研究等提供數(shù)據(jù)支持。

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

北回歸線(23°26′21″N)又名夏至線，是北半球熱帶和亞熱帶的分界線。本研究以北回歸線在云南省穿過的縣市作為研究區(qū)(簡稱“北回歸線(云南段)”)(圖1)，包括臨滄、玉溪、普洱、紅河和文山5個州(市)的富寧、西疇、麻栗坡、硯山、文山、蒙自、個舊、石屏、紅河、建水、元江、墨江、寧洱、景谷、雙江、滄源、耿馬17個縣(市)，總面積達54 074.95 km2。研究區(qū)自西向東跨越滇西縱谷區(qū)(屬橫斷山脈南延帶)、滇中高原和滇東高原(云貴高原組成部分)，西部山河相間、縱列分布，無量山脈、哀牢山脈、瀾滄江水系、元江水系等分布在此；中部高原地區(qū)地勢起伏平緩，屬滇中大面積塊斷式上升區(qū)，分布有蒙自壩、平遠壩等；東部高原位于云南高原邊緣地帶，絕對海拔低、喀斯特地貌明顯。豐富的地貌特征為驗證本文方法在不同地貌分類中的適應(yīng)性提供了條件支持。

圖1 北回歸線(云南段)位置及高程分布

研究數(shù)據(jù)采用先進星載發(fā)射和反射輻射儀全球數(shù)字高程模型(ASTER GDEM)數(shù)據(jù)，由美國航天局(NASA)與日本經(jīng)濟產(chǎn)業(yè)省(METI)利用NASA新一代對地觀測衛(wèi)星TERRA的觀測數(shù)據(jù)制作而成(http://www.gscloud.cn/)，水平誤差為30 m，垂直誤差為20 m[26]。

2 研究方法

本文方法流程(圖2)為：1)基于DEM數(shù)據(jù)進行地形因子提取，利用相應(yīng)指標進行最佳地形因子組合計算；2)進行雙尺度流域單元劃分，以劃分的小流域為單元對各因子進行分區(qū)統(tǒng)計；3)經(jīng)流域單元統(tǒng)計后的各個最佳地形因子與以上因子的非監(jiān)督分類結(jié)果(CART決策樹分類的必需波段)合并，得到本文的待分類數(shù)據(jù)；4)將選取的訓(xùn)練樣本數(shù)據(jù)(地貌類別)和多源地形因子數(shù)據(jù)(地貌屬性)代入CART工具中進行自動運算，得到分類規(guī)則決策樹；5)利用該決策樹對待分類數(shù)據(jù)進行自動劃分，可得研究區(qū)地貌類型分布圖，并對分類結(jié)果進行檢驗。

2.1 最佳地形因子提取

綜合考慮地貌的宏觀和微觀特征，并結(jié)合前人經(jīng)驗[7,19,27]，選定高程(E)、地勢起伏度(RF)、地表切割深度(SCD)、高程變異系數(shù)(EVC)4個宏觀地形因子和坡度(S)、坡度變率(SOS)、地表曲率(SCV)3個微觀地形因子以及可間接反映地勢起伏度的光照模擬值(LSV)[19]，各因子計算公式見式(1)-式(4)。以上地形因子均由DEM數(shù)據(jù)衍生得到，會存在信息冗余，在分類前需通過因子標準差和相關(guān)性計算進行篩選[28](標準差可反映地形因子對于均值的離散程度，值越大則灰度級分布越分散，包含的信息量越大；相關(guān)性能篩選出攜帶重復(fù)信息的地形因子，減少信息冗余，提高分類精度)。提取宏觀地形因子地勢起伏度時，需確定鄰域統(tǒng)計單元的大小，本文采用均值變點法[29]計算得到其最佳分析窗口為14×14(1.76 km2)。

圖2 技術(shù)流程

(1)

式中：Estd、Emean、Emax、Emin分別為統(tǒng)計窗口內(nèi)的高程標準差、平均高程、最大高程和最小高程。

(2)

式中：H為垂直高差；L為水平距離；ΔS為坡度差。

SCV=CC×CP

(3)

式中：CC為平面曲率；CP為剖面曲率。

LSV=255×(cos(Z_r)×cos(S_r)+sin(Z_r)×sin(S_r)×cos(Azi_r-Asp_r))

(4)

式中：Z_r為太陽天頂角弧度，本文設(shè)定為45°；S_r為某點坡度的弧度；Azi_r為太陽光線方向角的弧度，本文設(shè)定為315°；Asp_r為某點坡向的弧度。

2.2 雙尺度流域單元劃分標準

最佳地形因子確定后，需利用流域單元對各因子進行分區(qū)統(tǒng)計。為避免地貌類型劃分結(jié)果破碎或出現(xiàn)漏分現(xiàn)象，本文采用兩種不同大小的流域單元(即雙尺度流域單元)作為地貌劃分基本單元，即先利用較小流域單元進行地貌類型初步劃分，再利用較大流域單元對首次劃分結(jié)果中的破碎地貌區(qū)域進行二次整合劃分。因此，兩次流域單元大小的確定是關(guān)鍵。結(jié)合前人經(jīng)驗[13-15]，本文首次劃分的較小流域單元采用定性方法確定，要求分割出的子流域能夠基本涵蓋山體的谷、脊線[15]。二次整合的較大流域單元采用定性和定量相結(jié)合方法確定，從定性角度講，要求其對破碎地表單元有明顯合并效果，且能兼顧地貌單元的完整性和內(nèi)部同質(zhì)性特征；從定量角度講，選取高程和地勢起伏度作為驗證因子，利用Moran′sI(式(5))驗證兩地形因子在不同流域單元的地貌分類結(jié)果，確定出二次整合的最佳流域單元。驗證標準為：兩個地形因子的Moran′sI均位于[0,0.3]區(qū)間，表明同一地貌類型中各單元間差異性過高，位于(0.3,0.6]區(qū)間，各單元間相關(guān)性最佳，位于(0.6,1]區(qū)間，各單元間同質(zhì)性過高。

(5)

式中：n表示研究空間的對象數(shù)目；wij為第i、j區(qū)域的空間權(quán)重數(shù)值，兩者相鄰為1，否則為0；xi和xj為對象i、j的屬性值，x為該屬性的平均值(即平均高程、平均地勢起伏度)。

2.3 CART決策樹構(gòu)建

CART決策樹[30]的核心思想是根據(jù)學(xué)習(xí)區(qū)域變量和目標變量，對原始數(shù)據(jù)集以二叉樹的結(jié)構(gòu)形式進行循環(huán)分析[24]，通過計算Gini系數(shù)選取屬性集中的某個屬性作為二叉樹的分類屬性，并利用該屬性將待分類樣本集分為兩個子集，重復(fù)執(zhí)行此步驟，直到當前分類的樣本集達到葉節(jié)點[25]。構(gòu)建決策樹需要合格的地貌類型訓(xùn)練樣本和待分類數(shù)據(jù)集：本文待分類數(shù)據(jù)為多源地形因子數(shù)據(jù)，它包括經(jīng)流域單元分區(qū)統(tǒng)計后的各個最佳地形因子及其非監(jiān)督分類結(jié)果，利用GIS軟件中的波段合成工具合并而得；下面進行地貌類型訓(xùn)練樣本選取和可分離度計算，為便于地貌樣本點選取，先對地貌類型進行初步確定。

2.3.1 地貌類型初步確定 目前關(guān)于中國地貌類型劃分尚無統(tǒng)一標準，本研究主要參考沈玉昌等[31，32]的劃分標準，并結(jié)合研究區(qū)的空間尺度和高原地形進行適當調(diào)整。綜合考慮最終地貌劃分結(jié)果的細致性、實用性及對劃分方法驗證的完備性，本文利用高程和地勢起伏度將研究區(qū)劃分為11種地貌類型(表1)。

表1 地貌類型劃分標準

2.3.2 樣本點選取 按照代表性、均勻性和一定數(shù)量原則對以上地貌類型進行訓(xùn)練樣本選取，樣本點選取標準為：各地貌類型樣本的選取同等重要，面積越大，選取的樣本點越多；在相似地貌類型(如中起伏中高山和大起伏中高山)樣本點選取中，兩者差異性明顯區(qū)域增加樣本點數(shù)量，相似性較高區(qū)域減少樣本點數(shù)量。最終選取的樣本點數(shù)量及面積如表2所示，其中小、中、大起伏低山3種地貌類型所占面積較小，可分離度太低，為不影響整體分類精度，將其合并為中起伏低山。

2.3.3 樣本可分離度計算 選取J-M距離(式(6))[33]和轉(zhuǎn)換離散度(式(7))確定不同地貌類型樣本的可分離度，二者分別表示類別的可分離性和離散情況，取值范圍均為[0，2]：當數(shù)值較小時，表明地物可分離性較差，需進行適當合并；數(shù)值越大，代表兩類地物可分離性和離散性越強，當數(shù)值大于1.8時，可進行分類計算[33]。

表2 各地貌類型樣本點數(shù)量及面積

(6)

式中：B為c1、c2兩種類別基于f特征的巴氏距離；m1、m2和σ1、σ2分別為c1、c2兩種類別基于f特征的均值和方差[33]。

(7)

式中：i、j為要比較的兩類；Ci為i類的協(xié)方差矩陣；μi為i類的均值向量；tr為求矩陣的跡。

2.4 分類結(jié)果檢驗

分類結(jié)果檢驗包括指標檢驗和人工判讀檢驗兩部分。指標檢驗中，除使用上述的Moran′sI外，還采用平均值和標準差評估結(jié)果的可靠性，前者用于檢測地貌類型劃分指標的準確性，后者可以檢測不同地貌的異質(zhì)性程度，海拔越高、地勢起伏度越大，標準差越大。

3 結(jié)果與分析

3.1 最佳地形因子選取結(jié)果

通過地形因子相關(guān)性計算結(jié)果(表3)得出，相關(guān)性最大的地形因子是地表切割深度(SCD)和地勢起伏度(RF)，后者已被多數(shù)研究[6,8,18]選定為地貌劃分的主要地形因子。由于不同因子的單位和值域不同，需先對其進行歸一化處理，而后計算出高程(E)、高程變異系數(shù)(EVC)、地勢起伏度(RF)、坡度(S)、坡度變率(SOS)、地表切割深度(SCD)、地表曲率(SCV)、光照模擬值(LSV)的歸一化標準差分別為0.125、0.023、0.099、0.135、0.115、0.099、0.040、0.388，其中EVC和SCV的標準差明顯低于其余地形因子，予以剔除。最終確定的最佳地形因子組合為：高程(E)、地勢起伏度(RF)、坡度(S)、坡度變率(SOS)、光照模擬值(LSV)。

表3 地形因子相關(guān)性系數(shù)

3.2 雙尺度流域單元地貌類型劃分結(jié)果

3.2.1 雙尺度流域單元劃分 由不同流量閾值下地貌劃分單元地形因子的Moran′sI(表4)可知，當流量閾值從500增至1 000時，除中起伏中山的高程和中起伏中高山的地勢起伏度外，其余地形因子的Moran′sI呈增加趨勢，表明相關(guān)性隨著流域單元的擴大而增大；閾值從1 000增至2 000時，除壩子和河谷平原外，其余地形因子的Moran′sI仍呈增加或不變趨勢，即相關(guān)性仍在提高；閾值從2 000增至3 000時，多數(shù)地形因子的Moran′sI出現(xiàn)下降(部分降至0.3以下)，少部分呈上升趨勢(大于0.6)，可見流量閾值為2 000時的流域單元劃分效果最好，此時Moran′sI多在(0.3,0.6]區(qū)間內(nèi)。綜上，結(jié)合定性和定量兩種標準，本研究確定流量閾值為500時劃分的流域單元為第一次地貌劃分統(tǒng)計單元，此單元足夠小且基本包含地貌單元的谷、脊線,流量閾值為2 000時劃分的流域單元為第二次地貌統(tǒng)計單元，歸并兩次分類結(jié)果得出北回歸線(云南段)地區(qū)的地貌類型劃分單元(圖3,彩圖見附錄4)。

表4 不同流量閾值下地貌劃分單元地形因子的Moran′s I

圖3 流域單元劃分結(jié)果

3.2.2 CART決策樹構(gòu)建及雙尺度流域單元地貌類型劃分

(1)由表5可知，各地貌類型樣本的可分離度均大于等于1.8，故可用該樣本構(gòu)建CART決策樹(圖4)。首先，基于流量閾值為500的小流域單元分區(qū)統(tǒng)計高程、地勢起伏度、坡度、坡度變率、光照模擬值及其非監(jiān)督分類結(jié)果(圖5，彩圖見附錄4)，合并為多源地形因子數(shù)據(jù);然后，運行決策樹模型對多源地形因子數(shù)據(jù)進行劃分(圖6，彩圖見附錄4)，可以看出，基于流量閾值500的流域單元分類結(jié)果基本可識別出研究區(qū)地貌類型特征，壩子、河谷平原、中起伏低山和部分小起伏中山的劃分效果較好，地貌單元集中連片，破碎化程度低，地貌類型間界線明顯，但剩余地貌類型明顯存在破碎化程度高、山地被過度分割情況，在滇中、滇東高原交界地帶和西部縱谷地區(qū)尤為明顯。

表5 不同地貌類型可分離度計算結(jié)果

(2)保留壩子、河谷平原、中起伏低山、部分小起伏中山等分類效果較聚集和地貌類型分界較明顯區(qū)域，進而對過于破碎的地貌單元區(qū)擴大統(tǒng)計單元進行二次劃分(圖7，彩圖見附錄4)。可以看出，第二次劃分結(jié)果中的破碎地貌類型區(qū)有明顯合并效果，且符合客觀實際，兩次劃分結(jié)果地貌類型分布位置基本一致。但受研究區(qū)本身地貌形態(tài)和算法分類誤差影響，小起伏中山和中起伏中山、小起伏中高山和中起伏中高山、中起伏中高山和大起伏中高山等相似地貌單元交錯分布，在西部縱谷地區(qū)尤為明顯。這也符合實際的地貌特征，即相鄰地貌類型之間存在交錯分布。

圖4 地貌分類規(guī)則決策樹

圖5 地形因子統(tǒng)計結(jié)果

圖6 地貌類型劃分結(jié)果(流量閾值為500的流域單元)

(3)由北回歸線(云南段)地區(qū)地貌類型劃分結(jié)果(圖7)可以看出，研究區(qū)地貌類型以山地為主，其中小起伏中高山占比最高，主要分布在滇中、東高原區(qū)，分屬滇中大面積塊斷式上升、滇東穹隆式上升區(qū)[34]，符合實際地貌特征；其次為中起伏低山，主要分布在東部的富寧縣和西部低谷地(如元江谷地、南汀河谷地等)兩側(cè)，兩地海拔低，地勢有一定起伏，地貌類型主要為中、小起伏低山，與劃分結(jié)果一致；中起伏中山占比居第三，主要分布在西部縱谷區(qū)和滇中、東高原過渡區(qū)(富寧縣西部、麻栗坡縣、西疇縣)，地勢變化較大，且海拔相對較高，緊鄰小起伏中山和中起伏低山，側(cè)面說明CART決策樹算法劃分結(jié)果較科學(xué)，相似地貌類型在空間上具有相鄰和過渡性質(zhì)；剩余山地地貌類型有大起伏中山、中起伏中高山和大起伏中高山，主要分布在西部縱谷區(qū)，該地是橫斷山脈帚形地帶的擴展部位，受新老構(gòu)造運動影響，地勢起伏最大[34]，地貌多樣且變化頻繁。由此可見，本文山地地貌類型劃分結(jié)果符合北回歸線(云南段)地區(qū)實際地貌特征。平原地貌(壩子和河谷平原)中，壩子主要分布在石屏、建水、個舊、蒙自、文山、硯山等中東部高原地區(qū)，較為明顯的蒙自壩、平遠壩、建水壩均被識別；河谷平原主要分布在西部河谷地區(qū)，受水流作用影響，呈狹長形縱列分布在山地中間。

圖7 北回歸線(云南段)地區(qū)地貌類型分布

3.3 分類結(jié)果評估

3.3.1 指標檢驗 進一步計算不同地貌類型高程和地勢起伏度的平均值、Moran′sI和標準差[9](表6)?？梢钥闯?，9種地貌類型的高程和地勢起伏度均值都在規(guī)定的劃分標準內(nèi)，表明對分類指標的劃分基本準確；9種地貌類型的Moran′sI均大于0，說明劃分結(jié)果在相同地貌類型內(nèi)具有空間同質(zhì)性，但中起伏中山、小起伏中高山的Moran′sI值較大，表明這兩類地貌單元可能存在一定過分割現(xiàn)象；隨著地勢起伏度和高程增大，對應(yīng)地貌的標準差逐漸增大，說明不同高程和地勢起伏度的地貌間具有異質(zhì)性特征。

表6 地貌類型劃分結(jié)果參數(shù)檢驗

3.3.2 人工判讀檢驗 從劃分的地貌單元中隨機抽取140個小流域，將人工判讀結(jié)果與CART決策樹分類結(jié)果相比較(表7、表8)，140個小流域有115個自動分類結(jié)果與人工判讀結(jié)果一致，總體精度達到82.1%，Kappa系數(shù)為0.793。 25個錯誤分類點中7個分布在中、東高原地區(qū)，18個分布在西部縱谷地區(qū)，可見該方法對不同地貌類型的判別精度略有差別，對破碎程度小、地勢起伏小、異質(zhì)性強(河谷平原、中起伏低山、小起伏中山、小起伏中高山)的地貌劃分效果較好，精度可達85%以上；對破碎程度高、異質(zhì)性差(中起伏中山、大起伏中山、中起伏中高山、大起伏中高山)的地貌劃分效果稍差，但精度基本也在70%以上。整體而言，本文方法是地貌類型劃分的一種可行方法。

4 結(jié)論與討論

基于ASTER GDEM 30 m空間分辨率數(shù)據(jù)，利用最佳地形因子組合、雙尺度流域單元劃分、CART決策樹算法對北回歸線(云南段)地區(qū)主要地貌類型進行劃分，并對結(jié)果進行檢驗，主要結(jié)論如下：1)研究區(qū)地貌劃分最佳地形因子組合為高程、地勢起伏度、坡度、坡度變率、光照模擬值；研究區(qū)主要地貌為山地(占比93.39%)，其中小起伏中高山占比最高，壩子和河谷平原僅占6.61%；研究區(qū)雙尺度流域單元劃分的最佳流量閾值分別為500、2 000。2)利用平均值、標準差、Moran′sI和人工判讀結(jié)果對分類結(jié)果進行檢驗，結(jié)果顯示：基于CART決策樹的雙尺度流域單元地貌類型劃分方法具有較高的自動化程度和精確度，在北回歸線(云南段)地區(qū)總體精度可達到82.1%，Kappa系數(shù)為0.793，能夠準確識別出地貌類型空間分布特征，是地貌類型劃分的一種可行方法。

表7 CART決策樹分類結(jié)果與人工判讀結(jié)果的混淆矩陣

表8 不同地貌類型分類精度

本文雙尺度流域單元綜合考慮地貌劃分的細分性和完整性，提高了兩者在地貌劃分中的兼容性；利用CART決策樹算法對多個地形因子數(shù)據(jù)進行深度挖掘，綜合考慮多個地形因子對地貌劃分的影響，減少分類節(jié)點選擇時主觀因素的影響[35]，在保證分類結(jié)果科學(xué)性的同時也具有操作簡便性。但研究中仍存在不足之處：1)CART決策樹的分類結(jié)果對訓(xùn)練樣本依賴性較強，樣本選取會直接影響分類結(jié)果的準確性，這也是監(jiān)督分類中難以避免的問題。2)在地貌變化復(fù)雜區(qū)域存在對少量相似地貌類型的誤判，除受地貌形態(tài)復(fù)雜影響外，也受流域單元劃分大小的影響，“中山”、“中高山”單元劃分可能存在過度分割。今后需深入研究如何為不同地貌類型尋找合適的地貌劃分單元，以及如何應(yīng)用更加完備的定量方法得到最佳地貌劃分單元。另外，本研究為雙尺度流域單元劃分，在其他地區(qū)也可根據(jù)實際情況進行多尺度流域單元劃分。