基于TRMM數(shù)據(jù)提取降雨侵蝕力的校正方法
——以薊運河上游地域為例

靳秋桐,史明昌,?,張舉濤,王珊,胡影

(1.北京林業(yè)大學水土保持學院,100083,北京；2.北京林業(yè)大學水土保持與荒漠化防治教育部重點實驗室,100083,北京；3.北京地拓科技有限責任公司,100084,北京；4.高德軟件有限公司,102200,北京)

基于TRMM數(shù)據(jù)提取降雨侵蝕力的校正方法
——以薊運河上游地域為例

靳秋桐1,2,史明昌1,2,3?,張舉濤1,2,王珊4,胡影3

(1.北京林業(yè)大學水土保持學院,100083,北京；2.北京林業(yè)大學水土保持與荒漠化防治教育部重點實驗室,100083,北京；3.北京地拓科技有限責任公司,100084,北京；4.高德軟件有限公司,102200,北京)

由于TRMM數(shù)據(jù)的時間分辨率較低,導致用TRMM數(shù)據(jù)計算某一區(qū)域指定時間段內的降雨侵蝕力存在較大的系統(tǒng)誤差。為解決該問題提出通過統(tǒng)計研究區(qū)氣象站實際降雨過程的數(shù)據(jù),或者通過模擬研究區(qū)典型降雨過程,確定不同降雨過程中不同時間分辨率下降雨侵蝕力的差異程度,并以此對每一場降雨的降雨侵蝕力進行修正,進而校正指定時間段的降雨侵蝕力。為驗證該方法的有效性,以薊運河上游地區(qū)為研究區(qū),在Matlab 2010b中對當?shù)氐湫徒涤赀^程進行模擬,計算當?shù)夭煌瑫r間分辨率下降雨侵蝕力差異程度,并利用典型降雨年份2008年的TRMM 3B42數(shù)據(jù)(0.25°×0.25°空間分辨率及3 h時間分辨率)通過ENVI-IDL語言程序進行校正計算。結果表明,經過校正后,研究區(qū)年降雨侵蝕力分布約為900～2 900 MJ·mm/(hm2·h)。該結果與這一地區(qū)多年平均降雨侵蝕力的量級及分布一致,說明該校正方法是一種切實可行的、能有效提高TRMM數(shù)據(jù)應用精度的方法；但是受TRMM數(shù)據(jù)空間分辨率的限制,該方法只適用于研究大中尺度區(qū)域降雨侵蝕力。

降雨侵蝕力；TRMM；系統(tǒng)誤差；校正方法；薊運河流域

近些年來全球氣候變化導致極端天氣頻繁發(fā)生,我國很多地區(qū)旱澇災害事件劇增,這加重了我國土壤侵蝕嚴重程度[1]。降雨是引起土壤侵蝕的重要因子,由降雨引起土壤侵蝕的潛在能力,即為降雨侵蝕力[2]5-8[3]。準確計算和評估降雨侵蝕力對各種土壤侵蝕模型的設計、模型精度的評價及準確定量預報土壤侵蝕工作都有重要的意義[4-6]。

降雨侵蝕力是用降雨的總動能和降雨強度估算的多年平均指標[2]50-51[6]。指定時間長度內,每場降雨的總動能E和對應場次最大降雨強度I的乘積EI的累計值作為該段時間內的降雨侵蝕力[3,7]。通常降雨資料從氣象站點獲取[8-9],但氣象站點數(shù)量有限且分布不均,這使獲取空間面域上的降雨侵蝕力分布情況變得極為困難[5,10]；此外,及時準確地評價和預報土壤侵蝕,需要及時準確的降雨數(shù)據(jù)去估算降雨侵蝕力,這些也是在現(xiàn)階段技術條件下氣象站難以滿足的[9,11]。1997年11月美國和日本合作開發(fā)的全球第一顆搭載降水雷達的熱帶降雨觀測衛(wèi)星(Tropical Rainfall Measuring Mission,TRMM)發(fā)射成功,使獲取準全球空間面域上的降雨分布即時數(shù)據(jù)成為現(xiàn)實,彌補了氣象站觀測數(shù)據(jù)的不足[11-13]。

學者們對TRMM數(shù)據(jù)在測算降雨強度和降雨量時的精度問題爭議較大。C.M.Mannaerts等[14]認為TRMM數(shù)據(jù)的精度可用于大尺度空間降雨的研究；而T.Y.Pham[15]認為TRMM數(shù)據(jù)由于空間分辨率太低,并不能得到較為準確的降雨分布。研究區(qū)尺度越小,用TRMM數(shù)據(jù)計算結果的誤差越大。TRMM數(shù)據(jù)不宜用于小尺度區(qū)域的研究,用其研究中等尺度空間范圍內降雨侵蝕力也很少,通常只用其研究大尺度空間范圍內的降雨量和降雨侵蝕力的分布[9-11,13,16]?，F(xiàn)有TRMM數(shù)據(jù)計算降雨侵蝕力的方法分為2種,第1種是根據(jù)TRMM數(shù)據(jù)估算的降雨量結合研究區(qū)的經驗公式間接計算,第2種是根據(jù)TRMM數(shù)據(jù)利用RUSLE的公式直接計算。第1種方法因計算過程簡單,被廣泛使用但是其可靠性差,誤差不能估計。如Fan Jianrong等[9]選用經驗公式利用TRMM計算的降雨量數(shù)據(jù)去估算降雨侵蝕力分布。第2種方法計算過程較為復雜,但是能夠準確卻反映研究區(qū)降雨侵蝕力變化趨勢；為解決TRMM數(shù)據(jù)精度帶來的誤差,Zhu Qiang等[11]利用TRMM數(shù)據(jù)3 h統(tǒng)計的平均降雨強度和最大降雨強度代替30 min統(tǒng)計的平均降雨強度和最大降雨強度。A.Vrieling等[10,13,16]也簡單認為2種時間分辨率下的最大降雨強度相等；但是這些處理方式與實際情況不符(詳見第2節(jié)),從而降低了研究的可信度。

筆者基于以往利用TRMM數(shù)據(jù)計算降雨侵蝕力方法的研究,找到了一種有效提高TRMM數(shù)據(jù)在中等尺度上應用精度的方法,這既擴大TRMM數(shù)據(jù)的應用范圍,且彌補部分地區(qū)地面氣象觀測的空白。此外,筆者在用TRMM數(shù)據(jù)計算降雨侵蝕力的方法基礎上分析了該方法的系統(tǒng)誤差,提出了一種可行的校正方法,并以薊運河流域上游地區(qū)為研究區(qū),利用典型降雨年份2008年的TRMM 3B42數(shù)據(jù),對該校正方法進行了驗證。

1 降雨侵蝕力計算方法

降雨強度是一個實時變化的物理量,時間分辨率越高,監(jiān)測到的降雨強度變化越快,變化的幅度越大,統(tǒng)計到的降雨強度極值越大；反之,時間分辨率越低,變化幅度越小,統(tǒng)計到的極值越小。降雨侵蝕力是在平降雨強度和最大降雨強度的基礎上計算得到。與降雨強度一樣,降雨侵蝕力也受到時間分辨率的影響。為此RUSLE手冊建議,用每30 min統(tǒng)計的平均降雨強度和最大降雨強度計算降雨侵蝕力[3,7],在降雨侵蝕力的研究計算中,只針對具有侵蝕性的降雨進行研究和計算(通常將降雨量＞12 mm的次降雨認定為侵蝕性降雨[17])。

1.1 地面氣象站降雨資料計算降雨侵蝕力

地面氣象站降雨監(jiān)測數(shù)據(jù)的時間分辨率通常為30 min。RUSLE手冊建議,利用這種數(shù)據(jù)計算降雨侵蝕力的公式[3,18]如下:

式中:R為年平均降雨侵蝕力,MJ·mm/(hm2·h)；n為計算R值得總年數(shù)；j為參與計算的年份序數(shù)；m為第j年內的暴雨次數(shù)；k為參與計算的次暴雨序數(shù)；E為次降雨總動能,MJ/hm2；I30為30 min時間分辨率資料統(tǒng)計的次降雨最大30 min降雨強度,mm/h；i30為30 min時間分辨率資料統(tǒng)計的次降雨平均降雨強度,mm/h；ΔV為對應時間段內的降雨量,mm。

1.2 TRMM數(shù)據(jù)計算降雨侵蝕力

TRMM3B42數(shù)據(jù)集由TRMM科學資料信息系統(tǒng)和TRMM辦公室聯(lián)合制作(在http:∥trmm.gsfc. nasa.gov/或 http:∥old-cdc.cma.gov.cn/免費獲取),是基于紅外亮溫資料采用3B-42算法得到的準全球的降雨估量數(shù)據(jù),其圖像的一個像素點代表0.25°×0.25°經緯度范圍,像素值為每3 h平均降雨強度資料,單位mm/h。即TRMM數(shù)據(jù)得到降雨強度資料的時間分辨率為180 min。

基于TRMM的降雨侵蝕力計算方法是根據(jù)RUSLE手冊建議采用的降雨侵蝕力計算公式推導,采用TRMM 3B42中的3 h降雨量資料,篩除次降雨量小于12 mm的非侵蝕性降雨,采用K.G.Renard等[3]和Zhu Qiang等[11]的計算方法計算次降雨侵蝕力、月降雨侵蝕力和年降雨侵蝕力,計算公式如下:

式中:e為次降雨動能,MJ/(mm·hm2)；Rs為次降雨侵蝕力,MJ·mm/(hm2·h)；i180為180 min時間分辨率下的次降雨平均降雨強度,mm/h；I180為TRMM資料中的次降雨最大180 min降雨強度,mm/h；p為次降雨持續(xù)3 h段數(shù),il為該次降雨第l個3 h段的平均降雨強度。計算時,將該次降雨中il的最大值作為次降雨最大180 min降雨強度,將il的平均值作為180 min時間分辨率下的次降雨平均降雨強度。

對于月降雨侵蝕力,則是將一個月內的次降雨侵蝕力疊加得到:

式中:Rm為月降雨侵蝕力,MJ·mm/(hm2·h)；q為該月次降雨的次數(shù)。

對于年降雨侵蝕力,則是將一年內的次降雨侵蝕力疊加而得到:

式中:Ry為年降雨侵蝕力,MJ·mm/(hm2·h)；r為該年次降雨的次數(shù)。

1.3 TRMM數(shù)據(jù)計算降雨侵蝕力分布

編寫ENVI-IDL語言程序一次性讀入研究范圍內的TRMM 3B42數(shù)據(jù),利用式(3)～(6)計算出每個網格(0.25°×0.25°)站點的降雨量以及降雨侵蝕力。這樣得到的降雨侵蝕力可以表述大尺度上的降雨分布,但結果粗糙。通常利用地統(tǒng)計學方法對格網上的降雨侵蝕力進行克里格插值,得到更精細平滑的降雨侵蝕力分布。

為了準確構建插值模型,利用ArcGIS中探索性空間數(shù)據(jù)分析(ESDA)工具對數(shù)據(jù)進行檢測及篩選。其中:使用自然間斷點的分類符號方法對TRMM數(shù)據(jù)初步計算結果進行繪制和檢查；利用趨勢面分析任務,將經緯度坐標(x,y)與降雨侵蝕力的值(z)構成的點對投影在南北和東西2個平面上,得到的投影點分別進行趨勢線回歸擬合,得到南北和東西2個方向平面上的降雨侵蝕力變化趨勢線,即得到趨勢面分析圖。根據(jù)趨勢分析的結果決定克里金插值是否進行趨勢剔除處理。

協(xié)同克里格方法可以給出有限區(qū)域內區(qū)域變量的最佳線性無偏估計量,在遙感圖像處理中也得到了越來越廣泛的應用[19]。因降雨受地形影響較大[20],所以結合研究區(qū)域的數(shù)字高程模型(DEM)可以對格網上降雨侵蝕力進行協(xié)同克里格插值。在對TRMM數(shù)據(jù)進行克里格插值之前,需要將上一節(jié)中計算出的每個網格的降雨侵蝕力值作為該網格中心點位置的降雨侵蝕力,然后進行插值[10-11,13,16]。

2 降雨侵蝕力計算方法的校正

2.1 TRMM數(shù)據(jù)獲取降雨侵蝕力時的誤差來源

由遙感影像得到的TRMM數(shù)據(jù)源計算的降雨侵蝕力結果必然存在一定程度的誤差。該誤差分為2類:隨機性誤差和系統(tǒng)性誤差。隨機性誤差由TRMM衛(wèi)星遙感獲取數(shù)據(jù)時降雨強度大小及取樣時間的隨機性引起,也由不同尺度凸顯的地學要素(地形、地貌及地理位置等)影響的隨機性引起。隨機性誤差不可避免,尤其地學要素對降雨侵蝕力總體趨勢的影響會隨著尺度的縮小被逐步放大。系統(tǒng)性誤差由TRMM衛(wèi)星搭載的測雨雷達的精度水平、監(jiān)測時間間隔,及TRMM數(shù)據(jù)計算降雨侵蝕力的計算方法差異導致。由設備精度帶來的誤差隨著技術水平的改進可以提高,由計算方法帶來的系統(tǒng)誤差可以估計并消除。

利用不同時間分辨率氣象站降雨資料得到的降雨侵蝕力大小不一樣[8,18]。研究TRMM數(shù)據(jù)計算的降雨侵蝕力精度及誤差值,需要選定評判的標準,以往研究中通常采用 30 min降雨資料作為基準[3,7]。因為氣象站數(shù)據(jù)只是針對特定空間點的降雨數(shù)據(jù),所以對于面域的降雨分布來說,氣象站定點數(shù)據(jù)作為TRMM數(shù)據(jù)計算結果的基準并不合適；故不能對由隨機性造成的誤差進行準確的估計。筆者只分析由計算方法引起的系統(tǒng)誤差。

2.2 不同分辨率數(shù)據(jù)計算降雨侵蝕力方法的差異

受TRMM 3B42數(shù)據(jù)時間分辨率3 h的限制,不能獲取式(1)和式(2)中所需的30 min時間分辨率統(tǒng)計的次降雨最大降雨強度I30和平均降雨強度i30；因此,在利用TRMM數(shù)據(jù)計算降雨侵蝕力時,采用180 min分辨率下統(tǒng)計的次降雨平均降雨強度i180和最大降雨強度I180來替代較難獲得的30 min時間分辨率統(tǒng)計的次降雨最大降雨強度I30和平均降雨強度i30,利用式(3)和式(4)來計算次降雨侵蝕力。

不同的計算方法應用于同一地區(qū)會存在不同(甚至相差數(shù)倍)的降雨侵蝕力值[6]。降雨量、歷時及降雨時程(雨量隨歷時分配過程)的不同對一場降雨的降雨強度有顯著的影響。對于一場歷時恰好為30 min的侵蝕性降雨來說,如果在降雨開始時刻計時統(tǒng)計,最大降雨強度I30是I180的6倍,而平均降雨強度i30也是i180的6倍。由此可見,氣象站數(shù)據(jù)和TRMM數(shù)據(jù)計算降雨侵蝕力所得結果的差異,其實質上是不同時間分辨率下(最大和平均)降雨強度的差異[18]。在這種情況下,確定30 min與180 min時間分辨率下降雨強度的關系是校正在TRMM數(shù)據(jù)提取降雨侵蝕力中系統(tǒng)誤差的前提。

確定同一場降雨在不同時間分辨率下降雨強度之間關系的途徑有2種:1)通過研究區(qū)氣象站大量的實際降雨過程數(shù)據(jù),統(tǒng)計推算出研究區(qū)不同類型降雨的降雨歷時、降雨量和降雨強度之間的關系式；2)通過研究區(qū)域的暴雨強度公式并結合設計雨型,進行研究區(qū)不同典型降雨過程的時程分配,以這些時程分配統(tǒng)計出不同時間分辨率下不同類型降雨的降雨歷時、降雨量和降雨強度之間的關系式[21][22]3-5[23]14-17。 顯然,在缺少精細降雨資料的情況下,第二種途徑更為可行。根據(jù)得到不同時間分辨率下降強度的關系,計算出對于同一場降雨不同時間分辨率下降雨侵蝕力之間的關系:

式中:μ為降雨侵蝕力差異系數(shù)；R30和R180分別為30 min和180 min時間分辨率降雨資料計算的次降雨侵蝕力,MJ/(mm·hm2)。

2.3 TRMM提取降雨侵蝕力校正方法

基于以上分析,對于1.2節(jié)中使用TRMM數(shù)據(jù)計算大中尺度區(qū)域的降雨侵蝕力的方法,提出了以下校正方法:

1)建立不同歷時和不同降雨量的典型降雨過程降雨侵蝕力之間的關系表(μ值表)；

2)對于實際降雨,根據(jù)該場降雨的歷時和降雨量,在步驟1得到的差異值表中找出對應的降雨侵蝕力差異值μ；

3)利用該差異值μ乘以TRMM數(shù)據(jù)計算的初步結果,得到每一場次降雨降雨侵蝕力的修正值R′s,

4)通過累加每日、每月、每年內的對應降雨場次的修正值,得到修正后日、月、年的降雨侵蝕力。

得到修正后的降雨侵蝕力之后,根據(jù)1.3節(jié)所述利用協(xié)同克里格插值計算降雨侵蝕力分布。

3 應用案例

3.1 研究區(qū)概況

薊運河流域上游地區(qū)位于河海流域北部,東臨灤河,西臨潮白河,地跨北京、天津及河北3省市, E116°50′～118°15′,N39°45′～40°45′之間(圖1)。薊運河主要支流有泃河、州河,泃、州兩支于天津市寶坻區(qū)九王莊匯流后稱薊運河。九王莊以上地區(qū)為薊運河的主要匯水區(qū)域,即為薊運河上游地區(qū)。該區(qū)域銜接燕山山脈與華北平原,其北部地區(qū)為起伏的中低山、丘陵地形,南部地區(qū)以平原區(qū)為主。該地區(qū)屬于溫帶半濕潤大陸性季風氣候,季風氣候顯著。

3.2 數(shù)據(jù)來源

根據(jù)中國氣象數(shù)據(jù)共享服務網年降水量數(shù)據(jù),薊運河上游地區(qū)所在的北京市、天津市和承德市2003—2013年的11年間的降水量相關性很強,且各市2008年的降雨量相當于其11年平均值(圖2)；因此選取2008年為薊運河上游地區(qū)的典型降雨年份,并以此計算薊運河上游地區(qū)降雨侵蝕力分布特征和分析計算結果的精度及誤差。

圖1 研究區(qū)域及區(qū)內主要氣象站點分布示意圖Fig.1 Location of the study area with distribution of meteorological stations

圖2 天津市、北京市和承德市近10年間年降雨量分布Fig.2 Annual precipitation distribution of Tianjin,Beijing and Chengde in the recent decade

研究區(qū)降雨數(shù)據(jù)采用2008年的TRMM 3B42數(shù)據(jù)(數(shù)據(jù)獲取見1.2節(jié)),在ArcGIS中裁剪范圍在E116°26′57″～118°8′10″,N39°16′35″～40°47′33″的數(shù)據(jù)進行研究。

研究區(qū)數(shù)字高程(DEM)數(shù)據(jù)來自日本METI與美國NASA于2011年聯(lián)合開發(fā)的ASTER GDEM數(shù)據(jù)(可在http:∥gdem.ersdac.jspacesystems.or.jp/ index.jsp注冊下載)。選用DEM數(shù)據(jù)格式為Geo-TIFF,圖像的空間分辨率為30 m。

3.3 薊運河上游地區(qū)μ值計算

因缺少研究區(qū)以往精細降雨資料,采用設計典型降雨過程去計算薊運河上游地區(qū)μ值表。典型降雨過程的設計需要確定雨型和暴雨強度公式。

國內外常用的雨型有芝加哥雨型、模式雨型、Huff雨型、三角形雨型、SCS雨型等等。各種雨型之間存在差異,目前還沒有一種公認的雨型作為設計的依據(jù)[23]46-57。本研究區(qū)內,缺乏唐山及承德雨型的相關研究,而天津市區(qū)用Haff雨型的研究比較多,但天津市區(qū)屬于海洋性氣候,降雨過程與內陸有別,故不采用天津市的雨型公式。北京市降雨雨型多用芝加哥雨型或在其基礎上改進的雨型,并結合北京市暴雨公式進行相關設計研究[21][22]20-21,筆者選用北京市用過的芝加哥雨型。

通過薊運河上游地區(qū)暴雨強度公式結合常用的芝加哥雨型計算μ值表,但缺少薊運河上游地區(qū)暴雨強度公式的研究和記載；而天津市相關暴雨強度公式是根據(jù)天津市區(qū)氣象站資料求得[24]。天津市區(qū)臨海,其降雨量受海洋影響較大,而薊運河上游地區(qū)距其較遠,且在內陸,地形地貌差距較大,選用天津市暴雨強度公式不合適；因此采用緯度相近、地理范圍相互重疊、地形地貌相似的北京市第二暴雨分區(qū)的暴雨強度公式作為替代。筆者所用北京市第二區(qū)設計暴雨強度公式[21][22]3為

式中:i為平均降雨強度,mm/min；td為降雨歷時, min；P為設計降雨重現(xiàn)期,a；在選定總降雨量下,迭代試算出P值。

降雨歷程的分配使用芝加哥雨型,雨峰系數(shù)r取0.25,則峰前時間序列i(tb)和峰后時間序列i (tc)分別用如下公式[21,25]計算:

為了簡化運算,在校正降雨侵蝕力時,不考慮每場降雨的降雨量大小不同,只考慮侵蝕性降雨的降雨歷時變化。選取降雨量為25 mm(國家氣象局規(guī)定的中雨和大雨的分界線),歷時分別為0.5、1、2、3、4、5、6 h的7場降雨,利用式(8)～(11)在Matlab 2010b中合成降雨過程,模擬的步長為1 min。

對上述擬合所得不同降雨過程曲線進行統(tǒng)計計算,以降雨開始時刻為統(tǒng)計起點,得到30 min與180 min分辨率下最大、平均降雨強度及對應的降雨侵蝕力(表1)。

表1 不同降雨過程特征物理量及μ值Tab.1 Characteristics of physical quantities including μ value in different rainfall processes

平均降雨強度隨時間分辨率的差異隨降雨歷時而不同(表1),但以此計算降雨動能E相差只有16.5%左右,這相比于最大降雨強度隨時間分辨率的差異較小,此結果與前人研究結果[26]一致。最大降雨強度隨時間分辨率的差異均在135%以上,即I30≈2.35I180,這表明最大降雨強度受時間分辨率的影響程度與降雨歷時的關系不緊密。同時證明Zhu Qiang等[11]和 A.Vrieling等[10,16]研究時簡單認為I30=I180是不對的。

TRMM數(shù)據(jù)只能反映出每3 h內是否發(fā)生降雨,降雨持續(xù)了幾個3 h段,并能得出實際降雨的歷時。根據(jù)李建等[27]研究,20世紀90年代以后北京市降雨以6 h以下降雨為主,且貢獻最大的為2～4 h降雨。并且筆者通過TRMM數(shù)據(jù)發(fā)現(xiàn)＞6 h的降雨非常少,只占到6.7%左右,沒有連續(xù)＞9 h的降雨。另外,根據(jù)表1,隨著降雨歷時的增加,μ值變化幅度較小。因此我們采用了2、4、6 h的降雨代表＜3 h、介于3～6 h的降雨、＞6 h的降雨。即校正時,歷時≤3 h的降雨,其降雨侵蝕力乘以μ2h(2.85)；降雨介于3～6 h的降雨,其降雨侵蝕力乘以 μ4h(2.75)；歷時＞6 h的降雨,其降雨侵蝕力乘以μ6h(2.62)。本文中所模擬的2、4、6 h的降雨時程過程如圖3所示。

3.4 結果與分析

表2 薊運河上游地區(qū)7個縣市的主氣象站點所在網格的降雨侵蝕力校正前后對比Tab.2 Rainfall erosivity on grids with main meteorological sites on seven counties of Jiyun River upriver before and after the calibration

3.4.1 降雨侵蝕力校正前后數(shù)據(jù)對比 表2說明,對每場降雨的降雨侵蝕力經過校正之后,不同網格處累計得到2008年降雨侵蝕力的值均提高了2倍以上,使研究區(qū)2008年的降雨侵蝕力的量級與多年平均降雨侵蝕力[28-29]一致；故這樣的校正顯著提高了TRMM數(shù)據(jù)使用時的精度。

圖3 設計25 mm降雨量不同歷時降雨過程Fig.3 Rainfall process in different designed durations at precipitation of 25 mm

3.4.2 降雨侵蝕力分布趨勢 圖4表明較大的侵蝕力出現(xiàn)在該區(qū)域南部地區(qū),中部地區(qū)降雨侵蝕力居中,北部地區(qū)較小,即該地區(qū)降雨侵蝕力自東南到西北逐漸減弱的趨勢。圖5中可以直觀地看出,東西、南北方向變化趨勢近似。自東向西,數(shù)據(jù)點分布較為分散,東西方向擬合趨勢線有微弱的先升再降的趨勢；而南北方向擬合趨勢線數(shù)據(jù)點自南向北有微弱的減弱趨勢。以上2方向上的趨勢與薊運河上游地區(qū)地形相符合,北部山勢影響山前平原降雨比較多。

3.4.3 降雨侵蝕力插值結果及分析 協(xié)同克里格方法是根據(jù)半變異函數(shù)分析所提供的變量空間自相關程度的信息來進行插值,插值準確程度與半變異函數(shù)擬合的好壞直接相關[19]。表3示出,指數(shù)模型在平均值、標準平均值、標準均方根3個指標上是誤差最低。在均方根指標上,指數(shù)模型的結果與該指標最低值差異不大。只有在平均標準誤差上,指數(shù)模型的結果高于其他模型。綜合其他2個模型的誤差結果來看,指數(shù)模型模擬效果最好,故本研究采用指數(shù)模型進行后面的插值參數(shù)改進。

圖4 2008年降雨侵蝕力自然間斷點分類符號圖Fig.4 Classification symbols of natural discontinuities of rainfall erosivity in 2008

塊金值與基臺值的比值C0/(C0+C)稱為塊金效應,它表示隨機部分引起的空間變異占系統(tǒng)總變異的比例[30]。如果該比值較低,表明結構性因素引起的空間變異起主要作用；如果該值較高,說明隨機部分引起的空間變異起主要作用。通常認為:塊金效應＜25%,空間相關性強；25%～75%之間,空間相關性中等；＞75%,空間相關性弱。本文采用指數(shù)模型計算的塊金值為0,基臺值為72 241,塊金效應＜25%,表明TRMM像素值具有強烈的空間相關性,其結構性因素引起的空間變異起主要作用。

圖5 2008年降雨侵蝕力趨勢面分析Fig.5 Trend surface analysis of rainfall erosivity in 2008

依據(jù)擬合的半變異函數(shù)模型協(xié)同DEM進行克里格插值。研究區(qū)年降雨侵蝕力分布約為900～2 900 MJ·mm/(hm2·h)。圖6表明,該區(qū)域降雨侵蝕力東南到西北逐漸減弱,南部地區(qū)大于北部地區(qū),該變化趨勢和章文波等[28]、劉斌濤等[29]研究結果一致。

表3 不同半變異模型的誤差評價結果Tab.3 Error analysis of different fitting models of semi-variance function

圖6 薊運河上游各縣域降雨侵蝕力分布圖Fig.6 Rainfall erosivity distribution on the upriver counties of Jiyun River

4 結論與討論

以典型降雨年份2008年薊運河上游地區(qū)的降雨侵蝕力計算為研究案例,結果顯示校正后的降雨侵蝕力接近了多年平均降雨侵蝕力的量級及分布。這說明通過確定不同降雨過程中不同時間分辨率下降雨侵蝕力差異的μ值表,并以此對每一場降雨進行修正,進而減少TRMM數(shù)據(jù)計算降雨侵蝕力是系統(tǒng)誤差的方法,是一種切實可行的、有效提高TRMM應用結果精度的方法。

受到TRMM數(shù)據(jù)時空分辨率的影響,該方法只能用于大中尺度區(qū)域降雨侵蝕力的研究,不能用于小尺度區(qū)域的研究。該方法并不限制如何選擇降雨雨型和下墊面條件的影響。降雨雨型不同選擇,會使影響校正結果的精度,但是這種影響遠比未校正之前的誤差小。下墊面條件引起的誤差屬于隨機性誤差,在現(xiàn)有研究的基礎上,暫時無法估計下墊面條件的影響規(guī)律和引起誤差大小。

確定不同降水過程對應的μ值是應用該方法的前提和基礎。這一步驟需要以實際降雨過程數(shù)據(jù)進行統(tǒng)計分析為基礎,或者以更多典型降水過程的模擬數(shù)據(jù)為基礎。由于各地降雨量不同,降雨過程隨機變化,應用該方法是需要因地制宜,確定出適合當?shù)氐摩讨当?。當然在誤差允許的條件下,也可以制訂各氣候區(qū)或降雨類型區(qū)統(tǒng)一的μ值表,方便以后隨時調取使用。從本文案例的校正過程來看,影響校正結果效果的因素有雨型確定的方法、雨型種類、不同降雨歷時的降雨所占的比例、不同降雨歷時及不同降雨量的降雨對應的雨型種類等。追求更高精度的校正效果時,須對這些問題進行深入的研究,這也是今后進一步研究和改進完善該校正方法途徑。

[1] 王占禮.中國土壤侵蝕影響因素及其危害分析[J].農業(yè)工程學報,2000,16(4):32- 36

[2] Wischmeier W H,Smith D D.Predicting rainfall erosion losses:a guide to conservation planning[J].Agriculture Handbook 537,1978:5- 8,50- 51

[3] Renard K G,Foster G R,Weesies G A,et al.Predicting soil erosion by water:aguide to conservation planning with the revised universal soil loss equation(RUSLE) [J].Agriculture Handbook 703,1997:22- 38

[4] 葉芝苗,劉寶元,章大波,等.北京市降雨侵蝕力及其空間分布[J].中國水土保持科學,2003,1(1): 16- 20

[5] 章文波,付金生.不同類型雨量資料估算降雨侵蝕力[J].資源科學,2003,25(1):35- 41

[6] 卜兆宏,董勤瑞,周伏建,等.降雨侵蝕力因子新算法的初步研究[J].土壤學報,1992,29(4):408- 417

[7] Yu Baofu,Rosewell C J.A robust estimator of the R-factor for the Universal Soil Loss Equation[J].Transactions of the ASAE,1996,39(2):559- 561

[8] 章文波,謝云,劉寶元.用雨量和雨強計算次降雨侵蝕力[J].地理研究,2002,21(3):384- 390

[9] Fan Jianrong,Chen Yang,Yan Dong,et al.Characteristics of rainfall erosivity based on tropical rainfall measuring mission data in Tibet,China[J].Journal of Mountain Science,2013,10(6):1008- 1017

[10]Vrieling A,Hoedjes J C B,van der Velde M.Towards large-scale monitoring of soil erosion in Africa:Accounting for the dynamics of rainfall erosivity[J].Global and Planetary Change,2014,115:33- 43

[11]Zhu Qiang,Chen Xiuwan,Fan Qixiang,et al.A new procedure to estimate the rainfall erosivity factor based on Tropical Rainfall Measuring Mission(TRMM)data[J]. Science China Technological Sciences,2011,54(9): 2437- 2445

[12]Kummerow C,Barnes W,Kozu T,et al.The tropical rainfall measuring mission(TRMM)sensor package[J]. Journal of atmospheric and oceanic technology,1998,15 (3):809- 817

[13]Vrieling A,Sterk G,de Jong S M.Satellite-based estimation of rainfall erosivity for Africa[J].Journal of hydrology,2010,395(3):235- 241

[14]Mannaerts C M,Saavedra C P.Regional scale erosion modeling and monitoring using remotely sensed data: some spatial data scale issues[C]∥Proceedings of the international symposium on 25 years of assessment of erosion.Ghent,Belgium:Ghent University,2003:433- 440

[15]Pham T Y.Soil erosion risk modeling within upland landscapes in Vietnam using remotely sensed data and the RUSLE model[D].Nova Scotia:Dalhousie University, 2008:79- 84

[16]Vrieling A,de Jong S M,Sterk G,et al.Timing of erosion and satellite data:a multi-resolution approach to soil erosion risk mapping[J].International Journal of Applied Earth Observation and Geoinformation,2008,10, 267- 281

[17]Xie Yun,Liu Baoyuan,Nearing M A.Practical thresholds for separating erosive and non-erosive storms[J]. Transactions of the ASAE,2002,45(6):1843- 1847

[18]Yin Shuiqing,Xie Yun,Nearing M A,et al.Estimation of rainfall erosivity using 5-to 60-minute fixed-interval rainfall data from China[J].Catena,2007,70(3): 306- 312

[19]馮益明,雷相東,陸元昌.應用空間統(tǒng)計學理論解譯遙感影像信息“缺失”區(qū)[J].遙感學報,2004,8(4): 317- 322

[20]傅抱璞.地形和海拔高度對降水的影響[J].地理學報,1992,59(4):302- 314

[21]張大偉,趙冬泉,陳吉寧,等.芝加哥降雨過程線模型在排水系統(tǒng)模擬中的應用[J].給水排水,2008,34 (增刊2):354- 357

[22]北京市城市規(guī)劃設計研究院.城市雨水系統(tǒng)規(guī)劃設計暴雨徑流計算標準:DB11/T 969—2013[S].北京:中國計劃出版社,2013:3- 5,20- 21

[23]范澤華.天津市降雨趨勢分析及設計暴雨研究[D].天津:天津大學,2012:14- 17,46- 57

[24]黃津輝,向文艷,戶超,等.天津市設計暴雨方法比較及公式修正[J].天津大學學報,2013,46(4): 354- 360

[25]Keifer C J,Chu H H.Synthetic storm pattern for drainage design[J].Journal of the hydraulics division,1957,83 (4):1- 25

[26]殷水清,謝云,王春剛.用小時降雨資料估算降雨侵蝕力的方法[J].地理研究,2007,26(3):541- 547

[27]李建,宇如聰,王建捷.北京市夏季降水的日變化特征[J].科學通報,2008,53(7):829- 832

[28]章文波,謝云,劉寶元.中國降雨侵蝕力空間變化特征[J].山地學報,2003,21(1):33- 40

[29]劉斌濤,陶和平,宋春風,等.1960—2009年中國降雨侵蝕力的時空變化趨勢[J].地理研究,2013,32(2): 245- 256

[30]Schabenberger O,Pierce F J.Contemporary statistical models for the plant and soil sciences[M].Boca Ration: CRC press,2002:563- 599

(責任編輯:宋如華)

Calibration of rainfall erosivity calculation based on TRMM data: A case study of the upriver basin of Jiyun River,North China

Jin Qiutong1,2,Shi Mingchang1,2,3,Zhang Jutao1,2,Wang Shan4,Hu Ying3
(1.School of Soil and Water Conservation,Beijing Forestry University,100083,Beijing,China；2.Key Lab.of Soil&Water Conservation and Desertification Combating,Beijing Forestry University,Ministry of Education,100083,Beijing,China；3.Beijing Datum Science and Technology Development Co.,Ltd.,100084,Beijing,China；4.AutoNavi Software Co.,Ltd.,102200,Beijing,China)

The low temporal resolution of TRMM data could cause systematic errors in the rainfall erosivity calculation within the specified time period.For this reason,we propose a method of calibration in order to reduce the errors in this paper.Firstly,a relationship table of rainfall erosivity at different events and duration with typical rainfall process is established with historical rainfall data or numerical simulation；then,the corresponding difference is calculated from the above table based on the precipitation and duration of each rainfall in the study area；afterwards,the adjusted value of rainfall erosivity in each rainfall is calculated by multiplying the figure found in the pervious step by the original value of rainfall erosivity correspondingly；finally,the calibrated value of daily,monthly or yearly rainfall erosivity isobtained by adding daily, monthly oryearly adjusted value ofrainfallerosivity correspondingly.To verify the effectiveness of this calibration method,the typical rainfall process in the study area,the upriver regions of Jiyun River,was simulated within Matlab 2010b to obtain therelationship table,and then the rainfall erosivity in the study area was calculated using TRMM 3B42 data in 2008 with the spatial resolution of 0.25°×0.25°and the time resolution of 3 h.The calculation was conducted in the ENVI-IDL Language program and calibrated by the relationship table.The results showed that:the rainfall erosivity in the upriver of Jiyun River ranged from 900 to 2 900 MJ·mm/(hm2· h)with the highest values in central area,followed by southern and northern parts,which was consistent with the local average annual values and distribution of rainfall erosivity.The results confirm that this method can effectively improve the precision of TRMM data's application.Due to the limitation of spatial resolution of TRMM data,this calibration method is only applied in the rainfall erosivity calculation at large and middle scales.In order to improve the precision of this calibration method,the typical rainfall process,rain type design and distribution of rainfall duration should be detailed.

rainfall erosivity；TRMM；systematic error；calibration method；Jiyun River

S157.1

A

1672-3007(2015)04-0094-09

2014- 11- 15

2015- 05- 03

項目名稱:林業(yè)公益性行業(yè)科研專項經費項目課題“基于生態(tài)安全的水土保持措施空間配置技術”(201404209)

靳秋桐(1990—),男,碩士研究生。主要研究方向:3S技術集成與系統(tǒng)開發(fā)。E-mail:jinqiutong@bjfu.edu.cn

?通信作者簡介:史明昌(1969—),男,博士,教授。主要研究方向:地理信息科學。E-mail:shimc@bjfu.edu.cn