CN104715255A - 一种基于sar图像的滑坡信息提取方法 - Google Patents

Abstract

本发明公开了一种基于SAR图像的滑坡信息提取方法，通过对获取的灾前、灾后两幅原始SAR图像进行预处理，得到基准SAR图像，再对基准SAR图像进行变化检测，得到差值图像，再对差值图像进行纹理灰度信息计算，包括熵值、角二阶矩、对比度、均匀性、均值和方差六种纹理灰度特征。利用巴氏距离选择出识别滑坡效果最好的三个特征，融合成一幅图像。最后，结合滑坡在SAR图像中的特点，利用EM算法确定比例系数和阈值，设定阈值范围内的信息为滑坡信息，从而提取出滑坡信息。本发明利用滑坡在SAR图像中的特点，结合SAR图像的纹理灰度特征及EM阈值分类算法改进了对于SAR图像滑坡信息提取的方法，能够从较复杂地形的SAR图像中快速的提取出滑坡信息。

Description

一种基于SAR图像的滑坡信息提取方法

技术领域

本发明属于图像处理技术领域，更为具体地讲，涉及一种基于SAR图像的滑坡信息提取方法。

背景技术

滑坡是全球性的一种地质灾害，具有数量大、分布广、发生频繁等特点。它不仅对人民生命财产的危害性非常大，还造成环境和资源等的严重破坏。随着人为的影响导致大量的滑坡发生，造成的人员伤亡、经济损失和对环境等的破坏越来越严峻。我国约有70％的山地覆盖，因此滑坡发生密度大，频率高，是世界上发生滑坡较频繁的国家之一，每年给国家和群众造成巨大的损失。尤其是汶川和玉树地震发生后，灾区发生了大量的滑坡，给灾区的人民生命财产带来了极大的威胁和损失。因此，对滑坡灾害进行深入的研究是一项非常有必要而且具有重要社会价值的工作。

如何快速有效地对滑坡进行监测评估工作，对于及时进行抢险救援、减少损失是非常重要的。普通的滑坡监测技术大多依赖于进入滑坡规模、面积等进行调查分析。但是滑坡发生时往往当地的道路被损毁，短时间内很难进入灾害现场对滑坡进行调查，而且很多时候画皮具有分布广泛的特点，传统的现场勘测的方法很难再短时间内对滑坡信息进行提取。因此新的对滑坡进行快速准确提取的方法得到发展。随着卫星遥感技术的发展，可以获取的卫星影像分辨率越来越高，且卫星遥感技术具有快速、准确覆盖范围大等特点，在这种条件下，通过遥感手段对滑坡信息进行提取的技术越来越受到中外学者的关注。然而影响滑坡的因素比较复杂，包括地形地貌、地层岩性、地质构造、水文地质条件和植被等，任何因素都可能引起改变斜坡外形的变化，导致岩土性质恶化。另外，滑坡灾害发生后往往伴随阴雨云雾等恶劣天气，给光学遥感观测带来了极大的不便，因此合成孔径雷达(Synthetic Aperture Radar，SAR)以它可以全天时、全天候的获取数据且不受云雨雾等自然条件影响的优势在对滑坡信息进行提取中得到高度的重视。

然而由于SAR影像中阴影、叠掩等特点，现阶段国内外对于基于SAR图像提取地形复杂区域的滑坡信息还有很多不足，现有的信息提取研究也主要集中于灾害信息提取，不能较明确的区分滑坡和非滑坡区域，如何利用SAR图像提取地形复杂区域的滑坡信息是现有技术需要解决的难题。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于SAR图像的滑坡信息提取方法，以实现地形复杂地区突发性大面积滑坡的提取。

为实现上述发明目的，本发明一种基于SAR图像的滑坡信息提取方法，其特征在于，包括以下步骤：

(1)、获取原始SAR图像

在同一监测点，分别获取滑坡灾害发生前、后两幅原始SAR图像；

(2)、对原始SAR图像进行预处理

对两幅原始SAR图像分别进行辐射校正、滤波、几何校正和配准，得到两幅基准SAR图像；

(3)、获取差值图像

对两幅基准SAR图像进行变化检测，获取到差值图像；

(4)、获取纹理灰度图像

先利用差值图像生成灰度共生矩阵，再利用共生矩阵计算出六种纹理灰度特征图像，最后通过巴氏距离算法计算六种纹理灰度特征图，选取其中三个巴氏距离最大的纹理灰度特征图像，并按照灰度值比例组合成一幅纹理灰度图；

(5)、提取滑坡信息

利用期望极大化算法计算纹理灰度图像，提取出纹理灰度图像中的滑坡区域和非滑坡区域。

本发明的发明目的是这样实现的：

本发明基于SAR图像的滑坡信息提取方法，通过机载或星载文星获取灾前、灾后两幅原始SAR图像，对原始SAR图像进行预处理后，得到基准SAR图像，再对基准SAR图像进行变化检测，得到差值图像，再对差值图像进行纹理灰度信息计算，包括熵值、角二阶矩、对比度、均匀性、均值和方差六种纹理灰度特征。利用巴氏距离(Bhattacharyya)选择出识别滑坡效果最好的三个特征，融合成一幅图像。最后，结合滑坡在SAR图像中的特点，利用EM算法确定比例系数和阈值，设定阈值范围内的信息为滑坡信息，从而提取出滑坡信息。本发明利用滑坡在SAR图像中的特点，结合SAR图像的纹理灰度特征及EM阈值分类算法改进了对于SAR图像滑坡信息提取的方法，能够从较复杂地形的SAR图像中快速的提取出滑坡信息。

附图说明

图1是本发明基于SAR图像的滑坡信息提取方法流程图；

图2是两幅原始SAR图像；

图3是对原始SAR图像进行预处理之后的基准SAR图像；

图4是差值图像；

图5是六种纹理灰度特征图像；

图6是滑坡信息提取效果图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于SAR图像的滑坡信息提取方法流程图。

在本实施例中，如图1所示，本发明一种基于SAR图像的滑坡信息提取方法，包括以下步骤：

S1、获取原始SAR图像

在同一监测点，通过机载或星载卫星分别获取滑坡灾害发生前、后两幅原始SAR图像。

在本实施例中，可以通过Radarsat-2卫星获取C波段的高分辨率SAR图像。如图2所示，两幅原始SAR图像均为分辨率为3m×5m的Radarsat-2C波段图像，其中，图2(a)是滑坡灾害发生前的原始SAR图像，图2(b)是滑坡灾害发生后的原始SAR图像；从图像上可以看出，此处地形非常复杂，包括山地、河流、建筑物、大片林地，存在大量噪声，增加了滑坡的识别与提取难度，但这正是本发明要解决的问题，即克服高分辨率SAR影像复杂地形的大规模快速滑坡信息提取，实现滑坡信息的准确提取。

S2、对原始SAR图像进行预处理

对两幅原始SAR图像分别进行辐射校正、滤波、几何校正和配准，得到两幅基准SAR图像；其中，图3(a)是滑坡灾害发生前的基准SAR图像，图3(b)是滑坡灾害发生后的基准SAR图像；在本实施例中，辐射校正消除了太阳辐射造成的影响；采用Gamma Map滤波方法抑制相干斑噪声，同时能够更好的保持图像的纹理信息；通过几何校正消除图像中存在的各种几何变形，最终得到具有准确位置信息的图像；配准是将两幅图像之间的相同区域对准，即两幅图像相同位置上的两个像素配准。

S3、对基准SAR图像进行变化检测，获取到差值图像

S3.1、计算两幅基准SAR图像中对应像素点的灰度差值

$f\left(\mathrm{\Δ}{x}_{\tilde{i}}\right)=\log (f_1\left(x_{\tilde{i}}\right)/f_2\left(x_{\tilde{i}}\right)),$

其中，为第个像素点的灰度差值，表示滑坡灾害前基准SAR图像中第个像素点的灰度值，表示滑坡灾害后基准SAR图像中第个像素点的灰度值；

S3.2、将所有像素点的灰度差值按照像素点原来位置排列成差值图像。

在本实施例中，差值图像如图4所示，图像中的白色表示变化区域，颜色越深表示变化越小，即差值越小。

S4、获取纹理灰度图像

S4.1、利用差值图像生成灰度共生矩阵

在差值图像中任取一点P，以及偏离它最近点Q，将这两个点组成一个灰度值对为(g₁,g₂)，其中g₁＝f(P)，f(P)为点P的灰度值，g₂＝f(Q)，f(Q)为点Q的灰度值；

令点P在差值图像上移动时，出现各种(g₁,g₂)的值，设差值图像的灰度级为k，则当点P在差值图像上移动时会出现k²种灰度值对(g₁,g₂)；

在整个差值图像中，统计出每一个灰度值对(g₁,g₂)值出现的次数，并排成一个矩阵，在用(g₁,g₂)出现的总次数k²将该矩阵归一化为出现的频率p(g₁,g₂)，最后由p(g₁,g₂)组成k×k的矩阵，即为灰度共生矩阵，并标记为G(i',j')(i'＝1,2,...,k；j'＝1,2,...,k)；

S4.2、利用共生矩阵计算出六种纹理灰度特征图像，如图5所示；

熵值纹理灰度特征图像： $\mathrm{ENT}=-\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}G(i^{\′},j^{\′})\log G(i^{\′},j^{\′});$

角二阶矩纹理灰度特征图像：

对比度度纹理灰度特征图像：

均匀性纹理灰度特征图像：

均值纹理灰度特征图像：

方差纹理灰度特征图像 $\mathrm{VAR}=\frac{1}{k\×k}\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}{(G(i^{\′},j^{\′})-\mathrm{MEA})}^2;$

S4.3、利用巴氏距离算法分别计算六种纹理灰度特征图像，得到每幅纹理灰度特征图像的巴氏距离(BD)；

$\mathrm{BD}=\frac{1}{4}\frac{{({\mathrm{\μ}}_{i_1^{*}}-{\mathrm{\μ}}_{i_2^{*}})}^2}{{{\mathrm{\σ}}_{i_1^{*}}}^2+{{\mathrm{\σ}}_{i_2^{*}}}^2}+\log \left(\frac{{{\mathrm{\σ}}_{i_1^{*}}}^2+{{\mathrm{\σ}}_{i_2^{*}}}^2}{{2\mathrm{\σ}}_{i_1^{*}}{\mathrm{\σ}}_{i_2^{*}}}\right)$

其中，i^*＝1，2，3，4，5，6，和分别为第i^*幅纹理灰度特征图像上滑坡区域和非滑坡区域的均值，和分别表示第i^*幅纹理灰度特征图像上滑坡区域和非滑坡区域的标准差；均值和标准差均可以通过预先训练得到，即：选取已知滑坡和非滑坡区域的图像，分别训练出6幅纹理灰度特征图像中对应位置的滑坡区域和非滑坡区域的均值及标准差

在本实施例中，对均值和标准差进行样本训练时，在第i^*幅纹理灰度特征图像上的滑坡区域和非滑坡区域选择3×3大小的窗口，进行如下计算：

${\mathrm{\μ}}_{i_1^{*}}=\frac{1}{3\×3}\underset{x=1}{\overset{3}{\mathrm{\Σ}}}\underset{y=1}{\overset{3}{\mathrm{\Σ}}}{f}_{i^{*}}(x,y)$

${{\mathrm{\σ}}_{i_1^{*}}}^2=\frac{1}{3^2}\underset{x=1}{\overset{3}{\mathrm{\Σ}}}\underset{y=1}{\overset{3}{\mathrm{\Σ}}}{(f_{i^{*}}(x,y)-{\mathrm{\μ}}_{i_1^{*}})}^2$

其中，f_i*(x,y)表示第i^*幅纹理灰度特征图像中坐标为(x,y)的像素点的灰度值；计算同上，在此不再赘述；

S4.4、选取步骤S4.3中BD值最大的三幅纹理灰度特征图像，并按照灰度值比例组合成一幅纹理灰度图。

本实施例中，根据BD值的大小，选出熵值、均值以及角二阶矩三幅纹理灰度特征图像，再根据BD值计算出灰度比例a_i；

$a_i=\frac{{\mathrm{BD}}_i}{{\mathrm{BD}}_1+{\mathrm{BD}}_2+{\mathrm{BD}}_3},$

其中，a_i表示比例系数，BD_i表示熵值，均值或是角二阶矩三幅纹理灰度特征的BD值，i＝1,2,3；

纹理灰度图中的每个像素点均按照该比例系数进行组合，则组合后的纹理灰度图的像素灰度值F(x)＝a₁f_ent(x)+a₂f_mea(x)+a₃f_asm(x)，其中，f_ent(x),f_mea(x),f_asm(x)表示熵值、均值以及角二阶矩三幅纹理灰度特征图像的像素灰度值。

S5、提取滑坡信息

S5.1、计算分类阀值T

$({{\mathrm{\σ}}_n}^2-{{\mathrm{\σ}}_c}^2)T^2+2({\mathrm{\μ}}_n{{\mathrm{\σ}}_c}^2-{\mathrm{\μ}}_c{{\mathrm{\σ}}_n}^2)T+{{\mathrm{\μ}}_c}^2-{{\mathrm{\μ}}_n}^2{{\mathrm{\σ}}_c}^2-{{2\mathrm{\σ}}_n}^2{{\mathrm{\σ}}_c}^2\ln \left(\frac{{\mathrm{\σ}}_c{\mathrm{\α}}_n}{{\mathrm{\σ}}_n{\mathrm{\α}}_c}\right)=0$

其中，α_n和α_c是比例系数；μ_n和μ_c分别为纹理灰度图像中滑坡和非滑坡区域的均值，σ_n和σ_c分别为纹理灰度图像中滑坡和非滑坡区的标准差，均可以通过训练样本纹理灰度图得到，其计算方法与步骤S4.3中相同；

在本实施例中，比例系数的取值为：α_n＝0.65，α_c＝0.35；分类阀值T为一元二次方程，因此样本训练后可以得到两个分类阀值T，即：T＝1.75和T＝1.95；

S5.2、根据分类阈值T来区分滑坡区域和非滑坡区域

将纹理灰度图中像素点的灰度值与分类阈值T进行比较，如果纹理灰度图中像素点的灰度值在1.75～1.95范围内，则该像素点视为滑坡，否则视为非滑坡。

在本实施例中，通过上述方法对原始SAR图像经过处理后，可以从中提取到滑坡信息，如图6所示，图像中的白色区域即为滑坡区域。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (4)

1.一种基于SAR图像的滑坡信息提取方法，其特征在于，包括以下步骤：

(1)、获取原始SAR图像

在同一监测点，分别获取滑坡灾害发生前、后两幅原始SAR图像；

(2)、对原始SAR图像进行预处理

对两幅原始SAR图像分别进行辐射校正、滤波、几何校正和配准，得到两幅基准SAR图像；

(3)、获取差值图像

对两幅基准SAR图像进行变化检测，获取到差值图像；

(4)、获取纹理灰度图像

先利用差值图像生成灰度共生矩阵，再利用共生矩阵计算出六种纹理灰度特征图像，最后通过巴氏距离算法计算六种纹理灰度特征图，选取其中三个巴氏距离最大的纹理灰度特征图像，并按照灰度值比例组合成一幅纹理灰度图；

(5)、提取滑坡信息

利用期望极大化算法计算纹理灰度图像，提取出纹理灰度图像中的滑坡区域和非滑坡区域。

2.根据权利要求1所述的基于SAR图像的滑坡信息提取方法，其特征在于，所述的步骤(3)中，获取差值图像的方法为：

(2.1)、计算两幅基准SAR图像中对应像素点的灰度差值

$f\left(\mathrm{\Δ}{x}_{\tilde{i}}\right)=\log (f_1\left(x_{\tilde{i}}\right)/f_2\left(x_{\tilde{i}}\right)),$

其中，为第个像素点的灰度差值，表示滑坡灾害前基准SAR图像中第i个像素点的灰度值，表示滑坡灾害后基准SAR图像中第个像素点的灰度值；

(2.2)、将所有像素点的灰度差值按照像素点原来位置排列成差值图像。

3.根据权利要求1所述的基于SAR图像的滑坡信息提取方法，其特征在于，所述的步骤(4)中，获取纹理灰度图像的方法为：

(3.1)、利用差值图像生成灰度共生矩阵

在差值图像中任取一点P，以及偏离它点Q，将这两个点组成一个灰度值对为(g₁,g₂)，其中g₁＝f(P)，f(P)为点P的灰度值，g₂＝f(Q)，f(Q)为点Q的灰度值；

令点P在差值图像上移动时，出现各种(g₁,g₂)的值，设差值图像的灰度级为k，则当点P在差值图像上移动时会出现k²种灰度值对(g₁,g₂)；

在整个差值图像中，统计出每一个灰度值对(g₁,g₂)值出现的次数，并排成一个矩阵，在用(g₁,g₂)出现的总次数k²将该矩阵归一化为出现的频率p(g₁,g₂)，最后由p(g₁,g₂)组成k×k的矩阵，即为灰度共生矩阵，并标记为G(i',j')(i'＝1,2,...,k；j'＝1,2,...,k)。

(3.2)、利用共生矩阵计算出六种纹理灰度特征图像

熵值纹理灰度特征图像： $\mathrm{ENT}=-\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}G(i^{\′},j^{\′})\log G(i^{\′},j^{\′});$

角二阶矩纹理灰度特征图像： $\mathrm{ASM}=\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}G{(i^{\′},j^{\′})}^2;$

对比度度纹理灰度特征图像： $\mathrm{CON}=\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}{(i^{\′}-j^{\′})}^{2}G(i^{\′},j^{\′});$

均匀性纹理灰度特征图像： $\mathrm{HOM}=\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}\frac{1}{1+{(i^{\′}-j^{\′})}^2}G(i^{\′},j^{\′});$

均值纹理灰度特征图像： $\mathrm{MEA}=\frac{1}{k\×k}\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}G(i^{\′},j^{\′});$

方差纹理灰度特征图像 $\mathrm{VAR}=\frac{1}{k\×k}\underset{i^{\′},j^{\′}}{\mathrm{\Σ}}{(G(i^{\′},j^{\′})-\mathrm{MEA})}^2;$

(3.3)、利用巴氏距离算法分别计算六种纹理灰度特征图像的巴氏距离

$\mathrm{BD}=\frac{1}{4}\frac{{({\mathrm{\μ}}_{i_1^{*}}-{\mathrm{\μ}}_{i_2^{*}})}^2}{{{\mathrm{\σ}}_{i_1^{*}}}^2+{{\mathrm{\σ}}_{i_2^{*}}}^2}+\log \left(\frac{{{\mathrm{\σ}}_{i_1^{*}}}^2+{{\mathrm{\σ}}_{i_2^{*}}}^2}{2{\mathrm{\σ}}_{i_1^{*}}{\mathrm{\σ}}_{i_2^{*}}}\right)$

其中，i^*＝1，2，3，4，5，6，和分别为第i^*幅纹理灰度特征图像上滑坡区域和非滑坡区域的均值，和分别表示第i^*幅纹理灰度特征图像上滑坡区域和非滑坡区域的标准差；均值和标准差均可以通过预先训练得到，即：选取已知滑坡和非滑坡区域的图像，分别训练出6幅纹理灰度特征图像中对应位置的滑坡区域和非滑坡区域的均值及标准差

(3.4)、选取步骤(3.3)中BD值最大的三幅纹理灰度特征图像，再根据BD值计算出灰度比例a_i；

$a_i=\frac{{\mathrm{BD}}_i}{{\mathrm{BD}}_1+{\mathrm{BD}}_2+{\mathrm{BD}}_3},$

其中，a_i表示比例系数，BD_i表示三幅纹理灰度特征的BD值，i＝1,2,3；

则组合后的纹理灰度图的像素灰度值F(x)＝a₁f_ent(x)+a₂f_mea(x)+a₃f_asm(x)，其中，f_ent(x),f_mea(x),f_asm(x)表示三幅纹理灰度特征图像的像素灰度值。

4.根据权利要求1所述的基于SAR图像的滑坡信息提取方法，其特征在于，所述的步骤(5)中，期望极大化算法计算纹理灰度图像的方法为：

(4.1)、计算分类阀值T

$({{\mathrm{\σ}}_n}^2-{{\mathrm{\σ}}_c}^2)T^2+2({\mathrm{\μ}}_n{{\mathrm{\σ}}_c}^2-{\mathrm{\μ}}_c{{\mathrm{\σ}}_n}^2)T+{{\mathrm{\μ}}_c}^2{{\mathrm{\σ}}_n}^2-{{\mathrm{\μ}}_n}^2{{\mathrm{\σ}}_c}^2-2{{\mathrm{\σ}}_n}^2{{\mathrm{\σ}}_c}^2\ln \left(\frac{{\mathrm{\σ}}_c{\mathrm{\α}}_n}{{\mathrm{\σ}}_n{\mathrm{\α}}_c}\right)=0$

其中，α_n和α_c是比例系数；μ_n和μ_c分别为纹理灰度图像中滑坡和非滑坡区域的均值，σ_n和σ_c分别为纹理灰度图像中滑坡和非滑坡区的标准差，均可以通过训练样本纹理灰度图得到；

(4.2)、根据分类阈值T来区分滑坡区域和非滑坡区域

将纹理灰度图中像素点的灰度值在分类阈值T的范围内的像素点视为滑坡，否则视为非滑坡。