基于高空间分辨率遥感数据的天山云杉树冠信息提取研究

刘玉锋，潘 英，李 虎

(1.滁州学院计算机与信息工程学院,滁州 239000;2.安徽省高分辨率对地观测系统数据产品与应用软件研发中心，滁州 239000;3.滁州学院学生处，滁州 239000)

0 引言

树冠是树木的重要组成部分，是提取各种森林参数的主要依据。人们常常通过对树冠的研究来监测树木的长势、预测树木的生物量，判定木材的材性等。树冠的大小也是模拟能量或质量通过树冠传递的众多模型中的重要输入参数之一，对树冠的研究是实现从叶片到林分不同尺度生态转换的关键，对于揭示林木产量及其生长机制、提高森林资源监测和现代化管理水平均具有重要意义[1]。

从遥感图像当中进行信息提取主要采用图像分类的方法，分类依据是不同地物之间的光谱差异，而高空间分辨率遥感图像光谱波段较少，传统基于光谱特征的分类方法不能适用于高空间分辨率遥感图像分类。基于高空间分辨率遥感图像，以对象取代像元为单位进行地物信息提取，可以大大提高提取效率和提取精度。但其中尺度大小的确定对林冠提取的影响很大，如果尺度过小，可能会把一个完整的树冠分得支离破碎；如果尺度太大，则可能出现多个树冠重叠在一起被认成一个树冠。可见，找到能够对应单个树冠恰当的尺度是一件困难的事情，特别是对于林冠层次结构比较复杂的林分，只确定一个尺度可以适应所有树冠是不太可能的。

天山云杉广泛分布于天山及邻近的山地，是中亚山地森林最主要的建群树种之一。天山云杉属乔木，树形高大，树冠窄长，呈中间高而四周低的尖塔状，在遥感图像上会形成树冠顶点位置较亮而树冠边缘较暗的特征。基于该特征，天山云杉立木树冠信息提取可以采用半自动或自动方法进行。目前常用的树冠信息提取方法通常分2步完成，首先是单木位置探测，获取树冠中心位置；然后是树冠轮廓描绘，自动获取树冠边界。单木位置探测方法主要有局部最大值法[2]、多尺度分析法[3]和模板匹配法[4]等；树冠轮廓描绘方法主要有谷地跟踪法[5]、区域生长法[6]、分水岭分割法[7]、局部射线法[8]等。由于森林中树冠大小参差不齐，再加上图像空间分辨率的限制、树冠内在的重叠交叉情况的存在，目前单木树冠自动提取的精度还不高[9]，最好提取精度也只能与人工目视解译精度相当，还无法真正替代实地调查中的每木检尺。

针对以上问题，本文从林木树冠的空间几何特征出发，通过分析天山云杉林木树冠的形态特征及其在遥感图像上的成像规律，结合多尺度斑点检测和梯度矢量流(gradient rector flow，GVF)Snake模型，提出了“树冠顶点探测——树冠轮廓绘制——树冠轮廓优化”的天山云杉树冠信息提取技术流程，以期解决中亚山地森林资源遥感监测和调查中因子信息提取的技术难题。

1 研究区概况及数据源

1.1 研究区概况

研究区位于我国新疆天山西部林区(图1)，地处天山山脉北麓、伊犁河上游，地理位置在E81°25′2.3″～81°38′37.2″，N42°49′12.0″～ 42°57′11.6″之间。该林区地势起伏大，地形结构复杂，上部多冰川，且有以终年积雪为水源的天然水库，中部为高山森林带，下部为伊犁河谷盆地，是天山山地最主要的地带性植被——天山云杉的集中分布区，同时也是中亚山地森林的主要分布区。

图1 研究区范围示意图Fig.1 Schematic diagram of research area

1.2 数据源及其预处理

研究选择美国数字地球公司(Digital Globe)的WorldView-2高空间分辨率卫星遥感数据作为主要数据源，采用标准等级的预正射校正产品(ORStandard2A级)，其中多光谱波段的空间分辨率为2 m，全色波段的空间分辨率为0.5 m，成像时间为2010年8月30日，图像中云的含量较少，图像清晰，质量较好。

为了充分利用高空间分辨率卫星遥感图像的光谱和空间信息，采用NNDiffuse方法[10]对全色波段和多光谱波段进行融合，融合后的图像空间分辨率为0.5 m，包含红、绿、蓝、近红外等4个光谱波段信息。

与此同时，还开展了天山云杉森林资源样地补充调查，即在研究区范围内选择不同林龄、不同立地条件和不同密度的林分，设置标准样地若干，并对样地作包括树高、胸径、冠幅、郁闭度等内容在内的每木检尺工作，获取的样本数据用作检验数据。

2 研究方法

天山云杉树形高大，树冠窄长，在高空间分辨率卫星遥感图像上的成像呈圆形或近圆形，且与背景相比，树冠具有较高的光谱反射值，背景(包含阴影在内)具有较低的光谱反射值。基于此，提出了“树冠顶点探测——树冠轮廓绘制——树冠轮廓优化”的技术流程。

2.1 基于多尺度斑点检测的树冠顶点探测

天山云杉树冠在高空间分辨率卫星遥感图像上的成像呈圆形或近圆形(星下点位置)，其光谱反射率呈现中间高而四周低，且在树冠顶点位置处达到最大值的特点。因此采用多尺度斑点检测的方法来确定树冠的位置[11]。即通过构造图像尺度空间，得到在不同尺度下的滤波图像，再对滤波后的图像进行极值点检测，得到同时在位置空间和尺度空间的极值点作为树冠顶点。

1)尺度空间构建。如果待检测物体的大小未知，为了能够准确识别该物体，需要在图像信息处理模型当中引入尺度因子，通过连续变化尺度参数获得不同尺度下的视觉处理信息序列，以此构建图像尺度空间。图像的尺度空间L(x,y,σ)被定义为高斯函数G(x,y,σ)与原图像I(x,y)的卷积，即

L(x,y,σ)=G(x,y,σ)×I(x,y)，

(1)

式中：(x,y)为空间坐标；σ为尺度参数。在图像多尺度空间当中，通过改变特定尺度的σ值，可检测得到不同大小的斑点。通过对研究区范围内天山云杉立木参数实测数据统计得知，天山云杉树冠半径最小值为0.75 m，最大值为4.66 m，由此确定σ值的大致范围在0.5～3.295之间。这样避免了在进行斑点检测时把尺寸过大或过小的其他地物误认为树冠。

2)极值点检测。为了从图像尺度空间中检测到天山云杉树冠顶点，需要将图像中的每一个像素点与它周围的点进行比较，即将中间的待检测点同与它同尺度的8个相邻点以及上下相邻尺度对应的9×2个点(共计26个点)进行比较，如图2所示，如果中间待检测点的高斯响应值比它周围的26个相邻点都大或者都小，那么这个点即为尺度空间中的极值点，也就是天山云杉树冠顶点。天山云杉的林下主要为稀疏草地或裸地，因此采用归一化植被指数(normalized difference vegetation index,NDVI)阈值法对检测到的极值点进行非云杉滤除。本文通过多次实验尝试，将NDVI的阈值设定为0.12，即滤除NDVI值小于0.12的极值点，消除林下植被对信息提取结果的影响。

图2 尺度空间极值点检测示意图Fig.2 Schematic of extreme points detection in scale space

3)极值点精确定位。斑点检测是基于离散尺度空间的极值点检测，检测到的极值点在连续空间中可能并不一定也是极值点，图3显示了二维函数离散空间得到的极值点与连续空间极值点的差别。因此，当检测到树冠顶点的所在位置后，仍需采用子像素插值法对树冠顶点进行精确定位，即通过对尺度空间采用高斯函数进行曲线拟合，从而确定极值点的确切位置和尺度[12]。将高斯函数在特征点处二阶Taylor展开，即

(2)

(3)

图3 离散空间与连续空间极值点的差别示意图Fig.3 Extreme points’ difference between discrete space and continuous space

2.2 基于标记分水岭变换的树冠轮廓确定

天山云杉树冠轮廓是天山云杉及其背景(阴影或下垫面草地)在图像上的像元灰度值(即光谱反射率值)发生突变而形成的局部特征不连续性，表现为沿着边界走向的像元灰度值变化平缓，而垂直于边缘方向的像元灰度值变化剧烈，最终形成沿着边界的圆形或近圆形的闭合曲线。本文采用基于标记分水岭变换的区域分割方法得到树冠的初始轮廓，具体步骤如下：

1)梯度信息计算。待分割图像的边缘像素往往具有较大的梯度值，对应于地表上的“分水岭脊线”，而每个区域的内部像素通常具有较小的梯度值，对应地表上的“集水盆地”。分水岭变换的目的就是要求出梯度图像的“分水岭脊线”，因此，基于分水岭变换的图像分割方法，其性能在很大程度上依赖于用于计算待分割图像梯度的算法。采用基于数学形态学的梯度计算，定义为

grad(f)=(f⊕b)-(f⊖b)，

(4)

式中：f为输入图像；⊕和⊖分别代表形态学中的“膨胀”运算和“腐蚀”运算；b为结构元素，本方法中选择圆盘状结构元素，它具有各向同性，可以消除梯度对边缘方向的依赖性。

2)梯度图像重构。采用形态学开闭重建运算对梯度图像进行重构，目的是消除梯度图像中由于非规则灰度扰动和噪声引起的局部极值，从而保留天山云杉轮廓位置的极值信息。首先采用形态学开重建运算，消除梯度图像中尺度比结构元素小的极大值噪声和非规则干扰，随后再进行形态学闭重建运算，去除比结构元素小的暗噪声及非规则干扰。通过以上方法，区域极大值和极小值均得到了修正。

3)目标地物标记。标记应与天山云杉树冠的特性相关，即将天山云杉的所在位置作为前景标记，这里直接采用多尺度斑点检测方法得到的树冠顶点探测结果作为前景标记。

4)梯度图像修正。树冠位置标记后，就可以使用强制最小技术，使局部最小区域仅出现在标记的位置，并把其他像素按需要进行“上推”，以便删除其他的局部最小区域[13]，使梯度图像上对应于目标地物(天山云杉)的谷值和对应于目标边界(树冠轮廓)的峰值被保留下来。

5)分水岭变换。采用浸没模型算法(V-S算法)[14]，即将图像看作地形表面，图像中每一个像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域代表一个积水盆地，而集水盆的边界即为分水岭。通过变换，可在每个树冠顶点的四周，得到连通的、封闭的、单像素宽的树冠边缘，即为天山云杉树冠的初始轮廓。

2.3 基于GVF Snake演化的树冠轮廓优化

通过标记分水岭变换得到的树冠轮廓，是变换结果图像中灰度值为0的像元集合，即以像元中心点连线构成的闭合曲线作为树冠轮廓，这一结果实际是以栅格形式表现的地物轮廓形状，由于受制于像元大小，还略显粗糙，不够精准，需要进一步优化。研究采用GVF Snake主动轮廓模型对树冠轮廓进行优化，将树冠轮廓看作是在能量泛函引导下的自由形式的变形曲线，在轮廓连续性、光滑性等简单约束条件下，通过能量泛函的最小化过程，使得树冠轮廓曲线逐渐变形，并向目标边界演化，最终收敛到目标位置[15]。具体步骤为：

1)GVF力场生成。将GVF力场定义为V(x,y)=[u(x,y),v(x,y)]，其中，u(x,y)和v(x,y)是它的2个分量，分别表示图像灰度在x,y方向上的变化，树冠轮廓曲线的能量函数可表示为

(5)

式中：(x,y)为图像的坐标；f(x,y)为梯度场；ux,uy,vx,vy分别为u和v对x和y的一阶偏导；μ为控制参数，其值取决于图像中噪声的大小，当图像中噪声比较大时就增加μ值，反之就减小μ值。能量函数由2项组成，第一项称为平滑项，作用是产生一个缓慢变化的矢量场V(x,y)；第二项|f(x,y)|2|V(x,y)-f(x,y)|2称为数据项，它使力场V(x,y)尽可能地接近图像的边缘，即f(x,y)。当f(x,y)较小时，说明处于图像的均质区，能量函数主要由平滑项主导，极小化泛函时产生一个变化缓慢、光滑的力场；当f(x,y)较大时，说明处于图像的边缘区，能量函数主要是由数据项来主导，当V(x,y)=f(x,y)时，取得最小值。

2)初始轮廓曲线的设定。基于标记分水岭变换的树冠轮廓绘制结果虽然比较粗糙，还不够精准，但是基本处在GVF力场的有效逼近域内，因此可将其作为GVF Snake的初始轮廓曲线。

3)能量泛函最小化。为使矢量的能量函数最小，必须满足欧拉方程，即

μ2V-(V-f)|f|2=0，

(6)

分解式为

(7)

(8)

(9)

式中：C(x)和C(y)分别为轮廓线上控制点的横、纵坐标矩阵；γ和K为计算矩阵；inv为参数矩阵。通过分析GVF力场可知，在有效逼近域内，V(x,y)总是指向能量最小处，即边缘的位置。因此，初始轮廓线会在V(x,y)的作用下，向真实边缘逼近，当轮廓的总能量达到最小值时停止逼近，得到的轮廓即为目标结果轮廓。

3 结果与分析

3.1 提取结果

天山云杉多为天然更新林木，立地条件类型的多样性，使其形态不一、结构多样。为了验证本论文提出的天山云杉树冠信息提取技术的适用性，分别将其应用到不同林分密度的天山云杉样地当中。并且为了便于后续将树冠信息提取的结果与森林资源样地补充调查的实测资料进行对比，对所选择的样地按照郁闭度的密、中、疏的程度进行了归类，从中选择了具有代表性的样地，样地大小为60像素×60像素，如图4所示。其中绿色线框表示样地的范围，由于边缘部位的树冠会出现不完整的情况，故影像范围在样地的基础上，向外扩展了40像素。当林分郁闭度较高时，如图4(a)，(b)，(c)所示的样地20120911，郁闭度达到了0.7，林木分布均匀、密集。采用基于斑点检测的树冠顶点探测时，尽管大多数的树冠得以探测，树冠轮廓也能较好地被绘制出来，但还是可以看出粘连在一起的多株树冠或是由于形状不规则而被漏探测，如图4(b)中白色虚线框所示；或是由于无法区分单株树而形成一个“大的树冠”被探测出来，如图4(b)中绿色虚线框所示，这些现象导致了被探测到的树冠数目要比实际的少。当林分郁闭度适中时，如图4(d)，(e)，(f)所示的样地20120909，郁闭度为0.6。从图上可以看出，样地内绝大多数的天山云杉的树冠轮廓可以被绘制出来，并且树冠轮廓清晰可辨，能够较好地反映实际情况。当林分郁闭度较低时，如图4(g)，(h)，(i)所示的样地20120910，郁闭度为0.4，尽管大多数的树冠得以探测，树冠轮廓也能较好地被绘制出来，但由于样地中的云杉分布不均，导致用于识别树冠的背景在遥感图像上表现不一，从图上可以看出，林下背景既有低覆盖度草地，又有高大林木所成的阴影，尤其在低覆盖度草地与阴影结合的区域，容易出现伪树冠顶点，如图4(h)中黄色虚线框所示，从而使得探测到的树冠数目多于实际量测的数目。

(a)密林原始图像 (b)密林树冠顶点 (c)密林树冠轮廓

(d)中林原始图像 (e)中林树冠顶点 (f)中林树冠轮廓

(g)疏林原始图像 (h)疏林树冠顶点 (i)疏林树冠轮廓

图4 不同郁闭度样地天山云杉树冠轮廓分布

Fig.4CrownoutlinedistributionofsparsePiceaschrenkianavar.tianschanica

3.2 精度分析

3.2.1 立木株数

在天山云杉树冠信息提取的技术流程中，至关重要的一个步骤就是树冠顶点的识别，其识别效果的好坏直接决定树冠轮廓绘制的准确性，并且影响到后续冠幅遥感估算的精度。为此，本文针对树冠顶点识别的精度做了分析，把采用多尺度斑点检测方法得到的树冠顶点探测结果与森林资源样地补充调查的现场实测株数进行对比分析，表1列出了把识别株数与每个地面调查样地株数对照得出的正确、遗漏、误判情况。从表1可以看出，采用多尺度斑点检测方法的天山云杉林木识别数量与实测数量的散点较接近，准确率为82%；但该方法估算得到的天山云杉立木株树仍然存在漏分和误判，占比分别为18%和12%，这部分误差的绝对值虽然不高，但降低了天山云杉单木树冠识别的精度。

表1 样地识别株数与地面调查株数对照表Tab.1 Table of tree numbers of auto identification and field work from different stands

3.2.2 林木冠幅

基于提取得到的树冠轮廓，通过计算树冠南北向和东西向的冠幅长，然后取平均值作为其冠幅值，将此参数作为树冠信息提取的结果进一步与现场实测值进行对比分析，发现差异不大，但估算值普遍低于地面实测值，只有极少数高于地面实测值。原因在于对天山云杉冠幅的估算是树冠表面未被其他树冠遮挡部分的尺寸，如果两株树相邻，树冠彼此都被遮挡，则得到的树冠边界介于二者之间，即只是估算了树冠的阳性冠，而不是整个树冠。树冠被遮挡部分在图像上很难被识别，但实际上，外业测量时往往是将树冠上的最远枝作为冠幅的起测点，而不管树冠是否被遮挡。进一步对不同林分密度的天山云杉冠幅估算结果进行分析发现：当林分郁闭度较高时，如图5(a)为郁闭度为0.7的标准地20120911天山云杉(密林)冠幅估算结果与实测值的散点图，可以看出，冠幅估算的平均绝对误差为0.523 m，均方根误差为0.598 m，平均误差为10.8%，其中，最小误差绝对值为2.6%，最大误差绝对值为21.5%。当林分郁闭度适中时，如图5(b)为郁闭度为0.6的标准地20120909天山云杉(中林)冠幅估算结果与实测值得散点图，可以看出，郁闭度为0.6的标准地20120909冠幅估算的平均绝对误差为0.144 m，均方根误差为0.171 m，平均误差为4.5%，其中，最小误差绝对值为0.2%，最大误差绝对值为11.4%。当林分郁闭度较低时，如图5(c)为标准地20120910天山云杉(疏林)冠幅估算结果与实测值得散点图，可以看出，郁闭度为0.3的标准地20120910冠幅估算的平均绝对误差为0.232 m，均方根误差为0.261 m，平均误差为6.4%，其中，最小误差绝对值为0.7%，最大误差绝对值为18.5%。总体而言，本文方法对高、中、低郁闭度的天山云杉树冠信息都有较好的提取效果。

(a)标准地20120911,郁闭度0.7 (b)标准地20120909，郁闭度0.6 (c)标准地20120910,郁闭度0.3

图5 天山云杉树冠信息提取估算精度分析

Fig.5AccuracyanalysisincrownwidthestimationofdensePiceaschrenkianavar.Tianschanica

4 结论与讨论

4.1 研究结论

针对基于光谱特征的分类方法不能很好地适用于高空间分辨率遥感图像分类的问题，本文从天山云杉林木树冠的空间几何特征入手，通过分析林木树冠的形态特征及其在遥感图像上的几何形态特征，将多尺度斑点检测和GVF Snake主动轮廓模型有序地结合在一起，提出了“树冠顶点探测——树冠轮廓绘制——树冠轮廓优化”的树冠信息提取技术流程，解决了标记分水岭变换目标标记难以准确设定、主动轮廓模型演化结果受制于轮廓线的初始位置的问题，从而得到了定位准确、边界简洁的树冠轮廓结果。经与调查样地中每一株树的实测冠幅值进行比较，结果表明：该方法对高、中、低郁闭度的天山云杉树冠信息都有较好的提取效果，平均误差分别为10.8%，4.5%和6.4%。结果表明本研究能够较好地解决森林资源监测中高空间分辨率遥感数据树冠判读的关键技术问题，对中亚山地森林调查因子的遥感识别与信息提取具有一定的应用推广价值。

4.2 问题及讨论

本研究仍存在以下不足有待于进一步完善和深入研究。

1)基于多尺度斑点检测方法探测树冠顶点时，采用了各向同性的检测算子，这对于天山云杉的遥感成像呈圆形或近圆形固然是有效的，即当卫星高度角较大时，该方法是可行的。但随着太阳高度角的逐渐减小，具有一定高度的天山云杉在图像上不再是从树冠正上方拍摄的圆形，而是对树冠侧面成像的椭圆形或者锥形，此时仍采用该方法进行树冠顶点的探测将无法得到理想的结果，所以该方法存在一定的局限性。

2)对于高空间分辨率卫星遥感图像而言，如果不能消除由于高度差影响而造成的像点位移误差，那么通过采用多尺度斑点检测方法探测得到的树冠顶点的位置只能是一个粗略的位置，要想得到准确的位置，仍需要考虑成像几何关系对树冠顶点位置偏移影响。

3)阴影对天山云杉树冠信息的提取带来了不小的困难，对天山云杉冠幅估算的精度影响也较大。一方面，处于图像阴影当中的地物光谱特征受阴影噪声模糊作用的影响而变得不明显，使其中一部分树冠不能被有效识别，从而降低了树冠信息提取的精度；另一方面，由于阴影区域表现为深色或暗色图斑，与周围较亮区域形成明显反差，容易被识别为伪树冠。然而目前阴影的消除仍是个难点，需进一步开展深入研究。