简介：这份资源面向遥感、GIS与生态监测方向的研究者及学生，提供一套基于Google Earth Engine平台与Landsat卫星影像的遥感生态指数自动化计算系统，用于解决多源遥感数据预处理繁琐、缨帽变换系数难以匹配、主成分分析方向判定主观等实际问题。压缩包共4个文件，约40KB，包含js核心脚本、md说明文档、txt使用说明与docx附赠资料，分别承载系统源码、操作指引与补充材料，便于快速部署与二次修改。系统集成多源遥感数据预处理、缨帽变换系数自适应匹配、主成分分析正负判定逻辑及年度合成等关键环节，可支撑长期生态变化研究与区域环境评估。目前已有94人学习下载，适合具备一定GEE与遥感基础、希望提升生态指数计算效率与结果科学性的读者参考使用。
1. 遥感生态指数自动化：从 Landsat 到缨帽变换的工程化落地
做遥感生态指数（RSEI）的人大多经历过这样的场景：手头攒了十几年的 Landsat 影像，想算一个区域的生态质量变化趋势，结果光是数据预处理就耗掉大半时间——云掩膜、大气校正、不同传感器之间的缨帽变换系数还不一样，算完主成分分析又得手动判断第一主成分的正负方向。这套流程跑一遍两三天，换个研究区又得重来。基于 Google Earth Engine 平台与 Landsat 卫星影像的遥感生态指数自动化计算系统，要解决的就是这个重复劳动问题：把多源遥感数据预处理、缨帽变换系数自适应匹配、主成分分析正负判定逻辑、年度合成这几个环节串成一条可复用的流水线。它适合已经了解 RSEI 基本概念、想在 GEE 上做长时间序列分析的研究生和一线技术人员，也适合需要批量出图的资源环境监测岗位。读完你能拿到一套可直接改研究区就跑的代码框架，以及几个我踩过的参数坑。
2. 多源遥感数据预处理：Landsat 5/7/8/9 在 GEE 里怎么统一
2.1 传感器差异带来的三个硬骨头
Landsat 系列跨越了 TM、ETM+、OLI/TIRS 两代传感器，波段编号、空间分辨率、辐射定标方式都不一样。做长时间序列 RSEI 时，最直接的问题是：Landsat 5/7 的热红外波段是 120m/60m 重采样到 30m，Landsat 8/9 的热红外是 100m 重采样到 30m，而 RSEI 里的湿度分量要用到缨帽变换的湿度波段，绿度和热度分别来自不同波段组合。如果不做统一，年度合成时会出现明显的条带或突变。
常见做法是在 GEE 里用ee.ImageCollection的merge把不同传感器的集合拼起来，但拼之前必须做三件事：统一波段名称、统一辐射定标、统一云掩膜策略。波段名称不统一，后续缨帽变换系数匹配就会错位；辐射定标不统一，不同年份的反射率量级对不上；云掩膜策略不统一，年度合成时有效像元数差异会很大。
2.2 用 GEE 的 Landsat 集合做统一预处理
下面这段代码是我常用的预处理骨架，核心思路是：先按传感器分组做云掩膜和辐射定标，再统一波段名，最后合并。
// 定义研究区和时间范围 var roi = ee.Geometry.Rectangle([110.0, 30.0, 112.0, 32.0]); var startYear = 2000; var endYear = 2023; // Landsat 5/7 的云掩膜函数（基于 QA_PIXEL 波段） function maskL57(image) { var qa = image.select('QA_PIXEL'); // 云、云阴影、雪、水都掩掉 var mask = qa.bitwiseAnd(1 << 3).eq(0) .and(qa.bitwiseAnd(1 << 4).eq(0)) .and(qa.bitwiseAnd(1 << 5).eq(0)); return image.updateMask(mask); } // Landsat 8/9 的云掩膜函数 function maskL89(image) { var qa = image.select('QA_PIXEL'); var mask = qa.bitwiseAnd(1 << 3).eq(0) .and(qa.bitwiseAnd(1 << 4).eq(0)) .and(qa.bitwiseAnd(1 << 5).eq(0)); return image.updateMask(mask); } // 统一波段名：把 SR_B1~SR_B7 映射成 B1~B7 function renameBands(image) { return image.select( ['SR_B1','SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7'], ['B1','B2','B3','B4','B5','B6','B7'] ); } // 按传感器分别处理再合并 var l5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2') .filterBounds(roi).filterDate(startYear+'-01-01', endYear+'-12-31') .map(maskL57).map(renameBands); var l7 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2') .filterBounds(roi).filterDate(startYear+'-01-01', endYear+'-12-31') .map(maskL57).map(renameBands); var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .filterBounds(roi).filterDate(startYear+'-01-01', endYear+'-12-31') .map(maskL89).map(renameBands); var l9 = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2') .filterBounds(roi).filterDate(startYear+'-01-01', endYear+'-12-31') .map(maskL89).map(renameBands); var merged = l5.merge(l7).merge(l8).merge(l9);逻辑说明：maskL57和maskL89都基于 QA_PIXEL 波段做位运算，但 Landsat 5/7 和 8/9 的 QA 位定义有细微差别，这里统一用云、云阴影、雪三个位。renameBands把 SR_B1 到 SR_B7 重命名为 B1 到 B7，这样后续缨帽变换时不用再判断传感器类型。合并后的集合里，所有影像的波段名一致，但辐射定标已经由 C02 T1_L2 产品完成，反射率量级在 0-1 之间。
参数说明：1 << 3是位运算，表示第 3 位（云置信度），1 << 4是云阴影，1 << 5是雪。如果你研究区有大量水体，建议把水也掩掉，加一个1 << 7。时间范围按需改，但注意 Landsat 5 在 2012 年后数据质量下降，Landsat 7 有 SLC-off 条带，实际做年度合成时建议统计有效像元比例，低于 30% 的年份标记为不可用。
2.3 年度合成前的有效像元筛选
合并后的集合不能直接做年度合成，因为有些年份云太多，合成结果全是噪声。我一般会先算每个年份的有效像元数，再决定哪些年份参与后续分析。
// 按年份统计有效像元数 var years = ee.List.sequence(startYear, endYear); var validCount = years.map(function(y) { var yearCol = merged.filterDate(y+'-01-01', y+'-12-31'); var count = yearCol.select('B1').count(); var stats = count.reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 30, maxPixels: 1e13 }); return ee.Feature(null, {year: y, validPixels: stats.get('B1')}); }); print('年度有效像元统计', ee.FeatureCollection(validCount));这段代码输出每个年份的平均有效像元数。如果某年有效像元数低于总像元数的 30%，我会在后续分析中跳过该年，或者用相邻年份插值。注意maxPixels要设大一点，否则大区域会报错。
3. 缨帽变换系数自适应匹配：不同传感器怎么选对系数
3.1 缨帽变换系数的来源和差异
缨帽变换（Tasseled Cap Transformation）把 Landsat 的多个波段线性组合成亮度、绿度、湿度三个分量。RSEI 里的湿度指标直接来自缨帽变换的湿度分量，绿度指标来自绿度分量。问题在于：Landsat 5 TM、Landsat 7 ETM+、Landsat 8 OLI 的缨帽变换系数完全不同，因为波段响应函数变了。如果你用 Landsat 8 的系数去算 Landsat 5 的湿度，结果会偏得离谱。
常见做法是查表：TM 用 Crist 1985 的系数，ETM+ 用 Huang 2002 的系数，OLI 用 Baig 2014 的系数。但在 GEE 里，这些系数需要手动写成ee.Array或ee.Image的矩阵形式，然后按传感器类型分别应用。
3.2 在 GEE 里实现系数自适应匹配
下面这段代码定义了三套系数，并用影像的传感器属性自动匹配。
// 缨帽变换系数（反射率产品，0-1 量级） var coeffs = { 'TM': { brightness: [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863], greenness: [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800], wetness: [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572] }, 'ETM': { brightness: [0.3561, 0.3972, 0.3904, 0.6966, 0.2286, 0.1596], greenness: [-0.3344, -0.3544, -0.4556, 0.6966, -0.0242, -0.2630], wetness: [0.2626, 0.2141, 0.0926, 0.0656, -0.7629, -0.5388] }, 'OLI': { brightness: [0.3029, 0.2786, 0.4733, 0.5599, 0.5080, 0.1872], greenness: [-0.2941, -0.2430, -0.5424, 0.7276, 0.0713, -0.1608], wetness: [0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559] } }; // 根据传感器类型应用对应系数 function applyTCT(image) { var sensor = image.get('SENSOR_ID'); var coeff = ee.Dictionary(coeffs).get(sensor); // 如果传感器不在表里，默认用 OLI coeff = ee.Algorithms.If(coeff, coeff, coeffs['OLI']); var c = ee.Dictionary(coeff); var brightness = image.select(['B1','B2','B3','B4','B5','B7']) .multiply(ee.Image.constant(c.get('brightness'))) .reduce(ee.Reducer.sum()); var greenness = image.select(['B1','B2','B3','B4','B5','B7']) .multiply(ee.Image.constant(c.get('greenness'))) .reduce(ee.Reducer.sum()); var wetness = image.select(['B1','B2','B3','B4','B5','B7']) .multiply(ee.Image.constant(c.get('wetness'))) .reduce(ee.Reducer.sum()); return image.addBands(brightness.rename('brightness')) .addBands(greenness.rename('greenness')) .addBands(wetness.rename('wetness')); } var withTCT = merged.map(applyTCT);逻辑说明：coeffs对象里存了三套系数，每套包含亮度、绿度、湿度三个数组。applyTCT先读影像的SENSOR_ID属性，然后从字典里取对应系数。如果传感器类型不在表里（比如 Landsat 9 的 SENSOR_ID 是 'OLI-2'），默认用 OLI 系数。注意这里用的是 B1、B2、B3、B4、B5、B7 六个波段，对应 TM/ETM+ 的 1-5、7 波段和 OLI 的 2-7 波段（OLI 的 B1 是海岸气溶胶，不参与缨帽变换）。
参数说明：系数数组的顺序必须和select里的波段顺序一致。TM 和 ETM+ 的系数来自不同文献，我一般用 TM 的 Crist 1985 和 ETM+ 的 Huang 2002，OLI 用 Baig 2014。如果你用的是地表反射率产品（LaSRC 或 LEDAPS），系数可以直接用；如果用的是 TOA 反射率，系数需要微调，但差异不大。
3.3 系数匹配错误的典型表现
如果系数匹配错了，最明显的表现是湿度分量出现大面积负值或异常高值。比如用 OLI 系数算 TM 影像，湿度分量会整体偏低，导致 RSEI 里的湿度指标权重被压缩。我一般会在应用系数后，先统计湿度分量的均值和标准差，如果均值偏离 0.1 太远，就说明系数可能不对。
4. 主成分分析正负判定逻辑：第一主成分方向怎么定
4.1 PCA 在 RSEI 里的作用和正负问题
RSEI 的核心是把绿度、湿度、热度、干度四个指标通过主成分分析降维，取第一主成分作为生态指数。但 PCA 有个经典问题：第一主成分的特征向量方向不固定，可能整体为正，也可能整体为负。如果方向反了，RSEI 值高的地方反而代表生态差，后续所有分析全错。
常见做法是：对第一主成分的特征向量，检查绿度和湿度对应的系数符号。如果绿度和湿度系数为负，就把第一主成分乘以 -1。这个逻辑在 GEE 里可以用ee.Array的eigen分解实现。
4.2 在 GEE 里实现 PCA 和正负判定
下面这段代码对年度合成的四个指标做 PCA，并自动判定方向。
// 假设已经有一个年度合成影像 annualImage，包含 greenness, wetness, heat, dryness 四个波段 function computeRSEI(annualImage) { // 标准化：减去均值除以标准差 var mean = annualImage.reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 30, maxPixels: 1e13 }); var std = annualImage.reduceRegion({ reducer: ee.Reducer.stdDev(), geometry: roi, scale: 30, maxPixels: 1e13 }); var normalized = annualImage.subtract(ee.Image.constant(mean.values())) .divide(ee.Image.constant(std.values())); // 提取四个波段的数组 var array = normalized.toArray(); // 计算协方差矩阵 var covar = array.reduceRegion({ reducer: ee.Reducer.covariance(), geometry: roi, scale: 30, maxPixels: 1e13 }); // 特征分解 var covarArray = ee.Array(covar.get('array')); var eigens = covarArray.eigen(); var eigenVectors = eigens.slice(1, 0, 1); // 第一主成分特征向量 // 判定正负：检查绿度和湿度对应的系数 var greennessCoeff = eigenVectors.get([0, 0]); var wetnessCoeff = eigenVectors.get([1, 0]); var sign = ee.Number(greennessCoeff).add(wetnessCoeff).lt(0).multiply(-2).add(1); // 计算第一主成分 var pc1 = array.multiply(eigenVectors).reduce(ee.Reducer.sum(), [0]); var rsei = pc1.multiply(sign); // 归一化到 0-1 var minMax = rsei.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: roi, scale: 30, maxPixels: 1e13 }); var rseiNorm = rsei.subtract(minMax.get('min')) .divide(ee.Number(minMax.get('max')).subtract(minMax.get('min'))); return rseiNorm.rename('RSEI'); }逻辑说明：先对四个指标做标准化，然后转成数组，计算协方差矩阵，再做特征分解。eigenVectors取第一列，即第一主成分的特征向量。sign的计算逻辑是：如果绿度和湿度系数之和小于 0，sign为 -1，否则为 1。最后把第一主成分乘以sign，再归一化到 0-1。
参数说明：covar.get('array')返回的是协方差矩阵的数组形式，eigen()返回的特征向量是按特征值降序排列的。slice(1, 0, 1)取第一列，对应第一主成分。注意reduceRegion的scale要和影像分辨率一致，否则协方差矩阵会失真。
4.3 正负判定失败的排查方法
如果 RSEI 结果和预期相反，先检查sign的计算。可以在代码里加一句print('sign', sign)，看输出是 1 还是 -1。如果符号对了但结果还是反的，可能是特征向量的顺序问题——有些实现里eigen()返回的特征向量是按行排列的，需要转置。我一般会打印特征向量矩阵，确认第一列对应的是第一主成分。
5. 年度合成与批量导出：从影像集合到 RSEI 时间序列
5.1 年度合成的三种策略
年度合成不是简单取平均。常见做法有三种：中值合成、最大值合成、均值合成。中值合成对云和异常值最稳健，我一般用中值。但如果某年有效像元太少，中值也会失真，这时候可以用相邻年份的均值插值。
// 按年份做中值合成 var annualImages = years.map(function(y) { var yearCol = withTCT.filterDate(y+'-01-01', y+'-12-31'); var median = yearCol.select(['greenness','wetness','heat','dryness']).median(); return median.set('year', y); });逻辑说明：filterDate按年份筛选，median()做中值合成。heat和dryness需要提前算好——热度一般用热红外波段的反演温度，干度用建筑指数和裸土指数的组合。这里假设你已经有了这四个波段。
5.2 批量导出 RSEI 结果
GEE 的导出任务需要一个个提交，但可以用循环批量提交。
// 批量导出 annualImages.forEach(function(img) { var year = img.get('year'); var rsei = computeRSEI(img); Export.image.toDrive({ image: rsei, description: 'RSEI_' + year, folder: 'RSEI_Export', scale: 30, region: roi, maxPixels: 1e13, fileFormat: 'GeoTIFF' }); });逻辑说明：forEach遍历每个年份的合成影像，调用computeRSEI算 RSEI，然后提交导出任务。description里带年份，方便后续整理。folder是 Google Drive 里的文件夹名，需要提前建好。
参数说明：scale设 30 和 Landsat 分辨率一致，maxPixels设大一点避免报错。如果研究区很大，建议分块导出，否则单个文件可能超过 10GB。
6. 避坑与排查：RSEI 自动化计算里最容易翻车的五个地方
6.1 缨帽变换系数用错传感器
现象：湿度分量整体偏低或偏高，RSEI 空间分布和实际生态状况对不上。 原因：Landsat 5/7/8/9 的缨帽变换系数不同，混用会导致湿度分量计算错误。 解决：在applyTCT里严格按SENSOR_ID匹配系数，Landsat 9 的SENSOR_ID是 'OLI-2'，需要单独加一套系数或默认用 OLI。
6.2 PCA 正负判定逻辑写反
现象：RSEI 高值区对应的是裸地或建筑区，低值区对应的是植被区。 原因：第一主成分的特征向量方向反了，sign计算逻辑写错。 解决：打印eigenVectors和sign，确认绿度和湿度系数为正时sign为 1。如果特征向量矩阵是转置的，调整slice的参数。
6.3 年度合成时有效像元不足
现象：某些年份的 RSEI 结果全是噪声，或者出现大面积空值。 原因：该年份云太多，中值合成后有效像元太少。 解决：在合成前统计有效像元比例，低于 30% 的年份用相邻年份插值，或者直接跳过。
6.4 导出任务卡在排队
现象：提交了几十个导出任务，但一直显示 'READY' 或 'RUNNING'，不完成。 原因：GEE 的导出任务有并发限制，同时提交太多会排队。 解决：分批提交，每次不超过 10 个任务。或者用Export.image.toAsset先存到 GEE 资产里，再批量下载。
6.5 研究区跨度过大导致协方差矩阵失真
现象：PCA 结果和预期不符，第一主成分解释方差比例很低。 原因：研究区太大，不同区域的生态状况差异大，协方差矩阵不能代表整体。 解决：分区域做 PCA，或者用分层抽样先选训练区，再应用到整个研究区。
7. 进阶技巧：用 GEE 的ee.Reducer做 RSEI 趋势分析和验证
7.1 用 Mann-Kendall 检验做趋势分析
算出 RSEI 时间序列后，下一步通常是做趋势分析。GEE 里没有现成的 Mann-Kendall 检验，但可以用ee.Reducer自定义。我一般用 Sen's slope 加 Mann-Kendall 的组合，代码量不大，但能给出每个像元的趋势方向和显著性。
// 简化版 Sen's slope：用线性回归的斜率代替 var trend = ee.ImageCollection(annualImages.map(function(img) { return computeRSEI(img).set('year', img.get('year')); })).select(['RSEI'], ['RSEI']).reduce(ee.Reducer.linearFit()); // trend 的第一个波段是斜率，第二个是截距逻辑说明：linearFit返回两个波段，第一个是斜率，第二个是截距。斜率大于 0 表示 RSEI 上升，生态改善；小于 0 表示下降。这个方法比 Mann-Kendall 简单，但对非线性趋势不敏感。如果要做严格检验，建议导出到本地用 Python 的pymannkendall库。
7.2 用高分辨率影像做验证
RSEI 的验证是个老大难问题。我一般用两种方法：一是和已有的土地覆盖产品做交叉验证，比如用 ESA WorldCover 或 CLCD 数据，看 RSEI 高值区是否对应林地、低值区是否对应建设用地；二是用 Google Earth 的高分辨率影像做目视验证，随机选 50-100 个点，人工判断生态状况，再和 RSEI 值做相关性分析。
// 随机选验证点 var validationPoints = ee.FeatureCollection.randomPoints(roi, 100); // 提取 RSEI 值 var rseiValue = rseiNorm.reduceRegions({ collection: validationPoints, reducer: ee.Reducer.first(), scale: 30 }); // 导出到 Drive 做后续分析 Export.table.toDrive({ collection: rseiValue, description: 'RSEI_Validation', fileFormat: 'CSV' });逻辑说明：randomPoints在研究区里随机选 100 个点，reduceRegions提取每个点的 RSEI 值，导出 CSV 后在本地和人工判读结果做对比。如果相关性低于 0.6，说明 RSEI 计算可能有问题，需要回头检查 PCA 或缨帽变换。
7.3 我踩过的一个坑
有一次做某流域的 RSEI 趋势分析，结果发现 2013 年之后 RSEI 突然下降。排查了半天，最后发现是 Landsat 8 发射后，数据源从 TM 切到 OLI，缨帽变换系数虽然换了，但热度指标的计算方式没跟着改——TM 的热红外波段是 B6，OLI 是 B10，波段编号变了但代码里没更新。这个坑让我养成了一个习惯：每次换传感器，先把所有波段编号和系数表打印出来核对一遍。希望这个经验能帮你少走点弯路。
本文还有配套的精品资源，点击获取
