集合卡爾曼濾波數(shù)據(jù)同化方法改進(jìn)土壤水分模擬效果

陳　鶴，楊大文，劉　鈺，張寶忠（. 中國(guó)水利水電科學(xué)研究院，流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 00048；. 清華大學(xué)水利水電工程系，水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 00084）

?

集合卡爾曼濾波數(shù)據(jù)同化方法改進(jìn)土壤水分模擬效果

陳鶴1，楊大文2，劉鈺1，張寶忠1
（1. 中國(guó)水利水電科學(xué)研究院，流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100048；2. 清華大學(xué)水利水電工程系，水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084）

摘要：陸面過程模型是連續(xù)模擬土壤水分的有效工具，然而輸入數(shù)據(jù)及模型結(jié)構(gòu)本身的不確定性會(huì)導(dǎo)致模擬誤差在模型運(yùn)行過程中不斷積累。數(shù)據(jù)同化技術(shù)可以考慮模型不確定性，實(shí)時(shí)修正模型狀態(tài)變量，進(jìn)而提高土壤水分的模擬精度。本研究構(gòu)建集合卡爾曼濾波（EnKF， ensemble Kalman filter）數(shù)據(jù)同化方法，將其集成到水文強(qiáng)化陸面過程模型HELP （hydrologically-enhanced land process）中，對(duì)模型中土壤水分及表面溫度等狀態(tài)變量進(jìn)行優(yōu)化。模型選取山東位山生態(tài)水文觀測(cè)站2006年的數(shù)據(jù)進(jìn)行驗(yàn)證，采用未經(jīng)同化的模型率定結(jié)果作為基準(zhǔn)值。結(jié)果表明，數(shù)據(jù)同化后表層、根層、深層土壤水分模擬結(jié)果相比基準(zhǔn)值均有提高，土壤含水量均方根誤差減小30%～50%，證明采用數(shù)據(jù)同化方法能夠有效提高土壤水分的模擬結(jié)果。

關(guān)鍵詞：土壤；遙感；溫度；數(shù)據(jù)同化；陸面過程模型；土壤含水量

陳鶴，楊大文，劉鈺，張寶忠. 集合卡爾曼濾波數(shù)據(jù)同化方法改進(jìn)土壤水分模擬效果[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2016，32（2）：99－104.doi：10.11975/j.issn.1002-6819.2016.02.015http://www.tcsae.org

Chen He， Yang Dawen， Liu Yu， Zhang Baozhong. Data assimilation technique based on ensemble Kalman filter for improving soil water content estimation[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE），2016， 32（2）: 99－104. （in Chinese with English abstract）doi：10.11975/j.issn.1002-6819.2016.02.015http://www.tcsae.org

0　引　言

非飽和帶的土壤含水量是聯(lián)結(jié)地表與地下水分與能量交換的紐帶，影響著降水入滲過程、地表蒸散發(fā)等重要水文環(huán)節(jié)[1]，對(duì)于氣象學(xué)、水文學(xué)及農(nóng)學(xué)都具有重要意義。在農(nóng)業(yè)研究中，準(zhǔn)確模擬土壤含水量是農(nóng)業(yè)水管理、灌溉制度確定，以及農(nóng)業(yè)增產(chǎn)的必要前提[2]。土壤含水量可以通過地面觀測(cè)、模型模擬以及遙感反演定量獲取。然而，土壤含水量空間變異性較大，點(diǎn)尺度地面觀測(cè)難以獲取大范圍土壤含水量空間分布情況[3]；而遙感數(shù)據(jù)雖然能獲取空間信息，但僅能反演衛(wèi)星過境時(shí)刻的瞬時(shí)狀態(tài)，在時(shí)間上是不連續(xù)的[4]；陸面過程模型能夠連續(xù)模擬土壤含水量空間分布，但一方面模擬結(jié)果對(duì)初值誤差非常敏感，另一方面模擬誤差在模型運(yùn)行過程中不斷積累，導(dǎo)致模擬精度無法保證[5]。數(shù)據(jù)同化技術(shù)能夠考慮模型結(jié)構(gòu)及模型輸入數(shù)據(jù)的不確定性，將其與遙感信息結(jié)合，引入遙感觀測(cè)數(shù)據(jù)，在衛(wèi)星過境時(shí)刻實(shí)時(shí)修正模型狀態(tài)變量，可以充分發(fā)揮陸面過程模型和遙感數(shù)據(jù)的優(yōu)勢(shì)，為土壤含水量的連續(xù)模擬提供新途徑[6]。數(shù)據(jù)同化的核心思想是誤差估算及誤差模擬，即在動(dòng)力學(xué)模型框架下，融合不同來源、不同時(shí)空分辨率、不同精度的觀測(cè)數(shù)據(jù)，根據(jù)不同觀測(cè)之間的誤差關(guān)系，通過數(shù)學(xué)算法對(duì)模型中的狀態(tài)變量進(jìn)行優(yōu)化，以期提高模擬精度[7]。數(shù)據(jù)同化的概念和算法最早由海洋學(xué)家和氣象學(xué)家提出[8]，近年來逐漸引入水文學(xué)研究中[9]，例如土壤水同化[7]，地表溫度同化[10]等。隨著遙感技術(shù)的發(fā)展，越來越多的地表參數(shù)可以由遙感觀測(cè)反演獲取，給同化系統(tǒng)提供了數(shù)據(jù)基礎(chǔ)[11]。本研究構(gòu)建集合卡爾曼濾波（EnKF， ensemble Kalman filter）數(shù)據(jù)同化方法，將其集成到水文強(qiáng)化陸面過程模型HELP （hydrologically-enhanced land process）中，采用SEBS （surface energy balance system）遙感蒸散發(fā)模型模擬的衛(wèi)星過境時(shí)刻瞬時(shí)潛熱通量為觀測(cè)值，對(duì)陸面過程模型中土壤水分及表面溫度進(jìn)行同化，驗(yàn)證數(shù)據(jù)同化方法在改進(jìn)土壤含水量模擬精度中的作用。

1　數(shù)據(jù)與方法

1.1研究區(qū)域及數(shù)據(jù)

驗(yàn)證站點(diǎn)選取山東位山生態(tài)水文觀測(cè)站（36°38′55.5″N，116°3′15.3″E，海拔高度30 m（黃海海平面））。觀測(cè)站位于山東位山灌區(qū)（36°8′～37°1′N，115°25′～116°31′E）中部一處典型農(nóng)田內(nèi)，在區(qū)域內(nèi)具有代表性。觀測(cè)站內(nèi)設(shè)置有高10 m的鐵塔，在鐵塔上及鐵塔周邊的兩個(gè)土壤剖面中安置觀測(cè)儀器，對(duì)農(nóng)田內(nèi)的水、熱、碳通量以及常規(guī)氣象及土壤、作物狀態(tài)進(jìn)行觀測(cè)。與本研究?jī)?nèi)容相關(guān)的觀測(cè)包括通量觀測(cè)、氣象觀測(cè)及植被參數(shù)觀測(cè)。

1）通量觀測(cè)：包括土壤熱通量，以及由渦度相關(guān)系統(tǒng)觀測(cè)的顯熱通量和潛熱通量，以及由大孔徑激光閃爍儀觀測(cè)的顯熱通量；

2）氣象觀測(cè)：包括降雨、日照時(shí)數(shù)、空氣溫度及濕度、氣壓、風(fēng)速及風(fēng)向、向上/向下長(zhǎng)、短波輻射，以及熱紅外地表溫度；

3）植被參數(shù)觀測(cè)：包括冠層高度及葉面積指數(shù)，其中冠層高度采用實(shí)地采樣的方式每2周觀測(cè)一次，葉面積指數(shù)采用在農(nóng)田內(nèi)隨機(jī)采樣的方式進(jìn)行觀測(cè)。

遙感數(shù)據(jù)來源于Terra和Aqua衛(wèi)星上搭載的MODIS傳感器觀測(cè)的陸面產(chǎn)品，其中Terra衛(wèi)星的過境時(shí)間為上午10:30左右，Aqua衛(wèi)星的過境時(shí)間為下午1:30左右。MODIS陸面產(chǎn)品的空間分辨率為250 m到1 km，時(shí)間分辨率為逐日到16d合成不等。以上數(shù)據(jù)產(chǎn)品的介紹可參見MODIS網(wǎng)站（http://modis-land.gsfc.nasa.gov/），遙感數(shù)據(jù)可以通過NASA數(shù)據(jù)平臺(tái)下載（http://reverb.echo.nasa.gov/reverb/）。下載的數(shù)據(jù)產(chǎn)品經(jīng)過坐標(biāo)轉(zhuǎn)換、重采樣、質(zhì)量控制、數(shù)據(jù)插補(bǔ)等預(yù)處理環(huán)節(jié)，得到最終的地表輸入數(shù)據(jù)集。

1.2數(shù)據(jù)同化系統(tǒng)

數(shù)據(jù)同化系統(tǒng)的核心組成部分是模型算子、觀測(cè)算子和同化算法。在陸面數(shù)據(jù)同化系統(tǒng)中，模型算子通常采用水文模型或者陸面過程模型，用來模擬地表水熱耦合的物理過程；觀測(cè)算子用來連接需要被優(yōu)化的模型狀態(tài)變量和用來輔助同化的觀測(cè)數(shù)據(jù)；數(shù)據(jù)同化算法通常分為變分算法和濾波算法兩種，陸面數(shù)據(jù)同化研究中通常采用濾波算法。陸面數(shù)據(jù)同化系統(tǒng)中用到的數(shù)據(jù)，除了隨時(shí)間變化的模型輸入數(shù)據(jù)和不隨時(shí)間變化的模型參數(shù)集之外，還需要一套獨(dú)立來源的觀測(cè)數(shù)據(jù)用于同化計(jì)算。本研究中用到的模型算子為水文強(qiáng)化陸面過程模型HELP，數(shù)據(jù)同化算法為集合卡爾曼濾波方法（EnKF），用于同化運(yùn)算的觀測(cè)數(shù)據(jù)為由遙感蒸散發(fā)模型SEBS計(jì)算得到的瞬時(shí)潛熱通量。

1.2.1水文強(qiáng)化陸面過程模型

水文強(qiáng)化陸面過程模型（HELP）以SiB2模型（simple biosphere model 2）為基礎(chǔ)，并強(qiáng)化了產(chǎn)流過程等水文過程模擬。HELP模型中，陸氣間輻射傳輸，能量以及碳交換的模擬沿用SiB2模型，采用水熱傳輸模型和阻抗網(wǎng)絡(luò)模擬水熱狀況及土壤-植被-大氣間的水熱傳輸[12]。對(duì)于SiB2模型的水文模擬部分，HELP模型做出了以下改進(jìn)：1）在包氣帶土壤水分模擬中，HELP模型將包氣帶細(xì)化為厚度為0.1 m的土層，采用一維Richards方程描述土層之間的土壤水分交換；2）在計(jì)算過程中，HELP模型采用van Genuchten公式[13]取代SiB2模型中的Campbell/ Clapp-Hornberger公式，描述非飽和導(dǎo)水率與飽和導(dǎo)水率之間的關(guān)系；3）此外，HELP模型在SiB2地表產(chǎn)流模型的基礎(chǔ)上，增加了壤中流和地下水出流的模擬，采用質(zhì)量守恒方程和達(dá)西定律描述土壤水與河道之間的水量交換[14]。

為了精確模擬土壤含水量，進(jìn)一步對(duì)HELP模型土壤層離散方法加以細(xì)化，將土壤水分剖面細(xì)化為10層。因此HELP模型中共有17個(gè)狀態(tài)變量，分別為：10層土壤含水量，植被冠層和地表層對(duì)降雨的截留及貯存量，冠層、地表層及土壤深層的溫度，地下水位以及河道水位。在分析模擬結(jié)果時(shí)，再將10層土壤含水量換算為表層、根層、深層土壤含水量。

HELP模型的輸入數(shù)據(jù)主要為常規(guī)氣象數(shù)據(jù)，包括向下短波輻射、空氣溫度、相對(duì)濕度、風(fēng)速、降水量及灌溉量，以上數(shù)據(jù)均來自于位山站觀測(cè)數(shù)據(jù)。除了輸入數(shù)據(jù)外，HELP模型采用一系列參數(shù)描述土壤的物理特性及植被的生理特性。植被參數(shù)由植被類型決定，不隨時(shí)間變化，包括植被形態(tài)以及植被光學(xué)和生理學(xué)參數(shù)。土壤參數(shù)的取值見表1，植被參數(shù)的取值來源于位山站試驗(yàn)數(shù)據(jù)。除了兩個(gè)表中列出的參數(shù)外，模型的其他參數(shù)與SiB2模型中的默認(rèn)參數(shù)一致，未進(jìn)行人工調(diào)參。

HELP模型的模擬時(shí)段為2006年1月1日至2006年12月31日，模擬步長(zhǎng)為半小時(shí)。模型的初始條件（10層土壤含水量，冠層、地表及深層土壤溫度）取自實(shí)測(cè)值。模型的輸出值包括地表能量通量，土壤含水量，植被冠層、地表層、土壤深層3層溫度，徑流量，地下水位，河道水位以及碳通量方面的計(jì)算結(jié)果。本研究中對(duì)土壤含水量，冠層溫度，凈輻射，土壤熱通量，顯熱通量以及潛熱通量進(jìn)行驗(yàn)證，驗(yàn)證數(shù)據(jù)均來自位山站觀測(cè)值。其中，10層土壤含水量重新插值成表層、根層、深層含水量進(jìn)行驗(yàn)證。

表1　位山站土壤水分特征參數(shù)Table 1　Soil water characteristic parameter of Weishan Station

1.2.2集合卡爾曼濾波方法

集合卡爾曼濾波方法（EnKF）方法采用蒙特卡洛隨機(jī)采樣，依靠隨機(jī)生成一系列樣本，通過計(jì)算樣本的隨機(jī)誤差，直接得到模型誤差分布。采用EnKF方法進(jìn)行數(shù)據(jù)同化分為兩個(gè)步驟：1）集合預(yù)報(bào)；2）狀態(tài)變量校正。在每一個(gè)時(shí)刻t，采用蒙特卡洛方法對(duì)模型的狀態(tài)變量隨機(jī)生成一系列樣本，每一個(gè)樣本寫作xbi,t，其中下標(biāo)i表示第i個(gè)樣本，上標(biāo)b表示同化前的狀態(tài)變量。向量x的維度是m，表示模型中有m個(gè)狀態(tài)變量。在初始時(shí)刻，模型的集合預(yù)報(bào)可以用公式表示為

其中xi，0為模型的初始變量；M為模型算子。

同化前的狀態(tài)變量樣本矩陣為Xb

其中n是樣本數(shù)，則Xb的維度是m行n列。數(shù)據(jù)同化前狀態(tài)變量的平均值可以表示為

對(duì)同化前狀態(tài)變量樣本矩陣進(jìn)行背景場(chǎng)誤差計(jì)算如下

其中Pb是m×m維的矩陣，X′b是狀態(tài)變量離均值樣本矩陣，維度是m×n。

采用濾波法對(duì)狀態(tài)變量進(jìn)行校正，計(jì)算公式如下：

式（5）中的K是卡爾曼增益值（Kalman gain），計(jì)算公式為

其中R是作為同化量的觀測(cè)值背景場(chǎng)誤差，維度為p×p。

根據(jù)式（5）對(duì)每一個(gè)樣本的狀態(tài)變量xbi進(jìn)行校正，得到同化后的狀態(tài)變量xɑi，再進(jìn)行下一個(gè)時(shí)間步長(zhǎng)的模型預(yù)測(cè)

根據(jù)上述EnKF算法原理及計(jì)算公式，即可采用EnKF方法對(duì)模型進(jìn)行同化。首先隨機(jī)生成一組狀態(tài)變量的集合樣本，狀態(tài)變量的背景場(chǎng)誤差由下式計(jì)算

同樣，對(duì)觀測(cè)變量也加入隨機(jī)誤差，生成一組集合樣本

其中y是實(shí)測(cè)值，yi是加入了隨機(jī)誤差ηi之后的觀測(cè)值樣本。同化的具體過程可按照下式計(jì)算

式中I是單位矩陣。

1.2.3遙感蒸散發(fā)模型

采用基于能量平衡原理的SEBS（surface energy balance system）模型[15]模擬的衛(wèi)星過境時(shí)刻瞬時(shí)潛熱通量作為同化系統(tǒng)的觀測(cè)值。SEBS模型包含以下4個(gè)模塊：1）基于遙感空間反照率和輻射率的地表物理參數(shù)反演；2）熱量粗糙長(zhǎng)度的計(jì)算；3）顯熱通量的計(jì)算；4）潛熱通量的計(jì)算。與其他基于能量平衡原理的單層模型相比，SEBS模型的優(yōu)點(diǎn)在于每一個(gè)網(wǎng)格都獨(dú)立計(jì)算，因此即使在某天某些網(wǎng)格因?yàn)殛幱昊蛟频扔绊懭狈b感數(shù)據(jù)，也并不影響其他網(wǎng)格的計(jì)算結(jié)果，可以最大化地利用遙感數(shù)據(jù)。SEBS模型在位山站驗(yàn)證的結(jié)果表明，小麥季和玉米季的模擬均方根誤差均小于20%[16]，證明SEBS模型在研究區(qū)具有良好的模擬精度。

1.2.4數(shù)據(jù)同化方案

在數(shù)據(jù)同化系統(tǒng)中，需要計(jì)算初始背景場(chǎng)誤差協(xié)方差矩陣和模型誤差協(xié)方差矩陣。本研究對(duì)模型輸入數(shù)據(jù)加入隨機(jī)擾動(dòng)，生成一系列隨機(jī)樣本集合，通過模型自由運(yùn)行得到模型狀態(tài)變量的隨機(jī)樣本，進(jìn)而計(jì)算背景誤差。分別對(duì)HELP模型的輸入數(shù)據(jù)（降雨、風(fēng)速、相對(duì)濕度、空氣溫度、向下短波輻射）進(jìn)行隨機(jī)采樣，參考已有的研究成果，本研究采取的輸入數(shù)據(jù)隨機(jī)誤差見表2。

表2　HELP模型輸入數(shù)據(jù)隨機(jī)擾動(dòng)Table 2　Random disturbation of model input in HELP model

由于HELP模型的模擬步長(zhǎng)為半小時(shí)，SEBS模型反演的潛熱通量為每天2次，因此在有潛熱通量模擬值時(shí)進(jìn)行數(shù)據(jù)同化，其他計(jì)算時(shí)段模型自由運(yùn)行，給定觀測(cè)值的隨機(jī)誤差為20%。另外，在HELP模型的17個(gè)狀態(tài)變量中，冠層和地表截留是隨降雨輸入變化的狀態(tài)變量，在時(shí)間上不存在連續(xù)變化規(guī)律，即誤差不隨時(shí)間傳播。而地下水位與河道水位還受到臨近網(wǎng)格的影響，在一維垂向模型處理中無法準(zhǔn)確模擬。因此僅對(duì)10層土壤含水量及3層溫度進(jìn)行優(yōu)化計(jì)算，狀態(tài)變量向量的維度為13。選擇樣本數(shù)為10倍狀態(tài)變量向量維度，即130。

2　結(jié)果與分析

2.1模型自由運(yùn)行結(jié)果

在評(píng)價(jià)數(shù)據(jù)同化計(jì)算效果時(shí)，需要選定作為參照的基準(zhǔn)值，即模型未經(jīng)數(shù)據(jù)同化的自由運(yùn)行的模擬結(jié)果。模型自由運(yùn)行與數(shù)據(jù)同化方案采用同一套輸入數(shù)據(jù)、參數(shù)集以及初始狀態(tài)變量，模擬的3層土壤含水量（soil water content， SWC）如圖1所示。

表層、根層、深層土壤含水量的模擬相對(duì)誤差分別3.3%、4.7%和?0.9%，盡管整體的相對(duì)誤差較小，但模擬值與實(shí)測(cè)值吻合程度較低，3層土壤含水量絕對(duì)平均誤差分別為17.1%、16.3%和14.7%，均方根誤差分別為0.055、0.053、0.053 m3/m3。從圖中可以看出，在降雨/灌溉前后，模擬值的變化非常劇烈，表明在HELP模型中，土壤含水量對(duì)降雨/灌溉的輸入非常敏感。其中表層土壤含水量的變化幅度最為劇烈，根層其次，即使深層土壤含水量也會(huì)隨著降雨過程發(fā)生比較明顯的升高。

圖1　HELP模型自由運(yùn)行3層土壤含水量模擬值與實(shí)測(cè)值比較Fig.1　Open loop run of estimated and observed three-layer soil water content

2.2數(shù)據(jù)同化后模擬結(jié)果與基準(zhǔn)值對(duì)比

對(duì)同化結(jié)果進(jìn)行驗(yàn)證，模型從2006年1月1日起運(yùn)行至2006年12月31日，分別截取小麥季（3月1日－5 月31日）和玉米季（7月1日－9月30日）的結(jié)果進(jìn)行分析。采用百分比相對(duì)誤差（biɑs），百分比絕對(duì)平均誤差（mɑe），均方根誤差（rmse），以及數(shù)據(jù)同化效率系數(shù)（Eff）對(duì)同化結(jié)果進(jìn)行驗(yàn)證，其中Eff的計(jì)算公式為

式中上標(biāo)u和b分別表示同化后和同化前；Oi為第i時(shí)刻的觀測(cè)值；Sij為第i時(shí)刻的模擬值；i1和i2為模擬起止時(shí)間。Eff指標(biāo)大于0，表示數(shù)據(jù)同化后，模擬結(jié)果有改善；反之，小于0，則表示經(jīng)過數(shù)據(jù)同化后，模擬結(jié)果反而不如同化前。Eff值越趨近于100%，表示同化效果越好。

2.2.1表面溫度同化結(jié)果

圖2是小麥季和玉米季同化前后的表面溫度模擬值與實(shí)測(cè)值對(duì)比散點(diǎn)圖。小麥季表面溫度的模擬基準(zhǔn)值存在較大的系統(tǒng)偏差，數(shù)據(jù)同化后模擬值比基準(zhǔn)值系統(tǒng)偏差明顯降低，相對(duì)誤差減小至?16.6%，絕對(duì)平均誤差為20.9%，均方根誤差為1.79℃，數(shù)據(jù)同化的效率系數(shù)達(dá)到94.8%，確定性系數(shù)由0.66提高到0.95。玉米季表面溫度模擬基準(zhǔn)值模擬較好，均方根誤差為2.47℃，確定性系數(shù)為0.94。數(shù)據(jù)同化后，各項(xiàng)指標(biāo)也有了顯著提升，相對(duì)誤差由6.7%減小至3.2%，均方根誤差降低至1.29℃，確定性系數(shù)提高到0.97，數(shù)據(jù)同化效率系數(shù)為73.0%。結(jié)果表明數(shù)據(jù)同化對(duì)于改進(jìn)表面溫度狀態(tài)變量的模擬具有明顯效果。

圖2　表面溫度模擬值與實(shí)測(cè)值比較Fig.2　Scatter plot of estimated and observed surface temperature

2.2.2土壤含水量同化結(jié)果

圖3和4分別是小麥季和玉米季3層土壤含水量觀測(cè)值、基準(zhǔn)值和同化后的模擬值對(duì)比圖，表3是同化后3層土壤含水量與基準(zhǔn)值誤差對(duì)比分析。統(tǒng)計(jì)結(jié)果表明小麥季和玉米季同化后的模擬結(jié)果對(duì)比基準(zhǔn)值都有了不同程度的改善，均方根誤差降低，Eff的范圍從9.31%到74.17%。

從誤差評(píng)價(jià)指標(biāo)來看，在小麥季表層土壤水同化方面，同化以后的模擬值相對(duì)誤差和絕對(duì)平均誤差比基準(zhǔn)值略有增高。原因在于小麥季表層土壤水的基準(zhǔn)模擬值在降雨/灌溉前后變化非常劇烈，在降雨/灌溉發(fā)生后明顯升高，而過了一段時(shí)間后又顯著降低，導(dǎo)致在整個(gè)模擬時(shí)段內(nèi)，相對(duì)誤差處在較低的水平。但從圖中分析同化結(jié)果的具體過程可以看出，經(jīng)過數(shù)據(jù)同化后的表層土壤含水量在時(shí)間變化過程上與實(shí)測(cè)值吻合良好，但是存在一定的系統(tǒng)偏低。同化后的模擬均方根誤差比基準(zhǔn)值降低0.01 m3/m3，而同化后小麥季表層土壤含水量的效率系數(shù)為9.31%，表明雖然改善程度較根層和深層土壤含水量為低，但經(jīng)過數(shù)據(jù)同化后，表層土壤含水量的模擬值仍有一定程度的改善。小麥季根層和深層土壤含水量的同化值與基準(zhǔn)值相比有了明顯的改善，系統(tǒng)誤差顯著減小，與實(shí)測(cè)值的吻合程度也更高。

圖3　小麥季土壤含水量模擬值與實(shí)測(cè)值比較Fig.3　Soil water content estimation and observation in wheat season

圖4　玉米季土壤含水量模擬值與實(shí)測(cè)值比較Fig.4　Soil water content estimation and observation in maize season

與小麥季相比，玉米季土壤含水量的同化效果更好，尤其在相對(duì)誤差方面，同化后的3層土壤含水量模擬相對(duì)誤差分別為3.70%，3.62%和?5.45%，比基準(zhǔn)值有了顯著改善，絕對(duì)平均誤差、均方根誤差也顯著降低，效率系數(shù)均在60%以上，表3所示。

表3　同化后土壤含水量與基準(zhǔn)值誤差對(duì)比分析Table 3　Error analysis of updated soil water content withbenchmark run

3　結(jié)　論

將數(shù)據(jù)同化算法結(jié)合遙感數(shù)據(jù)集成到陸面過程模型中，為土壤含水量的連續(xù)模擬提供了新途徑。本研究將集合卡爾曼濾波算法集成到陸面過程模型HELP中，構(gòu)建陸面數(shù)據(jù)同化系統(tǒng)。采用對(duì)輸入數(shù)據(jù)加入隨機(jī)擾動(dòng)的方法生成隨機(jī)樣本和背景場(chǎng)誤差。利用遙感反演的潛熱通量模擬值作為同化系統(tǒng)的觀測(cè)值，對(duì)HELP模型狀態(tài)變量進(jìn)行優(yōu)化，并分別驗(yàn)證了同化結(jié)果對(duì)表面溫度和土壤含水量模擬的改進(jìn)程度。結(jié)果表明數(shù)據(jù)同化后的土壤含水量模擬結(jié)果明顯優(yōu)于未經(jīng)過同化的模擬基準(zhǔn)值，模擬相對(duì)誤差、絕對(duì)平均誤差以及均方根誤差都有不同程度的降低，在整個(gè)生育期內(nèi)的過程模擬與實(shí)測(cè)值吻合程度更高。綜合來看，同化后的土壤含水量均方根誤差降低幅度在30%～50%之間，表明對(duì)遙感反演的潛熱通量進(jìn)行同化，能夠改善HELP模型土壤水狀態(tài)變量模擬精度。

[參考文獻(xiàn)]

[1] Entekhabi D， Nakamura H， Njoku E G. Solving the inverse problem for soil moisture and temperature profiles by sequential assimilation of multifrequency remotely sensed observations[J]. IEEE Transactions on Geoscience and Remote Sensing， 1994， 32（2）: 438－448.

[2] Jackson T J， Hawley M E， O’Neill P.E. Preplanting soil moisture using passive microwave sensors[J]. Water Resources Bulletin， 1987， 23（1）: 11－19.

[3] Jackson T J， Schmugge T J， Nicks A D， et al. Soil moisture updating and microwave remote sensing for Hydrological simulation[J]. Hydrological Sciences Bulletin， 1981， 26（3）: 305－319.

[4] Lin S D， Wood E F， Troch Peter A， et al. Comparisons of remotely sensed and model-simulated soil moisture over a heterogeneous watershed[J]. Remote Sensing of Environment，1994， 48: 159－71.

[5] Feddes R A， Menenti M， Kabat P， et al. Is large-scale inverse modelling of unsaturated flow with areal average evaporationand surface soil moisture as estimated from remote sensing feasible[J]. Journal of Hydrology， 1993， 143: 125－52.

[6] Daley R. Atmospheric data analysis[M]. New York: Cambridge University Press， 1991， 457.

[7] Margulis S A， Mclaughlin D， Entekhabi D， et al. Land data assimilation and estimation of soil moisture using measurements from the Southern Great Plains 1997 field experiment[J]. Water Resources Research， 2002， 38（12）: 1－18.

[8] Salamon P， Feyen L. Assessing parameter， precipitation， and predictive uncertainty in a distributed hydrological model using sequential data assimilation with the particle filter[J]. Journal of Hydrology ， 2009， 376: 428－42.

[9] Mclaughlin D. An integrated approach to hydrologic data assimilation: interpolation， smoothing， and filtering[J]. Advances in Water Resourses， 2002， 25: 1275－86.

[10] Bosilovich M G， Radakovich J D， Silva A， et al. Skin Temperature Analysis and Bias Correction in a Coupled Land-Atmosphere Data Assimilation System[J]. Journal of the Meteorological Society of Japan. Ser.II， 2007， 85A: 205－28.

[11] Walker J P， Houser P R， Reichle R. New Technologies Require Advances in Hydrologic Data Assimilation[J]. EOS，2003， 84: 545－51.

[12] Sellers P J， Randall D A， Collatz G J， et al. A revised land surface parameterization （SiB2） for atmospheric GCMs. Part I: Model Formulation[J]. Journal of Climate， 1996， 9: 676－705.

[13] Genuchten M T. A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal， 1980， 44: 892－898.

[14] Yang Dawen， Tamagawa K. Hydrological improvement of the land surface process scheme using the CEOP observation data[C]//Proceeding of CEOP/IGWCO Joint Meeting. Tokyo，2005.

[15] Su Zhongbo. The surface energy balance system （SEBS） for estimation of turbulent heat fluxes[J]. Hydrological Earth System Science， 2002， 6: 85－99.

[16] Yang Dawen， Chen He， Lei Huimin. Estimation of evapotranspiration using a remote sensing model over agricultural land in the North China Plain[J]. International Journal of Remote Sensing， 2010， 31（14）: 3783－98.

Data assimilation technique based on ensemble Kalman filter for improving soil water content estimation

Chen He1, Yang Dawen2, Liu Yu1, Zhang Baozhong1
（1.Stɑte Key Lɑborɑtory of Simulɑtion ɑnd Regulɑtion of Wɑter Cycle in River Bɑsin， Chinɑ Institute of Wɑter Resources ɑnd Hydropower Reseɑrch， Beijing 100048， Chinɑ;2. Depɑrtment of Hydrɑulic Engineering, Stɑte Key Lɑborɑtory of Hydroscience ɑnd Engineering，Tsinghuɑ University， Beijing 100084， Chinɑ）

Abstract:Soil moisture in unsaturated zone connects the water and energy exchange process between surface water and groundwater， which has great influence to rainfall infiltration and surface evapotranspiration and hence has important meaning to agriculture， hydrology and meteorology. In agricultural study， accurate estimation of soil water content has significant importance to agricultural water management， irrigation regime determination and agricultural output increase. Soil water content can be quantified by surface observation， model estimation and remote sensing retrieval. Due to the large heterogeneity of soil property， surface observation at small scale can be hardly extended to large scale; while spatial distribution retrieved by remote sensing data can only obtain instantaneous value at satellite over-passing time. Land surface model is treated as a powerful tool in continuous estimation of soil water content， which is continuous in spatial and temporal dimension. However，the error tends to accumulate in the process of model simulation due to the inevitable uncertainty of forcing data and intrinsic error in model. Data assimilation technique can consider the uncertainty of the model and observation， update model states during the simulation period， and thus improve the accuracy of soil water content estimation， and exploit the advantages of both land surface model and remote sensing measurement. The concept and algorithm of data assimilation were first proposed by oceanologists and meteorologists， and have been gradually introduced to hydrology in recent years， such as soil water data assimilation and surface temperature data assimilation. As the development of remote sensing technique， more and more surface parameters can be obtained by remote sensing retrieval， which provides the available data sources for data assimilation system. The purpose of this study was to validate the data assimilation technique in improving soil water content estimation. To this end， an ensemble Kalman filter （EnKF） technique was coupled to a hydrologically-enhanced land process （HELP）model to update model states including soil water content and surface temperature. Random disturbance was added to the input data to generate ensemble model states and background error covariance matrix. The latent heat flux derived by MODIS data and surface energy balance system （SEBS） was used as the observation value of assimilation system to update the model states in HELP model. We chose a typical cropland in Weishan irrigation area （36°8′-37°1′N， 115°25′-116°31′E） as the study area，where located an eco-hydrological station （36°38′55.5″N， 116°3′15.3″E， average sea altitude of 30 m） with long series of flux data and meteorological measurements. The observation data used in this study were composed of flux observation data including soil heat flux， sensible and latent heat flux， meteorological observation data including rainfall， sunshine duration， air temperature and humidity， wind direction and speed， upward/downward longwave/shortwave radiation and infrared surface temperature， and vegetation observation data including canopy height and leaf area index. The model was firstly validated by the observation data in 2006， in which the open-loop estimation without state updating was treated as the benchmark run. The root mean square error （RMSE） of soil water content in surface， root and deep layer was 0.055， 0.053 and 0.053 m3/m3respectively. After data assimilation update， the surface temperature estimation of both wheat season and maize season was improved to a large extent， with an effectiveness coefficient of 94.8% and 73.0% respectively. Data assimilation also improved the estimation accuracy of soil water content， with a reduction of RMSE by 30%-50% compared to the benchmark run. In wheat season， the effectiveness coefficient of soil water content estimation of data assimilation in surface， root and deep layer ranged from 9.31% to 74.17%. Compared to wheat season， data assimilation showed better results in maize season， the relative error of soil water content in surface， root and deep layer was reduced to 3.70%， 3.62%， and -5.45%， respectively， and the effectiveness coefficient of all 3 layers was over 60%. These results demonstrate that the effect of data assimilation on improving soil water states is positive， which provides a new approach in continuous estimation of soil water content.

Keywords:soils; remote sensing; temperature; data assimilation; land surface model; soil water content

作者簡(jiǎn)介：陳鶴，女，遼寧大連人，中國(guó)水利水電科學(xué)研究院工程師，博士。北京 中國(guó)水利水電科學(xué)研究院，流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，100048。Email：chenhe@iwhr.com

基金項(xiàng)目：國(guó)家自然科學(xué)基金（51409277， 51379217）；“十二五”國(guó)家科技計(jì)劃（2012BAD08B01）

收稿日期：2015-08-03

修訂日期：2015-12-17

中圖分類號(hào)：S152.7

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1002-6819（2016）-02-0099-06

doi：10.11975/j.issn.1002-6819.2016.02.015