R/S和Mann-Kendall法在濟(jì)南市地下水管理模型中的應(yīng)用

齊 歡

(山東省地質(zhì)礦產(chǎn)勘查開發(fā)局八〇一水文地質(zhì)工程地質(zhì)大隊(duì)；山東省地質(zhì)礦產(chǎn)勘查開發(fā)局地下水資源與環(huán)境重點(diǎn)實(shí)驗(yàn)室，濟(jì)南 250014)

0 引 言

濟(jì)南的保泉始于20世紀(jì)80年代初，歷經(jīng)“采外補(bǔ)內(nèi)”保泉、“節(jié)水”保泉、“引黃”保泉、“封井”保泉等階段[1]，雖然目前泉群已連續(xù)噴涌15年，但都未從根本上解決泉水長期連續(xù)壯觀噴涌和城市居民飲用優(yōu)質(zhì)地下水問題。地下水管理模型是利用系統(tǒng)分析原理，為達(dá)到既定管理目標(biāo)所建立的優(yōu)化決策數(shù)學(xué)模型，由地下水流數(shù)值模型和優(yōu)化模型耦合而成[2-4]。地下水管理模型主要計(jì)算各種約束條件下的不同降水保證率的地下水最大可開采量，工作量較大，而準(zhǔn)確預(yù)測未來降水量的趨勢(shì)特征不僅可以減少工作量，而且針對(duì)性更強(qiáng)。

Mann-Kendall法是一種被廣泛用于分析趨勢(shì)變化特征的檢驗(yàn)方法，該方法不僅可以檢驗(yàn)時(shí)間序列趨勢(shì)上升與下降，而且還可以說明趨勢(shì)變化的程度，能很好地描述時(shí)間序列的趨勢(shì)特征[5]?；赗/S分析的Hurst指數(shù)可以根據(jù)時(shí)間序列過去變化趨勢(shì)來預(yù)測將來，該法已被廣泛用于水文，氣象要素時(shí)間序列的定性預(yù)測[6]。本文綜合利用Mann-kendall趨勢(shì)檢驗(yàn)和Hurst指數(shù)分析，預(yù)測未來濟(jì)南市降水量的趨勢(shì)特征，并將結(jié)果作為降水條件應(yīng)用于濟(jì)南市地下水管理模型中，分析不同保泉水位和補(bǔ)源條件下的水源地最優(yōu)開采量，做到有的放矢。

1 研究區(qū)概況

濟(jì)南市地處魯中低山丘陵與魯西北沖積平原交接帶上，南部為泰山山脈，北部為黃河平原，地勢(shì)南高北低，平原向東北緩傾。平均氣溫為12.8 ℃，多年平均降水量為649.15 mm(1956-2017年)，降水具有明顯的季節(jié)性，汛期6-9月份降水量占全年降水量的70%。地表水主要包括大汶河水系、小清河水系和徒駭河水系。

研究區(qū)包含長孝巖溶水子系統(tǒng)、趵突泉泉域和白泉泉域，面積2 766 km2。東至文祖斷裂；西至黃山巖脈牛角店斷裂，牛角店斷裂；北側(cè)以一系列近東西向隔水?dāng)嗔?，奧陶系灰?guī)r頂板埋深600 m一線為界；南至泰山群變質(zhì)巖地表分水嶺(見圖1)。

圖1 研究區(qū)地質(zhì)構(gòu)造略圖Fig.1 Sketch map of geological structure in the study area

2 研究數(shù)據(jù)和方法

2.1 研究數(shù)據(jù)

選取的濟(jì)南市1956-2017年降水量資料來自中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/site/index.html)，站號(hào)為54823，對(duì)于個(gè)別缺失數(shù)據(jù)，采用插值法確定。在數(shù)據(jù)處理方面，首先將逐日的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，分別得到各月、各年的降水量數(shù)據(jù)。

2.2 研究方法

2.2.1 Mann-kendall趨勢(shì)檢驗(yàn)

Mann-Kendall檢驗(yàn)法是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)法，一般取顯著性水平α=0.05，則臨界值UF0.05=±1.96，即95%的置信區(qū)間在±1.96之間。若UFk>0，則表明序列呈上升趨勢(shì)，UFk<0，則表明序列呈下降趨勢(shì)。當(dāng)它們超過臨界線時(shí)，表明上升或下降趨勢(shì)顯著，超過臨界線的范圍確定為出現(xiàn)突變的時(shí)間區(qū)域。如果UFk和UBk兩條曲線出現(xiàn)交點(diǎn)，且交點(diǎn)在臨界線之間，那么交點(diǎn)對(duì)應(yīng)的時(shí)刻便是置信度的突變點(diǎn)[7-9]。

2.2.2 Hurst指數(shù)分析

R/S分析是Hurst H.E.在大量實(shí)證研究基礎(chǔ)上提出的一種時(shí)間序列統(tǒng)計(jì)方法，該方法屬于非參數(shù)分析法，對(duì)考察的對(duì)象幾乎不作任何假設(shè)。通過Hurst指數(shù)H隨時(shí)間尺度的變化規(guī)律，可預(yù)測時(shí)間序列變化發(fā)展的趨勢(shì)性[10, 11]。當(dāng)H>0.5時(shí)，序列未來趨勢(shì)與過去呈正相關(guān)，或稱持續(xù)性，即過去一段時(shí)間內(nèi)的增長趨勢(shì)意味著未來也有增長趨勢(shì)，反之亦然；且H越接近1，持續(xù)性越強(qiáng)。當(dāng)H<0.5時(shí)，為負(fù)相關(guān)性，或稱反持續(xù)性，即過去一段時(shí)間內(nèi)的增長趨勢(shì)意味著未來減少的趨勢(shì)，反之亦然；H越接近0，反持續(xù)性越強(qiáng)。

2.2.3 Mann-kendall與Hurst方法分析

R/S法可以定性的認(rèn)識(shí)序列過去與未來是否存在相同或相反的變化特征，著重從定性的角度揭示未來的變化特征；Mann-kendall法利用線性回歸得出趨勢(shì)方程，揭示了序列的趨勢(shì)特征，著重從定量的角度分析序列在某一時(shí)間段內(nèi)的趨勢(shì)特征。因此綜合應(yīng)用R/S和Mann-Kendall法可以得出濟(jì)南市降水量未來的趨勢(shì)特征。

3 結(jié)果與分析

3.1 年降水量變化特征

對(duì)濟(jì)南市1956-2017年降水量序列進(jìn)行Mann-kendall趨勢(shì)檢驗(yàn)和Hurst指數(shù)分析，如圖2所示。

圖2 年降水量Mann-kendall趨勢(shì)檢驗(yàn)和Hurst指數(shù)分析圖Fig.2 Mann-kendall trend test and Hurst index chart of annual precipitation

從圖2可以看出：濟(jì)南降水序列UF的波動(dòng)情況較復(fù)雜，1956-1960年UF小于0，處于降水較少時(shí)期，UF與UB線在1961年出現(xiàn)交點(diǎn)，且位于置信區(qū)間內(nèi)，降水發(fā)生突變，降水量由1960年的460.5 mm猛增至1961年的983.6 mm。1967至1993年，僅1973-1980年UF大于0，處于降水較多時(shí)期，其余19年UF均小于0，處于降水較少時(shí)期。1994年以后，UF均大于0，且數(shù)值不斷增大，該時(shí)期降水量呈增加的趨勢(shì)；UF與UB線在1990年出現(xiàn)交點(diǎn)，且位于置信區(qū)間內(nèi)，降水發(fā)生突變，降水量由1989年的365.0 mm增至1990年的779.7 mm，變化顯著。濟(jì)南市降水的H值為0.67，具有很好的正持續(xù)性，未來降水趨勢(shì)與過去一致。綜上可知，未來濟(jì)南市處于降水較多時(shí)期，以豐水年為主。

3.2 地下水?dāng)?shù)值模型

采用地下水模擬軟件GMS10.0建立覆蓋研究區(qū)的地下水流數(shù)值模型。模型分為3層，第一層為潛水含水層，該目的層主要為第四系全新統(tǒng)及上更新統(tǒng)地層，概化為二維流；第二層為越流層，為潛水含水層底板以下的黏土層，石炭系、二迭系等弱透水層，只考慮垂向一維流；第三層為承壓水含水層，該目的層是寒武系和奧陶系灰?guī)r地層[12]。模擬期為2012年10月至2014年10月，采用200×200 m網(wǎng)格剖分，通過模型識(shí)別、驗(yàn)證過程，使模擬水位可以較好的擬合實(shí)測水位(見圖3)，水文地質(zhì)參數(shù)可以較為真實(shí)地反映研究區(qū)的實(shí)際水文地質(zhì)條件。

3.3 地下水管理模型

3.3.1 總體方案

在地下水?dāng)?shù)值模型的基礎(chǔ)上，利用脈沖響應(yīng)矩陣法建立地下水管理模型。在保持泉水持續(xù)噴涌的前提下，最大限度的開采地下水資源量，保泉水位分別取27.60和28.00，其余地段水位不低于有觀測資料以來的最低水位。

根據(jù)Mann-kendall和Hurst方法對(duì)濟(jì)南市降水量趨勢(shì)特征的分析，設(shè)置水文年為豐水年(P=25%)，降水量為764.4 mm。

濟(jì)南市已開展玉符河、興濟(jì)河、歷陽湖生態(tài)補(bǔ)源，為更好地與實(shí)際情況相結(jié)合，做到補(bǔ)源的常態(tài)化，本次共設(shè)計(jì)了3種管理方案，分別為：①方案1：無生態(tài)補(bǔ)源。②方案2：全年生態(tài)補(bǔ)源，入滲量為：玉符河10 萬m3/d，興濟(jì)河1 萬m3/d，歷陽湖1 萬m3/d。③方案3：1-6月生態(tài)補(bǔ)源，入滲量為：玉符河10 萬m3/d，興濟(jì)河1 萬m3/d，歷陽湖1 萬m3/d。

圖3 典型水位觀測孔水位擬合圖Fig.3 Fitting chart of typical water level observation well

3.3.2 管理時(shí)段劃分

研究區(qū)降水具有明顯的季節(jié)性，地下水位受到降水的影響，表現(xiàn)為汛期水位升高，枯水季節(jié)水位緩慢下降。一般6月份水位最低，10月水位最高。根據(jù)以上特點(diǎn)，將一個(gè)水文年劃分出三個(gè)管理時(shí)段，即10、11、12、次年1月為第一管理時(shí)段，時(shí)間為123 d；2、3、4、5月為第二管理時(shí)段，時(shí)間為120 d；6、7、8、9為第三管理時(shí)段，時(shí)間為122 d[13, 14]。

3.3.3 約束條件

地下水開采量數(shù)據(jù)主要來自濟(jì)南市城鄉(xiāng)水務(wù)局，部分?jǐn)?shù)據(jù)通過實(shí)際調(diào)查取得。農(nóng)業(yè)全年開采的巖溶地下水折合為8～10 萬m3/d，自備井開采的巖溶水量15 萬m3/d左右(包含張馬屯鐵礦排水4 萬m3/d、黃臺(tái)電廠自備井2.9 萬m3/d、長清水源地4 萬m3/d、齊河水源地4 萬m3/d等)，共計(jì)約24 萬m3/d?？紤]到白泉泉域未開展河流補(bǔ)源，故白泉泉域內(nèi)水源地參考最近三年的實(shí)際開采量，作為開采常量來處理，其中宿家水源地和李莊水源地開采量為0，白泉水源地開采量為1.1 萬m3/d，武家水源地開采量1.8 萬m3/d，黃土崖水源地開采量10 萬m3/d。第四系沖洪積扇孔隙水，主要用于農(nóng)業(yè)灌溉，開采點(diǎn)分散，因此，該開采量不作為決策變量，在管理方案中，只對(duì)孔隙水水位施加約束，以保證農(nóng)業(yè)供水。

水源地開采量不超過現(xiàn)有開采井的供水能力，即上限為各水源地的最大開采能力，下限根據(jù)濟(jì)南市區(qū)供水管網(wǎng)能力確定其供水量的約束值，綜上，確定8個(gè)水源地為優(yōu)化管理對(duì)象(見圖4和表1)。

圖4 主要水源地及河流分布圖Fig.4 Distribution map of main water sources and rivers

開采量/ (萬m3·d-1)曹樓橋子里冷莊古城峨眉山大楊莊臘山東源上限5.008.004.008.008.0010.005.005.00 下限00000003.00

3.3.3 結(jié)果分析

根據(jù)地下水管理模型運(yùn)行結(jié)果，得到豐水年各優(yōu)化管理水源地的最大開采量，如表2和表3所示。

從模擬可以看出：方案1時(shí)，在保泉水位為27.6 m時(shí)，各水源地的最優(yōu)開采量為平均8.37 萬m3/d；保泉水位為28.0時(shí)，各水源地的最優(yōu)開采量為平均5.43 萬m3/d，各水源地的最優(yōu)開采量隨著生態(tài)補(bǔ)源量的增加而增大，生態(tài)補(bǔ)源對(duì)濟(jì)南市開采巖溶地下水起到積極的作用。保泉水位取27.6 m時(shí)，相對(duì)于方案1，方案3時(shí)水源地的平均開采量增加了40.22%，方案2增加了45.68%；保泉水位取28.0 m時(shí)，相對(duì)于方案1，方案3時(shí)水源地的平均開采量增加了88.04%，方案2增加了96.50%。從“保泉供水”的角度，全年常態(tài)化補(bǔ)源效果更好，但從經(jīng)濟(jì)方面考慮，生態(tài)補(bǔ)源存在用水量大，費(fèi)用高的問題。同方案1相比，方案2比方案3補(bǔ)源入滲量增加一倍，而增采量僅增加了5.46%～8.47%?？紤]到河流的實(shí)際入滲能力，應(yīng)盡量不增采巖溶地下水，同時(shí)可在豐水年豐水期，利用水庫放水對(duì)河流的補(bǔ)給，科學(xué)合理的適當(dāng)開采巖溶地下水，滿足濟(jì)南市人民飲用優(yōu)質(zhì)地下水的需求。

表2 各水源地優(yōu)化開采量(約束水位27.6 m) 萬m3/d

表3 各水源地優(yōu)化開采量(約束水位28.0 m) 萬m3/d

選擇2013年數(shù)據(jù)對(duì)優(yōu)化結(jié)果進(jìn)行檢驗(yàn)，該年濟(jì)南市降水量723.90 mm，略小于豐水年的764.4 mm。2010-2012年降水量分別為762.20、605.30和537.50 mm，平均635.00 mm，接近濟(jì)南市多年平均降雨量值，可以避免前一年地下水位過高或過低對(duì)模擬結(jié)果的影響。將不同保泉水位、不同管理方案的水源地優(yōu)化開采量代入地下水?dāng)?shù)值模型，以趵突泉水位為例，通過對(duì)比模擬水位和管理時(shí)段節(jié)點(diǎn)最低水位約束值來驗(yàn)證管理模型的可靠性。

從表4可以看出，保泉水位取27.6 m時(shí)，水源地開采量主要受10月1日最低約束水位的控制，水源地開采主要在1月份和6-9月豐水期進(jìn)行；保泉水位取28.0 m時(shí)，水源地開采量受2月1日、6月1日和10月1日最低約束水位的綜合控制，方案1時(shí)2月1日水位略高于最低約束值，考慮到泉水位對(duì)補(bǔ)源響應(yīng)的滯后性，為維持6月1日28.0 m的保泉水位，水源地開采主要在6-9月豐水期進(jìn)行。

表4 水源地優(yōu)化開采時(shí)各個(gè)管理時(shí)段節(jié)點(diǎn)趵突泉水位值Tab.4 Water level value of the Baotuquan spring in each management period when the water source is optimized for mining

保泉水位27.6 m時(shí)，方案3在10月1日的模擬水位低于水位最低約束值0.05 m；保泉水位28.0 m時(shí)，方案2和方案3在10月1日的模擬水位均低于最低約束值0.01 m。該誤差主要由于2013年實(shí)際降水量小于豐水年的764. 4mm，使得降水入滲補(bǔ)給量相對(duì)偏小，從而模擬水位相對(duì)偏低。同時(shí)由于降水量存在年際與年內(nèi)分配不勻的現(xiàn)象，也會(huì)對(duì)地下水管理模型的結(jié)果造成一定影響。綜上可知，濟(jì)南市地下水管理模型的整體模擬效果較好。

5 結(jié) 論

(1)綜合應(yīng)用Mann-Kendall和R/S法對(duì)濟(jì)南市降水量趨勢(shì)特征進(jìn)行預(yù)測。Mann-kendall趨勢(shì)檢驗(yàn)得出，1994年以后，UF均大于0，且數(shù)值不斷增大，表明濟(jì)南市降水量為增加的趨勢(shì)；Hurst指數(shù)分析的H值為0.67，說明濟(jì)南市降水量具有很好的正持續(xù)性，未來降水趨勢(shì)與過去一致。綜上可知，未來濟(jì)南市處于降水較多時(shí)期，以豐水年為主。

(2)建立地下水管理模型，模擬豐水年保泉水位為27.60 m和28.00 m時(shí)，不同補(bǔ)源條件下水源地的最大可開采量。結(jié)果表明，生態(tài)補(bǔ)源對(duì)濟(jì)南市保泉起到積極的作用，水源地的最優(yōu)開采量隨著生態(tài)補(bǔ)源量的增加而增大，但從經(jīng)濟(jì)和地下水轉(zhuǎn)換利用率的方面考慮，在現(xiàn)有開采條件下，應(yīng)該盡量不增采巖溶地下水。

(3)Mann-Kendall和R/S法可以定性定量的預(yù)測濟(jì)南市未來降水量，使地下水管理模型目的性更強(qiáng)，顯著減少重復(fù)工作量。本文主要從年際變化上對(duì)降水量進(jìn)行評(píng)估，而年內(nèi)降水量分配不均也會(huì)對(duì)計(jì)算結(jié)果造成影響，后期可以嘗試以季降水量、月降水量為單位，揭示未來降水量的變化規(guī)律和趨勢(shì)，使管理模型評(píng)估效果更加精細(xì)化。對(duì)于趵突泉各管理時(shí)段約束水位的取得主要基于2003年泉水復(fù)噴以來的水位日觀測值，序列相對(duì)較短，隨著監(jiān)測時(shí)間的延長和數(shù)據(jù)的不斷更新，關(guān)鍵時(shí)間節(jié)點(diǎn)水位約束值的確定也將更加準(zhǔn)確。

□