基于GIS的黑龍江省降雨侵蝕力時空變化分析

周 寧，李 超，滿秀玲

（1.東北林業(yè)大學 林學院，黑龍江 哈爾濱150040；2.黑龍江省水土保持科學研究院，黑龍江 哈爾濱150070；3.北京林業(yè)大學，北京100083）

黑龍江省全境屬廣義上的東北黑土區(qū)，土壤侵蝕面積1.12×105km2，占總土地面積的24.74%，主要侵蝕類型包括水力侵蝕、風力侵蝕、凍融侵蝕和工程侵蝕，總體上以輕度水力侵蝕為主［1］。按照水土保持區(qū)劃劃分為大興安嶺山地水源涵養(yǎng)生態(tài)維護區(qū)，小興安嶺山地丘陵生態(tài)維護保土區(qū)，三江平原—興凱湖生態(tài)維護農(nóng)田防護區(qū)，長白山山地水源涵養(yǎng)減災區(qū)，長白山山地丘陵水質(zhì)維護保土區(qū)，東北漫川漫崗土壤保持區(qū)，松遼平原防沙農(nóng)田防護區(qū)和大興安嶺東南低山丘陵土壤保持區(qū)等8個分區(qū)。黑龍江省已成為我國重要商品糧基地中土壤侵蝕最嚴重的地區(qū)之一，正在逐步喪失其作為重要農(nóng)業(yè)生產(chǎn)基地的基礎［1］。鑒于水土流失對黑龍江省生態(tài)環(huán)境、糧食生產(chǎn)安全和社會經(jīng)濟發(fā)展造成的嚴重影響，針對黑龍江省黑土區(qū)土壤侵蝕展開深入的研究十分必要。

通用土壤流失方程［2－4］中的降雨侵蝕力因子作為衡量區(qū)域降雨侵蝕潛勢的量度，既受氣候條件的影響，也受時空變化的影響。張憲奎［5］等通過分析認為，黑龍江省的最佳降雨侵蝕力指標是E60I30，建立了降雨侵蝕力的簡化計算方法。謝云［6］等比較了用次降雨與日降雨量計算降雨侵蝕力的關(guān)系，得出在我國可以用日降雨來計算降雨侵蝕力。章文波［7］等比較了多種雨量資料估算降雨侵蝕力，得出日降雨量計算多年平均降雨侵蝕力精度最高。我國學者［8－12］對黃土高原區(qū)降雨侵蝕力變化規(guī)律和時空演變特征開展了研究，而對黑龍江省降雨侵蝕力時空研究還鮮見報道，劉燕玲［13］等僅對黑龍江省16氣象站點的降雨侵蝕力進行了研究。本研究采用黑龍江省近40a的日降雨量數(shù)據(jù)，依托GIS平臺，分析省域范圍內(nèi)降雨侵蝕力的時間和空間特異性。與以往同類研究相比，本研究著重于GIS工具的運用，時間序列劃分半月、10a和40a，從年際和多年角度進行變化分析；空間跨越尺度上以橫、縱山脈及松嫩和三江兩大平原為分界線，從反映垂直分布的地形角度進行趨勢分析。研究成果可應用于黑龍江省水土保持規(guī)劃、預報、監(jiān)測和治理工作，對加強水土流失治理，保護自然生態(tài)環(huán)境，保障糧食生產(chǎn)安全和社會經(jīng)濟發(fā)展具有積極意義。

1 材料與方法

1.1 數(shù)據(jù)來源和處理

中國氣象科學數(shù)據(jù)共享服務網(wǎng)公布的中國地面氣候資料日值數(shù)據(jù)集收錄了黑龍江省自1951年以來的基本、基準地面氣象觀測站及自動站日值數(shù)據(jù)集，其原始數(shù)據(jù)文件經(jīng)過較嚴格的質(zhì)量控制和檢查，并對1971—2000年的部分統(tǒng)計值進行了極值檢驗和時間一致性檢驗，并對已查出的錯誤記錄進行更改。選取地面氣候資料日值數(shù)據(jù)集中黑龍江省全部32個氣象站點日降雨量數(shù)據(jù)，通過初步統(tǒng)計，1972年以后數(shù)據(jù)序列較為完整，同時保證數(shù)據(jù)分析的時段性及時新性，確定選取1972—2011年32個氣象臺站的日降水量數(shù)據(jù)，除肇州氣象站數(shù)據(jù)時限為1988年1月1日至2011年12月31日、鶴崗氣象站數(shù)據(jù)時限為1972年1月1日至2008年12月31日外，其余各氣象臺站數(shù)據(jù)時限為1972年1月1日至2011年12月31日。保留降雨量數(shù)據(jù)，篩除雪、霧、露和霜等形式的降水量數(shù)據(jù)。

1.2 研究方法

從每年的1月1日至12月31日，每半月為一個計算區(qū)間，一年共有24個區(qū)間，1月上半月記為1-1，下半月記為1-2，以此類推，12月上半月記為12-1，下半月記為12-2。采用日降雨量計算半月降雨侵蝕力［14］。

式中：Ri——第i個半月的降雨侵蝕力值〔MJ·mm／（hm2·h）〕；n——半月內(nèi)天數(shù)，每月上半月天數(shù)均為15，下半月按實際計?。籔j——第i個半月第j日的日侵蝕性降雨量（mm）；α，β——參數(shù)值；Pd——日侵蝕性降雨量的日平均值（mm）；Py——日侵蝕性降雨量的年平均值（mm）。日侵蝕性降雨量取≥12mm 的日降雨量［6］。

逐年計算各氣象臺站的半月降雨侵蝕力，將結(jié)果導入ArcGIS構(gòu)建基礎數(shù)據(jù)庫。采用Kriging方法對平均降雨侵蝕力進行插值［15］，插值柵格大小為250m×250m。采用趨勢分析法［16］，從時間和空間上對黑龍江省降雨侵蝕力演變特征進行分析。

2 結(jié)果與分析

2.1 半月平均降雨侵蝕力時空分布特征

黑龍江省32個氣象臺站40a半月平均降雨侵蝕力呈現(xiàn)出相似的時間分布特征。圖1為12個主要氣象站半月平均降雨侵蝕力分布情況。由圖1可以看出，因篩除了雪、霧、露和霜等形式的降水量，降雨侵蝕力值在冬季、秋末及春初為零，呈單峰型集中分布在4月下半月至10月上半月，先逐漸增加，7—8月出現(xiàn)最大值，隨后逐漸降低。4月下半月至10月上半月平均降雨侵蝕力極差和期望值（圖2）表明，半月平均降雨侵蝕力在空間序列上存在較大差異，7和8月差異最為顯著。

圖1 黑龍江省1972－2011年主要氣象臺站半月平均降雨侵蝕力

圖2 黑龍江省4月下半月至10月上半月平均降雨侵蝕力極差及期望值

由黑龍江省半月平均降雨侵蝕力空間分布圖（圖略）可知，4月下半月平均降雨侵蝕力，伊勒呼里山以北地區(qū)出現(xiàn)最小值，小興安嶺、張廣才嶺、老爺嶺地區(qū)呈較小值，西部松嫩平原呈較大值，東部三江平原出現(xiàn)最大值。5月上半月平均降雨侵蝕力出現(xiàn)多個極大值中心，最小值出現(xiàn)在新林地區(qū)，黑龍江全省范圍由南向北、自西至東呈增大趨勢，最大值出現(xiàn)在海倫地區(qū)。5月下半月平均降雨侵蝕力，伊勒呼里山北部出現(xiàn)最小值，嫩江平原西南部、小興安嶺北部呈較小值，中部、東部地區(qū)呈較大值，張廣才嶺—完達山南部出現(xiàn)最大值。6月上半月平均降雨侵蝕力空間分布差異開始增大，最小值出現(xiàn)在漠河地區(qū)，北部、南部呈較小值，自西向東呈增大趨勢，最大值出現(xiàn)在鶴崗地區(qū)。6月下半月平均降雨侵蝕力，最小值出現(xiàn)在漠河地區(qū)，除北部和東部局部地區(qū)呈較小值外，全省范圍大部分地區(qū)呈較大值，最大值出現(xiàn)在海倫地區(qū)，并向周邊地區(qū)呈輻射分布減小。7月上半月平均降雨侵蝕力空間分布差異最為明顯，最小值出現(xiàn)在漠河地區(qū)，北部和東部大部分地區(qū)呈較小值，最大值出現(xiàn)在綏化地區(qū)，并向周邊地區(qū)呈遞減分布。7月下半月平均降雨侵蝕力，呈全年最大值，除大興安嶺地區(qū)呈現(xiàn)由北向南遞增分布外，全省范圍自西向東遞減趨勢明顯。8月上半月平均降雨侵蝕力，最小值出現(xiàn)在塔河地區(qū)，最大值出現(xiàn)在伊春地區(qū)，除小興安嶺—張廣才嶺地區(qū)、三江平原東部呈較大值，其他地區(qū)呈較小值。8月下半月平均降雨侵蝕力，最小值出現(xiàn)在塔河地區(qū)，三江平原呈較小值，東部山區(qū)呈較大值，最大值出現(xiàn)在鶴崗地區(qū)，小興安嶺向松嫩平原呈遞減趨勢。9月上半月平均降雨侵蝕力，松嫩平原出現(xiàn)最小值，張廣才嶺地區(qū)呈較小值，小興安嶺北部、老爺嶺嶺東部、三江平原呈較大值。9月下半月平均降雨侵蝕力，最小值出現(xiàn)在伊勒呼里山北部，小興安嶺北部、松嫩平原呈較小值，小興安嶺—張廣才嶺東部、三江平原呈較大值，最大值出現(xiàn)在鶴崗地區(qū)。10月上半月平均降雨侵蝕力，松嫩平原出現(xiàn)最小值，北部、西部呈較小值，小興安嶺—張廣才嶺東部出現(xiàn)3個極大值中心，最大值出現(xiàn)在鶴崗地區(qū)。

2.2 黑龍江省10a平均降雨侵蝕力時空分布特征

如圖3所示，以10a為刻度劃分時間序列進行分析，呼瑪、富裕和齊齊哈爾地區(qū)降雨侵蝕力變化趨勢相同，1972—1981年平均降雨侵蝕力最小，1982—1991年平均降雨侵蝕力最大值，1992—2001年平均降雨侵蝕力呈現(xiàn)波谷，2002—2011年平均降雨侵蝕力略有增大。大興安嶺、孫吳和佳木斯地區(qū)降雨侵蝕力變化趨勢相同，1972—1981年平均降雨侵蝕力最小，1992—2001年平均降雨侵蝕力最大。泰來、綏化和肇州地區(qū)降雨侵蝕力變化趨勢相同，2002—2011年平均降雨侵蝕力最小，1982—1991年平均降雨侵蝕力最大。海倫、伊春和雞西地區(qū)降雨侵蝕力變化趨勢相同，1992—2001年平均降雨侵蝕力最小，1982—1991年平均降雨侵蝕力最大。塔河、鶴崗、明水、鐵力、哈爾濱、牡丹江和綏芬河地區(qū)降雨侵蝕力變化趨勢相同，1972—1981年平均降雨侵蝕力最小，1982—1991年平均降雨侵蝕力出現(xiàn)最大值后呈持續(xù)降低趨勢。

圖3 黑龍江省10a平均降雨侵蝕力空間分布

采用Kriging方法對多年平均降雨侵蝕力進行插值。結(jié)果表明，黑龍江省范圍1972—1981年平均降雨侵蝕力最小，為1 207.55 〔MJ·mm／（hm2·h）〕，1982—1991年平均降雨侵蝕力最大，為1 436.99〔MJ·mm／（hm2·h）〕，1992—2001年、2002—2011年較1982—1991年平均降雨侵蝕力呈持續(xù)遞減趨勢。由黑龍江省10a平均降雨侵蝕力空間分布圖（圖略）可知，40a間每10a平均降雨侵蝕力最小值均出現(xiàn)在伊勒呼里山以北地區(qū)。1972—1981年平均降雨侵蝕力，小興安嶺北部、老爺嶺東南部、松嫩平原西部地區(qū)呈較小值，松嫩平原中南部、三江平原中部出現(xiàn)兩個極大值中心，并呈輻射狀向周邊遞減。1982—1991年平均降雨侵蝕力，三江平原呈較小值，黑龍江全省范圍出現(xiàn)4個極大值中心，小興安嶺—張廣才嶺向松嫩平原呈遞減趨勢。1992—2001年平均降雨侵蝕力，松嫩平原和三江平原呈較小值，南部出現(xiàn)了明顯的最大值中心，向周邊遞減趨勢顯著。2002—2011年平均降雨侵蝕力，張廣才嶺、老爺嶺地區(qū)呈較小值，三江平原南部出現(xiàn)局部較大值，小興安嶺東北部出現(xiàn)最大值中心，向松嫩平原遞減分布。

2.3 黑龍江省40a平均降雨侵蝕力時空分布特征

通過SPSS計算，1972—2011年平均降雨侵蝕力P值等于0.2＞0.05（W 檢驗），且在ArcGIS上得到的QQPlot分布圖中數(shù)據(jù)呈一條直線，服從正態(tài)分布，無需進行數(shù)據(jù)轉(zhuǎn)換即可進行插值，即可分析得出為40a平均降雨侵蝕力空間分布情況?；诓逯到Y(jié)果進行半方差分析，逐年分別計算塊金值和基臺值，塊金值平均值為667 932.91，基臺值平均值為1 681 401.34，塊金效應［17］（塊金值與基臺值的比值）為39.72%，介于25%～75%之間，表明多年平均降雨侵蝕力具有中等強度的空間相關(guān)性。分析32個站點1972—2011年平均降雨侵蝕力空間分布趨勢，可見全省范圍多年平均降雨侵蝕力，自西向東、由北至南呈明顯的拋物線狀分布趨勢。伊勒呼里山以北地區(qū)呈極小值，并出現(xiàn)最小值，以中部地區(qū)形成的最大值中心向四周呈遞減分布，在小興安嶺東部、張廣才嶺南部有極大值出現(xiàn)（圖4）。1972—2011年平均降雨侵蝕力總體上呈微弱的上升趨勢。降雨侵蝕力等值線（圖5）顯示，伊勒呼里山向北、松嫩平原自東向西、三江平原自西向東、小興安嶺由南至北、張廣才嶺、老爺嶺由北至南，多年平均降雨侵蝕力遞減趨勢明顯。

圖4 黑龍江省1972－2011年平均降雨侵蝕力空間分布

3 結(jié)論

（1）黑龍江省半月平均降雨侵蝕力在時間分布上呈單峰型，從4月下半月至10月上半月表現(xiàn)為先逐漸增加，7和8月出現(xiàn)最大值，隨后逐漸降低。半月平均降雨侵蝕力在空間分布上存在較大差異，7月差異最為顯著，7月下半月平均降雨侵蝕力為全年的最大值。

圖5 黑龍江省氣象臺站分布及1972－2011年平均降雨侵蝕力等值線

（2）以10a平均降雨侵蝕力為時間序列，黑龍江全省范圍1972—1981年平均降雨侵蝕力最小，1982—1991年平均降雨侵蝕力最大，1992—2001年、2002—2011年較1982—1991年平均降雨侵蝕力呈持續(xù)遞減趨勢。近40a平均降雨侵蝕力，伊勒呼里山向北、松嫩平原自東向西、三江平原自西向東、小興安嶺由南至北、張廣才嶺、老爺嶺由北至南，呈明顯遞減趨勢分布。

（3）多年平均降雨侵蝕力具有較強的空間相關(guān)性，全省范圍多年平均降雨侵蝕力，自西向東、由北至南呈明顯的拋物線狀分布趨勢，時間演進上呈微弱的上升趨勢。

（4）本研究采用日降雨量計算降雨侵蝕力的方法仍可優(yōu)化，在取得降雨侵蝕力實際值的基礎上選擇適合黑龍江省的計算公式，并通過增加氣象臺站數(shù)量，提高計算精度，將能夠取得更準確的黑龍江省降雨侵蝕力分布數(shù)據(jù)。就降雨侵蝕力在時空分布上表現(xiàn)出來的特異性，仍需進一步分析，今后將探究其受土地利用、氣候條件、植被類型和地形地貌等因素的影響程度。

［1］ 水利部，中國科學院，中國工程院.中國水土流失防治與生態(tài)安全：東北黑土區(qū)卷［M］.北京：科學出版社，2010：254-269.

［2］ Wischmeier W H.Prdicting Rainfall-erosion losses from cropland east of the Rocky Mountains［M］∥USDA.Agricultural Handbook，1965.

［3］ Wischmeier W H，Smith D D.Prdicting Rainfall Erosion Losses：A Guide to Conservation Planing［M］.Washington：Department of Agriculture，1978.

［4］ Renard K G，F(xiàn)oster G R，Weesies G A，et al.Predicting Soil Erosion by Water：A Guide to Conservation Planning with the Revised Universal Soil Loss Equation（RUSLE）［M］.Washington：Department of Agriculture，1997.

［5］ 張憲奎，許靖華.黑龍江省土壤流失方程的研究［J］.水土保持通報，1992，12（4）：1-9，18.

［6］ 謝云，劉寶元.用日雨量和雨強計算降雨侵蝕力［J］.水土保持通報，2001，21（6）：53-56.

［7］ 章文波，付金生.不同類型雨量資料估算降雨侵蝕力［J］.資源科學，2003，25（1）：35-41.

［8］ 殷水清，謝云.黃土高原降雨侵蝕力時空分布［J］.水土保持通報，2005，25（4）：29-33.

［9］ 李靜，劉志紅，李銳.黃土高原不同地貌類型降雨侵蝕力時空特征研究［J］.水土保持通報，2008，28（3）：124-127.

［10］ 信忠保，許炯心，馬元旭.近50a黃土高原侵蝕性降水的時空變化特征［J］.地理科學，2009，29（1）：98-104.

［11］ 戴海倫，苑爽，張科利，等.貴州省降雨侵蝕力時空變化特征研究［J］.水土保持研究，2013，20（1）：37-41.

［12］ 劉建祥，李智廣，張曉萍，等.最近30a河龍區(qū)間降雨侵蝕力的時空演變特征［J］.中國水土保持科學，2013，11（2）：6-13.

［13］ 劉燕玲，劉濱輝，王立剛，等.黑龍江省降雨侵蝕力的變化規(guī)律［J］.中國水土保持科學，2010，8（2）：24-29.

［14］ 章文波，謝云，劉寶元.利用日雨量計算降雨侵蝕力的方法研究［J］.地理科學，2002，22（6）：705-711.

［15］ 萬龍，馬芹，張建軍，等.黃土高原降雨量空間插值精度比較：Kriging與TPS法［J］.中國水土保持科學，2011，9（3）：79-87.

［16］ 陳晶，楊帆，田寶星，等.黑龍江省黑土區(qū)近50a來四季氣溫時空變化特征分析［J］.水土保持研究，2013，20（3）：150-154.

［17］ 張敏，賀鵬飛，陳偉強.基于GIS和地統(tǒng)計學的土壤養(yǎng)分空間變異分析［J］.東北農(nóng)業(yè)大學學報，2010，41（3）：53-58.