数值模式与统计模型相耦合的近岸海浪预报方法

毛科峰，萧中乐，王亮，季卫海

（1.解放军理工大学气象海洋学院，江苏南京 211101；2.解放军96631部队，北京 102208；3.解放军92899部队13分队，江苏宁波 315200）

数值模式与统计模型相耦合的近岸海浪预报方法

毛科峰1，萧中乐2，王亮1，季卫海3

（1.解放军理工大学气象海洋学院，江苏南京 211101；2.解放军96631部队，北京 102208；3.解放军92899部队13分队，江苏宁波 315200）

针对数值模式和统计模型预报近岸海浪存在的局限性，构建了数值模式和统计模型相耦合的近岸海浪预报框架，在模式计算格点和近岸预报目标点之间定义一个海浪能量密度谱传递系数，通过经验正交函数分解和卡尔曼滤波方法建立传递系数的统计预报模型并与数值模式进行耦合。经过对近岸波浪观测站1 a的预报试验表明：该方法能够提高近岸海浪有效波高预报精度，有效波高的均方根误差降低了约0.16 m，平均相对误差降低约9%。进一步试验和分析发现，该方法的预报有效时间小于24 h，将海浪能量密度谱经过分解后得到的基本模态反映了近岸波侯的主要特征，海浪能量密度谱传递系数的变化体现了波侯的季节变化特点。

近岸海浪；数值模式；统计模型；耦合；卡尔曼滤波

1 引言

波浪由外海传播至近岸过程中，除了包含波浪生成、耗散、波波非线性相互作用等物理过程外，由于受到水深、地形变化、岛屿、岛礁或近岸建筑物阻挡，会发生折射、绕射、反射、浅水变形、底摩擦、破碎等一系列复杂的变形。因此在近岸海浪预报模型中，准确描述这些物理过程是预报的关键。以WAM、SWAN、WAVEWATCH-Ⅲ为代表的第三代海浪能谱数值模式采用基于组成波谱能量平衡方程中能量源汇包含的物理过程有：风能输入、波浪破碎和白帽破碎时能量的耗散、波与波之间的非线性相互作用、底摩擦等，以及地形导致折射和变浅效应、流场导致折射、三波相互作用等过程，但对波浪折绕射、反射等物理过程的刻画还存在一定的局限［1］。数值模式中也还存在外强迫场输入误差，非优化的模式参数设置以及模式的数值求解方法等因素导致的误差等，这些因素制约着近岸海浪的预报准确度［2］。

构建统计模型也是进行近岸海浪预报的一条有效途径，众多学者［3—5］试图采用统计的方法去构建预报模型，如多元回归，神经网络、卡尔曼滤波等方法，尤其是卡尔曼滤波方法是一种具有严格数学理论基础的统计优化方法，以历史观测数据和当前观测数据递推当前时刻的状态，在气象和海洋业务预报中应用广泛，且效果较好。早期的研究主要是基于这些统计方法，选择海面风、海面气压、预报前若干时次的海浪实测值等因子去直接构建近岸海浪要素的预报模型。后来的研究基于这种思路，采用了新的统计方法或人工智能等算法，如遗传算法、支持向量机、决策树、遗传编程等构建预报模型［6—7］。ArıGüner［8］利用黑海西南部海岸气象观测站的风速、风向等作为预报因子构建了Karaburun海岸的海浪统计预报模型，并利用观测数据验证了预报效果。Breivik等［9—10］提出了近岸海浪预报的动力和统计降尺度方法，将数值模式的计算结果和实测波浪资料结合起来，采用了多种统计方法去建模，这主要是用统计手段作为数值模式结果的一种后处理。Reikard和Rogers［11］比较了以SWAN模式为代表的近岸海浪动力模式和统计模型的预报性能，指出了动力模式和统计模型的各自的优势和统计模型预报时效性不如动力模式长的特点。本文基于将数值模式和统计模型优势互补的思想，构建一种近岸海浪数值模式（SWAN）和统计方法相互耦合的预报框架，将数值模式的预报结果和近岸的观测资料通过优化的方法结合起来提高预报准确率，并利用实测海浪数据对这个预报框架模型进行试验和检验，分析模型对近岸海浪的预报性能。

2 海浪模式和统计模型耦合的预报思路

2.1 波浪能量密度谱传递系数的提出

本文研究的预报对象是近岸、港口某目标点的有效波高，如图1所示，O点为海浪预报的近岸目标点，当采用数值模式进行预报时，M点为模式计算网格点中距离O点最近的点，往往将该点的模式输出值作为O点的预报结果，但是两者之间常常有一定误差。主要是由于受到计算网格点分辨率的限制和数值模式本身的缺陷，图1中M点的模式结果和O点的海浪观测资料之间的误差反映了波浪模式对近岸过程中波浪发生变形等物理过程的刻画存在的不足。故从数值模式在M点较长时间段预报结果和O点的相应观测资料中挖掘出这两者之间的关系，可建立起耦合预报的技术途径。

海浪模式求解时，在空间网格点上求解海浪能量密度谱守恒方程，计算出格点上每个预报时次的波浪能量密度谱F（f，θ），然后根据波高、周期的统计分布关系式，计算出波浪要素，如有效波高、周期等，在数值模式SWAN中，有效波高的计算式为：

图1 模式计算格点和预报点Fig.1 Computational grids and forecasting points of the model

将近岸预报目标O点的有效波高记为Hso，数值模式计算得到M点离散的波浪能量密度谱值记为FiM，，由于海浪数值模式在近岸物理过程的局限和数

j值计算等方面的原因导致了直接利用式（2）从数值模式计算得到M点离散的波浪能量密度谱值计算O点的有效波高必然存在一定误差，故为了描述M点波浪能量密度谱值和O点有效波高之间的差别和关系，定义一个模式计算格点到近岸目标点的波浪能量密度谱传递系数矩阵，记为Βi，j其维数与FiM，j相同，则可借鉴式（2）得到一个通过修正M点的模式计算波浪能量密度谱值而得到近岸的预报目标O点有效波高的关系式：

2.2 预报模型建立

式中，βt是t时刻的前nk阶特征向量对应的波浪能量密度谱传递系数矩阵，它的维数为（nk×1）这一线性表达式，可以被视为卡尔曼滤波中的量测方程，εt为量测噪音，它是一维随机向量。

将波浪能量密度谱传递系数矩阵βt视为卡尔曼滤波系统中的状态向量，用方程式（10）描述其变化和量测方程一起构成卡尔曼滤波系统［13］：

式中，δt－1是动态噪音，δt－1和量测噪音εt都被认为是随机向量，假设两个量互不相关，且均值为0，方差为W和V，则应用广义最小二乘法，可以得到式（11），构成递推系统，其中为了简化表达式，记Xt＝AEtCt，由于H的下标s仅表示有效波高的含义，为了使公式简洁，将s略去，加入表示时刻的下标t。

2.3 预报流程

图2 耦合预报框架图Fig.2 Framework of coupling forecasting

基于上述的卡尔曼滤波系统和海浪数值模式可构成耦合预报框架，如图2所示，利用数值模式输出的波浪能量密度谱的时间序列进行经验正交分解，得到特征向量和特征值（本文也称为波浪能谱主模态和强度系数），然后利用卡尔曼滤波方法，根据近岸站点的观测有效波高值递推的波浪能谱主模态对应的从模式到近岸的能量密度谱传递系数，在图中被称为耦合预报的统计分析阶段；在实时预报阶段，对于某预报时次新得到的模式波浪谱，可利用统计分析阶段得到的波浪谱特征向量将波浪谱分解得到波浪能谱主模态，然后根据卡尔曼滤波系统得到的近岸能量密度谱传递系数进行波浪能量密度谱重构，进而进行近岸海浪有效波高预报。

在实时预报中，数值模式预报出计算格点上t时刻波浪能量密度谱，利用获取的近岸预报目标点的实测波高值，可以启动上述耦合预报流程，如图3所示，进行t时刻后的Δt时间段内的预报，Δt是耦合预报的时效，因此将上述的递推系统式（11）可写为：

式中，Δt可取为1，2，3…，取值的上限就是这个预报方法有效预报时效，具体情况下文将继续讨论它的取值，因此能量密度谱传递系数，既具有了空间尺度上模式格点和预报点之间能量传播变化的信息，还包含了Δt时间内的能量密度谱随时间演变的信息。

图3 预报流程示意图Fig.3 Schematic diagram of forecasting operation procedure

3 预报试验

3.1 试验方案

3.1.1 预报试验设计

为了分析预报模型的效果，选取琉球群岛海域中近岸站点NAHA站进行预报试验，从2008年1月1日1时起进行耦合预报，分别进行3 h、6 h、12 h和24 h预报，从1月1日1时至12月31日23时，每两小时进行一次四个时效的预报试验。模式预报是采用WAVEWATCH-Ⅲ模式和SWAN模式嵌套运行，其中WAVEWATCH-Ⅲ数值模式计算范围为18°～38°N，116°～138.75°E，格距为15′，并为SWAN模式提供嵌套波谱边界条件，SWAN模式计算范围为24°～30°N，126°～131°E，格距为3′；模式包括风浪相互作用、波波相互作用、白帽耗散、底摩擦作用等物理过程，风场采用CCMP资料［14］。图4中给出了海域内正方形标出的近岸海浪NAHA测站，该处波浪观测资料由日本NOWPHAS计划提供，为每两小时一次有效波高值。

3.1.2 预报模型参数的计算

图4 近岸数值模式计算区域Fig.4 Computational area of coastal numerical model

C0的取值：C0是的误差方差阵，认为是从样本资料直接计算得到，与理论值相等，故C0取为nk阶零方阵。

V的取值：由于量测量只有Hs一个量，故V为1 ×1矩阵，数值为V＝q／（k－m－1），q为第一组正交分解后线性回归的残差，预报与实际的差的平方和，k为取得样本个数，即数据的时间个数，m为选取主成分的个数。

3.2 试验结果分析

3.2.1 预报效果分析

首先对3 h预报的结果进行分析，试验中取m＝15 d的时间序列方向谱进行经验正交分解，取nk为9个主分量。图5为连续一年的NAHA站耦合预报结果和模式预报结果（距NAHA站最近的SWAN模式计算格点上的预报值）与观测资料的时间序列比较图，时间间隔为2 h。该图表明模式预报的有效波高值较实测值偏大，耦合预报结果与实测值非常接近，有效波高的偏度、均方根误差、平均相对误差均分别明显低于模式预报结果，如表1所示，其中耦合预报结果较模式预报的有效波高的均方根误差降低了约0.16 m，平均相对误差降低约9%。为了进一步分析耦合预报中两个参数设置对预报效果的影响，分别对正交经验分解的样本天数m和主分量个数nk的取值进行讨论，如表2和表3所示，给出了耦合预报均方根误差随选择样本天数和主分量个数的变化情况，通过试取值分析后发现，预报误差对两个参数的选择比较敏感，由表2可知，以15 d的模式计算波浪谱作为波谱分解的样本，以15 d为窗口构建递推系统使得预报误差较小，选择更多的天数后预报误差基本趋于稳定，故本文将m取为15 d。同理对主分量个数nk通过试取值分析后发现，如表3所示，经验正交分解后重构波谱时，收敛较快，12个主分量的方差贡献均在98%以上，预报误差较小，取更多的主分量，误差减小甚微，综合考虑计算效率，故nk取值为12，甚至发现nk大于16以后随着取值增加，误差也反而会增大，这也许和特征向量的方差贡献有关。

分析2008年1月2日14时之前14天168个时次分解的12个主模态与模式计算波谱分解得到对应的12个特征向量，如图6所示为2008年1月2日14时的方向谱，图7a、7b、6c、6d为谱分解后前4个主模态，方向谱的结构表明，这个时次主要是西北向波浪传播至该点，主模态的成分相对比较简单，也反映了这一主要特点，能量分布主要集中在150°至180°的方向。

表1 预报结果与浮标观测值统计分析结果Tab.1 Statistical analysis of forecasting results and buoy observations

表2 预报结果与浮标的观测值误差随m取值的变化Tab.2 The root-mean-square error between the forecasting results and the buoy observations with the change of the m values

表3 预报结果与浮标的观测值误差随取值的变化Tab.3 The root-mean-square error between the forecasting results and the buoy observations with the change of the nkvalues

3.2.2 预报的时效

对式（12）中Δt分别取3、6、12、24 h进行耦合预报，从2008年1月1日1时至12月31日23时进行预报试验，如图8为6 h预报结果的时间序列对比，（其他时效的结果图省略），不同预报时效的预报误差如表4所示，随着预报时效延长，预报误差逐渐增大，24 h的预报均方根误差为0.363 m较数值模式本身的预报误差0.324 m更大，这说明，本方法的预报时效在24 h以上不具备预报能力。

表4 预报误差随时效的变化Tab.4 The forecast error with the change of forecasting period

3.3 预报方法的机理讨论

图5 NAHA站1年的有效波高预报与实测值的比较Fig.5 Comparison of forecasted and observed significant wave heights of NAHA station during one year

图6 2008年1月2日14时的方向谱Fig.6 Directional spectrum at 14 o'clock on January 2nd，2008

图7 2008年1月2日14时的方向谱分解的主模态Fig.7 Principal component of directional spectrum decomposition at 14 o'clock on January 2nd，2008

图8 6 h预报结果的比较Fig.8 Comparison of 6-hour forecasting results with the observations

图9 波谱主模态Fig.9 Principal components of the spectrum

图10 能量密度谱传递系数Fig.10 Transfer coefficients of the wave energy density spectrum

由上述分析可见，本文的耦合预报方法效果较好，该方法的核心是将模式计算的波浪谱进行分解，定义了一个模式计算格点到近岸目标点的波浪能量密度谱传递系数，通过对系数的预报，描述模式的波浪能量密度谱值和预报点的有效波高之间的误差关系。为了探讨这种预报方法的机理，对数值试验的1年的波谱主模态和波浪能量密度谱传递系数的特点进行分析，如图9为12个主模态，第一、二模态体现了NAHA站点最主要的波浪特征：由于该站点位于琉球群岛西南部，常年受东亚季风影响，且东南面受岛屿遮挡，地形效应很明显。在西北风作用下，波浪有外海向岛屿近岸传播，以东南向为主，故第一、二模态的能量集中在150°～200°方向内，且以低频能量为主，从第一、二模态对应的能量密度谱传递系数的变化，如图10所示可以证明这个特点，在冬半年月份，这两个模态的系数值较大，在夏半年这两个模态的强度就非常的弱。夏半年盛行南风时，波浪由岛屿向外海传播，从第三模态可见波浪能量集中在320°～360°和0°～45°方向内，由于东南面的岛屿遮挡，这个模态的强度较之第一、二模态弱，且从能量密度谱传递系数的变化可见，主要集中在8－10月。其他模态和能量密度谱传递系数也体现了季节转换期的波浪特点。总之，海浪能量密度谱经过经验正交分解得到的主模态反映了近岸的波侯主要特征，海浪能量密度谱传递系数的变化体现了波侯的季节变化，因此本文的预报方法的机理就是通过描述近岸波浪能量密度谱的主模态的变化来预报近岸的有效波高值。

4 小结

（1）本文构建了近岸波浪数值模式和统计模型相耦合的预报框架，在计算格点和近岸预报点之间定义一个海浪能量密度谱传递系数矩阵，通过经验正交分解和卡尔曼滤波方法建立了传递系数的统计预报模型并与数值模式进行耦合。

（2）经过NAHA站1 a的预报试验表明，该方法能够较好地提高海浪有效波高预报精度，有效波高的均方根误差降低了约0.16 m，平均相对误差降低约9%。但是试验表明该方法的预报有效时间不大于24 h。

（3）通过试验说明了正交经验分解的样本天数m和主分量个数nk的取值对预报效果存在一定的影响。

（4）在耦合模型中，海浪能量密度谱经过经验正交分解得到的主模态能够反映近岸的波侯主要特征，海浪能量密度谱传递系数的变化体现了波侯的季节变化。

［1］Tolman H L，Banner M L，Kaihatu J M.The NOPP operational wave model improvement project［J］.Ocean Modelling，2013，70：2－10.

［2］Dietrich J C，Zijlema M，Allier P E，et al.Limiters for Spectral Propagation Velocities in SWAN［J］.Ocean Modelling，2012，70：85.

［3］Özger M.Significant wave height forecasting using wavelet fuzzy logic approach［J］.Ocean Engineering，2010，37（16）：1443－1451.

［4］Tseng C M，Jan C D，Wang JS，et al.Application of artificial neural networks in typhoon surge forecasting［J］.Ocean Engineering，2007，34（11／12）：1757－1768.

［5］Zamani A，Solomatine D，Azimian A，et al.Learning from data for wind-wave forecasting［J］.Ocean Engineering，2008，35（10）：953－962.

［6］Reikard G.Forecasting ocean wave energy：Tests of time-series models［J］.Ocean Engineering，2009，36（5）：348－356.

［7］Roulston M S，Ellepola J，Hardenberg J，et al.Forecasting wave height probabilities with numerical weather prediction models［J］.Ocean Engineering，2005，32（14／15）：1841－1863.

［8］ArıGüner H A，Yüksel Y，Özkan Cevik E.Estimation of wave parameters based on nearshore wind-wave correlations［J］.Ocean Engineering，2013，63：52－62.

［9］Breivik R，Gusdal Y，Furevik B R，et al.Nearshore wave forecasting and hindcasting by dynamical and statistical downscaling［J］.Journal of Marine Systems，2009，78：S235－S243.

［10］Camus P，Mendez F J，Medina R，et al.High resolution downscaled ocean waves（DOW）reanalysis in coastal areas［J］.Coastal Engineering，2013，72：56－68

［11］Reikard G，Rogers W E.Forecasting ocean waves：Comparing a physics-based model with statistical models［J］.Coastal Engineering，2011，58（5）：409－416.

Coastal wave forecasting by dynamical model coupled with a statistical method

Mao Kefeng1，Xiao Zhongle2，Wang Liang1，Ji Weihai3

（1.Institute of Meteorology and Oceanography，People's Liberation Army of China，University of Science and Technology，Nanjing 211101，China；2.The 96631 Unit of People's Liberation Army of China，Beijing 102208，China；3 The No.13 Branch Unit of the 92899 Unit of People's Liberation Army of China，Ningbo 315122，China）

The coupled coastal wave prediction scheme，which is a combination of a multi-scale numerical model and a statistical method，is proposed in order to avoid the limitations of one single scheme.By ocean wave model，the wave energy density spectrum of the computational grid in the coastal model is forecasted.We have defined a transfer coefficient matrix for thewave energy density spectrum between the computational grid and the coastal forecasting point.A statistical model for the prediction of this transfer coefficient is established using empirical orthogonal function（EOF）and Kalman filtering method.This statistical model is then coupled with the numerical model.The wave energy density spectrum of computational grid is optimized using the observed coastal buoy data.The coastal wave forecasting are validated by the observations of NAHA station for one year，indicating that this coupled method significantly improved the prediction power compared with the numerical model on its own.The rootmeansquare error of the significant wave height reduces about 0.16mand the average relative error is reduced by about 9%.It is also found that the forecasting accuracy of this method is limited within 24 hours；the principal components decomposed from the wave energy density spectrum reflect the main characteristics of local wave climate；and the change transfer coefficient of the spectrum reflects the seasonal variation of the wave climate.

coastal wave；numerical model；statistical method；coupled scheme；Kalman filtering method

P731.33

A

0253-4193（2014）09-0018-12

毛科峰，萧中乐，王亮，等.数值模式与统计模型相耦合的近岸海浪预报方法［J］.海洋学报，2014，36（9）：18—29，

10.3969／j.issn.0253-4193.2014.09.003

Mao Kefeng，Xiao Zhongle，Wang Liang，et al.Coastal wave forecasting by dynamical model coupled with a statistical method［J］.Acta Oceanologica Sinica（in Chinese），2014，36（9）：18—29，doi：10.3969／j.issn.0253-4193.2014.09.003

2013-04-01；

2014-01-26。

国家自然科学基金项目（41331174，11102232）。

毛科峰（1981—），男，湖南省常德市人，主要从事物理海洋和海洋环境调查技术研究。E-mail：maomaopla＠163.com