郭 藝,王 楓,甘甫平,閆柏琨
1.中國自然資源航空物探遙感中心,北京 100083;2.生態環境部環境規劃院,北京 100012
時間序列分析是一種基于系統的觀點研究各種相依有序的離散數據集合的方法,包括一般統計分析和數學模型的建立與預測[1,2]。其中基于時間序列數據的數學規律建立擬合度較高的數學模型并進行預測具有原理簡單、可操作性強、預測結果合理等優點[3,4],被廣泛應用于水動態特征分析及預測方面。常用的水動態預測模型包括季節分解模型[5,6]、移動平均模型[7]、指數平滑模型[8]、多層階梯模型[9],混合回歸模型[10,11]和小波神經網絡模型[12]等。
巖溶地貌占地球陸地表面的7%~12%,巖溶水提供了全球超過25%的人口的日常用水[13,14]。中國巖溶區面積達344×104km2,巖溶水是重要的供水來源[15]。巖溶泉可視為巖溶水系統的輸出,其水動態變化受多種輸入激勵的影響,是反映巖溶水系統水動力場的一個重要參數[16,17]。因此,客觀的分析巖溶泉水動態特征,并科學的預測泉流量,對巖溶水資源的開發利用與保護至關重要[18,19]。巖溶泉水位、流量表現出隨機性、非線性和多時間尺度變化等復雜特征[20,21]。過去對地下水水位或泉流量預測的研究大多數是利用單個預測模型對單個時間尺度的預測,缺乏不同時間尺度預測結果和不同模型預測結果的對比研究,本文利用移動平均模型和指數平滑模型預測濟南城區巖溶泉日尺度、月尺度和年尺度的流量,為當地巖溶水資源的開發利用和保護提供參考。
濟南泉域位于山東省濟南市,泉域地形南高北低,屬于半干旱大陸性季風氣候,多年平均氣溫為13.8℃,平均降水量為685 mm,主要集中在7—9月份。濟南市主要發育的河流有黃河和小清河兩大水系,與巖溶水系統有水力聯系的地表河流主要為黃河水系的玉符河和北大沙河。
濟南巖溶含水層主要由寒武系和奧陶系的碳酸鹽巖組成。巖溶水的主要補給來源為降水的入滲補給以及河道的滲漏補給,巖溶水系統接受補給后,水流沿著溶孔和裂隙由西南、南部和東南匯向城區泉群,在2.6 km2的城區出露100多個巖溶泉,其中以趵突泉、黑虎泉、五龍潭和珍珠泉最為聞名。
本文所用的2012年5月2日至2019年10月31日的趵突泉水位數據由“濟南市城鄉水務局”網站的“水位”專題資訊獲取,1959年到2002年濟南泉域地下水水位和巖溶泉流量數據由《濟南巖溶水系統研究》[22]中收集。
移動平均模型是基于移動平均法而建立的預測模型,又稱為滑動平均法,其基本思想是根據時間序列資料逐項推移,依次計算包含一定項數的序時平均值,以反映長期趨勢。使用移動平均法可以消除時間序列由于受周期變動和隨機波動的影響,顯示出數據的發展方向與趨勢(即趨勢線),然后依趨勢線分析預測序列的長期趨勢,其表達式為:
(1)

指數平滑模型是移動平均模型的改良模型,其特點是數據的權重按照由遠及近遞增,即對較近的數據給與更大的權重。根據平滑次數的不同,可分為一次指數平滑法、二次指數平滑法和三次指數平滑法。對時間序列數據建立指數平滑模型時,當時間序列無明顯趨勢變化的情況下,一次指數平滑模型即可滿足預測要求,一次指數平滑模型是用第t期的指數平滑值作為預測值進行預測的,其表達式為
(2)
式中:Et為第t期平滑值(E0為初始值,通常設E0=y0);α為平滑系數(0<α<1),通過試算法確定。
對建立的移動平均模型和指數平滑模型,本文采用均方根誤差(RMSE)來判別模型預測的精度,RMSE越小,預測結果的可靠性越大,反之預測結果不大可靠。其計算公式為:
(3)

1959—2008年濟南地下水水位和泉流量(包括西部的長孝巖溶區,中部的濟南城區泉域和東部白泉泉域)的長期監測結果表明泉流量與地下水水位之間存在線性正相關[23]。周娟[24]通過測量2015年2月23日至5月21日期間趵突泉日時間尺度的水位(H)和流量(Q),建立了兩者之間的回歸方程:Q=7.019H-179.45(R2=0.898)。由圖1可知,1959年到2002年城區巖溶泉年時間尺度的流量(Q)與地下水水位(H)之間存在正相關,回歸方程為:Q=7.071H-181.176(R2=0.899)。該方程與周娟觀測所得的方程系數一致,利用該方程將泉水位轉為泉流量。利用上述方程得到的2012年—2019年泉流量變化曲線見圖2。

圖1 1959—2002年地下水水位與泉流量關系
泉群出流和巖溶水人為開采為濟南巖溶水系統的主要排泄途徑。泉流量的年際變化曲線見圖3,濟南泉域1959—2002年平均流量為14.20萬m3/d,最大為1962年的51.15萬m3/d,最小流量出現在2000年,此時泉群全部斷流。由圖3可知,20世紀60年代,濟南市地下水開采量較小,泉群的流量高達30萬~50萬m3/d,20世紀70年伊始,隨著地下水開采量的增加,巖溶泉流量顯著降低,甚至出現斷流的現象。趵突泉于1972年6月出現了歷史上首次斷流現象,之后在1982年、1989年以及2000年又出現了不同程度的斷流,其中最長的斷流時間長達926天(1999年3月2日至2001年9月17日)。
21世紀初,隨著巖溶泉保護力度的加大,泉域巖溶水開采量減小,泉水位基本維持在28 m左右,流量保持在10萬m3/d。圖2顯示,2012年以后泉流量的周期性變化反映了降水量的周期項變化。2012年、2013年和2016年屬于豐水年,泉流量較高,年最高泉流量高于25萬m3/d;2014年、2015年和2017年屬于枯水年,泉流量較低,年最高泉流量低于20萬m3/d;2018年和2019年屬于平水年,年最高泉流量介于20萬~25萬m3/d。
由圖4可知,受降水的影響,2012—2019年泉流量呈現典型的季節性變化,1—5月持續下降,最小泉流量出現在6月,為15.50萬m3/d,從7月份開始流量上升,最高流量出現在9月,為22.23萬m3/d,多年月平均流量為18.63萬m3/d。

圖2 2012—2019年濟南泉域降水量、泉水位與泉流量變化曲線Fig.2 The daily variation of precipitation,spring water level and spring discharge in Jinanspring area from the year of 2012 to 2019
3.3.1 移動平均模型
利用移動平均模型預測泉流量的關鍵是確定移動平均的階數,本文采用試算法確定移動平均的階數。對于日時間尺度的泉流量時間序列,利用式(1)分別求取泉流量時間序列2階、5階、10階、15階、20階、25階和30階移動平均數,并依據式(3)計算其標準誤差。由圖5可知,隨著移動平均階數的增加,預測模型的標準誤差呈線性增大趨勢,故采用2階移動平均模型來對日時間尺度的泉流量進行預測,其預測結果如圖6所示。由圖6可知,基于日時間尺度的2階移動平均模型所得到的泉流量預測值與觀測值變化趨勢一致,即該模型能很好地擬合泉流量。但在利用該模型對泉水位進行預測時,三個時序之后(即當t=n+3時),泉流量的預測值保持在19.05萬m3/d,不再呈現時間變化。

圖3 1959—2002年濟南泉域降水量、開采量與泉流量變化曲線Fig.3 The year variation of precipitation,exploitation and spring discharge in Jinan spring area from the year of 1959 to 2002

圖4 2012—2019年泉流量月平均變化曲線

圖5 移動平均階數、平滑系數與標準誤差的關系曲線Fig.5 The relationship curve of moving average order,smoothing coefficient and standard error

圖6 日時間尺度泉流量的觀測值與預測值Fig.6 The daily prediction value and observation value of spring discharge
將2012年5月2日至2019年10月31日的日時間尺度的泉流量時間序列按月求取平均值得到月時間尺度的泉流量時間序列。對于月時間尺度的泉流量時間序列,分別求取2階至12階移動平均數,并計算其標準誤差,結果也如圖5所示。由圖5可知,隨著移動平均階數的增加,預測模型的標準誤差呈增大趨勢,但誤差增加速度小于階數增加速度,故也采用2階移動平均模型來對月時間尺度的泉流量進行預測,其預測結果如圖7所示。由圖7可知,基于月時間尺度下的2階移動平均模型所得到的泉流量預測值與觀測值變化趨勢一致,即該模型能很好地擬合泉流量。但在利用該模型對泉水位進行預測時,6個時序之后(即當t=n+6時),泉流量的預測值保持在19.88萬m3/d,不再呈現時間變化。

圖7 月時間尺度泉流量的觀測值與預測值Fig.7 the monthly prediction value and observation value of spring discharge
年時間尺度的泉流量時間序列可由上述日時間尺度的泉流量序列按年求取平均值獲得,其中2012年泉流量為5月2日至12月31日的平均值,2019年泉流量為1月1日至10月31日的平均值,其余年份為全年數據平均值。對于年時間尺度的泉流量時間序列,求取泉流量序列2階至4階移動平均數,并計算其標準誤差。由圖5可知,與日時間尺度和月時間尺度不同,隨著模型階數的增加,標準誤差先增加再減小,2階時誤差最小,故本次采用2階移動平均模型來對年時間尺度的泉流量進行預測,其預測結果如圖8所示。由圖8可知,2020年泉流量為17.57萬m3/d。
3.3.2 指數平滑模型
利用指數平滑模型預測泉流量的關鍵是確定平滑系數α,本文采用試算法確定平滑系數。由于平滑系數α的變化范圍為0~1,因此分別取α=0.1,0.2,…,0.9,代入式(2)求取不同時間尺度的泉流量時間序列的預測值,并依據式(3)計算其標準誤差,結果如圖5所示。由圖5可知,不同時間尺度條件下,隨著平滑系數α的增大,預測模型的標準誤差都呈減小趨勢,故采用α=0.9對泉流量進行預測。

圖8 年時間尺度泉流量的觀測值與預測值Fig.8 the yearly prediction value and observationvalue of spring discharge
平滑系數為0.9的指數平滑模型的預測結果如圖6~圖8所示。由圖6可知,基于日時間尺度下的指數平滑模型所得到的泉流量預測值與觀測值變化趨勢一致,即該模型能很好地擬合泉流量,但16個時序后(即當t=n+16時),泉流量的預測值保持在19.03萬m3/天,不再呈現時間變化。由圖7可知,基于月時間尺度下的指數平滑模型所得到的泉流量預測值與觀測值變化趨勢一致,即該模型能很好地擬合泉流量,在利用該模型對泉水位進行預測時,可預測14個時序,14個預測時序的預測值在19.67~20.17萬m3/天之間波動變化。由圖8可知,基于年時間尺度的指數平滑模型所得到的2020年泉流量為17.06萬m3/d。
由上述預測結果可知,利用兩種模型進行預測時,可預測的結果具有時間尺度效應。這種時間尺度效應主要表現在:(1)可預測的時序長度方面;(2)預測值準確度方面;(3)不同時間尺度模型預測值的差異方面。
對移動平均模型來說,日時間尺度模型可預測在3個時序,月時間尺度模型可預測6個時序。對指數平滑模型來說,日時間尺度模型可預測16個時序,月時間尺度模型可預測14個時序。
對移動平均模型和指數平滑模型,日時間尺度、月時間尺度和年時間尺度泉流量的預測值與測量值之間的均方根誤差(RMSE)分別為0.07、0.76和0.99以及0.55、2.74和3.64,表明隨著時間尺度的增加,預測值的誤差增加。
日時間尺度、月時間尺度和年時間尺度的移動平均模型預測的2020年泉流量分別為19.05萬m3/d,19.88萬m3/d和17.57萬m3/d;日時間尺度、月時間尺度和年時間尺度的指數平滑模型預測的2020年泉流量分別為19.03萬m3/d,19.82萬m3/d和17.06萬m3/d。無論是移動平均模型還是指數平滑模型,2020年泉流量月時間尺度的預測值最大,日時間尺度預測值次之,年時間尺度的預測值則顯著小于其他時間尺度的預測值。
造成這種現象的原因主要為無論移動平均模型還是指數平滑模型的原理都是用一組最近的實際數據值來預測未來一期或幾期內序列數值的常用方法,即適用于即期預測,只能預測較少時序的數據,對于多期的數據預測則存在滯后偏差[7]。另一方面,在利用時間序列分析法進行預測時需要大量的過去數據的記錄,本文中由于年泉流量的數據值較少,可能會造成利用年時間尺度進行泉流量預測時產生誤差。
日時間尺度、月時間尺度和年時間尺度的移動平均模型的均方根誤差(RMSE)均小于相同時間尺度的指數平滑模型的均方根誤差,表明對于本文來說基于移動平均模型得到的模擬值更接近觀測值,其預測結果優于指數平滑模型。楊曉俊[25]通過對比上述兩種模型對山西巖溶泉泉流量預測結果發現指數平滑模型的預測效果比移動平均模型好。造成上述這種差異的主要原因在于楊曉俊的時間序列數據選取的是無季節性變化的泉流量數據,而本文中泉流量展現出明顯的季節性變化。因此,即使指數平滑模型是移動平均模型的改良模型,在巖溶泉流量的時間序列存在較明顯的季節性波動時,移動平均模型的預測效果更好。
在對濟南城區巖溶泉水位和泉流量動態特征分析的基礎上,建立了巖溶泉流量移動平均模型和指數平滑模型,并開展了不同時間尺度的泉流量預測,主要結論如下:
(1)濟南城區地下水水位與泉流量之間呈現線性關系,可通過測量泉水位換算泉流量數據。
(2)受人為因素的影響,泉流量1959—2002年呈現出一定的下降趨勢;受自然因素的影響,泉流量2012—2019年呈現出一定的季節性波動;
(3)利用試算法確定出移動平均模型的階數為2,指數平滑模型的平滑系數為0.9,兩種模型的計算結果均較好地擬合了泉流量。
(4)在巖溶泉流量的時間序列存在較明顯的季節性波動時,移動平均模型的預測效果更好。
(5)兩種模型的預測結果顯示了相同的時間尺度效應,即基于日時間尺度的擬合效果最好,月時間尺度次之,年時間尺度誤差最大。