CN103049916B - 一种基于光谱斜率差异检测地表覆盖变化的方法 - Google Patents

Abstract

一种基于光谱斜率差异检测地表覆盖变化的方法，所述方法包括如下步骤，A、将同一地点的Landsat图像的六个波段的光谱数据依据波谱连接成五个光谱段，分别针对T1、T2时刻的Landsat图像的每个点位，计算每个光谱段的斜率，B、分别计算T1、T2时刻的Landsat图像的每个点位的各波段的斜率，形成斜率向量，并通过斜率向量的绝对距离计算斜率变化强度，C、根据步骤B中计算出的斜率变化强度计算结果通过设置阈值确定所述Landsat图像中的地表覆盖变化区域和地表覆盖不变区域。本发明所提供的一种基于光谱斜率差异检测地表覆盖变化的方法，解决了传统的变化检测方法很难控制变化检测中内源噪声的影响的问题。

Description

一种基于光谱斜率差异检测地表覆盖变化的方法

技术领域

本发明涉及一种对遥感影像进行处理的方法，特别是一种通过比对不同时期的Landsat图像来测定地表覆盖变化的方法。

背景技术

一直以来，在遥感成像领域，Landsat TM/ETM+图像有着大量的应用，美国陆地卫星历经三代，第一代是Landsat-1卫星、Landsat-2卫星和Landsat-3卫星，第二代是Landsat-4卫星和Landsat-5卫星，第三代是Landsat-6卫星和Landsat-7卫星。TM是指美国陆地卫星4～5号专题制图仪Thematic Mapper所获取的多波段扫描影像，TM图像有7个波段，其波谱范围是：TM-1为0.45～0.52微米，TM-2为0.52～0.60微米，TM-3为0.63～0.69微米，以上为可见光波段；TM-4为0.76～0.90微米，为近红外波段；TM-5为1.55～1.75微米，TM-7为2.08～2.35微米，为中红外波段；TM-6为10.40～12.50微米，为热红外波段。因TM影像具较高空间分辨率、波谱分辨率、极为丰富的信息量和较高定位精度，从20世纪80年代中后期开始成为世界各国广泛应用的重要的地球资源与环境遥感数据源。

LANDSAT-7卫星的改进型主题测绘仪Enhanced Thematic Mapper Plus，也就是ETM+，是在TM基础上改进的。ETM+相对TM增加了1个金色谱段和2个增益区域，增加了太阳定标器，并提高了红外谱段的分辨率。

在对遥感成像图片结果的分析工作中，科学准确地测定全球地表覆盖的动态变化，对于研究地球系统的能量平衡、碳循环及其他生物地球化学循环、气候变化等有着十分重要的意义。遥感影像变化检测就是是利用不同时期的遥感影像分析并确定发生变化的区域及变化类别的处理过程，及时、准确的变化信息对地表覆盖制图、资源管理、应急救灾和动态信息服务等具有重要的作用。

地表覆盖(Land Cover)是指地球表面各种物质类型及其自然属性与特征的综合体(参见Chen,J.,Chen,X.,Cui,X.,Chen,J.,2011.Change vectoranalysis in posterior space:a new method for land cover change detection.IEEEGeoscience and remote sensing letters 8(2),317-321.)。

地表覆盖变化检测的基本前提是由地表覆盖变化引起的有意义的变化要大于其它因素所引起的次要变化，否则次要的变化会干扰和影响有意义的变化(参见Singh,A.,1989.Digital change detection techniques usingremotely-sensed data.International Journal of Remote Sensing 10(6),989–1003.)。影像地表覆盖的变化检测的其它因素包括大气状况、太阳高度角和土壤湿度等，可将这些因素分为两大类，即外源因素和内源因素。外源因素是由大气状况、太阳高度角等外在条件引起的，选取同源同时期的影像并通过预处理能降低外源因素的影响。内源因素是同种地物内部构造在不同时期发生的变化，如不同时期土壤湿度的不同、水体浑浊度的不同和建成区新旧程度的不同等，这些因素并不能通过数据选取和预处理解决，因此要求变化检测模型中能够对内在因素的干扰具有鲁棒性。

目前有多种变化检测的算法，如变化向量分析、差值法、比值法、相关系数法和夹角余弦法等(参见Coppin,P.,Jonckheere,I.,Nackaerts,K.,Muys,B.,2004.Digital change detection methods in ecosystem monitoring:a review.International Journal of Remote Sensing 25(9),1565-1596.,Lu,D.,Mausel,P.,Brondizio,E.,Moran,E.,2004.Change detection techniques.International Journalof Remote Sensing 25(12),2365-2407.)，这些算法的特点是从光谱量值方面计算变化强度。内源噪声通常影响的是光谱曲线的量值，传统的变化检测方法在光谱空间计算光谱值变化强度的大小，这样就很难区分真变化和伪变化。

发明内容

本发明提供了一种基于光谱斜率差异检测地表覆盖变化的方法，该方法可减少或避免前面所提到的问题。

为解决上述问题，本发明提出了一种基于光谱斜率差异检测地表覆盖变化的方法，所述方法用于对同一地点的两个不同时刻的Landsat图像进行地表覆盖变化检测，从而测定所述地点发生地表覆盖变化的区域范围，所述两个不同时刻分别为T1、T2时刻，所述方法包括如下步骤，

A、将同一地点的Landsat图像的六个波段的光谱数据依据波谱连接成五个光谱段，分别针对所述T1、T2时刻的Landsat图像的每个点位，计算每个光谱段的斜率，光谱段(n,n+1)的斜率计算公式如下：

$k_{(n,n+1)}=\frac{\mathrm{\Δref}}{\mathrm{\Δ\λ}}=\frac{{\mathrm{ref}}_{n+1}-{\mathrm{ref}}_n}{{{\mathrm{\λ}}^{\′}}_{n+1}-{{\mathrm{\λ}}^{\′}}_n}$

其中，Δref是波段n+1和波段n的反射率差值，Δλ是波段n+1和波段n的波段长度之差，λ'_n是归一化后无量纲的第n波段波长，λ'_n的计算公式如下：

${\mathrm{\λ}}_n^{\′}=\frac{{\mathrm{\λ}}_n-\min \left(\mathrm{\λ}\right)}{\max \left(\mathrm{\λ}\right)-\min \left(\mathrm{\λ}\right)}$

其中，λ_n是第n波段的波长，min(λ)是所述Landsat图像的所有波段中波长的最小值，max(λ)是所述Landsat图像的所有波段中波长的最大值，λ_n、min(λ)、max(λ)的单位均为微米；

B、根据步骤A中的斜率计算公式，分别计算所述T1、T2时刻的Landsat图像的每个点位的各波段的斜率，形成斜率向量，以此表示光谱的形状，T1时刻的斜率向量是K_i＝(k_i(1,2),k_i(2,3),...,k_i(n,n+1))^T，T2时刻的斜率向量是K_j＝(k_j(1,2),k_j(2,3),…,k_j(n,n+1))^T，

通过斜率向量的绝对距离计算斜率变化强度，计算公式如下：

$\begin{array}{c}{d}_K\underset{n=1}{\overset{5}{\mathrm{\Σ}}}\left|\mathrm{\Δ}{k}_{(n,n+1)}\right|=\underset{n=1}{\overset{5}{\mathrm{\Σ}}}|k_{i(n,n+1)}-k_{j(n,n+1)}|=|\frac{{\mathrm{ref}}_{i,n+1}-{\mathrm{ref}}_{i,n}}{{{\mathrm{\λ}}^{\′}}_{n+1}-{{\mathrm{\λ}}^{\′}}_n}-\frac{{\mathrm{ref}}_{j,n+1}-{\mathrm{ref}}_{j,n}}{{{\mathrm{\λ}}^{\′}}_{n+1}-{{\mathrm{\λ}}^{\′}}_n}\\ =\underset{n=1}{\overset{5}{\mathrm{\Σ}}}|\frac{({\mathrm{ref}}_{i,n+1}-{\mathrm{ref}}_{j,n+1})-({\mathrm{ref}}_{i,n}-{\mathrm{ref}}_{j,n})}{{{\mathrm{\λ}}^{\′}}_{n+1}-{{\mathrm{\λ}}^{\′}}_n}\end{array}$

其中，ref_i,n+1和ref_i,n是T1时刻第n+1波段和第n波段的反射率，ref_j,n+1和ref_j,n是T2时刻第n+1波段和第n波段的反射率，对应波段的归一化波长分别是λ'_n+1和λ'_n；d_K是两斜率向量的变化强度；

C、根据步骤B中计算出的所述T1、T2时刻的Landsat图像的每个点位的斜率变化强度计算结果d_K，通过设置阈值确定所述Landsat图像中的地表覆盖变化区域和地表覆盖不变区域。

优选地，所述方法进一步包括如下步骤D，用于确定地表覆盖变化的区域的地物变化类别，所述步骤D为，根据已知的各种地物特征点及其所对应变化类型的斜率变化特征数据，利用地表覆盖变化区域的像素和上述已知的变化类型斜率变化特征数据的匹配，确定出所述像素的变化类型，从而确定出步骤C中所述Landsat图像中的地表覆盖变化区域的每个像素的地物变化类型。

本发明所提供的一种基于光谱斜率差异检测地表覆盖变化的方法，其用于利用Landsat图像进行地表覆盖变化检测，通过光谱斜率的计算，将变化检测从传统的光谱空间转化到斜率空间，通过比较光谱斜率的强度差异判断地表覆盖是否发生变化，解决了传统的变化检测方法在光谱空间计算光谱值变化强度的大小，这样就很难区分真变化和伪变化，也即是很难控制变化检测中内源噪声的影响的问题。此外，本发明所提供的一种基于光谱斜率差异检测地表覆盖变化的方法，其利用光谱斜率走势的变化规律确定地表覆盖的变化类型，通过对斜率走势差异的分析，构建斜率变化特征数据模型，即可容易的进一步确定地表覆盖的变化类型。

附图说明

以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范围。其中，

图1a为根据本发明的一个具体实施例的一种基于光谱斜率差异检测地表覆盖变化的方法的一个地点的T1时刻的Landsat影像示意图；

图1b为图1a所示的一个地点的T2时刻的Landsat影像示意图；

图1c为图1a、1b所示的点位1的光谱曲线变化示意图；

图1d为图1a、1b所示的点位2的光谱曲线变化示意图；

图1e为图1c、1d所示的点位1、点位2的光谱量值变化之和的示意图图；

图1f为图1a、1b所示的点位1的各光谱段斜率变化示意图；

图1g为图1a、1b所示的点位2的各光谱段斜率变化示意图；

图1h为图1c、1d所示的点位1、点位2的光谱段斜率差异之和的示意图；

图2是根据本发明的一个具体实施例的一种基于光谱斜率差异检测地表覆盖变化的方法的斜率变化特征链模型示意图；

图3a是根据本发明的一个具体实施例的陕西省礼泉县及其周边区域2000年的Landsat影像；

图3b是图3a所示区域的2009年的Landsat影像；

图3c是图3b所示Landsat影像的分类示意图；

图3d是图3a、图3b所示的Landsat影像的斜率变化强度影像示意图；

图3e是图3d所示的斜率变化强度影像的灰度直方图；

图3f为图根据图3d所示的Landsat影像的斜率变化强度影像所确定的变化区域示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。

下面参照附图详细说明根据本发明的一种基于光谱斜率差异检测地表覆盖变化的方法的实施步骤及其原理。所述方法用于对同一地点的两个不同时刻的Landsat图像进行地表覆盖变化检测，从而测定所述地点发生地表覆盖变化的区域范围，所述两个不同时刻分别为T1、T2时刻，所述方法包括如下步骤，

步骤A、将同一地点的Landsat图像的六个波段的光谱数据依据波谱连接成五个光谱段，分别针对所述T1、T2时刻的Landsat图像的每个点位，计算每个光谱段的斜率，光谱段(n,n+1)的斜率计算公式如下：

$k_{(n,n+1)}=\frac{\mathrm{\Δref}}{\mathrm{\Δ\λ}}=\frac{{\mathrm{ref}}_{n+1}-{\mathrm{ref}}_n}{{{\mathrm{\λ}}^{\′}}_{n+1}-{{\mathrm{\λ}}^{\′}}_n}$

其中，Δref是波段n+1和波段n的反射率差值，Δλ是波段n+1和波段n的波段长度之差，λ'_n是归一化后无量纲的第n波段波长，λ'_n的计算公式如下：

${\mathrm{\λ}}_n^{\′}=\frac{{\mathrm{\λ}}_n-\min \left(\mathrm{\λ}\right)}{\max \left(\mathrm{\λ}\right)-\min \left(\mathrm{\λ}\right)}$

其中，λ_n是第n波段的波长，min(λ)是所述Landsat图像的所有波段中波长的最小值，max(λ)是所述Landsat图像的所有波段中波长的最大值，λ_n、min(λ)、max(λ)的单位均为微米；

因为斜率是一个无量纲值，Δref是波段n+1和波段n的反射率差值，无量纲值；Δλ是波段n+1和波段n的波段长度之差，单位是微米，分子和分母的量纲不统一，所以需要对波长做无量纲处理。

光谱斜率k_(n,n+1)反映了光谱段(n,n+1)内光谱值的增减走势。如果k_(n,n+1)>0，表示从波段n到n+1光谱值增加；如果k_(n,n+1)<0，表示从波段n到n+1光谱值减少；如果k_(n,n+1)＝0，表示从波段n到n+1光谱值不变。|k_(n,n+1)|越大，表示光谱值的变化越大，反之则变化越小。

斜率差异包含两方面的内容：斜率变化强度差异和斜率走势差异。其中斜率变化强度差异是对两种地物光谱形状差异的总体量化，因此通过计算斜率变化强度可以很容易的判断地表覆盖是否发生变化。

图1a为根据本发明的一个具体实施例的一种基于光谱斜率差异检测地表覆盖变化的方法的一个地点的T1时刻的Landsat影像示意图；图1b为图1a所示的一个地点的T2时刻的Landsat影像示意图；图1c为图1a、1b所示的点位1的光谱曲线变化示意图；图1d为图1a、1b所示的点位2的光谱曲线变化示意图；图1e为图1c、1d所示的点位1、点位2的光谱量值变化之和的示意图图；图1f为图1a、1b所示的点位1的各光谱段斜率变化示意图；图1g为图1a、1b所示的点位2的各光谱段斜率变化示意图；图1h为图1c、1d所示的点位1、点位2的光谱段斜率差异之和的示意图；参照图1a-1h所示，如图1a、1b所示是土壤湿度不同对地物的影响。其中点位1是地物发生变化的区域，由裸地变化为植被，光谱曲线的形状变化较大，如图1c所示；点位2是不变区域，均为裸地，但是由于不同时期土壤湿度不同的影响，使得光谱曲线有所变化，如图1d所示，光谱曲线的整体形状不变，但是各波段的光谱值均有所变化。传统的变化检测方法在光谱空间计算光谱值变化强度的大小，图1a、1b中点位1和点位2的变化强度都较大，如图1e的累积柱状图所示，这样就很难区分真变化和伪变化，也就是说，传统的变化检测方法在判断区域是否发生地物类型变化时，存在很大的误差。

根据步骤A中的斜率计算方法，计算图1a、1b中两个点位不同时期光谱曲线的斜率，将其转化到斜率空间，如图1f、1g、1h所示。和图1c、1d、1e相比，地表真实变化的点位1在斜率空间光谱斜率的变化仍然很大，各个光谱段的斜率差异之和是6.925，而由于土壤湿度不同影响的点位2在斜率空间的变化差异量明显减小，其斜率差异之和是0.517，远小于点位1的斜率差异，因此变化和不变化的点位能明显的区分。由此可见，在斜率空间，不仅能有效的表示出地物的变化，也能有效规避噪声对变化检测的影响。

B、根据步骤A中的斜率计算公式，分别计算所述T1、T2时刻的Landsat图像的每个点位的各波段的斜率，形成斜率向量，以此表示光谱的形状，当地物发生变化时，光谱曲线的形状发生变化，因此通过计算光谱斜率向量的差异，即斜率变化强度，即可很容易的度量地物是否发生变化。

T1时刻的斜率向量是K_i＝(k_i(1,2),k_i(2,3),...,k_i(n,n+1))^T，T2时刻的斜率向量是K_j＝(k_j(1,2),k_j(2,3),…,k_j(n,n+1))^T，

通过斜率向量的绝对距离计算斜率变化强度，计算公式如下：

$\begin{array}{c}{d}_K\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}\left|\mathrm{\&Delta;}{k}_{(n,n+1)}\right|=\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}|k_{i(n,n+1)}-k_{j(n,n+1)}|=|\frac{{\mathrm{ref}}_{i,n+1}-{\mathrm{ref}}_{i,n}}{{{\mathrm{\&lambda;}}^{\&prime;}}_{n+1}-{{\mathrm{\&lambda;}}^{\&prime;}}_n}-\frac{{\mathrm{ref}}_{j,n+1}-{\mathrm{ref}}_{j,n}}{{{\mathrm{\&lambda;}}^{\&prime;}}_{n+1}-{{\mathrm{\&lambda;}}^{\&prime;}}_n}\\ =\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}|\frac{({\mathrm{ref}}_{i,n+1}-{\mathrm{ref}}_{j,n+1})-({\mathrm{ref}}_{i,n}-{\mathrm{ref}}_{j,n})}{{{\mathrm{\&lambda;}}^{\&prime;}}_{n+1}-{{\mathrm{\&lambda;}}^{\&prime;}}_n}\end{array}$

其中，ref_i,n+1和ref_i,n是T1时刻第n+1波段和第n波段的反射率，ref_j,n+1和ref_j,n是T2时刻第n+1波段和第n波段的反射率，对应波段的归一化波长分别是λ'_n+1和λ'_n；d_K是两斜率向量的变化强度；d_K越大表示不同时期曲线形状的整体差异越大，变化的可能性就越大，反之则越小。

C、根据步骤B中计算出的所述T1、T2时刻的Landsat图像的每个点位的斜率变化强度计算结果d_K，通过设置阈值确定所述Landsat图像中的地表覆盖变化区域和地表覆盖不变区域。

传统的变化阈值设置多基于经验或人工试错的方法，需要较长时间的摸索，主观性较强，效率不高。虽然目前也有多种自动阈值算法，如Maximumbetween-class variance method(参见Otsu,N.,1979.A Threshold SelectionMethod from Gray-Level Histograms.IEEE Transactions on Geoscience andRemote Sensing 9(1),62-66.)，Maximum Correlation method(参见Yen,J.,Change,F.,&Change,S.,1995.A New Criterion for Automatic MultilevelThresholding.IEEE Transaction on Geoscience and Remote Sensing 4(3),370-378.),Minimum error method(参见Kittler,J.,Illingworth,J.,1986.Minimum Error Thresholding.Pattern Recognition 19(1),41-47.)等，但是这些算法适用于直方图是双峰的灰度影像。在对Landsat图像的地表覆盖变化检测中，变化区域所占的比例可能较少，且变化区域的变化强度可能较大，这样就会集中在直方图的尾部，因此很难形成双峰，因此这些自动的阈值设置方法不适合用于确定Landsat图像的地表覆盖变化检测的变化阈值。

本发明中的阈值设置方法可以利用Chen等2003年提出的双窗口变步长阈值搜索法(Double-Windows Flexible Pace Search,DFPS)设置变化阈值(参见Chen,J.,Gong,P.,He,C.,Pu,R.,Shi,P.,2003.Land-use/land-cover changedetection using improved change-vector analysis.Photogrammetric Engineeringand Remote Sensing 69(4),369-379.)，该方法通过选取典型变化训练样区，在一定范围内利用变步长搜索能使检验成功率达到最大的阈值。

在步骤C中通过DFPS阈值设置法，即可将所述Landsat图像中的地表覆盖变化区域和地表覆盖不变区域区分开。

在一个优选实施例中，为了进一步确定所述Landsat图像中的地表覆盖变化类型，所述方法进一步包括步骤D，所述步骤D为，根据已知的各种地物特征点及其所对应变化类型的斜率变化特征数据，利用地表覆盖变化区域的像素和上述已知的变化类型斜率变化特征数据的比对，确定出所述像素的变化类型，从而确定出步骤C中所述Landsat图像中的地表覆盖变化区域的每个像素的地物变化类型。

目前有多种变化类型的确定方法，可分为两大类，一类是从不变区域上随机采集样本对变化区域进行分类(参见Xian,G.,Collin,H.,Fry,J.,2009.Updating the 2001National Land Cover Database land cover classification to2006by using Landsat imagery change detection methods.Remote Sensing ofEnvironment 113(6),1133-1147.)，另一类是采用变化向量分析法确定变化类别，常见的CVA类别提取算法如二维向量空间的三角函数法(参见Malila,W.A.,1980.Change vector analysis:an approach for detecting forest changes withLandsat,Proceedings of the 6th Annual Sym-posium on Machine Processing ofRemotely Sensed Data,03-06June,Purdue University,West Lafayette,Indiana,pp.326–335.)，高于二维空间的波段符号组合法(参见Michalek,J.L.,Wagner,T.W.,Luczkovich,J.J.,and Stoffle R.W.,1993.Multispectral change vectoranalysis for monitoring coastal marine environments,PhotogrammetricEngineering&Remote Sensing,59:381–384.；Jensen,J.,R.,2004.IntroductoryDigital Image Processing(3rd Edition).New Jersey:Prentice Hall press.)等。分类对训练样本的要求较高，随机选取的样本很难保证其可靠性和准确性，CVA的三角函数法和波段符号组合法适用于单一要素地物的变化检测(如森林的变化检测)或地物类较单一的区域。

借鉴Chen等提出的线目标间拓扑关系细化计算方法的思路(参见Chen,J.,Liu,W.,Li,Z.,Zhao,R.,Cheng,T.,2007,Detection of spatial ConflictsBetween Rivers and Contours in digital map updating,International Journal ofGIS,21(10),1093-1114.)，本发明在已有变化检测结果和参考图像分类的基础上，提出一种基于斜率变化特征链模型的变化类型确定方法。也即是基于已知的各种地物特征点及其所对应变化类型的斜率变化特征数据来对变化类型进行确定。

斜率变化特征数据，也即是斜率差异，包含两方面的内容：斜率变化强度差异和斜率走势差异。斜率变化强度是对两种地物光谱形状差异的总体量化，因此通过计算斜率变化强度判断地物是否发生变化。光谱走势包含了光谱斜率方向和斜率大小两方面的内容，通常一种地物变为另一种地物其光谱走势是有规律可循。因此，本发明从斜率方向和斜率大小两方面的变化分析光谱走势差异，可以建立基于已知的各种地物特征点及其所对应变化类型的斜率变化特征数据模型，也即是斜率变化特征链模型，进而确定变化类型。

假设已知各种地物特征点及其所对应变化类型的斜率变化特征链，则可以利用变化像素和变化类型斜率变化特征链的匹配，确定出变化像素的变化类型。因此，如何计算斜率变化特征链以及变化像素和特征点的匹配是该模型的关键。

对于需要进行地表覆盖变化检测的两期Landsat影像，和历史影像相比，当前的影像更容易获得地面真实的样本。因此，将当前影像作为参考影像并进行分类，结合基于步骤A、B得出的光谱的变化检测结果，确定出历史影像上不变区域的分类结果，然后求不变区域的每个类别的光谱求均值，就可以以此作为各地物类型的特征点。

图2是根据本发明的一个具体实施例的一种基于光谱斜率差异检测地表覆盖变化的方法的斜率变化特征链模型示意图。斜率变化特征链是对不同地物特征点的光谱变化细化描述与计算，具体来说，将整条光谱曲线分割为若干个局部的光谱段，在每个局部光谱段内从斜率走势变化和斜率大小变化两方面刻画光谱的不同，然后将其组合形成斜率变化特征链，如图2所示。假设图中植被特征点的光谱向量是(x₁，x₂，x₃，…，x_n)，裸地特征点的光谱向量是(x₁'，x₂'，x₃'，…，x_n')，按照相邻的光谱将其分为5个光谱段，对于每个光谱段，计算其斜率方向和斜率大小的变化。斜率值的正负反映了光谱段的方向，根据光谱斜率值的正负定性地将斜率方向定义为正向、负向和水平向。正向斜率表示在光谱段(n-1,n)内光谱处于上升趋势，负向斜率表示在光谱段(n-1,n)内光谱处于下降趋势，水平向斜率表示在光谱段(n-1,n)内光谱不发生变化。用t_vb(n-1,n)表示植被和裸地在光谱段(n-1,n)内斜率方向的变化，如下列公式所示，

$t_{\mathrm{vb}(n-1,n)}=\left\{\begin{array}{cccc}1,& \mathrm{if}{x}_{(n-1,n)}\&times;{x^{\&prime;}}_{(n-1,n)}>0& \mathrm{or}& x_{(n-1,n)}={x^{\&prime;}}_{(n-1,n)}=0\\ -1,& \mathrm{if}{x}_{(n-1,n)}\&times;{x^{\&prime;}}_{(n-1,n)}\&le;0& \mathrm{and}& x_{(n-1,n)}\&NotEqual;{x^{\&prime;}}_{(n-1,n)}\end{array}\right.$

如果两种地物的斜率方向一致，则t_vb(n-1,n)的值为1；如果两种地物的斜率方向相反或仅有一个斜率是水平向，则t_vb(n-1,n)的值为-1。

在斜率方向不变的情况下，斜率的大小也存在不同，如光谱段(1,2)(2,3)和(5,6)所示，因此需要进一步细化光谱斜率的差异。定义c_vb(n-1,n)描述植被和裸地两类地物在光谱段(n-1，n)内斜率大小的变化：

c_vb(n-1,n)＝|x_(n-1,n)-x′_(n-1,n)|

按照光谱段的顺序，将斜率方向和大小的差异组成斜率变化特征链：

f_vb(t_vb,c_vb)＝t_vb(1,2)c_vb(1,2)t_vb(2,3)c_vb(2,3)…t_vb(n-_1,n)c_vb(n-1,n)

该特征链不仅反映了植被和裸地局部光谱段斜率的变化，同时也反映了两种地物整体光谱的变化。

按照上述方法计算各种变化类型的斜率变化特征链。同样，对于变化像素，也可以按照图2所示的方法计算斜率变化特征链：

f′(t′,c′)＝t′_(1,2)c′_(1,2)t′_(2,3)c′_(2,3)…t′_(n-1,n)c′_(n-1,n)

根据变化像素在参考影像上的类别，选取对应的变化类型特征链向量。假设变化像素在参考影像上属于第i类，则选取i类与其它类型发生变化时的斜率变化特征链和变化像素做匹配。首先，评估每个变化像素的链单元和变化类型的差异度，形成匹配链，然后计算总体匹配系数，以此作为判断变化像素类型的依据。链单元匹配计算是求各变化像素链单元和变化类型链单元的绝对差，如下公式所示：

n_ij(n-1,n)＝|t_ij(n-1,n)-t′_(n-1,n)|，h_ij(n-1,n)＝|c_(n-1),n-c′_(n-1),n|

其中n_ij(n-1,2)是斜率走势变化t'_(n-1,n)和t_ij(n-1,n)的匹配值，h_ij(n-1,2)是斜率走势变化c'_(n-1,n)和c_ij(n-1,n)的匹配值。然后将其按照光谱段的顺序形成匹配链(如公式(11)所示)，通过匹配链的匹配值可看出变化像素每个光谱链斜率方向和斜率大小与各类型光谱特征链的匹配程度，各光谱段的匹配值越小，说明越接近。

m_ij(n_ij,h_ij)＝n_ij(1,2)h_ij(1,2)n_ij(2,3)h_ij(2,3)…n_ij(n-1,n)h_ij(n-1,n)

根据匹配链计算总体匹配系数，如下公式所示，总体匹配值越小，表明和该变化类型越相似，最后取总体匹配系数的最小值所对应的类型作为变化像素的变化类型。

M_ij＝∑n_ij×∑h_ij

这样就可以得出Landsat影像中变化像素的变化类型，也就得到了Landsat图像中的地表覆盖变化区域的变化类型。

在另一个具体实施例中，选取陕西省礼泉县及其周边为研究区域，该区域位于陕西省关中西北部(34°29′N，108°25′E)，地处黄土高原与渭河谷地交汇地带。图3a是该区域2000年的Landsat影像；图3b是图3a所示区域的2009年的Landsat影像；近十年间该区域经济快速发展，人口膨胀和城市扩张等因素引起了地表覆盖的变化。利用支持向量机(Support Vector Machine，SVM)对2009年影像进行分类，然后参考高分影像和其它辅助资料对分类结果做人工修改，分类结果如图3c所示，该区域地表覆盖共分五类，即农田，建成区，水体，湿地和裸地。

在做地表覆盖变化检测之前，为了消除不同时刻影像的非地表显著变化，可以对数据进行预处理。结合原始影像，本发明可以选用的预处理工作包括几何纠正和大气纠正。也即是首先对图3a、3b所示的两期影像做几何纠正，精度达到半个像素以内，然后利用ATCOR2(A Spatially-AdaptiveFastAtmospheric Correction-2)进行大气纠正，将影像的DN值转化为地表反射率以消除大气状况不同对变化检测的影响。

在对图3a、图3b所示的Landsat影像进行预处理后，可根据本发明所提供的步骤A、步骤B的方法计算图3a、图3b所示的Landsat影像的斜率，并将变化检测从光谱空间转化到斜率空间，然后计算两期影像的斜率变化强度，图3d是图3a、图3b所示的Landsat影像的斜率变化强度影像示意图；参照图3d可知，图3d是斜率变化强度影像，图中灰度越亮的区域表示变化的可能性越大，反之则越小。这样就可以进一步测定Landsat图像中的地表覆盖变化和地表覆盖不变区域。

根据本发明所提供的步骤C，要对上述从步骤B计算出的Landsat图像的每个点位的斜率变化强度计算结果进一步设置阈值，图3e是图3d所示的斜率变化强度影像的灰度直方图；参照图3e可知，图3e是斜率变化强度影像的灰度直方图，呈单峰分布，因此不适合用Otsu等自动阈值的方法。在本具体实施例中，可选取典型变化区域训练样本，利用DFPS算法设置阈值，这样就可得到该区域的变化阈值是4.369，光谱差异度大于该阈值的即为变化区域，否则不变。图3f为图根据图3d所示的Landsat影像的斜率变化强度影像所确定的变化区域示意图；阈值确定后，变化区域的确定结果如图3f所示，从图中可看出变化的区域主要集中在城区周边，和实际情况相符合。

为了确定地表覆盖变化区域的变化类型，根据本发明所提供的步骤D，可根据已知的各种地物特征点及其所对应变化类型的斜率变化特征链，利用变化像素和变化类型斜率变化特征链的匹配，确定出变化像素的变化类型，从而确定出步骤C中所述Landsat图像中的地表覆盖变化区域的变化类型。

在本具体实施例中，对于不变区域，2000年影像的分类结果和2009年保持不变。结合2000年原始影像，即可计算每个类别的特征点，然后分析变化类别的斜率差异并构建斜率变化特征链模型。通过变化区域和斜率变化特征链的匹配，即可确定变化区域的变化类型。

为检测在本具体实施例中根据本发明提出了一种基于光谱斜率差异检测地表覆盖变化的方法的变化检测效果，分别对变化区域确定和变化类型分类两部分做精度评价。对于变化区域确定，通过采用随机抽样的方法选取样本点，然后辅助参考资料确定样本点的变化和不变，其中变化区域检验样本842个，不变区域检验样本1537个。经检验，变化检测结果的总体精度是96.595％，Kappa系数是0.924；

对于变化类型分类确定，通过采用随机抽样和误差矩阵的方法评价分类精度，共选取1500个覆盖各类型的检验样本，变化类型分类确定检测结果的总体精度是89.266％，Kappa系数是0.858。

由上述检测效果可见，本发明所提供的一种基于光谱斜率差异检测地表覆盖变化的方法，在地表覆盖变化检测的准确率上有着很高的精度。

本发明所提供的一种基于光谱斜率差异检测地表覆盖变化的方法，其用于利用Landsat图像进行地表覆盖变化检测，通过光谱斜率的计算，将变化检测从传统的光谱空间转化到斜率空间，通过比较光谱斜率的强度差异判断地表覆盖是否发生变化，解决了传统的变化检测方法在光谱空间计算光谱值变化强度的大小，这样就很难区分真变化和伪变化，也即是很难控制变化检测中内源噪声的影响的问题。此外，本发明所提供的一种基于光谱斜率差异检测地表覆盖变化的方法，其利用光谱斜率走势的变化规律确定地表覆盖的变化类型，通过对斜率走势差异的分析，构建斜率变化特征链模型，即可容易的进一步确定地表覆盖的变化类型。

本领域技术人员应当理解，虽然本发明是按照多个实施例的方式进行描述的，但是并非每个实施例仅包含一个独立的技术方案。说明书中如此叙述仅仅是为了清楚起见，本领域技术人员应当将说明书作为一个整体加以理解，并将各实施例中所涉及的技术方案看作是可以相互组合成不同实施例的方式来理解本发明的保护范围。

以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的等同变化、修改与结合，均应属于本发明保护的范围。

Claims (2)

1.一种基于光谱斜率差异检测地表覆盖变化的方法，所述方法用于对同一地点的两个不同时刻的Landsat图像进行地表覆盖变化检测，从而测定所述地点发生地表覆盖变化的区域范围，所述两个不同时刻分别为T1、T2时刻，其特征在于，所述方法包括如下步骤，

A、将同一地点的Landsat图像的六个波段的光谱数据依据波谱连接成五个光谱段，分别针对所述T1、T2时刻的Landsat图像的每个点位，计算每个光谱段的斜率，光谱段(n,n+1)的斜率计算公式如下：

$k_{(n+n+1)}=\frac{\mathrm{\&Delta;ref}}{\mathrm{\&Delta;\&lambda;}}=\frac{{\mathrm{ref}}_{n+1}-{\mathrm{ref}}_n}{{\mathrm{\&lambda;}}_{n+1}^{\&prime;}-{\mathrm{\&lambda;}}_n^{\&prime;}}$

其中，Δref是波段n+1和波段n的反射率差值，Δλ是波段n+1和波段n的波段长度之差，λ'_n是归一化后无量纲的第n波段波长，λ'_n的计算公式如下：

${\mathrm{\&lambda;}}_n^{\&prime;}=\frac{{\mathrm{\&lambda;}}_n-\min \left(\mathrm{\&lambda;}\right)}{\max \left(\mathrm{\&lambda;}\right)-\min \left(\mathrm{\&lambda;}\right)}$

其中，λ_n是第n波段的波长，min(λ)是所述Landsat图像的所有波段中波长的最小值，max(λ)是所述Landsat图像的所有波段中波长的最大值，λ_n、min(λ)、max(λ)的单位均为微米；

B、根据步骤A中的斜率计算公式，分别计算所述T1、T2时刻的Landsat图像的每个点位的各波段的斜率，形成斜率向量，T1时刻的斜率向量是K_i＝(k_i(1,2),k_i(2,3),...,k_i(n,n+1))^T，T2时刻的斜率向量是K_j＝(k_j(1，2)，k_j(2，3)，…，k_j(n，n+1))^T，

通过斜率向量的绝对距离计算斜率变化强度，计算公式如下：

$\begin{array}{c}{d}_K=\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}\left|{\mathrm{\&Delta;k}}_{(n,n+1)}\right|=\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}|k_{i(n,n+1)}-k_{j(n,n+1)}|=\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}|\frac{{\mathrm{ref}}_{i,n+1}-{\mathrm{ref}}_{i,n}}{{\mathrm{\&lambda;}}_{n+1}^{\&prime;}-{\mathrm{\&lambda;}}_n^{\&prime;}}-\frac{{\mathrm{ref}}_{j,n+1}-{\mathrm{ref}}_{j,n}}{{\mathrm{\&lambda;}}_{n+1}^{\&prime;}-{\mathrm{\&lambda;}}_n^{\&prime;}}|\\ =\underset{n=1}{\overset{5}{\mathrm{\&Sigma;}}}\left|\frac{({\mathrm{ref}}_{i,n+1}-{\mathrm{ref}}_{j,n+1})-({\mathrm{ref}}_{i,n}-{\mathrm{ref}}_{j,n})}{{\mathrm{\&lambda;}}_{n+1}^{\&prime;}-{\mathrm{\&lambda;}}_n^{\&prime;}}\right|\end{array}$

其中，ref_i,n+1和ref_i,n是T1时刻第n+1波段和第n波段的反射率，ref_j,n+1和ref_j,n是T2时刻第n+1波段和第n波段的反射率，对应波段的归一化波长分别是λ'_n+1和λ'_n；d_K是两斜率向量的变化强度；

C、根据步骤B中计算出的所述T1、T2时刻的Landsat图像的每个点位的斜率变化强度计算结果d_K，通过设置阈值确定所述Landsat图像中的地表覆盖变化区域和地表覆盖不变区域。

2.根据权利要求1所述的方法，其特征在于，所述方法进一步包括如下步骤D，用于确定地表覆盖变化区域的地物变化类别，所述步骤D为，根据已知的各种地物特征点及其所对应变化类型的斜率变化特征数据，利用地表覆盖变化区域的像素和上述已知的变化类型斜率变化特征数据的比对，确定出所述像素的变化类型，从而确定出步骤C中所述Landsat图像中的地表覆盖变化区域的每个像素的地物变化类型。