基于EnKF 和ETKF 的比較研究

陳 麗，卞康旭

（中國民用航空飛行學(xué)院理學(xué)院,四川 廣漢 618307）

1 前言

預(yù)測(cè)物理系統(tǒng)通常需要系統(tǒng)的時(shí)間演變模型和系統(tǒng)當(dāng)前狀態(tài)的估計(jì)。在一些應(yīng)用中，可以高精度地直接測(cè)量系統(tǒng)的狀態(tài)。在其他應(yīng)用中，如天氣預(yù)報(bào)，直接測(cè)量全球系統(tǒng)狀態(tài)是不可行的。相反，必須從可用數(shù)據(jù)中推斷出狀態(tài)。雖然基于當(dāng)前數(shù)據(jù)的狀態(tài)估計(jì)是合理的，但通?？梢酝ㄟ^使用當(dāng)前和過去的數(shù)據(jù)來獲得更好的估計(jì)[1]?！皵?shù)據(jù)同化”在持續(xù)的基礎(chǔ)上提供這種估計(jì)，在預(yù)測(cè)步驟和狀態(tài)估計(jì)步驟之間迭代交替；后一個(gè)步驟通常被稱為“分析”。分析步驟結(jié)合來自當(dāng)前數(shù)據(jù)和先前短期預(yù)測(cè)（基于過去數(shù)據(jù)）的信息，產(chǎn)生當(dāng)前狀態(tài)估計(jì)。該估計(jì)用于初始化下一個(gè)短期預(yù)測(cè)，隨后用于下一次分析等[2]。

同化方法主要有兩種類型：變分同化方法和序列同化方法。序列數(shù)據(jù)同化方案顯式求解一系列方程，以找到系統(tǒng)的分析狀態(tài)。順序方法的例子是卡爾曼濾波器（KF）[3]以及從KF 的基礎(chǔ)上導(dǎo)出的各種濾波器。例如，集合卡爾曼濾波器（EnKF）[4]、集成變換卡爾曼濾波器（ETKF）[5]都是順序方案。

集合卡爾曼濾波器由Geir Evensen 于1994 年提出，并于1998 年進(jìn)行了修正[6]。它的半經(jīng)驗(yàn)性論證可能令人不安。盡管如此，集合卡爾曼濾波器已被證明在大量的學(xué)術(shù)和操作數(shù)據(jù)同化問題上是非常有效的。集合卡爾曼濾波器（EnKF)是一個(gè)降階卡爾曼濾波器，它只處理二階以下的誤差統(tǒng)計(jì)信息。由于其在計(jì)算過程中對(duì)引入的觀測(cè)信息進(jìn)行了隨機(jī)擾動(dòng)，又稱為隨機(jī)性集合卡爾曼濾波。

隨機(jī)EnKF[7]能夠單獨(dú)跟蹤集合中的每個(gè)成員，并通過分析步驟使它們相互作用。這個(gè)好的性質(zhì)是有代價(jià)的：其需要獨(dú)立地?cái)_動(dòng)每個(gè)成員的觀測(cè)向量。雖然這看起來易于實(shí)現(xiàn)，但在隨機(jī)擾動(dòng)時(shí)也引入了數(shù)值噪聲。這可能會(huì)影響算法的性能，特別是當(dāng)單個(gè)分析的觀測(cè)數(shù)量有限時(shí)，可能會(huì)造成算法發(fā)散。為了解決這個(gè)問題，發(fā)展出了另一種類型的濾波算法，稱為確定性集合卡爾曼濾波。在本文中，將重點(diǎn)關(guān)注確定性EnKF 的集合變換變體，通?？s寫為ETKF。該算法不是在狀態(tài)空間或觀測(cè)空間中執(zhí)行計(jì)算，而是在集合空間中執(zhí)行。

目前，對(duì)于隨機(jī)性和確定性集合卡爾曼濾波算法的比較研究較少，本文將基于不同實(shí)驗(yàn)?zāi)Ｐ捅容^隨機(jī)EnKF 與ETKF 的算法性能。

2 研究方法

2.1 隨機(jī)性集合卡爾曼濾波（EnKF)

2.1.1 分析步驟

在EnKF 算法框架中[8]，先從狀態(tài)集合開始分析，記狀態(tài)集合為：

將觀測(cè)向量進(jìn)行擾動(dòng)yi=y+ui，其中ui來自高斯分布ui～N(0,R)。首先，確保以避免偏差。其次，定義經(jīng)驗(yàn)誤差協(xié)方差矩陣：

當(dāng)m→∞時(shí)，Ru的極限應(yīng)趨于R。

其中：

那么分析誤差協(xié)方差矩陣為：

2.1.2 預(yù)報(bào)步驟

在預(yù)測(cè)步驟中，在分析中獲得的更新集合通過一個(gè)時(shí)間步長由模型傳播到下一個(gè)瞬間：

預(yù)報(bào)估計(jì)是預(yù)報(bào)集合的平均值：

而預(yù)測(cè)誤差協(xié)方差矩陣則由：

重要的是要注意，在算法過程中已經(jīng)避免使用切線性算子，或任何線性化。

2.2 集合變換卡爾曼濾波(ETKF)

ETKF 與EnKF 的不同之處在于分析集合的構(gòu)造。在下面的討論中，定義集合空間如下：

對(duì)于分析估計(jì)，利用卡爾曼增益的形式得到均值的最優(yōu)系數(shù)向量wa：

為了生成一個(gè)后驗(yàn)集合，下面需要分解。

2.3 乘法膨脹

為了防止濾波發(fā)散，引入乘法膨脹，其簡(jiǎn)單地將集合誤差協(xié)方差Pf進(jìn)行處理，以近似真實(shí)誤差協(xié)方差Pf：Pf←λPf，其中λ是膨脹因子。

3 數(shù)值實(shí)驗(yàn)

本節(jié)提供了EnKF、ETKF 在三種模型上的一些實(shí)驗(yàn)說明，一共進(jìn)行了6 個(gè)實(shí)驗(yàn)，從LA 模型的2 個(gè)實(shí)驗(yàn)開始，其次是Lorenz3 模型的2 個(gè)實(shí)驗(yàn)，最后是Lorenz96 模型的2個(gè)實(shí)驗(yàn)。這些實(shí)驗(yàn)的設(shè)置與文獻(xiàn)[9]中的設(shè)置相同。對(duì)于實(shí)驗(yàn)性能用以下“分析rmse”進(jìn)行考量：

其中，x和x真分別為待估計(jì)的模型狀態(tài)和真實(shí)的模型狀態(tài)。

3.1 LA 模型的實(shí)驗(yàn)

在本小節(jié)中，基于線性平流模型（LA）進(jìn)行了2個(gè)實(shí)驗(yàn)，線性平流模型表達(dá)如下：

狀態(tài)樣本xi由具有隨機(jī)幅度和相位的nk正弦波組成，表示為：

rand()函數(shù)返回一個(gè)均勻分布的隨機(jī)數(shù)，0≦rand()≦1,exp()是e的指數(shù)函數(shù)。實(shí)驗(yàn)基本設(shè)置都是相同的。唯一的區(qū)別在于觀測(cè)誤差。系統(tǒng)運(yùn)行700 個(gè)同化周期，前100個(gè)循環(huán)被丟棄。超過600個(gè)循環(huán)的平均結(jié)果見表1。

表1 LA 模型

從表1 中可知，對(duì)于線性模型，不論觀測(cè)誤差大小，ETKF 算法性能都優(yōu)于EnKF。

3.2 Lorenz3 模型的實(shí)驗(yàn)

在本小節(jié)中，基于Lorenz3 模型進(jìn)行了2 個(gè)實(shí)驗(yàn)，Lorenz3 模型方程如下所示：

其中為σ=10、ρ=28 和β=8/3。該系統(tǒng)用四階龍格-庫塔方案積分，采用固定的時(shí)間步長為0.01。

與LA 模型一樣，系統(tǒng)運(yùn)行700 個(gè)同化周期，前100 個(gè)循環(huán)被丟棄。超過600 個(gè)循環(huán)的平均結(jié)果見表2。

表2 Lorenz3 模型

從表2 中可知，在Lorenz3 模型中，觀測(cè)誤差較小時(shí)，ETKF 算法性能優(yōu)于EnKF,但觀測(cè)誤差較大時(shí)，EnKF 算法性能卻優(yōu)于ETKF。

3.3 Lorenz96 模型的實(shí)驗(yàn)

Lorenz96 表達(dá)式為：

在Lorenz96 模型中，依然進(jìn)行了2 個(gè)實(shí)驗(yàn)。系統(tǒng)依然運(yùn)行700 個(gè)同化周期，前100 個(gè)循環(huán)被丟棄。超過600 個(gè)循環(huán)的平均結(jié)果見表3。

表3 Lorenz96 模型

如表3 所示，在Lorenz96 模型中，和LA 模型結(jié)果相同，不論觀測(cè)誤差大小，ETKF 算法性能都優(yōu)于EnKF，不過EFKF 同化時(shí)間始終高于EnKF。

4 結(jié)論

隨機(jī)性集合卡爾曼濾波與集合變換卡爾曼濾波都是經(jīng)典的數(shù)據(jù)同化算法，在其他實(shí)驗(yàn)條件相同的情形下，兩者算法的性能表現(xiàn)出些許不同。對(duì)于線性模型（如LA 模型）和強(qiáng)非線性模型（如Lorenz96 模型）來說，不論觀測(cè)誤差大小，ETKF 的算法性能始終優(yōu)于EnKF，但其同化時(shí)間比EnKF 較多，不過在算法性能的優(yōu)勢(shì)下，時(shí)間成本可忽略不計(jì)。對(duì)于弱非線性模型（如Lorenz3模型），觀測(cè)誤差較小時(shí)，ETKF 算法性能優(yōu)于EnKF，但觀測(cè)誤差較大時(shí)，EnKF 算法性能卻優(yōu)于ETKF。