CN107527014A - 县级作物种植面积遥感统计抽样调查方案设计方法 - Google Patents

Abstract

本发明公开一种县级作物种植面积遥感统计抽样调查方案设计方法，其以多期现势高分一号16米宽幅影像为数据源，通过构建决策树分类方法自动解译得到研究区作物空间分布情况；然后以作物空间分布数据为抽样总体，结合研究区耕地面积及作物种植规模，采取PPS抽样方法抽选普查区，在抽中的普查区设计多种样本单元大小、多种样本单元数量，并将样方单元尺寸、样方单元数量和抽样方法结合推算方法进行随机组合，并将这些不同组合的统计抽样调查方案推算结果与高分二号0.8米分辨率影像数据解译得到的目标作物种植面积作比较；从结果数据的总体抽样精度、抽样稳定性以及最小样本量进行分析，并得出最适用于研究区目标作物遥感统计抽样调查方案。

Description

县级作物种植面积遥感统计抽样调查方案设计方法

技术领域

本发明涉及测绘遥感技术领域，具体涉及一种县级作物种植面积遥感统计抽样调查方案设计方法。

背景技术

传统的统计抽样已经被成熟的运用到业务化的运行中，但是在很多方面并不能满足现代统计业务的需求。一是传统的统计抽样业务是以基准数据为总体的，但基准数据往往来源于农业普查等大型国情国力调查，数据滞后性较强，影响统计抽样调查结果；二是传统的统计抽样大多采用目录抽样的方式进行抽选样本，样本本身并不具备地理空间信息属性，通过该方法得出的作物面积数据也不具备空间信息，更无法直观的呈现出作物的空间分布情况。

遥感技术作为获取农作物种植信息的空间技术，可以在很短的时间内获得大范围的农作物空间分布信息，可直观的对农作物空间分布地理空间信息进行展示，所获取的农作物空间分布信息也是当季最新数据。但是，单纯利用遥感技术来获取农作物面积信息也存在多方面的缺陷。受影像分辨率、天气因素以及解译方法等多方面的限制，遥感技术对于获取大的地物类别，如水体、建准用地、林地等可达到较高的精度，但对于具体作物类别的识别，如水稻、玉米等，这些作物类别往往表现为季节动态性，并且各种作物之间还存在同谱异物、同物异谱以及混合像元的现象，利用遥感技术进行高精度识别的难度比较大，获取的精度通常在70％左右，远远不能满足统计业务中95％以上精确度的需求，因此单纯利用遥感手段获取农作物种植面积也不能成熟的进行业务化运行。

虽然遥感所获取的农作物分布信息准确度并不是很高，但是通过遥感手段获取的农作物空间分布信息可以作为统计抽样的基准数据，因为这些分布数据是当年最新的分布数据，且具有空间分布信息，可直观的呈现出各种作物的空间分布情况，对抽样调查提供抽样依据以及便于抽选出具有空间属性的抽样单元，在很大程度上可减轻野外调查的工作强度。所以，将抽样技术理论与遥感技术结合进行农作物面积测量可进行很好的优势互补。

基于遥感与统计抽样相结合的农作物面积测量将是未来统计业务发展的趋势。但是目前，很多基于遥感与统计抽样相结合的农作物面积测量的研究都是以全国、省市为研究范围，尺度太大，往往得出的统计数据不够精确，无法获取县及县以下的农作物种植面积信息，以县级尺度范围进行的研究较少。另一方面，基于遥感的统计抽样研究并不系统化，很多研究都是以单一尺寸的样方单元大小或者是单一数量的样方个数或者是单一的抽样方法进行抽样研究，抽样外推方法也很单一。且由抽样反推所得出来的农作物种植面积数据也缺少系统的分析，单纯以总体反推误差进行分析，并未从总体反推精度、抽样稳定性以及抽样效率等多方面进行评价。

发明内容

本发明针对传统的统计抽样方法能够成熟的开展业务化运行的优势与缺少最新的基准数据和缺少样方空间属性的劣势，结合遥感技术能够全面短时间内获取含有空间信息的当年基准数据的优势，以各种尺寸的样本单元大小，不同数目的样本单元数量，结合简单随机抽样、系统抽样及分层抽样三种抽样方法，并借助直接推算以及回归模型推算方法，全面、系统的研究出一种县级作物种植面积遥感统计抽样调查方案设计方法，为及时、精准获取县级作物种植面积统计数据提供技术支持。

为达到上述技术目的，本发明的技术方案提供一种县级作物种植面积遥感统计抽样调查方案设计方法，所述县级作物种植面积遥感统计抽样调查方案设计方法包括以下步骤：

S1、基于自资源卫星遥感应用中心获取的高分二号卫星0.8米分辨率的中分遥感影像数据为准真值数据，并人工目视解译目标作物并对图斑属性赋值得到研究区目标作物空间分布情况；

S2、基于自资源卫星遥感应用中心获取的高分一号卫星16米分辨率的中分遥感影像数据为准真值数据，通过构建决策树分类方法自动解译得到研究区目标作物空间分布情况；

S3、采取PPS抽样方法抽选出多个普查区，在抽中的普查区设计多种样本单元大小、多种样本单元数量来进行目标作物种植面积统计抽样调查研究，分别采用简单随机抽样、系统抽样及分层抽样3种不同的抽样方法对样本单元进行抽选；以及分别采用直接推算法、回归模型推算法推算得到目标作物种植面积；

S4、将不同的样方单元尺寸、不同的样方单元数量和不同的抽样方法结合不同的推算方法进行随机组合，分析每一个因素对推算结果精度的影响；并将这些不同组合的统计抽样调查方案推算结果与基于高分二号卫星0.8米分辨率的中分遥感影像数据融合影像目视解译得到的目标作物种植面积作比较；从结果数据的总体抽样精度、抽样稳定性以及最小样本量分析各种统计抽样调查方案的优缺点以及抽样调查方案的适用条件，并得出最适用于研究区目标作物遥感统计抽样调查方案。

优选的，所述步骤S1包括以下分步骤：

S11、自资源卫星遥感应用中心获取高分二号卫星0.8米分辨率的中分遥感影像数据为准真值数据，对获取的遥感影像进行解压、辐射定标、大气校正、正射校正以及影像融合处理；

S12、通过野外实地调查工作，为亚米级影像数据人工目视解译工作建立真实可靠的目标作物解译知识库。

S13、根据影像中地块表现出的颜色、形状、大小、纹理、阴影等特征以及测量指标来确定所勾绘地块的边界，勾绘出的地块图斑矢量数据需要建立属性表。

S14、人工目视观察地块矢量边界内影像，并结合建立的解译知识库以及历年的农业统计数据综合进行解译识别，确定矢量地块内目标作物的具体类别并赋予属性，直至所有的地块完成解译识别，通过对解译得出的数据进行统计得到目标作物种植面积的准真值数据。

优选的，所述步骤S2包括以下分步骤：

S21、自资源卫星遥感应用中心获取高分一号卫星16米分辨率的中分遥感影像数据为准真值数据，对获取的遥感影像进行解压、辐射定标、大气校正以及正射校正；

S22、结合公安县的物候历数据以及不同时期影像中各种地物的形状、大小、颜色、纹理、阴影、图形及立体外貌等影像解译标志，初步确定各个时期遥感影像中所包含的典型地物的特征形态；

S23、通过构建决策树分类方法自动解译得到研究区目标作物空间分布情况。

优选的，所述步骤S23包括以下分步骤：

S231、选取目标作物抽穗时期的16米分辨率的中分遥感影像数据，根据影像中的NDVI和EVI值初步将影像中植被区与非植被区进行区分；

S232、选取目标作物移栽期的16米分辨率的中分遥感影像数据，根据近红外的光谱值将抽穗期影像中的山地林地地物进行剔除；

S233、选取目标作物抽穗期的16米分辨率的中分遥感影像数据，根据可见光红波段反射率值将目标作物与居民地、道路、裸露地进行区分并剔除；

S234、根据目标作物的移栽期、抽穗期、收获期三个时期遥感影像的EVI值初步提取出目标作物的大致分布；

S235、根据移栽期的遥感影像的近红外波段，将混杂在目标作物中的一些干扰目标作物进行剔除。

优选的，所述步骤S4中所述的将不同的样方单元尺寸、不同的样方单元数量和不同的抽样方法结合不同的推算方法进行随机组合，分析每一个因素对推算结果精度的影响；包括以下几种分析情况：

①当样方尺寸、样方数量以及抽样方法固定时，推算方式对抽样调查精度的影响：

②当样方尺寸、样方数量以及推算方法固定时，抽样方法对抽样调查精度的影响：

③当样方数量、抽样方法以及推算方法固定时，样方尺寸对抽样调查精度的影响：

④当样方尺寸、抽样方法以及推算方法固定时，样方数量对抽样调查精度的影响：

本发明所述的县级作物种植面积遥感统计抽样调查方案设计方法，其相较现有技术具有以下优势：

1、为及时、精准获取县级作物种植面积统计数据提供技术支持。

2、统计业务更加精细化。以往的统计业务是以省级、区域甚至是全国范围作为统计对象，尺度很大，不够精细化。本发明以县级尺度作为研究范围，在确保抽样调查方案更加合理的同时，实现了统计业务向县级尺度的精细化转变。

3、抽样调查方案的设计更加系统化、全面化。以往的调查方案的设计往往只简单地考虑到了样本单元大小或者抽样方法等一两个因素，本发明设计抽样调查方案则是从样本大小、样本数量、抽样方法以及推算方法四个维度对抽样调查方案进行全面系统化的设计。

4、抽样调查方案的评估更加合理化。以往调查方案的评估简单地以总体估算误差作为评估指标，本发明评估抽样调查方案则是以总体估算误差、抽样稳定性以及抽样效率三个指标进行综合评估，更加合理化。

5、减少野外实地调查工作量。传统的统计抽样缺少最新的基准数据和缺少样方空间属性的劣势，本发明结合遥感技术能够全面短时间内获取含有空间信息的当年基准数据的优势，在一定程度上可以减少野外实地调查工作量。

附图说明

图1为本发明所述县级作物种植面积遥感统计抽样调查方案设计方法的步骤流程图；

图2为本发明所述作物高分准真值数据获取的步骤流程图；

图3为本发明所述抽样方案设计与评估筛选的步骤流程图；

图4为本发明所述总体估算误差随样方尺寸与数量变化的曲线图；

图5为本发明所述标准差随样方尺寸与数量变化的曲线图；

图6为本发明所述估算误差绝对值总和与标准差总和随样方尺寸的变化曲线图；

图7为本发明所述估算误差绝对值总和与标准差总和随样方数量的变化曲线图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明涉及一种县级作物种植面积遥感统计抽样调查方案设计方法，如图1所示，所述县级作物种植面积遥感统计抽样调查方案设计方法包括以下步骤：

S1、基于自资源卫星遥感应用中心获取的高分二号卫星0.8米分辨率的中分遥感影像数据为准真值数据，并人工目视解译目标作物并对图斑属性赋值得到研究区目标作物空间分布情况；

S2、基于自资源卫星遥感应用中心获取的高分一号卫星16米分辨率的中分遥感影像数据为准真值数据，通过构建决策树分类方法自动解译得到研究区目标作物空间分布情况；

S3、采取PPS抽样方法抽选出多个普查区，在抽中的普查区设计多种样本单元大小、多种样本单元数量来进行目标作物种植面积统计抽样调查研究，分别采用简单随机抽样、系统抽样及分层抽样3种不同的抽样方法对样本单元进行抽选；以及分别采用直接推算法、回归模型推算法推算得到目标作物种植面积；

S4、将不同的样方单元尺寸、不同的样方单元数量和不同的抽样方法结合不同的推算方法进行随机组合，分析每一个因素对推算结果精度的影响；并将这些不同组合的统计抽样调查方案推算结果与基于高分二号卫星0.8米分辨率的中分遥感影像数据融合影像目视解译得到的目标作物种植面积作比较；从结果数据的总体抽样精度、抽样稳定性以及最小样本量分析各种统计抽样调查方案的优缺点以及抽样调查方案的适用条件，并得出最适用于研究区目标作物遥感统计抽样调查方案。

下面以湖北省公安县中稻种植面积遥感统计抽样调查方案设计为例做详细的阐述，具体包括以下步骤：

(1)作物高分准真值数据获取，如图2所示，包括以下步骤：

影像获取及预处理；将来源于资源卫星遥感应用中心的高分二号卫星0.8米分辨率的中分遥感影像数据作为公安县中稻准真值数据并进行解译。根据公安县中稻整个生长期在五月中下旬至十月中旬左右，获取公安县对应时间段的高分二号亚米级遥感影像数据，同时还获取了公安县的县级、乡镇级以及村级最新行政区划数据。由于从中国资源卫星应用中心获取的遥感影像数据为原始的压缩包影像数据，因此需要对获取的遥感影像进行解压、辐射定标、大气校正以及正射校正处理。辐射定标是将传感器获取的DN值转化成具有实际应用意义的大气顶层反射率；大气校正是为了减弱甚至消除大气因素对传感器获取的地物反射率的影响；正射校正是为了纠正因地形起伏或传感器因素引起的几何误差；辐射定标运用的是ENVI5.3中提供的Apply Gain and Offset工具进行定标的，需要设定相应的Gain Values和Offset Values参数；大气校正运用的是ENVI 5.3中提供的FLAASHAtmospheric Correction工具；正射校正使用的是ENVI 5.3中提供的RPCOrthorectification Workflow工具。全色遥感影像数据及多光谱影像数据做完辐射定标、大气校正以及正射校正处理后，需要进行影像融合处理。ENVI 5.3中的NNDiffuse融合算法是最新的影像融合算法，融合后的影像效果相对于其他几种融合算法效果更好，所以，应用ENVI 5.3中提供的NNDiffuse Pan Sharpening融合工具对全色以及多光谱影像数据进行融合处理。将融合后的遥感影像数据结合公安县县级行政区划矢量数据进行裁剪处理，裁剪后的影像数据和基准影像数据做影像配准处理。遥感影像配准是为了是使不同时间或者不同条件下获取的遥感影像中的同名点的地理空间坐标一致，也可以进一步消除遥感影像中畸变导致的地理空间误差。遥感影像配准可以应用ENVI 5.3中提供的ImageRegistration Workflow工具，配准时可以自动生成配准点也可以自行添加删除配准点，但是需要将配准点误差要控制在一个像元范围以内并尽量接近零值。至此，遥感影像数据的准备工作完成，可以直接应用于遥感影像的解译工作中。

外业调查建立解译知识库；野外实地调查工作主要是为了建立解译知识库而进行的实地调查采样的过程，结合调查所采集的中稻实地样本数据，可方便的识别出对应地理坐标处遥感影像中中稻所呈现出的颜色、形状、大小、纹理、阴影等特征，从而为亚米级影像数据人工目视解译工作建立真实可靠的中稻解译知识库。样方实地调查时，需要携带好相应的配套设备，如调查所需的PAD需安装外业电子数据包，包括调查地块以及周边的遥感影像和电子地图数据等。实地调查时，各工作人员按照野外调查前规划好的路线，到达指定的样本村范围，根据调查PAD具备的导航及定位功能导航到需要调查的具体样方中，到达调查样方时，需要根据样方的实地情况，将野外调查专题图以及PAD任务包中的样方地块分布图对应到实地中的田块。在当地有经验的人员指导下，对PAD上每一个矢量地块赋属性并且拍照，如果调查PAD中的矢量地块的边界与实地地块的边界不符合，则可以使用调查PAD软件中的编辑功能对矢量地块的边界进行编辑，编辑完成后进行保存并继续录入相应地块的属性，对于一些不便于人工调查的特殊样方，可借助无人机来辅助外业调查的整个过程。野外实地调查数据首先要保证数据的准确性，需要仔细检查每一个实地样方地块的属性是否准确，是否有遗漏以及属性填写错误的情况。其次，将野外调查获取的矢量数据叠加到遥感影像底图上，验证各个地块中的作物是否符合影像的光谱特征以及随时间变化的特征。对野外实地调查数据进行检验修改后，需要将野外调查数据进行整理以及归纳，以便建立解译知识库。解译知识库的建立可以为人工目视解译中稻提供很好的先验知识与训练样本。解译知识库中不仅包含中稻，还涵盖了许多影像中的其他地物类别，如水体、林地、建筑用地、裸露地、玉米等地物。解译知识库文件包含两类文件，一类是带有属性字段的记录了地块属性的矢量shp格式的文件，这类文件可以导入到ArcGIS软件中并与遥感影像进行叠加显示，具备很强的参考意义。另一种是依据第一类文件叠加显示后，进行数据整理，用图片以及文字对相应的地物类别进行描述的文件。

地块图斑的勾绘；以高分二号卫星0.8米分辨率的中分遥感影像数据作为图斑勾绘的基础影像，勾绘图斑的过程中，需要综合根据影像中地块表现出的颜色、形状、大小、纹理、阴影等特征以及测量指标来确定所勾绘地块的边界。在ArcGIS软件中勾绘地块边界时，要求勾绘的地块图斑闭合，公共边只需矢量化一次，地块之间不能存在重叠、自相交、空隙等拓扑错误；勾绘的地块矢量边界应尽量平滑，不能出现明显的尖角，除耕地地块外，房屋、水系、道路、林地等其他要素也必须进行勾绘。其中，当道路、水渠的宽度大于等于2米时，需要用多边形单独勾绘出来；当道路水渠等地块之间的边界小于2米时，只需要用单线条表示，作为地块间的分割线即可；当需要勾绘的地块图斑面积小于0.1亩时，该地块图斑不需要进行勾绘。勾绘出的地块图斑矢量数据需要建立属性表。

人工目视解译目标作物并对图斑属性赋值；解译知识库文件以及地块图斑矢量文件准备好后，便可进行人工目视解译，将高分二号亚米级影像数据和勾绘的公安县地块图斑数据都加载到ArcGIS软件中进行叠加显示，此时地块图斑样式需要调整为空心样式以便人工目视解译。按顺序选中其中某个地块矢量边界，人工目视观察该地块矢量边界内影像所表现出的地物的颜色、大小、形状、纹理等特征并结合建立的解译知识库以及公安县历年的农业统计数据综合进行解译识别，确定矢量地块内作物的具体类别。确认完成后，打开该地块的属性表，在属性表中的地块图斑类型中填写该地块的对应属性值并保存。按照地块的顺序依次进行解译识别并赋予属性直至所有的地块完成解译识别。保存解译识别后的地块图斑矢量数据，之后对解译得出的数据进行统计，得到公安县中稻种植面积的准真值数据。之后，利用ArcGIS软件对亚米级影像中稻目视解译得出的公安县中稻种植面积准真值数据进行统计，得出公安县2016年中稻种植面积准真值数据为54.3693万亩。

(2)基于中分遥感数据的研究区目标作物分布提取。

影像获取及预处理；基于中分遥感数据的公安县中稻分布提取解译应用的是从资源卫星遥感应用中心获取的高分一号卫星16米分辨率的中分遥感影像数据。要对获取的遥感影像进行解压、辐射定标、大气校正以及正射校正处理。辐射定标运用的是ENVI 5.3中提供的Radiometric Calibration工具进行定标的；大气校正运用的是ENVI 5.3中提供的FLAASH Atmospheric Correction工具；正射校正使用的是ENVI 5.3中提供的RPCOrthorectification Workflow工具。通过公安县县级行政区划矢量数据对经过辐射定标、大气校正以及正射校正处理的遥感影像数据进行裁剪处理。裁剪后的影像数据还需要和基准影像数据做影像配准处理。遥感影像配准是为了是使不同时间或者不同条件下获取的遥感影像中的同名点的地理空间坐标一致，也可以进一步消除遥感影像中畸变导致的地理空间误差。遥感影像配准可以应用ENVI 5.3中提供的Image Registration Workflow工具，配准时可以自动生成配准点也可以自行添加删除配准点，但是需要将配准点误差要控制在一个像元范围以内并尽量接近零值，保证每景影像数据的没有几何偏差，同名点能够重合。至此，通过对中分遥感数据进行解译后提取的中稻种植空间分布情况可直接使用。公安县中稻的移栽期大致在五月下旬至六月上旬，抽穗期在七月下旬至八月上旬，乳熟期在九月中下旬左右。其在移栽期、抽穗期及乳熟期三个关键时期中，在遥感影像中所呈现的各种特征参量以及在各个波段上的DN值都会产生规律性的变化。结合多种特征参量以及中稻和其他农作物物候期的不同，可以逐步剔除非中稻作物，再结合中稻本身的光谱值及组合参量随时间的时序变化规律，可以较高精度的提取出中稻的种植分布。基于专家知识的决策树分类方法也是目前最为先进的遥感影像分类方法，具有灵活、直观、效率高等优势，结合各种参量时序变化的规律，可以进一步提高中稻分类解译的精度。

中稻解译分类的基本步骤有：

初步判读；结合公安县的物候历数据以及不同时期影像中各种地物的形状、大小、颜色、纹理、阴影、图形及立体外貌等影像解译标志，初步确定各个时期遥感影像中所包含的典型地物的特征形态[46]。

区分植被与非植被区域；选取中稻抽穗时期的GF116米影像，分别得到影像的NDVI和EVI值。由于植被在可见光红波段反射率很低，在近红外波段反射率很高，所以植被相对于非植被，其NDVI与EVI值较高，对NDVI与EVI设定合理阈值(由于土壤中水分以及制备自身长势较弱等因素，EVI值或NDVI值会比正常的偏低，所以此时区分植被与非植被事，EVI或NDVI值可以稍微取小一些)，可以初步将影像中植被区与非植被区进行区分。

剔除山地林地；选取中稻移栽期影像，此时中稻处于刚移栽或者尚未移栽时期，植被信息很微弱，近红外光谱值相对于山地林地等植被信息较强的地物明显要偏低很多，此时对近红外的光谱值设定一个合理的阈值，可将抽穗期影像中的山地林地地物进行剔除，从而解译提取出种植有农作物的耕地。

进一步剔除居民地、道路、裸露地；第一步中区分植被与非植被地物只利用到了EVI或NDVI这个条件，但由于影像质量不佳或下雨积水等因素，植被与非植被的EVI或NDVI的值部分有交集，单纯对NDVI或EVI设定阈值并不能完全区分出植被以及非植被区域，故还需进一步剔除混在植被中的一些非植被地物(如居民地、道路、裸露地)，可通过中稻抽穗期影像上各种地物的光谱反射率值来区分，居民地、道路、裸露地的光谱反射率值在各个波段上都较高，而中稻的光谱反射率在可见光红波段上较低。所以，居民地、道路、裸露地与中稻的光谱反射率值在可见光红波段差异较大，对可见光红波段反射率值设定阈值，可以将中稻与居民地、道路、裸露地进行区分并剔除。

对种植农作物的耕地区域的作物品种进行进一步的分类解译；通过上述步骤，对中稻抽穗期影像中非植被及植被中的山地林地、长绿地等地物进行了剔除，抽穗期影像中提取出的地物就是种植农作物的耕地区域。结合公安县的物候历等数据，此时遥感影像中的耕地作物主要有棉花、大豆、玉米、中稻等。

初步提取出中稻；中稻移栽期，由于此时中稻处在刚移栽不久或即将移栽的时期，中稻的植被信息很弱，中稻地的植被覆盖度也很低，所以在中稻移栽期的遥感影像中，记录的更多是中稻的土壤背景，即含有大量水分的土壤。中稻抽穗期，这是中稻生长最为旺盛的一个时期，此时中稻的植被信息也最为强烈，光合作用最为旺盛，中稻种植地块上植被覆盖度也较高，所以在中稻抽穗期的遥感影像中，记录的更多是中稻本身的植被信息。中稻收获期，由于此时中稻处于即将收获或已经收获的时期，相比于中稻抽穗期，此时中稻的植被信息明显在减弱，光合作用也明显在减弱，中稻的整体颜色都呈现出偏黄的颜色，所以在中稻的收获期的遥感影像中，记录的更多是成熟期中稻的植被信息或者是中稻收获后的土壤信息。根据不同时期中稻在影像中呈现出的不同特征，可以对三个时期遥感影像的EVI值设定阈值，即中稻移栽期的EVI值小于中稻抽穗期的EVI值，中稻抽穗期的EVI值大于中稻成熟期的EVI值。通过对三个时期的EVI值设定相应的条件，可以初步提取出中稻的大致分布。

剔除干扰作物中的玉米、棉花以及大豆；通过步骤(5)中的参量，可以初步将中稻的分布提取出来，但是此时提取出的中稻地里还混杂着部分非中稻作物，比如玉米、棉花及大豆等。选取中稻移栽时期的遥感影像，由于中稻处于移栽期，其土壤中必定含有大量的水分，水分对近红外波段具有较强的吸收作用，而玉米、棉花以及大豆，都属于旱地作物，其土壤中水分含量相对于中稻种植地块的土壤水分含量明显要少，对近红外波段的吸收作用也不强。对移栽期的遥感影像的近红外波段设定合理的阈值，可以进一步将混杂在中稻中的一些干扰作物进行剔除，进一步提高提取中稻的精度。

(3)重点村的抽选。

抽样单元的选取首先需要构建总体抽样框架，明确构建抽样框架的对象。第一阶段，以公安县所有的行政村作为抽样对象，采取PPS抽样方法抽选出公安县内20个样本村，这20个样本村将作为抽样框架构建的总体对象。

PPS抽样方法即以耕地面积规模成比例的抽样方法。基本思想是先统计出公安县内各个行政村的耕地面积，计算出各个行政村的耕地面积占公安县总耕地面积的比重，在进行PPS抽选样本村时，每个行政村耕地面积占比就代表着每个村在抽样时能够被抽中的抽样概率，即耕地面积占比越大，那么它被抽中的概率也就越大，耕地面积占比越小，那么它被抽中的概率也就越小。

公安县共有16个乡镇、328个村委会。将公安县国土调查数据中的耕地图斑及各村的行政边界都加载到ArcGIS软件中，利用ArcGIS软件中的分区统计工具统计出公安县各个行政村的耕地面积，并计算出各行政村的耕地面积占公安县总耕地面积的比重，按与耕地面积规模成比例的抽样方法计算各行政村的抽中概率，并抽选出20个样本村作为总体抽样框架构建的研究对象，这20个样本村分别为金龙桥村、东风村、毛家台村、同强村、保恒垸村、榨岭新村、红安寺村、和祥村、国胜村、合兴村、黄堤村、跃进村、勇合村、国庆村、石牌村、欣荣村、中伏桥村、碾子沟村、杨家码头村及甘家铺村。

(4)抽样框架的构建。

以抽选出的20个样本村作为构建总体抽样框架的对象。公安县耕地虽然比较平整，但是耕地地块的破碎度较大，大多数地块面积就在几亩左右甚至不足一亩。结合公安县耕地实际情况，设计出5种具有不同大小样本单元的抽样框架来进行公安县中稻种植面积统计抽样调查研究，5种抽样单元大小分别为100米×100米、150米×150米、200米×200米、250米×250米、300米×300米。

将抽选出的20个样本村的矢量边界加载到ArcGIS中，以20个村的行政边界为总体抽样范围，应用ArcGIS软件中的构建渔网工具，构建100米×100米、150米×150米、200米×200米、250米×250米、300米×300米五种不同大小的渔网网格，这些网格便是总体抽样框架。在应用ArcGIS软件构建抽样框架时，样本村行政边界边缘处会有部分网格的局部在行政边界之外，为了保持抽样单元的完整性，需要保留整个的抽样单元，不需要对行政边界外的网格做裁剪舍弃处理，不同样方单元尺寸的抽样框架所包含的抽样单元个数也不一致。

(5)样本单元的抽选。

分别采用简单随机抽样、系统抽样及分层抽样3种不同的抽样方法对样本单元进行抽选，抽中的样本单元数量设计有5种，分别是60个、80个、100个、120个和140个。

简单随机抽样选取样本单元：简单随机抽样是按照每一个抽样单元被抽中的概率一致的原则进行抽样的。在本发明中，采用了系统性的简单随机抽样方法对样方单元进行抽选。即对每一个所选取的样本村，都采取简单随机抽样的方式抽选相同的样方单元数量。如共需抽选100个样方单元，那么每个村都抽选相同数量即5个样方单元。所以，该抽样方法属于系统性的简单随机抽样。采取该抽样方法分别对不同样方单元尺寸的抽样框架抽选出相应数量的样方单元，实际操作中，对每个村中的抽样单元按照顺序进行从1到n的编号，利用生成随机数的软件生成相应的随机数，与生成随机数对应编号的样方单元即被抽中，按照这种方式对所需的样方单元进行抽样。

系统抽样是将抽样框架中所有的样方单元按照顺序进行编号，再按照某一固定的规则对样方单元进行抽选的。其中，最常用的规则就是对按顺序进行编号的样方单元进行等距抽样，通常需要根据样方单元的总数和需要抽选的样方单元的数量来确定抽样距离，然后按照相同的抽样距离对样本单元进行抽选。具体操作时，首先在ArcGIS软件中对抽样框架中所有的样方单元按1到N进行编号，根据需要抽选的60个、80个、100个、120个、140个五种样方数量，计算出抽样距离分别为int(N/60)、int(N/80)、int(N/100)、int(N/120)、int(N/1400)。然后，以随机的方式分别在int(N/60)、int(N/80)、int(N/100)、int(N/120)、int(N/1400)范围内抽选出第一个样方单元，在第一个样方单元的基础上根据抽样距离抽选出剩余的需要抽选的样方单元。

分层抽样是先将抽样单元总体按照某种特定的标志或者特征进行分层，层与层之间将具有较大的差异性且并不会相互交叉重叠，层内的抽样单元将具备较大的同质性。在各层中，分别独立的采取简单随机抽样或者系统抽样的方式抽选出不定数量的抽样单元，从每层中抽样选取的样方单元合起来便是总的需要抽选的样方单元。简单来说，就是将总的抽样单元按照相似性进行聚类，原理类似于非监督分类算法中的IsoData算法，聚类完成后，再在每一个类别中抽选相应数量的样方单元。

(6)研究区作物面积的推算。

直接推算法：简单随机抽样与系统抽样的权数设计，以及中稻种植面积的推算：权数的设计一般需要与目标变量即公安县中稻种植面积有高度的相关性。以中分遥感解译所得出的公安县中稻种植面积作为直接推算权数的设计依据。

对简单随机抽样与系统抽样的权数进行设计时，对同一种抽样形式，若为每个村设计不一样的权数，由于每个村所抽选出的样本单元数量存在较大差异，会导致推算误差较大，所以以所有的样本村作为一个整体设计权数可以避免这种差异性带来的抽样误差。设计权数时，首先需要设计出20个样本村的权数。方法为统计出不同抽样尺寸下抽样框架中中分解译所得出的中稻种植面积S_c。由中分遥感解译所得出的公安县中稻种植面积S与抽样框架中的中稻面积的比值S_c作为所有样本村的设计的权数W，如公式(1)所示。

样本村的权数设定后，样方单元的权数w也需要进行设定。样方单元权数设计是以抽样框架包含的样方单元的总量N与抽中样方单元数量n的比值N/n作为抽中样方单元的权数的。按照权数设计的原理，简单随机抽样以及系统抽样所设计的权数如表1所示。完成权数的设定后，公安县中稻的种植面积可以依据抽中样方单元中的中稻准真值数据进行推算。实际操作过程中，需要先统计出抽中样方单元内的中稻准真值面积，然后乘以相应的样方单元权数和样本村权数便可以推算出公安县的中稻种植面积。

表1直接推算设计的权数

分层抽样权数的设计以及中稻种植面积的推算：分层抽样权数的设计中样本村的权数与简单随机抽样和系统抽样的权数设计是一样的。分层抽样中各个样方单元的权数与该样方单元所在的层数的样方单元总容量Nh相关。层内样方单元容量占总的样方单元容量比例越大，样方单元的权数也就越大，两者成正比例。分层抽样推算公安县中稻种植面积如公式(2)所示：

公式中，W代表样本村的权数，S代表推算出的公安县的中稻种植面积，N代表抽样框中总的样方单元的个数，h代表分层的层数，N_h代表第h层抽样单元的容量，Y_h代表第h层抽样单元内中稻准真值面积的平均值。

回归模型推算法：回归模型推算研究区作物种植面积，首先以样方单元内中分解译得出的目标作物种植面积作为自变量，以样方单元内亚米级数据解译得到的目标作物种植面积准真值数据作为因变量，由于每次抽样的样方单元数量为60、80、100、120或140个，所以需要建立相应个数的自变量向量与因变量向量，由这两者建立线性模型，如公式(3)所示。

y＝xβ+ε 式(3)

这里，x是样方单元内中分解译得出的目标作物种植面积构成的自变量向量，y是对应样方单元内亚米级数据解译得到的目标作物种植面积准真值数据构成的自变量向量，β为回归系数构成的列向量。通过最小二乘法可以得到回归系数参数。得出回归系数后，按照公式便可以推算出研究区目标作物种植面积。

其中，W代表样本村的权数，与直接推算中的样本村权数一致，N代表抽样框架的样本单元容量，代表抽中的样方单元作物准真值的平均值，代表所有样方单元中分解译所得作物种植面积的平均值，代表抽中的样方单元内中分解译所得作物种植面积的平均值。

(7)推算结果精度评估与抽样方案筛选

如图3所示，不同的样方单元尺寸、不同的样方单元数量和不同的抽样方法结合不同的推算方法对公安县中稻种植面积遥感抽样方案进行了系统性的研究。不同的遥感抽样方案推算出了不一样的成果。对于成果的精度评估，要从不同的角度进行分析。

总体估算精度，即由抽样方案推算出的公安县中稻种植面积与公安县中稻准真值面积的差值的绝对值，代表着多次抽样推算出的公安县中稻种植面积与实际的公安县中稻种植面积的误差绝对值。

抽样稳定性，即对同一种抽样方案进行反复多次的抽样，多次抽样所得的成果之间的差异性，一般由多次抽样所得成果数据的标准差进行度量。标准差越小，说明抽样稳定性越好。抽样稳定性是衡量抽样方案好坏的非常重要的指标，因为在实际的遥感抽样统计调查中，不可能会反复多次的进行抽样，只会按照单次抽样数据进行统计抽样。抽样稳定性越好，对统计抽样误差的纠正也会更加的简便。

最少样本量，在能够达到同等精度的条件下，使样本量尽量少一些，能够一定程度上减少抽样调查所需成本以及外业调查强度，大大提高抽样调查的效率。

所以，对抽样方案精度的评估，将从总体估算精度、抽样稳定性以及最少样本量三个角度进行分析。如表2数据所示，记录了由不同抽样方法通过直接推算的方法所得成果的总体估算精度以及抽样稳定性(标准差)。表3记录了由不同抽样方法通过回归模型推算的方法所得成果的总体估算精度以及抽样稳定性(标准差)。图4所示为平均估算误差随样方尺寸与数量变化的曲线图，图5所示为标准差随样方尺寸与数量变化的曲线图。接下来就依据实验数据从不同的角度来分析不同的抽样方案的精度。

①当样方尺寸、样方数量以及抽样方法固定时，推算方式对抽样调查精度的影响：由表2与表3中的数据可以分析出，通过回归模型推算出的数据无论是在总体估算精度方面还是在抽样稳定性方面都不如通过直接推算得出的数据。回归模型推算出的公安县中稻种植面积与公安县中稻种植面积准真值数据相差较大。精度方面完全不能满足统计业务的需求。通过对回归模型推算的各个过程进行分析，结果表明，由于所设计的抽样单元尺寸最大值为300米，而中分遥感解译所应用的影像分辨率为16米，所以单个样方单元内中分解译所得的中稻种植面积准确性很低也很不稳定，样方单元中的中分解译面积与准真值数据之间的回归关系也并不明显，推算的面积误差也就相对较大。所以，通过回归模型推算的方式推算公安县中稻种植面积并不适合该方案。所以，在推算方式的选择上，直接推算的方式是更加适用的，关于不同抽样方案的精度的分析数据将直接采用直接推算得出的数据进行对比分析。

②当样方尺寸、样方数量以及推算方法固定时，抽样方法对抽样调查精度的影响：以直接推算所得数据进行分析。由表2的数据表明，通过简单随机的抽样方法抽选不同尺寸、不同数量的样方单元，25种抽样方法推算结果的总体估算误差绝对值总和是38.492万亩；通过系统抽样抽选不同尺寸、不同数量的样方单元，25种抽样方法推算结果的总体估算误差绝对值总和是28.031万亩；通过分层抽样抽选不同尺寸、不同数量的样方单元，25种抽样方法推算结果的总体估算误差绝对值总和是16.53万亩；由图4中的折线数据可以分析出，5种样方尺寸与5种样方数量总共构成的25种抽样条件，其中，在10种抽样条件下分层抽样的总体估算误差绝对值是最小的，在8种抽样条件下系统抽样的总体估算误差绝对值是最小的，在7种抽样条件下简单随机抽样的总体估算误差绝对值是最小的。所以，在总体估算精度这一方面，分层抽样的精度是优于系统抽样的，系统抽样的精度是优于简单随机抽样的。

同样，通过简单随机的抽样方法抽选不同尺寸、不同数量的样方单元，25种抽样方法推算结果的标准差总和是178.79万亩；通过系统抽样抽选不同尺寸、不同数量的样方单元，25种抽样方法推算结果的标准差总和是164.91万亩；通过分层抽样抽选不同尺寸、不同数量的样方单元，25种抽样方法推算结果的总体估算误差绝对值总和是51.22万亩。同样，由图5中的折线数据可以分析出，5种样方尺寸与5种样方数量构成的25种抽样条件中，在21种抽样条件下分层抽样的标准差是最小的，在4种抽样条件下系统抽样的标准差是最小的，在0种抽样条件下简单随机抽样的标准差是最小的。所以，在抽样稳定性这一方面，分层抽样的稳定性是远远优于系统抽样和简单随机抽样的，系统抽样的稳定性略好于简单随机抽样。

③当样方数量、抽样方法以及推算方法固定时，样方尺寸对抽样调查精度的影响：以直接推算所得数据进行分析。图6记录的是当样方尺寸一定时，在所有的抽样条件下平均估算误差的总和以及标准差的总和随抽样尺寸的变化。从图6中可以分析出，总体估算精度随着样方尺寸的增大并没有明显的线性提高，标准差随着样方尺寸的增大在逐渐较小，也就意味着随着样方尺寸的增大，稳定性方面越来越好。样方尺寸与抽样稳定性方面的相关性较高。

④当样方尺寸、抽样方法以及推算方法固定时，样方数量对抽样调查精度的影响：以直接推算所得数据进行分析，图7记录的是当样方数量一定时，在所有的抽样条件下平均估算误差的总和以及标准差的总和随抽样数量的变化。从图7中可以分析出，总体估算精度随着样方数量的增大并没有明显的线性提高，标准差随着样方数量的增大在逐渐较小，也就意味着随着样方数量的增大，稳定性方面越来越好。所以，与样方数量类似，样方尺寸与抽样稳定性方面的相关性较高。

表2成果数据的总体估算误差绝对值与标准差(直接推算)

位：万亩

表3成果数据的总体估算误差绝对值与标准差(回归模型推算)

单位：万亩

通过精度评估，可以分析出，对于推算方法的选择，直接推算好于回归模型推算；对于抽样方法的选择，分层抽样好于系统抽样以及简单随机抽样；对于样方尺寸与样方数量的选择，尺寸越大或者数量越多，抽样稳定性也相对更好，但总体估算精度与这两个因素相关性不大。公安县中稻遥感抽样调查方案的筛选，也需要从总体估算精度、抽样稳定性以及最小样本量三个角度进行综合筛选。首先，分层抽样在总体估算精度以及抽样稳定性两个方面明显优于系统抽样以及简单随机抽样，所以抽样方案首选分层抽样。其次，需要筛选出在某一精度条件下的最小样本量。由图4和图5中数据可以分析出，分层抽样条件下，样方尺寸150米、样方数量80个，样方尺寸150米、样方数量120个以及样方尺寸150米、样方数量140个这三种条件下标准差依次为0.576万亩、0.670万亩以及0.656万亩，总体估算误差依次为0.152万亩、0.038万亩、0.058万亩，

综合上述分析，样方尺寸设定为150米×150米，样方大小设定为80个再结合分层抽样以及直接推断的方式可获得95％以上的推算精度，且抽样稳定性较高(即多次抽样推算结果的标准差很小，仅为公安县2016年中稻播种面积54.37万亩的1.06％)，因此该方案最适用于公安县中稻种植面积遥感抽样调查。

以上显示和描述了本发明的基本原理、主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。

Claims (5)

1.一种县级作物种植面积遥感统计抽样调查方案设计方法，其特征在于，所述县级作物种植面积遥感统计抽样调查方案设计方法包括以下步骤：

S1、基于自资源卫星遥感应用中心获取的高分二号卫星0.8米分辨率的中分遥感影像数据为准真值数据，并人工目视解译目标作物并对图斑属性赋值得到研究区目标作物空间分布情况；

S2、基于自资源卫星遥感应用中心获取的高分一号卫星16米分辨率的中分遥感影像数据为准真值数据，通过构建决策树分类方法自动解译得到研究区目标作物空间分布情况；

S3、采取PPS抽样方法抽选出多个普查区，在抽中的普查区设计多种样本单元大小、多种样本单元数量来进行目标作物种植面积统计抽样调查研究，分别采用简单随机抽样、系统抽样及分层抽样3种不同的抽样方法对样本单元进行抽选；以及分别采用直接推算法、回归模型推算法推算得到目标作物种植面积；

S4、将不同的样方单元尺寸、不同的样方单元数量和不同的抽样方法结合不同的推算方法进行随机组合，分析每一个因素对推算结果精度的影响；并将这些不同组合的统计抽样调查方案推算结果与基于高分二号卫星0.8米分辨率的中分遥感影像数据融合影像目视解译得到的目标作物种植面积作比较；从结果数据的总体抽样精度、抽样稳定性以及最小样本量分析各种统计抽样调查方案的优缺点以及抽样调查方案的适用条件，并得出最适用于研究区目标作物遥感统计抽样调查方案。

2.根据权利要求1所述县级作物种植面积遥感统计抽样调查方案设计方法，其特征在于，所述步骤S1包括以下分步骤：

S11、自资源卫星遥感应用中心获取高分二号卫星0.8米分辨率的中分遥感影像数据为准真值数据，对获取的遥感影像进行解压、辐射定标、大气校正、正射校正以及影像融合处理；

S12、通过野外实地调查工作，为亚米级影像数据人工目视解译工作建立真实可靠的目标作物解译知识库。

S13、根据影像中地块表现出的颜色、形状、大小、纹理、阴影等特征以及测量指标来确定所勾绘地块的边界，勾绘出的地块图斑矢量数据需要建立属性表。

S14、人工目视观察地块矢量边界内影像，并结合建立的解译知识库以及历年的农业统计数据综合进行解译识别，确定矢量地块内目标作物的具体类别并赋予属性，直至所有的地块完成解译识别，通过对解译得出的数据进行统计得到目标作物种植面积的准真值数据。

3.根据权利要求1所述县级作物种植面积遥感统计抽样调查方案设计方法，其特征在于，所述步骤S2包括以下分步骤：

S21、自资源卫星遥感应用中心获取高分一号卫星16米分辨率的中分遥感影像数据为准真值数据，对获取的遥感影像进行解压、辐射定标、大气校正以及正射校正；

S22、结合公安县的物候历数据以及不同时期影像中各种地物的形状、大小、颜色、纹理、阴影、图形及立体外貌等影像解译标志，初步确定各个时期遥感影像中所包含的典型地物的特征形态；

S23、通过构建决策树分类方法自动解译得到研究区目标作物空间分布情况。

4.根据权利要求3所述县级作物种植面积遥感统计抽样调查方案设计方法，其特征在于，所述步骤S23包括以下分步骤：

S231、选取目标作物抽穗时期的16米分辨率的中分遥感影像数据，根据影像中的NDVI和EVI值初步将影像中植被区与非植被区进行区分；

S232、选取目标作物移栽期的16米分辨率的中分遥感影像数据，根据近红外的光谱值将抽穗期影像中的山地林地地物进行剔除；

S233、选取目标作物抽穗期的16米分辨率的中分遥感影像数据，根据可见光红波段反射率值将目标作物与居民地、道路、裸露地进行区分并剔除；

S234、根据目标作物的移栽期、抽穗期、收获期三个时期遥感影像的EVI值初步提取出目标作物的大致分布；

S235、根据移栽期的遥感影像的近红外波段，将混杂在目标作物中的一些干扰目标作物进行剔除。

5.根据权利要求1所述县级作物种植面积遥感统计抽样调查方案设计方法，其特征在于，所述步骤S4中所述的将不同的样方单元尺寸、不同的样方单元数量和不同的抽样方法结合不同的推算方法进行随机组合，分析每一个因素对推算结果精度的影响；包括以下几种分析情况：

①当样方尺寸、样方数量以及抽样方法固定时，推算方式对抽样调查精度的影响：

②当样方尺寸、样方数量以及推算方法固定时，抽样方法对抽样调查精度的影响：

③当样方数量、抽样方法以及推算方法固定时，样方尺寸对抽样调查精度的影响：

④当样方尺寸、抽样方法以及推算方法固定时，样方数量对抽样调查精度的影响。