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

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 水利與土木建筑工程學(xué)院, 呼和浩特 010018; 2.內(nèi)蒙古自治區(qū)水資源保護(hù)與利用重點(diǎn)實(shí)驗(yàn)室, 呼和浩特 010018)

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

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

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況及數(shù)據(jù)來(lái)源

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

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

1.2 研究方法

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

式中:Ai為第i個(gè)子區(qū)域的面積;αi為第i個(gè)雨量站的面積權(quán)重;Pi為第i個(gè)子雨量站的降水量;n為雨量站數(shù)目。

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

小波方差的計(jì)算公式為：

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

1.2.3 EMD經(jīng)驗(yàn)?zāi)B(tài)分解 EMD經(jīng)驗(yàn)?zāi)B(tài)分解可對(duì)一個(gè)時(shí)間信號(hào)將其不同尺度即頻率的波動(dòng)或趨勢(shì)逐級(jí)分解開來(lái)，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，稱為本征模函數(shù)(Intrinsic Mode Function, IMF)，它是目前處理非平穩(wěn)、非線性信號(hào)，特別是分析時(shí)間序列趨勢(shì)的最好方法[27-28]，IMF分量可以反映出降水量序列中不同特征尺度的振蕩。其主要的計(jì)算原理如下：

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

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

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

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

(4) 若所得新序列I1(t)滿足兩個(gè)條件：極值點(diǎn)與過(guò)零點(diǎn)數(shù)目相等或相差一個(gè)，包絡(luò)線的均值為0，則I1(t)為固有模態(tài)函數(shù)(IMF)。

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

EMD方法在篩分過(guò)程中體現(xiàn)出了強(qiáng)大的直接性及自適應(yīng)性，可根據(jù)實(shí)際序列的數(shù)理特征自主形成特殊的基函數(shù)，每一個(gè)分離出來(lái)的 IMF都具有一定的物理意義，都是對(duì)原始序列物理信息的真實(shí)反映[29]。

1.2.4 Holt-Winters模型 對(duì)于未來(lái)發(fā)生的事情，最新觀察值較早期觀察值包含更多的信息，因而在預(yù)測(cè)時(shí)，最新觀測(cè)值較早期觀察值具有更大的權(quán)重，Holt-Winters模型具有周期性特征，可對(duì)時(shí)間序列進(jìn)行短期模擬預(yù)測(cè)[30]。Holt-Winters模型本質(zhì)上是一種高級(jí)指數(shù)平滑形式模型，具有周期調(diào)整與長(zhǎng)期趨勢(shì)調(diào)整特性，能對(duì)兼有長(zhǎng)期趨勢(shì)和季節(jié)模式的數(shù)據(jù)進(jìn)行預(yù)測(cè)[31-33]。本文利用的Holt-Winters乘性模型的基本原理方程如下：

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

Holt-Winters模型預(yù)測(cè)值計(jì)算為：

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

式中:St為平滑值即水平分量；α為水平權(quán)重;bt為長(zhǎng)期趨勢(shì)值；γ為趨勢(shì)權(quán)重;It為季節(jié)分量；β為季節(jié)權(quán)重;L為季節(jié)長(zhǎng)度(每年的月數(shù)或季數(shù));t為當(dāng)前時(shí)間;Xt為實(shí)際觀測(cè)值;yt+s為預(yù)測(cè)值；k為預(yù)測(cè)超前期數(shù);其中的γ,α,β范圍為0～1。

Holt-Winters乘性模型中的γ，α，β取值依賴于已知時(shí)間序列的性質(zhì)，通常為0.1～0.3的數(shù)值并產(chǎn)生一個(gè)依賴于大量過(guò)去觀測(cè)資料的預(yù)測(cè)[34]。

2 結(jié)果與分析

2.1 降水量年內(nèi)特征分析

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

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

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

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

2.2 降水量年際特征分析

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

圖3 降水量小波實(shí)部圖及小波方差

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

表1 氣象站各模態(tài)分量的方差貢獻(xiàn)率 %

以上3個(gè)氣象站的降水序列震蕩周期表明，各站方差貢獻(xiàn)率最高的IMF1分量波動(dòng)周期類似，短期預(yù)測(cè)其降水趨勢(shì)均為短幅下降后上升，中長(zhǎng)周期下IMF2和IMF3的波動(dòng)周期和未來(lái)趨勢(shì)性具有差異性，三站趨勢(shì)分量RES對(duì)未來(lái)的降水展望均呈增加態(tài)勢(shì)。

2.3 降水量預(yù)測(cè)

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

表2 Holt-Winters模型多年平均模擬結(jié)果

2.4 流域降水特征分析

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

表3 Holt-Winters模型未來(lái)12 a預(yù)測(cè)值

3 結(jié) 論

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

本文首次基于周期分析后利用Holt-Winters模型對(duì)多年降水量進(jìn)行了模擬預(yù)測(cè)，得到鄂托克旗站的預(yù)測(cè)值呈波動(dòng)平穩(wěn)態(tài)勢(shì)，橫山站和榆林站的預(yù)測(cè)值呈波動(dòng)增長(zhǎng)態(tài)勢(shì)，流域整體降水量呈波動(dòng)增長(zhǎng)態(tài)勢(shì)，這與小波分析和EMD分析的結(jié)論相同；盡管Holt-Winters模型在該流域的模擬預(yù)測(cè)效果不太理想，但對(duì)于氣候周期性明顯和氣象資料缺乏地區(qū)的降水量統(tǒng)計(jì)預(yù)測(cè)提供了新思路。