-
森林结构参数的准确获取是森林资源规划调查的重要工作之一[1-4]。林分密度与林分的树冠、胸径、树高的生长有显著的相关性[5-6]，是森林资源二类调查中一项重要的参数和指标[7-9]。同时，林分密度还与人工林的木材材性、生物量、蓄积量等有着密切的关系[10-11]，林分密度是表征森林生态结构，影响生态系统特征，以及林业可持续发展的关键因子。森林测计学中林分密度可分为疏密度、株数密度和郁闭度。本研究所讨论的林分密度仅指株数密度，为单位面积株数/单位样地面积。通常林分密度的获取方式主要包括样地测量和遥感手段估测[12]。常规的测量林分密度的方式一般是采用样地调查方法以布设样地为主，需要耗费较多的人力物力且不能实现大区域的估测，长期复查容易导致误差和重复测量[13-14]。遥感技术的推广应用给地区尺度进行林分密度的估测提供了有力的工具[15]。随着高分辨率遥感应用越来越广，探讨高分遥感提取森林植被参数就具有重要意义。基于高分遥感提取林分密度的方法有：二维各向同性高斯核函数平滑滤波[8]、分水岭方法[16]、发射或辐射的反演模型方法[17]、局部变化模型方法[18]和光谱最大值滤波法[19-20]。光谱最大值滤波是一种可以用来检测单木位置的技术，它是基于针叶树冠的顶点会存在较大遥感影像像元亮度值（DN值）的假设[21]，优势是能快速有效地探测到林地单木的位置[22]。因此，光谱最大值法可以用来间接估算林分密度。光谱最大值法的普适性需要进一步研究，特别是随着地域不同，林分类型不同时究竟采用多大的窗口以及应该采用何种方法精炼树冠点都是值得深入研究的问题，而本研究从此出发探究了不同的窗口大小和归一化植被指数（INDVI）阈值的光谱最大值法应用于不同的林分类型提取林分密度的精度问题。
HTML
-
鹫峰森林公园坐落于北京市海淀区西北部苏家坨镇境内，是距离京城最近的国家森林公园之一。鹫峰经纬度大约为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]。
-
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）进行了正射校正，覆盖了主要的研究区域。以匹配地面实测样地的位置，保证建模的精度。
-
研究小组收集了具有代表性的72个标准样地，样地大小为20 m × 20 m，样地基本上均匀分布在研究区域中。调查时间分别是2006年和2011年，调查方法采用全站仪和手持差分全球定位系统（GPS）相互配合的方式精确定位样地的4个角的位置具体流程可参见样地定位方法[26]，并且记录下每株树的树种、胸径、树高、冠幅等信息。根据实际调查样地的树种组成，针叶树种大于0.65的划分为针叶林，阔叶树种比例大于0.65的划分为阔叶林[27]，最终划分为25个针叶林样地，43个阔叶林样地，样地的冠幅描述性统计信息参见表 1和表 2。
平均冠幅/m 冠幅最小值/m 冠幅最大值/m 南北冠幅 东西冠幅 南北冠幅 东西冠幅 南北冠幅 东西冠幅 2.67 2.96 1.27 1.50 6.74 6.48 Table 1. Crown size statistiques in coniferous sample plot
平均冠幅/m 冠幅最小值/m 冠幅最大值/m 南北冠幅 东西冠幅 南北冠幅 东西冠幅 南北冠幅 东西冠幅 3.35 3.56 1.52 1.87 6.92 6.59 Table 2. Crown size statistics in broadleaved sample plot
1.1. 研究区域
1.2. 数据获取与处理
1.3. 外业调查数据
-
首先将QuickBird全色波段和多光谱波段分别进行预处理操作，首先进行基于控制点的正射校正，由于本研究获取了鹫峰区域的1:2 000大比例尺的地形图，由北京市测绘部门提供，地形图已经栅格化且经过了精校正，因此用地形图直接校正快鸟影像。DEM的获取采用如下方式：将等高线的线图层在ArcGIS地形分析模块中将其转换成不规则三角网（TIN），再将TIN转换成1 m精度的DEM。以二类调查数据为依据先确定林地区域，在QuickBird影像上裁剪出研究区。
-
当一个像素的直径（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语言中完成。
Figure 1. Flow Chart
2.1. 图像预处理
2.2. 窗口选择与INDVI阈值的选择
-
随着滤波窗口大小的变化，图像逐渐变得模糊，结果如图 2所示。图A′，图 2B′，图 2C′为利用全色图像减去滤波后图像的效果示意图。由于截取的区域土壤和树冠的差异较为明显，可以从图 2中比较清晰地分辨出树冠以及道路的边缘，有些较大的单木树冠清晰可辨。图 3以B11样地为例显示了将样地准确叠加在经过处理过后的光谱最大值点的图层上，可以清楚地统计出样地范围内所包含的单木数量。
Figure 2. Sketch map of local maximum filtering using different window size
Figure 3. Schematic map of B11 sample plot and local maximum points
-
根据研究方法中的介绍，分别选取了不同的窗口大小及INDVI阈值提取林分密度，其中根据分析处理结果见表 3~表 5。
所有林分 INDVI≥0.1 INDVI≥0.2 INDVI≥0.3 INDVI≥0.4 INDVI≥0.5 3×3窗口 0.523 4*** 0.521 5*** 0.523 1*** 0.485 3*** 0.435 1*** 5×5窗口 0.543 8*** 0.545 5*** 0.539 6*** 0.483 0*** 0.421 2*** 7×7窗口 0.475 1*** 0.482 9*** 0.484 6*** 0.427 7*** 0.342 8*** 说明：***代表模型在0.001水平上差异显著。 Table 3. Correlation analysis of all stand using the based on the local maximum of different INDVI and windows sizes (all stand)
针叶林 INDVI≥0.1 INDVI≥0.2 INDVI≥0.3 INDVI≥0.4 INDVI≥0.5 3×3窗口 0.740 2*** 0.741 5*** 0.743 0*** 0.718 5*** 0.656 3*** 5×5窗口 0.651 6*** 0.653 1*** 0.650 7*** 0.617 6*** 0.535 1*** 7×7窗口 0.645 1*** 0.649 8*** 0.644 8*** 0.612 2*** 0.525 4*** 说明：***代表模型在0.001水平上差异显著。 Table 4. Correlation analysis of all stand using the based on the local maximum of different INDVI and windows sizes (coniferous stand)
阔叶林 INDVI≥0.1 INDVI≥0.2 INDVI≥0.3 INDVI≥0.4 INDVI≥0.5 3×3窗口 0.366 1*** 0.364 2*** 0.338 7*** 0.140 6** 0.204 5** 5×5窗口 0.439 6*** 0.442 2*** 0.400 7** 0.165 5** 0.213 1** 7×7窗口 0.381 6*** 0.333 9*** 0.375 5** 0.120 5* 0.073 6* 说明：***代表模型在0.001水平上差异显著，**代表在0.01水平上差异显著，*代表在0.05水平上差异显著。 Table 5. Correlation analysis of all stand using the based on the local maximum of different INDVI and windows sizes (broadleaved stand)
由表 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的窗口大小获得的决定系数最高。由此可见，窗口大小是影响林分密度的估计最主要的因素。
-
根据上述分析，分别选取了较高的相关系数所对应的INDVI阈值和窗口大小，建立所有林分、阔叶林、针叶林的统计模型。统计结果如表 6所示。
INDVI 窗口大小选择 林分类型 决定系数 P值 F分布统计量 均方根误差 自由度 多0.2 5×5窗口 阔叶林 0.442 2 6.285E-07*** 34.08 10.97 1, 43 多0.2 3×3窗口 针叶林 0.741 5 1.648E-08*** 68.83 14.45 1, 24 多0.2 5×5窗口 所有林分 0.545 5 8.758E-14*** 85.23 13.97 1, 71 说明：***表示拟合结果在0.001水平上差异显著。 Table 6. Correlation analysis in different forest types
图 4显示了针叶林样地和阔叶林样地最佳INDVI阈值和最佳窗口选择下的林分密度提取散点图，纵坐标代表了实际林分密度，横坐标代表了基于光谱局部最大值滤波提取出的林分密度，从而建立线性回归模型。以针叶林和阔叶林范围生成的格网为单位，将模型应用于针叶林和阔叶林的光谱最大值点数量图层的栅格图，从而得到相应的林分密度分布图，其中针叶林和阔叶林的范围可以根据二类调查小班数据来确定。分析结果见图 5，其中红色区域代表相应林分的林分密度为0。
Figure 4. Scatter diagram between local maximum and real tree density in different forest type
Figure 5. Stand density estimation result
3.1. 窗口大小和INDVI阈值的选择
3.2. 林分密度估测模型及估测结果
-
本研究利用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）。研究结果进一步表明了窗口大小的选择对光谱最大值法提取林分密度有重要的影响，本研究分林分类型进一步提高了林分密度的估测精度。
-
本研究通过利用光谱最大值滤波方法应用于高空间分辨率遥感影像QuickBird试图实现数字化提取研究区内的林分密度流程，取得了较好的结果，其中针叶林的林分密度提取效果最好，建立的针叶林林分密度提取模型可以用来估测研究区域的林分密度，为进一步实现大尺度林分密度提供基础，而阔叶林的估测精度稍低，但模型仍然有利用价值。光谱最大值滤波方法在提取树冠中心点是假设每个树冠只有一个光谱反射值最大点为前提来推测的，因此，可以推测由于针叶林树冠具有规则的树冠形状，通常针叶树种的树冠往往只有一个光谱最大值点，而阔叶树由于具有较大面积的树冠和复杂的树冠结构往往不止一个树冠反射率最大点，这是造成针叶林提取精度高而阔叶林提取精度偏低的主要原因，而我们的研究结果跟前人的研究结论是一致的[19, 21]。注意到本研究得出结论阔叶林的林分密度选用5×5的窗口效果最优，这刚好对应于3 m × 3 m的树冠范围，这与表 2中阔叶林中样地的平均树冠大小较为接近，而在估测针叶林林分密度时选用3×3的窗口大小效果最优，也跟表 1中针叶林中样地的平均树冠大小一致，这说明窗口大小的选择应该要贴近于树冠的真实大小才能获得较高的精度。
本研究仍然有几点值得探讨：① 利用0.61 m分辨率的全色波段提取的光谱最大值点能够反映出亚米级的树冠位置，其结果跟影像的太阳高度角、太阳天顶角、卫星高度角、卫星天顶角等参数关系密切，因此，相同类型的传感器在不同时相得到的影像采用相同的方法结果是否有偏差还需进一步验证，因此可以利用不同时相的影像进行实验，分析验证本方案的可行性，进一步修正本研究方法。这是未来的研究内容之一。② 由于近红外波段与植被的关系密切，还应该充分挖掘近红外波段的潜力，例如采用近红外波段提取光谱局部最大值是否会得到更好的结果还有待研究。③ 本研究所提出的窗口大小的选择和INDVI过滤值的设定方案是否具有普适性还有待验证，研究相同的研究区不同的数据源是否需要设置不同的INDVI阈值和不同的窗口大小，及可变的窗口大小进行光谱最大值滤波进而探测树冠的位置应该是今后另一个提高精度的研究方向。
-
在论文完成之际，向北京林业大学林学院国家高技术研究发展计划（“863”计划）项目组所有人员在野外样地调查中所付出的努力表示衷心感谢！感谢相关部门提供的数据支持！
DownLoad: