何中政，方 麗，劉 萬，遲英凡，王永強，付吉斯，熊 斌

（1. 南昌大學(xué)工程建設(shè)學(xué)院，江西 南昌 330031； 2. 長江水利委員會長江科學(xué)院，湖北 武漢 4 30010；3. 南昌大學(xué) 鄱陽湖環(huán)境與資源利用教育部重點實驗室，江西 南昌 330031）

0 引 言

流域徑流過程受眾多因素的相互作用與影響，使得徑流時間序列表現(xiàn)出復(fù)雜的非線性特性［1］，具有顯著的非平穩(wěn)、高維度和模糊性等復(fù)雜特征［2，3］，研究預(yù)測精度高和穩(wěn)定性好的徑流預(yù)測模型仍是當(dāng)前水文水資源研究領(lǐng)域熱點問題之一。而目前水文預(yù)測方法主要可分為過程（機理）驅(qū)動和數(shù)據(jù)驅(qū)動模型兩大類：過程驅(qū)動模型以水文過程的物理機制或規(guī)律為基礎(chǔ)［4］，需要大量且可靠的水文數(shù)據(jù)來建立模型，從而進行流域徑流模擬或預(yù)測，但由于水文模型對水循環(huán)過程的簡化處理，往往難以精準(zhǔn)刻畫非線性的徑流過程，且由于部分水文數(shù)據(jù)資料缺失或可靠性差等因素影響，過程驅(qū)動模型在實際徑流預(yù)報中存在局限性［5］；而隨著計算機硬件、軟件和理論技術(shù)的發(fā)展，數(shù)據(jù)驅(qū)動模型依托統(tǒng)計學(xué)或其他人工智能技術(shù)，無需考慮水文過程物理機制［6］，通過對時間序列進行分析從而建立對預(yù)測對象和預(yù)測因子間復(fù)雜非線性關(guān)系刻畫更為完備的模型，使得數(shù)據(jù)驅(qū)動模型在水文數(shù)據(jù)資料缺少或復(fù)雜非線性的徑流過程將具有更好的適用性［2］。近年來，學(xué)者們將多種常見的人工智能模型引入到水文預(yù)測中，如前饋式神經(jīng)網(wǎng)絡(luò)［7］（Back-Propagation neural network，BP）、多元線性回歸［8］（Multiple Linear Regression，MLR）、支持向量機［9］（Support Vector Machine，SVM）等，相關(guān)研究均取得了較好的效果［10］，如何應(yīng)用人工智能技術(shù)開展流域徑流預(yù)測已經(jīng)成為未來發(fā)展的重要研究方向之一［11］。

高斯過程回歸（Gaussian Process Regression，GPR）作為一種新型的機器學(xué)習(xí)方法，采用高斯過程對先驗數(shù)據(jù)進行回歸分析，以統(tǒng)計學(xué)理論為基礎(chǔ)，在處理高維、非線性、小樣本等復(fù)雜的回歸問題時都具有很強的泛化能力，因而在各類研究領(lǐng)域得到應(yīng)用［12］。黃亞等人［3］采用GPR 模型在廣西天峨水文站進行短期徑流預(yù)測；孫斌等人［13］采用GPR 模型在進行短期風(fēng)速預(yù)測；李成宇［14］采用GPR 模型進行了土石壩沉降預(yù)測研究；周靖楠等人［15］采用GPR 模型與MI-KPCA 結(jié)合在北汝河進行中長期徑流預(yù)測。劉龍龍等人［16］采用GPR 模型與蜻蜓算法結(jié)合對居民社區(qū)時用水量動態(tài)實時區(qū)間預(yù)測。GPR 模型已在非線性回歸分析問題中得到成功應(yīng)用，GPR 回歸分析能力不僅取決于模型參數(shù)，還受核函數(shù)影響，但目前核函數(shù)的影響還未得到進一步研究。Scholkopf 和Smola［17］指出支持向量機和高斯過程兩種算法的核函數(shù)是等價的，并提出幾種高斯過程中可選擇的核函數(shù)類型。因此，本文分析了不同核函數(shù)作用下GPR 預(yù)測模型效果，并提出了一種基于指數(shù)核函數(shù)GPR 的流域短期徑流預(yù)測模型，以江西省贛江流域吉安水文站徑流過程為對象例開展實例分析。

1 研究方法

1.1 高斯過程回歸預(yù)測模型基本原理

高斯過程（Gaussian Process，GP）是概率論和數(shù)理統(tǒng)計中隨機過程的一種，是一系列服從正態(tài)分布的隨機變量的組合。GP繼承了正態(tài)分布的諸多性質(zhì)，其特征主要由數(shù)學(xué)期望和協(xié)方差函數(shù)決定。GPR 則是對數(shù)據(jù)進行回歸分析，來建立服從GP 先驗信息的非參數(shù)模型［18］。假設(shè)一個含有n對相互獨立觀測數(shù)據(jù)的學(xué)習(xí)樣本{，y}，其中是由n個輸入向量組成的輸入集，而y={y1，y2，…，yn}是由n個相對應(yīng)的一維輸出組成的輸出集。那么GP 的均值μ(x)和方差k(x，x′)可由下式確定：

式中：x，x′代表樣本中任意隨機變量，E代表數(shù)學(xué)期望。

因此GP可定義為：

而在實際應(yīng)用中，通常會進行數(shù)據(jù)預(yù)處理使得均值函數(shù)為0，故GP即也表示為：

對于實際工程中的回歸問題，還應(yīng)考慮輸出y受噪聲影響：

式中：y是受到噪聲污染的觀測值，f為函數(shù)值，x為輸入向量，噪聲ε服從均值為0、方差為σ2n的高斯分布。

從而可以推求得到y(tǒng)的先驗分布，見公式（6）；以及y和f(x*)的聯(lián)合先驗分布，見公式（7）。

式中：x*為測試輸入向量，I是n維單位矩陣，是n×n階的GP 協(xié)方差矩陣，是輸入集與測試輸入向量x*之間的n× 1階協(xié)方差矩陣，K(x*，)同理可知，k(x*，x*)為x*自身的協(xié)方差。

由此可知預(yù)測值f(x*)的后驗分布：

其中：

1.2 高斯過程回歸的核函數(shù)與超參數(shù)優(yōu)化

GPR 使用GR 作為先驗，其訓(xùn)練結(jié)果與預(yù)先設(shè)定的核函數(shù)有關(guān)。在隨機樣本均值為0 的前提下，核函數(shù)就與協(xié)方差函數(shù)的實際意義相近，可進行樣本間相關(guān)性的描述。不同核函數(shù)作用下的高斯過程回歸特征各不相同，高斯過程模型常用的核函數(shù)主要有以下幾種：

（1）徑向基核函數(shù)：

式中：r=‖X1-X2‖，σ為超參數(shù)，表征了學(xué)習(xí)樣本間的相似度。徑向基核函數(shù)也被稱為高斯核，常被應(yīng)用于SVM和GPR等各類機器學(xué)習(xí)算法中，參數(shù)簡單，但在處理樣本數(shù)量大或特征多時效果一般。

（2）二次有理核函數(shù)：

式中：α，l為超參數(shù)。二次有理核函數(shù)可以理解為是無窮個徑向基核函數(shù)的線性疊加，當(dāng)α→∞時，二次有理核函數(shù)等價于σ=l的徑向基核函數(shù)。與徑向基核函數(shù)相比，更為省時，作用域廣，但是對參數(shù)十分敏感。

（3）馬頓核函數(shù)：

式中：l，ν為超參數(shù)，Kν為修正貝塞爾函數(shù)。當(dāng)ν→∞時，相當(dāng)于以l為特征尺度的徑向基核函數(shù)。

（4）指數(shù)核函數(shù)：

指數(shù)核函數(shù)是馬頓核在ν= 0.5的特殊形式，這樣的變動減少對參數(shù)的依賴性降低，使得模型的訓(xùn)練學(xué)習(xí)難度大大降低，適用情景相對狹窄。

而核函數(shù)中的超參數(shù)一般可依據(jù)極大似然法，建立基于學(xué)習(xí)樣本條件概率的對數(shù)似然函數(shù)，見式（15）；然后可通過非線性規(guī)劃方法或群智能搜索方法尋找出超參數(shù)的近似最優(yōu)解。

1.3 預(yù)報因子篩選

由于徑流預(yù)測往往是采用輸入多個預(yù)測因子對預(yù)測量進行輸出，研究采用多重相關(guān)系數(shù)來篩選預(yù)報因子。多重相關(guān)系數(shù)是測量單個因變量與其他多個自變量之間線性相關(guān)程度的指標(biāo)。多重相關(guān)系數(shù)無法直接計算，可采用預(yù)測量y和其他多個預(yù)測因子x1，x2，…，xk線性組合的簡單相關(guān)系數(shù)進行測算。

首先，根據(jù)預(yù)測量y和其他多個預(yù)測因子x1，x2，…，xk進行回歸，建立如式（16）的回歸方程，然后進一步計算簡單相關(guān)系數(shù)，即為預(yù)測量y和多個預(yù)測因子x1，x2，…，xk之間的多重相關(guān)系數(shù)R。

式中：β0，β1，β2，…，βk為回歸系數(shù)。

2 實例研究

2.1 研究區(qū)域與資料概況

贛江是江西省內(nèi)第一大河流，地理位置為113.58°～116.63°E，24.52°～28.75°N。其發(fā)源于贛閩邊界武夷山西麓，主河道長約823 km，干支流自南向北貫穿整個江西省，地貌多以山地、丘陵為主，豐水期為3-8月，枯水期為9月-次年2月。

圖1 贛江中下游控制斷面分布示意圖Fig.1 Distribution diagram of control section in the middle and lower reaches of Ganjiang River

研究選取了贛江中游棟背、吉安關(guān)鍵斷面開展實例分析。棟背站位于江西省萬安縣棟背村，控制流域面積約4.02 萬km2；吉安站是贛江中游主要控制站，控制流域面積約5.62 萬km2；棟背站到吉安站的時滯約為18～24 h，棟背站對吉安站的流量過程有較大的影響。研究收集整理了兩個水文站點從2010 年1月1 日0∶00 至2018 年4 月14 日18∶00 的實測徑流數(shù)據(jù)，時段間隔為6 h，共計20 000 個數(shù)據(jù)。

2.2 基于高斯過程回歸的徑流預(yù)測流程

依托項目研究收集得到的數(shù)據(jù)資料，開展基于GPR 模型的徑流預(yù)測研究。具體實施步驟為：① 整編收集得到的徑流數(shù)據(jù)資料，將前80%作為預(yù)測模型訓(xùn)練期，后20%作為預(yù)測模型檢驗期；②根據(jù)率定期徑流資料，計算預(yù)報因子-預(yù)測量的多重相關(guān)系數(shù)，開展相關(guān)性分析，從而選定預(yù)報因子；③選取評價指標(biāo)，便于直觀表現(xiàn)模型的預(yù)測效果；④采用GPR 對預(yù)報因子和預(yù)報變量進行建模，用訓(xùn)練期資料進行參數(shù)率定；⑤將訓(xùn)練好的GPR模型用于檢驗期的徑流預(yù)報，檢驗?zāi)Ｐ皖A(yù)報效果。

2.3 評價指標(biāo)選取

為評價預(yù)測結(jié)果的綜合表現(xiàn)，選取了均方根誤差（Root Mean Squared Error，RMSE）、平均絕對誤差（Mean Absolute Error，MAE）、確定性系數(shù)（Deterministic coefficient，DC）和合格率（Qualified Rate，QR）評價指標(biāo)，具體計算公式如下：

圖2 基于高斯過程回歸的徑流預(yù)測流程圖Fig.2 Flow chart of runoff prediction based on Gaussian process regression

式中：Si為模擬值；Oi為實測值；為實測值的均值；k為樣本中相對誤差小于20%的樣本個數(shù)；n為測試樣本個數(shù)。評價指標(biāo)中，RMSE和MAE取值范圍為(0，∞)，其值越小表明模型預(yù)測誤差越小；DC取值范圍為[0，1]，其值越大表明預(yù)測值和真實值的相關(guān)性越大；QR取值范圍為[0，100%]，其值越大表明模型預(yù)測精度越高。

2.4 預(yù)測因子篩選

建立短期徑流預(yù)測模型首先需要篩選預(yù)測因子，記T為當(dāng)前時段，OT和IT分別表示本站點和上游站點T時段徑流量。采用IBM SPSS Statistics 26.0 軟件開展見2.3 小節(jié)中的多重相關(guān)系數(shù)分析，多重相關(guān)性分析結(jié)果如表1和表2所示。

表1 僅考慮吉安站的多重相關(guān)性分析結(jié)果Tab.1 The multiple correlation analysis results of Ji'an Station

表2 考慮棟背站和吉安站的多重相關(guān)性分析結(jié)果Tab.2 The multiple correlation analysis results of Dongbei station and Ji'an station

表1 和表2 僅列出了與吉安徑流多重相關(guān)性系數(shù)最大，且預(yù)測因子時段最短的組合。對比表1 和表2 可知，增加棟背徑流信息對“預(yù)測因子—預(yù)測量”復(fù)相關(guān)性系數(shù)增益較為明顯，因此研究選擇表2中預(yù)測因子開展短期徑流預(yù)測。

2.5 徑流預(yù)測實例分析

采用數(shù)據(jù)留出法［19］進行評估，將樣本數(shù)據(jù)的前80%作為模型訓(xùn)練樣本，后20%作為模型檢驗樣本，即訓(xùn)練期為2010 年1月1 日0∶00 至2016 年11 月19 日18∶00 的實測徑流數(shù)據(jù)，時段間隔為6 h；檢驗期為2016 年11 月19 日18：00 至2018 年4 月14日的實測徑流數(shù)據(jù)，時段間隔為6 h。根據(jù)表2 中的預(yù)測因子，建立28 個基于GPR 的短期徑流預(yù)測模型，分別對面臨T+1～T+28 時段的吉安徑流進行預(yù)測。為分析不同核函數(shù)對GPR 短期徑流預(yù)測模型的作用效果，研究分別選用了有理二次、徑向基、馬頓和指數(shù)核函數(shù)建立不同的GPR 短期徑流預(yù)測模型，還加入了MLR、回歸樹（Regression Tree，RT）、SVM、BP 等模型方法的預(yù)測結(jié)果作為對比。以上模型的訓(xùn)練以及檢驗均采用MATLAB R2020b 進行處理，GPR、MLR、RT 和SVM 模型訓(xùn)練基于Regression Learner 模塊，BP 模型訓(xùn)練基于采用的Neural Net Fitting模塊。

圖3 給出了吉安站不同模型短期徑流預(yù)測在檢驗期的結(jié)果，見T+1～T+27 時段4 種評價指標(biāo)的熱力圖。4 張熱力圖均為顏色越深，檢驗期預(yù)測效果越差。相關(guān)實驗結(jié)果表明：①從4種指標(biāo)綜合來看，研究建立的4 種GPR 模型和RT 模型均優(yōu)于MLR、SVM、BP 模型，在DC指標(biāo)上差異最為顯著，4 種GPR 模型和RT 模型的DC均在0.9以上，而MLR、SVM、BP 模型DC前2個時段大于0.9，隨預(yù)見期增長DC銳減至不足0.4；②從4 種GPR模型對比來看，有理二次、指數(shù)GPR 優(yōu)于徑向基和馬頓GPR，結(jié)合核函數(shù)表達(dá)式（11）～（14）來看，徑向基和馬頓核函數(shù)明顯比有理二次、指數(shù)核函數(shù)更復(fù)雜，這表明復(fù)雜的核函數(shù)形式不一定能提升GPR 模型非線性回歸分析能力；③隨著預(yù)測期的增加，MLR、SVM、BP 模型的4種評價指標(biāo)均呈現(xiàn)劣化的趨勢，4種GPR 模型和RT 模評價指標(biāo)有略微波動，但整體較為穩(wěn)定；④綜合各項評價指標(biāo)來看，有理二次、指數(shù)GPR 模型結(jié)果最為突出，而指數(shù)GPR 模型的RMSE、MAE、QR明顯優(yōu)于有理二次GPR，指數(shù)GPR 模型在28個時段的QR中僅有1處合格率不足100%，明顯少于有理二次GPR 的22 處。綜上所述，本文提出的指數(shù)GPR在吉安站短期徑流預(yù)測中顯著優(yōu)于其他6種模型方法。

圖3 吉安站檢驗期各模型評價指標(biāo)熱力圖Fig.3 Heat map of each model evaluation index for the testing period at Ji'an Station

圖4和表3分別給出了吉安站不同模型短期徑流預(yù)測T+28時段的預(yù)測結(jié)果和4 種評價指標(biāo)統(tǒng)計結(jié)果，7 種模型方法在訓(xùn)練期和檢驗期4 項評價指標(biāo)變化不大。此外，研究發(fā)現(xiàn)4 種GPR 模型檢驗期4項評價指標(biāo)在T+1到T+9時段逐漸劣化，4項評價指標(biāo)隨后在T+10到T+28時段變好并趨于穩(wěn)定。結(jié)合多重相關(guān)性分析結(jié)果，可知T+10 時段以后預(yù)測模型輸入信息激增，有助于數(shù)據(jù)驅(qū)動的GPR 模型開展回歸預(yù)測分析。這也說明應(yīng)用多重相關(guān)性系數(shù)分析篩選出的“預(yù)測因子—預(yù)測量”關(guān)系，并不能完全反映自身和上游站點徑流的非線性關(guān)系。

表3 吉安站T+28時段訓(xùn)練期和檢驗期結(jié)果評價指標(biāo)Tab.3 Evaluation indexes for the results of the training and testing period of the T+28 period at Ji'an Station

圖4 吉安站檢驗期T+28時段預(yù)測結(jié)果Fig.4 Forecast results for the T+28 period of the testing period at Ji'an Station

3 結(jié) 論

（1）以贛江流域吉安水文站為對象的預(yù)測步長為6 h 預(yù)見期為7 d 的實例分析為例，應(yīng)用不同核函數(shù)的GPR 模型預(yù)測結(jié)果表現(xiàn)存在明顯差異，不同方法預(yù)測表現(xiàn)由好到差分別為指數(shù)GPR、有理二次GPR、RT、馬頓GPR、徑向基GPR、SVM、MLR、BP。指數(shù)GPR 預(yù)測效果最好，28 個時段的DC和QR均分別接近1和100%。

（2）從4 種GPR 模型對比來看，徑向基和馬頓核函數(shù)表達(dá)式比有理二次、指數(shù)核函數(shù)更復(fù)雜，但預(yù)測結(jié)果不如有理二次、指數(shù)核函數(shù)的GPR 模型，說明復(fù)雜的核函數(shù)形式不一定能提升GPR模型非線性回歸分析能力。

（3）4 種GPR 模型檢驗期預(yù)測精度在T+1 到T+9 時段逐漸降低，但在T+10 到T+28 時段，隨預(yù)測模型輸入因子信息激增，預(yù)測精度提高并趨于穩(wěn)定。

（4）多重相關(guān)性系數(shù)分析篩選出的“預(yù)測因子—預(yù)測量”關(guān)系，主要依據(jù)線性關(guān)系，并不能反映自身預(yù)測量和上游站點預(yù)測因子的非線性關(guān)系。在未來的研究中，需探究能夠精準(zhǔn)描述預(yù)測因子與預(yù)測量復(fù)雜非線性關(guān)系的分析方法，從而更為準(zhǔn)確的篩選預(yù)測因子。