基于匯流累積計(jì)算的溝壑密度分析方法

吳秉校, 侯 雷, 宋敏敏, 吳發(fā)啟

(西北農(nóng)林科技大學(xué) 資源環(huán)境學(xué)院, 陜西 楊凌 712100)

基于匯流累積計(jì)算的溝壑密度分析方法

吳秉校, 侯 雷, 宋敏敏, 吳發(fā)啟

(西北農(nóng)林科技大學(xué) 資源環(huán)境學(xué)院, 陜西 楊凌 712100)

溝壑密度是評價地表侵蝕影響，水土流失情況，進(jìn)行地貌類型分析等的重要指標(biāo)。文章基于研究區(qū)填洼后的DEM，利用地表匯流累積分析原理和柵格重分類方法，通過研究溝道覆蓋區(qū)柵格數(shù)與匯流累積柵格二值化閾值間的關(guān)系，選擇合適的閾值用于提取地表溝道。然后結(jié)合自然地表高程，分析合理的高程和相對高程柵格重分類閾值并進(jìn)行高程柵格二值化，將二值化結(jié)果作為乘積模板并配合柵格細(xì)化算法優(yōu)化溝道提取結(jié)果。最終利用1 km2單位面積規(guī)則格網(wǎng)對溝道提取結(jié)果進(jìn)行目標(biāo)柵格數(shù)統(tǒng)計(jì)計(jì)算，得到地表溝壑密度值。結(jié)果顯示，該方法基于30 m分辨率DEM提取的西安市地表溝壑密度均值為1.29 km/km2。除雁塔區(qū)，其余各行政區(qū)統(tǒng)計(jì)均值與前人研究結(jié)果基本吻合。試驗(yàn)結(jié)果證明，基于本方法可以高效較精確地計(jì)算該市地表溝壑密度，分析結(jié)果可以量化研究區(qū)溝壑密度特征值以及通過柵格像元灰度變化體現(xiàn)溝壑密度變化趨勢。

溝壑密度; 匯流累積量; 柵格重分類; 柵格細(xì)化

溝壑密度(gully density)是地表每1 km2內(nèi)侵蝕溝(或水文網(wǎng))的總長度，是表征內(nèi)外營力對地表侵蝕影響的重要指標(biāo)，對研究區(qū)域水土流失狀況和等級、土壤侵蝕狀況、地貌演變與地質(zhì)災(zāi)害分析、水土保持區(qū)劃等具有重要意義。溝壑密度的大小與地表的切割破碎程度呈正相關(guān)關(guān)系，對于區(qū)域氣候、土壤、地質(zhì)條件、植被、生態(tài)適宜性研究亦有重要參考價值。作為一個重要的地貌特征表征指標(biāo)，前人在進(jìn)行相關(guān)研究的過程中利用了多種手段和方法進(jìn)行分析。在機(jī)助分析技術(shù)應(yīng)用以前，傳統(tǒng)的溝壑密度量測與計(jì)算主要基于外業(yè)調(diào)查和地形圖。隨著遙感影像和數(shù)字地形處理與分析技術(shù)的成熟和推廣，憑借其快捷、高效、高精度、適用范圍廣的優(yōu)勢逐漸取代了傳統(tǒng)的人為量測方法，成為溝壑密度分析的主要技術(shù)手段。

現(xiàn)在應(yīng)用最廣的溝壑密度分析方法主要基于數(shù)字等高線圖和數(shù)字高程模型(Digital Elevation Model，DEM)來進(jìn)行。利用數(shù)字等高線圖分析是基于矢量圖形分析原理，依據(jù)地形特征，運(yùn)用幾何形態(tài)學(xué)理論，分析矢量線劃的拓?fù)潢P(guān)系和幾何特征，從而進(jìn)行地表溝道提取，完成溝壑密度分析[1-2]。該類方法能夠自動化識別地形的特征點(diǎn)、線，提取結(jié)果連續(xù)性較好，適合于小范圍高精度溝谷提取和密度計(jì)算，但完全依據(jù)幾何理論分析，分析結(jié)果與實(shí)際地形特征有一定偏差。對于DEM分析提取的等高線，由于本身存在提取誤差和鋸齒狀畸形，該方法適用性不佳。利用DEM分析溝壑密度，前人主要依據(jù)兩種方法進(jìn)行：基于水文分析中的匯流累積算法提取出地表徑流路徑作為溝谷線進(jìn)行分析以及基于地表曲率的分析提取。基于地表曲率的方法是通過DEM進(jìn)行地表曲率計(jì)算和求導(dǎo)，通過導(dǎo)數(shù)變化趨勢，分析地形特征點(diǎn)、線，從而提取溝谷[3]。該方法可以避免水文分析中由于匯流閾值給溝道提取結(jié)果帶來的誤差，但處理結(jié)果受DEM自身分辨率和誤差影響較大，且需要人為剔除非溝谷線，溝谷提取結(jié)果連續(xù)性不佳，不利于地形破碎區(qū)的溝壑密度計(jì)算?；趨R流累積計(jì)算的方法是依據(jù)地表水文分析原理，通過一定的河網(wǎng)提取算法，利用DEM計(jì)算流向，匯流累積量并設(shè)置閾值提取地表徑流路徑，即區(qū)域水文網(wǎng)分布。該方法提取的溝谷特征線連續(xù)性較好，分析效率較高[4-7]。但該類方法在河網(wǎng)提取過程中由于DEM自身誤差和算法帶來誤差，造成了偽溝道和低洼平坦地區(qū)不合理平行水系的出現(xiàn)，在一定程度上影響了提取結(jié)果。單純憑借匯流累積量提取溝道沒有充分考慮到區(qū)域的自然地理背景，所以對于分析地形復(fù)雜區(qū)域溝壑密度存在缺陷。且傳統(tǒng)的溝壑密度分析方法多考慮的是大范圍內(nèi)平均溝壑密度，對區(qū)域內(nèi)部溝壑密度的變化狀況不能實(shí)現(xiàn)統(tǒng)計(jì)單元尺度的較精確分析。本文利用DEM，借助于水文分析算法，綜合考慮地形低洼區(qū)域的剔除和自然地理?xiàng)l件因素，優(yōu)化原有的溝壑密度分析算法，實(shí)現(xiàn)較大范圍的地表溝壑密度基于統(tǒng)計(jì)單元尺度的快速較精確分析，并以西安市為研究區(qū)，對該市的30 m分辨率DEM進(jìn)行分析計(jì)算，驗(yàn)證提取結(jié)果的合理性。

1 分析原理與方法

1.1 分析原理

(1) 溝壑密度計(jì)算。溝壑密度也稱為溝谷密度，是評價單位面積內(nèi)侵蝕溝(水文網(wǎng))長度的定量指標(biāo)，一般采用km/km2為單位進(jìn)行度量。自然地表形態(tài)分析中，溝壑可以視為地表徑流路徑。因此溝壑密度的分析可以借助于空間分析方法中的水文分析原理，提取出地表徑流路徑，計(jì)算徑流路徑長度，即可視為侵蝕溝長。

基于規(guī)則格網(wǎng)DEM的溝壑密度計(jì)算在提取溝壑長度時，結(jié)果為柵格數(shù)據(jù)格式，因此運(yùn)用矢量線劃度量必然存在矢柵轉(zhuǎn)化誤差，為了使計(jì)算效率提高，故可以對這一指標(biāo)做簡單轉(zhuǎn)化，以便于進(jìn)行柵格數(shù)據(jù)分析。以面積為1 km2的單元格網(wǎng)進(jìn)行統(tǒng)計(jì)，可以較精確地表現(xiàn)一區(qū)域內(nèi)溝壑密度的變化趨勢，而不僅僅是區(qū)域的溝壑密度均值。如DEM有固定的空間分辨率(CS，Cellsize)，單位為m，以單元面積為1 km2格網(wǎng)作為統(tǒng)計(jì)單元，N為水文網(wǎng)覆蓋區(qū)的目標(biāo)柵格像元數(shù)，則每一統(tǒng)計(jì)單元的溝壑密度(GD)可以通過下式進(jìn)行計(jì)算：

(1)

(2) 溝壑提取。基于水文分析方法的地表溝壑提取與地表徑流路徑提取類似，這里進(jìn)行概述。原始DEM存在誤差和不合理凹陷，自然水流在形成匯流路徑過程中，會填充洼地，并選擇洼地的最低出流點(diǎn)流出。因此在分析前需進(jìn)行地形填洼處理，形成一個平坦的地形，以減少水流路徑分析誤差。

水流路徑分析的關(guān)鍵在于確定在DEM表面，水流流經(jīng)哪些柵格單元，即每個柵格單元的流出方向,這需要在流向分析的基礎(chǔ)上進(jìn)行。在常規(guī)的流向計(jì)算中，常用“八方向法”進(jìn)行度量?！鞍朔较颉彼惴ㄔ谝欢ǖ南裨直媛蕳l件下利用“3×3”模板進(jìn)行逐像元掃描計(jì)算，以中心像元為基點(diǎn)構(gòu)建等分辨率寬度向八方向擴(kuò)展的高程擬合平面，并利用中心像元8方向的高程變化率構(gòu)建算法，分析鄰域范圍內(nèi)值變化率最大的下坡方向[8]。當(dāng)一像元的鄰域存在多個相同的下坡方向值時，會擴(kuò)大搜索，直到找到更大值。具體計(jì)算式如下：

MD=Zmax/distance×100

(2)

其中Zmax為中心像元與周圍8鄰域像元的高程差最大值，distance為相鄰像元間距離。一區(qū)域的流向分析完成后，需要進(jìn)一步分析地表匯流特征。匯流即從地表水流路徑源開始，每經(jīng)過一個像元，匯流累積量增加1，計(jì)算每個下坡像元的流入累積量，即對于每個像元，有多少的像元是其流入像元。最終匯流的分析結(jié)果可以理解為對于一個目標(biāo)像元，水流流經(jīng)時攜帶的上游像元數(shù)。在分析結(jié)果中，匯流值較大的像元必然處于坡降較大，地勢較低處，即為溝壑覆蓋的像元。

(3) 匯流路徑細(xì)化。從式(1)中可以得出，最終匯流提取結(jié)果包含的目標(biāo)柵格數(shù)將直接影響溝壑密度分析結(jié)果。分析上述匯流提取結(jié)果，在地勢起伏較大區(qū)域，匯流路徑的提取結(jié)果多為單像元寬度，在地勢低洼區(qū)，匯流柵格會存在團(tuán)塊狀集聚，增大了目標(biāo)柵格數(shù)和該區(qū)的溝壑密度，引起分析誤差。因此需要對提取結(jié)果進(jìn)行再處理才可以參與最終的溝壑密度計(jì)算，特別是團(tuán)塊區(qū)。這里使用改進(jìn)的柵格細(xì)化算法進(jìn)行處理，在基本不破壞柵格拓?fù)溥B通性的前提下，采用柵格剝離的方法刪除邊緣點(diǎn)，提取中心單像元寬度的特征骨架。假設(shè)該處理區(qū)為多像元寬度提取區(qū)，P為中心像元，對于該像元進(jìn)行判斷分析時，采用“3×3”模板，當(dāng)計(jì)算得到P1+P3+P5+P7=4，P2+P4+P6+P8=4(或0)，說明P被目標(biāo)像元包圍，則該像元保留。當(dāng)P3+P7=0，P2+P3+P4>0，P8+P7+P6>0時，說明P點(diǎn)為線狀柵格的轉(zhuǎn)折點(diǎn)或延續(xù)點(diǎn)，該點(diǎn)保留。當(dāng)P1+P3=0，P2=1，P4+P5+P6+P7+P8>0時，說明P點(diǎn)為線交叉點(diǎn)，進(jìn)行保留[9-10]。通過迭代式的剝離去除，對保留點(diǎn)進(jìn)行賦值區(qū)別處理后累計(jì)分析，最終保證提取結(jié)果低洼區(qū)團(tuán)塊狀集聚的柵格細(xì)化為單像元寬度的較連續(xù)線狀柵格，提升分析的合理性。

1.2 分析方法

原始DEM存在誤差和不合理凹陷，故在進(jìn)行坡向域計(jì)算之前，必須對原始DEM進(jìn)行填洼處理。利用上述方法對填洼后DEM進(jìn)行流向計(jì)算，并計(jì)算匯流累積量。對匯流累積計(jì)算結(jié)果進(jìn)行柵格細(xì)化，保證每條水流路徑為單像元寬度。將處理結(jié)果進(jìn)行像元二分處理，結(jié)合區(qū)域高程特征和DEM分辨率，選擇合理閾值，將小于閾值的匯流累積柵格賦值為“0”，大于閾值的賦值為“1”，“1”值柵格覆蓋區(qū)即為地表徑流流路，則可以得到表征溝壑分布的二值化柵格。

對于地形起伏較大區(qū)域來說，山谷低洼平坦區(qū)和平原區(qū)計(jì)算溝壑密度會影響分析結(jié)果的合理性，因此需要利用相對高程和高程的重分類與柵格像元的邏輯計(jì)算來輔助消除這些不合理區(qū)。高程即填洼后DEM的柵格屬性值，相對高程可以通過柵格鄰域統(tǒng)計(jì)的方法實(shí)現(xiàn)。利用一定的統(tǒng)計(jì)窗口，逐像元計(jì)算統(tǒng)計(jì)中心像元在搜索窗口范圍內(nèi)的高程變化值，其結(jié)果與窗口大小密切相關(guān)。本文中采用“3×3”窗口進(jìn)行分析。對于高程和相對高程分析結(jié)果進(jìn)行柵格重分類，高程重分類結(jié)果主要用于消除平原的不合理提取區(qū)，相對高程重分類結(jié)果主要用于消除地形起伏較小區(qū)和谷地低洼平坦區(qū)的溝壑提取結(jié)果。兩者的重分類是將屬性為一定閾值以上的提取結(jié)果賦值為“1”，閾值以下的賦值為“0”，將重分類的結(jié)果與溝壑二值化柵格分別進(jìn)行柵格相乘運(yùn)算。類似于Intersection運(yùn)算，在相乘過程中，高程和相對高程值低于閾值的溝壑像元值會賦值計(jì)算為“0”，達(dá)到自動去除的目的。

也就能夠得出教師的相關(guān)評價數(shù)值，根據(jù)每個教師所得出的數(shù)值的不同對教師加以判斷和評價。而這一過程在具體計(jì)算的時候相對比較復(fù)雜，這就需要將這些程序能夠和計(jì)算機(jī)相互配合使用，使得數(shù)據(jù)輸入之后，計(jì)算機(jī)能夠在經(jīng)過其中具體程序計(jì)算之后，得出必要的數(shù)據(jù)結(jié)果，根據(jù)相關(guān)結(jié)果內(nèi)容對教師群體加以評價，使得教師們能夠?qū)ψ陨淼脑u價更加公平公正對待，促使教師的綜合素質(zhì)不斷加以提升。

本文中方法計(jì)算溝壑密度是根據(jù)統(tǒng)計(jì)單元中所包含的“1”值像元數(shù)來計(jì)算各單元中特征線所覆蓋面積，再計(jì)算出溝壑密度。由于通過上述的多重重分類二值化，提取出的溝壑區(qū)目標(biāo)柵格屬性值已全部賦值為“1”，因此通過區(qū)域統(tǒng)計(jì)工具，創(chuàng)建一定單位面積的統(tǒng)計(jì)格網(wǎng)，采用單元求和法，計(jì)算每一單元中的特征柵格數(shù)，采用下式計(jì)算即可獲得該統(tǒng)計(jì)單元的溝壑密度。

(3)

其中Gsum為1 km2統(tǒng)計(jì)格網(wǎng)中每一單元所包含的“1”值柵格像元數(shù)。 在匯流累積二值化閾值，高程和相對高程重分類閾值的選擇上，需要根據(jù)實(shí)地的地形地貌條件和DEM分辨率進(jìn)行試驗(yàn)分析。下文將以西安市為試驗(yàn)區(qū)域，利用上述方法，分析這兩者取值對于該區(qū)域溝壑密度計(jì)算結(jié)果的影響。整體的算法流程如下圖所示(圖1)。在ArcGIS 10.2版本軟件中利用Python語言調(diào)用Arcpy函數(shù)包，編寫溝壑密度分析的腳本工具，該工具通過用戶交互輸入DEM柵格數(shù)據(jù)，匯流累積二值化閾值，高程重分類閾值，相對高程重分類閾值，目標(biāo)區(qū)域掩膜文件和統(tǒng)計(jì)格網(wǎng)矢量文件，運(yùn)行程序計(jì)算結(jié)果并保存[11]。統(tǒng)計(jì)格網(wǎng)采用漁網(wǎng)工具(Fishnet)生成覆蓋全域的以1 km2為單元的規(guī)則格網(wǎng)。使用的DEM數(shù)據(jù)為ASTER GDEM V2版，空間分辨率為30 m，投影為UTM/WGS84，經(jīng)拼接裁切處理后能夠滿足西安市全域溝壑密度分析需要。

1.3 試驗(yàn)區(qū)概述

本文以西安市市域范圍作為試驗(yàn)研究區(qū)。西安市位于關(guān)中平原中部，整體地貌形勢呈南高北低，南部為秦嶺山地，北部和西北部分布著大面積渭河及其支流的洪沖積平原，東部區(qū)為驪山低山丘陵，境內(nèi)海拔高度差異懸殊位列全國各城市之冠[12]。秦嶺山地與渭河平原界限分明，構(gòu)成西安市的地貌主體。在該市東南方向淺山與平原交匯區(qū)分布著面積較大的黃土臺塬，西部淺山區(qū)北麓亦有小面積臺塬分布。該市的平原高程多在海拔450 m以下，山地海拔在1 000 m以上居多，臺塬和低山海拔500 ～1 300 m。該市的地貌總體特征主要表現(xiàn)為南部秦嶺山體高大，多深溝峽谷地形，山體被南北向多條河流呈梳狀切割。山前分布有連片的洪(沖)積扇群，呈波浪形向北傾斜。低山與臺塬區(qū)溝壑縱橫，相對高程變化劇烈，地形最為破碎，其中溝壑發(fā)育最密集區(qū)主要是在驪山北麓、南麓臺塬與秦嶺淺山交界區(qū)以及藍(lán)田縣西部臺塬區(qū)。綜上所述，該區(qū)域的地形類型眾多，地貌形勢復(fù)雜多變，中低山、丘陵、臺塬、平原皆有分布，適合于對文中的溝壑密度分析算法合理性進(jìn)行檢測。

1.4 關(guān)鍵閾值分析

本文參考試驗(yàn)區(qū)自然地理背景條件，對于“匯流累積重分類閾值”，“高程重分類閾值”和“相對高程重分類閾值”三者取值方法進(jìn)行分析。在利用地表匯流算法進(jìn)行溝壑密度分析過程中，匯流累積重分類閾值的設(shè)置對于最終分析結(jié)果有巨大影響。在考慮DEM分辨率和區(qū)域地形特征的基礎(chǔ)上，選擇不同的閾值，對提取結(jié)果進(jìn)行溝壑覆蓋柵格數(shù)對比。閾值最小值設(shè)置為100，最大值設(shè)置為1 000，每隔50取值一次，統(tǒng)計(jì)分析結(jié)果如圖2所示。分析圖中統(tǒng)計(jì)曲線和趨勢函數(shù)可得出，當(dāng)閾值大于230時，曲線斜率開始小于趨勢函數(shù)斜率，曲線走勢趨緩，溝壑覆蓋柵格數(shù)受匯流累積重分類閾值影響開始降低[13]。因此，在試驗(yàn)區(qū)域的溝壑密度分析中，匯流累積重分類閾值取230。

圖1 溝壑密度分析方法示意圖

圖2 溝壑覆蓋柵格數(shù)和閾值關(guān)系曲線

綜合考慮DEM分辨率以及區(qū)域地表高程特征，將高程分級閾值設(shè)置為450 m。在計(jì)算過程中，海拔低于450 m的地表將被認(rèn)為不存在溝壑，這主要考慮到該市的平原平坦區(qū)和地形起伏區(qū)的臨界高程。如果研究區(qū)存在大面積地形平坦區(qū)或海拔較低，該值可以近似設(shè)定為0 m，即高程不做重分類處理。在分析中綜合高程重分類結(jié)果和整體地貌特征，將相對高程的重分類閾值設(shè)定在10 m以下進(jìn)行研究。相對高程重分類閾值取0時，說明相對高程不做重分類處理。對相對高程設(shè)置不同閾值后分析結(jié)果例舉如表1所示，表中還相應(yīng)分析了不同取值條件下溝壑密度分級后所各等級占據(jù)的區(qū)域面積比。參考前人的研究成果，西安市的平均溝壑密度為1.29 km/km2[14]。當(dāng)相對高程的重分類閾值取0時，該分析方法計(jì)算出西安市平均溝壑密度為1.35 km/km2。因此，進(jìn)行相對高程重分類二值化是必要的，對于地形平坦區(qū)進(jìn)行分析時，該值設(shè)為0 m。對于本研究區(qū)域，當(dāng)該值取3 m時，計(jì)算結(jié)果與前人研究的區(qū)域溝壑密度均值基本吻合，因此較為合理。將上述重要參數(shù)帶入分析程序中，計(jì)算該區(qū)的溝壑密度。

2 結(jié)果與分析

通過DEM處理分析和分析參數(shù)的設(shè)置研究，將匯流累積重分類閾值設(shè)置為230，高程重分類閾值設(shè)置為450 m，相對高程重分類閾值設(shè)置為3 m。將上述參數(shù)與區(qū)域統(tǒng)計(jì)格網(wǎng)，區(qū)域矢量圖帶入分析程序中，分析得到的溝壑分布情況和溝壑密度分析結(jié)果(圖3)。

由于西安市范圍較大，地形多樣，全域溝壑提取結(jié)果線劃過于復(fù)雜，故按照地貌類型例舉展示低山丘陵區(qū)、秦嶺淺山區(qū)和洪(沖)擊扇區(qū)的局部溝壑提取結(jié)果，見圖4。結(jié)果顯示西安市全域平均溝壑密度為1.29 km/km2，極大值為3.38 km/km2，極小值為0 km/km2。本分析方法提取結(jié)果除雁塔區(qū)，其余各行政區(qū)統(tǒng)計(jì)均值與前人研究結(jié)果基本吻合[14]，而且本方法在溝壑密度極值分析上有突破。對比與雁塔區(qū)具有相近地貌形態(tài)的灞橋區(qū)分析結(jié)果，并參考該行政區(qū)自然地貌統(tǒng)計(jì)數(shù)據(jù)可知，雁塔區(qū)中有近1/3的地貌類型屬于黃土臺塬，塬面破碎、丘洼相間，東部被河流深切，塬畔多沖溝發(fā)育，水土流失嚴(yán)重[15]。故可以推斷，前人基于人為量算，對于塬區(qū)溝壑密度分析結(jié)果可能偏小。參考前人研究的溝壑密度分級評價標(biāo)準(zhǔn)，該市溝壑密度大于2 km/km2的區(qū)域約占總面積的36%，其中溝壑密度<1 km/km2的區(qū)域面積最大，占據(jù)總區(qū)域面積的36.78%，≥3 km/km2的區(qū)域面積最小，約占據(jù)總區(qū)域面積的1.21%。

按照區(qū)域行政單元，將溝壑密度按照“<1 km/km2”，“1～2 km/km2”，“2～3 km/km2”，“≥3 km/km2”分級統(tǒng)計(jì)出西安市各區(qū)縣分布面積比和平均溝壑密度，結(jié)果見表2。通過分析可知，該方法計(jì)算得到西安市溝壑密度均值最大縣區(qū)為藍(lán)田縣，處于秦嶺與平原交界帶的周至縣、戶縣和長安區(qū)溝壑密度也較大。這些區(qū)域主要包含了復(fù)雜的地貌類型，存在大面積的淺山、臺塬區(qū)和洪沖積扇，溝壑縱橫密布。參考自然地理?xiàng)l件，這些地區(qū)的溝壑密度較大，該方法分析結(jié)果合理。臨潼區(qū)和灞橋區(qū)的溝壑密度均值較小，但區(qū)內(nèi)的溝壑密度變化范圍大，這是由于兩區(qū)位于臺塬低山與平原過渡區(qū)，存在大面積平坦地形，溝壑密度均值較小，但臺塬低山區(qū)溝壑密度大，因此極大值較大。由于高程重分類因子的加入，位于平原及城市中心地帶的未央?yún)^(qū)、蓮湖區(qū)、高陵區(qū)等區(qū)縣的溝壑密度參照文中方法分析均為0 km/km2，符合實(shí)地條件，避免了單純的匯流累積分析帶來的結(jié)果歧義，保證了分析合理性。北部閻良區(qū)的小部分地區(qū)位于渭北高原南緣的丘陵溝壑區(qū)，溝壑密度較大，通過文中分析方法亦得到較好的分析體現(xiàn)。

表1 相對高程重分類閾值計(jì)算結(jié)果對比

圖3 西安市溝壑密度計(jì)算結(jié)果

圖4 特征地貌類型區(qū)溝壑提取結(jié)果

縣區(qū)名比例/%<1km/km21～2km/km22～3km/km2≥3km/km2溝壑密度/(km·km-2)最小值最大值平均值前人研究結(jié)果藍(lán)田縣7.6844.9446.730.650.003.231.902.00周至縣24.3936.4838.640.500.003.381.571.49戶縣42.8727.9328.370.820.003.221.281.36長安區(qū)32.9732.6634.020.350.003.171.351.29臨潼區(qū)66.9817.8314.770.430.003.220.720.78灞橋區(qū)65.1919.8314.670.310.003.330.580.46雁塔區(qū)77.1515.956.900.000.002.400.430.06閻良區(qū)99.180.820.000.000.001.470.020.02新城區(qū)100.000.000.000.000.000.000.000.00碑林區(qū)100.000.000.000.000.000.000.000.00高陵區(qū)100.000.000.000.000.000.000.000.00蓮湖區(qū)100.000.000.000.000.000.000.000.00未央?yún)^(qū)100.000.000.000.000.000.000.000.00

3 結(jié) 論

本文利用水文分析原理，基于DEM分析地表匯流累積量，通過對匯流計(jì)算結(jié)果進(jìn)行柵格重分類二值化處理，選擇合理閾值，提取出溝道。并考慮研究區(qū)實(shí)際地貌特征，增加高程和相對高程二值化因子，對溝道提取結(jié)果進(jìn)行優(yōu)化。并利用構(gòu)建1 km2單元統(tǒng)計(jì)格網(wǎng)的方法，逐格統(tǒng)計(jì)目標(biāo)柵格數(shù)并計(jì)算出研究區(qū)的溝壑密度。該方法有如下特征：

(1)能夠?qū)崿F(xiàn)大范圍溝壑密度的快速較精確計(jì)算，同時可以體現(xiàn)區(qū)域地表與統(tǒng)計(jì)格網(wǎng)同尺度的地表溝壑密度變化趨勢，方便地獲取適合于特定分析需要的溝壑密度均值，極值等量值信息。在計(jì)算過程中，不需要對于分析數(shù)據(jù)進(jìn)行過多人為干預(yù)和刪改，實(shí)現(xiàn)效率較高，同時具有客觀性。

(2)在溝壑密度計(jì)算過程中，利用公式的簡單變形，將原先基于矢量圖形的計(jì)算變?yōu)槔酶窬W(wǎng)統(tǒng)計(jì)特征柵格數(shù)量方法，避免了溝道提取結(jié)果的柵矢轉(zhuǎn)化誤差和轉(zhuǎn)化后歧異圖形的出現(xiàn)。

(3)運(yùn)用邏輯判斷原理，基于一定的分析模板，將高程和相對高程進(jìn)行重分類二值化，并對溝道提取結(jié)果進(jìn)行柵格求積處理。通過該步驟可以實(shí)現(xiàn)低海拔平坦區(qū)溝壑密度計(jì)算異常值的剔除，考慮了區(qū)域的實(shí)際地貌特征，提升計(jì)算合理性。

該方法的處理過程涉及到分析窗口大小的選擇，高程和相對高程重分類閾值的分析以及溝道提取柵格結(jié)果的細(xì)化。其中重分類閾值的選擇對于最終溝壑密度分析結(jié)果具有較大影響。如何優(yōu)化這一分析過程仍需進(jìn)一步研究。由于格網(wǎng)是規(guī)則的，而研究區(qū)邊界會被格網(wǎng)切割。因此在利用該方法進(jìn)行溝壑密度分析時還需注意參與分析的DEM數(shù)據(jù)覆蓋范圍應(yīng)當(dāng)適當(dāng)大于統(tǒng)計(jì)格網(wǎng)范圍，如果與目標(biāo)區(qū)等同，則在后期的格網(wǎng)統(tǒng)計(jì)時，邊緣部分所在的格網(wǎng)統(tǒng)計(jì)結(jié)果會明顯小于鄰域結(jié)果，造成分析結(jié)果不合理。參考研究區(qū)試驗(yàn)結(jié)果，該方法在合理的參數(shù)設(shè)置條件下，分析得到的溝壑密度結(jié)果符合區(qū)域?qū)嶋H，可以較為準(zhǔn)確地體現(xiàn)樣區(qū)的溝壑密度變化趨勢，具有高效性和適用性。

[1] 熊漢江,李秀娟.一種提取山脊線和山谷線的新方法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2015,40(4):498-502.

[2] 靳海亮,康建榮,高井祥.利用等高線數(shù)據(jù)提取山脊(谷)線算法研究[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2005,30(9):809-812.

[3] 趙鵬,張立朝,張文詩,等.一種利用地形曲率提取地表溝道算法[J].測繪科學(xué)技術(shù)學(xué)報(bào),2014(1):13-17.

[4] 方磊,劉文兆,李懷有.基于GIS的黃土高塬溝壑區(qū)硯瓦川流域地形特征提取與分析[J].水土保持研究,2010,17(3):7-11.

[5] 曾紅偉,李麗娟,柳玉梅,等. Arc Hydro Tools及多源DEM提取河網(wǎng)與精度分析:以洮兒河流域?yàn)槔齕J]. 地球信息科學(xué)學(xué)報(bào),2011,13(1):22-31.

[6] 宋曉猛,張建云,占車生,等.基于DEM的數(shù)字流域特征提取研究進(jìn)展[J]. 地理科學(xué)進(jìn)展,2013,32(1):31-40.

[7] Turcotte R, Fortin J P, Rousseau A N, et al. Determination of the drainage structure of a watershed using a digital elevation model and a digital river and lake network[J]. Journal of Hydrology, 2001,240(3):225-242.

[8] Smith M J D, Goodchild M F, Longley P. Geospatial analysis: A Comprehensive Guide to Principles, Techniques and Software Tools[M]. Winchelsea, UK:The Winchelsea Press,2015.

[9] 張欣,陳國雄,鐘耳順.優(yōu)化柵格細(xì)化算法的線狀地物提取[J].地球信息科學(xué)學(xué)報(bào),2007,9(3):25-27.

[10] 郭邦梅,王濤,趙榮,等.面狀要素骨架線提取算法的研究[J].測繪通報(bào),2010(12):17-19.

[11] Pimpler E. Programming ArcGIS10.1with Python Cookbook [M]. Birmingham U K: Packt Publishing,2013.

[12] 宋德明,陰治農(nóng),馬正林,等.西安市地理志[M].西安:陜西人民出版社,1988.

[13] 李俊,湯國安,張婷,等.利用DEM提取陜北黃土高原溝谷網(wǎng)絡(luò)的匯流閾值研究[J].水土保持通報(bào),2007,27(2):75-78.

[14] 馬鐘彥.西安市水利志[M].西安:陜西人民出版社,1999.

[15] 西安市雁塔區(qū)地方志編纂委員會.雁塔區(qū)志[M].西安:三秦出版社,2003.

GullyDensityAnalysisBasedonFlowAccumulationMethod

WU Bingxiao, HOU Lei, Song Minmin, WU Faqi

(CollegeofNaturalResourcesandEnvironment,NorthwestA&FUniversity,Yangling,Shaanxi712100,China)

Gully density is an important indicator to assess the impact of surface erosion, soil erosion, landform analysis etc. Based on DEM filling, using surface flow accumulation and raster reclassification analysis principles, we researched the relationship between the number of grids in the channel coverage area and the binarization threshold of the accumulation grids, and selected an appropriate threshold to extract the surface channel. According to the natural surface elevation, the reasonable thresholds were analyzed. Reclasscifying the elevation and relative elevation rasters and thining, we used the binarization results as the product template to optimize the channel extraction results. Finally, we calculated the number of target grids in each 1-km2unit area to qualify the grid raster density of research zone. The results showed that the mean value of surface gully density in Xi′an City is 1.29 km/km2based on DEM with 30 m resolution. In addition to the Yanta District, the results were consistent with previous studies in average gully density of the remaining administrative regions. It certificated that this analytical method can be used to calculate the density of the urban gully effectively and accurately. The results can be used to quantify the characteristic values of the gully density in the study area and reflect the variation trend of the gully density through the gray-scale variation of the grid pixels.

gully density; flow accumulation; raster reclassified; grid refinement

2016-05-02

：2016-07-07

國家自然科學(xué)基金資助項(xiàng)目“黃土坡耕地地表糙度對產(chǎn)匯流的影響研究”(41271288)

吳秉校(1991—),男,安徽六安人,碩士,主要從事GIS空間分析方法與系統(tǒng)設(shè)計(jì)研究。E-mail:18291936701@163.com

吳發(fā)啟(1957—),男,陜西黃陵人,教授,博士生導(dǎo)師，主要從事水土保持與GIS應(yīng)用方面的研究。E-mail:wufaqi@263.net

P931; TP79

：A

：1005-3409(2017)03-0039-06