CN106845806A - 农田种植状态的遥感监测方法和系统 - Google Patents

农田种植状态的遥感监测方法和系统 Download PDF

Info

Publication number
CN106845806A
CN106845806A CN201710010876.8A CN201710010876A CN106845806A CN 106845806 A CN106845806 A CN 106845806A CN 201710010876 A CN201710010876 A CN 201710010876A CN 106845806 A CN106845806 A CN 106845806A
Authority
CN
China
Prior art keywords
vegetation index
period
farmland
threshold value
pixel
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201710010876.8A
Other languages
English (en)
Inventor
张淼
吴炳方
曾红伟
张鑫
李名勇
郑阳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Remote Sensing and Digital Earth of CAS
Original Assignee
Institute of Remote Sensing and Digital Earth of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Remote Sensing and Digital Earth of CAS filed Critical Institute of Remote Sensing and Digital Earth of CAS
Priority to CN201710010876.8A priority Critical patent/CN106845806A/zh
Publication of CN106845806A publication Critical patent/CN106845806A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00Administration; Management
    • G06Q10/06Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063Operations research, analysis or management
    • G06Q10/0639Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393Score-carding, benchmarking or key performance indicator [KPI] analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/02Agriculture; Fishing; Forestry; Mining

Landscapes

  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Engineering & Computer Science (AREA)
  • Economics (AREA)
  • Strategic Management (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Theoretical Computer Science (AREA)
  • Educational Administration (AREA)
  • Marketing (AREA)
  • General Physics & Mathematics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Development Economics (AREA)
  • Physics & Mathematics (AREA)
  • Quality & Reliability (AREA)
  • Operations Research (AREA)
  • Game Theory and Decision Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Agronomy & Crop Science (AREA)
  • Animal Husbandry (AREA)
  • Marine Sciences & Fisheries (AREA)
  • Mining & Mineral Resources (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种农田种植状态的遥感监测方法和系统,涉及遥感技术领域。所述方法包括:获取遥感数据,经过预处理后得到像元的植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值;提取像元植被指数曲线的极值点,所述极值点包括极大值点和极小值点;和逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态为休耕或有作物种植。本发明利用中低分辨率遥感数据即可有效识别、监测某一时间段的农田种植状态的动态变化,既提高了识别精度,又降低了成本。

Description

农田种植状态的遥感监测方法和系统
技术领域
本发明涉及遥感技术领域,具体地说,涉及一种农田种植状态的遥感监测方法和系统。
背景技术
在农作物生长季内,掌握农田种植状态的动态变化信息意义重大。首先,农田种植状态的动态变化信息能够反映出不同作物的轮作模式以及不同时期农田得到有效利用的比例与强度等信息。连续多年的高强度耕作会导致农田土壤养分、土壤有机质等的含量逐渐下降。另外,在作物收获之后秸秆还田量较少的情况下,农田土壤肥力和产出也会不断下降。因而不同作物的轮流种植,甚至于适时适度的休耕对保持土壤肥力,提高作物单产十分必要。及时、客观、准确的农田种植状态信息,可以科学地指导农业生产。
其次,作物种植面积是农业生产中一项必要的管理信息,对作物种植面积估算的准确性对农业生产有着重要意义。目前,通常依赖农作物精细分类的方法来估算作物种植面积。但是,这种方法不足以解决作物种植面积中存在的问题,且众多农作物精细分类方法往往仅适用于对小范围或者农田田块较大的区域作物分类,尤其是在当前生长季尚未结束时,由于缺少全生长季的数据,从而使得作物类型的精细识别是一个短期内无法解决的难题。
再有,除了作物种植面积的估算,作物单产预测以及作物长势监测也是农业生产管理中必须要收集的信息。而及时、客观、准确的农田种植状态信息,能够在一定程度上改善相应的估算、预测与监测精度。
美国国家宇航局(NASA)于1991年发起了一个称为地球科学事业(ESE)的综合性项目,在该项目中的地球观测卫星系列(EOS)部分,有两颗重要的卫星:Terra卫星和Aqua卫星。其中,Terra卫星每天上午从北向南通过赤道,因此又被称为地球观测第一颗上午星(EOS-AM1)。Aqua卫星每天下午从南向北通过赤道,因此被称为地球观测第一颗下午星(EOS-PM1),两颗星均为太阳同步极轨卫星,在数据采集时间上相互配合。
中分辨率成像光谱仪(MODerate-resolution Imaging Spectroradiometer,简称MODIS)是Terra卫星和Aqua卫星上搭载的主要传感器之一,两颗星相互配合,每天可重复观测整个地球表面,得到36个波段的观测数据。MODIS自2000年4月开始正式发布数据,MODIS传感器获取的遥感数据因其在时空监测尺度上的优越性,被广泛地用于植被、土地利用状况的监测。
通常,高分辨率的遥感数据能够较中低分辨率遥感数据获取更高精度的农田种植状态识别及监测结果,但是也带来了高成本的问题,无法有效开展大范围的农田种植状态的识别和监测。中低分辨率遥感数据为实现大范围农田种植状态的识别和监测提供了可能性,降低了监测成本,但监测精度也较高分辨率数据有所降低。
发明内容
本发明要解决的技术问题在于,针对现有技术中采用高分辨率遥感数据监测作物时的成本高、而采用中低分辨率遥感数据时又达不到足够监测精度的问题,提供了一种农田种植状态的遥感监测方法和系统,利用中低分辨率遥感数据有效识别和监测作物生长季内的农田种植状态的动态变化,提高识别精度。
为解决上述技术问题,根据本发明的一个方面,本发明提供了一种农田种植状态的遥感监测方法,其中,包括:
获取遥感数据,经过预处理后得到像元的植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值;
提取像元植被指数曲线的极值点,所述极值点至少包括极大值点;
逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态。
优选地,所述植被指数阈值包括休耕农田阈值和种植作物农田阈值,在所述时期的植被指数值大于或等于种植作物农田阈值时,确定所述时期的农田种植状态为有作物种植;
在所述时期的植被指数值小于所述休耕农田阈值时,确定所述时期的农田种植状态为休耕;
在所述时期的植被指数值小于种植作物农田阈值、且大于或等于休耕农田阈值时,根据所述时期的植被指数数据的变化趋势,参考所述植被指数极大值与种植作物农田阈值的比较结果,确定所述时期的农田种植状态。
优选地,参考所述植被指数极大值与种植作物农田阈值的比较结果,确定所述时期的农田种植状态具体包括:
判断所述时期的植被指数数据的变化趋势,如果所述时期的植被指数数据变化趋势为增大,获取所述时期之后最邻近的极大值;比较所述像元植被指数极大值与所述种植作物农田阈值的大小;
如果所述时期的植被指数数据变化趋势为减小,获取所述时期之前最邻近的极大值;比较所述像元植被指数极大值与所述种植作物农田阈值的大小;
如果所述植被指数极大值大于或等于所述种植作物农田阈值,确定所述时期的农田种植状态为有作物种植;如果所述植被指数极大值小于种植作物农田阈值,确定所述时期的农田种植状态为休耕;
如果所述时期的植被指数数据变化趋势不变,确定所述时期的植被指数值为极大值或极小值,当所述时期的植被指数值为极大值时,则所述时期的农田种植状态为休耕农田;
当所述时期的植被指数值为极小值时,在所述极小值所在时期之前和所在时期之后分别获取与其相邻的极大值,分别对比前、后相邻的极大值与所述种植作物农田阈值的大小,如果前、后极大值中至少一个大于或等于所述种植作物农田阈值,所述时期的农田种植状态为有作物种植,如果所述前、后极大值均小于所述种植作物农田阈值,所述时期的农田种植状态为休耕农田。
优选地,通过以下步骤判断所述时期的植被指数数据的变化趋势:
在所述时期对所述像元植被指数曲线求导,当所述时期的像元植被指数导数大于0时,所述时期的植被指数数据变化趋势为增大;
当所述时期的像元植被指数导数小于0时,所述时期的植被指数数据变化趋势为减小;
当所述时期的像元植被指数导数等于0时,所述时期的植被指数数据变化趋势不变;
当所述时期植被指数数据变化趋势不变,则所述时期植被指数处于极大值或极小值点,利用所述时期前一时期与后一时期像元植被指数导数的正负判断,如果前一时期像元植被指数导数大于0,后一时期像元植被指数导数小于0,则所述时期植被指数处于极大值点,如果前一时期像元植被指数导数小于0,后一时期像元植被指数导数大于0,则所述时期植被指数处于极小值点。
优选地,在逐像元、逐时期比较植被指数值与植被指数阈值之前还包括:
根据训练样本训练得到所述植被指数阈值,所述植被指数阈值至少包括休耕农田阈值和种植作物农田阈值。
优选地,在获得特定区域的农田种植状态后,还包括:
统计特定时期、特定区域范围内休耕农田和有作物种植农田像元数量;
计算有作物种植农田像元数量占休耕农田与有作物种植农田像元数量总和的比例,确定所述特定时期内的农田种植率。
为解决上述技术问题,根据本发明的另一个方面,本发明提供了一种农田种植状态的遥感监测系统,其中,包括:
预处理模块,用于对获取的遥感数据进行预处理,获得植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值;
极值提取模块,用于提取像元植被指数曲线的极值点,所述极值点至少包括极大值点;和
数据处理模块,用于逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态。
优选地,所述数据处理模块包括:
数据读取单元,用于读取植被指数阈值和指定时期的植被指数值,其中,所述植被指数阈值包括休耕农田阈值和种植作物农田阈值;
比较单元,用于比较指定时期的植被指数值和休耕农田阈值的大小,比较指定时期的植被指数值和种植作物农田阈值的大小,及比较植被指数极大值与种植作物农田阈值的大小;和
识别单元,根据所述比较单元的比较结果,确定所述指定时期的农田种植状态;
优选地,所述数据处理模块包括:
导数计算单元,用于在指定时期对所述像元植被指数曲线求导获得所述指定时期的导数;
所述比较单元包括比较所述导数与0的大小。
优选地,所述系统还包括:
植被指数阈值训练模块,根据训练样本训练得到所述植被指数阈值,所述植被指数阈值至少包括休耕农田阈值和种植作物农田阈值;
和/或验证模块,用于根据验证样本验证数据处理模块识别得到的农田种植状态的精度。
优选地,所述系统还包括:
数据统计模块,用于统计特定时期、特定区域的休耕农田、有作物种植农田的像元数量,计算有作物种植农田像元数量占休耕农田与有作物种植农田像元数量总和的比例,作为所述特定时期、特定区域的农田利用率。
本发明充分利用时间序列NDVI数据的高时间频率特性,准确地获取了不同时期休耕农田与种植作物的农田反映在遥感影像上的显著差异,仅根据免费的中低分辨率的遥感数据即可实现高精度的农田种植状态的动态监测,成本低、精度高,且实施的流程标准,简单易行,具有推广性;识别出的作物生长季内的农田种植状态动态变化情况,为全球尺度作物种植面积估算提供了解决方案和数据支持,可推动作物种植面积估算创新性方法的诞生。
附图说明
图1为本发明实施例一提供的一种农田种植状态的遥感监测方法流程示意图;
图2为本发明实施例一逐像元、逐时期识别农田种植状态的流程示意图;
图3为本发明实施例一利用植被指数阈值建立决策树、识别像元Ni的时期Tj的农田种植状态的具体流程示意图;
图4为本发明利用休耕农田训练样本多边形对与观测时间相对应的NDVI数据进行分区统计,获得休耕农田样本的NDVI累计直方图;
图5为本发明利用作物种植农田的训练样本多边形对同时期的NDVI数据进行分区统计,获得种植作物的农田样本的NDVI累计直方图;
图6为本发明实施例二所述农田种植状态的遥感监测系统的原理框图;
图7为本发明实施例二所述数据处理模块的原理框图;
图8为本发明实施例二中另一所述农田种植状态的遥感监测系统的原理框图。
图9为本发明实施例三提供一种农田种植率的计算方法流程示意图;
图10为本发明实施例三中遥感监测系统的另一原理框图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例一
参见图1,为本发明实施例一提供的一种农田种植状态的遥感监测方法流程示意图。所述方法包括:
步骤S1,获取遥感数据。
步骤S2,对取得的遥感数据预处理,从而得到像元的植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值。
步骤S3,提取像元植被指数曲线的极值点,所述极值点至少包括极大值点。
步骤S4,逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态。
具体地,在步骤S1中,所述遥感数据来源于中分辨率成像光谱仪(MODerateresolution Imaging Spectro-radiometer,简称MODIS),包括MODIS上午星(Terra)和下午星(Aqua)16日合成归一化植被指数(The Normalized Difference Vegetation Index,简称NDVI)产品,产品编号为MOD13Q1和MYD13Q1,空间分辨率250米,时间范围是从2011年10月1日至2013年9月30日,共获得了五个不同空间块(h26v04,h26v05,h27v04,h27v05和h28v05)的长时间序列植被指数产品数据。该数据是从美国国家宇航局(NASA)的Reverb网络工具(http://reverb.echo.nasa.gov)上下载获取。MODIS植被指数产品是由16天内每天的植被指数采用最大值合成方式合成而来,其目标是优先选择近星下点无云像元,尽可能减小残存云、暗影、大气气溶胶和BRDF效应的影响,在仪器特性和地表特性的限制条件下尽可能增加空间和时间的覆盖度,同时保证合成资料的质量及一致性。
在步骤S2中,对取得的遥感数据预处理包括数据拼接、重投影、数据类型转换和数据格式转换等操作,原始数据采用HDF(Hierarchical Data Format)科学数据集(Scientific Data Sets,简称SDSs)的方式分块存储。首先利用MODIS重投影工具对五个不同空间块进行拼接,并从HDF文件中读取NDVI波段数据,直接读取出来的数据是16位有符号整型数据,需要将直接读取的数据除以10000以转换为NDVI实际值,有效范围为-1至1之间。最后利用MODIS重投影工具将NDVI数据投影为Albers110大地坐标系,采用WGS84椭球体,并将转换为GEOTIFF数据格式。此时获得了为识别使用的像元的植被指数曲线数据,由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值。
在步骤S3中,利用预处理后得到的像元NDVI曲线提取极值点,具体包括:
采用S-G平滑算法(Savitzky-Golay filter)逐像元对时间序列NDVI曲线进行平滑处理,具体地,将平滑窗口的大小设置为7,高阶多项式的阶数设置为3。
而后,采用牛顿插值多项式对时间序列NDVI曲线进行逐点拟合,选用二阶牛顿插值多项式对相邻三个NDVI点进行插值曲线拟合,依据拟合出的插值曲线计算该时期NDVI导数值,以NDVI曲线前三个点为例,牛顿插值多项式拟合和导数计算过程如下:
假设前三个NDVI点满足某个二次多项式函数N,已知它在三个点上的取值为:
Nn(xi)=NDVIi,i=0,1,2;····································(1)
由于横坐标对应的为时间且时间间距固定,因此为了简化起见,可将x0,x1,x2设置为0,1,2。
则二阶牛顿插值多项式可记为:
N2(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)·······(2)
上式中,包含a0,a1,a2三个未知系数,而利用连续三个点上的NDVI值便可求解出各系数:
a0=x0······························(3)
由此便获得了前三个NDVI点的二阶牛顿插值多项式。本实施例优选了二阶牛顿插值多项式,有效避免了龙格现象对插值结果的影响。
对如上获得的二阶牛顿插值多项式求解一阶导数,作为前3个点中第2个点的导数;除第一期和最后一期NDVI数据外,其他各期NDVI数据的导数均按上述方式计算;对于NDVI序列上的第一期数据,将最后一期NDVI前移作为三点插值的第一点,对于NDVI序列上的最后一期数据,将第一期NDVI后移作为三点插值的最后一点,利用同样的计算过程获得相应位置的NDVI导数。
对二阶牛顿多项式求导之后可得其一阶导数可记为:
N′2(x)=a1+2a2*x-a2(x0+x1)·······················(6)
将x0,x1,x2的值(分别为0,1,2)代入(6)式中,则获得了第2个点的导数只与a1,a2有关,而a1,a2可以直接由(4)和(5)式计算,因此大大简化导数计算过程,提升导数计算效率。
最后,利用二阶牛顿插值多项式的导数计划结果进行极值点的判定。当二阶牛顿插值多项式的导数由正值过渡为负值时,该转折点为对应像元的NDVI曲线的极大值点,由负值过渡为正值的转折点为NDVI曲线的极小值点,从而对应地获得该处的NDVI值。
在步骤S4中,所述植被指数阈值包括休耕农田阈值NDVIuncropped和种植作物农田阈值NDVIcropped。其中,本发明利用植被指数阈值建立决策树,逐像元、逐时期识别农田种植状态,具体流程如图2所示:
步骤S41,选取像元Ni的时期Tj,其中,i为像元编号,i=1,2,3……m,其中m为像元最大值;j=1,2,3……n,其中,n为像元时期的最大值,即前述NDVI曲线中时间序列的总数。
步骤S42,识别像元Ni的时期Tj的农田种植状态,识别完像元Ni的时期Tj的农田种植状态后,执行步骤S43。
步骤S43,判断是否识别完像元Ni的所有时期,如果没有识别完,则转到步骤S42,继续识别时期Tj的下一个时期Tj+1的农田种植状态,如果已经识别完,则执行步骤S44。
步骤S44,判断是否识别完所有的像元,如果都已经识别完,则完成了农田种植状态的识别,结束本流程。如果还有像元没有识别,则取像元Ni的下一个像元Ni+1,返回步骤S41继续上述流程,直到识别完所有的像元。
其中,识别像元Ni的时期Tj的农田种植状态的具体流程如图3所示。
步骤S420,比较时期Tj的NDVI值与NDVIcropped的大小,其中,通过NDVI曲线可以得到对应时期Tj的NDVI值。NDVIcropped为植被指数阈值,用于代表种植作物的农田,为预选设置的一个参数,该参数可以通过一些训练样本数据得到。
步骤S421,判断时期Tj的NDVI值是否大于或等于NDVIcropped,如果时期Tj的NDVI值大于或等于NDVIcropped,说明此时期作物种植;如果NDVI值小于NDVIcropped,转到步骤S422。
步骤S422,判断时期Tj的NDVI值是否小于或等于NDVIuncropped。NDVIuncropped为植被指数阈值,用于代表休耕农田,如果时期Tj的NDVI值小于或等于NDVIuncropped,说明此时期没有作物种植,为休耕农田;如果时期Tj的NDVI值大于NDVIuncropped,即时期Tj的NDVI值介于NDVIuncropped和NDVIcropped之间,此时,执行步骤S423。
步骤S423,在时期Tj对所述像元Ni的NDVI曲线求导,得到时期Tj的导数。
步骤S424,判断所述导数是否等于0,如果所述导数等于0,则说明时期Tj对应的NDVI处于极大值点或极小值点,则转到步骤S425,如果不等于0,执行步骤S429。
步骤S425,判断当前时期Tj对应的NDVI是否处于极大值点,如果是极大值点,则说明此时期的农田没有作物种植,为休耕农田。如果是极小值点,则执行步骤S426。
步骤S426,获取时期Tj前后相邻的极大值。
步骤S427,分别对比所述前、后相邻的极大值与NDVIcropped的大小。
步骤S428,根据对比结果,判断是否至少有一个极大值大于或等于NDVIcropped。如果至少有一个极大值大于或等于NDVIcropped,则确定所述时期的农田有作物种植。如果这两个极大值都小于NDVIcropped,则确定所述时期的农田为休耕农田。
步骤S429,判断所述导数是否大于0,如果大于0,说明所述时期的NDVI值在NDVI曲线上的处于上升阶段,即所述像元NDVI曲线在时间上呈上升趋势,此时的极大值应位于所述时期Tj之后,寻找该时期Tj之后的第一个极大值,此时转到步骤S431。如果所述导数小于0,说明所述时期的NDVI值在NDVI曲线上处于下降阶段,即所述像元NDVI曲线在时间上呈现减小趋势,此时的极大值应位于所述时期Tj之前,应到该时期Tj之前寻找时间上最接近的极大值,此时转到步骤S430。
步骤S430,获取时期Tj前最邻近的极大值,而后执行步骤S432。
步骤S431,获取时期Tj后最邻近的极大值,而后执行步骤S432。
步骤S432,比较所述极大值与NDVIcropped的大小。
步骤S433,判断所述极大值是否大于或等于种植作物农田阈值NDVIcropped,如果是,则确定所述时期Tj的农田种植状态为有作物种植,如果不是,所述极大值小于所述种植作物农田阈值NDVIcropped,所述时期的农田种植状态为休耕农田。
在前述方案中,所述植被指数阈值NDVIcropped和NDVIuncropped为一预设值,为提高识别准确率,基于地面观测获取的农田种植状态样本数据提取不同农田种植状态下植被指数的累计直方图,并确定所述植被指数阈值。
例如:在2012年至2013年每年4月初、5月底、7月初、8月底先后八次在黄淮海平原开展地面观测,利用手持GPS记录了各观测时期农田是否有作物种植以及种植的作物类型数据,并在实验室内将点状数据结合高分辨率遥感影像拓展为面状数据,用于农田种植状态监测的训练样本和验证样本,基于对实验区不同类型的农田种植状态实际分布情况,选择每种类型的50%作为训练样本,其余50%的样本作为分类后验证样本。
利用休耕农田训练样本多边形对与观测时间相对应的NDVI数据进行分区统计,获得休耕农田样本的NDVI累计直方图与累计频率(见图4);利用种植作物的农田训练样本多边形对同时期的NDVI数据进行分区统计,获得作物种植区样本的NDVI累计直方图与累计频率(见图5)。在该图4、图5中,横轴表示NDVI值,纵轴表示累计频率。
结合专家知识确定休耕农田和种植作物农田阈值,截取休耕农田NDVI累计频率为96%处的值作为休耕农田阈值(NDVIuncropped),在本实施例中为0.25(如图4所示);截取种植作物的耕地NDVI累计频率为4%处的值作为种植作物农田阈值(NDVIcropped),在本实施例中为0.6(如图5所示)。
利用如上所述的步骤监测了黄淮海平原的农田种植状态动态变化过程,该变化过程反映了2012年10月冬小麦播种之后至2013年9月底玉米开始收获期间黄淮海平原农田种植状态动态变化规律特征。2012年11月至2013年3月期间,黄淮海平原北部部分农田一直处于休耕状态,表明该时期内北部部分地区农田无作物种植,主要原因是北部地区冬季气温较低,同时耕地肥力不足,部分农户放弃种植冬小麦,每年仅种植一季玉米。从2013年5月开始,由于春玉米、棉花等单季作物均开始生长,无作物种植的耕地逐渐减少,2013年6月下旬开始,部分地区夏粮作物收获完成,秋粮作物尚未生长,因此约1/2的耕地上这一时期无作物种植,但河南东部、安徽江苏北部等地区秋粮播种较早,6月下旬秋粮作物已经开始生长,因此这些区域为种植作物的农田,之后黄淮海平原的农田几乎全部种植作物。
空间上,河南、江苏、安徽北部和山东(除黄河三角洲地区)2012年10月至2013年9月期间,除冬小麦-玉米轮作换季期间,农田无作物生长外,其余时期农田几乎得到100%的利用,全年种植两季作物,且第一季作物收获至第二季作物播种的时间间隔很短;在河南北部、西部山区和河北南部的衡水地区农田种植状态较低,冬季无作物种植,每年仅5月之后耕地上才种植单季作物。
实施例二
如图6所示,为本发明所述农田种植状态的遥感监测系统的原理框图。所述遥感监测系统包括预处理模块1a、极值点提取模块2a和数据处理模块3a。其中,所述预处理模块1a用于对获取的遥感数据进行预处理,获得植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值。所述极值提取模块2a用于提取像元植被指数曲线的极值点,所述极值点包括极大值点和极小值点。所述数据处理模块3a利用植被指数阈值建立决策树,通过逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态。
具体地,所述数据处理模块如图7所示,包括:数据读取单元31、比较单元32、分类单元33和导数计算单元34。其中,所述数据读取单元31用于读取植被指数阈值和指定时期的植被指数值,其中,所述植被指数阈值为预设值,包括休耕农田阈值和种植作物农田阈值。
所述比较单元32用于进行比较运算,例如:比较指定时期的植被指数值和休耕农田阈值的大小,比较指定时期的植被指数值和种植作物农田阈值的大小,以及,比较植被指数极大值与种植作物的耕地阈值的大小,并将比较的结果发送给所述分类单元33。
所述分类单元33根据所述比较单元的比较结果,确定所述指定时期的农田种植状态。具体地,参考图3,所述识别单元33根据植被指数值与种植作物农田阈值或休耕农田阈值的大小,确定所述时期的农田为休耕或有作物种植状态;而在所述时期的植被指数值小于种植作物农田阈值、且大于或等于休耕农田阈值时,根据所述时期的植被指数数据的变化趋势,也就是该时期植被指数曲线的导数是大于0、小于0还是等于0来确定所述时期的农田种植状态。因而,所述数据处理模块3包括导数计算单元34,用于在指定时期对所述像元植被指数曲线求导获得所述指定时期的导数。
所述比较单元32还用于比较所述导数与0的大小,在所述导数不等于0时,将比较结果发送给分类单元33。所述分类单元33接收到所述比较结果后,根据所述植被指数极大值与种植作物农田阈值的比较结果,如果所述植被指数极大值大于或等于种植作物农田阈值,则确定该农田有作物种植,否则确定该农田为休耕状态。
所述比较单元32在比较得到所述导数等于0时,通过比较像元的植被指数曲线数据,确定当前所述植被指数值为极大值还是极小值,如果为极大值,将该消息发送给分类单元33。所述分类单元33在接收到该消息后,确定所述时期的农田种植状态为休耕农田。
所述比较单元32通过比较确定所述时期的植被指数值为极小值时,发送指令给所述数据读取单元31。所述数据读取单元31在所述极小值所在时期之前和所在时期之后分别获取与其相邻的极大值,并发送给所述比较单元32,所述比较单元32分别对比前、后相邻的极大值与所述种植作物农田阈值的大小,并将对比结果发送给分类单元33。分类单元接收到该结果后,根据对比结果,如果前、后极大值中至少一个大于或等于所述种植作物农田阈值,所述时期的农田种植状态为有作物种植,如果所述前、后极大值都小于所述种植作物农田阈值,所述时期的农田种植状态为休耕农田。
进一步地,在另一种实施方式中,如图8所示,为了提高识别准确率,该实施方式中的系统除了包括预处理模块1b、极值提取模块2b和数据处理模块3b外,还包括植被指数阈值训练模块4b,为了验证本系统的识别精度,还包括验证模块5b,并包括用于训练所述植被指数阈值的训练样本数据和进行验证的验证样本数据,所述训练样本数据和验证样本数据没有重复和交集,相互独立。所述植被指数阈值利用训练样本数据利用累计直方图及累计频率确定,具体如前述实施例一中的方法所述,在此不再重复说明。
关于验证,在一个具体实施方式中,利用采样采集的验证样本数据对前述实施例中的黄淮海平原农田种植状态遥感监测结果进行精度验证,具体过程是,将分类结果与地面实测获取的地物类型进行比较,比较后导出分类结果的混淆矩阵,而后使用总体分类精度和Kappa系数对分类结果的精度进行评价。其中,总体分类精度与Kappa系数的计算公式如下:
总体分类精度:
Kappa系数计算公式:
式中,pii为混淆矩阵第i行,第i列对应的像元数,n代表分类结果的类别数量,N为所有用于精度评价的像元总数,pi+和p+i和分别是混淆矩阵第i行和第i列的总像元数。
举例来说,表1是根据某次观测时获取的地物真实分布数据导出的分类结果混淆矩阵:
表1 黄淮海平原农田种植状态混淆矩阵
其中,地面观测获得的实际农田种植状态样本中,共1890个像元为种植作物的农田,2255个像元为休耕农田;1890个种植作物的农田中,有1882个像元正确的识别为种植作物农田,8个像元被错误的归类为休耕农田,生产者精度高达99.6%;2255个休耕农田像元中,2062个像元被正确的识别为休耕农田,193个像元被错误的归类为种植作物的农田,休耕农田的生产者精度为91.4%。
从农田种植状况监测结果的角度,4145个像元范围内(1890+2255),共2075个像元被归类为种植作物的农田,其中1882个像元实际为种植作物的农田像元,用户精度为1882/2075=90.7%;共2070个像元被归类为休耕农田,其中2062个像元实际为休耕农田,用户精度为2062/2075=99.6%。
利用表1中的数据,计算出利用本发明所述方法对黄淮海平原的农田种植状态识别结果的总体精度为(1882+2062)/(1882+8+193+2062)=95.2%,Kappa系数为90.5%,从而可见,尽管本发明使用中低分辨率的遥感数据作为源数据,但是仍然获得了较高的识别精度。
实施例三
如图9所示,为本实施例中农田种植率的计算方法流程示意图。
步骤S1',获取遥感数据,其为免费的中低分辨率的遥感数据。
步骤S2',识别农田种植状态,具体如前面实施例一中的方法,根据所述方法逐像元、逐时期地获取农田种植状态。
步骤S3',统计特定时期、特定区域范围内休耕农田与有作物种植的农田像元数量。
步骤S4',计算有作物种植的农田像元数量占休耕农田与有作物种植的农田像元数量总和的比例,从而得到所述特定时期内的农田种植率。
对应于系统,如图10所示,除了包括与实施例二所述的监测系统相同的模块,如预处理模块1a、极值提取模块2a和数据处理模块3a,或预处理模块1b、极值提取模块2b、数据处理模块3b、植被指数阈值训练模块4b和验证模块5b外,还包括数据统计模块4c,用于统计特定区域内休耕农田像元数量、种植作物农田像元数量及计算农田种植状态。
与现有技术相比,本发明具有以下明显的技术优势:
1.充分利用时间序列NDVI数据的高时间频率特性,准确的获取了不同时期休耕农田与种植作物的农田反映在遥感影像上的显著差异,且方法流程标准,简单易行,具有一定的推广性;
2.采用了结合样本数据统计获取的植被指数累计直方分布与决策树的农田种植状态遥感分类方法,仅需要免费的中低分辨率的遥感数据即可实现高精度的农田种植状态与农田种植状态的动态监测。
以上实施方式仅用于说明本发明,而并非对本发明的限制,有关技术领域的普通技术人员,在不脱离本发明的精神和范围的情况下,还可以做出各种变化和变型,因此所有等同的技术方案也属于本发明的范畴,本发明的专利保护范围应由权利要求限定。

Claims (10)

1.一种农田种植状态的遥感监测方法,其中,包括:
获取遥感数据,经过预处理后得到像元的植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值;
提取像元植被指数曲线的极值点,所述极值点包括极大值点和极小值点;和
逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态。
2.如权利要求1所述的农田种植状态的遥感监测方法,其中,所述植被指数阈值包括休耕农田阈值和种植作物农田阈值,在所述时期的植被指数值大于或等于所述种植作物农田阈值时,确定所述时期的农田种植状态为有作物种植;
在所述时期的植被指数值小于所述休耕农田阈值时,确定所述时期的农田种植状态为休耕;
在所述时期的植被指数值小于所述种植作物农田阈值且大于或等于休耕农田阈值时,根据所述时期的植被指数数据的变化趋势,参考所述植被指数极大值与所述种植作物农田阈值的比较结果,确定所述时期的农田种植状态。
3.如权利要求2所述的农田种植状态的遥感监测方法,其中,参考所述植被指数极大值与所述种植作物农田阈值的比较结果,确定所述时期的农田种植状态具体包括:
判断所述时期的植被指数数据的变化趋势,如果所述时期的植被指数数据变化趋势为增大,获取所述时期之后最邻近的极大值;如果所述时期的植被指数数据变化趋势为减小,获取所述时期之前最邻近的极大值;
比较所述像元植被指数极大值与所述种植作物农田阈值的大小;
如果所述植被指数极大值大于或等于所述种植作物农田阈值,确定所述时期的农田种植状态为有作物种植;如果所述植被指数极大值小于种植作物农田阈值,确定所述时期的农田种植状态为休耕;
如果所述时期的植被指数数据变化趋势不变,确定所述时期的植被指数值为极大值或极小值,当所述时期的植被指数值为极大值时,确定所述时期的农田种植状态为休耕;
当所述时期的植被指数值为极小值时,在所述极小值所在时期之前和所在时期之后分别获取与其最相邻的极大值,分别对比前、后相邻的极大值与所述种植作物农田阈值的大小,如果前、后极大值中至少一个大于或等于所述种植作物农田阈值,确定所述时期的农田种植状态为有作物种植,如果所述前、后极大值均小于所述种植作物农田阈值,确定所述时期的农田种植状态为休耕。
4.如权利要求3所述的农田种植状态的遥感监测方法,其中,通过以下步骤判断所述时期的植被指数数据的变化趋势:
在所述时期对所述像元植被指数曲线求导,当所述时期的像元植被指数导数大于0时,所述时期的植被指数数据变化趋势为增大;
当所述时期的像元植被指数导数小于0时,所述时期的植被指数数据变化趋势为减小;
当所述时期的像元植被指数导数等于0时,所述时期的植被指数数据变化趋势不变。
5.如权利要求1-4任一所述的农田种植状态的遥感监测方法,其中,在逐像元、逐时期比较植被指数值与植被指数阈值之前还包括:
根据训练样本训练得到所述植被指数阈值,所述植被指数阈值包括休耕农田阈值和种植作物农田阈值。
6.如权利要求1-4任一所述的农田种植状态的遥感监测方法,其中,还包括:
统计特定时期、特定区域范围内休耕农田与有作物种植农田像元数量;
计算有作物种植农田像元数量占休耕农田与有作物种植农田像元数量总和的比例,确定所述特定时期内的农田种植率。
7.一种农田种植状态的遥感监测系统,其中,包括:
预处理模块,用于对获取的遥感数据进行预处理,获得植被指数曲线数据,所述植被指数曲线数据包括由多个特定时期顺序组成的时间序列和与所述特定时期对应的植被指数值;
极值提取模块,用于提取像元植被指数曲线的极值点,所述极值点包括极大值点和极小值点;和
数据处理模块,用于逐像元、逐时期比较植被指数值与预设的植被指数阈值,根据所述时期的植被指数值与植被指数阈值的比较结果确定所述像元、所述时期的农田种植状态。
8.如权利要求7所述的农田种植状态的遥感监测系统,其中,所述数据处理模块包括:
数据读取单元,用于读取植被指数阈值和指定时期的植被指数值,其中,所述植被指数阈值包括休耕农田阈值和种植作物农田阈值;
比较单元,用于比较指定时期的植被指数值和休耕农田阈值的大小,比较指定时期的植被指数值和种植作物农田阈值的大小,及比较植被指数极大值与种植作物农田阈值的大小;和
识别单元,根据所述比较单元的比较结果,确定所述指定时期的农田种植状态。
9.如权利要求8所述的农田种植状态的遥感监测系统,其中,所述数据处理模块包括:
导数计算单元,用于在指定时期对所述像元植被指数曲线求导获得所述指定时期的导数;
所述比较单元还包括比较所述导数与0的大小。
10.如权利要求7所述的农田种植状态的遥感监测系统,其中,还包括:
植被指数阈值训练模块,根据训练样本训练得到所述植被指数阈值,所述植被指数阈值包括休耕农田阈值和种植作物农田阈值;和/或
验证模块,用于根据验证样本验证数据处理模块识别得到的农田种植状态的精度。
CN201710010876.8A 2017-01-06 2017-01-06 农田种植状态的遥感监测方法和系统 Pending CN106845806A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710010876.8A CN106845806A (zh) 2017-01-06 2017-01-06 农田种植状态的遥感监测方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710010876.8A CN106845806A (zh) 2017-01-06 2017-01-06 农田种植状态的遥感监测方法和系统

Publications (1)

Publication Number Publication Date
CN106845806A true CN106845806A (zh) 2017-06-13

Family

ID=59118290

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710010876.8A Pending CN106845806A (zh) 2017-01-06 2017-01-06 农田种植状态的遥感监测方法和系统

Country Status (1)

Country Link
CN (1) CN106845806A (zh)

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107437262A (zh) * 2017-07-19 2017-12-05 中国科学院遥感与数字地球研究所 作物种植面积预警方法和系统
CN109471131A (zh) * 2018-11-15 2019-03-15 中科禾信遥感科技(苏州)有限公司 通过遥感卫星照片统计监测轮作休耕情况的方法和装置
CN110869744A (zh) * 2017-07-18 2020-03-06 索尼公司 信息处理设备、信息处理方法、程序和信息处理系统
CN110909679A (zh) * 2019-11-22 2020-03-24 中国气象科学研究院 冬小麦历史种植区的休耕轮作信息遥感识别方法及系统
CN112308753A (zh) * 2020-11-24 2021-02-02 中国科学院东北地理与农业生态研究所 基于生态参量的生态系统属性组分组成结构描述方法
CN112347978A (zh) * 2020-11-25 2021-02-09 中国科学院东北地理与农业生态研究所 基于遥感植被指数的生态系统属性组分组成结构描述方法
CN112364302A (zh) * 2020-11-11 2021-02-12 中国科学院东北地理与农业生态研究所 融合属性分级信息的生态系统属性组分组成结构描述方法
CN112507858A (zh) * 2020-12-03 2021-03-16 中国科学院东北地理与农业生态研究所 基于遥感植被指数的生态系统属性组分组成结构描述方法
CN112840348A (zh) * 2019-10-11 2021-05-25 安徽中科智能感知产业技术研究院有限责任公司 一种基于时序遥感数据和卷积神经网络的作物种植分布预测方法
US20210158527A1 (en) * 2019-11-22 2021-05-27 Farmers Edge Inc. Dynamic Area Thresholding for Automatic Crop Health Change Detection and Alerting System
CN113822360A (zh) * 2021-09-24 2021-12-21 中化现代农业有限公司 一种作物的物候期确定方法及装置
CN114067158A (zh) * 2021-11-17 2022-02-18 江苏天汇空间信息研究院有限公司 应用多源遥感数据的农田使用状态的监测系统及方法
CN114463647A (zh) * 2021-12-22 2022-05-10 广州极飞科技股份有限公司 一种作业方法、装置、设备及存储介质
CN114708490A (zh) * 2021-12-14 2022-07-05 深圳先进技术研究院 水稻种植提取及复种指数监测方法、系统、终端及存储介质
CN118521969A (zh) * 2024-07-25 2024-08-20 广东省科学院广州地理研究所 一种水稻退种风险的监测方法

Cited By (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110869744A (zh) * 2017-07-18 2020-03-06 索尼公司 信息处理设备、信息处理方法、程序和信息处理系统
CN107437262A (zh) * 2017-07-19 2017-12-05 中国科学院遥感与数字地球研究所 作物种植面积预警方法和系统
CN107437262B (zh) * 2017-07-19 2021-06-04 中国科学院遥感与数字地球研究所 作物种植面积预警方法和系统
CN109471131B (zh) * 2018-11-15 2022-05-24 中科禾信遥感科技(苏州)有限公司 通过遥感卫星照片统计监测轮作休耕情况的方法和装置
CN109471131A (zh) * 2018-11-15 2019-03-15 中科禾信遥感科技(苏州)有限公司 通过遥感卫星照片统计监测轮作休耕情况的方法和装置
CN112840348B (zh) * 2019-10-11 2024-05-03 安徽中科智能感知科技股份有限公司 一种基于时序遥感数据和卷积神经网络的作物种植分布预测方法
CN112840348A (zh) * 2019-10-11 2021-05-25 安徽中科智能感知产业技术研究院有限责任公司 一种基于时序遥感数据和卷积神经网络的作物种植分布预测方法
CN110909679A (zh) * 2019-11-22 2020-03-24 中国气象科学研究院 冬小麦历史种植区的休耕轮作信息遥感识别方法及系统
US20210158527A1 (en) * 2019-11-22 2021-05-27 Farmers Edge Inc. Dynamic Area Thresholding for Automatic Crop Health Change Detection and Alerting System
US11688069B2 (en) * 2019-11-22 2023-06-27 Farmers Edge Inc. Dynamic area thresholding for automatic crop health change detection and alerting system
CN112364302A (zh) * 2020-11-11 2021-02-12 中国科学院东北地理与农业生态研究所 融合属性分级信息的生态系统属性组分组成结构描述方法
CN112308753A (zh) * 2020-11-24 2021-02-02 中国科学院东北地理与农业生态研究所 基于生态参量的生态系统属性组分组成结构描述方法
CN112347978A (zh) * 2020-11-25 2021-02-09 中国科学院东北地理与农业生态研究所 基于遥感植被指数的生态系统属性组分组成结构描述方法
CN112347978B (zh) * 2020-11-25 2022-06-03 中国科学院东北地理与农业生态研究所 基于遥感植被指数的生态系统属性组分组成结构描述方法
CN112507858B (zh) * 2020-12-03 2023-03-24 中国科学院东北地理与农业生态研究所 基于遥感植被指数的生态系统属性组分组成结构描述方法
CN112507858A (zh) * 2020-12-03 2021-03-16 中国科学院东北地理与农业生态研究所 基于遥感植被指数的生态系统属性组分组成结构描述方法
CN113822360A (zh) * 2021-09-24 2021-12-21 中化现代农业有限公司 一种作物的物候期确定方法及装置
CN114067158A (zh) * 2021-11-17 2022-02-18 江苏天汇空间信息研究院有限公司 应用多源遥感数据的农田使用状态的监测系统及方法
CN114708490A (zh) * 2021-12-14 2022-07-05 深圳先进技术研究院 水稻种植提取及复种指数监测方法、系统、终端及存储介质
CN114708490B (zh) * 2021-12-14 2024-08-16 深圳先进技术研究院 水稻种植提取及复种指数监测方法、系统、终端及存储介质
CN114463647A (zh) * 2021-12-22 2022-05-10 广州极飞科技股份有限公司 一种作业方法、装置、设备及存储介质
CN118521969A (zh) * 2024-07-25 2024-08-20 广东省科学院广州地理研究所 一种水稻退种风险的监测方法

Similar Documents

Publication Publication Date Title
CN106845806A (zh) 农田种植状态的遥感监测方法和系统
Boschetti et al. PhenoRice: A method for automatic extraction of spatio-temporal information on rice crops using satellite data time series
CN105760978B (zh) 一种基于温度植被干旱指数(tvdi)的农业旱灾等级监测方法
Yu et al. Automatic image-based detection technology for two critical growth stages of maize: Emergence and three-leaf stage
CN110472184A (zh) 一种基于Landsat遥感数据的多云雨雾地区水稻识别方法
CN102194127B (zh) 一种多频率sar数据农作物遥感分类方法
CN106372592B (zh) 一种基于冬小麦面积指数的冬小麦种植面积计算方法
CN109389049A (zh) 基于多时相sar数据与多光谱数据的作物遥感分类方法
CN109614891A (zh) 基于物候学和遥感的农作物识别方法
CN107014753A (zh) 作物长势监测方法和系统
Wu et al. Characterizing spatial patterns of phenology in cropland of China based on remotely sensed data
CN111209871B (zh) 一种基于光学卫星影像的油菜种植地块遥感自动识别方法
CN104951772B (zh) 一种基于ndvi时间序列曲线积分的冬小麦提取方法
Pazhanivelan et al. Rice crop monitoring and yield estimation through COSMO Skymed and TerraSAR-X: A SAR-based experience in India
CN105372672A (zh) 基于时间序列数据的南方冬种作物种植面积提取方法
Liu et al. Optimal MODIS data processing for accurate multi-year paddy rice area mapping in China
CN108805079A (zh) 冬小麦的识别方法及装置
CN103971199A (zh) 一种大范围农作物长势的遥感评级方法
CN110298322A (zh) 一种基于遥感数据的耕地提取方法及系统
CN107437262B (zh) 作物种植面积预警方法和系统
Dalposso et al. Spatial autocorrelation of NDVI and GVI indices derived from Landsat/TM images for soybean crops in the western of the state of Paraná in 2004/2005 crop season
Raza et al. Comparative geospatial approach for agricultural crops identification in interfluvial plain-a case study of Sahiwal district, Pakistan
CN112001809A (zh) 一种农林区退耕地信息获取方法
CN108280410A (zh) 一种基于二进制编码的农作物识别方法及系统
CN114219227A (zh) 一种精准施肥决策及地块级展示方法和系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20170613