面向高分辨率极化SAR四分量分解的相干矩阵自适应估计

杨　帅, 陈启浩, 刘修国, 陈　奇

(中国地质大学(武汉)信息工程学院, 湖北 武汉 430074)



面向高分辨率极化SAR四分量分解的相干矩阵自适应估计

杨帅, 陈启浩, 刘修国, 陈奇

(中国地质大学(武汉)信息工程学院, 湖北 武汉 430074)

摘要：提出了一种相干矩阵的自适应估计方法,基于估计结果矩阵进行传统四分量分解,以提高高分辨率极化合成孔径雷达(synthetic aperture radar, SAR)四分量分解的准确性。该方法根据像元主要散射类型对估计样本做自适应选取,消除异质像元对估计结果的影响;并通过对地物异质度的判别实现窗口平均估计和不动点迭代估计的自适应选择,以适应不同异质程度的地物分布情况。实验采用ESAR以及UAVSAR数据进行验证,结果表明,采用自适应矩阵估计的四分量分解结果更接近地物的实际散射机理。

关键词：极化合成孔径雷达; 四分量分解; 高分辨率; 矩阵估计; 自适应

0引言

目标分解是极化合成孔径雷达(synthetic aperture radar, SAR)图像解译、地物散射机理分析的一种重要手段[1]。基于模型的目标分解方法能得到符合地物实际散射机理的模型描述,广泛用于地表反演、目标识别以及图像分类等方面[2-3],因此一直以来是极化雷达领域的研究焦点之一[4-5]。

文献[6]基于反射对称假设提出的三分量分解,通过建立表面散射、二面角散射和体散射3种散射机制的模型,分别用于模拟粗糙表面、直角反射器以及森林冠层等地物。在此基础上文献[7]提出了四分量分解,通过引入螺旋体散射分量,来适应复杂人造地物反射不对称的情况。为了进一步提高四分量分解的准确性,国内外学者通过替换散射模型、优化模型参数、去取向等进行综合改进[8-12]。然而这些方法主要针对模型和流程进行修改,对四分量分解的源头对象——相干矩阵本身,未进行准确估计的考虑。

对中低分辨率极化SAR图像而言,由于受表面散射和体散射的随机矢量散射效应以及相干斑噪声的综合影响,传统的四分量分解方法需要采用窗口集合平均法(ensemble average, EA)估计目标的相干矩阵[13],再分解之。然而,高空间分辨率极化SAR图像呈现出新的特点:其一是分辨单元内子散射体数目明显减少,地物目标的异质性增大[14],使得大量地物目标的统计特性不再满足高斯模型,窗口平均法得到的结果不再是最优估计;其二是地物目标的空间信息高度细节化[15],导致矩形窗口边界与地物目标边界差异更为显著,此时窗口平均的矩阵估计方式会造成目标细节信息的严重损失[16-17]。因此,传统四分量分解采用的窗口平均估计不再适用于高分辨率SAR影像。

为此,本文提出了一种面向高分辨率极化SAR四分量分解的相干矩阵自适应估计方法。首先在矩形窗口的基础上,剔除部分与中心像元散射类型不一致的像元,将其作为自适应的估计样本取代原有的矩形窗口样本;然后对估计方式做自适应的选择,即通过对地物异质度的判别,在均质区保留窗口平均估计,在异质区引入文献[18]提出的不动点迭代估计(fixed point, FP);之后用自适应估计的结果进行四分量模型分解。

1极化SAR相干矩阵估计方法

传统四分量分解一般采用窗口平均的矩阵估计方法,对窗口内所有像元进行平均运算获得相干矩阵估计结果,即EA估计。这种估计方法适用于符合高斯模型的均质区,但用于高分辨率极化SAR图像中的异质区,会损失细节信息。

为此,文献[19]引入了球不变随机相量(spherically invariant random vector, SIRV)模型来对异质区进行建模。该模型中纹理随机变量的概率分布函数并不是被明确定义的,相较于其他模型,其通用性更强。文献[18]结合SIRV模型以及高分影像存在的特点,进一步提出了相干矩阵的FP估计,使用功率P和归一化相干矩阵M的组合来表示原来的相干矩阵。归一化矩阵M仅仅表示极化信息,而功率部分P则表示纹理随机变量。

归一化矩阵M的估算公式为

(1)

(2)

式中,Ti表示目标的相干矩阵;tr()表示矩阵的迹;“-1”表示矩阵求逆;N是样本像元的数目;m是归一化相干矩阵的维度。对M矩阵的计算采用迭代的方式进行。

对于功率部分P的计算,有两种方法:只考虑中心像元计算时,极化白化滤波(polarimetricwhiteningfilter,PWF)的功率为PPWF,如式(3)所示;利用窗口内所有像元计算时,多视极化白化滤波(multilookpolarimetricwhiteningfilter,MPWF)的功率为PMPWF,如式(4)所示。

(3)

(4)

由式(1)~式(4),得到相干矩阵估计的最终结果:

(5)

2面向四分量分解的相干矩阵自适应估计

针对高分辨率极化SAR影像,本文从估计样本和估计方式两个方面对原有的矩阵估计方法进行修改,提出一种新的相干矩阵自适应估计方法。并把估计结果矩阵用于传统四分量分解,从而获得更加准确的分解结果。

2.1估计样本的自适应选取

在矩阵估计时,无论是窗口平均估计还是FP估计,对图像的处理方式都是矩形估计窗口滑动遍历整幅影像。然而在高分辨率SAR影像中,地物微观信息更加丰富,导致原有的矩形估计窗口边界与地物目标边界差异更为显著,使用窗口内所有像元进行估计得到的很可能不再是最优的结果。

为此,在矩阵估计前对估计的样本进行一次自适应选取。考虑到分解结果的一致性,选用预分解的散射类型结果作为估计样本选择的参考依据。为了保留较多的估计样本以及对估计结果尽量小的改动,采用一种新的样本选择方式:在窗口范围内剔除异质像元,异质像元定义为与中心待估计像元主要散射类型不一致的像元。

首先,按照散射类型对像元进行分类:使用传统Yamaguchi分解方法进行分解,求出每一种功率所占像元总功率的比值。当功率最大的散射机制功率百分比超过50%时,认为该散射机制就是对应像元的主要散射机制,否则认为该像元是混合散射机制。处理后,整幅影像被分为表面散射像元、二面角散射像元、体散射像元、螺旋体散射像元和混合散射像元。

之后根据像元主要散射类型进行异质像元的剔除,具体包括:①提取中心像元的主要散射类型,对矩形窗口内所有像元进行遍历,并把与中心像元散射类型不匹配的像元剔除。②对于遍历到的混合像元,认为与中心像元类型基本相同,不做处理。③对于中心像元为混合散射机制像元的情况,不进行异质像元的剔除。④对于剔除像元超过50%的特殊情况,为避免估计样本不足导致估计结果不稳定,采用原始矩形窗口进行估计。

2.2估计方式的自适应选择

在均质区使用窗口平均估计能够得到比较准确的相干矩阵估计结果,但在异质区则是FP估计更优[18]。为了兼顾高分辨率极化SAR影像中的均质区与异质区,本文提出一种自适应的矩阵估计方式:通过待估计样本的异质程度来判别均质区和异质区,在均质区保留原本的窗口平均估计以及PMPWF,在异质区则引入FP估计以及PPWF。

实现上述自适应的矩阵估计方式,关键在于确定一个能够判别待估计像元异质程度的指标。对于单极化SAR图像,地物目标的异质度可以通过变差系数即标准差与均值的比值来反映。本文提出使用指标I作为异质度判别的指标,计算公式为

(6)

式中,var表示求方差。I是由变差系数推广至全极化SAR的,可以根据极化白化滤波[20],将矩阵变量C映射为标量x:

(7)

式中,Σ表示总体的协方差矩阵,可通过渐进最大似然估计法估算。变换后的变量x与原矩阵变量具有一致的统计特性,且x的均值为矩阵C的维度。

2.3基于自适应相干矩阵估计的四分量分解

采用动力学求解器进行仿真，总仿真时长为0.000 2 s。计算结束后分析导线仿真模型对于方案2设定的外部加载的响应，提取芯棒中心单元绘制应力时程曲线，如图10所示。由图10可以看出：在加载过程开始后，芯棒、外层铝股的应力快速上升，很快达到最大值，即达到最大应力状态；而后随着时间推移，铝股间出现相互作用，应力出现震荡并逐渐达到稳定。

基于自适应相干矩阵估计的四分量分解流程如图1所示。主要过程包括:①用Yamaguchi分解方法对图像进行一次预分解,使用第2.1节的方法自适应地确定每个待估计像元的估计样本。②把待估计的相干矩阵分拆成两部分表示:归一化相干矩阵M和功率P。③计算待估计像元的I值,针对不同的I值做自适应选择判别:如果I大于自适应阈值I0,使用FP估计方法来计算归一化矩阵MFP,同时由式(3)计算得到功率PPWF;否则保留原来的窗口平均估计方法获得归一化矩阵MEA,然后由式(4)计算得到功率PMPWF。④把估计后的归一化相干矩阵M和功率P合并,并把最终的相干矩阵估计结果用于传统四分量分解。

图1　利用自适应相干矩阵估计的四分量分解流程

3实验与分析

3.1实验数据说明

实验数据采用L波段ESAR的德国Oberpfaffenhofen地区的全极化影像,影像大小为1 200像素×1 000像素,单视数据,空间分辨率为1.5m×0.89m(距离向×方位向)。实验区的PauliRGB影像如图2(a)所示,图2(b)是该区域的光学影像,成像时间与极化SAR成像时间较为接近。可以看出该区域主要包含有密集的林地、建筑、农田、机场跑道等地物,跑道附近分布有一些零散的人工建筑物。

3.2异质度判别指标性能验证

计算实验区域影像的I值,结果如图3所示,图中颜色越趋近于白色表示I值越大。可以发现,整体而言在城区等异质度较高的区域,I值明显要比异质度较低的农田大,因此可以将I值作为异质度的判别指标。

图2　实验区域PauliRGB及光学影像图

图3　实验区域I值图

进一步,为了确定异质程度的判别阈值I0,在影像中选取了一些样本比较纯的农田、跑道、森林和城区等典型地物并计算I值,如表1所示。

表1　典型地物I值统计结果(EASR)

由表1可以看出,对于不同的地物类型,其I值的均值不同。农田、跑道等异质度较低的区域I值基本在1.6以下,而森林等I值则在1.75以上。对于极端异质的城区,I值则显著增大到2.5以上。一般认为森林、城区等是异质区,而跑道、农田等则是均质区。本文选择1.7作为区分均质区和异质区的自适应阈值I0。

为了验证不同的矩阵估计方法对四分量分解结果的影响,分别用以下4种估计方法进行处理:①窗口平均估计；②FP估计；③自适应矩阵估计(矩形窗口);④自适应矩阵估计(自适应估计窗口),即本文方法。在估计窗口大小的选择上,过小的窗口会导致估计算法不稳健,估计效果不明显;过大的窗口又会丢失图像的细节信息,影响估计结果,实验选定7×7窗口。

图4(b)~图4(d)分别为方法①、方法②和方法④的分解结果三分量RGB合成图,其中红色表示二面角散射、绿色表示体散射,蓝色表示表面散射。整体而言,引入FP估计之后,方法②、方法④的分解结果明显减少了窗口平均的影响,很多区域细节更为清晰。

图4　实验区域影像及四分量分解结果

进一步放大观察几个典型区域对比如图5所示,可以看出在农田、道路等均质区域,方法①能有效去除相干斑的影响,分解结果比较平滑,方法②的分解结果则趋于斑点状,与相对单一均匀的实际地物分布不符;在城区、森林等异质区域,方法①的分解结果比较模糊,丢失了很多细节信息,方法②的结果则较好地保留了建筑物边界、结构等信息。

方法④在均质区域能够保留窗口平均估计的平滑滤波效果,但又不至于过度平滑,如图5(d)中农田内部有一些异质目标(绿色部分)的分解细节得到了保留;而在异质区域,方法④由于引入了FP估计的优良特性,分解结果与实际地物类型更为接近,细节信息保留较好。

图5　ESAR分解结果局部对比图

为进一步验证本文第2.1节提出的自适应估计样本对分解结果的影响,增加方法③、方法④的定量对比实验。选择如图6(a)所示点目标区域作切片,该区域中心位置为一个亮点目标,主要散射类型是二面角,点周围存在着草地、机场跑道等地物像元。首先选择以点目标像元为主体的矩形区域,统计该区域内分解结果功率的平均值如表2所示;然后在该矩形区域周围选择如图6(b)所示的环形区域,同样统计该区域内的分解结果功率的平均值如表3所示。

表2　点目标区域分解结果功率平均值统计表

表3　环形区域分解结果功率平均值统计表

图6　点目标统计区域示意图

由表2可以看出,方法④相比方法③,二面角散射分量得到了增强,其他3个分量则略有下降。这与点目标的主要散射机制相符,说明经过估计样本的自适应选取,该区域的主要散射机制得到了合理的保留和增强。然后对比表2与表3中的总功率变化情况,环形区域的总功率略微下降而中心区域的总功率提高。这表明在采取了自适应的估计样本之后,周围像元对亮点区域像元进行估计的干扰减少,中心像元的功率被合理增强;而周围像元在估计时也减少了中心亮点像元的干扰,使得整体功率下降。本文采用的估计样本自适应途径是异质像元的剔除,整体而言剔除像元的数量较少,因此在总功率变化方面应该是比较小的,这与表2、表3中总功率变化的幅度较小也是一致的。

为了进一步验证本文方法的有效性,选取跑道附近的点目标区域做剖面分析。该区域如图7所示:中心位置为一个亮点目标,点目标右侧为异质度较高的黑色区域,其余部分则是跑道等均质区域。分别做如下三组分解结果功率剖面差值图:③-①、③-②、④-③,即把对应的分解结果各分量求差,然后从差值图上对比两种方法对分解结果的影响,结果如图8所示。红色表示二面角分量,绿色表示体散射分量,蓝色表示表面散射分量,黑色表示螺旋体散射分量。

图7　剖面区域示意图(画线处为剖面)

由图8(a)、图8(b)可以看出,方法③的分解结果在跑道等均质区与方法①的比较一致,在点目标等异质区与方法②的基本一致。这说明通过本文的自适应选择,方法③能够兼顾均质区和异质区,分别保留两种估计各自的优点。比较方法④与方法③的分解结果,发现在点目标像元处,二面角散射分量被显著加强,而这与该区域的主要散射类型是相符的。说明通过估计样本的自适应选择,该区域的主要散射功率得到了合理增强。

图8　跑道附近点目标功率剖面对比(功率差值)

3.4方法通用性验证

为了验证本文方法对不同平台极化SAR影像的有效性,选择L波段UAVSAR美国Haywrd地区的全极化影像进行实验。影像大小为2 500像素×2 500像素,空间分辨率为4.997 m×7.2 m(距离向×方位向),成像时间为2014年2月,如图9(a)所示。图9(b)为该地区对应时期的光学影像。对实验影像计算I值,如图10所示,图中颜色越趋近于白色表示I值越大。

选取一些典型地物样本统计其I值分布,结果如表4所示。可以看出,对于UAVSAR而言,其I值分布与ESAR略有不同,但总体趋势是一致的,这说明I值仍能较好地反映地物的异质情况。选择1.5作为I0,并采用4.3节中方法①、方法②、方法④进行实验。UAVSAR的影像成像质量较好,因此矩阵估计的窗口大小定为5×5。

图9　UAVSAR实验区域及分解结果图

图10　UAVSAR实验区域I值分布图

地物编号最小值最大值均值标准偏差河流0.360.700.510.07荒山0.333.190.770.36农田0.421.680.820.20森林0.653.941.700.48城区1.088.892.941.02

对分解结果局部细节进行对比分析,如图11所示。可以看出,本文方法对于异质度不同的区域适应性更好,能够在保留一些异质对象的分解结果细节的同时,降低均质区域的相干斑对分解结果的影响。这一结果与ESAR数据分解所得到的结论是类似的。

图11　UAVSAR分解结果局部对比图

4结论

本文提出的面向高分辨率极化SAR四分量分解的自适应相干矩阵估计方法,根据像元主要散射类型对样本做自适应的选取,避免了少数异质像元对分解结果的影响;同时通过对地物异质度的判别实现窗口平均估计和FP估计的自适应选择。实验表明,本文方法能够增强高分影像,尤其是异质度较高区域四分量分解的准确性。把自适应矩阵估计方法和模型分解的综合改进方法进行结合是下一步工作的重点。

参考文献：

[1] Sugimoto M, Ouchi K, Nakamura Y. On the novel use of model-based decomposition in SAR polarimetry for target detection on the sea[J].RemoteSensingLetters, 2013, 4(9): 843-852.

[2] Alberga V. A study of land cover classification using polarimetric SAR parameters[J].InternationalJournalofRemoteSensing, 2007, 28(17): 3851-3870.

[3] Han P, Xu J S, Zhao A J. PolSAR image runways detection based on multi-stage classification[J].SystemsEngineeringandElectronics, 2014, 36(5): 866-871. (韩萍, 徐建飒, 赵爱军. 基于多级分类的 PolSAR 图像机场跑道检测[J].系统工程与电子技术, 2014, 36(5): 866-871.)

[4] Van Beijma S, Comber A, Lamb A. Random forest classification of salt marsh vegetation habitats using quad-polarimetric airborne SAR, elevation and optical RS data[J].RemoteSensingofEnvironment, 2014, 149: 118-129.

[5] Lee J S, Ainsworth T L, Wang Y T. Generalized polarimetric model-based decompositions using incoherent scattering models[J].IEEETrans.onGeoscienceandRemoteSensing,2014,52(5):2474-2491.

[6] Freeman A, Durden S L. A three-component scattering model for polarimetric SAR data[J].IEEETrans.onGeoscienceandRemoteSensing, 1998, 36(3): 963-973.

[7] Yamaguchi Y, Yajima Y, Yamada H. A four-component decomposition of POLSAR images based on the coherency matrix[J].IEEEGeoscienceandRemoteSensingLetters,2006,3(3):292-296.

[8] An W T, Xie C H. An improvement on the complete model-based decomposition of polarimetric SAR data[J].IEEEGeoscienceandRemoteSensingLetters, 2014, 11(11): 1926-1930.

[9] Cui Y, Yamaguchi Y, Yang J, et al. On complete model-based decomposition of polarimetric SAR coherency matrix data[J].IEEETrans.onGeoscienceandRemoteSensing, 2014, 52(4): 1991-2001.

[10] Chen S W, Wang X S, Xiao S P, et al. General polarimetric model-based decomposition for coherency matrix[J].IEEETrans.onGeoscienceandRemoteSensing,2014,52(3):1843-1855.

[11] Liu X G, Huang X D, Chen Q H, et al. Three component model-based decomposition integrating de-orientation and generalized volume scattering model[J].JournalofElectronics&InformationTechnology,2012,34(10):2451-2457.(刘修国,黄晓东,陈启浩,等.综合去取向和广义体散射的三分量极化目标分解模型[J].电子与信息学报,2012,34(10):2451-2457.)

[12] Shan Z L, Wang C, Zhang H, et al. Improved four-component model-based target decomposition for polarimetric SAR data[J].IEEEGeoscienceandRemoteSensingLetters,2012,9(1):75-79.

[13] Wang C, Zhang H, Chen X, et al.Imageprocessingofpolarimetricsyntheticapertureradar[M].Beijing:Science Press, 2008.(王超,张红,陈曦,等.全极化合成孔径雷达图像处理[M].北京:科学出版社,2008.)

[14] 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].IEEETrans.onGeoscienceandRemoteSensing, 2004,42(10):2046-2057.

[15] Vasile G, Pascal F, Ovarlez J, et al. Optimal parameter estimation in heterogeneous clutter for high-resolution polarimetric SAR data[J].IEEEGeoscienceandRemoteSensingLetters, 2011, 8(6): 1046-1050.

[16] Alberga V, Krogager E, Chandra M, et al. Potential of coherent decompositions in SAR polarimetry and interferometry[C]∥Proc.oftheGeoscienceandRemoteSensingSymposium,2004:1792-1795.

[17] Lee J S, Pottier E.Polarimetricradarimaging:frombasicstoapplications[M].CRC Press, 2009.

[18] Vasile G, Ovarlez J, Pascal F, et al. Coherency matrix estimation of heterogeneous clutter in high-resolution polarimetric SAR images[J].IEEETrans.onGeoscienceandRemoteSen-sing, 2010, 48(4): 1809-1826.

[19] Yao K. A representation theorem and its applications to spherically-invariant random processes[J].IEEETrans.onInformationTheory, 1973, 19(5): 600-608.

[20] Novak L M, Burl M C. Optimal speckle reduction in polarimetric SAR imagery[J].IEEETrans.onAerospaceandElectronicSystems, 1990, 26(2): 293-305.

杨帅(1991-),男,博士研究生,主要研究方向为极化SAR信息提取与应用。

E-mail: cug_ys@163.com

陈启浩(1982-),男,讲师，博士,主要研究方向为极化SAR信息提取。

E-mail: cugcqh@163.com

E-mail: liuxg318@163.com

陈奇(1989-),男,硕士研究生,主要研究方向为极化SAR信息提取。

E-mail: chenqi106@foxmail.com

网络优先出版地址：http://www.cnki.net/kcms/detail/11.2422.TN.20150910.1130.012.html

Adaptive coherent matrix estimation method for high resolution

polarimetric SAR four-component decomposition

YANG Shuai, CHEN Qi-hao, LIU Xiu-guo, CHEN Qi

(CollegeofInformationEngineering,ChinaUniversityofGeosciences,Wuhan430074,China)

Abstract:To improve the accuracy of high resolution polarimetric synthetic aperture radar (SAR) four-component decomposition, an adaptive coherent matrix estimation method is proposed, which is applied before decomposition. First, in order to eliminate the effect of the extreme heterogeneity cell, proper samples are chosen adaptively according to the main types of cell scattering. Then, the ensemble average estimate and the fixed point estimate are applied separately according to the degree of feature heterogeneous to accommodate different feature distribution. This method is validated by using the ESAR data and the UAVSAR data. The results show that the decomposition results are closer to real scattering types by applying the adaptive matrix estimation.

Keywords:polarimetric synthetic aperture radar (SAR); four-component decomposition; high resolution; matrix estimation; adaptive

通讯作者刘修国(1969-),,男,教授,博士研究生导师,主要研究方向为遥感影像信息提取与3S集成。

作者简介：

中图分类号：TN 957.52

文献标志码：A

DOI:10.3969/j.issn.1001-506X.2016.01.09

基金项目：国家自然科学基金(41471355，41301477)；中国博士后科学基金面上项目(2012M521497)资助课题

收稿日期：2015-03-11；修回日期：2015-06-11；网络优先出版日期：2015-09-10。