做遥感生态监测的人应该都绕不过一个指数RSEI遥感生态指数。这名字听起来高大上但说白了就是用一个数字告诉别人“这片区域的生态环境到底好不好”。这两年我在ENVI里用Landsat影像给郑州市做了2000年到2019年共20年的RSEI变化分析从数据下载、大气校正到主成分分析、结果出图整条流水线跑通之后才发现这个指数在长时序城市生态监测里确实好用但坑也不少。今天把这套流程从0到1完整拆开包括每一步的参数设置、工具选择、踩过的雷以及为什么这么做的底层逻辑希望能给正在做RSEI或者准备做长时序遥感生态评价的朋友省点时间。1. 项目整体设计与思路拆解1.1 为什么要用RSEI而不是单看NDVI或植被覆盖度做城市生态评价最常见的一个误区是只看植被指数。绿度确实重要但城市生态系统里还有水体、裸土、水泥建筑、热岛效应单看NDVI等于“盲人摸象”。RSEI的聪明之处在于它把四个维度压到一起绿度NDVI、湿度Wet、热度LST、干度NDBSI然后通过主成分分析PCA自动分配权重。这样一来生态好坏不再靠人拍脑袋定权重而是让数据自己说话。我当时接手这个项目时业主方就问了一个很实际的问题“你们能不能用一个指标告诉我郑州过去20年生态是在变好还是变差”RSEI恰好就是干这个的。2000年到2019年正好跨越Landsat 5 TM和Landsat 8 OLI两个传感器时代数据连续性没问题30米分辨率对市域尺度也够用。1.2 技术路线和时相选择整个技术路线分五步数据获取、预处理、四指标计算、PCA合成、结果分析。听起来简单但每一步都有讲究。先说时相选择。我最终选了2000年、2005年、2010年、2015年、2019年五期影像基本是每5年取一个节点。为什么不是每年都做一是有些年份云量太大找不到好影像二是做长时序趋势分析5年间隔足够看出变化规律而且工作量可控。如果真要逐年做建议用GEEGoogle Earth Engine批量处理ENVI里逐年手动跑太折磨人。影像时相尽量选同季节。我全部选的7-9月份的夏季影像因为夏季植被生长旺盛NDVI区分度最高而且郑州夏天热岛效应明显LST的对比也更强烈。这个细节直接影响RSEI结果的可比性。1.3 数据源选型Landsat为什么是首选市域尺度生态监测可用的数据源其实不少MODIS、Sentinel-2都能做。但RSEI这个指数最早就是基于Landsat设计的四个指标的计算公式在TM/OLI传感器上已经验证得很成熟。另外2000年的时候Sentinel-2还没上天MODIS虽然有数据但分辨率太粗250米起步对郑州市区这种精细化城市生态评价来说不够看。Landsat从1999年Landsat 7开始就保证了对地连续观测Landsat 5更是从1984年硬扛到2011年长时序优势无法替代。郑州的位置大约在Landsat轨道行列号P124R36附近这个信息在USGS和地理空间数据云里都能查到。我下载数据时习惯先在USGS EarthExplorer里框选研究区范围再限定时间和云量比直接输行列号更稳妥。2. 数据获取与预处理实操2.1 影像下载的筛选标准下载看起来是最没技术含量的环节实际上是决定成败的一步。我筛选影像有几个硬性标准云量低于10%且研究区范围内无云。千万不要只看影像整体云量有时候云恰好压在城区上RSEI算出来湿度信号全乱了。传感器优先选Landsat 5 TM2000、2005、2010Landsat 8 OLI/TIRS2015、2019。2012年前后Landsat 5已经老化严重2013年Landsat 8接棒之后数据质量明显提升。绝对避开Landsat 7的条带问题。2003年后SLC故障导致L7影像出现大量条带虽然可以后期修复但RSEI需要用到短波红外波段条带区域的数值全是错的修复效果也一般。我的原则是能不用就不用。一个真实经验2010年的好影像特别难找。郑州那年夏季云雨比较多我在USGS上翻了半天最终选了一景2010年8月下旬云量8%的影像虽然整体质量不错但城市东南角的薄云还是让NDBSI出现了一点异常。后期只能用掩膜方式处理。所以做长时序项目建议第一件事先把所有年份的可选影像列表拉出来逐个检查别等预处理做完了才发现数据不行。2.2 辐射定标和大气校正的细节拿到原始Landsat影像一般是L1级别产品第一件事不是裁剪而是辐射定标。这一步把DN值转换为传感器入瞳处的辐射亮度值是后面所有计算的基石。ENVI里操作用Radiometric Calibration工具需要特别注意两个设置定标类型选Radiance输出数据类型选Float浮点型。如果选了Bias文件或保持整型输出后面大气校正极容易报错。大气校正我用的FLAASH模块。这一步很多人觉得可做可不做实际影响巨大。FLAASH校正之后植被波谱曲线会变得正常NDVI和Wet计算值更接近真实地表反射率。如果不做LST反演和NDBSI都会有系统性偏差长时序对比就更不可信了。FLAASH参数设置是实操里最容易出问题的地方我列一下要点影像中心经纬度在FLAASH界面里直接从影像头文件读取不要手动乱填。传感器类型Landsat 5 TM选TMLandsat 8 OLI选OLI选错会导致波段中心波长错位。飞行日期和时间从头文件里复制GMT时间会自动换算。大气模型夏季中纬度地区我一般选Mid-Latitude Summer。气溶胶模型城市区域通常选Urban能更好处理人为污染带来的气溶胶影响。初始能见度默认40公里如果当天影像里明显有霾改成20或更小否则校正结果发白。校正完成之后用Z Profile工具检查一下典型地物波谱比如植被在近红外处反射率高、红波段处吸收明显如果曲线已经是“正常”的植被形态说明校正到位。2.3 矢量裁剪与坐标系统一郑州市行政边界矢量我提前准备好投影坐标系用WGS84 UTM 50N郑州经度约113.6°E属于UTM 50N分带。这里有个容易忽略的坑下载的Landsat影像自带的地理坐标系是WGS84但很多边界矢量是CGCS2000或者Xian 1980投影不一致会导致裁剪出来错位好几公里最后还是对不上。正确做法是先把影像和矢量都统一到同一坐标系我习惯统一到UTM投影因为后面做面积统计时会方便很多面积单位直接是平方米不用再做度数换算了。裁剪工具用ENVI的Subset Data from ROIs确保勾选“Mask pixels outside of ROI”选项这样研究区之外的背景值会被自动设为0。不勾的话边缘会有一圈黑边后面做归一化时这些0值会干扰极值统计。3. 四大生态指标的计算方法与原理3.1 绿度指标NDVINDVI是四个指标里最“皮实”的一个公式简单、物理意义明确计算公式为NDVI (NIR - Red) / (NIR Red)对应的Landsat波段TM传感器用(B4-B3)/(B4B3)OLI传感器用(B5-B4)/(B5B4)。在ENVI里直接用Band Math写公式或者用Spectral Index工具一键生成。这里我有个小习惯生成NDVI之后会先看一下均值和标准差。城市区域夏季NDVI均值一般在0.3-0.5区间如果均值超过0.6说明大气校正可能没做干净反射率偏高。3.2 湿度指标Wet湿度指标来自缨帽变换K-T变换的湿度分量它能有效反映土壤和植被的含水量信息对城市水体、植被蒸腾和土壤湿润程度都很敏感。不同传感器的Wet系数差异很大一定要区分开Landsat 5 TM的Wet公式Wet 0.0315×B1 0.2021×B2 0.3102×B3 0.1594×B4 - 0.6806×B5 - 0.6109×B7Landsat 8 OLI的Wet公式Wet 0.1511×B2 0.1973×B3 0.3283×B4 0.3407×B5 - 0.7117×B6 - 0.4559×B7这两个公式我从网上和文献里核过好几遍确认无误。实际操作中最容易犯的错误是搞混波段顺序比如OLI数据的B2是蓝波段而不是TM里的B1蓝。我用Band Math输入公式之前会先打开波段列表核对一遍各波段名称宁可多花两分钟也别等算完才发现波段选错了。3.3 热度指标LSTLST地表温度的获取要比前两个复杂不少。我采用的是大气校正法辐射传输方程法原理是先从热红外波段获取辐射亮度值再估算大气上行辐射、下行辐射和大气透射率最终反演地表温度。ENVI里没有一键LST工具需要分步计算但逻辑其实很清晰第一步热红外波段辐射定标得到辐射亮度Lλ。第二步用Planck公式计算亮温T。亮温公式可以写成T K2 / ln(K1 / Lλ 1)TM热红外波段K1607.76K21260.56OLI的TIRS Band10K1774.8853K21321.0789。第三步计算植被覆盖度再估算地表比辐射率ε。这里用NDVI阈值法简单实用当NDVI小于0.2时视为裸土ε取0.97。当NDVI大于0.5时视为全植被ε取0.99。介于中间时按植被覆盖度线性插值ε 0.004×FVC 0.986。第四步用辐射传输方程反演地表温度公式为LST T / [1 (λ×T/ρ)×ln(ε)] - 273.15其中λ是热红外波段中心波长ρ取1.438×10⁻² m·Kε是比辐射率。这串公式看着头大但实际在ENVI的Band Math里一步步写下来也就每条一行的事。如果有条件也可以用ENVI扩展工具里的LST模块能省一半时间。但无论如何最后算出来的LST数值要检查一下郑州夏季地表温度一般40-55℃之间如果出现70℃这种离谱值肯定是参数填错了。3.4 干度指标NDBSI干度指标代表裸土和建筑用地的“干燥程度”是城市生态评价里特别重要的负向指标。计算公式为NDBSI (IBI SI) / 2其中SI是土壤指数公式SI [(SWIR1 Blue) - (NIR Green)] / [(SWIR1 Blue) (NIR Green)]TM传感器对应(B5B1-B4-B2)/(B5B1B4B2)OLI传感器对应(B6B2-B5-B3)/(B6B2B5B3)。IBI是基于建筑用地的指数计算公式IBI {2×SWIR1/(SWIR1NIR) - [NIR/(NIRRed) Green/(GreenSWIR1)]} / {2×SWIR1/(SWIR1NIR) [NIR/(NIRRed) Green/(GreenSWIR1)]}一句话概括NDBSI把“裸土建筑”的信号综合在一起。城市里高楼大厦、柏油路越多NDBSI值越高生态环境就越差。我第一次算完NDBSI之后用彩色渲染看了一眼郑州市中心明显是亮色高值区跟实际情况完全对得上这就是好结果。3.5 指标归一化处理四个指标算出来之后它们的量纲完全不同NDVI在-1到1之间LST是几十度的数值必须统一归一化到0-1区间。公式NI (I - I_min) / (I_max - I_min)ENVI里用Band Math做需要先对每个指标做一次统计记下最小值和最大值。注意统计时务必忽略背景0值否则背景值会被当成真实数据归一化结果会整体偏移。我一般用掩膜文件把研究区外的区域排除后再统计。另外因为LST和NDBSI是负向指标数值越高生态越差做主成分之前要把它们反转成1-LST、1-NDBSI确保四个指标方向一致都是数值越大代表生态越好。4. 主成分分析与RSEI构建4.1 主成分分析为什么能自动定权重四个指标都归一化之后接下来就是RSEI的核心——主成分分析PCA。为什么要用PCA因为绿度、湿度、热度、干度这四者之间存在信息重叠和相关性PCA能通过线性变换提取它们的主要共同特征把四个波段压缩成一个主成分主成分的方差越大承载的信息越多。最终我们取第一主成分PC1它的特征值贡献率一般能达到80%以上这就是RSEI的雏形。这个方法不需要人为设定权重PC1会自动把载荷高的指标作为主导因素理论上更客观。这也是RSEI相比其他生态指数最大的卖点。4.2 ENVI主成分分析实操ENVI里的操作路径是Toolbox - Transform - Principal Components - Forward PC Rotation。关键设置是Input File选择四个归一化后的指标图层按顺序组合成一个多波段文件Layer Stacking。协方差矩阵选择“Covariance Matrix”而非相关矩阵。输出选PC1同时勾选输出特征值和特征向量。计算完成后第一步看特征值贡献率。PC1贡献率如果低于70%说明四个指标的独立性太强主成分的代表性不足。这时候可以回头检查归一化是否到位或者考虑是不是某一期影像质量问题被混进来了。我做的五期数据里PC1贡献率基本都在75%-88%之间2015年最低原因是当时影像西北角有一片云影残留导致Wet指标异常。4.3 初始RSEI的计算和方向校验PC1计算出来后并不直接就是RSEI。常规处理是RSEI0 1 - PC1原因很简单在正确的计算方向下PC1里NDVI和Wet的载荷为正LST和NDBSI的载荷为负也就是说PC1越高生态越好。如果方向反了PC1会变成“生态差”的象征这时候RSEI0 PC1而不是1-PC1。方向校验的具体做法是在PC1结果图里用ROI工具在城市建成区和高植被覆盖区各采几个点看数值高低。如果城市建成区PC1值高、植被区PC1值低说明方向反了就需要用1-PC1。另外也可以直接查特征向量表判断四个指标在PC1上的载荷正负。最后对RSEI0再做一次归一化就得到最终的RSEI值范围在0-1之间数值越高代表生态越好。4.4 RSEI等级划分RSEI算出来之后是连续栅格直接看趋势不方便我一般会按0.2的间隔分成五级等级RSEI区间生态评价优0.8-1.0植被覆盖高生态质量好良0.6-0.8生态质量较好中等0.4-0.6生态质量一般较差0.2-0.4生态质量偏差差0.0-0.2生态质量差多为建设用地或裸地然后通过ENVI的Class Raster工具对RSEI分级再统计各等级像元数和面积占比就能精确地说出郑州从2000年到2019年生态优良区是变多了还是变少了变化区域主要在哪些方位。5. 长时序变化分析与影响解读5.1 五期RSEI结果对比我把2000年、2005年、2010年、2015年、2019年五期RSEI结果逐像元对比做了差值图后一期减前一期正负值直接反映生态变好还是变差。结论和郑州市实际的城市扩张方向高度吻合。市区核心区二七区、金水区、管城区RSEI整体偏低且2005年之后持续下滑主要原因当然是城市建成区扩张硬质地面增加NDVI和Wet下降、LST和NDBSI上升。但有意思的是郑东新区是一个明显的“先降后升”案例。2005年前后大规模建设期RSEI一度跌到谷底2015年以后随着龙子湖、象湖等水系公园的陆续建成区域RSEI又有回升。这说明RSEI不只是反映“开发破坏”也能捕捉到生态修复的积极信号。西郊和南部的山地、农田区域RSEI相对稳定偏优尤其是嵩山余脉一带2000年到2019年RSEI基本维持在0.7以上是整个郑州市的生态压舱石。5.2 如果要做GIS空间统计RSEI栅格本身不能直观说明“哪个街道生态最好”需要配合GIS做空间统计。我在ArcGIS或QGIS里把RSEI按行政区区县做分区统计算出每个区的RSEI均值再用柱状图对比。这里有个细节分区统计前建议把RSEI重采样成100米分辨率避免30米像元在边界处被大幅裁剪影响面积占比精度。转移矩阵部分我也顺手做了。把两期RSEI分级结果叠加用ENVI的Confusion Matrix或直接用ArcGIS的Tabulate Area可以得到各个等级之间的转移情况。比如2005年到2010年郑州市“优”级转“良”级的面积有多少“较差”级转“差”级的面积有多少。这些数据写报告的时候特别加分。5.3 数据可信度自查长时序分析最容易被人质疑的就是不同年份传感器不一样结果能比吗我做结果校验时主要从两个角度交叉验证一是看趋势合理性。郑州2000-2019年整体处于快速城市化和生态修复并行的阶段RSEI均值呈微弱下降但局部显著变化的趋势符合实际。如果结果出现全区域RSEI暴涨或暴跌的极端情况基本可以判定数据处理有问题。二是做相对验证。我拿了2015年的MODIS NDVI产品跟Landsat NDVI做了相关性分析相关系数约0.82说明Landsat估算结果在趋势上与独立数据源一致。这不能证明RSEI绝对准确但至少说明长时序变化规律是可靠的。6. 数据准备、参数选择与避坑经验6.1 不同传感器波段差异对照长时序项目最大的坑之一就是Landsat 5和Landsat 8的波段设置不同很多公式不能直接套用。我从项目一开始就整理了一张对照表波段功能Landsat 5 TMLandsat 8 OLI蓝波段B10.45-0.52B20.45-0.51绿波段B20.52-0.60B30.53-0.59红波段B30.63-0.69B40.64-0.67近红外B40.76-0.90B50.85-0.88短波红外1B51.55-1.75B61.57-1.65短波红外2B72.08-2.35B72.11-2.29热红外B6B10/B11尤其是热红外波段TM只有一个B6OLI有两个热红外波段B10和B11做LST时通常只用B10因为B11的定标参数和大气影响更不稳定。这个细节很多教程不会说但我实际测试过用B11反演出来的地表温度整体偏高2-3℃明显不靠谱。6.2 坐标系和掩膜的坑长时序栅格分析最怕坐标系不统一。2000年的影像下载下来可能是WGS842019年的可能是UTM投影USGS部分产品自带的投影方式不一致叠加对比时直接错开几百米。所以我每期数据做完大气校正之后第一件事就是统一投影和像元大小。用ENVI的Resize Data把像元统一为30米同时在UTM投影下做图像配准。不夸张地说这个步骤能省掉后面80%的对齐问题。6.3 影像质量的三个硬指标我总结了三个判断影像能不能用的硬指标拿不准的时候直接套云量小于10%并且研究区内无云影和阴影。影像获取日期在目标月份±15天以内尽量避开季相差异过大的时段。从Landsat头文件里看太阳高度角夏季影像太阳高度角低于50°的慎用光照角度太低会放大地形阴影和城市建筑阴影。这三个指标同时满足的影像才是“放心数据”。我曾经为了省事用过一景云量6%但太阳高度角只有42°的影像结果城市密集区阴影特别重NDVI明显偏低最终RSEI结果跟前后年份断档只能重新下载数据再做一遍。7. 常见问题与排查技巧实录7.1 ENVI报错排查方法很多人在ENVI里做RSEI会卡在奇奇怪怪的报错上。我把实际遇到频率最高的几类问题和排查思路整理成了一张速查表方便朋友们直接对照问题现象可能原因排查和解决方案FLAASH大气校正在输入数据处报错输入文件没有做辐射定标或数据类型不是Float重新做Radiometric Calibration输出Float型数据大气校正结果整体偏白或偏暗能见度参数不合适调整初始能见度多试几次参考影像实际天气PCA计算报错“Input stack has different... ”多波段文件里像元大小或范围不一致先Layer Stacking时统一投影、范围和像元大小归一化后影像全黑或全白统计极值时包含了背景0值裁剪边界生成掩膜只在掩膜范围内统计极值波段计算输出的结果分布诡异公式中波段号与传感器不匹配核对TM/OLI传感器的波段编号确认选的是B1-B7对应关系LST反演结果整体偏高热红外波段选错或大气参数不匹配改用OLI/TIRS B10检查比辐射率取值范围这些报错信息里最坑的不是报错本身而是ENVI有时候压根不报错输出一张看起来正常但数值不对的图。比如波段选错后Band Math照样能跑结果却是错的。所以每次在做完一个指标后我习惯随手打开直方图看数值范围是否符合预期然后再进入下一步。7.2 RSEI方向校验收到的教训有一次做2010年那期数据PC1特征向量里NDVI和Wet的载荷是负值LST和NDBSI是正值这意味着PC1数值越高生态越差方向完全反了。第一次算RSEI时我没注意直接1-PC1归一化结果发现RSEI空间分布和实际情况正好相反市中心显示生态良好山区反而生态差。后来养成了固定习惯每次算完PC1先看特征向量表确认方向再用散点图采集城市中心和植被区的值做对比。这个动作30秒就够但能避免整个结果推倒重来。7.3 关于工作效率的三个建议做这种长时序项目数据量不算大但步骤繁琐效率低纯属磨洋工。我建议第一用ENVI的Batch/Modeler批处理功能。五期数据的辐射定标、大气校正、波段计算完全可以用同一个流程模板批量跑设置好之后去吃个饭回来就全处理完了。第二保存中间结果。每做完一步就把结果另存为独立文件不要一路覆盖。后面想回头排查是哪一步出了问题有中间文件才能定位。第三记录每期数据的参数日志。哪期影像用了哪个公式、LST的大气透射率填的是多少、归一化的最大最小值是多少全部记在一个Excel里。长时序项目周期长过两个月你自己都会忘记当时怎么处理的数据有个日志比什么都强。8. 从RSEI到决策这套结果能拿来做什么RSEI做出来不是放在PPT里好看用的它在城市生态规划里有很多实际应用场景。郑州市的案例里我重点画出了2005-2010年间三环以内RSEI从“中等”跌到“较差”的区域这些区域高度集中在老工业区和城中村改造片区。绿化规划部门拿这个结果去排查这些区域的具体用地性质反向指导绿地布局比单纯看土地利用分类图更直接因为RSEI直接反映的是生态综合质量而不只是土地覆盖类型。另外RSEI还能用来做生态质量监测的时间预警。郑州近郊几个大型居住区建成后RSEI从0.6降到0.35左右基本用两年时间就能稳定下来。对这个下降幅度和时间规律我建议后续可以和气象站点数据结合分析热岛强度与RSEI的定量关系这又是另一个有意思的方向了。如果你准备在自己的城市复刻这套流程我强烈建议从Landsat 8 OLI的两期影像开始练手先跑通流程再回头补Landsat 5的历史数据。毕竟现在绝大多数在线教程都基于OLI传感器报错信息也更容易搜到答案新手直接上手TM老数据容易在第一步就被FLAASH参数劝退。最后再分享一个习惯RSEI模型的最优解未必是PC1如果某期影像PC1贡献率特别低可以尝试PC1与PC2加权合成或者检查四个指标是否需要做主成分旋转。但这属于进阶玩法了常规研究以PC1为主就完全可以支撑结论。数据质量过关、操作流程严谨RSEI这套方法绝对能成为你手里分析城市生态最顺手的工具之一。