泰安市蒸發量變化趨勢分析與基于神經網絡的預測

于小鴿, 王世超, 李巖, 錢麗麗

(1.山東科技大學資源學院, 泰安 271000; 2.山東科技大學地球科學與工程學院, 青島 266500;3.泰安市水文中心, 泰安 271000)

蒸發量是水資源調查與評價、水利工程、農田灌溉、城市用水等方面必須要考慮到的一個很重要的水文要素,是水文領域一直在研究的一個重要課題[1]。蒸散發是水循環重要的一環,約66%的降水通過蒸散發的途徑返回大氣,蒸發量的變化分析及預測對于合理利用研究區的水資源,提高研究區水資源的利用率有著重要意義[2]。

關于蒸發量的分析和預測的研究在中外已經取得了很大的進展,如劉蕊蕊等[3]應用M-K突變檢驗(Mann-Kendall檢驗法)、小波分析和灰色關聯度法,對石羊河流域蒸發量變化特征及影響因素進行了分析,為石羊河流域水資源變化及農業灌溉提供一定的理論依據;柴小輝[4]應用M-K突變檢驗和均值生成函數預測模型法,對宕昌縣蒸發量觀測資料特征值的變化趨勢進行了分析研究,為宕昌縣水資源評價、干旱分析等提供了技術支持和參考;近些年來隨著深度學習火熱,更是有基于深度學習的蒸發量預測研究不斷涌出,如史俊才等[5]應用帶有注意力機制的卷積神經網絡模型對北京地區蒸發量進行了預測;董力銘等[6]建立了分類梯度提升算法與蝙蝠算法耦合建模預測中國西北部地區水面蒸發量等,然而基于深度學習的預測模型受氣象因子限制較大,現實中這些氣象因子往往難以獲取,且無法提前對未來的水文工作起到指導作用,因此時間序列預測模型仍是目前的主流模型之一。不過,由于水文時間序列數據的多樣性,目前在模型的泛化性上仍存在有一定的空白。

現結合NeuralProphet算法和Optuna算法,創建一個可以適應不同數據特征的蒸發量預測模型。這種模型不僅具有高度的精確性,而且對各種不同的數據都具有很強的魯棒性,可以為水資源管理者提供一種預知水文信息變化的新途徑。

1 研究區概括與研究方法

1.1 研究區概括

泰安市(116°20′E～117°59′E,35°38′N～36°28′N)位于山東省中部,泰山山腳以南的平原上,總面積約7 762 km2[7]。作為山東省與黃河流域的重要紐帶,泰安市承擔著水資源管理、生態環境保護以及地區協同發展等多重任務。泰安市有著典型的溫帶大陸性季風型氣候,雨熱同期,歷史年平均降水量為690.6 mm,其中6—9月為528.6 mm,占全年的76.5%,3—5月為97.7 mm,占全年的14.1%,具有春旱、夏澇、深秋旱、年降雨量高低交替明顯的特點[8]。泰安市各代表性水文觀測站年平均水面蒸發量一般為1 000～1 220 mm,空間上東部丘陵地區小于西部平原地區。

利用東周水庫、黃前水庫、大汶口和戴村壩4個水文觀測站(圖1)1985—2021年蒸發量逐月觀測資料進行分析,并按3—5月、6—8月、9—11月、12—次年2月依次分為春、夏、秋、冬4個季節。四水文觀測站建站年份早,集水面積大,檢測資料齊全,地理位置上自東向西分布,可以從空間上將其作為泰安市蒸發量特點的代表性水文觀測站。使用泰山觀測站2005—2022年蒸發量日值數據進行預測,并按時間整理成周值和月值數據。泰山觀測站位于泰安市泰山區水文局院內,數據由山東省水文局采集和處理,可以代表泰安市整體蒸發量特征。蒸發量數據在3—11月使用E601型蒸發器,其他時間使用20 cm口徑蒸發器進行測量,所有蒸發量觀測數據均由泰安市水文中心提供,具有較高的可靠性。在使用數據之前,數據經過了質量控制和處理,包括刪除異常數據,對缺失數據進行插值處理等。泰安市水文觀測站分布如圖2所示。

圖1 泰安市代表性水文觀測站Fig.1 Representative hydrological observation station in Tai’an City

圖2 泰安市水文站點分布Fig.2 Distribution of hydrological stations in Tai’an City

1.2 研究方法

1.2.1 Mann-Kendall檢驗法

氣候的突變現象指的是氣候從一種穩定的狀態跨越式的跳轉到另一種穩定的狀態,本文研究使用目前最廣泛使用的M-K突變檢驗對泰安市4個代表性水文觀測站蒸發量觀測數據進行分析[9-12],對于具有n個樣本量的時間序列x,構造一秩序列Sk,即

(1)

式(1)中:

(2)

式(2)中:xi、xj分別為時間序列x在第i、j時刻的值;j=1,2,…,n?？梢?秩序列Sk是第i時刻數值大于j時刻數值個數的累計值。

在時間序列隨機獨立的假定下,定義統計量為

(3)

式(3)中:UF1=0;E(Sk)、var(Sk)分別累計數Sk的均值和方差。按時間序列x逆序,再重復上述過程,同時使UBk=-UFk,UB1=0。若UF、UB曲線在置信區間內有交點,該交點對應的時間為可能的突變開始時間點。當出現歧義點時,結合滑動t檢驗[13-15]確定真正的突變點。兩種方法均使用α=0.05的顯著性檢驗。

1.2.2 R/S分析法

R/S分析法(rescaled range analysis)通常用來分析時間序列的分形特征和長期記憶過程,是一種基于長程相關思想的時間序列分析方法[16-18]。通過R/S分析法得出的Hurst指數的大小可以判斷當前氣候特征是具有持續性還是反持續性。其Hurst指數計算步驟如下。

(1)設有時間序列x1,x2,…,xn。將數據分為長度為r的g組互不重疊的子序列,對每一組子序列,進行以下運算。

(4)

式(4)中:xij為第i個子序列中第j個值。

計算離差yij,即

(5)

計算累計離差Zij,即

(6)

式(6)中:yik為第i個子序列中第k個值的離差。

計算極差Ri,即

Ri=max(Zij)-min(Zij)

(7)

式(7)中:max(Zij)和min(Zij)分別為累計離差Zij的最大值和最小值。

計算標準差Si,即

(8)

計算RS值RSi,即

(9)

(2)將求出來的各子序列的RS求平均并求出其標準差,即為Hurst指數。Hurst指數介于0～1,若大于0.5,則表明該時間序列具有長期相關的特征。

1.2.3 Prophet算法

Prophet算法是一種時間序列模型中的GAM模型(加法模型),將時間序列分解成了趨勢項g(t)、季節項s(t)、假期項h(t)以及誤差項∈t,公式[19-21]為

y(t)=g(t)+s(t)+h(t)+∈t

(10)

式(10)中:∈t為非線性趨勢殘差,也叫隨機波動項。

趨勢項g(t)組成為

g(t)=(k+A?)t+[m+A(-s?)]

(11)

式(11)中:k、m為初始增長率。式(11)包括一個趨勢基準和一個增長變化,?為增長率的變化量,A為指示函數,(k+A?)t組成增長變化,s為當前時間戳,[m+A(-s?)]組成趨勢基準。

季節項s(t)使用傅里葉級數來表示,即

(12)

式(12)中:T為趨勢周期,周期為年時T=365.5,周期為月時T=30;un和vn為在模型逐漸擬合數據的過程中可確定的參數;對于以年為周期的序列,N一般取值為10,對于以周為周期的序列,N一般取值為3。

在本文模型中,假期項h(t)不予考慮。Prophet算法會對以上幾項分別進行擬合并等權相加,輸出最終結果。

1.2.4 NeuralProphet-Optuna模型設計

NeuralProphet基于神經網絡對時序數據進行預測[22-24],增加一個AR-Net(自回歸網絡)到分解的分量上進行時序自相關。Optuna是一個自動超參優化算法,對于給定的搜索范圍和搜索空間,Optuna可以不斷地對其進行采樣并縮小范圍,通過不斷地試錯給出最優的超參數值,Optuna基于貝葉斯算法中的Parzen樹優化,通過Parzen樹來模擬高斯過程中的p(x|y),即

(13)

式(13)中:y*為觀察后找到的最佳值;ξ(x)為觀察不同觀察值時形成的密度;g(x)為剩余觀察值形成的密度。

對于NeuralProphet模型來說,正則項的調整改變的是參數的分布的形狀,往往一個參數的改變都會很大的限制季節分解項靈活程度,進而對預測結果產生很大的影響,為解決此問題,本文研究加入了Optuna算法對最優參數解進行調整,以此來使模型的預測效果更精準,模型結構如圖3所示。

圖3 模型構建流程Fig.3 Model building process

1.2.5 模型評價指標

使用平均絕對誤差(MAE)、均方根誤差(RMSE)和確定系數(R2)作為模型評價指標。RMSE和MAE的定義[25-27]為

(14)

(15)

(16)

RMSE和MAE用于衡量模型預測值與實際觀測值之間的平均差異。MAE計算每個預測值與對應實際值之間的絕對誤差,然后將它們求和并取平均值。MSE是預測值與真實值之差的平方的平均值,RMSE是MSE的平方根,表示預測誤差的平均值,與MSE具有相同的單位。

此外,R2用于衡量模型的擬合度。其定義為

(17)

所有指標都是在測試集上計算的,以評估模型的泛化性能。RMSE和MAE越低以及R2越高,模型的性能越好。

2 分析結果與模型驗證

2.1 基于1985—2021年4站點數據的蒸發量變化特征分析

在分析總結東周水庫、黃前水庫、大汶口和戴村壩4個站點1985—2021年蒸發量檢測數據的基礎上,繪制成四季蒸發量變化圖、月際蒸發量變化圖以及年際蒸發量變化圖,如圖4～圖6所示。

圖4 泰安市四季蒸發量變化特征Fig.4 Seasonal evaporation variation characteristics in Tai’an City

圖5 泰安市年內蒸發量變化特征Fig.5 Annual evaporation variation characteristics in Tai’an City

圖6 泰安市年總蒸發量變化特征Fig.6 Annual total evaporation variation characteristics in Tai’an City

泰安市各季蒸發量均呈減少趨勢。尤以雨季減少趨勢最為顯著,其中夏季減少趨勢較明顯,春季次之;干季尤以冬季為主,蒸發量幾乎無顯著變化。泰安市月際蒸發量變化屬于單峰型,5、6月份達到最大值,分別為107.9、137.5、134.6、123.8 mm,占全年比例的12.5%、14.6%、14.8%和14.2%;1月份達到全年數值最小,分別為17.7、20.2、15.9、23.4 mm,占全年比例數為2.3%、2.1%、1.7%和2.7%[28]。

泰安市地區整體蒸發量呈十分顯著的減少趨勢。但由于4個水文觀測站不同的地理位置差異、地形地貌差異、植被覆蓋和人為活動影響,從空間分配上分析,無論是季節還是年總蒸發量,黃前水庫和大汶口呈顯著的下降趨勢,而東周水庫和戴村壩呈微弱的上升趨勢。其中黃前水庫的整體蒸發量要比東周水庫、大汶口和戴村壩明顯少得多。

2.2 基于M-K突變檢驗和滑動t檢驗的蒸發量突變特征分析

根據M-K突變檢驗分析的方法得到了4個水文站四季蒸發量和年總蒸發量的MK突變檢驗圖,如圖7～圖11所示。由圖7分析可得,對于黃前水庫,春季在1988年正序列UF和逆序列UB出現交點,并于2001年突破-1.96(α=0.05)的臨界線,意味著春季在1988年出現了明顯的突變減少;夏季UF曲線和UB曲線在1995年出現交點,并于2009年突破-1.96(α=0.05)的臨界線,這表明在1995年夏季發生了明顯的突變減少;秋季UF曲線和UB曲線有多個交點,突變點分別為:1988年,1990年,1997年,并于2009年突破-1.96(α=0.05)的臨界線,結合滑動t檢驗可知真正的突變點為1997年,這說明秋季在1997年發生顯著突變減少;冬季與夏季同理,沒有發生突變。

圖7 黃前水庫MK突變檢驗結果Fig.7 MK abrupt change test results for Huangqian Reservoir

由圖8分析可得,對于東周水庫,春季分別在2008年和2011年UF曲線和UB曲線交于一點,并于2018年突破1.96(α=0.05)的臨界線,結合圖3可明顯確認2008年已經開始突變,這說明春季從2008年發生顯著突變增多;夏季的UF曲線和UB曲線有多個交點,并于1999年突破-1.96(α=0.05)的臨界線,但交點均位于端點附近,結合滑動t檢驗可見夏季并沒有發生突變;秋季沒有發生突變;冬季則是在1990年發生顯著突變減少。

圖8 東周水庫MK突變檢驗結果Fig.8 MK abrupt change test results for Dongzhou Reservoir

由圖9分析可得,對于大汶口,春季從1988年發生顯著突變減少;夏季UF曲線和UB曲線有多個交點,并于2000年突破-1.96(α=0.05)的臨界線,結合滑動t檢驗確定1989年為真實的突變點并發生顯著突變減少;秋季沒有突破1.96(α=0.05)的臨界線;冬季與秋季情況相同均未發生突變。

圖9 大汶口MK突變檢驗結果Fig.9 MK abrupt change test results for Dawenkou

由圖10分析可得,對于戴村壩,春季與夏季都有正逆序列交點并突破臨界線,分別從1988年和1993年發生顯著突變減少;秋季UF曲線和UB曲線于2012出現交點,但是沒有突破1.96(α=0.05)的臨界線,秋季沒有發生突變;冬季與秋季情況相同。

圖10 戴村壩MK突變檢驗結果Fig.10 MK abrupt change test results for Daicunba

由圖11分析可得,對于年總蒸發量,黃前水庫的正序列UF和逆序列UB在1989年和1991年出現了兩個交點,并于2000年突破-1.96(α=0.05)的臨界線,結合滑動t檢驗確認1989年為真實的突變點;東周水庫年蒸發量并未突破-1.96(α=0.05)的臨界線,沒有發生突變;大汶口在1987年發生顯著突變減少;戴村壩分別在1986年和1989年UF曲線和UB曲線出現交點,并于1997年突破-1.96(α=0.05)的臨界線,由于1986年位于端點處不取,認為戴村壩從1989年發生顯著突變減少。并給出前面分析所用到滑動t檢驗圖,如圖12所示。

圖11 年總蒸發量MK突變檢驗結果Fig.11 MK abrupt change test results for annual total evaporation

圖12 滑動t檢驗Fig.12 Sliding t-test

2.3 基于R/S分析法的蒸發量變化趨勢預測

使用R/S分析法預測未來蒸發量變化趨勢。R/S分析法是一種在水文時間序列中被廣泛采用的用于分析持續性方法。它的原理基于統計物理學,已經在多個研究中證明了其可靠性。這種方法不僅可以提供對未來蒸發量的準確預測,而且其預測結果對于實際的水資源管理和規劃具有高度的實用性。泰安市4個觀測站hurst指數如表1所示,結果表明,4個觀測站均出現了明顯的hurst現象(H>0.5),表明年和四季蒸發量變化將會繼續維持之前的整體變化趨勢。由前述可知,泰安市蒸發量在研究時段內整體呈下降趨勢,可以預測,在未來一段時間內泰安市蒸發量將繼續下降,這將有利于農業生產活動,對于實際的水資源管理和規劃具有一定的實用性。

表1 4觀測站hurst指數Table 1 Hurst index of four evaporation stations

2.4 模型驗證與分析

Optuna算法對NeuralProphet模型的5個重要參數進行最優搜索,包括趨勢變化點范圍(changepoints_range)、年周期階數(year_season_order)、學習率(learning_rate)、趨勢變化點數量(n_changepoints)和滯后值(n_lags),最優參數搜索過程如表2所示。

表2 最優參數搜索過程Table 2 Optimal parameter search process

Optuna算法對表2中5個重要參數都給出了調整,這些參數的改變調整了模型的趨勢性組件和季節性組件,如圖13所示。通過對趨勢變化點范圍和數量的改變控制了時間序列中趨勢變化的可能位置,從而影響模型的趨勢組件;年周期階數的設置決定了年季節性組件的復雜度,使模型能夠靈活地適應年度循環變化;而不同大小的滯后值定義了自回歸組件中模型考慮的歷史觀測值的數量,從而影響模型的自回歸組件。

圖13 模型組件變化Fig.13 Model component changes

利用泰山站2005—2022年蒸發量逐日觀測數據,以7∶3的比例劃分訓練集和測試集,分別訓練Prophet模型、NeuralProphet模型與NeuralProphet-Optuna模型并預測,取部分預測結果進行可視化,如圖14所示。從圖14中可以看出,3種模型的預測精度差異明顯,其中NeuralProphet-Optuna模型的預測結果與真實值最為接近,其次為NeuralProphet算法。3種模型在時間序列的長短期趨勢性和季節性上均表現出了不錯的學習效果,但由于泰安市日值蒸發量復雜的水文特征,加入了Optuna算法的模型預測效果表現出了明顯的提升。特別是在蒸發量偏高且變化復雜的夏季,蒸發量數據經常發生短期的突變和波動,憑借著NeuralProphet算法中的自回歸組件,模型仍能取得較好的預測效果。為更明確比較模型的擬合效果,圖15給出了3種模型的預測分布和最優擬合直線。

圖14 模型預測結果對比Fig.14 Comparison of model prediction results

圖15 模型擬合效果對比Fig.15 Model fitting effect comparison

可見,NeuralProphet-Optuna模型的擬合優度高于其他兩個模型。尤其是在蒸發量值偏高的時候,NeuralProphet-Optuna模型給出了明顯優于其他兩個模型的預測結果,這與圖14所給出的信息是一致的,說明Optuna算法和自回歸組件的加入使模型更能學習到蒸發量數據中復雜的季節性和趨勢性特征。

為研究不同時間增量下模型的預測精度,按時間將日值數據整理成周值、月值數據分別對模型進行訓練并預測,并與經典的時間序列預測模型ARIMA、LSTM、Holt-Winters模型作比較,分別給出各個模型的評價指標,如表3所示。由于ARIMA模型主要關注序列的平穩性,對突變點的捕捉能力較弱,在預測日值的蒸發量數據時模型的預測效果很差,故不考慮。

表3 模型評價指標對比Table 3 Comparison of model evaluation metrics

可見,NeuralProphet-Optuna模型在3種不同的時間增量下都表現出了最優的預測精度,3種時間增量下NeuralProphet-Optuna模型的確定系數R2平均比改進前的NeuralProphet模型和Prophet模型提升了3.8%和7.6%,相比傳統的LSTM模型、Holt-Winters模型和ARIMA模型平均提升了12.8%、7.6%和5%。其中以周為時間增量時預測效果最佳,其次為月,然后是日。相比較日值和月值數據,周值數據下模型的確定系數R2分別提升了6.2%和4.6%。這是因為日值蒸發量數據中存在大量噪聲和短期波動,模型學習難度高,將數據整理成周值數據可以減小這些噪聲和短期波動的影響,提升了模型的預測精度。而將數據整理成月值數據時,樣本數量大幅降低,模型魯棒性降低,故不如周值數據下的預測。

3 結論

依據泰安市4個代表性水文觀測站的蒸發量觀測數據,對泰安市蒸發量變化過程進行了分析及預測,具體結論如下。

(1)通過M-K突變檢驗結合滑動t檢驗法的方法檢測泰安市蒸發量的突變特征,得到泰安市4個觀測站的四季蒸發量情況以及總蒸發量情況。從時間上來講,給出了4個蒸發站點年和四季蒸發量的突變特征及其準確的突變年份;從空間上來講,黃前水庫和大汶口蒸發總量將維持穩定下降的趨勢,其他兩個觀測站的變化幅度不大但也會維持目前的發展趨勢。

(2)基于NeuralProphet算法和Optuna算法提出了一種蒸發量預測模型,通過與其他經典的預測模型作比較,可以看出NeuralProphet-Optuna模型的預測結果更精準,穩定性更高。然而,NeuralProphet-Optuna模型在預測日值蒸發量時效果不如周值和月值數據,盡管依賴著其趨勢變化點識別能力和自回歸組件能夠很好地捕捉一些序列中的突變點和短期趨勢,但在學習一些極端的蒸發量日值時效果仍有限,未來可以考慮使用多種水文參數訓練深度學習類預測模型來解決這一問題。

(3)總體來說,本文研究采用了前沿的方法對泰安市水面蒸發量進行了全面分析和預測。所提出的模型不僅具有較高的精度和較好的魯棒性,而且為未來的水文研究和水資源管理者提供了一種有力的技術支持,具有一定的實用性和推廣價值。