GR4J模型在錢塘江流域的適用性研究

張徐杰，張發(fā)鴻，富 強，趙建鋒

(中國電建集團(tuán)華東勘測設(shè)計研究院有限公司，浙江 杭州 311122)

0 引 言

水文模型是模擬自然界水文循環(huán)的數(shù)學(xué)物理模型，模型的輸入數(shù)據(jù)一般包括降雨等，模型輸出結(jié)果一般為徑流。水文模型根據(jù)對流域的空間離散程度可以分為集總式水文模型和分布式水文模型。由于集總式水文模型具有結(jié)構(gòu)簡單、輸入數(shù)據(jù)要求較低等優(yōu)點，其應(yīng)用程度并不亞于分布式水文模型。常見的集總式水文模型有斯坦福模型、水箱模型、新安江模型、薩克拉門托模型等[1]。

GR4J模型是來自法國的一個集總式水文模型，是由Perrin等人于2003年在GR3J模型的基礎(chǔ)上開發(fā)而來的[2]。由于其簡單實用的特性，GR4J模型在國外已被廣泛地應(yīng)用在水文模擬、洪水預(yù)報等相關(guān)研究中，并且取得了較好的效果[3-6]。而在國內(nèi)，GR4J模型的應(yīng)用還相對較少，主要有鄧鵬鑫等[7]在贛江流域日徑流模擬中的應(yīng)用，管曉祥等[8]在黃河流域徑流過程模擬的應(yīng)用，王強等[9]在遼寧干旱半干旱流域的適用性研究等。

本文主要以錢塘江流域為研究對象，采用GR4J模型進(jìn)行逐日徑流模擬，對模型的模擬效果和適用性進(jìn)行分析和評價，以進(jìn)一步拓展GR4J模型在國內(nèi)流域水文模擬中的應(yīng)用。

1 GR4J模型原理

GR4J模型通過使用兩個非線性水庫進(jìn)行產(chǎn)流匯流計算，分別稱為產(chǎn)流水庫和匯流水庫，模型結(jié)構(gòu)如圖1所示[2]。模型中共有4個參數(shù)，分別是：產(chǎn)流水庫最大蓄水量X1(mm)、地下水交換系數(shù)X2(mm)、前一天匯流水庫最大蓄水量X3(mm)、單位線UH1的匯流時間X4(d)。下面對模型的產(chǎn)流和匯流過程作簡單介紹。

圖1 GR4J模型結(jié)構(gòu)示意

(1)產(chǎn)流階段。GR4J模型中，首先根據(jù)模型的輸入即降雨和蒸發(fā)(分別記為P和E，單位mm，下同)，確定凈降雨量Pn和蒸散發(fā)能力En。若P≥E，則Pn=P-E，En=0；若P0時，凈降雨量Pn中有一部分進(jìn)入產(chǎn)流水庫，記為Ps，其大小為

(1)

式中，S為產(chǎn)流水庫的水量，mm；X1為產(chǎn)流水庫最大蓄水量，mm。剩余部分為Pn-Ps，直接進(jìn)入?yún)R流階段。當(dāng)En>0時，產(chǎn)流水庫的蒸發(fā)量Es由下式?jīng)Q定。即

(2)

然后，產(chǎn)流水庫的水量S更新為

S=S-Es+Ps

(3)

則產(chǎn)流水庫的產(chǎn)流量

Perc=s(1-[1+(4S/9X1)4]-0.25)

(4)

因此，產(chǎn)流階段的總產(chǎn)流量

Pr=Perc+Pn-Ps

(5)

(2)匯流階段。GR4J模型的匯流過程分為兩部分，總產(chǎn)流量Pr的90%采用單位線UH1進(jìn)行匯流演算，并經(jīng)過匯流水庫進(jìn)行調(diào)節(jié)，剩余的10%采用單位線UH2進(jìn)行匯流演算。兩條單位線

UH1(i)=SH1(i)-SH1(i-1)

(6)

UH2(i)=SH2(i)-SH2(i-1)

(7)

式中，i為整數(shù)，表示第i天；SH1和SH2分別為單位線UH1和UH2對應(yīng)的S累積曲線，計算式分別如下

(8)

(9)

式中，t為時間，d。模型中還考慮到這兩部分匯流水量之間存在交換過程，交換量

F=X2(R/X3)3.5

(10)

式中，R為匯流水庫的水量，mm。然后，匯流水庫的水量更新為

R=max(0,Q9+F+R)

(11)

式中，Q9為經(jīng)過單位線UH1演算后進(jìn)入?yún)R流水庫的水量，mm。匯流水庫的出流量

Qr=R{1-[1+(R/X3)4]-0.25}

(12)

經(jīng)過單位線UH2演算和地下水交換水量匯合后的出流量

Qd=max(0,Q1+F)

(13)

式中，Q1為經(jīng)過單位線UH2演算后的水量，mm。因此，流域出口總匯流水量

Q=Qr+Qd

(14)

2 模型應(yīng)用

2.1 研究區(qū)域

錢塘江流域位于中國東部，東臨東海，流域范圍介于東經(jīng)117°～122°與北緯28°～31°之間，流域面積約為5.56萬km2。整個流域大部分位于浙江省西部，有南、北兩源。其中，正源為北源新安江，發(fā)源于安徽省休寧縣；南源蘭江也發(fā)源于安徽省休寧縣，源頭海拔約為810 m。南源匯流后，向東南流入浙江省，至馬金鎮(zhèn)后稱馬金溪，再下行至常山港、衢州后稱為衢江。衢江繼續(xù)下行，在蘭溪與金華江匯流后稱為蘭江。蘭江干流長約303 km，流域面積約1.95萬km2。由于新安江上有新安江水庫影響天然徑流過程，因此本次主要選取蘭江流域作為研究區(qū)域，如圖2所示。

圖2 蘭江流域示意

2.2 基礎(chǔ)數(shù)據(jù)

本文共選取蘭江流域24個水文氣象站點1980年～1990年的逐日觀測資料進(jìn)行分析計算。其中，流量站有6個，其逐日流量資料作為GR4J模型率定和驗證的實測徑流數(shù)據(jù)。氣象代表站有2個，主要用于計算GR4J模型所需的潛在蒸發(fā)量。雨量站有23個，通過反距離權(quán)重插值后分別得到6個子流域?qū)?yīng)的逐日降雨資料，作為模型輸入數(shù)據(jù)。水文氣象站點分布情況和基本信息見圖2和表1。

表1 水文氣象站點信息

3 結(jié)果與分析

采用蘭江流域1980年～1990年的逐日水文氣象數(shù)據(jù)，對GR4J模型的適用性進(jìn)行研究和分析。其中，1980年作為模型預(yù)熱期，1981年～1986年作為模型率定期，1987年～1990年為模型驗證期。由于GR4J模型中只有4個參數(shù)，一般采用簡單的率定方法就可以達(dá)到較理想的效果。因此，本文采用拉丁超立方體抽樣法進(jìn)行10 000組參數(shù)抽樣，然后進(jìn)行模型模擬，目標(biāo)函數(shù)選用Nash-Sutcliffe系數(shù)。另外，還選用相對誤差RBIAS來一起評價模型的模擬效果。即

(15)

(16)

式中，NS為Nash-Sutcliffe系數(shù)；RBIAS為模擬徑流總量的相對誤差；Qsim為模擬徑流量，m3/s；Qobs為實測徑流量，m3/s。

圖3～圖8是6個站點率定期和驗證期的逐日模擬流量與實測流量對比圖，表2和表3分別是GR4J模型在6個站點的參數(shù)率定結(jié)果和評價指標(biāo)結(jié)果。結(jié)合以上圖表可以看出，率定期和驗證期NS系數(shù)均達(dá)到0.7以上，模擬徑流總量的相對誤差總體控制在±10%以內(nèi)(除了金華站驗證期為-11.26%)。對于流域面積較大的衢州、金華和蘭溪3個站點，其NS系數(shù)均達(dá)到0.9以上，說明模型模擬效果總體較好。但從圖3也可以看出，流域面積最小的密賽站，部分年份的峰值徑流沒有較好地模擬，說明GR4J模型在流域面積小的站點模擬效果不如流域面積大的站點。這可能是因為，在流域面積較小的山區(qū)，徑流更加陡漲陡落，隨時間的變異性更大，水文模型更加難以模擬。

表2 GR4J模型參數(shù)率定結(jié)果

表3 模型率定期和驗證期評價指標(biāo)結(jié)果

圖3 密賽站模擬和實測流量對比

圖9所示的是6個站點10 000次模擬結(jié)果的NS系數(shù)分布情況。由圖9可見，較大NS值(比如0.8以上)的頻數(shù)并不是很多，這主要是因為本文對于模型參數(shù)率定采用了拉丁超立方體抽樣法，在參數(shù)范圍內(nèi)進(jìn)行了比較均勻的抽樣，并沒有參數(shù)迭代尋優(yōu)的過程，導(dǎo)致較理想模擬結(jié)果的次數(shù)所占比例較小。對于流域面積較大的衢州、金華、蘭溪3個站點，NS系數(shù)超過0.6的模擬次數(shù)基本超過50%，說明GR4J模型在蘭江流域總體具體較好的適用性。

圖4 雙塔底站模擬和實測流量對比

圖5 義烏佛堂站模擬和實測流量對比

圖6 衢州站模擬和實測流量對比

圖7 金華站模擬和實測流量對比

圖8 蘭溪站模擬和實測流量對比

圖9 各站水文模擬NS系數(shù)分布情況

4 結(jié)論與展望

GR4J是一個結(jié)構(gòu)簡單的4參數(shù)集總式水文模型，目前在國內(nèi)的應(yīng)用不是特別普遍。本文采用GR4J模型，在錢塘江流域的蘭江子流域進(jìn)行應(yīng)用。1981年～1990年的逐日徑流模擬結(jié)果表明，6個水文站點率定期和驗證期的最大NS系數(shù)均達(dá)到0.7及以上。其中，流域面積較大的3個站點，大概50%模擬次數(shù)的NS系數(shù)超過0.6，率定期和驗證期的最大NS系數(shù)達(dá)到0.9及以上。模擬結(jié)果表明，GR4J模型在蘭江流域總體具體較好的適用性。

本次研究的不足主要包括以下兩個方面：一是參數(shù)率定的方法較簡單，沒有迭代尋優(yōu)的過程；二是目前的參數(shù)率定結(jié)果表明，模型參數(shù)X3的推薦范圍(推薦范圍出自參考文獻(xiàn)[2])可能需要進(jìn)一步擴大(密賽站X3的率定結(jié)果在推薦范圍的下限)。后續(xù)研究中，將進(jìn)一步擴大部分參數(shù)范圍，探究不同地區(qū)的參數(shù)差異性，并在參數(shù)率定方法上尋找其他更優(yōu)方法的可能性。