基于网格平台的 MODIS 自动化处理方法研究

2012-11-20张东映孟令奎夏辉宇成建国

水利信息化 2012年6期

张东映，孟令奎，夏辉宇，刁伟，成建国

（1.武汉大学遥感信息工程学院，湖北武汉 430079；2.水利部水利信息中心，北京 100053）

0 引言

近年来，遥感在水灾害监测、水资源管理、水土保持等水利业务中的应用卓有成效。随着水利工作的深入开展，急需为不断积累的大量卫星遥感信息资源提供统一的处理服务环境，提高数据加工处理与共享能力。基于以上需求，从提高卫星遥感数据处理的效率着手开展卫星遥感数据高效处理的研究工作，为水利部水利信息中心奠定良好的技术和数据资源基础。目前，水利部信息中心每日从国家气象中心接收的 MODIS L1B 级数据总量大约在80GB，该卫星每天覆盖我国境内至少2次，需要将实时的单景影像处理为植被指数、地表温度、假彩色合成等基础数据产品[1]。

基于以上需求，针对每日接收的 MODIS 数据开展数据处理流程及网格化平台建设的研究，以期为水灾害监测、水土保持、水环境监测等水利业务提供及时、可靠的基础数据支撑。

1 MODIS 数据自动化预处理技术研究

MODIS 原始数据在气象卫星中心已经完成了由0级到 L1B 级产品的处理，因此数据处理主要完成2部分内容。

1.1 MODIS L1B 级数据预处理

MODIS L1B 数据经过了辐射和部分几何校正，但是还没有去除扫描成像时产生的 Bow-Tie 效应现象[2]。此外，数据使用的正弦投影坐标并非通用的坐标系统，并且数据格式也还没有转换成通用的遥感数据处理格式。所以 MODIS 数据处理的主要目标为1）消除数据中的 Bow-Tie 效应；2）根据 HDF（层次型数据格式）文件包含的经纬度信息对消除了Bow-Tie 效应的 MODIS L1B 影像进行几何校正，并将影像纳入 WGS-84坐标体系中；3）将处理得到的数据转换为 GeoTIFF 格式；4）根据辐射定标参数，将记录DN值（影像的像元数值）的原始数据文件转换为反射率和热辐射率数据。

传统预处理采取手工处理方式，用 ENVI（遥感图像处理软件）和 ERDAS 等遥感软件实现，步骤繁杂重复，费时费力，且容易出错。而通过自动化的批处理模式可以实现大量数据的预处理工作，不需要人工干预。应用 IDL 语言对 ENVI 进行二次开发是有效而快速的解决方案。图1是 MODIS 数据预处理流程图。

图1 250与500m 分辨率数据融合流程图

MRTSwath（The MODIS Reprojection Tool Swath）用来处理 MODIS L1B 数据，是由 MODIS 数据陆地研究小组利用 Java 语言开发的一个简易软件平台。MRTSwath 软件的核心功能就是实现 MODIS 数据从Swath 到 Grid 的格式转变，从而完成 MODIS L1B 数据的预处理工作（包括去除 Bow-Tie 效应和投影转换等几何校正过程）。辐射定标结果包括表观反射率和热辐射产品。MODIS L1B 产品一共36个波段，包括2种信息：波段 1～9和26为地物反射信息；波段 20～25和 27～35为地物发射信息。用 ENVI 打开 MODIS 数据，就是反射率（大气外层表观反射率）、辐射亮度及发射率3个数据类型。

由于 MODIS 信号数据的精度很高，用浮点数据保存文件占用空间较大，因此，MODIS L1B 采用16-bit 整数和尺度转换的方法，将浮点数据通过偏移量（offset）和尺度因子（scale）这2个参数转换为16-bit 整数，以减小文件大小。这些无符号的 16-bit数据没有实际的物理意义，因此，必须通过尺度转换的方法将其整数型数据重新转换为反射率值（反射通道）和热辐射度（发射通道）[3]。辐射亮度计算公式为LB1T2FS=RadS1B（SIB1T2FS-Rad_OB）；表观反射率计算公式为RB1T2FS=RefS1B（SIB1T2PS-Ref_OB）。式中：SIB1T2PS为某波段某像素点的计数值，B为相应波段号；RadS1B为辐射能量缩放比；Rad_OB为辐射能量偏移量；RefS1B为反射率缩放比；Ref_OB反射率偏移量；LB1T2FS为某波段某像素点的辐射亮度；RB1T2FS为某波段某像素点的表观反射率。上述参数均可以在相应波段科学数据集的属性域中读取。

综上所述，预处理数据产品包括3部分：MODIS L1B 预处理、热辐射和反射等产品。热辐射和反射产品是由 MODIS L1B 预处理数据根据热辐射率和反射率公式生成的。具体的说，热辐射产品由1000m 的热辐射波段生成，反射产品包括 250，500和1000m 的反射波段数据。

1.2 MODIS 产品再加工处理

根据现有数据和灾害监测业务需求，MODIS 数据自动化处理程序将提供处理和展示部分水利信息应用的影像产品。具体地说，由 MODIS 预处理数据通过云检测生成3个层次的再处理产品，第1层次产品为单幅的辐射亮温、地表温度、NDVI和假彩色合成产品；第2层次产品为全国日拼接地表温度、全国日拼接NDVI和全国日拼接假彩色合成产品；第3层次产品为旬NDVI和地表温度合成产品，然后再进行归一化合成。图2为 MODIS 预处理产品再加工处理逻辑流程。

图2 250与500m 分辨率数据融合流程图

1）云检测生成云二值图。在遥感图象处理中，云覆盖是常遇到的图像噪声之一，不仅为图像的处理带来许多困难，而且使后续的图象识别、分类精度难以保证，有时甚至无法进行，所以在生成第2层次产品之前，要进行云检测去除云的影响。云在可见光波段具有高反射率，遥感图像能较好地区分陆地和云的边界；云在近红外波段的波谱特征主要与大气中的含水量有关，主要反映大气中的水汽特征，即为水汽吸收谷。由于云与各种地物波谱特征形成明显反差，因此将其归一化处理，不仅可以突出云的信息，而且可以部分消除太阳高度角、卫星扫描角及大气程辐射的影响，由此通过设置适当的阈值将很容易区分云和其他地物。但是归一化处理后云与水的值差不多[2]。

2）产品分类。a）假彩色合成产品。考虑到近红外波段范围内植被、土壤、水分反射率差异较大，用近红外波段代替绿波段进行假彩色合成，图像包含的信息更加丰富，即 MODIS 影像的 1，2，3波段组合。因为 1，2两个波段的分辨率都是250m，而3波段的分辨率是500m，为了保留影像的高分辨率，将波段3重采样成250m，合成 1，2，3波段为 R，G，B 通道的假彩色图。b）全国 MODIS 影像拼接产品。Terra 卫星是太阳同步极轨卫星，我国白天过境的数据加上晚间过境数据，就可以得到每天最少白天和夜晚各1次更新数据，如果结合 Aqua卫星，可获得4次/d 的数据。每天接收到的 MODIS数据基本可以覆盖全国区域。实时获取当天全国MODIS3种分辨率的多轨合成图，应用于洪涝灾害等应急监测处理，为水资源利用和水环境调查提供动态监测。通过对同一地区多幅少云清晰影像进行光谱融合处理后，获取某时期全国的 MODIS 多天合成图。c）地表温度产品。地表温度产品首先需要由MODIS 数据处理生成辐射亮温数据。用普朗克方程可以反演出某个波段影像的亮度温度。首先采用最邻近法将云二值图重采样成1000m 的分辨率；其次提取辐射产品中的热红外波段（20～36波段，分辨率1000m），代入普朗克反演公式生成辐射亮温图。将辐射亮温数据与重采样后的云二值图做叠加处理，去除云的干扰；最后去云后的辐射亮温图按每6h 做时段拼接，生成4幅/d 亮温拼接图（上午、下午、前半夜和后半夜）。数据同时入库和做文件存储。地表温度反演在辐射亮温产品基础上，采用分裂窗算法进行反演。适用于 MODIS 数据的地表温度反演的分裂窗算法计算[3]。MODIS 数据的地表温度反演的分裂窗算法公式为Ts=A0+A1T31-A2T32。

式中：Ts是地表温度；T31和T32分别是 MODIS 第 31和32波段的辐射亮度温度；A0，A1，A2是分裂窗算法的参数。

电台发布的气象预警信息主要是由气象部门从广电中心直接推出的，其主要是以声音为传播媒介，电台的优势是能定时发布天气预报，还能随时插播预警信息，方便广大人民群众及时了解气象灾害的程度及预防措施。目前的收音机虽然体积较小、便于携带，但缺点同电视节目一样，需要打开收音机才能获知信息，并且电台节目的受众人群比较有限，缺乏视觉感染力，所以对农村居民的影响比较有限。

3）归一化植被指数合成。如果一个波段在可见光区域（VIS），另一个波段在近红外（NIR），则NDVI计算公式为NDVI= （NIR-VIS）/ （NIR+VIS）。

对于 MODIS 数据来说，计算NDVI所需要的可见光和近红外波段分别对应 MODIS 的第1和2波段。NDVI的合成包括日和旬合成，都遵循 MVC 原则，即最大值合成方法。对同一地理位置坐标点，选取最大的NDVI值作为合成最终结果。

2 遥感数据网格化并行处理平台建设研究

遥感数据网格化并行处理平台建设主要是基于网格调度软件 Platform 对遥感数据的预处理模块进行拆分，实现遥感数据网格平台下的并行计算处理，从而提高海量遥感数据的处理速度，构建高效的遥感数据网格化并行处理平台[4-5]。

1）基于 Platform 的网格平台搭建。Platform 作业调度解决方案主要包括 Platform Process Manager和 Platform LSF2个模块。Platform Process Manager能提供卫星数据处理流程的定义、执行、实时监视和控制。Process Manager 的服务将实际作业提交到LSF 管理的机群中[6]。

Platform 的软件环境与网络结构如图3所示。

图3 Platform 的软件环境与网络结构

其中的计算节点服务器配置如表1所示。

2）MODIS 数据处理的流程分解。面向水利应用的 MODIS 数据处理的整个流程比较复杂，需要经过几何校正、辐射校正、植被指数计算、地表温度计算等过程。根据 Platform 特点，将原有 MODIS 预处理流程以相同粒度划分为不同处理模块，如图4所示。

表1 网格节点服务器配置

图4 预处理流程任务分解图

模块 PreRef1km，PreRad1km代表1000m 分辨率数据几何校正，其处理流程如图5-a 所示。在此模块中，将反射率和发射率波段分开来处理，便于后续处理，且节省了处理时间。

模块 Prehkm和 Preqkm分别是500m 分辨率数据几何校正和250m 分辨率数据几何校正，处理流程和1000m 类似，处理流程如图5-b 所示。模块Ref1km和 Rad1km是1000m 分辨率数据辐射校正，其中Ref1km对应反射率波段，Rad1km对应发射率波段。

模块 Refhkm，Refqkm别是500m 分辨率数据辐射校正和250m 分辨率数据辐射校正，处理流程和Ref1km类似。模块 NDVI1km，NDVIqkm是1000和250m 分辨率的NDVI计算，其处理流程相同，流程如图5-c 所示。

模块 RBT 是辐射亮温处理，处理流程如图5-d所示。

模块 LST 是地表温度计算，处理流程如图5-e所示。

模块 Merge 是250与500m 分辨率数据融合，处理流程如图5-f 所示。

2）流程重组。对于 1000，500和250m 这3种分辨率影像的处理可以相互独立到不同节点中并行处理。只在 Merge 模块中需要同时用到500和250m影像处理的产品。另外，为了减少节点之间的数据传输，将 Ref1km和 NDVI1km处理节点设置为与PreRef1km相同，将 Rad1km，RBT1km，LST1km处理节点设置为与 PreRad1km相同，将 Refhkm，Merge 处理节点设置为与 Prehkm相同，将 Refqkm，NDVIqkm处理节点设置为与 Preqkm相同，流程重组图如图6所示。

3）网格服务封装与调用。所有模块都封装为带参数的 .exe 可执行程序，组成 MODIS 数据处理模块库。该模块库部署于每台执行节点中，且程序、输入输出路径都保持一致。利用 Process Manager 提供的可视化流编辑工具 Flow Editor，编辑可在 LSF 上执行的工作流。Flow Editor 定义的工作流为 XML格式，PlatForm LSF 调度程序读取该 XML 文件，获取传递任务信息及执行任务所需系统资源配置信息，并根据请求提交相应的处理任务，PlatForm 流程编辑图发图7所示。