氣候變化和區(qū)外調(diào)水對米山水庫水質(zhì)的影響

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

(1. 威海市水文中心，山東 威海 264209；2. 濟南大學(xué) a. 水利與環(huán)境學(xué)院，b. 山東省地下水?dāng)?shù)值模擬與污染控制工程技術(shù)研究中心，山東 濟南 250022；3. 山東省地質(zhì)測繪院，山東 濟南 250013；4. 威海市文登區(qū)水利事務(wù)服務(wù)中心，山東 威海 264400)

米山水庫是山東省威海市境內(nèi)最大的水庫，也是威海市最重要的城市供水水源。米山水庫作為南水北調(diào)工程山東段末端的調(diào)蓄水庫，主要承擔(dān)著威海市主城區(qū)和文登區(qū)的供水任務(wù)，其水質(zhì)對城市供水安全和區(qū)域水生態(tài)穩(wěn)定至關(guān)重要[1]。湖庫水質(zhì)主要受氣候變化和人類活動的影響，其動態(tài)變化表現(xiàn)出復(fù)雜的趨勢性、季節(jié)性及隨機性等特征。對于米山水庫來說，一方面隨著全球性氣候的變化，近年來膠東半島連續(xù)多年降雨量減少，氣候變暖，導(dǎo)致米山水庫入庫水量減少，對米山水庫的水質(zhì)產(chǎn)生影響；另一方面，米山水庫自2015年12月21日開始區(qū)外調(diào)水，調(diào)水水源包括黃河水、長江水和當(dāng)?shù)厮?，外調(diào)水同樣對米山水庫的水質(zhì)產(chǎn)生了顯著影響。文獻(xiàn)[2]中的研究已經(jīng)指出，米山水庫2018年水體中的電導(dǎo)率、鉀、鈉、鈣、氯化物、硫酸鹽、重碳酸鹽等指標(biāo)與2015年均值都有較大幅度的增長[2]。通過定量評估氣候變化和區(qū)外調(diào)水對米山水庫水質(zhì)演變的貢獻(xiàn)程度并評價在不同時期對水質(zhì)變化的作用過程，可以為米山水庫建立更加全面的水質(zhì)安全風(fēng)險管理體系，實施更加科學(xué)的水量-水質(zhì)聯(lián)合調(diào)控決策以及維護供水安全提供重要借鑒和參考。

前人關(guān)于調(diào)水工程對湖庫水質(zhì)影響的研究大多集中在水質(zhì)評價[3]、風(fēng)險評估[4-5]和調(diào)控政策[6]等方面，從定量角度評價氣候變化和區(qū)外調(diào)水對湖庫水質(zhì)演變貢獻(xiàn)度的研究尚未見報道。氣候變化和區(qū)外調(diào)水貢獻(xiàn)度的定量評價需要科學(xué)選擇水體水質(zhì)的度量指標(biāo)。在水體的眾多指標(biāo)中，水的電導(dǎo)率主要用來表征水體傳導(dǎo)電流能力的強弱，通過電導(dǎo)率值的大小可以推測水中帶電荷物質(zhì)的總濃度，因此可以作為反映水體污染情況的一個重要度量[7]。掌握電導(dǎo)率在水體中的時空分布對進行湖庫生態(tài)環(huán)境變化趨勢分析具有重要的參考價值。

本文中首先利用局部加權(quán)回歸的時間序列分解(STL)法[8-10]對降雨及電導(dǎo)率時間序列進行時序分解，揭示它們的變化趨勢、周期長度、隨機波動范圍等，分別提取出電導(dǎo)率的自然狀態(tài)(受氣候變化影響)和混合狀態(tài)(同時受氣候變化和區(qū)外調(diào)水影響)時間序列。利用自然狀態(tài)時間序列通過灰色系統(tǒng)GM(1,1)分解模型建立電導(dǎo)率模擬模型，得到自然狀態(tài)下的電導(dǎo)率模擬值，通過與混合狀態(tài)下實測電導(dǎo)率的比較，定量評價氣候變化和區(qū)外調(diào)水對米山水庫水質(zhì)演變的貢獻(xiàn)程度。

1 數(shù)據(jù)序列和研究方法

米山水庫位于山東省威海市文登區(qū)米山鎮(zhèn)，是以防洪、城鎮(zhèn)供水、發(fā)電等綜合利用的大(二)型水庫，控制流域面積為440 km2，總庫容為2.8×108m3，死庫容為5.07×106m3。2015年12月22日南水北調(diào)東線黃河水入威開始通水，2016年3月10日，長江水首次調(diào)入米山水庫，威海市實現(xiàn)長江水、黃河水和當(dāng)?shù)厮嗨绰?lián)合調(diào)度。米山水庫流域及雨量站分布如圖1所示。

圖1 米山水庫流域及雨量站分布

1.1 數(shù)據(jù)序列

1.1.1 數(shù)據(jù)來源

受制于月數(shù)據(jù)收集的困難，本研究中STL時間序列分解采用1980—2019年月均降雨量數(shù)據(jù)和2012—2019年月均電導(dǎo)率數(shù)據(jù)(見圖2)，通過對時間序列進行統(tǒng)計分析來提取歷史數(shù)據(jù)信息，以此找出其動態(tài)變化規(guī)律，揭示行為的未來演變趨勢；其中2013年電導(dǎo)率數(shù)據(jù)部分缺失，采用樣條插值進行插補(見圖2)。GM(1,1)分解模型采用1987—2019年年均電導(dǎo)率數(shù)據(jù)，通過建立電導(dǎo)率模擬模型，對電導(dǎo)率進行模擬計算。降雨數(shù)據(jù)來自于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)，電導(dǎo)率數(shù)據(jù)來自于威海市水文中心監(jiān)測數(shù)據(jù)。

1.1.2 降雨及電導(dǎo)率統(tǒng)計特征

對米山水庫1980—2019年年均降雨量和1987—2019年年均電導(dǎo)率長時間序列進行統(tǒng)計分析，結(jié)果如表1所示。從表中可以看出，研究區(qū)多年年均降雨量為811.67 mm，降雨較為豐富；年均降雨量的標(biāo)準(zhǔn)差為205.66 mm，說明研究區(qū)降雨的波動性相對較為平穩(wěn)；年均電導(dǎo)率為374.98 μS/cm，說明水體中離子含量較多；年均電導(dǎo)率的標(biāo)準(zhǔn)差為215.02 μS/cm，說明電導(dǎo)率波動性較大。

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

(b)2012—2019年月均電導(dǎo)率圖2 米山水庫月均降雨量和月均電導(dǎo)率及樣條插值數(shù)據(jù)

表1 米山水庫年均降雨量和電導(dǎo)率時間序列數(shù)據(jù)統(tǒng)計特征

1.2 研究方法

1.2.1 STL法

STL法是一種非參數(shù)統(tǒng)計方法，它以具有魯棒性的局部加權(quán)回歸散點平滑(LOESS)作為平滑方法，能夠?qū)r間序列以加法模型為原理分解為趨勢變化、季節(jié)性變化及不規(guī)則的隨機波動項。主要計算步驟參見文獻(xiàn)[8-9]，計算公式為

Yt=Tt+St+Rt，

(1)

式中：Yt是t時刻的觀測值；Tt是t時刻的趨勢變化，反映時間序列的變化趨勢，如增加或者減少；St是t時刻的季節(jié)性變化，反映時間序列的波動周期，如12個月；Rt是t時刻的隨機波動項，反映受外部因素影響的強弱。

1.2.2 GM(1,1)分解模型

目前常用的水質(zhì)模擬預(yù)測方法主要有數(shù)理統(tǒng)計法[11-12]、機器學(xué)習(xí)模型法[13-16]、數(shù)值模型法[17-18]和灰色模型法[19-22]等，其中灰色模型法中的GM(1,1)模型因用法簡單且所需實測數(shù)據(jù)少而被廣泛應(yīng)用，但是存在灰色偏差和抗干擾能力弱的局限性，只適合于對呈指數(shù)規(guī)律變化且增長速率較低的短時間序列進行預(yù)測[23]。水質(zhì)、水量等水文要素的增長速率有快有慢，并且數(shù)據(jù)本身存在周期變化和隨機擾動，造成傳統(tǒng)GM(1,1)模型對于水文要素的預(yù)測精度不高。為了解決此類問題，Hao等[24]提出了GM(1,1)分解模型，并將其應(yīng)用于水文時間序列的模擬預(yù)測，取得了非常理想的預(yù)測效果。GM(1,1)分解模型將時間序列分解為趨勢項、周期波動項和隨機變動項，因此能較好地擬合水文要素的動態(tài)變化特征。李華敏等[25]將GM(1,1)分解模型應(yīng)用于山西柳林泉水流量的預(yù)測中，并把預(yù)測結(jié)果與自回歸差分移動平均 (autoregressive integrated moving average, ARIMA)模型進行對比，結(jié)果顯示通過周期修正能顯著提高GM(1,1)分解模型預(yù)測精度，其預(yù)測結(jié)果要好于ARIMA模型的。

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

1)時間序列的趨勢項。利用GM(1,1)模型模擬計算時間序列的趨勢項。首先對原始時間序列x(0)(t)采用累加的方式生成一組趨勢明顯的新時間序列x(1)(t)，然后建立微分方程

(2)

式中：a為系統(tǒng)的發(fā)展系數(shù)；u為灰色作用量。

(3)

然后，進行累減操作恢復(fù)原始數(shù)據(jù)序列，得到預(yù)測結(jié)果如下：

(4)

2)時間序列的周期波動項。周期波動項通過利用正弦函數(shù)對殘差序列進行擬合獲得。殘差序列ε(t)計算公式如下：

(5)

文獻(xiàn)[24]中對殘差序列采用分段正弦曲線擬合的方法，本文中為了避免人為選擇分段周期和振幅所帶來的誤差，采用帶有截距的正弦曲線進行擬合，

(6)

(7)

3)時間序列的隨機波動項。時間序列的隨機波動項R(t)通過以下公式計算獲得：

(8)

合并趨勢項、周期波動項和隨機波動項，最終得到對時間序列動態(tài)變化的模擬預(yù)測值Q(t)，

(9)

2 實驗結(jié)果與分析

2.1 研究方法

利用STL法對1980—2019年月降雨量及2012—2019年月電導(dǎo)率值進行時序分解，揭示它們的變化趨勢、周期長度、隨機波動范圍等，并進一步提取出電導(dǎo)率的自然狀態(tài)(受氣候變化影響)和混合狀態(tài)(同時受氣候變化和區(qū)外調(diào)水影響)時間序列。

對降雨時間序列進行STL分解，獲取降雨量趨勢變化、季節(jié)性變化及隨機波動項，結(jié)果見圖3。從圖中的趨勢變化可以看出，降雨量從2012年開始有減少的趨勢，整體呈現(xiàn)穩(wěn)定中略有波動的狀態(tài)；降雨時間序列表現(xiàn)出顯著的周期性(周期為12個月)，且從2009年開始周期內(nèi)的季節(jié)性變化均呈遞減的趨勢；從隨機波動項可以看出，殘差波動程度在1998—2013年的時間序列要大于其他時間序列。

圖3 1980—2019年米山水庫月均降雨量時間序列分解結(jié)果

對月均電導(dǎo)率時間序列進行STL分解，獲取其趨勢變化、季節(jié)性變化及隨機波動項，結(jié)果見圖4。從圖4(a)2012—2019年月均電導(dǎo)率STL分解結(jié)果可以看出，電導(dǎo)率值在2016年之前呈現(xiàn)穩(wěn)定中略有波動的態(tài)勢，電導(dǎo)率波動范圍為250～350 μS/cm，此時電導(dǎo)率主要受氣候變化影響，即處于自然狀態(tài)。在水庫2016年開始調(diào)水后電導(dǎo)率迅速增大并在2017年達(dá)到最大值，波動范圍為857～977 μS/cm，之后處于穩(wěn)定波動狀態(tài)，此時電導(dǎo)率同時受氣候變化和區(qū)外調(diào)水的影響，即電導(dǎo)率處于混合狀態(tài)。電導(dǎo)率值的波動呈現(xiàn)出明顯的周期性，周期為12個月。從隨機波動項可以看出，殘差的波動程度在2016年后的時間序列要明顯大于之前的時間序列，說明米山水庫水質(zhì)在2016年實施區(qū)外調(diào)水后受人類活動等外部因素的干擾有所增加。為了更清楚地了解自然狀態(tài)下電導(dǎo)率的變化特征，對2012—2015年月電導(dǎo)率時間序列進行STL分解，如圖4(b)所示。由圖可以看出，自然狀態(tài)下電導(dǎo)率處于平穩(wěn)波動狀態(tài)數(shù)值的波動同樣呈現(xiàn)12個月的周期性，且隨機波動項的波動程度在2013—2014年之間的時間序列要大于其他的時間序列。

(a)2012—2019年月均電導(dǎo)率

(b)2012—2015年月均電導(dǎo)率圖4 米山水庫月均電導(dǎo)率的時間序列分解結(jié)果

2.2 自然狀態(tài)下電導(dǎo)率與降雨關(guān)系分析

對自然狀態(tài)下電導(dǎo)率與降雨的關(guān)系進行分析。1987—2015年年均電導(dǎo)率與年均降雨量變化關(guān)系圖如圖5所示。由圖可以看出，電導(dǎo)率值與降雨量呈現(xiàn)明顯的負(fù)相關(guān)關(guān)系，當(dāng)年均降雨量較大時，因降雨對水體的帶電荷物質(zhì)的稀釋作用而導(dǎo)致電導(dǎo)率年均值相應(yīng)的較?。?1999年研究區(qū)降雨量出現(xiàn)較小值(360.6 mm)，米山水庫電導(dǎo)率達(dá)到347 μS/cm，為歷年最大值；1990年降雨量出現(xiàn)較大值，為1 180.9 mm，電導(dǎo)率僅為232 μS/cm，為歷年最小值。由此可見，降雨對米山水庫電導(dǎo)率的稀釋作用較為明顯，進一步說明自然狀態(tài)下電導(dǎo)率主要受氣候變化的影響。

圖5 1987—2015年米山水庫年均電導(dǎo)率與降雨量變化關(guān)系

2.3 米山水庫電導(dǎo)率模擬計算

通過上述分析可以看出，米山水庫電導(dǎo)率在2016年前處于自然狀態(tài)，因此可以利用1987—2015年電導(dǎo)率時間序列通過GM(1,1)分解模型建立自然狀態(tài)電導(dǎo)率模型；運用該模型外推，獲得2016—2019年自然狀態(tài)下的電導(dǎo)率，利用同時間序列實測電導(dǎo)率減去自然狀態(tài)電導(dǎo)率，獲得區(qū)外調(diào)水對電導(dǎo)率增加的貢獻(xiàn)率，從而分別獲得氣候變化和區(qū)外調(diào)水對電導(dǎo)率增加的影響。

2.3.1 電導(dǎo)率的趨勢項

利用1987—2015年自然狀態(tài)下電導(dǎo)率值作為原始數(shù)據(jù)，建立GM(1,1)模型并求得模擬模型為

(10)

2.3.2 電導(dǎo)率的周期波動項

GM(1,1)模型的殘差值可以由式(5)求得，如圖7所示。利用式(6)帶截距的正弦曲線對殘差值進行擬合，得到殘差周期擬合曲線的函數(shù)表達(dá)式如下：

圖6 灰色系統(tǒng)GM(1,1)模型對米山水庫電導(dǎo)率模擬值與實測值

(11)

由圖可以看出，所求得的正弦殘差周期擬合曲線與殘差值整體擬合較好，擬合精度較高，能夠反映殘差的周期性變化。

2.3.3 電導(dǎo)率的隨機波動項

根據(jù)式(9)，合并電導(dǎo)率趨勢項、周期波動項和隨機波動變化項，可以對自然狀態(tài)下米山水庫電導(dǎo)率動態(tài)變化進行模擬預(yù)測，結(jié)果見圖8。由圖可以看出：米山水庫自然狀態(tài)下電導(dǎo)率值從1987年開始不斷增大，到2002年達(dá)到最高值；2002—2015年電導(dǎo)率值呈現(xiàn)減小的趨勢，從2016年開始電導(dǎo)率值又有增大的趨勢。同時，電導(dǎo)率模擬值與實際觀測值整體擬合較好，絕大部分觀測值均在模擬值的上、下限范圍內(nèi)。此外，模擬值與實測值一樣，都體現(xiàn)了明顯的上升和回落的周期性變化，說明GM(1,1)分解模型的模擬精度較高。

圖7 米山水庫電導(dǎo)率殘差周期擬合曲線、隨機波動曲線及殘差值

2.3.4 氣候變化和區(qū)外調(diào)水貢獻(xiàn)度定量評價

利用米山水庫2016—2019年電導(dǎo)率實測值減去自然狀態(tài)下的電導(dǎo)率值，即可得到區(qū)外調(diào)水造成的電導(dǎo)率變化值，結(jié)果見表2。從表中數(shù)據(jù)可以看出，2016、2017、2018、2019年自然狀態(tài)下米山水庫的電導(dǎo)率分別為269.3～317.26、270.94～318.9、273.4～321.36、276.63～324.59 μS/cm，處于正常的波動范圍內(nèi)；區(qū)外調(diào)水造成的電導(dǎo)率的增加量分別為540.41～588.37、643.70～691.66、656.18～704.14、629.74～677.70 μS/cm，分別是自然狀態(tài)下電導(dǎo)率的1.85～2.01、2.17～2.38、2.19～2.40、2.09～2.28倍，說明區(qū)外調(diào)水對米山水庫電導(dǎo)率的增加起著主要控制作用。通過對造成米山水庫電導(dǎo)率變化的主控因素進行定量研究表明，影響米山水庫水質(zhì)變化的主要原因不僅僅是氣候變化，這一結(jié)論對于米山水庫建立更加全面的水質(zhì)安全風(fēng)險管理體系，實施更加科學(xué)的水量-水質(zhì)聯(lián)合調(diào)控決策，維護供水安全具有重要意義。

圖8 灰色系統(tǒng)GM(1,1)分解模型對米山水庫電導(dǎo)率的模擬值與實測值

表2 氣候變化和區(qū)外調(diào)水對米山水庫電導(dǎo)率的影響

2.3.5 分析討論

本文中利用GM(1,1)分解模型建立了天然狀態(tài)下電導(dǎo)率的模擬預(yù)測模型，分別對氣候變化和區(qū)外調(diào)水對米山水庫電導(dǎo)率的影響進行了定量評價，研究結(jié)果具有重要的參考價值，但仍然存在著以下不足：

1)米山水庫從2016年才開始實施區(qū)外調(diào)水，較短的區(qū)外調(diào)水時間序列不能充分反映區(qū)外調(diào)水影響下電導(dǎo)率的變化。

2)GM(1,1)分解模型雖然具有用法簡單、所需實測數(shù)據(jù)少和比傳統(tǒng)的GM(1, 1)模型模擬精度更高的優(yōu)勢；但湖庫水體中的污染物濃度受到多種因素的影響(如降雨、氣溫和調(diào)水量等)，因此水質(zhì)預(yù)測模型涉及多種因素的非線性關(guān)系。GM(1, 1)分解模型無法充分考慮多種因素的影響，造成一定誤差。

3)在今后的研究中一方面要增加電導(dǎo)率時間序列的長度，另一方面要采用長短期記憶(LSTM)神經(jīng)網(wǎng)絡(luò)等多參數(shù)機器學(xué)習(xí)方法，通過輸入多個影響變量(如降雨量、氣溫和調(diào)水量等)進行水質(zhì)預(yù)測，由于輸入的不同變量之間可以相互驗證、相互作用，因此能夠有效提高預(yù)測精度。

3 結(jié)論

本文中首先利用STL法對米山水庫降雨量與電導(dǎo)率進行時間序列分解，從宏觀角度揭示它們的趨勢變化、周期波動長度和隨機波動范圍，并提取出電導(dǎo)率的自然狀態(tài)和混合狀態(tài)時間序列；在此基礎(chǔ)上，通過GM(1,1)分解模型建立自然狀態(tài)下電導(dǎo)率模擬模型，并將其外推，定量評價氣候變化和區(qū)外調(diào)水對米山水庫水質(zhì)演變的貢獻(xiàn)，得出以下結(jié)論：

1)通過STL分解結(jié)果可以看出，米山水庫降雨時間序列雖然從2012年開始有減少的趨勢，但整體上仍然呈現(xiàn)穩(wěn)定中略有波動的狀態(tài)；降雨表現(xiàn)出顯著的周期性。電導(dǎo)率在2016年之前同樣呈現(xiàn)穩(wěn)定中略有波動的狀態(tài)，此時電導(dǎo)率主要受氣候變化影響；2016年區(qū)外調(diào)水后電導(dǎo)率快速增大并在2017年達(dá)到最大值，之后處于穩(wěn)定波動狀態(tài)，此時電導(dǎo)率同時受氣候變化和區(qū)外調(diào)水的影響。電導(dǎo)率值的波動同樣呈現(xiàn)明顯的周期性。

2)GM(1,1)分解模型的電導(dǎo)率模擬值與實測值整體擬合較好，絕大部分實測值均在模擬值的上、下限范圍內(nèi)，且模擬值與實測值都能反映明顯的上升和回落的周期性變化，說明GM(1,1)分解模型的模擬精度較高。

3)模擬分析結(jié)果顯示，2016—2019年區(qū)外調(diào)水造成米山水庫電導(dǎo)率的增加量分別為540.41～588.37、643.7～691.66、656.18～704.14、629.74～677.7 μS/cm，分別是自然狀態(tài)下電導(dǎo)率的1.85～2.01、2.17～2.38、2.19～2.40、2.09～2.28倍，說明區(qū)外調(diào)水對米山水庫電導(dǎo)率的增加起到主要控制作用。

4)由于米山水庫水質(zhì)的演變主要受區(qū)外調(diào)水量、客水水質(zhì)影響，因此應(yīng)建立更加全面的區(qū)外調(diào)水水質(zhì)安全風(fēng)險管理體系，實施更加科學(xué)的水量-水質(zhì)聯(lián)合調(diào)控決策，以此維護米山水庫的供水安全和水生態(tài)穩(wěn)定。