集合變換卡爾曼濾波數據同化方法中的協方差局地化

王際朝， 王 玥， 臧紹東， 楊俊鋼， 紀艷菊， 阮宗利

(1. 中國石油大學(華東) 理學院， 山東 青島 266580； 2. 自然資源部第一海洋研究所， 山東 青島 266061)

數據同化是通過對觀測數據與預測值進行分析，尋找一個當前物理狀態的最優估計值(分析值)， 作為數值預測模型初始條件的方法。作為連接數值模式和觀測數據的技術手段， 數據同化在海洋科學領域得到了廣泛關注。當前， 集合數據同化方法得以快速發展， 已經成為業務化數值預測的一個可行選擇。集 合卡爾 曼 濾波(ensemble Kalman filter， EnKF)[1]，可以將卡爾曼濾波應用到強非線性模型， 擴展了卡爾曼濾波的適用性， 是數據同化領域的研究熱點之一。EnKF中的預測誤差協方差矩陣通過對預測集合成員進行統計計算得到， 預測誤差對于集合數據同化方法而言至關重要。然而， 產生足夠大的預測集合的計算成本令人望而卻步。當使用的集合數目過小時， 則不能在統計意義上充分表示模型的變化， 從而會引起欠采樣(under-sampling)。通常， 欠采樣會導致濾波發散、預測誤差協方差低估和虛假相關等負面問題。由于這些問題的影響， EnKF可能會產生一個次優的分析結果。

在進行數據同化時， 為了克服欠采樣帶來的負面影響， 有關學者已經提出了很多解決方法。當前基于集合的數據同化研究中， 最著名的解決方法是協方差膨脹(covariance inflation， CI)[2]、協方差局地化(covariance localization， CL)[3-4]和 局 地 化 分 析(local analysis， LA)[5-6]。CI方法通過一個膨脹因子校正預測誤差協方差的低估問題， 主要思想是通過膨脹集合均值和集合成員之間的誤差來增大預測誤差協方差。CI方法用于分析過程之前， 具體表示為：其中“←”表示對于之前狀態值的替代，r是所謂的膨脹因子。盡管CI方法克服了預測誤差協方差的低估問題， 但是協方差膨脹因子r解決不了虛假相關問題。CL方法在預測誤差協方差矩陣和局地化函數之間實施舒爾積運算[7]， 通過截斷預測誤差協方差矩陣中預先指定距離之外的誤差相關性以解決虛假相關問題。此外， 舒爾積運算能夠提高預測誤差協方差矩陣的秩[8]。LA方法以待更新狀態變量為中心， 建立一個虛擬的局部窗口， 得到該狀態變量預測誤差協方差的局部近似值， 在分析更新過程中， 將局部窗口外的集合擾動設定為零。例如， 通常可以以某個狀態變量為中心， 僅僅同化與其固定距離為M之內的觀測數據， 固定距離的長度決定了同化過程中觀測值的數量。顯然， 分別對狀態變量進行局部同化的特點使LA方法可以利用并行方法進行快速計算。

當前， 越來越多的學者[9-11]關注的是CL方法和LA方法。相比于CL方法， LA方法是一種獨立的方法， 它可以用于任何數據同化的框架， 但是它的實際同化效果并不如CL方法。Miyoshi等[12]指出LA方法的同化效果和CL方法類似， 但是LA方法通常會導致弱局地化影響。之后， Janji?等[9]通過對Lorenz96模型進行實驗發現： 如果觀測誤差方差小于初始狀態的預測誤差方差， CL方法將會導致比較小的估計誤差； 如果觀測誤差方差大于初始狀態的預測誤差方差， CL方法和LA方法有相同的局地化影響。因此， 從之前的研究和實驗可以得出， CL方法可能是相對較好的一個局地化方法。由于CL方法中的舒爾積運算僅僅用于預測誤差協方差矩陣，而集合轉換卡爾曼濾波(ensemble transform Kalman filter， ETKF)方法中的預測誤差協方差矩陣并沒有顯式計算， 而是傳遞預測集合擾動矩陣， 這導致ETKF方法中進行舒爾積運算的矩陣維度不一致，所以CL方法不能用于ETKF方法[13]。Janji?等[9]的文章給出了CL方法不適用于ETKF方法的具體討論。

為了克服CL方法不能用于ETKF方法的限制，在一元模型中， 通過對局地化函數的平方根矩陣和預測集合擾動矩陣進行舒爾積運算， Petrie[14]提出了一種用于ETKF方法的CL近似方法。對于矩陣維度不一致的問題， 通過在預測集合擾動矩陣中添加n-N列零向量(n表示狀態變量個數，N表示集合個數)， 對矩陣進行擴展來修正。然而， 韓培等[15]對此方法進行探究， 表明這種近似的CL方法是一種弱近似并產生了較差的同化效果。

基于集合樣本相關性和矩陣因式分解， 近幾年提出了動態計算局地化函數的方法， 使得CL方法可以應用到ETKF方法中。Bishop等[16]提出了一種生成局地化函數的新方法， 該函數能夠隨真實誤差相關函數移動并且還適應真實誤差相關函數的寬度。但是此方法的計算代價比較大， 鑒于此， Bishop等[17]提出了一種新的協方差自適應局地化方法(CALECO)的具體實施方案， 成功地將CL方法應用到ETKF方法， 但是此方法仍然具有較高的計算復雜度?；谏鲜龅难芯炕A， Bishop等[18]在ETKF方法中提出了較為可行的CL方法的實施方案， 并命名為增益集合轉換卡爾曼濾波(gain form of ensemble transform Kalman filter， GETKF)。GETKF中使用的是Gaspari等[19]提出的基于距離相關的一個局地裁剪函數作為局地化函數(GC局地化函數)， 此方法解決了在一元數值預測模型中CL方法不能用于ETKF方法的問題。

當前， 由于CL方法的實現存在限制， 國內外對于CL方法在ETKF等平方根濾波中的研究較少。本文將基于GETKF方法進行研究， 提出一種局地化效果較好的局地化方法， 并通過數值模擬實驗比較本文基于GETKF修改的方法與GETKF方法的同化效果。

1 理論背景與方法改進

Campbell等[20]比較了EnKF方法中分別使用模型空間CL方法和觀測空間CL方法的同化效果， 認為模型空間CL方法通常優于觀測空間CL方法。因此， 本文以下的研究基于模型空間CL方法。本節首先介紹GETKF方法， 然后對該方法進行改進。

1.1 增益集合轉換卡爾曼濾波

在ETKF方法中， 應用CL方法的難點在于預測誤差協方差矩陣并沒有顯式表示， 而是通過預測集合擾動矩陣隱式表達， 因而在預測集合擾動矩陣和局地化函數之間無法進行舒爾積運算， 當前主要的目的是在預測集合擾動矩陣和局地化函數的平方根矩陣中進行舒爾積運算來近似預測誤差協方差矩陣和局地化函數之間的舒爾積運算。假定預測集合擾動矩陣表示為

局地化函數ρ的平方根矩陣表示為W， 即

其中， 平方根矩陣W通過對局地化函數ρ進行特征值分解得到， 并且按特征值大小降序排列， 取前10個特征值和特征向量對構成， 即

一般情況， CL方法需要進行如下形式舒爾積運算，

由于之前討論的ETKF方法的局限性， GETKF方法中采用了集合擴展技術定義一個n× (N×L)維度的矩陣f

Z，

定義矩陣f Z表示局地化后的預測誤差協方差矩陣的平方根矩陣， 即

Bishop等[18]已經證明式(6)是成立的。因此成功地在預測集合擾動矩陣實現了CL方法?，F在使用代替作為當前的預測集合擾動矩陣， 即局地化后的預測誤差協方差矩陣的平方根矩陣， 重新構造集合擴展后的預測集合和預測誤差協方差矩陣。

令M=N×L， 則擴展預測集合表示為其中每個列向量如下定義：

其中，zk表示矩陣的第k列。很顯然擴展預測集合的均值與原始預測集合均值相同。擴展預測集合誤差協方差矩陣表示為：

GETKF方法的分析過程與標準的ETKF方法的分析過程類似， 區別是GETKF方法使用擴展預測集合以及擴展預測集合擾動矩陣代替原始預測集合和原始預測集合擾動矩陣定義標準化的預測-觀測集合擾動矩陣如公式(7)所示：

其中，

矩陣H表示觀測算子。為了提高ETKF方法同化的計算效率， 對標準化的預測-觀測集合擾動矩陣進行奇異值分解：

其中， 矩陣U和V是正交矩陣，的奇異值由矩陣∑的對角元素給出。接下來， 分析誤差協方差矩陣的平方根矩陣(分析集合擾動矩陣)表示為，

其中，

以上展示了GETKF方法的主要計算過程， 可以看出此方法存在改進的空間： 首先， GETKF方法利用集合擴展技術在分析過程中產生M個集合(N＜M)， 然后使用集合擴展之后的分析誤差協方差矩陣計算N個集合的方法較為復雜， 如何更加簡便地將M個集合轉換為N個集合進行預測過程是值得商榷的； 其次， GETKF方法中計算擴展后的集合誤差協方差矩陣采用有偏形式， 本文對此做出了修改，采用無偏形式的矩陣進行計算； 此外， GETKF方法中使用的GC局地化函數是一個基于距離相關的局地裁剪函數， 該函數是一元局地化函數， 所以擴展到多元模型存在限制； 最后， 由于GETKF方法將參數L固定為10， 所以在系統狀態變量個數較大的情況下， 擴展后的預測集合擾動矩陣不能很好地表示誤差變化?？梢钥闯鯣ETKF方法確實存在的不足之處， 對于GETKF方法的具體改進將在下一節進行說明。

1.2 增益集合轉換卡爾曼濾波的改進

基于GETKF方法， 本節從局地化函數平方根矩陣的選取、預測誤差協方差的無偏估計、隨機子采樣方法和GC局地化函數的限制4個方面對原方法進行改進， 以提高同化方法的效果。

1.2.1 局地化函數平方根的選取

對于公式(3)中局地化函數平方根矩陣列數L的選取， GETKF方法按照特征值降序排列， 取前十個特征值和對應的特征向量構成平方根矩陣W，即固定地取L為10。很顯然， 在應對狀態變量個數n很大時， 前10個特征值和特征向量不能完全表示系統變量的誤差變化， 所以本文重新定義了L的選?。?/p>

其中，E10%表示前10%的特征值的數量。即給定L的范圍， 通過進行多次數值實驗選取最優的數值L。

1.2.2 預測誤差協方差的無偏估計

通常， 我們假定卡爾曼濾波中的誤差協方差矩陣是無偏的。在統計學上， 通過對集合成員與集合均值的偏差除以集合自由度(樣本個數減一)來實現。但是GETKF方法中的公式(5)～(7)采用的是有偏形式的誤差協方差矩陣， 這導致誤差的產生， 造成對誤差協方差的低估， 所以本文通過將公式(5)～(7)中的來修正協方差公式， 構造無偏形式的誤差協方差矩陣， 從而在進行分析時減少對同化結果的影響。

1.2.3 隨機子采樣方法

由于擴展后的分析集合數量M大于原始預測集合數量N， 導致集合數量不能在預測模型中進行傳遞。GETKF方法中利用一個膨脹因子[式(12)]來克服此問題， 該方法需要計算擴展后的分析誤差協方差矩陣以及原始預測誤差協方差矩陣， 如果變量的個數n很大， 計算難度大大增加。本文使用隨機子采樣方法進行N列分析集合的選取， 計算簡便。隨機子采樣方法如下：

1.2.4 GC局地化函數的限制

GETKF方法中使用GC函數作為局地化函數，但是GC函數僅僅是一個一元變量的局地化函數， 多元模型中不同變量之間的相關性應該是有區別的，GC函數并沒有體現這一性質， 所以在多元模型中使用GETKF方法存在限制。本文選用一個基于距離的多元函數Askey函數[21]作為局地化函數， 當前此函數已經用于EnKF方法， 并且在多元模型中進行了實驗驗證， 得到了較好的同化效果。利用Askey函數將本文修改的方法擴展到多元模型變量中， 探究方法在多元情況下的適用性。

2 實驗方案

2.1 模型

2.1.1 KS模型

KS模型的表達式為，

KS模型是一個一元無量綱的非線性偏微分方程模型， 方程中二階項和四階項會增加模型的復雜性并導致混沌行為。在本文中， KS模型的偏微分方程通過一個指數時間的四階龍格庫塔數值格式進行離散， 時間步為Δt= 0.25， 終止時間T為250。

2.1.2 淺水模型

淺水模型作為一個多元模型， 在當前的數據同化研究中廣泛使用。其忽略摩擦效應、科氏力以及非線性項的方程表示為：

其中，u和v分別表示水平和垂直速度， 水面高度由h表示，g為取值為9.8 m/s2的重力加速度。模型區域選擇矩形區域取2 200 km。淺水模型的偏微分方程組使用Lax-Wendroff有限差分方法計算， 時間步為Δt= 0.01， 終止時間T為8。

2.2 局地化函數

對于KS方程等一元模型， 我們可以使用GC函數作為CL方法中的局地化函數ρ， 其表達式如下：

式中，z表示網格點之間或網格點與觀測點之間的距離。由于本文使用模型空間CL方法， 所以z表示物理空間中網格點之間的距離。局地化長度尺度c定義為為局地化半徑， 是一個可選參數。注意， 因子可以調整局地化函數達到最優[22]。

相反， 多元模型中使用Askey函數的一個多元擴展[21]作為局地化函數ρ， 形式為，

其中，m表示不同的狀態變量的總數。例如本文使用的淺水方程中m取3， 則Askey函數表示為：

即ρ的每個元素ijρ對應于第i個變量和第j個變量的局地化函數矩陣。由于多元模型引入了一些不同變量之間的聯系， 增加了數據同化的復雜度與計算量， 因此多元局地化函數的構造需要考慮這些不同變量之間的聯系以及聯系的強弱。uij表示第i個和第j個變量之間的局地化影響參數，

Askey函數中的參數z和c與GC函數中定義相同， 分別表示網格點之間的距離和局地化長度尺度。B是一個Beta函數，s表示狀態變量所在的歐式空間的維度， 而且要滿足條件v≥(s+ 1)/2。

2.3 實驗設計

其中， 矩陣H是觀測算子， 將真實狀態變量映射到觀測空間。假定觀測誤差之間不相關， 則矩陣R為對角觀測誤差協方差矩陣， 且對角元素取值為1。初始集合通過對真實值t x添加N次隨機誤差構造。詳細的實驗參數設置如下。

KS模型具有一維周期性， 設定狀態向量u的維度n為256， 初始真實值可以通過定義在周期域0≤x≤32π上的函數給出。實驗中分別取集合數N為5或N為10， 觀測數p是256或p是235。數值實驗的觀測頻率(同化間隔)不同， 假定每5或10個時間步長存在可用的觀測值。

對于多元淺水模型， 假定模型定義在矩形網格區域， 每個網格點定義3個變量： 水面高度h， 水平速度u和垂直速度v。假定水平速度u和垂直速度v的初始真實狀態為零， 即u=v= 0， 水面高度h的初始真實狀態由如下方式構造，

類似于KS模型中的參數設置， 在淺水方程中，集合數分別取N是5或N是10， 觀測數p為300或p為240， 每5或10個時間步長存在可用的觀測值。表1中給出了上述具體實驗參數的配置， 分別在兩

表1 不同實驗條件的參數Tab. 1 Parameters for different experimental conditions

個模型中進行8組不同的實驗， 其中對于KS方程全觀測p為256， 部分觀測p為235； 對于淺水模型全觀測p為300， 部分觀測p是240。

3 實驗結果與分析

首先， 分析GETKF方法和本文修改的方法(簡稱GCL方法)在一元KS模型中的實驗結果。為了研究兩種方法對預測誤差協方差矩陣的影響， 以實驗1為例， 取前兩個同化時刻(存在觀測數據的時刻)，比較同化過程中兩種方法對于預測誤差協方差矩陣的影響。在圖1至圖4的每一幅圖中， 從左到右，從上到下， 依次展示了GC局地化函數、預測誤差協方差矩陣、局地化函數和預測誤差協方差矩陣做舒爾積運算得到的矩陣、使用局地化函數平方根矩陣與預測集合擾動矩陣做舒爾積運算得到的預測誤差協方差矩陣。

圖1 實驗1中， 第一個同化時刻GCL對預測誤差協方差矩陣的影響Fig. 1 Effect of GCL on the forecast error covariance matrix at the first assimilation time in experiment 1

圖3 實驗1中， 第一個同化時刻GETKF對預測誤差協方差矩陣的影響Fig. 3 Effect of GETKF on the forecast error covariance matrix at the first assimilation time in experiment 1

圖4 實驗1中， 第二個同化時刻GETKF對預測誤差協方差矩陣的影響Fig. 4 Effect of GETKF on the forecast error covariance matrix at the second assimilation time in experiment 1

由于KS模型具有一維周期性， 所以GC局地化函數ρ是一個對稱矩陣， 且僅在主對角線附近和矩陣邊緣存在較強的相關性。很顯然， 對于GETKF和GCL這兩種方法， 圖1至圖4中的第4幅子圖均顯示了使用局地化函數之后協方差矩陣遠距離的偽相關得到消除， 并且和第3幅子圖相比效果類似， 說明兩種方法對于虛假相關性的消除是有效果的。但是，比較圖2和圖4， 可以看出第2個時刻兩種方法的局地化效果是存在偏差的(圖例中數值區間不同)， 這可能是源于兩種方法的具體實施不同。但是總體上兩者對于處理偽相關的效果較好， 即截斷距離以外的相關性被消除， 鄰近狀態變量和邊界上的相關性得到了保留。

圖2 實驗1中， 第二個同化時刻GCL對預測誤差協方差矩陣的影響Fig. 2 Effect of GCL on the forecast error covariance matrix at the second assimilation time in experiment 1

以實驗3和實驗7為例， 分別對GETKF方法和GCL兩種方法的同化效果進行評估。如圖5和圖6中不同方法的曲線軌跡所示， 大體上可以看出相比于GETKF方法的同化效果， 本文修改的GCL方法更加接近于真實值的軌跡曲線， 并且減少了濾波發散的發生， 展示了良好的同化效果， 說明在一元數值預測模型下， GCL方法的局地化效果優于GETKF方法。

圖5 實驗3， KS模型第一個變量的同化效果Fig. 5 Comparison of the assimilation effect of the first variable of the KS model in experiment 3

圖6 實驗7中， KS模型第一個變量的同化效果Fig. 6 Comparison of the assimilation effect of the first variable of the KS model in experiment 7

均方根誤差(root mean square error， RMSE)可以量化不同局地化方法的效果， 直觀表示不同方法的優劣程度。模型狀態變量的整體均方根誤差表示為，

其中，l表示同化時間的長度，表示第i個時刻第j個變量的分析值，表示第i個時刻第j個變量的真實值。此外， 由預測集合通過統計意義計算得到預測誤差協方差矩陣的秩至多為N-1， 該矩陣是一個低秩的矩陣， 這在傳遞誤差時易造成矩陣的病態， 因此需要提升矩陣的秩。理論上CL方法對于矩陣秩的提升是有效果的， 但是由于CL方法應用于ETKF方法的局限性， 使得效果大打折扣， 所以除了計算同化過程中的RMSE外， 也會對比GETKF和GCL方法使用前后預測誤差協方差矩陣秩的變化。具體如表2所示， 首先， 兩種方法的預測誤差協方差矩陣的秩與不使用局地化方法相比， 顯然均有了提高。GCL方法對于矩陣秩的提高更為明顯。其次， GCL方法的RMSE明顯小于GETKF方法的RMSE(實驗6有偏差)， 說明在一元模型中， GCL方法的同化效果要優于GETKF方法的同化效果。

表2 KS模型中不同實驗條件下RMSE和矩陣秩的對比Tab. 2 Comparison of the RMSE and matrix rank under different experimental conditions in the KS model

續表

表3 中展示了主流的數據同化方法(EnKF)的RMSE， 采用實驗1中的觀測數和同化間隔， 同時使用變化的集合數進行數值實驗。EnKF方法在集合數N是5的情況下， RMSE的數值遠遠大于GCL方法(實驗1)， 隨著集合數N增大到100， RMSE的數值在降低， 與GCL方法(實驗1)的RMSE相當， 但是程序的運行時間也在不斷增加。當集合數為N取200時，EnKF方法的RMSE數值低于GCL方法(實驗1)， 表現出了良好的數據同化效果， 但是程序運行時間是GCL方法使用集合數N為5時的數倍(實驗1條件下，GCL方法運行時間為4.83 s)。EnKF通過使用大集合數， 犧牲程序計算時間來獲取低RMSE， 并不是一種很好的處理方式， 并且在真實的數據同化方法應用中， 選取超過100個集合是不合適的。上述實驗結果表明兩種數據同化方法各有優劣， 但是可以看出在集合數較小的情況下， 選擇GCL方法進行數據同化的效果較好。

表3 KS模型中EnKF方法的RMSETab. 3 RMSE of the EnKF method in the KS model

最后， 我們將本文提出的GCL方法中的GC函數替換為Askey函數， 擴展方法的適用范圍， 在多元淺水模型中檢驗該方法。表4展示了在不同實驗條件下， GCL方法與未使用局地化方法的RMSE以及預測誤差協方差矩陣的秩的對比。類似于一元數值模型中的實驗結果， 使用GCL方法后的預測誤差協方差矩陣的秩與未使用局地化方法相比， 有了明顯的提升， 保證了誤差矩陣傳遞的可靠性， 從而GCL方法的RMSE與未使用局地化方法相比也有了明顯的降低， 說明利用Askey函數將GCL方法擴展到多元模型是成功的， 克服了GETKF方法只能用于一元模型的局限性。

表4 淺水模型中不同實驗條件下RMSE和矩陣秩的對比Tab. 4 Comparison of the RMSE and matrix rank under different experimental conditions in the shallow water model

4 結論

為解決基于集合的卡爾曼濾波中集合數目過少導致的欠采樣問題， 以及CL方法對于ETKF等平方根卡爾曼濾波方法的不適用， 本文對GETKF方法進行了改進， 提出GCL方法， 并與GETKF方法在相同實驗條件下進行了實驗對比。結果表明， GCL方法改善了在一元模型情況中的同化效果， 明顯提高了預測誤差協方差矩陣的秩， 降低了分析值的均方根誤差。此外， 利用一個Askey函數實現局地化方法在多元模型中的擴展， 使得提出的GCL方法可以用于多元模型中， 克服當前多元模型中局地化方法的欠缺?？傮w來說， 雖然GCL方法在處理欠采樣問題以及對于多元模型局地化具有一定的效果， 但是GCL方法應用到多元模型時， 為設計多元變量之間的相互作用關系及相關性， Askey函數需要預先指定的參數較多。在下一步的研究中應找到合適的參數確定依據，使得提出的新方法具有更廣泛的適用性。