多/高光谱遥感数据的类立体纹理特征

赵 巍，崔淑梅，吴 锐，刘家锋

(1.哈尔滨工业大学计算机科学与技术学院，150001哈尔滨;2.哈尔滨工业大学电气工程及自动化学院，150001哈尔滨)

源于灰度共生矩阵(Gray Level Co-occurrence Matrix，GLCM)的纹理特征30多年前已被证明其在分析和辨识遥感影像多种不同地物方面具有有效性［1-2］.近年来研究者［3］发现对卫星图像分类过程中采用纹理特征是克服光谱分析缺失的有效措施.如今遥感图像研究者［4-7］正以各种方式使用并改善着纹理这一经典的、可塑性高的遥感数据基元描述工具.

目前针对遥感数据的处理通常是对原始图像提取边缘、纹理等基本几何元素或特征.多/高光谱遥感数据的同源同点多幅图像同时获取的性质使其具有非常丰富的协同数据可利用;同时，分析遥感图像数据的有效特征可能就是获取图像时所设定的多个光谱的波长或波长范围［8］，因此研究者［9］对纹理等图像特征的提取很早就开始向多维化、立体化的方向发展.但目前的解译方法仍缺乏充分挖掘影像各种空间特征和地物其他属性特征等方面的探索［10-11］.文献［9］在进行纹理变换和分析时采用通道信息组合的方式计算3×3像素窗的中心像素点的局部属性，认为多通道信息融合对边缘检测和图像增强的效果要好于在地物主题映射的效果.Mounir等［12］从16个子波段、3个通道中分别抽取熵、纹理均值、3阶矩、4阶矩和能量5个参数组成特征向量.这种纹理特征构造方法可以增加信息量并提高分类的精度［13］.此外，采用均值法和主成分分析法，将子图像各个波段基于GLCM的纹理统计特征向量作为多波段遥感图像的特征向量［14］，以挖掘同源不同尺度信息，这样虽增加了纹理分析的复杂性，但分类性能均有不同程度的增益.上述同源遥感数据的纹理研究中较少有直接考虑波段之间像素灰度变化的关联关系的做法，多是分别提取单波段数据的纹理特征，再在特征层或决策层实现各波段信息融合的方式.比较而言，在多源遥感数据分析中的像素级融合的情况则非常普遍［15-17］.

由于GLCM是对单幅图像遍历得到的统计矩阵，这可能限制其在多/高光谱遥感解译分析能力.本文根据多光谱遥感数据的同源同点多幅同时获取的特点，提出了一种多光谱遥感图像底层数据综合方法，即灰度级差关联概率矩阵(Gray Level Difference Associated Possibility Matrix，GLDAP)描述.目前，能够采集立体信息的成像方式大致有两种:双目成像和全息成像.对于两幅同源不同波段的遥感图像来说，显然具有比单幅图像更丰富的视觉信息，当将其中一幅作为参考时，其与另一幅图像的差异即被显著化，灰度级差关联概率矩阵GLDAP以统计的方式记录了这些差异.GLDAP方法是对两幅同源图像的底层数据进行视觉差异统计，是对立体成像过程的模拟.模拟过程中所获得的信息被称为类立体信息，在类立体数据的基础上得到的纹理特征即为类立体纹理.

1 灰度级差关联概率矩阵

事实上多/高光谱遥感数据除了二维平面像素点的值之外，没有提供第3维度量信息，并非真正意义上的立体数据.每个波段图像可看作是透过不同透镜观察到的景物，只是无法可靠地指出哪个波段的图像更能准确、清晰地给出图像内容的全部细节信息，只是依赖于单波段图像提供的局部可靠信息来解译图像.人们已经根据地物的波谱特性获得了对特定类型地物有较好可分性的波段范围，从而有针对性地获取特定地物在相应单波段图像上的信息.

本文依据地物谱特征分析，参考区分特定地物的波段信息，从遥感数据中选择多幅有效波段图像，统计每两幅图像的灰度的协同变化规律，进行视觉差异统计分析，通过计算差异来刻画和提取有效纹理特征，以克服GLCM方法对单幅图像统计描述的局限性，实现数据解译分析的高精度目标.进行同源多幅图像处理时，对每次选定的两幅图像进行像素信息的协同统计，建立两图像在同一地理位置上对应像素点之间灰度值变化的关联关系，计算像素点对的灰度差值改变的概率，构成GLDAP概率矩阵;基于此矩阵计算统计纹理特征，如能量、对比度及熵等.视觉差异分析方法以类似于立体视觉的方式来处理多/高光谱图像.

Haralick等［18］提出经典的GLCM纹理分析方法设置窗口遍历一幅图像，统计差异并记录在灰度共生矩阵中，在灰度共生矩阵上提取包括能量、熵在内的14种纹理特征.一幅图像在0°，45°，90°和135°等4个方向上的灰度共生矩阵的计算为

式中:P为灰度共生矩阵;#为集合中元素的数目;(k，l)、(m，n)分别为像素点坐标;i，j分别为灰度级;d为两像素之间的距离.

本文在GLCM方法的基础上提出的GLDAP纹理分析方法，统计同源两目标波段图像F和F'中在相同方向、相同像素点距离d上灰度级差分别为Δi与Δi'的像素对的数目;以此来反映两图像在相同方向和像素距离上的灰度变化相关性.GLDAP矩阵的定义为

式中:~P为相应方向的灰度级差关联概率矩阵;#为集合中元素数目;(k，l)、(m，n)和(k'，l')、(m'，n')分别为F和F'中的相应两点坐标;Δi、Δi'分别为F和F'中相应两点灰度差值;d为窗口内两像素间距.

按照GLDAP的定义，应当选定方向和点间距d，根据式(1)同时遍历两个波段的图像并归一化便可获得GLDAP矩阵.

GLDAP矩阵是非对称的，其行列值分别是两幅图像的灰度差的绝对值，很好地反映了在两幅同源图像的相同距离和相同方向上像素灰度差的相关变化的统计信息.当两幅图像灰度分布一致时，主对角线和为1，各非主对角线和为0;随着灰度分布差异增大和纹理复杂程度增加，主对角线和将随之减小，各非主对角线和相应增加.因此，对角线和也是一种重要的统计特征.GLDAP矩阵主对角线元素表示两幅图的相邻像素灰度差相同的概率、各非主对角线元素表示两幅图的相邻像素灰度差不同的概率，计算公式为

式中:GLDAP矩阵大小为K×K;i，j分别为矩阵的行、列数;t为所求取的对角线;t=0为主对角线.

2 实验结果与分析

用GLDAP进行多光谱遥感图像中的水体、城区、山区3种类型的地物进行分类，并同GLCM的分类结果进行比较.为使实验具有可比性，实验过程中为GLCM矩阵方法提供与GLDAP方法相同的两幅单波段图像.设提取D种纹理特征，则GLCM方法采用文献［12］的常规处理方式，分别获得两幅图的纹理特征再构成特征向量，维数为2D;GLDAP方法除了这些纹理特征外还使用了主对角线和特征来反映两波段具有相同灰度差的统计量，其特征维数为D+1.

设图像的大小为N×N，预设窗口尺寸、灰度级及特征数等参数为常量，两种方法的时间复杂度均为O(N2).

实验系统的决策模块采用支持向量机(Support Vector Machine，SVM)方法，选择径向基函数(Radical Basis Function，RBF)作为核函数.软件在装有VC、Open CV环境的PC机上实现.

2.1 样本与前期实验准备

采用Landsat-7的增强型专题制图仪ETM+多光谱数据.第1～5及7波段的图像如图1所示，其中B为波段;图中不包括波段6(热红外)数据.单波段图像灰度级为256，为减少计算量和避免信息损失，将灰度级压缩为64级.借助ERDAS Imagine遥感图像处理系统为ETM+多光谱数据准备训练和测试样本.

图1 ETM＋多光谱波段数据图

水、陆两类分类的预备性实验中，进行了第1～5、7波段两两组合以及与GLDAP矩阵的不同计算窗口尺寸的组合测试;d取值为1，即只比较相邻两像素点的灰度变化;选取能量、对比度、均一性和相关性4个特征;选择了近1 500个样本点对支持向量机SVM进行训练，对近350 km2的水陆共存的遥感图像进行分类实验测试.两类分类的实验结果显示，GLDAP有着与GLCM相当的分类能力，将GLDAP应用于遥感图像的分类是完全可行的;选用第4、7波段组合和第1、7波段组合时，与5×5大小的窗口配合使用，地物边界效果理想.

基于上述实验参数条件和实验平台，用6个波段的10 973个样本对SVM分类器进行训练，截选若干区域进行了水体、城区和山区3类地物样本的分类测试.选择两幅有以上3种地物的区域作为本文实验数据处理的演示图像，其合成ETM+数据的第3、2、1波段组合为模拟真彩色图像如图2，3所示.肉眼观测该图可分辨城区纹理细密，山区纹理相对粗糙，水体灰度分布均匀、纹理单一.

图2 ETM＋多光谱波段数据第3、2、1波段合成示意图

图3 ETM＋多光谱波段数据第3、2、1波段合成示意图

2.2 分类实验结果样例与分析

对遥感图像进行水体、城区和山区的3类分类，选用第4、7波段组合和第1、7波段组合，配合5×5的滑动窗口.分类结果用黑、白、灰3色标记:将水体标记为白色，城区(人工迹地)为黑色，山区为灰色.

2.2.1 第4、7波段图像分类的结果

对图2的第4、7波段数据处理结果如图4所示.图4(a)是用GLCM方法的类别标记结果，能够大略看出山区的轮廓，但误识多.如:将部分山区地貌误识为城区，部分区域因纹理较单一被误识为水体;城区部分绿地与山区纹理很相近而被误识为山区;水体中间部分河床较高而水体很浅，被误识为山区.

图4(b)是用GLDAP方法对山区附近图像的分类结果，可以清晰地看出山区轮廓，山沟之间的居民区能正确地标记为城区，水体中间的浅滩部分也能正确识别.城区的绿地等部分由于纹理较粗糙，也存在误分类情况，但明显少于GLCM方法的误分量.这两种方法对水体部分的识别性能差异不大.

图4 图2的第4、7波段采用GLCM方法与GLDAP方法分类处理结果

对图3的第4、7波段数据处理结果如图5所示.图5(a)为GLCM方法的分类结果，山区和城区两类地物都互相存在大量的误分类，几乎无法分辨;水体两岸较为接近的区域误分成山区，而右下角灰度值较单一的山区部分则被误分为水体.使用GLCM的误识情况更为严重.

图5(b)为GLDAP分类结果，从图中可明显看到，山区和城区能够很清晰地分辨，且分类准确性很高.

图5 图3的第4、7波段采用GLCM方法与GLDAP方法分类处理结果

2.2.2 第1、7波段图像分类的结果

采用第1、7波段图像分类的结果分别如图6、7所示.从图6、7可以看出，使用GLCM纹理分类时水体和山区的区分度很差，水体的大部分被误分成了山区，而山区的大部分识别成了水体.造成山区和水体互有误识的情况主要是:1)水体的反射主要在第1、2波段，由于水体的反射率较高，第1波段图像水体的颜色较浅，而且受水体的深浅和悬浮物质的影响第1波段的纹理较第4或第7波段的复杂;水体对其他的波段吸收都很强，特别是第4波段，因此第4波段的水体图像呈黑色，且纹理单一;2)山区地表主要是绿色植被，绿色植被在蓝光区域(第1波段)有一个吸收谷，因此山区呈现的颜色也较深，纹理也较单一，与水体较接近;导致互相误识的情况大量发生.而由于GLDAP方法统计了两个波段间的灰度变化的关系，因此受单波段图像的影响较小，第1、7波段组合3种类型的区域基本上能够正确识别，山区的正确识别率较GLCM方法有大幅提高.

可见，选用GLDAP处理多波段图像时能够避免因单波段中某两类地物的可分性差而造成的误识.

图6 图2的第1、7波段采用GLCM方法与GLDAP方法分类处理结果

图7 图3的第1、7波段采用GLCM方法与GLDAP方法分类处理结果

2.3 统计数据及性能分析

两种方法分类效果的定量比较及平均分类正确率统计如表1所示.平均分类正确率是从样本图像中随机选取地物的若干典型样本区，计算每类地物正确分类像素数与实际该类像素数的百分比，将各类正确分类率取平均得到的.

表1 地物识别平均正确率%

实验表明，采用第1、7波段组合和第4、7波段组合，GLDAP对多光谱遥感图像进行分类比传统的GLCM方法更准确;第4、7波段组合的GLCM和GLDAP纹理方法优于第1、7波段组合的结果.

实验还分别统计了3种地物互有误识时所对应的6种平均误识率(如水体误识为城区、城区误识为山区的比率等).表1中第1、7波段组合实验GLCM方法的水体(14.4%)、山区(13.3%)的识别率极低，而城区识别率却高达98.3%.表2给出了两种方法的被误识为城区的平均统计量，从中可看出，GLCM将水体和山区误识为城区的比率之和高达80.8%，而GLDAP方法在误识为城区的比率则较该方法其他类型的误识差异不大.

表2 第1、7波段组合误识为城区的统计量比较%

同时还看到，与第4、7波段组合相比，第1、7波段组合时水体被误识得最为严重.这与特定波段反映地物谱信息的能力不同直接相关，导致了所构建的波段组合形成的视觉差别和基于此视差获得的统计信息在反映地物能力上的强弱有别.在对多/高光谱数据进行研究时，配合波段相关性分析，采用GLDAP方法充分测试不同波段组合，将对高效、完整获取类立体数据的可利用信息起到至关重要的作用.

3 结论

1)本文从各波段数据的内在关联出发，根据地物的波谱特征，选择双波段统计灰度协同变化规律，建立视觉差异统计方法，基于GLDAP矩阵提取类立体纹理特征，以克服GLCM方法对单波段图像统计描述的局限性;用类立体视觉方式处理多/高光谱图像，以提高遥感数据解译分析精度.方法同样可用于高光谱数据分析.

2)实验表明，GLDAP纹理提取方法具有良好的分类性能，3种地物分类效果明显优于GLCM方法，并能在一定程度上避免GLCM因单波段中两类地物可分性差而导致误识的情况.两种方法的时间复杂度相同，为多/高光谱遥感数据解译人员在相同时间条件下获得更优良的性能提供了新的可操作性强的分析途径.

［1］WESZKA J S，DYER C R，ROSENFELD A.A comparative study of texture measures for terrain classification［J］.Comparative and General Pharmaclology，1976，1(4):269-285.

［2］SHANMUGAN K S，NARAYANAN V，FROST V S，et al.Textural features for dadar image analysis［J］.IEEE Transactions on Geoscience and Remote Sensing，1981，19(3):153-156.

［3］CARLEER A P，WOLFF E.Urban land cover multilevel region-based classification of VHR data by selecting relevant features［J］.International Journal of Remote Sensing，2006，27(6):1035-1051.

［4］OHANIAN P P，DUBES R C.Performance evaluation for four class of texture features［J］.Pattern Recognition，1992，25(8):819-833.

［5］NYOUNGUI A N，TONYE E，AKONO A.Evaluation of speckle filtering and texture analysis methods for land cover classification from SAR images［J］.International Journal of Remote Sensing，2002，23(9):1895-1925.

［6］RAJESH K，JAWAHAR C V，AENGUPTA S，et al.Performance analysis of textural features for characterization and classification of SAR images［J］.International Journal of Remote Sensing，2001，22(8):1555-1569.

［7］舒宁.关于多光谱和高光谱影像的纹理问题［J］.武汉大学学报(信息科学版)，2004，29(4):292-295.

［8］RICHARDS J A，JIA X P.Remote sensing digital image analysis，an introduction［M］.Berlin，Heidelberg:Springer-Verlag，2006.

［9］JAMES R I，GARY W P.Texture transforms of remote sensing data［J］.Remote Sensing of Environment，1981，11:359-370.

［10］舒宁.关于遥感影像处理分析的理论与方法之若干问题［J］.武汉大学学报(信息科学版)，2007，32(11):1007-1015.

［11］骆剑承.遥感影像智能图解及其地学认知问题探索［J］.地理科学进展，2000，19(4):289-296.

［12］MOUNIR A K，AHMED H，DRISS A.Textural feature selection by joint mutual information based on gaussian mixture model for multispectral image classification［J］.Pattern Recognition Letters，2010，31(10):1168-1174.

［13］TANG X，STEWART W K.Optical and sonar image classification:wavelet packet transform vs fourier transform［J］.Computer Vision and Image Understanding，2000，79(1):25-46.

［14］田艳琴，郭平，卢汉清.基于灰度共生矩阵的多波段遥感图像纹理特征的提取［J］.计算机科学，2004，31(12):162-163.

［15］卢焱.多源遥感影像融合方法研究［J］.吉林大学学报(地球科学版)，2007，37(11):176-179.

［16］贾永红，李德仁.多源遥感影像像素级融合分类与决策级分类融合法的研究［J］.武汉大学学报(信息科学版)，2001，26(5):430-434.

［17］许君一，卢秀山，卿熙宏，等.基于热传导模型的像素级遥感图像融合［J］.测绘学报，2010，39(6):566-571.

［18］HARALICK R M，SHANMUGAM K，DINSTEIN I.Textural features for image classification［J］.IEEE Transactions on System，Man and Cybernetics，1973，SMC(3):610-621.

［19］BARALDI A，PARMIGGIANI F.An investigation of the textural characteristics associated with gray level Cooccurrence matrix statistical parameters［J］.IEEE Transaction on Geosciences and Remote Sensing，1995，33(2):293-304.