基于CSA-SIM-LSSVM的不同時間尺度參考作物蒸散發估算研究

王文川， 趙釗， 張磊

(1.華北水利水電大學 水資源學院，河南 鄭州 450046； 2.華北水利水電大學 水利學院，河南 鄭州 450046)

自2012年我國實行最嚴格水資源管理制度以來，如何更高效地管理、利用好水資源一直是熱門問題[1]。農業用水作為用水大戶，其中針對作物進行需水量精準把握是決策者的必修功課。參考作物蒸散發(Reference Evapotranspiration,ET0)是計算作物需水量的關鍵環節，以往該值的估算大多依賴于各種基于實驗的經驗公式，如Blaney-Criddle公式[2]、Irmak-Allen公式[3]、Ritchie模型[4]、Priestley-Taylor公式[5]和Penman-Monteith公式[6]。其中的Penman-Monteith公式在全球范圍內得到了廣泛應用，且能在不同類型地區取得較好的計算效果[7]。但該公式需要較為全面的氣象參數，而相關參數在有些區域是難以獲得的[8]。隨著人工智能和機器學習的發展，諸多新型技術得以在ET0估算上應用并取得了一定的效果。侯志強等[9]將ET0視作為各種氣象因素影響下的復雜非線性問題，利用最小二乘支持向量機(Least Square SVM,LSSVM)對河套地區日ET0進行了模擬研究。該研究顯示，在同等輸入氣象因子條件下，LSSVM的表現好于傳統的Priestley-Taylor公式和Hargreaves公式。但是LSSVM模型的性能高度依賴其控制參數，在此基礎上，一些智能優化算法展示了較高的客觀性和效率。鞠彬和王嘉毅[10]利用粒子群算法(Particle Swarm Optimization,PSO)與LSSVM耦合進行月尺度ET0的預測模擬，但是在全氣象要素輸入情況下其精度仍然不高，平均相對誤差高達13.52%。張育斌等[11]則采用耦合模擬退火算法(Coupled Simulated Annealing，CSA)優化LSSVM模型參數形成CSA-LSSVM模型，該模型的日尺度ET0預測的平均均方根誤差RMSE為0.39 mm/d，精度得到較高提升。與此同時，其他類型的機器學習模型諸如極限學習機[12]、極限梯度提升算法[13]、高斯回歸過程[14]等在ET0估算模擬領域也取得了一定的研究成果。

綜合上述研究成果可以看出，模型自身的結構對ET0模擬精度起著基礎性作用，對于模型參數的調整也是至關重要的一步。上述文獻中關于模型調整參數全局優化的PSO算法和CSA算法能力有限，十分有必要進一步平衡算法以提高全局尋優和局部挖掘能力。鑒于此，本文采用LSSVM為基礎預報模型，從優化算法的角度首先采用CSA算法為全局探索算法，利用多個起始點進入第一階段的尋優確定多個可能極值點，第二階段再利用擅長局部挖掘的單純形法進一步提供更可靠的候選解。此外，針對日尺度和月尺度ET0估算工作的不同特性，嘗試引入其他新興技術和特殊農時知識，例如LASSO回歸法篩選因子[15]、小波包分解降噪技術[16]和旬序數[17]。以此降低ET0估算對氣象數據的要求和進一步提高ET0估算的模擬精度，為水資源管理的決策提供更多的方法參考。

1 基本原理

1.1 耦合模擬退火-單純形(CSA-SIM)組合算法

1.1.1 耦合模擬退火算法

經典模擬退火算法(Simulated Annealing,SA)在尋優過程中引入了Metropolis準則，即以一定的概率接受較差目標函數值，從而避免算法掉入局部最優。伴隨著溫度的不斷收縮，SA尋優半徑逐漸減小趨于0以達到全局最優。但是，該算法在使用過程中受初始溫度和溫度收縮規則的影響，往往尋優精度不高。由此，Souza等提出了耦合模擬退火算法(CSA)[18]。不同于傳統的模擬退火——每次只生成一個解的尋優過程，CSA在解決方案空間的探測過程中引入了多次啟動初始化以提供更多的先驗信息，在各個退火過程中利用一個耦合能量(代表接受概率)和溫度等信息的公式實現信息的交互以提高多樣性，這種改進策略在一定程度上減少了尋求全局最優過程的迭代計算量。耦合模擬退火算法主要包括搜索向量的生成和接受兩個主要過程，公式表達如下:

(1)

式中：x為啟動器;Ω為搜索域;θ為當前超參數的集合;n為CSA初始啟動器數量;rand為[0,1]均勻分布上產生的隨機數;Tk為k時刻的溫度;k為迭代次數。因此，y的分布半徑會隨著Tk越大而更大，反之變小，最終值也會逐漸收斂于0。本文采用下式的溫度更新規則:

(2)

CSA算法可看作為多個SA算法的耦合算法。對最小化問題而言，單個啟動器的候選解的接受概率函數令為A(xi→yi)，即:

(3)

式中：xi為當前狀態;yi為新狀態;Tac為當前接受溫度，其更新規則參閱文獻[18];E指整個系統的能量。

利用Boltzmann參數求解僅有兩個狀態(i=1、2)的系統在第i個狀態的概率Pi:

(4)

式中：kB為Boltzmann常數；Ei為i狀態下當前系統的能量值；T為i狀態下的溫度；Z為當前系統所有狀態的能量和。

(5)

在狀態yi和溫度T已給定的情況下，狀態yi被接受的概率值由式(6)近似表示。為了實現耦合模擬退火，先初始化一個多狀態系統，x為狀態的集合，xi為當前的第i個狀態，yti為第i個當前狀態將要轉移的新狀態。設x∈(x1,x2,…,xn)，式(3)轉換為式(6):

A(xi→yti)=

(6)

此時，當前狀態x∈(x1,x2,…,xn)對應新狀態yt∈(yt1,yt2,…,ytn)的接受概率為A(x→yt)∈[A(x1→yt1)，A(x2→yt2),…,A(xn→ytn)]。狀態集合x內各個狀態接受對應的轉移狀態yt的概率，除了考慮自身外，還要考慮其它狀態的耦合。特殊情況下，當狀態總數n=1時，方法將退化為傳統的模擬退火求解問題。因此，CSA在接受精度高的解的同時能夠接受擬合較差的候選解，并以此實現信息的共享，促進增強單啟動SA算法。詳細的求解過程可參閱文獻[18]。

1.1.2 單純形法

單純形算法(Simplex Algorithm)是Nelder和Mead提出的一種擅長局部優化的方法[19]。單純形是具有n+1個頂點的幾何圖形，算法步驟主要有初始、準備、反射、延伸、收縮、棱長減半和程序出口。對于給定經度ε的最小化問題，單純形法過程簡便，詳細的求解步驟可參閱文獻[19]。

1.1.3 耦合模擬退火—單純形優化算法

首先利用耦合模擬退火算法進行全局尋優以獲得質量較好的數個局部最優點，再將其結果當作單純形優化的初始點進行局部挖掘，進而減少計算資源，提高效率。

1.2 最小二乘支持向量機(LSSVM)

最小二乘支持向量機(Least Squares Support Vector Machine,LSSVM)是一種在經典機器學習模型——支持向量機基礎上發展的改進算法[20]。原始支持向量機(Support Vector Machine,SVM)可用來解決分類和回歸問題[21]。對于回歸問題，SVM的基本思想是通過核函數把低維空間中的非線性問題變換到高維特征空間中的線性回歸，對應地在這個高維特征空間中尋求一個最優超平面或最優超曲面來尋求問題的解。假設某訓練樣本的集合為(xi,yi),i=1、2、…、N;xi∈RM,yi∈R，其中xi為輸入向量，yi為對應xi的輸出向量，N為樣本個數，M為維數。設經過變換后的非線性回歸函數為f(x)=+b，其中w為系數向量，b為偏置常量，Φ為非線性映射函數，為w與Φ(x)的內積。根據結構風險化最小化原則，通過引入相關損失函數和松弛變量，SVM的擬合過程將轉化為求解一個凸二次優化問題。LSSVM則是將該優化問題中的不等式約束改為等式約束，利用解方程組的方法避開了復雜耗時的二次規劃問題，LSSVM構造以下最優化目標函數:

(7)

式中:γ為正則化參數;δi為第i個松馳變量。

引入拉格朗日因子αi和KKT條件， Lagrange函數L可寫為:

(8)

分別對w、b、δ和α求偏微分可得:

(9)

消去(9)中的w和δ,可得:

(10)

式中:E=(1,1,…,1)T；Z=[Φ(x1),Φ(x2),…,Φ(xN)]T；A=(α1,α2,…,αN)T；Y=(y1,y2,…,yN)T；I為單位矩陣。經高維空間映射后計算量巨大，選用滿足Mercer條件的核函數K(x,xi)=<Φ(x),Φ(xi)>，那么LSSVM模型的回歸方程最終確定為:

(11)

由以上可以看出，合適的核函數是保障LSSVM模型的預測精度的必要條件。常用的核函數有:

(12)

式中：g、r、d均為核函數的參數，即需要調整的參數。

一般而言，線性核函數適用于線性可分的情況，比較適合特征數量很多，與樣本數相差不大的情形;多項式核函數和徑向基核函數均能將低維問題映射到高維特征空間，由前述公式不難看出，多項式核函數參數更多。對于參考作物蒸散發預測問題，特征數與樣本數均不大，綜合考慮，本文采用徑向基核函數。

確定了核函數之后，由公式(11)可知LSSVM模型最終的回歸方程與正則化參數γ和核函數參數g有關，也即算法欲優化的目標函數為:

error=f(x|γ,g)。

(13)

式中error為在參考作物蒸散發預測模擬問題中的各種誤差形式。由公式(13)知，在確定的輸入訓練矩陣x下，優化算法的目的在于尋找合適的γ和g使得error最小。

2 CSA-SIM-LSSVM估算模型實現

按照上述的CSA-SIM算法以及LSSVM模型的原理，本文提出的CSA-SIM-LSSVM模型旨在:搜尋一組最優的向量(γ,g)，在模型訓練過程中采用十折交叉檢驗法[22]，使得模型訓練代表誤差的目標函數值最小。針對確定的實驗樣本、輸入因子集合和輸出，CSA-SIM-LSSVM的實現流程如圖1所示。

圖1 CSA-SIM-LSSVM模型流程圖

3 應用實例

3.1 研究區域和數據特征

研究區域為華北平原(112°30′E～119°30′E,34°46′N～40°25′N)，北起燕山，南至黃河，東臨渤海，西以太行山為界，南北方向長約630 km，東西方向長約640 km，總面積約為13.92 km2。該地區是我國重要的農業基地，主要種植作物有小麥、玉米、棉花等。多年平均降水量為539.0 mm，降水季節分布不均勻，多集中在7—9 月份，且年際變化較大;年平均氣溫約為13.0 ℃;年平均日照時數約為2 430 h;年平均風速可達2.27 m/s。選取平原外圍4個不同方位的氣象站點作為研究對象(如圖2所示)，分別是密云(116°52′E,40°23′N)、靈壽(114°23′E,38°18′N)、延津(114°11′E,35°9′N)和沾化(118°7′E,37°41′N)。這4個站點具有較好的代表性，包含有1970年1月至2019年11月期間的每日氣象資料，氣象資料來源于中國氣象局國家氣象信息中心(http://data.cma.cn)地面氣象資料，數據質量可靠，總體正確率接近100%。收集的日尺度常規氣象因素包括相對濕度RH、日照時數n、2 m高處平均風速u2(由10 m高處平均風速換算得到)、平均氣溫Tmean、最高氣溫Tmax和最低氣溫Tmin，少數缺失的數據采用線性插補或往年同期數據進行移植，使其完整。

圖2 研究站點分布圖

3.2 研究準備

3.2.1 對照模型設置和相關參數

分別從算法和模型2個角度設置對照模型，包括PO-LSSVM模型、LSTM模型和RF模型。其中，PO算法是一種新型算法，已被相關實驗證實其具有較強的競爭力[23]。LSSVM模型γ和g的上下限范圍設為[e-10, e10]，e為自然對數。其余相關模型算法參數設置見表1。

表1 各模型有關參數設置

3.2.2 日尺度ET0模擬

選擇2015年1月份至2019年11月份(共1 795份樣本)的4個氣象站的日度氣象數據，包括相對濕度RH、日照時數n、2 m高處平均風速u2、平均氣溫Tmean、最高氣溫Tmax和最低氣溫Tmin，以及旬序數共7種影響因子作為輸入樣本。將Penman-Monteith

公式計算ET0值作為標準值，并將樣本數據的前2/3(1 197 份)作為訓練樣本，剩余1/3(598 份)作為檢驗樣本。利用LASSO回歸法確定不同要素組，可設置為:①全氣象數據和旬序數;②RH、n、u2、Tmean和旬序數;③n、u2、Tmean和旬序數。

3.2.3 月尺度ET0模擬

以密云站為例，利用自相關系數和偏自相關系數確定前期ET0序列對當前時間ET0的影響程度。圖3給出了原始月度ET0序列和8個基于db4小波函數的小波包分解子序列，8個子序列的直接疊加等于原始序列。

圖3 密云站ET0序列及其基于小波包分解的結果

每個子序列將利用提出的CSA-SIM-LSSVM模型獨立運行10次取平均值，8個子序列分別預報完畢后直接疊加得到預報結果。本次模型采用1970.01—2019.12序列末尾的100個月ET0作為檢驗樣本，訓練樣本由總樣本數和延遲預報月數確定。

3.2.4 評價標準

(14)

(15)

(16)

3.3 日尺度ET0估算結果對比及分析

以FAO Penman-Monteith公式為標準，按照3種不同影響因子輸入組合，為減小隨機性對本次研究的影響，對本節涉及的4種模型進行10次獨立運行輸出平均檢驗估算結果，記錄帶有模型優化過程的CSA-SIM-LSSVM模型、PO-LSSVM模型及LSTM模型的平均運行時間。3個模型的平均運行時間分別為9.16、25.85、14.99 s。4個站點上的相關平均評價結果見表2。從表2中可以得到以下結論:

1)從整體表現看，各個模型檢驗期的表現相比訓練期絕大多數有略微的下降，總的來看，沒有出現“過擬合”現象。

2)從不同類型模型之間來看，CSA-SIM-LSSVM模型表現最好，LSTM模型次之，RF模型表現最差。從算法角度進行對比，PO算法雖然更加新穎，但是本次研究證明其在優化LSSVM模型參數問題上相對于組合優化算法CSA-SIM較弱，不僅直觀上表現較差，而且運行時間更長，效率較低。

3)從影響因子組合來看，CSA-SIM-LSSVM模型、PO-LSSVM模型及LSTM模型的表現從好到差為因子組合①>因子組合②>因子組合③。但值得注意的是，對于RF模型，表現次序則為因子組合②>因子組合①>因子組合③?；贑SA-SIM-LSSVM模型，LASSO回歸法剔除系數為0或影響程度極低項后誤差增加不明顯，相較而言可縮短模型訓練時間。對于LSSVM模型和LSTM模型整體而言，加入“旬序數”這一具有農業指導意義的因子項是有效合理的。

4)從其他文獻提出模型的日尺度ET0模擬結果來看：張育斌等提出的CSA-LSSVM模型的誤差為RMSE=0.39 mm/d[11];陳晟提出的基于多氣象特征融合的RF模型的誤差為RMSE=0.17 mm/d[17];邢立文等利用LSTM對華北地區ET0模擬的日度平均誤差為RMSE=2.753 mm/d[24]。相比較而言，本文提出的CSA-SIM-LSSVM模型在前兩種要素組合下的平均誤差RMSE僅為0.061 mm/d,體現了該模型的優越性。

表2 各模型日度ET0平均擬合與檢驗結果評價

3.4 月尺度ET0估算結果及分析

圖4給出了基于小波包分解的CSA-SIM-LSSVM(WPD-CSA-SIM-LSSVM)模型在4個站點上測試樣本的預報結果曲線圖，用圓圈表示Penman-Monteith公式計算值，模型的預測結果用星號表示。整體上4個站點的預報結果較為可靠，模型很好地捕捉了月ET0序列的局部變化趨勢，在4個站點的ET0突變處、“雙峰”處均能提供較準確的預報值。

圖4 各站點測試樣本預報結果

表3記錄了WPD-CSA-SIM-LSSVM模型在4個站點上檢驗與測試的平均評價結果。與其他文獻提出模型的月尺度ET0模擬結果對比來看：侯志強等對于河套地區臨河氣象站的測試誤差為MRE=11.69%[9];鞠彬等對于哈巴河氣象站的模擬最小誤差為MRE=13.52%[10];邢立文等利用LSTM對華北地區ET0模擬的月度平均誤差為RMSE=1.460 mm/d[24]。本文提出的WPD-CSA-SIM-LSSVM模型月度模擬平均誤差MRE=1.13%，RMSE=0.028 mm/d，進一步體現了提出模型的優越性。

表3 各模型月度ET0平均擬合與檢驗結果評價

4 結論

利用耦合退火算法的全局尋優能力，配合單純形法良好的局部挖掘能力，并利用兩者耦合算法進行LSSVM模型的參數尋優以進行華北平原的日尺度和月尺度ET0的預測估算。評價結果顯示，提出的CSA-SIM-LSSVM模型在日尺度ET0估算上的表現明顯好于其他對照模型，展示了其提出的合理性。配合小波包分解技術，后續月尺度ET0的預報估算進一步展示了其可行性。本文提出的CSA-SIM-LSSVM模型具有較好的模擬效果，可為實際生產工作提供參考。