秦偉,朱清科,左長清,趙磊磊,鄺高明,王昭艷,張京鳳,單志杰
(1.中國水利水電科學研究院泥沙研究所,100048;2.北京林業大學 水土保持與荒漠化防治教育部重點實驗室,100083:北京)
地表過程(earth surface processes),即地球表面大氣圈、生物圈與巖石圈之間因物質與能量的交互作用而形成的一系列物理與化學過程,主要包括巖土物質的侵蝕、搬運和沉積,水熱養分的運移、轉化和循環,以及土地利用演變及其環境效應等。地表過程是地球表面固有的自然過程,由此引發地貌格局演變,水系發育變遷,并影響生態環境,同時,地表過程與人類活動密切相關、相互影響。在全球氣候變化與人類活動加劇的背景下,地表過程將更加深刻地影響人類的可持續發展。全面揭示地表過程規律,有效地向避害趨利的方向調控,以促進人與自然和諧相處,已為世界普遍關注,諸如土壤侵蝕、山地災害、生態水文等典型地表過程也因此長期成為學界研究焦點。由于地表過程的復雜性,通常需確立一定標準,據此將研究范圍劃分為一系列基本單元,先就基本單位完成模擬或分析,再經累加或串聯實現整個研究范圍的模擬以及相關特征的揭示。這種化繁為簡、化整為零的思路是解決地理學、生態學、水文學等眾多學科復雜問題的重要途徑。其中,所劃分的基本單元,即為本文的地表過程響應單元。地表過程響應單元劃分是否合理、獲取是否高效很大程度上決定地表過程研究結果是否可靠、過程是否順利。現有的劃分主要包括規則單元和不規則單元2 類。規則單元是將研究區域按一定分辨率劃分為均勻的方格單元,即柵格單元。這種方法簡單、快速、易實現,且在地理信息系統技術廣泛應用的背景下,能更好地與遙感影像等多源數據結合分析。目前采用最多,但其僅考慮了地形等環境因素影響地表過程的空間差異性,忽略了一定空間范圍內下墊面影響地表過程的相似性。不規則單元是根據地形地貌特點、水文運移路徑或土地利用類型等確定特征邊界將研究區域劃分為不規則的面狀班塊。這種方法在特定角度更真實反映了研究對象的固有屬性,更利于后期分析研究,但相比規則單元劃分,其實現難度更大,尤其針對較大空間尺度時,往往難以快捷地獲取預期結果。比較有代表性的不規則相應單元包括水文響應單元(Hydrological Response Unit,HRU)[1]、水文相似單元(Hydrological Similar Unit,HSU)[2]和分組相應單元(Grouped Response Unit,GRU)[3]等。
黃土丘陵溝壑區是世界上地形最破碎復雜的地區之一,加之干旱與暴雨交疊的氣候條件,造成該區土壤侵蝕、滑坡坍塌、生態水文等為代表的地表過程具有顯著的特殊性和復雜性。針對該區特點,建立一種合理、高效的地表過程響應單元劃分方法對推動地表過程研究具有重要的基礎意義。筆者以黃丘陵溝壑區典型流域為例,擬提出針對土壤侵蝕為代表的地表過程模擬研究的響應單元劃分方法及其在GIS 平臺中的實現技術,以期為促進該區地表過程研究提供有益參考。
四面窯溝流域位于陜北吳起縣鐵邊城鎮、北洛河上游南岸,E 107°45'47″~107°49'24″,N 36°54'08″~36°59'58″,屬北洛河水系,面積32.68 km2。流域內海拔1 369 ~1 730 m,溝壑縱橫,屬典型黃土高原丘陵溝壑地貌類型(圖1)。土壤主要為地帶性黑壚土剝蝕后廣泛發育在黃土母質上的黃綿土,質地為輕壤。多年平均降雨量478.3 mm,降雨年際變化大,季節分配不均,7—9 月降雨量占年均降雨量的62.4%,且多為暴雨。受耕作等人類活動影響,流域內天然植被退化嚴重,現存植被主要為人工次生植被,以沙棘(Hippophae rhamnoides)、檸條(Caragana mocrophylla)等灌木,河北楊(Populus hopeiensis)、小葉楊(Populus simonii)、山杏(Prunus sibirica)、杜梨(Pyrus betulaefolia)等闊葉喬木以及紫花苜蓿(Medicago sativa)等人工牧草為主,林草覆蓋率約80%。

圖1 研究流域區位圖Fig.1 Location of the study watershed
除氣候條件外,地表過程主要受下墊面影響。劃分響應單元的核心在于針對下墊面影響地表過程的空間異質性,分類、合并地表過程差異顯著的下墊面,形成內部地表過程特征相同或相似的下墊面單元。土地利用與地形地貌是下墊面的主要內容。其中,土地利用分類標準和劃分邊界比較明確,即按不同土地利用類型及其相連的界線進行劃分。地形地貌所包含的具體內容更為豐富,不同地形地貌指標均能不同程度影響地表過程,且以不同地形地貌指標為依據的劃分邊界也不盡清晰;因此,響應單元劃分的難點主要在于地形地貌特征解析及其邊界提取。就黃土丘陵溝壑區的土壤侵蝕而言,坡度、坡向、坡位是影響土壤侵蝕的主要地形地貌特征。坡面作為組成該區地貌結構的基本組份,其內部的坡度和坡向雖存在空間分異,但并不顯著,獨立坡面在整體上具有基本一致的坡向和坡度。坡位對土壤侵蝕的影響集中表現在梁坡和溝坡間土壤侵蝕強度和機制的差異,即梁坡以薄層徑流的剝蝕為主,侵蝕強度較低;溝坡以含沙股流的沖蝕為主,侵蝕強度較高。因此,以小流域為對象,參照水文研究中所提及的水文響應單元[1],將土壤侵蝕響應單元定義為由流域溝谷線、山脊線、分水線及溝緣線交疊構成的具有相對均一坡向、坡度與坡位的獨立坡面,并在Arc-GIS9.2 中,以數字高程模型為基礎,實現響應單元劃分及其地形特征提取。
數字高程模型(Digital Elevation Model,DEM)作為真實地貌的離散數字表達,是水文、土壤、地貌、生態等領域的重要基礎數據,分為規則格網(如grid、raster 等)和不規則三角網(Triangulated Irregular Network,TIN)[4]。通常多基于數字化等高線生成不規則三角網模型,然后再轉化為規則格網模型;然而,因內插計算誤差或少數真實地形影響,由此獲得的DEM 表面將存在局部凹陷洼地,阻斷匯水路徑從而影響后期的地表過程分析[5]。為此,采用陜西省1∶1萬紙質地形圖,掃描后在R2V 軟件中完成數字化,進行平滑、賦值等編輯處理,形成矢量等高線;然后在ArcGIS9.2 中,利用Spatial Analysis 模塊的Topo to Raster 工具,基于ANUDEM 算法生成水文關系正確的數字高程模型(Hydrologically Correct DEM,Hc-DEM)。Hc-DEM 消除了低洼等局部干擾地形,能形成通暢水文路徑,可為后期分析奠定基礎。
2.2.1 溝谷線、山脊線及分水線提取 溝谷線、山脊線及分水線不僅是流域地形地貌的分界線,也是流域水文循環過程中的特征線。在ArcGIS9.2 中,提取溝谷線、山脊線及分水線均以數字高程模型為基礎數據,在Hydrology 模塊中完成。
1)溝谷線提取。利用Hydrology 模塊的Flow Direction 工具和Flow Accumulation 工具,以Hc-DEM 為基礎數據,先后獲得網格單元的匯流方向和累積匯流量;再利用Raster Calculator 工具,輸入不同累積匯流閾值選擇生成不同級別的溝谷線。輸入的累積匯流閾值越小,獲得的溝谷密度越大,長度越長。流域內的溝谷系統是線狀水流作用形成的不同等級和規模的溝谷總稱。根據其形狀、大小和發育階段等特征,通常分為細溝、淺溝、切溝、沖溝和坳溝等5 類[6],內部相應存在以細溝侵蝕、淺溝侵蝕、切溝侵蝕、沖溝侵蝕和河道侵蝕為主的侵蝕類型。由于細溝、淺溝多存在于獨立坡面內,且難由一般精度的數字高程模型中提取,切溝、沖溝多存在于溝緣線上下,結構不穩定,不宜用于劃分響應單元;因此,以溝道和河道為提取對象,通過將不同累積匯流閾值對應的提取結果與高分辨率遙感影像(SPOT-5,2.5 m×2.5 m 分辨率,下同)進行對比,最終將累積匯流閾值0.6 km2時提取的溝谷作為研究流域的輸沙主溝道和永久河道。
2)山脊線與分水線提取。山脊線與分水線的提取方法與溝谷線一致。由于是匯流起點,因此,累積匯流閾值為零的柵格所組成的線即山脊線與分水線。其中,分水線是流域最外圍的連續閉合山脊線。2.2.2 溝緣線提取 黃土丘陵溝壑區從分水線到溝底,地貌形態存在明顯地帶分異,依次包括梁峁頂部、梁峁坡、溝坡和溝底。其中,梁峁頂部和梁峁坡構成溝間地(梁坡),溝坡和溝底構成溝谷地(溝坡),溝間地與溝谷地之間以溝緣線為界(圖2)。溝緣線上下的土壤、地形和植被存在差異,并由此導致土壤侵蝕等地表過程具有不同的特點。作為正負2大地形分區的界線,溝緣線是黃土丘陵溝壑區地表過程研究中十分重要的地貌特征線。

圖2 黃土丘陵溝壑區典型坡面橫剖面示意圖Fig.2 Structure feature of typical slope profile in loess hilly and gully region
溝緣線提取包括利用等高線或高分辨率航片、衛片手工勾繪和基于數字化高程模型或遙感影像的編程提取。手工勾繪精度較高,但主觀性大、效率低,難以滿足大、中空間范圍的研究需要;因此,基于特定數據源的編程自動提取應用越來越廣。以Hc-DEM 為數據源,在ArcGIS9.2 中將其轉為txt 文本格式作為編程計算對象。編程思路為:根據柵格單元間的匯水路徑,逐個計算從分水線出發到河道的每條匯水路徑上單元的坡度變化率,將坡度變化率最大的單元作為溝緣線組成點;在此基礎上,排除局部破碎點和獨立點,形成一條聯通、完整的溝緣線;將結果與高分辨率遙感影像進行比對,確定提取結果精度。具體步驟如下。
1)溝緣線上下的地形差異主要表現為坡度陡變,溝緣線位置通常對應坡面坡度變化率最大的部位,因此,首先以3×3 的窗口為范圍,采用三次加權差分法計算中心單元坡度,再以同樣方法獲得其與周圍8 個單元的坡度變化率(單元格間位置編號見圖3),以作為溝緣線提取依據:

式中:S 為單元格坡度,(°),S128為周邊8 個單元中某單元的坡度,其他同;fx為x 方向高程變化率;fy為y 方向高程變化率;Z 為高程,m,Z128為周邊8 個單元中某單元的高程,其他同;Sr0為單元格坡度變化率;g 為單元格大小,即柵格分辨率,m。
2)采用ArcGIS9.2 中Hydrology 模塊的Flow Direction 工具,輸入Hc-DEM,獲得各單元匯流方向;然后由分水線出發,按流向至河道結束的每條匯流路徑上選取坡度變化率最大的柵格,作為該路徑(坡面縱斷面)上的溝緣線位置點;將該點上、下的單元格分別標記為溝間地單元和溝谷地單元,并相應賦值(1、2)。

圖3 DEM 3×3 局部移動窗口Fig.3 3×3 moving part window of DEM
3)坡度變化率篩選和單元類型標記后可初步形成溝緣線,但因局部高地的存在,使個別位置存在單個或幾個孤立單元格構成的孤島,影響流域溝緣線的聯通性。為此,對上述孤立點進行排除處理,獲得完整、聯通的溝緣線以及由其劃分的流域溝間地和溝谷地分區。采用ArcGIS9.2 中Conversion 模塊的From Raster 工具,將柵格數據轉換為矢量線、面數據,并與高分辨率遙感影像比對,判斷提取結果的準確性。結果顯示,以上算法提取的溝緣線與實際地貌吻合較好,與遙感影像上清楚可辨的溝緣線相比,誤差不超過±10 m,即提取誤差在1 個單元格內,誤差較小,結果可靠(見圖4)。

圖4 流域溝緣線提取結果Fig.4 Extraction of shoulder line of gully
2.2.3 地表過程響應單元劃分 在ArcGIS9.2 中,利用Hydrology 模塊,將基于Hc-DEM 提取的流域山脊線和溝谷線按水流方向反向延長與分水線相交;疊加編程提取的溝緣線,并將各地貌特征線進行合并,形成地表過程響應單元[7](圖5)。

圖5 流域地表過程響應單元劃分結果Fig.5 Division of response units to earth surface processes
完成地表相應單元劃分后,可按以下方法提取其地形特征。
1)響應單元平均坡度提取。在ArcGIS9.2 中,直接提取的坡度是單元格與周圍單元格間的最大坡降,屬局部坡度,不能作為響應單元平均坡度。因此,采用3D Analysis Tool 的剖面線Interpolate Line功能,以Hc-DEM 為數據源,在響應單元內以分水線/山脊線為起點,溝緣線/溝谷線為終點,沿順坡方向創建帶有高程信息的直線,并利用Create Profile Graph 功能生成直線對應的剖面圖;根據剖面圖中所包含的直線長度和直線兩端的高程差計算坡面坡度[8]。根據規格,每個響應單元創建5 至10 條直線,將平均值作為對應響應單元的平均坡度。
2)響應單元整體坡向提取。采用3D Analysis Tool 的Suface Analysis 功能,提取坡向,并劃分為陰坡(337.5°~360°和0°~67.5°)、半陰坡(67.5°~112.5°和292.5°~337.5°)、陽坡(157.5°~247.5°)和半陽坡(112.5°~157.5°和247.5°~292.5°)4 個坡向類別,分別記為1、2、3、4。將響應單元內多數柵格單元的坡向類別作為其整體坡向,以消除局部地形影響。
3)響應單元順坡長度提取。在響應單元內,采用Measure 功能量算分水線/山脊線與溝緣線/溝谷線間沿垂直等高線方向的距離,除以所在響應單元平均坡度的余弦,獲得順坡實際長度。根據規格,每個響應單元創建5 至10 條直線,將平均值作為對應響應單元的順坡長度。
4)響應單元面積與周長提取。采用Conversion Tool 的To Coverage 功能將獲得的響應單元圖層轉化為Coverage 格式,則可在屬性表中自動增加各單元的周長與面積。
與傳統規則單元劃分法相比,針對黃土丘陵溝壑區土壤侵蝕的地形地貌分異特點,提出的在流域內以溝谷線、山脊線、分水線和溝緣線相互交疊構成的,具有相對均一坡向、坡度與坡位的地表過程響應單元劃分方法,在反映黃土丘陵溝壑區不同地形地貌分區侵蝕空間異質性的基礎上,保留了相對完整的地形地貌斑塊,能更好地與針對坡面的原位試驗觀測、基于地塊的水土保持措施配置相結合,在水土保持工程、林業生態工程和山地災害防治等領域具有較廣闊的應用前景。
針對所提出的響應單元,在ArcGIS9.2 中,以Hc-DEM 為基礎,分別確定的溝谷線、山脊線、分水線和溝緣線的提取方法,以及平均坡度、整體坡向、面積與周長等地形特征的測算技術,可實現地表過程響應單元的快速劃分及其地形特征高效提取,從而可克服不規則單元劃分通常不易實現的問題。
土壤侵蝕是黃土丘陵溝壑區的典型地表過程,土壤侵蝕模擬與預報是該學科的主要技術難點和熱點。基于不規則單元的模型構建和模擬分析,能拓展土壤侵蝕預報、模擬的研究思路,有助于推動該學科領域的發展。
[1] Beven K.Linking parameters across scales: sub-grid parameterization and scale dependent hydrological models[J].Hydrological process,1995,9:507-525
[2] 周大良.基于遙感信息的林冠截流模型[D].北京:北京大學,1997:28-30
[3] 任立良.數字流域模擬研究[J].河海大學學報,2000,28(4):1-7
[4] Ao T Q,Takeuchi K,Ishidaira H,et al.Development and application of a new algorithm for automated pits removal for grid DEMs[J].Hydrological Sciences Journal,2003,48(6):985-997
[5] 周云貴,劉瑜,鄔倫.基于數字高程模型的水系提取算法[J].地理學與國土研究,2000,16(4):77-81
[6] 左大康.現代地理學辭典[M].北京:商務印書館,1990:222
[7] 秦偉.北洛河上游土壤侵蝕特征及其對植被重建的響應[D].北京:北京林業大學,2009:96-97
[8] 秦偉,朱清科,趙磊磊,等.基于RS 和GIS 的黃土丘陵溝壑區淺溝侵蝕地形特征研究[J].農業工程學報,2010,26(6):58-64