董凡，張曉影，陳廣泉,，戴振學(xué)，王延誠, 3

(1.吉林大學(xué) 建設(shè)工程學(xué)院，吉林 長春 130026；2.自然資源部第一海洋研究所 自然資源部海岸帶科學(xué)與綜合管理重點實驗室，山東 青島 266061；3.青島海洋科學(xué)與技術(shù)試點國家實驗室 海洋地質(zhì)過程與環(huán)境功能實驗室，山東 青島 266237)

1 引言

近年來，受氣候變化和人類活動的雙重影響，全球海岸地區(qū)都面臨著海水入侵的威脅。海水入侵可引發(fā)地下水水質(zhì)惡化、土壤退化等生態(tài)環(huán)境問題，嚴(yán)重制約社會的可持續(xù)發(fā)展[1–4]。科學(xué)管理沿海地區(qū)地下水，預(yù)測近岸地下水位和電導(dǎo)率，掌握地下水位及鹽度的動態(tài)變化，對合理開發(fā)利用海岸帶地下水資源，防止海水入侵加劇，保護(hù)海岸帶生態(tài)環(huán)境有十分重要的指導(dǎo)意義。然而，受自然條件和人為因素的影響，海水入侵在時間序列上表現(xiàn)出時空性、不確定性及動態(tài)變化等特性，多種因素對海水入侵的影響并非多因子影響的線性疊加，從而使預(yù)測海水入侵的演化過程變得十分困難。

隨著人工智能的興起，由于人工神經(jīng)網(wǎng)絡(luò)（ANN）可以較好地映射非線性變化，因此被廣泛應(yīng)用在水文變量預(yù)測中。其中BP 神經(jīng)網(wǎng)絡(luò)模型的應(yīng)用范圍最廣[5]，但其神經(jīng)網(wǎng)絡(luò)訓(xùn)練的時間長，收斂速度慢，輸入信號由輸入層單向傳遞給輸出層，網(wǎng)絡(luò)模型只能使用有限步的歷史信息。而循環(huán)神經(jīng)網(wǎng)絡(luò)（RNN）憑借其循環(huán)結(jié)構(gòu)的特點，通過反饋連接可以使用長期的歷史信息，因此更加適合時間序列的預(yù)測[6]。具有外部輸入的非線性自回歸神經(jīng)網(wǎng)絡(luò)（Nonlinear Autoregressive Models with Exogenous Inputs，NARX）作為循環(huán)神經(jīng)網(wǎng)絡(luò)的一種，在地下水位的預(yù)測研究領(lǐng)域備受青睞。Guzman 等[7]使用NARX 模型預(yù)測密西西比河流域含水層，Wunsch 等[8]運用NARX 模型預(yù)測德國西南部數(shù)口井的地下水水位，Di Nunno 和Granata[9]應(yīng)用NARX網(wǎng)絡(luò)預(yù)測位于普利亞地區(qū)的76 口監(jiān)測井的日地下水水位。研究結(jié)果表明，NARX 神經(jīng)網(wǎng)絡(luò)在預(yù)測時間序列上具有顯著優(yōu)勢。

盡管ANN 在時間序列預(yù)測方面應(yīng)用廣泛，具有靈活性，但是如果不對輸入和輸出數(shù)據(jù)進(jìn)行預(yù)處理，它可能無法處理非平穩(wěn)數(shù)據(jù)[10–11]。小波變換作為一種數(shù)據(jù)預(yù)處理技術(shù)，可以將非平穩(wěn)時間序列的小波分解成不同的尺度，提供對序列結(jié)構(gòu)的解釋，并僅用幾個系數(shù)來提取有關(guān)其歷史的重要信息。該技術(shù)被廣泛應(yīng)用于非平穩(wěn)信號的時間序列分析[12]，通過將非平穩(wěn)信號分解成不同時間水平的子信號，可以較好地預(yù)測非平穩(wěn)信號[13–18]。將小波變換和人工神經(jīng)網(wǎng)絡(luò)模型相結(jié)合，發(fā)揮這兩種模型的優(yōu)勢，可以明顯提高預(yù)測精度，增加有效預(yù)測時段。目前，混合模型已廣泛應(yīng)用于部分水文要素的預(yù)測，如降雨、徑流和蒸發(fā)[19]。張建峰等[20]建立小波–BP 神經(jīng)網(wǎng)絡(luò)的混合模型預(yù)測華北平原城市地下水位。Seo 等[21]將小波分析（WA）分別結(jié)合BP 神經(jīng)網(wǎng)絡(luò)、自適應(yīng)神經(jīng)模糊推理系統(tǒng)，預(yù)測水庫水位并對結(jié)果進(jìn)行對比。Zhang 等[22]構(gòu)建了小波分析和NARX 神經(jīng)網(wǎng)絡(luò)混合模型對舟山填海區(qū)地下水位進(jìn)行預(yù)測，預(yù)測精度較高。

由于受全新世海侵、氣候干旱、海平面上升以及不合理的人類活動的影響，萊州灣地區(qū)是我國海水入侵災(zāi)害最為嚴(yán)重和典型的地區(qū)[23]。近年來，在上述因素綜合影響下，萊州灣地區(qū)地下咸水體分布面積已超過22 000 km2[24]，嚴(yán)重污染地下淡水資源，土壤次生鹽漬化面積不斷擴(kuò)大，對海岸帶區(qū)域的工農(nóng)業(yè)生產(chǎn)和生態(tài)環(huán)境造成嚴(yán)重影響，海水入侵已成為該地區(qū)嚴(yán)重的環(huán)境問題之一，嚴(yán)重制約經(jīng)濟(jì)社會發(fā)展。因此，開展萊州灣地區(qū)海水入侵動態(tài)特征分析，建立預(yù)測模型，對提高海岸帶地區(qū)防災(zāi)減災(zāi)能力和經(jīng)濟(jì)、社會可持續(xù)發(fā)展具有重要意義。在本研究選擇萊州灣南岸典型監(jiān)測剖面，首先分析地下水位和地下水電導(dǎo)率的動態(tài)變化特征，分析對比不同輸入?yún)?shù)對模型預(yù)測精度的影響，建立基于WA-NARX 混合的海水入侵預(yù)測模型，對近岸地下水位和電導(dǎo)率進(jìn)行預(yù)測模擬，研究結(jié)果可為海岸帶海水入侵的實時監(jiān)測、預(yù)測、預(yù)警提供理論與技術(shù)支撐。

2 研究區(qū)與研究方法

2.1 研究區(qū)與資料

2.1.1 研究區(qū)概況

研究區(qū)位于萊州灣南部（圖1），濰坊北部的濱海區(qū)。研究區(qū)內(nèi)設(shè)有6 口監(jiān)測井（BH01–BH06），分別沿垂直于海岸帶的剖面由靠近海岸帶的監(jiān)測井BH06 向內(nèi)陸的監(jiān)測井BH01 方向分布。研究收集了降雨、蒸發(fā)、潮汐、電導(dǎo)率、地下水位和農(nóng)業(yè)排灌用電量數(shù)據(jù)，見圖2。其中，潮汐的變化是潮位站監(jiān)測的海水位每小時的變化，降雨量與蒸發(fā)量是日測數(shù)據(jù)，排灌用電量是研究區(qū)每小時灌溉所用電量。研究選 取2011 年10 月1 日 至2012 年10 月31 之 間 的 監(jiān)測數(shù)據(jù)作為研究對象。

圖1 研究區(qū)位置圖及監(jiān)測井分布Fig.1 Location of the study area and monitoring wells

圖2 各影響因素時間序列Fig.2 Time series of influencing factors

2.1.2 研究區(qū)水文地質(zhì)特征

萊州灣南岸地區(qū)屬于魯北平原水文地質(zhì)區(qū)、濰北濱海平原水文地質(zhì)亞區(qū)。該地區(qū)地表水系較發(fā)達(dá)，地下水水文地質(zhì)特征在水平和垂向上均變化明顯。地貌類型由南向北依次為山前沖洪積平原、沖積海積平原和濱海海積平原，海岸類型為粉砂淤泥質(zhì)海岸。此處的沖洪積平原由多個扇緣相互交接的沖洪積扇群組成，沖洪積平原的前緣被海積物或黃河沖積物所掩埋，含水層厚度大、顆粒粗，具有多層結(jié)構(gòu)。含水層顆粒自上而下呈現(xiàn)由粗變細(xì)的趨勢，而在水平方向上則具有沖洪積扇的水文地質(zhì)特征[25]。

如圖3 所示，研究區(qū)普遍存在一層中更新統(tǒng)砂質(zhì)黏土弱透水層，埋深一般為 40～70 m，層厚約20 m。有3 層海侵層，分別對應(yīng)晚更新世以來發(fā)生的3 次大海侵。在3 次海侵中都滯留了濃縮于地層中的古海水（咸–鹵水），與淡水含水層在水平方向上相連，是該區(qū)咸–鹵水入侵的來源，主要出現(xiàn)在晚更新世和全新世地層中。

圖3 研究區(qū)剖面圖（修改自文獻(xiàn)[26]）Fig.3 Hydrogeology of the cross section from south to north in the study area (modified from reference [26])

2.2 研究方法

2.2.1 小波分析

在本研究中，運用小波分析的方法來識別顯著的振蕩周期，從原始數(shù)據(jù)中提取有用的信息，并將這些信息分解成一系列小波函數(shù)，小波函數(shù)由母波函數(shù)生成。本文采用的母波函數(shù)是Morlet 小波：

式中，w0為無量綱頻率；t為時間。通常，選擇w0=6的Morlet 小波[27–28]。連續(xù)小波變換（CWT）用于識別信號在時頻域的主振蕩。連續(xù)小波變換的f(t)如下所示：

式中，a為尺度或頻率參數(shù)；b為時間或位移參數(shù);(t) 是 復(fù)共軛函數(shù) ψ(t)。由式（2）可知，小波變換在不同分辨率或尺度下呈現(xiàn)出不規(guī)則和不對稱的特點，將f(t)分解為一系列低分辨率分量。

在實際應(yīng)用中，連續(xù)小波變換常被離散小波變換（DWT）所代替。令a0>1,b0∈R,k和j為整數(shù)，則f(t)的離散小波變換可表示為

在式（3）中，參數(shù)a0和b0最常用的選擇是使用對數(shù)縮放[29]，即a0=2,b0=1。原始時間序列在頻率（a或j）和時間（b或k）域的特征被連續(xù)小波變換的W(a,b)或離散小波變換的W(j,k)所反映。

將經(jīng)典的連續(xù)小波轉(zhuǎn)換為離散小波有多種算法，如Mallat 算法和á Trous 算法。本文采用了á Trous 算法[30]。á Trous 的分解算法定義如下：

式中，g(n)為離散低通濾波器；Ci(t) 和Wi(t) (i=1,2,···)分別是分辨率為i級的尺度系數(shù)和小波系數(shù)。W1(t),W2(t),···,Wp(t) 和Cp(t)組成分辨率為P的離散小波變換。

在小波分析中經(jīng)常使用近似特征和細(xì)節(jié)分量，近似特征表示信號的低尺度，即低頻信息；細(xì)節(jié)分量表示信號的高尺度，即高頻信息。離散小波變換使用低通濾波器g(n)和高通濾波器h(n)將f(t)分解成不同尺度上的低頻成分（C1，C2，C3··· ）和高頻成分（W1，W2，W3···）。因此，f(t)可以分解為包含其近似分量（低頻成分）和細(xì)節(jié)分量（高頻成分）。高頻信號捕獲原始時間序列數(shù)據(jù)中的小尺度特征，見圖4 小波分解部分。

圖4 WA-NARX 模型建模流程Fig.4 Flow chart of WA-NARX neural network hybrid model

2.2.2 NARX 神經(jīng)網(wǎng)絡(luò)

NARX 神經(jīng)網(wǎng)絡(luò)是一種有外部輸入的非線性自回歸網(wǎng)絡(luò)模型[31]，表示為

式中，y(t)代表輸出信號；x(t)代表外部輸入信號；f為非線性函數(shù)。即輸出信號y(t)可由輸出信號前ny個值和外部輸入信號x(t)前nu個值通過f非線性映射得到。

NARX 神經(jīng)網(wǎng)絡(luò)由網(wǎng)絡(luò)的輸出反饋和動態(tài)的神經(jīng)元構(gòu)成，是一種動態(tài)神經(jīng)網(wǎng)絡(luò)。它的結(jié)構(gòu)清晰，通常包含輸入層、輸出層、隱含層及輸入延時和輸出延時這些網(wǎng)絡(luò)結(jié)構(gòu)（圖5）。NARX 神經(jīng)網(wǎng)絡(luò)具有開環(huán)和閉環(huán)兩種模式，閉環(huán)網(wǎng)絡(luò)可以進(jìn)行多步預(yù)測。與一般神經(jīng)網(wǎng)絡(luò)相比，NARX 神經(jīng)網(wǎng)絡(luò)的反饋結(jié)構(gòu)和梯度下降算法使其擁有更好的收斂速度和學(xué)習(xí)效率，是一種研究較為廣泛的動態(tài)神經(jīng)網(wǎng)絡(luò)[32–33]，因此，NARX 神經(jīng)網(wǎng)絡(luò)既可以用于時間序列建模，也可以用于非線性濾波[34]。

圖5 輸入到輸出的NARX 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.5 Nonlinear auto-regression with exogenous input NARX

2.2.3 小波與NARX 混合模型

WA-NARX 混合模型是利用多分辨率小波分解得到的子時間序列作為輸入的NARX 模型。利用小波變換將原始信號依次分解成不同尺度上的信號成分，á Trous 算法產(chǎn)生低頻的近似分量（C）和高頻的細(xì)節(jié)分量（W）[35]。本研究利用L=int來確定小波的分解層數(shù)[36]，其中L表示小波層的個數(shù)，n表示時間序列中的數(shù)據(jù)點個數(shù)，int[ ]為取整函數(shù)。研究所包含時間序列的數(shù)據(jù)點個數(shù)為9 528，經(jīng)計算需要4 個分辨率級別。

本研究中，利用á Trous 算法首先將水位、潮汐和降雨等參數(shù)的小時時間序列分解為一定數(shù)量的子時間分量。根據(jù)潮汐水位隨時間的變化規(guī)律，運用小波分析的方法來識別顯著的振蕩周期，并從原始數(shù)據(jù)中提取潮汐水位變化的近似特征值和細(xì)節(jié)分量，然后將時間分量作為混合模型的輸入。選擇最重要的時間分量作為輸入，這些時間分量與原始數(shù)據(jù)具有很高的相關(guān)性。為了判別預(yù)測的準(zhǔn)確性，本文將均方根誤差（RMSE）和擬合度（R2）作為性能指標(biāo)來衡量模型的預(yù)測精度。同時，利用matlab 中的perform 函數(shù)結(jié)合網(wǎng)絡(luò)結(jié)構(gòu)與均方誤差（MSE）計算網(wǎng)絡(luò)性能來評價網(wǎng)絡(luò)的學(xué)習(xí)及預(yù)測能力。

WA-NARX 混合模型在建模過程中，數(shù)據(jù)被分為兩個子集：訓(xùn)練集和測試集。訓(xùn)練數(shù)據(jù)集（75%）用于優(yōu)化網(wǎng)絡(luò)配置，測試數(shù)據(jù)集（25%）用于測試網(wǎng)絡(luò)性能。此外，WA-NARX 混合模型的隱含層神經(jīng)元數(shù)是通過反復(fù)試驗確定的。在本次研究中，最終確定10 個隱含層神經(jīng)元。WA-NARX 模型建模流程見圖4，如圖所示：首先是參數(shù)和樣本的選擇，利用小波變換將原始時間序列分解成不同信號成分的小波分量，然后將分解的小波分量輸入到NARX 神經(jīng)網(wǎng)絡(luò)模型中進(jìn)行訓(xùn)練預(yù)測，最終得到預(yù)測結(jié)果。

3 結(jié)果和討論

3.1 海水入侵動態(tài)變化特征

3.1.1 地下水位動態(tài)變化特征

根據(jù)研究區(qū)監(jiān)測井2011–2012 年間的地下水水位監(jiān)測數(shù)據(jù)，結(jié)合研究區(qū)內(nèi)相應(yīng)時段的降雨量繪制監(jiān)測點的地下水位動態(tài)變化曲線（圖6）。根據(jù)研究區(qū)2011–2012 年各監(jiān)測井采集的地下水位動態(tài)數(shù)據(jù)，計算各監(jiān)測井地下水位的最大值、最小值、平均值以及標(biāo)準(zhǔn)差等空間變化統(tǒng)計特征，計算結(jié)果見表1。圖6 顯示，各監(jiān)測井地下水位的變化趨勢基本一致。監(jiān)測井在2011–2012 年間的地下水位的變幅為1.175～4.358 m，各監(jiān)測井地下水位在2011 年12 月、2012 年4 月、2012 年8 月這3 個時間段內(nèi)均有突變點，與降雨量變化相對應(yīng)。在降雨后，多數(shù)監(jiān)測井地下水位有明顯上升的趨勢。2011 年的11 月份到2012年的2 月份水位逐漸上升后趨于平緩，因降雨補給存在滯后效應(yīng)，隨后水位開始下降，2012 年3–4 月份水位隨著降雨上升，在4 月末水位達(dá)到全年最高值，因春灌開采影響，5 月份雖有降雨補給，但地下水位也開始下降，7–8 月是降雨最頻繁月份，受強降雨的影響，地下水位隨降雨量波動變化。在降雨過程中，間歇性的小雨只能濕潤土壤表層，不能有效補給地下水；集中式暴雨入滲系數(shù)低，大部分轉(zhuǎn)化為地表徑流。而連續(xù)小雨可以充分入滲到地下，因此有利于地下水補給[37]?？梢耘袛啾O(jiān)測井2011–2012 年地下水年內(nèi)動態(tài)變化特征為降雨入滲–開采型。

圖6 監(jiān)測井2011–2012 年間地下水位變化曲線Fig.6 Groundwater level variation in monitoring wells from 2011 to 2012

表1 監(jiān)測井地下水位統(tǒng)計特征分析（單位：m）Table 1 Analysis on statistical characteristics of groundwater level in monitoring wells (unit:m)

地下水超采是海水入侵最直接的原因，地下水位的持續(xù)下降會破壞水壓力平衡，當(dāng)沿海含水層地下水位低于海平面時，就會發(fā)生海水入侵，因此在分析區(qū)域內(nèi)海水入侵狀況時，可以根據(jù)地下水位的動態(tài)變化近似推斷海水入侵狀況。對海水入侵的監(jiān)測，通常是對地下水中氯化物和礦化度進(jìn)行分析[38]，也有相關(guān)研究[39–40]用電導(dǎo)率替代Cl–作為監(jiān)測海水入侵的指標(biāo)，指出水中Cl–含量與電導(dǎo)率成正比。如圖7 所示，地下水位和電導(dǎo)率在一定程度上呈負(fù)相關(guān)，即地下水位下降，電導(dǎo)率上升。當(dāng)?shù)叵滤幌陆禃r，淡水入海流量減小，海水入侵加劇，地下水中電導(dǎo)率升高。因此，從圖中可以看出地下水位動態(tài)變化特征在年內(nèi)時間尺度可較為準(zhǔn)確表征海水入侵的變化。

圖7 監(jiān)測井BH4 水位與電導(dǎo)率變化Fig.7 Groundwater level and electrical conductivity variation in BH4 monitoring well

3.1.2 電導(dǎo)率動態(tài)變化特征

根據(jù)研究區(qū)2011–2012 年各監(jiān)測井采集的電導(dǎo)率動態(tài)數(shù)據(jù)，計算各監(jiān)測井電導(dǎo)率的最大值、最小值、平均值、標(biāo)準(zhǔn)差以及標(biāo)準(zhǔn)離散率等空間變化統(tǒng)計特征，結(jié)果見表2。

表2 監(jiān)測井電導(dǎo)率統(tǒng)計特征分析（單位：mS/cm）Table 2 Analysis on statistical characteristics of electrical conductivity in monitoring wells (unit:mS/cm)

標(biāo)準(zhǔn)離散率θ的計算公式為

式中，SD為電導(dǎo)率的標(biāo)準(zhǔn)差；K為電導(dǎo)率的平均值。

溫度為20℃時，天然水的電導(dǎo)率介于0.5～1.5 mS/cm。在濱海地區(qū)，電導(dǎo)率可以作為一項評判海水入侵的指標(biāo)，電導(dǎo)率的值越小說明地下水受海水入侵程度越低[41–42]。由上表可以得出，監(jiān)測井BH1、BH2 和BH3 的電導(dǎo)率值大致在天然水電導(dǎo)率的范圍附近，沒有出現(xiàn)電導(dǎo)率異?，F(xiàn)象，因此這幾口監(jiān)測井的地下水基本未受海水入侵的影響；而BH4、BH5 及BH6 的電導(dǎo)率值遠(yuǎn)超出該范圍，說明地下水在很大程度上受到海水入侵的影響。

從電導(dǎo)率的平均值來看，BH5 的值最高，水質(zhì)相對較差，其主要原因可能與該處的地下水開采所導(dǎo)致的海水入侵有關(guān)，距離它較近的BH4 和BH6 的值也較高；BH3 的平均電導(dǎo)率最低，為1.20 mS/cm，說明此處的地下水質(zhì)相對較好。監(jiān)測井BH4 的標(biāo)準(zhǔn)離散率最高，為47.91%，說明該點地下水電導(dǎo)率相差較大；BH2 的標(biāo)準(zhǔn)離散率最低，為1.27%，說明該點的電導(dǎo)率的穩(wěn)定性較高。

3.1.3 小波相干性分析

為了理解潮汐對地下水位和地下水電導(dǎo)率的影響，引入小波相干性識別潮汐與地下水位和電導(dǎo)率之間的統(tǒng)計關(guān)系。潮汐與地下水位和電導(dǎo)率之間的小波相干性分析結(jié)果如圖8 所示（圖中以監(jiān)測井BH2和BH6 為例），譜之間的相位關(guān)系用箭頭的方向表示。影響錐（白色虛線）表示小波譜的區(qū)域，在該區(qū)域由于基礎(chǔ)數(shù)據(jù)的有限長度特性而產(chǎn)生的邊緣效應(yīng)不能被忽略。

圖8 小波相干性分析Fig.8 Wavelet coherence analysis

（1）地下水位和潮汐之間的小波相干性分析

小波相干性分析結(jié)果表明，每口監(jiān)測井的地下水位與潮汐數(shù)據(jù)之間存在不同的相干性。監(jiān)測井BH2 和BH6，在整個時間序列中頻譜值有相應(yīng)變化（圖8a 和圖8b），這些頻率代表地下水位對潮汐變化的響應(yīng)。對于監(jiān)測井BH2 和BH6，地下水位和潮汐在0.25 d 頻率上的相干性很小，0.5 d 相干性較高，1 d的相干性相對較弱，但相關(guān)分布在不同時間是間歇性的。在0.5 d 周期內(nèi)，與監(jiān)測井BH6 同期地下水位和潮汐資料相比，監(jiān)測井BH2 中的地下水位對潮汐的響應(yīng)時間不同，相干性更弱，這可能與BH2 的位置及地質(zhì)條件有關(guān)。小波相干性分析表明，海岸帶含水層的地下水位受到海洋潮汐變化的影響。潮汐信號對地下水位的影響具有時空差異?？拷０兜牡貐^(qū)更容易受潮汐影響。同時，某些地區(qū)的潮汐效應(yīng)一般與它們相對于海洋的位置和地質(zhì)條件的綜合效應(yīng)相關(guān)[43]。

（2）地下水電導(dǎo)率和潮汐之間的小波相干性分析

圖8c 和圖8d 顯示，將小波相干性分析應(yīng)用于監(jiān)測的地下水電導(dǎo)率數(shù)據(jù)和潮位時間序列時，未檢測到連續(xù)的短期時間相關(guān)性。監(jiān)測井BH2 在2012 年7 月1 日至2012 年10 月1 日時間段內(nèi)電導(dǎo)率和潮汐之間存在16～32 d 的長周期相關(guān)性。監(jiān)測井BH6 中電導(dǎo)率和潮汐之間存在短期的相關(guān)性，可以看出不同井位的電導(dǎo)率對潮汐變化的響應(yīng)不同。潮汐對地下水電導(dǎo)率的影響要弱于對地下水位的影響，因為地下水位的波動是受應(yīng)力變化的影響，而電導(dǎo)率的變化需要含水層儲水量的體積變化。更重要的是，地質(zhì)條件對控制潮汐對地下水位的影響至關(guān)重要。

3.2 模擬結(jié)果分析

3.2.1 水位模擬結(jié)果分析

WA-NARX 神經(jīng)網(wǎng)絡(luò)混合模型在NARX 神經(jīng)網(wǎng)絡(luò)模型的基礎(chǔ)上，引入小波分析，對數(shù)據(jù)進(jìn)行小波分解以減小預(yù)測誤差，提高預(yù)測精度。在進(jìn)行海岸帶地下水位預(yù)測時，把降雨和潮汐作為輸入?yún)?shù)，經(jīng)WANARX 混合模型模擬得到地下水位預(yù)測結(jié)果，同水位的觀測值進(jìn)行對比，得到預(yù)測結(jié)果的散點圖和時間序列圖。圖9 顯示了BH1–BH6 這6 口監(jiān)測井地下水位的預(yù)測結(jié)果。

圖9a、圖9c、圖9e、圖9g、圖9i、圖9k 為預(yù)測結(jié)果散點圖，圖9b、圖9d、圖9f、圖9h、圖9j、圖9l 為預(yù)測結(jié)果時間序列圖。通過預(yù)測結(jié)果散點圖可得：訓(xùn)練數(shù)據(jù)和預(yù)測數(shù)據(jù)集中分布一條直線上，模型在水位變幅較大時預(yù)測效果有所降低。如圖所示，監(jiān)測井BH3 的散點出現(xiàn)小幅度的偏離。通過預(yù)測結(jié)果的時間序列圖可得：訓(xùn)練和測試過程水位變化趨勢與觀測水位變化趨勢一致?；旌夏Ｐ驮谟?xùn)練和測試過程中均表現(xiàn)出較低的誤差，在地下水位波動處誤差相對較高，但模型總體上表現(xiàn)出良好的模擬效果。WANARX 神經(jīng)網(wǎng)絡(luò)混合模型在訓(xùn)練階段和測試階段的RMSE 和R2性能指標(biāo)值見表3，從表3 可見，所有監(jiān)測井在訓(xùn)練和測試階段的R2均大于0.99，RMSE 均小于0.03，表明模型性能良好。

圖9 WA-NARX 模型地下水位預(yù)測與觀測對比Fig.9 Comparison of groundwater level prediction based on WA-NARX and observation

表3 地下水位預(yù)測時WA-NARX 模型的均方根誤差（RMSE）和擬合度（R2）統(tǒng)計表Table 3 The WA-NARX model performance metrics of RMSE and R2 for groundwater level prediction

3.2.2 電導(dǎo)率模擬結(jié)果分析

BH1–BH6 6 口監(jiān)測井電導(dǎo)率的預(yù)測結(jié)果散點圖和時間序列圖如圖10 所示，從混合模型模擬結(jié)果可以看出：混合模型的訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)基本集中分布在一條直線上，電導(dǎo)率的變化趨勢與觀測電導(dǎo)率的變化趨勢一致，在訓(xùn)練和測試階段模型的誤差都相對較低。因此WA-NARX 神經(jīng)網(wǎng)絡(luò)混合模型適用于地下水電導(dǎo)率的模擬預(yù)測。表4 中所有監(jiān)測井在訓(xùn)練和測試階段的R2均大于0.98，RMSE 均小于0.009，說明該模型在進(jìn)行電導(dǎo)率的模擬預(yù)測時也表現(xiàn)出良好的性能。

表4 電導(dǎo)率預(yù)測時WA-NARX 模型的均方根誤差（RMSE）和擬合度（R2）統(tǒng)計Table 4 The WA-NARX model performance metrics of RMSE and R2 for electrical conductivity prediction

圖10 WA-NARX 模型電導(dǎo)率預(yù)測的性能圖Fig.10 Performance of WA-NARX model for electrical conductivity prediction

3.3 蒸發(fā)和排灌用電量對混合模型預(yù)測的影響

3.3.1 蒸發(fā)對水位預(yù)測的影響

在對比不同輸入?yún)?shù)對地下水位預(yù)測結(jié)果的影響時，在模型訓(xùn)練完成后，關(guān)閉WA-NARX 網(wǎng)絡(luò)的回路，輸出反饋進(jìn)行多步預(yù)測。以60 h 為時間步長進(jìn)行多步預(yù)測。表5 為僅有降雨、潮汐參數(shù)的地下水位預(yù)測結(jié)果與新增蒸發(fā)參數(shù)后（降雨、潮汐和蒸發(fā)）的地下水位預(yù)測結(jié)果的對比。由預(yù)測誤差可知，輸入?yún)?shù)僅為降雨和潮汐時，模型的預(yù)測精度較高；當(dāng)在此基礎(chǔ)上新增蒸發(fā)作為輸入?yún)?shù)時，模型預(yù)測結(jié)果發(fā)生了變化。與無新增參數(shù)的結(jié)果相比，監(jiān)測井BH1、BH3 和BH6 的預(yù)測精度有所提高，其余監(jiān)測井的預(yù)測精度有所降低。經(jīng)過整體對比發(fā)現(xiàn)，只有降雨和潮汐作為模型的輸入?yún)?shù)時，模型的預(yù)測效果略高于新增蒸發(fā)參數(shù)后的預(yù)測效果。說明在該地區(qū)進(jìn)行地下水位預(yù)測時，降雨和潮汐對地下水位變化的影響更為顯著。

圖11 是新增蒸發(fā)參數(shù)后得到的預(yù)測水位與實測水位的散點圖。圖中實線是觀測水位數(shù)據(jù)的連線，把降雨和潮汐作為WA-NARX 模型的輸入?yún)?shù)，預(yù)測60 h 時間步長的地下水位時，模型的預(yù)測值基本分布在實測的水位線上，預(yù)測效果良好。新增蒸發(fā)參數(shù)后，監(jiān)測井BH1、BH3 和BH6 的預(yù)測散點值更集中于實測水位線，表現(xiàn)出更好的預(yù)測效果，而其余觀測井的預(yù)測效果相比之前有所下降，這與表5 的分析結(jié)果一致。

圖11 新增蒸發(fā)參數(shù)的預(yù)測和實測地下水位散點圖Fig.11 Scatter diagrams of predicted and observed groundwater level with newly added evaporation parameters

表5 新增蒸發(fā)參數(shù)的地下水位預(yù)測誤差及網(wǎng)絡(luò)性能統(tǒng)計表Table 5 The prediction error of groundwater level and network performance with newly added evaporation parameters

3.3.2 排灌用電量對水位預(yù)測的影響

農(nóng)業(yè)排灌用電量與農(nóng)業(yè)抽水灌溉有關(guān)，灌溉則與地下水位密切相關(guān)，所以把排灌用電量作為一組參數(shù)，研究其對地下水位預(yù)測的影響。表6 為新增排灌用電量參數(shù)與只有降雨和潮汐作為參數(shù)時地下水位預(yù)測結(jié)果的對比。由預(yù)測誤差可知，只有降雨和潮汐作為輸入?yún)?shù)時模型有較高的預(yù)測精度；當(dāng)新增排灌用電量作為輸入?yún)?shù)后，模型的預(yù)測精度同之前相差甚小。其中監(jiān)測井BH5、BH6 的預(yù)測效果有細(xì)微的提高，其余監(jiān)測井的預(yù)測效果有所降低。說明在研究區(qū)內(nèi)農(nóng)業(yè)排灌用電量數(shù)據(jù)可以作為預(yù)測地下水位的一項參數(shù)，在使用時要結(jié)合研究區(qū)監(jiān)測井的實際情況，在農(nóng)業(yè)抽水灌溉量較大處的影響會更顯著。同時從電量時間序列可以看出，因統(tǒng)計關(guān)系，每個月的灌溉用電量基本為定值，而地下水位的預(yù)測周期是小時變化，所以，用電量的觀測與地下水預(yù)測頻率不完全一致，對模擬結(jié)果也有較大影響。如果用電量的統(tǒng)計頻率與預(yù)測相符，模型的預(yù)測精度應(yīng)有所提高。

表6 新增用電量參數(shù)的地下水位預(yù)測誤差及網(wǎng)絡(luò)性能統(tǒng)計Table 6 The prediction error of groundwater level and network performance with newly added power consumption parameters

圖12 是新增排灌用電量參數(shù)得到的預(yù)測水位與實測水位的散點圖。從圖中可以看出，每口井對新增排灌用電量參數(shù)后的預(yù)測效果響應(yīng)并不盡相同。BH1 與BH6 的預(yù)測結(jié)果相對近似；而BH2、BH3 及BH4 的預(yù)測水位分布及表6 中RMSE 共同顯示其預(yù)測值較只應(yīng)用降水與潮汐時有所下降。對于BH5 來說，加入用電量后預(yù)測結(jié)果整體向上偏移，精度有一定提高。對比結(jié)果說明井位地理位置分布及其區(qū)域地質(zhì)條件對用電量的響應(yīng)有較大關(guān)聯(lián)。

圖12 新增用電量參數(shù)的預(yù)測和實測地下水位散點圖Fig.12 Scatter diagrams of predicted and observed groundwater level with newly added power consumption parameters

3.3.3 蒸發(fā)對電導(dǎo)率預(yù)測的影響

圖13 是新增蒸發(fā)參數(shù)后得到的預(yù)測電導(dǎo)率與觀測電導(dǎo)率的曲線圖，表7 為不同輸入?yún)?shù)地下水電導(dǎo)率預(yù)測結(jié)果的對比統(tǒng)計表。由分析結(jié)果可知，在降雨、潮汐基礎(chǔ)上增加蒸發(fā)數(shù)據(jù)作為輸入?yún)?shù)時，無論在變化平緩區(qū)還是劇烈區(qū)，監(jiān)測井中電導(dǎo)率的預(yù)測精度都有顯著提高。特別對監(jiān)測井BH5，未加入蒸發(fā)數(shù)據(jù)時預(yù)測值略有偏高，加入蒸發(fā)對該井有明顯的校正作用。預(yù)測結(jié)果說明在該地區(qū)預(yù)測地下水電導(dǎo)率時，蒸發(fā)是重要影響因素。

圖13 新增蒸發(fā)參數(shù)的預(yù)測和實測電導(dǎo)率曲線圖Fig.13 Curve diagrams of predicted and observed electrical conductivity with newly added evaporation parameters

3.3.4 排灌用電量對電導(dǎo)率的影響

為了對比不同輸入?yún)?shù)對地下水電導(dǎo)率預(yù)測結(jié)果的影響，關(guān)閉WA-NARX 網(wǎng)絡(luò)的回路，輸出反饋進(jìn)行多步預(yù)測。把收集到的農(nóng)業(yè)抽水排灌用電量數(shù)據(jù)作為一項參數(shù)，對比其對地下水電導(dǎo)率預(yù)測效果的影響。由表7 可知，只有降雨和潮汐作為輸入?yún)?shù)時，模型有較好預(yù)測效果；新增農(nóng)業(yè)排灌用電量參數(shù)，監(jiān)測井BH1、BH2 和BH5 的預(yù)測效果有所提高，其余監(jiān)測井的預(yù)測效果基本保持一致，總體上新增排灌用電量參數(shù)后混合模型的預(yù)測效果有所提升。

表7 不同輸入?yún)?shù)電導(dǎo)率預(yù)測誤差及網(wǎng)絡(luò)性能評價統(tǒng)計Table 7 The prediction error of electrical conductivity and network performance with different input parameters

圖14 是新增排灌用電量參數(shù)得到的預(yù)測電導(dǎo)率與觀測電導(dǎo)率的曲線圖。預(yù)測值的變化趨勢和實測電導(dǎo)率的變化趨勢基本吻合，表現(xiàn)出較好的預(yù)測效果；當(dāng)新增農(nóng)業(yè)排灌用電量參數(shù)時，監(jiān)測井BH1、BH2 和BH5 的預(yù)測曲線與觀測電導(dǎo)率曲線的變化趨勢更加貼合，與表7 的分析結(jié)果一致。

圖14 新增用電量參數(shù)的預(yù)測和實測電導(dǎo)率曲線圖Fig.14 Curve diagrams of predicted and observed electrical conductivity with newly added power consumption parameters

4 結(jié)論與展望

萊州灣南岸是受海水入侵影響嚴(yán)重地區(qū)。通過對沿海典型剖面含水層地下水位和電導(dǎo)率的變化分析可知，該區(qū)海水入侵季節(jié)性變化趨勢明顯。主要表現(xiàn)在地下水電導(dǎo)率與水位存在一定的負(fù)相關(guān)關(guān)系，在降雨補給充足的情況下，地下水水位有較為明顯的上升趨勢，地下水電導(dǎo)率降低，海水入侵程度減弱。但由于受地下水開采的影響，水位和電導(dǎo)率均存在突變現(xiàn)象，水位降低打破原有的咸淡水平衡，使海水入侵加劇。同時，本文建立了基于小波分析的非線性自回歸神經(jīng)網(wǎng)絡(luò)混合模型（WA-NARX），分別對地下水位和地下水中的電導(dǎo)率進(jìn)行預(yù)測，通過增加模型的輸入?yún)?shù)，對比分析不同輸入?yún)?shù)對模型預(yù)測效果的影響。研究發(fā)現(xiàn)，地下水位和電導(dǎo)率的變化與降雨和潮汐的變化密切相關(guān)。當(dāng)降雨和潮汐作為模型的輸入?yún)?shù)預(yù)測地下水位與電導(dǎo)率時，模型預(yù)測效果良好。蒸發(fā)參數(shù)和農(nóng)業(yè)排灌用電量對地下水位與電導(dǎo)率也有一定影響。模型預(yù)測依然具有較高精度，個別觀測點精度有所提高。通過上述研究表明，WA-NARX 模型可以作為海水入侵預(yù)測的有效工具，不僅可以為水資源的高效利用提供參考，還可以掌握海岸帶地下水鹽度的變化，為海水入侵的監(jiān)測預(yù)測預(yù)警、保護(hù)地下淡水資源及海岸帶環(huán)境保護(hù)提供科學(xué)依據(jù)。