基于无人机影像的天山云杉树冠截留雪水当量估算分析∗

米尔扎提·依明，比拉力·依明，刘素红，伊加提·亚尔买买提，玛丽亚木·玛木提

(1.新疆大学 资源与环境科学学院 绿洲生态教育部重点实验室，新疆 乌鲁木齐 830017；2.北京师范大学 地理学部，北京 100875)

0 引言

雪是地球表面最为活跃的自然要素之一，雪面积、雪深、雪的反射率等特征是全球能量平衡模型中的重要输入参数[1].地球陆地表面的60%以上被季节性积雪所覆盖[2]，雪在水循环过程中占有重要地位，是陆地水资源尤其是干旱区水资源的重要组成部分，影响着区域甚至是全球尺度的气候变化过程，对区域流域的水文和生态具有极大的贡献意义[3−4].

天山是我国西北河流的重要发源地之一，天山云杉（Picea schrenkiana）森林在天山水文生态过程中起着非常关键的作用，但是冬季天山云杉的树冠截留造成了很大的水分流失.森林冠层截留降雪的再分配作用深刻影响森林内积雪的空间分布.在北方针叶林区域，树冠对降雪的截留比例很大，日本北部针叶林截留雪盖损失26%[5]，芬兰东部和北部雪盖截留损失约11%～30%[6]，北美山地针叶林的截留损失高达60%，雪水当量达到40 mm[7].森林冠层和大气之间强烈的能量交换对森林积雪、冠层截留以及融化等过程有显著的影响，季节性积雪覆盖区域是水文循环的重要过程之一[8−10].深度开展森林冠层对降雪截留研究，对于西北干旱区森林生态系统保护和流域水资源优化管理具有重要的意义.

基于遥感技术的雪盖面积检测起步较早并较成熟，被动微波遥感用于积雪监测始于1978年搭载多通道扫描微波辐射计（SMMR）的雨云-7卫星发射[11].随着MODIS、Landsat、风云（FY-3A）、环境减灾（HJ-1B）气象卫星等更多可见光近红外传感器卫星的发射，基于遥感的积雪监测业务能力得到了大幅度提高和发展.黄镇等利用MODIS数据以新疆为研究区，初步探讨了利用归一化雪盖指数NDSI识别积雪时的阈值问题[12].曹云刚和刘闯试图对比TM影像和MODIS影像，寻找NDVI、NDSI、雪盖比例三者之间的曲线关系，并制作了青藏高原的雪盖图[13].基于卫星遥感影像对积雪盖度监测建立在下垫面均一的理想平面上，但实际遥感监测中因下垫面各异，其积雪变化过程及积雪深度存在差异[14−15].近年来，国内外相关学者在雪水当量的微波遥感方面做了大量工作，如用被动微波遥感[16−18]、合成孔径雷达[19]进行了多种算法的反演[20]，但被动微波遥感影像的空间分辨率普遍较低，以至于难以准确获取雪水当量信息[11].近年来无人机遥感技术高速发展，其分辨率高、数据获取简单、低成本、时效性高、低风险等优势凸显[21−22].已广泛应用于土地利用分类、农业资源调查[23−24]、林业资源调查、林业有害生物防治、森林防火等方面[25−26]，正逐渐成为传统航天遥感及卫星遥感的优异补充，因此需要在卫星遥感手段的基础上结合无人机低空遥感数据进行补充和验证，从而提高定量遥感反演雪水当量精度.

本研究使用无人机对位于乌鲁木齐县鹰沟风景区的一片试验区进行成像并找出多个不同角度的单株云杉，对其进行无人机和云杉之间的角度计算、云杉面积计算，并得出单株云杉降雪截留面积以及不同角度下的变化.结合人工实验和无人机低空遥感数据，利用本团队通过模拟降雪的方式得出的云杉树苗雪盖面积和树冠截留雪水当量二者之间的关系[27]，计算得出单株云杉树冠截留雪水当量.为基于无人机研究森林积雪、天山云杉冠层对降雪截留和雪水当量相关影响提供参考和指导.

1 研究对象与研究方法

1.1 研究区概况

乌鲁木齐县地处亚欧大陆腹地，中国西北部，新疆中部，坐落于天山北麓，准噶尔盆地南部，介于86°37′33′′～88°58′24′′E，42°45′32′′～45°00′00′′N（见图1）.地形以山区、盆地和平原为主，地势南高北低，高差达4 000 m，流域内植被和土壤的垂直地带性显著.该流域属中温带大陆性干旱气候，全年干旱少雨，年降水量500～600 mm，降雪量32.4 mm，约占年降水量的12.3%[28].冬夏两季气温变化剧烈，年均气温3°C.鹰沟风景区位于乌鲁木齐县中部，天山云杉是林区主要的乔木，平均高度达15 m，是天山云杉分布比较典型的区域.

图1 试验区位置图

1.2 影像数据获取与处理

实验数据为无人机低空遥感影像，无人机为大疆精灵4 Pro V2.0，搭载CMOS相机，有效像素为2 000万，成像时间为2018年11月8日，无人机飞行稳定，高度为50 m，设置横向和纵向的飞行航线，各平行航线距离为10 m，整个拍摄过程中镜头垂直向下，进行多次航拍.一次飞行拍摄影像约80张.从无人机影像中选取树冠信息完整且未被其它云杉遮盖的、多个不同角度的单株云杉影像.为了提高分类的精度，借助图像处理软件对筛选出来的三组无人机影像进行裁剪，去除云杉树冠阴影、地面雪等背景信息.

本实验利用最大似然法对无人机影像提取树冠雪盖面积.采用最大似然分类方法进行监督分类，对类别样本对应的影像数据进行分类处理，记录影像数据隶属各类别的数量.根据训练样本提取特征向量，构造出各类别的多维正态分布模型，即概率密度函数或者概率分布函数.在得到各类的多维分布模型后，对未知类别的数据向量，通过贝叶斯公式计算概率值，属于哪一类的概率大（最大值选择器），则将该数据向量或者像元归为此类[29].

不考虑先验概率的最大似然分类法的判别函数为[30]:

式中：ωi代表类别i；g1i(X)代表ωi的判别函数；X为待分类的特征向量；P(x/ωi)为ωi的概率密度函数.

考虑先验概率后，最大似然分类法的判别函数为[31]:

式中：P(ωi)为ωi的先验概率，即ωi类在影像中的面积比例，当假定各类别具有相等的分布面积时，式(2)等同于式(1) ；g2i(X)表征着特征向量X对ωi类的后验概率P(x/ωi)的大小.因此，式(2)按照后验概率的大小确定像元的归属，不仅考虑了特征向量与预先给定的类别特征的相似程度，也考虑了类别在影像中的面积比例.

1.3 研究方法

1.3.1 无人机与云杉实际距离计算

本文通过最大似然法对单株云杉图像进行分类后，计算出每个分类特征的像元总数和占总像元数的百分比.无人机影像的大小为193.04 cm×128.69 cm，通过影像面积和各类特征的像元百分比乘积，得到各类特征在无人机图像上所占比例.为得到单株云杉的实际面积，在影像上选择一个已知尺寸的参考地物，以无人机成像时人工测量得到的两柱电线杆间的距离为参考，推算图像的比例尺信息.人工测量得到的两柱电线杆实际距离为15 m，在影像上的距离为68.85 cm，按此比例可以推算出影像的实地长和宽（长42.05 m，宽28.04 m）.最终计算出三组单株云杉的各类特征的实际面积.

1.3.2 单株云杉表面积及各类特征面积占比估算

假设单株云杉为圆锥体，计算圆锥的表面积需要底部圆的周长和侧面三角形的高.通过单株云杉成像距离和不同距离的云杉投影面积的相关性，构建二者之间的拟合方程.如图2所示，无人机拍摄过程中，当处于云杉正上方时Ds=0，通过拟合方程可估算云杉底部圆的面积；当与云杉达到一定距离后，云杉的侧面积达到最高值.把侧面积的最高值当作云杉侧面三角形的面积，借助云杉底面积信息推算云杉底部直径，用侧面积和底部直径来推算三角形的高.通过圆锥表面积公式S云杉表面积=1/2hC（S为面积，C为云杉底部周长，h为云杉侧面三角形的高）计算三组单株云杉的表面积.

图2 无人机成像原理图

1.3.3 云杉树冠截留雪水当量估算

雪水当量是指当积雪完全融化后所得到的水量.以本研究团队模拟降雪实验所构建的云杉雪盖面积与树冠截留的雪水当量模型为基础，估算了三组云杉树冠截留雪水当量[27].公式如下：

式中：W为雪水当量（cm），S为雪盖面积（m2）.

2 结果与分析

2.1 基于不同角度、不同距离的降雪截留面积与特征

分别对三组单株云杉进行不同角度无人机成像（共29张），并对第一组和第三组的10个不同角度、第二组的9个不同角度的单株云杉图片进行最大似然法分类，分类结果和位置示意见图3，黄点1～10为无人机拍摄的位置.分类结果基于像素的降雪截留面积，为得到实际面积，借助地面参照物得到了云杉实际面积和无人机成像的距离.

图3 三组单株云杉分类和位置示意图

三组云杉东西逆时针方向的不同角度信息对单株雪盖面积和不同成像角度下的变化趋势、雪盖面积变化趋势见图4.可以看出，三组云杉雪盖面积在不同角度下的变化趋势基本一致.由于一天中太阳方位角度变化，云杉向北面光照时间和光照强度比向南面大，所以出现云杉向北面和向南面的雪盖面积大小差异，向北面从东向西方向先减小后增大，三组单株云杉雪盖面积最小值出现在90°附近，面积分别为1.59 m2、0.68 m2和1.41 m2.三组云杉向南面雪盖面积变化幅度较小，雪盖面积最大值出现在270°～320°之间，分别为4.62 m2、7.81 m2和5.47 m2.单株云杉在不同角度的雪盖面积曲线中可能因植株的不同，雪盖面积变化幅度和变化曲线的曲率略有不同，但总体上仍保持基本一致.

图4 角度与雪盖面积变化曲线

2.2 云杉树冠雪盖面积及树冠截留雪水当量

图5为第一组云杉的成像距离和不同距离投影面积拟合曲线，拟合回归方程为y=-0.012 82x2+0.513 51x+2.939 93，拟合曲线的判定系数为0.948 6，有显著的相关性.通过拟合方程得出，当云杉处于无人机正下方（Ds=0）时，云杉冠幅面积为2.94 m2，利用圆锥表面积公式得出第一组云杉表面积为24.69 m2.

图5 成像距离和投影面积拟合曲线

对云杉的成像距离和云杉侧面积进行拟合，得出第二、第三组云杉的表面积分别为36.67 m2和27.05 m2，并通过三组云杉影像的每种特征平均值来推算雪（光照）、非雪（光照）、非雪（阴影）、雪（阴影）四个分类特征在云杉表面积中所占的特征实际面积（见表1），三组不同云杉降雪截留面积占总表面积的比例分别为42.5%、50.5%和79%.利用公式（3）得出三组云杉树冠截留雪水当量分别为2.324 L、7.267 L和9.683 L.

表1 三组云杉表面积所占的特征实际面积和雪水当量信息/m2

3 讨论

天山云杉树冠截留在降雨方面研究比较充分[32]，而对树冠截留降雪分配的研究目前较少.从观测尺度上来看，对森林树冠截留的研究主要集中在两个尺度：单株/小区尺度和卫星遥感尺度.Dozier[33]利用TM影像的2、5波段建立归一化积雪指数（NDSI）后，Harrison等[34]利用NDSI指数进行30 m/1 000 m尺度的雪盖面积制图研究，并发展了二值分类算法、混合像元分解算法、监督和非监督分类算法等，快速获取区域雪盖面积.但卫星遥感负责较大面积的积雪覆盖检测，很多区域研究的模型模拟分辨率较低，难以分出地面和树冠上截留的雪，需要地表验证，而利用无人机低空遥感技术以较低的成本可以获取高分辨率低空遥感数据.基于此，本实验利用最大似然法，对不同角度下的无人机影像进行云杉降雪截留面积提取，精度达到95.7%.以本研究团队模拟降雪实验所构建的云杉雪盖面积与树冠截留的雪水当量模型为基础，估算了三组云杉树冠截留的雪水当量，达到较好的提取效果.

云杉树冠截留雪与云杉冠层的承载能力及云杉的叶面积指数在很大程度上有关联，冠层截留雪取决于冠层截留雪和枝条的结合度、冠层截留雪晶体间的黏结强度以及枝条的支撑强度[35].森林降雪截留的影响因素很多，影响树冠截留雪盖蒸发和升华损失的因素包括林冠特征、空气温度[36−37]、风速[9]、湿度和地形[38−40]和太阳辐射等，这些因素相互耦合，共同作用于融雪过程[41−42].因此，天山云杉降雪截留面积的研究需要把握好时间因素，从降雪停止到开始收集数据之间的时间间隔会影响树冠降雪截留面积的大小.本实验收集数据的时间是降雪停止后的第二天，树冠上部分积雪已经蒸发或者滴落.

无人机影像估算云杉降雪截留研究中，因受太阳高度角和山体遮挡影响，在云杉林密度大的情况下，出现云杉冠层遮挡以及重叠，外围的云杉遮挡了太阳光照，从而会影响到高郁闭度区域的云杉降雪截留面积精度，由此可见，对于相互遮挡的高聚集分布林分仍需拓展思路，进一步开展研究.

4 结论

（1）利用最大似然法对天山云杉树冠截留的雪进行分类提取，提取精度达到95.7%.

（2）通过单株云杉多角度雪盖面积变化趋势可知，向北雪盖面积呈现出先减小后增大的趋势，最小面积出现在90°附近，最大面积出现在270°～320°之间；向南雪盖面积变化幅度比向北变化幅度小.由于太阳方位角的变化，云杉向北光照时间长，从而导致正北方位雪盖面积最小.云杉向南受太阳照射影响小，融化速度比较慢，雪盖面积大于向北.

（3）通过成像距离和不同距离的云杉投影面积拟合曲线，借助拟合方程和几何数学原理计算出了三组云杉表面积，分别为24.69 m2、36.67 m2和27.05 m2，雪盖面积在表面积中所占的比例分别为42.5%、50.5%和79%.

（4）将人工实验结果和无人机低空遥感数据相结合，利用本研究团队通过模拟降雪方式得出的云杉树苗雪盖面积和树冠截留的雪水当量二者之间的关系[27]，计算得出单株云杉树冠截留的雪水当量分别为2.324 L、7.267 L和9.683 L.