丁駿,呂忻,姚雅倩,孟鑫,姜雪敏,吳旭云,葛建忠
(1.上海市海洋監(jiān)測預報中心,上海 200062;2.國家海洋局東海預報中心,上海 200136;3.華東師范大學 河口海岸學國家重點實驗室,上海 200062)
風暴潮是世界上最嚴重的海洋災害之一[1-2],也是上海城市安全的主要威脅之一,風暴增減水[3]是一種非周期性的對眾多因素敏感的又備受關注的復雜現(xiàn)象。目前,國內外廣泛采用海洋站和衛(wèi)星等監(jiān)測到的海面高度與正常天文潮的代數(shù)差的方法來計算風暴增減水,即風暴潮大小[4-6]。由于風暴潮與天文潮的非線性耦合效應存在[7],實際增減水并非兩者簡單的線性相減。受海洋站空間分辨率限制,又常常以點代線、以點代面評估近岸海域風暴潮大小,其誤差不斷累積。雖然衛(wèi)星高度計運用到風暴潮觀測的研究[8]已經取得技術突破,但其時空分辨率和精度仍需提升。隨著理論技術發(fā)展與實測數(shù)據(jù)的豐富,風暴潮數(shù)值預報[7]的準確率和時效性進一步提高,可以直觀反映風暴增減水、漫灘過程,然而由于數(shù)值模式反映的物理過程是有限的,初始場和強迫場是近似的,模式參數(shù)是經驗的,預測結果不可避免存在誤差。在實踐中,用實測潮位分離法計算的增水值通常偏大,該值可達15~30 cm,而數(shù)值模擬的增水值往往偏小[9]。
集合卡爾曼濾波(Ensemble Kalman Filter,EnKF)[10-11]可以在受噪聲影響的測量值與數(shù)值模擬的預測值之間尋找最優(yōu)狀態(tài)估計。它基于集合預報思想,源于卡爾曼濾波理論和蒙特卡洛估計方法[12],克服了卡爾曼濾波技術計算量大、只適用線性模式的缺點。由Evensen[10]在1994 年提出并應用到海洋同化領域后,日漸成為國際上資料同化研究的熱點。在風暴潮資料同化應用方面,Butler 等[13]和Altaf 等[14]提出了一種基于EnKF 和ADCIRC(Advanced Circulation)風暴潮模型耦合的風暴潮預報方法,Colle 等[15]運用En-KF 對11 個風場資料同化驅動ADCIRC 模型來模擬“桑迪”颶風風暴潮,Etala 等[16]展現(xiàn)了局地集合變換卡爾曼濾波(Local Ensemble Transform Kalman Filter,LETKF)在衛(wèi)星遙感與驗潮站數(shù)據(jù)同化方面的優(yōu)勢,這些研究主要集中在EnKF 與風暴潮模型、氣象模型的耦合同化,EnKF 與衛(wèi)星、驗潮站等數(shù)據(jù)融合,其方法的改進和拓展等3 個方面。在風暴增減水特征分析方面,近30 年來風暴潮理論[17]沒有重大突破,也基本能夠滿足對風暴潮物理機制的解釋,但在不同臺風影響下,長江口區(qū)域風暴增減水變化趨勢不一[18],人們關注的近岸海域風暴潮的時空分布及規(guī)律仍以驗潮站、臺風等數(shù)據(jù)統(tǒng)計分析和風暴潮數(shù)值模擬為主[18-19],通過風暴增減水數(shù)據(jù)同化反演并進行分析鮮有報道。
本文選擇201810 號臺風“安比”登陸上海過程,采用EnKF 方法,對不同風暴潮資料進行同化融合,并進行結果驗證和時空分布特征分析。同化中增加了Schur 乘積算法[20-21],以消除運算中產生的濾波發(fā)散和偽相關問題。以期為風暴潮機理分析、風暴增減水再分析、風暴潮監(jiān)測預報等研究提供技術支撐,為海洋防災減災工作提供服務保障。
研究區(qū)域為受風暴潮影響明顯的長江口內、口外和杭州灣北海域。潮位監(jiān)測站包括9 個驗潮站(崇明南門、堡鎮(zhèn)、高橋、吳淞、蘆潮港、金山嘴、大戢山、灘滸島和佘山),其中,高橋、大戢山和灘滸島數(shù)據(jù)用于驗證,不參與同化,其余6 個站為參與同化的站,見圖1。
為計算風暴增減水值,搜集了2016?2018 年9 個海洋站的水位資料,資料符合《海洋觀測規(guī)范(GB/T14914—2018)》規(guī)定。
模型計算所需的風、氣壓強迫場資料來自上海市氣象臺訂正的歐洲中期天氣預報中心(European Centre for Medium-Range Weather Forecast,ECMWF)預報產品,預報時長為72 h,時間分辨率為3 h,空間分辨率為0.125°×0.125°,覆蓋區(qū)域0°~55°N,70°~140°E。徑流數(shù)據(jù)來自長江水文網(wǎng)(www.cjh.com.cn)大通站流量數(shù)據(jù)。岸線地形資料來源為美國國家地理數(shù)據(jù)中心(National Geophysical Data Center,NGDC)的分辨率為1′的ETOPO1 全球地形數(shù)據(jù)和精度為1∶250 000的岸線數(shù)據(jù)、航海保障部最新海圖水深和岸線資料、長江口航道管理局實測地形資料。開邊界調和常數(shù)基于全球潮汐模型TPXO8[22]和中國海洋圖集資料[23]調整得到。
圖1 數(shù)據(jù)同化區(qū)域、監(jiān)測站點分布和部分臺風路徑Fig.1 Data assimilation region,monitoring stations and parts of typhoon track
臺風是“安比”1949 年以來第3 個直接登陸上海且影響較大的臺風[24-25]。于2018 年7 月22 日12:30前后在上海市崇明島沿海登陸,登陸時中心附近最大風力達10 級(28 m/s),中心最低氣壓為982 hPa,部分臺風中心路徑(圖1)來自中央氣象臺臺風網(wǎng)(http://typhoon.nmc.cn/web.html)。
2.2.1 天文潮調和分析
風暴增減水由觀測水位(總水位)減去天文潮水位得到,其中,天文潮采用調和分析法[26]計算,公式表示為
式中,Zt表示天文潮位;A0表示平均水位;fi表示分潮i的交點因子;Hi表示分潮i振幅;qi表示分潮i角速度;(V0+u)i表示分潮i初相角;gi表示分潮i遲角;fi和(V0+u)i與時間t有關;n表示分潮數(shù)目。采用最小二乘法[26]分別對9 個站3 年逐時的實測潮位逼近得到68 個分潮調和常數(shù)Hi和gi,進而實現(xiàn)臺風“安比”過境期間天文潮位計算。
2.2.2 風暴潮數(shù)值模擬
FVCOM(Finite-Volume Coastal Ocean Model)[27]采用有限體積數(shù)值離散方法求解三維水動力原始控制方程組,運用不規(guī)則三角形網(wǎng)格和干濕網(wǎng)格處理技術[28],可以精細模擬復雜邊界和地形,在近岸海域的暴風潮數(shù)值預報中有明顯優(yōu)勢[29]。
本文利用FVCOM 對長江口風暴潮進行模擬,計算區(qū)域為28°21′~34°17′N,117°33′~124°30′E,水平網(wǎng)格分辨率最高達到100 m,垂向 σ坐標均分為10 層(圖2)。外模計算的時間步長為10 s,內外模時間步長比為10∶1,開邊界潮位設置13 個分潮調和常數(shù)(M2、S2、K1、O1、N2、K2、P1、Q1、Mf、Mn、M4、MS4、MN4)。風暴增減水場等于加風場模擬的總水位減去不加風場模擬的天文潮位,結果中包含了天文潮和風暴潮的非線性耦合效應。
2.2.3 集合卡爾曼濾波同化方法
圖2 長江口海域計算網(wǎng)格和水深Fig.2 Computing grids and water depth in the Changjiang River Estuary
采用雙線性插值方法將數(shù)值模擬值和觀測值插值到研究區(qū)域網(wǎng)格(0.008°×0.008°)進行空間分辨率匹配;利用蒙特卡洛隨機擾動法(Monte Carlo Forecast,MCF)[30-31]生成狀態(tài)變量和觀測值隨機變量集合。狀態(tài)變量為FVCOM 計算的逐時風暴增減水場,觀測值為6 個海洋站的逐時增減水值。狀態(tài)變量預測和觀測集合公式為式中,Xt,i、Xt?1,i為t時刻、t?1 時刻第i個樣本的狀態(tài)變量;F為狀態(tài)轉移矩陣,表示前后狀態(tài)值之間的關系;Wt?1為t?1 時刻的系統(tǒng)噪聲;Yt,i為t時刻第i樣本的觀測值;H為觀測矩陣,表示測量值與狀態(tài)值之間的關系,由于測量了狀態(tài)中的一個變量,且同一時刻風暴潮測量值與“真實”狀態(tài)值之間接近,故取值為1;Vt,i為t時刻的觀測噪聲;系統(tǒng)噪聲和觀測噪聲的方差為Q和R。
在狀態(tài)變量的誤差協(xié)方差矩陣計算中,采用Schur乘積算法[20-21]對其進行局地化[32]以提高估計質量,逐步濾除距觀測點較遠且不敏感的網(wǎng)格上的偽相關,減輕求逆矩陣不滿秩的問題。誤差協(xié)方差矩陣和Schur 乘積算法公式為
式中,Xt為狀態(tài)變量預測集合的平均值;n為樣本數(shù),取值為200,由下文討論給出;Et為狀態(tài)變量預測集合中的擾動值;Pt為狀態(tài)變量預測協(xié)方差矩陣;T 為矩陣轉置運算;Ω為局地相關函數(shù),是基于距離的相關系數(shù)矩陣;b為狀態(tài)值到觀測點的距離;a為Schur 半徑,表示與觀測點數(shù)據(jù)相關長度或區(qū)域的大小,a取0.8°,由下文討論給出。
根據(jù)觀測值和預測值,結合預測協(xié)方差Pt,計算卡爾曼濾波增益Kt,得到更新后的狀態(tài)變量Xt,i和Pt,并以此循環(huán)上述步驟??柭鼮V波增益、狀態(tài)變量和誤差協(xié)方差矩陣更新公式為
式中,Kt為卡爾曼增益矩陣;I為單位矩陣;為更新后的狀態(tài)變量;為更新后的誤差協(xié)方差矩陣。
本文選擇1949 年以來第3 個直接登陸上海的臺風“安比”風暴潮過程,將2018 年7 月21 日00 時至23 日23 時海洋站和FVCOM 數(shù)值模擬的風暴潮數(shù)據(jù)進行同化。為檢驗同化效果,設計了2 組實驗,見表1。
方案A 為控制實驗,不同化任何觀測資料,在FVCOM 處于穩(wěn)定預報狀態(tài)時,從2018 年7 月20 日15 時(世界時)熱啟動進行1 h 預報,得到16 時(世界時)預報結果,依次類推至2018 年7 月23 日15 時(世界時)??刂茖嶒炦€包括9 個驗潮站實測水位調和分析計算風暴潮的部分。
方案B 為同化實驗,在2018 年7 月20 日16 時(世界時)的預報場的基礎上,通過EnKF 同化2018年7 月20 日16 時(世界時)的6 個站的實測水位計算得到的風暴增減水,得到2018 年7 月20 日16 時(世界時)分析場,將其作為FVCOM 模型新的初始場進行1 h 預報得到17 時(世界時)的預報場,依次類推至2018 年7 月23 日15 時(世界時)。
圖3 給出了6 個站實測水位計算的風暴增減水序列(實測)、6 個站數(shù)值模擬的風暴增減水序列(模型)和6 個站同化的風暴增減水序列(同化)對比結果。模型與實測對比可以發(fā)現(xiàn),模型已能表現(xiàn)出風暴增水過程的大體特征,長江口內的崇明南門站、堡鎮(zhèn)站和吳淞站增水轉減水的時刻基本相同,但模型并未表現(xiàn)出實際增水的“雙峰”型特點和減水后“回彈”趨勢且增水高度低于實測;長江口門處的佘山站,模型最大減水與實際較一致,但最大增水誤差較大;處于杭州灣北岸海域的蘆潮港站和金山嘴站,模型和實測的最大增水較一致,對減水的刻畫不夠準確。盡管模型與實測給出的風暴增水存在差異,但經過同化后,其結果和實測符合的很好,具有較高的一致性。為了進一步分析同化結果的精度,本文統(tǒng)計了它們之間的平均絕對誤差(Mean Absolute Error,MAE)和均方根誤差(Root Mean Square Error,RMSE)??刂茖嶒炛?,實測與模型的MAE 為0.16 m,RMSE 為0.20 m;同化后,實測與同化的MAE 為0.06 m,比控制實驗降低了0.10 m,RMSE 為0.07 m,與控制實驗相比準確度提高了65%。
表1 同化檢驗實驗設計Table 1 Experimental design of assimilation verification
圖3 同化實驗的6 個站風暴增減水時間序列Fig.3 Storm surge time series of 6 stations participating in assimilation
圖4 展示了未參與同化的3 個站點實測、模型和同化計算的風暴增減水對比結果。由圖6 可見,模型與實測的風暴增減水相比仍存在一定差距,長江口內高橋站的特征對比基本同崇明南門站、堡鎮(zhèn)站和吳淞站,長江口門處的大戢山站和杭州灣海域的灘滸島站特征對比分別同佘山站、蘆潮港站,兩者MAE 和RMSE 分別為0.16 m、0.19 m。同化與實測相比,同化得到的結果較為準確,二者僅存在細微的偏差,其MAE 為0.07 m,RMSE 為0.09 m,說明了同化結果有效性和同化方法的正確性。
圖5 為同化實驗中集合離散度和集合均方根誤差的對比。兩者是否接近是檢驗集合預報系統(tǒng)可靠程度的重要指標[34]。集合離散度是同化實驗中集合成員與集合平均值之間的標準差(Standard Deviation,SD),集合均方根誤差則是與觀測值的差異[35]。除了集合離散度與均方根誤差保持一致外,集合離散度還應較大,以保證真實狀態(tài)包含在集合中[36],過低的離散度會使集合過度自信,易導致濾波發(fā)散[37]。本同化實驗中,集合離散度在0.04~0.22 m 之間,平均值為0.10 m,均方根誤差在0.05~0.25 m 之間,平均值為0.11 m,集合離散度與均方根誤差比值為0.90。在圖5中,集合離散度整體上略低于均方根誤差,可能由于數(shù)值模擬的結果偏低所致,但二者已展現(xiàn)出較好的一致性,進一步驗證了EnKF 的可靠性。
圖4 未參加同化的3 個站風暴增減水時間序列Fig.4 Storm surge time series of 3 stations not participating in assimilation
圖5 集合離散度和均方根誤差對比Fig.5 The comparison of set dispersion and root mean square error
圖6 和圖7 分別為典型時刻(臺風登陸前7 月22 日12 時)數(shù)值模擬的風暴增減水場、同化實驗的增減水場。沿著臺風路徑前進方向,數(shù)值模擬和同化實驗的增水場均呈現(xiàn)增水值逐漸升高趨勢,數(shù)值模擬的增減水范圍為?0.15~0.12 m,同化實驗的增減水范圍為?0.07~0.53 m;通過典型時刻的9 個站增減水對比發(fā)現(xiàn),實測水位計算的增減水與數(shù)值模擬的增減水之間MAE 為0.26 m,RMSE 為0.31 m;同化后,實測水位計算的與同化實驗的MAE 為0.10 m,比不同化的降低了0.16 m,RMSE 為0.12 m,與不同化的相比降低了0.19 m??梢?,同化作用明顯(表2)。另外,典型時刻集合離散度為0.17 m,均方根誤差0.19 m,兩者較一致。
地形是控制和制約增減水分布的重要條件,風和氣壓變化是造成長江口增減水的主要動力因子[38]。本文按照長江口特定的地形條件,劃分長江口斷面、杭州灣北斷面,分析同化后的臺風“安比”風暴潮特征;結合臺風“安比”登陸前風和氣壓狀況,進一步分析口內“雙峰”增水、臺風眼區(qū)域增水和增水鋒面特征。
圖8 和圖9 分別為臺風“安比”登陸前長江口斷面、杭州灣北斷面各站同化后的增水值。長江口斷面增水值在0.33~0.40 m 之間,東西方向的增水梯度不是很明顯,崇明南門站增水最大,高橋站增水最小??赡苡捎诘剞D流的影響,長江口南支北岸的增水略高于南岸,崇明南門站、堡鎮(zhèn)站增水高于吳淞站和高橋站。佘山站可能受臺風眼影響,增水略有升高。杭州灣北斷面增水值在0.04~0.13 m 之間,由于處于臺風第三象限,增水從灣內到灣外呈遞增趨勢,且梯度較明顯,金山嘴站增水最小為0.04 m,蘆潮港和大戢山站增水最大,為0.13 m。由于東南向風將水體吹送口內,西北風將水體輸運灣外,長江口斷面增水明顯高于杭州灣北斷面。
圖6 臺風登陸前數(shù)值模擬的增減水場Fig.6 Storm surge field calculated by model before the typhoon landing
圖7 臺風登陸前同化實驗的增減水場Fig.7 Storm surge field of assimilation before the typhoon landing
表2 典型時刻9 個海洋站增減水對比Table 2 Storm surge comparisons of the 9 stations as typical moment
圖3 和圖4 中長江口內崇明南門、堡鎮(zhèn)、吳淞和高橋站增水曲線呈現(xiàn)“雙峰”特征,即在最大增水前后較短時間(約12 h)內出現(xiàn)另一個較大的增水峰,兩峰高度基本一致。這與張文舟等[39]發(fā)現(xiàn)的福建沿海部分站點雙增水峰現(xiàn)象較相同,他認為在臺風橫穿臺灣海峽正面登陸福建時,臺灣海峽特殊地形和臺風大風區(qū)中心區(qū)域的經過共同造成。根據(jù)“水量平衡原理”[40],在風暴潮“初振”階段,向口內輸入的水量不斷增加,由于風的頂托、降雨等影響,總向海輸出的增水量較小,因此形成了風暴增水的“前峰”。臺風登陸前,強烈的氣旋風應力使得水量繼續(xù)向上游輸送,導致口內增水突然升高,由此產生了增水的“后峰”。
圖8 2018 年7 月22 日12:00 時長江口斷面各站增水值Fig.8 Storm surge of the Changjiang River Estuary section stations at 12:00 July 22,2018
圖9 2018 年7 月22 日12:00 時杭州灣北斷面各站增水值Fig.9 Storm surge of the Hangzhou Bay section stations at 12:00 July 22,2018
在臺風登陸前,各站增水曲線均出現(xiàn)大小不一的波峰,如長江口內4 個站的“后峰”,蘆潮港、金山嘴、灘滸等站在7 月22 日12 時的“小峰”??梢?,臺風登陸會對整個上海海域造成不同程度的增水。臺風登陸前,佘山站處于臺風眼附近位置,距離臺風中心約21 km(圖6)。按靜力學計算,氣壓降低1 hPa,水位約升高1 cm[38]。臺風中心氣壓值為980 hPa,以平均氣壓1 013 hPa 計算,約有33 cm 增水,這與佘山站增水36 cm 接近。圖7 中斷面a 處于臺風眼附近位置,其增水梯度相對較大,進一步表明臺風眼對風暴增水的影響明顯。
圖7 中區(qū)域b 是杭州灣梯度增水和長江口梯度增水的匯合鋒面,由于杭州灣梯度增水由西南向東北堆積,長江口梯度增水自東南向西北,此處存在增水“切變”,風暴增水呈現(xiàn)渦的特征,對于風暴增水研究有一定意義。囿于篇幅,其他時刻風暴潮特征未盡分析。
集合卡爾曼濾波同化質量與使用的集合大小有關[41],集合數(shù)目越多越可真實地描述系統(tǒng)狀態(tài)的空間分布,但集合數(shù)目過多會增加系統(tǒng)運行的時間。目前,還沒有很好的方法來確定最合適的集合樣本數(shù)[42],在大多數(shù)EnKF 的應用中,典型的集合大小數(shù)量是100 左右[41]。本文通過固定4 組不同的Schur 半徑a,分別選擇集合樣本數(shù)n為50、100、200、400、500,在臺風登陸前進行EnKF 同化,比較實測與同化之間的MAE 和RMSE 以及集合離散度和集合均方根誤差,討論集合樣本數(shù)和Schur 半徑對同化的影響,并給出適合的集合樣本數(shù)和Schur 半徑。
從圖10 可以看出,集合樣本數(shù)量增加,對同化結果有很顯著的積極影響,當集合數(shù)量繼續(xù)增大時,誤差基本維持,同化結果卻沒有顯著提升。如在a=0.6°時,集合樣本n≥200,MAE 和RMSE 維持在0.07 m、0.08 m 左右;在a為1.0°和1.4°時,集合樣本n≥400,MAE 和RMSE 維 持在0.07~0.10 m、0.09~0.12 m。但在a=0.2°時,趨勢卻相反,這可能與過小的Schur 半徑有關,對于較小的Schur 半徑,較小的集合樣本數(shù)有更好的同化結果[43]。另外,隨著集合樣本數(shù)增加,集合離散度與集合均方根誤差越趨于一致,同化系統(tǒng)可靠性加強。在a=0.6°和a=1.0°時,集合樣本數(shù)從50 增加到500 時,集合離散度與集合均方根誤差的比值分別從0.77 增加到0.95,從0.77 增加到0.94;在a=1.4°時,集合離散度與集合均方根誤差的比值從0.74 增加到0.89;在a=0.2°時,兩者一致性更高。
統(tǒng)計同化的單個步長時間,n=200 時為24 min,n=400 時為99 min,集合樣本數(shù)增加1 倍,計算成本增加了3 倍多。除了計算效率外,集合樣本數(shù)的選擇還需考慮到實測和模型計算的步長等。由于本文的計算步長為60 min,再綜合上述MAE 和RMSE 分析以及集合離散度與集合均方根誤差的比值,故風暴潮同化的集合數(shù)在100~200。
隨著Schur 半徑的增大,即相關長度的增大,受觀測數(shù)據(jù)影響的區(qū)域范圍越來越大,同化過程就會保持越來越多的信息量,進而會產生越來越大的偽相關性。如在集合樣本數(shù)為200 時,對比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;從集合離散度與集合均方根誤差的大小來看,當Schur 半徑從0.6°增大到1.4°時,兩條曲線越來越分開,在集合樣本數(shù)為200 時,隨著Schur 半徑的增加,比值從0.89 減為0.85,再到0.81,系統(tǒng)穩(wěn)定性有所降低。減小Schur 半徑,可以提高同化的效果,但Schur 半徑過小會導致部分信息丟失,使得風暴增水場出現(xiàn)一定誤差,如a=0.2°時,集合樣本數(shù)達到了100 時,誤差已超出預期,并不能提高估計的質量。綜合上述分析,當允許誤差在0.10 m 以下時,風暴增水同化的Schur 半徑可在0.6°~1.0°。實際上,在協(xié)方差局地化公式(7)中,Schur 半徑的選擇還可能受海洋要素類型、觀測點位置、網(wǎng)格大小等多種因素的制約,Schur 半徑與這些因素的關系還需要進一步研究。
圖10 不同集合樣本數(shù)及Schur 半徑同化誤差Fig.10 The assimilation results error of different set samples and Schur radius
本文針對實測計算和數(shù)值模擬的風暴潮尋優(yōu)問題,以典型臺風“安比”登陸上海的風暴增減水為研究對象,通過協(xié)方差局地化的集合卡爾曼濾波方法進行風暴潮數(shù)據(jù)同化融合與特征分析。得到結論如下:
(1)同化獲得了逐72 h 的風暴增減水的最優(yōu)解,同化與實測計算的增減水對比,MAE 為0.07 m,RMSE 為0.09 m,集合離散度與均方根誤差比值為0.90,同化效果較好且可信,為風暴潮、數(shù)值模擬結果訂正等研究提供了新思路;(2)同化后的風暴增減水能夠較好地刻畫臺風登陸前雙峰增水、臺風眼增水、增水鋒面等特征,對風暴潮研究有一定意義;(3)協(xié)方差局地化技術可以提高風暴潮估計質量,降低濾波發(fā)散和偽相關性,減輕求逆矩陣不滿秩的問題;(4)集合樣本數(shù)和Schur 半徑設置對同化效果影響較大,降低相關長度可以提高同化的效果,但Schur 半徑過小會導致部分信息丟失準確度降低。