雨量站網(wǎng)分布對雨量插值算法及徑流響應(yīng)的影響

汪青靜，楊欣玥，陳 華，許崇育，3，曾 強，徐 堅

（1.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2.河海大學(xué)水文與水資源學(xué)院，南京 210098；3.奧斯陸大學(xué)地學(xué)系，挪威奧斯陸 0317）

1 研究背景

降雨的空間分布對于氣象學(xué)、氣候?qū)W、水文學(xué)、地質(zhì)分析、環(huán)境監(jiān)控等都必不可少。現(xiàn)有的降雨數(shù)據(jù)主要通過雨量站點進(jìn)行收集，雨量站點能比較精確地測量點降雨，但難以獲得大范圍的降水分布，時空連續(xù)的降雨分布，需要根據(jù)已知雨量站點的資料通過空間插值得到。研究表明，降雨的空間分布不均勻性考慮得越充分，其水文過程模擬精度越高［1］。其中不同插值算法、不同雨量站密度會影響降雨空間分布的不均勻性［2-4］。

現(xiàn)有的空間插值方法可以被概括為確定性和地理統(tǒng)計的插值方法［5］，在比較不同插值算法時，大多選擇交叉驗證來比較其降雨插值結(jié)果。Driks等［6］用13個雨量站的降雨數(shù)據(jù)比較了泰森多邊形法、反距離權(quán)重法、區(qū)域平均值和克里金法，他們得出計算復(fù)雜的克里金法在結(jié)果上并沒有明顯的提高，建議在雨量站密集的地方使用反距離權(quán)重法進(jìn)行空間插值。高歌等［7］對反距離權(quán)重法和普通克里金法的日降雨插值結(jié)果進(jìn)行比較，發(fā)現(xiàn)普通克里金法的插值結(jié)果略好于反距離權(quán)重法。Ly等［8］在分析不同插值算法的日降雨時得到，考慮了高程的具有外部漂移的克里金法和協(xié)同克里金法對于日降雨的插值結(jié)果并沒有提高，普通克里金法和反距離權(quán)重法被認(rèn)為是最好的日降雨插值方法。也有少量研究將插值結(jié)果代入水文模型中，根據(jù)徑流模擬結(jié)果比較不同的插值算法。Ruelland等［9］在比較不同空間插值算法對集總式和半分布式模型的敏感性時得到，反距離權(quán)重法對于水文模型的模擬較好。綜上分析，研究區(qū)資料不同，評價指標(biāo)不同，不能統(tǒng)一得出某種最好的空間插值方法。

通常認(rèn)為，雨量站數(shù)目越多，對降雨的估計誤差就越小。然而，由于地形和經(jīng)濟(jì)因素的限制，很多地區(qū)的雨量站數(shù)目和分布都受到了限制。Ruelland等［9］指出在雨量站數(shù)目太少時，其插值結(jié)果對水文模型的輸入誤差很大，建議使用雨量站數(shù)目較多的區(qū)域進(jìn)行研究。Dong等［10］在研究水文模型對徑流模擬最合適的雨量站數(shù)目時得到，當(dāng)雨量站數(shù)目到達(dá)到一個閾值后，隨著雨量站數(shù)目的增加，其徑流模擬的效果不再增加。Anctil等［11］指出，相比使用所有雨量站的計算結(jié)果，合理布設(shè)的少量雨量站數(shù)目也可以得到較好的計算結(jié)果。Xu等［12］研究表明，在雨量站數(shù)目較少但空間分布合理時，其水文模型的表現(xiàn)會更好。因此，雨量站網(wǎng)的結(jié)構(gòu)不僅取決于雨量站數(shù)目，其分布對于降雨和徑流的估計也至關(guān)重要。然而現(xiàn)有研究主要分析不同雨量站網(wǎng)情況下的一種或兩種插值算法，無法給出不同插值算法在雨量站數(shù)目相同時的特定優(yōu)勢。

為比較雨量站網(wǎng)對不同空間插值算法雨量及其徑流響應(yīng)，本文將從以下3個方面展開研究：①在不同雨量站網(wǎng)密度下，不同空間插值算法的雨量誤差分析；②在不同雨量站網(wǎng)密度下，不同空間插值算法的徑流響應(yīng)分析；③在雨量站數(shù)目相同，但是空間分布不同時，不同空間插值算法的降雨-徑流響應(yīng)分析。

2 研究區(qū)域及資料

瀟水是湘江上游的一級支流，本文的研究區(qū)域為湘江上游雙牌站以上瀟水子流域。該流域位于24°39′N—26°4′N，111°6′E—112°10′E之間，整個子流域面積是10 434 km2。該子流域高程范圍80～1 998 m（圖1），流域的西南到東北為沿河平原，東南和西北為山地，高程超過1 000 m。該流域是亞熱帶季風(fēng)氣候，年均氣溫在17.6℃，年均降雨為1 600 mm，東部山區(qū)降雨較大，可達(dá)1 900 mm以上，北部年均降雨較小，在1 300 mm左右。

圖1 研究區(qū)域高程及氣象站點分布Fig．1 Drainage basin delineation and distribution of meteorological stations in the study area

本文研究區(qū)內(nèi)一共有44個分布均勻的雨量站，1個蒸發(fā)站，1個流域出口流量雙牌站（圖1）。這些站點具備流域2006—2014年完整的日尺度數(shù)據(jù)，均來源于湘江的水文年鑒。

3 主要研究方法

本文先將研究區(qū)域離散成1 km×1 km的網(wǎng)格，使用了4種空間插值算法：泰森多邊形法、反距離權(quán)重法、普通克里金法和協(xié)同克里金法，得到研究區(qū)2006—2014年所有1 km網(wǎng)格的日降雨序列，對比不同插值算法的降雨誤差，再將子流域內(nèi)所有1 km網(wǎng)格雨量取平均值，得到研究區(qū)每日的面均降雨，作為新安江模型的輸入，與實測徑流對比。

3.1 雨量空間插值方法

（1）泰森多邊形法（Thiessen，簡稱 THI法）［13］將雨量站兩兩相連并作連線的中垂線，中垂線相交形成與雨量站數(shù)目相同的多邊形，從而將流域劃分成多個多邊形（本文共形成44個多邊形），每個多邊形里面包含一個雨量站，則該區(qū)域內(nèi)的雨量站均用該雨量站的實測值替代。泰森多邊形法的優(yōu)點是計算簡單并且使用廣泛，缺點是空間插值的結(jié)果在多邊形內(nèi)均勻，在多邊形邊界上突變，這個與實際不符。

（2）反距離權(quán)重法（Inverse Distance Weighted，簡稱IDW法）以插值點與樣本點間的距離為權(quán)重進(jìn)行加權(quán)平均，離插值點越近的樣本賦予的權(quán)重越大［14］，本文的距離權(quán)重指數(shù)為2，選擇的是距離插值點最近的16個雨量站作為樣本點，不足16個樣本點的情況則取其所有值。反距離權(quán)重法的優(yōu)點是簡單易行，但缺點是易受極值的影響，出現(xiàn)“牛眼”現(xiàn)象［15-16］。

（3）普通克里金法（Ordinary Kriging，簡稱 OK法）提供了一個在有限區(qū)域內(nèi)對空間變量進(jìn)行無偏最優(yōu)估計的方法，其計算公式為

式中：Z*是雨量估算值；Zi（i＝1，2，…，n）是實測站點的雨量值；權(quán)重λi是根據(jù)克里金法的無偏性和估計方差最小計算得到，本文選擇指數(shù)模型作為擬合半方差的理論模型［17］。

（4）協(xié)同克里金法（Cokriging，簡稱 CK法）將區(qū)域化變量的無偏最優(yōu)估計方法從單一屬性發(fā)展到兩個或者兩個以上的協(xié)同區(qū)域化屬性［18］?？紤]到降雨和高程的影響，本文將高程作為第2變量考慮在其中，其計算公式為

式中：Z*是估算降雨值；Zi（i＝1，2，…，n）是初始變量，即實測站點的雨量值；Ej（j＝1，2，…，m）是二級變量，即實測站點的高程值，其權(quán)重λi和αj代表初始和二級變量的權(quán)重，也是根據(jù)克里金法的無偏性和估計方差最小得到。與普通克里金法一樣，本文選擇指數(shù)模型作為半變異函數(shù)的理論模型。

3.2 新安江模型

新安江模型［19］是一種概念性降雨徑流模型，本文選用的三水源新安江模型在濕潤及半濕潤地區(qū)得到了廣泛的認(rèn)可和應(yīng)用［12，19-21］。該模型按三水源劃分為地表徑流、壤中流和地下徑流，產(chǎn)流方式為蓄滿產(chǎn)流，蒸發(fā)分為上層、下層和深層，匯流分為坡地匯流和河網(wǎng)匯流2個階段，地面匯流采用納什單位線，壤中流和地下徑流采用的線性水庫。

本文將研究區(qū)的面均雨量站作為雨量輸入驅(qū)動新安江模型，模擬雙牌流量站的徑流過程。因為本文目的是比較雨量站網(wǎng)對不同插值算法的徑流響應(yīng)，所以只計算了率定期2006—2014年的日徑流量。參數(shù)率定采用SCE-UA優(yōu)化算法［22］，以納什效率系數(shù)（NSE）作為目標(biāo)函數(shù)自動搜尋水文模型參數(shù)最優(yōu)解。

3.3 雨量站網(wǎng)密度及分布

為了比較雨量站網(wǎng)密度及分布對不同插值算法的影響，本文選取了6種不同的雨量站網(wǎng)密度（見表1）。其中雨量站最少的情況也滿足了WMO在丘陵地區(qū)建議的最小雨量站密度［23］。對于前6種雨量站網(wǎng)密度，每次都從44個雨量站中隨機(jī)選取100種不同的雨量站網(wǎng)分布。Xu等［12］已經(jīng)證明過，100次隨機(jī)取樣可以代表不同雨量站密度下的區(qū)域降雨變化，因此對于每一種插值算法都要計算601次1 km×1 km網(wǎng)格插值結(jié)果。

表1 不同雨量站網(wǎng)等級的雨量站數(shù)目選取Table 1 Number of rain gauges corresponding to different density levels

對于每一種雨量站密度下的100種不同的空間分布，以空間統(tǒng)計中的最鄰近距離指數(shù)（NNI）作為選取標(biāo)準(zhǔn)［24-25］，NNI常用來刻畫點要素的分布格局是集聚、隨機(jī)還是均勻分布，是根據(jù)每個點與最近鄰點之間的平均距離計算，其計算公式為：

式中：NNI為雨量站分布最鄰近距離系數(shù)；d（NN）為雨量站平均最鄰近距離（m）；d（ran）為期望平均距離（m）；di為觀測點到第i個雨量站的距離（m）；A為流域面積（m2）；n為流域內(nèi)雨量站數(shù)目。一般認(rèn)為，如果NNI＜1，則認(rèn)為雨量站網(wǎng)在空間分布上聚類；如果NNI＞1，則認(rèn)為雨量站網(wǎng)在空間上分布均勻。

3.4 評價指標(biāo)

降雨和徑流是水文模型的兩個重要指標(biāo)，為了比較雨量站網(wǎng)對不同插值算法的降雨徑流影響，本文將從降雨和徑流2個方面對其插值結(jié)果進(jìn)行比較。

評價降雨精度，一般使用均方根誤差（RMSE，Root Mean Square Error），表示插值結(jié)果偏離實際值的大小。本文將從點雨量和面雨量2個方面進(jìn)行分析。點雨量誤差分析主要針對于雨量站點實測值與估算值的誤差計算；對于面雨量，本文假設(shè)44個雨量站插值的面均雨量結(jié)果為基準(zhǔn)值，而由不同的雨量站密度及分布計算的面均雨量結(jié)果作為估算值。基準(zhǔn)值與估算值的RMSE越大，則說明插值結(jié)果與實測降雨誤差大，RMSE越小，則認(rèn)為插值結(jié)果與實測降雨越小，計算公式為

式中：Zest為估測值；Zobs為實測值；n為降雨天數(shù)。

評價徑流時，將不同雨量站密度及分布的面均降雨代入新安江模型中計算，根據(jù)納什效率系數(shù)NSE比較其徑流模擬精度。

4 結(jié)果及分析

4.1 雨量空間分布規(guī)律分析

本文選用了4種經(jīng)常使用的空間插值算法，因為這4種插值算法代表了不同的降雨空間分布，也考慮了不同的降雨屬性。圖2給出了4種插值算法在使用44個雨量站時的年均降雨插值結(jié)果，從圖中可以看到，反距離權(quán)重法、普通克里金法、協(xié)同克里金法的插值結(jié)果相似，泰森多邊形法的降雨插值結(jié)果在空間上分布不連續(xù)。但是不論對于哪一種插值算法，都能從其年均降雨的空間分布中得出：在流域的東邊和西邊，其年均降雨量較大。

4.2 不同雨量站密度下不同插值算法的雨量誤差分析

對于不同插值算法的雨量誤差分析，本文將從點雨量和面雨量2個方面進(jìn)行分析。點雨量誤差分析是計算雨量站點實測值與估算值的RMSE。根據(jù)前人的研究，交叉驗證正被廣泛地使用［26-28］。針對本文使用的4種插值算法，在使用不同雨量站密度和分布插值時，將剩下來的雨量站（以研究區(qū)總共44個雨量站減去使用的插值雨量站）實測數(shù)據(jù)作為基準(zhǔn)值，與4種插值算法所計算的點雨量估算值進(jìn)行比較，使用RMSE作為評價標(biāo)準(zhǔn)，然后對剩下來雨量站的RMSE取均值，得到每種雨量站分布的點雨量插值誤差。圖3是這6種不同雨量站密度100次隨機(jī)分布下4種插值算法所計算RMSE的箱線圖。從圖3中可以看到，對于這4種插值算法，其RMSE都是隨著雨量站數(shù)目的增加而減小，說明隨著雨量站數(shù)目的增加，其降雨的插值誤差也越來越小。同時，反距離權(quán)重法、普通克里金法、協(xié)同克里金法的點雨量誤差范圍及趨勢基本相同，但是明顯的小于泰森多邊形的點雨量誤差范圍，這是受泰森多邊的插值原理所致，說明泰森多邊形的插值結(jié)果與實測值偏離較大。

圖2 4種不同插值算法在2006—2014年的年均降雨插值結(jié)果Fig．2 Interpolated average annual rainfall for 2006-2014 obtained by four spatial interpolation methods

圖3 不同雨量站密度下100次隨機(jī)分布的雨量站點間實測值與估算值的RMSE的箱線圖Fig．3 Box plot of RMSE of observed and estimated daily rainfall in the presence of different densities of rain gauges in 100 stochastic samplings for each number of stations respectively

區(qū)域的面均雨量是新安江模型的一個重要輸入項，而不同插值算法在雨量站密度及分布相同時也會得到不同的面均降雨。圖4給出了這4種不同插值算法的面均雨量基準(zhǔn)值與估算值的RMSE。

圖4 4種插值算法2006-2014年的日均降雨在不同雨量站密度下100次隨機(jī)分布的RMSE箱線圖Fig．4 Box plot of RMSE of observed and estimated average daily rainfall in 2006-2014 in the presence of different densities of rain gauges in 100 stochastic samplings by four interpolation methods

從圖4中可以看到，對于這4種插值算法，隨著雨量站數(shù)目的增加，其RMSE的值越小，其RMSE的誤差范圍也越小。

綜合點雨量和面雨量的分析得到：對于任一種空間插值算法，隨著參與計算的雨量站數(shù)目的增加，其降雨插值結(jié)果誤差越小。

4.3 不同雨量站密度下不同插值算法的徑流響應(yīng)分析

將不同雨量站密度及分布計算的面均雨量代入到新安江模型中，計算其徑流模擬的納什效率系數(shù)NSE。從圖5可以看到：①隨著雨量站數(shù)目的增加，模型的NSE范圍逐漸減小，并且NSE也會隨著參與計算的雨量站數(shù)目的增加而增大，但當(dāng)雨量站數(shù)目增加到一個閾值時，4種插值算法的NSE范圍逐漸減小并趨于一個穩(wěn)定值（即使用44個雨量站時所達(dá)到的NSE）；②泰森多邊形法在使用4個雨量站插值時，其NSE的范圍明顯大于其他3種插值方法，當(dāng)參與計算的雨量站數(shù)目增加時，這4種插值方法的NSE結(jié)果差別不大。

圖5 不同雨量站密度下100次隨機(jī)分布的插值結(jié)果代入新安江模型中的NSE箱線圖Fig．5 Box plot of Nash-Sutcliffe coefficients of Xinanjiang modelling calculated from different rain gauges’density for four different interpolation methods

由于本文是將研究區(qū)離散成空間上1 km網(wǎng)格計算日降雨，再取平均值得到面均雨量后代入新安江模型，當(dāng)雨量站數(shù)目較少時，不同插值算法的面均雨量波動均較大，其中泰森多邊形法的面均雨量波動范圍更大（見圖4（a）），新安江模型受不同插值算法的水量輸入影響較大，所以NSE波動范圍也較大；隨著雨量站數(shù)目的增加，4種插值算法計算的面均雨量差異逐漸減小，新安江模型受不同插值算法的水量輸入影響較小，因此在雨量站數(shù)目達(dá)到某個閾值后，其NSE變化范圍很小。Xu等［12］也得到了相似的規(guī)律，他指出當(dāng)雨量站數(shù)目達(dá)到一定后，新安江模型徑流模擬的NSE變化范圍較小。

4.4 相同雨量站的不同空間分布對插值算法的降雨-徑流響應(yīng)分析

本文選取了前面6種不同雨量站密度（表1），每一種密度都有100種不同的分布，根據(jù)NNI指數(shù)的最大和最小值，選出這100種隨機(jī)分布下最均勻（Best，以B表示）和最不均勻（Worst，以 W表示）的2種分布，并列出這12種不同的雨量站空間分布（圖6）。從圖6中可以看到，無論在哪種雨量站密度下，NNI指數(shù)最大時，其得到的雨量站在整個研究區(qū)空間分布均勻；而NNI指數(shù)最小時，其得到的雨量站在整個研究區(qū)空間分布聚集。所以，NNI指數(shù)大的雨量站空間分布會明顯好于NNI指數(shù)小的空間分布，因此將NNI指數(shù)作為本文定量描述雨量站空間分布的一個標(biāo)準(zhǔn)值。

圖6 6種不同雨量站密度下由NNI指數(shù)得到的雨量站空間分布Fig．6 Spatial distribution of rain gauges selected by NNI of six density levels

針對上述12種雨量站的空間分布，計算本文4種插值算法間的點雨量誤差RMSE（圖7（a））。從直方圖中的RMSE可以看到：①對于這4種插值算法，雨量站均勻的RMSE明顯要小于雨量站分布不均的RMSE，這說明，雨量站的空間分布越均勻，其降雨插值的誤差也越?。虎趯τ谌魏我环N雨量站分布情景，泰森多邊形法的RMSE明顯的要高于另外3種插值方法，而其中普通克里金的RMSE和協(xié)同克里金的RMSE結(jié)果非常相似，稍好于反距離權(quán)重法的RMSE，說明對于點雨量的估計，考慮了空間變量的克里金法（普通克里金和協(xié)同克里金）對日降雨的估計較為準(zhǔn)確。

圖7 不同雨量站分布下4種插值算法計算的RMSE和NSE的柱狀圖Fig．7 Histogram of RMSE and NSE in different distributions of rain gauges using four different interpolation methods

再將4種插值算法計算的12種雨量站分布的面均雨量代入新安江模型計算其徑流模擬系數(shù)NSE（圖7（b）），從直方圖中的NSE可以看到：①隨著雨量站數(shù)目的增加，當(dāng)雨量站數(shù)目達(dá)到某個閾值后，不論站點分布是否均勻，這4種插值方法所得到的NSE差異逐漸減小，這也從雨量站密度分布的個例分析中證明了圖5的結(jié)論；②當(dāng)雨量站點分布均勻時（圖7（b）中橫坐標(biāo)B表示的部分），這4種插值方法的徑流模擬系數(shù)NSE差別較小，說明在雨量站網(wǎng)布設(shè)均勻時，各空間插值算法插值結(jié)果差異較?。划?dāng)雨量站點分布不均勻時（圖7（b）中橫坐標(biāo)W表示的部分），在雨量站數(shù)目較少時，不同插值方法所得到的NSE差異較大，而隨著雨量站數(shù)目的增加，不同插值算法的NSE差異減小，說明在雨量站網(wǎng)布設(shè)不均勻時，站點數(shù)目越少，各空間插值算法插值結(jié)果差異越大。

綜合考慮圖7，針對這4種不同的插值算法，在計算點雨量時，考慮空間變量的克里金（包括普通克里金和協(xié)同克里金）插值算法要稍好于反距離權(quán)重法，明顯好于泰森多邊形法，所以對于估算單個雨量站點雨量或者分布式水文模型的雨量時，泰森多邊形法不能很好地估算實際降雨，克里金法相較于反距離權(quán)重法能更為準(zhǔn)確的估算日實際降雨；在計算面均雨量時，由于對流域所有網(wǎng)格取面均值后，不同插值算法間差異較小，因此可以選用計算簡便的插值算法，比如泰森多邊形、反距離權(quán)重法。

5 結(jié) 論

本文在雨量站密度及分布不同情況下，比較了4種常用的空間插值算法的日降雨結(jié)果，并將面均降雨帶入水文模型中比較其徑流的納什效率系數(shù)，得出結(jié)論如下：

（1）對于任一種空間插值算法，隨著參與計算的雨量站數(shù)目的增加，降雨插值結(jié)果誤差也減小，其徑流模擬的納什效率系數(shù)會先增大，到達(dá)一個閾值后不再明顯變化。

（2）在雨量站數(shù)目不變時，雨量站空間分布越均勻，雨量插值的結(jié)果誤差也越小，其徑流模擬的納什效率系數(shù)也越大。

（3）在雨量站網(wǎng)布設(shè)均勻時，各空間插值算法插值結(jié)果差異較?。辉谟炅空揪W(wǎng)布設(shè)不均勻時，站點數(shù)目越少各空間插值算法插值結(jié)果差異越大。

（4）在計算點雨量時，考慮空間變量的克里金（包括普通克里金和協(xié)同克里金）法要稍好于反距離權(quán)重法，明顯好于泰森多邊形法，因此克里金法能更準(zhǔn)確地估算日降雨；在計算面均雨量時，不同插值算法間差異較小，可以選用計算簡便的插值算法。