林馨貝，周崗，鄭澤濤，趙璐,2*，梁川

（1.四川大學(xué)，成都610065；2.南方丘區(qū)節(jié)水農(nóng)業(yè)研究四川省重點實驗室，成都 610066）

0 引言

【研究意義】農(nóng)田作物蒸散發(fā)（ET）包括土壤蒸發(fā)和植株蒸騰2 個部分[1]，是農(nóng)田能量平衡和水分平衡的重要組成部分[2]。精確估算作物ET可為制定農(nóng)田灌溉方案、評估水資源規(guī)劃管理提供依據(jù)?！狙芯窟M展】如今ET測量方法較多，如可通過渦度相關(guān)方法、大型蒸滲儀法、水文學(xué)法等進行直接測定，也可通過遙感法、波文比能量平衡法、空氣動力學(xué)法等間接測定[3-5]。隨著測量精度不斷提高，多種ET估算模型被相繼提出，其中應(yīng)用較多的是單源 Penman-Monteith 模型（ P-M ）、雙源Shuttleworth-Wallance 模型（S-W）[6-9]、多源Clumping模型（C）[10-12]。

從理論上分析，多源復(fù)雜模型能更明確地反映ET過程和ET水分的來源，但多源復(fù)雜模型參數(shù)較多，容易出現(xiàn)“過度擬合”的情況導(dǎo)致參數(shù)擬合精度下降。同時，完整獲得多源復(fù)雜模型所需的全部數(shù)據(jù)具有一定難度，數(shù)據(jù)中測量誤差的累積也會影響模型精度，因此多源復(fù)雜模型的估算精度不一定高于單源模型[13-14]。在單源模型中，P-M 模型由于反映了各氣候要素的綜合影響，具有明確的物理意義，是ET研究中應(yīng)用最為普遍的簡單模型[15]。

P-M 模型估算農(nóng)田ET的精度主要取決于冠層阻力（rc）的準確性。rc是P-M 模型中的重要參數(shù)，其受到葉片位置、土壤濕潤狀況及冠層內(nèi)空氣動力學(xué)特性等因素的綜合影響，無法通過儀器直接測定[16]。目前rc計算的有效方法有3 種，第一種是利用實測ET值按P-M 或S-W 模型反推rc[17-18]；第二種是利用單葉氣孔阻力，結(jié)合作物群體葉面積指數(shù)的空間垂直分布，計算得到群體的rc[19-21]；第三種是利用模型計算rc[22-23]。3 種方法中，利用模型法計算rc的應(yīng)用較多，其中Jarvis 模型（JA 模型）全面考慮了太陽輻射、水汽壓差、氣溫、土壤含水率及葉面積指數(shù)等因子與冠層阻力的綜合關(guān)系，包含了各環(huán)境響應(yīng)函數(shù)[24]，應(yīng)用較為廣泛。Srivastava 等[25]采用4 種方法估算rc，然后根據(jù)P-M 模型計算半濕潤地區(qū)灌溉玉米的ET，得出JA 模型具有較高的可靠性；李彩霞等[26]基于JA 模型，引入冠層上、下部溫差變量（ΔTc），提出了交替隔溝非充分供水下玉米rc的計算方法；周惠萍等[27]考慮氮元素對冠層導(dǎo)度的影響，建立了基于氮的JA 修正模型，提高了P-M 模型在不同施氮水平下估算番茄ET的準確性；李思恩等[28-29]基于JA 模型和流體運動的阻力定律，增加了土壤阻力，構(gòu)建了耦合的表層阻力模型（CO 模型），CO 模型所求rc值可被看作水汽通過冠層和土壤的所克服的阻力平均值，該模型在作物生長季前期（LAI＜2）時模擬效果較好[30]。

【切入點】JA 模型和CO 模型中都采用了一些經(jīng)驗參數(shù)，由于不同區(qū)域環(huán)境因素和作物品種存在一定差別，因此在利用模型計算rc時，需先進行參數(shù)優(yōu)化，以提高模擬精度。目前關(guān)于玉米ET研究中，大部分研究只選擇多種rc的計算模型進行比較，或采用單一參數(shù)優(yōu)化方法對rc的計算模型進行參數(shù)優(yōu)化[25-26]，采用不同參數(shù)優(yōu)化方法對同時對多種模型進行參數(shù)優(yōu)化的研究較少，因此無法通過結(jié)果比較來篩選地區(qū)適用性最好的模型和參數(shù)優(yōu)化方法的組合。

【擬解決的關(guān)鍵問題】懷來位于半干旱氣候區(qū)，光熱充足，適合夏玉米的生長[31]。夏玉米生長過程中耗水較多[32]，在干旱條件下易因缺水而減產(chǎn)[33-34]，因此需精確計算夏玉米ET以通過灌溉來保證玉米的正常生長發(fā)育。為此，本文基于懷來站點2013年的渦度相關(guān)數(shù)據(jù)和氣象數(shù)據(jù)，以P-M 模型反推的實際rc為標準，分別采用最小二乘法和蟻群算法率定JA 模型和CO 模型中的相關(guān)參數(shù)，進一步估算2014年玉米的rc和ET，篩選出最優(yōu)模型和參數(shù)優(yōu)化方法的組合，為提高rc的計算模型在懷來地區(qū)的適用性和計算精度提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

中國科學(xué)院懷來遙感綜合試驗站位于河北省懷來縣東花園鎮(zhèn)（115.7880°E，40.3491°N，海拔480 m）。站點所在區(qū)域處于半干旱地區(qū)，年均氣溫為10.1℃，最高氣溫39℃，最低氣溫達-20℃以下，年均降水總量約為370 mm，年均風(fēng)速為3.4 m/s，試驗站所在區(qū)域以沙質(zhì)沖積土為主，主要種植夏玉米，品種為鄭單958，每年5月上旬播種，9月中旬成熟[35-36]。

懷來站有渦度相關(guān)系統(tǒng)和自動氣象站各2 套和1套土壤水分無線傳感器網(wǎng)絡(luò)，可對地表通量和氣象要素進行實時觀測[37]。

1.2 數(shù)據(jù)處理

懷來站點2013年和2014年的氣象數(shù)據(jù)和渦度相關(guān)數(shù)據(jù)來自國家青藏高原科學(xué)數(shù)據(jù)中心提供的海河流域多尺度地表通量與氣象要素觀測數(shù)據(jù)集（http://data.tpdc.ac.cn/zh-hans/），內(nèi)容包括凈輻射Rn、土壤熱通量G、顯熱通量H、潛熱通量LE、風(fēng)速W、降水量P、土壤含水率θ、空氣濕度RH和氣溫T。儀器由于故障維修等因素，部分時間段出現(xiàn)數(shù)據(jù)缺失，日尺度數(shù)據(jù)由當(dāng)日有效數(shù)據(jù)取均值進行插補[38]，連續(xù)多日缺失的熱通量數(shù)據(jù)可以由非線性回歸方法進行插補[39]。

在裸露地面或者植被較低矮處，由熱力學(xué)第一定律可推導(dǎo)出地表能量平衡方程[40]：

式中：Rn為凈輻射（W/m2）；G為土壤熱通量（W/m2）；LE為潛熱通量（W/m2）；H為顯熱通量（W/m2）。式（1）左右二端相等時，代表能量閉合，反之則不閉合。能量不閉合的現(xiàn)象在測量系統(tǒng)中普遍存在，本文采用波文比能量強制閉合法使能量閉合[41]。

式中：D為能量不閉合差（W/m2）；β為波文比；ΔLE和ΔH為修正值，分別將其添加到對應(yīng)項即可得到LE和H修正后的數(shù)值。分析能量閉合度能夠有效提高測量精度[42-43]。

本文在進行參數(shù)擬合時，在觀測數(shù)據(jù)的基礎(chǔ)上剔除了降水量大于0 的觀測日，以提高后續(xù)參數(shù)優(yōu)化精度。

1.3 計算方法

1.3.1 Penman-Monteith 模型

1965年，Monteith 從能量平衡理論和大氣擴散理論出發(fā)，在 Penman 等的研究基礎(chǔ)上得到了Penman-Monteith（P-M）模型[7,38]，計算式為：

式中：λ為汽化潛熱（MJ/kg）；ET為實際蒸散量（mm）；Δ代表飽和水汽壓與溫度關(guān)系曲線的斜率；Rn為凈輻射（W/m2）；G為土壤熱通量（W/m2）；ρ為空氣密度（kg/m3）；Cp為空氣定壓比熱（J/（kg·℃））；VPD為飽和水汽壓差（kPa）；ra為空氣動力學(xué)阻力（s/m）；rc為冠層阻力（s/m）；γ為干濕表常數(shù)（kPa/℃）。

空氣動力學(xué)阻力ra可以通過計算[40]為：

式中：zm、zh分別為測風(fēng)速和測濕度的高度（zm=2.5m，zh=5m）；d為0 平面位移高度（m）；hc為作物平均高度（m）；zom為控制動量傳遞的粗糙度長度（m）；zoh為控制熱通量和水汽傳輸?shù)拇植陂L度。K為卡曼常數(shù)，取0.41，uz為參考高度的風(fēng)速（m/s）。

1.3.2 Jarvis 冠層阻力模型

Jarvis 模型（JA 模型）中包含了環(huán)境因素的響應(yīng)函數(shù)，考慮了輻射、溫度、濕度和土壤含水率對rc的影響。模型中各函數(shù)相互獨立，該模型是一個多變量的綜合性模型[24,44]：

式中：θF為田間持水率（m3/m3）；θW為凋萎含水率（m3/m3）。a1、a2、a3是經(jīng)驗參數(shù)，3 個經(jīng)驗參數(shù)和rcmin通過實測數(shù)據(jù)擬合得出。各參數(shù)擬合初值取a1=100[45]，a2=-0.12[24]，a3=0.0016[45]，rcmin=30 s/m[45]。

1.3.3 耦合冠層阻力模型

耦合冠層阻力模型（CO 模型）同時考慮了植物蒸騰和土壤蒸發(fā)的過程，所求rc指水汽通過冠層和土壤時克服的平均阻力，計算式為[28-31]：

1.3.4 利用Penman-Monteith 模型反算冠層阻力

利用站點觀測的ET，通過P-M 模型反算冠層阻力：

式中：λETEC為修正后的潛熱通量（W/m2），其余參數(shù)含義同式（1）。

1.4 敏感度分析

采用BP 神經(jīng)網(wǎng)絡(luò)模型分析rc與各影響因子的敏感度。BP 神經(jīng)網(wǎng)絡(luò)模型是一種按照誤差反向傳播的多層神經(jīng)網(wǎng)絡(luò)模型。本文中，模型輸入值為各影響因子（Rn、T、LAI、VPD和θ），輸出值為rc。使用BP神經(jīng)網(wǎng)絡(luò)模型控制輸入值，排除某一影響因子，進而改變模型的輸出值，通過輸出值的變化比較rc對各影響因子的敏感程度。排除某一因子后輸出值變化越大，則rc對被排除的因子敏感程度越高，反之越低[46-47]。

1.5 參數(shù)優(yōu)化算法

1.5.1 最小二乘法

最小二乘法（LSM）是一種進行曲線擬合的方法，在處理大量數(shù)據(jù)和分析誤差中具有廣泛的應(yīng)用。其基本原理[48]為：求一擬合目標函數(shù)φ(x)，使其在各點的偏差平方和達到最小：

式中：δ為擬合函數(shù)與原始數(shù)據(jù)差值；φ(xi)為擬合目標函數(shù)；yi為原始數(shù)據(jù)值。

1.5.2 蟻群算法

蟻群算法（ACO）通過模擬螞蟻覓食過程中個體間依靠信息素進行交流的過程，進而得到從蟻群到食物源的最短路徑的過程來進行優(yōu)化[49-50]。采用蟻群算法求優(yōu)化解的流程為：①預(yù)估模型參數(shù)范圍，初始化蟻群A(k)。預(yù)設(shè)蟻群規(guī)模、信息素蒸發(fā)系數(shù)、信息素強度系數(shù)和螞蟻爬行速度參數(shù)。②為每一只螞蟻選擇下一個節(jié)點。遍歷所有節(jié)點后計算最佳路徑。③更新信息素矩陣，目標至越優(yōu)的方格信息素越強。④檢查終止條件，看是否達到最大迭代次數(shù)或者最優(yōu)螞蟻是否達到要求精度。滿足條件則進行下一步，不滿足則返回步驟②，重復(fù)執(zhí)行步驟②、③、④。⑤輸出最優(yōu)值。

1.6 模型精度評價參數(shù)

本文采用決定系數(shù)R2，均方根誤差RMSE和一致性指數(shù)d來評價模型擬合的精度。

式中：Emi和EOi分別為模擬值和實測值；是實測值的平均值；n為樣本容量。本文中，利用λETEC帶入P-M 模型中反算所得的為實測值，利用模型擬合出的rc為模擬值；用相關(guān)儀器測得ETEC的為實測值，利用模型得到的ET為模擬值。

2 結(jié)果與分析

2.1 氣象數(shù)據(jù)的逐日變化

2013、2014年夏玉米全生育期的凈輻射Rn，風(fēng)速W，溫度T，飽和水汽壓差VPD，降水量P和土壤含水率θ的逐日變化如圖1所示。2013年和2014年玉米全生育期的Rn基本在50～250 W/m2間變化，2013年Rn的最大值為245.06 W/m2，最小值為23.04 W/m2，2014年Rn的最大值為218.95 W/m2，最小值為25.04 W/m2，6月末至7月初的Rn值在一個較高水平，之后隨時間呈下降趨勢；風(fēng)速W的值在2014年波動較大，2013年的最大值出現(xiàn)在5月19日，為5.28 m/s，2014年的最大值出現(xiàn)在5月25日，為5.97 m/s，隨后W呈下降趨勢；在夏玉米的生育期內(nèi)，溫度T的變化范圍主要在15～30℃之間，趨勢大致呈拋物線形式，在7、8月達到峰值，2013年的峰值為27.36℃，2014年的峰值為28.96℃；飽和水汽壓差VPD整體隨時間呈下降趨勢，2013年的最大值和最小值分別為2.31kPa 和0.29 kPa，在2014年波動較大，最大值和最小值分別為3.30 kPa 和0.15 kPa；2013年生育期內(nèi)總降水量P為374.5 mm，2014年生育期內(nèi)總P為293.4mm。P在時間上有明顯的分布不均現(xiàn)象，主要集中于6月下旬和7月下旬，有時會出現(xiàn)較大的日降水量。2013年的日P峰值出現(xiàn)在7月15日，為34.8 mm，2014年的日降水量峰值出現(xiàn)在8月28日，為24.6 mm；土壤含水率θ的變化趨勢明顯受到T和P的影響。每年7月T的水平較高，θ的值保持在一個較低的水平，在有降水時變化不明顯。進入8月后T開始下降，此時當(dāng)有降水后，θ的值會有明顯上升。

圖1 2013年和2014年玉米生育期內(nèi)氣象數(shù)據(jù)的逐日變化Fig.1 Daily variation of meteorological factors during growing season of maize in 2013 and 2014

2.2 實際冠層阻力的逐日變化及敏感程度分析

2013年、2014年的夏玉米全生育期的實際rc的逐日變化如圖2所示。實際rc在夏玉米生育前期時最高，接近5000 s/m，之后明顯減少，到6月下旬低于1000 s/m。結(jié)合圖1 與圖2，各氣象因子和LAI的變化會導(dǎo)致rc的值變化，如6、7月的T、Rn在一個較高水平，此時rc值下降較快。每年6月下旬有較為集中的降水，此時rc出現(xiàn)明顯的上升。

圖2 實際冠層阻力的逐日變化Fig.2 Daily variation of measured canopy resistance

利用BP 神經(jīng)網(wǎng)絡(luò)模型分析rc對各影響因子的敏感程度，結(jié)果如表1所示。在保留全部影響因子的條件下，輸出值與實測值間的相關(guān)性較高，R2=0.92，RMSE=346.61 s/m，d=0.90。當(dāng)控制Rn輸入后，R2=0.23，RMSE=1108.20 s/m，d=0.62，說明輸出值和實際值的相關(guān)關(guān)系明顯下降，rc對Rn的敏感程度較高；當(dāng)控制T輸入后，R2=0.92，RMSE=361.33 s/m，d=0.89；控制VPD輸入后，R2=0.92，RMSE=350.40 s/m，d=0.90；控制θ輸入后，R2=0.91，RMSE=374.41 s/m，d=0.89?？梢妑c對T、VPD和θ的敏感性程度都較低。當(dāng)控制LAI輸入后，R2=0.42，RMSE=959.76 s/m，d=0.63，rc對LAI的變化較為敏感。綜合可知rc對這5 個因子的敏感程度為：Rn＞LAI＞θ＞T＞VPD。

表1 敏感程度分析結(jié)果Table 1 Results of sensitivity analysis

2.3 模型參數(shù)擬合結(jié)果

以2013年實測的渦度相關(guān)數(shù)據(jù)通過P-M 模型反算得到的實際rc為標準，分別采用最小二乘法和蟻群算法對JA 模型和CO 模型中的參數(shù)進行擬合，擬合結(jié)果見表2。2 種優(yōu)化算法對參數(shù)擬合的效果較好，模型估計值和實測值的變化趨勢相同，R2的值均超過0.8?？傮w上CO 模型比JA 模型的擬合效果更優(yōu)，其中使用蟻群算法優(yōu)化的 CO 模型擬合結(jié)果最好（R2=0.89，RMSE=410.90 s/m，d=0.88），參數(shù)擬合結(jié)果為b1=2.20，b2=0.18，b3=0.01，b4=-0.18，b5=-1.39，b6=-0.45。用蟻群算法優(yōu)化的CO 模型在6月中旬以前的估計值比較接近實測值，之后的估計值和實測值之間有一定差距。用蟻群算法優(yōu)化的JA 模型估計rc的效果略好（R2=0.83，RMSE=515.52 s/m，d=0.85），參數(shù)擬合值分別為rcmin=7.43 s/m，a1=360.80，a2=0.25，a3=-0.003。

將4 種方法所得參數(shù)分別帶入對應(yīng)模型，利用2014年的氣象數(shù)據(jù)對夏玉米全生育期ET進行估計，并用渦度相關(guān)系統(tǒng)實測的ET值進行驗證，結(jié)果見表3。在JA 模型中，采用最小二乘法得到的R2=0.65，RMSE=1.20 mm，d=0.76，采用蟻群算法得到的R2=0.67，RMSE=1.16 mm，d=0.76，二者對JA 模型的擬合效果均較好；在CO 模型中，采用最小二乘法得到的R2=0.70，RMSE=1.11 mm，d=0.75，采用蟻群算法得到的R2=0.72，RMSE=1.07 mm，d=0.75，因此在計算ET時，誤差最小的是使用蟻群算法優(yōu)化后的CO 模型，總體上CO 模型的計算結(jié)果要優(yōu)于JA 模型，且在2 種模型中采用蟻群算法得到的結(jié)果精度更高。

3 討論

研究表明，rc在夏玉米生育前期達到最大值，之后隨著夏玉米生長逐漸下降，由于夏玉米生育期是5—9月，進入5月后，溫度開始上升，日照時間延長，加速了水分的蒸發(fā)，同時玉米葉面積指數(shù)變大，氣候變化和夏玉米生長使氣孔開度變大，最終導(dǎo)致了rc值變小。rc對Rn的敏感程度最高，因為Rn是地表能量的重要組成部分，在地表能量平衡和水分平衡中有直接參與[50]；LAI對rc的影響來自于作物本身的生理參數(shù)變化；其他氣候因素（θ，T和VPD）的變化對rc有影響，但是影響程度遠不及Rn和LAI明顯。這與Li 等[46]、Ding 等[52]和Zhang 等[53]研究結(jié)果一致。

優(yōu)化后CO 模型在估算rc和ET時的精度優(yōu)于JA模型，JA 模型對rc估算值偏小。CO 模型中不僅考慮了水分通過植物冠層時克服的阻力，還考慮了水分從土壤表面蒸發(fā)時克服的阻力[51]，無論對于稀疏還是密閉的植被都具有一定應(yīng)用性。夏玉米生育前期，由于植被稀疏，土壤表面的水分蒸發(fā)量占夏玉米ET的絕大部分，因此忽略土壤阻力的作用將會帶來一定的誤差。JA 模型沒有考慮土壤表面的阻力，因此估計的rc值偏小，可能導(dǎo)致估算的ET值偏大[29-30]。2 種優(yōu)化方法中，蟻群算法的魯棒性好，全局搜索能力強，能夠收斂于全局最優(yōu)解，在求解過程中能夠不依賴人工調(diào)整，適用于優(yōu)化求解問題[54]。而最小二乘法是線性估計，當(dāng)變量間為明顯的非線性關(guān)系時，最小二乘法得到的結(jié)果誤差較大，因此最小二乘法在使用上局限性較大。不論是對參數(shù)進行優(yōu)化，還是對2014年夏玉米的ET進行估算，蟻群算法的精度均優(yōu)于最小二乘法，可成為參數(shù)優(yōu)化的推薦算法。

由于數(shù)據(jù)質(zhì)量的影響，本文僅選取了夏玉米2a生長期渦度相關(guān)數(shù)據(jù)和氣象數(shù)據(jù)。在以后的研究中應(yīng)選擇時間跨度更長的數(shù)據(jù)，以提高夏玉米rc的參數(shù)擬合精度，進一步提高ET估計值的可靠度，為懷來地區(qū)的夏玉米作物需水量估計和農(nóng)業(yè)水管理提供更為合理的科學(xué)參考。

4 結(jié)論

1）rc在玉米的生育期呈下降趨勢，在生育前期rc值較大，接近5000 s/m。rc與各氣象因子和LAI存在明顯的非線性關(guān)系，采用BP 神經(jīng)網(wǎng)絡(luò)模型進行敏感分析，得到rc對各影響因子敏感程度為：Rn＞LAI＞θ＞T＞VPD。

2）最小二乘法和蟻群算法對rc的計算模型中參數(shù)擬合的效果較好，模型估計值和實測值的變化趨勢相同?？傮w上CO 模型比JA 模型的擬合效果更優(yōu)，其中使用蟻群算法優(yōu)化的CO 模型擬合結(jié)果最好，R2=0.89，RMSE=410.90 s/m，d=0.88，且6月中旬以前的估計值比較接近實測值，之后的估計值和實測值之間有一定差距。

3）利用2014年氣象數(shù)據(jù)，采用優(yōu)化后的rc的計算模型計算2014年的玉米ET，CO 模型的計算ET值精度高于JA 模型，其中采用蟻群算法優(yōu)化后的CO模型計算ET精度最高，R2=0.70，RMSE=31.44 mm，d=0.75，采用最小二乘法優(yōu)化后的JA 模型計算ET的精度較低，R2=0.70，RMSE=31.44 mm，d=0.76。