目录一、研究区与数据集加载二、四大生态指标计算函数一波段缩放校正二四大生态指标逐一声明与计算1. 湿度指数Wetness2. 绿度指数Greenness3. 干度指数Dryness4. 热度指数Heat三函数返回值三、年尺度指标均值合成四、分析参数与投影设置五、指标标准化六、指标融合与 PCA 分析一指标融合二PCA 前置处理三协方差矩阵与特征值分解四特征向量方向调整五PC1 提取六PCA 执行七、RSEI 最终计算与标准化八、结果导出九、结果可视化十、代码核心设计亮点与注意事项十一、运行结果若觉得代码对您的研究 / 项目有帮助欢迎点击打赏支持需要完整代码的朋友打赏后可在后台私信复制文章标题发给我我会尽快发您完整可运行代码感谢支持本代码是基于 Google Earth EngineGEE平台利用 Landsat 8/9 卫星影像计算年尺度遥感生态指数RSEI的实现代码核心采用均值合成适配年尺度分析需求通过绿度、湿度、干度、热度四大生态指标结合主成分分析PCA完成 RSEI 构建以下从代码逻辑分层、核心函数、关键计算步骤等方面逐段详细解析。相关遥感指数计算公式指数计算公式相关波段NDVI(NIR - Red) / (NIR Red)B5, B4 (Landsat8/9)WetTasseled Cap变换固定系数法多波段线性组合LST热红外波段反演单窗算法B10/B11 (Landsat8/9)NDBSI(SI IBI) / 2SI: (B5B3-B6)/(B5B3B6)IBI: [2B5/(B5B4) - (B2B6)/(B2B6)] / [2B5/(B5B4) (B2B6)/(B2B6)]一、研究区与数据集加载定义分析的空间范围加载并筛选符合年尺度分析的 Landsat 8/9 影像集是后续指数计算的基础。// 导入研究区 var geometry table; // 加载 Landsat 8/9 数据集 var landsat ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .merge(ee.ImageCollection(LANDSAT/LC09/C02/T1_L2)) .filterDate(2023-01-01, 2023-12-31) .filterBounds(geometry) .filter(ee.Filter.lt(CLOUD_COVER, 1));研究区导入geometry table表示从 GEE 资产中导入已上传的研究区矢量面需提前将研究区 SHP 等矢量文件上传至 GEE 资产并命名为table影像集加载与合并分别加载 Landsat 8LC08和 Landsat 9LC09的 C02/T1_L2 级产品地表反射率和地表温度产品已做辐射定标、大气校正等预处理通过.merge()合并为一个影像集消除单传感器的时间 / 空间限制影像筛选条件filterDate限定 2023 年全年影像适配年尺度分析filterBounds仅保留研究区内的影像减少无效计算filter(ee.Filter.lt(CLOUD_COVER, 1))筛选云量小于 1% 的影像最大限度保证影像质量避免云污染对指数计算的影响。二、四大生态指标计算函数定义自定义函数对影像集中每幅影像进行波段缩放校正和四大生态指标计算并将计算结果作为新波段添加到原影像是整个代码的核心计算模块。函数内分波段预处理和指数计算两部分。一波段缩放校正Landsat C02/T1_L2 产品的光学波段和热红外波段均有量化缩放系数需还原为真实的地表反射率和地表温度值代码如下function calculateIndices(image) { // 应用缩放因子光学波段还原为地表反射率0-1范围 var optical image.select([SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7]) .multiply(0.0000275).add(-0.2); // 热红外波段还原为星上辐射亮度后转换为地表温度初始为开尔文 var thermal image.select([ST_B10]) .multiply(0.00341802).add(149.0); // 光学波段重命名方便后续指数计算的表达式调用 optical optical.select( [SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7], [Blue, Green, Red, NIR, SWIR1, SWIR2] ); // 后续指数计算代码... }光学波段选择蓝、绿、红、近红外、短波红外 1、短波红外 2 波段通过multiply(0.0000275).add(-0.2)还原为地表反射率热红外波段选择 ST_B10 波段Landsat 8/9 的地表温度波段通过专属缩放因子还原为初始温度值波段重命名将光学波段按波长命名为 Blue/Green/Red/NIR/SWIR1/SWIR2避免后续表达式中波段名混淆。二四大生态指标逐一声明与计算四大指标是 RSEI 的核心构成绿度、湿度为正相关指标指标值越高生态状况越好干度、热度为负相关指标指标值越高生态状况越差代码中分别通过缨帽变换、归一化差值、指数组合等方法计算。1. 湿度指数Wetness采用Landsat 8/9 OLI 传感器专属的缨帽变换湿度分量公式计算是表征地表水分含量土壤湿度、植被冠层水分等的核心指标var wetness optical.expression( Blue * 0.1511 Green * 0.1973 Red * 0.3283 NIR * 0.3407 SWIR1 * (-0.7117) SWIR2 * (-0.4559), {波段映射} ).rename(wetness);原理缨帽变换通过线性组合将光学波段转换为具有生态意义的分量湿度分量的系数为 Landsat 8/9 专属校准值不可直接套用至 Landsat 5/7输出单波段影像命名为wetness。2. 绿度指数Greenness采用归一化植被指数NDVI表征是最经典的植被生长状况指标取值范围为 [-1,1]var ndvi optical.normalizedDifference([NIR, Red]).rename(greenness);原理normalizedDifference([NIR, Red])即(NIR−Red)/(NIRRed)植被覆盖度越高NDVI 值越接近 1输出单波段影像命名为greenness。3. 干度指数Dryness采用归一化建筑裸土指数NDBSI表征综合土壤指数SI和建筑指数IBI计算能同时反映地表裸土和人工建筑的干化程度// 1. 土壤指数SI var si optical.expression(((SWIR1 Red) - (NIR Blue)) / ((SWIR1 Red) (NIR Blue)),{波段映射}); // 2. 建筑指数IBI var ibi optical.expression((2*SWIR1/(SWIR1NIR) - (NIR/(NIRRed)Green/(GreenSWIR1))) / (2*SWIR1/(SWIR1NIR) (NIR/(NIRRed)Green/(GreenSWIR1))),{波段映射}); // 3. 干度指数NDBSISI和IBI取均值 var ndbsi si.add(ibi).divide(2).rename(dryness);原理SI 侧重表征裸土的干化特征IBI 侧重表征人工建筑的不透水面特征二者均值能更全面反映地表干度输出单波段影像命名为dryness。4. 热度指数Heat采用地表温度LST表征将热红外波段的开尔文温度转换为摄氏度℃更符合实际分析习惯var lst thermal.select([ST_B10]).subtract(273.15).rename(heat);原理开尔文温度转摄氏度的公式为℃输出单波段影像命名为heat。三函数返回值return image.addBands([wetness, ndvi, ndbsi, lst]);将计算得到的四大指标作为新波段添加到原 Landsat 影像中返回带新波段的影像保证后续能对指标进行批量合成分析。三、年尺度指标均值合成将上述自定义函数应用到筛选后的 Landsat 影像集对四大指标分别做年尺度均值合成得到研究区 2023 年全年的四大指标均值影像适配年尺度 RSEI 分析需求月 / 季尺度常用中值合成年尺度用均值更能反映全年平均生态状况。// 对影像集应用指数计算得到带四大指标波段的影像集 var withIndices landsat.map(calculateIndices); // 分别提取四大指标波段做均值合成得到年尺度均值影像 var greenness withIndices.select(greenness).mean(); var wetness withIndices.select(wetness).mean(); var heat withIndices.select(heat).mean(); var dryness withIndices.select(dryness).mean();.map(calculateIndices)对影像集中每一幅影像批量执行指标计算GEE 中map函数是对影像集 / 特征集进行批量处理的核心方法.select(波段名).mean()先提取指定指标的波段再通过mean()计算影像集的像素级均值得到单幅年尺度均值影像四大指标各生成一幅独立的均值影像。四、分析参数与投影设置定义后续计算和导出的空间参数保证结果的空间分辨率和投影一致性是 GEE 中空间分析的基础配置。// 研究区边界直接使用导入的矢量面 var region geometry; // 投影信息WGS84 UTM 48N需根据研究区经纬度调整 var projection EPSG:32648; // 空间分辨率Landsat 8/9的光学/热红外波段均为30m重采样后 var scale 30;projectionUTM 投影带需根据研究区的经度调整如东经 105-111° 为 UTM 48N避免投影变形scale设置为 30m与 Landsat 8/9 的原始分辨率一致保证计算精度。五、指标标准化定义标准化函数将四大指标的原始值归一化至 [0,1] 范围解决四大指标量纲不同的问题如 NDVI [-1,1]、LST [0,50℃]、湿度分量为任意浮点数为后续 PCA 分析消除量纲干扰PCA 对量纲敏感未标准化会导致方差大的指标主导分析结果。function standardizeForPCA(image) { var stats image.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: region, scale: scale, maxPixels: 1e13 }); var min ee.Number(stats.values().get(0)); var max ee.Number(stats.values().get(1)); return image.subtract(min).divide(max.subtract(min)); } // 对四大指标均值影像分别执行标准化 var greenness_std standardizeForPCA(greenness); var wetness_std standardizeForPCA(wetness); var heat_std standardizeForPCA(heat); var dryness_std standardizeForPCA(dryness);极值计算通过reduceRegion结合ee.Reducer.minMax()计算研究区内指标的像素级最小值和最大值maxPixels:1e13用于解除 GEE 的像素计算数量限制避免大研究区计算报错归一化公式采用最大 - 最小值标准化公式为标准化值原始值最小值最大值最小值将指标值映射至 [0,1]批量标准化对四大指标的均值影像分别执行标准化得到标准化后的影像命名后缀加_std区分。六、指标融合与 PCA 分析先将标准化后的四大指标融合为一幅多波段影像再通过主成分分析PCA对多波段影像进行降维提取第一主成分PC1作为 RSEI 的基础PC1 能解释四大指标的绝大部分方差最能反映整体生态状况并对 PCA 的特征向量进行方向调整保证 PC1 与生态状况正相关。一指标融合将四个标准化后的单波段影像拼接为一幅四波段影像波段顺序为绿度、湿度、热度、干度是 PCA 分析的输入数据var compositeImage ee.Image.cat([greenness_std, wetness_std, heat_std, dryness_std]);GEE 中ee.Image.cat()是影像拼接的核心方法波段顺序决定后续 PCA 特征向量的对应关系不可随意调整。二PCA 前置处理对融合后的多波段影像做中心化处理PCA 的必要步骤消除数据偏移的影响// 计算各波段的研究区均值 var meanDict image.reduceRegion({reducer: ee.Reducer.mean(), geometry: region, scale: scale, maxPixels: 1e13}); var means ee.Image.constant(meanDict.values(bandNames)); // 中心化像素值 - 波段均值 var centered image.subtract(means);三协方差矩阵与特征值分解PCA 的核心数学步骤通过协方差矩阵描述四大指标间的相关性再通过特征值分解提取特征值和特征向量// 影像转数组适配矩阵计算 var arrays centered.toArray(); // 计算像素级协方差矩阵 var covar arrays.reduceRegion({reducer: ee.Reducer.centeredCovariance(), geometry: region, scale: scale, maxPixels: 1e13}); // 特征值分解得到特征值和特征向量 var covarArray ee.Array(covar.get(array)); var eigens covarArray.eigen(); // 提取特征向量PCA的核心参数 var eigenVectors eigens.slice(1, 1);特征向量决定 PCA 各主成分的构成第一主成分PC1的特征向量系数对应四大指标的贡献度.slice(1,1)对特征值分解结果进行切片仅提取特征向量部分剔除特征值。四特征向量方向调整原始 PCA 的 PC1 方向可能与生态状况负相关需根据四大指标与生态的相关性对特征向量系数进行符号调整保证最终 PC1 值越高生态状况越好// 提取PC1对应四大指标的特征向量系数 var ndviCoef eigenVectors.get([0, 0]); // 绿度系数 var wetCoef eigenVectors.get([1, 0]); // 湿度系数 var lstCoef eigenVectors.get([2, 0]); // 热度系数 var ndbsiCoef eigenVectors.get([3, 0]); // 干度系数 // 构建调整矩阵对系数符号进行调整 var adjustMatrix ee.Array([ [ee.Number(ndviCoef).lt(0).multiply(2).subtract(1), 0, 0, 0], // 绿度系数为负则取反 [0, ee.Number(wetCoef).lt(0).multiply(2).subtract(1), 0, 0], // 湿度系数为负则取反 [0, 0, ee.Number(lstCoef).gt(0).multiply(2).subtract(1), 0], // 热度系数为正则取反 [0, 0, 0, ee.Number(ndbsiCoef).gt(0).multiply(2).subtract(1)] // 干度系数为正则取反 ]); // 调整特征向量 var adjustedEigenVectors eigenVectors.matrixMultiply(adjustMatrix);调整规则绿度 / 湿度正相关指标若特征向量系数为负将系数取反保证其对 PC1 的贡献为正热度 / 干度负相关指标若特征向量系数为正将系数取反保证其对 PC1 的贡献为负调整逻辑lt(0).multiply(2).subtract(1)是 GEE 中实现条件取反的常用方法系数满足条件时返回 - 1取反不满足时返回 1不变。五PC1 提取将调整后的特征向量与中心化后的影像数据做矩阵乘法提取第一主成分PC1同时输出 PC2-PC4备用一般仅用 PC1 计算 RSEIvar arrayImage arrays.toArray(1); var principalComponents eigenImage.matrixMultiply(arrayImage); return principalComponents.arrayProject([0]).arrayFlatten([[PC1, PC2, PC3, PC4]]);六PCA 执行将融合后的多波段影像传入 PCA 函数得到包含 PC1-PC4 的主成分影像var principalComponents calculatePCA(compositeImage);七、RSEI 最终计算与标准化以调整后的 PC1 为基础构建初始 RSEIRSEI₀并将其归一化至 [0,1] 范围得到最终的 RSEI 影像RSEI 值越接近 1研究区生态状况越好越接近 0生态状况越差。// 提取PC1作为初始RSEI var rsei0 principalComponents.select(PC1); // 最大-最小值标准化将RSEI₀映射至[0,1] var rsei rsei0.unitScale( rsei0.reduceRegion({reducer: ee.Reducer.min(), geometry: region, scale: scale, maxPixels: 1e13}).values().get(0), rsei0.reduceRegion({reducer: ee.Reducer.max(), geometry: region, scale: scale, maxPixels: 1e13}).values().get(0) );rsei0 principalComponents.select(PC1)因前期已对特征向量做方向调整PC1 可直接作为初始 RSEI无需额外计算unitScale(min, max)GEE 中内置的最大 - 最小值标准化函数等价于(x−min)/(max−min)将 RSEI₀的原始值映射至 [0,1]得到最终可直接解读的 RSEI 影像。八、结果导出将四大指标原始均值影像和最终 RSEI 影像导出至 Google 云端硬盘方便后续在 ArcGIS、ENVI 等软件中进行空间分析和可视化是 GEE 结果落地的关键步骤。以绿度指标和 RSEI 最终结果为例其余指标导出逻辑一致// 导出绿度原始指标 Export.image.toDrive({ image: greenness.clip(region), description: RSEI_greenness_raw, folder: RSEI_Results, scale: scale, crs: projection, region: region.geometry().bounds(), maxPixels: 1e13 }); // 导出最终RSEI结果 Export.image.toDrive({ image: rsei.clip(region), description: RSEI_final, folder: RSEI_Results, scale: scale, crs: projection, region: region.geometry().bounds(), maxPixels: 1e13 });导出参数说明image: 影像.clip(region)先对影像做研究区裁剪仅导出研究区内的结果减少文件大小description导出文件的名称需唯一folderGoogle 云端硬盘中存储结果的文件夹名会自动创建scale/crs与前期设置的分辨率和投影一致保证空间一致性region: region.geometry().bounds()以研究区的外接矩形作为导出范围GEE 导出影像要求为矩形范围maxPixels:1e13解除像素导出数量限制。注代码中对四大原始指标绿度、湿度、干度、热度和最终 RSEI 分别做了导出共 5 幅影像。九、结果可视化在 GEE 地图界面加载四大指标和最终 RSEI 的影像图层实现交互式可视化预览方便快速检查计算结果的合理性代码中对不同指标设置了专属的配色方案和显示范围。湿度指数采用动态显示范围通过reduceRegion计算研究区内湿度的实际极值避免固定范围导致的可视化失真wetness.reduceRegion({reducer: ee.Reducer.minMax(), geometry: region, scale: scale, maxPixels: 1e13}).evaluate(function (result) { Map.addLayer(wetness.clip(region), {min: result.wetness_min, max: result.wetness_max, palette: [#FFE4B5, #0000FF]}, Wetness); });.evaluate()将 GEE 服务器端的计算结果传回客户端实现动态参数赋值配色浅橙色低湿度→蓝色高湿度符合生态直观认知。绿度 / 干度 / 热度采用固定合理显示范围适配指标的常规取值区间配色与指标生态意义匹配绿度NDVImin:-1, max:1白色→绿色绿色越深表示植被覆盖度越高热度LSTmin:0, max:50蓝色→黄色→红色红色越深表示地表温度越高干度NDBSImin:-1, max:1深绿色→棕色棕色越深表示裸土 / 建筑占比越高。RSEI 最终结果min:0, max:1红色→黄色→绿色绿色越深表示生态状况越好是结果可视化的核心图层。地图视图定位Map.centerObject(region, 8)将 GEE 地图界面自动定位到研究区缩放级别为 8可根据研究区大小调整。十、代码核心设计亮点与注意事项代码核心设计亮点年尺度适配采用均值合成替代常规的中值合成更能反映年尺度的平均生态状况适合长时间序列的年际对比分析数据质量控制云量筛选阈值设为 1%最大限度保证影像质量避免云污染对指数计算的干扰标准化与 PCA 优化先对四大指标做 [0,1] 标准化再做 PCA 中心化消除量纲和数据偏移的双重影响同时对特征向量做方向调整保证 RSEI 与生态状况正相关符合 RSEI 的定义要求模块化设计将指标计算、标准化、PCA 分析封装为自定义函数代码可读性和复用性高可通过修改参数适配不同研究区、不同时间尺度如季尺度需调整filterDate和合成方式可视化与导出兼顾既实现 GEE 端的交互式预览又支持将结果导出至本地满足后续深入分析的需求。代码使用注意事项需提前将研究区矢量文件上传至 GEE 资产并命名为table否则会报table is not defined错误投影信息EPSG:32648需根据研究区的经度调整避免投影变形云量筛选阈值1%可根据研究区的云量特征调整如高云量地区可适当提高至 5%但需保证影像质量运行代码前需在 GEE 中开启 Google 云端硬盘的授权否则无法导出结果大研究区计算时可适当降低scale如重采样至 60m减少计算量避免 GEE 超时报错。十一、运行结果湿度指数Wetness绿度指数Greenness干度指数Dryness热度指数Heat遥感生态指数RSEI 计算结果点击RUN即可下载数据若觉得代码对您的研究 / 项目有帮助欢迎点击打赏支持需要完整代码的朋友打赏后可在后台私信复制文章标题发给我我会尽快发您完整可运行代码感谢支持
基于 GEE 实现年尺度遥感生态指数（RSEI）计算
目录一、研究区与数据集加载二、四大生态指标计算函数一波段缩放校正二四大生态指标逐一声明与计算1. 湿度指数Wetness2. 绿度指数Greenness3. 干度指数Dryness4. 热度指数Heat三函数返回值三、年尺度指标均值合成四、分析参数与投影设置五、指标标准化六、指标融合与 PCA 分析一指标融合二PCA 前置处理三协方差矩阵与特征值分解四特征向量方向调整五PC1 提取六PCA 执行七、RSEI 最终计算与标准化八、结果导出九、结果可视化十、代码核心设计亮点与注意事项十一、运行结果若觉得代码对您的研究 / 项目有帮助欢迎点击打赏支持需要完整代码的朋友打赏后可在后台私信复制文章标题发给我我会尽快发您完整可运行代码感谢支持本代码是基于 Google Earth EngineGEE平台利用 Landsat 8/9 卫星影像计算年尺度遥感生态指数RSEI的实现代码核心采用均值合成适配年尺度分析需求通过绿度、湿度、干度、热度四大生态指标结合主成分分析PCA完成 RSEI 构建以下从代码逻辑分层、核心函数、关键计算步骤等方面逐段详细解析。相关遥感指数计算公式指数计算公式相关波段NDVI(NIR - Red) / (NIR Red)B5, B4 (Landsat8/9)WetTasseled Cap变换固定系数法多波段线性组合LST热红外波段反演单窗算法B10/B11 (Landsat8/9)NDBSI(SI IBI) / 2SI: (B5B3-B6)/(B5B3B6)IBI: [2B5/(B5B4) - (B2B6)/(B2B6)] / [2B5/(B5B4) (B2B6)/(B2B6)]一、研究区与数据集加载定义分析的空间范围加载并筛选符合年尺度分析的 Landsat 8/9 影像集是后续指数计算的基础。// 导入研究区 var geometry table; // 加载 Landsat 8/9 数据集 var landsat ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .merge(ee.ImageCollection(LANDSAT/LC09/C02/T1_L2)) .filterDate(2023-01-01, 2023-12-31) .filterBounds(geometry) .filter(ee.Filter.lt(CLOUD_COVER, 1));研究区导入geometry table表示从 GEE 资产中导入已上传的研究区矢量面需提前将研究区 SHP 等矢量文件上传至 GEE 资产并命名为table影像集加载与合并分别加载 Landsat 8LC08和 Landsat 9LC09的 C02/T1_L2 级产品地表反射率和地表温度产品已做辐射定标、大气校正等预处理通过.merge()合并为一个影像集消除单传感器的时间 / 空间限制影像筛选条件filterDate限定 2023 年全年影像适配年尺度分析filterBounds仅保留研究区内的影像减少无效计算filter(ee.Filter.lt(CLOUD_COVER, 1))筛选云量小于 1% 的影像最大限度保证影像质量避免云污染对指数计算的影响。二、四大生态指标计算函数定义自定义函数对影像集中每幅影像进行波段缩放校正和四大生态指标计算并将计算结果作为新波段添加到原影像是整个代码的核心计算模块。函数内分波段预处理和指数计算两部分。一波段缩放校正Landsat C02/T1_L2 产品的光学波段和热红外波段均有量化缩放系数需还原为真实的地表反射率和地表温度值代码如下function calculateIndices(image) { // 应用缩放因子光学波段还原为地表反射率0-1范围 var optical image.select([SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7]) .multiply(0.0000275).add(-0.2); // 热红外波段还原为星上辐射亮度后转换为地表温度初始为开尔文 var thermal image.select([ST_B10]) .multiply(0.00341802).add(149.0); // 光学波段重命名方便后续指数计算的表达式调用 optical optical.select( [SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7], [Blue, Green, Red, NIR, SWIR1, SWIR2] ); // 后续指数计算代码... }光学波段选择蓝、绿、红、近红外、短波红外 1、短波红外 2 波段通过multiply(0.0000275).add(-0.2)还原为地表反射率热红外波段选择 ST_B10 波段Landsat 8/9 的地表温度波段通过专属缩放因子还原为初始温度值波段重命名将光学波段按波长命名为 Blue/Green/Red/NIR/SWIR1/SWIR2避免后续表达式中波段名混淆。二四大生态指标逐一声明与计算四大指标是 RSEI 的核心构成绿度、湿度为正相关指标指标值越高生态状况越好干度、热度为负相关指标指标值越高生态状况越差代码中分别通过缨帽变换、归一化差值、指数组合等方法计算。1. 湿度指数Wetness采用Landsat 8/9 OLI 传感器专属的缨帽变换湿度分量公式计算是表征地表水分含量土壤湿度、植被冠层水分等的核心指标var wetness optical.expression( Blue * 0.1511 Green * 0.1973 Red * 0.3283 NIR * 0.3407 SWIR1 * (-0.7117) SWIR2 * (-0.4559), {波段映射} ).rename(wetness);原理缨帽变换通过线性组合将光学波段转换为具有生态意义的分量湿度分量的系数为 Landsat 8/9 专属校准值不可直接套用至 Landsat 5/7输出单波段影像命名为wetness。2. 绿度指数Greenness采用归一化植被指数NDVI表征是最经典的植被生长状况指标取值范围为 [-1,1]var ndvi optical.normalizedDifference([NIR, Red]).rename(greenness);原理normalizedDifference([NIR, Red])即(NIR−Red)/(NIRRed)植被覆盖度越高NDVI 值越接近 1输出单波段影像命名为greenness。3. 干度指数Dryness采用归一化建筑裸土指数NDBSI表征综合土壤指数SI和建筑指数IBI计算能同时反映地表裸土和人工建筑的干化程度// 1. 土壤指数SI var si optical.expression(((SWIR1 Red) - (NIR Blue)) / ((SWIR1 Red) (NIR Blue)),{波段映射}); // 2. 建筑指数IBI var ibi optical.expression((2*SWIR1/(SWIR1NIR) - (NIR/(NIRRed)Green/(GreenSWIR1))) / (2*SWIR1/(SWIR1NIR) (NIR/(NIRRed)Green/(GreenSWIR1))),{波段映射}); // 3. 干度指数NDBSISI和IBI取均值 var ndbsi si.add(ibi).divide(2).rename(dryness);原理SI 侧重表征裸土的干化特征IBI 侧重表征人工建筑的不透水面特征二者均值能更全面反映地表干度输出单波段影像命名为dryness。4. 热度指数Heat采用地表温度LST表征将热红外波段的开尔文温度转换为摄氏度℃更符合实际分析习惯var lst thermal.select([ST_B10]).subtract(273.15).rename(heat);原理开尔文温度转摄氏度的公式为℃输出单波段影像命名为heat。三函数返回值return image.addBands([wetness, ndvi, ndbsi, lst]);将计算得到的四大指标作为新波段添加到原 Landsat 影像中返回带新波段的影像保证后续能对指标进行批量合成分析。三、年尺度指标均值合成将上述自定义函数应用到筛选后的 Landsat 影像集对四大指标分别做年尺度均值合成得到研究区 2023 年全年的四大指标均值影像适配年尺度 RSEI 分析需求月 / 季尺度常用中值合成年尺度用均值更能反映全年平均生态状况。// 对影像集应用指数计算得到带四大指标波段的影像集 var withIndices landsat.map(calculateIndices); // 分别提取四大指标波段做均值合成得到年尺度均值影像 var greenness withIndices.select(greenness).mean(); var wetness withIndices.select(wetness).mean(); var heat withIndices.select(heat).mean(); var dryness withIndices.select(dryness).mean();.map(calculateIndices)对影像集中每一幅影像批量执行指标计算GEE 中map函数是对影像集 / 特征集进行批量处理的核心方法.select(波段名).mean()先提取指定指标的波段再通过mean()计算影像集的像素级均值得到单幅年尺度均值影像四大指标各生成一幅独立的均值影像。四、分析参数与投影设置定义后续计算和导出的空间参数保证结果的空间分辨率和投影一致性是 GEE 中空间分析的基础配置。// 研究区边界直接使用导入的矢量面 var region geometry; // 投影信息WGS84 UTM 48N需根据研究区经纬度调整 var projection EPSG:32648; // 空间分辨率Landsat 8/9的光学/热红外波段均为30m重采样后 var scale 30;projectionUTM 投影带需根据研究区的经度调整如东经 105-111° 为 UTM 48N避免投影变形scale设置为 30m与 Landsat 8/9 的原始分辨率一致保证计算精度。五、指标标准化定义标准化函数将四大指标的原始值归一化至 [0,1] 范围解决四大指标量纲不同的问题如 NDVI [-1,1]、LST [0,50℃]、湿度分量为任意浮点数为后续 PCA 分析消除量纲干扰PCA 对量纲敏感未标准化会导致方差大的指标主导分析结果。function standardizeForPCA(image) { var stats image.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: region, scale: scale, maxPixels: 1e13 }); var min ee.Number(stats.values().get(0)); var max ee.Number(stats.values().get(1)); return image.subtract(min).divide(max.subtract(min)); } // 对四大指标均值影像分别执行标准化 var greenness_std standardizeForPCA(greenness); var wetness_std standardizeForPCA(wetness); var heat_std standardizeForPCA(heat); var dryness_std standardizeForPCA(dryness);极值计算通过reduceRegion结合ee.Reducer.minMax()计算研究区内指标的像素级最小值和最大值maxPixels:1e13用于解除 GEE 的像素计算数量限制避免大研究区计算报错归一化公式采用最大 - 最小值标准化公式为标准化值原始值最小值最大值最小值将指标值映射至 [0,1]批量标准化对四大指标的均值影像分别执行标准化得到标准化后的影像命名后缀加_std区分。六、指标融合与 PCA 分析先将标准化后的四大指标融合为一幅多波段影像再通过主成分分析PCA对多波段影像进行降维提取第一主成分PC1作为 RSEI 的基础PC1 能解释四大指标的绝大部分方差最能反映整体生态状况并对 PCA 的特征向量进行方向调整保证 PC1 与生态状况正相关。一指标融合将四个标准化后的单波段影像拼接为一幅四波段影像波段顺序为绿度、湿度、热度、干度是 PCA 分析的输入数据var compositeImage ee.Image.cat([greenness_std, wetness_std, heat_std, dryness_std]);GEE 中ee.Image.cat()是影像拼接的核心方法波段顺序决定后续 PCA 特征向量的对应关系不可随意调整。二PCA 前置处理对融合后的多波段影像做中心化处理PCA 的必要步骤消除数据偏移的影响// 计算各波段的研究区均值 var meanDict image.reduceRegion({reducer: ee.Reducer.mean(), geometry: region, scale: scale, maxPixels: 1e13}); var means ee.Image.constant(meanDict.values(bandNames)); // 中心化像素值 - 波段均值 var centered image.subtract(means);三协方差矩阵与特征值分解PCA 的核心数学步骤通过协方差矩阵描述四大指标间的相关性再通过特征值分解提取特征值和特征向量// 影像转数组适配矩阵计算 var arrays centered.toArray(); // 计算像素级协方差矩阵 var covar arrays.reduceRegion({reducer: ee.Reducer.centeredCovariance(), geometry: region, scale: scale, maxPixels: 1e13}); // 特征值分解得到特征值和特征向量 var covarArray ee.Array(covar.get(array)); var eigens covarArray.eigen(); // 提取特征向量PCA的核心参数 var eigenVectors eigens.slice(1, 1);特征向量决定 PCA 各主成分的构成第一主成分PC1的特征向量系数对应四大指标的贡献度.slice(1,1)对特征值分解结果进行切片仅提取特征向量部分剔除特征值。四特征向量方向调整原始 PCA 的 PC1 方向可能与生态状况负相关需根据四大指标与生态的相关性对特征向量系数进行符号调整保证最终 PC1 值越高生态状况越好// 提取PC1对应四大指标的特征向量系数 var ndviCoef eigenVectors.get([0, 0]); // 绿度系数 var wetCoef eigenVectors.get([1, 0]); // 湿度系数 var lstCoef eigenVectors.get([2, 0]); // 热度系数 var ndbsiCoef eigenVectors.get([3, 0]); // 干度系数 // 构建调整矩阵对系数符号进行调整 var adjustMatrix ee.Array([ [ee.Number(ndviCoef).lt(0).multiply(2).subtract(1), 0, 0, 0], // 绿度系数为负则取反 [0, ee.Number(wetCoef).lt(0).multiply(2).subtract(1), 0, 0], // 湿度系数为负则取反 [0, 0, ee.Number(lstCoef).gt(0).multiply(2).subtract(1), 0], // 热度系数为正则取反 [0, 0, 0, ee.Number(ndbsiCoef).gt(0).multiply(2).subtract(1)] // 干度系数为正则取反 ]); // 调整特征向量 var adjustedEigenVectors eigenVectors.matrixMultiply(adjustMatrix);调整规则绿度 / 湿度正相关指标若特征向量系数为负将系数取反保证其对 PC1 的贡献为正热度 / 干度负相关指标若特征向量系数为正将系数取反保证其对 PC1 的贡献为负调整逻辑lt(0).multiply(2).subtract(1)是 GEE 中实现条件取反的常用方法系数满足条件时返回 - 1取反不满足时返回 1不变。五PC1 提取将调整后的特征向量与中心化后的影像数据做矩阵乘法提取第一主成分PC1同时输出 PC2-PC4备用一般仅用 PC1 计算 RSEIvar arrayImage arrays.toArray(1); var principalComponents eigenImage.matrixMultiply(arrayImage); return principalComponents.arrayProject([0]).arrayFlatten([[PC1, PC2, PC3, PC4]]);六PCA 执行将融合后的多波段影像传入 PCA 函数得到包含 PC1-PC4 的主成分影像var principalComponents calculatePCA(compositeImage);七、RSEI 最终计算与标准化以调整后的 PC1 为基础构建初始 RSEIRSEI₀并将其归一化至 [0,1] 范围得到最终的 RSEI 影像RSEI 值越接近 1研究区生态状况越好越接近 0生态状况越差。// 提取PC1作为初始RSEI var rsei0 principalComponents.select(PC1); // 最大-最小值标准化将RSEI₀映射至[0,1] var rsei rsei0.unitScale( rsei0.reduceRegion({reducer: ee.Reducer.min(), geometry: region, scale: scale, maxPixels: 1e13}).values().get(0), rsei0.reduceRegion({reducer: ee.Reducer.max(), geometry: region, scale: scale, maxPixels: 1e13}).values().get(0) );rsei0 principalComponents.select(PC1)因前期已对特征向量做方向调整PC1 可直接作为初始 RSEI无需额外计算unitScale(min, max)GEE 中内置的最大 - 最小值标准化函数等价于(x−min)/(max−min)将 RSEI₀的原始值映射至 [0,1]得到最终可直接解读的 RSEI 影像。八、结果导出将四大指标原始均值影像和最终 RSEI 影像导出至 Google 云端硬盘方便后续在 ArcGIS、ENVI 等软件中进行空间分析和可视化是 GEE 结果落地的关键步骤。以绿度指标和 RSEI 最终结果为例其余指标导出逻辑一致// 导出绿度原始指标 Export.image.toDrive({ image: greenness.clip(region), description: RSEI_greenness_raw, folder: RSEI_Results, scale: scale, crs: projection, region: region.geometry().bounds(), maxPixels: 1e13 }); // 导出最终RSEI结果 Export.image.toDrive({ image: rsei.clip(region), description: RSEI_final, folder: RSEI_Results, scale: scale, crs: projection, region: region.geometry().bounds(), maxPixels: 1e13 });导出参数说明image: 影像.clip(region)先对影像做研究区裁剪仅导出研究区内的结果减少文件大小description导出文件的名称需唯一folderGoogle 云端硬盘中存储结果的文件夹名会自动创建scale/crs与前期设置的分辨率和投影一致保证空间一致性region: region.geometry().bounds()以研究区的外接矩形作为导出范围GEE 导出影像要求为矩形范围maxPixels:1e13解除像素导出数量限制。注代码中对四大原始指标绿度、湿度、干度、热度和最终 RSEI 分别做了导出共 5 幅影像。九、结果可视化在 GEE 地图界面加载四大指标和最终 RSEI 的影像图层实现交互式可视化预览方便快速检查计算结果的合理性代码中对不同指标设置了专属的配色方案和显示范围。湿度指数采用动态显示范围通过reduceRegion计算研究区内湿度的实际极值避免固定范围导致的可视化失真wetness.reduceRegion({reducer: ee.Reducer.minMax(), geometry: region, scale: scale, maxPixels: 1e13}).evaluate(function (result) { Map.addLayer(wetness.clip(region), {min: result.wetness_min, max: result.wetness_max, palette: [#FFE4B5, #0000FF]}, Wetness); });.evaluate()将 GEE 服务器端的计算结果传回客户端实现动态参数赋值配色浅橙色低湿度→蓝色高湿度符合生态直观认知。绿度 / 干度 / 热度采用固定合理显示范围适配指标的常规取值区间配色与指标生态意义匹配绿度NDVImin:-1, max:1白色→绿色绿色越深表示植被覆盖度越高热度LSTmin:0, max:50蓝色→黄色→红色红色越深表示地表温度越高干度NDBSImin:-1, max:1深绿色→棕色棕色越深表示裸土 / 建筑占比越高。RSEI 最终结果min:0, max:1红色→黄色→绿色绿色越深表示生态状况越好是结果可视化的核心图层。地图视图定位Map.centerObject(region, 8)将 GEE 地图界面自动定位到研究区缩放级别为 8可根据研究区大小调整。十、代码核心设计亮点与注意事项代码核心设计亮点年尺度适配采用均值合成替代常规的中值合成更能反映年尺度的平均生态状况适合长时间序列的年际对比分析数据质量控制云量筛选阈值设为 1%最大限度保证影像质量避免云污染对指数计算的干扰标准化与 PCA 优化先对四大指标做 [0,1] 标准化再做 PCA 中心化消除量纲和数据偏移的双重影响同时对特征向量做方向调整保证 RSEI 与生态状况正相关符合 RSEI 的定义要求模块化设计将指标计算、标准化、PCA 分析封装为自定义函数代码可读性和复用性高可通过修改参数适配不同研究区、不同时间尺度如季尺度需调整filterDate和合成方式可视化与导出兼顾既实现 GEE 端的交互式预览又支持将结果导出至本地满足后续深入分析的需求。代码使用注意事项需提前将研究区矢量文件上传至 GEE 资产并命名为table否则会报table is not defined错误投影信息EPSG:32648需根据研究区的经度调整避免投影变形云量筛选阈值1%可根据研究区的云量特征调整如高云量地区可适当提高至 5%但需保证影像质量运行代码前需在 GEE 中开启 Google 云端硬盘的授权否则无法导出结果大研究区计算时可适当降低scale如重采样至 60m减少计算量避免 GEE 超时报错。十一、运行结果湿度指数Wetness绿度指数Greenness干度指数Dryness热度指数Heat遥感生态指数RSEI 计算结果点击RUN即可下载数据若觉得代码对您的研究 / 项目有帮助欢迎点击打赏支持需要完整代码的朋友打赏后可在后台私信复制文章标题发给我我会尽快发您完整可运行代码感谢支持
