森林是陆地上最大的碳储存库，储存了陆地生态系统有机碳地上部分的76%~98%(刘华等，2005)，准确估测森林地上生物量(above-ground biomass，AGB)及其动态变化不仅是当前陆地碳循环研究的核心内容之一，同时也是森林生态学与全球气候变化研究的重要课题(Houghton et al., 2009)。遥感技术能够快速、准确、实时、宏观地获取森林植被信息，近年来，森林AGB的遥感估测研究受到国内外同行的广泛关注。光学遥感数据提供的光谱特征、植被指数特征以及纹理特征等能够较好地反映森林冠层信息，但对冠层下部枝干生物量的指示作用较弱，在森林AGB估测方面存在一定局限性(汤旭光等，2012；刘茜等，2015)。激光雷达(LiDAR)能获取高精度的森林垂直结构信息，在森林AGB估测方面起到了重要作用(曹林等，2013；穆喜云等，2015)，但LiDAR飞行成本较高，限制了其在大尺度范围内的应用(庞勇等，2005)。SAR具有较强的穿透性能，不仅可与树冠表层发生相互作用，还可与冠层下部枝干发生相互作用，对森林AGB具有较强的指示作用，以其全天时、全天候的工作能力和相对较低的数据获取成本，在区域及全球森林AGB估测中有着传统光学遥感和LiDAR无法比拟的优势(陈尔学，1999)。

目前，基于SAR数据的森林AGB估测方法主要有后向散射强度估测法、干涉相干性估测法、极化干涉SAR异速生长方程估测法和层析技术估测法等。相较于后向散射强度及干涉相干性估测法对生物量敏感性在高生物量值区间将会降低(Imhoff，1995；Sandberg et al., 2009)和异速生长方程估测法受树高估测精度以及模型误差的影响(Cloude et al., 2003；李春萍等，2007)，层析技术估测法能有效表征森林垂直方向空间结构信息，基于层析技术的森林AGB估测法已成为该领域前沿研究课题之一(李文梅等，2014)。层析技术作为一种三维成像技术，可实现森林内部散射体在垂直方向的有效分离，有助于深刻理解SAR信号与森林内部散射体相互作用的物理机制，在森林垂直结构参数反演领域具有巨大的应用潜力(Treuhaft，2002；Cloude，2006；庞蕾等，2010；庞礴等，2013)。

Reigber等(2000)首次利用多基线机载数据对森林试验区进行层析成像，成功获取了森林的垂直结构信息。Minh等(2014；2015；2016)利用P-波段多基线数据对热带雨林进行层析成像，发现30 m高度处层析相对反射率与森林AGB有较高的相关性，进而利用该高度处层析相对反射率进行森林AGB估测。由于多基线SAR数据获取难度较大，目前利用多基线InSAR层析技术提取森林AGB的研究在国内尚未见报道。国外虽已有多基线InSAR层析技术用于森林AGB估测的初步探索，但层析相对反射率三维分布提供较为精细的森林垂直结构信息，其与森林AGB相关的层析特征尚未得到充分挖掘和有效利用。因此，本文基于多基线InSAR层析技术，提取与森林AGB相关的特征参量，进一步构建森林AGB估测模型，从而提高其估测精度，为区域和全球森林生物量估测和碳储量研究提供关键技术支撑。

1 研究区概况和数据 1.1 研究区概况

研究区位于法属圭亚那热带雨林的巴拉库(Paracou)研究基地，中心纬度为5°16′N，经度为52°56′W(如图 1右上所示)。该区域常年炎热多雨，年均气温26 ℃，年均降水量2 980 mm，分雨季和旱季，旱季从8月中旬到11月中旬。地形以丘陵为主，海拔0~50 m。森林群落包括原始林和部分经不同程度采伐、干扰与自然因素影响后的天然次生林，树种繁多，每公顷森林140~160个树种(胸径>10 cm)，主要有金壳果科(Chrysobalanaceae)、大戟科(Euphorbiaceae)、橄榄科(Burseraceae)、桑科(Moraceae)、山榄科(Sapotaceae)、豆科(Leguminosae)等。森林结构复杂，树高在20~45 m之间，森林AGB在200~500 t·hm^-2之间。

1.2 多基线机载SAR数据获取与处理

多基线极化干涉SAR数据来自欧洲空间局(ESA)2009年热带林机载SAR遥感试验(TropiSAR 2009)。TropiSAR 2009是ESA BIOMASS地球探测项目计划的一部分，由ESA和法国国家空间研究中心(CNES)组织实施。机载数据采用法国国家航空航天研究中心(ONERA)研制的SETHI雷达系统，于2009年8月在法属圭亚那巴拉库研究基地获取试验数据。该数据由6轨重复飞行的P-波段全极化SAR数据组成，其中一轨极化SAR数据总功率影像如图 1左所示。主影像航高为3 962 m，空间基线以15 m间隔在垂直方向上均匀分布，时间基线为2 h，斜距向分辨率为1.0 m，方位向分辨率为1.2 m，入射角近距19°，远距52°。本研究采用数据已由ONERA进行辐射和极化定标、配准、平地相位去除等预处理，选取图 1所示红色方框区域开展研究工作，该区域主影像PauliRGB显示结果如图 1右下所示，影像大小为2 300行×1 500列。

1.3 地面数据获取与处理

法国农业发展国际合作研究中心(CIRAD)在该研究区内设立了16个固定大样地(如图 1红框范围内)，包括15个250 m×250 m样地(编号1~15)和1个500 m×500 m样地(编号16)，其中9个大样地(250 m×250 m)经历了3次不同程度采伐以进行伐后森林更新研究(Gourlet-Fleury et al., 2004)，其他样地为原始林，未经过人工干扰。对固定样地每木胸径检尺(起测径阶为10 cm)，测量部分代表性林木树高，建立树高-胸径相对生长方程(Réjou-Méchain et al., 2015)，计算获得每株树木的树高，进而以胸径、树高和密度(Chave et al., 2009)为自变量利用Chave等(2005)异速生长方程计算每木AGB，累加获得样地水平森林AGB。研究表明，Chave等(2005)提出的异速生长方程适用于本研究区：1)本研究区位于法属圭亚那热带雨林，试验区数据已参与该异速生长方程的建模研究；2)该异速生长方程已通过法属圭亚那热带雨林试验区数据验证；3)该异速生长方程基于热带雨林多个试验区(27个研究区)数据建立，估测模型较为稳健，适用于热带雨林众多森林类型的生物量估测。为增加样地数量，同时考虑样地尺度引起的建模误差(Chave et al., 2004)，本研究按子样地(共85块)进行评价，在编号1~15的大样地内选用子样地大小为125 m×125 m，在编号16的大样地内选用子样地大小为100 m×100 m。

1.4 LiDAR数据获取与处理

法国农业发展国际合作研究中心提供了覆盖研究区的机载LiDAR DEM和DSM产品，用于辅助分析多基线InSAR层析成像结果。该数据由ALTOA系统于2009年4月飞行获取，航高120~220 m，数据获取范围较小，仅覆盖SAR影像部分区域(图 1)。所获取的DEM和DSM产品由CIRAD从原始LiDAR点云数据中提取，采用WGS84坐标系，投影为通用横轴墨卡托投影(UTM)，空间分辨率为1 m，利用地面控制点对其精度进行检验，其高程平均误差为0.02 m(Vincent et al., 2012)。

2 研究方法

多基线InSAR通过多个高度不同的天线在法向(与雷达视线和飞行方向相垂直的方向)形成合成孔径而具备高程分辨能力，基于多基线InSAR层析技术可获取以层析相对反射率表达的森林垂直结构剖面，该剖面分布形状与森林空间分布形态具有较强的相似性，在一定程度上反映了森林垂直结构信息，其不同高度处层析相对反射率以不同形式表征森林AGB大小。多基线InSAR层析森林AGB反演具体流程如图 2所示。首先，基于HH极化多基线InSAR数据估计地形相位，并以此为参考对HV极化多基线InSAR数据的地形相位进行去除；然后，对HV极化多基线InSAR数据进行层析成像，并对层析结果进行地理编码；最后，提取不同高度处层析相对反射率，对其与地面样地森林AGB数据进行相关性分析，选择合适的特征建立森林AGB估测模型，并对其估测精度进行精度评价。

2.1 多基线InSAR层析成像

首先对多基线InSAR数据进行配准、平地相位去除等预处理，估计得到干涉协方差矩阵，根据成像几何计算导向矢量，采用常规波束形成法(beamforming)进行空间谱估计，得到沿高程方向分布的目标散射函数(即层析相对反射率垂直分布函数)(Krim et al., 1996)。其谱估计公式为：

$ {P_{{\rm{CBP}}}}\left(z \right) = a\left(z \right)Ra\left(z \right)。$

(1)

式中：P_CBP(z)表示待估算的层析相对反射率垂直分布函数；a(z)表示高度为z的导向矢量；R表示干涉协方差矩阵。

对各极化通道的多基线InSAR数据进行层析成像，得到层析相对反射率垂直分布信息。图 3a红线位置上像元的HH极化、HV极化和VV极化层析相对反射率垂直分布分别如图 3b-d所示，其中LiDAR DEM和DSM已由地理坐标投影转换到SAR几何坐标。由图 3可见，HH极化的后向散射能量主要集中在地表(图 3b)，HV极化的后向散射能量主要集中在冠层(图 3c)，VV极化的后向散射能量垂直分布介于HH极化和HV极化之间(图 3d)。为此，本研究基于HH极化多基线InSAR层析结果估计地形相位，基于HV极化多基线InSAR层析结果提取相关特征用于森林AGB估测。

2.2 特征提取

由图 3b可知，HH极化层析相对反射率在垂直方向的峰值高程位于林下地表，提取该峰值高程可反演得到林下地表高度。图 4为HH极化峰值高程与LiDAR DEM的剖面对比结果，其中蓝色实线为直接提取的峰值高程，红色点划线为高斯滤波后的峰值高程。由图 4可见，滤波后的峰值高程与LiDAR DEM较为吻合，因此，本研究以滤波后的峰值高程作为待反演的林下地表高度。

由提取的林下地表高度，根据式(2)估计得到地形相位，并以此为参考对HV极化多基线InSAR数据的地形相位进行去除，进而对去除地形相位的HV极化多基线InSAR数据进行层析成像：

$ {\varphi _0} = {k_z}{z_0}。$

(2)

式中：φ₀为待估计的地形相位；k_z为垂直有效波束；z₀为林下地表高度。

图 5为对应图 3a红线位置上像元的成像结果，其中黑色实线为去除地形相位后的林下地表高度(高程为0 m)，白色实线为LiDAR DSM与LiDAR DEM差分得到的冠层高度模型(CHM)。由图 5可知，地形相位去除后得到的层析相对反射率的高程以林下地表高度为高程起算面，据此可提取地表以上不同高度处层析相对反射率，作为森林AGB估测模型的输入特征。

2.3 森林AGB估测模型 2.3.1 相关性分析

对不同高度处层析相对反射率与地面样地森林AGB之间的相关性进行分析，以选择合适的特征用于森林AGB建模。本研究按5 m间隔提取地表以上不同高度处的层析相对反射率，分别对地表以上0、5、10、15、20、25、30、35、40 m高度处的层析相对反射率与森林AGB的相关性进行分析。

已有研究表明，在北方及温带森林地区(生物量＜250 t·hm^-2)，P-波段SAR后向散射强度对森林AGB有较高的敏感性(Le Toan et al., 2011)，但对于高生物量的热带雨林地区(生物量＞250 t·hm^-2)，利用后向散射强度进行森林AGB估测的研究较少。为了比较层析相对反射率相对于雷达后向散射强度对森林AGB的敏感性程度，本研究同时分析了各极化SAR后向散射强度与森林AGB的相关性。

2.3.2 建立森林AGB估测模型

选择与森林AGB相关性较强的特征，分别以其为自变量，采用一元线性回归法建立森林AGB估测模型。为充分发挥层析相对反射率在森林AGB估测中的作用，同时以多个特征为自变量，采用多元线性回归法建立森林AGB估测模型，以提高森林AGB的估测精度。

以试验区内85块样地调查数据为参考，采用留一交叉验证法(骆社周等，2015)对森林AGB估测模型进行精度分析。留一交叉验证法假设有N个样本，将每个样本作为测试样本，剩余N-1个样本作为训练样本，是评价回归模型泛化能力的一个有效方法，尤其在样本数量相对较少的情况下非常有效(Brovelli et al., 2008)。以均方根误差(RMSE)/相对均方根误差(RRMSE)、决定系数(R²)和绝对平均精度作为精度评价指标，其中RRMSE表示RMSE相对于平均森林AGB的百分比，绝对平均精度表示为：(1-|反演值-参考值|/参考值)×100%。

3 结果与分析 3.1 相关性分析

分别对HH极化、HV极化和VV极化后向散射强度与地面实测森林AGB的相关性进行分析，结果如图 6所示。由图 6可知，各极化后向散射强度与森林AGB之间的相关性很低，后向散射强度在该生物量范围内存在“饱和”现象，表明后向散射强度信息对热带雨林森林AGB的估测能力较弱。

图 7所示为以5 m为间隔，不同高度处层析相对反射率与森林AGB的相关性。由图 7可知，各高度处层析相对反射率与森林AGB相关性较高，20 m以下各高度处层析相对反射率与森林AGB呈不同程度的负相关(图 7a-d)，其中，5 m高度处层析相对反射率与森林AGB负相关系数最高(r=-0.58)；20 m以上各高度处层析相对反射率与森林AGB呈不同程度的正相关(图 7e-i)，以25 m高度处层析相对反射率与森林AGB正相关系数最高(r=0.63)。更为重要的是，本研究发现森林AGB高达500 t·hm^-2时，层析相对反射率仍未出现“饱和”现象，表明利用层析相对反射率有望实现热带雨林森林AGB的有效估测。

3.2 森林AGB估测

以上相关性分析表明，相较于雷达后向散射强度对森林AGB的敏感性，多基线InSAR层析技术得到的某一特定高度处层析相对反射率对森林AGB的敏感性更高。不同高度处层析相对反射率从不同层面间接反映了森林AGB大小，其中5 m和25 m高度处层析相对反射率对森林AGB具有较好的指示作用，分别以其作为输入特征变量，采用一元线性回归法建立森林AGB估测模型(5 m：y=-1 237.71-630.60x；25 m：y=2 231.72+744.68x)，交叉验证结果如图 8a、b所示。

可简单认为5 m和25 m高度处层析相对反射率分别代表来自森林底部和顶部的后向散射信号，同时利用森林底部和顶部的信息联合建模有望提高森林AGB的估测精度，故同时以5 m和25 m高度处层析相对反射率为输入特征变量，采用多元线性回归法建立森林AGB估测模型(y=957.19-286.45x₁+526.83x₂，其中，x₁为5 m高度处层析相对反射率，x₂为25 m高度处层析相对反射率)，交叉验证结果如图 8c所示。

由图 8可见，采用5 m高度处层析相对反射率进行建模的估测精度为88.44%，RMSE为49.85 t·hm^-2 (RRMSE为13.56%)，R²为0.31；采用25 m高度处层析相对反射率进行建模的估测精度为88.82%，RMSE为47.30 t·hm^-2(RRMSE为12.87%)，R²为0.37；进一步同时采用5 m和25 m高度处层析相对反射率进行联合建模，其估测精度为89.17%，RMSE为46.45 t· hm^-2(RRMSE为12.63%)，R²为0.40，优于仅采用5 m或25 m高度处层析相对反射率进行建模的精度。由此可见，利用特定高度处层析相对反射率可反演得到较高精度的森林AGB，且利用不同高度处层析相对反射率联合估测可进一步提高森林AGB估测精度。

4 讨论

层析相对反射率与森林AGB的相关性受波长、地形、森林类型、树种组成、空间结构、生物量水平等因素影响，针对P-波段SAR在热带雨林的多基线InSAR层析结果，本研究区域利用5 m和25 m高度处层析相对反射率建模效果较好，但尚未对其他区域进行对比分析与评估。不同高度处层析相对反射率与森林AGB的相关性也会因林分状况(如稀疏、低矮植被等)而有所差异，方法适用性仍需进一步深入研究。此外，样地间的空间相关性、大小以及外业调查的不确定性因素等也会影响森林AGB的估测精度。尽管如此，本研究所发展的森林AGB多基线InSAR层析估测法有效解决了热带雨林森林AGB遥感估测方法的信号“饱和”问题，可满足热带雨林区域森林生物量制图的应用需求，为区域及全球森林生物量估测与碳储量研究提供关键技术支撑。

本研究更多的是对层析技术和森林地上生物量估测方法的探讨，仅选用了冠层顶部和底部各一层的相对反射率参数特征进行生物量建模，有待利用更多层垂直结构信息进行多元非线性建模。不同于SAR散射机制的激光雷达可能会进一步完备表达森林垂直结构，其探测参数与层析反演参数的联合也可能进一步提高森林地上生物量的估测精度。另外，层析方法本身的优化仍有待进一步研究，以期建立更加精细(“高分辨率”)和完整的森林空间结构参数表达。

5 结论

层析技术可提取森林垂直方向空间结构信息，有助于深刻理解SAR信号与森林内部散射体相互作用的物理机制，同时森林的垂直结构可有效表征森林AGB的大小。本研究利用多基线InSAR层析技术对热带雨林进行三维成像，通过分析不同高度处层析相对反射率与森林AGB之间的相关性，提出了联合不同垂直高度处多基线InSAR层析特征估测森林AGB的方法，利用实测样地AGB数据验证了其有效性，结论如下：

1) 相较于雷达后向散射强度对森林AGB的敏感性，多基线InSAR层析技术得到的某一特定高度处层析相对反射率对森林AGB的敏感性更高。20 m以下各高度处层析相对反射率与森林AGB存在负相关，20 m以上各高度处层析相对反射率与森林AGB存在正相关。该现象为P-波段SAR信号在热带雨林条件下与森林散射体相互作用的结果，可能由于热带雨林森林散射体对SAR信号具有较强的消光作用，且消光系数随着森林AGB的增大而增大，使得森林AGB越大，P-波段SAR信号与森林上层部分的相互作用越强，而与森林下层部分的相互作用越弱。

2) 层析技术提取的不同层相对反射率与森林中的不同散射机制相对应，不同层相对反射率也以不同形式表征森林AGB大小。分别利用5 m和25 m高度处层析相对反射率可反演得到较高精度的森林AGB，且在森林AGB高达500 t·hm^-2时没有出现“饱和”现象，说明多基线InSAR层析技术得到的特定高度处层析相对反射率对热带雨林森林AGB具有较好的指示作用。

3) 5 m和25 m高度处层析相对反射率分别代表来自森林底部和顶部的后向散射信号，同时利用5 m和25 m高度处层析相对反射率联合估测进一步提高了森林AGB的估测精度，说明不同层相对反射率联合建模可充分发挥森林垂直结构在森林AGB估测中的作用，充分利用不同层次森林垂直结构信息可进一步提高复杂空间结构条件下森林AGB的估测精度。