多源海洋遥感叶绿素数据融合

崔建勇，刘晓东，岳增友，李连伟

(1.中国石油大学(华东) 地球科学与技术学院，山东 青岛 266580；2.海洋科学与技术国家实验室 海洋矿产资源评价与探测技术功能实验室，山东 青岛 266071)

0 引言

地球表面70%的面积被海洋所覆盖，海洋为人类生活提供了丰富的鱼类资源。随人类社会技术的不断进步，特别是遥感技术的快速发展，人类与海洋已融合为一个生命共同体。1978年美国宇航局(National Aeronautics and Space Administration，NASA)发射了Nimbus-7(雨云)号卫星，揭开了利用遥感技术探测海洋的新篇章。之后，NASA和其他各国卫星发射机构相继发射了多颗海洋卫星，实现了大区域、高频率、高精度的海洋观测，为科研人员研究海洋提供了海量遥感数据，对研究海洋生态环境演变、生物地球化学循环、气候变化、渔业捕捞等具有重要的意义[1-2]，也为探测海洋叶绿素α提供了坚实的数据基础。

不同的海洋卫星遥感传感器具有不同的空间分辨率和时间分辨率，光谱和覆盖率也各有不同，故单个海洋叶绿素α浓度数据源存在数据覆盖率、分辨率和数据可利用率等方面的缺陷。为弥补该缺陷，科研人员提出了多源数据融合技术。多源数据融合技术可实现多遥感数据信息互补，可扩展单个海洋遥感的空间覆盖率和时间分辨率[3]，增加整体的信息量，提高数据产品的时空连续性、一致性和可靠性。

多源数据融合技术最早被广泛应用于海表温度、卫星高度计测高数据等方面，随海洋遥感卫星的发展，该技术被用于海色数据的融合。1997年，NASA正式提出了SIMBIOS计划，同时欧空局(European Space Agency，ESA)提出了GLOBcolour计划[4]，相关计划大多使用平均法、GSM(Garver-Siegel-Maritorena)生物光学模型等方法，将多颗海色传感器的数据进行融合，建立了高质量、高精度、高覆盖率的海色数据集。李新星等[1]指出，在2005年GSM 方法提出后，Maritorena等人在2010年将该方法用于SeaWIFS、MODIS、MERIS的数据融合。国内学者在数据融合方面稍落后于国外学者，较早的研究是张灵凯等[5-6]用小波分析方法对SeaWIFS和MODIS的叶绿素α数据进行融合。

通过研究分析，国内外相关科研机构生产的海洋叶绿素α融合产品存在精度低、覆盖率低、时间跨度短等问题。本文基于前人研究成果，利用SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5个传感器的叶绿素α数据，实现了小波变换与Kalman滤波技术相结合的融合方法，构建了1998—2017年全球海洋表面的叶绿素α数据集，并与实测数据进行了对比分析。

1 数据获取与预处理

1.1 多源遥感数据

本文使用的遥感数据为SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5个传感器的叶绿素α浓度数据，时间范围为1998—2017年，数据产品为Level 3标准网格数据，数据格式为Netcdf(network common data form)和HDF4(hierarchical data format)，空间分辨率为4 km和9 km，如表1所示。相关数据均可从美国NASA的oceancolor网站免费下载。

表1 5个传感器叶绿素α浓度数据参数

因5颗卫星发射时间不同，因此各个卫星数据的时间覆盖范围也不同。考虑到融合产品的质量和精度，需要对同一时间段内的数据进行融合。

1.2 实测数据

实测叶绿素α浓度数据用于检验融合产品的精度和可靠性。实测数据共获取了5 710个观测文件。观测文件记录了站点编号、日期、时间、纬度、经度、水深、叶绿素浓度等，数据量共1.53 GB。下载网址为http：//seabass.gsfc.nasa.gov。

1.3 数据预处理

为方便后期融合算法的流程化计算，需要将5个传感器的数据转换成相同的格式和相同的空间分辨率。SeaWIFS、Terra-MODIS、Aqua-MODIS和VIIRS传感器的叶绿素α浓度数据格式为Netcdf格式，而MEIRS传感器的叶绿素α浓度数据是HDF4格式。空间分辨率方面也存在差异，SeaWIFS传感器的叶绿素α浓度数据的空间分辨率是9 km，其他4个传感器的叶绿素α浓度数据空间分辨率是4 km。

为解决上述2项问题，采用ENVI-IDL语言和Python语言编写数据预处理程序，分别对叶绿素α浓度数据进行处理。处理完成后，5个传感器的叶绿素α浓度数据的格式均为TIFF，空间分辨率均为4 km，地理空间投影一致。

2 数据融合方法研究

融合算法以小波变换为基础，并结合Kalman滤波技术。小波变换算法中正交或双正交多分辨率小波融合方法已用于多传感器图像信息的融合，该方法最大限度地保留了原多光谱图像的光谱信息[7-8]，并且能将图像分解为具有不同空间分辨率、频率特性和方向特性的子图像序列。其分频功能相当于一对共轭镜像正交滤波器组，把图像解构为1个低频带子图和3个高频带子图(包括水平、垂直和对角线3个方向)[9-11]。低频部分反映的是图像的整体视觉信息，高频部分反映的是图像的细节特征，具有多尺度、多分辨率、方向性和无冗余数据等优点[12-13]。Kalman滤波是一种基于协方差递归思想的预测器，能够寻求最优的预测结果，效率较高，将其应用于对高频信息的处理，可以生成质量更好的高频信息。

经数据预处理，可得到进行融合处理的标准图像。根据时间范围，将其划分为6个时间段进行融合，如表2所示。

表2 融合时间分配表

从表2可知，1998—1999年遥感数据只有SeaWIFS一种数据源，无法采用本文提出的融合方法进行处理，可采用查找表法进行融合；2000—2017年均有多个数据源，根据数据源个数的不同，研制不同的小波变换融合算法。

2.1 融合权重确定

不同数据源的数据精度不同，融合时会赋予不同的权重值。融合系数计算原理是：数据源对应的误差越小权重越大；一致性测度数值越大，权重越大。结合这2个信息，计算一致性测度，来确定低频系数矩阵的融合权重值。

假设有n个传感器，这些传感器获取同一观测范围的地物属性。设传感器i在K时刻的采集值为Xi(k)(k=1，2，…，s；k表示像素编号；s为观测范围像素个数)，传感器j在K时刻的采集值为Xj(k)，时间为K时，传感器i与j在k位置测量值的支撑力度矩阵如式(1)所示。

(1)

时间为K时，像素位置k的各传感器观测结果总支持度如式(2)所示。

(2)

总支持度矩阵中每一行的支持度之和，代表了对应传感器对最终观测结果的支持度。因此，时间为K时，像素位置k各传感器观测结果支持度如式(3)所示。

(3)

式中：一致性测度ri(k)仅反映了在k像素位置，传感器i的测量值与所有传感器测量值的相近水平。单个像素ri(k)的值，并不代表在整个观测区域上传感器的可靠性，应考虑到整个观测范围上传感器节点的可靠性。

第K时刻第i个传感器观测支持度的均值和方差分别如式(4)、式(5)所示。

(4)

(5)

若某传感器观测范围内像素平均值较大，且灰度的波动较小，则该传感器较稳定，具有较高的支持度，权重较大。因此，可用信噪比描述此传感器的支持度，信噪比的值为平均值除以信号波动的方差。时间为K时，传感器i的权重wi(K)如式(6)所示。

(6)

归一化后值如式(8)所示。

(7)

2.2 小波变换融合与Kalman滤波技术

对数据源计算权重后，利用小波变换对数据进行加权融合。融合算法对2个或多个数据源进行小波变换。首先采用小波多分辨率分解和Mallat快速算法，将原始图像分解成近似图像和细节图像，其分别代表了图像的不同结构；然后在各层的特征域上进行有针对性的融合。因为比较容易提取原始图像的结构信息和细节信息，所以融合效果较好，并且由于小波变换具有完善的重建能力，故保证了信号在分解重建过程中，很少有信息损失和信息冗余产生[14]。

小波变换需要指定小波基函数，本文采用Haar小波基。Haar函数是小波分析中最早用到的一个具有紧支撑正交小波函数。经过小波变换后，图像被分解为4个维度的子图像，由1个低频信息、3个高频信息组成，图像分辨率降为原来图像的一半。首先对每一个数据源小波变换的低频信息按权重进行加权处理，权重值根据式(7)计算得到，经过加权处理后的信息直接作为小波反变换的低频信息参与融合。接下来对高频信息进行Kalman滤波处理。

首先对n(n为传感器数量)个数据源的高频信息进行Kalman滤波计算，生成3×n对高频信息；然后按照取最大值的原则，对n组高频信息进行筛选，最终由一个加权运算得到的低频信息和3个新的高频信息进行小波反变换，实现多源图像加权小波融合过程，算法对每幅图像都作上述处理；最后得到相应的8 d融合产品，再根据8 d融合产品进行月季年产品的合成。

3 结果与讨论

3.1 全球叶绿素α浓度融合产品的分析

全球叶绿素α浓度融合产品可从数据均值、方差、信息量3个方面进行分析评价。因融合产品时间跨度大，为能充分反映不同时间段内融合数据的可利用率，本文将数据分3组进行对比分析：第1组数据的时间段为1998—2002年；第2组数据的时间段为2005—2009年；第3组数据的时间段为2012—2016年。

1)均值分析。均值可以反映整体变化情况。通过分析全球叶绿素α融合产品均值(图1)，可以发现融合产品的均值变化有一定的规律。每年的图像都是46幅，图中刚好有5个波峰，说明一年的数据中叶绿素值总会达到峰值，峰值出现的时间在一年的中间，在年初和年末的产品中均值较低。3组融合产品的均值曲线变化规律相似，峰谷时间大致吻合，说明融合产品在时间尺度上有很大的一致性和可靠性。

图1 多年融合产品均值对比

2)方差分析。通过分析全球叶绿素α融合产品方差(图2)，可知融合产品的方差曲线拟合性很好，曲线变化趋势一致，说明融合产品质量稳定，没有缺失数据存在，其中方差最小的是第1组数据，说明该组数据稳定性最好。

图2 多年融合产品方差对比

3)信息量分析。信息量反映的是图像本身所含信息的丰富程度。对融合产品信息量进行分析(图3)，3组数据中信息量的曲线变化趋势一致，说明在全年或者长时间的数据中，信息量的变化是一致的，在特定的时间信息量会增大，而在其他时间就会减少。信息量最大的是第3组数据，其次是第2组数据，最后是第1组数据。

图3 多年融合产品信息量对比

3.2 全球叶绿素α浓度融合产品的验证

1)融合产品与实测值的对比分析。从实测数据集中选择每8 d最大值和对应的坐标，从融合产品中对应选择相同坐标和相同时间的叶绿素值，将2组数据绘制成图(图4)。从图4可知，二者的点匹配度较高，超过60%的点拟合度很好。

图4 融合产品与实测值相关性

2)融合产品与已有产品的对比分析。

①可利用率对比分析。将2008年的融合产品与欧空局对应时间的GSM叶绿素产品进行可利用率分析，如图5所示。可知，融合产品的可利用率高于GSM产品，以编号9产品为例，可利用率差为37.46%，其中22号产品的可利用率差最小，为28.20%，故本文的融合产品可利用率较已有产品有很大提高，为数据的更深层次应用提供了充分的保障。

图5 融合产品与GSM产品可利用率对比

②精度对比分析。将本文融合产品与2008年的实测值、欧空局的GSM产品进行对比分析，如图6所示。可知，本文融合产品、GSM产品与实测值拟合性均较好，但有部分点叶绿素浓度值相差较大，融合产品与GSM产品的叶绿素浓度值均高于实测值。

图6 融合产品和GSM产品、实测值对比

选取2008年的融合产品、GSM产品与实测最大值对比，共匹配了24对点，分别以实测-融合、实测-GSM作为横纵坐标值画出相关性图，如图7所示。根据数据的相关性，可知融合产品与实测值的相关性为0.792 2，GSM与实测值的相关性为0.349 4，证明本文融合产品质量较高，与实测点匹配较好。

图7 融合产品、GSM产品与实测值相关性对比

4 结束语

由于在轨运行的海色卫星大多是极轨卫星，叶绿素α浓度产品受到空间覆盖率和空间分辨率不同程度的影响，从而限制了卫星产品的实际使用。本文提出了基于小波变换和Kalman滤波技术的融合算法，将SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5个传感器的叶绿素α浓度数据进行融合，完成了融合产品的均值、方差和信息量的分析，并完成了融合产品与实测值、欧空局的GSM产品的对比分析。通过对比分析可知，本文的融合产品在空间覆盖率、数据可利用率、叶绿素α浓度值的精度方面均优于欧空局的GSM产品。本文构建的小波变换方法与Kalman滤波技术相结合的融合算法有效提高了融合产品的质量，在未来海色数据的应用方面具有较大的潜力。

本文中因有部分图像数据存在缺值现象，并没有在融合之前进行补值处理，因此影响了融合产品的质量，对融合产品的覆盖率和精度有一定的影响。另外，融合产品只与欧空局的GSM产品进行了对比，没有在实际项目中应用，对实际工作的贡献有待后续考察。