基于Lorenz-96 的數(shù)據(jù)同化基礎算法數(shù)值實驗系統(tǒng)設計與實現(xiàn)

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

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

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

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

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

1 系統(tǒng)設計

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

1.1 算法模塊

早期同化方法是指出現(xiàn)時間較早、實現(xiàn)簡單的方法。變分同化算法和卡爾曼濾波算法是當前業(yè)務化應用的主流方法。粒子濾波算法由于不需要高斯分布假設，近些年越來越受到關注。混合方法是把前述的算法進行混合，充分利用不同算法的優(yōu)點，也是當前同化方法的一種趨勢[7]。由于方法名字本身的中文翻譯沒有統(tǒng)一的標準，本文同時給出各同化方法的英文。

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

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

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

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

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

該實驗系統(tǒng)的主要功能模塊如圖1 所示。

圖1 數(shù)據(jù)同化基礎算法實驗系統(tǒng)功能模塊

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

1.2 數(shù)學模型

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

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

2 系統(tǒng)實現(xiàn)

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

2.1 實驗系統(tǒng)界面

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

圖2 數(shù)據(jù)同化基礎算法數(shù)值實驗系統(tǒng)主界面

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

2.2 實驗系統(tǒng)功能解析

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

2.2.1 算法流程

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

（4）分析步和預測步循環(huán)進行。

2.2.2 程序?qū)崿F(xiàn)

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

Function [xa,mean_xf,mean_xa] =

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

%% EnKF 的分析步

%% 輸入:

%% xf0 -預報狀態(tài)矩陣

%% y - 觀測向量

%% H - 觀測算子

%% R - 觀測誤差協(xié)方差矩陣

% %n -集合數(shù)

%% 輸出:

% xa -分析的集合.

% mean_xf - 預報集合的均值.

% mean_xa - 分析集合的均值.

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

%預報集合均值

mean_xf = mean(xf0,2);

% 預報誤差協(xié)方差矩陣

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

% 觀測數(shù)

p = size(y,1);

% 觀測集合隨機誤差

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

% 觀測集合

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

% 觀測誤差協(xié)方差矩陣

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 所示），實驗系統(tǒng)彈出EnKF 算法子模塊的界面。

圖3 EnKF 二級菜單

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

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

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

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

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

function [H] = Generate_H(p,n)%%產(chǎ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 子模塊也可修改各輸入?yún)?shù)以查看算法的同化效果，分析各參數(shù)的敏感性。甚至，用戶可以直接修改Maltab 的“.m” 文件來實現(xiàn)定制的個性化算法，以滿足教學和科研的需要。本文僅給出了EnKF 子模塊的功能展示，根據(jù)算法原理的不同，其他子模塊的界面、輸入?yún)?shù)、輸出結(jié)果的展示方式等略有不同。培養(yǎng)學生的能力是大學教學的核心[13]，本系統(tǒng)的所有參數(shù)、代碼均可修改，可培養(yǎng)學生的自主學習和動手能力，為數(shù)據(jù)同化課程的研究性實驗教學提供了很好的實驗平臺。

3 結(jié)語

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

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