1 背景

快速、准确获取农作物的种植面积信息对于粮食估产和种植结构调整具有重要意义，有助于相关部门及时制定农业政策和确保粮食安全[1]。冬小麦是世界主要粮食作物之一。及时了解冬小麦种植面积，对开展冬小麦长势监测和估产、区域粮食安全评估等工作具有重要的社会意义。遥感技术具有探测面积大、数据更新周期短、真实客观等特点，是快速、准确提取冬小麦种植面积的有效手段之一[2]。山东省是农业大省，冬小麦播种面积在六千万亩左右，仅次于河南省。因此对山东省冬小麦种植分布提取具有十分重要的意义。

2 区域概况

本次冬小麦种植分布的区域为山东省。山东省位于中国东部沿海、黄河下游，陆域位于东经114°48′～122°42′、北纬34°23′～38°17′，东西长721.03公里，南北长437.28公里。境域包括半岛和内陆两部分，山东半岛突出于渤海、黄海之中，同辽东半岛遥相对峙；内陆部分自北而南与河北、河南、安徽、江苏4省接壤。全省陆域面积15.58万平方公里，海洋面积15.96万平方公里。山东省境内中部山地突起，西南、西北低洼平坦，东部缓丘起伏，形成以山地丘陵为骨架、平原盆地交错环列其间的地形大势。泰山雄踞中部，主峰海拔1532.7米，为全省最高点。黄河三角洲一般海拔2～10米，为全省陆地最低处。境内地貌复杂，大体可分为平原、台地、丘陵、山地等基本地貌类型。山东省光照资源充足，光照时数年均2290-2890小时，热量条件可满足农作物一年两作的需要。年平均降水量一般在550-950毫米之间，由东南向西北递减。

enter description here

3 材料与方法

3.1 影像处理平台

影像处理使用航天宏图自主研发的PIE-Engine 云服务平台（https://engine.piesat.cn/engine/home）。

PIE-Engine是一个集实时分布式计算、交互式分析和数据可视化为一体的在线遥感云计算开放平台，主要面向遥感科研工作人员、教育工作者、工程技术人员以及相关行业用户。它基于云计算技术，汇集遥感数据资源和大规模算力资源，通过在线的按需实时计算方式，大幅降低遥感科研人员和遥感工程人员的时间成本和资源成本。用户仅需要通过基础的编程就能完成从遥感数据准备到分布式计算的全过程，这使广大遥感技术人员更加专注于遥感理论模型和应用方法的研究，在更短的时间产生更大的科研价值和工程价值。

PIE-Engine是国内遥感行业一款跨时代产品，与GEE相比，在自主可控、国产数据适配等方面具有突出优势。PIE-Engine不但提供国外的Landsat系列、Sentinel系列等卫星遥感数据，还提供了国内的高分系列、环境系列、资源系列等卫星遥感数据，同时集成了大量的遥感通用算法和专题算法。如基于多时相的Landsat和Sentinel数据，可以快速开展作物长势监测、地区旱情分析、水体变化分析、城镇变化监测等分析处理。

PIE-Engine构建了“开放+共建+共享”的遥感云生态，在云上打通了卫星应用上、中、下游，一站式提供数据获取、数据分析、产品服务。用户可在平台上传自己的数据、注册自己的算法，并与他人实现成果共享，实现资源的高效率流转与利用，激发遥感应用产业活力，推动遥感应用技术进步。如今，PIE-Engine已在自然资源、应急管理、生态环境、气象海洋等行业与领域得到广泛应用。

3.2 影像获取

使用PIE-Engine自带的覆盖山东省的Sentinel-2 L2A的10米分辨率反射率影像，时间为2020年10月至2021年7月。

Sentinel-2号是高分辨率多光谱成像卫星，携带一枚多光谱成像仪(MSI)，分为2A和2B两颗卫星，其中一颗卫星的重访周期为10天，两颗互补，重访周期为5天（纬度较高的欧洲地区，仅需3天），常用于陆地监测，可提供植被、土壤和水覆盖、内陆水路及海岸区域等图像，还可用于紧急救援服务。

Sentinel-2数据集中包含的数据为L2A级产品数据，是经过大气校正的大气底层反射率数据（Bottom-of-Atmosphere corrected reflectance）。L2A级产品数据包含12个UINT16光谱带(无卷云波段B10)，3个QA频段（其中QA60是具有云掩码信息的位掩码频段），还包含了场景分类产品、气溶胶光学厚度（AOT）、水蒸汽及部分云和降雪概率产品。

3.3 提取方法

本次使用结合物候信息的阈值法，根据冬小麦的物候特征对不同时期的影像设置阈值进行提取。

华北地区冬小麦的NDVI时序曲线见下图[3]。冬小麦出苗以后NDVI逐渐增大，进入越冬期后NDVI逐渐减小，越冬期后NDVI开始增大，到四月份达到最大，随着小麦的渐渐成熟，NDVI也相应减小。根据这些特点，可以采用下面的思路来进行小麦提取：即10月NDVI较小，12月NDVI大于10月，2月NDVI小于12月，4月NDVI大于2月，6月NDVI小于4月（此处的大于或小于指两个时期的NDVI之差大于或小于一定的阈值）。但是采用这种方法需要多个时期的影像，特别是越冬期的时候高质量影像比较少，在大范围尺度时略有困难。因此有必要进行简化。经过不断的试验，发现使用10月、4月和6月能达到较好的效果。

enter description here

实践中发现，仅用NDVI植被指数进行提取会误提大蒜等播种时间比小麦早的作物。但是，在10月份，大蒜在短波红外、近红外、红波段假彩色合成影像上呈现出墨绿色[4]，如下图所示。在RGB为SWIR1、NIR、R的假彩色合成影像上，大蒜是墨绿色。具体表现为近红外波段大，短波红外小，因此可以使用NBR植被指数（归一化燃烧指数，Normalized Burnt Ratio），即NBR=（NIR-SWIR）/（NIR+SWIR）。

在影像上选取一定的样本点，取平均值计算两种作物的NDVI和NBR植被指数。10月的小麦和大蒜两个植被指数比较见下图，可以发现NBR差异更明显。

enter description here

根据影像情况根据试验最后确定阈值为播种10月的NDVI小于0.3，NBR小于0.07，4月的NDVI大于0.32，6月的NDVI小于4月。

3.4 PIE-Engine代码

PIE-Engine的小麦提取代码如下（链接：https://engine.piesat.cn/engine-share/shareCode.html?id=f294e2fbba2c43f4ab1dde434b8d4d72）。

代码内容为：

//加载范围
var aoi = pie.FeatureCollection("user/wcgiser/haha/China_Admin2")
                           .filter(pie.Filter.eq("省","山东省"))
//添加边界  
Map.addLayer(aoi,{ color: "ffff00", fillColor: "00000000" },"boundary")
var aoi = aoi.union().first().geometry();
//加载Sentinel-2 MSI数据集并按时间、云量筛选
var cloud_free_fun = function (image) {
    var cloudMask = pie.Algorithm.Sentinel2.cloudMask(image,"Sentinel2");
    return image.select(["B11","B8","B4","B3","B2"]).updateMask(cloudMask.eq(0));
}
//播种期影像
var img1 = pie.ImageCollection('S2/L2A')
    .filterDate("2020-10-01", "2020-11-11")
    .filterBounds(aoi)
    .filter(pie.Filter.lt("cloudyPixelPercentage", 60))
    .select(["B11","B8","B4","B3", "QA60"])
    .map(cloud_free_fun)
    .median();
//拔节期影像
var img2 = pie.ImageCollection('S2/L2A')
    .filterDate("2021-03-21", "2021-05-11")
    .filterBounds(aoi)
    .filter(pie.Filter.lt("cloudyPixelPercentage", 60))
    .select(["B11","B8","B4","B3", "QA60"])
    .map(cloud_free_fun)
    .median();
//成熟收获期影像
var img3 = pie.ImageCollection('S2/L2A')
    .filterDate("2021-05-30", "2021-07-05")
    .filterBounds(aoi)
    .filter(pie.Filter.lt("cloudyPixelPercentage", 60))
    .select(["B11","B8","B4","B3","B2", "QA60"])
    .map(cloud_free_fun)
    .median();
//提取
var red1 = img1.select("B4");
var nir1 = img1.select("B8");
var swir1 = img1.select("B11");
var red2 = img2.select("B4");
var nir2 = img2.select("B8");
var red3 = img3.select("B4");
var nir3 = img3.select("B8");
var ndvi1 = (nir1.subtract(red1)).divide(nir1.add(red1)).rename("NDVI").select("NDVI");
var nbr1 = (nir1.subtract(swir1)).divide(nir1.add(swir1)).rename("NBR").select("NBR");
var ndvi2 = (nir2.subtract(red2)).divide(nir2.add(red2)).rename("NDVI").select("NDVI");
var ndvi3 = (nir3.subtract(red3)).divide(nir3.add(red3)).rename("NDVI").select("NDVI");
//大蒜在10月份的近红外波段更大，短波红外波段更小
//条件1：播种期NDVI小，NBR小
//条件2：拔节期抽穗期 NDVI大
//条件3：成熟期NDVI小于拔节期NDVI
var wheat = ndvi1.lt(0.3)
            .and(nbr1.lt(0.07))
            .and(ndvi2.gt(0.32))
            .and(ndvi3.lt(ndvi2))
            .clip(aoi)
// 结果可视化
var visParam = {
    min: 1,
    max: 1,
    palette: '000088'
};
Map.addLayer(wheat.updateMask(wheat.eq(1)),visParam, "wheat")
//加载显示影像
Map.addLayer(img1, {min: 0, max: 4000}, "img1");
Map.addLayer(img2, {min: 0, max: 6000}, "img2");
Map.addLayer(img3, {bands:["B11","B8","B4"],min: 0, max: 6000}, "img3");
//统计面积
var areaImage = pie.Image().pixelArea().multiply(wheat.eq(1));
var data = areaImage.reduceRegion(pie.Reducer.sum(), aoi, 400);
print(pie.Number(data.get("constant")).multiply(0.0015*0.0001).round(),"万亩")