在林业碳汇量估算中遥感方法的应用分析

Application Analysis of Remote Sensing Method in Forestry Carbon Sink Estimation

摘要: 研究针对没有样地数据的林业碳汇量估算，将遥感手段引入其中，作为现行林业碳汇量估算方式的补充。方法选择参考森林碳汇量估算方法，结合林业碳汇量特征，加以改进使之更符合用林业碳汇量估算。以三种方法对研究区碳储量进行估算，生物量方程法的估算结果最小值为1.97 MgC到最大值7.91 MgC；光能利用率法的估算结果最小值为0.57 MgC到最大值9.75 MgC；NDVI反演法的估算结果最小值为1.36 MgC到最大值8.68 MgC。光能利用率法估算的碳储量空间分布为由北向南、由西向东逐渐递增，NDVI反演法估算的碳储量空间分布为由南向北、由东向西逐渐递增。将生物量方程法的估算结果作为真实数据，分别用光能利用率法和NDVI反演法与真实数据进行结果拟合，NDVI反演法的最终拟合效果最好，并且碳储量空间分布也更符合实际。

Abstract: In this study, for the estimation of forest carbon sink without sample plot data, remote sensing is introduced into the estimation of forest carbon sink as a supplement to the current estimation method of forest carbon sink. Methods The reference forest carbon sink estimation method was selected, combined with the characteristics of forest carbon sink, and improved to make it more in line with the estimation of forest carbon sink. Three methods were used to estimate the carbon storage in the study area. The minimum value of biomass equation method was 1.97 MgC to the maximum value of 7.91 MgC. The minimum value of the light energy utilization method is 0.57 MgC to the maximum value of 9.75 MgC; the minimum value of NDVI inversion method is 1.36 MgC to the maximum value of 8.68 MgC. The spatial distribution of carbon storage estimated by the light energy utilization method is gradually increasing from north to south and from west to east. The spatial distribution of carbon storage estimated by NDVI inversion method is gradually increasing from south to north and from east to west. The estimation results of biomass equation method were used as real data, and the results were fitted by light energy utilization method and NDVI inversion method with real data. The final fitting effect of NDVI inversion method was the best, and the spatial distribution of carbon storage was more realistic. 

1. 引言

随着双碳理论的提出，林业碳汇作为碳中和的重要环节，其计算方法和适用条件也在逐渐发展完善。林业碳汇是指通过人为方式，增强林分质量加强碳汇能力，使之可以吸收更多的温室气体 [1]，助力完成2060年碳中和目标。林业碳汇是区别于森林碳汇的，森林碳汇是指森林生态系统在碳循环过程中固定多少碳，整个森林生态系统中的乔木、灌木、草本植物、微生物和无机环境等均参与其中。而林业碳汇只算增量不算存量，只计算能通过间伐修枝、补植补造等抚育手段人为干预的那部分森林碳汇，所以林业碳汇的计算范围是小于森林碳汇的。现阶段国家认定的林业碳汇计算方法只有生物量方程法和生物量扩展因子法 [2]，这两种方法都需要样地实际的测树因子，然而并不是所有的林业部门都设置样地，并且在北方地区每年只有五月到十月适合到样地进行数据调查，不管是碳汇计算还是后续的监测，通过样地来获得实际的测树因子会消耗大量人力物力。我们可以类比森林碳汇量计算方法，找出符合林业碳汇特征，成本较低而精准度较高的方法，来计算林业碳汇量。

现有的森林碳汇量计算方法有很多，其中包括样地清查法、微气象学法、箱式法、模型模拟法和遥感法五大类 [3]，其中样地清查法即为现行的林业碳汇计算方法，优点是精度高，但需要设置长期样地监测；微气象学法和箱式法需要仪器设备长期监测，时间和经济成本较高；模型模拟法通过模拟气候、植被和土壤在碳循环过程中的动态作用过程，进而计算碳汇量，缺点是影响因素过多，需要考虑温度降水日照等等环境因素，无论从数据获取还是数据处理方面，难度都极大。遥感法对比传统方法，具有直观、快速、对环境无破坏以及没有时空限制等优点，但由于遥感仅能监测地上范围，所以利用遥感法计算出的森林碳汇仅局限在地上碳储量，而林业碳汇主要计算的是乔木碳储量，故遥感法在林业碳汇计算和监测中应用潜力很大。

本研究旨在将遥感手段应用在林业碳汇计算中，包括纯遥感方法和遥感与地面信息结合的方法这两种，纯遥感方法这里选取了光能利用率法，遥感与地面信息结合的方法这里选取了植被归一化指数(NDVI)反演法。本文以黑龙江省海林林业局为研究区利用遥感数据和研究区2018年的林地二类调查数据，分别用只需要地面信息的生物量方程法、纯遥感的光能利用率法、遥感与地面信息结合的NDVI反演法计算研究区全域的乔木碳储量，从碳储量的空间分布上探究三种方法的精准度，并随机选取110个样本点，对这110个样本点的碳储量进行精度分析，进而对比出性价比最高的林业碳汇计算方法。

2. 材料与方法

2.1. 研究区概况

海林林业局位于长白山脉张广才岭东坡，黑龙江省海林县的中部。地理坐标为东经128˚42'~129˚26'，北纬44˚36'~45˚02'。北部和东部与柴河林业局毗连，东部和东南部分别同牡丹江市郊区林业局、海林市林业局相接：南部与大海林林业局和山市种奶牛场为邻；西部以张广才岭主脉大平岭山脊与亚布力林业局分界，林区面积156,621公顷。地貌特征属山区和丘陵浅山区，地貌特征为“九山半水半分田”山峰连绵，地形复杂。地处中温带东部大陆性季风气候区。属寒温带大陆性季风气候，四季分明，气候宜人，无租期短，为85~130天。春秋季短、气候多变，夏季高温多雨，冬季漫长而寒冷，年平均气温42℃。土壤多为暗棕壤。海林林业局域内地表水系发育主要为横道河子水系和头道河子水系。横道河子全长110千米，流域面积71,766公顷，其主要支流为大荒沟和大西沟。头道河全长77 Km，流域面积62,682公顷，头道河大小支流十几条，其中较大的支流为大石沟和夹皮沟。此外，区内较大河流还有密江河，流域面积2150公顷。具体位置如图1所示。

Figure 1. Study location map

图1. 研究区位置图

2.2. 数据来源

本研究包所需要的数据包括林地数据、遥感数据、气象数据。其中，林地资源二类调查数据来自当地林业局；遥感数据包括来自地理国情监测云平台landsat8卫星遥感影像和多种卫星遥感数据反演光合有效辐射数据、来自NASA的包含气溶胶和水汽数据的MOD08-D3遥感数据、来自中国气象数据网的气象数据(主要以温度为主)。所有数据年份均为2018年。

2.3. 研究方法

1) 生物量方程法

生物量方程法是一种相对准确地估算乔木生物量的方法，主要通过建立乔木的胸径或树高等测树因子与生物量间的关系，以测树因子为变量，估算乔木整体生物量 [4]。由于研究区没有样地数据，所以根据二类调查数据，将各林班小班的优势树种统计出来，共分为红松、水曲柳、云杉、樟子松等12类树种，并统计各林班小班对应的平均胸径，根据代海军等 [5] - [10] 的研究，通过生物量与胸径间关系，估算出单株乔木生物量。由于树木生长受立地条件制约，故不同研究区内的不同地区活立木密度也不同，为了最终估算的结果更加精准，故再根据公顷株树，求出公顷生物量，将最终的公顷生物量与各树种对应含碳量相乘，得到公顷碳储量。为了不同方法最终结果进行精度对比时，各方法取样点具有可比性，故按照遥感数据像元格面积与一公顷做比，得到一个像元格内的单位面积碳储量，遥感数据分辨率均选择30 m。

单位面积碳储量计算方法如下：

$C=B\times CF\times N\times k$

式中：

C——单位面积碳储量；

B——单株乔木生物量；

CF——树种含碳率；

N——公顷株数；

k——像元格面积与一公顷的比值(0.09)。

12类树种具体生物量方程及含碳率公式见表1。

Table 1. Biomass equation and carbon content of tree species

表1. 树种生物量方程及含碳率

注：B-单株乔木生物量；D-优势树种平均胸径。

2) 光能利用率法

根据Patil [11] 等人的研究，用植被归一化指数、光能利用率和光合有效辐射的乘积得出地上生物量。其中植被归一化指数在遥感监测植被特征，例如地上生物量、叶绿素含量等有广泛的用途 [12]，本研究利用Landsat8卫星遥感数据中近红外波段、远红外波段和红外波段来计算研究区的植被归一化指数；光能利用率是估算植被生产力的关键参数，光能利用率的计算包括CASA模型、MODIS17计算方法、VPM模型、EC-LUE模型等 [13]，本研究选用参数更少精准度更高的MODIS17的计算方法 [14]；光合有效辐射是影响光合作用的关键因子，有助于碳循环和碳驱动机制的研究，其敏感性对全球气候系统有着重要的影响，而且它在不同的陆地生态系统模型中，都是重要的输入参数。这里直接使用地理国情监测云平台推出的气象或气候环境类系列数据，多种卫星遥感数据反演光合有效辐射(PAR)产品。

由上述过程得到研究区的地上生物量，根据树种的地下生物量与地上生物量的比值，得到研究区整体乔木生物量，再与树种含碳率相乘，最终得到研究区乔木碳储量，由于本方法为纯遥感法，故假定无地面数据，取针阔混交林的地下生物量与地上生物量的比值和含碳率，分别为0.248和0.498 [2]。

地上生物量计算公式如下：

$AGB=NDVI\times PAR\times LUE$

式中：

AGB——地上生物量；

NDVI——植被归一化指数；

PAR——光合有效辐射；

LUE——光能利用率。

植被归一化指数计算公式如下：

$NDVI=\frac{NIR-Red}{NIR+Red}$

式中：

NDVI——植被归一化指数；

NIR——近红外波段；

Red——红外波段。

光能利用率计算公式如下：

$\epsilon ={\epsilon }_{max}\times T_{MIN\text{\_}scalar}\times VPD\text{\_}scalar$

$T_{MIN\text{\_}scalar}=\frac{T_{min}-T_{min\text{\_}MIN}}{T_{min\text{\_}MAX}-T_{min\text{\_}MIN}}$

$VP{D}_{scalar}=\frac{VPD-VP{D}_{max}}{VP{D}_{min}-VP{D}_{max}}$

式中：

$\epsilon$ ——光能利用率；

${\epsilon }_{max}$ ——植被最大辐射转换效率，取值0.389g C·MJ⁻¹ [15] [16] [17]；

$T_{min\text{\_}MAX}$ —— $\epsilon ={\epsilon }_{max}$ 时的每日最低温度(适用于最佳VPD)；

$T_{min\text{\_}MIN}$ —— $\epsilon =0$ 时的每日最低温度(适用于任何VPD)；

$VP{D}_{max}$ —— $\epsilon ={\epsilon }_{max}$ 时的昼夜平均水汽压差(适用于最佳T_min)；

$VP{D}_{min}$ —— $\epsilon =0$ 时的昼夜平均水汽压差(适用于任何T_min)。

乔木碳储量计算公式如下：

$C=\left[AGB\times \left(1+R\right)\times CF\right]$

式中：

C——乔木碳储量；

AGB——地上生物量；

R——地下生物量/地上生物量；

CF——树种含碳率。

3) NDVI反演法

现在应用比较广泛的NDVI反演法，是将大气校正过的植被归一化指数(NDVI)，通过回归方法与地面数据相结合 [18]，反演出研究区的乔木碳储量。相较于样地数据反推出的区域碳储量而言，其特点是当研究区内树种单一，混交方式单一时，仅需对几个实测点进行碳储量测算，再与这几个实测点对应的NDVI建立回归方程，通过NDVI反演出研究区的碳储量，避免了大量人力物力资金的投入的同时，还能生成碳储量空间位置分布图，更加直观地观测到碳储量分布情况，但是当研究区树种繁杂时，为了保证回归方程的精准性，还是避免不了对不同树种搭配方式下、不同混交方式下大量的样地调查。

由于本研究不涉及样地实测数据，故对该方法进行了调整，首先使用6s模型将NDVI数据进行大气校正，NDVI可以反映植物长势，将NDVI转换成植被覆盖度，就可以代表植被分布情况 [19] [20]，那么我们可以将植被覆盖度作为加权值对研究区内碳储量的空间分布和碳储量的多少进行矫正，使结果更加准确。该方法中各树种生物量方程和含碳率同生物量方程法。

大气校正过程如图2所示：

Figure 2. Atmospheric correction process

图2. 大气校正过程

植被覆盖度计算公式如下：

$FVC=\frac{NDVI-NDV{I}_{soil}}{NDV{I}_{veg}-NDV{I}_{soil}}$

式中：

FVC——植被覆盖度；

NDVI——植被归一化指数；

$NDV{I}_{soil}$ ——完全是裸地或无植被覆盖区域像元的NDVI值；

$NDV{I}_{veg}$ ——完全被植被覆盖像元的NDVI值。

乔木碳储量计算公式如下：

$C=B\times CF\times FVC$

式中：

C——乔木碳储量；

B——单株乔木生物量；

CF——树种含碳率；

FVC——植被覆盖度。

3. 结果与分析

使用生物量方程法计算出研究区碳储量的空间分布图如图3所示。

Figure 3. Spatial distribution of carbon storage by biomass equation method

图3. 生物量方程法碳储量空间分布图

由图3可知，研究区整体以针阔混交林为主，估算的12类树种中水胡黄和落叶松分布较广、面积较大，其面积分别为24776.24公顷和21780.10公顷；柳树面积最小，面积为90.91公顷。大部分树种分布比较分散，其他阔叶乔木分布较为集中，集中在研究区东南地区。12类树种中单位碳储量最高的为冷杉，为5.69 MgC到7.91 MgC之间，单位面积碳储量最低的为杨树，为1.97 MgC到3.14 MgC之间，具体经过生物量方程法估算的各树种对应单位面积碳储量如表2所示。

使用光能利用率法估算出研究区碳储量的空间分布图如图4所示，从图中可以看出碳储量高的地区大部分集中在研究区的东南位置，还有极少部分分布在研究区东北位置，与其他阔叶类树种所在方位一致。由于光能利用率法是与乔木冠层的长势息息相关的，故叶面积越大或乔木冠层长势越好的地区，用光能利用率法估算出的碳储量越高。使用光能利用率法估算的研究区碳储量大小呈现由北至南，由西向东逐渐递增的趋势，最小值为0.57 MgC，最大值为9.75 MgC。

Table 2. Carbon storage per unit area of tree species

表2. 各树种单位面积碳储量

Figure 4. Spatial distribution of carbon storage by light energy utilization method

图4. 光能利用率法碳储量空间分布图

使用NDVI反演法计算出研究区碳储量的空间分布图如图5所示。

Figure 5. Spatial distribution of carbon storage by NDVI inversion

图5. NDVI反演法碳储量空间分布

如图5所示，NDVI反演法与光能利用率法估算出的碳储量空间分布情况恰恰相反，NDVI反演法估算出碳储量高的地区主要集中在研究区北部，呈现从北至南逐渐递减、由东向西逐渐递增的趋势。其碳储量的最小值为1.36 MgC，最大值为8.68 MgC。从最大值和最小值上看，NDVI反演法与生物量方程法更加接近。

在Arcgis中采用生成随机点的方法，在三张结果图中各取相同位置的110个随机点，由于研究区无样地数据支撑，故以生物量方程法计算得出的碳储量作为真实数据，将110个样本点对应的碳储量进行精度分析，结果如图6所示。

从图中可以看出光能利用率法计算出的碳储量极值差异较大，数据波动较大，而用NDVI反演法计算出的碳储量只有个别值较真实数据差距较大，整体趋势与真实数据接近。光能利用率法的RMSE和R²分别为23.35和−1.72，拟合效果很差，NDVI反演法的RMSE和R²分别为8.031和0.6784，拟合效果很好，故NDVI反演法的精准度较高。

4. 结论与讨论

在碳中和的大背景下，发展林业碳汇是根本趋势，本研究根据林业碳汇量估算特征，将碳汇量的估算范围从整个生态系统，聚焦到乔木碳储量上，减少对土壤、枯落物等无法人工干预改善部分的估算，节约人力物力等成本；对于无样地数据支撑的地区，也可以依靠林地二类调查数据和遥感数据，进行不限时空的横纵向碳汇能力对比；使用遥感方法作为林业碳汇量估算的扩充，可以使最终估算结果以空间分布图的形式呈现，可以更加直观地看出碳储量地分布情况，后续想要提高和改善碳汇能力，对于需要改善的区域位置的选择具有指导意义，也大大方便了对后续碳汇项目的监测范围管理。

Figure 6. Accuracy analysis of carbon storage corresponding to sample points

图6. 样本点对应碳储量精度分析

本研究选用的两种遥感手段，一个是全部依靠遥感数据的光能利用率法，一个是结合地面信息，将植被覆盖度作为加权的NDVI反演法。最终结果无论从空间分布还是碳储量的拟合程度，都是NDVI反演法更精准。光能利用率法与NDVI反演法差异大的原因，主要是由于光能利用率法的碳储量估算是根据叶面吸收辐射值和反射辐射值的不同，来区分植被位置分布和植被光合作用的强弱，故光能利用率法估算出的碳储量多少与植被叶面积大小有着直接关系。由于研究区地处寒温带，以针阔混交林为主，阔叶乔木密集的地区，使用光能利用率法估算的碳储量就越多，但针叶乔木的含碳率是普遍大于阔叶乔木的，所以光能利用率法在针阔混交林或者其他树种繁杂的地区，估算精度是不准确的。

本研究由于没有样地数据支撑，只能以基于林地二类调查数据、通过生物量方程法估算的碳储量作为真实数据，进行回归拟合精度分析。所以，NDVI反演法的精确度也仅限于对于没有样地数据支撑，或是需要对建立样地之前年份的碳储量估算的地区。