齐亚琳 林明森

（1 中国空间技术研究院，北京 100094） （2 国家卫星海洋应用中心，北京 100081）

1 引言

我国第一颗海洋动力环境探测卫星——海洋二号（HY-2）于2011年8月成功发射，可执行业务化监测和全球海洋动力环境参数探测。在HY-2卫星数据产品中，为最大限度地保存观测数据的信息，并节省存储空间，数据产品采用沿轨存储的方式，即沿卫星轨道方向存储风场、海面温度等海洋参数，并同时记录对应的观测时间等信息；并且HY-2卫星搭载的微波散射计和扫描微波辐射计的刈幅宽度虽然大于1 400km，但一天沿轨数据仍不能完全覆盖全球海域。而海洋环境预报等部门在实际使用卫星遥感数据中，通常需要用到全球化网格化的标准数据，并要求数据具有较高的空间和时间分辨率。为了解决上述问题，数据融合被认为是非常有效的手段之一。通过数据融合可以优化遥感信息资源利用，实现多源数据的优势互补。

本文以HY-2卫星获取的海面风场和海面温度作为数据源，结合美国国家环境预报中心（National Centers for Environmental Prediction，NCEP）的再分析数据，分别进行了海面温度和海面风场的数据融合，并得到了数据融合后的结果。通过定性分析，结果表明文中使用的融合方法对数据改善明显，并且可用于探测中、大尺度的海洋现象，证明数据融合技术可为HY-2卫星高分辨率数据产品的开发提供借鉴。

数据融合算法主要包括逐步订正算法、混合分析算法、截断滤波窗的变分分析算法、客观分析算法、三维变分同化分析算法、卡尔曼滤波算法和最优插值算法。逐步订正算法由于仅使用当前的观测数据，因此时间上的记忆性差；混合分析算法也没有时间上的记忆性，且时空分辨率都很低，分析场被严重平滑；统计最优估计算法和使用截断滤波窗的变分分析算法均选用了包括红外与微波在内的多种卫星遥感数据，空间分辨率高，但在有多种卫星遥感数据时，仅选取其中一种数据则可能给融合结果带入较大的噪声；三维变分同化分析算法对计算机存储要求低，计算速度快，三维变分同化分析和最优插值算法分析一样，融合结果受到数据分布稀疏的限制；最优插值算法可看作是卡尔曼滤波法的简化版，与卡尔曼滤波法相比，最优插值算法有以下几个方面的特点：在误差的计算上，从卡尔曼滤波法的公式看，它的误差是发展的，这是卡尔曼滤波法最吸引人的地方，但这一直也是制约着卡尔曼滤波法应用的致命之处，因为它的计算量很大，而最优插值算法完全省略了这点。因此最优插值算法减少了很大的计算量，使得其误差变为静止的；在最优插值算法的资料选择方面，它通常只是选取分析点附近的资料来进行分析。这样，大大减少了计算量；最优插值算法通常只进行单变量分析（例如海面温度）；时间空间距离加权插值的优点是公式比较简单，特别适用于结点散乱、不是网格点的问题。

总的来说，最优插值算法由于其计算量少，融合的“性价比”高，因此可以在单变量数据融合——海面温度融合业务上得到广泛的应用［1］。基于以上考虑，本文选用最优插值算法对海面温度数据进行融合，并且可利用该方法进行工程化和业务化。

根据散射计数据和插值方法本身的研究发现，时间空间权重插值法最适合结点散乱、非网格点的散射计数据融合，并且实现起来也是比较容易，适宜工程化和业务化操作。

2 海洋二号卫星海面温度数据融合

2．1 海洋二号卫星海面温度数据融合方法

2．1．1 最优插值算法

在最优插值算法中，空间网格点上的分析值是由网格点的背景场加上修订值而确定的，其修订值由周围各观测点的观测值与背景场值的偏差加权求得，其权重系数（即最优插值系数）不是任意选择的，应该使得网格点分析值的误差达到最小。对于海面温度，最优插值的表达式为［2］

式中：vam代表海面温度在空间网格点上的分析值；vem代表海面温度在空间网格点上的背景场值；K为权重系数矩阵；v0s代表海面温度在观测点的观测值；ves代表海面温度在观测点的背景场值。当观测点和网格点不重合时

式中：H为插值算子。

合并式（1）和式（2）得vam＝vem＋K（v0s-Hvem）。如果背景场和观测值的空间分辨率不相同，则背景场和观测值所属的权重不同，也即H不同。所以本研究分别用H1，H2作为双线性插值算子，分别对应背景场和观测点，通过该算子分别将网格点的背景场值和观测值插值到分析点上，转换为式（3）。

经过数据预处理后，可直接进行式（3）的运算，式（3）中H1、H2及K分别为双线性插值算子和权重系数矩阵。

权重系数矩阵K的表达形式为

式中：B为背景场误差协方差矩阵，R为观测误差协方差矩阵。欲计算权重系数K，必须先估算背景场误差协方差矩阵B和观测误差协方差矩阵R。

2．1．2 背景场误差及观测误差协方差矩阵计算

对B的确定基于以下两个假设条件：（1）B是定常的；（2）背景场误差的水平相关性满足水平距离的增加相关性呈指数递减的规律。则

式中：ρ为预报背景场水平相关系数矩阵，D为背景场误差组成的对角线矩阵。对于观测误差协方差矩阵R采用与背景场误差协方差矩阵B相同的处理方法。

对于相关系数矩阵ρ，鉴于以往的经验和观测站点横向比纵向稀疏分布的特点［3-5］，对相关系数矩阵做了相应的改进，使其更符合相关尺度的物理规律。相关系数矩阵ρ定义为

式（6）的相关性范围呈椭圆形状。式中：a为经度方向的相关尺度；b为纬度方向的相关尺度；，分别为经度、纬度方向上的距离。在本文中a取200km，b取150km。

2．1．3 矩阵的求解

由于观测点的不规则分布，在求解权重系数矩阵K时，（HBHT＋R）-1的求解过程通常出现误差，因此，应在实际求解时应尽可能避免求逆过程。本文采用求解方程的方法来求解权重系数矩阵K。将式（4）进行改写，两边同时右乘（HBHT＋R），得

对式（7）两边再同时进行矩阵转置，得

因此，将对K的求解转化为对KT的求解，且避免了对矩阵的求逆。但是因为稀疏矩阵H数据元素分布的不规则性（由于观测数据分布不规则导致），得到的转置矩阵方程组可能是病态的。对于病态方程组的求解，通常有三种方法求解：第一种方法是采用病态迭代算法，即用全选主元高斯法求解，得到方程组的一组近似解，将近似解带回原方程，得到剩余向量方程，然后求解剩余向量方程，直到满足精度为止。第二种方法是采用QR 分解与Household转换相结合，转化为求解线性方程组的最小二乘问题，得到线性方程组的最小二乘解。第三种方法是采用奇异值分解（SVD），得到线性方程组的最小二乘解。本文采用奇异值分解方法来处理可能出现的病态方程。将K分解为3个矩阵

式中：U为下三角矩阵；S为对角线矩阵；VT为上三角矩阵。SVD 法可以有效地消除不正常解，并满足方程求解精度。

为了解决计算机内存不足和观测数据数量稀少的问题，采用逐点逐次插值过程。其基本思路是以点为基本计算单位，寻求相关半径内的相关点，并计算其相关性，求得该点的融合值，同时在空间上逐行逐列推进。

2．2 海洋二号卫星海面温度数据融合结果

图1所示为HY-2 卫星扫描微波辐射计2011年12月3日升轨（白天）和降轨（夜间）的逐日全球海面温度（SST）分布，图2为使用最优插值算法制作的HY-2卫星扫描微波辐射计全覆盖的2011年12月3日白天和夜间全球SST 分布。

从图2可以发现，融合后的全球海面温度，不再存在由于轨道覆盖范围不够导致的轨道间空隙，而且全球海面温度的空间分辨率不降低。

融合前，HY-2卫星扫描微波辐射计测量的海面温度不能完全反应全球海面温度的分布情况，特别是不能反应出如黑潮这样的西边界洋流（图1）。经数据融合后，全球海面温度分布非常明显地分为热带、亚热带、寒带等，并且能反应出诸如黑潮等高温海流区域，如图3所示。

图1 HY-2卫星扫描微波辐射计数据2011年12月3日的全球SST 分布Fig．1 SST of the world ocean from HY-2satellite radiometer data on December 3，2011

图2 融合后HY-2卫星扫描微波辐射计数据2011年12月3日全覆盖的全球SST 分布Fig．2 SST of the world ocean from the fused HY-2radiometer data of December 3，2011

图3 融合后的HY-2卫星扫描微波辐射计数据西北太平洋区域SST 分布Fig．3 SST of the west north Pacific from the fused HY-2radiometer data

大洋中的中尺度涡旋有冷涡和暖涡之分，顾名思义冷涡在海面表现为海面温度相对临近海域低、暖涡在海表面表现为海面温度较临近海域高。融合后的海面温度可用来探测全球海洋的中尺度涡旋，而融合前不具备这样的能力。除此之外，融合后的海面温度可用来研究厄尔尼诺与南方涛动（El Niño-Southern Oscillation，ENSO）等海洋灾害，而融合前的数据由于不完整，不具有真正意义上的使用价值。

综上所述，针对HY-2卫星扫描微波辐射计获得的海面温度数据，本文使用最优插值算法进行单变量融合，融合结果在不降低精度的基础上提高了空间分辨率。完全可以在现有的基础上，融合国内外同类或红外卫星获取的海面温度数据，以提高空间分辨率和精度，并实现业务化和工程化。

3 海洋二号卫星海面风场数据融合

3．1 海洋二号卫星海面风场数据融合方法

采用时间空间权重插值法对多源卫星遥感海面风场进行融合。该方法实质上是一种比较简单的时

式中：Uest表示最终获得的估计值；下标k（k＝1，2，…，n）表示卫星观测数据点，Wk表示在k点的权重，Uk表示在k点的观测值，（xk，yk，tk）表示k点的时空坐标；下标o表示网格点，（xo，yo，to）表示其时空坐标；n表示插值点个数；R为距离影响半径，T为时间影响半径。

本文选定的时间空间权重插值法为普通Kriging插值法进行插值。Kriging 插值法是法国科学家Matheron［6］提出的，是对空间分布的数据求线性最优、无偏内插估计的一种方法，并在获得预测结果同时可以获得预测误差。其基本步骤如下：

1）求变异函数

变异函数r（h）反映了区域化变量的空间自相关性。

式中：h为观测点之间的空间间隔距离；N（h）为距离等于h的点对数；Z（xi）为处于点xi处变量的观测值；为与点xi偏离h处变量的实测值。在实际计算中，也可以将所有观测点的相对距离划分为若干等级，计算每个等级内观测点的个数，然后对每个等级内的所有点数取距离的平均值及r（h）的平均值。将计算所得的（h，r（h））点连接后就可以得到实验变异函数，再以最小二乘法计算出理论变异函数及其参数。理论变异函数模型有线性模型、指数模型、球面模型、高斯模型等。通过比较选取的高斯模型作为理论变异函数。

式中：C0称为块金常数；C0＋C为基台值；，C为拱高。

2）利用变异函数求权重系数

假设为未观测点，xi（i＝1，2，…，n）为其周围的观测点，Z表示计算的区域化变量。则对处某个区域化变量的估计值为

式中：λi表示权重，可由理论变异函数求得，为了保证是无偏估计，要求＝1。式（13）也可以改写为矩阵形式

其中

式中：r（x1，xn）为点x1和xn之间的变异量，可由理论变异函数根据两点间距离求出。同时该点估计值误差的平方可由下式求得δ2（x0）＝BTW-1B。

3．2 海洋二号卫星海面风场数据融合结果

图4和图5 分别是同时间内的HY-2 卫星微波散射计风场和美国国家环境预报中心（NCEP）再分析数据。该轨道的HY-2 卫星微波数据不能覆盖东海和南海，而NCEP 数据可覆盖全球范围内的海域。故在实际应用中需将高分辨率的HY-2风场数据和不受任何轨道等影响的NCEP数据进行融合。图6是图4 和图5 所示的HY-2 微波散射计风场数据和NCEP数据的融合结果。图6 显示的融合产品中，风场的空间分辨率达到与HY-2风场空间分辨率相同，并且覆盖东海和南海等海域。HY-2卫星微波散射计一天最多可覆盖全球90%的海域，而融合后的产品可100%覆盖全球海域。

图4 用于融合的海洋二号散射计风场Fig．4 HY-2satellite scatterometer wind field to be fused

图5 用于融合的NCEP风场Fig．5 NCEP wind field to be fused

图6 西北太融合海面风场Fig．6 Fused wind field in the west north Pacific

从图6可以看出，融合后的海面风场不仅弥补了高分辨率的HY-2卫星遥感数据不能全覆盖的不足，同时提高了NCEP 数据的空间分辨率，可为两种数据更好地使用，尤其是为业务化海洋监测应用提供便利。

另外，从图6可以看出，融合后的海面风场在东海黑潮上空加速，风向平行于东海黑潮主轴，这与已有的研究结果［7-10］一致。虽然海面温度对海面风速影响机制目前尚不完全清楚，但融合后的海面风场可为研究这一机制提供非常好的数据源。

目前国际上同类卫星微波散射计数据非常稀少，但在不久的将来同类卫星数据会逐渐增多。同时，目前国际上有较为成熟的风场再分析数据，如NCEP，欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）的数据等，可以用来和HY-2卫星散射计数据进行融合，并实现业务化和工程化。

4 结束语

本文针对HY-2卫星获取的海面温度和海面风场数据，利用不同的融合方法对HY-2卫星海面温度和海面风场数据进行融合。通过定性分析，融合后的海面温度和海面风场数据在空间分辨率不降低的前提下，弥补了自身的不足。

本文的研究结果证明数据融合方法应用到HY-2卫星数据处理系统中，可以较好地解决卫星观测数据覆盖范围不足的问题。这为今后卫星数据的高效使用和海洋环境预报部门更有效地使用卫星数据，提供了一种技术途径。

从本文的研究结果可以发现，HY-2 卫星数据经融合后得到的海面温度、海面风场数据，可为全球海气相互作用研究提供同时间内相同区域的数据源，对进一步找到海面温度和海面风场的定量物理关系奠定基础。同时，HY-2卫星海面温度和海面风场数据融合的成功实现和广泛应用，可为HY-2卫星其它数据的高效使用提供技术基础，并实现HY-2卫星与国内外同类卫星数据融合的业务化和工程化。

（References）

［1］奚萌．基于最优插值算法的红外和微波遥感海面温表温度数据融合［D］．北京：国家海洋环境预报中心，2011

Xi Meng．Merging infrared radiometer and microwave radiometer sea surface temperature data based on the optimum interpolation［D］．Beijing：National Marine Environmental Forecasting Center，2011（in Chinese）

［2］Reynolds R W，Smith T M．Improved global sea surface temperature analyses using optimum interpolation［J］．J．Climate，1994，7：929-948

［3］冯天祥，李世宏．矩阵QR 分解［J］．西南民族学院学报（自然科学版），2001，27（4）：418-422

Feng Tianxiang，Li Shihong．QR factorization of the metrix［J］．Journal of Southwest University for Nationalities（Natural Science Edition），2001，27（4）：418-422（in Chinese）

［4］马寨璞．海洋流场数据同化方法与应用的研究［D］．杭州：浙江大学，2002

Ma Zhaipu．Study of data assimilation method and its application to sea flow field［D］．Hangzhou：Doctoral dissertation of Zhejiang university，2002（in Chinese）

［5］张识．力法中病态方程的产生与改进方法［J］．力学与实践，1989，11（4）：57-59

Zhang Shi．Generation and improving method of the morbid equation in the force method［J］．Mechanics in Engineering，1989，11（4）：57-59（in Chinese）

［6］Matheron G．Principles of geostatistics［J］．Economic Geology，1963，58：1246-1266

［7］Pan Jiayi，Yan XiaoHai，Zheng Quanan，et al．Observation of western boundary current atmospheric convergence zones using scatterometer winds［J］．Geophys．Res．Lett．，2002，29（17）：1832

［8］Chelton D B，Schlax M G，Freilich M H，et al．Satellite measurements reveal persistent small-scale features in ocean winds［J］．Science，2004，303（5660）：978-983

［9］Chelton D B，Xie S P．Coupled ocean-atmosphere interaction at oceanic mesoscales［J］．Oceanography，2010，23（4）：52-69

［10］Xu Haiming，Xu Mimi，Xie Shangping，et al．Deep atmospheric response to the spring Kuroshio over the East China Sea［J］．J．Climate，2011，24：4959-4972