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