數據分解技術與若干智能算法優化的高斯過程回歸總氮預測

王永順,崔東文

(1.云南省水文水資源局文山分局,云南 文山 663000;2.云南省文山州水務局,云南 文山 663000)

總氮(TN)是指可溶性以及懸浮物中所含氮的總量,是反映水體污染程度和衡量湖庫營養化狀態的重要指標之一,其濃度超標時,將導致水體中微生物大量繁殖,浮游生物生長旺盛,出現水體富營養化狀態,嚴重威脅到居民生活、生產用水安全和社會安定。提高TN時間序列預測精度,對于保障飲水安全、富營養防治、水環境保護治理等具有重要意義。近年來,隨著機器學習與智能算法的快速發展,基于數據驅動的非機理組合預測模型逐漸興起,已在TN預測研究中得到應用,如BP神經網絡[1-2]、極限學習機(ELM)[3]、支持向量回歸機(SVR)[4]、廣義回歸神經網絡(GRNN)[5]、AdaBoost機器學習[6]、循環神經網絡(RNN)[7]、隨機森林(RF)[8]等。高斯過程回歸(Gaussian Process Regression,GPR)是基于貝葉斯理論與統計學習理論的非參數監督學習方法,適于處理高維數、非線性等復雜的回歸問題,已在時間序列預測分析、動態系統模型辨識等多個領域得到廣泛應用,但鮮見于TN時間序列預測研究。在實際應用中,GPR預測效果對超參數取值的依賴程度較高,目前主要采用實驗試湊、經驗選擇等方法對GPR超參數進行選取,存在隨機性大、容易陷入局部最優等問題。為改善其不足,粒子群優化(PSO)算法[9-10]、北方蒼鷹優化(NGO)算法[11]、鯨魚優化算法(WOA)[12]、人工蟻群優化(ACO)算法[13]等群體智能算法償試用于GPR超參數優化,并取得較好的優化效果。

TN時間序列廣泛存在于水質預測領域,主要根據已有的TN歷史數據對未來進行預測,但由于TN時間序列具有隨機因素影響大、非線性、非平穩性和多尺度等特征,傳統單一預測方法難以滿足TN的預測精度要求。近年來,基于“分解-預測-重構”的預測模型已廣泛應用于各行業領域的時間序列預測,并取得較好的預測效果。其中,分解技術是決定此類模型預測精度高低的關鍵,當前的主要分解技術有經驗模態分解(EMD)、集合經驗模態分解(EEMD)、互補集合經驗模態分解(CEEMD)、變分模態分解(VMD)、小波變換(WT)、奇異譜分解(SSA)等。但上述分解方法存在不足:EMD易出現模態混疊現象,導致分解效果不理想;EEMD通過加入白噪聲改善了模態混疊現象,但乃存在重構信號噪聲殘留問題;CEEMD在一定程度上抑制了模態混疊問題,但出現了偽分量,降低了模型的預測精度;VMD可以有效緩解信號的模態混疊,但存在分解模態數量不確定性問題;WT只對低頻信號再次分解,而忽視了高頻信號,易導致分解不徹底,影響預測精度;SSA存在窗口長度難以確定以及分解分量過多等不足。

為有效提高TN時間序列預測精度,改進GPR預測性能,本文分別基于經驗小波變換(Empirical Wavelet Transform,EWT)、小波包變換(Wavelet Packet Transform,WPT)分解技術,提出魚鷹優化算法(Osprey Optimization Algorithm,OOA)、霧凇優化算法(Rime Optimization Algorithm,ROA)、禿鷹搜索(Bald Eagle Search,BES)算法、黑寡婦優化算法(Black Widow Optimization Algorithm,BWOA)優化GPR的EWT-OOA-GPR、EWT-ROA-GPR、EWT-BES-GPR、EWT-BWOA-GPR、WPT-OOA-GPR、WPT-ROA-GPR、WPT-BES-GPR、WPT-BWOA-GPR預測模型(簡稱EWT-OOA-GPR等8種模型),并構建基于WT分解的WT-OOA-GPR、WT-ROA-GPR、WT-BES-GPR、WT-BWOA-GPR模型,基于支持向量機(SVM)的EWT-OOA-SVM、EWT-ROA-SVM、EWT-BES-SVM、EWT-BWOA-SVM、WPT-OOA-SVM、WPT-ROA-SVM、WPT-BES-SVM、WPT-BWOA-SVM模型,未經優化的EWT-GPR、WPT-GPR模型和未經分解的OOA-GPR、ROA-GPR、BES-GPR、BWOA-GPR模型作對比分析模型。通過全國重要飲用水水源地暮底河水庫2008年1月至2022年12月月監測TN濃度時序數據對各模型進行檢驗,旨在驗證EWT-OOA-GPR等8種模型用于TN預測的可行性。

1 數據來源與研究方法

1.1 數據來源

暮底河水庫位于云南省文山市,是全國重要飲用水水源地之一,屬于紅河流域瀘江水系。水庫設計壩高67.6 m,壩頂高程1 340.10 m,庫容5 784.9萬m3,正常蓄水位1 338.00 m,是一座以防洪、城市供水為主,兼顧灌溉、發電及改善生態環境用水的綜合性多功能中型水庫,水庫年供水能力1.01億m3,供水人口30余萬人。近年來,受徑流區居民生活、生產影響,水庫面臨農業面源、畜禽養殖以及生活點源等污染威脅,水庫TN時有超標,對生產生活造成嚴重影響。因此,開展暮底河水庫TN預測研究對于水源地的保護治理具有重要意義。

本文數據采用暮底河水庫2008年1月至2022年12月TN月尺度監測數據,其過程線見圖1—3原始序列。從原始序列可以看出,暮底河水庫TN月監測濃度序列波動性較大,復雜程度較高,呈現出較強的非線性和非平穩性,不利于直接預測。本文選取2008年1月至2022年12月暮底河水庫月監測TN濃度時序數據的70%作為訓練樣本,剩余30%作為預測樣本。

圖1 EWT分解

1.2 研究方法

1.2.1經驗小波變換(EWT)

EWT是Gilles于2013年提出的非平穩信號處理方法,它融合了EMD的自適應分解理念和WT理論的緊支撐框架,為信號處理提供了一種全新的自適應時頻分析思路。相比EMD,EWT能夠自適應選擇頻帶,完美地解決因信號時頻尺度不連續引發的模態混疊問題,同時EWT具備可靠的數學理論基礎,計算復雜度低,還能夠克服EMD分解中“過包絡”和“欠包絡”問題。EWT分解原理簡述如下[14-16]。

a)對分析信號f(t)進行傅里葉變換,求得傅立葉頻譜F(ω)。

b)將頻譜范圍歸一化到ω∈[0,π],并分割成N個連續區間Λn和N-1個分界頻率。

(1)

(2)

(3)

(4)

e)原始信號f(t)通過經驗小波變換最終被分解為若干經驗模態分量,描述為:

(5)

(6)

1.2.2小波包變換(WPT)

WPT能同時對信號低頻部分和高頻部分進行分解,更適用于TN時間序列分解。WPT對TN原始信號進行分解的公式[17-18]為:

(7)

重構算法為:

(8)

1.2.3魚鷹優化算法(OOA)

OOA是Trojovsky P等[19]于2023年受自然界魚鷹捕魚行為啟發而提出一種新型元啟發式優化算法。該算法模擬了魚鷹在海上識別魚群位置、獵捕魚兒(勘探)和將魚搬運安全位置(開發)2種行為,并通過對2種智能行為進行建模實現位置更新來求解待優化問題。OOA數學描述如下。

a)初始化。與其他元啟發式優化算法類似,OOA也是從種群初始化開始。在m維搜索空間中,種群規模為N的魚鷹初始化位置描述為:

xi,j=lbj+r·(ubj-lbj),i=1,2,…,N;j=1,2,…,m

(9)

式中xi,j——第i只魚鷹第j維空間位置;ubj、lbj——優化問題的上、下限值;r——介于0和1之間的隨機數;N——魚鷹種群規模;m——優化問題維度。

b)位置識別和獵捕魚兒(勘探階段)。魚鷹具有強大的視覺能力,能夠探測水下魚兒位置。在確定水下魚兒位置后,魚鷹潛入水下攻擊并捕獵魚兒。在OOA中,魚鷹隨機選擇其中一條魚兒位置并攻擊它。這一階段魚鷹位置更新描述如下:

(10)

c)將魚搬運至安全位置(開采階段)。獵獲一條魚后,魚鷹會將它運至一個安全的位置,并在那里進食。將魚運至安全位置的行為將導致魚鷹在搜索空間中的位置變化,這一行為提升了OOA在局部搜索中的開采能力。對于每只魚鷹,OOA通過選擇新的隨機位置作為“安全吃魚位置”。這一階段魚鷹位置更新描述為:

(11)

1.2.4霧凇優化算法(ROA)

ROA是Su H等[20]于2023年受自然界霧凇冰物理現象啟發而提出一種新型元啟發式優化算法。該算法模擬了霧凇冰的軟霧凇和硬霧凇生長過程,并通過對軟霧凇搜索策略、硬霧凇穿刺機制、貪婪選擇機制進行數學建模實現位置更新來求解待優化問題。ROA數學描述如下[20]。

a)初始化。與其他元啟發式優化算法類似,ROA也是從種群初始化開始。設在d維搜索空間中,種群規模為R的霧凇粒子初始化位置描述為:

xij=Lbij+r·(Ubij-Lbij),i=1,2,…,R;j=1,2,…,d

(12)

式中xij——第i個霧凇粒子第j維空間位置;Ubij、Lbij——優化問題的上、下限值;r——介于0和1之間的隨機數;R——霧凇種群規模;d——優化問題維度。

b)軟霧凇搜索策略。在微風環境中,軟霧凇的生長具有很強的隨機性,霧凇粒子可以自由地覆蓋物體的大部分表面,且在同一方向上生長緩慢。受軟霧凇生長的啟發,ROA提出軟霧凇搜索策略,該策略使得ROA能夠在早期迭代中快速覆蓋整個搜索空間,且不易陷入局部最優。軟霧凇粒子位置更新描述如下:

(13)

c)硬霧凇穿刺機制。在強風條件下,硬霧凇生長比軟霧凇更簡單、更規律。受霧凇穿刺生長的啟發,ROA提出硬霧凇穿刺機制。該機制用于實現穿刺霧凇粒子的位置更新,從而提高算法的收斂性,跳出局部極值尋優。硬霧凇粒子位置更新描述如下:

(14)

d)貪婪選擇機制。ROA提出了一種積極的貪婪選擇機制來參與種群更新,以提高全局搜索效率,即利用更新后的霧凇粒子適應度值與更新前霧凇粒子的適應度值進行比較,如果更新后的適應度值比更新前更好,則替換霧凇粒子更新前位置,否則保留原霧凇粒子位置。

1.2.5禿鷹搜索(BES)算法

BES是馬來西亞學者Alsattar等[21]于2019年提出的一種新型群體智能仿生算法。該算法模擬了禿鷹在狩獵過程中選擇搜索空間、搜索空間獵物和俯沖捕獲獵物3個階段,并通過迭代更新禿鷹捕獲獵物位置,即待優化問題最優解。BES算法描述簡述如下[21-22]。

a)選擇搜索空間階段。在該階段,禿鷹隨機選擇搜索區域,并通過判斷獵物數量來確定最佳搜索位置。該階段數學描述為:

Pi,new=Pbest+α×r×(Pmean-Pi)

(15)

式中Pi,new——當前禿鷹搜索位置;α——控制位置變化參數,取值范圍[1.5,2];r——[0,1]范圍內隨機數;Pbest——當前禿鷹最佳搜索位置;Pmean——先前搜索結束后禿鷹的平均分布位置;Pi——第i只禿鷹搜索位置。

b)搜索空間獵物階段。在該階段,禿鷹在選定搜索空間內以螺旋形狀飛行搜索獵物,加速搜索進程,尋找最佳俯沖捕獲位置。該階段數學描述為:

Pi,new=Pi+y(i)×(Pi-Pi+1)+x(i)×(Pi-Pmean)

(16)

c)俯沖捕獲獵物階段。在俯沖階段,禿鷹從搜索空間最佳位置快速俯沖飛向目標獵物,禿鷹種群其他個體也同時向最佳位置移動攻擊獵物。該階段數學描述為:

Pi,new=rand·Pbest+x1(i)×(Pi-c1·Pmean)+y1(i)×(Pi-c2·Pbest)

(17)

1.2.6黑寡婦優化算法(BWOA)

BWOA是Pea-Delgado A F等人于2020年通過模仿黑寡婦蜘蛛不同的求偶移動策略而提出的一種新型群體優化算法。BWOA優化原理簡述如下[23]:

a)移動策略。黑寡婦蜘蛛通過線性和螺旋形方式實現在蛛網中的移動,數學描述如下:

(18)

b)性信息素。性信息素在黑寡婦蜘蛛求偶交配中起著重要作用。雄性蜘蛛對來自營養充足的雌性蜘蛛的性信息素反應更靈敏,而對饑餓和缺乏食物的雌性蜘蛛的性信息素反應遲鈍。在BWOA中,雌性黑寡婦蜘蛛的性信息素描述如下:

(19)

式中 pheromone(i)——第i只雌性黑寡婦蜘蛛性信息素;fitnessmax、fitnessmin——當前最佳和最差適應度值;fitness(i)——第i只蜘蛛適應度值。

c)取代策略。在BWOA中,具有低性信息素的雌性蜘蛛代表同類相食類雌性蜘蛛,此類蜘蛛不受雄性蜘蛛的青睞;如果此類蜘蛛存在,其將被另一只蜘蛛取代:

(20)

1.2.7高斯過程回歸(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′)決定,因此高斯過程可定義為[9-13]:

f(x)～GP[m(x),k(x,x′)]

(21)

考慮噪聲后,GPR回歸模型可表示為:

yi=f(xi)+ε

(22)

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

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

(23)

式中K(X,X)——n×n階協方差矩陣;In——單位矩陣。

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

(24)

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

(25)

1.3 建模流程

步驟一分別利用EWT、WPT對暮底河水庫月監測TN濃度時序數據進行分解,分別得到6個和4個模態分量E1—E6、W1—W4,其中E1、W1分量為波動項,反映了原始時序數據的隨機變化情況,是TN預測精度提升的關鍵分量;E2—E5、W2—W3為周期項,反映了原始時序數據的微小震蕩變化;E6、W4為趨勢項,大致反映了原始時序數據的趨勢性,見圖1、2。作為對比分解方法,利用WT對TN濃度時序數據進行分解,得到6個模態分量D1—D6,見圖3。

步驟二采用文獻[17,24]中Cao方法確定圖1—3中各模態分量的輸入步長a,并利用前a月的TN分解分量來預測未來1月的分量值,即輸入層節點數為a,輸出層節點數為1。經計算,圖1模態分量E1—E6的輸入步長a分別為11、15、18、16、8、17;圖2模態分量W1—W4的輸入步長a分別為19、23、16、21;圖3模態分量D1—D6的輸入步長a分別為13、22、18、23、15、13;利用同樣的方法確定TN原始序列的輸入步長a為12。

圖2 WPT分解

圖3 WT分解

因此,預測模型的輸入、輸出可表述為:

(26)

式中M——樣本數量;a——輸入步長。

步驟三利用各模態分量訓練樣本的擬合值與實際值構建均方誤差(MSE)作為OOA、ROA、BES、BWOA優化GPR超參數的適應度函數:

(27)

步驟四設置OOA、ROA、BES、BWOA種群規模為30,最大迭代次數為100。隨機初始化魚鷹、霧凇粒子、禿鷹和蜘蛛個體空間位置。

GPR、SVM參數設置如下:基于二次有理式協方差函數(RQ)構建GPR模型,設置GPR核函數、噪聲方差、協方差的搜索空間為[0.01,100];SVM選擇高斯核函數,懲罰系數、核函數參數搜索范圍均設置為[0.01,100],不敏感損失系數設置為0.001;為驗證優化效果,WPT-GPR模型超參數采用試算法確定;所有模型的輸入數據均采用[0,1]進行歸一化處理。

步驟六魚鷹、霧凇粒子、禿鷹和蜘蛛個體的位置更新。

a)魚鷹位置更新。通過隨機確定所選魚兒位置,計算適應度值FP1,若FP1若優于F,則更新魚鷹位置并保存為最佳魚鷹位置XBest,否則保留原魚鷹最佳位置XBest;通過選擇新的隨機位置作為“安全吃魚位置”,并計算適應度值FP2,若FP2若優于F,則更新魚鷹位置并保存為最佳魚鷹位置XBest,否則保留原魚鷹最佳位置XBest。

b)霧凇粒子位置更新。計算附著系數E和隨機數r2,若r2>E,則更新軟霧凇粒子位置;計算隨機數r3,若r3

c)禿鷹位置更新。分別執行禿鷹隨機搜索區域、搜索空間獵物、俯沖捕獲獵物3個位置更新算子,更新禿鷹位置。

d)蜘蛛位置更新。生成隨機數rand(),若rand()<0.3,則通過線性方式更新蜘蛛位置,否則,通過螺旋形方式更新蜘蛛位置;計算每只雌性蜘蛛的信息素,并更新具有低性信息素的雌性蜘蛛位置。

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

步驟十選取平均絕對百分比誤差MAPE、平均絕對誤差MAE、均方根誤差RMSE、決定系數DC對模型進行評估。

(28)

2 預測結果及對比

構建EWT-OOA-GPR等22種模型對暮底河水庫TN分解分量進行訓練及預測,將預測結果加和重構后得到最終預測結果,結果見表1,預測誤差見圖4;同時構建未經分解的OOA-GPR、ROA-GPR、BES-GPR、BWOA-GPR模型對原始TN序列進行訓練及預測,結果見表1。

表1 暮底河水庫TN預測結果及對比

圖4 暮底河水庫TN預測相對誤差對比

依據表1、圖4可以得出以下結論。

a)EWT-OOA-GPR等8種模型對暮底河水庫TN預測的MAPE、MAE、RMSE分別在0.161%～0.219%、0.001 5～0.002 0 mg/L、0.001 9～0.002 4 mg/L,DC均為0.999 9,預測效果優于EWT-GPR、WPT-GPR模型,遠優于WT-OOA-GPR、WT-ROA-GPR、WT-BES-GPR、WT-BWOA-GPR、EWT-OOA-SVM、EWT-ROA-SVM、EWT-BES-SVM、EWT-BWOA-SVM、WPT-OOA-SVM、WPT-ROA-SVM、WPT-BES-SVM、WPT-BWOA-SVM模型,大大優于OOA-GPR、ROA-GPR、BES-GPR、BWOA-GPR模型,具有更高的預測精度和更好的泛化能力,將其用于水庫TN預測是可行的。其中,EWT-OOA-GPR、EWT-ROA-GPR、EWT-BES-GPR、EWT-BWOA-GPR模型的預測效果要優于WPT-OOA-GPR、WPT-ROA-GPR、WPT-BES-GPR、WPT-BWOA-GPR模型。

b)與EWT-GPR、WPT-GPR相比,EWT-OOA-GPR等8種模型預測的MAPE提高了42.7%以上,表明OOA、ROA、BES、BWOA能有效優化GPR超參數,提高模型預測性能,通過人工試算調試GPR超參數繁瑣、費時。

c)與EWT-OOA-SVM等8種模型相比,EWT-OOA-GPR等8種模型預測的MAPE提高了88.5%以上,表明在相同分解及優化情形下,GPR在處理高維度、小樣本、非線性等復雜回歸問題中表現出色。

d)與WT-OOA-GPR、WT-ROA-GPR、WT-BES-GPR、WT-BWOA-GPR相比,EWT-OOA-GPR等8種模型預測的MAPE提高了84.2%以上,表明EWT、WPT能將總氮時間序列分解為更具規模的分量,分解效果優于WT方法。

e)與OOA-GPR、ROA-GPR、BES-GPR、BWOA-GPR相比,EWT-EVO-GPR、EWT-CDO-GPR模型預測的MAPE提高了98.9%以上,表明未經分解直接采用OOA-GPR、ROA-GPR、BES-GPR、BWOA-GPR模型進行預測,雖然計算規模小、復雜度低,但由于原始TN時序數據具有較強的非線性和非平穩性特征,預測誤差最大,顯然不能滿足TN預測精度需求。

3 結論

基于EWT、WPT兩種時序數據處理技術,分別提出OOA、ROA、BES、BWOA算法與GPR相融合的TN時間序列預測模型,并構建WT-OOA-GPR等18種對比模型,以暮底河水庫2008年1月至2022年12月月監測TN濃度數據作為實例進行驗證,得到以下結論。

a)EWT-OOA-GPR等8種模型的預測效果優于其他18種對比模型,具有更高的預測精度,將其用于水庫TN預測是可行的。其中,基于EWT分解的EWT-OOA-GPR等4種模型的預測精度要高于基于WPT分解的WPT-OOA-GPR等4種模型。

b)EWT融合了WT以及EMD優勢,可以通過時間序列頻譜自適應地確定分解濾波器,能夠完美地解決分解頻率的混疊問題;WPT能同時對低頻、高頻信號再次分解,且對突變信號分解效果更好,二者的分解效果均優于WT方法。

c)OOA、ROA、BES、BWOA能有效優化GPR超參數,克服了人工試算調試GPR超參數繁鎖等不足的同時,有效提高了GPR預測性能。

d)在相同分解及優化情形下,GPR在處理高維度、小樣本、非線性等復雜回歸問題中表現較好,預測性能優于SVM。

e)由于EWT-OOA-GPR等8種模型僅利用了暮底河水庫月監測TN的歷史數據作為輸入,未考慮氨氮、硝酸鹽、亞硝酸鹽以及降水等因素的影響,因此,該模型及方法仍有進一步改進的空間。