基于時間序列的海流兔河流域降水量模擬預測分析

周生輝, 劉廷璽,2, 段利民,2, 張文瑞, 冀 如, 孫 龍

(1.內蒙古農業大學 水利與土木建筑工程學院, 呼和浩特 010018; 2.內蒙古自治區水資源保護與利用重點實驗室, 呼和浩特 010018)

大氣降水是水文循環的重要組成部分，是流域水資源的主要來源，近年來，隨著全球氣候變暖日趨明顯、極端氣候事件頻發，探究區域的降水循環規律成為了國內外研究的熱點[1-8]。降水量作為水文循環的重要驅動因子，直觀地反映了區域水資源的豐富程度，對區域水資源的量化及生態政策的制定至關重要，因此研究區域降水量的變化規律和未來趨勢具有現實意義。近年來，黃生志等[9-10]利用降水量與徑流的演變規律對渭河流域和哥倫比亞河流域未來干旱趨勢進行了預測分析，得到了具有數理統計意義的流域干濕特征規律。吳昊等[11]研究了陜西省14個氣象站的年、季降水量變化趨勢，結果發現渭河流域和黃土高原北部地區年降水量呈明顯下降趨勢。袁定波等[12]在泰森多邊形雨量法的基礎上研究了地理空間要素對降雨空間分布的影響，得到了雅礱江流域不同時間尺度的降雨量變化趨勢。韓知明等[13]利用小波分析對呼倫湖流域的降水序列進行了變化特征分析，揭示了該流域的降水變化主周期并對未來降水量的枯豐趨勢進行了預測。總之，眾多學者利用各種方法對大尺度區域降水量的規律進行了分析預測，而對于降水數據較少且需要重點關注的中小流域，相關研究卻并不常見。因此本文以毛烏素沙地中人類活動較為劇烈的海流兔河流域為研究對象，通過氣象站資料進行降水變化的周期分析和趨勢預測，為該區域以流域為單元的煤水協調開發提供氣候指導。

毛烏素沙地是中國四大沙地之一，經過多年的風沙治理及水土保持工作，其生態環境有了大幅改善，但近年來隨著域內農田開墾和煤炭資源開發，水資源平衡被重新調配，生態環境問題仍舊十分突出[14-19]。作為毛烏素沙地環境敏感區和未來蒙陜能源開發中心區的海流兔河流域受到了多方學者的關注，但基于該流域的降水量變化研究尚不多見[20-22]。本文收集流域周邊氣象站點的降水量資料，采用Moral小波分析、集合經驗態分解(EMD)、holt-winters模型對數據進行時間序列分析，最后利用泰森多邊形法加權以流域為整體單元，系統研究其多年降水的周期特征、變化規律以及未來降水趨勢的模擬預測，豐富海流兔河流域的降水研究成果。

1 數據與方法

1.1 研究區概況及數據來源

海流兔河流域位于陜蒙交界地區，地跨陜西省榆林市榆陽區、內蒙古自治區鄂爾多斯市烏審旗兩個縣市，處于毛烏素沙地邊緣與黃土高原接壤區。流域北起烏審旗烏蘭烏都，南至榆陽區王圪堵水庫，西起烏審旗木胡灘，東至榆陽區大海則，流域面積約2 600 km2，屬于典型的干旱、半干旱沙地和灘地相間分布的草原氣候環境，為溫帶大陸性季風氣候，大氣降水是域內主要的水資源補給項。區內多年降水量為370 mm，多年平均蒸發量為2 000 mm，降水多集中在7—9月份，風向以西北風為主。

本文所用的降水氣象資料來源于中國氣象數據網(http:∥data.cma.cn)，根據數據的完整性及可靠度，選取了海流兔河流域周邊最近的鄂托克旗站、榆林站、橫山站為研究對象，截取1959—2019年的逐月降水量同期數據，然后通過中國氣象數據網1961—2019年的0.5×0.5降水格點月數據集對3個站點的數據進行了對比檢驗。

1.2 研究方法

1.2.1 泰森多邊形法 各氣象站測得的降水量值只能代表點雨量，流域上的降水量值估算需要大量的氣象站通過面積加權進行計算。泰森多邊形法是將流域的所有雨量站之間連接直線，盡可能組成銳角三角形，對每個三角形求重心，利用得到的重心和三角形的邊垂直平分線將流域劃分成若干子區域，用每個子區域的雨量站代表降水量，加權面積后得到流域的平均降水量值[23]。具體計算過程如下：

式中:Ai為第i個子區域的面積;αi為第i個雨量站的面積權重;Pi為第i個子雨量站的降水量;n為雨量站數目。

1.2.2 Morlet小波分析 小波分析在時域和頻域兩個方面都具有良好的局部化功能，能揭示時間序列的多尺度變化特征，識別不同時間尺度的主要變化周期[24]。小波分析不僅能夠反映降水序列的局部變化特征，而且小波變換的結果可以顯示出氣候序列變化的尺度，以及顯示出序列變化的時間位置，能清晰地揭示出隱藏在時間序列中的多種變化周期，反映系統在不同時間尺度中的變化趨勢以估計未來的系統趨勢[25-26]。Morlet小波定義如下：

小波方差的計算公式為：

var(a)=∑[Wf(a,b)]2

1.2.3 EMD經驗模態分解 EMD經驗模態分解可對一個時間信號將其不同尺度即頻率的波動或趨勢逐級分解開來，產生一系列具有不同特征尺度的數據序列，稱為本征模函數(Intrinsic Mode Function, IMF)，它是目前處理非平穩、非線性信號，特別是分析時間序列趨勢的最好方法[27-28]，IMF分量可以反映出降水量序列中不同特征尺度的振蕩。其主要的計算原理如下：

(1) 提取原始時間序列X(t)中的極大值與極小值點，使用三次樣條插值函數對序列的上包絡線序列值Xmax(t)及下包絡線序列值Xmin(t)進行擬合，

(2) 計算上、下包絡線的均值m(t)：

(3) 使用原始數據序列X(t)減去(1) 中包絡線的均值得到一個新的序列I1(t)：

I1(t)=X(t)-m(t)

(4) 若所得新序列I1(t)滿足兩個條件：極值點與過零點數目相等或相差一個，包絡線的均值為0，則I1(t)為固有模態函數(IMF)。

(5) 用原始序列減去第一個固有模態函數I1(t)，得到剩余序列r1(t)：r1(t)=X(t)-I1(t)，將r1(t)作為新的原始序列，按照以上步驟，依次提取I2(t)，I3(t)，I4(t)，…，In(t)，直到rn(t)為一個單調序列。把分解后的各分量疊加，就得到原序列X(t)：

EMD方法在篩分過程中體現出了強大的直接性及自適應性，可根據實際序列的數理特征自主形成特殊的基函數，每一個分離出來的 IMF都具有一定的物理意義，都是對原始序列物理信息的真實反映[29]。

1.2.4 Holt-Winters模型 對于未來發生的事情，最新觀察值較早期觀察值包含更多的信息，因而在預測時，最新觀測值較早期觀察值具有更大的權重，Holt-Winters模型具有周期性特征，可對時間序列進行短期模擬預測[30]。Holt-Winters模型本質上是一種高級指數平滑形式模型，具有周期調整與長期趨勢調整特性，能對兼有長期趨勢和季節模式的數據進行預測[31-33]。本文利用的Holt-Winters乘性模型的基本原理方程如下：

bt=γ(St-St-1)+(1-γ)bt-1

Holt-Winters模型預測值計算為：

yt+s=(St+btk)It+k-L。

式中:St為平滑值即水平分量；α為水平權重;bt為長期趨勢值；γ為趨勢權重;It為季節分量；β為季節權重;L為季節長度(每年的月數或季數);t為當前時間;Xt為實際觀測值;yt+s為預測值；k為預測超前期數;其中的γ,α,β范圍為0～1。

Holt-Winters乘性模型中的γ，α，β取值依賴于已知時間序列的性質，通常為0.1～0.3的數值并產生一個依賴于大量過去觀測資料的預測[34]。

2 結果與分析

2.1 降水量年內特征分析

海流兔河流域北側的鄂托克旗站距流域約65 km，1959—2019年平均降水量為274 mm；流域南側的橫山站距流域約12 km，1959—2019年平均降水量為383 mm；流域東側的榆林站距流域約38 km，1959—2019年平均降水量為420 mm，流域從北向南，從西向東降水量逐漸增大。三站61 a線性降水量趨勢均呈上升狀態，氣象站記錄的降水量越大，則線性上升趨勢越明顯(圖1)。

圖1 氣象站1959-2019年逐月降水量變化

通過分析三站61 a的月平均降水量可以發現(圖2)，海流兔河流域多年平均月降水量呈明顯的單峰型，年內降水主要集中在7—9月份，分別占年均降水量的64.16%(鄂托克旗站)，61.19%(橫山站)和64.21%(榆林站)，基本占全年一半以上，具體表現為春冬枯，夏秋豐的降水時間分布格局。三站年內降水分布基本類似，降水峰值都集中在8月份，其中鄂托克旗站5月份降水趨勢明顯增大，其他站點5月降水趨勢相對穩定，直到進入7月份，三站降水趨勢顯著增大，流域進入汛期。海流兔河流域降水在年內分布集中，使年際間降水特征區分明顯，將有助于年際降水規律的分析。

圖2 氣象站1959-2019年月平均降水量

2.2 降水量年際特征分析

本文利用Morlet小波分析以16 a為時間尺度對3個氣象站1959—2019年的年降水量進行了周期分析(圖3)。鄂托克旗站的顯著性周期比橫山站和榆林站明顯，而橫山站與榆林站的差異性比較一致。鄂托克旗站在5～6 a的時間尺度和11～12 a的時間尺度信號能量變化較為強烈，干濕變化明顯；在5～6 a的尺度下共有8次干濕交替，2019年后的交替循環還未結束，表明該站未來的降水量是呈增加狀態的；在11～12 a的尺度下共有4次干濕交替，2019年后的交替循環還未閉合，大時間尺度上顯示該站還處于降水量少的區間內。榆林站和橫山站的小波實部和小波方差圖基本類似，干濕尺度在16 a以下的時間尺度上不甚明顯，在6～7 a的尺度下，主要發生在1959—1983年和1998—2019年度，其余信號能量變化較弱；在16 a以上的尺度下，橫山站和榆林站的干濕交替循環還未閉合，大時間尺度下這兩站未來還處在降水量少的時期內。

圖3 降水量小波實部圖及小波方差

通過經驗模態分解法(EMD)，對3個氣象站61 a的年降水量序列進行了分解，為了保證降水量信息的信號強度，均得到了方差貢獻率最大的3個IMF分量和1個趨勢分量(RES)，各分量表示的是不同時間尺度下的震蕩周期(表1)。鄂托克旗站IMF1分量的波動周期為3～6 a，1977年前的波動幅度較大，1978—1982年波動幅度較小， 1983—2019年波動幅度總體穩定，2019年之后未來2 a的降水趨勢是短幅下降后上升；IMF2分量的波動周期為7～12 a，1988年前的波動幅度較為明顯，1989—2015年波動幅度衰弱，2016—2019年波動幅度略有增大，2019年之后幾年的降水趨勢呈明顯下降狀態；IMF3分量的波動周期為30 a左右，2019年后的波動幅度大于前期水平，未來多年的降水量值將維持在2019年水平左右；RSE趨勢分量從20世紀60年代，降水量幅度處于歷史高位，未來該站降水整體仍將處于同位水平。橫山站IMF1分量的波動周期為3～7 a，1968年前的波動幅度較大，1969—2000年波動幅度逐漸衰弱，波動周期較短，2010—2019年振幅增大，波動周期較長，2019年之后的2～3 a的降水趨勢將會是短幅下降后上升；IMF2分量的波動周期為4～15 a，1988年前波動幅度明顯，1989—2000年波動幅度衰弱，2001年后波動周期增大，2019年之后的幾年將在高降水量持續一段時間后開始下降；IMF3分量的波動周期為35 a左右，1978—1988年處于波峰，1998—2008年處于波谷，2018年后期開始進入波峰時期，預計未來多年的降水量呈增長狀態；RSE趨勢分量在1988年處于波谷最低點，預計2019年之后該站的降水整體將處于高位。榆林站IMF1分量的波動周期為3～6 a，1968年前的波動幅度較大，1969—2001年波動幅度減弱，波動周期較短，2002年后波動幅度平穩，波動周期增加2 a左右，2019年后的2～3 a的降水趨勢也將是在短幅度下降后上升；IMF2分量的波動周期為5～10 a，1975年前的波動幅度具有一致性，1978—1988年波動幅度明顯，1989—1998年和2006—2014年期間波動幅度衰弱，2015年后的波動幅度明顯增強，2019年之后幾年間的降水趨勢將在短幅下降后上升；IMF3分量的波動周期為15～20 a，1988年之前的波動幅度具有一致性，1989年之后波動衰弱，一直下降至2008年、2009年后振幅明顯增大，2016年達到波峰后開始下降，預計未來多年的降水量將處于下降狀態；RSE趨勢分量同橫山站類似，預計2019年之后該站的降水整體也是處于高位。

表1 氣象站各模態分量的方差貢獻率 %

以上3個氣象站的降水序列震蕩周期表明，各站方差貢獻率最高的IMF1分量波動周期類似，短期預測其降水趨勢均為短幅下降后上升，中長周期下IMF2和IMF3的波動周期和未來趨勢性具有差異性，三站趨勢分量RES對未來的降水展望均呈增加態勢。

2.3 降水量預測

根據以上3個氣象站降水量數據的周期性分析可以看到，鄂托克旗站的小波分析顯著性周期分別為12 a和6 a，橫山站和榆林站的小波分析顯著性周期為6 a；EMD分析表明在中長周期下，三站的波動周期均為12 a左右，基于此本文利用乘性Holt-Winters模型以12 a為周期對歷史降水量數據提取信息模擬后，進行了未來12 a的降水量預測分析(表2)。在預測之前本文將61 a歷史數據劃分為49 a的識別期和12 a的驗證期，對該模型進行了適用性評估，結果表明鄂托克旗站12 a的驗證模擬值相對誤差為4.2%，橫山站為16.45%，榆林站為8.73%，震蕩周期類似，模擬效果較為理想。最終對61 a的數據進行模擬后得到鄂托克旗站Holt-Winters模型歷史降水數據估計下的參數α，γ，β分別為0，0.52，0.22，模擬值的平均相對誤差為28.61%，模擬值比實際值的震蕩強度較為劇烈，表明基于現有的人類活動及全球氣候變化下，未來幾年的預測趨勢為輕微下降后上升。橫山站的α，γ，β分別為0.12，0.08，0.2，平均相對誤差為22.62%，預測值在1985—1998年的震蕩強度普遍較大，未來幾年的趨勢預測為上升。榆林站的α，γ，β分別為0.05，0.38，0.24，平均相對誤差為22.27%，預測值的極端低值明顯，未來降水量的預測趨勢為顯著上升。整體上3個氣象站的模擬值多集中在平均值附近，極值明顯程度高于實際值，模擬情況系統性偏低，預測精度還有很大的提升空間。

表2 Holt-Winters模型多年平均模擬結果

2.4 流域降水特征分析

海流兔河流域是鄂爾多斯剝蝕高原向陜北黃土高原過渡的洼地小流域，整個流域處在毛烏素沙地之上，為蒙陜煤炭開采區水土流失嚴重的典型小流域，為更好地描述及預測該流域的降水量情況，本文利用泰森多邊形法對氣象站點數據進行加權分配，分割得到鄂托克旗站對該流域的控制面積為10.41%，橫山站為89.58%，榆林站為0.01%，其中橫山站的降水數據為主要流域控制項。通過數據加權后的小波及EMD分析得到海流兔河流域的基礎干濕交替周期尺度為6 a，震蕩以1988—1998年最為不明顯，16 a以上的尺度均顯示流域處在降水量上升期，未來流域整體降水趨勢預測呈顯著上升狀態。乘性Holt-Winters模型模擬預側值與橫山站類似，α，γ，β分別為0.12，0.09，0.21，平均相對誤差為22.28%，模擬精度不太理想，但結合小波及DEM可預測基于現有人類活動及全球氣候變化的影響下未來12 a流域降水量整體趨勢呈增加態勢(表3)。

表3 Holt-Winters模型未來12 a預測值

3 結 論

1961—2019年的降水量數據表明在現有人類活動強度及全球氣候趨勢下，鄂托克旗站、橫山站和榆林站的降水特征呈上升趨勢，線性上升傾向率為：榆林站(1.192 2)>橫山站(0.073 1)>鄂托克旗站(0.003)，流域整體的上升傾向率為0.065 6，未來降水展望為增長態勢。根據小波分析、EMD分析和乘性Holt-Winters模型分析得到3個氣象站具有明顯的周期性和趨勢性，干濕周期以6 a和12 a為主；通過降水量變化的主周期外推得到，3個氣象站的短期降水預測均呈上升趨勢，但在中長期的時間尺度下略有不同，其中榆林站的波動程度顯著于鄂托克旗站和橫山站；綜合加權分析得到，目前海流兔河流域整體處于降水增長期，預計在現有人類活動及全球氣候變化的影響下未來12 a的降水量將繼續呈波動增長趨勢。

本文首次基于周期分析后利用Holt-Winters模型對多年降水量進行了模擬預測，得到鄂托克旗站的預測值呈波動平穩態勢，橫山站和榆林站的預測值呈波動增長態勢，流域整體降水量呈波動增長態勢，這與小波分析和EMD分析的結論相同；盡管Holt-Winters模型在該流域的模擬預測效果不太理想，但對于氣候周期性明顯和氣象資料缺乏地區的降水量統計預測提供了新思路。