聯(lián)合EEMD 與BP 神經(jīng)網(wǎng)絡(luò)的灌區(qū)水源情勢(shì)預(yù)測(cè)研究

李雯晴 ，劉 招, ，王麗霞，李 強(qiáng)，吳小宏

（1.長(zhǎng)安大學(xué) 水利與環(huán)境學(xué)院，西安 710054；2.長(zhǎng)安大學(xué) 旱區(qū)地下水文與生態(tài)效應(yīng) 教育部重點(diǎn)實(shí)驗(yàn)室，西安 710054；3.長(zhǎng)安大學(xué) 水與發(fā)展研究院，西安 710054； 4.長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，西安 710054；5.涇惠渠灌溉管理局，陜西 三原 713800）

0 引 言

水是農(nóng)業(yè)的命脈。水源情勢(shì)對(duì)依賴引水灌溉的大型灌區(qū)而言，是直接關(guān)系到農(nóng)業(yè)生產(chǎn)和糧食安全的重大問(wèn)題，歷來(lái)受到地區(qū)和國(guó)家的高度重視。人類社會(huì)進(jìn)入21 世紀(jì)以來(lái)，氣候條件變化和人類活動(dòng)加劇，對(duì)原有的水文規(guī)律帶來(lái)了巨大沖擊，河川徑流銳減、洪水頻發(fā)、區(qū)域地下水大面積疏干[1]。【研究意義】水資源規(guī)律的改變給工農(nóng)業(yè)用水及區(qū)域水資源管理帶來(lái)了諸多問(wèn)題，尤其對(duì)干旱半干旱地區(qū)的廣大灌區(qū)而言，深入探討水源情勢(shì)與灌溉需水規(guī)律，確保糧食生產(chǎn)安全是一項(xiàng)緊迫而重要工作。

水源形勢(shì)的改變對(duì)區(qū)域水資源利用帶來(lái)新的問(wèn)題和挑戰(zhàn)?！狙芯窟M(jìn)展】近年來(lái)，許多專家學(xué)者對(duì)區(qū)域水源形勢(shì)與水資源可利用量進(jìn)行了大量研究[2-5]，包括對(duì)水沙變化規(guī)律[6]、流量與水沙過(guò)程響應(yīng)關(guān)系[7]、水沙規(guī)律預(yù)測(cè)[8]、地表水資源可利用量潛力[9]等問(wèn)題的研究。采用了滑動(dòng)平均[10]、Mann-Kendall 秩次相關(guān)檢驗(yàn)法[11]、BP 神經(jīng)網(wǎng)絡(luò)[12]、Copula 函數(shù)[13]等方法?！厩腥朦c(diǎn)】在前人對(duì)水源情勢(shì)的研究中，多集中于水沙演變規(guī)律這類淺層意義上的研究；對(duì)水源情勢(shì)的預(yù)測(cè)多采用BP 神經(jīng)網(wǎng)絡(luò)這種單一的方法，與其他方法如EEMD、小波分析等結(jié)合研究較少；在預(yù)測(cè)未來(lái)水源情勢(shì)并同時(shí)分析水資源可利用量方面研究較少。【擬解決的關(guān)鍵問(wèn)題】本研究嘗試聯(lián)合總體經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）與BP（Back Propagation）神經(jīng)網(wǎng)絡(luò)對(duì)灌區(qū)地表水源及灌溉需水進(jìn)行預(yù)測(cè)，分析未來(lái)情景下渠首引水和灌溉需水量間的可能關(guān)系，為緩解目前灌區(qū)水資源矛盾、高效利用和管理水資源、保障農(nóng)業(yè)生產(chǎn)安全提供技術(shù)支撐及參考。

1 材料與方法

1.1 研究區(qū)概況

1.1.1 灌區(qū)概況

涇惠渠灌區(qū)位于陜西省關(guān)中平原中部，其前身即為我國(guó)古代三大水利工程之一的鄭國(guó)渠。涇惠渠灌區(qū)灌溉面積9.69 萬(wàn)hm2，轄西安、咸陽(yáng)及渭南地區(qū)的6 縣48 個(gè)鄉(xiāng)鎮(zhèn)，總?cè)丝?18 萬(wàn)人。糧食作物以小麥、玉米為主，經(jīng)濟(jì)作物以蔬菜、果樹(shù)為主，是陜西省重要的農(nóng)產(chǎn)品生產(chǎn)區(qū)[14]。據(jù)涇惠渠渠首——涇河張家山站實(shí)測(cè)年徑流量數(shù)據(jù)顯示，灌區(qū)地表水源的涇河徑流泥沙呈現(xiàn)出明顯的遞減趨勢(shì)[15]，多年平均實(shí)測(cè)徑流量從20 世紀(jì)60 年代的17.15 億m3下降到20 世紀(jì)70、80 年代的13 億m3左右，直至2000 年來(lái)的10 億m3左右，徑流衰減迅速。與此同時(shí)，由于地下水的過(guò)量開(kāi)采，地下水位埋深也由1982 年的4.74 m 下降為2016 年的20 m，局部地下水埋深甚至降至 50 m 以下[16]。

1.1.2 灌區(qū)渠首涇河張家山站水文情勢(shì)

張家山站歷史徑流數(shù)據(jù)序列如圖1 所示。據(jù)涇河張家山站1970—2017 年年徑流量序列分析表明，灌區(qū)渠首年徑流量波動(dòng)性較大，多年變化幅度為3.22 億～20.87 億m3，多年均值為11.85 億m3。年徑流量在20 世紀(jì)90 年代前相對(duì)平穩(wěn)，但從20 世紀(jì)90 年代中期起出現(xiàn)明顯下降趨勢(shì)。進(jìn)行MK 趨勢(shì)檢驗(yàn)得其Z 值為-1.52，表明年徑流量為下降趨勢(shì)。

圖1 灌區(qū)渠首涇河張家山站歷史徑流泥沙數(shù)據(jù)序列 Fig.1 Historical runoff and sediment data sequence of Zhangjiashan station, canal head, irrigation area

1.1.3 灌區(qū)取用水現(xiàn)狀

灌區(qū)渠首引水主要用于灌溉、發(fā)電和補(bǔ)充西郊水庫(kù)。統(tǒng)計(jì)涇惠渠渠首多年實(shí)際引水及主要用水構(gòu)成如圖2 所示。由圖2 可見(jiàn)，20 世紀(jì)70 年代，渠首總引水量約為4.60 億m3，20 世紀(jì)80 年代后開(kāi)始出現(xiàn)下降趨勢(shì)，2005 年后維持在3.6 億m3左右。地表水有效灌溉量均值亦呈下降趨勢(shì)，從20 世紀(jì)80 年代的2.50 億m3下降到2010 年的1.44 億m3，約占總引水量的40%，多年平均值為1.80 億m3，約占據(jù)總引水量49.18%。以1985 和2015 年為例，統(tǒng)計(jì)不同時(shí)期涇惠渠灌區(qū)地表水資源利用情況如表1 所示。

發(fā)電用水量多年平均為1.29 億m3，占總引水量的35.24%，其余水量充蓄西郊水庫(kù)，若仍有盈余時(shí)將退泄水量排至石川河。西郊水庫(kù)為典型的灌區(qū)內(nèi)水庫(kù)，其除調(diào)蓄灌溉水資源外還供給三原縣等地的工業(yè)、生活用水，有效庫(kù)容0.20 億m3。

因來(lái)水與作物需水規(guī)律不盡一致，加之水價(jià)及農(nóng)民主管等因素，致使當(dāng)前灌區(qū)地表水有效灌溉水量占總引水量的份額較低，引取的相當(dāng)部分水量并非直接用于灌溉。一般年份，灌區(qū)灌溉總需水量均在3 億m3以上，這就意味著灌區(qū)消耗較多地下水的同時(shí)地表水未能得以充分利用，這種情形十分不利于灌區(qū)水資源利用。同時(shí)，若涇河來(lái)水和地表水可引水量進(jìn)一步降低，必然對(duì)維持良性的灌區(qū)水文生態(tài)帶來(lái)較大沖擊。

圖2 涇惠渠渠首多年實(shí)際引水及主要用水構(gòu)成 Fig.2 Actual water diversion and main water use composition of Jinghui canal head for many years

表1 不同時(shí)期涇惠渠灌區(qū)地表水資源利用情況 Table 1 Utilization of surface water resources in Jinghui Irrigation Area in different periods

1.2 基于EEMD-BP 的預(yù)測(cè)方法

為了給灌區(qū)未來(lái)水資源管理提供支撐，研究聯(lián)合經(jīng)驗(yàn)?zāi)B(tài)分解及神經(jīng)網(wǎng)絡(luò)來(lái)預(yù)測(cè)涇河未來(lái)10 a 來(lái)水。河川徑流受氣象、水文、水力等多方面因素影響，具有明顯的隨機(jī)性和復(fù)雜性。徑流序列是一個(gè)典型的非線性過(guò)程，預(yù)測(cè)難度較大[17]。EEMD 是一種自適應(yīng)的數(shù)據(jù)處理或挖掘方法，非常適合非線性、非平穩(wěn)時(shí)間序列，可有效地提取序列中變化特征、趨勢(shì)等隱含信息；BP 神經(jīng)網(wǎng)絡(luò)是基于人工神經(jīng)網(wǎng)絡(luò)所提出的誤差反向傳播算法，其對(duì)因果關(guān)系復(fù)雜、變化規(guī)律模糊的問(wèn)題有出色的分析及遞推能力[18]。將二者結(jié)合可有效發(fā)揮各自優(yōu)勢(shì)，大大增強(qiáng)對(duì)非平穩(wěn)時(shí)間序列的預(yù)測(cè)能力[19]。EEMD 為解決模態(tài)混淆問(wèn)題而提出，通過(guò)加入白噪音再進(jìn)行分解，利用白噪音均值為0 的特性進(jìn)行篩選，即可消除采用EMD方法時(shí)IMF（Intrinsic Mode Function）里夾雜2 個(gè)不同時(shí)間范圍訊號(hào)的現(xiàn)象。EEMD 分解出的IMF 減少模態(tài)混合現(xiàn)象，使不同時(shí)間尺度的訊號(hào)正確分布在對(duì)應(yīng)的IMF 中[20-21]。

EEMD-BP 時(shí)間序列預(yù)測(cè)模型建立思路如下：

1）利用EEMD 方法對(duì)水文時(shí)間序列進(jìn)行分解，得到不同頻率的IMF，并將其進(jìn)行歸一化處理。

2）對(duì)BP 神經(jīng)網(wǎng)絡(luò)的隱含層個(gè)數(shù)、訓(xùn)練速率、訓(xùn)練次數(shù)、控制誤差等參數(shù)進(jìn)行初始化設(shè)置，1973—2012 年為訓(xùn)練期，2013—2017 年實(shí)測(cè)數(shù)據(jù)為驗(yàn)證數(shù)據(jù)。將經(jīng)EEMD 分解處理后的輸入層與輸出層數(shù)據(jù)導(dǎo)入并進(jìn)行訓(xùn)練。

3）訓(xùn)練完成后，用訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)對(duì)灌區(qū)渠首2018—2027 年徑流泥沙量進(jìn)行預(yù)測(cè)，得到不同頻率預(yù)測(cè)值。

4）將不同頻率的預(yù)測(cè)值進(jìn)行疊加重構(gòu)，最終得到預(yù)測(cè)結(jié)果。

涇河作為知名的多泥沙河流，其高含沙水流及汛期洪水無(wú)法直接利用，在對(duì)灌區(qū)渠首水源形勢(shì)進(jìn)行預(yù)測(cè)和分析時(shí)，這部分不能作為可用水資源來(lái)對(duì)待。因此，考慮到含沙率影響，本研究中，不僅要預(yù)測(cè)徑流量，同時(shí)還應(yīng)同步預(yù)測(cè)河道逐年輸沙量。

1.3 模型的訓(xùn)練及驗(yàn)證

一般而言，同條河流同一站點(diǎn)的徑流及泥沙序列自相關(guān)性較強(qiáng)，除此之外，含沙量同時(shí)與徑流量也有較強(qiáng)的相依關(guān)系，即所謂“水大沙大”。研究基于上述EEMD-BP時(shí)間序列預(yù)測(cè)模型，對(duì)歷史數(shù)據(jù)系列進(jìn)行EEMD分解，再輸入BP神經(jīng)網(wǎng)絡(luò)對(duì)分解序列進(jìn)行預(yù)測(cè)，然后重構(gòu)，對(duì)模型反復(fù)訓(xùn)練，最終將EEMD-BP神經(jīng)網(wǎng)絡(luò)與BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值進(jìn)行對(duì)比，結(jié)果見(jiàn)圖3。EEMD-BP神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)和主要參數(shù)見(jiàn)表2。由圖3可以直觀看出，聯(lián)合EEMD的BP神經(jīng)網(wǎng)絡(luò)與實(shí)測(cè)值擬合效果比BP神經(jīng)網(wǎng)絡(luò)擬合效果好，能夠較為準(zhǔn)確地反映真實(shí)值。

表2 EEMD-BP神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)與參數(shù) Table 2 EEMD-BP neural network model structure and parameters

圖3 徑流泥沙預(yù)測(cè)結(jié)果 Fig.3 Prediction results of runoff and sediment

EEMD-BP 模型驗(yàn)證期年徑流及泥沙的標(biāo)準(zhǔn)均方差NMSE 分別為0.265 3 和0.307 8；平均絕對(duì)誤差MAE 分別為0.9 億m3和0.13 億t，相較多年平均值，誤差在可接受范圍內(nèi)；預(yù)測(cè)結(jié)果的方向變差對(duì)稱值DVS 均為100%，表明預(yù)測(cè)系列的變化趨勢(shì)與實(shí)測(cè)系列完全一致；決定系數(shù)R2分別為0.980 1 和0.986 2。故模型預(yù)測(cè)效果良好。

2 結(jié)果與分析

2.1 渠首徑流及輸沙量預(yù)測(cè)

基于1970—2017 年徑流泥沙數(shù)據(jù)，利用建立好的EEMD-BP 預(yù)測(cè)模型，對(duì)2018—2027 年的灌區(qū)渠首張家山站徑流泥沙進(jìn)行預(yù)測(cè)，結(jié)果見(jiàn)圖4。由圖4可知，2018—2027 年的模擬期內(nèi)，水沙變化較為同步，趨勢(shì)均為短暫的上升后再下降。年均徑流量約為11.86 億m3，最大為16.25 億m3，最小為8.14 億m3；輸沙量多年均值為1.00 億t，最大2.58 億t，最小0.29 億t。與1996—2017 年相比，平均徑流量略有上升，回升值為1.65 億m3；平均輸沙量略有下降，年均下降了0.17 億t。

對(duì)2008—2017 年實(shí)測(cè)數(shù)據(jù)和此次預(yù)測(cè)的2018—2027 年數(shù)據(jù)的特征參數(shù)進(jìn)行比對(duì)，若特征參數(shù)較為吻合，則說(shuō)明預(yù)測(cè)結(jié)果基本可信。檢驗(yàn)結(jié)果見(jiàn)表3。

圖4 涇河張家山站未來(lái)10 年徑流量輸沙量序列預(yù)測(cè) Fig.4 Prediction of runoff and sediment transport series of Zhangjiashan station in the next 10 years

表3 預(yù)測(cè)徑流及泥沙系列與實(shí)測(cè)系列特征對(duì)比 Table 3 Characteristics comparison between predicted runoff and sediment series and measured series

由表3 可知，徑流及輸沙量預(yù)測(cè)結(jié)果與歷史數(shù)據(jù)系列的特征值總體上并無(wú)明顯差異。且由于徑流泥沙受氣候變化和人類活動(dòng)影響，水文時(shí)間序列特征本身也會(huì)有一定波動(dòng)。同時(shí)考慮到?jīng)芎又猩嫌瘟饔騼?nèi)水土保持和退耕還林還草工作將持續(xù)展開(kāi)，流域水保及人類活動(dòng)程度基本穩(wěn)定，所以總體認(rèn)為預(yù)測(cè)結(jié)果基本合理可信。

2.2 渠首地表水資源可利用量分析

涇河流經(jīng)黃土高原，泥沙量很高。對(duì)灌區(qū)而言，這些高泥沙水流及汛期洪水一般無(wú)法取用，故其均為不可利用水量，在水資源分析時(shí)應(yīng)予以剔除。在灌區(qū)實(shí)際生產(chǎn)中，考慮到渠道沖淤等限制，實(shí)際灌溉引水控制含沙量不應(yīng)超7%，在持續(xù)干旱等特殊情況下可有所放寬，但含沙量仍不可超過(guò)10%。此外，在非汛期需要優(yōu)先保證河道生態(tài)基流，在計(jì)算灌區(qū)可利用水資源量時(shí)，這部分亦應(yīng)予以扣除。采用式（1）來(lái)估算渠首水資源可利用量[9]。計(jì)算式為：

式中：Wa為水資源可利用量；Wt為年徑流量；We為非汛期河道最小生態(tài)需水量；Wu為汛期不可利用水量。單位均為108m3。

計(jì)算非汛期河道最小生態(tài)需水量采用常用的基于水文學(xué)的Tenant 方法，在枯水季，即10 月—次年5 月，最小生態(tài)流量按多年平均流量的15%來(lái)計(jì)算[9]，求得We為0.734 億m3；汛期不可利用水量包括2 部分，即汛期洪水和高泥沙徑流。而由于水文不確定性，這2 部分在各年份情況差異較大，較難進(jìn)行準(zhǔn)確計(jì)算。事實(shí)上，如前述分析，涇河河道洪水與泥沙含量相關(guān)關(guān)系較為密切，高泥沙水流多由上游暴雨洪水所致，因此這二者可以合一考慮[22]。本文統(tǒng)一采用汛期水流含沙量這一指標(biāo)來(lái)推算汛期不可利用水量。為此，需先求得各年份汛期徑流量，基于年徑流及汛期徑流歷史數(shù)據(jù)，進(jìn)行回歸分析，最終求得相關(guān)關(guān)系較好的二次多項(xiàng)式回歸方程，其相關(guān)系數(shù)r 為0.89：

式中：Wf為汛期（每年6—9 月）徑流量。

因涇河汛期來(lái)沙占據(jù)全年來(lái)沙量的97%以上，故本文不考慮非汛期輸沙。根據(jù)汛期徑流量及預(yù)測(cè)的沙量即可求得汛期水流平均含沙量Cs，如表4 所示。再依據(jù)歷史數(shù)據(jù)，統(tǒng)計(jì)汛期棄水與平均含沙量關(guān)系，通過(guò)回歸分析求得二者相關(guān)關(guān)系較好的關(guān)系式，相關(guān)系數(shù)r 為0.86：

求得We和Wu后，即可利用式（1）計(jì)算得到渠首逐年水資源可利用量，結(jié)果如表4 所示。

表4 基于預(yù)測(cè)的涇惠渠灌區(qū)渠首地表水資源可利用量分析 Table 4 Analysis of available quantity of surface water resources at the head of canal in Jinghui irrigation area based on prediction

3 討 論

EEMD 非常適合非線性、非平穩(wěn)時(shí)間序列，可有效地提取序列中變化特征、趨勢(shì)等隱含信息，BP神經(jīng)網(wǎng)絡(luò)對(duì)因果關(guān)系復(fù)雜、變化規(guī)律模糊的問(wèn)題有出色的分析及遞推能力，將二者結(jié)合可有效發(fā)揮各自優(yōu)勢(shì)，大大增強(qiáng)對(duì)非平穩(wěn)時(shí)間序列的預(yù)測(cè)能力。未來(lái)涇惠渠灌區(qū)渠首可利用水資源量相對(duì)比較平穩(wěn)，平均值約為7.79 億m3，這較目前略有上升。主要是由于涇河泥沙呈下降趨勢(shì)，來(lái)水變清，在一定程度上提高了水資源可利用率。事實(shí)上，近年來(lái)涇惠渠灌區(qū)渠首年均引水量約為3.6 億m3，從預(yù)測(cè)結(jié)果看來(lái)，若灌區(qū)引水需求無(wú)大的變化，基本可以滿足未來(lái)幾年的地表水資源需求。2019—2027 年最大可利用量為2021 年的11.25 億m3，最小為2023 年的4億m3左右，在2023 年應(yīng)當(dāng)注意灌區(qū)水資源的調(diào)蓄和有效利用，謹(jǐn)防發(fā)生缺水減產(chǎn)問(wèn)題。此外輸沙量的預(yù)測(cè)值對(duì)水資源可利用量影響較為顯著。盡管年預(yù)測(cè)值最大的2019 年徑流總量達(dá)16.26 億m3，但其泥沙預(yù)測(cè)值也較高，致使不可利用水量達(dá)7.3 億m3，而可利用水資源量為8.22 億m3，尚不及徑流量相對(duì)較小的2020、2021 和2025 年。

由于渠首及引水工程有確定的設(shè)計(jì)標(biāo)準(zhǔn)，上述計(jì)算的可利用水資源量并非完全可為工程所用。涇惠渠灌區(qū)渠首的設(shè)計(jì)引水能力為50 m3/s[14]，若不考慮水庫(kù)調(diào)蓄，當(dāng)河道流量超出這一值時(shí)無(wú)法被工程利用。因此實(shí)際工程可利用量并非河道水資源可利用量，只有不超出引水能力的河道水量才是渠首工程的可引水量。

4 結(jié) 論

1）建立的EEMD-BP 模型預(yù)測(cè)年徑流及泥沙效果良好。二者標(biāo)準(zhǔn)均方差NMSE 分別為0.265 3 和0.307 8，平均絕對(duì)誤差MAE 分別為0.9 億m3和0.13 億t，誤差在可接受范圍內(nèi)，預(yù)測(cè)系列的變化趨勢(shì)與實(shí)測(cè)系列完全一致，決定系數(shù)R2分別為0.980 1 和0.986 2；

2）灌區(qū)渠首未來(lái)2019—2027 年均徑流量約為11.87 億m3，較1996—2017 年略有上升，年均輸沙量約為1.00 億t，較1996—2017 年略有下降，預(yù)測(cè)結(jié)果延續(xù)了涇河水沙多年變化的大致趨勢(shì)；

3）灌區(qū)渠首2019—2027 年可利用水資源量相對(duì)比較平穩(wěn)，平均水資源可利用量約為7.79 億m3，這較1996—2017 年可利用量略有上升但相差不大，在灌區(qū)引水需求無(wú)大的變化條件下，基本可滿足未來(lái)引水灌溉需求，但應(yīng)注意一些干旱年份水資源的調(diào)蓄和有效利用，謹(jǐn)防發(fā)生缺水減產(chǎn)問(wèn)題。