朱大運, 楊 倩, 陳 海, 陳 靜, 李少男
(1.貴州師范大學 喀斯特研究院, 貴州 貴陽 550001; 2.國家喀斯特石漠化防治工程技術研究中心, 貴州 貴陽 550001)
雨滴濺蝕和徑流剝蝕是降雨作用于土壤地表侵蝕的主要途徑[1],可通過降雨侵蝕力指標進行表達,與降雨量、降雨強度、降雨歷時、地形植被、人類活動等因子要素密切相關[2-3]。當前全球氣候變化大背景下,極端氣候事件發生的概率明顯增多[4],某種程度上增加了土壤侵蝕的風險。科學量化評估區域降雨侵蝕力,對水土流失定向阻控和水保措施制定具有重要指導意義。
以貴州高原為中心的西南喀斯特地區是中國水土流失最嚴重的區域之一[5],其降雨侵蝕問題已引起廣泛關注。基于EI30模型[6]、次降雨模型[7]、日降雨模型[8]、月降雨模型等[9]多種降雨侵蝕力計算模型,學者對貴州降雨侵蝕力時空特征開展了不同程度的研究。許月卿等[10]首先利用19個站點的氣象數據對貴州降雨侵蝕力時空特征進行了初步分析,指出降雨侵蝕力在空間上由南到北遞減,時間上呈增加趨勢,且年內分配不勻主要集中在夏季。這與戴海倫等[11]基于貴州羅甸小區實測數據和多站點雨量資料得到的貴州降雨力時空特征相一致。阮歐等[12]基于正交函數、Mann-Kendall檢驗、小波分析等分析了貴州降雨侵蝕力的空間分布、突變性和周期性規律。Zhu等[13]研究了貴州省長時間尺度降雨侵蝕力與厄爾尼諾—拉尼娜事件的關系,結果表明ENSO持續時間與降雨侵蝕力呈正比。上述研究成果在一定程度上豐富了對貴州地區降雨侵蝕力時空分布特征的認知。但是,關于降雨侵蝕力水系分布特征卻鮮有報道,貴州喀斯特生境特征脆弱,石漠化廣泛分布且橫跨長江、珠江2個流域,局地氣候與土壤侵蝕特征較為復雜,因此揭示長時間尺度不同水系降雨侵蝕力的內在差異十分必要。本文以日降雨資料為基礎,從水系的角度對貴州省1960—2017年降雨侵蝕力時空分布特征開展對比研究,以期為區內水土流失防治、生態保護等工作提供科學參考。
貴州省地處云貴高原東部,位于103°36′—109°35′E,24°37′—29°13′N之間,屬亞熱帶濕潤季風氣候,受西南季風和東亞季風的雙重影響。年平均氣溫15 ℃左右;多年平均降水量1 100~1 400 mm,其中47%的降水集中在夏季;植被類型豐富,以亞熱帶常綠闊葉林為主。水資源分布以苗嶺為分水嶺,南部屬珠江流域,境內面積6.00×104km2,下轄南盤江水系、北盤江水系、紅水河水系和柳江水系;北部屬長江流域,境內面積1.16×105km2,下轄牛欄江水系、烏江水系、赤水河水系和沅江水系(圖1)。貴州地貌以高原山地為主,山地和丘陵面積占92.5%。地勢西高東低,起伏較大,巖溶發育典型,石漠化與水土流失十分嚴重。據水利部全國水土流失遙感調查結果和貴州省水土保持規劃(2016—2030),20世紀80年代貴州省水土流失面積占土地總面積的43.52%,經過多方持續不斷的努力,水土流失得到有效控制,到2015年,水土流失面積占比下降到27.71%,防治形勢依然嚴峻。

圖1 貴州省境內水系與氣象站點分布
貴州省1960—2017年33個國家地面標準觀測站的逐日降水量數據來源于中國氣象科學數據共享服務網(http:∥data.cma.cn/)。遵循歐洲氣候評估數據審訂標準(數據不得少于40 a,單個站點數據總缺失量不得超過10%,年數據缺失量不得超過20%或超過連續3個月)對原始數據進行了檢驗[14],各站點氣象數據均通過嚴格的質量控制和均一性檢驗。對個別缺失數據采用相鄰站點線性回歸方法進行插補,經過處理修正后的數據序列具有很好的連續性。
1.3.1 降雨侵蝕力計算 采用章文波等[15]修正的Richardson日降雨侵蝕力模型計算降雨侵蝕力,該算法已經在國內大范圍尺度和區域尺度上得到廣泛應用。具體公式如下:
(1)
式中:Ri表示第i個半月時段的侵蝕力值〔MJ·mm/(hm2·h)〕;k表示該半月時段內的侵蝕性降雨日數(d);Pj表示半月時段內第j天≥12 mm的日雨量;α與β為模型參數,分別根據以下公式計算:
α=21.586β-7.198 1
(2)
(3)
式中:Pd12表示日雨量≥12 mm的日平均雨量(mm);Py12表示日雨量≥12 mm的年平均雨量;不同區域的模型參數α與β存在差異,經計算研究區參數α的數值范圍為0.219 1~3.879 1,平均值為1.146 0;參數β的數值范圍為1.269 7~1.893 7,平均值為1.528 4。式(1)降雨侵蝕力是以逐年各半月為基本統計單元,經過累加計算得到年降雨侵蝕力值。
1.3.2 時空變化分析方法 降雨侵蝕力年際、年代際變化趨勢通過線性傾向估計來進行表達,使用最小二乘法進行傾向率估計,研究其變化率。線性傾向估計法是一種非常實用有效的表征序列變化趨勢和傾向性的分析方法,在氣候變化研究中應用廣泛。采用ArcGIS地統計模塊的反距離加權插值法,分析多年平均降雨侵蝕力空間分布特征。該插值法基于相近相似的原理,主要依賴反距離的冪值和插值點與樣本點的距離,與任何實際的物理過程均不關聯,是土壤侵蝕點狀數據到連續面狀數據轉換的一種有效方法[16]。
Mann-Kendall突變分析是世界氣象組織推薦的一種非參數的突變檢驗法,具有人為干擾性少、較驗范圍寬的優點,可以對時間序列變化趨勢的顯著性和突變性進行檢驗。它主要基于統計量序列的分析,正序列與逆序列曲線在顯著性水平邊界線內的交點即為統計序列發生突變的時間[17]。運用變異系數(CV)量化描述研究區各水系降雨侵蝕力空間分布的不均勻性,衡量序列離散程度,其公式如下:
(4)

降雨侵蝕力年內空間變化采用重心模型來表達,它基于牛頓力學,動態權衡各個地區作用力的大小并指示變量變化方向,地理位置和屬性變化是決定重心的主要因素。降雨侵蝕力重心模型參考鐘業喜等[18],計算公式如下:
(5)
(6)

1.3.3 聚類與相關性分析方法 運用聚類分析法對降雨侵蝕力空間區域進行劃分。聚類分析是對樣本或變量進行分析的一種統計方法,在多個領域的區域研究中均有應用。它根據事物本身的特性將相似的事物歸類,被歸為一類的事物具有較高的相似性,而不同類間的事物有著很大差異。本研究借助主成分分析法獲取分類數目,通過isodata算法獲取像元分組特征,從而實現對研究區降雨侵蝕力的聚類劃分。
降雨侵蝕力與地形區位因子的相關性分析以像元點為基本單位。在ArcGIS 10.2中,將降雨侵蝕力空間插值數據轉化為矢量像元點,并提取各個水系每一個像元點的降雨侵蝕力值。同時,基于SRTM-DEM和坡度、坡向數據賦予每一個像元點經緯度坐標、高程、坡度、坡向值,而后以水系為分類單元,采用Person相關性分析法計算降雨侵蝕力與有關因子的相關系數,涉及矢量像元點累計47 300多個。
2.1.1 不同流域降雨侵蝕力時間變化特征 1960—2017年貴州省境內長江流域降雨侵蝕力范圍3 049.77~8 112.10 MJ·mm/(hm2·h),多年平均值5 823.42 MJ·mm/(hm2·h);珠江流域降雨侵蝕力范圍5 933.61~7 747.38 MJ·mm/(hm2·h),多年平均值6 893.22 MJ·mm/(hm2·h),珠江流域整體降雨侵蝕能力高于長江流域。由圖2可知,過去58 a貴州省境內長江流域年降雨侵蝕力呈現出線性上升趨勢,上升速率為4.13 MJ·mm/(hm2·h·a);在年代際變化上較為復雜,經歷了升高—下降—升高—下降—升高的變化過程。珠江流域降雨侵蝕力線性趨勢以下降為主,速率為-10.66 MJ·mm/(hm2·h·a);年代際變化上表現出波動特征。然而,長江流域與珠江流域年降雨侵蝕力線性趨勢均未達到顯著性水平。距平分析表明,珠江流域降雨侵蝕力在時間序列變化上比長江流域更加劇烈。

圖2 貴州省境內珠江流域及長江流域年均降雨侵蝕力變化與距平值
2.1.2 不同水系降雨侵蝕力時間變化特征 利用線性回歸法對比分析了8個水系降雨侵蝕力變化趨勢,除沅江水系和紅水河水系表現出統計意義上的上升趨勢,其余水系均為下降趨勢,反映出貴州大部分區域降雨侵蝕風險呈降低態勢,其中沅江水系上升速率最快達18.44 MJ·mm/(hm2·h·a),北盤江水系下降率最快為-24.44 MJ·mm/(hm2·h·a)(圖3)。各水系降雨侵蝕力波動變化特征明顯且差異顯著,牛欄江水系、烏江水系、赤水河水系、沅江水系、南盤江水系、北盤江水系、紅水河水系、柳江水系極差比分別為:6.56,3.32,4.49,3.11,9.01,3.60,3.95,5.26。通過5 a滑動平均線和極差比結果可知,珠江流域下轄水系降雨侵蝕力波動幅度明顯高于長江流域,并且整體在趨勢線起伏變化上更為一致,在特定時段尤其明顯,如1965—1970,2005—2015年等。長江流域下轄的4個水系降雨侵蝕力波動幅度相對較小,但在趨勢變化的同向性上差異較大。

圖3 貴州省境內不同水系降雨侵蝕力5 a滑動變化趨勢(k為趨勢系數)
2.1.3 降雨侵蝕力時間序列突變分析 利用Mann-Kendall法對不同水系降雨侵蝕力時間序列突變特征進行檢測,并基于移動T檢驗法和Yamamoto指數法綜合對比剔除虛假突變點。檢測結果表明,過去58 a各水系降雨侵蝕力突變特征不明顯,烏江、南盤江等水系雖然檢測到突變節點,但是均未通過顯著性檢驗,僅牛欄江水系在1976年達到p<0.05顯著性水平;而沅江水系、紅水河水系以及境內長江流域則未檢測到可靠的突變點,說明其歷年降雨侵蝕力不存在自然突變(表1)。

表1 貴州省境內不同水系降雨侵蝕力時間序列突變分析
2.2.1 不同水系降雨侵蝕力空間變化特征 圖4為貴州省境內不同水系年均降雨侵蝕力空間分布圖。從圖4可知,長江流域與珠江流域范圍內降雨侵蝕力在空間變化上表現出截然不同的趨勢,前者降雨侵蝕力從西北向東南遞增,后者降雨侵蝕力則從西向東遞減。長江流域以沅江水系降雨侵蝕最為嚴重,明顯高于其他水系,3個降雨侵蝕高值中心分別出現在都勻、松桃、織金;珠江流域降雨侵蝕最嚴重的區域為南盤江水系,興義、望謨、安順是該流域降雨侵蝕高值中心。

圖4 貴州省境內不同水系年均降雨侵蝕力空間分布
各水系年均降雨侵蝕力呈現出強烈的空間變異性,其中珠江流域下轄水系變異系數差異較小,流域平均值為0.22略高于長江流域。變異系數值最大的牛欄江水系,最小的烏江水系均屬長江流域,二者相差近1.8倍,反映出該流域強烈的區域空間變異性(表2)。此外,對比發現年均降雨侵蝕力的大小與空間變異性的強弱不存在顯著的對應關系,雖然在空間變異系數最高、年均降雨侵蝕力最低的牛欄江水系表現出反向對應關系,但在其他水系均未發現這種現象,這可能是受地形、土壤等多種因素的影響。

表2 貴州省境內不同水系降雨侵蝕力變異性分析
2.2.2 基于聚類分析分區的降雨侵蝕力對比 水系是根據河道集水范圍而形成的區域地理單元,主要受地形地貌的影響。為了減少下墊面的影響,進一步分析降雨侵蝕力在自然降雨條件下空間集聚變化特征,基于聚類分析法對貴州降雨侵蝕力時空變異特征進行區域劃分和對比。范俊甫等[19]在ArcGIS 10中對貴州省58 a降雨侵蝕力進行空間插值,運用主成分分析方法實現對空間插值數據的降維,第一主成分至第六主成分量累積貢獻率達90.16%,本文以此作為isodata聚類分析中類數量劃分的依據。最后基于最大似然法將isodata聚類分析計算的特征值進行像元歸類,劃分為6個空間區域。
由圖5可知,基于聚類分析分區與流域分區的多年平均降雨侵蝕力在空間界線上呈現截然不同的特征,聚類分區更加傾向于各中心點降雨侵蝕力高低的類別特征集聚。年均降雨侵蝕力最低值為Ⅰ區3 907.59 MJ·mm/(hm2·h),最高值為Ⅳ區7 188.76 MJ·mm/(hm2·h),極差值為3 281.87 MJ·mm/(hm2·h)(表3)。與水系分類統計結果相比,雖然聚類分區的年均降雨侵蝕力極差值小于水系分區,但是兩種分類法極值倍差相同均為1.8倍,進一步表明貴州省降雨侵蝕力較強的空間差異性。變異系數除了一區值較低為0.19,其他幾個區域則相對一致,介于0.24~0.26之間(表3)。總體上,基于聚類分區的各區域變異系數比基于水系劃分的變異系數要小。

圖5 基于聚類分區的貴州省年均降雨侵蝕力空間分布

表3 1960-2017貴州省境內不同聚類分區降雨侵蝕力變異性分析
進一步對6個區域降雨侵蝕力的時間變異性進行了趨勢分析,在長時間尺度上第Ⅲ區域、第Ⅴ區域降雨侵蝕力呈現不顯著增加趨勢,增加速率分別為10.87 MJ·mm/(hm2·h·a),7.64 MJ·mm/(hm2·h·a);其余4個區域均為不顯著減少趨勢,第Ⅳ區域遞減速率最高達25.68 MJ·mm/(hm2·h·a)(圖6)。1960—2017年各區域年均降雨侵蝕力變化曲線走向上整體上較為一致,在急驟轉換節點表現尤為明顯,如1990,2015年前后的降雨侵蝕力快速增加趨勢,但是波動幅度存在差別,反映出各區域對氣候響應的差異性。

圖6 貴州省基于聚類分區的降雨侵蝕力時間變化(k為趨勢系數)
基于重心模型計算了貴州境內各月降雨侵蝕力中心和大雨中心的位置,以及重心遷移方向(圖7)。如圖7所示,從地理空間上來看,降雨侵蝕力重心主要分布在貴陽轄區與周邊的龍里、貴定縣,降雨侵蝕力重心滯留在龍里縣境內時間最長達5個月,其次為貴陽市轄區。1—12月份降雨侵蝕力重心位置在這3個區域呈不規則方向遷移,累計遷移距離191 km,平均每月移動15.9 km,9—10月份重心遷移距離最長,單月向東平移了31.5 km,7—8月份重心遷移緩慢僅2.8 km。降雨侵蝕最為強烈的6月至9月其重心坐標全部位于貴陽市境內,與大雨(≥25 mm)的月重心位置一致。大雨月重心位置分布相對集中,除了降水較少的12,1,2月,全部集聚在貴陽東北的烏當區境內。

圖7 貴州降雨侵蝕力與大雨月重心遷移
通過降雨量推演的降雨侵蝕力僅是雨水對土壤侵蝕潛在能力的反映,它與地形地貌、人類活動等要素共同作用決定了區域水土流失狀況。貴州是中國典型喀斯特區,地貌特征脆弱而復雜,相關性分析表明地形、區位等要素對局地降雨侵蝕力的影響作用明顯(表4)。各水系年均降雨侵蝕力與經度、緯度相關性均達到顯著水平(p<0.05),與經度的相關性主要表現為正相關,其中柳江水系相關系數最高為0.89;與緯度的相關性主要表現為負相關,北盤江水系相關系數最高為-0.79。海拔與降雨侵蝕力相關性密切,大部分水系表現為顯著負相關,其中區域海拔較高的牛欄江水系和南、北盤江水系與降雨侵蝕力的負相關系數最高,這可能是受下墊面的影響。李維杰等[17]研究發現降雨侵蝕力與海拔負相關主要受地形以及地勢起伏的影響,海拔較高、坡面較陡、地勢起伏較大的區域較為明顯,這與牛欄江、南、北盤江水系下墊面地貌特征基本吻合。
由表4可知,坡向與降雨侵蝕力的相關性不顯著,但是坡度與降雨力的相關性全部通過顯著性水平檢驗,其中西北部牛欄江水系相關程度最高,存在較為嚴重的降雨侵蝕風險。

表4 不同水系降雨侵蝕力與影響因素相關性分析
降雨侵蝕力不僅直接影響土壤侵蝕的時空分布,而且作為重要因子參與侵蝕模型的計算,為水土流失預測提供關鍵支撐[3,15]。本文基于貴州省33個氣象站點逐日降水資料,分析了貴州省境內8個水系近60 a降雨侵蝕力時空變化特征,結果表明各水系年降雨侵蝕力呈波動變化,在長時間尺度上以線性減少趨勢為主。蘆鑫[20]研究發現長江流域上游年均降雨侵蝕力呈下降趨勢;賴成光等[21]對珠江流域降雨侵蝕力的分析發現,南、北盤江水系呈減少趨勢,紅水河水系呈增加趨勢,與本研究結果趨于一致。
全球氣候變化與各水系降雨侵蝕力時空演變關系密切,當氣溫上升,蒸發旺盛,大氣環流與降水空間格局在全球尺度上發生重新調整,進而導致侵蝕外營力降雨量與降雨強度發生改變[22]。貴州位于中國第二階梯向第三階梯過渡地帶,高海拔、低緯度的區位特征,使其氣候變化復雜且受季風影響明顯。趙志龍等[23]分析結果表明,1960—2016年貴州省降水量年際變化劇烈,呈不顯著減少趨勢,與本文降雨侵蝕力變化特征十分吻合,說明降雨量減少是引起侵蝕力呈遞減變化的關鍵因素。然而,降水量的減少只是區域氣候變化的一種直觀反映。
從深層次原因來看,過去幾十年中國近海海溫升高和青藏高原熱源作用減弱,季節性海陸溫差變化異常,導致東亞季風環流減弱、夏季西太平洋副熱帶高壓位置偏南,來自東部海洋的水汽無法深入長江流域西部地區[24];同時期西南季風的減弱現象也被研究發現[25]。因此,東亞季風與西南季風減弱的疊加影響,導致貴州地區降水減少,可能是引起貴州降雨侵蝕力呈線性遞減的主要原因。
此外,關于降雨侵蝕力在長時間尺度上的波動變化特征,厄爾尼諾—南方濤動(ENSO)的影響不容忽視。有學者研究發現貴州年降雨侵蝕力與ENSO指代因子海洋尼諾指數(ONI)、多變量ENSO指數(MEI)、南方濤動指數(SOI)關系密切,并且降雨侵蝕力大小與ENSO的持續時間呈正比[13]。然而,降雨侵蝕力變化是多種影響要素共同作用的結果,后續還需要從氣候變化、人類活動、地理因子等綜合影響的角度進行更多的分析。
(1) 貴州省大部分水系降雨侵蝕力呈現小幅波動下降趨勢,年際變化趨勢系數介于-24.44~18.44 MJ·mm/(hm2·h·a)之間,季風變化是導致降雨侵蝕力差異化的主要原因。珠江流域下轄水系降雨侵蝕力普遍高于長江流域,但各水系突變特征不明顯,僅牛欄江水系檢測到可靠突變點。
(2) 珠江流域范圍內各水系降雨侵蝕力空間變異性高于長江流域水系,在空間分布上從西向東遞減,高值中心在興義和望謨。與珠江流域不同,長江流域降雨侵蝕力從東南向西北遞減,高值中心在織金、都勻。降雨侵蝕力和地形區位關系密切,與經緯度、海拔、坡度均表現出較好的相關性,但是與坡向的相關關系不顯著。
(3) 聚類分析分區與水系分區的區域降雨侵蝕力在變化趨勢上相對一致,以遞減為主;從變異系數來看,基于聚類分區的各區域變異性普遍小于各水系。年內降雨侵蝕力重心主要集中在貴陽及周邊的龍里、貴定地區,其遷移路徑與區域大雨重心較為相似。