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/自定义 /自定义");
var Year='2019';
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) {
    // 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);
    var saturationMask = image.select('QA_RADSAT').eq(0);
    // 用缩放后的波段替换，并应用云掩膜
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
    var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
    // 用缩放后的波段替换，并应用云掩膜
    return image.addBands(opticalBands, null, true).addBands(thermalBands, null, true).updateMask(qaMask).updateMask(saturationMask);
}
/*************************************************************************************/
/******************************处理一年的数据_L8***************************************/
var collection2 = L8_SR.filterDate(star_date, end_date).filterBounds(roi).map(maskL8sr);
var vizParams = {
    bands: ['SR_B4', 'SR_B3', 'SR_B2'],
    min: 0,
    max: 0.3,
    gamma: [0.95, 1.1, 1]
}
var L8_median = collection2.median();

/*************************************************************************************/
/***************************************LST_L8****************************************/
//LST直接就是SR的红外波段，单位是开尔文
var Lst_8 = L8_median.select('ST_B10');
/*************************************************************************************/
/***************************************NDVI_L8****************************************/
var getNDVI8 = function(image) {
    var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']);
    return (ndvi);
};
var NDVI_8 = collection2.map(getNDVI8).median();

/*************************************************************************************/
/***************************************Wet_L8****************************************/
var getWet8 = function(image) {
    var wet = image.expression('B*(0.1511) + G*(0.1973) + R*(0.3283) + NIR*(0.3407) + SWIR1*(-0.7117) + SWIR2*(-0.4559)', {
        'B': image.select(['SR_B2']),
        'G': image.select(['SR_B3']),
        'R': image.select(['SR_B4']),
        'NIR': image.select(['SR_B5']),
        'SWIR1': image.select(['SR_B6']),
        'SWIR2': image.select(['SR_B7'])
    })
    return (wet);
};
var Wet_8 = collection2.map(getWet8).median();
/*************************************************************************************/
/***************************************NDBSI_L8****************************************/
var getNDBSI8 = function(image) {
    var ibi = image.expression('(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR + RED) - GREEN / (GREEN + SWIR1))) / (2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1)))', {
        'SWIR1': image.select('SR_B6'),
        'NIR': image.select('SR_B5'),
        'RED': image.select('SR_B4'),
        'GREEN': image.select('SR_B3')
    });
    var si = image.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', {
        'SWIR1': image.select('SR_B6'),
        'NIR': image.select('SR_B5'),
        'RED': image.select('SR_B4'),
        'BLUE': image.select('SR_B2')
    });
    var ndbsi = ibi.subtract(si).divide(2);
    return (ndbsi);
}
var NDBSI_8 = collection2.map(getNDBSI8).median();
/*************************************************************************************/
/*************************************归一化-MaxMin***********************************/
//归一化函数
var img_norm_1 = function(image) {
    var minMax = image.reduceRegion({
        reducer: ee.Reducer.minMax(),
        geometry: roi,
        scale: 30,
        maxPixels: 10e13,
        //tileScale: 16
    })
    var normalize = ee.ImageCollection.fromImages(image.bandNames().map(function(name) {
        name = ee.String(name);
        var band = image.select(name);
        return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))));
    })).toBands().rename(image.bandNames());
    return normalize;
}
/*************************************************************************************/
/***********************************标准化-中心均值化*********************************/
var img_norm_2 = function(image) {
    var mean_std = image.reduceRegion({
        reducer: ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true),
        geometry: roi,
        scale: 30,
        maxPixels: 10e9,
    });
    // use unit scale to normalize the pixel values
    var unitScale = ee.ImageCollection.fromImages(image.bandNames().map(function(name) {
        name = ee.String(name);
        var band = image.select(name);
        var mean = ee.Number(mean_std.get(name.cat('_mean')));
        var std = ee.Number(mean_std.get(name.cat('_stdDev')));
        var max = mean.add(std.multiply(3))
        var min = mean.subtract(std.multiply(3))
        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))))
        var result_band = band1.subtract(min).divide(max.subtract(min));
        return result_band;
    })).toBands().rename(image.bandNames());
    return unitScale
}
/*************************************************************************************/

/*************************************归一化波段，归并**********************************/
//归一化
var norm_NDVI= img_norm_1(NDVI_8);
var norm_Wet=img_norm_1(Wet_8);
var norm_NDBSI= img_norm_1(NDBSI_8);
var norm_Lst= img_norm_1(Lst_8);

var L8_median=L8_median.addBands(norm_NDVI.rename('NDVI').toFloat());
var L8_median=L8_median.addBands(norm_Wet.rename('Wet').toFloat());
var L8_median=L8_median.addBands(norm_NDBSI.rename('NDBSI').toFloat());
var L8_median=L8_median.addBands(norm_Lst.rename('LST').toFloat());

//掩膜-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')
             .filterDate(star_date,end_date).filterBounds(roi).select('waterClass').toBands();
var visualization = {
  min: 0.0,
  max: 3.0,
  palette: ['cccccc', 'ffffff', '99d9ea', '0000ff']
};
Map.addLayer(dataset2,visualization,'water');
var mask0 = dataset2.clip(roi).eq(2);
var mask1 = dataset2.clip(roi).eq(3);
var mask_ori=mask0.add(mask1).unmask().clip(roi).eq(0);

var L8_median = L8_median.updateMask(mask_ori)
/*************************************************************************************/
var bands = ["NDVI","Wet","NDBSI","LST"];
var bandImage =L8_median.select(bands)


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

/****************************************主成分分析************************************/
//Input-使用合成波段
var region = roi;
var pre_image =  bandImage.select(bands);
var scale = 30;
var bandNames = pre_image.bandNames();

// 数据平均
// 标准化拉伸-principal components.
var meanDict = pre_image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = pre_image.subtract(means);

// 重命名
var getNewBandNames = function(prefix) {
  var seq = ee.List.sequence(1, bandNames.length());
  return seq.map(function(b) {
    return ee.String(prefix).cat(ee.Number(b).int());
  });
};

//主成分分析函数
var getPrincipalComponents = function(centered, scale, region) {
  // 图像转为一维数组
  var arrays = centered.toArray();

  // 协方差矩阵
  var covar = arrays.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
  });

  // 获取“数组”协方差结果并转换为数组。
  // 波段与波段之间的协方差
  var covarArray = ee.Array(covar.get('array'));

  // 执行特征分析
  var eigens = covarArray.eigen();

  // 分割特征值
  var eigenValues = eigens.slice(1, 0, 1);

 //计算主成分贡献度
    var eigenValuesList = eigenValues.toList().flatten()
    var total = eigenValuesList.reduce(ee.Reducer.sum())
    var percentageVariance = eigenValuesList.map(function(item) {
      return (ee.Number(item).divide(total)).multiply(100).format('%.2f')
    })
      print('特征值',eigenValuesList )
    print("贡献率（%）", percentageVariance)


  // 分割出特征向量，其特征向量为行
  var eigenVectors = eigens.slice(1, 1);
   print('eigenVectors',eigenVectors);

  // 将图像转换为二维阵列
  var arrayImage = arrays.toArray(1);

  // 使用特征向量矩阵左乘图像阵列
  var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);

  // 将特征值的平方根转换为P波段图像。
  var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);

  // 将PC转换为P波段图像，通过SD标准化。
  return principalComponents
    // 抛出一个不需要的维度，[[]]->[]。
    .arrayProject([0])
    // 使单波段阵列映像成为多波段映像，[]->image。
    .arrayFlatten([getNewBandNames('pc')])
    // 通过SDs使PC标准化
    .divide(sdImage);

};

//进行主成分分析，获得分析结果
var pcImage = getPrincipalComponents(centered, scale, region);

//基于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')});

var RSEI = img_norm_1(RSEI_0).clip(roi);

print("RSEI"+Year, ui.Chart.image.histogram(RSEI, roi, 100, 258))

var visParam = {
    palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
        '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
 };
//添加图层
 Map.addLayer(RSEI, visParam, 'RSEI_'+Year)
/********************************导出 *****************************************************/
Export.image.toDrive({
      image: RSEI,
      description: 'RSEI_'+Year,
      folder: "Annual_RSEI",
      scale:30,
      region:roi,
      fileFormat:'GeoTIFF',
      crs: 'EPSG:4326',
    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.