趙然杭,甘 甜,逄曉騰,王興菊,茍偉娜,齊 真
(1.山東大學(xué)土建與水利學(xué)院,濟(jì)南250100;2.青島市供水事業(yè)發(fā)展中心,山東青島266071)
降雨是水文過程的重要環(huán)節(jié),能直接反映區(qū)域氣候變化及水文過程,降雨量變化直接關(guān)乎區(qū)域水資源開發(fā)利用乃至影響社會(huì)經(jīng)濟(jì)發(fā)展。近年來,由于氣候變化與人類活動(dòng)的影響,生態(tài)問題嚴(yán)峻、洪旱災(zāi)害頻發(fā),充分挖掘降雨時(shí)間序列所含信息能更好揭示其自然規(guī)律,為水資源的規(guī)劃管理提供科學(xué)的指導(dǎo),促進(jìn)人與自然和諧共處。對降雨量數(shù)據(jù)準(zhǔn)確預(yù)測能為防洪抗旱工作提供科學(xué)的建議,進(jìn)而改善居民生活質(zhì)量,促進(jìn)區(qū)域政治與經(jīng)濟(jì)的發(fā)展,是事關(guān)民生的大事。
數(shù)據(jù)挖掘是指從大量的、有噪聲的、模糊的、隨機(jī)的數(shù)據(jù)中,提取隱含其中的有用的信息和知識的過程[1]。國內(nèi)外學(xué)者對降雨時(shí)間序列數(shù)據(jù)挖掘進(jìn)行了相關(guān)研究:李宜軒運(yùn)用線性擬合的方法,綜合分析了桓臺縣1981-2018年的降水趨勢變化情況,結(jié)論與實(shí)際情況相符[2];衡彤等采用小波變換方法,對新安江流域黃山地區(qū)降水量時(shí)間序列的多時(shí)間尺度變化及突變特征進(jìn)行了探討,研究分析了主汛期降水與年降水的主周期,兩者較為接近,結(jié)果合理[3];李平蘭等分析了1970-2015年會(huì)東縣的降水變化特征,采用了累積距平法與Mann-Kendall 法的突變點(diǎn)綜合判定方法,彌補(bǔ)了單個(gè)方法檢驗(yàn)?zāi)芰Φ牟蛔悖?];Moravej等基于線性時(shí)間序列模型,利用ADF 函數(shù)、Mann-Kendall 檢驗(yàn),對伊朗West Azerbaijan 省的降雨時(shí)間序列進(jìn)行了趨勢分析,模型結(jié)論準(zhǔn)確且適用性較強(qiáng)[5]。
國內(nèi)外學(xué)者對降雨時(shí)間序列數(shù)據(jù)預(yù)測進(jìn)行了相關(guān)研究:王曉鵬等應(yīng)用基于Box-Jenkins方法的時(shí)間序列分析技術(shù),針對青南高原的四個(gè)典型地區(qū)1961-2005年降水量序列,建立了自回歸滑動(dòng)平均模型對未來降水量進(jìn)行了預(yù)測,模型因能反復(fù)識別修改,適用于復(fù)雜情況的數(shù)據(jù)擬合與預(yù)測[6];賈海峰等提出了灰色-時(shí)序組合預(yù)測模型,能對年降雨量進(jìn)行中長期預(yù)測,并以北京市大興縣黃村氣象站35年降雨量資料為例進(jìn)行了驗(yàn)證,模型能滿足精度要求,且對其他既具有擺動(dòng)又有趨勢數(shù)據(jù)的模擬和預(yù)測同樣適用[7];高瑞華利用煙臺地區(qū)1981-2010年降水資料,建立了干旱的灰色預(yù)測GM(1,1)模型,對2010-2030年的干旱災(zāi)害進(jìn)行了預(yù)測,模型通過了殘差與后驗(yàn)差檢驗(yàn)分析,在短期預(yù)測時(shí)精度能滿足要求,但不適用于長期預(yù)測[8];Kolachian 等采用分位數(shù)映射與平均貝葉斯模型,通過提高觀測值與集合預(yù)報(bào)平均值之間的相關(guān)系數(shù),對歐洲中心2020年各月份不同降水模式的若干氣象站的觀測降水序列進(jìn)行了預(yù)報(bào),結(jié)果較為準(zhǔn)確[9]。
以上研究分別對降雨時(shí)間序列進(jìn)行了數(shù)據(jù)挖掘或預(yù)測,沒有將數(shù)據(jù)挖掘與預(yù)測建立緊密聯(lián)系,而數(shù)據(jù)挖掘能為數(shù)據(jù)預(yù)測提供信息支撐,數(shù)據(jù)預(yù)測結(jié)果也能驗(yàn)證挖掘的準(zhǔn)確性,兩者相輔相成。由于降雨時(shí)間序列受趨勢、周期與突變等因素影響,復(fù)雜多變,需要在充分挖掘數(shù)據(jù)信息的基礎(chǔ)上,進(jìn)行降雨時(shí)間序列的預(yù)測研究。因此開展基于時(shí)間序列分解的降雨數(shù)據(jù)挖掘與預(yù)測研究,為雨(洪)水資源管理工作提供決策技術(shù)支持。
降雨時(shí)間序列因受多種因素影響,存在著趨勢性、周期性、突變性等多種變化特征,將時(shí)間序列分解為趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)與隨機(jī)項(xiàng)并對各項(xiàng)進(jìn)行分析與研究能更好揭示降雨量數(shù)據(jù)的變化規(guī)律。設(shè)X為時(shí)間序列,T為趨勢項(xiàng)、P為周期項(xiàng)、B為突變項(xiàng)、R為隨機(jī)項(xiàng),則X可以表示為:
常見的趨勢項(xiàng)分析法包括:累積距平[10]、Mann-Kendall 趨勢分析[11]、Hurst 指數(shù)[12]、特征點(diǎn)[13],周期項(xiàng)分析法包括:小波分析[14],突變項(xiàng)分析法包括:Mann-Kendall 突變檢驗(yàn)[15]、Pettitt 突變檢驗(yàn)[16],隨機(jī)項(xiàng)檢驗(yàn)法包括:自相關(guān)圖和單位根[17,18]。
提出了基于時(shí)間序列分解的降雨數(shù)據(jù)挖掘-預(yù)測模型(Data Mining And Forecasting model,DMAF 模型),不僅能挖掘出降雨時(shí)間序列的內(nèi)在信息,還能對其進(jìn)行準(zhǔn)確預(yù)測。
(1)數(shù)據(jù)挖掘方面:①數(shù)據(jù)分解:對于降雨時(shí)間序列X,采用特征點(diǎn)法分解其趨勢項(xiàng)T;小波分析法分解其周期項(xiàng)P;Mann-Kendall法、Pettitt法進(jìn)行突變性分析,無特殊情況突變項(xiàng)B不予考慮;計(jì)算余項(xiàng)X-T-P-B得隨機(jī)項(xiàng)R。②數(shù)據(jù)挖掘:分別采用累積據(jù)平法、Mann-Kendall趨勢分析法、Hurst指數(shù)法、特征點(diǎn)法進(jìn)行趨勢分析,采用小波分析法進(jìn)行周期分析,采用Mann-Kendall 法、Pettitt 法進(jìn)行突變性分析,采用自相關(guān)圖法、單位根法對隨機(jī)項(xiàng)進(jìn)行檢驗(yàn)。
(2)數(shù)據(jù)預(yù)測方面:利用線性延長方法將趨勢項(xiàng)T基于特征點(diǎn)法的末端直線進(jìn)行延長得到Ty,利用Matlab 小波工具箱Extend 函數(shù)將周期項(xiàng)P進(jìn)行延長得到Py,無特殊情況下不考慮突變項(xiàng),采用NAR 神經(jīng)網(wǎng)絡(luò)法對隨機(jī)項(xiàng)進(jìn)行預(yù)測。然后,將延長后的趨勢項(xiàng)Ty、周期項(xiàng)Py與預(yù)測所得隨機(jī)項(xiàng)Ry相加得預(yù)測結(jié)果,流程如圖1 所示。為提高預(yù)測精度,將NAR 神經(jīng)網(wǎng)絡(luò)經(jīng)遺傳算法進(jìn)行優(yōu)化,其流程圖如圖2。
圖1 研究流程圖Fig.1 Research flow chart
圖2 基于遺傳算法的神經(jīng)網(wǎng)絡(luò)優(yōu)化流程圖Fig.2 Neural network optimization flow chart based on genetic algorithm
桓臺縣是淄博市下轄縣,位處魯中山區(qū)和魯北平原的結(jié)合地帶,人口50.15 萬,縣域面積498.25 km2,屬暖溫帶季風(fēng)氣候[19]?;概_縣河網(wǎng)較為密布,烏河、東西豬龍河、孝婦河、澇滯河、小清河、預(yù)備河等10余條河流流經(jīng)縣內(nèi),但因處于魯山北麓山前洪沖積平原和黃泛平原迭交地帶,南受魯中山區(qū)洪水沖積,北受黃河泛濫淤淀,形成湖洼沼地,地形地貌復(fù)雜,降雨量預(yù)測難度大。加之季風(fēng)氣候的影響,降雨年際及年內(nèi)分布不均,旱澇災(zāi)害頻發(fā),嚴(yán)重制約著當(dāng)?shù)氐慕?jīng)濟(jì)發(fā)展。因此選取桓臺縣為研究區(qū)域,數(shù)據(jù)來源于桓臺縣水文站。
利用matlab 進(jìn)行編程,以桓臺縣1979-2018年降雨數(shù)據(jù)為例,時(shí)間步長為月,將時(shí)間序列分解為趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)與隨機(jī)項(xiàng)進(jìn)行數(shù)據(jù)挖掘。
(1)趨勢項(xiàng):分別采用累積距平、Mann-Kendall、Hurst 指數(shù)等方法進(jìn)行趨勢性分析,特征點(diǎn)法進(jìn)行趨勢項(xiàng)提取。累積距平圖如圖3所示,Hurst指數(shù)圖如圖4所示,Mann-Kendall趨勢分析結(jié)果如表1所示。
由圖3可得,累積距平圖總體較為平穩(wěn),存在先降后升的趨勢,但不明顯。由表1 可得,數(shù)據(jù)總體呈現(xiàn)微弱上升趨勢,但趨勢不顯著。圖4 中,Hurst 指數(shù)為0.067 7,H值<0.5,表明該時(shí)間序列具有長期相關(guān)性,但未來總體趨勢與過去相反。
表1 Mann-Kendall 趨勢檢驗(yàn)結(jié)果Tab.1 Mann Kendall trend test results
圖3 累積距平圖Fig.3 Cumulative anomaly
圖4 Hurst指數(shù)圖Fig.4 Hurst index
綜合分析可得,桓臺縣1979-2018年月降雨數(shù)據(jù)趨勢變化較為平穩(wěn),存在微弱上升趨勢,預(yù)測下一時(shí)段將呈現(xiàn)微弱下降趨勢。
以上為定性分析,現(xiàn)根據(jù)特征點(diǎn)法定量分解趨勢項(xiàng)。取閾值C=12,將雨量數(shù)據(jù)進(jìn)行分段,分段后的特征點(diǎn)利用插值法首尾相連,結(jié)果如圖5所示。
圖5 特征點(diǎn)分段圖Fig.5 Feature point segmentation
由圖5可得,桓臺縣1979-2018年的月降雨量大致呈現(xiàn)5次豐枯交替現(xiàn)象,為便于直觀描述,將各特征點(diǎn)的取值前30%標(biāo)注為high,中間40%標(biāo)注為mid,后30%標(biāo)注為low,如圖6所示。對各段斜率進(jìn)行描述,斜率小于-1標(biāo)注為驟減,-1~-0.1標(biāo)注為下降,-0.1~0.1 標(biāo)注為平穩(wěn),0.1~1 標(biāo)注為增加,大于1 標(biāo)注為激增。如圖7所示。
圖6 特征點(diǎn)取值描述圖Fig.6 Description of feature point value
圖7 各段斜率描述圖Fig.7 Slope description of each segment
經(jīng)計(jì)算插值法誤差為79.50,用回歸法連接后效果如圖8所示(圖8中紅色數(shù)字為各段回歸系數(shù))。
圖8 趨勢項(xiàng)分解圖Fig.8 Trend term decomposition
插值法更為直觀,便于描述序列變化情況?;貧w法誤差為53.53,明顯小于插值法,因此,選用回歸法的結(jié)果為時(shí)間序列的趨勢項(xiàng)。
(2)周期項(xiàng):使用Matlab 中的小波工具箱得小波方差圖,如圖9所示。由圖9可以看出,桓臺縣1979-2018年月降雨量的第一主周期為19(月)。以第一主周期時(shí)間尺度繪制小波系數(shù)圖,如圖10所示。選取第一主周期小波系數(shù)為周期項(xiàng)。
圖9 小波方差圖Fig.9 Wavelet variance
圖10 第一主周期小波系數(shù)圖Fig.10 First principal period wavelet coefficients
(3)突變項(xiàng):分別采用Mann-Kendall 法、Pettitt 法進(jìn)行突變性檢驗(yàn),結(jié)果如圖11、表2所示:
由圖11 可得,1979年-2018年的月降雨量突變情況發(fā)生少。表2中,突變點(diǎn)檢驗(yàn)結(jié)果不顯著,不能認(rèn)為產(chǎn)生突變。為避免單個(gè)檢驗(yàn)方法決策失誤,采用文獻(xiàn)[4]中突變點(diǎn)判別方法,認(rèn)為在兩種判別方法中均為突變點(diǎn)的月份才發(fā)生顯著突變。綜合分析,認(rèn)為1979-2018年月降雨不發(fā)生顯著突變情況,不進(jìn)行突變項(xiàng)分解。
表2 Pettitt法檢驗(yàn)結(jié)果Tab.2 Results of Pettitt method
圖11 Mann-Kendall突變檢驗(yàn)圖Fig.11 Mann-Kendall mutation test
(4)隨機(jī)項(xiàng):將時(shí)間序列減去趨勢項(xiàng)、周期項(xiàng)與突變項(xiàng)可得隨機(jī)項(xiàng)。若隨機(jī)項(xiàng)的自相關(guān)性較小且具有平穩(wěn)性,說明趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)信息被充分提取。利用Eviews 軟件對隨機(jī)項(xiàng)進(jìn)行自相關(guān)(AC)、偏自相關(guān)(PAC)檢驗(yàn),結(jié)果如圖12所示。
圖12 自相關(guān)圖Fig.12 Autocorrelation diagram
由圖12得,AC值與PAC值始終圍繞零軸小范圍波動(dòng),進(jìn)一步進(jìn)行單位根檢驗(yàn),得P值 為0、AIC為10.418 02、SC為10.54212、LogL為-2 423.818。由自相關(guān)與單位根檢驗(yàn)結(jié)果可知,隨機(jī)項(xiàng)不存在單位根,為平穩(wěn)時(shí)間序列。趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)信息已被較為充分提取。
2.3.1 隨機(jī)項(xiàng)預(yù)測
為便于驗(yàn)證預(yù)測結(jié)果,以1979-2014年432 組月降雨量數(shù)據(jù)為率定數(shù)據(jù),2015-2016年月降雨量數(shù)據(jù)為驗(yàn)證數(shù)據(jù),對2017-2018年月降雨數(shù)據(jù)進(jìn)行預(yù)測。
隨機(jī)項(xiàng)預(yù)測分別采用NAR和NARX[20]神經(jīng)網(wǎng)絡(luò)預(yù)測法。
(1)NAR 神經(jīng)網(wǎng)絡(luò)預(yù)測:構(gòu)建經(jīng)遺傳算法優(yōu)化的NAR 神經(jīng)網(wǎng)絡(luò),經(jīng)率定,滯后階數(shù)取18,中間層個(gè)數(shù)取8,訓(xùn)練集、測試集和驗(yàn)證集比例取70∶15∶15。所構(gòu)建神經(jīng)網(wǎng)絡(luò)模型的回歸系數(shù)圖、誤差圖、預(yù)測結(jié)果圖如圖13~15所示。
圖13 神經(jīng)網(wǎng)絡(luò)回歸系數(shù)圖Fig.13 Regression coefficient diagram of neural network
由圖13 知,神經(jīng)網(wǎng)絡(luò)模型回歸系數(shù)較高,總回歸系數(shù)0.793 26,訓(xùn)練集回歸系數(shù)0.815 96。由圖14 知,誤差集中分布于0附近,多為小誤差,大誤差較少。
圖14 神經(jīng)網(wǎng)絡(luò)誤差圖Fig.14 Neural network error chart
(2)NARX 神經(jīng)網(wǎng)絡(luò)預(yù)測:為了驗(yàn)證DMAF模型預(yù)測結(jié)果的合理性,采用NARX神經(jīng)網(wǎng)絡(luò)隨機(jī)項(xiàng)預(yù)測法,結(jié)果見圖16。
2.3.2 時(shí)間序列的總體預(yù)測
將圖15 和圖16 的預(yù)測結(jié)果,分別與所對應(yīng)的趨勢項(xiàng)和周期項(xiàng)相加,得到基于NAR 和NARX 神經(jīng)網(wǎng)絡(luò)的2017-2018年月降雨的預(yù)測結(jié)果,見圖17。其中,模型預(yù)測誤差16.79%,NARXANN預(yù)測誤差20.05%,直接預(yù)測誤差41.62%。
圖15 基于NAR神經(jīng)網(wǎng)絡(luò)隨機(jī)項(xiàng)預(yù)測結(jié)果圖Fig.15 Random term prediction result chart based on NAR neural network
圖16 基于NARX神經(jīng)網(wǎng)絡(luò)隨機(jī)項(xiàng)預(yù)測結(jié)果圖Fig.16 Random term prediction result chart based on NARX neural network
圖17 預(yù)測結(jié)果對比圖Fig.17 Comparison chart of prediction results
(1)時(shí)間序列挖掘:由圖3和表1可知,桓臺縣1979-2018年月降雨量數(shù)據(jù)有微弱的上升趨勢,此結(jié)論與文獻(xiàn)[2]相符。由圖3 可預(yù)測未來將呈微弱下降趨勢,由圖9 可知,其第一主周期是19(月),由圖11與表2可知,數(shù)據(jù)不存在明顯的突變情況。
(2)時(shí)間序列預(yù)測精度:由圖17 可知,與1979-2018年月降雨量的實(shí)測數(shù)據(jù)相比,基于NAR 和NARX 神經(jīng)網(wǎng)絡(luò)的預(yù)測誤差分別為16.79%和20.05%。利用NAR 神經(jīng)網(wǎng)絡(luò)直接預(yù)測法進(jìn)行時(shí)間序列的總體預(yù)測(結(jié)果見圖17),預(yù)測結(jié)果誤差較大(為41.62%)。因此,上述3 種預(yù)測方法,基于NAR 神經(jīng)網(wǎng)絡(luò)隨機(jī)項(xiàng)預(yù)測方法獲得月降雨數(shù)據(jù)的精度最高。
(1)關(guān)于降雨時(shí)間序列的數(shù)據(jù)挖掘與預(yù)測的關(guān)系:降雨時(shí)間序列受趨勢、周期與突變等因素影響,復(fù)雜多變,數(shù)據(jù)挖掘能為數(shù)據(jù)預(yù)測提供信息支撐,數(shù)據(jù)預(yù)測結(jié)果也能驗(yàn)證挖掘信息的準(zhǔn)確性。而現(xiàn)有水文數(shù)據(jù)挖掘與預(yù)測方面的文獻(xiàn)[2-5],未見有基于水文數(shù)據(jù)挖掘的預(yù)測研究,即沒有建立數(shù)據(jù)挖掘與預(yù)測的緊密聯(lián)系,而是直接預(yù)測。創(chuàng)新性地提出了DMAF 模型,將降雨時(shí)間序列分解為趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)及隨機(jī)項(xiàng),并對各項(xiàng)進(jìn)行數(shù)據(jù)挖掘,同時(shí)建立NAR 神經(jīng)網(wǎng)絡(luò)模型對隨機(jī)項(xiàng)進(jìn)行預(yù)測,再將預(yù)測結(jié)果與數(shù)據(jù)挖掘后的趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)相加得降雨時(shí)間序列預(yù)測結(jié)果,相較于直接預(yù)測法,結(jié)果更為準(zhǔn)確。因此DMAF 模型能在充分挖掘降雨時(shí)間序列的基礎(chǔ)上,對其進(jìn)行準(zhǔn)確預(yù)測。
(2)關(guān)于利用DMAF 模型對降雨時(shí)間序列長短期預(yù)測問題:對研究區(qū)域的降雨時(shí)間序列進(jìn)行分解后,利用線性延長方法將趨勢項(xiàng)T基于特征點(diǎn)法的末端直線進(jìn)行延長得到Ty;利用Matlab 小波工具箱Extend 函數(shù)將周期項(xiàng)P進(jìn)行延長得到Py;無特殊情況不考慮突變項(xiàng);利用經(jīng)遺傳算法優(yōu)化的神經(jīng)網(wǎng)絡(luò)預(yù)測隨機(jī)項(xiàng)得Ry;然后,將延長后的趨勢項(xiàng)Ty、周期項(xiàng)Py與神經(jīng)網(wǎng)絡(luò)預(yù)測所得隨機(jī)項(xiàng)Ry相加得預(yù)測結(jié)果。因降雨時(shí)間序列變化復(fù)雜,線性延長特征點(diǎn)法末端直線只能在短期滿足精度,DMAF模型適用于降雨時(shí)間序列的短期預(yù)測。
(3)關(guān)于利用DMAF 模型對降雨時(shí)間序列預(yù)測的步長問題:桓臺縣1979-2018年的480 組月降雨數(shù)據(jù),時(shí)間步長為月,若降雨數(shù)據(jù)能滿足一定的時(shí)間步長,研究提出的DMAF 模型還能對實(shí)時(shí)數(shù)據(jù)(30 min 或1 h 等步長)、周數(shù)據(jù)、年數(shù)據(jù)等進(jìn)行挖掘及預(yù)測。
(4)DMAF模型對流域的選擇沒有嚴(yán)格要求,因不考慮流域溫度、下墊面條件、海陸位置、流域河流結(jié)構(gòu)等因素,不同流域水文數(shù)據(jù)均能采用本模型進(jìn)行數(shù)據(jù)挖掘與預(yù)測。
(5)周期項(xiàng)由第一主周期小波系數(shù)確定,隨機(jī)項(xiàng)為原始數(shù)據(jù)減去趨勢項(xiàng)、周期項(xiàng)與突變項(xiàng)所得余項(xiàng),兩者存在負(fù)值是合理的,單位為mm。
以桓臺縣1979-2018年月降雨量數(shù)據(jù)為例,構(gòu)建了DMAF模型進(jìn)行數(shù)據(jù)挖掘與預(yù)測研究,主要結(jié)論如下。
(1)DMAF模型能在充分挖掘降雨時(shí)間序列的基礎(chǔ)上,對其進(jìn)行準(zhǔn)確預(yù)測,模型將時(shí)間序列分解為趨勢項(xiàng)、周期項(xiàng)、突變項(xiàng)與隨機(jī)項(xiàng),相較直接預(yù)測準(zhǔn)確度較高?;概_縣1979-2018年月降雨量數(shù)據(jù)挖掘與預(yù)測的結(jié)果表明,數(shù)據(jù)有微弱的上升趨勢,預(yù)測未來將呈微弱下降趨勢,其第一主周期是19(月),數(shù)據(jù)不存在明顯的突變情況。對桓臺縣2018年月降雨數(shù)據(jù)進(jìn)行預(yù)測,預(yù)測值與真實(shí)值較為接近,預(yù)測結(jié)果誤差較小,為16.79%。
(2)DMAF模型短期預(yù)測精度較高,不適宜長期預(yù)測。
(3)DMAF模型的推廣性較強(qiáng),能對實(shí)時(shí)數(shù)據(jù)、周數(shù)據(jù)、月數(shù)據(jù)、年數(shù)據(jù)等進(jìn)行挖掘與預(yù)測。
(4)DMAF模型的流域通用性較廣,對流域的選擇沒有嚴(yán)格要求,可適用于其他流域數(shù)據(jù)。 □