基于无人机影像天山云杉林主伐迹地提取研究

王雅佩，王振锡，刘梦婷，李 擎，杨勇强

(新疆农业大学林学与园艺学院，新疆教育厅干旱区林业生态与产业技术重点实验室，乌鲁木齐 830052)

0 引 言

【研究意义】天山云杉(PiceaSchrenkianavartianshanica)是新疆山地森林的主要优势针叶树种，是北方森林重要的组成部分，天山云杉林在水源涵养、水土保持以及维持新疆绿洲生态平衡等方面发挥着不可替代的重要生态功能[1-2]。天然林资源及其更新对森林的演变及发展具有重要的作用，是长期以来森林生态系统主要研究领域之一，这个过程中其受环境条件、自然干扰和人为干扰类型等因素及其相互作用的影响[3-4]。在林业生产中，森林主伐及其更新对森林资源的有效利用和恢复有着不可分割的作用，并且与森林的结构、功能和动态之间存在着密切的关系[5]。该方法耗费大量的人力物力财力、周期长、时效性低，特别是对主伐迹地而言，采伐地点分散并且不明确，致使野外调查时间延长、工作量增加，不足以满足当代林业发展的需要[6-7]。遥感技术是森林资源高效快速监测调查的重要技术手段。无人机遥感作为新兴遥感技术，能够在短时间内获取所需数据，通过自主设置航线更加灵活获取遥感数据，具有环境适应力强、分辨率高、低成本等特点，使无人机遥感技术成为获取高分辨率遥感数据的重要技术手段之一[8]。开展森林资源调查和统计对及时获取林木的生长状况，和森林未保育有实际意义。【前人研究进展】利用无人机影像对地物信息提取的应用研究已越来越广泛，由于其高分辨率的特性，使其分类方法主要是以面向对象方法为主。Chenari A[9]基于无人机正射影像，使用面向对象分类方法提取研究区的野生开心果和野生杏仁，并估算平均树冠面积，结果表明该方法对野生开心果估计精度为0.92，野生杏仁估计精度为0.89。吴金胜等[10]基于无人机影像，采用面向对象对估算研究区水稻面积，结果表明总体分类精度达到93%以上，该方法能够替代目视解译法，有效提高了遥感面积估算效率。彭佳忆等[11]基于无人机可见光影像采用面向对象方法提取荒漠地表类型信息，结果表明面向对象方法精度优于最大似然法，面向对象多尺度分割方法提取荒漠地表类型信息具有高精度。【本研究切入点】利用中高分辨率遥感影像在采伐监测中的研究有很多，但基于无人机可见光影像对山区主伐迹地提取及其分布研究方面还未见报道。研究基于无人机可见光影像，采用面向对象多尺度分割方法，利用最邻近分类方法与SEaTH算法相结合的分类方法提取不同等级的主伐迹地的面积。【拟解决的关键问题】基于无人机可将光影像，利用eCognition 9.0软件，采用面向对象多尺度分割技术，确定各层次最优分割尺度，通过最邻近分类以及SEaTH算法相结合的分类方法筛选分类特征，并建立类别层次结构及分类规则集，以实现提取不同等级主伐迹地面积，为新疆天保工程实施后山区森林资源更新恢复及评价提供科学依据。

1 材料与方法

1.1 材 料

1.1.1 研究区概况

研究区位于天山山脉中段中山带的新疆农业大学实习林场内，区域地处天山北麓乌鲁木齐西南方向110 km处，位于头屯河上游，地理坐标范围为43°16′～43°26′N，86°46′～86°57′E。林场年平均降水量约600 mm，其海拔高度在1 700～3 200 m间。经营面积为1.044 4×104hm2，其中林业用地面积为5.116×103hm2。地形多样，南部地势高，北部地势低，地形切割较为剧烈，坡度主要在10°～40°间。土壤发育较好，林下土壤为普通灰褐色森林土。林区以天山云杉纯林为主。

1.2.1 影像数据

采用的无人机搭载的相机型号是SONY DSC-RX1RM2(35 mm)相机，搭载镜头焦距为35 mm。为了避免夏季背景植被(草本、灌木等)的影响，无人机数据拍摄于2017年11月15日上午，天气晴朗，能见度高，风速较小。无人机起飞前，预先设定好飞行的方向、航线和高度载入飞行控制系统，此次飞行路线的航向重叠为80%，旁向重叠为60%，飞行高度为923 m。获取的无人机影像地面分辨率为0.118 m。

1.2.2 调查数据

数据包括新疆农业大学实习林场二类调查数据，样地调查数据以及林场伐区设计资料等其他辅助资料。2017年8月，在新疆农业大学实习林场内进行实地调查。根据天然林保护工程实施前的伐区设计材料确定主伐迹地，在主伐迹地范围内选取典型样地，设置30 m×30 m的典型样方，对样地内天山云杉进行每木检尺，实测并记录样方内全部天山云杉的树高、胸径、年龄、冠幅、株数等信息以及样地的郁闭度。利用亚米级的集思宝GPS记录样方的中心点及4个角点的位置信息，以便于将实测样地与无人机影像进行一一对应。

1.2 方 法

1.2.1 无人机影像处理

采用PIX4D软件对无人机数据进行处理，主要过程包括：(1)首先筛掉航拍中成像效果差、模糊等不合格的照片；(2)建立测区，导入有效的照片、POS数据和控制点文件，设置参数后提交，软件会自动提取照片及相机的相关信息；(3)数据的自动拼接，调整拼接线、投影切换和混合影像，最终生成正射影像(DOM)。

1.2.2 多尺度分割

在eCognition Developer 9.0软件，采用面向对象的方法提取不同等级的主伐迹地。面向对象信息提取中影像最优分割尺度的确定是关键技术之一。采用ESP算法[12]对最优分割尺度进行计算。ESP算法完成后，得到横坐标为尺度参数，纵坐标分别为局部方差(local variance，LV)和局部方差变化率(rate of change，ROC)的折线图。当ROC呈波峰状态即出现局部最大值时，则表示此尺度为某种地物对象的最佳分割尺度。

在多尺度分割技术中，影像分割结果的质量与均质性因子的设置密切相关。均质性因子包括颜色因子(Color)和形状因子(Shape)，其中形状因子又包含平滑度(Smoothness)以及紧致度因子(Compactness)。其中颜色与形状因子值是相对的，权值和为1；平滑度与紧致度因子权值和为1。颜色因子提供了对象的光谱信息，形状因子反映了分割对象的大小，平滑度因子是优化分割对象边界的平滑程度，紧致度因子使用形状标准来考虑整体紧致性，并优化影像对象结果的紧凑度。

1.2.3 分类

最邻近分类算法是一种监督分类方法，是通过选择训练样本，构建特征空间，计算相关特征之间的距离，根据距离确定对象所属类别[13]，进行分类。如果距离值越小，则表示该对象越可能隶属于该样本。为了避免由于采用过多参与分类的特征所带来的计算量大、特征冗余和分类低等问题，eCognition软件提供了特征空间优化的模块(Feature Space Optimazation)，选取预先选择的特征，通过该模块计算得到特征优选的结果，可以确定在几维时类别与类别之间的可分性最大，特征优选的结果可以直接应用于分类体系中执行分类。计算公式如下：

(1)

表1 初选特征因子
Table 1 Primary selection factor

特征Factor特征因子Texturefactor特征Factor特征因子Texturefactor纹理特征Texturefactor均值(Mean)方差(Variance)熵(Entropy)对比度(Contrast)均值性(Homogeneity)相异性(Dissimilarity)二阶距(SecondMoment)相关性(Correlation)光谱特征Spectralfactor形状因子Shapefactor平均值(Mean)亮度(Brightness)标准差(StdDev)波段比(Ratio)模拟植被指数(VI')长宽比(Length/Width)形状指数(ShapeIndex)

其中，VDVI表示可见光波段差异植被指数，VI' 表示模拟归一化植被指数提取植被的一个参数，公式如下：

(2)

VI'=(2×G'-R'-B')-(1.4×R'-G') .

(3)

SEaTH(Seperability and Thresholds)算法[14]是基于各类别的样本特征值，利用J-M(Jeffries-Matusita)距离实现两两类别间最优特征的自动选择，及其最优特征的阈值确定。SEaTH算法主要包括两部分：

(1)对区分两两类别特征进行优选

首先，选择各个类别的影像对象代表性样本，选择预先考虑到的特征，并导出样本的特征值。SEaTH算法采用J-M距离计算两两类别之间在某个特征的分离度，以评估两个类别的可分离性。计算J-M距离公式如下：

J=2(1-e-B) .

(4)

(5)

(2)最优特征阈值的确定

该算法根据高斯概率分布公式计算两两类别在某一个特征的最佳阈值。其公式如下：

p(x)=p(x|C1)p(C1)+p(x|C2)p(C2)

(6)

当p(x|C1)=p(x|C2)时，这两种类别间的混分最少。C1、C2表示两个类别，当特征值为X1时，C1和C2这两个类别的分离性最好，此时X1即为最佳特征阈值T，计算方式如下：

(7)

(8)

式中，T为该特征的阈值，T介于m1和m2之间，m1和m2分别为两个类别在某个特征上的均值，σ1σ2和σ2σ2分别为两个类别在某个特征上的标准差，n1和n2为类别C1和C2的样本数量。图1

图1 C1和C2最佳分离特征阈值
Fig.1 The schematic diagram of separation threshold between C1 and C2

1.2.4 分类精度评价

研究选取生产者精度、用户精度、Kappa系数和总体精度四个指标对不同等级主伐迹地提取的分类结果进行精度评价。

2 结果与分析

2.1 主伐迹地等级划分

根据实地主伐迹地调查数据包括林分平均胸径、平均树高、年龄及林分密度均和目视解译法勾绘获得的样地郁闭度及平均冠幅，利用R软件，采用中间距离法进行聚类分析，将13块典型样地划分为三个等级。主伐迹地天山云杉林的林分密度随着主伐迹地等级的增加逐渐提高。主伐迹地林木相较于周围林木树冠较小，且数量相对较多，且林木排列较为规律伴随着周围有天然更新的林木。不同等级的主伐迹地在图像上表现为等级Ⅰ主伐迹地相对林木树冠小，林木株数少以及裸地面积大；等级Ⅱ主伐迹地较等级Ⅰ的树冠大小适中，林木株数较多；等级Ⅲ主伐迹地树冠林木树冠较大，林木株数数量多以及裸地面积较小，且林木大小不一。表2,图2,图3

表2 主伐迹地分级及指标
Table 2 Classification and indicators of main cutting areas

等级Grade平均胸径AverageDBH(cm)平均树高Averageheight(m)平均年龄Averageage(a)平均冠幅Averagecrownwidth(m)郁闭度Crowndensity林分密度(株/hm2)Standdensity(plant/hm2)Ⅰ5.33～5.883.66～4.3615～221.68～1.860.20～0.54833～2100Ⅱ1.81～6.591.78～5.7517～261.01～2.480.28～0.872267～3367Ⅲ6.67～9.566.93～9.7326～291.13～1.990.41～0.653833～4632

图2 聚类结果
Fig.2 Clustering result graph

注:(a)(b)为Ⅰ等级，(c)(d)为Ⅱ等级，(e)(f)为Ⅲ等级

Annotation：(a) (b) is grade I, (c) (d) is grade II, (e) (f) is grade III

图3 不同等级主伐迹地部分样地示意
Fig.3 Images of parts of main cutting sites of different grades

2.2 最优分割尺度的确定

在eCognition软件支持下，采用ESP算法对研究区无人机影像数据进行最佳分割尺度计算。一次分割主要目的是为了区分植被和非植被，本次ESP处理计算结果中最佳分割尺度是43、82，113、188、245、274、316、343、404、467和504，使用这些最佳分割尺度分别对影像进行分割，通过观察其视觉效果，当分割尺度为43时，不论是道路、地面、建筑物还是林地等地物都被分割的十分破碎，未能达到此次分割的目的；当分割尺度为113时，发现不同地物均有合并现象，但分割效果仍然十分破碎；当分割尺度在113至245之间时，发现地面以及林地地物逐渐被合并，而道路和建筑物地物合并效果并不明显；当分割尺度在245之后开始，可以发现道路地物逐渐被合并，林地地物中也有部分区域进行进一步的合并；当分割尺度为343时，道路、地面、建筑物和林地地物中基本完成了对象的合并，此时分割的效果较为接近第一次分割的目的；当分割尺度为504时，道路地物的线状信息表达较为清晰，地面以及林地地物中的对象大面积的合并，几乎不存在破碎现象，达到此次分割目的，因此将分割尺度为504作为第一次分割的最优分割尺度，通过不同形状因子以及紧凑度因子的测试，最终确定形状因子为0.4，紧致度为0.5。

在eCognition软件支持下，采用ESP算法进行二次分割，其目的为在植被类别中提取主伐迹地，经过对比最后确定提取主伐迹地的最优分割尺度为292。对形状因子以及紧致度因子参数进行调节，最终确定形状因子以及紧致度因子分别为0.5、0.3。图4、图5,表3

图 4 一次分割ESP最佳分割尺度
Fig. 4 The optimal segmentation scale line chart in the first Segmentation 表3 多层次分割参数Table 3 Scale parameters of the multi-level segmentation

层次等级Hierarchicallevel目标地类Targetlocationclass背景地类Backgroundgroundclass分割尺度Divisionscale形状因子Shapefactor紧致度因子CompactnessfactorLevel1植被Vegetation非植被Non-vegetation5040.40.5Level2主伐迹地Finalfellingarea非主伐迹地Notfinalfellingarea2920.50.3

图5 一次分割不同分割尺度对比
Fig.5 The comparison chart of different segmentation scales in the first Segmentation

2.3 面向对象主伐迹地提取

2.3.1 level 1分类

在level 1中，主要区分植被与非植被，提取目标地类为植被。由于研究所使用的影像为无人机可见光影像，通过查阅文献最终确定VI' 和VDVI两个植被提取指标，经过多次试验对比观察，最终以VI' 大于等于-0.14将植被与非植被区分开。图6

图6 level 1层分类结果
Fig.6 The classification result of level 1

2.3.2 level 2分类

在level 2中，对level 1植被类别进行继承，在对植被类别二次分割前对一次分割结果进行合并，主要是在植被类别中提取主伐迹地，由于主伐迹地是2000年以前进行主伐过后的，相对于其他未主伐过的林地，林木相对于较整齐且树冠相对于较小，且单元面积内的林木株数相对于较多。通过SEaTH算法，计算主伐迹地和非主伐迹地类别间各种特征的J-M距离及其阈值的确定，选择J-M距离最大的前5个特征进行分类。表4、图7

图7 level 2层分类结果
Fig.7 The classification result of level 2表4 目标地类与背景地类特征J-M距离值及其阈值Table 4 The J-M and threshold value of features between target land covers and background land covers

目标地类Targetlocationclass背景地类Backgroundgroundclass特征FactorJ-M距离J-Mdistance阈值Threshold主伐迹地Finalfellingarea非主伐迹地NotfinalfellingareaHIS颜色特征亮度(Intensity,I)0.97≤0.497R波段均值(Mean_R)0.97≥88.5R波段波段比(Ratio_R)0.93≤0.167灰度共生矩阵均值(GLCM_Mean)0.83≤122.4亮度(Brightness)0.83≤112.79

2.3.3 Ⅰ、Ⅱ、Ⅲ等级的主伐迹地提取

在level 2层次提取的主伐迹地类别中，由于SEaTH算法适用于两两类别间的区分，而此次是为了提取Ⅰ、Ⅱ、Ⅲ等级的主伐迹地，根据主伐迹地等级，采用之前筛选的遥感特征因子与SEATH算法计算每两个等级间的J-M距离最大的特征因子相结合的方式，最终确定提取各等级主伐迹地的特征因子，经过多次试验对比观察确定阈值，建立其分类规则，提取Ⅰ、Ⅱ、Ⅲ等级的主伐迹地。表5、图8

表5 不同等级的主伐迹地分类规则
Table 5 Classification rules for main cutting areas of different grades

目标地类Targetlocationclass特征Factor规则表达式RegularexpressionⅠ色饱和度(Saturation,S)≤0.7476亮度(Brightness)≥133.313G波段平均值(MeanG)≥102.1649ⅡR波段平均值(MeanR)≤157.2225灰度共生矩阵均值(GLCMMean)>141.6125R波段波段比(RatioR)≤0.2273ⅢHIS颜色特征亮度(Intensity,I)≥0.4595亮度(Brightness)≤117.1661灰度共生矩阵方差(GLCMVariance)>21.437

图8 不同等级的主伐迹地分类结果
Fig.8 Classification results of main cutting areas of different grades

2.4 Ⅰ、Ⅱ、Ⅲ等级的主伐迹地面积精度评价及其分布

所使用无人机影像分辨率较高，达到了0.118 m，因可根据目视解译结果对分类结果进行验证，并提取不同等级的主伐迹地面积与目视解译结果作对比。利用eCognition软件中提供的基于样本的误差矩阵的精度评价方法，通过选取目视解译与分类结果图像中对应位置的提取类别进行比较计算，采用混淆矩阵的统计方法对分类结果进行分析以及对比。

从生产者精度来看，主伐迹地Ⅰ、Ⅱ和Ⅲ精度分别为80.77%、83.33%和64.10%；从用户精度来看，主伐迹地Ⅰ、Ⅱ和Ⅲ精度分别为84.00%、80.81%和71.43%。从总体看，利用最邻近分类与SEaTH算法相结合的分类方法总体分类精度达到81.82%，Kappa系数为0.74，实验分类质量取得了较好的效果。表6

表6 分类精度评价结果
Table 6 Classification accuracy evaluation results

地类CategoryⅠⅡⅢ非主伐迹地NotfinalfellingareaⅠ42602Ⅱ48096Ⅲ06254非主伐迹地notfinalfellingarea64587列总和Columnsum52963999生产者精度ProducerAccuracy,PA80.77%83.33%64.10%87.88%用户精度UserAccuracy,UA84.00%80.81%71.43%85.29%总体精度OverallAccuracy,OA81.82%Kappa系数Kappacoefficient0.74

对于不同等级的主伐迹地面积的提取结果与目视解译结果对比，主伐迹地Ⅰ的面积相对误差为12.91%，面积吻合度为87.09%；主伐迹地Ⅱ的面积相对误差为20.14%，面积吻合度为79.86%；主伐迹地Ⅲ的面积相对误差为33.67%，面积吻合度为66.33%；提取主伐迹地总面积为72.574 3 hm2，目视解译面积为92.174 9 hm2，面积相对误差为21.26%，面积吻合度为78.74%。该方法能够提取研究区主伐迹地。表7

表7 不同等级主伐迹地面积精度
Table 7 Precision scale of main cutting area of different grades

等级Grade提取面积Extractionarea(hm2)目视解译面积Visualinterpretationarea(hm2)面积相对误差ArearelativeError(%)面积吻合度Areacoincidencedegree(%)Ⅰ11.086112.728912.9187.09Ⅱ51.906865.000620.1479.86Ⅲ9.581414.445333.6766.33合计Total72.574392.174921.2678.74

3 讨 论

在面向对象分类多尺度分割技术中，最优尺度直接决定影像对象的大小及信息提取的精度[15]。通常情况下最优尺度的选择需要反复多次对分割尺度进行人工调整，然后通过目视判断法确定，但运用EPS算法可有效减少尝试次数，能够准确、快速的确定最优尺度。

在分类特征优化方面，研究采用最邻近和SEaTH算法相结合的分类方法提取研究区内不同等级的主伐迹地，其总体分类精度达到了81.82%，取得了相对较好的分类效果，可有效避免了分类特征较多而导致的特征冗余、计算量大和分类精度不高等问题，并为最优特征及其阈值的确定提供了较为便捷、有效的途径。这与刘鑫[16]2017年基于无人机影像提取四川低丘区耕地信息的研究结论一致。

在分类精度方面，不同等级(Ⅰ、Ⅱ、Ⅲ)的主伐迹地均存在不同程度的漏分和错分情况，主要原因可能是不同等级(Ⅰ、Ⅱ、Ⅲ)主伐迹地之间的图像阈值具有一定的重叠度，同时本研究使用的无人机图像仅是可见光范围的地物成像，植被信息提取所需的光谱信息量不足。在今后的研究中可以尝试采用无人机多光谱图像，持续改进和更新分类算法开展研究，进一步提高分类精度。

从不同等级主伐迹地面积提取精度来看，主伐迹地Ⅰ提取面积吻合度最高，这是由于主伐迹地Ⅰ内的天山云杉林分密度低于主伐迹地Ⅱ和Ⅲ，其林分平均冠幅、平均年龄等因子也都相对较小，因此在影像上更好区分。

4 结 论

4.1 通过调查，研究区主伐迹地分为三个等级，HIS颜色特征亮度(Intensity，I)、R波段均值(Mean_R)、R波段波段比(Ratio_R)、灰度共生矩阵均值(GLCM_Mean)和亮度(Brightness)是主伐迹地提取的重要特征。色饱和度(Saturation，S)、亮度(Brightness)和G波段平均值(Mean G)是主伐迹地Ⅰ提取的重要特征；R波段平均值(Mean R)、灰度共生矩阵均值(GLCM Mean)和R波段波段比(Ratio R)是主伐迹地Ⅱ提取的重要特征；HIS颜色特征亮度(Intensity，I)、亮度(Brightness)和灰度共生矩阵方差(GLCM Variance)是主伐迹地Ⅲ提取的重要特征。

4.2 基于面向对象多尺度分割技术，采用最邻近以及SEaTH算法相结合的分类方法，借助eCognition软件获取最优特征空间，可有效避免“椒盐”现象，总体分类精度为81.82%，该方法在不同等级主伐迹地分类中具有较高的分类精度，相比仅用单一方法能够更便捷、更有效的建立类别层次结构及分类规则集，是一种行之有效的方法。