非參數(shù)化蒸散發(fā)估算方法在華北灌溉農(nóng)田的適用性評(píng)價(jià)*

張曉龍,張玉翠,石嘉麗,2,王 妍,2,3,沈彥軍,2**

(1.中國(guó)科學(xué)院遺傳與發(fā)育生物學(xué)研究所農(nóng)業(yè)資源研究中心/中國(guó)科學(xué)院農(nóng)業(yè)水資源重點(diǎn)實(shí)驗(yàn)室/河北省節(jié)水農(nóng)業(yè)重點(diǎn)實(shí)驗(yàn)室 石家莊 050022; 2.中國(guó)科學(xué)院大學(xué) 北京 100049; 3.河北師范大學(xué)資源與環(huán)境科學(xué)學(xué)院/河北省環(huán)境演變與生態(tài)建設(shè)實(shí)驗(yàn)室 石家莊 050024)

蒸散發(fā)是水分從地表進(jìn)入大氣的過程,是土壤蒸發(fā)、冠層截留和植物蒸騰的總和。蒸散發(fā)是能量、水分和碳循環(huán)過程中的重要紐帶,也是區(qū)域水資源重要的消耗途徑。在全球范圍內(nèi)超過60%的降水通過蒸散發(fā)返回大氣,而在干旱區(qū)、農(nóng)田這一比例甚至超過90%。蒸散發(fā)以潛熱通量(latent heat flux,LE)的形式消耗了大約50%的地表吸收的太陽(yáng)輻射量。因此,精準(zhǔn)估算蒸散發(fā)已成為水量平衡計(jì)算、農(nóng)業(yè)水資源管理、水文模擬等眾多領(lǐng)域共同關(guān)注的核心問題。

盡管蒸散發(fā)可通過渦度相關(guān)系統(tǒng)、蒸滲儀、Bowen 比、大孔徑閃爍儀等方法直接觀測(cè),但由于這些原位技術(shù)的局限性,估算大尺度、長(zhǎng)期蒸散發(fā)仍然是一項(xiàng)重大的挑戰(zhàn)。因此,許多直接或間接的估算方法被提出,包括Penman 公式、Penman-Monteith 公式、經(jīng)驗(yàn)統(tǒng)計(jì)方法、地表能量平衡法、SPAC 系統(tǒng)模擬方法、互補(bǔ)相關(guān)理論方法、陸面過程與數(shù)據(jù)同化方法等。然而,大多數(shù)估算方法都比較復(fù)雜或需要對(duì)地表性質(zhì)有詳細(xì)了解,從而限制了這些方法的應(yīng)用,降低了估算精度。例如,Penman-Monteith 方程需要輸入冠層阻力,Shuttleworth-Wallace 方法需要空氣動(dòng)力特性、冠層粗糙度和土壤阻力,但這些參數(shù)在區(qū)域尺度上很難測(cè)量和校準(zhǔn)。近年來,Liu 等提出了一種基于哈密頓廣義能量原理的非參數(shù)化地表蒸散發(fā)估算(the nonparametric approach,NP)方法。NP 方法只需要地表凈輻射、空氣溫度、地表溫度和土壤熱通量作為輸入,克服了傳統(tǒng)蒸散發(fā)估算方法中存在經(jīng)驗(yàn)參數(shù)化的物理缺陷,降低了計(jì)算過程的不確定性。該方法所有必要的輸入均為可測(cè)量的,為實(shí)際應(yīng)用提供了一個(gè)新穎而簡(jiǎn)單的方法。目前,該方法已在部分站點(diǎn)進(jìn)行了驗(yàn)證,并結(jié)合遙感數(shù)據(jù)和氣象驅(qū)動(dòng)數(shù)據(jù)進(jìn)行了區(qū)域蒸散發(fā)適用性檢驗(yàn)。Liu 等利用全球26 個(gè)渦度相關(guān)站點(diǎn)觀測(cè)結(jié)果對(duì)NP 方法進(jìn)行了評(píng)價(jià),其中23 個(gè)點(diǎn)性能表現(xiàn)良好; 并進(jìn)一步表明NP 方法與Penman-Monteith 方法在性能上是兼容的。Yang 等發(fā)現(xiàn)NP 方法在干燥條件下會(huì)高估蒸散發(fā),應(yīng)避免在干燥條件下應(yīng)用,而在濕潤(rùn)條件下表現(xiàn)良好。王寧等對(duì)比發(fā)現(xiàn)在濕潤(rùn)下墊面NP 方法會(huì)低估蒸散發(fā),在干旱下墊面會(huì)高估蒸散發(fā),且夏季估算精度優(yōu)于冬季; 并利用遙感數(shù)據(jù)估算了黑河上中游區(qū)域蒸散發(fā),結(jié)果表明雖然存在一定誤差,但能夠反映區(qū)域蒸散發(fā)分布特征。Pan 等在黑河流域6 個(gè)站點(diǎn)驗(yàn)證了NP 方法的有效性,并基于該方法和遙感數(shù)據(jù)估算了黑河流域蒸散發(fā),結(jié)果表明遙感反演的蒸散發(fā)時(shí)空分布基本可靠,但在非植被站點(diǎn),估算值與地面觀測(cè)值一致性較差,而相對(duì)誤差較小?？偟膩碚f,NP 方法有必要在不同區(qū)域或生態(tài)系統(tǒng)進(jìn)行更細(xì)致的適用性分析,以識(shí)別出有助于提高估算精度的誤差源。

華北平原是我國(guó)主要的灌溉農(nóng)業(yè)區(qū)和重要的糧食生產(chǎn)基地,對(duì)保障國(guó)家糧食安全具有舉足輕重的作用。該地區(qū)種植結(jié)構(gòu)以冬小麥()和夏玉米()為主,其種植面積約占糧食總播種面積的80%左右。該地區(qū)降水量少且集中于夏季,無法滿足作物生長(zhǎng)需求,灌溉成為維持和提高糧食產(chǎn)量不可缺少的手段。農(nóng)田灌溉改變了區(qū)域水熱條件與氣候,直接和間接地影響了蒸散發(fā)過程。精準(zhǔn)地估算灌溉農(nóng)田蒸散發(fā)不僅對(duì)華北平原干旱監(jiān)測(cè)和水資源管理具有重要理論意義,而且對(duì)區(qū)域種植結(jié)構(gòu)調(diào)整和生態(tài)環(huán)境建設(shè)具有重要實(shí)踐意義。目前關(guān)于NP 方法應(yīng)用研究大多集中于干旱區(qū)流域,關(guān)于半濕潤(rùn)區(qū)灌溉農(nóng)田的適用性分析報(bào)道較少。因此,本文主要利用華北灌溉農(nóng)田區(qū)的欒城農(nóng)業(yè)生態(tài)系統(tǒng)試驗(yàn)站(簡(jiǎn)稱欒城站)、禹城綜合試驗(yàn)站(簡(jiǎn)稱禹城站)和館陶試驗(yàn)站(簡(jiǎn)稱館陶站)3 個(gè)通量站點(diǎn)的觀測(cè)數(shù)據(jù),應(yīng)用NP 方法估算這3 個(gè)站點(diǎn)30 min 和日尺度蒸散發(fā)(以LE 的形式),并以渦度相關(guān)系統(tǒng)修正的LE 為參考,分析NP 方法在華北平原灌溉農(nóng)田不同季節(jié)和不同時(shí)間尺度的適用性,該研究不僅為NP 方法改進(jìn)提供參考,而且也有助于加深對(duì)蒸散發(fā)理論的認(rèn)識(shí)。

1 研究站點(diǎn)及數(shù)據(jù)處理

1.1 研究區(qū)及觀測(cè)站點(diǎn)概況

華北平原位于太行山以東,黃河以北,燕山以南,渤海以西,面積約14.1 萬km。該地區(qū)土層較深,海拔低于100 m,跨越北京、天津、河北、河南、山東等多個(gè)地區(qū)。該地區(qū)屬于暖溫帶半濕潤(rùn)半干旱大陸性季風(fēng)氣候,年平均氣溫10.0～14.2 ℃,光照充足。主要糧食作物為小麥、玉米,主要經(jīng)濟(jì)作物為棉花(spp.)、果樹和蔬菜(圖1)。年平均降水量為400～500 mm,70%發(fā)生在夏季,農(nóng)業(yè)灌溉消耗了地下水總量的80%以上。本研究以華北灌溉農(nóng)田為研究對(duì)象,選取欒城站(37°53′N,114°41′E,海拔50 m)、禹 城 站(36°49′N,116°34′E,海拔22 m)和館陶站(36°30′N,115°07′E,海拔30 m) 3 個(gè)通量觀測(cè)站點(diǎn),欒城站和禹城站下墊面均為充分灌溉管理下的典型冬小麥-夏玉米一年兩熟輪作制農(nóng)田生態(tài)系統(tǒng),而館陶站下墊面為冬小麥/夏玉米和棉花綜合下墊面。各通量站下墊面均勻且地勢(shì)平坦,滿足通量觀測(cè)源區(qū)要求,觀測(cè)樣地具有較強(qiáng)的代表性。一般情況下,冬小麥10月初播種,次年6月中旬收獲; 夏玉米6月中旬播種,9月底收獲; 棉花4月底播種,10月底收獲。根據(jù)農(nóng)田生態(tài)系統(tǒng)水分、土壤、氣象和生物監(jiān)測(cè)規(guī)范和野外田間試驗(yàn)要求,各站配備有國(guó)際先進(jìn)的儀器設(shè)備。渦度相關(guān)系統(tǒng)主要由紅外開路式CO/HO氣體分析儀和三維超聲風(fēng)速儀組成,同時(shí)輔助有輻射分量、常規(guī)氣象要素、土壤熱通量的同步觀測(cè),為該研究提供所需的觀測(cè)數(shù)據(jù)。

圖1 華北平原作物類型空間分布及通量站點(diǎn)位置Fig.1 Spatial distribution of agricultural land-use types and location of flux stations in the North China Plain

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

本研究使用的通量數(shù)據(jù)是各站點(diǎn)提供的經(jīng)過質(zhì)量控制的二級(jí)產(chǎn)品,包括30 min 平均LE 和感熱通量(sensible heat flux,)。由于儀器故障、不利的微氣象條件、異常值剔除和質(zhì)量控制等原因造成數(shù)據(jù)缺失,通常一年中約有17%～50%的觀測(cè)數(shù)據(jù)缺測(cè)或被剔除。在水資源管理或水文模型研究中,通常需要得到日、月甚至年的累加值,本研究使用德國(guó)耶拿大學(xué)提供的R 軟件ReddyProc 包進(jìn)行數(shù)據(jù)插補(bǔ),并進(jìn)一步統(tǒng)計(jì)為日平均值。

渦度相關(guān)技術(shù)目前已成為獲得生態(tài)系統(tǒng)LE 和的標(biāo)準(zhǔn)觀測(cè)手段之一,但是該方法普遍存在能量不閉合的問題?，F(xiàn)有能量閉合校正方法主要有兩種類型,即Bowen 比閉合修正(BR)方法和能量殘差閉合修正(ER)方法。雖然兩種閉合修正方法都可能對(duì)系統(tǒng)內(nèi)的能量和水平衡檢測(cè)中有重大或潛在的不利影響,但對(duì)于大多數(shù)地表?xiàng)l件,采用ER 修正后的潛熱通量與各估算方法模擬結(jié)果的一致性更高。因此,采用ER 方法修正的LE 作為評(píng)價(jià)本研究蒸散發(fā)估算方法的基礎(chǔ)。ER 方法表達(dá)公式如下:

式中:LE為能量殘差法修正后LE (W?m),為渦動(dòng)相關(guān)系統(tǒng)測(cè)量的(W?m),為地表凈輻射(W?m),為土壤熱通量(W?m)。

為保障數(shù)據(jù)一致性,將輻射分量、空氣溫度、地表溫度、土壤熱通量等10 min 頻率觀測(cè)數(shù)據(jù)統(tǒng)計(jì)為30 min 和日平均值; 而降水統(tǒng)計(jì)為累加值。

2 研究方法

2.1 非參數(shù)化蒸散發(fā)估算方法

在NP 方法中,在宏觀層面上假設(shè)下墊面為均勻的,哈密頓量(勢(shì)能與動(dòng)能的總和)為該系統(tǒng)的總能量。作為勢(shì)能,則、和LE 為動(dòng)能; 地表溫度()作為這個(gè)系統(tǒng)的廣義坐標(biāo)。該方法計(jì)算了哈密頓量關(guān)于的偏微分方程,詳細(xì)推導(dǎo)過程見文獻(xiàn)[21]。計(jì)算公式如下:

式中:LE為NP 估算得到的LE (W·m);為地表空氣溫度(K); ?為地表比輻射率,由于數(shù)據(jù)的限制,本研究假設(shè) ?為常數(shù),在農(nóng)田下墊面站取0.95; σ為Stefan-Boltzman 常數(shù)(5.67× 10W?m?K); γ為干濕表常數(shù)(kPa?℃),可用近地表壓力來估算; ?為時(shí)的飽和水汽壓梯度(kPa?℃),其計(jì)算公式如下:

2.2 精度評(píng)價(jià)方法

本研究選取平均偏差(Bias)、相對(duì)誤差(RE)、決定系數(shù)(R)、均方根誤差(RMSE)和Nash-Sutcliffe效率系數(shù)(NSE)來量化估算方法模擬精度。其計(jì)算公式如下:

2.3 敏感性分析方法

為了檢測(cè)模型各驅(qū)動(dòng)變量對(duì)輸出結(jié)果的貢獻(xiàn),本研究采用簡(jiǎn)單的相對(duì)敏感系數(shù)法進(jìn)行敏感性分析。在這個(gè)方法中,通過比較改變輸入變量產(chǎn)生的LE(LE)相對(duì)于參考LE (LE)的比例來計(jì)算每個(gè)變量的敏感性。本研究中所有輸入變量的變化范圍設(shè)為±30%,取其中一種輸入?yún)?shù)逐步變化5%,而其他輸入?yún)?shù)保持不變。敏感系數(shù)(S)計(jì)算公式如下:

3 結(jié)果與分析

3.1 通量站點(diǎn)能量閉合分析

本文用線性回歸方法和能量平衡比法分析了欒城、禹城和館陶3 個(gè)站點(diǎn)的能量平衡狀況(圖2)。其中,線性回歸方法中斜率代表能量平衡程度,理想狀況下斜率為1,截距為0,但實(shí)際情況下截距通常不為0; 能量平衡比(EBR)即較長(zhǎng)時(shí)間段(如1年)熱通量之和(LE+)與可利用能量(?)的比值,通常EBR 范圍為0.34～1.69。從圖中可以看出,在日尺度上,欒城、禹城和館陶3 個(gè)站點(diǎn)的能量閉合率分別為83.7%、87.4%和51.7%,EBR 分別為0.95、1.06和0.94; 在30 min 尺度上這3 個(gè)站點(diǎn)的能量閉合率分別為81.1%、84.5%和73.7%,EBR 分別為1.13、1.06 和0.95?？芍?3 個(gè)站點(diǎn)的EBR 接近于1,說明這3 個(gè)站點(diǎn)數(shù)據(jù)公開前已經(jīng)過數(shù)據(jù)質(zhì)量控制; 欒城站和禹城站能量閉合率更好一些,且日尺度和30 min 尺度能量閉合率相差不大,表明這兩個(gè)臺(tái)站數(shù)據(jù)質(zhì)量相對(duì)更好一些。需說明的是,理論上每個(gè)站點(diǎn)的EBR 在不同時(shí)間尺度上應(yīng)是基本一致的,但本研究中欒城站2月份因缺少部分?jǐn)?shù)據(jù),故欒城站日尺度和30 min 尺度EBR 存在較大差異。一般認(rèn)為,如果能量閉合觀測(cè)誤差超過10%,其觀測(cè)結(jié)果對(duì)大氣模式特別是陸面過程模式的檢驗(yàn)是不能接受的。本研究發(fā)現(xiàn),即便這3 個(gè)站點(diǎn)經(jīng)過嚴(yán)格的質(zhì)量控制,依舊有超過12.6%以上的能量丟失,這樣會(huì)造成LE 低估現(xiàn)象。另外,隨著蒸散發(fā)估算理論和遙感技術(shù)的發(fā)展,需要提供較可靠的地面驗(yàn)證數(shù)據(jù),所以對(duì)通量觀測(cè)數(shù)據(jù)的校正十分有必要。

圖2 日尺度和30 min 尺度上3 個(gè)站點(diǎn)通量站點(diǎn)的能量平衡特征Fig.2 Energy balance characteristics of flux stations at daily scale and 30 min scale in the three stations

3.2 灌溉農(nóng)田蒸散發(fā)估算

基于NP 方法,利用通量站觀測(cè)的、、和等數(shù)據(jù),估算了2010年的欒城、禹城和館陶3 個(gè)華北平原典型灌溉農(nóng)田的日尺度和30 min 尺度的LE,將估算值LE與修正后的通量觀測(cè)值LE進(jìn)行對(duì)比驗(yàn)證分析。需說明的是,在30 min 尺度上,欒城站在2月缺失部分?jǐn)?shù)據(jù),剔除不用。表1 顯示了2010年3 個(gè)通量站點(diǎn)觀測(cè)的能量通量和環(huán)境參數(shù)均值。從表中可知,3 個(gè)站點(diǎn)的各輸入?yún)?shù)值相差不大,LE從大到小依次為欒城站、禹城站和館陶站;而降水量與LE呈反比,這可能是由于灌溉活動(dòng)保證了下墊面較穩(wěn)定的水分條件,在一定灌溉保證率的農(nóng)田中降水并不是蒸散發(fā)的控制因素; 欒城站、禹城站和館陶站潛熱比(LE/)分別為81.6%、82.3%和72.3%。欒城站和禹城站的潛熱比很接近,這可能與兩站下墊面條件更接近有關(guān)。農(nóng)田蒸散發(fā)在很大程度上受可利用能量控制,欒城站和禹城站潛熱比類似,欒城站可利用能量大于禹城站,故欒城站蒸散量大于禹城站; 另外,下墊面作物類型對(duì)農(nóng)田蒸散發(fā)也有較大影響,如禹城站和館陶站具有接近的可利用能量,但由于禹城站下墊面為冬小麥/夏玉米輪作農(nóng)田,館陶站則是冬小麥/夏玉米和棉花混合下墊面,導(dǎo)致館陶站蒸散量明顯低于禹城站。因?yàn)檫@3 個(gè)站點(diǎn)土壤類型、質(zhì)地差別不大,可忽略其影響。綜上,由于灌溉活動(dòng)保證了下墊面較穩(wěn)定的水分條件,所以華北平原灌溉農(nóng)田蒸散發(fā)的主要影響因素為可利用能量和作物類型。

表1 3 個(gè)通量站點(diǎn)能量通量和環(huán)境參數(shù)平均值Table 1 Average of energy fluxes and environment parameters at the three stations

3.3 非參數(shù)化方法估算結(jié)果精度評(píng)價(jià)

圖3 為2010年的欒城、禹城和館陶LE 逐日過程及與觀測(cè)值的對(duì)比結(jié)果?？傮w而言,欒城站和禹城站蒸散發(fā)在年內(nèi)明顯呈雙峰型,館陶站則呈不明顯的雙峰型,模擬效果可以較好地反映蒸散發(fā)年內(nèi)變化過程。從箱體圖結(jié)果可知,在灌溉農(nóng)田中NP 方法總體上低估LE,欒城、禹城和館陶站LE的均值較LE分別低12.88 W?m、16.18 W?m和15.02 W?m,平均低估14.69 W?m; 另外,LE年內(nèi)波動(dòng)范圍大于LE。在日尺度上綜合所有站點(diǎn)的LE與LE之 間 的Bias、RE、、RMSE 和NSE分別為?14.73 W?m、 27.5%、 0.81、 27.75 W?m和0.71(表2)。總之,在日尺度上NP 方法在3 個(gè)農(nóng)田站點(diǎn)總體表現(xiàn)較好,且欒城站模擬精度優(yōu)于館陶站和禹城站。

表2 日尺度和30 min 尺度3 個(gè)站點(diǎn)的蒸散發(fā)估算精度評(píng)價(jià)Table 2 Accuracy evaluation of evapotranspiration estimates at daily scale and 30 min scale in the three stations

圖3 日尺度上3 個(gè)站點(diǎn)的潛熱通量(LE)年內(nèi)變化過程及估算結(jié)果對(duì)比Fig.3 Comparison of latent heat flux (LE) between estimation and observation at daily scale in the three stations

因30 min 尺度數(shù)據(jù)量龐大,不便完整展示,本研究 以2010年9月1日0 點(diǎn) 至9月15日23 點(diǎn)30 分完整時(shí)間段為例,展示欒城、禹城和館陶3 站的LE日內(nèi)變化過程及與觀測(cè)值的對(duì)比結(jié)果(圖4)。該時(shí)間段夏玉米為乳熟期、棉花為吐絮期,作物生長(zhǎng)趨于成熟; 并且該時(shí)間段內(nèi)3 站有明顯的降水過程,降水量存在較大差異。從圖中可知,在示例時(shí)間段內(nèi)30 min 尺度上,LE與LE高度一致,估算結(jié)果能較理想地反映LE日內(nèi)變化過程,即使是陰雨天氣精度依舊很好。經(jīng)計(jì)算,3 個(gè)站點(diǎn)該時(shí)間段內(nèi)R可達(dá)0.97～0.99,NSE 可達(dá)0.96。在30 min 尺度上綜合2010年全年所有站點(diǎn)的LE與LE之間的Bias、RE、、RMSE 和NSE 分別為?10.34 W?m、21.5%、0.88、40.37 W?m和0.87(表2)。30 min 尺度上,非參數(shù)方法在欒城站模擬精度優(yōu)于館陶站和禹城站。相比于日尺度,非參數(shù)方法在3 個(gè)農(nóng)田站點(diǎn)總體表現(xiàn)除RMSE 指標(biāo)外均有所提升,表明該方法在30 min 這樣短時(shí)間尺度上表現(xiàn)更好?？傊?NP 方法在灌溉農(nóng)田低估LE,但是可以較好地捕捉LE 的年內(nèi)和日內(nèi)變化過程,具有較高的應(yīng)用價(jià)值。

圖4 30 min 尺度上3 個(gè)站點(diǎn)的潛熱通量日內(nèi)變化過程及估算結(jié)果對(duì)比(9月1—15日)Fig.4 Comparison of latent heat flux between estimation and observation at 30 min scale (from September 1st to 15th) in the three stations

3.4 非參數(shù)化方法在不同季節(jié)的適用性

從3 個(gè)站點(diǎn)不同季節(jié)的估算精度分析指標(biāo)(圖5)可以看出,在日尺度和30 min 尺度上NP 方法呈現(xiàn)季節(jié)變化的相似性,但在RE 和NSE 的量級(jí)上存在差異。Bias 整體上呈負(fù)值,最小值出現(xiàn)在4、5月份; RE 在冬季誤差較大,在作物生長(zhǎng)季表現(xiàn)穩(wěn)定且良好;年內(nèi)變化呈雙峰型,與年內(nèi)各月LE 變化過程一致,即當(dāng)LE 較高時(shí)可達(dá)0.8 以上,而在LE較低時(shí)(冬季和夏玉米苗期)普遍低于0.4; RMSE整體上呈單峰型,在日尺度上,RMSE 最大值出現(xiàn)在4、5月份,而在30 min 尺度上,RMSE 最大值出現(xiàn)在5、6月份; NSE 在各月的表現(xiàn)類似,作物生長(zhǎng)季表現(xiàn)優(yōu)于非生長(zhǎng)季,且30 min 尺度表現(xiàn)優(yōu)于日尺度。總體上,NP 方法在生長(zhǎng)季估算精度高于非生長(zhǎng)季,玉米估算精度高于小麥,夏季估算精度高于冬季。

圖5 日尺度和30 min 尺度3 個(gè)站點(diǎn)的蒸散發(fā)估算精度在各季節(jié)上的評(píng)價(jià)Fig.5 Accuracy evaluation of evapotranspiration estimates at daily scale and 30 min scale in the three stations

3.5 非參數(shù)化方法敏感性分析

為了分析各輸入?yún)?shù)對(duì)估算結(jié)果的影響,本研究分別在日尺度和30 min 尺度上對(duì)NP 方法各參數(shù)進(jìn)行敏感性分析。NP 方法輸入?yún)?shù)包括、、和,選取其中一種參數(shù)從變化±5%到±30%,而其余參數(shù)保持不變,得到各參數(shù)敏感性(圖6)。從圖中可知,和的敏感度與參數(shù)變化同向,和的敏感度與參數(shù)變化異向。日尺度上、、和參數(shù)變化5%時(shí)的敏感性系數(shù)在3 個(gè)站點(diǎn)平均分別為11.3%、?10.4%、6.0%和?0.24%; 30 min 尺度上這4 個(gè)參數(shù)變化5%時(shí)的敏感性系數(shù)分別為11.75%、?10.94%、6.24%和?0.38%。相對(duì)于站點(diǎn)來說,在日尺度上,參數(shù)、和敏感程度館陶站>禹城站>欒城站,而敏感程度禹城站>館陶站>欒城站; 在30 min 尺度上,參數(shù)和敏感程度館陶站>欒城站>禹城站,敏感程度館陶站>禹城站>欒城站;敏感程度禹城站>館陶站>欒城站,影響可忽略不計(jì),例如當(dāng)在敏感度最高的站點(diǎn)(禹城站)和時(shí)間尺度上(30 min)變化?30%時(shí),LE 僅變化2.66%?？傮w而言,無論在日尺度上還是在30min 尺度上,各參數(shù)的敏感性在各站點(diǎn)上沒有明顯差別,相較而言館陶站敏感性最高; 各參數(shù)敏感度依次為、、和。鑒于此,雖然是對(duì)通量站地表能量平衡的重要分量,但當(dāng)缺乏土壤熱通量觀測(cè)儀器埋深和土壤質(zhì)地等數(shù)據(jù)缺乏時(shí),可不考慮埋深(儀器以上土壤熱儲(chǔ)存量)對(duì)的影響。

圖6 非參數(shù)化方法在3 個(gè)站點(diǎn)的參數(shù)敏感性分析Fig.6 Sensitivity analysis of the parameters of the nonparametric approach in the three stations

4 討論

4.1 潛熱通量修正的影響

渦度相關(guān)方法目前已被廣泛應(yīng)用于陸地生態(tài)系統(tǒng)水碳循環(huán)過程研究,其測(cè)定的LE 結(jié)果常作為實(shí)際蒸散發(fā)“驗(yàn)證值”來完善或改進(jìn)蒸散發(fā)估算模型的參數(shù)化方案。然而,該方法的測(cè)量誤差可能對(duì)模型精度評(píng)價(jià)產(chǎn)生直接影響,其中30 min 尺度上LE 的標(biāo)準(zhǔn)誤差約為5%～20%,從30 min 到日的時(shí)間尺度上升會(huì)導(dǎo)致出現(xiàn)5%～10%的偏差。另外,基于渦度相關(guān)技術(shù)的能量平衡研究一直受到能量不閉合問題的挑戰(zhàn),該問題是近30年來困擾地-氣交換試驗(yàn)研究的主要難點(diǎn)之一。大量觀測(cè)結(jié)果表明,幾乎每個(gè)觀測(cè)站都存在地表能量平衡不閉合的問題,一般湍流通量(包括LE 和)僅占地表可利用能量的70%～90%。

Bowen 比閉合修正(BR)方法和能量殘差閉合修正(ER)方法是目前最常用的兩種方法。BR 方法是根據(jù)Bowen 比分配剩余能量到LE 和中,假設(shè)LE和具有相似的觀測(cè)精度,然而根據(jù)原位觀測(cè)數(shù)據(jù),渦度相關(guān)系統(tǒng)對(duì)的觀測(cè)精度要高于LE 的精度;此外,比LE 的隨機(jī)誤差的季節(jié)差異更小。因此,當(dāng)對(duì)和LE 測(cè)量的準(zhǔn)確性水平存在顯著差異時(shí),應(yīng)謹(jǐn)慎使用BR 方法。ER 方法是基于能量平衡原理,這種方法只依賴于、和的測(cè)量精度,但是這種修正方式將不閉合的能量全部歸于LE 值中,不可避免的會(huì)高估LE。當(dāng)空氣儲(chǔ)熱、冠層儲(chǔ)熱、地表到土壤熱通量板之間的土壤儲(chǔ)熱過高時(shí),ER 方法也會(huì)引起較大誤差。此外,使用Bias、RE、、RMSE和NSE 等指標(biāo)來表示模擬值與觀測(cè)值之間的誤差,通常假設(shè)將所有誤差都?xì)w結(jié)于模擬值中,即認(rèn)為觀測(cè)值無誤差。然而,這種假設(shè)在實(shí)際上幾乎不可能發(fā)生,觀測(cè)值的選取及校正方法的使用會(huì)直接影響精度評(píng)估結(jié)果,所以只有固定“真值”(觀測(cè)值)時(shí)對(duì)比精度評(píng)估結(jié)果才有意義。本研究以日尺度為例,分別將LE對(duì)比未修正的LE (LE)、BR 方法修正的LE (LE)和ER 方法修正的LE(LE)進(jìn)行驗(yàn)證,精度評(píng)價(jià)相關(guān)參數(shù)見表3。從表中可知,兩種修正方法使得LE與LE和LE一致性指標(biāo)(、RMSE 和NSE)均有所提升,偏差性指標(biāo)(Bias 和RE)略有下降;其中3 個(gè)站點(diǎn)平均的分別從0.57 提升到0.80 和0.82,RMSE 分別從33.41 W·m下降到26.80 W·m和27.50 W·m,NSE 分別從0.42 提升到0.69 和0.71,Bias 分別從?13.76 W·m下降到?14.86 W·m和?14.73 W·m,RE 分別從26.1%下降到28.1%和27.9%。相較而言,ER 修正后一致性優(yōu)于BR 修正,該結(jié)果與上述文獻(xiàn)結(jié)果一致。本研究目的是驗(yàn)證蒸散發(fā)估算方法的適用性,需提供一個(gè)較穩(wěn)定可靠的驗(yàn)證值,相較而言,ER 方法更適合灌溉農(nóng)田的LE 修正,故本文選擇利用ER 方法修正潛熱通量作為評(píng)價(jià)本研究蒸散發(fā)估算方法的基礎(chǔ)。但值得注意的是,也許因?yàn)長(zhǎng)E也是基于能量殘差進(jìn)行計(jì)算,LE和LE并非完全獨(dú)立,導(dǎo)致ER 方法修正結(jié)果表現(xiàn)更強(qiáng)的一致性,故最佳的能量閉合修正方法的選擇仍有待進(jìn)一步探索。

表3 日尺度上3 個(gè)站點(diǎn)的不同修正方法的潛熱通量對(duì)模擬結(jié)果的精度評(píng)價(jià)Table 3 Accuracy evaluation of simulation results by latent heat flux of different correction methods at daily scale in the three stations

4.2 非參數(shù)化方法適用性檢驗(yàn)

實(shí)際蒸散發(fā)精準(zhǔn)估算一直是水文、微氣象學(xué)研究領(lǐng)域的研究難點(diǎn)之一。NP 方法提供了一種新穎而簡(jiǎn)單的估算方法。雖然NP 方法基于簡(jiǎn)化的空氣動(dòng)力阻力框架,并結(jié)合平衡蒸發(fā)和的變化來估算LE,但該方法因減少了輸入變量的累積誤差,反而具有其他估算模型(Penman-Monteith 公式)相當(dāng)?shù)木? 實(shí)際上,由于不同作物的空氣動(dòng)力學(xué)阻力不同,采用通用穩(wěn)定函數(shù)計(jì)算的空氣動(dòng)力學(xué)阻力可能會(huì)產(chǎn)生較大的誤差。由于觀測(cè)數(shù)據(jù)限制,本研究未檢驗(yàn)NP 方法在月、年這樣的長(zhǎng)時(shí)間尺度的表現(xiàn),有待將來進(jìn)一步研究。本研究基于華北3 個(gè)灌溉農(nóng)田站點(diǎn)在日尺度和30 min 尺度上進(jìn)行了NP 方法適用性檢驗(yàn),結(jié)果表明,NP 方法可以較好地反映季節(jié)及日內(nèi)變化特征,在日尺度和30 min 尺度上總體表現(xiàn)較好,且在不同站點(diǎn)沒有產(chǎn)生顯著的預(yù)測(cè)誤差。另外,NP方法在玉米生長(zhǎng)季估算精度高于小麥生長(zhǎng)季,主要原因可能是NP 方法對(duì)灌溉活動(dòng)反應(yīng)不敏感,并不能充分反映灌溉對(duì)農(nóng)田蒸散發(fā)的影響。玉米生長(zhǎng)季為7月至10月,華北平均降水多集中于該時(shí)期,基本可保證作物生長(zhǎng)需求,而小麥生長(zhǎng)季為11月至次年6月,期間降水較少,無法滿足作物生長(zhǎng)需求,特別是在返青?拔節(jié)期,需要進(jìn)行大量灌溉; 因NP 方法對(duì)灌溉事件響應(yīng)能力不足,使得小麥季低估較玉米季更嚴(yán)重,且模擬精度更差,總體而言,NP 方法在華北灌溉農(nóng)田不同作物類型、不同時(shí)間尺度與參考值具有較高的一致性,但在水分供應(yīng)較充足時(shí)較嚴(yán)重低估LE,而在較干燥時(shí)輕度低估或不低估LE; 另外該方法對(duì)灌溉活動(dòng)的響應(yīng)考慮不足,有待進(jìn)一步改進(jìn)以提高灌溉農(nóng)田模擬精度。本研究結(jié)果與已有文獻(xiàn)結(jié)論基本一致。

研究發(fā)現(xiàn),NP 方法不僅可以較好地反映LE 的日內(nèi)變化、季節(jié)變化,而且基本不受陰雨天氣影響(圖4),這不僅表明NP 方法的穩(wěn)健性,更表明因不受云層的影響而展示出在遙感應(yīng)用中巨大潛力; 這與Liu 等研究結(jié)果一致。另外,值得注意的是,蒸散發(fā)估算值雖然整體被低估,但存在季節(jié)性差異(圖3)。為進(jìn)一步分析不同月份LE與LE散點(diǎn)關(guān)系的差異,本研究分析了日尺度和30 min 尺度上欒城站、禹城站和館陶站LE(軸)與LE(軸)之間的散點(diǎn)關(guān)系(圖7、圖8)。從圖7 可以看出,欒城、禹城和館陶3 站潛熱通量計(jì)算值和觀測(cè)值有較好的相符性,決定系數(shù)為0.80～0.83,3 站綜合為0.81; 一元線性擬合斜率低于1,且大部分點(diǎn)低于1∶1 線,表明蒸散發(fā)估算值整體被低估。值得注意的是,欒城站和禹城站偏離出95%預(yù)測(cè)帶的值集中于6月和7月,僅低估?1.38 W?m,并且主要表現(xiàn)為高于預(yù)測(cè)區(qū)間; 而館陶站6、7月份低估?23.5 W?m。這種現(xiàn)象的主要原因可能是欒城站和禹城站下墊面為冬小麥、夏玉米輪作,而館陶站是冬小麥、夏玉米輪作和棉花的混合下墊面。6月中旬到7月上旬是冬小麥?zhǔn)斋@后,而夏玉米七葉期之前的階段,該階段作物葉面積指數(shù)較低,尚不能覆蓋地表,下墊面性質(zhì)類似裸地到稀疏草地之間。在植被生長(zhǎng)茂盛時(shí),NP 方法低估LE,在裸地或植被稀疏時(shí)高估LE 的現(xiàn)象與上述文獻(xiàn)結(jié)果可相互驗(yàn)證。

圖7 日尺度上3 個(gè)站點(diǎn)的潛熱通量(LE)估算結(jié)果與觀測(cè)結(jié)果的比較Fig.7 Comparison of latent heat flux (LE) between estimation and observation at daily scale in the three stations

圖8 30 min 尺度上3 個(gè)站點(diǎn)的潛熱通量(LE)估算結(jié)果與觀測(cè)結(jié)果的比較Fig.8 Comparison of latent heat flux (LE) between estimation and observation at 30 min scale in the three stations

NP 方法所產(chǎn)生的殘差可能受到若干不確定來源的影響,包括NP 方法固有局限性、能量閉合通量校正、土壤熱通量校正以及氣象變量的測(cè)量誤差。能量閉合修正的影響已在上文討論。對(duì)測(cè)量深度的進(jìn)行校正可以改善2%～5%的能量閉合,進(jìn)而NP方法估算精度評(píng)價(jià)。地表溫度1 K 的誤差大約會(huì)導(dǎo)致13 W·m的誤差。的測(cè)量誤差通常小于5%。測(cè)量精度通常較高,但因?yàn)闇y(cè)量高度的差異,換算成2 m 的標(biāo)準(zhǔn)高度時(shí),會(huì)產(chǎn)生一定誤差。而從NP 方法固有局限性分析,其原因可能(但不完全)包含:1)該方法假設(shè)在局部熱平衡(=)情況下,蒸散發(fā)等于平衡蒸發(fā),未考慮平流項(xiàng)的影響。本研究發(fā)現(xiàn)公式(2)的第3 項(xiàng)數(shù)值極小,對(duì)估算結(jié)果幾乎無影響。在灌溉農(nóng)田站,第2 項(xiàng)所占數(shù)值比例也較低,該區(qū)域LE 低估主要是因?yàn)榈? 項(xiàng)低估。也許NP 方法可以參考Priestley ‐ Taylor 公式加入一個(gè)系數(shù),但這樣又改變了無參數(shù)化方法的初衷。因此,可能需要進(jìn)一步研究NP 方法與其他方法(如Priestley ‐ Taylor 公式、Penman-Monteith 公式)之間建立物理和數(shù)學(xué)上的聯(lián)系。2)可能與輸入數(shù)據(jù)和模型參數(shù)化中的錯(cuò)誤有關(guān),例如該方法使用代替了空氣動(dòng)力學(xué)阻力、表面阻力和土壤濕度的影響,這種空氣動(dòng)力學(xué)阻力和表面阻力的簡(jiǎn)化可能會(huì)削弱模型的性能。因?yàn)樵诓煌诤妥魑镱愋捅砻孀枇εc的關(guān)系也許不是固定線性關(guān)系,同時(shí)在氣象條件不穩(wěn)定、下墊面變化劇烈(灌溉、作物收獲等)情況下,這種替代關(guān)系可能并不總是適用。發(fā)展溫差(?)與地表溫度的線性關(guān)系可能是在不久的將來克服這一限制的一種可行方法?？偟膩碚f，由于所有輸入?yún)?shù)少,且易獲取/計(jì)算,NP 方法有可能成為一種常規(guī)監(jiān)測(cè)地表熱通量的方法; 利用常規(guī)衛(wèi)星反演或者再分析資料,可能在大區(qū)域尺度甚至全球尺度提供一個(gè)蒸散發(fā)估算方法; 也有助于開發(fā)廣泛應(yīng)用于氣象、水文、氣候模式的陸面方案。目前,NP 方法用于其他下墊面類型(如城市混合下墊面、森林、濕地等)的適用性仍不明確,有待進(jìn)一步研究。另外,NP 方法未來的研究將集中于計(jì)算方法、輸入精度和殘差的改進(jìn)上。

5 結(jié)論

本文利用非參數(shù)化蒸散發(fā)估算方法估算了日尺度和30 min 尺度華北灌溉農(nóng)田3 個(gè)典型站點(diǎn)的潛熱通量,基于通量觀測(cè)數(shù)據(jù)分析了該方法在不同季節(jié)及不同時(shí)間尺度的適用性,闡明了各輸入?yún)?shù)的敏感性,探索了能量閉合修正對(duì)蒸散發(fā)估算結(jié)果驗(yàn)證的重要影響。研究主要結(jié)論:

1)非參數(shù)化估算方法在華北灌溉農(nóng)田不同作物類型、不同時(shí)間尺度具有可靠和穩(wěn)健的表現(xiàn),可以較好地反映季節(jié)及日內(nèi)變化特征,在日尺度和30 min 尺度上總體表現(xiàn)較好。在日尺度上NP 估算的LE (LE)與能量殘差法修正后的渦度相關(guān)技術(shù)測(cè)定 的LE (LE)之 間 的Bias、RE、、RMSE 和NSE 分 別 為?16.18～?12.88 W·m、 20.6%～30.2%、0.80～0.83、22.45～31.06 W·m和0.66～0.75,在30 min尺度上分別為?13.30～?7.68 W·m、17.7%～24.8%、0.88、39.22～42.15 W·m和0.86～0.87。

2)非參數(shù)化估算方法在水分供應(yīng)較充足時(shí)較嚴(yán)重低估潛熱通量,而在較干燥時(shí)輕度低估或不低估潛熱通量; 該方法對(duì)蒸散發(fā)對(duì)灌溉活動(dòng)的響應(yīng)考慮不足。從整體上看,非參數(shù)化估算方法在灌溉農(nóng)田中普遍低估潛熱通量,欒城、禹城和館陶3 個(gè)站點(diǎn)年尺度上平均低估14.69 W?m(RE:27.5%),30 min尺度上平均低估10.34 W?m(RE:21.5%)。從作物類型上看,非參數(shù)化估算方法在玉米季(7?10月,供水以降水為主)估算精度高于小麥季(11月至次年6月,需大量灌溉補(bǔ)充降水的不足),3 個(gè)站點(diǎn)玉米生長(zhǎng)季平均為0.80,小麥生長(zhǎng)季平均為0.42。從生長(zhǎng)周期上看,作物生長(zhǎng)茂盛時(shí),非參數(shù)化估算方法低估潛熱通量,裸地或植被稀疏時(shí)幾乎不再低估潛熱通量,在欒城和禹城站6、7月份(裸地或作物稀疏)僅低估?1.38 W?m,同時(shí)館陶站(作物生長(zhǎng)良好)低估?23.5 W?m。

3)非參數(shù)化估算方法在華北灌溉農(nóng)田中參數(shù)敏感性從高到低依次為、和和,日尺度上這4 個(gè)參數(shù)變化5%時(shí)的敏感性系數(shù)分別為11.3%、?10.4%、6.0%和?0.24%; 30 min 尺度上這4 個(gè)參數(shù)變化5%時(shí)的敏感性系數(shù)分別為11.75%、?10.94%、6.24%和?0.38%。

4) Bowen 比修正方法和能量殘差修正方法可顯著提高非參數(shù)估算方法模擬的一致性，其中可分別提高0.23 和0.25,RMSE 可分別降 低6.61 W·m和5.91 W·m,NSE 可分別提高0.27 和0.29。