Google Earth Engine——基于改进的RSEI评估生态环境（水体掩膜后）

未经允许，禁止随意转载，尊重他人版权，仅供学习参考，欢迎交流。

背景介绍

   遥感生态指数(Remote Sensing Ecological Index)的获得，是使用主成分分析法耦合了绿度、湿度、干度、热度指标，以1-PC1（PC1，主成分分析结果的第一分量）标准化的结果作为遥感生态指数（徐涵秋，2013）。以上四个指标中，绿度反映了植被覆盖度，是生态环境好坏的重要指标；湿度和干度反映了生态环境中地表与地上的水分含量多少。热度指标即地表温度是区域和全球范围地表物理过程的重要因子，也是研究地表和大气物质交换和能量交换的关键参数。因此，使用主成分分析法耦合这四个指标，避免主观确定权重，为评估生态环境质量提供一种快速有效的方法。

   在使用RSEI模型时，一些学者根据2013年提出的1-PC1来计算RSEI（徐涵秋，2013），但也有其他学者直接使用PC1作为区域生态遥感指数。究其原因，是学者们只是简单地应用模型，对模型的运作和模型机制缺乏了解。由于主成分分析方法中特征向量方向的非唯一性，这两个模型会带来相反的结果。在实际研究中，学者们通常直接使用该方法，很少有人注意到特征向量及其方向，大大限制了遥感生态环境评估方法的发展。同时，由于缺乏对模型机制的研究，盲目地应用和改进模型，有时会误导学者做出错误的评价。

1. 通过研究主成分分析方法中特征向量方向的改变对RSEI的影响，只有当NDVI和Wet的特征向量为负值，NDSI和LST的特征向量为正值时，使用1-PC1计算的RSEI才是正确的。而直接用PC1的模型只有在NDVI和Wet的特征向量为正值，NDSI和LST的特征向量为负值时才能得到正确的结果（Li Ning等，2020），提出了改进模型如式1，这里的Vndvi、Vwet是指NDVI和Wet对PC1的特征向量值。
2. 在上述研究基础上，有学者通过大样本测试了RSEI模型特征向量在时间序列中的演变。结果发现，改变输入波段的顺序都可能导致RSEI是相反的，每个生态因子的特征向量都会影响相应的RSEI的方向。因此提出一个改进的模型，通过选择受季节性变化影响较小的Wet湿度分量来判断第一主成分PC1的特征向量方向。该模型的方向可以根据主观经验自动修改，无需研究人员干预，改进后的模型（如下式）可以适应不同时期的RSEI计算，同时无论输入指标的顺序如何变化，最终结果的方向是正确的（Zhang Yaqiu等，2021）。

   除了主成分分析来进行权重赋值，还有学者利用粒化熵方法来确定RSEI的权重。例如基于Modis数据建立了中国遥感生态指数（RSEI）模型,进行RSEIs的知识粒化，提出了根据指标的特征确定指标知识粒化熵权的方法（(Liao Weihua等，2020），避免了主成分分析特征向量方向不唯一性对RSEI的影响。

这里采用第二种改进方法，对研究区域基于改进的RSEI评估遥感生态环境。

• 首先是对数据和方法的说明。

  使用USGS Landsat 8 Level 2, Collection 2, Tier 1 (google.com)的30m分辨率数据，通过Google Earth Engine——基于新的Landsat SR数据集去云处理，对一年内的每景影像计算四个指标然后中值合成，进行归一化。在指标合成之前，如果研究区域内含有水体，还需要进行水体掩膜。

• 水体掩膜 

   可以根据MNDWI提取水体阈值，也可以通过现有的全球水体数据。

   一是MOD44W.006 Terra Land Water Mask Derived From MODIS and SRTM Yearly Global 250m  |  Earth Engine Data Catalog  |  Google Developers，源自 MODIS 和 SRTM 的陆地Water Mask，全球 250 米分辨率，提供2000年-2015年数据。

//源自 MODIS 和 SRTM 的陆地水体，全球 250 米 //0： 土地 //1： 水
var dataset = ee.ImageCollection('MODIS/006/MOD44W') .filterDate(start_date, end_date) .select('water_mask')

       二是JRC Yearly Water Classification History, v1.3  |  Earth Engine Data Catalog  |  Google Developers，从Landsat 5、7 和 8 获取，全球30米分辨率，实际提供1984年-2019年数据。

//JRC-年度水分类历史-30m  //0 No data //1 Not water //2 Seasonal water //3 Permanent water
var dataset2 = ee.ImageCollection('JRC/GSW1_3/YearlyHistory') .filterDate(start_date, end_date) .filterBounds(roi)

• GEE代码-指标计算+水体掩膜+PCA

/*目录*/
//L8去云
//年度数据预处理
//LST
//NDVI
//Wet
//NDBSI
//归一化-MaxMin
//标准化-中心均值化
//水体掩膜 
//归一化波段，归并
//主成分分析
/*************************************************************************************/
//导入自己的研究区，将其定义为roi
var roi = ee.FeatureCollection("users/自定义 /自定义");  16 var Year='2019';  17 var star_date = Year+'-01-01' //定义起始时间
var end_date = Year+'-12-31' //定义终止时间
var L8_SR = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") //加载L8_SR影像
/*************************************************************************************/
/****************************************L8 去云****************************************/
// 使用Landsat8 Collection 2,Level 2 QA_PIXEL波段（CFMask）去云
function maskL8sr(image) {  24     // Bit 0 - Fill
    // Bit 1 - Dilated Cloud
    // Bit 2 - Cirrus
    // Bit 3 - Cloud
    // Bit 4 - Cloud Shadow
    var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);  30     var saturationMask = image.select('QA_RADSAT').eq(0);  31     // 用缩放后的波段替换，并应用云掩膜
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);  33     var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);  34     // 用缩放后的波段替换，并应用云掩膜
    return image.addBands(opticalBands, null, true).addBands(thermalBands, null, true).updateMask(qaMask).updateMask(saturationMask);  36 }  37 /*************************************************************************************/
/******************************处理一年的数据_L8***************************************/
var collection2 = L8_SR.filterDate(star_date, end_date).filterBounds(roi).map(maskL8sr);  40 var vizParams = {  41     bands: ['SR_B4', 'SR_B3', 'SR_B2'],  42     min: 0,  43     max: 0.3,  44     gamma: [0.95, 1.1, 1]  45 }  46 var L8_median = collection2.median();  47 
/*************************************************************************************/
/***************************************LST_L8****************************************/
//LST直接就是SR的红外波段，单位是开尔文
var Lst_8 = L8_median.select('ST_B10');  52 /*************************************************************************************/
/***************************************NDVI_L8****************************************/
var getNDVI8 = function(image) {  55     var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']);  56     return (ndvi);  57 };  58 var NDVI_8 = collection2.map(getNDVI8).median();  59 
/*************************************************************************************/
/***************************************Wet_L8****************************************/
var getWet8 = function(image) {  63     var wet = image.expression('B*(0.1511) + G*(0.1973) + R*(0.3283) + NIR*(0.3407) + SWIR1*(-0.7117) + SWIR2*(-0.4559)', {  64         'B': image.select(['SR_B2']),  65         'G': image.select(['SR_B3']),  66         'R': image.select(['SR_B4']),  67         'NIR': image.select(['SR_B5']),  68         'SWIR1': image.select(['SR_B6']),  69         'SWIR2': image.select(['SR_B7'])  70  })  71     return (wet);  72 };  73 var Wet_8 = collection2.map(getWet8).median();  74 /*************************************************************************************/
/***************************************NDBSI_L8****************************************/
var getNDBSI8 = function(image) {  77     var ibi = image.expression('(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR + RED) - GREEN / (GREEN + SWIR1))) / (2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1)))', {  78         'SWIR1': image.select('SR_B6'),  79         'NIR': image.select('SR_B5'),  80         'RED': image.select('SR_B4'),  81         'GREEN': image.select('SR_B3')  82  });  83     var si = image.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', {  84         'SWIR1': image.select('SR_B6'),  85         'NIR': image.select('SR_B5'),  86         'RED': image.select('SR_B4'),  87         'BLUE': image.select('SR_B2')  88  });  89     var ndbsi = ibi.subtract(si).divide(2);  90     return (ndbsi);  91 }  92 var NDBSI_8 = collection2.map(getNDBSI8).median();  93 /*************************************************************************************/
/*************************************归一化-MaxMin***********************************/
//归一化函数
var img_norm_1 = function(image) {  97     var minMax = image.reduceRegion({  98  reducer: ee.Reducer.minMax(),  99  geometry: roi, 100         scale: 30, 101         maxPixels: 10e13, 102         //tileScale: 16
 }) 104     var normalize = ee.ImageCollection.fromImages(image.bandNames().map(function(name) { 105         name = ee.String(name); 106         var band = image.select(name); 107         return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max')))); 108  })).toBands().rename(image.bandNames()); 109     return normalize; 110 } 111 /*************************************************************************************/
/***********************************标准化-中心均值化*********************************/
var img_norm_2 = function(image) { 114     var mean_std = image.reduceRegion({ 115         reducer: ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true), 116  geometry: roi, 117         scale: 30, 118         maxPixels: 10e9, 119  }); 120     // use unit scale to normalize the pixel values
    var unitScale = ee.ImageCollection.fromImages(image.bandNames().map(function(name) { 122         name = ee.String(name); 123         var band = image.select(name); 124         var mean = ee.Number(mean_std.get(name.cat('_mean'))); 125         var std = ee.Number(mean_std.get(name.cat('_stdDev'))); 126         var max = mean.add(std.multiply(3)) 127         var min = mean.subtract(std.multiply(3)) 128         var band1 = ee.Image(min).multiply(band.lt(min)).add(ee.Image(max).multiply(band.gt(max))).add(band.multiply(ee.Image(1).subtract(band.lt(min)).subtract(band.gt(max)))) 129         var result_band = band1.subtract(min).divide(max.subtract(min)); 130         return result_band; 131  })).toBands().rename(image.bandNames()); 132     return unitScale 133 } 134 /*************************************************************************************/

/*************************************归一化波段，归并**********************************/
//归一化
var norm_NDVI= img_norm_1(NDVI_8); 139 var norm_Wet=img_norm_1(Wet_8); 140 var norm_NDBSI= img_norm_1(NDBSI_8); 141 var norm_Lst= img_norm_1(Lst_8); 142 
var L8_median=L8_median.addBands(norm_NDVI.rename('NDVI').toFloat()); 144 var L8_median=L8_median.addBands(norm_Wet.rename('Wet').toFloat()); 145 var L8_median=L8_median.addBands(norm_NDBSI.rename('NDBSI').toFloat()); 146 var L8_median=L8_median.addBands(norm_Lst.rename('LST').toFloat()); 147 
//掩膜-MNDWI
// var getMNDWI8 = function(image){
// var mndwi = image.normalizedDifference(['SR_B3', 'SR_B6']);
// return (mndwi);
// }
// var MNDWI_8 = collection2.map(getMNDWI8).median().clip(roi); 
// //Map.addLayer(MNDWI_8);
// //Treshhold-0.2
// //print("MNDWI_8", ui.Chart.image.histogram(MNDWI_8, roi, 100, 258))
// var mask = MNDWI_8.lt(0.2);
/*************************************水体掩膜 **********************************/
//JRC年度水分类历史-30m 
//0 cccccc No data
//1 ffffff Not water
//2 99d9ea Seasonal water
//3 0000ff Permanent water
var dataset2 = ee.ImageCollection('JRC/GSW1_3/YearlyHistory') 165              .filterDate(star_date,end_date).filterBounds(roi).select('waterClass').toBands(); 166 var visualization = { 167   min: 0.0, 168   max: 3.0, 169   palette: ['cccccc', 'ffffff', '99d9ea', '0000ff'] 170 }; 171 Map.addLayer(dataset2,visualization,'water'); 172 var mask0 = dataset2.clip(roi).eq(2); 173 var mask1 = dataset2.clip(roi).eq(3); 174 var mask_ori=mask0.add(mask1).unmask().clip(roi).eq(0); 175 
var L8_median = L8_median.updateMask(mask_ori) 177 /*************************************************************************************/
var bands = ["NDVI","Wet","NDBSI","LST"]; 179 var bandImage =L8_median.select(bands) 180 

/*************************************************************************************/

/****************************************主成分分析************************************/
//Input-使用合成波段
var region = roi; 187 var pre_image = bandImage.select(bands); 188 var scale = 30; 189 var bandNames = pre_image.bandNames(); 190 
// 数据平均
// 标准化拉伸-principal components.
var meanDict = pre_image.reduceRegion({ 194  reducer: ee.Reducer.mean(), 195  geometry: region, 196  scale: scale, 197     maxPixels: 1e9
}); 199 var means = ee.Image.constant(meanDict.values(bandNames)); 200 var centered = pre_image.subtract(means); 201 
// 重命名
var getNewBandNames = function(prefix) { 204   var seq = ee.List.sequence(1, bandNames.length()); 205   return seq.map(function(b) { 206     return ee.String(prefix).cat(ee.Number(b).int()); 207  }); 208 }; 209 
//主成分分析函数
var getPrincipalComponents = function(centered, scale, region) { 212   // 图像转为一维数组
  var arrays = centered.toArray(); 214 
  // 协方差矩阵
  var covar = arrays.reduceRegion({ 217  reducer: ee.Reducer.centeredCovariance(), 218  geometry: region, 219  scale: scale, 220     maxPixels: 1e9
 }); 222 
  // 获取“数组”协方差结果并转换为数组。
  // 波段与波段之间的协方差
  var covarArray = ee.Array(covar.get('array')); 226 
  // 执行特征分析
  var eigens = covarArray.eigen(); 229 
  // 分割特征值
  var eigenValues = eigens.slice(1, 0, 1); 232 
 //计算主成分贡献度
    var eigenValuesList = eigenValues.toList().flatten() 235     var total = eigenValuesList.reduce(ee.Reducer.sum()) 236     var percentageVariance = eigenValuesList.map(function(item) { 237       return (ee.Number(item).divide(total)).multiply(100).format('%.2f') 238  }) 239       print('特征值',eigenValuesList ) 240     print("贡献率（%）", percentageVariance) 241 

  // 分割出特征向量，其特征向量为行
  var eigenVectors = eigens.slice(1, 1); 245    print('eigenVectors',eigenVectors); 246 
  // 将图像转换为二维阵列
  var arrayImage = arrays.toArray(1); 249 
  // 使用特征向量矩阵左乘图像阵列
  var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage); 252 
  // 将特征值的平方根转换为P波段图像。
  var sdImage = ee.Image(eigenValues.sqrt()) 255     .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]); 256 
  // 将PC转换为P波段图像，通过SD标准化。
  return principalComponents 259     // 抛出一个不需要的维度，[[]]->[]。
    .arrayProject([0]) 261     // 使单波段阵列映像成为多波段映像，[]->image。
    .arrayFlatten([getNewBandNames('pc')]) 263     // 通过SDs使PC标准化
 .divide(sdImage); 265 
}; 267 
//进行主成分分析，获得分析结果
var pcImage = getPrincipalComponents(centered, scale, region); 270 
//基于Wet对PC1的特征向量判断，选择合适的RSEI计算公式
//Vwet<0
//var RSEI_0 = L8_median.expression('1-b1',{b1:pcImage.select('pc1')});

//Vwet>=0
var RSEI_0 = L8_median.expression('0+b1',{b1:pcImage.select('pc1')}); 277 
var RSEI = img_norm_1(RSEI_0).clip(roi); 279 
print("RSEI"+Year, ui.Chart.image.histogram(RSEI, roi, 100, 258)) 281 
var visParam = { 283     palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
        '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
 }; 286 //添加图层
 Map.addLayer(RSEI, visParam, 'RSEI_'+Year) 288 /********************************导出 *****************************************************/
Export.image.toDrive({ 290  image: RSEI, 291       description: 'RSEI_'+Year, 292       folder: "Annual_RSEI", 293       scale:30, 294  region:roi, 295       fileFormat:'GeoTIFF', 296       crs: 'EPSG:4326', 297     formatOptions:{cloudOptimized: true}
    });

•  参考文献

[1]徐涵秋. 2013. 区域生态环境变化的遥感评价指数[J]. 中国环境科学,33(05): 889-897.

[2]Li N,Wang J Y, Qin F. 2020. The improvement of ecological environment index model RSEI. Arab J Geosci 13, 403. [DOI 10.1007/s12517-020-05414-7]

[3]Zhang Y Q, Fang J. 2021. Developing a remote sensing-based ecological index based on improved biophysical features. Journal of Applied Remote Sensing 16(1), 012008.[DOI 10.1117/1.JRS.16.012008]

[4]LIAO W H, JIANG W G.2020. Evaluation of the Spatiotemporal Variations in the Eco-environmental Quality in China Based on the Remote Sensing Ecological Index[J]. Remote Sensing,12(15): 2462.