劉 茜,王明玉,趙英時

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

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

劉 茜,王明玉＊,趙英時

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

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

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

0 引言

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

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

1 算法介紹

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

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

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

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

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

(4）進入下一時刻,返回步驟(2）。

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

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

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

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

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

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

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

2 同化變權(quán)實驗

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

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

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

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

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

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

2.3 水文模型輸入和同化實驗

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

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

3 實驗結(jié)果與分析

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

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

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

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

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

4 結(jié)論

本文提出將集合預(yù)報成員不等權(quán)重思想與 En-KF相結(jié)合的方法,以En KF算法更新后集合成員的離散度作為權(quán)重因子,對集合成員權(quán)重進行調(diào)整,以集合變權(quán)平均值作為最后的預(yù)報值。研究利用水文過程模型(DHSVM）和實測數(shù)據(jù)進行土壤水分同化實驗,將等權(quán)平均和變權(quán)平均的預(yù)報結(jié)果進行比較,結(jié)論如下:1）基于En KF算法,將實測數(shù)據(jù)與水文過程模型(DHSVM）進行同化,利用等權(quán)平均法和變權(quán)平均法,均可以明顯改進單獨水文模型模擬的土壤水分結(jié)果;2）對于En KF優(yōu)化后的集合預(yù)報結(jié)果,集合的離散度和 En KF預(yù)報誤差具有較好的相關(guān)性,隨著集合離散度的增加,預(yù)報誤差也增加,集合離散度可以在一定程度上反映預(yù)報的效果;3）利用集合的離散度作為權(quán)重因子,將En KF濾波與變權(quán)重方法相結(jié)合,將不同權(quán)重的集合平均值替代原算術(shù)平均,可進一步提高集合的平均效果和 En KF的同化效果。本方法對集合數(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] 段明鏗,王盤興.一種新的集合預(yù)報權(quán)重平均方法[J].應(yīng)用氣象學(xué)報,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] 康麗莉,王守榮,顧駿強.分布式水文模型DHSVM對蘭江流域徑流變化的模擬試驗[J].熱帶氣象學(xué)報,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

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

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