一种融合近海多颗卫星海表高度异常数据的方法
【技术领域】
[0001] 本发明公开一种融合近海多颗卫星海表高度异常数据的方法，属于海洋卫星资料 处理技术领域。
【背景技术】
[0002] 海表高度异常（Sea Level Anomaly)是表征海洋状况的一个重要指标。从 1992 年 Topex/Poseidon 卫星发射成功后（Fu et al.，1994, TOPEX/POSEIDON mission overview, J. Geophys. Res.，99, 24369-24381)，卫星海表高度计的观测技术越来越成熟，目 前可以同时得到多颗卫星海表高度异常的观测。图1显示了 2012年1月1日至1月11日 11天内三颗卫星观测的中国近海（10° S-50° N，90° E-160。E)海表高度异常资料的空 间分布，这三颗卫星是Jason-1、Jason-2和CryoSat-2。从图1可以注意到卫星海表高度 异常观测在空间上采样不均匀，沿卫星轨道（沿轨）通常可以达到6km，但轨道间的距离为 200-300km。卫星海表高度资料在时间采样上也有其自身特点，比如Jason-1及Jason-2卫 星轨道重复采样的间隔为近10天，CryoSat-2时间采样上约为35天重复。
[0003] 在科学研宄及社会应用中，时间及空间上采样均匀的海表高度异常资料 会对应用有极大的益处并且方便使用。DUACS (Developing Use of Altimetry for Climate Studies)及 AVISO (Archiving, Validation and Interpretation of Satellite Oceanographic data)是法国的机构，专门从事海表高度异常资料在全球 范围内的处理，产生规则网格点（0.25° X0.25° )、时间分辨率为1天的全球海表 高度异常资料产品 ?ibarboure et al. ，2011，Jason_2in DUACS:Updated system description, first tandem results and impacts on processing and products, Marine Geodesy, 34, 214-241 ；Le Traon et al.,1998,An improved mapping method of multisatellite altimeter data, J. of Atmos, and Oceanic Technology, 15, 522-534)〇 但AVISO资料在产生全球产品时，只用了部分（1/3)的沿轨资料，并且在近海因为使用的 Lanzcos空间滤波方法会造成大的误差。这样便降低了海表高度资料在中国近海的精度，给 资料应用带来一些不便。
【发明内容】
[0004] 本发明针对AVISO全球海表高度异常资料处理中近海资料精度降低的问题，提供 一种新颖、有效、实用的融合中国近海多颗卫星海表高度异常数据的方法。
[0005] 为解决上述技术问题，本发明提供一种融合近海多颗卫星海表高度异常数据的方 法，其特征在于，包括以下步骤：
[0006] 1)获取多颗近海卫星沿轨海表高度异常观测值数据；
[0007] 2)利用经验模分解（Empirical Mode Decomposition)方法进行沿轨海表高度异 常观测值的空间虑波；
[0008] 3)对空间虑波后的数据，利用最优插值方法将不规则空间、时间采样点上海表高 度异常观测值插值到预设空间分辨率、时间分辨率的规则时间、空间网格上。
[0009] 本发明所达到的有益效果：本发明的方法针对AVISO的方法进行了改进，提出了 一种融合中国近海多颗卫星海表高度异常数据的方法，利用三颗卫星的沿轨海表高度异常 观测产生规则空间、时间网格上中国近海海表高度异常产品，与沿轨资料的均方根误差为 3. 49cm，而AVISO产品在中国近海海域与沿轨资料的均方根误差为3. 89cm，。因而本方法产 生的近海海表高度异常数据更接近沿轨资料，将误差降低了 10. 3%。与49个独立的测潮站 相比，本方法产生的近海海表高度异常数据与测潮站资料的均方根误差为10.01cm，AVISO 产品的均方根误差为11. 〇〇cm，因而与测潮站资料相比，本方法的误差降低了 9%。
【附图说明】
[0010] 图1为2012年1月1日至1月11日三颗卫星在中国近海海表高度异常的观测采 样点示意图；
[0011] 图2为经验模分解方法的示意图；
[0012] 图3a为模拟的海表高度异常真实值加上一个正态分布随机误差后的示意图；
[0013] 图3b为用Lanczos滤波方法及经验模分解滤波方法进行空间滤波后的结果示意 图；
[0014] 图4为2012年AVISO海表高度异常资料产品与沿轨资料的均方根误差及本方法 产生的规则网格点上的海表高度异常资料产品与沿轨资料的均方根误差对比图；
[0015] 图5为评价海表高度异常资料产品所用的49个测潮站的地理位置示意图；
[0016] 图6为比较AVISO及本方法生成的海表高度异常资料产品与49个测潮站的均方 根误差示意图。
【具体实施方式】
[0017] 下面结合附图对本发明做进一步的说明。
[0018] 本发明的一种融合近海多颗卫星海表高度异常数据的方法，包括以下步骤：
[0019] 1.获取三颗近海卫星沿轨海表高度异常观测值数据。
[0020] 2.利用经验模分解方法进行沿轨资料的空间虑波。
[0021] 在使用沿轨海表高度异常资料进行空间、时间最优插值前，要对沿轨资料进行空 间上的滤波以去掉一部分观测误差及小尺度噪音（在中炜度为60公里以下）。AVISO中使 用了传统的Lanzcos滤波方法，在中炜度滤掉空间尺度小于40-60km的扰动。这种方法将 沿轨资料进行FFT变换，在谱空间进行滤波，然后将滤波后的谱进行叠加以产生滤波后的 资料。因为这种方法的基函数是正弦及余弦函数，对一般性的信号并不能有效地表达，在沿 轨资料的两端易出现大的误差，这是AVISO资料处理中已被注意到的问题。
[0022] 为解决这一技术问题，本方法采用了经验模分解方法（Empirical Mode Decomposition,也称为 Hilbert-Huang Transform, Huang and ffu,2008, A review on Hilbert-Huang Transform:Method and its applications to geophysical studies, Rev. Geophys.，46, RG2006, doi : 10. 1029/2007RG000228)对沿轨资料进行空间滤 波。
[0023] 本方法将一段连续的沿轨海表高度异常观测值依照数据特点分解为多个经验模 态，然后去掉空间尺度最小的第一模态（在经验模分解方法中，第一模态就是空间尺度最 小的模态），将其它空间尺度大的模态叠加起来求和得到滤波后的沿轨海表高度异常观测 值数据。
[0024] 求取经验模的方法为筛选法，以下仅描述与沿轨海表高度异常观测处理有关的部 分，细节及数学原理可参考Huang和Wu (2008)。图2为经验模分解方法的示意图，x-轴为 空间距离，y_轴为信号的大小。a)显示任意一个随空间变化的信号；b)找出其极大值（菱 型）及极小值（圆圈）；c)用三次样条插值将极大值及极小值连起来，得到极大值包络线及 极小值包络线，求出极大值包络线及极小值包络线的平均，图中以虚线表示；d)求出原始 信号与极大值和极小值包络线平均的差；e)对d)中的差重复步骤b)至d) ;f)当以上过程 结束后，求出原始信号去掉第一经验模后的余项。图a)至图f)表示了求取第一经验模的 过程，重复以上过程可以用几个经验模的叠加来表示原始信号，并且余项为一常数或者只 有一个极值，无法定义极大值和极小值的包络线。
[0025] 将原始沿轨海表高度异常观测值表示为x(t)，t = 0,…，T为空间变量（沿卫星 轨道的距离），按照如下步骤可求得x (t)的经验模：
[0026] (1)找出原始沿轨海表高度异常观测值x(t)中的极大值和极小值，用三次样条 函数将极大值和极小值分别连接起来，得到极大值包络线函数M (t)和极小值包络线函 数m(t),求出所述极大值包络线函数M(t)和极小值包络线函数m(t)的平均值叫，11^ = 0? 5 [M (t) +m⑴]，如图2a，图2b，图2c所示；
[0027] (2)将原始沿轨海表高度异常观测值x(t)减去平均值叫，得到两者差值V
[0028] h1=x(t)_m1
[0029] 如图2d所示；
[0030] (3)对差值h重复以上步骤（1)、步骤（2)，得到第二差值h2，
[0031] h2=h「叫
[0032] 其中m2为差值h i极大值和极小值包络线的平均值，如图2e所示；
[0033] (4)重复以上步骤⑴至步骤（3)，可得到第三差值至第k差值h3，h4，…，匕+匕 等，当如下依据h k_JPhk定义的判据SDk满足，
[0035] 时停止，此时得到原始沿轨海表高度异常观测值x(t)的第一经验模Q，即Qihk， 其中e= 0. 1 ;
[0036] (5)求取原始沿轨海表高度异常观测值x(t)与Ci的余项r i，
[0037] r1=x(t)_C1，
[0038] 如图2f所示；
[0039] (6)按照步骤（1)至步骤⑷求取x(t)的第二经验模C2
[0040] r2=rfC。
[0041] (7)重复（5)、（6)求取原始沿轨海表高度异常观测值x(t)的其他经验模，最终将 原始沿轨海表高度异常观测值 x(t)分解为：
[0043] 其中Q，i= 1，???]!为x(t)的第i个经验模，rn为将x(t)表达成n个经验模叠 加后的余项；
[0044] (8)当rn为常数或者只有一个极值时，步骤（6)至（7)停止。
[0045] 在本方法中，将沿轨海表高度异常观测值进行经验模分解后，将第一模去掉，用其 他模的和来实现海表高度异常的空间滤波。在某些特殊情况下，当连续沿轨海表高度异常 观测值太少时（小于10个点），不进行空间滤波。
[0046] 图3用模拟的沿轨海表高度异常真实值xt_(t)来比较经验模分解滤波及 Lanczos滤波方法。对以下模拟的海表高度异常观测，
[0047]
[0048] 其中 t 为沿轨距离（单位，km)，t = (i-1) X6, i = 1，…，150 ;I\= 90km，T2 = 144km，T3= 600km。上式中的第一项为一个非线性相互作用项，第二项为一个大尺度变化 项，这种形式的选择是为了代表典型的海洋过程，并且x tnK(t)的均方差为7. 84cm，与海洋 中尺度涡旋的空间变率可以比较。
[0049] 在xtrue(t)上叠加一个平均值为0,均方差为2cm的正态分布随机数N(0,2)来模 拟沿轨海表高度资料中的观测误差。这些参数的选择与卫星海表高度资料的误差相一致。 图33表示了义_(1:)与1_(1:)+1^(0,2)随沿轨距离的变化。
[0050] 图3b比较了 Lanzcos及经验模分解法两种滤波方法的结果。使用AVISO中所用 的1^1^(3〇8方法滤掉1_(；(1:)+1^(0,2)中波数高于1/601〇^ 1的小尺度变化。结果表明，使用 Lanzcos方法滤波后，滤波后的信号xlan(x)与x tnM(t)的均方根误差为3. 76cm。使用经验 模分解方法进行空间滤波后，滤波后的信号