趙子敬，魏國孝，劉紅娟，田強龍，邱中齊，王福兵，楊澤偉

（蘭州大學資源環(huán)境學院，蘭州 730000）

0 引 言

幾乎所有的陸地-大氣耦合過程都受到陸地水儲存和通量的顯著調(diào)節(jié)，這些儲存和通量也會影響碳、養(yǎng)分、生態(tài)系統(tǒng)結(jié)構(gòu)和功能，并限制社會的淡水供應[1]。蒸散發(fā)（ET）包括土壤蒸發(fā)（ETs）和植物散發(fā)（ETc）是氣候系統(tǒng)的中心過程，是水、能量、循環(huán)的紐帶[2]。ET的準確評估對有效灌溉和實現(xiàn)水資源精細化管理有重要作用[3]。ET估算的準確性與模型的選擇和參數(shù)有關，并且取決于對氣候和生物因素的準確評估。因此，了解氣候和植被對ET的相對影響對于預測水循環(huán)將如何響應未來的氣候和生物變化具有重要意義[4]。

ET受多種因素控制，有研究表明，ET的變化是由太陽輻射、空氣溫度、水汽壓、風速等變量驅(qū)動的[5-7]。不同氣候區(qū)域?qū)T影響最大的氣候因素不同，影響ET過程的不同因素之間復雜的相互依賴關系使得準確量化ET具有挑戰(zhàn)性[8]。有些研究對蒸散發(fā)變化趨勢主要集中在氣候因素上，而沒有考慮植被的作用。然而，植被在蒸散發(fā)中起著至關重要的作用，尤其是通過植物控制ETc的過程。目前有很多方法可以分離氣候和植物驅(qū)動因素對ET的影響[4,9,10]。Stoy等[4]將渦動相關方差測量與線性擾動分析方法相結(jié)合，分離出美國東南部3個生態(tài)系統(tǒng)中物理和生物因素對ET的相對貢獻。YANG 等[10]對PenPan 蒸發(fā)模型進行了修正，推導出一階偏微分形式的蒸發(fā)皿蒸發(fā)方程，以量化氣候因子的貢獻。LI 等[9]利用Penman-Monteith模型、修正的作物系數(shù)法、Priestly-Taylor 模型、線性回歸模型解釋了2007-2013年中國西北玉米地ET年際變化。在該研究中，他們證實了植被驅(qū)動因子在調(diào)控ET過程中起著重要作用，強調(diào)準確估算冠層導度在ET建模和預測中的重要性，并為區(qū)分氣候和植被對ET變化的貢獻提供了一種新方法。

為了明確考慮通量相互作用，Shuttleworth 和Wallace[11]提出了一個雙源模型(以下稱為“SW 模型”)，該模型結(jié)合了兩個單源模型:一個是土壤蒸發(fā)(ETs)，另一個是植物散發(fā)(ETc)。SW 模型旨在捕捉稀疏植被條件下葉片下的裸土狀況，考慮冠層和土壤對總能量的貢獻。在模型中，植物冠層阻力(rcs)和土壤阻力(rss)和冠層高度與參考水平之間的空氣動力阻力(raa)以及葉片到冠層高度和土壤表面到冠層高度之間的空氣動力阻力(rca和rsa)分別通過控制地表與大氣之間的能量傳遞來調(diào)節(jié)植物的能量傳遞[12]。然而，SW 模型中復雜參數(shù)和模型結(jié)構(gòu)的不確定使得SW 模型模擬ET具有挑戰(zhàn)性。近年來，在開發(fā)用于驅(qū)動和處理參數(shù)和模型結(jié)構(gòu)不確定性的系統(tǒng)框架方面取得了重大進展[13,14]。特別是，一種新的自適應馬爾可夫鏈蒙特卡羅(MCMC)算法，又稱為差分進化自適應都市(Differential Evolution Adaptive Metropolis )DREAM 算法，被用來在貝葉斯框架內(nèi)有效估計驅(qū)動數(shù)據(jù)和參數(shù)不確定性[15,16]。所以根據(jù)貝葉斯原理利用DREAM 算法優(yōu)化SW 模型參數(shù)，這不僅簡化了參數(shù)的獲取，還可以降低參數(shù)的不確定性。偏導數(shù)方法可以量化驅(qū)動因子對ET變化的相對貢獻，所以本研究將SW 模型中的ETs和ETc分開，定量評估了冠層上方和土壤表面的可用能量(A和As)、飽和水汽壓差(es-ea)和阻力raa、rsa、rca、rcs、rss在ET變化中的貢獻，從而確定ET變化的主要驅(qū)動因子。

目前，利用SW 模型估算ET時很少考慮阻力的變化，將ET的年際變化歸因于這些阻力的研究也很少受到關注。本研究旨在探討大滿站玉米生態(tài)系統(tǒng)蒸散發(fā)變化的潛在驅(qū)動因素，利用偏導數(shù)方法量化氣候因素和植被因素對ETs和ETc變化的相對貢獻。而阻力因子主要受到氣候和植物因素的影響。所以考慮了阻力驅(qū)動因子的變化，并對SW 模擬的ETs和ETc的變化進行歸因分析。這一結(jié)果為了解玉米ET的年際變化提供了有用的見解，并強調(diào)了阻力在調(diào)控ET中的重要作用。

1 材料與方法

1.1 研究區(qū)概況

為了評估不同驅(qū)動因子對ETs和ETc年際變化的影響，本研究選擇了2016-2018年7-9月份大滿通量站的數(shù)據(jù)。玉米生長試驗是在中國西北部甘肅省張掖市的大滿站進行的。大滿綠洲位于黑河流域中游，是中國西北干旱區(qū)第二大內(nèi)陸河流域。黑河流域中游地區(qū)以綠洲灌溉農(nóng)業(yè)為特征，是生活用水和農(nóng)業(yè)用水量大的地區(qū)。它的年平均降水量和氣溫分別為125 mm 和7.2 °C（1960-2000年），年平均潛在蒸發(fā)量約為2 290 mm，年平均日照時數(shù)3 106 h，無霜期148 d。土壤類型以粉質(zhì)黏土為主，凍土層深度約143 mm。研究區(qū)為典型的灌溉農(nóng)業(yè)區(qū)，主要水源為祁連山融雪。玉米和春小麥是該地區(qū)種植的主要農(nóng)作物。玉米一般在4月下旬播種，9月中旬收獲，種植行距40 cm，株距30 cm，研究區(qū)作物密度約為6.6 萬株/hm2。大滿站的實驗位置和研究地點見圖1。

圖1 大滿站的實驗位置和研究地點Fig.1 Experiment locations and study sites at Daman Station

1.2 測量與處理

大滿站的觀測數(shù)據(jù)是從黑河流域聯(lián)合遙測實驗研究（HiWATER）項目的野外觀測系統(tǒng)中收集的，詳細描述參考李等[17]。渦度（EC）相關系統(tǒng)的架設高度為4.5 m。渦度數(shù)據(jù)中包含潛熱通量和感熱通量，均為半小時尺度，數(shù)據(jù)處理包括：野值點剔除，延遲時間校正，三維旋轉(zhuǎn)[18]，頻率響應修正，超聲虛溫修正和密度校正等，頻率響應校正等[19]。并刪除由于降水、水汽凝結(jié)、系統(tǒng)故障導致的虛假數(shù)據(jù)。站點儀器測量的數(shù)據(jù)主要有降雨量、氣溫、風速、大氣壓強、土壤水分、土壤熱通量、凈輻射等數(shù)據(jù)。在生長季，葉面積指數(shù)(LAI)大約每10 d測量一次。

1.3 模型與方法

1.3.1 Shuttleworth-Wallace(SW)模型

SW[13]模型包括植物蒸騰的一維模型和土壤蒸發(fā)的一維模型。在模型中，地表和冠層阻力調(diào)節(jié)土壤和植物表面之間的質(zhì)量和熱量傳遞，而空氣動力阻力調(diào)節(jié)兩表面和大氣邊界層之間的阻力[13,14]。土壤蒸發(fā)和植物蒸騰用以下公式[13]計算：

式中：λET為來自植物（λT）和土壤（λE）的潛熱通量之和，W/m2；Δ 為飽和水汽壓曲線斜率，kPa/K；γ為濕度常數(shù)，kPa/K；Cp為干燥空氣的比熱容，1 013 J/(kg?K-1)；es和ea分別為飽和水汽壓和實際水汽壓，kPa；raa為植物冠層高度與參考高度之間的空氣動力學阻力，s/m；rca為冠層內(nèi)邊界阻力，s/m；rsa為地面與冠層高度間的空氣動力阻力，s/m；rcs為植物群體冠層阻力，s/m；rss為下墊面土壤表面阻力，s/m。

A和As(W/m2)分別為冠層和土壤表面以上的可用能量，定義為：

式中：Rn和Rns分別為冠層上方和地面接收的太陽凈輻射，W/m2；G為土壤熱通量，W/m2；Ka為消光系數(shù)；LAI為葉面積指數(shù)。阻力模型的計算參考魏等[1]。

1.3.2 貝葉斯推理和DREAM 算法

參數(shù)的后驗概率分布由貝葉斯定理計算：

式中：π(θ|M)表示模型M下的θ的先驗密度；p(D|θ,M)是模型M及參數(shù)θ的聯(lián)合似然性。邊際似然或貝葉斯模型證據(jù)是：

用于參數(shù)估計的似然函數(shù)p(D|θ,M)是根據(jù)觀測誤差的分布指定的。在時間t的每個觀測值D（t）中的誤差有e( )t表示為：

假設e(t)服從均值為零的高斯分布，似然函數(shù)可以表示為：

式中：n是觀測次數(shù)；σ表示誤差方差。

本文使用DREAM 算法來探索ET模型的參數(shù)空間，并優(yōu)化SW 模型中的參數(shù)。DREAM 采樣方案是對混洗復雜進化Metropolis (SCEM-UA)的全局優(yōu)化算法的改編。Vrigt 等[15,16]更詳細地描述了該算法。

1.4 評估模型性能的傳統(tǒng)統(tǒng)計指標

用于評估模型性能的傳統(tǒng)誤差度量包括決定系數(shù)（R2）、一致性指數(shù)（IA）、相對誤差度量（EF）、均方根誤差（RMSE）、平均偏差誤差（MBE）。

1.5 不同驅(qū)動因子對ET變化趨勢的貢獻

本文采用偏微分方程研究氣候和植被因素對ET變化趨勢的貢獻。通過這種方法，可以分離出驅(qū)動生長季蒸散發(fā)變化的氣候和植物因素的相對貢獻。李等[9]使用的Penman-Monteith模型將ET表示為4個氣候變量和一個植物變量的連續(xù)函數(shù):

本研究把SW 模型中的ETs和ETc分別表示為兩個不同的連續(xù)函數(shù)，并假設各個變量是相互獨立的。用偏導數(shù)計算ETs和ETc的變化，并分析各驅(qū)動因子對ETs和ETc的相對貢獻。并用作物潛熱通量(λT)和土壤潛熱通量(λE)(W/m2)來表示ETc和ETs。

λE可以表示為A、As、es-ea、raa、rsa、rss的連續(xù)函數(shù)：

同樣地，λT可以表示為A、As、es-ea、raa、rca、rcs的連續(xù)函數(shù)：

λE和λT的相對變化可以用泰勒展開級數(shù)表示為一階偏微分形式。

λE的相對變化：

類似的,λT的相對變化可以表達為：

λET是通過渦度觀測系統(tǒng)測得的，并利用SW 模型和實測氣象資料計算λE和λT驅(qū)動因子的偏導數(shù)。為方便表述，下文用LE表示潛熱通量。

2 結(jié)果與分析

2.1 能量閉合分析

所有表面能量和質(zhì)量交換模型都基于基本守恒定律；即能量守恒和質(zhì)量守恒。能量守恒方程的主要組成部分，通常稱為“能量平衡閉合”，可以描述為：

式中：H為顯熱湍流通量，LE為潛熱湍流通量，Rn為凈輻射，G為地表熱通量。能量平衡比（EBR），定義為H和LE之和與Rn和G之差的比值。在大多數(shù)站點EBR為70%～90%，這稱之為能量不閉合問題。

2016-2018年基于半小時白天湍流通量（H+LE）總和與可用能量（Rn-G）的線性回歸見圖2。由于晚間風速小，或者湍流不發(fā)生，以平流的方式傳遞通量，導致EC測得的通量嚴重不閉合。因此，本研究只選擇白天的數(shù)據(jù)，并定義“Rn-G>0”為白天。結(jié)果表明，2016-2018年的EBR分別為0.75、0.89、0.84，決定系數(shù)（R2）也都在0.7 以上，說明擬合效果良好。

圖2 2016-2018年能量平衡閉合Fig.2 Energy balance closure from 2016 to 2018

2.2 參數(shù)優(yōu)化及模型性能

SW 模型中有7 個參數(shù)gmax、Q50、D50、Kq、Ka、b1和b2（表1）。每個參數(shù)的先驗概率密度被指定為均勻分布。2016-2018年7-9月的數(shù)據(jù)用于估計每個參數(shù)的后驗概率密度函數(shù)，并使用DREAM 算法生成了50 000N個樣本來實現(xiàn)。在計算中，鏈數(shù)N等于相關模型中的參數(shù)數(shù)，對于SW 模型，N等于7。對于每個鏈，前10 000 個樣本作為老化數(shù)據(jù)丟棄，其余40 000個樣本用于設置后驗密度函數(shù)。優(yōu)化后的SW 模型參數(shù)值的最大似然估計見表1。

表1 模型參數(shù)優(yōu)化后的最大似然估計Tab.1 Maximum likelihood estimation of model parameters after optimization

使用2016-2018年白天每半小時數(shù)據(jù)評估了優(yōu)化后的SW模型性能。使用模型的校準參數(shù)，并運行SW 模型來估計每半小時的LE值。SW模型模擬與測量LE值的線性回歸見圖3。一般來說，SW 模型與測量的LE有較好的擬合效果(R2>0.77)。圖3可以看出，雖然模型和觀測值仍存在一定的差異，但SW 模型能夠捕捉到蒸散發(fā)的動態(tài)變化。

圖3 2016-2018年SW 模型和實測LE線性擬合Fig.3 The SW model was fitted with the observed LE from 2016 to 2018

半小時蒸散量實測值與模擬值的R2、RMSE、IA、EF、MBE見表2?；趯@些傳統(tǒng)誤差度量的分析，這三年的IA和EF比較接近，而MBE分別為-9.61、13.82、-19.11。從總體上來看經(jīng)過參數(shù)優(yōu)化后的SW模型能夠準確模擬ET。

表2 模型性能的傳統(tǒng)統(tǒng)計指標Tab.2 Traditional statistical metrics of model performance

2.3 ET及驅(qū)動因子的年平均變化

2016-2018年SW 模擬的ET、ETs、ETc、A、As、es-ea、raa、rsa、rss、rca、rcs白天每半小時的平均值見表3。2016-2018年ET的實測值分別為195.62、232.17、213.66 W/ m2，從表3可

表3 2016-2108年白天每半小時的均值Tab.3 Mean values for every half hour during daytime from 2016 to 2108

以看出2016年和2017年SW 模型模擬的ET均值與實測ET相差不大，2018年則模擬的偏大。2018年ETs為138.07 W/m2，明顯比前兩年偏大，主要是因為As模擬值偏大，導致分配給冠層的可用能量偏小，所以ETc的值就偏小。阻力方面，其中raa、rsa、rca三年變化不明顯。2017年rss值為317.13 s/m，明顯比2016年134.37 s/m 和2018年85.17 s/m 大；2018年rcs值為217.73 s/m，與其余兩年相比偏大。這些差異可能是導致玉米農(nóng)田ET年際變化的潛在因素。

2.4 驅(qū)動ETs和ETc變化的因素

由A、As、es-ea、raa、rsa、rss對ETs變化的貢獻率和A、As、es-ea、raa、rca、rcs對ETc變化的貢獻率見圖4。它由方程式（14）和（15）得到，貢獻率的總和應該等于1，然而由于各種不確定的因素導致結(jié)果并不嚴格等于1。從圖4可以看出ETs與A、As、es-ea、raa呈正相關與rsa和rss呈負相關。與2016年相比，ETs分別增加了2.56 W/m2、72.9 W/m2（表3），這主要是因為As的增加。2017年和2018年As分別增加了55.09 W/m2和156.93 W/m2。ETs的增加主要歸因于As、es-ea、A的變化（圖4），并且rsa和rss的變化抵消了部分對ETs的正貢獻。其中，對ETs變化最大的正相關驅(qū)動因子為As，最大負相關驅(qū)動因子為rss，分別貢獻了124%、95.4%和-99.1%、-26.9%。對于影響ETc變化的驅(qū)動因子對2017 和2018年的貢獻有較大區(qū)別。其中，2017年ETc與A、es-ea呈正相關，與As呈負相關；而在2018年，ETc與A、es-ea呈負相關，與As呈正相關。rcs對于ETc的變化都是正貢獻，并且占了很大比例，分別為68.8%和71.9%，間接說明植物因素在調(diào)節(jié)ETc中起著不可忽視的作用。

圖4 2016-2018年不同驅(qū)動因子對ETS和ETC變化的相對貢獻Fig.4 The relative contributions of the different driving factors to changes in ETs and ETc,from 2016 to 2018

3 討 論

3.1 能量閉合分析

渦度相關（EC）法是基于對垂直速度波動（w′）和標量濃度波動（c′）乘積的直接測量，并假設平均垂直速度可以忽略不計，從而直接估計H和LE；Rn和G分別是通過凈輻射儀和土壤熱通量儀測量的。表面能量收支不閉合的原因有很多，Mahrt[20]討論了幾個原因，其中包括：①在非常接近表面的地方測量的各種通量分量之間的源區(qū)域缺乏重合，例如葉子的蒸發(fā)和來自炎熱干燥土壤表面的顯熱；②測量時間為典型的30 min可能會錯過由非常低的頻率波動引起的協(xié)方差；③儀器位置的平均垂直速度不可能系統(tǒng)的為零，從而產(chǎn)生垂直平流通量；④與傳感器分離、頻率響應、對準問題以及來自塔或儀器安裝結(jié)構(gòu)的干擾有關的測量誤差等。而Foken[21]表明以前關于測量誤差或存儲項是未閉合能量平衡的原因的假設站不住腳，因為即使是來自校準傳感器的湍流通量、凈輻射和地面熱通量也無法閉合能量平衡。作為一個假設，他將能量平衡閉合問題假設為一個尺度問題，并且閉合只可能在景觀尺度上發(fā)生。由于能量不閉合的原因至今沒有一個明顯的定論，并且在極端情況下，整個能量不閉合可能歸因于潛熱通量[22]。研究結(jié)果表明，H+LE的真實表面通量可能比EC 測量值高20%～30%。能量不閉合可能會對SW 模型的不確定性產(chǎn)生影響，模型使用的H值比實際值高，從而導致模擬的LE偏高。

3.2 不確定性和誤差分析

關于DREAM 算法的效率，SW 模型的接受率遠高于之前研究中使用的一些MCMC 算法獲得的接受率[23]。DREAM 算法可以有效地處理涉及高維、多模態(tài)、非線性的問題。使用DREAM 算法優(yōu)化的參數(shù)Ka，可以明確地將總可用能量劃分為SW 模型中冠層和土壤吸收的能量。方程（5）的分析發(fā)現(xiàn)Ka的變化不僅可以解釋消光效應，還可以糾正能量驅(qū)動數(shù)據(jù)的錯誤。這也意味著使用標定數(shù)據(jù)估算的Ka值實際上不僅僅是真正的消光系數(shù)，還包括SW 模型中的能量不平衡校正。從這個分析中，可以看到，Ka不僅涉及到冠層和土壤表面之間的能量分布，還涉及到能量的不平衡。因此，參數(shù)Ka對SW 模型的性能有很大的影響。此外，由于Ka影響冠層和土壤的可用能量分配，Ka值較低導致分配給土壤表面可用能量偏大，從而導致模擬的ETs偏大。

試驗結(jié)果可能存在不確定因素。首先，對ETs和ETc估計的誤差可能歸因于在泰勒展開式中考慮的不同驅(qū)動因子獨立性的假設。事實上，這些變量并不是完全獨立的，它們確實在一定程度上相互影響，這導致試驗分析存在一定的不確定性。其次，這些變量并不是固定的值，而是隨著環(huán)境的變化而變化，同時也會受到測量誤差的影響。再次，泰勒偏微分方程適用于兩個比較相近的值，由于所考慮的變量的年變異性造成的差異導致了較大的誤差。在本研究中，沒有考慮人為因素和下墊面條件的影響，因此對每個變量的分析引入了一些不確定性因素。除了測量和計算誤差外，模型的選擇還存在一定的不確定性。以往的研究表明，冠層阻力和空氣動力阻力模型的不同會導致不同的模擬結(jié)果[24]。此外，模型的結(jié)構(gòu)也是不確定性的進一步來源。本研究只是用DREAM 算法優(yōu)化SW 模型參數(shù)使模型達到良好的擬合效果，并沒有土壤蒸發(fā)和植物散發(fā)的實測值，這也是導致結(jié)果不確定的重要因素。

3.3 歸因分析

本文利用泰勒展開法求解偏微分方程，分析了氣候和植被變量對SW 模型中ETs和ETc的相對貢獻。結(jié)果表明，土壤表面可用能量(As)是影響ETs變化最主要的外因，這或許與分配到土壤表面的凈輻射有關。這與ZHANG 等[25]揭示凈輻射下降是影響青藏高原東部蒸發(fā)皿蒸發(fā)的主要因素相符。而LIU 等[5]認為風速是影響是青藏高原蒸發(fā)皿蒸發(fā)的主要因素，這可能與風速影響了空氣動力阻力和飽和水氣壓差有關；并且李等[26]發(fā)現(xiàn)空氣動力學項是造成東北氣象站參考作物蒸散發(fā)量下降趨勢的主要原因。由此可見，ETs年際變化主要受氣候變化的影響，尤其是As的影響。對于ETc而言，rcs變化是導致ETc變化的主要植物變量。Wilson 等[27]研究了北美溫帶闊葉林ET的年際變化，結(jié)果表明不同年份ET的變化主要是由于冠層導度造成的，其次是凈輻射的變化。LI等[28]表明冠層導度主導了能量分配過程，然后發(fā)現(xiàn)ET的年際變化主要歸因于冠層導度的變化而不是氣候因素的變化[9]。主要是因為冠層導度是調(diào)節(jié)ETc的重要因素。這些結(jié)果表明，植物驅(qū)動因素對ET有重要影響，在未來氣候條件下模擬蒸散發(fā)應更加重視冠層導度的作用。

4 結(jié) 論

（1）本文利用DREAM 算法優(yōu)化了SW 模型參數(shù)，降低了大部分參數(shù)的假設先驗不確定性，并達到良好的模擬效果。

（2）泰勒偏微分方程法不僅可以得到主導因子，還可以評估各因子的貢獻率。本研究應用此方法分析驅(qū)動因子對玉米ETs和ETc年際變化的貢獻，歸因分析表明As是影響ETs變化最主要因素，除了As外，es-ea和rss對ETs變化也起著很重要的作用。

（3）相對于氣候驅(qū)動變量，植物控制在調(diào)節(jié)ETc方面發(fā)揮了更大的作用。關于阻力對ETc的影響，由rcs變化引起的ETc的變化最大，說明相較于其他阻力而言，rcs對ETc的影響最大。

（4）一般來說，ET的變化受氣候和植物變量的復雜相互作用的影響，As和rcs分別是控制玉米ETs和ETc年際變化的兩個主要因素。