丁駿，呂忻，姚雅倩，孟鑫，姜雪敏，吳旭云，葛建忠

(1.上海市海洋監(jiān)測(cè)預(yù)報(bào)中心，上海 200062；2.國(guó)家海洋局東海預(yù)報(bào)中心，上海 200136；3.華東師范大學(xué) 河口海岸學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200062)

1 引言

風(fēng)暴潮是世界上最嚴(yán)重的海洋災(zāi)害之一[1-2]，也是上海城市安全的主要威脅之一，風(fēng)暴增減水[3]是一種非周期性的對(duì)眾多因素敏感的又備受關(guān)注的復(fù)雜現(xiàn)象。目前，國(guó)內(nèi)外廣泛采用海洋站和衛(wèi)星等監(jiān)測(cè)到的海面高度與正常天文潮的代數(shù)差的方法來計(jì)算風(fēng)暴增減水，即風(fēng)暴潮大小[4-6]。由于風(fēng)暴潮與天文潮的非線性耦合效應(yīng)存在[7]，實(shí)際增減水并非兩者簡(jiǎn)單的線性相減。受海洋站空間分辨率限制，又常常以點(diǎn)代線、以點(diǎn)代面評(píng)估近岸海域風(fēng)暴潮大小，其誤差不斷累積。雖然衛(wèi)星高度計(jì)運(yùn)用到風(fēng)暴潮觀測(cè)的研究[8]已經(jīng)取得技術(shù)突破，但其時(shí)空分辨率和精度仍需提升。隨著理論技術(shù)發(fā)展與實(shí)測(cè)數(shù)據(jù)的豐富，風(fēng)暴潮數(shù)值預(yù)報(bào)[7]的準(zhǔn)確率和時(shí)效性進(jìn)一步提高，可以直觀反映風(fēng)暴增減水、漫灘過程，然而由于數(shù)值模式反映的物理過程是有限的，初始場(chǎng)和強(qiáng)迫場(chǎng)是近似的，模式參數(shù)是經(jīng)驗(yàn)的，預(yù)測(cè)結(jié)果不可避免存在誤差。在實(shí)踐中，用實(shí)測(cè)潮位分離法計(jì)算的增水值通常偏大，該值可達(dá)15～30 cm，而數(shù)值模擬的增水值往往偏小[9]。

集合卡爾曼濾波（Ensemble Kalman Filter，EnKF）[10-11]可以在受噪聲影響的測(cè)量值與數(shù)值模擬的預(yù)測(cè)值之間尋找最優(yōu)狀態(tài)估計(jì)。它基于集合預(yù)報(bào)思想，源于卡爾曼濾波理論和蒙特卡洛估計(jì)方法[12]，克服了卡爾曼濾波技術(shù)計(jì)算量大、只適用線性模式的缺點(diǎn)。由Evensen[10]在1994 年提出并應(yīng)用到海洋同化領(lǐng)域后，日漸成為國(guó)際上資料同化研究的熱點(diǎn)。在風(fēng)暴潮資料同化應(yīng)用方面，Butler 等[13]和Altaf 等[14]提出了一種基于EnKF 和ADCIRC（Advanced Circulation）風(fēng)暴潮模型耦合的風(fēng)暴潮預(yù)報(bào)方法，Colle 等[15]運(yùn)用En-KF 對(duì)11 個(gè)風(fēng)場(chǎng)資料同化驅(qū)動(dòng)ADCIRC 模型來模擬“桑迪”颶風(fēng)風(fēng)暴潮，Etala 等[16]展現(xiàn)了局地集合變換卡爾曼濾波（Local Ensemble Transform Kalman Filter，LETKF）在衛(wèi)星遙感與驗(yàn)潮站數(shù)據(jù)同化方面的優(yōu)勢(shì)，這些研究主要集中在EnKF 與風(fēng)暴潮模型、氣象模型的耦合同化，EnKF 與衛(wèi)星、驗(yàn)潮站等數(shù)據(jù)融合，其方法的改進(jìn)和拓展等3 個(gè)方面。在風(fēng)暴增減水特征分析方面，近30 年來風(fēng)暴潮理論[17]沒有重大突破，也基本能夠滿足對(duì)風(fēng)暴潮物理機(jī)制的解釋，但在不同臺(tái)風(fēng)影響下，長(zhǎng)江口區(qū)域風(fēng)暴增減水變化趨勢(shì)不一[18]，人們關(guān)注的近岸海域風(fēng)暴潮的時(shí)空分布及規(guī)律仍以驗(yàn)潮站、臺(tái)風(fēng)等數(shù)據(jù)統(tǒng)計(jì)分析和風(fēng)暴潮數(shù)值模擬為主[18-19]，通過風(fēng)暴增減水?dāng)?shù)據(jù)同化反演并進(jìn)行分析鮮有報(bào)道。

本文選擇201810 號(hào)臺(tái)風(fēng)“安比”登陸上海過程，采用EnKF 方法，對(duì)不同風(fēng)暴潮資料進(jìn)行同化融合，并進(jìn)行結(jié)果驗(yàn)證和時(shí)空分布特征分析。同化中增加了Schur 乘積算法[20-21]，以消除運(yùn)算中產(chǎn)生的濾波發(fā)散和偽相關(guān)問題。以期為風(fēng)暴潮機(jī)理分析、風(fēng)暴增減水再分析、風(fēng)暴潮監(jiān)測(cè)預(yù)報(bào)等研究提供技術(shù)支撐，為海洋防災(zāi)減災(zāi)工作提供服務(wù)保障。

2 數(shù)據(jù)和方法

2.1 研究區(qū)域與數(shù)據(jù)來源

研究區(qū)域?yàn)槭茱L(fēng)暴潮影響明顯的長(zhǎng)江口內(nèi)、口外和杭州灣北海域。潮位監(jiān)測(cè)站包括9 個(gè)驗(yàn)潮站（崇明南門、堡鎮(zhèn)、高橋、吳淞、蘆潮港、金山嘴、大戢山、灘滸島和佘山），其中，高橋、大戢山和灘滸島數(shù)據(jù)用于驗(yàn)證，不參與同化，其余6 個(gè)站為參與同化的站，見圖1。

為計(jì)算風(fēng)暴增減水值，搜集了2016?2018 年9 個(gè)海洋站的水位資料，資料符合《海洋觀測(cè)規(guī)范（GB/T14914—2018）》規(guī)定。

模型計(jì)算所需的風(fēng)、氣壓強(qiáng)迫場(chǎng)資料來自上海市氣象臺(tái)訂正的歐洲中期天氣預(yù)報(bào)中心（European Centre for Medium-Range Weather Forecast，ECMWF）預(yù)報(bào)產(chǎn)品，預(yù)報(bào)時(shí)長(zhǎng)為72 h，時(shí)間分辨率為3 h，空間分辨率為0.125°×0.125°，覆蓋區(qū)域0°～55°N，70°～140°E。徑流數(shù)據(jù)來自長(zhǎng)江水文網(wǎng)（www.cjh.com.cn）大通站流量數(shù)據(jù)。岸線地形資料來源為美國(guó)國(guó)家地理數(shù)據(jù)中心（National Geophysical Data Center，NGDC）的分辨率為1′的ETOPO1 全球地形數(shù)據(jù)和精度為1∶250 000的岸線數(shù)據(jù)、航海保障部最新海圖水深和岸線資料、長(zhǎng)江口航道管理局實(shí)測(cè)地形資料。開邊界調(diào)和常數(shù)基于全球潮汐模型TPXO8[22]和中國(guó)海洋圖集資料[23]調(diào)整得到。

圖1 數(shù)據(jù)同化區(qū)域、監(jiān)測(cè)站點(diǎn)分布和部分臺(tái)風(fēng)路徑Fig.1 Data assimilation region,monitoring stations and parts of typhoon track

臺(tái)風(fēng)是“安比”1949 年以來第3 個(gè)直接登陸上海且影響較大的臺(tái)風(fēng)[24-25]。于2018 年7 月22 日12:30前后在上海市崇明島沿海登陸，登陸時(shí)中心附近最大風(fēng)力達(dá)10 級(jí)（28 m/s），中心最低氣壓為982 hPa，部分臺(tái)風(fēng)中心路徑（圖1）來自中央氣象臺(tái)臺(tái)風(fēng)網(wǎng)（http://typhoon.nmc.cn/web.html）。

2.2 研究方法

2.2.1 天文潮調(diào)和分析

風(fēng)暴增減水由觀測(cè)水位（總水位）減去天文潮水位得到，其中，天文潮采用調(diào)和分析法[26]計(jì)算，公式表示為

式中，Zt表示天文潮位；A0表示平均水位；fi表示分潮i的交點(diǎn)因子；Hi表示分潮i振幅；qi表示分潮i角速度；(V0+u)i表示分潮i初相角；gi表示分潮i遲角；fi和(V0+u)i與時(shí)間t有關(guān)；n表示分潮數(shù)目。采用最小二乘法[26]分別對(duì)9 個(gè)站3 年逐時(shí)的實(shí)測(cè)潮位逼近得到68 個(gè)分潮調(diào)和常數(shù)Hi和gi，進(jìn)而實(shí)現(xiàn)臺(tái)風(fēng)“安比”過境期間天文潮位計(jì)算。

2.2.2 風(fēng)暴潮數(shù)值模擬

FVCOM（Finite-Volume Coastal Ocean Model）[27]采用有限體積數(shù)值離散方法求解三維水動(dòng)力原始控制方程組，運(yùn)用不規(guī)則三角形網(wǎng)格和干濕網(wǎng)格處理技術(shù)[28]，可以精細(xì)模擬復(fù)雜邊界和地形，在近岸海域的暴風(fēng)潮數(shù)值預(yù)報(bào)中有明顯優(yōu)勢(shì)[29]。

本文利用FVCOM 對(duì)長(zhǎng)江口風(fēng)暴潮進(jìn)行模擬，計(jì)算區(qū)域?yàn)?8°21′～34°17′N，117°33′～124°30′E，水平網(wǎng)格分辨率最高達(dá)到100 m，垂向 σ坐標(biāo)均分為10 層（圖2）。外模計(jì)算的時(shí)間步長(zhǎng)為10 s，內(nèi)外模時(shí)間步長(zhǎng)比為10∶1，開邊界潮位設(shè)置13 個(gè)分潮調(diào)和常數(shù)（M2、S2、K1、O1、N2、K2、P1、Q1、Mf、Mn、M4、MS4、MN4）。風(fēng)暴增減水場(chǎng)等于加風(fēng)場(chǎng)模擬的總水位減去不加風(fēng)場(chǎng)模擬的天文潮位，結(jié)果中包含了天文潮和風(fēng)暴潮的非線性耦合效應(yīng)。

2.2.3 集合卡爾曼濾波同化方法

圖2 長(zhǎng)江口海域計(jì)算網(wǎng)格和水深Fig.2 Computing grids and water depth in the Changjiang River Estuary

采用雙線性插值方法將數(shù)值模擬值和觀測(cè)值插值到研究區(qū)域網(wǎng)格（0.008°×0.008°）進(jìn)行空間分辨率匹配；利用蒙特卡洛隨機(jī)擾動(dòng)法（Monte Carlo Forecast，MCF）[30-31]生成狀態(tài)變量和觀測(cè)值隨機(jī)變量集合。狀態(tài)變量為FVCOM 計(jì)算的逐時(shí)風(fēng)暴增減水場(chǎng)，觀測(cè)值為6 個(gè)海洋站的逐時(shí)增減水值。狀態(tài)變量預(yù)測(cè)和觀測(cè)集合公式為式中，Xt,i、Xt?1,i為t時(shí)刻、t?1 時(shí)刻第i個(gè)樣本的狀態(tài)變量；F為狀態(tài)轉(zhuǎn)移矩陣，表示前后狀態(tài)值之間的關(guān)系；Wt?1為t?1 時(shí)刻的系統(tǒng)噪聲；Yt,i為t時(shí)刻第i樣本的觀測(cè)值；H為觀測(cè)矩陣，表示測(cè)量值與狀態(tài)值之間的關(guān)系，由于測(cè)量了狀態(tài)中的一個(gè)變量，且同一時(shí)刻風(fēng)暴潮測(cè)量值與“真實(shí)”狀態(tài)值之間接近，故取值為1；Vt,i為t時(shí)刻的觀測(cè)噪聲；系統(tǒng)噪聲和觀測(cè)噪聲的方差為Q和R。

在狀態(tài)變量的誤差協(xié)方差矩陣計(jì)算中，采用Schur乘積算法[20-21]對(duì)其進(jìn)行局地化[32]以提高估計(jì)質(zhì)量，逐步濾除距觀測(cè)點(diǎn)較遠(yuǎn)且不敏感的網(wǎng)格上的偽相關(guān)，減輕求逆矩陣不滿秩的問題。誤差協(xié)方差矩陣和Schur 乘積算法公式為

式中，Xt為狀態(tài)變量預(yù)測(cè)集合的平均值；n為樣本數(shù)，取值為200，由下文討論給出；Et為狀態(tài)變量預(yù)測(cè)集合中的擾動(dòng)值；Pt為狀態(tài)變量預(yù)測(cè)協(xié)方差矩陣；T 為矩陣轉(zhuǎn)置運(yùn)算；Ω為局地相關(guān)函數(shù)，是基于距離的相關(guān)系數(shù)矩陣；b為狀態(tài)值到觀測(cè)點(diǎn)的距離；a為Schur 半徑，表示與觀測(cè)點(diǎn)數(shù)據(jù)相關(guān)長(zhǎng)度或區(qū)域的大小，a取0.8°，由下文討論給出。

根據(jù)觀測(cè)值和預(yù)測(cè)值，結(jié)合預(yù)測(cè)協(xié)方差Pt，計(jì)算卡爾曼濾波增益Kt，得到更新后的狀態(tài)變量Xt,i和Pt，并以此循環(huán)上述步驟?？柭鼮V波增益、狀態(tài)變量和誤差協(xié)方差矩陣更新公式為

式中，Kt為卡爾曼增益矩陣；I為單位矩陣；為更新后的狀態(tài)變量；為更新后的誤差協(xié)方差矩陣。

3 數(shù)據(jù)同化檢驗(yàn)

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

本文選擇1949 年以來第3 個(gè)直接登陸上海的臺(tái)風(fēng)“安比”風(fēng)暴潮過程，將2018 年7 月21 日00 時(shí)至23 日23 時(shí)海洋站和FVCOM 數(shù)值模擬的風(fēng)暴潮數(shù)據(jù)進(jìn)行同化。為檢驗(yàn)同化效果，設(shè)計(jì)了2 組實(shí)驗(yàn)，見表1。

方案A 為控制實(shí)驗(yàn)，不同化任何觀測(cè)資料，在FVCOM 處于穩(wěn)定預(yù)報(bào)狀態(tài)時(shí)，從2018 年7 月20 日15 時(shí)（世界時(shí)）熱啟動(dòng)進(jìn)行1 h 預(yù)報(bào)，得到16 時(shí)（世界時(shí)）預(yù)報(bào)結(jié)果，依次類推至2018 年7 月23 日15 時(shí)（世界時(shí)）?？刂茖?shí)驗(yàn)還包括9 個(gè)驗(yàn)潮站實(shí)測(cè)水位調(diào)和分析計(jì)算風(fēng)暴潮的部分。

方案B 為同化實(shí)驗(yàn)，在2018 年7 月20 日16 時(shí)（世界時(shí)）的預(yù)報(bào)場(chǎng)的基礎(chǔ)上，通過EnKF 同化2018年7 月20 日16 時(shí)（世界時(shí)）的6 個(gè)站的實(shí)測(cè)水位計(jì)算得到的風(fēng)暴增減水，得到2018 年7 月20 日16 時(shí)（世界時(shí)）分析場(chǎng)，將其作為FVCOM 模型新的初始場(chǎng)進(jìn)行1 h 預(yù)報(bào)得到17 時(shí)（世界時(shí)）的預(yù)報(bào)場(chǎng)，依次類推至2018 年7 月23 日15 時(shí)（世界時(shí)）。

3.2 同化效果檢驗(yàn)

圖3 給出了6 個(gè)站實(shí)測(cè)水位計(jì)算的風(fēng)暴增減水序列（實(shí)測(cè)）、6 個(gè)站數(shù)值模擬的風(fēng)暴增減水序列（模型）和6 個(gè)站同化的風(fēng)暴增減水序列（同化）對(duì)比結(jié)果。模型與實(shí)測(cè)對(duì)比可以發(fā)現(xiàn)，模型已能表現(xiàn)出風(fēng)暴增水過程的大體特征，長(zhǎng)江口內(nèi)的崇明南門站、堡鎮(zhèn)站和吳淞站增水轉(zhuǎn)減水的時(shí)刻基本相同，但模型并未表現(xiàn)出實(shí)際增水的“雙峰”型特點(diǎn)和減水后“回彈”趨勢(shì)且增水高度低于實(shí)測(cè)；長(zhǎng)江口門處的佘山站，模型最大減水與實(shí)際較一致，但最大增水誤差較大；處于杭州灣北岸海域的蘆潮港站和金山嘴站，模型和實(shí)測(cè)的最大增水較一致，對(duì)減水的刻畫不夠準(zhǔn)確。盡管模型與實(shí)測(cè)給出的風(fēng)暴增水存在差異，但經(jīng)過同化后，其結(jié)果和實(shí)測(cè)符合的很好，具有較高的一致性。為了進(jìn)一步分析同化結(jié)果的精度，本文統(tǒng)計(jì)了它們之間的平均絕對(duì)誤差（Mean Absolute Error，MAE）和均方根誤差（Root Mean Square Error，RMSE）?？刂茖?shí)驗(yàn)中，實(shí)測(cè)與模型的MAE 為0.16 m，RMSE 為0.20 m；同化后，實(shí)測(cè)與同化的MAE 為0.06 m，比控制實(shí)驗(yàn)降低了0.10 m，RMSE 為0.07 m，與控制實(shí)驗(yàn)相比準(zhǔn)確度提高了65%。

表1 同化檢驗(yàn)實(shí)驗(yàn)設(shè)計(jì)Table 1 Experimental design of assimilation verification

圖3 同化實(shí)驗(yàn)的6 個(gè)站風(fēng)暴增減水時(shí)間序列Fig.3 Storm surge time series of 6 stations participating in assimilation

圖4 展示了未參與同化的3 個(gè)站點(diǎn)實(shí)測(cè)、模型和同化計(jì)算的風(fēng)暴增減水對(duì)比結(jié)果。由圖6 可見，模型與實(shí)測(cè)的風(fēng)暴增減水相比仍存在一定差距，長(zhǎng)江口內(nèi)高橋站的特征對(duì)比基本同崇明南門站、堡鎮(zhèn)站和吳淞站，長(zhǎng)江口門處的大戢山站和杭州灣海域的灘滸島站特征對(duì)比分別同佘山站、蘆潮港站，兩者M(jìn)AE 和RMSE 分別為0.16 m、0.19 m。同化與實(shí)測(cè)相比，同化得到的結(jié)果較為準(zhǔn)確，二者僅存在細(xì)微的偏差，其MAE 為0.07 m，RMSE 為0.09 m，說明了同化結(jié)果有效性和同化方法的正確性。

圖5 為同化實(shí)驗(yàn)中集合離散度和集合均方根誤差的對(duì)比。兩者是否接近是檢驗(yàn)集合預(yù)報(bào)系統(tǒng)可靠程度的重要指標(biāo)[34]。集合離散度是同化實(shí)驗(yàn)中集合成員與集合平均值之間的標(biāo)準(zhǔn)差（Standard Deviation，SD），集合均方根誤差則是與觀測(cè)值的差異[35]。除了集合離散度與均方根誤差保持一致外，集合離散度還應(yīng)較大，以保證真實(shí)狀態(tài)包含在集合中[36]，過低的離散度會(huì)使集合過度自信，易導(dǎo)致濾波發(fā)散[37]。本同化實(shí)驗(yàn)中，集合離散度在0.04～0.22 m 之間，平均值為0.10 m，均方根誤差在0.05～0.25 m 之間，平均值為0.11 m，集合離散度與均方根誤差比值為0.90。在圖5中，集合離散度整體上略低于均方根誤差，可能由于數(shù)值模擬的結(jié)果偏低所致，但二者已展現(xiàn)出較好的一致性，進(jìn)一步驗(yàn)證了EnKF 的可靠性。

圖4 未參加同化的3 個(gè)站風(fēng)暴增減水時(shí)間序列Fig.4 Storm surge time series of 3 stations not participating in assimilation

圖5 集合離散度和均方根誤差對(duì)比Fig.5 The comparison of set dispersion and root mean square error

圖6 和圖7 分別為典型時(shí)刻（臺(tái)風(fēng)登陸前7 月22 日12 時(shí)）數(shù)值模擬的風(fēng)暴增減水場(chǎng)、同化實(shí)驗(yàn)的增減水場(chǎng)。沿著臺(tái)風(fēng)路徑前進(jìn)方向，數(shù)值模擬和同化實(shí)驗(yàn)的增水場(chǎng)均呈現(xiàn)增水值逐漸升高趨勢(shì)，數(shù)值模擬的增減水范圍為?0.15～0.12 m，同化實(shí)驗(yàn)的增減水范圍為?0.07～0.53 m；通過典型時(shí)刻的9 個(gè)站增減水對(duì)比發(fā)現(xiàn)，實(shí)測(cè)水位計(jì)算的增減水與數(shù)值模擬的增減水之間MAE 為0.26 m，RMSE 為0.31 m；同化后，實(shí)測(cè)水位計(jì)算的與同化實(shí)驗(yàn)的MAE 為0.10 m，比不同化的降低了0.16 m，RMSE 為0.12 m，與不同化的相比降低了0.19 m?？梢?，同化作用明顯（表2）。另外，典型時(shí)刻集合離散度為0.17 m，均方根誤差0.19 m，兩者較一致。

4 風(fēng)暴潮特征分析

地形是控制和制約增減水分布的重要條件，風(fēng)和氣壓變化是造成長(zhǎng)江口增減水的主要?jiǎng)恿σ蜃覽38]。本文按照長(zhǎng)江口特定的地形條件，劃分長(zhǎng)江口斷面、杭州灣北斷面，分析同化后的臺(tái)風(fēng)“安比”風(fēng)暴潮特征；結(jié)合臺(tái)風(fēng)“安比”登陸前風(fēng)和氣壓狀況，進(jìn)一步分析口內(nèi)“雙峰”增水、臺(tái)風(fēng)眼區(qū)域增水和增水鋒面特征。

圖8 和圖9 分別為臺(tái)風(fēng)“安比”登陸前長(zhǎng)江口斷面、杭州灣北斷面各站同化后的增水值。長(zhǎng)江口斷面增水值在0.33～0.40 m 之間，東西方向的增水梯度不是很明顯，崇明南門站增水最大，高橋站增水最小?？赡苡捎诘剞D(zhuǎn)流的影響，長(zhǎng)江口南支北岸的增水略高于南岸，崇明南門站、堡鎮(zhèn)站增水高于吳淞站和高橋站。佘山站可能受臺(tái)風(fēng)眼影響，增水略有升高。杭州灣北斷面增水值在0.04～0.13 m 之間，由于處于臺(tái)風(fēng)第三象限，增水從灣內(nèi)到灣外呈遞增趨勢(shì)，且梯度較明顯，金山嘴站增水最小為0.04 m，蘆潮港和大戢山站增水最大，為0.13 m。由于東南向風(fēng)將水體吹送口內(nèi)，西北風(fēng)將水體輸運(yùn)灣外，長(zhǎng)江口斷面增水明顯高于杭州灣北斷面。

圖6 臺(tái)風(fēng)登陸前數(shù)值模擬的增減水場(chǎng)Fig.6 Storm surge field calculated by model before the typhoon landing

圖7 臺(tái)風(fēng)登陸前同化實(shí)驗(yàn)的增減水場(chǎng)Fig.7 Storm surge field of assimilation before the typhoon landing

表2 典型時(shí)刻9 個(gè)海洋站增減水對(duì)比Table 2 Storm surge comparisons of the 9 stations as typical moment

圖3 和圖4 中長(zhǎng)江口內(nèi)崇明南門、堡鎮(zhèn)、吳淞和高橋站增水曲線呈現(xiàn)“雙峰”特征，即在最大增水前后較短時(shí)間（約12 h）內(nèi)出現(xiàn)另一個(gè)較大的增水峰，兩峰高度基本一致。這與張文舟等[39]發(fā)現(xiàn)的福建沿海部分站點(diǎn)雙增水峰現(xiàn)象較相同，他認(rèn)為在臺(tái)風(fēng)橫穿臺(tái)灣海峽正面登陸福建時(shí)，臺(tái)灣海峽特殊地形和臺(tái)風(fēng)大風(fēng)區(qū)中心區(qū)域的經(jīng)過共同造成。根據(jù)“水量平衡原理”[40]，在風(fēng)暴潮“初振”階段，向口內(nèi)輸入的水量不斷增加，由于風(fēng)的頂托、降雨等影響，總向海輸出的增水量較小，因此形成了風(fēng)暴增水的“前峰”。臺(tái)風(fēng)登陸前，強(qiáng)烈的氣旋風(fēng)應(yīng)力使得水量繼續(xù)向上游輸送，導(dǎo)致口內(nèi)增水突然升高，由此產(chǎn)生了增水的“后峰”。

圖8 2018 年7 月22 日12:00 時(shí)長(zhǎng)江口斷面各站增水值Fig.8 Storm surge of the Changjiang River Estuary section stations at 12:00 July 22,2018

圖9 2018 年7 月22 日12:00 時(shí)杭州灣北斷面各站增水值Fig.9 Storm surge of the Hangzhou Bay section stations at 12:00 July 22,2018

在臺(tái)風(fēng)登陸前，各站增水曲線均出現(xiàn)大小不一的波峰，如長(zhǎng)江口內(nèi)4 個(gè)站的“后峰”，蘆潮港、金山嘴、灘滸等站在7 月22 日12 時(shí)的“小峰”?？梢?，臺(tái)風(fēng)登陸會(huì)對(duì)整個(gè)上海海域造成不同程度的增水。臺(tái)風(fēng)登陸前，佘山站處于臺(tái)風(fēng)眼附近位置，距離臺(tái)風(fēng)中心約21 km（圖6）。按靜力學(xué)計(jì)算，氣壓降低1 hPa，水位約升高1 cm[38]。臺(tái)風(fēng)中心氣壓值為980 hPa，以平均氣壓1 013 hPa 計(jì)算，約有33 cm 增水，這與佘山站增水36 cm 接近。圖7 中斷面a 處于臺(tái)風(fēng)眼附近位置，其增水梯度相對(duì)較大，進(jìn)一步表明臺(tái)風(fēng)眼對(duì)風(fēng)暴增水的影響明顯。

圖7 中區(qū)域b 是杭州灣梯度增水和長(zhǎng)江口梯度增水的匯合鋒面，由于杭州灣梯度增水由西南向東北堆積，長(zhǎng)江口梯度增水自東南向西北，此處存在增水“切變”，風(fēng)暴增水呈現(xiàn)渦的特征，對(duì)于風(fēng)暴增水研究有一定意義。囿于篇幅，其他時(shí)刻風(fēng)暴潮特征未盡分析。

5 討論

集合卡爾曼濾波同化質(zhì)量與使用的集合大小有關(guān)[41]，集合數(shù)目越多越可真實(shí)地描述系統(tǒng)狀態(tài)的空間分布，但集合數(shù)目過多會(huì)增加系統(tǒng)運(yùn)行的時(shí)間。目前，還沒有很好的方法來確定最合適的集合樣本數(shù)[42]，在大多數(shù)EnKF 的應(yīng)用中，典型的集合大小數(shù)量是100 左右[41]。本文通過固定4 組不同的Schur 半徑a，分別選擇集合樣本數(shù)n為50、100、200、400、500，在臺(tái)風(fēng)登陸前進(jìn)行EnKF 同化，比較實(shí)測(cè)與同化之間的MAE 和RMSE 以及集合離散度和集合均方根誤差，討論集合樣本數(shù)和Schur 半徑對(duì)同化的影響，并給出適合的集合樣本數(shù)和Schur 半徑。

從圖10 可以看出，集合樣本數(shù)量增加，對(duì)同化結(jié)果有很顯著的積極影響，當(dāng)集合數(shù)量繼續(xù)增大時(shí)，誤差基本維持，同化結(jié)果卻沒有顯著提升。如在a=0.6°時(shí)，集合樣本n≥200，MAE 和RMSE 維持在0.07 m、0.08 m 左右；在a為1.0°和1.4°時(shí)，集合樣本n≥400，MAE 和RMSE 維 持在0.07～0.10 m、0.09～0.12 m。但在a=0.2°時(shí)，趨勢(shì)卻相反，這可能與過小的Schur 半徑有關(guān)，對(duì)于較小的Schur 半徑，較小的集合樣本數(shù)有更好的同化結(jié)果[43]。另外，隨著集合樣本數(shù)增加，集合離散度與集合均方根誤差越趨于一致，同化系統(tǒng)可靠性加強(qiáng)。在a=0.6°和a=1.0°時(shí)，集合樣本數(shù)從50 增加到500 時(shí)，集合離散度與集合均方根誤差的比值分別從0.77 增加到0.95，從0.77 增加到0.94；在a=1.4°時(shí)，集合離散度與集合均方根誤差的比值從0.74 增加到0.89；在a=0.2°時(shí)，兩者一致性更高。

統(tǒng)計(jì)同化的單個(gè)步長(zhǎng)時(shí)間，n=200 時(shí)為24 min，n=400 時(shí)為99 min，集合樣本數(shù)增加1 倍，計(jì)算成本增加了3 倍多。除了計(jì)算效率外，集合樣本數(shù)的選擇還需考慮到實(shí)測(cè)和模型計(jì)算的步長(zhǎng)等。由于本文的計(jì)算步長(zhǎng)為60 min，再綜合上述MAE 和RMSE 分析以及集合離散度與集合均方根誤差的比值，故風(fēng)暴潮同化的集合數(shù)在100～200。

隨著Schur 半徑的增大，即相關(guān)長(zhǎng)度的增大，受觀測(cè)數(shù)據(jù)影響的區(qū)域范圍越來越大，同化過程就會(huì)保持越來越多的信息量，進(jìn)而會(huì)產(chǎn)生越來越大的偽相關(guān)性。如在集合樣本數(shù)為200 時(shí)，對(duì)比a為0.6°、1.0°和1.4°發(fā) 現(xiàn)，MAE 和RMSE 從0.07 m 和0.08 m 增加到0.11 m 和0.13 m，再到0.13 m 和0.16 m；從集合離散度與集合均方根誤差的大小來看，當(dāng)Schur 半徑從0.6°增大到1.4°時(shí)，兩條曲線越來越分開，在集合樣本數(shù)為200 時(shí)，隨著Schur 半徑的增加，比值從0.89 減為0.85，再到0.81，系統(tǒng)穩(wěn)定性有所降低。減小Schur 半徑，可以提高同化的效果，但Schur 半徑過小會(huì)導(dǎo)致部分信息丟失，使得風(fēng)暴增水場(chǎng)出現(xiàn)一定誤差，如a=0.2°時(shí)，集合樣本數(shù)達(dá)到了100 時(shí)，誤差已超出預(yù)期，并不能提高估計(jì)的質(zhì)量。綜合上述分析，當(dāng)允許誤差在0.10 m 以下時(shí)，風(fēng)暴增水同化的Schur 半徑可在0.6°～1.0°。實(shí)際上，在協(xié)方差局地化公式（7）中，Schur 半徑的選擇還可能受海洋要素類型、觀測(cè)點(diǎn)位置、網(wǎng)格大小等多種因素的制約，Schur 半徑與這些因素的關(guān)系還需要進(jìn)一步研究。

圖10 不同集合樣本數(shù)及Schur 半徑同化誤差Fig.10 The assimilation results error of different set samples and Schur radius

6 結(jié)論

本文針對(duì)實(shí)測(cè)計(jì)算和數(shù)值模擬的風(fēng)暴潮尋優(yōu)問題，以典型臺(tái)風(fēng)“安比”登陸上海的風(fēng)暴增減水為研究對(duì)象，通過協(xié)方差局地化的集合卡爾曼濾波方法進(jìn)行風(fēng)暴潮數(shù)據(jù)同化融合與特征分析。得到結(jié)論如下：

（1）同化獲得了逐72 h 的風(fēng)暴增減水的最優(yōu)解，同化與實(shí)測(cè)計(jì)算的增減水對(duì)比，MAE 為0.07 m，RMSE 為0.09 m，集合離散度與均方根誤差比值為0.90，同化效果較好且可信，為風(fēng)暴潮、數(shù)值模擬結(jié)果訂正等研究提供了新思路；（2）同化后的風(fēng)暴增減水能夠較好地刻畫臺(tái)風(fēng)登陸前雙峰增水、臺(tái)風(fēng)眼增水、增水鋒面等特征，對(duì)風(fēng)暴潮研究有一定意義；（3）協(xié)方差局地化技術(shù)可以提高風(fēng)暴潮估計(jì)質(zhì)量，降低濾波發(fā)散和偽相關(guān)性，減輕求逆矩陣不滿秩的問題；（4）集合樣本數(shù)和Schur 半徑設(shè)置對(duì)同化效果影響較大，降低相關(guān)長(zhǎng)度可以提高同化的效果，但Schur 半徑過小會(huì)導(dǎo)致部分信息丟失準(zhǔn)確度降低。