集合變換卡爾曼濾波數(shù)據(jù)同化方法中的協(xié)方差局地化

王際朝， 王 玥， 臧紹東， 楊俊鋼， 紀(jì)艷菊， 阮宗利

(1. 中國(guó)石油大學(xué)(華東) 理學(xué)院， 山東 青島 266580； 2. 自然資源部第一海洋研究所， 山東 青島 266061)

數(shù)據(jù)同化是通過(guò)對(duì)觀測(cè)數(shù)據(jù)與預(yù)測(cè)值進(jìn)行分析，尋找一個(gè)當(dāng)前物理狀態(tài)的最優(yōu)估計(jì)值(分析值)， 作為數(shù)值預(yù)測(cè)模型初始條件的方法。作為連接數(shù)值模式和觀測(cè)數(shù)據(jù)的技術(shù)手段， 數(shù)據(jù)同化在海洋科學(xué)領(lǐng)域得到了廣泛關(guān)注。當(dāng)前， 集合數(shù)據(jù)同化方法得以快速發(fā)展， 已經(jīng)成為業(yè)務(wù)化數(shù)值預(yù)測(cè)的一個(gè)可行選擇。集 合卡爾 曼 濾波(ensemble Kalman filter， EnKF)[1]，可以將卡爾曼濾波應(yīng)用到強(qiáng)非線性模型， 擴(kuò)展了卡爾曼濾波的適用性， 是數(shù)據(jù)同化領(lǐng)域的研究熱點(diǎn)之一。EnKF中的預(yù)測(cè)誤差協(xié)方差矩陣通過(guò)對(duì)預(yù)測(cè)集合成員進(jìn)行統(tǒng)計(jì)計(jì)算得到， 預(yù)測(cè)誤差對(duì)于集合數(shù)據(jù)同化方法而言至關(guān)重要。然而， 產(chǎn)生足夠大的預(yù)測(cè)集合的計(jì)算成本令人望而卻步。當(dāng)使用的集合數(shù)目過(guò)小時(shí)， 則不能在統(tǒng)計(jì)意義上充分表示模型的變化， 從而會(huì)引起欠采樣(under-sampling)。通常， 欠采樣會(huì)導(dǎo)致濾波發(fā)散、預(yù)測(cè)誤差協(xié)方差低估和虛假相關(guān)等負(fù)面問(wèn)題。由于這些問(wèn)題的影響， EnKF可能會(huì)產(chǎn)生一個(gè)次優(yōu)的分析結(jié)果。

在進(jìn)行數(shù)據(jù)同化時(shí)， 為了克服欠采樣帶來(lái)的負(fù)面影響， 有關(guān)學(xué)者已經(jīng)提出了很多解決方法。當(dāng)前基于集合的數(shù)據(jù)同化研究中， 最著名的解決方法是協(xié)方差膨脹(covariance inflation， CI)[2]、協(xié)方差局地化(covariance localization， CL)[3-4]和 局 地 化 分 析(local analysis， LA)[5-6]。CI方法通過(guò)一個(gè)膨脹因子校正預(yù)測(cè)誤差協(xié)方差的低估問(wèn)題， 主要思想是通過(guò)膨脹集合均值和集合成員之間的誤差來(lái)增大預(yù)測(cè)誤差協(xié)方差。CI方法用于分析過(guò)程之前， 具體表示為：其中“←”表示對(duì)于之前狀態(tài)值的替代，r是所謂的膨脹因子。盡管CI方法克服了預(yù)測(cè)誤差協(xié)方差的低估問(wèn)題， 但是協(xié)方差膨脹因子r解決不了虛假相關(guān)問(wèn)題。CL方法在預(yù)測(cè)誤差協(xié)方差矩陣和局地化函數(shù)之間實(shí)施舒爾積運(yùn)算[7]， 通過(guò)截?cái)囝A(yù)測(cè)誤差協(xié)方差矩陣中預(yù)先指定距離之外的誤差相關(guān)性以解決虛假相關(guān)問(wèn)題。此外， 舒爾積運(yùn)算能夠提高預(yù)測(cè)誤差協(xié)方差矩陣的秩[8]。LA方法以待更新?tīng)顟B(tài)變量為中心， 建立一個(gè)虛擬的局部窗口， 得到該狀態(tài)變量預(yù)測(cè)誤差協(xié)方差的局部近似值， 在分析更新過(guò)程中， 將局部窗口外的集合擾動(dòng)設(shè)定為零。例如， 通常可以以某個(gè)狀態(tài)變量為中心， 僅僅同化與其固定距離為M之內(nèi)的觀測(cè)數(shù)據(jù)， 固定距離的長(zhǎng)度決定了同化過(guò)程中觀測(cè)值的數(shù)量。顯然， 分別對(duì)狀態(tài)變量進(jìn)行局部同化的特點(diǎn)使LA方法可以利用并行方法進(jìn)行快速計(jì)算。

當(dāng)前， 越來(lái)越多的學(xué)者[9-11]關(guān)注的是CL方法和LA方法。相比于CL方法， LA方法是一種獨(dú)立的方法， 它可以用于任何數(shù)據(jù)同化的框架， 但是它的實(shí)際同化效果并不如CL方法。Miyoshi等[12]指出LA方法的同化效果和CL方法類似， 但是LA方法通常會(huì)導(dǎo)致弱局地化影響。之后， Janji?等[9]通過(guò)對(duì)Lorenz96模型進(jìn)行實(shí)驗(yàn)發(fā)現(xiàn)： 如果觀測(cè)誤差方差小于初始狀態(tài)的預(yù)測(cè)誤差方差， CL方法將會(huì)導(dǎo)致比較小的估計(jì)誤差； 如果觀測(cè)誤差方差大于初始狀態(tài)的預(yù)測(cè)誤差方差， CL方法和LA方法有相同的局地化影響。因此， 從之前的研究和實(shí)驗(yàn)可以得出， CL方法可能是相對(duì)較好的一個(gè)局地化方法。由于CL方法中的舒爾積運(yùn)算僅僅用于預(yù)測(cè)誤差協(xié)方差矩陣，而集合轉(zhuǎn)換卡爾曼濾波(ensemble transform Kalman filter， ETKF)方法中的預(yù)測(cè)誤差協(xié)方差矩陣并沒(méi)有顯式計(jì)算， 而是傳遞預(yù)測(cè)集合擾動(dòng)矩陣， 這導(dǎo)致ETKF方法中進(jìn)行舒爾積運(yùn)算的矩陣維度不一致，所以CL方法不能用于ETKF方法[13]。Janji?等[9]的文章給出了CL方法不適用于ETKF方法的具體討論。

為了克服CL方法不能用于ETKF方法的限制，在一元模型中， 通過(guò)對(duì)局地化函數(shù)的平方根矩陣和預(yù)測(cè)集合擾動(dòng)矩陣進(jìn)行舒爾積運(yùn)算， Petrie[14]提出了一種用于ETKF方法的CL近似方法。對(duì)于矩陣維度不一致的問(wèn)題， 通過(guò)在預(yù)測(cè)集合擾動(dòng)矩陣中添加n-N列零向量(n表示狀態(tài)變量個(gè)數(shù)，N表示集合個(gè)數(shù))， 對(duì)矩陣進(jìn)行擴(kuò)展來(lái)修正。然而， 韓培等[15]對(duì)此方法進(jìn)行探究， 表明這種近似的CL方法是一種弱近似并產(chǎn)生了較差的同化效果。

基于集合樣本相關(guān)性和矩陣因式分解， 近幾年提出了動(dòng)態(tài)計(jì)算局地化函數(shù)的方法， 使得CL方法可以應(yīng)用到ETKF方法中。Bishop等[16]提出了一種生成局地化函數(shù)的新方法， 該函數(shù)能夠隨真實(shí)誤差相關(guān)函數(shù)移動(dòng)并且還適應(yīng)真實(shí)誤差相關(guān)函數(shù)的寬度。但是此方法的計(jì)算代價(jià)比較大， 鑒于此， Bishop等[17]提出了一種新的協(xié)方差自適應(yīng)局地化方法(CALECO)的具體實(shí)施方案， 成功地將CL方法應(yīng)用到ETKF方法， 但是此方法仍然具有較高的計(jì)算復(fù)雜度。基于上述的研究基礎(chǔ)， Bishop等[18]在ETKF方法中提出了較為可行的CL方法的實(shí)施方案， 并命名為增益集合轉(zhuǎn)換卡爾曼濾波(gain form of ensemble transform Kalman filter， GETKF)。GETKF中使用的是Gaspari等[19]提出的基于距離相關(guān)的一個(gè)局地裁剪函數(shù)作為局地化函數(shù)(GC局地化函數(shù))， 此方法解決了在一元數(shù)值預(yù)測(cè)模型中CL方法不能用于ETKF方法的問(wèn)題。

當(dāng)前， 由于CL方法的實(shí)現(xiàn)存在限制， 國(guó)內(nèi)外對(duì)于CL方法在ETKF等平方根濾波中的研究較少。本文將基于GETKF方法進(jìn)行研究， 提出一種局地化效果較好的局地化方法， 并通過(guò)數(shù)值模擬實(shí)驗(yàn)比較本文基于GETKF修改的方法與GETKF方法的同化效果。

1 理論背景與方法改進(jìn)

Campbell等[20]比較了EnKF方法中分別使用模型空間CL方法和觀測(cè)空間CL方法的同化效果， 認(rèn)為模型空間CL方法通常優(yōu)于觀測(cè)空間CL方法。因此， 本文以下的研究基于模型空間CL方法。本節(jié)首先介紹GETKF方法， 然后對(duì)該方法進(jìn)行改進(jìn)。

1.1 增益集合轉(zhuǎn)換卡爾曼濾波

在ETKF方法中， 應(yīng)用CL方法的難點(diǎn)在于預(yù)測(cè)誤差協(xié)方差矩陣并沒(méi)有顯式表示， 而是通過(guò)預(yù)測(cè)集合擾動(dòng)矩陣隱式表達(dá)， 因而在預(yù)測(cè)集合擾動(dòng)矩陣和局地化函數(shù)之間無(wú)法進(jìn)行舒爾積運(yùn)算， 當(dāng)前主要的目的是在預(yù)測(cè)集合擾動(dòng)矩陣和局地化函數(shù)的平方根矩陣中進(jìn)行舒爾積運(yùn)算來(lái)近似預(yù)測(cè)誤差協(xié)方差矩陣和局地化函數(shù)之間的舒爾積運(yùn)算。假定預(yù)測(cè)集合擾動(dòng)矩陣表示為

局地化函數(shù)ρ的平方根矩陣表示為W， 即

其中， 平方根矩陣W通過(guò)對(duì)局地化函數(shù)ρ進(jìn)行特征值分解得到， 并且按特征值大小降序排列， 取前10個(gè)特征值和特征向量對(duì)構(gòu)成， 即

一般情況， CL方法需要進(jìn)行如下形式舒爾積運(yùn)算，

由于之前討論的ETKF方法的局限性， GETKF方法中采用了集合擴(kuò)展技術(shù)定義一個(gè)n× (N×L)維度的矩陣f

Z，

定義矩陣f Z表示局地化后的預(yù)測(cè)誤差協(xié)方差矩陣的平方根矩陣， 即

Bishop等[18]已經(jīng)證明式(6)是成立的。因此成功地在預(yù)測(cè)集合擾動(dòng)矩陣實(shí)現(xiàn)了CL方法?，F(xiàn)在使用代替作為當(dāng)前的預(yù)測(cè)集合擾動(dòng)矩陣， 即局地化后的預(yù)測(cè)誤差協(xié)方差矩陣的平方根矩陣， 重新構(gòu)造集合擴(kuò)展后的預(yù)測(cè)集合和預(yù)測(cè)誤差協(xié)方差矩陣。

令M=N×L， 則擴(kuò)展預(yù)測(cè)集合表示為其中每個(gè)列向量如下定義：

其中，zk表示矩陣的第k列。很顯然擴(kuò)展預(yù)測(cè)集合的均值與原始預(yù)測(cè)集合均值相同。擴(kuò)展預(yù)測(cè)集合誤差協(xié)方差矩陣表示為：

GETKF方法的分析過(guò)程與標(biāo)準(zhǔn)的ETKF方法的分析過(guò)程類似， 區(qū)別是GETKF方法使用擴(kuò)展預(yù)測(cè)集合以及擴(kuò)展預(yù)測(cè)集合擾動(dòng)矩陣代替原始預(yù)測(cè)集合和原始預(yù)測(cè)集合擾動(dòng)矩陣定義標(biāo)準(zhǔn)化的預(yù)測(cè)-觀測(cè)集合擾動(dòng)矩陣如公式(7)所示：

其中，

矩陣H表示觀測(cè)算子。為了提高ETKF方法同化的計(jì)算效率， 對(duì)標(biāo)準(zhǔn)化的預(yù)測(cè)-觀測(cè)集合擾動(dòng)矩陣進(jìn)行奇異值分解：

其中， 矩陣U和V是正交矩陣，的奇異值由矩陣∑的對(duì)角元素給出。接下來(lái)， 分析誤差協(xié)方差矩陣的平方根矩陣(分析集合擾動(dòng)矩陣)表示為，

其中，

以上展示了GETKF方法的主要計(jì)算過(guò)程， 可以看出此方法存在改進(jìn)的空間： 首先， GETKF方法利用集合擴(kuò)展技術(shù)在分析過(guò)程中產(chǎn)生M個(gè)集合(N＜M)， 然后使用集合擴(kuò)展之后的分析誤差協(xié)方差矩陣計(jì)算N個(gè)集合的方法較為復(fù)雜， 如何更加簡(jiǎn)便地將M個(gè)集合轉(zhuǎn)換為N個(gè)集合進(jìn)行預(yù)測(cè)過(guò)程是值得商榷的； 其次， GETKF方法中計(jì)算擴(kuò)展后的集合誤差協(xié)方差矩陣采用有偏形式， 本文對(duì)此做出了修改，采用無(wú)偏形式的矩陣進(jìn)行計(jì)算； 此外， GETKF方法中使用的GC局地化函數(shù)是一個(gè)基于距離相關(guān)的局地裁剪函數(shù)， 該函數(shù)是一元局地化函數(shù)， 所以擴(kuò)展到多元模型存在限制； 最后， 由于GETKF方法將參數(shù)L固定為10， 所以在系統(tǒng)狀態(tài)變量個(gè)數(shù)較大的情況下， 擴(kuò)展后的預(yù)測(cè)集合擾動(dòng)矩陣不能很好地表示誤差變化?？梢钥闯鯣ETKF方法確實(shí)存在的不足之處， 對(duì)于GETKF方法的具體改進(jìn)將在下一節(jié)進(jìn)行說(shuō)明。

1.2 增益集合轉(zhuǎn)換卡爾曼濾波的改進(jìn)

基于GETKF方法， 本節(jié)從局地化函數(shù)平方根矩陣的選取、預(yù)測(cè)誤差協(xié)方差的無(wú)偏估計(jì)、隨機(jī)子采樣方法和GC局地化函數(shù)的限制4個(gè)方面對(duì)原方法進(jìn)行改進(jìn)， 以提高同化方法的效果。

1.2.1 局地化函數(shù)平方根的選取

對(duì)于公式(3)中局地化函數(shù)平方根矩陣列數(shù)L的選取， GETKF方法按照特征值降序排列， 取前十個(gè)特征值和對(duì)應(yīng)的特征向量構(gòu)成平方根矩陣W，即固定地取L為10。很顯然， 在應(yīng)對(duì)狀態(tài)變量個(gè)數(shù)n很大時(shí)， 前10個(gè)特征值和特征向量不能完全表示系統(tǒng)變量的誤差變化， 所以本文重新定義了L的選取：

其中，E10%表示前10%的特征值的數(shù)量。即給定L的范圍， 通過(guò)進(jìn)行多次數(shù)值實(shí)驗(yàn)選取最優(yōu)的數(shù)值L。

1.2.2 預(yù)測(cè)誤差協(xié)方差的無(wú)偏估計(jì)

通常， 我們假定卡爾曼濾波中的誤差協(xié)方差矩陣是無(wú)偏的。在統(tǒng)計(jì)學(xué)上， 通過(guò)對(duì)集合成員與集合均值的偏差除以集合自由度(樣本個(gè)數(shù)減一)來(lái)實(shí)現(xiàn)。但是GETKF方法中的公式(5)～(7)采用的是有偏形式的誤差協(xié)方差矩陣， 這導(dǎo)致誤差的產(chǎn)生， 造成對(duì)誤差協(xié)方差的低估， 所以本文通過(guò)將公式(5)～(7)中的來(lái)修正協(xié)方差公式， 構(gòu)造無(wú)偏形式的誤差協(xié)方差矩陣， 從而在進(jìn)行分析時(shí)減少對(duì)同化結(jié)果的影響。

1.2.3 隨機(jī)子采樣方法

由于擴(kuò)展后的分析集合數(shù)量M大于原始預(yù)測(cè)集合數(shù)量N， 導(dǎo)致集合數(shù)量不能在預(yù)測(cè)模型中進(jìn)行傳遞。GETKF方法中利用一個(gè)膨脹因子[式(12)]來(lái)克服此問(wèn)題， 該方法需要計(jì)算擴(kuò)展后的分析誤差協(xié)方差矩陣以及原始預(yù)測(cè)誤差協(xié)方差矩陣， 如果變量的個(gè)數(shù)n很大， 計(jì)算難度大大增加。本文使用隨機(jī)子采樣方法進(jìn)行N列分析集合的選取， 計(jì)算簡(jiǎn)便。隨機(jī)子采樣方法如下：

1.2.4 GC局地化函數(shù)的限制

GETKF方法中使用GC函數(shù)作為局地化函數(shù)，但是GC函數(shù)僅僅是一個(gè)一元變量的局地化函數(shù)， 多元模型中不同變量之間的相關(guān)性應(yīng)該是有區(qū)別的，GC函數(shù)并沒(méi)有體現(xiàn)這一性質(zhì)， 所以在多元模型中使用GETKF方法存在限制。本文選用一個(gè)基于距離的多元函數(shù)Askey函數(shù)[21]作為局地化函數(shù)， 當(dāng)前此函數(shù)已經(jīng)用于EnKF方法， 并且在多元模型中進(jìn)行了實(shí)驗(yàn)驗(yàn)證， 得到了較好的同化效果。利用Askey函數(shù)將本文修改的方法擴(kuò)展到多元模型變量中， 探究方法在多元情況下的適用性。

2 實(shí)驗(yàn)方案

2.1 模型

2.1.1 KS模型

KS模型的表達(dá)式為，

KS模型是一個(gè)一元無(wú)量綱的非線性偏微分方程模型， 方程中二階項(xiàng)和四階項(xiàng)會(huì)增加模型的復(fù)雜性并導(dǎo)致混沌行為。在本文中， KS模型的偏微分方程通過(guò)一個(gè)指數(shù)時(shí)間的四階龍格庫(kù)塔數(shù)值格式進(jìn)行離散， 時(shí)間步為Δt= 0.25， 終止時(shí)間T為250。

2.1.2 淺水模型

淺水模型作為一個(gè)多元模型， 在當(dāng)前的數(shù)據(jù)同化研究中廣泛使用。其忽略摩擦效應(yīng)、科氏力以及非線性項(xiàng)的方程表示為：

其中，u和v分別表示水平和垂直速度， 水面高度由h表示，g為取值為9.8 m/s2的重力加速度。模型區(qū)域選擇矩形區(qū)域取2 200 km。淺水模型的偏微分方程組使用Lax-Wendroff有限差分方法計(jì)算， 時(shí)間步為Δt= 0.01， 終止時(shí)間T為8。

2.2 局地化函數(shù)

對(duì)于KS方程等一元模型， 我們可以使用GC函數(shù)作為CL方法中的局地化函數(shù)ρ， 其表達(dá)式如下：

式中，z表示網(wǎng)格點(diǎn)之間或網(wǎng)格點(diǎn)與觀測(cè)點(diǎn)之間的距離。由于本文使用模型空間CL方法， 所以z表示物理空間中網(wǎng)格點(diǎn)之間的距離。局地化長(zhǎng)度尺度c定義為為局地化半徑， 是一個(gè)可選參數(shù)。注意， 因子可以調(diào)整局地化函數(shù)達(dá)到最優(yōu)[22]。

相反， 多元模型中使用Askey函數(shù)的一個(gè)多元擴(kuò)展[21]作為局地化函數(shù)ρ， 形式為，

其中，m表示不同的狀態(tài)變量的總數(shù)。例如本文使用的淺水方程中m取3， 則Askey函數(shù)表示為：

即ρ的每個(gè)元素ijρ對(duì)應(yīng)于第i個(gè)變量和第j個(gè)變量的局地化函數(shù)矩陣。由于多元模型引入了一些不同變量之間的聯(lián)系， 增加了數(shù)據(jù)同化的復(fù)雜度與計(jì)算量， 因此多元局地化函數(shù)的構(gòu)造需要考慮這些不同變量之間的聯(lián)系以及聯(lián)系的強(qiáng)弱。uij表示第i個(gè)和第j個(gè)變量之間的局地化影響參數(shù)，

Askey函數(shù)中的參數(shù)z和c與GC函數(shù)中定義相同， 分別表示網(wǎng)格點(diǎn)之間的距離和局地化長(zhǎng)度尺度。B是一個(gè)Beta函數(shù)，s表示狀態(tài)變量所在的歐式空間的維度， 而且要滿足條件v≥(s+ 1)/2。

2.3 實(shí)驗(yàn)設(shè)計(jì)

其中， 矩陣H是觀測(cè)算子， 將真實(shí)狀態(tài)變量映射到觀測(cè)空間。假定觀測(cè)誤差之間不相關(guān)， 則矩陣R為對(duì)角觀測(cè)誤差協(xié)方差矩陣， 且對(duì)角元素取值為1。初始集合通過(guò)對(duì)真實(shí)值t x添加N次隨機(jī)誤差構(gòu)造。詳細(xì)的實(shí)驗(yàn)參數(shù)設(shè)置如下。

KS模型具有一維周期性， 設(shè)定狀態(tài)向量u的維度n為256， 初始真實(shí)值可以通過(guò)定義在周期域0≤x≤32π上的函數(shù)給出。實(shí)驗(yàn)中分別取集合數(shù)N為5或N為10， 觀測(cè)數(shù)p是256或p是235。數(shù)值實(shí)驗(yàn)的觀測(cè)頻率(同化間隔)不同， 假定每5或10個(gè)時(shí)間步長(zhǎng)存在可用的觀測(cè)值。

對(duì)于多元淺水模型， 假定模型定義在矩形網(wǎng)格區(qū)域， 每個(gè)網(wǎng)格點(diǎn)定義3個(gè)變量： 水面高度h， 水平速度u和垂直速度v。假定水平速度u和垂直速度v的初始真實(shí)狀態(tài)為零， 即u=v= 0， 水面高度h的初始真實(shí)狀態(tài)由如下方式構(gòu)造，

類似于KS模型中的參數(shù)設(shè)置， 在淺水方程中，集合數(shù)分別取N是5或N是10， 觀測(cè)數(shù)p為300或p為240， 每5或10個(gè)時(shí)間步長(zhǎng)存在可用的觀測(cè)值。表1中給出了上述具體實(shí)驗(yàn)參數(shù)的配置， 分別在兩

表1 不同實(shí)驗(yàn)條件的參數(shù)Tab. 1 Parameters for different experimental conditions

個(gè)模型中進(jìn)行8組不同的實(shí)驗(yàn)， 其中對(duì)于KS方程全觀測(cè)p為256， 部分觀測(cè)p為235； 對(duì)于淺水模型全觀測(cè)p為300， 部分觀測(cè)p是240。

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

首先， 分析GETKF方法和本文修改的方法(簡(jiǎn)稱GCL方法)在一元KS模型中的實(shí)驗(yàn)結(jié)果。為了研究?jī)煞N方法對(duì)預(yù)測(cè)誤差協(xié)方差矩陣的影響， 以實(shí)驗(yàn)1為例， 取前兩個(gè)同化時(shí)刻(存在觀測(cè)數(shù)據(jù)的時(shí)刻)，比較同化過(guò)程中兩種方法對(duì)于預(yù)測(cè)誤差協(xié)方差矩陣的影響。在圖1至圖4的每一幅圖中， 從左到右，從上到下， 依次展示了GC局地化函數(shù)、預(yù)測(cè)誤差協(xié)方差矩陣、局地化函數(shù)和預(yù)測(cè)誤差協(xié)方差矩陣做舒爾積運(yùn)算得到的矩陣、使用局地化函數(shù)平方根矩陣與預(yù)測(cè)集合擾動(dòng)矩陣做舒爾積運(yùn)算得到的預(yù)測(cè)誤差協(xié)方差矩陣。

圖1 實(shí)驗(yàn)1中， 第一個(gè)同化時(shí)刻GCL對(duì)預(yù)測(cè)誤差協(xié)方差矩陣的影響Fig. 1 Effect of GCL on the forecast error covariance matrix at the first assimilation time in experiment 1

圖3 實(shí)驗(yàn)1中， 第一個(gè)同化時(shí)刻GETKF對(duì)預(yù)測(cè)誤差協(xié)方差矩陣的影響Fig. 3 Effect of GETKF on the forecast error covariance matrix at the first assimilation time in experiment 1

圖4 實(shí)驗(yàn)1中， 第二個(gè)同化時(shí)刻GETKF對(duì)預(yù)測(cè)誤差協(xié)方差矩陣的影響Fig. 4 Effect of GETKF on the forecast error covariance matrix at the second assimilation time in experiment 1

由于KS模型具有一維周期性， 所以GC局地化函數(shù)ρ是一個(gè)對(duì)稱矩陣， 且僅在主對(duì)角線附近和矩陣邊緣存在較強(qiáng)的相關(guān)性。很顯然， 對(duì)于GETKF和GCL這兩種方法， 圖1至圖4中的第4幅子圖均顯示了使用局地化函數(shù)之后協(xié)方差矩陣遠(yuǎn)距離的偽相關(guān)得到消除， 并且和第3幅子圖相比效果類似， 說(shuō)明兩種方法對(duì)于虛假相關(guān)性的消除是有效果的。但是，比較圖2和圖4， 可以看出第2個(gè)時(shí)刻兩種方法的局地化效果是存在偏差的(圖例中數(shù)值區(qū)間不同)， 這可能是源于兩種方法的具體實(shí)施不同。但是總體上兩者對(duì)于處理偽相關(guān)的效果較好， 即截?cái)嗑嚯x以外的相關(guān)性被消除， 鄰近狀態(tài)變量和邊界上的相關(guān)性得到了保留。

圖2 實(shí)驗(yàn)1中， 第二個(gè)同化時(shí)刻GCL對(duì)預(yù)測(cè)誤差協(xié)方差矩陣的影響Fig. 2 Effect of GCL on the forecast error covariance matrix at the second assimilation time in experiment 1

以實(shí)驗(yàn)3和實(shí)驗(yàn)7為例， 分別對(duì)GETKF方法和GCL兩種方法的同化效果進(jìn)行評(píng)估。如圖5和圖6中不同方法的曲線軌跡所示， 大體上可以看出相比于GETKF方法的同化效果， 本文修改的GCL方法更加接近于真實(shí)值的軌跡曲線， 并且減少了濾波發(fā)散的發(fā)生， 展示了良好的同化效果， 說(shuō)明在一元數(shù)值預(yù)測(cè)模型下， GCL方法的局地化效果優(yōu)于GETKF方法。

圖5 實(shí)驗(yàn)3， KS模型第一個(gè)變量的同化效果Fig. 5 Comparison of the assimilation effect of the first variable of the KS model in experiment 3

圖6 實(shí)驗(yàn)7中， KS模型第一個(gè)變量的同化效果Fig. 6 Comparison of the assimilation effect of the first variable of the KS model in experiment 7

均方根誤差(root mean square error， RMSE)可以量化不同局地化方法的效果， 直觀表示不同方法的優(yōu)劣程度。模型狀態(tài)變量的整體均方根誤差表示為，

其中，l表示同化時(shí)間的長(zhǎng)度，表示第i個(gè)時(shí)刻第j個(gè)變量的分析值，表示第i個(gè)時(shí)刻第j個(gè)變量的真實(shí)值。此外， 由預(yù)測(cè)集合通過(guò)統(tǒng)計(jì)意義計(jì)算得到預(yù)測(cè)誤差協(xié)方差矩陣的秩至多為N-1， 該矩陣是一個(gè)低秩的矩陣， 這在傳遞誤差時(shí)易造成矩陣的病態(tài)， 因此需要提升矩陣的秩。理論上CL方法對(duì)于矩陣秩的提升是有效果的， 但是由于CL方法應(yīng)用于ETKF方法的局限性， 使得效果大打折扣， 所以除了計(jì)算同化過(guò)程中的RMSE外， 也會(huì)對(duì)比GETKF和GCL方法使用前后預(yù)測(cè)誤差協(xié)方差矩陣秩的變化。具體如表2所示， 首先， 兩種方法的預(yù)測(cè)誤差協(xié)方差矩陣的秩與不使用局地化方法相比， 顯然均有了提高。GCL方法對(duì)于矩陣秩的提高更為明顯。其次， GCL方法的RMSE明顯小于GETKF方法的RMSE(實(shí)驗(yàn)6有偏差)， 說(shuō)明在一元模型中， GCL方法的同化效果要優(yōu)于GETKF方法的同化效果。

表2 KS模型中不同實(shí)驗(yàn)條件下RMSE和矩陣秩的對(duì)比Tab. 2 Comparison of the RMSE and matrix rank under different experimental conditions in the KS model

續(xù)表

表3 中展示了主流的數(shù)據(jù)同化方法(EnKF)的RMSE， 采用實(shí)驗(yàn)1中的觀測(cè)數(shù)和同化間隔， 同時(shí)使用變化的集合數(shù)進(jìn)行數(shù)值實(shí)驗(yàn)。EnKF方法在集合數(shù)N是5的情況下， RMSE的數(shù)值遠(yuǎn)遠(yuǎn)大于GCL方法(實(shí)驗(yàn)1)， 隨著集合數(shù)N增大到100， RMSE的數(shù)值在降低， 與GCL方法(實(shí)驗(yàn)1)的RMSE相當(dāng)， 但是程序的運(yùn)行時(shí)間也在不斷增加。當(dāng)集合數(shù)為N取200時(shí)，EnKF方法的RMSE數(shù)值低于GCL方法(實(shí)驗(yàn)1)， 表現(xiàn)出了良好的數(shù)據(jù)同化效果， 但是程序運(yùn)行時(shí)間是GCL方法使用集合數(shù)N為5時(shí)的數(shù)倍(實(shí)驗(yàn)1條件下，GCL方法運(yùn)行時(shí)間為4.83 s)。EnKF通過(guò)使用大集合數(shù)， 犧牲程序計(jì)算時(shí)間來(lái)獲取低RMSE， 并不是一種很好的處理方式， 并且在真實(shí)的數(shù)據(jù)同化方法應(yīng)用中， 選取超過(guò)100個(gè)集合是不合適的。上述實(shí)驗(yàn)結(jié)果表明兩種數(shù)據(jù)同化方法各有優(yōu)劣， 但是可以看出在集合數(shù)較小的情況下， 選擇GCL方法進(jìn)行數(shù)據(jù)同化的效果較好。

表3 KS模型中EnKF方法的RMSETab. 3 RMSE of the EnKF method in the KS model

最后， 我們將本文提出的GCL方法中的GC函數(shù)替換為Askey函數(shù)， 擴(kuò)展方法的適用范圍， 在多元淺水模型中檢驗(yàn)該方法。表4展示了在不同實(shí)驗(yàn)條件下， GCL方法與未使用局地化方法的RMSE以及預(yù)測(cè)誤差協(xié)方差矩陣的秩的對(duì)比。類似于一元數(shù)值模型中的實(shí)驗(yàn)結(jié)果， 使用GCL方法后的預(yù)測(cè)誤差協(xié)方差矩陣的秩與未使用局地化方法相比， 有了明顯的提升， 保證了誤差矩陣傳遞的可靠性， 從而GCL方法的RMSE與未使用局地化方法相比也有了明顯的降低， 說(shuō)明利用Askey函數(shù)將GCL方法擴(kuò)展到多元模型是成功的， 克服了GETKF方法只能用于一元模型的局限性。

表4 淺水模型中不同實(shí)驗(yàn)條件下RMSE和矩陣秩的對(duì)比Tab. 4 Comparison of the RMSE and matrix rank under different experimental conditions in the shallow water model

4 結(jié)論

為解決基于集合的卡爾曼濾波中集合數(shù)目過(guò)少導(dǎo)致的欠采樣問(wèn)題， 以及CL方法對(duì)于ETKF等平方根卡爾曼濾波方法的不適用， 本文對(duì)GETKF方法進(jìn)行了改進(jìn)， 提出GCL方法， 并與GETKF方法在相同實(shí)驗(yàn)條件下進(jìn)行了實(shí)驗(yàn)對(duì)比。結(jié)果表明， GCL方法改善了在一元模型情況中的同化效果， 明顯提高了預(yù)測(cè)誤差協(xié)方差矩陣的秩， 降低了分析值的均方根誤差。此外， 利用一個(gè)Askey函數(shù)實(shí)現(xiàn)局地化方法在多元模型中的擴(kuò)展， 使得提出的GCL方法可以用于多元模型中， 克服當(dāng)前多元模型中局地化方法的欠缺?？傮w來(lái)說(shuō)， 雖然GCL方法在處理欠采樣問(wèn)題以及對(duì)于多元模型局地化具有一定的效果， 但是GCL方法應(yīng)用到多元模型時(shí)， 為設(shè)計(jì)多元變量之間的相互作用關(guān)系及相關(guān)性， Askey函數(shù)需要預(yù)先指定的參數(shù)較多。在下一步的研究中應(yīng)找到合適的參數(shù)確定依據(jù)，使得提出的新方法具有更廣泛的適用性。