林木年龄是指林木种子萌发后生长的年数，对林分常以组成林分的林木平均年龄表示林分年龄(宁杨翠, 2012)。通过林龄信息能知道森林所处的生长时期，帮助制定合理的森林经营措施。从森林经营角度讲，林龄是森林经营和森林资源调查工作中的重要因子(胡安泰, 1984)。另外，林龄也被认定为影响森林生态系统净初级生产力(NPP)的主要驱动因子(Bradford et al., 2008)，在森林碳储量、森林生物量及生态系统净生产力(NEP)等相关研究中均扮演着重要角色。

在已有研究中, 确定林龄空间格局特征通常是基于森林样点实测的林龄数据, 建立林龄与其他替代指标的经验关系, 并以此来估算林龄空间分布特征。戴铭等(2011)以第5次全国森林详查获取的省级优势树种的平均林龄及分布面积为基础，以同期生长季节的NOAA/AVHRR的NDVI遥感数据为辅助，在空间降尺度统计技术的支持下，得出了全国8 km分辨率的林龄定量分布。Li等(2014)使用生长期和枯叶期2个物候期的Landsat ETM+数据估算了我国东北地区落叶松(Larix)林的林龄。以6个反射率波段、9个植被指数和6个纹理数据作为原始数据源，分别通过一元线性回归、多元线性回归和多层神经网络模型估算了落叶松林林龄。结果表明，落叶松林林龄估计的最小和最大误差分别为0.47年和21.3年；而白桦(Betula platyphylla)林林龄估计误差在0.6年和10.1年。Zhang等(2014)利用森林资源调查数据拟合“树高-生物量”和“生物量-林龄”之间的函数关系，通过变量替换建立“树高-林龄”的统计函数，再利用全球l km空间分辨率的树高数据(Duncanson et al., 2010)，绘制中国1 km林龄空间分布图。

综上所述，我国关于林龄的研究中存在以下3点不足：首先，传统方法没有考虑观测样点的空间分布，估算结果受林龄空间异质性程度、观测样点数及分布、替代指标的精度与可获取性及经验关系的普适性等多种因素限制。其次，我国大尺度的林龄空间分布数据没考虑干扰影响。以大兴安岭地区为例，1965—2010年46年间发生火灾1 614起，森林总过火林地面积达3 523 011.86 hm²(杜春英等, 2010；胡海清等, 2012a)。仅2001—2010年10年间该区域发生火灾367起，总过火面积487 469.78 hm²(胡海清等, 2012b)。因此，林龄的遥感估算不可忽略林火干扰影响。最后，在确定干扰林龄时，已往的方式是将干扰发生年林龄设置为0，以此类推干扰后的林龄。而现实中，林火干扰后并非所有树木完全死亡，林木死亡情况受干扰强度、树种组成等因素影响。因此，将干扰发生年林龄设为0并不合理，综合考虑干扰强度、树种等因素的影响才能使得干扰林龄估算的更合理。

黑龙江省拥有丰富的林木资源，需发展有效可行的方法来对林龄分布进行评价。由于本研究估算林龄时使用的样地数据没被林火干扰，不能反映林火干扰后的林龄，而研究区域林火过火面积较大，在估测林龄时不能忽略这部分过火区域。因此，本研究首先采用GWR模型估算非干扰林龄；其次，采用dNBR将火烧迹地按未过火、轻度过火、中度过火、重度过火分为4个级别，分别讨论不同林火烈度下不同树种的更替情况，并确定干扰林龄。本研究成果可为进一步研究森林生物量、碳储量等打下基础，同时也可为黑龙江省林区森林经营管理提供科学依据。

1 研究区域概况与研究数据 1.1 研究区域概况

研究区域为黑龙江省林区。黑龙江省位于中国东北部，地理坐标为121°11′E—135°05′E；43°25′N—53°33′N。南北相距约为1 120 km，东西相距约930 km，总面积454 000 km²。黑龙江省地形起伏大，地貌复杂多山，大部分地区海拔超过300 m，最高处达到1 500 m，其地势为东高西低。该地主要为大陆性季风气候，年均气温在0 ℃左右，降水量700 mm左右，且南北、东西差别很大，分布不均。土壤以暗棕壤为主，辅以黑土、黑钙土、白浆土、草甸土。主要树种为兴安落叶松(L. gmelinii)、樟子松(Pinus sylvestris var. mongolica)、红皮云杉(Picea koraiensis)、白桦、蒙古栎(Quercus mongolica)、山杨(Populus davidiana)等。

1.2 研究数据

本研究基础数据包括研究区2010年生长季(6—9月)Landsat-5 TM影像，2000―2010年林火发生区域火前与火后的Landsat-5 TM影像；黑龙江省2010年国家森林资源一类清查数据；黑龙江省2000—2010年火灾记录数据。

火灾记录来自于黑龙江省林业厅，该记录包括起火地点、地理坐标、发现时间、灭火时间、过火面积和起火原因等信息。辅助数据包括由黑龙江省1:5万地形图矢量化得到的10 m分辨率DEM、2000年东北地理与农业生态研究所出版的1:10万东北植被资源分布图。样地与火点分布如图 1所示。

2 研究方法 2.1 非干扰林龄模型构建

1) 提取变量 本研究林龄实测值来自2010年国家森林资源一类清查数据，是林分优势树种的平均年龄。样地共计2 216块样地，覆盖了黑龙江省4个区域：大兴安岭、小兴安岭、三江平原及长白山。随机抽取444块用作模型的精度检验，剩余1 772块参与建模。变量分为林分变量、遥感变量和地形变量。

林分变量来自一类清查数据，包括平均胸径、平均树高、断面积、郁闭度等。遥感变量包括Landsat-5 TM影像的6个波段灰度值；Landsat-5 TM影像6个波段在对角线方向的纹理信息；Landsat-5 TM影像穗帽变换得到的亮度(Brightness)、绿度(Greeness)和湿度(Wetness)波段。地形变量来自黑龙江省1:5万地形图矢量化得到的10 m分辨率DEM，包括海拔、坡度、坡向等。TM影像反射率穗帽变换矩阵系数如表 1所示。

2) 筛选最优变量 通过逐步回归在显著性水平α=0.05下筛选最优变量, 通过检验的变量分别为林分平均高(ASH)、海拔(Altitude)、林分平均胸径(ADBH)、绿度指数(Greeness)与湿度指数(Wetness)(P＜0.001，VIF＜10)。逐步回归结果如表 2所示。ASH和ADBH与林龄存在显著正相关，因生长是年龄的单调非减函数。Greeness、Wetness和林龄存在显著的负相关，这与Wulder等(2004)的研究相符，同时Sivanpillai(2006)研究发现穗帽变换湿度指数与林龄有显著相关关系，表明穗帽变换湿度指数较Landsat ETM+原始波段对林龄有更强的预测能力(Sivanpillai et al., 2006)。所以，ASH、Altitude、ADBH、Greeness、Wetness可作为影响林龄分布的变量进行研究。因变量与自变量的描述性统计如表 3所示。

3) GWR模型 GWR(Geographically weighted regression)模型是由Brunsdon等(1998)基于局域光滑理论提出的一种局部模型，它将样本的地理位置信息参与建模，能够有效地消除残差的空间自相关性。传统的回归分析建立在假设参数估计与数据的地理位置无关的基础上，即所有点的估计皆为无偏估计。GWR模型是传统回归方法的扩展，为了反映局域参数的变化，在传统线性回归统计模型的基础上将空间的影响以距离权重的方式加入到模型中(李明泽 2017)。GWR模型的基本形式如下：

$ Y_{i}=\beta_{0}\left(u_{i}, v_{i}\right)+\sum\limits_{k=1}^{p} \beta_{k}\left(u_{i}, v_{i}\right) x_{i k}+\varepsilon_{i}, i=1, 2, \cdots, n。$

式中：Y_i为响应变量林龄，(u_i, v_i)为第i个样本点的坐标，β_k(u_i, v_i)是第i个样本点上的第k个回归参数，是关于地理位置的函数，ε_i为误差项。

4) 模型残差的检验和评价统计量 (1)残差空间自相关性检验 空间自相关(Spatial autocorrelation)是指一些变量在同一个分布区内的观测数据之间潜在的相互依赖性。Tobler(1970)曾指出“地理学第一定律：任何东西与别的东西之间都是相关的，但近处的东西比远处的东西相关性更强”。传统统计模型的基本假设是要求模型残差之间相互独立，本研究采用空间相关性指数莫兰指数(Moran I)来评价残差的空间自相关性，其表达式如下：

$ I=\frac{n \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} w_{i j}(d)\left(x_{j}-\overline{x}\right)}{\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} w_{i j}(d) \sum\limits_{i}^{n}\left(x_{i}-\overline{x}\right)^{2}}。$

式中：x_i和x_j分别为样点i和j的观测值(i≠j)，x为观测值的平均值，w_ij(d)为根据样地点i和j之间的距离估算的权重。

本研究采用标准正态检验(standard normal test)来判定模型残差的空间分布是否随机。计算GWR模型残差的全局Moran I，Moran I通常假设其值为零时残差分布完全空间随机。如果正态检验得到结果－1.96≤Z≤1.96，则P＞0.05(α=0.05)，表示接受零假设，即模型残差的空间分布模式很可能是随机空间产生的结果，若Z>1.96或者Z < －1.96，则P＜0.05(α=0.05)，则拒绝零假设，即模型残差的空间分布模式为具有统计显著性的聚类或分散模式。

(2) 评价统计量 本研究使用调整决定系数R_adj²和均方根误差RMSE对GWR模型进行拟合效果评价。

5) 精度评价 本研究随机抽取444块样地作为独立检验样本，对模型拟合值与实测值进行对比分析。

6) 方法实现 本研究使用GWR4.0软件建立林龄GWR模型，核函数使用adaptive bi-square，最佳带宽的选择方法为golden section search，其评价标准为AICc；使用Excel宏文件ROOKCASC计算全局Moran I；使用ENVI5.3及ARCGIS10.3软件进行空间可视化分析。

2.2 林火烈度分级与干扰林龄的确定

1) 提取火场边界 首先对2000―2010年火前与火后的TM影像进行几何配准，校正误差严格控制在0.5个像元以内。根据2000―2010年火灾记录表中的经纬度和过火面积记录，在ENVI5.3软件中对TM影像进行火场目视解译，确定火场边界。

2) 林火烈度分级 林火烈度分级在火场边界内进行。林火烈度的估算分级采用dNBR。差值归一化火烧指数(dNBR)由Key等(2006)提出，是火前NBR与火后NBR的差值，研究表明dNBR所利用的多时相技术，相较单时相的NBR能更准确地提取火烧迹地, 李明泽等(2017)采用dNBR对呼中林区火烧迹地林火烈度分级并取得了较好效果。dNBR的取值范围-2~2。为便于分析和处理，采用1 000作为dNBR计算时转化为整数的相乘系数。其取值范围则变为-2 000~2 000。采用Mazuelas(2012)在西班牙火烧迹地制图时制定的dNBR估算火烧烈度的阈值范围(表 4)。NBR与dNBR公式如下：

$ \mathrm{NBR}=\frac{\rho_{\mathrm{nir}}-\rho_{\mathrm{swir}}}{\rho_{\mathrm{nir}}+\rho_{\mathrm{swir}}}, \\ \mathrm{d} \mathrm{NBR}=\mathrm{NBR}_{\mathrm{pre-fire}}-\mathrm{NBR}_{\mathrm{post-fire}}。$

3) 干扰林龄的确定 由于本研究所使用的林龄实测值来自国家森林资源一类清查数据，是林分中优势木的平均年龄，因此确定干扰林龄时，过火区域林龄是否发生改变取决于该地区优势树种是否因林火干扰而发生变化。使用ArcGIS的空间分析模块，将林火烈度分级图与植被资源分布图空间叠加，分析林火烈度与植被类型的相互作用，分别讨论不同火烈度下火烧迹地内各个植被类型的更替情况，并在非干扰林龄的基础上修正干扰森林的林龄。

3 结果与分析 3.1 非干扰林龄拟合结果

1) GWR模型拟合结果与系数 GWR模型能够较为有效地估算黑龙江省林区非干扰林龄，R_adj²和RMSE分别为0.68和16.171 7。GWR模型系数的几种统计量见表 5。

利用反距离权重法(Inverse distance weighted，IDW)对模型中各变量对应参数分别插值得到GWR模型各个变量对应参数的空间分布结果。图 2为GWR模型在预测林龄时各变量对应参数的空间分布，可以看出，海拔对林龄的影响主要是正向的，其中大兴安岭西南部和小兴安岭中部的回归系数最高，呈现向双侧递减的现象，表明在其他条件不变时，较高海拔处具有较高的林龄。在图 2b中，大兴安岭地区的回归系数最大且为正值，小兴安岭地区的回归系数最小。

对整个研究区域而言，林分平均树高对林龄的影响呈现东高西低。在图 2c中，除大兴安岭东南部一部分外其余大部回归系数均为正值，说明林分平均胸径对林龄的影响是正向的，在其他条件不变的情况下，林分平均胸径增大会使林龄增大。

从图 2d中可看出，绿度指数对林龄的影响主要是负向的，大兴安岭地区的回归系数最小且为负，小兴安岭地区回归系数呈现向双侧递增的现象。总体上绿度指数对林龄呈负相关关系，说明绿度指数较小的林分对应着较高的林龄。

在图 2e中，大兴安岭西部、小兴安岭中部、三江平原西部回归系数最小且为负值。大兴安岭地区东部、三江平原东部回归系数最大。综合整个研究区域，湿度指数对林龄的影响是负向的。

2) GWR模型残差的空间自相关性分析 本研究计算了GWR模型的全局Moran I并通过基于随机假设的标准正态检验来评价残差的空间自相关性。国家森林资源一类清查样地的布设间隔为8 km，因此以8 km带宽计算模型的全局Moran I及基于基本假设的标准正态检验结果(Z值)。GWR模型的全局Moran I为-0.013 39，Z值为-0.013 (大于-1.96且小于1.96)，通过了标准正态检验，这说明GWR模型能有效地消除模型残差的空间自相关性。

3) 模型精度验证结果 使用独立样本检验了GWR模型预测效果。GWR模型的林龄模拟和实测结果相关关系见图 3，观测值与预测值R²和RMSE分别为0.683和15.006。

4) 非干扰林龄拟合结果 利用GWR模型定量分析了研究区林龄，其在2010年的非干扰空间分布见图 4。

由图 4可知，黑龙江省森林林龄的变化范围为0~186年，平均林龄48年，标准差为16年，与现地调查结果相符。大兴安岭地区林龄普遍高于研究区平均水平，主要分布在60~80年。对整个研究区域而言，林龄峰值集中在大兴安岭西南部、三江平原东部及小兴安岭东部。小兴安岭地区林龄主要分布在20~60年。大兴安岭东北部地区林龄主要分布在0~20年，为研究区域最小值，主要原因为该地区林火干扰较为频繁，森林多以次生林为主。

3.2 林火烈度分析

1) 黑龙江森林林火烈度分析 通过对比分析火灾前后两景TM影像，目视解译提取火场118处。依据dNBR对2000—2010年黑龙江林区过火区域进行林火烈度分级(图 5)，并计算出各烈度等级的过火面积(表 6)。结果表明，黑龙江林火主要发生在大兴安岭及小兴安岭地区，主要原因是该区森林易燃、气候干旱、地广人稀，对林火控制能力薄弱。总过火面积为527 932 hm²，其中重度火烧29 157 hm²，中度火烧180 268 hm²，轻度火烧318 507 hm²。

2) 林火烈度与植被 通过林火烈度分布图与植被资源分布图的空间叠加，得到不同植被类型上不同林火烈度等级的过火面积(表 7)，可以看出，兴安落叶松林和蒙古栎林在整个研究区中过火面积最大，分别占总面积的28.63%和47.23%。兴安落叶松是大兴安岭地区的代表性植被，约占该地区森林面积的55%(孙家宝 2010)；另外，兴安落叶松林多为单层林，林木稀疏，林下喜光杂草滋生，增加了易燃性。蒙古栎在大兴安岭、小兴安岭地区均有大量分布，因蒙古栎林叶子干燥易燃，不易腐烂，非常容易引起火灾。不同烈度等级的林火在各植被类型上的分布均表现为：轻度火＞中度火＞重度火。

3) 树种更替与干扰林龄估算结果 本研究所讨论的干扰林龄在定义上和非干扰林龄一致，是林分优势树种的林龄。优势树种的改变意味着林龄发生改变，所以有必要讨论各森林类型在不同林火烈度下的树种更替情况，从而确定林火干扰发生后火烧迹地的优势树种是否发生变化。以下是研究区域各森林类型在不同林火烈度下的树种更替情况。

兴安落叶松林：孙家宝(2010)根据大量的火烧迹地调查数据，统计了不同火烈度下兴安落叶松林火后各类物种的重要值。其研究表明，重度林火干扰下，70%活立木被烧死，火烧迹地内仅存留大径级兴安落叶松；火烧1年后，白桦、山杨依靠根茎的萌生能力迅速繁殖占据裸地，成为最先侵入的先锋树种，3年后形成白桦和山杨幼苗群落，10年后白桦、山杨萌生苗迅速生长后开始郁闭，此时白桦的重要值最高，占有绝对优势。中度林火干扰下，30%~50%活立木被烧死，火烧后10年内形成白桦幼苗和兴安落叶松混交林，由于白桦重要值最高，因此认为10年之内白桦为优势树种。轻度林火干扰下，兴安落叶松林受到的影响不大，火烧后数年间仍然处于优势地位。

兴安落叶松、白桦林：根据孙家宝(2010)研究，重度林火干扰后，70%以上林木被破坏，郁闭度降低。由于火烧迹地白桦种源丰富，白桦大量繁殖，火后19年间白桦的重要值都保持在最高水平，具有绝对优势。中度火干扰后，小径级的兴安落叶松被烧死，原有白桦的地上部分被烧死，地下根基部重新萌发，其数量增加，上升为优势树种。轻度林火干扰下，兴安落叶松、白桦林树种组成未受到影响，只有部分小径级的幼树和幼苗被烧死，促进了白桦的更新。

蒙古栎林及蒙古栎、黑桦林：研究(刘广菊，2008)表明，在重度林火干扰下，耐寒、耐旱、萌生能力较强的蒙古栎、黑桦幼树占据迹地。由于立地条件、气候条件的影响，树木生长缓慢，甚至难以恢复，靠天然更新多形成火偏途顶级群落蒙古栎矮曲林或退化草坡。7年内火烧迹地内群落结构和组成变化比较剧烈，7年后火烧迹地内群落结构趋于稳定，蒙古栎逐渐萌发并成为优势树种。中度林火干扰下，火烧迹地内群落结构物种组成变化不大，多形成以蒙古栎为主的偏顶极演替，仍然是蒙古栎、黑桦幼树占据火烧迹地。轻度火烧后，乔木层组成结构没被破坏，植被演替过程只有小乔木、下木和地被物发生变化。

白桦林及白桦、山杨林：研究(刘广菊，2008)表明，山杨更新繁殖能力很强，但树木的抗火、适应能力相对较弱，难以占有优势，白桦实生、萌生能力都很强。在重度林火干扰后演替初期，迹地上白桦幼苗主要依靠萌生呈簇状分布，迅速占据火烧迹地，4~9年后逐渐郁闭，在乔木层成为优势树种。中度火干扰对天然次生林群落草本、灌木影响较大，而对乔木物种数的影响不大，火烧迹地内树种组成比较稳定。轻度火烧干扰下，原有林分中的主林层未受到或轻微受到破坏，故林分基本上保持火烧前水平，林下更新仍以原生地带植被为主，整体上不能改变原有的树种组成。

根据上述讨论，当林火干扰不足以使火烧迹地发生树种更替时，则认为林火干扰对林龄不造成影响。当林火干扰足以使火烧迹地内发生树种更替时，则林龄在林火发生当年记为0，并在新的优势树种萌发时从1开始累加。表 8列出了各森林类型在不同林火烈度下对应的干扰林龄计算方法。图 6为2010年黑龙江森林干扰林龄分布情况。

4 讨论

本研究使用地理加权回归模型估算林龄，与Li等(2014)研究相比建模过程更科学合理，消除了模型残差的空间自相关，更加符合传统统计模型残差间相互独立的基本假设，这是因为地理加权回归模型在模型拟合时加入了样地位置信息，并通过距离衰减函数进行加权。另外，相比前人的研究本研究在林龄估算过程中引入了林火干扰。在确定干扰林龄时，以往的处理方式是将干扰发生年森林林龄设置为0，以此类推干扰后森林的林龄。而现实中发生林火干扰后，森林并非完全死亡，火烧迹地内部林木的死亡情况与演替过程取决于当地的树种组成、生境条件、林木发育期和林火烈度。本研究在确定干扰林龄时推翻了以往研究“干扰即归零”的假设，讨论了火烧迹地的林火烈度和树种组成，细化了定义干扰林龄的过程。在未来的研究中可以结合林木发育期、生境条件，对林火干扰后火烧迹地内的演替过程做出更深入的讨论。由于受到研究数据限制，本研究在估算林龄时仅讨论了林火干扰，未来干扰林龄研究中可增加森林病虫害干扰及人为干扰，使干扰林龄估算更合理。关于火烧迹地内树种演替的分析，本研究仅讨论了2000—2010年短时间内的森林演替情况。近年来随着各种森林生态系统模型的出现，对森林景观长期、大范围的研究变为可能(奚为民，2016)。干扰林龄估算的研究应结合诸如LANDIS(Mladenoff et al., 2004)等人工智能耦合模型，更加真实地还原森林群落演替过程。

5 结论

本研究利用森林资源一类清查数据、Landsat-5 TM数据，使用GWR模型估算了黑龙江省2010年非干扰林龄；利用黑龙江省2000—2010年火灾记录数据、Landsat-5 TM数据，获取了黑龙江省2000—2010年的林火范围，对火烧迹地林火烈度分级，并将林火烈度分布图与植被资源分布图相叠加，分析了不同林火烈度下森林的演替情况，估算了黑龙江省2010年干扰林龄，得出以下结论。

1) GWR模型适合进行林龄拟合，能有效消除模型残差的空间自相关性。黑龙江省非干扰森林平均年龄为48年，标准差为16年，林龄整体空间分布状态不均匀，大兴安岭地区林龄主要分布在60~80年，普遍高于黑龙江林区平均水平。

2) 黑龙江省2000—2010年总过火面积527 932 hm²，其中重度火烧29 157 hm²，中度火烧180 268 hm²，轻度火烧318 507 hm²。林火主要集中在大兴安岭和小兴安岭地区。兴安落叶松林和蒙古栎林过火面积占其总面积的比例最大，分别为28.63%和47.23%，这是因为兴安落叶松以单层林为主，林木稀疏，导致林下喜光杂草滋生，增加了易燃性，而蒙古栎林叶子干燥易燃，不易腐烂，幼林叶子不易脱落，造成火灾周期性发生。

3) 本研究在林龄估算中加入了林火干扰，为林龄估算方法提供了新思路，同时本研究成果也能为黑龙江地区森林NPP、NEP及森林碳储量、生物量等相关研究提供数据支持。