前年接了一个县域生态质量评估的项目甲方要的指标很简单这片区域的生态是在变好还是在变差。翻了半天规范最后落在徐涵秋老师提出的遥感生态指数RSEI上。这个指数不依赖统计年鉴数据完全用遥感影像反演的四个分量来刻画生态状况——绿度、湿度、热度、干度——再通过主成分分析把它们集成为一个0到1的值高就代表生态好低就代表生态差。整套流程跑下来最繁琐、也最容易出错的恰恰就是开头这“四个指数的计算”。NDVI、WET、NDBSI、LST每一个单独拎出来都有对应的成熟公式但把它们拼到同一套流程里波段顺序、传感器型号、大气校正方式、归一化区间这些细节一叠加问题就来了。这篇博文就围绕这四个指数的计算展开把每一步的原理、公式、实操写法、常见坑位一次讲透适合正在用ENVI或者GEE做区域生态评价的同学参考。1. RSEI在算什么四个指数背后的一盘棋1.1 为什么偏偏是绿度、湿度、热度、干度这四个分量RSEI的全称是遥感生态指数英文Remote Sensing Ecological Index。这个指数的核心思路很朴素人在自然界中感知一个地方生态好不好通常就看几件事——植被茂不茂盛、空气湿润还是干燥、体感热不热、地面是土石裸露还是建筑密集。徐涵秋把这四个维度对应到了遥感可量化的四个指数上绿度用NDVI归一化植被指数代表植被覆盖越密、长势越旺值越高。湿度用WET缨帽变换湿度分量代表反映土壤和植被冠层的水分含量。热度用LST地表温度代表地表越热生态舒适度越差。干度用NDBSI归一化差值裸土与建筑指数代表裸土和建筑用地越多地表越“干化”。这四个分量之所以被组合到一起是因为它们能从不同侧面捕捉人类活动对生态的干扰。一片区域如果被开发建设直观表现就是植被减少、湿度下降、地表增温、裸土和硬化地面增加恰好对应NDVI和WET下降、LST和NDBSI上升。所以RSEI本质上是在用遥感波段响应来量化“人类活动干扰程度”。生态维度遥感指数生态含义与生态质量的关系绿度NDVI植被覆盖度与活力正向湿度WET地表水分含量正向热度LST地表热环境负向干度NDBSI裸土与建筑占比负向1.2 从影像到RSEI的完整技术链路在动手计算四个指数之前一定先在心里把整条流程过一遍否则很容易在中间某一步迷路。完整的RSEI计算流程是影像预处理辐射定标、大气校正、云掩膜、研究区裁剪。计算四个指数NDVI、WET、NDBSI、LST。四个指数各自归一化到0到1区间。对归一化后的四个指数做主成分分析PCA。取第一主成分PC1根据载荷方向得到初始RSEI。对初始RSEI再做一次归一化得到最终0到1的生态指数。这里要特别提醒一点很多人以为RSEI就是“四个指数算完加个权”实际上RSEI的权重不是人为设定的而是由主成分分析自动确定的。这个设计的好处是避免主观赋权带来的争议也是RSEI跟其他生态评价指数最大的区别。后续的归一化和PCA都会围绕这个逻辑展开。2. 数据准备影像选型和预处理决定成败2.1 用哪个传感器Landsat系列的波段差异RSEI前面这“四个指数的计算”全部基于Landsat系列影像因为30米分辨率和波段设置对这几个指数最友好。现在常用的有Landsat 5 TM、Landsat 7 ETM、Landsat 8 OLI、Landsat 9 OLI-2。波段数量虽然不同但计算四个指数需要的核心波段是齐全的蓝波段Landsat 5/7是Band 1Landsat 8/9是Band 2。绿波段Landsat 5/7是Band 2Landsat 8/9是Band 3。红波段Landsat 5/7是Band 3Landsat 8/9是Band 4。近红外Landsat 5/7是Band 4Landsat 8/9是Band 5。短波红外1Landsat 5/7是Band 5Landsat 8/9是Band 6。短波红外2Landsat 5/7是Band 7Landsat 8/9是Band 7。热红外Landsat 5 TM是Band 6Landsat 7是Band 6Landsat 8/9是Band 10或Band 11。这是最容易栽跟头的地方。Landsat 8把原来的近红外和短波红外波段整体往前提了一位如果沿用Landsat 5/7的波段号去套公式算出来的指数全是错的。我见过不少初学者把Landsat 8的Band 4当近红外用结果NDVI算出来一片异常——本质就是把红波段和近红外波段整反了。2.2 预处理哪几步不能省四个指数虽然只是波段运算但输入影像是辐射定标和大气校正后的地表反射率产品时结果才有意义。第一辐射定标。将DN值转为辐射亮度或反射率这是所有定量反演的基础。如果用的是USGS的Collection 2 Level 2产品地表反射率波段已经做好了大气校正可以跳过这一步如果用Level 1原始数据就必须做。第二大气校正。LST反演和NDVI计算都非常依赖大气校正。ENVI里常用FLAASH模块输入影像、传感器类型、成像时间、中心经纬度、大气模型、气溶胶模型这几个参数就能跑出来。GEE用户更简单直接用Collection 2 Level 2的SR波段即可。第三云掩膜。Landsat影像上的云和云影会导致NDVI和LST出现异常值尤其LST在云覆盖区会显示为低温在RSEI结果里就会表现为一片“虚假的生态良好区”这是非常隐蔽的错误。GEE里可以使用QA_PIXEL波段做CFMask云掩膜ENVI里则需要结合目视检查和阈值方法手动剔除。第四研究区裁剪和水体掩膜。水体在WET上表现为高值、在LST上表现为低温直接参与PCA会干扰第一主成分的载荷方向。大多数RSEI研究都会用改进的归一化差值水体指数MNDWI把水体掩掉只对陆地区域计算。这一点后面踩坑清单里还会再提。2.3 多时相影像的季节一致性如果项目要做的是多年变化分析那影像时相的选择要尽量一致最好都选在7月到9月植被生长旺盛期。因为NDVI受物候影响非常大春季和秋季的植被指数差异足以淹没真实的生态变化信号。我自己的习惯是优先选8月前后、云量小于10%的影像如果某一年实在没得选再放宽到7月或9月但会在结果里说明这个限制。3. 绿度和湿度NDVI与WET的计算方法3.1 NDVI最熟悉也最容易马虎的指数NDVI的计算公式是NDVI (NIR - Red) / (NIR Red)Landsat 8/9NDVI (Band 5 - Band 4) / (Band 5 Band 4)Landsat 5/7NDVI (Band 4 - Band 3) / (Band 4 Band 3)这个公式本身没有任何难度但在ENVI Band Math里要注意数据类型。如果直接拿整型影像做减法结果会自动截断为整型NDVI会变成一堆0和1后面根本没法用。正确写法是在表达式里强制转浮点(float(b5) - b4) / (float(b5) b4)GEE里使用归一化差值函数可以避免这个问题var ndvi img.normalizedDifference([SR_B5, SR_B4]).rename(NDVI);NDVI的取值范围是-1到1理论上裸土接近0、浓密植被接近0.8甚至更高。但在后续RSEI计算中我们并不直接用原始NDVI而是用它做归一化和植被覆盖度估算所以不需要对NDVI本身做额外处理计算精度达到float就够。3.2 WET缨帽变换湿度分量系数表必须查对WET来自缨帽变换Tasseled Cap Transformation它是通过对原始波段做固定系数线性组合得到的代表地表湿度信息。国内做RSEI普遍采用徐涵秋团队论文中的WET系数不同传感器的系数完全不同传感器BlueGreenRedNIRSWIR1SWIR2Landsat 5 TM0.03150.20210.31020.1594-0.6806-0.6109Landsat 7 ETM0.26260.21410.09260.0656-0.7629-0.5388Landsat 8 OLI0.15110.19730.32830.3407-0.7117-0.4559以Landsat 8为例WET在ENVI Band Math里的表达式是0.1511 * b2 0.1973 * b3 0.3283 * b4 0.3407 * b5 - 0.7117 * b6 - 0.4559 * b7这里b2到b7依次对应Blue、Green、Red、NIR、SWIR1、SWIR2。GEE里用expression按波段名引用同样可行var wet img.expression( 0.1511 * B2 0.1973 * B3 0.3283 * B4 0.3407 * B5 - 0.7117 * B6 - 0.4559 * B7, { B2: img.select(SR_B2), B3: img.select(SR_B3), B4: img.select(SR_B4), B5: img.select(SR_B5), B6: img.select(SR_B6), B7: img.select(SR_B7) } ).rename(WET);这里要特别强调WET系数是跟着传感器走的不同传感器的系数不能互换。把Landsat 8的系数套到Landsat 5影像上WET值的分布会完全失真后续PCA载荷方向的判断也会跟着出错。3.3 计算前的数据质量自查算出NDVI和WET后先不要急着往下走用统计工具看一下最大值、最小值、均值和标准差。NDVI应该有正值也有负值水体区域接近-1WET在城市建成区应该偏低在水体和林区偏高。如果发现WET在全图都是负值或者整体量级不对优先检查用的是不是对应传感器的系数以及波段顺序是不是对应对了。这一步自查能省下后面排查的一堆时间。4. 干度和热度NDBSI与LST的完整计算链路4.1 NDBSI裸土指数SI与建筑指数IBI的合成干度指标用NDBSI表示它由裸土指数SI和建筑指数IBI取平均得到NDBSI (SI IBI) / 2SISoil Index裸土指数的公式为SI [(SWIR1 Red) - (NIR Blue)] / [(SWIR1 Red) (NIR Blue)]IBIIndex-based Built-up Index建筑指数的公式为IBI {2 × SWIR1 / (SWIR1 NIR) - [NIR / (NIR Red) Green / (Green SWIR1)]} / {2 × SWIR1 / (SWIR1 NIR) [NIR / (NIR Red) Green / (Green SWIR1)]}以Landsat 8为例ENVI Band Math里可以分两步写先分别算出SI和IBI; SI ((b6 b4) - (b5 b2)) / ((b6 b4) (b5 b2)) ; IBI (2 * b6 / (b6 b5) - (b5 / (b5 b4) b3 / (b3 b6))) / (2 * b6 / (b6 b5) (b5 / (b5 b4) b3 / (b3 b6)))然后再做一次平均得到NDBSI。这里有个经验点IBI公式里绿波段对SWIR1、近红外对红波段的比值方向不能写反否则建筑区的响应会从高值变成低值。我建议算完之后在城市建成区取几个样本点核对数值NDBSI在建成区应该明显高于农田和林地。GEE里对应的写法是把表达式直接拼接注意波段引用名保持一致var si img.expression( ((SWIR1 Red) - (NIR Blue)) / ((SWIR1 Red) (NIR Blue)), { SWIR1: img.select(SR_B6), Red: img.select(SR_B4), NIR: img.select(SR_B5), Blue: img.select(SR_B2) } ).rename(SI); var ibi img.expression( (2 * SWIR1 / (SWIR1 NIR) - (NIR / (NIR Red) Green / (Green SWIR1))) / (2 * SWIR1 / (SWIR1 NIR) (NIR / (NIR Red) Green / (Green SWIR1))), { SWIR1: img.select(SR_B6), NIR: img.select(SR_B5), Red: img.select(SR_B4), Green: img.select(SR_B3) } ).rename(IBI); var ndbsi si.add(ibi).divide(2).rename(NDBSI);4.2 LST从辐射定标到地表温度反演的完整步骤热度指标LST的严格反演要经过辐射定标、亮度温度、植被覆盖度、地表比辐射率、大气参数代入这五步。第一步辐射定标。把热红外波段的DN值转为大气顶辐射亮度LλLλ ML × Qcal AL其中Qcal是热红外波段的DN值ML和AL是元数据里的增益和偏置参数。在USGS的Collection 2元数据里Landsat 8/9的ML一般存为RADIANCE_MULT_BAND_10或RADIANCE_MULT_BAND_11AL同理。第二步计算亮度温度BTBT K2 / ln(K1 / Lλ 1)K1和K2是热红外波段的定标常数不同传感器取值不同这个必须查对。传感器K1 (W/(m²·sr·μm))K2 (K)Landsat 5 TM607.761260.56Landsat 7 ETM666.091282.71Landsat 8 TIRS Band 10774.88531321.0789Landsat 9 TIRS-2 Band 10774.88531321.0789第三步计算植被覆盖度Pv。常用的方法是像元二分模型Pv (NDVI - NDVI_min) / (NDVI_max - NDVI_min)这里的NDVI_min和NDVI_max取累计百分比5%和95%对应的值而不是直接取全图极值。因为影像上难免有异常像元直接取极值会把植被覆盖度整体带偏。第四步估算地表比辐射率ε。用植被覆盖度做线性加权ε 0.004 × Pv 0.986这个公式的经验来源是完全裸土的自然地表比辐射率约为0.986完全植被覆盖的地表约为0.99中间按植被覆盖度线性内插。第五步代入大气参数。这一步最关键也最容易忽略。地表温度不是把BT简单换算就能得到的因为大气上行辐射L↑、大气下行辐射L↓、大气透过率τ会干扰传感器接收的信号。完整的辐射传输方程是Lλ [ε × B(Ts) (1 - ε) × L↓] × τ L↑解出B(Ts)B(Ts) [Lλ - L↑ - τ × (1 - ε) × L↓] / (τ × ε)最后反演地表温度LST K2 / ln(K1 / B(Ts) 1)这三个大气参数L↑、L↓、τ可以从NASA的大气校正参数查询网站获取按成像时间和中心经纬度查询。这个网站是公开的输入坐标和时间就能拿到对应的大气剖面参数取当天过境时刻的值即可。4.3 更省事的做法直接用Collection 2 Level 2的温度产品现在USGS的Landsat Collection 2 Level 2产品里已经直接提供了地表温度波段ST_B10不需要自己从头反演。GEE里处理起来非常简单var lst img.select(ST_B10) .multiply(0.00341802) .add(149.0) .subtract(273.15) .rename(LST);这个表达式是把ST_B10原始值换算成摄氏温度。0.00341802和149.0是L2产品的定标系数USGS官方文档里写得很明确。不过要提醒的是如果你希望跟文献里的RSEI结果保持可比性很多已发表论文用的是辐射传输方程法直接用L2温度产品得到的结果跟那些论文之间可能会有零点几度的系统性偏差。我的建议是做单期静态评价时用哪种方法都可以做多年动态对比时要么全部用L2产品要么全部自己反演别混着用。5. 归一化处理四指数量纲统一的必经一步5.1 为什么必须归一化四个指数的物理单位完全不同NDVI在-1到1之间WET可能是负几十到正几十LST是摄氏度NDBSI理论上在-1到1之间但实际分布也很分散。如果不做归一化直接丢进PCA数值范围大的变量会主导第一主成分的方向RSEI就变成了WET或LST的独角戏失去了多指标综合的意义。归一化就是把每个指数的原始值线性映射到0到1区间NI (I - I_min) / (I_max - I_min)I_min和I_max的选取非常讲究。如果直接取全图最小值和最大值云影、水体、异常像元会把区间拉得很宽导致正常地物被压缩到一个很窄的区间里。绝大多数RSEI研究会取累计概率2%到98%或者5%到95%作为区间端点把两端离群值裁掉。这一步对整个结果影响很大同一个研究区用全图极值还是用5%到95%分位最终RSEI的空间分布可能会有明显差异。5.2 ENVI和GEE里的实操写法在ENVI里先对每个指数做统计读取2%和98%分位数然后在Band Math里用裁剪归一化的方式写; 以NDVI为例 (float(b1) - 0.15) / (0.75 - 0.15)这里0.15和0.75只是示例实际值要根据统计结果替换。如果b1小于I_min公式会算出负值需要再用条件语句把结果约束在0到1之间。ENVI里有foran expression可以写条件判断但更简单的做法是接受部分像元落在0到1之外因为后续PCA本质上不关心取值范围是否严格只要所有输入的量纲一致即可。GEE里用表达式可以更方便地截断var ndviNorm ndvi.subtract(ndviMin).divide(ndviMax.subtract(ndviMin)) .clamp(0, 1) .rename(NDVI_norm);clamp函数会把小于0的置为0、大于1的置为1正好对应端点截断的效果。5.3 一个常见的分位数取值陷阱同一期的四景影像拼接成一幅研究区图时归一化的统计值应该基于整幅拼接影像计算而不是基于单景影像计算。否则每景影像各自归一化接缝处会出现明显的色调差后续PCA也会在接缝处产生伪边界。这个细节在GEE里尤其要注意统计分位数时把研究区内所有像元放在一个reducer里算不要分块统计。6. 主成分分析集成把四个指数变成一个生态值6.1 PCA在这里扮演什么角色四个指数分别刻画了生态的不同侧面但如果把它们并列展示决策者很难直接得出结论。主成分分析在这里的作用是在保持信息损失最小的前提下把四个相关变量压缩成一个综合变量。RSEI用的是第一主成分PC1它在理想情况下应该解释四个指数总方差的70%以上。如果PC1的解释率偏低说明四个指数之间的相关性较弱此时RSEI的可靠性也会打折扣。这个解释率在输出的本征值信息里可以直接读取。还有一个关键判断PC1的载荷方向。理论上NDVI和WET在PC1上的载荷应为正LST和NDBSI应为负。因为生态好的地方绿度高、湿度高、温度低、干度低。如果实际算出来的载荷方向跟这个预期相反说明取反了。常规处理是若PC1与生态质量正相关NDVI、WET载荷为正则RSEI PC1。若PC1与生态质量负相关则RSEI 1 - PC1。我在实操中习惯不管方向统一用“正向指标”对齐的思路先把四个归一化指数中预期方向一致的变量确认好再看PC1载荷如果方向反了就取1减PC1。千万不要想当然地认为所有的PC1都可以直接当RSEI用。6.2 ENVI里的操作要点ENVI里做PCA一般走两条路一是用Forward PC Rotation工具直接对四波段影像做主成分变换二是先从Compute Statistics里算出协方差矩阵再做主成分变换。推荐用第二种因为可以明确看到本征值和载荷矩阵。操作上把四个归一化指数合成为一个四波段文件执行PCA变换后第一波段就是PC1。检查PC1的载荷方向决定是否取反向然后把PC1按0到1做最终归一化得到的就是RSEI。6.3 GEE里的主成分分析思路GEE里没有现成的“一键PCA”工具需要手动实现。核心步骤是用reduceRegion计算四个归一化指数的协方差矩阵用eigen分解得到特征向量再用矩阵乘法把四波段影像投影到第一主成分轴上。代码大概长这样var image ee.Image.cat([ndviNorm, wetNorm, ndbsiNorm, lstNorm]).rename([NDVI, WET, NDBSI, LST]); var region table.geometry(); var covar image.reduceRegion({ reducer: ee.Reducer.centeredCovariance(), geometry: region, scale: 30, maxPixels: 1e10 }); var covarArray ee.Array(covar.get(covariance)); var eigens covarArray.eigen(); var eigenVectors eigens.slice(1, 0, 4); var pc1Vector eigenVectors.slice(0, 0, 1); var pc1Image image.select([NDVI, WET, NDBSI, LST]) .toArray() .matrixMultiply(pc1Vector.transpose()) .arrayProject([0]) .arrayFlatten([[PC1]]) .rename(PC1);这段代码的要点是centeredCovariance计算的是去中心化后的协方差矩阵eigen函数返回的特征向量第一行对应第一主成分的载荷。跑完之后同样要检查载荷方向再决定是否需要取1减PC1。7. 实操总结四个指数计算中我踩过的坑7.1 云掩膜没做干净LST和NDVI同时失真有一年在GEE里算长江中游某县的RSEI结果图上出现了一大片极高值位置恰好跟当天的薄云区域重合。原因就是云的低温让LST偏低薄云区域的NDVI又因云层反射而偏高两个“正向”信号叠加后PC1在那个区域被顶到了接近1看起来就像生态极好实际上是被云骗了。从那以后我所有的时间序列RSEI在计算前都会用QA_PIXEL做云和云影掩膜掩膜后的像元不参与任何统计。7.2 水体不掩膜PCA载荷方向可能直接被带偏水体在WET上是极高值在LST上是极低值如果研究区内水体面积超过5%水体的强信号会主导协方差矩阵导致PC1从“生态综合指标”变成“水体识别器”。最明显的特征是NDVI在第一主成分上的载荷突然变成负值跟理论预期完全相反。遇到这种情况第一步不是怀疑公式而是把MNDWI大于0的区域掩掉再看载荷方向是否回到正轨。7.3 Landsat 7条带区要不要保留Landsat 7的SLC-off条带是历史问题如果早期年份只能用Landsat 7条带区域的像元值会是空值或0。这类像元在计算NDVI和WET时会表现为异常低值在PCA里又会形成一条条规则的线性噪声。我的建议是做RSEI长期趋势分析时尽量避开Landsat 7优先用Landsat 5补历史、Landsat 8/9补近期如果实在绕不开至少把条带区在最终结果里标为NoData不要参与统计平均。7.4 归一化区间取全图极值结果被两个“极端像元”绑架有一版RSEI结果整幅图看起来灰蒙蒙一片空间差异很小。排查到最后发现某一个像元的WET值异常大很可能来自云影边缘把归一化的分母拉得极大导致其他所有像元被压到0.2到0.6的窄区间里。改成2%到98%分位数截断后影像的层次感立刻出来了。现在我做任何RSEI实验归一化区间都会明确写进方法论不会含糊地用“min-max”。7.5 多年对比项目要统一所有参数做2000年到2020年的RSEI变化分析最容易出的问题就是不同年份用了不同的处理参数。比如2000年用的是Landsat 52020年用的是Landsat 8WET系数变了LST反演的大气参数也不一样连归一化区间都可能不同。这些参数如果不统一最后算出来的变化趋势可能就是参数差异造成的假象。应对办法是在项目启动时就把各年份的处理参数写死成一张表Landsat 5的年份全部用同一组系数、Landsat 8的年份用另一组LST尽量全部用L2产品或者全部用统一流程自己反演归一化区间统一用5%到95%分位数。7.6 一个提高效率的小技巧如果研究区范围大或者要做多个年份的RSEI建议直接在GEE里搭一个批量处理函数把“计算四个指数 归一化 PCA”封装成一个函数输入Image对象输出RSEI影像。这样一年一景的影像只需要循环调用就行不需要下载到本地反复操作。我自己现在的主力流程就是GEE批量算RSEIENVI只用来做小范围精细验证两者结合效率最高。算RSEI这几年我最大的体会是四个指数的计算表面上就是套公式难度不大真正拉开差距的是你对每个公式背后的物理含义、每个参数适用条件的理解。NDVI不能拿错波段WET不能混用系数LST不能跳过大气参数归一化不能糊弄区间PCA不能忽略载荷方向——这些细节单独看都很小叠在一起就决定了结果能不能用。如果看完这篇还是要自己跑一遍建议先拿一景小范围的影像、选一个年份把全流程从NDVI到RSEI完整走一遍每个中间结果都用统计值和图层目视两个方式检查一遍。跑通一次之后后面再换研究区、换年份就只是重复劳动加微调参数的事。最后再分享一个习惯每次跑完四个指数我都会随手导出一张四联图把NDVI、WET、NDBSI、LST并列放在一起检查。这张图能在五分钟内帮我看出80%的波段错位问题推荐你也试试。