大尺度时间序列遥感图像变化区域的检测方法

摘要

本发明提供了一种大尺度时间序列遥感图像变化区域的检测方法。在该检测方法中，通过将大尺度高分辨率时间序列遥感图像变化区域检测问题转化为图像区域时间序列变化信息标注问题，实现了大尺度时间序列图像中变化区域的快速定位，降低大尺度遥感图像变化检测的难度，增强大尺度时间序列图像判读解译的实用性。

说明书

技术领域

本发明涉及遥感技术领域，尤其涉及一种大尺度时间序列遥感图像变 化区域的检测方法。

背景技术

大尺度高分辨率时间序列遥感图像变化区域检测是实现土地类型长 期变化跟踪分析的关键基础，特别是随着近年来遥感数据爆炸式增加，高 分辨率遥感图像数据量极大，直接对大尺度高分辨率时间序列遥感图像进 行有监督变化检测处理困难较大：需要对每个时间序列图像均提供训练样 本，耗时费力，可操作性不强。

多示例学习方法(Multi-InstanceLearning，简称MIL)是近年来提出 的一种方便的图像分类、标注方法。与传统的图像标注方法相比，多示例 学习方法在分类器训练过程中不需要对每一个图像区域进行标记，只需要 对包的标记进行标记，指示包中是否包含示例标记，因此，在分类器训练 中比较简单方便，实际可操作性比较强。

从多示例学习思想出发，将大尺度高分辨率时间序列遥感图像变化区 域检测问题转化为图像区域时间序列变化信息标注问题，降低了有监督变 化检测对训练样本的要求：只需要标记是否发生了变化，不需要标记哪里 发生了变化，从而极大地提高了针对大尺度高分辨率时间序列遥感图像的 有监督变化检测方法的可操作性。

发明内容

(一)要解决的技术问题

鉴于上述技术问题，本发明提供了一种大尺度时间序列遥感图像变化 区域的检测方法，以提高了针对大尺度高分辨率时间序列遥感图像的有监 督变化检测方法的可操作性。

(二)技术方案

本发明大尺度时间序列遥感图像变化区域的检测方法包括：步骤A， 将大尺度时间序列遥感图像，即在同一区域不同时间获取的n幅大尺度遥 感图像，进行配准处理，得到配准后的n幅大尺度遥感图像{X₁，X₂，...，X_n}； 步骤B，将经配准处理后的在同一区域不同时间获取的n幅大尺度遥感图 像{X₁，X₂，...，X_n}按照时间序列(t₁，t₂，…，t_n)排列，按照预设大小的窗口 将每个时相上的图像均匀分割成不重叠的k个图像区域将同一位置上的图像区域按照时间序列排列其中，i＝1、 2、……，n；j＝1、2、……、k；步骤C，以每个时间序列图像区域为多示 例学习中词袋，以词袋中连续两个时间点上的图像区域为多示例学习中的 示例，通过“图像区域-示例，图像序列-词袋”的映射，在大尺度时间序 列遥感图像上完成多示例学习框架的构造；步骤D，以图像区域为单位对 每个图像区域提取颜色特征、纹理特征，SIFT特征以及熵特征，以每个图 像区域为单位构造特征向量；对词袋中的每个示例计算特征向量；以及步 骤E，利用训练样本对多示例学习框架中的分类器进行训练，得到最佳分 类器参数设置；利用已训练的多示例分类器，对测试集进行测试，给出k 个词袋的标记，从而将时间序列上发生过变化的图像区域标记出来，并以 概率形式给出每个区域发生变化的百分比。

(三)有益效果

从上述技术方案可以看出，本发明大尺度时间序列遥感图像变化区域 的检测方法具有以下有益效果：

(1)降低了有监督变化检测对训练样本的要求，极大地提高了针对 大尺度高分辨率时间序列遥感图像的有监督变化检测方法的可操作性，提 高了大尺度时间序列图像判读解译的实用性。

(2)以区域百分比形式给出每个区域发生变化的概率，方便后续进 一步的细化变化分析处理。

附图说明

图1为根据本发明实施例大尺度时间序列遥感图像变化区域的检测方 法的流程图；

图2为真实测试数据集；

图3为应用本发明实施例的检测方法对图2所示的真实测试数据集进 行检测的结果。

具体实施方式

本发明通过将大尺度高分辨率时间序列遥感图像变化区域检测问题 转化为图像区域时间序列变化信息标注问题，实现了大尺度时间序列图像 中变化区域的快速定位，降低大尺度遥感图像变化检测的难度，增强大尺 度时间序列图像判读解译的实用性。

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实 施例，并参照附图，对本发明进一步详细说明。

在本发明的一个示例性实施例中，提供了大尺度时间序列遥感图像变 化区域的检测方法。图1为根据本发明实施例大尺度时间序列遥感图像变 化区域的检测方法的流程图。如图1所示，本实施例包括以下步骤：

步骤A，图像配准预处理：将大尺度时间序列遥感图像，即在同一区 域不同时间获取的n幅大尺度遥感图像，进行配准处理，得到配准后的n 幅大尺度遥感图像{X₁，X₂，...，X_n}；

其中，大尺度时间序列遥感图像优选高分辨率光学遥感图像，其时间 序列至少包括5个以上时相，每个时相上图像的大小不小于1000*1000像 素，图像至少包含RGB三个波段。此外，此处的大尺度是指单幅图像大 小不小于1000*1000像素，高分辨率是指图像空间分辨率在2.0米～0.5米/ 像素范围内。

步骤B，时间序列图像区域划分：将经配准处理后的在同一区域不同 时间获取的n幅大尺度遥感图像{X₁，X₂，...，X_n}按照时间序列(t₁，t₂，…，t_n) 排列，按照预设大小的窗口将每个时相上的图像均匀分割成不重叠的k个 图像区域将同一位置上的图像区域按照时间序列排列其中，i＝1、2、……，n；j＝1、2、……、k；

本步骤中，窗口的大小依据感兴趣变化目标的大小确定，一般情况下， 取最小感兴趣目标面积大小的二倍以上，本发明中取50*50像素。k由图 像自身大小和窗口大小确定，当窗口大小和图像大小确定时，k为常数。

步骤C，时间序列图像的多示例框架构建：以每个时间序列图像区域 为多示例学习中词袋，以词袋中连续两个时间点上的图像区域为多示例学 习中的示例，通过“图像区域-示例，图像序列-词袋”的映射，在大尺度 时间序列遥感图像上完成多示例学习框架的构造；

其步骤为：将步骤B产生的同一位置上的按时间序列排列的图像区域 视为多示例学习框架中的一个词袋；将词袋中连续两个时 间点(t_i，t_i+1)上的图像区域视为多示例学习框架中的一个示例；如果一个 词袋中所有的示例均为不变化的区域，则该词袋标记为不变化词袋，反之， 如果词袋中有一个示例为变化的区域，则该词袋标记为变化词袋。

步骤D，以图像区域为单位对每个图像区域提取颜色特征、纹理特征， SIFT特征以及熵特征，以每个区域为单位构造特征向量；对词袋中的每个 示例计算特征向量；

该步骤又可以分为以下步骤：

分步骤D1，颜色特征主要是提取CIELab颜色特征，其步骤为：对原 始输入图像进行平滑滤波，图像中的每个像素p(m，n)，以坐标(m，n) 为中心，取3*3大小的图像区域，计算图像区域中像素的光谱均值作为当 前像素p(m，n)的滤波后光谱值，得到滤波处理后的两时相遥感图像； 将图像RGB颜色空间投影CIELab颜色空间，得到图像中每个像素在 CIELab颜色空间中的颜色特征值；以前一步获取的过分割图像为基础， 每个过分割区域中所有像素的Lab颜色特征的均值作为该区域的CIELab 特征，得到每个区域的3维CIELab颜色特征；其中，RGB颜色空间到 CIELab颜色空间的转换计算公式为：

$\begin{array}{c}X=\frac{0.412453*R+0.357580*G+0.180423*B}{0.950456}\\ Y=\frac{0.212671*R+0.715160*G+0.072169*B}{1}\\ Z=\frac{0.019334*R+0.119193*G+0.950227*B}{1.088754}\\ L=116*f\left(Y\right)-16\\ a=500*[f\left(X\right)-f\left(Y\right)]\\ b=200*[f\left(Y\right)-f\left(Z\right)]\end{array}---\left(1\right)$

其中，X、Y、Z仅为计算中间值，没有具体的物理含义。R、G和B分别 为相应像素的RGB颜色体系的R/G/B分量值。L、a和b分别为转换后的 相应像素的CIELab颜色体系的L、a和b各分量的值。

f(·)为转换函数，无具体物理含义，以t为变量的f(t)函数的定义如下：

$f\left(t\right)=\left(\begin{array}{cc}{t}^{\frac{1}{3}},& \begin{array}{cc}if& t>{\left(\frac{6}{29}\right)}^3\end{array}\\ \frac{1}{3}*{\left(\frac{29}{6}\right)}^2*t+\frac{4}{29},& other\end{array}\right)---\left(2\right)$

分步骤D2，纹理特征主要是提取Gabor纹理特征，其步骤为：将原 始输入遥感图像转化为灰度图像，转化方法为直接将图像R、G、B三个 通道上每个像素的光谱值相加求均值，作为灰度图像的光谱值；构造Gabor 滤波器如下：

其中，λ为波长，θ为Gabor滤波器的条带方向，为Gabor滤波器的相 位参数，γ为长宽比，控制Gabor滤波器的椭球率。Gabor能量函数定义 如下形式：

$e_{\lambda ,\theta }(x,y)=\sqrt{{\gamma }_{\lambda ,\theta ,0}^2(x,y)+{\gamma }_{\lambda ,\theta ,-\frac{\pi }{2}}^2(x,y)}---\left(4\right)$

其中，和分别为g_λ，θ，0(x，y)和与灰度图像卷积的 结果；本实施例中设定Gabor滤波器的尺度参数λ分别取[12345]，方向 参数θ分别取[0π/4(2*π)/4(3*π)/4(4*π)/4(5*π)/4(6* π)/4(7*π)/4](π≈3.1415926)，可以构造40中Gabor滤波器，使 用每个滤波器对灰度图像滤波，即将每个Gabor滤波器与灰度图像进行卷 积；以图像过分割区域为单位，区域中所有像素的Gabor纹理特征的均值 作为该区域的Gabor特征，得到每个区域的40维Gabor纹理特征；

分步骤D3，SIFT特征(D.Lowe，“Distinctiveimagefeaturesfromscale invariantkeypoints，”Int.J.ComputerVision，vol.60，no.2，pp.91-110，Dec. 2004.)提取步骤为：

首先，对图像进行归一化预处理，再将图像放大两倍并预滤波剔除噪 声。定义图像的尺度空间为L(m，n，σ)，它是由可变尺度高斯G(m，n，σ)和输 入图像I(m，n)卷积得到：

$L(m,n,\sigma )=G(m,n,\sigma )\otimes I(m,n)---\left(5\right)$

式中，为m和n的卷积算子，σ为尺度因子，(m，n)为图像坐标；

其次，在图像二维平面空间和高斯差分(DOG，DifferenceofGaussian) 尺度空间中同时检测局部极值以作为特征点，DOG算子如下：

D(m，n，σ)＝L(m，n，kσ)-L(m，n，σ)(6)

式中，k为常数因子；

再次，得到精确特征点后，去除其中对比度低和稳定性差的边缘特征 点；

最后，当检测完图像的特征点后，这时每个特征点包含了三个信息： 位置、所在尺度、方向。SIFT特征是在此关键点及其邻近的4×4的小块 中，计算8个方向的梯度直方图，得到一个128维的向量。此时SIFT特 征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征 向量的长度归一化，则可以进一步去除光照变化的影响。

对全部用来测试的图像提取了SIFT特征。

分步骤D4，熵特征提取步骤为：将原始输入遥感图像转化为灰度图 像，转化方法为直接将图像R、G、B三个通道上每个像素的光谱值相加 求均值，作为灰度图像的光谱值；以过分割图像中每个区域为单位，计算 图像区域的熵，得到每个区域的1维熵特征。熵特征的计算公式如下：

$entry=\underset{i=1}{\overset{255}{\Sigma }}{p}_i{\mathrm{logp}}_i.---\left(7\right)$

其中，p_i表示图像区域中灰度值为i的像素个数占整个图像区域中像素个 数的比例。

分步骤D5，计算每个图像区域内像素各特征均值作为该待处理图像 区域的特征；

分步骤D6，分别对分步骤D1～分步骤D5中获取的每个图像区域的颜 色特征、纹理特征、SIFT特征和熵特征进行归一化处理，将归一化之后的 颜色特征、纹理特征、SIFT特征和熵特征合并为一个特征向量，作为该图 像区域的特征向量。其中，归一化处理公式表示如下：

f′＝(f-min(min(f)))/(max(max(f))-min(min(f)))(8)

其中，f分别为上述的CIELab颜色特征、Gabor纹理特征、SIFT特征、熵 特征，f′为归一化处理后的特征向量。

分步骤D7，按照步骤C中示例构建方法，对每个示例计算特征向量。

对于每个示例，利用分步骤D6中获取的归一化特征向量，计算每个 示例的特征向量如下：

$f_{ins\tan ce}=abs(f_{t_i}^{\prime }-f_{t_{i+1}}^{\prime }).---\left(9\right)$

其中，和分别为时间点t_i和t_i+1上图像区域的经公式(8)归一化处理 后的特征向量。

步骤E，利用训练样本对多示例学习框架中的分类器进行训练，得到 最佳分类器参数设置；利用已训练的多示例分类器，对测试集进行测试， 给出k个词袋的标记，从而将时间序列上发生过变化的图像区域标记出来， 并以概率形式给出每个区域发生变化的百分比，以供后续变化分析使用

该步骤又可以分为以下步骤：

分步骤E1，训练样本数据集准备：其中，X_i^*为某一大 尺度高分辨率遥感图像序列{X₁^*，X₂^*，...，X_l^*}按照时间序列(t₁，t₂，…，t_l)排 列，按照步骤B和步骤C将X_i^*构造为词袋和示例；按照步骤D计算每个 示例的特征向量Y_i^*为训练样本X_i^*中对应每个词袋的类别概率标记， 指示训练样本X_i^*中每个区域属于变化类或者属于不变化类的概率；如果 一个词袋中所有的示例均为不变化的区域，则该词袋标记为不变化词袋 Y_i^*＝-1，反之，如果词袋中有一个示例为变化的区域，则该词袋标记为变 化词袋Y_i^*＝+1；

分步骤E2，分类器训练：本发明中采用支持向量机(SupportVector Machine，SVM)分类器，分类器使用高斯径向基核函数(GaussianRadial BasisFunction，RBF)，表达式如下：

$k(f_{ins\tan ce\left(i\right)}^{*},f_{ins\tan ce\left(j\right)}^{})={\theta }_0\exp (-\frac{|f_{ins\tan ce\left(i\right)}^{*}-f_{ins\tan ce\left(j\right)}^{*}{|}^2}{2{\sigma }^2})---\left(10\right)$

其中，θ₀为高斯径向基核函数的幅值控制参数，σ为高斯径向基核函数的 尺度控制参数；分别表示步骤D中获取的第i个示例和第j 个示例的特征向量集合，θ₀和σ即为需要训练设置的两个高斯核参数；利 用给定的训练样本集对应的特征向量集合和标记集合 Y_i^*，对支持向量机分类器进行训练，得到最佳的分类器参数设置；

分步骤E3，利用步骤E2训练所得的支持向量机分类器，对待测试数 据X_i＝{X₁，X₂，...，X_n}由步骤D所获取的特征向量组成的测试特征集f_instance进行分类，得到每个词袋的标注类别概率Y_i。

以下进一步说明本发明的应用实例：实例给出了一组QuickBird遥感 图像序列。图2所示，给出了时间序列数据样本集合。该图像序列包含8 幅图像，每幅图像的大小为16476*16386像素，图像序列的时间跨度为 2002年至2013。将本实施例的检测方法应用于真实数据集上的部分测试 应用实例如图3所示。

需要说明的是，上述对各元件的定义并不仅限于实施方式中提到的各 种具体结构或形状，本领域的普通技术人员可对其进行简单地熟知地替换， 例如：

(1)步骤C中图像区域均匀分割的窗口大小的选择可以根据应用目 的进行经验选取；

(2)步骤D中特征提取方法可以使用其他类型图像特征描述方法， 对图像进行特征提取；

(3)步骤E中分类器的选取可以使用其他分类器，例如Boosting分 类器、高斯过程分类器、KNN分类器等，这些分类器均为本领域内公知 的方法，此处不再详细描述。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行 了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而 已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修 改、等同替换、改进等，均应包含在本发明的保护范围之内。