氣候變化和區外調水對米山水庫水質的影響

張明芳，黃林顯，劉 揚，張 杰，劉建光，宋協良，李文虎

(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)由于米山水庫水質的演變主要受區外調水量、客水水質影響，因此應建立更加全面的區外調水水質安全風險管理體系，實施更加科學的水量-水質聯合調控決策，以此維護米山水庫的供水安全和水生態穩定。