GEE实现遥感生态指数RSEI自动化计算与PCA正负判定

简介：一套基于 Google Earth Engine 与 Landsat 影像的遥感生态指数自动化计算系统，面向从事生态监测、环境遥感与 GIS 分析的研究者，用于解决 RSEI 计算流程繁琐、参数依赖人工调试的问题。资源包仅 4 个文件、约 40KB，以 js 脚本、md 说明、txt 说明与 docx 附赠资料组成，其中 RSEI.js 为核心计算脚本，README 与说明文件用于指导部署调用，docx 文件补充操作与案例背景。内容覆盖多源遥感数据预处理、缨帽变换系数自适应匹配、主成分分析正负判定逻辑及年度合成等关键环节，便于研究者快速掌握自动化流程，也可按需修改代码适配不同研究区。当前已有 94 人学习，适合具备一定遥感基础、希望借助 GEE 提高生态指数计算效率的入门及进阶用户。

1. 遥感生态指数自动化计算：先想清楚四个分量，再谈系统

每年做区域生态质量评价，遥感生态指数（RSEI）都是绕不开的产出项。用 Google Earth Engine 平台加 Landsat 卫星影像搭一套自动化计算系统，听起来就是把缨帽变换、主成分分析、年度合成这些环节串起来跑一遍，但真正落地时，麻烦全藏在“自动化”三个字里。最典型的就是主成分分析正负判定：同一套代码，换一个年份或换一个区域，PC1 的符号就翻过来，RSEI 的生态含义直接反转。这篇笔记就按数据预处理、缨帽变换系数自适应匹配、PCA 判定逻辑、年度合成这条线，把能直接复现的 GEE 代码和参数选择讲清楚。适合正在做遥感生态指数、生态遥感指数评价的工程师和研究生，也适合想把单景结果扩展到逐年时序的人。

2. 年度合成与数据底座：Landsat Collection 2 在 GEE 上的预处理取舍

2.1 RSEI 四分量在 Landsat 影像上的遥感表达

RSEI 不是随便把几个指数平均，它用主成分分析把四个生态相关分量压缩成一个综合值。四个分量分别是绿度、湿度、热度、干度，对应遥感变量是 NDVI、缨帽变换的 Wet 分量、地表温度 LST、以及由裸土指数 SI 和建筑指数 IBI 合成的 NDBSI。这里有一个容易被忽视的点：湿度分量不是土壤含水量，而是缨帽变换的第三分量，它反映植被冠层和土壤的总体湿润程度，在 PCA 里通常对生态质量起正贡献。干度则相反，裸土和建设用地越多，NDBSI 越高，生态越差。

四个分量量纲差异很大，NDVI 在 -1 到 1，LST 是开尔文，NDBSI 是比值，直接加总没有意义。RSEI 的做法是先把四个分量标准化，再做 PCA，取第一主成分 PC1 作为基础，最后归一化到 0 到 1。高值代表生态质量好，低值代表生态压力大。换句话说，RSEI 的价值在于把多维遥感信息压成一个可比较的指数，而 PCA 的载荷符号决定了这个指数是否真的“高值即优”。这也是后面主成分分析正负判定逻辑存在的根本原因。

2.2 年度合成与云掩膜：中位数合成的取舍

单景 Landsat 很难覆盖完整研究区，尤其像长江中下游这种多云地区，一景影像里可能有三分之一被云挡住。年度合成的基本思路是把生长季内所有影像做时间维压缩，用中位数代替单景像元，这样既保留地表反射率的稳定状态，又规避云和阴影。常见做法是用imageCollection.median()，它对异常值不敏感，比均值更稳，也不需要假设数据分布。

云掩膜放在合成之前还是之后，结果差异很大。放在合成之前，能保证进中位数的像元都是有效观测；放在合成之后，就只能靠多波段规则去猜哪些像元是被云污染过的。我的习惯是用 Landsat Collection 2 Level 2 的 QA_PIXEL 波段，把云、云阴影、雪和填充像元全部掩掉，再做中位数合成。Landsat 5、7、8、9 的 QA 波段定义一致，这在自动化系统里非常重要，否则每换一个传感器就要写一套掩膜逻辑。

function maskLandsat(img) { var qa = img.select('QA_PIXEL'); // Collection 2 QA_PIXEL: bit3=cloud, bit4=cloud shadow, bit5=snow var mask = qa.bitwiseAnd(1 << 3).eq(0) .and(qa.bitwiseAnd(1 << 4).eq(0)) .and(qa.bitwiseAnd(1 << 5).eq(0)); return img.updateMask(mask); } var start = '2023-06-01'; var end = '2023-09-30'; var roi = ee.FeatureCollection('你的研究区').geometry(); var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')) .filterDate(start, end) .filterBounds(roi) .map(maskLandsat); var composite = l8.median().clip(roi);

这段代码里，1 << 3是把 1 左移三位，得到对应 QA_PIXEL 第 3 位的掩码值。.eq(0)表示该位不为 1，也就是没有云。合并 LC08 和 LC09 是因为 Landsat 8 和 9 同时运行时，时间分辨率能加密到 8 天，合成样本更充足。median()的参数不需要额外调，真正要调的是日期窗口：如果研究区在北方温带，6 到 9 月可以；如果研究区是热带雨林，这个窗口内云还是太多，就得放宽到全年，甚至改成逐月合成再取年最大。

2.3 自动化系统的模块划分：从影像到指数的数据流

把整个 RSEI 计算拆成模块，自动化才有意义。我的划分是五层：数据获取层负责按年份过滤 Landsat Collection 2，并做云掩膜和波段重命名；预处理层做缨帽变换、NDVI、LST、NDBSI，生成四个分量；合成层把四个分量各自合成年度影像；分析层做 PCA 和正负判定；输出层做归一化、调色板和导出。这五层对应到 GEE 代码里就是五个函数，每个函数只做一件事，方便单独调试。

数据流上最常出问题的是波段重命名。Landsat 5/7 和 8/9 的原始波段名不一样，同样叫“红波段”，L5 是 SR_B3，L8 是 SR_B4。如果不在预处理层统一成red,green,blue,nir,swir1,swir2，后面缨帽变换系数表就会对错位置。我一般会写一个renameBands函数，按SPACECRAFT_ID判断传感器，再统一映射。这个函数是整个系统的地基，地基歪了，后面 PCA 出来的载荷再好看也是错的。

3. 多源遥感数据预处理与缨帽变换系数自适应匹配

3.1 不同 Landsat 传感器的缨帽变换系数表

缨帽变换（Tasseled Cap）是把六个反射率波段投影到亮度、绿度、湿度三个新轴上。问题在于，Landsat 5 的 TM 传感器、Landsat 7 的 ETM+ 和 Landsat 8/9 的 OLI 光谱响应函数不同，同一套系数不能跨传感器通用。如果拿 TM 的湿度系数去算 OLI 影像，Wet 分量会带上明显的亮度信息，生态含义就乱了。自适应匹配的核心是一张按SPACECRAFT_ID索引的系数表，我整理过最常用的一组：

OLI 的系数来自 Baig 等人对 Landsat 8 的重算结果，Landsat 9 与 8 同为 OLI 传感器，系数可以直接复用。这张表里 TM 和 ETM+ 的系数基本一致，可以直接合并，但代码里我仍会分开存，方便以后单独修正。实际项目中，如果研究区同时横跨 L5、L7、L8、L9，自适应匹配就不只是调用一个系数数组，而是要保证每景影像在完成缨帽变换后再进入年度合成。

3.2 在 GEE 里做缨帽变换：波段顺序与参数说明

代码实现时，最容易犯的错是波段顺序。缨帽变换要求输入波段严格按蓝、绿、红、近红外、SWIR1、SWIR2 排列，系数表和影像波段必须一一对应。我在早期项目里把 SWIR1 和 SWIR2 写反过一次，结果湿度分量在城区出现大面积高值，排查了很久才意识到是波段顺序问题。下面这段代码先把波段统一命名，再按顺序选择，最后乘系数并求和。

function renameBands(img) { var sat = img.get('SPACECRAFT_ID'); if (sat === 'LANDSAT_5' || sat === 'LANDSAT_7') { return img.select( ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'], ['blue', 'green', 'red', 'nir', 'swir1', 'swir2'] ); } else { return img.select( ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'], ['blue', 'green', 'red', 'nir', 'swir1', 'swir2'] ); } } function tasseledCap(img) { var ren = renameBands(img); var sat = img.get('SPACECRAFT_ID'); var coef; if (sat === 'LANDSAT_8' || sat === 'LANDSAT_9') { coef = oliTC; } else { coef = tmTC; } var b = ren.select(['blue', 'green', 'red', 'nir', 'swir1', 'swir2']); var wetness = b.multiply(coef.wetness) .reduce(ee.Reducer.sum()) .rename('WET'); var greenness = b.multiply(coef.greenness) .reduce(ee.Reducer.sum()) .rename('GREEN'); var brightness = b.multiply(coef.brightness) .reduce(ee.Reducer.sum()) .rename('BRIGHT'); return ren.addBands([wetness, greenness, brightness]); }

这里b.multiply(coef.wetness)是逐波段乘系数，随后reduce(ee.Reducer.sum())把六个波段的加权和压成一个值。coef 数组长度必须是 6，顺序和 b 的 select 顺序一致。注意 Landsat Collection 2 的 SR 产品已经是反射率，不需要再做大气校正，但需要确认尺度因子：Collection 2 的 SR 波段已经乘以 10000，直接拿来算缨帽没问题；如果数据源是 TOA 或 DN 级影像，必须先换算成 0 到 1 的反射率，再乘系数，否则湿度分量会整体偏移。

提示：在 GEE 里调试缨帽变换时，不要直接 print 整景影像的像素值，而是用reduceRegion取一个均质区域的均值来对比。比如水域的 Wet 分量应该明显高于裸地，如果反过来，先检查波段顺序而不是系数表。

3.3 用 NDVI 阈值做水体裸地掩膜：年度合成前的最后一件事

四个分量算完后，直接丢进 PCA 之前，最好做一次生态无效像元掩膜。水体对缨帽湿度分量贡献极高，但水体本身不是“生态质量好”，如果大面积水域进入 PCA 样本，PC1 会被水体的高湿度主导，陆地生态信息反而被压缩。常见做法是用 NDVI 阈值把纯水体像元置为无效，或者单独生成一个水体掩膜波段。

function maskWater(img) { var ndvi = img.normalizedDifference(['nir', 'red']).rename('NDVI'); var waterMask = ndvi.lt(0); return img.updateMask(waterMask.not()); }

NDVI 小于 0 的像元通常是水体、积雪或极端裸地，在年度合成影像里，这个阈值对清澈水体效果很好，但对浑浊水体或阴影区域会漏掉一部分。更稳妥的做法是同时用 mNDWI 和 NDVI 做决策：mNDWI 大于 0.3 且 NDVI 小于 0 判为水体。不过每加一个阈值就多一个参数，自动化系统里建议把它们设计成可配置项，而不是写死在代码里。

我习惯在进入 PCA 之前，把这个掩膜影像的像素均值打印出来做 sanity check。比如湿润区的湿地项目，WET 分量均值应当在绿度分量均值的一倍以上；如果 WET 均值接近 0 甚至负值，多半是合成影像里云掩膜过度，把高湿度像元全删掉了。

4. 主成分分析正负判定逻辑：让 PCA 从黑匣子变成可解释的生态指数

4.1 在 GEE 里对四个分量做主成分分析的最小实现

GEE 里做 PCA 分三步：采样、协方差矩阵、特征分解。最小实现可以用reduceRegion直接算四个波段的协方差矩阵，再用eigen()做特征分解。这里的核心参数是采样像元数、尺度和采样区域。采样区域必须是研究区内部的非 NoData 区域，最好先用掩膜像元统计一遍，确保四个分量都有有效值。

var bands = ['GREEN', 'WET', 'LST', 'NDBSI']; var covar = image4.select(bands).reduceRegion({ reducer: ee.Reducer.covariance(), geometry: roi, scale: 30, maxPixels: 1e10, tileScale: 4 }).get('covariance'); var eigens = ee.Array(covar).eigen(); var eigenVectors = eigens.slice(1, 1); var pcImage = image4.select(bands) .toArray() .matrixMultiply(eigenVectors.transpose()) .arrayFlatten([['PC1', 'PC2', 'PC3', 'PC4']]); var pc1 = pcImage.select('PC1');

这段代码里，reduceRegion的reducer: ee.Reducer.covariance()返回一个数组，maxPixels设得很大是为了避免大区域采样时因为像元数超限报错。tileScale: 4是一个容易忽略的参数：它在计算协方差时把任务切分成更小的块，能有效避免 GEE 的 computation time out。eigens.slice(1, 1)取特征向量矩阵：eigen() 返回的数组第一列是特征值，后面的列是特征向量，切片要去掉第一列。

特征向量矩阵的行顺序和 bands 数组的列顺序是一致的，也就是第一行对应 GREEN 的载荷，第二行对应 WET，第三行 LST，第四行 NDBSI。这个顺序在下一步正负判定时非常关键。载荷载荷的绝对值越大，说明该分量对 PC1 的贡献越大；符号决定方向。

4.2 正负判定逻辑：什么时候该用 1-PC1

PCA 是纯数学变换，数据稍微变化，PC1 的方向就可能翻转。同一套代码，今年算出 PC1 与 WET 正相关，明年可能变成负相关。RSEI 的生态语义要求高值代表生态好，所以必须做正负判定。判定的原则是：看 PC1 与绿度 GREEN、湿度 WET 的载荷符号。如果 GREEN 和 WET 载荷为正、LST 和 NDBSI 载荷为负，PC1 直接就是生态质量的正向表达；如果反过来，PC1 就需要取反。

var pc1Row = eigenVectors.slice(1, 0, 1); // 第 1 行，即 PC1 的 4 个载荷 var greenLoading = ee.Number(pc1Row.get([0, 0])); var wetLoading = ee.Number(pc1Row.get([0, 1])); // 生态正向指标（绿度+湿度）载荷为负，说明 PC1 方向反了 var needFlip = greenLoading.add(wetLoading).lt(0); var rseiRaw = ee.Image(ee.Algorithms.If(needFlip, pc1.multiply(-1), pc1)) .rename('RSEI_raw');

这段代码的精髓在greenLoading.add(wetLoading).lt(0)。只要绿度和湿度两个正向分量的载荷之和小于零，就判定 PC1 方向反了。用两者之和而不是只看其中一个，是因为某些区域 WET 分量对 PC1 的贡献很小，单看 WET 可能发生误判。加和之后，即便某个分量贡献弱，另一个也能拉住方向。

判定逻辑写进自动化系统后，还有一个容易出问题的地方：阈值是“小于 0”还是“小于某个负数”。数据噪声可能导致载荷在 0 附近抖动，比如 GREEN 载荷 0.01，WET 载荷 -0.02，加和 -0.01，被判定为需要取反，但实际结果和上一年差别不大。我一般会加一个很小的缓冲阈值，比如载荷之和小于 -0.05 才取反。具体值要看研究区，但至少不要用裸的 0，否则边界年份容易出现 RSEI 空间分布突然整体反转。这个参数建议写进配置文件，标成signBuffer。

4.3 从 PC 影像到 RSEI 的归一化与出图参数

正负判定之后，RSEI 原始值还不是 0 到 1 的指数。GEE 里可以用unitScale做最小最大归一化，但要注意：最大最小值应该从研究区统计得到，而不是全图默认。如果研究区里有一片云残留的高亮像元，最大最小值会被它们拉偏，导致整个 RSEI 空间对比度变差。

var minMax = rseiRaw.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: roi, scale: 30, maxPixels: 1e10, bestEffort: true }); var rsei = rseiRaw.unitScale( ee.Number(minMax.get('RSEI_raw_min')), ee.Number(minMax.get('RSEI_raw_max')) ).rename('RSEI'); var palette = ['#d7191c', '#fdae61', '#ffffbf', '#a6d96a', '#1a9641']; Map.addLayer(rsei.clip(roi), {min: 0, max: 1, palette: palette}, 'RSEI');

调色板从红到绿，红色代表生态差、绿色代表生态好，这是生态评价最常用的视觉表达。输出参数上，scale: 30和原始 Landsat 分辨率一致，不建议为了省计算量直接降到 30 米以下的重采样。如果研究区很大，bestEffort: true可以让 GEE 自动调整采样范围，避免内存爆炸；如果 RSEI 要用于定量分析而不是出图，我建议关掉 bestEffort，改用固定 scale 和 maxPixels，确保结果可复现。

5. RSEI 计算的五条踩坑记录：从黑屏到跨年跳变

5.1 PCA 结果全黑，先查 NoData 混进采样

现象：PC1 影像大部分区域是黑色，统计直方图只有 0 或极少数值。原因：四个分量影像里有大量 NoData 像元，比如云掩膜把某块区域全部删掉，而reduceRegion在协方差计算时把 NoData 当作数字参与运算，协方差矩阵出现 NaN，特征分解后 PC1 全是无效值。解决：采样前先用掩膜把无效像元统一设为 NoData，并在 reduceRegion 里加maxPixels和tileScale；如果是用 sample 方式做 PCA，要加dropNulls: true。我最开始做自动化时只在一景影像上测试，没遇到这个问题，后来换到大范围年度合成时黑屏，排查半天才发现是云掩膜留下的空洞在作怪。

5.2 Landsat 7 SLC-off 条带让湿度分量翻车

现象：Landsat 7 参与年度合成后，WET 分量影像上出现整齐的斜向条带。原因：Landsat 7 的 SLC-off 故障导致扫描行校正器失效，影像存在永久条带，QA 波段里没有单独标识这些缝隙像元，中位数合成虽然能削弱条带，但在条带密集区域仍然会残留低值异常。解决：在云掩膜函数里，把 QA_PIXEL 的 bit 0 填充像元一并排除，同时用多时相影像做中位数时把观测次数阈值设高，比如少于 3 次有效观测的像元直接设为 NoData。这样能保证每个像元至少来自多次观测，条带区域不会被单次异常值主导。

5.3 跨年 RSEI 跳变：日期窗口和云掩膜都不一致

现象：逐年 RSEI 均值变化曲线在某一年突然跃升或骤降，图形上像台阶。原因：自动化系统里逐年合成窗口不统一，某年用的是 6 到 9 月，另一年用的是 5 到 10 月，物候差异被当成生态变化。解决：把日期窗口作为固定参数，例如startDay = '06-01'、endDay = '09-30'，年份只替换前缀；云掩膜参数也统一。另一个隐蔽原因是相邻年份用了不同等级的 Landsat 产品，比如前一年是 Collection 1，后一年是 Collection 2，两个版本的 SR 定标有差异，尽量全部使用 Collection 2 Tier 1。

5.4 PC1 符号反转：不同年份的正负判定必须固定

现象：某一年 RSEI 高值区在城市，低值区在森林，空间分布完全反转。原因：PCA 在不同年份的样本统计量不同，PC1 的特征向量符号自动翻转，而代码里写死了RSEI = 1 - PC1，没有做符号判定。解决：用第 4 章的needFlip逻辑，在每一年都重新计算 PC1 载荷，按绿度和湿度载荷方向决定取反还是保留。我曾经遇到过一种更隐蔽的情况：同一年里分两个行政区分区计算，一个区判定取反，另一个区不取反，边界处 RSEI 出现明显跳变。自动化系统如果要分区计算，正负判定的样本区域必须统一，不能每个小区域单独判定。

5.5 裸地出现高 RSEI：NDBSI 的 SWIR 波段映射错了

现象：荒漠或裸土区 RSEI 反而比林地高。原因：干度指标 NDBSI 里的 SWIR1 和 SWIR2 波段映射错误，Landsat 8 的 SR_B6 对应 SWIR1、SR_B7 对应 SWIR2，但有人沿用 Landsat 5 的习惯把 SR_B5 当 SWIR1，计算出的 IBI 在裸地区域变成负数，干度分量失去意义。解决：在统一波段命名时，Landsat 8/9 的swir1必须来自SR_B6，swir2来自SR_B7，并打印 NDBSI 在裸地和水体处的均值做验证。裸地 NDBSI 应该明显偏高，水体 NDBSI 应该偏低，如果相反，优先查波段映射而不是公式。

提示：以上任何一个坑，都可能让 RSEI 看起来合理但实际不可用。建议每次跑完都导出四分量影像，在 GEE 里叠加目视检查一轮。RSEI 是压过一次主成分的结果，原始分量如果错了，指数层面的数字再漂亮都救不回来。

6. 进阶：用 RSEI 时间序列做生态区划验证

6.1 用 linearFit 算五年 RSEI 斜率并分级

年度 RSEI 算完后，最直接的进阶用法是计算逐年斜率，识别生态变好和变差的区域。GEE 的linearFit可以快速算每个像元的时间趋势，输入是时间和 RSEI 的多波段影像。先把逐年 RSEI 合成一个 ImageCollection，再转成多波段影像，时间轴放 1 到 5。

var rseiCol = ee.ImageCollection.fromImages([ rsei2020, rsei2021, rsei2022, rsei2023, rsei2024 ]); var years = ee.List.sequence(1, 5); var rseiStack = rseiCol.toBands().rename(years.map(function(y) { return 'y' + y; }).getInfo()); var slope = rseiStack.select(years.map(function(y) { return 'y' + y; })).reduce(ee.Reducer.linearFit());

slope 影像的scale波段就是像元尺度的年均变化量，阈值按研究区实际情况定。我常用的分级是：大于 0.02 为明显变好，小于 -0.02 为明显变差，中间为稳定区。这个阈值不是通用的，但至少能区分出趋势信号和年际噪声，再结合土地利用数据看驱动因素。

6.2 抽样验证：把 RSEI 和实测样方做相关

RSEI 如果没有地面验证，始终只是影像统计结果。做生态区划前，我会在野外样方里记录植被覆盖度、土壤含水量、地表温度，然后和同年的 RSEI 做 Pearson 相关。如果相关系数低于 0.4，问题通常不在统计方法，而在四个分量里的某一个是错的，比如 LST 没做地表比辐射率校正，或者 NDBSI 在低植被覆盖区失效。验证不是把 R² 算完就结束，而是把 RSEI 的误差空间分布图叠加到原始影像上看，找到系统偏离的区域类型。

我做 RSEI 自动化系统最大的教训，是太相信 PCA 的“无监督”能力，认为它自己会找到生态质量的潜在结构。实际上，无监督意味着符号不确定，正负判定逻辑才是整个系统里最容易翻车但也最值得深挖的环节。现在我的代码里会显式输出每年的特征向量载荷，把needFlip结果写进导出文件名，比如RSEI_2023_flip_true.tif，这样后续任何一次结果异常都能回溯到当年的 PCA 方向。希望这套流程能帮你少踩几个坑，也帮你在 RSEI 这条路上走得更稳。

本文还有配套的精品资源，点击获取