基于多角度卫星遥感的胡杨林结构参数获取

杨雪峰 木尼热·买买提

（新疆师范大学地理科学与旅游学院 新疆干旱区湖泊环境与资源实验室 乌鲁木齐 830054）

森林是一个复杂的统一体，由多种生物和非生物环境组成，其彼此之间相互影响、相互制约，表现出一定的结构特征和相互作用规律。森林结构特征是森林植物与周围环境条件之间以及森林植物彼此之间相互作用的表现形式，可直观反映森林生长状况，既是林学、生态学和地学领域的主要研究对象，也是众多陆地生态模型的重要输入参数，获取大尺度森林结构特征具有重要科学价值（Pedroni, 1997；Schröteret al.,2019）。

森林结构特征通常采用覆盖度（fractional vegetation cover，FVC）、林分密度（density）、树高（height）、冠幅（crown width）、胸径（diameter at breast height）、叶面倾角分布（leaf angle distribution，LAD）和叶面积指数（leaf area index，LAI）等指标表示（Rogasset al., 2012）。目前，基于遥感手段获取森林结构参数已取得较大发展，能够获取大面积森林结构参数的遥感技术主要有激光雷达（light detection and ranging，LiDAR）、雷达、光学遥感等。机载LiDAR是测量植物冠层结构、形状和地上生物量的有效探测手段（董立新, 2016；Lefsky, 2010；Omasaet al., 2002；van Leeuwenet al., 2010）；但因LiDAR机载设备观测范围较小，使用成本较高，星载LiDAR设备地面采样点又十分稀疏，限制了LiDAR的应用。干涉合成孔径雷达（interferometric synthetic aperture radar，InSAR）和极化干涉合成孔径雷达（polarimetric SAR interferometry，PolInSAR）等技术发展出随机体散射体/地表二层（random volume over ground，RVoG）、极化相干层析（polarization coherence tomography，PCT）等模型方法，可用来获取森林结构的精细结构特征（Cloudeet al.,2003；Cloude, 2006；Parket al., 2015；谈璐璐等, 2010）；然而，现阶段星载极化干涉数据面临时间失相干、空间失相关等影响，其应用受一定限制，甚至无法取得令人信服的反演结果（李震等, 2014；刘晓萌等, 2007；李廷伟等, 2009；Leeet al., 2012）。随着高分遥感技术的发展，利用影像高分辨率特点获取树木结构参数的研究日益增多；但由于森林阴影的长度、太阳高度角、方位角以及坡度、坡向等多个因子具有不确定性，使用高分影像在测量中容易造成较大的人为误差，无法推广到大范围森林参数测量中（崔少伟等, 2011；韩学锋, 2008；周志强等, 2012；Keet al., 2011）。与单一方向光学遥感相比，多角度光学遥感不仅利用多光谱观测提取地物组分的波谱信息, 而且利用卫星前向、后向和星下点不同角度的观测，既富含地表结构信息，同时较多的观测角度也为反演求解地表结构参数提供了便利。随着多角度光学遥感的发展, 多种二向反射分布函数（bidirectional reflectance distribution function，BRDF）模型可用来反演并估计地表植被类型（杨雪峰等，2017），采用物理模型结合多角度观测数据提取植被结构参数具有很大潜力（Suet al., 2007；王亚楠,2010；金晟业, 2009)。目前的星载多角度传感器中，多角度光谱辐射计（multi-angle imaging spectro-radiometer，MISR）观测角度最多，通常采用MISR数据结合几何光学模型（geometric-optical model，GOM）和基于核驱动模型的AMBRALS （algorithm for model bidirectional reflectance anisotropies of the land surface）算法对目标进行区分并反演推测相应的植被结构参数（Dineret al., 1998）。鉴于MISR独特的大范围多角度观测能力以及免费获取政策，与 LiDAR和 PolInSAR 相比，采用MISR进行森林结构参数提取，在数据获取能力、后期数据处理便利性和使用成本方面更具优势。

在干旱区环境下，当林上冠层郁闭度不高时，直射光通过孔隙直接到达林下灌、草植被，此时林下灌、草以及较亮的土壤背景反射很大程度影响传感器接收的信号，传感器接收的信号是可见光与植被（乔、灌、草、土壤）之间非常复杂的相互交互过程；如果不能正确分离树木、林下植被和土壤的光学信号，则反演森林结构参数的任务很难实现（黄健熙等, 2005；Rahmanet al., 2004）。GOM结构未能表达背景反射的二向异性特性，也未充分利用 MISR 的多角度观测优势，因此一些研究对GOM进行修改以使其适用于干旱区环境，如简单几何光学模型（simple geometric-optical model，SGM）（Choppinget al., 2005；2008；2009）。SGM在模型中加入Walthall土壤反射模型（Walthallet al.,2000）描述林下土壤背景二向异性反射，Walthall模型参数由定标站点回归计算得到，SGM理论上更适合森林郁闭度较低的环境下使用，但模型适用性未得到广泛验证。

对塔里木河下游而言，胡杨（Populus euphratica）的高度、密度、冠幅和群落覆盖度既是植被恢复的标志性指标，也是研究胡杨群落生态过程与演化机制、改善优势种群保护措施的途径（陈海燕等，2015）。目前，该区域有关胡杨结构参数的研究工作主要包括：利用人工地面调查方式建立胡杨林植被动态变化GIS数据库（玉米提·哈力克等, 2008）；使用QuickBird、WorldView等高分卫星遥感影像数据，采用人机交互方法提取胡杨冠幅、株数等森林结构数据（牛婷等,2008；万红梅等, 2011)；基于Landsat、Hyperion等中等分辨率卫星遥感影像数据，利用像元二分模型分解反演塔河下游植被覆盖度（黄粤等, 2013）等。总体来看，以往研究评价指标较单一，除覆盖度外，研究范围局限于样点和样带等小范围区域，尚缺乏大范围、系统的胡杨林结构参数提取研究。鉴于此，本研究以MISR为主要数据源，以二向反射模型为理论基础，采用SGM模型，在无人机摄影测量技术支持下对塔里木河下游胡杨林主要结构参数（树高、冠幅、林分密度和覆盖度）进行反演和评价。

1 研究区概况

塔里木河地处我国西部内陆干旱荒漠地区，既是干旱区生态环境变化的敏感地区，也是我国生物多样性保护和全球变化研究的关键区域之一。在人类对自然水资源时空格局改变为主要形式的扰动影响下，塔里木河下游以胡杨林为主要建群种的自然植被受到严重影响，生物多样性受损，沙漠化过程加剧。自2000年以来，有关部门实施向塔里木河下游应急输水工程，以恢复下游生态系统，胡杨是应急输水生态恢复的目标植物之一。

研究区位于若羌县英苏附近的一处河段，87°59′36″—88°11′48″ E，40°22′51″—40°27′32″ N，面积约100 km2（图1），平均海拔840 m，地形较为平坦。属暖温带大陆性干旱气候，平原地区年均气温10 ℃，年均降水量约55 mm，降水匮乏，植被稀少。河岸附近发育有河岸林，主体为胡杨，表层土壤为粉砂土，且盐碱化程度较高。地下水水位较高区域分布有柽柳（Tamarix chinensis）灌丛、黑枸杞（Lycium ruthenicum）灌丛等植被。

图1 研究区Fig. 1 Study area

2 数据与方法

2.1 研究数据

2.1.1 MISR数据 MISR提供9个角度的观测信息，分别为4个前向观测角AF（26.1°）、BF（45.6°）、CF（60.0°）、DF（70.5°），4个后向观测角AA（26.1°）、BA（45.6°）、CA（60.0°）、DA（70.5°）以及1个天顶角AN（0°）。每个传感器均有蓝光、绿光、红光和近红外4个波段。MISR 9个相机中，AN相机4个波段为275 m空间分辨率数据，其余8个角度相机红光波段为275 m空间分辨率数据。

MI1B2T、MI1B2GEOP、MIL2ASAE 3种产品是目前最常用的数据。MI1B2T为TOA辐射数据，MIL2ASAE为17.6 km空间分辨率气溶胶数据，MI1B2GEOP为17.6 km空间分辨率的太阳天顶角、太阳方位角以及观测天顶角和观测方位角数据。本研究用数据获取时间为2014年8月25日，数据轨道号078111。使用MI1B2T 9相机红光波段和MI1B2GEOP数据构建多角度观测数据集。

2.1.2 无人机影像 为了对MISR影像数据提取的胡杨对象进行高度定标和数据验证，在研究区设置3个无人机验证区。2018、2019、2020年7月，在观测区使用大疆精灵4 Pro/RTK进行无人机倾斜摄影测量，相机搭载1英寸2 000万像素CMOS传感器。基本过程为：规划航线，设置飞行高度110 m， 航向重叠度80%、旁向重叠度75%，通过5向倾斜摄影方式获取研究区地面影像数据。

2.1.3 地面实测数据 2018年7月，在研究区进行胡杨结构参数实地调查，位置如图2。

图2 地面采样示意Fig. 2 Schematic diagram of field-measurement

树高和冠幅调查方法：选取49株不同生态型胡杨，应用激光测距仪的测高功能从3个不同方向测量树冠最高点，获取胡杨树高，取平均值作为实测树高；使用长卷尺分别获取相隔60°的3个方向树木冠幅，取平均值；利用OruxMaps地图软件记录每株树坐标，并在地面分辨率0.5 m的WorldView-2多光谱影像地图上进行位置标注，方便后期与无人机影像比对。

林分密度调查方法：利用OruxMaps地图软件，以WorldView-2 影像为参照，在地图中设置 25个100 m×100 m 采样格网，实地对每个格网中的胡杨进行株树统计，并在地图软件中标记。

树高、冠幅和林分密度实测值描述统计如表1 所示。

表1 采样点的描述统计Tab. 1 Descriptive statistics of sampling points

2.2 数据处理与分析

2.2.1 无人机影像处理 采用Pix4d Mapper摄影测量软件对无人机影像进行自动配准、空三加密等处理后生成点云数据（点云密度为74点·m-3），点云数据经分类得到地表和植被点云，再插值生成分辨率为0.5 m的数字表面模型（digital surface model，DSM）和数字高程模型（digital elevation model，DEM），进而得到冠层高度模型（canopy height model，CHM）（图3）。CHM数据首先使用eCognition软件的“Min/Max Filter”最小/最大值滤波算法（主要参数：Mode: Diff. brightest to center；2D kernel Size: 7）生成种子点，根据种子点位置提取树高，然后使用“Pixel-Based Object Resizing”增长算法（主要参数：mode: Growing；Pixel Layer Constraint: 2）生成植被冠层矢量多边形对象，继而通过ArcGIS计算、统计出冠幅和密度。提取与实测49株胡杨样本位置对应的胡杨对象，二者树高、冠幅比较发现，树高R2为0.90，RMSE为0.45 m；冠幅R2为0.84，RMSE为0.68 m；对与25个采样格网重叠CHM部分分别提取胡杨对象数量，与格网内实测统计株数比较，R2为0.94，RMSE为4.25株·hm-2。这表明，无人机获取的胡杨林结构数据可代替地面调查数据，作为SGM的定标参数以及最终精度验证。

图3 无人机CHM及MISR样地Fig. 3 CHM obtained by UAV and MISR samples extenta. 无人机采样区1 1st sample area of UAV b. 无人机采样区2 2nd sample area of UAV c. 无人机采样区3 3rd sample area of UAV

2.2.2 MISR影像处理 1） 预处理 基于MI1B2T TOA BRF辐射数据，对9个相机的红光波段TOA反射率采用SMAC（simplified method for the atmospheric correction）大气校正方法（Rahmanet al., 1994）转换为地面反射率。通过提取研究区影像数据并进行投影转换，与无人机影像配准后提取重叠区域的MISR像元大小（275 m2）样本47个（图3），随机选择30个作为训练集，剩余17个作为测试集。

2） AMBRALS核驱动模型反演 AMBRALS（Wanneret al., 1995）是一种半经验核驱动线性模型，其将场景反射视为各向同性反射、体散射和表面散射的组合，数学形式表示为：

式中：R为反射率；kvol为体散射核；kgeo是几何光学核；fiso、fvol和fgeo分别为各向同性反射、体散射和表面散射的权重系数；i、v、φ分别为光照天顶角、视天顶角和相对方位角。

fiso、fvol和fgeo是模型的核心参数，可反映出对应场景基本的反射特征，本研究使用这3个系数构建场景土壤背景回归方程。

根据研究区大部分地区林地灌木稀疏、草本稀薄的特点，反演选择RossThin核为体散射核，LiSparseModis核为几何光学核。

构建多角度反射数据集，使用AMBRALS反演得到47个样本区域对应的fiso、fvol和fgeo数据（图4）。

图4 SGM模型反演流程Fig. 4 The inversion flow chart based on SGM model

3） SGM 模型反演 SGM作为一种物理模型，使用不同视角观测的土壤背景反射和植被反射2部分的组合表示整个场景反射，其数学表达为：

式中：R为反射率；Gwalthall表示采用Walthall模型模拟土壤背景反射；CRoss表示采用RossThin冠层体散射模型模拟植被反射；kG和kc分别表示可见光照地面和可见光照树冠比例；ϑi、ϑv、φ表示光照天顶角、视天顶角和相对方位角。

其余部分数学表示如下：

式中：c1～c4为Walthall模型系数；ϑi、ϑv、φ为变换后的光照天顶角、视天顶角和相对方位角；ξ为光照天顶角、视天顶角和相对方位角的三角函数。

式中：λ为林分密度；γ为树冠水平半径；ε′为变换后的散射相位角；O为视野和照明阴影部分重叠区域。

SGM模型反演分为3个步骤：

1） 训练集样本的多角度反射率观测数据和森林几何参数作为已知固定变量输入SGM模型，使用GRG（Facóet al.,1989）优化算法反演得到Walthall模型系数（图4A）；

2） 使用训练集的AMBRALS模型系数和AN相机多光谱反射率数据，建立Walthall模型系数的回归方程（图4B）；

3） 使用Walthall模型系数回归方程获得测试集土壤背景反射，SGM和GRG算法反演获取测试集的几何参数（图4C）。

3 结果与分析

3.1 背景反射模拟及精度分析

塔里木河河岸林郁闭度低，高亮土壤背景与暗深河水对比反差较大，再加上场景内灌丛和草本植被，SGM模型能否准确模拟土壤背景反射尤为重要。

反演第一阶段使用SGM模型将观测反射强度分解为土壤背景反射和冠层反射2部分，背景反射模拟采用Walthall模型，冠层反射模拟采用RossThin模型。输入数据包括无人机获取的训练集样地结构参数（树高、覆盖度、冠幅、林分密度）、模型固定参数LAI（固定值为2.3）（巴比尔江·迪力夏提等，2012）、HB（固定值为2.0）、BR（固定值为1.0）（图5）。

图5 SGM模型树冠形态参数Fig. 5 The crown shape parameters of SGM model

反演以观测反射值与模拟反射值的均方根误差（root mean square error，RMSE）最小化为目标函数（式7）（图6a），最终输出Walthall的c1～c4参数值：

图6 背景反射模拟与评价Fig. 6 Simulation and evaluation of Walthall model parametersa. 背景反射模拟Simulation of background reflection；b. 背景反射评价Evaluation of background reflection.

以上处理过程使用随机选择的30块训练样地，模拟反射值与观测反射值R2最大为0.99，最小为0.72，均值为0.92，大部分样地的反射值R2为0.90以上（图6b）。R2较低和RMSE较大的多为河床附近、覆盖度较大的样地，原因可能因河水背景影响，以及胡杨呈现出聚生、丛生状态，在此情形下，SGM模拟误差较大。

对训练样地反演获取的树高和覆盖度与参考值进行比较发现，覆盖度与参考值的一致性较好，R2达0.99，RMSE为0.73%，树高的一致性稍差，R2为0.66，RMSE为0.61 m（图6b）。总体来看，SGM模型能较好模拟样地各向异性反射特征。在反射率、结构参数等均在适当精度的情况下，获取的c1～c4参数值描述的Walthall背景反射模型参数也视为能正确反映场景背景反射的各向异性特征。

3.2 Walthall模型参数敏感性分析

Walthall模型的4个参数中，c1描述背景反射的基本强度水平、c2描述前后向反射间的差异程度、c3描述天底角反射强度、c4描述天底角反射与前后向反射的差异程度。为评价Walthall模型参数对最终反演结果的影响，使用无人机获取的样地结构参数反演得到的c1～c4值为参考，分别对其中1个参数值进行±20%范围变化，同时固定其他3个参数，反演得到覆盖度、树高、冠幅和林分密度，与无人机测得的结构参数进行平均相对误差（mean relative error，MRE）计算比较结果见图7。

可以观察到，斜率越大的参数对变量影响越大。c1值偏离参照值±10%，覆盖度误差增加20%～40%；偏离参照值±20%，覆盖度误差增加60%～80%；其次是c4，c2和c3影响最小（图7a）。由于模型中树高由冠幅计算得到，c1～c4变化对树高和冠幅的影响相同。相对覆盖度而言，树高、冠幅对c1～c4在正负变化方向上的响应比较复杂（图7b、c）。总的来说，影响程度最大的还是c1，其次是c4，c2和c3影响最小。c1～c4对林分密度值变化影响相对最小（图7d）。实际结构参数变化受到c1～c4变化的叠加影响。

3.3 Walthall模型参数回归及精度分析

通过建立训练集AN相机的多光谱（Blue、Green、Nir）反射率、AMBRALS模型参数fiso、fvol、fgeo和c1～c4参数的多元线性回归方程，用于对测试集样地的背景反射进行预测。

使用多元线性回归方程预测的训练集样地c1～c4参数，与反演值进行对比，c1～c4参数的调整R2分别为0.78、0.98、0.59和0.75（图8）。

图8 Walthall参数回归建模Fig. 8 Regression of Walthall model parametersa. c1；b. c2；c. c3；d.c4.

3.4 测试集森林结构参数反演及精度分析

输入17个测试集AN相机的多光谱（Blue、Green、Nir）反射率、AMBRALS模型参数fiso、fvol、fgeo，使用Walthall模型参数多元线性回归方程，计算出测试集样地的背景反射，作为SGM模型的输入信息，其他输入数据还包括模型固定参数LAI、HB、BR值。

以测试集观测反射值与模拟反射值的RMSE最小化为目标函数（式7），使用GRG优化算法计算输出17个测试集样地的树高、覆盖度、林分密度和冠幅等参数。

17个测试集样地反演结果与无人机测量获取的数值比较如图9。

图9 测试集结构参数评价Fig. 9 Evaluation of structural parameters of test dataseta. 覆盖度FVC；b. 树高Height；c. 林分密度Density；d. 冠幅Crown width.

在测试样地（275 m2）尺度上，与无人机获取的结构参数相比，反演获取的测试集结构参数使用线性模型拟合时，覆盖度、树高、林分密度、冠幅的R2分别为0.54、0.47、0.41、0.24，RMSE分别为3%、0.76 m、112株、0.64 m，MRE分别为24.7%、8.9%、22.5%、10%；使用幂函数模型时，覆盖度、树高、林分密度、冠幅的R2分别为0.80、0.53、0.55、0.30，RMSE分别为1%、0.4 m、33 株、0.32 m，MRE分别为10%、5%、6%、6%（图9）。

4 结论与讨论

基于MISR多角度遥感观测数据和SGM模型，通过模拟塔里木河下游河岸胡杨林的二向反射特性，分解后获得地面背景反射特征，最终反演得到样地主要森林结构信息。

研究得到主要结论如下：

1） UAV获取树高、冠幅、林分密度与实测值的R2分别为0.90、0.84、0.94，RMSE为0.45 m、0.68 m、4.25株·hm-2，无人机获取的胡杨林结构数据可用作SGM模型的定标参数以及最终精度验证数据。

2） SGM对训练集样地的模拟反射值R2最大值为0.99、最小值为0.72、均值为0.92，模拟覆盖度R2为 0.99，RMSE为0.73%，树高R2为0.66，RMSE为0.61m。SGM模型能较好模拟胡杨林样地的各向异性反射特征。

3） 对Walthall模型参数进行敏感性分析发现，c1参数对覆盖度和树高的影响最大，其次是c4参数，c2和c3参数影响最小，c1～c4参数对林分密度值变化影响很小。

4） 多元线性回归方程预测的c1～c4参数的调整R2分别为0.78、0.98、0.59和0.75，说明使用AN多光谱反射率和AMBRALS模型参数可以较好描述Walthall模型参数。

5） SGM模型反演获取的Walthall模型参数精度以及由此建立的参数回归方程精度对于最终的结构参数反演结果有较大影响。

6） 使用多角度卫星遥感（MISR）反演获取的胡杨覆盖度、树高、林分密度和冠幅的线性模型R2分别为0.54、0.47、0.41和0.24，幂函数模型R2分别为0.80、0.53、0.55和0.30，说明使用多角度卫星遥感（MISR）和反演方法可以实现区域尺度上的森林调查。

研究中存在的问题：

1） SGM反演需要获取与卫星影像像元尺度对应的森林结构数据作为模型定标数据，这需使用地面实测或无人机LiDAR、无人机摄影测量等技术获取高精度的森林调查数据。由于MISR像元与森林调查数据存在空间定位误差，加之分辨率相差较大，造成二者之间空间配准误差难以确定。干旱区环境下，河岸林植被覆盖度本身较低且空间分布异质性较大，较小的空间误差可能造成对应结构参数的大幅度改变，也会间接造成背景反射模型误差增大，最终影响到反演精度。对于这一问题，目前尚无较好的解决办法。

2） 通过对Walthall模型参数的敏感性和回归方程效果分析，发现背景反射模拟精度对反演结果有直接影响。研究区较为复杂的背景反射（高亮土壤、暗色河流和冠层反射与胡杨接近的灌木等）均增大背景反射模拟的难度和不确定性。如何获得更准确的Walthall模型标定参数以及在此基础上建立更有效精确的回归模型，是提高最终反演结构数据精度的关键。目前情况下，有限的定标站点不能代表所有实际场景，因此，基于样地提取的Walthall参数代表性也受到限制，在此基础上产生的参数回归方程会进一步增大该误差。一个可能的改进办法是采用DART（discrete anisotropic radiative transfer model）等仿真模型，根据场景组分实测光谱，使用不同组分组合方案，重建MISR尺度大小的胡杨林场景。模拟该尺度场景下中等分辨率光学卫星（如Landsat、Sentinel-2）多光谱反射( Mortonet al., 2015; Gastellu-Etchegorryet al., 2015），建立与背景反射模型参数的回归关系，生成Walthall模型参数的参数矩阵查找表（look up table, LUT），然后使用真实Landsat、Sentinel-2多波段反射率回归生成背景反射，提高对背景的模拟精度。

3）SGM模型反演结构参数的一个重要部分是优化算法，本研究使用的GRG是处理非线性约束问题的一种有效算法，但是 GRG 无法保证最后的寻优结果是全局最优值。当针对多极值问题时，其最终校准结果是一个临近初始给定值的局部最优值（李川等，2014）。与常规优化演算法相比，遗传演算法等近些年来出现的一些模仿自然选择与进化的随机搜索演算法只要能够提供合适的突变概率和运算时间，就能够确保其很好地遍历整个参数空间，最后寻优结果也往往能够落在一个较优的局部参数区域内。因此，使用更好的优化算法是提高反演精度的可行办法之一。