大中流域尺度分布式土壤侵蝕學(xué)坡長(zhǎng)提取的數(shù)據(jù)分片方法

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

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

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

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

1 研究區(qū)概況

本研究分別在地形比較平緩的東北漫崗丘陵區(qū)和地形較陡的黃土丘陵區(qū)展開(kāi)研究。東北漫崗丘陵位于大小興安嶺和長(zhǎng)白山的山前臺(tái)地,由低平漫川和波狀起伏的丘陵組成,本研究樣區(qū)選在黑龍江克山和拜泉一帶。黃土丘陵區(qū)位于黃河中游流域,是我國(guó)乃至世界上土壤侵蝕最嚴(yán)重的地區(qū),本研究樣區(qū)選在位于典型黃土丘陵區(qū)的陜西安塞。據(jù)25 m分辨率坡度統(tǒng)計(jì),東北樣區(qū)和黃土樣區(qū)域的平均坡度分別是2.0°和21.5°。2個(gè)研究區(qū)地勢(shì)特征見(jiàn)圖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個(gè)研究區(qū)各選覆蓋9個(gè)1∶25萬(wàn)標(biāo)準(zhǔn)圖幅的DEM為基礎(chǔ)。所用DEM數(shù)據(jù)系根據(jù)1∶5萬(wàn)DLG(包括等高線、高程點(diǎn)和水系3個(gè)專題層),在ANUDEM支持下建立的水文地貌關(guān)系正確的DEM(hydrological correctly DEM,Hc-DEM)[9-10],分辨率為25 m[11]。

2 研究方法

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

2.1 流域坡長(zhǎng)概念模型

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

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

圖3 坡長(zhǎng)與分水線Fig.3 Slope length and watershed boundary

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

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

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

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

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

3)基于坡長(zhǎng)實(shí)際值推求:以流域?yàn)閱卧崛〉钠麻L(zhǎng),對(duì)坡長(zhǎng)做統(tǒng)計(jì)分布分析,用累計(jì)頻率99.99%作為最大值(LENmax/m)。

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

2.3 坡長(zhǎng)提取與對(duì)比

以流域?yàn)閱卧鸵詧D幅為單元(邊界緩沖后)提取坡長(zhǎng),并將結(jié)果進(jìn)行對(duì)比。主要步驟如下。

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

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

對(duì)比分析:對(duì)流域?yàn)閱卧蛨D幅為單元的坡長(zhǎng),從提取結(jié)果的圖形、統(tǒng)計(jì)特征和工作效率等方面進(jìn)行對(duì)比分析。

3 結(jié)果與分析

3.1 以流域?yàn)閱卧钠麻L(zhǎng)

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

圖5 東北漫崗丘陵樣區(qū)(左)和西北黃土丘陵樣區(qū)(右)地貌中等流域單元(面積分別為11.09萬(wàn)km2和13.42萬(wàn)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 坡長(zhǎng)特征分析 以下從空間格局和統(tǒng)計(jì)分布(圖6、表1)2個(gè)方面,對(duì)樣區(qū)坡長(zhǎng)做簡(jiǎn)單分析。

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

圖7 以流域?yàn)閱卧钠麻L(zhǎng)提取工作流程Fig.7 Working flowchart for extracting slope length by a watershed as unit

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

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

3.2 以圖幅為單元的坡長(zhǎng)

3.2.1 坡長(zhǎng)不完整帶的寬度 不完整流域帶寬度(Dw):根據(jù)匯流面積—河流密度關(guān)系分析,并結(jié)合對(duì)河流—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個(gè)研究區(qū)坡長(zhǎng)基本統(tǒng)計(jì)特征Tab.1 Statistical characteristics of two sites'slope lengths

坡長(zhǎng)實(shí)際值的統(tǒng)計(jì)分布:利用前述以流域?yàn)閱卧崛〉钠麻L(zhǎng),用流域邊界切割剔除流域周遍不完整坡長(zhǎng)部分,用累計(jì)頻率99.99%作為最大值(LENmax),結(jié)果在黃土和東北樣區(qū)坡長(zhǎng)最大值分別為491.4和2 124.5 m(表3)。

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

表3 實(shí)際坡長(zhǎng)統(tǒng)計(jì)表Tab.3 Statistics of actual slope length

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

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

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

圖8 標(biāo)準(zhǔn)圖幅坡長(zhǎng)提取流程Fig.8 Flowchart of extracting slope length from standard map sheet

3.3 2種方法比較

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

圖9 流域?yàn)閱卧蜆?biāo)準(zhǔn)圖幅為單元坡長(zhǎng)直方圖Fig.9 Slope length histogram while a watershed as an unit and a standard map sheet as an unit

據(jù)我們工作記錄,主要工作步驟的消耗時(shí)統(tǒng)計(jì)見(jiàn)表5。以流域?yàn)閱卧钠麻L(zhǎng)提取總耗時(shí)46 min,以標(biāo)準(zhǔn)圖幅為的單元提取坡長(zhǎng)總耗時(shí)8 min,為常規(guī)方案的17.4%；所以,以標(biāo)準(zhǔn)圖幅為單元提取坡長(zhǎng),效率遠(yuǎn)大于以中等流域?yàn)閱卧墓ぷ鞣绞健?/p>

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

4 結(jié)論與討論

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

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

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

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

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

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

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

5 參考文獻(xiàn)

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

[2] 楊勤科,郭偉玲,張宏鳴,等,基于DEM的流域坡度坡長(zhǎng)因子計(jì)算方法研究初報(bào)[J].水土保持通報(bào),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] 王程,陳正江,楊勤科,等.流域分布式坡長(zhǎng)不確定性的初步分析[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].中國(guó)水土保持科學(xué),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] 國(guó)家測(cè)繪局.基礎(chǔ)地理信息數(shù)字產(chǎn)品1∶10 000、1∶50 000數(shù)字高程模型:CH/T 1008—2001[S].北京:中國(guó)標(biāo)準(zhǔn)出版社,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ū)域坡度坡長(zhǎng)因子提取算法[J].計(jì)算機(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.