CN102073869A - 基于点位数据和遥感影像数据的区域耕地质量监测方法 - Google Patents

Abstract

本发明公开了一种基于点位数据和遥感影像数据的区域耕地质量监测方法，包括：S1、获取NDVI时间序列数据，并进行平滑处理；S2、提取分类特征，并根据分类特征对作物种植模式分类；S3、在不同种植模式下进行NDVI时间序列主成分分析，并提取多个主成分；S4、利用点位数据，对反映耕地质量的指标进行主成分分析，提取主成分，计算耕地质量指数，进行监测点的耕地质量的定量评价；S5、利用点位数据，对不同种植模型下耕地质量指数与NDVI时间序列主成分分别进行相关性分析；S6、利用点位数据建立耕地质量指数与不同种植模式下的第一主成分建立回归模型，得到耕地质量等级分布图。该方法能从本质上全面反映耕地质量等级分布，且能获得好的区分度。

Description

基于点位数据和遥感影像数据的区域耕地质量监测方法

技术领域

本发明涉及遥感监测技术和农业资源应用技术领域。尤其涉及一种基于点位数据和遥感影像数据的区域耕地质量监测方法。

背景技术

随着经济发展、城市化进程加快以及农业利用结构调整，耕地资源的数量和质量变化都受到了影响。决策者不仅要了解耕地资源的数量，更要了解耕地质量的变化特征及规律，以保证耕地资源的可持续利用，确保粮食生态安全。因此迫切需要进行长期、稳定、全面、科学的耕地质量监测。

自我国1984年开始耕地质量监测以来的26年内，耕地质量监测都是通过设立监测点，持续定期地通过科学而合理的空间取样，来调查影响耕地质量的主要指标或属性，并结合其他数据的分析，对耕地质量的变化作出评估。虽然利用定位监测资料可以发现耕地土壤肥力、土壤污染等指标的变化，但是受监测点制约，一方面监测范围较小，另一方面监测周期长，不能实时地实现对耕地质量的监测。而且对于耕地质量的空间分布特征，则利用区域定位监测点数据通过地统计学的插值方法来预测未采样点的属性值，从而揭示区域耕地质量的空间分布特征。但是这种将点的监测结果扩展到区域的方法其精度无法保证，不能代表耕地质量大区域的空间分布特征。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何提供一种能够从本质上全面地反映耕地质量的等级分布，且能获得较好的区分度的区域耕地质量监测方法。

(二)技术方案

为解决上述技术问题，本发明提供了一种基于点位数据和遥感影像数据的区域耕地质量监测方法，包括以下步骤：

S1、获取归一化植被指数NDVI时间序列数据，并对该数据进行平滑处理；

S2、利用所述数据提取分类特征，并根据该分类特征采用决策树分类方法对作物种植模式进行分类；

S3、在不同种植模式下进行NDVI时间序列主成分分析，并提取多个主成分；

S4、利用耕地质量监测点的点位数据，对反映耕地质量的指标进行主成分分析，并提取所述主成分，利用所提取的主成分计算耕地质量指数，进行监测点的耕地质量的定量评价；

S5、利用所述点位数据，对不同种植模式下的耕地质量指数与进行主成分分析所得到的几个主成分分别进行相关性分析，若得到各种植模式下的耕地质量指数均与一主成分有一定程度的正相关关系，则将该主成分作为第一主成分进行耕地质量的监测；

S6、利用所述点位数据，建立耕地质量指数与不同种植模式下的第一主成分建立回归模型，得到耕地质量等级分布图。

其中，步骤S4中，用下式将提取的主成分进行加权平均后得到每个耕地质量监测点的耕地质量指数，其中将各主成分的特征值作为权重：

$\mathrm{Tint}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\×{\mathrm{pc}}_i/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i$

其中，Tint是利用主成分的特征值对提取的主成分加权平均后所得的反映耕地质量的指数，λ_i是提取的主成分pc_i的特征值，i表示第i个主成分，n表示所提取的主成分的个数。

其中，步骤S3中的所述主成分分析针对各种植模式的NDVI时间序列分别进行，分析时首先针对每种种植模式提取6个主成分，这6个主成分从不同方面反映了农用地生产力的状况，然后对各种植模式下第1个主成分的特征向量进行分析来认识第1个主成分代表的物理意义，主成分的名次按特征值或贡献率大小的顺序从大到小排列。

其中，步骤S1中，利用权重滑动平均滤波法Savizky-Golay对NDVI时间序列数据进行平滑滤波。

其中，利用中分辨率成像光谱仪MODIS获得NDVI时间序列数据。

其中，所述反映耕地质量的指标包括气候因素和土壤因素。

其中，所述回归模型为线性回归模型。

(三)有益效果

本发明通过点位数据和遥感影像数据相结合进行区域耕地质量监测，试验证明，该方法能够从本质上全面地反映耕地质量的等级分布，且能获得较好的区分度。

附图说明

图1是本发明的方法流程图；

图2是NDVI时间序列平滑前后对比图；

图3是决策树分类中所用到的NDVI组合图；

图4是北京市大兴区2007年作物种植模式分类图；

图5是不同种植模式的主成分分析第一主成分的特征向量曲线图；

图6是不同种植模式主成分一与耕地质量指数散点图；

图7是北京市大兴区2007年耕地质量等级分布图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细说明。以下实施例用于说明本发明，但不用来限制本发明的范围。

植被遥感信息中包含了叶子及其生长状况等综合信息，这些信息直接指示作物干物质的积累以及生物量的多少。一些研究指出NDVI与作物产量、生物量有较好的相关性，并且用于表征农田生产力。农田生产力取决于三方面因素：作物种植模式、耕地质量和农田管理措施。因此，利用遥感数据的高时间分辨率归一化植被指数(Normalized Difference Vegetation Index，NDVI)年序列与耕地质量构成关系，实现区域耕地质量的遥感监测，对区域耕地质量的实时监测具有重要的科学意义和实用价值。

如图1所示，本发明包括以下步骤：

在本发明的优选实施例中，首先利用MODIS NDVI时间序列进行作物种植模式的分类。包括MODIS NDVI时间序列数据获取及处理、作物种植模式的分类和分类精度评价。

S1、MODIS(中分辨率成像光谱仪)NDVI时间序列数据的获取与处理

首先，获取2007年MODIS250米分辨率的8d(8天)合成地表反射率产品，该数据来源于NASA的MODIS数据产品分发网站(http://modis.gsfc.nasa.gov/)，包括了第一、第二波段两个波段。利用MODIS产品专门处理工具MRT(Modis Reprojection Tool)软件对下载的产品进行文件格式、空间镶嵌和地图投影的转换，最后进行掩模得到研究区域的影像图。

其次，利用NDVI计算公式对该影像图进行计算得到MODIS NDVI时间序列数据。

NDVI是近红外与红色通道反射比值的一种变换形式。许多研究表明NDVI与LAI、绿色生物量、植被覆盖度、光合作用等植被参数有关，是植被生长状况及植被覆盖度的最佳指示因子。其计算公式为：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{NIR}}-{\mathrm{\ρ}}_R}{{\mathrm{\ρ}}_{\mathrm{NIR}}+{\mathrm{\ρ}}_R}$

式中：ρ_NIR代表近红外波段的反射率，ρ_R代表红光波段的反射率；由公式可以看出NDVI的取值范围在(-1 1)开区间。负值代表地面覆盖为云、水、雪等，对可见光高反射；0代表有岩石或裸土等，ρ_NIR和ρ_R近似相等；正值表示有植被覆盖，且随覆盖度增大而增大。在ENVI(The Environment for Visualizing Images)软件中采用band math功能对所有时期的反射率影像进行计算，并进行波段叠加，可以得到研究区每个像元的NDVI时间序列影像数据集。

最后，对时间序列NDVI进行平滑处理。本发明利用现有技术Savitzky-Golay滤波法对MODIS-NDVI时序数据进行去噪和重构处理。

本发明在ENVI中的IDL(Interactive Data Language)平台中应用Savizky-Golay滤波方法，编写程序对北京北部山区2007年NDVI时序数据逐点进行重构处理，平滑前后时序数据效果见图4，图中是某单个像元NDVI值按照一定时间顺序排列，以时间为横坐标轴排，得到的该像元NDVI时序的动态曲线。从图4中平滑前后曲线对比可以看出：平滑前时序曲线有严重的锯齿状波动，如果直接用于时间序列分析很难比较准确地捕捉作物生长趋势，对建立分类模型也造成一定的难度；重构后的生长曲线有效消除了云和缺失数据的影响，但同时基本上保持了原有曲线的基本形状，较为真实地恢复了植被指数曲线。因此，重构后的生长曲线可以更清楚地反应作物生长曲线变化趋势和规律，有利于分类模型的建立。

S2、作物种植模式分类

本发明主要利用决策树分类器对北京市大兴区的作物种植模式进行了研究，详细步骤论述如下：

根据北京市大兴区统计年鉴，确定大兴区主要的作物类型，建立大兴区作物种植模式分类体系(表1)。其中一年一熟、蔬菜瓜类、果树和一年两熟这4种种植模式都作为耕地上的种植模式。

表1北京市大兴区作物种植模式分类体系

在决策树分类中，通过对NDVI组合设置分类规则进行分类，通过边分类边学习的方法建立并改进分类规则。NDVI组合是由不同波段的NDVI影像通过加减乘除等运算而得到的一景NDVI影像，NDVI组合是为做分类所用，所以尽量选取可以明显区分出不同土地利用类型的组合。图3为分类所用到的部分NDVI组合。

大兴区属于近郊区，地块较零碎，作物类型多样。在250m的尺度下，每个象元中一般会有三种及以上类型的地类。因此导致象元的NDVI时间序列曲线多种多样，增加了分类难度。作物类型特征识别为：

1)一年两熟(小麦-玉米)：全年有植被覆被；冬季返青较早；种植较密；并且NDVI曲线有两个波峰。

2)“果树”与“蔬菜瓜类”的共同特点：NDVI曲线都表现为单峰；生长期比“玉米或花生”长；全年中NDVI在高值的时间较长。

3)“果树”与“蔬菜瓜类”的区别：在39时相(2007-11-01---2007-11-08)，瓜地菜地已经收获收割。而果树此时的NDVI值还相对较高。

4)一年一熟(玉米或花生等)：单波峰；全年生长期较短；全年NDVI平均值低。

5)其他：5月到9月的NDVI平均值远低于有作物生长的地类。

根据上述作物类型特征识别，建立分类规则，得到分类结果图，从中提取4种种植模式(图4)。

接下来还可以进行分类精度评价：

本发明抽取北京市大兴区高分辨率遥感影像SPOT5对低分辨率影像分类结果作精度检验。

用2007年2.5m分辨率的SPOT5遥感影像对2007年北京市大兴区的分类图做精度检验。通过目视解译的方法选取不同土地利用类型的训练样区，选取时量均匀选取，覆盖全区；不仅选取土地利用类型单一的训练区，并要选取多种土地利用类型混杂的训练区，当其中某土地利用类型的面积占此混杂训练区面积的1/2时即分为此类型；此外由于融合后SPOT5影像的分辨率为2.5米，而MODIS为250米，所以在SPOT5上选取的训练样区的长、宽都要大于100个栅格，才可以使SPOT5训练样区的地表实际面积与MODIS一个象元的地表实际面积相匹配。在ENVI软件中作精度检验，总体分类精度为82.78％，Kappa系数为0.77。

S3、在本发明的优选实施例中，针对不同种植模式，进行NDVI时间序列主成分分析，提取主成分。

采用主成分变换可将NDVI年序列中初级生产力的有用信号压缩到少量的前几个主分量，在信息损失最小的前提下，用较少的分量替代原来的高维数据，生成的主成分影像既包含了时间序列信息，也包含了空间信息，每个主成分表达的主题需要利用专业领域知识解释。

本发明假设这些主成分从不同的方面反映了同一作物种植模式下初级生产力的变化特征，它包含了土地自然质量、化肥投入、灌水措施等管理措施的影响。经过主成分分析，这些影响因素所造成的区域生产力的差异保持相对独立。在同一种植模式下，通常情况下土地自然质量是造成初级生产力空间差异最大的主导因素。

经过主成分分析，每种种植模式的前6个主成分在总方差中的累计贡献率如表2所示。根据方差的大小排序，每种种植模式的前6个主成分的方差最大，包含了信息量的大部分，累计达到86％以上，占有绝对优势。不同作物种植模式下，第一主成分(PC 1)(也称为主成分一)反映了初级生产力的32.77％～53.66％，蔬菜瓜类种植模式最大，其次是果园，一年两熟种植模式最小。

表2不同种植模式下主成分变换前6个主成分累积贡献率

各种植模式的PC 1特征向量如图5所示，4种种植模式下的PC 1特征向量在1年中大部分时间为正值，而在作物成熟期到衰败期期间表现为负的特征值，通常这些时期是作物经济产量形成的关键时期，一年一熟(春玉米或花生或大豆或薯类)表现为20-22个8d(薯类，6月份)或35-37个8d(春玉米或花生或大豆，10月份)，蔬菜瓜类表现在26-35个8d(7月下旬到10月初)，果树表现在21-25个8d(桑葚和桃等果树，6月中旬到7月中旬)或32-38个8d(梨和苹果等果树，9月中旬到10月下旬)，一年两熟(冬小麦和夏玉米)表现在19-25个8d(5月底到7月中旬)和33-39个8d(9月中旬到10月底)。这些时期的NDVI以负权重表征在1年的NDVI合成上，使得关键生长期的累积NDVI得到对比突出，增强了作物的成熟期与其他时间的生物量对比情况，可见4种种植模式下的第一主成分增强了作物的生长关键期对最终产量(尤其经济产量)的影响。作为耕地质量会持续稳定地支持一年中作物的生长，因此假设耕地质量的差异将会显著地在第一主成分体现。

S4、在本发明的优选实施例中，利用耕地质量监测点的点位数据，对反映耕地质量的指标进行主成分分析，并提取所述主成分，利用所提取的主成分计算每个监测点的耕地质量指数，进行监测点的耕地质量的定量评价。步骤如下：

首先是耕地质量监测点采集与分布情况。监测点采集区为大兴区六环路以南约800km²区域中的农田、园地、林区耕层(0～20cm)、亚耕层(20～40cm)。基本农田区(约300km²)按3件/km²的采样与分析密度，一般农田、园地、林区按1件/km²的采样与分析密度，共采集1240个样点。

其次构建耕地质量评价指标体系。本发明从自然环境要素功能及其彼此不能替代的原理出发，经过环境辨识、系统分析，重点选取对平原区耕地质量起决定作用的2种评价因素及其9个评价因子，构成耕地质量评价指标体系，反映耕地质量的指标包括：

(1)气候因素。水热因子对农业生产而言，对生物的正常发育和生长起着决定性的作用，温度与水分及其对比状况不仅决定作物的种类、熟制、产量和品质，而且在很大程度上也决定土地用途种类和发展农业生产中应采取的方向性措施，它是农用地评价的基本因素。反映土地质量的气候因素指标：温度生长期和水分生长期。温度是能量的重要表现形式，它决定着作物整个生长发育阶段的生命活力，温度生长期是联合国粮农组织在农业生态区研究中提出的反应土地质量的新概念，它不针对具体作物而言，只反映某一地区土地资源的热量状况。水分生长期像温度生长期一样，是描述土地特性的指标，它是当地主导土壤、作物种植制度下，不靠灌溉的参照生长期长度，与具体作物和土壤无关。

(2)土壤因素。对农用地来说土壤是其重要组成因素之一，农用地的生产力主要取决于土壤的肥力水平，即土壤能够供应和协调植物生长发育所必须的水分、养分、空气和热量以及其他生活条件的能力，它直接影响作物的生命活动，决定作物的产量和质量，形成耕地质量的差异，它是耕地评价的最基本因素。选用的反映土壤因素的指标包括：土壤全氮、土壤水解性氮、有机质、有效磷、速效钾、含水量和PH值。

然后，对上述指标进行主成分分析，解释从变量中观察到的大部分变化。基于变量相关矩阵变换，得到主成分的特征值、方差贡献率、累积贡献率和特征向量，提取特征值超过1.0的成分，建立主成分方程，计算这些主成分的得分。在1240个样点中，从评价指标中提取了3个主成分，它们的特征值和特征向量见表3和表4。前3个主成分共同说明了耕地质量的73.86％，尤其是第一主成分解释了耕地质量的48.36％。第一和第二主成分反映了耕地质量的土壤养分状况；第三主成分反映了耕地质量的水热状况。

表3耕地质量评价指标主成分分析特征值

主成分

特征值

方差贡献率％

累计贡献率％

1	4.353	48.36	48.36
				2	1.159	12.88	61.24
3	1.136	12.62	73.86
				4	0.696	7.73	81.56
5	0.599	6.65	88.25
				6	0.449	4.99	93.25
7	0.328	3.64	96.89
				8	0.236	2.63	99.51
9	0.044	0.49	100.000

表4耕地质量评价指标主成分分析前三个成分的特征向量

指标	第一特征向量	第二特征向量	第三特征向量
				全氮	0.450	-0.057	0.054
水解性氮	0.388	-0.206	0.107
				有机质	0.429	-0.087	-0.017
有效磷	0.184	-0.576	0.449
				速效钾	0.297	0.321	0.107
含水率	0.370	0.364	-0.113
				PH	-0.394	0.076	0.048
十度温度生长期	0.033	0.597	0.618
				水分生长期	-0.217	-0.136	0.612

最后，用下式将提取的主成分进行加权平均后得到每个采样点(即耕地质量监测点)的耕地质量指数，其中把成分的特征值作为权重。

$\mathrm{Tint}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\×{\mathrm{pc}}_i/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i---\left(6\right)$

这里，Tint是利用主成分的特征值对提取的主成分加权平均后所得到的反映耕地质量的指数(即耕地质量指数)，λ_i是提取的主成分pc_i的特征值，i表示第i个主成分，n表示所提取的主成分的个数。

在本发明的优选实施例中，利用耕地质量指数与不同种植模式下进行主成分分析得到的几个主成分的关系计算区域耕地质量指数，实现区域耕地质量监测。步骤如下：

S5、首先是不同种植模式下的耕地质量指数与进行主成分分析所得到的几个主成分的相关性分析。将耕地质量的点位数据与各类种植模式叠加，确定每个耕地点位的种植模式，其中，有589个点属于一年一熟种植模式，102个点属于蔬菜瓜类种植模式，110个点属于果园种植模式，136个点属于一年两熟种植模式。然后，在各类种植模式下，提取各点位数据所在栅格像元的每个主成分的值，对点位数据的耕地质量与每个主成分分值进行相关性分析，得到相关系数(表5)。

表5耕地质量与不同种植模式主成分的相关系数

表中，**表示显著水平0.01；*表示显著水平0.05。

在4类种植模式下，每类种植模式NDVI时间序列的PC1都与耕地质量呈0.01水平的显著正相关关系，这说明第一主成分可以反映耕地质量差异影响的生产力。因此可以用第一主成分进行耕地质量的监测。

S6、利用所述点位数据，建立耕地质量指数与不同种植模式PC1的回归模型。在4类种植模式下，绘制点位数据PC 1与耕地质量指数的散点图(图6)，各种植模式下的耕地质量与由耕地质量所造成的生产力呈显著地线性相关关系。通过各种植模式下的线性关系，建立不同种种模式的线性回归模型(表6)，计算每个像元(250m*250m)的耕地质量指数，利用表7对耕地质量指数进行分级，还可以得到耕地质量等级图(图7)，等级I到VIII表示耕地质量越来越差。

表6不同种植模式耕地质量指数的线性回归模型

表7耕地质量等级划分

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (7)

1.一种基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，包括以下步骤：

S1、获取归一化植被指数NDVI时间序列数据，并对该数据进行平滑处理；

S2、利用所述数据提取分类特征，并根据该分类特征采用决策树分类方法对作物种植模式进行分类；

S3、在不同种植模式下进行NDVI时间序列主成分分析，并提取多个主成分；

S4、利用耕地质量监测点的点位数据，对反映耕地质量的指标进行主成分分析，并提取所述主成分，利用所提取的主成分计算耕地质量指数，进行监测点的耕地质量的定量评价；

S5、利用所述点位数据，对不同种植模式下的耕地质量指数与进行主成分分析所得到的几个主成分分别进行相关性分析，若得到各种植模式下的耕地质量指数均与一主成分有一定程度的正相关关系，则将该主成分作为第一主成分进行耕地质量的监测；

S6、利用所述点位数据，建立耕地质量指数与不同种植模式下的第一主成分建立回归模型，得到耕地质量等级分布图。

2.如权利要求1所述的基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，步骤S4中，耕地质量定量评价时利用主成分分析提取主成分，用下式将提取的主成分进行加权平均后得到每个耕地质量监测点的耕地质量指数，其中将各主成分的特征值作为权重：

$\mathrm{Tint}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\×{\mathrm{pc}}_i/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i$

其中，Tint是利用主成分的特征值对提取的主成分加权平均后所得到的反映耕地质量的指数，λ_i是提取的主成分pc_i的特征值，i表示第i个主成分，n表示所提取的主成分的个数。

3.如权利要求1所述的基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，步骤S3中的所述主成分分析针对各种植模式的NDVI时间序列分别进行，分析时首先针对每种种植模式提取6个主成分，这6个主成分从不同方面反映了农用地生产力的状况，然后对各种植模式下第1个主成分的特征向量进行分析来认识第1个主成分代表的物理意义，主成分的名次按特征值或贡献率大小的顺序从大到小排列。

4.如权利要求1所述的基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，步骤S1中，利用权重滑动平均滤波法Savizky-Golay对NDVI时间序列数据进行平滑滤波。

5.如权利要求1所述的基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，利用中分辨率成像光谱仪MODIS获得NDVI时间序列数据。

6.如权利要求1所述的基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，所述反映耕地质量的指标包括气候因素和土壤因素。

7.如权利要求1～6任一项所述的基于点位数据和遥感影像数据的区域耕地质量监测方法，其特征在于，所述回归模型为线性回归模型。

Images