多元自適應(yīng)回歸樣條算法模擬川中丘陵區(qū)參考作物蒸散量

陳宣全，崔寧博,2,3，李繼平，徐浩若，劉雙美，麻澤龍，樂進(jìn)華，王 軍

多元自適應(yīng)回歸樣條算法模擬川中丘陵區(qū)參考作物蒸散量

陳宣全1，崔寧博1,2,3※，李繼平1，徐浩若1，劉雙美4，麻澤龍4，樂進(jìn)華5，王 軍5

（1. 四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室水利水電學(xué)院，成都 610065；2. 南方丘區(qū)節(jié)水農(nóng)業(yè)研究四川省重點實驗室，成都 610066；3. 西北農(nóng)林科技大學(xué)旱區(qū)農(nóng)業(yè)水土工程教育部重點實驗室，楊凌 712100；4. 四川省水利科學(xué)研究院，成都 610072；5. 北京東方潤澤生態(tài)科技股份有限公司，北京 100083）

參考作物蒸散量（reference crop evapotranspiration, ET0）是作物精準(zhǔn)灌溉管理與農(nóng)業(yè)高效用水的核心參數(shù)。為提高川中丘陵區(qū)氣象資料缺省下的ET0預(yù)報精度，利用不同的氣象因子組合，建立15種基于多元自適應(yīng)回歸樣條算法（multivariate adaptive regression splines, MARS）的ET0預(yù)報模型。選取11個代表性氣象站點1961—2016年逐日氣象資料進(jìn)行分析，將其與其他ET0預(yù)報模型進(jìn)行對比，并利用可移植性分析評價MARS模型在川中丘陵區(qū)的適用性。結(jié)果表明：基于溫度和風(fēng)速項輸入的MARS5（輸入大氣頂層輻射、最高氣溫、最低氣溫、2 m處風(fēng)速）、MARS9（輸入最高氣溫、最低氣溫、2 m處風(fēng)速）和MARS13（輸入最高氣溫、2 m處風(fēng)速）模型，以及僅基于風(fēng)速項輸入的MARS15模型都具有良好的模擬精度；大氣頂層輻射和風(fēng)速是決定機(jī)器學(xué)習(xí)模型地域性適應(yīng)能力的關(guān)鍵；引入大氣頂層輻射后，MARS6（輸入大氣頂層輻射、最高氣溫、最低氣溫、相對濕度）、MARS7（輸入大氣頂層輻射、最高氣溫、最低氣溫、日照時長）、MARS8（輸入大氣頂層輻射、最高氣溫、最低氣溫）模型均優(yōu)于相同氣象因子依賴下的Irmak-Allen、Irmak、Hargreaves-M4模型；通過可移植性分析發(fā)現(xiàn)，在訓(xùn)練站點和測試站點的隨機(jī)交叉組合下，MARS5模型保持了較高的精度（納什效率系數(shù)和決定系數(shù)均大于0.985），且輸出較為穩(wěn)定的模擬結(jié)果，均方根誤差變化范圍為0.121～0.193 mm/d，平均相對誤差變化范圍為2.7%～4.2%。因此，基于多元自適應(yīng)回歸樣條算法的ET0預(yù)報模型可作為川中丘陵區(qū)ET0預(yù)報的推薦模型。

蒸散；算法；模型；多元自適應(yīng)回歸樣條；川中丘陵區(qū)；可移植性

0 引 言

參考作物蒸散量（以下簡稱ET0），又稱潛在蒸散量，是指在水分充足、生長情況良好的條件下，冠層蒸散阻力為70 s/m、反照率為0.23、高度為12 cm的草本植被完全覆蓋地面時的蒸散量[1-2]。ET0可以表征大氣蒸散能力，是農(nóng)田水資源優(yōu)化配置、農(nóng)作物水量需求補償?shù)葐栴}的重要參考量[3-4]。

ET0可以通過實際測量法和模型估算法得到，但實測難度較大且成本較高。目前，針對不同地區(qū)的ET0已有很多估算模型，其中具有較強普適性的模型是FAO 56 Penman-Monteith模型（以下簡稱P-M模型）。P-M模型充分考慮了各種氣象因素的影響[1]，基于完整的氣象數(shù)據(jù)，便能得到較為精確的ET0結(jié)果。此外，還發(fā)展出60余種簡化估計模型，包括溫度法中的Hargreaves-M4（H-M）模型[5]和McGuinness-Bordne模型[6]等，輻射法中的Irmak-Allen（I-A）模型[7]和Irmak（IK）模型[8]等。由于ET0主要受到氣溫和太陽輻射的影響[9]，因此基于氣溫和太陽輻射的Priestley-Taylor等ET0預(yù)報模型迅速發(fā)展。

中國現(xiàn)有功能較為全面的國家基本氣象站僅有756個，隨精準(zhǔn)農(nóng)業(yè)迅速發(fā)展，對精細(xì)化氣象服務(wù)的需求不斷增強，粗放型農(nóng)業(yè)氣象服務(wù)對農(nóng)業(yè)生產(chǎn)的指導(dǎo)性逐漸減弱，對低成本投入下基于較少氣象參數(shù)輸入的ET0簡化精準(zhǔn)模擬模型生產(chǎn)需求愈來愈突出。隨著機(jī)器學(xué)習(xí)算法迅速發(fā)展，較少氣象參數(shù)輸入驅(qū)動的多種ET0模擬模型相繼被提出。崔寧博等[10]將思維進(jìn)化（mind evolutionary algorithm, MEA）算法與誤差反向傳播神經(jīng)網(wǎng)絡(luò)（back-propagation neural network, BPNN）模型運用于中國西北地區(qū)的ET0預(yù)測，發(fā)現(xiàn)預(yù)報精度高于相同輸入下的經(jīng)驗?zāi)Ｐ?，?dāng)氣象數(shù)據(jù)缺省時能作為西北旱區(qū)的ET0預(yù)報模型。Abdullah等[11]將極限學(xué)習(xí)機(jī)（extreme learning machine, ELM）模型運用于伊拉克地區(qū)的ET0預(yù)測，結(jié)果表明ELM模型運行高效且精度較高，在干旱、半干旱地區(qū)具有較高的泛化能力。Tabari等[12]利用支持向量機(jī)（support vector machine, SVM）和自適應(yīng)神經(jīng)模糊推理系統(tǒng)（adaptive neural fuzzy inference system, ANFIS）進(jìn)行了ET0預(yù)報模型優(yōu)化，提高了預(yù)報精度。Feng等[13-14]將極限學(xué)習(xí)機(jī)、隨機(jī)森林（random forests, RF）模型和廣義回歸神經(jīng)網(wǎng)絡(luò)（generalized regression neural networks, GRNN）模型用于四川的ET0預(yù)測，基于溫度資料和地外輻射數(shù)據(jù)獲取了較好的模擬精度。

目前常用的時間序列分析模型（例如線性回歸模型）和神經(jīng)網(wǎng)絡(luò)模型等，不能體現(xiàn)自變量間的交互作用，且前者精度較低；后者盡管對部分因子的處理優(yōu)于線性回歸模型，但物理意義模糊，不能得到顯式表達(dá)。本文采用多元自適應(yīng)回歸樣條（multivariate adaptive regression splines，MARS）模型進(jìn)行ET0預(yù)報。MARS模型結(jié)合了樣條回歸、累加回歸等多個回歸函數(shù)的優(yōu)點，考慮了各個變量之間的交互作用，能甄別氣象因子與ET0間復(fù)雜的非線性映射關(guān)系，并對每個變量的函數(shù)關(guān)系實現(xiàn)顯式表達(dá)，有很強的自適應(yīng)性和模型解讀能力，具有極佳的推廣價值。Krzemień等[15]基于MARS提出了活動礦井的溫度預(yù)測模型，能更好地防止煤氣化火災(zāi)的發(fā)生。Rezaie-Balf等[16]通過改良MARS算法，提出了精度較高的W-MARS模型來預(yù)測地下水資源的動態(tài)變化。

本文以P-M模型計算的ET0作為標(biāo)準(zhǔn)值，建立了15種不同的氣象因子輸入組合，利用MARS模型對川中丘陵區(qū)的ET0進(jìn)行預(yù)報，將其與常用簡化經(jīng)驗?zāi)Ｐ蛯Ρ?，并分析其在川中丘陵區(qū)的可移植性，為氣象資料缺損地區(qū)提供一種ET0預(yù)報的可靠方法。

1 材料與方法

1.1 研究區(qū)域概況

四川盆地中部以丘陵地貌為主，常稱川中丘陵區(qū)。川中丘陵區(qū)是典型的方山丘陵區(qū)，東至華鎣山，西達(dá)龍泉山，北通大巴山南麓，南臨長江，約8.4萬km2；海拔在250～600 m以內(nèi)，相對高差小于100 m。因此，該地區(qū)低山廣布，河谷縱橫，地表起伏較多[17]。川中丘陵區(qū)地處濕潤、半濕潤帶，面積占四川盆地近50%，是四川重要的農(nóng)業(yè)生產(chǎn)區(qū)。其北部是山地到淺山丘的過渡區(qū)，中部是淺山丘與平原區(qū)，南部是山地與丘陵?？紤]川中丘陵區(qū)復(fù)雜的地理氣候因素，選取11個代表性站點[18]，具體分布情況見圖1，站點基本地理、氣候信息見表1。

圖1 氣象站點分布圖

表1 川中丘陵區(qū)各氣象站點基本數(shù)據(jù)

Table 1 Basic data of each meteorological station in hilly area of central Sichuan

1.2 數(shù)據(jù)獲取

在數(shù)據(jù)獲取過程中，選用了綿陽站、樂山站、宜賓站等11個代表性站點的氣象數(shù)據(jù)，具體包括日值平均溫度（mean）、最低氣溫（min）、最高氣溫（max）、10 m處風(fēng)速（10）、日照時長（sun）、相對濕度（relative humidity, RH）等。數(shù)據(jù)均來自國家氣象信息中心（https://data.cma. cn/），其中，內(nèi)江站資料為1999—2016年日值數(shù)據(jù)；其余10個站點資料均為1961—2016年日值數(shù)據(jù)（針對廣元站及達(dá)縣站缺失的部分?jǐn)?shù)據(jù)，通過插值法[19]補全）。

參考FAO推薦模型，距地面2 m處風(fēng)速由風(fēng)廓線[1]關(guān)系推出：

式中為測點到地面的垂直距離，m；u為距地面高度處風(fēng)速，m/s.

Samani[9]和Hargreaves等[20]認(rèn)為太陽輻射是引起白晝溫差的主要原因，故利用大氣頂層輻射R和溫差來推演凈輻射，彌補缺失日照時長sun的計算誤差。此外，馮禹等[21]發(fā)現(xiàn)，將大氣頂層輻射R加入廣義回歸神經(jīng)網(wǎng)絡(luò)（generalized regression neural network, GRNN）模型和小波神經(jīng)網(wǎng)絡(luò)（wavelet neural network, WNN）模型后，能提高模型預(yù)測精度。因此，本文將R作為模型的輸入項，R是與站點緯度和時間相關(guān)的量[1]。

式中Latitude為站點緯度，rad；為日序數(shù)，每年1月1日起，至12月31日結(jié)束，從1到365（或366）循環(huán)。

1.3 模型準(zhǔn)備

本文選取了6個與ET0相關(guān)的氣象因子——max、min、2、RH、sun、R，除R可以直接計算得到外，其余需觀測得到，分析ET0與所有因子之間的相關(guān)性，比較各個因子對ET0的影響大小，選取其中2~6個因子進(jìn)行組合，得到15種輸入?yún)?shù)組合，見表2。

表2 MARS模型與輸入?yún)?shù)

注：MARS(=1, 2, 3,…,15)表示15種不同輸入下的MARS模型；R、max、min、2、RH和sun各表示大氣頂層輻射、日最高氣溫、日最低氣溫、離地面2 m處風(fēng)速、相對濕度和日照時長，下同；R可通過站點緯度和日序數(shù)計算得到，不用觀測，有別于其他氣象數(shù)據(jù)。

Notes: MARS(=1, 2, 3,…,15) represents MARS models with 15 different inputs;R,max,min,2, RH andsunrepresent aerodynamic resistance, the highest daily temperature, the daily minimum temperature, wind speed 2 m away from the ground, relative humidity and sunshine duration, the same below;Rcan be calculated by the latitude and date sequence, without observation, different from other meteorological parameters.

考慮到川中丘陵區(qū)南北跨度較大，故將研究區(qū)域分為3部分，北部：廣元、萬源、巴中、閬中、綿陽；中部：達(dá)縣、遂寧、樂山；南部：宜賓、敘永。并將1961—2016年的逐日氣象數(shù)據(jù)按7:3的比例分為2部分。前39 a作為訓(xùn)練組建立MARS模型，后17 a作為測試組，驗證不同氣象因子輸入下MARS模型的預(yù)報精度。

1.4 參考作物蒸散量計算模型

P-M模型基于能量平衡方程和水蒸氣擴(kuò)散理論，考慮了各種氣象因素，物理意義明確，本文使用P-M模型的計算結(jié)果作為ET0的標(biāo)準(zhǔn)值。其表達(dá)式[1]為

按照P-M模型規(guī)定，數(shù)據(jù)周期為1~10 d時，可忽略土壤熱通量的影響[1]，而本文獲得的氣象數(shù)據(jù)以天為單位，故取土壤熱通量為0。

為驗證MARS模型在缺省氣象數(shù)據(jù)時的預(yù)測精度，本文還選取了一些在川中丘陵區(qū)精度較高的ET0經(jīng)驗?zāi)Ｐ妥鳛閷Ρ?。各模型計算公式及相關(guān)氣象參數(shù)見表3，其中R為太陽（或短波）輻射[1]，MJ/(m2·d).

表3 ET0經(jīng)驗?zāi)Ｐ图坝嬎愎?/p>

注：R為太陽（或短波）輻射，MJ·m-2·d-1.

Notes:Ris solar (or short-wave) radiation, MJ·m-2·d-1.

1.5 MARS模型

MARS是由美國的統(tǒng)計學(xué)家Jerome Friedman于1991年提出的數(shù)據(jù)分析方法。MARS模型通過把整體數(shù)據(jù)劃分為許多小區(qū)域，從而將遞歸分割和樣條擬合結(jié)合起來，得出變量間的非線性關(guān)系[22]，進(jìn)而通過廣義交叉驗證（generalized cross-validation，GCV）準(zhǔn)則，并根據(jù)擬合對象的動態(tài)特征及變量間的相互作用調(diào)整擬合路徑，可以充分?jǐn)M合不同維度的函數(shù)。利用Matlab工具開發(fā)基于氣象因子的MARS模型用以預(yù)報ET0，它源自于Jekabsons[23]開發(fā)的程序包，可從http://www.cs.rtu.lv/ jekabsons/下載。

若考慮個基函數(shù)，一階MARS模型可表示為

1）開始為只含1個常數(shù)項的基礎(chǔ)模型，借助直傳截斷過程對樣本函數(shù)進(jìn)行分割處理，并考慮變量的交互作用，不斷增加基函數(shù)數(shù)量，提高模型的精度，直至殘差平方和達(dá)到最小值或者基函數(shù)個數(shù)達(dá)到最大值，得到1個過度擬合的模型[16, 22]；

2）通過向后剪枝過程刪除貢獻(xiàn)較低的基函數(shù)，并對剩余各項的系數(shù)不斷進(jìn)行修正計算，若能保證模型精度，則刪除多余基函數(shù)[25]，否則保留基函數(shù)；

3）最后對向后剪枝過程得到的一系列模型進(jìn)行對比，通過GCV準(zhǔn)則選擇出最優(yōu)模型，當(dāng)模型的精確度上升時，GCV值下降[16, 26]。

3類輸入?yún)?shù)，最大交互數(shù)目為2時，直傳截斷過程生成基函數(shù)的簡化拓?fù)鋱D見圖2。

圖2 MARS模型基函數(shù)導(dǎo)出過程拓?fù)鋱D

1.6 模型評價方法

用納什效率系數(shù)（Nash efficiency coefficient，NSE）、均方根誤差（root mean square error，RMSE）、平均相對誤差（mean relative error，MRE）和決定系數(shù)（2）4個統(tǒng)計參數(shù)度量2個序列之間的差異大小，然后，綜合考慮每個參數(shù)，提出綜合評價指標(biāo)（comprehensive performance indicator, CPI），對模型的預(yù)報精度和可行性進(jìn)行評估[29-30]，計算方法如下：

1）參數(shù)歸一化處理

式中為原參數(shù)序列，由同一站點在同一種輸入下，運用不同模型計算的某種參數(shù)值構(gòu)成，min為序列的最小值，max為序列的最大值，為歸一化序列。

2）基于序列中位數(shù)得到每個參數(shù)的“得分”

3）計算模型的綜合評分CPI

2 結(jié)果與分析

2.1 不同因子相關(guān)性分析

6項輸入因子與ET0的相關(guān)系數(shù)見表4，所有變量的線性相關(guān)度見圖3，max、min、2和R與ET0的相關(guān)性較強，其中R與ET0相關(guān)性最大，這與馮禹等[21]在四川盆地采用機(jī)器學(xué)習(xí)模型研究ET0預(yù)報模型時得到的結(jié)論一致，R加入模型輸入層，能提高模型精確度。

圖3 模型中所有變量的線性相關(guān)系數(shù)

表4 川中丘陵區(qū)各站點6種輸入因子與ET0的相關(guān)度分析

注：相關(guān)系數(shù)絕對值介于0.8～1.0，極強相關(guān)；介于0.6～0.8， 強相關(guān)；介于0.4～0.6，中等程度相關(guān)；介于0.2～0.4，弱相關(guān)；介于0～0.2，極弱相關(guān)或無相關(guān)；**表示在1%的水平上顯著相關(guān)，下同。

Notes: Pearson correlation coefficient between 0.8-1.0, very strong correlation; 0.6-0.8, strong correlation; 0.4-0.6, moderate correlation; 0.2-0.4, weak correlation; 0-0.2, very weak correlation or no correlation; ** represents significant correlation at 1%, the same below.

MARS模型考慮了氣象因子的交互作用，所以能實現(xiàn)輸入項的替代和刪除，為保持模型的精度，需要特別注意具有特殊意義的氣象因子。從圖3中可以看出6個輸入項的相互關(guān)系，不難發(fā)現(xiàn)2較為特殊，其與ET0的相關(guān)系數(shù)較大，而與其他5個氣象因子的相關(guān)系數(shù)都極低，不同于max、min、sun、RH和R等（至少與1～3個氣象因子的相關(guān)系數(shù)較大）。2在所有輸入項中相對獨立，與其他氣象因子的關(guān)系較為復(fù)雜，這說明模型響應(yīng)對2較為敏感，推測2的缺失將大幅度降低MARS模型的ET0模擬精度。

2.2 不同氣象因子輸入下MARS模型的預(yù)報精度

15種不同輸入下，川中丘陵區(qū)不同地區(qū)MARS模型的ET0預(yù)報結(jié)果及模型精度排名見表5，模型誤差分析見圖4的箱線圖。

表5 川中丘陵區(qū)3個分區(qū)內(nèi)不同氣象因子輸入下MARS模型ET0預(yù)報精度分析

注(Note)：<0.01.

據(jù)表5可知，在川中丘陵區(qū)3個分區(qū)內(nèi)，MARS1和MARS4（缺失RH）的NSE、RMSE、MRE和2分別為0.999、0.035～0.039 mm/d、1.05%～1.15%和0.999，CPI排名并列第1，說明MARS模型能準(zhǔn)確擬合各氣象參數(shù)與ET0之間非線性的復(fù)雜映射關(guān)系；MARS9在輸入3個參數(shù)（max、min、2）時，NSE和2均高于0.92，CPI排名為第5；MARS13（輸入max、2）和MARS15（輸入R、2）在輸入2個參數(shù)時，NSE和2均大于0.9，CPI排名分別為第6和第7。

另外，從圖4可見，MARS1和MARS4具有一致的預(yù)報精度，證明僅缺失RH時，對MARS模型的預(yù)報精度沒有影響；MARS5（缺失RH和sun）和MARS2（缺失sun）精度也較高，僅次于MARS1和MARS4，預(yù)報效果極佳。因此，MARS5可作為川中丘陵區(qū)缺少氣象參數(shù)輸入下的ET0預(yù)報模型。

對比MARS3和MARS7、MARS6和MARS8、MARS11和MARS12，在較少氣象數(shù)據(jù)輸入時，RH對模型的預(yù)報精度貢獻(xiàn)并不一定是積極的，在川中丘陵區(qū)北部，增加RH作為輸入項在一定程度上提高了模型預(yù)報精度，而在中部和南部RH的引入降低了模型預(yù)報精度，主要體現(xiàn)在2的減小，說明RH對川中丘陵區(qū)中、南部ET0的影響小于北部，輸入項的多余造成了模型的過擬合，預(yù)測結(jié)果精度降低。

圖4 2000—2016年川中丘陵區(qū)MARS模型ET0預(yù)報誤差箱線圖

對比MARS8、MARS9、MARS10、MARS11和MARS12，在僅有溫度作為輸入時，模型的預(yù)報精度降低（排名為15），增加R、sun、RH對精度提升不大（排名在9以后），而引入風(fēng)速項（MARS9）后，精度明顯提高（排名為5）。這是因為其他很多氣象因子間都有一定的非線性函數(shù)關(guān)系，模型能通過學(xué)習(xí)它們的交互過程，起到一定的彌補作用，而風(fēng)速在所有氣象因子中較為獨立，風(fēng)速的缺失會導(dǎo)致映射關(guān)系中關(guān)鍵部分的缺失，所以2對模型預(yù)報精度的影響大于其他因子，這與張皓杰等[29]在西北地區(qū)的ET0預(yù)報研究中所得結(jié)論一致。同樣，2對模型的重要性在5參數(shù)輸入的3個模型中也有體現(xiàn)，MARS4和MARS2的預(yù)報精度遠(yuǎn)大于MARS3。

MARS9中只保留max時，即MARS13，發(fā)現(xiàn)預(yù)報精度變化不大，略有降低。在MARS9中去除所有溫度數(shù)據(jù)，并添加大氣頂層輻射R后，精度又回到與MARS9相當(dāng)?shù)乃剑虼?，缺失氣象?shù)據(jù)時，可以考慮將R加入模型的輸入層，提高精度。這是因為大氣頂層輻射雖僅為緯度與日序數(shù)的函數(shù)，但日照時長的影響[20]在R中也得到了體現(xiàn)。同樣，這一點在組合MARS5和MARS8、MARS6和MARS11、MARS7和MARS10中也有體現(xiàn)。

2.3 MARS模型月尺度誤差分析

表6為川中丘陵區(qū)不同地區(qū)各月日均ET0結(jié)果，3個地區(qū)的ET0計算結(jié)果在全年內(nèi)的變化趨勢接近，月分布曲線呈拋物線狀，ET0最大值出現(xiàn)在5—7月（不同站點間略有差別），ET0最小值出現(xiàn)在12月和1月。

通過2.2節(jié)的分析發(fā)現(xiàn)，氣象因子觀測數(shù)（表2）為3、2、1時，ET0預(yù)報精度最高的分別是MARS5、MARS13和MARS15，另外，氣象因子觀測數(shù)為4時，MARS模型預(yù)報結(jié)果普遍較高，精度較低的是MARS3模型。故選取MARS3、MARS5、MARS13和MARS15，比較其在川中丘陵區(qū)北部、中部、南部的ET0預(yù)報精度，各月預(yù)報結(jié)果見圖5。1）整體上MARS5模型的精度最高，月相對誤差在?1%～3%，且在ET0較小的月份（1—3月、9—11月）精度較高，在4—10月（主要作物生長期）的相對誤差為0～2.5%。2）MARS3模型的預(yù)報精度低于MARS5模型，但其各月ET0預(yù)報的相對誤差變化不大，且整體上中部地區(qū)精度更高，相比之下，北部偏大，南部偏小。3）MARS15模型的ET0預(yù)報值在1—6月整體偏大，在7—12月整體偏小，除春季（3—6月）南部地區(qū)ET0預(yù)報精度略高于北部地區(qū)以外，其他時間3條曲線較接近。4）MARS13模型的預(yù)報結(jié)果不同于其他3個模型，月相對誤差不具有明顯的地域差別，3條曲線在各月都十分接近，這是因為MARS13模型相比于MARS3、MARS5和MARS15模型，缺失了大氣頂層輻射R作為模型輸入變量。

表6 川中丘陵區(qū)不同地區(qū)各月日均ET0

圖5 川中丘陵區(qū)不同地區(qū)ET0模擬月相對誤差對比圖

MARS模型在川中丘陵區(qū)不同區(qū)域的月尺度誤差分析說明大氣頂層輻射R是機(jī)器學(xué)習(xí)模型識別地域性差別（地理位置）的重要因子，輸入變量中R的引入能提高模型的適應(yīng)能力，是模型甄別不同地理環(huán)境差異的關(guān)鍵。

2.4 MARS模型與其他經(jīng)驗?zāi)Ｐ蛯Ρ?/h3>

將MARS模型與輸入?yún)?shù)相同的經(jīng)驗?zāi)Ｐ瓦M(jìn)行對比，結(jié)果見表7。由表7可知，僅基于溫度資料預(yù)報時MARS12優(yōu)于H-M模型，基于溫度和輻射資料預(yù)報時MARS10優(yōu)于IK模型，NSE和2均有所提高而RMSE和MRE均降低；基于溫度和大氣濕度資料預(yù)報時MARS11優(yōu)于I-A模型，雖然MRE和2差別不大，但NSE提高而RMSE降低。相同氣象因子輸入下MARS模型的預(yù)報效果優(yōu)于其他經(jīng)驗?zāi)Ｐ停也煌军c的ET0模擬精度穩(wěn)定，體現(xiàn)為MARS模型的4個統(tǒng)計參數(shù)變化范圍較小。

表7 ET0 MARS模型與經(jīng)驗?zāi)Ｐ捅容^

表7中MARS模型的對比再次證明輸入層中引入R后，ET0模擬精度得到很大提升。此外，R可由測點地理位置和時間經(jīng)簡單計算得到，R的引入不對氣象觀測工作增加任何技術(shù)或經(jīng)濟(jì)上的困難，故推薦將引入R的MARS模型作為川中丘陵區(qū)的ET0預(yù)報模型

2.5 MARS模型可移植性分析

綜上，MARS5（輸入R、max、min、2）和MARS9（輸入max、min、2）都能在缺省較多氣象參數(shù)時具有較高的ET0預(yù)報精度，且MARS5考慮了基于站點位置的大氣頂層輻射影響，模型更加精確。為進(jìn)一步論證2.3所得結(jié)論，考察大氣頂層輻射對模型地域性差異的影響，并檢驗MARS5模型在川中丘陵區(qū)不同地區(qū)間的適應(yīng)能力，從10個站點中隨機(jī)選取4個（加上內(nèi)江站）作為預(yù)測組，再隨機(jī)選取5組（每組3個站點）作為訓(xùn)練組，將訓(xùn)練組的氣象資料打亂重分配，構(gòu)建15組MARS5模型，預(yù)報結(jié)果分析見表8。訓(xùn)練樣本與測試樣本比例為1∶1，均為56 a日值氣象數(shù)據(jù)。結(jié)果表明，跨站點組合時，MARS5模型NSE>0.985，RMSE為0.121～0.193 mm/d，MRE為2.7%～4.2%，2>0.985，ET0預(yù)報精度較高，且對站點的地理位置表現(xiàn)出極強的穩(wěn)定性。說明MARS5模型在氣候條件相似的地區(qū)內(nèi)具有極高的適用性，其預(yù)報能力高且穩(wěn)定，可作為川中丘陵區(qū)預(yù)報ET0的理想模型。

表8 不同站點間MARS5 ET0模型可移植性分析

3 結(jié) 論

本文選取川中丘陵區(qū)11個氣象站點的日值氣象數(shù)據(jù)作為輸入，建立了基于多元自適應(yīng)回歸樣條（multivariate adaptive regression splines，MARS）的ET0預(yù)報模型，結(jié)果表明：

1）減少1個輸入?yún)?shù)（相對濕度）的MARS4預(yù)報精度極高，與全部6個參數(shù)輸入的模型綜合評價指標(biāo)并列排名第1；繼續(xù)減少1～2個參數(shù)（日照時長、輻射），模型（MARS5、MARS9）精度稍微降低，但納什效率系數(shù)和2仍不低于0.92；甚至只有2個參數(shù)的模型（MARS13和MARS15）的納什效率系數(shù)和2均大于0.9，均可作為缺少氣象數(shù)據(jù)時川中丘陵區(qū)預(yù)報ET0的推薦模型，尤其是MARS15僅需要1個待觀測氣象因子（風(fēng)速）（另一參數(shù)輻射可計算獲得）作為輸入就能達(dá)到較高的ET0預(yù)報精度。

2）僅基于溫度資料預(yù)報時MARS12優(yōu)于H-M模型；基于溫度和大氣濕度資料預(yù)報時MARS11優(yōu)于I-A模型；基于溫度和輻射資料預(yù)報時MARS10優(yōu)于IK模型，而將這些MARS模型分別引入輻射參數(shù)后，模擬的精度得到了進(jìn)一步提升?？梢?，當(dāng)氣象數(shù)據(jù)缺失時，相較于經(jīng)驗?zāi)Ｐ?，MARS模型更適用于川中丘陵區(qū)的ET0預(yù)報。

3）大氣頂層輻射和風(fēng)速是決定MARS模型地域性適應(yīng)能力的關(guān)鍵。在訓(xùn)練站點和預(yù)測站點的組合測試下，MARS5不僅預(yù)報精度高，還具有極強的可移植性，整體上納什效率系數(shù)≥0.985，均方根誤差<0.20 mm/d，平均相對誤差<4.2%，2>0.985。當(dāng)某站點氣象數(shù)據(jù)較少或時間尺度存在缺失時，可以采用附近站點的氣象數(shù)據(jù)用于MARS模型的學(xué)習(xí)過程，仍可獲得精度較高的ET0預(yù)報結(jié)果。

本文通過對比分析找到了影響ET0的關(guān)鍵氣象參數(shù)，并分析各參數(shù)的相互作用，實現(xiàn)了氣象因子的刪除和替換，降低了模型對氣象數(shù)據(jù)的依賴性。為氣象數(shù)據(jù)缺省地區(qū)的ET0預(yù)報研究提供了一種思路。但本文設(shè)置的輸入組合有限，關(guān)于模型輸入項的對照研究還可以進(jìn)一步改善，主要體現(xiàn)在為溫度、風(fēng)速、相對濕度等氣象因子設(shè)計不同的輸入組合，進(jìn)行氣象因子在機(jī)器模型中對ET0貢獻(xiàn)率的對比研究。

[1] Allen R G, Pereira L S, Raes D, et al. Crop evapotranspiration: Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56[J]. Rome: FAO, 1998, 300(9): D05109.

[2] 龔元石. Penman-Monteith 公式與 FAO-PPP-17 Penman 修正式計算參考作物蒸散量的比較[J]. 北京農(nóng)業(yè)大學(xué)學(xué)報, 1995, 21(1)：68－75. Gong Yuanshi. Comparison of the reference evapotranspiration estimated by the Penman-Monteith and FAO-PPP-17 Penman Methods[J]. Acta Agriculturae Universitatis Pekinensis, 1995, 21(1): 68－75. (in Chinese with English abstract)

[3] 左德鵬, 徐宗學(xué), 李景玉, 等. 氣候變化情景下渭河流域潛在蒸散量時空變化特征[J]. 水科學(xué)進(jìn)展, 2011, 22(4)：455－461. Zuo Depeng, Xu Zongxue, Li Jingyu, et al. Spatiotemporal characteristics of potential evapotranspiration in the Weihe River basin under future climate change[J]. Advances in Water Science, 2011, 22(4): 455－461. (in Chinese with English abstract)

[4] 梁麗喬, 閆敏華, 鄧偉, 等. 松嫩平原西部參考作物蒸散量變化過程[J]. 地理科學(xué)進(jìn)展, 2006, 25(3)：22－31. Liang Liqiao, Yan Minhua, Deng Wei, et al. Change of reference crop evapotranspiration from west Songnen plain[J]. Progress in Geography, 2006, 25(3): 22－31. (in Chinese with English abstract)

[5] Trajkovic S. Hargreaves versus Penman-Monteith under humid conditions[J]. Journal of Irrigation and Drainage Engineering, 2007, 133(1): 38－42.

[6] Oudin L, Hervieu F, Michel C, et al. Which potential evapotranspiration input for a lumped rainfall–runoff model?: Part 2—Towards a simple and efficient potential evapotranspiration model for rainfall–runoff modelling[J]. Journal of Hydrology, 2005, 303(1/2/3/4): 290－306.

[7] Irmak S, Irmak A, Allen R, et al. Solar and net radiation-based equations to estimate reference evapotrans-piration in humid climates[J]. Journal of Irrigation and Drainage Engineering, 2003, 129(5): 336－347.

[8] Irmak S, Allen R, Whitty E. Daily grass and alfalfa-reference evapotranspiration estimates and alfalfa-to-grass evapotrans-pir-ation ratios in Florida[J]. Journal of Irrigation and Drainage Engineering, 2003, 129(5): 360－370.

[9] Samani Z. Estimating solar radiation and evapotranspiration using minimum climatological data[J]. Journal of Irrigation and Drainage Engineering, 2000, 126(4): 265－267.

[10] 崔寧博, 魏俊, 趙璐, 等. 基于 MEA-BPNN 的西北旱區(qū)參考作物蒸散量預(yù)報模型[J]. 農(nóng)業(yè)機(jī)械學(xué)報, 2018, 49(8)：228－236. Cui Ningbo, Wei Jun, Zhao Lu, et al. Reference crop evapotranspiration prediction model of arid areas of northwest China based on MEA-BPNN[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(8): 228－236. (in Chinese with English abstract)

[11] Abdullah S S, Malek M A, Abdullah N S, et al. Extreme learning machines: A new approach for prediction of reference evapotranspiration[J]. Journal of Hydrology, 2015, 527: 184－195.

[12] Tabari H, Kisi O, Ezani A, et al. SVM, ANFIS, regression and climate based models for reference evapotranspiration modeling using limited climatic data in a semi-arid highland environment[J]. Journal of Hydrology, 2012, 444: 78－89.

[13] Feng Y, Cui N, Gong D, et al. Evaluation of random forests and generalized regression neural networks for daily reference evapotranspiration modelling[J]. Agricultural Water Management, 2017, 193: 163－173.

[14] Feng Y, Peng Y, Cui N, et al. Modeling reference evapotranspiration using extreme learning machine and generalized regression neural network only with temperature data[J]. Computers and Electronics in Agriculture, 2017, 136: 71－78.

[15] Krzemień A. Fire risk prevention in underground coal gasification (UCG) within active mines: Temperature forecast by means of MARS models[J]. Energy, 2019, 170: 777－790.

[16] Rezaie-Balf M, Naganna S R, Ghaemi A, et al. Wavelet coupled MARS and M5 model tree approaches for groundwater level forecasting[J]. Journal of Hydrology, 2017, 553: 356－373.

[17] 王朕, 梁川, 趙鵬, 等. 川中丘陵區(qū)地表干濕長程相關(guān)性及影響因素研究[J]. 四川大學(xué)學(xué)報：工程科學(xué)版, 2016, 48(S1)：61－68. Wang Zhen, Liang Chuan, Zhao Peng, et al. Long-range correlation of surface dry/wet condition and its influential factors in hilly area of central Sichuan[J]. Journal of Sichuan University: Engineering Science Edition, 2016, 48(S1): 61－68. (in Chinese with English abstract)

[18] 詹存, 梁川, 趙璐. 川中丘陵區(qū)季節(jié)性干旱時空分布特征及成因分析[J]. 農(nóng)業(yè)工程學(xué)報, 2013, 29(21)：82－90. Zhan Cun, Liang Chuan, Zhao Lu. Temporal and spatial distribution characteristics and causes analysis of seasonal drought in hilly area of central Sichaun [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(21): 82－90. (in Chinese with English abstract)

[19] 杜加強, 舒儉民, 劉成程, 等. 黃河上游參考作物蒸散量變化特征及其對氣候變化的響應(yīng)[J]. 農(nóng)業(yè)工程學(xué)報, 2012, 28(12)：92－100. Du Jiaqiang, Shu Jianmin, Liu Chengcheng, et al. Variation characteristics of reference crop evapotranspiration and its responses to climate change in upstream areas of Yellow River basin[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(12): 92－100. (in Chinese with English abstract)

[20] Hargreaves G H, Allen R G. History and evaluation of Hargreaves evapotranspiration equation[J]. Journal of Irrigation and Drainage Engineering, 2003, 129(1): 53－63.

[21] 馮禹, 崔寧博, 龔道枝. 機(jī)器學(xué)習(xí)算法和 Hargreaves 模型在四川盆地ET0計算中的比較[J]. 中國農(nóng)業(yè)氣象, 2016, 37(4)：415－421. Feng Yu, Cui Ningbo, Gong Daozhi. Comparison of machine learning algorithms and hargreaves models for reference evapotranspiration estimation in Sichuan basin[J]. Chinese Journal of Agrometeorology, 2016, 37(4): 415－421. (in Chinese with English abstract)

[22] Friedman J H. Multivariate adaptive regression splines[J]. The Annals of Statistics, 1991, 19(1): 1－67.

[23] Jekabsons G. ARESLab: Adaptive regression splines toolbox for Matlab/Octave[J/OL]. 2016-05-15[2018-01-01]. http://www.cs.rtu.lv/jekabsons/regression.html.

[24] Friedman J H, Roosen C B. An introduction to multivariate adaptive regression splines[M]. Thousand Oaks, CA: Sage Publications Sage CA, 1995.

[25] 張亞軍. 基于MARS的車用鋰離子電池荷電狀態(tài)估計的研究[D]. 長春：長春工業(yè)大學(xué), 2015. Zhang Yajun. Research on Estimation of Lithium-ion Battery SOC for Electric Vehicle based on Multivariate Adaptive Regression Splines[D]. Changchun: Changchun University of Technology, 2015. (in Chinese with English abstract)

[26] Nieto P G, Torres J M, De Cos Juez F J, et al. Using multivariate adaptive regression splines and multilayer perceptron networks to evaluate paper manufactured using Eucalyptus globulus[J]. Applied Mathematics and Computation, 2012, 219(2): 755－763.

[27] Friedman J, Hastie T, Tibshirani R. The elements of statistical learning[M]. Springer Series in Statistics. New York: Springer, 2001.

[28] Feng Y, Gong D, Mei X, et al. Estimation of maize evapotranspiration using extreme learning machine and generalized regression neural network on the China Loess Plateau[J]. Hydrology Research, 2017, 48(4): 1156－1168.

[29] 張皓杰, 崔寧博, 徐穎, 等. 基于ELM的西北旱區(qū)參考作物蒸散量預(yù)報模型[J]. 排灌機(jī)械工程學(xué)報, 2018, 36(8)：779－784. Zhang Haojie, Cui Ningbo, Xu Ying, et al. Prediction for reference crop evapotranspiration in arid northwest China based on ELM[J]. Journal of Drainage and Irrigation Machinery Engineering, 2018, 36(8): 779－784. (in Chinese with English abstract)

[30] 馮禹, 崔寧博, 龔道枝, 等. 基于極限學(xué)習(xí)機(jī)的參考作物蒸散量預(yù)測模型[J]. 農(nóng)業(yè)工程學(xué)報, 2015, 31(增刊1)：153－160. Feng Yu, Cui Ningbo, Gong Daozhi, et al. Prediction model of reference crop evapotranspiration based on extreme learning machine[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(Supp.1): 153－160. (in Chinese with English abstract)

Simulation of reference crop evapotranspiration in hilly area of central Sichuan based on MARS

Chen Xuanquan1, Cui Ningbo1,2,3※, Li Jiping1, Xu Haoruo1, Liu Shuangmei4, Ma Zelong4, Le Jinhua5, Wang Jun5

(1.610065,; 2.,610066,; 3.,,712100,; 4.610072,; 5.100083,)

The reference crop evapotranspiration (ET0) is a kernel parameter for precise irrigation management of crops and agriculture efficient water use. In order to improve the accuracy of the ET0prediction in the hilly area of central Sichuan with missing meteorological data in some area, 15 different prediction models based on multivariate adaptive regression splines (MARS) were established by using different meteorological factors. The daily meteorological data of 11 representative meteorological stations from 1961 to 2016 were analyzed by the MARS models. These data were divided into training set and test set in a ratio of 7:3, and the simulation results of the MARS models were statistically evaluated using the calculation results of the FAO 56 Penman-Monteith model as a standard. In the statistical evaluation, 4 statistical parameters were obtained by the prediction sequence and the calculation result of the FAO 56 P-M model. They were root mean square error (RMSE), mean relative error (MRE), Nash efficiency coefficient (NSE), and2. The value of the index above were used to calculate a score for evaluating the prediction accuracy of the models, and rank the models based on the scores. Then the results were compared with other ET0prediction models and the applicability of the models in the hilly area of central Sichuan was evaluated by the portability analysis. The results showed that the full MARS model with 6 input parameters had the highest accuracy. Decreasing 1 input of relative humidity, the model still had the higher accuracy, ranking No 1 based on comprehensive performance indicator (CPI), which was same with the full model ranking. Reducing continually 1 input of sunshine duration still yielded the high simulation accuracy with NSE and2higher than 0.985. Further decreasing 1-2 input, the model NSE and2still were higher than 0.9. Among these models, the model with 2 inputs of radiation and wind speed was the most easy to use since the radiation could be calculated and only wind speed was required to measure. Radiation and wind speed were the keys to determine the regional adaptability of machine learning models. Radiation contained the geographic and temporal information of the site, which made it a key factor in the MARS models to distinguish the differences in geographical environment. On the other hand, radiation could compensate for the negative impact caused by the lack of sunshine duration on the prediction accuracy of the MARS models. The wind speed was more important than the other meteorological factors because the response of MARS models were more sensitive to it. Compared with the Irmak-Allen Model, the Irmak Model, and the Hargreaves-M4 Model, the MARS6, MARS7, and MARS8improve the accuracy. Under the same meteorological factors input, the MARS models had a stronger simulation ability for ET0than the existing empirical models; Through the portability analysis, the MARSmodel with 4 input parameters of radiation, maximum and minimum temperature and wind speed maintained high precision with NSE and2both higher than 0.985, RMSE 0.121-0.193 mm/d and MRE 2.7%-4.1%. In sum, the MARS model realized the deletion and replacement of meteorological factors, reduced dependence of ET0forecasting on meteorological data, and maintained a relatively high forecasting accuracy and wide applicability. The MARS was recommended as a reliable ET0prediction model in the hilly area of central Sichuan.

evapotranspiration; algorithm; models; multivariate adaptive regression splines; hilly area of central Sichuan; portability

2019-03-29

2019-07-10

“十三五”國家重點研發(fā)計劃項目（2016YFC0400206）；國家自然科學(xué)基金資助項目（51779161）；中央高?；究蒲袠I(yè)務(wù)費（XSHZ201604）；2016年四川省級財政項目（平壩丘陵區(qū)經(jīng)濟(jì)作物高效節(jié)水灌溉技術(shù)集成研究與示范）

陳宣全，主要從事節(jié)水灌溉理論與技術(shù)研究。Email：im_chenxq@126.com

崔寧博，教授，博士生導(dǎo)師，主要從事節(jié)水灌溉理論與技術(shù)研究。Email：cuiningbo@126.com

10.11975/j.issn.1002-6819.2019.16.017

P426.2

A

1002-6819(2019)-16-0152-09

陳宣全，崔寧博，李繼平，徐浩若，劉雙美，麻澤龍，樂進(jìn)華，王 軍. 多元自適應(yīng)回歸樣條算法模擬川中丘陵區(qū)參考作物蒸散量[J]. 農(nóng)業(yè)工程學(xué)報，2019，35(16)：152－160. doi：10.11975/j.issn.1002-6819.2019.16.017 http://www.tcsae.org

Chen Xuanquan, Cui Ningbo, Li Jiping, Xu Haoruo, Liu Shuangmei, Ma Zelong, Le Jinhua, Wang Jun. Simulation of reference crop evapotranspiration in hilly area of central Sichuan based on MARS[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(16): 152－160. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2019.16.017 http://www.tcsae.org