1 引 言
合成孔径雷达(synthetic aperture radar,SAR)因具有全天时、全天候和对某些地物的穿透性能等对地观测能力，成为当前遥感发展的前沿技术之一。与光学遥感相比，雷达遥感具有主动式、多极化、可变观测角度、宽幅成像等特点，其对探测目标的几何特征和物理特征也非常敏感[1]。随着ALOS PALSAR、TerraSAR-X及Radarsat 2等星载SAR系统的运行，多波段SAR图像被用于许多地理要素的提取、反演和变化监测，如全球森林监测、农作物长势监测、灾害监测等[2, 3, 4, 5, 6, 7, 8]，这些应用都需要自动而可靠的多时相SAR图像非监督变化检测技术和信息提取方法。
目前国内外针对光学遥感图像，根据地物的光谱特征、纹理特征以及上下文关系等，对传统的非监督变化检测方法——变化向量分析法进行改进，提高了土地覆盖非监督变化检测的精度[9, 10, 11, 12, 13, 14, 15]。SAR图像的成像机理和地物的影像特征不同于光学图像，由于受SAR图像斑点噪声和分类精度等因素的影响，目前SAR图像在地理要素变化监测方面的应用明显不足。多时相SAR图像的非监督变化检测方法的标准探测器是基于局部平均值的比值法[16]，该方法探测地物的突变效果较理想，但几乎探测不到地物的弱变化信息。文献[17]基于广义高斯模型，文献[16]基于局部统计量相似性测度实现SAR图像土地覆盖的非监督变化检测；文献[18]基于特征和像素值利用多时相SAR图像对城市地区进行了非监督变化检测，这些变化检测技术均是针对区分变化区域和非变化区域而提出的，因而仅能提取变化与非变化这两种信息。然而在复杂的地理环境下，地物的后向散射特征易受多种因素的影响：地物的生物、物理参数、介电常数、环境因素(如天气状况)、雷达参数(雷达波入射角、波长、极化方式)等[19]，因而即使是未变的地物，在不同时相SAR图像上的灰度值及其统计特征也存在着不同程度的差异，这些差异增加了变化类型的复杂性和变化检测技术的难度。由于地物的变化非常复杂，为了从SAR差异图像上提取出更丰富的信息，对于变化区域，需要进一步识别出后向散射增强区域和后向散射减弱区域，并且这两种区域中像素灰度值的统计特征也不相同，因此需要选取双阈值对差异图像进行分割。
用于SAR差异图像分割的阈值选取方法主要有最小误差法(KI)和期望最大化法(EM)，这两种方法都需要假定差异图像上变化与非变化类型灰度值的统计分布，统计分布模型不同，阈值选取的方法也就不同。目前已经发展了基于高斯分布模型、广义高斯分布模型[20]和广义Gamma模型的KI准则单阈值选取方法[21]，仅能识别变化与非变化这两种类型。文献[20]提出了基于高斯模型的EM双阈值选择方法，假定不同变化类型的灰度值服从高斯分布，采用EM算法产生双阈值，在变化的区域上进一步区分出变化区域增强类和变化区域减弱类。尽管广义高斯模型更适于描述差异图像上各变化类型的概率分布[17]，但很难构造出基于广义高斯模型的EM双阈值解算方程。本文基于广义高斯分布模型，将KI准则进行推广，提出仅利用差异图像的灰度直方图自动选取最优双阈值的方法，以实现未发生变化类、后向散射减弱类和后向散射增强类这3种类型的非监督变化检测信息提取。
2 基于广义高斯分布模型和KI准则的双阈值选择 2.1 SAR差异图像的提取假设X1和X2分别是t1、t2时刻所获得的，经过辐射校正、配准及滤波后的同一区域的SAR图像，图像的大小均为I像素×J像素，X1(i,j)和X2(i,j)分别是相应图像上第i行、第j列像素的灰度值，其中1≤i≤I，1≤j≤J。由于SAR图像乘性斑点噪声的特点，常用对数比值法来构建SAR差异图像[17]，其计算方法如下
h(Xl)(Xl=0，1，…，L-1)是差异图像的灰度直方图，表示灰度值Xl的概率。假设差异图像上有3种类型，分别是：未发生变化类ωu、变化类别ωc1(对应后向散射减弱类)和变化类别ωc2(对应后向散射增强类)。根据贝叶斯决策理论，任一像素Xl的概率密度计算方法如下
式中，P(ωu)、P(ωc1)和P(ωc2)分别表示3种类型的先验概率；p(Xl/ωu)、p(Xl/ωc1)和p(Xl/ωc2)表示条件概率密度函数。 2.2 KI准则的双阈值选取
KI准则的单阈值选取方法是在最优化预先定义的准则函数J(T)的基础上选择合理的阈值T∈{0,1,…,L－1}，准则函数平均了直方图上的价值函数c(Xl,T)，使得图像的总体分类效果最好。价值函数c(Xl,T)通过对灰度水平和阈值的比较，测度了分类像素的价值。本文将KI准则单阈值的选择方法推广到双阈值的选择，定义准则函数J(T1,T2)如下
式中，T1、T2∈{0,1,…，L－1}是所选择的双阈值。代价函数的计算如下 式中，P(ωi/Xl,T1,T2)(i=c1,u,c2)分别是3种类型在给定的灰度值Xl和双阈值T1、T2下的后验概率。选取最佳阈值要使得图像的分类误差最小，即使J(T1,T2)取得最小值
根据贝叶斯准则，后验概率可由先验概率和条件类的概率密度函数计算出来，其计算方法如下
其中条件类的概率密度函数则取决于所构建的类型统计模型。 2.3 基于广义高斯分布模型的KI双阈值选取
研究表明，SAR信号由于具有显著的尖峰和拖尾现象，采用广义高斯分布模型描述地物的概率密度分布函数要优于高斯模型[17]。与高斯模型相比，构建广义高斯模型最多需一个形状参数。本文在基于广义高斯分布的单阈值准则函数的基础之上，构造出双阈值准则函数。基于广义高斯分布模型，3种类型的概率密度函数为
式中
mi、σi2和βi分别表示分布函数的均值、方差和形状参数；Γ(·)指的是Gamma函数，以未变化类型ωu为例，βu的估计值计算如下[17]：
(1) 计算广义高斯比率函数r(βu)，以产生查找表
(2) 计算|Xl－mu|的均值
(3) 计算比率
(4) 根据查找表，计算
。
将基于广义高斯模型的单阈值准则函数进行推广，定义双阈值准则函数如下
式中，H(Ω,T1,T2)表示类别集合Ω∈{c1,u,c2}的熵，其计算公式如下
式(7)～(11)中各类型的先验概率Pi(T1,T2)、均值mi(T1,T2)和方差σi2(T1,T2)(i=c1,u,c2)可根据差异图像的灰度直方图h(Xl)计算出，具体如下
假设各类型的条件概率密度服从广义高斯分布，所选取的最优分割阈值T1、T2，要使式(10)的取值最小。
3 试验与分析 3.1 试验数据描述本文以徐州市某高校新校区为试验区，选取2个时相的SAR图像进行非监督变化检测试验，图像的获取时间分别是1999年4月18日和2009年7月30日，分别来源于加拿大Radardat-1和Radardat-2卫星，均采为C波段的波束、HH极化方式。其中1999年的波束入射角为37°～41°，图像的距离向分辨率为8.3 m，方位向分辨率为8.4 m；2009年的波束入射角为20°～41°，图像的距离向分辨率为16.5～8.0 m，方位向分辨率为8.0 m。首先对SAR图像进行几何校正、配准、重采样、滤波等预处理，选取201像素×201像素的图像为试验数据，图像像素的实际尺寸是8.3 m×8.3 m(图 1)。
|
| 图 1 试验区SAR图像 Fig. 1 SAR image in the study area |
1999年试验区内的地物主要包括山林、麦田和农村居民点，如图 1(a)中的A、B、C处所示地物，由于新校区的建设，2009年该区发生了复杂的变化，图 1(b)中D、E、F、G处的地物分别是草地、操场、建筑物、人工湖。由于两个年份图像获取的时间跨度较大、季节不同，加之传感器的波束入射角、分辨率、斑点噪声等因素的影响，同一地理位置的像素灰度值呈现出复杂的变化。特别是新建的教学楼、体育馆等建筑物，有与雷达波束相垂直的平面，其后向散射回波会产生角反射效应，在雷达图像上表现为一定形状的亮线，而其他平面由于表面光滑，则表现为较暗的区域，因此建筑物区域在图像上呈现出明暗相间的影像特征，其亮线具有一定的方向性，目视解译很容易识别，但自动变化检测的难度较大。
3.3 试验结果分析计算图 1(b)与图 1(a)对应像素的比值，获得比值差异图像，对差异图像的像素值取对数后，再进行线性增强，获得增强后的对数比值差异图像Xl(图 2)。计算出差异图像中各灰度值像素数占总像素数的比例，从而得到差异图像的灰度直方图h(Xl)。
|
| 图 2 对数比值差异图像 Fig. 2 Difference image of log-ratio |
由第2节，基于广义高斯分布模型KI准则对图像Xl进行双阈值选取，获得最佳分割阈值T1=71、T2=93，根据这两个阈值对图像Xl的像素逐个进行判别，获得变化类型图像(图 3)，判别准则如下
式中，ωu是指未发生变化类(图 3中的灰色调区域)；ωc1是发生变化的后向散射减弱类(图 3中的黑色调区域)；ωc2为发生变化的后向散射增强类(图 3中的白色调区域)。
|
| 图 3 变化类别图像 Fig. 3 Change type image |
图 3中变化类型ωc1主要包括农村居民点变为人工湖(图 1中的C变为F)、麦田变为操场(图 1中的B变为E)、麦田、草地变为道路等变化类型；变化类型ωc2主要是具有中等灰度值的麦田、草地等地物变为具有较大灰度值的建筑物(图 1中的B变为G)，从图 3中可以看出建筑物的位置和布局；未发生变化类ωu主要是山林。
根据所选取的阈值，计算出这3种类型的广义高斯分布模型参数(表 1)，绘制出各类型的条件概率分布密度曲线图(图 4)，类别ωc1与类别ωu概率密度曲线的交点所对应的灰度值为73.3，类别ωu与类别ωc2的交点所对应的灰度值为92.2。根据类别ωc1与ωu的先验概率P、阈值71在这2类中的条件概率p(71/ωc1)和p(71/ωu)，可计算出类别在阈值71处的后验概率p(ωc1/71)和p(ωu/71)分别为0.002 1和0.002 9，同样的，计算出类别ωu与ωc2在阈值93处的后验概率p(ωu/93)和p(ωc2/93)分别为0.004 5和0.005 4。根据EM算法，在阈值处两相邻类别的后验概率应相等，但由于像素的个数有限以及灰度值是正整数的原因，两类别的后验概率并不相等，但据此可以验证出与其他灰度值相比，这两个阈值处两类别的后验概率相差最小，从而说明了本文提出的SAR图像非监督变化检测方法是有效可行的。
| 参数 | ωc1 | ωu | ωc2 |
| a | 0.187 3 | 0.279 9 | 0.287 9 |
| b | 0.561 8 | 0.909 3 | 0.253 7 |
| β | 0.601 4 | 0.567 5 | 0.616 8 |
| m | 60.394 1 | 84.301 5 | 116.677 2 |
| σ | 9.100 8 | 6.992 5 | 18.434 2 |
| P | 0.207 7 | 0.408 8 | 0.383 5 |
|
| 图 4 各类型的概率分布密度图 Fig. 4 Probability distribution density of types |
通过分析比较1999年、2009年的土地利用现状图、高分辨率光学遥感图像，并结合遥感图像目视判读结果以及部分实地调查资料，本文提出的变化检测方法的总体误差为482像素。为进行对比分析，采用基于高斯模型的EM双阈值选取方法对图 2进行阈值选取，阈值为T1=87、T2=101，其总体误差为785像素，因此基于广义高斯模型的KI双阈值自动分割SAR图像非监督变化检测精度有显著的提高。
为了进一步验证本文提出的SAR图像非监督变化检测方法的可行性和有效性，取1999年、2007年两个时相的徐州市新城区Radardat-1卫星SAR图像进行试验，通过查阅当时的气象资料，这两个时相的图像是在相似的大气条件和地表湿度条件下获取的，尺寸为513像素×513像素，如图 5所示。根据本文提出的方法，计算出最佳分割阈值为T1=65、T2=105，变化检测结果如图 6所示，图中的灰色调区域是指未发生变化类，黑色调区域是发生变化的后向散射减弱类，白色调区域是发生变化的后向散射增强类。图 5(a)中选取A、B、C、D 4个样区，共计11 296像素，对照地面资料，其中未发生变化的像素个数为2397，后向散射减弱类像素个数为4794，后向散射增强类像素个数为4105。1999年样区的地物分别是麦田、大龙口水库的水体、居民区和鱼塘，在图 5(b)中，淮海食品城的建设，使得A变成了建筑区；由于徐州市新城区的建设，在大龙口水库的基础上开挖修建了大龙湖，湖的周边种植了林木、草地等进行绿化，B区域有的水体部分未发生变化，图 6中呈现灰色调，依然可以看出原来的水库形状，有的水体部分则发生了变化，其中白色调的区域则是由水体变为绿地的部分；由于新城区的建设，许多居民区被拆迁，C区域由原来的居民区变成了绿地，色调由亮变暗；位于大龙口水库边的区域D原来是鱼塘，随着大龙湖的开挖及其周边环境的变化，也随之发生变化，在图 6中能够被检测出来，呈现出白色调的小格网状。除此之外，在大龙湖周围新修的几条道路也被检测出来，呈现出黑色调的条带状。经检验，样区的变化检测率达到92.6%，虚警率仅为2.3%。
|
| 图 5 Radardat-1卫星SAR图像 Fig. 5 SAR image from Radardat-1 satellite |
|
| 图 6 变化类别图像 Fig. 6 Change type image |
在复杂的地理环境下，多时相SAR图像非监督变化检测技术不仅要检测出变化的区域，还需要对变化区域进行细分，进而识别出后向散射增强类和后向散射减弱类。与传统的高斯模型相比，构建广义高斯模型只需增加1个形状参数，它更适于描述多时相SAR对数比值差异图像上各变化类型的概率分布。本文采用广义高斯分布模型描述未变化类、后向散射减弱类和后向散射增强类这3种类型在SAR差异图像上的概率密度分布，各模型参数均根据差异图像的灰度直方图自动计算，并基于KI准则，构建了双阈值准则函数，提出了基于广义高斯分布模型的KI准则最优双阈值自动选取方法。试验结果表明该方法可行、有效，与基于高斯模型的EM双阈值选取方法相比，变化检测精度显著提高。如何充分利用空间上下文信息，提高SAR图像非监督变化检测技术的精度和速度是下一步研究的重点。
| [1] | XU Yibo, LAI Xijun, ZHOU Chunguo. Study on the Remote Sensing Monitoring of Wetland Vegetation in East Dongting Lake Using ENVISAT ASAR Data [J]. Resources and Environment in the Yangtze Basin, 2010,19(4):452-459.(徐怡波,赖锡军,周春国. 基于ENVISAT ASAR数据的东洞庭湖湿地植被遥感监测研究[J]. 长江流域资源与环境,2010,19(4):452-459.) |
| [2] | LNNQVIST A, RAUSTE Y, MOLINIER M, et al. Polarimetric SAR Data in Land Cover Mapping in Boreal Zone [J]. IEEE Transactions on Geoscience and Remote Sensing, 2010,48(10): 3652-3662. |
| [3] | ZHENG Wei, LIU Chuang, WANG Zhengxing. ENVISAT ASAR Data Based Extraction of Flood Spatial Distribution Information [J]. Journal of Natural Disasters, 2009,18(4):120-124.(郑伟,刘闯,王正兴. 基于ENVISAT-ASAR数据的洪涝水体空间分布信息提取[J]. 自然灾害学报,2009,18(4):120-124.) |
| [4] | BOUVET A, TOAN T L, LAMDAO N. Monitoring of the Rice Cropping System in the Mekong Delta Using ENVISAT/ASAR Dual Polarization Data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2009,47(2):517- 526. |
| [5] | AMINI J, SUMANTYO J. Employing a Method on SAR and Optical Images for Forest Biomass Estimation [J]. IEEE Transactions on Geoscience and Remote Sensing, 2009,47(12):4020-4026. |
| [6] | CHAKRABORTY M, MANJUNATH K R, PANIGRAHY S, et al. Rice Crop Parameter Retrieval Using Multi-temporal, Multi-incidence Angle Radarsat SAR Data [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2005,59(5):310-322. |
| [7] | SANTORO M, FRANSSON J, ERIKSSON L B. Signatures of ALOS PALSAR L-band Backscatter in Swedish Forest[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009,47(12):4001-4019. |
| [8] | TANASE M A, SANTORO M, RIVA J, et al. Sensitivity of X-, C-, and L-band SAR Backscatter to Burn Severity in Mediterranean Pine Forests [J]. IEEE Transactions on Geoscience and Remote Sensing, 2010,48(10):3663-3675. |
| [9] | ZHANG Jinshui, PAN Yaozhong, HAN Lijian, et al. Land Use/Cover Change Detection with Multi-source Data [J]. Journal of Remote Sensing, 2007, 11(4):500-510. (张锦水,潘耀忠,韩立建, 等. 光谱与纹理信息复合的土地利用/覆盖变化动态监测研究[J].遥感学报,2007, 11(4):500-510.) |
| [10] | GHOSH S, BRUZZONE L, PATRA S, et al. A Context-sensitive Technique for Unsupervised Change Detection Based on Hopfield-type Neural Networks [J]. IEEE Transactions on Geoscience and Remote Sensing, 2007,45(3):778-789. |
| [11] | BOUCHER A, KAREN C S, JOURNEL A G. A Novel Method for Mapping Land Cover Changes: Incorporating Time and Space with Geostatistics [J]. IEEE Transactions on Geoscience and Remote Sensing, 2006,44(11):3427-3435. |
| [12] | NEMMOUR H, CHIBANI Y. Multiple Support Vector Machines for Land Cover Change Detection: an Application for Mapping Urban Extensions [J]. ISPRS Journal of Photogrammetry & Remote Sensing, 2006,61(5):125-133. |
| [13] | BOVOLO F, BRUZZONE L. A Theoretical Framework for Unsupervised Change Detection Based on Change Vector Analysis in the Polar Domain [J]. IEEE Transactions on Geoscience and Remote Sensing, 2007,45(1):218-236. |
| [14] | CHEN J, GONG P, HE C R, et al. Land-use/Land-cover Change Detection Using Improved Change-vector Analysis [J]. Photogrammetry Engineering and Remote Sensing, 2003, 69(4): 369-379. |
| [15] | WANG Min, ZHANG Xingyue. Change Detection Using High Spatial Resolution Remotely Sensed Imagery by Combining Evidence Theory and Structural Similarity [J]. Journal of Remote Sensing, 2010, 14(3):558-570. (汪闽,张星月.多特征证据融合的遥感图像变化检测[J].遥感学报,2010,14(3):558-570.) |
| [16] | INGLADA J, MERCIER G. A New Statistical Similarity Measure for Change Detection in Multitemporal SAR Images and Its Extension to Multiscale Change Analysis [J]. IEEE Transactions on Geoscience and Remote Sensing, 2007,45(5):1432-1445. |
| [17] | BAZI Y, BRUZZONE L. 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. |
| [18] | GAMBA P, DELLACQUA F, TRIANNI G. Rapid Damage Detection in the Bam Area Using Multitemporal SAR and Exploiting Ancillary Data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2007,45(6):1582-1589. |
| [19] | SHAO Yun, LIAO Jingjuan, FAN Xiangtao, et al. Analysis on Rice Backscatter Signatures in Time Domain: Comparison between Radarsat-SAR Observation and Simulated Model Results [J]. Journal of Remote Sensing, 2002, 6(6):440-449. (邵芸,廖静娟,范湘涛,等.水稻时域后向散射特性分析:雷达卫星观测与模型模拟结果对比[J].遥感学报,2002,6(6):440-449.) |
| [20] | 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. (黄世奇, 刘代志, 胡明星,等.基于小波变换的多时相SAR图像变化检测技术[J]. 测绘学报,2010,39(2):180-186.) |
| [21] | 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):718-724. (高丛珊,张红,王超,等.广义Gamma模型及自适应KI阈值分割的SAR 图像变化检测[J].遥感学报,2010,14(4):718-724.) |