近30年来, 中国是全球森林资源增长最多的国家, 森林覆盖率已由20世纪80年代初的12%提高到目前的23.04%, 森林蓄积量提高到175.6亿立方米, 人工林面积稳居全球第一^[1—2]。从生态系统稳定性来看, 人工林生态系统较天然林更易受到有害生物的攻击^[3], 导致我国森林生物灾害表现为发生种类的多样性, 发生时间与空间的复杂性与不均衡性^[4], 使生物灾害预测评估困难, 尤其缺少定量化指标与方法进行全局分析^[5]。我国森防测报机构近些年一直在推广建立不同类型的生物灾害的发生监测体系^[6—8], 组建了1000多个国家级测报中心对300多种重要林业生物灾害进行持续的监测^[9—10], 对不同种类建立标准化灾害评估方法, 形成了不同程度发生面积的统计报告^[11—12], 但如何利用这些系列发生面积数据产生精准化的全局性比较分析还尚无统一认识。

时间序列数据分析是统计学中一类几乎在所有领域都广泛应用的分析手段。针对数据的平稳性、波动幅度、周期性、随机性等特征, 分析检验方法众多, 经济统计领域尤甚。而气象、灾害等领域时间序列数据, 针对其波动幅度大、非正态数据分布等特性, 目前对于其趋势分析主要方法有曼-肯德尔(Mann Kendall, M-K)趋势检验、创新性趋势分析(ITA)、稳健回归^[13—14]等；对于其突变分析主要方法有M-K突变检验、累计距平法以及滑动T检验^[15]等；对于其稳定性分析常用变异系数^[16]；而对于未来持续性分析则常用Hurst指数^[17]。在不同区域和种类的农林害虫灾害趋势和突变预测中得到应用。

对于我国森林生物灾害的整体趋势, 常采用病情指数、虫情指数^[18]、发生率等来反映有害生物危害的严重程度, 后相继提出了森林病虫指数^[19]、为害指数^[20]、发生指数^[21—22]、森林病虫基指数^[23]、成灾率^[24]等指标, 但均在局部、单种灾害应用或仅是概念的提出, 未能建立可在全局不同尺度时空分析应用的统一系列时间数据集。由于全局性历史灾害数据调查获取与系列时间数据建立的难度, 如何利用现有森林有害生物发生面积统计数据, 建立时间序列数据集, 验证其在全局灾害趋势与突变分析的应用价值, 促进今后的相关持续数据产出和应用是本文研究的出发点。本研究有效利用森林灾害发生面积统计数据, 计算衍生了森林有害生物发生指数(FPOI)、灾害发生率(FPOAR)两个系列指标集。并选用非参数M-K检验、滑动T检验、变异系数、Hurst指数方法, 分析评估我国不同省区森林生物灾害的趋势特征、突变特征、稳定性和持续性特征在时间和空间上的同步与差异, 为未来森林有害生物灾害全局性合理分析评估与持续监测提供了借鉴。

1 数据来源与处理方法 1.1 灾害指数集构建方法

1998—2019年森林有害生物发生面积和森林面积数据源于《中国林业统计年鉴》和《中国林业和草原统计年鉴》^[25—27], 有害生物发生面积包括轻度、中度、重度及总计面积, 在森林有害生物的划分上, 根据林业有害生物发生及成灾标准(LY/T 1681—2006)^[28], 将森林病虫害危害程度划分为轻度、中度、重度三个等级。本文结合危害程度与发生面积, 计算了森林有害生物发生指数(FPOI)、灾害发生率(FPOAR)两个指标集。其中FPOI权衡了有害生物灾害的轻、中、重度发生的不同影响, FPOAR为同期森林生物灾害发生总计面积与森林面积的比值, 不考虑灾害程度的权重。FPOI参考文献^[21]及专利^[22], 将公式中该树种面积替换为森林面积, 以FPOI和FPOAR比对全国森林有害生物在时间层面、空间格局上的发生及变化情况。计算公式如下：

(1)

(2)

1998—2019年全国森林有害生物数据中, 对于异常数据采取删除处理。由于行政区域的划分, 1998—2003年四川的数据中包含重庆, 因此计算1998—2003年四川、重庆FPOI、FPOAR时均除以四川省森林面积；西藏缺失1998—2001、2003—2004年数据；2001—2002年全国合计不包括龙江集团, 2003—2005年不包括龙江集团和新疆兵团。由于本文主要针对省区级有害生物, 因此未将龙江集团、新疆兵团、大兴安岭数据列入分析中。森林资源数据(森林面积)源于第四至九次全国森林资源清查(1989—1993、1994—1998、1999—2003、2004—2008、2009—2013、2014—2018)结果。将单位统一为：公顷。本研究尚缺中国港澳台统计数据。

1.2 分析方法 1.2.1 灾害趋势分析

Sen′s斜率估计与M-K趋势检验已成为时间序列数据变化趋势及显著性趋势检验的重要方法^[29—30], 本文利用Sen′s斜率估计和M-K检验判断森林生物灾害FPOI、FPOAR时间序列的变化趋势, 当Sen′s斜率估计斜率大于0时, 表示森林生物灾害呈上升趋势, 反之呈下降趋势。以M-K方法进行变化显著性检验, 对统计检验量Z, 在0.05置信水平, 若|Z|≥1.96, 表示序列变化趋势显著；在0.01的置信水平, 若|Z|≥2.58时表明序列变化极显著。在RStudio调用trend包采用sens.slope、mk.test函数计算。

1.2.2 灾害突变分析

M-K检验是已广泛使用的非参数的统计检验方法^[31—32], 文中将其用于时段内森林有害生物变化趋势和突变出现时间研究。M-K突变检验通过计算统计量时间序列UF和UB判断所检验时间序列是否存在突变点, 当UF_k和UB_k交点位于临界线内时, 交点对应的时刻便是突变的时间(P < 0.05), 若其超过临界线, 即存在明显的变化趋势。UF_k < 0或UF_k>0表示有害生物发生呈下降或上升趋势。为增强突变结果可信度, 辅以滑动T检验法对突变点进行可信度检验^[33]。在R中编写M-K突变检验代码, 在Python环境下编写滑动T检验代码分别对全国生物灾害逐一检验，以确定具体的突变年份。

1.2.3 灾害稳定性分析

变异系数可用于表示全国各省区森林生物灾害的相对波动程度, 其值越大, 说明该省区有害生物受到的干扰多、不稳定, 反之, 则处于相对稳定状态^[34—35]。其计算公式为：

(3)

式中, CV是变异系数, σ为标准差, x为各省区FPOI、FPOAR的平均值。

1.2.4 灾害未来持续性分析

Hurst指数是预测时间序列数据相对于过去未来发展趋势的一个重要指数。在M-K检验的基础上, 常采用重标极差(R/S)计算序列的Hurst指数, 广泛应用于气象、水文等领域^[36—37]。运用Hurst指数来定量描述FPOI、FPOAR时间序列变化的持续性或反持续性。

Hurst指数H的取值范围为[0, 1], 当0≤H≤0.5时, 具有反向持续性, H越接近0, 反持续性越强；当H=0.5时, 表明独立随机状态；当0.5≤H≤1, 表明该时间序列未来变化趋势具有正向持续性, 即森林有害生物未来变化趋势与过去22年的趋势一致, H值越接近1, 正向持续性越强。根据文献^[38]将Hurst指数定义如下等级：反持续性(H < 0.5), 弱持续性(0.5 < H≤0.65), 强持续性(0.65 < H < 1)。在R中调用pracma包采用hurstexp函数进行计算。

1.2.5 统计学分析

使用Excel 2016整理森林有害生物数据；在R 4.0.2、RStudio中计算森林有害生物的FPOI、FPOAR, 应用trend包进行灾害趋势分析, pracma包进行Hurst指数分析, ggplot2包对数据进行可视化；Python 3.9用于计算变异系数、滑动T检验。

2 结果与分析 2.1 指数集FPOI、FPOAR的构成与特征

本研究计算了1998至2019年森林有害生物FPOI、FPOAR两个系列指标集, 数据集由有害生物发生面积包括轻度、中度、重度及总计面积和森林面积构成, 包含698条数据。全国合计FPOI数据集波动范围：0.24—0.36；FPOAR数据集波动范围：4.82%—6.91%。省区尺度指数波动范围明显大于全国尺度。1998—2019各省年均FPOI、FPOAR按从低到高排列见图 1, 全国尺度FPOI与省区尺度年均FPOI时间序列变化见图 2。

2.2 森林生物灾害历史趋势特征

根据趋势检验结果(表 1), 近22年全国森林生物灾害FPOI和FPOAR分别有21、20个省区以及全国合计通过趋势检验(95%置信区间), 在以下计算中出现的百分比指占省区个数(31个省)的百分比, 下同。

以FPOI指标, 全国尺度呈极显著下降趋势(P < 0.01)。省区尺度森林生物灾害呈显著增加趋势占比9.68%, 主要分布在新疆、西藏、天津等地；呈显著下降的省区为58.06%, 其中极显著下降的省区为48.39%, 集中分布在华北(北京、河北、山西)、东北(黑龙江)、华东(上海、江苏、福建、山东)、中南(河南、广东、海南)、西南(四川、贵州、云南)以及西北(陕西、甘肃、青海、宁夏)地区。近三分之一的省区(安徽、广西、湖北、湖南、吉林、江西、辽宁、内蒙古、浙江、重庆, -1.96≤Z < 1.96)变化趋势不明显。

以FPOAR指标, 全国尺度呈显著上升趋势(P < 0.05)。省区尺度森林生物灾害有19.35%的省区呈极显著上升趋势, 具体表现为天津、湖北、湖南、西藏、新疆、重庆等地区极显著增加；呈极显著下降趋势的省区为41.94%, 分布于华北(北京、山西)、东北(黑龙江)、华东(上海、江苏、山东)、西北(陕西、甘肃、青海、宁夏)以及广东、河南、四川等地；无显著变化趋势省区占38.71%(安徽、福建、广西、贵州、海南、河北、吉林、江西、辽宁、内蒙古、云南、浙江, -1.96≤Z < 1.96)。

2.3 森林生物灾害突变分析

为进一步诊断全国森林生物灾害序列的突变特征, 采用M-K突变检验和滑动T检验法对各省区FPOI、FPOAR序列进行突变性检验, 检验结果见表 2, 给定显著性水平α=0.05。

以FPOI指标, 综合两检验结果共8个省份发生突变。在2003、2007、2012、2017年天津、新疆、西藏、湖北等分别发生了突变；湖南、广西、江西、重庆等地在各自临界线内存在多个交点, 湖南在1998—2014年内有多个交点, 表明1999、2004、2009、2011、2013年是潜在的突变点, 根据滑动T检验显示其在2009年出现突变, 综合两种检验其突变时间为2009年。同理, 在2001年广西发生了突变；江西的突变点为2002年；重庆2014年发生突变；天津、新疆、西藏地区呈上升趋势；广西、湖北、湖南、江西、重庆等地变化趋势不显著。

以FPOAR指标, 全国合计突变点在2001年呈上升趋势, 此外M-K检验显示13个省份发生了突变, 根据滑动T检验辽宁、内蒙古、浙江地区统计量未超出临界线, 不存在突变点, 综合两检验结果共10个省份发生突变。在2007、2009、2012年湖北、湖南、西藏分别发生了突变；2004年为天津、重庆的突变点；2010年为安徽、新疆的突变点；广西、江西、吉林等地在各自临界线内存在多个交点, 借助滑动T检验, 广西的突变点为2001年和2002年；江西和吉林的突变时间均为2014年。湖北、湖南、天津、西藏、新疆、重庆地区表现为上升趋势；安徽、广西、江西、吉林变化不明显。

FPOI、FPOAR突变结果有所差别, 就突变省份而言, FPOI指标表明天津、新疆、西藏具突变年份有害生物呈上升趋势, 而广西、湖北、湖南、江西、重庆等虽具突变年份但无显著性；FPOAR指标显示湖北、湖南、天津、西藏、新疆、重庆地区及全国合计均发生突变灾害呈上升趋势, 安徽、广西、江西、吉林等虽在不同年份具突变点但无显著性。这与M-K趋势检验结果一致。就突变年份而言, 全国合计FPOI未发生突变, 湖北、江西、重庆三个省份FPOI、FPOAR突变点差异较大。

2.4 森林生物灾害稳定性分析

基于变异系数法分析中国森林生物灾害空间格局稳定性(图 3), 可以看出这22年间森林生物灾害整体情况较为稳定, FPOI变异系数的变化范围为0.095—1.55, 均值为0.45±0.3；FPOAR变异系数的变化范围为0.078—1.43, 均值为0.38±0.29。变异系数低于0.5以下FPOI占64.52%, 主要分布在华北、华东、西南、西北等部分地区；FPOAR低于0.5以下的区域占83.87%, 除海南、青海、浙江、重庆、上海地区外其他省区较为稳定。海南、青海地区FPOI、FPOAR变异系数较高, 说明该地区森林有害生物灾害交替发生概率大、不稳定。FPOI、FPOAR变异系数存在差异, FPOAR指标衡量海南、青海、浙江、重庆、上海等5省区表现不稳定, 但以FPOI指标, 增加到除以上省区的甘肃、西藏、山西、河南、江苏、吉林11省区, 其他地区为低值区。

2.5 森林生物灾害未来变化趋势

1998—2019年, 全国合计FPOI的Hurst指数为0.67；全国合计FPOAR的Hurst指数为0.69, 全国FPOI的Hurst指数值域0.52—0.74, 均值为0.67±0.06；FPOAR的Hurst指数值域0.55—0.73, 均值为0.67±0.05, FPOI、FPOAR的Hurst值大于0.5的省区均占比为100%, 说明全国森林生物灾害未来的变化呈正向持续。FPOI、FPOAR的Hurst指数强弱持续性不同, 主要表现在湖南、黑龙江、广西、湖北、浙江、山东、云南等地, FPOI显示前三个省份为弱持续性, 后四个省份表现为强持续性, FPOAR与之相反。

为揭示我国森林有害生物变化趋势与持续性, 将1998—2019年Sen+Mann-Kendall趋势分析与Hurst指数结果叠加, 将不显著增加-强持续、不显著增加-弱持续、不显著减少-强持续、不显著减少-弱持续4个等级定义为不显著区域, 得到变化趋势与持续性的耦合结果(图 4)。1998—2019年间, FPOI各省份变化趋势所占的比例从大到小依次是：显著减少-强持续(45.16%)>不显著(32.26%)>显著减少-弱持续(12.90%)>显著增加-强持续(9.68%)；FPOAR变化趋势比例：不显著(38.71%)>显著减少-强持续(35.48%)>显著增加-强持续(16.13%)>显著减少-弱持续(6.45%)>显著增加-弱持续(3.23%)。

3 讨论

森林生物灾害目前已受到学者的广泛关注^[39—40]。许多研究表明, 森林生物灾害面积居高不下呈高发频发态势, 局地成灾^[41—42], 与以往有关森林生物灾害趋势研究结论有所不同, 研究显示全国有害生物呈下降趋势省区数量明显多于上升省区数量。主要原因有：(1)森林生物灾害面积在不断增长的同时, 同期森林面积也在增长, 森林面积由第四次(1989—1993)全国森林清查的133703500公顷提高到第九次(2014—2018)的220446200公顷, 灾害面积增加相对于森林面积呈下降态势, 而天津、西藏、新疆等显著增加首先是数据本身反映的, 其次可能发生新的病虫害或环境因子所导致；(2)时间序列的差异性, 不同时间序列下的生物灾害发生程度不同；(3)所用方法的差异性, 方法间的原理机制、计算步骤不同, 对灾害趋势的估算结果各异。(4)数据源、调查方式以及所用指标的不同, 会影响灾害统计分析结果, 因此在以后的研究中需结合地域特征使数据标准化、规范化以减小误差。

FPOI、FPOAR指标适用于森林生物灾害发生趋势和突变特征的分析, 其结果略有差异。从森林生物灾害历史趋势角度, FPOAR显示重庆、湖北、湖南呈显著上升趋势, 而FPOI则不显著；FPOI显示福建、贵州、海南、河北、云南呈下降态势, 而FPOAR则不显著, 主要源于两个指标概念的差异, FPOI强调有害生物灾害的轻度、中度、重度发生的影响, FPOAR则不考虑灾害程度的权重, 权重赋值不同也会导致结果不同, 两指标均表明我国森林生物灾害区域化特征明显。从森林生物灾害突变角度, 森林生物灾害是生物和非生物因素共同作用的结果, 环境因素是诱发森林灾害的原因之一^[43—44], 特别是温度、降水、林分结构等因子, 这可能是造成各省区森林生物灾害发生突变的原因。全国尺度上趋势突变发生在2001年, 而在2001年全国只有广西发生突变, 这可能是由于全国1998—2000年有害生物累加结果导致, 根据突变结果突变年份恰好是厄尔尼诺、拉尼娜年或其次年, 厄尔尼诺-南方涛动(ENSO)现象所形成的几种气象灾害都与有害生物的生存、繁殖、蔓延和流行有关^[45]。研究表明ENSO事件与马尾松毛虫、稻飞虱以及萎蔫病等的发生有显著关联^[46—47], 因此如何定量分析这些因素对森林生物灾害的影响等也需进一步深入研究。

森林生物灾害的数据源于历史资料, 由于调查者专业受限这类数据往往会受到质疑^[48]。由于部分省区数据缺失、数据异常等, 导致结果会存在一定误差。另外本文仅31个省市区数据进行如果以地区、县或1000个测报中心的数据进行将更准确, 但所应用的Sen+Manna-Kendall检验、M-K突变检验、变异系数、Hurst指数等研究方法准确把握了全国森林生物灾害的动态变化被证明是切实可行的。在开发沿用以往灾害数据利用基础上, 今后加强原始数据产生的规范性、准确性, 采用数据点更均匀及更密集尺度等前提下, 此指标及方法将可以更直观准确显示我国森林生物灾害整体趋势。以便更为深刻地认识森林有害生物的发生规律。

4 结论

通过森林有害生物FPOI、FPOAR分析了中国22年来森林生物灾害变化趋势与突变特征, 探讨了灾害发生的可能影响因素, 得出以下结论：FPOI、FPOAR指标用于森林生物灾害分析切实可行, 显著程度有所差别, 从灾害变化显示度, 强调灾害权重因子FPOI优于FPOAR；1998—2019年森林生物灾害华北、东北、华东、西北等大部分省区表现为下降态势, 而新疆、西藏、天津等地森林生物灾害呈上升趋势；突变检验中局部区域突变时间节点具有一致性, 而大部分省区并未展现突变；华北、华东、西南、西北等部分地区生物灾害形势较稳定, 海南、青海地区稳定性差；Hurst指数显示未来森林有害生物灾害变化以正向持续为主, 与历史趋势保持一致。