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.