简介这份资源面向遥感、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 重采样到 30mLandsat 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 2002OLI 用 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的计算逻辑是如果绿度和湿度系数之和小于 0sign为 -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自定义。我一般用 Sens slope 加 Mann-Kendall 的组合代码量不大但能给出每个像元的趋势方向和显著性。// 简化版 Sens 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 的热红外波段是 B6OLI 是 B10波段编号变了但代码里没更新。这个坑让我养成了一个习惯每次换传感器先把所有波段编号和系数表打印出来核对一遍。希望这个经验能帮你少走点弯路。本文还有配套的精品资源点击获取 
长河君 · 2026/9/29 1:06:21 · 阅读 1.2k GEE遥感生态指数自动化：Landsat多源数据与缨帽变换工程化实践
简介这份资源面向遥感、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 重采样到 30mLandsat 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 2002OLI 用 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的计算逻辑是如果绿度和湿度系数之和小于 0sign为 -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自定义。我一般用 Sens slope 加 Mann-Kendall 的组合代码量不大但能给出每个像元的趋势方向和显著性。// 简化版 Sens 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 的热红外波段是 B6OLI 是 B10波段编号变了但代码里没更新。这个坑让我养成了一个习惯每次换传感器先把所有波段编号和系数表打印出来核对一遍。希望这个经验能帮你少走点弯路。本文还有配套的精品资源点击获取