Landsat长时间序列数据格式统一与反射率转换方法实现

沈文娟， 李明诗

(南京林业大学林学院，南京 210037)

0 引言

Landsat卫星已积累了超过40 a的连续对地观测档案数据，这些数据在区域生态环境监测方面的应用优势明显[1]。已有研究将其用于重构过去几十年的森林干扰及恢复历史，并评价现存森林经营方法的有效性[2-3]。陆地卫星生态系统干扰自适应处理系统(landsat ecosystem disturbance adaptive processing system，LEDAPS)是美国航空航天局(NASA)的一个项目，主要研究长时间跨度的稠密时间序列Landsat GeoCover影像数据集的预处理算法。该影像数据集获取的真实地表反射率能够用于多种对象的长时间尺度评价。鉴于大区域范围土地覆盖变化的复杂性(包括火、森林采伐、城市化及农业领域变化等)，干扰模式及其与气候、生物多样性和生态系统健康的定量化关系，实现长时间序列遥感数据的预处理和分析是科学分析生态系统变化的前提条件[4-5]。

现有用于大气纠正的6S模型[6]，对云层较厚区域的反演精度较差[7]，应用时往往假设整个研究区的大气条件均一，太阳和传感器的观测角度一致，地表具有相同的海拔高度[8]，但是对反演区的气候和气象资料以及地面状况的了解信息不足，很多影像获取时刻的气象数据将很难准确获得，所以常使用一些标准模式来估算相关参数，特别是水汽、O3含量[9]。本研究以ETM+和TM数据为例，实现了LEDAPS算法可处理数据的格式统一，并验证长时间序列地面反射率处理产品的精度。在国内推广使用LEDAPS算法形成长时间序列可用的Landsat科学产品，对于推动多学科深入挖掘Landsat数据的科学价值具有重要的示范作用。

1 数据预处理

1.1 LEDAPS算法的原理

1.1.1 数据选择与定标

LEDAPS的定标方法首先使用日平均定标灯亮度历史数据移除先前错误的定标，再将修改后的定标系数用于相应日期的影像[4]。主要基于太阳方位、日地距离、TM或ETM+波段以及太阳辐照度，定标影像为表观(top-of-otmosphere,TOA)反射率。 Landsat TOA反射产品的输出格式为HDF-EOS，LEDAPS网站提供数据的下载。

1.1.2 大气校正原理

大气校正旨在补偿大气成分的散射和吸收效应，获取精确的地表反射率值。表观反射率ρTOA与地表反射率ρs的关系为

(1)

式中：Tg为其后面括号内所含气体的传输值；TR+A是瑞利和气溶胶传输值；ρR+A是瑞利和气溶胶大气固有的反射值；SR+A是瑞利和气溶胶球面反照率。传输、固有反射值和球面反照度由6S辐射传输模型计算所得[10]； O3含量来源于搭载在Nimbus-7，Meteor-3以及地球探测平台上的O3总量绘图系统(total ozone mapping spectrometer，TOMS)； 水蒸气来源于NOAA国家环境预报中心(national centers for environmental prediction，NCEP)重分析数据，分辨率为2.5°×2.5°(http: //rda.ucar.edu/datasets/ds090.0/)。数字地形和NCEP表面压力数据用于调整瑞利散射为实地的状况。

Kaufman等[11]充分利用影像中浓密植被(DDV)信息直接提取气溶胶光学厚度(aerosol optical thicrnecs,AOT)，在暗目标之间使用样条算法空间插值获取每个波段的气溶胶光学厚度值[10,12-13]。将插值得到的AOT和通过相关资料得到的O3、大气压以及水汽值用于6S传输模型算法，获取地表反射率。在处理过程中，还使用了Landsat7 自动云量评估算法进行云掩模[14]。在Landsat7项目中该方法提供以景为单位的云覆盖元数据具有优势[4]。

1.2 地表反射率及其验证

Masek等[4]总结了LEDAPS产品质量评价的3种方法： ①将LEDAPS得到的地面反射数据与Terra MODIS日反射产品作比较； ②Landsat 影像气溶胶光学厚度值与AERONET(aerosol robotic network)AOT作比较； ③Landsat地面反射值与高分辨率影像或者机载辐射数据的替代值作比较。本研究尝试获取LEDAPS处理前后的Landsat影像上典型地物反射率，并与标准波谱作比较，以此来判断处理结果的精度。

1.3 LEDAPS模型实现

LEDAPS预处理代码作为自动化的预处理和分析工具，是Landsat 1T数据定标、TOA反射率转换、云掩模和大气校正预处理的单独版本，并且该软件适用的数据为Landsat5和Landsat7标准EROS L1T产品，结果储存在HDF-EOS格式中。HDF-EOS数据图层包括所有反射波段的图层(不含热红外波段)、蓝光波段的气溶胶光学厚度图层以及一个“QA”图层(含有云的分离标记，缺失数据和来源于Landsat4的数据)。

LEDAPS模型(ftp://hrsl.arsusda.gov/pub/fgao/ Ledaps/)(图1)含4个具体操作命令，分别为lndpm(提取标准Landsat产品元数据信息，生成后续处理使用文件)、lndcal(进行TM和ETM+数据的定标，获取大气上界反射率并得到大气上界亮温)、lndcsm(使用自动云量评估算法得到云掩模)、lndsr(使用MODIS/6S辐射传输模型和黑暗密集植被算法所得气溶胶进行Landsat 反射波段的大气校正)。模型输出文件包含lndcal.*.hdf，lndth.*.hdf，lndcsm.*.hdf及lndsr.*.hdf，每个文件均带有地理信息的ENVI头文件。

图1 LEDAPS模型运行界面

2 数据格式的统一

在国内推广LEDAPS方法，可用的数据总结为以下3部分，即Landsat CDR(climate data records)产品、标准数据的LEDAPS处理数据以及非标准数据转换后的LEDAPS处理数据。综合以上数据，能够形成完整的长时间序列数据产品。

2.1 标准数据准备

2.1.1 LEDAPS处理后的地表反射率数据

Landsat地表反射率CDR数据由美国地质调查局(USGS)地球资源观测和科学(earth resources observation and science，EROS)中心科学处理部门(ESPA)处理所得。ESPA提供了地表反射率数据以及相关的源数据和生成数据。数据由LEDAPS软件处理所得。具体信息可参考http: //landsat.usgs.gov/documents/cdr_product_guide.pdf[15]。

2.1.2 LEDAPS可用的标准数据

标准数据的定义来源于USGS EDC或GeoCover的Landsat TM和ETM+数据，由LPGS或者NLAPS系统创建，具备有限的元数据格式(*_MTL.txt/*.met/*.H1/*.hdr)，其他数据格式或LPGS和NLAPS的其他类型均不能使用。含元数据(*_MTL.txt)格式的TM和ETM+数据部分可下载，数据包中含有Landsat单波段数据、GCP文件信息及元数据文件信息等。在LEDAPS处理中，主要使用数据包中的3种数据类型，其中元数据主要用于自动生成后续定标、掩模及大气校正所需的信息。

2.2 非标准数据格式处理

LEDAPS只能运行标准数据。在国内获取到的非标准Landsat数据(相对于标准数据格式)需要统一成标准数据格式后才能运行LEDAPS，得到地表反射率系列产品。FAST B数据包含band1—7.dat和header.dat以及产品信息文件共9个文件。在数据格式统一之前需先将国内地面站获取的*.dat数据与标准数据进行配准和重采样，转换后投影为UTM，基准面为WGS-84，波段5和波段7的分辨率为30 m，波段6的分辨率为60 m。然后将单波段数据转换成单波段GeoTiff格式，再模仿生成元数据信息。元数据信息主要参考头文件和产品信息文件，影像对应日期、太阳方位角和高度角、定标参数、角点经纬度坐标和投影坐标(分别在ENVI和ERDAS中读取)等需设置，控制点文件信息可以省略。数据格式统一后，后续产品实现参照1.3节。

3 方法试验与效果分析

LEDAPS在国内的成功使用案例较少，可参照的相关信息不多。中国科学院净月潭遥感实验站有关于旧版本(ledapsSrc.20110204.tar.gz)源数据程序编译安装的信息。本文推广标准数据和非标准数据的方法试验用于Landsat影像p120r38不同年份的数据转换和格式统一，完成了长时间序列(1987—2011)数据预处理及校正前后的评价。表1所列的影像信息除了包含多个时相的标准数据和非标准数据，还包含多个时相的Landsat CDR数据，含地表反射率数据，可供直接使用。

表1 长时间序列影像信息(1987—2011)

3.1 方法试验

3.1.1 标准数据转换方法

实现本试验采用长时间序列影像数据中成像日期为2003-01-28的ETM+数据，数据包中内容(参照2.1.2)有 L71120038_03820030128_B10-B80.TIF.L71120038_03820030128_GCP.txt及L71120038_03820030128_MTL.txt。经过LEDAPS程序转换得到的数据文件(参照1.3)有cloud_mask.log，L71120038_03820030128.carbon_met.txt,lndcal.L71120038_03820030128.hdf，lndcal.L71120038_03820030128.hdf.hdr，lndcal.L71120038_03820030128.txt，lndcsm.L71120038_03820030128.hdf，lndcsm.L71120038_03820030128.hdf.hdr，lndcsm.L71120038_03820030128.txt，lndsr.L71120038_03820030128.hdf，lndsr.L71120038_03820030128.hdf.hdr，lndsr.L71120038_03820030128.txt，lndth.L71120038_03820030128.hdf及lndth.L71120038_03820030128.hdf.hdr。

3.1.2 非标准数据转换和格式统一方法

首先将国内地面站获取的*.dat数据转换成标准数据格式，之后采取LEDAPS处理。本试验采用长时间序列影像数据中成像日期为2011-05-18的TM数据，数据包中内容(参照2.2.2)有band1.dat，band2.dat，band3.dat，band4.dat，band5.dat，band6.dat，band7.dat，header.dat，及ProductDescription.self。配准和重采样数据后得到新的数据包，投影为UTM，基准面为WGS-84。后续的处理参照3.1.1。

3.2 效果分析

LEDPAS处理的长时间序列地表反射率影像，其校正前后的影像视觉效果和反射辐射亮度有明显变化，校正后的影像更加清楚。

限于篇幅，本文在标准数据和非标准数据中各选一个年份的影像LEDDAS校正和评价(图2,3)。

图2 LEDAPS校正前(左)后(右) ETM+3,2,1真彩色合成图像(2003-01-28)

图3 LEDAPS校正前(左)后(右) TM3,2,1真彩色合成图像(2011-05-18)

从图2(左)可以看出，校正前表观反射率图像上地物边缘模糊，整幅图像总体呈偏蓝色调，这种辐射失真是由于在可见光波段，大气窗口内的太阳辐射受大气分子吸收作用的影响较少，蓝光波段受大气瑞利散射影响较大造成的。校正后影像(图2(右))气溶胶基本被消除，地物轮廓清晰，地物可识别性高，图像清晰度明显提高，基本还原了地物真实反射率，因而图像的亮度明显提高，无论是视觉效果还是所反映物理信息的准确性都得到较大改善。

从图3(左)可以看出，校正前表观反射率影像整体颜色偏暗，色彩不分明，清晰度不高； 校正后影像(图3(右))清晰度增加，能够识别出的地物种类增多，层次更加丰富，影像质量、视觉效果均得到较大改善，反演的物理信息更接近实况。特别是校正后地物边缘变得清晰，影像对比度提高，其中植被的颜色更鲜明，可识别性更高。

将图2，3中的局部区域(红框)放大，能够更好地比较校正前后的影像效果。

(a) 原始影像(b) TOA反射率影像(c) 地表反射率影像

(a) 原始影像(b) TOA反射率影像(c) 地表反射率影像

为了更好地理解本研究中标准数据和非标准数据经LEDAPS校正前后各波段反射率的变化，在LEDAPS处理得到的TOA和地面反射率影像上选择典型地物100个，将反射率取平均后与标准波谱曲线比较。图4和图5局部影像图分别对应选取的典型地物水体和植被。由图6(a)可知，LEDAPS处理后水体的波谱曲线形态与标准波谱曲线的相似，不同波段对应不同图像波谱的具体趋势比较与姚薇等[16]研究结果相似。由图6(b)可知，经过LEDAPS处理后的植被波谱曲线显示了典型绿色植被波谱特征，而未处理之前的植被波谱曲线特征不明显。此外，除处理后的近红外波段的反射率值低于标准值外，其他各波段的反射率值与标准值基本相同，与姚薇等[16]研究结果基本相同。分析校正前后各波段的动态变换范围，发现前3个波段的反射率值在校正后显著减小，而后3个波段的反射率值在大气校正后有所增大，这与大气瑞利散射和气溶胶的散射作用以及水汽的吸收密切关联，与杨静学等[9]的研究结果基本吻合。说明标准数据和非标准数据经过LEDAPS处理后得到的地面反射率产品均能有效地降低大气中臭氧、水汽及气溶胶等对影像波段反射率的影响，从而更有助于地表真实光谱信息的提取与研究，为遥感分类及建模应用研究提供更精确的数据源。

(a) 水体(2003-01-28)(b) 植被(2011-05-18)

分别对比图3和图2以及图5和图4可以看出，非标准数据的图像处理效果不及标准数据的处理效果。一方面，由于LEDAPS预处理方法的使用初期仅处理标准数据，本研究为了实现长时间序列反射率图像统一转换，将国内地面站获取的非标准数据统一至标准数据的格式，数据统一过程中不可避免存在误差； 另一方面，被处理的Landsat数据来源于国内外不同的生产机构，其初级产品的处理级别有差异。虽存有以上不足，考虑到Landsat连续时相数据成功应用的案例，以及可用标准数据的有限性，为了满足后续连续时相遥感影像相关的应用研究工作，现有的处理结果具有使用价值。并且在实现长时间序列的影像转换中TM和ETM+影像产品对于真实地物的光谱反应符合波谱规律，表明两者具备可比性。

4 结论

本研究介绍了一种长时间序列Landsat影像预处理程序LEDAPS，详细描述了该程序的数据选择、格式统一以及算法的实现过程，并在试验中成功地运用LEDAPS实现了不同来源的Landsat影像数据格式的处理，为土地覆盖变化和干扰等长时间序列的监测、生物物理参数遥感反演提供了科学产品。本研究将有助于促进国内形成处理长时间序列影像数据的准则。主要结论如下：

1)LEDAPS反射率转换方法提供了较为全面的参数用于6S辐射传输模型生成地表反射率数据。

2)基于LEDAPS 的Landsat影像校正产品的生成比单一的6S模型校正结果使用更加完备的参数，结果精度更有说服力。前者能自动化生成长时间序列产品，后者模型方法不断成熟并改进，为前者在方法的自动化、模块化、使用化及工程化的发展提供详尽的基础。

3)结果评价指出标准和非标准Landsat影像数据经LEDAPS处理前后反射波谱与标准波谱具有相似的形态，表明处理后的波谱更接近真实波谱，能有效地降低大气中O3、水汽及气溶胶等对影像波段反射率的影响。

Landsat影像以其能提供长时间序列产品的优势，为连续时相动态监测提供了可能性。但该方法在非标准数据参数设置、算法验证以及结果精度等方面尚有不足之处，是今后研究的重点。

参考文献(References)：

[1] Cohen W B,Goward S N.Landsat’s role in ecological applications of remote sensing[J].Bioscience,2004,54(6):535-545.

[2] Kennedy R E,Cohen W B,Schroeder T A.Trajectory-based change detection for automated characterization of forest disturbance dynamics[J].Remote Sensing of Environment,2007,110(3):370-386.

[3] Huang C Q,Goward S N,Masek J G,et al.Development of time series stacks of Landsat images for reconstructing forest disturbance history[J].International Journal of Digital Earth,2009,2(3):195-218.

[4] Masek J G,Vermote E F,Saleous N,et al.A Landsat surface reflectance dataset for North America,1990-2000[J].Geoscience and Remote Sensing Letters,2006,3(1):68-72.

[5] Wolfe R,Masek J,Saleous N,et al.LEDAPS:Mapping North American disturbance from the Landsat record[C]//Proceedings of the 2004 IEEE International Geoscience and Remote Sensing Symposium.IEEE,2004.

[6] 亓雪勇,田庆久.光学遥感大气校正研究进展[J].国土资源遥感,2005，17(4):1-6.

Qi X Y,Tian Q J.The advances in the study of atmospheric correction for optical remote sensing[J].Remote Sensing for Land and Resources,2005，17(4):1-6.

[7] 赵春江,宋晓宇,王纪华,等.基于6S模型的遥感影像逐像元大气纠正算法[J].光学技术,2007,33(1):11-15.

Zhao C J,Song X Y,Wang J H,et al.An algorithm based on 6S model for removing atmospheric effects form satellite imagery pixel by pixel[J].Optical Technique,2007,33(1):11-15.

[8] 徐永明,覃志豪,陈爱军.基于查找表的MODIS逐像元大气校正方法研究[J].武汉大学学报:信息科学版,2010,35(8):959-962.

Xu Y M,Qin Z H,Chen A J.A pixel-by-pixel atmospheric correction algorithm for MODIS data based on look-up table[J].Geomatics and Information Science of Wuhan University,2010,35(8):959-962.

[9] 杨静学,王云鹏,杨勇.基于高程或气溶胶厚度与6S模型校正参数回归方程的遥感图像大气校正模型[J].遥感技术与应用,2009,24(3):331-340.

Yang J X,Wang Y P,Yang Y.Atmospheric correction model of remote sensing image by regression analysis between elevation/aerosol optical thickness and correction parameters of 6S model[J].Remote Sensing Technology and Application,2009,24(3):331-340.

[10]Vermote E F,Saleous N E,Justice C O,et al.Atmospheric correction of visible to middle-infrared EOS-MODIS data over land surfaces:Background,operational algorithm,and validation[J].Journal of Geophysical Research:Atmospheres(1984-2012),1997,102(D14):17131-17141.

[11]Kaufman Y J,Wald A E,Remer L A,et al.The MODIS 2.1μm channel-correlation with visible reflectance for use in remote sensing of aerosol[J].IEEE Transactions on Geoscience and Remote Sensing,1997,35(5):1286-1298.

[12]Vermote E F,Saleous N,Justice C O.Atmospheric correction of MODIS data in the visible to middle infrared:First results[J].Remote Sensing of Environment,2002,83:97-111.

[13]Ouaidrari H,Vermote E F.Operational atmospheric correction of Landsat TM data[J].Remote Sensing of Environment,1999,70(1):4-15.

[14]Irish R R.Landsat7 automatic cloud cover assessment:Algorithms for multispectral,hyperspectral,and ultraspectral imagery[J].Proceedings of SPIE,2000,4049:348-355.

[15]Department of the Interior U S.Geological Survey.Landsat climate data record(CDR)surface reflectance product guide version 2.0[EB/OL].[2013-03-27].http://landsat.usgs.gov/documents.

[16]姚薇,李志军,姚珙,等.Landsat卫星遥感影像的大气校正方法研究[J].大气科学学报,2011,34(2):251-256.

Yao W,Li Z J,Yao G,et al.Atmospheric correction model for Landsat images[J].Transactions of Atmospheric Sciences,2011,34(2):251-256.