灰度共生矩阵纹理特征与马尾松林测树因子关系初探

李 进，王 鑫，张 园

（浙江农林大学环境科技学院，浙江 临安 311300）

近年来，高空间分辨率的遥感图像（卫片或航片）作为森林资源调查的有效手段，发挥着越来越重要的作用。由于空间分辨率很高（一般在5 m以内），因此包含了很多空间纹理方面的信息[1]。纹理是复杂的视觉实体或者子模式的组合，有亮度、色彩、陡度、大小等特征，因而纹理可以认为是在局部窗口内，影像灰度级之间的空间分布及空间相互关系[2]。灰度共生矩阵（空间灰度相关方法）是目前最常见和广泛应用的一种纹理统计分析方法，它通过对影像灰度级之间联合条件概率密度的计算表示纹理[3]。

国内结合灰度共生矩阵纹理特征进行森林测树因子建模的文章较少，李明诗、谭莹、潘洁等利用SPOT5 HRG影像进行纹理指标提取，结合研究区地形特征，分别对5个森林类型实现其生物量估算模型辨识与验证，结果发现，少数纹理特征ME（均值）、VA（方差）在森林生物量估算上是有效的并且是重要的，大多数纹理指标与生物量的关系不紧密；纹理指标对针叶林空间形态的表达能力优于阔叶林[4]。国外学者Christine分别通过建立平均胸径与原始图像灰度值、方差纹理指标的平均值以及第一主成分的图像灰度值的回归关系进行估计[5]，结果发现用7×7窗口的方差纹理图像拟合的平均胸径，其相关系数的平方达到了0.623；Ibrahim Ozdemir等则利用分辨率为15 m的ASTER图像提取纹理信息，建立了一致性纹理值与树冠大小的回归关系，其相关系数达到了 0.61[6]。本文初步探讨了森林资源二类调查资料中马尾松林的年龄、平均胸径、平均树高与郁闭度等测树因子与纹理指标之间的关系，为研究灰度共生矩阵纹理特征与其他测树因子的关系提供了参考依据。

1 方法

1.1 研究区概况

研究区为浙江省台州市仙居县，位于浙江省东南部，地处括苍山脉中段北麓，台州市的西部，东连临海、黄岩，南邻永嘉县，西接缙云县，北靠东阳市、磐安县和天台县，位于 120° 9′ ～ 121° 40′ E，28° ～ 28° 5′ N。该地区属浙东盆地低山区，海拔 700 ～ 1200 m[7]，气候常年温和湿润，森林覆盖率高达 77.2%。马尾松（Pinus massoniana）是我国分布面积最广的针叶树种之一[8]，同时也是仙居县的优势树种，面积7万hm2，占仙居乔木林总面积的62%[9]。

1.2 数据预处理

本文所采用的遥感数据是2001年1月到2003年11月期间拍摄的多幅仙居县单波段航片图，空间分辨率为1.27 m。图像采用了西安80坐标系进行配准。

本文所用的样地数据是仙居县森林资源二类调查资料，对其中的马尾松树种相关资料进行了筛选。为了准确提取马尾松的纹理信息，我们将马尾松的样地中心坐标进行了矢量化，以便于与航片图像进行连接。由于每块样地面积为800 m2，我们在图像上样地中心坐标所在的位置周围选取了与中心点纹理相似的位置进行感兴趣区域（AOI）截图，每个子图的大小为边长15 ～ 25个像素不等的正方形。需要注意的是，马尾松所在的样地的树种结构很多，为了便于研究马尾松的纹理特征，本文提取的样地树种结构大多为针叶纯林，少部分为针叶相对纯林。由于航片的阴影会造成一定的影响，本文在进行图像提取时，尽量避免山地阴影和航片制图过程中产生的阴影，选取了符合条件的15个马尾松纯林样地的子图进行纹理计算。

1.3 方法

1973年，Haralick首先提出灰度共生矩阵（GLCM）[3]，此后，它成为最常见和广泛应用的一种纹理统计分析方法。灰度共生矩阵（空间灰度相关方法）通过对影像灰度级之间联合条件概率密度P（i,j/d,θ）的计算表示纹理。通常θ方向为：0、45、90、135° 4个方向。这样P（i,j/d,θ）为一对称矩阵。如果d相对纹理的粗糙度小，共生矩阵的元素值将集结在对角线附近；反之，如果d较大，共生矩阵的元素值将离开主对角线向外散开分布。

Haralick定义了14种纹理特征指标[10]，其中常用的用于提取遥感图像中纹理信息的特征统计量主要有：均值（Mean）、方差（VAR）、一致性（HOM）、对比度（CON）、非相似度（DIS）、熵（ENT）、角二阶矩（ASM）以及灰度相关（COR）等，相应的计算公式如下：

本文基于Matlab7.1平台编程进行图像的灰度共生矩阵和纹理指标的计算。Matlab是一款方便、强大的人机交互式编程工具软件，上手容易，简单易学，而且里面的内核程序允许用户根据实际需要进行修改，这给用户解决实际问题带来了极大的便利。

灰度共生矩阵一般要求4个纹理参数：窗口大小、步长、方向和灰度等级。本文的研究对象是马尾松林，其纹理图像的方向特征不是很明显，为了方便计算，取方向为 45°。根据原始图像（子图）的大小，窗口大小范围为3×3，5×5…，13×13，由于步长随窗口变化而变化，本文暂将步长设为1，2，…，窗口大小－1。灰度等级设为4，8，16，…，256。

为了探求测树因子（平均胸径、平均树高、郁闭度、年龄等）与纹理指标（均值、方差、一致性、对比度等）的相关程度，我们将纹理参数的组合与纹理指标结合起来计算，生成相应的纹理图像，然后取其图像灰度值的平均值作为纹理因子（即自变量），用森林二类调查数据的测树因子作为因变量，分别进行线性回归，统计相关系数的平方，选择每组测树因子中相关系数平方最大值的组合。表1为仙居县森林资源二类调查数据。

表1 仙居县森林资源二类调查数据Table 1 Forest management inventory of Xianju, Zhejiang province

2 结果

根据纹理计算公式和一元线性回归方法，得到了纹理指标与测树因子的最优估算模型（表2，图1）。从中可以看出，用角二阶矩这一纹理指标拟合得到的年龄模型，其相关系数的平方达到了0.5566，表明年龄与角二阶矩的相关程度较高；与之相比，郁闭度与方差的相关程度则并不明显，相关系数的平方只有0.3914。

表2 最优拟合结果Table 2 The optimal fitting result

3 结论与讨论

本文根据高分辨率的航片初步建立了纹理图像与测树因子之间的线性关系，结果表明，用角二阶矩和方差这两个纹理指标拟合年龄、平均树高这两个测树因子效果较好，其相关系数的平方达到了 0.5以上。另外，通过实验发现，灰度等级越高（即图像压缩程度小），对纹理图像与测树因子的关系贡献程度越高，因此我们推测图像的纹理信息越丰富对研究纹理与测树因子的关系越有帮助。

图1 最优模型结果Figure 1 The optimal modeling result

值得注意的是，在航片图像预处理时，需要考虑航片的几何校正，以便于与地面样地数据进行连接，而调查的地面数据同样需要矢量化，这样才方便准确样地的方位并弄清周围的环境。为了确保马尾松纹理图像的信息丰富，我们选择树种类型时，考虑到马尾松纯林含有更多的马尾松信息，而其他的树种类型则有不同程度的树种信息干扰。

由于所用的航片数据含有不少噪声（如等高线、阴影等），这对数据的提取多少产生了一些影响，同时也限制了窗口大小的选择。我们预计，随着窗口的增大，纹理将出现更加规律的周期性变化，纹理与测树因子的关系将更加明显，这对研究它们的关系很有帮助。由于数据量的限制，本文还没做模型的精度验证，这将是以后要做的工作。

[1]宫鹏，黎夏，徐冰，等. 高分辨率影像解译理论与应用方法中的一些研究问题[J]. 遥感学报，2006（1）：1－5.

[2]Rosenfeld A，Kak A. Digital Picture Processing (2nd edit-ion)[M]. Academic Press，1982.

[3]Haralick R M，Shanmugam K，Dinstein I. Texture features for image classification. [J]. IEEE Trans Sys Man Cyber，1973，3（6）：610－621.

[4]李明诗，谭莹，潘洁，等. 结合光谱、纹理及地形特征的森林生物量建模研究[J]. 遥感信息，2006（6）：6－9.

[5]Blinn C E. Estimation of Important Scenic Beauty Covariates from Remotely Sensed Data[D]. Blacksburg, Virginia, 2000.

[6]Ozdemir I, Norton D A, Ulas Ozkan Y,et al. Estimation of Tree Size Diversity Using Object Oriented Texture Analysis and Aster Imagery[J]. Sensors 2008（8）：4709－4724.

[7]沈兵明，金艳. 基于GIS的山地人居环境自然要素综合评价——以浙江省仙居县为例[J]. 经济地理，2006，26（增刊）：305－311.

[8]项云飞，吴卫阳. 马尾松林二种阔叶化改造模式对比试验[J]. 浙江林业科技，2006，26（5）：46－48.

[9]张国庆，黄从德，郭恒，等. 不同密度马尾松人工林生态系统碳储量空间分布格局[J]. 浙江林业科技，2007，27（6）：10－13.

[10]Haralick R M. Statistical and Structural Approaches to Texture[J]. Proc IEEE, 1979, 167（5）：786－804.