朱傳華，李化平

(安徽建筑大學(xué) 環(huán)境與能源工程學(xué)院，安徽 合肥 230601)

0 引言

水土流失是我國(guó)主要環(huán)境問(wèn)題之一，水土流失生態(tài)過(guò)程發(fā)生的潛在可能性及其程度的研究，對(duì)合理制定區(qū)域水土保持措施，規(guī)劃區(qū)域可持續(xù)發(fā)展具有重要意義[1]。對(duì)于土壤侵蝕風(fēng)險(xiǎn)評(píng)價(jià)，國(guó)外比較知名的模型分為經(jīng)驗(yàn)?zāi)Ｐ?如RUSLE)與物理模型(SWAT)兩大類(lèi)，目前我國(guó)的水土流失模型大多為借鑒國(guó)外模型適應(yīng)于本土的土壤侵蝕狀況[1，2]。國(guó)內(nèi)土壤侵蝕模型的研究也取得了一定成果[3-5]，但以上模型都需要相當(dāng)多的時(shí)間來(lái)收集模型建立所需的數(shù)據(jù)，如降雨和土壤理化信息等。隨著RS和GIS技術(shù)的廣泛應(yīng)用，參照我國(guó)水利部部頒標(biāo)準(zhǔn)《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》(SL190-2007)，可相對(duì)快速地度量一個(gè)區(qū)域的水土侵蝕強(qiáng)度[6-8]。本文借鑒該思路，在GIS和RS技術(shù)的支持下，綜合考慮自然因素(坡度和植被覆蓋度)和人為因素(土地利用類(lèi)型)，參考《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》，應(yīng)用Python語(yǔ)言編寫(xiě)水土流失風(fēng)險(xiǎn)度分類(lèi)模型，對(duì)董鋪水庫(kù)庫(kù)區(qū)的水土流失風(fēng)險(xiǎn)度進(jìn)行評(píng)價(jià)，并使用轉(zhuǎn)移矩陣和動(dòng)態(tài)度等方法定量分析2004-2018年間4個(gè)年份董鋪水庫(kù)庫(kù)區(qū)水土流失風(fēng)險(xiǎn)度動(dòng)態(tài)變化特征，以期為水庫(kù)水土流失防治和生態(tài)環(huán)境保護(hù)提供參考。

1 數(shù)據(jù)來(lái)源與研究方法

1.1 研究區(qū)概況

合肥市位于北緯31o52'、東經(jīng)117o17'，安徽省省會(huì)城市，2017年末全市常住人口796.5萬(wàn)人。董鋪水庫(kù)和大房郢水庫(kù)位于市區(qū)西北近郊，是合肥市飲用水主要來(lái)源。庫(kù)區(qū)地形西北高東南低，高程介于0-254 m，均屬第四紀(jì)更新世松散堆積層覆蓋，河谷為不對(duì)稱(chēng)河谷，南岸高北岸平緩。土壤以黃棕壤、水稻土兩類(lèi)為主要土壤。研究區(qū)位于亞熱帶之北邊緣，屬濕潤(rùn)季風(fēng)氣候區(qū)，平均氣溫15.7oC攝氏度，流域多年平均年降雨量975 mm，多年平均年徑流量0.6288×103m3。[9]水庫(kù)周邊濕地植被常見(jiàn)假稻、陌上菅、蘆葦和香蒲等。庫(kù)區(qū)人工林主要為楊樹(shù)，其他喬木樹(shù)種較少，林下灌木也很少，但草本層蓋度大，野大豆分布范圍廣。

1.2 數(shù)據(jù)來(lái)源與預(yù)處理

Landsat 5 TM和Landsat 8 OLI遙感影像數(shù)據(jù)和數(shù)字高程模型(DEM)數(shù)據(jù)來(lái)源于中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)(http：//www.gscloud.cn)。其中4個(gè)時(shí)期的Landsat遙感影像成像時(shí)間分別為2004年、2008年、2013年和2018年，月份在4-6月之間，軌道號(hào)為121/38，圖像質(zhì)量良好且已做過(guò)系統(tǒng)輻射較正和幾何粗校正，空間分辨率為30 m。DEM數(shù)據(jù)為ASTER GDEM數(shù)字高程數(shù)據(jù)產(chǎn)品，空間分辨率為30 m。通過(guò)DEM數(shù)據(jù)，使用ArcGIS軟件的水文分析工具提取董鋪水庫(kù)(包含大房郢水庫(kù))流域邊界，確定研究區(qū)范圍，如圖1所示。

圖1 研究區(qū)

1.3 坡度分析

地形是水土流失模型中關(guān)鍵的地表特征，坡度對(duì)地表徑流和土壤侵蝕有著重大的影響[6]。在ArcGIS軟件中對(duì)DEM數(shù)據(jù)進(jìn)行坡度計(jì)算并重分類(lèi)。坡度范圍為 0°～58.6°，其中 0°～8°的區(qū)域約占整個(gè)研究區(qū)的98%，參考《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》(SL190-2007)，將坡度分成 6類(lèi)：＜5 °，5 °-8 °，8 °-15 °，15 °-25 °，25 °-35 °，＞35 °，如圖 2 所示。

圖2 坡度圖

1.4 土地利用分類(lèi)

不同土地利用類(lèi)型對(duì)應(yīng)不同的植被覆蓋度和人類(lèi)干擾程度，從而影響水土流失的動(dòng)力和抗侵蝕阻力系統(tǒng)，在區(qū)域水土流失發(fā)展中起重要作用[1]。土地利用分類(lèi)圖采用4個(gè)年份的Landsat影像解譯結(jié)果，在ENVI軟件中執(zhí)行監(jiān)督分類(lèi)，解譯精度符合要求。參考《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》(SL190-2007)，將土地利用分類(lèi)分成3類(lèi)：水體、居民地為第1類(lèi)，農(nóng)用地為第2類(lèi)，除上述以外土地利用類(lèi)型(其他)為第3類(lèi)，如圖3所示。

1.5 植被覆蓋度估算

植被的冠層截留、削弱雨滴對(duì)地面的打擊力度、增加雨水下滲以及根系增強(qiáng)土壤的抗侵蝕力，這些有益作用成為影響水土流失敏感性的重要因素[1]。植被覆蓋度是衡量地表植被的一個(gè)最重要指標(biāo)，而像元二分法是一種廣泛應(yīng)用于植被覆蓋度的遙感估算模型，它假設(shè)一個(gè)像元的地表由有植被覆蓋和無(wú)植被覆蓋混合組成，或完全無(wú)植被覆蓋(裸土)，或完全被植被覆蓋三種情況[10]。歸一化植被指數(shù)(NDVI)能夠較好地反映植被生長(zhǎng)狀態(tài)及覆蓋信息[11]，由一個(gè)像元的NDVI值可以計(jì)算出像元中有植被覆蓋的面積比例fc。具體公式如式(1)：

式(1)中，NDVIsoil為完全是裸土(無(wú)植被覆蓋)像元的NDVI值，NDVIveg則代表完全被植被所覆蓋像元的NDVI值。

在ENVI軟件中對(duì)Landsat影像數(shù)據(jù)進(jìn)行輻射校正、大氣校正和裁剪等預(yù)處理，計(jì)算NDVI，通過(guò)土地利用類(lèi)型提取掩膜數(shù)據(jù)，運(yùn)用像元二分法估算植被覆蓋度。參考《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》(SL190-2007)，結(jié)合董鋪水庫(kù)庫(kù)區(qū)植被覆蓋的實(shí)際情況，將植被覆蓋度分成5級(jí)：1級(jí)為高度植被覆蓋(fc≥75%)，2 級(jí)為中高度植被覆蓋(60%≤fc＜75%)，3級(jí)為中度植被覆蓋(45%≤fc＜60%)，4級(jí)為中低度植被覆蓋(30%≤fc＜45%)，5級(jí)為低度植被覆蓋(fc＜30%)。4個(gè)年份的植被覆蓋度如圖4所示。

圖4 植被覆蓋度圖

1.6 水土流失風(fēng)險(xiǎn)度分類(lèi)模型

在ArcGIS軟件中，將坡度分類(lèi)、土地利用分類(lèi)和植被覆蓋度分類(lèi)三種柵格數(shù)據(jù)，按空間位置進(jìn)行空間連接，并按表1將三種影響因子交叉形成水土流失風(fēng)險(xiǎn)度分類(lèi)，參考《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》(SL190-2007)，結(jié)合董鋪水庫(kù)庫(kù)區(qū)水土流失的實(shí)際情況，分成6類(lèi)：微度，輕度，中度，強(qiáng)烈，極強(qiáng)烈和劇烈。為提高數(shù)據(jù)處理效率，采用Python語(yǔ)言編寫(xiě)分類(lèi)模型，其函數(shù)原型為 calcuValue(!lan‐duse!，!slope!，!landcover!)，其中 calcuValue為函數(shù)名稱(chēng)，括號(hào)內(nèi)為輸入?yún)?shù)，分別為土地利用分類(lèi)值，坡度分類(lèi)值和植被覆蓋度分類(lèi)值。函數(shù)體內(nèi)由42條類(lèi)似“if土地利用分類(lèi)值==3 and坡度分類(lèi)==1 and植被覆蓋度分類(lèi)==1：return微度”的決策語(yǔ)句構(gòu)成。在ArcGIS中建立分類(lèi)模型(圖 5)，輸入上述同一時(shí)期的坡度分類(lèi)、土地利用分類(lèi)和植被覆蓋度分類(lèi)數(shù)據(jù)后，即可執(zhí)行水土流失風(fēng)險(xiǎn)度分類(lèi)(圖6)。4個(gè)年份的水土流失風(fēng)險(xiǎn)度分類(lèi)結(jié)果如圖7所示。

表1 水土流失風(fēng)險(xiǎn)度分類(lèi)表

圖5 分類(lèi)模型建模

圖6 分類(lèi)模型界面

圖7 水土流失風(fēng)險(xiǎn)度

1.7 轉(zhuǎn)移矩陣

采用馬爾科夫鏈的轉(zhuǎn)移矩陣，將水土流失風(fēng)險(xiǎn)度的6種狀態(tài)兩兩結(jié)合，考察它們?cè)?004-2018年這段時(shí)期內(nèi)的變遷特征。轉(zhuǎn)移矩陣動(dòng)態(tài)變化指標(biāo)主要包括水土流失風(fēng)險(xiǎn)度類(lèi)型轉(zhuǎn)換概率(P)和水土流失風(fēng)險(xiǎn)度類(lèi)型動(dòng)態(tài)度(D)。P的計(jì)算公式為[12，13]：

式(2)中：Pij表示從一個(gè)水土流失風(fēng)險(xiǎn)度狀態(tài)轉(zhuǎn)化為另一個(gè)狀態(tài)的概率，fij是從狀態(tài)i到狀態(tài)j變遷的觀測(cè)數(shù)據(jù)，6是狀態(tài)的總數(shù)。根據(jù)式(2)計(jì)算，獲得水土流失風(fēng)險(xiǎn)度轉(zhuǎn)換概率矩陣(見(jiàn)表3)。

水土流失風(fēng)險(xiǎn)度類(lèi)型動(dòng)態(tài)度可以定量描述區(qū)域水土流失風(fēng)險(xiǎn)度變化的速度，它對(duì)比較水土流失風(fēng)險(xiǎn)度變化的區(qū)域差異和預(yù)測(cè)未來(lái)水土流失風(fēng)險(xiǎn)度變化趨勢(shì)具有積極的作用。D的計(jì)算公式為[12，13]：

式(3)中：D表示研究時(shí)段某一水土流失風(fēng)險(xiǎn)度類(lèi)型動(dòng)態(tài)度；Aa表示研究時(shí)期初某一水土流失風(fēng)險(xiǎn)度類(lèi)型的數(shù)量；Ab表示研究期末某一水土流失風(fēng)險(xiǎn)度類(lèi)型的數(shù)量；T表示研究期時(shí)間間隔，當(dāng)T的時(shí)間單位為年時(shí)，D值就是該研究區(qū)某種水土流失風(fēng)險(xiǎn)度類(lèi)型的年變化率。

2 結(jié)果與分析

根據(jù)4個(gè)時(shí)期水土流失風(fēng)險(xiǎn)度的評(píng)價(jià)結(jié)果，在Excel中統(tǒng)計(jì)了董鋪水庫(kù)庫(kù)區(qū)2004-2018年4個(gè)時(shí)期水土流失風(fēng)險(xiǎn)度數(shù)量特征(見(jiàn)表2)。由表2可以看出，盡管微度風(fēng)險(xiǎn)是研究區(qū)主要的水土流失風(fēng)險(xiǎn)度級(jí)別，土地面積總量占比達(dá)99%左右，但2004-2018年間董鋪水庫(kù)庫(kù)區(qū)水土流失風(fēng)險(xiǎn)變化較明顯，主要表現(xiàn)為水土流失風(fēng)險(xiǎn)度有減少的趨勢(shì)：極強(qiáng)烈風(fēng)險(xiǎn)面積由2.7 hm2減少到0 hm2；強(qiáng)烈風(fēng)險(xiǎn)面積由84.6 hm2減少到5.4 hm2；中度風(fēng)險(xiǎn)面積由495 hm2減少到 451.8 hm2；輕度風(fēng)險(xiǎn)面積由2976.3 hm2減少到2 254.5 hm2。微度風(fēng)險(xiǎn)面積由288 606.6 hm2增加到289 453.5 hm2。水土流失風(fēng)險(xiǎn)動(dòng)態(tài)度結(jié)果顯示：?jiǎn)我凰亮魇эL(fēng)險(xiǎn)度在2004-2008年間多發(fā)生劇烈變化，如輕度風(fēng)險(xiǎn)年變化率-13.15%，中度風(fēng)險(xiǎn)年變化率-20.36%，強(qiáng)烈風(fēng)險(xiǎn)年變化率-18.62%，極強(qiáng)烈風(fēng)險(xiǎn)年變化率-25%。2008-2018年間，微度和強(qiáng)烈風(fēng)險(xiǎn)呈加速減少趨勢(shì)，而輕度和中度呈加速增加趨勢(shì)。2004-2018年間水土流失風(fēng)險(xiǎn)總體上呈減少趨勢(shì)，微度、輕度、中度、強(qiáng)烈和極強(qiáng)烈風(fēng)險(xiǎn)年際變化率分別為0.02%，-1.73%，-0.62%，-6.69%和-7.14%。

表2 董鋪水庫(kù)2004-2018年水土流失風(fēng)險(xiǎn)度強(qiáng)烈的變化

研究區(qū)域水土流失風(fēng)險(xiǎn)度的動(dòng)態(tài)變化不僅表現(xiàn)在各種類(lèi)型的數(shù)量特征上，還可通過(guò)各種類(lèi)型之間的轉(zhuǎn)化來(lái)描述。本文通過(guò)轉(zhuǎn)移矩陣來(lái)分析各種水土流失風(fēng)險(xiǎn)度轉(zhuǎn)化的趨勢(shì)(見(jiàn)表3)。由表3可以看出，合肥市董鋪水庫(kù)庫(kù)區(qū)水土流失風(fēng)險(xiǎn)轉(zhuǎn)化主要是微度風(fēng)險(xiǎn)和其他風(fēng)險(xiǎn)度之間的轉(zhuǎn)化。2004-2008年間微度風(fēng)險(xiǎn)有99.87%沒(méi)有轉(zhuǎn)化，轉(zhuǎn)化為輕度風(fēng)險(xiǎn)和中度風(fēng)險(xiǎn)分別為0.11%和0.01%，同時(shí)輕度風(fēng)險(xiǎn)和中度風(fēng)險(xiǎn)大部分轉(zhuǎn)化為微度風(fēng)險(xiǎn)，分別為62.47%和91.09%；極強(qiáng)烈風(fēng)險(xiǎn)和強(qiáng)烈風(fēng)險(xiǎn)全部轉(zhuǎn)化為微度風(fēng)險(xiǎn)，水土保持措施成效顯著。2008-2013年間微度風(fēng)險(xiǎn)有99.71%沒(méi)有轉(zhuǎn)化，轉(zhuǎn)化為輕度風(fēng)險(xiǎn)和中度風(fēng)險(xiǎn)分別為0.25%和0.04%，同時(shí)輕度風(fēng)險(xiǎn)、中度和強(qiáng)烈風(fēng)險(xiǎn)少部分轉(zhuǎn)化為微度風(fēng)險(xiǎn)，分別為36.54%，48.04%和16.67%，輕度和強(qiáng)烈風(fēng)險(xiǎn)大部分沒(méi)有轉(zhuǎn)化，分別為60.14%和75%，水土流失風(fēng)險(xiǎn)稍微增加。調(diào)查結(jié)果表明，為保護(hù)飲用水安全，2012年間水庫(kù)周邊地區(qū)有大范圍的村民搬遷導(dǎo)致用地類(lèi)型變化。2013-2018年間微度風(fēng)險(xiǎn)有99.46%沒(méi)有轉(zhuǎn)化，轉(zhuǎn)化為輕度風(fēng)險(xiǎn)和中度風(fēng)險(xiǎn)分別為0.47%和0.07%，同時(shí)除中度風(fēng)險(xiǎn)大部分(52.92%)轉(zhuǎn)化為微度風(fēng)險(xiǎn)外，輕度風(fēng)險(xiǎn)和強(qiáng)烈風(fēng)險(xiǎn)少部分轉(zhuǎn)化為微度風(fēng)險(xiǎn)，分別為35.02%，和17.39%，水土保持有所提升，搬遷后的濕地建設(shè)、造林、綠化長(zhǎng)廊和苗木花卉基地等措施效果初顯。

表3 董鋪水庫(kù)2004-2018年水土流失風(fēng)險(xiǎn)度轉(zhuǎn)換概率矩陣

3 結(jié)論

以合肥市董鋪水庫(kù)庫(kù)區(qū)為研究區(qū)，參照中國(guó)水利部《土壤侵蝕分類(lèi)分級(jí)標(biāo)準(zhǔn)》，應(yīng)用Python編寫(xiě)水土流失風(fēng)險(xiǎn)度分類(lèi)模型，借助GIS和RS技術(shù)處理空間數(shù)據(jù)，運(yùn)用轉(zhuǎn)移矩陣和動(dòng)態(tài)度等方法定量分析了2004-2018年董鋪水庫(kù)庫(kù)區(qū)水土流失風(fēng)險(xiǎn)及其動(dòng)態(tài)變化。研究結(jié)果表明：(1)研究區(qū)域水土流失以微度風(fēng)險(xiǎn)為主，水土流失風(fēng)險(xiǎn)轉(zhuǎn)化主要是微度風(fēng)險(xiǎn)和其他風(fēng)險(xiǎn)度之間的轉(zhuǎn)化，多年間水土流失風(fēng)險(xiǎn)度沒(méi)有重大變化，但總體水土流失風(fēng)險(xiǎn)程度呈現(xiàn)下降趨勢(shì)，水土保持措施有效。(2)Python語(yǔ)言編寫(xiě)的分類(lèi)模型具有一定創(chuàng)新性，盡管分類(lèi)模型只考慮了3種水土流失因子，但可較快地評(píng)價(jià)某個(gè)區(qū)域的水土流失風(fēng)險(xiǎn)度狀況。