耦合模拟退火优化最小二乘支持向量机的日参照蒸散量模拟计算

张育斌，魏正英，张 磊，张 帅，蔚磊磊

(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.