降雨輸入對青龍河流域HSPF模型模擬結(jié)果的影響

劉興坡，陳心能，盧木子，夏澄非，李永戰(zhàn)

(1.上海海事大學 海洋科學與工程學院，上海 201306； 2.上海海事大學 海洋環(huán)境與生態(tài)模擬研究中心，上海 201306； 3.河北省桃林口水庫管理局，河北 秦皇島 066400)

當前，BASINs/HSPF模型在國內(nèi)外流域水文水質(zhì)模擬方面得到了較為廣泛的應(yīng)用[1].評價模型模擬結(jié)果不確定性、提升模型模擬精度開始成為模型應(yīng)用的重要問題.一般而言，模型模擬結(jié)果的不確定性取決于模型結(jié)構(gòu)的不確定性、模型參數(shù)的不確定性以及模型輸入的不確定性等[2].首先，模型結(jié)構(gòu)不確定性一般需要通過模型結(jié)構(gòu)的重新研發(fā)、組合(或集合)模型[3]等途徑進行優(yōu)化，前者涉及流域水文水質(zhì)系統(tǒng)模型的系統(tǒng)開發(fā)，而后者則超出了單個流域模型研究的范疇.其次，HSPF模型參數(shù)不確定性是當前模型不確定性研究的熱點，主要包括用于不確定性分析的貝葉斯理論方法、用于靈敏度分析的方差分解方法以及用于校準優(yōu)化的啟發(fā)式智能搜索算法等[2,4-5]，典型方法如靈敏度分析[6]、PEST自動校準[7]、正交極差分析[8]等.在降雨輸入研究方面，劉潔等應(yīng)用降雨隨機模擬模型分析了3種降雨情景對東江流域徑流的影響[9].該文從降雨情景分析的角度考慮了降雨隨機性變化，但并未就其對HSPF模型精度影響作出評估.ZHOU Zhang等則針對氣候變化導致的降雨模式變化，應(yīng)用BASINS-HSPF-CAT建模系統(tǒng)評估其對中國海南尖峰嶺熱帶雨林流域水文過程的影響[10].綜上發(fā)現(xiàn)，對于HSPF模型降雨輸入不確定性對其模擬結(jié)果影響的研究還未見相關(guān)報道.鑒于降雨是HSPF模型最為重要的輸入條件，具有非常大的波動性和隨機性.本文針對降雨輸入較大的波動性和隨機性，應(yīng)用降雨隨機模擬方法來考查其對HSPF模型模擬結(jié)果的影響，為HSPF模型的深入應(yīng)用提供參考和借鑒.

1 研究數(shù)據(jù)與方法

1.1 研究流域及數(shù)據(jù)

以青龍河流域HSPF模型為研究對象考查降雨輸入的影響.青龍河干流總長246 km，流域面積6 340 km2，流域范圍大部分在河北省青龍滿族自治縣、平泉縣和寬城縣境內(nèi)，部分隸屬遼寧省凌源市和建昌縣，其中河北境內(nèi)3 776 km2，遼寧省境內(nèi)1 284 km2.受太平洋副熱帶高壓的影響，該流域暴雨多發(fā)生在7月下旬—8月上旬，兩月雨量約占全年降雨量的60%～70%，最大可達80%.

本文研究數(shù)據(jù)主要包括：1)精度為90 m的SRTM3文件，經(jīng)ArcGIS10.2處理后可獲得流域DEM圖；2)應(yīng)用ArcGIS10.2軟件對2009年0.5 m空間分辨率的航拍影像進行目視解譯，并按照GlobCover2009數(shù)據(jù)的分類體系，可獲得流域12種土地利用類型數(shù)據(jù)；3)基于流域DEM數(shù)據(jù)，應(yīng)用ArcGIS10.2軟件進行流域水文系統(tǒng)分析，包括水流方向提取、水流長度計算、河網(wǎng)分級、河流網(wǎng)絡(luò)獲取、匯流累積量計算以及流域分割(子流域獲取)等步驟，從而可獲得HSPF建模所需要的河網(wǎng)和子流域.本研究中，將整個流域劃分為35個子流域，分別對應(yīng)35個河段，如圖1所示.4)氣象數(shù)據(jù)來自中國氣象科學數(shù)據(jù)共享中心網(wǎng)站河北省青龍站(東經(jīng)118.950°，北緯40.400°)，該站點數(shù)據(jù)包括1957年—2013年逐日降雨、相對濕度、氣溫、風速、光照、蒸發(fā)皿蒸發(fā)和氣壓等.應(yīng)用WDMUtil程序可將氣象數(shù)據(jù)計算轉(zhuǎn)換為小時數(shù)據(jù)，生成wdm格式文件可為HSPF模擬使用.本文研究所采用的即為1957年1月1日—2013年9月30日的日降雨數(shù)據(jù)，數(shù)據(jù)來源為中國氣象科學數(shù)據(jù)共享中心網(wǎng)站.針對從該網(wǎng)站下載的數(shù)據(jù)文件，首先利用Surfer軟件和Excel對降雨數(shù)據(jù)進行格式編排以及異常值的修正，然后應(yīng)用Surfer軟件進一步將降雨數(shù)據(jù)轉(zhuǎn)化為WDM時間序列格式，作為HSPF模型的輸入數(shù)據(jù).

圖1 青龍河流域HSPF模型界面

1.2 建模方法

根據(jù)隨機水文學原理，降雨時間序列包括趨勢成分、周期成分、相依隨機成分和獨立隨機成分(白噪聲)等.其中，趨勢成分和周期成分為降雨時間序列的確定性成分.當降雨時間序列扣除確定性成分后便成為平穩(wěn)隨機序列.該平穩(wěn)隨機序列包括相依隨機成分和獨立隨機成分兩個部分.為此，分別對上述各個時間序列成分進行建模，建立降雨時間序列的隨機模擬模型，采用蒙特卡洛方法，獲得多組降雨隨機模擬序列，逐個作為HSPF模型的降雨輸入進行模擬，獲得每次模擬的Nash-Sutcliffe效率系數(shù)(ENS，下同)，據(jù)此來評價降雨輸入對HSPF模型模擬效果的影響.研究技術(shù)路線參見圖2.

1.2.1 趨勢成分檢驗與建模方法

首先分別采用非參數(shù)Mann-Kendall趨勢檢驗法、Mann-Kendall秩次相關(guān)檢驗法、Spearman秩次相關(guān)檢驗法、線性回歸法以及滑動平均法評估實際降雨時間序列(Z(t))的趨勢顯著性.然后采用如下方法對降雨時間序列的趨勢成分進行建模.降雨時間序列趨勢成分可選取為如下形式[11]：

c8t-1/2+c9e-t+c10lnt.

(1)

式中：t為時間變量，ci(i=0,1,2,…,10)為待定系數(shù).

然后采用逐步回歸分析方法確定函數(shù)形式.當樣本數(shù)目很大時，可選用相對的時間單位，即有

(2)

式中：N為樣本數(shù)目，τ為采樣間隔.

圖2 研究路線

1.2.2 周期成分建模方法

降雨時間序列的周期項p(t)可表達為d個諧波函數(shù)的組合，即

(3)

(4)

(5)

式中，Y(t)=Z(t)-f(t)，為實際序列去除趨勢項之后的殘差序列.

確定諧波個數(shù)d可采用累積周期圖的方法[12].具體方法如下：

1)計算序列Y(t),t=1,2,…,N的方差：

(6)

2)計算第j個諧波對序列方差S2的貢獻：

(7)

(8)

(9)

1.2.3 平穩(wěn)時間序列建模方法

1)相依隨機成分建模方法.實際降雨時間序列去除趨勢成分和周期成分后，就可按照平穩(wěn)時間序列模型進行建模.其中的相依隨機成分建模一般采用自回歸滑動平均模型(ARMA)，包括AR(p)和MA(q)兩部分.本文采用最小信息準則(AIC準則和BIC準則)來確定ARMA(p,q)模型階數(shù)p和q，并采用矩法來估計模型參數(shù).相依隨機成分的建模表達式一般可以表示為

s(t)=β0+β1s(t-1)+β2s(t-2)+…+βns(t-n)+
e(t)+α1e(t-1)+α2e(t-2)+…+αde(t-d).

(10)

式中：s(t)和e(t)分別代表自回歸模型和滑動平均模型，α和β分別代表相應(yīng)系數(shù).

2)獨立隨機成分(白噪聲)建模方法.根據(jù)建立的平穩(wěn)時間序列模型，可以反求殘差序列.殘差成分的均值應(yīng)為0，方差可經(jīng)過統(tǒng)計計算得出.由此就可以確定獨立隨機成分(白噪聲成分)的數(shù)學表達形式.殘差成分的正態(tài)獨立性需要在建模后驗證.若已知獨立隨機成分的均值和方差，可按照正態(tài)分布對其進行隨機模擬.

3)殘差成分的正態(tài)獨立性檢驗.進行上述建模后，需要進行殘差成分的正態(tài)獨立性檢驗，驗證所建立的模型是否滿足要求，能否反映實際降雨時間序列的時序特性和隨機變化規(guī)律.

由于平穩(wěn)序列模型ARMA(p,q)是建立在具有正態(tài)概率分布特性的獨立隨機成分(白噪聲成分)基礎(chǔ)之上.在模型建立之后，必須對殘差成分的獨立正態(tài)性進行檢驗.一般采用自相關(guān)函數(shù)檢驗法，首先構(gòu)造統(tǒng)計量：

(11)

若殘差序列為獨立序列，則統(tǒng)計量R服從自由度為(m-p-q)的χ2分布.公式中m為最大滯時(一般可選取為N/4左右)，N為序列長度，rk為殘差序列的k階自相關(guān)系數(shù).

1.2.4 降雨時間序列的蒙特卡洛隨機模擬方法

2 降雨隨機序列模擬結(jié)果

2.1 趨勢成分檢驗與建模結(jié)果

2.1.1 趨勢成分檢驗結(jié)果

根據(jù)常用的5種時間序列趨勢性檢驗方法，對青龍河流域的實際降雨時間序列進行分析(見表1).根據(jù)檢驗結(jié)果，Mann-Kendall秩次相關(guān)檢驗法的檢驗結(jié)果為顯著，其余4種方法均為不顯著.為了保證趨勢分析的可靠度，繼續(xù)通過定量分析的方法識別青龍河流域降雨時間序列的趨勢成分.

表1 降雨時間序列趨勢檢驗結(jié)果比較

2.1.2 趨勢成分確定

為了定量確定青龍河流域?qū)嶋H降雨時間序列的趨勢成分，基于Matlab平臺，采用逐步回歸分析方法確定其趨勢成分，如圖3所示.根據(jù)式(1)，設(shè)定95%置信限，獲得趨勢成分的擬合公式(見式(11)).根據(jù)式(12)，獲得降雨時間序列的趨勢成分，如圖4所示.

f(t)=0.026 7t2+0.053 4t-1+0.133lnt.

(12)

圖3 Matlab平臺上應(yīng)用逐步回歸分析法優(yōu)化趨勢成分

Fig.3 Stepwise regression analysis for optimizing trend components on Matlab

圖4 趨勢成分建模結(jié)果

2.2 周期成分建模結(jié)果

根據(jù)累積周期圖建模方法，可以發(fā)現(xiàn)實際降雨序列存在約13年的變化周期，與小波分析、功率譜方法獲得的周期結(jié)果相同.式(13)為實際降雨序列周期成分表達式.周期成分的建模結(jié)果如圖5所示.

(13)

圖5 周期成分結(jié)果

2.3 相依隨機成分建模結(jié)果

首先根據(jù)AR模型以及MA模型的偏相關(guān)函數(shù)判斷ARMA模型的階次，如圖6所示.然后應(yīng)用矩法進行參數(shù)估計，從而獲得3個ARMA模型，如表2所示.分別計算上述3個模型的AIC值和BIC值，并進行比較.由表2可知，ARMA(1,1)模型的AIC值、BIC值最小，故根據(jù)最小準則選擇其作為相依隨機成分的建模形式.

2.4 獨立隨機成分(白噪聲)建模結(jié)果

圖6 ARMA模型偏相關(guān)函數(shù)

表2 ARMA模型優(yōu)選

2.5 隨機模擬模型實用性檢驗

為了檢驗隨機模擬序列是否能夠反映實際降雨時間序列的大部分統(tǒng)計特性.此處采用長序法對該模型的實用性進行初步分析.檢驗內(nèi)容包括均值、離差系數(shù)Cv、偏差系數(shù)Cs和一階自相關(guān)系數(shù).按相對誤差≤10%統(tǒng)計參數(shù)通過率，結(jié)果見表3.

表3 統(tǒng)計參數(shù)通過率

根據(jù)表3中各項統(tǒng)計數(shù)據(jù)，降雨隨機模擬序列通過了長序檢驗，驗證了降雨隨機模擬方法及其結(jié)果的可信度.這表明，采用趨勢成分、周期成分與ARMA建模以及正態(tài)隨機模擬獲得降雨隨機模擬序列是量化降雨輸入隨機性的可行方法.

2.6 降雨隨機模擬序列生成

獲得降雨時間序列的趨勢成分、周期成分、相依隨機成分以及獨立隨機成分之后，就可以應(yīng)用蒙特卡羅方法，在計算機上可實現(xiàn)降雨時間序列的隨機模擬，獲得與實際降雨時間序列統(tǒng)計特性相似的隨機模擬序列.圖7所示為200個降雨隨機模擬序列與實際降雨時間序列的比較.

3 HSPF模型模擬結(jié)果與分析

為了分析降雨時間序列隨機性對于青龍河流域HSPF模型模擬效果的影響，將200組降雨隨機模擬序列分別輸入到HSPF模型中，計算HSPF模型驗證期(2011年)實測徑流與模擬徑流的ENS值，評價模擬效果.

圖7 降雨隨機模擬序列與實際降雨序列的比較

Fig.7 Comparison of rainfall random simulation series and actual rainfall series

3.1 青龍河流域HSPF模型靈敏參數(shù)優(yōu)化

為了合理評價降雨時間序列對HSPF模型模擬結(jié)果的影響，首先需要優(yōu)化HSPF模型參數(shù)，保證HSPF模型精度.為此，應(yīng)用PEST自動校準程序識別出9個最為靈敏的模型參數(shù)，然后采用正交極差分析方法優(yōu)化獲得9個參數(shù)的校準結(jié)果[8]，即LZSN為2，INFILT為0.167，AGWRC為0.95，DEEPFR為0.209，BASETP為0.133，AGWETP為0.1，CEPSC為0.27，UZSN為1.35，IRC為0.483.其余不靈敏參數(shù)則采用HSPF模型設(shè)定默認值，并進行相應(yīng)參數(shù)不確定性的影響分析.

3.2 降雨隨機模擬序列對HSPF模型模擬結(jié)果的影響

將上述獲得的200組降雨隨機模擬序列分別輸入到HSPF模型中運行模擬，可以獲得200組模擬徑流序列(如圖8所示).并根據(jù)2011年實測徑流數(shù)據(jù)與HSPF模型的模擬徑流數(shù)據(jù)，計算獲得200組ENS值，如表4所示.根據(jù)表4，200次HSPF模擬獲得的ENS值的變化范圍為[71.09%, 74.96%]，ENS平均值為73.03%，ENS波動幅度為3.87%.由此可以發(fā)現(xiàn)，在模型參數(shù)得到較好優(yōu)化的前提下，降雨時間序列隨機性仍會導致HSPF模擬ENS值3.87%的起伏.這表明降雨輸入的隨機性是影響HSPF模型模擬效果的顯著因素，是改善HSPF模型模擬效果的重要考量.

圖8 模擬徑流序列與實測徑流序列比較

Fig.8 Comparison of the simulated runoff series and the actual runoff series

表4 HSPF隨機模擬的ENS值

Tab.4 Nash-Sutcliffe efficiency (ENS) values by HSPF stochastic simulations %

3.3 降雨時間序列中極值點對HSPF模擬結(jié)果的影響

由于一般的隨機建模方法無法充分還原實際降雨時間序列極值點或?qū)O值點作為奇異點濾除，從而導致隨機模擬序列較實際時間序列更為平滑.因此，為了評估降雨時間序列中極值點對HSPF模擬結(jié)果的影響，此處在隨機模擬序列中考慮將實際降雨時間序列的若干極值點還原，從而考查極值點對HSPF模擬結(jié)果的影響.選用3.2節(jié)中ENS最高值(74.96%)所對應(yīng)的降雨隨機模擬序列，按照每年中日降雨量由大到小順序，依次還原1個、2個、3個、4個和5個極值點，分別獲得5個對應(yīng)的降雨隨機模擬序列，然后依次輸入到HSPF模型中，評估其對模擬結(jié)果的影響,結(jié)果如表5所示.由表5可知，經(jīng)過5次還原極大值后，HSPF模擬ENS的最小值為75.35%，均大于沒有還原極值的降雨隨機模擬序列(其ENS最大值為74.96%)，這表明降雨時間序列極大值對于HSPF模擬效果具有顯著影響.因此，提高HSPF模擬精度需要關(guān)注降雨時間序列極大值點的影響.此外，隨著每年還原的日降雨極大值數(shù)量從1增長到3時，ENS值逐漸增大，每年還復極大值數(shù)量為3時，ENS值達到最大，而后逐漸下降.這說明降雨時間序列還原3個最大極值點時模擬結(jié)果最佳.這也啟示，在HSPF模擬的降雨情景優(yōu)選時，還原的極大值點并非越多越好，而是要重點考慮若干最大極值點的影響.

表5 還原的日降雨極值個數(shù)與相應(yīng)的ENS值

Tab.5 Number of modified daily rainfall extremes and the corresponding ENS values %

同理，為了考查降雨時間序列極小值點對模擬結(jié)果的影響，經(jīng)過5次還原極小值后，HSPF模擬ENS的最小值為75.03%，也均大于沒有還原極值的降雨隨機模擬序列(其ENS最大值為74.96%)，這表明降雨時間序列極小值點對于模擬結(jié)果具有優(yōu)化作用.然而，還原極小值之后，HSPF模擬的ENS值的波動幅度僅為0.16%，這表明將降雨模擬序列中每年日降雨量的極小值還原后，對于HSPF模型的影響并不顯著.此外，隨著每年還原的日降雨極小值數(shù)量從1增長到4時，ENS值逐漸增大，每年還復極小值數(shù)量為4時，ENS值達到最大，而后下降.這說明降雨時間序列還原4個最大極值點時模擬結(jié)果最佳.這也啟示，在HSPF模擬的降雨情景優(yōu)選時，還原的極小值點并非越多越好，而是要重點考慮若干最小極值點的影響.

通過降雨輸入隨機性和極值點對HSPF模型模擬結(jié)果的影響分析，發(fā)現(xiàn)在模型參數(shù)相對優(yōu)化的條件下，兩個因素對于HSPF模型模擬結(jié)果影響具有7.33%的貢獻率，這表明降雨輸入的不確定性是HSPF模型模擬結(jié)果不確定性的重要來源，也是提高HSPF模型精度的重要方面.

4 結(jié) 論

1)采用趨勢成分、周期成分與ARMA建模以及正態(tài)隨機模擬獲得降雨隨機模擬序列是量化降雨輸入隨機性的可行方法.

2)在模型參數(shù)優(yōu)化的條件下，降雨隨機模擬序列HSPF模擬ENS值的變化區(qū)間為[71.09%,74.96%]，波動幅度達3.87%，表明降雨輸入隨機性對于HSPF模擬結(jié)果具有顯著影響.

3)當考慮每年的日降雨量極大值時，ENS值變化區(qū)間為[75.35%,78.81%]，波動幅度達3.46%，且結(jié)果均優(yōu)于沒有考慮降雨極值點的降雨隨機模擬序列，表明降雨時間序列的極大值對于HSPF模擬效果具有顯著影響.

4)隨著降雨時間序列中所考慮極大值點數(shù)量的逐漸增多，HSPF模擬效果出現(xiàn)下降趨勢，表明HSPF模擬應(yīng)特別關(guān)注若干最大極值點的影響.

5)當考慮每年的日降雨量極小值時，ENS值變化區(qū)間為[75.03%,75.19%]，波動幅度僅為0.16%，表明降雨時間序列的極小值對于HSPF模擬效果的影響并不顯著.

6)降雨輸入的不確定性是HSPF模型模擬結(jié)果不確定性的重要來源，改善HSPF模擬效果需要考慮降雨時間序列隨機性和極值點因素的影響.

本研究可為量化降雨輸入對HSPF模型模擬的影響以及HSPF模擬降雨情景的選擇提供借鑒.