Reconstruction and fusion correction method of multi-source LAI products in Maoershan forest based on ensemble Kalman filter

研究区位于黑龙江省尚志市帽儿山实验林场(45°14′~45°29′ N，127°29′~127°44′ E，图1)，平均海拔为428.0 m，属于温带季风气候，年均气温为2.8 ℃，土壤为典型暗棕壤。植被属于长白山植物区系，主要针叶树种有红松Pinus koraiensis、落叶松Larix gmelinii、樟子松Pinus sylvestris var. mongolica等，主要阔叶树种有白桦Betula platyphylla、山杨Populusdavidiana、色木槭Acer mono、水曲柳Fraxinus mandshurica、胡桃楸Juglans mandshurica等^[17]。

Figure 1.  Geographical location and plot distribution of the study area

在研究区内根据植被类型选取2块1 km×1 km大样地(图1)，每个大样地内包含5块100 m×100 m小样地。使用LAI-2200和TRAC冠层分析仪测量LAI^[18−19]。使用GPS记录位置坐标。

考虑到研究区数据的可获取性，本研究选取了3种应用广泛的LAI产品(表1)。利用HEG软件，对MOD15A2H和VNP15A2H数据进行格式转换和投影变换。利用ArcGIS 10.2对CGLS 300m V1数据进行投影变换。

选取2017年云层覆盖率低于10%的11幅Sentinel-2影像，利用Sen2cor工具进行大气校正得到L2A级数据，利用SNAP工具箱将数据重采样至10 m分辨率，并转换为与LAI产品一致的UTM WGS84投影。

采用2017年8月14日的Sentinel-2 MSIL2A影像，利用ENVI 5.3进行最大似然法分类^[20]。定义了农田、阔叶林、针叶林、草地和其他5类训练样本，通过实地调查数据和30 m地表覆盖分类产品GLC_FCS30相结合的方式目视判读各样本点的类别，采用混淆矩阵的方法进行精度检验，计算得到的总体精度为79.34%，Kappa系数为0.7517，整体分类质量高。将分类后的影像掩膜提取生成10 m分辨率的阔叶林和针叶林土地覆盖图(图1)。

利用2014—2016年的MOD15A2H和VNP15A2H的质量控制数据集进行像元质量标识，采用时间加权插值算法建立两者的LAI背景库，基于背景库对2017年的MODIS LAI和VIIRS LAI进行数据重建，计算公式参考文献[15]。最后使用三次样条帽盖算法(LACC)^[21]对初修正的2个数据集和PROBA-V LAI数据集进行平滑处理，去除异常值。

根据解混理论，假设每个粗分辨率像元的反射率是构成混合像元的每个土地覆盖类型反射率的线性组合，同理混合像元的LAI由不同类别的LAI及其在像元中所占的面积比例(即丰度)构成，基于Sentinel-2的10 m土地覆盖类型图来计算丰度，利用约束最小二乘法计算混合像元内不同类别的平均LAI值，并根据土地覆盖类型图分配到对应类别的像元上，从而得到10 m分辨率的高精度LAI数据^[22]。

利用降尺度后的3种LAI数据集分别构建经验性的LAI动态模型，将模型预测的LAI作为参数，连同其他参数一起输入PROSAIL模型模拟森林冠层反射率。同时，将Sentinel-2的红光波段(Band 4)和近红外波段(Band 8)影像作为观测数据，通过EnKF算法比较观测数据集合和模拟反射率集合，对预测LAI进行校正。结合精度评价结果对同化后的3种LAI数据进行赋权融合。

LAI动态模型能模拟叶面积指数^[23]，描述LAI随时间的动态变化过程，计算公式参考文献[24]。PROSAIL模型是由叶片光学特性模型(PROSPECT)和冠层二向反射率模型(SAIL)组合成的耦合模型，参数设置参考文献[25−27]。PROSPECT模型输入参数为：叶片结构参数(N=1.4)、叶绿素含量(Cab=58 μg·cm⁻²)、叶片含水量(Cw=0.025 g·cm⁻²)和干物质含量(Cm=0.009 g·cm⁻²)；SAIL模型输入参数为：叶面积指数(LAI)﹑土壤亮度参数(Psoil=0.25)、热点参数(Hspot=0.037)、太阳天顶角(tts=30)、观测天顶角(tto=0)和相对方位角(psi=0)，叶片反射率和透射率采用PROSPECT模型输出数据，tts、tto和psi可通过遥感影像头文件获取。

EnKF是基于蒙特卡罗的集合预报方法估计预报误差统计量的数据同化技术^[11]，计算公式如下：

式(1)~(4)中：${\boldsymbol{A}}=({x}_{1}, {x}_{2}, \cdots , {x}_{N})\in {{{{{\boldsymbol{R}}}}}}^{n\times N}$，是包含N个成员的状态变量集合，x表示n维状态向量(LAI)；M代表LAI动态模型，$ \varepsilon $为模型误差；${{\boldsymbol{A}}}_{t}^{\mathrm{f}}$表示t时刻模型预测的状态变量集合；a表示分析状态变量，T表示转置矩阵。${\boldsymbol{D}}\in {{{{{\boldsymbol{R}}}}}}^{m\times N}$代表扰动观测集合矩阵，是通过给定的m维观测向量${\boldsymbol{d}}\in {{{{{\boldsymbol{R}}}}}}^{m}$加上N个扰动构成，本研究中m=2。${\boldsymbol{P}}\in {{{{{\boldsymbol{R}}}}}}^{n\times n}$代表预测误差协方差矩阵，${\boldsymbol{R}}\in {{{{{\boldsymbol{R}}}}}}^{m\times m}$代表观测误差协方差矩阵，其算法与P相似，如式(3)所示。式(4)中，$\overline{{{\boldsymbol{A}}}_{t}^{\mathrm{a}}}$为状态最优估计值。H代表PROSAIL模型，因为该模型是非线性观测算子，不适用于标准EnKF，所以构建了扩展集$\widehat{{\boldsymbol{A}}}=[{{\boldsymbol{A}}}_{t}^{\mathrm{f}},\; h({{\boldsymbol{A}}}_{t}^{\mathrm{f}}\left)\right]$，将H转化成线性算子$\widehat{{\boldsymbol{H}}}$，此时P为扩展集协方差矩阵，故式(2)改写如下：

本研究观测集合和第1次的LAI预测集合是通过添加N个扰动(符合均值为0，方差为1的标准正态分布的随机数)生成，其中集合数N取值为100。

对3种同化LAI数据的精度检验结果归一化处理后分别赋予权重，将3种同化LAI数据分别与各自的权重相乘得到融合LAI数据。除了相关系数外的指标数值的绝对值越大，精度越差，需要在归一化处理前对它们取倒数，计算公式如下：

式(6)~(7)中：n表示评价指标的个数，m表示LAI同化数据的个数，i=1$, \cdots , $ m，$ {Y}_{i,j} $表示第i个LAI同化数据对应的第j个指标的值，$ {Y}_{i} $表示第i个LAI同化数据在所有指标进行归一化后的值。

利用实测数据对3种同化LAI数据和融合LAI数据进行精度评价，评价指标选用相关系数(R)、均方根误差(E_RMSE)、平均绝对误差(E_MAE)和相对偏差(E_RE)，计算公式参考文献[28]。

由图2可见：原始MODIS和VIIRS LAI时序曲线存在较为剧烈的跳跃现象，影响了其对植被物候特征的表征，而PROBA-V LAI时序曲线较为平滑，非生长季出现较小的跳跃现象。各产品都清晰表现了植被LAI夏季高、冬季低的生长特点。在植被LAI关键点缺失的情况下，MODIS和VIIRS产品基于背景库修正的方法，能够更好地重建LAI时序曲线，修正后的LAI与PROBA-V LAI差异较小，且时序曲线在形态上与PROBA-V也更加接近。修正后的LAI使用LACC算法进行平滑处理，平滑后的非生长季LAI均值都有所提高。在生长季的阔叶林VIIRS和PROBA-V LAI变化趋势相似；在生长季的针叶林MODIS和PROBA-V LAI变化趋势相似。在生长季，PROBA-V阔叶林的LAI高于针叶林。总体上，PROBA-V LAI均值高于MODIS和VIIRS。

Figure 2.  Changes and reconstruction effects of 3 LAI time series curves in broad-leaved forest and coniferous forest areas

由表2可见：3种LAI产品与实测数据的相关系数均较低，均方根误差、平均绝对误差和相对偏差值均较高，说明3种LAI产品与实测数据之间的一致性不高，没有明显的线性相关关系。经过EnKF同化后，3种LAI产品的相关系数均提升到0.49以上，MODIS LAI、VIIRS LAI和PROBA-V LAI的相关系数在阔叶林分别提高了0.57、0.52和0.57，在针叶林分别提高了0.52、0.32和0.40。3种同化后LAI数据的均方根误差均在0.40以下。其中，PROBA-V LAI在阔叶林和针叶林的均方根误差最小。VIIRS LAI在阔叶林的均方根误差小于MODIS LAI，而MODIS LAI在针叶林的均方根误差小于VIIRS LAI。相比同化前，VIIRS LAI在针叶林的均方根误差减小程度最低，仅减小了0.67，而VIIRS LAI在阔叶林的均方根误差减小程度最高，减小了1.69。对比3种同化后LAI数据的平均绝对误差，MODIS和VIIRS LAI在阔叶林的平均绝对误差最小，均为0.60，其中，VIIRS LAI在阔叶林的平均绝对误差相比同化前减小了1.11，减小程度最高。对比相对偏差，同化后除了MODIS LAI在阔叶林的相对偏差为负值，其他相对偏差都为正值。相比同化前，除了MODIS LAI在针叶林的偏差增大(同化前相对偏差=−0.72%，同化后相对偏差=8.34%)，其他数据的偏差均有不同程度的降低。

由图3可见：同化后LAI数据在阔叶林和针叶林从大到小依次为PROBA-V、VIIRS、MODIS，PROBA-V LAI具有明显的高估现象，与表2的相对偏差情况一致，与PROBA-V LAI原始数据就具有高估情况有关。同化后MODIS和VIIRS LAI数值相近，与两者原始数据就具有较高一致性有关。

Figure 3.  Scatter verification diagram of assimilated LAI and measured value

由图4可见：相比同化后LAI，融合后LAI与实测数据的相关性较高。融合LAI在阔叶林和针叶林的相关系数分别为0.83和0.76，均方根误差均在0.20以下，平均绝对误差和相对偏差趋势相同。在阔叶林，相比同化后的MODIS、VIIRS和PROBA-V LAI，融合后的LAI相关系数分别提高了0.24、0.27和0.21；对比同化后的MODIS和VIIRS LAI，融合后的LAI均方根误差分别降低了0.22和0.16；融合后的LAI平均绝对误差相比同化后的降低了0.02~0.28；融合后的LAI 相对偏差相比同化后的PROBA-V LAI，偏差降低程度最大，降低了8.59%。在针叶林，相比同化后的MODIS、VIIRS和PROBA-V LAI，融合后的LAI相关系数分别提高了0.17、0.27和0.20；均方根误差分别降低了0.11、0.15和0.06；对比同化后的PROBA-V LAI，融合后的LAI平均绝对误差降低了0.29；融合后的LAI相对偏差相比同化后的PROBA-V LAI，偏差降低程度最大，降低了9.88%。

Figure 4.  Scatter verification diagram of fused LAI and measured values

一般而言，LAI时序曲线是光滑且连续的，因此，当严重锯齿状的变化曲线中出现与正常趋势差距较大的拐点时，认为不符合植被正常生长规律^[29]。MODIS和VIIRS LAI时序曲线存在剧烈的跳跃现象，PROBA-V LAI即使利用时间平滑和间隙填充技术，使得时序曲线较为平滑，但仍存在较多缺失值^[30]，因此在产品应用前进行重建至关重要。本研究对3种LAI产品进行重建，初步改善了数据集质量较差的问题，使得三者LAI时序曲线变化趋势趋于一致。由于LAI背景库值反映的是研究区LAI生长变化的平均趋势，受空间和时间因素的影响非常大，因此与实际值有一定的偏差。本研究利用混合像元分解法进行降尺度处理，由于粗分辨率像元存在异质性，会存在一定的不确定性，经过降尺度处理，样地内同质像元具有相同的LAI值。本研究表明：经过同化方法校正了LAI值，使得LAI产品与实测值具有较高的相关性，但存在低值高估和高值低估的现象。经过EnKF同化，3种LAI产品精度和一致性均有所提高，经过同化融合后的LAI数据，降低了单一LAI数据误差，可有效弥补地面实测数据缺乏的不足，但能否满足植被监测和地表模型输入参数的要求还需进一步验证。本研究通过将LAI数据降尺度后作为同化算法的输入数据，使融合LAI的精度得以提高，生成10 m空间分辨率LAI数据，然而由于缺少研究区长期有效的观测数据，生成数据的时间分辨率受限制，在后续研究中可以考虑加入其他反射率数据，以期得到高精度长时序的LAI产品。

本研究结果表明：3种原始LAI产品与实测数据的一致性均较差，精度较差；EnKF同化后的LAI更接近实测数据；相比针叶林，阔叶林LAI同化效果较好；总体上MODIS LAI同化效果最好，其次是PROBA-V LAI，VIIRS LAI在阔叶林和针叶林的同化效果存在较大差异。融合后的LAI精度优于同化后单一的LAI产品，其相关系数为0.83和0.76，比同化后LAI提升了0.17~0.27。同化融合后的LAI数据提高了单一LAI产品的精度，减小了不同源数据的差异，达到了充分利用多源遥感数据优势和提高产品精度的目的，可以为开展高精度LAI数据生成方法的研究提供参考。