遥感生态指数RSEI模型原理与GEE实现全解析

这些年做生态环境遥感评估，我周围很多同行一开始都习惯用NDVI单指标来看一个区域的生态状况，但说实话，单看绿度太片面了。一个城市周边可能有大量农田，NDVI很高，但地表温度、建筑裸土占比这些生态压力因素完全没体现出来。后来接触到徐涵秋老师提出的遥感生态指数（RSEI）模型，才真正觉得找到了一个既科学严谨、又能快速出结果的多指标评价框架。这篇文章我就把这套方法从原理到实操完整梳理一遍，给正在做生态环境评价、城市规划研究、或者毕业论文需要出生态质量图的朋友做个参考。

RSEI的核心思路其实很直接：把影响生态环境的四大要素——绿度、湿度、热度、干度——用遥感数据分别量化，再用主成分分析把它们客观地压缩成一个综合指数。整个流程在ENVI、ArcGIS或者Google Earth Engine上都能跑通，关键是理解每一步为什么要这么做。

1. RSEI模型到底在评价什么：四大分量的生态语义

在动手算之前，我们得先把RSEI的四个输入分量吃透。这决定了你后面选数据、调参数的时候心里有没有底。

1.1 绿度：植被覆盖的"健康表"

绿度指标用的是NDVI，即归一化植被指数。它的原理是利用植被在红光波段强吸收、近红外波段强反射的光谱特性，把植被覆盖信息从地表信号里分离出来。公式是：

NDVI = (NIR - Red) / (NIR + Red)

对Landsat 8/9来说，就是（Band5 - Band4）/（Band5 + Band4）。NDVI的取值范围在-1到1之间，理论上裸土接近0，茂密森林可以到0.8以上，水体为负值。

我在实际处理中最想提醒的一点是：NDVI虽然看起来简单，但不同季节的植被状态差异极大。如果是做多时期对比，影像的获取月份必须尽量一致，比如都用每年8-9月份的影像，否则你很难区分生态变化和物候变化。

1.2 湿度：土壤含水与地表湿润度

湿度分量来自缨帽变换（Tasseled Cap Transformation）的湿度分量WET。缨帽变换是一种针对植被、土壤、水体等典型地物设计的线性变换，它能把原始多波段数据压缩成几个有明确生态含义的组分：亮度、绿度和湿度。

不同传感器的WET系数差别很大，这点非常容易踩坑。Landsat 8 OLI地表反射率产品的WET系数是：

WET = 0.1511 * Blue + 0.1973 * Green + 0.3283 * Red + 0.3407 * NIR - 0.7117 * SWIR1 - 0.4559 * SWIR2

注意，这里的系数是基于地表反射率数据推导的，如果你用的是大气表观反射率（TOA）产品，系数必须换一套。另外，有些教程里给的是Landsat 5 TM的系数，直接套到Landsat 8上，结果出来全是乱的，这个坑后面我会详细讲。

1.3 热度：地表温度的生态压力

热度分量就是地表温度LST。城市热岛效应是生态质量恶化的最直接表现之一，所以在RSEI里热度被当作重要压力指标。

LST的反演方法有很多种，比如辐射传输方程法、劈窗算法、单窗算法等。实操中我们最常用的是辐射传输方程法（大气校正法），它的基本思路是：先把热红外波段的DN值转换为辐射亮度，再转换为亮温，然后结合地表比辐射率做发射率校正，最终得到真实的地表温度。

这个指标的计算链路最长，涉及到好几个中间参数，很多新手就是在这步被劝退的。别急，第3章我会把每一步的推导和代码都给你列出来。

1.4 干度：建筑与裸土的"人工印记"

干度指标NDBSI（Normalized Difference Bare Soil and Building Index）是裸土指数SI和建筑指数IBI的均值：

NDBSI = (SI + IBI) / 2

为什么要用两个指数的均值？因为单一指数很难同时兼顾自然裸土和人工建筑。比如传统的裸土指数容易把建筑也识别进来，建筑指数又容易漏掉裸土。把两者平均，能在一定程度上相互弥补，将这两类"不透水面"和"裸露地表"统一表征为生态干度。

这个指标是RSEI里最能反映人类活动干扰强度的分量。一个区域如果NDBSI高，通常意味着城市化程度高、地表硬化严重或者水土流失明显，这些都属于生态退化的重要信号。

理解了这四大分量的生态含义之后，你会发现RSEI其实构建了一个"压力-状态-响应"的简化框架：植被状态（绿度、湿度）代表生态系统的基础状况，人类活动干扰（热度、干度）代表生态压力。综合起来，就得到了一个能同时反映自然禀赋和人为影响的生态质量指数。

2. 数据源选型与预处理：决定RSEI成败的第一道关口

方向不对，努力白费。数据选型和预处理的质量，直接决定了后面所有计算结果的可靠性。

2.1 Landsat 8/9还是Sentinel-2：怎么选

目前做RSEI用得最多的是Landsat系列，主要原因是它拥有从1980年代至今的连续存档，非常适合做几十年的生态变化分析。具体到传感器，我推荐以下选型原则：

数据源	分辨率	优势	注意点
Landsat 8/9 OLI+TIRS	30m（热红外100m）	2013年至今，存档连续，有成熟的热红外波段	重访周期16天，受云影响大
Landsat 5 TM	30m（热红外120m）	1984-2011年历史数据	数据较老，辐射定标参数需仔细处理
Landsat 7 ETM+	30m	1999年至今	2003年后有条带，需要插值修复
Sentinel-2 MSI	10/20m	分辨率高，重访周期5天	没有热红外波段，LST需要单独处理

有一点必须澄清：Sentinel-2本身没有热红外传感器，直接用Sentinel-2做完整的RSEI是算不出LST的。虽然有些研究用Landsat的LST降尺度配合Sentinel-2做高分辨率RSEI，但那超出了本文基础教程范畴。我建议新手先老老实实用Landsat 8/9的Collection 2 Level 2产品。

2.2 云掩膜与影像合成

云是光学遥感最大的敌人。一片薄云或者云阴影覆盖的区域，反射率信号完全失真，如果没做掩膜就参与计算，那几个像元的NDVI、LST都会变成离谱的异常值，直接影响后续PCA的特征向量。

在Google Earth Engine里，我们可以直接用Landsat Collection 2 Level 2产品自带的QA_PIXEL波段做云掩膜，这是目前最省事也最靠谱的做法。我常用的JavaScript代码如下：

function maskL8sr(image) { // 获取QA波段 var qa = image.select('QA_PIXEL'); // 设置要掩膜的位标志：云、云阴影、冰雪 var cloudBitMask = 1 << 1; // Cloud var shadowBitMask = 1 << 3; // Cloud Shadow var snowBitMask = 1 << 4; // Snow // 创建掩膜 var mask = qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(shadowBitMask).eq(0)) .and(qa.bitwiseAnd(snowBitMask).eq(0)); // 去掉缩放系数，保留SR波段和热红外波段 return image.updateMask(mask) .select(['SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7','ST_B10'], ['Blue','Green','Red','NIR','SWIR1','SWIR2','Thermal']); }

这里要特别注意热红外波段ST_B10的定标单位。Collection 2 Level 2的地表温度产品单位是开尔文乘以10，在后续计算时需要先乘以0.00341802并加上149.0来还原成真实的开尔文温度，这个细节很坑，后面演示代码里我会带上。

如果你用的是ENVI/ArcGIS路线，那么下载数据时应优先选Landsat Collection 2 Level 2科学产品，而不是Level 1原始DN产品。Level 2已经帮你完成了大气校正，可以直接拿SR波段做指数计算，省掉了最耗时的大气校正环节。

2.3 投影、裁剪与统一分辨率

RSEI的四个分量中，NDVI、WET、NDBSI是反射率比值或线性组合，不受投影影响；但LST是热红外波段计算来的，原始分辨率是100米。如果在ENVI里直接做波段合成，需要先把热红外重采样到30米，否则最后PCA的输入图层分辨率不一致，会产生错位。

在GEE中处理非常简单，用select和resample即可：

// 重采样到30米并统一投影 var thermal = image.select('Thermal').resample('bilinear'); var sortedImage = image.addBands(thermal, null, true);

统一投影时，建议全部转到UTM投影。如果研究区跨多个UTM带，可以用Albers等积投影，保证面积计算不变形。

另外别忘了裁剪到研究区边界。用clip接研究区的FeatureCollection即可：

var roi = ee.FeatureCollection('projects/your-project/assets/roi'); var imageClipped = sortedImage.clip(roi);

这一步能显著减少后续PCA的计算量，尤其是研究区只占影像一小部分的时候。

3. 四大指标计算实操：公式、参数与GEE代码

预处理做完，就到了最核心的指标计算环节。我直接以Landsat 8/9 Collection 2 Level 2为输入，给出完整的GEE实现。

3.1 NDVI与WET的计算

NDVI前面已经给了公式，在GEE里用normalizedDifference一行就能解决：

var ndvi = image.normalizedDifference(['NIR', 'Red']).rename('NDVI');

WET的计算需要用到缨帽变换系数。这里我给出Landsat 8 OLI地表反射率对应的完整公式：

var wet = image.expression( '0.1511 * B2 + 0.1973 * B3 + 0.3283 * B4 + 0.3407 * B5 - 0.7117 * B6 - 0.4559 * B7', { 'B2': image.select('Blue'), 'B3': image.select('Green'), 'B4': image.select('Red'), 'B5': image.select('NIR'), 'B6': image.select('SWIR1'), 'B7': image.select('SWIR2') } ).rename('WET');

我把这一段单独列出来，是因为系数搞混的概率实在太高了。网上流传的缨帽变换系数有很多版本，有基于TOA反射率的、有基于地表反射率的、有TM的、有ETM+的，用错一个版本，湿度分量就可能出现大面积负数，后面PCA的结果也会跟着错。

3.2 Landsat 8地表温度（LST）反演的完整链路

LST是RSEI四个分量里计算链路最长的，我按顺序把每一步拆开讲清楚。

第一步：读取热红外波段的DN值并定标。GEE里ST_B10已经是Level 2处理过的地表温度产品，但单位是"开尔文×10"，需要还原：

var thermalDN = image.select('Thermal'); var lstKelvin = thermalDN.multiply(0.00341802).add(149.0).rename('LST_K');

如果用的是Collection 2 Level 1的原始数据，则需要先做辐射定标得到辐射亮度L，再用热红外波段的K1、K2常数换算亮温。Landsat 8 Band 10的K1=774.8853、K2=1321.0789（这是当前Collection 2版本的官方常数，不同版本有差异，务必查询影像MTL文件）：

// 如果从L1数据出发，DN转辐射亮度 var ML = 0.0003342; // 从MTL文件获取 var AL = 0.1; // 从MTL文件获取 var radiance = image.select('B10').multiply(ML).add(AL); // 辐射亮度转亮温 var brightnessTemp = radiance.expression( 'K2 / log(K1 / rad + 1)', {'rad': radiance, 'K1': 774.8853, 'K2': 1321.0789} );

第二步：计算地表比辐射率。这一步需要用NDVI估算植被覆盖度FVC，再根据植被和裸土的比辐射率加权平均：

// 植被覆盖度 var ndviMin = 0.05; var ndviMax = 0.75; var fvc = ndvi.subtract(ndviMin).divide(ndviMax - ndviMin).clamp(0, 1).rename('FVC'); // 地表比辐射率：水体取0.99，城镇和裸地取0.97，自然地表用FVC加权 var emissivity = fvc.expression( '(FVC >= 0) && (NDVI < 0.05) ? 0.995 : ' + '((NDVI >= 0.05 && NDVI <= 0.75) ? (0.004 * FVC + 0.986) : 0.99)', {'FVC': fvc, 'NDVI': ndvi} ).rename('EM');

第三步：用单窗算法反演真实地表温度。经典的简化公式是：

LST = BT / [1 + (λ * BT / ρ) * ln(ε)]

其中λ是热红外波段中心波长（Landsat 8 Band 10取10.895微米），ρ = h * c / σ ≈ 1.438 * 10^-2 m·K，ε为比辐射率。代码实现如下：

var lambda = 10.895; // 微米 var rho = 1.438e-2; // 米·开尔文 var lst = brightnessTemp .divide(brightnessTemp.multiply(lambda).divide(rho).log().multiply(-1).add(1))

等等，这个公式我写反了。正确的单窗算法简化式是：

LST = BT / (1 + (λ * BT / ρ) * ln(ε))

因为对数项是负值（ε < 1），所以分母小于1，LST会比亮温略高。GEE里正确写法是：

var lst = brightnessTemp.divide( ee.Image(1).add( brightnessTemp.multiply(lambda).divide(rho).multiply(emissivity.log()) ) ).rename('LST');

最后别忘了把开尔文转摄氏度，方便制图时标注：

var lstC = lst.subtract(273.15).rename('LST_C');

3.3 NDBSI建筑与裸土指数计算

NDBSI包含SI和IBI两个指数，先分别算再取平均。SI的公式为：

SI = [(SWIR1 + Red) - (NIR + Blue)] / [(SWIR1 + Red) + (NIR + Blue)]

IBI的公式为：

IBI = {2SWIR1/(SWIR1+NIR) - [NIR/(NIR+Red) + Green/(Green+SWIR1)]} / {2SWIR1/(SWIR1+NIR) + [NIR/(NIR+Red) + Green/(Green+SWIR1)]}

从结构上可以看出，IBI利用的是建筑在短波红外波段的反射率明显高于植被的特点，以及植被在近红外高反射、建筑在绿光波段反射率相对高的光谱差异。代码实现如下：

var si = image.expression( '((SWIR1 + Red) - (NIR + Blue)) / ((SWIR1 + Red) + (NIR + Blue))', { 'SWIR1': image.select('SWIR1'), 'Red': image.select('Red'), 'NIR': image.select('NIR'), 'Blue': image.select('Blue') } ).rename('SI'); var ibi = image.expression( '(2*SWIR1/(SWIR1+NIR) - (NIR/(NIR+Red) + Green/(Green+SWIR1))) / ' + '(2*SWIR1/(SWIR1+NIR) + (NIR/(NIR+Red) + Green/(Green+SWIR1)))', { 'SWIR1': image.select('SWIR1'), 'NIR': image.select('NIR'), 'Red': image.select('Red'), 'Green': image.select('Green') } ).rename('IBI'); var ndbsi = si.add(ibi).multiply(0.5).rename('NDBSI');

到这里，NDVI、WET、LST、NDBSI四个分量影像就全部准备好了。接下来的PCA才是RSEI的灵魂。

4. 主成分分析（PCA）与RSEI合成：从四个指标到一个指数

为什么要用PCA而不是简单地给四个指标各赋一个权重求平均？这是RSEI模型最有技术含量的地方。

4.1 为什么先做归一化

四个分量的量纲和数值范围差异巨大：NDVI大致在-1到1之间，LST可能在280到320开尔文，WET的范围因传感器而异，NDBSI也在-1到1附近。如果不归一化直接做PCA，数值范围大的LST会主导主成分的方差贡献，导致特征向量严重偏向热度分量，其余分量的信息被压制。

归一化公式是每个指标按0-1区间拉伸：

NI = (I - I_min) / (I_max - I_min)

在GEE里实现需要注意：这里的I_min和I_max应该基于研究区影像的实际像元分布取，而不是理论值。但直接取影像的最小最大值容易被异常像元干扰。我习惯用2%到98%的分位数来截断归一化，这样能有效抑制极少数云残留或水体异常值的影响：

function normalize(img) { var min = img.reduceRegion({ reducer: ee.Reducer.percentile([2]), geometry: roi, scale: 30, maxPixels: 1e10 }).values().get(0); var max = img.reduceRegion({ reducer: ee.Reducer.percentile([98]), geometry: roi, scale: 30, maxPixels: 1e10 }).values().get(0); return img.subtract(ee.Image.constant(min)).divide(ee.Image.constant(max).subtract(ee.Image.constant(min))) .clamp(0, 1); }

这里get(0)取到的值是数组形式，在GEE里要小心处理。更稳妥的做法是用reduceRegion返回的字典，按波段名分别取min和max。另外注意maxPixels要设得足够大，否则大范围研究区会报错。

4.2 PCA的输入与主成分贡献率判断

归一化完成后，将四个分量合成一个多波段影像，然后在GEE里用reduceRegion配合ee.Reducer.principalComponentAnalysis()做PCA。PCA的核心逻辑是：求四个波段构成的协方差矩阵的特征值和特征向量，按特征值大小排序，提取第一主成分PC1。

var image4bands = ndvi.addBands([wet, lst, ndbsi]).select(['NDVI','WET','LST','NDBSI']).float(); var pca = ee.Reducer.principalComponentAnalysis(4); var pcImage = image4bands.reduceNeighborhood(pca, ee.Kernel.square(1));

等等，reduceNeighborhood不是我们想要的方案。它是在滑动窗口里做PCA，计算量和结果的解释方式都不对。正确的做法应该是先在全影像范围内采样构建协方差矩阵，再对影像做线性变换：

// 先采样得到均值与协方差 var mean = image4bands.reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 30, maxPixels: 1e10 }); var centered = image4bands.select(['NDVI','WET','LST','NDBSI']) .subtract(ee.Image.constant(mean.values())); var covar = centered.toArray().reduceRegion({ reducer: ee.Reducer.covariance(), geometry: roi, scale: 30, maxPixels: 1e10 }); // 从协方差做PCA var eigen = ee.Array(covar.get('array')).eigen(); var eigenVectors = eigen.slice(1, 0, 1).slice(0, 0, 4); var pc1 = centered.toArray().matrixMultiply(eigenVectors.transpose()).arrayProject([0]) .arrayFlatten([['PC1']]);

这一段代码有几个技术细节值得展开：eigen()返回的是一个二维数组，每一行格式是[特征值, 特征向量分量...]，按特征值降序排列。所以slice(1, 0, 1)取出的是第一行后面的部分，即最大特征值对应的特征向量。这个特征向量就是PC1对应的权重，它的四个分量分别代表NDVI、WET、LST、NDBSI在PC1中的载荷方向。

拿到了PC1之后，我们要看两个关键指标：

第一是特征值贡献率。PC1对应的特征值占所有特征值之和的比例，一般要求大于80%。如果贡献率偏低，说明四个分量的信息没有很好地收敛到一个主轴上，模型解释力不够。

第二是特征向量的符号。在原始的RSEI框架中，我们期望NDVI和WET的载荷为正（生态好的指标贡献正向），LST和NDBSI的载荷为负（生态压力指标贡献负向）。如果你算出来的PC1载荷方向完全相反，即NDVI为负、LST为正，说明第一主成分表达的是"生态退化"而不是"生态质量"。这种情况下需要把RSEI反过来，用1减去PC1来校正方向，否则高值代表恶劣、低值代表良好，语义全反了。

4.3 RSEI合成、再归一化与方向校正

无论特征向量方向如何，最终都要做一个关键步骤：把PC1再归一化到0-1区间，然后判断是否需要方向修正。

// 先归一化PC1到0-1 var pc1Norm = normalize(pc1); // 判断方向并校正 // 如果NDVI、WET在特征向量中的系数之和为正，说明PC1越大生态越好，直接用pc1Norm // 如果为负，则需要用 1 - pc1Norm var rsei = pc1Norm; // 或者 rsei = ee.Image(1).subtract(pc1Norm);

这里我强烈建议你打印出特征向量来人工检查一下，别全指望自动判断。具体做法是print(eigenVectors)，在Console里查看四分量各自的载荷值。我处理过不少区域，大部分时候NDVI和WET的载荷为正，但也有例外——比如在大范围荒漠地区，PC1可能完全由LST主导，这时候就要特别小心。

合成之后的RSEI理论上取值范围是0到1，数值越高代表生态质量越好。但实际像元分布往往不均匀，可能集中在0.4到0.7之间。这时候我会做一次直方图拉伸或分位数拉伸，方便后续分级和可视化。

5. 结果分级、变化检测与制图输出

RSEI算出来是个连续变量，但研究上最常用的呈现方式是分级。分级既方便空间对比，也好做面积统计。

5.1 五级分类阈值的设定

目前学界用得最多的是等间距五级分类法，把0-1区间均分为五等份：

RSEI区间	等级	生态含义
0.0 - 0.2	差	生态质量极差，以高强度开发或严重退化区域为主
0.2 - 0.4	较差	生态质量偏差，人为干扰明显
0.4 - 0.6	中等	生态质量一般，自然与人为活动交织
0.6 - 0.8	良	生态质量较好，植被覆盖较高
0.8 - 1.0	优	生态质量优异，基本为茂密植被区

等间距分类最大的优势是可对比性。不同年份、不同研究区的RSEI结果都能用同一套标准比较。不过也有学者用自然断点法（Jenks）分级，让每个类别的内部差异最小化。我自己做多时期对比时，坚持用等间距法，因为自然断点法每年的断点不同，变化统计的结果会被分类方法本身干扰，很难说清是生态真变了还是分级变了。

GEE里的分级可以用expression或者where实现：

var rseiClassified = rsei.expression( '(rsei > 0 && rsei <= 0.2) ? 1 : ' + '(rsei <= 0.4) ? 2 : ' + '(rsei <= 0.6) ? 3 : ' + '(rsei <= 0.8) ? 4 : 5', {'rsei': rsei} );

5.2 多时期变化检测

RSEI最有价值的应用场景其实是时间序列变化分析。比如评估一个城市十年的生态环境演变，或者一个生态修复工程实施前后的效果对比。

操作思路很简单：分别计算2005年、2010年、2015年、2020年等时间节点的RSEI，然后做差值或矩阵转移分析。差值法就是后一期RSEI减前一期RSEI，正值代表生态改善，负值代表生态退化。矩阵转移分析则需要统计不同等级之间转移的面积。

GEE里做差值非常方便：

var rsei2020 = ...; // 2020年RSEI var rsei2010 = ...; // 2010年RSEI var change = rsei2020.subtract(rsei2010);

然后设置阈值判断改善区域和退化区域。但我建议差值分析用连续RSEI值做，分类转移矩阵则用分级后的结果做。两条技术路线互补，一个评估"程度"，一个评估"类别的转变轨迹"。

5.3 制图输出的规范

RSEI制图有几个容易忽略但影响观感的地方：

第一是配色方案。生态质量从差到优，我推荐使用从红色经黄色、绿色到深绿的连续色带。具体来说，差用暗红（#8B0000），较差用橙红（#FF4500），中等用黄色（#FFD700），良用浅绿（#32CD32），优用深绿（#006400）。这套配色有很直观的生态认知感，审稿人也习惯这种表达。千万别用蓝色当"优"，因为蓝色通常被读者理解为水体。

第二是统计信息。图上一定要附上各等级面积和百分比统计表。很多人在ArcGIS里割了图就走，评委一问"中等生态区占比多少"就答不上来，这就尴尬了。

第三是空间位置信息。研究区示意图、行政区边界、比例尺、指北针一个都不能少，这是学术制图的基本盘。

GEE里用Export.image.toDrive导出整幅栅格，再拿到ArcGIS或QGIS里做正式制图。导出时注意scale设置成30米对应Landsat分辨率，crs设置成UTM投影，maxPixels适当调大。

6. 我踩过的坑和给你的建议

RSEI的坑说多不多，但每一个都能让结果偏得离谱。我把这几年实际处理中踩过、以及看别人踩过的坑集中说一下。

6.1 坑一：归一化时被0除

这是最常见也最隐蔽的问题。归一化的分母是max - min，如果某个波段的max和min相等（比如研究区内全是均匀地物，或者掩膜后有效像元极少），分母为0，结果全面变成NoData或者Infinity。

解决办法有两个：一是归一化之前先做掩膜统计，确保有效像元数量充足；二是在公式里加一个极小值偏移量防止除零：

// 加一个极小值防止除零 var range = max.subtract(min).add(0.0001); var normalized = img.subtract(ee.Image.constant(min)).divide(ee.Image.constant(range));

6.2 坑二：PCA主成分被单一指标主导

有次我帮人复查数据，发现他的PC1特征向量里LST载荷达到0.9以上，其他三个分量几乎可以忽略。查了一圈原因：一是他直接用原始LST没有归一化，把整个分析带偏了；二是他的研究区有大量火烧迹地，LST异常高，而且没做掩膜。后来我把归一化修好、火烧迹地掩掉，PC1的贡献结构就正常了。

所以每次算出特征向量后，我都建议顺手把载荷打印出来看一眼，四个分量载荷量级应该相差不多，如果某一个分量载荷超过0.8，基本可以断定预处理出了问题。

6.3 坑三：大水体干扰严重

LST在大水面上会异常偏低，NDVI在水体上是负值，WET在水体上反而极高。如果不处理，整个PCA会被水体"带节奏"，尤其是南方湖泊密集的研究区，这个问题特别突出。

常规做法是先用水体指数（如MNDWI）把水体掩膜掉再算RSEI，最后制图时把水体单独叠加回来。这样既避免了水体对统计的干扰，图面上也能看到水体的空间分布。公式是：

var mndwi = image.normalizedDifference(['Green', 'SWIR1']).rename('MNDWI'); var waterMask = mndwi.gt(0.2).not(); // MNDWI > 0.2 视为水体，掩膜掉

6.4 经验和流程建议

最后总结一套我在多个项目里验证过的稳定流程，你可以直接照抄：

1. 明确研究区和时间范围，确定影像年份和季节（夏季晴天优先）。
2. 在GEE里批量筛选目标年份、云量低于10%的Landsat 8/9影像，做云掩膜后取中值合成。
3. 计算NDVI、WET、LST、NDBSI四个分量，并做水体掩膜。
4. 每个分量分别做2%-98%分位数截断归一化。
5. 合成后做PCA，检查PC1贡献率（需大于80%）和特征向量方向。
6. 对PC1归一化得RSEI，按需做方向校正。
7. 五级分类、统计各等级面积占比、导出栅格。
8. 在ArcGIS/QGIS里出图，配上等级统计表和比例尺。

做RSEI这个模型，公式本身不复杂，难点全在预处理和参数细节上。如果你能在自己的研究区把上面这八步完整跑通一遍，你会发现它不仅是一个评价工具，更是一套理解"自然-人类耦合系统"的思维方式。后续如果你想扩展，还可以把RSEI和人口密度、GDP等社会经济数据结合做相关性分析，或者用它来筛选生态修复的优先区域，这都是一旦用熟了RSEI之后很自然的延伸方向。