海水表面温度是一种重要的海洋环境参数, 被广泛应用于海洋动力过程、海气相互作用和全球气候变化等研究领域中. 但就传统原位测量方式而言, 其设备成本代价高、站点分布离散、监测范围有限, 无法满足较大区域范围、较长时间序列和连续性观测的需求. 针对以上不足, McMillin^[1]最早提出使用遥感影像劈窗算法反演地表温度. 随着遥感技术和卫星反演算法的成熟, 覃志豪等^[2]提出采用MODIS遥感影像反演海表温度产品, 并应用于相关大尺度、长周期的研究中. 同时伴随遥感数据的海量增长和应用算法复杂度的急剧增加, 传统计算方式难以满足日常的影像计算需求. 朱义明、付天新等^[3,4]使用MapReduce编程模型对遥感影像实现了并行处理. Almeer、常生鹏等^[5,6]使用Hadoop实现了图像处理中的滤波、锐化、增强和分割等算法. 刘震^[7]提出了一种基于五元组的遥感图像并行处理方法. 以上遥感数据并行计算研究均采用离线存储的方式, 数据下载、人工预处理程度高; 算法设计的参数变量较单一, 对于中间输入参数较多的劈窗算法并没有提供相应的解决方案参考; 没有充分利用RDD窄依赖特性提高计算的效率. 本文综合相关研究和算法特性, 提出了一种基于Apache Spark的MODIS海表温度快速反演方法, 融合在线获取、实时解析、弹性分布式数据集处理和适于快速计算的反演算法提高海量遥感数据处理效率.

1 方法

本文研究涉及的主要方法有MODIS影像在线获取与动态解析、适于快速计算的改进劈窗算法和基于Apache Spark的遥感影像计算过程优化. MODIS影像在线获取与动态解析使用自动化的影像查询机制提高影像定制的效率; 适于快速计算的改进劈窗算法通过减少输入参数, 降低模型复杂度的形式提高算法的时空效率; 基于Apache Spark的遥感影像计算过程优化通过RDD窄依赖关系减少计算过程中的网络数据交换延迟.

1.1 MODIS影像在线获取与动态解析

针对人工影像数据获取的低效性, 本文采用自动化影像条件查询与多线程下载程序实现海量数据自定义批量下载, 避免了中间过程的人工干预, 提高了影像获取的效率. 美国宇航局以RESTful API方式提供了查询MODIS影像元数据和服务器存储路径的网络请求方法. 如图1, 用户可以通过不同的网页请求参数获取不同的影像元数据信息. 如通过searchForFiles页面传递产品名称、数据集编号、时空范围和影像状态信息等参数查询研究区影像ID列表; 通过getFileProperties页面传递影像ID列表查询具体影像的状态信息, 并据此下载影像或提交订单; 通过getFileUrls页面获取影像的服务器路径, 并据此下载相应的影像数据文件. 本文为了提高处理效率, 避免磁盘IO, 在数据处理前将影像文件以字节数组的形式缓存在内存中.

针对影像数据提取与解析过程, 本文采用一种适于并行化处理的第三方处理方案. MODIS L1B数据文件存储格式是一种自描述多对象超文本文件格式——HDF4. HDF4通过数据描述符和数据元素来管理不同的数据对象. 数据描述符包含了标识符、参照数、数据偏移长度和数据长度. 通过相应的数据描述符可以将文件中不同对象的数据提取出来, 达到解析的目的. Spark兼容所有Hadoop支持的文件格式, 但是对于如HDF4这样的遥感数据并不提供原生的支持, 需要通过自定义的方式进行解析. 目前常用的第三方解析方案有美国国家超级计算机应用中心基于C语言的解析方案和美国大气研究中心基于Java的解析方案. 本文考虑到Spark运行于Java虚拟机支持的开发环境和平台可移植性的需求, 采用基于Java的解决方案. 通过自定义的解析程序, 可以将Spark不识别的HDF4文件数据, 提取转化为其支持并行化的操作数据.

1.2 适于快速计算的改进劈窗算法

劈窗算法, 是一种利用相邻两个热红外通道对大气水汽吸收作用的差异特性, 通过波段组合消除大气影响, 最终将海表温度表达为波段亮度温度函数的方法. 本文利用影像元数据、影像值和大气参数拟合公式, 减少反演过程中参数输入和拟合的效率损耗, 避免计算节点间数据交换. 其具体算法流程如图2.

1.2.1 辐射定标

MODIS L1B数据文件中保存的波段像元值为16位量化级数. 根据MODIS用户手册^[8]所述, 反演计算中所需的辐射亮度值需要通过以下公式进行计算:

${{Radiance = Scale}} \times \left( {{{DN}} - {{Offsets}}} \right)$

(1)

其中, Radiance为辐射校正后的波段辐射亮度, Scale和Offsets为MODIS元数据提供的定标参数, DN为波段量化级数.

1.2.2 大气水分含量

大气透过率通常是根据与大气水分含量的关系计算而得. Kaufman和Gao^[9]提出利用MODIS相应光谱范围的波段函数组合获取大气水分含量的方法, 相较地面观测数据, 其误差在5%~10%之间.

$\omega = {\left[ {\frac{{\alpha - {\rm{ln}}\left( {\frac{{{\rho _{19}}}}{{{\rho _2}}}} \right)}}{\beta }} \right]^2}$

(2)

其中, $\omega $ 为大气水分含量(g·cm^–2), $\alpha$ =0.02, $\beta = 0.632$ , ${\rho _i}$ 为波段 $i$ 的地面反射率.

1.2.3 大气透过率

大气透过率与大气水分含量的关系一般通过大气传输模型模拟确定. 毛克彪等^[10]针对大气传输模型参数难以实时获取的问题, 利用神经网络拟合MODTRAN模型, 获得了了以下统计回归表达式, 其波段反射率平均误差不大于0.008.

${\tau _{31}} = 2.9 - 1.88 \times {e^{\frac{\omega }{{21.23}}}}$

(3)

${\tau _{32}} = 4.6 \times {e^{- \textstyle\frac{\omega }{{ 32.71}}}} - 3.59$

(4)

其中, ${\tau _i}$ 为 ${i}$ 波段的大气透过率, $\omega $ 为上一步所求大气水分含量. 考虑到反演计算过程中外部模型参数计算与输入会引起节点间数据交换的效率损耗, 本文采用以上统计回归表达式计算大气透过率.

1.2.4 线性组合系数^[11]

地表比辐射率主要取决于地表的物质组成, 就海表温度反演而言, 其地表类型单一, 因此直接取相应的水体比辐射率.

${\varepsilon _{31 {\rm{水体}}}} = 0.996\;83,{\rm{ }}{\varepsilon _{32 {\rm{水体}}}} = 0.992\;324$

(5)

${C_i} = {\varepsilon _i}{\tau _i}\left( \theta \right)$

(6)

${D_i} = \left[ {1 - {\tau _i}\left( \theta \right)} \right]\left[ {1 + \left( {1 - {\varepsilon _i}} \right){\tau _i}\left( \theta \right)} \right]$

(7)

其中, $i$ 为31或32波段, ${C_i}$ 、 ${D_i}$ 为计算线性组合系数 ${A_j}\left( {{j} = 1,2,3} \right)$ 的中间参数, ${\varepsilon _i}$ 是波段 $i$ 的地表比辐射率, ${\tau _i}\left( \theta \right)$ 为波段 $i$ 视角为 $\theta $ 的大气透过率.

$ \begin{array}{c} \displaystyle{A_0} = \left[ {\frac{{{D_{32}}\left( {1 - {C_{31}} - {D_{31}}} \right)}}{{{D_{32}}{C_{31}} - {D_{31}}{C_{32}}}}} \right]{a_{31}} - \left[ {\frac{{{D_{31}}\left( {1 - {C_{32}} - {D_{32}}} \right)}}{{{D_{32}}{C_{31}} - {D_{31}}{C_{32}}}}} \right]{a_{32}}\displaystyle \\ {A_1} = 1 + \frac{{{D_{31}}}}{{{D_{32}}{C_{31}} - {D_{31}}{C_{32}}}} + \left[ {\frac{{{D_{32}}\left( {1 - {C_{31}} - {D_{31}}} \right)}}{{{D_{32}}{C_{31}} - {D_{31}}{C_{32}}}}} \right]{b_{31}}\displaystyle \\ {A_2} = \frac{{{D_{31}}}}{{{D_{32}}{C_{31}} - {D_{31}}{C_{32}}}} + \left[ {\frac{{{D_{31}}\left( {1 - {C_{32}} - {D_{32}}} \right)}}{{{D_{32}}{C_{31}} - {D_{31}}{C_{32}}}}} \right]{b_{32}} \end{array}$

其中, ${a_{31}}$ = –64.603 63 , ${a_{32}}$ = –68.725 75, ${b_{31}}$ = 0.440 817, ${b_{32}}$ = 0.473 453.

1.2.5 亮度温度

亮度温度是与物体在相同频率下具有相同辐射率的黑体的绝对温度, 是海表温度反演的重要参数之一. 利用普朗克黑体辐射方程可以直接求解相应波段的亮度温度.

${T_i} = \frac{{{k_{i,2}}}}{{{\rm{ln}}\left( {1 + \frac{{{k_{i,1}}}}{{{B_i}}}} \right)}}$

(11)

其中, ${T_i}$ 为波段的亮度温度, $ {k_{31,1}} = 729.541\; 636$ W·m^–2·sr^–1·μm^–1, $ {k_{32,1}} = 474.684\;7799$ W·m^–2·sr^–1·μm^–1, $ {k_{31,2}} = 1304.413\;871$ K, $ {k_{32,2}} = 1196.978\;785$ K

1.2.6 海表温度反演

覃志豪等^[2]综合研究国内外17种劈窗算法, 通过比较其参数复杂性和算法精度, 在保证较高反演精度的情况下提出了线性反演公式. 其具有良好的反演精度, 已被应用于农业旱灾监测. 本文基于简化输入参数, 减少像元间计算引起的节点数据交换的需要, 采用以下线性反演公式进行海表温度的快速反演.

${T_s} = {A_0} + {A_1} \times {T_{31}} - {A_2} \times {T_{32}} - 273.15$

(12)

其中, ${T_s}$ 是海表温度(℃), ${T_{31}}$ 和 ${T_{32}}$ 分别是MODIS第 31 和 32 波段的亮度温度, $ {A_0}$ 、 ${A_1}$ 和 ${A_2}$ 是劈窗算法参数.

1.3 基于Apache Spark的遥感影像计算过程优化

针对基于Spark的海表温度反演, 本文综合考虑了算法时间复杂度和空间复杂度, 采用如图3所示计算流程.

网络下载的影像文件以独立影像的形式保存在分布式文件系统中(本文采用的是HDFS). 为减少影像处理时节点间的数据传输, 数据块的大小一般按照影像文件的大小分块. 如数据容量最大为169 MB, 分块大小应至少设置为169 MB才能保证不会因文件分块导致跨节点访问. Spark作业开始时, 影像文件以绝对路径集合的形式传入, 经并行化操作构建包含有n个分区的键值对弹性分布式数据集 (以下简称RDD). 其中, 键名为每个文件的绝对路径字符串, 键值为相应的影像文件. 利用第三方解析方案提取出每个影像文件中的波段数据. 针对每个影像的波段数据, 按照算法设计, 结合Spark提供的转换算子和动作算子按照劈窗算法执行相同的波段组合函数, 进行海表温度的反演. 为了避免数据节点间的数据块传输和RDD间的组合计算, 从创建键值对RDD-0到生成键值对RDD-2的过程中, 通过对影像文件分块优化、算法过程优化和计算过程优化实现RDD间的窄依赖关系. 最终的反演结果保存到分布式文件系统中.

2 实验 2.1 实验数据

本文实验数据采用MODIS L1B Terra星载数据(https://ladsweb.modaps.eosdis.nasa.gov/search/), 数据空间范围涵盖山东半岛、渤海湾一带, 时间范围为2015年1月1日至2015年12月25日, 共计6634景, 总计1000 GB. 由于本文采用的劈窗算法需要可见光波段参与计算, 因此实验数据均为白昼期间获取的影像.

2.2 硬件环境

整个集群构建在相对封闭的千兆局域网中. 集群采用主从式架构, 使用商用台式机作为节点硬件基础.

主节点负责整个集群的任务提交、分配、管理与汇总, 不参与数据计算过程. 从节点负责数据的分布式存储与实际计算. 集群及节点详细硬件配置如图4和表1.

2.3 软件环境

为了保持节点行为的一致性, 所有节点均采用相同的软件配置. 操作系统选用的是CentOS,在其上构建程序运行所需的Java和Scala运行环境. Hadoop框架提供计算平台所需的分布式文件系统HDFS和统一的资源调度机制YARN. 软件环境配置如表2.

3 结果及分析

实验结果按照是否使用Spark集群计算框架, 分为单机结果和集群结果. 单机结果在单一计算节点上本地运行获取, 集群结果在YARN Cluster模式下获取. 每种结果均采用不同的数据量作为变量进行测试: 10 GB (66景)、50 GB (329景)、100 GB (658景)、500 GB (3347景)和1000 GB (6634景). 每项数据均在相同实验条件下重复三次, 取平均值. 具体实验结果如表3.

在单机模式下, Scala程序能够实现约每秒0.04 GB (或每秒0.29个影像文件)的处理能力. 在加入Spark并行处理框架后, Scala程序能够实现约每秒0.40 GB (或每秒2.67个影像文件)的处理能力. 总体上, 在集群模式下影像处理的效率约为单机模式的10倍. 就不同的数据量而言, 当数据量较小时, 由于集群启动、并行化和任务调度等的开销, 处理效率低于平均水平. 当数据量在50~500 GB时, 影像处理效率呈现略微上升的态势. 但当数据量达到1000 GB时, 影像的处理效率出现小幅度回落. 可能在集群硬件资源达到饱和时, 在相同分区数量下, 单一分区的处理粒度加大, 处理任务间的间隙、资源调度与任务分发延迟加大, 导致了效率的下降.

4 结论

本文提出了一种将Spark大数据并行计算技术应用于海表温度反演的方法. 通过将自动化影像获取过程、简化算法流程和RDD窄依赖并行计算结合, 实现了一种海表温度快速反演的机制. 通过对比单机模式和集群模式的处理效率, 验证了基于Apache Spark的海表温度反演方法的快速计算能力. 通过将大数据并行计算技术应用于遥感图像处理, 可以有效提高数据产品的更新频率, 为大范围、长周期和高分辨率的海洋、气象观测提供良好的支持. 本文采用的改进简化劈窗算法和RDD窄依赖特性也可以应用于其他遥感产品生产和海量数据并行计算过程.