合成孔径雷达(synthetic aperture radar,SAR)具有穿透云雨雾霾，全天时、全天候获取地面信息的能力，自20 世纪50年代以来，已被广泛用于地表沉降监测、地形测绘、资源勘探、环境遥感以及军事等领域。近年来随着SAR技术的迅速发展，利用SAR影像进行变化检测[1]逐渐成为研究的热点[2, 3, 4, 5]。
现有SAR影像变化检测方法多使用灰度信息[6, 7, 8, 9, 10]，由于SAR影像中存在明显的斑点噪声[11]，变化检测结果通常具有较高的虚警率和漏检率，此时可以采用图像处理方法减弱斑点噪声的影响[1, 12]。随着SAR影像分辨率的提高，纹理特征作为一种可以减小斑点噪声影响的图像信息在变化检测中得到应用[13, 14, 15]。然而，目前基于纹理的变化检测研究，一般仅选择单一纹理或几个适用于试验区域的纹理特征[14, 15]。不同纹理对同一地物的敏感性不同，由于地物变化的不确定性，仅使用一个或几个有限纹理进行变化检测的普适性不高。由于某一种描述纹理的特征，对差异的描述结果有两种：①可以区分出地物的变化，在计算纹理特征向量的欧氏距离时，这一分量较大；②无法区分出地物变化，在计算纹理特征向量的欧氏距离时，这一分量接近0。本文针对4个方向的灰度共生矩阵(gray-level co-occurrence matrix,GLCM)，分别计算8种具有代表性的纹理特征，组成32维特征向量，用于构造考虑空间邻域信息的差异影像。
阈值自动选取是实现非监督变化检测的关键之一，由于阈值选取中常用的KI算法、EM算法等一般都要假设差异图像中变化类和未变化类服从某种分布[8, 9, 10, 12, 16]，如果分布函数选择不恰当就会得到不合理的阈值。因此，本文使用不依赖分布模型的最大化熵法自动选取阈值[17, 18]，从而实现SAR影像非监督变化检测，方法流程见图 1。
|
| 图 1 本文方法的流程图 Fig. 1 Flowchart of the proposed method |
图像纹理由于形式上的广泛性和多样性，到目前还没有为众人所公认的定义。文献[19]分析了几个具有代表性的定义，总结出两点共识：①纹理不同于灰度和颜色等图像特征，它通过像素及其周围空间邻域的灰度分布来表现，即局部纹理信息；②局部纹理信息不同程度的重复性，即全局纹理信息。
文献[20]将纹理特征按照提取方法归为5类，即统计方法、几何方法、结构方法、模型方法和信号处理方法。SAR是主动成像，灰度信息与地物的后向散射强度相关，虽然具有丰富的纹理信息，但是不具有常规纹理不规则的局部模式和简单的重复周期，这使得SAR影像纹理通常只具有统计学的意义。文献[21]在利用陆地卫星图像研究美国加利福尼亚海岸带的土地利用问题时,首次提出GLCM方法，该统计方法对纹理的细节性和随机性描述较好，具有较强的适应能力与稳健性，是分析SAR图像纹理的有效方法。
GLCM用两个位置像素的联合概率密度来定义，它不仅反映亮度的分布特性，也反映具有同样亮度或接近亮度的像素之间的位置分布特性，是有关图像亮度变化的二阶统计特征，可作为分析图像基元和排列结构的基础。
设f(x,y)为一幅二维数字图像，其大小为M×N，灰度级别为Ng,则满足一定空间关系的GLCM为G
式中，F(i，j)是f(x，y)中满足{(x1，y1)，(x2，y2)∈M×N|f(x1，y1)=i，f(x2，y2)=j}的点对数，显然G为Ng×Ng的矩阵，若(x1，y1)与(x2，y2)间距离为d，两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的GLCM:G(i，j，d，θ)。
图像纹理是由GLCM计算出来的特征量来描述的，文献[26]根据GLCM定义了14 种纹理特征，本文选择了角二阶矩(ASM)、熵(ENT)、惯性矩(CON)、逆差矩(HOM)、异质性(DIS)、均值(u) 、标准差(σ)、相关性(COR) 8种常用纹理统计量进行试验。
1.2 纹理特征向量计算GLCM时，在d=1，θ∈{0°，45°，90°，135°}时，得到4个方向上的GLCM，分别计算8个纹理特征值，按顺序排列，得到一个32维的特征向量，用于描述局部图像中心像元位置的纹理信息。由于纹理特征值与窗口大小有关，为了分析窗口大小对变化检测结果的影响，本文对3×3至25×25之间的12个奇数窗口分别进行试验。
1.3 最大化熵法阈值熵[22]原理本质上是“高概率的事物容易出现”这一朴素公理的一个推论。最大熵原理[23]指出，当需要对一个随机事件的概率分布进行预测时，应当满足全部已知的条件，而对未知的情况不要做任何主观假设。匈牙利著名数学家、香农奖得主Csiszar证明，对任何一组不自相矛盾的信息，最大熵模型不仅存在，而且是唯一的，都具有指数函数形式。
通常认为差异影像包含了变化与未变化两类信息，确定把两类信息分开的合理阈值是变化检测的难点之一，本文采用阈值求取与信息熵相结合的策略，认为使差异影像信息熵最大的阈值是要确定的合理阈值，方法如下：
设二维数字图像f(x，y)中灰度值为i的概率为p(i)，则熵H0可以通过式(2)计算
设f(x，y)的灰度分布直方图为h，阈值t把灰度分布直方图分割为h1、h2两部分，按照式(2)分别计算h1、h2的熵H1、H2，则此时整幅图像的熵Ht可通过式(3)计算
设熵向量为H，对∀t∈(0，255)，H[t]的值可通过式(4)计算
设t0为最大化熵法选取的阈值，则Ht0为H的无穷范数，可通过式(5)计算
2 试 验
本文基于Visual Studio 2012和OpenCV 2.4.2 环境，采用C++语言编程实现以上算法并进行试验。
2.1 试验数据试验区域为兖州某矿开采沉陷区，试验数据为2011-12-25、2012-04-02、2012-09-14获取的Terra-SAR影像，进行了两组试验，影像参数见表 1。
| 试验 | 数据来源 | 数据时相 | 波段数量 | 分辨率/m |
| 试验1 | Terra-SAR | 2012-04-02 | 1 | 2.47 |
| Terra-SAR | 2012-09-14 | 1 | 2.47 | |
| 试验2 | Terra-SAR | 2011-12-25 | 1 | 1.96 |
| Terra-SAR | 2012-04-02 | 1 | 1.96 |
图 2为试验1两个时相Terra-SAR影像，大小为449×549像素，影像中主要地物为水体、植被和农田。比较图 2(a)和图 2(b)可以看出两个时相影像的变化情况，为了便于比较变化检测的结果，结合遥感影像和地形图资料，人工描绘了真实的变化图，如图 3所示。
|
| 图 2 Terra-SAR影像(试验1) Fig. 2 Terra-SAR image(experiment 1) |
|
| 图 3 真实变化图(试验1) Fig. 3 The real change(experiment 1) |
图 4为试验2两时相Terra-SAR影像，大小为300×300像素，影像中主要地物为水体、植被和农田。与试验1相比，该地区变化检测难度较高，主要表现在3个方面：①两幅影像间地物类型变化(植被变为水体、水体变为植被、无水的坑塘变为水体、中间一块农田因开采沉陷，由农田变为水体等)更多；②灰度差异影像信噪比较低(农田与水体的灰度差异较小)，受噪声影响严重；③地物变化斑块更多，存在多个小面积变化区域。图 5(a)为灰度差异影像，为了便于比较变化检测的结果，比较图 4(a)和图 4(b)，人工描绘了真实的变化图，如图 5(b)所示。
|
| 图 4 Terra-SAR影像(试验2) Fig. 4 Terra-SAR image(experiment 2) |
|
| 图 5 差异影像与真实变化图(试验2) Fig. 5 The real change(experiment 2) |
由于SAR影像具有明显的斑点噪声，试验时选择Enhanced-Frost方法[24]对原始影像作滤波处理，用SIFT算法[25]对两时相影像进行配准，为了节省内存、提高计算效率并使SAR影像纹理特征明显，对两时相SAR影像作了灰度级压缩处理，SAR影像的初始灰度级为256，压缩之后的灰度级为16。
2.3 构造差异影像对两时相SAR影像中的每一个像元计算纹理特征向量，用两个特征向量间的欧氏距离描述两幅SAR影像间的差异信息，将差异信息归化到0～255之间，得到纹理差异影像。灰度信息采用差值法获得差异影像。
设两幅SAR影像中位置对应的两个像元的纹理特征向量分别为F1，F2，则欧氏距离D 可用向量范数表示为
试验共获得了26幅差异影像，限于篇幅，这里仅列出试验1的4幅差异影像。图 6(a)为利用灰度信息和差值法获得的差异影像，图 6(b)—图 6(d)是窗口大小为5×5、9×9、13×13时获得的纹理差异影像。由图 6可以看出，灰度差异影像中斑点噪声较多，大小为5×5的差异影像纹理特征明显，由不规则豆状块体构成，随着窗口的增大，纹理特征逐渐消失，变化区域逐渐变亮，与不变化区域的灰度差逐渐增大，可区分性增强，细节变化则趋于消失。
|
| 图 6 差异影像(试验1) Fig. 6 The difference image(experiment 1) |
试验中得到的差异影像均采用式(5)计算最大化熵法阈值。通过阈值分割，得到仅含有变化信息和未变化信息的二值影像，其中白色代表变化信息，黑色代表未变化信息。
以试验1中纹理窗口为7×7和9×9的差异影像为例，说明阈值选取的有效性。计算的最大化熵法阈值分别为74和99，以最大化熵法阈值为中心，在双邻域内对称等间距的选择4个阈值，对差异影像进行阈值分割，对得到的变化检测结果用虚警率与漏检率进行精度分析，检测精度与阈值的关系见图 7。
|
| 图 7 纹理差异影像检测精度与阈值关系图 Fig. 7 Detection precision of texture difference image changing with threshold |
从图 7中可见，在两种纹理窗口下，随着阈值的增大，虚警率在减小，漏检率在增大。在最大化熵法阈值附近，两条曲线相交，说明采用最大化熵法阈值可综合考虑虚警率和漏检率，在SAR变化检测中用于阈值分割是合理有效的。
2.5 试验结果本文的试验中，对26幅差异影像以及看作马尔科夫随机场(二阶邻域系统)的2幅灰度差异影像按照最大化熵法阈值分割，得到28幅变化检测结果。试验1和试验2基于马尔科夫随机场和灰度信息的变化检测结果见图 8和图 9，在不同窗口大小下基于纹理特征向量的变化检测结果见图 10和图 11。
|
| 图 8 对比方法的变化检测结果(试验1) Fig. 8 The change detection result of comparative method(experiment 1) |
|
| 图 9 对比方法的变化检测结果(试验2) Fig. 9 The change detection result of comparative method (experiment 2) |
|
| 图 10 基于纹理特征向量的不同窗口大小下变化检测结果(试验1) Fig. 10 The change detection result based on texture feature vector at different window size(experiment 1) |
|
| 图 11 基于纹理特征向量的不同窗口大小下变化检测结果(试验2) Fig. 11 The change detection result based on texture feature vector at different window size (experiment 2) |
纹理窗口大小会影响变化检测的虚警率、漏检率和总体误差。为了建立窗口大小与3个精度指标之间的关系，假设两组试验数据的权值相等，对加权后的数据进行回归分析。以窗口大小为自变量，精度指标为因变量，根据离散点的客观分布趋势，同时为避免复杂模型导致的过拟合现象，选择二阶多项式模型用于回归分析。3个精度指标随窗口大小变化的回归曲线与离散点分布见图 12。
|
| 图 12 精度指标随窗口变化的回归分析 Fig. 12 The regression analysis of window size and precision index |
图 12中虚警率、漏检率和总体误差拟合时的相关系数分别为0.884 6、0.729 8、0.696 0，对应的拟合公式分别为式(7)、式(8)和式(9)
在图 12的回归曲线中，虚警率与漏检率曲线的交点位于窗口9×9和11×11之间，根据回归公式求出窗口大小为9×9和11×11的虚警率和漏检率如表 2所示，其中窗口大小1对应于加权后基于灰度的精度检测指标。
| (%) | |||
| 窗口大小 | 虚警率 | 漏检率 | 总体误差 |
| 1 | 37.51 | 42.23 | 12.57 |
| 9×9 | 34.07 | 30.84 | 12.55 |
| 11×11 | 27.94 | 31.35 | 9.63 |
根据图 12和表 2得到窗口大小与3个精度指标之间的关系如下。
(1) 随着窗口的增大，虚警率和总体误差先减小后增大，在窗口为17×17时均取得最小值。
(2) 漏检率随着窗口的增大而增大，增长率在窗口为11×11之前较小，当窗口大于11×11时，增长率变大。
(3) 基于GLCM纹理特征向量的变化检测结果中，虚警率曲线与漏检率曲线的交点位于窗口9×9和11×11之间。由表 2可知，基于GLCM纹理特征向量的变化检测结果与基于灰度的变化检测结果相比，纹理窗口为9×9时，虚警率减小了3.44%，漏检率减小了11.39%，总体误差减小了0.02%；纹理窗口为11×11时，虚警率减小了9.57%，漏检率减小了10.88%，总体误差减小了2.94%。因此在综合考虑虚警率与漏检率的影响时，最优的奇数纹理窗口为11×11。
基于马尔科夫随机场的方法与本文方法(窗口11×11)变化检测结果的精度比较见表 3。从表 3中可见，基于马尔科夫随机场的变化检测结果与本文方法相比有最小的虚警率，但同时有最大的漏检率，而且总体误差大于本文方法中纹理窗口为11×11时的试验结果。
| % | ||||
| 试验 | 方法 | 总体误差 | 虚警率 | 漏检率 |
| 试验1 | MRF | 8.69 | 8.36 | 42.72 |
| 本文方法(窗口11×11) | 8.08 | 21.63 | 23.44 | |
| 试验2 | MRF | 10.69 | 6.22 | 70.81 |
| 本文方法(窗口11×11) | 9.71 | 32.12 | 35.55 | |
基于GLCM的纹理特征是一种可以减小SAR影像斑点噪声影响的空间邻域信息，用不同方向的多个纹理特征描述差异信息可以提高变化检测的精度和普适性。本文使用基于GLCM的32维纹理特征向量描述两时相SAR影像的差异信息，使用不依赖分布模型假设的最大化熵法自动选取阈值。对3×3至25×25之间的12个奇数窗口，进行SAR影像变化检测试验，通过回归分析得到适合SAR影像变化检测的纹理窗口为11×11。窗口太小，GLCM纹理特征不明显，总体误差较大；窗口太大，变化区域与未变化区域的边缘容易被误检测。试验表明，本文方法优于马尔科夫随机场法，受斑点噪声的影响较小，可以提高高分辨率SAR影像变化检测的精度。接下来的研究内容有两个，一是通过融合其他信息，进一步提高本文方法的稳健性；二是对本文方法中窗口的自动选择问题进行研究，进一步提高本文方法的易用性。
致谢： 特别感谢陈炳乾和杨俊凯在论文写作与修改过程中给予的竭诚帮助。
| [1] | HACHICHA S, CHAABANE F. On the SAR Change Detection Review and Optimal Decision[J]. International Journal of Remote Sensing, 2014, 35(5):1693-1714. |
| [2] | MARINO A, HAJNSEK I. A Change Detector Based on an Optimization with Polarimetric SAR Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(8):4781-4798. |
| [3] | AGHABABAEE H, AMINI J, TZENG Y C. Improving Change Detection Methods of SAR Images Using Fractals[J]. Scientia Iranica, 2013, 20(1):15-22. |
| [4] | AGHABABAEE H, AMINI J, TZENG Y C, et al. Unsupervised Change Detection on SAR Images Using a New Fractal-based Measure[J]. Photogrammetrie-Fernerkundung-Geoinformation, 2013(3):209-220. |
| [5] | XIONG Boli, CHEN J M, KUANG Gangyao. A Change Detection Measure Based on a Likelihood Ratio and Statistical Properties of SAR Intensity Images[J]. Remote Sensing Letters, 2012, 3(3):267-275. |
| [6] | 郝洪美, 张永红, 石海燕, 等. 统计假设检验方法在全极化SAR变化检测中的应用[J]. 遥感学报, 2012, 16(3):520-532. HAO Hongmei, ZHANG Yonghong, SHI Haiyan, et al. Application of Test Statistic Method in Fully Polarimetric SAR Change Detection[J]. Journal of Remote Sensing, 2012, 16(3):520-532. |
| [7] | AGHABABAEE H, TZENG Y C, AMINI J. Swarm Intelligence and Fractals in Dual-pol Synthetic Aperture Radar Image Change Detection[J]. Journal of Applied Remote Sensing, 2012, 6(1):63596-63596. |
| [8] | 黄世奇, 刘代志, 胡明星, 等. 基于小波变换的多时相SAR图像变化检测技术[J]. 测绘学报, 2010, 39(2):180-186. HUANG Shiqi,LIU Daizhi,HU Mingxing,et al.Multi-temporal SAR Image Change Detection Technique Based on Wavelet Transform[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(2):180-186. |
| [9] | 高丛珊, 张红, 王超, 等. 广义Gamma模型及自适应KI阈值分割的SAR图像变化检测[J]. 遥感学报, 2010, 14(4):710-724. GAO Congshan, ZHANG Hong, WANG Chao, et al. SAR Change Detection Based on Generalized Gamma Distribution Divergence and Auto-threshold Segmentation[J]. Journal of Remote Sensing, 2010, 14(4):710-724. |
| [10] | 胡召玲. 广义高斯模型及KI双阈值法的SAR图像非监督变化检测[J]. 测绘学报, 2013, 42(1):116-122. HU Zhaoling. An Unsupervised Change Detection Approach Based on KI Dual Thresholds under the Generalized Gauss Model Assumption in SAR Images[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(1):116-122. |
| [11] | 韩春明,郭华东,王长林.SAR图像斑点噪声抑制的本质[J]. 遥感学报, 2002, 6(6):470-474. HAN Chunming, GUO Huadong, WANG Changlin. The Essence of SAR Image Speckle Suppression[J]. Journal of Remote Sensing, 2002, 6(6):470-474. |
| [12] | MOSER G,SERPICO S B.Unsupervised Change Detection from Multichannel SAR Data by Markovian Data Fusion[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(7):2114-2128. |
| [13] | WANG Guangxue, HUANG Xiaotao, ZHOU Zhimin, et al. A New SAR Image Change Detection Algorithm Based on Texture Feature[C]//Proceedings of the 3rd International Asia-Pacific Conference on Synthetic Aperture Radar. Seoul, Korea:IEEE, 2011. |
| [14] | GONG Maoguo,LI Yu,JIAO Licheng,et al.SAR Change Detection Based on Intensity and Texture Changes[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 93:123-135. |
| [15] | 陈志鹏, 邓鹏, 种劲松, 等. 纹理特征在SAR图像变化检测中的应用[J]. 遥感技术与应用, 2002, 17(3):162-166. CHEN Zhipeng, DENG Peng, CHONG Jinsong, et al. Application of Textural Features to Change Detection in SAR Image[J]. Remote Sensing Technology and Application, 2002, 17(3):162-166. |
| [16] | 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. |
| [17] | WONG A K C, SAHOO P K. A Gray-level Threshold Selection Method Based on Maximum Entropy Principle[J]. IEEE Transactions on Systems, Man and Cybernetics, 1989, 19(4):866-871. |
| [18] | KAPUR J N, SAHOO P K, WONG A K C. A New Method for Gray-level Picture Thresholding Using the Entropy of the Histogram[J]. Computer Vision, Graphics and Image Processing, 1985, 29(3):273-285. |
| [19] | 刘丽, 匡纲要. 图像纹理特征提取方法综述[J]. 中国图象图形学报, 2009, 14(4):622-635. LIU Li, KUANG Gangyao. Overview of Image Textural Feature Extraction Methods[J]. Journal of Image and Graphics, 2009, 14(4):622-635. |
| [20] | TUCERYAN M, JAIN A K. Texture Analysis[M]//CHEN C H, PAU L F, WANG P S P. The Handbook of Pattern Recognition and Computer Vision.[S.l.]:World Scientific Publishing Co., 1998:207-248. |
| [21] | HARALICK R M, SHANMUGAM K, DINSTEIN I H. Textural Features for Image Classification[J]. IEEE Transactions on Systems, Man and Cybernetics, 1973, SMC-3(6):610-621. |
| [22] | RÉNYI A. On Measures of Entropy and Information[C]//Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability. Berkeley:University of California Press, 1961. |
| [23] | KAPUR J N, KESAVAN H K. Entropy Optimization Principles with Applications[M]. Boston:Academic Press, 1992. |
| [24] | 贾伟洁, 张继贤, 杨景辉, 等. 基于滤波的SAR影像去噪方法分析[J]. 微计算机信息, 2012, 28(6):24-26. JIA Weijie, ZHANG Jixian, YANG Jinghui, et al. The Analysis of Noise Reduction for SAR Image Based on Filters[J]. Microcomputer Information, 2012, 28(6):24-26. |
| [25] | LOWE D G. Distinctive Image Features From Scale-invariant Keypoints[J]. International Journal of Computer Vision, 2004, 60(2):91-110. |