Forest above ground biomass estimation using airborne P band polarimetric SAR data

研究区为位于内蒙古呼伦贝尔盟的根河市大兴安岭森林生态系统国家野外科学观测研究站(大兴安岭生态站)，50°49′~50°51′N，121°30′~121°31′E。该研究站是中国纬度最高的森林生态研究站，面积为102 km²。研究区气候类型为典型的寒温带大陆性季风气候。研究区地势相对平缓，区域内80%的坡度小于15°，海拔高度分布在800~1 200 m，森林覆盖率大于75%，主要树种为兴安落叶松Larix gmelinii、白桦Betula platyphylla、樟子松Pinus sylvestris var. mongolica等。

采用的P波段SAR数据由机载CAMSAR系统获得^[15]。通过协议方式获取根河实验区机载P波段全极化(HH/HV/VV/VH) SLC数据。该数据是2013年9月以“奖状Ⅱ”飞机平台，以CASMSAR系统右视观测获取的全极化SAR数据。获取过程中飞行方向为自西向东，飞行高度为5 807 m，获取时间为2013年9月13—16日，极化方式为HH/HV/VH/VV，产品模式为SLC，幅宽为6 km × 7 km，中心入射角为55.058°，距离向分辨率为0.666 m，方位向分辨率为0.625 m。机载P波段SAR数据的预处理关键步骤包括多视处理、正射校正、入射角校正和极化方位角校正。多视视数在距离向和方位向均为3，数据的正射校正、入射校正参考文献[10]；极化方位角校正过程和算法见文献[16]。P波段SAR数据预处理结果、森林AGB抽样点及林分概况见图1。

Figure 1.  P band SAR data (A), the distribution of samples (B) and stand examples (C)

LiDAR获取的数字表面模型(DSM)、数字高程模型(DEM)用于P波段SAR数据的地理编码，由冠层高度模型(CHM)获取的LiDAR森林AGB数据用于反演模型的训练及验证。本研究获取的机载LiDAR数据是将Leica机载雷达系统荷载于“运-5”飞机平台上，于2012年8—9月在根河实验区开展飞行任务。该原始数据密度为5.6个·m⁻²的点云数据，激光中心波谱值为1550 nm，在初始点云数据的基础上，提取了研究区域高精度的DEM (图2A)、CHM (图2B)和森林AGB(图2C)等衍生产品。DSM、DEM、CHM数据的详细生成方法参考文献[17]，高精度LiDAR森林AGB的详细提取过程与方法参考文献[10,18]。

Figure 2.  Lidar derived data

为保证建模样本和验证样本能够代表整个研究区的森林AGB水平，以高精度LiDAR森林AGB图为基础，按照750 m的空间采样间隔，在ArcGIS中采用交互人工干预(去除道路及裸露地物)的方法选取113个样点(图1B)作为森林AGB反演模型的训练与验证。机载P波段SAR数据覆盖区森林林相不均匀，平均AGB较低，约46.7 t·hm⁻²，且大于100 t·hm⁻²的采样点仅有5个。113个样点AGB以10 t·hm⁻²间隔得分布情况(图1B)：0~10 t·hm⁻²，3个；＞10~20 t·hm⁻²，14个；＞20~30 t·hm⁻²，17个；＞30~40 t·hm⁻²，17个；＞40~50 t·hm⁻²，18个；＞50~60 t·hm⁻²，14个；＞60~70 t·hm⁻²，10个；＞70~80 t·hm⁻²，8个；＞80~90 t·hm⁻²，7个；＞90 t·hm⁻²，5个。

极化目标分解方法是从全极化数据中提取地物极化信息的有效方法，目前多种极化分解参数在森林类型识别中表现出巨大的潜力。本研究采用目前常用的极化后向散射系数、常用的3种极化分解方法提取极化SAR特征，并基于此分析森林的极化散射特征。提取P波段HH、HV和VV等3个极化的后向散射系数，基于3个后向散射系数的雷达植被指数(RVI)、极化辨别率参数(PDR)；基于Freeman-Durden三分量分解的体散射分量(FVOL)、单次散射分量(FODD)、二次散射分量(FDBL)、体地散射比分量(FD1/FD2，1表示地散射分量是ODD和DBL的和，2表示地散射分量仅为ODD)；基于Yamaguchi的体散射分量(YVOL)、单次散射分量(YODD)、二次散射分量(YDBL)、螺旋体散射分量(YHLX)；基于H-A-ALPHA极化分解的极化散射熵(entropy)、反熵(anisotropy)、散射角(alpha)、目标方位向角(beta)、相位差角1(gamma)、相位差角2(delta)。3种极化分解方法参考文献[19-20]。由于本次飞行试验时并未布设角反射器用于定标，因此本研究使用的后向散射系数值仅有相对含义。

森林在P波段极化特征响应分析的目的是为了确定P波段对森林AGB动态变化敏感的极化特征参数，从而确定有效的极化特征参数进行森林AGB的估测。将研究区的森林AGB划分为A (表示生物量在0~30 t·hm⁻²变化时对应后向散射系数的箱线变化，均值约20 t·hm⁻²)，B (31~50 t·hm⁻²，40 t·hm⁻²)，C (51~70 t·hm⁻²，60 t·hm⁻²)，D (71~90 t·hm⁻²，80 t·hm⁻²)，E (＞91 t·hm⁻²，100 t·hm⁻²)等5个变化等级，分别制作各等级相应P波段SAR提取参数值的箱线图，分析研究区各极化特征参数对森林AGB动态变化的响应，进而分析其极化散射特征。为了定量的分析各极化特征与森林AGB变化的关系，计算了它们之间的皮尔逊相关系数及显著性水平。

多元线性逐步回归模型(MLR)是森林AGB估测中最常用、最经典的方法之一。与常规的线性回归模型相比，MLR可同时完成模型输入参数的优选，进而提高森林AGB估测的效率和精度。MLR是将自变量逐个引入，每次判断自变量对因变量影响的显著性，并对模型中的自变量进行检验，逐个从模型中剔除不显著的变量，从而得到最优模型。即在保证显著性值在0.05以下的情况下，筛选相关性高的特征变量，进而得到因变量的最优估计。MLR的实现算法详见文献[21]。

KNN^[22-25]、SVR^{[22, 24]}和RF^{[24, 26]}是森林AGB估测中常用的非参数模型，与参数模型相比，无固定的模型结构，通常通过数据驱动的方法来确定模型结构并用于森林AGB的估测，因此也称为机器学习方法。在森林AGB估测中，这3种方法各有优势，因此选取这3种方法来探究非参数方法在P波段森林AGB估测中的潜力。

反演结果精度的定量评价通过反演结果与真值之间的决定系数(R²)、均方根误差(RMSE)、平均绝对误差(MAE)、估测精度(Acc)来表征。这4个参数的计算公式参考文献[27]。

本研究获取的P波段中心波长为64 cm，该波段在森林中的穿透性较强。由图3可知：P波段20个SAR特征均表现出对研究区森林AGB变化的敏感性。P波段3种极化方式后向散射系数对森林AGB在0~110 t·hm⁻²内动态变化的响应与已有研究^[28]相似，即后向散射系数的值均随着森林AGB的升高而缓慢增长。在各森林AGB水平中，HH、HV和VV散射能量强度相差并不明显，并且从箱线图(图3)的中值可以看出：HH的单调增长趋势较其他2个极化明显，且异常值较少。这可能是由于P波段波长较长，比较粗树干等圆柱体形状的散射体为森林中的散射体，使得去极化特征明显降低，因此同极化的HH和VV后向散射能量对AGB变化的敏感性明显高于HV^[29]。此外，由后向散射系数计算的RVI也随着森林AGB增加而单调增加且未出现饱和现象，并且通过箱线图的宽度变化可以看出：各水平RVI箱宽变化不明显，说明该值对森林结构的变化不敏感，而对森林AGB的变化比较敏感。相比RVI，PDR非单调增长，在不同森林AGB水平，占主导散射的散射体会发生明显变化^[30-31]，该现象也可由Freeman-Durden中各个参数随森林AGB变化的响应加以说明。P波段Freeman-Durden分解相关特征中，除FODD和FDBL外，均随着森林AGB的增加而增加，并且2种体地散射比的值均远远大于1。Yamaguchi分解特征对森林AGB的响应趋势与Freeman-Durden分解相关特征的趋势基本一致。P波段H-A-ALPHA分解提取的参数中，散射熵(entropy)随着森林AGB增加而单调增加，而alpha则敏感性不明显。但beta、gamma和delta对森林AGB变化的敏感性则高于散射角alpha。P波段森林的极化散射特征随森林AGB的变化与已有研究中短波长的研究差异明显，说明现有的极化分解建模方法与该森林AGB水平P波段的散射特征有明显差异。

Figure 3.  Forest scattering mechanism analysis at P band

由表1可知：极化分解中表征体散射分量的极化特征与森林AGB变化的显著性较低或不存在显著性，而二次散射特征分量则与森林AGB变化具有强相关性且显著性水平最高。这说明在本研究区，树干与地表形成的二次散射对森林AGB的变化最为敏感，且目前极化分解中的体散模型并不适合P波段研究区森林AGB水平下森林的散射机制。表1的定量分析结果与图4的定性分析结果对比可知：由于各极化参数的值域范围差异明显，图3对森林AGB敏感性较明显的参数在定量分析时相关性并不一定最高，因此图4的分析仅可作为初步参考，对于森林AGB敏感的极化特征参数的分析，仍需要采用皮尔逊系数进行相关性定量分析。

Figure 4.  Scatter plots of estimated and LiDAR forest AGB of four models

从表2可知：4种模型中，MLR模型估测结果精度最低，RF精度最高；而KNN和SVR模型的估测精度则相近，R²相同，RMSE、MAE和Acc也基本相同，比RF估测结果的最高值仅低约2%。

从图4可看出：在整个森林AGB分布中，MLR估测结果分布较其他3种非参数方法分散，在森林AGB大于80 t·hm⁻²时出现了明显的饱和低估现象。其他3种方法尽管也有低估的现象，但是饱和现象并不明显。此外，本研究中4种方法估测结果的相对误差约30%，而以往区域性森林AGB的估测误差为37%~67%^[32]，采用P波段HV后向散射系数的估测结果中，同质性森林地区的相对误差约13%，而异质性地区相对误差则约60%^[28]。在采用L-波段极化分解参数进行森林AGB反演研究中，在森林AGB水平低于120 t·hm⁻²时，后向散射的特征优于3种极化分解的特征，这与本研究中P波段的研究结果基本一致^[33]。4种方法的估测结果中均出现了低值高估和高值低估的现象，这与以往采用不同遥感数据源进行森林AGB估测的结果一致^[32]。图5中，RF的残差分布最接近高斯分布，尽管峰值出现在10 t·hm⁻²左右，由于其值分布较窄，因此具有较高的估测精度。KNN和SVR的残差分布图均较为连续，并且分布区间明显高于MLR，但由于正值占比较大，因此总体上呈现高估现象。MLR估测结果的残差分布图在10和20 t·hm⁻²出现了2个明显的峰值，且所有残差值的分布范围较宽，解释了其估测精度低于其他3种方法的原因。在ENGHART等^[34]基于C和L-波段参数和非参数方法森林AGB估测的研究中，也发现MLR的估测结果要略低于非参数模型。然而参数和非参数模型的适用性还受到训练样本大小的影响，因此在训练样本较大时可选择非参数模型进行估测，当样本较小时则可优先选择MLR方法进行森林AGB估测^[32]。

Figure 5.  Histograms of difference between forest AGB estimation results and LiDAR sampling points

为了进一步分析P波段极化特征对森林AGB估测的潜力，本研究将研究区森林AGB划分为不同的等级，然后采用4种估测模型中估测结果最优的RF模型进行估测，估测结果见表3。在本研究中将森林AGB划分为3组，其中第1组分别以30、60 t·hm⁻²为边界划分为3个子组；而第2组和第3组分别以40和50 t·hm⁻²为界划分为2个子组。由表3可知：分组界限不同，估测精度有明显的差异，分组划分越详细，估测精度(Acc)越高；此外，3种分组情况的估测结果均表明：在森林AGB平均值约45 t·hm⁻²，最高值不超过120 t·hm⁻²时，P波段在森林AGB水平较高的分组估测精度较高，如在第1组中，森林AGB大于30 t·hm⁻²分组的估测精度比小于30 t·hm⁻²分组的估测精约高5%；而在以50 t·hm⁻²为分组界限的2组中，森林高AGB组的估测精度比低AGB组的估测精度约高出6%。表3的结果分析表明：待估森林的AGB水平对P波段极化特征进行森林AGB估测的估测精度有明显影响，且P波段极化信息更适合森林AGB较高区域森林AGB的估测。

针对P波段极化SAR数据在森林AGB估测中的潜力，提取了20个P波段极化SAR参数，探索了森林AGB估测中常用的MLR、KNN、SVR和RF方法在使用P波段极化SAR数据进行森林AGB估测的潜力。结果表明：P波段极化SAR信息在森林AGB估测中具潜力，但估测精度受到待估区域森林AGB水平高低的影响；4种估测方法中，非参数方法的估测结果明显优于MLR估测。P波段森林AGB估测结果中，同样存在AGB低值高估和高值低估的现象，其原因仍需要进一步探索。此外由于本研究区的森林AGB均值为45 t·hm⁻²，最高值低于120 t·hm⁻²，所以P波段极化信息在AGB高于120 t·hm⁻²的森林覆盖区中对森林AGB的估测能力仍有待进一步研究。