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

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

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

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

1 研究區(qū)概況

沂蒙山區(qū)位于淮河流域北部、山東省的中南部,屬北方土石山區(qū),面積3萬(wàn)2 571.29 km2,涉及27個(gè)縣(市、區(qū))。屬半濕潤(rùn)暖溫帶季風(fēng)氣候,年均氣溫12～14℃,年均降水量為700～900 mm,主要集中在夏季汛期(6—9月),汛期降雨量為592 mm,占全年降水量的71.3%;河流水系為沂沭泗河水系,總長(zhǎng)度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ù)來(lái)源與處理

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

年均降雨侵蝕力計(jì)算公式為:

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

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

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

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

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

2.3 不同方法預(yù)測(cè)精度評(píng)價(jià)

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

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

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

2.4 空間分布差異性分析

根據(jù)降雨侵蝕力插值結(jié)果,按照同一分級(jí)標(biāo)準(zhǔn)分類( ＜3 300 MJ·mm/(hm2·h·a)(1 級(jí)),≥3 300 ～3 700 MJ·mm/(hm2·h·a)(2 級(jí)), ＞ 3 700 ～ 4 100 MJ·mm/(hm2·h·a)(3 級(jí)), ＞ 4 100 ～ 4 500 MJ·mm/(hm2·h·a)(4 級(jí)), ＞ 4 500 ～ 4 900 MJ·mm/(hm2·h·a)(5 級(jí)), ＞4 900 ～5 300 MJ·mm/(hm2·h·a)(6級(jí))和 ＞5 300 MJ·mm/(hm2·h·a))(7 級(jí)),從小到大依次編號(hào)為1～7級(jí)?；贏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ù)達(dá)13.73%。根據(jù)反映離散程度變異系數(shù)的分級(jí)標(biāo)準(zhǔn)[22]可知,沂蒙山降雨侵蝕力屬中等變異。

表1 沂蒙山區(qū)降雨侵蝕力統(tǒng)計(jì)特征值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個(gè)方向逐漸遞減的特征。

3.2 普通克里格插值分析

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

圖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ù)達(dá)0.917,表明降雨侵蝕力隨年均降雨量、年均汛期雨量和年均侵蝕性降雨量的增加而增大。年均降雨量與降雨侵蝕力的相關(guān)系數(shù)小于年均汛期雨量和年均侵蝕性降雨量,表明年均降雨量對(duì)降雨侵蝕力的影響小于年均汛期降雨量和年均侵蝕性降雨量。

3.3.2 回歸預(yù)測(cè)模型 采用多元線性逐步回歸的方法,基于上述75個(gè)插值站點(diǎn)各降雨特征值擬合計(jì)算降雨侵蝕力變化特征。方程擬合結(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ù)測(cè)值;x1為年均汛期降雨量;x2為年均侵蝕性降雨量。

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

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

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

表4 降雨侵蝕力預(yù)測(cè)殘差半方差函數(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é)果的均值、最大值和最小值與實(shí)算值分別相差42.63、24.52 和 32.04 MJ·mm/(hm2·h·a),說(shuō)明回歸克里格插值結(jié)果與實(shí)算值相近。

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

3.4 不同插值方法對(duì)比分析

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

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

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

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

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

圖6 不同方法降雨侵蝕力插值分級(jí)結(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個(gè)方向遞減的特征,這與區(qū)域內(nèi)降雨特征值分布規(guī)律一致。同時(shí),在鄒城東部和莒縣東部存在2個(gè)高值中心,在沂南東北部和山亭區(qū)北部存在2個(gè)低值中心。由各縣降雨侵蝕力分布特征(表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ū)、棗莊市臺(tái)兒莊區(qū)、費(fèi)縣等6縣,降雨侵蝕力均值在4 500 ～5 000 MJ·mm/(hm2·h·a)之間;臨沂市蘭山區(qū)、蒼山縣和羅莊區(qū)降雨侵蝕力均值大于5 000 MJ·mm/(hm2·h·a) 。

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

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

4 結(jié)論與討論

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

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

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

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