一种基于新差异算子和纹理的SAR图像水体变化检测算法

2020-04-23李玲玉

计算机与现代化 2020年4期

李玲玉,张毅

(中国科学院电子学研究所,北京 100190)

0 引言

合成孔径雷达(Synthetic Aperture Radar, SAR)是一种高分辨率主动式微波成像系统，不受复杂气候等条件的限制，可全天时全天候地获取影像，甚至可透过地表、植被、云层获取其掩盖的信息。且SAR具有成本低、覆盖广的特点，使其在水体变化检测等领域具有独特的优势[1]。利用SAR进行水体变化的检测可优化水体资源配置，完善防洪减灾体系，完善湖泊等生态环境建设，进行水体附近综合管理。各种应用的难点在于如何提高SAR图像水体变化检测精度[2]。

为精确得到SAR图像的变化信息，研究人员对其进行了大量的研究。变化检测方法有直接比较法和分类后比较法。研究人员通常用直接比较法，首先利用变化信息构建差异图，然后针对差异图进行分析，区分出变化类与非变化类。

传统的差异图生成法主要包含如下3类：1)基于像素的算子，如Dekker[3]提出的比值法(LR)、Celik[4]提出的对数比值法(LR-D)等；2)基于邻域的算子，如Inglada等[5]提出的均值比算子(MR)、Gong等[6]提出的邻域比值法(NR)，以及Xiong等[7]提出的最大似然比算子(LLR)等；3)几种方法的结合，因为一般单一的差异图生成算子不能很好地反映出图像之间的差异。如Hou等[8]提出的LR与高斯对数比结合的方法，Zheng等[9]提出的差值法与比值法结合的方法，Yan等[10]提出的LR与对数均值比结合的方法。上述方法均没有结合单像素算子和邻域算子的优点。

差异图分析法包括阈值法、主动轮廓法和聚类法等。阈值法能快速区分变化区域和非变化区域，然而当变化类和不变类之间存在较强重叠或无法准确建模时，使用阈值法提取结果较差。主动轮廓法可以在一定程度上降低相干斑噪声的影响，但在初始轮廓线与真实轮廓线相距较远时误差较大。聚类法，如Bezdek等[11]提出的模糊C均值(FCM)方法，它较依赖于人工给出的参数，没有考虑邻域信息；Krinidis等[12]提出的局部信息模糊C均值聚类(FLICM)方法，解决了上述问题，但仍容易被孤立噪声点影响；Gong等[13]提出的基于改进MRF能量函数的FCM(MRFFCM)算法，降低了散斑的影响，但不能保证噪声与细节的平衡。而空间、纹理等信息可降低孤立噪声的影响，引入更多的SAR图像特征。

传统方法精度有限，为了提高水体变化检测精度，本文结合单像素算子和邻域算子的优点，将空间、纹理等信息与适应调节抑制噪声和细节平衡的模糊因子结合运用到算法中，提出一种基于新差异算子和纹理的SAR图像变化检测方法，该方法分为3步:

1)根据SAR图像的特征，结合LR和LLR这2种算子，提出一种新的差异算子。

2)提出一种新的FLICM方法，将目标分为3类。

3)根据差异图像分割阈值将过渡区域再次分类，得到最终的变化检测结果。

本文算法的具体流程如图1所示。

图1 本文的SAR图像水体变化检测算法流程

1 本文的SAR图像水体变化检测算法

1.1 差异图像的生成

针对SAR图像的噪声特征，比值法的比率运算符被广泛使用于SAR图像的变化检测中，它通常以对数刻度表示，可将SAR图像的乘性噪声变为加性噪声，增强了SAR图像对辐射误差的鲁棒性。本文基于比值法的LR算子与LLR算子，分析其特性，提出将两者结合的变化检测新算子。

LR是应用非常广泛的比率检测算子，可以很好地反映单像素差，但其忽视了邻域信息，LR的计算公式为：

(1)

其中，X1、X2为变化前后的图像。DLR最小值为1，代表不变的像素。而变化像素的值通常较大，因此采用对数运算符将其缩小。采用对数运算符后，最小值对应为0，采用绝对值运算符把正负像素视为变化像素。因此LR变为：

(2)

水体一般为镜面散射，因此水体在SAR图像中近似为均匀区域，伽马分布能很好地拟合水体的SAR图像。LLR算子假设SAR图像服从伽玛分布，并构造了似然比检验[7]，采用3×3的矩形窗代表局部信息，既能反映邻域差异，又具有较好的抗噪声能力。仍然将对数运算符及绝对值运算符运用到简化后的LLR算子中计算：

(3)

其中，Ω为图像像素(m,n)在3×3窗口内的像素集。上述的2个算子，不变像素的值在0附近形成了一个尖峰，变化的像素形成一条长尾。最后将差异图像归一化到[0,255]。本文利用1997年7月和8月加拿大渥太华的2幅RADARSAT SAR图像作为数据集1来展示本文算法的优势，该图像分辨率为12 m，大小为350×290，如图2(a)和图2(b)所示，图2(a)为洪水后的SAR图像，图2(b)为洪水前的SAR图像，参考差异图像如图2(c)所示，DLR、DLLR归一化直方图如图3(a)和图3(b)所示。

(a) 7月原始图像

(b) 8月原始图像

(a) DLR归一化直方图

(b) DLLR归一化直方图

本文通过将LR算子与LLR算子相乘构成SAR图像水体变化检测新算子。在单像素差和领域差同时都很强时，新算子的值将会变大，在单像素差和领域差同时都很弱时，新算子的值将会变小。而在有噪声造成单像素差弱但邻域差强，或单像素差强但邻域差弱时，新算子的值介于两者之间。最终的差异图生成新算子为：

DI=DLR×DLLR

(4)

上述DLR、DLLR、DI这3种算子生成的差异图像如图4(a)～图4(c)所示。

(a) DLR差异图

(b) DLLR差异图

DI归一化直方图如图5(a)所示。本节利用差异图像的相邻直方图比值进行阈值分割，得到初始分割图像，如图5(b)所示。

(a) DI归一化直方图

(b) 差异图像DI的相邻直方图比值图

差异矩阵直方图中，峰值表示无变化区域，震荡部分为变化区域，可以得到变化与非变化的检测门限在其中间，具体值可用相邻直方图比值图得出。差异图像的相邻直方图比值计算公式为：

(5)

其中，hDI为差异图像直方图，i∈[i0,255)，i0为hDI最大值对应的灰度。自动选择比率小于1的点作为过渡点的第一个点为阈值,即：

k={i|η(i)1}

(6)

分析DLR、DLLR、DI算子的直方图，所提出的DI的直方图在不变部分显示出更尖锐的峰值，在变化部分显示出更长的尾部。该算子放大了变化区域与不变区域的特征，同时降低了噪声效果。

分析差异图像的相邻直方图比值图，如图6(a)所示，得到差异图的初始分割。DLLR与DI的初始分割如图6(b)、图6(c)所示。其中，DLLR的阈值为文献[14]中所提出的阈值。与参考图对比，可以得到DLLR的PCC为0.9611，kappa系数为0.8392；DI的PCC为0.9746，kappa系数为0.9038。对差异图进行阈值分割，DI的变化检测结果好于DLLR。

(a) 相邻直方图比值图

(b) DLLR差异图初始分割

1.2 基于纹理的差异图获取

现今，有许多纹理分析技术，如灰度共生矩阵(GLCM)、傅里叶功率谱、马尔科夫随机场(MRF)和Gabor滤波器[14～17]等。其中，Gabor滤波器在分类等问题上应用广泛，文献[18]验证了Gabor与其他纹理特征相比，可以在计算成本低的情况下达到同样的性能。二维高斯滤波器公式如下：

(7)

利用尺度和方向的不同，将图像与二维Gabor核进行卷积，得到了DI(x，y)的N组Gabor响应。本文选择尺度和方向个数分别为4、6。Gabor滤波器及其对应的响应如图7(a)和图7(b)所示。

(a) Gabor滤波器

(b) 对应的响应

将每组响应变为一列，得到Gabor差分图像矩阵为：

Gabor=(gbi1,gbi2,…,gbii,…,gbin)

(8)

其中，n=24，gbi为每组Gabor滤波器对应的响应。除纹理信息外，本节还加入了2个空间特征，它们为归一化后的图像x坐标和y坐标。

因此，差分图像的像素可以由这26维特征向量表示，为了减少冗余，采用PCA进行降维，留下特征值超过最大特征值的5%的特征向量。最后，采用公式(9)将其变为一维进行分析[19]。

(9)

将其整形为原始图像大小，得到的结果如图7(c)所示。

1.3 改进的FLICM方法

改进的FLICM方法称为结合纹理的FLICM方法(FLICM_texture)，该方法分为3个步骤：估计初始聚类中心、FLICM_texture的聚类和阈值分割的再次分类。

1.3.1 初始聚类中心的获得

由于变化类与非变化类之间有一定的过渡区域，该区域有灰度级低的植被，或者受到了噪声影响。因此本文将差异图分为3类：变化区域、非变化区域和过渡区域。

由1.1节中提出的DI差异图像的相邻直方图比值阈值法提取出核心变化区域，即：

(10)

其中，k为1.1节中对应的阈值。经过对比分析，该阈值分割的精度较高，因此可利用它得到3个初始聚类中心：

(11)

1.3.2 FLICM_texture方法

2010年，Stelios在FCM目标函数中引入一个新的因子，提高了对噪声和异常值的鲁棒性，不需提前计算参数，不再需要人工选择参数。该模糊因子计算公式如下：

(12)

其中，k为聚类中心，j为第i个像素局部窗口的像素集合，m为每个模糊隶属度的加权指数。dij为像素i与j之间的空间欧氏距离，ukj为第j个像素在第k个聚类中心的隶属度，vk为聚类k的中心。可以看出，因子Gki没有使用任何控制图像噪声和图像细节之间平衡的参数，且该因子包含了局部信息。

将局部空间和灰度信息应用于目标函数中，并将纹理信息添加到目标函数中，得到最后的目标函数如下式：

(13)

其中，yi、wk分别为纹理图像中像素与聚类中心的值。使上式最小化的隶属度及聚类中心表达式为：

(14)

(15)

FLICM_texture算法流程如下：

1)设置类别数c、模糊参数m和停止条件ε。

2)将1.3.1节得到的初始聚类中心设为该算法的初始聚类中心，计算纹理聚类中心，设置循环数b=0。

3)用式(14)计算隶属度。

4)用式(11)修正聚类中心，计算新的纹理聚类中心。

5)当max{u(b)-u(b+1)}<ε时停止，否则，令b=b+1，进入步骤3。

当算法收敛后，进行去模糊化处理，将满足下面条件的像素i划分给第k类。

Ci=argk{max {uki}},k=1,2,…，c

(16)

由式(16)可以看出，该方法并未增加FLICM的时间复杂度。

生成纹理时，二维Gabor滤波器中的行数m的选取过程：利用FLICM_texture算法检测图像，当m=9时，PCC为0.9539，kappa系数为0.8047，FA为0.0083，MA为0.2855。当m=39时，PCC为0.9579，kappa系数为0.8238，FA为0.0083，MA为0.2599。m=39时，其PCC、kappa系数比m=9时高，FA、MA比m=9时低，显然，m=39时的检测效果好于m=9，因此令二维Gabor滤波器中的行数为39，其中PCC、kappa、FA、MA系数说明见第2章，生成的纹理图像如图8所示。

(a) m=9时纹理图像

(b) m=39时纹理图像图8 纹理图像

在一切参数相同的条件下，对FLICM和FLICM_texture将差异图分为2类，得到的检测结果如图9(a)和图9(b)所示。在本文的方法中，去除FLICM_texture的纹理及完整的本文方法得到的结果如图9(c)和图9(d)所示。

(a) FLICM结果 (b) FLICM_texture结果

对于图9(a)和图9(b)，与参考图像对比，可以得到FLICM的PCC为0.9499，kappa系数为0.7843，MA为0.3152；FLICM_texture的PCC为0.9580，kappa系数为0.8238，MA为0.2599。PCC增长了约0.85%，kappa系数增长了约5.0%，MA下降了约17.5%。对于图9(c)和图9(d)，图9(c)的PCC为0.9761，kappa系数为0.9067，MA为0.1226；图9(d)的PCC为0.9800，kappa系数为0.9203，MA为0.0776。这验证了FLICM_texture方法的有效性。

若只是使用FLICM_texture方法将差异图等分为2类则精度较低，因此用FLICM_texture方法将其分为3类：水体、背景和过渡区域，使用阈值法将过渡区域再次划分为2类。

1.3.3 阈值分类

由于变化类与非变化类之间有一定的过渡区域，该区域有灰度级低的植被，或者受到了噪声影响。本文采用阈值法对其进行二次分类。判断公式为：

(17)

其中，Ωg为差分图像过渡区域的像素集合，Dm为该区域的像素，k为1.1节中利用相邻直方图比值得到的阈值。

2 定量分析变化精度

将本文提出的方法得到的变化检测结果与FCM、模糊局部C均值法(FLICM)、基于马尔科夫随机场的ICM算法(MRF&ICM)[20]、阈值法做对比，利用以下4个定量指标分析算法的有效性：1)虚警(FA)表示错误检测到的未变化像素的百分比;2)漏警(MA)表示未检测到的已变化像素的百分比;3)kappa系数表示2幅图像的一致性程度;4)正确分类率(PCC)表示百分比正确分类[10]。

表1 各种方法在加拿大渥太华上空数据集的变化检测结果

1)数据集1为渥太华数据集，它分别是RADARSAT在1997年5月和8月在加拿大渥太华上空拍摄的SAR图像的切片(290×350)，它们为曾经遭受洪水侵袭的地区。各种方法在渥太华数据集的变化检测精度如表1所示，各种方法在渥太华数据集的变化检测结果图如图10所示。

(a) FCM+DLR (b) FLICM+DLR

(e) FLICM+DLLR (f) 阈值+DLLR

(g) MRF&ICM+DLLR (h) DI+阈值

(i) 本文方法 (j) 参考图

2)数据集2为伯尔尼数据集，它分别是ERS-2在1999年4月和5月在瑞士伯尔尼市上空拍摄的SAR图像的切片(301×301)，如图11(a)和图11(b)所示，DI差异图如图11(c)所示。在这2个日期之间，部分城市、伯尔尼机场遭到了阿雷河的淹没。各种方法在伯尔尼数据集的变化检测精度如表2所示，各种方法在伯尔尼数据集的变化检测结果如图12所示。

(a) 4月原始图像

(b) 5月原始图像

表2 各种方法在伯尔尼数据集的变化检测结果

(a) FCM+DLR (b) FLICM+DLR

(e) FLICM+DLLR (f) 阈值+DLLR

(g) MRF&ICM+DLLR (h) DI+阈值

(i) 本文方法 (j) 参考图

3)数据集3为印度金奈数据集，它分别是哨兵在2018年8月4日和2019年7月30日于金奈上空拍摄的切片(1373×1927)，如图13(a)和图13(b)所示，DI差异图如图13(c)所示。其中金奈2019年遭受严重干旱，造成湖泊等水体大面积缩减[21]。各种方法在金奈数据集的变化检测精度如表3所示，各种方法在金奈数据集的变化检测结果如图14所示。比较水体变化面积可以得到，金奈附近水体相比于去年减少了73.24%，其中金奈最大的湖Puzhal湖面积减少了近83.88%。

(a) 2018年8月原始图像

(b) 2019年7月原始图像

表3 各种方法在金奈数据集的变化检测结果

(a) FCM+DLR (b) FLICM+DLR

(e) FLICM+DLLR (f) 阈值+DLLR

(g) MRF&ICM+DLLR (h) DI+阈值

(i) 本文方法 (j) 参考图

综上可以看出，由于差异图中引入局部平均，差异图像中一些细小的变化丢失了，但变化丢失相对较少。比较由DLLR和DI得到阈值分割后的结果的PCC与kappa系数，发现本文提出的新算子具有良好的性能。将本文方法与其他方法相比，本文方法提高了检测率，降低了虚警率，提高了kappa系数，因此本文方法具有较高的鲁棒性。

3 结束语

对SAR图像水体变化进行检测具有完善防洪减灾体系、进行水体附近综合管理等作用。本文针对多时相SAR变化检测中孤立噪声点、需人工输入部分参数、信息应用不全等问题，提出了基于新差异算子和纹理的变化检测算法。

将本文方法应用于加拿大渥太华上空、瑞士伯尔尼市上空、印度金奈上空的SAR水体变化检测，可以看出，本文方法在3个地区的PCC和kappa系数上均优于FCM、FLICM、MRF&ICM、阈值法，因此本文方法提高了SAR图像水体变化检测精度。但是，本文提出的新差异算子对低后向散射区域较为敏感，如果SAR图像中有大的低后向散射区域，虚警率可能会变大。下一步将针对低后向散射区域进行处理，使其更具有通用性。