Google Earth Engine——基于改進(jìn)的RSEI評(píng)估生態(tài)環(huán)境（水體掩膜后）
未經(jīng)允許，禁止隨意轉(zhuǎn)載，尊重他人版權(quán)，僅供學(xué)習(xí)參考，歡迎交流。
背景介紹
遙感生態(tài)指數(shù)(Remote Sensing Ecological Index)的獲得，是使用主成分分析法耦合了綠度、濕度、干度、熱度指標(biāo)，以1-PC1（PC1，主成分分析結(jié)果的第一分量）標(biāo)準(zhǔn)化的結(jié)果作為遙感生態(tài)指數(shù)（徐涵秋，2013）。以上四個(gè)指標(biāo)中，綠度反映了植被覆蓋度，是生態(tài)環(huán)境好壞的重要指標(biāo)；濕度和干度反映了生態(tài)環(huán)境中地表與地上的水分含量多少。熱度指標(biāo)即地表溫度是區(qū)域和全球范圍地表物理過程的重要因子，也是研究地表和大氣物質(zhì)交換和能量交換的關(guān)鍵參數(shù)。因此，使用主成分分析法耦合這四個(gè)指標(biāo)，避免主觀確定權(quán)重，為評(píng)估生態(tài)環(huán)境質(zhì)量提供一種快速有效的方法。
在使用RSEI模型時(shí)，一些學(xué)者根據(jù)2013年提出的1-PC1來計(jì)算RSEI（徐涵秋，2013），但也有其他學(xué)者直接使用PC1作為區(qū)域生態(tài)遙感指數(shù)。究其原因，是學(xué)者們只是簡單地應(yīng)用模型，對(duì)模型的運(yùn)作和模型機(jī)制缺乏了解。由于主成分分析方法中特征向量方向的非唯一性，這兩個(gè)模型會(huì)帶來相反的結(jié)果。在實(shí)際研究中，學(xué)者們通常直接使用該方法，很少有人注意到特征向量及其方向，大大限制了遙感生態(tài)環(huán)境評(píng)估方法的發(fā)展。同時(shí)，由于缺乏對(duì)模型機(jī)制的研究，盲目地應(yīng)用和改進(jìn)模型，有時(shí)會(huì)誤導(dǎo)學(xué)者做出錯(cuò)誤的評(píng)價(jià)。
- 通過研究主成分分析方法中特征向量方向的改變對(duì)RSEI的影響，只有當(dāng)NDVI和Wet的特征向量為負(fù)值，NDSI和LST的特征向量為正值時(shí)，使用1-PC1計(jì)算的RSEI才是正確的。而直接用PC1的模型只有在NDVI和Wet的特征向量為正值，NDSI和LST的特征向量為負(fù)值時(shí)才能得到正確的結(jié)果（Li Ning等，2020），提出了改進(jìn)模型如式1，這里的Vndvi、Vwet是指NDVI和Wet對(duì)PC1的特征向量值。
![]()
在上述研究基礎(chǔ)上，有學(xué)者通過大樣本測試了RSEI模型特征向量在時(shí)間序列中的演變。結(jié)果發(fā)現(xiàn)，改變輸入波段的順序都可能導(dǎo)致RSEI是相反的，每個(gè)生態(tài)因子的特征向量都會(huì)影響相應(yīng)的RSEI的方向。因此提出一個(gè)改進(jìn)的模型，通過選擇受季節(jié)性變化影響較小的Wet濕度分量來判斷第一主成分PC1的特征向量方向。該模型的方向可以根據(jù)主觀經(jīng)驗(yàn)自動(dòng)修改，無需研究人員干預(yù)，改進(jìn)后的模型（如下式）可以適應(yīng)不同時(shí)期的RSEI計(jì)算，同時(shí)無論輸入指標(biāo)的順序如何變化，最終結(jié)果的方向是正確的（Zhang Yaqiu等，2021）。
除了主成分分析來進(jìn)行權(quán)重賦值，還有學(xué)者利用粒化熵方法來確定RSEI的權(quán)重。例如基于Modis數(shù)據(jù)建立了中國遙感生態(tài)指數(shù)（RSEI）模型,進(jìn)行RSEIs的知識(shí)粒化，提出了根據(jù)指標(biāo)的特征確定指標(biāo)知識(shí)粒化熵權(quán)的方法（(Liao Weihua等，2020），避免了主成分分析特征向量方向不唯一性對(duì)RSEI的影響。
這里采用第二種改進(jìn)方法，對(duì)研究區(qū)域基于改進(jìn)的RSEI評(píng)估遙感生態(tài)環(huán)境。
-
首先是對(duì)數(shù)據(jù)和方法的說明。
使用USGS Landsat 8 Level 2, Collection 2, Tier 1 (google.com)的30m分辨率數(shù)據(jù)，通過Google Earth Engine——基于新的Landsat SR數(shù)據(jù)集去云處理，對(duì)一年內(nèi)的每景影像計(jì)算四個(gè)指標(biāo)然后中值合成，進(jìn)行歸一化。在指標(biāo)合成之前，如果研究區(qū)域內(nèi)含有水體，還需要進(jìn)行水體掩膜。
-
水體掩膜
可以根據(jù)MNDWI提取水體閾值，也可以通過現(xiàn)有的全球水體數(shù)據(jù)。
一是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年數(shù)據(jù)。
//源自 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米分辨率，實(shí)際提供1984年-2019年數(shù)據(jù)。
//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代碼-指標(biāo)計(jì)算+水體掩膜+PCA
1 /*目錄*/
2 //L8去云
3 //年度數(shù)據(jù)預(yù)處理
4 //LST
5 //NDVI
6 //Wet
7 //NDBSI
8 //歸一化-MaxMin
9 //標(biāo)準(zhǔn)化-中心均值化
10 //水體掩膜
11 //歸一化波段，歸并
12 //主成分分析
13 /*************************************************************************************/
14 //導(dǎo)入自己的研究區(qū)，將其定義為roi
15 var roi = ee.FeatureCollection("users/自定義 /自定義");
16 var Year='2019';
17 var star_date = Year+'-01-01' //定義起始時(shí)間
18 var end_date = Year+'-12-31' //定義終止時(shí)間
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 // 用縮放后的波段替換，并應(yīng)用云掩膜
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 // 用縮放后的波段替換，并應(yīng)用云掩膜
35 return image.addBands(opticalBands, null, true).addBands(thermalBands, null, true).updateMask(qaMask).updateMask(saturationMask);
36 }
37 /*************************************************************************************/
38 /******************************處理一年的數(shù)據(jù)_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 //歸一化函數(shù)
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 /***********************************標(biāo)準(zhǔn)化-中心均值化*********************************/
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 // 數(shù)據(jù)平均
192 // 標(biāo)準(zhǔn)化拉伸-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 //主成分分析函數(shù)
211 var getPrincipalComponents = function(centered, scale, region) {
212 // 圖像轉(zhuǎn)為一維數(shù)組
213 var arrays = centered.toArray();
214
215 // 協(xié)方差矩陣
216 var covar = arrays.reduceRegion({
217 reducer: ee.Reducer.centeredCovariance(),
218 geometry: region,
219 scale: scale,
220 maxPixels: 1e9
221 });
222
223 // 獲取“數(shù)組”協(xié)方差結(jié)果并轉(zhuǎn)換為數(shù)組。
224 // 波段與波段之間的協(xié)方差
225 var covarArray = ee.Array(covar.get('array'));
226
227 // 執(zhí)行特征分析
228 var eigens = covarArray.eigen();
229
230 // 分割特征值
231 var eigenValues = eigens.slice(1, 0, 1);
232
233 //計(jì)算主成分貢獻(xiàn)度
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("貢獻(xiàn)率（%）", percentageVariance)
241
242
243 // 分割出特征向量，其特征向量為行
244 var eigenVectors = eigens.slice(1, 1);
245 print('eigenVectors',eigenVectors);
246
247 // 將圖像轉(zhuǎn)換為二維陣列
248 var arrayImage = arrays.toArray(1);
249
250 // 使用特征向量矩陣左乘圖像陣列
251 var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
252
253 // 將特征值的平方根轉(zhuǎn)換為P波段圖像。
254 var sdImage = ee.Image(eigenValues.sqrt())
255 .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
256
257 // 將PC轉(zhuǎn)換為P波段圖像，通過SD標(biāo)準(zhǔn)化。
258 return principalComponents
259 // 拋出一個(gè)不需要的維度，[[]]->[]。
260 .arrayProject([0])
261 // 使單波段陣列映像成為多波段映像，[]->image。
262 .arrayFlatten([getNewBandNames('pc')])
263 // 通過SDs使PC標(biāo)準(zhǔn)化
264 .divide(sdImage);
265
266 };
267
268 //進(jìn)行主成分分析，獲得分析結(jié)果
269 var pcImage = getPrincipalComponents(centered, scale, region);
270
271 //基于Wet對(duì)PC1的特征向量判斷，選擇合適的RSEI計(jì)算公式
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 /********************************導(dǎo)出 *****************************************************/
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 });
-
參考文獻(xiàn)
[1]徐涵秋. 2013. 區(qū)域生態(tài)環(huán)境變化的遙感評(píng)價(jià)指數(shù)[J]. 中國環(huán)境科學(xué),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.
浙公網(wǎng)安備 33010602011771號(hào)