海浪资料同化是将数值模式与实测资料相结合, 用于给出真实状态(又称为分析值)的最佳估计的一种有效方法, 对改善海浪模拟和预报效果具有至关重要的作用(Lionello et al, 1992)。相对于大气资料同化, 海浪数据同化的研究起步较晚, 所采用的同化方法主要为变分方法(de las Heras et al, 1992, 1994; Bauer et al, 1996)、最优插值(Lionello et al, 1992; Young et al, 1996; Hasselmann et al, 1997; Voorrips et al, 1997; Greenslade, 2001; Waters et al, 2013)、集合最优插值(Qi et al, 2016)和基于静态集合样本的滤波方法(孙盟等, 2014)等。然而, 在海洋环流模式和大气模式中广泛开展的集合Kalman滤波(Kalman, 1960; Kalman et al, 1961; Evensen, 1994; Burgers et al, 1998)同化方法在海浪数据同化中的应用相对较少。如何构造稳定可信的海浪模式集合运行成为制约海浪模式集合滤波同化应用的一个关键问题。

已知, 对于空间单点的某个变量进行多次采样, 当采样次数足够多时, 该变量集合样本通常满足一定的统计分布(高斯分布)。针对同一变量, 对上述单点(又称中心点)附近的空间点进行多次采样, 同样的, 当采样次数足够多时, 邻点的变量集合样本也满足一定的统计分布, 且邻点与中心点的变量集合样本存在一定的相关性, 这也是海浪同化中背景误差相关结构的信息来源(Kalnay, 2003)。根据上述分析, 全球尺度范围内的空间点均满足上述两条特征。由此, 为保证海浪模式集合滤波同化应用的合理性和正确性, 所构造的海浪模式集合需要满足上述两条基本特征, 这也是衡量海浪模式集合构造是否合理的重要指标。

由于海浪模式对初始场的敏感性较弱(Lionello, 1992), 通过初始场扰动构造模式集合难以保证模式集合运行的稳定性。已知, 海浪模式由风场驱动, 若对风场进行持续有效的集合扰动, 可满足海浪模式集合运行的稳定性要求(Greenslade, 2004)。如何合理地对风输入进行扰动, 根据集合滤波方法的基本思想, 引入风场扰动量, 该扰动量的量级应与风场模拟误差的量级相当。本文提出三种风场集合扰动方案, 分别为纯随机数、随机场和时间滞后的风场扰动方法。纯随机数方案是对当前时刻每个空间位置的风场矢量乘以满足某分布的随机数, 该集合扰动量必然满足一定的统计分布。随机场(Evensen, 1994)扰动方案是预先按照蒙特卡罗方法生成三维伪随机集合样本, 该集合样本总体均满足正态分布, 且单个样本在空间上具有局地性和平滑性, 然后将该集合样本作为扰动量, 叠加给当前风场, 构成风场集合。时间滞后方案是指利用6h间隔的风场偏差近似风场模拟误差, 由于当前时刻仅存在一个风场偏差样本, 因此对单个样本乘以满足某分布的随机数生成扰动集合。

基于上述分析, 本文利用2014年1月ECMWF全球风场, 基于MASNUM-WAM海浪模式, 采用纯随机数、随机场和时间滞后方案, 构造风场集合, 由风场集合驱动生成海浪模式集合, 开展海浪数值模拟实验, 对海浪要素(有效波高)和二维波数谱集合样本进行统计分析, 对比三种风场集合扰动方案对于海浪模式的影响。上述为引言, 第一节介绍本文所采用的海浪模式和数据, 详细阐述了三种风场集合扰动方案, 第二节为集合海浪模式实验及结果分析, 最后一节为结论与展望。

1 风场集合扰动方案 1.1 模式与数据

本文采用球坐标系下的第三代海浪模式MASNUM-WAM (Marine Science and Numerical Modeling) (Yuan et al, 1992a, b; 杨永增等, 2005)。该模式应用了基于破碎波统计理论发展的海浪破碎耗散源函数(Yuan et al, 1986), 并采用复杂特征线嵌入计算格式。波数谱被离散成24个方向和25个波数, 对应频率范围是0.042-0.413Hz。本文模式计算区域为: 79°S-65°N; 0°-360°E, 空间分辨率为1.0°×1.0°, 时间步长为15min, 模式每1h输出一次。地形数据为ETOPO5。

本文模式驱动风场数据采用2014年1月ECMWF风场(ERA-Interim), 该数据由欧洲中期天气预报中心提供。风场数据时间间隔为6h, 空间分辨率为1.0°×1.0°, 覆盖范围为90°S-90°N; 0°-360°E。此外, 本文中采用的三种风场集合扰动方案是针对全球海域提出的, 具有普遍适用性, 不局限于某海域或某时段, 研究中采用的全球风场数据具有足够时长即可。

1.2 风场扰动方案

为保证海浪集合模式运行的稳定性, 本文对驱动风场进行持续集合扰动, 即对当前时刻风场叠加扰动集合。首先, 需要确定风场的扰动幅度。若扰动幅度过大, 由于风场变化过大, 则容易出现海浪模式结果跳跃不连续、集合发散等问题; 若扰动幅度过小, 则容易出现集合收敛的情况, 与构造海浪集合模式的初衷相违背。由此, 风场扰动幅度应与风场模拟误差的量级相当。

确定风场模拟误差。已知在基于繁殖向量法的集合预报研究中, 数值模拟的误差可以分为三种分量:增长分量、恒定分量和衰减分量。其中, 增长分量将随着模式的积分不断成长, 恒定分量不随模式的积分改变, 衰减分量则随着模式的积分不断减小。与NMC方法(Parrish et al, 1992)类似, 采用固定时间间隔的模拟结果的差异可代替模式误差(孙盟等, 2014)。具体做法如式(1)所示。

其中, Range表示模式误差的标准差, , , K表示采样个数, W表示风速。由于天气时间变化尺度较小, 此处固定时间间隔Δt=6h。

根据上述分析, 本文利用2014年1月ECMWF全球风场数据, 统计6h间隔风速偏差的月平均值及标准差, 如图 1所示。从图 1可以看出, 6h间隔风速偏差月平均的浮动范围介于±0.2m/s之间, 全球平均约为0m/s; 6h间隔风速偏差的标准差取全球平均为1.6212m/s, 西风带区域标准差较大, 介于3-4m/s, 其他区域标准差普遍低于2m/s。

构造风场集合的基本思路是将扰动集合叠加到当前时刻的风场, 生成该时刻的风场集合。上文已确定风场扰动幅度(=1.6212m/s), 下面具体介绍三种风场集合扰动方案。

(1)方案A:纯随机数扰动法

方案A纯随机数扰动法中, 扰动集合由(-1, 1)区间内满足平均分布的随机数与风矢量的乘积构成, 如式(2)所示。若扰动幅度过小, 模式集合易出现收敛的趋势, 违背了风场集合扰动的初衷; 若扰动幅度过大, 可能会导致海浪模式的溢出或崩溃。为了避免上述问题, 引入控制扰动幅度的参数α_A, 将α_A·β_{i, j, n}作为最终的风场扰动系数。

其中, W_{i, j}^b和W_{i, j, n}^a分别表示扰动前和扰动后的风场, i和j分别表示经度和纬度方向的模式空间网格点, n(n=1, 2, …, N)表示集合样本数。β_{i, j, n}表示(-1, 1)区间内满足平均分布的随机数。

为了确定扰动幅度参数α_A, 本文首先对2014年1月ECMWF全球风场数据进行月平均及标准差统计, 如图 2所示。从图 2a可以看出, 西风带附近海域风速可达10m/s以上, 全球风速月平均为6.2678m/s; 从图 2b可以看出, 西风带附近海域标准差偏大, 普遍大于3m/s, 全球风速标准差的均值为2.3841m/s。根据上文的分析, 已确定风场扰动幅度为1.6212m/s, 结合式(2), 故本文取。

(2)方案B:随机场扰动法

方案B随机场扰动法中, 采用Evensen (1994)的方法生成随机场, 作为扰动集合, 该三维伪随机集合样本总体服从正态分布, 单个集合样本空间上具有局地性和平滑性, 将该扰动集合叠加到风矢量场, 构成风场集合, 如式(3)所示。与方案A类似, 引入控制扰动幅度的参数。

其中, 表示采用Evensen (1994)的方法所生成的随机场元素。

Evensen (1994)采用蒙特卡罗(Monte Carlo)方法所构造的随机场空间分布较为平滑, 单点空间相关性随距离的增加而减小, 具有较好的局地性, 将该随机场作为扰动形成驱动风场集合不会增加额外的空间梯度或破坏海浪模式本身的稳定性。构造随机场需要预先给出该随机场的空间相关距离尺度, 随机场作为扰动叠加到风场, 故其相关距离尺度应与风场相关距离尺度(本文取5°)一致。

为了确定扰动幅度参数, 对每一个随机场样本计算其标准差。考虑到方案A中包含平均分布(已知β的标准差为), 为保持扰动幅度的量级一致, 根据上文的分析, 结合式(3), 故本文取。

(3)方案C:时间滞后扰动法

方案C时间滞后扰动法中, 采用固定时间间隔的风场偏差与(-1, 1)区间内满足平均分布的随机数γ的乘积构造扰动集合, 如式(4)所示。驱动风场为大气模式结果, 其本身已包含大气模式的模拟误差, 利用固定时间间隔Δt的风场偏差构造扰动样本, 该扰动样本保留了风场误差的演变及空间相关信息。由于两个时刻的风场偏差仅足以构造一个扰动样本, 由此, 本文对该扰动样本乘以满足平均分布的随机数γ_n, 从而生成扰动集合。与方案A类似, 考虑到风场扰动幅度对海浪模式的影响, 此处引入系数α_C, 用于调整扰动幅度。

其中, W_{i, j}^{b, t}和W_{i, j, n}^{a, t}分别表示t时刻扰动前后的风场。考虑到天气时间变化尺度较小, 此处取Δt=6h。γ_n表示(-1, 1)区间内满足平均分布的随机数(已知γ的标准差为)。根据上文的分析, 结合式(4), 故本文取α_C=1。

2 扰动实验与结果分析

利用2014年1月ECMWF风场数据, 按照1.2节中提出的扰动集合构造方案, 分别构造集合风场, 用于驱动海浪模式集合, 开展为期一个月的海浪集合数值模拟实验, 集合样本数为100。为了确保风场扰动的正确性, 本文首先对扰动集合样本进行抽样检查, 详见2.1节。驱动风场误差是海浪数值模拟误差的重要来源, 海浪模式本身对输入的驱动风场误差存在调整和自适应过程, 为考察不同的风场集合扰动方案对海浪模式的影响, 本文以海浪特征要素(有效波高)和二维波数谱为统计对象, 对比分析三种风场集合扰动方案对海浪模式的影响, 详见2.2和2.3节。

2.1 扰动集合样本分析

对1.2节中的三种扰动方案, 以2014年1月1日0时扰动集合样本为例, 提取第一个样本, 分析其空间分布特征, 如图 3所示。从图 3可以看出, 方案A的扰动样本分布呈现沙盘化, 该现象由随机数对风场空间结构的瓦解造成; 方案B的扰动样本呈现碎片状特征, 具有局地性和平滑性; 方案C的扰动样本保留了较多的风场空间结构信息(尺度较大), 这也是前两种方案所不具备的。

2.2 有效波高集合样本统计分析

(1)发散度时空分布

有效波高集合样本发散度可以反映模式集合发散的程度。为了分析有效波高集合样本发散度的时空分布特征, 采用式(5), 对减掉有效波高集合平均的三维变量S(样本发散度)进行统计分析。

其中, H表示有效波高(空间二维场), 上划线表示集合平均, n=1, 2, …, N(N=100)。发散度S为具有时空属性的三维变量, 月平均的空间分布如图 4所示; 对发散度S进行空间(全球)平均后, 可得到随时间演变的曲线, 如图 5所示。

从图 4可以看出, 三种扰动方案的发散度在中高纬度地区普遍偏高(0.4m以上), 低纬度区域基本低于0.15m;与其他两种方案相比, 方案B发散度空间整体分布更均匀。三种扰动方案的发散度在低纬度地区普遍偏小, 初步分析其原因为低纬度地区海浪有效波高随时间变化较小, 对集合扰动敏感性较弱。

从图 5可以看出, 三种扰动方案的发散度曲线随时间变化不大, 纯随机数和随机场扰动法的发散度介于0.1-0.2m, 时间滞后扰动法的发散度较前两者偏大, 在0.2-0.3m浮动。

已知, 风场集合扰动方案优劣的重要判据是海浪模式集合是否具有均匀特征、稳定性和持续性。综上所述, 以有效波高作为海浪模式集合发散度的检验指标, 随机场扰动法所生成的海浪模式集合的发散度最佳; 其发散度空间分布平滑且基本均匀, 其发散度的全球平均随时间变化不大, 基本保持在0.1-0.2m, 即风场扰动对有效波高的影响范围为0.1-0.2m。

(2)离散概率密度分布

海浪模式集合样本发散度反映的是集合的整体特征, 良好的时空发散度是证明集合构造方案合理的必要条件。已知, 集合滤波方法中的集合需要满足一定统计分布, 即海浪模式集合样本的统计分布特征越明显, 风场扰动方案越合理。

选取某地理空间位置, 若其有效波高集合样本满足一定的统计特征, 那么, 该集合样本的概率密度分布可以反映模式集合的离散程度。由此, 将有效波高的值域, 划分成较小的网格区间, 统计该地理空间位置上某时刻的有效波高集合在该网格区间内的样本数, 得到相应的离散概率密度分布, 由不同时刻的有效波高集合样本, 可以得到该点离散概率密度随时间的演变过程。由于集合样本扰动量相对于有效波高是个小量, 故采用式(6), 对减掉有效波高集合平均的集合样本H^*进行统计分析, 结果如图 6所示。

图 6为空间点(60°E, 30°S)、(180°W, 0°)、(30°W, 40°N)位置处集合样本H^*的离散概率密度分布随时间的变化情况。从图 6可以看出, 三种扰动方案的集合样本分布均接近正态分布, 但近似程度不一。由图 6b、e、h可以看出, 虽然空间点不同, 但随机场扰动方案的集合分布特征保持稳定, 而纯随机数和时间滞后扰动方案在不同空间点的集合分布变化较大, 时而发散时而收敛。即从随机误差统计分布的角度来讲, 随机场扰动为最佳方案。

(3)相关系数月平均分布

利用有效波高集合样本可以计算某地理空间点与周围网格点的相关系数, 该相关系数反映了海浪同化调整过程中所需的权重信息, 是衡量海浪模式集合构造是否合理的关键指标之一。选定某地理空间点作为参考点, 以有效波高为统计量, 计算其他模式网格点相对于参考点的相关系数。每个时刻存在一个参考点相关系数的全场分布, 由此可以得到相关系数月平均的空间分布, 如图 7所示。

图 7a中存在三个相关系数0.4等值线区域, 即纯随机数扰动方案的相关系数分布存在虚假相关; 随机场扰动方案的相关系数分布具有较好的局地性, 相关系数随着与参考点距离而增大, 相关尺度大小适宜; 时间滞后扰动方案的相关系数(大于0.4)空间范围涉及到全球海域, 且存在虚假相关(如图 7c、f所示), 其相应的同化调整计算量过大, 不适于实际应用。时间滞后扰动方案中, 同一时刻的扰动样本仅相差一个系数, 即同时刻的扰动集合样本间的相关性较大, 导致该方案有效波高空间相关尺度较大。从同化调整的角度来讲, 随机场扰动方法为最佳方案。

2.3 海浪谱集合样本统计分析

上述2.2节中, 以海浪要素有效波高为例, 统计分析了三种风场集合扰动方案对海浪模式的影响, 但海浪模式运行过程中是以二维海浪谱进行积分计算的, 有效波高仅为海浪谱的积分量, 难以细致地反映海浪的内部结构信息。由此, 为了了解三种风场集合扰动方案对海浪谱的影响, 针对海浪谱集合样本进行统计分析。

选取空间某点, 针对该点某时刻(此处取2014年1月31日24时)的二维波数谱(记为P), 按照公式(7)统计其发散度S_P(二维空间变量), 如图 8所示。与上述发散度对应, 单点海浪谱集合样本平均P的空间分布如图 9所示。

对比图 8和图 9可以看出, 海浪谱集合发散度主要集中于谱峰的位置, 且随与谱峰距离的增大而减小, 说明海浪模式集合扰动反映了海浪谱的主要部分, 并没有破坏海浪谱的主体结构。图 8和图 9中, 红色散点代表集合风场的方向(以下风、浪方向均指去向), (180°W, 0°)和(30°W, 40°N)位置处的集合风场方向与波向一致; (60°E, 30°S)位置处, 三种方案都存在两个谱峰(东北向和西南向), 但三种方案的集合风场方向都集中分布于西南向, 东北方向的谱峰可能是西风带位置的涌浪传播导致的, 这点也可以从相关系数分布图 7a-c中看出, (60°E, 30°S)位置与西风带位置相关性较强。

3 结论和展望

现今, 集合Kalman滤波方法在海洋和大气资料同化领域得到广泛研究和应用, 但具体到海浪资料同化相关研究较少。模式集合样本是否具有代表性是决定集合滤波同化效果的关键因素之一。考虑到海浪模式对初始场敏感性较弱, 本文通过对风场进行持续扰动, 构造海浪模式集合。所构造的海浪模式集合应当具有良好的时空发散度, 集合样本作为统计量应满足一定的统计分布, 此外, 由海浪模式集合统计得到背景误差相关信息应合理可信, 满足上述要求的模式集合样本方可应用于集合Kalman滤波同化研究。

根据集合滤波的基本思想, 本文提出三种风场扰动方案:纯随机数扰动、随机场扰动和时间滞后扰动。利用2014年1月ECMWF风场, 基于MASNUM-WAM海浪模式, 构造风场集合, 开展为期一个月的海浪模式集合运行实验, 以检验风场扰动方案对MASNUM-WAM海浪模式的影响。针对海浪要素(有效波高)和海浪谱进行集合样本统计分析, 结果表明, 随机场扰动方案为最佳方案。随机场扰动方案的发散度在空间上分布均匀, 发散度随时间变化不大, 浮动范围为0.1-0.2m;有效波高集合样本的离散概率密度分布随时间变化的统计表明, 该方案的离散概率密度接近正态分布, 且不随空间点位置的改变而出现较大形变; 有效波高集合样本的相关距离尺度统计分析表明, 该方案的相关系数分布具有较好的局地性, 相关距离尺度适中, 适用于下一步的集合Kalman滤波同化研究。此外, 二维波数谱的集合样本发散度统计分析表明, 风场集合扰动方案没有破坏海浪谱的基本结构, 海浪谱集合发散度主要集中于谱峰附近。

综上所述, 随机场扰动方案为最佳风场集合扰动方案, 即利用Evensen (1994)提出的方法生成随机场集合, 将该随机场集合作为扰动量叠加到风矢量场, 构成集合风场, 驱动海浪模式集合, 可用于后续的集合滤波同化实验。集合Kalman滤波方法是一种较为成熟的同化方法, 下一步工作将把随机场扰动方案与集合滤波方法相结合, 基于MASNUM-WAM海浪模式, 开展海浪资料同化实验, 以检验集合滤波同化方法在海浪资料同化过程中的效果。