梁 曉鄭小谷戴永久師春香

1）（國(guó)家氣象信息中心，北京100081）2）（北京師范大學(xué)，北京100875）

En KF中誤差協(xié)方差優(yōu)化方法及在資料同化中應(yīng)用

梁 曉1）*鄭小谷2）戴永久2）師春香1）

1）（國(guó)家氣象信息中心，北京100081）2）（北京師范大學(xué)，北京100875）

集合卡爾曼濾波（the Ensemble Kalman Filter，簡(jiǎn)稱En KF）中將預(yù)報(bào)集合的統(tǒng)計(jì)協(xié)方差作為預(yù)報(bào)誤差協(xié)方差，但該估計(jì)可能嚴(yán)重偏離真實(shí)的預(yù)報(bào)誤差協(xié)方差，影響同化精度?；跇O大似然估計(jì)理論，發(fā)展了一種優(yōu)化預(yù)報(bào)誤差協(xié)方差矩陣的實(shí)時(shí)膨脹方法，即MLE（the Maximum Likelihood Estimation）方法。利用蒙古國(guó)基準(zhǔn)站Delgertsgot（簡(jiǎn)稱DGS站）觀測(cè)資料，基于En KF方法和 MLE方法，在通用陸面模式（the Common Land Model，簡(jiǎn)稱Co LM）中同化了地表溫度和10 cm土壤溫度觀測(cè)資料，建立了土壤溫度同化系統(tǒng)。結(jié)果表明：MLE方法對(duì)地表溫度和各層土壤溫度（尤其深層土壤溫度）的估計(jì)比En KF方法準(zhǔn)確。考慮到淺層和深層土壤溫度的差別，在實(shí)施MLE方法時(shí)對(duì)淺層和深層土壤溫度采用了不同的膨脹因子。對(duì)比膨脹因子為單一標(biāo)量時(shí)的結(jié)果，多因子膨脹能緩解深層土壤溫度的不合理膨脹，改善同化效果。

數(shù)據(jù)同化；集合卡爾曼濾波；誤差協(xié)方差膨脹

引 言

集合卡爾曼濾波（the Ensemble Kalman Filter，簡(jiǎn)稱En KF）是目前常用的一種數(shù)據(jù)同化方法，自Evensen［1］提出和Burgers等［2］發(fā)展后，已應(yīng)用于很多同化研究。理論上，通過將觀測(cè)資料融入預(yù)報(bào)模型，數(shù)據(jù)同化能夠?qū)崿F(xiàn)模型輸出和觀測(cè)資料的最佳融合。但這種“最佳”依賴于預(yù)報(bào)（背景）誤差協(xié)方差和觀測(cè)誤差協(xié)方差的精度。預(yù)報(bào)誤差協(xié)方差和觀測(cè)誤差協(xié)方差的估計(jì)是否準(zhǔn)確，是決定En KF方法同化精度高低的關(guān)鍵［3］。然而，實(shí)際應(yīng)用中很難精確得到這些統(tǒng)計(jì)量信息［4］。因此，優(yōu)化誤差協(xié)方差的估計(jì)，對(duì)En KF方法發(fā)展十分必要［5-6］。

En KF方法中，通常將預(yù)報(bào)集合的統(tǒng)計(jì)協(xié)方差作為真實(shí)的預(yù)報(bào)誤差協(xié)方差。但研究表明，由于樣本量有限，該估計(jì)中存在的采樣誤差會(huì)導(dǎo)致預(yù)報(bào)集合的統(tǒng)計(jì)協(xié)方差低估真正的預(yù)報(bào)誤差協(xié)方差，集合擴(kuò)散度迅速減小，使得觀測(cè)對(duì)同化的影響越來越小，甚至完全失效，發(fā)生濾波發(fā)散［7-8］。為此，一種旨在增大集合方差的協(xié)方差膨脹技術(shù)逐漸發(fā)展起來。目前常用的一種方法是給預(yù)報(bào)集合成員偏離平均值的偏差乘一個(gè)略大于1的常數(shù)［7］，該方法中的膨脹因子通過反復(fù)試驗(yàn)得到，需多次運(yùn)行同化系統(tǒng)。對(duì)于大型復(fù)雜的實(shí)際模型，這種試驗(yàn)方法的計(jì)算代價(jià)十分巨大，且膨脹因子在同化時(shí)段內(nèi)不變也不合理。因此，有必要發(fā)展無(wú)需反復(fù)試驗(yàn)并能在同化循環(huán)中實(shí)時(shí)更新膨脹因子的協(xié)方差膨脹技術(shù)。

Dee［9］及其相關(guān)工作［10-11］發(fā)展了誤差協(xié)方差的極大似然估計(jì)方法。首先對(duì)誤差協(xié)方差矩陣建立參數(shù)化模型，再通過觀測(cè)減預(yù)報(bào)得到的殘差的極大似然估計(jì)來估計(jì)模型中的未知參數(shù)。但其工作沒有提出可靠的誤差協(xié)方差矩陣參數(shù)化模型，因此一直未得到普及。近年來，根據(jù)觀測(cè)減預(yù)報(bào)得到的殘差的統(tǒng)計(jì)特征，逐漸提出了一些實(shí)時(shí)估計(jì)膨脹因子的方法［12-14］。然而，這些方法均基于矩估計(jì)，而非極大似然估計(jì)。

作為極大似然估計(jì)方法［9-11］的擴(kuò)展，Zheng［15］提出了預(yù)報(bào)誤差協(xié)方差矩陣的多變量實(shí)時(shí)膨脹方法，將膨脹因子從常數(shù)標(biāo)量擴(kuò)展到實(shí)時(shí)變化的對(duì)角矩陣。Liang等［16］基于Zheng［15］的工作，提出了同時(shí)膨脹預(yù)報(bào)和觀測(cè)誤差協(xié)方差矩陣的方法，但上述工作只在簡(jiǎn)單的小模型中進(jìn)行試驗(yàn)。作為文獻(xiàn)［16］工作的延續(xù)，本文使用陸面過程模型，將Zheng［15］和Liang等［16］的膨脹方法用于陸面同化。本文的膨脹因子已擴(kuò)展為對(duì)角矩陣，即非標(biāo)量。

土壤溫度是陸面過程中的一個(gè)重要變量，可以通過模型模擬或儀器觀測(cè)得到，但兩種手段各有利弊。站點(diǎn)觀測(cè)反映了觀測(cè)時(shí)間和地點(diǎn)上的真實(shí)狀態(tài)，但由于土壤的空間異質(zhì)性，它所能代表的空間范圍有限，無(wú)法滿足區(qū)域性研究的要求；遙感觀測(cè)可以得到衛(wèi)星過境時(shí)刻大范圍的地表溫度信息，但遙感反演過程引入了很大的不確定性，而且只能獲得瞬時(shí)信息。陸面過程模型可以給出土壤溫度在時(shí)間和空間上的連續(xù)演進(jìn)，但模型參數(shù)、驅(qū)動(dòng)場(chǎng)、初始條件等的不完備都會(huì)影響模擬精度。如何有效結(jié)合觀測(cè)和陸面過程模型模擬得到更真實(shí)且時(shí)空連續(xù)的信息，尤為重要，因此陸面數(shù)據(jù)同化研究越來越多［17-20］。本文將利用通用陸面模式（the Common Land Model，簡(jiǎn)稱Co LM）和兩種同化方法發(fā)展土壤溫度同化方案，改進(jìn)對(duì)土壤溫度廓線的估計(jì)。為便于表達(dá)，未使用和使用本文膨脹方法的En KF方法分別簡(jiǎn)稱為En KF方法和MLE方法。

1 MLE方法

一個(gè)非線性時(shí)間離散的預(yù)報(bào)系統(tǒng)和線性觀測(cè)系統(tǒng)可表示為

其中，i是積分步數(shù)，是真值向量，Mi－1是非線性模型算子，yoi是觀測(cè)向量，Hi是將狀態(tài)變量映射到觀測(cè)空間的觀測(cè)算子，ηi和εi是背景誤差和觀測(cè)誤差向量，假定二者相互獨(dú)立且時(shí)間上不相關(guān)，各自服從平均值為0、協(xié)方差矩陣分別為背景誤差協(xié)方差矩陣Pi和觀測(cè)誤差協(xié)方差矩陣Ri的高斯分布。同化目標(biāo)是找到盡可能接近真值的分析狀態(tài)

En KF方法在很多文獻(xiàn)中均有介紹［2，21-22］，不再贅述。En KF方法中通常采用預(yù)報(bào)集合的采樣協(xié)方差矩陣Pi（式（3））作為真正的預(yù)報(bào)誤差協(xié)方差矩陣

有限的樣本量導(dǎo)致該估計(jì)中存在較大的采樣誤差，集合協(xié)方差會(huì)低估真正的誤差協(xié)方差，導(dǎo)致濾波發(fā)散［23］。最初，抑制濾波發(fā)散的方法是給預(yù)報(bào)集合中每個(gè)成員偏離平均值的離差乘略大于1的因子相當(dāng)于將預(yù)報(bào)誤差協(xié)方差矩陣Pi膨脹為原來的λ倍。但該膨脹因子λ隨時(shí)間不變，且需要多次反復(fù)試驗(yàn)來確定［7］。

基于觀測(cè)減預(yù)報(bào)得到的殘差di的極大似然估計(jì)，Zheng［15］和 Liang等［16］提出一種隨時(shí)間變化的膨脹方法，通過極小化di的似然函數(shù)－2Li（λi）（即對(duì)di的概率密度函數(shù)?。?倍的自然對(duì)數(shù)）實(shí)時(shí)估計(jì)膨脹因子［λi］：

其中，det表示矩陣的行列式。該方法根據(jù)文獻(xiàn)［9-11］提出的估計(jì)誤差協(xié)方差方法發(fā)展而來，膨脹因子［λi］為對(duì)角矩陣，不是標(biāo)量。

估計(jì)得到［λi］后，預(yù)報(bào)誤差協(xié)方差矩陣假定為［λi］Pi［λi］，相當(dāng)于用膨脹后的預(yù)報(bào)集合重新計(jì)算式（3）中定義的預(yù)報(bào)誤差協(xié)方差矩陣：

結(jié)合預(yù)報(bào)和觀測(cè)信息，可計(jì)算出i時(shí)刻同化更新后的分析集合

式（7）中，ε′i，j為服從平均值為0、協(xié)方差矩陣為Ri的正態(tài)分布隨機(jī)擾動(dòng)。

本文將使用同化結(jié)果的均方根誤差來評(píng)價(jià)同化效果。通常來說，均方根誤差越小，同化結(jié)果越好。

2 試驗(yàn)設(shè)計(jì)

本文選取國(guó)際協(xié)同強(qiáng)化觀測(cè)期計(jì)劃（Coordinated Energy and Water Cycle Observations Project，簡(jiǎn)稱CEOP）蒙古國(guó)基準(zhǔn)站Delgertsgot（簡(jiǎn)稱DGS站），使用 En KF和 MLE兩種同化方法，在Co LM［24-25］中以12 h 1次的頻率同化DGS站地表溫度和10 cm土壤溫度觀測(cè)資料，比較兩種方法對(duì)土壤溫度廓線的估計(jì)效果。同化系統(tǒng)狀態(tài)向量包括地表溫度、10層土壤溫度和10層土壤濕度。本文還分別使用經(jīng)驗(yàn)算法得到的葉面積指數(shù)和MODIS LAI產(chǎn)品，用來分析葉面積指數(shù)精度對(duì)土壤溫度的影響。

2.1 研究區(qū)概況及觀測(cè)數(shù)據(jù)

蒙古國(guó)基準(zhǔn)站位于蒙古高原，地表覆蓋為均質(zhì)矮草，觀測(cè)數(shù)據(jù)較好［18，26］。其中 DGS站（46.1°N，106.4°E）于2002年10月1日—2003年9月30日進(jìn)行了地表溫度和土壤溫度及氣象要素的連續(xù)觀測(cè)，本文選取該站開展同化試驗(yàn)。土壤溫度的觀測(cè)深度為3，10，40 cm和100 cm，相應(yīng)還有土壤濕度觀測(cè)，30 min 1次。在該站點(diǎn)已開展了一些陸面同化研究［17-18］，本文在選取站點(diǎn)觀測(cè)資料和研究時(shí)段等方面均予以考慮。本研究中，在Co LM中同化了DGS站每日08：00（世界時(shí)，下同）和20：00的地表溫度和10 cm土壤溫度的觀測(cè)數(shù)據(jù)，并使用逐時(shí)連續(xù)的各層土壤溫度和土壤濕度觀測(cè)數(shù)據(jù)來驗(yàn)證同化結(jié)果。

Co LM運(yùn)行所需的大氣驅(qū)動(dòng)數(shù)據(jù)中，除入射太陽(yáng)短波輻射和大氣長(zhǎng)波輻射沒有觀測(cè)資料外，其他變量均來自自動(dòng)氣象站資料。日本氣象廳利用其全球數(shù)據(jù)同化系統(tǒng)為CEOP的基準(zhǔn)站提供了始于2002年10月1日的時(shí)間分辨率為1 h的大氣數(shù)據(jù)集，Huang等［17］從該數(shù)據(jù)集獲取了蒙古國(guó)基準(zhǔn)站的太陽(yáng)短波輻射和大氣長(zhǎng)波輻射，本研究所需的輻射信息來源于此。

為分析葉面積指數(shù)（leaf area index，簡(jiǎn)稱LAI）對(duì)土壤溫度的影響，本研究還使用MODIS LAI產(chǎn)品代替原模式經(jīng)驗(yàn)算法計(jì)算的葉面積指數(shù)，比較使用經(jīng)驗(yàn)葉面積指數(shù)和MODIS LAI產(chǎn)品得到的結(jié)果。MODIS LAI產(chǎn)品提供了全球的1 km分辨率的葉面積指數(shù)（每8 d 1次），本文從中提取DGS站的2002年10月1日—2003年9月30日的葉面積指數(shù)和質(zhì)量控制說明數(shù)據(jù)。

2.2 MLE方法的兩種情況

本文使用的膨脹因子［λi］是對(duì)角矩陣，而不再局限為標(biāo)量。當(dāng)［λi］的主對(duì)角線上含有不同數(shù)值時(shí)，即可對(duì)不同的狀態(tài)變量使用不同系數(shù)進(jìn)行不同程度的膨脹；而當(dāng)［λi］的主對(duì)角線上均為相同數(shù)值時(shí)，膨脹因子［λi］便相當(dāng)于標(biāo)量，即對(duì)所有的狀態(tài)變量都使用同一系數(shù)進(jìn)行相同程度的膨脹。

Co LM中地表溫度和深層土壤溫度的性質(zhì)差異大，如地表溫度隨時(shí)間變化快且變化幅度大，而深層土壤溫度隨時(shí)間變化平緩。如果對(duì)整個(gè)土壤溫度廓線均采用同一膨脹因子，可能無(wú)法兼顧淺層和深層土壤溫度膨脹的合理性，最終影響土壤溫度廓線的同化效果。因此，本文對(duì)淺層和深層土壤溫度采用不同的膨脹因子，即膨脹因子矩陣［λi］的主對(duì)角線上含有兩個(gè)待估計(jì)的標(biāo)量膨脹因子λi，1和λi，2。其中［λi］主對(duì)角線上的前3個(gè)元素值均為λi，1，用于膨脹地表溫度和上兩層土壤溫度（2.8 cm以上）；后8個(gè)元素值均為λi，2，用于膨脹下8層土壤溫度（2.8 cm以下）；與土壤濕度對(duì)應(yīng)的10個(gè)元素值則為1.0，即 不 膨 脹 土 壤 濕 度。由 式 （5）可 知，－2Li（λi）實(shí)際上是關(guān)于Hi［λi］的函數(shù)，即無(wú)論膨脹因子矩陣［λi］中包含多少未知參數(shù)，通過極小化－2Li（λi）只能估計(jì)出Hi［λi］能包含的那些參數(shù)。如果被同化的觀測(cè)信息只有地表溫度，則Hi［λi］只包含用于膨脹地表溫度和上兩層土壤溫度的參數(shù)λi，1，即 MLE方法通過極小化－2Li（λi）只能估計(jì)出λi，1，無(wú)法估計(jì)出λi，2。為此，本文在同化系統(tǒng)中又引入10 cm土壤溫度觀測(cè)信息，這樣Hi［λi］中就同時(shí)包含了λi，1和λi，2，從而保證 MLE方法可以同時(shí)估計(jì)得到這兩個(gè)參數(shù)。

為探討對(duì)淺層和深層土壤溫度進(jìn)行不同程度膨脹的必要性，本文將比較淺層和深層土壤溫度使用同一系數(shù)和不同系數(shù)進(jìn)行膨脹時(shí)的結(jié)果。為便于表達(dá)，淺層和深層土壤溫度使用同一膨脹系數(shù)和不同膨脹系數(shù)的MLE方法分別簡(jiǎn)稱為MLE1方法和MLE2方法。

3 同化結(jié)果

3.1 確定觀測(cè)算子和觀測(cè)誤差

本研究同化的觀測(cè)變量為地表溫度和10 cm土壤溫度，狀態(tài)變量為地表溫度、10層土壤溫度和10層土壤濕度。

本研究同化的地表溫度觀測(cè)是站點(diǎn)觀測(cè)數(shù)據(jù)，觀測(cè)變量與模型中的地表溫度變量一致，故地表的觀測(cè)算子為1。本研究還同化了站點(diǎn)觀測(cè)的10 cm土壤溫度，而這一深度無(wú)法與Co LM中定義的任一層土壤節(jié)點(diǎn)的深度直接對(duì)應(yīng)。因此，利用相鄰兩層模擬深度上的土壤溫度經(jīng)線性內(nèi)插得到10 cm土壤溫度，該插值算法即為10 cm深度觀測(cè)算子：

其中，y1和y2分別為地表和10 cm土壤溫度的模擬值，d3和d4分別為Co LM中第3層和第4層土壤節(jié)點(diǎn)的深度（6.2 cm和11.9 cm），x0為Co LM 的地表溫度模擬值，x3和x4分別為Co LM中第3層和第4層的土壤溫度模擬值。該觀測(cè)算子將Co LM輸出的模型狀態(tài)變量與觀測(cè)變量聯(lián)系起來。

在DGS站，地表溫度使用4000.4G紅外溫度傳感器觀測(cè)，土壤溫度使用溫度計(jì)觀測(cè)。查閱相關(guān)儀器手冊(cè)，最終確定地表溫度和10 cm土壤溫度的觀測(cè)誤差分別為1.0 K和0.5 K。

3.2 同化結(jié)果

使用DGS站2002年10月1日—2003年9月30日的觀測(cè)數(shù)據(jù)進(jìn)行試驗(yàn)。Co LM從2002年10月1日積分到2003年5月29日，作為spin-up使模型達(dá)到平衡，2003年5月30日同化開始，積分步長(zhǎng)為1 h。基于En KF方法和MLE方法，Co LM每12 h同化1次地表溫度和10 cm土壤溫度觀測(cè)數(shù)據(jù)（08：00和20：00），集合樣本量為100。選用研究時(shí)段為2003年9月1—30日不同深度的土壤溫度來驗(yàn)證同化效果。

淺層土壤溫度的時(shí)間變化較快且變幅較大，很難區(qū)分其在30 d內(nèi)的觀測(cè)場(chǎng)、模擬場(chǎng)和同化場(chǎng)。因此，對(duì)地表溫度和3 cm土壤溫度在30 d內(nèi)的結(jié)果進(jìn)行平均，得到其日變化（圖1）。En KF方法和MLE方法對(duì)地表溫度（圖1a）和3 cm土壤溫度（圖1b）的估計(jì)都比模擬值準(zhǔn)確，尤其在09：00—23：00時(shí)段。與模擬相比，En KF方法在01：00—07：00有所低估，MLE方法在此時(shí)段內(nèi)的結(jié)果則可以與模擬場(chǎng)相當(dāng)甚至更優(yōu)。

圖1 2003年9月1—30日地表溫度（a）和3 cm土壤溫度（b）觀測(cè)場(chǎng)、模擬場(chǎng)和同化場(chǎng)平均日變化Fig.1 The diurnal variation of observed，simulated and assimilated soil temperature at 0 cm（a）and 3 cm（b）averaged from 1 Sep to 30 Sep in 2003

圖2為2003年9月1—30日10 cm和40 cm土壤溫度的觀測(cè)場(chǎng)、模擬場(chǎng)和同化場(chǎng)。在10 cm深度（圖2a），En KF方法和MLE方法均比模擬更準(zhǔn)確地估計(jì)出了土壤溫度日變化的趨勢(shì)和幅度，但MLE方法相比En KF方法的優(yōu)勢(shì)不大。在40 cm深度（圖2b），與觀測(cè)場(chǎng)相比，模擬場(chǎng)明顯低估，而兩種同化方法均有效改善了這種低估，且MLE方法同化得到的土壤溫度比En KF方法更接近觀測(cè)。100 cm的結(jié)果與40 cm相似（圖略）。

圖2 2003年9月1—30日10 cm（a）和40 cm（b）土壤溫度觀測(cè)場(chǎng)、模擬場(chǎng)和同化場(chǎng)變化Fig.2 The observed，simulated and assimilated soil temperature at 10 cm（a）and 40 cm（b）from 1 Sep to 30 Sep in 2003

表1為各層土壤溫度模擬場(chǎng)和同化場(chǎng)在2003年9月1—30日平均的均方根誤差。En KF方法和MLE方法同化得到的均方根誤差均低于模擬結(jié)果的均方根誤差，而且膨脹因子為對(duì)角矩陣的MLE2方法比En KF方法更好。表1還表明，同化對(duì)土壤溫度估計(jì)效果的改善在10 cm，40 cm和100 cm深度比0 cm和3 cm深度明顯。淺層的土壤溫度變化較快且變幅較大，因此，雖然同化時(shí)刻的同化結(jié)果會(huì)接近觀測(cè)值，但同化時(shí)刻后，在持續(xù)較長(zhǎng)的預(yù)報(bào)時(shí)間內(nèi)，變化很快的淺層土壤溫度會(huì)以較快的速度向模型模擬結(jié)果靠攏，導(dǎo)致上一步同化對(duì)淺層土壤溫度的改善在下一步同化到來時(shí)影響不大。這樣，同化時(shí)刻對(duì)淺層土壤溫度造成的影響就很難在整個(gè)預(yù)報(bào)過程中得到有效保持和傳遞。為進(jìn)一步分析該問題，本文單獨(dú)針對(duì)2003年9月1—30日每日08：00和20：00的同化時(shí)刻，統(tǒng)計(jì)了膨脹因子為對(duì)角矩陣的MLE2方法同化得到的地表溫度、3 cm和10 cm土壤溫度在同化更新前的背景場(chǎng)均方根誤差和同化更新后的分析場(chǎng)均方根誤差?？偟膩碚f，時(shí)間平均的地表溫度、3 cm和10 cm土壤溫度同化場(chǎng)均方根誤差（分別為0.402，1.920 K和0.157 K）均顯著低于背景場(chǎng)均方根誤差（分別為3.028，3.509 K和1.110 K），這表明同化更新后的同化場(chǎng)比背景場(chǎng)更加接近觀測(cè)。這與前面的分析一致：在每日兩次的同化時(shí)刻，同化的確對(duì)淺層土壤溫度產(chǎn)生了較大影響，但這種影響很難在整個(gè)預(yù)報(bào)過程中得到有效保持。

考慮到淺層和深層土壤溫度的性質(zhì)不同，本文在實(shí)施MLE方法時(shí)，對(duì)淺層和深層土壤溫度采用了不同的膨脹因子。從2003年9月1—30日平均結(jié)果看，MLE1方法（淺層和深層土壤溫度使用同一系數(shù)進(jìn)行膨脹）中地表溫度和10層土壤溫度的膨脹因子值均為2.2，而MLE2方法（淺層和深層土壤溫度使用不同系數(shù)進(jìn)行膨脹）中地表溫度和上兩層（0.7～2.8 cm）土壤溫度的膨脹因子值為3.2，下8層（6.2～286.5 cm）土壤溫度的膨脹因子值為1.3?？梢姡琈LE1方法和MLE2方法都通過膨脹因子而對(duì)預(yù)報(bào)誤差協(xié)方差矩陣的離散度發(fā)揮了作用。但MLE1方法沒有對(duì)淺層和深層土壤溫度區(qū)別膨脹，導(dǎo)致深層土壤溫度膨脹過度；而MLE2方法對(duì)淺層和深層土壤溫度采用了不同的膨脹因子，使得深層土壤溫度的膨脹相對(duì)合理。表1中分別列出了MLE1和MLE2方法得到的土壤溫度的均方根誤差?？梢钥吹?，MLE1方法對(duì)各層土壤溫度的估計(jì)均沒有MLE2方法準(zhǔn)確，尤其對(duì)于深層土壤溫度（40 cm和100 cm），MLE1方法同化得到的均方根誤差比MLE2方法高50%以上，甚至比En KF方法同化得到的均方根誤差還高。圖3為En KF方法、MLE2方法和MLE1方法同化得到的40 cm和100 cm土壤溫度。由圖3可見，由于MLE1方法對(duì)整個(gè)土壤溫度廓線都采用同一因子進(jìn)行膨脹，將適用于變化很快的淺層土壤溫度的膨脹因子，也用于膨脹變化較慢的深層土壤溫度，造成深層土壤溫度被不合理膨脹，進(jìn)而影響了深層的同化效果。

表1 2003年9月1—30日土壤溫度模擬場(chǎng)和同化場(chǎng)平均的均方根誤差（單位：K）Table 1 Root mean square error of the simulated and assimilated soil temperature from 1 Sep to 30 Sep in 2003（unit：K）

圖3 MLE1方法和MLE2方法分別同化得到2003年9月1—30日40 cm（a）和100 cm（b）土壤溫度Fig.3 Soil temperature at 40 cm（a）and 100 cm（b）assimilated by MLE1 and MLE2 from 1 Sep to 30 Sep in 2003

圖4為2003年9月1—30日3 cm近地表土壤濕度的觀測(cè)場(chǎng)、模擬場(chǎng)和同化場(chǎng)。與觀測(cè)相比，Co LM的模擬普遍高估了土壤濕度，而En KF方法和MLE方法均對(duì)此有所緩解，且MLE方法的改善更明顯。這表明，雖然同化過程中沒有引入土壤濕度的觀測(cè)信息，同化過程仍可以通過誤差協(xié)方差矩陣和模式傳輸而影響土壤濕度。

圖4 2003年9月1—30日3 cm土壤濕度觀測(cè)場(chǎng)、模擬場(chǎng)和同化場(chǎng)Fig.4 The observed，simulated and assimilated soil moisture at 3 cm from 1 Sep to 30 Sep in 2003

4 結(jié)論和討論

本文基于極大似然估計(jì)理論，發(fā)展了一種優(yōu)化預(yù)報(bào)誤差協(xié)方差的實(shí)時(shí)膨脹方法——MLE方法。當(dāng)然，除了協(xié)方差膨脹方法外，國(guó)內(nèi)外在預(yù)報(bào)誤差協(xié)方差的分析和優(yōu)化方法上還有很多其他研究［27-30］，這些超出了本文目前的研究范圍，有待于下一步的跟蹤和關(guān)注。本研究表明：

1）總的來說，MLE方法對(duì)地表溫度和各層土壤溫度的估計(jì)都比En KF方法準(zhǔn)確，對(duì)深層土壤溫度（40 cm和100 cm）尤其如此。

2）考慮到淺層和深層土壤溫度的差別，本文對(duì)淺層和深層土壤溫度采用不同的膨脹因子，即MLE2方法。與膨脹因子為單一標(biāo)量值時(shí)的MLE1方法比較表明，由于對(duì)整個(gè)土壤溫度廓線采用相同的膨脹因子，MLE1方法對(duì)各層（尤其深層）土壤溫度的估計(jì)沒有MLE2方法準(zhǔn)確。這表明，多因子膨脹可以緩解深層土壤溫度的不合理膨脹問題。

3）為分析葉面積指數(shù)對(duì)土壤溫度的影響，本文在Co LM中引入MODIS LAI產(chǎn)品代替經(jīng)驗(yàn)算法計(jì)算的葉面積指數(shù)。結(jié)果表明，與使用經(jīng)驗(yàn)葉面積指數(shù)相比，使用MODIS LAI可以改善Co LM對(duì)地表溫度和3 cm土壤溫度的模擬，對(duì)地表溫度和各層土壤溫度的同化也更加準(zhǔn)確。

本文在估算模型誤差協(xié)方差矩陣時(shí)使用對(duì)角矩陣（協(xié)方差項(xiàng)為0），沒有考慮不同土壤層變量之間的相關(guān)。Jin等［31］研究表明，當(dāng)合理給定模型誤差協(xié)方差矩陣中的協(xié)方差項(xiàng)時(shí)，同化系統(tǒng)能將表層優(yōu)化后的信息快速傳遞給深層，從而快速顯著改善整個(gè)土壤廓線狀態(tài)變量的估計(jì)。但其工作是單獨(dú)優(yōu)化土壤溫度或濕度，未涉及土壤溫度與濕度之間的誤差協(xié)相關(guān)。在同步優(yōu)化土壤溫度和土壤濕度時(shí)，如何確定土壤溫度和濕度之間的協(xié)相關(guān)，值得研究。

本文膨脹因子隨時(shí)間變化，但在空間上不變，這實(shí)際上不合理，尤其是觀測(cè)資料在空間上分布不均勻時(shí)。適用于觀測(cè)資料濃密區(qū)的膨脹因子，若持續(xù)用于所有地區(qū)，可能導(dǎo)致觀測(cè)資料稀疏區(qū)過度膨脹。為對(duì)淺層和深層土壤溫度進(jìn)行不同程度的膨脹，本文引入了10 cm土壤溫度的觀測(cè)信息。如果隨時(shí)空變化的膨脹方法得以實(shí)現(xiàn)，則僅同化地表溫度就能針對(duì)不同深度的土壤溫度估計(jì)出不同的膨脹因子。Anderson［32］提出的方法中膨脹因子時(shí)空上均可變化，但其工作只考慮了觀測(cè)誤差空間獨(dú)立的情況，不利于遙感資料在同化中的應(yīng)用。因此，有必要發(fā)展適用于空間相關(guān)觀測(cè)系統(tǒng)的、時(shí)空變化的膨脹方法。

［1］ Evensen G.Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte-Carlo methods to forecast error statistics.J Geophys Res，1994，99（C5）：10143-10162.

［2］ Burgers G，Leeuwen P J V，Evensen G.Analysis scheme in the Ensemble Kalman Filter.Mon Wea Rev，1998，126：1719-1724.

［3］ Evensen G.The ensemble Kalman filter：Theoretical formulation and practical implementation.Ocean Dynam，2003，53（4）：343-367.

［4］ Senegas J，Wackernagel H H，Rosenthal W，et al.Error covariance modeling in sequential data assimilation.Stoch Env Res Risk A，2001，15：65-86.

［5］ Daley R.Atmospheric Data Analysis.Cambridge：Cambridge University Press，1991.

［6］ Kalnay E.Atmospheric Modeling，Data Assimilation，and Predictability.Cambridge：Cambridge University Press，2002.

［7］ Anderson J L，Anderson S L.A Monte Carlo implementation of the non-linear filtering problem to produce ensemble assimilations and forecasts.Mon Wea Rev，1999，127：2741-2758.

［8］ Constantinescu E M，Sandu A，Chai T，et al.Ensemble-based chemical data assimilation.I：General approach.Q J R Meteorol Soc，2007，133：1229-1243.

［9］ Dee D P.On-line estimation of error covariance parameters for atmospheric data assimilation.Mon Wea Rev，1995，123（4）：1128-1145.

［10］ Dee D P，da Silva A M.Maximum-likelihood estimation of forecast and observation error covariance parameters.PartⅠ：Methodology.Mon Wea Rev，1999，127：1822-1834.

［11］ Dee D P，Gaspari G，Redder C，et al.Maximum-likelihood estimation of forecast and observation error covariance parameters.PartⅡ：Applications.Mon Wea Rev，1999，127：1835-1849.

［12］ Wang X，Bishop C H.A comparison of breeding and ensemble transform Kalman filter ensemble forecast schemes.J Atmos Sci，2003，60：1140-1158.

［13］ Li H，Kalnay E，Miyoshi T.Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter.Q J R Meteorol Soc，2009，135（639）：523-533.

［14］ Miyoshi T.The gaussian approach to adaptive covariance inflation and its implementation with the Local Ensemble Transform Kalman Filter.Mon Wea Rev，2011，139：1519-1535.

［15］ Zheng X G.An adaptive estimation of forecast error covariance parameters for Kalman filtering data assimilation.Adv Atmos Sci，2009，26：154-160.

［16］ Liang X，Zheng X G，Zhang S P，et al.Maximum likelihood es-timation of inflation factors on error covariance matrices for ensemble Kalman filter assimilation.Q J R Meteorol Soc，2011，138：263-273，doi：10.1002／qj.912.

［17］ Huang C L，Li X，Lu L.Retrieving soil temperature profile by assimilating MODIS LST products with ensemble Kalman filter.Rem Sens Environ，2008，112：1320-1336.

［18］ Yang K，Koike T，Kaihotsu I，et al.Validation of a dual-pass microwave land data assimilation system for estimating surface soil moisture in semiarid regions.J Hydrometeorology，2009，10（3）：780-793.

［19］ 楊曉峰，陸其峰，楊忠東.基于AMSR-E土壤濕度產(chǎn)品的LIS同化試驗(yàn).應(yīng)用氣象學(xué)報(bào)，2013，24（4）：435-445.

［20］ 吳統(tǒng)文，宋連春，劉向文，等.國(guó)家氣候中心短期氣候預(yù)測(cè)模式系統(tǒng)業(yè)務(wù)化進(jìn)展.應(yīng)用氣象學(xué)報(bào)，2013，24（5）：533-543.

［21］ 王莉，黃嘉佑.Kalman濾波的試驗(yàn)應(yīng)用研究.應(yīng)用氣象學(xué)報(bào)，1999，10（3）：276-282.

［22］ 趙曉琳，朱國(guó)富，李澤椿.基于TIGGE資料識(shí)別適應(yīng)性觀測(cè)敏感區(qū)的應(yīng)用研究.應(yīng)用氣象學(xué)報(bào)，2010，21（4）：405-415.

［23］ Whitaker J S，Hamill T H.Ensemble data assimilation without perturbed observations.Mon Wea Rev，2002，130：1913-1924.

［24］ Dai Y J，Zeng X B，Dickinson R E，et al.The common land model.Bull Amer Meteor Soc，2003，84：1013-1023.

［25］ 孟春雷，張朝林.路面氣象數(shù)值預(yù)報(bào)模型及性能檢驗(yàn).應(yīng)用氣象學(xué)報(bào)，2012，23（4）：451-458.

［26］ Koike T.Coordinated Enhanced Observing Period （CEOP）-An initial step for integrated global water cycle observation.World Meteorological Organization Bulletin，2004，53（2）：115-121.

［27］ 龔建東，趙剛.全球資料同化中誤差協(xié)方差三維結(jié)構(gòu)的準(zhǔn)確估計(jì)與應(yīng)用：背景誤差協(xié)方差調(diào)整與數(shù)值試驗(yàn)分析.氣象學(xué)報(bào)，2006，64（6）：669-682.

［28］ 曹小群，黃思訓(xùn)，張衛(wèi)民，等.區(qū)域三維變分同化中背景誤差協(xié)方差的模擬.氣象科學(xué)，2008，28（1）：8-14.

［29］ 馬旭林，莊照榮，薛紀(jì)善，等.GRAPES非靜力數(shù)值預(yù)報(bào)模式的三維變分資料同化系統(tǒng)的發(fā)展.氣象學(xué)報(bào)，2009，67（1）：50-60.

［30］ 王曼，李華宏，段旭，等.WRF模式三維變分中背景誤差協(xié)方差估計(jì).應(yīng)用氣象學(xué)報(bào)，2011，22（4）：482-492.

［31］ Jin R，Li X.Improving the estimation of hydrothermal state variables in the active layer of frozen ground by assimilating in situ observations and SSM／I data.Sci China Ser D：Earth Sci，2009，39（9）：1220-1231.

［32］ Anderson J L.Spatially and temporally varying adaptive covariance inflation for ensemble filters.Tellus，2009，61：72-83.

A Method of Improving Error Covariances in En KF and Its Application to Data Assimilation

Liang Xiao1）Zheng Xiaogu2）Dai Yongjiu2）Shi Chunxiang1）

1）（National Meteorological Information Center，Beijing100081）
2）（Beijing Normal University，Beijing100875）

In the ensemble Kalman filter（En KF），the forecast error covariance matrix is estimated as the sampling covariance matrix of the forecast ensemble.However，previous studies suggest that the sampling error resulting from finite-size ensembles may make such estimations far from the true forecast error covariance，and finally degrade the performance of En KF.A common way to address this problem is covariance inflation with a time-constant inflation factor.A time-dependent infiation approach on forecast error covariance matrix（i.e.，MLE method）is developed based on the maximum likelihood estimation theory，so as to improve estimates of forecast error covariances.At Delgertsgot（DGS）Station in the Mongolian Plateau reference site，point observations of ground temperature and soil temperature at the depth of 10 cm are assimilated into the Common Land Model（Co LM）with a frequency of every 12 hours，using two assimilation algorithms（En KF method and MLE method），in order to test the effectivity of MLE in practical assimilation.In this way，a soil temperature assimilation system is constructed on the point scale.

Results indicate that MLE method performs better than En KF method for ground temperature and soil temperatures at most depths（especially for soil temperatures at deeper depths）.Moreover，considering differences between soil temperatures at shallower depths and those at deeper depths，different inflation factors are adopted to them when implementing MLE method.Compared to results of MLE using a single scalar inflation factor，it shows that multiple-factor inflation is able to alleviate the unreasonable inflation of soil temperatures at deeper depths and therefore get better assimilation results.In addition，the leaf area index（LAI）in the Co LM is updated dynamically by MODIS LAI products，and results derived using MODIS LAI are compared to those derived using LAI computed by experiential formula，so as to study the effect of the LAI accuracy on simulated and assimilated soil temperatures.It shows that using MODIS LAI can get better simulation of soil temperature at depths of 0 cm and 3 cm，as well as more accurate assimilation of soil temperature at most depths.

The inflation factor is set to be variable in time，but constant in space.However，variables such as soil temperature and soil moisture behave quite differently at shallow surfaces and deep depths，and observations may be unevenly distributed in space in regional assimilation researches.Therefore，it is necessary to adopt different inflation factors to different variables or to the same variable at different locations.In the future，it is necessary to develop a time-and-space dependent inflation method and test its capability in real applications.

data assimilation；En KF；error covariance inflation

梁曉，鄭小谷，戴永久，等.En KF中誤差協(xié)方差優(yōu)化方法及在資料同化中應(yīng)用.應(yīng)用氣象學(xué)報(bào)，2014，25（4）：397-405.

2014-01-06收到，2014-05-05收到再改稿。

公益性行業(yè)（氣象）科研專項(xiàng)（GYHY201206013，GYHY201306045），國(guó)家國(guó)際科技合作專項(xiàng)（2011DFG23150）

*email：liangx＠cma.gov.cn