褚永彬，朱利東，楊文光，卞玉霞，崔林林

(1.成都信息工程大學資源環(huán)境學院，成都 610224；2.成都理工大學沉積地質(zhì)研究院，成都 610051)

地形是構造、營力和時間的函數(shù)[1]，其中的營力指的是外營力作用。作為一種地表分布廣泛的外營力，流水作用具有顯著的削高填洼效果。在一個流域中，流水總是將高處的物質(zhì)搬運到低處。水流的侵蝕、搬運、堆積，塑造了流水地貌。隨著時間的推移，流水地貌從幼年期到老年期的演化，被搬運的物質(zhì)越來越多，殘留的物質(zhì)越來越少。如何定量地描述殘留物質(zhì)的多少，反映流水地貌的演化階段，成為地貌研究的一個科學問題。針對這一問題，Strahler[2]于1952年提出了面積高程積分(hypsometric integral，HI)的概念，記錄流域某一條等高線的相對高程與流域高程差之比以及此條等高線的水平橫截面積與流域面積之比，組成一條曲線來代表集水盆地在侵蝕作用下殘留的物質(zhì)比例，用曲線的凹凸形狀表達流域演化階段。由于河流對巖性、構造、氣候的敏感響應，使得河流的地貌過程發(fā)生瞬時改變，面積高程積分也相應改變?；诖?，通過HI的分析來探討巖性、構造和氣候的變化，尤其是構造活動，成為河流地貌研究的熱點。在構造分析方面，Parisa等[3]利用包括HI在內(nèi)的地貌指標討論了阿爾伯茲山的構造活動性；Siddiqui等[4]分析了HI對研究意大利北部亞平寧山脈活動構造的意義，揭示了高HI值與構造活躍和區(qū)域構造正相關；Gong等[5]利用HI再現(xiàn)了梅江流域地貌特征與構造響應之間的關系；楊帆等[6]、高澤民等[7]分別利用HI分析了墨脫地區(qū)、河套盆地北緣大青山地區(qū)的構造地貌特征，為當?shù)氐姆勒饻p災提供參考；邵崇建等[8]、洪艷等[9]基于HI計算結果，討論了龍門山地區(qū)的構造地貌特征以及構造活動強弱。在巖性分析方面，張威等[10]計算并分析了他念他翁山玉曲流域的面積高程積分值與構造活動以及巖性之間的關系。在氣候?qū)Φ匦斡绊懛矫?，辛聰聰?shù)萚11]討論了雅魯藏布江下游的構造活動、氣候與面積高程積分值之間的關系；Taloor等[12]則討論了構造活動、氣候變化對喜馬拉雅南麓Durung Drung盆地地貌演化的影響?？梢哉f，面積高程積分已成為定量地貌的重要方法之一，其應用研究文獻較多，但有關計算方法的報道則很少。目前應用最多的是Pike等[13]在1971年提出的高程起伏比法，此方法通過簡易的公式計算面積高程積分值，無法繪制曲線，不能直觀展示流域演化階段。

基于此，利用常用軟件ArcGIS和Excel，設計了一種簡易方法，實現(xiàn)了面積高程曲線繪制和積分值計算，并以祁連山大通河流域為實驗區(qū)進行了方法實用性的驗證。

1 面積高程積分曲線的概念、原理

“水往低處流”這是自然規(guī)律。在一個流域中，受侵蝕、搬運作用，地表物質(zhì)被剝蝕，高程逐漸降低。伴隨著流域從幼年期到老年期的不斷演化，低海拔的地表面積占流域面積的比值逐漸增大，高海拔的地表面積占比則逐漸減少。基于此，采用面積-高程曲線將流域中某一條等高線的水平橫截面積與等高線到流域出水口的相對高程聯(lián)系起來，如圖1[2]所示。為了得到無量綱值便于進行比較，面積高程曲線中縱坐標采用等高線的相對高度h與流域地勢高差H之比，即y=h/H；橫坐標采用等高線所切的水平斷面面積a與流域總面積A之比，即x=a/A。面積高程曲線下方與坐標軸之間的面積值即為面積高程積分值。如圖2[1]所示。面積高程積分可以反映流域內(nèi)物質(zhì)相對總量、流域發(fā)育進程以及流域勢能[14]。

圖1 面積高程積分的定義[2]Fig.1 Definition of Hypsometric Integral[2]

圖2 面積高程積分曲線示意圖[1]Fig.2 Sketch map for Hypsometric Integral curve[1]

2 數(shù)據(jù)和方法

2.1 所需數(shù)據(jù)

數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)是地貌分析重要的數(shù)據(jù)來源。常用的DEM數(shù)據(jù)，如SRTM、ASTER-GDEM的數(shù)據(jù)結構為柵格數(shù)據(jù)結構，在ArcGIS中，柵格數(shù)據(jù)具有對應的柵格屬性表。當柵格屬性表生成時，該表中會創(chuàng)建3個默認字段：OID、VALUE和COUNT，其中VALUE為柵格數(shù)據(jù)集中列出的唯一單元值(以格網(wǎng)形式表示，均為整數(shù))，COUNT為柵格數(shù)據(jù)集中具有某個“值”列單元值的單元數(shù)量[15]，OID(Object ID)為系統(tǒng)定義的針對表格中每行的唯一對象標識符編號，如圖3所示。

圖3 柵格數(shù)據(jù)屬性表Fig.3 The graphic of Raster dataset attribute table

為了更好地進行說明，選擇青海省大通河某次級流域為例。圖4為該次級流域的DEM。圖5為此流域的DEM局部放大?？芍车雀呔€所切的水平斷面面積可以轉(zhuǎn)換成大于等于此等高線高程的柵格投影面積，即柵格個數(shù)柵格單元長柵格單元寬。

圖4 大通河某次級流域DEM三維顯示Fig.4 3D map of a certain sub watershed of Datong River

圖5 大通河某次級流域DEM局部放大圖Fig.5 Localenlarged drawing of a certain sub watershed of Datong River

2.2 提取方法及流程

由面積高程積分的概念可知，繪制面積高程曲線需要某一等高線Hi的相對高程h及此高程以上的流域面積a，流域的地勢高差H以及流域面積A。提取面積高程積分曲線的關鍵在于得到高程及位于此高程以上的流域面積。

由圖3可知，從屬性表中可以獲得高程值Hi對應的柵格個數(shù)。柵格單元格的大小是確定的，則高程值對應的面積計算公式為

SHi=R2CHi

(1)

式(1)中：R為柵格大小，R2為單個柵格的面積；CHi為高程值Hi對應的柵格個數(shù)；SHi為高程值Hi對應的面積。那么位于高程值Hi以上的柵格面積則為

(2)

a=Si

(3)

相對高程可表示為

h=Hi-Hmin

(4)

式(4)中：Hmin為流域內(nèi)最小高程值。

在ArcGIS中，可以通過水文分析提取流域以及流域范圍內(nèi)的DEM數(shù)據(jù)，構建其柵格屬性表，已知屬性表中記錄了高程唯一值及對應的柵格個數(shù)(圖6)，將屬性表導出為文本文件并用Excel表格打開，然后按VALUE字段進行降序排序(圖7)。增加a、A、h、H、x、y，共計6個字段(圖8)。利用Excel公式對6個字段分別賦值：

圖6 示例流域的柵格屬性表(部分)Fig.6 A partial view of the Raster dataset attribute table of an example watershed

圖7 示例流域的柵格屬性EXCEL表(部分)Fig.7 A partial view of the Raster dataset attribute excel table of an example watershed

a為高程值Hi及其以上部分的面積，需要獲取柵格個數(shù)，Excel函數(shù)為“=SUM($B$2:B2)”，即統(tǒng)計COUNT所在列的高程值Hi當前單元格與第一個高程值之間的個數(shù)總和；A為流域面積，需要獲取全流域的柵格個數(shù)，Excel函數(shù)為“=SUM($B$2:$B$1756)”；h為高程值Hi的相對高度，即高程值Hi與所有高程值中最小值之間的差值，Excel函數(shù)為“=A2-MIN($A$2:$A$1756)”；H為流域地勢高差，即所有高程值中最大值與最小值之間的差值，Excel函數(shù)為“=MAX($A$2:$A$1756)-MIN($A$2:$A$1756)”；則x=a/A，y=h/H，得到橫坐標和縱坐標的值(圖8)，即可繪制面積高程曲線，如圖9所示。曲線呈典型的S形，根據(jù)曲線形態(tài)可以判斷該流域地形屬于壯年期(均衡)階段。

圖8 面積高程積分參數(shù)計算結果Fig.8 Calculation results of Hypsometric Integral

圖9 示例流域的面積高程曲線Fig.9 Hypsometric Curve of an example watershed

2.3 面積高程曲線積分值的簡易計算

常直楊等[16]比較了3種HI計算方法：起伏比法、積分曲線法、體積比例法，對比結果顯示不同方法計算結果幾乎一致，其中起伏比法是計算HI值最高效最便捷的方法。起伏比法是Pike等[13]于1971年推導出的HI簡易計算方法，其簡化公式為

(5)

ArcGIS分區(qū)統(tǒng)計工具可以方便地統(tǒng)計出流域內(nèi)高程的最大值Hmax、最小值Hmin、均值Hmean，由式(5)計算HI。在示例流域中，最大值Hmax=4 274，最小值Hmin=2 191，均值Hmean=3 112，則計算出的HI值為0.44，同樣顯示流域地形屬于均衡的壯年期地形。

3 祁連山地貌分析應用實例

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

為了驗證上述方法的實用性，祁連山大通河流域作為實驗區(qū)，提取了12個次級流域的面積高程積分值并繪制了積分曲線。大通河發(fā)源于祁連山東段大通山和托勒南山之間，是湟水最大的支流，也是祁連山東段最大的縱向河流之一。大通河流域受若干斷裂控制，所在區(qū)域具有低速率的抬升，HI為0.49，剛剛進入壯年期階段。

考慮到HI對DEM分辨率不具有依賴性[17]，故而采用高程精度高于ASTER-GDEM的SRTM3數(shù)據(jù)[18]。SRTM3是研究人員廣泛使用的地形數(shù)據(jù)，其空間分辨率90 m。

3.2 方法及結果

利用ArcGIS水文分析提取大通河次級流域。選擇了其中12個流域。利用分割工具將此12個流域面要素分割為12個單獨的面圖層文件，通過按掩膜提取批處理的方式，完成流域DEM的提取，并分別對流域DEM進行屬性表構建，屬性表導出及計算操作，最后完成HI的計算和曲線繪制，結果如圖10、圖11、表1所示。

表1 大通河12個次級流域HI值及演化階段Table 1 Values and evolution stages of 12 sub basins of Datong River

1～12為次級流域編號圖10 祁連山大通河12個次級流域分布Fig.10 Distribution of 12 sub basins of Datong River in Qilian Mountains

圖11 祁連山大通河12個次級流域HI曲線Fig.11 Hypsometric Curves of 12 sub basins of Datong River

大通河上游流經(jīng)大通河盆地，其次級流域面積高程積分曲線下凹，流域演化已進入老年期階段(流域編號1～4)，與整個大通河流域處于壯年期的演化階段不協(xié)調(diào)。其余河段次級流域雖處于壯年期階段，但不同河段的積分值依然存在差異，相較于流域7～12，流域5、流域6曲線較平直。次級流域面積高程積分的差異，說明上游河段演化早于中下游，結合斷裂的分布，也表明中下游受構造運動影響較大，河流演化較為復雜。流域7和流域9具有相同的HI值，但曲線卻略有差異。曲線的差異表明在次級小流域的上游，流域9的侵蝕程度高于流域7；在下游，流域9堆積的物質(zhì)量大于流域7。

4 結論

(1)DEM柵格屬性表具有高程值對應的柵格數(shù)目，為面積高程積分的計算提供了基礎。

(2)基于DEM的特點，設計了基于ArcGIS和Excel的面積高程積分提取技術流程，可方便地完成面積高程積分曲線的繪制及積分值的計算。計算結果與常用的高程起伏比法一致。

(3)對青海省大通河流域SRTM數(shù)字高程模型數(shù)據(jù)進行了實驗研究和定量分析。研究表明，該方法可較為方便地完成面積高程積分曲線的繪制及積分值的計算，準確反映流域發(fā)育程度及對構造活動的響應。

(4)該方法無需借助其他軟件，無需編程，簡便易操作，便于面積高程積分方法的廣泛使用，對地形地貌信息提取和定量分析具有一定價值。