遥感生态指数（RSEI）计算教程：绿度、湿度、热度、干度四指标合成

1. 引言：为什么需要RSEI

区域生态环境质量评估长期面临一个核心矛盾：单一指标（如植被覆盖度、地表温度）只能反映生态系统的某个侧面，而综合指标体系又往往面临权重赋值主观、量纲不统一、跨区域可比性差等问题。2013年，福州大学徐涵秋教授在《生态学报》发表论文，提出完全基于遥感信息、以主成分分析（PCA）自动确定权重的遥感生态指数（RSEI），为这一困境提供了优雅的解决方案。

RSEI的核心思想可以概括为一句话：用四个与人类生存最密切相关的自然要素——绿度、湿度、热度、干度——来表征生态系统的综合状态，并通过PCA让数据自己"说话"，避免人为赋权的主观性。这一设计思路与联合国可持续发展目标（SDGs）中关于陆地生态、气候行动、清洁水源的多个目标高度契合，也与"山水林田湖草沙"生命共同体理念形成呼应。

本文评述：RSEI之所以在十余年间被广泛应用于城市生态、矿区修复、流域治理、自然保护区评估等场景，根本原因在于它找到了"综合性"与"客观性"之间的平衡点。但笔者认为，RSEI并非万能钥匙——它对数据质量、季节一致性、PCA符号稳定性都有较高要求，工程实践中稍有不慎就可能得出与实际情况相悖的结论。本文的写作目标，正是把这条从原始影像到可信结论的完整链路拆解清楚。

2. RSEI理论基础与四维指标体系

2.1 生态系统的四维表征逻辑

RSEI选择绿度、湿度、热度、干度四个维度，并非随意拼凑，而是基于生态系统对人类生存环境影响的优先级排序。绿度代表植被生产力与碳汇能力，湿度反映水循环与土壤墒情，热度关联城市热岛与能量平衡，干度则指向土地退化与荒漠化风险。这四个维度恰好覆盖了生态系统的"生产—调节—支撑"三大功能。

从信息论角度看，四个指标之间存在一定的相关性（例如植被茂盛区域通常湿度较高、温度较低），但并非完全冗余。PCA的作用正是在保留主要信息的前提下，将这些相关指标压缩为一个综合指数。本文评述：笔者认为，四维指标的选择体现了"最小充分集"的思想——用尽可能少的指标覆盖尽可能多的生态信息维度，这一点比许多动辄十几个指标的综合评价体系更为务实。

2.2 RSEI的数学定义

RSEI的标准计算流程如下：首先对四个指标进行归一化处理，消除量纲差异；然后对归一化后的四波段影像进行主成分变换；取第一主成分（PC1）作为初始生态指数；最后对PC1进行归一化，得到取值在[0,1]区间的RSEI值。RSEI越接近1，表示生态质量越好。

其中f(·)表示归一化处理，PCA表示主成分分析。需要特别注意的是，当PC1对绿度和湿度的载荷为负、对热度和干度的载荷为正时，需要对PC1取反，确保RSEI值越大生态越好。这一符号判断是工程实现中最容易出错的环节之一。

3. 绿度指标：NDVI及其变体

3.1 NDVI的计算与物理意义

归一化植被指数（NDVI）是RSEI中绿度分量的标准选择。其计算公式为：

NDVI利用植被在近红外波段的强反射和红光波段的强吸收特性，有效表征植被覆盖度与生长状态。Landsat 8 OLI传感器对应波段为Band 5（NIR）和Band 4（Red），Landsat 9延续了相同波段设置，Sentinel-2则对应Band 8（NIR）和Band 4（Red）。

3.2 NDVI的局限性与替代方案

NDVI在高植被覆盖区容易饱和，在低覆盖区受土壤背景影响较大。工程实践中，可考虑以下替代或补充方案：

• EVI（增强植被指数）：引入蓝波段进行大气校正，减少气溶胶影响，适合高密度植被区
• SAVI（土壤调整植被指数）：引入土壤调节因子L，适合稀疏植被区
• kNDVI（核归一化植被指数）：2021年Camps-Valls等提出，用核函数替代比值形式，对饱和问题有显著改善

本文评述：笔者认为，对于大多数RSEI应用场景，NDVI仍然是性价比最高的选择——计算简单、物理意义明确、与Landsat/Sentinel系列兼容性好。但在干旱区或高密度植被区，建议进行NDVI与EVI的对比实验，选择与研究区匹配度更高的指标。

4. 湿度指标：WET与缨帽变换

4.1 缨帽变换原理

湿度分量（WET）通常通过缨帽变换（Tasseled Cap Transformation, TCT）提取。缨帽变换是一种线性正交变换，将原始多光谱波段投影到几个具有明确物理意义的特征空间：亮度（Brightness）、绿度（Greenness）和湿度（Wetness）。其中湿度分量反映了土壤水分和植被含水量的综合信息。

不同传感器的缨帽变换系数不同。Landsat 8 OLI的湿度分量计算公式（Baig等，2014）为：

其中B1至B7分别对应Landsat 8 OLI的沿海/气溶胶、蓝、绿、红、近红外、短波红外1、短波红外2波段。注意该公式要求输入反射率数据（已进行辐射定标和大气校正）。

4.2 Sentinel-2的湿度计算

Sentinel-2 MSI传感器没有直接的缨帽变换系数，但可以通过波段对应关系进行近似计算，或使用Shi和Xu（2019）提出的Sentinel-2缨帽变换系数。工程实践中，更稳妥的做法是将Sentinel-2波段重采样至Landsat波段宽度后套用Landsat系数，或直接使用NDWI（归一化水体指数）作为湿度代理指标。

本文评述：笔者认为，缨帽变换的湿度分量虽然物理意义明确，但其系数依赖于传感器特性，跨传感器使用时必须格外谨慎。在多源数据融合场景下，建议优先统一到同一传感器体系，或采用NDWI/MNDWI等对传感器差异不敏感的指数作为替代。

5. 热度指标：地表温度LST反演

5.1 单通道算法

地表温度（LST）是RSEI中热度分量的核心指标。对于Landsat系列数据，常用的反演方法包括辐射传输方程法、单窗算法和单通道算法。其中Jiménez-Muñoz和Sobrino（2003）提出的单通道算法（SC算法）应用最为广泛。

SC算法的核心步骤包括：首先将热红外波段的DN值转换为辐射亮度，然后利用大气水汽含量计算大气透过率，最后结合地表比辐射率计算地表温度。对于Landsat 8 TIRS Band 10，公式可简化为：

其中γ和δ为与亮温相关的中间变量，ψ₁、ψ₂、ψ₃为大气函数，ε为地表比辐射率，L_sensor为传感器处辐射亮度。

5.2 劈窗算法与温度产品

对于Landsat 8/9 TIRS的两个热红外波段，可采用劈窗算法（Split-Window）提高反演精度。此外，USGS提供的Landsat Collection 2 Level-2产品已包含地表温度波段（ST_B10），经过大气校正，可直接使用，大幅降低工程实现难度。

本文评述：笔者认为，对于RSEI应用而言，Landsat Collection 2 Level-2的ST_B10产品已经足够——其精度在大多数场景下优于自行反演的结果，且省去了大气参数获取的麻烦。如果研究区跨越多个Landsat行列号，建议优先使用GEE中的LANDSAT/LC08/C02/T1_L2数据集，其ST_B10波段已做统一处理。

6. 干度指标：NDBSI构建

6.1 裸土指数与建筑指数

干度指标（NDBSI）由裸土指数（SI）和建筑指数（IBI）取平均得到，分别反映自然裸土和人工不透水面的干燥程度。

裸土指数SI的计算公式（Rikimaru等，2002）为：

建筑指数IBI的计算公式（Xu，2008）为：

最终NDBSI = (SI + IBI) / 2。

6.2 干度指标的工程考量

NDBSI对城市扩张和土地退化极为敏感，是RSEI中变化最剧烈的分量之一。在矿区、城市边缘带等区域，NDBSI的上升往往先于NDVI的下降，具有早期预警价值。

本文评述：笔者认为，NDBSI的计算涉及多个波段比值，对大气校正精度和云掩膜质量要求较高。工程实践中，建议对NDBSI结果进行3×3中值滤波，消除孤立异常像元，同时结合研究区特点调整SI和IBI的权重（默认等权，但在纯自然区域可适当降低IBI权重）。

7. 主成分分析与RSEI合成

7.1 PCA的数学本质

主成分分析是一种正交线性变换，将原始变量投影到方差最大的方向上。对于RSEI而言，PCA的作用有三：降维、去相关、自动赋权。第一主成分（PC1）通常解释了原始四个指标60%—80%的方差，是生态信息的最大公约数。

PCA的数学推导可简述为：对归一化后的数据矩阵X（n×4），计算协方差矩阵C = XᵀX/(n-1)，对C进行特征分解，取最大特征值对应的特征向量作为PC1的载荷向量。

7.2 符号判断与归一化

PCA结果的一个关键工程问题是符号不确定性：特征向量的方向可以整体翻转而不影响方差解释率。因此，必须根据PC1对四个指标的载荷符号来判断是否需要取反。判断规则是：如果PC1对NDVI和WET的载荷为负，对LST和NDBSI的载荷为正，则PC1取反。

归一化采用极差标准化：RSEI = (PC1 - PC1_min) / (PC1_max - PC1_min)。归一化范围通常取研究区内的2%—98%分位数，以避免极端值影响。

本文评述：笔者认为，PCA符号判断是RSEI计算中最容易出错的环节。建议在代码中显式打印PC1的载荷向量，并设置自动检查逻辑——如果载荷方向与预期不符，自动取反并输出警告信息。此外，对于多时相分析，应确保各时相PCA的符号方向一致，否则时序对比将失去意义。

8. Google Earth Engine实现

8.1 数据准备与预处理

GEE平台提供了Landsat和Sentinel系列数据的在线访问与计算能力，是RSEI工程实现的首选平台。以下代码展示了基于Landsat 8 Collection 2 Level-2数据的RSEI计算核心流程。

// 定义研究区与时间范围
var roi = ee.Geometry.Rectangle([116.0, 39.5, 117.0, 40.5]);
var startDate = '2023-06-01';
var endDate = '2023-09-30';

// 加载Landsat 8 Collection 2 Level-2数据
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterBounds(roi)
  .filterDate(startDate, endDate)
  .filter(ee.Filter.lt('CLOUD_COVER', 20));

// 云掩膜函数
function maskL8(image) {
  var qa = image.select('QA_PIXEL');
  var cloud = qa.bitwiseAnd(1 << 3).eq(0);
  var shadow = qa.bitwiseAnd(1 << 4).eq(0);
  return image.updateMask(cloud).updateMask(shadow);
}

// 中值合成
var composite = l8.map(maskL8).median().clip(roi);

// 反射率缩放（Collection 2 Level-2需乘以0.0000275再减0.2）
var opticalBands = composite.select(['SR_B.*']).multiply(0.0000275).add(-0.2);
var thermalBand = composite.select('ST_B10').multiply(0.00341802).add(149.0);

8.2 四指标计算

// NDVI计算
var ndvi = opticalBands.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');

// WET计算（Landsat 8 OLI缨帽变换湿度分量）
var wet = opticalBands.expression(
  '0.1511*B1 + 0.1973*B2 + 0.3283*B3 + 0.3407*B4 - 0.7117*B5 - 0.4559*B7', {
    'B1': opticalBands.select('SR_B1'),
    'B2': opticalBands.select('SR_B2'),
    'B3': opticalBands.select('SR_B3'),
    'B4': opticalBands.select('SR_B4'),
    'B5': opticalBands.select('SR_B5'),
    'B7': opticalBands.select('SR_B7')
  }).rename('WET');

// LST计算（Collection 2 Level-2 ST_B10已为开尔文温度）
var lst = thermalBand.rename('LST');

// NDBSI计算
var si = opticalBands.expression(
  '((SWIR1 + Red) - (NIR + Blue)) / ((SWIR1 + Red) + (NIR + Blue))', {
    'SWIR1': opticalBands.select('SR_B6'),
    'Red': opticalBands.select('SR_B4'),
    'NIR': opticalBands.select('SR_B5'),
    'Blue': opticalBands.select('SR_B2')
  });

var ibi = opticalBands.expression(
  '(2*SWIR1/(SWIR1+NIR) - (NIR/(NIR+Red) + Green/(Green+SWIR1))) / ' +
  '(2*SWIR1/(SWIR1+NIR) + (NIR/(NIR+Red) + Green/(Green+SWIR1)))', {
    'SWIR1': opticalBands.select('SR_B6'),
    'NIR': opticalBands.select('SR_B5'),
    'Red': opticalBands.select('SR_B4'),
    'Green': opticalBands.select('SR_B3')
  });

var ndbsi = si.add(ibi).divide(2).rename('NDBSI');

8.3 PCA合成与RSEI输出

// 归一化函数（基于研究区2%-98%分位数）
function normalize(image, region) {
  var percentiles = image.reduceRegion({
    reducer: ee.Reducer.percentile([2, 98]),
    geometry: region,
    scale: 30,
    maxPixels: 1e13
  });
  var min = ee.Number(percentiles.values().get(0));
  var max = ee.Number(percentiles.values().get(1));
  return image.subtract(min).divide(max.subtract(min));
}

// 四指标归一化
var ndviNorm = normalize(ndvi, roi);
var wetNorm = normalize(wet, roi);
var lstNorm = normalize(lst, roi);
var ndbsiNorm = normalize(ndbsi, roi);

// 组合为四波段影像
var fourBands = ndviNorm.addBands(wetNorm).addBands(lstNorm).addBands(ndbsiNorm);

// PCA变换
var pca = fourBands.toArray().principalComponents(4);
var pc1 = pca.arraySlice(0, 0, 1).arrayProject([0]).arrayFlatten([['PC1']]);

// 符号判断与归一化
var pc1Norm = normalize(pc1, roi);
var rsei = pc1Norm.rename('RSEI');

// 输出
Map.addLayer(rsei, {min: 0, max: 1, palette: ['#8B0000', '#FF4500', '#FFD700', '#90EE90', '#006400']}, 'RSEI');
Export.image.toDrive({
  image: rsei,
  description: 'RSEI_2023',
  region: roi,
  scale: 30,
  maxPixels: 1e13
});

本文评述：笔者认为，GEE实现RSEI的最大优势在于算力与数据获取的便利性，但需要注意两个工程细节：一是PCA变换在GEE中是基于整个研究区计算的，如果研究区跨越多个气候带，建议分区计算；二是归一化的分位数选择会影响RSEI的绝对数值，多时相分析时应保持分位数设置一致。

9. Python本地实现

9.1 技术栈与数据读取

对于需要精细控制或离线处理的场景，Python本地实现是更灵活的选择。推荐技术栈：rasterio用于栅格读写，numpy用于数值计算，scikit-learn用于PCA，matplotlib用于可视化。

import rasterio
import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# 读取Landsat 8波段数据
def read_band(path):
    with rasterio.open(path) as src:
        return src.read(1).astype(np.float32), src.profile

# 读取各波段
b2, profile = read_band('LC08_B2.tif')  # Blue
b3, _ = read_band('LC08_B3.tif')        # Green
b4, _ = read_band('LC08_B4.tif')        # Red
b5, _ = read_band('LC08_B5.tif')        # NIR
b6, _ = read_band('LC08_B6.tif')        # SWIR1
b7, _ = read_band('LC08_B7.tif')        # SWIR2
b10, _ = read_band('LC08_B10.tif')      # Thermal

# 辐射定标与反射率转换（Collection 2 Level-2）
def to_reflectance(dn):
    return dn * 0.0000275 - 0.2

b2, b3, b4, b5, b6, b7 = [to_reflectance(b) for b in [b2, b3, b4, b5, b6, b7]]
lst = b10 * 0.00341802 + 149.0  # 开尔文温度

9.2 四指标计算与PCA

# NDVI
ndvi = (b5 - b4) / (b5 + b4 + 1e-10)

# WET
wet = (0.1511*b2 + 0.1973*b3 + 0.3283*b4 +
       0.3407*b5 - 0.7117*b6 - 0.4559*b7)

# NDBSI
si = ((b6 + b4) - (b5 + b2)) / ((b6 + b4) + (b5 + b2) + 1e-10)
ibi = (2*b6/(b6+b5+1e-10) - (b5/(b5+b4+1e-10) + b3/(b3+b6+1e-10))) / \
      (2*b6/(b6+b5+1e-10) + (b5/(b5+b4+1e-10) + b3/(b3+b6+1e-10)))
ndbsi = (si + ibi) / 2

# 归一化（2%-98%分位数）
def normalize(arr):
    p2, p98 = np.nanpercentile(arr, [2, 98])
    return np.clip((arr - p2) / (p98 - p2 + 1e-10), 0, 1)

ndvi_n = normalize(ndvi)
wet_n = normalize(wet)
lst_n = normalize(lst)
ndbsi_n = normalize(ndbsi)

# 堆叠为四波段
stack = np.stack([ndvi_n, wet_n, lst_n, ndbsi_n], axis=-1)
h, w, c = stack.shape
flat = stack.reshape(-1, c)

# 去除NaN
mask = ~np.isnan(flat).any(axis=1)
valid = flat[mask]

# PCA
pca = PCA(n_components=4)
pc = pca.fit_transform(valid)

# 符号判断
loadings = pca.components_[0]
if loadings[0] < 0 and loadings[1] < 0 and loadings[2] > 0 and loadings[3] > 0:
    pc1 = -pc[:, 0]
else:
    pc1 = pc[:, 0]

# 归一化
rsei_valid = normalize(pc1)

# 还原为栅格
rsei_flat = np.full(flat.shape[0], np.nan)
rsei_flat[mask] = rsei_valid
rsei = rsei_flat.reshape(h, w)

# 保存
profile.update(dtype=rasterio.float32, count=1, compress='lzw')
with rasterio.open('RSEI_2023.tif', 'w', **profile) as dst:
    dst.write(rsei.astype(np.float32), 1)

本文评述：笔者认为，Python实现的最大价值在于可复现性和可定制性。建议将上述代码封装为函数或类，便于批量处理多时相数据。同时，PCA的符号判断逻辑应显式输出载荷向量，便于人工核查。

10. 时序监测与变化检测

10.1 多时相RSEI计算策略

时序RSEI分析的核心挑战在于保证各时相之间的可比性。建议采用以下策略：

• 统一季节窗口：各年份使用相同的月份范围（如6—9月），避免物候差异干扰
• 统一归一化基准：使用所有年份的全局2%—98%分位数，而非逐年计算
• 统一PCA基准：将所有年份数据合并后进行PCA，或逐年PCA后检查载荷一致性
• 统一传感器：尽量使用同一传感器（如全部用Landsat 8），或进行交叉定标

10.2 变化检测与趋势分析

时序RSEI的变化检测常用方法包括：差值法（后一时相减前一时相）、Theil-Sen中位数趋势分析和Mann-Kendall显著性检验。Theil-Sen方法对异常值稳健，适合长时序分析。

import pymannkendall as mk
from scipy.stats import theilslopes

# 假设rsei_series为某像元多年的RSEI值序列
# Theil-Sen斜率
slope, intercept, lo_slope, up_slope = theilslopes(rsei_series, years)

# Mann-Kendall检验
result = mk.original_test(rsei_series)
print(f"趋势: {result.trend}, p值: {result.p:.4f}")

本文评述：笔者认为，时序RSEI分析中最容易被忽视的是"基准漂移"问题——如果各年份独立归一化，RSEI的绝对值将失去可比性，趋势分析可能得出完全错误的结论。建议在工程实践中始终采用全局归一化策略，并在论文中明确说明归一化方法。

11. 应用案例与工程洞察

11.1 城市扩张监测

以某中部省会城市为例（模拟数据，基于Landsat 8 2023年影像），城市核心区RSEI均值约为0.35，城市边缘带约为0.52，远郊农田区约为0.68，山区林地约为0.85。RSEI的空间梯度清晰反映了城市化对生态质量的影响。

值得注意的是，城市公园和绿地的RSEI值可达到0.6—0.7，显著高于周边建成区，说明城市内部生态修复措施的有效性。本文评述：笔者认为，RSEI在城市生态评估中的价值不仅在于"打分"，更在于识别生态质量的热点与冷点，为城市绿地规划提供空间依据。

11.2 矿区生态修复评估

矿区是RSEI应用的重要场景。采矿活动导致NDVI下降、NDBSI上升、LST升高，RSEI显著降低。生态修复后，RSEI逐步回升，但恢复速度因修复方式而异。根据笔者对某矿区（模拟数据）的分析，自然恢复区的RSEI年均回升约0.008，人工复垦区约0.015，说明人工干预可加速生态恢复。

11.3 流域生态质量评估

流域尺度的RSEI分析需考虑上下游差异。上游水源涵养区RSEI通常较高（0.7—0.9），中游农业区中等（0.5—0.7），下游城市区较低（0.3—0.5）。RSEI的沿程变化可用于识别流域生态退化热点，为水土保持和生态补偿提供依据。

12. 前沿趋势与未来展望

12.1 高光谱与多源数据融合

随着高分五号、ZY1-02D等高光谱卫星的投入使用，RSEI的指标构建有望从多光谱扩展到高光谱。高光谱数据可提供更精细的植被生化参数（叶绿素、水分、木质素），使绿度和湿度指标更加精准。此外，SAR数据（Sentinel-1）可穿透云层，为多云雨区的RSEI计算提供补充。

12.2 深度学习与RSEI

近年来，深度学习在遥感领域的应用日益广泛。对于RSEI而言，深度学习可用于：云检测与修复、超分辨率重建、指标预测与时空插值。例如，基于ConvLSTM的RSEI时空预测模型，可利用历史RSEI序列预测未来生态质量变化趋势。

12.3 多尺度耦合与尺度效应

RSEI的尺度效应是一个尚未充分研究的问题。同一区域在30m、100m、250m、1km等不同尺度下计算的RSEI可能存在显著差异。未来研究需要建立RSEI的尺度转换模型，实现多尺度生态评估结果的兼容与融合。

本文评述：笔者认为，RSEI的未来发展应重点关注三个方向：一是与生态系统服务价值评估的耦合，将RSEI从"状态指标"扩展为"功能指标"；二是与气候变化模型的嵌套，评估未来气候情景下的生态质量演变；三是与社交媒体、手机信令等社会感知数据的融合，构建"自然—社会"耦合的生态评估框架。

13. 总结与最佳实践清单

RSEI作为一个成熟且持续演进的遥感生态评估工具，其价值已在大量研究中得到验证。但工具的价值最终取决于使用者的理解深度。笔者认为，掌握RSEI的关键不在于记住公式，而在于理解每个指标背后的生态学含义，以及PCA合成过程中的工程细节。只有把"数据—指标—合成—验证"这条链路走通，RSEI才能真正成为可信的生态评估工具。

14. 主要参考文献