基于海浪譜分解與重構(gòu)的資料同化試驗(yàn)

毛科峰,蕭中樂(lè),宋海波,鐘奕飛

（1 .解放軍理工大學(xué)氣象海洋學(xué)院,江蘇南京211101；2 .解放軍96631部隊(duì),北京102208；3 .解放軍92721部隊(duì),浙江舟山316000；4 .解放軍95857部隊(duì),湖南衡陽(yáng)421000）

數(shù)值計(jì)算時(shí)表達(dá)式為:

?

基于海浪譜分解與重構(gòu)的資料同化試驗(yàn)

毛科峰1,蕭中樂(lè)2,宋海波3,鐘奕飛4

（1 .解放軍理工大學(xué)氣象海洋學(xué)院,江蘇南京211101；2 .解放軍96631部隊(duì),北京102208；3 .解放軍92721部隊(duì),浙江舟山316000；4 .解放軍95857部隊(duì),湖南衡陽(yáng)421000）

摘要:針對(duì)有效波高資料提出一種海浪譜分解與重構(gòu)的資料同化方案:利用歷史時(shí)段內(nèi)的有效波高觀測(cè)資料和模式計(jì)算波高場(chǎng),采用最優(yōu)插值方法得到分析波高場(chǎng)；在WAV E WAT C H -Ⅲ模式的波浪能量密度譜和有效波高分析值之間引入一個(gè)變異系數(shù)矩陣,描述模式的誤差,以此為狀態(tài)向量構(gòu)建卡爾曼濾波系統(tǒng),對(duì)分解過(guò)的海浪譜進(jìn)行修正和重構(gòu),得到同化后的海浪譜初始場(chǎng)。利用美國(guó)阿拉斯加灣北部海域的7個(gè)浮標(biāo)站進(jìn)行同化和72 h預(yù)報(bào)試驗(yàn),對(duì)連續(xù)1個(gè)月的預(yù)報(bào)結(jié)果進(jìn)行統(tǒng)計(jì)表明:采用該同化方案后24 h預(yù)報(bào)結(jié)果的有效波高均方根誤差比未同化的結(jié)果降低了0.13 m；同化方案對(duì)預(yù)報(bào)效果的影響可持續(xù)36 h左右,隨著預(yù)報(bào)時(shí)效延長(zhǎng),同化的效果減弱。

關(guān)鍵詞:海浪譜；有效波高；同化；WAV E WAT C H -Ⅲ模式；卡爾曼濾波

1 引言

為了提高海浪模式預(yù)測(cè)能力,將模式結(jié)果和實(shí)時(shí)資料進(jìn)行同化,能夠改善模式預(yù)報(bào)的準(zhǔn)確度［1］。目前同化的海浪資料包括實(shí)測(cè)波高和實(shí)測(cè)海浪譜,由于海浪模式的控制方程是建立在波譜能量平衡方程基礎(chǔ)上,觀測(cè)得到的有效波高對(duì)應(yīng)著波譜能量平衡方程中的波浪能譜的積分,如何用有效波高去同化波浪能譜是資料同化時(shí)面臨的一個(gè)重要問(wèn)題。Komen［2］首先嘗試將實(shí)測(cè)波高應(yīng)用于波浪數(shù)值模式資料同化的研究,之后Hasselmann等［3］和Bauer等［4］以實(shí)測(cè)波高與模式計(jì)算值之比,作為修正二維波譜的參數(shù)去調(diào)整海浪譜。這種方法相對(duì)簡(jiǎn)單,但是對(duì)波浪譜的調(diào)整,比較“粗糙”,沒(méi)有考慮實(shí)測(cè)資料對(duì)波浪譜方向和頻率分布的影響,后來(lái)Hasselmann等［5］在波高資料同化研究中改善了波譜修正的尺度因子,也有學(xué)者提出將波浪譜分離成風(fēng)浪和涌浪部分分別同化。另外一方面,在同化高度計(jì)波高資料時(shí)采用考慮空間效應(yīng)的修正因子去修正波浪譜,Greenslade和Young［6］采用最優(yōu)插值方法同化E RS-2高度計(jì)資料,研究了觀測(cè)誤差方差和背景誤差方差的不同模型和背景誤差相關(guān)尺度隨緯度變化等因素導(dǎo)致的同化效果。E m manouil等［7］發(fā)展了一種基于集合卡爾曼濾波方法的波浪資料同化方法,能夠結(jié)合誤差協(xié)方差的空間分布特點(diǎn),使得24 h內(nèi),波高的均方根誤差減小30 %～60 %。國(guó)內(nèi)的海浪資料同化研究起步較晚,張志旭等［8］通過(guò)在臺(tái)風(fēng)浪預(yù)報(bào)中同化高度計(jì)資料后改進(jìn)預(yù)報(bào)效果。王毅［9］、王毅和余宙文［10］用最優(yōu)插值方法將衛(wèi)星高度計(jì)的波高資料同化到S WA N模式中,并通過(guò)業(yè)務(wù)預(yù)報(bào)檢驗(yàn)了效果,表明預(yù)報(bào)效果有一定程度改善。Eeng等［11］采用最優(yōu)插值方法在S WA N模式中首次同化了臺(tái)灣島臨近海域的實(shí)測(cè)資料用于臺(tái)風(fēng)浪的預(yù)報(bào)研究,發(fā)現(xiàn)資料同化效果的影響范圍取決于臺(tái)風(fēng)發(fā)展?fàn)顟B(tài)及涌浪的傳播方向。Guo等［12］利用蒙特卡羅方法模擬計(jì)算有效波高的背景誤差協(xié)方差,并發(fā)展了一種新的基于平均波向的有效波高的背景誤差協(xié)方差模型,用于改進(jìn)波高資料的同化效果。

目前,利用有效波高資料對(duì)波浪模式同化的方法研究是為了將觀測(cè)波高的信息在物理意義上和空間相關(guān)性上去修正模式的波浪能譜。其中風(fēng)浪和涌浪分離的方案在物理意義上具有優(yōu)勢(shì),然而涌浪的判據(jù)和分離本身也是一個(gè)較復(fù)雜的問(wèn)題；通過(guò)修正尺度因子去調(diào)整波浪譜取得了較好的效果,但方法也依賴于經(jīng)驗(yàn)。因此,本文基于海浪模式中離散的海浪能量密度譜與有效波高的線性關(guān)系式,提出一種基于海浪譜分解與重構(gòu)的波高資料同化方法,并在數(shù)值模式中進(jìn)行預(yù)報(bào)效果試驗(yàn)分析。

2 有效波高資料的同化方案

2.1 波浪能量密度譜變異系數(shù)的引入

海浪模式求解海浪能量密度譜平衡方程時(shí),預(yù)報(bào)出計(jì)算格點(diǎn)上,每個(gè)預(yù)報(bào)時(shí)次的波浪能量密度譜E（f,θ）,然后根據(jù)波高、周期的統(tǒng)計(jì)分布關(guān)系式,計(jì)算出波浪要素,如平均波高、周期等,其中,有效波高的計(jì)算式為:

數(shù)值計(jì)算時(shí)表達(dá)式為:

式中,Ei,j是頻率為fi、方向?yàn)棣萰時(shí)離散的波浪能量密度譜值,方向和頻率離散的維數(shù)是nf和nθ,它表示波浪能量不僅分布在一定的頻率范圍內(nèi),而且分布在相當(dāng)寬的方向范圍內(nèi),這描述了海浪內(nèi)部結(jié)構(gòu),也能間接地表征海浪外在表現(xiàn)特征。值得注意的是式（1）被離散成式（2）后,與Ei,j構(gòu)成了一個(gè)線性關(guān)系式。為了便于區(qū)分,本文用表示模式計(jì)算的有效波高Hs的真實(shí)值,數(shù)值模式計(jì)算得到的離散波浪能量密度譜值記為,為了描述模式計(jì)算的波浪能量密度譜值和真實(shí)有效波高之間的誤差和關(guān)系,引入一個(gè)模式計(jì)算的波浪能量密度譜和有效波高真值之間的變異系數(shù)矩陣,記為Bi,j其維數(shù)與相同,則有:

2.2 資料同化的思路

本文的資料同化思路就是基于這個(gè)變異系數(shù),描述模式計(jì)算的能量密度譜和實(shí)測(cè)波高之間的誤差,以此為狀態(tài)向量構(gòu)建卡爾曼濾波系統(tǒng),對(duì)分解過(guò)的海浪譜進(jìn)行修正和重構(gòu)。為此首先對(duì)模式計(jì)算得到的較長(zhǎng)一段時(shí)間內(nèi)的波浪能量密度譜EMi,j進(jìn)行經(jīng)驗(yàn)正交分解,記ek為特征值,對(duì)分解出的特征向量進(jìn)行顯著性檢驗(yàn),并判定其物理意義,當(dāng)式（4）成立時(shí),

ek和ek+ 1是可以分離的,特征值是有意義的。然后將特征值按從大到小的順序進(jìn)行排序,這里不妨設(shè)下式成立:

選取特征值較大前nk階特征向量來(lái)表示實(shí)際的。對(duì)t時(shí)刻模式得到的波浪能量密度譜Ei,j|t= t簡(jiǎn)記為Et,它的維數(shù)為（nfnθ）,分解成:

ck,t和ek,t是t時(shí)刻的k階的特征向量和特征值,也有:

則式（2）可寫為:

式中,Et=（e1,t…ek,t…enk,t）,

因此對(duì)于式（3）,可以寫成:

式中,βt是前nk階特征向量對(duì)應(yīng)的波浪能量密度譜變異系數(shù)矩陣,它的維數(shù)為（nk×1）。這一線性表達(dá)式可以被視為卡爾曼濾波中的量測(cè)方程,引入量測(cè)噪音εt,它是一維隨機(jī)向量,則寫為:

將波浪能量密度譜變異系數(shù)矩陣βt視為卡爾曼濾波系統(tǒng)中的狀態(tài)向量,用狀態(tài)方程式（12）描述其變化,和量測(cè)方程一起構(gòu)成卡爾曼濾波系統(tǒng),

δt- 1是動(dòng)態(tài)噪音,與量測(cè)噪音εt一樣都是隨機(jī)向量,假設(shè)這兩個(gè)量滿足互不相關(guān),均值為0,方差分別為W 和V ,因此構(gòu)成如下遞推系統(tǒng):

為了簡(jiǎn)化表達(dá)式,令Xt= DEtCt,則:

利用初始時(shí)刻t以前的歷史資料通過(guò)上述遞推系統(tǒng)得到t時(shí)刻的變異系數(shù)矩陣βt后,就可以利用式（6）重構(gòu)t時(shí)刻的波浪譜,記Ea為同化后的波浪譜,則有:

2.3 有效波高的最優(yōu)插值方法

式中,Nobs為影響該模式格點(diǎn)的實(shí)測(cè)浮標(biāo)站點(diǎn)的數(shù)目,Wij則為各觀測(cè)站點(diǎn)相對(duì)于計(jì)算格點(diǎn)的最佳權(quán)重,下標(biāo)i表示計(jì)算格點(diǎn),下標(biāo)j表示浮標(biāo)觀測(cè)站。由于實(shí)測(cè)波高值HO、模式的起始猜測(cè)值HP及分析值Ha與真值HT之間有誤差存在,要求取各測(cè)站對(duì)于計(jì)算格點(diǎn)的最佳權(quán)重,必須求解Nobs組聯(lián)立方程組,若僅有少量的觀測(cè)點(diǎn),則用高斯消元法可求解。故:

若已知真值的條件下,易求得實(shí)測(cè)值與模式值的均方誤差,但實(shí)際上,真值未知,故采用已有研究成果中的處理方法［6］,觀測(cè)誤差協(xié)方差取為Okj=,其中δkj為kronecher符號(hào),R為實(shí)測(cè)均方誤差與模式均方誤差之比。Pkj是觀測(cè)點(diǎn)k 和j上的模式計(jì)算波高的誤差協(xié)方差,Pik是計(jì)算格點(diǎn)i和觀測(cè)點(diǎn)k的模式計(jì)算波高誤差協(xié)方差。采用Greenslade和Young［6］在海浪同化研究中的方案,波高背景場(chǎng)誤差協(xié)方差取為:

式中,dkj是觀測(cè)點(diǎn)k和j之間的距離；σk和σj是觀測(cè)點(diǎn)k,j上的有效波高背景場(chǎng)的均方根誤差；L為影響尺度半徑；α、β、γ都是經(jīng)驗(yàn)常數(shù)。

2.4 資料同化的計(jì)算步驟

首先選擇同化某時(shí)刻t之前連續(xù)m天的模式計(jì)算波高場(chǎng)和方向譜矩陣序列,利用相應(yīng)時(shí)刻的波高觀測(cè)資料通過(guò)最優(yōu)插值方法得到分析波高場(chǎng),對(duì)方向譜序列進(jìn)行第一次正交經(jīng)驗(yàn)分解。選擇前nk個(gè)主成分,記錄前nk個(gè)特征向量；然后選取同樣時(shí)間長(zhǎng)度的方向譜樣本,用第一次分解得到的nk個(gè)空間函數(shù),將之進(jìn)行第二次正交經(jīng)驗(yàn)分解,選取同樣個(gè)數(shù)的特征向量。利用第一次分解得到的nk個(gè)空間函數(shù)和分析波高,進(jìn)行線性擬合,得到第一組波浪能量密度譜傳遞系數(shù)值；利用第二次正交經(jīng)驗(yàn)分解得到的nk空間函數(shù)和預(yù)報(bào)點(diǎn)的浪高,進(jìn)行線性擬合,得到第二組波浪能量密度譜傳遞系數(shù)值。記C0為的誤差方差陣,認(rèn)為是從樣本資料直接計(jì)算得到,與理論值相等,故C0取為nk階零方陣。W矩陣主對(duì)角線上的值是:Δ/Δt,i = 1,2,3,4 ,其他元素為0,其中,Δ、Δ分別是第一組波浪能量密度譜傳遞系數(shù)值與第二組波浪能量密度譜傳遞系數(shù)對(duì)應(yīng)的差的平方, Δt分母為分解所用數(shù)據(jù)的時(shí)間樣本個(gè)數(shù)。由于量測(cè)量只有Hs一個(gè)量,故V為1×1矩陣,數(shù)值為V = q/（k - m - 1）,q為第一組正交分解后線性回歸的殘差,預(yù)報(bào)與實(shí)際的差的平方和,k為取得樣本個(gè)數(shù),即數(shù)據(jù)的時(shí)間個(gè)數(shù),m為選取主成分的個(gè)數(shù)。確定上述參數(shù)的取值后,可以利用卡爾曼遞推系統(tǒng)得到t時(shí)刻的變異系數(shù)矩陣βt,然后重構(gòu)t時(shí)刻的波浪譜,得到同化后的初始場(chǎng)。

3 同化預(yù)報(bào)試驗(yàn)

3.1 資料和海浪模式

為了開(kāi)展同化預(yù)報(bào)的數(shù)值試驗(yàn),將波浪計(jì)算的區(qū)域選為美國(guó)阿拉斯加灣北部海域,一個(gè)重要的原因是因?yàn)樵搮^(qū)域能夠獲得相當(dāng)數(shù)量的波浪實(shí)測(cè)資料。

（1）同化的浮標(biāo)資料

浮標(biāo)資料由美國(guó)國(guó)家數(shù)據(jù)浮標(biāo)中心（N DBC）提供,浮標(biāo)主要集中在美國(guó)阿拉斯加灣北部海域,地理位置如圖1示意,浮標(biāo)的編號(hào)和該點(diǎn)的水深情況如表1,采用的是“波浪騎士”浮標(biāo),采樣頻率為1.28 Hz,浮標(biāo)實(shí)測(cè)有效波高資料的間隔為1 h。

圖1 同化資料的浮標(biāo)位置Eig .1 The buoys used for assimilation

表1 浮標(biāo)的位置和水深Tab .1 Buoy coordinates and water depth

（2）WAV E WAT C H -Ⅲ模式

本文的研究基于WAV E WAT C H -Ⅲ第三代海浪數(shù)值模式開(kāi)發(fā)了新的資料同化模塊,并進(jìn)行數(shù)值試驗(yàn),模式基于波浪能譜平衡方程:

式中,N即是波作用密度譜,它是頻率f、傳播方向θ、時(shí)間t和地理空間位置的函數(shù),cg為群速,S是能量源匯,關(guān)于模式的具體內(nèi)容見(jiàn)參考文獻(xiàn)［13］。

3.2 試驗(yàn)方案

采用不同的初始場(chǎng)方案進(jìn)行72 h預(yù)報(bào)設(shè)計(jì)了兩個(gè)試驗(yàn):

（1）試驗(yàn)I:采用模式上次預(yù)報(bào)的波譜條件初始化后,進(jìn)行預(yù)報(bào)72 h；

（2）試驗(yàn)Ⅱ:采用本同化方案得到的初始場(chǎng)進(jìn)行72 h預(yù)報(bào)。

海浪模式的計(jì)算區(qū)域?yàn)?7.5°～62°N,155°～142.5°W ,如圖2。模式的網(wǎng)格分辨率為0.25°× 0.25°,網(wǎng)格總數(shù)為51×19；波譜頻率范圍為0.042～0.413 Hz,離散為25個(gè),波譜方向分辨率為15°,離散為24個(gè)方向。模式計(jì)算的時(shí)間步長(zhǎng)為300 s,模式預(yù)報(bào)結(jié)果每間隔1 h輸出1次。模式風(fēng)場(chǎng)驅(qū)動(dòng)采用CCM P風(fēng)場(chǎng)數(shù)據(jù)（本文中預(yù)報(bào)所用的風(fēng)場(chǎng)都不是實(shí)際的預(yù)報(bào)風(fēng)場(chǎng),故在分析預(yù)報(bào)誤差時(shí)可忽略風(fēng)場(chǎng)預(yù)報(bào)誤差導(dǎo)致的海浪預(yù)報(bào)誤差）,CC M P風(fēng)場(chǎng)的空間分辨率是0.25°×0.25°,每次時(shí)間間隔為6 h。水深資料采用美國(guó)國(guó)家地球物理數(shù)據(jù)中心（N G D C）提供的E T O P O-1數(shù)據(jù)。模式物理過(guò)程包括深度誘導(dǎo)波破碎、底摩擦耗散作用、白浪耗散作用、非線性4項(xiàng)波相互作用等主要物理過(guò)程。這些物理過(guò)程的參數(shù)化表達(dá)式都采用模式的默認(rèn)設(shè)置。

試驗(yàn)過(guò)程分為同化得到初始場(chǎng)和進(jìn)行72 h預(yù)報(bào)兩個(gè)過(guò)程。以2008年10月21日08時(shí)起預(yù)報(bào)72 h為例說(shuō)明,圖3表示這個(gè)同化和預(yù)報(bào)的過(guò)程,預(yù)報(bào)時(shí)間是2008年10月21日08時(shí),同化過(guò)程是2008年10月20日08時(shí)開(kāi)始同化至10月21日08時(shí)獲得同化后的初始場(chǎng),開(kāi)始預(yù)報(bào),這24 h內(nèi)同化5個(gè)時(shí)次的波高觀測(cè)資料,每個(gè)同化時(shí)次,利用該時(shí)次之前15 d（即m取值為15）的每小時(shí)實(shí)測(cè)資料,對(duì)模式計(jì)算結(jié)果進(jìn)行有效波高的最優(yōu)插值,然后根據(jù)上述的方案,對(duì)格點(diǎn)上有效波高的最優(yōu)插值波高和模式計(jì)算的海浪譜進(jìn)行波譜分解和重構(gòu),之后,數(shù)值模式繼續(xù)同化,具體流程如圖3。

3.3 結(jié)果分析

（1）初始場(chǎng)對(duì)比

用浮標(biāo)的有效波高實(shí)測(cè)值對(duì)模式計(jì)算的波高場(chǎng)（下稱背景場(chǎng)）進(jìn)行最優(yōu)插值,上式中對(duì)L和α、β、γ等經(jīng)驗(yàn)常數(shù)的取值基于前人的研究成果［2］,α= 0,β= 2, γ= 0.5,L = 40 k m,圖2中給出了浮標(biāo)站點(diǎn)的位置和同化浮標(biāo)資料的計(jì)算格點(diǎn),從圖中可見(jiàn)共有51個(gè)計(jì)算格點(diǎn)的波高經(jīng)過(guò)最優(yōu)插值訂正,形成新的波高分析場(chǎng),如圖4所示。圖4a是2008年10月20日08時(shí)的模式背景場(chǎng),圖4b是該時(shí)刻的分析場(chǎng),背景場(chǎng)波高偏小的現(xiàn)象,在同化格點(diǎn)上得到了修正。圖4c和4d分別為2008年10月20日14時(shí)的模式波高背景場(chǎng)和波高分析場(chǎng),可見(jiàn)到類似的現(xiàn)象,最優(yōu)插值對(duì)波高的影響也由局部范圍擴(kuò)大到整個(gè)計(jì)算區(qū)域。

圖2 海浪模式計(jì)算區(qū)域Eig .2 The waves mode calculation area

圖3 同化和預(yù)報(bào)的流程Eig .3 The flowchart of assimilation and forecast

圖4 模式計(jì)算有效波高場(chǎng)和分析場(chǎng)（m）Eig .4 The significant wave height of model result and analysis fields（m）a .2008年10月20日08時(shí)模式背景場(chǎng),b .2008年10月20日08時(shí)分析場(chǎng),c .2008年10月20日14時(shí)模式背景場(chǎng),d .2008年10月20日14時(shí)分析場(chǎng)a .Significant wave height background field at 08:00 of October 20th,2008,b .significant wave height analyticalfield at 08:00 of October 20th, 2008,c .significant wave height background field at 14:00 of October 20th,2008,d .significant wave height analyticalfield at 14:00 of October 20th,2008

得到有效波高的分析場(chǎng)后,構(gòu)建卡爾曼濾波系統(tǒng)對(duì)分解過(guò)的波譜進(jìn)行修正和重構(gòu)得到模式預(yù)報(bào)所需的海浪譜初始場(chǎng)。對(duì)每個(gè)同化計(jì)算的格點(diǎn),用同化時(shí)刻之前的15 d的分析場(chǎng)和模式輸出的海浪譜,進(jìn)行波譜修正和重構(gòu),得到同化時(shí)刻的海浪譜。以2008 年10月21日08時(shí)計(jì)算格點(diǎn)（57.75°N,144.75°W）為例,圖5a是同化之前的海浪方向譜,圖5b是同化后的海浪方向譜,該點(diǎn)同化后的方向譜在能量的大小上有一定增大,計(jì)算格點(diǎn)（57.75°N,144.50°W）和計(jì)算格點(diǎn)（59.50°N,144.50°W）同化后的方向譜在能量的大小和主要頻率及譜形狀上發(fā)生了改變,分別如圖6和圖7所示。

圖5 57.75°N,144.75°W格點(diǎn)的方向譜的比較Eig .5 Directional spectra for the grid 57.75°N,144.75°W

圖6 57.75°N,144.50°W格點(diǎn)的方向譜的比較Eig .6 Directional spectra for the grid 57.75°N,144.50°W

（2）預(yù)報(bào)結(jié)果對(duì)比分析

試驗(yàn)中m取值為15,即用15 d的方向譜時(shí)間序列進(jìn)行譜分解和重構(gòu),取nk為9個(gè)主分量。如圖8所示,在預(yù)報(bào)時(shí)刻之前進(jìn)行24 h的資料同化,形成預(yù)報(bào)時(shí)刻的初始場(chǎng),然后進(jìn)行預(yù)報(bào),以24 h預(yù)報(bào)結(jié)果為例分析同化后的預(yù)報(bào)效果。圖8為連續(xù)1個(gè)月的浮標(biāo)站實(shí)測(cè)有效波高和兩個(gè)試驗(yàn)方案預(yù)報(bào)結(jié)果的時(shí)間序列比較圖,時(shí)間間隔為1 h。該圖和有效波高均方根誤差的比較表明未同化的模式預(yù)報(bào)的有效波高值較實(shí)測(cè)值普遍偏小,同化后的預(yù)報(bào)結(jié)果與實(shí)測(cè)值接近,特別是在實(shí)測(cè)波高較大時(shí)同化后效果改進(jìn)比較明顯, 表2給出了這7個(gè)測(cè)站1個(gè)月連續(xù)24 h預(yù)報(bào)效果的統(tǒng)計(jì)結(jié)果。

圖7 59.50°N,144.50°W格點(diǎn)的方向譜的比較Eig .7 Directional spectra for the grid 59.50°N,144.50°W

分別進(jìn)行了24 h、36 h、48 h、72 h同化預(yù)報(bào)試驗(yàn),考察不同預(yù)報(bào)時(shí)效的效果,結(jié)果如表2所示。由表2給出的7個(gè)測(cè)站1個(gè)月的預(yù)報(bào)誤差可見(jiàn),采用同化方法進(jìn)行預(yù)報(bào)的結(jié)果在24 h內(nèi)預(yù)報(bào)效果得到了明顯改進(jìn)。將7個(gè)測(cè)站的結(jié)果平均后,發(fā)現(xiàn)同化后有效波高均方根誤差較未同化的效果降低了0.13 m；48 h預(yù)報(bào)效果是同化后有效波高均方根誤差較未同化降低了0.05 m；72 h的預(yù)報(bào)結(jié)果中同化后和未同化的有效波高均方根誤差沒(méi)有差別,因此可以認(rèn)為同化后初始場(chǎng)對(duì)預(yù)報(bào)效果的改善作用已經(jīng)消失了。

表2 預(yù)報(bào)結(jié)果與浮標(biāo)觀測(cè)值的均方根誤差統(tǒng)計(jì)（單位:m）Tab .2 Statistical parameters R M SE of observed and calculated significant wave height（Unit:m）

通過(guò)兩個(gè)測(cè)站同化后和未同化的預(yù)報(bào)結(jié)果的比較來(lái)分析同化作用對(duì)預(yù)報(bào)影響的時(shí)效。圖9a為46077站從10月20日08時(shí)至10月21日08時(shí)每6 h同化一次后獲得10月21日08時(shí)的初始場(chǎng)后進(jìn)行預(yù)報(bào)的有效波高時(shí)間序列,圖9b為46105站的結(jié)果,在10月21日08時(shí)之后的30 h內(nèi)預(yù)報(bào)波高和實(shí)測(cè)符合較好,此后隨著時(shí)效增長(zhǎng),效果變差,從46077站的結(jié)果看,50 h后同化和未同化的結(jié)果完全一致。值得注意的是在10月21 日08時(shí)之前的24 h內(nèi)在每間隔6 h的時(shí)刻同化了觀測(cè)資料后,波高會(huì)出現(xiàn)“跳躍式變化”,模式計(jì)算值會(huì)出現(xiàn)突然大于或小于實(shí)測(cè)波高的現(xiàn)象,波高變化較快。這個(gè)現(xiàn)象可能是由于間斷的將實(shí)測(cè)資料與模式進(jìn)行同化,同化后的波浪譜是在接近實(shí)測(cè)波高的條件下重構(gòu)得到的,與模式的平衡關(guān)系被打破了,所以模式的計(jì)算結(jié)果也體現(xiàn)了模式對(duì)同化資料后的平衡調(diào)整過(guò)程。

圖8 2010年10月有效波高預(yù)報(bào)結(jié)果時(shí)間序列比較Eig .8 Time series of observed and calculated significant wave height for the assimilation experiments on October 2010

從整個(gè)區(qū)域的有效波高等值線也可發(fā)現(xiàn)這個(gè)現(xiàn)象。圖10a為同化后的初始場(chǎng),圖10b為未同化的初始場(chǎng),圖11a和圖11b分別為未同化和同化后6 h預(yù)報(bào)結(jié)果,波高的差值分布從同化的格點(diǎn)區(qū)域擴(kuò)散到了整個(gè)計(jì)算區(qū)域。圖12a、12b、12c、12d分別為6 h、12 h、36 h、48 h預(yù)報(bào)結(jié)果的差值,從整個(gè)區(qū)域有效波高預(yù)報(bào)結(jié)果的差值變化圖可見(jiàn),差值從同化的計(jì)算格點(diǎn)逐漸擴(kuò)散到近岸,并減小為0,這也可以認(rèn)為同化后的初始場(chǎng)對(duì)模式的作用是以波浪能量的形式傳播影響整個(gè)計(jì)算區(qū)域,因此同化影響的時(shí)效和范圍也可能與波浪傳播的速度有關(guān)。

圖9 2010年10月有效波高結(jié)果和實(shí)測(cè)比較Eig .9 The time series of observed and calculated significant wave height on October 2010

圖10 2010年10月20日08時(shí)有效波高（m）預(yù)報(bào)初始場(chǎng)Eig .10 The initialfields of significant wave height（m）on October 20th,2010

4 小結(jié)和討論

（1）提出一種基于海浪譜分解與重構(gòu)的海浪有效波高資料同化方案:該方采用最優(yōu)插值方法得到分析波高場(chǎng)；在模式的波浪能量密度譜和有效波高分析值之間引入一個(gè)變異系數(shù)矩陣,描述模式的波浪能量密度譜值和有效波高之間的誤差關(guān)系,以此為狀態(tài)向量構(gòu)建卡爾曼濾波系統(tǒng),對(duì)分解過(guò)的波譜進(jìn)行修正,進(jìn)而重構(gòu)海浪譜初始場(chǎng)。

（2）通過(guò)利用美國(guó)阿拉斯加灣的7個(gè)浮標(biāo)站進(jìn)行同化和預(yù)報(bào)的數(shù)值試驗(yàn),對(duì)連續(xù)1個(gè)月的24 h預(yù)報(bào)結(jié)果進(jìn)行統(tǒng)計(jì)表明:采用同化方案后預(yù)報(bào)結(jié)果的有效波高均方根誤差比未同化的結(jié)果降低了0.13 m；同化后初始場(chǎng)對(duì)預(yù)報(bào)效果的影響可持續(xù)36 h左右。

圖11 6 h有效波高（m）預(yù)報(bào)結(jié)果Eig .11 The significant wave height（m）fields of 6 h forecast

圖12 未同化和同化后有效波高預(yù)報(bào)結(jié)果之差（m）Eig .12 Difference of significant wave height for experiments with and without assimilation（m）

（3）文中對(duì)最優(yōu)插值的幾個(gè)參數(shù)的設(shè)置主要采用了前人研究成果；對(duì)同化方法中波浪方向譜時(shí)間序列的天數(shù)選取和譜分解的主分量個(gè)數(shù)設(shè)置沒(méi)有詳細(xì)討論,而是直接給出了經(jīng)驗(yàn)取值,這一點(diǎn)還需進(jìn)一步研究。

（4）本文的研究只直接同化了有效波高資料,對(duì)觀測(cè)波周期等波浪要素的直接同化也將開(kāi)展進(jìn)一步研究,并采用更多源的觀測(cè)資料來(lái)檢驗(yàn)同化預(yù)報(bào)的效果。

參考文獻(xiàn):

［1］ Tolman H L,Banner M L,Kaihatu J M . The N O PP operational wave modelimprovement project［J］. Ocean M odelling,2013,70（5）:2 - 10 .

［2］ Komen G J .Introduction to wave models and assimilation of satellite data in wave models（ocean waves）［M］. Cambridge:Cambridge U niversity Press,1996:554 .

［3］ Hasselmann K,Hasselmann S,Bauer E,et al. Development of a satellite SA R image spectra and altimeter wave height data assimilation system for E RS - 1［R］. Technical Report N ASA - C R - 182685,1988 .

［4］ Bauer E,Hasselmann S,Hasselmann K,et al. Validation and assimilation of Seasat altimeter wave heights using the W A M wave model［J］.Journal of Geophysical Research:Oceans（1978 - 2012）,1992,97（C8）:12671 - 12682 .

［5］ Hasselmann S,Lionello P,Hasselmann K . An optimalinterpolation scheme for the assimilation of spectral wave data［J］.Journal of Geophysical Research:Oceans（1978 - 2012）,1997,102（C7）:15823 - 15836 .

［6］ Greenslade D J M ,Young I R . The impact ofinhomogenous background errors on a global wave data assimilation system［J］.Journal of Atmospheric & Ocean Science,2005,10（2）:61 - 93 .

［7］ E m manouil G,Galanis G,Kallos G . A new methodology for using buoy measurementsin sea wave data assimilation［J］. Ocean Dynamics,2010,60 （5）:1205 - 1218 .

［8］ 張志旭,齊義泉,施平,等.衛(wèi)星高度計(jì)波高資料的同化試驗(yàn)分析［J］.海洋學(xué)報(bào),2003,25（5）:21 - 28 . Zhang Zhixu,Qi Yiquan,Shi Ping,et al. Preliminary study on assimilation of significant wave heights from T/P altimeter［J］. Haiyang Xuebao, 2003,25（5）:21 - 28 .

［9］ 王毅.S W A N模式及數(shù)據(jù)同化技術(shù)在海浪預(yù)報(bào)中的試驗(yàn)研究和應(yīng)用［D］.青島:中國(guó)海洋大學(xué),2011 . W ang Yi.Study and application of S W A N model and data assimilation technology in wave forecast［D］. Qingdao:China Ocean U niversity,2011 .

［10］ 王毅,余宙文.衛(wèi)星高度計(jì)波高數(shù)據(jù)同化對(duì)西北太平洋海浪數(shù)值預(yù)報(bào)的影響評(píng)估［J］.海洋學(xué)報(bào),2009,31（6）:1 - 8 . W ang Yi,Yu Zhouwen . Validation ofimpact of assimilation of altimeter satellite significant wave height on wave forecastin the northwest Pacific ［J］. Haiyang Xuebao,2009,31（6）:1 - 8 .

［11］ Eeng X,Zheng J,Yan Y . W ave spectra assimilation in typhoon wave modeling for the East China Sea［J］. Coastal Engineering,2012,69（11）:29 - 41 .

［12］ Guo Y,H ou Y,Zhang C,et al. A background error covariance model of significant wave height employing M onte Carlo simulation［J］. Chinese Journal of Oceanology and Limnology,2012,30（5）:814 - 821 .

［13］ Tolman H L . User manual and system docu mentation of W AV E W AT C HⅢT M version 3.14［S］. Technical Note,M M A B Contribution,2009: 276 .

中圖分類號(hào):P731.22

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):0253-4193（2016）03-0015-12

收稿日期:2014-12-05；

修訂日期:2015-04-02。

基金項(xiàng)目:江蘇省自然科學(xué)基金（B K20150711）；國(guó)家自然科學(xué)基金資助項(xiàng)目（11102232）。

作者簡(jiǎn)介:毛科峰（1981—）,男,湖南省常德市人,講師,研究方向?yàn)楹Ｑ蟓h(huán)境預(yù)報(bào)與海洋調(diào)查。E-mail:maomaopla @ 163 .com

毛科峰,蕭中樂(lè),宋海波,等.基于海浪譜分解與重構(gòu)的資料同化試驗(yàn)［J］.海洋學(xué)報(bào),2016,38（3）:15 - 26,doi:10.3969/j.issn . 0253-4193.2016.03.002

M ao Kefeng,Xiao Zhongle,Song Haibo,et al. Assimilation experiment based on wave spectrum decomposition and reconstruction［J］. Haiyang Xuebao,2016,38（3）:15 - 26,doi:10.3969/j.issn .0253-4193.2016.03.002

Assimilation experiment based on wave spectrum decomposition and reconstruction

M ao Kefeng1,Xiao Zhongle2,Song Haibo3,Zhong Yifei4
（1 .Instituteof Meteorology and Oceanography,PLA University of Science and Technology,Nanjing 211101,China；2 . The 96631 Unitof PLA ,Beijing 102208,China；3 .The 92721 Unitof PLA ,Zhoushan 316000,China；4 . The 95857 Unitof PLA ,Hengyang 102208,China）

Abstract:In the near-shore wave forecasting model,the ocean buoy data assimilation method based on wave spectral decomposition and optimization is proposed . The calculated wave energy spectru m before the initialtime is studied with orthogonal decomposition,and the result,combined with synchronous buoy observations of significant wave height value is used to construct a Kalman filtering system,and the initial wave energy density spectru m of wave modelis revisesd by the multi-time significant wave height value . This method has been applied to 72 h wave forecast experiments with assimilation 7 buoy's significant wave height data in the Gulf of Alaska . One month experiments show that method can improve the forecasts of significant wave height at different degree of different prediction time . The mean square errorfor 24 h forecast of significant wave heightreduces 0.13 m . M eanwhile,the effect ofinitialfield after assimilation to forecast will extend 36 h or so,but the assimilation effectis weakened by extend the prediction time .

Key words:ocean wave spectru m；significant wave height；assimilation；WAV E WAT C H -Ⅲmodel；Kalman filtering