基于EOF法的長江流域1960—2015年降雨侵蝕力時空變化特征研究

龐延杰，董林垚，丁文峰，張平倉，朱秀迪

(1.長江水利委員會 長江科學院，湖北 武漢 430010；2.水利部 山洪地質災害防治工程技術研究中心，湖北 武漢 430010)

降雨侵蝕力反映了降雨引起土壤侵蝕的潛在能力，是一種對引起土壤侵蝕的氣候因子的定量評價指標[1]。降雨侵蝕力通常用R指標來表示，它首次出現在WISCHMEIER et al.[2]提出的建立在美國東部降雨-徑流觀測數據基礎上的通用土壤侵蝕模型(USLE)中[3]。隨后，R指標受到廣大學者的推崇，并應用于各國土壤侵蝕模型的構建和水土保持規劃工作中。R指標與降雨強度、降雨量、降雨歷時等降雨參數有關[4]，國內外學者通常使用不同時間段的次降雨總動能(E)和雨強(I)組合來計算降雨侵蝕力。考慮到這種算法是建立在基于小時或分鐘雨量信息觀測數據基礎上的，資料的獲取有一定難度，章文波等[4-5]等根據日降雨量資料構建了更簡易的算法。由于日降雨量資料相對容易獲得，且能夠提供豐富的降雨特征信息，可以近似估算降雨侵蝕力，因此被廣泛應用于降雨侵蝕力時空變化特征分析[6-8]。

區域降雨侵蝕力時空變化特征的研究對于區域水土流失量預測預報、水土保持措施制定和農業發展戰略研究具有重要的現實意義。當前國內關于區域尺度降雨侵蝕力時空變化特征的研究多側重于區域均值或單個站點降雨侵蝕力的趨勢、突變及周期分析，而對降雨侵蝕力時空變化的信息特征提取不足[9-13]。經驗正交函數(EOF)分析法在大氣、海洋等地球物理科學的數據分析領域應用廣泛，能夠將復雜的地學系統分解為若干個相互正交模態的正交組合，并提取其主要時空變化特征。EOF法被廣泛地應用于長江流域氣溫、降雨和洪澇災害時空變化特征及其對大氣環流因子的響應研究[9-11]。考慮到基于EOF法對長江流域降雨侵蝕力時空變化特征的研究成果不多，本研究以長江流域內186個雨量站的日降雨資料為基礎，應用日雨量的降雨侵蝕力模型計算長江流域1960—2015年多年平均降雨侵蝕力，使用EOF法提取降雨侵蝕力的時空變化特征，以期為長江流域水土流失研究和水土保持措施制定提供科學依據。

1 研究區域與研究數據

長江流域位于90°33′～122°25′E、24°30′～35°45′N(圖1)，流域面積180萬km2，涉及19個省、自治區和直轄市，流域人口約占全國總人口的1/3。長江流域上游為深切割高原區，中上游為切割山地區，中下游為低山丘陵與平原區。流域內年均氣溫為-4～20 ℃，并呈東高西低、南高北低的分布趨勢。流域內地域遼闊、地形復雜，季風氣候特征十分典型，降水量和暴雨的時空分布很不均勻，除江源地區年均降水量小于400 mm，四川盆地西部和東部邊緣、江西和湖南、湖北部分地區年均降水量大于1 600 mm外，大部分地區年均降水量在800～1 600 mm之間。長江流域水豐沙富，據大通站資料，年徑流總量為9 240億m3，汛期(5—10月)徑流量占年徑流總量的71.7%，年均含沙量為0.544 kg/m3，年均輸沙量為4.86億t。

本研究采用1960—2015年長江流域186個雨量站的日降雨量，數據來源于中國氣象數據網(http://data.cma.cn/site/index.html)，各雨量站數據完備，無缺測漏測。根據長江流域地理及水系分布特征，將流域劃分為11個區域：區域一為金沙江水系，包含雨量站27個；區域二為沱江水系，包含雨量站11個；區域三為嘉陵江水系，包含雨量站11個；區域四為上游干流區間，包含雨量站25個；區域五為烏江水系，包含雨量站14個；區域六為漢江水系，包含雨量站11個；區域七為中游干流區間，包含雨量站39個；區域八為洞庭湖水系，包含雨量站24個；區域九為下游干流區間，包含雨量站16個；區域十為鄱陽湖水系，包含雨量站8個；區域十一為太湖水系，無雨量站。各雨量站分布詳見圖1。

2 研究方法

2.1 降雨侵蝕力計算方法

本研究采用基于日雨量資料的降雨侵蝕力計算模型，其計算公式[4-5]為

圖1 長江流域區域劃分及雨量站分布

(1)

(2)

α=21.586β7.189 1

(3)

(4)

上式中：R為年降雨侵蝕力，MJ·mm/(hm2·h)；Ri為第i個半月時段的降雨侵蝕力，MJ·mm/(hm2·h)；k為半月時段內的降雨天數；Pj為半月時段內第j天的侵蝕性降雨量，mm，取日降雨量≥12 mm作為侵蝕性降雨的判別標準；α、β為模型參數；Pd12為日雨量≥12 mm的日平均降雨量，mm；Py12為日雨量≥12 mm的年降雨量，mm。

2.2 EOF方法

經驗正交函數(EOF)分析法也稱主成分分析法(PCA)，是一種提取主要數據特征量的方法，其計算步驟[14-15]為：

(1)對要分析的數據進行預處理，通常處理為距平形式。標準化處理后得到一個數據矩陣Xm×n，其中m為站點數，n為年數。

(2)計算矩陣Xm×n與其轉置矩陣Xn×m的交叉積，結果為相關系數矩陣Cm×m。

(3)計算Cm×m的特征根矩陣Em×m和特征向量Vm×m，其中Em×m是由特征根組成的對角矩陣，即

(5)

每個特征根λi對應特征向量Em×m的第i列，是對應的第i個EOF模態。

(4)將EOF投影到原始資料矩陣上，得到空間特征向量對應的時間系數(主成分)，即

(6)

式中：PCm×n的每行數據是對應特征向量的時間系數。

(5)第k個模態對總方差的貢獻率為

(7)

(6)EOF結果的顯著性檢驗一般通過與隨機或者虛假數據EOF分析結果進行比較，在95%置信度水平下，特征根的誤差[12]為

(8)

式中：N*為數據的有效自由度。

將特征根λi按順序依次檢查，標上誤差范圍。如果前后兩個特征根λi之間的誤差范圍無重疊，則通過顯著性檢驗。

3 結果與分析

3.1 降雨侵蝕力計算結果

利用公式(1)、(2)、(3)、(4)對長江流域186個雨量站1960—2015年歷年降雨侵蝕力進行計算，求均值后得到各個站點的多年平均降雨侵蝕力，通過Kriging插值法計算得到長江流域年均降雨侵蝕力空間分布，見圖2(a)。長江流域年均降雨侵蝕力值變化范圍為114.423～16 233.136 MJ·mm/(hm2·h)，平均值為5 855.716 MJ·mm/(hm2·h)，標準差為3 057.963 MJ·mm/(hm2·h)。相比黃土高原地區[6]和珠江流域[8]，長江流域年均降雨侵蝕力均值雖介于兩者之間，但其空間差異卻更大。長江流域降雨侵蝕力總體上從西北向東南遞增：降雨侵蝕力低值主要分布在金沙江水系、岷沱江水系北部、嘉陵江水系西北部和漢江水系東部；降雨侵蝕力高值主要分布在下游干流區間、鄱陽湖水系和洞庭湖水系。上述結果與HUANG et al.[7]對長江流域1965—2005年降雨侵蝕力分析結果大致相同。

圖2 1960—2015年長江流域多年平均降雨侵蝕力、降水量分布

由圖2可知，長江流域多年平均降雨侵蝕力、降水量的空間分布規律基本一致，說明降水量的空間分布結果對于研究降雨侵蝕力具有很好的借鑒作用。采取線性回歸方法構建長江流域年均降雨侵蝕力與經度、緯度、海拔和年均降水量之間的相關關系，結果表明(圖3)：年均降雨侵蝕力總體上隨著經度的增加而增加(R2=0.606，p<0.01)，隨著緯度的增加而減小(R2=0.096，p<0.01)，隨著海拔升高而減小(R2=0.403，p<0.01)，與年均降水量呈極顯著的正相關關系(R2=0.901，p<0.01)。

圖3 長江流域多年平均降雨侵蝕力與經度、緯度、海拔和年均降水量的相關關系

3.2 空間分布特征

應用EOF法對長江流域1960—2015年降雨侵蝕力進行時空分解，并進行顯著性檢驗，得到降雨侵蝕力的主要空間分布模態。EOF分解的特征向量貢獻率及顯著性檢驗結果表明，前4個特征向量特征值的誤差范圍不重疊，通過North顯著性檢驗，累積貢獻率達46.39%，因此前4個特征向量可以較好地解釋長江流域1960—2015年降雨侵蝕力的空間分布特征。

模態第一特征向量的方差貢獻率為24.73%，高于其他模態的貢獻率，是長江流域降雨侵蝕力的主要空間分布形式。圖4(a)顯示，除金沙江上游部分區域外，其他站點的特征值均為正值，表明1960—2015年長江流域大部分地區降雨侵蝕力變化趨勢具有高度一致性，即同時呈現高值或低值的狀況。第一模態的特征值呈現由西向東遞增的趨勢，反映長江下游地區降雨侵蝕力波動程度比上中游地區大。高值中心位于下游干流區間東部、鄱陽湖水系東北部和太湖水系，反映該區域降雨侵蝕力變化量較大。

圖4 長江流域降雨侵蝕力EOF分析第一、第二、第三和第四模態特征值空間分布

模態第二特征向量的方差貢獻率為10.75%，也是長江流域降雨侵蝕力的主要空間分布形式。圖4(b)顯示，除洞庭湖水系、鄱陽湖水系、四川盆地和長江上游局部地區外，長江流域其他地區均為正值區，因此洞庭湖水系、鄱陽湖水系、四川盆地和長江上游局部地區降雨侵蝕力會呈現與長江流域其他地區不同的震蕩類型。正值高值中心位于長江中下游河段以北區域，負值高值中心處于鄱陽湖流域東南，表明上述區域降雨侵蝕力變化較大。

模態第三特征向量和第四特征向量的方差貢獻率分別為6.18%和4.73%，反映了降雨侵蝕力空間分布的局部特征。模態第三特征向量和第四特征向量在金沙江流域、洞庭湖流域和鄱陽湖流域均呈正值，即具有相近的震蕩類型，而在其他地區正負值不同，呈現互補的變化特征。

3.3 時間變化特征

根據降雨侵蝕力空間分布特征分析結果，長江流域降雨侵蝕力場主要有4種表現類型，其中第一模態和第二模態為降雨侵蝕力的主要空間變化特征。時間系數代表了所對應特征向量空間分布模態的時間變化特征，系數的符號表示與模態震蕩類型的異同，絕對值的大小表示模態的典型程度。為了分析降雨侵蝕力場的時間變化特征，分別采用線性趨勢分析和小波變換方法分析降雨侵蝕力的變化趨勢和震蕩周期。

降雨侵蝕力前四個模態時間變化系數及線性趨勢見圖5。第一模態時間趨勢斜率大于零，其顯著性檢驗p值小于0.1，在一定程度上說明長江流域降雨侵蝕力總體上隨時間呈增加趨勢，每年的增加量約為290 MJ·mm/(hm2·h)。第一模態時間系數極大值對應年份為1998年，與長江流域強降雨引發全面洪災相對應；極小值對應年份為1978年，與長江流域中下游地區特大干旱事件對應。第二、三、四模態時間系數趨勢分析的顯著性檢驗p值均大于0.1，表明這三個模態的線性變換趨勢不明顯。該擬合結果說明利用時間變化系數可有效反映年降雨侵蝕力極值點出現的時間，同時對于預測未來幾年降雨侵蝕力的變化趨勢有一定的借鑒作用。

圖5 長江流域降雨侵蝕力第一、第二、第三和第四模態時間系數

降雨侵蝕力前四個模態的時間系數均呈周期性波動，為了探究降雨侵蝕力變化的顯著周期，使用小波變換工具對降雨侵蝕力前四個模態時間系數進行功率譜計算。第一模態時間系數波動周期集中在2～5 a；第二模態時間系數波動周期集中在2～5 a；第三模態時間系數波動周期集中在1～11 a，第四模態時間系數波動周期不顯著。第一模態時間系數小波功率譜分析結果表明長江流域降雨侵蝕力總體上呈2～5 a周期波動，在1995—2005年間波動最為明顯，與此期間全流域降雨在2～5 a周期變化相關；第二模態時間系數小波功率譜分析顯示，1990—2010年間長江流域中下游河段以北與鄱陽湖流域東南在2～5 a的周期呈干濕交替變化，引起局地降雨侵蝕力呈周期變化。

4 結 論

使用EOF法對長江流域186個雨量站1960—2015年日降雨量計算得到的降雨侵蝕力進行分析，提取其主要時空變化特征向量，采用統計方法分析其時空變化的主要特征，得出主要結論如下：

(1)長江流域1960—2015年多年平均降雨侵蝕力變化范圍為114.423～16 233.136 MJ·mm/(hm2·h)，介于黃河流域和珠江流域之間，但其空間變化比這兩者更為明顯。降雨侵蝕力總體上從西北向東南遞增，與經度和年均降水量呈正相關關系，與緯度和海拔呈負相關關系，與年均降水量相關性最強。

(2)EOF分析的前4個特征向量可以較好地解釋長江流域1960—2015年降雨侵蝕力的空間分布特征，其中前2個為主要變化特征：模態第一特征向量分析結果表明1960—2015年間長江流域大部分地區降雨侵蝕力變化趨勢具有高度一致性，下游地區降雨侵蝕力波動程度比上中游地區大；模態第二特征向量分析結果表明洞庭湖水系、鄱陽湖水系、四川盆地和長江上游局部地區降雨侵蝕力呈現與長江流域其他地區不同的震蕩類型，長江中下游河段以北區域和鄱陽湖流域東南降雨侵蝕力變化較大。

(3) 長江流域降雨侵蝕力在1960—2015年間總體上以每年290.246 MJ·mm/(hm2·h)的速率增加，第一模態時間系數的極值與長江流域旱澇事件相對應，在1990—2010年呈現2～5 a的周期性波動，與區域干濕交替變化相關。