張 學(xué) 儒 ，官 冬 杰，牟 鳳 云，陳 春

基于ASTER GDEM數(shù)據(jù)的青藏高原東部山區(qū)地形起伏度分析

張 學(xué) 儒 ，官 冬 杰，牟 鳳 云，陳 春

（重慶交通大學(xué)河海學(xué)院，重慶 400074）

以青藏高原東部山區(qū)為研究區(qū)，基于高空間分辨率的ASTER GDEM數(shù)據(jù)，通過AML語言編程調(diào)用ArcGIS中用于鄰域分析的focal函數(shù)，計(jì)算不同鄰域尺度單元下地形起伏度。研究表明：鄰域尺度單元大小對(duì)地形起伏度計(jì)算至關(guān)重要，起伏度值先隨鄰域尺度單元面積增大而快速增大，當(dāng)鄰域尺度單元面積達(dá)一定閾值后，其增大速度開始減緩并趨于平穩(wěn)，且在增速減慢過程中存在一明顯拐點(diǎn)，即最佳鄰域尺度單元。通過高差顯著性變化檢驗(yàn)法，確定最佳鄰域尺度單元為5.0625 km2，據(jù)此制作地形起伏度分級(jí)圖，發(fā)現(xiàn)研究區(qū)自西北向東南地形起伏度逐步增加，地勢(shì)以中度起伏（200～500 m）為主。

地形起伏度;ASTER GDEM;鄰域分析;青藏高原東部山區(qū)

0 引言

地形起伏度作為反映地表起伏程度的一個(gè)指標(biāo)，在地貌制圖、資源環(huán)境評(píng)價(jià)等領(lǐng)域有廣泛應(yīng)用，是地貌類型提取和劃分的重要指標(biāo)，也是區(qū)域地貌對(duì)比研究的客觀依據(jù)［1］。近年來隨著遙感、地理信息系統(tǒng)以及全球定位系統(tǒng)等技術(shù)的快速發(fā)展，基于數(shù)字高程模型（Digital Elevation Model，DEM）的地形起伏度計(jì)算研究工作得以開展。國際地理學(xué)會(huì)地貌調(diào)查制圖委員會(huì)基于DEM數(shù)據(jù)采用16 km2的統(tǒng)計(jì)面積，計(jì)算地形起伏度并編制了歐洲1∶250萬地貌圖［2］。涂漢明等利用國家數(shù)字地形數(shù)據(jù)庫，結(jié)合中國地貌類型的基本特征，提出地形起伏度最佳統(tǒng)計(jì)單元應(yīng)滿足山體完整性與區(qū)域普適性兩條原則，論證了中國地形起伏度最佳統(tǒng)計(jì)單元為21 km2［3］;湯國安等基于DEM數(shù)據(jù)，使用窗口分析法計(jì)算多尺度地形起伏度［1］;劉新華等基于我國1∶100萬DEM數(shù)據(jù)，在全國采集6個(gè)樣區(qū)，運(yùn)用窗口分析法得出我國地形起伏度最佳統(tǒng)計(jì)窗口為25 km2［4］;唐飛等通過對(duì)準(zhǔn)噶爾盆地及其西北山區(qū)地形起伏度的研究，得出新疆克拉瑪依地區(qū)地形起伏度的最佳統(tǒng)計(jì)單元面積為4 km2［5］;郎玲玲等比較了福建低山丘陵區(qū)不同尺度DEM提取的地形起伏度，認(rèn)為基于1∶10萬、1∶25萬DEM計(jì)算地形起伏度的最佳統(tǒng)計(jì)單 元 分 別 為 0.41 km2、4.41 km2［6］;龍 恩 基 于SRTM-DEM數(shù)據(jù)對(duì)長白山區(qū)地形起伏度的最佳斷點(diǎn)進(jìn)行了判斷［7］;張磊基于1∶25萬DEM數(shù)據(jù)，采用鄰域分析和Focal函數(shù)計(jì)算了多尺度地形起伏度，實(shí)現(xiàn)了京津冀地區(qū)地貌類型的自動(dòng)劃分［8］;王玲等采用變點(diǎn)分析法處理1∶25萬DEM數(shù)據(jù)，實(shí)現(xiàn)了新疆地區(qū)地貌類型的自動(dòng)劃分［9］。傳統(tǒng)方法由于地形數(shù)據(jù)量大、計(jì)算繁瑣，使定量的地表形態(tài)、地貌研究發(fā)展緩慢，DEM的快速精確獲取和GIS技術(shù)的發(fā)展為基于DEM進(jìn)行地形起伏度提取與分析提供了有效手段。本文基于ASTER GDEM數(shù)據(jù)，探討一種高效率計(jì)算不同鄰域尺度單元下地形起伏度的方法。

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

研究區(qū)地處青藏高原東南邊緣與四川盆地西部邊緣山地交接地帶（北緯31°04′～35°64′、東經(jīng)96°92′～104°40′），橫跨青海、四川兩省，涉及瑪多、達(dá)日、班瑪、色達(dá)、壤塘、馬爾康、金川7縣，面積74 865 km2，呈西北東南走向的條帶狀，屬于高原與高山峽谷過渡地帶［10］，也是我國第一個(gè)階梯與第二個(gè)階梯的過渡地區(qū)，按地貌分為高原區(qū)、山原區(qū)和高山峽谷區(qū)，地勢(shì)自西北向東南逐漸降低。西北部高原區(qū)多為各水系上游切割地區(qū)，嶺平谷寬，山體較大，山坡較緩，地形起伏度較小;東南部高山峽谷區(qū)多為各水系切割劇烈地區(qū)，山脊多呈鋸齒狀，地勢(shì)險(xiǎn)要，地形垂直差異明顯［11，12］。研究區(qū)平均海拔 4 245 m，其中4 200～4 400 m的高程所占比重最高，最高峰是位于瑪多縣南部的巴顏喀拉山，海拔5 286 m［10］。

本研究采用ASTER GDEM數(shù)據(jù)，該數(shù)據(jù)是美國航空航天管理局（NASA）和日本經(jīng)貿(mào)及工業(yè)部（METI）2009年7月2日共同發(fā)布的全球空間分辨率為30 m的最新柵格型數(shù)字地形數(shù)據(jù)，是根據(jù)NASA新一代對(duì)地觀測(cè)衛(wèi)星TERRA的詳盡觀測(cè)結(jié)果制作完成的，數(shù)據(jù)覆蓋北緯83°至南緯83°之間的所有陸地區(qū)域。DEM的空間分辨率對(duì)地形起伏度的計(jì)算精度至關(guān)重要，在ASTER GDEM推出以前，可以獲取到的大范圍DEM數(shù)據(jù)是空間分辨率為90 m的SRTM DEM數(shù)據(jù)，該數(shù)據(jù)在小尺度地形提取時(shí)一直受到空間分辨率過低的影響。ASTER GDEM數(shù)據(jù)最大的優(yōu)勢(shì)在于其擁有較高的空間分辨率，在小尺度精細(xì)地形特征提取時(shí)明顯優(yōu)于SRTM DEM數(shù)據(jù);另外ASTER GDEM數(shù)據(jù)空間分辨率與TM遙感影像一致，為兩者的結(jié)合使用提供了便利條件。

2 地形起伏度計(jì)算

2.1 地形起伏度計(jì)算原理

地形起伏度（Relief Amplitude）是表征地表單元地勢(shì)起伏復(fù)雜程度的一個(gè)指標(biāo)，也是定量描述地貌形態(tài)、劃分地貌類型的重要指標(biāo)，其表達(dá)式如下［13］：

式中：RA表示地形起伏度，Hmax和Hmin分別是單位面積內(nèi)最大和最小高程值?？梢岳肁rcGIS用于鄰域分析的focal函數(shù)實(shí)現(xiàn)上式的計(jì)算，采用Focalmax和Focalmin分別計(jì)算單位面積內(nèi)的最大高程值和最小高程值，二者相減即可求得地形起伏度。

2.2 鄰域分析

鄰域分析又稱窗口分析，是一種局部運(yùn)算的空間分析，涉及一個(gè)中心柵格單元和一組環(huán)繞柵格單元，根據(jù)中心單元和環(huán)繞單元的柵格數(shù)值為中心單元位置生成一個(gè)新的函數(shù)值或鄰域統(tǒng)計(jì)量，即輸出柵格圖中每個(gè)格子的值均是輸入柵格圖中對(duì)應(yīng)格子鄰域單元的函數(shù)。在ArcGIS中基于鄰域函數(shù)的統(tǒng)計(jì)方法有10種，在地形起伏度計(jì)算中僅涉及Maximum和Minimum兩種鄰域分析方法。

focal函數(shù)的兩個(gè)重要參數(shù)（即鄰域形狀和鄰域單元大?。?duì)地形起伏度的影響至關(guān)重要，因此需要對(duì)鄰域分析問題展開深入探討。在ArcGIS軟件中基于鄰域分析函數(shù)包括空間分析模塊下Neighborhood Statistic工具和GRID模塊下的focal函數(shù)。由于focal函數(shù)可以通過AML語言編程反復(fù)調(diào)用，更適用于運(yùn)算量較大的計(jì)算，因此本文采用focal函數(shù)計(jì)算地形起伏度。運(yùn)用focal函數(shù)展開鄰域計(jì)算時(shí)，先審查當(dāng)前單元格鄰域內(nèi)的所有單元，此時(shí)鄰域相當(dāng)于一個(gè)濾波窗口。根據(jù)所選統(tǒng)計(jì)類型開展鄰域計(jì)算，并將統(tǒng)計(jì)結(jié)果賦給輸出單元，然后濾波窗口平移到下一個(gè)格子，重復(fù)上述過程，直到所有柵格計(jì)算完畢，生成新柵格圖層，鄰域計(jì)算完成。Arc GIS軟件提供了矩形、圓形、環(huán)形、扇形及自定義的不規(guī)則多邊形5種鄰域，在地形起伏度計(jì)算中多采用矩形或圓形鄰域［2］，本研究采用矩形鄰域。

2.3 多尺度地形起伏度計(jì)算

鄰域范圍直接影響地形起伏的計(jì)算精度，一般情況下，地形起伏度開始會(huì)隨著鄰域統(tǒng)計(jì)范圍的增大而顯著增大，當(dāng)增加到一定閾值后，其隨鄰域統(tǒng)計(jì)范圍增大不會(huì)發(fā)生顯著變化。因此基于鄰域分析方法計(jì)算地形起伏度時(shí)，應(yīng)將鄰域統(tǒng)計(jì)范圍逐步擴(kuò)大，計(jì)算不同尺度下的地形起伏度。目前最被認(rèn)可的方法是規(guī)則窗口遞增法［14］。根據(jù)式（1），以n×n（n＝3，7，11，…，167）增加步長為4的柵格矩形作為模版算子，3×3為起始窗口，對(duì)整個(gè)研究區(qū)進(jìn)行遍歷計(jì)算，到167×167終止。由于運(yùn)算量很大，因此采用AML語言編程調(diào)用focal函數(shù)來實(shí)現(xiàn)，程序設(shè)計(jì)流程如圖1所示。

圖1 地形起伏度計(jì)算程序設(shè)計(jì)流程Fig.1 Program design flow of calculating relief amplitude

3 結(jié)果分析

3.1 地形起伏度與鄰域單元變化分析

鄰域統(tǒng)計(jì)單元不同所計(jì)算出的地形起伏度差異顯著，而且多數(shù)鄰域單元的計(jì)算效果并不理想。因此定量分析地形起伏度與鄰域統(tǒng)計(jì)范圍關(guān)系是確定最佳鄰域尺度單元的基礎(chǔ)。將地形起伏度與鄰域統(tǒng)計(jì)范圍關(guān)系綜合成表1，在EXCEL中對(duì)表1中數(shù)據(jù)繪制散點(diǎn)圖并進(jìn)行Logarithmic對(duì)數(shù)方程擬合（圖2）。

表1 地形起伏度與鄰域面積關(guān)系Table 1 The relationship between relief amplitude and neighborhood area

圖2 地形起伏度與鄰域尺度擬合曲線Fig.2 The fitting curve of relief amplitude and neighborhood scales

圖2的擬合統(tǒng)計(jì)結(jié)果經(jīng)擬合優(yōu)度檢驗(yàn)，擬合度較好。擬合曲線總體變化規(guī)律是：起伏度值先隨統(tǒng)計(jì)格網(wǎng)單元面積增大而快速增大，當(dāng)格網(wǎng)單元面積達(dá)一定閾值后，其增大速度開始減緩，最后趨于平穩(wěn)，且在增速開始減慢過程中存有一明顯拐點(diǎn)，即為曲線由陡變緩的閾值，該值所對(duì)應(yīng)的鄰域統(tǒng)計(jì)單元即為最佳尺度單元。

3.2 最佳鄰域尺度單元

常規(guī)的確定最佳尺度單元方法是人工作圖判定法和最大高差法。人工作圖判定法雖然可以快速地判定曲線上的拐點(diǎn)，但其判定受人為主觀因素的影響較大，往往精度不高;最大高差法容易產(chǎn)生多個(gè)拐點(diǎn)，多作為人工作圖判定法的輔助方法。針對(duì)常規(guī)方法的不足，本文采用基于計(jì)算相鄰尺度最大起伏度差與鄰域面積差的比值的高差顯著性變化檢驗(yàn)法來確定最佳尺度單元。

基于正方形鄰域獲得的多尺度地形起伏度高差顯著性變化檢驗(yàn)公式為∶

式中：Ij表示地形起伏度隨面積變化的強(qiáng)度，Rj表示j×j尺度地形起伏度的最大值，Sj表示j×j尺度地形起伏度計(jì)算時(shí)的鄰域面積，（Sj＋4－Sj）/Sj＋4是將面積變化幅度標(biāo)準(zhǔn)化。Ij越大，表示單位鄰域增長面積對(duì)應(yīng)的地形起伏度變化越大;當(dāng)Ij最大時(shí)，即出現(xiàn)對(duì)應(yīng)的地形起伏度拐點(diǎn)，該點(diǎn)是地形起伏度所對(duì)應(yīng)的最佳面積單元。

依據(jù)上述公式，得到高差顯著性值隨尺度變化圖（圖3）?？梢钥闯?，正方形鄰域的最佳尺度單元位于第19個(gè)點(diǎn)，對(duì)應(yīng)鄰域單元大小為75×75個(gè)柵格單元，面積為5.0625 km2。采用同樣的計(jì)算與檢驗(yàn)方法，基于90 m空間分辨率的SRTM-DEM數(shù)據(jù)再次計(jì)算研究區(qū)地形起伏度，得到最佳鄰域尺度單元為14.29 km2，由此可見，DEM的精度對(duì)地形起伏度的計(jì)算至關(guān)重要。

圖3 變化強(qiáng)度Ij隨鄰域尺度變化Fig.3 The change of Ij with neighborhood scales

3.3 地形起伏度分級(jí)

參考中國1∶100萬數(shù)字地貌制圖規(guī)范中的分級(jí)標(biāo)準(zhǔn)［15］，對(duì)最佳尺度下的地形起伏度進(jìn)行分級(jí)，生成地形起伏度分級(jí)圖（圖4）。研究區(qū)地形起伏度的總體趨勢(shì)是：從西北向東南方向起伏度逐步增加?，敹唷⑦_(dá)日、色達(dá)3縣地勢(shì)較平坦，地形起伏度小于500 m，壤塘、馬爾康、金川3縣地勢(shì)起伏劇烈，絕大部分地區(qū)地形起伏度在500 m以上。從各級(jí)地形起伏度所占比例結(jié)構(gòu)看，200～500 m起伏度所占比例高達(dá)40.12%;其次是500～1 000 m，比例為26.05%;起伏度小于30 m的區(qū)域所占比例最小，僅有1.53%;起伏度大于1 000 m的區(qū)域有3 924.12 km2，占5.01%。由此可見，研究區(qū)地勢(shì)以中度起伏為主。

圖4 地形起伏度分級(jí)Fig.4 The class of relief amplitude

4 結(jié)論

基于ASTER GDEM數(shù)據(jù)和ArcGIS平臺(tái)AML編程，調(diào)用GRID模塊中用于鄰域分析的focal函數(shù)，可以高效計(jì)算不同鄰域尺度統(tǒng)計(jì)單元下地形起伏度，是一種有效的窗口遞增分析方法。鄰域統(tǒng)計(jì)單元大小對(duì)地形起伏度計(jì)算至關(guān)重要。兩者之間的擬合關(guān)系表明：起伏度值先隨鄰域統(tǒng)計(jì)單元面積增大而快速增大，當(dāng)鄰域統(tǒng)計(jì)單元面積達(dá)一定閾值后，其增大速度開始減緩，最后趨于平穩(wěn)，且在增速減慢過程中存在一明顯拐點(diǎn)。經(jīng)高差顯著性變化檢驗(yàn)法確定最佳鄰域統(tǒng)計(jì)單元面積為5.06 km2。研究區(qū)自西北向東南地形起伏度逐步增加，地勢(shì)以中度起伏為主。由于研究區(qū)面積較大，且ASTER GDEM數(shù)據(jù)空間分辨率高，導(dǎo)致地形起伏度計(jì)算過程運(yùn)算量很大，所以窗口分析的步長設(shè)置為4，如果將步長縮小，可能會(huì)得到更精確的結(jié)果。

［1］湯國安，楊瑋瑩，楊昕，等.對(duì)DEM地形定量因子挖掘中若干問題的探討［J］.測(cè)繪科學(xué)，2003，28（1）：28－32.

［2］湯國安，劉學(xué)軍，閭國年.數(shù)字高程模型及地學(xué)分析的原理與方法［M］.北京：科學(xué)出版社，2005.

［3］涂漢明，劉振東.中國地勢(shì)起伏度最佳統(tǒng)計(jì)單元的求證［J］.湖北大學(xué)學(xué)報(bào)（自然科學(xué)版），1990，12（3）：266－271.

［4］劉新華，張曉萍，楊勤科，等.不同尺度下影響水土流失地形因子指標(biāo)的分析與選?。跩］.西北農(nóng)林科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2004，32（6）：107－110.

［5］唐飛，陳曦，程維明，等.基于DEM的準(zhǔn)噶爾盆地及其西北山區(qū)地勢(shì)起伏度研究［J］.干旱區(qū)地理，2006，29（3）：388－392.

［6］郎玲玲，程維明，朱啟疆，等.多尺度DEM提取地勢(shì)起伏度的對(duì)比分析——以福建低山丘陵區(qū)為例［J］.地球信息科學(xué)，2007，9（6）：1－6.

［7］龍恩.基于Srtm-DEM與遙感的長白山基本地貌類型提取方法［J］.山地學(xué)報(bào)，2007，25（5）：557－565.

［8］張磊.基于地形起伏度的地貌形態(tài)劃分研究——以京津冀地區(qū)為例［D］.河北師范大學(xué)，2009.

［9］王玲，同小娟.基于變點(diǎn)分析的地形起伏度研究［J］.地理與地理信息科學(xué)，2007，23（6）：65－67.

［10］張學(xué)儒.基于地貌和植被的青藏高原東部樣帶土地類型研究［D］.中國科學(xué)院地理科學(xué)與資源研究所，2011.

［11］劉林山.黃河源地區(qū)高寒草甸退化研究——以達(dá)日縣為例［D］.中國科學(xué)院地理科學(xué)與資源研究所，2006.

［12］閻建忠，張鐿鋰，擺萬奇，等.大渡河上游生計(jì)方式的時(shí)空格局與土地利用/覆被變化［J］.農(nóng)業(yè)工程學(xué)報(bào)，2005，21（3）：83－89.

［13］朱紅春，陳楠，劉海英，等.自1∶10000比例尺DEM提取地形起伏度——以陜北黃土高原的實(shí)驗(yàn)為例［J］.測(cè)繪科學(xué)，2005，30（4）：86－88.

［14］周侗，龍毅，湯國安，等.面向DEM地形復(fù)雜度分析的分形方法研究［J］.地理與地理信息科學(xué)，2006，22（1）：26－30.

［15］中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室.中華人民共和國1∶100萬數(shù)字地貌制圖規(guī)范［Z］.2005.

Analysis on the Relief Amplitude Based on ASTER GDEM Data in Mountain Area of Eastern Tibetan Plateau

ZHANG Xue－ru，GUAN Dong－jie，MOU Feng－yun，CHEN Chun
（CollegeofRiverSea，ChongqingJiaotongUniversity，Chongqing400074，China）

Relief amplitude is an indicator that can reflect change of rolling surface topography.It is widely used in terrain mapping，evaluation of resources and the environment and so on.This paper takes mountain area of the eastern Tibetan Plateau as study area.The relief amplitude in different neighborhood scale were calculated by AML programming invoking focal functions in ArcGIS base on ASTER GDEM data with high spatial resolution.The researches show that：Neighborhood scales have a critical effect on relief amplitude.In beginning，the values of relief amplitude increase rapidly with growing in area of neighborhood statistics units.But when the area of neighborhood units reaches a certain threshold value，the growth rate begins to slow down and stabilize.And there is a significant inflection point that is best neighborhood-scale unit when the growth rate begins to slow down.The best neighborhood-scale unit for measuring relief amplitude is 5.0625 km2by significance test.Class of relief amplitude is got，and it is found that the relief amplitude gradually increases from northwest to southeast in study area.The relief amplitude between 200 m to 500 m has the highest proportion，about 40.12%，and level of 500～1 000 m accounts for 26.05%.Therefore fluctuation of terrain is mainly moderate in study area.

relief amplitude;ASTER GDEM;neighborhood analysis;mountain area of eastern Tibetan Plateau

P 931;P208

A

1672－0504（2012）03－0011－04

2011－11- 30;

2012－02－13

國家自然科學(xué)基金項(xiàng)目（40901057）;國家科技支撐計(jì)劃項(xiàng)目（2006BAB15B02-04）;河北省重點(diǎn)基礎(chǔ)研究項(xiàng)目（08966712D）

張學(xué)儒（1982－），男，博士，講師，主要從事土地利用與GIS應(yīng)用研究。E－mail：zhangxueru5＠126.com