用于地市及县两级粮食补贴核算的种植面积遥感估算方法

摘要

本发明公开用于地市及县两级粮食补贴核算的种植面积遥感估算方法，包括如下步骤：步骤P1：对省域内粮食作物的种植面积抽样调查；步骤P2：遥感影像选择与预处理；步骤P3：粮食作物遥感识别；步骤P4：混合像元的处理；步骤P5：建立粮食作物种植面积估算模型，对区县级行政单元的大宗粮食作物种植面积进行估算。

说明书

技术领域

本发明涉及遥感监测方法，特别涉及一种用于地市及县两级粮食补贴核 算的粮食种植面积遥感估算方法。

背景技术

粮食直接补贴政策是近年来我国最为重要的政府支农政策之一，旨在 “进一步加大国家惠农、支农政策力度，更好地调动农民种粮积极性，保护 种粮农民利益”。随着对农民补贴资金的不断增加，政府对补贴资金的管理 也需要加强。目前我国各省均启动了以粮食种植面积基础的补贴政策。为了 实施粮食种植的补贴，需要对地市及县级行政单元的粮食种植面积进行核算。 目前主要是采用统计上报方式，由各县上报粮食种植面积，然后由省级财政 部门根据上报的粮食种植面积核算补贴金额。

通过统计方式估算农作物种植面积是目前统计部门的主要手段，但是由 于样本布设、地面调查和统计误差，以及主观方面因素，数据质量难以保证。 而且由于该数据是用于计算各县自己的粮食补贴金额，因此存在明显的虚报 现象，使国家财政资金存在较大的浪费现象。因此必须发展一种客观、可靠 的监测方法，而遥感技术因其在时间与空间尺度上的优越性，在快速准确监 测大宗粮食作物种植面积方面具备优势。然而由于粮食补贴一般在区县级行 政单元上进行，因而要求具有在区县级尺度具有较高的准确性（优于85%）， 使利用遥感方法监测大区域内精确到区县级的大宗粮食作物种植面积成为一 项难度较大又急需解决的任务。

发明内容

针对现有技术中存在的不足，本发明的目的在于提供一种用于地市及县 两级粮食补贴核算的种植面积遥感估算方法，能够有效降低在大范围内靠地 面调查获取作物种植面积信息的成本，同时也能够有效地避免地面调查过程 中因人工参与而产生的主观影响，为我国粮食补贴政策的执行提供客观、准 确的基础信息。

本发明的技术方案是这样实现的：一种用于地市及县两级粮食补贴核算 的种植面积遥感估算方法，包括如下步骤：

步骤P1：对目标省域内粮食作物的种植面积进行抽样调查；

步骤P2：遥感影像选择与预处理；

步骤P3：粮食作物遥感识别；

步骤P4：混合像元的处理；

步骤P5：建立粮食作物的种植面积估算模型，对地市及县两级行政单元 的粮食作物的种植面积进行估算。

上述用于地市及县两级粮食补贴核算的种植面积遥感估算方法，在步骤 P1中：以研究区历史农作物种植面积和种植比例结构为基础，建立省域范围 的面积抽样框架，并根据粮食作物的种植结构、期望精度与可靠性计算样本 量，并采用系统布样方式进行布样，具体包括如下步骤：

步骤P101：利用地理信息系统工具建立覆盖省域范围的矩形格网，矩形 格网尺寸大小为：长度和宽度均为5公里，删除跨越省界的矩形格网单元， 按照自上而下、自左而右的顺序对每个矩形格网单元进行编号，使得每个矩 形格网单元具有唯一的编号；

步骤P102：利用研究区最近5年的粮食作物的种植面积和种植比例结构， 计算样本量，样本量计算公式如下：

$n_0={\left(\frac{t}{r}\right)}^2\frac{S^2}{{\bar{Y}}^2},$

n=n₀+n₀*5%；

式中：n₀是初始样本量，n是修正样本量；t是抽样概率，当置信水平为 95%时，t等于1.96；r是相对误差，研究中取r=5%；是总体均值，即研究区 最近5年的粮食作物的种植面积平均；S²是总体方差，即研究区最近5年的 粮食作物的种植面积方差；

步骤P103：采用系统布样方法，选择地面调查样方，具体包括如下步骤： （1）用总的矩形格网数除以样本量，得到布样间隔；（2）随机抽取一个介于 0和布样间隔数值之间的随机数，作为布样的起始值，即抽取第一个样本；（3） 根据上一个样本的矩形格网编号，累加一个布样间隔，得到一个新的样本的 矩形格网单元编号；（4）重复步骤（3），直至结束；

步骤P104：人工调查每个样本矩形格网内的每块耕地上的农作物种植类 别；

步骤P105：统计推断每种作物的种植成数，并根据全省耕地面积统计估 算全省的种植面积，即抽样方法估算的全省各农作物种植面积；推断公式如 下：

农作物总种植面积=农作物种植成数*总耕地面积。

上述用于地市及县两级粮食补贴核算的种植面积遥感估算方法，在步骤 P2中：根据作物物候期，确定合适时相的环境星CCD数据作为数据源，并对 选定的环境星CCD数据进行预处理，具体包括如下步骤：

步骤P201：参考研究区历史物候数据，分别在作物播种期后30天、拔 节期、开花期和收获前期获取4景影像；

步骤P202：以1:10万地形图为参考影像，应用二次多项式方法对环境 星CCD影像进行几何精校正，误差控制在1个像元以内；

步骤P203：对环境星CCD影像进行辐射定标处理，定标采用如下公式进 行：

$\rho =(\frac{\mathrm{DN}}{A}+L_0)*\frac{\pi *d^2}{\mathrm{ESUN}*\cos \theta },$

公式中，DN为遥感影像的灰度值，A和L₀为中国资源卫星中心公布的环境 减灾星座A/B星各载荷在轨绝对辐射定标系数，d为日地天文单位距离，ESUN 为大气顶层太阳辐照度，θ为太阳天顶角；

步骤P204：计算光学遥感影像的比值植被指数RVI、归一化植被指数 NDVI、增强型植被指数EVI：

$\mathrm{RVI}=\frac{{\rho }_{\mathrm{Nir}}}{{\rho }_{\mathrm{Red}}};$

$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{Nir}}-{\rho }_{\mathrm{Red}}}{{\rho }_{\mathrm{Nir}}+{\rho }_{\mathrm{Red}}};$

$\mathrm{EVI}=(1.5+L)\times \frac{{\rho }_{\mathrm{Nir}}-{\rho }_{\mathrm{Red}}}{L+{\rho }_{\mathrm{Nir}}+C_1\times {\rho }_{\mathrm{Red}}+C_2\times {\rho }_{\mathrm{Blue}}}$

其中，ρ_Blue为环境星CCD数据蓝波段反射率，ρ_Nir为近红外反射率，ρ_Red为 红波段反射率。L是调整参数，C₁、C₂为大气修正参数，它们均可以减少背景 和大气的作用，对于环境星影像来说，L、C₁、C₂的经验值分别为1、6和7.5；

步骤P205：将NDVI数据与环境星CCD数据的4个波段叠加在一起，生 成具有5个波段的合成图像；

步骤P206：耕地范围提取：利用全省耕地数据库提取出全省的耕地分布 数据，以全省耕地分布数据为掩膜，对合成的环境星图像做掩膜处理，切割 出耕地，作为后续处理的基础数据。

上述用于地市及县两级粮食补贴核算的种植面积遥感估算方法，在步骤 P3中：利用最大似然分类或支持向量机分类方法，以步骤P104所述的地面 样方调查结果为训练样本，对环境星CCD影像进行分类，识别不同的粮食作 物，具体包括如下步骤：

步骤P301：调入地面调查样本矢量数据，利用地理信息系统工具将其分 成两个面积相等的部分，即将样本分成两个样本库；分别将不同的粮食作物 多边形提取出来，并对每个多边形进行编号，使得每个多边形具有唯一标识 号；按照多边形的唯一标识号，将粮食作物的多边形分为2个部分，奇数记 录第1部分，偶数记录为第2部分，同时删除面积小于8100平方米的多边形；

步骤P302：选取步骤P301所得的两个样本库中的任何一部分，作为农 作物遥感分类的训练样本库，从中选取地块较大且色彩稳定的多边形，作为 光谱训练集，然后利用光谱训练工具，训练出不同农作物的光谱特征；

步骤P303：采用最大似然法或支持向量机进行影像分类，将每个像元赋 为某种农作物，最终得到影像的农作物分类结果，并对分类结果进行重编码， 把一种作物对应的像元值设为1，其他作物设为0，得到该作物纯像元分布数 据集；对每种作物重复此步骤，分别得到每种作物纯像元分布数据集，从而 得到不同农作物的分布图。

上述用于地市及县两级粮食补贴核算的种植面积遥感估算方法，在步骤 P4中：在步骤P3分类的基础上，对于多种农作物交错区域的混合像元，建 立基于光谱特征的分解模型，计算不同像元内各种作物的面积百分比，具体 包括如下步骤：

步骤P401：对P206所获影像进行非监督分类，获取不同农作物的种植 区域，并去除步骤P303所获得的纯象元区域，得到混合像元分布区；

步骤P402：从步骤P3的光谱训练特征中读取不同作物的光谱特征（平 均值），计算每种农作物和裸地RVI、NDVI、EVI的均值，建立线性混合像元 分解模型，计算公式为：

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i+p_l=1$

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i*{\mathrm{RVI}}_i+p_l*{\mathrm{RVI}}_l={\mathrm{RVI}}_p$

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i*{\mathrm{NDVI}}_i+p_l*{\mathrm{NDVI}}_l={\mathrm{NDVI}}_p$

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i*{\mathrm{EVI}}_i+p_l*{\mathrm{EVI}}_l={\mathrm{EVI}}_p$

其中，n为农作物种植数量，p_i为一种农作物所占的像元比例，p_l为该 裸地所占的像元比例，RVI_i、NDVI_i、EVI_i为该农作物的训练样本的平均RVI 值、NDVI值、EVI值，RVI_l、NDVI_l、EVI_l裸地的训练样本的平均RVI值、NDVI 值、EVI值，RVI_p、NDVI_p、EVI_p为对应像元的RVI值、NDVI值、EVI值；

步骤P403：采用线性光谱分解模型，求解不同像元的不同作物种植面积 比例。

上述用于地市及县两级粮食补贴核算的种植面积遥感估算方法，在步骤 P5中：以步骤P1、步骤P3和步骤P4为基础，建立粮食作物种植面积估算模 型，对区县级行政单元的粮食作物种植面积进行估算，具体包括如下步骤：

步骤P501：根据不同像元的面积比例，统计不同地市及县两级行政单元 内不同粮食作物的像元数据，并根据像元大小计算不同地市及县两级行政单 元内各种作物的毛种植面积；

步骤P502：利用不同地市及县两级行政单元内的耕地系数（耕地占图斑 面积的比例），估算每个地市及县两级行政单元内各种作物的净种植面积：

净种植面积=毛种植面积*耕地系数；

步骤P503：用回归估算模型，抽样得到的作物的种植成数与分类得到的 作物的种植成数相互修正，计算各地市及县两级行政单元各种粮食作物的种 植面积；

步骤P504：统计汇总得到全省各种粮食作物的种植面积。

本发明具有以下有益效果：

1.充分利用遥感技术在大范围作物识别方面的能力，从而有效降低了在 大范围内靠地面调查获取作物种植面积信息的成本，同时也有效避免了地面 调查过程中因人工参与而产生的主观影响；

2.以我国自主研发的环境1号卫星数据为数据源，结合地面样方调查数 据，通过CCD数据各波段反射率和归一化植被指数，分析大宗粮食作物的光 谱、纹理等特征，并利用最大似然分类、支持向量机等方法识别大宗粮食作 物，在此基础上，结合国土资源第二次调查耕地数据库和田埂系数数据库， 建立大宗粮食作物种植面积估算模型，并对模型进行标定，发展出的基于遥 感的作物种植面积监测方法，将报告单元由省级行政单元细化到县级行政单 元，填补了在区县级报告单元的大宗粮食作物种植面积估算领域的空白，为 粮食补贴政策的基础性种植面积估算提供重要的技术支撑，从而为我国粮食 补贴政策的执行提供了客观、准确的基础信息；

3.环境1号卫星具有较高的空间和时间分辨率，在空间分辨率上明显优 于MODIS数据，在时间分辨率上明显优于Landsat TM数据，同时其幅宽较宽， 能够保证作物种植面积业务化监测的需要。

4.适用于小麦、玉米、大豆、水稻、马铃薯等粮食作物的种植面积估算。

附图说明

图1早稻样方调查面积与环境星分类面积回归分析；

图2本发明用于粮食补贴的作物种植面积遥感估算方法流程图；

图3中稻样方调查面积与环境星分类面积回归分析；

图4晚稻样方调查面积与环境星分类面积回归分析。

具体实施方式

结合附图对本发明做进一步的说明：

如图2所示，本发明用于地市及县两级粮食补贴核算的种植面积遥感估 算方法包括如下步骤：

步骤P1：对目标省域内粮食作物的种植面积进行抽样调查。

以研究区历史农作物种植面积和种植比例结构为基础，建立省域范围的 面积抽样框架，并根据粮食作物的种植结构、期望精度与可靠性计算样本量， 并采用系统布样方式进行布样，具体包括如下步骤：

步骤P101：利用ArcGIS等地理信息系统（GIS）工具产生矩形格网，矩 形格网尺寸大小为：长度和宽度均为5公里，要求覆盖省域范围。删掉横跨 省界的矩形格网单元，然后对每个矩形格网单元按照自上而下，自左而右的 顺序进行编号，使得每个矩形格网单元具有唯一的编号。

步骤P102：利用研究区最近5年的粮食作物的种植面积和种植比例结构 数据，计算样本量，样本量计算公式如下：

$n_0={\left(\frac{t}{r}\right)}^2\frac{S^2}{{\bar{Y}}^2},$

n=n₀+n₀*5%；

式中：n₀是初始样本量，n是修正样本量；t是抽样概率，当置信水平为 95%时，t等于1.96；r是相对误差，研究中取5%；是总体均值，即5年的 粮食作物的种植面积平均；S²是总体方差，即研究区最近5年的粮食作物的 种植面积方差。为了保证抽样精度，样本量在实际执行过程中上调5%。

步骤P103：采用系统布样方法，选择地面调查样方，具体操作是：（1） 用总的矩形格网数除以样本量，得到布样间隔；（2）随机抽取一个介于0和 布样间隔数值之间的随机数，作为布样的起始值，即抽取第一个样本；（3） 根据上一个样本的矩形格网编号，累加一个布样间隔，得到一个新的样本的 矩形格网单元编号；（4）重复步骤（3），直至结束。

步骤P104：人工调查每个样本矩形格网内的每块耕地上的作物种植类 别。具体操作是：利用地理信息系统（GIS）处理样本格网内的图斑图，用差 分GPS测量勾绘出每个地块的边界，然后根据地块内的作物类别填写相应的 编号，具体编码如下：

1-冬小麦；2-一季稻/中稻；3-玉米；4-春小麦；5-早稻；6-晚稻；7- 大豆；9-高粱；10-马铃薯；11-其它豆类；12-休耕；13-弃耕。

步骤P105：统计推断每种作物的种植成数，并根据全省耕地面积统计估 算全省的种植面积，即抽样方法估算的全省各农作物种植面积；推断公式如 下：

农作物总种植面积=农作物种植成数*总耕地面积。

步骤P2：遥感影像选择与预处理。

根据作物物候期，确定合适时相的环境星CCD数据作为数据源，并对选 定的环境星CCD数据进行预处理，具体包括：

步骤P201：参考研究区历史物候数据，分别在作物播种期后30天、拔 节期、开花期和收获前期获取4景影像；

根据目标估算区（监测区或省份）的作物物候期选择合适的环境星CCD 影像，自播种后1个月至收割前的影像均可使用，为了提高分类精度，要求 至少2期影像，相隔时间应超过20天。

步骤P202：以1:10万地形图为参考影像，对环境星CCD影像进行几何 精校正，校正方法采用二次多项式方法，重采样方法采用最邻近距离法；误 差控制：平原地区要求误差不超过1个像元，山区误差不超过1.5个像元。

步骤P203：对环境星CCD影像进行辐射定标处理，定标采用如下公式进 行：

$\rho =(\frac{\mathrm{DN}}{A}+L_0)*\frac{\pi *d^2}{\mathrm{ESUN}*\cos \theta },$

公式中，DN为遥感影像的灰度值，A和L₀为中国资源卫星中心公布的环境 减灾星座A/B星各载荷在轨绝对辐射定标系数，d为日地天文单位距离，ESUN 为大气顶层太阳辐照度，θ为太阳天顶角；

步骤P204：计算比值植被指数RVI、归一化植被指数NDVI、增强型植被 指数EVI，计算公式如下：

$\mathrm{RVI}=\frac{{\rho }_{\mathrm{Nir}}}{{\rho }_{\mathrm{Red}}};$

$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{Nir}}-{\rho }_{\mathrm{Red}}}{{\rho }_{\mathrm{Nir}}+{\rho }_{\mathrm{Red}}};$

$\mathrm{EVI}=(1.5+L)\times \frac{{\rho }_{\mathrm{Nir}}-{\rho }_{\mathrm{Red}}}{L+{\rho }_{\mathrm{Nir}}+C_1\times {\rho }_{\mathrm{Red}}+C_2\times {\rho }_{\mathrm{Blue}}}$

其中，ρ_Blue为环境星CCD数据蓝波段反射率，ρ_Nir为近红外反射率，ρ_Red为 红波段反射率。L是调整参数，C₁、C₂为大气修正参数，它们均可以减少背景 和大气的作用，对于环境星影像来说，L、C₁、C₂的经验值分别为1、6和7.5；

步骤P205：将归一化植被指数（NDVI）数据与环境星CCD数据的4个波 段叠加在一起，生成具有5个波段的合成图像。

步骤P206：耕地范围提取：利用目标估算区（如全省）耕地数据库提取 出全省的耕地分布数据，以目标估算区（如全省）耕地分布数据为掩膜，对 合成的环境星图像做掩膜处理，切割出耕地，作为后续处理的基础数据。

步骤P3：粮食作物遥感识别。

利用最大似然分类或支持向量机分类方法，以地面样方调查为训练样本， 对环境星CCD影像进行分类，识别不同的粮食作物，具体包括如下步骤：

步骤P301：调入地面调查样本矢量数据，利用地理信息系统工具将其分 成两个面积相等的部分（即分成两个样本库），分别将不同粮食作物的多边形 提取出来，并对每个多边形进行编号，使得每个多边形具有唯一标识号。

按照多边形的唯一标识号，将粮食作物的多边形分为2个部分（如：奇 数记录为第1部分，偶数记录为第2部分）；同时删除面积小于8100平方米 （9个环境星像元）的多边形【面积太小而没有代表性，没法使用】。

步骤P302：选取步骤P301所得的两个样本库中的任一部分，作为农作 物遥感分类的训练样本库，从中选取地块较大且对应的环境星数据色彩稳定 的多边形，作为光谱训练集，然后利用光谱训练工具，计算多边形对应像元 的各波段的最大值、最小值、均值、方差等，作为不同农作物的光谱特征。

步骤P303：采用最大似然法或支持向量机进行影像分类，将每个像赋为 某种农作物。

对分类结果进行重编码，把一种作物对应的像元值设为1，其他作物设 为0，得到该作物纯像元分布。对每种作物重复此步骤，分别得到每种作物 纯像元分布，从而得到不同农作物的分布图。

步骤P4：混合像元的处理。

对多种农作物交错区域进行混合像元分解，建立基于光谱特征的线性分 解模型，计算不同像元内各种作物的面积百分比，具体包括如下步骤：

步骤P401：对合成的环境星数据进行非监督分类，分类后处理，把像元 集群组分别划分为各种地类，以相应的代码表示。非监督分类方法获取农作 物种植区域，去除监督分类方法获取的纯象元区域，得到混合像元分布区。

步骤P402：根据步骤P3中的光谱特征，计算每种农作物和裸地RVI、NDVI、 EVI的均值，利用线性混合像元分解模型，建立方程，计算像元中该农作物 所占的比例，计算公式为：

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i+p_l=1$

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i*{\mathrm{RVI}}_i+p_l*{\mathrm{RVI}}_l={\mathrm{RVI}}_p$

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i*{\mathrm{NDVI}}_i+p_l*{\mathrm{NDVI}}_l={\mathrm{NDVI}}_p$

$\underset{i=1}{\overset{n}{\Sigma }}{p}_i*{\mathrm{EVI}}_i+p_l*{\mathrm{EVI}}_l={\mathrm{EVI}}_p$

其中，n为农作物种植数量，p_i为一种农作物所占的像元比例，p_l为该 裸地所占的像元比例，RVI_i、NDVI_i、EVI_i为该农作物的训练样本的平均RVI 值、NDVI值、EVI值，RVI_l、NDVI_l、EVI_l裸地的训练样本的平均RVI值、NDVI 值、EVI值，RVI_p、NDVI_p、EVI_p为对应像元的RVI值、NDVI值、EVI值；

步骤P403：采用线性光谱分解模型，求解不同像元的不同作物种植面积 比例。

步骤P5：以步骤P1、步骤P3和步骤P4为基础，建立粮食作物种植面积 估算模型，对地市及县两级行政单元的粮食作物种植面积进行估算，具体包 括如下步骤：

步骤P501：对步骤P3和步骤P4中得到的同类地物数据进行镶嵌,重叠 区域像元值选择最大值，得到包含纯像元和混合像元的作物分布。根据不同 像元的面积比例，统计不同地市及县两级行政单元内不同粮食作物种植的像 元数据，并根据像元大小计算不同地市及县两级行政单元内各种作物的毛种 植面积。

步骤P502：利用各地市及县两级行政单元内的耕地系数（耕地占图斑面 积的比例），估算地市及县两级行政单元内各种作物的净种植面积，即：

净种植面积=毛种植面积*耕地系数，

步骤P503：将各地市及县两级行政单元农作物面积求和，利用遥感方法 得到全省不同农作物种植面积，并结合全省耕地面积，计算全省范围内不同 农作物种植成数

统计各样方内分类得到的农作物面积，计算样方范围内农作物种植成数 然后对各个样方内地面调查得到的农作物面积与分类得到的农作物面积 进行回归分析，计算回归方程，得到回归系数b。

如在图1中，b=1.105；如在图3中，b=0.856；如在图4中，b=0.917。 图1、3、4中，验证样方的大小为1*1KM，因此根据这个验证结果，我们可 以得出我们的估算方法在1*1KM的尺度已经具有很高的精度，针对更粗的县 级尺度、地级尺度、省级尺度来说，估算精度显然是可以充分保证的。

在遥感识别估算领域，由于南方地区地块比较破碎，种植结构复杂，且 水稻种植分为早、中、晚三季，中稻和早晚稻都有生长时间重合的情况，其 种植状况较其他作物明显复杂，因此，水稻种植面积是粮食作物中公认的最 难以进行识别估算的，由图1、图3和图4所示，本实例技术方案能够对早 稻、中稻和晚稻的种植面积进行很准确地估算，因而本实例的技术方案也能 够对除水稻之外的其它粮食作物种植面积进行更加准确的估算。

利用抽样方法和分类方法估算的全省各农作物种植面积，建立不同农作 物种植面积回归估算模型：

$\overline{y_{\mathrm{reg}}}=\bar{y}+b(\overline{p_{\mathrm{pop}}}-\bar{p}),$

其中，回归估算模型修正的作物种植成数，为地面调查得到的作物 种植成数，b为修正系数，为分类方法统计的全省范围内农作物种植成 数，为分类方法统计的样方范围内农作物种植成数。

步骤P504：得到修正的作物种植成数后，由各地市及县两级行政单 元耕地面积数据计算各地市及县两级行政单元不同农作物种植面积。

上述实施例仅仅是为清楚地说明本发明创造所作的举例，而并非对本发 明创造具体实施方式的限定。对于所属领域的普通技术人员来说，在上述说 明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有 的实施方式予以穷举。凡在本发明的精神和原则之内所引伸出的任何显而易 见的变化或变动仍处于本发明创造权利要求的保护范围之中。