賈 薛,唐彥君
(1.河海大學地理信息科學與工程研究所,南京 211100;2. 寧波弘泰水利信息科技有限公司,浙江 寧波 315192)
準確地獲取降雨的時空分布信息是進行流域水資源管理,特別是洪水預報、徑流模擬、水質預測的重要環節[1],然而由于經濟和人力原因,雨量觀測站點數量極其有限且其分布極為不均勻[2],因此,將空間分布不均勻的降水數據精確插值成規則分布的空間數據集對于解決水文關鍵問題具有重要的實際應用價值[3]。對于降水插值,國內外學者進行了大量研究[4, 5]。由于降水過程的復雜性,單純考慮降水數據本身的插值可能有較大誤差,因此加入高程因素可能對插值精度的提高有所幫助[6, 7]。岳文澤[8]、曾懷恩[9]等進行了基于地統計方法的空間插值研究認為研究區的地理特征對降水插值有較大影響;趙登忠[10]等進行了基于DEM的地理要素PRISM空間內插研究,認為要充分的將地理要素考慮到降水插值中來。Daly[11]等認為對于山區或者降水站點不是很密集的地區, 反距離加權法(IDW)有助于提高預測數據的精度;而IDW沒有考慮到高程等因素,有很大的缺點。本文在綜合國內外研究的基礎上,以瀏陽河流域為例,選取反距離權重法、克里金法,對流域暴雨進行空間分布的研究,然后分別采用加入高程因素的方法,改進現有2種插值方法,并對比插值精度,分析各種插值方法對瀏陽河流域的適用性。
瀏陽河是湖南省最大河流——湘江的最主要支流,地處湖南省東部偏北,全長234.8 km,流域總面積4 665 km2。流域屬于中亞熱帶季風濕潤氣候區,夏季炎熱多雨,冬季寒冷濕潤。多年平均降雨量1 400~1 800 mm,各地年平均氣溫為16~18 ℃,其中夏季是降雨最多的季節,一年中絕大部分暴雨也分布于夏季。
圖1是本研究的降雨站點分布圖。根據瀏陽河上中下游的劃分,流域共有10個雨量站點在上游(張坊、光明、寒婆坳等),4個雨量站在中游(雙江口、淮川、幸福、目蓮渡),4個在下游(江背、黃興、同升湖、樹木嶺)。圖2為流域30 m分辨率的DEM數據經過分層設色后的結果。DEM數據總體反應了流域的高程的一般特征,整個流域海拔從東南向西北呈現出遞減趨勢,東北部主要為山地丘陵,地勢險峻,海拔較高,最大高度超過1 500 m,中部主要為低山丘陵,大部分在1 000 m以下,西南部主要為沖擊平原地貌,地勢較為平坦,海拔較低大部分在500 m以下。

圖1 研究區雨量站點分布Fig.1 Distribution of rainfall station in research area

圖2 瀏陽河流域DEMFig.2 DEM of Liuyang River basin
流域降水數據是本研究的重要數據。根據中國氣象上暴雨的定義,每小時降雨量16 mm以上、或連續12 h降雨量30 mm以上、24 h降雨量50 mm以上的雨稱為“暴雨”,利用VBA嵌套語句分別檢測出研究區18個站點2013-2014年的逐日逐時降雨量數據中滿足上述3種定義的降雨,并篩選出2場具有代表性的場次暴雨(從降雨開始至降水結束,時間跨度可能超過24 h)作為本研究的原始數據(A場次2013年6月27左右;B場次2014年5月24日左右),選擇這2個場次主要是考慮到該時段內區域降雨相對比較集中,盡量避免因無降雨而反映不出降雨的空間變異,如圖3所示。

圖3 瀏陽河流域場次暴雨與高程數據(以站點高程升序排列)Fig.3 Data of Liuyang River basin
反距離權重法是一種基于相近相似原理的確定性插值方法,即2個物體離得越近就越相似,反之,離得越遠則相似性越小[10]。其表達式為:
(1)
式中:P為待估數值;pi為i點的實測數值;di為待估值點距離實測數值點i的距離;n為計算待估值點的實測數據點的個數。
克里金方法是一種建立在地質統計學基礎上的插值方法。該方法最早由法國地理學家Matheron和南非礦山工程師Krige[12]提出并用于礦山勘探。這種方法充分吸收了地理統計的思想,認為任何在空間連續變化的屬性是非常不規則的,不能用簡單的平滑函數進行描述,只能用隨機表面函數進行模擬[13]??死锝鸱椒ǖ年P鍵在于確定權重系數,該方法在插值過程中根據某種優化準則函數來動態地決定變量的數值,從而使內插函數處于最佳狀態[14]。其表達式為:
(2)
式中:Z(x0)為x0處的估計值;Z(xi)為xi處的觀測值;λi為克里金權重系數;n為觀測點個數。
考慮高程的反距離權重法是對反距離權重方法的改進,他不僅考慮待插值點與樣本點的水平距離,還考慮到他們的垂直高程,即樣本點距離插值點越近,與待插值點高程差越小,對插值點的影響越大,相應的權重系數也越大。其表達式為:

(3)


(4)
式中:WD、WZ分別為距離和高程對降雨的影響權重,它們之和為1;Pi為觀測站i的降雨量;di為待估點與觀測站i之間的距離,km;zi為待估點與觀測站i之間的高程差,m;zi+10是為了避免待插值點與某些實測點高程很接近,或是相同,而使w(z)i過大的情況。
在加入一個輔助變量的協克里金中,對任意一待插位置x0,其估值可用以下公式求得:
(5)
式中:Z1是主變量;Z2是輔助變量;n為待插值點x0鄰域內的變量Z1的實測點個數;m為待插值點x0鄰域內的變量Z2的實測點個數;ai是權重系數,表示空間樣本點xi處的觀測值Z1(xi)對估計值Z(x0)的貢獻程度;bj是權重系數,表示空間樣本點xj處的觀測值Z(xj)對估計值Z(x0)的貢獻程度。
在無偏條件下,用拉格朗日乘數法得到協克里金方程組:
(6)

由此可見,求解協克里金插值的問題的關鍵是互變異函數的求解。
考慮到雨量站點的數量情況,采用交叉驗證[15]的方法來驗證插值的結果。關于插值效果的標準選擇,本研究選擇平均絕對誤差(MAE)、平均相對誤差(MRE)這2個指標。MAE估量模擬值可能的誤差范圍,MRE則反映模擬值相對于測量值的準確度[3]。其表達式為:
(8)

使用4種方法分別對瀏陽河流域18個站點,A、B 2個場次的暴雨雨量進行空間插值。由圖4可以看出不考慮高程的反距離權重法和普通克里金插值結果較為平滑,反距離權重法主要依賴空間距離,距離的遠近決定了權重的大小,因此插值結果表現出以實測點為中心,向外依次變化的趨勢;而克里金插值考慮距離和空間數據結構,因此插值效果與反距離權重法有一些差別。由于本研究插值點數據較少,插值易受極值點的影響,“牛眼”現象較為明顯,結果中某些值域區間的斑塊菱角分明,變化較為劇烈。考慮高程的反距離權重法插值結果相對于不考慮高程的反距離權重法,結果更為不平滑。在按每個實測點為中心向外依次變化的同時,加入了高程信息的影響,使得降雨分布發生了一定變化,預測結果反映出一定的高程信息,變得更加粗糙??紤]高程的協克里金方法的插值結果趨勢與克里金法基本一致,協克里金法運用高程信息作為輔助變量,插值結果在局部地區與克里金法稍有差異。


圖4 瀏陽河流域場次暴雨空間插值結果 Fig.4 The results of rainfall spatial interpolation of Liuyang River basin
由表4可以看出,對于瀏陽河流域A、B 2個場次的暴雨雨量插值,加入高程因素后,插值結果均有所提升,4種插值方法MAE的排序為IDW>OK>CK>IDEW,MRE的排序與MAE相同。其中,不考慮高程的反距離權重和普通克里金插值精度相差不大,而協克里金法由于考慮了高程因素,精度略有提高。4種方法中,考慮高程的反距離權重取得了略好的效果。

表4 精度評價Tab.4 Accuracy assessment
經過對瀏陽河流域A、B 2場暴雨的空間插值結果進行分析,并用交叉驗證的方法,利用平均相對誤差、平均絕對誤差2種精度評價模型對插值結果進行統一評價,結果表明:
(1)在山丘區流域,由于地形復雜性的影響,加入高程影響因素后,降雨插值結果的精度有所提高,且較為真實地反映了瀏陽河流域降雨的空間分布,其中考慮高程的反距離權重方法略好于考慮高程的協克里金方法,在瀏陽河流域的降水插值中效果較好。
(2)當研究區采樣點(雨量站)個數較少時,對于克里金插值法來說,要想獲得其中所必須的半變異函數是比較困難的,對于反距離權重法插值來說,插值結果“牛眼”現象較為明顯,這在一定程度上也會影響插值效果。因此,保證一定密度的站網分布很有必要。
(3)由于瀏陽河流域上下游高程差異較大,高程是影響降雨的重要因素之一。除了高程, 降水的空間分布影響因素還有很多, 如降雨雨型、站點所處的坡度及坡向等其他地理環境要素, 將這些可能影響降雨的因素合理引進到降雨插值中來將進一步提高降水插值的精度。
□
[1] 于 野, 王 闖, 王 錚. 地理信息系統支持下的降雨時空統計分析[J]. 測繪通報, 2003,(2):44-46.
[2] Willmott C J, Robeson S M, Feddema J J. Estimating continental and terrestrial precipitation averages from rain-gauge networks[J]. International Journal of Climatology, 1994,14(4):403-414.
[3] 江善虎, 任立良, 雍 斌, 等. 老哈河流域降水的空間插值方法比較[J]. 干旱區資源與環境, 2010,(1):80-84.
[4] Houghton J T. Climate change 1995: the science of climate change[J]. Climatic Change, 1996,67(5 520):1 261-1 261.
[5] Lam N S-N. Spatial interpolation methods: a review[J]. American Cartographer, 1983,10(2):129-150.
[6] Hevesi J A, Istok J D, Flint A L. Precipitation estimation in mountainous terrain using multivariate geostatistics, part I: structural analysis[J]. Japplmeteor, 1992,31(7):661-676.
[7] Goovaerts P. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall[J]. Journal of Hydrology, 2000,228(1):113-129.
[8] 岳文澤, 徐建華, 徐麗華. 基于地統計方法的氣候要素空間插值研究[J]. 高原氣象, 2005,(6):974-980.
[9] 曾懷恩, 黃聲享. 基于Kriging方法的空間數據插值研究[J]. 測繪工程, 2007,(5):5-8,13.
[10] 趙登忠, 張萬昌, 劉三超. 基于DEM的地理要素PRISM空間內插研究[J]. 地理科學, 2004,(2):205-211.
[11] Daly C. Guidelines for assessing the suitability of spatial climate data sets[J]. International Journal of Climatology, 2006,26(6):707-721.
[12] 侯景儒. 地質統計學的理論與方法[M]. 北京: 地質出版社, 1990:69-78.
[13] 李俊曉, 李朝奎, 殷智慧. 基于ArcGIS的克里金插值方法及其應用[J]. 測繪通報, 2013,(9):87-90,97.
[14] 劉登偉, 封志明, 楊艷昭. 海河流域降水空間插值方法的選取[J]. 地球信息科學, 2006,(4):75-79,83.
[15] 方書敏, 錢正堂, 李遠平. 甘肅省降水的空間內插方法比較[J]. 干旱區資源與環境, 2005,(3):47-50.