本发明涉及一种遥感数字图像处理技术，尤其是涉及一种变化检测方法。
背景技术：
通过对地表进行重复观测所获得的多期遥感数据，快速、准确地检测出地表覆盖发生变化的区域，是遥感领域一直以来着力解决的技术问题(参见文献：HANSEN M C,LOVELAND T R.A review of large area monitoring of land cover change using Landsat data[J].Remote Sensing of Environment,2012,122(66-74))。随着国内外光学遥感卫星的相继投入使用，已经基本具备了稳定的对地重复观测能力，形成了多个时期的多光谱遥感影像数据集和归一化差异植被指数(Normalized Difference Vegetation Index，NDVI)时间序列数据集。然而，由于遥感多光谱传感器对地表成像过程中的大气条件等随机噪声会使地表覆盖的光谱量值产生误差，当较多量值均受误差影响时光谱就会带有噪声。同时，“同物异谱，异物同谱”现象会导致同类地表覆盖的光谱具有较大系统性偏差，而各类地表覆盖光谱间的区分度较差。受到物候、季相和地表粗糙度等因素影响，NDVI时间序列在不同年份也常常具有一定系统性波动差异。因此，当前的地表覆盖变化检测项目是在这种数据背景下，实施对地表覆盖发生变化的区域进行检测的。
目前已经发展出了多种变化检测方法，用于变化检测过程。影像差值法和比值法(参见文献：LU D,MAUSEL P,BROND ZIO E,et al.Change detection techniques[J].International Journal of Remote Sensing,2004,25(12):2365-401.)利用地表覆盖在两期影像光谱间的量值差异或者相似度作为判别发生变化与否的依据。它们的检测效率高，但容易受部分光谱值噪声的干扰，使得检测精度不高。变化向量分析法(参见文献：MALILA W A.Change vector analysis an approach for detecting forest changes with Landsat[M].LARS Symposia.1980:385.)将光谱的多个量值组成向量，以两个向量的绝对偏差作为判别变化的依据。该方法应用广泛，但难以克服地表覆盖光谱内源噪声的不利影响。光谱角制图法(参见文献：KRUSE F A,LEFKOFF A B,BOARDMAN J W,et al.The spectral image processing system(SIPS)—interactive visualization and analysis of imaging spectrometer data[J].Remote sensing of environment,1993,44(2-3):145-63.)利用两期光谱向量间的夹角作为判别发生变化的依据，在遇到光谱内源噪声时仍然难以克服光谱内源噪声的影响。相关系数法(参见文献：周启明.多时相遥感影像变化检测综述[J].地理信息世界,2011,9(2):28-33.)利用两期光谱间的相关性确定地表覆盖的变化强度，当光谱间存在系统误差时其判别依据得出的变化和未变化区域区分度较差，检测精度较低。光谱梯度差异法(参见文献：陈军,陆苗,陈利军.一种基于光谱斜率差异检测地表覆盖变化的方法[P].北京：CN103049916A,2013-04-17.)利用两期光谱斜率间的差异来检测地表覆盖变化，避免了光谱内源噪声对变化检测的影响，提高了检测精度。然而，该方法难以克服“同物异谱”对检测结果的不利影响。NDVI变化向量分析法(参见文献：LAMBIN E F,STRAHLER A H.Change-vector analysis in multitemporal space-a tool to detect and categorize land-cover change processes using high temporal-resolution satellite data[J].Remote Sensing of Environment,1994,48(2):231-44.陈云浩,李晓兵,陈晋,史培军.1983—1992年中国陆地植被NDVI演变特征的变化矢量分析[J].遥感学报,2002(01):12-18+81.)利用两个年际的NDVI时间序列组成向量，以两个向量间NDVI量值的总体差异来确定变化强度，进而检测地表变化。该方法容易受到两个NDVI时间序列系统性差异的影响，依据其得出的变化强度不能准确地判断变化。NDVI梯度差异检测法(参见文献：Lu,M.；Chen,J.；Tang,H.；Rao,Y.；Yang,P.；Wu,W.Land cover change detection by integrating object-based data blending model of landsat and modis[J].Remote Sensing of Environment 2016,184,374-386.)利用两个年际NDVI时间序列曲线的梯度差异，作为计算地表覆盖变化强度的依据，提高了无NDVI时间序列系统差异情况下的地表覆盖检测精度。基于兰氏距离NDVI的变化检测法(参见文献：赵威成,祁向前,叶欣.基于兰氏距离NDVI的土地覆盖变化检测[J].安徽农业科学,2017,45(06):221-223.)在两个年份NDVI量值的基础上，以兰氏距离度量地表覆盖的变化强度。因此，该方法仍然难以避免两个年际NDVI系统偏差的影响。李月臣等(参见文献：李月臣,陈晋,宫鹏,岳天祥.基于NDVI时间序列数据的土地覆盖变化检测指标设计[J].应用基础与工程科学学报,2005(03):44-58.)设计了表征NDVI时间序列曲线数值差异的值指数和形状差异的形指数。该值指数是采用兰氏距离公式，分别计算各个时间节点上NDVI的差与和的比值，并累加这些比值得出；形指数是依据由每个时间节点计算出的交叉相关系数公式得出。其地物变化强度是由值指数和形指数的算数平均值得出。该方法由于采用了各时间节点上的NDVI量值计算变化强度，因此难以克服两个年际NDVI时间序列系统差异的不利影响。殷守敬(参见文献：殷守敬.基于时序NDVI的土地覆盖变化检测方法研究[D].武汉大学,2010.)在李月臣的基础上，通过对时间节点增加时间偏移量，计算出新的形指数，并将原值指数和新的形指数的加权平均值作为计算得出的变化强度。该方法仍然难以克服两个年际NDVI时间序列系统差异的不利影响，变化强度上的变化和未变化的区域对比度较低，检测精度有待提高。
总体上，现有基于多光谱影像或者NDVI时间序列的变化检测方法，还不能克服两个时期的光谱间或者两个年际的NDVI时间序列间系统性差异的不利影响，检测精度不高，难以对地表变化实现准确地发现。造成上述问题的原因在于，现有变化检测方法从光谱或者NDVI时间序列的量值、比值、角度、斜率、相关性等方面确定地表的变化强度，而光谱或者NDVI时间序列间的系统差异会影响它们的值，使得变化和未变化地表区域的变化强度值较接近，进而难以准确检测出变化。
技术实现要素：
本发明的目的在于针对上述现有技术的不足，提供一种基于多光谱影像和NDVI时间序列的地表覆盖变化检测方法。
本发明的思想是：以光谱相关系数减弱光谱间系统差异的影响，以NDVI时间序列曲线的相角累积量、基线累积量、相对累积率和相对基线过零率四种形状参数减弱NDVI时序曲线间系统差异的影响；结合光谱相关系数和四种形状参数共同增大变化与未变化地表区域变化强度的对比度，并根据阈值判断出地表变化。
基于多光谱影像和NDVI时间序列的地表覆盖变化检测方法，包括以下步骤：
a、叠加两个时期的多光谱影像，对每个像元分别计算两个光谱间的皮尔森相关系数,并由该系数确定出光谱变化强度；
b、叠加两个时期上各年际的NDVI时间序列数据，对每个像元分别计算各曲线的相角累积量、基线累积量、相对累积率和相对基线过零率四种形状参数值；由形状参数的差异确定出相应的时序变化强度；
c、将每个像元的光谱变化强度和时序变化强度分别进行归一化处理，得到相应的归一化后的相对变化强度，并加权集成为一种光谱时序变化强度；
d、设定变化阈值，将光谱时序变化强度大于该阈值的像元划分为变化像元，将其余像元划分为未变化像元，由全部的变化像元组成地表覆盖变化区域。
步骤a所述的光谱变化强度由公式
MC＝1-p (1)
计算得出，其中p为皮尔森相关系数。
步骤b所述的两个时期上各年际的NDVI时间序列数据，是空间分辨率与步骤a中的多光谱影像相同，时间分辨率不低于30天的NDVI数据；相角累积量由公式(2)计算，
其中n为相角个数，V(i,s)和V(i,e)分别表示时间序列曲线上第i个NDVI值，t(i,s)表示相角的开始季相时间，t(i,e)表示相角的结束季相时间；
基线累积量由公式(3)计算，
其中ΔVi表示第i个NDVI曲线值和基线值间的差异，B表示基线区间内的所有NDVI集合，基线为NDVI曲线上前三个值的平均值和后三个值的平均值的连线，由平面两点坐标公式可求出该连线上的NDVI值，进而计算出ΔVi；
相对累积率由公式(4)计算，
其中V0表示生长季开始的NDVI值，取前三个NDVI值的平均值，Vi表示NDVI曲线上的第i个NDVI值；
相对基线过零率由公式(5)和由公式(6)计算，
其中T是总体时间间隔，V是NDVI值，Vm表示该时间段内所有NDVI的平均值；
四种形状参数的时序变化强度由公式(7)计算，
MΔf＝|Δf|p (7)
其中Δf表示相同的两个形状参数的值的偏差，对于相对累积率，p为2，对于相角累积量，基线累积量和相对基线过零率，p为1；
步骤c所述的归一化处理由公式(8)计算，
其中M是变化强度，Mmin和Mmax分别是变化强度的最小值和最大值；
加权集成由公式(9)实现，
w是加权因子，I是相对变化强度，i为0时表示光谱相对变化强度，i从1到4依次表示四种时序相对变化强度。
步骤d所述的将光谱时序变化强度大于该阈值的像元划分为变化像元，是通过运算操作中的大于运算符实现的。
所述的设定变化阈值，可以是人工根据变化强度设定，也可以是由计算机自动判别设定；运算操作中的大于运算符，实现将光谱时序变化强度大于该阈值的像元划分为变化像元。
有益效果：
本发明与现有技术共有的技术特征是：
均是使用两期多光谱影像或者两个年际的NDVI时间序列数据进行变化强度计算，计算结果均是形成一种变化强度；均利用了统计学中常用的皮尔森相关系数作为确定变化强度的依据之一；均采用了设定阈值的方式，将地表划分为变化区域和无变化区域。
本发明与现有技术的区别特征是：
本发明同时采用多光谱影像和NDVI时序曲线进行地表覆盖变化检测，特别是通过两个年际的NDVI时间序列曲线的相角累积量、基线累积量、相对累积率和相对基线过零率四种形状参数值间的差异作为时序变化强度的依据。
本发明通过归一化处理统一不同变化强度的量纲，并结合光谱变化强度和时序变化强度增大变化和未变化地表区域的对比度。
本发明克服了两个时期的光谱间以及两个年际的NDVI时间序列间系统性差异的不利影响，以较高的精度实现对地表覆盖变化的准确发现。同时，本发明所依据的数学基础清晰、检测步骤简洁，便于地表覆盖变化检测的自动化开展。
附图说明
附图1基于多光谱影像和NDVI时间序列的地表覆盖变化检测方法原理，
其中R1和R2分别指时期1和时期2的多光谱影像，T1和T2分别指时期1和时期2的NDVI时间序列数据，c指皮尔森相关系数，Mc指光谱变化强度；f1到f4顺次分别指T1的四种形状参数值，f’1到f’4顺次分别指T2的四种形状参数值，M1到M4顺次分别指四种形状参数的时序变化强度；I0指相对光谱变化强度，I1到I4顺次分别指四种相对时序变化强度，M指光谱时序变化强度，H指设定的变化阈值，A指变化地表，B指未变化地表
附图2NDVI时间序列曲线的相角示意图
附图3NDVI时间序列曲线的基线累积示意图
附图4NDVI时间序列曲线计算得出的相对累积曲线示意图
附图5NDVI时间序列曲线的相对基线过零点示意图
附图6本发明的地表变化检测实施例示意图
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例对本发明作进一步详细说明。
本实施例中所使用的多光谱影像为美国地质调查局网站上公开的，空间分辨率为30米的两景Landsat影像如附图6-1和附图6-2。经过ENVI软件中的辐射定标和大气校正功能，将影像的DN值转换为地表反射率，并将两幅影像进行了配准。其T1年的光谱如附图6-5，其T2年的光谱如附图6-6。
实施例中所使用的NDVI时间序列数据为由美国国家航空航天局网站下载的MOD13Q1产品，经过降尺度处理得到的空间分辨率为30米，时间分辨率为16天的NDVI时序数据(如附图6-3和6-4)，其T1年的NDVI时序曲线如附图6-7，其T2年的NDVI时序曲线如附图6-8。
假设需要检测某地区在T1年到T2年之间发生的变化，则需要经过以下步骤：
a、叠加两幅多光谱影像，由ENVI软件的波段运算功能，获得每个像元上两个光谱间的皮尔森相关系数(如附图6-9),并由公式
MC＝1-p (1)
计算得出光谱变化强度Mc；
b、叠加T1年和T2年的NDVI时间序列数据，由ENVI软件波段运算功能，对每个像元的各个曲线，由公式
计算曲线的相角累积量f1(如附图6-10)，由公式
计算曲线的基线累积量f2(如附图6-11)，由公式
计算曲线的相对累积率f3(如附图6-12)，由公式
计算曲线的相对基线过零率f4(如附图6-13)。根据T1年和T2年的NDVI时序曲线的四种形状参数值f1、f2、f3和f4，由公式
MΔf＝|Δf|p (7)
计算四种形状参数的时序变化强度，得到与形状参数值顺次相对应的M1、M2、M3和M4；
c、将光谱变化强度Mc和四种形状参数的时序变化强度M1、M2、M3和M4，根据由公式
分别进行归一化处理，得到顺次相对应的归一化后的相对变化强度I0(如附图6-14)、I1(如附图6-15)、I2(如附图6-16)、I3(如附图6-17)和I4(如附图6-18)。取加权因子均为1，将光谱相对变化强度I0和四种时序相对变化强度I1、I2、I3和I4由公式
加权集成处理，得出最终的光谱时序变化强度(如附图6-19)；
d、通过目视判别光谱时序变化强度，人工设定变化阈值为常数H，利用ENVI软件的波段运算功能，将大于H的像元判别为变化像元，标记为1，将小于或者等于H的像元判别为未变化像元，标记为0；将全部变化像元作为地表变化区域(如附图6-20中白色区域)。