CN108984803B - 一种农作物产量空间化的方法及*** - Google Patents
一种农作物产量空间化的方法及*** Download PDFInfo
- Publication number
- CN108984803B CN108984803B CN201811226264.3A CN201811226264A CN108984803B CN 108984803 B CN108984803 B CN 108984803B CN 201811226264 A CN201811226264 A CN 201811226264A CN 108984803 B CN108984803 B CN 108984803B
- Authority
- CN
- China
- Prior art keywords
- ndvi
- yield
- crop
- statistical
- area
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 42
- 238000012417 linear regression Methods 0.000 claims abstract description 38
- 238000000638 solvent extraction Methods 0.000 claims abstract description 21
- 238000004458 analytical method Methods 0.000 claims abstract description 15
- 238000007781 pre-processing Methods 0.000 claims description 16
- 238000003066 decision tree Methods 0.000 claims description 10
- 238000001914 filtration Methods 0.000 claims description 9
- 238000010606 normalization Methods 0.000 claims description 9
- 238000005192 partition Methods 0.000 claims description 8
- 238000011160 research Methods 0.000 claims description 8
- 238000012952 Resampling Methods 0.000 claims description 7
- 238000000605 extraction Methods 0.000 claims description 7
- 238000009499 grossing Methods 0.000 claims description 7
- 230000001419 dependent effect Effects 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 5
- 238000012216 screening Methods 0.000 claims description 5
- 238000005457 optimization Methods 0.000 abstract description 6
- 241000209140 Triticum Species 0.000 description 81
- 235000021307 Triticum Nutrition 0.000 description 81
- 238000009826 distribution Methods 0.000 description 16
- 238000012795 verification Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 235000013339 cereals Nutrition 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 239000008267 milk Substances 0.000 description 3
- 210000004080 milk Anatomy 0.000 description 3
- 235000013336 milk Nutrition 0.000 description 3
- 239000002689 soil Substances 0.000 description 3
- 238000011497 Univariate linear regression Methods 0.000 description 2
- 238000012271 agricultural production Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 235000017060 Arachis glabrata Nutrition 0.000 description 1
- 244000105624 Arachis hypogaea Species 0.000 description 1
- 235000010777 Arachis hypogaea Nutrition 0.000 description 1
- 235000018262 Arachis monticola Nutrition 0.000 description 1
- 229920000742 Cotton Polymers 0.000 description 1
- 244000068988 Glycine max Species 0.000 description 1
- 235000010469 Glycine max Nutrition 0.000 description 1
- 241000219146 Gossypium Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 240000008042 Zea mays Species 0.000 description 1
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 description 1
- 235000002017 Zea mays subsp mays Nutrition 0.000 description 1
- 230000009418 agronomic effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 235000005822 corn Nutrition 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 235000020232 peanut Nutrition 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 238000013316 zoning Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION 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/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Business, Economics & Management (AREA)
- Computational Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Evolutionary Biology (AREA)
- Tourism & Hospitality (AREA)
- Primary Health Care (AREA)
- General Business, Economics & Management (AREA)
- Marketing (AREA)
- Human Resources & Organizations (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Economics (AREA)
- Health & Medical Sciences (AREA)
- Strategic Management (AREA)
- Mining & Mineral Resources (AREA)
- Marine Sciences & Fisheries (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Animal Husbandry (AREA)
- Algebra (AREA)
- Agronomy & Crop Science (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种农作物产量空间化的方法及***。该方法包括:通过遥感影像获取随时间变化的NDVI;采用CART分类回归树算法确定种植范围;对种植范围进行分区;获取各区域农作物的统计产量;确定农作物的生长季,生长季包括农作物的多个单时相以及由单时相组合而成的不同物候期;根据生长季各单时相和物候期的NDVI,采用线性回归分析法分析农作物的统计产量与生长季的NDVI的相关性,确定与农作物的统计产量相关性最高、均方根误差最小的生长季的NDVI变量,记为目标生长季NDVI变量;根据目标生长季中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量。本发明能够反映出农作物产量的空间格局和动态变化,为农作物种植结构的优化提供参考。
Description
技术领域
本发明涉及农作物产量统计领域,特别是涉及一种农作物产量空间化的方法及***。
背景技术
社会统计数据是反映一个国家或地区社会发展状况的重要指标之一,是国家或地区制定法规政策和宏观控制管理的重要依据。传统的社会统计数据已经不能满足人们的需求。传统社会统计数据一般是以行政区划(省级、地市级、县级、乡镇等)为统计单元收集、整理统计资料,具有空间分辨率低,缺乏对数据空间分布特征的描述,数据不完整等缺点,不能详细的表达行政区划内的空间差异,造成社会统计数据与其他数据的叠置分析比较困难,影响了有效信息的采集,严重阻碍了有效信息的提取。社会统计数据的空间化是解决这一问题的有效方法。
20世纪90年代初期,学者们提出了“空间化”的概念,空间化是利用一定的方法或参数构建模型,然后利用模型反演出数据在一定时间和空间尺度上的分布情况,其实质就是把数据反演到公里格网上。社会经济统计数据空间化成为诸多学科中的热点,目前人们对于社会统计数据空间化研究方法主要分为空间插值模型、土地利用/土地覆盖影响模型、多源数据融合模型和遥感反演模型。
目前社会统计数据空间化的内容主要集中在人口统计数据空间化;GDP统计数据空间化和其他属性数据空间化等。关于农业经济统计空间化的研究很少。农业生产活动受地理位置、气候特征、土壤水文等自然条件的影响,对农业统计数据的空间化难度较大。在农业经济统计数据的空间化中,研究者主要集中在农作物种植面积的空间化和农业生产投入方面的空间化方面。对粮食产量统计数据的空间化研究更少。
发明内容
本发明的目的是提供一种农作物产量空间化的方法及***,能够反映出农作物产量的空间格局和动态变化,为农作物种植结构的优化提供参考。
为实现上述目的,本发明提供了如下方案:
一种农作物产量空间化的方法,所述方法包括:
通过遥感影像获取随时间变化的NDVI;
根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围;
对所述种植范围进行分区;
获取各区域的所述农作物的统计产量;
确定所述农作物的生长季,所述生长季包括所述农作物的多个单时相以及由所述单时相组合而成的多个物候期;
根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与各所述生长季的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,记为目标生长季NDVI变量;
根据所述目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量。
可选的,所述对所述种植范围进行分区,具体包括:
结合地势地貌、DEM数字高程数据和行政单元矢量数据对所述种植范围进行分区。
可选的,所述根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与所述生长季NDVI变量的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,具体包括:
计算各区域在各单时相的NDVI;
根据各区域在各单时相的NDVI计算各区域在各物候期的NDVI;
以各单时相的NDVI以及各物候期的NDVI为自变量,以农作物的统计产量为因变量构建线性回归方程;
选取与所述农作物的统计产量相关性最高且拟合方程的均方根误差最小的回归方程,记为目标回归方程,将所述目标回归方程中自变量所代表的生长季的NDVI记为目标生长季NDVI变量。
可选的,所述根据所述目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的产量,具体包括:
根据计算各像元的农作物空间化产量,其中,Y′j表示第j(j=1,...n)个像元的所述农作物的空间化产量,n表示待空间化的区域中的像元个数,Y表示待空间化的区域的统计产量,i表示目标生长季,k表示目标生长季的个数,NDVIi,j表示第i个目标生长季第j个像元的NDVI,表示待空间化的区域中第i个目标生长季的NDVI总和,pi表示第i个目标生长季的NDVI与统计产量之间的相关系数,αi表示第i个时期的相关系数归一化结果。
可选的,所述根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围之前,还包括:
对所述NDVI进行预处理,所述预处理包括分辨率重采样、滤波去噪和平滑处理。
本发明还提供了一种农作物产量空间化的***,所述***包括:
数据获取模块,用于通过遥感影像获取随时间变化的NDVI;
种植范围确定模块,用于根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围;
分区模块,用于对所述种植范围进行分区;
统计产量获取模块,用于获取各区域的所述农作物的统计产量;
生长季确定模块,用于确定所述农作物的生长季,所述生长季包括所述农作物的多个单时相以及由所述单时相组合而成的多个物候期;
目标生长季确定模块,用于根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与各所述生长季的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,记为目标生长季NDVI变量;
产量空间化模块,用于根据所述目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量。
可选的,所述分区模块,具体包括:
分区单元,用于结合地势地貌、DEM数字高程数据和行政单元矢量数据对所述种植范围进行分区。
可选的,所述目标生长季确定模块,具体包括:
第一计算单元,用于计算各区域在各单时相的NDVI;
第二计算单元,用于根据各区域在各单时相的NDVI计算各区域在各物候期的NDVI;
线性回归方程构建单元,用于以各单时相的NDVI以及各物候期的NDVI为自变量,以农作物的统计产量为因变量构建线性回归方程;
目标生长季确定单元,用于选取与所述农作物的统计产量相关性最高且拟合方程的均方根误差最小的回归方程,记为目标回归方程,将所述目标回归方程中自变量所代表的生长季的NDVI记为目标生长季NDVI变量。
可选的,所述产量空间化模块,具体包括:
产量空间化单元,用于根据计算各像元的农作物空间化产量,其中,Y′j表示第j(j=1,...n)个像元的所述农作物的空间化产量,n表示待空间化的区域中的像元个数,Y表示待空间化的区域的统计产量,i表示目标生长季,k表示目标生长季的个数,NDVIi,j表示第i个目标生长季第j个像元的NDVI,表示待空间化的区域中第i个目标生长季的NDVI总和,pi表示第i个目标生长季的NDVI与统计产量之间的相关系数,αi表示第i个时期的相关系数归一化结果。
可选的,所述***还包括:
预处理模块,用于对所述NDVI进行预处理,所述预处理包括分辨率重采样、滤波去噪和平滑处理。
根据本发明提供的具体实施例,本发明公开了以下技术效果:本发明提供的农作物产量空间化方法及***,基于农作物生长季的NDVI数据,采用CART决策树分类方法得到农作物种植面积分布图,提取农作物生长季的不同物候期的NDVI均值和各单时相的NDVI值,通过线性回归分析筛选与农作物统计产量相关性最高、拟合方程的均方根误差最小的NDVI变量,建立农作物产量空间化模型,实现农作物统计产量的空间化。农作物产量统计数据的空间化可以反映出地区内部的粮食产量空间格局和动态变化,为农作物种植结构的优化提供参考。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例农作物产量空间化方法的流程示意图;
图2为本发明实施例研究区的位置图;
图3为本发明实施例山东省冬小麦空间分布图;
图4为本发明实施例回归标准化残差直方图;
图5为本发明实施例山东省冬小麦县级产量空间分布图;
图6为本发明实施例山东省冬小麦市级产量空间分布图;
图7为本发明实施例山东省冬小麦产量空间化结果对比;
图8为本发明实施例农作物产量空间化***的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的目的是提供一种农作物产量空间化的方法及***,能够反映出农作物产量的空间格局和动态变化,为农作物种植结构的优化提供参考。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明实施例农作物产量空间化方法的流程示意图,如图1所示,本发明提供的农作物产量空间化方法包括以下步骤:
步骤101:通过遥感影像获取随时间变化的NDVI;
步骤102:根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围;
步骤103:对种植范围进行分区;
步骤104:获取各区域的农作物的统计产量;
步骤105:确定所述农作物的生长季,所述生长季包括所述农作物的多个单时相以及由所述单时相组合而成的多个物候期;
步骤106:根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与各所述生长季的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,记为目标生长季NDVI变量;
步骤107:根据所述目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量。
其中,步骤101具体为:从美国地质勘探局USGS的LP DAAC(land processdistributed active archive center)下载待空间化年份的250m分辨率的MOD13Q1植被指数时间序列数据。
在步骤102之前,对植被指数时间序列数据进行预处理,具体为:利用MODIS投影工具(MRT)软件将MOD13Q1时间序列数据转为阿尔伯斯(Albers)等面积投影,影像的分辨率重采样为250m,利用Savitzky-Golay滤波方法对MODIS NDVI时间序列数据进行滤波去噪、平滑处理,提取农作物生长季的MODIS NDVI数据。
步骤102具体包括:
基于分类回归树算法的农作物种植范围提取,具体包括结合高分辨率遥感数据的农作物样本选取、基于CART算法的农作物识别的决策树建立和农作物识别的精度验证。本专利所记载的农作物指的是某一种农作物,例如,该农作物为冬小麦,并不是多种农作物的总和。
步骤103具体包括:
结合地势地貌、DEM数字高程数据和行政单元矢量数据对种植范围进行分区。分区的原则包括:1)行政单元的完整性;2)不同行政单元间高程差异大;3)不同行政单元间冬小麦地块的破碎程度不同。
步骤104具体为:获取各区域的农作物的统计产量,例如,整理冬小麦种植面积和产量统计数据,检查数据的完整性和规范性,利用ArcGIS软件将统计数据和矢量边界数据进行空间链接,确保每个分区都有冬小麦种植面积和产量统计数据。分区可以以行政单元县进行划分。
步骤106具体包括:
计算各区域在各单时相的NDVI;
根据各区域在各单时相的NDVI计算各区域在各物候期的NDVI;
式中,NDVI′i表示农作物第i个单时相的NDVI均值,m表示待空间化的分区(如县)中的农作物像元个数,NDVIi表示第i个单时相的NDVI,NDVI″表示农作物物候期内第j个物候期的NDVI,n表示第j个物候期内所包含的农作物单时相个数。该农作物可以为冬小麦,单时相的NDVI均值可以为:3月6日、3月22日、4月7日、4月23日、5月9日、5月25日、6月10日和6月26日的冬小麦NDVI的均值,物候期的NDVI可以为返青期、拔节期、抽穗期和乳熟期的冬小麦的NDVI均值。
以各单时相的NDVI以及各物候期的NDVI为自变量,以农作物的统计产量为因变量构建多元线性回归方程;例如,首先,建立各个因子(单时相的NDVI均值和物候期的NDVI均值)与冬小麦统计产量之间的一元回归方程,确定相关系数并进行显著性检验;然后,采用自变量进入的方法建立单时相的NDVI、物候期的NDVI与冬小麦统计产量之间的多元回归方程,确定输入因子的个数和相关系数;最后,对比分析各因子与冬小麦产量之间的相关系数、回归标准化残差直方图以及均方根误差(RMSE),选取与冬小麦产量相关性最高、拟合方程的RMSE最小的回归方程的NDVI变量作为冬小麦产量空间化的因子。一元回归方程和多元回归方程公式如下:
Y=aXi+b
Y=β0+β1X1+...+βnXn
式中,Y表示冬小麦的统计产量值,Xi表示不同时相的NDVI,a,b和β0…βn均为常数。
选取与所述农作物的统计产量相关性最高、拟合方程的均方根误差最小的回归方程,记为目标回归方程,将所述目标回归方程中自变量所代表的生长季的NDVI记为目标生长季NDVI变量。
步骤107具体包括:
根据计算各像元的农作物空间化产量,其中,Y′j表示第j(j=1,...n)个像元的所述农作物的空间化产量,n表示待空间化的区域中的像元个数,Y表示待空间化的区域的统计产量,i表示目标生长季,k表示目标生长季的个数,NDVIi,j表示第i个目标生长季第j个像元的NDVI,表示待空间化的区域中第i个目标生长季的NDVI总和,pi表示第i个目标生长季的NDVI与统计产量之间的相关系数,αi表示第i个时期的相关系数归一化结果。
下面举例说明:
本案例选取山东省作为研究区(如图2所示)。山东省位于中国东部沿海、黄河下游,北纬34°22.9′-38°24.01′、东经114°47.5′-122°42.3′之间,陆地面积15.71万km2,17个地级市,县级单位137个。山东省的气候属暖温带季风气候类型,降水集中,雨热同期,全年无霜期由东北沿海向西南递增,光照资源充足,热量条件可满足农作物一年两作的需要。山东省是全国的产量大省之一,主要农作物有小麦、玉米、水稻、大豆、棉花和花生等。
本案例所用数据包括:1)来源于美国地质勘探局USGS的LP DAAC(land processdistributed active archive center)的2016年250m分辨率的MOD13Q1时间序列数据,MOD13Q1为植被指数产品,用于提取冬小麦的种植面积和产量空间化;2)来源于美国地质勘探局(USGS)(http://glovis.usgs.gov/)90m分辨率的DEM数字高程数据,辅助研究区域的分区;3)来源于山东省***(http://xxgk.stats-sd.gov.cn/)的2016年冬小麦种植面积和产量数据,用于冬小麦面积提取的精度验证和冬小麦产量的空间化;4)来源于地理空间数据云(http://www.gscloud.cn/)的2016年小麦生长季内的LandsatTM、GF-1、GF-2等高分辨率遥感数据,用于冬小麦样本的选取。
本案例主要流程包括:(1)研究区MODIS NDVI时间序列数据的获取和预处理;(2)冬小麦种植面积和产量统计数据的整理;(3)基于分区和分类回归树算法的冬小麦种植范围提取;(4)冬小麦主要生长季NDVI变量计算;(5)分析NDVI变量与产量间的关系,筛选最优空间化因子;(6)冬小麦产量的空间化模型构建与应用。
本案例研究区数据的预处理包括:首先,利用MRT软件将2016年的MOD13Q1时间序列数据转为Albers等面积投影,分辨率重采样为250m。由于MOD13Q1产品本身存在噪声和云污染现象,利用Savitzky-Golay滤波方法对2016年的MODIS-NDVI时间序列数据进行滤波去噪、平滑处理,提取冬小麦生长季的MODIS-NDVI数据作为待分类数据。
本案例冬小麦种植面积和产量统计数据的整理是检查数据的完整性和规范性,将统计数据和矢量边界数据进行空间链接,确保每个地区都有冬小麦种植面积和产量统计数据。
本案例在基于分区和分类回归树算法的冬小麦种植范围提取中考虑到山东省的地势地貌以及冬小麦的种植地块复杂程度,结合DEM数字高程数据和行政单元矢量数据将山东省按行政单元(县)划分为平原-简单地物区域和平原/丘陵-复杂地物区域。通过高分辨率遥感数据(LandsatTM、GF-1、GF-2等)分区域选取冬小麦样本点(图2),基于冬小麦生长季的MODIS-NDVI数据,采用CART分类方法自动建立冬小麦种植面积识别的决策树,分别提取平原-简单地物区域和平原/丘陵-复杂地物区域的冬小麦种植面积,利用县级矢量数据统计提取冬小麦种植面积,与冬小麦种植面积的统计数据进行比较,计算每个区县的种植面积提取精度,筛选出分类精度低于60%的区县,针对筛选出区县的影像重新建立决策树进行分类,以提高分类精度。最终得到的山东省的冬小麦种植分布如图3所示。同时通过生成随机点的方法在山东省选取了595个冬小麦样点(图2),其中平原-简单地物区域冬小麦样本点为403,平原/丘陵-复杂地物区域冬小麦样本点为192,采用混淆矩阵方法对提取的种植面积进行精度验证,验证结果如表1所示。从表1可以看出简单地物区域分类结果总体精度高于复杂地物区域的分类精度,简单地物区域分类结果中的漏分误差高于复杂地物区域的分类结果,复杂地物区域的分类结果中的错分误差高于简单地物区域分类结果;山东省冬小麦种植面积识别的总体精度为82.51%,制图精度为83.26%,用户精度为78.91%。另外,通过统计提取的冬小麦种植面积与实际冬小麦种植面积相比,发现简单地物区域的冬小麦种植面积识别的平均精度为92.88%,复杂地物区域的冬小麦种植面积识别的平均精度为81.28%,山东省冬小麦种植面积识别的平均精度为87.64%。两种精度验证表明简单地物区域分类结果优于复杂地物区域的分类结果。
表1冬小麦种植面积精度验证
本案例在冬小麦主要生长季NDVI变量计算中提取了2016年山东省的冬小麦主要生长季(3月-6月)的合计8期的16天合成的NDVI数据(每期数据的起始时间分别为:3月6日、3月22日、4月7日、4月23日、5月9日、5月25日、6月10日和6月26日),用NDVImd表示,m表示月份,d表示天数(如:NDVI0509表示5月9日至5月24日期间的16天合成的NDVI)。结合山东省的冬小麦物候信息,将冬小麦生长季的8个单时相数据分配到各个物候期,以县级为单元统计冬小麦每个物候期内单时相的NDVI均值(式1),以及返青期(NDVIrgs,为NDVI0306和NDVI0322的均值)、拔节期(NDVIjs,为NDVI0407和NDVI0423的均值)、抽穗期(NDVIhs,为NDVI0509和NDVI0525的均值)和乳熟期(NDVImrs,为NDVI0610和NDVI0626的均值)的NDVI(式2)。
式中,NDVI′i表示冬小麦第i个单时相的NDVI均值,m表示县域中的冬小麦像元个数,NDVIi表示第i个单时相的NDVI,NDVI表示冬小麦物候期内第j个物候期的NDVI,n表示第j个物候期内所包含的冬小麦单时相个数。
本案例分析NDVI变量与产量间的关系,筛选最优空间化因子的具体做法是在SPSS软件中,建立各个因子与冬小麦统计产量之间的一元回归方程,确定相关系数并进行显著性检验;然后采用自变量进入的方法建立单时相的NDVI、物候期的NDVI与冬小麦统计产量之间的多元回归方程,确定输入因子的个数和相关系数。一元回归方程和多元回归方程表示方法如下:
式3:Y=aXi+b
式4:Y=β0+β1X1+...+βnXn
式中,Y表示冬小麦的统计产量值,Xi表示不同时相的NDVI,a,b和β0…βn均为常数。最后对比分析各因子与冬小麦产量之间的相关系数、回归标准化残差直方图以及拟合方程的RMSE,选取与冬小麦产量相关性最高、RMSE最小的回归方程的NDVI变量作为冬小麦产量空间化的因子。通过一元线性回归和多元回归分析得到不同NDVI变量和冬小麦产量之间的相关关系(表2)。由表2可知,在一元线性回归分析中,当自变量为单时相的NDVI时,NDVI0407与冬小麦产量的相关性最高,R2为0.901,P≤0.05,通过显著性检验;当自变量为物候期的NDVI时,拔节期的NDVI(NDVIjs)与冬小麦产量的相关性最高,R2为0.900,P≤0.05,通过显著性检验。
在多元线性回归分析中,自变量输入方法为“Enter”,当自变量输入为单时相的NDVI变量时,筛选出参与线性回归的因子为:NDVI0306、NDVI0423、NDVI0525、NDVI0626;当自变量输入为物候期的NDVI变量时,筛选出参与线性回归的因子为:NDVIrgs、NDVIjs、NDVIhs、NDVImrs,所有物候期NDVI都参与线性回归模型;当输入自变量为单时相NDVI和物候期的NDVI变量时,筛选出参与线性回归的因子为:NDVI0306、NDVI0423、NDVI0525、NDVImrs。3种不同自变量的多元线性回归中,R2均为0.903,P≤0.05,均通过显著性检验。
表2线性回归相关性
注:NDVIrgs表示返青期的NDVI,NDVIjs表示拔节期的NDVI,NDVIhs表示抽穗期的NDVI,NDVImrs表示乳熟期的NDVI。*表示P≤0.05,**表示P≤0.01。
为了更进一步筛选参与冬小麦产量空间化的因子,分别对每个一元线性回归和多元线性回归的残差进行分析,得到回归标准化残差直方图(图4)。在回归标准化残差直方图中,正态曲线是判断回归标准化残差直方图是否符合正态分布的标准,由图4可知,在一元线性回归中,所有线性回归方程的回归标准化残差直方图都符合正态分布趋势,其中多元线性回归的回归标准化残差直方图的标准差比较小,更符合正态分布趋势,表明多元线性回归模型优于一元线性回归模型。
通过比较3个多元线性回归的拟合效果,分别计算拟合方程的RMSE,发现自变量为单时相的NDVI变量时,RMSE最小。单时相的NDVI变量为:NDVI0306、NDVI0423、NDVI0525、NDVI0626,因此确定参与冬小麦产量空间化为多元线性回归模型中的NDVI0306、NDVI0423、NDVI0525、NDVI06264个单时相NDVI变量。
本案例中空间化模型的构建是在确定回归因子的基础上,采用权重分配(面积插值法)方法,按县级行政单元计算每个像元中各时期NDVI占县级行政单元中NDVI总量的权重,将县级行政单元总产量按照因子权重分配到每个像元,得到山东省冬小麦产量空间分布图,具体空间化模型如式5和式6所示,生成的县级产量空间分布图如图5所示。由图5可知,山东省西部地区和西南部分地区的冬小麦产量较高,中部丘陵地区和沿海边缘地区冬小麦产量较低,产量分布呈现西高东低的现象。其原因可能是山东省西部地区和西南地区主要以平原为主,雨水充足,耕种条件良好,适应于冬小麦的种植;中部丘陵地区地形复杂,耕种条件差,土壤湿度较低,相应冬小麦产量较低;沿海边缘地区主要以第二三产业为主,工业的快速发展迫使土壤质量的下降以及农民对冬小麦种植的重视程度,导致冬小麦产量较低。
式中,Y′j表示第j(j=1,...n)个像元的冬小麦空间化产量,n表示县域中冬小麦的像元个数,Y表示县级冬小麦统计产量,i(i=1,2,3,4)表示3月6日、4月23日、5月25日和6月26日4个单时相时期,NDVIi,j表示第i个时期第j个像元的NDVI变量,表示县级单元中第i个时期的NDVI总和,pi表示第i个时期的NDVI与统计产量之间的相关系数,αi表示第i个时期的相关系数归一化结果。
本发明提供的农作物产量空间化方法,基于农作物生长季的NDVI(MODIS NDVI)数据,采用CART决策树分类方法得到农作物种植面积分布图,提取农作物生长季的不同物候期的NDVI均值和单时相的NDVI值,通过线性回归分析筛选与农作物统计产量相关性最高、拟合方程的均方根误差最小的NDVI变量,建立农作物产量空间化模型,实现农作物统计产量的空间化。农作物产量统计数据的空间化可以反映出地区内部的粮食产量空间格局和动态变化,为农作物种植结构的优化提供参考。
本发明一种农作物产量空间化的***,图8为本发明实施例农作物产量空间化***的结构示意图,如图8所示,本发明提供的***包括:
数据获取模块801,用于通过遥感影像获取随时间变化的NDVI;
种植范围确定模块802,用于根据NDVI,采用CART分类回归树算法确定农作物的种植范围;
分区模块803,用于对种植范围进行分区;
统计产量获取模块804,用于获取各区域的农作物的统计产量;
生长季确定模块805,用于确定农作物的生长季,生长季包括农作物的多个单时相以及由单时相组合而成的多个物候期;
目标生长季确定模块806,用于根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析农作物的统计产量与各生长季的相关性,确定与农作物的统计产量相关性最高且均方根误差最小的NDVI变量,记为目标生长季NDVI变量;
产量空间化模块807,用于根据目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量。
其中,分区模块803,具体包括:
分区单元,用于结合地势地貌、DEM数字高程数据和行政单元矢量数据对种植范围进行分区。
目标生长季确定模块806,具体包括:
第一计算单元,用于计算各区域在各单时相的NDVI;
第二计算单元,用于根据各区域在各单时相的NDVI计算各区域在各物候期的NDVI;
线性回归方程构建单元,用于以各单时相的NDVI以及各物候期的NDVI为自变量,以农作物的统计产量为因变量构建线性回归方程;
目标生长季确定单元,用于选取与农作物的统计产量相关性最高且拟合方程的均方根误差最小的回归方程,记为目标回归方程,将目标回归方程中自变量所代表的生长季的NDVI记为目标生长季NDVI变量。
产量空间化模块807,具体包括:
产量空间化单元,用于根据计算各像元的农作物空间化产量,其中,Y′j表示第j(j=1,...n)个像元的农作物的空间化产量,n表示待空间化的区域中的像元个数,Y表示待空间化的区域的统计产量,i表示目标生长季,k表示目标生长季的个数,NDVIi,j表示第i个目标生长季第j个像元的NDVI,表示待空间化的区域中第i个目标生长季的NDVI总和,pi表示第i个目标生长季的NDVI与统计产量之间的相关系数,αi表示第i个时期的相关系数归一化结果。
本发明提供的农作物产量空间化***还包括:
预处理模块,用于对NDVI进行预处理,预处理包括分辨率重采样、滤波去噪和平滑处理。
本发明提供的农作物产量空间化***,基于农作物生长季的NDVI(MODIS NDVI)数据,采用CART决策树分类方法得到农作物种植面积分布图,提取农作物生长季的不同物候期的NDVI均值和单时相的NDVI值,通过线性回归分析筛选与农作物统计产量相关性最高、拟合方程的均方根误差最小的NDVI变量,建立农作物产量空间化模型,实现农作物统计产量的空间化。农作物产量统计数据的空间化可以反映出地区内部的粮食产量空间格局和动态变化,为农作物种植结构的优化提供参考。
本说明书中各个实施例采用递进的方式描述,每个实施例重点说明的都是与其他实施例的不同之处,各个实施例之间相同相似部分互相参见即可。对于实施例公开的***而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处。综上所述,本说明书内容不应理解为对本发明的限制。
Claims (10)
1.一种农作物产量空间化的方法,其特征在于,所述方法包括:
通过遥感影像获取随时间变化的NDVI;
根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围;
对所述种植范围进行分区;
获取各区域的所述农作物的统计产量;
确定所述农作物的生长季,所述生长季包括所述农作物的多个单时相以及由所述单时相组合而成的多个物候期;
根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与各所述生长季的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,记为目标生长季NDVI变量;
根据所述目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量;
其中,所述根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围,具体包括:
结合DEM数字高程数据和行政单元矢量数据将研究区域进行分区,选取各区域的农作物样本点,基于农作物生长季的NDVI数据,采用CART分类方法建立农作物种植面积识别的决策树,分别提取各区域的农作物种植面积,并与农作物种植面积的统计数据进行比较,计算各区域的种植面积提取精度,筛选出分类精度低于设定阈值的区域,对筛选出的区域重新建立决策树进行分类。
2.根据权利要求1所述的农作物产量空间化的方法,其特征在于,所述对所述种植范围进行分区,具体包括:
结合地势地貌、DEM数字高程数据和行政单元矢量数据对所述种植范围进行分区。
3.根据权利要求1所述的农作物产量空间化的方法,其特征在于,所述根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与所述生长季NDVI变量的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,具体包括:
计算各区域在各单时相的NDVI;
根据各区域在各单时相的NDVI计算各区域在各物候期的NDVI;
以各单时相的NDVI以及各物候期的NDVI为自变量,以农作物的统计产量为因变量构建线性回归方程;
选取与所述农作物的统计产量相关性最高且拟合方程的均方根误差最小的回归方程,记为目标回归方程,将所述目标回归方程中自变量所代表的生长季的NDVI记为目标生长季NDVI变量。
5.根据权利要求1所述的农作物产量空间化的方法,其特征在于,所述根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围之前,还包括:
对所述NDVI进行预处理,所述预处理包括分辨率重采样、滤波去噪和平滑处理。
6.一种农作物产量空间化的***,其特征在于,所述***包括:
数据获取模块,用于通过遥感影像获取随时间变化的NDVI;
种植范围确定模块,用于根据所述NDVI,采用CART分类回归树算法确定所述农作物的种植范围,具体为:结合DEM数字高程数据和行政单元矢量数据将研究区域进行分区,选取各区域的农作物样本点,基于农作物生长季的NDVI数据,采用CART分类方法建立农作物种植面积识别的决策树,分别提取各区域的农作物种植面积,并与农作物种植面积的统计数据进行比较,计算各区域的种植面积提取精度,筛选出分类精度低于设定阈值的区域,对筛选出的区域重新建立决策树进行分类;
分区模块,用于对所述种植范围进行分区;
统计产量获取模块,用于获取各区域的所述农作物的统计产量;
生长季确定模块,用于确定所述农作物的生长季,所述生长季包括所述农作物的多个单时相以及由所述单时相组合而成的多个物候期;
目标生长季确定模块,用于根据生长季各单时相和物候期的NDVI,在各区域内,采用线性回归分析法分析所述农作物的统计产量与各所述生长季的相关性,确定与所述农作物的统计产量相关性最高且均方根误差最小的NDVI变量,记为目标生长季NDVI变量;
产量空间化模块,用于根据所述目标生长季NDVI变量中各像元的NDVI占所属区域NDVI总和的比重,确定各像元的农作物空间化产量。
7.根据权利要求6所述的农作物产量空间化的***,其特征在于,所述分区模块,具体包括:
分区单元,用于结合地势地貌、DEM数字高程数据和行政单元矢量数据对所述种植范围进行分区。
8.根据权利要求6所述的农作物产量空间化的***,其特征在于,所述目标生长季确定模块,具体包括:
第一计算单元,用于计算各区域在各单时相的NDVI;
第二计算单元,用于根据各区域在各单时相的NDVI计算各区域在各物候期的NDVI;
线性回归方程构建单元,用于以各单时相的NDVI以及各物候期的NDVI为自变量,以农作物的统计产量为因变量构建线性回归方程;
目标生长季确定单元,用于选取与所述农作物的统计产量相关性最高且拟合方程的均方根误差最小的回归方程,记为目标回归方程,将所述目标回归方程中自变量所代表的生长季的NDVI记为目标生长季NDVI变量。
10.根据权利要求6所述的农作物产量空间化的***,其特征在于,所述***还包括:
预处理模块,用于对所述NDVI进行预处理,所述预处理包括分辨率重采样、滤波去噪和平滑处理。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811226264.3A CN108984803B (zh) | 2018-10-22 | 2018-10-22 | 一种农作物产量空间化的方法及*** |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811226264.3A CN108984803B (zh) | 2018-10-22 | 2018-10-22 | 一种农作物产量空间化的方法及*** |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108984803A CN108984803A (zh) | 2018-12-11 |
CN108984803B true CN108984803B (zh) | 2020-09-22 |
Family
ID=64544524
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811226264.3A Active CN108984803B (zh) | 2018-10-22 | 2018-10-22 | 一种农作物产量空间化的方法及*** |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108984803B (zh) |
Families Citing this family (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109800921B (zh) * | 2019-01-30 | 2019-12-20 | 北京师范大学 | 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 |
CN111160781B (zh) * | 2019-12-30 | 2023-09-08 | 内蒙古鄂尔多斯资源股份有限公司东昊厂 | 一种羊绒纺织生产计算投入量的复合动态制成率模型 |
CN112084839B (zh) * | 2020-07-21 | 2024-03-05 | 沈阳农业大学 | 一种天空地一体解析小地块玉米受非生物胁迫成因的方法 |
CN111882573B (zh) * | 2020-07-31 | 2023-08-18 | 北京师范大学 | 一种基于高分辨率影像数据的耕地地块提取方法及*** |
CN113052433B (zh) * | 2021-02-22 | 2024-02-06 | 中国科学院空天信息创新研究院 | 基于关键时相和农田景观特征参量的作物单产估算方法 |
CN113269464B (zh) * | 2021-06-10 | 2024-04-23 | 中国科学院地理科学与资源研究所 | 一种生态恢复评估方法和生态恢复评估装置 |
CN114492987A (zh) * | 2022-01-24 | 2022-05-13 | 浙江大学 | 一种资产存量空间化方法、***及存储介质 |
CN114529826B (zh) * | 2022-04-24 | 2022-08-30 | 江西农业大学 | 一种基于遥感影像数据的水稻估产方法、装置及设备 |
CN116579521B (zh) * | 2023-05-12 | 2024-01-19 | 中山大学 | 产量预测时间窗口确定方法、装置、设备及可读存储介质 |
CN116542403B (zh) * | 2023-07-06 | 2023-10-20 | 航天宏图信息技术股份有限公司 | 作物产量预测方法、装置、电子设备及可读存储介质 |
Family Cites Families (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101858971A (zh) * | 2010-06-02 | 2010-10-13 | 浙江大学 | 一种基于modis数据的水稻单产遥感估算方法 |
CN101937079B (zh) * | 2010-06-29 | 2012-07-25 | 中国农业大学 | 基于区域相似度的遥感影像变化检测方法 |
CN103335953B (zh) * | 2013-06-04 | 2016-03-23 | 中国科学院遥感与数字地球研究所 | 一种个体与群体特征相结合的农作物长势遥感评价方法 |
CN103310197A (zh) * | 2013-06-13 | 2013-09-18 | 山东省农业可持续发展研究所 | 一种利用中分辨率卫星数据提取黄淮海平原区大蒜种植面积的方法 |
KR20150096103A (ko) * | 2014-02-14 | 2015-08-24 | 한국전자통신연구원 | 농산물 수확량 예측 장치 및 방법 |
CN104134095B (zh) * | 2014-04-17 | 2017-02-15 | 中国农业大学 | 一种基于尺度转换和数据同化的农作物产量估测方法 |
CN104615977B (zh) * | 2015-01-26 | 2018-02-06 | 河南大学 | 综合关键季相特征和模糊分类技术的冬小麦遥感识别方法 |
CN104794336B (zh) * | 2015-04-17 | 2017-06-27 | 武汉大学 | 一种农田秸秆资源空间分布估算方法 |
US10628895B2 (en) * | 2015-12-14 | 2020-04-21 | The Climate Corporation | Generating digital models of relative yield of a crop based on nitrate values in the soil |
CN107330801A (zh) * | 2017-06-07 | 2017-11-07 | 北京师范大学 | 一种冬小麦种植比例的计算方法及装置 |
CN107274297A (zh) * | 2017-06-14 | 2017-10-20 | 贵州中北斗科技有限公司 | 一种土地作物种植适宜性评估方法 |
CN108205718B (zh) * | 2018-01-16 | 2021-10-15 | 北京师范大学 | 一种粮食作物抽样测产方法及*** |
-
2018
- 2018-10-22 CN CN201811226264.3A patent/CN108984803B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN108984803A (zh) | 2018-12-11 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108984803B (zh) | 一种农作物产量空间化的方法及*** | |
Zhang et al. | Land cover classification of the North China Plain using MODIS_EVI time series | |
CN111598045B (zh) | 一种基于对象图谱和混合光谱的遥感耕地变化检测方法 | |
Kiptala et al. | Land use and land cover classification using phenological variability from MODIS vegetation in the Upper Pangani River Basin, Eastern Africa | |
CN107256451A (zh) | 农业大数据基于gis的农业适宜性评价分析方法 | |
Göttlicher et al. | Land‐cover classification in the Andes of southern Ecuador using Landsat ETM+ data as a basis for SVAT modelling | |
Gong et al. | Spatiotemporal dynamics of urban forest conversion through model urbanization in Shenzhen, China | |
Li et al. | Assimilating spatiotemporal MODIS LAI data with a particle filter algorithm for improving carbon cycle simulations for bamboo forest ecosystems | |
Li et al. | Mapping population density distribution at multiple scales in Zhejiang Province using Landsat Thematic Mapper and census data | |
CN108205718B (zh) | 一种粮食作物抽样测产方法及*** | |
Yang et al. | Role of the countryside landscapes for sustaining biodiversity in karst areas at a semi centennial scale | |
Gessesse et al. | Land use dynamics and base and peak flow responses in the Choke mountain range, Upper Blue Nile Basin, Ethiopia | |
CN117611993B (zh) | 一种基于遥感实际蒸散发估算植被分类的方法 | |
Cai et al. | Consistency Assessments of the land cover products on the Tibetan Plateau | |
Wu et al. | Phenology-based cropland retirement remote sensing model: A case study in Yan’an, Loess Plateau, China | |
CN113793023A (zh) | 一种基于卫星遥感的区域绿色发展指标评价方法 | |
Xing et al. | Mapping irrigated, rainfed and paddy croplands from time-series Sentinel-2 images by integrating pixel-based classification and image segmentation on Google Earth Engine | |
CN113570273A (zh) | 一种灌溉耕地统计数据的空间化方法及*** | |
Rimal et al. | Crop cycles and crop land classification in Nepal using MODIS NDVI | |
Shifaw et al. | Farmland dynamics in Pingtan, China: understanding its transition, landscape structure and driving factors | |
CN115830464A (zh) | 基于多源数据的高原山地农业大棚自动提取方法 | |
Aryastana et al. | Irrigation Water Management by Using Remote Sensing and GIS Technology to Maintain the Sustainability of Tourism Potential in Bali | |
CN113610294A (zh) | 一种中小城市色彩规划方法、存储介质和装置 | |
Siljander et al. | A predictive modelling technique for human population distribution and abundance estimation using remote-sensing and geospatial data in a rural mountainous area in Kenya | |
Shen et al. | Capturing landscape changes and ecological processes in Nikko National Park (Japan) by integrated use of remote sensing images |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |