基于集合卡爾曼濾波的土壤水分同化變權(quán)實(shí)驗(yàn)研究

劉 茜,王明玉,趙英時(shí)

(中國(guó)科學(xué)院研究生院,北京 100049）

基于集合卡爾曼濾波的土壤水分同化變權(quán)實(shí)驗(yàn)研究

劉 茜,王明玉＊,趙英時(shí)

(中國(guó)科學(xué)院研究生院,北京 100049）

將集合預(yù)報(bào)成員不等權(quán)重思想與集合卡爾曼濾波(En KF）同化方法相結(jié)合,利用集合成員的離散度作為權(quán)重因子,對(duì)En KF算法優(yōu)化后的集合成員采用不等權(quán)重取平均值,作為同化后的預(yù)報(bào)值。首先檢驗(yàn)了集合離散度和預(yù)報(bào)誤差的相關(guān)性,證明將集合離散度作為權(quán)重因子的可靠性;利用一個(gè)水文過(guò)程模型(DHSVM）和實(shí)測(cè)數(shù)據(jù)進(jìn)行了土壤水分的同化變權(quán)實(shí)驗(yàn),對(duì)EnKF分析和更新后產(chǎn)生的土壤水分集合,分別采用算術(shù)平均和變權(quán)平均的方法,計(jì)算土壤水分預(yù)報(bào)結(jié)果并進(jìn)行比較。實(shí)驗(yàn)表明,集合變權(quán)平均法可以進(jìn)一步提高同化的預(yù)報(bào)效果。

數(shù)據(jù)同化;集合預(yù)報(bào);權(quán)重;En KF;土壤水分

0 引言

集合卡爾曼濾波(En KF）算法目前被廣泛應(yīng)用于大氣、海洋、陸地等數(shù)據(jù)同化研究和應(yīng)用中[1-3],許多學(xué)者將 En KF與其他同化方法進(jìn)行了比較研究[4,5],進(jìn)一步證明了 En KF方法的有效性和優(yōu)越性。En KF基于集合預(yù)報(bào)的思想發(fā)展起來(lái),它通過(guò)一組集合預(yù)報(bào)得到預(yù)報(bào)誤差協(xié)方差矩陣,即從集合預(yù)報(bào)估計(jì)狀態(tài)變量與觀測(cè)變量之間的協(xié)方差,利用觀測(cè)資料和協(xié)方差,通過(guò)卡爾曼濾波方程更新集合預(yù)報(bào)結(jié)果。為了提高計(jì)算精度,En KF需要很大的集合數(shù)目,但是在利用En KF同化之后,只是將所有集合成員的算術(shù)平均值作為預(yù)報(bào)結(jié)果,集合中每個(gè)預(yù)報(bào)結(jié)果具有相同的權(quán)重,這并不能夠完全解釋集合分布中的信息。如果集合中存在與大多數(shù)成員差別非常大的成員,其將對(duì)集合平均結(jié)果有較大的影響。集合預(yù)報(bào)的結(jié)果更傾向于多數(shù)集合預(yù)報(bào)成員的值,高預(yù)報(bào)成功率往往對(duì)應(yīng)于較多成員比較一致的結(jié)果[6],基于變權(quán)集合平均法的集合預(yù)報(bào)可以獲得更好的預(yù)報(bào)結(jié)果[7,8]。變權(quán)集合平均法的出發(fā)點(diǎn)就是通過(guò)對(duì)集合成員權(quán)重的調(diào)整,加大預(yù)報(bào)結(jié)果較為一致的集合成員的權(quán)重,強(qiáng)調(diào)這種一致性對(duì)預(yù)報(bào)結(jié)果的影響。在En KF的集合數(shù)據(jù)同化研究中尚未運(yùn)用該方法。

本文將集合預(yù)報(bào)成員不等權(quán)重的思想與 En KF同化方法相結(jié)合,在不影響 En KF算法對(duì)誤差協(xié)方差矩陣計(jì)算,即不影響整個(gè)同化方法優(yōu)化分析和更新的情況下,以 En KF算法更新后集合成員的離散度作為權(quán)重因子,采用不等權(quán)重的集合平均值替代簡(jiǎn)單算術(shù)平均值作為最后的預(yù)報(bào)值,以期獲得更好的預(yù)報(bào)效果。

1 算法介紹

1.1 集合卡爾曼濾波(En KF）

Evensen根據(jù) Ep stein的隨機(jī)動(dòng)態(tài)預(yù)報(bào)理論提出集合卡爾曼濾波(En KF）算法[9],Burgers對(duì)該算法進(jìn)行了改進(jìn),提出擾動(dòng)觀測(cè)的 En KF算法[10],其思想和計(jì)算如下:

(1）初始化背景場(chǎng)。給定 N個(gè)符合高斯分布的狀態(tài)變量Xi(i=1,…,N）,計(jì)算每個(gè)狀態(tài)變量在 k+ 1時(shí)刻的預(yù)報(bào)值 Xfi,k+i:

式中:Xif,k+i、Xai,k分別為第i個(gè)狀態(tài)變量在k+1時(shí)刻的預(yù)報(bào)值和在 k時(shí)刻的分析值;Mk(■）表示非線性模型算子;Wk為模型誤差方差矩陣;wi,k表示期望為0、方差為Wk的高斯白噪聲。

(2）計(jì)算 k+1時(shí)刻的卡爾曼增益矩陣 Kk+1:

(4）進(jìn)入下一時(shí)刻,返回步驟(2）。

1.2 集合變權(quán)平均法

在En KF算法優(yōu)化更新后,會(huì)產(chǎn)生一組優(yōu)化后的狀態(tài)變量集合。為了提高集合預(yù)報(bào)的效果,考慮到集合成員之間的差異,選用集合成員之間離散度w作為權(quán)重因子,對(duì)集合平均值進(jìn)行調(diào)整。

集合成員的離散度w可以定義為:

集合中每個(gè)元素的權(quán)重可以認(rèn)為與其離散度成反比,離散度越大,賦權(quán)重越小。集合變權(quán)平均預(yù)報(bào)結(jié)果′可表示為:

1.3 變權(quán)集合同化效果檢驗(yàn)方法

利用En KF算法對(duì)狀態(tài)變量進(jìn)行更新,產(chǎn)生優(yōu)化后的狀態(tài)變量集合,對(duì)集合分別采用等權(quán)平均法和變權(quán)平均法計(jì)算預(yù)報(bào)結(jié)果;然后利用實(shí)測(cè)的土壤水分?jǐn)?shù)據(jù),采用平均誤差(M B E）和均方根誤差(RM SE）對(duì)預(yù)報(bào)效果進(jìn)行檢驗(yàn)。公式為:

式中:P為預(yù)報(bào)值;O為實(shí)地觀測(cè)值;n為觀測(cè)的數(shù)量。

2 同化變權(quán)實(shí)驗(yàn)

為檢驗(yàn)此方法的有效性,基于En KF算法,將實(shí)測(cè)的土壤水分?jǐn)?shù)據(jù)與分布式水文土壤植被模型(DHSVM）模擬的土壤水分進(jìn)行同化,對(duì)優(yōu)化后的土壤水分集合分別采用等權(quán)平均和變權(quán)平均的方法,計(jì)算預(yù)報(bào)的土壤水分,在黑河流域進(jìn)行了土壤水分同化變權(quán)實(shí)驗(yàn)。

2.1 研究區(qū)與數(shù)據(jù)

黑河流域是我國(guó)第二大內(nèi)河流域,流域面積11.6萬(wàn)km2,年均降水量108 mm。研究區(qū)位于黑河流域中游(38°45′～39°15′N,100°00′～100°45′E）,地面高程在1 000～2 000 m,年均氣溫6℃。研究區(qū)內(nèi)的自動(dòng)氣象觀測(cè)站記錄了水文模型所需的6個(gè)氣象驅(qū)動(dòng)變量:氣溫、風(fēng)速、相對(duì)濕度、短波輻射、長(zhǎng)波輻射、降水;在每個(gè)氣象站放置了 TDR儀器,用于測(cè)量10 cm、20 cm、40 cm、80 cm、120 cm、160 cm深度的土壤水分。實(shí)驗(yàn)在 HUAZHA IZI站點(diǎn)進(jìn)行,此處主要為低矮稀疏灌草覆蓋的荒漠,土壤以砂礫土為主。

2.2 分布式水文土壤植被模型(DHSVM）

DHSVM模型詳細(xì)描述了地形和植被對(duì)水流動(dòng)的影響,模型主要包括植被截留、蒸散發(fā)、積雪、融雪和徑流等部分。DHSVM在流域數(shù)字高程模型(DEM）的網(wǎng)格尺度上對(duì)水文過(guò)程進(jìn)行動(dòng)態(tài)描述,以DEM的節(jié)點(diǎn)為中心,流域被分成若干計(jì)算網(wǎng)格單元,每個(gè)網(wǎng)格被賦予不同的土壤和植被特性。每一計(jì)算時(shí)段內(nèi),模型對(duì)流域內(nèi)各網(wǎng)格的能量平衡方程和質(zhì)量平衡方程提供聯(lián)立解,各網(wǎng)格之間通過(guò)坡面匯流和壤中流的匯流演算發(fā)生水文聯(lián)系。

非飽和土壤水運(yùn)動(dòng)采用兩層模式進(jìn)行模擬,透過(guò)植被下滲的降水以及雪蓋融水都進(jìn)入土壤。第一層的深度是下冠層的根系深度,上層土壤的厚度 d1等于下冠層植被根系平均深度;第二層的深度是上冠層的根系深度,下層土壤厚度 d2從 d1處伸展到上冠層植被根系的平均深度。下冠層一般只從第一層中提取水分,但上冠層能從兩層中提取水分;如果土壤達(dá)到飽和狀態(tài),多余的水分將形成地面徑流。土壤蒸散發(fā)的計(jì)算一般只限制在第一層。該模型的詳細(xì)描述可參考文獻(xiàn)[11,12]。

2.3 水文模型輸入和同化實(shí)驗(yàn)

DHSVM模型的輸入數(shù)據(jù)包括DEM、土壤類型、土層厚度、植被類型、氣象數(shù)據(jù)、河網(wǎng)/路網(wǎng)分布圖和道路/河流分級(jí)文件,其中植被和土壤類型數(shù)據(jù)由中科院資源環(huán)境科學(xué)數(shù)據(jù)中心提供。以30 m× 30 m的DEM網(wǎng)格作為所有輸入數(shù)據(jù)的基本網(wǎng)格單元,并利用反距離加權(quán)插值法對(duì)輸入數(shù)據(jù)進(jìn)行重采樣。將DHSVM模型的輸入變量作為柵格圖層加入A rcGIS,并將這些柵格圖層Grid格式轉(zhuǎn)換成Binary格式,便于DHSVM模型讀取和處理。在每個(gè)時(shí)間步長(zhǎng)和每個(gè)單元網(wǎng)格上模型會(huì)自動(dòng)運(yùn)行。DHSVM模型運(yùn)用許多刻畫(huà)植被和土壤特征的參數(shù)描述地表過(guò)程,這些參數(shù)主要通過(guò)地面實(shí)測(cè)和文獻(xiàn)[13]獲得。在本研究中,土壤水分的模擬時(shí)間步長(zhǎng)為3 h。

基于En KF算法,每隔3 d加入一次實(shí)測(cè)表層土壤水分?jǐn)?shù)據(jù)(0～10 cm）,與DHSVM模型模擬進(jìn)行數(shù)據(jù)同化。對(duì) En KF而言,模型誤差和觀測(cè)誤差的準(zhǔn)確估計(jì)十分重要,它直接影響著同化效果。然而,觀測(cè)誤差和模型誤差隨著時(shí)間不斷變化,而且其影響因素非常多,所以對(duì)觀測(cè)誤差和模型誤差的量化非常困難。假設(shè)本研究中實(shí)測(cè)土壤水分是精確的,先將DHSVM模型在同化實(shí)驗(yàn)之前運(yùn)行一個(gè)月,在初始條件并不是很穩(wěn)定的條件下,模型自動(dòng)進(jìn)行調(diào)整而達(dá)到平衡,同時(shí)獲得模型的初始場(chǎng)和初始場(chǎng)誤差方差。根據(jù)初始場(chǎng)誤差方差對(duì)土壤水分的初始場(chǎng)進(jìn)行擾動(dòng),產(chǎn)生水文模型的初始狀態(tài)集合。此外,對(duì)站點(diǎn)觀測(cè)的表層土壤水分進(jìn)行擾動(dòng),增加均值為0、方差為2E-5的隨機(jī)噪聲,生成觀測(cè)數(shù)據(jù)的集合,集合數(shù)目設(shè)為100。通常在En KF算法中,在每個(gè)時(shí)間步長(zhǎng),最后的狀態(tài)變量取值由同化預(yù)測(cè)的集合算術(shù)平均值確定。在 En KF同化變權(quán)實(shí)驗(yàn)中,模型運(yùn)行的每個(gè)時(shí)間步長(zhǎng)最后的狀態(tài)變量取值由同化預(yù)測(cè)的集合變權(quán)平均值確定。在進(jìn)行同化實(shí)驗(yàn)的同時(shí),單獨(dú)運(yùn)行DHSVM模型,并將其運(yùn)行結(jié)果與等權(quán)平均預(yù)測(cè)、變權(quán)平均預(yù)測(cè)的同化結(jié)果相比較。

3 實(shí)驗(yàn)結(jié)果與分析

首先基于等權(quán)平均方法,利用 En KF算法進(jìn)行同化實(shí)驗(yàn),在不同的時(shí)間步長(zhǎng)進(jìn)行集合預(yù)報(bào)。為了進(jìn)一步驗(yàn)證將集合離散度作為權(quán)重因子的可靠性,對(duì)集合離散度和預(yù)報(bào)誤差進(jìn)行相關(guān)分析(圖1）。由圖1可知,同化后的預(yù)報(bào)誤差隨著集合離散度的增大而增大,預(yù)報(bào)誤差與集合離散度之間具有較好的線性相關(guān)性(R2=0.7843）。因此,為了提高預(yù)報(bào)精度,獲得更多集合預(yù)報(bào)的有效信息,將集合離散度作為En KF集合預(yù)報(bào)的權(quán)重因子取平均值是合理的。

圖1 預(yù)報(bào)誤差與集合離散度的相關(guān)性Fig.1 The relationship between forecast errorsand the ensemble spread

實(shí)驗(yàn)采用均方根誤差RM SE和平均誤差MB E,對(duì)2008年5月27日(148 d）—7月10日(190 d）HUAZHA IZI站點(diǎn)模型單獨(dú)模擬、等權(quán)平均預(yù)測(cè)、變權(quán)平均預(yù)測(cè)結(jié)果進(jìn)行分析,結(jié)果如圖2和表1。由表1可知,利用等權(quán)平均預(yù)報(bào)和變權(quán)平均預(yù)報(bào),同化后表層和根區(qū)土壤水分的模擬值均方根誤差和平均誤差均明顯降低。此外,對(duì)于表層土壤水分而言,一般等權(quán)平均和變權(quán)平均預(yù)報(bào)的模擬值的RM SE分別為0.0047和0.0045,MB E分別為-0.0018和-0.0012;對(duì)于根區(qū)土壤水分而言,等權(quán)平均預(yù)報(bào)和變權(quán)平均預(yù)報(bào)的模擬值的RM SE分別為0.0029和0.0014,M B E分別為-0.0019和-0.0014,變權(quán)平均預(yù)報(bào)的優(yōu)化效果明顯優(yōu)于等權(quán)平均預(yù)報(bào)。研究表明,在En KF同化更新后,取集合的變權(quán)平均值,將實(shí)地觀測(cè)與水文模型同化可以明顯改善土壤水分的估測(cè)結(jié)果,同化使得DHSVM模型輸出更接近于實(shí)測(cè)值。

圖2 實(shí)地觀測(cè)、DHSVM模型模擬、等權(quán)平均預(yù)報(bào)和變權(quán)平均預(yù)報(bào)同化結(jié)果比較Fig.2 Comparisons of field measurement,DHSVM simulation and assim ilation results

表1 土壤水分模擬和同化結(jié)果的誤差統(tǒng)計(jì)Table 1 The soilmoisture errors statistic of simulation and assim ilation results

4 結(jié)論

本文提出將集合預(yù)報(bào)成員不等權(quán)重思想與 En-KF相結(jié)合的方法,以En KF算法更新后集合成員的離散度作為權(quán)重因子,對(duì)集合成員權(quán)重進(jìn)行調(diào)整,以集合變權(quán)平均值作為最后的預(yù)報(bào)值。研究利用水文過(guò)程模型(DHSVM）和實(shí)測(cè)數(shù)據(jù)進(jìn)行土壤水分同化實(shí)驗(yàn),將等權(quán)平均和變權(quán)平均的預(yù)報(bào)結(jié)果進(jìn)行比較,結(jié)論如下:1）基于En KF算法,將實(shí)測(cè)數(shù)據(jù)與水文過(guò)程模型(DHSVM）進(jìn)行同化,利用等權(quán)平均法和變權(quán)平均法,均可以明顯改進(jìn)單獨(dú)水文模型模擬的土壤水分結(jié)果;2）對(duì)于En KF優(yōu)化后的集合預(yù)報(bào)結(jié)果,集合的離散度和 En KF預(yù)報(bào)誤差具有較好的相關(guān)性,隨著集合離散度的增加,預(yù)報(bào)誤差也增加,集合離散度可以在一定程度上反映預(yù)報(bào)的效果;3）利用集合的離散度作為權(quán)重因子,將En KF濾波與變權(quán)重方法相結(jié)合,將不同權(quán)重的集合平均值替代原算術(shù)平均,可進(jìn)一步提高集合的平均效果和 En KF的同化效果。本方法對(duì)集合數(shù)據(jù)同化方法具有一定的借鑒作用。

[1] 張學(xué)峰,黃大吉,章本照,等.集合數(shù)據(jù)同化方法的發(fā)展與應(yīng)用概述[J].海洋學(xué)研究,2007,25(3）:89-96.

[2] ALLEN J I,EKNES M,EVENSEN G.An ensemble Kalman filter with a complex marine ecosystem model:Hindcasting phytoplankton in the Cretan Sea[J].Annals Geophysicae,2002,20: 1-13.

[3] DONG Y N,GU Y Q,OL IVER D S.Sequential assimilation of 4D seismic data for reservoir description using the ensemble Kalman filter[J].Journal of Petroleum Science and Engineering,2006,53(1）:83-99.

[4] WALKER J P,W ILLGOOSE G R,KALMA J D.One-dimensional soilmoisture p rofile retrieval by assimilation of near-surface observations:A comparison of retrieval algorithms[J].Advances in Water Resources,2001,24(6）:631-650.

[5] HAN X J,L I X.An evaluation of the nonlinear/non-Gaussian filters for the sequential data assimilation[J].Remote Sensing of Environment,2008,112(4）:1434-1449.

[6] TOTH Z,ZHU Y J,MARCHOK T.The use of ensemble to identify forecasts with small and large uncertainty[J].Weather Forecasting,2001,16(8）:463-477.

[7] HOU TEKAM ER P L.Methods fo r ensemble p rediction[J]. Monthly Weather Review,1995,123:2181-2196.

[8] 段明鏗,王盤(pán)興.一種新的集合預(yù)報(bào)權(quán)重平均方法[J].應(yīng)用氣象學(xué)報(bào),2006,17(4）:488-493.

[9] EVERNSEN G.Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte-Carlo methods to fo recast erro r statistics[J].Journal of Geophysical Research,1994,99 (C5）:10143-10162.

[10] BERGERSG,LEEUWEN PJ,EVENSEN G.Analysis scheme in the ensemble Kalman filter[J].Monthly Weather Review, 1998,129:2884-2903.

[11] W IGMOSTA M S,VA IL L,LETTENMA IER D P.A distributed hydrology-vegetation model for complex terrain[J].Water Resource Research,1994,30(6）:1665-1679.

[12] W IGMOSTA M S,N IJSSEN B,STORCK P,et al.The Distributed Hydrology Soil Vegetation Model,in Mathematical Models of Small Watershed Hydrology and Applications[M]. Water Resource Publications Littleton CO.,2002.7-42.

[13] 康麗莉,王守榮,顧駿強(qiáng).分布式水文模型DHSVM對(duì)蘭江流域徑流變化的模擬試驗(yàn)[J].熱帶氣象學(xué)報(bào),2008,24(2）:176 -182.

A Weighted Average Soil M oisture Assim ilation Experiment Based on Ensemble Kalman Filter

L IU Qian,WANGM ing-yu,ZHAO Ying-shi
(Graduate University of Chinese Academ y of Sciences,Beijing 100049,China）

The arithmetic average of the ensemblemembers is usually taken as the final assimilation result of Ensemble Kalman Filter(EnKF）.In this study,theweighted averageof the forecasting ensembleswas used as the assimilated resultof EnKF.The dispersion of the ensemblememberswas considered as the weighting facto r.After the EnKF has analyzed and updated each ensemblemember,the weighted average of the ensemblemembers was adop ted as the final fo recasting results.Since the relationship between the dispersion of ensemblemembers and the forecast errors was significant,taking the dispersion of the ensemble members as the weighting factor was reasonable.Thisapp roach wasexamined based on a soilmoisture assimilation experiment, w hich assimilated the in situ measurements of the surface soil moisture w ith a distributed hydrologic soil vegetation model (DHSVM）using EnKF.The resultsof the weighted averagewere compared w ith the resultsof arithmetic average based on En-KF.This study indicates that the assimilation can be imp roved significantly using the weighted average of ensemble members.

data assimilation;ensemble fo recasting;weight;EnKF;soilmoisture

P208

A

1672-0504(2010）01-0094-04

2009-10-09;

2009-12-10

國(guó)家973項(xiàng)目“陸表生態(tài)環(huán)境要素主被動(dòng)遙感協(xié)同反演理論與方法”(2007CB714407）;中國(guó)科學(xué)院百人計(jì)劃擇優(yōu)支持項(xiàng)目“流域尺度復(fù)雜地質(zhì)介質(zhì)地下水污染與水資源安全”

劉茜(1983-）,女,博士研究生,主要從事遙感與水文模型的同化研究。＊通訊作者E-mail:mwang@gucas.ac.cn