1. 写在前面

  今天分享一个有意思的文章，用于进行农作物分类。文章提出了一个灵活的物候辅助监督水稻(PSPR)制图框架。主要是通过提取植被物候，并自动对物候数据进行采样，获得足够多的样本点，再使用随机森林等机器学习方法进行分类。这种方法有效解决了样本量不足或者样本位置不够精确的问题，并且分类结构相较于之前的方法更高。我认为这是一种比较有意思的文章，当然这种方法还可以用到其他植被类型分类中。

灵活的物候辅助监督水稻(PSPR)制图框架：

2. 北京作物分类

   首先，使用Landsat5、7、8数据获取植被物候信息，再提取随机采样点。
  各植被类型的物候：

/**************************使用PSPR方法生成随机采样点*************************/
// PSPR: 灵活的物候辅助监督水稻制图框架
//设置研究区位置: 北京
var table = ee.FeatureCollection("users/cduthes1991/boundry/China_province_2019");
var BJGrid4 = table.filter(ee.Filter.eq('provinces','beijing'));
var roi = BJGrid4;
Map.addLayer(roi,{'color':'grey'},'roi',false);
Map.centerObject(roi,7);var GridTest = GridRegion(BJGrid4.geometry(),6,6).filterBounds(roi);
print("GridTest size:",GridTest.size());
var color = {'color':'0000FF','fillColor':'FF000000'};
Map.addLayer(GridTest.style(color),null,'GridTest');/********* cropland mask**************************************************/
//这个数据为中科院的30LUCC数据，其中11和12分别表示水田和旱地
var CAS_LULC_2018 = ee.Image("users/chengkangmk/Global-LULC-China/LULC_CAS/CAS_LULC_30m_2018");
var cropland = CAS_LULC_2018.eq(11).or(CAS_LULC_2018.eq(12)).clip(roi);
Map.addLayer(cropland.randomVisualizer(),null,'cropland',false);/*******************************自定义函数*********************************/
// Landsat 4, 5 and 7 去云
function rmL457Cloud(image) {var qa = image.select('pixel_qa');// If the cloud bit (5) is set and the cloud confidence (7) is high// or the cloud shadow bit is set (3), then it's a bad pixel.var cloud = qa.bitwiseAnd(1 << 5).and(qa.bitwiseAnd(1 << 7)).or(qa.bitwiseAnd(1 << 3));// Remove edge pixels that don't occur in all bandsvar mask2 = image.mask().reduce(ee.Reducer.min());var mask3 = image.select("B1").gt(2000);return image.updateMask(cloud.not()).updateMask(mask2).updateMask(mask3.not()).toDouble().divide(1e4).copyProperties(image).copyProperties(image, ["system:time_start",'system:time_end']);
}// Landsat-8 去云
function rmL8Cloud(image) { var cloudShadowBitMask = (1 << 3); var cloudsBitMask = (1 << 5); var qa = image.select('pixel_qa'); var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); var mask2 = image.select("B2").gt(2000);                 return image.updateMask(mask).updateMask(mask2.not()).toDouble().divide(1e4).copyProperties(image).copyProperties(image, ["system:time_start",'system:time_end']);
}// Sentinel-2 去云
function rmS2cloud(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));var mask2 = image.select("B2").lte(2000);return image.updateMask(mask).updateMask(mask2).toDouble().divide(1e4).copyProperties(image).copyProperties(image, ["system:time_start", "system:time_end"]);
}// 计算相关指数，包括NDVI、LSWI（植被水分含量指数）、EVI
function addIndex(image){// original bandsvar blue = image.select('blue'); var red = image.select('red');var green  = image.select('green');var nir = image.select('nir');var swir1 = image.select('swir1');var ndvi = image.normalizedDifference(["nir", "red"]).rename("NDVI").toDouble();var lswi = image.normalizedDifference(["nir", "swir1"]).rename("LSWI").toDouble();var lswi2ndvi = lswi.subtract(ndvi).rename("LSWI2NDVI").toDouble();var evi = image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {'NIR': nir,'RED': red,'BLUE':blue}).rename("EVI");return image.addBands(ndvi).addBands(lswi).addBands(lswi2ndvi).addBands(evi);
}// 为影像的特定波段指定名称
var LC8_BANDS = ['B2',   'B3',    'B4',  'B5',  'B6',    'B7']; //Landsat 8
var LC7_BANDS = ['B1',   'B2',    'B3',  'B4',  'B5',    'B7']; //Landsat 7
var LC5_BANDS = ['B1',   'B2',    'B3',  'B4',  'B5',    'B7']; //Llandsat 5
var S2_BANDS  = ['B2',   'B3',    'B4',  'B8',  'B11',   'B12']; // Sentinel-2
var STD_NAMES = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2'];// 定义时间变量、LSWI阈值、采样点数量
var year = '2015';
var th_lswi2NDVI = 0;
var pointNum = 1000;// 导入数据
// landsat 8
var l8Col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').map(rmL8Cloud).filterBounds(roi).filterDate(year+'-04-01',year+'-11-01')// .filter(ee.Filter.lte('CLOUD_COVER',10)).select(LC8_BANDS, STD_NAMES); 
// landsat 7 
var l7Col = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').map(rmL457Cloud).filterBounds(roi).filterDate(year+'-04-01',year+'-11-01').select(LC7_BANDS, STD_NAMES); 
// landsat 5  
var l5Col = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').map(rmL457Cloud).filterBounds(roi).filterDate(year+'-04-01',year+'-11-01').select(LC5_BANDS, STD_NAMES); var L578COl = ee.ImageCollection(l8Col.merge(l7Col).merge(l5Col)).sort("system:time_start").map(addIndex);
print("L578COl: ", L578COl);//****************************************************************************
//***************************** 基于LSWI的水稻制图****************************
//****************************************************************************
var palette1 = {min: 0, max: 1.0, palette: ['000000', '00FF00']};// 1.水淹阶段，定义植被的时间窗口
var imgCol_flood = L578COl.filterDate(year+'-05-01',year+'-06-30') .filterBounds(roi);// 使用 LSWI-NDVI 提取水稻
// "rice_0_lswi"： 表示初步产生水稻的时间
var rice_0_lswi =  imgCol_flood.select("LSWI2NDVI").map(function(image){return image.select("LSWI2NDVI").gt(th_lswi2NDVI);  
});// 从rice_0_lswi图像集合中获取每个像素位置上的LSWI指数最大值(找到最大的 LSWI 值)。
// LSWI是一种用于检测土地表面水体的指数，其数值通常与水体的存在程度成正比。在稻田中，LSWI 的最大值对应于稻田灌溉时水体充足的情况。因此这样可以最大程度的保证水稻位置的获取。
var rice_0_lswi = rice_0_lswi.reduce(ee.Reducer.max()).clip(roi).updateMask(cropland);
Map.addLayer(rice_0_lswi, palette1, "rice_0_lswi candidates",false);//2.生长阶段，水稻幼苗被移植到稻田中，直到水稻成熟
var imgCol_growth =  L578COl.filterDate(year+'-07-01',year+'-10-31').filterBounds(roi);// 不同阶段图像
var PalettePanel = {bands:["swir1","nir","red"],min:0,max:0.3};var imgCol_flood_qmosaic = imgCol_flood.qualityMosaic('LSWI2NDVI').clip(roi);
Map.addLayer(imgCol_flood_qmosaic, PalettePanel, 'imgCol_flood_qmosaic');var imgCol_growth_qmosaic = imgCol_growth.qualityMosaic('NDVI').clip(roi);
Map.addLayer(imgCol_growth_qmosaic, PalettePanel, 'imgCol_growth_qmosaic', false);// 使用 0 值来填充rice_0_lswi的区域
var rice_0_map = rice_0_lswi.unmask(0).clip(roi);
Map.addLayer(rice_0_map.selfMask(), {"palette":'#FF0000'}, "rice_0_map",false);//****************************************************************************
//*********************基于CCVS 方法提取水稻物候******************************
//****************************************************************************
var visParam = {min: -0.2,max: 0.8,palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};// 使用LSWI_max, LSWI_min, NDVI_max, NDVI_min来推导RCLN，其中RCLN表示LSWI相对于与NDVI的变化幅度
var RCLE = imgCol_growth_qmosaic.select("LSWI").subtract(imgCol_flood_qmosaic.select('LSWI')).abs().divide(imgCol_growth_qmosaic.select("EVI").subtract(imgCol_flood_qmosaic.select('EVI')).abs()).rename("RCLE");
Map.addLayer(RCLE, visParam,'RCLE', false);// 创建一个布尔类型的图像，其中大于 0.1 的像素被设置为 true，小于或等于 0.1 的像素被设置为 false
var LSIW_min = imgCol_flood.select('LSWI').min().gt(0.1).clip(roi).updateMask(cropland);
Map.addLayer(LSIW_min.selfMask(), {"palette":'#FF0000'}, "LSIW_min",false);// 
var RCLE_rice = RCLE.updateMask(RCLE.gt(0)).updateMask(LSIW_min).lt(0.6).unmask(0).clip(roi); 
Map.addLayer(RCLE_rice.selfMask(), {palette: 'green'}, "CCVS_rice",false);var riceCombine = RCLE_rice.add(rice_0_map);// 将以上两种方法得到的数据图像进行合并求交
var stableMask = riceCombine.eq(0).or(riceCombine.eq(2));
var LULU_mutual = riceCombine.updateMask(stableMask).where(riceCombine.eq(2),1).rename("RiceMap").updateMask(cropland);
Map.addLayer(LULU_mutual,palette1,'LULU_mutual',false);var riceMapCol = GridTest.toList(GridTest.size()).map(function(ROIFea){// set study areavar roi = ee.FeatureCollection([ROIFea]);var ricemap = riceTrainData(roi);return ee.FeatureCollection(ricemap);
});
var samplePoint = ee.FeatureCollection(riceMapCol).flatten();
Map.addLayer(samplePoint,null,'samplePoint',false);// check the generated sample data
print("samplePoint:",samplePoint.limit(10));var ricePoint_1 = samplePoint.filter(ee.Filter.eq('landcover',1));
Map.addLayer(ricePoint_1,{'color':'#FFA500'},'ricePoint_1');
print("ricePoint_1 size",ricePoint_1.size());var NonricePoint_1 = samplePoint.filter(ee.Filter.eq('landcover',0));
print("NonricePoint_1 size",NonricePoint_1.size());Export.table.toAsset(samplePoint,"PSPRSampleGeneration"+year,"PSPRSampleGeneration"+year)function riceTrainData(roiRegion){var roi = roiRegion;var LULU_mutual2 = LULU_mutual.clip(roi);/***********************************************************************Define neighboor function* and generate the samples***********************************************************************/function neighFun(img,kernalRadius,roi){var kernel = ee.Kernel.square(kernalRadius,'pixels',false);var kernelArea = (ee.Number(kernalRadius).multiply(2).add(1)).pow(2);var imgNeibor = ee.Image(img).convolve(kernel).eq(kernelArea).set("system:footprint",roi.geometry());return img.updateMask(imgNeibor);}var samplePoint = ee.List([]);for(var i=0;i<=1;i++){var class_Num = i;var class_i_mask = neighFun(LULU_mutual2.eq(class_Num),1,roi);var class_i = LULU_mutual2.updateMask(class_i_mask);samplePoint = samplePoint.add(class_i);}var samplePoint = ee.ImageCollection(samplePoint).mosaic().rename("landcover").updateMask(cropland);var pointSample =  samplePoint.stratifiedSample({numPoints:pointNum,classBand:"landcover",region:roi.geometry(),scale:30,seed:0,// tileScale:8,geometries:true});return pointSample;
}/**************************************************************************
generate the grid
***************************************************************************/
function generateGrid(xmin, ymin, xmax, ymax, dx, dy) {var xx = ee.List.sequence(xmin, ee.Number(xmax).subtract(0.0001), dx);var yy = ee.List.sequence(ymin, ee.Number(ymax).subtract(0.0001), dy);var cells = xx.map(function(x) {return yy.map(function(y) {var x1 = ee.Number(x);var x2 = ee.Number(x).add(ee.Number(dx));var y1 = ee.Number(y);var y2 = ee.Number(y).add(ee.Number(dy));var coords = ee.List([x1, y1, x2, y2]);var rect = ee.Algorithms.GeometryConstructors.Rectangle(coords); return ee.Feature(rect);});}).flatten(); return ee.FeatureCollection(cells);
}function GridRegion(roiRegion,xBlock,yBlock){//roiRegion: area of interest in the form of geometry// compute the coordinatesvar bounds = roiRegion.bounds();var coords = ee.List(bounds.coordinates().get(0));var xmin = ee.List(coords.get(0)).get(0);var ymin = ee.List(coords.get(0)).get(1);var xmax = ee.List(coords.get(2)).get(0);var ymax = ee.List(coords.get(2)).get(1);var dx = (ee.Number(xmax).subtract(xmin)).divide(xBlock); //4var dy = (ee.Number(ymax).subtract(ymin)).divide(yBlock);var grid = generateGrid(xmin, ymin, xmax, ymax, dx, dy);  grid = grid.filterBounds(roiRegion); return grid;
}

植被水分含量指数（LSWI） 在其公式中使用了NIR和SWIR通道，SWIR波段对植被含水量的变化较敏感：

结果展示：