拓冰网站定制 
把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 TM2011年退役、Landsat 7 ETM2003年后有条带故障、Landsat 8 OLI/TIRS2013年发射。你自己跑的时候年份和传感器的对应关系大概是年份可选传感器注意点2000Landsat 5 TM数据质量稳定大气校正参数好定2005Landsat 5 TM同上2010Landsat 5 TM或Landsat 7如果用Landsat 7要处理SLC条带2015Landsat 8 OLI/TIRS波段编号与TM不同注意转换2019Landsat 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×B7Landsat 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) / 2SI裸土指数的公式 SI ((B5 B1) - (B4 B3)) / ((B5 B1) (B4 B3))IBI建筑指数的公式 IBI (2×B5/(B5B4) - (B4/(B4B3) B2/(B2B5))) / (2×B5/(B5B4) (B4/(B4B3) B2/(B2B5)))注意这里波段编号同样是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 6Landsat 8是Band 10TIRS不要搞混。第二步计算地表比辐射率ε。徐涵秋在RSEI论文里推荐了简化方法先用NDVI计算植被覆盖度FVCFVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中NDVI_soil一般取0.05NDVI_veg取0.70。然后比辐射率 ε 0.004×FVC 0.986这个简化公式对混合像元有不错的精度在ENVI里用Band Math一条表达式就能算出来。第三步用大气校正法反演LST。公式是LST T / (1 (λ×T/ρ)×ln(ε))其中T是亮度温度Kλ是热红外波段中心波长TM是11.45μmOLI 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. 主成分分析合成RSEIPC1方向、贡献率与结果翻转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 载荷方向判断什么时候RSEI01-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×30m900平方米换算成平方公里。另一个实用技巧把五期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的最小值会被水体拉得很低归一化后建成区和裸地的干度区分度会受到压制。建议的处理方式是计算一个简单的NDWIGreen - NIR/Green NIRNDWI大于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 TM2011年退役、Landsat 7 ETM2003年后有条带故障、Landsat 8 OLI/TIRS2013年发射。年份和传感器的对应需要注意年份可选传感器注意点2000Landsat 5 TM数据质量稳定大气校正参数好定2005Landsat 5 TM同上2010Landsat 5 TM或Landsat 7如果用Landsat 7要处理SLC条带2015Landsat 8 OLI/TIRS波段编号与TM不同注意转换2019Landsat 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×B7Landsat 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) / 2SI裸土指数的公式 SI ((B5 B1) - (B4 B3)) / ((B5 B1) (B4 B3))IBI建筑指数的公式 IBI (2×B5/(B5B4) - (B4/(B4B3) B2/(B2B5))) / (2×B5/(B5B4) (B4/(B4B3) B2/(B2B5)))这里同样注意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 6Landsat 8是Band 10TIRS不要搞混。第二步计算地表比辐射率ε。徐涵秋在RSEI论文里推荐了简化方法先用NDVI计算植被覆盖度FVCFVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中NDVI_soil一般取0.05NDVI_veg取0.70。然后比辐射率 ε 0.004×FVC 0.986这个简化公式对混合像元有不错的精度在ENVI里用Band Math一条表达式就能算出来。第三步用大气校正法反演LST公式是LST T / (1 (λ×T/ρ)×ln(ε))其中T是亮度温度Kλ是热红外波段中心波长TM是11.45μmOLI 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. 主成分分析合成RSEIPC1方向、贡献率与结果翻转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 载荷方向判断什么时候RSEI01-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的最小值会被水体拉得很低归一化后建成区和裸地的干度区分度会下降。建议的做法是计算NDWIGreen - NIR/Green NIRNDWI大于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历史文件高效得多。这算是我跑遥感数据这几年养成的一个笨习惯但每次帮上忙的时候都觉得值。 RSEI遥感生态指数实战：用ENVI和Landsat分析郑州20年生态变化
把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 TM2011年退役、Landsat 7 ETM2003年后有条带故障、Landsat 8 OLI/TIRS2013年发射。你自己跑的时候年份和传感器的对应关系大概是年份可选传感器注意点2000Landsat 5 TM数据质量稳定大气校正参数好定2005Landsat 5 TM同上2010Landsat 5 TM或Landsat 7如果用Landsat 7要处理SLC条带2015Landsat 8 OLI/TIRS波段编号与TM不同注意转换2019Landsat 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×B7Landsat 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) / 2SI裸土指数的公式 SI ((B5 B1) - (B4 B3)) / ((B5 B1) (B4 B3))IBI建筑指数的公式 IBI (2×B5/(B5B4) - (B4/(B4B3) B2/(B2B5))) / (2×B5/(B5B4) (B4/(B4B3) B2/(B2B5)))注意这里波段编号同样是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 6Landsat 8是Band 10TIRS不要搞混。第二步计算地表比辐射率ε。徐涵秋在RSEI论文里推荐了简化方法先用NDVI计算植被覆盖度FVCFVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中NDVI_soil一般取0.05NDVI_veg取0.70。然后比辐射率 ε 0.004×FVC 0.986这个简化公式对混合像元有不错的精度在ENVI里用Band Math一条表达式就能算出来。第三步用大气校正法反演LST。公式是LST T / (1 (λ×T/ρ)×ln(ε))其中T是亮度温度Kλ是热红外波段中心波长TM是11.45μmOLI 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. 主成分分析合成RSEIPC1方向、贡献率与结果翻转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 载荷方向判断什么时候RSEI01-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×30m900平方米换算成平方公里。另一个实用技巧把五期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的最小值会被水体拉得很低归一化后建成区和裸地的干度区分度会受到压制。建议的处理方式是计算一个简单的NDWIGreen - NIR/Green NIRNDWI大于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 TM2011年退役、Landsat 7 ETM2003年后有条带故障、Landsat 8 OLI/TIRS2013年发射。年份和传感器的对应需要注意年份可选传感器注意点2000Landsat 5 TM数据质量稳定大气校正参数好定2005Landsat 5 TM同上2010Landsat 5 TM或Landsat 7如果用Landsat 7要处理SLC条带2015Landsat 8 OLI/TIRS波段编号与TM不同注意转换2019Landsat 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×B7Landsat 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) / 2SI裸土指数的公式 SI ((B5 B1) - (B4 B3)) / ((B5 B1) (B4 B3))IBI建筑指数的公式 IBI (2×B5/(B5B4) - (B4/(B4B3) B2/(B2B5))) / (2×B5/(B5B4) (B4/(B4B3) B2/(B2B5)))这里同样注意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 6Landsat 8是Band 10TIRS不要搞混。第二步计算地表比辐射率ε。徐涵秋在RSEI论文里推荐了简化方法先用NDVI计算植被覆盖度FVCFVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中NDVI_soil一般取0.05NDVI_veg取0.70。然后比辐射率 ε 0.004×FVC 0.986这个简化公式对混合像元有不错的精度在ENVI里用Band Math一条表达式就能算出来。第三步用大气校正法反演LST公式是LST T / (1 (λ×T/ρ)×ln(ε))其中T是亮度温度Kλ是热红外波段中心波长TM是11.45μmOLI 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. 主成分分析合成RSEIPC1方向、贡献率与结果翻转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 载荷方向判断什么时候RSEI01-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的最小值会被水体拉得很低归一化后建成区和裸地的干度区分度会下降。建议的做法是计算NDWIGreen - NIR/Green NIRNDWI大于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历史文件高效得多。这算是我跑遥感数据这几年养成的一个笨习惯但每次帮上忙的时候都觉得值。