集合卡爾曼濾波同化估算地表水熱通量的研究

劉翔舸，王鵬新，唐伯惠，黃健熙，陶 欣

(1. 中國農業大學信息與電氣工程學院 北京 海淀區 100083；2. 中國科學院地理科學與資源研究所 北京 朝陽區 100101；3. 美國馬里蘭大學帕克分校地理科學系 美國 馬里蘭帕克大學城 20742)

地表水熱通量在地面表層范圍內的大氣、水文和生態等諸多領域研究中具有重要的地位，其影響著地表水汽和能量的輸送過程。在農業大環境中，地表水熱通量的估計對于農業管理、作物估產以及農業生態系統的研究具有十分重要的意義。目前，在能量平衡原理、空氣動力學原理、土壤-植被-大氣傳輸理論的不斷發展下，已經建立了多種估算地表水熱通量的模型。如研究者利用建立相關的能量和物質交換過程方法的單層或多層模型，進行水熱通量的估算研究[1-2]。

另一方面，陸面數據同化方法的發展使數據同化技術在陸面過程研究中得到較多的應用[3]。其中，數據同化技術在實際應用中可同時利用多源觀測數據結合動態模型的多層次架構的特點，能夠將模型與觀測通過數據同化算法有機結合在一起，共同作用于對目標參數的估算反演中，從而獲得高分辨率的估算結果，使最終估算的目標數據質量得到提升。因此，數據同化技術在估算地表水熱通量的相關研究也不同程度地開展起來，文獻[4]應用變分同化方法，將遙感反演的地表溫度數據同化入陸面能量平衡與地表溫度為約束條件的動態方程組中，得到了潛熱和地表蒸散參數的有效估算。文獻[5]應用擴展卡爾曼濾波方法，并通過設計方案中考慮觀測和模型誤差的不確定性，將衛星反演溫度誤差與實測表層土壤溫度進行擬合，以擬合作為觀測誤差同化入地表陸面模型中，估計下層土壤溫度與地表蒸散。文獻[6]以能量平衡方程的方法提出了一種近地表陸面模式，并應用變分同化方法與遙感反演的地表溫度耦合，實現對地表熱通量的估算。文獻[7]應用弱約束的四維變分方法，對簡單陸面模式進行了數值試驗，將地面站實測地表溫度與陸面模式等關鍵參數共同建立代價函數，實現地表水熱通量的估計。文獻[8]采用通用陸面模式結合EnKF對地表水熱通量進行估算。文獻[9]綜合考慮了能量傳遞過程中土壤表層與植被冠層溫度之間的差異，并應用變分同化將地表溫度同化入一個雙層模型，估算了地表水熱通量數據。

由于多數地表水熱通量的研究工作是基于國外通量觀測站數據的支持，對國內研究的參考價值有限。因此本文嘗試利用國內試驗站點的觀測數據，通過搭建一個以同化地表溫度改進陸面過程模型中，地表水熱通量估計的數據同化系統，為數據同化方法在農田環境中的地表水熱通量估計的進一步研究提供借鑒與基礎。

1 研究站點概況與數據預處理

本文選用的試驗數據來自山東禹城農業綜合試驗站(36.83N ，116.57E )，如圖1所示，站點位于山東省禹城市，是中國科學院首批的開放式試驗站，并于1999年由國家科技部設立為首批國家重點試驗站。該站地理環境處于黃淮海平原，地處暖溫帶半濕潤季風氣候區。年平均氣溫13.1℃，年均降水量582 mm，光熱資源豐富，雨熱同期，是主要的農業生產區，站點農作物主要是小麥、玉米等，該站點同時具有農業生產的有利環境和農業作物類型的代表性。

禹城農業綜合試驗站通過不斷發展，陸續建立了多種類型儀器的監測與研究系統，可具備同時監測氣象數據、土壤數據、植物信息等能力。為農業的持續發展和陸面過程中的土壤-植物-大氣間的能量循環監測研究，提供完備的科技手段和觀測數據支持。

2 EnKF同化系統的構建

圖2所示為本文的研究數據同化系統的具體流程圖。從圖中可以大致看出，數據同化系統的主要結構可分為：外部觀測、同化算法、動態模型3個主要部分組成。數值試驗研究中的外部觀測數據與氣象驅動數據由上述國內禹城試驗站點的實測數據提供，下面將對數據同化系統的同化算法和陸面過程模型分別進行描述。

圖2 同化系統流程圖

2.1 EnKF同化方法

為了有效地應對非線性問題，從而合理得到誤差的估計，一種集合的卡爾曼濾波(ensemble Kalman filter，EnKF)發展起來[10-18]。EnKF運用Monte-Carlo方法與集合的思想，主要解決了一般同化方法在實際應用中背景誤差協方差矩陣估計和預報困難的問題；解決非線性系統的近似問題；避免使用變分模式下的伴隨模式；較易考慮模式誤差問題；并可以實現并行化計算[19]。因此，本文選擇EnKF方法作為數據同化系統中的優化算法。

EnKF算法作為應用最廣泛的非線性濾波同化方法，最初由文獻[12]提出并應用于海洋數據同化領域。但是隨著近年來陸面數據同化研究的不斷發展，EnKF方法也逐步應用于陸面過程研究中，其EnKF算法由KF算法推衍而來，主要結構仍分為預報階段和分析階段兩個組成部分不變，具體算法描述如下：將動態模型M在每個隨機變量第i時刻的模型的預報值xi作為背景場的狀態集合向量bX ，

2.2 動態模型描述

式中，Rn是地表凈輻射；H和LE分別是顯熱通量和潛熱通量；P是地表的熱慣量； 為每日頻率(1/86 400 s)；Tsi表示前一天第i時刻的地表溫度值；表示土壤溫度；DN表示一天內的時間步長。式(11)為該模式的主體部分，其中熱通量數據(潛熱、顯熱與凈輻射)依賴于地表溫度Ts的狀態估計。因此，準確估算與優化當前時刻的地表溫度Ts，對下一刻的地表水熱通量數據的估計至關重要。所以，可認為通過改變地表溫度的預報，進而可改變地表水熱通量的估計精度。模型中的參數具體求解與設定參考文獻[6-7, 20]。

2.3 同化系統中誤差的確定

在數據同化系統的研究中，由于系統中模型表述狀態變量的復雜性，在實際應用中同化系統的誤差估計與處理是受到廣泛關注的問題。數據同化的目的是在動態模型的動力框架內，通過綜合不同來源的不同時空分辨率的常規和非常規觀測，將動態模型和各種觀測信息集合成為不斷依靠觀測而調整優化模型的軌跡，減小模型預報誤差而獲得更高精度的預報值。因此，誤差的研究應始終貫穿于數據同化的研究之中[21]。

通常情況認為動態模型在運行一段時間后將自動達到一種穩定和平衡的狀態，此時的模型模擬結果在一定程度上具有模型自身誤差的代表性。本文根據該方法，利用運行一段時間后的模型模擬結果與地面試驗站點的實測值進行比較，可獲得統計上的模型誤差。

如圖3所示，利用禹城試驗站點提供的氣象數據使模型從第91天開始，持續運行90天的模擬結果。圖中所示從第100天左右開始，模型模擬的LST結果高于地表觀測的LST結果，之后在一段時間內(第100天～160天)的觀測值與模擬值的差別較大，此外，圖中顯示在第165天左右，試驗站點實測的地表溫度數值有部分超過320 K，認為該部分數據存在較大的測量誤差，在以后計算中以320 K為閾值，對觀測站點實測的地表溫度數據進行質量控制。

為了研究并得到模型的誤差，認為模型自身運行10天后的模擬結果即可達到穩定。將模型在一段時間內的(第100天～160天)對應試驗站點的觀測值與模型獨立運行的模擬值進行以天為單位的差值平均，得到該時段一天內的模型誤差變化趨勢，如圖4所示。圖中的結果表示，地表溫度的誤差主要體現在白天階段，并且白天的誤差相比夜間變化的波動較大。綜合考慮，模型地表溫度平均高于地面實測地表溫度約為1.5 K，本文在試驗中假設模型運行時的模型誤差為定值，不隨其他變量相關，通過上述統計方式得出平均每天模型誤差的均值，作為同化系統中模型誤差的初始設定。

圖3 山東禹城站點模擬地表溫度結果與觀測結果比較圖

圖4 一天內模擬值與觀測值的差值變化

3 結果分析與討論

3.1 集合大小的確定

EnKF集合的多少決定最終整個數據同化系統的效率和精度。集合數目過多影響系統執行效率，過少卻影響最終的同化精度。因此在本文試驗中，考慮利用調整集合m的大小，通過系統運行結果的均方根誤差(root mean square error，RMSE)變化，試確定適合本同化系統的集合數目m。

在確定集合大小的試驗中，同化步長，模型誤差、背景場誤差、模型驅動參數等均取值相同，只調整集合m的大小，將得到的同化后地表溫度結果與試驗站觀測的地表溫度對比，利用誤差RMSE變化分析EnKF的集合數目的影響，如圖5所示。

由圖5可以看出，RMSE變化趨勢隨著集合數目的增大顯著減小，當集合大小為30時RMSE最小，之后逐漸趨于穩定。因此取m=30。

3.2 同化結果

利用已設定的數據同化系統，對山東禹城2009年4～6月內進行EnKF的數據同化試驗。考慮數據的連續性和演示效果的清晰性，選擇試驗站點的第91天開始之后30天的數據同化結果進行展示。

圖5 集合大小對同化系統的影響

圖6所示為模型模擬的LST和同化后的LST的趨勢變化，兩者總體趨勢與實測的LST變化相同，但是從圖中可以明顯看出，同化后的LST趨勢更接近實測的LST曲線，說明相較于模型模擬的LST，同化后的LST精度得到了有效的改善。其中，模型模擬LST與地面實測LST的RMSE為4.85 K，經過同化后LST的RMSE則下降到1.57 K。證明在同化過程中對于LST估計精度的改善較為明顯。因此，伴隨在整個同化的過程中，改善的LST將使模型最終估計的地表水熱通量的精度同時得到提高。

圖6 山東禹城站點地表溫度、地表水熱通量的同化結果

另一方面，由于試驗中缺少對地面觀測通量數據的支持，為了驗證數據同化后的通量結果，利用MODIS的數據產品進行輔助驗證。將同化后得到的水熱通量中的潛熱通量與MODIS的ET產品(MOD16A2)進行對比分析。MODIS的ET產品分為8天和月的ET/LE的全球產品，分辨率為1 km′1 km。驗證試驗選擇MODIS對應禹城站點的8天ET產品與陸面模型獨自模擬的ET(未同化)和同化系統估算的ET進行對比驗證，驗證試驗將分別從ET變化趨勢和RMSE的結果兩方面做進一步的判斷與分析。

MOD16的ET/LE產品采用8天累加的ET值。文中的數據均為30 m ins，因此需要對數據進行時間尺度上的匹配。根據MODIS產品的定義，本文將模型模擬的ET與同化結果的ET，在8天內的數據累加生成8天的數據產品(mm/8 d)，與MOD16產品對比結果。

圖7所示為2009年第89天～177天(4～6月)間的8天的ET結果隨時間的變化的趨勢情況，圖中的3個方法得到的ET的總體變化趨勢相同。由圖中可以明顯看出，模型模擬得到的ET值總體偏高(與MOD16的ET之間的RMSE=4.18 mm)。

圖7 同化后ET結果對比圖

而經過同化后的ET整體趨勢與模型相比略低，但ET值仍高于MODIS的ET產品(與MOD16的ET之間的RMSE=2.99 mm)。綜合LST與ET結果的RMSE對比試驗，認為經過數據同化后的模型模擬結果(LST、H和LE)均得到了有效改善。因此，本文構建的利用地表溫度改進地表水熱通量估計數據同化系統具有可應用性與適用性。

4 結 論

本文利用山東禹城試驗站點數據通過EnKF同化方法，改進模型地表水熱通量估算精度的嘗試是成功的，具有較高的可重復性與穩定性，并具備一定的適用性。在EnKF的數據同化系統中，探討了系統中陸面模型誤差的設定問題，認為一定時間內模型誤差的變化可用統計的方法標定，進而運用于整個數據同化過程。同時，探討了EnKF算法中集合大小對同化結果的影響，通過對集合大小變化的誤差協方差RMSE的比較分析，確定已建立同化系統的集合數目。最后，利用該同化系統同化地表溫度數據，提高了模型的地表水熱通量的最終估計精度，認為該同化系統方案可有效地改善陸面過程模型對地表水熱通量的估計。該研究方法和結果在農田地表水熱通量的估算方面具有參考價值。

感謝中國科學院地理科學與資源研究所提供山東禹城試驗站觀測數據。

[1] 劉紹民, 李小文, 施生錦, 等. 大尺度地表水熱通量的觀測、分析與應用[J]. 地球科學進展, 2010, 25(11): 1113-1127.LIU Shao-m in, LI Xiao-wen, SHI Sheng-jin, et al.Measurement, analysis and application of surface energy and water vapor fluxes at large scale[J]. Advances in Earth Science, 2010, 25(11): 1113-1127.

[2] 辛曉洲, 田國良, 柳欽火. 地表蒸散定量遙感的研究進展[J]. 遙感學報, 2003(3): 233-240.

XIN Xiao-zhou, TIAN Guo-liang, LIU Qin-huo. A review of researches on remote sensing of land surface evapotranspiration[J]. Journal of Remote Sensing, 2003(3):233-240.

[3] LIANG S. Advances in land remote sensing, system,modeling, inversion and application[M]. New York, USA:Springer, 2008.

[4] BONI G, ENTEKHABI D, CASTELLI F. Land data assimilation w ith satellite measurements for the estimation of surface energy balance components and surface control on evaporation[J]. Water Resources Research, 2001, 37(6):1713-1722.

[5] KUMAR P, KALEITA A L. Assimilation of near-surface temperature using extended Kalman filter[J]. Advances in Water Resources, 2003, 26(1): 79-93.

[6] CAPARRINI F, CASTELLI F, ENTEKHABI D. Estimation of surface turbulent fluxes through assimilation of radiometric surface temperature sequences[J]. Journal of Hydrometeorology, 2004, 5(1): 145-159.

[7] JUN Q, LIANG S L, LIU R G, et al. A weakconstraint-based data assimilation scheme for estimating surface turbulent fluxes[J]. Geoscience and Remote Sensing Letters, 2007, 4(4): 649-653.

[8] XU T R, LIANG S L, LIU S M. Estimating turbulent fluxes through assimilation of geostationary operational environmental satellites data using ensemble Kalman filter[J]. Journal of Geophysical Research-Atmospheres,2011: 116(D09): 109-125.

[9] BATENI S M, LIANG S. Estimating surface energy fluxes using a dual-source data assimilation approach adjoined to the heat diffusion equation[J]. Journal of Geophysical Research-Atmospheres, 2012, 117(D17): 118-132.

[10] EVENSEN G. Sequential data assimilation w ith a nonlinear quasi-geostrophic model using monte-carlo methods to forecast error statistics[J]. Journal of Geophysical Research-Oceans, 1994, 99(C5): 10143-10162.

[11] EVENSEN G. Advanced data assimilation for strongly nonlinear dynamics[J]. Monthly Weather Review, 1997,125(6): 1342-1354.

[12] EVENSEN G. The ensemble Kalman filter: theoretical formulation and practical implementation[J]. Ocean Dynamics, 2003, 53(4): 343-367.

[13] BURGERS G, VAN LEEUWEN P J, EVENSEN G.Analysis scheme in the ensemble Kalman filter[J].1998(126): 1719-1724.

[14] BISHOP C H, TOTH Z. Ensemble transformation and adaptive observations[J]. Journal of the Atmospheric Sciences, 1999, 56(11): 1748-1765.

[15] BOCQUET M. Ensemble Kalman filtering w ithout the intrinsic need for inflation[J]. Nonlinear Processes in Geophysics, 2011, 18(5): 735-750.

[16] WANG X G, BISHOP C H, JULIER S J. Which is better,an ensemble of positive-negative pairs or a centered spherical simplex ensemble?[J]. Monthly Weather Review,2004, 132(7): 1590-1605.

[17] WHITAKER J S, HAM ILL T M. Ensemble data assim ilation w ithout perturbed observations[J]. Monthly Weather Review, 2006, 134(6): 1722-1722.

[18] WANG X G, HAM ILL T M, WHITAKER J S, et al. A comparison of the hybrid and EnSRF analysis schemes in the presence of model errors due to unresolved scales[J].Monthly Weather Review, 2009, 137(10): 3219-3232.

[19] 黃健熙, 武思杰, 劉興權, 等. 基于遙感信息與作物模型集合卡爾曼濾波同化的區域冬小麥產量預測[J]. 農業工程學報, 2012, 28(4): 142-148.

HUANG Jian-xi, WU Si-jie, LIU Xing-quan, et al.Regional w inter wheat yield forecasting based on assim ilation of remote sensing data and crop grow th model w ith ensemble Kalman method[J]. Transactions of the Chinese Society of Agricultural Engineering, 2012, 28(4):142-148.

[20] CAPARRINI F, CASTELLI F, ENTEKHAB D. Variational estimation of soil and vegetation turbulent transfer and heat flux parameters from sequences of multisensor imagery[J].Water Resources Research, 2004, 40(12): 5-20.

[21] 擺玉龍, 李新, 韓旭軍. 陸面數據同化系統誤差問題研究綜述[J]. 地球科學進展, 2011, 26(8): 795-804.

BAI Yu-long, LI Xin, HAN Xu-jun. A review of error problems for land data assim ilation systems[J]. Advances in Earth Science, 2011, 26(8): 795-804.

編 輯 黃 莘