一、植被覆盖度计算方法
植被覆盖度FVC(FractionalVegetationCover)定义为单位面积内绿色植被冠层垂直投影面积所占比例。FVC是衡量地表植被状况的重要指标之一,也是区域生态系统环境变化的重要指示,对水文、生态、区域变化等都具有重要意义。
像元二分模型主要对规则、均匀、单一的像元进行植被覆盖度提取。其原理是判定一个单位像元由植被与非植被两部分组成,而遥感影像中光谱波段组合由这两个组成部分权重线性合成,各因子的权重等于总像元中所占的比例,总像元中植被权重的集合与植被覆盖度具有很强的相关性。
基于遥感影像得到的全部像元信息S由地表植被覆盖信息Sveg和地表土壤覆盖信息Ssoil两部分组成,公式:
S= Sveg+ Ssoil
该像元单位中植被覆盖的像元比例大小公式:
Sveg= fc× Sveg(端元)
式中,fc表示在此遥感影像的一个像元中地表植被覆盖信息所占的比例大小。
该单位像元中非植被覆盖度信息公式:
Ssoil= (1-fc)× Ssoil(端元)
式中,1-fc表示单位像元中的非植被覆盖信息比例。
S= fc×Sveg+(1-fc)×Ssoil
即有像元二分模型中的植被覆盖度(FVC)公式:
FVC = (S- Ssoil) / (Sveg - Ssoil)
通常,像元二分法的植被覆盖度提取会基于NDVI数据计算,即NDVI代替S,公式:
FVC = (NDVI - NDVIsoil) / (NDVIseg- NDVIsoil)
通常情况下,NDVIsoil取值与NDVI累计5%值最接近,NDVIveg取值与NDVI累计95%值最接近。不同区域针对不同植被状况可适当调整置信区间。
二、FVC计算代码
// import Puer SHPvar roi = ee.FeatureCollection("users/chuanbaojiangpr/puer_merge");Map.centerObject(roi,10);// Area Displayvar style ={color:"red",fillColor:"00000000"};Map.addLayer(roi.style(style),{},"geometry");// remove cirrus cloud shadow and cloudfunction maskS2clouds(image){var qa =image.select("QA60");var cloudBitMask = 1<<10;var cirrusBitMask = 1<<11;var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));return image.updateMask(mask).divide(10000);}// S2 SR_datavar startTime = "-11-01";var endTime = "-04-30";var S2 = ee.ImageCollection("COPERNICUS/S2_SR").filterDate(startTime,endTime).filterBounds(roi).filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE",10)).map(maskS2clouds).select("B2","B3","B4","B8").median().clip(roi);// min\max:experience thresholdvar RGBvis = {min:0,max:0.5,bands:["B4","B3","B2"]};Map.addLayer(S2,RGBvis,"S2");// ndvivar NDVI = S2.expression("(NIR-R)/(NIR+R)",{"R":S2.select(["B4"]),"NIR":S2.select(["B8"])});// sort pixels(downSample,because maxPixels allows only 10000000)var threshold = NDVI.reduceRegion({reducer:ee.Reducer.percentile([5,95]),geometry:roi.geometry().bounds(),// must downsample , in order to caculate onlinescale:110// maxPixels:1e13});// debug and get the valName-B8// print(threshold)var NDVI_5 = ee.Number(threshold.get("B8_p5"));var NDVI_95 = ee.Number(threshold.get("B8_p95"));// var RGBvis = {min:-0.5,max:1};// Map.addLayer(NDVI,RGBvis,"NDVI");// FVCvar imgFVC = ((NDVI.subtract(ee.Image(NDVI_5))).divide(ee.Image(NDVI_95).subtract(ee.Image(NDVI_5)))).float();// normalize FVC 0-1var FVC = (imgFVC.lt(0).multiply(0)).add((imgFVC.gt(0).and(imgFVC.lte(1))).multiply(imgFVC)).add(imgFVC.gt(1).multiply(1));// rank and colormapvar FVC_rank = FVC.where(FVC.lt(0.2),1).where((FVC.gte(0.2).and(FVC.lt(0.4))),2).where((FVC.gte(0.4).and(FVC.lt(0.6))),3).where((FVC.gte(0.6).and(FVC.lt(0.8))),4).where(FVC.gt(0.8),5);// Map.addLayer(FVC_rank,{min:1,max:5,palette:["DCDCDC","FFEBCD","99B718","529400","011301"]},"FVC_rank")// Export FVC_rankExport.image.toDrive({image:FVC_rank,description:"FVC_dry_",folder:"users/chuanbaojiangpr/DataGET/Image",fileFormat:"GeoTIFF",region:roi,scale:110,crs:"EPSG:4326"});
三、计算过程遇到的问题与解决方法
(一)计算像素数超限
解决方法:①在统计整幅NDVI影像的NDVIsoil、NDVIveg近似值时,采用分块计算处理方法[4];②采用重采样统计处理,降低计算像素的空间分辨率,使其重采样后最大像素在1e13数量内(存在较大误差)。
(二)导出影像数据像素大小(空间分辨率)无法与Sentinel-2 MSI 可见光-近红外10m保持一致
解决方法:重采样,将影像10m空间分辨率降低至110m。
参考资料:
[1]GEE项目代码整理篇(GEE快速入门) - 知乎
[2]/earth-engine/apidocs
[3] 百度百科
[4]GEE:内存超限?将研究区划分成规则的小块运算_gee分块计算__养乐多_的博客-CSDN博客