基于时间序列的海流兔河流域降水量模拟预测分析

周生辉, 刘廷玺,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模型在该流域的模拟预测效果不太理想，但对于气候周期性明显和气象资料缺乏地区的降水量统计预测提供了新思路。