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

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

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

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

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

1 研究方法

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

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

(1)

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

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

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

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

圖1 參照蒸散量的LSSVM模型

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

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

(2)

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

考慮復雜的函數(shù)擬合誤差,給出的回歸問題可以根據(jù)結(jié)構(gòu)最小化原則轉(zhuǎn)換，相當于求解下面的最小值問題：

(3)

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

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

(4)

根據(jù)KKT條件，最優(yōu)條件可以通過分別對w,b,ξi和αi求偏導數(shù)得到。

(5)

將公式(4)簡化記成：

(6)

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

(7)

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

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

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

超參數(shù)生成過程如下式：

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

(8)

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

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

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

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

(9)

式中：Tk是k時刻生成溫度；k是T的調(diào)整次數(shù)。

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

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

(10)

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

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

(11)

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

(12)

(13)

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

2 研究區(qū)域和評價方法

2.1 地區(qū)案例研究

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

2.2 評價方法

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

表1 作物參照蒸散量ET0經(jīng)典計算模型

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

(14)

(15)

(16)

(17)

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

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

3 應用結(jié)果與分析

3.1 氣象數(shù)據(jù)缺失下CSA-LS-SVM模型的模擬精度

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

表2 原始氣象數(shù)據(jù)與ET0的相關(guān)性分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表5 驗證期模型總ET0性能統(tǒng)計表

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

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

圖5 CSA-LS-SVM模型在決策系統(tǒng)運行情況

4 結(jié) 語

基于LSSVM算法研究基礎(chǔ)上，提出耦

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

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

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

(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].農(nóng)業(yè)工程學報，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].中國農(nóng)村水利水電，2014,(12):64-68.

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

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