劉明錦，張智涌?，王 賓，陳萬林

(1.成都水生態(tài)文明建設(shè)研究重點基地，四川 都江堰，611830；2.四川水利職業(yè)技術(shù)學(xué)院，四川 崇州，611231)

0 引言

山洪災(zāi)害是世界上危害最大的自然災(zāi)害之一，具有突發(fā)性強、破壞力大、難以預(yù)測等特點，容易造成大量人員傷亡、嚴重的財產(chǎn)損壞和環(huán)境災(zāi)難。由于自然地理環(huán)境、極端災(zāi)害天氣以及社會經(jīng)濟活動等各類因素的共同影響，造成了我國的山洪災(zāi)害問題呈現(xiàn)出頻發(fā)、多發(fā)的態(tài)勢。1950－2016年間，我國山洪災(zāi)害造成的直接經(jīng)濟損失高達1.72×105億元，2000年以來，每年因為山洪死亡的人數(shù)超過1000人[1]。山洪災(zāi)害的防御工作是我國防汛減災(zāi)工作的重點與難點。及時、準確地進行山洪災(zāi)害的預(yù)警預(yù)報工作，有利于指導(dǎo)受災(zāi)群眾快速撤離、減輕災(zāi)害損失、保障人民的生命財產(chǎn)安全，是目前最為有效可行的防災(zāi)減災(zāi)非工程措施，是山洪災(zāi)害防御的重要研究方向。

為了對山洪進行精準預(yù)測，不同領(lǐng)域的國內(nèi)外學(xué)者近年來從多個不同的角度進行了很多有意義的研究[2]。由于山洪和強對流降水的密切關(guān)系，氣象類的學(xué)者更傾向于對觸發(fā)山洪的強降水檢測進行研究[3]；水文地質(zhì)防汛等領(lǐng)域類的學(xué)者則傾向于采用水文類臨界雨閾值以及一些較為穩(wěn)定的水文匯流模型來進行計算，其本質(zhì)上也依賴于氣象信息[4]；隨著計算機建模的興起，也有學(xué)者利用需要進行山洪預(yù)測的計算機進行建模，考慮地形、地質(zhì)條件等因素，結(jié)合氣象信息對山洪進行預(yù)測分析。

總的來說，目前主流的山洪預(yù)測的核心思想還是結(jié)合各地區(qū)影響暴雨洪水的各類相關(guān)資料(歷史水文數(shù)據(jù)、氣象類、地形地貌數(shù)據(jù)等)，研究其規(guī)律并建立的相關(guān)經(jīng)驗關(guān)系。在以往的研究中，大多以山洪成因的各類影響因素為初始條件，用歷史數(shù)據(jù)作為主要參考對象，從而根據(jù)不同的水文地質(zhì)環(huán)境建立山洪預(yù)測模型。該類模型針對性較強，有歷史數(shù)據(jù)支撐，但本身存在模型靈活度不高、歷史數(shù)據(jù)準確性等問題，因此，山洪數(shù)據(jù)預(yù)測值誤差較大。同時因為各地區(qū)環(huán)境不一樣，山洪成因不同，因此存在模型遷移難度大的問題。

區(qū)別于以往側(cè)重于氣象資料或者計算機建模模擬山區(qū)構(gòu)造等進行的山洪預(yù)測，本研究的意義在于引入了機器學(xué)習模型，利用傳感器實時采集的降雨量作為輸入(采集時間可以設(shè)置較小的時間片)和最終形成的洪水量進行實時比對，形成大數(shù)據(jù)，通過機器學(xué)習模型不斷校正不同情況下降雨量與山洪之間的關(guān)系，最終能夠生成一個趨于理想條件的洪水預(yù)測模型。本研究期望利用機器學(xué)習和大數(shù)據(jù)技術(shù)為山洪預(yù)測提供技術(shù)支撐。

1 研究區(qū)域情況及數(shù)據(jù)來源

研究區(qū)域位于四川省阿壩州小金川流域。小金川河是大渡河的一級支流，流域面積5275km2，有北、東兩源。北源為撫邊河，在猛固橋處與東源沃日河匯合后，始稱小金川河，該處流域面積3681.95km2。猛固橋往下，河谷逐漸拓寬，水流趨緩，水量大增，兩岸集鎮(zhèn)、耕地、人口增多。其間分布有小金縣城(美興鎮(zhèn))、新格鄉(xiāng)以及丹巴縣的太平、半扇門、岳扎和中路四個鄉(xiāng)，總長60.8km，落差472.6m，平均比降7.8‰，小金川河分屬小金縣和丹巴縣管轄。分界線位于三叉溝口，上段屬小金縣，該處集雨面積4751.8km2，猛固橋到三叉溝，河段長26.51km，總落差143m，平均比降5.39‰，河道較為平緩。河谷多呈“U”形，河道彎度小，耕地和居民多分布于兩岸階地、灘地。本段匯入的支溝眾多，其中集雨面積大于100km2的有左岸的美沃溝、沙龍溝以及右岸的崇德溝和水卡子溝。其流域內(nèi)共有小金水文站、達維雨量站和兩河雨量站，其中小金水文站有完整的雨量數(shù)據(jù)和流量數(shù)據(jù)。研究區(qū)域基本情況見圖1。

圖1 研究區(qū)域基本情況

根據(jù)阿壩州水文局提供的采集樣本中的歷年雨量數(shù)據(jù)和山洪數(shù)據(jù)，取2013年5月1日后開始的數(shù)據(jù)。其中，小金水文站的雨量數(shù)據(jù)5343條，達維雨量站雨量數(shù)據(jù)5961條，兩河雨量站數(shù)據(jù)2339條，小金水文站流量數(shù)據(jù)有26759條，如表1所示。其中由于前期時間采集不一致和流域內(nèi)降雨時間和降雨量不相同，因此，雨量數(shù)據(jù)和流量數(shù)據(jù)產(chǎn)生條目不一樣，雨量站之間多是因為雨量不同而導(dǎo)致雨量條目不一致。

表1 研究區(qū)域基礎(chǔ)數(shù)據(jù)

經(jīng)過數(shù)據(jù)分析，發(fā)現(xiàn)部分數(shù)據(jù)中出現(xiàn)異常流量，排除掉因為梯級電站放水等異常流量，取當日中高段雨量為基準，篩選出研究數(shù)據(jù)，如表2所示。

表2 研究區(qū)域雨量

2 系統(tǒng)設(shè)計思路

本研究的設(shè)計思路是建立一套具有適用性廣、不限制區(qū)域的山洪預(yù)警模型。建立模型首先需要考慮山洪形成的原因。

根據(jù)國內(nèi)相關(guān)學(xué)者的研究，山洪的成因與降雨量有著密切的關(guān)系，隨著降雨量的不斷增多，不同下滲率的土壤濕度逐漸飽和，達到山洪爆發(fā)的臨界雨量，根據(jù)地質(zhì)地形的不同，山洪的形成時間會有所不同[5－6]。根據(jù)流域面積的不同，降雨與最終匯流成山洪有一定的時間差。

因此，可以看出，降雨是山洪爆發(fā)的最主要成因[7]，考慮到山區(qū)存在著資料匱乏等問題，在系統(tǒng)設(shè)計上盡量簡化以往需要考慮的其他各類因素影響[8]，設(shè)傳感器數(shù)據(jù)關(guān)系參數(shù)W1、W2…Wn來替代各類以往山洪計算中需要考慮的地形數(shù)據(jù)、地質(zhì)數(shù)據(jù)、土壤濕度指數(shù)、下滲率、山洪與降雨的時間差以及其他孕災(zāi)環(huán)境指標等山洪影響系數(shù)，降雨量為x和山洪水量y從傳感器中實時獲取并存入大數(shù)據(jù)平臺。系統(tǒng)設(shè)計思路如圖2所示。

圖2 系統(tǒng)設(shè)計思路

模型中通過NTP服務(wù)器進行數(shù)據(jù)時間同步，保證所有平臺的數(shù)據(jù)時間步調(diào)保持一致。傳感器采集的實時數(shù)據(jù)x通過采集器傳輸?shù)酱髷?shù)據(jù)中心進行存儲，本研究利用山洪預(yù)警模型機器學(xué)習算法，來計算各類傳感器中獲取的降雨量的實時數(shù)據(jù)x與最終匯流的山洪流量y之間的關(guān)系參數(shù)W。在不斷的采集過程中，大數(shù)據(jù)平臺會積累大量傳感器數(shù)據(jù)，山洪預(yù)警模型機器學(xué)習算法利用這些大數(shù)據(jù)來不斷優(yōu)化關(guān)系參數(shù)W，使其精度越來越高，從而形成一個不斷診斷與改進的數(shù)據(jù)優(yōu)化螺旋上升模型。

3 山洪預(yù)測模型

3.1 山洪預(yù)測模型參數(shù)

根據(jù)山洪形成的主要因素，給定一個大小為n的數(shù)據(jù)集Yi，xi1，xi2…xid(i＝1，…，n)，其中xi1…，xid是在山洪預(yù)測中第i個傳感器屬性的取值，yi為待預(yù)測的山洪值，即模型的預(yù)測目標，預(yù)測目標yi與傳感器之間的關(guān)系為線性假設(shè)組合，則之間的關(guān)聯(lián)公式為:

該式也是山洪預(yù)警模型中機器算法公式，其中w1，w2…，wd，b均為影響因素，w表示權(quán)重，b表示偏重。

3.2 山洪預(yù)測優(yōu)化

用yi表示為預(yù)測的山洪值，即模型的預(yù)測目標，y表示真實值。預(yù)測目標值與實際值是不一樣的，在這個模型中未知的影響因素為w1，w2…，wd，b，即要計算和優(yōu)化的目標。通過計算和優(yōu)化影響因素w1，w2…，wd，b，使其無限接近實際值y。

為了計算和優(yōu)化影響因素w1，w2…，wd，b，引入損失函數(shù)(Loss Function，或Cost Function)概念。通過傳感器采集的實際值與模型的預(yù)測值進行比對，損失函數(shù)輸出一個非負數(shù)的實值，這個實值表示模型誤差的大小，公式如下:其中，y＾

i表示預(yù)測值，yi表示實際值，即對于一個大小為n的數(shù)據(jù)集，MSE是n個數(shù)據(jù)預(yù)測結(jié)果誤差平方的均值。損失函數(shù)計算出來之后，需要通過不斷地優(yōu)化。優(yōu)化方式采用梯度下降法:設(shè)如果fx在點xn有定義且可微，則認為fx在點xn沿著梯度的負方向－f(xn)下降是最快的。反復(fù)調(diào)節(jié)xn，使得fx接近最小值或者極小值，調(diào)節(jié)的方式為:

上式中，λ代表學(xué)習率。

4 機器學(xué)習模型訓(xùn)練

定義了山洪預(yù)測模型結(jié)構(gòu)之后，通過以下幾個步驟(見圖3)進行機器學(xué)習的山洪預(yù)測模型訓(xùn)練:

圖3 機器學(xué)習流程

(1)采集傳感器數(shù)據(jù)初始值，初始化權(quán)重wi和偏置b參數(shù)值；

(2)利用網(wǎng)絡(luò)正向傳播來計算網(wǎng)絡(luò)輸出和損失函數(shù)；

(3)根據(jù)損失函數(shù)進行反向誤差傳播(backpropagation)，將網(wǎng)絡(luò)誤差從輸出依次向前傳遞，并更新網(wǎng)絡(luò)中的參數(shù)；

(4)重復(fù)2－3步驟，直至網(wǎng)絡(luò)訓(xùn)練誤差達到預(yù)期值，即預(yù)測值趨近目標值為止。

訓(xùn)練程序的代碼配置如下:

＃準備樣本特征X，訓(xùn)練初期用dawei，lianghe，xiaojin作為樣本特征

x＝data[[′dawei′，′lianghe′，′xiaojin′]]

＃接著用xiaojinshuiliang作為樣本輸出y

y＝data[[′xiaojinshuiliang′]]

＃劃分訓(xùn)練集和測試集

from sklearn.model＿selection import train＿test＿split

x＿train，x＿test，y＿train，y＿test＝train＿test＿split(x，y，test＿size＝0.2，random＿state＝0)

＃查看訓(xùn)練集和測試集的維度

＃print(′x＿train.shape:′，x＿train.shape)

＃print(′x＿＿test.shape:′，x＿test.shape)

＃print(′y＿train.shape:′，y＿train.shape)

＃print(′y＿test.shape:′，y＿test.shape)

＃從sklearn庫中導(dǎo)入函數(shù)

from sklearn.linear＿model import LinearRegression

＃執(zhí)行函數(shù)獲得一個模型

LR＝LinearRegression()＃未經(jīng)訓(xùn)練的機器學(xué)習模型

＃對模型傳入輸入數(shù)據(jù)x＿train和輸出數(shù)據(jù)y＿train

LR.fit(x＿train，y＿train)＃經(jīng)過訓(xùn)練的機器學(xué)習模型

＃輸出截距和各個系數(shù)

print(′LR.intercept＿:′，LR.intercept＿)

print(′LR.coef＿:′，LR.coef＿)

＃評價模型。這里使用MSE和RMSE來評價模型的好壞

y＿pred＝LR.predict(x＿test)

＃引入sklearn模型評價工具庫

from sklearn import metrics

print(″RMSE:″，np.sqrt(metrics.mean＿squared＿error(y＿test，y＿pred)))

import matplotlib.pyplot as plt

fig，ax＝plt.subplots()

＃畫散點圖

ax.scatter(y＿test，y＿pred)

ax.plot([y.min()，y.max()]，[y.min()，y.max()]，′k－－′，lw＝4)

＃設(shè)置標題

ax.set＿title(′Plot′)

＃設(shè)置X軸標簽

ax.set＿xlabel(′Measured′)

＃設(shè)置Y軸標簽

ax.set＿ylabel(′Predicted′)

＃顯示所畫的圖

plt.show()

為了更好地進行機器學(xué)習訓(xùn)練，將數(shù)據(jù)集分割為兩份:一份用于調(diào)整模型的參數(shù)，即進行模型的訓(xùn)練，模型在這份數(shù)據(jù)集上的誤差被稱為訓(xùn)練誤差；另外一份被用來測試，模型在這份數(shù)據(jù)集上的誤差被稱為測試誤差。訓(xùn)練模型的目的是通過從訓(xùn)練數(shù)據(jù)中找到規(guī)律來預(yù)測未知的新數(shù)據(jù)，所以測試誤差是更能反映模型表現(xiàn)的指標。

分割數(shù)據(jù)的比例主要考慮以下兩個因素:更多的機器學(xué)習訓(xùn)練數(shù)據(jù)會降低參數(shù)估計的方差，從而得到更可信的模型，即傳感器采集的數(shù)據(jù)越多，計算越精準；而更多的測試數(shù)據(jù)會降低測試誤差的方差，從而得到更可信的測試誤差。

5 結(jié)果驗證

結(jié)果驗證是利用python的sklearn深度學(xué)習框架平臺來進行數(shù)據(jù)的驗證和仿真。在仿真過程中，輸入了一系列傳感器數(shù)據(jù)集和相應(yīng)的山洪數(shù)據(jù)集來進行模型的訓(xùn)練和預(yù)測。下面的散點圖展示了使用模型后對山洪數(shù)據(jù)的預(yù)測。其中，在XY軸中橫坐標X表示當前情況下的山洪實際情況，縱坐標Y表示根據(jù)山洪預(yù)測模型預(yù)測的結(jié)果，當二者值完全相等的時候就會落在虛線上。所以模型預(yù)測得越準確，則點離虛線越近。

模型在運行時間T1之后的結(jié)果變化如圖4所示，程序運行結(jié)果如表3所示。

圖4 模型在運行時間T1之后結(jié)果變化

表3 模型運行時間T1數(shù)據(jù)

LR.intercept＿:[51.43598567]

LR.coef＿:[[2.66256014 7.81232586 1.88804097]]

RMSE:8.15

將小金水文站5343條數(shù)據(jù)，達維雨量站5961條數(shù)據(jù)，兩河雨量站2339條數(shù)據(jù)以及小金水文站流量數(shù)據(jù)26759條數(shù)據(jù)放入模型中進行機器學(xué)習，最終模型運行時間T2的結(jié)果變化如圖5所示，程序運行結(jié)果如表4所示。

圖5 模型在運行時間T2之后結(jié)果變化

表4 模型運行時間T2數(shù)據(jù)

(1096，4)

LR.intercept＿:[125.79028717]

LR.coef＿:[[6.13921908 2.48457902 2.91936099]]

RMSE:61.35

在運行時間T2之后，隨著采集數(shù)據(jù)量的不斷增多，發(fā)現(xiàn)整體數(shù)據(jù)偏差過大，離散程度較高，RMSE數(shù)據(jù)有其他影響因素。進行原因分析后針對現(xiàn)有模型中的參數(shù)進行了優(yōu)化，加入了連續(xù)下雨量、下雨時間、連續(xù)下雨天等參數(shù)，進行了模型的進一步細化，離散變化如圖6所示，程序優(yōu)化輸出后結(jié)果如下:

圖6 模型參數(shù)優(yōu)化后(在運行時間T2之后結(jié)果變化)

LR.intercept＿:[53.75315188]

LR.coef＿:[[0.33490087－0.62412436 2.64343157 1.22031563 1.06905883－1.25737151－0.18086965 1.16870534 0.94289013 1.27620644 0.46292934－2.19127078 1.84307365 1.98384269 0.70598535]]

RMSE:34.96

在運行T3時間之后，模型得到了優(yōu)化，離散程度更小，如圖7所示，程序運行結(jié)果如下:

圖7 模型在運行時間T3之后結(jié)果變化

LR.intercept＿:[53.84795419]

LR.coef＿:[[0.30243946－0.55534329 2.67885757 1.28878908 1.07345404－1.26625626 －0.1710308 1.15262515 0.8902583 1.26158854 0.45054248－2.13878691 1.84879842 1.94619908 0.70365609]]

RMSE:33.25

由此可以看出，隨著時間的不斷向前推移，該模型對于山洪的預(yù)測離散度會逐漸減小，真實值與預(yù)測值之間的誤差也會越來越小，導(dǎo)致產(chǎn)生誤差的原因分析，主要為:

一是山洪的成因機制較為復(fù)雜，具有非常明顯的非線性特征，山洪形成的影響因素眾多，降雨是其中最為重要的變量，其他如地形地質(zhì)、土壤、植被等因素對于山洪的形成也有一定程度的影響。對于面域較大的水文環(huán)境，可能受到上游梯級電站不定期放水影響從而導(dǎo)致數(shù)據(jù)偏移較大。

二是在機器學(xué)習過程中，需要大量數(shù)據(jù)，由于傳感器數(shù)量偏少，不能覆蓋整個山洪區(qū)域，因此，在機器學(xué)習中會導(dǎo)致參數(shù)偏移較大。隨著時間的推進，采集數(shù)據(jù)量不斷增多，機器學(xué)習會不斷優(yōu)化和調(diào)整參數(shù)，盡最大可能性提高數(shù)據(jù)精度。

三是小金流域上游小水電站影響，缺乏水電站放水量數(shù)據(jù)導(dǎo)致山洪流量無法修正。

以下是四川省阿壩州龍溪溝利用該模型進行的數(shù)據(jù)盲測，龍溪鄉(xiāng)位于汶川縣的最北部，距國道317線3.5km，西與理縣桃坪鄉(xiāng)相連，東北同茂縣接壤，南與克枯鄉(xiāng)相鄰，全鄉(xiāng)幅員面積214.3km2。龍溪溝流域面積187km2，設(shè)有龍溪和木扎兩個雨量站，上游無梯級電站。利用本洪水預(yù)測模型進行數(shù)據(jù)盲測，最終的程序運行結(jié)果如表5所示。

表5 模型盲測數(shù)據(jù)

LR.intercept＿:[1614.4154718]

LR.coef＿:[[－0.02352098 0.03386048 0.01865036－0.00408701 0.04258365 0.04459187－0.00338166 0.02211623]]

MSE:0.02

RMSE:0.12

由此可以看出，在面域較小，且上游無梯級電站的山區(qū)水文環(huán)境下，該模型具有較小的離散型，數(shù)據(jù)的分散程度小。通過模擬盲測，說明該模型具有一定的遷移性。

6 結(jié)論與討論

本文建立的基于機器學(xué)習和大數(shù)據(jù)的山洪預(yù)測模型，通過Scikit－Learn(sklearn)機器學(xué)習平臺對其進行了模擬，模擬有一定的效果，說明該模型適用于小規(guī)模山洪的預(yù)測，主要結(jié)論包括:

(1)山洪的成因機制較為復(fù)雜，具有非常明顯的非線性特征，山洪形成的影響因素眾多，降雨是其中最為重要的變量，其他如地形地質(zhì)、土壤、植被等因素對于山洪的形成也有一定程度的影響。對于面域較大的水文環(huán)境，可能受到上游梯級電站放水影響從而導(dǎo)致數(shù)據(jù)偏移較大。

(2)對于面域較小的水文環(huán)境，由于受到降雨的影響會更大，因此，數(shù)據(jù)準確度會隨著數(shù)據(jù)量的增多而提高，更適合于小規(guī)模部署。

(3)預(yù)測值與真實值之間存在一定的誤差，但隨著時間的逐漸向前推移，RMSE的離散程度在縮小。

本文中所涉及的洪水預(yù)測模型經(jīng)過前期的虛擬仿真和實地盲測，根據(jù)模型特性，需要去不斷尋找山洪與山洪成因的各種因素之間的聯(lián)系和影響關(guān)系，還需要不斷地學(xué)習和大量數(shù)據(jù)的對比計算，因此，采集器的數(shù)量和學(xué)習時間是一個重要不可或缺的條件，未來還將通過增加傳感器類型、增加傳感器數(shù)量等方式進一步提高計算精度。在算法改進方面，還將引入卷積神經(jīng)網(wǎng)絡(luò)，卷積神經(jīng)網(wǎng)絡(luò)是將來模型的一個重要方向。