MASNUM-WAM海浪模式集合Kalman濾波同化研究*
——II. 集合樣本對(duì)同化效果的影響

孫 盟 尹訓(xùn)強(qiáng), 楊永增,① 吳克儉 孫寶楠

(1. 中國海洋大學(xué)海洋與大氣學(xué)院 青島 266100; 2. 國家海洋局第一海洋研究所海洋環(huán)境與數(shù)值模擬研究室青島 266061; 3. 海洋國家實(shí)驗(yàn)室區(qū)域海洋動(dòng)力學(xué)與數(shù)值模擬功能實(shí)驗(yàn)室 青島 266071)

海浪數(shù)據(jù)同化是改善海浪模式和提高海浪預(yù)報(bào)水平的重要途徑之一。影響同化效果的關(guān)鍵因素之一是背景誤差的確定。背景誤差是模擬與真值的偏差。由于真值是未知的, 通常需要對(duì)背景誤差進(jìn)行近似構(gòu)造和估計(jì), 不同的同化方法對(duì)此有各自的處理方式。最優(yōu)插值方法(Lorenc, 1981)將背景誤差作為一個(gè)常值, 同化過程中僅對(duì)其進(jìn)行一次估計(jì), 即靜態(tài)統(tǒng)計(jì)。集合Kalman濾波方法(Kalman, 1960; Kalman et al,1961; Evensen, 1994; Burgers et al, 1998)是當(dāng)前較為先進(jìn)的同化方法, 同化過程中, 多個(gè)模式同時(shí)運(yùn)行,模式集合可以對(duì)背景誤差進(jìn)行實(shí)時(shí)更新和動(dòng)態(tài)估計(jì)。由此, 集合Kalman濾波(Ensemble Kalman Filter, 簡稱EnKF)對(duì)計(jì)算資源的要求較高?；趥鹘y(tǒng)集合濾波方法, Anderson(2001, 2003)提出了集合調(diào)整 Kalman濾波(Ensemble Adjustment Kalman Filter, 簡稱EAKF)。 該方法無需對(duì)觀測(cè)進(jìn)行擾動(dòng),其引入了一個(gè)線性算子, 用于替代傳統(tǒng)的增益矩陣。當(dāng)集合數(shù)目較小時(shí)(10—20個(gè)), EAKF方法的表現(xiàn)優(yōu)于 EnKF(Evensen, 2003)。

在海洋和大氣領(lǐng)域, 已有采用EAKF方法的相關(guān)研究工作(Zhang et al, 2005; Karspeck et al, 2007;Zhang et al, 2007; Yin et al, 2010, 2011; Karspeck et al,2013; Chen et al, 2016; Hill et al, 2016; Shen et al,2016)。 但具體到海浪資料同化, 相關(guān)研究較少。孫盟等(2014)提出基于靜態(tài)樣本的集合 Kalman濾波同化方法, 采用24h間隔有效波高模擬偏差構(gòu)造誤差集合, 用于近似背景誤差, 將該誤差集合與模式變量相疊加, 構(gòu)成模式狀態(tài)變量集合, 結(jié)合兩步濾波方法得到再分析的波高場(chǎng), 同化效果顯著。靜態(tài)樣本集合Kalman濾波方法僅需運(yùn)行一個(gè)海浪模式, 大大降低了對(duì)計(jì)算資源的需求, 適用于業(yè)務(wù)化海浪預(yù)報(bào)。為了量化該方法與傳統(tǒng)集合濾波方法的差異, 本文分別采用靜態(tài)樣本集合Kalman濾波和EAKF方法, 開展海浪資料同化實(shí)驗(yàn)。EAKF同化方案需要構(gòu)造具有代表性的模式集合。 考慮到海浪模式對(duì)于初始場(chǎng)的敏感性較弱, 本文對(duì)風(fēng)場(chǎng)進(jìn)行隨機(jī)場(chǎng)集合擾動(dòng)(孫盟等,2016), 由風(fēng)場(chǎng)集合驅(qū)動(dòng)生成海浪模式集合。

本文基于MASNUM-WAM海浪模式, 針對(duì)2014年全球海域, 分別采用靜態(tài)樣本集合 Kalman濾波和EAKF同化方法, 利用Jason-2衛(wèi)星高度計(jì)資料, 開展海浪資料同化實(shí)驗(yàn)。上述為引言, 第一節(jié)為模式與數(shù)據(jù)介紹, 第二節(jié)為同化方案與實(shí)驗(yàn)設(shè)計(jì), 其中詳細(xì)介紹了靜態(tài)樣本集合Kalman濾波同化方案、EAKF同化方案及實(shí)驗(yàn)設(shè)計(jì), 第三節(jié)為實(shí)驗(yàn)結(jié)果與分析, 最后一節(jié)為結(jié)論。

1 模式與數(shù)據(jù)

本文采用球坐標(biāo)系下的第三代海浪模式MASNUM-WAM(Marine Science and Numerical Modeling) (Yuan et al, 1991, 1992; 楊永增等, 2005)。該模式應(yīng)用了基于破碎波統(tǒng)計(jì)理論發(fā)展的海浪破碎耗散源函數(shù)(Yuan et al, 1986), 并采用復(fù)雜特征線嵌入計(jì)算格式。球坐標(biāo)系下, 波譜能量平衡方程和復(fù)雜特征線方程如式(1)和式(2)(3)(4)所示,

其中, 波數(shù)譜 E =E（ K ,λ , φ,t）為波數(shù)、經(jīng)度λ和緯度φ的函數(shù);表示背景流場(chǎng);）表示群速度。等式(1)右側(cè)包含風(fēng)輸入項(xiàng)、破碎耗散項(xiàng)、底摩擦耗散項(xiàng)、非線性波-波相互作用項(xiàng)和波流相互作用項(xiàng)。θ1為波數(shù)矢量輻角(指向東為 0, 逆時(shí)針旋轉(zhuǎn)為正),）表示波矢量的單位矢量,）表示頻散關(guān)系。

波數(shù)譜被離散成24個(gè)方向和25個(gè)波數(shù), 對(duì)應(yīng)頻率范圍是 0.042—0.413Hz。本文模式計(jì)算區(qū)域?yàn)?60°S—60°N; 0°—360°, 空間分辨率為 0.5°×0.5°, 時(shí)間步長為 15min, 模式輸出為 1h一次。地形數(shù)據(jù)為ETOPO5。風(fēng)場(chǎng)數(shù)據(jù)采用2014年全球ECMWF風(fēng)場(chǎng)(ERA-Interim), 該數(shù)據(jù)由歐洲中期天氣預(yù)報(bào)中心提供。風(fēng)場(chǎng)數(shù)據(jù)時(shí)間間隔為 6h, 空間分辨率為 0.5°×0.5°。

海浪同化資料采用 Jason-2衛(wèi)星高度計(jì)數(shù)據(jù)。Jason-2為Topex/Poseidon和Jason-1的后繼衛(wèi)星, 于2008年 6月發(fā)射升空, 工作至今, 重復(fù)周期約為 10天, 測(cè)高精度為2.5cm。海浪同化檢驗(yàn)資料采用Saral衛(wèi)星高度計(jì)數(shù)據(jù)。Saral衛(wèi)星隸屬于印度空間研究組織(ISRO), 由法國國家空間研究中心(CNES)研制,其軌道路線沿襲了ERS-1/2和Envisat衛(wèi)星的軌道設(shè)置。該衛(wèi)星于2013年2月發(fā)射升空, 工作至今, 重復(fù)周期約為35天, 測(cè)高精度為8mm。上述衛(wèi)星資料由AVISO(Archiving, Validation and Interpolation of Satellite Oceanographic)提供。

2 同化方案與實(shí)驗(yàn)設(shè)計(jì)

本文分別采用靜態(tài)樣本集合 Kalman濾波和EAKF方法, 針對(duì) 2014年全球海域開展海浪資料同化實(shí)驗(yàn)。靜態(tài)樣本集合Kalman濾波方案介紹詳見2.1節(jié); EAKF方案介紹詳見2.2節(jié)。海浪同化實(shí)驗(yàn)設(shè)計(jì)介紹詳見2.3節(jié)。

2.1 靜態(tài)樣本集合Kalman濾波方案

靜態(tài)樣本集合Kalman濾波(孫盟等, 2014)同化過程可分為兩個(gè)部分。首先利用海浪模式歷史數(shù)據(jù), 由24h間隔有效波高模擬偏差構(gòu)造誤差集合, 如式(5)所示, 該誤差集合可用于近似背景誤差, 將此誤差集合與同化時(shí)刻的模式變量相疊加, 構(gòu)成模式狀態(tài)變量集合。然后, 采用兩步濾波方法對(duì)模式狀態(tài)變量集合進(jìn)行更新, 生成再分析的波高場(chǎng)。該方法僅需運(yùn)行一個(gè)海浪模式, 可以大幅度降低計(jì)算成本, 適用于業(yè)務(wù)化海浪預(yù)報(bào)。本文利用2013年12月份海浪模式結(jié)果,構(gòu)造靜態(tài)樣本集合, 已知模式結(jié)果輸出頻率為1h, 故集合樣本數(shù)目為720個(gè)。

其中, Cov(x,y)表示空間兩點(diǎn)x和y的協(xié)方差, 上角標(biāo)t和t–24h表示時(shí)間, E表示期望。

2.2 EAKF方案

(1) EAKF方法簡介

根據(jù)上述定義, 聯(lián)合向量的概率密度函數(shù)可表示為

其中,tkY 表示t時(shí)刻以前的所有觀測(cè)及t時(shí)刻的k個(gè)觀測(cè)子集。結(jié)合貝葉斯公式, 得到

上式表示, 一個(gè)新的觀測(cè)資料子集合將如何改變基于已有觀測(cè)資料更新得到的初猜聯(lián)合向量。

當(dāng)假設(shè)觀測(cè)資料誤差和先驗(yàn)概率密度分布為高斯分布時(shí), 式(8)可視為兩個(gè)高斯過程的卷積,結(jié)果仍為高斯分布, 新高斯過程的協(xié)方差矩陣和均值為

其中, R代表觀測(cè)誤差協(xié)方差矩陣。引入線性算子A(滿足 ∑u= A∑pAT), 用于更新聯(lián)合向量的每個(gè)樣本。

EAKF方法的主要優(yōu)點(diǎn)在于無需對(duì)觀測(cè)進(jìn)行樣本化處理, 盡可能的保留經(jīng)驗(yàn)分布中的非線性高階信息, 此外, 當(dāng)集合數(shù)目相對(duì)較小時(shí), 該方法能夠獲得合理可靠的結(jié)果(Anderson, 2001)。

(2) 模式集合的構(gòu)造

在EAKF方案中, 多個(gè)海浪模式同時(shí)運(yùn)行, 該方法的關(guān)鍵是構(gòu)造具有代表性的海浪模式集合?？紤]到海浪模式對(duì)于初始場(chǎng)的敏感性較弱, 此處對(duì)風(fēng)場(chǎng)進(jìn)行持續(xù)地集合擾動(dòng), 由風(fēng)場(chǎng)集合驅(qū)動(dòng)生成海浪模式集合。風(fēng)場(chǎng)集合擾動(dòng)采用隨機(jī)場(chǎng)擾動(dòng)方案(孫盟等,2016), 如式(12)所示, 風(fēng)場(chǎng)集合由當(dāng)前時(shí)刻風(fēng)場(chǎng)矢量和擾動(dòng)集合疊加而成, 擾動(dòng)集合為三維隨機(jī)場(chǎng), 該隨機(jī)場(chǎng)根據(jù)蒙特卡羅方法生成(Evensen, 1994), 其整體滿足正態(tài)分布, 單個(gè)隨機(jī)場(chǎng)樣本的空間分布具有局地性和平滑性。

2.3 實(shí)驗(yàn)設(shè)計(jì)

對(duì) 2014 年全球海域(60°S—60°N, 0°—360°)開展海浪集合濾波同化實(shí)驗(yàn), 實(shí)驗(yàn)分為三組, 分別為未加入同化的控制實(shí)驗(yàn)、靜態(tài)樣本集合Kalman濾波同化實(shí)驗(yàn)和 EAKF同化實(shí)驗(yàn)(集合數(shù)目為 10個(gè))。同化數(shù)據(jù)為 Jason-2衛(wèi)星高度計(jì)資料, 同化結(jié)果檢驗(yàn)數(shù)據(jù)為Saral衛(wèi)星高度計(jì)資料。

海浪同化系統(tǒng)運(yùn)行流程如圖1所示。海浪模式初始化后, 每間隔6h同化一次。同化過程中, 由單個(gè)海浪模式(靜態(tài)樣本集合 Kalman濾波方案)或海浪模式集合(EAKF方案)輸出控制文件和海浪譜, 調(diào)用同化模塊, 待控制文件被同化模塊刪除后, 輸入再分析的海浪譜, 海浪模式繼續(xù)向前積分, 至此完成一次海浪資料同化過程。同化模塊啟動(dòng)后, 讀入控制文件和海浪譜, 判斷同化方法的類型, 若為靜態(tài)樣本集合Kalman濾波方案, 則利用靜態(tài)樣本集合構(gòu)造模式狀態(tài)變量集合; 若為EAKF方案, 則直接利用多個(gè)模式的變量構(gòu)造模式狀態(tài)變量集合。隨后同化模塊利用兩步濾波方法, 對(duì)模式狀態(tài)變量集合進(jìn)行更新, 輸出同化后的海浪譜并刪除控制文件。

海浪模式運(yùn)行過程中對(duì)二維海浪譜進(jìn)行積分計(jì)算, 有效波高為海浪譜的積分量。同化過程中, 本文采用 Esteva(1988)的方法對(duì)海浪譜進(jìn)行調(diào)整, 即僅對(duì)波譜總能量進(jìn)行比例系數(shù)縮放, 不改變譜形, 如式(13)所示。

其中, 上標(biāo)a和b分別表示同行后和同化前, 下標(biāo)i和j分別表示經(jīng)度和緯度方向的空間網(wǎng)格, E(f, θ)表示二維波譜, f和θ分別表示頻率和方向, H表示有效波高。

圖1 海浪同化系統(tǒng)運(yùn)行流程Fig.1 Flowchart of wave data assimilation system

3 同化實(shí)驗(yàn)結(jié)果與分析

利用Saral衛(wèi)星高度計(jì)資料檢驗(yàn)海浪數(shù)據(jù)同化實(shí)驗(yàn)結(jié)果。由于衛(wèi)星高度計(jì)觀測(cè)為沿軌數(shù)據(jù), 海浪模式輸出結(jié)果為網(wǎng)格數(shù)據(jù), 對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行誤差統(tǒng)計(jì)前,本文利用雙線性插值方法, 對(duì)模擬結(jié)果和觀測(cè)資料進(jìn)行了時(shí)空匹配處理。實(shí)驗(yàn)結(jié)果檢驗(yàn)分為兩個(gè)部分,整體誤差統(tǒng)計(jì)和誤差時(shí)空分布特征分析。

3.1 整體誤差統(tǒng)計(jì)

為檢驗(yàn)同化實(shí)驗(yàn)的整體效果, 采用絕均差和均方根誤差公式, 對(duì)有效波高模擬誤差進(jìn)行分區(qū)統(tǒng)計(jì)分析。根據(jù)通常意義的地理緯度劃分標(biāo)準(zhǔn)及風(fēng)浪和涌浪全球空間分布特征(Chen et al, 2002), 考慮了高緯度地區(qū)海冰的季節(jié)性變化和衛(wèi)星高度計(jì)的實(shí)際觀測(cè)范圍, 劃分誤差統(tǒng)計(jì)分區(qū)如圖2所示, 統(tǒng)計(jì)結(jié)果見表1和表 2。從表 1中全球統(tǒng)計(jì)結(jié)果來看, 與控制實(shí)驗(yàn)相比, 兩組同化實(shí)驗(yàn)的有效波高模擬絕均差和均方根誤差減少了約 25%和 20%, 靜態(tài)樣本集合 Kalman濾波方法和EAKF方法的同化效果顯著, 兩者同化效果相差不大。在統(tǒng)計(jì)分區(qū)Ⅰ和 Ⅲ , EAKF 方案的有效波高模擬誤差明顯小于靜態(tài)樣本集合Kalman濾波方案; 在統(tǒng)計(jì)分區(qū)Ⅱ, 情況相反, 后者表現(xiàn)更佳。

為了檢驗(yàn)?zāi)M與觀測(cè)的線性相關(guān)程度, 繪制模擬-觀測(cè)空間散點(diǎn)圖。具體做法為: 將有效波高的模擬-觀測(cè)空間劃分為分辨率 0.05m×0.05m 的網(wǎng)格, 統(tǒng)計(jì)落在每個(gè)網(wǎng)格內(nèi)的觀測(cè)點(diǎn)數(shù)目, 統(tǒng)計(jì)結(jié)果如圖 3所示。散點(diǎn)分布越集中于黑色對(duì)角虛線, 說明模擬與觀測(cè)的線性相關(guān)程度越高, 與控制實(shí)驗(yàn)相比, 兩組同化實(shí)驗(yàn)的散點(diǎn)分布更集中于對(duì)角線。與控制實(shí)驗(yàn)相比,靜態(tài)樣本集合Kalman濾波和EAKF方案的模擬-觀測(cè)相關(guān)系數(shù)提高了約4.7%和4.9%。

圖2 有效波高模擬誤差統(tǒng)計(jì)分區(qū)示意圖Fig.2 Schematic diagram of statistical regions for SWH simulation error

表1 2014年有效波高模擬絕均差統(tǒng)計(jì)結(jié)果Tab. 1 Mean absolute error comparing simulated SWH with observations from satellite Saral

表2 2014年有效波高模擬均方根誤差統(tǒng)計(jì)結(jié)果Tab.2 Root mean square error comparing simulated SWH with observations from satellite Saral

此外, 從圖 3中可以看出, 當(dāng)有效波高較小時(shí)(4m 以下), 模擬較觀測(cè)偏大, 當(dāng)有效波高較大時(shí)(4m以上), 模擬較觀測(cè)偏小。為了進(jìn)一步了解兩組同化實(shí)驗(yàn)對(duì)于不同區(qū)間內(nèi)有效波高的修正情況, 將有效波高劃分三個(gè)區(qū)間, 0.5—2m、2—4m和4m以上, 針對(duì)上述三個(gè)區(qū)間, 分別統(tǒng)計(jì)模擬誤差的概率密度分布,如圖4所示。

圖3 2014年全球有效波高模擬-觀測(cè)空間散點(diǎn)圖Fig.3 Scatter diagrams comparing simulated SWH with global observations in 2014

圖4 2014年全球有效波高模擬誤差概率密度分布Fig.4 Probability density distribution of global SWH simulation error in 2014

圖 4中, 曲線峰值越接近于 0, 說明模擬誤差的均值越小, 曲線的聚攏程度表示模擬誤差相對(duì)于其均值的集中程度。當(dāng)有效波高觀測(cè)值在0.5—2m區(qū)間內(nèi), 模擬較觀測(cè)偏大, 靜態(tài)樣本集合 Kalman濾波方案對(duì)控制實(shí)驗(yàn)的改善優(yōu)于EAKF方案; 當(dāng)有效波高觀測(cè)值在2—4m區(qū)間內(nèi), 模擬與觀測(cè)相差不大, 兩組同化方案對(duì)控制實(shí)驗(yàn)的改善程度相當(dāng); 當(dāng)有效波高觀測(cè)值大于4m時(shí), 模擬較觀測(cè)偏小, EAKF方案對(duì)控制實(shí)驗(yàn)的修正優(yōu)于靜態(tài)樣本集合Kalman濾波方案。

通過上述整體統(tǒng)計(jì)分析, 可以得到以下結(jié)論, 靜態(tài)樣本集合Kalman濾波和EAKF方案同化效果顯著,總體來講, 兩者相差不大; 在波高較小的低緯度地區(qū),前者占優(yōu); 在波高較大的中高緯度地區(qū), 后者占優(yōu)。

3.2 誤差時(shí)空分布特征分析

為進(jìn)一步檢驗(yàn)和評(píng)估靜態(tài)樣本集合Kalman濾波和EAKF方案的同化效果, 本節(jié)將對(duì)有效波高模擬絕均差的空間分布和時(shí)間序列進(jìn)行統(tǒng)計(jì), 分析兩組同化實(shí)驗(yàn)的誤差時(shí)空分布特征。

(1) 誤差空間分布

為了解兩組同化實(shí)驗(yàn)的有效波高模擬誤差空間分布特征, 分別統(tǒng)計(jì)控制實(shí)驗(yàn)和兩組同化實(shí)驗(yàn)的有效波高模擬絕均差空間分布情況。與控制實(shí)驗(yàn)相比,靜態(tài)樣本集合Kalman濾波和EAKF同化實(shí)驗(yàn)的有效波高模擬絕均差減小量的空間分布如圖5所示。從圖中可以看出, 兩組同化實(shí)驗(yàn)效果顯著, 且相差不大。靜態(tài)樣本集合具有背景誤差的特征信息。靜態(tài)樣本集合采用24h間隔有效波高模擬偏差近似背景誤差, 單個(gè)時(shí)刻僅存在一個(gè)誤差樣本, 為獲取足夠多的樣本,該方法對(duì)有效波高歷史模擬結(jié)果進(jìn)行長時(shí)間序列采樣, 由此, 靜態(tài)樣本集合包含了每個(gè)空間位置在某時(shí)段內(nèi)的背景誤差的特征信息, 這也是靜態(tài)樣本集合Kalman濾波方法的優(yōu)勢(shì)之一。

圖5 2014年有效波高模擬絕均差同化減小量的空間分布(m)Fig.5 Spatial distribution of MAE reduction of simulated SWH with assimilation in 2014 (unit: m)

此外, 從圖5可以看出, 兩組同化實(shí)驗(yàn)有效波高模擬絕均差減小量的空間分布結(jié)構(gòu)類似, 低緯度地區(qū)的同化效果優(yōu)于中高緯度地區(qū), 這與涌浪和風(fēng)浪在某空間位置的比例有關(guān)。大洋東岸的熱帶和副熱帶海域, 存在明顯的涌浪區(qū), 特別是太平洋和大西洋,涌浪區(qū)顯著(Chen et al, 2002)。海浪同化過程中, 對(duì)波譜進(jìn)行了整體縮放。風(fēng)浪受風(fēng)強(qiáng)迫的影響, 變化較快,同化信息難以長時(shí)間保留; 與之相反, 涌浪則能夠較好的保留海浪同化信息, 因此低緯度地區(qū)的同化效果優(yōu)于中高緯度地區(qū)。

從有效波高模擬絕均差同化減小量的空間分布情況來看, 靜態(tài)樣本集合Kalman濾波和EAKF同化效果相近, 但仍然存在一定的差異。為更直觀地分析兩組同化方案效果的差異, 將兩者的有效波高模擬絕均差做差, 后者減前者, 得到兩者同化效果差異的空間分布, 如圖6所示。

圖6 2014年靜態(tài)樣本集合Kalman濾波和EAKF有效波高模擬絕均差的差異的空間分布(m)Fig.6 Spatial distribution of difference of MAE of simulated SWH between static sample ensemble Kalman filter and EAKF(2014, units: m)

圖 6中, 藍(lán)色代表 EAKF同化方案的優(yōu)勢(shì)區(qū)域,紅色代表靜態(tài)樣本集合Kalman濾波同化方案的優(yōu)勢(shì)區(qū)域, 白色代表兩者同化效果相當(dāng)?shù)膮^(qū)域。從圖中可以看出, 在中高緯度地區(qū), EAKF方案表現(xiàn)較好, 在低緯度地區(qū), 靜態(tài)樣本集合 Kalman濾波方案表現(xiàn)較好, 兩者的差異幅度約為 0.02m。EAKF方案具有傳統(tǒng)集合濾波方法的天然優(yōu)勢(shì), 包含背景誤差的瞬時(shí)信息, 在風(fēng)場(chǎng)變化較大的中高緯度地區(qū), 特別是西風(fēng)帶區(qū)域, 其優(yōu)勢(shì)顯著。在EAKF方案中, 風(fēng)場(chǎng)擾動(dòng)對(duì)風(fēng)浪的影響較大, 對(duì)涌浪的影響則較小; 在靜態(tài)樣本集合 Kalman濾波方案中, 靜態(tài)樣本的構(gòu)造方式?jīng)Q定了其能夠包含更多的涌浪信息; 由此, 前者的空間相關(guān)系數(shù)較后者相對(duì)局地化, 該情況在涌浪占主導(dǎo)的低緯度地區(qū)更加凸顯, 這也是后者同化效果在低緯度區(qū)域略優(yōu)于前者的原因。

(2) 誤差時(shí)間序列

為了解兩組同化實(shí)驗(yàn)的有效波高模擬誤差時(shí)間分布特征, 統(tǒng)計(jì)控制實(shí)驗(yàn)和兩組同化實(shí)驗(yàn)的有效波高模擬絕均差隨時(shí)間的變化情況。對(duì)每個(gè)誤差統(tǒng)計(jì)分區(qū)(如圖 2所示), 逐日統(tǒng)計(jì)有效波高模擬絕均差,得到模擬誤差時(shí)間序列如圖 7所示。從圖 7(b)(c)和(d)藍(lán)色曲線變化情況可以看出, 控制實(shí)驗(yàn)在低緯度地區(qū)(30°S—30°N)及中高緯度地區(qū)(30°N—60°N, 30°S—60°S)的夏季, 有效波高模擬絕均差在0.6m上下浮動(dòng), 在其他區(qū)域和時(shí)段內(nèi), 模擬誤差位于 0.8m 上下, 即風(fēng)場(chǎng)變化較小時(shí), 海浪模式的模擬效果更好。從圖7(b)(c)和(d)綠色和紅色曲線變化情況可以看出, 在風(fēng)場(chǎng)變化較小的區(qū)域和時(shí)段內(nèi),靜態(tài)樣本集合Kalman濾波和EAKF方案的同化效果顯著且兩者差異不大; 兩者的差異主要體現(xiàn)在中高緯度地區(qū)的冬季, 后者具有明顯優(yōu)勢(shì)。從圖 7(a)全球范圍內(nèi)的模擬誤差統(tǒng)計(jì)結(jié)果來看, 兩組同化實(shí)驗(yàn)的有效波高模擬絕均差變化曲線幾乎重合, 基本穩(wěn)定在0.5m上下。

圖7 2014年有效波高模擬絕均差隨時(shí)間變化曲線(m)Fig.7 Curves of MAE of simulated SWH in 2014 (unit: m)

為更直觀地分析靜態(tài)樣本集合 Kalman濾波和EAKF方案同化效果的差異, 在不同統(tǒng)計(jì)分區(qū)內(nèi), 將兩者的有效波高模擬絕均差逐日做差, 后者減前者,得到兩者差異隨時(shí)間的變化情況, 如圖8所示。

圖8 2014年靜態(tài)樣本集合Kalman濾波和EAKF有效波高模擬絕均差的差異隨時(shí)間變化曲線(m)Fig.8 Curves of difference of MAE of simulated SWH between static sample ensemble Kalman filter and EAKF in 2014 (unit: m)

在圖8中, 虛線(零線)以下代表EAKF方案占優(yōu),虛線以上代表靜態(tài)樣本集合Kalman濾波方案占優(yōu)。從全球統(tǒng)計(jì)結(jié)果(綠色曲線)來看, 兩種方案的同化效果相當(dāng), 兩者差異在零線上下浮動(dòng), 幅度約為0.02m。在低緯度地區(qū)(30°S—30°N), 靜態(tài)樣本集合Kalman濾波同化方案優(yōu)勢(shì)明顯, 兩者差異在 0.03m上下浮動(dòng); 在中高緯度地區(qū)(30°N—60°N, 30°S—60°S), 風(fēng)速較大, 且風(fēng)場(chǎng)變化較快, EAKF同化方案優(yōu)勢(shì)明顯。從圖中藍(lán)色曲線(30°N—60°N)的變化情況可以看出, 北半球的冬季, 風(fēng)場(chǎng)變化較大, EAKF方案優(yōu)勢(shì)顯著; 北半球的夏季, 風(fēng)場(chǎng)變化相對(duì)冬季更為平穩(wěn),靜態(tài)樣本集合Kalman濾波方案占優(yōu)。從圖中紅色曲線(30°S—60°S)的變化情況, 可以得到上述類似的結(jié)論。

4 結(jié)論

背景誤差是影響海浪同化效果的關(guān)鍵因素之一,最優(yōu)插值和三維變分方法僅對(duì)背景誤差進(jìn)行一次估計(jì), 即“靜態(tài)”統(tǒng)計(jì), 集合 Kalman濾波方法可以對(duì)背景誤差進(jìn)行“動(dòng)態(tài)”統(tǒng)計(jì), 實(shí)時(shí)更新背景誤差, 是更為合理的方法。但集合Kalman濾波方法對(duì)計(jì)算資源的要求較高, 在實(shí)際應(yīng)用中難以實(shí)現(xiàn)。孫盟等(2014)提出基于靜態(tài)樣本的集合 Kalman濾波同化方法, 該方法利用24h有效波高模擬偏差構(gòu)造靜態(tài)樣本集合, 用于近似背景誤差, 在同化時(shí)刻, 將靜態(tài)樣本集合疊加到有效波高場(chǎng), 實(shí)現(xiàn)模式狀態(tài)變量集合的構(gòu)造, 該方法僅需運(yùn)行一個(gè)海浪模式, 可大大降低對(duì)計(jì)算資源的需求。為量化靜態(tài)樣本集合Kalman濾波同化方案與傳統(tǒng)集合 Kalman濾波同化方案的差異, 本文基于MASNUM-WAM 海浪模式, 采用靜態(tài)樣本集合Kalman濾波和EAKF方法, 利用Jason-2衛(wèi)星高度計(jì)資料, 針對(duì) 2014年全球海域開展了海浪資料同化實(shí)驗(yàn)。EAKF同化方案需要構(gòu)造具有代表性的海浪模式集合, 本文采用隨機(jī)場(chǎng)擾動(dòng)方案(孫盟等, 2016), 對(duì)風(fēng)場(chǎng)進(jìn)行集合擾動(dòng), 再由風(fēng)場(chǎng)集合驅(qū)動(dòng)生成海浪模式集合。同化實(shí)驗(yàn)的檢驗(yàn)資料為Saral衛(wèi)星高度計(jì)資料, 結(jié)果表明, 兩組同化實(shí)驗(yàn)的效果顯著, 有效波高模擬絕均差和均方根誤差較未加入同化的控制實(shí)驗(yàn)減小了約25%和20%。根據(jù)2014年全球的誤差統(tǒng)計(jì)結(jié)果, 靜態(tài)樣本集合Kalman濾波和EAKF方案的同化效果相差不大, 但兩者仍然存在一定差異。在風(fēng)場(chǎng)變化較小的低緯度地區(qū),靜態(tài)樣本集合Kalman濾波方案較EAKF方案略優(yōu); 在風(fēng)場(chǎng)變化較大的中高緯度地區(qū), EAKF方案優(yōu)勢(shì)顯著。

靜態(tài)樣本集合Kalman濾波同化方案在同化過程中, 僅需運(yùn)行一個(gè)海浪模式, 大大降低了計(jì)算成本,雖然其在中高緯度地區(qū)的表現(xiàn)不如 EAKF同化方案,但從整體統(tǒng)計(jì)結(jié)果來看, 靜態(tài)樣本集合 Kalman濾波方案的同化效果已經(jīng)十分接近EAKF同化方案的同化效果, 是一種性價(jià)比較高的海浪數(shù)據(jù)同化方案,可應(yīng)用于海浪實(shí)時(shí)業(yè)務(wù)化預(yù)報(bào)和歷史再分析數(shù)據(jù)的重構(gòu)。

孫 盟, 尹訓(xùn)強(qiáng), 楊永增, 2014. 靜態(tài)集合樣本的構(gòu)造及其在全球海浪濾波同化中的應(yīng)用. 海洋與湖沼, 45(5):918—927

孫 盟, 尹訓(xùn)強(qiáng), 楊永增等, 2016. MASNUM-WAM海浪模式集合Kalman濾波同化研究—I. 風(fēng)場(chǎng)擾動(dòng)對(duì)海浪模擬影響.海洋與湖沼, 47(6): 1091—1100

楊永增, 喬方利, 趙 偉等, 2005. 球坐標(biāo)系下 MASNUM 海浪數(shù)值模式的建立及其應(yīng)用. 海洋學(xué)報(bào), 27(2): 1—7

Anderson J L, 2001. An ensemble adjustment Kalman filter for data assimilation. Monthly Weather Review, 129(12):2884—2903

Anderson J L, 2003. A local least squares framework for ensemble filtering. Monthly Weather Review, 131(4): 634—642

Burgers G, Jan van Leeuwen P, Evensen G, 1998. Analysis scheme in the ensemble Kalman filter. Monthly Weather Review, 126(6): 1719—1724

Chen G, Chapron B, Ezraty R et al, 2002. A global view of swell and wind sea climate in the ocean by satellite altimeter and scatterometer. Journal of Atmospheric and Oceanic Technology, 19(11): 1849—1859

Chen H, Yin X Q, Bao Y et al, 2016. Ocean satellite data assimilation experiments in FIO—ESM using ensemble adjustment Kalman filter. Science China Earth Sciences,59(3): 484—494

Esteva D C, 1988. Evaluation of preliminary experiments assimilating Seasat significant wave heights into a spectral wave model. Journal of Geophysical Research, 93(C11):14099—14105

Evensen G, 1994. Sequential data assimilation with a nonlinear quasi—geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research:Oceans, 99(C5): 10143—10162

Evensen G, 2003. The Ensemble Kalman Filter: theoretical formulation and practical implementation. Ocean Dynamics,53(4): 343—367

Hill A J, Weiss C C, Ancell B C, 2016. Ensemble sensitivity

analysis for mesoscale forecasts of dryline convection initiation. Monthly Weather Review, 144(11): 4161—4182 Kalman R E, 1960. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1):35—45

Kalman R E, Bucy R S, 1961. New results in linear filtering and prediction theory. Journal of Basic Engineering, 83(1):95—108

Karspeck A R, Anderson J L, 2007. Experimental implementation of an ensemble adjustment filter for an intermediate ENSO model. Journal of Climate, 20(18): 4638—4658

Karspeck A R, Yeager S, Danabasoglu G et al, 2013. An ensemble adjustment Kalman filter for the CCSM4 ocean component. Journal of Climate, 26(19): 7392—7413

Lorenc A C, 1981. A global three-dimensional multivariate statistical interpolation scheme. Monthly Weather Review,109(4): 701—721

Shen Z Q, Zhang X M, Tang Y M, 2016. Comparison and combination of EAKF and SIR-PF in the Bayesian filter framework. Acta Oceanologica Sinica, 35(3): 69—78

Yin X Q, Qiao F L, Shu Q, 2011. Using ensemble adjustment Kalman filter to assimilate Argo profiles in a global OGCM.Ocean Dynamics, 61(7): 1017—1031

Yin X Q, Qiao F L, Yang Y Z et al, 2010. An ensemble adjustment Kalman filter study for Argo data. Chinese Journal of Oceanology and Limnology, 28(3): 626—635

Yuan Y L, Hua F, Pan Z D et al, 1991. LAGFD-WAM numerical wave model——I. basic physical model. Acta Oceanologica Sinica, 10(4): 483—488

Yuan Y L, Hua F, Pan Z D et al, 1992. LAGFD-WAM numerical wave model—Ⅱ. Characteristics inlaid scheme and its application. Acta Oceanologica Sinica, 11(1): 13—23

Yuan Y L, Tung C C, Huang N E, 1986. Statistical characteristics of breaking waves. In: Phillips O M, Hasselmann K eds.Wave Dynamics and Radio Probing of the Ocean Surface.New York: Springer, 265—272

Zhang S, Harrison M J, Rosati A et al, 2007. System design and evaluation of coupled ensemble data assimilation for global oceanic climate studies. Monthly Weather Review, 135(10):3541—3564

Zhang S, Harrison M J, Wittenberg A T et al, 2005. Initialization of an ENSO forecast system using a parallelized ensemble filter. Monthly Weather Review, 133(11): 3176—3201