肖麗瓊，周汝良**，陳永剛，蔡志勇，黃 彥

(1. 西南林業(yè)大學(xué) 地理與生態(tài)旅游學(xué)院，云南 昆明 650224；2. 中航通飛研究院有限公司 中國特種飛行器研究所，廣東 珠海 519000)

地貌區(qū)劃及制圖是地球科學(xué)系統(tǒng)研究與應(yīng)用中一項(xiàng)基礎(chǔ)性工作，對國土空間規(guī)劃、資源開發(fā)利用、自然災(zāi)害管理、生態(tài)環(huán)境保護(hù)、生態(tài)恢復(fù)治理等都具有重要意義. 中國地貌學(xué)家開展了多種地貌區(qū)劃方法研究，同時編制了各種地貌區(qū)劃圖[1-2]，如省級地貌分區(qū)[3]及部分地區(qū)的地貌分區(qū)圖等[4]. 李炳元等[5]利用半定量與定性方法將中國劃分為兩級地貌區(qū)，程維明等[6]基于2013 年提出的地貌區(qū)劃結(jié)果構(gòu)建了全國五級地貌區(qū)劃方案，但各種區(qū)劃方案結(jié)果都有所差別. 其輸出結(jié)果不確定性主要受分類原則、分類指標(biāo)、分類方法以及數(shù)據(jù)分辨率的影響，由于缺乏定量化數(shù)據(jù)標(biāo)準(zhǔn)和模型化的計(jì)算機(jī)分類方法，多使用定性或半定量化劃分地貌，造成了各種分類方法無法相容和統(tǒng)一.

在地學(xué)分類指標(biāo)方面有指標(biāo)多樣、等級較少，且較多研究趨向于利用海拔高度和地勢起伏度進(jìn)行劃分[5]. 然而中國地貌復(fù)雜多樣，區(qū)域差異明顯，僅利用海拔高度和地勢起伏度數(shù)據(jù)進(jìn)行地貌分類不能完全反映地表形態(tài)的整體性，必須尋找地貌特征明顯的指標(biāo)量化其地貌形態(tài)，將其細(xì)分為更易于管理與應(yīng)用的組成部分[7]. 在分類方法方面，隨著數(shù)字高程模型數(shù)據(jù)[8]、ArcGIS 軟件技術(shù)和遙感技術(shù)的發(fā)展，為地貌自動劃分提供了更為廣泛的數(shù)據(jù)和技術(shù)支持[9-12]. 現(xiàn)有的基于聚類、規(guī)則知識及典型樣點(diǎn)等地貌自動劃分方法，都受限于語義定義模糊、多尺度分割困難、空間信息難以量化的影響[13].聚類是利用統(tǒng)計(jì)學(xué)原理將地貌屬性相似的地理柵格歸并為空間斑塊，再將其定義為地形的某個部位[14-16]. 事實(shí)上，在山地地貌區(qū)域，空間變異性強(qiáng)，聚類生成的斑塊大小不一且極其瑣碎，因而不能保證地貌的完整性[17]，也無法產(chǎn)出符合認(rèn)知概念和規(guī)則知識的地貌類型，難以進(jìn)行地學(xué)解釋[18]. 其根本原因是地貌區(qū)劃是定義在空間等模式變量上的函數(shù)，而不是像素上的測度. 現(xiàn)有研究中面向?qū)ο蟮挠跋穹椒朔藮鸥穹诸惖钠扑閳D斑和“椒鹽現(xiàn)象”[19]，而分割尺度的確定與面向?qū)ο蟮目臻g聚類是關(guān)鍵[20]. 其次，由于衛(wèi)星影像存在包含地表覆蓋物信息，難以有效融合地貌特征和尺度差異等問題而制約著地貌分類的進(jìn)一步研究[21-22]. 數(shù)字高程模型(Digital Elevation Models，DEM)的優(yōu)勢在于可“過濾”地表覆蓋物信息而只表達(dá)地形自身信息，而“紋理”特征又是區(qū)分不同地貌形態(tài)的重要依據(jù)[23-24]，成為宏觀尺度上地形特征量化與識別研究的突破點(diǎn). 在數(shù)據(jù)分辨率方面，數(shù)字高程模型(DEM)分辨率越高，地形模擬越成功，但地貌形態(tài)具有分形幾何特征，隨分辨率的增加，數(shù)據(jù)量呈幾何級數(shù)遞增.選擇符合應(yīng)用要求又能兼顧計(jì)算機(jī)容量與處理能力的分辨率數(shù)據(jù)，制約著DEM 應(yīng)用范圍及地貌區(qū)劃研究[25].

基于以上分析，本文提出了能夠表達(dá)地貌特征的新地學(xué)指標(biāo)，并結(jié)合傳統(tǒng)的地貌指標(biāo)構(gòu)建分類體系. 為克服一次聚類產(chǎn)生的瑣碎圖斑及山地地貌完整性被分裂的問題，引入了“核顆粒”和“紋理指標(biāo)”特征概念，并提出了二次瑣碎圖斑再聚合的方法，進(jìn)行中國全陸域范圍的地貌區(qū)劃研究，該方法對分異性較強(qiáng)的地理空間分類制圖具有參考應(yīng)用價值.

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

1.1 研究區(qū)概況 中國位于亞洲東部、太平洋西岸，地勢西高東低呈三級階梯分布. 第一級階梯內(nèi)有號稱“世界屋脊”的青藏高原，其間聳立著多條著名的高大山系，其中喜馬拉雅山脈是地球上最新隆起的年輕山系之一. 第二階梯包含若干高原和盆地，從青藏高原北緣、東緣到大興安嶺、太行山、巫山、雪峰山. 第三級階梯地勢低平主要以高原、盆地和平原為主[26]. 中國地勢起伏較大，分布著許多高大的山脈及被這些山脈所圍繞隔開的大型地貌單元，由西向東逐級下降[27]. 其中，山地和高原是構(gòu)成中國地貌基本輪廓的主體，山地、丘陵和高原約占中國總面積的2/3[28].

1.2 數(shù)據(jù)源及預(yù)處理 下載覆蓋中國區(qū)域的SRTM3 數(shù)據(jù)集，柵格分辨率為90 m（數(shù)據(jù)來源于https://earthdata.nasa.gov）. 運(yùn)用ArcGIS 軟件底層函數(shù)對數(shù)據(jù)文件拼接運(yùn)算，將1°×1°的分塊地圖拼接成中國的整體地圖. 由于接縫邊緣存在像素灰階不一致，構(gòu)成肉眼可見的異?？p，采用濾波進(jìn)行消除.對于某些超過正常海拔值的區(qū)域，用鄰域灰階進(jìn)行填充或條帶過濾處理.

像元降尺度處理：數(shù)字高程模型是描繪地形地貌的基礎(chǔ)數(shù)據(jù)，經(jīng)過提取地形信息和再演算，可實(shí)現(xiàn)各種地表形態(tài)的數(shù)字化模擬表達(dá). 隨著DEM 分辨率的提高，可分辨和提取流域的尺度不斷提高，同時數(shù)據(jù)量也呈幾何級數(shù)增長. 本研究以全國陸域范圍為對象，考慮到研究區(qū)及大尺度地貌劃分對數(shù)據(jù)的需求，分別利用90 m、500 m、1 000 m 分辨率的DEM 數(shù)據(jù)提取地形指標(biāo)，經(jīng)測算比對分析1 000 m尺度柵格數(shù)據(jù)能滿足研究應(yīng)用要求，地貌碎塊的表達(dá)也較好，故以其為地理指標(biāo)的柵格計(jì)算尺度.

2 原理與方法

2.1 地貌新變量 地貌要素特征的量化表達(dá)一直是難點(diǎn)，傳統(tǒng)的海拔、坡度、坡向等不足以表達(dá)宏觀地貌特征，要實(shí)現(xiàn)地貌類的劃分需尋找新的特征指標(biāo).

2.1.1 地形隆起高度與切割深度 對于發(fā)育于大傾斜面上的山地地貌，在水平和垂直方向上具有復(fù)雜的向上隆起和向下切割的變化特征，如圖1 所示. 傳統(tǒng)的海拔難以表達(dá)地貌隆起、溝谷切割或盆地塌陷的狀態(tài). 基于此，提出以每個流域的平均海拔為樣本點(diǎn)構(gòu)建地表趨勢面，以此為起算點(diǎn)測算地形隆起高度與切割深度的變化特征[29]. 該指標(biāo)以山體基部為起算點(diǎn)，度量山體隆起或切割下沉程度.

圖1 沿24°N 云南地形的隆起與切割剖面Fig. 1 Uplift and cutting profile along 24°N mountain in Yunnan

2.1.2 皺褶多度 褶皺是指巖石因受力發(fā)生一系列的波狀彎曲而顯示的變形，使原來近于平直的面表現(xiàn)出向上隆起或塌陷的曲面形態(tài). 皺褶多度與地貌類型密切相關(guān)，平原地貌幾乎無皺褶，山地地貌破碎，皺褶明顯. 王艷霞等[30]基于滑動窗口方法，利用地表表面積與地物投影的平面面積比值提取褶皺多度（公式（1）），但該種方法采用面積線性對比反映研究區(qū)的褶皺變化程度，而中國山體坡度多介于20°～50°中等大小的山體發(fā)育，坡度差異小，故投影面積區(qū)分也不顯著，使得褶皺量化對比不明顯. 基于此，本文以具體坡度數(shù)值對褶皺多度指標(biāo)進(jìn)行改進(jìn)（公式（2）），以DEM 為數(shù)據(jù)源，利用ArcGIS軟件構(gòu)造褶皺地貌單元模型.

其中，R為褶皺多度，為一個常數(shù)；S為不同地貌單元內(nèi)的地形坡度（單位：（°））.

2.1.3 地表破碎度 地表破碎度描述了單位面積上的山體數(shù)量或密度，如中國南部的云貴高原和南嶺地區(qū)山體多而密，地形破碎明顯. 地表破碎程度需要從水平和垂直方向進(jìn)行刻畫，水平方向上揭示山體水平密度變化，可用地圖投影面上的河流切割密度表示[31]. 垂直方向上，選取山麓到山頂海拔的起伏變化表征垂直破碎形態(tài)，同時用區(qū)域內(nèi)地形的相對高差、地形隆起高度與切割深度及地形差異性表征該地學(xué)含義. 現(xiàn)實(shí)中垂直尺度上的地形隆起高度與切割也是垂直破碎度的主要影響因素，故將水平破碎及垂直破碎聯(lián)合起來，更能體現(xiàn)地貌整體的破碎程度. 考慮到不同地形指標(biāo)表現(xiàn)特征各異，取值范圍存在差異，為了權(quán)衡不同指標(biāo)的影響，對其進(jìn)行歸一化處理. 計(jì)算方法如下：

其中，A為地表破碎度，為一個常數(shù)；T為歸一化后的地形隆起高度與切割深度；St表達(dá)歸一化后的窗口的海拔變差，即計(jì)算窗口標(biāo)準(zhǔn)差；D為歸一化后的地表切割密度，即溝壑密度.

2.2 地貌類的柵格定義與二次聚合 利用

ArcGIS 軟件提供的ISO 聚類方法進(jìn)行多指標(biāo)統(tǒng)計(jì)聚類，將地貌劃分為碎片單元. 因各碎塊面積大小不一，有大面積的片區(qū)，也有小面積的椒鹽，通過消除將小面積逐步歸并到大片區(qū)中，被歸并的小面積碎片稱之為“核顆?！? 如果碎片面積大小相近，多類碎片可利用重復(fù)模式結(jié)構(gòu)表達(dá)地貌類，其集合體即為“紋理指標(biāo)”. 依據(jù)“核顆?！睂⒖臻g相鄰且相互關(guān)系密切的小碎片聯(lián)通為大片區(qū)“紋理”，可劃分出地貌類別.

2.2.1 地貌類定義和空間再聚合算法 假設(shè)集合Ω={A1，A2···，An}，其中元素Ai(i=1，2，···，n)為地貌類別的地學(xué)指標(biāo)，Ai由柵格像元組成并用數(shù)字矩陣描述. 用Ω中的元素進(jìn)行多元聚類或分類，可將Ω劃分為k類，表示為Ω1={1，2，···，k}，每一類代表距離度量相近的碎片集合.

依據(jù)統(tǒng)計(jì)學(xué)原理，類序號標(biāo)識相鄰的碎塊更為相近，即j+1 與j、j+2 相近，與j-1、j+3 相離(j=1，2，···，k)，按序號相近原則為優(yōu)先秩序，碎片再聚合聯(lián)通為大片區(qū).

2.2.2 “核顆?！钡拿枋龊吞幚?假設(shè)Ω1中k成為地貌類最小面積Sk，“核顆?！倍x為第k類中面積小于Sk的所有碎塊. 遍歷k中每個面積小于Sk的元素，并搜索其相鄰或有包含關(guān)系的碎塊，按類序號相近者優(yōu)先聚合，逐步將其歸并為片區(qū). 片區(qū)將繼承面積大的“核顆?！钡念悓傩?，面積較小者做為椒鹽被歸并而忽略其屬性. 通過不斷搜索和迭代，當(dāng)某片區(qū)的面積大于Sk時，該片區(qū)形成地貌單元，并停止再聚合運(yùn)算.

2.2.3 “紋理指標(biāo)”的描述和處理 若存在一種山地地貌類，由多種不同類編碼但面積近似的碎片單元組成，即該地貌類不存在“核顆?！焙徒符}可歸并關(guān)系時，需利用斑塊大小、斑塊之間的鄰近及可聯(lián)通性等拓?fù)潢P(guān)系，將碎塊圖斑迭代聚合為片區(qū)稱為“紋理指標(biāo)”.

假設(shè)Ω1的子集{j1，j2，···，jm}(j1，j2，···，jm分別取值于1，2，···，k)形成“紋理指標(biāo)”，給定“核顆?！钡拿娣e上下限Smax和Smin為約束條件，對滿足Smin

2.2.4 聚合迭代運(yùn)算 以Ω1={1，2，···，k；k+1，k+2，···}為新的初始地貌類，不斷增大Sk、Smax和Smin的面積約束值，重復(fù)進(jìn)行片區(qū)與椒鹽形成“核顆?！薄ⅰ凹y理指標(biāo)”的聚合迭代運(yùn)算，即可實(shí)現(xiàn)地貌對象的區(qū)域劃分.

3 結(jié)果與分析

3.1 地貌分類指標(biāo)分析 基于數(shù)字地貌信息提取的地形因子是一個定量化指標(biāo)，能客觀地反映地貌形態(tài)特征與發(fā)生發(fā)育演化的規(guī)律，揭示各指標(biāo)與地形地貌的因果聯(lián)系. 因此，地貌指標(biāo)的選擇成為地貌分類的關(guān)鍵. 不同的地形指標(biāo)從不同維度刻畫地貌特征，其間存在相關(guān)關(guān)系，為消除指標(biāo)相關(guān)性及信息冗余的影響，借助ArcGIS 軟件柵格代數(shù)運(yùn)算工具，對9 個指標(biāo)進(jìn)行中國全陸域范圍的相關(guān)性統(tǒng)計(jì)分析，表1 得到各指標(biāo)之間的相關(guān)系數(shù)矩陣. 經(jīng)對各指標(biāo)進(jìn)行可視化及相關(guān)性分析，剔除表1 中相關(guān)性在0.6 以上的指標(biāo)，但在進(jìn)行單個地貌指標(biāo)分析時發(fā)現(xiàn)，地表破碎度、地形起伏度雖相關(guān)性較高，但兩者從不同方面刻畫地貌形態(tài)特征，故都保留運(yùn)用于最佳地貌區(qū)劃指標(biāo)體系的篩選構(gòu)建中. 經(jīng)對各個指標(biāo)組合對比分析，本文最終選擇地形隆起高度與切割深度、皺褶多度、地表破碎度、海拔高度、地形起伏度共5 個指標(biāo)研究樣區(qū)地貌區(qū)劃.

表1 地貌指標(biāo)相關(guān)系數(shù)Tab. 1 Correlation among the geomorphology factors

3.1.1 地形隆起高度與切割深度 地表形態(tài)指標(biāo)是地貌區(qū)劃的主要依據(jù)之一. 地形隆起高度與切割深度能夠度量山體向上隆起、溝谷向下切割的程度，體現(xiàn)地貌在高度上的真正垂直變化，刻畫大傾斜面上大地域間構(gòu)造狀況及整體變化趨勢. 由圖2(a)可知，地貌隆起最大高度為2 755 m，切割最大深度為-2 051 m. 從青藏高原到周邊區(qū)域過渡帶，隆起數(shù)值較大，其中以橫斷山區(qū)隆起最為劇烈、隆起面積最大. 其次，青藏高原北面和柴達(dá)木盆地西面、南面的昆侖山，北面的天山地形變化呈窄條帶發(fā)育，海拔下降大、溝谷深切，隆起與向下切割形態(tài)顯著. 塔里木盆地、準(zhǔn)噶爾盆地、四川盆地、東北平原、華北平原和長江中下游平原等地為地貌切割平緩的典型代表單元，隆起數(shù)值較小.

3.1.2 皺褶多度 地貌形態(tài)特征主要受內(nèi)外營力作用影響反映研究區(qū)地貌差異性，如花崗巖長期受外營力作用侵蝕，在中國北方地區(qū)常呈現(xiàn)高大險峻的山地地貌如華山、泰山等. 研究定量化的物質(zhì)組成分異指標(biāo)對地貌區(qū)劃意義重大，褶皺多度能夠衡量特定區(qū)域內(nèi)地表整體形態(tài)，度量地貌坡面組合差異. 由圖2(b)可知，中國地貌單元內(nèi)褶皺度最低值集中分布于平原和中低山. 褶皺多度較明顯區(qū)域地形起伏劇烈、坡度大，主要分布于昆侖山脈、天山山脈、橫斷山區(qū)、祁連山區(qū)、秦嶺、臺灣山脈等高山區(qū).

3.1.3 地表破碎度 地表破碎度能從水平垂直方向上揭示研究區(qū)地表組合完整程度，刻畫空間上的分布規(guī)律. 從外營力方面表征地貌分帶現(xiàn)象，成為地貌區(qū)劃單元轉(zhuǎn)折部位劃分的依據(jù). 由圖2(c)可知，中國地表破碎度呈現(xiàn)較高空間異質(zhì)性，如東部的東北平原、華北平原、長江中下游平原等大部分地區(qū)；中西部的準(zhǔn)噶爾盆地、塔里木盆地、柴達(dá)木盆地、四川盆地等典型盆地地貌區(qū)，地表完整度良好，破碎程度小；中部的秦嶺、太行山、云貴高原等地地表相對崎嶇，存在一定程度的破碎；西部的天山、阿爾金山、昆侖山、祁連山、橫斷山等區(qū)域多高山峽谷分布，高差變化大，地表破碎顯著.

3.1.4 海拔高度 海拔高度是度量地表形態(tài)特征的最基本定量化指標(biāo). 以海平面為參照，用于測算地面上某一點(diǎn)與大地水準(zhǔn)面之間的高差，是地貌區(qū)劃的重要依據(jù). 由圖2(d)可知，中國地勢西高東低呈階梯狀分布，地勢按海拔大約可劃分為4 個梯次，低海拔和中海拔地貌占主導(dǎo). 其中，低海拔主要分布在第三級階梯上即東部平原及丘陵地區(qū)，以及準(zhǔn)噶爾盆地、四川盆地和塔里木盆地部分地區(qū)；中海拔主要分布于中國的第二級階梯上的云貴高原、內(nèi)蒙古高原、黃土高原、塔里木盆地等地區(qū)；青藏高原地區(qū)、橫斷山區(qū)、天山山脈等為中國高海拔地區(qū)；喜馬拉雅及昆侖山山脈主體部分為中國極高海拔區(qū).

圖2 主要地學(xué)指標(biāo)特征Fig. 2 Characteristics of main geoscience indexe

3.1.5 地形起伏度 地表起伏高度從地表形態(tài)方面研究地貌類型，能體現(xiàn)增強(qiáng)或均勻化地貌的復(fù)雜程度，區(qū)分區(qū)域尺度地形地貌的變化. 本文根據(jù)徐漢明等[32]研究確定中國最佳統(tǒng)計(jì)單元為22 km2進(jìn)行地形起伏度計(jì)算. 由圖2(e)可知，中國地形起伏度在空間上東西、南北差異明顯，西部多為高原及高大山區(qū)起伏大，東部地貌類型多為平原及丘陵起伏微小，南部主要表現(xiàn)為中等起伏，北部大部分地區(qū)地貌起伏和緩. 在利用窗口分析法提取指標(biāo)的過程中，因不同區(qū)域具有不同的最佳分析窗口，在進(jìn)行大尺度地形研究時難以兼顧不同大小的地貌實(shí)體單元，但該指標(biāo)能夠更好地表達(dá)窗口地形起伏變化，也能夠更清晰地表達(dá)山體地貌特征，成為地貌區(qū)劃研究最重要的指標(biāo)之一.

3.2 地貌類型劃分 用ArcGIS 軟件系統(tǒng)中ISO聚類和最大似然法劃分地貌其結(jié)果表現(xiàn)為“核顆?！苯徊娴幕ò呓符}結(jié)構(gòu)，難以形成真正的地貌類(圖3).其次，聚類過程中類別數(shù)目不同結(jié)果也存在差距，因此本文分別聚類為6、9、12、18、36、42 類進(jìn)行實(shí)驗(yàn)分析. 如當(dāng)聚類為6 類時，在塔里木盆地、四川盆地、柴達(dá)木盆地、準(zhǔn)噶爾等典型盆地地貌區(qū)存在較多小面積的椒鹽；在內(nèi)蒙古高原、黃土高原、云貴高原等地的高原地貌區(qū)由面積近似的多類碎塊組成，地貌分類結(jié)果較為破碎. 當(dāng)聚類為9 類時，因類別設(shè)置差異較小，故與聚類為6 類時結(jié)果較為相似，但典型地貌單元特征基本顯現(xiàn)，如四川盆地區(qū)基本結(jié)構(gòu)框架雖已初見雛形，但還摻雜較多小面積的椒鹽. 當(dāng)聚類為18 類時，太行山、秦嶺、大別山、昆侖山與祁連山脈地區(qū)多面積近似的碎塊. 當(dāng)聚類為36 類時（圖3），各個具有代表性的地形單元已經(jīng)具有典型的地貌特征. 當(dāng)聚類為12 類時，其結(jié)果與聚類為18 時差異較?。划?dāng)聚類為42 類時，破碎圖斑較多，地貌類型破碎. 因此，經(jīng)試驗(yàn)初步聚類為36 類最接近于真實(shí)地貌類型.

圖3 初始聚類地貌Fig. 3 Initial cluster geomorphic

初始聚類結(jié)果在平原和盆地地區(qū)表現(xiàn)最為突出，而在復(fù)雜山地地貌區(qū)表現(xiàn)較差，在此過程中地形屬性的有用信息與干擾信息同時增長. 如符合地理認(rèn)知的橫斷山區(qū)應(yīng)為獨(dú)立類，但聚類結(jié)果中具有多個顆粒類別，使得山地等復(fù)雜地貌區(qū)精確化難度大，無法產(chǎn)出符合認(rèn)知概念和具有現(xiàn)實(shí)意義的地貌類型. 因此，需進(jìn)行二次聚合處理實(shí)現(xiàn)中國地貌的計(jì)算機(jī)自動區(qū)劃與制圖.

3.2.1 “核顆?！钡纳?“核顆?！笔堑乩碜兞肯嘟南裨酆象w，在中國初始聚類中Sk定義為地貌類別中面積取值小于50 km2的碎塊. 按類序號相近者在中國范圍內(nèi)搜索相鄰或有包含關(guān)系的碎塊，將所有面積小于50 km2的小碎塊椒鹽進(jìn)行消除形成片區(qū)，但該片區(qū)還未形成獨(dú)立地貌“紋理”，故又將面積取值小于100 km2的椒鹽“核顆?！边M(jìn)行消除. 經(jīng)多次試驗(yàn)以倍數(shù)遞增聚合成片區(qū)，對中國地貌區(qū)劃研究具有普適性. 結(jié)合碎塊圖斑迭代聚合成“紋理”類地貌的方法進(jìn)行消除，當(dāng)最終聯(lián)通面積Sk為9×103km2時，二級地貌區(qū)初見雛形，即該區(qū)域完成地貌類的“核顆粒”聚合，形成新的地貌區(qū).

3.2.2 “紋理”類地貌的生成 華北、華東、塔里木盆地等平坦地形區(qū)，無“核顆粒”和椒鹽可歸并元素，即初始聚類就形成大面積片區(qū)“紋理” (圖3).除此之外，其它地貌類幾乎無法用聚類分析的輸出類確認(rèn)，需設(shè)置規(guī)則逐步迭代形成地貌類“紋理”.山地地貌區(qū)由多類碎塊組成且各碎塊的面積近似，不存在明顯的“核顆?！焙徒符}，優(yōu)先歸并關(guān)系無法聚合迭代運(yùn)算，需設(shè)置單獨(dú)算法. 圖4(a)顯示，喜馬拉雅高山、極高山，類別8、9 碎塊面積近似，若以特定“核顆粒”約束面積聯(lián)通為片區(qū)，西北部將被歸并為類別9，東南部歸并為類別8，而類別9 與其北部羌塘高原為一個地貌類別，會造成喜馬拉雅山區(qū)山體不完整. 若想將其進(jìn)行定量化單獨(dú)分區(qū)而不以人工勾繪方式，需設(shè)置單獨(dú)算法，即在保證其它區(qū)域地貌類別代碼不變的前提下. 按類序號相近者歸并形成片區(qū). 新形成的“紋理”看成“核顆?！崩^續(xù)合并，如此循環(huán)，直至該區(qū)域界線明顯后將其融合到中國地貌矢量圖中. 圖4(b)顯示，優(yōu)先保證類別8連成片區(qū)，并將該圖層融合到中國地貌圖內(nèi)形成新“紋理”進(jìn)行運(yùn)算. 當(dāng)逐步迭代到300 km2時，可通過人機(jī)交互識別設(shè)置“核顆粒”面積的上下限為350 km2和300 km2，進(jìn)行多類多斑塊同時聯(lián)通，形成如阿爾泰高山區(qū)等新地貌類，如此迭代實(shí)現(xiàn)地貌區(qū)域劃分.

圖4 初始聚類地貌和紋理類地貌的處理與生成Fig. 4 Processing and generation of initial clustering landform and texture classes

3.2.3 自動區(qū)劃及制圖 依據(jù)實(shí)際需要重復(fù)進(jìn)行碎片與椒鹽及“紋理”聚合的迭代運(yùn)算，完成中國地貌空間劃分與制圖. 當(dāng)聯(lián)通面積Sk為1.2×105km2時，一級地貌區(qū)基本框架形成，圖5(a)中將中國地貌劃分為6 個一級地貌區(qū). 三級地勢階梯特征是中國一級地貌區(qū)劃分的基礎(chǔ)，第一級階梯內(nèi)大致以3 500 m海拔高度所圈定的青藏高原為1 個地貌區(qū)；第二級階梯內(nèi)將其劃分為3 個地貌區(qū)，秦嶺南北地貌組合明顯不同，以南地區(qū)以中山、亞高山廣布，以北地區(qū)為海拔較低的盆地及高平原相間分布，而賀蘭山以西地區(qū)大型盆地與亞高山占主導(dǎo)，與賀蘭山以東的地貌區(qū)相比地形特征差異顯著，故單獨(dú)劃分為1個大區(qū). 第三級地勢階梯南北以大別山為界，東西以武夷山、雪峰山及苗嶺為界將其劃分為2 個地貌區(qū). 在一級地貌區(qū)骨架范圍制約下，再次進(jìn)行二次聚合迭代，依據(jù)地形自然紋理和已有的地貌區(qū)劃研究，聯(lián)通面積Sk為3×104km2時，大部分二級地貌區(qū)形成，并將其劃分為33 個二級地貌區(qū)(圖5(b))，能夠突出分區(qū)區(qū)域性. 東部低山平原大區(qū)（I）框架制約下劃分為7 個二級地貌區(qū)；東南低中山大區(qū)（Ⅱ）下劃分5 個二級地貌區(qū)；中北部中山高原大區(qū)（Ⅲ）范圍內(nèi)分5 個二級地貌區(qū)；西北部中高山盆地大區(qū)(Ⅳ)劃分為5 個二級地貌區(qū)；西南部亞高山、中山大區(qū)(Ⅴ)劃分為5 個二級地貌區(qū)；青藏高原大區(qū)(Ⅵ)制約下劃分為6 個二級地貌區(qū). 將本文結(jié)果與程維明等[5-6]的研究成果比對分析發(fā)現(xiàn)，一級地貌區(qū)邊界與之基本吻合，二級地貌區(qū)輪廓源于等高線的自然形態(tài).

圖5 中國地貌兩級分區(qū)的區(qū)劃結(jié)果Fig. 5 Results of two-grade regionalization of geomorphology in China

4 討論與結(jié)論

4.1 討論

（1）山地地貌類是柵格地理單元組成的“核顆粒”和“紋理指標(biāo)”的空間再聚合集合體. 基于統(tǒng)計(jì)聚類結(jié)果為基礎(chǔ)，利用斑塊大小、紋理特征、斑塊之間的鄰近及可聯(lián)通性等拓?fù)潢P(guān)系，將碎塊圖斑迭代聚合為片區(qū)可實(shí)現(xiàn)地貌類型的計(jì)算機(jī)自動區(qū)劃及制圖. 同時，自動區(qū)劃和制圖受綜合自然地理要素及地貌分類指標(biāo)體系的影響，需考慮不同地貌區(qū)邊界. 在進(jìn)行山地地貌區(qū)劃時，應(yīng)注重整個區(qū)域垂直地形特征并結(jié)合溝谷發(fā)育和大地構(gòu)造的差別，如秦嶺北坡山地邊界的劃分，以分水嶺為界從東向西根據(jù)垂直地形差異關(guān)系劃分地貌區(qū). 在以低山為主，且有個別中山分布區(qū)的界線確定時，可根據(jù)地形起伏狀況把突出的中山或與周邊區(qū)別較大區(qū)域單獨(dú)劃分出來. 在以平原為主，有低山及丘陵地貌區(qū)，可根據(jù)實(shí)際地貌制圖的平面結(jié)構(gòu)及各指標(biāo)的“紋理”特征進(jìn)行劃分，如魯東低山丘陵周邊各指標(biāo)在該區(qū)域表現(xiàn)平緩，低山與周圍其余地貌相比較突出，結(jié)合其所處地理位置故將其劃分為魯東低山丘陵區(qū).

（2）由于地貌指標(biāo)選取與分類方法不同，使其部分區(qū)劃邊界和屬性判斷存在差異. 而又缺少可對比的矢量數(shù)據(jù)，只能通過可視分析進(jìn)行比對，其結(jié)果與現(xiàn)實(shí)會存在些許偏差. 如在類別上，程維明等[5-6]提出江河源丘狀山原、羌塘高原湖盆區(qū)2 個二級地貌區(qū)，但因該區(qū)受地形影響各指標(biāo)表現(xiàn)出高度相似性，聚類后可形成大面積片區(qū)，因此將其合并劃分為羌塘高原區(qū)；橫斷山和江河上游高山峽谷區(qū)多高大山體發(fā)育，其形態(tài)相似聚合后由多類面積相似的碎塊組成，不存在優(yōu)先歸并順序故將其合并劃分；喀喇昆侖山極高山與昆侖極高山各指標(biāo)表現(xiàn)特征相似而合并為一個地貌類；寧鎮(zhèn)山脈（江蘇）面積較小，地形起伏變化和緩且與周圍地貌區(qū)分不顯著，故將其歸并到華北平原地貌區(qū). 在分界線上，程維明等[5-6]以伏牛山中部為界劃分秦嶺-大巴亞高山區(qū)，本文依據(jù)定量化的指標(biāo)表達(dá)為保證山體完整性仍以伏牛山為界. 除此之外，中國南部地區(qū)地形破碎，聚類后多面積近似多類碎塊鑲嵌形成的模式結(jié)構(gòu)，設(shè)置算法時具有不確定性，故在一級、二級地貌區(qū)及界線劃分上與已有研究結(jié)果存在差異. 但與中國地形數(shù)據(jù)對比分析發(fā)現(xiàn)，本文的地貌分類結(jié)果與原始數(shù)據(jù)及地貌的科學(xué)認(rèn)知體系相吻合，具有一定的現(xiàn)實(shí)意義.

4.2 結(jié)論

（1）本文構(gòu)建了地形隆起高度與切割深度、褶皺多度、地表破碎度、地形起伏度、海拔高度組合的中國地貌計(jì)算機(jī)自動劃分的指標(biāo)體系. 在新地貌指標(biāo)特征表達(dá)方面，地形隆起高度與切割深度指標(biāo)更能表征地貌的垂直變化，實(shí)現(xiàn)宏觀大地域間的地貌構(gòu)造差異；地表破碎度指標(biāo)用水平溝壑密度、垂直方向的隆起與切割高度狀況以及窗口地形變化特征共同描述地貌破碎程度，比以往的指標(biāo)更為全面、客觀；褶皺多度指標(biāo)從外營力方向表達(dá)了地貌起伏變化特征，對多種地理科學(xué)現(xiàn)象的研究都可提供借鑒作用.

（2）以聚類分析輸出的地貌類碎片斑塊為基本單元，構(gòu)建基于“核顆?！焙汀凹y理指標(biāo)”逐步迭代的二次空間聚合算法，可實(shí)現(xiàn)中國地貌的計(jì)算機(jī)自動區(qū)劃和制圖研究. 最終，發(fā)現(xiàn)聯(lián)通面積Sk為9×103km2時，二級地貌區(qū)開始形成；聯(lián)通面積Sk為3×104km2時，大部分二級地貌區(qū)形成；聯(lián)通面積Sk為1.2×105km2時，一級地貌區(qū)形成. 據(jù)此，將中國地貌劃分為6 個一級地貌區(qū)和33 個二級地貌區(qū). 二次聚合解決了山地地貌的計(jì)算機(jī)分類問題，該方法也適合于地理分異性較強(qiáng)的土壤制圖、植被分類制圖及遙感圖像分類研究等.