一种基于MODISNDVI时序数据的耕地损失评估方法与流程

一种基于modis ndvi时序数据的耕地损失评估方法
技术领域
1.本发明属于遥感图像处理以及遥感变化检测领域，利用modis ndvi时序数据进行由于城市扩张导致的耕地损失变化检测，获取由耕地转变为建设用地的具体时间，特别指一种基于modis ndvi时序数据的耕地损失评估方法。


背景技术：

2.2001
‑
2005年间，城市扩张占用的农业用地达到了每年2178km2，且占用的大多是平原地区或居民点周围的优质耕地；中国发展报告2010预测，未来20
‑
30年，中国人口规模达到顶峰，城市化率将达到65％。快速的城市化进程使得农业用地将会进一步大面积缩减，对我国粮食安全构成潜在威胁。准确、客观的获取每年由于城市扩张导致的农业用地损失状况，能为城市发展、土地管理等政策的制定提供有效的科学依据。然而，常规的统计农业用地现状的方法为土地调查，然而土地调查耗费大量财力人力，不可能频繁进行；同时土地调查持续周期长，数据具有不同程度的滞后性。遥感则具有大范围、周期短、连续性等优点，能够很好的弥补土地调查的缺点。
3.当前，基于遥感影像分析耕地变化的方式主要是变化检测方法，许多研究人员提出了遥感影像变化检测模型，总结如下：
4.1)基本的土地利用变化检测方法是多时相遥感影像分类法，即将多期遥感影像进行地表覆盖分类，然后提取变化区域，分类方法包括支撑向量机、最大似然分类、深度学习方法等，然而，由于遥感影像场景复杂，目前分类算法精度仍然不高，且单一分类模型无法适用于所有地区；
5.2)第二种变化检测方法是基于多时相遥感影像的主成分分析(pca)，以减少因分类误差的影响，该方法不依赖于图像分类，但受物候影响，相同地表覆盖在影像上的特征会随着季节变化，所以该方法经常时效；
6.3)目前，随着深度学习的发展，许多研究人员提出了基于深度学习的遥感影像变化检测方法，如基于resnet的变化检测方法，该类方法收到样本数量限制，无法适用于所有地区，且精度有限。
7.变化检测方法仅能获知耕地的用地类型是否发生变化，无法确定耕地发生变化的具体时间。


技术实现要素：

8.本发明所要解决的技术问题在于提取由于城市扩张导致的耕地损失时空变化情况，提出一种基于modis ndvi时序数据的耕地损失评估方法，该方法不仅能提取由耕地转变为建设用地的像素，还能获取用地类型变化的时间。
9.本发明采用的技术方案为：
10.一种基于modis ndvi时序数据的耕地损失评估方法，包括以下步骤：
11.(1)获取modis ndvi时序数据，分别转换投影至地理坐标系下，并按照日期叠加在
一起生成时序数据；
12.(2)利用现有土地利用产品，提取出第一期土地覆盖类型为耕地，后一期土地覆盖类型为建设用地的像素；
13.(3)对于步骤(2)中提取出的像素，利用时序模型逐像素进行拟合；
14.(4)根据拟合结果，利用aic准则判断modis ndvi时序数据中是否存在断点以及断点时刻；若存在断点，则判断地表覆盖类型是否从耕地转变为建筑用地。
15.进一步的，步骤(3)中的时序模型为：
16.y
v
＝t
v
s
v
r
v
ꢀꢀ
(1)
17.其中：
[0018][0019]
t
v
＝α βt
ꢀꢀ
(3)
[0020]
式中，t
v
代表趋势性成分，由线性函数表达；s
v
代表季节性成分，由三角函数的傅里叶级数表达；r
v
代表残差，α、β为线性函数的系数，γ
k
和θ
k
为系数，f为每年的modis ndvi观测数量，k为确定三角函数周期的系数，t为自变量，代表第t个观测值。
[0021]
进一步的，步骤(3)中拟合的方式具体为：
[0022]
若地表覆盖不发生变化，则modis ndvi时序数据不存在断点，故：
[0023][0024]
拟合结果得到残差平方和rss1；
[0025]
若地表覆盖发生变化，则modis ndvi时序数据中存在1个断点，设断点位置为t
*
，f＜t
*
＜l
‑
f，故：
[0026][0027]
式中，f取值为一个周期，l为modis ndvi时序数据观测总数量，α1、β1、α2、β2分别为两个线性函数的系数，γ
i
、θ
i
、γ
j
、θ
j
为系数，i、j为断点前后确定三角函数周期的系数，t
*
为断点位置，t为自变量，代表第t个观测值；
[0028]
利用公式(5)遍历t
*
的所有值，求取残差平方和的最小值rss2，并记录下此时t
*
的取值即为断点时刻。
[0029]
进一步的，步骤(4)中利用aic准则判断modis ndvi时序数据中是否存在断点以及断点时刻，具体为：
[0030]
aic准则表达式如下：
[0031]
aic＝2l lln(rss/l)
ꢀꢀ
(6)
[0032]
式中，l为公式(4)或公式(5)中的参数个数，将公式(4)或公式(5)中的参数个数代入公式(6)分别对应得到两个aic值aic1和aic2，若aic1>aic2，则时序数据中存在断点，断点时刻为断点时刻即为地表覆盖类型发生变化的时刻。
[0033]
进一步的，步骤(4)中判断地表覆盖类型是否从耕地转变为建筑用地，具体为：
[0034]
若地表覆盖类型从耕地转变为建筑用地，则断点前后的季节性成分的振幅会减小，即：
[0035][0036]
若公式(7)成立，则地表覆盖类型从耕地转变为了建筑用地。
[0037]
本发明与背景技术相比具有如下有益效果：
[0038]
本发明能够获得每个像素由耕地转变为建设用地的时间，并以此统计出每年由耕地转化为建设用地的面积。
附图说明
[0039]
图1是本发明实施例耕地损失评估流程图。
[0040]
图2是本发明实施例modis ndvi时序数据分解示意图。
[0041]
图3是本发明实施例山东省逐年耕地损失面积。
具体实施方式
[0042]
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
[0043]
具体来说，该方法包括以下步骤：
[0044]
(1)构造研究区modis ndvi时序数据。下载时间分辨率高于1个月的modis ndvi数据，分别转换投影至地理坐标系下，然后将所有modis ndvi数据按照日期叠加到一起生成时序数据；利用研究区矢量数据裁剪出研究区modis ndvi时序数据。
[0045]
(2)利用现有土地利用产品，提取出第一期土地覆盖类型为耕地，后一期土地覆盖类型为建设用地的像素；
[0046]
(3)对于(2)中提取出的像素，逐像素按照下式进行建模和拟合：
[0047]
y
v
＝t
v
s
v
r
v
ꢀꢀ
(1)
[0048]
其中，t
v
代表趋势性成分，由线性函数表达；s
v
代表季节性成分，由三角函数的傅里叶级数表达；r
v
代表残差。
[0049]
其中：
[0050][0051]
t
v
＝α βt
ꢀꢀ
(3)
[0052]
其中α、β为线性函数的系数，γ
k
和θ
k
为系数，f为每年的modis ndvi观测数量，k为确定三角函数周期的系数，t为自变量，代表第t个观测值。
[0053]
1)若地表覆盖不发生变化，则modis ndvi时序数据不存在断点，故：
[0054]
[0055]
拟合结果得到残差平方和rss1。
[0056]
2)若地表覆盖发生变化，根据一般情况，地表覆盖变化不会频繁发生，则modis ndvi时序数据中存在1个断点，假设断点位置为t
*
(f＜t
*
＜l
‑
f)(一般而言，f取值为一个周期；l为modis ndvi时序数据观测总数量)，故：
[0057][0058]
式中，f取值为一个周期，l为modis ndvi时序数据观测总数量，α1、β1、α2、β2分别为两个线性函数的系数，γ
i
、θ
i
、γ
j
、θ
j
为系数，i、j为断点前后确定三角函数周期的系数，t
*
为断点位置，t为自变量，代表第t个观测值；
[0059]
利用公式(5)遍历t
*
的可能值，求取残差平方和最小值rss2，并记录下此时t
*
的取值为
[0060]
(4)利用aic准则判断modisndvi时序中是否存在断点。
[0061]
aic准则表达式如下：
[0062]
aic＝2l lln(rss/l)
ꢀꢀ
(6)
[0063]
其中，l为公式(4)(modis ndvi时序数据中不存在断点)或公式(5)(modis ndvi时序中存在一个断点)的参数个数。
[0064]
将rss1及公式(4)的参数个数代入公式(6)得到aic1，将rss2及公式(5)的参数个数代入公式(6)得到aic2；若aic1<aic2，则该时序数据中不存在断点，反之，若aic1>aic2，则该时序数据中存在断点，断点时刻为该时刻为地表覆盖类型发生变化的时刻。
[0065]
确认地表覆盖变化是否属于从耕地转变为建筑物。若地表覆盖类型从耕地转变为建筑用地，则断点前后的季节性成分的振幅会减小，即：
[0066][0067]
若地表覆盖变化从耕地转变为建筑，则掩模图设置为1，否则设置为0。
[0068]
下面是一个更具体的例子：
[0069]
如图1所示，一种基于modis时序数据的耕地损失评估模型包括以下步骤：
[0070]
1、构造研究区modisndvi时序数据
[0071]
下载2001年1月
‑
2009年12月行列号为h27v05的月合成modisndvi数据，每年观测数量为23，共有207个观测值，空间分辨率为250米；将投影转换到地理坐标系下，然后将其按照时间叠加起来，再按照山东省行政区划进行裁剪。
[0072]
2、利用globeland30
‑
2000和globeland30
‑
2010产品，提取出2000年土地覆盖类型为耕地，但2010年土地覆盖类型为建设用地的像素，并构建土地利用变化掩模图，使用众数聚合法将掩模图重采样到250米。
[0073]
3、逐像素按照公式(1)进行建模和拟合，其中f＝23。图2给出某个像素的拟合结果。
[0074]
4、对每个像素的建模和拟合过程，利用aic准则判断时序数据中是否存在断点，确
定最终的拟合方式。图2中，利用aic准则判断为断点存在，且断点位置为第94个观测值，即为2005年2月份。
[0075]
5、对每个地表覆盖发生变化的像素，判断断点前后的季节性成分振幅变化，若振幅减小，则确认像素的地表覆盖由耕地转变为建设用地。图2中的季节性成分振幅在地表覆盖发生变化后减小了，则可确定为由耕地转变为建设用地；即该像素在2005年2月份由耕地转变为建设用地。
[0076]
6、按照年份统计山东省每年由耕地转变为建设用地的面积，如图3。
[0077]
总之，本发明利用modisndvi时序数据提取了由耕地转变为建设用地的像素，并获得了由耕地转变为建设用地的时间。
[0078]
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。