李科，苑福利，劉廠

（1.海軍研究院，天津 300061；2.哈爾濱工程大學(xué) 智能科學(xué)與工程學(xué)院，黑龍江 哈爾濱 150001）

數(shù)據(jù)同化算法已被廣泛應(yīng)用于大氣和海洋環(huán)境的預(yù)報(bào)研究中，依靠模式狀態(tài)的概率密度函數(shù)（沈哲奇等，2016），對(duì)大氣和海洋信息進(jìn)行準(zhǔn)確估計(jì)。當(dāng)前數(shù)據(jù)同化主要使用的模式都是非線性的（白楊等，2020）。針對(duì)模式的非線性問題，數(shù)據(jù)同化方法一般采用兩種方案來解決（Chorin et al，2009），一種方法是使用高斯線性假設(shè)，如集合卡爾曼濾波數(shù)據(jù)同化方法（許立兵等，2020）。然而，實(shí)際應(yīng)用的模式都是高維非線性的，針對(duì)該模式傳統(tǒng)的線性化假設(shè)會(huì)產(chǎn)生截?cái)嗾`差，同時(shí)在海氣耦合模式中也存在局限性（Evensen et al，1994）。其中最有代表性的是傳統(tǒng)的EAKF 方法，使用蒙特卡洛方程積分控制的概率密度函數(shù)，進(jìn)一步描述先驗(yàn)的假設(shè)和誤差統(tǒng)計(jì)。如果模式是弱非線性時(shí)，EAKF 可以假設(shè)先驗(yàn)和似然密度分布都是高斯分布假設(shè)（Anderson et al，2001）。但如果系統(tǒng)是高度非線性的，這些迭代映射很難收斂，會(huì)影響協(xié)方差的準(zhǔn)確性。因此，EAKF 在高度非線性的模式中具有一定的局限性（Freitas et al，2001）。針對(duì)該方法的局限性，近幾年研究人員更多地將該方法與變分同化方法結(jié)合，構(gòu)建集合變分方案，以En4Dvar為代表的混合變分方法（Tian et al，2018）可以較好地適應(yīng)非高斯模式的要求，對(duì)于Kalman 的另一種改進(jìn)方法是針對(duì)實(shí)際應(yīng)用的需求構(gòu)建局地化方案，其中Penny 等（2013）成功將LETKF 應(yīng)用于在小尺度建模（COSMO）框架中。

另一種解決非線性問題的方法是粒子濾波方法（PF）。該方法基于統(tǒng)計(jì)和貝葉斯公式發(fā)展而來（Van Leeuwen et al，1996），同時(shí)結(jié)合先驗(yàn)和后驗(yàn)概率密度分布來更新粒子狀態(tài)。該方法在同化過程中每個(gè)粒子將根據(jù)觀測(cè)信息被賦予一個(gè)權(quán)重，該權(quán)重與觀測(cè)的概率成比例，這是粒子濾波獨(dú)有的特點(diǎn)（Van Leeuwen，2011）。但該方法也存在局限性，其中最為主要的就是“粒子退化”和“維數(shù)災(zāi)難”問題。針對(duì)粒子濾波自身存在的問題，對(duì)其改進(jìn)研究也主要集中于兩個(gè)方面：第一是采用混合方法，使用粒子濾波過濾狀態(tài)空間的強(qiáng)非高斯部分，同時(shí)在高維狀態(tài)使用卡爾曼濾波，該方法同時(shí)結(jié)合卡爾曼濾波與粒子濾波的優(yōu)點(diǎn)，最有代表性的是Slivinskiet 等（2015）提出的混合粒子集合卡爾曼濾波。另一種方法是通過提議密度調(diào)整粒子防止粒子崩潰的發(fā)生，該方法的主要代表是Van Leeuwen（2011）提出了均權(quán)重粒子濾波（EWPF），該方法使用提議密度引導(dǎo)所有粒子靠近觀測(cè)，同時(shí)提高粒子濾波的同化質(zhì)量。而且Ades 等（2015） 和Chen 等（2020）分別將該方法進(jìn)一步改進(jìn)局地化方案并應(yīng)用于高維模式中。

目前數(shù)據(jù)同化方法主要應(yīng)用于耦合模式中。因此本文在一個(gè)簡單的非線性耦合模型上比較兩種算法的同化特點(diǎn)，進(jìn)一步研究它們的同化特性。

1 試驗(yàn)方法

1.1 簡單耦合模式介紹

在實(shí)際應(yīng)用中，數(shù)據(jù)同化方法主要應(yīng)用于耦合模式中，如CGCM 等，然而該動(dòng)力學(xué)模式需要大量的計(jì)算資源支持，因此本文并沒有使用復(fù)雜的CGCM 模式進(jìn)行實(shí)驗(yàn)。Zhang 等（2013）在Lorenz 63 模型的基礎(chǔ)上開發(fā)了簡單耦合預(yù)測(cè)模型，該模型具有CGCM（Zhang et al，2012） 的基本特征，方程為：

該模型有五個(gè)變量，x1，x2和x3代表大氣，w代表上層海洋，濁代表深層海洋。Zhang 等（2014）和Han 等（2014）詳細(xì)介紹了參數(shù)的意義。大氣模式誤差為0.5，上層海洋模式誤差為0.13，深海模式誤差為0.01。

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

本研究使用“孿生”實(shí)驗(yàn)框架設(shè)計(jì)實(shí)驗(yàn)，根據(jù)給定初始條件進(jìn)行10 000 次迭代進(jìn)行spin-up 擾動(dòng)，然后積分10 000 次迭代產(chǎn)生真值。觀測(cè)信息是由真值疊加白噪聲近似，兩個(gè)算法的集合成員數(shù)都是20。

實(shí)驗(yàn)框架設(shè)計(jì)參考了Zhang 等（2013）和Ades等（2015）的研究。該耦合模型初始條件為（0，1，0， 0，0），觀測(cè)標(biāo)準(zhǔn)差，大氣為2，上層海洋為0.5，由于深海缺乏觀測(cè)，因此本實(shí)驗(yàn)假設(shè)深海沒有觀測(cè)值。根據(jù)Ades 等（2015）的實(shí)驗(yàn)設(shè)計(jì)，確定大氣的同化觀測(cè)間隔為45，上層海洋為2 356。實(shí)驗(yàn)積分100 000 次迭代，再根據(jù)最后的10 000 次迭代研究兩種方法的同化特性。

均方根誤差（root-mean-square-error，RMSE）是評(píng)價(jià)同化可靠性的重要標(biāo)準(zhǔn)，計(jì)算集合均值（x，y和z）與真值（xt，yt和zt）的差值（RMSET）：

本文采用兩種均方根誤差對(duì)同化結(jié)果進(jìn)行評(píng)價(jià)。

1.3 集合調(diào)整卡爾曼濾波

在本文中，假設(shè)模式在時(shí)刻t 的狀態(tài)為xt，觀測(cè)狀態(tài)為yt。基于模型方程產(chǎn)生觀測(cè)值（Zhang et al，2011），其中觀測(cè)誤差假設(shè)為高斯：

1.4 均權(quán)重粒子濾波

對(duì)于數(shù)據(jù)同化的最優(yōu)估計(jì)問題，可以近似理解為對(duì)于后驗(yàn)概率密度分布的估計(jì)。貝葉斯定理提供了一種在先驗(yàn)概率基礎(chǔ)上，計(jì)算后驗(yàn)概率的方法。粒子濾波方法就是基于貝葉斯定理而來，它利用條件概率密度分布函數(shù)表示同化結(jié)果（Ades et al，2015）。根據(jù)觀察，條件概率密度函數(shù)可以表示為（Van Leeuwen，2015）：

其中，p（x）為先驗(yàn)概率密度函數(shù)，p（y|x）為條件概率密度，p（x|y）為后驗(yàn)概率密度。

EWPF 由Van Leewen（2015）詳細(xì)闡述，該方法使用提議密度保證集合中大多數(shù)粒子可以靠近觀測(cè)信息，根據(jù)調(diào)整后的集合粒子確定后驗(yàn)概率密度：

其中子=（tj-t0）/（tn-t0）表示分析時(shí)間到觀測(cè)時(shí)間的歸一化時(shí)間尺度，t0和tn分別表示前一次觀測(cè)時(shí)刻和后一次觀測(cè)時(shí)刻。松弛項(xiàng)的函數(shù)B（子）是一個(gè)比例因子，控制對(duì)觀測(cè)的松弛度。第二步是等權(quán)重調(diào)整，EWPF 利用松弛提議密度來調(diào)整靠近觀測(cè)的粒子，保證集合中粒子的權(quán)重幾乎相等，粒子權(quán)重可以表示為：

2 簡單耦合模式實(shí)驗(yàn)結(jié)果

2.1 確定膨脹因子、松弛系數(shù)和時(shí)間尺度子

由于模態(tài)誤差和采樣誤差的存在，對(duì)于EAKF方法中存在濾波發(fā)散的問題，使用膨脹因子避免濾波發(fā)散。而EWPF 使用時(shí)間尺度子和比例因子b控制對(duì)粒子趨向觀測(cè)的松弛調(diào)整，從而調(diào)節(jié)提議密度調(diào)整粒子。因此，確定合適的膨脹因子和松弛系數(shù)可以有效地提高兩種方法的同化質(zhì)量。

根據(jù)上文介紹，從圖1（a） 中可以看出當(dāng)EAKF 中的大氣膨脹因子為4 580.5 時(shí)，大氣RMSE 最低。而從圖1（b）可以得出，根據(jù)大氣膨脹系數(shù)決定海洋膨脹系數(shù)，當(dāng)EAKF 海洋膨脹系數(shù)為43.5 時(shí)，RMSE 最低。

圖1 耦合模式中EAKF 中（a）大氣膨脹因子和（b）海洋膨脹因子變化的RMSE 值

使用RMSE 評(píng)價(jià)EWPF 中的松弛比例因子對(duì)于同化結(jié)果的影響。由圖2 可知，當(dāng)EWPF 中的大氣松弛參數(shù)為4，EWPF 中的海洋松弛參數(shù)為6時(shí)，大氣和海洋同時(shí)RMSE 處于最低。

圖2 在耦合模式中使用EWPF，松弛比例因子與RMSE 之間的關(guān)系

子的選取決定了在觀測(cè)間隔中使用提議密度的時(shí)刻，針對(duì)該耦合模式中大氣和海洋的不同觀測(cè)頻率，分別選擇獲得最優(yōu)大氣海洋狀態(tài)RMSET的時(shí)間尺度子，繪制觀測(cè)頻率與時(shí)間尺度子的關(guān)系，如圖3 所示。

圖3 在耦合模式中使用EWPF 不同的觀測(cè)頻率與子的關(guān)系

由圖3 可以看出，隨著觀測(cè)頻率的增加，子的值逐漸降低。由公式可知，子越小表示開始使用提議密度調(diào)整粒子越早，造成該結(jié)果的主要原因是隨著同化觀測(cè)之間的間隔增加，為了保證粒子更好地靠近觀測(cè)，需要在觀測(cè)間隙更多地獲取觀測(cè)信息，更好地調(diào)整粒子狀態(tài)，因此子的取值變小，是為了更好地調(diào)整粒子狀態(tài)。

2.2 耦合模式比較實(shí)驗(yàn)結(jié)果

在海氣耦合模式中，針對(duì)不同大氣和海洋觀測(cè)頻率對(duì)于同化結(jié)果的影響，比較分析不同觀測(cè)頻率下大氣和海洋狀態(tài)的RMSET，結(jié)果如圖4 所示。圖4（a）和圖4（b）給出了大氣觀測(cè)頻率為20，海洋觀測(cè)頻率為1 500 時(shí)的RMSET，圖4（c） 和圖4（d）給出了大氣觀測(cè)頻率為35，海洋觀測(cè)頻率為2 000 的結(jié)果，圖4（e）和圖4（f）給出了大氣觀測(cè)頻率為45，海洋觀測(cè)頻率為2 356 的RMSET，由圖中可以看出兩種方法在不同大氣觀測(cè)頻率下結(jié)果較為接近，隨著觀測(cè)頻率增加，EWPF結(jié)果略優(yōu)于EAKF 方法。對(duì)于上層海洋狀態(tài)，EWPF 的同化結(jié)果在三種不同海洋觀測(cè)條件下，其同化結(jié)果明顯優(yōu)于EAKF。在后續(xù)實(shí)驗(yàn)中將選擇大氣觀測(cè)頻率為45，海洋觀測(cè)為2 356 繼續(xù)實(shí)驗(yàn)，因?yàn)樵撚^測(cè)頻率可以更好地體現(xiàn)模式的非線性。

圖4 兩種算法的耦合模型根據(jù)不同觀測(cè)頻率的RMSET

在海氣耦合模式中，圖5（a） 描述了在海氣耦合模式中選擇在最后10 000 次迭代的統(tǒng)計(jì)結(jié)果，比較兩種算法對(duì)于大氣狀態(tài)的估計(jì)，其中黑色菱形表示觀測(cè)值，黑色點(diǎn)虛線表示真值。從圖中可以看出EWPF（點(diǎn)劃線）總優(yōu)于EAKF（實(shí)線），而且EWPF 的軌跡更加趨近于觀測(cè)值。圖5（b） 給出了在60 000～100 000 模式迭代過程中兩種算法對(duì)于上層海洋的狀態(tài)估計(jì)，同樣的EWPF 結(jié)果優(yōu)于EAKF 方法。從兩組實(shí)驗(yàn)中可以看出，EAKF 結(jié)果明顯偏離真值，僅在同化時(shí)刻快速變化趨近真值，其原因是在實(shí)驗(yàn)中大氣和海洋觀測(cè)間隔選擇較大，導(dǎo)致有效觀測(cè)減少。而EWPF 可以使用提議密度來調(diào)整粒子，使其接近同化時(shí)刻前的觀測(cè)，但EAKF 只在同化時(shí)刻使用觀測(cè)值來調(diào)整當(dāng)時(shí)的集合。因此，EAKF 狀態(tài)估計(jì)在所有迭代中都有明顯的突變，并在同化時(shí)迅速接近真值。為了進(jìn)一步分析這一問題及其同化特征，圖6 描述了兩種同化方法的RMSET和RMSEO。

圖5 在耦合模式下（a）給出大氣x1 和（b）給出上層海洋兩種方法的同化結(jié)果

根據(jù)模式狀態(tài)統(tǒng)計(jì)結(jié)果，圖6 給出了對(duì)于不同狀態(tài)估計(jì)的RMSET和RMSEO，其中圖6（a） 和（b）給出全時(shí)間序列的RMSET，EWPF（點(diǎn)劃線）整體低于EAKF（實(shí)線），在同化時(shí)刻（黑點(diǎn)）EAKF 的RMSET明顯下降，這一結(jié)果也印證了先前對(duì)EAKF 的描述，也證明了EAKF 狀態(tài)估計(jì)與真值接近。圖6（c）和（d）給出了兩種同化方法在耦合模式下同化結(jié)果的RMSEO，可以清晰地看出EAKF 和EWPF 有相似的結(jié)果。EWPF 上層海洋低于EAKF，基于以上結(jié)果可以得出EWPF 結(jié)果更接近觀測(cè)。

基于圖6 對(duì)RMSE 描述，使用概率密度分布直方圖進(jìn)一步研究兩種方法集合粒子的分布，結(jié)果如圖7 所示。圖7（a）和（c）給出了EAKF 中20個(gè)集合成員對(duì)于大氣和海洋狀態(tài)的概率分布，圖7（b）和（d）給出了EWPF 中20 粒子在大氣和海洋的概率密度分布，其中大氣的實(shí)驗(yàn)選擇在98 100時(shí)間步長的結(jié)果，海洋的實(shí)驗(yàn)選擇在87 172 時(shí)間步長的結(jié)果。從圖中可以得出，大部分EAKF 集合成員都分布在真值和觀測(cè)值之間，而EWPF 粒子與觀測(cè)值更接近。原因是EWPF 使用建議密度來調(diào)整粒子接近觀測(cè)，而觀測(cè)值對(duì)于最終確定粒子權(quán)重有著重要的意義。EAKF 的最終結(jié)果由觀測(cè)增量決定，因此其估計(jì)值與真實(shí)值更接近。

圖6 使用這兩種算法在耦合模型進(jìn)行不同實(shí)驗(yàn)條件下大氣海洋對(duì)于真值和觀測(cè)的RMSE 結(jié)果

圖7 卡爾曼濾波和粒子濾波使用20 集合粒子數(shù)的概率密度分布圖

在海氣耦合模型中，假設(shè)深海沒有可用的觀測(cè)值，但通過公式（4）可知，深海的狀態(tài)可以根據(jù)其他狀態(tài)的觀測(cè)來調(diào)整。圖8 給出了兩種同化方法在深海狀態(tài)下，時(shí)間步長為92 000～100 000 時(shí)的RMSET。由于深海沒有觀測(cè)數(shù)據(jù)用于同化實(shí)驗(yàn)，因此最終結(jié)果僅由積分公式（1），通過同化調(diào)整后的大氣和海洋狀態(tài)來獲得，從圖中可以清晰地看出，EWPF 的RMSET明顯低于EAKF 方法。說明在缺少觀測(cè)的情況下，單純依靠其他狀態(tài)觀測(cè)通過積分方程調(diào)整時(shí)，EWPF 方法可以更好地獲得同化結(jié)果。

圖8 深層海洋在時(shí)間序列92 000 到100 000 時(shí)間步長的RMSET

3 結(jié)論與展望

本文選擇在非線性海氣耦合模式上，比較分析EAKF 和EWPF 兩種同化方法的同化特點(diǎn)。實(shí)驗(yàn)結(jié)果表明，在非線性耦合模式中，EWPF 模式狀態(tài)估計(jì)明顯優(yōu)于EAKF，特別是在觀測(cè)時(shí)間間隔較大的情況下，EWPF 可以有效避免EAKF 在同化時(shí)刻突然靠近模式真值，同化結(jié)果更加穩(wěn)定。其次，基于兩種不同的RMSE 和集合粒子分布直方圖，可以看出EWPF 的同化結(jié)果更加接近觀測(cè)值，而EAKF方法的同化結(jié)果更加接近狀態(tài)真值。因此，EWPF同化結(jié)果的準(zhǔn)確性與觀測(cè)的質(zhì)量有關(guān)。最后，在海氣耦合中，假設(shè)深海中沒有可用觀測(cè)值，深海狀態(tài)只受到積分方程中其他變量同化結(jié)果的影響，通過比較無觀測(cè)深海狀態(tài)RMSE 結(jié)果表明，EWPF 對(duì)于耦合模式有著更好的應(yīng)用，對(duì)于缺少觀測(cè)的狀態(tài)，可以根據(jù)已有觀測(cè)較好地調(diào)節(jié)。這是EWPF 的另一個(gè)特點(diǎn)，這些實(shí)驗(yàn)也證明了基于粒子濾波改進(jìn)的EWPF 在該耦合模型中具有很好的應(yīng)用潛力。

盡管EWPF 在簡單模型上得到有效證明，但該方法缺少在中尺度耦合模式中與EAKF 方法的進(jìn)一步驗(yàn)證，其特性和同化質(zhì)量還需要更多的實(shí)驗(yàn)研究。