CN113868826B - 一种基于栅格尺度耕地农田生产能力重心变化的测度方法 - Google Patents

Abstract

本申请公开了一种基于栅格尺度耕地农田生产能力重心变化的测度方法，包括：步骤1，将待测度地区所在区域栅格化，并基于各年度的土地利用数据，进行耕地区域标注，生成栅格耕地数据；步骤2，获取待测度地区各个年度的累计植被指数数据，并结合对应年度的栅格耕地数据，基于地理位置存在的交集区域，进行数据提取，生成待测度地区各年度的累计耕地植被栅格数据；步骤3，以年度为单位，依次计算每一年度的累计耕地植被栅格数据的耕地重心经纬度坐标，基于每一年度的耕地重心经纬度坐标，对待测度地区的农田生产能力重心变化的速度和趋势进行预测分析。通过本申请中的技术方案，优化耕地或粮食生产重心经纬度坐标计算结果，使其与现实结果更符合。

Description

一种基于栅格尺度耕地农田生产能力重心变化的测度方法

技术领域

本申请涉及农业地理学、土地资源管理、耕地变化预测的技术领域，具体而言，涉及一种基于栅格尺度耕地农田生产能力重心变化的测度方法。

背景技术

耕地资源是农业生产最基本的物质条件，其质量和数量的变化会直接影响到粮食的供应格局和安全水平。随着城市化速度的不断加快，部分地区，尤其是经济较为发达的东部地区，耕地非农化趋势明显且短期内难以逆转。耕地在不同地区的增减变化，直接影响到我国耕地的空间分布格局，进而影响粮食的种植与生产格局。因此，掌握耕地的空间分布格局及其变化，有助于帮助我们了解粮食生产格局的状况，对于及时制定相关应对策略具有重要的现实意义。

传统的耕地研究往往利用统计年鉴得到的耕地面积和粮食产量数据，无法知晓行政单元内部耕地和产量的具体空间分布。

而现有的耕地或粮食生产重心，通常是基于地理信息系统专业制图软件ArcGIS中的Feature to point工具计算得到的，基于将所有代表各个地区的矢量多边形合并为一个，以计算全国或者某一指定区域耕地或粮食生产重心。

在实际计算过程中，往往由于矢量多边形过多、过细碎，计算耕地或粮食生产重心时费时费力，且计算人员必须有专门的经验知识，才能保证最终计算结果的准确性。

另外，在研究耕地或粮食生产重心时，目前使用的统计数据的最小统计单元往往是县域，一方面，无法对县域级别下的地区进行耕地或粮食生产重心的计算；另一方面，这种计算结果通常是假设所统计的各个县域内耕地或粮食生产能力一致而得出的，导致计算出的结果与实际的情况存在一定的误差，并且，当县域面积越大，这种计算结果的误差越大。

发明内容

本申请的目的在于：提高耕地或粮食生产重心经纬度坐标计算结果的准确度和计算效率，使其与现实统计结果更符合，进而优化待测度地区的农田生产能力重心变化的速度和趋势进行预测分析。

本申请的技术方案是：提供了一种基于栅格尺度耕地农田生产能力重心变化的测度方法，该方法包括：步骤1，将待测度地区所在区域栅格化，并基于各年度的土地利用数据，进行耕地区域标注，生成栅格耕地数据；步骤2，获取待测度地区各个年度的累计植被指数数据，并结合对应年度的栅格耕地数据，基于地理位置存在的交集区域，进行数据提取，生成待测度地区各年度的累计耕地植被栅格数据，其中，累计植被指数数据为栅格型数据；步骤3，以年度为单位，依次计算每一年度的累计耕地植被栅格数据的耕地重心经纬度坐标，基于每一年度的耕地重心经纬度坐标，对待测度地区的农田生产能力重心变化的速度和趋势进行预测分析。

上述任一项技术方案中，进一步地，步骤2中，获取待测度地区各个年度的累计植被指数数据，具体包括：

步骤201，获取待测度地区各年度的连续5日平均温度数据；步骤202，根据预设积温阈值，采用累加求和的方式，计算待测度地区各年度以日为单位的归一化植被指数NDVI的第一和值，将第一和值记作累计植被指数数据，累计植被指数数据的计算公式为：

式中，N_r为说是累计植被指数数据，NDVI_i为第i天的归一化植被指数NDVI的取值，SDT为连续5日平均温度数据大于预设积温阈值的起始日期，EDT为连续5日平均温度数据大于预设积温阈值的截至日期；

上述任一项技术方案中，进一步地，栅格耕地数据中耕地区域对应栅格内数据的取值被标注为1，栅格耕地数据中非耕地区域对应栅格内数据的取值被标注为0。

上述任一项技术方案中，进一步地，步骤2中，生成待测度地区各年度的累计耕地植被栅格数据，具体包括：步骤211，选取累计植被指数数据和栅格耕地数据中地理位置上存在交集的区域作为掩膜，并生成初始耕地植被栅格数据，其中，初始耕地植被栅格数据中各个栅格对应的数据初始值为0；步骤212，基于栅格耕地数据中被标注为耕地区域的栅格标号，选取累计植被指数数据中对应栅格标号的指数数据；步骤213，根据指数数据和栅格标号，对初始耕地植被栅格数据中各个栅格对应的数据初始值进行更新，生成累计耕地植被栅格数据。

上述任一项技术方案中，进一步地，步骤3中，计算累计耕地植被栅格数据的耕地重心的方法，具体包括：步骤31，采用累加求和方式，计算累计耕地植被栅格数据中各个栅格对应数据的第二和值，将第二和值记作耕地总面积；步骤32，逐行计算累计耕地植被栅格数据中每一行栅格对应数据的第三和值与对应行序号的乘积，记作单行耕地面积；步骤33，计算各个单行耕地面积的第四和值，以及第四和值与耕地总面积的第一商，将第一商记作耕地重心行坐标；步骤34，逐列计算累计耕地植被栅格数据中每一列栅格对应数据的第五和值与对应列序号的乘积，记作单列耕地面积；步骤35，计算各个单列耕地面积的第六和值，以及第六和值与耕地总面积的第二商，将第二商记作耕地重心列坐标；步骤36，根据耕地重心行坐标与耕地重心列坐标，生成耕地重心经纬度坐标。

上述任一项技术方案中，进一步地，耕地重心行坐标的计算公式为：

式中，X为耕地重心行坐标，i为累计耕地植被栅格数据中的行标号，i＝1,2,...,rows，rows为累计耕地植被栅格数据中行数最大值，j为累计耕地植被栅格数据中的列标号，j＝1,2,...,cols，cols为累计耕地植被栅格数据中列数最大值，为累计耕地植被栅格数据中第i行第j列对应的数据，cropland_all为耕地总面积。

上述任一项技术方案中，进一步地，耕地重心列坐标的计算公式为：

式中，Y为耕地重心列坐标，i为累计耕地植被栅格数据中的行标号，i＝1,2,...,rows，rows为累计耕地植被栅格数据中行数最大值，j为累计耕地植被栅格数据中的列标号，j＝1,2,...,cols，cols为累计耕地植被栅格数据中列数最大值，为累计耕地植被栅格数据中第i行第j列对应的数据，cropland_all为耕地总面积。

本申请的有益效果是：

本申请中的技术方案，通过栅格化的方式，基于各年度的土地利用数据，生成栅格耕地数据，并与累计植被指数数据相结合，将粮食产量数据在栅格尺度上进行体现，可以准确地获得待测度地区的耕地或粮食生产能力的具体地域空间分布状况，进而准确地掌握粮食生产在地域空间上的分布特点，解决了传统方法在研究耕地中无法确定耕地的空间分布格局及变化趋势的问题，并提高了耕地或粮食生产重心的计算效率和准确度，降低了计算过程中的人力成本，且对操作人员的要求较矢量数据处理相对更低。

本申请中采用栅格化的方式，相对于传统统计县域尺度而言，可以更加详细的了解耕地及粮食生产的位置信息和数量特征，其计算出的重心变化结果与现实结果更符合。

附图说明

本申请的上述和/或附加方面的优点在结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是根据本申请的一个实施例的基于栅格尺度耕地农田生产能力重心变化的测度方法的示意流程图；

图2是根据本申请的一个实施例的基于Matlab实现的示意流程图；

图3是根据本申请的一个实施例的农田生产力的示意图；

图4是根据本申请的一个实施例的累计耕地植被栅格数据计算的示意图；

图5是根据本申请的一个实施例的耕地重心计算的示意流程图；

图6是根据本申请的一个实施例的某区域耕地重心变化的仿真示意图。

具体实施方式

为了能够更清楚地理解本申请的上述目的、特征和优点，下面结合附图和具体实施方式对本申请进行进一步的详细描述。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互结合。

在下面的描述中，阐述了很多具体细节以便于充分理解本申请，但是，本申请还可以采用其他不同于在此描述的其他方式来实施，因此，本申请的保护范围并不受下面公开的具体实施例的限制。

如图1所示，本实施例提供了一种基于栅格尺度耕地农田生产能力重心变化的测度方法，该方法包括：

步骤1，将待测度地区所在区域栅格化，并基于各年度的土地利用数据，进行耕地区域标注，生成栅格耕地数据，其中，栅格耕地数据中耕地区域对应栅格内数据的取值被标注为1，栅格耕地数据中非耕地区域对应栅格内数据的取值被标注为0。

具体的，为了实现本实施例中的测度方法，以Matlab软件为基础，搭建本实施例的实现框架，如图2所示。

由于Matlab软件能够快速地处理矩阵数据，因此，将待测度地区的相关数据进行栅格(网格)化，之后输入至Matlab软件时，即可转化为相应的矩阵形式，便于Matlab软件进行处理。

欧洲航天局全球土地覆盖产品(ESA CCI-LC)、美国宇航局和美国地质调查局全球土地覆盖产品(NASA/USGS MODIS-LC)每年都会对全球土地使用情况据进行统计，同样的，中国科学院资源环境科学与数据中心每2-3年也会对中国土地使用情况据进行统计。

本实施例以ESA CCI-LC逐年数据为基础，作为待测度地区的土地利用数据，将待测度地区根据设定的栅格边长进行栅格化，若某栅格对应区域为耕地区域时，将栅格耕地数据中对应的数值置为1，反之，若为非耕地区域，则将数值置为0。

需要说明的是，栅格耕地数据为矩阵形式。

步骤2，获取待测度地区各个年度的累计植被指数数据，并结合对应年度的栅格耕地数据，基于地理位置存在的交集区域，进行数据提取，生成待测度地区各年度的累计耕地植被栅格数据，其中，累计植被指数数据为栅格型数据；

具体的，如图3所示，本实施例将归一化植被指数NDVI作为反映农作物长势和营养信息的重要信息，用来反映农田或者耕地的生产能力，将其与栅格数据相结合，以便用Matlab软件进行处理，掌握粮食生产在空间上的分布特点。

进一步的，本实施例中采用的归一化植被指数NDVI为逐日进行统计的数据，而基于对农作物生长所需温度因素的分析，当环境温度较低时(如温度低于10℃)，农作物基本不会生长，因此，可以基于该地区的环境温度，剔除归一化植被指数NDVI中的部分数据，以减少计算量，剔除异常数据、提高计算精度。因此，该步骤2中，获取待测度地区各个年度的累计植被指数数据，具体包括：

步骤201，获取待测度地区各年度的连续5日平均温度数据；

步骤202，根据预设积温阈值，采用累加求和的方式，计算待测度地区各年度以日为单位的归一化植被指数NDVI的第一和值，将第一和值记作累计植被指数数据，累计植被指数数据的计算公式为：

式中，N_r为说是累计植被指数数据，NDVI_i为第i天的归一化植被指数NDVI数值，SDT为连续5日平均温度数据大于预设积温阈值的起始日期，EDT为连续5日平均温度数据大于预设积温阈值的截至日期。

需要说明的是，本实施例设定预设积温阈值的取值为10℃，起始日期SDT、截至日期EDT为气象学中的界限温度起止日期，即起始日期SDT为在日平均积温数据中，连续5天的滑动平均温度值大于预设积温阈值10℃的第一天；截至日期EDT为在日平均积温数据中，连续5天的滑动平均温度值大于预设积温阈值10℃的最后一天。

该NDVI_i的取值由对16天时间间隔的归一化植被指数NDVI进行线性插值计算所得，具体线性插值计算过程本实施例并不限定。

进一步的，步骤2中，生成待测度地区各年度的累计耕地植被栅格数据，具体包括：

步骤211，选取累计植被指数数据和栅格耕地数据中地理位置上存在交集的区域作为掩膜，并生成初始耕地植被栅格数据，其中，初始耕地植被栅格数据中各个栅格对应的数据初始值为0；

具体的，在地理信息系统软件ArcGIS中，所使用的地理数据都带有地理坐标系、投影坐标系，该软件ArcGIS中，默认的是栅格耕地数据的第i行j列对应的地理区域与累计植被指数数据的第i行j列对应的地理区域相同。

并且，如果出现数据区域不对应的情况，如某地区对应的栅格耕地数据为5*4，而累计植被指数数据为3*6，则可以使用软件ArcGIS自身具有的Arctoolbox工具箱中功能Extract by mask工具，直接基于地理位置上的交集区域，提取两个数据的交集区域，保持两类数据的行列数一致，以便后续Matlab对数据进行读取。

步骤212，基于栅格耕地数据中被标注为耕地区域的栅格标号，选取累计植被指数数据中对应栅格标号的指数数据；

步骤213，根据指数数据和栅格标号，对初始耕地植被栅格数据中各个栅格对应的数据初始值进行更新，生成累计耕地植被栅格数据。

如图4所示，在Matlab中建立空数据立方体或三维数组，以三维数组为例，其行列数为掩膜的行列数。将栅格耕地数据的栅格标号与第一层三维数组的行列号一一对应，同样的，将累计植被指数数据的行列号与第二层三维数组的行列号一一对应，之后，使用find函数，对第一层三维数组中耕地区域(即取值被标注为1的区域)进行标记，根据标记结果，选取第二层三维数组中相同位置的数据，对第三层三维数组中的数据进行更新，将耕地区域与农田或者耕地的生产能力进行对应，得到累计耕地植被栅格数据，即累计耕地植被栅格数据中每一个栅格的数值代表该区域耕地的生产能力。

步骤3，以年度为单位，依次计算每一年度的累计耕地植被栅格数据的耕地重心经纬度坐标，基于每一年度的耕地重心，对待测度地区的农田生产能力重心变化的速度和趋势进行预测分析。

进一步的，步骤3中，计算累计耕地植被栅格数据的耕地重心经纬度坐标的方法，具体包括：

步骤31，采用累加求和方式，计算累计耕地植被栅格数据中各个栅格对应数据的第二和值，将第二和值记作耕地总面积cropland_all；

步骤32，逐行计算累计耕地植被栅格数据中每一行栅格对应数据的第三和值与对应行序号的乘积，记作单行耕地面积；

步骤33，计算各个单行耕地面积的第四和值，以及第四和值与耕地总面积cropland_all的第一商，将第一商记作耕地重心行坐标，对应的计算公式为：

式中，X为所述耕地重心行坐标，i为累计耕地植被栅格数据中的行标号，i＝1,2,...,rows，rows为累计耕地植被栅格数据中行数最大值，j为累计耕地植被栅格数据中的列标号，j＝1,2,...,cols，cols为累计耕地植被栅格数据中列数最大值，为累计耕地植被栅格数据中第i行第j列对应的数据；

步骤34，逐列计算累计耕地植被栅格数据中每一列栅格对应数据的第五和值与对应列序号的乘积，记作单列耕地面积；

步骤35，计算各个单列耕地面积的第六和值，以及第六和值与耕地总面积cropland_all的第二商，将第二商记作耕地重心列坐标，对应的计算公式为：

式中，Y为所述耕地重心列坐标，i为累计耕地植被栅格数据中的行标号，i＝1,2,...,rows，rows为累计耕地植被栅格数据中行数最大值，j为累计耕地植被栅格数据中的列标号，j＝1,2,...,cols，cols为累计耕地植被栅格数据中列数最大值，为累计耕地植被栅格数据中第i行第j列对应的数据；

步骤36，根据耕地重心行坐标与耕地重心列坐标，生成耕地重心经纬度坐标(X，Y)。

具体的，如图5所示，利用Matlab软件中的geotiffread函数，得到的耕地上NDVI数据(累计耕地植被栅格数据)导入Matlab中，生成二维数组a，并使用geotiffinfo函数读取数据的投影信息，方便后续输出为TIFF数据；之后，使用size函数读取二维数组a列数(row,col)，并建立一个值为0，行列数为row和col的二维数组b。

基于累计耕地植被栅格数据，计算对应的耕地重心行坐标X和耕地重心列坐标Y，将二维数组b中(X,Y)位置的数值赋值为1，最后，使用geotiffwrite函数，将二维数组输出为TIFF文件。

为便于查看，以5年为时间间隔制图，如图6所示，展示了2000-2020年中国农田生产能力重心分布图。通过本实施例中的上述过程，将2000、2005、2010、2015和2020年耕地重心经纬度坐标借助软件ArcGIS可视化在地图中，如点A至点E所示。

通过耕地重心经纬度坐标的变化可以发现，近20年耕地农田生产力的重心呈现出向西北迁移的趋势，这可能主要是由以下两方面的原因引起：

一是直接原因：北方地区，特别是西北的新疆地区，耕地面积持续增加，与此同时节水灌溉技术在西北地区不断推广使用，提高了西北地区耕地农田生产力；

二是间接原因：由于农业生产的比较效益较低，传统从事农业活动的人们转而进城从事务工劳动，引起部分东南方地区出现耕地撂荒或原本一年种植两季的地区仅种植一季，引起耕地单位面积农田生产量下降，这间接的导致了耕地农田生产力向西北移动。

此外，耕地重心经纬度坐标出现年际波动，可能是由于东南地区本身拥有相对好的水热资源且适宜农作物一年多季生长，粮食价格以及老百姓的种植意愿(种植一季或种植两季)直接会影响到耕地在一年的种植和使用频率。当东南地区耕地被充分利用时，耕地重心会出现小范围的向南波动。但由于近年来北方地区耕地持续开垦，因此耕地农田生产重心总体向西北地区移动的趋势不会出现大的变化。

然而，必须注意的是，耕地粮食生产重心的西北移动趋势也需要引起相关的关注和重视，尽管中国北方地区后备耕地资源相对充裕，但水资源相对不足且生态系统相对脆弱，耕地开发与利用的同时需要兼顾保护生态与环境，才能维持北方地区土地的可持续性利用与发展。

以上结合附图详细说明了本申请的技术方案，本申请提出了一种基于栅格尺度耕地农田生产能力重心变化的测度方法，包括：步骤1，将待测度地区所在区域栅格化，并基于各年度的土地利用数据，进行耕地区域标注，生成栅格耕地数据；步骤2，获取待测度地区各个年度的累计植被指数数据，并结合对应年度的栅格耕地数据，基于地理位置存在的交集区域，进行数据提取，生成待测度地区各年度的累计耕地植被栅格数据；步骤3，以年度为单位，依次计算每一年度的累计耕地植被栅格数据的耕地重心经纬度坐标，基于每一年度的耕地重心经纬度坐标，对待测度地区的农田生产能力重心变化的速度和趋势进行预测分析。通过本申请中的技术方案，优化耕地或粮食生产重心经纬度坐标计算结果，使其与现实结果更符合。

本申请中的步骤可根据实际需求进行顺序调整、合并和删减。

本申请装置中的单元可根据实际需求进行合并、划分和删减。

尽管参考附图详地公开了本申请，但应理解的是，这些描述仅仅是示例性的，并非用来限制本申请的应用。本申请的保护范围由附加权利要求限定，并可包括在不脱离本申请保护范围和精神的情况下针对发明所作的各种变型、改型及等效方案。

Claims (6)

1.一种基于栅格尺度耕地农田生产能力重心变化的测度方法，其特征在于，所述方法包括：

步骤1，将待测度地区所在区域栅格化，并基于各年度的土地利用数据，进行耕地区域标注，生成栅格耕地数据；

步骤2，获取所述待测度地区各个年度的累计植被指数数据，并结合对应年度的栅格耕地数据，基于地理位置存在的交集区域，进行数据提取，生成所述待测度地区各年度的累计耕地植被栅格数据；

步骤3，以年度为单位，依次计算每一年度的累计耕地植被栅格数据的耕地重心经纬度坐标，基于每一年度的耕地重心经纬度坐标，对所述待测度地区的农田生产能力重心变化的速度和趋势进行预测分析，

所述步骤3中，计算累计耕地植被栅格数据的耕地重心的方法，具体包括：

步骤31，采用累加求和方式，计算所述累计耕地植被栅格数据中各个栅格对应数据的第二和值，将所述第二和值记作耕地总面积；

步骤32，逐行计算所述累计耕地植被栅格数据中每一行栅格对应数据的第三和值与对应行序号的乘积，记作单行耕地面积；

步骤33，计算各个单行耕地面积的第四和值，以及所述第四和值与所述耕地总面积的第一商，将所述第一商记作耕地重心行坐标；

步骤34，逐列计算所述累计耕地植被栅格数据中每一列栅格对应数据的第五和值与对应列序号的乘积，记作单列耕地面积；

步骤35，计算各个单列耕地面积的第六和值，以及所述第六和值与所述耕地总面积的第二商，将所述第二商记作耕地重心列坐标；

步骤36，根据所述耕地重心行坐标与所述耕地重心列坐标，生成所述耕地重心经纬度坐标。

2.如权利要求1所述的基于栅格尺度耕地农田生产能力重心变化的测度方法，其特征在于，所述步骤2中，获取所述待测度地区各个年度的累计植被指数数据，具体包括：

步骤201，获取所述待测度地区各年度的连续5日平均温度数据；

步骤202，根据预设积温阈值，采用累加求和的方式，计算所述待测度地区各年度以日为单位的归一化植被指数NDVI的第一和值，将所述第一和值记作所述累计植被指数数据，所述累计植被指数数据的计算公式为：

式中，N_r为说是累计植被指数数据，NDVI_i为第i天的归一化植被指数NDVI的数值，SDT为连续5日平均温度数据大于所述预设积温阈值的起始日期，EDT为连续5日平均温度数据大于所述预设积温阈值的截至日期。

3.如权利要求1所述的基于栅格尺度耕地农田生产能力重心变化的测度方法，其特征在于，所述栅格耕地数据中耕地区域对应栅格内数据的取值被标注为1，所述栅格耕地数据中非耕地区域对应栅格内数据的取值被标注为0。

4.如权利要求1所述的基于栅格尺度耕地农田生产能力重心变化的测度方法，其特征在于，所述步骤2中，生成所述待测度地区各年度的累计耕地植被栅格数据，具体包括：

步骤211，选取所述累计植被指数数据和所述栅格耕地数据中地理位置上存在交集的区域作为掩膜，并生成初始耕地植被栅格数据，其中，所述初始耕地植被栅格数据中各个栅格对应的数据初始值为0；

步骤212，基于所述栅格耕地数据中被标注为耕地区域的栅格标号，选取所述累计植被指数数据中对应栅格标号的指数数据；

步骤213，根据所述指数数据和所述栅格标号，对所述初始耕地植被栅格数据中各个栅格对应的数据初始值进行更新，生成所述累计耕地植被栅格数据。

5.如权利要求1所述的基于栅格尺度耕地农田生产能力重心变化的测度方法，其特征在于，所述耕地重心行坐标的计算公式为：

式中，X为所述耕地重心行坐标，i为所述累计耕地植被栅格数据中的行标号，i＝1,2,…,rows，rows为所述累计耕地植被栅格数据中行数最大值，j为所述累计耕地植被栅格数据中的列标号，j＝1,2,…,cols，cols为所述累计耕地植被栅格数据中列数最大值，为所述累计耕地植被栅格数据中第i行第j列对应的数据，croland_all为所述耕地总面积。

6.如权利要求1所述的基于栅格尺度耕地农田生产能力重心变化的测度方法，其特征在于，所述耕地重心列坐标的计算公式为：

式中，Y为所述耕地重心列坐标，i为所述累计耕地植被栅格数据中的行标号，i＝1,2,…,rows，rows为所述累计耕地植被栅格数据中行数最大值，j为所述累计耕地植被栅格数据中的列标号，j＝1,2,…,cols，cols为所述累计耕地植被栅格数据中列数最大值，为所述累计耕地植被栅格数据中第i行第j列对应的数据，cropland_all为所述耕地总面积。