馮曉蕾,何錫君,胡孫堅
(1.浙江水文新技術(shù)開發(fā)經(jīng)營公司,浙江 杭州 310009;2.浙江省水文局,浙江 杭州 310009;3.浙江省工業(yè)環(huán)保設(shè)計研究院有限公司,浙江 杭州 310000)
水文預報就是根據(jù)前期或現(xiàn)時已出現(xiàn)的水文、氣象等信息,運用水文學、氣象學、水力學的原理和方法,對河流、湖泊等水體未來一定時段內(nèi)的水文情勢做出定量或定性的預報[1]。按其預見期的長短,可分為短期水文預報與中長期水文預報[2]。中長期水文預報的定義為[3]:根據(jù)前期水文氣象要素,用成因分析與數(shù)理統(tǒng)計的方法,對未來較長時間的水文要素進行科學的預測,其預見期一般為3 d ~ l a。
本文以衢州水文站的年平均徑流量為預報對象,22項氣象環(huán)流因子為預報因子,分析對衢州水文站年平均徑流量有重要影響的因子。使用逐步回歸分析模型,對各年的平均流量分別進行擬合、預報,以期為流域水資源的合理調(diào)度提供重要依據(jù)。
回歸分析法是從分析影響預報對象的因子中挑選出一批預報因子,建立其統(tǒng)計規(guī)律作為預報的依據(jù)。為了預報某對象(如某一水文要素)未來時刻的變化,選擇預報量前期已發(fā)生的多個有關(guān)的水文、氣象要素或其他物理要素,即預報因子為研究對象,利用回歸方法去分析多個預報因子與這個預報量之間的關(guān)系,建立它們統(tǒng)計關(guān)系的方程式。從而,對預報對象做出估計,其主要步驟為:① 確定預報量并選擇適當?shù)囊蜃樱虎诟鶕?jù)數(shù)據(jù)計算回歸系數(shù)標準方程組所包含的有關(guān)統(tǒng)計量;③解方程組,定出回歸系數(shù);④建立回歸方程進行統(tǒng)計顯著性檢驗;⑤利用已出現(xiàn)的因子值代入回歸方程做出預報量的估計。
利用統(tǒng)計方法來挑選預報因子,統(tǒng)計方法中最常用的就是相關(guān)分析,通過計算2個序列y和x之間的線性相關(guān)系數(shù),從而判斷它們之間的線性相關(guān)程度。預報因子與預報對象的相關(guān)性計算有很多種方法,本文中采用單相關(guān)計算的方法。相關(guān)系數(shù)的計算公式為:
式中:xi、yi分別為影響因子和待預報水文要素的實測值;x、y 分別為影響因子和待預報水文要素的多年平均值;r為樣本的相關(guān)系數(shù),其顯著程度常用檢驗來說明:
式中:r為樣本的相關(guān)系數(shù);n為樣本容量(個)。
選定信度α后可從t分布表中查出相應的tα,當t≥tα時,則認為在這一信度下兩者是線性相關(guān)的,否則認為是不相關(guān)的。
2.2.1 基本思路
先定義一個衡量因子對預報對象重要性的指標,以便從中挑選出對其影響顯著的因子。因子的挑選應逐步進行,在建立回歸方程的過程中,每一步只挑選出一個因子,要求當步選出的因子能滿足殘差平方和下降最多并且能通過指定的信度檢驗。直至還未引入回歸方程的因子中不存在對作用顯著的因子為止。隨著因子的引入,由于因子之間的相互配合關(guān)系,可能產(chǎn)生當后面的因子引入以后會引起前面的因子對的作用顯著變小,甚至不顯著,要將不顯著的因子剔除。因此,在逐步回歸中每一步都要做剔除和引進因子的檢驗。直至方程中既不能引入也不能剔除因子為止。最后的方程中只包含了對預報對象影響顯著的因子。而沒有進入方程的因子,添加任何一個,都不會對回歸效果有顯著的改進。
2.2.2 計算步驟
2.2.2.1 建立標準化的正規(guī)方程組
假設(shè)預報因子X是由p個組成,每個因子的長度為n;預報對象y,其長度為n。在逐步回歸中采用的是標準化的正規(guī)回歸方程組:
其系數(shù)與常數(shù)項組成的零步增廣矩陣為:
記作R(0)= Rij(0)。
2.2.2.2 因子的剔除和引入
假設(shè)逐步回歸已經(jīng)進行了l步,回歸方程引入了l個預報因子,則第l + 1步的計算步驟如下:
(1)計算全部待選因子的方差貢獻:
式中:為第l步時矩陣R(l)中的元素,表示已引入因子的方差貢獻,表示未引入因子的方差貢獻。
(2)剔除因子。在已引入的l個因子中挑選出方差貢獻最小者,記為min{Vk(l)},其對應的因子為xk,定義方差比為:
給定顯著性水平α,在F分布表中查出自由度為1,(n - l - 1)時的上分位點Fα,若Fk≤Fα,則認為因子xk的方差貢獻不顯著,可剔除;否則不剔除。
(3)引入因子。在為引進方程的因子中挑選出方差貢獻最大者,記為max{},其對應的因子為xk,定義方差比為:
式中:n為樣本容量;l為當前因子數(shù)。
給定顯著性水平α,在F分布表中查出自由度為1,(n - l - 2)時的上分位點Fα,若Fk≥Fα,則引進因子xk,否則不引進。
在逐步回歸建模過程中,因子xk的剔除和引進都需要進行消去變換,消去變換的方法為:
式中: i表示矩陣的行數(shù);j表示矩陣的列數(shù),k是被選入的因子的序號。
2.2.2.3 計算結(jié)果分析
循環(huán)因子引進和剔除步驟,若可供挑選的因子全被引進,或者既無因子可引進又無因子可剔除時,則終止逐步回歸分析。
假設(shè)最后共引入了m個因子進入回歸方程,則應用逐步回歸分析法建立的標準化回歸方程:
根據(jù)標準化前后數(shù)據(jù)之間的關(guān)系,即:
則標準化回歸方程可轉(zhuǎn)化為下式:
式中:m為因子數(shù)(個);b0,b1,…bm為回歸系數(shù)。
衢州水文站屬國家級重點水文站,是錢塘江上游的一個重要控制站,流域?qū)偕降氐匦?,多山溪河流,坡陡流急,容易誘發(fā)洪澇災害。流域內(nèi)降雨充沛,多年平均降雨量1 815.2 mm,但年內(nèi)分配極不均勻,降雨主要集中在3 —6月,占全年總量的57%,特別是梅雨季節(jié),常常出現(xiàn)大暴雨、特大暴雨天氣。據(jù)實際資料統(tǒng)計,流域平均匯流時間為8 h左右。衢州水文站的水文預報,不僅關(guān)系到浙贛鐵路的安全,還關(guān)系到碗窯、白水坑、湖南鎮(zhèn)、銅山源、黃壇口等大型水庫的錯峰調(diào)度,在錢塘江流域防洪體系中,處于非常重要的地位。
以衢州水文站年平均流量為預報對象,模型編制依據(jù)的水文資料為衢州水文站1951 — 1999年的年平均流量。
根據(jù)各參數(shù)物理特性,對原始資料及衢州水文站資料進行合理性、一致性分析,數(shù)據(jù)合理,資料系列完整,無需插補延長,可直接用于計算。由于衢州水文站從1951年起有完整的水文資料,故資料系列取1951 — 1994年共44 a同步資料進行預報分析及擬合檢驗;1995 — 1999年共5 a同步資料進行預報方程預報檢驗。
根據(jù)前述預報因子分析方法,采用單相關(guān)系數(shù)法計算衢州水文站年平均流量與22個預報因子之間的相關(guān)系數(shù),并對求得的相關(guān)系數(shù)進行t檢驗,篩選出與預報對象相關(guān)性強的因子。
具體計算時,取衢州水文站1952 — 1994年的長系列年平均流量資料,與某一預報因子1951 — 1993年的系列資料分別計算相關(guān)系數(shù),即在提前1 a的范圍內(nèi)進行挑選,在滿足r>rα的基礎(chǔ)上,選擇相關(guān)系數(shù)最大的系列為該因子相關(guān)月份系列選?。蝗绻A報因子1 — 12月系列與預報對象相關(guān)系數(shù)都小于rα,則表明該因子與預報對象線性相關(guān)不顯著,不予引進;依次計算該預報對象與22個預報因子1 — 12月系列的相關(guān)系數(shù),即可挑選出相應的預報因子。
由可信度(α = 0.05)及樣本數(shù)(n = 43),查表得最低相關(guān)系數(shù)rα = 0.301;由可信度(α = 0.01)及樣本數(shù)(n = 43),查表得最低相關(guān)系數(shù)rα = 0.389,具體預報因子挑選結(jié)果見表1。
表1 預報因子及相關(guān)系數(shù)表
首先,依次挑選出該預報對象與22個預報因子系列中相關(guān)系數(shù)最大的系列,再利用預報因子及衢州水文站年平均流量實測資料,根據(jù)上文逐步回歸模型原理,通過程序計算方程的各項系數(shù)。
逐步回歸每一步的回歸方程系數(shù),過程一共進行了11步,從最后第11步的計算結(jié)果可知:22個因子最終有9個因子進入方程中,分別是因子7、因子2、因子11、因子4、因子10、因子19、因子21、因子20及因子16?!胺菢藴驶貧w系數(shù)”中的“B”列系數(shù)代入,得衢州水文站年平均流量的預報方程:
y = 1 049.377 + 1.221X1- 0.727X2- 2.117X3- 4.439X4+0.451X5+ 1.618X6- 17.557X7+ 5.894X8- 4.013X9(11)
式中:y為衢州水文站年平均流量(m3/s);X1為上一年 12 月 H 500 mb(25 ~ 35 N、110 ~ 130 E)長江中下游區(qū)7點合計;X2為上一年3月 H 500 mb(50 ~ 55 N、70 ~ 90 E + 40 ~ 45 N + 65 ~ 85 E)巴爾喀什湖區(qū);X3為上一年10月H 500 mb(30 ~ 40 N、80 ~ 90 E)西安高原子6點合計;X4為上一年1月H 500 mb(120 E、20 ~ 40 N)高度差(沿120 E線20 ~ 40 N);X5為上一年4月H 500 mb(50 ~ 60 N、100 ~ 120 E)貝加爾湖區(qū)8點合計;X6為上一年10月鄂海平均高度H(135 ~ 150 E、45 ~ 60 N);X7為上一年1月C102 102站西風指數(shù)(115 E、25 ~ 30 N);X8為上一年1月西風風速V27.5 N(m/s)(105 E、25 ~ 35 N);X9為上一年8月西風風速V42.5 N(m/s)(105 E、40.5 N)。
(1)復相關(guān)系數(shù)(R)。復相關(guān)系數(shù)為:
式中:R為復相關(guān)系數(shù);U為回歸平方和;Q為殘差平方和;Syy為總和平方和。
(2)剩余標準差(Sy)。剩余標準差為:
式中:Sy為剩余標準差;n為資料年限(a);Q為殘差平方和; m為挑選的因子(個)。
(3) 回歸效果的F檢驗?;貧w效果的F檢驗為:
由計算程序求得上述回歸方程的復相關(guān)系數(shù)R、剩余標準差Sy和F檢驗值。
根據(jù)本例樣本數(shù)n = 43,因子個數(shù)m = 9,在信度α =0.01的條件下,查復相關(guān)系數(shù)表有臨界值Rα= 0.671,而衢州水文站逐步回歸方程的R = 0.852,R>Rα。因此在α =0.01的水平下回歸的效果是顯著的。該方程剩余標準差Sy=34.04 m3/s。
根據(jù)本例fU = 9,fQ = 33,在信度α = 0.01的條件下,查F分布表得臨界值Fα= 3.01,而衢州水文站逐步回歸方程的F = 9.72,F(xiàn)>Fα。因此在α = 0.01的水平下回歸效果是顯著的。
根據(jù)已建立的衢州水文站年均流量預報模型所用的1952 — 1994年共43 a的實測年平均流量對預報模型進行歷史擬合檢驗,計算相應年份的預報誤差及許可誤差,其中許可誤差采用實測流量的20%來計算,如預報誤差<許可誤差,則該年預報值合格,否則為不合格。該預報方程在43 a的平均流量歷史擬合中預報合格次數(shù)為37次,合格率為86.0%,根據(jù)GB 22482 — 2008 — T《水文情報預報規(guī)范》規(guī)定,中長期預報合格率≥85%時,預報方程等級為甲等。具體擬合結(jié)果見圖1和表2。
圖1 衢州水文站年平均流量逐步回歸方程擬合圖
表2 衢州水文站年平均流量逐步回歸方程擬合表
續(xù)表2
根據(jù)已建立的預報方程,以1995 — 1999年共5 a的實際年平均流量資料對預報模型進行預報檢驗,合格總數(shù)為5,合格率為100%。具體預報結(jié)果見表3。
表3 衢州水文站年平均流量逐步回歸方程預報檢驗表
由表3可知,年平均流量回歸方程預報檢驗效果較好,達到了預報精度的要求。
(1)本文以衢州水文站的年平均徑流量為預報對象,22項氣象環(huán)流因子為預報因子,分析對衢州水文站年平均徑流量有重要影響的因子,使用逐步回歸分析模型,對衢州水文站各年的平均流量進行擬合、預報。結(jié)果表明,逐步回歸模型的擬合、預報效果達到預期精度要求,可用于衢州水文站年均流量的預報。
(2)使用逐步回歸預報方程在43 a的平均流量歷史擬合中預報合格次數(shù)為37次,合格率為86.0%,根據(jù)GB 22482 — 2008 — T《水文情報預報規(guī)范》規(guī)定,預報方程等級為甲等預報方案。預報檢驗合格率為100%。
(3)采用氣象因子與水文要素建立的中長期水文預報模型,能夠獲得較高的預報精度及較長的預見期,可以為地方政府及防汛部門提前提供決策依據(jù),有效地在防災減災中降低人民的生命財產(chǎn)損失,也可以對水資源進行優(yōu)化配置,合理調(diào)配。
參考文獻:
[1]中華人民共和國水利部.水文基本術(shù)語和符號標準:GB/T 50095 — 98[S].北京:中國計劃出版社,1999.
[2] 詹道江,葉守澤.工程水文學[M].北京:中國水利水電出版社,2000.
[3] 范鐘秀.中長期水文預報[M].南京:河海大學出版社,1999.