基于GIS的夷陵区松材线虫入侵松林时空分布特征分析

周 丹，肖正利，赵 勇，田静波，李祖梅，敬顺华

（1.夷陵区自然资源和规划局，湖北 宜昌 443100；2.远安县林业局，湖北 宜昌 444200；3.夷陵区森林植物检疫站，湖北 宜昌 443100；4.夷陵区分乡镇林业管理站，湖北 宜昌 443100；5.夷陵区黄花镇林业管理站，湖北 宜昌 443100）

松材线虫病，又称松材线虫萎蔫病、松树萎蔫病、松树枯萎病［1］，是国际重大检疫性病害，是目前最为严重的自然灾害之一，该病害在热带和亚热带地域迅速蔓延造成严重危害［2］，也是危害较大的外来入侵物种之一［3］，松树一旦被松材线虫感染，基本无法幸存，被感染的松林成片枯死，对松林生态系统造成了严重的破坏。松材线虫有关的研究涉及从分类到防治、病原、媒介、基因、行为以及共生菌等诸多领域，总体上没有太大突破［4］。

松林线虫病是由滑刃科（Aphelenchoididae）伞滑刃属（Bursaphelenchus）的松材线虫（Bursaphelenchus xylophilus）引起的松树上的一种毁灭性病害［1，5］；在全国范围内，几乎所有气候带的绝大多数松林生态系统都有可能被松材线虫侵染［2］，其传播方式主要是自然传播和人为传播，自然传播是借助媒介昆虫天牛携带来完成，可作为媒介的昆虫有13种［1］，国内发现能携带松材线虫的昆虫有6种［6］，松材线虫病的媒介昆虫除松褐天牛外，国内外已经报道了林间存在其他昆虫可携带松材线虫而扩散传播松材线虫病的现象［7］；该病的传播扩散趋势还未得到有效遏制。

松材线虫病疫情传播和扩散具有一定的地理空间位置分布和发生时间规律，并具有关联性，分析松材线虫病地理空间位置数据信息，其目的在于准确了解疫情发生的空间分布状况，挖掘其空间规律；量化松材线虫病疫情地理空间过程并构建相应的时空模型，为科学监测、预测、调查、评估等发挥辅助决策，通过地理信息系统分析研究，对受害松林密度、分布以及变化进行了时空模拟，助力松材线虫病疫情防控提供辅助科技支撑。

1 材料与方法

1.1 研究区域概况

夷陵区位于湖北省长江中、上游，地理位置处于110.512 258 33°—111.392 230 33°E，30.322 233 33°—31.282 230 33°N，举世瞩目的三峡大坝座落在夷陵区境内；地处大巴山余脉，主体地貌以中低山为主，地势由东南向西北逐渐升高，地形起伏多变，错综复杂，高度相差悬殊，海拔51～2 005 m，平均海拔650 m，从南到北兼用中亚热带、北亚热带和暖温带的气候特点，气候适宜，自然条件优越，日照充足，年平均气温16.6℃，极端最高气温41.4℃，最低气温-12℃，≥10℃的年积温5 408℃，无霜期278 d，年平均日照时数为1 669.2 h；年平均相对湿度78.5%。全区平均年降水量在997～1 370 mm。降雨季节和降雨量极不均匀，5～9月雨量集中，降雨量在800 mm；自然气候温暖湿润，适合外来入侵生物松材线虫的生长与繁殖。

该区域森林面积24.838×104hm2，现有森林植被以天然次生植被为主，植物种类繁多，森林植被类型属于亚热带落叶阔叶林，以落叶阔叶与常绿针叶混交林等多种类型的植物结构为主，同时兼有中亚热带常绿阔叶林和北亚热带常绿阔叶林，森林覆盖率74.87%，地段性森林植被类型主体为落叶阔叶混交林，以壳斗科（Fagaceae）物种组成植被为主体森林群落居多。松科植物马尾松占27.90%，马尾松林呈纯林和混交分布；阔叶混交林占45.33%，其他林占26.77%，马尾松在该区域森林植被类型上占有一定的生态空间与分布范围。

松材线虫病于2006年首次在小溪塔街办梅子垭村松林检测出［8］，2009年成为疫区［9］。该区域人口相对密集，交通发达、经济活动频繁，松木在不同区域间的调运频繁，大规模流动为松材线虫疫情的传播和扩散创造了条件。松材线虫入侵马尾松林致病力强，死亡速度快、数量大、扩散快，治理难度大，给该区域造成了严重的经济和生态损失，该病给国民经济造成重大损失，破坏了自然景观和生态环境，对松林资源构成严重威胁。

1.2 资料数据来源

资料数据来源于夷陵区2019—2020年2个年度松材线虫病枯死木采伐管理小班数据库（*.SHP格式），基于每一株病枯木空间位置，以二清为基础底图建立的数据库，用于外业调查、疫木（枯死木）采伐管理、监测、检查验收，ArcGIS软件平台运用空间连接（Spatial Join）工具将病枯木位置与株数汇集成面图层数据库。

1.3 研究方法

为了便于数据分析，松材线虫病枯死木采伐管理小班数据库点面图层之间相互转换采用点模式分析，点模式分析可以用来描述任何类型的事件数据［10］，每个小班（点）数据视为一个疫点（或侵染源、疫源）。由于松树感染松材线虫导致死亡具有发生空间位置和发生时间，导致松树死亡时间段主要集中在秋季，难以获取地面上每个疫点（小班）每一个月（周）时间段疫情数据，时间分析采用年度比较法；运用数据描述性分析、趋势分析；时空聚类分析采用空间自相关、最近邻距离指数、核密度分析［10-14］；空间信息图制作与计算采用ArcGIS 10.8软件分析、统计。

1.3.1 描述性分析 描述松材线虫病疫情发生基本情况、比较时间与空间趋势演变，用数据量化松材线虫病疫情空间传播分布特征，绘制疫情图（包括松林危害程度空间分布、疫点分布、疫点三维趋势）。

1.3.2 趋势面分析 采用线性模型与数学曲面拟合样本数据，建立二元多项式回归模型统计，运用趋势分析工具（Trend analysis）将疫点松林受害株数值投影到一个由东西向和南北向平面构成的正交平面上，投影平面上拟合曲线形成趋势图，显示疫情在空间的总体分布规律，以此反映区域的系统变异，建立趋势面模型。

1.3.3 时空聚类分析

1）空间自相关（Global Moran’s I）。探索性空间数据分析（ESDA）是检验是否存在一个统计学意义上显著的空间分布，并进一步了解生成这一分布空间过程的一种方法［15］，空间自相关是ESDA的重要方法，空间自相关分为全局和局部自相关，可揭示空间变量的区域结构形态，是检验具有某变量在特定空间位置上的属性值是否与其相邻空间位置上的属性值显著相关［13-16］，相邻对象属性值相似，则为正相关；相邻对象属性值相异，则为负相关；相邻对象属性值随机，则相邻对象无相关性。

对松材线虫疫点空间传播的空间自相关性进行实证分析，可描述疫点空间整体分布状况，检验疫点是否与其相邻空间点上的属性值相关联的重要指标，其结果正相关表明属性值变化与其相邻空间单元具有相同的变化趋势，代表空间现象有集聚性的存在；负相关则为随机现象。当一个区域具有高的受害松林密度，且和周围区域为正相关，则该区域就可以认为是个疫点。

该分析采用全局自相关、空间自相关指标莫兰指数（Moran’sI）表示，是从整体上反映研究区域的空间相关性，其取值范围［-1，1］，当I＞0时，表示正相关，且越接近于1聚集性越大；当I＜0，负相关，I＝0，是随机现象［10，11］。局部空间自相关是在全局分析的基础上，没有聚集出现时，用来检查是否有小的、可能被忽视的聚集，可反映子区域间的相关程度。全局空间自相关的Moran’sI计算公式为：

式中，zi为疫点i的属性与其平均值偏差，wi,j为疫点i和j之间的空间权重，n为疫点的总数，So为所有空间权重的聚合。E[I]、V[I]分别表示Moran’sI的期望值和方差，当p很小时，所观测到的空间要素属性不太可能产生于随机过程（小概率事件），因此可以拒绝零假设，z得分和p值都与标准正态分布相关联。

2）最近邻距离指数（NNI）。最近邻距离指数是测量每个要素的质心与其最近邻要素的质心之间的距离，计算最近邻距离的平均值，用最邻近点之间的距离描述空间分布模式，是判断点要素空间分布特征和空间格局特征的空间计量方法，将疫点作为要素点分析，可计算范围区域内任意一疫点与最邻近疫点的距离d min，取平均值得到实际最近距离dˉmin，与随机分布模式下的理论最邻近距离E（d min）比较，得到最邻近指数R，如果R=1，说明事件分布模式属于随机分布；R＜1，说明事件分布模式属于集聚分布；R＞1，则表现的模式趋向于扩散。计算见式（6）。

式中，d i为疫点i与最邻近疫点之间的距离；A为区域面积（采用默认面积，数据最小外接矩形决定）；n为疫点个数；d min为最邻近距离；dˉmin平均距离。

3）核密度分析（KDE）。核密度分析工具（Kernel density）用于计算疫点受害松林值在其周围邻域中的密度，以样本疫点为圆心，搜索圆半径，离圆心越近，栅格单元的密度值就越高，能直观地表现研究对象的分布密度，核密度值的高低代表研究对象在空间分布上的集聚程度［16］。该方法可计算任意疫点区域半径范围内受害松林值的密度，可反映受害松林空间位置分布的相对集中程度，可以据此判断受害松林整体空间分布特点。计算见式（7）。

预测疫点位置（x，y）范围内的受害松林株数值密度，i=1，…，n，是输入疫点，位于（x，y）位置的半径距离内，包括总和中的疫点；popi是i点的领域数量（population）字段值参数，该字段参数为全部小班内受害死亡松值，disti是疫点i和（x，y）位置之间的距离。

2 结果与分析

2.1 基本情况

该区域松材线虫病疫情2009年经国家林业局鉴定确认疫区以来，疫情扩散蔓延迅速，10年来累计死亡松树208.97×104株，发生面积达7.49×104hm2，自然传播媒介为松褐天牛成虫，一年发生1代，羽化周期长，躲藏在松树皮夹缝或树干中，难以防治；主要危害松科植物马尾松、华山松，给该区域造成经济与生态损失2.1亿元。由表1分析可知，2020年度，发生范围扩散到11个乡镇108个村，传播途径为自然传播与人为活动传播；受害松林面积1.126×104hm2，发生疫点（小班）个数7 360个，清理与采伐枯死木36.2×104株，小班受害程度1～5 000株，处置方式为现场焚烧，目前，森林成灾率＜3‰，危害程度处于轻度；该区域疫区、疫点（小班）较多，存在较高的自然传播扩散风险。2020年度危害程度与2019年同期相比，在防治措施下，发生形势和危害程度减轻，但受害松林值基数仍然较大，存在加重或扩散趋势，当前，应继续加大资金投入，采取了防治措施，遏制了危害程度。

表1 松材线虫疫情基础数据

2.2 可视化分析

运用ArcGIS软件处理和输出进行可视化，对疫情图层分析，以自然间断点分级法，将受害松林程度值从低到高的顺序分成10个等级，用由浅入深的颜色和点的大小以地图可视化形式来表示，得到受害松林程度分布图和疫点分布，详见图1至图4。松材线虫病致松林受害程度高发区域分布在分乡镇、黄花镇、雅雀岭、龙泉镇、邓村乡，其受害松林平面分布特点聚集性和集中连片暴发，疫点较多，受害株数值程度较高，分布面存在一定的区域，属于传播高风险区。

图1 2019年度松林危害程度空间分布

图4 2020年度疫点分布

传播时间上，2019年度发生疫点（小班）5 542个，危害程度区间值1～4 917株；2020年度发生疫点（小班）7 360个，危害程度区间值1～2 209株，在防控措施下，受害松林值有所下降，但降幅小；传播空间上，随着时间的推移，2020年度疫情的传播范围向外围发生了扩散，由中部向西部和北部区域发生了扩散，在防控措施下，疫点受害程度从高度聚集分布趋向稀疏分布，且分布较为分散。2个年度疫情整体时空特征为疫点（小班）较多，分布存在一定的区域，难以防控到全部疫点，疫情短期内难以根除。

图2 2020年度松林危害程度空间分布

2.3 趋势面分析

以疫点地理位置坐标（x，y），松林受害株数值为自变量，利用空间探索数据分析，将疫点进行趋势面分析，三维数据点（X，Y，Z），X和Y代表各区域几何中心处的经纬度方向，X轴为方向东，Y轴为方向北，Z代表自变量松林受害株数值（i）；疫点杆的高低代表受害株数值，i分别在平面XY和YZ上投影散点图，通过散点做二项式拟合，得到松林受害株数值在经纬度方向上的某种曲线变化趋势（图5、图6）。结果表明，2019—2020 2个年度松线虫病疫情趋势为东西部，东部较高，2020年向中部进行了扩散，从投影趋势线值分析，趋势线近似于直线，总体较平、低，疫点呈现锥形向中心点聚集性扩散。

图3 2019年度疫点分布

图5 2019年度疫点三维趋势

2.4 疫点空间聚集特征分析

根据疫点的分布和受害松林株数值密度分布特征，通过多次设定不同距离阈值，发现距离半径（R）达到5 km时效果最为理想，生成的密度栅格越平滑且概化程度越高，反映松林株数值密度的相对集中程度。搜集每个疫点之间距离阈值半径5 km受害松林株数，得到受害松林株数发生区域密度分布（图7、图8），不同层次颜色栅格分别代表受害松林株数密度变化规律、区域分级、扩散趋势、空间分布聚集变化特征。核密度显示，2019年度聚集程度以290株为基数8个自然间断点分级，聚集程度区间0～2 321株；2020年度聚集程度以152株为基数8个自然间断点分级，聚集程度区间0～1 361株；2个年度疫区大体相同，2019年发生形势和危害程度偏重；2020年疫情范围出现了扩张趋势，偏重区域发生了转移；松林受害高密集区分布在乐天溪镇、分乡镇、黄花镇；密度图类似于热力图，由中心点高度聚集，向四周扩散，与最邻近距离指数结果一致，松林受害特征为集聚分布；密度分布图可作为松材线虫入侵疫情危害程度分区的依据。

图7 2019年受害松林密度空间分布

图8 2020年受害松林密度空间分布

2.5 全局自相关分析

由表2可知，2个年度时间段所有疫点样本分别计算，所涉距离阀值分别取值1与2 km，搜索疫点周围半径1与2 km范围疫点松林值，检验空间位置上的疫点属性值是否与其相邻空间位置上的疫点属性值相关，结果显示，2个年度所选距离阀值莫兰指数（Moran’s I）取值范围大于0，P值（概率）为0，则随机产生此聚类模式的可能性小于1%。区域范围内疫点存在正相关关系，整体格局上呈现空间聚集性，即疫源传播致松林受害株数高度聚集性。

表2 疫点空间自相关指数

2.6 最邻近距离指数分析

空间自相关分析仅能分析数据的空间聚集性，但无法揭示空间分布格局；空间分布分为三种类型：聚集的、均匀的和随机。均匀分布的点状要素最邻近距离最大，随机分布次之，聚集分布最小；最邻近距离是表示点状事物在地理空间中相互邻近程度的地理指标，最邻近点指数能很好地反映点状要素的空间分布特征；运用平均最近邻分析工具（Average Nearest Neighbor），根据每个疫点预期最近疫点的平均距离计算最近邻指数，预期平均距离是假设随机分布模式中领域间的平均距离；所涉面积以系统默认面积、疫点最小外接矩形来决定，计算结果见表3，2个年度样本疫点分别计算，疫点的平均观测值为173.34和169.82 m，小于预期平均值433.4和348.9 m，最近邻比率指数为0.399 9和0.486 7（0～1），R＜1，对应Z检验值得分-85.46和-84.24，P值为0，即随机分布概率为0，表明分布模式呈现聚集分布。

表3 疫点最近邻指数

3 小结与讨论

3.1 小结

1）空间可视化分析了该区域松材线虫近2年疫情数据，通过时空对比，2019年度松树枯死木达到历年峰值，在防治措施下，2020年受害松林呈现下降趋势，疫点（小班）下降幅度偏小，总体受害松林株数值较高；其空间分布特点为聚集性和集中连片暴发，分布面较大，疫点较多，存在一定的区域性，具有传播扩散高风险，自然传播扩散暴发将长期存在；松材线虫病疫情传播空间上，随着时间发展推移，疫情的传播范围向外围发生扩散，疫点受害程度从高度聚集分布趋向稀疏分布，新发疫点呈跳跃式发生，且分布较为分散，疫情难以防控到全部疫点，疫点短期内难以根除。

2）三维趋势面直观展示松材线虫疫情在三维尺度上的变化趋势，2个年度松材线虫病疫情趋势为东西部，东部较高，2020年中部发生了扩散，从投影趋势线值分析，趋势线近似直线，总体较平、低，疫点呈现锥形向中心点聚集性扩散。

3）时空聚类结果。①2个年度疫点样本空间自相关取距离阀值半径1与2 km，莫兰指数值大于0，为正相关关系，整体格局呈现空间聚集性，疫源传播致松树受害株数值高度聚集。②最近邻距离指数平均观测值为173.34和169.82 m，小于预期平均值433.4和348.9 m，最邻近比率指数为0.399 9和0.486 7，对应检验值得分为-85.46和-84.24，P为0，疫点的分布模式为聚集分布。

4）疫点空间聚集特征核密度显示，距离阈值半径（R）5 km 2个年度受害松林株数以8个自然间断点分级，聚集程度区间为0～2 321株和0～1 361株，2019年危害程度较为偏重，2020年出现了扩张趋势，2个年度发生疫区大体相同，2020年度偏重区域发生了转移；密度图类似于热力图，由中心点高度聚集，向四周扩散，与最邻近距离指数结果一致，松林受害特征为聚集分布；密度分布图可作为松材线虫入侵疫情危害程度分区的依据。

3.2 讨论

松材线虫病疫情扩散、蔓延发生的因子很多，本研究存在一定局限性，疫情的传播因素与媒介昆虫的密度、海拔、气候、松林的分布、松林的密度、环境因子、人为活动干扰等相关，以及历年疫点数据不完善，需要叠加研究，将松材线虫病疫情进行空间管理和空间分析是一种科学有效的管理方法。

1）从疫情2个年度发生时间分析，在采取防治措施下，松林受害株数总体呈现下降趋势，但分布范围发生了扩散，疫情整体发展较快，局部区域疫情呈加快扩散趋势，局部受害松树在增加，存在集中暴发趋势，跳跃式分布；该区域疫情分布广、疫点分散，一旦到达一定疫点数量和分布区域，在短期内难以根除。

2）从疫情2个年度发生空间分析，大部分区域分布集中，局部零星分布，新发疫情呈跳跃式分布，疫情分布广、疫点分散，空间分布与空间格局呈现高度聚集性；从疫情空间发展趋势上看，疫情整体快速发展，部分区域疫情呈加快扩散趋势。新区域受害松树逐年增加，受害面积不断扩大，部分区域存在集中暴发趋势。认为松材线虫疫点应由最初离散分布到聚类分布发展，是一个长期过程，达到聚集性，松林危害程度极重。

3）松材线虫病防治措施，目前国内还没有更好的防治方法，地面活体松林无法判断是否感染；松材线虫不会主动感染松树，需要松褐天牛携带才能传播入侵健康松树。综合防治松材线虫病的传播媒介松褐天牛，减少传播媒介虫口密度，空间上阻挡疫点，通过逐年减少疫点，加大疫木的清理力度、及早发现和清除受感染的树木、阻止疫情扩散，调节松树密度，建立阔叶林生态系统。