基于隨機(jī)森林算法和氣象因子的碭山酥梨始花期預(yù)報(bào)

李 德，陳文濤，樂(lè)章燕，范孝玲，孫 義，孟雅婷，楊 健

（1. 安徽省宿州市氣象局，宿州 234000；2. 河北省廊坊市氣象局，廊坊 065000；3. 安徽省碭山縣氣象局，碭山 235300）

0 引 言

碭山酥梨（Pyrusbretshneidericv. Dangshansu pear）原產(chǎn)安徽省碭山縣，是目前中國(guó)栽培面積最大的梨品種和最重要的梨樹(shù)資源之一[1-2]。截至2018年，碭山縣境內(nèi)碭山酥梨栽培面積已超過(guò)50 000 hm2。開(kāi)花是梨樹(shù)最關(guān)鍵的物候期，花期的早晚和長(zhǎng)短及開(kāi)花的質(zhì)量直接影響生產(chǎn)和觀賞的效果。碭山酥梨花期一般在早春氣溫波動(dòng)較大的3月底到4月上、中旬，始花期的早遲受氣溫影響顯著，且期間有時(shí)會(huì)發(fā)生霜凍害，既影響人工授粉、疏花等管理活動(dòng)，也影響梨花觀賞活動(dòng)的制定與準(zhǔn)備[3]。為此，準(zhǔn)確預(yù)測(cè)梨花開(kāi)始日期，是梨花觀賞活動(dòng)舉辦和花期霜凍防范等活動(dòng)的重要前提，已成為氣象服務(wù)重要內(nèi)容。目前，國(guó)內(nèi)外關(guān)于樹(shù)木花期等植物重要物候期的預(yù)報(bào)方法，除少數(shù)學(xué)者利用其他指示植物的物候期[4]和依據(jù)芽生長(zhǎng)量[5]預(yù)報(bào)目標(biāo)植物花期外，多數(shù)是基于一個(gè)或多個(gè)氣象因子與物候期之間的相關(guān)關(guān)系建立統(tǒng)計(jì)回歸模型，開(kāi)展迎春花[6]、蘋(píng)果[7]、梨[8-9]和桂花[10]花期預(yù)測(cè)并取得一定成效。然而，統(tǒng)計(jì)回歸模型由于存在經(jīng)驗(yàn)特征強(qiáng)、易出現(xiàn)過(guò)擬合以及因子之間共線(xiàn)性等問(wèn)題，限制了預(yù)報(bào)精度的提高。同時(shí)，植物的始花期與氣象因子之間是一非線(xiàn)性問(wèn)題[11-12]，由線(xiàn)性回歸模型量化具有較大的局限性。因此，尋找新的預(yù)報(bào)方法尤為重要。

近年來(lái)，隨著計(jì)算機(jī)技術(shù)的發(fā)展，機(jī)器學(xué)習(xí)方法被應(yīng)用到預(yù)測(cè)研究中，特別是隨機(jī)森林算法（Random Forest,RF）在城市需水量[13]、空氣質(zhì)量[14]、森林火災(zāi)[15]、大氣降水量[16]和小麥產(chǎn)量預(yù)測(cè)[17]以及強(qiáng)對(duì)流天氣分類(lèi)預(yù)報(bào)[18]、小麥葉片SPAD值估算[19]等不同領(lǐng)域已取得較好應(yīng)用效果。白琳等[20]和Zhang等[21]研究均揭示RF比傳統(tǒng)的多元線(xiàn)性回歸的結(jié)果更為理想，處理非線(xiàn)性和分級(jí)關(guān)系更具優(yōu)勢(shì)。Fernandez-delgado等[22]評(píng)估了179 種機(jī)器學(xué)習(xí)算法在121個(gè)數(shù)據(jù)集上的性能后認(rèn)為，RF性能最好等。

然而，目前，應(yīng)用隨機(jī)森林算法進(jìn)行植物物候期預(yù)測(cè)的研究相對(duì)較少。為此，本文篩選影響碭山酥梨花期早遲的關(guān)鍵氣象因子作為特征變量，以始花期為目標(biāo)變量，根據(jù)不同的起報(bào)時(shí)間，利用RF訓(xùn)練構(gòu)建碭山酥梨始花期預(yù)報(bào)模型進(jìn)行花期氣象預(yù)報(bào)，以期提升始花期預(yù)報(bào)時(shí)效和精度，為梨花花期管理活動(dòng)實(shí)施和觀光產(chǎn)業(yè)發(fā)展等提供技術(shù)支撐。

1 材料與方法

1.1 研究區(qū)概況

安徽省碭山縣位于蘇魯豫皖 7縣交界處，境內(nèi)年平均氣溫 14.4 ℃，年日照時(shí)數(shù) 2 120.0 h，年降水量746.7 mm，全年≥10 ℃的有效積溫4 864.4 ℃·d。土壤為近代黃泛沉積物發(fā)育而成的潮土，土層深厚，通透性好，pH值8.2左右，適宜酥梨栽培。在中國(guó)梨樹(shù)氣候區(qū)劃中，碭山縣屬于碭山酥梨種植的適宜和較適宜氣候區(qū)[23]。

1.2 數(shù)據(jù)來(lái)源

1983－2018年碭山酥梨逐年始花期資料、逐年梨樹(shù)越冬期即上年 12月到當(dāng)年梨樹(shù)始花期間的逐日日照時(shí)數(shù)、日平均氣溫、日最高氣溫、日最低氣溫和降水量等氣象數(shù)據(jù)，均來(lái)自安徽省碭山縣國(guó)家氣象觀測(cè)站。其中，始花期數(shù)據(jù)為碭山縣氣象觀測(cè)站的碭山酥梨長(zhǎng)期定位觀測(cè)地段資料；觀測(cè)地段位于碭山縣國(guó)家氣象站觀測(cè)場(chǎng)東南方向2.5 km處的安徽省碭山縣園藝場(chǎng)果園內(nèi)，果園土壤肥力中等。觀測(cè)植株為1966年定植，梨樹(shù)栽植密度為241株/hm2。觀測(cè)地段與觀測(cè)植株選擇標(biāo)準(zhǔn)和梨樹(shù)始花期觀測(cè)方法按照《農(nóng)業(yè)氣象觀測(cè)規(guī)范（下卷）》[24]執(zhí)行，其中始花期的判定標(biāo)準(zhǔn)是花序上第一批花朵開(kāi)放。

1.3 研究方法

1.3.1 隨機(jī)森林回歸模型構(gòu)建

1）特征集及目標(biāo)變量：開(kāi)花是果樹(shù)生長(zhǎng)發(fā)育進(jìn)程中重要的物候期，果樹(shù)的開(kāi)花早遲與前期氣象條件關(guān)系十分密切，特別是多數(shù)果樹(shù)在盛花期前40 d左右的平均氣溫或平均最高氣溫與開(kāi)花期早遲有較好的相關(guān)性[12,25-31]。為此，首先采用Pearson相關(guān)系數(shù)法，篩選梨樹(shù)上年度越冬期間（開(kāi)始期為12月上旬）至開(kāi)花之前（3月下旬）逐旬及其跨旬氣象要素與始花期相關(guān)程度較高（即通過(guò)信度水平P＜0.050檢驗(yàn)）的氣象因子，作為基本特征因子。逐年的始花期數(shù)據(jù)，采用日序法（即儒略日序數(shù)）轉(zhuǎn)換為數(shù)值型數(shù)據(jù)，即1月1日、1月2日、1月3日……分別為1、2、3……，其余類(lèi)推，并作為目標(biāo)集。

其次，為反映不同界限溫度的積溫效應(yīng)對(duì)碭山酥梨花芽分化進(jìn)程的影響，本文依據(jù)碭山酥梨的生物學(xué)特性[32-34]、結(jié)合環(huán)境氣候特點(diǎn)和生產(chǎn)服務(wù)經(jīng)驗(yàn)，參考相關(guān)研究成果[25-26,35-37]，選取≥0 、≥3 ℃、≥5 ℃、≥7.2 ℃4個(gè)界限溫度指標(biāo)，分別對(duì)應(yīng)梨樹(shù)越冬結(jié)束后樹(shù)液開(kāi)始流動(dòng)、花芽萌動(dòng)、打破休眠的最佳溫度和生理休眠終止的上限溫度指標(biāo)[28-30,32-33]。同時(shí)，兼顧所選預(yù)報(bào)因子距離實(shí)際始花日期應(yīng)有一定提前量，以提升預(yù)報(bào)結(jié)論的實(shí)際應(yīng)用價(jià)值，選擇各界限溫度指標(biāo)計(jì)算的終止日期為較常年始花期早6 d的3月25日。即從2月11日開(kāi)始，每個(gè)指標(biāo)分別計(jì)算逐日平均氣溫累積到3月10日、3月11日、…、3月25日，計(jì)16個(gè)不同終止日期的統(tǒng)計(jì)量；其中，≥0、≥3 ℃、≥5 ℃、≥7.2 ℃ 4個(gè)指標(biāo)分別統(tǒng)計(jì)其到16個(gè)終止日期的活動(dòng)積溫及其累積日數(shù)，計(jì) 8個(gè)特征指標(biāo)；≥3 ℃、≥5 ℃、≥7.2 ℃ 3個(gè)溫度指標(biāo)分別統(tǒng)計(jì)其到16個(gè)終止日的有效積溫，計(jì)3個(gè)特征指標(biāo)，共計(jì)11×16組變量，作為積溫效應(yīng)因子。然后，采用Pearson相關(guān)系數(shù)法，篩選出與始花期相關(guān)程度較高的特征變量作為積溫效應(yīng)特征變量。

2）逐日滾動(dòng)氣象預(yù)報(bào)模型構(gòu)建：為建立逐日氣象預(yù)報(bào)模型，實(shí)現(xiàn)始花期滾動(dòng)預(yù)報(bào)，滿(mǎn)足實(shí)際氣象服務(wù)需求。本文將開(kāi)始預(yù)報(bào)日期確定在3月10日、終止預(yù)報(bào)日期選擇在3月25日，自3月10日開(kāi)始至3月25日，期間每增加1 d、分別選取不同的特征變量進(jìn)入特征集進(jìn)行模型訓(xùn)練，累計(jì)訓(xùn)練構(gòu)建16個(gè)逐日始花期預(yù)報(bào)模型，實(shí)現(xiàn)始花期逐日滾動(dòng)氣象預(yù)報(bào)。

隨機(jī)森林是由多棵分類(lèi)回歸樹(shù)（Classification And Regression Tree， CART）構(gòu)成的組合分類(lèi)模型[38]，選定的特征變量作為特征數(shù)據(jù)與始花期數(shù)據(jù)進(jìn)行集成共同構(gòu)成隨機(jī)森林的樣本數(shù)據(jù)集，采用隨機(jī)抽樣的方法，抽取16%的樣本構(gòu)成測(cè)試集，剩余的84%的樣本構(gòu)成訓(xùn)練集。本文所用資料為1983－2018年、共計(jì)36個(gè)樣本，隨機(jī)抽取1988、1990、1994、2003、2016年5 a作為測(cè)試集，剩余的31 a作為訓(xùn)練集。然后，通過(guò)自助法（bootstrap）從原始樣本集采樣得到構(gòu)建N棵樹(shù)所需的N個(gè)子集，每次未被抽到的數(shù)據(jù)稱(chēng)為袋外數(shù)據(jù)（Out-of-Bag，OOB），用來(lái)進(jìn)行內(nèi)部誤差估計(jì)和特征變量重要性評(píng)價(jià)；生成每棵樹(shù)時(shí)，從規(guī)模為M的特征變量集中隨機(jī)選擇m個(gè)變量（m＜M），對(duì)于回歸，采用均方差作為節(jié)點(diǎn)分裂標(biāo)準(zhǔn)，遞歸執(zhí)行選取最優(yōu)分枝的操作。由于隨機(jī)森林采用樣本和特征的雙重隨機(jī)抽樣構(gòu)建決策樹(shù)，因此即使不對(duì)決策樹(shù)進(jìn)行剪枝操作也不會(huì)出現(xiàn)傳統(tǒng) CART決策樹(shù)過(guò)擬合的現(xiàn)象[39]。最后將這些樹(shù)的預(yù)報(bào)結(jié)果取平均值，作為目標(biāo)變量的預(yù)測(cè)值。本文在逐日始花期的RF算法中，最大節(jié)點(diǎn)數(shù)、最大樹(shù)深度、最小子節(jié)點(diǎn)數(shù)、模型數(shù)量分別選取為1 000、10、5、100。

3）各特征變量重要性量化評(píng)價(jià)方法：本文采用精度下降率來(lái)表征各特征變量的重要性，即對(duì)所有訓(xùn)練樣本來(lái)說(shuō)，隨機(jī)打亂某一特征變量取值，采用訓(xùn)練成的模型對(duì)測(cè)試樣本進(jìn)行再次預(yù)報(bào)，并計(jì)算其預(yù)報(bào)誤差率；若這一特征變量的袋外數(shù)據(jù)（OOB）的擬合誤差越大，說(shuō)明該特征變量越重要[38-39]。本文使用 R語(yǔ)言中的程序包來(lái)計(jì)算分析各特征變量的重要性。

1.3.2 模型精度評(píng)價(jià)與驗(yàn)證

采用 3個(gè)指標(biāo)作為評(píng)價(jià)模型擬合程度優(yōu)劣，即決定系數(shù)（Coefficient of Determination，R2）、均方根誤差（Root Mean Square Error，RMSE）和預(yù)報(bào)準(zhǔn)確率/預(yù)報(bào)誤差率（Nd）。

式中n為樣本數(shù)量，ymi表示實(shí)際值，ysi表示模擬值，和分別為實(shí)際和模擬樣本的平均值，Nr表示采用訓(xùn)練集和測(cè)試集數(shù)據(jù)利用預(yù)報(bào)模型進(jìn)行預(yù)報(bào)的始花日期與實(shí)際觀測(cè)始花日期相差在±1 d和±2 d（不包含±1 d）及其誤差在±3 d以上的年份數(shù)，Nf為采用訓(xùn)練集和測(cè)試集數(shù)據(jù)進(jìn)行預(yù)報(bào)的總年份數(shù)；當(dāng)誤差日數(shù)為±1 d和±2 d（不包含±1 d）時(shí)，計(jì)算得到的Nd稱(chēng)為預(yù)報(bào)準(zhǔn)確率（單位：%）；當(dāng)計(jì)算得到的誤差為±3 d以上時(shí)，Nd稱(chēng)為預(yù)報(bào)誤差率（單位：%）。

2 結(jié)果與分析

2.1 碭山酥梨始花期的變化特征

1983－2018年，碭山酥梨始花期平均出現(xiàn)在4月2日，最早出現(xiàn)在3月24日（2004年）、最遲出現(xiàn)在4月12日（1988年），最早與最遲日期最大相差19 d。采用峰度和偏度檢驗(yàn)法[40]，對(duì)歷年始花期出現(xiàn)的時(shí)間序列進(jìn)行正態(tài)分布性檢驗(yàn)。結(jié)果表明，始花期時(shí)間序列的峰度和偏度絕對(duì)值分別為 0.695、0.068，均小于理論峰度（2.450）和理論偏度（1.225），表明，1983－2018年碭山酥梨始花期序列呈正態(tài)分布（圖 1a），可以采用線(xiàn)性趨勢(shì)法進(jìn)行演變趨勢(shì)分析[40]。

線(xiàn)性趨勢(shì)分析結(jié)果表明，1983－2018年，碭山縣碭山酥梨始花期呈極顯著（P＜0.001）提早發(fā)生趨勢(shì)（圖1b），每10 a約提早2.750 d。

2.2 影響始花期的氣象因子篩選及其評(píng)價(jià)

2.2.1 影響始花期的旬尺度氣象要素

經(jīng)過(guò)相關(guān)性分析，在所計(jì)算的旬以及不同旬組合的氣象因子與始花期之間相關(guān)性中，通過(guò)P＜0.050信度水平檢驗(yàn)的特征變量有14個(gè)。其中，平均氣溫5個(gè)、最高氣溫4個(gè)、最低氣溫3個(gè)、日照時(shí)數(shù)2個(gè)。各特征變量的物理意義見(jiàn)表1。

圖1 1983—2018年碭山酥梨始花期分布特征與演變趨勢(shì)Fig.1 Distribution characteristics and evolution trend of Dangshansu pear in its first flowering dates from 1983 to 2018

表1 旬及跨旬尺度氣象要素與始花期相關(guān)系數(shù)Table 1 Correlation coefficients between meteorological elements and first flowering dates at ten-day scale and interten-day scale

由表1可見(jiàn)，不同旬及跨旬尺度溫度要素與始花期之間均呈極顯著（P＜0.001）負(fù)相關(guān)關(guān)系，各溫度要素的開(kāi)始時(shí)間為2月中旬，此時(shí)正值梨樹(shù)越冬中后期，表明越冬中后期及之后的溫度條件與梨樹(shù)花芽完成休眠并開(kāi)始萌芽、開(kāi)花密切相關(guān)，其中以2月中旬—3月中旬、2月中旬—2月下旬、2月中旬—3月上旬、3月下旬的平均氣溫和3月下旬平均最高氣溫與始花期相關(guān)程度最高，相關(guān)系數(shù)在?0.750～?0.650之間；3月中旬的平均最低氣溫及其他時(shí)段溫度要素的相關(guān)程度次之、相關(guān)系數(shù)在?0.600～?0.550之間。旬及跨旬尺度降水量與始花期早遲之間未有通過(guò)顯著性檢驗(yàn)的因子，表明在碭山縣梨樹(shù)自越冬以后到萌芽開(kāi)花期間，降水量的多少對(duì)梨樹(shù)開(kāi)花早遲無(wú)顯著性影響。光照條件有2個(gè)時(shí)段，即3月中旬和3月下旬的日照時(shí)數(shù)與始花期之間呈弱相關(guān)關(guān)系，相關(guān)系數(shù)分別為?0.343和?0.444。考慮到梨樹(shù)始花期預(yù)報(bào)工作的實(shí)際應(yīng)用價(jià)值，和實(shí)際預(yù)報(bào)時(shí)所選的預(yù)報(bào)因子應(yīng)為實(shí)況值以降低采用預(yù)報(bào)值或估測(cè)值進(jìn)行預(yù)報(bào)所產(chǎn)生的誤差，為此，本文篩選始花期預(yù)報(bào)模型所需特征變量的原則，一是應(yīng)為實(shí)況值，二是距離實(shí)際始花期要有一定的提前量，以利于提前部署花期低溫霜凍害防御工作和開(kāi)展授粉、疏花等管理活動(dòng)，同時(shí)，也可為梨花觀賞活動(dòng)的適時(shí)開(kāi)展提供一定的籌備時(shí)間；因此，對(duì)于部分相關(guān)程度較高的3月下旬氣象因子，不作為特征變量入選；同時(shí)，優(yōu)先選擇通過(guò)較高信度水平（P＜0.001）檢驗(yàn)的因子作為預(yù)報(bào)特征變量，共計(jì)9個(gè)，分別為 Tav-1、Tav-2、Tav-3、Tav-4、Tmin-1、Tmin-2、Tmax-1、Tmax-2、Tmax-3。

2.2.2 不同界限溫度積溫與始花期相關(guān)性

Pearson相關(guān)系數(shù)分析表明，不同積溫量和累積日數(shù)與梨樹(shù)始花期相關(guān)性，僅有6個(gè)特征指標(biāo)通過(guò)P＜0.010信度水平的顯著性檢驗(yàn)，分別為逐日平均氣溫≥0與≥3 ℃活動(dòng)積溫、逐日平均氣溫≥3 ℃、≥5 ℃和≥7.2 ℃有效積溫和逐日平均氣溫≥7.2 ℃的累積日數(shù)；每個(gè)指標(biāo)自2月11日累計(jì)到3月10日為1個(gè)變量，即3月10日為 6個(gè)特征變量，參于訓(xùn)練每日始花期氣象預(yù)報(bào)模型；自3月10日開(kāi)始，向后每增加1 d，6個(gè)特征指標(biāo)即增加1組，至3月25日，共計(jì)6×16組變量。各個(gè)特征指標(biāo)的表示符號(hào)及其意義見(jiàn)表2。

圖2為不同日期6個(gè)特征變量與碭山酥梨始花期之間的相關(guān)系數(shù)。

表2 積溫效應(yīng)特征變量Table 2 Characteristic variable of accumulated temperature effect

圖2 不同界限溫度的積溫及累積日數(shù)與始花期之間相關(guān)系數(shù)Fig.2 The correlation coefficient between accumulated temperature ,cumulative days and first flowering dates

由圖2可見(jiàn)，6個(gè)特征變量總體上呈現(xiàn)3大特征，1）各特征變量自3月10日開(kāi)始至3月25日的逐日相關(guān)程度越接近梨樹(shù)始花期相關(guān)程度越高。2）ΣTa1i、ΣTa2i、ΣTa3i三組特征變量與始花期的相關(guān)程度普遍高于 ΣTa5i、ΣTa6i、ΣTni三組特征變量；同時(shí)，ΣTa1i、ΣTa2i、ΣTa3i和 ΣTa5i四組特征變量的16個(gè)時(shí)段與始花期之間的相關(guān)程度均通過(guò)P＜0.001的顯著性檢驗(yàn)，表明，越冬后到始花前之間≥0與≥3 ℃活動(dòng)積溫及其3 ℃以上有效積溫 3個(gè)變量，相較于其他 3個(gè)變量，對(duì)碭山酥梨始花期早遲的有著重要影響。3）每個(gè)特征變量與始花期之間的相關(guān)程度均呈波動(dòng)性變化的，且自3月15日開(kāi)始到3月25日隨著距離實(shí)際花期日期的臨近其相關(guān)程度亦逐漸增強(qiáng)。為此，選擇 3月 10日－3月25日之內(nèi)，5個(gè)積溫量和1個(gè)累計(jì)日數(shù)，計(jì)6個(gè)變量，累計(jì) 16日共計(jì)96個(gè)特征變量，分別按照逐日預(yù)報(bào)模型的預(yù)報(bào)的時(shí)間先后順序，作為特征變量入選預(yù)報(bào)模型。

2.3 逐日始花期氣象預(yù)報(bào)模型及其特征變量重要性分析

2.3.1 逐日始花期氣象預(yù)報(bào)模型

以3月10日為始花期開(kāi)始預(yù)報(bào)日，進(jìn)行首個(gè)氣象預(yù)報(bào)模型構(gòu)建，之后，每增加1 d訓(xùn)練構(gòu)建1個(gè)始花期氣象預(yù)報(bào)模型的方法，累計(jì)構(gòu)建16個(gè)預(yù)報(bào)模型。

由表1和圖2，分別選取預(yù)報(bào)日之前的氣象因子作為特征變量進(jìn)行模型訓(xùn)練，其中，Tav-3、Tav-4、Tmin-2、Tmax-3四個(gè)氣象因子的計(jì)算終止日期是3月20日，雖然其相關(guān)性強(qiáng)，但在訓(xùn)練構(gòu)建3月10日—3月19日的10個(gè)逐日氣象預(yù)報(bào)模型時(shí)，為避免氣象因子為估算值或預(yù)報(bào)值，故未選擇這4個(gè)氣象因子，僅選用計(jì)算終止日在3月20日前的11個(gè)特征變量；而訓(xùn)練構(gòu)建3月20日—3月25日的6個(gè)逐日氣象預(yù)報(bào)模型時(shí)，則選擇所有通過(guò)顯著性檢驗(yàn)的15個(gè)特征變量。各逐日氣象預(yù)報(bào)模型入選的具體特征變量情況見(jiàn)表3。

根據(jù)參加各日始花期預(yù)報(bào)模型訓(xùn)練構(gòu)建的氣象因子，結(jié)合表1和圖2中的各氣象因子的相關(guān)系數(shù)，比較后發(fā)現(xiàn)，16個(gè)逐日氣象預(yù)報(bào)模型中，共計(jì)200個(gè)氣象因子，其中相關(guān)系數(shù)最小的為0.469（即因子ΣTa61）、最大的為0.789（即因子ΣTa216）。

2.3.2 逐日氣象預(yù)報(bào)模型中各特征變量的重要性評(píng)價(jià)

經(jīng)計(jì)算，得到3月10日、3月11日、…、3月25日16個(gè)預(yù)報(bào)日氣象預(yù)報(bào)模型中各特征變量的重要性排序。比較發(fā)現(xiàn)：

表3 不同模型所選特征變量Table 3 Characteristic variables from different models

1）不同預(yù)報(bào)模型的特征變量重要性排在前2位和最后2位的并不相同。3月10日－3月17日的8個(gè)預(yù)報(bào)日預(yù)報(bào)模型的特征變量重要性排在前 2位的是 Tav-2和ΣTa2i（i=1，2，…，8）；3月18日－3月25日8個(gè)預(yù)報(bào)日預(yù)報(bào)模型的特征變量重要性，除3月22日與3月23日 2個(gè)預(yù)報(bào)日預(yù)報(bào)模型排在重要性首位和次位的分別為ΣTa313、ΣTa314和Tav-2、ΣTa214外，其他日期預(yù)報(bào)模型排在前 2 位的分別是 ΣTa19、ΣTa110、ΣTa111、ΣTa112、ΣTa115、ΣTa116、Tav-2和ΣTa316。各預(yù)報(bào)日預(yù)報(bào)模型特征變量重要性排在最后2位的2個(gè)特征變量中，出現(xiàn)次數(shù)最多的是Tmax-i和ΣTni，分別出現(xiàn)12和11次；其次是ΣTa3i和Tav-i，各出現(xiàn)3次；而ΣTa5i、和Tmin-i出現(xiàn)次數(shù)最少，分別為2次和1次。

2）不同預(yù)報(bào)模型中重要性所占比例較大的特征變量個(gè)數(shù)存在差異，且臨近預(yù)報(bào)終止日、占比較大的特征變量個(gè)數(shù)趨于增多。3月15日前6個(gè)預(yù)報(bào)日預(yù)報(bào)模型中，各特征變量重要性排在前 2位的所占比例較大，自第 3位特征變量開(kāi)始其重要性快速下降；而3月16日及之后的10個(gè)預(yù)報(bào)模型中，占比例較大的特征變量個(gè)數(shù)逐漸增多，如第16預(yù)報(bào)日的預(yù)報(bào)模型各特征變量重要性在第8個(gè)特征變量才開(kāi)始迅速下降。

另外，比較各特征變量與始花期的相關(guān)系數(shù)及其在各預(yù)報(bào)模型中的重要性（表1和圖2）可以發(fā)現(xiàn)，a）各預(yù)報(bào)模型中特征變量重要性排在前 2位的氣象因子，并非是其相關(guān)性最強(qiáng)的。如3月10日－3月17日8個(gè)預(yù)報(bào)模型中，特征變量排在前2位的氣象因子Tav-2和ΣTa2i（i=1，2，…，8），其與始花期的相關(guān)系數(shù)分別為?0.667和在?0.640～?0.704之間，相關(guān)性明顯弱于重要性排在后面的特征變量ΣTa1i和ΣTa3i；而與始花期相關(guān)性最強(qiáng)的特征變量ΣTa3i，僅在3月22日、23日的預(yù)報(bào)模型中排在特征變量重要性首位。b）各預(yù)報(bào)模型中特征變量重要性排在最后2位的氣象因子，并非是相關(guān)系最弱的，如模型16中的ΣTn16的相關(guān)系數(shù)為?0.691。說(shuō)明不同時(shí)期的氣象條件對(duì)梨樹(shù)花芽發(fā)育的影響有明顯差異性和復(fù)雜性。

綜合上面分析可發(fā)現(xiàn)，3月10日－3月17日，決定梨樹(shù)始花早遲的最重要?dú)庀髼l件是2月中旬－3月上旬平均氣溫；3月18日－3月25日則是日平均氣溫≥0 活動(dòng)積溫。在具體預(yù)報(bào)實(shí)踐中，應(yīng)對(duì)重要性排在前面的特征變量給予特別關(guān)注，以避免其較小變化而引起預(yù)報(bào)精度的顯著下降。

2.4 預(yù)報(bào)模型評(píng)價(jià)及其2019年預(yù)報(bào)試驗(yàn)

2.4.1 預(yù)報(bào)模型預(yù)報(bào)準(zhǔn)確率和誤差率評(píng)價(jià)

按照公式（3）分別計(jì)算 16個(gè)預(yù)報(bào)日預(yù)報(bào)模型的訓(xùn)練集和測(cè)試集誤差范圍在±1 d、±2 d（不包含±1 d）的準(zhǔn)確率及其誤差在±3 d以上的誤差率。

從圖 3可以發(fā)現(xiàn)，各預(yù)報(bào)模型的預(yù)報(bào)準(zhǔn)確率（即誤差在2 d以?xún)?nèi)的頻率）訓(xùn)練集為92.9%、測(cè)試集為75.5%；同時(shí)，各預(yù)報(bào)模型的訓(xùn)練集與測(cè)試集誤差±2 d之內(nèi)的準(zhǔn)確頻率，呈現(xiàn)隨著預(yù)報(bào)時(shí)間距離實(shí)際開(kāi)花日期的臨近、其準(zhǔn)確率逐漸增加的特征；尤其是模型15和模型16，即預(yù)報(bào)日期在3月24日和3月25日時(shí)，測(cè)試集未出現(xiàn)±3 d以上的誤差，可見(jiàn)各氣象預(yù)報(bào)模型的預(yù)報(bào)精度較高。

圖 4給出了各預(yù)報(bào)模型的訓(xùn)練集與測(cè)試集的均方根誤差（RMSE）和決定系數(shù)（R2）。

由圖4a可見(jiàn)發(fā)現(xiàn)，各預(yù)報(bào)模型的RMSE，訓(xùn)練集在1.693～2.870之間，且呈逐漸減少趨勢(shì)；測(cè)試集的RMSE在2.240～7.237之間，除模型7、模型8和模型9三個(gè)預(yù)報(bào)模型的RMSE有所增大外，其余13個(gè)預(yù)報(bào)模型，特別是自模型8開(kāi)始，測(cè)試集的RMSE逐漸減小趨勢(shì)明顯。

圖3 逐日預(yù)報(bào)模型訓(xùn)練集與測(cè)試集誤差分布Fig.3 Error distribution of the training set and test set of the daily forecast model

圖4 各預(yù)報(bào)模型訓(xùn)練集與測(cè)試集的均方根誤差（RMSE）和決定系數(shù)（R2）Fig.4 Root mean square error (RMSE)and determination coefficient (R2) of training set and test set of each forecast model

16個(gè)預(yù)報(bào)模型訓(xùn)練集與測(cè)試集的決定系數(shù)（圖4b）平均值分別為0.891和0.701且呈現(xiàn)相似特征，除測(cè)試集在模型7、模型8和模型9時(shí)的決定系數(shù)（R2）減小外，其他時(shí)間的預(yù)報(bào)模型，訓(xùn)練集和測(cè)試集的預(yù)報(bào)結(jié)果的決定系數(shù)均呈增大趨勢(shì)，表明越接近實(shí)際開(kāi)花日期，預(yù)報(bào)模型的預(yù)報(bào)精度越高。

2.4.2 2019年預(yù)報(bào)試驗(yàn)準(zhǔn)確率

利用2019年安徽省碭山縣國(guó)家氣象站觀測(cè)到的實(shí)際氣象數(shù)據(jù)，計(jì)算得到入選各預(yù)報(bào)模型的特征變量值，驅(qū)動(dòng)各時(shí)期的預(yù)報(bào)模型，即得到2019年碭山酥梨始花期的預(yù)報(bào)日期，并與2019年碭山縣國(guó)家氣象站實(shí)際觀測(cè)的碭山酥梨始花日期相比較，得到始花期預(yù)報(bào)值與實(shí)際值的誤差（圖5）。

圖5 2019年碭山酥梨始花期各預(yù)報(bào)模型的預(yù)報(bào)值與實(shí)際值Fig.5 Predicted value and the observed value of each forecast model of Dangshansu pear at the first flowering date in 2019

由圖5發(fā)現(xiàn)，各預(yù)報(bào)模型對(duì)2019年碭山酥梨始花期的預(yù)報(bào)精度均較高，其中，模型1～模型5的預(yù)報(bào)誤差在?2 d之內(nèi)，模型6和模型7的誤差在?1 d之內(nèi)，然而自3月17日開(kāi)始（模型8）及以后8個(gè)預(yù)報(bào)模型的預(yù)報(bào)值與實(shí)況值完全一致，即在15 d前準(zhǔn)確預(yù)報(bào)出當(dāng)年始花期，可見(jiàn)，本文建立的氣象預(yù)報(bào)方法具有較高預(yù)報(bào)精度。這樣的預(yù)報(bào)結(jié)論對(duì)指導(dǎo)梨花低溫霜凍防御準(zhǔn)備工作和賞花節(jié)活動(dòng)安排等具有很高的實(shí)用價(jià)值。

3 討 論

各特征變量的重要性度量表明，16個(gè)預(yù)報(bào)模型排在前2位和最后2位的特征變量在不同預(yù)報(bào)日不并相同，其中，2月中旬－3月上旬平均氣溫出現(xiàn)最多（14次），其次是≥0活動(dòng)積溫（11次）、≥3 ℃活動(dòng)與有效積溫（各10次），出現(xiàn)次數(shù)最少的是≥5 ℃與≥7.2 ℃有效積溫和2月下旬平均最低氣溫（各1次），表明對(duì)梨樹(shù)花期早遲有著重要影響的氣象條件隨著距離實(shí)際花期的臨近是變化的，且是多個(gè)氣象條件疊加影響的結(jié)果。利用RF算法能夠處理多維特征數(shù)據(jù)并能將多因素的疊加影響綜合反映出來(lái)的特點(diǎn)[38-39]，本文基于11～15個(gè)氣象因子，實(shí)現(xiàn)了梨樹(shù)花期逐日滾動(dòng)氣象預(yù)報(bào)，經(jīng)2019年預(yù)報(bào)試驗(yàn)，表現(xiàn)出一定的預(yù)報(bào)潛力。這與傳統(tǒng)上建立單個(gè)或多個(gè)氣象因子回歸方程[4-7,9-10]的預(yù)報(bào)方法有明顯區(qū)別。

特征變量的篩選是訓(xùn)練花期預(yù)報(bào)模型的基礎(chǔ)和提升RF模型輸出精度的有效方法[17,19]。本文依據(jù)相關(guān)研究成果、生產(chǎn)服務(wù)經(jīng)驗(yàn)和碭山酥梨生物學(xué)特性及環(huán)境氣候特點(diǎn)，確定不同時(shí)間段的日平均氣溫、最高氣溫、最低氣溫和部分界限溫度積溫及其累積日數(shù)為初選特征變量，以相關(guān)程度高為原則，篩選氣象因子并進(jìn)行模型訓(xùn)練；其中，3月10－19日10個(gè)模型入選特征變量11個(gè)、3月20－25日6個(gè)模型入選特征變量15個(gè)；這為解決目前對(duì)影響碭山酥梨花期的氣象條件研究尚較欠缺的條件下，開(kāi)展始花期氣象預(yù)報(bào)提供了思路和方法。

本研究發(fā)現(xiàn)，不同時(shí)段的前期氣溫與梨樹(shù)始花期早遲密切相關(guān)，尤其是2月中旬－3月中旬平均氣溫與始花期呈極顯著負(fù)相關(guān)，而且影響始花期早遲的溫度與不同界限溫度的積溫及其累積日數(shù)越臨近始花期相關(guān)程度越高，這與部分學(xué)者揭示蘋(píng)果[7]、桃樹(shù)[27]花期與前期氣溫相關(guān)性的結(jié)論相吻合。不同時(shí)段的日照時(shí)數(shù)中，除 3月下旬與花期呈顯著相關(guān)外，其他時(shí)間僅與始花期存在弱相關(guān)關(guān)系，這與一些學(xué)者認(rèn)為的果樹(shù)花期早遲與日照條件關(guān)系相對(duì)較弱[9,35-37]的結(jié)論一致。不同時(shí)段降水量與花期之間相關(guān)性均未通過(guò)顯著性檢驗(yàn)，表明降水量多少對(duì)梨樹(shù)開(kāi)花早遲無(wú)顯著性影響；分析認(rèn)為這與碭山縣地處北亞熱帶向暖溫度帶的氣候過(guò)渡帶，淺層地下水位較高，年降水量相對(duì)豐富[3,23]，碭山酥梨越冬向萌芽開(kāi)花轉(zhuǎn)變期自身需水量較少的生物學(xué)特性有關(guān)，也佐證了果樹(shù)開(kāi)花主要受前期溫度條件所支配的結(jié)論[25-26,33-37]。本研究中入選16個(gè)逐日預(yù)報(bào)模型的特征變量均為氣溫及其氣溫基礎(chǔ)上的統(tǒng)計(jì)量，也印證了一些學(xué)者的研究與實(shí)驗(yàn)結(jié)論，即溫度是決定植物物候期最重要的氣象因子，且完成某發(fā)育期需要一定數(shù)量的積溫[28,30]。

另外，有研究指出土壤溫度[8]、需冷量[32-33]和時(shí)積溫[41]等環(huán)境參量以及樹(shù)齡和施肥與灌水等管理措施[42]也與樹(shù)木花期有關(guān)。本文是從氣象預(yù)報(bào)服務(wù)角度對(duì)區(qū)域內(nèi)總體花期進(jìn)行預(yù)報(bào)而并非單一果樹(shù)或果園，考慮到資料的易取性，并未涉及這些參量，這是本研究未來(lái)需要改進(jìn)之處。

4 結(jié) 論

1）在全球氣候變暖影響下，1983－2018年碭山酥梨始花期呈極顯著提前發(fā)生趨勢(shì)，每 10 a約提前 2.750 d（P＜0.001）。

2）碭山酥梨始花期早遲與前期平均氣溫、最低氣溫、最高氣溫以及不同界限溫度的積溫量等氣象因子密切相關(guān)，相關(guān)系數(shù)在 0.469～0.789之間；且距離始花期越近、其相關(guān)性越強(qiáng)。

3）從3月10日起報(bào)到3月25日終報(bào)，由16個(gè)氣象預(yù)報(bào)模型實(shí)現(xiàn)逐日滾動(dòng)氣象預(yù)報(bào)，各預(yù)報(bào)模型的預(yù)報(bào)精度均較高，其訓(xùn)練集與測(cè)試集誤差在±2 d以?xún)?nèi)的平均正確率分別為92.9%、75.5%；在2019年試驗(yàn)預(yù)報(bào)中，提前15 d正確預(yù)報(bào)出當(dāng)年始花期。