基于地球系統(tǒng)模式的局地化粒子濾波器與集合卡爾曼濾波器同化實(shí)驗(yàn)

張鈺婷，沈浙奇,2,3*，伍艷玲,2,3

(1.自然資源部第二海洋研究所 衛(wèi)星海洋環(huán)境動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310012；2.河海大學(xué) 海洋學(xué)院 資料同化與預(yù)測(cè)研究所，江蘇 南京 210098；3.南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室（珠海），珠海 廣東 519080)

1 引言

作為研究大氣、海洋科學(xué)熱門課題之一，資料同化技術(shù)不僅可以為海洋數(shù)值預(yù)報(bào)模式提供初始場(chǎng)，還可以構(gòu)造海洋再分析資料集，為海洋觀測(cè)計(jì)劃和數(shù)值預(yù)報(bào)模式物理量及參數(shù)等提供設(shè)計(jì)依據(jù)，近年來得到了廣泛的研究和應(yīng)用[1]。集合卡爾曼濾波器（EnKF）是一種有效的資料同化方法，自Evensen[2]于1994年首次提出以來經(jīng)過了20多年的發(fā)展和改進(jìn)，已經(jīng)在海洋資料同化中得到了廣泛的研究和應(yīng)用。近年來，隨著動(dòng)力模式的不斷發(fā)展和計(jì)算能力的提高，粒子濾波器（PF）由于不受模型線性和誤差高斯分布假設(shè)的約束，也逐漸成為了當(dāng)前資料同化算法研究的熱門[3]。

EnKF是集合預(yù)報(bào)同卡爾曼濾波器（KF）的結(jié)合，它使用集合表示模式變量的概率密度分布，并采用KF的更新公式提供線性模型高斯分布假設(shè)下的同化最優(yōu)解。由于EnKF避免了卡爾曼濾波器和擴(kuò)展卡爾曼濾波器（EKF）中的協(xié)方差更新模型，使得它能夠被應(yīng)用于大型地球物理模式。集合調(diào)整卡爾曼濾波器（EAKF）是在EnKF基礎(chǔ)上發(fā)展起來的衍生方法，它一般被視為一種確定性的EnKF格式[4]。EAKF從濾波理論出發(fā)推導(dǎo)了一個(gè)用于模式變量更新的算子，取代了傳統(tǒng)EnKF中的增益矩陣，且無需增加額外的觀測(cè)擾動(dòng)，在計(jì)算量上有一定的優(yōu)勢(shì)[4–5]。作為目前海洋資料同化所采用的主流方法之一—集合卡爾曼濾波器以及衍生方法隱含了預(yù)報(bào)集合高斯分布的假設(shè)，適用于線性系統(tǒng)的同化[6]。

相較于集合卡爾曼濾波，粒子濾波不含高斯假設(shè)，對(duì)非線性非高斯同化系統(tǒng)能產(chǎn)生更好的同化效果。粒子濾波算法基于貝葉斯估計(jì)理論，是貝葉斯公式的蒙特卡羅算法近似?！傲Ｗ印迸c集合卡爾曼濾波器中的集合成員相同，是用于表示模式變量的概率分布的集合樣本。如果計(jì)算資源充足，隨著粒子數(shù)目的增加，粒子的概率密度分布會(huì)逐漸趨向于真實(shí)狀態(tài)場(chǎng)的概率密度分布，粒子濾波器能夠?qū)崿F(xiàn)最優(yōu)貝葉斯估計(jì)的同化效果[3,7]。經(jīng)典的粒子濾波器使用似然函數(shù)計(jì)算每個(gè)粒子的標(biāo)量權(quán)重，因此當(dāng)狀態(tài)場(chǎng)空間的維數(shù)較大時(shí)，狀態(tài)場(chǎng)數(shù)值的微小變化會(huì)引起對(duì)應(yīng)權(quán)重在量級(jí)上的巨大變化，多數(shù)集合成員會(huì)因權(quán)重過小而失效，從而導(dǎo)致粒子濾波器的退化[8]。相當(dāng)多的方法已經(jīng)被提出來處理粒子濾波器中的退化[9]，如最優(yōu)重要性粒子濾波器[10]、等權(quán)重粒子濾波器[11]、集合卡爾曼粒子濾波器[12]、局地化粒子濾波器[8,13]等。本文主要考察的局地化粒子濾波器是最近才被提出的一種同化方法，它通過在經(jīng)典粒子濾波器中引入EnKF中常用的局地化方法來解決粒子濾波器的退化問題，得到了廣泛關(guān)注。

局地化方法自21世紀(jì)初在集合卡爾曼濾波器的同化中被提出以來[14]，已經(jīng)被廣泛地應(yīng)用于各種業(yè)務(wù)化集合同化系統(tǒng)中，得到了普遍認(rèn)可。由于模式集合成員的數(shù)量非常有限，在計(jì)算背景誤差協(xié)方差和增益矩陣的過程中往往會(huì)出現(xiàn)由于樣本不足而造成的虛假遠(yuǎn)距離相關(guān)，造成錯(cuò)誤的同化更新。通過引入局地化可以使得同化的更新過程在一個(gè)較小子空間中進(jìn)行，從而抑制這種虛假相關(guān)，并大大降低計(jì)算量。利用相似的思想，最近的一些工作開始在粒子濾波器中使用局地化，例如，Poterjoy[13]、Shen等[15]、Penny和Miyoshi[8]改善的局地化粒子濾波器已經(jīng)被初步驗(yàn)證能夠使用與集合卡爾曼濾波器相當(dāng)?shù)募铣蓡T數(shù)來避免粒子退化現(xiàn)象，也逐漸地被用于地球物理模式的資料同化中。

本文在耦合的通用地球系統(tǒng)模式（Community Earth System Modal,CESM）中開展了集合卡爾曼濾波器和局地化粒子濾波器的觀測(cè)系統(tǒng)模擬試驗(yàn)。通過同化模擬的衛(wèi)星海表溫度（SST）資料，考察不同局地化參數(shù)對(duì)于兩種濾波器方法的不同影響。在此基礎(chǔ)上，進(jìn)一步比較了兩種濾波器方法的同化效果，探討了兩類方法的優(yōu)缺點(diǎn)以及發(fā)展前景。

2 數(shù)據(jù)與方法

2.1 模式和同化系統(tǒng)

實(shí)驗(yàn)使用的模式為CESM，它是由美國國家大氣研究中心（NCAR）于2010年發(fā)布的新一代地球系統(tǒng)模式，是目前最先進(jìn)、使用最廣泛的地球系統(tǒng)模式之一。該模式采用模塊化框架，主體由大氣、海洋、陸地、海冰、河流等幾大模塊組成，并由耦合器（CPL7）管理模塊間的數(shù)據(jù)信息交換和模式運(yùn)行。實(shí)驗(yàn)使用了CESM1.2.1版本，它的海洋模式為POP2，大氣模式為CAM4。本實(shí)驗(yàn)使用全耦合的模式設(shè)置（B compset）和當(dāng)代（present day）的外強(qiáng)迫。采用的模式分辨率為0.9×1.25_gx1v6，即大氣模式水平分辨率為 0.9°×1.25°，垂向26層；海洋模式水平分辨率接近1°，在赤道區(qū)域緯向加密約為0.5°，垂向60層。

實(shí)驗(yàn)所涉及的兩種同化方法都是借助NCAR開發(fā)的資料同化 研究 平臺(tái)（Data Assimilation Research Testbed，DART）實(shí)現(xiàn)。DART是由NCAR的數(shù)據(jù)同化研究部門開發(fā)和維護(hù)的一個(gè)開源軟件。它提供了多種確定性和隨機(jī)性濾波器算法，包括集合調(diào)整卡爾曼濾波器、集合卡爾曼濾波器、核濾波器和粒子濾波器等[16]。我們建立了DART和POP2模式的接口，實(shí)現(xiàn)了耦合模式框架下的海洋觀測(cè)資料的弱耦合同化。鑒于計(jì)算量的考慮，我們使用了20個(gè)集合成員。在所有的實(shí)驗(yàn)中，我們也采用了固定系數(shù)為1.02的協(xié)方差膨脹和針對(duì)20個(gè)集合成員的樣本誤差訂正（Sampling Error Correction）方案[17]。

2.2 數(shù)據(jù)

本文采用的實(shí)驗(yàn)方案為觀測(cè)系統(tǒng)模擬試驗(yàn)（OSSE），其基本思想是使用一組真值實(shí)驗(yàn)作為參考值，然后根據(jù)實(shí)際存在的觀測(cè)系統(tǒng)對(duì)真值取樣，并且疊加上具有給定方差的隨機(jī)誤差來模擬觀測(cè)。在實(shí)驗(yàn)中同化上述模擬的觀測(cè)資料，并利用真值實(shí)驗(yàn)得到的參考值來評(píng)估不同同化方法得到的分析場(chǎng)。

圖1 觀測(cè)系統(tǒng)模擬試驗(yàn)流程圖設(shè)計(jì)Fig.1 Flow chart design of observation system simulation experiment

圖1顯示了OSSE的流程圖，其中我們利用真值實(shí)驗(yàn)生成每周一次的全球海表溫度觀測(cè)，其空間分辨率為1°×1°，并疊加上符合高斯分布的隨機(jī)誤差作為觀測(cè)誤差。為簡(jiǎn)化討論，我們針對(duì)所有觀測(cè)點(diǎn)使用相同的觀測(cè)誤差標(biāo)準(zhǔn)差。參考OISST[18]中的SST觀測(cè)誤差的全球平均標(biāo)準(zhǔn)差，假設(shè)所有位置的誤差標(biāo)準(zhǔn)差均為0.3℃。

此外，本文使用了哈德雷中心（Met Office Hadley Center）提供的長(zhǎng)度為 100 a（1919–2018）的月平均HadISST (Hadley Centre Sea Ice and Sea Surface Temperature data set)再分析資料[19]，主要用于計(jì)算 SST 資料對(duì)不同距離的變量的相關(guān)系數(shù)，并與模式20個(gè)集合成員計(jì)算的樣本相關(guān)系數(shù)進(jìn)行比較，為局地化系數(shù)的選取提供理論依據(jù)（見3.1.2節(jié)）。

2.3 同化方法和局地化

2.3.1 EAKF 中的局地化

集合調(diào)整卡爾曼濾波器（EAKF）[4]是在集合卡爾曼濾波器[2]的基礎(chǔ)上發(fā)展起來的。EAKF避免了集合卡爾曼濾波器中對(duì)觀測(cè)資料的擾動(dòng)，且在集合樣本數(shù)目較少時(shí)仍能得到較為滿意的結(jié)果。EAKF基于線性回歸理論逐個(gè)順次地同化觀測(cè)資料，使用觀測(cè)空間的局地化后可以將每個(gè)觀測(cè)點(diǎn)的更新范圍限制在一個(gè)較小的子空間中進(jìn)行，節(jié)省了計(jì)算量。

EAKF的同化步驟可以表示如下，若用x表示狀態(tài)向量，用yo表示誤差方差為的觀測(cè)標(biāo)量，h表示觀測(cè)算子，則EAKF公式的第一步是使用觀測(cè)算子h將狀態(tài)空間的每個(gè)集合成員投影到觀測(cè)點(diǎn)上，作為每個(gè)成員對(duì)應(yīng)觀測(cè)的先驗(yàn)估計(jì)，如下

式中，下標(biāo)p代表先驗(yàn)值；n代表集合成員

式中，下標(biāo)u代表后驗(yàn)值

接著使用以下公式計(jì)算觀測(cè)空間的每個(gè)成員的后驗(yàn)估計(jì)

EAKF最后利用相關(guān)系數(shù)將每個(gè)觀測(cè)的增量回歸到模式變量中

式中，xm,n表示模式第m個(gè)分量的第n個(gè)集合成員；σxm,y為xm與yp的協(xié)方差。

由于集合模式成員數(shù)有限，因此所模擬的背景誤差協(xié)方差會(huì)存在一定的虛假信息。這些虛假信息隨著距離的增加而增加，且會(huì)在同化過程中造成虛假的更新量，影響結(jié)果的準(zhǔn)確性。為了避免這些問題，我們?cè)贓AKF中引入了局地化方案，使用一個(gè)數(shù)值介于0和1之間且隨距離的大小單調(diào)遞減的因子 ρ，并將式（5）替換為

公式（6）中的 ρ被稱為局地化因子，是一個(gè)依賴于距離的函數(shù)，一般可以由以下公式[20]計(jì)算得到

式中，dn代表模式點(diǎn)xm,n和觀測(cè)點(diǎn)yp,n的距離；c是一個(gè)局地化參數(shù)，可以控制 ρ的去相關(guān)距離。由式（7）可知，當(dāng)dn大 于2c時(shí)， ρ =0，從而觀測(cè)點(diǎn)yp,n和模式點(diǎn)xm,n完全不 相關(guān)。

2.3.2 LPF 中的局地化

相比于EAKF，局地化的粒子濾波器（LPF）最近才被提出來解決傳統(tǒng)粒子濾波器中的退化問題[13]，因此關(guān)于它的研究結(jié)果相對(duì)較少。在作者所了解的文獻(xiàn)中，目前尚沒有將LPF應(yīng)用于CESM模式的研究也鮮有對(duì)LPF的局地化參數(shù)進(jìn)行的討論。

從原理上說，在粒子濾波器中采用局地化是為了使用較小的計(jì)算成本來避免濾波退化問題。粒子濾波器的原理是給每個(gè)集合成員配給一個(gè)權(quán)重，用集合成員的加權(quán)組合來表示同化變量的完整概率分布密度函數(shù)，并隨著同化的進(jìn)行不斷更新權(quán)重。經(jīng)典的粒子濾波器使用一個(gè)標(biāo)量的數(shù)來代表權(quán)重，因此當(dāng)模式的狀態(tài)變量維數(shù)巨大時(shí)，很容易造成其中一個(gè)成員占據(jù)大部分權(quán)重，而其他成員權(quán)重都幾乎為0的現(xiàn)象這就是濾波退化[2]。

在經(jīng)典粒子濾波器中，狀態(tài)場(chǎng)Xn的標(biāo)量權(quán)重可由以下公式計(jì)算

式中，wn表示第n個(gè)粒子的標(biāo)量權(quán)重；m是所有觀測(cè)點(diǎn)的總數(shù)；觀測(cè)算子hi用 于將模式預(yù)報(bào)場(chǎng)投影到其中的某個(gè)觀測(cè)點(diǎn)yi上；是 yi的 方差，這里的正比例符號(hào)“∝”意味著計(jì)算等式右邊之后還需要進(jìn)行一次標(biāo)準(zhǔn)化使得所有權(quán)重的和為1。當(dāng)m很大時(shí)，多個(gè)指數(shù)函數(shù)相乘導(dǎo)致不同粒子的權(quán)重很容易有量級(jí)上的差異。

在粒子濾波器中引入局地化的基礎(chǔ)是將粒子的權(quán)重?cái)U(kuò)展成為矢量，即不同的模式網(wǎng)格點(diǎn)使用不同的權(quán)重。在此基礎(chǔ)上應(yīng)用粒子濾波器可以將狀態(tài)分析過程轉(zhuǎn)移到一個(gè)較小子空間中進(jìn)行，由于子空間的差異性會(huì)增加集合成員的多樣性，從而降低集合退化的可能。粒子濾波器中的局地化主要用于計(jì)算矢量權(quán)重，使其適用于大型地球物理模式。但是由于粒子濾波器的同化原理與卡爾曼濾波器有本質(zhì)的不同，局地化的引入方式也有所不同[21]。局地化粒子濾波器的第一步是將式（8）中的標(biāo)量權(quán)重?cái)U(kuò)展到矢量權(quán)重，也就是說狀態(tài)場(chǎng)中各個(gè)不同的分量采用不同的局部權(quán)重，與此同時(shí)，集合成員的重取樣也都在局地進(jìn)行。Poterjoy[13]借助局地化公式，利用局地化因子 ρ給出不同分量的權(quán)重公式如下

2.4 同化實(shí)驗(yàn)設(shè)計(jì)

本文首先通過敏感性實(shí)驗(yàn)考察局地化參數(shù)對(duì)于EAKF和LPF的同化效果的影響。局地化的參數(shù)數(shù)值對(duì)應(yīng)了觀測(cè)點(diǎn)和模式網(wǎng)格點(diǎn)的去相關(guān)距離：假設(shè)局地化參數(shù)的值為c，那么觀測(cè)點(diǎn)與距離2c以上的模式網(wǎng)格點(diǎn)的相關(guān)系數(shù)為0。CESM模式中的距離使用弧度制單位表示，例如d=0.1 rad，那么它實(shí)際對(duì)應(yīng)的赤道經(jīng)度為為了處理海洋在垂直和水平方向的不同尺度，同化算法使用一個(gè)垂向歸一化尺度系數(shù)（Vertical Normalization Factor，以下記做ν，單位為m/rad）來進(jìn)行垂向距離的轉(zhuǎn)化。假設(shè)水平方向的距離為dh，單位為（°），垂直方向的距離為dv，單位為m，那么兩點(diǎn)之間的模式距離如下計(jì)算

根據(jù)d和c的比值，可以使用式（7）計(jì)算局地化因子 ρ 。顯然ν的數(shù)值越大，相同的實(shí)際垂向距離dv條件下式（10）右端的第二項(xiàng)越小，相同的c可以影響到越遠(yuǎn)的水平距離。而當(dāng)ν為無窮大的時(shí)候，實(shí)際上就關(guān)閉了垂直方向的局地化，局地化因子僅由dh決定。

實(shí)驗(yàn)所選取的同化數(shù)據(jù)為海表溫度，集合成員數(shù)為20，同化頻率為每7 d同化一次，同化實(shí)驗(yàn)的時(shí)間為12個(gè)月。同化的初始集合由如下方法產(chǎn)生：我們首先對(duì)CESM模式進(jìn)行了100 a的自由積分，然后在得到的初始場(chǎng)上層30層溫度變量上疊加偽隨機(jī)場(chǎng)[22]的擾動(dòng)，構(gòu)成20個(gè)集合成員。我們對(duì)該集合進(jìn)行2 a的模式積分，使得每個(gè)集合成員各變量之間保持一定的動(dòng)力平衡，這些積分的結(jié)果則被用來作為同化實(shí)驗(yàn)的初始集合。

設(shè)計(jì)的敏感性實(shí)驗(yàn)如表1所示：我們首先考察EAKF中的局地化方案，比較了不同局地化參數(shù)c和垂向歸一化尺度系數(shù)ν對(duì)于同化效果的影響。我們先固定c為 0.1 rad，然后分別設(shè)定ν=1 000 m/rad,1 500 m/rad 2 000 m/rad，以及+∞（即垂直方向不采用局地化）來考察垂向局地化方案對(duì)于同化效果的影響。并且分別將實(shí)驗(yàn)名稱記作Kc0.1v1000、Kc0.1v1500、Kc0.1v2000以及Kc0.1vinf。根據(jù)CESM氣候模式的距離設(shè)置Kc0.1v1000試驗(yàn)中，觀測(cè)點(diǎn)在海表與水平距離超過大約11.5°的模式格點(diǎn)相關(guān)系數(shù)為0，同時(shí)與其正下方垂直距離超過200 m的模式格點(diǎn)相關(guān)系數(shù)也為0。然后，我們固定最優(yōu)垂向歸一化尺度系數(shù)ν，來考察水平局地化方案對(duì)同化效果的影響。由于前一組實(shí)驗(yàn)得出的結(jié)論為EAKF的最優(yōu)垂向局地化方案為垂向關(guān)閉局地化（見3.1.1節(jié)），在此基礎(chǔ)上，分別設(shè)置局地化參數(shù)c為 0.05 rad、0.1 rad、0.2 rad、0.3 rad，分別將對(duì)應(yīng)的實(shí)驗(yàn)名稱記為Kc0.05vinf、Kc0.1vinf、Kc0.2vinf Kc0.3vinf。類似地，我們也針對(duì)LPF探討了系數(shù)ν和參數(shù)c的相關(guān)問題，利用敏感性實(shí)驗(yàn)尋找最優(yōu)的局地化方案，揭示LPF和EAKF對(duì)于局地化的不同要求最后，基于使用最優(yōu)的局地化方案，我們比較了LPF和EAKF的同化效果，揭示LPF的潛在優(yōu)勢(shì)與不足同時(shí)，為了顯示同化效果，本實(shí)驗(yàn)設(shè)置了不進(jìn)行任何同化的控制試驗(yàn)(下文簡(jiǎn)稱FREERUN)進(jìn)行對(duì)照。

表1 實(shí)驗(yàn)列表Table 1 Experimental list

3 實(shí)驗(yàn)結(jié)果與討論

3.1 EAKF 局地化對(duì)同化效果影響

3.1.1 EAKF 的垂向局地化方案

我們首先固定局地化參數(shù)c為0.1 rad，考察不同垂向歸一化尺度系數(shù)ν對(duì)于同化效果的影響。實(shí)驗(yàn)采用的系數(shù)ν為 1 000 m/rad、1 500 m/rad、2 000 m/rad與局地化參數(shù)c相乘以后，對(duì)應(yīng)的距離分別為100 m 150 m/200 m，也就是說，局地化系數(shù)ρ在超過 200 m 300 m/400 m 的數(shù)值為 0，因此，SST 觀測(cè)最多能影響的垂直深度為 200 m/300 m/400 m。實(shí)驗(yàn)選取了均方根誤差（RMSE）作為評(píng)判同化效果優(yōu)劣的標(biāo)準(zhǔn)，其計(jì)算公式為

我們首先根據(jù)RMSE的垂直分布討論系數(shù)ν對(duì)同化效果的影響。圖2為設(shè)置不同ν時(shí)，EAKF同化SST資料后的垂向溫鹽RMSE。由于同化資料為SST，因此溫度的同化效果要明顯強(qiáng)于鹽度的，在深度較淺時(shí)其優(yōu)勢(shì)更加顯著。對(duì)于溫度變量，使用不同的系數(shù)ν同化SST得到的RMSE在100 m以淺并沒有顯著差異，且都遠(yuǎn)小于控制實(shí)驗(yàn)的RMSE。而隨著深度的增加，變量與觀測(cè)之間的距離增大，一方面所有實(shí)驗(yàn)的同化效果逐漸削弱，另一方面不同實(shí)驗(yàn)的RMSE出現(xiàn)差異。特別地，當(dāng)深度約為150～300 m時(shí)，在關(guān)閉垂向局地化的Kc0.1vinf實(shí)驗(yàn)產(chǎn)生的RMSE是最小的。根據(jù)式（6）可知，EAKF利用SST更新深層溫度變量的原理是使用SST和給定層溫度之間的相關(guān)性將表層的觀測(cè)增量回歸到深層。因?yàn)樯顚訙囟鹊碾x散度（標(biāo)準(zhǔn)差）很小，僅有O（10?2），所以即使只使用20個(gè)集合成員計(jì)算相關(guān)系數(shù)，表層和深層的相關(guān)系數(shù)也能夠正確表達(dá)。所以，雖然Kc0.1vinf實(shí)驗(yàn)沒有引入垂向局地化，較深層的溫度也能夠得到正確更新。

圖2 不同垂向局地化方案 EAKF 實(shí)驗(yàn)中區(qū)域平均（60°S～60°N，環(huán)地球）垂向均方根誤差Fig.2 Regional mean (60°S?60°N,ring the earth) root mean square error in EAKF experiments with different vertical localization schemes

另一方面，同化SST對(duì)于鹽度的改進(jìn)基于溫鹽相關(guān)得到。在150 m以淺，仍然可以發(fā)現(xiàn)鹽度變量的RMSE小于控制實(shí)驗(yàn)的RMSE，但是不同垂直局地化方案對(duì)于鹽度同化效果的差異主要在較淺層而非深層出現(xiàn)。同時(shí)，圖2b的結(jié)果表明同化SST時(shí)ν=1 500 m/rad/2 000 m/rad的局地化能更好地改進(jìn)淺層的鹽度，由式（10）可知，在局地化參數(shù)c固定的條件下，對(duì)于淺層的變量，尺度系數(shù)ν越大，dh起到的作用也相對(duì)越大。體現(xiàn)在鹽度上就是相同深度的鹽度變量獲取海表觀測(cè)的信息就越多。所以我們可以斷言，淺層的鹽度同化差異實(shí)際上受到水平局地化的影響更大綜合兩者的討論，在垂直方向不采用局地化是相對(duì)更優(yōu)的方案。

3.1.2 EAKF 局地化參數(shù)

對(duì)EAKF垂直局地化方案的敏感實(shí)驗(yàn)結(jié)果顯示，EAKF的同化效果在關(guān)閉垂直方向的局地化時(shí)相對(duì)較好。因此，我們?cè)陉P(guān)閉垂直局地化的基礎(chǔ)上，分析不同參數(shù)c對(duì)應(yīng)的水平局地化方案對(duì)EAKF同化效果的影響，第2組實(shí)驗(yàn)采用的局地化參數(shù)分別為0.05 rad、0.1 rad、0.2 rad、0.3 rad。圖 3 從時(shí)間尺度上分析各物理量RMSE的變化情況。首先分析表層溫度、鹽度（SSS）和海表高度（SSH）變量的RMSE變化情況，如圖3a所示，SST從啟用同化的第2個(gè)月開始就一直保持較小的RMSE，且使用不同參數(shù)實(shí)驗(yàn)的結(jié)果沒有明顯差異。這是由于SST的觀測(cè)資料比較密集，即使在某個(gè)點(diǎn)的同化中引入虛假更新，也很容易被其他點(diǎn)的正確更新所抵消。而對(duì)其他兩個(gè)變量來說，同化SST并不能立刻減小誤差，需要進(jìn)行一段時(shí)間的持續(xù)同化和模式積分，誤差才能顯著減小，這點(diǎn)對(duì)于SSS特別明顯。圖3b顯示，SSS的RMSE尚處于下降階段，沒有達(dá)到穩(wěn)定，其同化效果與局地化參數(shù)的大小并沒有明顯的相關(guān)關(guān)系，在最后3個(gè)月，幾種參數(shù)實(shí)驗(yàn)的結(jié)果沒有顯著差別。而SSH在局地化參數(shù)稍大時(shí)同化效果較好，這說明SST和SSH之間的相關(guān)性較好，不容易產(chǎn)生虛假相關(guān)。

圖3 不同局地化參數(shù) EAKF 實(shí)驗(yàn)中區(qū)域平均（60°S～60°N，環(huán)地球）的均方根誤差時(shí)間序列Fig.3 RMSE time series of regional mean (60°S?60°N,ring the earth) in EAKF experiments with different local parameters

我們進(jìn)一步分析較深層的同化效果發(fā)現(xiàn)，當(dāng)深度為200 m時(shí)，圖3d和圖3e表現(xiàn)出最優(yōu)同化效果的局地化參數(shù)c=0.1 rad。這是因?yàn)槲覀冴P(guān)閉了垂向局地化功能，因而使得一定水平距離之外的海表觀測(cè)資料能夠無差別地影響整個(gè)水柱，隨著深度的增加，實(shí)際的距離增大，遠(yuǎn)距離觀測(cè)帶來的虛假相關(guān)也會(huì)出現(xiàn)。從溫度和鹽度兩個(gè)變量，都可以看出c=0.1 rad的同化效果好于c=0.05 rad，這是因?yàn)楦嗟挠^測(cè)資料被用來更新200 m深度的變量；而它的同化效果也好于局地化參數(shù)更大的另外兩個(gè)實(shí)驗(yàn)，這是因?yàn)楹髢烧邥?huì)因?yàn)樘摷傧嚓P(guān)而帶來虛假的更新。

由于集合濾波器的集合成員數(shù)量有限，在相距較遠(yuǎn)的兩個(gè)點(diǎn)之間會(huì)產(chǎn)生虛假的相關(guān)，造成錯(cuò)誤的更新，這是在EAKF中引入局地化的主要原因。因此，如果可以量化表示給定集合成員數(shù)可能造成的虛假相關(guān)，就可以相應(yīng)地選擇最優(yōu)的局地化參數(shù)。針對(duì)某個(gè)給定點(diǎn)上的觀測(cè)，以1月的控制實(shí)驗(yàn)預(yù)報(bào)值為例，我們使用以下公式定義對(duì)應(yīng)局地化參數(shù)c的虛假相關(guān)占比（Spurious Correlation Ratio,SCR）

圖4為SST相關(guān)系數(shù)的空間分布，我們分別選取了位于太平洋的兩個(gè)點(diǎn)（0°，180°）、（20°S，120°W），位于印度洋的一個(gè)點(diǎn)（20°S，60°E），以及位于大西洋的一個(gè)點(diǎn)（20°N，40°W），計(jì)算了使用不同局地化參數(shù)下的有效相關(guān)系數(shù)，樣本相關(guān)系數(shù)和分析相關(guān)系數(shù)，分別對(duì)應(yīng)圖4中的前4列、第5列以及第6列，圖中的白色區(qū)域表示相關(guān)系數(shù)小于0部分。從圖中可以看出，幾乎所有位置的樣本相關(guān)系數(shù)都大于分析相關(guān)系數(shù)，因此在某些較遠(yuǎn)區(qū)域，樣本會(huì)夸大SST的相關(guān)性，進(jìn)而引入虛假相關(guān)。我們也隨機(jī)選取了其他多個(gè)格點(diǎn)進(jìn)行相同的計(jì)算分析，得到的結(jié)論類似。通過引入局地化算法，設(shè)置恰當(dāng)?shù)木值鼗霃娇梢栽诒ＷC近距離相關(guān)的前提下去除遠(yuǎn)距離虛假相關(guān)，避免遠(yuǎn)距觀測(cè)引起的錯(cuò)誤更新[23]。

圖4 海表溫度相關(guān)系數(shù)Fig.4 Correlation coefficient of sea surface temperature

我們進(jìn)一步計(jì)算了幾種不同局地化方案下的虛假相關(guān)占比（表2），發(fā)現(xiàn)局地化方案的引入大大減小了虛假相關(guān)的占比。其中當(dāng)局地化參數(shù)c=0.1 rad時(shí)，既不會(huì)因?yàn)榘霃竭^大而引入一些不必要的虛假相關(guān)，也不會(huì)因?yàn)榫值鼗霃竭^小而濾去大部分正確相關(guān)。因此c=0.1 rad相對(duì)而言是更優(yōu)的局地化參數(shù)，與前面敏感性實(shí)驗(yàn)的結(jié)果一致。

表2 SST相關(guān)系數(shù)虛假相關(guān)占比Table 2 Proportion of false correlation in SST correlation coefficient

3.2 LPF 局地化對(duì)同化效果影響

3.2.1 LPF 垂向局地化方案

第3組實(shí)驗(yàn)考慮了LPF垂向局地化方案對(duì)同化效果的影響，類似EAKF局地化實(shí)驗(yàn)，我們首先固定c為 0.1 rad，設(shè)置垂向歸一化系數(shù)ν=500/1000/1500/m/rad。首先，我們從海洋溫度和鹽度RMSE的垂向分布分析LPF垂直局地化對(duì)同化效果的影響。如圖5所示，在淺層，當(dāng)ν=1000/1500/∞m/rad，即對(duì)應(yīng)最大垂向同化距離為200 m/300 m時(shí)，LPF的同化效果較好。由式（10）可知，當(dāng)局地化參數(shù)c固定，垂直距離dv較小且固定時(shí)，觀測(cè)點(diǎn)到它能同化到的最遠(yuǎn)模式網(wǎng)格點(diǎn)的水平距離dh隨著v的增大而增大。在較淺的深度，v的值越大就能在水平方向取得越多的觀測(cè)資料來計(jì)算權(quán)重，如果不發(fā)生退化，就能獲得越好的同化效果。而當(dāng)深度較大時(shí)，無論對(duì)于溫度還是鹽度來說，都可以發(fā)現(xiàn)關(guān)閉垂向局地化的LPF的效果較差，甚至?xí)斐韶?fù)面同化效果—即同化后的誤差反而大于控制實(shí)驗(yàn)（例如大于500 m的溫度誤差）。這是由于粒子濾波器依賴權(quán)重的計(jì)算和重采樣來更新變量，關(guān)閉垂向局地化意味著所有層都會(huì)根據(jù)表層觀測(cè)計(jì)算的權(quán)重來重分配樣本，這會(huì)過度夸張表層觀測(cè)的影響范圍，從而給深層的變量帶來錯(cuò)誤的更新[24]。綜合溫度和鹽度的實(shí)驗(yàn)結(jié)果，可以發(fā)現(xiàn)當(dāng)ν=1 000 m/rad，也就是垂向局地化距離為100 m時(shí)，LPF的同化效果最好。

圖5 不同垂向局地化方案 LPF 實(shí)驗(yàn)中區(qū)域平均（60°S～60°N，環(huán)地球）垂向均方根誤差Fig.5 Regional mean (60°S?60°N,ring the earth) root mean square error in LPF experiments with different vertical localization schemes

3.2.2 LPF 局地化參數(shù)

接下來，我們固定LPF最優(yōu)垂向局地化距離為100 m，分析LPF的不同局地化參數(shù)c對(duì)同化效果的影響。由式（10）可知，考慮垂向局地化的情況下，為了保持垂直方向的局地化距離為100 m，針對(duì)不同的c需要采用不同的ν。圖6為采用不同局地化參數(shù)的LPF同化SST資料后的垂向RMSE，從圖中可以很明顯的看出，LPF的效果對(duì)局地化參數(shù)c非常敏感。當(dāng)c=0.2 rad時(shí)，同化效果非常不理想，這主要是由濾波退化造成的。c的數(shù)值越大，越多的觀測(cè)能同時(shí)影響給定點(diǎn)的權(quán)重，從而導(dǎo)致該點(diǎn)的權(quán)重更容易退化。集合退化會(huì)大大降低粒子濾波器的效率，特別地，我們可以發(fā)現(xiàn)c=0.2 rad時(shí)的鹽度沒有同化效果。另一方面，當(dāng)c的數(shù)值太小時(shí)，由于獲取資料過少，無法正確表達(dá)狀態(tài)變量的概率分布密度函數(shù)，也會(huì)導(dǎo)致同化效果不佳。圖7進(jìn)一步展示了不同局地化參數(shù)得到的分析場(chǎng)的表層變量月平均RMSE演變，不同于圖3，即使只考慮表層變量，局地化參數(shù)對(duì)于同化效果也有很大的影響。我們可以進(jìn)一步確認(rèn)c=0.2 rad會(huì)導(dǎo)致LPF的退化。綜合圖6和圖7的結(jié)果，本實(shí)驗(yàn)中最優(yōu)的LPF局地化參數(shù)c=0.1 rad。

圖6 不同局地化參數(shù) LPF 實(shí)驗(yàn)中區(qū)域平均（60°S～60°N，環(huán)地球）垂向均方根誤差Fig.6 Regional mean vertical root mean square error(60°S?60°N,ring the earth) in LPF experiments with different local parameters

圖7 不同局地化參數(shù) LPF 實(shí)驗(yàn)中區(qū)域平均（60°S?60°N，環(huán)地球）均方根誤差時(shí)間序列Fig.7 Regional mean (60°S?60°N,ring the earth) root mean square error time series in the experiments of LPF with different local parameters

3.3 EAKF 與 LPF 的同化效果比較

最后我們采用最優(yōu)局地化參數(shù)的EAKF與LPF進(jìn)行對(duì)比。圖8比較了EAKF與LPF的RMSE垂向分布。從溫度RMSE可以看出，當(dāng)深度小于100 m時(shí)，LPF（藍(lán)線）同化效果要優(yōu)于 EAKF（紫線），但當(dāng)深度在100～300 m時(shí)，EAKF同化效果要優(yōu)于LPF濾波器。從前面的討論我們很容易理解這個(gè)結(jié)果，這是因?yàn)長(zhǎng)PF需要啟用垂向局地化來避免深層的錯(cuò)誤更新而EAKF在沒有垂向局地化的條件下可以通過相關(guān)系數(shù)（即使很?。﹣碛行Ц律顚幼兞?。如果EAKF也使用ν=1 000 m/rad 的局地化，其在 100 m 以下的同化效果和LPF是一樣的（紅線）。對(duì)于鹽度變量，LPF的同化效果整體要略優(yōu)于EAKF，并且在淺層的效果更加明顯，這從一方面也說明溫度和鹽度的關(guān)系具有一定的非線性，因而粒子濾波器的表現(xiàn)更好。

圖8 最優(yōu)局地化參數(shù) EAKF、LPF 對(duì)比實(shí)驗(yàn)中區(qū)域平均（60°S～60°N，環(huán)地球）垂向均方根誤差Fig.8 Regional mean (60°S?60°N,ring the earth) vertical root mean square error in the experiments of EAKF and LPF with best local parameters

圖9為最優(yōu)局地化參數(shù)的EAKF和LPF實(shí)驗(yàn)的RMSE之差在海表的空間分布，其中差異由以下公式得到

圖9 最優(yōu)局地化參數(shù) EAKF、LPF 實(shí)驗(yàn)均方根誤差之差空間分布Fig.9 Spatial distribution of the difference between the root mean square error of the EAKF and LPF experiments

由定義，若它們之差大于0，則表示LPF的同化效果優(yōu)于EAKF；若他們之差小于0，則相反。圖中可以明顯的看出LPF對(duì)于SST的同化效果要優(yōu)于EAKF，除了北太平洋和赤道東太平洋的小部分區(qū)域外，LPF的誤差都相對(duì)較小。對(duì)于SSS和SSH來說，兩種濾波器同化效果互相有優(yōu)劣。在一些典型海區(qū)（如印度洋、大西洋），LPF對(duì)鹽度同化略有優(yōu)勢(shì)。而對(duì)于赤道太平洋，LPF對(duì)于SSH的同化效果更好。我們?cè)诒?進(jìn)一步計(jì)算了不同區(qū)域的平均RMSE以及同化影響（Influence of Assimilation,IOA）并進(jìn)行比較。其中IOA定義如下：

表3 最優(yōu)局地化參數(shù)EAKF、LPF不同區(qū)域平均均方根誤差（RMSE）、同化影響（IOA）對(duì)比Table 3 Comparison table of mean root mean square error (RMSE) and influence of Assimilation (IOA) of EAKF and LPF in different regions

顯然LPF的SST同化效果在所有海區(qū)都優(yōu)于EAKF。LPF對(duì)SSS的同化效果在印度洋和大西洋的優(yōu)勢(shì)也非常顯著。而對(duì)于SSH來說，EAKF的同化效果全面優(yōu)于LPF。實(shí)際上，我們發(fā)現(xiàn)使用集合成員有限的EAKF同化SST對(duì)于SSH的改進(jìn)效果本身就不是非常顯著，其結(jié)果受到了局地化以外的很多其他因素（如平均動(dòng)力地形MDT，集合離散度，跨變量相關(guān)系數(shù)等）的影響，因此還需要更多的研究來進(jìn)一步解釋兩種不同同化格式對(duì)于SSH的同化效果差異。

4 討論與展望

本文在大型地球系統(tǒng)模式CESM中開展了SST資料的同化實(shí)驗(yàn)?？疾炝司值鼗桨负蛥?shù)對(duì)于集合調(diào)整卡爾曼濾波器和局地化粒子濾波器兩種同化方法的效果的影響，并對(duì)兩種方法進(jìn)行了比較。

其中，集合卡爾曼濾波器使用變量之間的相關(guān)性來傳遞觀測(cè)點(diǎn)的更新量。為此我們提出了一個(gè)虛假相關(guān)占比公式，對(duì)集合樣本的相關(guān)和真實(shí)的相關(guān)進(jìn)行比較，判別局地化對(duì)于抑制遠(yuǎn)距離虛假相關(guān)的貢獻(xiàn)，并以此為依據(jù)選擇局地化參數(shù)。在垂直方向，我們發(fā)現(xiàn)由于深層變量的集合離散度較小，在同化SST資料的時(shí)候不需要特別采用垂向局地化就能達(dá)成效果。而局地化粒子濾波器基于使用局部觀測(cè)來決定不同位置變量的權(quán)重，并利用重采樣方法進(jìn)行更新，因此局地化參數(shù)對(duì)于粒子濾波器的影響更大。我們必須謹(jǐn)慎選擇局地化參數(shù)來避免粒子濾波器的退化，否則會(huì)產(chǎn)生負(fù)面的同化效果。此外我們也驗(yàn)證了垂直方向的局地化對(duì)于LPF是必要的。

選用最優(yōu)局地化方案的兩種濾波器的對(duì)比表明，LPF在淺層的同化效果比EAKF好，在較深層，由于垂向局地化的影響，LPF的效果略差于EAKF。由于粒子濾波器在理論上不含高斯分布假設(shè)，比較適用于非線性的系統(tǒng)和誤差，在實(shí)驗(yàn)中LPF對(duì)于鹽度的變量的同化優(yōu)勢(shì)也比較明顯。

綜上，本文驗(yàn)證了粒子濾波器在全球耦合氣候模式中的可行性，強(qiáng)調(diào)了有效的局地化策略對(duì)于避免濾波退化和改進(jìn)同化效果方面的重要性。通過與集合卡爾曼濾波器的同化效果對(duì)比，我們進(jìn)一步揭示了局地化粒子濾波器這一種新興算法在未來具有非常大的潛力被應(yīng)用于復(fù)雜模式的業(yè)務(wù)化同化系統(tǒng)中。

當(dāng)前的實(shí)驗(yàn)僅考慮了同化SST資料過程中的局地化方案，對(duì)于其他的觀測(cè)資料類型，如溫鹽廓線資料和海表面高度異常資料，其結(jié)論不一定相同。另一方面，當(dāng)前的工作僅考慮使用20個(gè)集合成員進(jìn)行計(jì)算，與局地化相關(guān)的結(jié)論也有一定的限制。我們將在下一階段的工作中，進(jìn)一步考慮使用不同的濾波器對(duì)于多源海洋觀測(cè)資料的同化效果差別以及局地化方案的選擇策略，也將進(jìn)一步嘗試不同集合成員數(shù)對(duì)于同化效果的影響。