基于時間序列分解的降雨數(shù)據(jù)挖掘與預(yù)測

趙然杭，甘 甜，逄曉騰，王興菊，茍偉娜，齊 真

（1.山東大學(xué)土建與水利學(xué)院，濟(jì)南250100；2.青島市供水事業(yè)發(fā)展中心，山東青島266071）

0 引 言

降雨是水文過程的重要環(huán)節(jié)，能直接反映區(qū)域氣候變化及水文過程，降雨量變化直接關(guān)乎區(qū)域水資源開發(fā)利用乃至影響社會經(jīng)濟(jì)發(fā)展。近年來，由于氣候變化與人類活動的影響，生態(tài)問題嚴(yán)峻、洪旱災(zāi)害頻發(fā)，充分挖掘降雨時間序列所含信息能更好揭示其自然規(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ù)據(jù)挖掘進(jìn)行了相關(guān)研究：李宜軒運(yùn)用線性擬合的方法，綜合分析了桓臺縣1981-2018年的降水趨勢變化情況，結(jié)論與實際情況相符［2］；衡彤等采用小波變換方法，對新安江流域黃山地區(qū)降水量時間序列的多時間尺度變化及突變特征進(jìn)行了探討，研究分析了主汛期降水與年降水的主周期，兩者較為接近，結(jié)果合理［3］；李平蘭等分析了1970-2015年會東縣的降水變化特征，采用了累積距平法與Mann-Kendall 法的突變點(diǎn)綜合判定方法，彌補(bǔ)了單個方法檢驗?zāi)芰Φ牟蛔悖?］；Moravej等基于線性時間序列模型，利用ADF 函數(shù)、Mann-Kendall 檢驗，對伊朗West Azerbaijan 省的降雨時間序列進(jìn)行了趨勢分析，模型結(jié)論準(zhǔn)確且適用性較強(qiáng)［5］。

國內(nèi)外學(xué)者對降雨時間序列數(shù)據(jù)預(yù)測進(jìn)行了相關(guān)研究：王曉鵬等應(yīng)用基于Box-Jenkins方法的時間序列分析技術(shù)，針對青南高原的四個典型地區(qū)1961-2005年降水量序列，建立了自回歸滑動平均模型對未來降水量進(jìn)行了預(yù)測，模型因能反復(fù)識別修改，適用于復(fù)雜情況的數(shù)據(jù)擬合與預(yù)測［6］；賈海峰等提出了灰色-時序組合預(yù)測模型，能對年降雨量進(jìn)行中長期預(yù)測，并以北京市大興縣黃村氣象站35年降雨量資料為例進(jìn)行了驗證，模型能滿足精度要求，且對其他既具有擺動又有趨勢數(shù)據(jù)的模擬和預(yù)測同樣適用［7］；高瑞華利用煙臺地區(qū)1981-2010年降水資料，建立了干旱的灰色預(yù)測GM（1，1）模型，對2010-2030年的干旱災(zāi)害進(jìn)行了預(yù)測，模型通過了殘差與后驗差檢驗分析，在短期預(yù)測時精度能滿足要求，但不適用于長期預(yù)測［8］；Kolachian 等采用分位數(shù)映射與平均貝葉斯模型，通過提高觀測值與集合預(yù)報平均值之間的相關(guān)系數(shù)，對歐洲中心2020年各月份不同降水模式的若干氣象站的觀測降水序列進(jìn)行了預(yù)報，結(jié)果較為準(zhǔn)確［9］。

以上研究分別對降雨時間序列進(jìn)行了數(shù)據(jù)挖掘或預(yù)測，沒有將數(shù)據(jù)挖掘與預(yù)測建立緊密聯(lián)系，而數(shù)據(jù)挖掘能為數(shù)據(jù)預(yù)測提供信息支撐，數(shù)據(jù)預(yù)測結(jié)果也能驗證挖掘的準(zhǔn)確性，兩者相輔相成。由于降雨時間序列受趨勢、周期與突變等因素影響，復(fù)雜多變，需要在充分挖掘數(shù)據(jù)信息的基礎(chǔ)上，進(jìn)行降雨時間序列的預(yù)測研究。因此開展基于時間序列分解的降雨數(shù)據(jù)挖掘與預(yù)測研究，為雨（洪）水資源管理工作提供決策技術(shù)支持。

1 研究理論與方法

1.1 時間序列法簡介

降雨時間序列因受多種因素影響，存在著趨勢性、周期性、突變性等多種變化特征，將時間序列分解為趨勢項、周期項、突變項與隨機(jī)項并對各項進(jìn)行分析與研究能更好揭示降雨量數(shù)據(jù)的變化規(guī)律。設(shè)X為時間序列，T為趨勢項、P為周期項、B為突變項、R為隨機(jī)項，則X可以表示為：

常見的趨勢項分析法包括：累積距平［10］、Mann-Kendall 趨勢分析［11］、Hurst 指數(shù)［12］、特征點(diǎn)［13］，周期項分析法包括：小波分析［14］，突變項分析法包括：Mann-Kendall 突變檢驗［15］、Pettitt 突變檢驗［16］，隨機(jī)項檢驗法包括：自相關(guān)圖和單位根［17，18］。

1.2 基于時間序列分解的降雨數(shù)據(jù)挖掘-預(yù)測模型

提出了基于時間序列分解的降雨數(shù)據(jù)挖掘-預(yù)測模型（Data Mining And Forecasting model，DMAF 模型），不僅能挖掘出降雨時間序列的內(nèi)在信息，還能對其進(jìn)行準(zhǔn)確預(yù)測。

（1）數(shù)據(jù)挖掘方面：①數(shù)據(jù)分解：對于降雨時間序列X，采用特征點(diǎn)法分解其趨勢項T；小波分析法分解其周期項P；Mann-Kendall法、Pettitt法進(jìn)行突變性分析，無特殊情況突變項B不予考慮；計算余項X-T-P-B得隨機(jī)項R。②數(shù)據(jù)挖掘：分別采用累積據(jù)平法、Mann-Kendall趨勢分析法、Hurst指數(shù)法、特征點(diǎn)法進(jìn)行趨勢分析，采用小波分析法進(jìn)行周期分析，采用Mann-Kendall 法、Pettitt 法進(jìn)行突變性分析，采用自相關(guān)圖法、單位根法對隨機(jī)項進(jìn)行檢驗。

（2）數(shù)據(jù)預(yù)測方面：利用線性延長方法將趨勢項T基于特征點(diǎn)法的末端直線進(jìn)行延長得到Ty，利用Matlab 小波工具箱Extend 函數(shù)將周期項P進(jìn)行延長得到Py，無特殊情況下不考慮突變項，采用NAR 神經(jīng)網(wǎng)絡(luò)法對隨機(jī)項進(jìn)行預(yù)測。然后，將延長后的趨勢項Ty、周期項Py與預(yù)測所得隨機(jī)項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

2 基于時間序列分解的降雨數(shù)據(jù)挖掘-預(yù)測應(yīng)用實例

2.1 研究區(qū)域概況與數(shù)據(jù)來源

桓臺縣是淄博市下轄縣，位處魯中山區(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ù)來源于桓臺縣水文站。

2.2 時間序列分解

利用matlab 進(jìn)行編程，以桓臺縣1979-2018年降雨數(shù)據(jù)為例，時間步長為月，將時間序列分解為趨勢項、周期項、突變項與隨機(jī)項進(jìn)行數(shù)據(jù)挖掘。

（1）趨勢項：分別采用累積距平、Mann-Kendall、Hurst 指數(shù)等方法進(jìn)行趨勢性分析，特征點(diǎn)法進(jìn)行趨勢項提取。累積距平圖如圖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，表明該時間序列具有長期相關(guān)性，但未來總體趨勢與過去相反。

表1 Mann-Kendall 趨勢檢驗結(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ù)測下一時段將呈現(xiàn)微弱下降趨勢。

以上為定性分析，現(xiàn)根據(jù)特征點(diǎn)法定量分解趨勢項。取閾值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)計算插值法誤差為79.50，用回歸法連接后效果如圖8所示（圖8中紅色數(shù)字為各段回歸系數(shù)）。

圖8 趨勢項分解圖Fig.8 Trend term decomposition

插值法更為直觀，便于描述序列變化情況。回歸法誤差為53.53，明顯小于插值法，因此，選用回歸法的結(jié)果為時間序列的趨勢項。

（2）周期項：使用Matlab 中的小波工具箱得小波方差圖，如圖9所示。由圖9可以看出，桓臺縣1979-2018年月降雨量的第一主周期為19（月）。以第一主周期時間尺度繪制小波系數(shù)圖，如圖10所示。選取第一主周期小波系數(shù)為周期項。

圖9 小波方差圖Fig.9 Wavelet variance

圖10 第一主周期小波系數(shù)圖Fig.10 First principal period wavelet coefficients

（3）突變項：分別采用Mann-Kendall 法、Pettitt 法進(jìn)行突變性檢驗，結(jié)果如圖11、表2所示：

由圖11 可得，1979年-2018年的月降雨量突變情況發(fā)生少。表2中，突變點(diǎn)檢驗結(jié)果不顯著，不能認(rèn)為產(chǎn)生突變。為避免單個檢驗方法決策失誤，采用文獻(xiàn)［4］中突變點(diǎn)判別方法，認(rèn)為在兩種判別方法中均為突變點(diǎn)的月份才發(fā)生顯著突變。綜合分析，認(rèn)為1979-2018年月降雨不發(fā)生顯著突變情況，不進(jìn)行突變項分解。

表2 Pettitt法檢驗結(jié)果Tab.2 Results of Pettitt method

圖11 Mann-Kendall突變檢驗圖Fig.11 Mann-Kendall mutation test

（4）隨機(jī)項：將時間序列減去趨勢項、周期項與突變項可得隨機(jī)項。若隨機(jī)項的自相關(guān)性較小且具有平穩(wěn)性，說明趨勢項、周期項、突變項信息被充分提取。利用Eviews 軟件對隨機(jī)項進(jìn)行自相關(guān)（AC）、偏自相關(guān)（PAC）檢驗，結(jié)果如圖12所示。

圖12 自相關(guān)圖Fig.12 Autocorrelation diagram

由圖12得，AC值與PAC值始終圍繞零軸小范圍波動，進(jìn)一步進(jìn)行單位根檢驗，得P值 為0、AIC為10.418 02、SC為10.54212、LogL為-2 423.818。由自相關(guān)與單位根檢驗結(jié)果可知，隨機(jī)項不存在單位根，為平穩(wěn)時間序列。趨勢項、周期項、突變項信息已被較為充分提取。

2.3 時間序列預(yù)測

2.3.1 隨機(jī)項預(yù)測

為便于驗證預(yù)測結(jié)果，以1979-2014年432 組月降雨量數(shù)據(jù)為率定數(shù)據(jù)，2015-2016年月降雨量數(shù)據(jù)為驗證數(shù)據(jù)，對2017-2018年月降雨數(shù)據(jù)進(jìn)行預(yù)測。

隨機(jī)項預(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，中間層個數(shù)取8，訓(xù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ù)測：為了驗證DMAF模型預(yù)測結(jié)果的合理性，采用NARX神經(jīng)網(wǎng)絡(luò)隨機(jī)項預(yù)測法，結(jié)果見圖16。

2.3.2 時間序列的總體預(yù)測

將圖15 和圖16 的預(yù)測結(jié)果，分別與所對應(yī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ī)項預(yù)測結(jié)果圖Fig.15 Random term prediction result chart based on NAR neural network

圖16 基于NARX神經(jīng)網(wǎng)絡(luò)隨機(jī)項預(yù)測結(jié)果圖Fig.16 Random term prediction result chart based on NARX neural network

圖17 預(yù)測結(jié)果對比圖Fig.17 Comparison chart of prediction results

3 結(jié)果分析與討論

3.1 結(jié)果分析

（1）時間序列挖掘：由圖3和表1可知，桓臺縣1979-2018年月降雨量數(shù)據(jù)有微弱的上升趨勢，此結(jié)論與文獻(xiàn)［2］相符。由圖3 可預(yù)測未來將呈微弱下降趨勢，由圖9 可知，其第一主周期是19（月），由圖11與表2可知，數(shù)據(jù)不存在明顯的突變情況。

（2）時間序列預(yù)測精度：由圖17 可知，與1979-2018年月降雨量的實測數(shù)據(jù)相比，基于NAR 和NARX 神經(jīng)網(wǎng)絡(luò)的預(yù)測誤差分別為16.79%和20.05%。利用NAR 神經(jīng)網(wǎng)絡(luò)直接預(yù)測法進(jìn)行時間序列的總體預(yù)測（結(jié)果見圖17），預(yù)測結(jié)果誤差較大（為41.62%）。因此，上述3 種預(yù)測方法，基于NAR 神經(jīng)網(wǎng)絡(luò)隨機(jī)項預(yù)測方法獲得月降雨數(shù)據(jù)的精度最高。

3.2 討 論

（1）關(guān)于降雨時間序列的數(shù)據(jù)挖掘與預(yù)測的關(guān)系：降雨時間序列受趨勢、周期與突變等因素影響，復(fù)雜多變，數(shù)據(jù)挖掘能為數(shù)據(jù)預(yù)測提供信息支撐，數(shù)據(jù)預(yù)測結(jié)果也能驗證挖掘信息的準(zhǔn)確性。而現(xiàn)有水文數(shù)據(jù)挖掘與預(yù)測方面的文獻(xiàn)［2-5］，未見有基于水文數(shù)據(jù)挖掘的預(yù)測研究，即沒有建立數(shù)據(jù)挖掘與預(yù)測的緊密聯(lián)系，而是直接預(yù)測。創(chuàng)新性地提出了DMAF 模型，將降雨時間序列分解為趨勢項、周期項、突變項及隨機(jī)項，并對各項進(jìn)行數(shù)據(jù)挖掘，同時建立NAR 神經(jīng)網(wǎng)絡(luò)模型對隨機(jī)項進(jìn)行預(yù)測，再將預(yù)測結(jié)果與數(shù)據(jù)挖掘后的趨勢項、周期項、突變項相加得降雨時間序列預(yù)測結(jié)果，相較于直接預(yù)測法，結(jié)果更為準(zhǔn)確。因此DMAF 模型能在充分挖掘降雨時間序列的基礎(chǔ)上，對其進(jìn)行準(zhǔn)確預(yù)測。

（2）關(guān)于利用DMAF 模型對降雨時間序列長短期預(yù)測問題：對研究區(qū)域的降雨時間序列進(jìn)行分解后，利用線性延長方法將趨勢項T基于特征點(diǎn)法的末端直線進(jìn)行延長得到Ty；利用Matlab 小波工具箱Extend 函數(shù)將周期項P進(jìn)行延長得到Py；無特殊情況不考慮突變項；利用經(jīng)遺傳算法優(yōu)化的神經(jīng)網(wǎng)絡(luò)預(yù)測隨機(jī)項得Ry；然后，將延長后的趨勢項Ty、周期項Py與神經(jīng)網(wǎng)絡(luò)預(yù)測所得隨機(jī)項Ry相加得預(yù)測結(jié)果。因降雨時間序列變化復(fù)雜，線性延長特征點(diǎn)法末端直線只能在短期滿足精度，DMAF模型適用于降雨時間序列的短期預(yù)測。

（3）關(guān)于利用DMAF 模型對降雨時間序列預(yù)測的步長問題：桓臺縣1979-2018年的480 組月降雨數(shù)據(jù)，時間步長為月，若降雨數(shù)據(jù)能滿足一定的時間步長，研究提出的DMAF 模型還能對實時數(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）周期項由第一主周期小波系數(shù)確定，隨機(jī)項為原始數(shù)據(jù)減去趨勢項、周期項與突變項所得余項，兩者存在負(fù)值是合理的，單位為mm。

4 結(jié) 論

以桓臺縣1979-2018年月降雨量數(shù)據(jù)為例，構(gòu)建了DMAF模型進(jìn)行數(shù)據(jù)挖掘與預(yù)測研究，主要結(jié)論如下。

（1）DMAF模型能在充分挖掘降雨時間序列的基礎(chǔ)上，對其進(jìn)行準(zhǔn)確預(yù)測，模型將時間序列分解為趨勢項、周期項、突變項與隨機(jī)項，相較直接預(yù)測準(zhǔn)確度較高。桓臺縣1979-2018年月降雨量數(shù)據(jù)挖掘與預(yù)測的結(jié)果表明，數(shù)據(jù)有微弱的上升趨勢，預(yù)測未來將呈微弱下降趨勢，其第一主周期是19（月），數(shù)據(jù)不存在明顯的突變情況。對桓臺縣2018年月降雨數(shù)據(jù)進(jìn)行預(yù)測，預(yù)測值與真實值較為接近，預(yù)測結(jié)果誤差較小，為16.79%。

（2）DMAF模型短期預(yù)測精度較高，不適宜長期預(yù)測。

（3）DMAF模型的推廣性較強(qiáng)，能對實時數(shù)據(jù)、周數(shù)據(jù)、月數(shù)據(jù)、年數(shù)據(jù)等進(jìn)行挖掘與預(yù)測。

（4）DMAF模型的流域通用性較廣，對流域的選擇沒有嚴(yán)格要求，可適用于其他流域數(shù)據(jù)。 □