基于氣象-生理的夏玉米作物系數(shù)及蒸散估算

王振龍，范 月，呂海深，許瑩瑩，董國強

（1. 安徽?。ㄋ炕次┧茖W研究院 水利水資源安徽省重點實驗室，蚌埠 233000；2. 河海大學理學院，南京 210098；3. 河海大學水文水資源學院，南京 210098）

0 引 言

蒸散是植物蒸騰量與植株間土壤蒸發(fā)量之和，是土壤-植物-大氣系統(tǒng)（Soil-Plant-Atmosphere Cotinuum，SPAC）中重要的水交換過程[1]。其對作物的生長發(fā)育情況有重要影響，在大氣循環(huán)中扮演重要角色，還能進一步調(diào)節(jié)氣候[2-3]。構(gòu)建精準的作物系數(shù)模型對準確預測蒸散量及進一步制定精準的灌溉計劃有重要意義。

目前針對作物蒸散量計算方法，主要方法是作物系數(shù)法。世界糧農(nóng)組織（Food and Agriculture Organization of the United Nations，F(xiàn)AO）提出了一種估算作物蒸散量的標準方法（參考作物蒸散量乘以作物系數(shù)），其中參考蒸散通過修正的彭曼（Penman-Menteith，P-M）模型計算得到，目前該方法已得到廣泛運用[4]，而其中的關(guān)鍵問題就是確定作物系數(shù)。然而，根據(jù)FAO 推薦的4 個生長階段的作物系數(shù)計算得到的蒸散量與實際蒸散可能存在較大偏差[4-5]，因此國內(nèi)外學者建立了多種作物系數(shù)模型。如Steele 等[6]利用回歸分析建立了作物系數(shù)與播后天數(shù)間的數(shù)學關(guān)系式，生成了玉米生長期內(nèi)作物系數(shù)的動態(tài)變化曲線；王娟等[7]利用氣象數(shù)據(jù)和P-M 模型構(gòu)建了夏玉米蒸散模型，并得出了不同生育期的作物系數(shù)；雷志棟等[8]采用氣象系列數(shù)據(jù)，通過P-M 模型和非水分脅迫條件下作物騰發(fā)（Non-Water Stress Evapotranspiration，NWSE）模型模擬計算了冬小麥返青-收獲期逐日作物系數(shù)；劉海軍和康躍虎[9]在冬小麥拔節(jié)抽穗期探討了農(nóng)田實際蒸散量和參考作物蒸散量的關(guān)系，利用回歸分析建立了葉面積指數(shù)和株高與作物系數(shù)的數(shù)學關(guān)系；王振龍等[10]基于溫度效應對作物生長發(fā)育的影響從而計算作物系數(shù)。盡管這些模型取得了一定成果，但綜合考慮作物系數(shù)影響要素的具體計算模型較少，且大多都是通過傳統(tǒng)方法計算擬合。本研究嘗試采用地溫及葉面積指數(shù)構(gòu)建了氣象-生理雙函數(shù)乘法模型，并融合了機器學習算法（涉及多領(lǐng)域交叉學科），進一步提高了作物系數(shù)計算精度。

本研究采用五道溝水文實驗站大型稱重式蒸滲儀實測蒸散數(shù)據(jù)，結(jié)合系列氣象資料計算出實際作物系數(shù)，將作物系數(shù)與葉面積指數(shù)利用米氏方程進行擬合，得出受生物因子影響的作物系數(shù)，進而利用多元回歸篩選出影響作物系數(shù)最密切的氣象因子，即地溫（0 cm），進行了指數(shù)擬合，最后通過梯度下降法優(yōu)化求取模型參數(shù)，采用實測數(shù)據(jù)對計算結(jié)果進行評估，對掌握夏玉米生育期內(nèi)作物系數(shù)的動態(tài)變化特征，精確估算作物實際蒸散量具有重要意義。

1 材料與方法

1.1 試驗區(qū)概況

試驗研究在安徽淮北平原五道溝水文實驗站進行。該地區(qū)（33°09’N，117°21’E）屬暖溫帶半濕潤季風氣候區(qū)，四季分明，雨熱同期。地下水位變化范圍為1～3 m，種植作物主要有玉米、小麥、花生和大豆。試驗區(qū)土壤主要包括2 種：砂姜黑土和黃潮土，分別約占淮北平原土壤的54%和33%。試驗區(qū)土壤田間持水率為28%～30%，適宜作物生長的土壤含水率18%～25%，凋萎含水率為10%～13%。試驗區(qū)設(shè)有氣象觀測場，可獲取地溫、空氣溫度及濕度、降雨量、風速及風向、水面蒸發(fā)量和日照時數(shù)等60 余年不間斷長系列氣象要素。其中地溫有0、5、10、15、20、40、80、160 和320 cm 長系列資料。另外還設(shè)有自動高精度氣象站，可獲取空氣溫度、空氣濕度、風速、凈輻射、太陽輻射、土壤熱通量等，可用于計算參考作物蒸散量。

1.2 試驗設(shè)計與方法

本研究以淮北平原分布較為廣泛的砂姜黑土區(qū)為研究對象，作物為夏玉米（登海618），于2018 年6 月22日播種，10 月8 日收獲。玉米實際蒸散量由大型稱重式蒸滲儀測得，蒸滲儀土柱口徑為4.0 m2，高度為4.0 m，蒸滲儀地下水埋深設(shè)為1 m，土壤為分層回填土，10、30、50 和100 cm 埋深處分別設(shè)有土壤水分、溫度、電導3 個參數(shù)傳感器，數(shù)據(jù)每10 min 獲取1 次。蒸滲儀自動采集數(shù)據(jù)資料時段選取2018 年6 月22 日至10 月8 日。選取同期氣象觀測場水文氣象要素數(shù)據(jù)。利用葉面積儀測定蒸滲儀內(nèi)玉米不同生育階段的葉面積，結(jié)合實際種植密度計算葉面積指數(shù)。

1.3 蒸滲儀數(shù)據(jù)處理

由于蒸滲儀會受天氣、儀器自身因素等的影響，需對數(shù)據(jù)進行預處理：1）對蒸滲儀蒸散量異常值進行剔除和插值處理；2）對1 h 內(nèi)實際蒸發(fā)筒重和10 min 差值分別平均再求和計算整點蒸發(fā)筒重，計算出時變化量，再通過移動加權(quán)平均法修正1 h 變化量，得出時蒸散量；3）對小時蒸散量累計計算日蒸散量。

1.4 參考作物蒸散量

本研究采用高精度氣象站數(shù)據(jù)，利用FAO 修正的P-M 公式計算逐日參考作物蒸散量（Reference Evapotranspiration，ET0，mm/d），如式（1）[11-14]所示

式中Rn為地表凈輻射，MJ/(m2·d)；G 為土壤熱通量，MJ/(m2·d)；U2為地面上方2 m 處風速，m/s；es為飽和水汽壓，kPa；ea為實際水汽壓，kPa；Δ 為飽和水汽壓隨溫度變化曲線的斜率，kPa/℃；γ 為干濕球常數(shù)，kPa/℃；T為2 m 高處日平均氣溫，℃，由日最高氣溫和最低氣溫平均值計算得出。

1.5 作物系數(shù)計算

本研究所用測筒控制地下水埋深為1 m，屬于不充分供水，這種情況下作物實際蒸散量（Evapotranspiration，ET，mm/d）除了受環(huán)境因素和作物自身特性的影響外，還受到土壤含水率的制約，因此，本研究引入了土壤水分脅迫系數(shù)Ks，并進行實際作物系數(shù)的計算[15]，如式（2）所示

式中Kc為作物系數(shù)。

1.6 土壤水分脅迫系數(shù)計算

土壤水分脅迫系數(shù)反映了土壤水分對蒸散量的影響，通過式（3）[16]確定為

式中θ(0-40)為0～40 cm 的平均土壤質(zhì)量含水率，由蒸滲儀內(nèi)10 cm 和30 cm 體積含水率平均計算所得，土壤容重為1.4 g/cm3；θw為凋萎含水率，θw= 10%；θf為田間持水率，θf= 28%。

2 作物系數(shù)模型建立

2.1 模型構(gòu)建

2.1.1 基本模型

作物系數(shù)可以反映作物生長狀況。當前作物系數(shù)的計算一般基于作物生理特性，但以往研究表明：作物系數(shù)受到作物生理特性及氣象因子的共同影響。因此，本研究綜合考慮兩種類型因素對作物系數(shù)（Crop Coefficient，Kc）的影響[17]，通過建立氣象-生理雙函數(shù)乘法模型計算作物系數(shù)，模型式如式（4）所示

式中Ke為由氣象因素構(gòu)成的作物系數(shù)函數(shù)，KLAI為由作物生理因素構(gòu)成的作物系數(shù)函數(shù)，其中LAI（Leaf Area Index)為葉面積指數(shù)，m2/m2。

當前，不同研究者針對作物系數(shù)與葉面積指數(shù)間關(guān)系采取了不同擬合方式，如高斯擬合、指數(shù)擬合等。在生長初期，作物快速生長，蒸散隨作物蒸騰量的增加而增加，作物系數(shù)也增大，當增大到一定程度時，LAI 與作物系數(shù)的關(guān)系不再明顯，作物系數(shù)主要受環(huán)境因子影響，本研究根據(jù)作物系數(shù)對葉面積指數(shù)的響應特征，遂采用米氏（Michaelis-Mentenn）方程擬合作物系數(shù)與LAI[18-19]，即KLAI為式（5）所示

式中m 和n 為待擬合參數(shù)；KLAI為僅受生物因子影響時的作物系數(shù)。

將實際Kc去除生物因素影響的處理，即用實際Kc除以擬合得到的KLAI，得Ke，分析其與同期各氣象因子包括日平均氣溫、日最高氣溫、日最低氣溫、日照時數(shù)、風速、絕對濕度、相對濕度、飽和差、水汽壓力差和0、5、10、20、40、80、160、320 cm 地溫的相關(guān)性，采用多元逐步回歸法篩選出解釋程度較高的因子，并對各要素進行共線性診斷，根據(jù)SPSS 分析結(jié)果可知各個系數(shù)的方差膨脹因子（Variance Inflation Factor，VIF）較?。?

式中a、b、c 為未知參數(shù)；D0為地表溫度，℃。

表1 各函數(shù)擬合優(yōu)度和均方根誤差比較

2.1.2 梯度下降法

針對上述模型中的未知參數(shù)，本研究利用梯度下降法確定參數(shù)，梯度下降算法是機器學習中較為廣泛的優(yōu)化算法[20-22]，本研究主要以線性回歸算法損失函數(shù)求極小值確定參數(shù)a、b、c，因此需將式（6）轉(zhuǎn)化為線性形式，其步驟如下：

1）首先對式（6）進行轉(zhuǎn)化，對其兩邊取對數(shù)，得式（7）為

式中θ0、θ1和θ2為未知參數(shù)；x0為特征常量設(shè)為1；x1為地溫，℃；x2為地溫的平方，℃2。

為了算法能更快收斂，要確保不同特征的取值在相近的范圍內(nèi)，需要對自變量進行特征歸一化。式（8）中D0和D02取值范圍相差較大，需要對其進行特征縮放，即為式（10）和式（11）所示

式中x1'為均值歸一化后的地溫，℃；x1av為地溫平均值，℃；x1max為地溫最大值，℃；x1min為地溫最小值，℃；x2'為均值歸一化后的地溫平方，℃2；x2av為地溫平方平均值，℃2；x2max為地溫平方最大值，℃2；x2min為地溫平方最小值，℃2。

因此將式（9）轉(zhuǎn)化為式（12）所示

2）對式（12）利用梯度下降法求得參數(shù)；在線性回歸算法中，目標損失函數(shù)定義為式（13）所示

式中J(θ)為損失函數(shù)；θ 為由θ0、θ1和θ2構(gòu)成的參數(shù)向量；m 為樣本總數(shù)；xi為由x0、x1'及x2'構(gòu)成的特征向量；yi為實際作物系數(shù)構(gòu)成的向量；hθ(xi)為對應預測值構(gòu)成的函數(shù)。具體過程如下：對θ 賦初值，這個值可以隨機選取，也可以讓其為一個全0 的向量，本研究令其初始值為0 向量；更新θ 的值，使得目標損失函數(shù)按梯度下降的方向減少，通過式（14）迭代

式中α 為學習率，需人為指定，本研究取值0.001，α 若過大會導致震蕩無法收斂，若過小收斂速度會很慢；對多元函數(shù)的參數(shù)求偏導，將求得的各個參數(shù)的偏導數(shù)以向量的形式表示就是對應的梯度，本式中為對應的梯度，梯度計算結(jié)果如式（15）所示

3）當每次更新后的結(jié)構(gòu)都能讓損失函數(shù)變小，最終達到最小則停止下降。

2.2 模型評價指標

本研究評價指標包括平均絕對誤差（Mean Absolute Deviation，MAE）、均方誤差（Mean-Square Error，MSE）、均方根誤差（Root Mean Square Error，RMSE）和相關(guān)系數(shù)（Correlation Coefficient，r），用來評價計算結(jié)果與實測數(shù)據(jù)間的誤差情況；采用準確率（Precision，P）對模型精度進行評價，對于作物系數(shù)，定義絕對誤差在0.2、0.3 和0.4 以內(nèi)的數(shù)據(jù)個數(shù)占總數(shù)據(jù)集的比例為相應準確率，分別記為P0.2、P0.3和P0.4；對于蒸散量，定義絕對誤差在2、3 和4 mm 以內(nèi)的數(shù)據(jù)個數(shù)占總數(shù)據(jù)集的比例為相應準確率，分別記為P2、P3和P4。各指標通過式（16）～式（19）計算

式中yi為實際值；i為模型計算值；i 為樣本序數(shù)，i=1,2,…m；為實際值均值；為計算值均值；m 為樣本個數(shù)。

3 結(jié)果與討論

3.1 作物系數(shù)模型

夏玉米全生育期地溫、葉面積指數(shù)及實際作物系數(shù)變化（圖1）。地溫、LAI 與Kc的r 分別為0.75 和0.52，相關(guān)性較高。作物根系生長對地溫變化非常敏感，地溫越高，作物的生長發(fā)育越快，作物蒸騰作用增強，土壤水的移動越頻繁；地溫低，根系的生長和代謝都會受到抑制，作物根系吸水緩慢，蒸騰作用較弱。且作物蒸騰主要通過葉片蒸騰作用，與葉面積關(guān)系密切[23-26]。

圖1 全生育期地溫和作物系數(shù)Kc 變化過程線圖 Fig.1 Process of changes in ground temperature and crop coefficient Kc during the whole growth period

利用蒸滲儀實測蒸散數(shù)據(jù)及氣象系列資料對作物系數(shù)模型中各參數(shù)進行估計。首先對于KLAI，選取2018 年玉米生長季（6 月22 日至10 月8 日）的實測葉面積指數(shù)數(shù)據(jù)進行研究，每隔1 d 抽取1 d 數(shù)據(jù)作為訓練樣本（55 d），其余作為檢驗樣本（54 d）。將根據(jù)式（2）計算得到的作物系數(shù)Kc與實測葉面積指數(shù)LAI 進行擬合，利用Matlab擬合為式（5）的形式，所得方程如式（20）所示

式中0.158 4 為米氏常數(shù)，數(shù)值越小，作物系數(shù)對LAI 變化敏感性越強；由于時間尺度為日尺度，故LAI 均對應逐日作物系數(shù)Kc，故R2較小，但兩者之間具有顯著相關(guān)性（P<0.01）。

其次針對Ke，利用梯度下降法（如式13～式15）解得式（12）中的參數(shù)θ0、θ1和θ2分別為-0.063 4、0.542 8和0.478 8，再根據(jù)所得θ0、θ1和θ2的值求得式（6）中a、b 及c 的值分別為0.804 7、-0.566 8 和1.445 2。所得Ke計算式如式（21）所示

最后，根據(jù)式（20）及式（21）建立的乘法模型式如式（22）所示

圖2 為訓練樣本根據(jù)式（22）所得作物系數(shù)計算值及實際值變化過程圖。兩者平均值分別為0.399 和0.405，相差小于0.01，MAE 為0.11，MSE 為0.02，RMSE 為0.14，對于作物系數(shù)估算準確率P0.2為80%、P0.3為98%、P0.4為100%，誤差較小，說明模型精度較高。

圖2 訓練樣本計算作物系數(shù)Kc 與實際值比較 Fig. 2 Calculated crop coefficient Kc compared with actual value in training sample

3.2 模型檢驗

圖3 為檢驗樣本根據(jù)式（22）所建模型計算得到的作物系數(shù)Kc及實際值變化過程圖。兩者平均值分別為0.40 和0.44，僅相差0.04，MAE 為0.12，MSE 為0.02，RMSE 為0.15，對于作物系數(shù)估算準確率P0.2為81%、P0.3為96%、P0.4為100%，誤差較小，模型預測精度較高。

圖3 檢驗樣本計算作物系數(shù)Kc 與實際值比較 Fig.3 Calculated crop coefficient Kc compared with actual value in testing sample

全生育期夏玉米實際作物系數(shù)與計算作物系數(shù)均值分別為0.42 和0.40，僅相差0.02，MAE 為0.12，RMSE為0.15，r 為0.91。對于作物系數(shù)估算準確率P0.2為81%、P0.3為97%、P0.4為100%。

表2 為以實際作物系數(shù)為依據(jù)，不同樣本范圍模型的計算精度評價指標值。結(jié)果表明本研究構(gòu)建的作物系數(shù)模型能以較高精度計算夏玉米生長期內(nèi)逐日作物系數(shù)，再利用聯(lián)合國糧農(nóng)組織（FAO）所推薦的估計作物實際蒸散量方法，可以進一步計算作物逐日蒸散量ET。

表2 不同樣本作物系數(shù)計算精度評價指標值 Table 2 Evaluation index value of crop coefficient calculated accuracy in different samples

3.3 蒸散量計算及精度評價

圖4 為依據(jù)式（2）計算所得訓練集蒸散值與實際值對比圖。兩者平均值分別為3.88 和3.86 mm/d，相差小于0.1 mm/d，MAE 為1.0 mm/d，RMSE 為1.4 mm/d，對于蒸散量估算準確率P2為80%、P3為95%、P4為98%，誤差較小，預測精度較高。

圖5 為依據(jù)式（2）計算所得驗證集蒸散值與實際值對比圖。兩者平均值分別為4.01 和4.16 mm/d，僅相差0.15 mm/d，MAE 為1.0 mm/d，RMSE 為1.3 mm/d，對于蒸散量估算準確率P2為81%、P3為96%、P4為100%，誤差較小，模型預測精度較高。

圖5 檢驗樣本計算蒸散量ET 與實際值比較 Fig.5 Calculated evapotranspiration ET compared with actual value in testing sample

表3 為以實際作物系數(shù)Kc為依據(jù)，不同樣本范圍模型的計算精度評價指標值。結(jié)果表明本研究構(gòu)建的作物系數(shù)模型能以較高精度計算夏玉米生長期內(nèi)逐日作物系數(shù)，并利用 FAO 所推薦的估計作物實際蒸散量（Evapotranspiration，ET）方法所計算的作物逐日蒸散量ET 精度也較高。

為防止出現(xiàn)過擬合問題，本研究選取夏玉米6 月22日至7 月22 日連續(xù)時間作模型驗證，所得蒸散量ET 實際值與估算值如圖6 所示。兩者平均值分別為2.33 和2.54 mm/d，僅相差0.21 mm/d，MAE 為0.83 mm/d，RMSE為2.82 mm/d，r 為0.85，精度較高。

表3 不同樣本蒸散量計算精度評價指標值 Table 3 Evaluation index value of evapotranspiration calculated accuracy in different samples

圖6 連續(xù)時段計算蒸散量ET 與實際值比較 Fig.6 Calculated evapotranspiration ET compared with actual value in continuous period

3.4 模型比較與評價

本研究采用相同數(shù)據(jù)并根據(jù)FAO 推薦作物系數(shù)及P-M 模型計算了全生長季的玉米蒸散量。FAO 建議將作物的生長期劃分為生長初期、發(fā)育期、生長中期和生長后期這4 個階段[27-30]。根據(jù)試驗地區(qū)作物實際生長情況，夏玉米各生長階段劃分及對應的推薦作物系數(shù)如表4，其中發(fā)育期作物系數(shù)通過插值得到。

表4 夏玉米生長階段劃分及對應作物系數(shù) Table 4 Summer maize growth stage division and corresponding crop coefficient

將計算結(jié)果與本研究作物系數(shù)模型計算的玉米蒸散量分別與其實際蒸散量進行了對比（圖7），發(fā)現(xiàn)三者的趨勢基本一致，但根據(jù)FAO 推薦作物系數(shù)計算的蒸散量ET 值整體偏高，計算精度明顯低于本研究所建作物系數(shù)模型，具體計算結(jié)果如表5 所示。

基于本研究模型計算所得蒸散值與實際值均值僅相差0.07 mm/d，MAE 為1.00 mm/d，RMSE 為4.50 mm/d，r 為0.75，對于蒸散量估算準確率P2為81%、P3為95%、P4為99%。根據(jù)FAO 推薦作物系數(shù)計算所得蒸散值與實際值存在較大偏差，兩者均值相差5.00 mm/d，MAE 為5.07 mm/d，RMSE 為6.09 mm/d，r 為0.75，對于蒸散量估算準確率P2為23%、P3為30%、P4為36%。本研究模型計算的蒸散值準確率相比FAO 推薦作物系數(shù)得到的計算值準確率（誤差在2～3 mm/d 以內(nèi)）提高了3 倍以上，且誤差幾乎都在4 mm/d 以內(nèi)。因此本研究構(gòu)建的氣象-生理雙函數(shù)模型，綜合反應了氣象和生物因子對玉米生長的共同影響，可以更高精度估算玉米作物系數(shù)和蒸散量。

表5 全生育期不同模型蒸散量計算精度評價指標值 Table 5 Evaluation index value of evapotranspiration calculated accuracy by different models during the whole growth period

圖7 全生育期不同方法蒸散量ET 計算值與實際值比較 Fig.7 Calculated evapotranspiration ET compared with actual value based on different ways during the whole growth period

4 結(jié) 論

1）本研究基于地溫和葉面積指數(shù)構(gòu)建了氣象-生理雙函數(shù)乘法模型，可以反映作物逐日作物系數(shù)變化。該模型在訓練和檢驗樣本中計算精度均較高。在訓練樣本中，作物系數(shù)計算值與實際值均值相差0.01，平均絕對誤差（Mean Absolute Deviation，MAE）為0.11，均方根誤差（Root Mean Square Error，RMSE）為0.14，相關(guān)系數(shù)（Correlation Coefficient，r）為0.94；在檢驗樣本中，兩者均值相差0.04，MAE 為0.12，RMSE 為0.15，r 為0.87。該模型簡化了夏玉米作物系數(shù)的計算，明確了地溫和葉面積指數(shù)對作物系數(shù)的綜合影響程度，提高了計算精度，可用于玉米作物系數(shù)的動態(tài)計算。

2）基于本研究作物系數(shù)模型計算的實際蒸散量精度較高，在訓練樣本中，蒸散量計算值與實際值均值相差0.1 mm/d，MAE 為1.0 mm/d，RMSE 為1.4 mm/d，r 為0.75；在檢驗樣本中，兩者均值相差0.15 mm/d，MAE 為1.0 mm/d，RMSE 為1.3 mm/d，r 為0.76。與聯(lián)合國糧農(nóng)組織（FAO）推薦作物系數(shù)計算結(jié)果相比，精度更高。

3）根據(jù)FAO 推薦作物系數(shù)計算蒸散值與實際值兩者均值相差5.0 mm/d，MAE 為5.07 mm/d，RMSE 為6.09 mm/d，r 為0.75，本研究模型計算的蒸散值準確率相比FAO推薦作物系數(shù)得到的計算值準確率（誤差在2～3 mm/d 以內(nèi)）提高了3 倍以上。

本研究作物系數(shù)模型是基于五道溝水文實驗站大型稱重式蒸滲儀實測數(shù)據(jù)構(gòu)建，本模型研究時間及空間尺度還有一定的局限性，通過蒸散發(fā)、潛水蒸發(fā)和包氣帶土壤水運移等多要素、多過程、面尺度、時尺度原型監(jiān)測試驗有待進一步完善，以大型蒸滲儀與大田地下水位同步變化的蒸散量試驗有待進一步研究和實施完善。