把2000年和2019年两张郑州的Landsat影像并排放在一起，不用算任何指数，肉眼就能看出这座城市的骨架拉大了多少。但“看起来变差了”和“到底差在哪、差了多少”是两回事，而要回答后者，就需要一个能同时反映绿度、湿度、热度、干度的综合指标，这就是RSEI（遥感生态指数）。这篇文章我用自己的实际处理流程，完整拆解如何在ENVI里基于Landsat影像构建郑州市2000-2019年的RSEI，从数据选择、预处理、四指标计算、主成分分析到变化检测，把每一步的关键参数和底层逻辑都讲清楚。适合正在做长时序生态遥感评价、毕业论文涉及城市生态方向、或者刚接触RSEI想少走弯路的朋友参考。
1. 为什么用RSEI盯一座城市，而不是只看NDVI或绿地率
1.1 RSEI在算什么：四个指标一个主成分
很多刚接触生态遥感的人会习惯性用NDVI代表生态状况，这当然有道理，但不够全面。城市生态系统受到的影响是多维度的：绿地减少了生态变差，但建筑和裸土增加、地表温度上升、湿地被填埋同样在恶化生态。只拿绿度一个维度说话，很容易得出“这城市绿化挺好所以生态没问题”的片面结论。
RSEI的思路很直接，由徐涵秋提出，用主成分分析把四个指标压缩成一个综合指数，这四个指标分别是：
- 绿度指标：NDVI（归一化植被指数）
- 湿度指标：WET（缨帽变换湿度分量）
- 干度指标：NDBSI（归一化裸土与建筑指数）
- 热度指标：LST（地表温度）
这四个指标分别对应城市生态的主要压力来源和改善来源。然后把它们做归一化，再用主成分分析取第一主成分PC1。PC1能解释原始数据里大部分信息，而且会自动根据数据本身的分布给四个指标分配权重，不需要人为打分，这正是RSEI相比其他综合评价方法更客观的地方。
1.2 郑州这个案例的典型性
选择郑州做长时序RSEI，有它天然的代表性。郑州是典型的平原城市，处于华北平原腹地，2000年到2019年这20年间经历了极其剧烈的城市扩张。郑东新区从一片鱼塘洼地变成建成区，三环、四环不断外扩，常西湖片区、航空港区相继开发，这种大规模的土地利用变化正好是RSEI最擅长捕捉的对象。
同时郑州也不是单一的城市化负面样本。西部有嵩山余脉的山地丘陵，北部有黄河滩区湿地，中心城区有大面积公园绿地和水系，这些生态本底让RSEI的空间异质性很明显，做出来的分级图不会是一片平铺的色块，而是有层次、有对比、能讲出故事的。用这样一座城市做研究区，RSEI的变化分析结果会更有说服力，也更容易在论文里落地。
2. 2000-2019年Landsat数据怎么选、去哪里拿
2.1 轨道号、时间窗口和云量要求
做长时序RSEI，数据源基本就是Landsat系列。郑州所在区域的WRS-2轨道号是124/036，这个固定不变，每年每季都能拿到影像。
时间窗口我建议统一选在6月到9月之间。原因有三个：第一，这个时期植被生长旺盛，NDVI对生态差异的区分度最高；第二，夏季地表温度空间差异明显，城市热岛效应突出，LST指标更有辨识度；第三，尽量保持20年间物候一致，否则你用5月的影像和9月的影像比NDVI，差值里有很大一部分是物候差异而不是真实生态变化。
云量是另一个硬性筛选条件。我一般要求整景影像云量小于10%，同时还要看一下研究区范围内是否有云。有时候整景云量很小，但恰好郑州市上空飘一朵云，那就不能用。USGS EarthExplorer和中国地理空间数据云都能在线预览缩略图，下载前一定要肉眼确认。
2.2 Landsat 5/7/8的差异与预处理链路
2000年到2019年跨越了三个传感器：Landsat 5 TM（2011年退役）、Landsat 7 ETM+（2003年后有条带故障）、Landsat 8 OLI/TIRS（2013年发射）。你自己跑的时候，年份和传感器的对应关系大概是：
| 年份 | 可选传感器 | 注意点 |
|---|---|---|
| 2000 | Landsat 5 TM | 数据质量稳定，大气校正参数好定 |
| 2005 | Landsat 5 TM | 同上 |
| 2010 | Landsat 5 TM或Landsat 7 | 如果用Landsat 7，要处理SLC条带 |
| 2015 | Landsat 8 OLI/TIRS | 波段编号与TM不同，注意转换 |
| 2019 | Landsat 8 OLI/TIRS | 数据最新，质量最好 |
预处理链路方面，如果下载的是Level-1产品，完整流程是：辐射定标 → FLAASH大气校正 → 研究区裁剪。如果直接下载Landsat Collection 2 Level-2表面反射率产品，那多光谱波段就不需要再做大气校正了，能省不少事。但热红外波段反演LST时，Level-2产品也有对应的地表温度产品（ST_B10），可以直接用，也可以自己用辐射传输方程反演，后面细说。
2.3 裁剪与坐标系统一
郑州市行政边界shapefile是必备材料。在ENVI里用Resize Data配合shp做裁剪，或者用Subset Data from ROIs按行政区范围裁剪，都可以。我的习惯是先把影像投影到WGS-84 UTM 50N，再做裁剪，这样后续算面积的时候单位是平方米，统计各生态等级面积时非常方便。
这里有个容易忽略的细节：20年间郑州市的行政区划发生过调整，早期下载的shp范围可能和现在的市界不一样。如果你要严格按行政区范围统计，建议统一使用同一版本的最新市界，或者用研究区统一范围，避免因边界变化导致面积统计对不上。
3. ENVI里跑四指标：绿度、湿度、干度、热度的完整计算过程
3.1 辐射定标和FLAASH大气校正（含Landsat 8的参数坑）
虽然Level-2表面反射率产品越来越普及，但很多情况下你手上只有Level-1数据，所以辐射定标和大气校正这个流程还是得会。
在ENVI工具箱里打开Radiometric Calibration，选择多光谱波段，输出类型选Float，定标类型选Reflectance。这里要注意，Landsat 8需要把多光谱和全色波段分开处理，全色波段不用参与后续计算。
大气校正强烈建议用FLAASH模块。操作上确实有几个坑：
第一，FLAASH要求输入辐射定标后的辐射亮度数据（单位是μW/(cm²·sr·nm)），不是反射率，所以定标时Output Interleave要选BIL或者BIP，单位要选对。
第二，传感器类型要选对。Landsat 5 TM、Landsat 7 ETM+和Landsat 8 OLI各自对应不同的参数文件。Landsat 8在FLAASH里需要手动输入中心波长和FWHM，这两个参数可以从MTL元数据文件里找到，或者在ENVI里用Metadata Viewer查看。
第三，大气模型和气溶胶模型的选择要结合时相和纬度。郑州属于中纬度地区，6-9月一般选Tropical或者Mid-Latitude Summer，气溶胶模型选Urban，初始能见度设个40km左右，具体可以根据影像实际情况调整。
如果实在不想折腾FLAASH，用QUAC快速大气校正也能凑合，但这个属于简化方案，精度上不如FLAASH，尤其在LST反演时会放大误差。能跑FLAASH就别偷懒。
3.2 NDVI：最没悬念但最关键的一步
NDVI的公式很简单：(NIR - Red) / (NIR + Red)。在ENVI Band Math里输入表达式，选对波段就行。注意Landsat 5 TM的近红外是Band 4，红光是Band 3；而Landsat 8 OLI的近红外是Band 5，红光是Band 4，波段编号整整差了一号。
NDVI虽是四指标里最简单的，但它的作用不只是当绿度指标。后面算LST的时候，植被覆盖度要用NDVI反推，所以NDVI算完之后不要删，保留下来备用。
还有一个经验：NDVI影像里如果有云或者水体，会出现异常低值。水体通常在-0.5以下，云在0.8以上，后续做归一化和PCA时会干扰结果。如果你不想单独做水体掩膜，至少要在算植被覆盖度时把NDVI小于0的像元设定为0，把大于0.8的像元设为0.8，做一次截断处理。
3.3 缨帽变换WET：注意不同传感器的系数表
WET分量需要用到缨帽变换系数。很多人在这一步翻车，因为Landsat 5和Landsat 8的系数表完全不同，不能混用。我列一下常用的：
Landsat 5 TM的WET系数（基于TM波段1-5、7）： WET = 0.0315×B1 + 0.2021×B2 + 0.3102×B3 + 0.1594×B4 - 0.6806×B5 - 0.6109×B7
Landsat 8 OLI的WET系数（基于OLI波段2-7）： WET = 0.1511×B2 + 0.1973×B3 + 0.3283×B4 + 0.3407×B5 - 0.7117×B6 - 0.4559×B7
用Band Math输入时，直接把对应波段乘上系数再求和。这里需要注意Landsat 8的B2是蓝波段、B3是绿波段、B4是红波段、B5是近红外、B6是短波红外1、B7是短波红外2，和TM时代的编号习惯完全不同。
WET分量的物理含义是土壤和植被的湿度状况，城市建成区、裸土、道路的WET偏低，水体、植被茂密区的WET偏高。这也是RSEI里体现“湿度生态效应”的指标。
3.4 NDBSI：裸土和建筑指数的叠加
干度指标NDBSI是由裸土指数SI和建筑指数IBI合成的，公式是：
NDBSI = (SI + IBI) / 2
SI（裸土指数）的公式： SI = ((B5 + B1) - (B4 + B3)) / ((B5 + B1) + (B4 + B3))
IBI（建筑指数）的公式： IBI = (2×B5/(B5+B4) - (B4/(B4+B3) + B2/(B2+B5))) / (2×B5/(B5+B4) + (B4/(B4+B3) + B2/(B2+B5)))
注意这里波段编号同样是Landsat 5和Landsat 8有区别。Landsat 5的B1蓝、B3红、B4近红外、B5短波红外；Landsat 8的B2蓝、B4红、B5近红外、B6短波红外1。
我建议分步算，先把SI和IBI分别算出来，再用Band Math做平均，这样中间可以检查每步结果的数值范围是否合理。SI和IBI的值域不会严格落在[-1,1]内，会出现一些超界像元，这是正常的，归一化之后会被拉回0-1区间。
这里有个实际体会：在郑州这种平原城市里，IBI对高楼密集区的响应很灵敏，但黄河水面、大型水体的IBI会出现明显负值，如果不做任何处理，后续归一化时水体会把干度指标的分布拉得很偏。所以更稳妥的做法是，如果研究区内有较大面积水体，建议先用NDVI阈值或水体指数做一个简单掩膜，把水域剔除后，再进入PCA。后面我会单独讲水体的取舍。
3.5 地表温度LST：大气校正法反演
地表温度反演是RSEI四指标里最绕的一步，但绕也得做，因为热度指标是城市生态评价的核心。
我用的是大气校正法，也叫辐射传输方程法，分三步。
第一步，把热红外波段辐射定标后的数据换算成亮度温度。ENVI里有现成的工具，也可以直接用辐射定标后得到的辐亮度数据，用Planck反函数手动算。Landsat 5 TM的热红外是Band 6，Landsat 8是Band 10（TIRS），不要搞混。
第二步，计算地表比辐射率ε。徐涵秋在RSEI论文里推荐了简化方法：先用NDVI计算植被覆盖度FVC：
FVC = (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)
其中NDVI_soil一般取0.05，NDVI_veg取0.70。然后比辐射率： ε = 0.004×FVC + 0.986
这个简化公式对混合像元有不错的精度，在ENVI里用Band Math一条表达式就能算出来。
第三步，用大气校正法反演LST。公式是：
LST = T / (1 + (λ×T/ρ)×ln(ε))
其中T是亮度温度（K），λ是热红外波段中心波长（TM是11.45μm，OLI Band 10是10.9μm），ρ=h×c/σ，约等于1.438×10⁻² m·K。
大气校正法需要大气剖面参数，如果手头没有同步气象数据，可以用NASA提供的在线大气参数查询工具获取，或者使用ENVI扩展工具中的大气校正参数估算工具。我实际操作下来，工具的估算精度足够RSEI这种相对比较分析用，结果以摄氏度或开尔文为单位都可以，因为后面归一化会抹掉单位差异。
不过这里必须强调一个细节：Landsat 8热红外波段原始分辨率是100米，和30米的多光谱波段不一致。在做LST之前，要把热红外波段用三次卷积法重采样到30米，否则后面PCA的时候影像尺寸对不上，会很麻烦。
3.6 四指标归一化的顺序问题
四个指标计算完之后，先做归一化再做PCA。归一化公式是：
NI = (I - I_min) / (I_max - I_min)
在ENVI里可以直接用Statistics查看每幅影像的Max和Min，然后把值带入Band Math。这里有三种做法：一是用全局最大最小值归一化，二是用直方图2%截断后的值归一化，三是用每个波段各自的最大最小值归一化。
我用的是第三种，也就是各期影像各自归一化。为什么？因为长时序分析里，如果2000年影像整体偏暗、NDVI偏低，用全局统计会把低值像元的差异拉平，不利于反映真实的空间梯度。每期影像各自归一化虽然会导致不同年份之间RSEI的绝对数值不完全可比，但它的等级划分和空间分布更符合该年份自身的生态特点。实际上多数RSEI论文也是这么做的。
如果你要用多期影像的最小最大值统一归一化，那就要在ENVI里先做影像镶嵌或者层叠分析，把多期数据放在一起统计算，操作更复杂，但跨年份的可比性会更好。两种方案没有绝对的对错，关键是在写论文时交代清楚。
4. 主成分分析合成RSEI：PC1方向、贡献率与结果翻转
4.1 PCA的操作步骤
四个归一化指标准备好之后，在ENVI工具箱里用Principal Components工具，也可以直接用Layer Stacking把NDVI、WET、NDBSI、LST按顺序合成一个多波段文件，然后做主成分分析。
这里有一个选择：用协方差矩阵还是相关系数矩阵。RSEI的标准做法是用协方差矩阵，因为四个指标都已经归一化到0-1，量纲一致了，协方差矩阵能保留各指标本身的方差信息，而不是把方差标准化。在ENVI的PCA对话框里选Covariance Matrix即可。
输出的PC1会是一个单波段影像，后台会生成一份特征值和特征向量的报告，这个报告很重要，一定要保存下来，里面包含了PC1的贡献率（Eigenvalue percentage）和四个指标在PC1上的载荷（Eigenvector）。
我跑郑州数据的经验是，PC1的贡献率普遍在65%-85%之间。你不需要过分纠结贡献率是不是一定要到90%，RSEI的标准做法就是只取PC1。如果PC1贡献率连60%都不到，那要回头检查四指标的计算和归一化是不是有问题。
4.2 载荷方向判断：什么时候RSEI0=1-PC1
这是RSEI流程里最容易被忽略的一步，直接影响最终结果的正负向。
PCA本身只是一种线性变换，它不关心“生态好”是正方向还是负方向。因此跑出来的PC1可能呈现出两种方向：一种是指标中绿度和湿度高的像元PC1值较大；另一种是反过来，干度和热度高的像元PC1值较大。
怎么判断？看PCA输出的特征向量表。如果PC1对应的特征向量里，NDVI和WET的载荷为正，NDBSI和LST的载荷为负，说明PC1越大生态越好，可以直接用RSEI0 = PC1。但如果NDVI和WET是负的，NDBSI和LST是正的，说明PC1越大生态越差，这时候就要翻转，用RSEI0 = 1 - PC1。
郑州近20年的数据里，我碰到的情况大都是PC1方向符合第二种模式（干度、热度为正，绿度、湿度为负），所以都用RSEI0 = 1 - PC1。但你做每一期影像都得单独检查方向，不能想当然用固定公式。
翻转之后再做一次0-1归一化，得到最终的RSEI值。RSEI越接近1，生态状况越好；越接近0，生态状况越差。
4.3 分级阈值怎么定，生态等级怎么解释
得到0-1连续值之后，为了分析和制图，一般会分成等级。最常用的是等间距分级，也就是：
| 等级 | RSEI范围 | 生态含义 |
|---|---|---|
| 差 | 0.0 - 0.2 | 生态质量差，以高密度建成区、裸地为主 |
| 较差 | 0.2 - 0.4 | 生态质量较差，以中低密度建成区为主 |
| 中等 | 0.4 - 0.6 | 生态质量中等，城郊过渡带、农业区 |
| 良 | 0.6 - 0.8 | 生态质量良好，以农田、林草覆盖区为主 |
| 优 | 0.8 - 1.0 | 生态质量优，以密林、湿地、水体为主 |
有人会用NDVI的等分间隔在RSEI上找断点做非等间距分级，但我在实际项目中觉得等间距分类的可读性和可解释性最好，后续做转移矩阵也方便。ENVI里用Density Slice做分级，导出RGB色带文件后可以用于制图。
5. 五期/六期结果的时间序列对比：变化检测与结果解读
5.1 差值法与等级转移矩阵
做完每一期的RSEI分级后，接下来是长时序研究里最有价值的部分：变化分析。我一般做两种分析。
第一种是差值法：把2019年的RSEI减去2000年的RSEI，做RSEI变化检测图。正差值表示生态变好，负差值表示生态变差。在ENVI里用Band Math两幅影像相减就行，然后设定变化阈值（比如±0.05以内算稳定，超过算变好/变差）。
第二种是等级转移矩阵：在ENVI里做分类后比较（Confusion Matrix Using Ground Truth Image的思路类似，或者用Change Detection Statistics工具），统计2000年各等级到2019年各等级的转移百分比。这张转移矩阵能非常直观地反映生态退化的主要路径：比如“良”变成了“中等”、“中等”变成了“较差”，这些转移方向对生态修复决策很有参考价值。
如果你用的是每期影像各自归一化再分级的结果，做等级转移矩阵没问题；但如果做RSEI数值相减的差值分析，建议所有年份都用同一套最小最大值归一化，否则数值差的语义会打折扣。
5.2 郑州生态格局的实际变化
我用这个流程跑完2000年、2005年、2010年、2015年和2019年五期数据之后，空间分布上能看到几个很明显的变化。
中心城区的RSEI等级这些年在核心位置始终偏低，但低值区范围一直在外扩，尤其是郑东新区CBD、高铁站周边、航空港区这些新建区域，生态等级从“中等”逐步跌成“较差”或“差”，这是城市扩张的直接代价。
西部沿嵩山余脉的山地丘陵区域，因为植被覆盖好、人口密度低，RSEI始终稳定在“良”和“优”的水平，这20年里波动很小，说明山体生态本底是郑州抵抗生态退化的压舱石。
北部黄河沿岸的滩区比较有意思，有些年份RSEI上升、有些年份下降，和当年黄河水量、滩区耕地利用方式有关。黄河湿地保护政策落地后，2015年到2019年这段时间，沿岸不少像元的生态等级明显好转。
整体来看，郑州2000-2019年的RSEI呈现出“中心恶化、边缘改善、总体略有下降但局部修复”的格局。这种细节丰富的研究结论，正是长时序RSEI相比单一时期生态评价最大的优势。
5.3 制图与输出建议
做生态评价，图比表格更重要。我的制图习惯是：RSEI分级图用从红到绿的渐变色（红=差，黄=中等，绿=优），加上郑州市界和主要水系，输出分辨率300dpi以上的GeoTIFF，最后在制图软件里加比例尺、指北针和图例。
ENVI里可以直接在Layer Manager里对RSEI分级图层做Density Slice配色，然后导出为矢量或直接出图。如果论文里需要各等级面积表，可以在ENVI里用Class Statistics统计每级像元数，乘以单个像元面积（30m×30m=900平方米），换算成平方公里。
另一个实用技巧：把五期RSEI影像用同一个色带和同一套分级阈值制图，排版在一个图版里，读者一眼就能看出空间格局的变化轨迹。这也是RSEI论文里的标准展示方式。
6. 长时序RSEI最容易翻车的几个坑
6.1 传感器系数不一致：Landsat 5和Landsat 8的缨帽系数不能混用
这是长时序RSEI里最隐蔽、最坑人的一个错误。因为TM和OLI的光谱响应范围不一样，WET系数表完全不是一套，你要是图省事复制粘贴同一套系数，算出来的WET在2000年和2019年之间会存在系统性的差异，导致主成分分析结果里出现虚假的变化趋势。
我自己的习惯是，在做数据预处理之前先建一个表格，把每个年份的传感器类型、波段编号对应关系、WET系数、SI/IBI的波段映射全部列清楚。这个表看起来费时间，但它能避免你在算到第三个年份时突然发现自己用错了波段。具体系数前面已经列出，这里不再重复。
同样的问题也存在NDBSI的波段引用上。Landsat 5的B1对应Landsat 8的B2，搞错一个波段，整个干度指标就废了。
6.2 Landsat 7条带和Landsat 5老数据的大气校正问题
如果你研究时期跨越了2003年到2012年，很可能会用到Landsat 7 ETM+的数据。2003年SLC故障之后，每景影像都有楔形条带缺失，大约损失22%的像元。在ENVI里可以用Landsat Gapfill工具做条带修复，用同区域前后时相的影像填充。但要注意，2010年前后如果可以用Landsat 5影像，优先用Landsat 5，毕竟条带修复是插值，不是原始观测。
Landsat 5的老数据也有坑。2000年左右的TM影像噪声水平偏高，FLAASH做完之后局部区域会出现条带状伪影。我碰到过右下角研究区内出现横向条纹的情况，后来发现是L0数据本身的扫描行校正问题，重新下载另一景过境影像才解决。遇到异常影像不要硬着头皮处理，直接换相邻日期、相同轨道号的影像是最省事的选择。
6.3 水体、负值和异质性像元的处理
RSEI的四指标里，水体在NDVI上是负值、在WET上是高值、在NDBSI上是负值、在LST上是低值，它在主成分分析里会变成一个极端的离群群体，把PC1的空间分布拉向水体主导的方向。
很多RSEI论文说不做水体掩膜，理由是RSEI本身对水体也有响应。但以我做郑州数据的经验来看，郑州北部有黄河、中部有贾鲁河、城区散布着不少人工湖，如果不做掩膜，干度指标NDBSI的最小值会被水体拉得很低，归一化后建成区和裸地的干度区分度会受到压制。
建议的处理方式是：计算一个简单的NDWI=（Green - NIR）/（Green + NIR），NDWI大于0的像元视为水体，在进入PCA之前用ENVI的Build Mask功能掩膜掉。如果你不想掩膜，那就至少要在论文方法部分说明你保留水体的理由，并且检查PC1载荷方向是否因为水体而发生了反转。
另外，归一化时如果影像里出现了异常负值或极端高值（比如云的边界像元），最大最小值会被污染。我通常先用直方图查看各指标的前后2%截断值，如果异常值非常明显，直接在Band Math里做一个条件判断，把小于0的设0、大于1的设1，再做归一化。虽然粗暴，但对PCA的稳健性很有帮助。
6.4 物候一致性和结果校验
另一个经常被忽略的问题，是不同年份影像的获取日期不完全一致。2000年你是8月15号的影像，2005年是7月3号，2010年是9月10号，虽然都在生长季，但作物物候和植被状态有差异，这会直接干扰NDVI和RSEI的年际对比。
碰到这种情况，我的做法是优先选择同旬或同月的影像，如果实在找不到，就在结果解释时明确指出物候差异可能带来的误差范围。严格的论文里可以考虑用像元二分模型或者谐波函数做物候校正，但对于大多数工程和研究场景，控制数据源在相邻月份的难度已经足够。
最后，结果输出之后一定要做校验。拿RSEI分级结果和高分影像或Google Earth历史影像做点位对比，在高德地图历史影像上随机选20-30个点，人工判读生态等级和RSEI等级是否一致，误差超过两个等级的样本要回头查原因。这个环节看起来费功夫，但它能帮你发现很多数据预处理的隐蔽问题。
我自己在项目里养成的一个习惯是：永远保留每一期的中间结果（辐射定标、大气校正、四指标归一化、PCA载荷表），这样如果后面发现某个环节出了偏差，不需要整条流水线推倒重来，只需要替换出问题的中间产物，重新跑下游步骤即可。20年的数据量不算大，但5期全流程重跑的代价也不小，保留中间结果是对自己时间和耐心的最大尊重。
长时序RSEI这个项目做到后面，你会发现自己对城市生态的理解会因为数据而变得具体：哪里在变差，哪里在恢复，哪里虽然建成区在扩张但绿地跟上了，这些都能在指数曲线上得到验证。把2000年和2019年两张郑州的Landsat影像并排摆在一起，不用算任何指数，肉眼就能看出这座城市的骨架拉大了多少。但“看起来变差了”和“到底差在哪、差了多少”是两回事，而要回答后者，就需要一个能同时反映绿度、湿度、热度、干度的综合指标，这就是RSEI（遥感生态指数）。这篇文章我用自己的实际处理流程，完整拆解如何在ENVI里基于Landsat影像构建郑州市2000-2019年的RSEI，从数据选择、预处理、四指标计算、主成分分析到变化检测，把每一步的关键参数和底层逻辑都讲清楚。适合正在做长时序生态遥感评价、毕业论文涉及城市生态方向、或者刚接触RSEI想少走弯路的朋友参考。
1. RSEI到底在算什么东西，以及为什么选郑州
1.1 RSEI的原理：四个指标压缩成一个指数
很多刚接触生态遥感的人会习惯性用NDVI代表生态状况，这当然有道理，但不够全面。城市生态系统受到的影响是多维度的：绿地减少了生态变差，但建筑和裸土增加、地表温度上升、湿地被填埋同样在恶化生态。只拿绿度一个维度说话，很容易得出“这城市绿化挺好所以生态没问题”的片面结论。
RSEI的思路很直接，由徐涵秋提出，用主成分分析把四个指标压缩成一个综合指数，这四个指标分别是：
- 绿度指标：NDVI（归一化植被指数）
- 湿度指标：WET（缨帽变换湿度分量）
- 干度指标：NDBSI（归一化裸土与建筑指数）
- 热度指标：LST（地表温度）
这四个指标分别对应城市生态的主要压力来源和改善来源。然后把它们做归一化，再用主成分分析取第一主成分PC1，作为生态质量的综合评价指标。PC1的贡献率越高，说明这四个指标对生态的联合解释能力越强，RSEI也就越可靠。
不人为设置权重是RSEI最大的优点之一。传统生态评价往往靠层次分析法或者专家打分，把绿度权重设多少、温度权重设多少，主观性很强。而RSEI让数据自己说话，各指标的权重由PC1的特征向量决定，这既简化了流程，也提高了研究结果的可重复性。
1.2 郑州作为研究区的独特价值
选择郑州做长时序RSEI，有它天然的代表性。郑州是典型的平原城市，处于华北平原腹地，2000年到2019年这20年间经历了极其剧烈的城市扩张。郑东新区从一片鱼塘洼地变成建成区，三环、四环不断外扩，常西湖片区、航空港区相继开发，这种大规模的土地利用变化正好是RSEI最擅长捕捉的对象。
同时郑州也不是单一的城市化负面样本。西部有嵩山余脉的山地丘陵，北部有黄河滩区湿地，中心城区有大面积公园绿地和水系，这些生态本底让RSEI的空间异质性很明显，做出来的分级图不会是一片平铺的色块，而是有层次、有对比、能讲出故事的。用这样一座城市做研究区，RSEI的变化分析结果会更有说服力，也更容易在论文里落地。
2. 2000-2019年Landsat数据怎么选、去哪里拿
2.1 轨道号、时间窗口和云量要求
做长时序RSEI，数据源基本就是Landsat系列。郑州所在区域的WRS-2轨道号是124/036，这个固定不变，每年每季都能拿到影像。
时间窗口我建议统一选在6月到9月之间。原因有三个：第一，这个时期植被生长旺盛，NDVI对生态差异的区分度最高；第二，夏季地表温度空间差异明显，城市热岛效应突出，LST指标更有辨识度；第三，尽量保持20年间物候一致，否则你用5月的影像和9月的影像比NDVI，差值里有很大一部分是物候差异而不是真实生态变化。
云量是另一个硬性筛选条件。我一般要求整景影像云量小于10%，同时还要看一下研究区范围内是否有云。有时候整景云量很小，但恰好郑州市上空飘一朵云，那就不能用。在USGS EarthExplorer和国内的地理空间数据云平台都能在线预览缩略图，下载前一定要肉眼确认研究区不超标。
2.2 Landsat 5/7/8的差异与预处理链路
2000年到2019年跨越了三个传感器：Landsat 5 TM（2011年退役）、Landsat 7 ETM+（2003年后有条带故障）、Landsat 8 OLI/TIRS（2013年发射）。年份和传感器的对应需要注意：
| 年份 | 可选传感器 | 注意点 |
|---|---|---|
| 2000 | Landsat 5 TM | 数据质量稳定，大气校正参数好定 |
| 2005 | Landsat 5 TM | 同上 |
| 2010 | Landsat 5 TM或Landsat 7 | 如果用Landsat 7，要处理SLC条带 |
| 2015 | Landsat 8 OLI/TIRS | 波段编号与TM不同，注意转换 |
| 2019 | Landsat 8 OLI/TIRS | 数据较新，质量稳定 |
如果下载的是Level-1产品，完整预处理链路是：辐射定标 → FLAASH大气校正 → 研究区裁剪。如果直接下载Landsat Collection 2 Level-2表面反射率产品，那多光谱波段就不需要再做大气校正了，能省不少事。但热红外波段反演LST时，建议还是从Level-1开始自己走一遍，或者使用Level-2自带的地表温度产品ST_B10。
2.3 裁剪与坐标系统一
郑州市行政边界shapefile是必备材料。在ENVI里用Resize Data配合shp做裁剪，或者用Subset Data from ROIs按行政区范围裁剪，都可以。我的习惯是先把整景影像投影到WGS-84 UTM 50N，再做裁剪，这样后面算面积的时候单位是平方米而不是度，统计各生态等级面积时非常方便。
这里有个容易忽略的细节：20年间郑州的行政区划发生过调整。如果你要严格按行政区范围统计，建议统一使用同一版本的最新市界，或者直接用一个固定的矩形研究区范围，避免因边界变化导致面积统计对不上。我采用的方案是下载一份郑州全市域的最新shp，所有年份统一用这个范围裁剪，保证可比性。
3. ENVI里跑四指标：绿度、湿度、干度、热度的完整计算过程
3.1 辐射定标和FLAASH大气校正（含Landsat 8的参数坑）
虽然Level-2表面反射率产品越来越普及，但很多情况下你手上只有Level-1数据，所以辐射定标和大气校正这个流程还是得会。
在ENVI工具箱里打开Radiometric Calibration，选择多光谱波段，输出类型选Float，定标类型选Reflectance。这里要注意，Landsat 8需要把多光谱和全色波段分开处理，全色波段不用参与后续计算。
大气校正强烈建议用FLAASH模块。有几个坑需要注意：
第一，FLAASH要求输入辐射定标后的辐射亮度数据，单位一般是μW/(cm²·sr·nm)，不是反射率。所以定标时输出类型要选Radiance而不是Reflectance，否则FLAASH会报错或者结果异常。
第二，传感器类型要选对。Landsat 5 TM、Landsat 7 ETM+和Landsat 8 OLI各自对应不同的参数文件。Landsat 8在FLAASH里需要手动输入中心波长和FWHM，这两个参数可以从MTL元数据文件里找。
第三，大气模型和气溶胶模型的选择要结合时相和纬度。郑州属于中纬度地区，6-9月一般选Tropical或者Mid-Latitude Summer，气溶胶模型选Urban，初始能见度设个40km左右，具体可以根据影像实际情况调整。
如果实在不想折腾FLAASH，用QUAC快速大气校正也能凑合，但这个属于简化方案，精度上不如FLAASH，尤其在LST反演时会放大误差。能跑FLAASH就别偷懒。
3.2 NDVI：最没悬念但最关键的一步
NDVI的公式很简单：(NIR - Red) / (NIR + Red)。在ENVI Band Math里输入表达式，选对波段就行。注意Landsat 5 TM的近红外是Band 4，红光是Band 3；而Landsat 8 OLI的近红外是Band 5，红光是Band 4，波段编号整整差了一号。
NDVI虽是四个指标里最简单的，但它的作用不只是当绿度指标。后面算LST的时候，植被覆盖度要用NDVI反推，所以NDVI算完之后不要删，保留下来备用。
还有一个经验：NDVI影像里如果有云或者水体，会出现异常低值或高值。水体通常在-0.5以下，云在0.8以上，后续做归一化和PCA时会干扰结果。如果你不想单独做水体掩膜，至少要在算植被覆盖度时把NDVI小于0的像元设定为0，把大于0.8的像元设为0.8，做一次截断处理。
3.3 缨帽变换WET：注意不同传感器的系数表
WET分量需要用到缨帽变换系数。很多人在这一步翻车，因为Landsat 5和Landsat 8的系数表完全不同，不能混用。我列一下常用的：
Landsat 5 TM的WET系数（基于TM波段1-5、7）： WET = 0.0315×B1 + 0.2021×B2 + 0.3102×B3 + 0.1594×B4 - 0.6806×B5 - 0.6109×B7
Landsat 8 OLI的WET系数（基于OLI波段2-7）： WET = 0.1511×B2 + 0.1973×B3 + 0.3283×B4 + 0.3407×B5 - 0.7117×B6 - 0.4559×B7
用Band Math输入时，直接把对应波段乘上系数再求和。注意Landsat 8的B2是蓝波段、B3是绿波段、B4是红波段、B5是近红外、B6是短波红外1、B7是短波红外2，和TM时代的编号习惯完全不同。
WET分量的物理含义是土壤和植被的湿度状况。城市建成区、裸土、道路的WET偏低，水体、植被茂密区的WET偏高。这也是RSEI里体现“湿度生态效应”的指标。
3.4 NDBSI：裸土和建筑指数的叠加
干度指标NDBSI是由裸土指数SI和建筑指数IBI合成的，公式是：
NDBSI = (SI + IBI) / 2
SI（裸土指数）的公式： SI = ((B5 + B1) - (B4 + B3)) / ((B5 + B1) + (B4 + B3))
IBI（建筑指数）的公式： IBI = (2×B5/(B5+B4) - (B4/(B4+B3) + B2/(B2+B5))) / (2×B5/(B5+B4) + (B4/(B4+B3) + B2/(B2+B5)))
这里同样注意Landsat 5和Landsat 8的波段差异。我建议分步算，先把SI和IBI分别算出来，再做平均，这样中间可以检查每步结果的数值范围是否合理。SI和IBI的值域不会严格落在[-1,1]内，会出现一些超界像元，这是正常的，归一化之后会被拉回0-1区间。
在郑州这种平原城市里，IBI对高楼密集区的响应很灵敏，但黄河水面、大型水体的IBI会出现明显负值。如果研究区内有较大面积水体，建议先用NDVI阈值做掩膜剔除，避免水体把干度指标的分布拉得很偏。
3.5 地表温度LST：大气校正法反演
地表温度反演是RSEI四指标里最绕的一步，但绕也得做，因为热度指标是城市生态评价的核心。
我用的是大气校正法，也叫辐射传输方程法，分三步。
第一步，从热红外波段辐射定标数据计算亮度温度T。ENVI里可以用Band Math手动代入Planck反函数，也可以直接用一些现成扩展工具。Landsat 5 TM的热红外是Band 6，Landsat 8是Band 10（TIRS），不要搞混。
第二步，计算地表比辐射率ε。徐涵秋在RSEI论文里推荐了简化方法：先用NDVI计算植被覆盖度FVC：
FVC = (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)
其中NDVI_soil一般取0.05，NDVI_veg取0.70。然后比辐射率： ε = 0.004×FVC + 0.986
这个简化公式对混合像元有不错的精度，在ENVI里用Band Math一条表达式就能算出来。
第三步，用大气校正法反演LST，公式是：
LST = T / (1 + (λ×T/ρ)×ln(ε))
其中T是亮度温度（K），λ是热红外波段中心波长（TM是11.45μm，OLI Band 10是10.9μm），ρ=h×c/σ，约等于1.438×10⁻² m·K。
大气校正法需要大气剖面参数，如果手头没有同步气象数据，可以用NASA提供的在线工具查询过境时刻的大气参数。对于RSEI这种相对比较分析，精度足够，结果以摄氏度或开尔文都可以，因为后面归一化会抹掉单位差异。
必须提醒一个细节：Landsat 8热红外波段原始分辨率是100米，和多光谱的30米不一致。在做LST之前，要把热红外波段用三次卷积法重采样到30米，否则后面PCA时影像尺寸对不上，会很麻烦。
3.6 四指标归一化的顺序问题
四个指标计算完之后，先做归一化再做PCA。归一化公式是：
NI = (I - I_min) / (I_max - I_min)
在ENVI里可以直接用Statistics查看每幅影像的Max和Min，然后带入Band Math。这里有两种思路：一是每期影像各自归一化，二是多期影像统一归一化。
我用的是前者，也就是各期影像各自归一化。为什么？因为长时序分析里，如果2000年影像整体偏暗、NDVI偏低，用全局统计会把低值像元的差异拉平，不利于反映真实的空间梯度。每期影像各自归一化虽然会导致不同年份之间RSEI的绝对数值不完全可比，但它的等级划分更符合该年份自身的生态特点。实际上多数RSEI论文也是这么做的。
如果你要用多期影像统一归一化，那就要在ENVI里先做影像层叠分析，把多期数据放在一起统计算，操作更复杂，但跨年份的可比性会更好。两种方案没有绝对的对错，关键是在写论文时交代清楚用了哪一种。
4. 主成分分析合成RSEI：PC1方向、贡献率与结果翻转
4.1 PCA的操作步骤
四个归一化指标准备好之后，在ENVI工具箱里用Principal Components工具，也可以直接用Layer Stacking把NDVI、WET、NDBSI、LST按顺序合成一个多波段文件，然后做主成分分析。
这里有一个关键选择：用协方差矩阵还是相关系数矩阵。RSEI的标准做法是用协方差矩阵，因为四个指标都已经归一化到0-1，量纲一致了，协方差矩阵能保留各指标本身的方差信息，而不是把方差标准化。在ENVI里选Covariance Matrix即可。
输出的PC1会是一个单波段影像，后台会生成一份特征值和特征向量的报告，这个报告很重要，一定要保存下来，里面包含了PC1的贡献率和四个指标在PC1上的载荷。
我跑郑州数据的经验是，PC1的贡献率普遍在70%-85%之间。你不需要过分纠结贡献率是不是一定要到90%，RSEI就是这样设计的，如果PC1贡献率连60%都不到，那就要回头检查四指标的计算和归一化是不是有问题。
4.2 载荷方向判断：什么时候RSEI0=1-PC1
这是RSEI流程里最容易被忽略的一步，直接影响最终结果的正负向。
PCA本身只是一种线性变换，它不关心“生态好”是正方向还是负方向。因此跑出来的PC1可能呈现出两种方向：一种是指标中绿度和湿度高的像元PC1值较大；另一种是反过来，干度和热度高的像元PC1值较大。
怎么判断？看PCA输出的特征向量。如果PC1对应的特征向量里，NDVI和WET的载荷为正，NDBSI和LST的载荷为负，说明PC1越大生态越好，可以直接用RSEI0 = PC1。但如果NDVI和WET是负的，NDBSI和LST是正的，说明PC1越大生态越差，这时候就要用RSEI0 = 1 - PC1。
郑州近20年的数据里，我碰到的情况大都是后者，即PC1中绿度和湿度是负贡献，干度、热度是正贡献，所以RSEI0 = 1 - PC1。但你做每一期影像都得单独检查方向，不能想当然用固定公式。
翻转之后再做一次0-1归一化，得到最终的RSEI值。RSEI越接近1，生态状况越好；越接近0，生态状况越差。
4.3 分级阈值怎么定，生态等级怎么解释
得到0-1连续值之后，为了分析和制图，一般会分成等级。最常用的是等间距分级，也就是：
| 等级 | RSEI范围 | 生态含义 |
|---|---|---|
| 差 | 0.0 - 0.2 | 生态质量差，以高密度建成区、裸地为主 |
| 较差 | 0.2 - 0.4 | 生态质量较差，以中低密度建成区为主 |
| 中等 | 0.4 - 0.6 | 生态质量中等，城郊过渡带、农业区为主 |
| 良 | 0.6 - 0.8 | 生态质量良好，以农田、林草覆盖区为主 |
| 优 | 0.8 - 1.0 | 生态质量优，以密林、湿地、水体为主 |
ENVI里可以用Density Slice做分级，设置好每一级的颜色范围后直接保存为色带文件。这个分级结果既可以用来出图，也可以用来统计各等级的面积。
5. 五期结果的时间序列对比：变化检测与生态解读
5.1 差值法与等级转移矩阵
做完每一期的RSEI分级后，接下来是长时序研究里最有价值的部分：变化分析。我一般做两种分析。
第一种是差值法：把2019年的RSEI减去2000年的RSEI，得到RSEI变化检测图。正差值表示生态变好，负差值表示生态变差。在ENVI里用Band Math两幅影像相减就行，然后设定变化阈值（比如±0.05以内算稳定，超过算变好或变差）。
第二种是等级转移矩阵：在ENVI里做分类后比较，统计2000年各等级到2019年各等级的转移百分比。这张矩阵能非常直观地反映生态退化的路径，比如“良”变成了“中等”、“中等”变成了“较差”，这对生态修复决策很有参考价值。
如果你用的是每期影像各自归一化再分级的结果，做等级转移矩阵没问题；但如果做RSEI数值相减的差值分析，建议所有年份都用同一套最小最大值归一化，否则数值差的语义会打折扣。实际操作中，两种分析可以互补使用。
5.2 郑州生态格局的变化细节
我用这个流程跑完2000年、2005年、2010年、2015年和2019年五期数据之后，空间分布上能看到几个明显特征。
中心城区的RSEI等级始终偏低，但低值区范围一直在外扩。尤其是郑东新区、高铁站周边、航空港区这些新建区域，生态等级从“中等”逐步跌成“较差”或“差”，这是城市扩张的直接代价。有意思的是，老城区核心位置因为绿化提升和老旧小区改造，部分像元从“差”变成了“较差”，说明局部修复是有效果的。
西部沿嵩山余脉的山地丘陵区域，因为植被覆盖好、人口密度低，RSEI稳定在“良”和“优”的水平，这20年里波动很小。北部黄河沿岸的滩区则表现出明显的年际波动，与黄河水量和滩区利用方式有关，黄河湿地保护政策落地后，2015到2019年间不少像元等级明显好转。
整体来看，郑州2000-2019年的RSEI呈现“中心恶化、边缘改善、总体略有下降但局部修复”的格局。这种空间细节丰富的结论，正是长时序RSEI相比单期生态评价最大的价值。
5.3 制图与结果输出建议
做生态评价，图比表格重要。我的制图习惯是：RSEI分级图用从红到绿的渐变色（红=差，黄=中等，绿=优），加上郑州市界和主要水系，输出300dpi以上的GeoTIFF，最后在制图软件里加比例尺、指北针和图例。
ENVI里可以直接在Layer Manager里对RSEI分级图层做配色，然后导出。如果论文里需要各等级面积表，可以在ENVI里用Class Statistics统计每级像元数，乘以900平方米（30m×30m）换算成平方公里。
一个实用技巧：把五期RSEI影像用同一个色带和同一套分级阈值制图，排版在一个图版里，读者一眼就能看出空间格局的变化轨迹。这也是RSEI论文里的标准展示方式。
6. 长时序RSEI最容易翻车的几个坑
6.1 传感器系数不一致：Landsat 5和Landsat 8不能混用
这是长时序RSEI里最隐蔽、最坑人的错误。TM和OLI的光谱响应范围不一样，WET系数表完全不是一套，你要是图省事复制粘贴同一套系数，算出来的WET在2000年和2019年之间存在系统性差异，导致主成分分析里出现虚假的变化趋势。
我自己的习惯是，在做预处理之前先建一个表，把每个年份的传感器类型、波段编号、WET系数、SI/IBI公式全部列清楚。这个表看起来费时间，但能避免你在算到第三个年份时突然发现自己用错了波段。
同样，NDBSI的波段引用也要格外小心。Landsat 5的B1对应Landsat 8的B2，差一个波段，干度指标就废了。
6.2 Landsat 7条带修复和Landsat 5老数据的大气校正问题
研究区间跨了2000-2019年，难免会用到一些中间年份的Landsat 7数据。Landsat 7在2003年SLC故障后，每景影像都有条带缺失。在ENVI里可以用Landsat Gapfill工具修复，但插值出来的数据在条带区域精度有限。如果同一时期能找到Landsat 5的影像，优先用Landsat 5。
Landsat 5的老数据问题主要出在FLAASH大气校正后的质量检查。我碰到过2010年左右一景TM影像在校正后出现明显的横向噪声带，和原始影像逐波段查看后，发现是源数据某一段扫描行出了问题。这种情况下不要犹豫，直接换相邻日期的影像重做。
6.3 水体和极端像元的处理：掩膜还是保留
水体在四个指标里的表现很特殊：NDVI是负值、WET是高值、NDBSI是负值、LST是低值。它在PCA里会成为一个极端的离群群体，把PC1的空间分布往水体方向拉。
很多RSEI论文说不做水体掩膜，理由是RSEI本身对水体也有响应，或者说PC1已经包含了水体的生态含义。但以郑州数据来看，黄河和市区人工湖面积不小，如果不做掩膜，NDBSI的最小值会被水体拉得很低，归一化后建成区和裸地的干度区分度会下降。
建议的做法是：计算NDWI=（Green - NIR）/（Green + NIR），NDWI大于0的像元视为水体，在进入PCA之前用ENVI的Build Mask功能抠掉。如果你选择保留水体，至少要在论文方法部分说明理由，并且检查PC1载荷方向有没有因水体发生反转。
6.4 物候一致性检查与异常年份复核
另一个容易被忽略的问题是：不同年份影像的获取日期不完全一致。2000年如果是8月15号，2005年可能是7月3号，虽然都在生长季，但作物物候有差异，可能造成RSEI年际变化中出现假信号。
我的处理方法是：优先选同一个旬或同一个月的影像；实在选不到，就在结果解释时明确说明物候差异的影响范围。更严格的做法是结合MODIS NDVI时间序列，对Landsat影像做物候订正，但对于大多数城市尺度的RSEI研究，控制数据源在相邻月份的精度已经足够。
最后，结果出来之后一定要复核。把RSEI分级图叠加到同区域的Google Earth历史影像或高分影像上，随机抽几十个点检查分级是否合理。尤其是那些生态等级突变明显的区域，比如从“优”直接跳到“差”的像元，很可能对应着云影或者条带修复的伪影。
我自己做一个长时序项目，每期影像跑完四指标和PCA之后，会顺手把PC1载荷、贡献率、各指标最大最小值都记录在一个Excel表格里。这个表格既是论文方法部分的素材，也是排查问题的关键线索。比如某一年PC1载荷方向突然和前后年份不一致了，多半就是归一化或者波段用错了，回去查表比重新翻ENVI历史文件高效得多。这算是我跑遥感数据这几年养成的一个笨习惯，但每次帮上忙的时候都觉得值。
