RSEI遥感生态指数计算——基于GEE与Landsat8/9卫星
汉祚历传数将终，
天下半壁归曹公。
识时务者为俊杰，
人力岂能抗天工。
目录
1.前言
2.数据准备
2.1.指数介绍
2.2.计算思路
2.3.代码
3.结果展示
1.前言
距离上次公众号发表已经有一年之久，InVEST模型产水量的部分还没有讲完，主要还是由于时间原因一直没有更新，临近期末的最后两周还是有些空余时间，讲讲RSEI（remote sensing based ecological index）生态遥感指数，InVEST模型B站已有教学视频，跟着学问题不大，但不得不说这个模型生态指标（Ecological Indicators）上发文没有1万，也有8千。开源的模型流通性确实强。诸如生态廊道，生态源地之类的生态安全格局更是老生常谈。
RSEI这个指数已经有很久了，最早是福州大学徐涵秋教授(徐涵秋, 2013a, 2013b)在2013年发表的两篇文章中首次提出，至今还在应用可见影响力不小。但在操作时其网上的疑难杂症不少：（1）RSEI计算判断的问题，用PC1计算RSEI0还是用1-PC1去计算得到RSEI0？（2）RSEI计算的详细步骤是什么，如何从f（NDVI,Wet,LST,NDBSI）中得出RSEI。这篇文章希望能尽量解答上述以及相关问题。
本文是基于GEE和Landsat8/9卫星来处理，所以像Landsat 5或者Modis卫星上在公式上就不适用，但是原理都是一样的，GEE能够快速地获取数据而Landsat8/9能够提供较高时空分辨率并且波段全的数据。
2.数据准备
所有的操作在GEE上完成。
2.1.指数介绍
RSEI需要用到四个指数：NDVI,Wet,LST,NDBSI(Table 1)。NDVI、Wet、NDBSI这三个指数都可以在Landsat卫星用波段得出，而LST可以由Landsat卫星的TIR热红外波段得到。这里在GEE处理时得注意两个问题：
（1）我这里用的是Landsat 8/9的C2L2数据，C2L2尽管经过了大气校正得到了地表反射率，但储存时用了缩放因子，NDVI和LST在计算时一定要通过缩放因子和偏移量计算后得到，才是正确的数值，详细见USGS官网说明： （usgs.gov/media/files/landsat-8-9-collection-2-level-2-science-product-guide)。同时注意在“地理空间数据云”上下载的C2L2数据一样也得通过缩放后才能计算NDVI等指数，直接在ArcGIS中用栅格计算器算也是不对的。
（2）Wet计算：Wet的计算实际上是通过Tasseled Cap Transformation (TCT)缨帽变换得到(Crist, 1985; Qingsheng Liu et al., 2014)。什么是缨帽变化这里得解释下，TCT是一种特殊的主成分变换，通过线性变换将原始波段数据转换为新的特征空间，变换系数是基于大量实验数据统计得到的固定系数。主成分变化也就是通过各个解释变量（X）降维处理，找出贡献最大的变量，然后写成y=aX1+bX2+cX3+…+nXn的线性形式。而这里的a，b一直到n这些系数就是通过大量的实验得到。
这里要注意的主要问题是，不同卫星不同传感器上的Wet的计算方式不一样，Landsat7ETM+和Landsat8/9OLI传感器的不一样，modis也不一样。这个得看具体得传感器和卫星。同时值得注意的是TCT不同季节也有差异，而计算Wet用的的是固定系数而没有考虑不同季节下的变化，这一点会对计算出的RSEI不利。
Table 1. 用于RSEI计算的四个指数
| | 绿度指数 | 湿度指数 | 热指数 | 干燥指数 |
| 遥感指数 | NDVI | Wet | LST | NDBSI |
| References | (Hu and Xu, 2018) | (Zhang et al., 2024) | | (Hu and Xu, 2018) |
NDVI=(ρNIR-ρRed)/(ρNIR+ρRed) （1）
Wet=0.1511ρblue+0.1973ρgreen+0.3283ρred +0.3407ρnir-0.7117ρswir1-0.4559ρswir2（2）
NDBSI=(IBI+SI)/2 （3）
IBI={2ρSWI1/(ρSWI1+ρNIR)-(ρNIR/(ρNIR+ρRed)+ρGreen/(ρGreen+ρSWI1)]&/2ρSWI1/(ρSWI1+ρNIR)+[ρNIR/(ρNIR+ρRed)+ρGreen/(ρGreen+ρSWI1)]} （4）
SI=[(ρSWI1+ρRed)-(ρNIR+ρBlue)]/[(ρSWI1+ρRed)+(ρNIR+ρBlue)] （5）
2.2.计算思路
RSEI的计算思路实际上很简单，参照(Fig. 1)：第一步从Landsat8/9中获得NDBI，Wet ，LST，NDBSI四个指标。对四个栅格做标准化处理，减少不同数量级上带来的误差。
第二步进行PCA主成分分析得到PC1,PC2,PC3,PC4四个结果(Table 2)。这四个结果的贡献度总和为1（100%），PC1中贡献度可以看到最大44.86%。
第三步判断PC1中的NDBI，Wet ，LST，NDBSI的正负情况，(Zheng et al., 2022)这篇文章实际上已经对具体计算的细节讲的很清楚。我们先从结果分析我们得到来的RSEI越高，反应得生态环境越好，既然是这样那么NDVI和Wet越是在正方向，越是促进，而LST和NDBSI应该是在负方向削弱。这样我们得出的RSEI才符合观察的标准。从(Table 2)中的PC1可以看出NDBI，Wet同为正（+）而LST，NDBSI同为负（-），所以这里就可以直接使用PC1，而我们来看看PC3，和PC1的正负情况相反，如果出现这种情况，那么就是用1-PC1来得到初始的RSEI也就是RSEI0。这里就要可以说出在引言部分的第二个问题：详细计算的步骤是什么？
我们前者已经说过TCT本来就是一种PCA主成分分析，是一种线性表达，故我们要计算的RSEI0实际上就是RSEI0=0.40025*NDVI-0.62451*LST+ 0.670142* WET-0.0265*NDBSI 。从这里我想各位应该就不难发现，系数就是PC1(Table 2)中的结果，到此上述的两个疑难全部得以解释。何时使用PC1或者1-PC1以及具体的计算原理和步骤。
Fig. 1.计算步骤
Table 2.主成分分析指标
| Indicator | PC1 | PC2 | PC3 | PC4 |
| NDVI | 0.40025 | 0.135958 | -0.88467 | -0.19668 |
| LST | -0.62451 | 0.035449 | -0.42287 | 0.65568 |
| WET | 0.670142 | -0.00901 | 0.139757 | 0.728901 |
| NDBSI | -0.0265 | 0.990039 | 0.1379 | 0.010166 |
| Feature value | 1.794494 | 1.002018 | 0.870283 | 0.333207 |
| Contribution rate（%） | 44.86233 | 25.05044 | 21.75706 | 8.330173 |
2.3.代码
代码说明：注意我这里用的是年尺度的RSEI计算，故采用的均值合成而不是中值合成。
javascript
// 导入研究区
var geometry = ee.FeatureCollection("users/RaySpaniare/jiujiang");
// 加载 Landsat 8/9 数据集
var landsat = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') // Landsat 8
.merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')) // Landsat 9
.filterDate('2023-01-01', '2023-12-31')
.filterBounds(geometry)
.filter(ee.Filter.lt('CLOUD_COVER', 1));
// 计算各个指数
function calculateIndices(image) {
// 应用缩放因子
var optical = image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])
.multiply(0.0000275).add(-0.2); // 光学波段的缩放
var thermal = image.select(['ST_B10'])
.multiply(0.00341802).add(149.0); // 热红外波段的缩放
// 重命名波段以便计算
optical = optical.select(
['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2']
);
// 湿度指数 - Tasseled Cap Wetness for Landsat 8/9 OLI
var wetness = optical.expression(
'Blue * 0.1511 + Green * 0.1973 + Red * 0.3283 + NIR * 0.3407 + SWIR1 * (-0.7117) + SWIR2 * (-0.4559)',
{
'Blue': optical.select('Blue'),
'Green': optical.select('Green'),
'Red': optical.select('Red'),
'NIR': optical.select('NIR'),
'SWIR1': optical.select('SWIR1'),
'SWIR2': optical.select('SWIR2')
}
).rename('wetness');
// 绿度指数 - NDVI (-1 到 1)
var ndvi = optical.normalizedDifference(['NIR', 'Red'])
.rename('greenness');
// 干度指数 - NDBSI
// 1. 计算 SI (Soil Index)
var si = optical.expression(
'((SWIR1 + Red) - (NIR + Blue)) / ((SWIR1 + Red) + (NIR + Blue))',
{
'SWIR1': optical.select('SWIR1'),
'Red': optical.select('Red'),
'NIR': optical.select('NIR'),
'Blue': optical.select('Blue')
}
);
// 2. 计算 IBI (Index-based Built-up Index)
var ibi = optical.expression(
'(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR + Red) + Green / (Green + SWIR1))) / ' +
'(2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR + Red) + Green / (Green + SWIR1)))',
{
'SWIR1': optical.select('SWIR1'),
'NIR': optical.select('NIR'),
'Red': optical.select('Red'),
'Green': optical.select('Green')
}
);
// 3. 计算 NDBSI (Normalized Difference Bare Soil Index)
var ndbsi = si.add(ibi).divide(2)
.rename('dryness');
// 热度指数 - 使用地表温度
var lst = thermal.select(['ST_B10'])
.subtract(273.15) // 转换为摄氏度
.rename('heat');
return image.addBands([wetness, ndvi, ndbsi, lst]);
}
// 对影像集应用指数计算并获取均值影像
var withIndices = landsat.map(calculateIndices);
// 分别获取四个指标的均值影像
var greenness = withIndices.select('greenness').mean();
var wetness = withIndices.select('wetness').mean();
var heat = withIndices.select('heat').mean();
var dryness = withIndices.select('dryness').mean();
// 直接使用研究区实际边界
var region = geometry;
// 设置投影信息和计算参数
var projection = 'EPSG:32648'; // WGS84 UTM Zone 48N
var scale = 30; // Landsat 的空间分辨率
// 将四个指标标准化用于PCA计算
function standardizeForPCA(image) {
var stats = image.reduceRegion({
reducer: ee.Reducer.minMax(),
geometry: region,
scale: scale,
maxPixels: 1e13
});
var min = ee.Number(stats.values().get(0)); // 获取最小值
var max = ee.Number(stats.values().get(1)); // 获取最大值
return image.subtract(min).divide(max.subtract(min));
}
// 标准化四个指标
var greenness_std = standardizeForPCA(greenness);
var wetness_std = standardizeForPCA(wetness);
var heat_std = standardizeForPCA(heat);
var dryness_std = standardizeForPCA(dryness);
// 组合标准化后的指标用于PCA
var compositeImage = ee.Image.cat([greenness_std, wetness_std, heat_std, dryness_std]);
// PCA 计算函数
function calculatePCA(image) {
var scale = 30;
var bandNames = image.bandNames();
// 计算均值
var meanDict = image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: region,
scale: scale,
maxPixels: 1e13
});
// 中心化
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);
// 转换为数组
var arrays = centered.toArray();
// 计算协方差矩阵
var covar = arrays.reduceRegion({
reducer: ee.Reducer.centeredCovariance(),
geometry: region,
scale: scale,
maxPixels: 1e13
});
// 特征值分解
var covarArray = ee.Array(covar.get('array'));
var eigens = covarArray.eigen();
// 获取特征向量
var eigenVectors = eigens.slice(1, 1);
// 获取PC1的系数（第一列）
var ndviCoef = eigenVectors.get([0, 0]);
var lstCoef = eigenVectors.get([1, 0]);
var wetCoef = eigenVectors.get([2, 0]);
var ndbsiCoef = eigenVectors.get([3, 0]);
// 创建调整矩阵
var adjustMatrix = ee.Array([
[ee.Number(ndviCoef).lt(0).multiply(2).subtract(1), 0, 0, 0],
[0, ee.Number(lstCoef).gt(0).multiply(2).subtract(1), 0, 0],
[0, 0, ee.Number(wetCoef).lt(0).multiply(2).subtract(1), 0],
[0, 0, 0, ee.Number(ndbsiCoef).gt(0).multiply(2).subtract(1)]
]);
// 调整特征向量
var adjustedEigenVectors = eigenVectors.matrixMultiply(adjustMatrix);
// 将调整后的特征向量转换为图像格式
var eigenImage = ee.Image(adjustedEigenVectors);
// 计算主成分
var arrayImage = arrays.toArray(1);
var principalComponents = eigenImage.matrixMultiply(arrayImage);
return principalComponents
.arrayProject([0])
.arrayFlatten([['PC1', 'PC2', 'PC3', 'PC4']]);
}
// 计算主成分（使用标准化后的数据）
var principalComponents = calculatePCA(compositeImage);
// 计算 RSEI₀（直接使用PC1，因为已经调整了方向）
var rsei0 = principalComponents.select('PC1');
// 标准化到0-1范围
var rsei = rsei0.unitScale(
rsei0.reduceRegion({
reducer: ee.Reducer.min(),
geometry: region,
scale: scale,
maxPixels: 1e13
}).values().get(0),
rsei0.reduceRegion({
reducer: ee.Reducer.max(),
geometry: region,
scale: scale,
maxPixels: 1e13
}).values().get(0)
);
// 导出四个原始指标
Export.image.toDrive({
image: greenness.clip(region),
description: 'RSEI_greenness_raw',
folder: 'RSEI_Results',
scale: scale,
crs: projection,
region: region.geometry().bounds(),
maxPixels: 1e13
});
Export.image.toDrive({
image: wetness.clip(region),
description: 'RSEI_wetness_raw',
folder: 'RSEI_Results',
scale: scale,
crs: projection,
region: region.geometry().bounds(),
maxPixels: 1e13
});
Export.image.toDrive({
image: heat.clip(region),
description: 'RSEI_heat_raw',
folder: 'RSEI_Results',
scale: scale,
crs: projection,
region: region.geometry().bounds(),
maxPixels: 1e13
});
Export.image.toDrive({
image: dryness.clip(region),
description: 'RSEI_dryness_raw',
folder: 'RSEI_Results',
scale: scale,
crs: projection,
region: region.geometry().bounds(),
maxPixels: 1e13
});
// 导出 RSEI 结果（原有代码）
Export.image.toDrive({
image: rsei.clip(region),
description: 'RSEI_final',
folder: 'RSEI_Results',
scale: scale,
crs: projection,
region: region.geometry().bounds(),
maxPixels: 1e13
});
// 计算湿度指数的实际范围用于显示
wetness.reduceRegion({
reducer: ee.Reducer.minMax(),
geometry: region,
scale: scale,
maxPixels: 1e13
}).evaluate(function(result) {
// 使用计算得到的实际范围来显示湿度
Map.addLayer(wetness.clip(region), {
min: result.wetness_min,
max: result.wetness_max,
palette: ['#FFE4B5', '#0000FF']
}, 'Wetness');
});
// 其他图层显示保持不变
Map.addLayer(greenness.clip(region), {
min: -1,
max: 1,
palette: ['white', 'green']
}, 'Greenness (NDVI)');
Map.addLayer(heat.clip(region), {
min: 0,
max: 50,
palette: ['blue', 'yellow', 'red']
}, 'Heat (LST °C)');
Map.addLayer(dryness.clip(region), {
min: -1,
max: 1,
palette: ['#006400', '#8B4513']
}, 'Dryness (NDBSI)');
//显示RSEI
Map.addLayer(rsei.clip(region), {
min: 0,
max: 1,
palette: ['red', 'yellow', 'green']
}, 'RSEI');
// 设置地图视图范围
Map.centerObject(region, 9);
|
3.结果展示
Fig. 2. 计算得出的四个因子真实值以及RSEI指数（a）NDVI(b)Wet（c）LST（d）NDBSI
(Fig. 2)为RSEI的计算结果和四个变量的真实值。到此RSEI的操作全部结束。
RaySpaniare 2024.12.21于雁栖湖，是日冬至
References
Crist, E.P., 1985. A TM Tasseled Cap equivalent transformation for reflectance factor data. Remote Sensing of Environment 17, 301–306. https://doi.org/10.1016/0034-4257(85)90102-6
Hu, X., Xu, H., 2018. A new remote sensing index for assessing the spatial heterogeneity in urban ecological quality: A case from Fuzhou City, China. Ecological Indicators 89, 11–21. https://doi.org/10.1016/j.ecolind.2018.02.006
Qingsheng Liu, Gaohuan Liu, Chong Huang, Suhong Liu, Jun Zhao, 2014. A tasseled cap transformation for Landsat 8 OLI TOA reflectance images, in: 2014 IEEE Geoscience and Remote Sensing Symposium. Presented at the IGARSS 2014 - 2014 IEEE International Geoscience and Remote Sensing Symposium, IEEE, Quebec City, QC, pp. 541–544. https://doi.org/10.1109/IGARSS.2014.6946479
Zhang, H., Ma, C., Liu, P., 2024. Dynamic evaluation of the ecological evolution and quality of arid and semi-arid deserts in the Aibugai River Basin based on an improved remote sensing ecological index. Ecological Informatics 82, 102727. https://doi.org/10.1016/j.ecoinf.2024.102727
Zheng, Z., Wu, Z., Chen, Y., Guo, C., Marinello, F., 2022. Instability of remote sensing based ecological index (RSEI) and its improvement for time series analysis. Science of The Total Environment 814, 152595. https://doi.org/10.1016/j.scitotenv.2021.152595
徐涵秋, 2013a. 区域生态环境变化的遥感评价指数. 中国环境科学 33, 889–897.
徐涵秋, 2013b. 城市遥感生态指数的创建及其应用. 生态学报 33, 7853–7862.
扫一扫
4865
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
