局地誤差子空間變換卡爾曼濾波方法的最優參數選取

王昊運，王輝，張宇，萬莉穎

（1. 國家海洋環境預報中心，北京100081；2. 中國海洋大學海洋環境與生態教育部重點實驗室，山東青島266100；3. 國家海洋局海洋災害預報技術研究重點實驗室，北京100081）

1 引言

在大氣和海洋模式中，初始狀態對于數值模擬的效果至關重要。作為一種順序同化方法，集合卡爾曼濾波（Ensemble Kalman Filter，EnKF）及其演變的集合同化方法，集諸多優點于一身，是現有同化方法中最具應用發展前景的一個[1-3]。EnKF 通過將多個擾動的初始樣本作為一個集合[4]，利用這個集合估計背景誤差協方差，對卡爾曼濾波方法進行了有效簡化。EnKF 的優勢還在于：其計算代價比卡爾曼濾波和擴展卡爾曼濾波小得多；不要求背景誤差協方差是線性演變的；不要求發展模式的線性和伴隨模式；可以給集合預報提供好的初始擾動[4-5]。

為了解決EnKF 在大氣和海洋數值模擬的應用問題，最重要的兩個改進就是背景誤差協方差的膨脹和局地化分析[6]。通常，模式的狀態向量維數很高（為107），遠遠大于集合樣本的維數，這會導致背景誤差協方差矩陣中的虛假相關。同時因為集合樣本離散度的問題，也會導致同化分析對背景誤差協方差的低估[7-8]。EnKF及其演變而來的局地化分析在大氣和海洋數值模擬中得到了廣泛的應用[9-12]。

集合同化的優勢在于：集合同化方法不僅給出模式狀態的最優估計，而且不需要建立預報模式的切線性和伴隨，背景誤差協方差“流依賴”。在En-KF同化方法獲得成功的基礎上，為了不對觀測進行擾動，發展了一系列演變、改進的同化方法。例如：集合調整卡爾曼濾波（Ensemble Adjustment Kalman Filter，EAKF）、集合變換卡爾曼濾波（Ensemble Transform Kalman Filter，ETKF）、奇異演變插值卡爾曼濾波（Singular Evolutive Interpolated Kalman Filter，SEIK）、集合平方根濾波（Ensemble Square Root Filter，EnSRF）等。由于順序同化公式簡明、應用相對容易，EnKF 及其演變版本得到了快速的發展和應用[13-14]。觀測上不加擾動進行同化，能夠解決EnKF 計算量大和集合成員少時收斂速度慢的問題，這一系列方法被統稱為集合均方根濾波[15-17]。例如比EnKF 更加高效的SEIK 方法，由于采用了二階取樣法，因此可以用更少的集合樣本數達到比EnKF 更好的效果，同時避免了EnKF 對觀測向量擾動帶來的人為誤差，所以更加節省計算資源并且更加有效。實際應用中決定集合同化方法計算代價最關鍵的就是集合樣本數，它直接決定了該同化方法需要在預報步時積分模式的次數。Nerger等[18]結合了ETKF在低維集合樣本展開的子空間對誤差協方差矩陣估計的優勢和SEIK 采用二階取樣法可以帶來的小樣本數的優勢，提出了誤差子空間變換卡爾曼濾波（Error Subspace Transform Kalman Filter，ESTKF）。ESTKF 被證明是一種更加高效的誤差子空間濾波，其局地化分析方案稱為局地誤差子空間變換卡爾曼濾波（Local Error Subspace Transform Kalman Filter，LESTKF）[19]。

LESTKF 在實際應用中，同樣需要解決背景誤差協方差矩陣的低估和因為計算能力不足導致集合樣本數過小引起的虛假相關。因此局地化方案中的重要參數“局地化半徑”和LESTKF分析方案中用來膨脹背景誤差協方差矩陣的“ 遺忘因子”（Forgetting Factor）這兩個重要參數對同化效果起著決定性的作用[19]。研究這兩個參數如何影響LESTKF 的同化性能和同化效果，以及如何選取這兩個參數才能最小化分析誤差，對于將該同化方法應用于實際的大氣海洋模式是十分關鍵的[20]。本文利用Lorenz96 模型結合LESTKF 同化方法，通過“局地化半徑”和“遺忘因子”設置孿生實驗，研究這兩個重要參數對于同化效果的影響。

本文將首先探究“遺忘因子”對背景誤差協方差膨脹效果及同化結果的影響；其次將研究“局地化半徑”和“ 遺忘因子”的共同影響；最后結合Lorenz96 模型孿生實驗的結果分析這兩個重要參數在集合同化分析中扮演的角色，并針對如何優化選取參數提出建議和總結。

2 同化方法及同化框架

2.1 ESTKF同化方法介紹

對于非線性海洋大氣系統而言，海洋大氣在tk時刻的n 維狀態向量為xk，和它對應的誤差協方差矩陣為PK。m 個集合成員組成的狀態向量集合可表示為x(α)k,α = 1,…,m。用集合均值來表示對tk時刻的狀態估計：

同化預報時首先將集合樣本分別通過模式積分至同化時刻。

在ESTKF分析步中，背景誤差協方差Pf在形式上用集合樣本Xf表示，即：

式（7）和（8）中：T?是一個m ×( m - 1)維的滿秩矩陣，并且每列元素的和為零，即

T?矩陣的作用是在計算矩陣L 時剔除集合樣本矩陣Xf的集合平均，即計算集合擾動。值得注意的是，矩陣L 是個n ×( m - 1)維矩陣，只存儲前m - 1個集合擾動。

分析場通過集合擾動矩陣L給出：

式（10）——（12）中：權重向量wˉ、變換矩陣A?分別為（m-1）維、（m-1）×（m-1）維。參數取值范圍為0＜ρ?≤1，稱作“遺忘因子”，用來放大背景誤差協方差。ρ?作為該同化方法的重要參數，其取值直接影響著ESTKF的同化效果，也是本文中的主要研究對象之一。

該方法在得到分析場后需要進行再取樣。由于集合樣本矩陣經過了變換用X?a表示，P?a同理。

根據之前的研究，在SEIK 中，矩陣C?-1由A?-1經過Cholesky 分解得到，即( C?-1)TC?-1= A?-1。Ω 是一個m ×( m - 1) 維的矩陣，該矩陣所有列向量相互正交，并且與(1,…,1)T也正交。矩陣Ω 的作用是與( m - 1)×( m - 1) 維 的 集 合 變 換 矩 陣A?-1重 新 生成新的集合擾動矩陣。在ESTKF 中，Ω 矩陣記為Ω?，即：

式中：Ω?為Householder 矩陣[18]。Ω?的作用為將Xf（n×m 維）集合樣本空間中的向量投影在矩陣L（m-1個集合擾動）所在的誤差子空間。同樣矩陣Ω?是滿秩矩陣，并且列向量和為零。將式（8）——（13）中的T?矩陣替換為Ω?矩陣，最終算法就是ESTKF。

除了傳統的EnKF 等濾波的優勢之外，ESTKF的本質是在SEIK 的基礎上將集合變換矩陣的計算通過Ω?矩陣投影在誤差子空間來完成。該方法相對于SEIK 來說：計算矩陣L 時，不必像之前一樣忽略最后一列的集合樣本，同時剔除了集合成員順序的影響；ESTKF 由于在誤差子空間計算，所以計算代價要小于SEIK。

2.2 局地化方案

為了消除集合樣本數不足導致的背景誤差協方差矩陣中的虛假相關，假設只有距離模式格點在一定范圍（即局地化半徑）內的觀測才會對模式格點產生影響。局地化不僅可以減小虛假相關對同化效果的影響，還可以減小計算量，尤其當觀測數量遠遠大于模式集合樣本數時。而且局地化可以保證集合同化方法在高維模式應用上的合理性，即對于每個模式格點而言，模式狀態的調整是在一個相對較大的集合空間內實現的。需要注意的是，由于局地化半徑是為了消除虛假相關，而不同模式可能有不同的虛假相關，所以局地化半徑隨模式不同而不同。同時，局地化可能會在局地區域的邊緣引起模式狀態的不連續。為了解決這一問題，通常引入平滑的、以局地化半徑為參數、以觀測和格點距離為自變量的局地化函數，使得觀測的影響隨距離逐漸衰減。本研究中選用的局地化函數是五階Gaspari-Cohn 相關函數[17]，即：

式中：d 表示觀測到模式格點的距離，rloc表示局地化半徑，計算結果表示觀測對分析點的權重。利用Gaspari-Cohn 相關函數按照權重結合局地化半徑中的觀測得到該格點處的分析值。

參數ρ?的選取對同化的影響主要發生在LESTKF 同化的分析步中，在協方差矩陣的計算中引入了一個膨脹因子來增加濾波穩定性。因為模式的強非線性，導致在相空間中對初值十分敏感，在同化中太過相信模式預報的背景場將容易導致濾波發散，因此該參數的選取對于同化效果有著重要的影響。

局地化半徑通過把同化的區域分解成小的子區域，并行地同時更新每個格點的分析，更新時僅用到這個格點某一半徑內的所有觀測。在分析局地化的基礎上，對該半徑內的不同觀測引入權重的局地化方法被稱為觀測局地化（Observation Localization，OL）。研究表明，在局地化分析時，針對不同的集合樣本數局地化半徑存在一個最優的選擇[18]。

2.3 并行數據同化框架（Parallel Data Assimilation Framework，PDAF）

集合同化方法作為一種順序同化方法，在海洋和大氣領域得到廣泛應用。EnKF、集合最優插值（Ensemble Optimal Interpolation，EnOI）提出后，為了解決計算代價等問題，集合平方根濾波等分析方法演變出一系列的集合同化方法，如ETKF、EAKF、減秩卡爾曼濾波（Singular Evolutive Extended Kalman Filter，SEEK）、SEIK、ESTKF 等。這些集合同化方法都是先由預報步分別積分集合樣本；分析步僅需要模式提供的部分信息，通常只依賴于狀態向量，而不是單個物理量場。例如在海洋模式的狀態向量中，存儲著U、V、z、S、T等模式變量場，或是需要估計的模式參數。對于觀測算子H 的計算，只需要知道觀測位置在狀態向量中的存放位置即可。以上屬性使得建立一個用通用方式實現集合同化方法的核心算法，并通過調用通用接口來進行同化的同化框架成為可能。這將極大地降低同化方案的實施。由于通用接口的設定，同化參數也便于系統化的管理和調整。

PDAF 同化框架就是一些集合同化方法的算法庫（網址：http://pdaf.awi.de）[21]。目前全球大部分同化系統都是離線進行的，也就是模式集合積分與同化分析步分兩個程序進行。這種同化方式雖穩定但低效。因為模式積分程序和同化程序之間需要用文件來傳遞信息，這種運行方式稱為“離線模式”。另一種運行方式為直接將同化程序寫進模式代碼中整合成為一個程序，這種運行方式更高效，但是需要調整很小一部分模式代碼，這種運行方式稱為“在線模式”。前者是將模式積分和同化分析步分兩個單獨程序運行，模式的輸出結果輸入到同化程序，同化的分析場又作為模式的啟動場；后者是將同化程序與模式耦合起來，需要對模式代碼做擴展，調用PDAF 核心函數，形成一個完整的運行程序。 本研究中采用PDAF 同化 框 架結合Lorenz96 以“ 在線 模 式”運行。實際情況中由于模式復雜，往往采用“離線模式”[21-23]。

3 Lorenz96-PDAF孿生實驗

3.1 實驗設置及PDAF參數選取

Lorenz96 模型作為大氣和海洋的低階近似，具有非線性，存在混沌吸引子。Lorenz96 己經被廣泛應用于大氣海洋預測和資料同化的研究中[24]。其非線性強，同時對初值非常敏感，其動力框架為：

式中：Xi為模型變量,F為定常強迫，F的大小決定著非線性的程度。本文的實驗設計中，模型維數等于dim_state=40，即i = 1,2,…,40，強迫項F 選為8，數值積分采用四階龍格－庫塔格式。

對于孿生實驗而言，將模式積分10 000 步得到的狀態向量的時間序列作為“真值”，并用第1 001步時的模式狀態作為初始場（前1 000 步為模式調整階段，即“spin up”階段），其狀態向量序列即為40×10 000維的矩陣。

觀測序列是在Lorenz96 模式“真值”序列的基礎上，加上不相關的隨機噪聲生成的。本節中，為了探索不同同化強度下遺忘因子和局地化半徑的影響，擾動了不同的觀測序列，兩個觀測序列的標準差分別取為1.0和0.1。為了去除模式調整階段的影響，本實驗只同化1 000 步之后的觀測，采用積分一步同化一步觀測的方式。

用于集合同化的初始集合樣本由二階取樣法（Second-order Exact Sampling）在10 000 步模式“歷史真值”的基礎上生成初始樣本。本節為了探究遺忘因子和局地化半徑在不同樣本數下的同化效果，分別用10個、30個樣本做了對照試驗。集合孿生同化實驗流程見圖1。

圖1 集合孿生同化實驗示意圖

由于本孿生實驗的Lorenz96模型和同化方法均在PDAF同化框架中進行，表1給出了主要同化參數的設置和選項。其中dim_ens 為集合樣本數，forget和local_range分別對應“遺忘因子”和“局地化半徑”。

表1 同化實驗主要參數

3.2 遺忘因子對同化效果的影響

為了探究不同遺忘因子對同化結果的影響，本實驗采用local_range為5個模式格點、rms_obs=1.0、dim_ens=30；但將遺忘因子取不同值。將同化結果、模式預報結果與模式真值、估計值（30 個樣本的集合平均）做均方根誤差（Root Mean Square Error，RMSE）對比，綜合探究遺忘因子對同化結果的影響。為了方便表述，以下遺忘因子均記為f值。

圖2 給出了同化至1 000 步，即模式第2 000 步的真值、觀測值、估計值和分析場。圖2a 和2b 分別為第2 000 步的真值和觀測值；圖2c 和2d 分別為f =1.0 時模式估計場和同化分析場；圖2e和2f分別為f =0.9時的模式估計場和同化分析場。

通過對比發現模式第40 個格點上的觀測值相對真值出現了較大擾動（見圖2b、2d、2f 的紅框部分）。而只有f= 0.9的分析場有一個抬升，向觀測值做出了調整，這正是因為f 對背景誤差協方差的膨脹作用導致的分析結果向觀測偏移。

為了研究f 對分析場和預報場的影響，本實驗又將f 的取值擴大為0.5～1.0，并計算其同化5 000步后與真值的RMSE。

圖3 分別為預報場和分析場的RMSE，橫坐標為同化步（因為時間步過于密集，因此每隔50 步填值，橫坐標100 即代表第5 000 步，以此類推）。f 從上到下由0.5（協方差膨脹最大）遞增至1.0（協方差不膨脹）。

圖3 預報場和分析場RMSE（橫坐標為縮放50倍的時間步，縱坐標為與真值的RMSE）

對于單個RMSE 序列來說，在同化開始時RMSE 最大。這是因為初始樣本是從模式之前的“歷史狀態”中提取的，與真值的誤差最大；但隨著同化預報步和分析步的交替，預報場和分析場均迅速向真值收斂。但是隨著同化實驗的進行會出現誤差的突然增長。該特征與Pham[8]的圖1一致。這一現象可能由于Lorenz 系統進入了強的非線性區域所導致。

通過圖3a、b 的對比，發現LESTKF 對于所有實驗而言均對預報場有所改善，但不同f 值的改善程度不同。通過縱向比較可以看出：f為0.95與f為1.0時，f 的使用明顯避免了同化至3 000 步左右（圖中橫坐標約60）的誤差爆發；f 為0.5 時誤差最大；從上到下RMSE有一個先減小再增大的過程。

為了探究遺忘因子整體對同化結果的影響，分別對5 000 步同化實驗的預報場和分析場的RMSE求時間平均（見圖4）。結果表明：f =0.95 時分析場和預報場的誤差最小；隨著f 的減小，誤差隨之增大。f =1.0 時分析誤差和背景誤差均存在著低估，隨著f 減小低估得到改善；但隨著f 繼續減小，則出現了分析誤差和背景誤差的高估，導致同化結果過于接近觀測從而偏離真值。藍線和黑線的交點則是最優f 值。這說明遺忘因子通過控制對背景誤差協方差的估計會顯著影響同化效果，高估和低估均不是最優f值。

圖4a、b 的對比表明對于固定的f 而言，分析場相對于真值的RMSE 明顯小于預報場相對于真值的RMSE。例如當f =0.95 時，預報場相對于真值的RMSE 約為0.25，而分析場相對于真值的RMSE 約為0.2。這說明f 通過對背景誤差協方差的放大的確可以改善LESTKF 的同化效果；但隨著f 越來越小，分析場對于預報場的提升依然存在，但提升效果開始變小。

3.3 局地化半徑與遺忘因子的共同影響

本實驗根據Nerger 等[18]的實驗設置，將total_steps 設置為60 000 步，得到的結果作為真值；觀測依然在真值的基礎上增加隨機擾動，觀測誤差分別以1.0 和0.1 的標準差生成。前1 000 步仍然作為“spin up”階段。本實驗將通過不同的觀測誤差控制同化實驗的強度。集合樣本分別取10 個和30個。因為在實際的集合同化實驗中，往往因為計算資源的限制，不會選取過多的集合樣本，因此選取10個集合樣本更加接近實際情況。

local_range 的范圍為0～20 個格點，這是因為：Lorenz96 模式共有40 個格點，并且是周期邊界條件；local_range 超過20 個格點后相當于沒有局地化分析。將不同f 和local_range 組合進行同化實驗，結果依然用RMSE 的時間平均來表示。當時間平均RMSE 大于設定的觀測誤差時，認為發生了濾波發散。

圖5a 給出了30 個集合樣本、觀測誤差為1.0 的組合實驗結果；5b為10個集合樣本、觀測誤差為1.0的組合實驗結果；5c為10個集合樣本但觀測誤差降低至0.1 的組合實驗結果。格點上的數字代表二者組合實驗分析場的時間平均RMSE，白色格點說明該參數組合實驗結果的誤差超過了設定的擾動的觀測誤差，發生了濾波發散，因此不填色。

圖5 可以看出當集合樣本數為30 個時，不論local_range 和f 如何選取，均沒有發生濾波發散，而且存在一個最優的參數組合區域。這表明當集合樣本數足夠大、同化方法對背景誤差協方差的估計很精確時，引入遺忘因子和局地化分析依然對提高同化結果有幫助，能夠找到一個最優的參數搭配區域。例如圖5a中local_range＞6、f ＞0.93的區域。

當集合樣本數為10個時，這種情形更加接近實際同化，如圖5b、c 中出現了大面積的濾波發散區域。值得注意的是，最優的參數選擇區域和濾波發散區域非常接近，兩者緊鄰。例如local_range = 7、f = 0.98 時，RMSE 約為0.2；但是當local_rangen = 8、f仍然為0.98時，則出現了RMSE為3.5的濾波發散，這個結果與Nerger等[18]的結果仍然一致。Nerger解釋是由初始樣本隨機選取引起的同化實驗前期的過渡階段導致的，差的隨機初始樣本會使該過渡階段變長，因此導致濾波發散。

圖4 時間平均RMSE（藍線代表與真值的誤差，黑線代表與估計真值的誤差）

圖5 組合實驗結果（表中數字代表RMSE，白色區域為濾波發散區）

為了探討觀測誤差對于同化實驗的影響，將觀測重新用0.1 的標準差隨機生成，并同樣用10 個集合樣本做了組合實驗，結果見圖5c。當通過減小觀測誤差來增加同化強度時，濾波發散的區域有所增大，但依然存在一個最優的參數選擇區域。

4 總結與討論

本文通過Lorenz96 模型的孿生實驗，分別驗證了局地化半徑和遺忘因子的不同選取對同化分析場的影響，并且詳細介紹了PDAF 同化中一些參數的設置。結論如下：

（1）局地化半徑對分析結果的空間分布影響明顯：局地化半徑過大，不能很好地濾去背景誤差協方差矩陣中的虛假相關；局地化半徑過小則分析太細節化，使得物理量場不符合實際。

（2）遺忘因子作為單獨影響同化實驗的因子時，分析誤差隨著f 值的增大有一個先減小再增大的過程。這是因為f 值控制著對分析誤差協方差和背景誤差協方差的估計，對真值的低估和高估均會引起誤差。圖4中可以清楚的看到，估計的RMSE與實際的RMSE 的交點就是理論上最優的遺忘因子。但在實際同化應用中，由于對真值的估計不準確，因此對于參數f的選取只能通過實驗方法去尋找。

（3）f的選取對于同化效果影響顯著。Lorenz96理想實驗結果表明：f 的作用是人為放大背景誤差協方差，防止濾波發散。f 取值并非越小越好，太小會使同化結果過于接近模式，從而減弱觀測信息對模式的調整。f 恰當選?。ㄈ≈禐?～1），則可以明顯提高同化效果。因此在實際同化中選取該參數時應格外注意。

（4）遺忘因子和局地化半徑作為共同因子影響同化時，的確存在一個最優的區域。但是最優區域的選取需要慎重，因為最優的參數組合的區域往往在濾波發散的臨界區域附近。一旦選取不當，則很容易出現濾波發散，這在實際同化中一定要注意。本文為今后進一步利用LESTKF實際同化業務應用中參數的最優化選取做了鋪墊。

作為仍在不斷發展的同化方案，ESTKF 及PDAF 同化框架已經得到越來越廣泛的應用。由于PDAF 同化框架接口化、開源的特性，可以更方便地研究分析步中的其他參數并實現優化。例如，如何選取自適應的局地化半徑（隨緯度、同化的物理量而改變），以及如何進一步選取自適應的遺忘因子，以便協方差膨脹隨著同化問題不同而變化，是值得進一步研究的問題。

資料同化是一個綜合問題，不僅需要同化方案的演進，同時還需要計算技術的進步，二者缺一不可。怎樣在不提高計算代價的前提下提高同化效果，是資料同化理論研究中的一個重要方向。是否能夠提出一個更有效率的同化方案，并且能擁有集合濾波中隨時間演變的背景誤差協方差矩陣，仍是資料同化研究的一個重要問題。