四維變分同化技術(shù)在風(fēng)暴潮數(shù)值模擬中的應(yīng)用

王宗辰，于福江,2，原野,2（.國家海洋環(huán)境預(yù)報(bào)中心，北京0008；2.國家海洋局海洋災(zāi)害預(yù)報(bào)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，北京0008）

四維變分同化技術(shù)在風(fēng)暴潮數(shù)值模擬中的應(yīng)用

王宗辰1，于福江1,2，原野1,2
（1.國家海洋環(huán)境預(yù)報(bào)中心，北京100081；2.國家海洋局海洋災(zāi)害預(yù)報(bào)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，北京100081）

摘要：利用區(qū)域海洋模式ROMS(Regional Ocean Modelling System)及其四維變分同化模塊，建立了具有資料同化能力的東中國海風(fēng)暴潮數(shù)值模式，通過將海洋站水位觀測(cè)資料同化到風(fēng)暴潮模式中，提高了模式對(duì)風(fēng)暴潮的模擬精度。四維變分同化技術(shù)能夠在整個(gè)同化時(shí)間窗口保持動(dòng)力協(xié)調(diào)，使模擬結(jié)果在該時(shí)間窗口內(nèi)最大程度的靠近觀測(cè)，同時(shí)，得到了最優(yōu)預(yù)報(bào)初始場(chǎng)。利用該模式，對(duì)兩次溫帶風(fēng)暴潮過程進(jìn)行了數(shù)值模擬，結(jié)果表明：在同化窗口內(nèi)，同化對(duì)模擬精度有明顯的提高；結(jié)束同化之后，得到的最優(yōu)預(yù)報(bào)初始場(chǎng)對(duì)臨近預(yù)報(bào)精度也有一定提高。

關(guān)鍵詞：區(qū)域海洋模式ROMS；四維變分同化；風(fēng)暴潮數(shù)值模擬；最優(yōu)預(yù)報(bào)初始場(chǎng)

1　引言

目前風(fēng)暴潮數(shù)值預(yù)報(bào)結(jié)果多是由確定性數(shù)值模式計(jì)算得到，但是采用確定性風(fēng)暴潮數(shù)值模式，有時(shí)模擬結(jié)果和實(shí)際情況相比仍然有較大的誤差[1]。這些誤差的來源主要包括3個(gè)方面：一是動(dòng)力模式對(duì)真實(shí)物理過程的各種近似處理，例如模式物理參數(shù)化過程；二是大氣風(fēng)應(yīng)力和模式開邊界等輸入資料的不確定性；三是由數(shù)值計(jì)算方法引起的誤差。解決這一問題比較有效的方法之一是將風(fēng)暴潮數(shù)值模式和寶貴的觀測(cè)資料通過最優(yōu)的方法融合起來，發(fā)展具有同化能力的數(shù)值模式。近年來，各種同化方法在風(fēng)暴潮預(yù)報(bào)領(lǐng)域開展了較為廣泛的應(yīng)用[1-10]，其中比較有代表性的是Kalman濾波在西北歐國家風(fēng)暴潮業(yè)務(wù)化中的應(yīng)用[2-4]，和變分同化在參數(shù)估計(jì)、狀態(tài)估計(jì)中的應(yīng)用[6-9]。

Kalman濾波方法的缺點(diǎn)是對(duì)觀測(cè)資料序列有嚴(yán)格限制，即觀測(cè)資料不能出現(xiàn)缺測(cè)，并且要求資料序列必須是等時(shí)間間隔[1]。本文所采用的四維變分同化算法突破了這些限制，在個(gè)別站點(diǎn)缺測(cè)或者資料序列時(shí)間間隔不等的情況下，模式仍可以順利運(yùn)行得到結(jié)果。Sasaki[11]第一次提出了滿足控制方程強(qiáng)約束條件的伴隨方法，完成了四維變分同化的基本理論。20世紀(jì)末，伴隨方法被引入到氣象數(shù)據(jù)分析和同化領(lǐng)域[12]，之后，又被引入到全球海洋狀態(tài)分析和海洋環(huán)流模式中。Zhang等[6-7]通過將模式得到的“偽觀測(cè)數(shù)據(jù)”（包括余水位和近岸潮汐，余水位指總水位減去潮汐）同化到二維POM（Princeton Ocean Model）中，估計(jì)了風(fēng)拖曳系數(shù)和側(cè)邊界的潮汐等。Peng等[8-9]在臺(tái)風(fēng)風(fēng)暴潮研究中應(yīng)用四維變分同化技術(shù)，將“偽觀測(cè)數(shù)據(jù)”（由高分辨率模式產(chǎn)生）同化到模式中，對(duì)臺(tái)風(fēng)風(fēng)暴潮進(jìn)行了后報(bào)試驗(yàn)，結(jié)果表明只同化水位就可以有效的改進(jìn)預(yù)報(bào)結(jié)果。但遺憾的是，超過同化的3 h窗口，模式結(jié)果的改進(jìn)大概只能持續(xù)5 h左右，這也是由臺(tái)風(fēng)風(fēng)暴潮的特點(diǎn)決定的，其風(fēng)場(chǎng)變化劇烈，增水速度快，持續(xù)時(shí)間短，一般在6 h到12 h之間。本文則選取了溫帶風(fēng)暴潮作為試驗(yàn)對(duì)象，主要考慮其風(fēng)場(chǎng)變化相對(duì)緩慢，增減水的時(shí)間尺度更長，可以持續(xù)12 h到36 h，這樣，同化窗口可以設(shè)置得更寬，最優(yōu)預(yù)報(bào)初始場(chǎng)就可以包含更多的水位信息。

本文選擇ROMS作為水動(dòng)力模式，第一個(gè)原因是其具備模擬風(fēng)暴潮的能力。韓國氣象廳的You等[13]通過3次影響韓國的臺(tái)風(fēng)風(fēng)暴潮過程對(duì)ROMS在西北太平洋的應(yīng)用進(jìn)行了檢驗(yàn)，模擬表明ROMS的預(yù)報(bào)結(jié)果和韓國現(xiàn)有的風(fēng)暴潮預(yù)報(bào)系統(tǒng)（RTSM, Regional Tide/Storm Model）結(jié)果相當(dāng)，而且有更大的空間變化，以后很有可能用ROMS取代現(xiàn)有的業(yè)務(wù)化預(yù)報(bào)系統(tǒng)；美國馬里蘭大學(xué)的Li等[14]和愛爾蘭氣象中心的Wang等[15]也都利用ROMS分別對(duì)發(fā)生在切薩皮克灣和愛爾蘭附近的風(fēng)暴潮過程進(jìn)行了模擬，結(jié)果令人滿意。另一個(gè)重要原因是其四維變分同化技術(shù)的有效性已被驗(yàn)證，Powell等[16-17]和Broquet等[18-20]在北美邊緣海開展了一系列基于ROMS的四維變分同化研究，取得了不錯(cuò)的效果。

本文基于ROMS建立了具有資料同化能力的東中國海風(fēng)暴潮數(shù)值模式，并進(jìn)行了溫帶風(fēng)暴潮個(gè)例試驗(yàn)。文章接下來介紹模式的基本配置，然后簡要說明強(qiáng)約束增量型四維變分同化的算法及其實(shí)現(xiàn)，最后分析試驗(yàn)結(jié)果并加以總結(jié)和討論。

2　模式配置

ROMS是近10年來新發(fā)展起來并被廣泛認(rèn)可的一個(gè)三維、自由海面和基于地形跟隨坐標(biāo)的非線性斜壓原始方程模式，具體可參考文獻(xiàn)[21-23]。

本文選取的模擬區(qū)域?yàn)?14.233°—133°E，22.2°—41°N，包含渤海、黃海、東海以及西日本海和西北太平洋海域，水平網(wǎng)格分辨率為（1/10）°× （1/10）°，格點(diǎn)數(shù)為188×189，垂直方向分為20層，這樣的分辨率設(shè)置既可以較好的刻畫風(fēng)暴潮過程又能夠兼顧計(jì)算效率。

水深資料來自英國海洋資料中心（BODC, British Oceanographic Data Centre）的格點(diǎn)數(shù)據(jù)gebco_08（0.5′×0.5′），同時(shí)，近岸水深（h<100 m）由國家海洋環(huán)境預(yù)報(bào)中心業(yè)務(wù)化水深數(shù)據(jù)（2′×2′）插值代替，最小水深5 m，最大水深約為6000 m（見圖1）。岸線數(shù)據(jù)采用美國地球物理資料中心（NGDC, National Geophysical Data Center）最高分辨率的GSHHG-2.2.2數(shù)據(jù)。南、東、北三個(gè)方向設(shè)為開邊界，二維平均流速采用Flather邊界條件，即正壓模態(tài)的偏差以重力外波的速度（傳播出去；與之配合，水位采用Chapman邊界條件，即水位以淺水波速（播出去；三維流速采用輻射邊界條件。

ROMS的大氣強(qiáng)迫方式有2種方法可供選擇，因風(fēng)暴潮過程在近岸主要受風(fēng)強(qiáng)迫影響，所以只采用風(fēng)應(yīng)力作為外強(qiáng)迫，風(fēng)應(yīng)力由歐洲中長期天氣預(yù)報(bào)中心的再分析風(fēng)場(chǎng)（ERA-Interim）根據(jù)Large & Pond[24]的公式計(jì)算得到，空間分辨率0.75°×0.75°，時(shí)間分辨率3 h，內(nèi)模積分時(shí)間步長取180 s，外模取9 s。

3　四維變分同化算法和實(shí)現(xiàn)

3.1四維變分同化算法

資料同化的過程是將確定性模式結(jié)果和觀測(cè)數(shù)據(jù)融合在一起，生成一個(gè)對(duì)海洋狀態(tài)的重新估計(jì)，這種估計(jì)應(yīng)該比確定性模式結(jié)果更加接近真值。四維變分同化技術(shù)通過動(dòng)力約束使模式結(jié)果與觀測(cè)之間距離達(dá)到最小化，動(dòng)力協(xié)調(diào)一致是這種技術(shù)的優(yōu)點(diǎn)；缺點(diǎn)是通過迭代的方法求解伴隨方程解的過程非常消耗計(jì)算資源，雖然可以通過具體的迭代次數(shù)來控制計(jì)算時(shí)間，但實(shí)際情況是，為了保證目標(biāo)函數(shù)的收斂，狀態(tài)變量的維數(shù)不宜太大，這也是水平網(wǎng)格選擇（1/10）°×（1/10）°的主要原因。

ROMS自帶的四維變分同化技術(shù)由Courtier[25]首先提出，現(xiàn)已在業(yè)務(wù)化數(shù)值天氣預(yù)報(bào)模式中展開了廣泛的應(yīng)用，下面簡單介紹其原理和方法。

根據(jù)貝葉斯理論，得到有觀測(cè)條件下“最優(yōu)初始場(chǎng)”（或稱“分析場(chǎng)”）的條件概率（表達(dá)式略）；再根據(jù)最大似然估計(jì)推出，若要使概率最大，則下面的表達(dá)式（目標(biāo)函數(shù)）取極小值，即

式中，δxk=xk-xk-1、db k-1=xb-xk-1、do k-1= yo(i)-G(xk-1)；x∈Rn表示控制變量；n代表控制變量的空間維數(shù)；xb∈Rn表示背景場(chǎng)估計(jì)值；δxk表示外循環(huán)的迭代增量；k代表外循環(huán)次數(shù)；yo(i)∈Rm表示觀測(cè)值；i代表觀測(cè)值的時(shí)間維數(shù)；m代表觀測(cè)值的空間維數(shù)；Gk-1是一個(gè)Rn→Rm的線性空間映射算子，用來實(shí)現(xiàn)切線性模式的向前積分和模式結(jié)果向觀測(cè)空間線性插值；B(x)∈Rn×n和R∈Rm×m分別代表背景誤差協(xié)方差和觀測(cè)誤差協(xié)方差的估計(jì)值，都是對(duì)稱正定矩陣。

另外循環(huán)次數(shù)為1，則k=1，一般的x0=xb，目標(biāo)函數(shù)簡化為：

對(duì)簡化了的目標(biāo)函數(shù)求導(dǎo)得：

使其等于0便可得到δx的解析解，但很遺憾，目前只能通過迭代算法使目標(biāo)函數(shù)向某一點(diǎn)收斂，ROMS所采用的是基于Lanczos方法的共軛梯度下降算法[26-28]。

3.1.1背景誤差協(xié)方差

到目前為止，如何準(zhǔn)確的構(gòu)造背景誤差協(xié)方差仍是一個(gè)非常具有挑戰(zhàn)性的問題，原因是無法獲得真值，同時(shí)又缺乏足夠的樣本資料。本文參照“NMC”方法[29]，選取12個(gè)不同時(shí)效的預(yù)報(bào)值對(duì)背景誤差進(jìn)行估計(jì)。具體做法是，從12個(gè)不同的時(shí)刻冷啟動(dòng)模式，然后獲得同一時(shí)刻后報(bào)結(jié)果作為樣本，把樣本的平均值看做“真值”，12個(gè)樣本離差的平均值近似的看做背景誤差，背景誤差協(xié)方差可以表達(dá)為

3.1.2觀測(cè)誤差協(xié)方差

考慮到風(fēng)暴潮過程的時(shí)空尺度，海洋站的水位數(shù)據(jù)是一種比較理想的同化資料。因此，選取了中國沿海24個(gè)海洋站的逐時(shí)水位數(shù)據(jù)作為同化資料。為了簡化問題，本文不考慮潮汐以及潮汐和風(fēng)增水的非線性相互作用，因此選擇T_TIDE[30]作為調(diào)和分析的工具，對(duì)24個(gè)海洋站的水位進(jìn)行調(diào)和分析。利用調(diào)和分析的結(jié)果推算潮汐，總水位減去潮汐，就得到了風(fēng)暴潮，風(fēng)暴潮誤差統(tǒng)一設(shè)定為0.02 m，假定R是一個(gè)對(duì)角矩陣，則R對(duì)角線上所有元素均為0.0004 m2。

3.2算法實(shí)現(xiàn)

設(shè)定同化窗口[t0,t1]，外循環(huán)設(shè)為1層，內(nèi)循環(huán)10層。根據(jù)本文的試驗(yàn)結(jié)果，這樣的設(shè)置可以保證目標(biāo)函數(shù)基本收斂。

（1）從t0開始熱啟動(dòng)非線性模式，向前積分到t1，得到[t0,t1]時(shí)刻內(nèi)的背景場(chǎng)xb，其中?表示非線性模式的轉(zhuǎn)移矩陣，f表示大氣強(qiáng)迫；

（2）從t0時(shí)刻開始，選擇一個(gè)δx（δx?xb），向前積分切線性模式到t1時(shí)刻，得到t1時(shí)刻的Gδx-d和J，其中?(t)=??/?x|xb；

（3）再從t1時(shí)刻開始，用R-1(Gδx-d)強(qiáng)迫伴隨模式向后積分到t0時(shí)刻，根據(jù)公式（3）可以得到?J/?δx；

（4）選擇共軛梯度下降算法不斷更新δx；

（5）重復(fù)步驟2—4，直到滿足收斂條件，即得到最優(yōu)初始場(chǎng)xb+δx。

4　數(shù)值模擬和結(jié)果

本文選擇了2011年初發(fā)生在渤海的兩次溫帶風(fēng)暴潮過程作為試驗(yàn)對(duì)象，分別是1月15日0時(shí)—1 月16日24時(shí)（世界時(shí)，后同）的一次減水過程和2月26日0時(shí)—2月27日24時(shí)一次增水過程，兩次過程的持續(xù)時(shí)間大約都為48 h，為了觀察后續(xù)變化，模式運(yùn)行的結(jié)束時(shí)間向后延長12 h。這兩次風(fēng)暴潮的影響范圍僅限于渤海區(qū)域，因此，雖然同化了24個(gè)海洋站水位資料，但每次過程只選出渤海區(qū)域風(fēng)暴潮位最大的4個(gè)海洋站作為分析對(duì)象。

4.1減水個(gè)例

4.1.1試驗(yàn)設(shè)計(jì)

控制試驗(yàn)（CR）：風(fēng)暴潮過程開始前48 h，冷啟動(dòng)模式運(yùn)行至風(fēng)暴潮過程結(jié)束；

同化試驗(yàn)1（DA1）：控制試驗(yàn)運(yùn)行48 h后，加入同化模塊，同化時(shí)間窗口設(shè)為向后12 h，然后停止同化，運(yùn)行至風(fēng)暴潮過程結(jié)束；

同化試驗(yàn)2（DA2）：控制試驗(yàn)運(yùn)行48 h后，加入同化模塊，同化時(shí)間窗口設(shè)為向后24 h，然后停止同化，運(yùn)行至風(fēng)暴潮過程結(jié)束；

同化試驗(yàn)3（DA3）：控制試驗(yàn)運(yùn)行48 h后，加入同化模塊，同化時(shí)間窗口設(shè)為向后36 h，然后停止同化，運(yùn)行至風(fēng)暴潮過程結(jié)束。

4.1.2結(jié)果分析

控制試驗(yàn)和3個(gè)同化試驗(yàn)?zāi)M結(jié)果與觀測(cè)水位的時(shí)間序列如圖2所示。可以看出，控制試驗(yàn)和同化試驗(yàn)都較好的模擬出了風(fēng)暴減水的趨勢(shì)，同化試驗(yàn)的模擬結(jié)果明顯優(yōu)于控制試驗(yàn)。在0—48 h期間，同化試驗(yàn)3的模擬結(jié)果和觀測(cè)最為接近，這表明36 h的同化窗口最為準(zhǔn)確的模擬出了本次風(fēng)暴減水過程。

圖2　1月15日0時(shí)—1月17日12時(shí)，4個(gè)海洋站的風(fēng)暴減水時(shí)間序列（0時(shí)是同化窗口的起始時(shí)刻，有色虛線表示對(duì)應(yīng)顏色同化試驗(yàn)的同化結(jié)束時(shí)刻）

同化試驗(yàn)3可以看做一個(gè)再分析過程，其模擬結(jié)果可以近似認(rèn)為是真值。計(jì)算了同化試驗(yàn)3和控制試驗(yàn)在減水過程中不同時(shí)刻模擬結(jié)果的絕對(duì)差值，其空間分布見圖3。可以看出，不同時(shí)刻的大值區(qū)分布不盡相同，但差值比較大的區(qū)域主要分布在靠近灣頂或近岸處，這表明確定性模式的預(yù)報(bào)水位絕對(duì)誤差在靠近灣頂或海灣處比較大。

對(duì)于風(fēng)暴潮數(shù)值預(yù)報(bào)，同化的目的是要得到一個(gè)比確定性模式更為接近實(shí)際的預(yù)報(bào)初始場(chǎng)，從而提高確定性模式的預(yù)報(bào)精度[1]。同化試驗(yàn)1同化了12 h的模擬結(jié)果可以看做最優(yōu)預(yù)報(bào)初始場(chǎng)，為了定量描述其對(duì)預(yù)報(bào)結(jié)果的改進(jìn)，定義水位改進(jìn)度為

式中，error表示誤差；abs表示絕對(duì)值算符。

利用同化試驗(yàn)1得到的最優(yōu)預(yù)報(bào)初始場(chǎng)對(duì)水位進(jìn)行預(yù)報(bào)，并計(jì)算了逐時(shí)的水位改進(jìn)度（見圖4），同時(shí)，設(shè)定當(dāng)改進(jìn)度不小于30%時(shí)，最優(yōu)預(yù)報(bào)初始場(chǎng)對(duì)預(yù)報(bào)精度有顯著地提高。可以看到，對(duì)于這次風(fēng)暴減水過程，從初始時(shí)刻同化12 h之后得到的最優(yōu)預(yù)報(bào)初始場(chǎng)，不同程度的提高了4個(gè)海洋站的預(yù)報(bào)精度，不同站點(diǎn)的改進(jìn)時(shí)長在15—17 h之間，但對(duì)于芝罘島預(yù)報(bào)水位改進(jìn)度偏低的現(xiàn)象，原因不詳。同時(shí)可以看到，4站的改進(jìn)度曲線存在不同程度的波動(dòng)，龍口和芝罘島表現(xiàn)得尤其明顯，這和圖2中水位曲線的波動(dòng)是對(duì)應(yīng)出現(xiàn)的，可能是因?yàn)檎军c(diǎn)水位里包含了和模式不兼容的物理過程，引起了模式的不穩(wěn)定。

圖3　同化試驗(yàn)3與控制試驗(yàn)水位差絕對(duì)值的空間分布(單位：cm)

4.2增水個(gè)例

4.2.1試驗(yàn)設(shè)計(jì)

控制試驗(yàn)（CR）：風(fēng)暴潮過程開始前48 h，冷啟動(dòng)模式運(yùn)行至風(fēng)暴潮過程結(jié)束；

同化試驗(yàn)1（DA1）：控制試驗(yàn)運(yùn)行60 h后，加入同化模塊，同化時(shí)間窗口設(shè)為向后6 h，然后停止同化，運(yùn)行至風(fēng)暴潮過程結(jié)束；

同化試驗(yàn)2（DA2）：控制試驗(yàn)運(yùn)行60 h后，加入同化模塊，同化時(shí)間窗口設(shè)為向后12 h，然后停止同化，運(yùn)行至風(fēng)暴潮過程結(jié)束；

同化試驗(yàn)3（DA3）：控制試驗(yàn)運(yùn)行60 h后，加入同化模塊，同化時(shí)間窗口設(shè)為向后24 h，然后停止同化，運(yùn)行至風(fēng)暴潮過程結(jié)束。

4.2.2結(jié)果分析

控制試驗(yàn)和3個(gè)同化試驗(yàn)?zāi)M結(jié)果與觀測(cè)水位的時(shí)間序列如圖5所示。可以看出，控制試驗(yàn)和同化試驗(yàn)都較好的模擬出了增水趨勢(shì)，控制試驗(yàn)極大值點(diǎn)附近出現(xiàn)了較大誤差，其中黃驊站極值誤差達(dá)到了70 cm，極值的相位也存在1—4 h的滯后。3個(gè)同化試驗(yàn)對(duì)極值水位的預(yù)報(bào)精度都有提高，同化試驗(yàn)3的模擬結(jié)果和觀測(cè)最為接近，黃驊站的最大增水誤差減少了約30 cm，對(duì)極值點(diǎn)位相的訂正也非常明顯。

圖4　同化試驗(yàn)1中4個(gè)海洋站預(yù)報(bào)結(jié)果的逐時(shí)水位改進(jìn)度（虛線代表30%的改進(jìn)度）

圖5　2月26日0時(shí)—2月28日12時(shí)，4個(gè)海洋站風(fēng)暴增水時(shí)間序列（黑色虛線表示同化窗口的起始時(shí)刻，其他有色虛線表示對(duì)應(yīng)顏色同化試驗(yàn)的同化結(jié)束時(shí)刻）

同化試驗(yàn)3同樣可以看做一個(gè)再分析過程，其模擬結(jié)果可以近似認(rèn)為是真值。計(jì)算了同化試驗(yàn)3和控制試驗(yàn)在4個(gè)臨近增水極大值時(shí)刻模擬結(jié)果的絕對(duì)差值，其空間分布見圖6。分析結(jié)果與減水試驗(yàn)相似，不再贅述。

圖6　同化試驗(yàn)3與控制試驗(yàn)水位差絕對(duì)值的空間分布(單位：cm)

同化試驗(yàn)1在6時(shí)的模擬結(jié)果作為最優(yōu)預(yù)報(bào)初始場(chǎng)，對(duì)水位進(jìn)行了預(yù)報(bào)。對(duì)于風(fēng)暴增水過程，我們更加關(guān)注增水的極值和極值出現(xiàn)的時(shí)間，因此根據(jù)觀測(cè)選擇了每個(gè)站點(diǎn)增水極大值出現(xiàn)時(shí)刻前后各3 h共7 h的預(yù)報(bào)結(jié)果進(jìn)行分析，表1列出了同化試驗(yàn)1中這7 h的逐時(shí)水位改進(jìn)度。從表中可以看到，大部分時(shí)刻的預(yù)報(bào)精度都有所提高，但同時(shí)，曹妃甸和龍口兩站的改進(jìn)度出現(xiàn)了一些負(fù)值和小值，這和圖5中兩站水位曲線波動(dòng)現(xiàn)象是對(duì)應(yīng)的，可能是因?yàn)檎军c(diǎn)水位里包含了和模式不兼容的物理過程，引起了模式的不穩(wěn)定。把同化試驗(yàn)2在12時(shí)的模擬結(jié)果作為最優(yōu)預(yù)報(bào)初始場(chǎng)，對(duì)水位進(jìn)行了預(yù)報(bào)，結(jié)果顯示預(yù)報(bào)改進(jìn)度有進(jìn)一步提高，顯然是因?yàn)轭A(yù)報(bào)初始場(chǎng)更加靠近增水極值時(shí)刻。

表1　同化試驗(yàn)1在4個(gè)海洋站增水極值時(shí)刻前后3 h的逐時(shí)水位改進(jìn)度（單位：%）

5　總結(jié)和討論

本文利用ROMS及其強(qiáng)約束增量型四維變分同化模塊建立了中國渤、黃、東海風(fēng)暴潮數(shù)值預(yù)報(bào)模式，通過將觀測(cè)余水位同化到模式中，改進(jìn)了同化窗口內(nèi)的風(fēng)暴潮數(shù)值模擬結(jié)果。同時(shí)，資料同化技術(shù)可以提供更為合理的預(yù)報(bào)初始場(chǎng)，對(duì)風(fēng)暴潮的臨近預(yù)報(bào)結(jié)果有一定改進(jìn)。對(duì)于風(fēng)暴減水過程，30%的預(yù)報(bào)水位改進(jìn)可以持續(xù)約15—17 h；對(duì)于風(fēng)暴增水過程，預(yù)報(bào)結(jié)果的改進(jìn)程度和最優(yōu)預(yù)報(bào)初始場(chǎng)時(shí)刻有關(guān)，其越接近增水極值時(shí)刻，對(duì)風(fēng)暴增水預(yù)報(bào)結(jié)果的改進(jìn)越大。

目前，國家海洋環(huán)境預(yù)報(bào)中心的溫帶風(fēng)暴潮業(yè)務(wù)化預(yù)報(bào)每天運(yùn)行兩次，分別是7時(shí)和15時(shí)，發(fā)布預(yù)報(bào)時(shí)間為8時(shí)和16時(shí)。以7時(shí)開始運(yùn)行為例，模式所采用的同化資料只能是7時(shí)之前的海洋站數(shù)據(jù)，風(fēng)速為預(yù)報(bào)中心業(yè)務(wù)化風(fēng)場(chǎng)預(yù)報(bào)結(jié)果，這就要求模式要在1 h之內(nèi)運(yùn)行結(jié)束。根據(jù)本文的試驗(yàn)結(jié)果，24 h的預(yù)報(bào)結(jié)果5 min便可計(jì)算得到，同化及預(yù)報(bào)共計(jì)需要約50 min，可以滿足8時(shí)發(fā)報(bào)的要求。實(shí)際上，預(yù)報(bào)中心可以實(shí)時(shí)接收海洋站數(shù)據(jù)，只要風(fēng)場(chǎng)在預(yù)報(bào)時(shí)效內(nèi)，風(fēng)暴潮模式隨時(shí)可以運(yùn)行。

風(fēng)暴潮過程基本可以看做一個(gè)外強(qiáng)迫問題，即初始場(chǎng)只在前期對(duì)預(yù)報(bào)結(jié)果影響比較大，隨著時(shí)間推移，風(fēng)強(qiáng)迫的作用變得更加重要。本文的一個(gè)隱含假設(shè)是風(fēng)應(yīng)力等其他資料是準(zhǔn)確的，模式也是準(zhǔn)確的，僅海洋的初始狀態(tài)存在誤差，但實(shí)際情況是，即便通過同化技術(shù)訂正了初始場(chǎng)，預(yù)報(bào)結(jié)果和觀測(cè)相比仍存在系統(tǒng)性誤差。因此，在模式?jīng)]有改動(dòng)的情況下，可以認(rèn)為風(fēng)應(yīng)力、近岸地形和水深是誤差主要來源，通過資料同化技術(shù)更加合理的估計(jì)風(fēng)應(yīng)力是下一步的研究方向。

致謝：本文所用海洋站資料均來自國家海洋環(huán)境預(yù)報(bào)中心信息室，在此表示感謝！

參考文獻(xiàn)：

[1]于福江,張占海,林一驊.一個(gè)穩(wěn)態(tài)Kalman濾波風(fēng)暴潮數(shù)值預(yù)報(bào)模式[J].海洋學(xué)報(bào), 2002, 24(5): 26-35.

[2] Heemink A W. Storm surge prediction using Kalman filtering[D]. Twente: Twente University of Technology, 1986.

[3] Heemink A W, Mouthaan E E A, Roest M R T, et al. Inverse 3D shallow water flow modeling of the continental shelf [J]. Continental Shelf Research, 2002, 22(3): 465-484.

[4] Sorensen J V T, Madsen H. Efficient Kalman filter techniques for the assimilation of tide gauge data in three-dimensional modeling of the North Sea and Baltic Sea system[J]. Journal of Geophysical Research, 2004, 109(C3), doi: 10.1029/2003JC002144.

[5] Butler T, Altaf M U, Dawson C, et al. Data Assimilation within the Advanced Circulation (ADCIRC) Modeling Framework for Hurricane Storm Surge Forecasting[J]. Monthly Weather Review, 2012, 140(7): 2215-2231.

[6] Zhang A J, Parker B B, Wei E. Assimilation of water level data into a coastal hydrodynamic model by an adjoint optimal technique[J]. Continental Shelf Research, 2002, 22(14): 1909-1934.

[7] Zhang A, Wei E, Parker B B. Optimal estimation of tidal open boundary conditions using predicted tides and adjoint data assimilation technique[J]. Continental Shelf Research, 2003, 23 (11-13): 1055-1070.

[8] Peng S Q, Xie L. Effect of determining initial conditions by four-dimensional variational data assimilation on storm surge forecasting[J]. Ocean Modelling, 2006, 14(1-2): 1-18.

[9] Peng S Q, Xie L, Pietrafesa L J. Correcting the errors in the initial conditions and wind stress in storm surge simulation using an adjoint optimal technique[J]. Ocean Modelling, 2007, 18(3-4): 175-193.

[10] Lionello P, Sanna A, Elvini E, et al. A data assimilation procedure for operational prediction of storm surge in the northern Adriatic Sea[J]. Continental Shelf Research, 2006, 26(4): 539-553.

[11] Sasaki Y. Some basic formalisms in numerical variational analysis [J]. Monthly Weather Review, 1970, 98(12): 875-883.

[12] Le Dimet F X, Talagrand O. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[J]. Tellus A, 1986, 38(2): 97-110.

[13] You S H, Lee W-J, Moon K S. Comparison of storm surge tide predictions between a 2Doperational forecasting system operational forecast system, RTSMand ROMS[J]. Ocean Dynamics, 2010, 60(2): 443-459.

[14] Li M, Zhong L J, Boicourt W C, et al. Hurricane-induced storm surges, currents and destratification in a semi-enclosed bay[J]. Geophysical Research Letters, 2006, 33(2): L02604. doi: 10.1029/ 2005GL024992

[15] Wang S Y, McGrath R, Hanafin J, et al. The impact of climate change on storm surges over Irish waters[J]. Ocean Modelling, 2008, 25(1-2): 83-94.

[16] Powell B S, Arango H G, Moore A M, et al. 4DVAR data assimilation in the Intra-Americas Sea with the Regional Ocean Modeling System (ROMS)[J]. Ocean Modelling, 2008, 25(3-4): 130-145.

[17] Powell B S, Moore A M. Estimating the 4DVAR analysis error of GODAE products[J]. Ocean Dynamics, 2009, 59(1): 121-138.

[18] Broquet G, Edwards C A, Moore A M, et al. Application of 4D-variational data assimilation to the California Current System [J]. Dynamics of Atmospheres and Oceans, 2009, 48(1-3): 69-92.

[19] Broquet G, Moore A M, Arango H G, et al. Ocean state and surface forcing correction using the ROMS-IS4DVAR data assimilation system[J]. Mercator Ocean Quarterly Newsletter, 2009, 34: 5-13.

[20] Broquet G, Moore A M, Arango H G, et al. Corrections to ocean surface forcing in the California Current Systemusing 4D-variational data assimilation[J]. Ocean Modelling, 2011, 36 (1-2): 116-132.

[21] Shchepetkin A F, McWilliams J C. A method for computinghorizontal pressure-gradient force in an oceanic model with a nonaligned vertical grid[J]. Journal of Geophysical Research, 2003, 108(C3), doi: 10.1029/2001JC001047

[22] Shchepetkin A F, McWilliams J C. The regional oceanic modeling system (ROMS): a split explicit, free-surface, topographyfollowing-coordinate oceanic model [J]. Ocean Modelling, 2005, 9 (4): 347-404.

[23] Haidvogel D B, Arango H, Budgell W P, et al. Ocean forecasting in terrain-following coordinates: formulation and skill assessment of the Regional Ocean Modeling System[J]. Journal of Computational Physics, 2008, 227(7): 3595-3624.

[24] Large W G, Pond S. Open ocean momentum flux measurements inmoderatetostrongwinds[J]. Journal of Physical Oceanography, 1981, 11(3): 324-336.

[25] Courtier P, Thépaut J-N, Hollingsworth A. A strategy for operational implementation of 4D-Var, using an incremental approach[J]. Quarterly Journal of the Royal Meteorological Society, 1994, 120(519): 1367-1387.

[26] Fisher M. Minimization algorithms for variational data assimilation[R]//Proceedings of ECMWF Seminar Recent Developments in Numerical Methods for Atmospheric Modelling. U.K.: ECMWF Publication, 1998, 364-385.

[27] Tshimanga J. On a class of limited memory preconditioners for large-scalenonlinearleast-squaresproblems[D]. Namur, Belgium: Facultes Universitaires Notre-Dame de la Paix, 2007.

[28] Tshimanga J, Gratton S, Weaver A T, et al. Limited-memory preconditioners with application to incremental variational data assimilation[J]. Quarterly Journal of the Royal Meteorological Society, 2008, 134(632): 751-769.

[29] Parrish D F, Derber J C. The national meteorological center's spectral statistical-interpolation analysis system[J]. Monthly Weather Review, 1992, 120(8): 1747-1763.

[30] Pawlowicz R, Beardsley B, Lentz S. Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE[J]. Computers and Geosciences, 2002, 28(8): 929-937.

Application of four-dimensional variational data assimilation technique in storm surge simulations

WANG Zong-chen1，YU Fu-jiang1,2，YUAN Ye1,2
(1. National Marine Environmental Forecasting Center, Beijing 100081 China; 2. Key Laboratory of Research on Marine Hazards Forecasting, State Oceanic Administration, Beijing 100081 China)

Abstract：Base on Regional Ocean Modeling System(ROMS) with four-dimensional variational (4D-Var) data assimilation modules, a model is established for storm surge simulation in the East China Sea, and the deterministic model output is corrected by assimilating the available tidal gauge station data. The 4D-Var technique is able to compensate for the errors between modeled outputs and observations by containing dynamically consistent and persistent during a period of time(also called assimilation window), and the optimal initial condition for forecasting is obtained. Two applications on extra tropical storm surges that occurred in Bohai are performed by the model. The results show that the procedure is capable of strongly improving simulation accuracy, and the optimal initial condition effectively improves the forecasting accuracy in a short period.

Key words：Regional Ocean Modeling System (ROMS); 4D-Var data assimilation; storm surge simulation; optimal initial conditions

作者簡介：王宗辰（1988-），男，碩士研究生，主要從事風(fēng)暴潮數(shù)值模擬及同化技術(shù)研究。E-mail：781004920@163.com

基金項(xiàng)目：海洋公益性行業(yè)專項(xiàng)（200905013）

收稿日期：2014-01-06

DOI:10.11737/j.issn.1003-0239.2015.01.001

中圖分類號(hào)：P731.23

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1003-0239（2015）01-0001-09