上海邻近海域风暴潮数据同化与特征分析

丁骏，吕忻，姚雅倩，孟鑫，姜雪敏，吴旭云，葛建忠

(1.上海市海洋监测预报中心，上海 200062；2.国家海洋局东海预报中心，上海 200136；3.华东师范大学 河口海岸学国家重点实验室，上海 200062)

1 引言

风暴潮是世界上最严重的海洋灾害之一[1-2]，也是上海城市安全的主要威胁之一，风暴增减水[3]是一种非周期性的对众多因素敏感的又备受关注的复杂现象。目前，国内外广泛采用海洋站和卫星等监测到的海面高度与正常天文潮的代数差的方法来计算风暴增减水，即风暴潮大小[4-6]。由于风暴潮与天文潮的非线性耦合效应存在[7]，实际增减水并非两者简单的线性相减。受海洋站空间分辨率限制，又常常以点代线、以点代面评估近岸海域风暴潮大小，其误差不断累积。虽然卫星高度计运用到风暴潮观测的研究[8]已经取得技术突破，但其时空分辨率和精度仍需提升。随着理论技术发展与实测数据的丰富，风暴潮数值预报[7]的准确率和时效性进一步提高，可以直观反映风暴增减水、漫滩过程，然而由于数值模式反映的物理过程是有限的，初始场和强迫场是近似的，模式参数是经验的，预测结果不可避免存在误差。在实践中，用实测潮位分离法计算的增水值通常偏大，该值可达15～30 cm，而数值模拟的增水值往往偏小[9]。

集合卡尔曼滤波（Ensemble Kalman Filter，EnKF）[10-11]可以在受噪声影响的测量值与数值模拟的预测值之间寻找最优状态估计。它基于集合预报思想，源于卡尔曼滤波理论和蒙特卡洛估计方法[12]，克服了卡尔曼滤波技术计算量大、只适用线性模式的缺点。由Evensen[10]在1994 年提出并应用到海洋同化领域后，日渐成为国际上资料同化研究的热点。在风暴潮资料同化应用方面，Butler 等[13]和Altaf 等[14]提出了一种基于EnKF 和ADCIRC（Advanced Circulation）风暴潮模型耦合的风暴潮预报方法，Colle 等[15]运用En-KF 对11 个风场资料同化驱动ADCIRC 模型来模拟“桑迪”飓风风暴潮，Etala 等[16]展现了局地集合变换卡尔曼滤波（Local Ensemble Transform Kalman Filter，LETKF）在卫星遥感与验潮站数据同化方面的优势，这些研究主要集中在EnKF 与风暴潮模型、气象模型的耦合同化，EnKF 与卫星、验潮站等数据融合，其方法的改进和拓展等3 个方面。在风暴增减水特征分析方面，近30 年来风暴潮理论[17]没有重大突破，也基本能够满足对风暴潮物理机制的解释，但在不同台风影响下，长江口区域风暴增减水变化趋势不一[18]，人们关注的近岸海域风暴潮的时空分布及规律仍以验潮站、台风等数据统计分析和风暴潮数值模拟为主[18-19]，通过风暴增减水数据同化反演并进行分析鲜有报道。

本文选择201810 号台风“安比”登陆上海过程，采用EnKF 方法，对不同风暴潮资料进行同化融合，并进行结果验证和时空分布特征分析。同化中增加了Schur 乘积算法[20-21]，以消除运算中产生的滤波发散和伪相关问题。以期为风暴潮机理分析、风暴增减水再分析、风暴潮监测预报等研究提供技术支撑，为海洋防灾减灾工作提供服务保障。

2 数据和方法

2.1 研究区域与数据来源

研究区域为受风暴潮影响明显的长江口内、口外和杭州湾北海域。潮位监测站包括9 个验潮站（崇明南门、堡镇、高桥、吴淞、芦潮港、金山嘴、大戢山、滩浒岛和佘山），其中，高桥、大戢山和滩浒岛数据用于验证，不参与同化，其余6 个站为参与同化的站，见图1。

为计算风暴增减水值，搜集了2016−2018 年9 个海洋站的水位资料，资料符合《海洋观测规范（GB/T14914—2018）》规定。

模型计算所需的风、气压强迫场资料来自上海市气象台订正的欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecast，ECMWF）预报产品，预报时长为72 h，时间分辨率为3 h，空间分辨率为0.125°×0.125°，覆盖区域0°～55°N，70°～140°E。径流数据来自长江水文网（www.cjh.com.cn）大通站流量数据。岸线地形资料来源为美国国家地理数据中心（National Geophysical Data Center，NGDC）的分辨率为1′的ETOPO1 全球地形数据和精度为1∶250 000的岸线数据、航海保障部最新海图水深和岸线资料、长江口航道管理局实测地形资料。开边界调和常数基于全球潮汐模型TPXO8[22]和中国海洋图集资料[23]调整得到。

图1 数据同化区域、监测站点分布和部分台风路径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）来自中央气象台台风网（http://typhoon.nmc.cn/web.html）。

2.2 研究方法

2.2.1 天文潮调和分析

风暴增减水由观测水位（总水位）减去天文潮水位得到，其中，天文潮采用调和分析法[26]计算，公式表示为

式中，Zt表示天文潮位；A0表示平均水位；fi表示分潮i的交点因子；Hi表示分潮i振幅；qi表示分潮i角速度；(V0+u)i表示分潮i初相角；gi表示分潮i迟角；fi和(V0+u)i与时间t有关；n表示分潮数目。采用最小二乘法[26]分别对9 个站3 年逐时的实测潮位逼近得到68 个分潮调和常数Hi和gi，进而实现台风“安比”过境期间天文潮位计算。

2.2.2 风暴潮数值模拟

FVCOM（Finite-Volume Coastal Ocean Model）[27]采用有限体积数值离散方法求解三维水动力原始控制方程组，运用不规则三角形网格和干湿网格处理技术[28]，可以精细模拟复杂边界和地形，在近岸海域的暴风潮数值预报中有明显优势[29]。

本文利用FVCOM 对长江口风暴潮进行模拟，计算区域为28°21′～34°17′N，117°33′～124°30′E，水平网格分辨率最高达到100 m，垂向 σ坐标均分为10 层（图2）。外模计算的时间步长为10 s，内外模时间步长比为10∶1，开边界潮位设置13 个分潮调和常数（M2、S2、K1、O1、N2、K2、P1、Q1、Mf、Mn、M4、MS4、MN4）。风暴增减水场等于加风场模拟的总水位减去不加风场模拟的天文潮位，结果中包含了天文潮和风暴潮的非线性耦合效应。

2.2.3 集合卡尔曼滤波同化方法

图2 长江口海域计算网格和水深Fig.2 Computing grids and water depth in the Changjiang River Estuary

采用双线性插值方法将数值模拟值和观测值插值到研究区域网格（0.008°×0.008°）进行空间分辨率匹配；利用蒙特卡洛随机扰动法（Monte Carlo Forecast，MCF）[30-31]生成状态变量和观测值随机变量集合。状态变量为FVCOM 计算的逐时风暴增减水场，观测值为6 个海洋站的逐时增减水值。状态变量预测和观测集合公式为式中，Xt,i、Xt−1,i为t时刻、t−1 时刻第i个样本的状态变量；F为状态转移矩阵，表示前后状态值之间的关系；Wt−1为t−1 时刻的系统噪声；Yt,i为t时刻第i样本的观测值；H为观测矩阵，表示测量值与状态值之间的关系，由于测量了状态中的一个变量，且同一时刻风暴潮测量值与“真实”状态值之间接近，故取值为1；Vt,i为t时刻的观测噪声；系统噪声和观测噪声的方差为Q和R。

在状态变量的误差协方差矩阵计算中，采用Schur乘积算法[20-21]对其进行局地化[32]以提高估计质量，逐步滤除距观测点较远且不敏感的网格上的伪相关，减轻求逆矩阵不满秩的问题。误差协方差矩阵和Schur 乘积算法公式为

式中，Xt为状态变量预测集合的平均值；n为样本数，取值为200，由下文讨论给出；Et为状态变量预测集合中的扰动值；Pt为状态变量预测协方差矩阵；T 为矩阵转置运算；Ω为局地相关函数，是基于距离的相关系数矩阵；b为状态值到观测点的距离；a为Schur 半径，表示与观测点数据相关长度或区域的大小，a取0.8°，由下文讨论给出。

根据观测值和预测值，结合预测协方差Pt，计算卡尔曼滤波增益Kt，得到更新后的状态变量Xt,i和Pt，并以此循环上述步骤。卡尔曼滤波增益、状态变量和误差协方差矩阵更新公式为

式中，Kt为卡尔曼增益矩阵；I为单位矩阵；为更新后的状态变量；为更新后的误差协方差矩阵。

3 数据同化检验

3.1 实验设计

本文选择1949 年以来第3 个直接登陆上海的台风“安比”风暴潮过程，将2018 年7 月21 日00 时至23 日23 时海洋站和FVCOM 数值模拟的风暴潮数据进行同化。为检验同化效果，设计了2 组实验，见表1。

方案A 为控制实验，不同化任何观测资料，在FVCOM 处于稳定预报状态时，从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.2 同化效果检验

图3 给出了6 个站实测水位计算的风暴增减水序列（实测）、6 个站数值模拟的风暴增减水序列（模型）和6 个站同化的风暴增减水序列（同化）对比结果。模型与实测对比可以发现，模型已能表现出风暴增水过程的大体特征，长江口内的崇明南门站、堡镇站和吴淞站增水转减水的时刻基本相同，但模型并未表现出实际增水的“双峰”型特点和减水后“回弹”趋势且增水高度低于实测；长江口门处的佘山站，模型最大减水与实际较一致，但最大增水误差较大；处于杭州湾北岸海域的芦潮港站和金山嘴站，模型和实测的最大增水较一致，对减水的刻画不够准确。尽管模型与实测给出的风暴增水存在差异，但经过同化后，其结果和实测符合的很好，具有较高的一致性。为了进一步分析同化结果的精度，本文统计了它们之间的平均绝对误差（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 可见，模型与实测的风暴增减水相比仍存在一定差距，长江口内高桥站的特征对比基本同崇明南门站、堡镇站和吴淞站，长江口门处的大戢山站和杭州湾海域的滩浒岛站特征对比分别同佘山站、芦潮港站，两者MAE 和RMSE 分别为0.16 m、0.19 m。同化与实测相比，同化得到的结果较为准确，二者仅存在细微的偏差，其MAE 为0.07 m，RMSE 为0.09 m，说明了同化结果有效性和同化方法的正确性。

图5 为同化实验中集合离散度和集合均方根误差的对比。两者是否接近是检验集合预报系统可靠程度的重要指标[34]。集合离散度是同化实验中集合成员与集合平均值之间的标准差（Standard Deviation，SD），集合均方根误差则是与观测值的差异[35]。除了集合离散度与均方根误差保持一致外，集合离散度还应较大，以保证真实状态包含在集合中[36]，过低的离散度会使集合过度自信，易导致滤波发散[37]。本同化实验中，集合离散度在0.04～0.22 m 之间，平均值为0.10 m，均方根误差在0.05～0.25 m 之间，平均值为0.11 m，集合离散度与均方根误差比值为0.90。在图5中，集合离散度整体上略低于均方根误差，可能由于数值模拟的结果偏低所致，但二者已展现出较好的一致性，进一步验证了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 时）数值模拟的风暴增减水场、同化实验的增减水场。沿着台风路径前进方向，数值模拟和同化实验的增水场均呈现增水值逐渐升高趋势，数值模拟的增减水范围为−0.15～0.12 m，同化实验的增减水范围为−0.07～0.53 m；通过典型时刻的9 个站增减水对比发现，实测水位计算的增减水与数值模拟的增减水之间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，两者较一致。

4 风暴潮特征分析

地形是控制和制约增减水分布的重要条件，风和气压变化是造成长江口增减水的主要动力因子[38]。本文按照长江口特定的地形条件，划分长江口断面、杭州湾北断面，分析同化后的台风“安比”风暴潮特征；结合台风“安比”登陆前风和气压状况，进一步分析口内“双峰”增水、台风眼区域增水和增水锋面特征。

图8 和图9 分别为台风“安比”登陆前长江口断面、杭州湾北断面各站同化后的增水值。长江口断面增水值在0.33～0.40 m 之间，东西方向的增水梯度不是很明显，崇明南门站增水最大，高桥站增水最小。可能由于地转流的影响，长江口南支北岸的增水略高于南岸，崇明南门站、堡镇站增水高于吴淞站和高桥站。佘山站可能受台风眼影响，增水略有升高。杭州湾北断面增水值在0.04～0.13 m 之间，由于处于台风第三象限，增水从湾内到湾外呈递增趋势，且梯度较明显，金山嘴站增水最小为0.04 m，芦潮港和大戢山站增水最大，为0.13 m。由于东南向风将水体吹送口内，西北风将水体输运湾外，长江口断面增水明显高于杭州湾北断面。

图6 台风登陆前数值模拟的增减水场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 中长江口内崇明南门、堡镇、吴淞和高桥站增水曲线呈现“双峰”特征，即在最大增水前后较短时间（约12 h）内出现另一个较大的增水峰，两峰高度基本一致。这与张文舟等[39]发现的福建沿海部分站点双增水峰现象较相同，他认为在台风横穿台湾海峡正面登陆福建时，台湾海峡特殊地形和台风大风区中心区域的经过共同造成。根据“水量平衡原理”[40]，在风暴潮“初振”阶段，向口内输入的水量不断增加，由于风的顶托、降雨等影响，总向海输出的增水量较小，因此形成了风暴增水的“前峰”。台风登陆前，强烈的气旋风应力使得水量继续向上游输送，导致口内增水突然升高，由此产生了增水的“后峰”。

图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

在台风登陆前，各站增水曲线均出现大小不一的波峰，如长江口内4 个站的“后峰”，芦潮港、金山嘴、滩浒等站在7 月22 日12 时的“小峰”。可见，台风登陆会对整个上海海域造成不同程度的增水。台风登陆前，佘山站处于台风眼附近位置，距离台风中心约21 km（图6）。按静力学计算，气压降低1 hPa，水位约升高1 cm[38]。台风中心气压值为980 hPa，以平均气压1 013 hPa 计算，约有33 cm 增水，这与佘山站增水36 cm 接近。图7 中断面a 处于台风眼附近位置，其增水梯度相对较大，进一步表明台风眼对风暴增水的影响明显。

图7 中区域b 是杭州湾梯度增水和长江口梯度增水的汇合锋面，由于杭州湾梯度增水由西南向东北堆积，长江口梯度增水自东南向西北，此处存在增水“切变”，风暴增水呈现涡的特征，对于风暴增水研究有一定意义。囿于篇幅，其他时刻风暴潮特征未尽分析。

5 讨论

集合卡尔曼滤波同化质量与使用的集合大小有关[41]，集合数目越多越可真实地描述系统状态的空间分布，但集合数目过多会增加系统运行的时间。目前，还没有很好的方法来确定最合适的集合样本数[42]，在大多数EnKF 的应用中，典型的集合大小数量是100 左右[41]。本文通过固定4 组不同的Schur 半径a，分别选择集合样本数n为50、100、200、400、500，在台风登陆前进行EnKF 同化，比较实测与同化之间的MAE 和RMSE 以及集合离散度和集合均方根误差，讨论集合样本数和Schur 半径对同化的影响，并给出适合的集合样本数和Schur 半径。

从图10 可以看出，集合样本数量增加，对同化结果有很显著的积极影响，当集合数量继续增大时，误差基本维持，同化结果却没有显著提升。如在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 半径，较小的集合样本数有更好的同化结果[43]。另外，随着集合样本数增加，集合离散度与集合均方根误差越趋于一致，同化系统可靠性加强。在a=0.6°和a=1.0°时，集合样本数从50 增加到500 时，集合离散度与集合均方根误差的比值分别从0.77 增加到0.95，从0.77 增加到0.94；在a=1.4°时，集合离散度与集合均方根误差的比值从0.74 增加到0.89；在a=0.2°时，两者一致性更高。

统计同化的单个步长时间，n=200 时为24 min，n=400 时为99 min，集合样本数增加1 倍，计算成本增加了3 倍多。除了计算效率外，集合样本数的选择还需考虑到实测和模型计算的步长等。由于本文的计算步长为60 min，再综合上述MAE 和RMSE 分析以及集合离散度与集合均方根误差的比值，故风暴潮同化的集合数在100～200。

随着Schur 半径的增大，即相关长度的增大，受观测数据影响的区域范围越来越大，同化过程就会保持越来越多的信息量，进而会产生越来越大的伪相关性。如在集合样本数为200 时，对比a为0.6°、1.0°和1.4°发 现，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°时，两条曲线越来越分开，在集合样本数为200 时，随着Schur 半径的增加，比值从0.89 减为0.85，再到0.81，系统稳定性有所降低。减小Schur 半径，可以提高同化的效果，但Schur 半径过小会导致部分信息丢失，使得风暴增水场出现一定误差，如a=0.2°时，集合样本数达到了100 时，误差已超出预期，并不能提高估计的质量。综合上述分析，当允许误差在0.10 m 以下时，风暴增水同化的Schur 半径可在0.6°～1.0°。实际上，在协方差局地化公式（7）中，Schur 半径的选择还可能受海洋要素类型、观测点位置、网格大小等多种因素的制约，Schur 半径与这些因素的关系还需要进一步研究。

图10 不同集合样本数及Schur 半径同化误差Fig.10 The assimilation results error of different set samples and Schur radius

6 结论

本文针对实测计算和数值模拟的风暴潮寻优问题，以典型台风“安比”登陆上海的风暴增减水为研究对象，通过协方差局地化的集合卡尔曼滤波方法进行风暴潮数据同化融合与特征分析。得到结论如下：

（1）同化获得了逐72 h 的风暴增减水的最优解，同化与实测计算的增减水对比，MAE 为0.07 m，RMSE 为0.09 m，集合离散度与均方根误差比值为0.90，同化效果较好且可信，为风暴潮、数值模拟结果订正等研究提供了新思路；（2）同化后的风暴增减水能够较好地刻画台风登陆前双峰增水、台风眼增水、增水锋面等特征，对风暴潮研究有一定意义；（3）协方差局地化技术可以提高风暴潮估计质量，降低滤波发散和伪相关性，减轻求逆矩阵不满秩的问题；（4）集合样本数和Schur 半径设置对同化效果影响较大，降低相关长度可以提高同化的效果，但Schur 半径过小会导致部分信息丢失准确度降低。