趙思遠(yuǎn)，郝振純，劉文斌

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098；2.黑龍江省水文局，黑龍江 哈爾濱 150001)

黑龍江干流水位預(yù)報方法研究*

趙思遠(yuǎn)1，郝振純1，劉文斌2

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098；2.黑龍江省水文局，黑龍江 哈爾濱 150001)

選用基于混合線性回歸模型的黑箱模型方法對黑龍江干流中上游段包括漠河、鷗浦、三道卡、長發(fā)屯、烏云、勤得利6個站點進(jìn)行水位過程模擬，探討我國江河河道的水位預(yù)報問題。模型結(jié)構(gòu)采用AIC準(zhǔn)則和DW檢驗綜合確定，參數(shù)通過最小二乘估計，模型以2006-2010年汛期(5-9月)日時間尺度的水文觀測資料進(jìn)行參數(shù)率定，以2011年汛期日時間尺度的水文觀測資料進(jìn)行模擬驗證；并對模型輸出不確定性，參數(shù)敏感性進(jìn)行分析研究。結(jié)果表明：① 所提出的模型方法在黑龍江干流地區(qū)應(yīng)用較好，6個站點的預(yù)報合格率均達(dá)到85%以上，確定性系數(shù)也都為0.97以上，預(yù)報方案精度均能達(dá)到甲等，可以用于今后干流作業(yè)預(yù)報。② 通過統(tǒng)計Monte Carlo模擬的輸出變量標(biāo)準(zhǔn)差量化了模型輸出不確定性，標(biāo)準(zhǔn)差越大，說明預(yù)報的洪水過程線越不確定，從而模型不確定性也就越大；運用標(biāo)準(zhǔn)回歸系數(shù)敏感性分析得出，模型中干流入流不同時間階數(shù)中最小的參數(shù)最為敏感，支流入流不同時間階數(shù)中最小的參數(shù)次之，時間階數(shù)大的參數(shù)一般不太敏感。

水位預(yù)報；混合線性回歸；黑龍江干流；不確定性分析；蒙特卡羅方法

近幾年我國江河河道水文水位監(jiān)測站點大幅增加，由于觀測條件的限制，許多站點為無人值守的自動水位觀測站，無水位流量關(guān)系，但在河道防汛中，水位是最主要的防洪依據(jù)，防洪決策也需要通過獲知實時水位信息以進(jìn)行有效的防洪調(diào)度和提前安排緊急防洪措施，來減少洪水帶來的損失。因此，研究建立黑龍江各站點水位預(yù)報工作對黑龍江干流及沿岸的防洪有著重要的意義。

國內(nèi)包括黑龍江地區(qū)在河道水位預(yù)報方面，多年來一直主要沿用古老的相應(yīng)水位法和流量演算-水位流量關(guān)系轉(zhuǎn)換法[1-2]，又稱為經(jīng)驗分析法。采用該法建立的預(yù)報模型，用于有回水頂托、沖淤變化和人工控制河道，預(yù)報精度不高，一般僅用于預(yù)報洪峰水位，不宜用于預(yù)報水位過程。自20世紀(jì)70年代計算水力學(xué)問世以來，國內(nèi)外提出了一系列基于水力學(xué)理論的水位計算模型[3]，但大多只適用于穩(wěn)定河床。實際應(yīng)用中需較詳細(xì)的河道斷面和糙率資料，且河床特征常受沖淤影響隨時間變化，故預(yù)報精度不高。系統(tǒng)理論和方法向洪水預(yù)報技術(shù)的滲透是從50年代開始的。它經(jīng)歷了一個從個別方法的借用逐步到系統(tǒng)理論引入的過程[4-6]。水文統(tǒng)計方法就是基于系統(tǒng)論的觀點和思維方式，針對水位、流量等預(yù)報對象，根據(jù)流域?qū)嶋H情況和物理成因，尋找主要的相關(guān)因子，建立預(yù)報對象和主要相關(guān)因子之間的相關(guān)關(guān)系，利用相關(guān)因子的實際取值來預(yù)測預(yù)報對象的取值。這種方法簡單實用，但模型參數(shù)沒有明確的物理意義，屬于黑箱模型。采用不同系統(tǒng)計算方法結(jié)果差異較大，所以預(yù)報精度很大程度上依賴于水文資料及模型的求解。

黑龍江是以雨水補(bǔ)給為主、積雪融水補(bǔ)給為輔的河流，全部徑流中雨水補(bǔ)給約占75%～80%，融雪水補(bǔ)給約占15%～20%，干流汛期5-9月降水量占年降水量的70%～85%，且全年徑流量的約70%也在這一期間[7]。干流周圍缺乏雨量站和溫度觀測點，且缺少詳細(xì)的河道工情信息，針對黑龍江干流的水文特性和資料情況，本文采用基于混合線性回歸模型的黑箱模型方法來建立包括漠河、鷗浦、三道卡、長發(fā)屯、烏云、勤得利6個站點的水位預(yù)報模型，并應(yīng)用蒙特卡羅方法來量化各個站點混合線性回歸模型輸出不確定性，對抽樣的輸入?yún)?shù)進(jìn)行敏感性分析，目的是將各參數(shù)對模型輸出的不確定性的影響進(jìn)行排序，找出顯著的影響參數(shù)。

1 混合線性回歸模型

設(shè)有mv個變量X1,X2,…,Xmv，其中X1是因變量。X1,X2，…，Xmv的觀察值序列為Xij，i=1,…,mv，j=1,…,nt。定義混合線性回歸模型

X1t=b0+b11X1,t-p1 0 -1+ … +b1p1X1,t-p1 0 -p1+…

+bmv,1Xmv,t - pmv 0 -1+ … +bmv,pmvXmv,t-pmv 0 -pmv

(1)

式中，pi0為各變量的滯后時段；pi為各變量的階數(shù)，pi0+pi，i=1,…,mv為各變量的時間滯后。

因變量X1t不僅可以依賴于其自身歷史值X1j，j=t-p10-1,…,t-p10-p1，而且還依賴于自變元的某些歷史值Xij，j=t-pi0-1,…,t-pi0-pi，i=2,…,mv。模型的預(yù)見期

可以用最小二乘估計法進(jìn)行模型的參數(shù)估計。經(jīng)典的線性回歸模型沒有體現(xiàn)出因變量自身前后及因變量與各自變量“過去”的依賴關(guān)系；混合線性回歸模型又與自回歸模型不同，自回歸模型沒有考慮其他變量對因變量的影響。

用近似最小信息準(zhǔn)則挑選回歸因子的思想是：從高階模型到低階模型逐步剔除不重要變量，使準(zhǔn)則函數(shù)為最小的模型便是近似最佳模型。

此外，本文還引進(jìn)DW統(tǒng)計量來檢驗?zāi)Ｐ褪欠翊嬖谧韵嚓P(guān)性。DW(Durbin-Waston)檢驗自相關(guān)性的一種經(jīng)典方法。DW統(tǒng)計量為：

2 模型不確定性分析方法

模型不確定性分析是度量不確定性和描述不確定性的過程。在模型的構(gòu)建與應(yīng)用過程中都存在大量的不確定性，其中模型參數(shù)的不確定性反映了模型結(jié)構(gòu)不完善、信息不完備以及參數(shù)估計方法不足[11]，故本文著重模型參數(shù)不確定性的分析定量。對于模型參數(shù)的不確定性分析，目前最常用的是Monte Carlo分析方法[12-13]。

用于不確定性分析的蒙特卡羅方法一般包括以下幾個步驟[14]：① 獲取輸入?yún)?shù)概率分布的隨機(jī)取樣；② 將所有參數(shù)取樣值的組合輸入模型，執(zhí)行模型模擬；③ 對模型輸出結(jié)果進(jìn)行統(tǒng)計分析。本文中，先將蒙特卡羅抽樣取得的輸入樣本點通過混合線性回歸模型得到輸出樣本點，進(jìn)行統(tǒng)計分析量化模型輸出不確定性，再利用標(biāo)準(zhǔn)回歸系數(shù)(SRC)來描述各輸入?yún)?shù)敏感性[15]，利用該系數(shù)的相對大小來確定各輸入?yún)?shù)對模型輸出的不確定性的影響。

其中S(Y)和S(Xi)分別代表取樣的輸出值Y和輸入?yún)?shù)Xi的標(biāo)準(zhǔn)差，bi為混合線性回歸模型中確定的輸入?yún)?shù)的系數(shù)。每個參數(shù)的敏感性就可以通過SRC的絕對值來評價。

3 預(yù)報模型建立

黑龍江中上游的水位觀測站共有15個，即漠河、開庫康、鷗浦、呼瑪、三道卡、黑河、長發(fā)屯、勝利屯、奇克、烏云、嘉蔭、蘿北、中興鎮(zhèn)、同江、勤得利，圖1為黑龍江中上游水系概況。在黑龍江干流范圍內(nèi)，根據(jù)《SL250-2000水文情報預(yù)報規(guī)范》[16]對資料的要求，按照易獲取、資料連續(xù)、代表性強(qiáng)等原則，選取漠河、鷗浦、三道卡、長發(fā)屯、烏云、勤得利6個代表站點建立水位預(yù)報模型。為了盡可能準(zhǔn)確地預(yù)報各站點水位，消除水位的非平穩(wěn)性影響，因此對水位時間序列作一階差分變換，同時考慮了上下邊界雨量因素對水位漲差的影響。整個預(yù)報系統(tǒng)流程如圖2所示。

圖1 黑龍江中上游水系概況Fig.1 Water generalization for the upper and middle reaches of the Heilongjiang River

圖2 黑龍江干流站點水位預(yù)報方案組成概化Fig.2 Water-level forecast scheme generalization for the stations of the Heilongjiang mainstream

設(shè)漠河、開庫康、鷗浦、呼瑪、三道卡、黑河、長發(fā)屯、勝利屯、奇克、烏云、嘉蔭、蘿北、中興鎮(zhèn)、同江、勤得利站點水位時間序列分別為Z1t,…,Z15t，表示第t日8時水位(單位m)；對各站點水位作相應(yīng)的一階差分記為X1t=Z1t-Z1t-1，X2t,…,X15t同前，稱為各站點第t日8時水位漲差(單位m)；P1t，P2t,...,P15t代表各站分區(qū)第t-1日8時至第t日8時雨量(單位mm)。

建立各站點混合線性回歸模型，計算步驟為：將該站點上邊界入流t日8時以前水位漲差及匯流分區(qū)t-1日8時以前雨量作為式(1)的輸入，計算得該站點t日8時水位漲差預(yù)報值，再加上該站點t-1日8時水位，可得該站點t日8時水位預(yù)報值。

漠河站水位過程預(yù)報模型為：

(2)

式中，X0t為洛古河站第t日8時水位漲差。

鷗浦站水位過程預(yù)報模型為：

(3)

三道卡站水位預(yù)報模型為：

(4)

式中，Xct為呼瑪橋站第t日8時水位漲差。

長發(fā)屯站水位過程預(yù)報模型為：

(5)

式中，Xdt為小沙贊卡站第t日8時水位漲差。

烏云站水位過程預(yù)報模型為：

(6)

式中，Xet為雙河屯站第t日8時水位漲差；Xft為卡緬卡站第t日8時水位漲差。

勤得利站水位過程預(yù)報模型為：

(7)

式中，Xgt為富錦站第t日8時水位漲差。

4 預(yù)報精度評定

采用2006-2010年黑龍江干流各水位站點汛期(5-9月)水文觀測資料，分別率定相應(yīng)的水位過程預(yù)報模型，采用2011年汛期的水位觀測資料驗證各站的水位過程預(yù)報模型，每個站點的有效預(yù)見期為24 h。根據(jù)《SL250-2000水文情報預(yù)報規(guī)范》[16]，可以采用許可誤差和確定性系數(shù)作為方案精度評定標(biāo)準(zhǔn)，變幅均方差σΔ和確定性系數(shù)DC按下式計算：

(8)

(9)

按上述評定準(zhǔn)則，所有6個站點的預(yù)報合格率均達(dá)到85%以上，確定性系數(shù)也都為0.97以上，因此各站點預(yù)報方案均評定為甲等預(yù)報方案。各站的驗證結(jié)果見表1，預(yù)報模型的擬合效果見圖3-4。其中鷗浦站擬合曲線在預(yù)報時段中期有突變現(xiàn)象，原因分析為鷗浦站上游入流站點開庫康2011年7月16號數(shù)據(jù)觀測或摘錄可能有誤。

表1 各站水位過程預(yù)報精度統(tǒng)計

圖3 2011年汛期(6-9月)水位過程預(yù)報檢驗Fig.3 Verifying the validity of water-level forecasting for flood season during June to September,2011

5 模型不確定性分析

敏感性分析也是一種動態(tài)不確定性分析，它用于分析輸入?yún)?shù)對模型輸出的不確定性的影響排序，找出顯著的影響參數(shù)。要對6個站點的混合線性回歸模型進(jìn)行不確定性分析，首先要獲取輸入?yún)?shù)先驗概率分布，由于本模型輸入?yún)?shù)可能包含多階歷史值，在對其統(tǒng)計分布特征進(jìn)行設(shè)定時，同一種輸入?yún)?shù)的分布不會因選取不同時間階數(shù)改變，如代表漠河站的逐日水位漲差時間序列X1t與X1t-1的分布相同。不同種輸入?yún)?shù)的統(tǒng)計分布特征由自然規(guī)律并結(jié)合各站點實際情況確定。本文中輸入?yún)?shù)總體分為水位漲差和雨量兩種，而根據(jù)以往經(jīng)驗和自然規(guī)律，本文假定日水位漲差X服從正態(tài)分布，日降雨量P服從指數(shù)分布，各站點輸入?yún)?shù)的先驗分布具體情況見表2。

根據(jù)上表中假定的先驗分布，對輸入?yún)?shù)進(jìn)行蒙特卡羅抽樣并進(jìn)行模型模擬。該步驟重復(fù)1000次以獲得模型輸出值的統(tǒng)計分布信息，見表3。各站點模型輸出經(jīng)過驗證近似服從正態(tài)分布，且輸出水位漲差序列近似無偏。從各站蒙特卡羅模擬得到的輸出樣本點統(tǒng)計可看，各站樣本點標(biāo)準(zhǔn)差同先驗分布十分接近，證明所用的蒙特卡羅抽樣方法可性，且統(tǒng)計輸出的標(biāo)準(zhǔn)差量化了模型的不確定性，標(biāo)準(zhǔn)差越大，說明預(yù)報的洪水過程線越不確定，從而模型不確定性也就大一些；相反，標(biāo)準(zhǔn)差越小就說明不確定性越小。

表2 各站點混合線性回歸模型參數(shù)的統(tǒng)計分布特征

表3 輸出值統(tǒng)計分布概況

圖4 2011年汛期(5-9月)水位過程預(yù)報檢驗Fig.4 Verifying the validity of water-level forecasting for flood season of 2011(May to September)

由于本文是以線性回歸模型來逼近具有非線性特性的水位漲差，所以，采用標(biāo)準(zhǔn)回歸系數(shù)方法(SRCs)來表示每個輸入?yún)?shù)對輸出結(jié)果的貢獻(xiàn)率，其絕對值越大則表明該輸入?yún)?shù)越敏感，該模型結(jié)果不確定性的貢獻(xiàn)率也越大，各站點輸入?yún)?shù)敏感性統(tǒng)計指標(biāo)SRC列于表4。

表4 輸入?yún)?shù)的標(biāo)準(zhǔn)回歸系數(shù)

從表4可看出，對于各個預(yù)報站點的混合線性回歸模型，其上邊界干流入流站點不同時間階數(shù)中最小的水位漲差(如X4t-1)始終是對模型模型輸出不確定性貢獻(xiàn)率最大的；若上邊界有支流入流，則支流站點不同時間階數(shù)中最小的水位漲差(如Xct-1)參數(shù)敏感性次之；時間階數(shù)大的參數(shù)一般不太敏感；這些也都符合上下游相應(yīng)水位關(guān)系及匯流時間規(guī)律。此外，有兩點需要注意的是：各站自身水位漲差歷史值并不一定是顯著的影響參數(shù)(如X1t-1,X10t-1)，這是因為當(dāng)通過上游入流站點的水位漲差已經(jīng)可以滿足模型模擬輸出時，那么該站自身水位漲差歷史值對模型輸出不確定性就基本沒有貢獻(xiàn)；雨量因素在混合線性回歸模型中相關(guān)系數(shù)一直是最小的，但理論上該變量是水文循環(huán)中最活躍的，也是造成水位漲落的直接因素，這里通過標(biāo)準(zhǔn)回歸系數(shù)(SRC)指標(biāo)就量化出了它對模型輸出不確定性的貢獻(xiàn)率大小，在預(yù)報因子數(shù)目較少的情況下，它(如P1t-1)甚至對模型輸出不確定性有關(guān)鍵性影響。

6 結(jié) 論

針對黑龍江干流汛期洪水的特點和防洪減災(zāi)決策的需求，本文采用基于混合線性回歸模型的統(tǒng)計相關(guān)方法，建立了黑龍江干流各站點的水位預(yù)報模型，然后對模型輸出不確定性，參數(shù)敏感性進(jìn)行分析研究。通過預(yù)報模型應(yīng)用與不確定性分析結(jié)果能夠更好的為實際預(yù)報工作決策提供依據(jù)，從本研究中主要得出以下結(jié)論：

1)從模型模擬的精度評定指標(biāo)與模擬過程線來看，考慮了多階自變元歷史值的混合線性回歸模型模擬結(jié)果較好，6個站點的預(yù)報合格率均達(dá)到85%以上，確定性系數(shù)也都為0.97以上，預(yù)報方案的預(yù)報精度均能達(dá)到甲等。因此，所提出的模型方法具有一定的精度和較好的適用性，可以用于作業(yè)預(yù)報。

2)將MonteCarlo抽樣樣本代入模型計算，通過統(tǒng)計輸出的標(biāo)準(zhǔn)差量化了模型的不確定性，標(biāo)準(zhǔn)差越大，說明預(yù)報的洪水過程線越不確定，從而模型不確定性也就大一些；相反，標(biāo)準(zhǔn)差越小就說明不確定性越小。再通過標(biāo)準(zhǔn)回歸系數(shù)法對各站混合線性回歸模型中輸入?yún)?shù)作敏感性分析，結(jié)果表明模型上邊界干流入流不同時間階數(shù)中最小的參數(shù)變化一般對模型輸出影響最大，最為敏感，如有支流，則支流入流不同時間階數(shù)中最小的參數(shù)敏感性次之；時間階數(shù)大的參數(shù)一般不太敏感。

[1] 黃國如,朱慶平,馬駿,等.水位演算模型及其在水位預(yù)報中的應(yīng)用[J].水文,1999(2):1-6.

[2] 包為民.水文預(yù)報[M].北京:中國水利水電出版社,2009:109-112.

[3]CUNGEJ,HOLLYF,VERWAYA.PracticalAspectsofComputationRiverHydraulics[M].London:PitmanPublishingLimited, 1990.

[4] 葛守西.現(xiàn)代洪水預(yù)報技術(shù)[M].北京:中國水電水電出版社,1999:1-23.

[5]PHIENH,KHAN.FloodforecastingfortheupperreachoftheRedRiverBasin,NorthVietnam[J].WaterSA,2003,29(3):267-272

[6]DOOGEJ.LinearTheoryofHydrologicSystems[R].USDeptAgrTechBull1948SuperintendentofDocuments.Washington,1973:327.

[7] 牟金玲,狄娟.黑龍江干流中上游段水文特性[J].黑龍江水利科技,2007,35(2):101-102.

[8] 胡煒,董增川.統(tǒng)計相關(guān)分析在太湖流域水位預(yù)報中的應(yīng)用[J].水利經(jīng)濟(jì),2008,26(3):55-58.

[9] 張顯揚,董增川,王建群,等.雅魯藏布江尼洋河流域洪水預(yù)報方法研究[J].河海大學(xué)學(xué)報(自然科學(xué)版), 2005,33(5):530-533.

[10] 梁忠民,欒承梅,李致家.線性動態(tài)系統(tǒng)模型在實時洪水預(yù)報中的應(yīng)用[J].水力發(fā)電,2003,29(6):16-18,25.

[11]HAANC.Parametricuncertaintyinhydrologicmodeling[J].TransactionofASAS,1989,32(1):137-146.

[12]ANNANJ.Modelingunderuncertainty:MonteCarlomethodsfortemporallyvaryingparameters[J].EcologicalModeling,2001,136(2-3):297-302.

[13] 裴洪平,汪勇.杭州西湖富營養(yǎng)化模型的不確定性分析[J].生物數(shù)學(xué)學(xué)報,2004,19(1):117-122.

[14]HANNAS.MonteCarloestimatesofuncertaintiesinpredictionsbyaphotochemicalgridmodel(UAM-IV)duetouncertaintiesininputvariables[J].AtmosphericEnvironment,1997,32(21):3619-3628.

[15] 邢可霞,郭懷成.環(huán)境模型不確定性分析方法綜述[J].環(huán)境科學(xué)與技術(shù),2006,29(5):112-114.

[16]SL250-2000水文情報預(yù)報規(guī)范[S].北京:中國水利水電出版社,2000.

The method of water-level forecasting for the Heilongjiang mainstream

ZHAOSiyuan1，HAOZhenchun1，LIUWenbin2

(1.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,Hohai University, Nanjing 210098, China；2.Hydrological Bureau of Heilongjiang Province, Harbin 150001, China)

The black-box model based on multivariate hybrid linear regression is applied for water-level simulation of the upper and middle reaches of the Heilongjiang mainstream, including the Mohe station, the Oupu station, the Sandaoka station, the Changfatun station, the Wuyun station, and the Qindeli station. The model structure is determined by AIC and DW method; parameters of the model are estimated by the least square method. The models are calibrated with day time-scale hydrologic data of 2006-2010 in the flood season(May to September), and verified with day time-scale hydrologic data in the flood season(May to September) in 2011. And the model’s uncertainty and parameters sensivitity are analyzed. The results are as follows: ① The model is applied well in the Heilongjiang mainstream and accurate within acceptable range. Forecast schemes’ pass rate reaches more than 85% and coefficients of determination are above 0.97 for six stations. The model can be used for water-level forecasting. ② The paper quantifies uncertainty of the Multivariate hybrid linear regression model results using the Monte Carlo methods. The larger the standard deviation, the more uncertain the flood hydrograph forecasting, and thus model is more uncertain. The parameters’ sensitivity analysis based on standardized regression coefficients shows that the parameters with minimum values in time order of mainstream are the most sensitive, the parameters with minimum values in time order of tributary take the second place, and the parameters with higher values in time order are less sensitive.

water-level forecasting; mixed linear regression; Heilongjiang mainstream; uncertainty analysis; Monte Carlo

10.13471/j.cnki.acta.snus.2016.02.25

2015-05-13

國家自然科學(xué)基金資助項目(41371047)；中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項基金資助項目(XDA05110102)；水文水資源與水利工程科學(xué)國家重點實驗室專項經(jīng)費資助項目(1069-514031112)

趙思遠(yuǎn)(1991年生)，男；研究方向：水文物理規(guī)律模擬及水文預(yù)報；通訊作者：郝振純;E-mail:hzchun@hhu.edu.cn

P338

A

0529-6579(2016)02-0151-08