黃忠民,肖提榮
(1.測繪出版社,北京100045;2.云南省地圖院,云南昆明650034)
龍川江、勐果河兩流域位于云南省北部、金沙江流域中游,行政區劃屬于楚雄州。綜合自然區劃為滇東高原區,地層為揚子地層區,地貌區劃為滇中紅色高原,氣候區劃為南亞熱帶半干旱氣候區,土壤區劃為山原山地紅壤地帶,植被區劃為滇中、滇東高原半濕潤常綠闊葉林、云南松林區[1]。龍川江發源于南華縣天子廟魯都拉山腳,經南華、楚雄、元謀等地于龍街匯入金沙江,全長257 km,流域面積約9293 km2;龍川江流域東、西、南三面均是高山環繞,來自東南和西南的暖濕氣流到此已是強弩之末,水汽含量很少,氣流越山后下沉,產生“焚風效應”,使空氣變得更為干燥,成為云貴高原雨量最少的地區,干濕季分明,植物蓋度低,土壤侵蝕十分嚴重,著名的“元謀土林”、“彎保土林”就出現在這里,是長江上游水土流失重災區。勐果河流域主要處于武定縣境內,集水面積約1872 km2。
考慮到云南的實際和模型的簡便、適用性,本方法以應用最為廣泛的USLE模型為藍本,即

(1)降雨侵蝕力因子(R)
就本流域而言,R值計算方法是:R=E60·I30。本研究區域內及附近各代表站降雨侵蝕力見文獻[2]中云南金沙江流域各站降雨侵蝕力計算表。
(2)土壤可蝕性因子(K)
楊子生教授曾在滇東北山區針對嚴重侵蝕的3類土壤——紅壤、黃壤和紫色土作了小區試驗,得到了實測K值,并修正了Wischmeier等諾模方程式中的系數,得到適合金沙江流域的K值計算方程式[3]:K=[2.737M1.14·(10-4) ·(12- a)+4.236(b-2)+3.259(c-3)]/100,并根據該流域土壤調查資料,按此式計算出了20個土類(亞類)的K值。
(3)坡度坡長因子(LS)
楊子生教授曾通過小區試驗的方法得到了適合金沙江流域的LS因子計算公式:LS=(L/20)0.24·(S/5°)1.32。
(4)植被覆蓋與作物經營管理因子(C)
C屬無量綱數,其值介于0~1之間[4]。楊子生教授以小區試驗方法加上參考USLE和國內一些地區的研究成果,得到了幾種自然植被類型不同地面覆蓋度的C值。
(5)水土保持因子(P)
可定義為:在其他條件相同的情況下,設某一水土保持措施的坡地土壤流失量(AP)與無任何水土保持措施的坡地土壤流失量(A)之比值,即P=AP/A,P值的大小介于0~1之間。
水土流失快速定量估算的實現涉及4個基本技術環節:科學實用的估算方法、規范的技術路線、快速的數據獲取手段、先進的數據分析處理技術。數據的獲取及更新一直是限制估算效率的瓶頸。本研究利用遙感影像來迅速獲取水土流失的動態變化信息,實現對各項估算數據的快速更新,并通過GPS檢驗其數據提取的精確性,從而解決了這一瓶頸問題。水土流失的快速估算涉及多種來源、多種類型的數據,通常包括遙感圖像(植被和土地利用)、氣候觀測資料(降雨)、專題地圖(土壤和地形)、調查資料,以及相關研究的數據產品。本文利用ASTER立體像對提取DEM,降低了經費,提高了工作效率,特別是對某些人員很難到達或經濟落后的地區,意義更加突出。整個工作流程如圖1所示。

圖1 水土流失3S定量估算方案流程
研究區域涉及ASTER圖像8景,其中L1A級的有5景,L1B級的有3景。由于L1A級數據是未經過處理的原始數據,因此在影像鑲嵌之前必須對其進行定標、輻射校正和幾何重采樣,以及與L1B級數據一起進行影像鑲嵌。
現在,隨著遙感技術的不斷發展,利用遙感影像立體像對提取DEM數據已經成為可能,如利用SPOT數據生產的DEM精度可達到10 m以內[5]。由于本項目中沒有實地的控制點(GCP),因此提取DEM只能利用ASTER DTM方式自動提取,但這樣只能提取每幅ASTER圖像的相對DEM。同時,項目又要求較為精確的DEM,因此采用SRTM的DEM來糾正提取的DEM,以使其更為精確。通過對每幅ASTER立體像對數據處理、裁剪、拼接后得到的DEM即為本文最終用來提取地形信息的DEM。
坡度坡長因子被用來衡量地形對于土壤侵蝕的影響。可利用DEM來實現基于像元的坡度與坡長的計算。在ArcGIS中,通過表面分析的坡度命令可以直接計算出每個像元的坡度值;ArcGIS軟件中沒有直接求坡長的功能,但可以通過先求負地形,再通過水文分析模塊,求出負地形的水流方向、水流長度等,再進一步利用柵格計算器把坡度、坡長代入本區域的LS因子計算公式,得到研究區的坡度坡長層。圖2是研究區坡長坡度綜合因子分級圖(注:圖中顏色越深表示地形越平坦,越不易發生水土流失;顏色淺則反之)。

圖2
將現有的土壤分布圖矢量化及與DEM配準后,根據土壤分類圖中每種類型的K值,在ArcGIS空間分析模塊中將矢量數據轉換為柵格,得到了研究區的土壤可蝕性(K)層,如圖3所示(注:圖中顏色越深表示K值越小,土壤被侵蝕的可能性越小,越不易發生水土流失;顏色越淺則反之)。

圖3
土地利用層和植被覆蓋層不直接參與水土流失量的計算。本文使用兩個時期的遙感數據來確定不同時期的土地利用狀況,這樣便可以計算出兩個不同時期的水土流失量,并通過對比分析來評價水土保持政策的實施狀況及成效。
(1)土地利用層
在對原始TM/ETM+影像進行幾何校正、影像增強、彩色合成、比值處理的基礎上,遵循以下判讀流程[6]:影像分層→監督分類→目視分析,進行分類。首先用閾值法將水體分成水體與陰影層,用比值法(NDVI)和閾值法區分植被覆蓋區與混合區;然后對照現有土地利用圖,利用ERDAS分別在TM和ETM+影像上建立土地利用類型(共10種,如林地、耕地、水體、裸巖及沙灘等)的訓練樣本,并分別對其進行監督分類及分類后處理(如聚類、殲滅等),得到兩個時期的土地利用圖;最后再與現有的土地利用圖進行疊加分析、計算。
(2)植被覆蓋層
本文在ETM+/TM上用標準化植被指數NDVI=(波段4-波段3)/(波段4+波段3)×100+127計算得到植被覆蓋層,并通過NDVI的差值了解兩時期植被動態變化信息。
根據云南金沙江流域各站降雨侵蝕力計算表中的值,結合年平均降水等線圖,使用ArcGIS的空間分析模塊中反距離權重插值函數實現降雨侵蝕力因子的空間插值,得到了研究區的降雨侵蝕力因子層(R)。
水土保持因子作為侵蝕動力的抑制因子,主要反映了植被或作物,以及管理措施對土壤流失量的影響。在土地利用圖的基礎上,依據不同作物類型及其覆蓋類型、地物類型的實測和觀察結果及已有的相關研究成果,對每一類土地利用分類中的不同植被覆蓋率進行不同的CP因子賦值。由于缺少研究區內相關的水土保持措施資料,因此此處P值取值為1。然后再在ArcGIS空間分析模塊中將矢量轉為柵格,這樣便可得到研究區在兩個時期(20世紀80年代末與90年代末)內的水土保持因子層(CP)。
RUSLE方程實際上就是將所有因子相乘從而獲得一個土壤流失值。因此,結合GIS強大的空間運算能力,便可計算研究區每個像元的土壤流失量。可利用ArcGIS空間分析模塊的柵格計算器執行疊置(A=R·K·LS·CP),然后得到像元土壤侵蝕模數分布圖,再根據像元面積計算像元年土壤侵蝕量。
由于本文中有兩時相的CP因子層,因此可以分別計算并統計出20世紀80年代末與90年代末研究區年土壤侵蝕量,并通過對比分析(見表1),實現了研究區內水土流失的動態監測。

表1 研究區兩時期年土壤侵蝕量 萬t
在像元土壤侵蝕模數分布圖的基礎上,根據國家行業標準《土壤侵蝕分類分級標準》(SL 190—2007)確定土壤侵蝕強度分級指標,對像元侵蝕量進行分類,將侵蝕量在同一侵蝕等級的像元進行歸并(重分類),進而繪制出兩時期土壤侵蝕強度圖(如圖4、圖5所示)。

圖4
由于本項目中大量應用遙感方式獲取研究區土壤侵蝕相關信息,為證實其可靠性,必須通過實地考察(路線如圖6所示)。GPS數據驗證:通過先后幾次的野外實地考察,證明由ASTER立體像對提取的DEM在實際工作中是適用的;通過3S技術監測的強度分級結果與實際情況較為相符。

圖5

圖6
本文提出的基于3S技術的水土流失定量遙感方法,首先依據監測模型計算出水土流失量,然后再按照國家相關標準分級,克服了先定級后計算流失量的人為定級的局限,是一種快速、動態、準確且較為先進的定量遙感方法。但是由于各方面原因,本文的結果將受到許多不規則因素(如工程建設)的影響,仍需進一步改進與完善。
[1]云南百科全書編纂委員會.云南百科全書[M].北京:中國大百科全書出版社,1999:15-173.
[2]楊子生.云南金沙江流域土壤流失方程研究[J].山地學報,2002,20(S1):1-9.
[3]楊子生.滇東北山區坡耕地土壤可蝕性因子[J].山地學報,1999,17(S0):10-15.
[4]WISCHMEIER W H,SMITH D D.Predicting Rainfall Erosion Losses:A Guide to Conservation Planning[M].USDA:Science and Education Adminitration,1978:40-62.
[5]黃玉琪.SPOT影像的 DEM 自動生成[J].測繪通報,1998(9):13-16.
[6]楊樹華,王寶榮,王崇云,等.流域生態系統的生態保護及其數字化管理[M].北京:科學出版社,2006:136-157.