劉博弈，王海渝，龔 嚴(yán),2，劉文豪，衛(wèi) 琦,3，徐俊增,3

(1.河海大學(xué)農(nóng)業(yè)科學(xué)與工程學(xué)院，南京 210098; 2.河海大學(xué)水利水電學(xué)院，南京 210098; 3.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，南京 210098)

0 引 言

中國農(nóng)業(yè)用水占總用水比例大，發(fā)展節(jié)水灌溉對(duì)農(nóng)業(yè)的可持續(xù)發(fā)展至關(guān)重要[1]。節(jié)水灌溉的關(guān)鍵是進(jìn)行實(shí)時(shí)的灌溉預(yù)報(bào)，而逐日作物騰發(fā)量的預(yù)測(cè)是實(shí)時(shí)灌溉預(yù)報(bào)的基礎(chǔ)[2]。參考作物騰發(fā)量(Reference crop evapotranspiration，ET0)是計(jì)算作物騰發(fā)量的重要參數(shù)之一，ET0的準(zhǔn)確預(yù)報(bào)對(duì)于灌溉計(jì)劃以及區(qū)域水資源分配具有重要的指導(dǎo)意義。

鑒于上述研究不足，本文以天氣預(yù)報(bào)和符號(hào)回歸算法為基礎(chǔ)，擬提出一種新的ET0預(yù)測(cè)模型，并在干旱區(qū)和濕潤區(qū)分別選擇典型站點(diǎn)對(duì)其預(yù)測(cè)精度進(jìn)行評(píng)價(jià)，其研究結(jié)果對(duì)于不同氣候區(qū)實(shí)現(xiàn)ET0準(zhǔn)確預(yù)測(cè)以及制定有效合理的灌溉計(jì)劃具有重要意義。

1 材料與方法

1.1 數(shù)據(jù)來源

本文選取位于濕潤區(qū)和干旱區(qū)的6個(gè)典型站點(diǎn)(各3個(gè))作為研究對(duì)象，從中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.gov.cn)分別收集了6個(gè)站點(diǎn)1989-2016年逐日歷史實(shí)測(cè)氣象資料，包括歷史逐日的最高氣溫、最低氣溫、平均氣溫、日照時(shí)數(shù)、平均風(fēng)速以及平均相對(duì)濕度。并在天氣網(wǎng)(http:∥www.tianqi.com)收集了6個(gè)站點(diǎn)在2013-2016年逐日歷史天氣預(yù)報(bào)，包括日最高氣溫、日最低氣溫、天氣類型和風(fēng)力等級(jí)。6個(gè)站點(diǎn)的基本信息如表1所示。

表1 典型站點(diǎn)基本信息Tab.1 Basical information of typical stations

1.2 ET0預(yù)報(bào)模型

1.2.1 HS方法及其地區(qū)率定

HS公式是一種基于溫度估算ET0的方法，由于其需要的氣象數(shù)據(jù)少，常用于ET0預(yù)報(bào)。其表達(dá)式[13]為：

ET0,HS=0.408×a×(Tmean+17.8)×TRc×Ra

(1)

式中：ET0,HS為采用HS公式計(jì)算得到的參考作物騰發(fā)量，mm/d；TR和Tmean分別為逐日溫差以及平均溫度，℃，可由日最高氣溫Tmax和最低氣溫Tmin計(jì)算得到，Tmean=(Tmax+Tmin)/2，TR=Tmax-Tmin；Ra為天頂輻射，MJ/(m2·d)；參數(shù)a和參數(shù)c需要根據(jù)當(dāng)?shù)貧庀筚Y料率定。

由于FAO56 PM公式[14]計(jì)算精度高，被廣泛用于評(píng)價(jià)其他ET0預(yù)報(bào)和估算模型。本研究以1989～2009年歷史實(shí)測(cè)氣象資料為率定樣本，以PM公式計(jì)算的ET0為對(duì)照，對(duì)HS公式的參數(shù)a和參數(shù)c進(jìn)行率定，并以2010-2016年氣象數(shù)據(jù)為驗(yàn)證樣本。將2013-2016年氣象預(yù)報(bào)數(shù)據(jù)輸入，對(duì)率定后的HS公式進(jìn)行ET0預(yù)報(bào)精度分析。

1.2.2 符號(hào)回歸算法

圖1所示為單相PWM整流的拓?fù)鋱D,Un和In分別為電路交流側(cè)的電壓和電流,Rn為線路電感,Ln為交流側(cè)電流,Udc為直流側(cè)電壓,C2為直流側(cè)儲(chǔ)能電容,單相PWM整流電路主要通過將交流能量通過IGBT的開關(guān)狀態(tài)轉(zhuǎn)換成直流能量,通過調(diào)整不同開關(guān)狀態(tài)的持續(xù)時(shí)間,可以控制直流側(cè)電壓的大小,同時(shí)保證交流側(cè)的輸入電壓和輸入電流為同相位,保證單位功率因數(shù)的能量傳輸。

符號(hào)回歸算法(SR)是在Eureqa軟件應(yīng)用上逐步發(fā)展完善[15, 16]。運(yùn)用符號(hào)回歸算法建立模型，需要確定輸入的自變量以及因變量，函數(shù)組成以及適應(yīng)度函數(shù)[12]。本文選用的函數(shù)組成包括C(常數(shù))、x(輸入變量)、ex、lg(x)、sqrt、冪、+、-、×、÷，選擇擬合優(yōu)度R2作為衡量標(biāo)準(zhǔn)。擬合優(yōu)度越大，模型估算精度越高。程序自動(dòng)運(yùn)行后將對(duì)產(chǎn)生的模型進(jìn)行反復(fù)修正評(píng)估，直到產(chǎn)生的模型穩(wěn)定且滿足衡量標(biāo)準(zhǔn)，即運(yùn)行結(jié)束。本文從產(chǎn)生的模型中挑選出擬合優(yōu)度R2最大的模型作為SR模型。

以最高氣溫Tmax，最低氣溫Tmin，相對(duì)日照時(shí)數(shù)(Relative sunshine duration，n/N)，天頂輻射Ra為自變量，以1989-2009年歷史實(shí)測(cè)氣象資料為率定樣本，以PM公式計(jì)算的ET0為對(duì)照，經(jīng)2010-2016年驗(yàn)證期驗(yàn)證后建立ET0符號(hào)回歸估算模型。其模型如下：

ET0,SR=f(Tmax,Tmin,n/N,Ra)

(2)

n/N是估算凈輻射的重要參數(shù)，一些學(xué)者將天氣預(yù)報(bào)中天氣類型解析為相對(duì)日照時(shí)數(shù)用于ET0預(yù)測(cè)[7, 17]，但天氣類型預(yù)報(bào)精度需要提高。也有相關(guān)研究利用溫度建立相對(duì)日照時(shí)數(shù)定量表示模型[18]。本文以最高氣溫Tmax，溫差TR為自變量，以1989-2009年逐月歷史實(shí)測(cè)資料為率定樣本，以實(shí)測(cè)資料計(jì)算的n/N為對(duì)照，經(jīng)2010-2016驗(yàn)證期驗(yàn)證后建立逐月n/N符號(hào)回歸估算模型，模型表達(dá)式如下：

n/N=f(Tmax,TR)

(3)

式中：Tmax和TR分別為日最高氣溫和溫差，℃；n/N為逐日相對(duì)日照時(shí)數(shù)，由日照時(shí)數(shù)n和白晝最大可能時(shí)長N來確定。具體公式見文獻(xiàn)[14]。

將2013-2016年天氣預(yù)報(bào)數(shù)據(jù)先代入到模型(3)得到n/N預(yù)報(bào)值，再將同時(shí)期溫度預(yù)報(bào)數(shù)據(jù)與n/N預(yù)報(bào)值代入到模型(2)進(jìn)行ET0預(yù)報(bào)。

1.3 誤差分析

采用均方根誤差(RMSE)、平均絕對(duì)誤差(MAE)和加權(quán)決定系數(shù)(WR2)等3個(gè)統(tǒng)計(jì)指標(biāo)來評(píng)價(jià)ET0預(yù)報(bào)精度，并采用準(zhǔn)確率分別評(píng)價(jià)氣溫預(yù)報(bào)和ET0預(yù)報(bào)，具體公式參見文獻(xiàn)[8, 12]。

2 結(jié)果與討論

2.1 氣溫預(yù)報(bào)精度評(píng)價(jià)

典型站點(diǎn)的最低和最高氣溫預(yù)報(bào)準(zhǔn)確率如圖1所示?？梢钥闯?，位于干旱區(qū)的典型站點(diǎn)的最高、最低氣溫預(yù)報(bào)的平均準(zhǔn)確率(72.39%和79.41%)均略低于位于濕潤區(qū)的典型站點(diǎn)的值(82.48%和84.98%)，且這一差異在最高氣溫的預(yù)報(bào)精度上表現(xiàn)更為明顯?？梢?，氣溫預(yù)報(bào)在濕潤區(qū)各站點(diǎn)的精度比在干旱區(qū)各站點(diǎn)更為穩(wěn)定，這可能與采用不同的氣溫預(yù)報(bào)模型相關(guān)[19]。整體上，不同氣候區(qū)的典型站點(diǎn)最低氣溫預(yù)報(bào)的準(zhǔn)確率均略高于最高氣溫預(yù)報(bào)準(zhǔn)確率，與李國翠等[20]的研究結(jié)果較為一致。

圖1 各站點(diǎn)氣溫預(yù)報(bào)準(zhǔn)確率的比較Fig.1 Comparison of accuracy of the forecasted temperature in different stations

2.2 HS公式率定與驗(yàn)證

表2給出了HS公式參數(shù)率定結(jié)果，可以看出，位于相同氣候區(qū)的站點(diǎn)參數(shù)取值大致落在相同范圍，且參數(shù)值范圍與胡慶芳等[21]研究結(jié)果較為一致，差別可能是由于所采用時(shí)間步長不同導(dǎo)致。

表2 HS公式參數(shù)修正結(jié)果Tab.2 Results of parameter correction of the HS model

表3 HS公式率定期和驗(yàn)證期統(tǒng)計(jì)指標(biāo)Tab.3 Statistical indices for HS models in calibration period and validation period

以PM公式計(jì)算結(jié)果為對(duì)照，進(jìn)一步分析率定后的HS公式在各站點(diǎn)率定期和驗(yàn)證期的誤差情況(表3)，可以看出，率定后的HS公式在各站點(diǎn)的率定期和驗(yàn)證期其MAE和RMSE均維持在較低水平，且HS公式在位于干旱區(qū)的站點(diǎn)精度略高于位于濕潤區(qū)的站點(diǎn)，但在總體上率定后的HS公式在各站點(diǎn)均保持較高的精度，可以應(yīng)用到預(yù)報(bào)中。

2.3 ET0符號(hào)回歸估算模型

以PM公式和1989～2009年歷史氣象資料計(jì)算的ET0作為基準(zhǔn)，通過符號(hào)回歸算法建立ET0與各氣象因素之間的定量關(guān)系，以決定系數(shù)(R2)最大作為篩選依據(jù)，篩選出各站點(diǎn)的ET0符號(hào)回歸估算模型，分別如下式(4)～(9)所示：

0.026Tmax(n/R)Ra-0.65

(4)

ET0,SR常州=0.171 3+0.045 69Ra+0.022 98Tmin+

(5)

0.002 994Tmax(n/R)Ra-0.310 5

(6)

ET0,SR西寧=0.11(n/N)Ra+0.002 4TmaxRa-

0.020 8 8Tmax-2.2 (n/R)+0.82

(7)

ET0,SR銀川=0.032 68Ra+0.002 093TmaxRa+

(8)

ET0,SR烏魯木齊=0.05Ra+9×10-4TmaxRa+

(9)

以2010-2016年實(shí)測(cè)氣象數(shù)據(jù)為驗(yàn)證樣本，以PM公式計(jì)算結(jié)果為對(duì)照，對(duì)ET0符號(hào)回歸估算模型的估算精度(圖2)進(jìn)行分析，可以看出，與率定后的HS公式相比，SR模型在各站點(diǎn)的ET0計(jì)算精度顯著提升，回歸直線方程斜率均接近1，加權(quán)決定系數(shù)WR2均超過了0.91，干旱區(qū)和濕潤區(qū)各站點(diǎn)采用SR模型的MAE值較率定后的HS公式分別平均降低了49.56%和28.40%，RMSE分別減少了48.90%和29%?？梢姡琒R模型在引入相對(duì)日照時(shí)數(shù)輸入變量后，ET0估算精度得到明顯提高。

2.4 基于天氣預(yù)報(bào)數(shù)據(jù)的ET0預(yù)測(cè)結(jié)果評(píng)價(jià)

以2013-2016年歷史逐日天氣預(yù)報(bào)數(shù)據(jù)為依據(jù)，以PM公式計(jì)算值(ET0,PM)為對(duì)照，分別對(duì)率定后的HS公式和SR模型的ET0預(yù)報(bào)值(ET0,HS和ET0,SR)進(jìn)行評(píng)價(jià)(如圖3所示)，可以看出，整體上，兩種方法在各站點(diǎn)的預(yù)報(bào)值均與PM公式計(jì)算值變化趨勢(shì)較為一致，與PM公式相比，引入相對(duì)日照時(shí)數(shù)的SR模型和率定后的HS模型對(duì)ET0的預(yù)報(bào)精度均有明顯提升，SR模型的ET0預(yù)報(bào)精度略高于率定后的HS模型，且其提升幅度在濕潤區(qū)各站點(diǎn)表現(xiàn)得更為明顯。

進(jìn)一步分析SR模型和HS模型的ET0預(yù)報(bào)精度(表4)，可以發(fā)現(xiàn)，與率定后的HS公式的計(jì)算結(jié)果相比，位于濕潤區(qū)的各站點(diǎn)采用SR模型的MAE、RMSE值分別平均減少了18.98%和20.97%，WR2和準(zhǔn)確率分別平均增加了24.65%和5.57%。而位于干旱區(qū)的各站點(diǎn)SR預(yù)報(bào)模型MAE、RMSE值分別減少了9.79%和7.53%，WR2和準(zhǔn)確率分別平均增加了0.49%和0.80%。

HS公式雖然預(yù)報(bào)精度整體較高，但仍存在一些缺陷。由于不考慮風(fēng)速以及相對(duì)濕度的影響，在極端風(fēng)速情況下，預(yù)報(bào)精度較低[9]。例如位于干旱區(qū)的烏魯木齊2015年7月22日ET0達(dá)到了11.6 mm/d，最高溫與最低溫度預(yù)報(bào)值只相差1.4 ℃和0.8 ℃，而實(shí)測(cè)風(fēng)速為5.9 m/s，使得HS預(yù)報(bào)值僅為7.7 mm/d。此外在位于不同氣候區(qū)的站點(diǎn)預(yù)報(bào)精度差異較大，位于濕潤區(qū)各站點(diǎn)HS預(yù)報(bào)的MAE值平均比干旱區(qū)增大了30.3%，HS公式在位于濕潤區(qū)的各站點(diǎn)預(yù)報(bào)效果較差，這也與Xu等[2]研究結(jié)論一致。

圖2 驗(yàn)證期符號(hào)回歸模型ET0估算值與PM公式ET0計(jì)算值的關(guān)系Fig.2 Relation in ET0 between those caculated by SR and PM models during validation period

圖3 2013-2016年ET0預(yù)報(bào)值和PM計(jì)算值過程變化圖Fig.3 Variation patterns for forecasted ET0 and ET0 calculated from PM model (2013-2016)

站點(diǎn)方法MAE/(mm·d-1)RMSE/(mm·d-1)WR2準(zhǔn)確率/%武漢SR0.520.720.7693.70HS0.690.910.5489.18常州SR0.550.700.7696.44HS0.660.890.5890.27杭州SR0.540.690.7596.29HS0.640.870.5990.28西寧SR0.320.480.8798.70HS0.370.520.8597.88銀川SR0.480.670.7995.28HS0.520.730.8493.98烏魯木齊SR0.450.700.9193.97HS0.490.750.8693.70

符號(hào)回歸算法通過引入相對(duì)日照時(shí)數(shù)，使得ET0預(yù)報(bào)精度較HS公式提高，但在不同氣候區(qū)提升幅度不同。濕潤區(qū)夏季雨量增多導(dǎo)致濕度增大，使得HS公式在夏季計(jì)算誤差較大，且濕潤區(qū)夏季天氣類型變化較大[9]，引入相對(duì)日照時(shí)候后即增加了天氣類型的修正，因而使得在濕潤區(qū)提升幅度較高。例如武漢站點(diǎn)HS模型在夏季AAE值為2.04 mm/d，而SR模型AAE值僅為1.32 mm/d。而干旱區(qū)風(fēng)速ET0對(duì)風(fēng)速較為敏感，SR模型沒有考慮風(fēng)速的影響，且HS公式在干旱區(qū)預(yù)報(bào)精度較高，因此在干旱區(qū)SR模型提升幅度較小。

3 結(jié) 論

本研究提出了一種基于符號(hào)回歸算法(SR)和天氣預(yù)報(bào)的ET0預(yù)報(bào)方法，并應(yīng)用于干旱區(qū)和濕潤區(qū)的典型站點(diǎn)，并與率定HS公式的ET0預(yù)報(bào)結(jié)果進(jìn)行了比較。得出的主要結(jié)論如下：

(1)與率定后的HS公式相比，引入相對(duì)日照時(shí)數(shù)的SR模型的ET0預(yù)報(bào)精度明顯提高，且在位于濕潤區(qū)典型站點(diǎn)其預(yù)報(bào)精度的提高幅度大于位于濕潤區(qū)的典型站點(diǎn)。

(2)對(duì)于干旱區(qū)的站點(diǎn)，由于率定后的HS公式預(yù)報(bào)精度較高，SR模型提升幅度較小，因此建議使用HS公式進(jìn)行ET0預(yù)報(bào)；對(duì)于天氣類型(尤其是夏季)變化較大的濕潤區(qū)站點(diǎn)，由于SR模型提升幅度較高，因此建議使用SR模型進(jìn)行預(yù)報(bào)。

□