利用纹理融合与广义高斯模型的高分辨率SAR影像变化检测

刘本强，赵 争，盛玉婷，张忠芳

(1.山东科技大学 测绘科学与工程学院，山东 青岛 266590；2.中国测绘科学研究院，北京 100830；3.北京航天世景信息技术有限公司，北京 100089)

合成孔径雷达(Synthetic Aperture Radar, SAR)具有全天时、全天候的监测能力。利用SAR影像进行变化检测是当前研究的热点问题。变化检测是通过对同一地区不同时期的遥感影像比对分析，获取地物的变化信息，在灾害监测[1-2]、土地利用[3]、城市化发展[4]等方面应用广泛。

近年来，Cosmo-SkyMed、TerraSAR-X、Radarsat-2和ALOS-2等商业合成孔径雷达卫星的发射，极大满足获取高分辨率SAR影像的需求。不同于中低分辨率SAR图像，高分辨率的SAR图像具有更加丰富的纹理、几何和细节等特征，得到研究人员的广泛关注[5-7]。文献[8-9]针对高分辨率的机载SAR数据，综合考虑灰度共生矩阵纹理、LBP纹理、Gabor纹理等纹理特征，描述震后建筑物的损坏等级。在变化检测中引入纹理信息不仅能够抑制斑点噪声，而且在检测性能上也会有所提高。陈志鹏等[10]采用纹理信息代替灰度信息，使用差值法获取城区变化检测结果，并比较不同纹理的检测性能。李坤等[11]针对地震灾害导致的地表变化问题，分析SAR图像灰度和纹理特征的敏感程度，结合主成分分析和K均值聚类技术提取震前震后的变化信息。Wang等[12]利用方向滤波和分形纹理进行SAR影像变化检测。赵凌君等[13]构建9种常见的纹理特征，然后利用Bhattacharyya距离进行特征选择，通过K近邻分类器实现建筑区的检测。

广义高斯模型(Generalized Gaussian Model，GGM)能够较好的描述SAR图像统计特性，文献[14-16]均利用广义高斯模型描述差异图的变化类像元和未变化类像元的概率分布，并基于KI阈值法实现多时相单极化SAR影像的非监督变化检测。Zhao等[17-18]结合统计模型的方法分别利用Kittler-Illingworth(KI)算法和改进的KI阈值算法对相应的差异影像进行阈值提取，实现多时相全极化SAR数据的非监督变化检测，获得较好的变化检测结果。

本文为了充分利用高分辨率SAR影像的灰度信息和纹理信息，提高变化检测精度，提出一种纹理信息融合与广义高斯模型相结合的SAR影像变化检测方法。该方法以灰度共生矩阵(gray level co-occurrence matrix, GLCM)计算的纹理特征构造纹理差值差异图，以对数比值运算构造灰度差异图，以小波变换为工具融合两类差异图。然后通过广义高斯模型对融合后的差异图统计建模，估计变化类像元和未变化类像元的概率分布，利用KI阈值准则获取最佳分割阈值，得到最终变化检测结果。

1 算法介绍

1.1 灰度共生矩阵与纹理特征

纹理描述像素的邻域灰度空间分布规律，是特征提取、影像分类、目标识别的重要手段。常用的纹理提取方法有统计方法、几何方法、结构方法、模型方法和信号处理方法。灰度共生矩阵是最常用的一种基于统计的纹理提取方法，由Haralick等[19]提出，定义14种基于GLCM的纹理特征。

灰度共生矩阵定义为

P(i,j,d,θ)=#{[(x,y),(x+Δx,y+Δy)|d,

θ,f(x,y)=i,f(x+Δx,y+Δy)=j]}.

(1)

归一化方式为

(2)

式中：#{•}表示大括号中成立的像元对数；(x,y)为像元坐标，且x=0,1,…,M-1,y=0,1,…,N-1;M，N为影像的高度和宽度；i,j=0,1,…,L-1表示灰度级；Δx，Δy为像素坐标偏移量；d为步长，θ为方向，一般取0°、45°、90°、135°；f(•)表示像素坐标的灰度值。

根据灰度共生矩阵可以定义如下纹理特征：角二阶矩(angular second moment, ASM)、熵(entropy)、逆差矩(inverse difference moment, IDM)。

角二阶矩：

(3)

熵：

(4)

逆差矩：

(5)

1.2 基于小波变换的差异图融合

小波变换作为一种数学工具，在图像融合处理领域有着很重要的地位。Gong等[20]和慕彩虹等[21]分别利用离散小波变换(discrete wavelet transform, DWT)和离散平稳小波变换(stationary wavelet transform, SWT)将对数比值差异图和均值比差异图融合，取得显著的效果。

首先利用式(6)和式(7)构造纹理差异图和灰度差异图，其中Ti(x,y)为i时刻SAR影像的纹理特征影像，Ii(x,y)为i时刻SAR强度影像，i=1,2表示两个时相，x和y表示影像的行列。

TDI(x,y)=|T1(x,y)-T2(x,y)|,

(6)

LR(x,y)=log(I1(x,y)/I2(x,y)).

(7)

然后，对两种差异图分别进行3层SWT小波分解，每一层分解分别得到一个低频系数矩阵和3个高频系数矩阵。对于低频系数，采用式(8)的加权平均融合规则，对于高频系数，采用式(9)的局部能量最小融合规则。

(8)

(9)

Ek(x,y)=∑Ω(Hk(x,y))2

(10)

表示以(x,y)为中心的局部能量，Ω为计算局部能量时窗口大小，本文中取3×3。

最后，将融合后的低频系数和高频系数执行SWT逆变换，得到融合后的差异影像。

1.3 广义高斯模型下的KI阈值(GGKI)

对于融合后差异图Fswt(x,y)上存在着的变化类ωc和未变化类ωu，根据贝叶斯决策理论，差异图上的任意一个像元Fswt(x,y)可由式(11)表示，其中Xf表示融合后差异图的像元灰度值。

p(Xf)=P(ωu)p(Xf|ωu)+P(ωc)p(Xf|ωc).

(11)

在广义高斯模型的假设下，两类像元的条件概率密度可以表示为

p(Xf|ωi)=aiexp(-(bi|Xf|-mi)βi),

i∈{u,c},

(12)

其中

(13)

参数mi，σi和βi分别表示广义高斯分布的均值、标准差和形状控制参数，Γ(•)为Gamma函数，表达式为

(14)

KI阈值选择法是一种以最小错误率贝叶斯理论为基础的阈值选择方法。图像直方图h(Xf)的最佳分割阈值T由一个能够描述整幅图像的准则函数J(T)确定。

(15)

式中c(Xf,T)为代价函数，其表达式为

(16)

其中，p(ωi|Xf,T)(i=u,c)为阈值T下的未变化类和变化类的后验概率密度。使得图像分类误差最小的阈值为最佳阈值，即T*为最佳阈值时需使得准则函数J(T)最小。

(17)

根据贝叶斯理论，后验概率密度函数可以由先验概率和条件概率密度表示为

p(ωi|Xf,T)=

(18)

综上所述，可以推导出准则函数的表达式为

H(Ω,T)-(Pu(T)lnau(T)+Pc(T)lnac(T)).

(19)

式中H(Ω,T)表示类别集合Ω∈{ωu,ωc}的熵，表达式为

(20)

对于广义高斯模型中相关参数的估计，可以采用

(21)

2 实验与分析

2.1 实验数据

实验选用3 m分辨率的TerraSAR-X数据，覆盖广州市和贵阳市某区域。为方便实验研究，本文选取其中1 500 m×1 300 m和1 000 m×1 000 m的区域进行实验，见图1和图2。其中图1(c)和图2(c)为人工标定的变化参考图，用来定量评价检测效果。

图1 广州市TerraSAR-X数据

图2 贵阳市TerraSAR-X数据

2.2 实验方法

图3是本文方法的流程图。主要包括3个部分内容：

1)对实验影像做预处理(辐射校正、几何校正、影像配准、影像滤波)，其中影像滤波采用3×3 Lee滤波，以削弱斑点噪声的影响；

2)GLCM计算中，纹理窗口选择11×11，d=1，θ分别取0°、45°、90°、135°，计算后的4个方向纹理特征取均值，再参照式(6)构造不同纹理特征的纹理差异图；

3)按照式(7)构造对数比值差异图，并与不同纹理特征的纹理差值差异图做基于SWT变换的融合处理；

4)采用广义高斯分布对融合后的差异图建模，利用KI阈值准则获得最佳分割阈值，提取变化信息并评价检测精度。

图3 本文方法流程

2.3 实验结果及分析

将对数比值图与纹理差值差异图经SWT变换，得到的融合差异图见图4和图5，其分别是对数比值图与角二阶距纹理(ASM)、熵纹理(ENT)、逆差距纹理(IDM)的融合结果、对数比值差异图。可以看出，在对数比值差异图中，背景中的噪点特别多，变化区域也“破破碎碎”。而融合纹理信息的差异图相对平滑，噪点明显减少，且变化区域更加完整。这是因为噪点为频率域的高频部分，而高频系数采用局部能量最小原则进行融合，对噪点有一定的抑制作用。同时，每一层的SWT分解也会降低噪声的影响。

采用广义高斯分布对融合差异图建模，拟合变化类像素和未变化类像素的概率分布，利用KI阈值准则获取最佳分割阈值提取变化信息，得到变化二值图。图6和图7分别是广义高斯模型对两组数据集融合差异图的拟合结果。可以看出对于融合后的差异图，广义高斯模型均能较好的拟合其变化类像元和未变化类像元的概率密度函数，特别注意的是，广义高斯模型对于融合逆差矩(IDM)纹理的差异图拟合效果最好。

图4 广州数据集差异图

图5 贵阳数据集差异图

图6 广州数据集广义高斯模型拟合结果

图7 贵阳数据集广义高斯模型拟合结果

基于广义高斯模型KI阈值算法(GGKI)的变化检测结果见图8和图9。可以直观地看出融合3种纹理信息均能够较好的提取出变化区域，但融合ASM纹理检测结果的背景中含有较多的噪点。

为了定量评价检测精度，验证所提方法的有效性，分别对不同方法得到的结果计算检测指标，进行精度评定，结果见表1。

错检率：实际未变化像元被错分为变化像元的个数与实际未变化像元总数的比率。

漏检率：实际变化像元被错分为未变化像元的个数与实际变化像元总数的比率。

总正确率：两类检测正确的像元和与总像元的比率。

Kappa系数：其值越大，检测精度越高。

图8 广州数据集变化检测结果

图9 贵阳数据集变化检测结果

数据集方法错检率/%漏检率/%总正确率/%Kappa系数广州融合ASM+GGKI2.3720.1096.730.6941融合ENT+GGKI1.8324.9197.000.7011融合IDM+GGKI1.7035.4897.300.7068贵阳融合ASM+GGKI4.5720.7292.940.7339融合ENT+GGKI3.3318.4394.340.7829融合IDM+GGKI2.5718.6894.950.8025

从表1中可以看出，针对高分辨率SAR影像：①融合纹理信息与广义高斯模型的变化检测方法其Kappa系数均达到70%左右，能够满足变化检测需求。②不同纹理信息的检测性能不同，在具体应用时可以针对特定需求选择不同的纹理特征参与融合。③融合IDM纹理信息的检测方法具有最低的错检率，并且其总正确率和Kappa系数最高，检测效果最好。

针对水域界限，本文算法能够较好的提取其变化信息。水体多为面状和条带状，在SAR影像上呈暗色调，灰度值较小，与周围其他地物具有很高的对比度，容易区分与识别。本文在计算GLCM时窗口大小取11×11，由于水体区域的灰度值变化幅度较小，所以计算出的纹理特征更加趋于同质，更加便于区分水体与非水体，通过纹理差值就能发现变化的区域。再通过小波变换融合纹理差异图和灰度差异图，进一步增强变化信息，抑制背景中的噪点，提高变化检测性能。

3 结束语

本文充分利用高分辨率SAR影像丰富的纹理信息，以SWT变换为工具，将纹理差异图与对数比值图融合。融合后的差异图有效地利用两类差异图的优势信息，充分反映变化信息的位置与范围，性能更加优异。同时，利用基于广义高斯模型的KI阈值分割方法，提取变化信息，得到变化检测二值图。本文选用两组3 m分辨率的TerraSAR-X影像，对比分析融合不同纹理信息的检测性能，实验结果证明广义高斯模型能够较好地拟合出融合纹理信息差异图的概率分布。同时纹理融合与广义高斯模型相结合的变化检测方法具有较可靠的检测精度，应用于高分辨率SAR影像变化检测可行且有效。

参考文献：

[1] ZHAO J,YANG J,LU Z,et al.A Novel Method of Change Detection in Bi-Temporal PolSAR Data Using a Joint-Classification Classifier Based on a Similarity Measure[J].Remote Sensing,2017,9(8): 846.

[2] 刘莹,陶超,闫培,等.图割能量驱动的高分辨率遥感影像震害损毁建筑物检测[J].测绘学报,2017,46(7):910-917.

[3] 杨育丽,马明国,周兆叶.基于地形梯度的山地流域土地利用格局及时空演变研究-以云南小江流域为例[J].测绘工程,2017,26 (12): 8-16.

[4] SCHLAFFER S,MATGEN P,HOLLAUS M,et al.Flood Detection from Multi-temporal SAR Data Using Harmonic Analysis and Change Detection[J].International Journal of Applied Earth Observation and Geoinformation,2015(38): 15-24.

[5] LEFEBVRE A,CORPETTI T,HUBERT-MOY L.Object-oriented Approach and Texture Analysis for Change Detection in Very High Resolution Images[C]//Geoscience and Remote Sensing Symposium,2008,IV-663-IV-666.

[6] TISON C,NICOLAS J M,TUPIN F,et al.A New Statistical Model for Markovian Classification of Urban Areas in High-resolution SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(10): 2046-2057.

[7] WANG G G,HUANG X T,ZHOU Z M,et al.A New SAR Image Change Detection Algorithm Based on Texture Feature [C]//Synthetic Aperture Radar (APSAR),2011 3rd International Asia-Pacific Conference on.IEEE,2011: 1-4.

[8] SUN W,SHI L,YANG J,et al.BuildingCollapse Assessment in Urban Areas Using Texture Information from Postevent SAR Data[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2016,9(8): 3792-3808.

[9] SHI L,SUN W,YANG J,et al.Building Collapse Assessment by the Use of Postearthquake Chinese VHR Airborne SAR[J].IEEE Geoscience and Remote Sensing Letters,2015,12(10): 2021-2025.

[10] 陈志鹏,邓鹏,种劲松,等.纹理特征在SAR图像变化检测中的应用[J].遥感技术与应用,2002,17(3): 162-166.

[11] 李坤,杨然,王雷光,等.基于敏感特征向量的SAR图像灾害变化检测技术[J].计算机工程与应用,2012,48(3):24-28.

[12] WANG G,HUANG X,ZHOU Z,et al.A New SAR Image Change Detection Algorithm Based on Texture Feature[C]//Synthetic Aperture Radar (APSAR),2011 3rd International Asia-Pacific Conference on.IEEE,2011: 1-4.

[13] 赵凌君,秦玉亮,高贵,等.利用 GLCM 纹理分析的高分辨率 SAR 图像建筑区检测[J].遥感学报,2009 (3): 475-490.

[14] BAZI Y,BRUZZONE L,MELGANI F.An Unsupervised Approach Based on The Generalized Gaussian Model to Automatic Change Detection in Multitemporal SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing,2005,43(4): 874-887.

[15] 袁增云,胡召玲.基于广义高斯模型的多时相 SAR 图像变化检测[J].江苏师范大学学报(自然科学版),2013,31(1): 74-78.

[16] 胡召玲.广义高斯模型及 KI 双阈值法的 SAR 图像非监督变化检测[J].测绘学报,2013 (1): 116-122.

[17] ZHAO J Q,YANG J,LI P X,et al.An Unsupervised Change Detection Based on Test Statistic and KI from Multi-temporal and Full Polarimetric SAR Images[J].ISPRS-International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2016: 611-615.

[18] ZHAO J Q,YANG J,LU Z,et al.An Unsupervised Method of Change Detection in Multi-Temporal PolSAR Data Using a Test Statistic and An Improved K&I Algorithm[J].Applied Sciences,2017,7(12): 1297.

[19] HARALICK R M,SHANMUGAM K.Textural Features for Image Classification[J].IEEE Transactions on systems,man,and cybernetics,1973 (6): 610-621.

[20] GONG M G,ZHOU Z Q,MA J J.Change Detection in Synthetic Aperture Radar Images Based on Image Fusion and Fuzzy Clustering[J].IEEE Transactions on Image Processing,2012,21(4):2141-2151.

[21] 慕彩红,霍利利,刘逸,等.基于小波融合和PCA-核模糊聚类的遥感图像变化检测[J].电子学报,2015,43(7):1375-1381.