方磊,劉文兆,2,李懷有

(1.西北農(nóng)林科技大學(xué)資源環(huán)境學(xué)院,陜西楊陵712100;2.中國(guó)科學(xué)院水利部水土保持研究所,陜西楊陵712100;3.黃河水利委員會(huì)西峰水土保持科學(xué)試驗(yàn)站,甘肅 慶陽745000)

地形形態(tài)分布特征對(duì)地表徑流和土壤侵蝕、植被分布和土地利用、水土保持措施布設(shè)等具有重要影響[1],它直接影響著地表物質(zhì)的遷移與能量的轉(zhuǎn)換,標(biāo)志著地表基本環(huán)境格局的形成與發(fā)展。

基于數(shù)字高程模型(Digital Elevation M odel,簡(jiǎn)稱DEM)的數(shù)字地形分析技術(shù)經(jīng)過40多年的發(fā)展,其諸多基礎(chǔ)理論問題得到了深入的研究[2],而利用GIS中DEM的自動(dòng)提取技術(shù)已經(jīng)成為地形因子的主要提取手段。同時(shí)DEM的應(yīng)用也從簡(jiǎn)單的地形要素提取擴(kuò)展到更廣的范圍,如輔助建立水文模型、土壤侵蝕模型和洪水預(yù)報(bào)模型等。目前國(guó)內(nèi)學(xué)者對(duì)基于DEM的數(shù)字地形分析研究主要是集中在地形因子提取與定量分析以及DEM信息容量等方面。如湯國(guó)安、趙牡丹等[3-4]系統(tǒng)研究了坡度分級(jí)方法并對(duì)陜北不同地貌區(qū)進(jìn)行坡度提取、分級(jí),最終對(duì)不同地面坡譜進(jìn)行了統(tǒng)計(jì)分析,劉學(xué)軍等[5-9]將DEM誤差分為DEM模型誤差以及空間要素提取算法誤差,并在考慮DEM誤差自相關(guān)的前提下對(duì)基于DEM 的空間因素(包括坡度、坡向、河網(wǎng)等)的不同提取算法進(jìn)行了系統(tǒng)的誤差分析,陳楠等對(duì)不同尺度DEM做了信息容量比較,定量化研究了信息損失量的規(guī)律性[10-11]。但以上諸多研究多集中在陜北丘陵溝壑區(qū),對(duì)高塬溝壑區(qū)研究較少,尤其是對(duì)位于隴東,有天下第一塬號(hào)稱的董志塬相關(guān)研究更是鮮有報(bào)道。高塬溝壑區(qū)地貌復(fù)雜,塬高、溝深、坡陡,地勢(shì)起伏較大,是水土流失嚴(yán)重地區(qū)之一,其中水蝕和重力侵蝕強(qiáng)烈,所以如何正確快速掌握地形特征規(guī)律,為當(dāng)?shù)叵嚓P(guān)部門尤其是水保部門提供決策依據(jù)是亟待解決的問題。

基于以上情況,本研究選擇典型高塬溝壑區(qū)——硯瓦川流域?yàn)檠芯繀^(qū)域,在眾多空間因子中選取了3個(gè)最能突出反映該地區(qū)地形特征的指標(biāo)因子,指標(biāo)選取基本原則就是這些地形因子能從不同角度反映該區(qū)域地形特征,獨(dú)立性強(qiáng),最終確定坡度、坡長(zhǎng)和剖面曲率三個(gè)地形因子作為該區(qū)域地形特征的主要指標(biāo),因?yàn)槠露却笮≈苯佑绊懼乇砦镔|(zhì)的流動(dòng)及能量轉(zhuǎn)換的規(guī)模與強(qiáng)度,坡長(zhǎng)能夠影響坡面徑流的流速和流量,進(jìn)而影響土壤侵蝕強(qiáng)度,而剖面曲率能夠反映出坡度的變化率。最后結(jié)合GIS,以數(shù)字地形分析技術(shù)為手段,建立相應(yīng)地形因子的數(shù)字地面模型(Digital Terrain M odel簡(jiǎn)稱DTM),并統(tǒng)計(jì)分析各因子分布特點(diǎn)與規(guī)律,望給當(dāng)?shù)亟?jīng)濟(jì)建設(shè)提供支持,為后續(xù)相關(guān)研究提供參考。

1 材料與方法

1.1 研究區(qū)域概況

硯瓦川流域系涇河水系,位于甘肅省慶陽市西峰區(qū)境內(nèi) ,107°37′-107°55′E,35°31′-35°44′N,流域中軸線呈NW-ES走向,總面積約368 km2,該區(qū)屬半濕潤(rùn)季風(fēng)氣候區(qū),年均溫為8.1℃,年均降水量為579 mm,最多年降水量為852.8 mm(1964年),最少年降水量336.7 mm(1986年),7-9月約占全年降水量的60%左右,并多以暴雨形式出現(xiàn),大量降水以地表經(jīng)流的形式匯入江河,發(fā)生暴洪的機(jī)率大。流域以黑壚土為主,占流域面積的70%,主要分布在塬面,土層厚約250 cm,其上松下粘,耕性良好,適耕期長(zhǎng)。流域地面海拔高度為 940～1 409m,平均高程 1 243 m,平均溝壑密度2.17 km/km2,溝谷是侵蝕最劇烈的部位,重力侵蝕頻繁發(fā)生,對(duì)區(qū)域經(jīng)濟(jì)發(fā)展、居住安全以及下游水利工程正常運(yùn)行帶來不利影響,水土保持任務(wù)艱巨。

1.2 數(shù)據(jù)基礎(chǔ)及預(yù)處理

目前制作DEM的主要方法有兩種,即通過遙感立體像得到DEM,或是基于數(shù)字地形圖通過內(nèi)插的方法得到DEM,后者是迄今為止最為成熟和經(jīng)濟(jì)實(shí)用的DEM建立方法,本文綜合考慮該區(qū)域的數(shù)據(jù)可獲性,選擇后者建立本研究區(qū)域的DEM,所使用到的基礎(chǔ)數(shù)據(jù)為覆蓋研究區(qū)域的4幅1∶5萬國(guó)家標(biāo)準(zhǔn)地形圖,基本參數(shù)為:54北京坐標(biāo)系,56黃海高程系,等高距為20 m,使用投影為高斯-克呂格投影。

一般來說,地形圖掃描后都有一定變形,必須進(jìn)行預(yù)處理,也就是掃描糾正,并對(duì)糾正后的柵格圖進(jìn)行檢查,才能得到矢量化準(zhǔn)確的數(shù)字等高線。本研究采用A rcGIS 9.2提供的Georeferencing工具進(jìn)行圖像糾正,先在每幅地形圖上選取12個(gè)均勻分布的控制點(diǎn),為方便計(jì)算,本研究選取的都是格網(wǎng)交點(diǎn),然后建立鏈接表,并逐個(gè)對(duì)控制點(diǎn)進(jìn)行誤差檢查,合格后對(duì)圖像進(jìn)行重采樣,最后得到經(jīng)過精確糾正的掃描圖。接下來使用Geoway 3.5對(duì)該掃描圖進(jìn)行等高線、高程點(diǎn)、河流、湖泊以及流域邊界線的矢量化,通過要素拼接最終得到硯瓦川流域完整的數(shù)字地形圖。

1.3 建立水文增強(qiáng)型DEM

DEM是一組有序數(shù)據(jù)陣列對(duì)地表高度空間變化的定量表現(xiàn)。是一定范圍內(nèi)規(guī)則格網(wǎng)點(diǎn)的平面坐標(biāo)(X,Y)及其高程(Z)的數(shù)據(jù)集,它主要是描述區(qū)域地貌形態(tài)的空間分布,可通過等高線或相似立體模型進(jìn)行數(shù)據(jù)采集(包括采樣和量測(cè)),然后進(jìn)行數(shù)據(jù)內(nèi)插而形成的。根據(jù)DEM對(duì)水文特征和地貌關(guān)系的表現(xiàn)能力,可將DEM分為2類,即水文增強(qiáng)型DEM(hydrologically enhanced DEM s/drainage enforced DEM s)[12],也稱為水文地貌關(guān)系正確DEM(hydrologically correct DEM s,H c-DEM)[13],以及非水文增強(qiáng)型DEM。本研究采用澳大利亞國(guó)立大學(xué)Hutchinson教授開發(fā)的ANUDEM軟件[13-14]來建立水文增強(qiáng)型DEM(H c-DEM),考慮到本研究區(qū)域?yàn)檐珔^(qū),相對(duì)于丘陵溝壑區(qū)地形平坦,用ANUDEM軟件插值時(shí)設(shè)置的三個(gè)關(guān)鍵參數(shù)[15]分別為:分辨率10m,第二糙度系數(shù)0.5,迭代次數(shù)35,最終得到水文增強(qiáng)的DEM:H c-DEM,圖1是基于H c-DEM建立的本研究區(qū)域的三維示意圖。

1.4 地形指標(biāo)計(jì)算

坡度是指地球表面任一點(diǎn)通過該點(diǎn)的切平面與水平地面的夾角,表示地表面在該點(diǎn)的傾斜程度,是一個(gè)重要的地形因子。坡度大小直接影響著地表物質(zhì)的流動(dòng)及能量轉(zhuǎn)換的規(guī)模與強(qiáng)度。坡度的計(jì)算如公式(1),其實(shí)現(xiàn)過程是在3×3的DEM 柵格窗口中進(jìn)行的,窗口在DEM數(shù)據(jù)矩陣中連續(xù)移動(dòng)后完成整幅圖的計(jì)算工作。

式中:fx,fy——南北和東西方向上的偏導(dǎo)數(shù)。目前對(duì) fx,fy的計(jì)算有二階差分、三階不帶權(quán)差分(Sharp nack算法)、三階反距離平方權(quán)差分(Horn算法)、三階反距離權(quán)差分、Frame差分、簡(jiǎn)單差分等6種常見算法,對(duì)于這6種算法的分析比較,國(guó)內(nèi)學(xué)者已經(jīng)做了大量的研究[5,16-17],最終,本研究確定使用三階反距離平方權(quán)差分(Horn算法),該算法已經(jīng)集成到A rcGIS中,利用該算法提取的坡度圖如附圖4。在實(shí)際應(yīng)用中,坡度有兩種表達(dá)方法,一種是坡度即水平面與地面表面之間的夾角;另一種是坡度百分比即高程增加量與水平增加量之比的百分?jǐn)?shù),本研究采用第一種方法來表達(dá)坡度。剖面曲率是對(duì)地面坡度沿最大坡降方向地面高程變化率的度量,通過在A rcGIS下對(duì)坡度求“slope”來實(shí)現(xiàn)提取(附圖5)。

坡長(zhǎng)通常是指在地面上沿水流方向到其流向起點(diǎn)間的最大地面距離在水平面上的投影長(zhǎng)度。坡長(zhǎng)因子是水土保持上的重要因子之一,坡長(zhǎng)通過影響坡面徑流的流速和流量,影響水流挾沙力,進(jìn)而影響土壤侵蝕強(qiáng)度,因而是定量計(jì)算土壤流失的重要指標(biāo),也是土壤侵蝕模型中的必要參數(shù)。在不考慮其他影響因素時(shí),坡長(zhǎng)越大,徑流速度就愈大,匯集的徑流也愈大,水力侵蝕也愈強(qiáng)[18]。研究坡長(zhǎng)對(duì)產(chǎn)流產(chǎn)沙過程對(duì)建立土壤侵蝕預(yù)報(bào)模型、制訂有效的水土保持體系有重要的理論價(jià)值。本研究對(duì)坡長(zhǎng)的提取是基于Hickey和Van Remortel等提出的基于DEM提取區(qū)域坡度坡長(zhǎng)因子的算法,H ickey和Van Rem ortel等編制了相應(yīng)AM L語言程序[19-21],后又改進(jìn)為利用C++語言編寫的程序,使運(yùn)行效率大為提高[22]。本研究正是通過運(yùn)行改進(jìn)后的程序來實(shí)現(xiàn)坡長(zhǎng)的提取(附圖6),其流程圖如圖2。

圖2 硯瓦川流域坡長(zhǎng)提取流程圖

2 結(jié)果與分析

2.1 流域坡度特征分析

將坡度圖屬性表導(dǎo)出,利用Visual Foxp ro進(jìn)行頻率統(tǒng)計(jì)得到圖6,可知流域的坡度分布主要集中在1°～5°,其在累積頻率曲線上表現(xiàn)為曲線增勢(shì)較陡,當(dāng)坡度＞5°時(shí),累積頻率曲線的增勢(shì)趨于緩和,且5°以上的各坡度頻率分布均小于3%。流域16%的面積坡度在 1°以下,坡度在9°以下的約占流域面積的50%,＞30°的地區(qū)不到流域總面積的20%。流域坡度分布在 0～78.5°,平均坡度為14.9°,坡度標(biāo)準(zhǔn)差為 14.6,變異系數(shù)為 0.98,說明研究區(qū)的地面較為平坦,溝壑發(fā)育還不完全。

圖3 硯瓦川流域坡度頻率曲線圖

根據(jù)文獻(xiàn)資料記錄,在黃土高原丘陵溝壑區(qū)縣南溝流域[23],基于5m分辨率DEM提取的坡度最小值為0°,最大值為 72.9°,平均坡度為28.9°;榆林地區(qū)典型樣區(qū)[24]10 m分辨率DEM提取的坡度最小值為0°,最大值為 59.8°,平均坡度為20.3°,坡度標(biāo)準(zhǔn)差為10.48,變異系數(shù)為0.52;延川樣區(qū)的平均坡度為 32.27°,坡度標(biāo)準(zhǔn)差為 11.94,變異系數(shù)為0.37。比較結(jié)果知,相對(duì)于丘陵溝壑區(qū),本研究區(qū)的坡度分布范圍較廣,平均坡度小,說明區(qū)域內(nèi)平緩坡面積比較大,而坡度標(biāo)準(zhǔn)差和變異系數(shù)均大于丘陵溝壑區(qū),反映了研究區(qū)域的地形起伏波動(dòng)程度要大于丘陵溝壑區(qū),具體統(tǒng)計(jì)數(shù)據(jù)見表1。

表1 研究區(qū)及丘陵溝壑區(qū)若干樣區(qū)地形特征值統(tǒng)計(jì)表[23-24]

2.2 流域剖面曲率的特征分析

導(dǎo)出流域的剖面曲率圖屬性表,并在V isual Foxp ro下進(jìn)行頻率統(tǒng)計(jì)(圖4)。本研究區(qū)的剖面曲率分布主要集中在1～3,其中曲率為1的地區(qū)面積最廣,約占流域的29.2%,其在累積頻率曲線上的表現(xiàn)為曲線增勢(shì)很陡,當(dāng)曲率＞1時(shí),累積頻率曲線的增勢(shì)趨于緩和,4以上各曲率的頻率分布均小于3%;由累積頻率曲線知,流域50%的區(qū)域坡面曲率小于7,＞24的地區(qū)不到流域總面積的20%。流域剖面曲率分布在0～75.4,平均剖面曲率為12,曲率標(biāo)準(zhǔn)差為13.6,變異系數(shù)為1.13。

圖4 硯瓦川流域剖面曲率頻率曲線圖

根據(jù)文獻(xiàn)資料記錄[24],在黃土高原丘陵溝壑區(qū),榆林地區(qū)典型樣區(qū)10 m分辨率DEM提取的剖面曲率最小值為0,最大值為67.8,平均剖面曲率為21.6,曲率標(biāo)準(zhǔn)差為13.2,變異系數(shù)為0.61。比較結(jié)果知,相對(duì)于丘陵溝壑區(qū),本研究區(qū)的剖面曲率分布稍廣,而平均剖面曲率不到丘陵溝壑區(qū)一半,說明陡坡區(qū)域相比丘陵溝壑區(qū)較少,曲率標(biāo)準(zhǔn)差非常接近,但變異系數(shù)是丘陵溝壑區(qū)的一倍,表明高塬溝壑區(qū)的坡度變化率的離散程度比較大,地貌起伏的波動(dòng)程度要大于丘陵溝壑區(qū)。

2.3 流域坡長(zhǎng)特征分析

通過導(dǎo)出坡長(zhǎng)圖的屬性表,在Visual Foxp ro下進(jìn)行頻率統(tǒng)計(jì)(圖5),本研究區(qū)的坡長(zhǎng)分布主要集中在10～150 m,其在累積頻率曲線上表現(xiàn)為曲線增勢(shì)較陡,其中坡長(zhǎng)為10 m的坡分布最廣,占流域的11.21%,150m以上各坡長(zhǎng)的頻率分布均小于2%;由累積頻率曲線知,流域50%的區(qū)域坡長(zhǎng)小于80m,當(dāng)坡長(zhǎng)大于200m時(shí),累積頻率曲線的趨勢(shì)變緩,坡長(zhǎng)大于240 m的坡只占流域總面積的20%。坡長(zhǎng)較大值均分布在海拔1 300 m以上的高塬區(qū),主要在流域的邊界處,溝壑發(fā)育不完全;坡長(zhǎng)較小值主要分布在溝谷區(qū),這些地區(qū)地表破碎。流域坡長(zhǎng)分布在0～9 583.1m,平均坡長(zhǎng)為181.8m,標(biāo)準(zhǔn)差為350.8,變異系數(shù)為1.93,說明研究區(qū)的地面較為平坦,溝壑發(fā)育還不完全。

圖5 硯瓦川流域坡長(zhǎng)頻率曲線圖

根據(jù)黃土高原丘陵溝壑區(qū)文獻(xiàn)資料記錄,在縣南溝流域[23],基于5 m分辨率DEM提取的坡長(zhǎng)最小值為2.5 m,最大值為485.2m,平均坡長(zhǎng)為33.7m;基于25 m分辨率DEM提取的坡長(zhǎng)最小值為12.8 m,最大值為443.8 m,平均坡長(zhǎng)為54.2 m;在延河流域10m分辨率的典型樣區(qū)[24],坡長(zhǎng)最小值為5m,最大值為1 642.25m,平均坡長(zhǎng)為64.16m,坡長(zhǎng)標(biāo)準(zhǔn)差為65.30,變異系數(shù)為1.02;榆林樣區(qū)基于10m分辨率DEM提取的坡長(zhǎng),最小值為5 m,最大值為571.3 m,平均坡長(zhǎng)為37.05 m,坡長(zhǎng)標(biāo)準(zhǔn)差為38.2,變異系數(shù)為1.03。比較結(jié)果知,相對(duì)于丘陵溝壑區(qū),本研究區(qū)的坡長(zhǎng)分布范圍非常廣,平均坡長(zhǎng)是丘陵溝壑區(qū)的3～6倍,坡長(zhǎng)標(biāo)準(zhǔn)差也是丘陵溝壑區(qū)的5～9倍,變異系數(shù)約是丘陵溝壑區(qū)的2倍,反映了本研究區(qū)地表在水平方向上的平均轉(zhuǎn)折變化水平較小,坡長(zhǎng)平均值較大,說明研究區(qū)的連續(xù)平坦區(qū)域面積較大,這也與本區(qū)域塬面占主導(dǎo)的實(shí)際地形特征相一致。

3 結(jié)論

(1)流域的坡度分布主要集中在1°～5°,且5°以上的各坡度的頻率分布均小于3%。

(2)剖面曲率分布主要集中在1～3,其中曲率1約占流域的29.2%,4以上各曲率的頻率分布均小于3%。

(3)坡長(zhǎng)分布主要集中在10～150 m,10 m的坡長(zhǎng)坡占流域的11.21%,150 m以上各坡長(zhǎng)的頻率分布均小于2%。

(4)相對(duì)丘陵溝壑區(qū),高塬溝壑區(qū)的坡度分布范圍較廣,平均坡度小,區(qū)域內(nèi)平緩坡面積比較大,但坡度變化率的離散程度比較大,地貌起伏的波動(dòng)程度要大于丘陵溝壑區(qū)。

這些結(jié)果的獲取有助于區(qū)域水保措施部門因地制宜,制定出合理的水保措施,為流域規(guī)劃方案的制訂提供理論依據(jù)和技術(shù)指導(dǎo),本研究?jī)H考慮了坡度、剖面曲率與坡長(zhǎng)因子,并未考慮其他地形因子如坡向、高程極差、高程曲線積分等,需要今后進(jìn)一步進(jìn)行研究。

[1] 朱顯謨.中國(guó)黃土高原土地資源[M].西安:陜西科學(xué)技術(shù)出版社,1986.

[2] 李志林,朱慶.數(shù)字高程模型[M].武漢:武漢大學(xué)出版社,2001.

[3] 湯國(guó)安,宋佳.基于DEM坡度圖制圖中坡度分級(jí)方法的比較研究[J].水土保持學(xué)報(bào),2006,20(2):157-160.

[4] 趙牡丹,湯國(guó)安,陳正江,等.黃土丘陵溝壑區(qū)不同坡度分級(jí)系統(tǒng)及地面坡譜對(duì)比[J].水土保持通報(bào),2002,22(4):33-36.

[5] 劉學(xué)軍,任志峰,王彥芳,等.基于DEM的任意方向坡度計(jì)算方法[J].地域研究與開發(fā),2009,28(4):139-141.

[6] 劉學(xué)軍,晉蓓,王彥芳.DEM流徑算法的相似性分析[J].地理研究,2008,27(6):1347-1357.

[7] 劉學(xué)軍,卞璐,盧華興,等.顧及DEM誤差自相關(guān)的坡度計(jì)算模型精度分析[J].測(cè)繪學(xué)報(bào),2008,37(2):200-206.

[8] 劉學(xué)軍,王永君,龔健雅,等.DEM流域網(wǎng)絡(luò)提取算法的誤差特性分析[J].測(cè)繪學(xué)報(bào),2007,36(2):224-230.

[9] 劉學(xué)軍,龔健雅,周啟鳴,等.基于 DEM 坡度坡向算法精度的分析研究[J].測(cè)繪學(xué)報(bào),2004,33(3):258-263.

[10] 陳楠,王欽敏.基于地形因子信息量的數(shù)字高程模型分辨率的選擇:以黃土高原的研究為例[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2009,34(6):692-695.

[11] 陳楠,林宗堅(jiān),湯國(guó)安,等.數(shù)字高程模型的空間信息不確定性分析[J].測(cè)繪通報(bào),2005(11):14-17.

[12] 楊勤科,M cVicar T R,VanN iel T G,等.ANUDEM和TIN兩種建立DEM方法的對(duì)比研究[J].水土保持通報(bào),2006,26(6):84-88.

[13] Hutchinson M F.A new p rocedure for gridding elevation and stream line data with automatic removal of spurious pits[J].Journal of Hyd rology,1989,106(3/4):211-232.

[14] 周買春,黎子浩.數(shù)值地形圖的生成及其水文地貌特征評(píng)價(jià)[J].水利學(xué)報(bào),2002(2):71-74.

[15] 楊勤科,M cVicar T R,VanNiel T G,等.用 ANUDEM建立水文地貌關(guān)系正確DEM的方法研究[J].測(cè)繪科學(xué),2006,31(6):155-157.

[16] 陳楠,王欽敏,湯國(guó)安,等.6種坡度提取算法的應(yīng)用范圍分析:以在黃土丘陵溝壑區(qū)的研究為例[J].測(cè)繪信息與工程,2006,31(4):20-22.

[17] 劉學(xué)軍,王葉飛,曹志東,等.基于 DEM 的坡度坡向誤差空間分布特征研究[J].測(cè)繪通報(bào),2004(12):11-13.

[18] 唐克麗,陳永宗.黃土高原地區(qū)土壤侵蝕區(qū)域特征及綜合治理途徑[M].北京:中國(guó)科學(xué)技術(shù)出版社,1990.

[19] Van Remortel R,H am ilton M,Hickey R.Estimating the LS fac tor for RUSLE through iterative slope length p rocessing of digital elevation data[J].Cartography,2001,30(1):27-35.

[20] H ickey R.Slope angle and slope length solutions for GIS[J].Cartography,2000,29(1):1-8.

[21] H ickey R,Sm ith A,Jankow ski P.Slope Length Calcu lations from a DEM within ARC/INFO GRID[J].Computers,Environment and U rban Systems,1994,18(5):365-380.

[22] Van Remortel H M,H ickey R.Calcu lations LS Factor Using DEM dataaccording to RUSLE-based C riteria.www.on linegeographer.com/slope/rusleV4pc.am l2004.

[23] 謝紅霞.延河流域土壤侵蝕時(shí)空變化及水土保持環(huán)境效應(yīng)評(píng)價(jià)研究[D].西安:陜西師范大學(xué),2008.

[24] 張婷.基于DEM的陜北黃土高原多地形因子空間關(guān)聯(lián)特征研究[D].西安:西北大學(xué),2005.