CN102542295B - 一种采用图像分类技术从遥感图像中进行滑坡检测的方法 - Google Patents

Abstract

本发明涉及一种采用图像分类技术从遥感图像中进行滑坡检测的方法，首先采用基于重叠的面积均分法对预处理后的遥感图像进行分块，得到面积相等的正方形图像块，并将得到的图像块分成两个集合——训练集和测试集；其次，提取训练集和测试集中所有图像块的SIFT特征，对训练集中的SIFT特征采用k-means聚类方法得到单词和词典；然后，用BoVW模型表示训练集和测试集中的每个图像块；最后，利用pLSA模型提取每个图像块的主题，并用KNN分类器将测试集中的图像块分成滑坡和非滑坡两类，从而实现对遥感图像的滑坡检测。该方法运算量小、检测效率高，同时又具有很高的检测正确率。

Description

一种采用图像分类技术从遥感图像中进行滑坡检测的方法

技术领域

本发明涉及一种采用图像分类技术从遥感图像中进行滑坡检测的方法，可以应用于各种遥感图像的滑坡检测。 

背景技术

滑坡是一类破坏力极强的地质灾害，它常常给工农业生产以及人民生命财产造成巨大损失。因此，滑坡研究也越来越受到人们的重视。早期的滑坡检测主要采用传统的地面调查方法，该方法野外工作强度很大，特别在危险及气候恶劣地区，工作效率十分低下。到20世纪90年代末，立体镜航空照片解译配合一定的地面验证成为滑坡检测和制图最常用的方法，然后将解译结果通过手工转绘到相应比例尺的地形图上，制作滑坡分布图，该方法使部分野外工作转移到室内，在一定程度上提高了工作效率，减轻了野外工作强度。近年来，开始利用ArcGIS、CoreIDRAW等软件平台结合数字高程模型进行滑坡检测和分析，但该方法需要获取三维立体地形学信息，计算量很大。 

遥感作为一种滑坡调查和监测手段，随着传感器的不断发展和遥感图像分辨率的不断提高，越来越受到地质灾害研究人员的关注。从遥感图像中进行滑坡检测是滑坡分析、评价、预测和监测的基础，它利用滑坡发生处的影像光谱信息、地形地貌和形态特征等，对遥感图像进行分析和判别，获取滑坡灾害发生范围，从而达到灾害调查和制图的目的。近年来，随着高分辨率遥感图像在资源环境监测、灾害管理等领域的广泛应用，高分辨率遥感图像信息提取方法以及遥感图像分类方法已成为近年来研究的热点。现有的采用遥感图像进行滑坡检测的方法主要有基于对象的高分辨率遥感图像滑坡检测方法和基于监督分类的遥感图像滑坡检测方法。但是，基于对象的高分辨率遥感图像滑坡检测方法需要对多期遥感图像进行多尺度、多层次分割，检测效果直接依赖于图像分割的结果，总体检测率不高；基于监督分类的遥感图像滑坡检测方法 只能实现将图像划分为滑坡解译所需要的简单几种类型，且检测过程中仍需借助于遥感专业处理软件ENVI，检测效率和检测精度都有待提高。 

发明内容

要解决的技术问题 

为了避免现有技术的不足之处，本发明提出了一种采用图像分类技术从遥感图像中进行滑坡检测的方法，能够从遥感图像中进行滑坡检测，且具有很高的检测正确率。 

本发明的思想在于：(1)采用基于重叠的面积均分法对预处理后的遥感图像进行分块，得到面积相等的正方形图像块，然后将得到的图像块分成两个集合——训练集和测试集；(2)提取训练集和测试集中所有图像块的SIFT特征，然后对训练集中的SIFT特征采用k-means聚类方法得到单词和词典；(3)用BoVW模型表示训练集和测试集中的每个图像块；(4)利用pLSA模型提取每个图像块的主题，最后用KNN分类器将测试集中所有的图像块分成滑坡和非滑坡两类，从而实现对遥感图像的滑坡检测。 

技术方案 

一种采用图像分类技术从遥感图像中进行滑坡检测的方法，其特征在于：步骤1：预处理：采用加权平均法对遥感图像的RGB三个分量进行加权平均得到灰度图像，然后利用灰度线性变换函数将灰度图像的灰度范围映射到 

的灰度区间，得到预处理后的图像； 

其中：加权平均法计算公式为f(i，j)＝0.3R(i，j)+0.59G(i，j)+0.11B(i，j)，f(i，j)为加权平均法得到的(i，j)像素点的灰度值，R(i，j)、G(i，j)和B(i，j)分别为(i，j)像素点的RGB三个分量值；灰度线性变换函数为 

为(i，j)像素点经过灰度映射后的灰 度值，f_min、f_max分别为灰度映射前的灰度图像的最小灰度值和最大灰度值；步骤2：图像分块：采用2m×2m的滑动窗口对预处理后的灰度图像进行重叠分块，相 

邻图像块的重叠区域为m个像素，得到图像块的集合D＝{I₁，…，I_N}；从集合D中 选取N_landslide个具有滑坡特征的滑坡类图像块和N_{non-landslide}个不具有滑坡特征的非滑坡类图像块组成训练集 

将集合D中剩余的图像块组成测试集  $B=\{b_1,L,b_{N_{\mathrm{test}}}\};$

其中：I表示集合D中的图像块；下标N为集合D中图像块的个数；m的取值范围为 

γ为遥感图像的空间分辨率，单位为米；a表示训练集中的图像块；下标N_training为训练集中图像块的个数，且N_training＝N_landslide+N_{non-landslide}，N_landslide/N_{non-landslide}≈0.5；b表示测试集中的图像块；下标N_test为测试集中图像块的个数； 

步骤3：单词和词典的形成：共分两个步骤，具体过程如下： 

步骤a1：提取集合D中每个图像块的SIFT特征，具体过程如下： 

首先计算每个图像块的图像卷积 

其中：σ_n＝0.5，*表示卷积运算； 

然后计算每个图像块的图像金字塔 

图像金字塔共O组，每组S层，下一组的第一幅图像由上一组的最后一幅图像降采样得到，采样因子为2；其中：σ＝σ₀·2^o+s/S，σ₀＝1.6·2^1/S，o＝0，…, O-1，s＝0，…,  S-1，S＝3，O＝log₂ 2m； 

再对每个图像块中同一组相邻的GSS_σ求差分得到DOG_σ，将DOG_σ每个像素点的值分别和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，如果该像素点的值为极小值或者极大值，则该像素点为图像显著点，显著点周围以σ为半径的区域为显著区域，由此可以得到一系列的图像显著点； 

最后计算每个图像块的梯度图像卷积 

其中： 

为 

的梯度图像；然后在 

上计算以图像显著点为中心、σ为半径的圆形区域的梯度方向直方图，梯度方向直方图共分36个区间，每个区间10度，梯度方向直方图峰值的方向区域为该显著点的主方向；最后在 上将以每个图像显著点为中心、σ为半径的圆形区域，按主方向及其垂直方向等分成16个区域，在每个小区域中分别统计梯度方向直方图，梯度方向直方图共分8个区间，每个区间45度，并将每个方向直方图的幅值量化到[0，255]区间，得到16×8＝128维的SIFT特征；由此，每一个图像块I可以得到一系列的SIFT特 征； 

步骤b1：设定聚类中心数目M，利用K-means聚类算法对训练集中所有图像块的SIFT特征进行聚类，得到含有M个单词的词典W＝{w₁，…，w_M}； 

其中：M表示词典中单词的个数，取值范围为500～2500；w表示词典中的单词； 

步骤4：用BoVW模型表示图像块：将训练集和测试集中所有图像块的每个SIFT特征映射为词典中的一个单词；然后统计每一个图像块中所有单词的频率，得到一个N_training×M的矩阵 和一个N_test×M的矩阵  $N_{\mathrm{test}}={\left(n{(b_k,w_j)}_{\mathrm{kj}}\right)}_{N_{\mathrm{test}}\×M};$

其中：n(a_i，w_j)为图像块a_i的BoVW模型表示，代表单词w_j在图像块a_i中出现的频率；下标ij表示元素n(a_i，w_j)在N_training中处于第i行、第j列，i＝1，2，…，N_training，j＝1，2，…，M；n(b_k，w_j)为图像块b_k的BoVW模型表示，代表单词w_j在图像块b_k中出现的频率；下标kj表示元素n(b_k，w_j)在N_test中处于第k行、第j列，k＝1，2，…，N_test； 

步骤5：用pLSA模型提取图像块的主题：共分两个步骤，具体过程如下： 

步骤a2：计算训练集图像块的P(w_j|z_t)和P(z_t|a_i)，具体过程如下： 

初始化 $P\left(z_t\right)=\frac{1}{T},$ $P\left(w_j\right|z_t)=\frac{1}{M}$ 和 $P\left(a_i\right|z_t)=\frac{1}{N_{\mathrm{training}}},$

计算 $P\left(z_t\right|a_i,w_j)=\frac{P\left(z_t\right)P\left(w_j\right|z_t)P\left(a_i\right|z_t)}{{\mathrm{\Σ}}_{m=1}^{T}P\left(z_m\right)P\left(w_j\right|z_m)P\left(a_i\right|z_m)},$

将P(z_t|a_i，w_j)代入公式 $\left\{\begin{array}{c}P\left(w_j\right|z_t)=\frac{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{r=1}^M{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}n(a_i,w_r)P\left(z_t\right|a_i,w_r)}\\ P\left(a_i\right|z_t)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{s=1}^{N_{\mathrm{training}}}n(a_s,w_j)P\left(z_t\right|a_s,w_j)}\\ P\left(z_t\right)=\frac{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}P\left(a_i\right|z_t)}{{\mathrm{\Σ}}_{m=1}^T{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}P\left(a_i\right|z_m)}\end{array}\right.$

计算P(w_j|z_t)、P(a_i|z_t)和P(z_t)，然后再次计算P(z_t|a_i，w_j)并重复上述过程，直到 $L={\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)\log [{\mathrm{\Σ}}_{t=1}^{T}P\left(z_t\right)P\left(w_j\right|z_t)P\left(a_i\right|z_t)]$ 的对数似然函数期望 值的增加量小于指定的阈值δ时停止迭代，得到P(w_j|z_t)；然后，得到 

$P\left(z_t\right|a_i)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)};$

其中：z表示主题；P(w_j|z_t)表示单词w_j在主题z_t中出现的概率；P(z_t|a_i)表示主题z_t在图像块a_i中出现的概率；P(z_t)表示主题z_t在所有主题中出现的概率；P(a_i|z_t)表示图像块a_i在主题z_t中出现的概率；P(z_t|a_i，w_j)表示主题z_t在图像块a_i和单词w_j条件下的后验概率；下标t和m表示主题的序号，t＝1，2，…，T，m＝1，2，…，T，T表示主题个数，T的取值范围为10～25；下标r表示单词的序号，r＝1，2，…，M；下标s表示训练集图像块的序号，s＝1，2，…，N_training；log表示对数运算符；δ为设定的阈值，取值范围为0.1～0.5； 

步骤b2：计算测试集图像块的P(z_t|b_k)，具体过程如下： 

初始化 $P\left(z_t\right)=\frac{1}{T},$ $P\left(b_k\right|z_t)=\frac{1}{N_{\mathrm{test}}},$

计算 $P\left(z_t\right|b_k,w_j)=\frac{P\left(z_t\right)P\left(w_j\right|z_t)P\left(b_k\right|z_t)}{{\mathrm{\Σ}}_{m=1}^{T}P\left(z_m\right)P\left(w_j\right|z_m)P\left(b_k\right|z_m)},$

将P(z_t|b_k，w_j)代入公式 $\left\{\begin{array}{c}P\left(b_k\right|z_t)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)P\left(z_t\right|b_k,w_j)}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{d=1}^{N_{\mathrm{test}}}n(b_d,w_j)P\left(z_t\right|b_d,w_j)}\\ P\left(z_t\right)=\frac{{\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}P\left(b_k\right|z_t)}{{\mathrm{\Σ}}_{m=1}^T{\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}P\left(b_k\right|z_m)}\end{array}\right.$

计算P(b_k|z_t)和P(z_t)，然后再次计算P(z_t|b_k，w_j)并重复上述过程，直到  $L={\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)\log \left[{\mathrm{\Σ}}_{t=1}^{T}P\left(z_t\right)P\left(w_j\right|z_t)P\left(b_k\right|z_t)\right]$ 的对数似然函数期望值的增加量小于指定的阈值δ时停止迭代，计算 

其中：P(z_t|b_k)表示主题z_t在图像块b_k中出现的概率；P(b_k|z_t)表示图像块b_k在主题z_t中出现的概率；P(z_t|b_k，w_j)表示主题z_t在图像块b_k和单词w_j条件下的后验概率；下标d表示测试集图像块的序号，d＝1，2，…，N_test； 

步骤6：KNN分类：计算P(z_t|b_k)与P(z_t|a_i)的欧式距离，得到一个N_test×N_training的欧式距离矩阵 

然后根据设定的近邻数目K，利用KNN分类器将 测试集中所有的图像块分成滑坡和非滑坡两类； 

其中： 

表示P(z_t|b_k)与P(z_t|a_i)的欧式距离，计算公式为 

下标ki表示元素 

在Θ中处于第k行、第i列；K为KNN分类器的近邻数目，取值范围为7～15。 

有益效果 

本发明提出的采用图像分类技术从遥感图像中进行滑坡检测的方法，无需获取待检测地点的三维立体地形学信息，直接采用基于BoVW和pLSA的图像分类方法将分块后的待检测遥感图像分为两类——滑坡类和非滑坡类，从而实现从遥感图像中进行滑坡检测。该方法运算量小、检测效率高，同时又具有很高的检测正确率。采用此检测方法得到的滑坡对于后续的滑坡分析、评价、预测和监测都具有重要意义。 

附图说明

图1：本发明方法的基本流程图 

图2：本发明方法的示意图 

具体实施方式

现结合实施例、附图对本发明作进一步描述： 

用于实施的硬件环境是：Intel(R)Core(TM)i5CPU 3.2G计算机、2GB内存、1G显卡，运行的软件环境是：Matlab R2008a和Windows XP。我们用Matlab软件实现了本发明提出的方法。该实验所用的遥感图像取自于新疆伊犁自治州的巩留县和新源县的三个滑坡重灾区，具体的地理坐标为：东经82°32′41″-82°39′33″，北纬43°06′13″-43°12′02″；东经82°44′00″-83°10′50″E，北纬43°08′57″-43°14′21″；东经83°17′41″-83°35′15″，北纬43°12′02″-43°31′32″。 

本发明具体实施如下： 

1、预处理：采用加权平均法对遥感图像的RGB三个分量进行加权平均得到灰度图像，然后利用灰度线性变换函数将灰度图像的灰度范围映射到[0，255]的灰度区间，得到预处理后的图像。 

其中，加权平均法计算公式为f(i，j)＝0.3R(i，j)+0.59G(i，j)+0.11B(i，j)；灰度线性变 换函数为 

f(i，j)为加权平均法得到的(i，j)像素点的灰度值；R(i，j)、G(i，j)和B(i，j)分别为(i，j)像素点的RGB三个分量值； 

为(i，j)像素点经过灰度映射后的灰度值，f_min、f_max分别为灰度映射前的灰度图像的最小灰度值和最大灰度值。 

2、图像分块：实验所用的遥感图像的空间分辨率为1米，选取m＝256，采用512×512的滑动窗口对预处理后的灰度图像进行重叠分块，相邻图像块的重叠区域为256个像素，得到图像块的集合D＝{I₁，…，I_N}；从集合D中选取66个具有滑坡特征的滑坡类图像块和134个不具有滑坡特征的非滑坡类图像块组成训练集 

将集合D中剩余的图像块组成测试集 

其中：N＝21000，N_training＝200，N_test＝20800。 

3、单词和词典的形成：共分两个步骤，具体过程如下： 

(1)提取集合D中21000个图像块的SIFT特征，具体过程如下： 

1)计算每个图像块的图像卷积 

其中：σ_n＝0.5，*表示卷积运算； 

2)计算每个图像块的图像金字塔 

图像金字塔共O组，每组S层，下一组的第一幅图像由上一组的最后一幅图像降采样得到，采样因子为2；其中：σ＝σ₀·2^o+s/S，σ₀＝1.6·2^1/S，o＝0，…, O-1，s＝0，…,  S-1，S＝3，O＝9； 

3)对每个图像块中同一组相邻的GSS_σ求差分得到DOG_σ，将DOG_σ每个像素点的值分别和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，如果该像素点的值为极小值或者极大值，则该像素点为图像显著点，显著点周围以σ为半径的区域为显著区域，由此可以得到一系列的图像显著点； 

4)计算每个图像块的梯度图像卷积 其中： 为 

的梯度图像；然后在 

上计算以图像显著点为中心、σ为半径的圆形区域的梯度方向直方图，梯度方向直方图共分36个区间，每个区间 10度，梯度方向直方图峰值的方向区域为该显著点的主方向；最后在 

上将以每个图像显著点为中心、σ为半径的圆形区域，按主方向及其垂直方向等分成16个区域，在每个小区域中分别统计梯度方向直方图，梯度方向直方图共分8个区间，每个区间45度，并将每个方向直方图的幅值量化到[0，255]区间，得到16×8＝128维的SIFT特征；由此，每一个图像块I可以得到一系列的SIFT特征； 

(2)设定聚类中心数目M＝800，利用K-means聚类算法对训练集中200个图像块的SIFT特征进行聚类，得到含有800个单词的词典W＝{w₁，…，w₈₀₀}。 

4、用BoVW模型表示图像块：将训练集和测试集中所有图像块的每个SIFT特征映射为词典中的一个单词；然后统计每一个图像块中所有单词的频率，得到一个200×800的矩阵N_training＝(n(a_i，w_j)_ij)_200×800和一个20800×800的矩阵N_test＝(n(b_k，w_j)_kj)_20800×800；其中：i＝1，2，…，200，j＝1，2，…，800，k＝1，2，…，20800。 

5、用pLSA模型提取每个图像块的主题：共分两个步骤，具体过程如下： 

(1)计算训练集图像块的P(w_j|z_t)和P(z_t|a_i)，具体过程如下：选取δ＝0.2，T＝20， 

初始化P(z_t)＝0.05、P(w_j|z_t)＝0.00125和P(a_i|z_t)＝0.005，计算  $P\left(z_t\right|a_i,w_j)=\frac{P\left(z_t\right)P\left(w_j\right|z_t)P\left(a_i\right|z_t)}{{\mathrm{\Σ}}_{m=1}^{20}P\left(z_m\right)P\left(w_j\right|z_m)P\left(a_i\right|z_m)},$ 将P(z_t|a_i，w_j)代入公式  $\left\{\begin{array}{c}P\left(w_j\right|z_t)=\frac{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{r=1}^M{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}n(a_i,w_r)P\left(z_t\right|a_i,w_r)}\\ P\left(a_i\right|z_t)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{s=1}^{n_{\mathrm{training}}}n(a_s,w_j)P\left(z_t\right|a_s,w_j)}\\ P\left(z_t\right)=\frac{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}P\left(a_i\right|z_t)}{{\mathrm{\Σ}}_{m=1}^T{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}P\left(a_i\right|z_m)}\end{array}\right.$ 计算P(w_j|z_t)、P(a_i|z_t)和P(z_t)，然后再次计算P(z_t|a_i，w_j)并重复上述过程，直到  $L={\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)\log \left[{\mathrm{\Σ}}_{t=1}^{T}P\left(z_t\right)P\left(w_j\right|z_t)P\left(a_i\right|z_t)\right]$ 的对数似然函数期望值的增加量小于指定的阈值δ时停止迭代，得到P(w_j|z_t)；然后，得到 

$P\left(z_t\right|a_i)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)};$

(2)计算测试集图像块的P(z_t|b_k)，具体过程如下：选取δ＝0.2，T＝20，初始化P(z_t)＝0.05、 $P\left(b_k\right|z_t)=\frac{1}{20800},$ 计算 $P\left(z_t\right|b_k,w_j)=\frac{P\left(z_t\right)P\left(w_j\right|z_t)P\left(b_k\right|z_t)}{{\mathrm{\Σ}}_{m=1}^{T}P\left(z_m\right)P\left(w_j\right|z_m)P\left(b_k\right|z_m)},$

将P(z_t|b_k，w_j)代入公式 $\left\{\begin{array}{c}P\left(b_k\right|z_t)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)P\left(z_t\right|b_k,w_j)}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{d=1}^{N_{\mathrm{test}}}n(b_d,w_j)P\left(z_t\right|b_d,w_j)}\\ P\left(z_t\right)=\frac{{\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}P\left(b_k\right|z_t)}{{\mathrm{\Σ}}_{m=1}^T{\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}P\left(b_k\right|z_m)}\end{array}\right.$ 计算P(b_k|z_t)和P(z_t)，然后再次计算P(z_t|b_k，w_j)并重复上述过程，直到  $L={\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)\log \left[{\mathrm{\Σ}}_{t=1}^{T}P\left(z_t\right)P\left(w_j\right|z_t)P\left(b_k\right|z_t)\right]$ 的对数似然函数期望值的增加量小于指定的阈值δ时停止迭代，计算 

6、KNN分类：计算P(z_t|b_k)与P(z_t|a_i)的欧式距离，得到一个20800×200的欧式距离矩阵 选取近邻数目K＝11，利用KNN分类器将测试集中所有的图像块分成滑坡和非滑坡两类；其中： 

的计算公式为  ${\mathrm{\ρ}}_{b_k,a_i}=\sqrt{{\mathrm{\Σ}}_{t=1}^T{[P\left(z_t\right|a_i)-P\left(z_t\right|b_k)]}^2}.$

采用本发明方法可获得91.23％的检测正确率。其中，检测正确率定义为正确检测的滑坡数量与总的滑坡数量之比。 

Claims (1)

1.一种采用图像分类技术从遥感图像中进行滑坡检测的方法，其特征在于：

步骤1：预处理：采用加权平均法对遥感图像的RGB三个分量进行加权平均得到灰度图像，然后利用灰度线性变换函数将灰度图像的灰度范围映射到

的灰度区间，得到预处理后的图像；

其中：加权平均法计算公式为f(i,j)＝0.3R(i,j)+0.59G(i,j)+0.11B(i,j)，f(i,j)为加权平均法得到的(i,j)像素点的灰度值，R(i,j)、G(i,j)和B(i,j)分别为(i,j)像素点的RGB三个分量值；灰度线性变换函数为

为(i,j)像素点经过灰度映射后的灰度值，f_min、f_max分别为灰度映射前的灰度图像的最小灰度值和最大灰度值；

步骤2：图像分块：采用2m×2m的滑动窗口对预处理后的灰度图像进行重叠分块，相邻图像块的重叠区域为m个像素，得到图像块的集合D＝{I₁,…，I_N}；从集合D中选取N_landslide个具有滑坡特征的滑坡类图像块和N_{ono-landslide}个不具有滑坡特征的非滑坡类图像块组成训练集

将集合D中剩余的图像块组成测试集 $B=\{b_1,\·\·\·,{b_N}_{\mathrm{test}}\};$

其中：I表示集合D中的图像块；下标N为集合D中图像块的个数；m的取值范围为

γ为遥感图像的空间分辨率，单位为米；a表示训练集中的图像块；下标N_train为训练集中图像块的个数，且N_training＝N_landslide+N_{ono-landslide}，N_landslide/N_{ono-landslide}≈0.5；b表示测试集中的图像块；下标N_test为测试集中图像块的个数；

步骤3：单词和词典的形成：共分两个步骤，具体过程如下：

步骤a1：提取集合D中每个图像块的SIFT特征，具体过程如下：

首先计算每个图像块的图像卷积

其中：σ_n＝0.5，*表示卷积运算；

然后计算每个图像块的图像金字塔

图像金字塔共O组，每组S层，下一组的第一幅图像由上一组的最后一幅图像降采样得到，采样因子为2；其中：σ＝σ₀·2^o+s/S，σ₀＝1.6·2^1/S，o＝0,...0-1，s＝0,.00S-1，S＝3，O＝log₂2m；

再对每个图像块中同一组相邻的GSS_σ求差分得到DOG_σ，将DOG_σ每个像素点的值分别和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，如果该像素点的值为极小值或者极大值，则该像素点为图像显著点，显著点周围以σ为半径的区域为显著区域，由此可以得到一系列的图像显著点；

最后计算每个图像块的梯度图像卷积 ${\hat{I}}_{{\mathrm{\σ}}_n}=I_{G_{{\mathrm{\σ}}_n}}*\frac{1}{2\mathrm{\π}\·{\left(1.5\mathrm{\σ}\right)}^2}{e}^{-(x^2+y^2)/2\·{\left(1.5\mathrm{\σ}\right)}^2},$ 其中：为

的梯度图像；然后在

上计算以图像显著点为中心、σ为半径的圆形区域的梯度方向直方图，梯度方向直方图共分36个区间，每个区间10度，梯度方向直方图峰值的方向区域为该显著点的主方向；最后在

上将以每个图像显著点为中心、σ为半径的圆形区域，按主方向及其垂直方向等分成16个区域，在每个小区域中分别统计梯度方向直方图，梯度方向直方图共分8个区间，每个区间45度，并将每个方向直方图的幅值量化到[0,255]区间，得到16×8＝128维的SIFT特征；由此，每一个图像块I可以得到一系列的SIFT特征；

步骤b1：设定聚类中心数目M，利用K-means聚类算法对训练集中所有图像块的SIFT特征进行聚类，得到含有M个单词的词典W＝{w₁,…，w_M}；

其中：M表示词典中单词的个数，取值范围为500～2500；w表示词典中的单词；

步骤4：用BoVW模型表示图像块：将训练集和测试集中所有图像块的每个SIFT特征映射为词典中的一个单词；然后统计每一个图像块中所有单词的频率，得到一个N_training×M的矩阵

和M一个

的矩阵

其中：n(a_i,w_j)为图像块a_i的BoVW模型表示，代表单词w_j在图像块a_i中出现的频率；下标ij表示元素n(a_i,w_j)在N_training中处于第i行、第j列，i＝1,2,…，N_training，j＝1,2,…，M；n(b_k,w_j)为图像块b_k的BoVW模型表示，代表单词w_j在图像块b_k中出现的频率；下标kj表示元素n(b_k,w_j)在N_test中处于第k行、第j列，k＝1,2,…，N_test;

步骤5：用pLSA模型提取图像块的主题：共分两个步骤，具体过程如下：

步骤a2：计算训练集图像块的P(w_j|z_t)和P(z_t|a_i)，具体过程如下：

初始化 $P\left(z_t\right)=\frac{1}{T},$ $P\left(w_j\right|z_t)=\frac{1}{M}$ 和 $P\left(a_i\right|z_t)=\frac{1}{N_{\mathrm{training}}},$

计算 $P\left(z_t\right|a_i,w_j)=\frac{P\left(z_t\right)P\left(w_j\right|z_t)P\left(a_i\right|z_t)}{{\mathrm{\Σ}}_{m=1}^{T}P\left(z_m\right)P\left(w_j\right|z_m)P\left(a_i\right|z_m)},$

将P(z_t|a_i,w_j)代入公式 $\left\{\begin{array}{c}P\left(w_j\right|z_t)=\frac{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{r=1}^M{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}n(a_i,w_r)P\left(z_t\right|a_i,w_r)}\\ P\left(a_i\right|z_t)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{s=1}^{N_{\mathrm{training}}}n(a_s,w_j)P\left(z_t\right|a_s,w_j)}\\ P\left(z_t\right)=\frac{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}P\left(a_i\right|z_t)}{{\mathrm{\Σ}}_{m=1}^T{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}P\left(a_i\right|z_m)}\end{array}\right.$

计算P(w_j|z_t)、P(a_i|z_t)和P(z_t)，然后再次计算P(z_t|a_i,w_j)并重复上述过程，直到 $L={\mathrm{\Σ}}_{i=1}^{N_{\mathrm{training}}}{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)\log \left[{\mathrm{\Σ}}_{t=1}^{T}P\left(z_t\right)P\left(w_j\right|z_t)P\left(a_i\right|z_t)\right]$ 的对数似然函数期望值的增加量小于指定的阈值δ时停止迭代，得到P(w_j|z_t)；然后，得到 $P\left(z_t\right|a_i)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)P\left(z_t\right|a_i,w_j)}{{\mathrm{\Σ}}_{j=1}^{M}n(a_i,w_j)};$

其中：z表示主题；P(w_j|z_t)表示单词w_j在主题z_t中出现的概率；P(z_t|a_i)表示主题z_t在图像块a_i中出现的概率；P(z_t)表示主题z_t在所有主题中出现的概率；P(a_i|z_t)表示图像块a_i在主题z_t中出现的概率；P(z_t|a_i,w_j)表示主题z_t在图像块a_i和单词w_j条件下的后验概率；下标t和m表示主题的序号，t＝1,2,…，T，m＝1,2,…，T，T表示主题个数，T的取值范围为10～25；下标r表示单词的序号，r＝1,2,…,M；下标s表示训练集图像块的序号，s＝1,2,…,N_training；log表示对数运算符；δ为设定的阈值，取值范围为0.1～0.5；

步骤b2：计算测试集图像块的P(z_t|b_k)，具体过程如下：

初始化 $P\left(z_t\right)=\frac{1}{T},$ $P\left(b_k\right|z_t)=\frac{1}{N_{\mathrm{test}}},$

计算 $P\left(z_t\right|b_k,w_j)=\frac{P\left(z_t\right)P\left(w_j\right|z_t)P\left(b_k\right|z_t)}{{\mathrm{\Σ}}_{m=1}^{T}P\left(z_m\right)P\left(w_j\right|z_m)P\left(b_k\right|z_m)},$

将P(z_t|b_k,w_j)代入公式 $\left\{\begin{array}{c}P\left(b_k\right|z_t)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)P\left(z_t\right|b_k,w_j)}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{d=1}^{N_{\mathrm{test}}}n(b_d,w_j)P\left(z_t\right|b_d,w_j)}\\ P\left(z_t\right)=\frac{{\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}P\left(b_k\right|z_t)}{{\mathrm{\Σ}}_{m=1}^T{\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}P\left(b_k\right|z_m)}\end{array}\right.$

计算P(b_k|z_t)和P(z_t)，然后再次计算P(z_t|b_k,w_j)并重复上述过程，直到 $L={\mathrm{\Σ}}_{k=1}^{N_{\mathrm{test}}}{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)\log \left[{\mathrm{\Σ}}_{t=1}^{T}P\left(z_t\right)P\left(w_j\right|z_t)P\left(b_i\right|z_t)\right]$ 的对数似然函数期望值的增加量小于指定的阈值δ时停止迭代，计算 $P\left(z_t\right|b_k)=\frac{{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)P\left(z_t\right|b_k,w_j)}{{\mathrm{\Σ}}_{j=1}^{M}n(b_k,w_j)};$ 其中：P(z_t|b_k)表示主题z_t在图像块b_k中出现的概率；P(b_k|z_t)表示图像块b_k在主题z_t中出现的概率；P(z_t|b_k,w_j)表示主题z_t在图像块b_k和单词w_j条件下的后验概率；下标d表示测试集图像块的序号，d＝1,2,…，N_test;

步骤6：KNN分类：计算P(z_t|b_k)与P(z_t|a_i)的欧式距离，得到一个N_test×N_training的欧式距离矩阵

然后根据设定的近邻数目K，利用KNN分类器将测试集中所有的图像块分成滑坡和非滑坡两类；

其中：

表示)与P(z_t|a_i)的欧式距离，计算公式为

下标ki表示元素

在Θ中处于第k行、第i列；K为KNN分类器的近邻数目，取值范围为7～15。

Images