未经允许，禁止随意转载，尊重他人版权，仅供学习参考，欢迎交流。
背景介绍
遥感生态指数(Remote Sensing Ecological Index)的获得，是使用主成分分析法耦合了绿度、湿度、干度、热度指标，以1-PC1（PC1，主成分分析结果的第一分量）标准化的结果作为遥感生态指数（徐涵秋，2013）。以上四个指标中，绿度反映了植被覆盖度，是生态环境好坏的重要指标；湿度和干度反映了生态环境中地表与地上的水分含量多少。热度指标即地表温度是区域和全球范围地表物理过程的重要因子，也是研究地表和大气物质交换和能量交换的关键参数。因此，使用主成分分析法耦合这四个指标，避免主观确定权重，为评估生态环境质量提供一种快速有效的方法。
在使用RSEI模型时，一些学者根据2013年提出的1-PC1来计算RSEI（徐涵秋，2013），但也有其他学者直接使用PC1作为区域生态遥感指数。究其原因，是学者们只是简单地应用模型，对模型的运作和模型机制缺乏了解。由于主成分分析方法中特征向量方向的非唯一性，这两个模型会带来相反的结果。在实际研究中，学者们通常直接使用该方法，很少有人注意到特征向量及其方向，大大限制了遥感生态环境评估方法的发展。同时，由于缺乏对模型机制的研究，盲目地应用和改进模型，有时会误导学者做出错误的评价。
- 通过研究主成分分析方法中特征向量方向的改变对RSEI的影响，只有当NDVI和Wet的特征向量为负值，NDSI和LST的特征向量为正值时，使用1-PC1计算的RSEI才是正确的。而直接用PC1的模型只有在NDVI和Wet的特征向量为正值，NDSI和LST的特征向量为负值时才能得到正确的结果（Li Ning等，2020），提出了改进模型如式1，这里的Vndvi、Vwet是指NDVI和Wet对PC1的特征向量值。
-
在上述研究基础上，有学者通过大样本测试了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
1 /*目录*/
2 //L8去云
3 //年度数据预处理
4 //LST
5 //NDVI
6 //Wet
7 //NDBSI
8 //归一化-MaxMin
9 //标准化-中心均值化
10 //水体掩膜
11 //归一化波段，归并
12 //主成分分析
13 /*************************************************************************************/
14 //导入自己的研究区，将其定义为roi
15 var roi = ee.FeatureCollection("users/自定义 /自定义");
16 var Year='2019';
17 var star_date = Year+'-01-01' //定义起始时间
18 var end_date = Year+'-12-31' //定义终止时间
19 var L8_SR = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") //加载L8_SR影像
20 /*************************************************************************************/
21 /****************************************L8 去云****************************************/
22 // 使用Landsat8 Collection 2,Level 2 QA_PIXEL波段（CFMask）去云
23 function maskL8sr(image) {
24 // Bit 0 - Fill
25 // Bit 1 - Dilated Cloud
26 // Bit 2 - Cirrus
27 // Bit 3 - Cloud
28 // Bit 4 - Cloud Shadow
29 var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
30 var saturationMask = image.select('QA_RADSAT').eq(0);
31 // 用缩放后的波段替换，并应用云掩膜
32 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 // 用缩放后的波段替换，并应用云掩膜
35 return image.addBands(opticalBands, null, true).addBands(thermalBands, null, true).updateMask(qaMask).updateMask(saturationMask);
36 }
37 /*************************************************************************************/
38 /******************************处理一年的数据_L8***************************************/
39 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
48 /*************************************************************************************/
49 /***************************************LST_L8****************************************/
50 //LST直接就是SR的红外波段，单位是开尔文
51 var Lst_8 = L8_median.select('ST_B10');
52 /*************************************************************************************/
53 /***************************************NDVI_L8****************************************/
54 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
60 /*************************************************************************************/
61 /***************************************Wet_L8****************************************/
62 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 /*************************************************************************************/
75 /***************************************NDBSI_L8****************************************/
76 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 /*************************************************************************************/
94 /*************************************归一化-MaxMin***********************************/
95 //归一化函数
96 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
103 })
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 /*************************************************************************************/
112 /***********************************标准化-中心均值化*********************************/
113 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
121 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 /*************************************************************************************/
135
136 /*************************************归一化波段，归并**********************************/
137 //归一化
138 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
143 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
148 //掩膜-MNDWI
149 // var getMNDWI8 = function(image){
150 // var mndwi = image.normalizedDifference(['SR_B3', 'SR_B6']);
151 // return (mndwi);
152 // }
153 // var MNDWI_8 = collection2.map(getMNDWI8).median().clip(roi);
154 // //Map.addLayer(MNDWI_8);
155 // //Treshhold-0.2
156 // //print("MNDWI_8", ui.Chart.image.histogram(MNDWI_8, roi, 100, 258))
157 // var mask = MNDWI_8.lt(0.2);
158 /*************************************水体掩膜 **********************************/
159 //JRC年度水分类历史-30m
160 //0 cccccc No data
161 //1 ffffff Not water
162 //2 99d9ea Seasonal water
163 //3 0000ff Permanent water
164 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
176 var L8_median = L8_median.updateMask(mask_ori)
177 /*************************************************************************************/
178 var bands = ["NDVI","Wet","NDBSI","LST"];
179 var bandImage =L8_median.select(bands)
180
181
182 /*************************************************************************************/
183
184 /****************************************主成分分析************************************/
185 //Input-使用合成波段
186 var region = roi;
187 var pre_image = bandImage.select(bands);
188 var scale = 30;
189 var bandNames = pre_image.bandNames();
190
191 // 数据平均
192 // 标准化拉伸-principal components.
193 var meanDict = pre_image.reduceRegion({
194 reducer: ee.Reducer.mean(),
195 geometry: region,
196 scale: scale,
197 maxPixels: 1e9
198 });
199 var means = ee.Image.constant(meanDict.values(bandNames));
200 var centered = pre_image.subtract(means);
201
202 // 重命名
203 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
210 //主成分分析函数
211 var getPrincipalComponents = function(centered, scale, region) {
212 // 图像转为一维数组
213 var arrays = centered.toArray();
214
215 // 协方差矩阵
216 var covar = arrays.reduceRegion({
217 reducer: ee.Reducer.centeredCovariance(),
218 geometry: region,
219 scale: scale,
220 maxPixels: 1e9
221 });
222
223 // 获取“数组”协方差结果并转换为数组。
224 // 波段与波段之间的协方差
225 var covarArray = ee.Array(covar.get('array'));
226
227 // 执行特征分析
228 var eigens = covarArray.eigen();
229
230 // 分割特征值
231 var eigenValues = eigens.slice(1, 0, 1);
232
233 //计算主成分贡献度
234 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
242
243 // 分割出特征向量，其特征向量为行
244 var eigenVectors = eigens.slice(1, 1);
245 print('eigenVectors',eigenVectors);
246
247 // 将图像转换为二维阵列
248 var arrayImage = arrays.toArray(1);
249
250 // 使用特征向量矩阵左乘图像阵列
251 var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
252
253 // 将特征值的平方根转换为P波段图像。
254 var sdImage = ee.Image(eigenValues.sqrt())
255 .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
256
257 // 将PC转换为P波段图像，通过SD标准化。
258 return principalComponents
259 // 抛出一个不需要的维度，[[]]->[]。
260 .arrayProject([0])
261 // 使单波段阵列映像成为多波段映像，[]->image。
262 .arrayFlatten([getNewBandNames('pc')])
263 // 通过SDs使PC标准化
264 .divide(sdImage);
265
266 };
267
268 //进行主成分分析，获得分析结果
269 var pcImage = getPrincipalComponents(centered, scale, region);
270
271 //基于Wet对PC1的特征向量判断，选择合适的RSEI计算公式
272 //Vwet<0
273 //var RSEI_0 = L8_median.expression('1-b1',{b1:pcImage.select('pc1')});
274
275 //Vwet>=0
276 var RSEI_0 = L8_median.expression('0+b1',{b1:pcImage.select('pc1')});
277
278 var RSEI = img_norm_1(RSEI_0).clip(roi);
279
280 print("RSEI"+Year, ui.Chart.image.histogram(RSEI, roi, 100, 258))
281
282 var visParam = {
283 palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
284 '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
285 };
286 //添加图层
287 Map.addLayer(RSEI, visParam, 'RSEI_'+Year)
288 /********************************导出 *****************************************************/
289 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}
298 });
-
参考文献
[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.