气候变化和区外调水对米山水库水质的影响

张明芳，黄林显，刘 扬，张 杰，刘建光，宋协良，李文虎

(1. 威海市水文中心，山东 威海 264209；2. 济南大学 a. 水利与环境学院，b. 山东省地下水数值模拟与污染控制工程技术研究中心，山东 济南 250022；3. 山东省地质测绘院，山东 济南 250013；4. 威海市文登区水利事务服务中心，山东 威海 264400)

米山水库是山东省威海市境内最大的水库，也是威海市最重要的城市供水水源。米山水库作为南水北调工程山东段末端的调蓄水库，主要承担着威海市主城区和文登区的供水任务，其水质对城市供水安全和区域水生态稳定至关重要[1]。湖库水质主要受气候变化和人类活动的影响，其动态变化表现出复杂的趋势性、季节性及随机性等特征。对于米山水库来说，一方面随着全球性气候的变化，近年来胶东半岛连续多年降雨量减少，气候变暖，导致米山水库入库水量减少，对米山水库的水质产生影响；另一方面，米山水库自2015年12月21日开始区外调水，调水水源包括黄河水、长江水和当地水，外调水同样对米山水库的水质产生了显著影响。文献[2]中的研究已经指出，米山水库2018年水体中的电导率、钾、钠、钙、氯化物、硫酸盐、重碳酸盐等指标与2015年均值都有较大幅度的增长[2]。通过定量评估气候变化和区外调水对米山水库水质演变的贡献程度并评价在不同时期对水质变化的作用过程，可以为米山水库建立更加全面的水质安全风险管理体系，实施更加科学的水量-水质联合调控决策以及维护供水安全提供重要借鉴和参考。

前人关于调水工程对湖库水质影响的研究大多集中在水质评价[3]、风险评估[4-5]和调控政策[6]等方面，从定量角度评价气候变化和区外调水对湖库水质演变贡献度的研究尚未见报道。气候变化和区外调水贡献度的定量评价需要科学选择水体水质的度量指标。在水体的众多指标中，水的电导率主要用来表征水体传导电流能力的强弱，通过电导率值的大小可以推测水中带电荷物质的总浓度，因此可以作为反映水体污染情况的一个重要度量[7]。掌握电导率在水体中的时空分布对进行湖库生态环境变化趋势分析具有重要的参考价值。

本文中首先利用局部加权回归的时间序列分解(STL)法[8-10]对降雨及电导率时间序列进行时序分解，揭示它们的变化趋势、周期长度、随机波动范围等，分别提取出电导率的自然状态(受气候变化影响)和混合状态(同时受气候变化和区外调水影响)时间序列。利用自然状态时间序列通过灰色系统GM(1,1)分解模型建立电导率模拟模型，得到自然状态下的电导率模拟值，通过与混合状态下实测电导率的比较，定量评价气候变化和区外调水对米山水库水质演变的贡献程度。

1 数据序列和研究方法

米山水库位于山东省威海市文登区米山镇，是以防洪、城镇供水、发电等综合利用的大(二)型水库，控制流域面积为440 km2，总库容为2.8×108m3，死库容为5.07×106m3。2015年12月22日南水北调东线黄河水入威开始通水，2016年3月10日，长江水首次调入米山水库，威海市实现长江水、黄河水和当地水多水源联合调度。米山水库流域及雨量站分布如图1所示。

图1 米山水库流域及雨量站分布

1.1 数据序列

1.1.1 数据来源

受制于月数据收集的困难，本研究中STL时间序列分解采用1980—2019年月均降雨量数据和2012—2019年月均电导率数据(见图2)，通过对时间序列进行统计分析来提取历史数据信息，以此找出其动态变化规律，揭示行为的未来演变趋势；其中2013年电导率数据部分缺失，采用样条插值进行插补(见图2)。GM(1,1)分解模型采用1987—2019年年均电导率数据，通过建立电导率模拟模型，对电导率进行模拟计算。降雨数据来自于中国气象科学数据共享服务网，电导率数据来自于威海市水文中心监测数据。

1.1.2 降雨及电导率统计特征

对米山水库1980—2019年年均降雨量和1987—2019年年均电导率长时间序列进行统计分析，结果如表1所示。从表中可以看出，研究区多年年均降雨量为811.67 mm，降雨较为丰富；年均降雨量的标准差为205.66 mm，说明研究区降雨的波动性相对较为平稳；年均电导率为374.98 μS/cm，说明水体中离子含量较多；年均电导率的标准差为215.02 μS/cm，说明电导率波动性较大。

(a)1980—2019年月均降雨量

(b)2012—2019年月均电导率图2 米山水库月均降雨量和月均电导率及样条插值数据

表1 米山水库年均降雨量和电导率时间序列数据统计特征

1.2 研究方法

1.2.1 STL法

STL法是一种非参数统计方法，它以具有鲁棒性的局部加权回归散点平滑(LOESS)作为平滑方法，能够将时间序列以加法模型为原理分解为趋势变化、季节性变化及不规则的随机波动项。主要计算步骤参见文献[8-9]，计算公式为

Yt=Tt+St+Rt，

(1)

式中：Yt是t时刻的观测值；Tt是t时刻的趋势变化，反映时间序列的变化趋势，如增加或者减少；St是t时刻的季节性变化，反映时间序列的波动周期，如12个月；Rt是t时刻的随机波动项，反映受外部因素影响的强弱。

1.2.2 GM(1,1)分解模型

目前常用的水质模拟预测方法主要有数理统计法[11-12]、机器学习模型法[13-16]、数值模型法[17-18]和灰色模型法[19-22]等，其中灰色模型法中的GM(1,1)模型因用法简单且所需实测数据少而被广泛应用，但是存在灰色偏差和抗干扰能力弱的局限性，只适合于对呈指数规律变化且增长速率较低的短时间序列进行预测[23]。水质、水量等水文要素的增长速率有快有慢，并且数据本身存在周期变化和随机扰动，造成传统GM(1,1)模型对于水文要素的预测精度不高。为了解决此类问题，Hao等[24]提出了GM(1,1)分解模型，并将其应用于水文时间序列的模拟预测，取得了非常理想的预测效果。GM(1,1)分解模型将时间序列分解为趋势项、周期波动项和随机变动项，因此能较好地拟合水文要素的动态变化特征。李华敏等[25]将GM(1,1)分解模型应用于山西柳林泉水流量的预测中，并把预测结果与自回归差分移动平均 (autoregressive integrated moving average, ARIMA)模型进行对比，结果显示通过周期修正能显著提高GM(1,1)分解模型预测精度，其预测结果要好于ARIMA模型的。

GM(1,1)分解模型的工作原理分述如下。

1)时间序列的趋势项。利用GM(1,1)模型模拟计算时间序列的趋势项。首先对原始时间序列x(0)(t)采用累加的方式生成一组趋势明显的新时间序列x(1)(t)，然后建立微分方程

(2)

式中：a为系统的发展系数；u为灰色作用量。

(3)

然后，进行累减操作恢复原始数据序列，得到预测结果如下：

(4)

2)时间序列的周期波动项。周期波动项通过利用正弦函数对残差序列进行拟合获得。残差序列ε(t)计算公式如下：

(5)

文献[24]中对残差序列采用分段正弦曲线拟合的方法，本文中为了避免人为选择分段周期和振幅所带来的误差，采用带有截距的正弦曲线进行拟合，

(6)

(7)

3)时间序列的随机波动项。时间序列的随机波动项R(t)通过以下公式计算获得：

(8)

合并趋势项、周期波动项和随机波动项，最终得到对时间序列动态变化的模拟预测值Q(t)，

(9)

2 实验结果与分析

2.1 研究方法

利用STL法对1980—2019年月降雨量及2012—2019年月电导率值进行时序分解，揭示它们的变化趋势、周期长度、随机波动范围等，并进一步提取出电导率的自然状态(受气候变化影响)和混合状态(同时受气候变化和区外调水影响)时间序列。

对降雨时间序列进行STL分解，获取降雨量趋势变化、季节性变化及随机波动项，结果见图3。从图中的趋势变化可以看出，降雨量从2012年开始有减少的趋势，整体呈现稳定中略有波动的状态；降雨时间序列表现出显著的周期性(周期为12个月)，且从2009年开始周期内的季节性变化均呈递减的趋势；从随机波动项可以看出，残差波动程度在1998—2013年的时间序列要大于其他时间序列。

图3 1980—2019年米山水库月均降雨量时间序列分解结果

对月均电导率时间序列进行STL分解，获取其趋势变化、季节性变化及随机波动项，结果见图4。从图4(a)2012—2019年月均电导率STL分解结果可以看出，电导率值在2016年之前呈现稳定中略有波动的态势，电导率波动范围为250～350 μS/cm，此时电导率主要受气候变化影响，即处于自然状态。在水库2016年开始调水后电导率迅速增大并在2017年达到最大值，波动范围为857～977 μS/cm，之后处于稳定波动状态，此时电导率同时受气候变化和区外调水的影响，即电导率处于混合状态。电导率值的波动呈现出明显的周期性，周期为12个月。从随机波动项可以看出，残差的波动程度在2016年后的时间序列要明显大于之前的时间序列，说明米山水库水质在2016年实施区外调水后受人类活动等外部因素的干扰有所增加。为了更清楚地了解自然状态下电导率的变化特征，对2012—2015年月电导率时间序列进行STL分解，如图4(b)所示。由图可以看出，自然状态下电导率处于平稳波动状态数值的波动同样呈现12个月的周期性，且随机波动项的波动程度在2013—2014年之间的时间序列要大于其他的时间序列。

(a)2012—2019年月均电导率

(b)2012—2015年月均电导率图4 米山水库月均电导率的时间序列分解结果

2.2 自然状态下电导率与降雨关系分析

对自然状态下电导率与降雨的关系进行分析。1987—2015年年均电导率与年均降雨量变化关系图如图5所示。由图可以看出，电导率值与降雨量呈现明显的负相关关系，当年均降雨量较大时，因降雨对水体的带电荷物质的稀释作用而导致电导率年均值相应的较小： 1999年研究区降雨量出现较小值(360.6 mm)，米山水库电导率达到347 μS/cm，为历年最大值；1990年降雨量出现较大值，为1 180.9 mm，电导率仅为232 μS/cm，为历年最小值。由此可见，降雨对米山水库电导率的稀释作用较为明显，进一步说明自然状态下电导率主要受气候变化的影响。

图5 1987—2015年米山水库年均电导率与降雨量变化关系

2.3 米山水库电导率模拟计算

通过上述分析可以看出，米山水库电导率在2016年前处于自然状态，因此可以利用1987—2015年电导率时间序列通过GM(1,1)分解模型建立自然状态电导率模型；运用该模型外推，获得2016—2019年自然状态下的电导率，利用同时间序列实测电导率减去自然状态电导率，获得区外调水对电导率增加的贡献率，从而分别获得气候变化和区外调水对电导率增加的影响。

2.3.1 电导率的趋势项

利用1987—2015年自然状态下电导率值作为原始数据，建立GM(1,1)模型并求得模拟模型为

(10)

2.3.2 电导率的周期波动项

GM(1,1)模型的残差值可以由式(5)求得，如图7所示。利用式(6)带截距的正弦曲线对残差值进行拟合，得到残差周期拟合曲线的函数表达式如下：

图6 灰色系统GM(1,1)模型对米山水库电导率模拟值与实测值

(11)

由图可以看出，所求得的正弦残差周期拟合曲线与残差值整体拟合较好，拟合精度较高，能够反映残差的周期性变化。

2.3.3 电导率的随机波动项

根据式(9)，合并电导率趋势项、周期波动项和随机波动变化项，可以对自然状态下米山水库电导率动态变化进行模拟预测，结果见图8。由图可以看出：米山水库自然状态下电导率值从1987年开始不断增大，到2002年达到最高值；2002—2015年电导率值呈现减小的趋势，从2016年开始电导率值又有增大的趋势。同时，电导率模拟值与实际观测值整体拟合较好，绝大部分观测值均在模拟值的上、下限范围内。此外，模拟值与实测值一样，都体现了明显的上升和回落的周期性变化，说明GM(1,1)分解模型的模拟精度较高。

图7 米山水库电导率残差周期拟合曲线、随机波动曲线及残差值

2.3.4 气候变化和区外调水贡献度定量评价

利用米山水库2016—2019年电导率实测值减去自然状态下的电导率值，即可得到区外调水造成的电导率变化值，结果见表2。从表中数据可以看出，2016、2017、2018、2019年自然状态下米山水库的电导率分别为269.3～317.26、270.94～318.9、273.4～321.36、276.63～324.59 μS/cm，处于正常的波动范围内；区外调水造成的电导率的增加量分别为540.41～588.37、643.70～691.66、656.18～704.14、629.74～677.70 μS/cm，分别是自然状态下电导率的1.85～2.01、2.17～2.38、2.19～2.40、2.09～2.28倍，说明区外调水对米山水库电导率的增加起着主要控制作用。通过对造成米山水库电导率变化的主控因素进行定量研究表明，影响米山水库水质变化的主要原因不仅仅是气候变化，这一结论对于米山水库建立更加全面的水质安全风险管理体系，实施更加科学的水量-水质联合调控决策，维护供水安全具有重要意义。

图8 灰色系统GM(1,1)分解模型对米山水库电导率的模拟值与实测值

表2 气候变化和区外调水对米山水库电导率的影响

2.3.5 分析讨论

本文中利用GM(1,1)分解模型建立了天然状态下电导率的模拟预测模型，分别对气候变化和区外调水对米山水库电导率的影响进行了定量评价，研究结果具有重要的参考价值，但仍然存在着以下不足：

1)米山水库从2016年才开始实施区外调水，较短的区外调水时间序列不能充分反映区外调水影响下电导率的变化。

2)GM(1,1)分解模型虽然具有用法简单、所需实测数据少和比传统的GM(1, 1)模型模拟精度更高的优势；但湖库水体中的污染物浓度受到多种因素的影响(如降雨、气温和调水量等)，因此水质预测模型涉及多种因素的非线性关系。GM(1, 1)分解模型无法充分考虑多种因素的影响，造成一定误差。

3)在今后的研究中一方面要增加电导率时间序列的长度，另一方面要采用长短期记忆(LSTM)神经网络等多参数机器学习方法，通过输入多个影响变量(如降雨量、气温和调水量等)进行水质预测，由于输入的不同变量之间可以相互验证、相互作用，因此能够有效提高预测精度。

3 结论

本文中首先利用STL法对米山水库降雨量与电导率进行时间序列分解，从宏观角度揭示它们的趋势变化、周期波动长度和随机波动范围，并提取出电导率的自然状态和混合状态时间序列；在此基础上，通过GM(1,1)分解模型建立自然状态下电导率模拟模型，并将其外推，定量评价气候变化和区外调水对米山水库水质演变的贡献，得出以下结论：

1)通过STL分解结果可以看出，米山水库降雨时间序列虽然从2012年开始有减少的趋势，但整体上仍然呈现稳定中略有波动的状态；降雨表现出显著的周期性。电导率在2016年之前同样呈现稳定中略有波动的状态，此时电导率主要受气候变化影响；2016年区外调水后电导率快速增大并在2017年达到最大值，之后处于稳定波动状态，此时电导率同时受气候变化和区外调水的影响。电导率值的波动同样呈现明显的周期性。

2)GM(1,1)分解模型的电导率模拟值与实测值整体拟合较好，绝大部分实测值均在模拟值的上、下限范围内，且模拟值与实测值都能反映明显的上升和回落的周期性变化，说明GM(1,1)分解模型的模拟精度较高。

3)模拟分析结果显示，2016—2019年区外调水造成米山水库电导率的增加量分别为540.41～588.37、643.7～691.66、656.18～704.14、629.74～677.7 μS/cm，分别是自然状态下电导率的1.85～2.01、2.17～2.38、2.19～2.40、2.09～2.28倍，说明区外调水对米山水库电导率的增加起到主要控制作用。

4)由于米山水库水质的演变主要受区外调水量、客水水质影响，因此应建立更加全面的区外调水水质安全风险管理体系，实施更加科学的水量-水质联合调控决策，以此维护米山水库的供水安全和水生态稳定。