曹佳云，楊勤科，2，王 程，丑述仁

（1.西北大學(xué) 城市與環(huán)境學(xué)院，西安710127；2.中國(guó)科學(xué)院 水利部 水土保持研究所，陜西 楊凌712100）

全國(guó)第四次土壤侵蝕調(diào)查已經(jīng)在全國(guó)范圍內(nèi)展開，區(qū)域土壤侵蝕因子是主要的調(diào)查內(nèi)容。由于實(shí)用型區(qū)域土壤侵蝕模型尚在開發(fā)中，USLE和RUSLE被廣泛應(yīng)用于土壤侵蝕調(diào)查制圖［1］。USLE和RUSLE中地形對(duì)土壤流失的影響用坡度坡長(zhǎng)因子（LS）表示［2－3］。LS計(jì)算是利用USLE等坡面模型進(jìn)行流域與區(qū)域土壤侵蝕評(píng)價(jià)與制圖的基礎(chǔ)［4］。國(guó)內(nèi)外研究者針對(duì)流域和區(qū)域土壤侵蝕調(diào)查制圖中的LS計(jì)算進(jìn)行了一系列探索。但是在這些研究中，多集中在一般方法的討論，且多以流域?yàn)閱卧鐓菛|亮、張照錄、汪邦穩(wěn)等人的研究［5－7］。國(guó)外 Moore，Wilson，Williams和Desmet等專門討論了流域LS計(jì)算方法［8－12］。Hickey 和 Van Remortel等 研 究 了 基 于DEM 提取區(qū)域LS的方法［13－16］；在此基礎(chǔ)上，張宏鳴等設(shè)計(jì)開發(fā)了區(qū)域尺度LS因子計(jì)算工具［17］。已有研究中多以流域?yàn)閱卧阅承姓^(qū)為單元提取LS因子的方法較少涉及；較多關(guān)注了丘陵區(qū)，對(duì)于高塬溝壑區(qū)注意不夠。本研究以黃土高原溝壑區(qū)為典型研究區(qū)，探討縣域地形因子提取和分析方法，為第四次全國(guó)土壤侵蝕普查及其后續(xù)工作提供技術(shù)支撐。

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

1.1 研究區(qū)概況

本次研究的區(qū)域選在長(zhǎng)武縣，東西跨越107°38′49″—107°58′02″E，長(zhǎng)度約為27.23km，南北跨越34°59′09″—35°18′37″N，長(zhǎng)度約為30.06km，面積567.1 km2，黃土塬、梁峁丘陵等各種地貌發(fā)育比較齊全。

1.2 數(shù)據(jù)基礎(chǔ)

本次研究的基礎(chǔ)數(shù)據(jù)來源為國(guó)家基礎(chǔ)比例尺1∶5萬地形圖，所需圖幅為長(zhǎng)武縣界緩沖2km后所包含的圖幅。地形圖等高距為10m，高斯克呂格6°投影。從測(cè)繪局購(gòu)買并自行按照國(guó)家技術(shù)規(guī)程數(shù)字化，并手工勾繪地形特征線。

2 LS因子提取流程

以長(zhǎng)武縣國(guó)家基礎(chǔ)比例尺1∶5萬地形圖為基礎(chǔ)數(shù)據(jù)，首先通過ANUDEM［18］軟件建立水文地貌關(guān)系正確的 DEM（Hydrology Correct DEM，簡(jiǎn)稱Hc—DEM），在此基礎(chǔ)上通過自行開發(fā)的LS工具［17］從DEM中提取LS因子，操作流程見圖1。

圖1 LS提取流程圖

2.1 Hc－DEM建立

通過地形圖建立DEM，首先要對(duì)地形數(shù)據(jù)進(jìn)行編輯處理，消除其中存在的錯(cuò)誤，并手工勾繪地形特征線；然后將等高線、高程點(diǎn)、河流、湖泊、流域邊界以及地形特征線輸入ANUDEM軟件，建立Hc—DEM。

2.1.1 地形數(shù)據(jù)編輯 采用ANUDEM軟件建立DEM，需要輸入的基礎(chǔ)地形數(shù)據(jù)包括：等高線、高程點(diǎn)、河流和湖泊。這些地形數(shù)據(jù)可以在1∶5萬地形圖中獲取。

要獲取全縣范圍的等高線、高程點(diǎn)、河流和湖泊數(shù)據(jù)，首先需要進(jìn)行標(biāo)準(zhǔn)分幅地形圖的拼接和投影變換。將數(shù)據(jù)拼接成一個(gè)圖幅后，需要對(duì)地形數(shù)據(jù)進(jìn)行檢查和錯(cuò)誤編輯，主要包括：檢查等高線和高程點(diǎn)錯(cuò)誤，保證高程有序變化；檢查和編輯河流流向，確保河流由高往低處流；湖泊位置檢查，使其位于低洼部位；格式轉(zhuǎn)換：對(duì)于編輯好的數(shù)據(jù)，用ungenerate命令生成AUNDEM可以識(shí)別的文本文件。

2.1.2 地形特征線提取 在等高線基礎(chǔ)上，手工勾繪主要的地形特征線，主要包括：塬邊線、溝沿線以及坡腳線。轉(zhuǎn)換生成地形特征線的gen文件。關(guān)于特征線提取，詳見文獻(xiàn)［19］。

2.1.3 DEM建立 利用編輯好的地形數(shù)據(jù)，運(yùn)行ANUDEM生成縣域DEM。輸入文件為等高線、高程點(diǎn)、河流、湖泊、緩沖的縣域邊界以及地形特征線。主要參數(shù)設(shè)置為，迭代次數(shù)20，第二粗糙度0.5。輸出文件為2進(jìn)制或文本格式，可在ArcGIS中轉(zhuǎn)化為grid格式。由于數(shù)據(jù)量比較大，如果運(yùn)行中斷，可分塊逐個(gè)運(yùn)行生成DEM，然后再用Mosaic功能合并得到縣域DEM。

2.1.4 流域劃分 流域劃分通過AML編程方式自動(dòng)實(shí)現(xiàn)，然后對(duì)流域邊界進(jìn)行編輯和輸出。根據(jù)本研究需要，本著各流域單元面積適中、形狀盡量接近圓形的原則，將長(zhǎng)武縣劃分為4個(gè)流域單元。將各流域邊界單獨(dú)存為一個(gè)文件。對(duì)各流域邊界建立1km的緩沖區(qū)，得到緩沖的流域邊界，并將緩沖邊界文件轉(zhuǎn)換生成gen文件。

2.2 LS因子提取

在生成的Hc－DEM的基礎(chǔ)上，利用自行開發(fā)的LS因子提取工具［17］，提取坡度、坡長(zhǎng)，根據(jù)第四次土壤侵蝕普查技術(shù)規(guī)程中的規(guī)定，計(jì)算LS因子。

2.2.1 坡度、坡長(zhǎng)和LS因子算法 自行開發(fā)的LS計(jì)算工具［17］，坡度提取采用D8算法，坡長(zhǎng)采用單流向下的徑流累計(jì)算法。其中坡長(zhǎng)的提取須以流域?yàn)閱卧?/p>

2.2.2 LS因子計(jì)算的實(shí)現(xiàn) 在ArcGIS環(huán)境下將grid格式DEM轉(zhuǎn)換為文本文件，然后輸入LSTools中，每個(gè)流域?qū)⑤敵?個(gè)文本文件，分別是坡度、坡長(zhǎng)、坡度因子、坡長(zhǎng)因子和坡度坡長(zhǎng)因子。這些文本文件可在ArcGIS環(huán)境下轉(zhuǎn)化為grid，并將各個(gè)流域的5個(gè)因子圖進(jìn)行拼接，生成縣域完整的LS因子系列圖，包括坡度、坡向、坡度因子、坡長(zhǎng)因子、坡度坡長(zhǎng)因子。

3 結(jié)果與分析

3.1 DEM與地形宏觀結(jié)構(gòu)特征

利用地形圖上的地形信息和地形特征線生成的DEM，從宏觀上表達(dá)該地區(qū)塬、塬坡、現(xiàn)代溝坡、川地等各種地貌結(jié)構(gòu)，較之常規(guī)方法建立的DEM，對(duì)侵蝕地形的表達(dá)能力有所改善［19］（圖2）。

圖2 長(zhǎng)武縣典型樣區(qū)Hc－DEM

3.2 坡度、坡長(zhǎng)和LS值的空間分布

從生成的坡度、坡長(zhǎng)、LS因子的分布圖（圖3）可以看出，在溝道、河灘、塬的頂部，坡度較小、坡長(zhǎng)較長(zhǎng)，LS值較??；在溝道兩邊的坡上，坡度較大、坡長(zhǎng)較小，LS值較大。LS值的分布與坡度的分布最相似，而與坡長(zhǎng)的分布差別較大。這是因?yàn)樵邳S土高原丘陵區(qū)，坡度對(duì)侵蝕的影響大于坡長(zhǎng)［20］。

圖3 坡度、坡長(zhǎng)、LS因子圖

3.3 坡度、坡長(zhǎng)和LS值的頻率分布

3.3.1 整個(gè)區(qū)域坡度、坡長(zhǎng)和LS值的頻率分布

將整個(gè)縣的坡度、坡長(zhǎng)和LS因子的特征值進(jìn)行統(tǒng)計(jì)，將坡度、坡長(zhǎng)和LS因子的屬性表導(dǎo)出，在Excel中制作各因子的頻率與累積頻率分布曲線，如圖4所示。從頻率分布曲線中可以看出，0～25°的坡度占總坡度的80%左右，其中坡度值為0～5°的坡度占了50%。坡度曲線在0.5°附近最高，先驟降，隨后在20°～25°出現(xiàn)小的峰值，最后開始下降，到60°附近基本持平，60°～80°的坡度只占總坡度的0.05%左右。這是因?yàn)樵搮^(qū)域存在比較多的塬、河灘等平坦地貌，有坡度較陡的區(qū)域存在是因?yàn)镈EM構(gòu)建時(shí)添加了地形特征線。坡長(zhǎng)的分布中，小于500m的坡長(zhǎng)占了95%，說明黃土高原丘陵區(qū)一方面地形比較破碎復(fù)雜，相對(duì)起伏度較大，另一方面說明該區(qū)域也包含較多的平坦地貌。LS的頻率曲線從值為1的最高值開始驟降，到值為5左右逐漸下降。它是坡度與坡長(zhǎng)共同影響的結(jié)果，與坡度頻率分布更加相近，這也驗(yàn)證了在黃土高原丘陵區(qū)，坡度對(duì)侵蝕的影響大于坡長(zhǎng)的影響。從坡度、坡長(zhǎng)和LS因子的特征值與頻率曲線可以得出，提取出的因子可以符合該地區(qū)的實(shí)際地形狀況。

圖4 坡度、坡長(zhǎng)與LS因子頻率分布曲線

圖5 各類地形的坡度分布

圖6 各類地形的坡長(zhǎng)分布

3.3.2 各地形類型區(qū)坡度、坡長(zhǎng)和LS值的頻率分布 在整個(gè)縣選擇4個(gè)典型樣區(qū)（塬面、塬坡、現(xiàn)代溝坡和川地），在提取的坡度、坡長(zhǎng)和LS圖基礎(chǔ)上裁剪出4個(gè)典型區(qū)，統(tǒng)計(jì)其地形特征值的平均值，結(jié)果表明，各個(gè)地形類型中，現(xiàn)代溝坡坡度最陡，頻率曲線偏向右（高值區(qū)）；塬面最緩，平均不到1°，頻率曲線偏向左（低值區(qū)）；塬坡和川地介于中間（圖5）。各類地形坡長(zhǎng)分布可分為3組（圖6），塬面和塬坡最長(zhǎng)（大于全區(qū)坡長(zhǎng)），現(xiàn)代溝坡最短（小于全區(qū)坡長(zhǎng)），而川地介于其間（大于全區(qū)坡長(zhǎng)）。

4 結(jié)論

（1）縣域行政單元LS因子提取的主要技術(shù)環(huán)節(jié)為基礎(chǔ)數(shù)據(jù)編輯、流域劃分、地形特征線提取、用ANUDEM軟件生成Hc—DEM以及用LS＿Tools提取LS因子。

（2）對(duì)于一個(gè)行政單元或流域而言，良好的基礎(chǔ)數(shù)據(jù)、先進(jìn)的DEM插值方法、科學(xué)合理的LS因子提取技術(shù)，是LS因子提取的基本要求，保證這幾個(gè)方面是LS因子提取的關(guān)鍵技術(shù)環(huán)節(jié)。

（3）以流域?yàn)閱卧腖S因子提取，是針對(duì)LS因子提取的科學(xué)原理而言的；以縣域?yàn)閱卧腖S因子提取，是針對(duì)LS因子提取的實(shí)用而言的；縣域LS因子的提取，也必須遵循LS因子提取的科學(xué)原理，注意避免因?yàn)镈EM可辨識(shí)的最小流域不完整而帶來的邊際效應(yīng)。

（4）利用我們的前期理論研究和開發(fā)的LS因子提取工具，可提取縣級(jí)行政單元的LS因子專題層和相關(guān)單項(xiàng)要素，如坡度、坡長(zhǎng)、流域邊界等。從提取結(jié)果的統(tǒng)計(jì)和空間分布特征看，基本可反映該地區(qū)的地貌特征。

（5）長(zhǎng)武塬區(qū)各不同地貌單元的坡度、坡長(zhǎng)特征有明顯不同，坡度表現(xiàn)為塬面坡度比較平緩，塬坡和川地介于中間，現(xiàn)代溝坡的坡度最陡；坡長(zhǎng)表現(xiàn)為塬面和塬坡最長(zhǎng)，川地其次，現(xiàn)代溝坡坡長(zhǎng)最短。

［1］ 楊勤科，李銳，劉詠梅.區(qū)域土壤侵蝕普查方法的初步討論［J］.中國(guó)水土保持科學(xué)，2008，6（3）：1-7.

［2］ Renard K G，F(xiàn)oster G R，Weesies G A.Predicting Rainfall Eosion by Water：A Guide to Conservation Planning with the Revised Universal Soil Loss Equation（RUSLE）［M］.Washington D C：USDA Agric.Handb，1997.

［3］ Wischmeier W H，Smith D D.Predicting Rainfall Eosion Losses from Cropland East of the Rocky Mountains：A Guide for Soil and Water Conservation Planning［M］.Washington D C：USDA Agric.Handb.，1978.

［4］ 楊勤科，郭偉玲，張宏鳴，等.基于GIS和DEM的流域坡度坡長(zhǎng)因子計(jì)算方法［J］.水土保持通報(bào)，2010，30（2）：203-206.

［5］ 汪邦穩(wěn)，楊勤科，劉志紅，等.基于DEM和ArcGIS的修正通用土壤流失方程的地形因子值提?。跩］.中國(guó)水土保持科學(xué)，2007，5（2）：18-23.

［6］ 吳東亮，劉鵬舉，唐小明.基于GIS的柵格化坡面徑流路徑模擬與LS值計(jì)算［J］.北京林業(yè)大學(xué)學(xué)報(bào)，2001，23（5）：10-14.

［7］ 張照錄，崔繼紅.基于柵格GIS土壤侵蝕地形因子的提取算法［J］.計(jì)算機(jī)工程，2006，32（5）：226-228.

［8］ Williams J R，Berndt H D.Determining the Universal Soil Loss Equation′s Length-slope Factor for Watersheds［M］Ankeny I A.Erosion and Solid Matter Transport in Inland Waters，IAHS-AISH Publication No.122，1977.

［9］ 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-428.

［10］ 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-184.

［11］ Desmet P，Grovers G.GIS-based simulation of erosion and deposition patterns in an agricultural landscape：A comparison of model results with soil map information［J］.Catena，1995，25（1／4）：389-401.

［12］ Desmet P J.A GIS procedure for automatically calculting the USLE LS factor on topographically complex landscape units［J］.Journal of Soil and Water Conservation，1996，51（5）：427-433.

［13］ Hickey R，Smith A，Jankowski P.Slope length calculations from a DEM within ARC／INFO grid［J］.Computers，Environment and Urban Systems，1994，18（5）：365-380.

［14］ Hickey R.Slope angle and slope length solutions for GIS［J］.Cartography，2000，29（1）：1-8.

［15］ Van Remortel R，Hamilton M，Hickey R.Estimating the LS factor for RUSLE through iterative slope length processing of digital elevation data［J］.Cartography，2001，30（1）：27-35.

［16］ Van Remortel R，Maichle R W，Hickey R J.Computing the LS factor for the Revised Universal Soil Loss Equation through array-based slope processing of digital elevation data using a C＋＋executable［J］.Computers ＆ Geosciences，2004，30（9／10）：1043-1053.

［17］ 張宏鳴，楊勤科，郭偉玲，等.基于GIS的區(qū)域LS因子算法及實(shí)現(xiàn)［J］.計(jì)算機(jī)工程，2010，36（9）：246-248.

［18］ Hutchinson M.ANUDEM Version 5.2 ［EB／OL］.［2011-06-15］.http：∥fennerschool.anu.edu.au／research／publications／software-datasets／anudem.

［19］ 古云鶴，楊勤科，羅儀寧，等.突變地形特征在DEM上的表達(dá)［J］.水土保持研究，2011，18（2）：174-179.

［20］ 唐克麗.中國(guó)水土保持［M］.北京：科學(xué)出版社，2004.