南方红壤区是我国重要的林木产区。由于侵蚀性降雨集中、土壤抗蚀性差、林下植被匮乏^[1]以及强烈的人为干扰^[2], 南方红壤区已成为仅次于黄土高原的中国第二大水土流失区^[3]。长汀县是我国南方水土流失最严重的地区之一^[4], 该地区通过人工植树种草^[5]、封山育林^[6]等措施, 水土流失势头得到了初步控制。森林碳储量不仅能够维持生态平衡^[7]、反映森林质量^[8], 也对水土流失治理成效有一定的指示作用。因此, 森林碳储量可为区域生态环境评价、区域水土流失治理成效评估提供辅助^[9]。

马尾松(Pinus massoniana)作为长汀水土流失区生态恢复和重建的先锋树种^[10], 具有耐贫瘠、成活率高的特点^[11]。马尾松林地上林木碳储量的变化直接反映区域森林总碳储量的状况。传统的碳储量测量方法多采用样地清查法^[12], 该方法成本高, 工作量大, 仅适用于较小的区域, 并且人为选择样地, 主观因素强, 导致结果具有不确定性。针对大面积的碳储量估算, 遥感技术具有明显的优势。遥感技术可以对大面积的森林进行动态监测^[13], 更好地了解森林碳储量的变化。众多学者利用遥感技术对森林碳储量进行估算。Fauzi等^[14]对马来西亚吉兰丹热带森林不同海拔梯度碳储量进行估算。覃连欢^[15]使用生物量换算因子法估算广西省森林植被的碳储量。马尾松单木生物量模型是准确测算马尾松林碳储量不可缺少的工具, 对马尾松林生物量、碳储量的研究多采用福建省通用模型估算^[16-17]。但众多研究表明, 不同起源^[18]、不同生态环境^[19]、不同林龄^[20]的同一树种的生长存在较大差异。长汀水土流失区马尾松所处的生态环境比福建其他地区马尾松林特殊, 土壤极度贫瘠, 造成该区域马尾松生长状况较差, 同年生马尾松的胸径、树高通常达不到正常生长水平, 导致通用模型并不适用于该区马尾松生长预测。

基于此, 研究以长汀县河田镇马尾松林为对象, 利用项佳等^[21]提出的长汀红壤侵蚀区马尾松林生物量估算模型, 结合地面样地碳储量数据和河田镇2017年11月1日Landsat 8 OLI遥感影像, 构建马尾松林地上林木碳储量遥感反演模型。基于伪不变特征原理的线性归一化法((pseudo-invariant features, PIF)), 实现反演模型在2003年、2010年同期影像上的适用性校正, 反演3期马尾松林地上林木碳储量, 进一步分析区域马尾松林地上林木碳储量时空分异特征, 以期为我国南方水土流失区下一步生态恢复重建工作提供重要依据。

1 材料与方法 1.1 研究区概况

河田镇(东经116°16′—116°34′, 北纬25°30′—25°44′)位于福建省长汀县中部地区, 总面积约296 km² (图 1)。该地区属亚热带海洋性季风气候, 年均温17—19.5℃, 年降雨量1700 mm^[22], 降雨强度大。全镇低山高丘陵环绕, 属于河谷盆地, 平均海拔为390 m, 以山地丘陵红壤为主, 侵蚀严重, 土壤贫瘠、地表含沙量大, 属于强水土流失区。以河田镇为中心的红壤侵蚀劣地, 原有地带性植被常绿阔叶林被次生马尾松林代替, 稀疏矮小、生长缓慢的“小老头马尾松”众多, 林下分布以芒萁(Dicranopteris dichotoma)为主, 形成以马尾松-芒萁为主的植被群落结构^[19]。在大力开展封山育林、种植水保林草、改造低效林等工程措施的综合治理下, 河田镇水土流失治理取得初步成效, 但林分结构不稳定、地表径流频发的背景下其生态环境仍十分脆弱。

1.2 数据的获取与预处理 1.2.1 样地的设置

根据典型取样原则, 于2017年7月在全镇范围内随机但空间上又具有一定均匀性和代表性设置45个20 m×20 m的样地, 采用UG905高精度手持GPS准确定位。考虑到研究区马尾松林木林分质量相对较弱, 故本次调查起测胸径为3 cm。对45个样地内达到起测胸径的马尾松进行每木检尺, 测量并记录胸径、树高和冠幅。

1.2.2 样地碳储量的计算

基于马尾松林样地每木检尺的胸径、树高信息, 利用项佳等^[21]提出的长汀红壤侵蚀区马尾松林生物量估算模型(见表 1), 计算各样地蓄积量、生物量。

采用李海奎^[23]提出的马尾松含碳系数(0.4596)计算各样地地上林木碳储量值：

(1)

式中, C为平均碳储量(t/hm²)；B为平均生物量(t/hm²)；C_C为含碳系数。

1.2.3 遥感数据获取及预处理

利用遥感影像定量反演时, 理论上外业调研时间应与影像获取时间保持一致或接近。但区域可获取的2017年7月份左右有效卫星影像缺失, 考虑到本研究时间跨度长、马尾松生长指标相同年份内变化相对较小, 最终选定2017年11月1日的Landsat 8 OLI影像作为同期模型构建的基础数据, 2003年、2010年则选择的2003年12月29日、2010年1月14日Landsat 5 TM影像, 以上数据均来源于地理空间数据云(http://www.gscloud.cn/)。在此基础上采用2017年分辨率更高的Spot 6影像和Eco Drone-UAS八旋翼无人机获取的无人机影像作为参考影像, 对3期基础影像进行匹配校准, 配准精度在0.5个栅格以内, 并对影像进行辐射定标、大气校正、影像裁剪等预处理。

1.3 模型自变量的选择

研究区马尾松林下以芒萁等耐旱瘠草为主。曾宏达^[24]在对比马尾松与芒萁的光谱曲线后, 发现芒萁反射率高于马尾松, 特别是在绿光、红光等波段差异明显。而植被指数是光谱反射率的组合运算, 可以突出影像中的植被信息。为此, 研究选取植被指数(表 2)、红光波段、绿光波段参与建模。

多元回归常用于多变量模型的构建, 多变量的引入可以提高模型的决定系数(R²)。通过对各植被指数进行相关性分析后发现, 各植被指数间均存在极显著正相关, 如果直接进行多元回归可能造成模型的“估计不稳定”。主成分分析可在保留原指标信息的同时避免了多元变量的共线性问题。基于此, 为减少变量相关性对模型估测的影响, 研究通过主成分分析(Principal Component Analysis, PCA)选取主成分因子(Principal Component Factor, PCF)(式2)参与建模。

(2)

1.4 基于伪不变特征原理的线性归一化法

研究区模型是基于2017年11月1日Landsat 8 OLI影像构建的, 为更好的将模型反演于2003年12月9日、2010年1月14日Landsat 5 TM影像, 研究基于伪不变特征原理(PIF), 实现三期影像的数据匹配校正。伪不变特征原理(PIF)是从多幅影像中选取不变特征点, 提取这些点的灰度值, 建立参考影像与目标影像特征间的回归方程, 达到目标影像校正的目的。由于研究时序长, 且研究区属于山区, 2003年至2017年间不变地物特征点相对难寻, 故通过分析在确定不变地物点上采取间接的方法更加合适。

2 结果与分析 2.1 2017年马尾松林地上林木碳储量反演模型的确定 2.1.1 2017年马尾松林地上林木碳储量遥感估测模型的构建

利用公式(1)计算了45个样地的马尾松林地上林木碳储量。在对45个样地点进行初步回归拟合和残差分析的基础上, 去掉3个异常点。按照建模集样本80%的选取原则, 随机选取34个样本作为建模集, 8个样本作为检验集。

利用SPSS软件对34个建模样本进行指数回归模型的拟合, 采用决定系数(R²)和F统计量、显著性水平(Sig.)、估算标准误差等指标比较模型拟合效果。各模型及拟合指标见表 3。

由表 3可知, 除绿光波段反射率和红光波段反射率外, 其余植被指数建模Sig.值均为0.00, 模型具有显著性。其中以GNDVI、NRI、PCF为自变量拟合的回归模型R²大于0.5, F统计值分别为42.08、41.60、34.60, 估算标准误差分别为0.42、0.42、0.44。初步选择这三组指数模型为较优模型。

2.1.2 模型精度检验

采用均方根误差、平均相对精度两个指标对上述3种较优模型的估算精度进行检验。由图 2可知, 三组模型对马尾松林地上林木碳储量整体估算效果均较好。结合表 4可知, 以GNDVI构建的指数模型(C₂₀₁₇=0.006e^{14.357GNDVI₂₀₁₇})均方根误差较低、平均相对精度较高, 最终选定该模型为估算长汀县河田镇马尾松林地上林木碳储量的反演模型。

2.2 2003年、2010年模型的确定及3期林木碳储量的反演

基于伪不变特征原理的线性归一化法(PIF), 在提取三期不变地物点GNDVI值基础上, 以校正之后的2010年GNDVI影像校正2003年的GNDVI影像, 最终得到三期影像GNDVI值间的校正关系(表 5)。两个模型拟合R²均达到0.85以上, 整体表明拟合效果可以满足不同时期影像校正需求。采用均方根误差检验待校正影像(2003年、2010年)与参考影像的相似性关系, 检验结果见表 6。校正后影像GNDVI值的均方根误差均小于校正前。

将表 5中2017年GNDVI与2003年、2010年GNDVI值转换关系式分别代入构建的马尾松林地上林木碳储量遥感反演模型y=0.006e^14.357x中, 得到三个时期的长汀县河田镇马尾松林地上林木碳储量遥感反演模型见表 7。利用构建的遥感估算模型反演2003年、2010年及2017年长汀县河田镇马尾松林地上林木碳储量。

根据区域碳储量数值分布范围, 将碳储量划分为0—10 t/hm²、10—20 t/hm²、20—30 t/hm²、30—40 t/hm²、40—50 t/hm²、>50 t/hm²这6个等级进行统计, 得到马尾松林地上林木碳储量等级分布专题图(见图 3)。

2.3 马尾松林地上林木碳储量的时空分异特征 2.3.1 马尾松林地上林木碳储量的时间分异特征

河田镇马尾松林地上林木碳储量统计结果见图 4。2003、2010及2017年地上林木碳储量分别为8.24 t/hm²、11.34t/hm²、16.14 t/hm², 整体呈上升趋势。

分别对2003、2010及2017年马尾松林各等级地上林木碳储量像元数进行统计, 由表 8可知, 2003—2017年地上林木碳储量在0—10 t/hm²的林分面积不断减少, 地上林木碳储量大于10t/hm²的林分面积不断增加。2010年地上林木碳储量在0—10 t/hm²的区域比2003年减少了7.43个百分点；2017年地上林木碳储量在0—10 t/hm²的区域比2010年减少了27.84个百分点。2010年地上林木碳储量在10—20 t/hm²、20—30 t/hm²、30—40 t/hm²、40—50 t/hm²和>50 t/hm²的区域比2003年分别增加了3.45、1.68、0.71、0.34和1.25个百分点；2017年地上林木碳储量在10—20 t/hm²、20—30 t/hm²、30—40 t/hm²、40—50 t/hm²和>50 t/hm²的区域比2010年分别增加了17.28、6.05、2.77、1.07和0.67个百分点。

2.3.2 马尾松林地上林木碳储量的空间分异特征

由图 5可知, 2003年、2010年和2017年马尾松林地上林木碳储量在海拔 < 300m的区域碳储量最少、在海拔>500m的区域碳储量最大, 碳储量随海拔的升高而增加。2003—2010年, 各海拔梯度上地上林木碳储量增长率分别为107.82%、59.94%、27.92%、37.78%、26.08%、5.77%和-9.90%;2010—2017年, 各海拔梯度上地上林木碳储量增长率分别为123.22%、84.35%、34.26%、16.51%、20.44%、24.85%和21.16%, 整体表现出随海拔增加碳储量增长率降低的特征。

由图 6可知, 2003年、2010年和2017年马尾松林地上林木碳储量在坡度0—5°的区域碳储量最少、在坡度>45°的区域碳储量最大, 碳储量随坡度的增加而增加。2003—2010年间, 各坡度地上林木碳储量增长率分别为53.71%、46.09%、37.99%、25.94%、16.36%、14.56%;2010—2017年, 各坡度地上林木碳储量增长率分别为81.52%、70.27%、41.86%、28.51%、8.31%、4.27%, 表现出随坡度增加碳储量增长率降低的特征。

图 7表明, 向阳坡(67.5°—247.5°)碳储量高于背阴坡(0—67.5°, 247.5°—360°)。2003—2010年间, 向阳坡和背阴坡地上林木碳储量增长率分别为35.06%、53.96%;2010—2017年间, 向阳坡和背阴坡地上林木碳储量增长率分别为28.02%、119.39%, 表现出背阴坡地上林木碳储量增长率高于向阳坡。

3 结论与讨论

本研究以长汀县河田镇马尾松林为研究对象, 基于2017年Landsat 8 OLI遥感影像, 充分提取与碳储量潜在相关的信息, 构建马尾松林地上林木碳储量遥感反演模型, 并基于伪不变特征原理的线性归一化法(PIF)实现该遥感反演模型在2003年、2010年Landsat 5 TM影像上的适用性校正转换。实现对长汀县河田镇马尾松林地上林木碳储量的遥感估测和时空分异特征研究。

2017年遥感估算的平均碳储量值与研究布设的42个样地实测碳储量平均值接近, 表明研究构建的遥感模型对碳储量的估算结果具有可靠性。2003、2010、2017年长汀县河田镇马尾松林地上林木碳储量超过50 t/hm²的面积分别占3.33%、4.58%、5.26%, 表明了以河田镇为代表的南方水土流失区虽然森林覆盖度恢复显著, 但林木生长质量并不佳、森林碳汇功能的发挥仍较弱, 这与刘政等^[3]的研究结果大体一致。2010—2017年这7年间各级碳储量变化幅度大于前7年, 说明研究区在水土流失治理工作中采取的封山育林、种植水保林草及低效林改造等措施是成功的。

研究区马尾松林碳储量空间分布与海拔、坡度、坡向的空间分布密切相关。马尾松林碳储量随海拔升高而增加, 这一结果与黄绍霖^[30]的研究结果大体一致。海拔变化间接造成局地气候条件变化, 且随海拔增高人为干扰强度减小, 马尾松分布相对密集, 而低海拔区除人类活动频繁外, 地面原生植被覆盖较少, 多以补植的幼龄林为主, 林分固碳能力差。平均碳储量增长率随海拔升高而降低, 产生这种空间分异规律的原因可能是：低海拔区域生态治理强度大, 马尾松林生长快速, 高海拔区域多以封山育林为主, 且部分结构不合理、立地质量差的林分存在倒退的可能。马尾松林地上林木碳储量随着坡度增加而增加, 分析原因可能有两方面：一方面是坡面增加了林冠层接受光照的面积, 有利于马尾松林木的生长和吸收固定CO₂；另一方面, 与人类活动密不可分, 随坡度上升, 虽然土壤养分流失增多, 但森林采伐利用难度也增加, 人为干扰活动较少, 使得这部分马尾松林分生长健康稳定。碳储量增长率随坡度升高而降低, 主要原因可能是早期水土流失治理主要选择在森林覆盖较低、生长质量较差的缓坡区, 该部分区域补植造林强度较大, 随着树木生长到一定年龄, 生物量碳储量会显著增加, 而陡峭区域更多采取封山育林等保护政策, 自然更新能力低于人为治理下林木生长更新。向阳坡的马尾松地上林木碳储量高于背阴坡, 主要原因可能有两方面：一是河田镇位于北半球, 向阳坡受到的太阳辐射强度和日照时长远高于背阴坡, 且马尾松具有喜光喜温的特性, 光照充足的条件有利于马尾松生长和干物质的积累；二是在水、肥、气、热的影响下, 阳坡和半阳坡土壤养分和土壤碳储量的积累增强了马尾松生长能力, 同时, 山地不同坡向迎风、背风效果不同, 迎风坡容易造成地形雨, 局部降水分配强度不同也会影响马尾松碳储量的积累。背阴坡碳储量增长率远高于向阳坡, 产生这种背阴坡碳储量总值低但增长率高的原因可能和水土流失治理密不可分, 补植的马尾松处于生长初期, 胸径、树高在一定年龄内生长较快, 而随着年龄的增长, 马尾松生长速度会减缓至平稳状态。

本研究的研究对象为马尾松地上乔木层部分的碳储量, 未进一步考虑林下芒萁对碳储量的影响。在选择植被指数作为遥感反演模型自变量时, 虽然参考了前人研究中能有效区分芒萁和马尾松的光谱信息, 但尚未实现芒萁和马尾松的剥离, 而研究时段内林下芒萁覆盖度的变化, 可能会对模型反演的结果产生误差。因此, 可以考虑从像元尺度剥离芒萁的影响。