基于Lorenz-96 的數據同化基礎算法數值實驗系統設計與實現

王際朝，臧紹東，阮宗利，吳國麗

（中國石油大學（華東） 理學院，山東 青島 266580）

數據同化（data assimilation）亦稱資料同化，根據一定的標準和方法，將不同時空、不同分辨率、不同手段獲取的觀測數據與動力模式有機結合，建立數據與模型相互協調的優化關系，從而給出當前模式狀態的變量分析場。作為一種技術，數據同化已廣泛應用于大氣[1]、海洋[2]、陸面[3]、石油工程[4]等諸多領域。尤其是當前各領域的觀測數據飛速積累，數值模擬技術也日趨成熟，數據同化作為連接觀測和模式的橋梁，其作用愈發重要。

國內很多高校，尤其是涉及大氣、海洋、石油等專業的院校，面向高年級本科生或研究生開設數據同化課程。數據同化課程的教學內容是一系列同化方法或者算法的理論講解和上機實驗。此課程的理論基礎是估計理論和控制理論[5]，教師在教學過程中易偏重理論講解，大量的復雜公式使學生感到抽象而缺乏興趣，很難發揮學生學習的積極性和主動性。算法是數據同化課程的本質，在完成算法詳細講解的同時，教師應適當利用計算機可視化技術，使得教學內容更形象和生動。

本課題組查閱了大量資料，發現鮮有對此課程進行教學研究的文獻，這可能是由于數據同化課程的“科研” 屬性較強。本文基于Lorenz-96 模型，設計并實現了數據同化課程基礎算法數值實驗系統，旨在為數據同化教學過程中基礎算法的講解提供一個可視化的載體，從而加深學生對算法的理解程度，培養學生的動手能力，提高學生的學習興趣[6]，進而提升教學質量。

1 系統設計

考慮到數據同化各種方法的理論基礎、出現時間、是否業務化應用主流、是否當前研究前沿等，本文將數據同化課程中涉及的基礎算法分為5 類：早期同化方法、變分同化算法、卡爾曼濾波算法、粒子濾波算法和混合方法。所有算法都在統一的數學模型下進行程序設計，以方便比較各種方法的同化效果。

1.1 算法模塊

早期同化方法是指出現時間較早、實現簡單的方法。變分同化算法和卡爾曼濾波算法是當前業務化應用的主流方法。粒子濾波算法由于不需要高斯分布假設，近些年越來越受到關注?；旌戏椒ㄊ前亚笆龅乃惴ㄟM行混合，充分利用不同算法的優點，也是當前同化方法的一種趨勢[7]。由于方法名字本身的中文翻譯沒有統一的標準，本文同時給出各同化方法的英文。

早期同化方法包括：多項式插值（polynomial interpolation）、逐步訂正法（successive corrections method）、松弛逼近法（nudging method）和最優插值（optimal interpolation，OI）法。

變分算法包括： 三維變分（ three-dimensional variational，3D-VAR）和四維變分（four-dimensional variational，4D-VAR）算法。

卡爾曼濾波算法包括：卡爾曼濾波（Kalman filter，KF）、擴展卡爾曼濾波（extended Kalman filter，EKF）、集合卡爾曼濾波（ensemble Kalman filter，EnKF）、集合調整卡爾曼濾波（ensemble adjustment Kalman filter，EAKF）、集合轉換卡爾曼濾波（ensemble transform Kalman filter，ETKF）和局地化集合轉換卡爾曼濾波（local ensemble transform Kalman filter，LETKF）算法。

粒子濾波算法包括：基本粒子濾波（particle filter，PF）、重采樣粒子濾波（resampling particle filter，RPF）、隱式粒子濾波（implicit particle filter，IPF）、等權重粒子濾波（equivalent-weights particle filter，EWPF）和映射粒子濾波（mapping particle filter，MPF）算法。

混 合 方 法 包 括 ： 4DVAR+EnKF 、 PF+EnKF 、4DVAR+KF 和EnVarPF 方法。

該實驗系統的主要功能模塊如圖1 所示。

圖1 數據同化基礎算法實驗系統功能模塊

數據同化的算法仍在不斷改進和發展中，部分新的算法或者非常復雜的算法并未包括在本實驗系統中。本文設計的模塊基本涵蓋了數據同化課程的基礎算法，同時兼括了部分近幾年發展起來的同化算法。

1.2 數學模型

數據同化算法的設計與實現需要數學模型作為基礎。經過多年數據同化課程教學經驗的積累，本文選擇Lorenz-96 模型作為所有算法的模型基礎，即本文的所有算法均基于此數學模型進行編程和設計。Lorenz-96 模型是高度非線性的，且對初值極端敏感，是驗證各種同化算法的理想模型[8]。此模型是Lorenz等于1996 年從動力模式中推導出的簡化的動力模型，被廣泛地用來驗證各類同化方法。其動力學方程為

其中：F為外強迫項，n為變量個數。Lorenz-96 系統的線性耗散項使得總能量迫項F能增大或降低總能量，二次平流項對貢獻，具有保存總能量的性質。

2 系統實現

Matlab 軟件為數據同化課程的實驗教學提供了強有力的設計平臺[9]。 GUI 是指圖形用戶界面（graphical user interface，GUI），由窗口、按鈕、文字、圖形、菜單等對象構成的用戶視窗，是實現人機交互的有效工具[10-11]。本文利用Matlab GUI，設計了數據同化基礎算法數值實驗系統。

2.1 實驗系統界面

Matlab GUI 形成包含GUI 組件屬性的擴展名為“.fig” 的文件和保存程序代碼的擴展名為“.m” 的文件[12]。GUI 程序主要包括各控件的初始化函數、輸出函數、布局以及回調函數。實驗系統的主界面采用經典的菜單模式，主菜單上的5 個部分對應5 個功能模塊，點擊可以彈出各子模塊對應的二級子菜單。單擊每個子菜單會調出各子模塊對應的界面。實驗系統主界面如圖2 所示。

圖2 數據同化基礎算法數值實驗系統主界面

點擊主菜單的下拉二級菜單，會調出各個算法的執行模塊。由于算法眾多，本文優選出集合卡爾曼濾波算法進行實驗系統的功能解析，并給出算法的原理和主要程序的說明。

2.2 實驗系統功能解析

集合卡爾曼濾波（EnKF）是由Geir Evensen 利用蒙特卡洛思想提出的，隨后涌現了基于此方法思想的各種集合卡爾曼濾波類同化方法。集合卡爾曼濾波在大量的理論和實際業務化數據同化問題上已經被證明是非常有效的。

2.2.1 算法流程

EnKF 算法中，計算預測誤差協方差矩陣是通過集合成員完成的，因此能夠巧妙地避開切線性算子的求解問題，對非線性系統有著更好的適用性。這種處理方式不僅解決了KF 方法應對非線性系統的局限，同時保留了誤差協方差矩陣隨時間演化的特點。由于EnKF 方法中對于除了卡爾曼增益矩陣的計算外其他所有在集合成員之間的操作都是互相獨立的，因此很容易運用計算機進行并行計算，這是EnKF 方法得到推廣的重要原因。EnKF 算法流程如下：

（4）分析步和預測步循環進行。

2.2.2 程序實現

在EnKF 算法中，數據同化分析步是整個算法的核心。同化分析步作為一個函數 “EnKF_analysis”，供實驗系統調用。用戶只需輸入預測狀態矩陣、觀測向量、觀測算子、觀測誤差協方差矩陣和集合數，即可得到同化分析值等。其Matlab 程序如下：

Function [xa,mean_xf,mean_xa] =

EnKF_analysis(xf0, y, H, R,n)

%% EnKF 的分析步

%% 輸入:

%% xf0 -預報狀態矩陣

%% y - 觀測向量

%% H - 觀測算子

%% R - 觀測誤差協方差矩陣

% %n -集合數

%% 輸出:

% xa -分析的集合.

% mean_xf - 預報集合的均值.

% mean_xa - 分析集合的均值.

%%----------------------------------------------%%

%預報集合均值

mean_xf = mean(xf0,2);

% 預報誤差協方差矩陣

Pf = ((xf0 - repmat(mean_xf,1,n)) * (xf0 -repmat(mean_xf,1,n))') / (n-1);

% 觀測數

p = size(y,1);

% 觀測集合隨機誤差

u = 0 + sqrtm(R) * randn(p,n);

% 觀測集合

y_ensemble = repmat(y,1,n) + u;

% 觀測誤差協方差矩陣

Ru = (u * u') / (n-1);

% 集合卡爾曼增益

Ku = Pf * H' * inv((H * Pf * H' + Ru));

% 分析集合

xa = xf0 + Ku * (y_ensemble - H * xf0);

% 分析均值

mean_xa = mean(xa,2);

end

2.2.3 功能演示

在主界面下點擊菜單“卡爾曼濾波算法”，在彈出的下拉菜單中點擊 “Ensemble Kalman Filter (EnKF)”（如圖3 所示），實驗系統彈出EnKF 算法子模塊的界面。

圖3 EnKF 二級菜單

在彈出的EnKF 算法界面中，輸入變量數、集合數、同化間隔、觀測算子等初始化參數后，點擊“Run_EnKF” 按鈕，如圖4 所示，實驗系統自動給出同化效果的結果圖形。效果圖展示了EnKF 同化觀測數據后與真值比較的結果，其中黑色虛線表示 “真值”，紅色實線表示EnKF 的分析值，藍色方框表示輸入參數規定的觀測數據。

圖4 第1 個變量的同化效果

圖4 展示的是第1 個變量的同化效果。在不改變其他輸入參數的情況下，實驗系統亦可展示其他變量的同化效果。圖5 給出的是第19 個變量的同化效果。從效果圖4 和5 中可以看出，EnKF 算法具有較好的同化效果。根據計算結果，亦可分析均方根誤差、相關系數等誤差統計指標值。

圖5 第19 個變量的同化效果

此EnKF 子模塊可對輸入參數進行任意符合規則的設置，提高了算法的人機交互性能。由于EnKF 算法輸入參數較多，有些參數是作為子程序進行輸入的，比如觀測算子H的構造，其本身就是一個子程序，輸入和輸出要符合實驗系統的要求。這里給出“參數”Generate_H 的Matlab 程序。實驗系統的用戶可在后臺修改這些參數，以滿足個性化的需求。

function [H] = Generate_H(p,n)%%產生觀測算子矩陣

H = zeros(p,n); %預分配空間

if rem(n,p) == 0; %如果能整除

j = 1:n/p: n;

for i = 1: p

H(i,j(i)) = 1;

end else

sd=1; %隨機種子

rng(sd);

k_randn = randperm(n,p);

[k_randn] = sort(k_randn,'ascend');

j = k_randn; for i = 1 : p

H(i,j(i)) = 1;

end

end

end

除了可以展示前述功能外，EnKF 子模塊也可修改各輸入參數以查看算法的同化效果，分析各參數的敏感性。甚至，用戶可以直接修改Maltab 的“.m” 文件來實現定制的個性化算法，以滿足教學和科研的需要。本文僅給出了EnKF 子模塊的功能展示，根據算法原理的不同，其他子模塊的界面、輸入參數、輸出結果的展示方式等略有不同。培養學生的能力是大學教學的核心[13]，本系統的所有參數、代碼均可修改，可培養學生的自主學習和動手能力，為數據同化課程的研究性實驗教學提供了很好的實驗平臺。

3 結語

基于Lorenz-96 數學模型，利用Matlab GUI，設計和實現了數據同化課程基礎算法數值實驗系統，并以EnKF 算法為例展示了實驗系統的部分功能。系統界面友好，交互簡單，內容豐富，流程清晰，演示和擴展方便。師生輸入算法所需參數即可獲得可視化實驗結果，并可對比分析不同條件下、不同算法的數據同化效果，使得原來算法的枯燥理論推導過程可在“所見即所得” 的環境下變得更加生動和形象。本實驗系統還應用到了線上教學活動中，豐富了教學手段，取得了理想的教學效果。

數據同化課程基礎算法數值實驗系統為課程理論教學和實驗教學的有機結合提供了抓手，為實驗教學提供了強有力的載體，豐富了教學手段，可用在數據同化課程的的輔助性教學中以提升教學質量。本系統亦可為學生提供一個同化算法的數值實驗平臺，有益于加深學生對算法的理解和提高學生的動手能力。