茹哲敏
(山西省水文水資源勘測總站,山西 太原 030001)
地下水是山西省重要的基礎資源和戰略資源,是城鄉生活和工農業用水的重要供水水源,是生態與環境的主要控制性要素。對地下水位進行預測和分析,是地下水監督管理的重要方式,也是水利行業強監管的重要內容。定襄縣盆地內地下水的補給源包括大氣降水垂直入滲補給、基巖裂隙水的側向滲透補給、農田灌溉入滲補給、河道滲漏補給等。地下水位變化特征復雜多變,需要對各影響因素進行相關性分析,確定主要的影響因素,建立模型,對地下水位進行預測。回歸模型和神經網絡模型是2 種較為常見的預測分析模型,2 種模型各有優缺點和適用條件,其預測結果可以相互對比和印證。2 種模型預測成果可指導和管理定襄縣地下水資源的合理開發利用,對其他區域地下水位預測具有參考價值。
定襄縣位于山西省中北部,東西長48 km,南北寬36 km,滹沱河由西向東橫貫全境,全縣地形由東向西呈簸箕形,北、東、南三面群山環繞,西部和中部為忻定盆地,地勢較為平坦,盆地面積約406 km2。境內主要河流為滹沱河、牧馬河、云中河等。盆地基底為起伏不平的塊狀陷落和隆起,地下水位埋深較淺,地下水位年內變化較為平緩[1],含水層厚度大、顆粒粗,地下水總體由盆地周邊山前傾斜平原區向沖積平原區匯集,并沿滹沱河谷區由上游向下游逐漸徑流。定襄縣年降水量450 mm左右,降水量較為穩定。本次研究共選取10處地下水監測站和1處雨量站的水文年鑒資料作為研究資料,地下水監測站均勻分布,監測類型全部為淺層孔隙水,開采量使用統計年報數據,采用1993—2016 年的系列資料,研究區域及監測站分布如圖1所示。

圖1 研究區域示意
地下水水位變動受降水、開采、蒸發以及補給能力等多因素影響。降水是地下水的主要補給來源[2],豐水年的地下水位一般會明顯升高,而枯水年的地下水位一般會下降。在盆地區的潛水含水層,降水時間與地下水位變化的時間較為接近,水位變化相對于降水時間上的滯后現象不太明顯。地表水位的變化也影響地下水位的變化,一般兩者升降方向一致但存在時間滯后。人類對地下水的開采、灌溉、人工補給等活動直接影響地下水位變化,如果過量開采地下水而缺少補給,會使地下水位持續下降。
地下水位變化特征復雜多變,但也有規律可循,經初步分析可知,影響定襄縣地下水位的主要因素包括上年末水位、本年降水量、上年降水量、開采量、蒸發量等。將本年與上年水位差(水位變幅)與降水量建立關系,如圖2 所示,水位變幅與降水量呈正相關,說明降水量為重要影響因素。經分析本年降水量和開采量的關系發現,降水量與開采量呈負相關,說明降水量大的年份開采量減少,比較符合當地實際情況。再利用此方法分析本年末地下水位與其他因素之間的關系。結果表明,本年末地下水位與上年降水量的關系為負相關,說明降水影響地下水位的滯后性并不明顯,故不考慮上年降水量;同時,蒸發對地下水位影響非常微小,故不考慮蒸發量;同時,經以上分析,本次研究將上年末水位、本年降水量、開采量等主要影響因素作為模型參數。

圖2 水位變幅與降水量關系
回歸分析是研究一個隨機變量Y對另一組變量X1,X2,…,Xn之間的相互關系的一種統計方法,通過建立統計模型來研究這種關系,它主要通過對變量的觀察所獲得的統計數據來確定反映變量間關系的經驗公式,并通過所得公式進行統計描述、分析和推斷,從而解決預測、控制和優化問題。用來進行回歸分析的數學模型方程為:
式中:Y為因變量;X1,X2,…,Xn為一組自變量。
若因變量和自變量之間是線性關系,則稱為線性回歸模型,只含有一個回歸變量的線性回歸模型為一元線性回歸模型,含有多個回歸變量的線性回歸模型為多元線性回歸模型[3]。線性回歸可以通過最小二乘法求出其方程。
相關系數r可以反映變量之間相關程度。判定本年地下水位與其他因素之間是否存在相關關系,需要計算相關系數。如果相關系數接近1表示正相關,如果接近-1 表示負相關,如果接近0 表示不相關。相關系數計算公式為:
式中:Cov(X,Y)為X與Y的協方差;Var(X)為X的方差;Var(Y)為Y的方差。
由于地下水位與多種因素存在關系,因此采用多元線性回歸模型,影響本年末地下水位的主要因素包括上年末水位、本年降水量、開采量等。選用定襄1993—2016 年的監測數據作為分析數據,將2017、2018 年作為預測年份。利用數據分析工具建立地下水動態多元線性回歸模型,根據資料情況,率定模型系數。確定后的模型方程為:
式中:Hb為本年末地下水位(m);Hs為上年末地下水位(m);Pb為本年降水量(mm);Q為開采量(萬m3);a,c,d,e為模型系數,應滿足以下條件:a>0,c>0,d<0。
回歸方程中的自變量都應該與因變量是顯著相關的,且各個自變量之間應該是相互獨立的。統計檢驗是利用統計學的理論檢驗回歸方程的可靠性,包括擬合優度檢驗、模型的顯著性檢驗(F檢驗)和模型參數的顯著性檢驗(t檢驗)等。
分析變量是否對本年地下水位有影響,可通過將本年地下水位的總離差平方和進行分解來確定各因素對本年地下水位的各樣本之間的差異做出的貢獻。總的離差平方和SST可分解為回歸平方和SSR及殘差平方和SSE,即SST=SSR+SSE。F檢驗統計量公式為:
式中:n為樣本數量;p為自變量數量;SSR為回歸平方和;SSE為殘差平方和。
當假設成立時,F服從自由度為p和n-p-1的F分布。可以利用F統計量對回歸方程的總體顯著性進行檢驗。如果通過檢驗,說明全體自變量在整體上對本年地下水位是有影響的。顯著性檢驗的方差分析,詳見表1。

表1 顯著性檢驗的方差分析
在通過回歸方程的顯著性檢驗后,還需要進一步對每個自變量的顯著性進行t檢驗。如果某個自變量未通過t檢驗,說明其對本年地下水位影響不顯著,應把此自變量從回歸方程中剔除。模型計算參數,詳見表2。

表2 模型計算參數
可決系數R2是擬合優度的度量,取值范圍為[0,1],R2越接近1,說明擬合值的離差平方和占總的離差平方和的比重越大,回歸擬合效果越好。可決系數分解公式為:式中:SSR為回歸平方和;SSE為殘差平方和;SST為總離差平方和。
復相關系數R是R2的平方根,表示回歸方程對原數據擬合程度的好壞,它反映了所有自變量與因變量之間線性關系及密切程度。經計算得出復相關系數R為0.857,可決系數R2為0.734。
根據模型檢驗情況,可決系數R2接近1,說明實際水位與預測水位擬合效果較好。參數a>0,表示本年末地下水位與上年末地下水位正相關;參數c>0,表示本年末地下水位與本年降水量正相關,降水量的增加使地下水位上升;參數d<0表示本年末地下水位與開采量負相關,地下水的開采會導致地下水位的下降。上年末地下水位、本年降水量、開采量的t檢驗的P值都<0.05,則表示顯著性較好。因此,這3個參數較為合理。通過各項分析,修正不合理數據,剔除不顯著的變量,對剩余因素重新建立回歸方程,直到通過各項檢驗為止。最終形成的回歸方程式為:
將上年末地下水位、本年降水量、開采量代入統計模型,計算出預測地下水位,再與實際地下水位相比較。2017 和2018 年預測地下水位分別為753.52和753.75 m,與實際地下水位相差分別為-0.01 和-0.14 m,預測結果較為良好。實際地下水位與預測地下水位對比情況,如圖3所示。

圖3 實際地下水位與回歸模型預測地下水位對比
采用誤差反向傳播(Error back propagation)算法的神經網絡稱為BP神經網絡。網絡由輸入層、中間隱含層和輸出層組成,中間層可以有一層或多層,層與層之間全連接,同一層之間的節點無連接。算法的基本過程由信號的正向傳播與誤差反向傳播2個過程組成。正向傳播時,輸入數據從輸入層傳入,經過各隱含層處理后逐層向后傳播,傳向輸出層。根據指定的初始化權重值和偏置值來計算輸出值,若與期望輸出相差過大,則轉入誤差的反向傳播階段。誤差反向傳播是將輸出誤差通過隱含層向輸入層逐層反傳,并將誤差分攤給各層的所有單元,從而獲得各層的誤差信號,此誤差信號即作為修正單元權值的依據。這種信號正向傳播與誤差反向傳播的各層權值調整過程周而復始地進行,權值不斷調整,一直進行到網絡輸出的誤差減少到可預先設定的值或進行到預先設定的迭代次數為止[4]。迭代結束后,得出最優參數即性能函數取最小值的參數,包括最終權重矩陣和偏置,即可用此參數進行預測。
BP 神經網絡優化算法有很多,包括動量BP 法、學習率可變的BP 算法、共軛梯度法、牛頓法和擬牛頓法、Levenberg-Marquardt 算法等。選擇Levenberg-Marquard 算法來計算性能函數的最小值,它類似擬牛頓算法,同時具有梯度法和牛頓法的優點,根據多次試驗對比,相對于其他算法,它的收斂速度較快,不容易陷入局部最小,算法結果更合理。
網絡由輸入層、中間隱含層和輸出層組成。將影響本年末地下水位的主要因素作為網絡的輸入層數據,輸入層節點數為3 個,分別為上年末地下水位、本年降水量、本年開采量。輸出層節點數為1個,輸出層數據為本年末地下水位。采用1993—2016 年系列資料作為輸入數據,將2017 和2018 年作為預測年份,分析掌握兩者之間潛在的規律,最終根據這些規律來推算2017 和2018 年輸出結果。通過Matlab?軟件構建BP 神經網絡模型,模型結構如圖4所示。

圖4 BP神經網絡模型結構
需要最優化的函數為性能函數,用來評價網絡的精確程度。選擇均方誤差MSE函數作為性能函數,均方誤差是指預測值與實際值之差的平方和的平均值,它包含了偏差和方差[5]。均方誤差函數公式為:
式中:X'i為預測輸出數據;Xi為實際監測數據;n為樣本數量。
由于模型因素的數值單位不統一,為避免大數值因子弱化小數值因子,要進行歸一化處理,將數據按比例縮放。當網絡訓練結束后,再反歸一化到原始數據范圍。本次采用最大最小法進行歸一化,統一映射到[-1,1]區間上,歸一化公式為:
式中:Xmax為原始數據的最大值;Xmin為原始數據的最小值;X為原始數據;Y為歸一化后數據。
構建BP 神經網絡,需要對權重進行初始化,選擇激活函數,設置隱含層的層數和節點數、最大迭代次數、學習率等參數。具體包括以下方面。
4.3.1 初始化權重
網絡通過迭代的方式確定權重,需要一個初始值。較大的初始權重會造成大的信號傳遞給激活函數,導致網絡飽和,從而降低網絡學習到更好權重的能力。初始權重通常為較小的非零隨機值,不能將權重初始化為相同的恒定值或零,否則會導致反向傳播誤差得到平分,使網絡無法更新權重。采用Nguyen-Widrow 算法對每個權重矩陣和偏置進行初始化,使得每層節點的活動區域大致均勻地分布在輸入空間上,由于很少有神經元被浪費,訓練速度更快,相對于純隨機權重和偏置,此算法更具有優勢[6]。
4.3.2 激活函數
激活函數可選擇Sigmoid 函數、線性函數、ReLU函數、Softplus 函數等函數。隱含層選用tanSigmoid函數,輸出層選用純屬線性函數。Sigmoid 函數是光滑、可微、單調的函數[7],具有非線性的放大功能,tanSigmoid 函數是雙曲正切Sigmoid 函數,函數將輸入從正負無窮的范圍映射在(-1,1)區間,當輸入的絕對值非常大時會出現飽和現象,對輸入的微小變化變得不敏感。tanSigmoid函數公式為:
式中:x為自變量;e為自然常數。
4.3.3 隱含層的層數和節點數
神經網絡可以包括一個或多個隱含層,并需要確定隱含層的層數。隱含層節點數對BP 神經網絡的性能影響較大,一般采用經驗公式來確定。設置1個隱含層和5個隱含層節點數,可以滿足要求。
4.3.4 最大迭代次數
網絡在訓練時,如果達到最大迭代次數,即使達不到誤差要求,也會終止計算。設置最大迭代次數為2 000,設置性能值為0.000 01。
4.3.5 學習率η
學習率決定了在一次學習中多大程度上更新參數。學習率過大,容易出現振蕩,可能導致學習過程不穩定而不能收斂,無法達到要求的精度;學習率過小,又可能會導致訓練周期增加、收斂慢。剛開始可先選擇較小的學習率訓練網絡以保持網絡穩定運行,同時密切觀察誤差下降曲線的變化情況,若下降較快,則說明學習率選取恰當;但若出現比較大的震蕩,則說明學習率選得略大,可適當調小。學習率的選取直接影響閾值和權值的調整量,也影響網絡收斂速度和精度。通過試驗,學習率設置為0.001較為合適。
模型參數設置完成后,將歸一化處理后的訓練數據輸入網絡,調用模型程序進行訓練,訓練期間沿著性能函數最快下降的方向,不斷地調整網絡的權重和偏置,直到性能函數值小于設置值。經過多次試算,得出最優參數即性能函數取最小值的參數,包括輸入權重矩陣和偏置。輸入層到中間層、中間層到輸出層的權重矩陣和偏置,詳見表3。

表3 權重矩陣和偏置
利用訓練好的模型參數,采用神經網絡算法,推算出2017、2018 年輸出結果。2017、2018 年預測地下水位分別為753.58 和753.67 m,與實際地下水位相差分別為0.05 和-0.21 m,預測結果較為良好。擬合相關性如圖5所示,最佳擬合由虛線表示,完美擬合(輸出等于目標)由實線表示,輸出和目標之間的相關系數接近1,這表明擬合較好。實際地下水位與預測地下水位對比情況,如圖6所示。

圖5 擬合相關性

圖6 實際地下水位與神經網絡模型預測地下水位對比
回歸模型可以建立因變量與自變量之間的關系,回歸分析可以相對準確地計量各個因素之間的相關程度,通過標準的統計方法可以計算出唯一的結果。運用回歸模型需要確定各影響因素,影響因素的選擇很重要,定襄縣地下水位與上年降水量的關系不大,而其他縣區地下水位可能與上年降水關系密切,如果選擇的因素不合適,將不能反映實際情況,需要根據不同研究區域的實際情況選擇合適的自變量。當地下水位與各影響因素之間的關系較為顯著時,采用回歸模型較為合適。
神經網絡模型實現了一個從輸入到輸出的映射過程,它可以實現任何復雜非線性映射的功能,適合于求解內部機制復雜的問題,在函數擬合、優化計算、模式識別及聚類等多種領域應用廣泛。但神經網絡容易得到局部最優解,如果參與訓練數據量過少,容易發生過擬合,即只能擬合訓練數據,但對其他數據不能很好擬合,可以通過增加數據量以防止過度擬合;網絡需要設置的參數較多,隱含層數和節點數的選擇無理論依據,只能通過經驗來給定一個粗略的范圍;對權重初始值敏感,由于初始權重是隨機給定的,網絡訓練具有不可重現性,每次訓練的結果會不同,但可保存權重矩陣和偏置,以便下次預測重復使用。當無法確定地下水位與各影響因素之間的公式時,采用BP神經網絡模型較為合適。
(1)監測資料的可靠性對2 種模型預測的精度至關重要。降水量和地下水位監測資料的時序較長、系列資料較完整,其預測精度相對較高,而地下水開采量以收集統計為主,開采量數據應進行進一步分析和復核。選用地下水監測站水位數據應與開采層位相對應,且監測質量較好,能夠代表附近水位變化情況,如果與開采層位不對應,計算結果將嚴重失真,不能反映區域水位變化規律。
(2)2 種模型樣本的數量應足夠多,否則可能會得出局部最優解,而不適用全局,對預測精度將產生不利影響。定襄縣的監測資料選取1993—2016 年,資料系列年限較長,不會出現過度擬合現象,模型預測精度能夠得到保證。
(3)本文建立了回歸模型和神經網絡模型,用這2 種模型對定襄縣忻定盆地的地下水位進行了預測,預測結果基本可靠,2 種模型預測結果可以相互對比和印證,若相差較大,應詳細分析原因并調整模型參數,直至結果合理。本文研究方法可指導并管理定襄縣地下水資源的合理開發利用,對其他區域地下水位預測具有參考價值。