基于光谱局部最大值滤波的林分密度估计

王书涵，张晓丽，朱程浩，瞿 帅

（北京林业大学 林学院 省部共建森林培育与保护教育部重点实验室，北京100083）

基于光谱局部最大值滤波的林分密度估计

王书涵，张晓丽，朱程浩，瞿 帅

（北京林业大学 林学院 省部共建森林培育与保护教育部重点实验室，北京100083）

林分密度对林分生长有很重要的影响，既是生态学中重要的研究内容，也是林分因子调查的重要参数之一。采用快鸟影像的全色波段利用局部最大值滤波方法提取了研究区的林分密度。采用了皮尔森相关系数衡量了实际林分密度同树冠点数量之间的相关性，重点探讨了3×3，5×5，7×7（以像素为单位）等3种不同窗口大小及不同的归一化植被指数（INDVI）阈值对提取树冠点数量的影响，选择出最佳的窗口和INDVI滤值的组合，并建立线性回归模型，将整个研究区划分成样地大小的格网，统计格网中光谱最大值点的数量并转换成林分密度栅格图层，运用建立的模型得到研究区林分密度。经实验发现：采用光谱局部最大值滤波方法提取出的树冠点数量确实同实际林分密度存在较强的相关性（R2＝0.545 5，ERMSE＝13.97，P＜0.001），特别是针叶林，经F检验采用3×3窗口大小，INDVI≥0.2作为阈值具有极显著的相关性并得到最高的相关系数（R2＝0.741 5，ERMSE＝14.45，P＜0.01）；阔叶林较针叶林相关系数略低（R2＝0.442 2，ERMSE＝10.97，P＜0.01），并采用5×5窗口大小以及INDVI≥0.2作为阈值达到最佳的效果；最后利用建立的模型生成了研究区的林分密度分布图。光谱最大值法能较好地提取林分密度。图5表6参27

森林测计学；快鸟影像；局部最大值滤波；林分密度

森林结构参数的准确获取是森林资源规划调查的重要工作之一［1-4］。林分密度与林分的树冠、胸径、树高的生长有显著的相关性［5-6］，是森林资源二类调查中一项重要的参数和指标［7-9］。同时，林分密度还与人工林的木材材性、生物量、蓄积量等有着密切的关系［10-11］，林分密度是表征森林生态结构，影响生态系统特征，以及林业可持续发展的关键因子。森林测计学中林分密度可分为疏密度、株数密度和郁闭度。本研究所讨论的林分密度仅指株数密度，为单位面积株数/单位样地面积。通常林分密度的获取方式主要包括样地测量和遥感手段估测［12］。常规的测量林分密度的方式一般是采用样地调查方法以布设样地为主，需要耗费较多的人力物力且不能实现大区域的估测，长期复查容易导致误差和重复测量［13-14］。遥感技术的推广应用给地区尺度进行林分密度的估测提供了有力的工具［15］。随着高分辨率遥感应用越来越广，探讨高分遥感提取森林植被参数就具有重要意义。基于高分遥感提取林分密度的方法有：二维各向同性高斯核函数平滑滤波［8］、分水岭方法［16］、发射或辐射的反演模型方法［17］、局部变化模型方法［18］和光谱最大值滤波法［19-20］。光谱最大值滤波是一种可以用来检测单木位置的技术，它是基于针叶树冠的顶点会存在较大遥感影像像元亮度值（DN值）的假设［21］，优势是能快速有效地探测到林地单木的位置［22］。因此，光谱最大值法可以用来间接估算林分密度。光谱最大值法的普适性需要进一步研究，特别是随着地域不同，林分类型不同时究竟采用多大的窗口以及应该采用何种方法精炼树冠点都是值得深入研究的问题，而本研究从此出发探究了不同的窗口大小和归一化植被指数（INDVI）阈值的光谱最大值法应用于不同的林分类型提取林分密度的精度问题。

1 研究区概况与数据

1.1 研究区域

鹫峰森林公园坐落于北京市海淀区西北部苏家坨镇境内，是距离京城最近的国家森林公园之一。鹫峰经纬度大约为40°3′1.618″N，116°2′25.722″E，横跨海淀和门头沟2个区，总面积约为832.04 hm2，主要的林分类型有侧柏Platycladus orientalis林、油松Pinus tabuliformis林、落叶松Larix gmelinii林、刺槐Robinia pseudoacacia林、栓皮栎Quercus variabilis林、栓皮栎与槲栎Quercus aliena混交林等，其中针叶林是主要的森林类型，面积为474.85 hm2，占整个林区的57.1%［23-24］。公园内地形多变，最高海拔为1 153 m，属于华北大陆季风气候，具有冬春季干燥多风、夏季凉爽多雨的特征［25］。

1.2 数据获取与处理

QuickBird影像是美国数字全球公司QuickBird卫星获取的图像，全色图像分辨率为0.61 m，多光谱图像为2.44 m。本研究采用QB LV2A数据，获取时间为2008年10月24日，位深度为16 Bit，太阳高度角为37.3°，太阳方位角为166.6°，卫星方位角为74.6°，卫星高度角为67.3°，图像本身已经经过了几何粗校正。由于实验区地形起伏明显，采用北京市测绘局提供的1∶2 000大比例尺地形图结合提取出的高精度数字高程模型（DEM）进行了正射校正，覆盖了主要的研究区域。以匹配地面实测样地的位置，保证建模的精度。

1.3 外业调查数据

研究小组收集了具有代表性的72个标准样地，样地大小为20 m×20 m，样地基本上均匀分布在研究区域中。调查时间分别是2006年和2011年，调查方法采用全站仪和手持差分全球定位系统（GPS）相互配合的方式精确定位样地的4个角的位置具体流程可参见样地定位方法［26］，并且记录下每株树的树种、胸径、树高、冠幅等信息。根据实际调查样地的树种组成，针叶树种大于0.65的划分为针叶林，阔叶树种比例大于0.65的划分为阔叶林［27］，最终划分为25个针叶林样地，43个阔叶林样地，样地的冠幅描述性统计信息参见表1和表2。

表1 针叶林样地树冠大小描述性统计Table 1 Crown size statistiques in coniferous sample plot

表2 阔叶林样地树冠大小描述性统计Table 2 Crown size statistics in broadleaved sample plot

2 研究方法

2.1 图像预处理

首先将QuickBird全色波段和多光谱波段分别进行预处理操作，首先进行基于控制点的正射校正，由于本研究获取了鹫峰区域的1∶2 000大比例尺的地形图，由北京市测绘部门提供，地形图已经栅格化且经过了精校正，因此用地形图直接校正快鸟影像。DEM的获取采用如下方式：将等高线的线图层在ArcGIS地形分析模块中将其转换成不规则三角网（TIN），再将TIN转换成1 m精度的DEM。以二类调查数据为依据先确定林地区域，在QuickBird影像上裁剪出研究区。

2.2 窗口选择与INDVI阈值的选择

当一个像素的直径（DN）值比其他给定窗口大小的周围像素的DN都要大的时候被定义为局部最大值［8］。光谱最大值提取采用了ERDAS中的聚集分析方法，能够对这些感兴趣的像素执行的操作包括标准差、和、平均值、中值、最小值、最大值等。可以通过选择窗口的大小评价围绕在感兴趣的中央像素周围的区域。

本研究首先对经过预处理的快鸟影像的全色波段分别用3×3，5×5，7×7等3种不同窗口大小（除特殊指出外，以下所有窗口均以像素为单位）进行光谱局部最大值滤波，接着用全色波段减去滤波后的图像，值为0的点则为局部光谱最大值点。将过渡图像中的值为0的点提出来并转换成点图层，利用INDVI对非树冠点进行剔除。本研究采用以下方式对树冠点进行筛选：以0.1为梯度将INDVI图层按照0.1～0.5划定为5个等级，当设定INDVI阈值为大于等于0.1时，若此时的光谱最大值点的INDVI值小于0.1则剔除，剩下的点进行下一步的统计。

将野外获取的精确的样地的位置叠加在不同窗口滤波得到的树冠点图层上，空间关联出落在每个样地中的树冠点的数量，以光谱最大值点作为自变量，实际林分密度作为因变量进行皮尔森相关性分析，这样不同的窗口和INDVI阈值组合方案会得到不同的相关系数值，由此评价不同林分的最佳的窗口和INDVI组合方式以建立林分密度估测模型。

经过过滤处理之后的树冠点图层按照20 m×20 m大小的方格栅格化，以每个方格为统计单位统计落入方格中的局部最大值点的数量，从而得到初始林分密度，将初始林分密度作为自变量输入建立的回归模型中，便能够计算出研究区域实际林分密度。图1显示了整个研究操作的具体流程。所有的空间分析和统计分析步骤均在Erdas，ArcGIS和R语言中完成。

3 研究结果

随着滤波窗口大小的变化，图像逐渐变得模糊，结果如图2所示。图A′，图2B′，图2C′为利用全色图像减去滤波后图像的效果示意图。由于截取的区域土壤和树冠的差异较为明显，可以从图2中比较清晰地分辨出树冠以及道路的边缘，有些较大的单木树冠清晰可辨。图3以B11样地为例显示了将样地准确叠加在经过处理过后的光谱最大值点的图层上，可以清楚地统计出样地范围内所包含的单木数量。

3.1 窗口大小和INDVI阈值的选择

根据研究方法中的介绍，分别选取了不同的窗口大小及INDVI阈值提取林分密度，其中根据分析处理结果见表3～表5。

图1 技术流程图Figure 1 Flow Chart

图2 不同窗口大小的局部最大值滤波的效果图Figure 2 Sketch map of local maximum filtering using different window size

图3 样地位置同获取的光谱最大值点叠加效果示意图Figure 3 Schematic map ofB11sample plot and local maximum points

由表3可以看出：对于所有林分，决定系数的最高值出现在5×5窗口大小，这跟阔叶林的分析结果一致，而针叶林的决定系数的最高值出现在3×3的窗口大小中。随着INDVI阈值的升高，其相关系数的数值主要呈现先上升后下降的趋势。以针叶林为例，无论采用何种窗口，INDVI阈值为0.1时相关性较高，随着INDVI增加到0.2，决定系数有所增加，当INDVI阈值为0.3，决定系数开始下降，当INDVI取值继续增加到0.5，决定系数继续下降。针叶林的决定系数最高值出现在3×3窗口，INDVI选择为≥0.3。另外，从决定系数的数值来看，无论选择哪种INDVI阈值，几乎总存在一个固定的窗口是拟合的最佳窗口选择。例如对于阔叶林，无论采用哪种INDVI阈值，5×5的窗口大小获得的决定系数最高。由此可见，窗口大小是影响林分密度的估计最主要的因素。

表3 不同窗口和INDVI阈值的光谱最大值滤波的林分密度与实际林分密度相关性比较（所有林分）Table 3 Correlation analysis of all stand using the based on the local maximum of different INDVIand windows sizes（all stand）

表4 基于不同窗口大小和INDVI阈值的光谱最大值滤波的林分密度与实际林分密度相关性分析（针叶林）Table 4 Correlation analysis of all stand using the based on the local maximum of different INDVIand windows sizes（coniferous stand）

表5 基于不同窗口大小和INDVI阈值的光谱最大值滤波的林分密度与实际林分密度相关性分析（阔叶林）Table 5 Correlation analysis of all stand using the based on the local maximum of different INDVIand windows sizes（broadleaved stand）

3.2 林分密度估测模型及估测结果

根据上述分析，分别选取了较高的相关系数所对应的INDVI阈值和窗口大小，建立所有林分、阔叶林、针叶林的统计模型。统计结果如表6所示。

表6 不同林分类型的相关性分析Table 6 Correlation analysis in different forest types

图4显示了针叶林样地和阔叶林样地最佳INDVI阈值和最佳窗口选择下的林分密度提取散点图，纵坐标代表了实际林分密度，横坐标代表了基于光谱局部最大值滤波提取出的林分密度，从而建立线性回归模型。以针叶林和阔叶林范围生成的格网为单位，将模型应用于针叶林和阔叶林的光谱最大值点数量图层的栅格图，从而得到相应的林分密度分布图，其中针叶林和阔叶林的范围可以根据二类调查小班数据来确定。分析结果见图5，其中红色区域代表相应林分的林分密度为0。

4 结论

本研究利用QuickBird的全色波段采用光谱局部最大值的方法针对不同的林分类型提取了研究区的林分密度，结论如下：不区分林分类型（将针叶林阔叶林同时考虑），采用5×5的滤波窗口以及采用INDVI≥0.2作为阈值过滤树冠点并拟合实际的林分密度，能达到最高的相关性（R2＝0.545 5,ERMSE＝13.97）。针叶林样地采用3×3窗口大小以及采用INDVI≥0.2作为阈值建立的模型精度最高（R2＝0.741 5,ERMSE＝10.97）。研究结果进一步表明了窗口大小的选择对光谱最大值法提取林分密度有重要的影响，本研究分林分类型进一步提高了林分密度的估测精度。

图4 不同林分类型的局部最大值林分密度提取散点图Figure 4 Scatter diagram between local maximum and real tree density in different forest type

图5 林分密度估测结果Figure 5 Stand density estimation result

5 讨论

本研究通过利用光谱最大值滤波方法应用于高空间分辨率遥感影像QuickBird试图实现数字化提取研究区内的林分密度流程，取得了较好的结果，其中针叶林的林分密度提取效果最好，建立的针叶林林分密度提取模型可以用来估测研究区域的林分密度，为进一步实现大尺度林分密度提供基础，而阔叶林的估测精度稍低，但模型仍然有利用价值。光谱最大值滤波方法在提取树冠中心点是假设每个树冠只有一个光谱反射值最大点为前提来推测的，因此，可以推测由于针叶林树冠具有规则的树冠形状，通常针叶树种的树冠往往只有一个光谱最大值点，而阔叶树由于具有较大面积的树冠和复杂的树冠结构往往不止一个树冠反射率最大点，这是造成针叶林提取精度高而阔叶林提取精度偏低的主要原因，而我们的研究结果跟前人的研究结论是一致的［19,21］。注意到本研究得出结论阔叶林的林分密度选用5×5的窗口效果最优，这刚好对应于3 m×3 m的树冠范围，这与表2中阔叶林中样地的平均树冠大小较为接近，而在估测针叶林林分密度时选用3×3的窗口大小效果最优，也跟表1中针叶林中样地的平均树冠大小一致，这说明窗口大小的选择应该要贴近于树冠的真实大小才能获得较高的精度。

本研究仍然有几点值得探讨：①利用0.61 m分辨率的全色波段提取的光谱最大值点能够反映出亚米级的树冠位置，其结果跟影像的太阳高度角、太阳天顶角、卫星高度角、卫星天顶角等参数关系密切，因此，相同类型的传感器在不同时相得到的影像采用相同的方法结果是否有偏差还需进一步验证，因此可以利用不同时相的影像进行实验，分析验证本方案的可行性，进一步修正本研究方法。这是未来的研究内容之一。②由于近红外波段与植被的关系密切，还应该充分挖掘近红外波段的潜力，例如采用近红外波段提取光谱局部最大值是否会得到更好的结果还有待研究。③本研究所提出的窗口大小的选择和INDVI过滤值的设定方案是否具有普适性还有待验证，研究相同的研究区不同的数据源是否需要设置不同的INDVI阈值和不同的窗口大小，及可变的窗口大小进行光谱最大值滤波进而探测树冠的位置应该是今后另一个提高精度的研究方向。

6 致谢

在论文完成之际，向北京林业大学林学院国家高技术研究发展计划（“863”计划）项目组所有人员在野外样地调查中所付出的努力表示衷心感谢！感谢相关部门提供的数据支持！

［1］ KAYITAKIRE F,HAMEL C,DEFOURNY P.Retrieving forest structure variables based on image texture analysis and ikonos-2 imagery［J］.Remote Sens Environ,2006,102（3/4）：390-401.

［3］ PASHER J,KING D J.Multivariate forest structure modelling and mapping using high resolution airborne imagery and topographic information［J］.Remote Sens Environ,2010,114（8）：1718-1732.

［4］ van LEEUWEN M,NIEUWENHUIS M.Retrieval of forest structural parameters using lidar remote sensing［J］.Eur J For Res,2010,129（4）：749-770.

［5］ BURKES E C,WILL R E,BARRON-GAFFORD G A,et al.Biomass partitioning and growth efficiency of intensively managed Pinus taeda and Pinus elliottii stands of different planting densities［J］.For Sci,2003,49（2）：224-234.

［6］ 段劼，马履一，贾黎明，等.北京地区侧柏人工林密度效应［J］.生态学报，2010，30（12）：3206-3214.

DUAN Jie,MA Lüyi,JIA Liming,et al.The density effect of Platycladus orientalis plantation in Beijing area［J］.Acta Ecol Sin,2010,30（12）：3206-3214.

［7］ COPENHAVER P E,TINKER D B.Stand density and age affect tree-level structural and functional characteristics of young,postfire lodgepole pine in Yellowstone National Park［J］.For Ecol Manage,2014,320（5）：138-148

［8］ DRALLE K,RUDEMO M.Stem number estimation by kernel smoothing of aerial photos［J］.Can J Fot Res,1996,26（7）：1228-1236

［9］ MCCOMBS J W,ROBERTS S D,EVANS D L.Influence of fusing lidar and multispectral imagery on remotely sensed estimates of stand density and mean tree height in a managed loblolly pine plantation［J］.For Sci,2003,49（3）：457 -466.

［10］ 贾全全，罗春旺，刘琪璟，等.不同林分密度油松人工林生物量分配模式［J］.南京林业大学学报（自然科学版），2015，39（6）：87-92.

JIA Quanquan,LUO Chunwang,LIU Qijing,et al.Biomass allocation in relation to stand density in Pinus tabuliformis plantation［J］.J Nanjing For Univ Nat Sci Ed,2015,39（6）：87-92.

［11］ 张华林，李天会，吴志华，等.不同林分密度对桉树幼林木材材性的影响［J］.中南林业科技大学学报，2010，30（6）：85-91.

ZHANG Hualin,LI Tianhui,WU Zhihua,et al.Effects on wood properties of young Eucalypt plantation in different stands density［J］.J Cent South Univ For Technol,2010,30（6）：85-91.

［12］ OZDEMIR I,KARNIELI A.Predicting forest structural parameters using the image texture derived from worldview-2 multispectral imagery in a dryland forest,Israel［J］.Int J Appl Earth Obs,2011,13（5）：701-710.

［13］ OMULE S A Y.Personal bias in forest measurements［J］.Forest Chron,1980,56（5）：222-224.

［14］ MCROBERTS R E,HAHN J T,HEFTY G J,et al.Variation in forest inventory field measurements［J］.Can J ForRes,1994,24（9）：1766-1770.

［15］ 杜文峰，王凤臻，李庆.提高郁闭度调查精度的几点建议［J］.林业资源管理，1999（3）：63-65.

DU Wenfeng,WANG Fengzhen,LI Qing.Some advices on how to improve the stand density［J］.For Resour Manage, 1999（3）：63-65.

［16］ PERRIN G,DESCOMBES X,ZERUBIA J.2D and 3D vegetation resource parameters assessment using marked point processes［C］//TANG Y Y,WANG S P,LORETTE G,et al.ICPR’06 Proceedings of the 18th International Conference on Pattern Recognition.Washington D C：IEEE Computer Society,2006：1-4.

［17］ SCARTH P,PHINN S.Determining forest structural attributes using an inverted geometric-optical model in mixed eucalypt forests,southeast Queensland,Australia［J］.Remote Sens Environ,2000,71（2）：141-157.

［18］ BA S,CAVAYAS F.Automated forest structure mapping from high resolution imagery based on directional semivariogram estimates［J］.Remote Sens Environ,1997,61（1）：82-95.

［19］ WULDER M,NIEMANN K O,GOODENOUGH D G.Error reduction methods for local maximum filtering of high spatial resolution imagery for locating trees［J］.Can J Remote Sens,2002,28（5）：621-628.

［20］ NELSON T,BOOTS B,WULDER M A.Techniques for accuracy assessment of tree locations extracted from remotesly sensed imagery［J］.J Environ Manage,2005,74（3）：265-271.

［21］ WULDER M A,WHITE J C,NIEMANN K O,et al.Comparison of airborne and satellite high spatial resolution data for the identification of individual trees with local maxima filtering［J］.Int J Remote Sens,2004,25（11）：2225-2232.

［22］ 韩学峰.基于高分辨率遥感林分调查因子的提取研究［D］.福州：福建师范大学，2008.

HAN Xuefeng.A study on the Extraction of the Stand Description Factors from High Resolution Remote Sensing Images［D］.Fuzhou：Fujian Normal University,2008.

［23］ 赵永泉，彭道黎.北京鹫峰公园主要人工林群落多样性研究［J］.西南林学院学报，2008，28（1）：17-22.

ZHAO Yongquan,PENG Daoli.Diversity of the main plantation communities in Jiufeng Forest Park,Beijing［J］.J Southwest For Coll,2008,28（1）：17-22.

［24］ SUN Fengbin,YIN Zhe,LUN Xiaoxiu,et al.Deposition velocity of PM2.5 in the winter and spring above deciduous and coniferous forests in Beijing,China［J］.PLoS One,2014,9（5）：e97723.doi：10.1371/journal.pone.0097723.

［25］ 王昆，张晓丽，王珊，等.鹫峰地区QuickBird影像纹理特征与生物量估测关系初探［J］.地理与地理信息科学，2013，29（3）：52-55.

WANG Kun,ZHANG Xiaoli,WANG Shan,et al.Study on the relationship between texture of QuickBird image and biomass estimation in area of Jiufeng［J］.Geogr Geo-Inf Sci,2013,29（3）：52-55.

［26］ 张晓丽，王书涵，陈江凌，等.应用于快鸟影像的山区单木快速绝对定位及坐标修正方法：CN201310564134.1［P］.2014-02-05.

ZHANG Xiaoli,WANG Shuhan,CHEN Jiangling,et al.Single wooden absolute positioning and coordinate rapid correctio CN201310564134.1［P］.2014-02-05.

［27］ 孟宪宇.测树学［M］.北京：中国林业出版社，2006.

Stand density estimates based on a local maximum spectral filter

WANG Shuhan,ZHANG Xiaoli,ZHU Chenghao,QU Shuai
（Key Laboratory for Silviculture and Conservation of Ministry of Education,Forestry College,Beijing Forestry University,Beijing 100083,China）

This study aims to estimate stand density for different forest types via local maximum （LM）filtering method from high-resolution remote sensing imagery.Stand density was extracted by the LM method to count the number of spectral maximum points extracted from a QuickBird （QB）panchromatic.Research was implemented in the Jiufeng National Forest Park.A high-accuracy digital elevation model（DEM）was used to perform precise ortho-rectification and topographic corrections to correct the images’geometric and spectral distortions.Precise positioning coordinates for the four corner points of a plot were obtained through a combination of differential GPS（DGPS）and Total Station.Spurious tree density calculated within each sample plot was extracted by counting the spectral maximum points with QB imagery.A linear regression model between the true tree density and spurious tree density was established.Spurious stand density was used as the independent variable and stand density was used as the dependent variable.Results showed that the final total correction of the multispectral images was controlled within one pixel at 0.99 Root Mean Square Error（ERMSE）,and the ERMSEof the full-color image correction was 5.86.For a broadleaf forest in Jiufeng National Forest Park,a 5×5 window size and Normalized Difference Vegetation Index（INDVI）≥0.2 achieved the best estimation results（R2= 0.442 2,ERMSE=10.97,P＜0.01）.For the coniferous,broadleaf,and whole area forest models,the coniferousforest had the best results using a 3×3 window size and INDVI≥0.2 （R2=0.741 5,ERMSE=14.45,P＜0.01）.The stand density planning map was also completed using the regression model and the inventory data.The accuracy of stand density estimations of coniferous forest was better than that of broadleaf forest via LM method.［Ch,5 fig.6 tab.27 ref.］

forest measuration;QuickBird image;local maximum filtering;stand density

S758.5；TP701

A

2095-0756（2017）03-0413-08

浙 江 农 林 大 学 学 报，2017，34（3）：413-420

Journal of Zhejiang A＆F University

10.11833/j.issn.2095-0756.2017.03.005

2016-05-30；

2016-07-07

国家高技术研究发展计划（“863”计划）项目（2012AA102001）；教育部北京市森林培育与保护省部共建重点实验室资助项目（2009GJSYS02）

王书涵，博士研究生，从事林业遥感与3S技术研究。E-mail：370926730@qq.com。通信作者：张晓丽，教授，博士后，从事林业遥感与3S技术等研究。E-mail：Zhang-xl@263.net