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

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

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

在進(jìn)行數(shù)據(jù)同化時， 為了克服欠采樣帶來的負(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方法通過一個膨脹因子校正預(yù)測誤差協(xié)方差的低估問題， 主要思想是通過膨脹集合均值和集合成員之間的誤差來增大預(yù)測誤差協(xié)方差。CI方法用于分析過程之前， 具體表示為：其中“←”表示對于之前狀態(tài)值的替代，r是所謂的膨脹因子。盡管CI方法克服了預(yù)測誤差協(xié)方差的低估問題， 但是協(xié)方差膨脹因子r解決不了虛假相關(guān)問題。CL方法在預(yù)測誤差協(xié)方差矩陣和局地化函數(shù)之間實施舒爾積運(yùn)算[7]， 通過截斷預(yù)測誤差協(xié)方差矩陣中預(yù)先指定距離之外的誤差相關(guān)性以解決虛假相關(guān)問題。此外， 舒爾積運(yùn)算能夠提高預(yù)測誤差協(xié)方差矩陣的秩[8]。LA方法以待更新狀態(tài)變量為中心， 建立一個虛擬的局部窗口， 得到該狀態(tài)變量預(yù)測誤差協(xié)方差的局部近似值， 在分析更新過程中， 將局部窗口外的集合擾動設(shè)定為零。例如， 通?？梢砸阅硞€狀態(tài)變量為中心， 僅僅同化與其固定距離為M之內(nèi)的觀測數(shù)據(jù)， 固定距離的長度決定了同化過程中觀測值的數(shù)量。顯然， 分別對狀態(tài)變量進(jìn)行局部同化的特點使LA方法可以利用并行方法進(jìn)行快速計算。

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

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

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

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

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

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

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

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

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

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

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

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

Z，

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

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

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

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

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

其中，

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

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

其中，

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

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

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

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

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

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

1.2.2 預(yù)測誤差協(xié)方差的無偏估計

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

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

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

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

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

2 實驗方案

2.1 模型

2.1.1 KS模型

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

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

2.1.2 淺水模型

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

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

2.2 局地化函數(shù)

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

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

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

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

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

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

2.3 實驗設(shè)計

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 實驗3， KS模型第一個變量的同化效果Fig. 5 Comparison of the assimilation effect of the first variable of the KS model in experiment 3

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

表2 KS模型中不同實驗條件下RMSE和矩陣秩的對比Tab. 2 Comparison of the RMSE and matrix rank under different experimental conditions in the KS model

續(xù)表

表3 中展示了主流的數(shù)據(jù)同化方法(EnKF)的RMSE， 采用實驗1中的觀測數(shù)和同化間隔， 同時使用變化的集合數(shù)進(jìn)行數(shù)值實驗。EnKF方法在集合數(shù)N是5的情況下， RMSE的數(shù)值遠(yuǎn)遠(yuǎn)大于GCL方法(實驗1)， 隨著集合數(shù)N增大到100， RMSE的數(shù)值在降低， 與GCL方法(實驗1)的RMSE相當(dāng)， 但是程序的運(yùn)行時間也在不斷增加。當(dāng)集合數(shù)為N取200時，EnKF方法的RMSE數(shù)值低于GCL方法(實驗1)， 表現(xiàn)出了良好的數(shù)據(jù)同化效果， 但是程序運(yùn)行時間是GCL方法使用集合數(shù)N為5時的數(shù)倍(實驗1條件下，GCL方法運(yùn)行時間為4.83 s)。EnKF通過使用大集合數(shù)， 犧牲程序計算時間來獲取低RMSE， 并不是一種很好的處理方式， 并且在真實的數(shù)據(jù)同化方法應(yīng)用中， 選取超過100個集合是不合適的。上述實驗結(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ò)展方法的適用范圍， 在多元淺水模型中檢驗該方法。表4展示了在不同實驗條件下， GCL方法與未使用局地化方法的RMSE以及預(yù)測誤差協(xié)方差矩陣的秩的對比。類似于一元數(shù)值模型中的實驗結(jié)果， 使用GCL方法后的預(yù)測誤差協(xié)方差矩陣的秩與未使用局地化方法相比， 有了明顯的提升， 保證了誤差矩陣傳遞的可靠性， 從而GCL方法的RMSE與未使用局地化方法相比也有了明顯的降低， 說明利用Askey函數(shù)將GCL方法擴(kuò)展到多元模型是成功的， 克服了GETKF方法只能用于一元模型的局限性。

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

4 結(jié)論

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