互花米草(Spartina alterniflora)起源于美洲大西洋沿岸和墨西哥湾, 是滨海盐沼常见的禾本科(Gramineae)多年生草本植物, 具有区域性入侵特征^[1]。互花米草对温度变化具有一定适应性, 分布纬度跨度大, 适宜生长在广阔的潮间带上^[2]。互花米草在许多引种地快速蔓延, 包括北美西海岸、欧洲、新西兰和中国沿海等^[3-4]。1990—2015年中国互花米草面积从4375.5 hm²增加到54579.7 hm^2[5], 成为中国东海岸分布最广的1个入侵种^[6-7], 改变了原来生境特征与生物群落生态结构^[8-9], 严重威胁生态系统的稳定性和功能^[10-12]。近年来, 如何防范互花米草入侵成为中国生物入侵问题研究的焦点之一。广西沿海是互花米草适宜分布区, 主要分布在丹兜海、铁山港湾、廉州湾等地^[13], 广西山口红树林保护区是互花米草扩散和蔓延面积最大的区域^[14], 其面积从2003年的167 hm²增加到2015年的413.00 hm^{2[5, 15]}已对当地自然生态系统产生严重的负面影响, 亦成为学术界关注的重点区域。

目前关于互花米草扩散的研究成果包括监测、扩散模拟、种群空间分布、繁殖方式、驱动力等方面。在监测和扩散模拟方面, 有学者利用遥感影像对互花米草进行监测^[16-18]或基于元胞自动机构建互花米草种群扩散模型^[19-20]。在种群分布方面, 有学者利用景观格局指数、质心变化等方法分析互花米草群落的景观空间分布特征^[21]。在繁殖方式方面, 前人从互花米草种群繁殖方式探讨其扩张能力^[22], 或从种群的入侵力、生态系统可入侵性和入侵通道分析互花米草的爆发机制^[2]。在影响因子及驱动力方面, 对于中国沿海互花米草扩散的影响因素有海温、海流纬向分布、海盐度、温度、降雨量^[23], 有学者以气候、水文、经济、社会和政策为驱动力因子分析湿地类型景观格局动态变化^[24-25]。研究驱动力的方法主要包括Pearson双变量相关性分析法^[26]、灰色关联度分析方法^[27]以及Logistic模型等^[28], 其中灰色关联度分析方法对于样本数量较小的情况下较为适宜^[27]。

综上, 迄今关于互花米草的研究成果颇为丰硕, 然而结合互花米草扩散动态并探讨其驱动力研究成果相对较少。本研究通过分析2009—2019年间广西山口红树林保护区互花米草空间分布、面积变化动态、质心的变化与景观格局, 探讨互花米草在红树林湿地中扩散特征与规律, 并通过灰色关联度分析法探讨互花米草扩散动态的驱动力, 旨在为互花米草监控、红树林资源保护和管理提供科学依据。

1 数据与方法 1.1 样地概况

广西山口国家级红树林生态自然保护区(中心位置地理坐标21°28'N, 109°43'E)距离广西北海市105 km, 东与广东省湛江红树林保护区接壤, 由广西合浦县东南部沙田半岛的东西两侧海岸及海域组成(图 1), 海岸线长约50 km, 面积8000 hm², 保护区于2002年1月被列入国际重要湿地名录。该保护区地属亚热带海洋性季风气候, 年均气温23.4℃, 年平均降水量为1887.54 mm。保护区滩涂主要为淤泥质, 分布着红树林建群种有红海榄(Rhizophora stylosa)、木榄(Bruguiera gymnorrhiza)、桐花树(Aegiceras corniculatum)、白骨壤(Aricennia marina), 互花米草主要分布在丹兜海潮滩、英罗潮滩, 呈单优群落生长或与红树林混生。

1.2 数据来源及处理 1.2.1 遥感影像数据

以2013年12月高分一号全色和多光谱高分辨率卫星影像为主要数据源, 分辨率分别为2 m与8 m。利用ERDAS9.2遥感影像处理平台做正射校正、数据配准与融合处理, 并组合2、4、3波段影像, 采用人机互译的方式进行解译。利用覆盖研究区域1 10000地形图进行遥感影像校正。以同期QuickBird影像为基本验证数据, 随机抽取100个互花米草斑块, 利用对比法进行精度验证^[16]。2013年研究区互花米草抽样斑块总体面积误差为1.1%, 斑块面积总体重叠率达到95.8%。

1.2.2 互花米草面积数据

研究区1979、2000、2003、2007、2008、2011、2015年互花米草面积数据来源于参考文献(表 1)。2009年10月、2019年12月互花米草数据来源于Google Earth影像, 影像分辨率为0.61—2.4 m ^[31-32]。利用ArcGIS的“线切割面”的功能目视解译研究区Google Earth影像。随机抽取100个斑块进行现场手持GPS绕测精度验证, 2019年研究区抽样斑块总体面积误差为1.3%, 斑块面积总体重叠率达到95.4%。

本次研究对象为互花米草, 主要分布在滩涂上, 因此将保护区滩涂用地进行分类: 互花米草、红树林、光滩、养殖滩涂。此外, 本研究还解译了海堤内废弃虾塘中的互花米草。辅助数据有ALOS卫星12.5 m Dem数据、潮位信息数据。本研究采用ArcGIS 10.4软件进行矢量数据处理及制图。

1.2.3 气象、经济、社会数据

年均最高温、最低温、年降雨量数据来源于北海历史天气数据(http://lishi.tianqi.com/beihai/201911.html)和国家气象科学数据中心地面资料数据(https://data.cma.cn/site/index.html)。地区国民生产总值、人口数据来源于北海市统计年鉴(http://xxgk.beihai.gov.cn/bhsdfzbzwyhbgs/zjtd/bhnj_98776/)；赶海人数数据来源于广西红树林研究中心开展的山口红树林保护区本底调查。

1.3 研究方法 1.3.1 土地类型转移矩阵

用土地类型转移矩阵来分析用地类型的结构特征及不同用地类型变化方向, 根据不同时期的用地类型, 按公式(1), 在GIS中进行地图代数计算^[33]。

(1)

式中, C_ij表示K时期至K+1时期用地类型变化, A_ij^k、A_ij^k+1表示两个时期的用地类型面积(hm²)。利用ArcGIS 10.4软件叠置的分析功能, 得到用地类型转移矩阵。

1.3.2 质心的变化

在ArcGIS 10.4中标出各互花米草斑块的质心点, 记录斑块质心的经纬度坐标。通过计算互花米草的面积加权质心变化, 揭示互花米草在空间上的扩散方向和速度, 计算方法按公式(2)^[21]:

(2)

式中, X_c和Y_c是按面积加权的互花米草质心坐标, X_i和Y_i是互花米草的第i个斑块的质心坐标, C_i为互花米草的第i个斑块的面积, n为互花米草斑块的总数目。

1.3.3 景观格局分析

对2009、2013、2019年山口红树林保护区互花米草景观格局指数开展分析, 以了解研究区互花米草格局及其变化特征。通过Fragstats4.1进行计算斑块数量(Number of patches, NP)、斑块密度指数(Patch Density, PD)、平均斑块大小(Patch Density, MPS)、分维数(Dimension, D)和最大斑块指数(Patch Density, LPI)等5个景观指数^[34]。

1.3.4 灰色关联度分析方法

灰色关联度分析方法, 是指两个系统或两个因素之间关联性的度量, 目的在于找到系统间影响因素的关系, 以反映各评价对象与参考序列的关联关系。本研究通过灰色关联度分析互花米草景观格局指数与人为因素、气候因素的关联性, 在具体研究中对指标数据利用初值法做无量纲化处理, 计算比较序列与参考序列对应元素的绝对差值, 分辨系数取0.5, 最终定量各因素之间的数值关系。关联度值大于0.5时, 说明关联度显著；关联度值越接近1, 关联度越强^[27]。

2 结果与分析 2.1 互花米草空间分布和扩散速度 2.1.1 互花米草空间分布

2009—2013年间, 广西山口红树林保护区滩涂主要分布着红树林和互花米草斑块, 互花米草在保护区丹兜海潮滩、英罗港海塘村、英罗村滩涂上蔓延, 与红树林形成竞争滩涂空间的群落(图 2)。2019年养殖滩涂斑块面积增加13.77 hm², 主要是由于当地居民刈割处理禾荣村潮滩(中心位置: 经度109°38′19", 纬度21°35′15")互花米草, 在滩涂上养殖红树蚬(Geloina erosa)。2009—2013年间互花米草面积最大的斑块分布位置不变(中心位置: 经度109°39′45″, 纬度21°33′44″), 但面积变小, 2009年、2013年与2019年依次为64.34 hm²、53.91 hm²、53.54 hm²。

丹兜村海堤坍塌(中心位置: 经度109°40′13", 纬度21°35′19"), 互花米草种子随着海水倒灌进堤内, 互花米草在废弃的虾塘蔓延。2009年堤内互花米草斑块24个, 面积3.57 hm², 2013堤内斑块53个, 面积37.38 hm²。2013—2019年间, 当地居民对废弃虾塘进行整治将堤内互花米草用地转为虾塘, 至2019年堤内互花米草斑块7个, 面积约7.06 hm²。

2.1.2 互花米草扩散速度及面积变化趋势

2009—2019年间, 保护区互花米草面积呈现增加趋势, 但增长速度下降。2009年互花米草面积为301.04 hm², 至2019年增加到454.24 hm²。互花米草年均增长率大幅下降, 2009—2013年年均增长率为7.60%, 2013—2019年则下降为1.99%(表 2)。

2009—2013年间红树林面积呈增加态势, 其年均增长率为3.42%, 但面积增长率小于互花米草, 而2013—2019年呈减少趋势, 主要是位于禾荣村马扒欄滩涂的红树林死亡(中心位置: 经度109°38′36″、纬度21°35′40″)。

以2009—2019年间广西山口自然保护区的面积数据为因变量(表 2), 通过Excel绘制互花米草面积散点图(图 3), 得到互花米草面积预测模型式(3):

(3)

式中, x为预估测的年份, y为预测的互花米草面积, R²=0.785, 利用SPSS 22.0软件进行方差分析、模型统计检验, F=23.013, P=0.001 < 0.01, 该模型具有统计学意义。结果表明互花米草面积呈线性递增态势。如果互花米草不进行控制, 预测模型显示2029年互花米草面积将达到555.387 hm², 2039年达682.417 hm², 界时将有更多的适林滩涂被侵占。这也印证了前人的研究结论^[14]。

2.1.3 用地类型转移分析

通过保护区10年间土地利用转移矩阵分析, 可进一步说明研究区互花米草面积变化情况。总体而言, 2009—2019年间研究区各地类中活跃程度位列前三的依次为养殖滩涂、互花米草、红树林(表 3)。互花米草侵占滩涂的空间速度比红树林快, 红树林将面临生长空间的限制。2009—2019年间, 光滩转化为互花米草是其转化为红树林面积的1.507倍。2019年红树林面积为868.98 hm², 增加面积主要来自互花米草(17.99 hm²)、光滩(137.52 hm²)转入。位于那潭村滩涂(经度109°39′56"、纬度21°36′23")互花米草转入为红树林, 主要是相关部门在此种植人工林。2019年互花米草面积为454.24 hm², 主要来自红树林(13.77 hm²)、光滩(207.26 hm²)、堤内虾塘(3.49 hm²)转入；2019年保护区滩涂养殖面积来自红树林(1.21 hm²)、互花米草(3.37 hm²)以及光滩(9.19 hm²)转入。

2.2 保护区互花米草扩散特征

2009、2013、2019年互花米草的质心坐标都位于丹兜海潮滩(图 4)。互花米草斑块质心整体移动方向: 2009—2013年向西北方向偏移, 2013—2019年向东南方向偏移。2009—2013年向北移动速度较快, 而2013—2019年向南移动速度较快。前一阶段质心变化, 主要是因为互花米草扩散进入保护区西北侧稀疏的红树林, 而后一阶段互花米草在英罗湾、烟楼下等潮滩平行海岸线向海一侧扩散。

保护区互花米草种群扩散表现为点源扩散、多点爆发的特点。以丹兜潮滩上互花米草斑块(中心位置: 经度109°39′8", 纬度21°34′53")扩散过程为例(图 5), 图中2009年互花米草斑块由破碎的多个点源构成, 随着时间变迁多点爆发, 至2013年多个点源连接形成面状的斑块, 至2019年形成片状的互花米草斑块。这种扩散的特点, 被邓自发等称之为短程“流”式扩张^[2]。互花米草对土壤要求不高, 在淤泥质的滩涂上生长最好^[35], 而山口红树林保护区为海湾淤泥滩涂生境, 有利于互花米草生长扩散, 因此在丹兜海潮滩上形成大面积互花米草优势群落。

2.3 互花米草景观格局变化

本研究用高度浓缩景观格局信息的景观格局指数反映斑块结构组成和景观变化特征。如表 4所示, 2009—2019年红树林平均斑块大小先上升后降低, 分维数数值和斑块密度值均先降低后上升, 表明期间斑块形状趋向复杂变化, 红树林斑块有破碎化趋势。

互花米草斑块数量2013年较2009年减少, 主要是小斑块连接成大斑块的缘故, 2013 —2019年斑块数量增加, 互花米草斑块整体质心偏移, 位于丹兜海潮滩的烟楼下、永安、上新村、以及英罗、海塘的滩涂上有新斑块生成。

互花米草分维数2019年 < 2009年 < 2013年, 表明2009—2013年互花米草爆发高峰期, 斑块复杂, 并且随着时间推移, 至2019年斑块逐步成型。相对2013年, 2019年的互花米草斑块形状由复杂变为简单。2009—2019年间互花米草斑块密度数值先减小后增大、平均斑块大小数值先增大后减小, 说明期间互花米草呈现破碎化趋势。

2.4 互花米扩散的驱动力因素

本研究利用社会环境因素结合灰色关联度分析法研究样地互花米草动态变化驱动力。引起互花米草扩散动态变化的社会环境因素分为人为因素和自然因素, 其中人为因素主要体现在地区国民生产总值变化、山口赶海人口比例；自然因素包括年均最低温、年均最高温(图 6)。因研究期间年均气温与年均降雨量变化趋势不明显, 未纳入到影响因素中。

表 5所示, 分别以2009—2019年间研究区互花米草分维数、最大斑块指数、破碎化指数、斑块密度、面积作为参考序列, 针对4个驱动因素赶海人口比例、年均最高温、年均最低温、地区生产总值, 进行关联度分析。各个影响因子关联度值均大于0.5, 各因子选择具有一定合理性^[27]。研究结果表明, 2009—2019年间互花米草动态变化受到自然因素与人为因素共同影响, 互花米草景观格局指数对年均最高温、赶海人口比例、年均最低温较敏感。其中, 年均最高温与互花米草分维数、破碎化指数、面积相关性最大；赶海人口比例对最大斑块指数、斑块密度、斑块数影响最大；年均最低温度与破碎化指数、互花米草面积、最大斑块指数、斑块密度、斑块数的关联度值均排在第二位；地区生产总值与互花米草景观格局指数的关联度相对较小。

2009—2019年间影响互花米草面积变化的各个因子关联度排序依次为: 年均最高温>年均最低温>赶海人口比例>地区生产总值。因此研究区域互花米草面积主要是受到气候因素的影响, 其次为人为因素。从气候因素上看, 随着全球温度的上升, 研究区域年均高温和年均最低温都有上升趋势(图 6), 互花米草面积与研究样地的年均最高温正相关(R²=0.9958), 与年均最低温正相关(R²=0.9304)(图 7)。从人为因素来看, 互花米草面积与地区生产总值正相关(R²=0.8448), 互花米草面积与研究样地赶海人口比例呈负相关(R²=0.9986)。2009—2019年间, 北海市地区生产总值为从335亿元上升到1300.80亿元, 年均增长率为14.53%(R²=0.9805)。2009—2019年间从事赶海活动的人口比例下降, 说明该地经济发展为当地居民提供了更多的就业岗位与多样化的经济收入来源, 赶海人口比例相应减少, 有利于互花米草扩散。

3 讨论

(1) 气候变化对互花米草的景观格局指数和面积变化有一定影响。互花米草分维数、破碎化指数、面积对年均最高温变化表现较为敏感, 互花米草面积与当地的年均最高温、年均最低温正相关。Chuine等认为增温显著提高了C4植物的竞争力^[36], 本研究中互花米草面积年均增长率大于红树林, 可能是因为在温度变化背景下, 互花米草的竞争力高于红树林的缘故。

(2) 在全球变暖的趋势下, 2009—2019年间山口保护区互花米草总面积呈现增长趋势, 但扩散速率呈现下降趋势。这与广西海岸互花米面积草整体变化趋势基本一致^{[14, 16]}。究其原因, 一方面是由于温度的影响, 祝振昌等研究结果认为互花米草萌发率随温度升高呈现先增大后减小的趋势^[37]。另一方面因为经过前期的互花米草爆发式扩散, 面积迅速增加, 但受保护区滩涂空间限制, 其扩散速率降低。在研究区域互花米草入侵的红树林空间主要是保护区西北部的红树林密度小的区域, 在红树林密度大、树形高大的英罗港区域, 互花米草主要是位于红树林斑块外缘^[23], 未能入侵到斑块内部空间中, 这与前人研究互花米草容易向低矮、稀疏的红树林内扩散的结论一致^{[14, 38]}。

(3) 从2009—2019年土地利用变化可以看出, 2019年局部滩涂上的互米草(面积3.37 hm²)被刈割, 互花米草滩涂转为养殖滩涂。当然, 养殖滩涂是否长期有效抑制互花米草生长, 需要后期进一步监测。此外通过滩涂养殖取代互花米草的方式, 其生态影响也需要进一步关注, 因养殖滩涂面积还包括红树林面积(1.21 hm²)、光滩面积(9.19 hm²)转入, 人为干扰加剧, 将对红树林生长造成一定影响。

(4) 本研究中驱动力分析的因素主要包括自然因素和人为因素, 在气候因子中选择年均值, 具体月份的气候因子影响差异还需要进一步研究。在后期研究中, 拟将互花米草入侵前后土壤的理化性质变化、地形、土地利用纳入驱动因素中进行深度分析, 为红树林湿地环境生态保护提供依据。

4 结论

(1) 广西山口红树林保护区互花米草呈现扩散趋势, 与红树林争夺生长空间。2009—2019年间互花米草面积呈现增长趋势, 在光滩上生长优势较红树林明显, 且年均增长率大于红树林, 生长占用滩涂空间, 将使宜林空间减少。期间, 红树林趋向破碎化, 其中主要是由于互花米草侵入红树林边缘, 局部入侵斑块内部的缘故, 特别是丹兜海西北侧稀疏的红树林, 使红树林斑块形状趋向复杂变化。近年互花米草质心向保护区东南方向偏移, 应注意防范。

(2) 2009—2019年间互花米草动态变化受到人为因素与自然因素共同影响。在人为因素方面, 赶海人口比例是景观格局变化的主要影响因子之一, 与最大斑块指数、斑块密度、斑块数关联度最大。在自然因素方面, 气候变化是互花米草景观格局及面积变化的重要影响因素。2009—2019年间, 年均最高温对互花米草分维数、破碎化指数变化影响相对较大。互花米草面积变化的主要影响因子亦为气候因素, 与年均最高温、年均最低温均呈正相关。随着全球气候变暖, 互花米草面积增加, 红树林生长将面临严峻挑战。

致谢: 感谢福建师范大学地理科学学院祁新华教授对论文提出宝贵的修改意见。