利用大面积无人机航拍影像研究松材线虫病病死树林间分布规律

胡可炎,邓创创,陆雪雷,谢伟龙,刘春燕,孙 思*

(1. 华南农业大学林学与风景园林学院,广州 510642;2. 河源市林业局,广东河源 517001;3. 广东省森林资源保育中心,广州 510130)

松材线虫病是由松材线虫Bursaphelenchusxylophilus引起的一种可导致多种松属植物迅速死亡的森林病害(叶建仁, 2019)。目前松材线虫病在全球多个国家都有分布,包括美国、中国、日本、韩国、葡萄牙、希腊、加拿大、墨西哥、西班牙等(Fernándezetal., 2006; Kimetal., 2020)。1982年中国首次出现松材线虫入侵的报道(孙永春, 1982)。截止2022年3月,全国共有19个省(区、市)731个县级行政区发生松材线虫病疫情(国家林业和草原局, 2022)。受该病害影响的松树多达16种,且大多对松材线虫高度或中度感病,而马尾松Pinusmassoniana是南方松林的常见树种,分布于南方15个省区,面积达2×107ha,因此对该病害的防治任务刻不容缓(叶建仁, 2019)。广东省自1988年首次在深圳发现松材线虫病后,一直都是我国松材线虫病危害较为严重的地区(黄焕华, 1990)。

针对松材线虫病治愈难特点,及时发现病死树并清理显得尤为重要,因此病害防控的重点环节之一就是监测。目前,监测的手段主要是地面调查、无人机遥感和卫星遥感3种(Wulderetal., 2006a, 2006b)。由于我国松林面积广阔,而且多数松林分布于陡峭的山区(李浩等, 2020),因此近年来,利用各种遥感技术进行森林资源保护和监测已经成为趋势(胡根生等, 2013)。就航天遥感技术而言,由于受分辨率、卫星周期和天气情况等因素的影响,难以准确及时地发现染病松树和病死树(李浩等, 2020)。随着近几年无人机技术、数字摄影数字处理平台和图像处理软件的快速发展,利用高分辨率航拍影像进行森林病虫害监测成为研究热点(Lehmannetal., 1998; Park &Kim, 2014; Näsietal., 2015),并逐渐发展成为一种低能耗、低成本、高效环保的新型森林病虫害监测方式(郭庆华等, 2016; 黄华毅等, 2021; 李凤迪等, 2021)。

目前利用无人机航拍调查病死松树的研究主要集中于在航拍影像上快速准确地识别和定位病死松树等内容,均属于调查监测技术本身,尚无松材线虫病病死树林间分布规律与空间分布格局的报道。而要进行病死树木分布空间规律的分析,必须要有大面积、足够数量的具有代表性的林地样本,因此本研究选择了固定翼无人机进行航拍作业(黄焕华等, 2018),旨在探索固定翼无人机大面积航拍的工作效率,同时又利用大面积的林地样本探究不同地形因子对松材线虫病病死树的分布规律。为以后松材线虫病的预报监测提供基础数据;也可以结合疫区的地形因子,提高今后对疫情的防治效率(马涛等, 2018; 温阿敏等, 2018)。

1 材料与方法

1.1 试验地概况

河源市(23°10′～24°47′N,114°14′～115°36′E)位于广东省东北部,地形以山地、丘陵为主,山地约占53%,丘陵约占36%,平原和谷地占11%。河源市丘陵、山地海拔较低,全市坡度大都在30°以下。河源市属亚热带季风气候,气温偏高, 年平均气温20.7℃, 年平均降水量为1 768.9 mm,雨量较为丰沛。试验地位于和平县阳明镇(YM)、紫金县紫城镇(ZC)、东源县义合镇(YH)3地,是主要以马尾松为优势种的针阔混交林(图1),面积约为3 500 ha。

1.2 无人机数据获取与影像处理

航拍无人机为北京韦加无人机公司的翔宇Ⅰ型固定翼电动无人机(表1),配备尼康D810航拍专用相机(表2)进行拍摄。拍摄时飞行高度为480 m,飞行速度为14.6 m/s,照片重叠度为航向75%,旁向65%,未携带差分GPS。共飞行了 3个架次,飞行时间为9∶30-15∶30,当天完成了外业数据采集工作。

利用Pix4Dmapper对航拍影像图进行拼接,生成DOM(Digital Ortho Model)数据,分辨率为12.5 cm/pix。接下来用ENVI(The Environment for Visualizing Images)软件对有多余的图像进行剪裁或者对图像进行纠正,导入矢量边界,剪裁出需要的栅格数据,再利用eCognition(易康)软件的规则树对图像进行分割,分类和树冠提取,筛选出病死树信息,并对病死树进行标记。

试验采用的是多尺度分割和光谱差异分割结合,分割尺度参数为30,shape为0.1,compactness为0.7。由于试验研究主体为病死树,故一定程度加大紧致度和颜色的权重。影像分类采用监督分类的方法,将颜色作为主要分类特征。由于软件功能和算法的局限,经过软件初步处理过的图像仍需要通过人工对软件无法识别的病死树进行标记。

图1 试验地航拍影像局部Fig.1 Part of aerial image of the experimental site

表1 无人机参数

1.3 样地设置及调查

在松材线虫病疫情严重的范围内设置连续样方,阳明镇36个,紫城镇30个,义合镇26个,共92个样方,每个样方的面积为2 500 m2(500 m×500 m)。为避免冬季其他阔叶树颜色干扰,调查时间设在松材线虫病死树高峰期9-10月。依据地形的划分标准,对各个地形要素进行分级,具体见表3,分级的标准依据国家林业和草原局颁布的《森林资源规划设计调查主要技术规定》(2003年)和《森林生态学》(李俊清,2006)。通过实地调查,确定样地及航拍的区域,由于疫区较大,还需要包括不同的立地因子,故选择500 m×500 m相邻地块作为样地进行研究。

表2 相机参数

表3 地形要素分级与赋值

1.4 病死树数量及地形因子统计方法

1.4.1病死树数量统计方法

影像分类完成后可通过人工目视结合样地调查时的了解对分割结果进行调整,然后将每棵病死树转为矢量点,再导入ArcGIS进行松材线虫病病死树数量统计。

1.4.2方位与坡向数据获取

在ArcGIS中分别加载广东省河源市和平县阳明镇、紫金县紫城镇、东源县义合镇的的矢量行政边界和DEM数据(均由河源市林业局提供,下同),选择地理处理工具[ArcToolbox]中,找到[3D Analyst]—[栅格表面]—[栅格表面],指定方位和坡向参数,执行命令后生成坡向栅格图(图2-I、II)。

1.4.3坡度数据获取

在ArcGIS中分别加载广东省河源市和平县阳明镇、紫金县紫城镇、东源县义合镇的矢量行政边界和DEM数据。坡度可以利用地表的曲率进行描述和量化,直线形和凸型斜坡在曲率上的对应为曲率大于或等于0,凹型坡和阶梯型坡在曲率上的对应为曲率小于0。在ArcGIS中,选择spatial analysis中的surface-curvature,提取曲率,最终获得图像(图2-III)。

1.4.4海拔数据获取

在91卫图助手上下载广东省河源市和平县阳明镇、紫金县紫城镇、东源县义合镇的高程图,导入ArcGIS。然后在[ArcToolbox]中,执行命令[3DAnalyst工具]—[栅格表面]—[等值线],按地形划分标准指定各参数,执行后生成等高线矢量图层,后对等高线图层进行详细标注(图2-IV)。

1.5 松材线虫病病死树分布格局计算及检验方法

种群分布格局采用以下指标来进行测定:扩散系数(C)、负二项指数(K)、丛生指标(I)、Cassie指标(CA)、聚块性指标(m*/m)和ArcGIS中的“平均最近邻空间统计”工具。

图2 地形因子获取Fig.2 Terrain factor acquisition注:I,病死树方位分布图;II,病死树坡向分布图;III,病死树坡度分布图;IV,病死树海拔分布图;A,阳明镇;B,紫城镇;C,义合镇。Note: I, Azimuth distribution of dead trees; II, Slope aspect distribution of dead trees; III, Slope position distribution of dead trees; IV, Altitude distribution of dead trees; A, YM; B, ZC; C, YH.

(1)扩散系数(C)

C=S2/x

式中:S2为松材线虫病病死树出现的频次的方差,x是松材线虫病病死出现频次的平均数。当C=1时,松材线虫病病死树分布格局为随机分布;C<1时,分布格局为均匀分布;C>1时,分布格局为聚集分布。

(2)负二项指数(K)

K=x2/(S2-x)

式中:S2为松材线虫病病死树出现的频次的方差,x是松材线虫病病死出现频次的平均数。K值越小,聚集度越大,反之,K值越大,聚集度越小,当K值达到无穷大时,则为随机分布。

(3)Cassie指标(CA)

CA=1/K

式中:K为负二项指数。当CA=0,分布格局为随机分布。当CA<0时,分布格局为均匀分布。当CA>0时,分布格局为聚集分布。

(4)丛生指标(I)

I=S2/x-1

式中:S2为松材线虫病病死树出现的频次的方差,x是松材线虫病病死出现频次的平均数。当I=0时,分布格局为随机分布。当I<0时,分布格局为均匀分布。当I>0时,分布格局为聚集分布。

(5)聚块性指标(m*/m)

m*/m= 1+1/K

式中:K为负二项指数。当m*/m=1时,分布格局为随机分布。当m*/m<1时,分布格局为均匀分布。当m*/m>1时,分布格局为聚集分布。

(6)抽样方法和数量的确定

在每个标准地内,用双对角线取样法、平行线取样法、五点法、“Z”字取样法,抽样,通过与3地的病死率进行比较,得到最佳的抽样方法。

双对角线取样法:双对角线取样法是在地块两条对角线上较为平均地分配调查样点进行取样。

平行线取样法:在样地中每隔数行或者一定距离进行调查。

五点取样法:随机选择五个点取样。

“Z”字取样法:取样的样点分布于样地边,中间较少。

1.6 数据统计及分析方法

数据采用SPSS 22.0软件进行分析,重复数据计算其平均数和标准差,并采用Duncan新复极差法比较数据之间的差异性。在P=0.05的水平下进行差异性检测。

2 结果与分析

2.1 松材线虫病病死树空间分布格局

计算得出和平县阳明镇、紫金县紫城镇、东源县义合镇3地松材线虫病病死树分布格局各项判定参数,详细数据见表4。经过各项分布格局判定指标,3个标准地皆为聚集分布。

利用“平均最近邻空间统计”工具分析YM、ZC、YH这3个样地病死树的分布规律,结果z值分别为-16.029755,-11.237614,-8.537191,均代表聚集分布。3份报告中的图形结果十分接近(图3)。

2.2 抽样方法的比较

在阳明镇、紫城镇及义合镇样地随机选择5块样地,分别采用五点法、双对角线法、平行线法、“Z”字法4种不同抽样方法进行样地抽样统计,分别计算4种抽样方法所得数量与总体平均数的标准差(见表5)。结果显示双对角线法、平行线法、“Z”字法所得平均数与总体平均数均有明显差异,五点法与总体平均值没有明显差异,故在调查病死树数量时可选择五点法进行抽样。

表4 松材线虫病病死树分布格局

图3 枯死松树“平均最近邻空间统计”工具分析报告Fig.3 Analysis report of“average nearest neighbor” tool for dead pines

2.3 广东省河源市病死树分布总体特征

经软件分析和人工补充标记,河源市3处试验地病死树共计5 794株。其中各试验地由软件分析得到的病死树均占本试验地病死树总数的60%～70%,需人工补充标记30%～40%。

3处试验地中,和平县阳明镇有病死树3 710株,占64.03%;紫金县紫城镇有病死树1 563株,占26.98%;东源县义合镇有病死树701株,占12.10%。病死树主要分布在西坡、南坡和东南坡,西坡最多为25.94%,其次是南坡23.57%(图4)。病死树主要分布在半阳坡和阳坡,半阳坡占36.54%,阳坡占34.09%(图5)。病死树主要分布海拔区间在300～350 m和250～300 m之间,其中300～350 m的病死树占30.43%,250～300 m的病死树占21.83%(图6)。

2.4 不同地形因子病死树分布规律

2.4.1方位

经统计病死树发生轻微林分(病死率小于1%,下同)共有66块样地,选取其中24块进行分析。分析后发现,病死树发生轻微林分中病死树主要分布在南坡和西坡,分别为总病死树数量的24.40%与20.37%,单位面积样地中同样是南坡病死树数量最多,为5.7株/ha。

表5 四种不同抽样方式病死树数量统计表

图4 不同方位病死树数量统计图Fig.4 Statistics on the number of dead trees in different azimuth

图5 不同坡向病死树数量统计图Fig.5 Statistical on the number of dead trees in different slope aspects

图6 不同海拔分级病死树数量统计图Fig.6 Statistical on the number of dead trees in different altitude

经统计病死树发生中等林分(病死率1%～3%,下同)共有24块样地。在发生中等林分中病死树主要分布也在南坡和西坡,分别为33.01%与15.98%,单位面积中南坡病死树数量最多,为11.25株/ha,位于南坡的病死树数量与除西坡外地其他方位病死树数量具有显著差异(P<0.05)(图7)。

图7 不同方位病死树数量分析图Fig.7 Analysis of the number of dead trees in different azimuth

病死树发生严重林分(病死率大于3%,下同)只有和平县阳明镇的9号和18号样地分别达到了5.80%和5.81%。经分析,病死树主要分布于西坡和东南坡,西坡病死树最多占44.67%,东南坡29.65%次之。单位面积下东南坡病死树数量20.4株/ha最多,略多于西坡。

2.4.2坡向

病死树发生轻微林分中,病死树主要分布在阳坡和半阳坡,其中阳坡数量最多为31.06%,半阳坡次之,占30.18%。单位面积样地中也是阳坡病死树数量3.45株/ha最多。

病死树发生中等林分中,病死树主要分布在半阳坡和阳坡,半阳坡病死树数量最多,达到36.90%,单位面积下同样是半阳坡病死树数量11.7株/ha最多。

病死树发生严重林分中,病死树主要分布于半阳坡和阳坡,半阳坡病死树与阳坡病死树数量之比约为1∶1,阴坡和半阴坡没有病死树分布。单位面积下半阳坡病死树数量22.61株/ha远多于阳坡。

2.4.3坡度

病死树发生轻微林分中,位于凸坡的病死树数量要略高于凹坡,占总数51.72%。凸坡与凹坡的面积比约为1.38∶1,单位面积下位于凹坡的病死树数量2.1株/ha要高于凸坡病死树数量。

病死树发生中等林分中,位于凸坡的病死树要远超位于凹坡的病死树,占总数62.95%,凸坡与凹坡的面积比约为1.31∶1,通过将病死树数量与面积相比较,发现中等林分中病死树发生位于凸坡的病死树数量4.8株/ha高于凹坡的病死树数量。

病死树发生严重林分中,松材线虫病病死树主要分布在凹坡,共约占总数76.75%,远超过凸坡的病死树数量。凸坡与凹坡的面积比约为1∶2.47,单位面积样地下的病死树分布在凹坡的数量20.1株/ha超过分布在凸坡的数量。

2.4.4海拔

病死树发生轻微林分中,病死树主要分布于100～150 m、300～350 m以及350～400 m,其中100～150 m病死树数量最多(P<0.05),占总数48.18%。海拔位于100～150 m的单位面积病死树数量远超其他海拔分级,为8.4株/ha。

病死树发生中等林分中,病死树主要分布于300～350 m与250～300 m,共占总数53.04%,其中300～350 m病死树数量最多(P<0.05),占28.93%。海拔位于100～150 m的单位面积病死树数量最多,为7.05株/ha(图8)。

图8 不同海拔病死树数量分析Fig.8 Analysis of the number of dead trees in differentaltitude

病死树发生严重林分中,病死树主要分布于300～350 m与250～300 m,共占总数74.25%,其中300～350 m病死树数量最多,占48.49%。位于200～250 m的单位面积病死树数量要远超过其他海拔分级,为33株/ha。

3 结论与讨论

3.1 松材线虫病病死树空间分布

嵇保中等人(2000)应用频次分析方法研究了林分中原有松树和松材线虫病病死树的空间分布型,发现林分内原有松树的分布为随机分布,而松材线虫病病死树则呈现为聚集分布。说明松材线虫萎蔫病病死木的聚集分布不是由林分内原有松树的分布所引起,更大程度是由媒介昆虫活动特性决定。松墨天牛Monochamusalternatus善于飞行,从海上岛屿到大陆隔海一次性能够飞行3.2 km,但是其飞行意愿较低,很少能观察到其在白天飞行,且有研究表明纯林中松墨天牛在100 m的范围内飞行和活动(朋金和等, 1997; 来燕学, 1998; 张心团等, 2004)。这导致了林间松墨天牛的卵、幼虫及树体内刚羽化的成虫都呈聚集分布(Shibata, 1984),并进一步导致了其寄主松树受害之后分布格局呈现聚集分布。本研究根据各项判定参数计算得出和平县阳明镇、紫金县紫城镇、东源县义合镇3地松材线虫病病死树分布均呈聚集分布,这与前人的研究结果一致,并且研究结果更加形象直观。

五点取样法、双对角线取样法、平行线法和“Z”形抽样法是包括昆虫学调查在内的实际调查中被广泛采用的方法,选择适当的取样方法,常可获得代表性较高的样本值,得到的数据更具有实际价值(金瑞华等, 1982)。例如,平行线抽样法进行绿鳞象甲Hypomecessquamosus(董丽丽和陈向阳, 2010)和枣飞象Scythropusyasumatsui(阎雄飞等, 2019)田间调查最理想,“Z”形取样法和对角线取样法进行异迟眼蕈蚊Bradysiadifformis(吕玉华等, 2022)田间调查时结果较理想。本研究在使用不同抽样方法计算松材线虫病病死树数量后发现,双对角线法、平行线法、“Z”字法等3种不同抽样方法调查平均数与总体平均数比较均有明显差异,五点取样法所得平均数与总体平均数没有明显差异。在以后调查病死树数量的工作中,可优先选择五点取样法。

目前此类以病死树空间分布格局为目标的研究较少,仅有嵇保中等(2000)以地面调查数据为基础的报道,更多的是对有害昆虫的空间分布型研究(陆俊姣等, 2021; 殷郑艳等, 2022)。而以无人机航拍图为基础开展病死树空间分布格局研究,可使研究对象更加直观,数据量更充足。但一个时间点的数据本质上还是静态数据,无法真正体现松材线虫病发生发展的流行病学趋势,因此下一步需要开展等时间间隔的长时序研究,如每月拍摄1次,持续1-2年。

3.2 不同立地因子的病死树分布规律

立地因子不仅是树木生长发育的基本条件(蔡飞等, 1997),对生物多样性的空间分布也有较大的影响。坡向、坡度、海拔等地形地貌因子在一定程度上控制着太阳辐射和地表水分的空间再分配,形成局部生态环境的小气候条件,是导致各种生态现象和过程发生变化的根本性因素,如空气湿度随着海拔的升高而下降(Hikosaka, 2002;Gong, 2008)。因此,空间尺度较小时,立地因子与病虫害的发生有密切关系(陈宏伟等, 2011)。例如,落叶松毛虫Dendrolimussuperans通常发生在海拔800 m以下的低山丘陵区,在背风朝阳,光照充足,气候温暖,土壤贫瘠的环境下发生严重,在向阳的三面环山山谷中的森林缓坡内则更为严重(李莉等, 2002);栗山天牛Massicusraddei在阳坡、坡顶发生较重(赫传杰等, 2008);昆嵛山腮扁叶蜂CephalciakunyushanicaXiao多数分布在海拔500 m以下,但极少分布在海拔较低处(<200 m),阳坡的林分内虫口密度较高(胡瑞瑞等, 2021)。

在日本,松材线虫病首先在低海拔区域发生并蔓延,随后向其他区域扩散,在400 m海拔以下的地区,松材线虫病发生最为严重;随着海拔的上升,疫情逐渐减轻;而在海拔超过700 m的地区,松林中几乎没有松材线虫病发生(Mamiya, 1983)。这是因为松材线虫发育的最适宜温度在10～28℃之间,一旦温度超出这个区间时就会不同程度抑制松材线虫的发育和繁殖(Kosaka &Ogura, 1995; 宋玉双等, 1996; 顾晓丽, 2016)。海拔高度影响局部生境中的水热状况,进而影响松材线虫的生长发育过程,最终导致松材线虫病表现出海拔差异性。同样地,方位和坡向也会影响到局部生境中光照、水热状况及土壤理化性质。在一定范围内,南坡所获得的总辐射能是北坡的1.6～2.3倍,因而南坡的温度较高。其次,由于马尾松是阳性树种,喜光喜温,南坡、阳坡的马尾松数量较多。坡度也是影响立地条件尤其是水分条件重要的因子,不同坡度对土壤含水量变化有显著影响(苏子龙等, 2013)。

本研究中,松材线虫病病死树数量在一定范围内随着海拔逐渐上升,在低于100 m和超过500 m的地方病死树极少分布;在同一个病死树聚集区域,松材线虫病病死树的分布呈现由南坡向西坡和东南坡逐渐扩散、由阳坡向半阳坡逐渐扩散的趋势。其分布规律和上述已有报道的松材线虫病发生规律相符合,且更为细化,有利于在生产实际中对防治工作的指导。综合本研究和前人研究结果,未来在防治松材线虫病时,当疫情未发生或者发生较轻时,可将100～150 m海拔和南坡作为重点进行监测,尽可能将疫情控制在小范围内。疫情开始扩散之后,需将150～500 m海拔,西坡和东南坡也纳入防治重点。

3.3 松树总量对病死树数量的影响

本研究主要分析了多个地理要素和松材线虫病疫情发生程度之间的关系,没有涉及调查地点的松树总量。广东省多地已观察到松树大量死亡后,残余松树中的病死树明显减少。这一现象说明松树总量与病死树数量可能存在相关性,即松树总量可能会影响到本试验的结果。但要通过试验证实这一相关性,必须首先获取相对精确的松树总量。目前从可见光航拍影像上没有开展这类工作的报道,主要原因可能是现有分析方法,无论是常规遥感分析还是深度学习等技术,都无法获得可达实用水平的准确率。这一问题的解决有待激光雷达和人工智能的结合与继续发展。