考慮地表反照率時域重建的海河流域蒸散發(fā)模擬研究

段 浩，韓 昆，趙紅莉，徐浩瑋，文 銘，靳曉輝，蔣云鐘

（1.中國水利水電科學研究院，北京 100048；2.黃河水利科學研究院，河南 鄭州 450003）

1 研究背景

地表蒸散發(fā)（ET）是聯(lián)系水、能量和碳循環(huán)的中間環(huán)節(jié)［1］，是陸面生態(tài)過程［2］和水資源管理［3］的關鍵參數(shù)。準確獲取ET時空分布，對理解陸氣間水熱交換、加強水資源管理具有重要意義。

基于遙感監(jiān)測技術，可獲取ET在大時空尺度上的變化特征，與傳統(tǒng)的站點觀測相比具有顯著優(yōu)勢。在Penman-Monteith（P-M）公式基礎上發(fā)展而來的Penman-Monteith-Leuning（PML）模型［4］，是基于遙感數(shù)據(jù)計算ET的方法之一，該方法通過引入地表導度Gs，可基于植被葉面指數(shù)和P-M公式直接計算ET，已在全球不同氣候條件下得到了應用。王海波等［5］基于PML模型，模擬了黑河流域ET；李紅霞等［6］和Zhang等［7］分別利用此模型對澳大利亞的ET進行了模擬和驗證；Li等［8］利用該模型分析了植被變化對ET的影響；Zhang等［9］還基于PML模型分析了全球范圍內氣象要素對土壤蒸發(fā)和植被散發(fā)的影響。這些研究使PML模型在ET模擬上逐漸成熟，成為開展大尺度ET監(jiān)測的重要手段［10］。

地表反照率（ac）是反映地表對太陽輻射反射能力的指標，它的細微變化都會對輻射（Rn）收支分配產生顯著影響［11］，最終影響ET的模擬精度。在利用PML模型進行ET模擬時，通常通過兩種方式對ac進行處理，一種是直接使用8天合成的ac產品，如Zhang等［12］對MODIS的8天合成ac產品進行了質量控制和平滑處理，在重采樣后得到0.05°分辨率的數(shù)據(jù)；該作者在計算全球ET時［9］同樣是使用了Global Land Surface Satellite（GLASS）8天0.05°的ac產品。這種處理方式只能得到每8天的ET結果。另一種方式是對8天合成的ac產品做線性插值，得到逐日數(shù)據(jù)，如Li等［8］對8天合成的ac進行線性插值，分析了子牙河流域ET對植被變化的響應；Mu等［13］在計算ET時，同樣使用了線性插值的方法，得到逐日的ac產品?；诰€性插值的方法獲取時間連續(xù)的ac產品，處理方式簡單，但該方法將ac在時間上的演變簡化為線性關系，機理性不足。ac的變化受到土壤含水量［14］、下墊面條件［15］、積雪［16］、溫度［17］等多種要素的綜合影響，具有顯著的時空異質性特征，線性插值難以體現(xiàn)其變化過程。張振宇等對SEBAL蒸散發(fā)模型的研究表明［18］，通過不同方法得到的ac間的不同，是導致日蒸散量差異的主要原因。

為強化PML模型模擬中ac產品的遙感機理，本文引入充分考慮遙感數(shù)據(jù)時空信息的諧波分析與泊松方程協(xié)同的時域重建算法，在充分考慮ac時空變異特征的基礎上對海河流域ac產品進行重建，并基于時域重建后的ac產品進行流域蒸散發(fā)反演。同時基于線性插值方法得到的ac，分析不同方法得到的ac對PML模型中ET結果的影響。

2 研究方法

2.1 基于葉面指數(shù)的PML蒸散發(fā)模型 PML模型計算ET的基本思想是，在潛在蒸發(fā)計算方法基礎上引入“表面阻抗”的概念［19］，得到非飽和下墊面蒸散發(fā)的計算形式，其計算過程可拆分為土壤蒸發(fā)和植被蒸騰兩部分，具體可表示為［4］：

其中：ET為蒸散發(fā)；λ為汽化潛熱；ε為溫度-飽和水汽壓曲線斜率與干濕表常數(shù)的比值；ρa為空氣的密度；Cp為空氣定壓比熱；Da為參考高度飽和水汽壓差；A為可用能量，代表凈輻射與土壤熱通量的差值；Ga為空氣動力學導度；Gs為地表導度；Gc為冠層導度；f為土壤蒸散發(fā)系數(shù)；Ac和As分別為可用能量中被冠層和土壤吸收的部分。

冠層導度的計算需要利用最大氣孔導度gsx和LAI的關系得到：

其中：kQ為短波輻射衰減系數(shù)，常取0.6；Qh為冠層上方的可見光輻射通量；Q50和D50為當氣孔導度gs=gsx/2（gsx是gs的最大值）時的可見光輻射通量和水汽壓差，通常取值分別為2.6 MJ/（m2/d）和0.8 kPa。

可用能量A為凈輻射Rn與土壤熱通量G之差，而Rn則由ac及太陽短波輻射和凈長波輻射得到，計算方法為：

其中Rs為太陽短波輻射（MJ·m-2·d-1），Rnl為凈長波輻射（MJ·m-2·d-1）。

在對f和gsx兩個參數(shù)的處理上，gsx的參數(shù)化參考Zhang等的研究，根據(jù)土地利用類型確定［9］；在f的參數(shù)化方面，本文采用Zhang［12］提出的基于累積降雨與平衡蒸發(fā)速率比值的方案來確定f的取值。

2.2 地表反照率時域重建方法 ac數(shù)據(jù)對ET的模擬有重要影響，為充分考慮ac的時空變異性，增強其處理過程的物理機理，本文引入諧波分析與泊松方程協(xié)同的ac時域重建算法［20］，計算得到海河流域逐日的ac產品。時域重建的主要過程如下：

（1）初始值確定。指定目標年（模擬年份）和參考年的ac數(shù)據(jù)序列，本文選取MODIS的ac產品進行分析，對每一個目標年數(shù)據(jù)，以其前三年ac數(shù)據(jù)為參考年數(shù)據(jù)。首先根據(jù)ac的數(shù)據(jù)質量信息對全部原始數(shù)據(jù)做質量控制，區(qū)分有效區(qū)域和數(shù)據(jù)缺失區(qū)域；然后利用矩匹配的方法對參考年數(shù)據(jù)做歸一化，并分析與目標數(shù)據(jù)的相關性，得到參考年的權值；最后通過多年參考數(shù)據(jù)的加權處理初步計算缺失區(qū)域的ac。

關于參考年數(shù)據(jù)的歸一化處理，又可分為以下三種情況：

任意年份參考數(shù)據(jù)與目標數(shù)據(jù)的有效區(qū)域均有重疊，此時以目標年有效區(qū)域重疊部分的均值和標準差為基準，進行歸一化；任意年份參考數(shù)據(jù)與目標數(shù)據(jù)的有效區(qū)域均無重疊網格，此時無法計算權值，基于多年均值進行數(shù)據(jù)缺失區(qū)域的填補；存在參考年數(shù)據(jù)與目標數(shù)據(jù)有效區(qū)域都不存在重疊的部分情況，根據(jù)缺失區(qū)域的質量標識和在上述兩種情況下得到的臨時像元值來確定。

在上述各情況下初始值的具體計算過程，可參見文獻［20］的相關描述。

（2）時域插值。當同一網格在參考年的不同年份均為空值，則通過時域插值來填補步驟（1）未能填補的區(qū)域，采用真值諧波分析（HANTS）法進行，其基本算法如下：

3 研究數(shù)據(jù)

3.1 遙感數(shù)據(jù) 研究所用到的遙感影像數(shù)據(jù)包括MODIS的ac數(shù)據(jù)MCD43A3、LAI數(shù)據(jù)MCD15A2H，覆蓋范圍為海河流域，數(shù)據(jù)時間為2005至2018年（其中2005—2007的ac數(shù)據(jù)作為2008年ac時域重建的參考年數(shù)據(jù)），以上數(shù)據(jù)均可從EARTHDATA（https：//search.earthdata.nasa.gov/）下載，所有影像均通過MRT工具轉換為0.01°的分辨率。MCD43A3數(shù)據(jù)在時間上雖然為逐日序列，但在空間上覆蓋不完整，可通過時域重建方式得到時空連續(xù)的逐日數(shù)據(jù)產品。此外，本研究中土地利用數(shù)據(jù)采用LUCC（Land-Use and Land-Cover Change）地物分類標準。

為檢驗本研究得到的ET結果，本文還使用了MODIS的蒸散發(fā)產品（MOD16）。MOD16數(shù)據(jù)由美國國家航空航天局（NASA）發(fā)布，是目前應用最廣泛的遙感ET產品之一［21-22］。

3.2 實測數(shù)據(jù) 由于ET的站點監(jiān)測成本較高，實測資料缺乏，本研究采集了大興站2008年、館陶站和密云站2010年和望都站2018年的渦度相關數(shù)據(jù)，作為對PML模型模擬結果精度的驗證數(shù)據(jù)。研究所用氣象數(shù)據(jù)來自國家氣象站的地面監(jiān)測數(shù)據(jù)（www.data.cma.gov.cn），包括用來計算蒸散發(fā)的氣溫、濕度、日照時數(shù)、風速和氣壓等數(shù)據(jù)，數(shù)據(jù)序列為2008、2010和2018年三年的逐日觀測資料。將國家氣象站點的觀測數(shù)據(jù)基于IDW（Inverse Distance Weighted）方法插值成柵格數(shù)據(jù)，空間分辨率統(tǒng)一為0.01°。海河流域的主要土地利用分類及氣象站、渦度相關站點的位置如圖1所示，各渦度相關站點的數(shù)據(jù)情況如表1所示。

表1 海河流域渦度相關測站情況

圖1 海河流域土地利用及氣象站、渦度相關站點分布示意圖

4 結果及討論

4.1 海河流域ac重建結果對比 根據(jù)時域重建方法和線性插值法可得到的海河流域時空連續(xù)的ac產品。圖2展示了基于時域重建方法得到的2018年各季節(jié)海河流域ac產品。從圖中可以看出，2018年ac具有明顯的時空變異特征。在空間分布上，2018年ac整體上呈現(xiàn)出北部高，南部低的趨勢，同時城市區(qū)域的ac相對較低。在季節(jié)變化上，夏季和冬季ac的高值區(qū)范圍相對較小，主要集中在海河流域北部地區(qū)；而春季和秋季ac的空間分布相對均衡。在ac的最值方面，冬季的最大值為0.49，同時也是全年的最大值，秋季的最大值較小，為0.31。

圖2 時域重建法獲取的海河流域2018年各季節(jié)ac均值分布圖

與時域重建方法相比，基于線性插值方法得到的ac結果整體偏高（見圖3）。在2018年各季節(jié)中，時域重建的ac結果比線性插值結果偏低0.2以內的地區(qū)分別占到88.46%、92.60%、88.31%和82.23%，覆蓋了海河流域的大部分地區(qū)。在空間分布上，時域重建結果偏高的地區(qū)主要集中在海河流域北部及西部山區(qū)，其余地區(qū)時域重建的ac結果總體偏低。形成這種差異的原因在于兩種方法對數(shù)據(jù)缺失日期插補方法的不同，當某個網格的ac數(shù)據(jù)出現(xiàn)缺失時，線性插值法基于缺失日期前后的數(shù)據(jù)建立線性關系，插值結果不會出現(xiàn)明顯的低值，當有效數(shù)據(jù)在時間上的分布較為集中時，還可能會插值出異常高的ac值，樊輝［23］使用不同插值方法對NDVI（Normalized Difference Vegetation Index）數(shù)據(jù)的重建結果也表明，線性插值的重建結果偏差最大。而基于HANTS的時域重建方法在插補數(shù)據(jù)時是通過選取均值和顯著諧波，利用最小二乘擬合來迭代求解的過程，對缺失數(shù)據(jù)序列的插補會出現(xiàn)小于缺失日期前后數(shù)據(jù)的情況。

圖3 時域重建法與線性插值法獲取的海河流域各季節(jié)ac差值分布及統(tǒng)計圖

4.2 模型精度評價 基于PML模型對海河流域蒸散發(fā)進行反演，將采用時域重建方法得到ac作為輸入記為方案“PML_HANTS”，并利用大興站、館陶站和密云站的渦度相關數(shù)據(jù)對模擬結果進行驗證（表2）。對模型精度的評價指標選用相關系數(shù)（r）、均方根誤差（RMSE）和偏離度（Bias）［24］。從圖中可以看到，PML_HANTS方案的模擬結果與觀測值在三個站點具有較好的相關性，相關系數(shù)r分別為0.73、0.77和0.76；在RMSE和Bias指標方面，RMSE指標在2.3 mm/d以下，Bias指標為負值，表明本文的模擬值與觀測值相比偏低。與其他基于分布式的f和gsx參數(shù)化的模擬結果相比，三個站點的模擬結果的相關系數(shù)與Morillas等［25］的研究結果較為一致，表明本文基于PML_HANTS方案在海河流域的ET模擬結果與觀測值有較好的相關性。

表2 PML_HANTS方案下大興、館陶和密云三個站點ET模擬值與觀測值對比

為進一步分析PML_HANTS方案在流域尺度上的精度，選用MOD16的全球蒸散發(fā)產品與PML_HANTS方案的模擬結果進行對比。由于MOD16產品是8天合成的ET，在對比時將實測數(shù)據(jù)與模擬結果統(tǒng)一到與MOD16產品相同的時間尺度上。圖4給出了海河流域各渦度相關站點（見表1）8天合成的觀測數(shù)據(jù)與PML_HANTS方案模擬結果及MOD16產品的對比情況，從圖中可以看出，PML_HANTS方案模擬結果與渦度相關數(shù)據(jù)的r是0.68，RMSE是14.23，而MOD16產品與渦度相關數(shù)據(jù)的r和RMSE分別為0.45和15.09，本文提出的PML_HANTS方案模擬結果要優(yōu)于MOD16產品?？傮w來看，PML_HANTS方法在8天尺度上與觀測值的RMSE仍有一定偏差，其原因可能有兩方面：一是使用基于降雨數(shù)據(jù)的f參數(shù)化方法可能帶來不確定性，Morillas等［25］的研究也表明f的參數(shù)化方法對ET的模擬精度有一定影響；二是本文采用的氣象數(shù)據(jù)輸入是通過對氣象站點數(shù)據(jù)插值后得到，氣象數(shù)據(jù)的插值過程中也會對ET的模擬結果帶來一定誤差。

圖4 海河流域8天尺度渦度相關數(shù)據(jù)與PML_HANTS方案模擬結果及MOD16產品對比圖

4.3 ac對ET的影響分析 為分析本文兩種ac結果對ET的影響，還利用線性插值方法得到的ac對海河流域ET進行了模擬，記為方案“PML_LINEAR”，并對比在兩種ac輸入條件下，PML模型對望都站模擬結果的差異（表3）。從表2可以看出，兩種方案的模擬結果與觀測值均有較好的相關性，PML_LINEAR方案的r值為0.52，與PML_HANTS方案差異較小。在其他指標方面，PML_HANTS方案的RMSE比PML_LINEAR方案降低了0.18 mm/d，偏離度也減少了7.61%。因此，對望都站而言，PML_HANTS方案的模擬值與渦度相關觀測值的準確性更好。

表3 望都站PML_HANTS與PML_LINEAR方案模擬精度對比

在基于PML模型的ET模擬中，ac主要參與凈輻射的計算，進而影響到被冠層和土壤吸收的能量Ac和As，并與凈輻射的大小成反比。為更加清楚地反映PML_HANTS和PML_LINEAR方案中ac的差異對ET模擬結果的影響，以望都站2018年5月份（有觀測值日期）兩種方案下ac變化過程與ET變化過程為例進行對比分析（圖5）。從圖中可以看出，5月份兩種方案得到的ac值變化均較為平緩，但PML_HANTS方案下的ac值比PML_LINEAR方案整體偏低，平均值偏低約0.14。在ET的差異方面，兩種方案的模擬結果均較好地反映了ET的變化過程，但由于PML_HANTS方案的ac值相對較小，計算得到的凈輻射偏大，即冠層和土壤吸收到的輻射更多，PML_HANTS方案得到的ET結果平均偏高0.32 mm/d。由此可以看出，ac的計算方案的不同，影響了望都站ET的模擬精度。

圖5 2018年5月不同方案下ac及ET模擬值與觀測值對比圖

為進一步分析兩種ac重建方案對ET結果在海河流域范圍內的影響，將PML_HANTS方案的年蒸散發(fā)量與PML_LINEAR方案的結果進行了對比，兩種方案條件下2018年海河流域ET總量的差值及頻率統(tǒng)計如圖6所示。從圖中可以看出，PML_HANTS方案得到的年ET總量整體高于PML_LINEAR方案，年ET總量與PML_LINEAR方案的差值在0～200 mm的地區(qū)占到92.79%，其中海河流域中部和南部地區(qū)偏高較為明顯，望都站PML_HANTS方案的結果比PML_LINEAR方案的模擬結果偏高，是這種空間差異特征的具體體現(xiàn)。此外，還有7.15%的地區(qū)PML_HANTS方法偏低0～200 mm，主要分布在海河流域北部等地區(qū)。

圖6 海河流域PML_HANTS與PML_LINEAR模擬結果差值分布及統(tǒng)計圖

5 結論

本文基于PML模型對海河流域逐日蒸散發(fā)進行了反演研究，通過引入時域重建方法，獲取了時空連續(xù)的地表反照率數(shù)據(jù)，并與通過線性插值方法得到ac來計算ET的方法進行對比，探討了ac對ET模擬結果的影響，主要得到以下結論：

（1）通過與不同年份、4個渦度相關站點的監(jiān)測數(shù)據(jù)進行對比，本文提出的PML_HANTS方案模擬結果與渦度相關數(shù)據(jù)的相關性系數(shù)在0.53至0.77之間，相關性較好，但比觀測值偏低，需在后續(xù)工作中持續(xù)改進；同時在8天的時間尺度上，PML_HANTS方案模擬結果與渦度相關數(shù)據(jù)的相關性也優(yōu)于MOD16產品。

（2）對海河流域2018年ac的時域重建結果顯示，ac整體上呈現(xiàn)出北部高，南部低的分布特征，同時夏季和冬季ac的高值區(qū)范圍相對較小，而春季和秋季ac的空間分布相對均衡。

（3）對望都站2018年5月的分析表明，基于時域重建的ac結果整體上低于基于線性插值得到的ac數(shù)據(jù)，蒸散發(fā)模擬結果的RMSE降低了0.18 mm/d，偏離度降低了7%以上，可為進一步改善蒸散發(fā)精度提供思路。