李科，苑福利，劉廠

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

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

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

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

1 試驗方法

1.1 簡單耦合模式介紹

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

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

1.2 實驗設(shè)計

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

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

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

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

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

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

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

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

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

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

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

2 簡單耦合模式實驗結(jié)果

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

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

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

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

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

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

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

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

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

2.2 耦合模式比較實驗結(jié)果

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

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

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

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

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

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

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

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

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

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

3 結(jié)論與展望

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

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