岳斌，牛最榮，曹志宏，王啟優(yōu)，朱詠

（1.甘肅省水文站，甘肅 蘭州 730030；2.甘肅農(nóng)業(yè)大學(xué)水利水電工程學(xué)院，甘肅，蘭州 730070）

由于人類活動改變了流域下墊面條件，導(dǎo)致入滲、徑流、蒸發(fā)等水平衡要素發(fā)生一定的變化，從而造成徑流的減少（或增加）。按第三次《全國水資源調(diào)查評價技術(shù)細(xì)則》要求，對河川年徑流系列需進(jìn)行還原和一致性處理，推薦使用降水徑流關(guān)系法。該方法本質(zhì)是建立降水與徑流之間的一元函數(shù)關(guān)系，方程經(jīng)最小二乘法概化，與每一個數(shù)據(jù)點(diǎn)均有一定誤差，計算修正值按系數(shù)比例縮小，修正結(jié)果誤差較大。國內(nèi)相關(guān)學(xué)者對一致性修正也有諸多報道，韓瑞光等［1］研究得出海河流域人類活動加劇，下墊面發(fā)生較大變化，利用天然徑流一致性修正方法可較好地分析出人類活動的間接影響；劉俊等［2］應(yīng)用降水徑流關(guān)系法修正唐山市地表水資源，結(jié)果表明修正后的結(jié)果更科學(xué)；李芳等［3］對渭河干流代表水文站進(jìn)行一致性分析，將過去的年徑流資料修正為近期下墊面條件下的徑流數(shù)據(jù)，以保證水文資料具有一致性；張波等［4］提出了基于典型解集模型的非一致性年徑流過程設(shè)計方法，該方法充分考慮變化環(huán)境的影響，設(shè)計年徑流過程的結(jié)果更加真實、可靠。張利茹等［5］指出變化環(huán)境下的水文序列資料不能直接用于水文模型參數(shù)的率定，應(yīng)首先對水文序列資料進(jìn)行突變檢測分析，然后再進(jìn)行還原分析計算。顧冉浩等［6］采用分項調(diào)查法對徑流還原計算，利用天然徑流構(gòu)建SWAT降雨徑流模型，模擬得到全研究區(qū)的天然徑流數(shù)據(jù)。以上文獻(xiàn)應(yīng)用不同方法對天然徑流進(jìn)行一致性修正，建立了不同的降雨徑流模型，但以上模型均僅考慮了降水因素，而未分析氣溫、蒸發(fā)等其他影響因素。

本文基于LM-BP神經(jīng)網(wǎng)絡(luò)法對天然徑流進(jìn)行一致性修正，模型輸入變量為降水、氣溫等5種影響因子的10種不同組合，引入誤差函數(shù)MAPE作為最優(yōu)模型的判定依據(jù)，該方法充分考慮了與天然徑流相關(guān)的不同影響因子，并利用神經(jīng)網(wǎng)絡(luò)任意非線性逼近的特征，建立天然徑流量預(yù)測模型，推求天然徑流修正值。

1 研究區(qū)概況

武山水文站位于渭河干流中上游，1974年7月建站，集水面積8 080 km2，河長138 km。測站以上較大支流有秦祁河、咸河、榜沙河、漳河等；流域以干流為界南北差異較大，南部植被較好，有森林、草場分布，水量較充沛；北部屬黃土高原溝壑區(qū)，植被差，水量小，多數(shù)支溝非雨期干涸。該站洪水一般由暴雨形成，主要發(fā)生在5~10月，洪水歷時一般1~3 d；測驗河段順直，左右岸均有石砌河堤，沙卵石河床，復(fù)式斷面［7-8］。采用低水流速儀、高水浮標(biāo)法施測流量，臨時曲線法推流整編。武山水文站位置示意圖見圖1。

圖1 武山水文站位置示意圖Figure 1 Location of Wushan hydrological station

受全球氣候變化和人類活動的影響，近年來渭河流域的徑流時空分布發(fā)生了很大變化，水資源持續(xù)減少，用水矛盾突出，生態(tài)環(huán)境脆弱，這對水資源合理并持續(xù)利用帶來了極大的挑戰(zhàn)［7］。

2 數(shù)據(jù)來源及分析方法

2.1 數(shù)據(jù)來源及資料“三性”分析

采用1956~2016年武山水文站及半陰坡、慶坪、堯甸等雨量站的資料來源于水文年鑒，氣象站的資料來源于中國氣象數(shù)據(jù)網(wǎng)，部分站點(diǎn)資料系列不同步的經(jīng)鄰近站插補(bǔ)延長獲得。

應(yīng)用泰森多邊形法計算流域平均面降水量；由測站年資料算數(shù)平均法求得流域年平均氣溫、水面蒸發(fā)；采用分項調(diào)查、分月計算還原量，并由實測徑流量和還原量之和得出天然徑流量。所用資料均依據(jù)相應(yīng)規(guī)范要求測驗、整編，經(jīng)一致性、完整性、代表性檢查，數(shù)據(jù)完整、準(zhǔn)確、可靠。

2.2 分析方法

2.2.1 徑流趨勢與突變檢驗法 采用水文時間序列趨勢與突變分析系統(tǒng)對天然年徑流序列變化趨勢和突變檢驗，該軟件集成了目前多種較成熟的水文檢驗方法，包括曼-肯德爾、斯波曼秩次相關(guān)、線性回歸、雙累積曲線法趨勢性檢驗法；R/S分析、累計距平、有序聚類分析、里和海哈林法分析系列突變點(diǎn)檢驗法［9-11］。首先對天然年徑流序列的變化趨勢進(jìn)行檢驗，若變化趨勢顯著，再進(jìn)行突變檢驗，并以突變點(diǎn)年份將徑流序列分為變前和突變后兩段。發(fā)生突變的主要原因是兩個時段內(nèi)的產(chǎn)匯流條件發(fā)生顯著變化，若要保持天然徑流序列的一致性，需根據(jù)突變后產(chǎn)匯流條件修正突變前的天然徑流序列。

本文在分析天然徑流趨勢和突變檢驗的基礎(chǔ)上，以天然徑流的影響因子為輸入變量，應(yīng)用LMBP神經(jīng)網(wǎng)絡(luò)模型建立突變后天然徑流的預(yù)測模型，將該模型移植到突變前天然徑流的推求中，用突變前的影響因子推求天然徑流，達(dá)到一致性修正目的。

2.2.2 LM-BP神經(jīng)網(wǎng)絡(luò)一致性修正法

2.2.2.1 模型原理 LM-BP神經(jīng)網(wǎng)絡(luò)法是指網(wǎng)絡(luò)結(jié)構(gòu)采用誤差反向傳播的BP神經(jīng)網(wǎng)絡(luò)，求解誤差函數(shù)極小值的算法采用LM（Levenberg-Marquardt）算法。BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)具有多層感知器，通過對人類神經(jīng)元的功能進(jìn)行模擬、儲存及學(xué)習(xí)輸入和輸出數(shù)據(jù)，不需對變量的映射關(guān)系進(jìn)行描述，利用樣本數(shù)據(jù)建模，實現(xiàn)從輸入到輸出的任意非線性映射［12-14］。BP網(wǎng)絡(luò)的學(xué)習(xí)過程是一種誤差反向傳播的過程，同時通過修正各層神經(jīng)元的權(quán)值、閾值，調(diào)整輸出使誤差信號最小。根據(jù)Kolmogorov定理，具有一個隱含層的三層BP神經(jīng)網(wǎng)絡(luò)能在閉集上任意精度逼近非線性連續(xù)函數(shù)［15-16］。其結(jié)構(gòu)見圖2。

圖2 神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)Figure 2 Structure of neural network model

本文采用LM算法，其基本思想是使其每次迭代不再沿著單一的負(fù)梯度方向，而是允許誤差沿著惡化的方向進(jìn)行搜索，同時通過在最速梯度下降法和高斯-牛頓法之間自適應(yīng)調(diào)整來優(yōu)化網(wǎng)絡(luò)權(quán)值，使網(wǎng)絡(luò)誤差函數(shù)能夠有效收斂于該函數(shù)的最小點(diǎn)，提高了網(wǎng)絡(luò)的收斂速度和泛化能力［17-20］。LM算法基于避免計算修正速率中Hessian矩陣而設(shè)計，根據(jù)下式修正網(wǎng)絡(luò)權(quán)值：

式中：J為包含誤差性能函數(shù)對網(wǎng)絡(luò)權(quán)值一階導(dǎo)數(shù)的雅克比矩陣；I為單位矩陣；e為誤差向量；Y k為正向計算的網(wǎng)絡(luò)輸出向量；Tk為實際的輸出樣本向量；p為樣本個數(shù)；w為神經(jīng)網(wǎng)絡(luò)權(quán)值組成的向量；μ為自適應(yīng)調(diào)整系數(shù)，為試探性參數(shù)。在實際操作過程中，開始時μ取小值，如果求得的w能使誤差函數(shù)指標(biāo)E（w）降低，則該值取μ/β（β>1）；反之，該值取μ×β。通常設(shè)置μ初始值0.01，β取10。

LM算法更適用于訓(xùn)練權(quán)值和閾值數(shù)目不超過幾百個的神經(jīng)網(wǎng)絡(luò)［21］。根據(jù)本文待解問題的特征，可選用LM算法作為BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型的訓(xùn)練方法［22-23］。

2.2.2.2 模型設(shè)計 采用經(jīng)典的三層BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，即輸入層、隱含層和輸出層，隱含層神經(jīng)元節(jié)點(diǎn)數(shù)過多或過少都會影響預(yù)測模型的精度［24］。輸入層神經(jīng)元個數(shù)等于輸入變量的個數(shù)，隱含層神經(jīng)元個數(shù)q和輸入層神經(jīng)元個數(shù)M之間一般具有以下近似關(guān)系［25］：

模型的隱含層神經(jīng)元傳遞函數(shù)采用purelin，輸出層神經(jīng)元傳遞函數(shù)采用tansig［26］。神經(jīng)網(wǎng)絡(luò)的輸入和輸出值需限制在一定的范圍內(nèi)，使較大的輸入值仍落在神經(jīng)元轉(zhuǎn)化函數(shù)梯度大的地方，可加快網(wǎng)絡(luò)的訓(xùn)練速度，輸入輸出值一般在［-1，1］之間，本文采用mapminmax函數(shù)進(jìn)行歸一化處理。

LM-BP神經(jīng)網(wǎng)絡(luò)在開始訓(xùn)練前，各層的連接權(quán)值和閾值隨機(jī)初始化為［0，1］之間的值，因此神經(jīng)網(wǎng)絡(luò)每次訓(xùn)練得到的結(jié)果不一樣。本文采用相對誤差平均值MAPE指標(biāo)函數(shù)，設(shè)定每次迭代次數(shù)為200，總共運(yùn)行100次，指標(biāo)函數(shù)最優(yōu)的結(jié)果作為天然徑流一致性修正模型的參數(shù)。MAPE指標(biāo)函數(shù)：

式中：yk為實際天然徑流值；ˉyk為預(yù)測天然徑流值。

天然徑流預(yù)測模型的輸入為與徑流相關(guān)的多個影響因子，輸出為天然徑流量。應(yīng)用突變點(diǎn)后的系列值進(jìn)行模型訓(xùn)練，根據(jù)誤差函數(shù)MAPE值最小原則確定最優(yōu)模型，將該模型應(yīng)用到突變點(diǎn)前，以突變點(diǎn)前的影響因子為輸入，求出修正后的天然徑流量。

3 結(jié)果與分析

3.1 天然徑流變化趨勢及突變點(diǎn)

基于武山水文站1956~2016年天然徑流序列，采用3種趨勢檢驗法和5種突變點(diǎn)檢驗法，置信水平α=0.05，檢驗結(jié)果見表1。綜合分析檢驗結(jié)果可得，天然年徑流序列存在明顯下降趨勢，傾向率為-0.758億m3/10·a；突變點(diǎn)為1993年，突變顯著。以1993年為分界點(diǎn)，突變前徑流序列為1956~1993年，突變后為1994~2016年。

表1 武山站天然年徑流系列趨勢及跳躍分析結(jié)果Tab1 Trend and jump analysis results of natural annual runoff series at Wushan station

3.2 天然徑流一致性修正

由3.1可知，受人類活動影響，流域內(nèi)下墊面發(fā)生變化，徑流系列于1993年發(fā)生突變，LM-BP神經(jīng)網(wǎng)絡(luò)模型以1994~2016年數(shù)據(jù)進(jìn)行訓(xùn)練，求出模型最優(yōu)參數(shù)，再用1956~1993年的影響因子推求出天然徑流量。

選用流域內(nèi)氣溫（T）、降水量（P）、水面蒸發(fā)（E）、太陽黑子相對數(shù)（R）及武山站6~9月降水量（P6-9）影響因子作為模型的輸入量。按照影響天然年徑流量的因子數(shù)量和類型，共建立不同輸入變量組合形式的10種組合模型，對10種模型進(jìn)行訓(xùn)練，并對1956~1993年天然徑流一致性修正，可得不同組合影響因子情況下修正結(jié)果（表2）。

由表2可得，修正后天然年徑流量范圍3.786～4.857億m3，其中BP4-1組合模型的MAPE值最小為-0.301，表明預(yù)測值與實際值誤差小，故采用BP4-1組合因子作為模型輸入量，1994~016年天然徑流預(yù)測值與實際值過程對比見圖3。模型訓(xùn)練過程收斂情況如圖4，天然值和預(yù)測值回歸分析見圖5。

表2 武山站天然年徑流一致性修正統(tǒng)計表Tab2 Statistics of consistency correction of natural annual runoff at Wushan station

圖3 1994～2016年天然徑流預(yù)測值與實際值對比圖Figure 3 Comparison between predicted value and actual value of natural runoff from 1994 to 2016

圖4 均方差目標(biāo)值收斂曲線圖Figure 4 Convergence curve of mean square deviation target value

圖5 回歸分析結(jié)果Figure 5 Regression analysis results

3.3 綜合分析

模型輸入變量組合為BP4-1的MAPE值為-0.301，修正后1956～1993年天然徑流量4.875億m3。經(jīng)計算分析，采用降水徑流法一致性修正后的值為4.434億m3，本文提出的方法較傳統(tǒng)方法天然徑流修正值增加9.95%。兩種方法的修正后天然年徑流過程結(jié)果對比見圖6。

圖6 武山水文站年徑流系列曲線圖Figure 6 Annual runoff series curve of Wushan hydrological station

由圖6可得，降水徑流相關(guān)法修正結(jié)果相當(dāng)于修正前天然徑流向下平移，兩者峰型變化完全一致，而LM-BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型修正后徑流過程與修正前峰型變化趨勢相似，局部差異變化呈現(xiàn)非線性擬合關(guān)系。

修正前1956~2016年天然年徑流5.761億m3，采用LM-BP神經(jīng)網(wǎng)絡(luò)模型修正后天然徑流量4.422億m3，降水徑流相關(guān)法修正后天然徑流量4.147億m3，較原天然年徑流量分別減少23.2%、28.0%。本文提出的方法修正結(jié)果有利于水資源利用、經(jīng)濟(jì)社會規(guī)劃和可持續(xù)發(fā)展，且BP神經(jīng)網(wǎng)絡(luò)具有任意非線性關(guān)系擬合的特征，可以更好的逼近流域內(nèi)影響因子-天然徑流關(guān)系，綜合分析確定武山水文站1956~2016天然年徑流量4.422億m3較為合理。

4 討論與結(jié)論

1）經(jīng)水文時間序列趨勢與突變分析系統(tǒng)分析可得，武山水文站1956~2016年天然徑流量呈波動變化，減少趨勢顯著，傾向率-0.758億m3/10 a；綜合分析確定天然徑流突變點(diǎn)為1993年。

2）天然徑流預(yù)測模型引入誤差函數(shù)MAPE，根據(jù)MAPE值最小原則分析可得，武山水文站天然徑流LM-BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型的輸入量確定為氣溫（T）、降水量（P）、太陽黑子相對數(shù)（R）及武山站6~9月降水量（P6-9），最優(yōu)模型的MAPE值為-0.301。

3）基于LM-BP神經(jīng)網(wǎng)絡(luò)最優(yōu)模型的天然徑流一致性修正方法，利用了模型任意非線性逼近特征，綜合了與天然徑流相關(guān)的多種影響因子，1994~2016年天然徑流預(yù)測值與實際值擬合度較高，1956~1993年修正后天然年徑流過程與修正前峰型變化趨勢相似，修正值更為合理，可作為武山站天然年徑流量一致性修正優(yōu)選方法，且對其他水文站年徑流一致性修正有參考意義。

4）武山水文站位于渭河上游，近年來水土保持措施不斷強(qiáng)化，致使該區(qū)域下墊面條件發(fā)生明顯變化，徑流系列的一致性也發(fā)生改變，該區(qū)域徑流一致性修正分析具有典型性和代表性。BP神經(jīng)網(wǎng)絡(luò)模型可以更好地擬合現(xiàn)狀條件下的產(chǎn)匯流關(guān)系，但模型結(jié)構(gòu)和參數(shù)的確定方法不確定，有待進(jìn)一步研究分析。