基于曲波变换的全极化SAR数据与TM影像融合方法

肖世忱，廖静娟

(1.中国科学院遥感与数字地球研究所数字地球重点实验室，北京 100094；2.中国科学院大学，北京 100049)

1 引 言

随着传感器技术、遥感技术等的迅猛发展，多平台、多传感器、多时相、多光谱和多分辨率的遥感影像构成了同一地区的多源遥感信息。在众多的信息源中，SAR数据对目标的几何特性非常敏感，而多光谱影像对目标的物理和化学属性则更为敏感，主要反映不同地物的轮廓与光谱信息[1]，这些不同的信息类型涉及到完全不同的目标质量并且极大地不相关[2]，因此SAR数据与多光谱影像的融合将能充分利用它们的互补性，弥补各自的不足，这对于更好地实现对被观测目标的理解是非常重要的。

SAR数据与多光谱影像进行融合的方法主要有IHS变换和小波变换等。其中基于IHS变换的经典融合方案是在IHS空间内，使用经过拉伸处理的单极化SAR数据替换多光谱影像的I分量的方式来产生融合影像[3]，该方式虽容易扭曲原始光谱特性，但可分离多光谱影像中的光谱信息和强度信息。小波变换因其具有良好的多尺度与时频局部特性在影像融合领域也得到广泛应用。近年来Candes等提出曲波变换的概念[4]，该变换不仅继承了小波变换的优良特性，还具有多方向和各向异性的特征，在给定相同的重构精度下可以近乎最优地表示影像中具有方向性的直线及曲线等奇异性特征，对函数的稀疏表达能力较小波变换也更强。基于Curvelet变换的经典融合方案是在粗尺度层取系数平均值，细尺度层取绝对值较大的系数。许星等采用这种方案进行SAR数据与TM影像的融合时，获得了相比于小波变换更好的融合效果[5]。

目前，研究者在进行SAR数据与多光谱影像的融合时，多采用单极化数据，主要因为SAR数据记录的后向散射信息对地面目标的类型、方向、同质性和空间关系更为敏感，相比于光学遥感影像呈现出有望能被更好使用的结构特征[6]。相比于单极化SAR数据，多极化数据在具备目标丰富结构信息的同时，还记录了目标在不同极化方式下的信息，尤其是全极化数据，以散射矩阵的方式记录了地面目标的信息。随着Radarsat-2等卫星的成功发射，全极化SAR数据逐渐丰富，其参与融合可使融合结果中包含全极化数据中蕴含的丰富地物信息。Xu等尝试在RGB空间采用相应通道以一定的权重相加来实现SAR全极化数据与多光谱影像的融合[7]，但由于RGB颜色空间是根据红、绿和蓝3原色的混合比例来定义不同的颜色，而人眼只能通过感知颜色的色调、亮度和饱和度来识别和区分物体[8]，因此选择在IHS空间进行融合将更加优越。

本文在IHS空间提出了一种基于第二代曲波变换融合全极化SAR数据与TM影像的方法。针对两种数据源的影像特征，设计了一种充分保留原影像重要特征的低频子带系数融合规则。实验结果表明，该方法所得融合影像较基于IHS变换的经典融合方案和基于曲波变换的经典融合方案所得影像在主观视觉效果和客观定量指标上均有所改善。

2 数据源与数据处理

2.1 数据源

SAR数据源是Radarsat-2卫星于2012年6月中旬在黑龙江省获得的全极化数据。产品类型是单视复数据，空间分辨率为8m。TM影像是由Landsat-5卫星于2008年6月30日获取的，产品等级为L4，包含7个波段，分辨率为30m。影像内大部分区域地势平坦，包括居民区、耕地、河流和道路等典型地类。

2.2 数据处理

Radarsat-2全极化数据与TM影像融合前的数据处理主要涉及SAR数据的相干斑抑制、影像配准以及直方图匹配等环节。由于TM影像的L4级产品是经过辐射校正、几何校正和几何精校正等处理后的高精度产品，因此数据处理环节中涉及到信息统一的操作将以TM影像为基准。

SAR数据的相干斑抑制步骤是先对SAR数据进行方位向：距离向等于2∶1的多视处理来解决方位向和距离向采样大小差异较大的问题，再采用窗口大小为3×3的Gamma Map自适应滤波方法分别对4种极化数据进行降斑处理，处理后地物内部平滑均匀，地物边缘锐化清晰。配准处理是以TM影像为基准来校正SAR数据，控制点的选取位于道路和田梗等的交叉点和拐点处，纠正方法选用的是一阶多项式。配准后X、Y方向的配准误差均控制在1个像元内，总体RMS为0.60个像元。由于从不同传感器采集到的影像会因平台角度、影像获取时间以及卫星波长扩展间存在差异等造成相同的场景具有不同照明条件[9]，所以融合前还需将SAR数据各通道同TM影像相应波段进行直方图匹配来调整SAR数据的亮度。

3 融合方法

曲波(Curvelet)变换理论可细分为两代，其中第二代Curvelet变换[10]是一种更易于数字实现、更便于理解的快速变换，与第一代Curvelet变换[4]相比，在构造上已经完全不同，其在实现过程中没有引入Ridgelet变换[11]，而是在频域中直接给出了Curvelet基的具体表示形式，可以说是一种真正意义上的Curvelet变换。第二代Curvelet变换算法包括两种具体实现：缠绕算法和非等间隔快速傅里叶算法[12]。相比于后者，缠绕算法的等间隔特性使其更容易实现且忠实于连续变换，再加上运算量和复杂度都相对较小[13]，因此本文将基于缠绕算法对数据进行Curvelet变换。SAR数据与TM影像融合方案如图1所示，I分量Curvelet变换流程如图2所示。

3.1 低频子带系数融合规则

影像经Curvelet变换后所得低频子带系数代表的是原影像的概貌信息，包含了影像的主要能量。由于SAR与TM影像的成像机理不同，同一目标在两种影像中的灰度分布特性往往存在较大差异，甚至极性完全相反。因此，对低频子带系数沿用传统加权平均的融合规则通常不能获得满意的结果。为此本文提出了一种基于SAR数据与TM影像特征的低频子带系数融合规则。

Radarsat-2数据的分辨率高于TM影像，在刻画目标边缘等细节信息方面的能力要优于TM影像，因此首先对SAR数据的低频子带系数求取梯度并作归一化处理，然后适当选取梯度阈值TG，当SAR低频子带系数的梯度值≥TG时，融合系数F(x，y)取SAR与TM低频子带系数的线性加权。由于SAR数据中亮度较大的目标往往也是需要融合的目标，而TM影像在描述影像纹理特征方面要优于SAR数据，因此当SAR低频子带系数的梯度

影像f(x，y)在像素点(x，y)处的梯度可定义为一个向量：

(1)

梯度的模为各分量的平方和再求平方根，即为：

(2)

对于数字图像，连续导数形式可以用求差来近似表示，则：

|G(x,y)|≅|f(x,y)-f(x+1,y)|+|f(x,y)-f(x,y+1)|

(3)

低频子带系数的融合规则为：

(4)

3.2 高频子带系数融合规则

经Curvelet变换后所得高频子带系数描述的是影像的细节和边缘等特征，而这些特征在曲波变换域中则表现为具有绝对值较大的高频子带系数。因此，高频子带系数采用取绝对值大的系数进行融合。

图1 SAR与TM影像的融合方案

图2 基于Curvelet变换的影像I分量融合流程图

4 融合结果与效果的分析评价

融合结果的总体效果如图3所示。由于试验区较大，总体效果无法详细描述细部信息，文中又选取道路和耕地结果放大图进行具体说明，如图4、图5所示。每幅图中又包括6景子图像，其中(a)为原SAR影像，(b)为原TM影像，(c)为HV极化数据基于曲波变换经典融合方案的结果，(d)为全极化SAR数据基于IHS变换经典融合方案的结果，(e)为全极化SAR数据基于曲波变换经典融合方案的结果，(f)为文中方法的融合结果。

图3 总体效果图

图4 道路效果图

图5 耕地效果图

从目视判读的效果上看，相比于影像(c)，影像(e)和(f)不仅影像内的信息更加丰富，而且影像效果更加平滑，这说明全极化SAR数据产生的I分量相比于SAR单极化通道不仅能有效综合各极化通道的信息，还能进一步抑制单极化通道中的相干斑噪声，从而产生更好的融合效果。对比影像(d)、(e)和(f)可以发现经典IHS变换所得结果(d)较原TM影像扭曲严重，光谱失真较大，而曲波变换的引入则能有效克服传统IHS变换方式的不足，实现对原TM影像光谱信息更好的保留。对比道路、居民区等地物的光谱信息可以发现影像(f)在保持原TM影像的光谱信息以及影像清晰度方面要优于影像(e)，这说明文中方法的融合规则较经典曲波变换融合规则更加有效。

SAR数据与TM影像融合效果的定量化评价应综合考虑对光谱信息的保持和空间细节信息的增强等两方面。其中，反映光谱信息的指标包括扭曲程度和相关系数等，而反映空间细节信息的指标则可以用信息熵和清晰度等进行衡量。

扭曲程度直接反映了融合影像相比于原影像的光谱失真程度。扭曲程度值越小，融合影像的光谱失真程度就越少。其定义为：

(5)

式中DNF(i，j)与DN(i，j)分别为融合后和原始影像上(i，j)处的灰度值。

相关系数揭示了原始影像和融合影像之间在小尺寸结构上的相似性，反映了原始影像和融合影像的相关程度。相关系数越大，融合后的影像与原始影像的相关性越好。其定义为：

(6)

式中M与N分别代表影像的行与列，D(X)、D(Y)分别为影像X和Y的方差。

信息熵是衡量影像信息丰富程度的指标，是将影像内频率的变化作为信息。融合影像的信息熵越大，说明融合影像的信息量增加的越多。其定义为：

(7)

式中pi为影像X中像素灰度值为i的概率。

清晰度又称为平均梯度，是描述影像清晰程度的指标，反映了影像中微小细节反差与纹理变化特征，清晰度越大，影像越清晰，边缘信息保留越多。

(8)

原SAR数据、原TM影像以及各方法融合结果相应指标的统计结果如表1所示。从信息熵指标上看，融合结果的信息量较原SAR数据和TM影像都有增加，这与目视判读效果相一致，而方法间无明显差别。从扭曲程度和相关系数两项指标上看，融合结果较原TM影像均产生一定程度的扭曲和失相关。但在程度上文中方法所产生的扭曲和失相关均低于另外3种方法，这说明文中方法在对TM影像光谱信息的保持方面性能更优。从清晰度指标上看，融合结果的清晰度均优于原始数据。在基于全极化数据产生的融合结果中，文中方法较另外两种融合方法所得融合结果的清晰度更高，这说明文中方法能更有效地在融合结果中注入全极化SAR数据的细节信息。所不足的是基于全极化数据的融合结果在影像清晰度方面要低于单极化SAR数据的融合结果，这主要因为全极化数据I分量的产生实质上是HH、HV和VV极化数据间的线性加权，而这种方式在一定程度上会降低原影像的对比度，从而导致融合结果影像清晰度的部分损失。从各方法执行时所耗费的时间看，IHS变换经典融合方案耗费的时间最少，原因在于没有对数据进行曲波变换；文中方法所耗费的时间要略多于另外两种基于曲波变换的融合方法，原因在于文中方法需要进一步计算低频子代系数的梯度等信息。

表1 影像融合结果的定量评价

5 结束语

本文在IHS空间提出了一种基于二代Curvelet变换融合全极化SAR数据与TM影像的算法。针对全极化SAR数据与TM影像的特征，对Curvelet变换后的低频子带系数依据SAR数据的梯度信息与亮度信息确定其融合权值，对高频子带系数采用取绝对值大的系数进行融合。实验结果表明，该方法既可以很好地保留TM影像中的光谱信息，又能有效地注入全极化SAR数据中的细节信息。在目视效果与定量化指标上均优于基于经典IHS变换和经典Curvelet变换的融合算法。但方法在综合全极化SAR数据中的信息时存在损失影像清晰度的不足。另外，文中仅将全极化SAR数据当作彩色图像使用，而全极化数据是可通过极化变换等方式进一步挖掘信息的，如何提取出全极化SAR数据中最有益于参与融合的I分量也是后续需要进一步研究的问题。

参考文献：

[1] 李晖晖，郭雷，刘航.基于互补信息特征的SAR与可见光图像融合研究[J].计算机科学，2006，33(04)：221-224.

[2] HELLWICH O，GUNZL M.Landuse classification by fusion of optical and multitemporal SAR imagery[C].Proceedings of the International Geoscience and Remote Sensing Symposium，Honolulu，Hawaii，July，2000：1921-1923.

[3] 黄登山，杨敏华，姚学恒，等.基于Contourlet与IHS变换的SAR与多光谱影像融合方法[J].宇航学报，2011，32(1)：187-192.

[4] CANDES E J，DONOHO D L.Curvelets：A surprisingly effective nonadaptive representation for objects with edges[C].Proceedings of International Conference on Curves and Surfaces，4th，Saint-Malo，France，July，1999：105-120.

[5] 许星，李映，孙瑾秋，等.基于Curvelet变换的SAR与TM图像融合研究[J].西北工业大学学报，2008，26(03)：395-398.

[6] ZENG Y，ZHANG J X，VAN GENDEREN J L，et al.Image fusion for land cover change detection[J].International Journal of Image and Data Fusion，2010，1(2)：193-215.

[7] XU M S，XIA Z S，ZHANG F L，et al.Multi-temporal polarimetric SAR and optical data fusion for land cover mapping in southwest China[C].International Conference on Multimedia Technology，2010：1-4.

[8] 张微，丁献文，陈汉林.基于改进的IHS方法的SAR图像增强与融合[J].海洋学研究，2007，25(01)：73-79.

[9] GUNGOR O，SHAN J.An optimal fusion approach for optical and SAR images[C].Proceedings of the ISPRS Commission VII Mid-term Symposium “Remote Sensing：From Pixels to Processes”，Enschede，The Netherlands，May，2006：111-116.

[10] CANDES E J，DONOHO D L.New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities[J].Communications on Pure and Applied Mathematics，2003，57(2)：219-266.

[11] CANDES E J.Ridgelets：Theory and applications[D].Stanford University，1998.

[12] CANDES E J，DEMANET D L，YING L.Fast discrete curvelet transforms[J].Multiscale Model & Simulation，2006，5(3)：861-899.

[13] SVEINSSON J R，BENEDIKTSSON J A.Combined wavelet and curvelet denoising of SAR images using TV segmentation[C].Proceedings of the International Geoscience and Remote Sensing Symposium，Barcelona，Spain，2007：503-506.