耦合模擬退火優化最小二乘支持向量機的日參照蒸散量模擬計算

張育斌，魏正英，張 磊，張 帥，蔚磊磊

(1.西安交通大學機械制造系統工程國家重點實驗室，西安 7100454；2.寧波大紅鷹學院，寧波 315175)

參照作物蒸散量(ET0)是計算作物需水量的關鍵，精準的計算及預測對作物的按需生長及區域水資源管理具有重要意義。而參照作物需水量ET0的計算方法有很多，其中ASCE使用分布在世界各大洲不同氣候條件下的實驗儀器所測數據作為參照，發現彭曼-蒙蒂斯公式(PM公式)計算的ET0與實測值都很接近[1]。ET0精準模擬計算對灌溉工程師和農戶掌握作物需水情況，做出灌溉分析和決策具有重要意義。

ET0的模擬計算過程看作是氣象參數間的復雜非線性回歸過程，伴隨智能算法的發展，該領域研究學者提出了基于智能算法研究的ET0模擬計算模型。Ozgur Kisi等[2]研究基于最小二乘支持向量機的日參照蒸散量計算，模型在美國加州應用中模擬精度高于ANN的模型，但其模擬精度隨著輸入氣象因子項的減少而相對降低。馮禹等[3]建立基于極限學習機(ELM)的ET0模擬模型，模型在川中丘陵地區ET0模擬精度較高。Shiri等[4]建立了基于基因表達式編程(GEP)的ET0模擬模型，將其應用于西班牙巴斯克地區ET0模擬，發現 GEP 模擬性能高于ANFIS、Priestley-Taylor和H-Si模型。Kisi等[5,6]使用進化的前饋神經網絡、小波分析與經典模型進行ET0模擬對比研究。Matri等[7]、Kumar等[8]分別建立了基于BP神經網絡的ET0預測模型，但 BP 神經網絡需要較多的學習樣本，很難求得全局最優解。

本文在研究最小二乘支持向量機在陜西地區應用基礎上，考慮LS-SVM模型求解過程在得到簡化的同時，缺失了支持向量機擁有的魯棒性以及稀疏性，為了較好的反映樣本分布特征與 LS-SVM 優化模型的性能，采用耦合模擬退火( CSA) 算法優化LS-SVM超參數，避免超參數陷入局部極小值。并通過實驗，分析驗證改進算法在日參照作物蒸散量模擬計算中具有更快學習速度和更好的泛化能力。

1 研究方法

1.1 作物參照蒸散量(ET0)計算

作物需水量(ETc)是植物生長過程中一個非常重要的數據，反映植物生長發育時所需要的水量，是灌溉過程中重要的參數，影響區域水資源的優化配置，其等于作物系數與參考作物蒸散量(ET0)之積，與植物的參考作物蒸散量息息相關。ET0的計算方法有很多，常見的有溫度法、輻射法(模型法)、經驗公式法和綜合法等。當前，國內外應用最多的是FAO-56 Penman-Monteith方法(PM56法)，該法基于能量平衡和空氣動力學原理，考慮了大量的氣象因素，計算結果較為準確[9]，常作為標準值來研究。

(1)

式中：ET0為參照蒸散量，mm/d；Δ為溫度-飽和水汽壓關系曲線在T處的斜率，kPa/℃；T為2 m處日平均溫度，℃；Rn為作物表面凈輻射，MJ/(m2·d)；G為土壤熱通量，MJ/(m2·d)；γ為濕度表常數，kPa/℃；es為飽和水汽壓，kPa；ea為實際水汽壓，kPa；u2為距離地面2 m處的平均風速，m/s。

這些參數不是直接得到，而是通過氣象原始的8項數據(①日最低氣溫；②日最高氣溫；③地理緯度弧度；④海拔高度；⑤日平均溫度；⑥平均相對濕度；⑦實際日照時長；⑧風標高度及風速)計算得來，且有些參數難以測得，如果有缺失數據就不能計算作物需水量，因而，采用部分氣象因子和農情信息作為模型輸入量來預測ET0的值，是現實灌溉指導意義。

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

LSSVM本質上是一種變形算法。變形算法主要是通過增加函數項、變量或系數等方法使公式變形，從而產生出各種具有某一方面優勢或者一定應用范圍的算法[10]。圖1是作物參照蒸散量LSSVM回歸算法模型。

圖1 參照蒸散量的LSSVM模型

考慮輸入量xi(日最低氣溫、日最高氣溫、地理緯度弧度、海拔高度、日平均溫度、平均相對濕度、實際日照時長、風標高度及風速)和輸出量yi(蒸散量)時間序列，根據最小二乘支持向量機方法，非線性函數可以表達為：

f(x)=wφ(x)+b

(2)

式中：f(x)指的是氣象變量與ET0的關系；w是m維權向量；φ(x)是映射x到m維特征向量的映射函數；b是偏移選項；w和b是待求的回歸參數。

考慮復雜的函數擬合誤差,給出的回歸問題可以根據結構最小化原則轉換，相當于求解下面的最小值問題：

(3)

式中：γ為正則化參數；ξi是xi的松弛變量。

為了求解公式(3)中的最優問題，目標函數可以通過改變約束問題轉化為一個無拘無束的問題,引入Lagrange函數為：

(4)

根據KKT條件，最優條件可以通過分別對w,b,ξi和αi求偏導數得到。

(5)

將公式(4)簡化記成：

(6)

這樣作物參考蒸散量LSSVM的模擬計算問題就變成：

(7)

由式(7) 可以看出, 確定回歸參數的關鍵在于計算核相關矩陣的逆計算。

1.3 耦合模擬退火算法(CSA)

模擬退火在常規迭代算法的基礎上引入 Metropolis 準則，以一定接受概率接受擬合精度偏低的參數，避免算法陷入局部最優，并且伴隨著溫度參數的迭代下降，SA算法逐漸趨向全局最優。但是SA算法收斂速度變慢，不利于在線優化；另一方面，SA 算法魯棒性受到接受溫度的初始賦值和接受溫度規則選取的影響。本文借鑒文獻[11]耦合模擬退火算法(CSA)用于優化 LS-SVM的超參數的選擇。CSA是通過反映接受概率函數的能量和狀態等信息的耦合項來實現各個模擬退火過程的相互耦合與信息共享，接受函數與耦合項的逐漸匹配引導 CSA算法獲得全局最優解。耦合模擬退火算法主要包括超參數(正則化參數γ、核寬度σ)的生成和接受的兩個主要過程：

超參數生成過程如下式：

yi=xi+εixi∈θ,yi∈Ω

(8)

i=1,2,3,…,n;θ∈Ω

式中：Ω是超參數域；θ是當前超參數的集合；n是CSA數量；εi是隨機變量。

通常εi選用正切函數作為分布，這樣對應的式(8)中生成超參數的表達式為：

yi=xi+tan[π(τ-0.5)]Tk

(9)

式中：Tk是k時刻生成溫度；k是T的調整次數。

同時可知，yi的分布主要受Tk的影響，Tk越大，則yi分布范圍廣，反之變窄，而最終值也會逐漸收斂。在本文，采用生成溫度更新規則U(Tk,k)：

U(Tk,k)←Tk=T0/k

(10)

在CSA算法的超參數接受過程中，生成超參數被接受的概率由接受溫度Tac值的分布控制，SA算法接受溫度更新規則采用V(Tack,k)：

V(Tack,k)←Tack=Tac0×e-k-1

(11)

CSA 算法將多個獨立 SA 算法的接受概率耦合。CSA 超參數的接受概率函數為Aθ(β,xi→yi)，耦合項β是所有超參數xi對應擬合精度的函數，分別選為：

(12)

(13)

CSA算法生成溫度更新的兩個階段，首先設定超參數生成迭代次數m，超參數接受過程中，依據留一交叉驗證求得獨立的各個超參數擬合精度，然后，結合接受概率耦合函數式(12)，CSA接受擬合精度高的超參數的同時，并行接受擬合效果差的超參數，從而促進獨立SA算法優化過程中，實現信息的共享。

2 研究區域和評價方法

2.1 地區案例研究

本文氣象資料來自中國氣象數據共享服務系統，數據經過嚴格質量控制，質量良好。選取陜西地區 3個氣象站點：榆林(緯度38°16′, 經度109°47′)、安康(緯度32°43′, 經度109°02′)和西安(緯度34°18′, 經度108°56′)的1971-2014年逐日中8項氣象資料進行訓練與測試，包括日平均氣溫(T)、最高氣溫(Tmax)、最低氣溫(Tmin)、日照時數(n)、距地面10 m 高處的風速(Uh)、2 m 高度相對濕度(RHm)，其中以 1971-2000年逐日氣象資料進行訓練，2001-2009年逐日氣象資料進行測試，2010-2014年逐日氣象資料進行驗證。

2.2 評價方法

本文主要集中在LSSVM相關模型與Hargreaves, Priestley-Taylor(P-T)這兩種典型模型對比研究。其中ET0計算是采用世界糧農組織推薦FAO-56 PM公式作為計算ET0的標準公式。

表1 作物參照蒸散量ET0經典計算模型

根據以上內容，本文首次使用原始氣象的8項數據作為輸入量，采用FAO-56 PM方法的ET0輸出量作為模型的校對值。而均方根誤差(RMSE),平均絕對誤差(MAE)和自相關系數(R2)的統計數據作為模型評價的主要因素。RMSE,MAE和R2定義如下：

(14)

(15)

(16)

(17)

式中：N表示數據的數量；x和y分別是預測值和FAO-56的ET0值。

RMSE和MAE越小，表明模型偏差越小；EF越接近于 1，表明模型效率越高；R2越接近1，表明吻合度越高。

3 應用結果與分析

3.1 氣象數據缺失下CSA-LS-SVM模型的模擬精度

選擇8項氣象參數，以不同組合的氣象因子做為輸入量，采用FAO 56 Penman-Montieth公式的計算結果作為模型的預期輸出值，建立相應的 CSA-LS-SVM 模型，考慮組合的有效性，對其8個原始氣象數據與ET0的相關性進行研究，結果見表2。

表2 原始氣象數據與ET0的相關性分析

從表2可知，與參照蒸散量ET0相關的氣象數據的重要性，最大的是日最高氣溫，其次日平均氣溫，日最低氣溫，日照時長，風速，海拔和緯度，最后與濕度關系系數呈負相關，且相關關系極顯著。

參考原始氣象數據與ET0的相關性，選擇安康地區的，針對性的氣象因子進行組合作為為輸入量，探討氣象資料缺失對CSA-LS-SVM模型模擬精度的影響，結果見表3。

表3 不同氣象參數組合下的CSA-LS-SVM的ET0模擬精度

當采用全部8個氣象因子進行模擬計算時，與其6個氣象因子進行模擬計算的效果一樣，符合表2所述海拔和緯度影響系數非常小，幾乎可以忽然不計。溫度對預測結果的影響非常大，尤其是最高溫度。只要有溫度，即便輸入量減少到兩個，也可以預測出非常好的結果，其R2和EF均在0.90以上，具有較好的相關性和模型有效性。只要帶有最高溫度的2個氣象因素組合進行模擬計算時，模擬精度都較高，R2和EF均在0.95以上；其他2個氣象因素組合進行模擬計算時，模擬精度明顯下降，其中缺少與溫度相關條件時，模擬結果較差，相對誤差較大。

3.2 CSA-LS-SVM模型與其他模型模擬精度對比

以FAO-56 PM 模型為標準, 將CSA-LS-SVM模型與最小二乘支持向量機(LSSVM)模型以及經典ET0計算方法的 Priestley-Taylor和Hargreaves模型模擬精度進行比較，參數選擇相同下模型誤差，見表4。

表4 CSA-LS-SVM模型與其他模擬模型比較

由表4可知，CSA-LS-SVM模型取得模擬計算效果，較其他模型的誤差小，總體都在0.45以內。

為進一步對 CSA-LS-SVM模型的可移植性及誤差成因進行分析，不同氣象站所有模型的蒸散量的模型計算結果見圖2～圖4中散點圖所示，我們可以清楚地看到6個氣象原始數據作為CSA-LS-SVM模型輸入所模擬計算結果是最接近 FAO-56 PM 計算得到ET0值，比其他幾個模型都要高。該模型的線性相關性(假設滿足公式y=ax+b)和自相關系數R2比其他模型都好。這里6參數模型系數a和b分別接近1和0，自相關系數R2比其他模型的值要高。

由圖2知，對于榆林地區CSA-LSSVM模型計算值比其他模型更接近FAO-56 PM的標準值，LSSVM 和 經典Hargreaves這兩模型幾乎相等，自相關系數分別為0.931 6和0.936 7。但P-T經典公式所得計算自相關系數最差。

由圖3可知，安康地區的CSA-LS-SVM模型模擬計算效果最好R2=0.988 9，其次經典的Hargreaves 模型模擬計算效果較好R2=0.983 3，而LSSVM模型模擬計算效果也很好，其R2=0.958 9，比Priestley-Taylor模型好。這幾種方式都可以適用，但是Priestley-Taylor模型最好做修正。

在圖4中，由于西安地區氣象數據缺失，只有統計到2010年的氣象數據，且訓練2006年缺失，模型測試是采用08-10數據作為測試，結果可以看出，測試中CSA-LS-SVM模型模擬計算效果最好R2=0.950 3，其次是LSSVM模型，R2=0.948 8；經典Priestley-Taylor和Hargreaves模型中模型模擬R2分別為0.880 4和0.833 5，在該地區可以看出缺失數據會對模型產生較大影響。

圖2 榆林模型模擬值與 FAO-56 PM的標準值比較

圖3 安康地區模型模擬值與 FAO-56 PM標準值比較

圖4 西安地區模型模擬值與FAO-56 PM標準值比較

3.3 CSA-LS-SVM模型的驗證分析

考慮ET0預測對于灌溉管理的重要性，進一步計算ET0的總值模擬計算精度，見表5。從表中相對誤差指標來看，智能算法CSA-LS-SVM和LSSVM模型在4個地區應用指標明顯優于經典模型，所模擬計算所得驗證期間的ET0總值是最接近FAO-56 PM模型的ET0總值，相對誤差都在10%以內。

對于榆林地區，CSA-LS-SVM和LSSVM模型是同一級精度，且CSA-LS-SVM預測相對誤差4.35%好于LSSVM的8.91%，是最好的計算精度，其次為LSSVM模型；經典模型中Hargreaves和Priestley-Taylor模型都相對誤差都較高，達到40% 以上，說明選擇的經典模型雖部分線性相關性較好，但誤差大，不適用。在安康地區，情況與榆林地區相似，CSA-LS-SVM預測模型計算精度很好，相對誤差都在1%以內；同時Hargreaves好于 Priestley-Taylor模型，但相對誤差達到20%以上。西安地區CSA-LS-SVM模型優于其他模型，其次LSSVM和Hargreaves模型都較好，相對誤差在5%以內，而Priestley-Taylor模型模擬計算精度最差。從總體來看，情況就西安還有點不同，原因估計還是跟數據缺失有關，造成部分經典公式相對誤差較小。

表5 驗證期模型總ET0性能統計表

3.4 CSA-LS-SVM模型在KBE決策模型中應用

將以上研究的CSA-LS-SVM模型嵌入于基于Delphi平臺開發基于知識工程灌溉決策系統中進行ET0決策分析，同時調用網頁上天氣預報，如西安從2016年4月26日未來幾天的天氣溫度和濕度分別為：15～25， 14～28， 15～29， 16～29， 17～31 ℃，56%。運行基于知識工程決策系統中決策ET0分別為2.066，2.125，2.511，2.418，2.651 mm，如圖5所示，這些預測結果，有利于指導灌溉。

圖5 CSA-LS-SVM模型在決策系統運行情況

4 結 語

基于LSSVM算法研究基礎上，提出耦

合模擬退火優化LSSVM超參數問題，CSA-LS-SVM算法在作物日參照蒸散量模擬計算，結果表明其能取得較高的模擬精度，均優于其他模型，有著較好的應用前景。

(1)CSA-LS-SVM模型有效解決了一般LSSVM訓練速度慢、參數選擇困難等缺陷問題。

(2)本文首次提出用8項氣象原始監測數據，建立CSA-LS-SVM模型，在陜西3個地區其模擬計算精度略高于LSSVM模型，明顯高于其他經典模型(Hargreaves和 Priestley-Taylor)模擬結果；在ET0的總計算值模型驗證中CSA-LS-SVM和LSSVM模型都較接近總FAO-56ET0標準值，經典模型只能部分地區適用，但應用時要對模型進行修正。

(3)在氣象資料缺乏情況下，只要獲取當日最高氣溫，再加其他一個氣象因子，CSA-LS-SVM模型模擬較FAO-56ET0標準模型的預測精度在92%以上，為作物需水量的智能灌溉決策提供參考值。

[1] Allen R G, Pereiral L S, Raes D, et al. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements[M]. Rome: FAO Irrigation and Drainage 56, 1998.

[2] 馮 禹,崔寧博.基于極限學習機的參考作物蒸散量預測模型[J].農業工程學報，2015,31(1):153-158.

[3] Ozgur Kisi.Least squares support vector machine for modeling daily reference evapotranspiration[J]. Irrig Sci, 2013.31:611-619.

[4] Shiri J, Kisi O, Landeras G, et al. Daily reference evapotranspiration modeling by using genetic programming approach in the Basque Country (Northern Spain)[J]. Journal of Hydrology, 2012,414/415:302-316.

[5] Kisi O. Evapotranspiration modeling using a wavelet regression model[J]. Irrig Sci,2011a,29(3):241-252.

[6] Kisi O.Modeling reference evapotranspiration using evolutionary neural networks[J]. ASCE J Irrig Drain Eng, 2011b,137(10):636-643.

[7] Marti P, Gonzalez-Altozano P, Gasque M. Reference evapotranspiration estimation without local climatic data[J]. Irrig Sci, 2011b,29:479-495.

[8] Kumar M, Raghuwanshi NS, Singh R.Artificial neural networks approach in evapotranspiration modeling: a review[J]. Irrig Sci, 2011,29:11-25.

[9] 張育斌，魏正英.極端天氣下作物參照蒸散量計算方法研究[J].中國農村水利水電，2014,(12):64-68.

[10] 郭振凱, 宋召青, 毛劍琴.一種改進的在線最小二乘支持向量機回歸算法[J].控制與決策，2009，1(24)：145-148.

[11] 衷路生，陳立勇等. 基于耦合模擬退火優化最小二乘支持向量機的車輪踏面磨耗量預測[J].計算機應用研究,2015,2(32):397-402.