沂蒙山區(qū)降雨侵蝕力空間分布推算方法

楊韶洋,劉霞,姚孝友,吳迪,張榮華,齊斐,張洪達

(1.山東省土壤侵蝕與生態(tài)修復(fù)重點實驗室,山東農(nóng)業(yè)大學(xué)林學(xué)院水土保持系,271018,山東泰安;2.美國密西西比州立大學(xué)森林資源系,39759,美國密西西比州斯塔克維爾;3.南京林業(yè)大學(xué)林學(xué)院,210000,南京;4.水利部淮河水利委員會水土保持監(jiān)測中心站,233001,安徽蚌埠)

降雨侵蝕力(rainfall erosivity)是評價降雨引起土壤侵蝕潛在能力的一個動力指標(biāo)[1],研究區(qū)域降雨侵蝕力空間插值方法,對提高預(yù)測精度,了解區(qū)域降雨侵蝕力空間分布規(guī)律具有重要意義。空間插值方法主要有確定性插值方法(反距離權(quán)重插值、樣條曲線插值)和地統(tǒng)計插值方法(克里格插值),其中,地統(tǒng)計插值方法既充分考慮了樣本點的方向、位置和距離,又能夠?qū)?shù)據(jù)中存在的趨勢和異向性進行處理,選出最優(yōu)、最適合的模型進行擬合[2]。因此,近年來國內(nèi)外諸多學(xué)者采用地統(tǒng)計插值方法中的普通克里格法(Ordinary Kriging)對不同區(qū)域降雨侵蝕力的空間分布特征進行了研究[3-8],并取得了一定的成果;但在土壤養(yǎng)分、溫度空間預(yù)測的相關(guān)研究探索中發(fā)現(xiàn),由多元線性回歸和普通克里格法結(jié)合而成的回歸克里格法(Regression Kriging)比普通克里格法更能提高空間預(yù)測精度[9-15],而目前應(yīng)用回歸克里格法對降雨侵蝕力空間分布特征的研究鮮見報道。

筆者借助沂蒙山區(qū)88個雨量站點1980—2010年日降水?dāng)?shù)據(jù)和降雨侵蝕力日雨量公式,計算站點年均降雨侵蝕力,采用回歸克里格法和普通克里格法對沂蒙山區(qū)降雨侵蝕力進行空間插值預(yù)測,探討不同插值方法對降雨侵蝕力空間預(yù)測精度的影響,并對沂蒙山區(qū)降雨侵蝕力及其空間分布特征進行研究,以期為北方土石山區(qū)尤其是魯中南低山丘陵區(qū)降雨侵蝕力研究及區(qū)域水土流失定量評價提供技術(shù)支撐。

1 研究區(qū)概況

沂蒙山區(qū)位于淮河流域北部、山東省的中南部,屬北方土石山區(qū),面積3萬2 571.29 km2,涉及27個縣(市、區(qū))。屬半濕潤暖溫帶季風(fēng)氣候,年均氣溫12～14℃,年均降水量為700～900 mm,主要集中在夏季汛期(6—9月),汛期降雨量為592 mm,占全年降水量的71.3%;河流水系為沂沭泗河水系,總長度2 586 km。地貌類型以山地丘陵為主,山丘區(qū)面積占63%;主要土壤類型有棕壤、褐土、潮土、砂漿黑土、粗骨土;地帶性植被屬暖溫帶闊葉林帶,森林植被主要有油松(Pinus tabulaeformis),赤松(Pinus densiflora)、側(cè)柏(Platycladus orientalis)刺槐(Robinia pseudoacacia)、麻櫟(Quercus acutissima)等,灌草植被主要有黃荊(Vitex negundo)、酸棗(Ziziphus jujuba)、結(jié)縷草(Zoysia japonica)等。

2 材料與方法

2.1 數(shù)據(jù)來源與處理

以沂蒙山區(qū)88個雨量站點(圖1)1980—2010年日降水資料為基礎(chǔ)數(shù)據(jù)源,采用章文波等[16]日雨量修正模型[17](式(1)和式(2)),計算得到年均降雨侵蝕力實算值,并統(tǒng)計各站點年均降雨量、年均汛期(6—9月)降雨量及年均侵蝕性降雨量(日降雨量≥12 mm)?；贏rcGIS地統(tǒng)計分析功能,采用普通克里格法獲得年均降雨量、年均汛期降雨量和年均侵蝕性降雨量插值圖。

年均降雨侵蝕力計算公式為:

式中:R半月k為第 k 半月的降雨侵蝕力,MJ·mm/(hm2·h·a);Pdij為第i年第k半月第j日≥12 mm的日降雨量;α、β為回歸系數(shù);Pd12為日降雨量≥12 mm的日平均值,mm;Pdl表示統(tǒng)計時段內(nèi)第l日≥12 mm的日降雨量;j=1,2,…,m為第i年第k半月日降雨量≥12 mm的時間,d;i=1,2,…,N為時間,a;l=1,2,…,n為統(tǒng)計時段內(nèi)所有日降雨量≥12 mm的時間,d。為多年平均年降雨侵蝕力,MJ·mm/(hm2·h·a);k=1,2,…,24,表示1年24個半月。

2.2 降雨侵蝕力空間插值方法

2.2.1 普通克里格插值 以年均降雨侵蝕力為區(qū)域化變量,半方差函數(shù)為分析工具,對插值點區(qū)域化變量的取值進行線性無偏最優(yōu)估計[18-19]。

2.2.2 回歸克里格插值 以年均降雨侵蝕力為目標(biāo)變量、降雨特征值為輔助變量建立回歸方程,利用回歸方程對趨勢項進行預(yù)測,計算趨勢項殘差并其進行普通克里格插值,最終將趨勢項與殘差項進行空間相加得到目標(biāo)變量預(yù)測值。

式中:z(so)為未知點so目標(biāo)變量,即降雨侵蝕力回歸克里格預(yù)測值;m(so)為趨勢項,即降雨侵蝕力回歸方程預(yù)測值,通過線性回歸進行擬合;ε(so)為殘差項(降雨侵蝕力實算值-回歸預(yù)測值),利用普通克里格插值預(yù)測;so為未知點,o=1,2,…,n。

2.3 不同方法預(yù)測精度評價

隨機抽取13個站點作為數(shù)據(jù)精度驗證子集,驗證剩余75個站點的插值預(yù)測結(jié)果精度。通過比較驗證站點的實算值與預(yù)測值,采用式(5)和(6)計算平均絕對誤差和均方根預(yù)測誤差[20],并以均方根誤差作為指標(biāo),利用式(7)計算相對精度改進值[21],評價預(yù)測精度。

式中:M為平均絕對誤差;m為驗證數(shù)據(jù)子集采樣點數(shù)目;Zoi、Zpi分別為雨量站點i的實算值和預(yù)測值;R為均方根預(yù)測誤差;RI為回歸克里格插值相對于普通克里格插值的相對精度改進值;RR為普通克里格插值的均方根預(yù)測誤差值;RRK為回歸克里格插值的均方根預(yù)測誤差值。

圖1 沂蒙山區(qū)雨量站點分布圖Fig.1 Distribution map of rainfall stations in Yimeng Mountain Area

2.4 空間分布差異性分析

根據(jù)降雨侵蝕力插值結(jié)果,按照同一分級標(biāo)準(zhǔn)分類( ＜3 300 MJ·mm/(hm2·h·a)(1 級),≥3 300 ～3 700 MJ·mm/(hm2·h·a)(2 級), ＞ 3 700 ～ 4 100 MJ·mm/(hm2·h·a)(3 級), ＞ 4 100 ～ 4 500 MJ·mm/(hm2·h·a)(4 級), ＞ 4 500 ～ 4 900 MJ·mm/(hm2·h·a)(5 級), ＞4 900 ～5 300 MJ·mm/(hm2·h·a)(6級)和 ＞5 300 MJ·mm/(hm2·h·a))(7 級),從小到大依次編號為1～7級?；贏rcGIS空間分析模塊,疊加生成圖譜和面積轉(zhuǎn)移矩陣,分析不同插值方法降雨侵蝕力空間分布差異性。

3 結(jié)果與分析

3.1 降雨侵蝕力與降雨特征值

由表1可知,年均降雨侵蝕力均值為4 133.92 MJ·mm/(hm2·h·a);變化范圍在 3 033.23 ～5 438.22 MJ·mm/(hm2·h·a)之間,標(biāo)準(zhǔn)差為567.59,變異系數(shù)達13.73%。根據(jù)反映離散程度變異系數(shù)的分級標(biāo)準(zhǔn)[22]可知,沂蒙山降雨侵蝕力屬中等變異。

表1 沂蒙山區(qū)降雨侵蝕力統(tǒng)計特征值Tab.1 Rainfall erosivity statistical characteristics of Yimeng Mountain Area MJ·mm/(hm2·h·a)

圖2 沂蒙山區(qū)降雨特征值空間分布圖Fig.2 Rainfall characteristics values distribution map of Yimeng Mountain Area

根據(jù)沂蒙山區(qū)降雨特征值空間分布圖(圖2)可知,年均降雨量、年均汛期降雨量以及年均侵蝕性降雨量的變化范圍分別在457.27～851.64、457.27～613.42和379.82～671.46 mm之間,各降雨特征值在空間分布上存在差異性,整體上遵循從南向西北和東北2個方向逐漸遞減的特征。

3.2 普通克里格插值分析

隨機選取75個雨量點的年均降雨侵蝕力實算值作為區(qū)域化變量,采用普通克里格法插值生成年均降雨侵蝕力空間分布圖(圖3),分析可知,降雨侵蝕力均值為4 214.85 MJ·mm/(hm2·h·a),變化范圍為3 293.30 ～5 271.71 MJ·mm/(hm2·h·a),其空間分布規(guī)律與降水特征值分布規(guī)律一致。但插值結(jié)果的均值、最大值和最小值與實算值分別相差80.93、166.51 和 260.07 MJ·mm/(hm2·h·a),表明普通克里格插值結(jié)果與實算值差異較大。

圖3 沂蒙山區(qū)普通克里格法降雨侵蝕力空間分布圖Fig.3 Rainfall erosivity map of the Ordinary Kriging method in Yimeng Mountain Area

3.3 回歸克里克插值分析

3.3.1 降雨侵蝕力與降雨特征值相關(guān)性分析 降雨侵蝕力與各降雨特征值呈極顯著相關(guān)(表2),其中與年均汛期降雨量相關(guān)性最高,相關(guān)系數(shù)達0.917,表明降雨侵蝕力隨年均降雨量、年均汛期雨量和年均侵蝕性降雨量的增加而增大。年均降雨量與降雨侵蝕力的相關(guān)系數(shù)小于年均汛期雨量和年均侵蝕性降雨量,表明年均降雨量對降雨侵蝕力的影響小于年均汛期降雨量和年均侵蝕性降雨量。

3.3.2 回歸預(yù)測模型 采用多元線性逐步回歸的方法,基于上述75個插值站點各降雨特征值擬合計算降雨侵蝕力變化特征。方程擬合結(jié)果如下:

模型1:m=13.892x1-3 426.750,R2=0.839(P＜0.05);

表2 沂蒙山區(qū)降雨侵蝕力與降雨特征值相關(guān)性Tab.2 Correlations between the rainfall erosivity and rainfall characteristics of Yimeng Mountain Area

模型2:m=9.906x1+2.822x2-2 819.957,R2=0.888(P ＜0.05)。

式中:m為年均降雨侵蝕力回歸方程預(yù)測值;x1為年均汛期降雨量;x2為年均侵蝕性降雨量。

從擬合方程來看,模型1與模型2均為逐步回歸結(jié)果,且模型2的決定系數(shù)大于模型1;因此選擇模型2作為最優(yōu)擬合結(jié)果,進行降雨侵蝕力擬合計算。同時所選降雨特征值中,年均降雨量被剔除,說明雖然年均降雨量與降雨侵蝕力極顯著相關(guān),但其不能很好反映降雨侵蝕力的空間分布特征,而年均汛期降雨量和年均侵蝕性降雨量均被保留,說明二者是影響降雨侵蝕力空間預(yù)測的主要因素。

3.3.3 殘差分析 由殘差統(tǒng)計特征值(表3)可以看出,回歸預(yù)測殘差均值為 -0.37 MJ·mm/(hm2·h·a),變化范圍在 -389.97 ～356.88 MJ·mm/(hm2·h·a)之間,殘差值變化范圍較大,存在較強的變異性。對預(yù)測殘差進行半方差函數(shù)分析,結(jié)果(表4)表明,指數(shù)模型基臺效應(yīng)和均方根誤差均小于其他2個模型,因此選取指數(shù)模型作為插值模型。降雨侵蝕力預(yù)測殘差基臺效應(yīng)為25.73%,屬中等空間相關(guān)性[23],表明年均汛期降雨量和年均侵蝕性降雨量作為結(jié)構(gòu)因素是引起降雨侵蝕力空間變異的主要原因,而隨機因素影響較小。

表3 沂蒙山區(qū)降雨侵蝕力回歸預(yù)測殘差統(tǒng)計特征值Tab.3 Regression prediction residual statistical characteristic values of rainfall erosivity of Yimeng Mountain Area

表4 降雨侵蝕力預(yù)測殘差半方差函數(shù)模型及其參數(shù)Tab.4 Prediction residual semi variance function models and parameters of rainfall erosivity

3.3.4 插值結(jié)果分析 分析年均降雨侵蝕力插值分布圖(圖4)可知,采用回歸克里格插值后降雨侵蝕力均值為4 176.55 MJ·mm/(hm2·h·a),變化范圍為3 065.27 ～5 413.70 MJ·mm/(hm2·h·a),其空間分布規(guī)律與降雨特征值空間分布規(guī)律一致。插值結(jié)果的均值、最大值和最小值與實算值分別相差42.63、24.52 和 32.04 MJ·mm/(hm2·h·a),說明回歸克里格插值結(jié)果與實算值相近。

表5 不同插值方法預(yù)測精度對比Tab.5 Precision comparation of different interpolation methods

3.4 不同插值方法對比分析

3.4.1 插值精度評價 采用剩余13個驗證站點來比較2種插值方法的預(yù)測精度,從預(yù)測精度(表5)可知,回歸克里格法的平均絕對誤差(108.69)小于普通克里格法的平均絕對誤差(250.22),同時各驗證點絕對誤差(圖5)表明回歸克里格法要明顯低于普通克里格法;對比均方根誤差可知,回歸克里格法均方根誤差為139.11,明顯小于普通克里格法的300.05,回歸克里格法相對精度改進值為53.64%,說明回歸克里格插值結(jié)果比較理想,其預(yù)測精度高于普通克里格插值。

3.4.2 空間分布差異性分析 根據(jù)降雨侵蝕力插值分級結(jié)果空間疊加圖譜(圖6)和降雨侵蝕力插值圖層空間疊加轉(zhuǎn)移矩陣(表6)所示,2種方法插值結(jié)果中,降雨侵蝕力等級一致的面積為2萬4 794.51 km2,占總面積的76.12%,即疊加圖譜中的白色區(qū)域,說明2種插值方法的降雨侵蝕力空間分布規(guī)律具有一致性;但進一步分析發(fā)現(xiàn),2種方法插值結(jié)果在各級別上存在差異,其中:普通克里格插值結(jié)果分級中缺少第7級,而回歸克里格插值結(jié)果中存在,且普通克里格插值結(jié)果的第6級中有1 285.58 km2在回歸克里格結(jié)果中屬于第7級,即疊加圖譜中的正方形網(wǎng)格區(qū);其次普通克里格插值結(jié)果的第1級面積僅為1.42 km2,且在回歸克里格中屬于第2級,即變化圖譜中的豎線區(qū),而回歸克里格插值結(jié)果的第1級面積為372.67 km2,在普通克里格中屬于第2級,即疊加圖譜中的斜網(wǎng)格區(qū)。

圖4 沂蒙山區(qū)回歸克里格法降雨侵蝕力空間分布圖Fig.4 Rainfall erosivity map of the Regression Kriging method in Yimeng Mountain Area

圖5 不同方法預(yù)測絕對誤差分析圖Fig.5 Prediction absolute error analysis chart of different interpolation methods

上述分析說明,雖然2種插值結(jié)果在空間分布規(guī)律上具有一致性,但普通克里格插值結(jié)果在各降雨侵蝕力級別中分布不完全,對局部變異特別是極值區(qū)體現(xiàn)不夠充分,而回歸克里格插值結(jié)果除了在數(shù)值上與降雨侵蝕力實算值接近外,其插值結(jié)果對于局部變異體現(xiàn)更為精確。從空間分布的情況來看,回歸克里格插值結(jié)果更為接近現(xiàn)實情況,說明其對于空間分布預(yù)測的效果更好。

圖6 不同方法降雨侵蝕力插值分級結(jié)果空間疊加圖譜Fig.6 Spatial overlay map of rainfall erosivity classify results of different interpolation methods

表6 不同方法降雨侵蝕力插值圖層空間疊加轉(zhuǎn)移矩陣Tab.6 Area transfer matrix of the rainfall erosivity interpolation maps of different methods km2

3.5 降雨侵蝕力空間分布特征

根據(jù)沂蒙山區(qū)降雨侵蝕力回歸克里格插值圖(圖4)所示:沂蒙山區(qū)年均降雨侵蝕力最小值出現(xiàn)在沂源縣北部,最大值出現(xiàn)在蒼山縣中部,降雨侵蝕力整體呈現(xiàn)由南向西北和東北2個方向遞減的特征,這與區(qū)域內(nèi)降雨特征值分布規(guī)律一致。同時,在鄒城東部和莒縣東部存在2個高值中心,在沂南東北部和山亭區(qū)北部存在2個低值中心。由各縣降雨侵蝕力分布特征(表7)可知:沂源縣降雨侵蝕力均值最小為3 449.60 MJ·mm/(hm2·h·a),臨沂市羅莊區(qū)降雨侵蝕力均值最大為5 324.29 MJ·mm/(hm2·h·a),除沂源縣外微山縣、沂水縣、莒縣和蒙陰縣等9縣降雨侵蝕力均值低于 4 000 MJ·mm/(hm2·h·a);鄒城市、新泰市、五蓮縣和泗水縣等9縣,降雨侵蝕力均值在 4 000 ～4 500 MJ·mm/(hm2·h·a)之間;日照市嵐山區(qū)、棗莊市臺兒莊區(qū)、費縣等6縣,降雨侵蝕力均值在4 500 ～5 000 MJ·mm/(hm2·h·a)之間;臨沂市蘭山區(qū)、蒼山縣和羅莊區(qū)降雨侵蝕力均值大于5 000 MJ·mm/(hm2·h·a) 。

沂蒙山區(qū)降雨侵蝕力變異系數(shù)為12.01%,說明沂蒙山區(qū)降雨侵蝕力在空間分布屬中等變異;但各縣(市、區(qū))除鄒城市屬中等變異外,其他各縣變異系數(shù)均小于10%,屬弱變異,說明鄒城市降雨侵蝕力的空間變異性大于其他各縣,主要原因是鄒城市存在降雨侵蝕力極值中心。

表7 沂蒙山區(qū)各縣降雨侵蝕力統(tǒng)計特征值Tab.7 Rainfall erosivity statistical characteristics of different counties in Yimeng Mountain Area MJ·mm/(hm2·h·a)

4 結(jié)論與討論

1)回歸克里格法插值預(yù)測結(jié)果與降雨侵蝕力實算值相近。與普通克里格法相比,回歸克里格法的相對精度改進值為53.64%;雖然2種方法降雨侵蝕力插值結(jié)果在空間分布規(guī)律上具有一致性,但回歸克里格插值結(jié)果對降雨侵蝕力分布規(guī)律的描述更為精確。

2)沂蒙山區(qū)降雨侵蝕力與降雨特征值多元線性逐步回歸結(jié)果表明,年均降雨量并不是影響降雨侵蝕力空間分布的主要因素,年均汛期雨量和年均侵蝕性雨量等結(jié)構(gòu)因素是影響降雨侵蝕力空間分布的主要因素,而隨機因素對降雨侵蝕力的影響較小。

3)沂蒙山區(qū)降雨侵蝕力在空間分布上遵循從南向西北和東北2個方向逐漸遞減的分布規(guī)律,與區(qū)域內(nèi)降雨分布規(guī)律一致。沂蒙山區(qū)27個縣(市、區(qū))中,除鄒城市外其他各縣降雨侵蝕力在空間分布上均屬于弱空間變異,且各縣降雨侵蝕力空間變異性均小于沂蒙山區(qū)。

本研究采用回歸克里格法僅對沂蒙山區(qū)降雨侵蝕力進行了空間預(yù)測,并與普通克里格法對比分析,但回歸克里格法是否在其他區(qū)域也具有相同的效果尚不明確;同時本研究僅對降雨特征值與降雨侵蝕力相關(guān)性進行分析,對于降雨特征值如何影響降雨侵蝕力還有待進一步研究。雖然沂蒙山區(qū)各縣(市、區(qū))降雨侵蝕力在空間分布上屬于弱變異性,但縣域間降雨侵蝕力存在明顯差異;因此各縣應(yīng)根據(jù)區(qū)域降雨侵蝕力分布特點及實際情況合理進行水土流失預(yù)防和治理。沂蒙山區(qū)各縣降雨侵蝕力空間變異性均小于沂蒙山區(qū),說明降雨侵蝕力空間變異可能與區(qū)域尺度有關(guān),如何確定不同區(qū)域尺度對降雨侵蝕力空間變異的影響還有待深入研究。