警告：本教程中介绍的用于协调 Landsat ETM+ 和 OLI 数据的程序已过时，在处理 Landsat Collection 2 地表反射率数据时，不建议使用这些程序，也没有必要使用。如需了解详情，请参阅有关 Landsat 协调的常见问题解答条目。Collection 1 数据已弃用，本教程中的代码将不再有效或更新。
本教程涉及将 Landsat ETM+ 地表反射率协调为 Landsat OLI 地表反射率。它提供：
- 一种频谱转换函数
- 用于创建随时可供分析的数据的函数
- 时序图表示例
本文旨在提供一份端到端指南，帮助您协调和直观呈现 35 年以上的 Landsat 区域时序数据，以便立即应用于您感兴趣的区域。
请注意，将 ETM+ 转换为 OLI 的系数也适用于 TM。因此，在本教程中，提及 ETM+ 时即指 TM。
Landsat 简介
Landsat 是一项卫星成像计划，自 1972 年以来一直在收集中等分辨率的地球图像。作为历时最长的太空地球观测计划，它为识别地貌变化的时空趋势提供了宝贵的时间记录。本教程中使用了 Thematic Mapper (TM)、Enhanced Thematic Mapper Plus (ETM+) 和 Operational Land Imager (OLI) 仪器数据。它们密切相关，并且相对容易整合为一致的时间序列，从而生成 1984 年至今的连续记录，每个传感器的频率为 16 天，空间分辨率为 30 米。多光谱扫描仪 (MSS) 仪器将 Landsat 记录的时间范围扩展到了 1972 年，但本教程中未使用该仪器。其数据差异很大，因此很难与后来的传感器集成。如需查看跨所有传感器的协调示例，请参阅 Savage 等人的著作（2018 年）和 Vogeler 等人的著作（2018 年）。
为何要进行协调
Roy 等人 (2016) 表明，根据应用的不同，Landsat ETM+ 和 OLI 的光谱特征之间存在细微但可能显著的差异。您可能需要协调数据集的原因包括：生成跨越 Landsat TM、ETM+ 和 OLI 的长时间序列；生成近期的年内合成数据，以减少 ETM+ SLC 关闭间隙和云/阴影掩码造成的观测缺失影响；或提高时间序列内的观测频率。如需了解详情，请参阅上文中的链接手稿。
操作说明
函数
以下是一系列将 ETM+ 协调到 OLI 并生成分析就绪数据所需的函数，这些数据将在本教程的应用部分中使用，以直观呈现像素的光谱时间历史记录。
协调
通过根据 Roy 等人 (2016) 表 2 中的 OLS 回归系数将 ETM+ 光谱空间线性转换为 OLI 光谱空间，实现了协调一致性。以下字典中定义了各个波段的系数，其中包含斜率 (slopes) 和截距 (itcps) 图像常量。请注意，y 轴截距值乘以了 10,000，以匹配 USGS Landsat 地表反射率数据的缩放比例。
var coefficients = {
itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172])
.multiply(10000),
slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
同一光谱响应窗口的 ETM+ 和 OLI 波段名称不相等，需要进行标准化。以下函数仅从每个数据集中选择反射率波段和 pixel_qa 波段，并根据它们所代表的波长范围对其进行重命名。
// Function to get and rename bands of interest from OLI.
function renameOli(img) {
return img.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
// Function to get and rename bands of interest from ETM+.
function renameEtm(img) {
return img.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
最后，定义转换函数，该函数会将线性模型应用于 ETM+ 数据，将数据类型转换为 Int16 以与 OLI 保持一致，并重新附加 pixel_qa 波段以供日后在云和阴影掩码中使用。
function etmToOli(img) {
return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
.multiply(coefficients.slopes)
.add(coefficients.itcps)
.round()
.toShort()
.addBands(img.select('pixel_qa'));
}
云和阴影遮盖
准备好进行分析的数据应已遮盖云和云阴影。以下函数使用 CFmask（Zhu 等人，2015 年）中包含的pixel_qa 波段，用于将识别为云和云阴影的像素设置为 null。
function fmask(img) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = img.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask)
.eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return img.updateMask(mask);
}
谱指数计算
即将推出的应用使用归一化燃烧比率 (NBR) 光谱指数来表示受野火影响的森林像素的光谱历史记录。之所以使用 NBR，是因为 Cohen 等人 (2018) 发现，在 13 个光谱指数/波段中，就美国各地的森林扰动（信号）而言，NBR 的信噪比最高。计算公式如下。
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
在您自己的应用中，您可以选择使用不同的光谱指数。您可以在此处更改或添加其他索引。
合并函数
为每个数据集定义一个封装函数，该函数可整合上述所有函数，以便方便地将其应用于各自的图片集合。
// Define function to prepare OLI images.
function prepOli(img) {
var orig = img;
img = renameOli(img);
img = fmask(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
// Define function to prepare ETM+ images.
function prepEtm(img) {
var orig = img;
img = renameEtm(img);
img = fmask(img);
img = etmToOli(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
在应用中，您可以决定包含或排除某些函数。根据需要更改这些函数。
时间序列示例
上述 prepOli 和 prepEtm 封装函数可以映射到 Landsat 地表反射率集合，以创建可用于分析的跨传感器数据，从而直观呈现像素或像素区域的光谱时间顺序。在此示例中，您将创建一个 35 年以上的时序，并显示单个像素的光谱历史记录。此特定像素与成熟的太平洋西北针叶林地块（图 1）的近期历史有关，该地块在 20 世纪 80 年代和 90 年代经历了一些扰动，并在 2011 年发生了一次高强度燃烧。
图 1. 位置和网站字符，例如感兴趣的区域。美国俄勒冈州胡德山北坡上成熟的太平洋西北针叶林。图片由以下机构提供：Google 地球、美国农业部林务局、Landsat 和 Copernicus。
定义感兴趣的区域 (AOI)
此应用的结果是像素的 Landsat 观测时间序列。ee.Geometry.Point 对象用于定义像素的位置。
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
在您的应用中，您可能会选择其他像素。为此，请替换上述经度和纬度坐标。或者，您也可以使用其他 ee.Geometry 对象定义（例如 ee.Geometry.Polygon() 和 ee.Geometry.Rectangle()）来总结一组像素的光谱历史记录。如需了解详情，请参阅开发者指南的几何图形部分。
检索 Landsat 传感器集合
获取 OLI、ETM+ 和 TM 的 Landsat USGS 地表反射率集合。
点击链接可详细了解每个数据集。
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
定义图片集过滤条件
定义一个过滤条件，以根据感兴趣区域的空间边界、光合作用高峰季节和质量来限制图片集。
var colFilter = ee.Filter.and(
ee.Filter.bounds(aoi), ee.Filter.calendarRange(182, 244, 'day_of_year'),
ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
ee.Filter.or(
ee.Filter.eq('IMAGE_QUALITY', 9),
ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));
您可以根据需要更改这些过滤条件。在上述链接的数据说明页面的“图片属性”标签页下，找到要按哪些其他属性进行过滤。
准备合集
合并集合，应用过滤条件，并对所有图片应用 prepImg 函数。以下代码段的结果是一个 ee.ImageCollection，其中包含来自 OLI、ETM+ 和 TM 传感器集合的符合过滤条件的图片，并且这些图片已处理为可用于分析的 NBR。
// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOli);
etmCol = etmCol.filter(colFilter).map(prepEtm);
tmCol = tmCol.filter(colFilter).map(prepEtm);
// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);
制作显示所有观测结果的时序图表
通过绘制所有观测结果的散点图（其中传感器通过颜色区分），可以定性评估传感器间的一致性。Earth Engine 的 ui.Chart.feature.groups 函数提供了此实用程序，但首先必须将观测到的 NBR 值作为属性添加到每个影像中。对图像集合应用可计算与 AOI 相交的所有像素的中位数并将结果设置为图像属性的 map-reduce 函数。
var allObs = col.map(function(img) {
var obs = img.reduceRegion(
{geometry: aoi, reducer: ee.Reducer.median(), scale: 30});
return img.set('NBR', obs.get('NBR'));
});
在您自己的应用中，如果 AOI 是大面积区域，您可以考虑增加 scale 并为 reduceRegion 函数指定 bestEffort、maxPixels、tileScale 实参，以确保操作不会超出像素数上限、内存或超时限制。您还可以将中位数归约函数替换为您偏好的统计函数。如需了解详情，请参阅《开发者指南》的图像区域的统计信息部分。
ui.Chart.feature.groups 函数现在可以接受集合。
该函数需要一个集合以及要映射到图表的集合对象属性的名称。在此示例中，“system:time_start”属性用作 x 轴变量，而“NBR”用作 y 轴变量。请注意，NBR 属性是由上述 reduceRegion 函数根据要缩减的图片中的波段名称命名的。如果您使用其他频段（而非“NBR”），请相应地更改 y 轴属性名称实参。最后，将分组（序列）变量设置为“SATELLITE”，并为其分配不同的颜色。
var chartAllObs =
ui.Chart.feature.groups(allObs, 'system:time_start', 'NBR', 'SATELLITE')
.setChartType('ScatterChart')
.setSeriesNames(['TM', 'ETM+', 'OLI'])
.setOptions({
title: 'All Observations',
colors: ['f8766d', '00ba38', '619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
pointSize: 6,
dataOpacity: 0.5
});
print(chartAllObs);
经过一段时间的处理后，控制台中会显示类似于图 2 的图表。请注意以下几点：
- 每年有多次观测。
- 该时间序列由三个不同的传感器组成。
- 没有可辨别的传感器偏差。
- 观测频率每年都会有所不同。
- 存在一些年内差异，这些差异在各年间保持一致且相对较小。如果使用 NDVI（img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');）尝试此操作，则会发现年内和年际响应变异性要大得多。
- 对于大多数时间序列，NBR 保持在较高水平，只有一些微小的扰动。
- 森林火灾导致 NBR 响应大幅快速下降，这在
- 不过，请注意，火灾（多拉湖）发生在 2011 年 9 月。检测年份出现错误的原因是，年度综合日期范围为 7 月至 8 月；在此期间之后发生的更改不会被纳入，直到下一次可用的非 null 观测结果（可能是下一年或更晚）出现为止。
- 在检测到 NBR 大幅下降后两年，植被开始恢复（NBR 响应增加）。
图 2. 单个像素的光谱响应时间序列图，其中包含来自 Landsat TM、ETM+ 和 OLI 的表示。TM 和 ETM+ 影像通过线性转换与 OLI 影像保持一致。图片拍摄于 7 月至 8 月，并经过滤选，确保画质上乘。
制作显示年度中位数的时序图表
为了简化时序并消除其中的噪声，可以应用年内观测缩减。此处使用的是中位数。
第一步是按年份对图片进行分组。通过映射每个图片的 ee.Image.Date 中的集合设置“year”，为每个图片添加新的“year”属性。
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
新的“year”属性可用于通过联接对同一年份的图片进行分组。 不同图片年份集合与完整集合之间的联接可提供同年列表，从中可以执行中位数缩减。将完整集合子集化为一组不同的年份代表。
var distinctYearCol = col.distinct('year');
定义一个过滤条件和一个联接，用于确定不同的年份集合 (distinctYearCol) 与完整集合 (col) 之间匹配的年份。匹配项将保存到名为“year_matches”的属性中。
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
应用联接并将生成的 FeatureCollection 转换为 ImageCollection。
var joinCol = ee.ImageCollection(join.apply(distinctYearCol, col, filter));
// Apply median reduction among matching year collections.
var medianComp = joinCol.map(function(img) {
var yearCol = ee.ImageCollection.fromImages(img.get('year_matches'));
return yearCol.reduce(ee.Reducer.median())
.set('system:time_start', img.date().update({month: 8, day: 1}));
});
如果您希望使用中位数以外的统计信息来表示特定年份的图片集合，请更改上方的 ee.Reducer。
最后，创建一个图表，用于显示来自年内夏季观测结果集的 NBR 中位数值。根据需要更改区域缩减统计信息。
var chartMedianComp = ui.Chart.image
.series({
imageCollection: medianComp,
region: aoi,
reducer: ee.Reducer.median(),
scale: 30,
xProperty: 'system:time_start',
})
.setSeriesNames(['NBR Median'])
.setOptions({
title: 'Intra-annual Median',
colors: ['619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
lineWidth: 6
});
print(chartMedianComp);
经过一段时间的处理后，控制台中会显示类似于图 3 的图表。
图 3. 单个像素的夏季中位数光谱响应时间序列图，其中包含来自 Landsat TM、ETM+ 和 OLI 的表示。TM 和 ETM+ 影像通过线性转换与 OLI 影像保持一致。图片拍摄于 7 月至 8 月，并经过高质量过滤。
其他信息
替代转换函数
Roy 等人 (2016) 的表 2 提供了 OLS 和 RMA 回归系数，用于将 ETM+ 地表反射率转换为 OLI 地表反射率，反之亦然。上述教程仅演示了通过 OLS 回归实现 ETM+ 到 OLI 的转换。您可以在代码编辑器脚本或 GitHub 源代码脚本中找到所有翻译选项对应的函数。
Dollar Lake 火灾
本教程中提到的火灾是 2011 年 9 月中旬在美国俄勒冈州胡德山北坡发生的 Dollar Lake 火灾。如需详细了解此功能，请参阅这篇博文，并参阅 Varner 等人 (2015)。