曹夏雨，楊勤科，2?，蘭 敏，王春梅，2

(1.西北大學城市與環(huán)境學院,710069,西安；2.中國科學院水利部水土保持研究所,712100,陜西楊凌；3.中煤西安設(shè)計工程有限公司,710054,西安)

土壤侵蝕預報研究中,坡長是指地表從徑流源點開始,沿徑流線到達到坡度減小直至有沉積出現(xiàn)地方之間或者到一個明顯的溝道(自然的溝道或人工渠道)之間的水平距離[1]。針對流域內(nèi)每個點(DEM柵格)提取的坡長稱為流域分布式土壤侵蝕學坡長(簡稱流域坡長)[2]。

流域坡長的計算,是基于USLE(或RUSLE、CSLE)進行流域和區(qū)域尺度土壤侵蝕評價制圖的最關(guān)鍵問題之一[3-5]。其計算受到諸多因素影響,主要包括DEM分辨率、流向算法和基礎(chǔ)數(shù)據(jù)的范圍等[6]。關(guān)于數(shù)據(jù)分片方面,研究者強調(diào)須以完整的流域為單元,以便能在完整的徑流路徑上計算坡長,以避免完整的坡面被數(shù)據(jù)邊界截斷,使數(shù)據(jù)邊界處的部分坡長被遺失[7-8]。事實上,這種認識是針對小流域或大中流域較粗分辨率(數(shù)據(jù)量較小)的坡長提取而言的；而對于較高分辨率、較大面積(數(shù)千到上萬個標準圖幅組成)的情況下,坡長提取的數(shù)據(jù)如何分片(數(shù)據(jù)單元),卻未見討論。筆者對流域坡長提取的數(shù)據(jù)單元劃分方法做出探討,以便避免或減少邊際效應(yīng)前提下,在較大區(qū)域范圍內(nèi)快速提取坡長,滿足區(qū)域土壤侵蝕評價與制圖的需要。

1 研究區(qū)概況

本研究分別在地形比較平緩的東北漫崗丘陵區(qū)和地形較陡的黃土丘陵區(qū)展開研究。東北漫崗丘陵位于大小興安嶺和長白山的山前臺地,由低平漫川和波狀起伏的丘陵組成,本研究樣區(qū)選在黑龍江克山和拜泉一帶。黃土丘陵區(qū)位于黃河中游流域,是我國乃至世界上土壤侵蝕最嚴重的地區(qū),本研究樣區(qū)選在位于典型黃土丘陵區(qū)的陜西安塞。據(jù)25 m分辨率坡度統(tǒng)計,東北樣區(qū)和黃土樣區(qū)域的平均坡度分別是2.0°和21.5°。2個研究區(qū)地勢特征見圖1。

圖1 東北漫崗丘陵樣區(qū)(左)和西北黃土丘陵樣區(qū)(右)地形Fig.1 Landform of the Northeast China undulated hill site(Left,NE site for short below)and the Northwest China loess hill site(Right,NW site for short below)

本研究以2個研究區(qū)各選覆蓋9個1∶25萬標準圖幅的DEM為基礎(chǔ)。所用DEM數(shù)據(jù)系根據(jù)1∶5萬DLG(包括等高線、高程點和水系3個專題層),在ANUDEM支持下建立的水文地貌關(guān)系正確的DEM(hydrological correctly DEM,Hc-DEM)[9-10],分辨率為25 m[11]。

2 研究方法

本研究基于坡長提取的一般原理,利用數(shù)字地形分析方法,結(jié)合對水系、流域和坡長統(tǒng)計特征的分析來完成。

2.1 流域坡長概念模型

流域坡長的計算,從局部高點(位于分水地帶)開始,順徑流路徑向下坡方向?qū)γ總€柵格單元的長度不斷累加,直到發(fā)生泥沙沉積的地方或明顯的溝道(圖2)[7,12]。從典型小流域坡長提取結(jié)果看,分水線往下,坡長的值逐漸增加(圖3,由白到黑坡長增加)；所以流域坡長提取的基本數(shù)據(jù)單元是DEM上可以辨認的、具有水文地貌意義的最小流域—最小有效流域(WME),只要數(shù)據(jù)以WME的邊界或相應(yīng)級別的流水線為邊界,即可提取完整的坡長。

圖2 坡長提取技術(shù)流程Fig.2 Flowchart of extracting slope length

圖3 坡長與分水線Fig.3 Slope length and watershed boundary

2.2 基于WME的坡長提取數(shù)據(jù)單元

據(jù)坡長概念模型,較大范圍矩形數(shù)據(jù)塊劃分的WME,其最外圍流域邊界(或流水線)的連線(BNDWME)將接近于矩形(圖4中粗黑色折線)。BNDWME之內(nèi)為一組完整小流域,之外為一個不完整小流域帶,也是不完整坡長帶。假定圖4中內(nèi)部黑色矩形為某標準圖幅的圖廓,則只要對標準圖幅向外緩沖寬度大于不完整坡長帶寬度(Dwr)后,即可以規(guī)則圖幅(包括緩沖帶)為數(shù)據(jù)單元計算得到完整的坡長。Dwr通過下面3種方法確定。

圖4 規(guī)則數(shù)據(jù)塊和不規(guī)則流域邊界Fig.4 Rectangular data and incomplete watershed zone

1)直接量算:即在一個工作區(qū)選擇典型樣區(qū),提取最小有效流域,勾繪近似矩形的BNDWME,量取其最大內(nèi)接和最小外接矩形間的最大距離(Dw/m)。

2)基于溝壑密度推求:在DEM上提取與WME適應(yīng)的河流,統(tǒng)計河流總長度,用式(1)推算平均坡長[13]。

3)基于坡長實際值推求:以流域為單元提取的坡長,對坡長做統(tǒng)計分布分析,用累計頻率99.99%作為最大值(LENmax/m)。

最終,規(guī)則數(shù)據(jù)向外擴充(緩沖)的數(shù)據(jù)寬度(Dwr/m)為上述各值的最大值。

2.3 坡長提取與對比

以流域為單元和以圖幅為單元(邊界緩沖后)提取坡長,并將結(jié)果進行對比。主要步驟如下。

流域劃分:經(jīng)過填洼、流向計算、累計流量計算等步驟,通過在ARC/INFO workstation下編程,選擇合適匯水面積閾值,提取最小有效流域,用以確定不完整坡長帶寬度；然后將樣區(qū)9個1∶25萬標準圖幅的DEM拼接為一個數(shù)據(jù)單元后,將其劃分為20個左右的中等流域,用以作為流域坡長提取的基本單元。

坡長提取:分別以流域和以標準圖幅(以Dwr為寬度向外緩沖)為單元,利用筆者開發(fā)的坡度坡長提取工具軟件提取坡長[2,14]。

對比分析:對流域為單元和圖幅為單元的坡長,從提取結(jié)果的圖形、統(tǒng)計特征和工作效率等方面進行對比分析。

3 結(jié)果與分析

3.1 以流域為單元的坡長

3.1.1 流域劃分與坡長提取 東北和黃土樣區(qū)1∶25萬標準分幅、25 m分辨率DEM劃分的較大流域單元如圖5示。為了提取正中間矩形圖框(1∶25萬標準圖幅)部分的坡長,將涉及的中等流域(圖5中的11—14,21—24)邊界,經(jīng)過緩沖后切割出4個流域的數(shù)據(jù)塊,然后逐一提取坡長,拼接后從中間切割出工作區(qū)的坡長(圖6,均為1∶25萬標準圖幅的左上角)。以流域為單元的坡長提取工作流程如圖7。有6個基本步驟。

圖5 東北漫崗丘陵樣區(qū)(左)和西北黃土丘陵樣區(qū)(右)地貌中等流域單元(面積分別為11.09萬km2和13.42萬km2)Fig.5 Medium watershed units in the NE site(Left)and NW site(Right)(Areas are 11.09×104km2and 13.42×104km2respectively)

3.1.2 坡長特征分析 以下從空間格局和統(tǒng)計分布(圖6、表1)2個方面,對樣區(qū)坡長做簡單分析。

圖6 東北漫崗丘陵樣區(qū)(上)西北黃土丘陵樣區(qū)(下)坡長Fig.6 Slope length of NE site(Up)and NW site(Down)

圖7 以流域為單元的坡長提取工作流程Fig.7 Working flowchart for extracting slope length by a watershed as unit

1)東北漫崗丘陵區(qū)坡長:坡長平均值479 m,中值320 m,最大值3 570 m。空間上坡度比較平緩的漫崗地坡度最長,最大值出現(xiàn)在漫崗與漫川的轉(zhuǎn)折部位；緩坡丘陵坡度較短。

2)黃土丘陵區(qū)的坡長:樣區(qū)平均坡長86.1 m,中值68.4 m,最大值1 223.9 m。空間分布是,從丘陵頂部向下增加,到坡面向溝道或川地交界處理最長達到最大。

3.2 以圖幅為單元的坡長

3.2.1 坡長不完整帶的寬度 不完整流域帶寬度(Dw):根據(jù)匯流面積—河流密度關(guān)系分析,并結(jié)合對河流—DEM表面關(guān)系的觀察,東北漫崗丘陵區(qū)和黃土丘陵區(qū)在25 m分辨率DEM上匯水面積閾值分別為1.0和0.5 km2。根據(jù)前述方法量算,在東北漫崗丘陵區(qū)和黃土丘陵區(qū),Dw的值分別為4 791.6 m和2 776.5 m。

表1 2個研究區(qū)坡長基本統(tǒng)計特征Tab.1 Statistical characteristics of two sites'slope lengths

坡長實際值的統(tǒng)計分布:利用前述以流域為單元提取的坡長,用流域邊界切割剔除流域周遍不完整坡長部分,用累計頻率99.99%作為最大值(LENmax),結(jié)果在黃土和東北樣區(qū)坡長最大值分別為491.4和2 124.5 m(表3)。

表2 基于溝壑密度的坡長統(tǒng)計Tab.2 Slope length statistics based on gully density

表3 實際坡長統(tǒng)計表Tab.3 Statistics of actual slope length

坡長數(shù)據(jù)緩沖寬度(Dwr):據(jù)以上統(tǒng)計,黃土樣區(qū)和東北樣區(qū)不完整坡長帶寬度分別是2 776.5和4 791.6 m(表4),實際應(yīng)用中,可設(shè)置為3和5 km。

表4 矩形圖幅周邊不完整流域帶參數(shù)統(tǒng)計Tab.4 The statistic parameters of the incomplete watershed around the rectangular map sheet

3.2.2 坡長提取流程 對標準圖幅向外緩沖3 000(黃土樣區(qū))和5 000 m(東北樣區(qū)),然后在LS_Tool系統(tǒng)中提取流域坡長,工作流程見圖8。與圖7相比,工作步驟由6步減少到3步。尤其重要的是,減少了流域劃分、按流域切割數(shù)據(jù)和提取坡長后的再拼接,從而使工作效率大為提高。

圖8 標準圖幅坡長提取流程Fig.8 Flowchart of extracting slope length from standard map sheet

3.3 2種方法比較

把相同樣區(qū)2種方法提取的坡長做差值運算,結(jié)果全部等于0,表明坡長的值完全一樣。由于表面一致,因而使其統(tǒng)計直方圖(比較上部位兩條線為黃土樣區(qū)坡長累計頻率)也完全一致(圖9)。有此可見,以標準圖幅為單元提取坡長是可行的。

圖9 流域為單元和標準圖幅為單元坡長直方圖Fig.9 Slope length histogram while a watershed as an unit and a standard map sheet as an unit

據(jù)我們工作記錄,主要工作步驟的消耗時統(tǒng)計見表5。以流域為單元的坡長提取總耗時46 min,以標準圖幅為的單元提取坡長總耗時8 min,為常規(guī)方案的17.4%；所以,以標準圖幅為單元提取坡長,效率遠大于以中等流域為單元的工作方式。

表5 流域坡長和規(guī)則數(shù)據(jù)塊坡長計算時間比較Tab.5 Comparison of calculation times between two calculating methods(Watershed slope length and rectangular data) min_

4 結(jié)論與討論

流域坡長提取的數(shù)據(jù)單元,是一個長期被忽視的問題,本研究結(jié)論及有關(guān)問題歸納如下。面向大區(qū)域、多各圖幅和海量數(shù)據(jù)的坡長提取數(shù)據(jù)單元研究,可得出如下初步結(jié)論。

1)根據(jù)流域分布式土壤侵蝕學坡長的含義和提取算法,本文提出流域坡長概念模型—流域坡長提取的基本數(shù)據(jù)單元是DEM上可以辨認的、具有水文地貌意義的最小流域。據(jù)此,對于一個由多個圖幅構(gòu)成的較大工作區(qū),只要對將每個圖幅范圍向外擴展一定距離(Dwr),即可直接用規(guī)則圖幅為單元逐幅提取坡長。

2)基于對不完整流域帶寬度(Dw)、溝壑密度和坡長均值、和坡長實際值的統(tǒng)計(LENmax),在地形較平緩東北樣區(qū)和地形較陡黃土樣區(qū),標準圖幅向外緩沖的寬度分別為2 776.5和4 791.6 m。

3)基于中等流域和標準圖幅的坡長提取,其結(jié)果完全一致。工作步驟分別為6步和3步,總耗時間分別為45和8 min。所以,基于標準圖幅提取大范圍坡長,工作步驟大為簡化,效率大為提高。

筆者論述了以標準圖幅或規(guī)則矩形塊為數(shù)據(jù)單元提取流域坡長的基本原理,也提出了標準圖幅為單元提取坡長應(yīng)該緩沖的寬度推算方法；然而以下問題有待深入研究:

1)坡長提取運算的數(shù)據(jù)量上限:經(jīng)試驗或者理論推導,推算出在給定的硬件條件下利用LS_Tool工具提取坡長時,可順利運行數(shù)據(jù)的DEM最大數(shù)據(jù)量(矩陣行列數(shù)),以便更有效劃分數(shù)據(jù)單元。

2)不完整坡長帶的推算:上述方法雖然可推算出不完整坡長帶的寬度,但是原則上要根據(jù)工作區(qū)域的地形特征做出率定和優(yōu)化。今后的工作中,應(yīng)在足夠樣本基礎(chǔ)上,選擇合適參數(shù),建立比較通用的模型以便能推算出不完整坡長帶寬度值。

5 參考文獻

[1] SMITH D D,WISCHMEIER W H.Factors affecting sheet and rill erosion[J].Trans.Am.Geophys.Union,1957.38(6):889.

[2] 楊勤科,郭偉玲,張宏鳴,等,基于DEM的流域坡度坡長因子計算方法研究初報[J].水土保持通報,2010.30(2):203.YANG Qinke,GUO Weiling,ZHANG Hongming,et al.Method of extracting LS factor at watershed scale based on DEM[J].Bulletin of Soil and Water Conservation,2010,30(2):203.

[3] MOORE I.D,BURCH G J.Physical basis of the lengthslope factor in the universal soil loss equation[J].Soil Science Society of America Journal,1986.50(5):1294.

[4] WILSON J P.Estimating the topographic factor in the universal soil loss equation for watersheds[J].Journal of Soil and Water Conservation,1986.41(3):179.

[5] MOORE I D,WILSON J P.Length-slope factors for the revised universal soil loss equation:Simplified method of estimation[J].Journal of Soil and Water Conservation,1992.47(5):423.

[6] 王程,陳正江,楊勤科,等.流域分布式坡長不確定性的初步分析[J].水土保持研究,2012,19(2):15.WANG Cheng,CHEN Zhengjiang,YANG Qinke,et al.Analysis on uncertainty of DEM derived watershed Dia[J].Research of Soil and Water Conservation,2012,19(2):15.

[7] HICKEY R,A SMITH,AND P JANKOWSKI.Slope length calculations from a DEM within ARC/INFO GRID[J].Computers,Environment and Urban Systems,1994.,18(5):365.

[8] GRUBER S,PECKHAM S.Land-surface parameters and objects in hydrology[M]//HENGL T,REUTER H I.Geomorphometry:Concepts,Software,Applications,Series Developments in Soil Science Vol.33.Amsterdam:Elsevier,2009:207.

[9] YANG Qinke,McVICAR T R,Van NIEL T G,et al.Improving a digital elevation model by reducing source data errors and optimising interpolation algorithm parameters:an example in the Loess Plateau,China[J].International Journal of Applied Earth Observation and Geoinformation),2007,9(3):235.

[10] 楊勤科,師維娟,McVicar T R,等.水文地貌關(guān)系正確的DEM建立方法[J].中國水土保持科學,2007,5(4):1.YANG Qinke,SHI Weijuan,McVICAR T R,et al.On constructing methods of hydrologically correct DEMs[J].Science of Soil and Water Conservation,2007,5(4):1.

[11] 國家測繪局.基礎(chǔ)地理信息數(shù)字產(chǎn)品1∶10 000、1∶50 000數(shù)字高程模型:CH/T 1008—2001[S].北京:中國標準出版社,2002:3.National Administration of Surveying.CH/T 1008- 2001 Digital products of fundamental geographic information 1∶10 000,1∶50 000 digital elevation models:CH/T 1008-2001[S].Beijing:Standards Press of China,2002:3.

[12] HICKEY R.Slope angle and slope length solutions for GIS[J].Cartography,2000.29(1):1.

[13] 陸中臣,賈紹鳳,黃克新,等著.流域地貌系統(tǒng)[M].大連:大連出版社,1991:32.LU Zhongchen,JIA Shaofeng,HUANG Kexin,et al.Watershed landform system[M].Dalian:Dalian Publishing House,1991:32.

[14] 張宏鳴,楊勤科,劉晴蕊,等,基于GIS的區(qū)域坡度坡長因子提取算法[J].計算機工程,2010,36(9):246.ZHANG Hongming,YANG Qinke,LIU Qingrui,et al,Regional slope length and slope steepness factor extraction algorithm based on GIS,[J].Computer Engineering,2010,36(9):246.