基于多源数据的台风风暴潮概率预报研究：数值预报系统

郭文云，安佰超，裘诚，李铖，李丕学，葛建忠，丁平兴

（1.上海海事大学海洋科学与工程学院，上海201300；2.上海市海洋监测预报中心，上海200062；3.华东师范大学河口海岸学国家重点实验室，上海200062）

1 引言

台风引起的风暴潮灾害每年都给我国造成大量经济损失[1-2]。开展风暴潮预报是防止（减少）风暴潮灾害并且保障人民生命财产安全的有效手段。随着SPLASH（Special Program to List Amplitudes of Surges from Hurricanes）、SLOSH（Sea，Lake，and Overland Surges from Hurricanes model）和ADCIRC（The ADvanced CIRCulation model）等模型的应用，数值模型已成为风暴潮预报的重要方法[3-9]。风暴增水强度对台风路径和强度等参数很敏感，但台风路径和强度预报存在很大的不确定性。据统计，2016年中央气象台台风24 h、48 h和72 h预报的平均路径误差分别为76.2 km、147.3 km和244.7 km，平均风速误差分别为5.4 m/s、7.3 m/s和7.6 m/s[10]。这些台风预报误差给风暴潮数值预报工作带来挑战。充分考虑台风预报中的误差，进行多个台风可能路径及可能强度的风暴潮集合数值预报，从而提供更多更全面的预报信息，可以有效弥补这一不足[11-15]。对集合预报中样本的发生概率进行仔细估计，称为概率预报。

目前，多个国家的风暴潮业务化预报系统已经开始采用集合（概率）预报技术。美国国家飓风中心（National Hurricane Center）使用的P_surge风暴潮概率预报模型考虑了台风路径、台风中心附近最大风速和台风大小的不确定性，根据历史预报误差统计生成27个台风样本来进行风暴潮概率预报[16]。澳大利亚气象局（Bureau of Meteorology Australia）则基于特定方法生成1 000条可能的台风路径，然后从中随机选取200条来进行风暴潮集合预报[17]。

我国在风暴潮集合预报研究上也取得不少进展。国家海洋环境预报中心已经开发了一套风暴潮概率预报模式，并实现了业务化运行。该模式主要考虑了台风路径的不确定性，计算最大可能情况下的风暴增减水极值。利用该集合预报模式已对0814号强台风“黑格比”进行数值预报实验[18]。王培涛等[12]基于预报路径衍生偏快、偏慢、偏左和偏右4条路径对福建沿海进行了风暴潮集合预报。陈永平等[19]在5条可能路径（包括预报路径）的基础上进一步考虑了3个可能的最大风速，使用15个台风样本进行风暴潮集合计算。以上结果都显示，采用集合预报技术可以有效降低风暴潮预报误差。目前我国的研究大多还停留在集合预报上，只是通过构造不同台风样本来进行多个情景的预报，并通过分析预报增水平均值来给出预报增水幅度，在统计学上缺乏实际概率意义，并未达到真正的概率预报。事实上，不同台风样本的发生概率应该有所不同。

郭文云等[1]建立了一个台风集合构造方案。该方案可构造27个台风样本（9条可能路径和3个可能最大风速），并基于历史预报误差的统计分布来确定不同台风样本的发生概率。该集合构建方案可以较为全面地覆盖可能的情景。本文基于这一台风集合构造方案，采用中尺度天气预报（Weather Research and Forecasting，WRF）模式和有限体积海岸海洋模式（Finite Volume Coastal Ocean Model，FVCOM）建立长江口台风风暴潮集合预报系统，并以1909号台风“利奇马”为例，详细说明了概率预报提供的丰富的风暴潮概率预报信息。这些信息可以为风暴潮防灾减灾工作提供详细的数据参考。

2 长江口风暴潮数值模型系统

台风是一种局地的强非线性气象系统，其最大风速半径一般为几十公里，对风暴潮进行数值模拟需要精细地刻画台风的风场。台风中心具有快速移动性，因此需要较大的模型区域来准确模拟风暴增水过程。本文采用高分辨率台风风场叠加大尺度背景风场的方法来提供风暴潮模型所需的表面风应力输入。其中台风风场使用经典的藤田公式来描述，大尺度背景风场由三重嵌套的WRF中尺度气象模型提供，而风暴潮模型则选用国际先进的FVCOM三维模型进行两重嵌套。该模型采用无结构网格，可以精确地刻画长江口复杂的岸线和地形条件。

2.1 WRF气象模型

WRF模式具有计算高效、可移植性强、方便且易于维护等优点，适合于大气动力学研究和业务化预报等工作。为提高WRF的计算效率，同时保证计算范围足够大，以涵盖台风风暴潮的生成、发展和消亡的全过程，因此采用三重嵌套和双向数据交换的方法。模式的计算区域如图1所示。WRF大区覆盖从15°～51°N，109°～149°E的区域；中区覆盖从23°～41°N，117°～129°E的区域；而小区覆盖长江口和杭州湾沿海地区及江浙沪部分地区。各区域的分辨率分别为81 km、27 km和9 km。

图1 长江口风暴潮预报系统中的三重嵌套WRF气象模式区域

WRF模型在进行后报时采用美国国家环境预报中心（National Centers for Environmental Prediction，NCEP）提供的经过同化后的FNL（Final Operational Global Analysis）数据，空间分辨率为1°×1°，时间间隔为6 h；预报时则采用NCEP提供的全球预报产品GFS（Global Forecast System），空间分辨率为0.5°×0.5°，时间间隔为3 h。

2.2 台风模型

本文采用对称台风气压场模型。模型需要用到的台风参数包括台风位置、台风中心气压P0、台风中心附近最大风速Vmax和最大风速半径Rmax。其中最大风速半径无法在预报台风中给出，因此使用经验公式[20]进行估计：

台风气压场使用藤田公式[21]来估计：

式中：P0为台风中心气压；P∞为台风外围气压；r是计算点至气压中心的距离；R0是表征台风系统特征的参数，可根据最大风速半径调整。

气旋的风速分布根据梯度风关系[22-23]参照气压分布给出：

式中：f=2Ωsinφ为柯氏参数；Ω为地球自转角速度；φ为纬度；ρa为空气密度。台风移行产生的风场采用宫崎正卫公式[24-25]：

式中：Vx和Vy分别是台风中心移动速度的正东分量和正北分量。

2.3 风场融合

合成后的台风风场可以表示为：

W→=c1W1[-sin(ϕ+β),cos(ϕ+β)]+c2W→2（5）

式中：ϕ是计算点与台风中心连线与正东方向的夹角；β是梯度风与海面风的夹角；c1和c2是订正系数。

由对称模型风场计算所得的风场与背景风场之间的合成，根据计算点与台风中心的距离r与判别距离R1和R2的关系进行插值合成：

式中：a=(r-R1)/(R2-R1)。

表面风应力的计算采用二次率参数化形式：

式中：ρa为空气密度，取1.23 kg/m3；W→为风矢量；CD为拖曳系数，随风速大小而变化：

2.4 风暴潮模型

长江口风暴潮模型采用葛建忠等[26-27]建立的东中国海-长江口FVCOM嵌套模型（见图2）。东中国海大区模型覆盖我国渤海、黄海、东海及其邻近海域。长江口精细化模型覆盖了长江各个入海口以及杭州湾和舟山群岛。长江口上界取至潮流界江阴以上200 km左右的大通；钱塘江上界取至老盐仓附近；外海开边界东至124.5°E，北至34.5°N，南至28.5°N。网格分辨率在外海陆架海域为0.5～5 km，在长江口则精细到200 m左右。模型开边界基于葛建忠等[26]的数据调整得到，在8个主要分潮（M2、S2、K1、O1、N2、K2、P1和Q1）的基础上，加入两个长周期分潮（Mf和Mm）和3个浅水分潮（M4、MS4和MN4）。所加入的分潮都是在OSU TOPEX/POSEIDON全球海洋潮汐模式[28]的基础上进行适当调整后的结果。

2.5 模型验证

2018年7月18日20时（北京时，下同），台风“安比”（1810）在西北太平洋洋面上生成（见图3a）。2018年7月22日12时30分前后，台风“安比”在上海市崇明岛沿海登陆，登陆时中心附近最大风力达10级（28 m/s），中心最低气压为982 hPa，成为1949年以来第三个直接登陆上海的台风。2018年7月24日夜间，台风“安比”（热带风暴级）减弱变性为温带气旋，之后强度持续减弱，中央气象台于25日02时对其停止编号。

图2 长江口风暴潮模型的东中国海-长江口嵌套网格

2019年8月4日，台风“利奇马”（1909）生成（见图3b），并于2019年8月10日01时45分左右在浙江省温岭市城南镇沿海登陆，登陆时中心附近最大风力16级（52 m/s），属超强台风，是1949年以来登陆我国的第五个超强台风。随后台风“利奇马”纵穿浙江和江苏两省并移入黄海海面，又于8月11日20时50分许在山东省青岛市黄岛区沿海再次登陆，此后其移入渤海海面并不断减弱，最终于8月13日14时被中央气象台停止编号。截至2019年8月14日10时，台风“利奇马”共造成我国57人死亡，直接经济损失537.2亿元。金山嘴潮位站的水位观测数据显示，台风“利奇马”造成该站最大增水超过1.0 m。

采用本文构建的风暴潮模型对台风“安比”和“利奇马”进行后报检验，所得南门、堡镇、芦潮港和金山嘴4站的增水过程对比如图4所示。本文构建的风暴潮模型不仅可以准确地刻画台风“安比”和“利奇马”引起的增水过程，对台风过后的减水过程也能准确模拟。南门、堡镇和芦潮港的平均增水误差都小于15 cm。根据观测数据，金山嘴站在两个台风过后都存在较强余振（振幅超过40 cm），本文的模拟结果虽然也能体现余振效果，但对余振的刻画还不够精细，使得该站的平均误差超过20 cm。这可能是由于台风登陆后风场发生变形，而本文所采用的对称台风模型难以刻画这种变形。模型模拟的最大增水时间与观测基本一致，最大增水幅度比观测稍强，从预报安全性角度上是可以接受的。此外，台风“利奇马”引起的长江口增水都带有双峰特征，本文的模型也能反映出类似的特征。

3 预报实验

本文的主要目的是建立长江口概率预报系统。相对于单一预报，概率预报系统能提供更多的风暴增水预报结果。如何从这些结果中挖掘出所需的风暴增水信息，也是风暴潮概率预报的一项重要工作。

仍然以台风“利奇马”为例，选取芦潮港站最大增水发生前24 h左右（2019年8月9日14时）为预报时刻进行虚拟预报实验。采用郭文云等[1]的台风集合构造方案，先根据中国中央台、中国香港台、中国台湾台、美国台、日本台和韩国台的预报路径生成一条误差更小的分析路径（见图3中红色线），然后基于分析路径及其统计误差生成两个概率圆上的8条衍生路径（见图3），路径集合覆盖了实测路径的运移范围。类似方法可以得到预报台风最大风速集合（见表1）。其中48 h和72 h预报时效最大风速都覆盖到了实测最大风速，但24 h预报最大风速没有覆盖到实测最大风速，实测风速弱于概率偏小的预报最大风速。

图3 两个台风的实测路径和台风路径样本构建（注：台风强度标注只用于实测路径）

以上路径集合和最大风速集合相互交叉形成27个台风样本，每个台风样本的发生概率都根据分析台风的误差分布确定。对以上27个台风样本分别进行风暴潮预报，可得到27个风暴潮预报结果。

图4 两个台风过程中4站的实测和后报增水过程对比

3.1 最大和最小可能增水分布

将以上27个预报中的水位数据输出，那么对于任意时刻任意空间点，都存在27个可能的风暴水位。同时进行只考虑天文潮汐而不考虑台风过程的模拟，相减即可得到27个不同的增水场。其中最大增水值属最坏情况，最小增水值属最好情况。从预报24 h后即2019年8月10日14时的最大和最小可能增水分布图上可以看到（见图5），此时长江口4条分支和杭州湾内的最大可能增水都超过2.0 m，但最小可能增水还不到0.5 m；长江口门附近甚至还有可能出现减水情况。最大和最小可能增水之间存在巨大差异，充分说明台风预报中的路径和风速等参数的误差会使风暴潮预报增水产生巨大偏差，也体现了建立风暴潮概率预报系统的重要性。若采用确定性的风暴潮单一预报系统，以上信息都无法提供。

3.2 平均增水及其标准偏差分布

结合以上27个台风样本的风暴潮预报结果及其发生概率，可以计算不同时刻不同空间点的平均增水值及其标准差：

式中：ηi为第i个数值预报实验的增水值；pi为该数值实验所用台风样本的发生概率。

图6为24 h预报时效（即2019年8月10日14时）的平均增水分布及其标准差。图中可以看出，此时东中国海的增水并不明显，一般在0.5 m以下；沿长江口向陆以及沿杭州湾向内，平均增水值不断增大，长江口4条分支和杭州湾内的平均增水都超过1.5 m，但实测增水为1.0 m左右，这与台风预报最大风速较实际最大风速显著偏大有关。即使是概率偏弱（20%）的台风最大风速都比实测最大风速更大（见表1）。该时刻的预报增水标准偏差也呈现外海较小近岸较大的空间分布特征，在长江口内和杭州湾内，大部分区域的标准偏差大于0.5 m，表明该时刻预报增水的不确定性较大。这些信息对于风暴潮的防灾减灾工作都具有重要的参考价值。

图5 台风“利奇马”24 h预报时效的可能增水（预报时间为2019年8月9日14时；单位：m）

图6 台风“利奇马”24 h预报时效结果（预报时间为2019年8月9日14时，单位：m）

为方便对比，本文同时也将分析路径的预报增水结果（代表确定性单一预报结果）呈现于图6c。其总体增水场分布与预报平均增水分布相似，但在长江口内和杭州湾内的增水值略微更大。只根据这单一预报结果，我们无从得知该时刻增水的不确定性。由于预报台风路径存在明显误差，势必导致预报风暴潮增水的不确定性。相对于单一预报系统，概率预报系统则可以提供这些信息，充分体现了建立风暴潮概率预报系统的必要性。

3.3 概率增水值和增水概率值

风暴潮概率预报系统可以得到某个站点在未来某个时刻的多个增水值及其发生概率。将这些增水值从大到小排列并计算其累积可得到该站点在某时刻的增水概率图。图7为24 h预报时效的芦潮港站台风增水累积概率图，在所有可能情况中，该时刻芦潮港站都处于增水状态。对于某个固定的概率，其增水值可以达到多大，此为概率增水值。芦潮港站在该时刻的20%概率增水值为1.65 m（见图7a）。

风暴潮可以根据其最大增水值来划分等级，一般以1.0 m、1.5 m、2.0 m和2.5 m为界将风暴潮分为一般、中等、较大、大和特大5个等级[29]。我们同样希望得到对于某个固定的增水值，其发生概率如何，此为增水概率值。根据图7b可以确定芦潮港站（高于）1.0 m增水概率为84.9%。

依据以上对概率增水值和增水概率值的定义，我们可以根据概率预报结果得到未来某个时刻的概率增水分布及增水概率分布。图8a为2019年8月10日14时预报的20%概率增水分布。图中可以看到外海的增水值与平均增水没有太大差别，表明台风预报误差对外海的风暴增水影响较小，而在长江口内及杭州湾内，20%概率增水超过2.0 m。

图8b为1.0 m增水概率分布，即预报该时刻增水高于1.0 m的概率。在外海，这一概率都为0，表明外海增水不可能超过1.0 m；一旦进入河道或杭州湾内，由于地形缩窄，潮能汇聚，1.0 m增水概率迅速增加，直至达到100%。

3.4 站点时间序列

将以上27个模拟中的重要站点数据输出，可以得到27条不同的风暴过程水位线（见图9a）。结合不同台风样本的发生概率，可以计算站点在指定发生概率范围内的增水范围。图9b呈现了芦潮港站预报的100%和70%概率范围的增水，其中垂向虚线指示预报时刻。图中可以看到随着预报时效的延伸，预报的增水100%范围先随着台风预报误差的增长逐渐加大；然后台风逐渐远离预报区域，即使台风预报误差更大，但由于台风对预报区域的水位影响逐渐减小，增水100%范围也逐渐减小。2019年8月11日00时，由于台风运动到离长江口较近的位置，台风预报误差也较为显著，此时的增水100%范围达到最大。芦潮港站最大可能增水几乎可以达到2 m，最小可能增水却只有-1.80 m，显示出显著的风暴潮预报的不确定性。预报增水70%范围同样也呈现先增大后减小的特征，但其最大范围比100%范围显著减小，其发生的时间也与100%范围的略有差异。

图7 2019年8月9日14时芦潮港站的台风增水累积概率图

图8 台风“利奇马”24 h预报时效结果（预报时间为2019年8月9日14时)

图9 台风“利奇马”影响下不同站点的预报结果（预报时间为图中虚线所示的2019年8月9日14时）

图9显示本次预报的最大增水值比实测的最大增水显著偏大，这是由于台风预报的最大风速显著偏大所致；特别是在24 h预报时效，即使是概率偏小的最大风速，依然大于实测的最大风速。

表2 4个潮位站概率预报结果统计检验（概率预报为平均值加1个标准差）

3.5 集合预报精度检验

采用图4中台风“利奇马”期间南门、堡镇、芦潮港和金山嘴4个潮位站的逐时观测数据对本文建立的概率预报系统进行检验。最大增水及其发生时刻是风暴潮预报的重要参数。本文对其误差进行统计，同时将分析路径作为确定性单一预报的代表，所得结果见表2。

台风“利奇马”活动期间，南门、堡镇、芦潮港和金山嘴4站的实测最大增水分别为0.84 m、0.66 m、0.84 m和1.28 m。9日14时对台风“利奇马”近中心最大风速的24 h预报风速显著偏大，对于24 h预报，即使是概率偏弱的风速（30.4 m/s）也比实测风速（28 m/s）更大（见表1），因此，预报实验所得最大增水的估计也较实测最大增水显著偏大。27个集合样本所得到的4站最大增水分别为0.80～3.33 m、0.73～2.88 m、1.24～2.47 m和1.73～3.11 m。集合平均增水过程的最大值（见图9红线）分别为1.87 m、1.50 m、1.90 m和2.56 m，都较实测最大增水偏大约1.0 m。相对于确定性单一预报，本文建立的概率预报系统对最大增水的估计都更接近实测最大增水，表明本文建立的概率预报系统可以通过集合消减来提高风暴潮的预报精度。不仅如此，概率预报系统还能够提供更多的可用信息，如表2中所示的各站最大增水的标准误差。误差越大，表明台风路径和强度等参数中的预报误差对该站的最大增水影响越显著。对于这些海域，实际业务中需要特别注意台风实际路径变化的影响。

台风“利奇马”期间，南门和堡镇站的实测增水在8月10日17时达到最大，芦潮港和金山嘴站则在14时达到最大。概率预报和确定性预报对最大增水时刻（约24 h预报时效）的估计误差为2 h左右（见表2）。概率预报和确定性预报的最大增水时刻误差没有显著差异，但采用概率预报技术可以得到更多的信息。例如：虽然预报南门和堡镇站的最大增水时刻与实测只相差约1.5 h，但其标准差分别为5.14 h和6.38 h，比芦潮港和金山嘴的标准误差明显更大。因此，如果台风的移动路径、移动速度和台风强度等参数出现偏差，这两站的最大增水时刻有可能变化更大。

4 结论与展望

本文采用郭文云等[1]的台风集合构造方案，基于WRF气象模型和FVCOM三维海洋模型，建立了一套适用于长江口的风暴潮概率预报系统。为兼顾台风风速空间变化强烈和移动范围大的特点，WRF模型采用三重嵌套模式，FVCOM海洋模型采用两重嵌套，在保证精度的前提下提高计算效率。系统采用9条可能路径和3个可能最大风速进行27个风暴过程模拟。每个模拟风暴过程的发生概率都根据历史预报误差分布进行确定。

风暴潮概率预报较单一预报的优点是提供了多个可能的预报增水值及各种丰富的预报信息。本文建立了一套详细分析风暴潮概率预报结果的产品，如最大和最小可能增水值、平均增水值及其标准偏差、概率增水值及增水概率值等。基于以上分析，我们可以从风暴潮概率预报结果中得到更全面的信息，为防灾减灾工作提供更多数据参考。当然，本文所建立的长江口概率预报方案还存在一些可改进之处。其一，预报方案只考虑了预报台风的路径及最大风速误差，没有考虑预报台风中心气压的误差，而台风中心气压直接影响到台风的最大风速半径（见式（1）），是决定台风大小的重要参数，因此，台风中心气压的预报误差也应该在风暴潮概率预报系统中加以考虑；其二，本文所建立的概率预报系统采用了9条可能路径，使得预报系统运行过程中需要用到较多的计算资源，需要进一步做不同台风集合方案的对比，在计算资源与预报精度上取得合理的平衡；其三，本文所采用的概率预报方案是基于历史统计得来的，实际预报过程中，不同台风的预报误差也有所不同（如对于转向台风及走向异常的台风，其预报误差往往更大）。所以，加强与台风预报部门的合作，争取采用实时的台风预报误差数据是改进风暴潮概率预报系统的根本途径。