珠海淇澳岛人工次生无瓣海桑纯林的植被碳储量变化

孙学超， 黄展鹏， 张琼锐， 文 嫱， 黄志琴， 王晨茜， 张若琪， 徐颂军*

(1. 华南师范大学地理科学学院， 广州 510631； 2. 洲际环境科学研究院(广州)有限公司， 广州 510660)

随着全球温室气体排放的增加，气候变暖越来越被国际社会所重视。湿地生态系统因其巨大的碳储量和较强的温室气体排放能力成为了全球变化科学中的前沿和热点。

湿地占陆地总面积的4%～5%，但其碳储量高达450 Gt，相当于陆地生态系统总碳储量的20%[1]。其中，红树林地上部分碳储量约为2.5 Gt，土壤的碳储量为4～20 Gt[2]。红树林覆盖了地球表面0.1%的面积，却固定了相当于大气中5%的碳，其单位面积的有机碳储量是陆地森林的5倍之多[3]。然而，20世纪70年代开始，人类活动使红树林遭受到严重的破坏，全球超过1/3面积的红树林消失。作为目前中国最大的红树林恢复区的珠海市淇澳岛红树林，也曾遭受过严重破坏。淇澳岛历史上生长有600 hm2的红树林，1950年开始的大规模开垦致使红树林遭受严重破坏，到1998年时，原生林仅剩大围湾最内侧 32.2 hm2的秋茄-桐花树群落[4]。1999年，珠海市政府开展了红树林修复计划，截至2017年已成功营造了约600 hm2的以无瓣海桑为主的人工红树林[5]。目前有关红树林碳库的研究主要集中在红树林植被各器官间的碳含量关系[6]、植被的固碳能力[4,7-10]和湿地土壤的固碳能力方面[11]。其研究目的多为计算单一年份的红树林生态系统碳库，而植被碳储量是该碳库中的重要部分，约占红树林生态系统总碳库的30%以下[12]。相较土壤碳库，红树林植被碳储量的年际变化更大[13]，这方面的研究是目前较少被关注到的。

本研究结合实地采样和多年高分辨率遥感影像数据，对1999—2019年间淇澳岛无瓣海桑纯林的生长范围、年际生物量和生物碳储量进行定量研究，同时分析了引起这些变化的原因及无瓣海桑的固碳能力，以期评价以无瓣海桑为主要树种的中国亚热带红树林生态修复方法的效果，为亚热带红树林生态系统的固碳研究提供基础数据，为今后生态修复的树种优选提供理论依据。

1 研究区、数据来源和方法

1.1 研究地概况

淇澳岛红树林自然保护区位于广东省珠海市淇澳岛西北部(113°36′40"E～113°39′15"E，22°23′40"N～22°27′38"N)、珠江内西侧，东与香港、深圳隔海相望，南联万山群岛及南海，西靠金星门，北与中山、番禺、东莞毗邻，总面积为23.9 km2。该自然保护区内地势落差较大，北高中低，具有陆地、岛屿、海岸、低山和丘陵等丰富多样的地理地貌资源。地质条件以花岗岩为主，表层是黄沙黏土，土壤类型是滨海盐渍草甸沼泽土，岛屿东侧湾多为沙岸，西北两侧湾多为淤积泥滩。该地属南亚热带海洋性季风气候区，年平均气温为23.5 ℃，降雨量为1 948.9 mm，日照时数为1 910.0 h。

1.2 样方设置及调查数据

样方位置的设置主要遵循以下原则[14]：

(1)代表性：根据研究对象的群落特征，尽可能少的点获得最具代表性的群落生长状态信息；

(2)兼顾性：考虑红树林湿地采样难度，样方位置尽可能归一或靠近布设；

(3)可达性：样地调查中的各项信息的获取需要有可靠的技术方法和仪器设备。

由于乔木的胸径对计算其生物量有重要影响，因此，采用《植物群落清查的主要内容、方法和技术规范》[15]为指南进行胸径测量。乔木样方采用标准样方(10 m×10 m)，用GPS定位仪记录样方四角坐标点，测定样方内乔木的树种、胸径、树高、数量和冠幅并拍照记录。

采样时间为2019年，共设置了20个样方(图1)。由于研究区域的林分为无瓣海桑纯林，样方内的乔木均为无瓣海桑。本研究沿淇澳岛红树林保护区内航河道选取了15个样方(A1～A15)，该区域的无瓣海桑林较茂密；在离河道较远处无瓣海桑较稀疏的区域又选取了5个样方(B1～B5)，共计20个样方(图1)。用随机拆分法将样方数据按照3∶1的比例拆分，其中 15个样方数据(编号1～15)用于构建生物量反演模型， 5个样方(编号16～20)数据用于模型的精度评价(表1)。

图1 无瓣海桑纯林样方的位置

表1 研究区设置的无瓣海桑纯林样方的基本信息

1.3 研究方法

1.3.1 植被碳储量及生物量计算 植被碳储量与生物量的关系较为密切，其估算公式如下[14]：

C植=B×0.5，

(1)

其中，B为生物量。

无瓣海桑中树干、树枝、树叶、树根、树皮、花果的生物量分别用以下异速生长方程确定[16]：

lgB树干=0.360 7+0.330 2×lg(DBH2×H)，

(2)

lgB树枝=0.159 0+0.387 9×lg(DBH2×H)，

(3)

lgB树叶=-0.756 0+0.435 5×lg(DBH2×H)，

(4)

lgB树根=-0.039 2+0.421 1×lg(DBH2×H)，

(5)

lgB树皮=-0.379 0+0.355 9×lg(DBH2×H)，

(6)

lgB花果=-2.345 6+0.379 1×lg(DBH2×H)，

(7)

其中，B为生物量(kg)，DBH为树胸径(cm)，H为树高(m)。

1.3.2 遥感影像预处理 本研究采用1999年的Landsat5卫星TM遥感影像、2003年和2009年的QuickBird影像、2013年的GeoEye-1影像以及2019年的高分二号遥感影像。用ENVI5.3对遥感影像进行辐射定标和大气校正后，基于2019年的遥感影像对其他年份的遥感影像进行正射校正，再计算用于生物模型的构建和反演的比值植被指数和纹理数据方差(Variance)，具体公式如下[15]：

(8)

Mean=∑i∑jp(i,j)×i，

(9)

Variace=∑i∑jp(i,j)×(i-Mean),

(10)

其中，SR为简单植被指数，ρNIR、ρRED分别为近红外波段、红光波段在遥感图像里的辐射亮度值；i、j分别为像元在影像中的行数、列数，p(i,j)为处理的模板矩阵中(i,j)位置的值；Mean和Variance分别是灰度共生矩阵统计量的均值和方差，Mean表现了影像纹理的均匀程度，Variance表现了影像像元值与均值的偏差，二者均可通过计算处理窗口中像元的方差值得到。

1.3.3 无瓣海桑范围提取 淇澳岛红树林湿地生态恢复的主要手段是引种无瓣海桑来抢夺互花米草入侵的区域，给本土红树林留下生存空间，再利用无瓣海桑和本土红树林树种混合进行物种恢复，因此，当前淇澳岛红树林的主要林分是无瓣海桑林[17]。本研究在2019年田野调查的基础上，采用人工神经网络的方法，对1999、2003、2009、2013、2019年的遥感影像进行解译，研究对象仅为无瓣海桑纯林，不包括少量的无瓣海桑与其他红树林的混交林。通过利用目视解译得到的验证集来计算混淆矩阵，以判定2003年和2009年的无瓣海桑纯林分类精度；通过2013、2019年实地考察所得的验证集来计算混淆矩阵，以判定2013、2019年无瓣海桑纯林的分类精度。

1.3.4 生物量反演 在利用遥感影像估算红树林生物量并折算成红树林植被碳储量方面已有诸多研究[18-20]，目前常用的有多光谱影像、高光谱影像、高分辨率影像、雷达影像以及以上影像的多源组合，其中，加入了雷达影像的多源组合影像有较高的准确率。考虑到研究区域的时空范围和实际需要，本研究选取包含红、绿、蓝和近红外波段的高分辨率遥感影像作为反演数据。在反演生物量的自变量选取方面，选取了比值植被指数(SR)及纹理特征中的方差特征(Variance)。

本研究在进行生物量反演时，遵循如下步骤：首先，根据2019年15个无瓣海桑纯林样方的数据，分别计算单株无瓣海桑树干、树枝、树叶、树根、树皮和花果的生物量后，加总得到单株无瓣海桑的生物量，将单个样方的所有单株生物量加总，除以样方面积，得到该样方的单位面积生物量；然后，将预处理好的2019年高分二号遥感图像的SR和Variance与所得样方的生物量进行多元回归建模，得到生物量反演模型；接着，将5个样方的SR和Variance代入模型，得到预测值，计算其相对误差；最后，如果模型预测值和真实值偏离不大，表明模型的准确度较高，则可根据得到的生物量反演模型对2003、2009、2013、2019年的遥感影像进行生物量反演。

2 结果与分析

2.1 无瓣海桑生物量模型构建的结果

将2019年11月17日在淇澳岛红树林保护区得到的野外采样数据代入异速生长方程(2)～(7)，从而分别得到20个样方的无瓣海桑的总生物量。由结果(表2)可知：20个样方的无瓣海桑的总生物量最高为4.08 t，最低为0.64 t。

表2 20个样方中无瓣海桑的生物量

2.1.1 单位面积生物量和SR、Variance的相关性 将2019年各个样方的无瓣海桑的单位面积生物量与样方对应的SR、Variance进行相关性分析，得到了三者之间的相关关系：无瓣海桑的单位面积生物量与Variance、SR的相关关系较为密切，分别呈显著相关(P<0.05)、极显著相关(P<0.01)，相关系数分别为-0.720、-0.605。可见无瓣海桑的单位面积生物量与这2个变量的相关性高，且为负相关，即随着SR和Variance值的增大，无瓣海桑的单位面积生物量减小；Variance比SR解释了更多的无瓣海桑的单位面积生物量变化。

(11)

(12)

(13)

由拟合模型的结果可知，模型(11)～(13)的P值均小于0.05，说明模型具有统计学意义。其中，模型(11)的相关系数为0.605，R2=0.367，RMSE=16.01 kg/m2；模型(12)的相关系数为0.720，R2=0.519，RMSE=19.03 kg/m2；模型(13)的相关系数为0.837，R2=0.700，RMSE=15.63 kg/m2。可见将SR和Variance同时作为自变量时，所建立的模型的精度更高。利用模型(13)对15个样方的无瓣海桑的单位面积生物量进行预测，预测值和实测值间的交叉检验结果(图2)显示：交叉检验决定系数为0.704 6，可见预测精度较高。将用于验证的5个样方代入模型(13)，预测值和实测值的相对误差均值为9.18%，预测值与实际值偏离不大。可见该模型不仅在样本集的拟合效果好，在测试集的拟合效果也较理想，能够满足研究需求。

图2 15个样方的单位面积生物量的实测值与预测值对比

2.2 无瓣海桑纯林面积、生物量和碳储量的变化

通过人工神经网络解译2003、2009、2013、2019年的遥感图像，得到无瓣海桑纯林的分布范围(图3，图4)。2003、2009、2013、2019年的无瓣海桑纯林的识别精度分别达到95.87%、91.68%、90.33%、90.04%，精度较高，满足分析需要。

图3 2019年淇澳岛无瓣海桑纯林的生长范围

1999年是淇澳岛红树林生态修复的初级阶段，且根据1999年的Landsat影像(图5)可以判断，后期引种无瓣海桑的地区基本上是裸地，此时的红树林仅大围湾附近有32 hm2的秋茄桐花树群落，因此无瓣海桑纯林面积计为0 hm2；由图4和表3可知：2003年，无瓣海桑纯林开始出现在研究区沿大坝向海逐渐延伸的位置，此时的密度较为稀疏，面积为29.52 hm2；2009年，无瓣海桑纯林的生长面积逐渐向北扩大(面积为114.25 hm2)，林分密度也有所增加，并延伸至西部的沿海的狭长地带；2013年，无瓣海桑纯林的林分密度进一步增加，并继续在西部沿海地带延伸生长，面积为140.28 hm2；2019年，淇澳岛西部沿海已被无瓣海桑纯林包围，中部核心地带的林分密度明显升高，并向东部沿海海岸延伸，面积为332.40 hm2。可见，2003—2019年，本研究区域的无瓣海桑纯林面积一直在增长，其增长速度在不同阶段有所差异：1999—2003年，无瓣海桑的面积年均增加5.90 hm2；2003—2009年，无瓣海桑的面积年均增加12.10 hm2；2009—2013年，无瓣海桑的面积年均增加5.21 hm2；2013—2019年，无瓣海桑的面积年均增加27.45 hm2。增长的无瓣海桑区域主要是从滩涂裸地地和海洋部分转变而来，且2013—2019年期间增长最快。

图4 2003—2019年淇澳岛无瓣海桑纯林的生长范围变化

表3 淇澳岛无瓣海桑纯林面积和生物量

图5 1999年淇澳Landsat卫星影像图

由淇澳岛不同时期无瓣海桑纯林的生物量变化情况(图6)可知：2003年，淇澳岛无瓣海桑纯林的单位面积生物量为21.71～48.48 kg/m2，此时无瓣海桑纯林较为稀疏，集中分布在研究区域的中部；2009年，淇澳岛无瓣海桑纯林的单位面积生物量为6.74～63.95 kg/m2，高值点集中分布在研究区域北部，低值点集中分布在研究区域东南部；2013年，淇澳岛无瓣海桑纯林的单位面积生物量为9.41～60.4 kg/m2，高值点集中分布在西侧沿海一带，低值点分布在东侧和南侧；2019年，淇澳岛无瓣海桑纯林的单位面积生物量为18.76～58.93 kg/m2，高值点集中分布在西侧沿海一带，低值点集中分布在东南部地区。

图6 2003—2019年淇澳岛无瓣海桑纯林的单位面积生物量变化

由研究期内无瓣海桑纯林的生物量、植被碳储量及其年际变化情况(表3、表4)可知：(1)1999—2013年，无瓣海桑纯林总生物量的年均增长率和年均固碳速率不断提升；(2)总体上看，1999—2019年，研究区无瓣海桑纯林的总生物量不断增加，各年份的无瓣海桑纯林的单位面积生物量也在同步增加。

表4 淇澳岛无瓣海桑纯林总生物量、面积和植被碳储量年际变化

3 讨论

3.1 无瓣海桑解译和反演模型精度对比

本研究区域的树种为无瓣海桑纯林，在遥感影像上的纹理和波段都具有典型特点，又配合野外调查核证，因此，基于人工神经网络模型解译各个年份遥感影像的精度较高(90.04%～95.87%)，Kappa系数也大于0.8。目前提取红树林面积和红树林种间分类的精度多集中于86.7%～95.0%，在使用基于监督分类的方法对红树林区域进行识别、分类时，若遥感图像的分辨率越高且纹理越清晰，则遥感影像的分类精度越高[18,21-24]。此外，在进行淇澳岛红树林的遥感影像解译时，唐焕丽等[19]基于WorldView-2影像的8波段数据，对淇澳岛红树林进行了种间分类，分别计算无瓣海桑、秋茄和桐花树等群落的面积，得到2010年11月时淇澳岛红树林中无瓣海桑的面积为146.41 hm2，大于本研究得到的2009、2013年的无瓣海桑的面积(分别为114.24、140.28 hm2)。该情况应是其提取的无瓣海桑林中包括了混交林中的无瓣海桑。

在构建无瓣海桑反演生物量模型，本文分别用了SR、Variance以及二者结合的方式作为自变量与生物量进行回归建模。由拟合模型的结果可知：仅用Variance为自变量的模型的R2、P、RMSE值均优于仅用SR指数为自变量的模型；以Variance和SR同时为自变量的模型的R2、P、RMSE值均优于单一自变量的模型，精度较高。因此，采用二者结合的方法作为研究区红树林生物量反演模型是一个较为可靠的反演方法。截至2020年，国内利用多光谱遥感影像对红树林生物量进行反演的研究较少，已有研究运用多元逐步回归方法反演生物量，其判定系数主要在0.30～0.78之间[25]，精度波动范围较大；而反演无瓣海桑生物量的研究更少，仅有王月如[14]基于WorldView-2、Landsat8和RadarSat-2遥感影像数据进行深圳湾红树林的生物量反演，报道了无瓣海桑的生物量反演效果。该研究所采用的指数和方法与本研究类似，所得生物量反演模型的相关系数为0.858，判定系数为0.737，与本研究的反演精度基本一致。

3.2 淇澳岛无瓣海桑纯林生物量和碳储量的时空差异

当红树林湿地遭受破坏之后，对其进行人工复种的碳汇潜力十分巨大[26]。淇澳岛于1999年开始引种无瓣海桑，进行红树林湿地的恢复[4]。2003年，淇澳岛无瓣海桑纯林的总生物量已经达到了9 745.47 t，碳储量达4 872.74 t，5年期间年均可固碳974.55 t。2003—2009、2009—2013、2013—2019年的碳储量年均增加量分别为1 454.06、2 168.28、4 296.67 t。尽管在上述这3个时间段里，无瓣海桑纯林在各个年份里的面积年均增长率下降，但在2009—2013、2013—2019年期间，由于无瓣海桑纯林的碳储量基数较大，碳储量年均增加量依然呈增加趋势。2003—2019年，无瓣海桑纯林的单位面积生物量持续上升。由此可见，1999—2019年淇澳岛的红树林生物量、植被碳储量和单位面积固碳能力均呈现上升态势，淇澳岛无瓣海桑的固碳能力尚未饱和。

淇澳岛无瓣海桑纯林生物量和植被碳储量在空间分布上是有差异的(图4)。2003年生物量处于较高水平的中、东部地区在2009年的生物量呈现下降趋势。导致该变化的主要原因可能是：2008年的南方冰冻灾害[4]导致无瓣海桑林大面积死亡，而靠近西部和北部沿海地区的红树林水热条件相对较好，寒害后无瓣海桑的生物量增加迅速，从而维持了寒害之前的水平。从2009、2013、2019年的生物量分布上，可以看到研究区域西侧的红树林的生物量和植被碳储量稳步上升，而中、东部地区的红树林生物量增长较慢。除了水热条件的差别以外，人工干扰的程度也是重要的影响原因之一:研究区西侧的红树林是缓冲保育区，受到的人工干扰较少；科普展示区和山林漫步区主要分布在红树林保护区的中部、东部，受到的人工干扰相对较多；研究区域的东南部有水坝，工程用土的填充对附近的红树林也产生了一定程度的影响。

3.3 无瓣海桑的快速固碳及碳汇能力

与同处于亚热带的其他生态系统相比，淇澳岛无瓣海桑纯林在种植初期(1999—2003年)的年均固碳速率已经超过了我国亚热带区域的阔叶林[27]，是阔叶林的1.5～2倍，2003—2009、2009—2013、2013—2019年的年均固碳速率递增，且固碳速率远高于我国其余亚热带地区树种的[27]。

但是，1999—2009年淇澳岛无瓣海桑纯林的年均固碳速率均低于全球红树林的年均固碳速率(11.10～13.63 t/(hm2·a)[28-29]，主要的原因有2点：一是淇澳岛的无瓣海桑纯林在2008年南方冻灾中大面积死亡，导致无瓣海桑的生物量锐减，影响了年均固碳速率的计算；二是2009年遥感影像研究区域的无瓣海桑纯林大部分是寒害后人工复种的幼龄林，而用于全球红树林年均固碳速率计算的红树林包括各个林龄阶段的红树林植被。2009—2013年，研究区域的无瓣海桑纯林的年均固碳速率迅速提升至19.48 t/(hm2·a)，2013—2019年的年均固碳速率的增长趋势显著放缓，这2个阶段的无瓣海桑纯林的年均固碳速率均超过了全球红树林的年均固碳速率。

在这20年时间，研究区无瓣海桑纯林的单位面积植被碳储量从0 t/hm2增至168.38 t/hm2，增长幅度、生长速度和固碳能力是其他森林生态系统难以达到的[30]。这主要是因为无瓣海桑作为大乔木，其胸径和树高都超过其他亚热带红树林，拥有更强的光合氮利用效率和水分利用率[31]。同时，无瓣海桑的地下部分也有很强的固碳能力，林下土壤的碳库也高于其他本土红树林树种[32-33]。这些优势使得无瓣海桑能够快速生长并抢占有利生态位[34-35]。综上可知，无瓣海桑在快速固碳方面具有较大的优势。

如表5所示，2019年无瓣海桑纯林的单位面积植被碳储量达到168.38 t/hm2，大于广东省广州市阔叶和针叶林地(30.93 t/hm2)[36]和广西省森林(55.37 t/hm2)[37]等亚热带地区森林的碳密度，说明无瓣海桑的储碳能力高于其他亚热带森林。与其他红树林的植被碳储量[1,18,25,38-40]相比，无瓣海桑纯林的植被碳储量也处于较高水平。显然，这与无瓣海桑的快速固碳能力息息相关。

表5 不同热带、亚热带生态系统的单位面积生物量、单位面积植被碳储量

由此可见：在亚热带沿海滩涂地区的裸地上恢复植被时，可将无瓣海桑作为先锋树种进行裸地改良，以获取更高的固碳效益。

4 结论

本文基于多期遥感影像和实测样方数据，对1999—2019年期间，淇澳岛无瓣海桑的生长范围、年际生物量、植被碳储量及固碳速率进行了定量计算，并分析了引起这些变化的原因。主要结论如下：

(1)无瓣海桑在1999—2003、2003—2009、2009—2013、2013—2019年间的固碳速率分别为3.66、7.55、19.48、28.06 t/(hm2·a)。其固碳速率远高于普通的亚热带常绿阔叶林，是亚热带沿海裸地恢复植被的优选树种。

(2)淇澳岛无瓣海桑纯林的植被碳储量快速增长，2003、2009、2013、2019年的植被碳储量分别为4 872.74、15 051.14、25 892.55、55 969.24 t，单位面积植被碳储量分别为14.66、45.28、77.90、168.38 t/hm2。20年间，淇澳岛种植无瓣海桑将沿海裸地生物量提升到168.38 t/hm2。每年新增的植被碳储量还在加速上升，显然，无瓣海桑纯林的固碳速率快、碳汇能力大，且淇澳岛无瓣海桑纯林的碳汇能力还在持续上升。

(3)无瓣海桑纯林的生长明显受到人类活动的影响：人工干扰程度较强的地方(如旅游区、工程建设区)，无瓣海桑的生长及固碳效果会受到一定程度上的抑制；人类活动或干扰较小的地区，无瓣海桑的生长条件较好，固碳效果较强。

本文提供了无瓣海桑在固碳能力方面的基础数据，为红树林生态系统碳循环研究提供一定的参考。在全球气候变化背景下，在更大尺度下研究红树林固碳能力是下一步研究的重点。