五種地下水埋深預(yù)測模型對比分析
——以肇州縣為例

張嗣路，李治軍，于博文，王濤

（1.黑龍江大學(xué)水利電力學(xué)院，黑龍江哈爾濱 150080；2.黑龍江大學(xué)寒區(qū)地下水研究所，黑龍江哈爾濱 150080）

地下水動態(tài)變化是一個受自然條件與人類活動影響的復(fù)雜的水文過程，在時間序列上具有高度的滯后性與隨機(jī)性。當(dāng)?shù)叵滤裆畛掷m(xù)增加，會造成地面沉降形成地下水漏斗；當(dāng)?shù)叵滤裆钭儨\，會出現(xiàn)土壤次生鹽漬化。因此，預(yù)測地下水埋深是區(qū)域水資源管理配置的重要基礎(chǔ)，對于水資源的可持續(xù)利用和生態(tài)環(huán)境的保護(hù)至關(guān)重要［1］。傳統(tǒng)地下水埋深預(yù)測通過數(shù)學(xué)物理方法建立確定性模型模擬地下水運(yùn)動機(jī)制，其中特定參數(shù)的識別和確定，需要精度較高的數(shù)據(jù)支撐，且計(jì)算量較大，假定條件復(fù)雜［2］。近年來隨著科技發(fā)展和隨機(jī)論理論的成熟，許多學(xué)者將隨機(jī)性模型應(yīng)用到地下水埋深預(yù)測。隨機(jī)性模型基于隨機(jī)理論描述地下水運(yùn)動過程，主要方法有回歸分析法、小波分析法、神經(jīng)網(wǎng)絡(luò)模型、時間序列分析模型等，各種模型應(yīng)用廣泛但仍具有一定局限性［3］。自回歸滑動平均模型等時間序列分析模型預(yù)測精確度較高，但對于非平穩(wěn)趨勢序列和周期性特征處理效果不好［4］；灰色模型所需參數(shù)少、容錯性小，對歷史數(shù)據(jù)有很強(qiáng)的依賴性，誤差偏大；人工神經(jīng)網(wǎng)絡(luò)具有較強(qiáng)的容錯性和非線性映射能力，但存在運(yùn)行時間過長，陷入局部最優(yōu)等問題。近年來，很多專家學(xué)者針對不同地區(qū)地下水埋深運(yùn)用了各類模型進(jìn)行合理預(yù)測模擬，徐強(qiáng)等［6］采用多元線性回歸模型對北京平原平谷地下水水位動態(tài)進(jìn)行預(yù)測；劉雷等［7］采用灰色GM（1，1）理論建立渭北旱塬區(qū)地下水埋深預(yù)測模型；付強(qiáng)等［8］將人工神經(jīng)網(wǎng)絡(luò)進(jìn)行改進(jìn)并應(yīng)用到井灌水稻區(qū)進(jìn)行模擬預(yù)測。文獻(xiàn)研究表明，專家學(xué)者主要運(yùn)用單一預(yù)測模型進(jìn)行地下水埋深預(yù)測研究，因此本文選擇五種常見的預(yù)測模型，對肇州縣地下水埋深進(jìn)行模擬預(yù)測與對比分析，以期了解五種模型在地下水動態(tài)預(yù)測中性能的優(yōu)劣，為地下水動態(tài)監(jiān)測預(yù)報(bào)和地下水資源開采規(guī)劃提供科學(xué)依據(jù)。

1 研究方法

1.1 多元線性回歸模型

多元線性回歸模型是通過建立多元線性回歸方程并進(jìn)行顯著性檢驗(yàn)來研究因變量和多個自變量之間是否存在線性關(guān)系、存在何種數(shù)量變化，根據(jù)多個自變量的變化來預(yù)測模擬因變量值［9］。多元線性回歸模型可以有效地綜合考慮各個自變量的關(guān)系，適用于地下水埋深預(yù)測的多因素性。其數(shù)學(xué)表表達(dá)式為：

式中：β為樣本偏回歸系數(shù)；ε為隨機(jī)誤差項(xiàng)。

1.2 灰色GM（1，1）模型

灰色GM（1，1）模型是基于灰色系統(tǒng)理論對單變量時間序列數(shù)據(jù)進(jìn)行預(yù)測的一種預(yù)測方法。其基本原理是首先對原始數(shù)據(jù)進(jìn)行累加，得到緊鄰均值序列和矩陣，運(yùn)用最小二乘法和求解微分方程對預(yù)測模型進(jìn)行計(jì)算，最后將計(jì)算值進(jìn)行累減得到預(yù)測結(jié)果［10］。具體方法如下：

設(shè)有原始非負(fù)數(shù)據(jù)序列為：

對其進(jìn)行一次累加處理，得到新的數(shù)列：

其中參數(shù)a，b可采用最小二乘法確定，即：

進(jìn)行累減得到灰色GM（1，1）預(yù)測值為：

1.3 基于馬爾科夫鏈優(yōu)化的灰色GM（1，1）模型

馬爾科夫鏈理論基于隨機(jī)過程理論研究系統(tǒng)從一個空間狀態(tài)轉(zhuǎn)移到另一個狀態(tài)的概率進(jìn)而根據(jù)所得到的下一時刻的狀態(tài)來預(yù)測系統(tǒng)未來狀態(tài)的發(fā)展趨勢［11］。通過馬爾科夫鏈將灰色GM（1，1）模型預(yù)測值與實(shí)際值的相對誤差進(jìn)行修正，達(dá)到優(yōu)化灰色GM（1，1）預(yù)測值效果［12］?；隈R爾科夫鏈優(yōu)化步驟如下：

首先對灰色GM（1，1）預(yù)測產(chǎn)生的相對誤差劃分為n個區(qū)間狀態(tài)，建立狀態(tài)轉(zhuǎn)移概率矩陣：

式中：pij表示從狀態(tài)Si轉(zhuǎn)移到狀態(tài)Sj的概率（pij≥0，=1），其中：

式中：Mij表示狀態(tài)Si出現(xiàn)的次數(shù)；Mi表示從狀態(tài)Si轉(zhuǎn)移到狀態(tài)Sj的次數(shù)。

假設(shè)初始狀態(tài)向量S0，在經(jīng)過K步轉(zhuǎn)移后的狀態(tài)向量定義為Sk=S0P(k)。在確定預(yù)測狀態(tài)所處轉(zhuǎn)態(tài)序列區(qū)間后，就能得到對灰色GM（1，1）優(yōu)化修正的預(yù)測結(jié)果［13］，即：

1.4 BP神經(jīng)網(wǎng)絡(luò)模型

BP 神經(jīng)網(wǎng)絡(luò)模型是一種通過誤差逆向傳播進(jìn)行算法訓(xùn)練的多層前饋型人工神經(jīng)網(wǎng)絡(luò)模型［14］。網(wǎng)絡(luò)結(jié)構(gòu)可分為輸入層、隱含層和輸出層，每層包含若干個神經(jīng)元，前后層之間通過網(wǎng)絡(luò)權(quán)值相連接，各層神經(jīng)元之間無連接。選擇Sigmoid 型函數(shù)作為隱含層的激活傳遞函數(shù)，如：

BP 神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)過程可分為信號正向傳播和誤差逆向傳播兩部分，工作信號沿輸入層、隱含層、輸出層的路線逐層計(jì)算傳播，當(dāng)輸出層結(jié)果與目標(biāo)輸出誤差較大時，誤差信號反向傳播對每個層之間的權(quán)值和閾值進(jìn)行修正，最后回到輸入層。兩種信號循環(huán)反復(fù)進(jìn)行訓(xùn)練，直到將誤差控制到期望范圍內(nèi)，學(xué)習(xí)過程結(jié)束。

1.5 基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)模型

遺傳算法按照自然界進(jìn)化理論模擬染色體遺傳、進(jìn)化、變異的過程形成的全局高效優(yōu)化搜索算法。遺傳算法通過對種群進(jìn)行選擇、交叉和變異操作，將適應(yīng)度高的個體保留并進(jìn)行繁衍，產(chǎn)生符合條件的個體即找到最優(yōu)解［15］。遺傳算法具有全局高效搜索能力，可以尋找到BP 神經(jīng)網(wǎng)絡(luò)權(quán)值閾值最優(yōu)解，提高模型的魯棒性，使得最終預(yù)測值更加接近實(shí)測值。

2 實(shí)例應(yīng)用

2.1 研究區(qū)概況

肇州縣位于黑龍江省西南方向，地理位置在東經(jīng)124°48'～125°48'、北緯45°35'～46°16'之間，地形為沖積型平原，地貌較為平坦，占地面積為2 445 km2。肇州縣作為農(nóng)業(yè)大縣，是全國重要的商品糧基地［16］。肇州縣境內(nèi)無江河，地表水資源匱乏，地下水是肇州縣生產(chǎn)生活主要用水來源。受大陸性季風(fēng)氣候影響，肇州縣地下水補(bǔ)給主要來源為降雨入滲和側(cè)向徑流補(bǔ)給，以人工開采、潛水蒸發(fā)和側(cè)向徑流方式排泄。

2.2 數(shù)據(jù)資料

本文選用肇州縣1980-2019 年共40 年的年平均地下水埋深資料。在應(yīng)用多元線性回歸模型、BP神經(jīng)網(wǎng)絡(luò)模型和基于遺傳算法優(yōu)化的BP 神經(jīng)網(wǎng)絡(luò)時，需要明確地下水位埋深變化的主要影響因素［17］。根據(jù)研究區(qū)地下水動態(tài)變化類型并結(jié)合相關(guān)參考文獻(xiàn)，選用降水量、蒸發(fā)量、地下水開采量和前期水位共四個因素作為輸入變量。五種模型都將1980-2009年地下水位埋深作為訓(xùn)練樣本，2010-2019 年地下水位埋深數(shù)據(jù)作為測試集。

2.3 模型建立及訓(xùn)練

2.3.1 多元線性回歸模型

通過SPSS26.0 統(tǒng)計(jì)軟件建立多元回歸線性模型并進(jìn)行顯著性檢驗(yàn)，模型建立過程見表1～3。由表1 可知，復(fù)變相關(guān)系數(shù)R為0.893，決定系數(shù)R2為0.797，調(diào)整后的R2為0.764，表示建立的多元線性回歸模型其所選用的自變量可以解釋76.4%的地下水埋深的變化［18］。表2中，F(xiàn)=24.524＞F0.05，P=0＜0.01，說明建立的回歸模型高度顯著具有良好的統(tǒng)計(jì)學(xué)意義。根據(jù)表3 以P＜0.05 為檢驗(yàn)水平，蒸發(fā)量、降水量和前期水位對于地下水位埋深呈正相關(guān)，差異顯著；地下水開采量差異不顯著［19］。以地水位埋深為因變量，蒸發(fā)量、降水量、前期水位為自變量建立多元線性回歸預(yù)測方程為：

表1 模型擬合匯總Tab.1 Model fitting summary

表2 方差分析Tab.2 Analysis of variance

表3 回歸系數(shù)及顯著性檢驗(yàn)Tab.3 Regression coefficient and significance test

式中：y為研究區(qū)地下水平均埋深；x1為蒸發(fā)量；x3為降水量；x4為前期水位。

2.3.2 灰色GM（1，1）模型

根據(jù)1980-2019 年肇州縣地下水為埋深數(shù)據(jù)和灰色GM（1，1）模型原理，使用Matlab 進(jìn)行編程計(jì)算，得到各參數(shù)值a=-0.010 781 24，u=5.618 878 842，帶入到微分方程中，得到模型預(yù)測公式為：

2.3.3 基于馬爾科夫鏈優(yōu)化的灰色GM（1，1）模型

在基于馬爾科夫鏈對灰色GM（1，1）模型進(jìn)行優(yōu)化時，需要將灰色GM（1，1）模型訓(xùn)練樣本的相對誤差進(jìn)行區(qū)間狀態(tài)劃分［22］。訓(xùn)練的相對誤差范圍為-15.35%～18.64%，分別以-20%、-5%、0、5%、10%、20%為閾值，將其劃分為5 個區(qū)間狀態(tài)。確定狀態(tài)轉(zhuǎn)移矩陣為：

2009 年訓(xùn)練數(shù)據(jù)處于［-5%～0）狀態(tài)，確定初始狀態(tài)向量S0=［0，1，0，0，0］，在經(jīng)過一步轉(zhuǎn)移后的狀態(tài)概率向量Sk=S0P(k)=［0，3/5，0，1/5，1/5］。由此可以得到2010 年預(yù)測數(shù)據(jù)相對誤差可能處于［-5%～0）狀態(tài)，即可采用公式對灰色GM（1，1）模型預(yù)測結(jié)果進(jìn)行修正。

2.3.4 BP神經(jīng)網(wǎng)絡(luò)模型

通過Matlab 編寫程序建立BP 神經(jīng)網(wǎng)絡(luò)預(yù)測模型模擬預(yù)測地下水埋深［23］。模型選用三層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，輸入層神經(jīng)節(jié)點(diǎn)設(shè)置為4個，由降水量、蒸發(fā)量、地下水開采量和前期水位組成；輸出層節(jié)點(diǎn)為1 個即研究區(qū)地下水位埋深；隱含層節(jié)點(diǎn)個數(shù)根據(jù)經(jīng)驗(yàn)公式設(shè)置為9，確定模型結(jié)構(gòu)為4∶9∶1。隱含層的傳遞函數(shù)選擇tansing 函數(shù)，輸出層傳遞函數(shù)選擇purelin 函數(shù)，學(xué)習(xí)算法采用梯度下降法trainlm 函數(shù)訓(xùn)練，設(shè)置學(xué)習(xí)速率和訓(xùn)練精度分別是0.01和0.000 01［24］。

2.3.5 基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)模型

基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)模型采用和BP神經(jīng)網(wǎng)絡(luò)相同的網(wǎng)絡(luò)結(jié)構(gòu)和參數(shù)設(shè)置［27］。GA-BP 神經(jīng)網(wǎng)絡(luò)模型遺傳算法部分采用實(shí)數(shù)編碼的方式，編碼長度為55，種群規(guī)模為40，進(jìn)化代數(shù)為100，選擇操作為輪盤賭法，變異概率0.7，交叉概率為0.1。

2.4 模擬結(jié)果對比分析

將5個模型得到的模擬預(yù)測值和實(shí)際地下水位埋深進(jìn)行擬合得到預(yù)測模型結(jié)果對比圖，如圖1 所示。由圖1 可知，多元線性回歸模型的預(yù)測值小于實(shí)際值；灰色GM（1，1）模型的預(yù)測值呈線性遞增；基于馬爾科夫鏈優(yōu)化的灰色GM（1，1）模型的預(yù)測值可以反映出實(shí)際值的波動變化情況；BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測值較為接近，但個別年份波動誤差較大；基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)預(yù)測值更加精準(zhǔn)，整體擬合情況最好。具體數(shù)值可見表4。

圖1 預(yù)測模型結(jié)果對比圖Fig.1 Comparison of prediction model results

表4 各模型模擬預(yù)測值表 mTab.4 Simulation and prediction values of each model

本文采用絕對誤差A(yù)E、相對誤差RE、平均絕對誤差MAE、平均絕對百分比誤差MAPE、均方誤差MSE和均方根誤差RMSE作為評價指標(biāo)對模型預(yù)測精度進(jìn)行評價，計(jì)算公式如下：

式中：Yi為實(shí)際地下水埋深值；為模擬預(yù)測值；n為預(yù)測個數(shù)。

評價指標(biāo)具體數(shù)值可見表4，預(yù)測模型的絕對誤差和相對誤差對比圖見圖2、3。由表5和圖2、3分析可知，GA-BP神經(jīng)網(wǎng)絡(luò)模型各項(xiàng)評價指標(biāo)都比其他四個模型效果好，且預(yù)測精度最高，預(yù)測結(jié)果最穩(wěn)定，說明GA-BP 神經(jīng)網(wǎng)絡(luò)模型可以很好地將遺傳算法和神經(jīng)網(wǎng)絡(luò)兩者優(yōu)勢相結(jié)合，既增強(qiáng)全局搜索能力又保持良好的容錯性和泛化能力。

表5 模型性能參數(shù)Tab.5 Performance parameters of the model

圖2 預(yù)測模型絕對誤差圖Fig.2 Absolute error diagram of prediction model

圖3 預(yù)測模型相對誤差圖Fig.3 Relative error diagram of prediction model

5 種預(yù)測模型精度和穩(wěn)定性按由低到高的順序依次為：多元線性回歸＜灰色GM（1，1）＜Markov-灰色GM（1，1）＜BP 神經(jīng)網(wǎng)絡(luò)＜GA-BP神經(jīng)網(wǎng)絡(luò)。多元線性回歸在考慮解決簡單的線性問題上效果較好，但影響地下水埋深條件較為復(fù)雜，相關(guān)因素具有非線性的特點(diǎn)。所以本文在考慮3 個主要因素的影響下，預(yù)測結(jié)果精度較差?；疑獹M（1，1）模型預(yù)測結(jié)果偏向于數(shù)據(jù)整體的長期發(fā)展趨勢［8］，對于均勻增減且無波動的序列預(yù)測精度效果較好，所以本次預(yù)測結(jié)果呈現(xiàn)線性遞增現(xiàn)象，誤差較大。基于馬爾科夫鏈理論優(yōu)化的灰色GM（1，1）模型通過轉(zhuǎn)移概率對原始數(shù)據(jù)序列波動進(jìn)行修正，相較于灰色GM（1，1）模型提高預(yù)測精度，更加可靠。BP神經(jīng)網(wǎng)絡(luò)對長序列資料的學(xué)習(xí)能力較強(qiáng)，整體擬合較好。GA-BP 神經(jīng)網(wǎng)絡(luò)的預(yù)測值與實(shí)際值擬合情況最好相較于其他4種預(yù)測模型可以更高效準(zhǔn)確地對地下水位埋深進(jìn)行預(yù)測。

由表4 可知BP 神經(jīng)網(wǎng)絡(luò)預(yù)測的相對誤差最大值為6.61%最小值為0.13%，而GA-BP 神經(jīng)網(wǎng)絡(luò)預(yù)測的相對誤差最大值為3.06%最小值為0.04%，表明經(jīng)過遺傳算法優(yōu)化的BP 神經(jīng)網(wǎng)絡(luò)模型可以顯著提升預(yù)測模擬精度和穩(wěn)定性；由圖4、5 所示，BP神經(jīng)網(wǎng)絡(luò)經(jīng)過五輪訓(xùn)練迭代得到最優(yōu)解，而GA-BP神經(jīng)網(wǎng)絡(luò)僅通過兩輪訓(xùn)練，均方根誤差就能達(dá)到0.024 8，快速收斂且穩(wěn)定，表明遺傳算法優(yōu)化的BP 神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練效率更高。由圖6可知，個體適應(yīng)度隨著遺傳迭代次數(shù)增加不斷減小，而個體適應(yīng)能力不斷增強(qiáng)，在經(jīng)過27次迭代后，個體適應(yīng)度保持穩(wěn)定。

圖4 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練圖Fig.4 BP neural network training diagram

圖5 GA-BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練圖Fig.5 GA-BP neural network training diagram

圖6 適應(yīng)度曲線圖Fig.6 Fitness curve

GA-BP 神經(jīng)網(wǎng)絡(luò)相較于BP 神經(jīng)網(wǎng)絡(luò)的優(yōu)勢體現(xiàn)在：預(yù)測模擬結(jié)果更加精準(zhǔn)、訓(xùn)練效率更加高效和適應(yīng)能力更強(qiáng)。遺傳算法在訓(xùn)練樣本階段通過全局尋優(yōu)和快速收斂的能力優(yōu)化BP神經(jīng)網(wǎng)絡(luò)權(quán)值和閾值，解決了BP 神經(jīng)網(wǎng)絡(luò)易陷入局部最優(yōu)、收斂能力差等問題，同時迭代優(yōu)選出的適應(yīng)力強(qiáng)個體，提高預(yù)測模型的穩(wěn)定性。

3 結(jié)論

（1）本文選用多元線性回歸、灰色GM（1，1）、基于馬爾科夫鏈優(yōu)化的灰色GM（1，1）、BP 神經(jīng)網(wǎng)絡(luò)和基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)五種預(yù)測模型，對肇州縣地下水位埋深進(jìn)行模擬預(yù)測，對比分析五種模型的預(yù)測結(jié)果和適用范圍。結(jié)果表明：GA-BP 神經(jīng)網(wǎng)絡(luò)模型相對誤差在0.04%～3.06%之間，絕對誤差在0.003～0.24之間，預(yù)測精度較高，穩(wěn)定性好。

（2）本文選擇降雨量、蒸發(fā)量、地下水開采量和前期水位作為模型的輸入層輸入應(yīng)用于肇州縣地下水埋深預(yù)測，預(yù)測結(jié)果較為準(zhǔn)確，可以為其他地區(qū)模型輸入?yún)?shù)的選擇提供參考。

（3）對于地下水位埋深預(yù)測，組合模型兩兩結(jié)合可以優(yōu)勢互補(bǔ)，提高預(yù)測結(jié)果的準(zhǔn)確性和穩(wěn)定性，達(dá)到較小誤差科學(xué)預(yù)測的效果，相較于單一方法具有較高的可行性和實(shí)用性。