沈浙奇,唐佑民,2 *,高艷秋

（1 .國家海洋局第二海洋研究所衛(wèi)星海洋環(huán)境動力學國家重點實驗室,浙江杭州310012；2 .北大不列顛哥倫比亞大學環(huán)境科學及工程系,加拿大喬治王子城V2 N 4Z9）

?

集合資料同化方法的理論框架及其在海洋資料同化的研究展望

沈浙奇1,唐佑民1,2 *,高艷秋1

（1 .國家海洋局第二海洋研究所衛(wèi)星海洋環(huán)境動力學國家重點實驗室,浙江杭州310012；2 .北大不列顛哥倫比亞大學環(huán)境科學及工程系,加拿大喬治王子城V2 N 4Z9）

摘要:在海洋動力系統(tǒng)的數(shù)值模擬中,海洋資料同化是一種能夠有效融合多源海洋觀測資料和數(shù)值模式的方法。它不僅可以顯著地提高數(shù)值模擬的效果,構(gòu)造海洋再分析資料場,還能有效減少海洋和氣候預報時模式初始條件的不確定性。因此,海洋資料同化對于海洋研究和業(yè)務化應用具有非常重要的意義。資料同化方法的研究一直是大氣、海洋科學的熱門課題之一。其中,集合卡爾曼濾波器（En K E）是一種有效的資料同化方法,自提出以來經(jīng)過了20多年的發(fā)展和改進,已經(jīng)在海洋資料同化中得到了廣泛的研究和應用。近年來,隨著動力模式的不斷發(fā)展和計算能力的提高,粒子濾波器由于不受模型線性和誤差高斯分布假設的約束,也逐漸成為了當前資料同化方法研究的熱點。本文分析和總結(jié)了目前關于集合卡爾曼濾波器和粒子濾波器的一些最新理論研究結(jié)果,在貝葉斯濾波理論的框架下討論了這兩類算法的關聯(lián)和區(qū)別,以及各自在資料同化實踐中的優(yōu)勢和不足。在此基礎上,我們探討了粒子濾波器應用于海洋模式資料同化的主要困難和目前可行的一些解決方法,展望了集合資料同化方法研究的新趨勢,為集合資料同化方法的進一步發(fā)展和應用提供理論基礎。

關鍵詞:資料同化；集合卡爾曼濾波器；粒子濾波器；非高斯噪聲；貝葉斯濾波

1 引言

資料同化（Data Assimilation,也稱數(shù)據(jù)同化）是一種能夠有機結(jié)合數(shù)值模式和觀測這兩種海洋學研究基本手段的方法。資料同化的核心思想是在動力學模式的運行過程中不斷地融合新的觀測信息,以達成改進模式初值、優(yōu)化模式參數(shù)、改進模式分析和預報效果的目的。通過引入資料同化,不僅能夠有效地提高海洋模擬效果,減少海洋和氣候預報初始條件的不確定性,還能夠為缺少海洋觀測的深海和邊緣海提供再分析資料［1—2］,為海洋觀測計劃和數(shù)值預報模式物理量及參數(shù)等提供設計依據(jù)。因此海洋資料同化在現(xiàn)代物理海洋學的研究中占據(jù)了非常重要的地位［3—4］。

目前,資料同化方法主要分為兩類:變分方法（Variational method）和序貫方法（Sequential method）。變分方法包括三維變分（3 D-VA R）和四維變分（4 D-VA R）［5—6］,其主要特點是利用數(shù)值最優(yōu)化算法求解目標函數(shù)的最小化問題。而以集合卡爾曼濾波器（EnkE）［7］為代表的序貫方法則是基于統(tǒng)計估計理論［8］的遞歸方法。所謂的序貫（一般也翻譯為順序）,指的是在每次有觀測的時刻就利用顯式的更新公式進行一次同化,然后再將所得的分析場作為下一次預報的初始場。這兩類方法在大氣、海洋等諸領域都已經(jīng)取得了巨大的成功［9—11］。這里我們主要關注集合卡爾曼濾波器［7,12］,相比于變分方法,其優(yōu)點在于無需發(fā)展四維變分同化算法要求的切線性和伴隨模式,在同化時也考慮預報誤差的動力演變,并且可以顯式地提供集合預報的初始擾動,因此近年來隨著計算條件的改善已得到了廣泛的應用,也成為當前資料同化方法研究的熱點［13—14］。

另一方面,粒子濾波器（Particle Eilter）也是一種集合資料同化方法,它在動力系統(tǒng)的狀態(tài)估計中有著非常廣泛的應用。相比于集合卡爾曼濾波器,粒子濾波器算法不受模型狀態(tài)變量和誤差高斯分布假設的約束,可以用于任意非線性非高斯的動力系統(tǒng)。它通過采樣方法近似表示模式狀態(tài)場的概率分布密度函數(shù),能更好地同化非高斯信息。但是受限于濾波退化等問題,粒子濾波器一直沒有在海洋資料同化領域得到較大關注。近年來,隨著粒子濾波器被成功應用于一些高維動力系統(tǒng)的資料同化［15—16］,粒子濾波器開始成為當前資料同化算法研究的一個熱門方向［17］。

本文在第二部分簡要介紹了集合卡爾曼濾波器和粒子濾波器的研究背景。第三部分在貝葉斯濾波的統(tǒng)一框架下推導了集合卡爾曼濾波器和粒子濾波器的算法,并從理論和實踐兩方面指出了兩者的優(yōu)勢和不足。第四部分主要介紹了粒子濾波器算法研究的最新進展和主要挑戰(zhàn),最后則是對集合資料同化方法的總結(jié)和展望。

2 集合卡爾曼濾波器和粒子濾波器的研究背景

集合卡爾曼濾波器是在卡爾曼濾波器的基礎上發(fā)展起來的。對于線性的模式和高斯分布的噪聲系統(tǒng),經(jīng)典的卡爾曼濾波器給出了精確的最優(yōu)解［18］。而對于非線性的模式,最直接的處理方法是將卡爾曼濾波器應用于該模式的一階線性化近似,得到相應的近似解,這種方法被稱為擴展卡爾曼濾波器（E K E）。實際上,無論是經(jīng)典的卡爾曼濾波器還是擴展卡爾曼濾波器,都需要在更新系統(tǒng)狀態(tài)場的同時,計算和保存背景場的誤差協(xié)方差矩陣,這在狀態(tài)場空間維數(shù)非常巨大的海洋模式資料同化中是難以實現(xiàn)的。同時,擴展卡爾曼濾波器中的線性化過程會引起很大的截斷誤差,造成同化結(jié)果的偏差。這幾點都極大地妨礙了擴展卡爾曼濾波器在實際海洋資料同化中的應用。

Evensen于1994年提出了集合卡爾曼濾波器［7］,通過采用蒙特卡洛的思想來避免模式線性化和直接計算預報誤差協(xié)方差的困難。集合卡爾曼濾波器在算法上的最大突破是引入了集合的概念,即通過多次積分模式來得到預報集合,并使用預報集合的經(jīng)驗協(xié)方差來更新所有的集合成員,得到分析集合。經(jīng)驗協(xié)方差是一個統(tǒng)計上的概念,指的是根據(jù)有限的樣本計算得到的協(xié)方差矩陣。當集合成員的數(shù)目趨于無窮的時候,可以證明預報集合的經(jīng)驗協(xié)方差趨向于預報誤差協(xié)方差的真值。一方面,集合卡爾曼濾波器通過積分完整的非線性模式來進行預報,避免了截斷誤差的產(chǎn)生；另一方面,集合卡爾曼濾波器使用實時計算的經(jīng)驗協(xié)方差矩陣來代替卡爾曼濾波器公式中出現(xiàn)的預報誤差協(xié)方差矩陣,避免了協(xié)方差矩陣的儲存和更新,以及伴隨而來的巨大運算量。由于集合卡爾曼濾波器的提出,卡爾曼濾波器的思想得以應用于大型海洋模式的資料同化。在此之后,基于相同的集合化思想,一系列集合卡爾曼濾波器的衍生算法也相繼被提出,如集合轉(zhuǎn)移卡爾曼濾波器（E T K E）［19］,集合調(diào)整卡爾曼濾波器（E A K E）［20］,集合平方根濾波器（En-SRE）［21］等。這些方法雖然基于不同的角度來推導狀態(tài)場和誤差協(xié)方差矩陣的更新公式,但歸根結(jié)底還是都使用了經(jīng)典卡爾曼濾波器。

在此基礎上,一些新的方法和技巧也陸續(xù)被引入來改善集合卡爾曼濾波器及其衍生方法在實際同化中的表現(xiàn)。比如,協(xié)方差松弛（covariance inflation,也被稱為協(xié)方差膨脹）技術被引入來改變集合的離散度,從而解決因集合成員數(shù)目較少或者模式誤差的影響而產(chǎn)生的協(xié)方差矩陣低估問題［22］。協(xié)方差的松弛系數(shù)最初使用的是根據(jù)經(jīng)驗給定的固定參數(shù)。但在近幾年,一系列選擇自適應松弛系數(shù)的方法也得到了很大的發(fā)展［23—24］。在集合卡爾曼濾波器的同化實踐中,模式的初值、大氣外強迫、模式開邊界和河流輸入等因素的擾動會在很大程度上影響集合的離散程度,一些研究也表明使用這些擾動構(gòu)成預報集合能大大改善集合卡爾曼濾波器的同化效果［25］。

集合卡爾曼濾波器及其衍生方法能成功應用于大型海洋模式,另一個很主要的原因是局地化策略［26］的應用。一般來說,對于空間上不穩(wěn)定的動力系統(tǒng),分析格式并不能修正距觀測距離較遠的格點處的誤差,相反,由于分析空間范圍過大,反而會在距離較大的格點之間產(chǎn)生虛假的相關性。局地化策略相當于在觀測周圍區(qū)域內(nèi)的一個低維子空間內(nèi)進行分析,使得觀測不再影響給定距離之外的格點。通過引入局地化,除了能夠避免上述的虛假相關性以外,還能大大減少預報誤差協(xié)方差矩陣的非零元素數(shù)目,降低分析的運算量。所以,局地化方案的選取和局地化參數(shù)的設置對于集合卡爾曼濾波器的同化效果具有很大的影響,其理論研究在同化領域仍然是熱門方向［27—28］。

集合卡爾曼濾波器目前在國際上得到了廣泛的研究和應用。加拿大氣象中心于2005年最先在其全球集合預報系統(tǒng)中采用了業(yè)務化的集合卡爾曼濾波器來提供預報的初始條件［29］。W W RP/T H O RPE X國際計劃于2008年在阿根廷的布宜諾斯艾利斯開展了“比較集合卡爾曼濾波器和四維變分”的專題研討。該研討會通過將集合卡爾曼濾波器與當前主流業(yè)務化預報系統(tǒng)中采用的四維變分方法進行比較,進一步地闡明了集合卡爾曼濾波器方法的業(yè)務化前景［30］。另外,挪威2010年基于H Y C O M海洋模式,利用En-K E方法建立了北大西洋海洋預報系統(tǒng)［31］。美國N O A A/G ED L基于M O M 4模式,利用En K E方法建立了全球氣候分析與預報系統(tǒng)［32］。

雖然集合卡爾曼濾波器及其衍生算法都得到了廣泛的關注和認可,但是有必要指出卡爾曼濾波公式本身采用的一些假定實際上限制了上述方法的同化效果。如前述,只有在保證數(shù)值模式為線性模式,并且預報誤差和觀測誤差都是高斯分布的前提下,卡爾曼濾波器才提供最優(yōu)的解。卡爾曼濾波器的公式本身只包含狀態(tài)場的平均值和協(xié)方差信息,它實際上默認了公式中涉及到的所有誤差（預報誤差和觀測誤差）都是呈高斯分布的。而一個簡單的事實是,模式的非線性會不可避免地將誤差的分布變成非高斯分布［33］。對集合卡爾曼濾波器的漸近分析也表明,即使同化系統(tǒng)中的預報初值誤差和觀測誤差都是高斯分布的,分析集合的經(jīng)驗分布也會趨向于一個錯誤的極限——即與一般的貝葉斯濾波器的極限不同（貝葉斯濾波理論將在第三部分展開討論）。換句話說,即使我們的計算資源無限,集合卡爾曼濾波器也無法保證獲得最可靠的分析。在此背景下,集合卡爾曼濾波器以外的非線性集合同化方法在近些年得到了迅速發(fā)展［33—34］,其中一類典型的方法就是粒子濾波器。

相比于集合卡爾曼濾波器,粒子濾波器不受模型狀態(tài)變量和誤差高斯分布假設的約束,可以用于任意非線性非高斯的動力系統(tǒng)。因此粒子濾波器從21世紀初開始就被廣泛應用于地球系統(tǒng)數(shù)據(jù)同化［35—36］,并成為了一個熱門的研究方向。實際上,粒子濾波器并不是一種很新的概念。早在20世紀50 - 60年代, Ham mersley等就提出了一種稱為序貫重要性采樣的蒙特卡洛方法［37］,它的基本原理是通過離散的隨機樣本逼近狀態(tài)場的概率分布密度函數(shù)。但是,由于計算復雜性和樣本退化的問題,該方法在相當長的一段時間內(nèi)沒有取得多大的進展。直到1993年Gordon等提出了重取樣的概念［38］,在一定程度上克服了算法的退化問題,才出現(xiàn)了第一個真正可行的序貫蒙特卡洛方法。該方法稱為自助式粒子濾波器,被公認為是現(xiàn)代粒子濾波器的基石［39］。所謂的粒子實際上就是用來近似表示概率分布密度函數(shù)的樣本,這和集合卡爾曼濾波器中的集合成員是相對應的概念。粒子濾波器利用觀測來更新粒子的權(quán)重,并通過加權(quán)來逼近后驗概率分布密度。最初的序貫重要性采樣使用轉(zhuǎn)移概率密度函數(shù)作為重要性采樣函數(shù),在每一次觀測更新的時刻利用遞推關系來改變各個粒子的權(quán)重,但不改變粒子本身的數(shù)值。隨著迭代次數(shù)的增加,幾乎肯定會出現(xiàn)粒子退化的問題。粒子退化指的是這樣一種情況:隨著同化的進行,少數(shù)粒子逐漸占據(jù)了幾乎所有的權(quán)重,而大部分粒子的權(quán)重則基本可以忽略不計。如果出現(xiàn)了粒子退化,大量的計算成本會被用于更新對概率分布密度函數(shù)的貢獻幾乎為零的粒子,造成無法承受的運算負擔。而過低的有效樣本率也會導致后驗估計的精度完全不可靠,無法改善預報的效果。自助式粒子濾波器實際上就是在序貫重要性采樣的基礎上加入了重取樣過程,因此也被稱為重取樣粒子濾波器（SIR-PE）。重取樣的基本思想是根據(jù)后驗分布中的各粒子權(quán)重重新分配粒子集,將權(quán)重較大的粒子復制多份,同時消除權(quán)重很小的粒子。在這個過程中粒子的總數(shù)保持不變,最后得到的是等權(quán)重的粒子集。重取樣的依據(jù)是粒子的權(quán)重,具體的實施方法有很多種,包括多項式重取樣、殘差重取樣、系統(tǒng)重取樣等［40］。總之,通過重取樣,可以消除權(quán)重較小的粒子,避免退化的發(fā)生。但是這樣的處理也會導致重取樣后的粒子幾乎都是少數(shù)幾個粒子的后代,粒子的多樣性明顯減弱,不能充分表征后驗概率密度函數(shù),該現(xiàn)象稱為粒子貧化問題。為解決貧化問題,Kotecha和Djuric提出了高斯粒子濾波器（Gaussian particle filter）及高斯和粒子濾波器（Gaussian su m particle filter）［41—42］,Xiong等提出了后驗高斯重取樣的粒子濾波器［43］,Nakano等提出了融合粒子濾波器［44］,這些方法對于抑制貧化都起到了一定的效果。

對于維數(shù)較小的模式和同化系統(tǒng),重取樣過程能夠解決粒子的權(quán)重退化問題。但對于海洋模式資料同化中涉及的高維模式,結(jié)果并沒有那么樂觀。Snyder等系統(tǒng)地研究了粒子濾波中的退化問題,并且指出:為了防止粒子的權(quán)重發(fā)生退化,粒子的數(shù)目必須隨著狀態(tài)空間的維數(shù)而指數(shù)增長,即使是引入了重取樣過程,也無法從根本上解決粒子數(shù)目需要隨系統(tǒng)維數(shù)增長的難題［45］。對于海洋模式資料同化中的高維模式來說,這樣的計算量顯然是難以接受的。簡而言之,“維數(shù)災難”（curse of dimensionality）［46］問題是粒子濾波器應用于實際大氣海洋研究中的高維模式的最大障礙。Van Leeuwen詳細總結(jié)了高維模式系統(tǒng)資料同化中的粒子濾波算法,在重申了維數(shù)困難的前提下,也列出了幾種已知可行的粒子濾波器,如引導粒子濾波器、邊際粒子濾波器、輔助粒子濾波器等［36］。大體上來看,目前主要的努力方向仍然在序貫重要性取樣和重取樣的結(jié)合方法上面。

最初提出的序貫重要性取樣使用轉(zhuǎn)移密度函數(shù)作為重要性取樣函數(shù)（也叫做建議分布密度）,也就是說直接利用前一步的粒子進行預報,將預報的結(jié)果不加處理地作為新的粒子。這種情況下,算法本身不改變粒子的數(shù)值,而只是利用粒子與觀測的相近程度更新權(quán)重。當狀態(tài)場空間的維數(shù)很大時,很容易由于粒子與觀測的差別太大而無法獲取有意義的權(quán)重,這就是退化現(xiàn)象。相對地,如果使用一個依賴于觀測值的建議分布來產(chǎn)生新的粒子,就可以人為地使粒子更接近觀測值,避免計算得到的權(quán)重過小,在一定程度上也能夠降低退化發(fā)生的幾率?；诖饲疤峥梢缘玫阶顑?yōu)建議分布。需要注意的是,這里的“最優(yōu)”并不是指同化的效果最優(yōu),而是指根據(jù)這種建議分布下產(chǎn)生的粒子計算出來的權(quán)重的方差最小,能最大限度地避免退化的發(fā)生［47］。在絕大多數(shù)情況下,這種最優(yōu)建議分布無法解析求得,目前熱門的方法都是在試圖去近似表示這個最優(yōu)建議分布。這些方法在第四部分具體闡述。

3 貝葉斯理論框架下的集合卡爾曼濾波和粒子濾波

資料同化在數(shù)學上是動力學系統(tǒng)的狀態(tài)估計問題,由狀態(tài)空間方法描述。狀態(tài)空間模型分為狀態(tài)預報模型和觀測模型兩部分,二者在資料同化系統(tǒng)中又分別被稱為模式算子和觀測算子。模式算子對應著數(shù)值模式,在不同的動力學系統(tǒng)中可以分別是大氣模式、海洋模式、陸面過程模式或水文模式等。模式算子描述的是狀態(tài)變量場隨時間的演變,一般情況下可以將它理解成為一個“黑匣子”,利用給定的初值和外強迫場得到具有一定準確性的預報場。觀測算子則描述了狀態(tài)場和觀測資料的對應關系:對于直接觀測來說,一般需要使用空間插值的方法將模式的網(wǎng)格點投影到觀測點；而對于遙感等非直接觀測來說,也可能需要使用非線性的觀測模型來將狀態(tài)變量和觀測變量相關聯(lián),比如遙感輻射傳輸模型。狀態(tài)空間模型可以用以下公式描述:

式中,xk代表離散時間tk上的狀態(tài)變量場,yk代表相同時刻的對該狀態(tài)場的觀測。η和ζ分別代表了隨機的模式誤差和觀測誤差,二者的概率分布密度函數(shù)可以是任意形式的（但在很多情況下兩者都被假設為高斯分布）,并且假設兩者互不相關。

貝葉斯理論是集合濾波資料同化方法的理論基礎。基于貝葉斯理論的貝葉斯濾波器從原理、方法和符號系統(tǒng)上為序貫資料同化方法提供了一個統(tǒng)一的框架［48—49］。貝葉斯濾波器以概率分布密度函數(shù)為對象,其目標是:在初始狀態(tài)場的概率分布密度函數(shù)p（x0）為已知的前提下,利用模式算子式（1）和觀測算子式（2）遞歸地更新狀態(tài)場的后驗概率分布密度函數(shù)p（xk|y1:k）,其中y1:k代表從時刻t1到tk的所有觀測,xk| y1:k意味著利用所有這些觀測信息對xk作出的估計。貝葉斯濾波器的基本前提有兩個:第一,模式狀態(tài)服從一階M arkov過程p（xk| x0:k- 1）= p（xk| xk- 1）,也就是說每一步預報的結(jié)果只依賴于預報初值,與歷史狀態(tài)場不相關；第二,觀測值與模式當前狀態(tài)場相互獨立。

遞歸貝葉斯濾波器的每個濾波循環(huán)分為兩步——預報和更新。其中的預報過程利用前一時刻狀態(tài)場的后驗概率分布密度函數(shù)p（xk- 1|y1:k- 1）得到當前時刻的先驗概率分布密度函數(shù)p（xk| y1:k- 1）。假設p（xk- 1| y1:k- 1）已知,根據(jù)M arkov過程轉(zhuǎn)移密度的Chap man-Kolmogorov方程,下一時刻狀態(tài)場的先驗概率分布密度為［50］:

式中,p（xk| xk- 1）稱為轉(zhuǎn)移概率密度函數(shù),等價于模式誤差的概率表達。當模式誤差取成式（1）中的加法噪聲形式時,可以簡單地認為:

這也是η的概率分布密度函數(shù)。先驗概率分布密度函數(shù)刻畫了我們僅根據(jù)模式而不利用觀測信息得到的知識,其分布密度函數(shù)僅與模式初值和模式誤差的分布相關。而更新過程（在同化術語中一般稱為分析）通過整合觀測信息來更新狀態(tài)場,得到后驗概率分布密度。根據(jù)貝葉斯定理,系統(tǒng)狀態(tài)的后驗概率分布密度函數(shù)為:

公式（3）和（5）完整地表達了序貫資料同化方法的貝葉斯遞歸濾波形式:即模式系統(tǒng)的動態(tài)演進和觀測信息的不斷融合,在考慮模型和觀測誤差的基礎上,不斷地更新對系統(tǒng)狀態(tài)的估計。貝葉斯濾波從原理上為各種資料同化方法勾勒了一個統(tǒng)一的框架,集合卡爾曼濾波器和粒子濾波器的具體算法都可以在此基礎上得到。

3.1 集合卡爾曼濾波器

如果假設式（5）中先驗概率分布密度函數(shù)、似然函數(shù)和后驗概率分布密度函數(shù)為協(xié)方差矩陣分別為Pf、R和Pa的多維正態(tài)分布,展開式（5）并對公式兩邊取對數(shù),再使得Pa最大化,可以得到如下的目標函數(shù)

這與從線性無偏估計得到的三維變分的目標函數(shù)相同［5］。為了方便表示,式（7）中使用xfk代替先驗的狀態(tài)場xk| y1:k- 1。由式（7）很容易推導出經(jīng)典卡爾曼濾波器的分析公式:

式中,H為h的線性化算子（公式（7）～（9）的具體推導見附錄1）。

集合卡爾曼濾波器中采用了集合的概念,預報集合中的每一個成員實際上是通過對不同的初值進行積分得到的。使用預報場的集合成員可以計算其經(jīng)驗協(xié)方差矩陣,即

3.2 粒子濾波器

與集合卡爾曼濾波器中采用的高斯分布假設不同,粒子濾波器使用一系列粒子的加權(quán)和來表示各階段的概率分布,并且直接使用貝葉斯濾波式（3）和（5）來進行預報和分析。假設前一時刻的后驗概率分布密度函數(shù)可以表示為如下形式:

這是一個以xk為變量的連續(xù)概率分布密度函數(shù),由,i = 1,2,…,N｝這一組基構(gòu)成。接下來使用式（5）進行分析,把式（16）代入,就可以得到后驗概率分布密度函數(shù)如下:

這依然只是一個連續(xù)的分布密度函數(shù),我們需要對其取樣來得到與式（15）對應的δ函數(shù)的有限組合形式。理論上,我們需要根據(jù)一個建議分布q（xk）來取樣得到新的粒子,也就是說令-）,那么式（17）可以改寫為:

假定每一個貝葉斯濾波循環(huán)中的權(quán)重之和都為1,那么因子A可以看作是一個歸一化因子。因此,權(quán)重的遞歸關系可以表示為:該式中的正比符號意味著需要一個歸一化過程使得等式的兩邊相等。

也就是說,僅使用似然函數(shù)的值就能夠更新每個粒子的權(quán)重。

基礎的粒子濾波器算法不需要使用類似集合卡爾曼濾波器中的預報場（xf）和分析場（xa）這兩套不同的符號來區(qū)別分析前后的粒子。這是因為先驗的粒子（xk| y1:k- 1）和后驗的粒子（xk| y1:k）在數(shù)值上并沒有發(fā)生變化,變化的僅僅是兩者的權(quán)重。這一點也是退化問題發(fā)生的根源。為了處理退化問題, Gordon等［38］引入了重取樣過程,將其應用于式（20）對應的權(quán)重更新之后的加權(quán)集合。借用Van Leeuwen［36］綜述論文中的示意圖（圖1）,可以很清楚地看出重取樣過程的實質(zhì)就是根據(jù)權(quán)重的大小重新分配樣本。圖中黑點的大小代表了粒子的權(quán)重, 圖1a由于沒有采用重取樣過程,很快就會發(fā)生所有權(quán)重集中在同一個粒子上的退化現(xiàn)象；而圖1b在求得權(quán)重之后,將權(quán)重較大的粒子復制多份,同時消除權(quán)重較小的粒子,使用這樣的粒子集進行預報和更新就不容易產(chǎn)生退化。重取樣算法有很多種,但其基本原則大致相同［40］,加入重取樣的基本粒子濾波器也常常被叫做重取樣粒子濾波器。

3.3 集合卡爾曼濾波器和粒子濾波器的關聯(lián)和比較

如前面所述,貝葉斯濾波理論提供了一個統(tǒng)一的框架,各種資料同化方法都可以由此推出。粒子濾波器通過蒙特卡洛的思想直接應用貝葉斯濾波理論,當粒子數(shù)目趨向于無窮的時候,它可以收斂到貝葉斯濾波器的結(jié)果［54］。而集合卡爾曼濾波器則不同,它采用了高斯近似來逼近貝葉斯濾波器的結(jié)果。也就是說,集合卡爾曼濾波器衍生的方法預先假定了貝葉斯式（3）和（5）中涉及到的所有概率分布密度都是高斯分布的。理論研究證明:在一個分析循環(huán)中,只要模式是非線性的,即使初始場、模式噪聲和觀測噪聲都服從高斯分布,當集合尺寸N趨于無窮時,預報集合和分析集合的分布密度函數(shù)也都會趨向于錯誤的極限［55］。也就是說兩個集合的極限分布會不同于貝葉斯預報場式（3）和貝葉斯分析場式（5）的分布。關于這一點有非常嚴密的數(shù)學證明,但是最直觀的解釋還是因為模式的非線性會不可避免地帶來分布的非高斯性。所以高斯分布的初始集合加上高斯噪聲會產(chǎn)生非高斯的預報集合,然后使用預報集合的高斯近似進行分析,也會導致后驗分布的偏差。該理論指出的是集合卡爾曼濾波器類的方法在同化效果上的極限,即在集合成員數(shù)目足夠多時,集合卡爾曼濾波器的同化效果會逐漸趨于一個不正確的穩(wěn)定結(jié)果。

圖1 粒子濾波器中的粒子退化與重取樣示意圖［36］Eig .1 The schematic diagram for filter degeneracy and resampling in SIR-PE［36］黑點代表粒子,點的大小代表權(quán)重大小Black dot represents the particle,the size of dot represents the weight value

目前在集合資料同化方法的研究中,集合卡爾曼濾波器相較于粒子濾波器得到了更多的關注和應用,其最主要的原因是由于集合卡爾曼濾波器及一系列的衍生方法所需要的計算成本較小。這里的計算成本主要體現(xiàn)在集合成員數(shù)量上。以前面介紹的集合卡爾曼濾波器和重取樣粒子濾波器為例,一般來說,即使對于狀態(tài)空間場維數(shù)巨大的數(shù)值模式,集合卡爾曼濾波器使用10～100個左右的集合成員就能保證較好的同化效果［56］。而相比較而言,重取樣粒子濾波器需要非常大的集合來防止發(fā)生退化。更糟的是,在重取樣粒子濾波器中這個數(shù)目還會隨著問題的維數(shù)而增長。Bocquet等［33］和Nakano等［44］使用了洛倫茨63（L63）模式［57］和洛倫茨96（L96）模式［58］對這兩種方法進行了一系列的比較。結(jié)果表明,對于只有3維的L63模式,重取樣粒子濾波器需要200個左右的集合成員來使得同化結(jié)果的均方根誤差（R M SE）明顯小于集合卡爾曼濾波器；而對于10維和40維的L96模式,則分別需要10 000和30 000以上的集合成員數(shù)才能確保重取樣粒子濾波器的優(yōu)勢。

兩者比較的結(jié)論:集合卡爾曼濾波器非常適用于狀態(tài)空間維數(shù)巨大的動力系統(tǒng),即使集合成員的數(shù)目比較小,該方法也能夠給出后驗分布的前兩階統(tǒng)計矩（均值和協(xié)方差）的合理估計。再加上集合卡爾曼濾波器可以采納局地化方案以抑制長距離的虛假相關,進一步地避免了發(fā)生類似粒子濾波器中的退化。與之相比,重取樣粒子濾波器更適合于維數(shù)較小的強非線性系統(tǒng),它能夠突破高斯假設的限制,達到更好的同化效果。一旦系統(tǒng)的維數(shù)增加,“維數(shù)災難”問題將直接導致重取樣粒子濾波器的計算成本無法負擔。因此盡管有不少文章對重取樣粒子濾波器與集合卡爾曼濾波器進行比較,但是在實際高維海洋模式的資料同化中,重取樣粒子濾波器是無法使用的［45］。如何在高維的模式系統(tǒng)中使用粒子濾波器仍然是一個極具挑戰(zhàn)性的問題,是當前集合資料同化方法研究的熱門方向［29,31］。目前的研究表明,在粒子濾波器算法中采用觀測依賴的建議分布,能夠極大地降低權(quán)重退化的幾率,降低運算成本,這也是當前粒子濾波器研究的熱門方向。

4 粒子濾波器算法的新進展和研究重點

由3.2可知,所有的粒子濾波器都基于使用Dirac-Delta函數(shù)的加權(quán)和形式來表示完整的概率分布密度函數(shù),但是權(quán)重的更新公式根據(jù)不同的算法又各有不同。這是因為更新公式（18）中涉及的建議分布密度函數(shù)q（xk）可以選擇不同的形式。為理解建議分布這個概念,假設兩個概率分布密度函數(shù)分別記為p（x）和q（x）,根據(jù)q（x）隨機取樣得到的樣本集合｛xi,i = 1,…,N｝顯然逼近于q（x）的分布。進一步地,假如定義權(quán)重為wi= p（xi）/q（xi）,那么這些加權(quán)樣本｛xi,wi｝就會逼近于p（x）的分布。這兩個分布密度函數(shù)p（x）和q（x）在統(tǒng)計上分別被稱為目標分布和建議分布。

重取樣粒子濾波器以轉(zhuǎn)移概率密度函數(shù)p（xk| xk- 1）為建議分布,推出形式上比較簡單的遞推公式（19）?？紤]到轉(zhuǎn)移密度函數(shù)等價于模式噪聲的概率表達,可以看到重取樣粒子濾波器實際上不對預報粒子進行任何的處理,而是直接使用以上一步的分析粒子為初值得到的預報粒子集合去進行分析。當狀態(tài)場的空間維數(shù)很大時,因為所有粒子的權(quán)重都太小而較容易產(chǎn)生退化。最近的研究結(jié)果都指出,使用依賴于觀測的建議分布q（xk,yk）可以在進行分析之前對預報粒子進行一定的處理,使其更接近觀測資料,減少權(quán)重退化的幾率。這一類的研究是當前粒子濾波器方法的研究熱點。

最理想的情況是令q（xk,yk）= p（xk| xk- 1,yk）,這個分布在相關文獻中被稱為最優(yōu)建議分布,它可以使得tk時間新產(chǎn)生的粒子集合｛xik｝對應各權(quán)重｛wi｝的方差最小。由貝葉斯定理可得到:

將式（22）帶入式（18）就得到:

對比式（19）可以看出,由式（21）更新的wi與從建議分布中隨機選取的樣本xik并不相關,因此不會增加權(quán)重集合的方差。在實際同化中,考慮到對觀測的依賴性,這個最優(yōu)建議分布的表達式很難通過計算得到。因此,目前大量工作的重點都在于逼近最優(yōu)分布或者在特定的情況下將其簡化。

M orzfeld等［59］利用一種依賴于有效平均值的方法求解得到了p（xk| xk- 1,yk）的解析表達式。該表達式適用于模式和觀測算子都含有加法高斯噪聲,且觀測算子為線性函數(shù)（記為H）的簡單情況。Snyder給出了具體的解析式［47］如下:假設同化系統(tǒng)為

式中,ηk～N（0,Q）,ζk～N（0,R）。這里的～p（x）代表樣本是依據(jù)概率分布密度函數(shù)p（x）隨機取樣得到的,N（0,Q）代表均值為0,協(xié)方差為Q的（多元）高斯分布。根據(jù)Doucet等［60］的理論,由于模式噪聲η服從高斯分布可得xk| xk- 1～N（f（xk- 1）, Q）,利用標準的卡爾曼濾波器公式,可得xk| xk- 1,yk～N（,P）,其中

也就是說,最優(yōu)建議分布p（xk| xk- 1,yk）可以寫成高斯分布N（,P）的形式。而相應的權(quán)重更新公式（21）也可以利用以下關系得到:

上面提到的是比較理想的情況,在實踐中鑒于計算量的考慮,一般使用集合成員來近似表示這些高斯分布。Papadakis等［61］通過比較總結(jié)了集合卡爾曼濾波器和重取樣粒子濾波器兩種方法的特點,并在此基礎上提出了加權(quán)集合卡爾曼濾波器（weighted ensemble Kalman filter）以結(jié)合兩者優(yōu)勢。雖然該方法的名字中含有“卡爾曼濾波器”,但它實際上是一種粒子濾波器。其基本原理相當于使用集合卡爾曼濾波器的后驗概率分布（由前面的論述很容易知道它服從高斯分布）作為最優(yōu)建議分布的一個近似,在此基礎上求得權(quán)重并用重取樣方法更新集合（加權(quán)集合卡爾曼濾波器的具體算法參考附錄2）。可以看出,加權(quán)集合卡爾曼濾波器在算法的實行上非常簡便,可以看作是依次執(zhí)行集合卡爾曼濾波器和重取樣粒子濾波器,只是在權(quán)重的更新公式上進行了一定的修改。這也是這種粒子濾波器會被叫做“加權(quán)集合卡爾曼濾波器”的主要原因。加權(quán)集合卡爾曼濾波器自提出以來得到了很多的關注和應用,具有進一步研究和應用的價值［62—63］。

同樣地,為了結(jié)合集合卡爾曼濾波器和重取樣粒子濾波器這兩類方法的特點,Erei和Kunsch提出了集合卡爾曼粒子濾波器（ensemble Kalman particle filter）［64］。與加權(quán)集合卡爾曼濾波器的思想（使用集合卡爾曼濾波器的分析作為建議分布）不同,集合卡爾曼粒子濾波器通過使用一個可調(diào)制的參數(shù)控制集合卡爾曼濾波器和重取樣粒子濾波器兩階段的比重,構(gòu)建了一個混合高斯模型（Gaussian Mixture M odel）,以達成可計算性和同化非高斯信息能力的平衡。我們最新的工作進一步發(fā)展了這種方法,使之能夠適用于非線性的觀測算子［65］。數(shù)值實驗表明,集合卡爾曼粒子濾波器能非常有效地改善集合卡爾曼濾波器在強非線性模式系統(tǒng)中的同化效果。但是因為自適應算法調(diào)制參數(shù)的計算量比較大,我們還需要一些更深入的研究才能將其應用到較高維的系統(tǒng)。

除了上述的方法以外,其他一些最新的粒子濾波器也大多集中于使用這種依賴于觀測的建議分布。例如,Van Leeuwen［66］通過充分改進建議分布提出了一種完全非線性的粒子濾波器,他們稱之為等權(quán)重粒子濾波器（Equivalent W eights Particle Eilter）。等權(quán)重粒子濾波器的基本思想是在模式運行中對每一個粒子求解一個目標為“使得各粒子的權(quán)重大致相等”的最優(yōu)化問題,從而避免退化的發(fā)生。這種方法也受到了大量的關注［67 - 68］。此外,Chorin等提出的隱式粒子濾波器（Implicit Particle Eilter）［59,69］,Luo等提出的帶殘量拖曳的粒子濾波器（Particle Eilter with Residual Nudging）［70］,這些方法也在一定程度上促進了粒子濾波器的發(fā)展。

5 總結(jié)和展望

隨著各種觀測數(shù)據(jù)的迅猛增多,以及大氣海洋數(shù)值模式的不斷發(fā)展,如何充分利用各種觀測數(shù)據(jù)以滿足數(shù)值模式和預報的需要是一個非常重要的課題。作為連接觀測和模式的橋梁,資料同化技術日益受到關注,并得到深入的研究與廣泛的應用。

本文著重探討了集合資料同化方法的一些研究背景和發(fā)展現(xiàn)狀,主要包括兩大類的方法:集合卡爾曼濾波器和粒子濾波器。集合卡爾曼濾波器是為數(shù)不多的由海洋學者率先提出的資料同化方法［7,12］,在海洋資料同化之中已經(jīng)有了非常廣泛的應用,其同化效果也得到了普遍認可。與另一種在業(yè)務化中普遍使用的變分方法相比,集合卡爾曼濾波器的主要優(yōu)勢如下:（1）使用流依賴的背景場協(xié)方差矩陣；（2）無需使用切線性模式和伴隨模式；（3）容易并行化使用。國內(nèi)外已經(jīng)有了非常多的關于集合卡爾曼濾波器方法本身的研究,探討了初始集合如何產(chǎn)生、如何使用協(xié)方差松弛技術、如何進行局地化以及如何進行參數(shù)估計等方面問題。在應用方面,集合卡爾曼濾波器也被廣泛用于諸如太平洋海洋高度計資料同化［71］、E NSO的預報［72］和印度洋海面高度異常的短期預報［73］等方面。預計在不久的將來,還會在參數(shù)估計以及耦合模式的資料同化方面得到更多的應用。然而,集合卡爾曼濾波器也有一些不足:比如需要一定的集合樣本數(shù)目來避免背景場協(xié)方差的低估；積分模式的計算成本較高等。從原理上看,還可以知道集合卡爾曼濾波器及其衍生的方法都隱含了高斯分布假設——即假定預報誤差、模式誤差和觀測誤差都服從高斯分布。在該假設下,更新公式只使用前兩階的統(tǒng)計矩（均值和協(xié)方差）來刻畫集合的概率分布密度函數(shù)。當實際預報集合不服從高斯分布時,同化效果會有一定的限制。

作為另一類的集合資料同化方法,粒子濾波器被更多地用于機械、計算機和信號處理等工程領域,只是在近幾年才逐漸被認識到可用于地球系統(tǒng)的資料同化［16］,并逐漸成為了研究的熱門方向之一。從理論上看,包括集合卡爾曼濾波器和粒子濾波器在內(nèi)的序貫同化方法都可以從貝葉斯濾波的框架推導得到。相應于集合卡爾曼濾波器的高斯分布假設,粒子濾波器對于預報集合的分布不作任何的先驗假定,而是采用大量的樣本去估計完整的概率分布密度函數(shù)。在貝葉斯框架內(nèi),我們非常容易理解這兩類方法的強項和缺陷:集合卡爾曼濾波器的分析場分布隨著N的增大能夠很快收斂,但是無法收斂到由貝葉斯濾波器給出的真實的極限；粒子濾波器能夠收斂到真實的極限,但是需要的集合成員數(shù)量N過于巨大,難以在實際同化中實現(xiàn)。一些簡單模式中的比較［33］也揭示了這一點。

如何在高維的模式系統(tǒng)中使用粒子濾波器一直是一個挑戰(zhàn)性的問題,還沒有得到徹底的解決,也是當前研究的一個熱門方向。理論已經(jīng)證明了基本的重取樣粒子濾波器無力解決高維系統(tǒng)中的濾波退化問題（“維數(shù)災難”）,需要更多的努力來避免粒子的效率過低。通過對一系列最新研究的總結(jié),本文指出解決粒子濾波器計算量大的一個主要手段是使用依賴于觀測的建議分布。理論上的最優(yōu)建議分布難以在實際同化中計算得到,所以大部分當前可行的粒子濾波器都采用它的近似來作為產(chǎn)生新樣本的建議分布。在當前比較受關注的新興粒子濾波器中,本文選擇加權(quán)集合卡爾曼濾波器做了較詳細的介紹。這不僅僅是因為加權(quán)集合卡爾曼濾波器的算法表述比較簡明,更重要的是,它在集合卡爾曼濾波器和粒子濾波器之間構(gòu)筑了一個橋梁,非常好地展現(xiàn)了兩者的區(qū)別和關聯(lián)。此外,本文也提到了一些當前熱門的粒子濾波器。

總的來說,集合資料同化方法目前正處于快速發(fā)展中。隨著當前數(shù)值模式復雜性的不斷提高,計算能力的不斷增長,集合卡爾曼濾波器和粒子濾波器的重要性正在逐漸凸顯。當然,兩類資料同化方法的研究中都還有非常多急需解決的難題,本文通過對一些理論工作的整理和總結(jié),旨在為將來集合資料同化方法的發(fā)展提供一些理論基礎。

參考文獻:

［1］ Yan C,Zhu J,Xie J . An ocean reanalysis system for thejoining area of Asia and Indian-Pacific ocean［J］. Atmospheric and Oceanic Science Letters, 2010,3（2）:81 - 86 .

［2］ Han G,Li W ,Zhang X,et al. A regional ocean reanalysis system for coastal waters of China and adjacent seas［J］. Advances in Atmospheric Sciences,2011,28（3）:682 - 690 .

［3］ 馬建文,秦思嫻.數(shù)據(jù)同化算法研究現(xiàn)狀綜述［J］.地球科學進展,2012,27（7）:747 - 757 . M a Jianwen,Qin Sixian . Recent advances and development of data assimilation algorith ms［J］. Advancesin Earth Science,2012,27（7）:747 - 757 .

［4］ 李宏,許建平.資料同化技術的發(fā)展及其在海洋科學中的應用［J］.海洋通報,2011,30（4）:463 - 472 . Li H ong,Xu Jianpin . Development of data assimilation and its application in ocean science［J］. M arine Science Bulletin,2011,30（4）:463 - 472 .

［5］ Courtier P,Andersson E,Heckley W ,et al. The EC M W E implementation of three-dimensional variational assimilation（3 D-Var）.I:Eormulation ［J］. Quarterly Journal of the Royal M eteorological Society,1998,124（550）:1783 - 1807 .

［6］ Le Dimet E X,Talagrand O . Variational algorith ms for analysis and assimilation of meteorological observations:theoretical aspects［J］. Tellus A, 1986,38（2）:97 - 110 .

［7］ Evensen G .Sequential data assimilation with a nonlinear quasi-geostrophic model using M onte Carlo methods to forecast error statistics［J］.Journal of Geophysical Research:Oceans（1978 - 2012）,1994,99（C5）:10143 - 10162 .

［8］ Gelb A . Applied optimal estimation［M］. Boston:The MIT Press,1974 .

［9］ Evensen G . The ensemble Kalman filter:Theoreticalformulation and practicalimplementation［J］. Ocean Dynamics,2003,53（4）:343 - 367 .

［10］ Kalnay E . Atmospheric modeling,data assimilation,and predictability［M］. Cambridge:Cambridge U niversity Press,2003 .

［11］ Park S K,Xu L . Data Assimilation for Atmospheric,Oceanic and H ydrologic Applications（Vol.Ⅱ）［M］. Berlin:Springer Science & Business M edia,2013 .

［12］ Burgers G,Van Leeuwen P J,Evensen G . Analysis scheme in the ensemble Kalman filter［J］. M onthly W eather Review,1998,126（6）:1719 -1724 .

［13］ 劉成思,薛紀善.關于集合Kalman濾波的理論和方法的發(fā)展［J］.熱帶氣象學報,2005,21（6）:628 - 633 . Liu Chengsi,Xue Jishan . The ensemble Kalman filter theory and method development［J］.Journal of Tropical M eteorology,2005,21（6）:628 -633 .

［14］ Gillijns S,M endoza O B,Chandrasekar J,et al. W hatis the ensemble Kalman filter and how well doesit work？［C］//Proceedings of the A merican Control Conference,2006 . Minneapolis,Minnesota,U SA,2006.

［15］ Bengtsson T,Snyder C,Nychka D . Toward a nonlinear ensemble filter for high-dimensional systems［J］.Journal of Geophysical Research,2003, 108（D24）:8775 .

［16］ M oradkhani H,Hsu K L,Gupta H,et al. U ncertainty assessment of hydrologic modelstates and parameters:Sequential data assimilation using the particle filter［J］. W ater Resources Research,2005,41（5）:W 05012 .

［17］ 畢海蕓,馬建文.粒子濾波算法在數(shù)據(jù)同化中的應用研究進展［J］.遙感技術與應用,2014,29（5）:701 - 710 . Bi Haiyun,M a Jianwen . Advancesin the study of partcilefilterin data assimilation［J］. Remote Sensing Technology and Application,2014,29（5）: 701 - 710 .

［18］ Kalman R E . A new approach to linear filtering and prediction problems［J］.Journal of basic Engineering,1960,82（1）:35 - 45 .

［19］ Bishop C H,Etherton B J,M aju mdar S J . Adaptive sampling with the ensemble transform Kalman filter .PartⅠ:Theoretical aspects［J］. M onthly W eather Review,2001,129（3）:420 - 436 .

［20］ Anderson J L . An ensemble adjustment Kalman filter for data assimilation［J］. M onthly W eather Review,2001,129（12）:2884 - 2903 .

［21］ Tippett M K,Anderson J L,Bishop C H,et al. Ensemble square root filters［J］. M onthly W eather Review,2003,131（7）:1485 - 1490 .

［22］ Li H,Kalnay E,Miyoshi T .Simultaneous estimation of covarianceinflation and observation errors within an ensemble Kalman filter［J］. Quarterly Journal of the Royal M eteorological Society,2009,135（639）:523 - 533 .

［23］ Miyoshi T . The Gaussian approach to adaptive covariance inflation and its implementation with the local ensemble transform Kalman filter［J］. M onthly W eather Review,2011,139（5）:1519 - 1535 .

［24］ Zheng X,W u G,Zhang S,et al. Using analysis state to construct a forecast error covariance matrix in ensemble Kalman filter assimilation［J］. Advances in Atmospheric Sciences,2013,30（5）:1303 - 1312 .

［25］ Shu Y,Zhu J,W ang D,et al. Assimilating remote sensing and in situ observationsinto a coastal model of northern South China Sea using ensemble Kalman filter［J］. Continental Shelf Research,2011,31（6）:S24 - S36 .

［26］ H unt B R,Kostelich E J,Szunyogh I . Efficient data assimilation for spatiotemporal chaos:A local ensemble transform Kalman filter［J］. Physica D:Nonlinear Phenomena,2007,230（1）:112 - 126 .

［27］ Elowerdew J . Towards a theory of optimallocalisation［J］. Tellus A,2015,67,25257 .

［28］ Perianez A,Reich H,Potthast R . Optimallocalization for ensemble Kalman filter systems［J］.気象集誌第2輯,2014,92（6）:585 - 597 .

［29］ H outekamer P L,Mitchell H L,Pellerin G,et al. Atmospheric data assimilation with an ensemble Kalman filter:Results with real observations ［J］. M onthly W eather Review,2005,133（3）:604 - 620 .

［30］ Buehner M ,Charette C,He B,et al.Intercomparison of 4-D Var and En K E systems for operational deterministic N W P［R］. W W RP/T H O RPE X W orkshop on 4 D-VA R and Ensemble Kalman Eilter Intercomparisons . Buenos Aires,2008 .

［31］ Samuelsen A,Bertino L,Hansen C .Impact of data assimilation of physical variables on the spring bloom from T O P A Z operational runs in the North Atlantic［J］. Ocean Science,2009,5（4）:635 - 647 .

［32］ Chassignet E P,H urlburt H E,M etzger E J,et al. U .S . G O D A E:Global Ocean Prediction With the H Ybrid Coordinate Ocean M odel（H Y C O M）［J］. Oceanography,2009,22（2）:64 - 75 .

［33］ Bocquet M ,Pires C A,W u L . Beyond gaussian statistical modeling in geophysical data assimilation［J］. M onthly W eather Review,2010,138（8）: 2997 - 3023 .

［34］ Han X,Li X . An evaluation of the nonlinear/non-Gaussian filters for the sequential data assimilation［J］. Remote Sensing of Environ ment,2008, 112（4）:1434 - 1449 .

［35］ Van Leeuwen P J . A variance-minimizing filter for large-scale applications［J］. M onthly W eather Review,2003,131（9）:2071 - 2084 .

［36］ Van Leeuwen P J . Particle filtering in geophysical systems［J］. M onthly W eather Review,2009,137（12）:4089 - 4114 .

［37］ Ham mersley J M ,Handscomb D C . M onte carlo methods［M］. London:M ethuen Publishing,1964 .

［38］ Gordon N J,Salmond D J,S mith A E . Novel approach to nonlinear/non-Gaussian Bayesian state estimation［J］. Proceedings of the IEEE,1993, 140（2）:107 - 113 .

［39］ Capp O,Godsill S J,M oulines E . An overview of existing methods and recent advances in sequential M onte Carlo［J］. Proceedings of the IEEE, 2007,95（5）:899 - 924 .

［40］ Douc R,Capp O . Comparison of resampling schemes for particle filtering［C］//Proceedings of the 4th International Symposiu m on Image and Signal Processing and Analysis,ISP A05 . Zagreb,Criatia,2005:64 - 69 .

［41］ Kotecha J H,Djuric P M . Gaussian su m particle filtering［J］.Signal Processing,IEEE Transactions on,2003,51（10）:2602 - 2612 .

［42］ Kotecha J H,Djuric P M . Gaussian particle filtering［J］.Signal Processing,IEEE Transactions on,2003,51（10）:2592 - 2601 .

［43］ Xiong X,Navon I M ,Uzunoglu B . A note on the particle filter with posterior Gaussian resampling［J］. Tellus A,2006,58（4）:456 - 460 .

［44］ Nakano S,Ueno G,Higuchi T . M erging particle filter for sequential data assimilation［J］. Nonlinear Processesin Geophysics,2007,14（4）:395 -408 .

［45］ Snyder C,Bengtsson T,Bickel P,et al. Obstacles to high-dimensional particle filtering［J］. M onthly W eather Review,2008,136（12）:4629 -4640 .

［46］ Dau m E,H uang J . Curse of dimensionality and particle filters［C］// Proceedings of the Aerospace Conference,IEEE . Newport,2003 .

［47］ Snyder C .Particlefilters,the“optimal”proposaland high-dimensionalsystems［C］//Proceedings ofthe EC M W E Seminar on Data Assimilation for Atmosphere and Ocean,2011 .

［48］ Diard J,Bessiere P,M azer E . A survey of probabilistic models using the bayesian program ming methodology as a unifying framework［C］//The Second International Conference on Computational Intelligence,Robotics and Autonomous Systems . Erance,2003 .

［49］ 李新,擺玉龍.順序數(shù)據(jù)同化的Bayes濾波框架［J］.地球科學進展,2010,25（5）:515 - 522 . Li Xin,Ba Yulong . A Bayesian filter framework for sequential data assimilation［J］. Advances in Earth Science,2010,25（5）:515 - 522 .

［50］ Chen Z . Bayesian filtering:Erom Kalman filters to particle filters,and beyond［J］.Statistics,2003,182（1）:1 - 69 .

［51］ H outekamer P L,Mitchell H L . A sequential ensemble Kalman filter for atmospheric data assimilation［J］. M onthly W eather Review,2001,129 （1）:123 - 137 .

［52］ Tang Y,A mbandan J,Chen D . Nonlinear measurement function in the ensemble Kalman filter［J］. Advances in Atmospheric Sciences,2014,31 （3）:551 - 558 .

［53］ Yang C,Min J,Tang Y .Evaluation oftwo modified Kalman gain algorith msfor radar data assimilation in the W R E model［J］. Tellus A,2015,67, 25950.

［54］ Crisan D,Doucet A . A survey of convergence results on particle filtering methods for practitioners［J］.Signal Processing,IEEE Transactions on, 2002,50（3）:736 - 746 .

［55］ Le Gland E,M onbet V,Tran V - D .Large sample asymptoticsforthe ensemble Kalman filter［M］//Crisan D . The Oxford Handbook of Nonlinear Eiltering . Oxford:Oxford U niversity Press,2009:598 - 634 .

［56］ Anderson J L . Ensemble Kalman filters for large geophysical applications［J］. Control Systems,IEEE,2009,29（3）:66 - 82 .

［57］ Lorenz E N . Deterministic nonperiodic flow［J］.Journal of the Atmospheric Sciences,1963,20（2）:130 - 141 .

［58］ Lorenz E N . Predictability:A problem partly solved［C］// Proceedings of the Seminar on Predictability . U K,1996 .

［59］ M orzfeld M ,Tu X,Atkins E,et al. A random map implementation ofimplicitfilters［J］.Journal of Computational Physics,2012,231（4）:2049 -2066 .

［60］ Doucet A,Godsill S,Andrieu C . On sequential M onte Carlo sampling methods for Bayesian filtering［J］. Statistics and Computing,2000,10（3）: 197 - 208 .

［61］ Papadakis N,M Min E,Cuzol A,et al. Data assimilation with the weighted ensemble Kalman filter［J］. Tellus A,2010,62（5）:673 - 697 .

［62］ Beyou S,Cuzol A,Gorthi S S,et al. W eighted ensemble transform kalman filter for image assimilation［J］. Tellus A,2013,65,18803 .

［63］ Nakano S Y . H ybrid algorith m of ensemble transform and importance sampling for assimilation of non-Gaussian observations［J］. Tellus A,2014, 66 .

［64］ Erei M ,Kunsch H R . Bridging the ensemble Kalman and particle filters［J］. Biometrika,2013,100（4）:781 - 800 .

［65］ Shen Z,Tang Y . A modified ensemble Kalman particle filter for non-Gaussian systems with nonlinear measurement functions［J］.Journal of Advances in M odeling Earth Systems,2015,7（1）:50 - 66 .

［66］ Van Leeuwen P J . Nonlinear data assimilation in geosciences:an extremely efficient particle filter［J］. Quarterly Journal ofthe Royal M eteorological Society,2010,136（653）:1991 - 1999 .

［67］ Ades M ,Van Leeuwen P J . An exploration of the equivalent weights particle filter［J］. Quarterly Journal of the Royal M eteorological Society, 2013,139（672）:820 - 840 .

［68］ Ades M ,Van Leeuwen P J . The equivalent-weights particle filterin a high-dimensional system［J］. Quarterly Journal of the Royal M eteorological Society,2015,141（687）:484 - 503 .

［69］ Chorin A J,M orzfeld M ,Tu X .Implicit particle filters for data assimilation［J］. Com municationsin Applied M athematics and Computational Science,2010,5（2）:221 - 240 .

［70］ Luo X,H oteit I . Efficient particle filtering through residual nudging［J］. Quarterly Journal of the Royal M eteorological Society,2014,140（679）: 557 - 572 .

［71］ 萬莉穎.集合同化方法在太平洋海洋高度計資料同化中的應用研究［D］.北京:中國科學院大氣物理研究所,2006 . W an Liying . Ensemble methods and applications to altimetry data assimilation in the Pacific［D］. Beijing:Institute of Atmospheric Physics,Chinese Academy of Science,2006 .

［72］ Deng Z,Tang Y . The retrospective prediction of E NSO from 1881 to 2000 by a hybrid coupled model:（Ⅱ）Interdecadal and decadal variationsin predictability［J］.Climate dynamics,2009,32（2/3）:415 - 428 .

［73］ Haugen V E,Evensen G . Assimilation of SL A and SST data into an O G C M for the Indian Ocean［J］. Ocean Dynamics,2002,52（3）:133 - 151 .

［74］ Arulampalam M S,M askell S,Gordon N,et al. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking［J］. Signal Processing,IEEE Transactions on,2002,50（2）:174 - 188 .

附錄1:從貝葉斯公式推導卡爾曼濾波器

在貝葉斯濾波理論中,狀態(tài)場的更新公式由貝葉斯定理給出:

式中,xk|y1:k- 1表示tk時間先驗的狀態(tài)變量（預報）,我們假設其滿足均值為xfk,協(xié)方差為Pf的多元高斯分布。同時,假設觀測誤差滿足均值為0,協(xié)方差為R的高斯分布,即

式中,d和m分別為狀態(tài)空間和觀測空間的維數(shù)。代入式（A1）,并在等式兩邊取對數(shù)得到

式中,A0是一個常數(shù)。

因此可以得到目標函數(shù)

通過最小化J（xk）可以使得后驗概率分布密度最大,在后驗概率分布密度函數(shù)為高斯分布的假定下這也相當于使后驗分布的協(xié)方差最小。容易發(fā)現(xiàn),式（A4）和從最佳線性無偏估計得到的三維變分的目標函數(shù)完全一致。

如果模式的觀測算子為線性的函數(shù)H ,不妨假設后驗狀態(tài)場（分析）是預報和觀測的線性組合

為簡化記號,我們這里忽略了時間下標k。根據(jù)預報誤差協(xié)方差和觀測誤差協(xié)方差的定義,可知

為使Pa最小,令

式中,trace代表矩陣的跡?？梢越獬?/p>

公式（A6）和（A7）是經(jīng)典的卡爾曼濾波器中的狀態(tài)場更新公式和協(xié)方差更新公式。而式（A9）則是卡爾曼增益矩陣的解析表達式。

附錄2:加權(quán)集合卡爾曼濾波器算法流程

（1）應用集合卡爾曼濾波器將觀測yk同化到預報集合中,得到分析集合并估算其協(xié)方差矩陣

其中

（3）應用標準的重取樣算法［74］:對于粒子每個,根據(jù)其權(quán)重大小計算重取樣指標s（i）,得到所有的重取樣指標之后,再令,i = 1,…,N。舉例來說,圖1中根據(jù)權(quán)重計算的重取樣指標為｛4,4,5,5,5,5,6, 6,6,7｝,意味著根據(jù)權(quán)重,前3個和最后3個粒子全部刪去,中間的第4到第7個粒子分別重復2、4、3、1份。具體求重取樣指標的多種算法可參考相關文獻［38］,重取樣后的粒子集合｛,i = 1,…,N｝再次變成等權(quán)重的。

然后就可以用這個等權(quán)重的粒子集合來積分模式,進入下一個“預報-分析”循環(huán)。

原論文中假設觀測算子H是線性的,實際上即使觀測算子是非線性的函數(shù)h ,我們只需要將加權(quán)集合卡爾曼濾波器按照第3.1節(jié)中集合卡爾曼濾波器的方法進行一定的擴展就可以拓展其適用范圍。

中圖分類號:P731

文獻標志碼:A

文章編號:0253-4193（2016）03-0001-14

收稿日期:2015-05-24；

修訂日期:2015-08-07。

基金項目:國家自然科學基金項目（41276029,41321004）；科技部國家基礎科研項目（2013CB430302）；衛(wèi)星海洋環(huán)境動力學國家重點實驗室自主課題（SO E DZZ1404,SO E DZZ1518）。

作者簡介:沈浙奇（1984—）,男,浙江省杭州市人,博士,從事數(shù)據(jù)同化方法的研究。E-mail:zqshen @ sio .org .cn

*通信作者:唐佑民,男,教授,研究員,主要從事大尺度海氣相互作用、可預報性理論、數(shù)據(jù)同化、E NSO和氣候預測等方面的研究。E-mail:ytang @ unbc .ca

沈浙奇,唐佑民,高艷秋.集合資料同化方法的理論框架及其在海洋資料同化的研究展望［J］.海洋學報,2016,38（3）:1 - 14,doi: 10.3969/j.issn .0253-4193.2016.03.001

Shen Zheqi,Tang Youmin,Gao Yanqiu . The theoreticalframework of the ensemble-based data assimilation method and its prospect in oceanic data assimilation［J］. Haiyang Xuebao,2016,38（3）:1 - 14,doi:10.3969/j.issn .0253-4193.2016.03.001

The theoretical framework of the ensemble-based data assimilation method and its prospectin oceanic data assimilation

Shen Zheqi1,Tang You min1,2,Gao Yanqiu1
（1 .State Key Laboratory of Satellite Ocean Environment Dynamics,Second Institute of Oceanography,State Oceanic Administration,Hangzhou 310012,China；2 .EnvironmentalScience and Engineering,University of Northern British Columbia,Prince George V2 N 4Z9,Canada）

Abstract:In the nu mericalsimulation ofthe ocean dynamic system,data assimilation is able to use thelimited observation data and nu merical modelto best estimate the ocean state,and effectively reduce the uncertainty from theinitial conditions . Therefore,data assimilation plays an important role in the study of modern physical oceanography . The ensemble Kalman filter（En K E）is an effective data assimilation method,which has attracted broad attention in oceanic data assimilation since it is proposed about twenty years ago .In recent years,the particle filter（PE）has become a hot research field,foritis not restricted by the linear and Gaussian assu mption of the model. This paper analyzes and su m marizes the current theories about the En K E and PE,in the framework of Bayesian filtering theory . The En K E and PE algorith ms are proposed and compared . On this basis,we further discuss the major obstacle for applying the particle filter in oceanic data assimilaiton at present. Some feasible solutions are also introduced . This paper is expected to provide theoretical basis for further develop ment and application of the ensemble-based data assimilation method in oceanic data assimilation .

Key words:data assimilation；ensemble Kalman filter；particle filter；non-Gaussian noise；Bayesian filter