若干新型群體智能算法優化高斯過程回歸的年降水量預測

李 杰，崔東文

（1.云南省水利水電勘測設計研究院，昆明 650000；2.云南省文山州水務局，云南 文山 663000）

0 引 言

降水量預測是水資源合理開發利用的重要基礎，同時是緩解區域水資源供需矛盾、提高防洪抗旱應急能力、保障區域水生態安全的重要支撐[1]。近年來，以神經網絡為代表的機器學習技術因非線性問題處理能力強、無需問題解析函數、泛化能力好等優點，已在降水量預測研究中得到廣泛應用，如BP 神經網絡[2,3]、支持向量機（SVM）[4,5]、長短時記憶神經網絡（LSTM）[6,7]等。然而，BP神經網絡雖然能夠以任意精度逼近函數，但存在模型結構復雜、易陷入局部最優、調節參數選取困難等不足；SVM 方法雖然需要的樣本較少，但是其預測精度受限于核函數及其參數的準確選取，且易陷入局部極值；LSTM 網絡雖然在時間序列預測方面具有較大優勢，但LSTM 存在內部參數多、收斂速度慢、系統消耗資源大等不足。高斯過程回歸（Gaussian process regression，GPR）是基于貝葉斯理論與統計學習理論的非參數監督學習方法，適于處理高維數、非線性等復雜的回歸問題，已在時間序列預測分析、動態系統模型辨識、系統控制等多個領域得到廣泛應用，并取得良好效果。然而，GPR 預測效果對超參數取值的依賴程度較高，目前主要采用實驗試湊、經驗選擇等方法對GPR 超參數進行選取，存在隨機性大、容易陷入局部最優等問題。針對這一問題，粒子群優化（PSO）算法[8,9]、北方蒼鷹優化（NGO）算法[10]、鯨魚優化算法（WOA）[11]、人工蟻群優化（ACO） 算法[12]等群體智能算法（swarm intelligence algorithms，SIA）償試用于GPR 超參數優化，并取得較好的優化效果。

由于降水受地理位置、地形條件、大氣環流、人類活動等諸多因素影響，其時間序列常表現出高噪聲、非線性、非平衡性等特性[13]，若直接預測會使結果產生較大誤差。為充分挖掘降水量時間序列所含信息，小波分解（WD）、變分模態分解（VMD）、完全集合經驗模態分解（CEEMDAN）、總體經驗模態分解（EEMD）、奇異譜分析（SSA）等分解技術已在降水量時間序列數據處理中得到應用，如王文川等[14]提出的基于小波分解（WD）和郊狼優化（COA）算法的長短期記憶神經網絡（LSTM）降水量預測模型，徐冬梅等[15]提出VMD-時間卷積網絡（TCN）月降水量組合預測模型，羅上學等[16]建立的CEEMDAN-LSTM 耦合模型，楊倩等[17]基于總體經驗模態分解（EEMD） 和長短時記憶神經網絡（LSTM） 構建的EEMDLSTM 耦合模型，張以晨[18]等建立的奇異譜分解（SSA）-支持向量回歸機（SVR）耦合模型。然而，在實際應用中，上述分解技術存在一定不足，如WD只對低頻部分進行再分解，缺乏對時序數據高頻部分的分析，降低了分解效果；VMD 在非線性、非平穩數據分解中具有較高的精確性，但分解數值K 的合理選取對預測結果影響較大；CEEMDAN、EEMD 方法雖然在一定程度上解決了EMD 模型混疊和虛假分量問題，但人為經驗添加噪聲、分解分量過多等不足制約了其應用；SSA技術雖然具有較好的分解效果，但存在分解分量多、預測復雜度高、計算規模大等缺點。小波包變換（wavelet packet transform，WPT）衍生于WD，與之不同的是，WPT 同時將低頻、高頻信號再次分解，能將原始信號分解為更具規律的子序列分量，在各行業領域具有廣泛的應用。

因此，為提高降水量預測精度，改進GPR 預測性能，同時拓展SIA 在GPR 超參數優化中的應用，本文基于WPT 分解技術、孔雀優化算法（peafowl optimization algorithm，POA）、沙貓優化（sand cat swarm optimization，SCSO）算法、獵豹優化（cheetah optimization，CO）算法和GPR 方法，建立WPTPOA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型，同時構建基于支持向量機（SVM）的WPT-POA-SVM、WPT-SCSOSVM、WPT-CO-SVM 模型、基于RBF 神經網絡的WPT-POARBF、WPT-SCSO-RBF、WPT-CO-RBF 模型和未經優化的WPT-GPR 模型作對比分析模型，并利用云南省文山州1956-2021 年年降水量數據對10 種模型進行率定與驗證，旨在檢驗WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型用于年降水量預測的可行性。

1 材料與方法

1.1 數據來源

本文選用文山州1956-2021年年尺度降水量數據（數據來源于文山州歷年水資源公報和文山州統計年鑒），過程線如圖1 所示。從圖1 可以看出，文山州年降水量序列波動性較大，復雜程度較高，呈現出較強的非線性和非平穩性，不利于直接預測。文山州位于云南省東南部，轄文山、硯山等7 縣1市，國土面積31 456 km2，分屬珠江、紅河兩大流域。珠江流域面積17 309 km2，占全州總面積的54.5%，主要有南盤江、清水江、馱娘江、西洋江等，多年平均水資源量75.36 億m3；紅河流域面積14 147 km2，主要有盤龍河、八布河、南利河、迷福河等，多年平均水資源量82.34 億m3。文山州水資源總量豐富，多年平均降水量1 208 mm，水資源總量157.7 億m3，屬相對豐水地區。區域降水高度集中，汛期（5-10 月）降水量約占全年降水量的80%，其中約60%的降水集中在7-9月，降水量年內高度集中導致水旱災害頻發和水資源供需矛盾突出。

圖1 文山州1956-2021年降水量Fig.1 Precipitation in Wenshan Prefecture from 1956 to 2021

1.2 研究方法

1.2.1 小波包變換（WPT）

WPT 能同時對信號低頻部分和高頻部分進行分解，更適用于降水量時間序列分解。WPT 對降水量原始信號進行分解的公式為[19-21]：

重構算法為：

式中各參數意見詳見文獻[19-21]。

1.2.2 孔雀優化算法（POA）

POA 是Jingbo Wang 等人于2022 年受綠孔雀群求偶、覓食和追逐行為啟發而提出一種新型元啟發式算法。該算法通過雄孔雀、雌孔雀和幼年孔雀來模擬種群在覓食過程中的群體行為，并基于孔雀種群角色劃分建立雄孔雀、雌孔雀和幼年孔雀位置更新算子來求解待優化問題[22]。POA數學描述如下：

（1）角色劃分。POA 將孔雀種群分為3 個角色：雄孔雀、雌孔雀和幼年孔雀。在優化過程中，所有個體都根據適應度值進行排名，并將適應度值最優的前五只孔雀視為雄孔雀Xpci（i=1,2,…,5），剩余的前30%的孔雀視為雌孔雀XPh，其余的孔雀視為幼年孔雀XPhC。

（2）雄孔雀位置更新。雄孔雀擁有華麗的羽毛，在其找到充足的食物后四處游走，并通過搖晃羽毛來求偶，雄孔雀越漂亮，吸引的雌孔雀就越多。POA 中，孔雀擁有的適應度值越高，其圍繞食物源繞圈的概率就越大，圈半徑越小；適應度值差的孔雀更容易原地旋轉，圓半徑較大。雄孔雀位置更新描述如下：

式中：Xpci為第i只雄孔雀位置，i=1,2,…,5；Rs為旋轉半徑；Xr為隨機向量，描述為Xr= 2 · rand(1,Dim) - 1；‖Xr‖為Xr的模數；Dim 為問題維度；r1、r2、r3、r4為均勻分布在[0,1]范圍內的隨機數；t為當前迭代次數。

（3）雌孔雀接近行為。雌孔雀看到雄孔雀求偶舞蹈時，往往會先靠近雄孔雀再觀察四周，雌孔雀被吸引的概率與雄孔雀的適應度值成正比。雌孔雀位置描述為：

式中：XPh為雌孔雀位置；r5為[0,1]范圍內均勻分布的隨機數；θ為雌孔雀接近系數，描述為θ=θ0+θ1-θ0(t/tmax)，t、tmax為當前迭代次數和最大迭代次數；θ0、θ1為接近系數初始值和最大值，分別取0.1和1。

（4）幼年孔雀搜索行為。除了接近具有良好食物資源（高適應度值）的雄性孔雀外，幼年孔雀還隨機搜索，希望在搜索空間中找到更高質量的食物資源。幼年孔雀位置描述為：

式中：XPcC為幼年孔雀位置；α、δ為隨迭代次數動態變化的因子；Levy 為萊維飛行，一種隨機游走策略；XSPc、XPcC為隨機選定的雄孔雀位置和幼年孔雀位置；其他參數意義同上。

（5）孔雀間互動行為。由于雄孔雀Xpc1擁有最好的食物資源，其余四只雄孔雀將被吸引逐漸向Xpc1移動。

1.2.3 沙貓優化（SCSO）算法

SCSO 是Amir Seyyedabbasi 等人于2022 年提出一種新型優化算法。該算法靈感來自于沙貓特殊的低頻噪聲檢測能力，該能力有助于沙貓在地面或地下快速找到并捕捉獵物。SCSO算法中，沙貓覓食行為分搜索獵物（勘探）和攻擊獵物（開發）兩個階段，并通過自適應機制在勘探和開發間保持平衡[23]。SCSO數學描述簡述如下：

（1）初始化。設置沙貓種群規模N，利用式（1）初始化沙貓個體位置Xi。

式中：Xi表示第i只沙貓個體初始位置；ub、lb 分別表示搜索空間上、下限值；rand(0,1)表示介于0 和1 之間均勻分布的隨機數。

（2）搜索獵物（探索）。在該階段，每只沙貓根據最佳候選位置、當前位置及其敏感度范圍r更新自己位置，使其能在搜索區域中獲得局部最優值。沙貓位置更新描述如下：

式中：Xi(t+ 1)為第i只沙貓第t+ 1 次迭代位置；Xbc(t)為第t次迭代沙貓最佳候選位置；Xc(t)為第t次迭代沙貓位置，即沙貓當前位置；r為沙貓敏感度范圍，描述為r=rG· rand(0,1)，其中rG為靈敏度控制參數，描述為rG=SM-(2SM·t)/(T+t)，SM為沙貓聽覺特性參數，取值為2；t、T分別為當前迭代次數和最大迭代次數；其他參數意義同上。

SCSO 通過定義探索與開發之間轉換的平衡參數R來保持探索與開發之間的平衡，數學描述如下：

式中：R為探索與開發之間轉換的平衡參數；其他參數意義同上。

（3）攻擊獵物（開發）。在該階段，SCSO 通過輪盤選擇算法為每只沙貓選擇一個隨機角度進行位置更新以接近獵物。沙貓位置更新描述如下：

式中：Xb(t)為迄今為止沙貓最佳位置；θ為隨機角度，θ介于0和360之間；其他參數意義同上。

1.2.4 獵豹優化（CO）算法

CO 算法是Mohammad AminAkbari 等人于2022 年受自然界獵豹狩獵啟發而提出一種新型群體智能優化算法。該算法通過模擬獵豹在狩獵過程中搜索、坐等和攻擊3種策略來實現位置更新，即待優化問題的求解[24]。CO算法數學描述如下：

（1）初始化。與其他群體智能優化算法類似，CO 算法也是從種群初始化開始。設在d維搜索空間中，獵豹初始化位置描述為：

式中：Xi,j為第i頭獵豹第j維位置；UBj、LBj為第j維搜索空間上、下限值；rand 為介于0 和1 之間的隨機數；n為獵豹種群規模；d為問題維度。

（2）搜索策略。獵豹在其領地（搜索空間）或周圍區域進行全范圍掃描或主動搜索，以找到獵物。該策略數學描述為：

（3）坐等策略。在搜索模式下，獵物可能會暴露在獵豹視野中，在這種情況下，獵豹的每一個動作都可能會導致獵物逃跑。為避免該情況發生，獵豹采取坐等伏擊策略（躺在地上或躲進灌木叢）以接近獵物。該策略數學描述為：

式（12）各參數意義同上。該策略不但提高狩獵成功率（獲得取優解），而且避免CO過早收斂。

（4）攻擊策略。在CO 算法中，每頭獵豹都可以根據逃跑獵物、領頭獵豹或附近獵豹的位置來調整自己的位置，以獲得最佳攻擊位置。該策略數學描述為：

1.2.5 高斯過程回歸（GPR）

GPR 是通過有限的高維數據來擬合出相應的高斯過程，從而來預測任意隨機變量下的函數值。設訓練集(X,y)={(Xi,yi)|i=1,2,…,n}，其中X=[x1,x2,…,xn]代表d維的輸入向量矩陣，y=[y1,y2,…,yn]表示輸出值。高斯過程（Gaussian process，GP）f(x)的有限維隨機變量都服從一個多元高斯分布，其性質由均值函數m(x)和協方差函數（核函數）k(x,x＇)決定，因此高斯過程可定義為[8-12]：

考慮噪聲后，GPR回歸模型可表示為：

式中：ε為高斯噪聲，且滿足ε～N(0,σ2)。

考慮噪聲觀測值的先驗分布為：

式中：K(X,X)為n×n階協方差矩陣；In為單位矩陣。

觀測值y和預測值f*的聯合先驗分布為：

測試過程中根據Bayes原理求得后驗分布為：

協方差函數k(x,x＇)（又稱“核函數”）體現了輸入樣本之間的相似關系，決定了GPR 的預測精度。GPR 核函數中超參數選取對模型精度有重要影響，本文基于二次有理式協方差函數（RQ）構建GPR，并采用上述POA、SCSO、CO 算法對GPR核函數參數、噪聲方差、協方差3個超參數進行調優，以改善GPR模型預測性能。

1.3 建模流程

WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型的預測步驟歸納如下（其他模型參考實現）：

步驟一：利用WPT 將實例年降水量時序數據進行2 層分解，得到1 個周期項分量[2,4]和3 個波動項分量[2,1]～[2,3]，見圖2。周期項分量可有效反映原始年降水量數據的周期性，波動項分量可反映原始年降水量數據的震蕩變化。本文選取1956-2015 年降水量時序數據作為訓練樣本，2016-2021 年降水量時序數據作為預測樣本。

圖2 文山州年降水量時序數據分解Fig.2 Decomposition of Annual Precipitation Time Series Data in Wenshan Prefecture

步驟二：采用Cao 方法[19-21]確定圖1 中周期項分量[2,4]和波動項分量[2,1]～[2,3]的輸入步長a，并利用前a年降水量預測當年降水量，即輸入層節點數為a，輸出層節點數為1。經Cao方法計算，圖1中周期項分量[2,4]和波動項分量[2,1]～[2,3]的輸入步長a分別為14、13、10、11。

因此，WPT-POA-GPR 等模型的輸入、輸出可表述為：

式中：M為樣本數量；a為輸入步長。

步驟三：利用周期項、波動項分量訓練樣本的預測值與實際值構建均方誤差（MSE）作為優化GPR 超參數的適應度函數：

步驟四：設置POA、SCSO、CO種群規模為30，最大迭代次數為100，其他采用算法默認值。初始化蜣螂、珍鲹、獵豹個體位置。

本文基于二次有理式協方差函數（RQ）構建GPR 模型，GPR 核函數、噪聲方差、協方差搜索空間分別設置為[0.01,10]、[0.01,10]、[0.1,100]；SVM 選擇RBF 徑向基核函數，懲罰因子、核函數參數、不敏感損失系數搜索空間分別設置為[0.01,100]、[0.01,100]、[0.000 1,0.1]；RBF 擴展速度搜索空間設置為[0.001,1 000]。為驗證優化效果，WPT-GPR模型超參數采用試算法確定；所有模型的輸入數據均采用[0,1]進行歸一化處理。

步驟五：基于式（20）計算孔雀、沙貓、獵豹個體適應度值，找到并保存最佳個體位置。令t= 1。

步驟六：分別利用上述POA、SCSO、CO 算法位置更新算子更新個體新位置。

步驟七：利用位置更新后的孔雀、沙貓、獵豹個體計算適應度值，比較并保存當前最佳個體位置。

步驟八：重復步驟六～步驟八直至滿足算法最大迭代次數。

步驟九：輸出最佳個體位置，該位置即為GPR 最佳超參數向量。利用該向量建立WPT-POA-GPR 等模型對周期項分量[2,4]和波動項分量[2,1]～[2,3]進行預測和加和重構。

步驟十：利用平均絕對百分比誤差MAPE、平均絕對誤差MAE、均方根誤差RMSE和決定系數DC對模型進行評價。其中，MAPE、MAE用于反映模型預測誤差大小，RMSE用于衡量觀測值同真值之間的偏差，MAPE、MAE、RMSE越小，說明模型性能越優，預測精度越高；DC反映變量之間相關關系的密切程度，其值越大，說明模型性能越好。

2 結果與分析

構建WPT-POA-GPR 等10 種模型對實例年降水量進行訓練及預測，結果見圖3 和圖4；預測相對誤差、絕對誤差見圖5。

圖3 年降水量訓練誤差Fig.3 Annual precipitation training error

圖4 年降水量預測誤差Fig.4 Annual precipitation prediction error

圖5 年降水量訓練-預測誤差效果對比圖（后6年為預測樣本）Fig.5 Comparison Chart of Annual Precipitation Training and Prediction Error Effects (the next 6 years are prediction samples)

依據圖3～圖5可以得出以下結論：

（1） WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR模型對實例年降水量預測的MAPE在0.46%～0.52%、MAE在5.25～5.80 mm、RMSE在7.72～8.20 mm 之間，確定性系數DC≧0.996 7，預測精度優于WPT-POA-RBF、WPT-SCSO-RBF、WPT-CO-RBF、WPT-GPR 模型，遠優化WPT-POA-SVM、WPT-SCSO-SVM、WPT-CO-SVM 模型；模型對實例年降水量訓練精度同樣優化其他7 種模型。可見，WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型具有更高的預測精度和更好的泛化能力，將其用于年降水量時間序列預測是可行的。

（2） 與WPT-POA-RBF、WPT-SCSO-RBF、WPT-CORBF 模型和WPT-POA-SVM、WPT-SCSO-SVM、WPT-COSVM 模型相比，WPT-POA-GPR、WPT-SCSO-GPR、WPTCO-GPR 模型預測的MAPE分別提高64.9%、76.7%以上，表明GPR 在處理高維度、小樣本、非線性等復雜回歸問題中表現出色，具有較好的擬合精度及泛化性能。

（3） 與WPT-GPR 模型相比，WPT-POA-GPR、WPTSCSO-GPR、WPT-CO-GPR 模型預測的MAPE 提高72.8%以上，表明POA、SCSO、CO能有效優化GPR超參數，提高GPR預測性能。

3 討 論

提高年降水量預測精度是水文預報研究的熱點和難點。然而，由于年降水量影響因素眾多，往往呈現出非線性、多尺度等特征，傳統單一模型難以獲得較好的預測精度。當前，基于“分解算法+預測模型”的預測方法已在年降水量預測研究中得到應用，并取得較好的預測效果，但也存在一定的不足，如文獻[14]建立的模型中，WD分解方法存在分解不徹底，缺乏對高頻數據的分析，LSTM 存在計算量大、耗時長等不足；文獻[15]中，VMD 分解方法存在分解數值K選取困難，TCN 存在運行速度慢、超參數選取依靠人工調試等不足；文獻[16]中，CEEMDAN 分解方法存在計算量大、復雜度高等問題；文獻[17]中，EEMD 分解方法存在人為經驗添加噪聲、分解分量過多等不足；文獻[18]中，SSA 分解方法存在分解過程復雜，分解分量多、計算規模大等不足，SVR 存在超參數選取困難、處理大規模數據時效果不佳等不足。

本文融合了WPT 分解方法、POA/SCSO/CO 3 種新型群體智能算法和GPR 預測器，建立了WPT-POA-GPR、WPTSCSO-GPR、WPT-CO-GP 年降水量時間序列預測模型，并構建 WPT-POA-SVM、WPT-SCSO-SVM、WPT-CO-SVM、WPT-POA-RBF、WPT-SCSO-RBF、WPT-CO-RBF、WPTGPR 對比分析模型，通過文山州年降水量預測實例驗證了所構建的3種模型具有較好的預測效果，模型及方法可在類似降水量預測研究中進一步研究與推廣。今后可在以下方面作進一步研究：

（1） 通過實例對比驗證WPT 與上述WD、VMD、CEEMDAN、EEMD、SSA分解方法的分解效果。

（2）通過與傳統粒子群優化算法、遺傳算法優化GPR 超參數的對比，進一步驗證POA/SCSO/CO 算法在調優GPR 超參數中的優勢。

（3）本文僅以年降水量數據作為模型輸入，未考慮氣象、地形等因素，因此，該模型及預測精度仍有進一步提升的空間。

4 結 論

為提高年降水量預測精度，本文融合了WPT 方法、POA/SCSO/CO 三種新型群體智能算法和GPR 模型，提出WPTPOA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 預測模型，并構建若干對比模型，結合具體算例，主要結論有：

（1） WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR模型預測誤差小于其他對比模型，具有較高的預測精度和較好的泛化能力，將其用于年降水量時間序列預測是可行的。模型及方法可為相關降水量時間序列預測研究提供參考。

（2）針對GRP 在超參數尋優時依賴參數初始值、容易陷入局部最優等問題，利用POA、SCSO、CO 優化GRP 超參數，不但改善了GPR 預測性能，而且有效提升了模型的智能化水平，具有較好的適用性。

（3）針對非線性、非平穩較強的年降水量時間序列，利用WPT 對其進行平穩化處理，可得到更具規律的周期項分量和波動項分量，顯著提高了年降水量的預測精度。