利用ESAR极化数据的复杂地形区森林地上生物量估算

张海波,汪长城,2，朱建军，付海强

1. 中南大学地球科学与信息物理学院，湖南 长沙 410083； 2. 中南大学有色金属成矿预测教育部重点实验室，湖南 长沙 410083

森林覆盖了地球陆地30%左右的面积，在维持陆地生态平衡、保护生态安全、防止生态危机方面起着决定性的作用[1-2]。同时，森林也是陆地上最经济的吸碳器和最大的碳储库[3],森林生物量表征着森林生态系统中的碳储量[4]，因此，精确估测生物量特别是大尺度范围估测，能使人们更好地理解全球碳循环和全球变暖。一般来说，生物量包括地上生物量(above-ground biomass，AGB)和地下生物量，由于收集地下生物量的实测数据比较困难，因此，当前的研究大多针对地上生物量[5]。

目前森林生物量估算的方法主要有传统的地面调查法和遥感监测手段。传统的地面调查方法可直接获得单木树高、胸径、冠幅等参数，通过异速生长方程可获取森林生物量[6]。这种方法需要耗费大量的人力、物力和财力，对森林破坏性大，并且只能局限于小范围，不适用于区域、全国乃至全球尺度的生物量估测与动态变化监测[7-8]。与之相比，具有大范围、快速、准确的遥感技术成为了获取地上生物量的主要途径，遥感数据已广泛应用于区域尺度植被动态变化监测研究[9-10]。其中，雷达遥感因其具有全天候、强穿透力、对地表粗糙度敏感等优于光学遥感的特性，更有利于生物量的估测[11]。

近年来，利用雷达遥感数据估测生物量取得了许多进展，一些研究表明长波长SAR(L、P波段)数据后向散射系数在评估地上生物量方面有着相当大的潜力[12-14]。但是，由于SAR特殊的侧视成像方式，其后向散射容易受复杂地形的影响，其中明显的影响是降低了雷达的频率[15],而且地形坡度在SAR影像中是很难处理的。目前在降低复杂地形对SAR影响方面的研究主要是根据面积积分[16]或利用投影角进行地形辐射校正[17-19]。该方法在以单次和二次散射为主的非植被区域能较好地校正地形引起的后向散射异常，但对于结构复杂、内部存在多次散射现象的森林区域，难以描述地形变化引起的多次散射的变化，因此对于森林参数定量分析改善不明显。在森林AGB估测模型中融入地形因子的方法，不但考虑了地形对后向散射的综合影响，而且不需要考虑地形对森林散射机制的具体影响[20]。目前，该方面的研究相对较少。

以往的研究中，缺乏地形因子对估算模型估算森林AGB的能力及稳定性影响的分析，地形变化对后向散射强度与森林AGB的响应关系的影响同样缺乏可靠分析。因此，本文利用机载P波段全极化SAR数据，分析地形变化与不同极化后向散射的关系，找到相关性最高的极化方式用于森林AGB估算，建立适用于地形复杂区森林AGB估算的模型，并通过不同模型的比较，验证模型的适用性与稳定性，同时利用该模型对研究区森林AGB进行估算。文中的分析方法除了能适用山地、丘陵地区森林AGB估测之外，对地形复杂区的农作物估产方面同样具有一定的借鉴意义。

1 材料与方法

1.1 研究区概况

研究区位于瑞典北部西博滕省温德恩市(Vindeln municipality)的Krycklan流域，处于64°11′E—64°17′E，19°41′N—19°53′N之间，是温德恩市森林实验场的一部分，同时也是瑞典大学野外森林研究基地，面积约为9390 hm2，属于亚寒带针叶林气候。林区地形起伏较大，海拔高度范围为100～400 m。森林类型为自然生长的混交林，其中大多数为针叶林树种，包括云杉(Spruce)，松树(Pine)，落叶阔叶林树种有桦树(Birch)，除此之外还有其他少量的落叶乔木如山杨(Aspen)、花楸(Rowan)等。研究区如图1所示。

图1 研究区及林分分布Fig.1 The location of test site and the distribution of stand

1.2 森林实测数据

森林实测数据来自于研究区24个林分采样区，样区分布见图1，面积为3.07～24.34 hm2，在每个林分中根据林分面积按照一定的间距(50～160 m)布置8到13个半径为10 m的样地，样地总数为310个，样地调查时间为2008年10月13日—17日。在这些样地中对胸径大于4 cm树木，逐一测定单木的树种、树高和胸径，并记录其他一些参数，如植被种类、土壤类型等。各样地中枯死的树木及落叶不参与生物量计算，因此不予记录。实地测量的单木数据根据不同林分区域的树种类型先汇总成样地尺度的森林参数，包括每个树种的断面积加权平均年龄、胸高断面积加权平均直径(cm)、胸高断面积加权平均树高(dm)。不同树种AGB信息参考H. Perersspm的计算方式，将每个树种的生物量分成3个组分(树干、树枝、树叶)分别计算。由于季节原因，落叶林不计算树叶生物量，落叶乔木归于桦树一类统计。

1.3 DEM数据

DEM数据来源于欧空局BIOSAR项目，不包含植被信息。首先对其进行插值和重采样，以获取规则的样地DEM数据，处理之后的数据分辨率为1 m×1 m，然后提取出研究区坡度。由于SAR主要受距离向的坡度影响[21]，因此参照飞行方向，规定坡面朝向雷达视线方向的坡度为正，背向雷达视线方向的坡度为负。雷达入射角、局部入射角和地形坡度的几何特征及关系如图2所示，其中θi是雷达入射角，θl是局部入射角，θs是坡度。

图2 坡度几何特征Fig.2 Geometry characteristics of terrain slope

局部入射角可由地形坡度与雷达入射角获得，其计算公式如下

θi=arccosH/SR

(1)

θl=θi-θs

(2)

式中,SR为斜距距离，是斜距方向上近距起点的斜距值加上像素数乘以斜距向采样间距，起点斜距为4 294.21 m；H是飞机的飞行高度。研究区坡度与局部入射角见图3。

图3 研究区坡度和局部入射角Fig.3 The slope and local incidence angel of the study area

1.4 SAR数据

SAR数据来源于欧空局BIOSAR项目，是由德国宇航局(DLR)应用机载E-SAR传感器获取的P波段全极化数据(HH、HV、VH、VV)，数据支持平台高度为4 090.1 m，飞行速度为92 m/s，中心频率为0.35 GHz，天线俯视角为40°，航向角为-47.2°，脉冲发射频率(PRF)为2000 Hz。影像成像时间为2008年10月14日，数据类型为单视复数影像(SLC)，影像的重叠度为50%。在经过多视处理(2×1)、滤波、地理编码，重采样等SAR影像基本处理之后，最终的影像坐标系统为UTM WGS-84，像元大小为1 m×1 m，SAR影像的基本处理由DLR完成。根据雷达系统的互易性，HV极化方式与VH极化方式近似相等，因此，本文从HH、HV、VV这3种极化方式的影像数据中获取后向散射系数，其获取公式如下

(3)

式中，i为P波段不同极化方式的SAR影像(HH、HV和VV)；m为林分中的像元；n为林分；k和c是校正常量，分别为32 769和10-6(文献[22])；sinθ是几何投影因子；θ为雷达视角，式(3)可通过10lgσ0将其转换为dB值。

1.5 估算模型

在对植被生物量的估算中，半经验模型往往能够提供清晰的结构和思路来理解不同后向散射相互作用的物理机制，能够对研究对象提供合理的解析[23-24]。本文采用的半经验模型是以改进的水云模型为前提，公式如下

(4)

(5)

式中，β为需通过经验定义的系数；V为蓄积量。将式(5)代入式(4)可得式(6)

(6)

对式(6)进一步参数化，由于文中估算的参数为森林AGB，因此，将公式中的V用AGB替换，如下

(7)

本文在上述分析模型的基础上引入局部入射角，得到最终模型(水云分析模型)如下

σ0cos-1(θi-θs)=β1+β2eβ3AGB

(8)

为验证文中水云分析模型在地形复杂区估算AGB的适用性，文中使用顾及地形因子的线性模型、对数模型、二次模型、指数模型与水云分析模型对比分析。为方便表述，令γ0=σ0cos-1(θi-θs)，各模型具体形式见表1。

表1 森林AGB估算模型

1.6 模型参数确定方法

构建森林AGB估算模型后，需要通过样地实测数据与后向散射系数值确定模型的参数值。由于使用的估算模型多为非线性模型，难以用传统方法确定模型的参数值，而遗传算法(GA)适用范围较广，无论非线性问题能不能转化为线性问题，都可以直接进行数据拟合，并且对于可以转化为线性问题的非线性数据问题，遗传算法要优于其他算法[27]，因此本文使用遗传算法来确定生物量估算模型的最优参数值。遗传算法采用了“适者生存”的原则来获得最佳解决方案，基本操作包括了选择、交叉和变异，染色体的基因对应于模型的待求参数，进化过程中的适应度函数为实测生物量值与模型计算生物量值之间的误差，具有最小误差的参数组合表示能很好地适应环境。本文方法流程见图4。

图4 方法流程Fig.4 Methodology flowchart

2 结果与分析

2.1 不同极化方式后向散射与森林AGB响应关系

研究区不同极化方式的后向散射系数值与对应生物量值的相关性分析结果见表2。由表2可知，在3种极化方式中HH与HV极化方式后向散射系数与生物量的相关性要高于VV极化方式，其中HV极化方式后向散射系数与生物量呈强相关性(0.6～0.8)，相关性达到0.696，而VV极化方式后向散射系数与生物量的相关性最弱，为-0.052。HH和HV极化方式后向散射与生物量呈正相关，表明在HH和HV极化方式下，后向散射系数值随着地面生物量的增加而增加，VV极化方式后向散射与生物量呈负相关，表明在VV极化方式下，后向散射系数随着生物量的增加呈减少趋势。这种现象与地面植被覆盖的情况相关，在植被覆盖较高且为大型植被的区域，交叉极化方式对生物量的敏感程度高于单极化方式，在植被覆盖较低的区域，则会出现相反的情况，尤其对于VV极化方式最为明显，其在高覆盖的大型植被区对生物量最不敏感，然而在稀疏低矮的植被区对生物量的敏感性较好。

表2 后向散射系数与森林AGB的相关性

注：* *表示在0.01水平显著相关。

从不同极化方式后向散射系数值的动态范围来看，极化方式不同，对应的后向散射系动态范围不同(图5)。HV极化方式与HH、VV极化方式相比要集中，动态范围为3.94 dB(-16.98～-13.04 dB)，而HH、VV极化方式分别为9.43 dB(-11.41～-1.98 dB)和5.85 dB(-10.79～-4.94 dB)。相同极化方式中不同生物量水平对应的后向散射范围不同，在HV极化方式中，生物量低于120 t/hm2，其后向散射系数的动态范围为3～4 dB，生物量高于120 t/hm2，后向散射系数动态范围仅为1～2 dB。由图5还可知，HH、HV和VV极化方式的后向散射系数变化趋势在森林AGB处于较低值的情况下能与AGB变化趋势保持一致，但随着AGB值的增大，这种一致性仅在HV极化方式下继续保持，HH和VV极化方式均出现高AGB值对应低后向散射系数值，并且比低AGB值对应的后向散射系数值还要低的情况。

通过以上对比分析可知在不同极化方式中，HV极化方式更适用于复杂地形区的森林AGB估算，因此本文后面的分析主要针对HV极化方式展开。

图5 后向散射系数与地上生物量散点图Fig.5 The scatter plot between backscatter coefficient and AGB

2.2 估算模型不同坡度的适用性

为了更深入了解不同地形对估算模型的影响，以及对比不同模型在不同坡度情况下的适用性，本文将坡度分为3个等级，分别为0～5°、5°～10°、≥10°。每个等级对应的林分样地数目分别为24、18和11个。分级后的样地具体情况见表3。

表3 不同坡度等级样地情况

由图6可知，在坡度小于5°的地区，5种估算模型决定系数分别为0.443、0.539、0.552、0.501、0.57，均为3种坡度等级中最高，随着坡度的增加，决定系数均有下降，在坡度≥10°的区域，5种模型的决定系数达到最低，分别为0.13、0.197、0.32、0.16、0.424。这表明在复杂地形区，地形因子对森林AGB的估算具有极大的影响，后向散射系数与森林AGB的相关程度随地形坡度的增加而减小。这可能是由于在地形起伏较小的区域，电磁波与植被层的相互作用构成的“水云层”在整个区域较稳定，没有发生实质的变化，而随着地形坡度的增加，降低了雷达频率，“水云层”结构发生改变，电磁波与地面的相互作用增加，这些易使电磁波与植被之间的相互作用出现异常。

在不同坡度等级中，5种模型估算AGB的决定系数大小排序是相同的，均为水云分析模型>二次模型>对数模型>指数模型>线性模型，这表明水云分析模型在复杂地形区估测生物量的适应能力最强，线性模型适应能力最弱。不同坡度等级变化中，各模型估算森林AGB的决定系数变化值见表4，坡度由0～5°变化至5°～10°中，水云分析模型的决定系数变化值最小为0.066，线性模型最高为0.153，变化值大小排序为线性模型>指数模型>对数模型>二次模型>水云分析模型。坡度由5°～10°至≥10°中，变化值最小的为水云分析模型，为0.08，最大值为对数模型，为0.206，不同模型变化值大小排序为对数模型>二次模型>线性模型>指数模型>水云分析模型。坡度由0～5°至≥10°中，变化值最小的为水云模型，变化值为0.146，变化值最大的为对数模型，为0.342，不同模型变化值大小排序为对数模型>指数模型>线性模型>二次模型>水云分析模型。这表明，无论在地形起伏较小的地区还是地形起伏较大的地区，水云分析模型估算森林AGB的可靠性和稳定性均高于其他4种模型。

表4 不同坡度过渡中决定系数变化值

注：(a)、(b)、(c)、(d)、(e)分别指代图6中对应的模型

2.3 森林AGB估算结果

由2.2节分析可知，不同坡度等级内水云分析模型是5种模型中最适用于复杂地形区森林AGB估算的。然而在实际应用中，如区域性森林AGB估算和生物量制图，一个林分样地往往包含了多种坡度等级，其后向散射系数值是多类坡度下后向散射的结果。因此为完成研究区生物量制图，先利用研究区24个完整林分的实测生物量值与对应的HV极化方式后向散射系数值，通过遗传算法得到估算模型最优系数值，结果见表5。

图6 不同坡度中模型拟合结果比较Fig.6 Comparison of models fitting results in different slopes

表5 森林AGB估算结果

由表5可知，在整个区域中，水云分析模型估算生物量的决定系数最高，为0.796，均方根误差(RMSE)最小，为0.466，说明水云分析模型相对其他4种模型最适用于整个研究区生物量反演和制图，为便于理解及生物量估算结果检验，水云分析模型可写成如下形式

(9)

确定适合复杂地形生物量估算模型之后，对模型精度检验是必要的。首先将林分对应的HV极化方式后向散射系数值代入模型中，估算出生物量值，然后计算出实测生物量值与模型估算生物量值的均方根误差，结果见图7。实测生物量与模型估算的生物量值决定系数为0.597，RMSE为30.876 t/hm2，拟合精度为77.40%。最终的AGB制图结果见图8，研究区生物量分布范围大致为：0～260 t/hm2，平均生物量为76.01 t/hm2。

图7 森林AGB估算结果精度评价Fig.7 Accuracy of the estimated forest AGB

3 结 论

本文分析了复杂地形区不同极化方式后向散射系数与森林AGB的相关性，首次将地形因子引入改进的水云模型，并通过其他模型与之对比分析，采用遗传算法确定模型的最优参数，并对模型在不同坡度情况下的可靠性与稳定性进行了对比分析，结果显示：

3种极化方式后向散射系数变化趋势在森林AGB处于较低值的情况下能与AGB变化趋势保持一致，但随着AGB值的增大，这种一致性仅在HV极化方式下继续保持，因此相比之下，HV极化方式更适用于复杂地形区生物量的估算。

图8 森林AGB分布图(t/hm2)Fig.8 Distribution maps of the forest AGB

地形对森林AGB的估算具有极大的影响，后向散射系数与AGB的相关性随着地形坡度的增加而减小。在地形坡度为0～5°、5°～10°和≥10°内，5种模型估算生物量的能力大小排序均为：水云分析模型>二次模型>对数模型>指数模型>线性模型。

模型在复杂地形区估算生物量的稳定性随坡度变化有所改变，当坡度由0～5°变化至5°～10°时，水云分析模型的决定系数变化值最小,线性模型最高，稳定性大小排序为：水云分析模型>二次模型>对数模型>指数模型>线性模型；坡度由5°～10°至≥10°中，变化值最小的为水云分析模型，最大值为对数模型，稳定性排序为：水云分析模型>指数模型>线性模型>二次模型>对数模型；坡度由0～5°至≥10°中，变化值最大的为对数模型，稳定性排序为：水云分析模型>二次模型>线性模型>指数模型>对数模型。虽然模型稳定性排序因坡度变化而不同，但对于水云分析模型而言，无论在坡度起伏较小的地区还是较大的地区，其稳定性都最高。

利用本文所述融入地形因子的水云分析模型估算的森林AGB与实测AGB值进行精度检验，其决定系数为0.597，RMSE为30.876 t/hm2，拟合精度为77.40%，比报告中使用Le Toan方法的精度有所提高(RMSE为35 t/hm2)，估算出研究区平均生物量为76.01 t/hm2。

本文采用林分实测数据用于生物量拟合及验证生物量估测模型，虽然林分样地的数量有限，尤其是将坡度分级后，在坡度较大的区域，实测数据只有11个，但是使用实测数据在一定程度上能够确保模型系数的正确性及检验的准确性。另外，文中试验数据为P波段数据，其波长长，穿透性强，在植被区信号能穿透植被直接与地面交互，这种情况下使用本文方法能够明显改善地形对后向散射的影响，提高生物量的估算精度。而对于短波波段，由于波段短，穿透能力弱，其信号主要与植被冠层交互，这种情况下，使用在模型中融入地形因子的方法对生物量估算精度的改善可能不明显。本文研究方法除适用于山地、丘陵地区森林AGB估测之外，对地形复杂区的农作物估产方面同样具有一定的借鉴意义。