时空过程是随时间进程而发展的对象的空间模式(亦称空间格局)^[1]。过程与模式的耦合并不十分完美, 两个独立的过程可能产生相同的模式。由于直接对过程进行研究往往是非常困难的, 因此常从模式入手来研究过程。作为一个生态过程, 森林害虫的发生必然表现出一定的时空模式。在不同尺度上解析这些模式的成因, 可为研究害虫种群暴发机制提供线索。理论上, 种群的时空模式可能是很复杂的, 如螺旋波、空间混沌和晶格等。西部古毒蛾Orgyia vetusta的种群就呈现一种“图灵斑”的模式^[2], 落叶松小卷蛾的暴发过程中出现周期波的模式^[3]。但在自然系统中, 这些复杂的模式是很少见的。相比之下, 一种常见的模式是种群在宽地理区域内表现出同步性的振荡, 即种群动态的空间同步性, 也称为空间自相关^[4]。研究宽尺度同步性一直是种群生态学的中心问题之一。因为随着地区相关性的增加, Meta种群的全球性留存将下降。害虫管理后, 其种群的恢复力也与种群动态中这种相关性的程度密切相关^[5]。本文主要介绍森林害虫的空间同步性的特点、成因及其研究方法方面的进展。

1 空间同步性及其检测

种群同步性通常指的是不同地点间种群密度(或者丰度)的同步性振荡, 种群的密度既可以是同种的, 也可以是不同种的。空间同步性检测最简单的方法莫过于计算两个空间上不相邻亚种群的时间序列的成对相关性, 然后用整个研究区域的相关性的总平均数作为地区水平上的同步性。Pearson相关系数常用于两个系列之间线性相关性的研究。Haynes^[6]还使用Spearman秩相关来度量同步性。这类度量同步性的方法被称为序列值之间相关性的检测^[7], 在20世纪90年代得到了广泛应用^[8]。但用相关系数来度量种群的同步性并不具有生态学上的意义, 且这类方法的结果不能排除同步性不均质的影响^[9]。

Hanski和Woiwod在1993年提出, 同步性应该指的是变化速率的协同变异^[10], 可用当年与来年丰度比率的时间序列(R_t = N_t/N_t-1或R_t = log N_t - log N_t-1)的相关性来表示^[8]。这类方法被称为变化上关联(如变化上的相关性、方向上变化的一致性)的检测^[7], 其结果最能代表同步性的这个定义, 因为它直接测度两个序列如何一起移动。Buonaccorsi等曾提出一个相对简单的度量方法, A_ij=序列i和j同向移动的次数/(T-1), 而2A_ij-1则为Kendall′s tau的变形^[7]。

对峰值或其他极端行为的相符性检测方法常用于研究害虫暴发在所关注地区是否是经常性的。它需要较长的数据序列, 卡方检验曾用于这类研究^[11-12]。这类方法发展的新技术为相位分析。该方法的度量范围已远超过空间同步性, 可对一些复杂的空间模式进行度量^[13]。

通过模拟的数据研究表明, 以上3类度量有时与其对同步性的评价相符, 有时却不相符^[7]。

空间同步性研究的另一个重要内容是同步性与距离的关系。对此研究的方法可分为参数和非参数两大类。参数检验假定同步性与距离之间存在一个特定的函数关系。由于非参数方法能够处理更加复杂的空间相关模式, 因而得到更加广泛的应用^[8]。Mantel检验曾作为这种研究的标准方法^[4]。该方法是Mantel于1967年提出的^[14], 检验的是样地间距离与相似性之间的总体关系。后来又在此方法上衍生出一些新方法来, 如基于距离的偏相关统计、Mantel相关图、距离矩阵多重回归、以及多距离矩阵间相合性检验(ADM)等^[15]。偏Mantel检验与Mantel相似, 不同之处在于考虑了第3个矩阵的影响(通常是地理距离)^[16]。Mantel相关图先将距离分成类别再检验类别内由成对值所计算出来的自相关程度是否大于或小于样地间总体平均自相关性。然后将标准化的Mantel z值对距离类别作图^[17]。另一种广泛使用的非参数统计方法是NCF法。该方法使用平滑样条来产生空间协方差作为距离延迟函数的直接、连续的估计, 并用bootstrap法产生估计值的置信区间^[18]。Gouveia^[19]将一些第二协方差(如纬度)也考虑在内, 对NCF进行了扩展, 称为偏NCF。

由于Mantel检验的H₀是两个由不相似性构成的矩阵之间没有关系, 而不是两个随机变量或数据表之间是独立的, 因而在生物学上用Mantel检验来检测数据的空间结构或在两个数据集关系研究中控制空间相关性的用法被认为是误用^{[15, 20]}。也就是说Mantel检验并不是对空间自相关的检验。常用的直接描述空间自相关的相关图有Moran′s I^[21]和Geary′s c^[22]。除了相关图外, 常用来描述空间结构的结构函数还有变异函数。变异函数来自于地统计学, 可用于空间模式的检测、建模和估计。变异函数计算的估计值实际上是Geary′s c的未标准化形式^[23]。空间相关分析可得出空间自相关的方向、程度及距离范围, 变异函数可反映出空间自相关性引起的空间异质性的程度和距离范围。如使用地统计学方法对松褐天牛Monochamus alternatus的研究表明, 在松褐天牛种群存在空间自相关的范围为93 m, 其强度随着距离的增加而减小。且空间自相关可解释总空间变异的99.9%^[24]。

以上方法均是针对整个研究区域的空间相关性, 即全局性的相关性, 其结果不能显示研究区域的细节。Anselin^[25]提出了一个将全局空间相关系数分解成为对每个取样单位的局部统计量的方法, 称为LISA。LISA的计算可用Moran′s I或Geary′s c来进行, 分别称为局部Moran′s I或局部Geary′s c。计算结果也可绘成图。

生态学家对空间统计学家和地统计学家的批评似乎并不认同。生态学家认为空间同步性涉及到时间, 是随时间变异的一种特殊的空间自相关, 是对同步性的绝对统计。而简单的自相关指的是一种静态的(非时间上的)依赖, 仅度量相对同步性^[8]。双方似乎也达成了一种“井水不犯河水”的默契, 空间同步性的文章中几乎不提及空间自相关性, 反之亦然。对于多时间点的空间数据, 只能对每个时间点单独分析其半差函数^[26]。

2 森林昆虫种群同步性的特点

空间同步性在昆虫中广泛存在^[27]。为期14年的监测结果表明, 英国的蝴蝶种群动态呈空间同步性^[28]。对美国10个分属于不同目(半翅目、鳞翅目、膜翅目和双翅目)昆虫的研究也证实了这一点^[29]。1938—2002年的数据表明, 森林天幕毛虫Malacosoma disstria在安大略省西北部、安大略省东部/魁北克西部和魁北克西南部3个区域的动态呈同步性^[30]。在法国, 松带蛾Thaumetopoea pityocampa在气候相似的邻近地区同步发生^[31]。此外, 空间同步性还在黄杉毒蛾Orgyia pseudotsugata^[32]、山毛榉舟蛾Quadricalcarifera punctatella^[33]、云杉色卷蛾Choristoneura fumiferana, 落叶松芽小卷蛾Zeiraphera diniana、山松小蠹Dendroctonus ponderosae和舞毒蛾Lymantria dispar^[9]等昆虫中发现。

不同昆虫种类的空间同步性不一样。在英国, 蚜虫的空间同步性要高于蛾类^[34]。

空间同步性与距离有关, 常随距离的增大而下降。在北极挪威北部的Troms县的小岛上, 冬尺蠖Operophtera brumata在其主要栖息地成熟白桦林内表现出同步的动态, 然而在不同小岛的样地间则表现出最大程度的不同^[35]。尽管空间同步性随距离的增大而下降, 空间同步性的范围可达1000 km^[8]。在英国, 蚜虫和蛾类的空间同步性达到了800 km^[34]。在北美东部, 云杉小卷蛾暴发的空间同步性随距离的增加而下降, 当距离接近2000 km时, 空间同步性趋于0^[36]。当然也有空间同步性先下降而后随距离上升的现象, 这取决于所采取的尺度^[37]。

空间同步性还与时间尺度有关。例如, 在1993年前后, 英国的蚜虫在长时间尺度上(>4年)的同步性下降, 而短时间尺度(≤4年)上的同步性却增加了^[38]。

3 森林害虫发生空间同步性的机制及研究方法 3.1 森林害虫发生空间同步性的机制

已有的理论研究表明地区随机性和扩散是解释空间同步性的较合理的两种主要机制。引起地区随机性最明显的因素是气象因素^[39]。如果受到地区化密度依赖因素(通常是气象因素)的相似影响, 空间上隔离的具有完全相同密度依赖动态的两个种群将被同步化。这就是由Moran^[40]提出, 并由Royama^[41]完善的Moran效应。扩散对空间同步性的影响是通过循环种群锁相实现的。当然扩散也会通过空间诱导的混沌引起空间异步性, 这取决于扩散的模式^[5]。这两种机制可单独发挥作用, 但通常通过非加性模式联合发挥作用。仅当环境相关性和扩散速率都很小时, 这两种机制的交互效应才较小。由于交互效应的符号与环境相关性相反, 因而在通常情况下当样地之间的环境相关性为正时, 通过简单相加得到的种群同步性比实际情况要小。且随着种群管理强度的增加, 种群同步性下降。因此在解释种群同步性的模式时, 必须将扩散和环境相关性联合考虑^[5]。第三种被广泛认可的空间同步机制是与其他物种(如天敌等)的生态作用。但这些营养级上的相互作用者很难对一些昆虫的大尺度同步性产生显著影响, 因为这些天敌(如动物、拟寄生物、病原微生物等)的活动能力基本上是相当于或小于昆虫自己的活动力^[9]。

目前大多数的研究结果表明, 相比扩散而言, 气象因素对空间同步性的影响要大。如英国, 在夜蛾和蚜虫中的结果就是这样^[34]。在新西兰, 舞毒蛾暴发的同步性也是由Moran效应所驱动的^[42]。对来自北半球(包括亚洲、北美洲和欧洲)的26种森林鳞翅目昆虫的140次暴发的资料分析表明, 气象因素, 特别是温度和降水与这些害虫的暴发相关。而种群暴发的周期及凉爽的天气则与太阳黑子的低谷相关^[43]。Peltonen对加拿大5种森林昆虫在景观尺度上的发生数据进行了研究。结果表明, 空间同步性并不直接与扩散能力相联系。相比之下, 气象因素的空间相关性可解释发生数据的同步性水平。因而对这些昆虫来说, Moran效应可能是影响空间动态的决定性过程^[9]。在中部阿巴拉契亚山(Appalachian)的落叶林中, 10种蝴蝶的空间同步性由气象因子驱动, 或至少部分由气象因子驱动^[44]。根据来自阿尔卑斯山脉拱处65个地点的3113根落叶松树轮宽度序列的分析表明, 落叶松小卷蛾Zeiraphera diniana的发生与两个世纪以来的气温变化相一致^[45]。在阿拉斯加的中南部和西南部, 气象因素(温度、降水等)是驱动红翅大小蠹Dendroctonus rufipennis暴发的主要因素^[46]。在英国, NAO(North Atlantic Oscillation)指数与云杉举蚜Elatobium abietinum的空间同步性相关^[47]。在PNW(Pacific Northwest)地区的中心地带, 重建的西部云杉卷蛾Choristoneura occidentalis三个世纪的发生数据表明, 其发生与干旱历史相吻合^[48]。对横跨全英的20种蚜虫的研究表明, 在长时间尺度上(>4年), 80%的同步性是由冬季的气候驱动的^[38]。草地螟Loxostege sticticalis是一种迁飞性昆虫, 其在中国、俄罗斯亚洲部分和蒙古国的大发生具有高度同步性, 高空气流可能是驱动这一同步性的主要原因^[49]。

除了具有迁飞性的昆虫之外, 扩散引起的同步性通常是短距离的(几公里)^[50]。但这并不意味着扩散的作用可以忽略。在北美东部, 云杉小卷蛾暴发的空间同步性就是由Moran效应和高的扩散率共同作用的结果^[36]。

但已有的研究表明, 有些害虫的发生不能仅由扩散和Moran效应来解释。两个相隔种群间即使是很少个体的交流也能导致同步性^[51]。当地其他的环境条件变异的同步性可能会掩盖气象因子驱动力的同步性效应。如舞毒蛾种群同步性下降的速度要快于气象因子同步性下降的速度^[9]。从1967—1998年的记录数据分析来看, 气象因素的空间自相关性几乎不影响云杉色卷蛾C. fumiferana在加拿大安大略省的空间分布及发生频率^[52]。还有, 森林昆虫暴发的空间同步性受到森林管理措施的影响。云杉色卷蛾C. fumiferana在不同强度营林措施(采伐)下的发生动态不能用不同取样地点的气候来解释。在保护区(几乎未采伐)其空间同步性要大于精细管理(采伐强度很大)的区域^[53]。栖息地特征(树木密度、物候上的不同步性、农药使用情况、地貌指数、非寄主植物的丰度等)可解释舞毒蛾所造成的失叶变异的21%—34%。在去掉栖息地特征的影响后, 剩余失叶变异的空间自相关性存在的距离范围为788 m—461 m, 而舞毒蛾幼虫的扩散范围为100 m—1 km^[54]。

3.2 空间同步性成因的研究方法

一个空间模式的成因来自两个方面, 一是其自身的空间结构, 称为真正的空间自相关；二是来自于具有空间自相关的外部因素(如环境因素), 称为诱导的空间依赖。通常情况下, 这两个方面是很难区分的。对森林害虫的暴发而言, 真正的空间自相关主要由扩散等因素引起, 而诱导的空间依赖主要由Moran效应引起。对森林保护工作来说, 对真正的空间自相关是无能为力的。森林保护工作者唯一能做的是了解诱导的空间依赖, 根据环境因素的变化, 制定相关的管理措施。

由于空间相关性的天然存在, 使研究者获取的数据不是随机的和相互独立的, 不满足一些经典的统计方法所要求的误差独立的假设。因而必须采用另外的统计方法。为了区分天气因素、空间邻近和森林类型对舞毒蛾暴发空间同步性的影响, Haynes^[6]采用了MRM(multiple regression on resemblance(or distance) matrices)法来研究一个因变量距离或相似矩阵与其他解释变量距离或相似矩阵的关系, 然后将空间邻近矩阵和其他环境因子相似矩阵都作为解释变量, 将因变量矩阵的方差分割为环境来源和空间来源(如扩散)。Foster^[54]先使用一个简单的线性模型和等级分割来研究寄主丰度、物候异步性、地貌特征及杀虫剂使用情况对舞毒蛾密度的影响, 然后利用变异函数研究残差的空间自相关性。

使用混合效应模型也是应对数据不独立的有效方法^[55]。随机效应的加入, 好比给不同组的数据赋予了一个不同的“基线”, 从而解决数据不独立的问题。可在GLM(generalized linear model)、GLS(generalized least squares)、LME(linear mixed-effects model)、GLMM(generalized linear mixed model)、GAMM(generalized additive mixed model)或GEE(generalized estimating equation)等模型中加入以变异函数表示的空间相关结构。Zuur在Mixed Effects Models and Extensions in Ecology with R一书对这些方法的应用有详细的讲解和举例^[56]。Dormann于2007年对这些方法进行了回顾和模拟研究^[57]。从Ⅰ类Ⅱ类错误率来看, 这些模型对正态数据和泊松数据的拟合都不错, 但对于二项分布的数据要差一些。在对这些方法进行选择时, 响应变量中误差分布的类型应是一个重要的标准。对于正态分布, 基于GLS的方法(GLS、SAR(simultaneous autoregressive models)、CAR(conditional autoregressive models))非常有效。对具有不同误差分布的空间自相关来说, 最灵活的方法是空间GLMM和GEE。尽管在作者的模拟研究中对系数的估计并不好, 自协变的方法被认为足以区别有意义的解释变量和随机协变量。自回归模型(SAR和CAR)需要对权重矩阵制定一个赋值的方案, 因而在研究中通常不是首选。GEE需要很高的存储能力来解未聚类GEE得分方程, 因此对巨大样本而言, 在固定模型中的应用就非常有限。但它特别适合于缺失值和非栅格数据。在实践中一种好的策略是使用不同的方法进行拟合, 然后使用类似于AIC(Akaike information criterion)的标准对这些模型进行比较。

近年来, 越来越多的研究将空间邻近做为影响种群结构的一些无法测量的过程(如扩散)的解释变量^[55]。在这些研究中多采用MEM(Moran′s eigenvector map)做为空间变量。MEM是一些来自于一个以一定距离截取的空间距离网络的特征向量。这些特征向量由样地的空间坐标来计算, 样地位置的排列可规则, 也可不规则。MEM的特点是可以在所有尺度上模拟空间变异, 而且这些MEM之间是正交的。在研究中, 这些MEM被分成一些子模型来对应不同的研究尺度。由于这些MEM之间是正交的, 子模型之间是线性独立的, 没有可共同解释的变异, 因而可用于线性模型, 包括变异分解^[58]。MEM用于去趋的数据, 对未去趋数据, 可用AME法。MEM建模也可用于时间序列的分析^[59], 只不过用于时间序列分析时, MEM是来自于以一定时间间隔(或延迟(lag))截取的时间距离网络的特征向量。

变异分解是将害虫发生的总变异分解成几个独立的组分, 一般包括4个组分：一个纯粹的环境组分(a), 一个由于邻近效应引起的纯粹的空间组分(b), 一个共同的组分(c)以及一个未解释的变异(d)。然后使用一个线性的模型就可以将总变异进行分解(V_t=∑a+b+c+d)^[60]。计算的结果通常以Venn图的形式来表现, 这样就可以直观地看出哪个组分是害虫发生的主要动驱动力。

4 森林虫害发生数据的获取

森林虫害发生的数据, 通常有3个来源。其一是利用观测数据或历史记录数据^[61]。其二是利用遥感或空中调查数据^[54]。其三是利用树轮学的原理进行重建^[62]。第一个来源需要观测者或者记录者具有一定的专业知识。从国内的情况来看, 基层林场的森林保护工作者的专业水平有待进一步提高, 而且从业者往往不固定, 因而这类数据经常受到诟病。第二类数据一方面受到年限的限制, 不易获得较长的历史数据, 另一方面图片经常受到云层等因素的干扰。第三类数据是利用害虫暴发影响树木径生长的科学理论进行历史数据的重建, 一方面, “树木不撒谎”, 因而这类数据的可靠性非常高；另一方面, 可获得较长的历史数据(几百年, 甚至上千年, 取决于树的年龄)。由于这两方面的优点, 这类数据已广泛用于考古、气象和生态学研究。利用树轮序列重建害虫发生的历史, 已是成熟的方法, 只是在中文文献中鲜有报道。由于在我国, 森林有害生物的管理工作中将有害生物作为灾害来管理, 而是否成灾的标准通常是是否对树木的生长造成影响。因而利用第三类数据进行森林生物灾害管理的话, 亦可避开争论多年的关于害虫“发生”与“成灾”的话题。当然也有研究者担心用生长锥取样会对树木造成伤害。事实上伤害是存在的。但物理性的伤害是几乎可以忽略的, 而较明显的后果主要是引起木材腐朽菌的感染, 导致伤口周围变色。总体上伤害是非常小的, 且对生长锥进行消毒或堵塞取样孔对减小伤害的效果并不显著^[63]。

重建的理论基础是生物灾害影响树木的径生长。采用树轮学的方法, 找出径生长受干扰的年份(年轮显著变窄的年份), 然后排除影响树木径生长的其他因素如自身的生长趋势、火灾、气象因素等的影响。利用树轮学原理重建森林害虫发生历史的方法从取样到构建地区年表阶段, 与在考古、气象和生态学方面的方法基本一样。不同之处在于害虫暴发年份的确定。通常采用RWI(ring-width index)法来确定害虫暴发年份。首先是RWI的计算, 公式为：RWI_t= R_t/G_t

式中, RWI_t是t年的年轮宽度指数, R_t为t年的实际年轮宽度, G_t为该样芯在t年的去趋势后的年轮宽度。RWI法又分非寄主RWI法和地区RWI法。这两种方法的目的都是去掉气象因素对年轮生长的影响, 得到害虫对年轮生长的影响。非寄主RWI法的理论基础是害虫对寄主和非寄主的危害不同, 因而其年轮的径生长不同。但要去掉气象因素的影响, 所选取的非寄主与寄主在对气候变化的响应(径生长)上应该相同或相似。地区RWI法的理论基础是害虫的暴发总是在一个足够大的地区内的一些小地区中。因此, 地区范围的划定是否合适将直接影响到判定结果的对错。气象因素之外, 火灾也是一个影响林木径生长的主要因素之一。非寄主RWI法不能排除火灾的影响, 地区RWI法在理论上是可以的。火灾的发生情况可通过调查样地内的火疤来确定或询问林场的经营者。而在国内, 绝大多数的森林都是天然次生林和人工林, 林龄都不大, 所有的火灾应该都有记录。当然没有一种方法是完美的, 研究结果的得出还需结合森林保护专业的知识。

非寄主RWI法多采用Swetnam^[64]的方法计算一个校正的指数, 即

式中, I代表RWI, c代表校正的意思, t意为第t年, h意为寄主, n意为非寄主, SD意为标准差, I意为平均RWI。

然后按一定的标准来判定暴发年。如Ryerson^[65]采用以下标准：(1)将I_c标准化后, 只少有8个负值；(2)发生年的前后各有1个正值；(3)发生区间内只少有1年的I_c低于-1.28SD。

而地区RWI法的判定标准通常是当年样芯的RWI小于1, 且比地区的RWI低于1.28SD判定为受干扰年。地区的RWI是该地区内所有样芯的平均值。1.28SD也是一个常用的标准, 在研究中可根据具体情况进行调整。Paritsis^[66]在发生年份的最终判定上并不只依赖于RWI值, 而是依据以下3个标准：(1)从年轮上看, 发生年的年轮明显窄于相邻的年份；(2)Pointer Years法判定为受干扰年；(3)RWI法判定为害虫发生年。对样芯而言, 满足以上3个条件的发生年为高置信度发生年, 满足其中2个条件的称为中置信度发生年, 满足第2和第3中的一个的称为低置信度发生年。对样地和地区而言, 同时满足第2和第3的发生年称为高置信度发生年, 只满足其中一个的发生年, 称为低置信度发生年。

5 小结

我国是森林生物灾害多发的国家, 而研究森林害虫的空间模式是揭示害虫成灾机制的关键内容之一。在诸多空间模式中, 同步性是常见的模式。引起同步性的主要因素是扩散和Moran效应, 即真实的空间自相关和诱导的空间依赖。然而这两种情况通常很难区分。本文综述了森林害虫同步性方面的研究进展以及研究空间同步性的方法, 特别是利用树轮学原理重建森林害虫发生历史以及近年来兴起的将dbMEM作为空间变量来解释因变量变异的方法, 以期推动我国森林害虫成灾机制方面的研究。