高度计波高数据同化对印度洋海域海浪模式预报影响研究

齐鹏，范秀梅

（1.中国科学院海洋研究所海洋环流与波动重点实验室，山东青岛 266071；2.中国科学院大学，北京 100039）

1 引言

近年，我国海洋环境业务预报部门已建立基于第三代海浪模式的全球海浪数值预报业务化系统。但是，目前在海浪预报工作中，仍未能有效利用卫星高度计波高测量数据进行海浪的同化预报。从另一方面看，越来越多的海洋卫星通过高度计提供全球大覆盖量、高精度和实时的海洋表面测量数据，使得将高度计实时波高观测数据同化到海浪模式以改善模式预报初始场，进而提高海浪数值预报精度成为可能。

如上所言，海浪观测方面，当前海洋卫星观测已成为海洋环境立体观测的主要手段之一。1978年，美国国家航空航天局（NASA）发射了具里程碑意义的海洋卫星SeaSat-A，它发回的数据通过处理可获得包括海面风速、风向、波高、波长、波谱、海面温度、大气水含量、海冰、海面地形、海洋水准面和高分辨率雷达图像，而且测量精度基本达到实用要求。其后，主要有欧空局（ESA）1991和1995年先后发射的资源卫星ERS-1和ERS-2；1992年NASA与法国国家空间研究中心（CNES）合作发射的TOPEX/Poseidon（简记为T/P），其后继卫星Jason-1于2001年底发射；2008年发射的Jason-2则是T/P和Jason-1的后继卫星，它的发射出于以下两个目的：一是为海洋研究提供高质量的观测，二是为同化和预报提供业务产品。另外，值得特别指出的是，2012年4月起我国第一颗海洋动力环境卫星海洋二号（HY-2）已开始向社会发布海洋环境监测数据产品。

海浪模式方面，WAVEWATCH IIITM是由美国NCEP（National Centers for Environmental Prediction）的 EMC(Environmental Modeling Center)下面的MMAB（Marine Modeling andAnalysis Branch）在WAVEWATCHⅠ和WAVEWATCHⅡ的基础上开发的全谱空间的第三代海浪模式。其不同之处主要体现在控制方程、程序结构、数值和物理方法方面，允许用户用自己的数值或物理方法发展一个基于WAVEWATCH IIITM框架的新的海浪模式；同时，WAVEWATCH IIITM框架下的优化方案、并行化算法、嵌套方案、输入和输出功能能方便地得到共享和使用。此外，本文所用版本WAVEWATCH IIITMversion 3.14提供多网格相互作用即网格嵌套的模式方案，使得在处理有开阔水域边界时，能使用网格嵌套方案为开边界提供边界条件。WAVEWATCH IIITM没有提供可用的海浪同化代码（虽然NCEP表示已在做海浪同化技术方面的工作，但目前并不打算发布海浪同化模块），但提供了海浪同化的接口，这就方便我们将自己开发的同化程序植入模式和随模式的积分运算而被调用，并将同化结果反馈到模式中，为模式预报所用。

只使用高度计波高数据的同化方法是由Esteva[1]和Janssen等[2]较早提出。Lionello等[3]认为Janssen等[2]的同化方法太复杂，提出了更简明的方法。其方法分两步：第一步，以实测有效波高（SWH）做最优插值分析；第二步，由浪高分析场求总波能，然后利用总波能和初估模式谱得到最终的海浪分析谱。他们是利用插入观测法或最优插值法（Optimal Interpolation，OI）将观测值同化到WAM（Wave Modelling）模式中。Mastenbroek等[4]，Breivik 和Reistad[5],Lionello 等[6]，Bender和 Glowacki[7]以 及Greenslade[8]应用类似方法开展同化海浪预报，均收到不同程度改进海浪预报的效果。国内方面，张志旭等[9]使用一种简单的逐时同化的方法对海浪模式进行高度计SWH数据同化试验；王毅和余宙文[10]将最优插值和海浪谱重构技术与SWAN模式结合，建立海浪最优插值同化预报模式。

总之，为提高海浪模式预报精度和延长其预报时效，当前，人们有条件大量使用高度计海浪观测资料，将其同化到海浪模式以提高模式预报初始场的质量，进而提高海浪预报效果。

根据我国海洋公益性行业科研专项的要求，结合当前我国海上能源运输安全和海军护航任务迫切需求，特别是围绕关于重点开展海洋环境预报中关键技术的研究，建立印度洋海域海洋环境数值预报系统，提升我国远洋环境保障能力，我们研制了基于OI方法的高度计SWH数据同化并行模块，并将其植入第三代海浪模式WAVEWATCH IIITMversion 3.14。通过同化Jason-2高度计SWH数据，实现模式边积分边同化的功能，进行印度洋目标区域海浪场短期预报，并对同化提高模式预报初始场的质量和同化对短期预报影响的效果和时效性进行检验。

2 高度计有效波高OI同化模块研制

最早，Eliassen[11]利用最小二乘法推导观测值和背景场向量的OI方程。在对背景场和观测误差性质充分了解的基础上，OI方法充分利用背景场和观测值的信息，给出一个方差最小意义下的最优线性估计。式（1）和（2）为OI同化表达式。

式中，各下标代表矩阵的维数，n表示所有点数，p表示观测点数；上标a代表分析场，也即加入同化后的模式波高场；b代表背景场，也即同化之前模式给出的（待同化的）波高场；o代表观测场；矩阵H代表有效波高；P是观测算子，它是一个插值算子（本文的插值方法为双线性插值），作用于背景场Hb，可以得到观测值的第一猜值PHb；W代表最优权重矩阵；R代表观测误差协方差矩阵；B代表背景误差协方差矩阵。

对于观测误差协方差矩阵R，与国内外同类研究的做法一样，这里R中元素的值取为：

式中，δkj是kronecker符号；σ0的取值，根据Jason-2产品手册中给出的高度计有效波高误差为5%或0.25 m，同时考虑文中所用观测时次的有效波高大都在5 m以内，故取σ0=0.25。易知，这里R事实上是一个p×p维的常值对角方阵，即不同点上观测值之间是不相关的。

对于背景误差协方差矩阵B，与国内外同类研究的做法一样，根据相关长度法，其元素的值采用下面的经验表达式：

式中，σp是模式预报的背景场的均方根误差（RMS），并参考已有的业务化预报结果，这里σp取为0.6 m；dkj是观测点k与 j之间的距离，以度为单位。L为经验常数，根据相关文献，取L=8Sx，Sx为模式的空间网格分辨率，本文取为1/3°，故L约为2.67°；换算成以km为单位，则L在300 km左右。通过试算同化效果，将L值调整到相对较为合理。

计算背景误差时，还需考虑两点的背景方差能相互影响的最大距离，即影响半径的问题。文中影响半径r的取法为r=11Sx，即约为400 km。最后，B矩阵中的元素由下式得到：

由式（5）可知，两点之间的距离如果大于影响半径r，就认为它们之间的背景方差不相关。在OI同化方法中，观测误差协方差和背景误差协方差对同化效果至关重要。因而，不排除为了同化的效果更好，不同海区在参数选取上可有不同，一般通过试算，可以获得一个较合理的参数设置。

在处理占用内存较大的背景误差协方差矩阵时，文中采用稀疏矩阵的表示方法，对稀疏矩阵的计算则采用专门处理稀疏矩阵的Sparse Blas函数库里的函数，以提高运算速度。

完成上面工作后，接下来要做的是使编写好的OI同化模块与WAVEWATCH IIITM模式实现对接。值得指出的是，这里开发的同化模块也是按并行设计的，这样就与要植入其中的海浪模式在程序并行风格上取得了一致。

WAVEWATCH IIITM由模块ww3_shel控制整个模式的运行（见图1），而w3wavemd是模式的核心计算模块；调用同化模块是在模式积分到目标时刻，输出结果之前，被调用；同化完有效波高数据后，调用输出模块，输出同化后的模式结果（见图 1）。

作为全谱空间模式，WAVEWATCH IIITM模式中的海浪是以二维波谱形式出现，而被同化的是高度计的有效波高，由此需要在每次同化过后利用同化后得到的有效波高分析场重构出相应的二维海浪谱。同化前与同化后的波能表达式如下：

式中，i、j为网格点坐标；上标b和a分别代表同化前和同化后；Hij代表有效波高，Fij(f,θ)代表海浪二维波谱；f是海浪频率，θ是波向。在有效波高的变化不会引起波能在频率和方向上的重新分布，而只会引起波能在谱空间中相同的线性变化的假定下，由上面同化前与同化后的波能表达式可以得到同化后的海浪谱调整形式为：

在使用Jason-2高度计有效波高（SWH）数据时首先需要对数据进行质量控制，即剔除不合理的异常数据。此外，我们在同化数据选取时规定，取计算范围内有观测数据的时刻作为同化时刻，以此时刻为中心，在其±0.5 h时间范围内取观测数据，作为该时刻同化所采纳的数据，这样基本上可保证同化时刻观测值与模式输出值在时间上取得一致，从而保证同化效果。本文同化2010年12月15—17日期间Jason-2在北印度洋海域上空轨道线上高度计测量SWH（见图5），得到同化初始场，并以此做2010年12月18—20日3天的海浪预报。前后6天对应的卫星轨道周期号是Cycle90。

图1 植入同化模块的WAVEWATCH IIITM流程图

图2 网格嵌套模式方案

3 高度计SWH同化初始场的检验

为检验上述高度计SWH数据同化模块用于WAVEWATCH IIITM得到的同化初始场对比无同化的情形，对模式初始场是否有改善以及改善程度，设计了以下数值试验。以5°S以北印度洋海域作为目标计算区域，即5°S—25°N，40°—105°E；空间网格分辨率取1/3°×1/3°；海浪谱网格分辨率取24个方向（即方向分辨率为15°）,和25个频谱段，即按公式fm=αfm-1，这里，m=1,2,···,25 ，初始频率 f0取0.0418，频率增长因子α取1.1。海面大气强迫采用业务预报单位提供的WRF（Weather Research and Forecasting）模式的海面10 m风场产品，1 h 1次的。

考虑目标计算区域南边界连接开阔大洋，极易受传入大洋涌浪信号影响，采用全球网格嵌套印度洋区域网格的方案（见图2），为其提供开边界条件。全球网格计算范围：78°S—78°N，180°W—180°E；所给空间网格分辨率相对较低：1.0°×1.25°；海面大气强迫采用6 h1次的NCEP再分析风场。全球网格模式中的频谱和方向谱参数设置与嵌套其中的区域网格的一致。

WAVEWATCH IIITMversion 3.14提供了多网格嵌套的模式运行功能。采用嵌套方案后，由ww3_multi模块替代ww3shel（单网格模块）管理和调用其它的函数，完成多网格嵌套后的模式积分。

2010年12月15—17日这3天在目标计算区域有Jason-2高度计测量数据的时刻列于表1，数据同化即发生在这3天中有观测数据的那些时刻。

首先，我们检验利用同化数据生成同化分析场的效果。我们将同一时刻同化模式的结果（SWH）、无同化的模式结果（SWH）分别与高度计测量SWH进行比较（见图3）。图3为沿轨比较，容易看出，同化模式的结果与高度计SWH更接近。

表1 用于同化的高度计波高数据所在时刻

表2 同化和无同化的模式SWH分别与高度计测量SWH之间的统计检验(2010年12月15—17日)

图3 不同时刻同化的和无同化的模式结果（SWH）与高度计测量（SWH）的比较

为定量评估同化改善模式预报初始场的效果，引入以下统计量，

均方根误差：

相关系数：

式中，Pi代表模式输出值，Oi代表观测值，N为样本个数。

为简单起见，我们将15日、16日和17日各日在观测时刻（见表1）的数据（包括同化和无同化的模式输出SWH及高度计SWH）按所在的日期分别并入一个统计组，进行统计特征量的计算，计算结果列于表2。由表2容易看出，同化模式输出的SWH与高度计SWH之间的均方根误差较无同化的情况明显地减小，而相关系数明显地增大。这表明，高度计SWH数据同化模式的初始场与观测更为接近，表明所研制的高度计波高数据同化模块是有效的，能够起到改善模式分析场的作用。

4 同化初始场影响预报结果的时效性

海浪模式预报对于时间跨度太长的初始场是不敏感的。为考察同化初始场影响短期海浪预报结果的时效性，取轨道上数据连续，数据较多的同化时刻，且其后不能再有与其相交或离它很近的同化轨道，以免影响到检测效果。我们选择2010年12月17日14时，以同化和无同化所形成的初始场（见图4中17日14时的同化初始场与高度计SWH的沿轨分布更接近）为起点，对其后6 h、12 h、24 h、36 h、48 h、60 h和72 h等各预报时间点上输出的预报结果进行比较（见图4），其中，还将6 h前（即17日14时）的沿轨高度计SWH拿来与6 h预报的同化预报结果和无同化的预报结果进行比较，以考察它们何者更接近于观测。

图4 同化初始场对预报影响时效性的试验结果

图5 Jason-2卫星轨道（2010年12月15日00时—20日23时时间段）

表3 预报阶段高度计波高数据所在时刻

显然，对于6 h预报（即17日20时），同化初始场的6 h预报结果较无同化的预报结果更接近6 h前的高度计观测。随着预报时段延长，对后期预报结果的影响，同化模式与无同化模式之间趋于减小；从60 h以后的预报结果来看，同化预报和无同化预报之间差别基本消失。另外，从这一组图来看，无同化的模式预报结果普遍偏大，特别是在波高较小区域。

5 高度计SWH同化对北印度洋海域海浪预报影响评估

为评估高度计SWH数据同化初始场对印度洋海域海浪预报影响，我们以连续同化12月15日、16日和17日3天的该海域高度计SWH数据得到的同化分析场作为预报初始场，进行12月18日、19日和20日3天（即从6 h直到72 h）的印度洋海域短期海浪预报。前3天的高度计测量数据被同化到预报模式，形成同化分析场和预报初始场，后3天（即预报阶段）的高度计数据则不再参加同化，而是用于对预报结果的检验评估。上述6天（即12月15—20日）有卫星观测值的时刻详见表1和表3。

图5粗实线表示用于同化的高度计数据的卫星轨道位置，细虚线表示用来作检验的数据的轨道位置。该图下方是各条轨道所对应的具体时刻，如15.01表示15日01时刻。这些轨道分别就是2010年12月15日、16日、17日、18日、19日和20日各日各时刻（以此时刻为中心，并在其±0.5 h时间范围内取观测数据）经过该目标区域的轨道。由于Jason-2卫星在这6天之内在该区域没有重复的轨道，一天之内又分为下行线和上行线，所以同化之后，在同一个同化位置没有新的高度计观测数据可用来检测同化后的预报效果，因而检验预报效果只能利用其附近的轨道上的观测数据。相邻的用来同化和用来检验的轨道如果同为上行线或下行线则它们是相互平行的；如果一根为上行线而另一根为下行线则它们是相交的；但是，它们之间的时间差却都不小。举例来说，2010年12月18日01时刻的轨道是下行线，而与其（在时间上）最近的下行线是15日01时刻和15日02时刻的，与其相交的上行线是15日14时刻的，故它们之间时间上相距都已超过72 h。由此，在同化时刻过去72 h之后再来检验同化对预报影响的效果其实是很困难的。因为，模式经过1 h一次的风场驱动，3天以后的同化与无同化的风浪场结果之间基本就没有什么差别了。另外，由于我们在该目标区域目前还难以获得浮标等其它来源的观测数据，所以，为检验同化预报的效果也只有利用高度计测量数据了。故在文中还是给出了同化预报和无同化预报在有高度计观测时刻的SWH比较图，但因为用于同化的与用于检验的二者时间上相距太长，同化虽有优势但已不明显。图6中，每幅图下方括号中所标注的是与该验证时刻较近的用来同化的观测数据其轨道所对应的时刻。可以看出，它们在时间间隔上都很长，都超过了48 h，有的还超过了72 h。但从这组图中还是可以看出加入同化的预报要比无同化的效果好，不过随着积分时间的延长，同化与无同化效果渐趋一致。

6 结论与展望

海浪的模式预报对于时间跨度太长的初始场并不敏感，而大气模式提供的即时风场又总是存在着一定的误差。因此，为提高海浪模式预报精度，及时更新和提高海浪模式预报初始场的质量是解决途径之一，即在海浪模式的积分运算过程中应充分利用当前卫星遥感技术所能提供的观测数据，尽可能大面积吸收高度计测量海浪数据并利用一定的同化手段及时更新海浪分析谱，如此可为模式预报提供更接近准确的初始场，从而提高海浪模式预报的精度。

图6 预报阶段同化、无同化预报结果比较（包括与最近的高度计观测的比较）

研制了基于OI方法的高度计波高数据同化并行模块，并将其植入WAVEWATCH IIITM，使同化模式仍能保证其以并行方式运行。为解决同化背景场误差协方差矩阵占用内存较大问题，采用稀疏矩阵的存储方法节省内存和使用Sparse Blas函数库中的相关函数，使该问题得到较好解决。

以5°S以北的印度洋海域作为目标计算区域，嵌套在WAVEWATCH IIITM的全球网格中，较好解决了目标区域的开边界条件问题。所建立的海浪同化预报模式由大气模式WRF输出的一小时一次的海面10m风场驱动运行。同化模式的结果（SWH）和无同化的模式结果（SWH）分别与高度计测量波高（SWH）比较，显示高度计波高数据同化可明显改善模式预报初始场。

进行印度洋目标区域海浪场72 h以上时段同化预报试验，初步分析表明，同化高度计SWH数据进入海浪模式进行边同化边预报，能起到改善海浪短期预报效果的作用；但同化影响模式预报的时效性，限于72 h内；预报时效越短，同化模式的结果与无同化的模式结果之间差别越显著。

总之，本文初步研究表明，连续同化高度计SWH数据进行边同化边预报可望明显提高海浪短期预报结果的精度。

在强调高度计波高数据同化改善海浪模式预报水平的同时，我们也注意到，由于卫星观测所获得的仅是星下点数据，因此，其轨道重现频率与空间覆盖率一直如同一对矛盾。若要追求空间上的高覆盖率，则其轨道重现频率必然受到影响。例如，Envisat卫星两相邻轨道之间最宽距离仅80 km，但其轨道完全重合一次的周期却需35天之久；Jason-2，Jason-1的完全重合周期为10天，但两相邻轨道间的距离却宽达315 km。总之，在本文的研究中不难发现，仅依靠单卫星提供的波高数据是相当有限的，而且无法做到时空同步（实时提供大范围海浪场的高度计测波资料），难以满足进行大洋区域海浪场同化预报时对大面积吸收实测海浪资料的渴求。因而，今后应考虑利用多卫星融合的高度计波高资料进行大洋海浪场的同化预报，应能至少在一定程度上弥补上述不足。例如，利用最近的短时段可获得的Cryosat-2，Jason-2，Jason-1以及Envisat等4颗卫星（随着新卫星上天，将不断有新的卫星高度计资料作为数据源加入）在内的最少两颗卫星的数据，通过一定的滤波方法进行空间平滑和外插处理之后，得到准实时的多卫星融合的测波资料用于大范围海浪场的同化预报。

[1]Esteva D C.Evaluation of preliminary experiments assimilating Sea sat significant wave heights into a spectral wave model[J].J GeophysRes,1988,93(C11):14099-14105.

[2]Janssen PA E M,Lionello P,Reistad M,et al.Hind cast and data assimilation studies with the WAM model during the Sea sat period[J].J Geophys Res,1989,94:973-993.

[3]Lionello P, GÜNTHER H, Janssen PAE M.Assimilation of altimeterdata in a global third-generation wave model[J]. J Geophys Res,1992, 97(C9): 14453-14474.

[4]Mastenbroek C,Makin V K,Voorrips A C,et al.Validation of ERS-1 altimeter wave height measurements and assimilation in a North Sea wave model[J].Global Atmos Ocean Syst,1994,2:143-161.

[5]Breivik L A,Reistad M.Assimilation of ERS-1 altimeter wave heights in an operational numerical wave model[J].Weather Fore-cast,1994,9:440-451.

[6]Lionello P,Guenther H,Hansen B.ASequential data Assimilation Scheme Applied to the Global Wave Analysis and Prediction[J].J MarSyst,1995,6:87-107.

[7]Bender LC,GlowackiT.The assimilation of altimeter data into the Australian wave model[J].Aust Meteorol Mag,1996,45:41-48.

[8]Greensladed J M.The assimilation of ERS-2 significant wave height data in the Australian region[J].J Mar Syst,2001,28:141-160.

[9]张志旭,齐义泉,施平，等.卫星高度计波高资料的同化实验分析[J].海洋学报,2003,25(5):21-28.

[10]王毅,余宙文.卫星高度计波高数据同化对西北太平洋海浪数值预报的影响评估[J].海洋学报,2009,31(6):1-6.

[11]EliassenA.Provisional report on calculation of spatial covariance and autocorrelation of the pressure field[R].Institute for Weather andClimateResearch,AcadSci,Oslo,ReportNo.5,1954.