基于時(shí)間序列模型的濟(jì)南趵突泉地下水位預(yù)測(cè)

張鄭賢，劉 藝，張鋒賢

（1.濟(jì)南大學(xué) 水利與環(huán)境學(xué)院/山東省地下水?dāng)?shù)值模擬與污染控制工程技術(shù)研究中心，山東 濟(jì)南 250022；2.中國(guó)科學(xué)技術(shù)大學(xué) 工程科學(xué)學(xué)院，安徽 合肥 230026）

1 研究背景

濟(jì)南市素有“泉城”的美譽(yù)，轄區(qū)內(nèi)分布多處名泉，其中尤以趵突泉最為著名。而今由于社會(huì)經(jīng)濟(jì)和城市化進(jìn)程的不斷推進(jìn)，使得泉水補(bǔ)給區(qū)面積驟減，大量新增的硬化面積，削弱了泉域巖溶水的補(bǔ)給量，破壞了區(qū)域地下水系統(tǒng)的自然平衡狀態(tài)，近年來(lái)，泉水分布區(qū)地下水位的持續(xù)下降，使得泉水噴涌受到嚴(yán)重威脅［1-2］。為了保障泉水的持續(xù)噴涌，緩解日益突出的水資源供需矛盾，濟(jì)南市在啟動(dòng)引黃供水工程的同時(shí)充分利用南水北調(diào)東線的水源，將外調(diào)水與當(dāng)?shù)氐乇硭嘟Y(jié)合，在泉域巖溶水的直接補(bǔ)給區(qū)開(kāi)展回灌補(bǔ)源工程。

隨著濟(jì)南市保泉供水工作陸續(xù)開(kāi)展，王慶兵等［3］設(shè)計(jì)了4類地下水資源開(kāi)采方案，運(yùn)用地下水流數(shù)值模型模擬了4種方案對(duì)泉群地下水水位及泉水流量動(dòng)態(tài)變化的影響，最終確定泉域地下水可持續(xù)開(kāi)采方案。高宗軍等［4］采用微量元素水文地球化學(xué)方法，對(duì)濟(jì)南泉域范圍內(nèi)的豐枯水期巖溶水化學(xué)特征及其變化進(jìn)行調(diào)查研究。邢立亭等［5］采用示蹤試驗(yàn)、泉水位動(dòng)態(tài)觀測(cè)等方法揭示濟(jì)南巖溶含水介質(zhì)特征。張展羽等［6］針對(duì)地下水位在時(shí)間序列上表現(xiàn)出高度的隨機(jī)性和滯后性，建立了基于主成分分析與多變量時(shí)間序列CAR模型，對(duì)泉域內(nèi)巖溶水水位進(jìn)行預(yù)測(cè)。齊歡等［7］利用地下水?dāng)?shù)值模擬GMS軟件對(duì)玉符河補(bǔ)源進(jìn)行評(píng)價(jià)，提出西郊水位的抬升可以緩解地下水開(kāi)采對(duì)趵突泉水位的影響，減小保泉的壓力。

本文借鑒已有的研究成果，從保泉供水、區(qū)域地下水資源保護(hù)的角度出發(fā)［8］，結(jié)合濟(jì)南泉域巖溶水補(bǔ)給區(qū)的水文地質(zhì)特征和現(xiàn)有的回灌補(bǔ)源工程，利用灰色動(dòng)態(tài)GM（1，1）和時(shí)間序列ARIMA、Holt-Winters三種模型分別對(duì)趵突泉?dú)v年的水位數(shù)據(jù)進(jìn)行評(píng)價(jià)分析，并選擇擬合效果最優(yōu)的模型對(duì)泉水位的波動(dòng)趨勢(shì)進(jìn)行預(yù)測(cè)，為濟(jì)南市水生態(tài)文明建設(shè)及泉水資源保護(hù)提供依據(jù)。

2 研究區(qū)概況

濟(jì)南市位于山東省中部，面積8154km2，屬暖溫帶季風(fēng)氣候，多年平均降雨量為647.9mm，在6—9月集中降雨，占全年降雨量的70%，大氣降水是濟(jì)南巖溶地下水最主要的補(bǔ)給來(lái)源，降雨量的多少直接決定區(qū)域內(nèi)巖溶水資源的總量［9］。巖溶泉域位于濟(jì)南中部，為一復(fù)雜的巖溶水系統(tǒng)，對(duì)濟(jì)南泉域的劃分尚存在爭(zhēng)議，目前普遍認(rèn)同的觀點(diǎn)為：南邊界為玉符河、北大沙河源頭的地表水分水嶺（即長(zhǎng)城嶺），北邊界為濟(jì)南北部的石炭、二疊系地層和巖漿巖體，東西兩側(cè)分別以濟(jì)南東部的東塢斷層和長(zhǎng)清附近的馬山斷層作為弱透水邊界與外界發(fā)生水力聯(lián)系，泉域總面積1500km2［10］。趵突泉的泉水主要來(lái)源于南部山區(qū)，大氣降水及地表徑流通過(guò)灰?guī)r出露和裂隙巖溶發(fā)育的地帶入滲形成裂隙巖溶水，該部分水流沿巖層傾斜方向由南向北運(yùn)動(dòng)，至城區(qū)受到火成巖巖體的阻擋而大量匯聚，在灰?guī)r和侵入巖體的接觸地帶及第四系沉積層較薄弱處形成泉水涌出地表［11］。濟(jì)南泉域范圍及趵突泉位置見(jiàn)圖1。

圖1 濟(jì)南泉域范圍圖

3 研究方法

3.1 灰色時(shí)間序列模型簡(jiǎn)介對(duì)灰色系統(tǒng)建立的預(yù)測(cè)模型稱為灰色模型（Grey Model），簡(jiǎn)稱GM模型，通過(guò)少量的、不完備的信息，建立灰色微分預(yù)測(cè)模型，對(duì)事物發(fā)展規(guī)律作出模糊性的長(zhǎng)期描述［12］。時(shí)間序列分析是根據(jù)系統(tǒng)觀測(cè)值得到的時(shí)間序列數(shù)據(jù)，通過(guò)擬合逼近和參數(shù)估計(jì)來(lái)建立數(shù)學(xué)模型對(duì)觀測(cè)序列進(jìn)行評(píng)價(jià)預(yù)測(cè)［13］。地下水位波動(dòng)成因復(fù)雜，是典型的灰色系統(tǒng)，本文選取了濟(jì)南泉域趵突泉近5年的地下水位進(jìn)行分析，并采用三種灰色時(shí)間序列模型對(duì)趵突泉的水位進(jìn)行擬合，通過(guò)比較三種模型的擬合結(jié)果，選擇最優(yōu)模型對(duì)濟(jì)南泉域的地下水位進(jìn)行預(yù)測(cè)、評(píng)價(jià)。3種模型適用條件及優(yōu)缺點(diǎn)對(duì)比見(jiàn)表1。

表1 3種灰色時(shí)間序列模型適用條件及優(yōu)缺點(diǎn)對(duì)比［14］

3.2 模型分類（1）GM（1，1）模型。已知觀測(cè)序列經(jīng)過(guò)一次累加生成序列（1-AGO），建立灰色微分方程：

記則由最小二乘法使得J(u)=(Y-Bu)T(Y-Bu)取最小值，u的估計(jì)值為：

預(yù)測(cè)方程為：

式中：x(0)(k)為觀測(cè)值序列，x(1)(k)為一次累加序列，u為參數(shù)向量，Y和B為參數(shù)矩陣，a為發(fā)展系數(shù)，b為灰色向量，t和k為時(shí)間響應(yīng)值，u*為u的估計(jì)值，x*(1)(k+1)為預(yù)測(cè)值。

（2）ARIMA模型。模型表達(dá)式為：

其中等式左邊的模型是自回歸AR（p）部分，非負(fù)數(shù)p表示自回歸階數(shù)，{φ1，φ2，…，φp}稱為自回歸系數(shù)；等式右邊是移動(dòng)平均MA（q）模型，非負(fù)數(shù)q稱為移動(dòng)平均階數(shù)，{θ1，θ2，…，θp}稱為移動(dòng)平均系數(shù)，εt為隨機(jī)干擾誤差項(xiàng)，將觀測(cè)序列進(jìn)行d次差分，作平穩(wěn)化處理。

（3）Holt-Winters模型

預(yù)測(cè)方程為：

式中：St為平滑值即水平分量，α為水平權(quán)重；b為長(zhǎng)期趨勢(shì)值，γ為趨勢(shì)權(quán)重；It為季節(jié)分量，β為季節(jié)權(quán)重；L為季節(jié)長(zhǎng)度（每年的月數(shù)或季數(shù)）；t為當(dāng)前時(shí)間；xt為是實(shí)際觀測(cè)值；yt+k為預(yù)測(cè)值，k為預(yù)測(cè)超前期數(shù)。

4 趵突泉地下水位預(yù)測(cè)

4.1 趵突泉水位波動(dòng)分析對(duì)2012年5月至2017年12月濟(jì)南泉域趵突泉的地下水位數(shù)據(jù)進(jìn)行分析，水位變化的時(shí)序圖見(jiàn)圖2，2010—2017年濟(jì)南市降雨量與地下水開(kāi)采量見(jiàn)表2。

圖2 趵突泉地下水位變化時(shí)序圖

由圖2可以看出趵突泉的水位整體呈現(xiàn)下降趨勢(shì)，在每一年中水位隨著季節(jié)的變化出現(xiàn)波動(dòng)并且表現(xiàn)出一定的規(guī)律。每年3月地下水位開(kāi)始下降，5—6月份下降幅度最為迅猛，主要原因是：3月份入春之后，天氣干燥，氣溫回升塊，降水少且蒸發(fā)量大，土壤失水量大。5—6月份隨著濟(jì)南市西郊春耕的陸續(xù)開(kāi)始，大量抽取地下水用于農(nóng)業(yè)灌溉，導(dǎo)致地下水位迅速下降，7月份之后由于春耕逐漸結(jié)束，大規(guī)模用水接近尾聲，趵突泉水位止跌回漲，7—9月份雨水較為豐富，降雨形成的地表徑流入滲補(bǔ)給深層巖溶水使得泉域地下水位逐漸抬升，9月地下水位波動(dòng)達(dá)到峰值。濟(jì)南位于魯中地區(qū)，冬小麥最佳的播種期為10—11月，因此，10月份之后由于降水減少加之冬灌對(duì)地下水的開(kāi)采使得水位出現(xiàn)回落［15］。全年來(lái)看，影響地下水位的主要因素是季節(jié)性開(kāi)采和降雨，每年4—7月的春耕季節(jié)性開(kāi)采，使得全年中這一時(shí)期的水位最低，該時(shí)期的泉水位是保泉供水的重點(diǎn)監(jiān)控水位［16］。在降雨量較小且泉域地下水開(kāi)采量較大的情況下，若不及時(shí)開(kāi)展回灌補(bǔ)源工程，濟(jì)南市的四大泉群將面臨停噴危機(jī)。

雙頂徑是指胎頭兩個(gè)頂骨隆突間的距離。位置約為兩側(cè)太陽(yáng)穴上方最寬處骨頭之間的長(zhǎng)度。隨著孕周的增加，胎兒雙頂徑逐漸增大，孕足月時(shí)平均值約為9.3厘米。在分娩的過(guò)程中胎兒的頭部通過(guò)產(chǎn)道的最大徑線就是雙頂徑。

分析表2可知，濟(jì)南市近4年中：僅2016年為平水年，連續(xù)3年遭遇枯水年，降雨量整體偏少使得巖溶水的補(bǔ)給量驟減。另一方面，目前濟(jì)南市的地下水開(kāi)采布局，主要還是集中開(kāi)采泉域內(nèi)的東部地下水，濟(jì)西巖溶水開(kāi)采量較小，據(jù)濟(jì)南市名泉辦統(tǒng)計(jì)，目前濟(jì)東地下水集中開(kāi)采區(qū)，日開(kāi)采量22萬(wàn)噸。同期地下水開(kāi)采量雖有降低但降幅較小，因此，在現(xiàn)狀開(kāi)采條件下，趵突泉水位呈現(xiàn)出下降趨勢(shì)。

4.2 數(shù)據(jù)的預(yù)處理及模型訓(xùn)練灰色時(shí)間序列數(shù)據(jù)與普通數(shù)據(jù)不同，它要求建模數(shù)據(jù)必須為光滑離散數(shù)據(jù)并且具有嚴(yán)格的順序，需要定義灰色時(shí)間變量以滿足模型要求，分別對(duì)觀測(cè)數(shù)據(jù)進(jìn)行1次累加和1次差分，使數(shù)據(jù)平滑化。GM（n，1）模型通過(guò)多次擬合及后驗(yàn)差檢驗(yàn)，最終選擇1階線性動(dòng)態(tài)GM（1，1）作為預(yù)測(cè)模型；對(duì)ARIMA模型進(jìn)行反復(fù)訓(xùn)練，其中自回歸階數(shù)p為1，移動(dòng)平均階數(shù)q為1，并進(jìn)行1次季節(jié)差分時(shí)，模型ARIMA（1，0，1）（1，1，1）對(duì)觀測(cè)值的擬合效果較好。Holt-Winters模型經(jīng)測(cè)試最終選擇Holt-Winters加法模型進(jìn)行優(yōu)化擬合時(shí)模擬精度較高。

表2 2010—2017年濟(jì)南市泉域降雨量與地下水開(kāi)采量

4.3 數(shù)學(xué)模型的比選借助MATLAB數(shù)據(jù)處理軟件，采用3種灰色時(shí)間序列模型，對(duì)趵突泉2012—2017年逐月地下水位數(shù)據(jù)進(jìn)行分析，并選取3項(xiàng)指標(biāo)對(duì)模型的擬合結(jié)果進(jìn)行檢驗(yàn)：利用確定系數(shù)（R2）檢驗(yàn)擬合值與觀測(cè)值的相關(guān)程度；均方根誤差（RMSE）檢驗(yàn)擬合數(shù)據(jù)的可靠性；效率系數(shù)（NSE）對(duì)水文模型的預(yù)測(cè)能力進(jìn)行評(píng)價(jià)。觀測(cè)值與擬合值時(shí)序圖見(jiàn)圖3—圖5，模型的檢驗(yàn)見(jiàn)表3。

由圖3—圖5可見(jiàn)，ARIMA與Holt-Winters模型與實(shí)際觀測(cè)序列的變化趨勢(shì)高度相近且模型較為成功的預(yù)測(cè)了水位波動(dòng)的峰值和谷值，可以判定使用該模型是較為合理的；而灰色動(dòng)態(tài)GM（1，1）模型僅成功的捕捉到了水位下降的趨勢(shì)，對(duì)地下水位的周期性、季節(jié)性波動(dòng)變化趨勢(shì)的擬合并不理想。

圖3 GM（1，1）模型擬合值與觀測(cè)值時(shí)序圖

圖4 ARIMA模型擬合值與觀測(cè)值時(shí)序圖

圖5 Holt-Winters模型擬合值與觀測(cè)值時(shí)序圖

表3 模型擬合情況檢驗(yàn)表

表3給出了三種模型的擬合、預(yù)測(cè)情況檢驗(yàn)表，可以看出GM（1，1）模型的擬合、預(yù)測(cè)情況均最差，決定系數(shù)R2和效率系數(shù)NSE數(shù)值均為最低，并且誤差項(xiàng)RMSE數(shù)值較大，說(shuō)明利用該模型分析趵突泉水位波動(dòng)規(guī)律精度較差，因此，不考慮使用GM（1，1）模型。時(shí)間序列ARIMA與Holt-Winters模型的擬合效果整體較好，3項(xiàng)檢驗(yàn)指標(biāo)數(shù)值較為接近，差異不明顯，需要對(duì)擬合情況進(jìn)行深入分析。綜合以上的分析結(jié)果，對(duì)趵突泉的地下水位分析、預(yù)測(cè)優(yōu)先使用ARIMA與Holt-Winters模型。

表4 ARIMA與Holt-Winters模型的擬合優(yōu)度指標(biāo)

為了探尋對(duì)趵突泉水位波動(dòng)規(guī)律擬合效果最佳的模型，分別對(duì)ARIMA與Holt-Winters模型選取了8項(xiàng)擬合優(yōu)度指標(biāo)進(jìn)行對(duì)比分析，由表4可知：兩個(gè)模型均無(wú)離群值（逸出值）且平穩(wěn)的R2數(shù)值較大，說(shuō)明建模數(shù)據(jù)為平滑的離散序列，數(shù)值較為可靠；Ljung-Box Q檢驗(yàn)結(jié)果顯示系統(tǒng)殘差為高斯白噪聲，表明模型的設(shè)定較為合理；為防止模型復(fù)雜度較高引起過(guò)擬合現(xiàn)象，根據(jù)“BIC貝葉斯信息量準(zhǔn)則”選擇最優(yōu)模型（BIC的值越小越好），結(jié)果表明：Holt-Winters模型對(duì)水位觀測(cè)值的擬合效果優(yōu)于ARIMA模型；此外，Holt-Winters模型的MAPE（平均相對(duì)誤差絕對(duì)值）、MaxAPE（最大平均相對(duì)誤差絕對(duì)值）、MAE（平均絕對(duì)誤差）、MaxAE（最大平均絕對(duì)誤差）評(píng)價(jià)指標(biāo)均為最小，說(shuō)明Holt-Winters相對(duì)于ARIMA模型的誤差項(xiàng)更小，對(duì)原始數(shù)據(jù)的擬合精度更高。因此，綜合考慮模型的合理性、復(fù)雜度及誤差項(xiàng)等優(yōu)度指標(biāo)，最終選擇Holt-Winters模型作為預(yù)測(cè)模型。

4.4 最優(yōu)模型的應(yīng)用通過(guò)對(duì)3種數(shù)學(xué)模型的深入比選，最終選擇擬合效果最優(yōu)的Holt-Winters模型對(duì)趵突泉的地下水位進(jìn)行分析，并利用該模型預(yù)測(cè)了2018—2019年趵突泉的水位波動(dòng)情況。Holt-Winters時(shí)間序列模型參數(shù)見(jiàn)表5，2018—2019年趵突泉水位預(yù)測(cè)值時(shí)序圖見(jiàn)圖6。

表5 Holt-Winters模型的參數(shù)估計(jì)表

圖6 觀測(cè)值與擬合值、預(yù)測(cè)值時(shí)序圖

隨著濟(jì)南市保泉工作的陸續(xù)開(kāi)展，地下水開(kāi)采量近年來(lái)已趨于穩(wěn)定，由表2可知濟(jì)南市近8年來(lái)降雨量持續(xù)偏枯，近10年中未出現(xiàn)豐水年，近20年中僅有2年為豐水年［17］。以近年的降雨與開(kāi)采現(xiàn)狀為依據(jù)，利用Holt-Winters模型預(yù)測(cè)短時(shí)期內(nèi)的泉水位的波動(dòng)情況，并以此為基礎(chǔ)評(píng)價(jià)保泉形勢(shì)。

圖6給出了Holt-Winters模型的預(yù)測(cè)值、預(yù)測(cè)值95%的置信區(qū)間值，并標(biāo)注出了水位波動(dòng)的峰值、谷值及保泉紅色警戒線、泉水停噴線。

4.5 預(yù)測(cè)結(jié)果分析利用Holt-Winters模型對(duì)2018年1月至2019年12月趵突泉的地下水位進(jìn)行預(yù)測(cè)，由圖6可知：趵突泉的地下水位在未來(lái)的24個(gè)月總體呈現(xiàn)出下降趨勢(shì)且水位標(biāo)高均值分別為27.734m和27.605m，2019年6月泉水位預(yù)測(cè)值為近7年最低水平，數(shù)值僅為27.124m，低于《濟(jì)南市保持泉水噴涌應(yīng)急預(yù)案》規(guī)定的27.6m的保泉紅色警戒線，逼近27.01m的停噴線。表明現(xiàn)狀降雨與開(kāi)采條件下趵突泉在2019年6月可能存在潛在停噴危機(jī)。由此可見(jiàn)，濟(jì)南市未來(lái)兩年內(nèi)“保泉”形勢(shì)依然嚴(yán)峻。

表6 2013—2019年趵突泉年平均水位及峰值谷值對(duì)比 （單位：m）

由表6可知：2014、2015年濟(jì)南市連續(xù)遭遇枯水年，農(nóng)業(yè)灌溉對(duì)地下水的季節(jié)性開(kāi)采量相對(duì)較大而有效補(bǔ)給量偏少，趵突泉水位下降明顯；近4年中僅2016年為平水年，降雨量總體偏枯。濟(jì)南市政府及時(shí)啟動(dòng)保泉應(yīng)急預(yù)案，加大了外調(diào)水源的使用力度，關(guān)停了泉域范圍內(nèi)的部分自備井，在一定程度上緩解了水位下降勢(shì)頭，使得泉水位出現(xiàn)小幅上漲，但由于降雨量并未明顯增加，深層巖溶水的有效補(bǔ)給量增長(zhǎng)并不顯著，地下水超采形成的局部降落漏斗在短時(shí)期內(nèi)未能得到完全修復(fù)，使得泉水位難以持續(xù)上漲。Holt-Winters模型的預(yù)測(cè)結(jié)果表明：2018和2019年的泉水位預(yù)測(cè)值與2013—2017年同期的多年平均地下水位觀測(cè)值相比分別下降了53.4cm和66.3cm，降幅為1.9%和2.4%，水位下降幅度較大，且仍然存在不斷下降的趨勢(shì)，采用必要的調(diào)蓄補(bǔ)源措施已迫在眉睫。

5 結(jié)論

通過(guò)對(duì)趵突泉2012—2017年逐月的地下水位變化規(guī)律進(jìn)行分析，得出泉水位既存在不斷下降的趨勢(shì)，又受季節(jié)變化的影響呈現(xiàn)出明顯的周期性波動(dòng)，因而，采用灰色時(shí)間序列模型對(duì)趵突泉的地下水位進(jìn)行分析預(yù)測(cè)較為合適。對(duì)3種模型水位觀測(cè)值擬合效果進(jìn)行檢驗(yàn)，初步確定ARIMA與Holt-Winters時(shí)間序列模型的預(yù)測(cè)效果要優(yōu)于灰色動(dòng)態(tài)GM（1，1）模型，再分別選取8種擬合優(yōu)度指標(biāo)對(duì)ARIMA與Holt-Winters模型進(jìn)行深入分析，最終選擇擬合效果最優(yōu)的Holt-Winters模型對(duì)趵突泉的水位進(jìn)行預(yù)測(cè)，該方法的擬合優(yōu)度、可信度及模型的預(yù)測(cè)能力均較好，適用于中、短期正常氣象條件的地下水位預(yù)測(cè)。

采用Holt-Winters模型預(yù)測(cè)了2018—2019年2年內(nèi)趵突泉的地下水位變化情況，預(yù)測(cè)結(jié)果表明2018、2019年的泉水位均值分別值為27.734m和27.605m，泉水位波動(dòng)的峰值為28.215m、谷值為27.124m且谷值將出現(xiàn)于2019年的6月份，是近7年來(lái)的最低水位，該水位低于27.6m的保泉紅色警戒線，逼近27.01m的泉水停噴線，照此趨勢(shì)發(fā)展，若無(wú)法有效緩解水位下降的勢(shì)頭，在2019年6月趵突泉泉水可能面臨潛在的停噴危機(jī)。