蒋馥根,孙 华,*,李成杰,马开森,陈 松,龙江平,任蓝翔
1 中南林业科技大学 林业遥感信息工程研究中心, 长沙 410004
2 林业遥感大数据与生态安全湖南省重点实验室, 长沙 410004
3 南方森林资源经营与监测国家林业与草原局重点实验室, 长沙 410004
4 广西壮族自治区森林资源与生态环境监测中心, 南宁 530000
5 国家林业和草原局中南调查规划设计院, 长沙 410004
森林是陆地生态系统的主体,也是生物圈中最大的碳储库,在全球碳循环和气候变化中起着至关重要的作用[1]。森林地上生物量(Aboveground Biomass,AGB)作为评价森林服务功能和森林质量的主要指标之一,越来越受到人们的重视[2]。准确地估测AGB是森林资源管理和生态系统动态监测的基础性工作[3]。与传统的人工调查相比,遥感技术具有快速获取大尺度植被生长状况的潜力,为森林资源的监测与管理提供了有效的方法参考[4- 5]。通过抽样调查的方式获取少量野外样地数据,结合遥感影像变量进行反演已成为区域范围内生物量估测的一种有效方式[6- 7]。
光谱反射率能体现地物之间的差异,是森林参数遥感反演的理论基础。光学遥感影像是世界上覆盖范围最广、类型最多、时间序列最丰富的数据[7]。光学遥感影像广泛应用于土地变化和森林资源监测,丰富的光谱信息能够有效反映植被的生长情况[8- 9]。然而,云层覆盖、数据条带缺失以及高植被覆盖引起的光谱饱和现象限制了光学遥感影像的进一步应用[6,10]。与传统光学遥感成像手段相比,合成孔径雷达(SAR,Synthetic Aperture Radar)和激光雷达(LiDAR,Light Detection and Ranging)可以穿透植被冠层以获得林分或森林详细的垂直结构信息,理论上可以获得更为准确的森林参数估测结果[11- 14]。然而,SAR必须在特定频率的波段中工作,Gleason等[7]发现这些波段不一定适合生物量估计。LiDAR利用冠层高度模型(CHM,Canopy Height Model)获取研究区域的平均树高,通过构建估测模型实现对整个区域的森林生物量估算[15]。但大范围的高分辨率数字高程模型(DEM,Digital Elevation Model)是获取准确的CHM的基础,这对于地形复杂的山区来说无疑是困难的[16]。此外,主动式遥感数据成本较高,目前大区域内没有完整的高质量周期性数据覆盖,不利于开展大范围的森林监测[10]。
目前,中分辨率光学数据仍然是大面积人工林AGB估测最常用的数据源[7]。Sentinel-2由于其携带的红边波段对植被生长状况较敏感,已被广泛应用于森林AGB估测[17-18]。GF-6作为携带红边波段的第一颗国产光学卫星,具有大幅提高对农业、林业、草原等植被资源监测的潜力[19-21]。梁继等[19]研究了红边特征在农作物识别中的表现,结果表明通过红边波段构建的红边归一化植被指数在作物样本类别区分度上比传统NDVI(Normalized Difference Vegetation Index)更显著,GF-6红边特征波段在作物识别中表现优越,使得识别精度显著提高。Yu等[21]证明了GF-6两个红边波段对于土地覆盖的变化检测可以达到89%的整体精度,与使用同期GF-1数据相比有显著提升。目前,GF-6在森林资源定量估测上的研究较少,对于森林AGB估测的有效性仍需进一步验证。此外,单个光学数据源由于其光谱信息的限制,应用于生物量反演的效果始终受到限制。联合多个数据源进行森林地上生物量反演的方式越来越受到欢迎[5,13]。
本研究以旺业甸林场和塞罕坝林场为研究区,基于GF-6和Sentinel-2影像分别提取遥感变量,探索联合数据源以及红边波段对AGB估测的效应。同时,结合野外实测样地数据构建落叶松和樟子松AGB线性和非线性估测模型,比较GF-6、Sentinel-2、联合GF-6和Sentinel-2三种数据源在AGB估测的效果,选择估测效果最优的模型进行AGB反演和空间分布制图。通过对最优模型进行适用性评价,以期为人工林AGB遥感估测提供方法参考。
旺业甸林场位于内蒙古自治区赤峰市西南部(118°09′—118°30′E,41°21′—41°39′N),总面积259.58 km2(图1)。林场属中温带大陆性季风气候,海拔范围为500—1890 m,温度和降水随海拔升高有明显变化。年均气温4.2℃,无霜期117 d,年平均日照2913.3 h,年平均降水量为400 mm。林场森林覆盖率约为93%,主要树种为落叶松(Larixgmelinii(Rupr.)Kuzen.)、油松(PinustabuliformisCarr.)、樟子松(Pinussylvestris)和云杉(PiceaasperataMast.)等。落叶松是林场主要的木材产品,也是北方地区重要的生态资源,及时掌握落叶松生物量分布和动态变化对于了解中国北方地区森林生态系统碳循环具有重要的价值。
塞罕坝林场成立于1962年,总面积27300 hm2,是全国最大的人工林林场。林场位于河北省承德市北部(116°51′—117°39′E,42°02′—42°36′N)(图1)。年平均气温为-1.4℃,冬季漫长,生长季节短,平均海拔为1700 m。塞罕坝属寒温带大陆性季风气候,极端最低气温-43.2℃,年平均积雪7个月。年平均降水量约490 mm,年平均无霜期68 d。林场目前森林覆盖率已达80%,林场树种以樟子松(Pinussylvestris)和落叶松(Larixgmelinii(Rupr.)Kuzen.)为主,林木总蓄积量达到1012万m3。
图1 研究区位置示意图Fig.1 Location of the study area
1.2.1野外AGB数据获取及处理
利用eCognition 8.7分类软件对覆盖旺业甸林场的GF-6 PMS(数据获取时间为2019年9月30日)影像(2 m分辨率)进行多尺度分割,结合2018年旺业甸林场植被分布图提取植被类型,最终确定林场内落叶松分布区域(图2)。根据不同立地条件和生长状况在林场内落叶松分布区域共设置42个25 m×25 m的样地,并于2019年9至10月进行野外实测调查。选用Trimble Geo 7x GPS全球定位仪记录样地中心坐标。对样地内胸径大于5 cm的样木进行每木检尺,测量树高和胸径数据,同时记录样地海拔、坡度、坡向等信息。利用塞罕坝林场2018年森林资源调查数据库,提取樟子松边界后计算、统计所有小班蓄积量结果,并在边界内随机设置100个25 m×25 m的样地,最终确定所有样地的森林蓄积量(图2)。
图2 落叶松和樟子松边界及样地分布Fig.2 Boundary and plot distribution of Larch and Scotch pine
在全国标准信息公共服务平台(http://std.samr.gov.cn/)获取2018年内蒙古自治区人工落叶松二元立木材积公式,并计算所有样地内落叶松单木材积,累加得到样地蓄积量(m3)后换算成每公顷蓄积量(m3/hm2)。樟子松单位面积蓄积量通过小班蓄积量值直接转换。利用生物量转换因子连续函数法[22]推算出每公顷地上生物量,参数a和b采用李海奎《中国森林植被生物量和碳储量评估》[23]中落叶松和樟子松树种/树种组参数统计结果,落叶松取0.6096和33.8060,樟子松取1.0945和2.0040。
V1=0.002466-0.000014×D2+0.001142×H+0.000036×H×D2
(1)
式中,V1为落叶松材积(m3),D为胸径(cm),H为树高(m)。
B=aV2+b
(2)
式中,B为每公顷生物量(t/hm2),V2为每公顷蓄积量(m3/hm2),a和b为常数。
落叶松和樟子松AGB分布范围分别为60.72—269.87 t/hm2和30.06—414.23 t/hm2(表1)。旺业甸林场落叶松样本标准差和变异系数分别为54.60 t/hm2和37.1%,说明样本呈一定程度的离散分布。塞罕坝林场樟子松样本标准差和变异系数分别为92.50 t/hm2和92.8%,样本变异程度较大。落叶松和樟子松AGB均值分别为147.11 t/hm2和99.65 t/hm2。
表1 研究区森林生物量统计结果Table 1 Forest aboveground biomass statistical results of the study area
1.2.2遥感影像数据处理
本研究获取的遥感影像数据包括覆盖研究区且与野外调查时间相近的Sentinel-2和GF-6影像。Sentinel-2卫星携带一枚多光谱成像仪(Multispectral imager,MSI),卫星重访周期为10 d,幅宽为290 km[17]。Sentinel-2是唯一拥有三个以上红边波段的光学卫星,这对准确获取植被生长状况非常有效[17]。为了获取更准确的植被信息,研究选择空间分辨率较高的Band 2—Band 8A,共计8个波段。GF-6是我国第一颗携带红边波段信息的光学遥感卫星,携带WFV(Wide Field of View)和PMS(Panchromatic and Multi-spectral)两个传感器,能够有效反映森林和植被生长特性(表2)[19]。GF-6幅宽达800 km,与其他高分系列卫星组网运行后重返时间减少至2天,能持续、高效地提供高分辨率影像,大幅提高对森林、草原等植被的监测能力[20]。
表2 研究采用的Sentinel-2和GF-6波段信息Table 2 The selected bands of Sentinel-2 and GF-6 in the study
利用Sen2cor 2.5.5和ENVI 5.3遥感图像处理软件分别对Sentinel-2和GF-6影像进行辐射定标、大气校正等预处理,最终将像元亮度值转换为反射率值[17-18]。为了将样地与遥感数据进行匹配,获得更准确的植被信息,将两种影像像元大小均重采样至与样地大小一致。
1.2.3遥感反演变量提取及筛选
光谱反射率和植被指数是建立森林地上生物量遥感反演模型的重要变量[10]。植被指数能对植被生长状况进行有效地度量,合适的植被指数参与建模能显著提高AGB估测精度[18]。作为常用的植被指数,归一化植被指数(Normalized Difference Vegetation Index,NDVI)[24]、增强型植被指数(Enhanced Vegetation Index,EVI)[25]、大气抗阻植被指数(Atmospherically Resistant Vegetation Index,ARVI)[25]、可见光抗大气指数(Visible Atmospheric Resistant Index,VARI)[25]和红绿植被指数(Red-Green Vegetation Index,RGVI)[25]与植物生长周期中重要的生理参数之间有良好的相关关系,已被证明对森林AGB估测较敏感。红边波段处于近红外与红光波段之间快速变化的区域,对植被叶绿素的微小变化较敏感[26]。研究将同时构建红边归一化植被指数(Red-Edge Normalized Difference Vegetation Index,RENDVI)[26]和红边叶绿素指数(Red Edge Chlorophyll Index,RECI)[27]作为遥感变量参与建模,以评价红边植被指数对于AGB估测的效应(表3)。
表3 植被指数计算公式Table 3 Formulas of the vegetation index
研究同时建立线性和非线性回归模型用于研究区森林AGB估测和模型精度比较。为了验证变量数量对于模型估测的影响,建立单变量线性模型、单变量非线性模型(三次多项式)、多元线性回归[28]和多变量非参数模型(随机森林)[29]用于比较模型估测效果。相比单变量回归,多元线性回归模型能定量地、更全面地衡量多个遥感变量对于AGB估测的影响。逐步回归分析用于检测并保留与AGB显著相关的变量,同时结合方差膨胀因子VIF(Variance Inflation Factor)指数消除变量间的共线性,最终保留的变量应当与AGB显著相关且变量之间相互独立[24]。此外,非线性模型中的非参数模型对目标函数的形式没有提前假设,更适合于复杂数据的预测[29]。随机森林作为非参数模型的代表,已成为AGB估测的常用方法[8,10,17]。随机森林通过快速构造大量的决策树来估计AGB,模型训练过程中,决策树之间相互独立,训练速度快。另外,随机森林可以评估模型中每个变量的重要性,从而有效地判断变量对模型的贡献[17,29]。研究将通过评价变量的重要性,通过重要性排序筛选出对模型贡献度较大的变量组合进行AGB估测。
为了验证GF-6和Sentinel-2影像对AGB估测的有效性,研究分为Sentinel-2变量组、GF-6变量组以及联合GF-6和Sentinel-2变量组,分别建立回归模型进行AGB估测。同时,为了充分利用样本以提高模型的可靠性,研究选择留一交叉法[25]对估测结果进行验证,即每次只留下1个样地做测试集,其他样地做训练集,重复n次后确定最终估测结果。选择决定系数(R2)[30]和均方根误差(RMSE)[30]进行模型评价,计算方式如下:
(3)
(4)
计算所有遥感变量与森林AGB之间的Pearson相关系数,并将相关系数矩阵可视化。由图3可知,从GF-6和Sentinel-2影像中提取的遥感变量中,光谱反射率与落叶松AGB存在较高的相关关系,且整体相关系数高于植被指数。其中,红边波段反射率均与AGB显著相关(P<0.05)且相关性系数较高。GF-6和Sentinel-2影像中相关系数较大的变量主要由红边波段、近红外波段、蓝色波段和黄色波段组成。与落叶松不同的是,樟子松植被指数相关系数整体大于光谱反射率(图4)。以GF-6为数据源时,与樟子松AGB相关性最高的变量为RECI,表现为极显著正相关(P<0.01)。
图3 Sentinel-2与GF-6遥感变量与落叶松AGB相关系数矩阵图Fig.3 Matrix diagram of correlation coefficient between remote sensing variables and AGB of Larch of Sentinel-2 and GF-6AGB:地上生物量Aboveground Biomass;ARVI: 大气抗阻植被指数Atmospherically Resistant Vegetation Index;EVI: 增强型植被指数Enhanced Vegetation Index;NDVI: 归一化植被指数Normalized Difference Vegetation Index;RGVI:红绿植被指数Red Green Vegetation Index;VARI: 可见光抗大气指数Visible Atmospheric Resistant Index;RENDVI: 红边归一化植被指数Red Edge Normalized Difference Vegetation Index;RECI: 红边叶绿素指数Red Edge Chlorophyll Index
图4 Sentinel-2与GF-6遥感变量与樟子松AGB相关系数矩阵图Fig.4 Matrix diagram of correlation coefficient between remote sensing variables and AGB of Scotch pine of Sentinel-2 and GF-6
利用R语言RandomForest函数包分别计算变量的重要性,按照重要性从高到低的顺序依次增加变量分别构建随机森林模型,分析不同变量数量下模型RMSE的变化趋势。由图5和图6可知,对于GF-6、Sentinel-2以及联合GF-6和Sentinel-2三类变量,当变量数量分别取3、10、2和5、10、10时,误差达到最小,将此时的变量组合最终用于AGB估测。
图5 旺业甸林场不同变量组的重要性排序和RMSE变化图Fig.5 Importance ranking under different variable group and RMSE change in Wangyedian forest farmRMSE: 均方根误差Root Mean Square Error
图6 塞罕坝林场不同变量组的重要性排序和RMSE变化图Fig.6 Importance ranking under different variable group and RMSE change in Saihanba forest farm
以GF-6和Sentinel-2作为数据源,结合旺业甸林场和塞罕坝林场AGB实测数据建立单变量线性模型、单变量三次多项式模型、多元线性回归和随机森林模型进行AGB估测和模型精度比较。结果如表4所示,单变量线性模型和三次多项式模型估测效果较差,不能有效地进行AGB估测。多变量模型中,多元线性回归在三组数据源中均优于随机森林模型,在落叶松和樟子松构建的模型中RMSE分别降低了8.0%、19.3%、28.5%和11.6%、5.0%、4.4%。
表4 AGB回归估测结果Table 4 AGB regression estimation results
利用逐步回归分别对三组变量进行变量筛选,最终获得的变量主要由红边波段、蓝色波段、近红外波段和黄色波段构成,说明这些波段信息对AGB具有较好的敏感性。相比Sentinel-2,基于GF-6建立的多元线性回归模型决定系数R2有一定提高,估测误差略有下降,估测效果较优。联合GF-6和Sentinel-2影像构建的多元线性回归模型估测效果最好,决定系数R2分别达到了0.66和0.65,RMSE为31.45 t/hm2和54.77 t/hm2。多元线性回归模型在所有模型中估测效果最好,将作为最优模型进行旺业甸林场落叶松和塞罕坝林场樟子松AGB反演和空间分布制图。
为了验证红边波段对落叶松、樟子松AGB估测的有效性,将三组数据源剔除红边波段信息之后,结合多元线性回归再次建立估测模型进行比较。结果如表5所示,在剔除红边波段之后,三组数据源在估测落叶松AGB模型精度显著下降,模型RMSE分别提高了9.6%、17.6%和28.7%,樟子松分别提高了6.2%、7.9%和8.2%,说明加入红边波段能显著提高AGB估测精度。此外,在剔除三组数据源的红边波段信息后,联合GF-6和Sentinel-2的数据源相比GF-6、Sentinel-2估测效果仍然有显著提高,RMSE显著降低。
表5 剔除红边波段后多元线性回归估测结果Table 5 Multiple linear regression estimation results after eliminating red edge band
通过Sentinel-2和GF-6影像分别构建的多元线性回归模型均存在一定程度的高估和低估,但将两种影像联合后建立的回归模型误差显著降低,拟合效果达到最优(图7、图8)。联合GF-6和Sentinel-2建立的估测模型预测AGB值与实测AGB值具有良好的线性关系,绝大部分预测值随机分布在回归线两侧。
图7 基于Sentinel-2、GF-6和联合数据源的落叶松AGB回归模型拟合图Fig.7 Fitting diagram of AGB regression models of Larch based on Sentinel-2, GF-6 and federated data source
图8 基于Sentinel-2、GF-6和联合数据源的樟子松AGB回归模型拟合图Fig.8 Fitting diagram of AGB regression models of Scotch pine based on Sentinel-2, GF-6 and federated data source
基于Sentinel-2、GF-6、联合GF-6和Sentinel-2筛选的变量通过构建多元线性回归分别得到了旺业甸林场落叶松AGB的空间分布图(图9)。图9所示的AGB空间分布相似。在三种AGB空间分布中,落叶松AGB均值分别为136.8 t/hm2、128.7 t/hm2和140.98 t/hm2,与使用单个数据源相比,联合GF-6和Sentinel-2数据源得到的AGB均值有一定的增加,且最接近实测样地的AGB均值147.11 t/hm2。较大的AGB预测值主要分布在西部和东南地区,而北部和西南地区的预测值相对较小。由于靠近建筑物和道路,人为干扰活动较多,中部和西北地区的AGB值偏小。对于塞罕坝林场,生物量值较大的值主要分布在西部和西北部地区,联合数据源构建的模型预测效果较好(图10)。以联合GF-6和Sentinel-2为数据源预测的AGB空间分布与林场实际AGB分布基本一致,可以作为有效的数据源进行AGB反演。
图9 旺业甸林场落叶松AGB空间分布图Fig.9 AGB Spatial distribution of Larch in Wangyedian forest farm
图10 塞罕坝林场樟子松AGB空间分布图Fig.10 AGB Spatial distribution of Scotch pine in Saihanba forest farm
研究通过建模变量与最优模型产生的残差之间的显著性检验以及残差的空间自相关检验进行模型的适用性评价。对于落叶松和樟子松,所有建模变量与残差之间均不存在显著相关(P>0.1)。以落叶松为例,构建AGB残差与建模变量之间的散点图(图11),结果表明变量对AGB估测的不确定性没有显著影响,变量筛选方法和结果是可行的。此外,残差之间的空间自相关通过莫兰指数(Moran′s I)来检验,最终落叶松和樟子松的结果分别为0.21和0.05(P>0.1),表明残差之间没有显著的空间自相关现象,呈随机分布。
图11 残差与建模变量散点图Fig.11 Scatter plots of residuals and modeling variables
光谱反射率和植被指数与森林生长状况密切相关,可以作为森林AGB估测基本的建模变量[10,17]。相比可见光波段,红边波段能对植被冠层结构和叶绿素含量的微小变化做出快速反应,有效提高AGB估测精度[18- 20]。同时,基于红边波段构建的植被指数可以减弱NDVI等传统植被指数的饱和现象,从而更准确地预测和评价植被的生长状况[26]。由旺业甸林场和塞罕坝林场GF-6和Sentinel-2提取的遥感变量中,红边波段的反射率均与AGB显著相关,且经过变量筛选之后保留的变量中也均包含红边波段反射率及红边波段构建的植被指数。联合GF-6和Sentinel-2建立的多元线性回归模型取得了最高的决定系数和最低的RMSE,相比单个数据源,模型RMSE分别得到了显著降低。Li等[31]以资源三号高分辨率影像为数据源结合多元线性回归实现了旺业甸落叶松AGB估测,模型RMSE相比本研究中的Sentinel-2和GF-6数据源分别降低了17.0%和13.6%,然而GF-6在联合Sentinel-2后建立的模型RMSE相比资源三号数据源下降了7.2%,说明联合GF-6和Sentinel-2影像结合多元线性回归能有效地对落叶松AGB进行估测。
相比Sentinel-2影像,GF-6拥有的黄色波段(590—630 nm)在植被生长监测中也展现了优越性。胡秀娟等[32]通过WorldView- 2卫星影像的黄色波段构建了新型植被指数,有效地对水土流失区的植被健康状况进行了监测。曾庆伟等[33]证明了基于GF-6近红外和黄边波段建立的YNDVI指数比传统的NDVI能更有效地反映森林扰动信息。但黄色波段对森林AGB估测的敏感性仍需进一步验证,研究发现GF-6黄色波段与旺业甸林场落叶松AGB 显著相关,且GF-6、联合GF-6和Sentinel-2变量组在多元线性回归中筛选的结果均保留了黄色波段,说明黄色波段对AGB估测较敏感。
旺业甸林场和塞罕坝林场均为人工林林场,森林状况分布较均匀,光学遥感获取的冠层光谱差异有限,这导致通过构建大量决策树进行预测的适用于更复杂情况的非线性随机森林模型估测效果有限。此外,样地数量和特征变量数量能显著影响模型的估测效果,对于小样本人工林,多元线性回归能取得比非线性模型更好的估测效果,这与Li等[31]在旺业甸林场人工林得到的研究结果相似。由于森林类型、树种、林龄的差异,遥感影像提取的地表植被信息存在同物异谱、同谱异物和其他能显著影响估测精度的因素。常用的遥感变量包括波段反射率、植被指数、地形因子、纹理指数等,这些变量与不同树种AGB之间的相关性及显著性并不一致。所以在进行AGB建模和估测时,需要同时考虑变量与AGB之间的线性和非线性关系,利用不同的特征变量选择方法构建相对应的模型探索更多的估测方式以提高估测精度和模型的可解释性。
SAR和LiDAR等主动遥感方式能穿透森林冠层,获取更准确的森林结构信息,从而提高AGB估测精度。但由于地形影响和数据覆盖有限等限制,利用主动遥感数据源进行大范围AGB估测效率有限[15]。此外,由于高分辨率光学影像数据量大、覆盖不完全、数据昂贵等因素,目前与国家森林资源连续清查样地大小基本一致的中分辨率光学数据仍是森林AGB反演的主流数据源。研究利用GF-6和Sentinel-2影像结合多元线性回归实现了旺业甸林场落叶松和塞罕坝林场樟子松AGB估测,能对森林AGB遥感监测研究提供一定的参考。然而,本研究所用的野外观测样地数量相对偏少且树种单一,采用线性回归所建立的生物量模型的适用性需要进一步验证。在今后的研究中,可以考虑通过增加实测样地数量和其他研究区对比进行方法验证。
以Sentinel-2和GF-6作为数据源,结合野外实测AGB数据建立线性和非线性回归模型,对旺业甸林场落叶松和塞罕坝林场樟子松AGB进行反演和空间分布制图。结果表明:(1)由Sentinel-2和GF-6影像中提取的红边波段信息均与AGB显著相关,且经过筛选用于建模的变量主要由红边、蓝色、近红外和黄色波段等对AGB较敏感的波段信息构成;(2)GF-6相比Sentinel-2拥有对植被参数较敏感的黄色波段,且作为数据源建立的多元线性回归模型RMSE相比Sentinel-2显著降低,表明国产GF-6可以代替Sentinel-2用于AGB估测。GF-6幅宽达到800 km,能更高效地监测大范围森林动态变化;(3)增加红边波段进行落叶AGB估测能显著提高模型估测精度。三组数据源分别加入红边波段信息后进行建模,落叶松估测模型RMSE分别下降了9.6%,17.6%和28.7%,樟子松则分别下降了6.2%、7.9%和8.2%;(4)联合GF-6和Sentinel-2建立的多元线性回归模型取得了最高的决定系数和最低的RMSE。相比单个数据源,模型RMSE得到了显著降低。联合GF-6和Sentinel-2数据源结合多元线性回归能有效地对森林AGB进行估测,为森林资源管理和遥感动态监测提供新的数据源参考。