2. 山东正元地理信息工程有限责任公司 潍坊分公司,山东 潍坊261021;
3. 辽宁工程技术大学 测绘与地理科学学院,辽宁 阜新 123000;
4. 山东科技大学 测绘科学与工程学院,山东 青岛 266510
2. Shandong Zhengyuan Geographic Information Engineering Company Limited(Weifang Branch),Weifang 261021,China;
3. School of Mapping and Geographic Science,Liaoning Technical University,Liaoning 123000,China;
4. College of Geomatics,Shandong University of Science and Technology,Qingdao 266510,China
由于成像不受天气及时间的影响,基于合成孔径雷达影像的变化检测技术在灾害管理(如地震、洪涝等自然灾害的灾情评估与监测)、资源调查(如农作物长势监测)等方面具有独特的优势。但是,SAR固有的斑点噪声增加了变化检测处理的复杂性和不确定性[1]。基于光学影像的变化检测在很多应用领域实现了实用化[2],相比之下,对SAR影像变化检测的研究还远未成熟。
早期的SAR变化检测主要是利用单极化SAR影像。由于斑点噪声为乘性噪声,为了减弱斑点噪声的影响,变化检测主要是利用前后时相SAR幅度或强度影像的比值[1, 3, 4, 5, 6],或者是前后时相SAR影像比值对数[7, 8],生成一个包含变化信息的差异影像,然后对差异影像进行阈值处理提取变化信息。利用小波变化对噪声不敏感的优势,国内学者提出了将比值对数图像进行小波变换并通过EM算法确定变化阈值的方法[9],以及将小波变换和区域级决策融合结合在一起的SAR变化检测技术[10, 11]。全极化SAR影像包含更为丰富的地物信息,因此在变化检测方面比单极化更有优势。国内外已经发表的全极化变化监测方法主要有3类。第1类是将传统变化检测技术移植到极化SAR影像中的,如文献[12]将主成份分析的思想应用于多时相极化SAR影像变化检测。这类方法的不足之处是对极化信息使用不充分。第2类是利用一些极化特征进行变化信息提取,例如文献[13]通过对目标协方差矩阵的特征值和特征向量提取广义标准差、相对比率和椭率角等参数的分析提取变化区域,文献[14]基于复相干因子发展了一个四维向量,通过对此向量的分析提取变化信息。第3类是基于极化SAR影像统计分布的变化检测[15, 16, 17, 18],其中文献[17, 18]提出的基于极化协方差矩阵的似然比检验的变化检测方法是其中的代表。当极化通道的后向散射系数或者通道间的相位差/相关性这两类指标中的任何一个发生变化时,该方法都有较好的效果。但是,这类方法的不足是需要对该极化协方差矩阵做统计分布假设,高分辨率SAR影像的统计分布比较复杂。同时由于似然比统计量的概率密度函数无法显式表达,在利用给定的显著性水平确定拒绝或接受原假设(前后时相协方差矩阵相等,即无变化)的临界值时,使用了大样本条件下的卡方分布近似替代,因此,该方法理论上只在使用了很高视数(大于30)进行多视处理的数据集上或者变化图斑很大时有效,对于一般多视的极化SAR数据或者小图斑的变化检测可靠性无法保证。
针对上述情形,本文从用于极化SAR影像分类的Wishart距离出发,提出一种新的适用于变化检测的极化距离测度,并发展了相应的极化变化检测方法,最后利用北京地区的多时相全极化RADARSAT-2影像进行试验分析。该方法无需考虑极化SAR影像的统计分布,计算简便,而且在像元级进行处理,对小的变化图斑也能检测。
1 极化距离测度全极化SAR影像记录了地面每个分辨单元在4种基本极化状态,即hh、hv、vh、vv(hh表示水平发射/水平接收状态,其他类推)的散射回波的幅度和相位,形成如下散射矩阵S
通常假定介质满足互易性,此时Svh=Shv。
自然界中的地物绝大多数为分布式目标,对于这类地物,通常用极化相干矩阵T3或者极化协方差矩阵C3来描述其散射信息[19]。C3矩阵和T3矩阵均为复对称矩阵,二者之间酉相似。因此,在下面的叙述中只以极化相干矩阵为例,结论同样适用于极化协方差矩阵。
C3及T3矩阵的迹,被称为span,表达了地物目标总的散射功率
由于span为4个通道的强度之和,因此比任一通道的强度影像都有更低的斑点噪声[19]。
多视处理的极化SAR影像的相干矩阵T3服从复Wishart分布[20],此时,根据最大似然分类原理,可以推导出如下的样本相干矩阵与类中心相干矩阵之间的Wishart距离
式中,Vm为第m类中心的相干矩阵,为该类中所有像元相干矩阵的平均值;|·|表示取矩阵的行列式,上标-1表示求矩阵的逆。任一像元可根据其与所有类别中心相干矩阵的Wishart距离最小值判断其归属于哪一类[20, 21]。
在此基础上,文献[22]提出了两个类别中心相干矩阵之间的距离测度
式(4)结果越小,表明两个相邻的类别极化散射特性越相近,则可以合并为同一类别,最终可实现极化SAR影像的非监督聚类。
2 基于极化距离测度的变化检测变化检测与分类之间具有一定的联系。极化分类是将极化特性相似的像元归为同一类别,而变化检测则是要从前后时相影像上找出极化特性相差较大的像元。由于极化相干矩阵包含了地物极化散射的完整信息,因此,可以通过前后时相的极化相干矩阵之间的距离大小来判断像元是否发生变化。
但是,式(4)定义的极化距离测度用于变化检测还有明显不足。主要是距离定义中的第1部分d1(Vi,Vj)(即12ln|Vi|+ln|Vj|)与像元是否变化无关,因为不同类别的极化协方差矩阵可能有很接近的行列式值,同时这部分的值与SAR影像是否经过辐射定标有关。但是,距离定义中的第2部分d2(Vi,Vj)(即12trVi-1Vj+Vj-1Vi)对变化敏感。当像元在前后时相完全一样时,该值为常数3,若有变化发生,则必然导致前后时相极化相干矩阵的不同,此时该值将明显大于3。
为了进一步说明这种情况,从本文所用的试验影像中(已经过辐射定标处理)选择了一些样本,这些样本包括在前后时相上地类没有发生变化的像元8对及地类改变的像元8对,针对前后时相的相干矩阵,分别计算了d1(Vi,Vj)及d2(Vi,Vj)的值,列于表 1。
样本点序号 | 样本点属性 | d1(Vi,Vj) | d2(Vi,Vj) |
1 | 草地→建筑 | -6.21 | 10.05 |
2 | 草地→建筑 | -6.15 | 14.35 |
3 | 裸地→建筑 | -9.64 | 35.60 |
4 | 裸地→建筑 | -7.61 | 27.13 |
5 | 裸地→草地 | -9.71 | 26.79 |
6 | 裸地→草地 | -9.37 | 15.17 |
7 | 水体→草地 | -13.85 | 19.71 |
8 | 水体→草地 | -13.25 | 13.34 |
9 | 草地→草地 | -7.68 | 5.45 |
10 | 草地→草地 | -7.20 | 4.94 |
11 | 裸地→裸地 | -9.44 | 4.30 |
12 | 裸地→裸地 | -9.39 | 3.55 |
13 | 水体→水体 | -15.78 | 6.10 |
14 | 水体→水体 | -14.62 | 6.10 |
15 | 建筑→建筑 | -4.23 | 5.54 |
16 | 建筑→建筑 | -2.57 | 5.59 |
从表 1看出,变化像元与非变化像元的d1(Vi,Vj)值相混合,表明其对变化完全不敏感。但是所有变化像元的d2(Vi,Vj)均显著高于未变化像元。
因此,本文提出如下的适用于变化检测的极化距离测度
式中,T3,1、T3,2分别表示前后时相像元的极化相干矩阵;减3是为了保证同质像元之间的距离为零,以满足距离测度的非负性。对此极化距离测度的进一步分析如下。相干矩阵为埃尔米特半正定矩阵,可表示为
式中,λi为特征值;ui为相应的特征向量,也是酉矩阵U的列向量;上标H表示共轭转置。经推导,上述极化距离测度可表示为
式中,(λi,ui)、(ρi,vi)分别为T3,1、T3,2的特征值及特征向量。从式(7)可以看到,极化距离测度是对两个相干矩阵所包含的全部散射机理的相关性及其幅度比值的综合比较,因此能敏感地反映相干矩阵之间的差别,是一种具有普适性的极化距离测度。
将上述极化距离作为变化检测的差异图像。对差异图像进行阈值分割从而提取变化图斑是遥感变化检测中另一关键问题,国内外研究者提出了很多方法,其中马尔科夫随机场模型方法较有优势。该方法将图像的变化提取问题转化为图像标记问题,通过利用马尔科夫随机场与吉布斯随机场之间的等价关系,将差异图像中的灰度信息与空间信息融入寻求全局最优的变化提取结果的过程中[23],因此能可靠、准确地提取变化信息。该方法引入了邻域信息,对SAR斑点噪声引起的伪变化有很好的抑制作用[24]。因此本文利用马尔科夫随机场模型法完成从极化距离图像中提取最后的变化图斑。相应的,基于极化距离测度的变化检测方法的流程如图 1所示。
3 试验及分析 3.1 数据预处理试验数据为北京西部石景山地区2009年9月与2010年7月获取的两景RADARSAT-2全极化单视复影像(SLC),均为降轨图像,波束模式为Fine Quad-Polarization。具体参数如表 2所示。
成像时间 | 像元大小(距离向×方位向) | 像元个数(距离向×方位向) | 中心入射角/(°) | 中心经纬度/(°) |
2009-09-23 | 4.733 m×4.765 m | 3248×6134 | 35.314 | (116.151,39.943) |
2010-07-15 | 4.733 m×4.742 m | 3644×6207 | 39.998 | (116.182,39.938) |
为了保证两幅影像的辐射一致性,首先对两幅SLC影像分别进行了如下辐射定标处理
式中,Iσ、Qσ分别为定标后像元的实部和虚部;DNI、DNQ为原始单视复影像像元的实部和虚部;Aj是与斜距有关的增益系数,由附带的查找表文件(lutSigma.xml)提供,每一列对应一个增益系数[25]。
辐射定标后对两幅SAR影像进行精确配准。由于两幅SAR影像入射角有近5°的差异,在地形起伏较大区域无法精确配准,因此将原始影像中相对平坦的子区(大小为1340像元×2140像元)作为试验区。为了消除入射角差异造成像元地面分辨率不一致对精确配准的影响,首先分别将其转换成具有方形像元的地距图像,且方位分辨率与原图一致。然后分别提取两幅地距影像的span图像,通过span图像得到两幅影像间的配准多项式。以2009年的span图像为参考图像,沿距离向和方位向均匀选择了20×25个窗口(窗口大小为64像元×64像元,称为目标窗口),对每一个目标窗口影像,在待配准影像(2010年span影像)的搜索窗口内利用最大互相关系数准则寻找同名点[26]。为了确定亚像元级的同名点位置,将目标窗口影像和搜索窗口影像均进行了2×2倍过采样。对所有匹配到的同名点对,利用二次多项式拟合其对应的几何位置关系,对拟合偏差大于1个像元的同名点对进行剔除,迭代求得配准多项式。最后,根据最终得到的配准多项式将2010年地距影像的每一个极化通道的实部和虚部均变换至参考影像的坐标系统内,利用双三次采样方式赋值,完成了两幅SAR影像每一个极化通道间的精确配准。
3.2 变化检测辐射定标和精确配准后的两幅地距极化散射SAR影像依然用极化散射矩阵表达,记为S2009、S2010,然后进行2×2的多视生成各自的极化相干矩阵i>T3,2009、T3,2010。多视处理不仅抑制了斑点噪声,而且使得配准精度再提高一倍。为进一步抑制斑点噪声,对相干矩阵施行了窗口大小为3×3的增强Lee滤波。至此完成了全部的预处理工作。预处理后的两幅极化SAR影像显示如图 2,大小为1077行×1140列。图中假彩色合成的RGB通道分别为12|Shh-Shv|2、2|Shv|2、12|Shh+Shv|2。
3.3 变化检测分析在预处理后影像的基础上,根据式(4)生成了2009年9月23日和2010年7月15日两幅影像间的极化距离图像dc(T3)。
作为对比,分别利用极化总功率和各个极化通道的回波强度图像生成了如下的差异影像
上述基于后向散射强度对数比值的差异性变量是SAR变化检测的常用度量。
为了定量比较本文提出的极化距离测度和上述基于后向散射强度的差异性度量的变化检测能力,从预处理后的影像中选取了4个感兴趣区域为样本。它们分别代表了研究区内4种典型的地表覆盖变化:水体变为草地、建筑物变为裸地、草地变为建筑物、草地变为裸地。样本的选择及地表覆盖类型的判断均参考了从Google Earth下载的高分辨率光学影像。在每个感兴趣区域中分别手动勾绘了一个地表覆盖发生变化的子区和一个与其相邻的未变化子区。然后统计了这些子区分别在dc(T3)、dspan、dshh、dshv、dsvv等5个差异图像上的均值和标准差,结果分别列于图 3~图 6。在图 3-图 6的图(a)及图(b)中,红色边框标记的为变化子区,蓝色或黄色边框标记的为非变化子区;可分离性统计图(c)表示在上述5个图像中的变化区域和非变化区域的均值及标准差,均值以带颜色的点标示,横坐标1、3、5、7、9分别表示变化区域的dc(T3)、dspan、dshh、dshv、dsvv;2、4、6、8、10分别表示非变化区域的dc(T3)、dspan、dshh、dshv、dsvv。
由图 3-图 6可以看出,本文提出的极化距离测度对上述全部4种类型的地表覆盖变化都具有较好的分离能力。dspan对水体变为草地、建筑物变为裸地、草地变为建筑物等3种变化类型具有较好的分离能力,但无法区分草地到裸地的变化。其原因是水体和建筑物有明显的后向散射特点,水体在每一个极化通道上均只有很弱的后向散射,建筑物由于其几何特点(墙面与地面构成二面角散射),在每一个极化通道上均呈现很强的后向散射,因此地物在前或后时相为水体或建筑物之一时,在dspan上都会呈现较大的值;裸地与草地的后向散射相近,差别主要表现在裸地的后向散射中以奇次散射占优,而草地的后向散射中漫散射占优,但在极化总功率上非常接近。dshv对水体变为草地、建筑物变为裸地等两种变化具有较好区分能力,Shv与回波中的去极化分量有关,反映了地物的去极化能力[27]。水体的后向散射是比较纯粹的镜面反射(奇次散射),裸地也以奇次散射占优,因此在水体变为草地、建筑物变为裸地这两种变化中均表现出hv通道回波强度的明显变化。从图 5中可以发现,草地到建筑物的变化在hv通道上也几乎可以分离开。图 5同时表明,草地到建筑物的变化在hh通道上可以得到很好的区分。这主要是因为建筑物后向散射的中的绝大部分为偶次散射(二面角散射),因此其hh通道的回波强度远大于vv及hv通道的回波强度。
进一步地,设定dc(T3)、dspan、dshh、dshv、dsvv等5个图像只分为两类:变化类和未变化类,在整幅影像中选取了40个像元,其中20个为变化像元,20个为未变化像元,按式(11)计算每幅图像的类别分离度(divergence)[28]。类别分离度越大说明两类之间的可分性高。
式中σi、σj分别为两类的方差;μi,μj为其均值。
表 3定量地显示出,与基于后向散射强度生成的各种差异图像相比,极化距离图像能更好地区分变化。以极化距离图像dc(T3)为基础(图 7(a)),利用基于Markov随机场模型方法提取了最终的变化像元,结果如图 7(b)所示。为了对最终变化检测结果进行定量评价,选择图中一块区域,在光学影像辅助下,手工提取了所有的变化图斑,作为参考影像。图 8(a)为图 7(b)白框区域的放大,与其对应的变化参考影像显示于图 8(b)。对分别由dc(T3)、dspan、dshh、dshv、dsvv生成的变化图斑(均通过Markov随机场模型方法提取变化图斑)进行了精度评价,结果列于表 4,可见由极化距离测度提取的变化图具有最低的虚警率和漏检率,总体变化检测精度达到94.54%(1-5.46%)。
(%) | |||
差异影像 | 虚警率 | 漏检率 | 总体误差 |
dc(T3) | 2.51 | 4.25 | 5.46 |
dspan | 6.47 | 7.67 | 11.47 |
dshh | 6.52 | 16.01 | 17.79 |
dshv | 8.21 | 8.77 | 13.77 |
dsvv | 4.37 | 22.93 | 21.52 |
在基于极化协方差矩阵分类的Wishart距离基础上,本文提出了适用于多时相极化SAR影像变化检测的计算简便的距离测度,并给出了相应的完整的极化变化提取技术流程。利用北京地区2009年和2010年获取的两景RadarSat-2全极化影像进行了试验,通过与基于SAR后向散射强度得到的多种差异图像进行对比,极化距离测度能反映全部不同类型的地物变化,而且在极化距离图像上变化与未变化区域的分离度更高。这说明,本文提出的极化距离测度在变化检测应用中具有很好的潜力。
[1] | RIGNOT E J M,VAN ZYL J J.Change Detection Techniques for ERS-1 SAR Data[J].IEEE Transactions on Geoscience and Remote Sensing,1993,31(4):896-906. |
[2] | LU D,MAUSEL P,BRONDIZIO E,et al.Change Detection Techniques[J].International Journal of Remote Sensing, 2004, 25(12): 2365-2407. |
[3] | CIHLAR J,PULTZ T J,GRAY A L.Change Detection with Synthetic Aperture Radar[J].International Journal of Remote Sensing, 1992, 13(3): 401-414. |
[4] | VILLASENOR J D, FATLAND D R, HINZMAN L D. Change Detection on Alaska’s North Slope Using Repeat Pass ERS-1 SAR Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 1993, 31(1): 227-236. |
[5] | MOSER G, SERPICO S B. Generalized Minimum Error Thresholding for Unsupervised Change Detection from SAR Amplitude Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(10): 2972-2982. |
[6] | KIM D J,MOON W M,KIM Y S.Application of TerraSAR-X Data for Emergent Oil-spill Monitoring[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(2): 852-863. |
[7] | DEKKER R J.Speckle Filtering in Satellite SAR Change Detection Imagery[J]. International Journal of Remote Sensing, 1998, 19(6): 1133-1146. |
[8] | BOVOLO F, BRUZZONE L.A Detail-preserving Scale-driven Approach to Change Detection in Multitemporal SAR Images[J]. IEEE Transactions on Geoscience and Remote Sensing,2005,43(12):2963-2972. |
[9] | HUANG Shiqi,LIU Daizhi,HU Mingxing,et al.Multitemporal 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.) |
[10] | WAN Honglin,JIAO Licheng,XIN Fangfang.Interactive Segmentation Technique and Decision Level Fusion Based Change Detection for SAR Images[J].Acta Geodaetica et Cartographica Sinica,2012,41(1):74-80.(万红林,焦李成,辛芳芳.基于交互式分割技术和决策级融合的SAR图像变化检测[J].测绘学报,2012,41(1):74-80.) |
[11] | WAN Honglin,JIAO Licheng,WANG Guiting,et al.A Region-of-interest Level Method for Change Detection in SAR Imagery[J].Acta Geodaetica et Cartographica Sinica,2012,41(2):239-245.(万红林,焦李成,王桂婷,等.在感兴趣的区域层面上进行SAR图像变化检测的方法研究[J].测绘学报,2012,41(2):239-245.) |
[12] | BARONTI S,CARLA R,SIGISMONDI S,et al.Principal Component Analysis for Change Detection on Polarimetric Multitemporal SAR Data[C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposium. Pasadena: IEEE, 1994:2152-2154. |
[13] | KERSTEN P R,LEE J S,AINSWORTH T L.A Comparison of Change Detection Statistics in POLSAR Images[C]//Proceedings of 2005 IEEE International Geoscience and Remote Sensing Symposium. New York: IEEE, 2005: 4836-4839. |
[14] | SABRY R. A New Coherency Formalism for Change Detection and Phenomenology in SAR Imagery:A Field Approach[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(3): 458-462. |
[15] | 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): 526-532. (郝洪美,张永红,石海燕,等.统计假设检验方法在全极化SAR变化检测中的应用[J].遥感学报, 2012, 16(3): 526-532.) |
[16] | WU Fan,WANG Chao,ZHANG Hong.Residential Areas Extraction in High Resolution SAR Image Based on Texture Features[J].Remote Sensing Technology and Application,2005,20(1):148-152. (吴樊,王超,张红.基于纹理特征的高分辨率SAR影像居民区提取[J]. 遥感技术与应用, 2005, 20(1): 148-152.) |
[17] | CONRADSEN K,NIELSEN A A,SCHOU J,et al.Change Detection in Polarimetric SAR Data and the Complex Wishart Distribution[C]//Proceedings of IEEE 2001 International Geoscience and Remote Sensing Symposium. Sydney: IEEE, 2001: 2628-2630. |
[18] | CONRADSEN K,NIELSEN A A,SCHOU J,et al.A Test Statistic in the Complex Wishart Distribution and Its Application to Change Detection in Polarimetric SAR Data[J].IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(1): 4-19. |
[19] | LEE J S,POTTIER E. Polarimetric Radar Imaging:From Basics to Applications[M].Boca Raton: CRC Press, 2009. |
[20] | LEE J S,GRUNES M R,KWOK R. Classification of Multi-look Polarimetric SAR Imagery Based on the Complex Wishart Distribution[J]. International Journal of Remote Sensing, 1994, 15(11): 2299-2311. |
[21] | FAMIL L F, POTTIER E, LEE J S. Unsupervised Classification of Multifrequency and Fully Polarimetric SAR Images Based on the H/A/Alpha-Wishart Classifier[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(11): 2332-2342. |
[22] | LEE J S.Unsupervised Terrain Classification Preserving Polarimetric Scattering Characteristics[J].IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(4): 722-731. |
[23] | MOSER G, SERPICO S, VERNAZZA G. Unsupervised Change Detection from Multichannel SAR Images[J]. IEEE Geoscience and Remote Sensing Letters, 2007, 4(2): 278-282. |
[24] | KASETKASEM T, VARSHNEY P K. An Image Change Detection Algorithm Based on MARKOV Radom Filed Models[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1815-1823. |
[25] | MDA. RADARSAT-2 Product Format Definition: RN-RP-51-2713: Issue 1/7[R]. Richmond: MDA, 2008. |
[26] | LIU Guoxiang,DING Xiaoli,LI Zhilin,et al.Co-registration of Satellite SAR Complex Images[J].Acta Geodaetica et Cartographica Sinica,2001, 30 (1): 60-66.(刘国祥,丁晓利,李志林,等.星载SAR复数图像的配准[J].测绘学报,2001, 30 (1): 60-66.) |
[27] | ZEBKER H A, VAN ZYL J J. Imaging Radar Polarimetry: A Review[J]. Proceedings of IEEE, 1991, 79 (11): 1583-1606. |
[28] | SWAIN P H,DAVIS S M.Remote Sensing:The Quantitative Approach[M].New York:McGraw Hill Book Company, 1978. |