卜丽静,吴 畅,何志博,张正鹏
(辽宁工程技术大学 测绘与地理科学学院,辽宁 阜新123000)
合成孔径雷达(Synthetic Aperture Radar,SAR)是一种具有全天候作业、识别伪装和穿透掩盖物等优点的主动式微波对地观测系统,能够提供高分辨率的SAR图像。光学图像具有地物光谱信息丰富的特点,但云雾雨雪天气成像困难。遥感图像融合,能够将单一或多个传感器的数据和信息进行关联、相关和合成,得到更精确的位置和特征估计,提高遥感数据在地物提取中的识别能力,扩大应用范围[1]。光学图像之间融合研究较多[2-3],近几年,SAR与光学图像的融合也逐渐成为研究热点,如,经典的 PCA,I HS[4],BT,HPF,WT 等融合方法[5-6],小 波 融 合 方 法[7-11];非 下 采 样 Contourlet变换和HIS变换相结合融合算法[12]等。这些方法具有各自的特点,如PCA方法能很好的保持多光谱信息,但SAR图像纹理结构信息保持相对较少;HPF和WT、小波、Contourlet变换方法纹理结构信息保持较好,但也存在损失SAR图像高频信息或空间细节信息的情况 I HS和BT方法视觉效果较好,色彩及光谱失真现象严重[6]。这些方法主要关注纹理或光谱特征,但存在某一特征保留较好,另一特征有不同程度损失的情况。而且,从SAR图像的散射特性角度考虑纹理特征较少,噪声抑制效果也较差。针对以上问题,本文提出基于纹理特征的SAR与光学图像融合方法。分析SAR的散射特性和图像特点,在考虑光谱特征和SAR的纹理特征基础上,达到尽量多的融合两种图像优点的目的。
SAR图像与光学图像比较,具有以下特性:①存在相干斑噪声,地物细节信息不丰富,灰度相似的地物间区分困难,但纹理特征明显。SAR图像中含有丰富的纹理,不同的地表粗糙度呈现出不同的纹理特征;纹理对图像区域和表面的感知、描述的独特作用可以区分某些目标,作为一种有力的边缘检测手段[13],可以提高分类精度[14]。它是雷达识别地物的一项关键技术[15];②图像中金属物体和建筑物等具有强散特性的地物回波较强,在图像中较亮,特征明显。由于建筑物后向散射由特殊反射与漫反射的后向散射组成,前者一般要强得多,建筑物的特殊反射可以分为来自倾斜屋顶的单反射、来自墙-地结构的二次反射、来自墙-墙-地结构的三次反射三种机制,一般来说反射次数越多,后向散射回波越弱,对于建筑物来说前两种反射机制占主导地位,因此,建筑物后向散射强,在图像中较亮[16]。光学图像中没有这些特性。因此,可以充分利用SAR图像中的这些信息,与光学图像融合,达到改善图像质量,增加信息的目的。
高分辨率SAR图像的纹理特性和强散射特性优于光学图像。因此,融合过程中尽量保留利用这些优势特性。首先,利用灰度共生矩阵方法,提取SAR 的 均 值 (Mean)、方 差 (Variance)、协 同 性(Ho mogeneity)、对比 度(Contrast)、相异 性(Dissi milarity)、信息熵(Entropy)、二阶矩(Second Moment)和相关性(Correlation)8种常用纹理统计量。这些纹理特征值经归一化后能够组合成类似于多光谱图像的多波段图像。由于各纹理特征之间具有一定的相关性,所以当纹理特征的数目增加时,图像的复杂性也随之增加,并不一定得到良好的组合效果[17]。因此,通过分析各个纹理特征之间的相关性矩阵,选择可分离性和重要的纹理,并要求强散射特征也被很好保留。通过上述分析,本文选择均值、方差、相异性、二阶矩四种纹理生成主要纹理特征图,该图像由式(1)计算。
式中:(i,j)表示图像行列值;SARTexture表示求得的综合纹 理 特 征 图 像;SARMean、SARVariance、SARDissimilarity、SARSecondMoment表示提取出来的均值、方差、相异性、二阶矩四种纹理特征;a1,a2,a3,a4表示计算重要纹理特征图时的阈值。Otsu[17]方法是一种全局化的灰度图像动态二值化方法。该算法的基本思想是:设使用某一个阈值,将灰度图像根据灰度值大小,分成目标部分和背景部分两类,在这两类的类内方差最小和类间方差最大的时候,得到的阈值是最优的二值化阈值。因此,该方法能够按图像的灰度特性将图像分成背景和目标两个部分。本文在处理各个纹理特征图像时,需要保留的就是纹理特征中的主要信息,即目标信息,该方法得到的阈值能够很好的分离出这一信息,因此采用该方法计算式(1)中的阈值。计算过程中,对四种纹理特征像素值与相应阈值逐一判断,只要满足式(1)中的情况,就保留这一位置对应的SAR图像像素值,否则这一位置值为0。通过上述计算,既能够保留SAR图像中主要纹理特征信息,又可以去除相对不重要的特征信息和噪声,为下一步的融合提供最佳信息。
与高分辨率SAR图像相比,低分辨率多光谱图像的色调、饱和度信息丰富,但V分量分辨率低、纹理信息较少。因此,可以充分利用光学的光谱信息和SAR的高分辨率信息改善融合效果。在HSV(Hue Saturation Value)彩色空间可以将色调、饱和度、亮度信息分离[17]。因此,对光学图像进行HSV变换,分离出H,S,V分量。通过改善V分量可以达到改善图像质量的效果。由于Contourlet变换作为一种“真正”的图像二维表示方法,能捕捉到图像信息内在的几何结构特征[8],具有高度方向性和各向各异性,对图像纹理细节信息表达具有很好的效果。因此,利用Contourlet多尺度变换,对SAR和光学图像的高低频信息进行融合是目前研究的重要方法之一。但以往方法中,存在多光谱图像色彩信息保留好,SAR图像信息处理单一的问题,如全部保留低频分量,舍弃高频分量[8];由于SAR的低频和高频分量中均含有纹理和强散射特征,这样做容易造成SAR信息损失。因此,综合考虑SAR和多光谱图像中的优势信息,本文提出改进Contourlet融合方法计算新的低频和高频信息方法,用主要信息占优原则计算融合后的低频和高频信息。新低频信息计算时,首先将光学和SAR图像的每一个像素值和先验阈值进行判断,如果纹理特征图像的灰度值小于阈值,新低频计算完全采用光学图像,如果大于阈值,则采用改进的加权平均法,算式如式(2)所示[12]。此处改进的目的是达到控制融合过程中SAR和光学信息的目的,融合时尽量保留最占优信息。
式中:Ca,j和Cb,j分别表示SAR和光学图像分解后的低频分量系数;Cc,j表示融合图像的低频分量系数;a,b为加权因子;(Ca,j+Cb,j)×a表示两幅图像的加权均值,对图像的亮度起着决定性的作用;(Ca,j-Cb,j)×b表示两幅图像的 加 权 差值的绝对值,包含两幅图像的边缘信息。新高频信息的计算采用绝对值较大法,如式(3)所示。
式中:i=1,2,…;DiA,JDiB,JDiF,J分别为各分解层上对应非下采样Contourlet变换系数[12]。
最后,将计算出最优的高频信息与低频信息进行Contourlet逆变换,得到新的强度分量V′。用原光学图像的H,S和新的亮度分量V′进行HSV逆变换,得到融合后的图像。算法流程如图1所示。
图1 融合方法流程
实验采用Cos mo-Sky Med高分辨率SAR图像和Landsat8光学图像,图像参数如表1所示。实验区域中包括公路、铁路、建筑物、植被、金属管线、水池等地物。SAR图像中铁路、管线、建筑物等地物具有强散射特征,强度值较高,复垦区、植被区、建筑物等的纹理特征明显。相比之下,光学图像分辨率低,这些信息稍弱,但光谱信息丰富。原始SAR和光学图像如图2所示。
实验过程是首先对SAR图像用Frost滤波去噪,为兼顾图像中地物的平滑效果,选择3×3的滤波窗口进行噪声抑制。然后,用灰度共生矩阵方法计算纹理特征,选择7×7窗口,64级的灰度量化级。分析纹理特征间相关性,选择均值、方差、相异性、二阶矩四种纹理,并按式(1)计算出重要纹理特征图像。实验中阈值的设置原则是尽量保持SAR图像中的重要纹理特征,阈值选择过大,会导致保留的纹理特征过少,下一步融合用到的纹理特征就会减少,导致融合结果中SAR的纹理特征少;但阈值设置过小又会导致保留的特征过多,主要特征不突出,影响后期Contourlet变换高低频信息选择结果,使最后的融合结果与直接Contourlet变换相比改善不大。实验中设置的阈值分别为a1=0.2,a2=0.1,a3=0.4,a4=0.1,得到主要纹理特征图像如图3(a)所示。并将光学图像进行HSV变换,亮度分量V如图3(b)所示,V中包含光学图像的灰度信息。对比图3中的(a)和(b),铁路、建筑物、复垦区、金属管线等地物在光学图像中表现不明显或边缘信息对比度差,SAR重要纹理图像中这些信息丰富。然后,对SAR重要纹理和光学V分量进行Contourlet分解。利用分解后的高低频信息计算新低频和新高频信息,为体现本文的主要信息占优原则对结果的影响,本文方法做了两组实验。实验1中,归一化SAR和光学的低频信息图,阈值设为0.7,综合利用SAR和光学的信息,按式(2)计算,新低频信息计算的权值a=0.7,b=0.3,高频计算使用式(3)。实验2中,以SAR信息占优为主,阈值设为0.3。最后,将H,S,V′分量进行 HSV逆变换,转换到RGB彩色空间,得到本文算法融合结果,如图3中(g h 所示。为了与其它方法进行对比,将滤波去噪后的SAR图像与光学图像进行小波融合、HSV融合、Br ovey融合、直接Contourlet变换融合。其中,小波融合中融合方法使用I HS方法。结果如图3中的(c)、(d)、(e)所示。本文将处理前后图像均归一化到0~255之间,选取均值、均方根误差(RMSE)、标准差、平均梯度、边缘保持指数(EKI)5个指标对融合后的图像做定量评价,评价结果见表2。
表1 SAR与光学图像数据参数
图2
均值表达图像平均亮度值。如果亮度值适中,则视觉效果良好。其值越大图像越亮,其值越小图像越暗。RMSE和标准差的计算见式(4)、式(5)。
式中:M,N 表示图像的长和宽;f(xi,yj)表示标准图像像素值;f′(xi,yj)表示处理后图像像素值;σ表示图像灰度的标准差;μ表示处理后图像灰度的统计均值。两幅图像的RMSE是以一幅图像为标准,另一幅图像与它之间的差值来计算。它反映处理后图像像元点偏离标准图像的程度,防止图像偏差过大失去原有信息过多,值越小说明处理效果越好。本文计算时以光学为标准图像。从表2中可以看出本文方法2的RMSE最小,说明融合后图像保留了光学图像的主要信息,与光学图像偏差不大。标准差表示处理后图像各个像素点处的灰度与灰度均值之间的离散情况。如果该值较大,图像灰度级的分布越为离散,图像的反差大,从图像上能够看出更多的信息。其值较小图像反差则较小,对比度不大,色调单一均匀,看不出太多的信息。本文方法的标准差最大,说明融合后图像能够看出更多信息。平均梯度能够很敏感地反映图像对微小细节反差表达的能力,可用来评判图像的清晰程度,同时还反映出图像中微小细节反差和纹理变换特征。平均梯度值越大,说明图像越清晰。EKI表示图像边缘的保持程度,在原图像上EKI值为1,EKI小于1表示边缘被模糊,EKI大于1表示边缘被增强。平均梯度和EKI值两个指标中Contourlet和本文方法2较好,这是因为它们保留了较多的SAR的信息。因此,从评价指标来看本文方法最好。
图3 实验结果
表2 融合效果评价指标
从主观评价角度对比实验结果,HSV和Br ovey融合后的图像基本上保留所有SAR图像的信息,但多光谱图像的纹理信息基本没有体现,仅保留一些色彩信息,整体效果类似与SAR图像上简单赋予光学的颜色,而且颜色信息与光学图像中的光谱信息有差异,没有达到综合两种数据源优势的目的。小波变换融合则是更多的保留了光学的纹理和色彩,对SAR的信息保留非常少,图像仍模糊不清,没有达到提高图像分辨率、增加地物细节的目的。对比原光学图像(图2中(a))、SAR图像(图2中(b))和本文融合结果图3中(f),可以看出,图2中(a)的1和2区域,实地地物为输煤铁路,SAR图像中信息明显,本文方法信息被很好地保留,3,4,5的复垦区阶梯纹理也都被保留,6和7对应的建筑物区域,强散射信息被保留,而且光学的色彩信息也被很好的保留。对比图3中的(f)、(g)、(h),可以看出上述的几个对比区域均比原图像和其它融合方法好,但三者之间又有区别。直接Contourlet融合的图像由于没有主要特征判断,导致两幅图像的信息都有,但结果并不一定有利于解译识别。如6区域的房屋信息,光学的信息更好些,图3中的(g)结果更好些。4,5区域中,SAR的信息更好些,边缘特征体现明显,图3中的(h)最好。因此,总的来说,从定量指标分析和融合后图像可以看出,本文融合方法能够根据关注的主要特征,更好地综合两种数据的优势,增加地物细节信息,提高图像的分辨率。
针对SAR与光学图像的融合问题,在分析SAR与光学图像特点基础上,提出一种基于SAR图像中纹理和强散射信息的融合方法。通过对纹理特征相关性分析得到重要纹理特征,并利用改进Contourlet变换融合方法综合光学和SAR的高低频信息,进而融合两种图像的优势。并将本文方法与小波、HSV、Br ovey融合方法对比分析,实验表明本文方法能够较好的保持高分辨率SAR图像强散射、纹理细节信息和光学图像的光谱特征,具有较好的融合效果,能够提高图像的解译能力。总的来说,SAR和光学融合是一个非常复杂的问题,本文方法还有一些问题,如实验中3 m SAR和30 m的光学图像分辨率差异有些大,如果融合数据源分辨率相差小些效果会进一步改善。而且,融合图像实质上是光学和SAR哪个信息占优的问题,如何能够自适应的选择需要保留的信息,也需进一步研究。同时,Contourlet分解的级数和最佳权重中的确定问题,低分辨率光学图像的光谱信息是否全部保留问题,还有待于更进一步的研究。
[1] 李军,林宗坚.基于特征的遥感影像数据融合方法[J].中国图像图形学报,1997(Z1):103-107.
[2] 姜芸,王军.多源遥感影像融合在哈大齐土地利用分类中的应用[J].测绘工程,2010,19(4):34-38.
[3] 朱蕊,邱茂,胡英男.面向空间数据更新的多源数据融合关键技术研究[J].测绘工程,2013,22(4):22-25.
[4] 贾永红,李德仁,刘继林.四种HIS变换用于SAR与T M 影像复合的比较[J].遥感学报,1997,2(2):103-106.
[5] 杨智翔,何秀凤,危来龙,等.一种多极化SAR图像融合方法研究[J].测绘科学,2014,39(2):14-18.
[6] 邹丽丽,崔海山,李颖,等.SAR与SPOT数据融合方法研究[J].遥感技术与应用,2010,25(6):836-841.
[7] 高文涛,汪小钦,凌飞龙.雷达与光学图像小波融合方法研究[J].地球信息科学,2007,9(4):129-132.
[8] HONG Gang,ZHANG Yun,MERCER Bryan.A wavelet and I HS integration method to f use high resolution SAR with moderate resolution multispectral i mages[J].Journal of the American Society for Photogrammetr y and Remote Sensing,2009,75(10):1213-1223.
[9] CHIBANI Y.Additive integration of SAR features into multispectral SPOT i mages by means of the a’trous wavelet decomposition[J].ISPRS Journal of Photogrammetr y and Remote Sensing,2006,60(5):306-314.
[10]ALPARONE L,BARONTI S,GARZELLI A,et al.Landsat ET M+and SAR i mage f usion based on generalized intensity modulation[J].IEEE Transactions on geoscience and remote sensing,2004,42(12):2832-2839.
[11]黄登山,杨敏华,姚学恒,等.基于a’trous小波与广义HIS变换的SAR与多光谱影像融合[J].遥感信息,2011(1):9-13.
[12]苏志渊,张蕾,普杰信,等.NSCT变换的SAR和可见光图像融合[J].计算机工程与应用,2011,47(5):166-168.
[13]CHEN Shu-peng,TONG Qing-xi,GUO Hua-dong.Study of Mechanism on Remote Sensing Infor mation[M].Beijing:Science Press,1998.(in Chinese).
[14]胡召玲,李海权,杜培军.SAR图像纹理特征提取与分类研究[J].中国矿业大学学报,2009,38(3):422-427.
[15]SHAO Yun,GUO Hua-dong,WEI Xiu-ping.Macroscopi-cal Texture Analysis and Geology Application Effect of Radar Image[A].Airborne Radar Remote Sensing and.
[16]郭华东.雷达对地观测理论与应用[M].北京:科学出版社,2000.
[17]R.C.冈萨雷斯,R.E.伍兹,著,数字图像处理[M].阮秋琦,阮宇智,译.北京:电子工业出版社,2007:229-235.