基于INFORM模型与GEE的根河火烧迹地冠层含水量反演研究

饶月明，刘凯军，黄华国，林起楠，王 川

（1.北京林业大学 省部共建森林培育与保护教育部重点实验室,北京 100083；2.内蒙古自治区巴彦淖尔市乌拉特后旗林草局，内蒙古 乌拉特后旗 015543）

水分含量是评价森林健康的关键参数，代表了森林对气候、干旱和虫害的综合反馈[1]。植物冠层中水分含量占40%～80%，并且是影响绿色植物光合作用和生物量的主要因素[2]，叶片含水量被认为是监测植物干旱胁迫的有效指标[3]，土壤水分在虫害初期由于疏导组织受阻无法顺利到达冠层[4]，森林火灾蔓延过程中产生的飞火与树冠火均与冠层可燃物含水量有密切联系[5]。森林火灾破坏了森林的结构和功能，降低了涵养水源的能力。冠层水分含量（Canopy water content，CWC）是叶片的等效水厚度（Equivalent water thickness，EWT）与叶面积指数（Leaf area index，LAI）的乘积，表示单位面积冠层水分总厚度。植物水分主要吸收近红外与短波红外的太阳辐射，植物的反射率在970、1 200、1 450、1 950 和2 500 nm 光谱范围内有不同程度的降低，因此利用光谱信息可定量反演植物的水分含量[6-9]。

遥感提取植物的水分含量主要包括统计模型和物理模型两种方法。统计模型通常使用光谱指数[10]、人工神经网络[11]和最小二乘回归将叶片水分参数与反射率关联提取植物水分含量。然而，统计模型的建立需要标定，对于大区域反演能力有限。相比而言，物理模型研究光线与植被生化组分之间相互作用的辐射传输过程，从机理上解释植被光谱与生化组分之间的关系，物理模型更有通用性[12-14]。Zhang 等[15]指出，高光谱数据有助于提高辐射传输模型反演植被水分含量的精度。Trombetti 等[11]使用人工神经网络与辐射传输模型从MODIS 中提取了植物的水分含量，并以高光谱数据验证。Jurdao 等[16]通过PROSPECT 模型和GeoSAIL 模型反演了西班牙部分林地的含水率，并证明林木含水率的反演有助于森林火险的评估。Quan 等[17]耦合不同的物理模型提取森林可燃物的含水率，达到森林火险预警的效果。但使用物理模型反演森林植被含水量的研究相对较少。原因可能是：1）物理模型在反演植被水分含量时面临病态反演的问题[18]；2）相对于多光谱数据，高光谱数据偏少，难以用于大面积反演研究[15]；3）植被含水量在物理模型中属于弱敏感参数，反演过程中容易受到其它因素干扰，具有一定反演难度[19]。如果能解决上述问题，就可以推广物理模型，获得更好的反演精度。选择合适的物理模型对于反演精度提升至关重要，混合模型（如GO-RT 模型[20]、FLIM 模型[21]、INFORM 模型[22]）结合几何光学模型在解释阴影投射面积和地物表面空间相关性上的基本优势，与辐射传输模型在解释均匀媒质中多次散射的优势，对于火烧迹地此类稀疏分布且具有一定冠层高度的林地，模拟效果更好。

遥感大数据的发展也进一步降低了病态反演的问题。Google Earth Engine（GEE）平台是Google 公司、卡耐基梅隆大学、美国地质调查局联合开发推出的遥感大数据云端运算平台，其提供了海量遥感数据，对高效检测森林植被变化有重要意义[23]。

内蒙古自治区根河市是我国重要的天然林保护区，但在春夏常发生森林火灾，特别是1987年5月和2003年5月发生了特大森林火灾，对森林生态环境造成了严重破坏，森林涵养水源的能力降低，面临干旱、虫害与再次火烧的风险。本研究使用INFORM 模型反演火烧迹地的冠层含水量，并基于GEE 对火烧迹地CWC 恢复进行时序性分析，最后绘制了2018年根河地区CWC 的分布图，为量化监测火后森林植被水分恢复提供技术支持。

1 研究区概况与研究方法

1.1 研究区概况

研究区位于内蒙古自治区根河大兴安岭林区，地处120°12′～122°55′E，50°20′～52°30′N，海拔多集中于700～1 300 m。该地区属于寒温带，具有大陆性季风气候，年降水量为450～500 mm，雨季集中在7—8月。根河地区地形由西南至东北逐渐升高，东北区域以坡度小于15°的山地为主。森林资源丰富，森林覆盖率为75%，树种类型以兴安落叶松Larix gmelinii与白桦Betula platyphylla为主。该区域常年发生森林火灾，造成了巨大的经济损失和生态损失，当地一般在春夏季节封山防火。图1是研究区样地分布。

1.2 地面实测数据

样地分布于根河大兴安岭林区，火灾发生年份分别是1987年、1998年、2003年、2010年和2015年，同时选取顶极群落作为对照样地，样地分布见图1。每个样地群选取2～5 块样地，结合森林火灾资料、历史影像数据和野外调查数据，分类别确定为中度火烧、重度火烧与对照样地。使用差分GPS 确定样地坐标，样地大小为30 m×30 m，由于火烧迹地分布广泛且地形复杂，实际调查样地总数有20 个。

图1 样地分布Fig.1 The location of plots

样地中针叶树种是兴安落叶松，阔叶树种是白桦。使用搭载鱼眼镜头的相机和CAN-EYE 软件提取叶面积指数作为样地的估计值[24]，提取过程见图2。通过每木检尺确定不同树种的标准木，并将标准木树冠分为上、中、下3 层，用高枝剪采集3 层叶片样本，每层取5～10 片（针叶每层取10～20 束）。扫描并提取阔叶表面积，针叶通过平铺测量边长确定表面积，并将所采集叶片带回实验室，使用烘干法确定含水量。

图2 基于鱼眼镜头的冠层LAI 提取Fig.2 The extraction of canopy LAI using spherical lens

叶片等效水厚度（EWT，g/cm2）指叶片的含水量与叶面积的比值，表示单位叶面积的含水量。

式（1）中：mf指叶片鲜质量，g；md指叶片干质量，g；A指叶片表面积，cm2。

冠层含水量（CWC，kg/m2）通过计算各树种叶片EWT、样地LAI 与各树种比例乘积得到[25]。

式（2）中：f为各树种在样地树种中的比例。

20 个样地的EWT、LAI、株树密度（SD）、叶绿素含量（Cab）的统计值如表1所示，表2为林分样地信息。由表1可知，针叶树种EWT为0.016 2～0.028 2 g/cm2，平均值为0.021 7 g/cm2；阔叶树种为0.006 7～0.011 3 g/cm2，均值为0.009 6 g/cm2；样地总EWT 为0.010 1～0.024 2 g/cm2，均值为0.017 7 g/cm2；LAI 为0.01～2.83，均值为1.40。

表1 样地调查数据统计Table 1 Summary of statistics of leaf and canopy variables

表2 林分样地信息Table 2 Information of plots

1.3 遥感数据与预处理

使用遥感数据分为地面调查同期卫星数据和历史数据两类。同期遥感数据采用Landsat 8 OLI影像Level 1 产品，此产品经过了几何校正和地形校正。使用2018年7—8月期间无降水影像，原始数据经过了辐射校正、FLAASH (Fast line-ofsight atmospheric analysis of spectral hypercubes) 大气校正[26]等预处理流程。历史数据使用Landsat 5、Landsat 7 和Landsat 8 地表反射率数据集，在GEE平台上通过Javascript 编译代码对影像批量预处理，并提取各样地1985—2018年6—8月无云影响高质量地表反射率数据。

根据Landsat 8 OLI 数据波段特征，对已有的水分指数进行分析，最终选用归一化水分指数（Normalized difference water index,NDWI）[27]验证冠层含水量的反演精度。

式（3）中：NIR 代表Landsat 8 第5 波段反射率；SWIR 代表Landsat 8 第6 波段反射率。

1.4 INFORM 模型

INFORM 模型是一种结合几何光学模型和辐射传输模型的混合模型，模拟400～2 500 nm 的森林二向反射率。该模型实际上结合了FLIM 模型、SAILH 模型和PROSPECT 模型[22,28]。对于特定波段，INFORM 模型反射率的公式为：

式（4）中：R是样地反射率；RC是无限厚度冠层的反射率；RG是背景反射率。C因子代表冠层因子，G因子代表地表因子。

C因子和G因子的计算公式为：

式（5）～（6）中：C因子主要与观测方向的冠层地表覆盖率Co与太阳方向的地面覆盖率Cs相关；而G因子主要与地面部分分量Fx相关。

Co和Cs具有一定的相关性，是由模型中冠层直径、光照和观测角度计算得到。对于给定的天顶角θo、Co的计算公式为：

式（7）中：SD 是林分密度（n·hm-2）；k值代表树冠的平均水平面积（hm2）。k的计算公式是：

式（8）中：CD 为冠层直径（m）。

INFORM 模型中的冠层透过率采用SAILH 模型计算得到，To与Ts的计算引入平均叶面积指数LAI、平均叶倾角ALA (Average leaf inclination angle)等参数。

在INFORM 模型中，无限冠层深度的冠层反射率（RC）和背景反射率（RG）采用SAILH 和适宜本研究区的叶片光学模型计算，背景反射率RG计算中涉及的土壤反射率光谱使用研究区土壤光谱曲线。植物叶片等效水厚度（EWT）、叶绿素含量（Cab）、干物质含量（Cm）等叶片生化参数使用PROSPECT 模型中的输入参数。国外对于耦合INFORM 与PROSPECT 模型反演针阔混交林EWT 已经有过相关研究[25]。因此，使用INFORM模型计算样地反射率（R），可以将其表示为内部冠层参数、叶片参数和外部参数的函数，即：

式（9）中：N为描述叶片内部细胞结构的参量；Cw为叶片水分含量；Cm是叶片干物质含量；Scale代表土壤反射率因子；lais为单株树LAI 值；laiu为地面冠层LAI 值；sd 为株树密度；h为树高；cd 为冠层直径；ala 为平均叶倾角；tetas代表太阳天顶角；tetao为观测天顶角；phi 为方位角；skyl为天空漫辐射参量；rsoil为土壤反射率。

需要特别指出的是，INFORM 模型中的LAI 分为单株树木lais(Single tree leaf area index)和地表灌木laiu（Leaf area index of understorey），本研究中林分的LAI 通过单株树木lais与覆盖率Co得到，即：

INFORM 模型中的EWT 等效水厚度为PROSPECT 叶片辐射传输模型的输入参数，样地对应的CWC 是EWT 与LAI 的乘积。

1.5 反演方法与策略

使用查找表（LUT）联合反演EWT 与LAI 得到CWC。LUT 实际上是建立不同的输入参数组合，得到不同的光谱数据组合，建立冠层含水量和光谱之间的关系，反演时用线性或非线性内插的方法计算出参数值。根据INFORM 模型EWT 与LAI在400～2 500 nm 之间参数敏感性结果与实际调查值确定查找表输入参数与变动范围，由此得到冠层含水量与冠层反射率的查找表。通过筛选查找表中的实测光谱值与模型模拟光谱值的RMSE值优化各参数值，筛选当RMSE 最小时为所选取参数的最优标定值。大量研究表明，INFORM 模型的EWT 仅在近红外波段与短波红外波段范围内有较高的敏感性，对应Landsat 8 波段为第5、第6 波段，而LAI 对应敏感的波段为2～6 波段[25,29-31]，LAI 敏感范围大于EWT 敏感范围。故本研究在筛选时，先使用LAI 敏感的波段筛选出RMSE 小的前10%参数值作为反演EWT 的查找表子集，再结合EWT 敏感的波段从子集中筛选出RMSE 最小的前5%参数值的平均值，得到LAI与EWT，相乘得到最终的CWC 值。

式中：Rmes为反射率实际值；Rmod为模型反射率模拟值。

1.6 反演结果的验证

使用实测数据与水分指数评价反演结果。其中样地反演结果结合实测LAI、EWT 和CWC，通过决定系数R2、均方根误差RMSE、标准化均方根误差nRMSE(normalized RMSE)评价反演精度。

式中：yi和yi′ 分别指实测和反演值；指样地测量的均值。

由于野外调查采集样本量偏少，在林区选择280 个样点，用以评价NDWI 与反演CWC 的相关关系，从另一个角度评价反演精度。其中70 个样点位于火烧迹地附近，70 个样点位于纯健康林地附近，其他为普通林地样点。

1.7 利用GEE 大数据平台时序性反演分析

使用GEE 在线数据，加载Landsat 8 OLI 地表反射率数据。编写代码预处理，筛选1987—2018年近30 a 的Landsat 系列数据，并选取6—8月同一时间段数据，预处理后，直接提取各样地火灾发生前3 a 至2017年的各波段反射率数据。根据不同的传感器，选择特定的光谱响应函数，得到适应不同传感器的查找表，以此反演得到不同时间段的火烧迹地LAI 与CWC，使用Mann-kendall时间序列非参数估计模型，确定LAI 与CWC 恢复速率。并基于2018年8月遥感数据，使用GEE重采样到1 km，使用同样反演方法得到整个区域的冠层含水量数据，并绘制分辨率为1 km 的森林冠层含水量的分布，图3为研究的技术路线。

图3 根河地区火烧迹地冠层含水量反演技术路线Fig.3 The flow chart about the canopy water content inversion of the fire scars in Genhe county

2 结果与分析

2.1 查找表的设立

结合INFORM 模型的参数敏感性分析[25]与样地实测数据，建立用于反演冠层含水量的查找表（表3）。INFORM 模型中所输入的参数包括叶片叶绿素含量Cab（μg·cm-2）、叶片等效水厚度Cw(g·cm-2)、叶片干物质含量Cm(g·cm-2)、叶片内部结构参数N等。N是描述叶片内部细胞结构的参数，与植物的种类和生长状态有关，设定N值为1.5。

表3 查找表Table 3 The information about the LUT

2.2 反演结果

联合反演结果显示，LAI 反演的R2为0.81，nRMSE为0.40，当LAI 高于1.5 时，模型模拟得到的LAI 值低于实际测量的LAI（图4a）。而叶片水分含量EWT 反演的R2为0.35，nRMSE为0.20，总体上模型模拟得到的EWT 明显高于实际测量的EWT（图4b）。CWC 反演的R2为0.79，nRMSE为0.52，植被水分在冠层尺度上的反演精度高于叶片尺度（图4c）。而实测的CWC 与NDWI 呈现指数关系，R2=0.78，精度较INFORM模型偏低（图4d），显示INFORM 模型反演冠层含水量的优势。造成反演单株树叶片含水量精度不如冠层含水量精度高的主要原因有：1）实际调查中混交林的叶片含水量难以精确确定；2）多光谱数据的局限性无法选择水分最敏感的光谱波段反演；3）冠层含水量CWC 由EWT 与LAI 共同组成，而LAI 的准确反演确实提高了冠层尺度植被水分反演精度。

2.3 NDWI 和CWC 对比

在研究区附近选取280 个象元，用同样方法反演CWC，并与NDWI 拟合。大范围反演结果（图4e）显示，NDWI 较小（如火烧迹地）时，CWC与NDWI 呈线性关系；当NDWI 区域饱和（健康林地）时，INFORM 模型仍可以继续反演CWC。结果显示，CWC 与NDWI 总体上呈明显的指数关系，R2为0.77。

图4 反演结果Fig.4 The result of inversion

NDWI 和CWC 的良好相关关系说明了局部地区反演的有效性。但在NDWI 饱和时，反演CWC仍在上升，反演CWC 与NDWI 总体呈现指数关系。一方面，CWC 与NDWI 相关性高，说明反演CWC 精度较好；另一方面，NDWI 达到饱和值时，而反演得到的CWC 未饱和，反映了物理模型反演可以提供更多量化的水分信息。对于统计模型，含水量反演通常需要精准的实测值标定，相比而言，INFORM 模型不需要标定，更适于大区域植被水分的反演。

2.4 基于Google Earth Engine 火烧迹地冠层含水量时序反演结果

基于GEE 大数据平台编写Java Script API 批量处理，通过INFORM 模型反演研究区1987—2018年生长季期间的LAI 与CWC。火灾影响后，其LAI 和CWC 都有明显下降，符合森林火灾的一般表现，结果如图5所示。火灾前3 a 到火灾后10 a的LAI 与CWC 的动态变化表明，部分火烧迹地在火灾发生前，其CWC 已呈下降趋势，可能是火灾发生的重要原因之一，如1987年样地、1998年样地；除火灾后经人为补植的2010年样地外，火灾发生后1～2 a 内，LAI 与CWC 均出现明显下降；火灾后10 a 的恢复过程中，部分样地（1987年样地，2003年样地）的LAI 与CWC 在火灾发生后一年内出现明显的恢复，2010年与2015年样地的CWC 与LAI 在火灾后均还未出现明显的恢复现象。

图5 火烧迹地LAI 与CWC 时序性变化分析（红线代表火灾发生年份）Fig.5 Analysis on temporal change of LAI and CWC in burned areas（The red line represents the burning year）

受不同火烧程度、气候和地理环境因素的影响，火烧迹地的恢复速率各不相同。基于Mannkendall 模型计算得到1987年、1998年、2003年3块火烧迹地10 a内LAI的年恢复速率分别为0.17、0.19、0.12，CWC年恢复速率分别为0.010 3、0.011 7、0.011 5 kg/m2，从而量化了LAI 与CWC的恢复速率。

续图5Continuation of Fig.5

其中1998年样地、2003年样地、2010年样地由于地形坡度较大，恢复速率并不理想，即使CWC 恢复至最高值，仍低于火前水平。其中2010年样地实地调查发现即使在火后人工补植，恢复情况仍不明显，且出现了病害影响。一方面由于森林火灾严重破坏了森林涵养水源的能力，另一方面由于补植树种没有符合“适地适树”原则，共同导致恢复结果不理想。

2.5 基于INFORM 模型根河地区冠层含水量反演结果

基于GEE 和INFORM 模型对2018年8月根河地区冠层含水量反演，得到分辨率为1 km 的冠层含水量的分布（图6），与实际调查结果吻合。火灾发生10 a 内火烧迹地冠层含水量大部分集中在第一、二等级（0.02～0.27 kg/m2），与其他未受到林火干扰的森林相比，冠层含水量明显降低，如图6所示红色窗口内为2003年火烧迹地，其冠层含水量存在明显的低值，而NDWI 表现并不明显。本森林冠层含水量反演产品可以为森林防火与森林病虫害防治提供监测预警信息。

图6 CWC 反演与NDWI 对比结果Fig.6 The comparison between CWC inversion result and NDWI

作为物理模型，INFORM 模型反演冠层含水量具有通用性，对于根河这类面积大且防火压力严峻的区域，植被指数通常会达到饱和值，无法测得更准确的冠层水分含量。且对于统计模型，不同区域冠层含水量与光谱指数往往呈现不同的关系，需要多次准确含水量实测值的标定，增加了大面积反演的难度。

3 结论与讨论

森林冠层含水量是评价森林健康的重要指标，本研究使用INFORM 模型建立了CWC 与森林冠层反射率的查找表，并结合Landsat 系列数据与Google Earth Engine 遥感大数据平台，实现了根河地区火烧迹地CWC 的定量反演，并进行时序性研究分析，为动态监测火烧迹地恢复速率和森林健康提供了新的方法。主要得到以下结论：

1）INFORM 模型中水分敏感的波段集中在近红外与短波红外，使用近红外与短波红外反演冠层含水量以评价植被水分，可以得到精确结果；2）基于INFORM 模型反演火烧迹地LAI 与CWC精度较高，适用于火烧迹地恢复监测；3）CWC反演精度比EWT 反演精度高，森林冠层平均叶片尺度的验证存在客观困难，难以精确估计林分平均叶片含水量，而引入森林结构参数LAI 可以提高冠层尺度上的水分反演精度。CWC 作为LAI 与EWT 的乘积，更能反映火烧迹地的恢复状况，更适用于本研究的评价；4）Google Earth Engine 作为前沿的遥感大数据平台，在解决森林生态遥感持续性监测相关研究时，具有方便快捷的特点；5）火烧迹地的冠层含水量在发生火灾后，与该地区其他森林相比，有明显下降。基于Mann-kendall模型计算得到火烧迹地LAI 与CWC 恢复的速率各不相同。虽然在林火发生后，部分地区补植过林木，但由于火烧迹地涵养水源的能力降低，出现过大面积的林木病害。由于火烧迹地中林木水分含量难以短时间恢复，火烧迹地有再次发生火灾与病害的风险；6）使用GEE 与INFORM 模型反演得到根河地区森林冠层含水量分布，冠层含水量可以作为检测森林健康的指标，对防止森林病虫害与森林火灾有一定意义。

植被含水量可以用于准确监测植被环境胁迫程度[32]、健康状况并进行火险评估[5,33]。本研究通过长时间监测与实地调查表明，冠层含水量低的区域林木生长状态普遍不佳（1998年样地、2003年样地、2010年样地），更易受病虫害影响（2010年样地），为森林健康监测提供了有效且便捷的手段。

国内外有关植被水分的遥感监测大多集中于农业，森林植被水分监测验证存在客观困难性。王长青[34]结合Sentinel-1B 和Landsat8 OLI 数据，通过经验模型反演小区域针叶林冠层叶片含水量，规则函数模型模拟最高精度为R2=0.629 9，与本研究结果（R2=0.79）相比，不但明显提高了反演精度，并且更适用于大范围反演。全兴文[19]使用物理模型反演草地冠层含水量R2=0.68，但通过考虑自由参数相关后，R2提高至0.84，但草地比森林更易于反演。Zhu 等[25]在反演混交林冠层含水量时，使用机载激光雷达数据作为结构参数数据，能够提高水分反演的精确度（R2=0.87），其部分结果与本研究相符。对于大范围反演研究，激光雷达和高光谱数据都较难获得，相比之下，本研究使用易获得的高光谱数据，推广性更好。

然而本研究还存在一些局限和不足。模型反演过程中未考虑参数之间的内在关系，存在病态反演问题，未来可以尝试通过考虑自由参数相关性来缓解病态反演问题。反演过程需要经过查找表多次迭代分析，计算量巨大对硬件有较高的要求。研究中的实测数据偏少，植被叶片水分的测量估计（特别是针叶）存在客观的瓶颈。GEE 作为大数据平台，提供了研究者持续性研究的快捷途径，但高光谱数据与激光雷达数据太少，难以提供更多准确的水分光谱信息与林分结构信息，是本研究的局限性。“高分五号”作为我国重要的高光谱卫星，能够为反演提供更多的光谱信息。这些问题的解决，有助于利用物理模型反演植被生化指数的进一步研究。