一种基于多特征优选的可见光无人机遥感影像森林树种分类
方法
技术领域
[0001] 本
发明涉及农林业资源遥感信息自动提取领域，特别是
一种基于多特征优选的可见光无人机遥感影像森林树种分类方法。
背景技术
[0002] 树种是森林资源信息中的一项重要参数，精准的树种信息提取对于森林资源调查、外来物种监测、生态健康评估等方面具有重要意义。传统的树种识别方法以实地调查为
基础，劳动强度大，耗时长，受空间范围的限制。遥感手段可对大范围内的目标地物进行周期性的可重复性的观测，目前已广泛应用于森林信息提取。其中，中等
分辨率的卫星遥感影像可对大区域范围的森林类型进行识别，但受限于分辨率，对于树种混交区域难以实现树种的分类识别，难以满足林分尺度树种资源调查的需求。高分辨率卫星遥感影像上不同树种类型的纹理差异明显，可进行林分尺度的树种分类识别，但是高分辨率卫星遥感数据获取难度大，受天气和回访周期的影响往往不能及时获取所需区域的影像。无人机遥感作为新兴的遥感技术，其时效性更强、成本低、数据易于获取，成为小区域林业遥感数据获取的新手段。
[0003] 目前利用无人机遥感数据进行树种分类研究主要使用的是多
光谱和高光谱数据，这两种类型的数据可以提供丰富的光谱特征，然而数据获取成本高处理流程复杂，实际林业生产和林业资源调查中使用最多的是仅包含红，绿，蓝三个波段可见光无人机遥感影像。可见光无人机遥感影像多用于林分株数、树高和蓄积量等参数的估算，在森林树种分类识别的研究还很少。行业和部
门积攒了大量的可见光无人机数据，没有充分发掘其应用价值。
可见光无人机遥感数据因其提供的光谱信息仅有三个波段，进行树种分类难度较大。特征的构建和特征选择是实现可见光无人机遥感数据树种分类识别的关键。
发明内容
[0004] 有鉴于此，本发明的目的是提供一种基于多特征优选的可见光无人机遥感影像森林树种分类方法，解决可见光无人机遥感影像光谱信息少，不能有效的区分多树种的问题。
[0005] 本发明采用以下方案实现：一种基于多特征优选的可见光无人机遥感影像森林树种分类方法，包括以下步骤：
[0006] 步骤S1：使用无人机航拍获取所需林区的可见光遥感影像，对可见光遥感影像进行预处理生成数字表面模型和数字正射影像图；
[0007] 步骤S2：建立可见光无人机遥感影像树种分类体系，并从可见光遥感影像上选择不同森林树种及地类样本数据；
[0008] 步骤S3：从步骤S1中生成的数字表面模型和数字正射影像图中提取多类型的特征构建特征空间；
[0009] 步骤S4：利用递归消除
随机森林特征选择
算法选择从步骤S3构建的特征空间中提取最佳分类特征子集；
[0010] 步骤S5：在在步骤S4获取的最佳分类特征子集上使用随机森林算法实现树种的分类，用以提取树种分布图，对分类结果进行
精度评价。
[0011] 进一步地，步骤S1中所述对可见光遥感影像进行预处理包括以下步骤：
[0012] 步骤SA：在无人机处理
软件中进行POS数据同名点寻找；
[0013] 步骤SB：通过空三测量计算步骤S1中获取的可见光遥感影像的真实
位置和拼接参数，生成点
云模型；
[0014] 步骤SC：利用步骤SB中获取的真实位置和拼接参数对可见光遥感影像进行自动镶嵌，输出数字表面模型和数字正射影像图；正射影像图空间分辨率为0.1～0.5m,其中数字正射影像图中仅包含红、绿和蓝三个波段信息，影像上树冠边界树种类型清晰可见。
[0015] 进一步地，步骤S2中所述建立树种分类体系的具体内容为：根据目视结果结合森林资源规划设计调查主要技术规定划分影像区域的树种类型。
[0016] 进一步地，步骤S3中所述提取多类型的特征构建特征空间具体包括以下步骤：
[0017] 步骤Sa：直接采用数字表面模型的像元值作为高度特征；
[0018] 步骤Sb：直接采用数字正射影像图的红绿蓝三波段的像元值为光谱特征；
[0019] 步骤Sc：利用红绿蓝三个波段构建常见的8种植被指数，植被指数名称和计算公式分别为：
[0020] 归一化红绿差分指数：NGRDI＝(G-R)/(G+R)；
[0021] 超绿指数：EXG＝2g-r-b；
[0022] 植被
颜色指数：CIVE＝0.441r-0.881g+0.385b+18.78745；
[0023] 植被指数：VEG＝g/rab(1-a),a＝0.667；
[0024] 超绿超红差分指数：EXGR＝EXG-1.4r-g；
[0025] Woebbecke指数：WI＝(g-b)/(r-g)；
[0026] 联合指数：COM＝0.25EXG+0.3EXGR+0.33CIVE+0.12VEG；
[0027] 联合指数2：COM(2)＝0.36EXG+0.47CIVE+0.17VEG；
[0028] 其中R为红波段、G为绿波段、B为蓝波段；r、g、b分别是R、G、B三波段的标准化结果，其中r＝R/(R+B+G)，g＝G/(R+G+B),b＝B/(R+G+B)；a为常数取值为0.667；步骤Sd：将数字正射影像图从RGB色彩空间转到HSV颜色空间，获得HSV三个颜色分量色相、
饱和度和
亮度；转换公式如下：
[0029]
[0030] V＝max
[0031]
[0032] 其中max和min分别对应红、绿、蓝波段的最大值和最小值；
[0033] 步骤Se:以HSV颜色分量为原始数据使用灰度共生矩阵提取纹理特征,其中步长设置为1，方向设置为45°，使用平均J-M距离表达样本间的可分性，用以确定最佳的纹理提取窗口，J-M距离的公式为：
[0034]
[0035] 其中：p(X/ωi)表示第i个像元属于第ωi个类别的几率；
[0036] 步骤Sf：选择特定的形态学属性
滤波器，根据各自设置的
阈值集合{k1,k2,...,kn}分别对数字正射影像的三个波段进行多次的细化和粗化操作，获得单个波段的属性剖面AP，计算公式如下所示：
[0037] AP(I)＝{φkn(I),...,φk1(I),I,γk1(I),...,γkn(I)}；
[0038] 其中φ和γ分别表示形态学中的细化和粗化操作，φkn(I)表示灰度影像I在形态学滤波器下根据阈值kn进行细化得到的特征影像，γkn(I)表示灰度影像I在形态学属性滤波器下根据阈值kn进行粗化得到的特征影像；
[0039] 步骤Sg:将三个波段的特定属性剖面AP联合到一起组合成扩展的形态学属性剖面：
[0040] EAP(S)＝{AP(B1),AP(B2),AP(B3)}
[0041] 其中B1、B2、B3分别代表红、绿、蓝三个波段灰度影像；
[0042] 步骤Sh：选择面积(a)、对
角线长度(d)、灰度值标准差(s)、
转动惯量(i)四种形态学属性滤波器分别构建EAP(S)；
[0043] 步骤Si：对每个形态学属性滤波器获得的EAP(S)进行主成分分析，选择累计信息量大于99％的前n个主分量，作为该属性滤波器下获得的EAP(S)。
[0044] 进一步地，所述步骤S4具体包括以下步骤：
[0045] 步骤S41：使用bootstrap
采样从步骤S2中采集的样本中选择训练样本；
[0046] 步骤S42：利用随机森林算法计算特征重要性；
[0047] 步骤S43：去除重要性最小的特征，更新特征集；
[0048] 步骤S44：利用更新的子集进行随机森林分类，使用总体精度和Kappa系数评价分类精度；
[0049] 步骤S45：判断更新后的特征集是否为空，若为空则统计每种特征组合的总体精度和Kappa系数值，值最大的特征集为最佳特征子集；否则返回步骤S41。
[0050] 进一步地，所述步骤S5的具体内容为在步骤S4获取的最佳特征子集上使用随机森林算法实现树种的分类，提取树种分布图，建立混淆矩阵使用总体精度和Kappa系数评价总体的分类情况，使用F统计值来评价单个树种的分类精度。其中总体精度等于被正确分类的像元总数除以总像元数；Kappa系数的计算公式如下：
[0051]
[0052] m表示分类的类别数,xii表示混淆矩阵第i列第i行上
像素的数量即正确分类的数目，xi+和x+i分别是第i行和第i列的总像素数量，N是用于精度评价的总像素数量；F统计值的计算公式如下：
[0053]
[0054] 其中，PAi表示类别i的生产精度，即从分类结果图中任取一个随机样本，其所具有的类型与地面实际类型相同的条件概率PAi＝xii/xi+，UAi表示类别i的用户精度,即实际任意一个随机样本与分类图上同一地点的分类结果一致的条件概率UAi＝xii/x+i。F统计值越高表明单个地物分类精度越高。
[0055] 与
现有技术相比，本发明具有以下有益效果：
[0056] 本发明从可见光无人机影像上提取多类型的特征，弥补光谱特征不足的问题，同时利用递归消除随机森林算法对多类型特征进行优选，综合不同类型特征中最有益于树种分类的特征实现树种的分类识别，提高了可见光无人机遥感影像森林树种类型识别的精度和自动化程度，为农林业资源调查中树种信息的提取提供了新的途径。
附图说明
[0058] 图2为本发明实施例的采集的无人机影像数据，图2中，(a)为数字正射影像图与样本分布图，(b)为数字高程模型图。
[0059] 图3为本发明实施例的各树种类型的植被指数均值图。
[0060] 图4为本发明实施例的不同窗口大小提取纹理的特征下样本的可分离度图。
[0061] 图5为本发明实施例的递归消除随机森林特征选择结果图。
[0062] 图6为本发明实施例的递归消除随机森林特征选择的特征图，其中，
[0063] 图6(a)为高度特征图，图6(b)为色相纹理均值图，图6(c)为饱和度纹理均值图，图6(d)为EAPa第三主成分图，图6(e)亮度纹理均值图，图6(f)WI指数图，图6(g)色
相图，图6(h)EAPd第五主成分图，图6(i)EAPs第二主成分图，图6(j)EAPi第十一主成分图，(k)EAPd第三主成分图。
[0064] 图7为本发明实施例的分类结果图。
具体实施方式
[0065] 下面结合附图及实施例对本发明做进一步说明。
[0066] 应该指出，以下详细说明都是例示性的，旨在对本
申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
[0067] 需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本
说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
[0068] 如图1所示，本实施例提供一种基于多特征优选的可见光无人机遥感影像森林树种分类方法，包括以下步骤：
[0069] 步骤S1：使用无人机航拍获取所需林区的可见光遥感影像，对可见光遥感影像进行预处理生成数字表面模型和数字正射影像图；
[0070] 步骤S2：建立可见光无人机遥感影像树种分类体系，并从可见光遥感影像上选择不同森林树种及地类样本数据；
[0071] 步骤S3：从步骤S1中生成的数字表面模型和数字正射影像图中提取多类型的特征构建特征空间；
[0072] 步骤S4：利用递归消除随机森林特征选择算法选择从步骤S3构建的特征空间中提取最佳分类特征子集；
[0073] 步骤S5：在步骤S4获取的最佳分类特征子集上使用随机森林算法实现树种的分类，用以提取树种分布图，对分类结果进行精度评价。
[0074] 在本实施例中，步骤S1中所述对可见光遥感影像进行预处理包括以下步骤：
[0075] 步骤SA：在专业无人机处理软件中进行POS数据同名点寻找；
[0076] 步骤SB：通过空三测量计算步骤S1中获取的可见光遥感影像的真实位置和拼接参数，生成点云模型；
[0077] 步骤SC：利用步骤SB中获取的真实位置和拼接参数对可见光遥感影像进行自动镶嵌，输出数字表面模型和数字正射影像图；正射影像图空间分辨率为0.1～0.5m,其中数字正射影像图中仅包含红、绿和蓝三个波段信息，影像上树冠边界树种类型清晰可见。
[0078] 在本实施例中，步骤S2中所述建立树种分类体系的具体内容为：根据目视结果结合《森林资源规划设计调查主要技术规定》划分影像区域的树种类型。由于主要研究对象为树种类型，因此其他地类统称为“其他”类。样本选择时尽可能全面。
[0079] 在本实施例中，数字表面模型表示地物表面的高度分布，因此使用数字表面模型的像元值作为高度特征。数字正射影像中的像元值作为光谱特征，另在数字正射影像图中计算提取植被指数、HSV颜色特征、HSV纹理特征、扩展的形态学多属性剖面特征。六个类型特征构建特征空间；
[0080] 步骤S3中所述提取多类型的特征构建特征空间具体包括以下步骤：
[0081] 步骤Sa：直接采用数字表面模型的像元值作为高度特征；
[0082] 步骤Sb：直接采用数字正射影像图的红绿蓝三波段的像元值为光谱特征；
[0083] 步骤Sc：利用红绿蓝三个波段构建常见的8种植被指数，植被指数名称和计算公式分别为：
[0084] 归一化绿红差分指数：NGRDI＝(G-R)/(G+R)；
[0085] 超绿指数：EXG＝2g-r-b；
[0086] 植被颜色指数：CIVE＝0.441r-0.881g+0.385b+18.78745；
[0087] 植被指数：VEG＝g/rab(1-a),a＝0.667；
[0088] 超绿超红差分指数：EXGR＝EXG-1.4r-g；
[0089] Woebbecke指数：WI＝(g-b)/(r-g)；
[0090] 联合指数：COM＝0.25EXG+0.3EXGR+0.33CIVE+0.12VEG；
[0091] 联合指数2：COM(2)＝0.36EXG+0.47CIVE+0.17VEG；
[0092] 其中R为红波段、G为绿波段、B为蓝波段；r、g、b分别是R、G、B三波段的标准化结果，其中r＝R/(R+B+G)，g＝G/(R+G+B),b＝B/(R+G+B)；a为常数取值为0.667；
[0093] 步骤Sd：将数字正射影像图从RGB色彩空间转到HSV颜色空间，获得HSV三个颜色分量色相、饱和度和亮度；转换公式如下：
[0094]
[0095] V＝max
[0096]
[0097] 其中max和min分别对应红、绿、蓝波段的最大值和最小值；
[0098] 步骤Se:以HSV颜色分量为原始数据使用灰度共生矩阵提取纹理特征其中步长设置为1，方向设置为45°，使用平均J-M距离表达样本间的可分性，用以确定最佳的纹理提取窗口，J-M距离的公式为：
[0099]
[0100] 其中：p(X/ωi)表示第i个像元属于第ωi个类别的几率；
[0101] 步骤Sf：选择特定的形态学属性滤波器，根据各自设置的阈值集合{k1,k2,...,kn}分别对数字正射影像的三个波段进行多次的细化和粗化操作，获得单个波段的属性剖面AP，计算公式如下所示：
[0102] AP(I)＝{φkn(I),...,φk1(I),I,γk1(I),...,γkn(I)}；
[0103] 其中φ和γ分别表示形态学中的细化和粗化操作，φkn(I)表示灰度影像I在形态学滤波器下根据阈值kn进行细化得到的特征影像，γkn(I)表示灰度影像I在形态学滤波器下根据阈值kn进行粗化得到的特征影像；
[0104] 步骤Sg:将三个波段的特定属性剖面AP联合到一起组合成扩展的形态学属性剖面：
[0105] EAP(S)＝{AP(B1),AP(B2),AP(B3)}
[0106] 其中B1、B2、B3分别代表红、绿、蓝三个波段灰度影像；
[0107] 步骤Sh：选择面积(a)、对角线长度(d)、灰度值标准差(s)、转动惯量(i)四种形态学属性滤波器分别构建EAP(S)；
[0108] 步骤Si：对每个形态学属性滤波器获得的EAP(S)进行主成分分析，选择累计信息量大于99％的前n个主分量作为该属性滤波器下获得的EAP(S)。
[0109] 在本实施例中，所述步骤S4具体包括以下步骤：
[0110] 步骤S41：使用bootstrap采样从步骤S2中采集样本中选择训练样本；
[0111] 步骤S42：利用随机森林算法计算特征重要性；
[0112] 步骤S43：去除重要性最小的特征，更新特征集；
[0113] 步骤S44：利用更新的子集进行随机森林分类，使用总体精度和Kappa系数评价分类精度；
[0114] 步骤S45：判断更新后的特征集是否为空，若为空则统计每种特征组合的总体精度和Kappa系数值，值最大的特征集为最佳特征子集；否则返回步骤S41。
[0115] 在本实施例中，步骤S5的具体内容为在步骤S4获取的最佳特征子集上使用随机森林算法实现树种的分类，提取树种分布图，建立混淆矩阵使用总体精度和Kappa系数评价总体的分类情况，使用F统计值来评价单个树种的分类精度。其中总体精度等于被正确分类的像元总数除以总像元数；Kappa系数的计算公式如下：
[0116]
[0117] m表示分类的类别数,xii表示混淆矩阵第i列第i行上像素的数量(即正确分类的数目)，xi+和x+i分别是第i行和第i列的总像素数量，N是用于精度评价的总像素数量。F统计值的计算公式如下：
[0118]
[0119] 其中，PAi表示类别i的生产精度，即从分类结果图中任取一个随机样本，其所具有的类型与地面实际类型相同的条件概率PAi＝xii/xi+，UAi表示类别i的用户精度即实际任意一个随机样本与分类图上同一地点的分类结果一致的条件概率UAi＝xii/x+i。F统计值越高表明单个地物分类精度越高。
[0120] 较佳的，在本实施例中，结合一实例对本实施例进行详细介绍。本实施例中选择的影像区域位于福建省建瓯市房道镇，属亚热带季
风气候，多年平均气温18.7℃，年降
水量1663.8cm，总面积0.16km2。研究区树种以天然生长的阔叶林和人工种植的毛竹
(Phyllostachys edulis)杉木(Cunninghamialanceolata)为主，地物类型较为典型。
[0121] 实验数据为2016年8月16日获取，所使用的无人机为瑞士SenseFly公司生产的小型固定翼eBee无人机，装备SONYWX220相机，1800万像素,等效焦距为24mm。航拍时高度为500m，旁向重叠率为60％、航向重叠率为80％，同时使用差分GPS获取四个地面控制点。
[0122] 将航拍的相片、POS数据和地面控制点输入Pix4Dmapper软件进行影像拼接生成正射影像(DOM)和数字表面模型(DSM)，其工作流程主要包括以下三个部分：首先利用飞行POS数据寻找相片中的同名点；然后通过空三测量计算航拍相片的真实位置和拼接参数，建立点云模型，生成摄影测量点云；最后基于影像内容优化和区域网平差技术，自动校准影像并生成DOM和DSM。影像的分辨率为0.107m，均方根误差为0.039m，摄影测量点云
密度为1.27/m2。研究区影像大小为4511行3260列，总面积为0.16km2，DOM和DSM分别如图2(a)和图2(b)所示。
[0123] 目视结果结合《森林资源规划设计调查主要技术规定》将影像区域的林分组成分为针叶林、阔叶林和竹林，因主要研究对应为森林树种信息，因此其他地类统称为“其他”类。其中阔叶林选择26803个像素、针叶林选择了39601个像元、竹林选择了30039个像元，其他类别选择42486个像元。样本分布如图2(a)所示。
[0124] 利用DOM和DSM提取多类型特征，以DSM的像元值作为高度特征，DOM三波段的像元值作为光谱特征，利用RGB三波段信息构建8中常见的可见光植被指数。统建样本在各个植被指数上的均值，如图3所示，不同植被指数对树种的区分程度不同，其中阔叶与其他树种的区分度较好，针叶和毛竹在多个植被指数上存在混淆。
[0125] 利用色彩空间转换将RGB空间的DOM转到HSV颜色空间，获得色相(Hue)、饱和度(Saturation)、亮度(Value)三个颜色分量。
[0126] 在HSV颜色分量的基础上使用灰度共生矩阵提取纹理特征。纹理特征提取的主要参数有窗口大小、步长和方向，其中步长和方向对分类结果的影响较小，步长设置为1，角度为45°即可。分别取3*3、5*5、7*7、…、25*25共12组纹理特征，计算每个组纹理特征下样本的J-M距离，如图4所示，当纹理窗口大于19之后样本的可分离性趋于稳定，各树种的平均J-M距离均大于1.95，具有较高的可分性，因此选择最佳的纹理窗口为19。
[0127] 提取RGB影像扩展的形态学多属性剖面特征，根据影像上各类型树种树冠的大小，设置属性面积阈值为(5、9、61、115、171、325、1027)、对角线长度阈值推荐设置为(10、20、30、40、50、60、70、80、90)、灰度值标准差阈值设置为(10、20、30、40、50、60、70、80、90)、转动惯量阈值设置为(0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9)。获得4组EPA属性特征，对每组特征进行主成分分析，选择累计信息量达99％的前N个主分量，
叠加生成EMAPs特征。因此EMAPs中分别包含EAPa前四个主成分，EAPs前五个主成分，EAPd的前五个主成分，以及EAPi的前十三个主成分，维度为27。
[0128] 共提取6种类型的特征属性，包括高度特征1个、光谱特征3个、植被指数8个、HSV颜色分量3个、HSV纹理特征24个，扩展的形态学多属性剖面27个，共66个特征，使用递归消除随机森林算法从中进行筛选，如图5所示，当参与分类的特征为11时分类精度最高，且当参与分类的特征继续增加时，分类精度趋于平稳，选择前11个特征作为最佳分类特征子集。这11个特征按重要性从大到小的顺序为DSM、Hue_GLCM_Mean(色相纹理均值)、Sat_GLCM_Mean(饱和度纹理均值)、EAPa_pca3(EAPa的第三主成分)、Val_GLCM_Correlation(亮度纹理相关性)、WI、色相、EAPd_pca5(EAPd的第五主成分)、EAPs_pca2(EAPs的第二主成分)、EAPi_pca11(EAPi的第十一主成分)、EAPd_pca3(EAPd的第三主成分)。特征图如图6所示。
[0129] 为了验证本文方法所优选特征的有效性，分别使用原始RGB影像和优选的11组特征进行分类。特征优选后树种分布图如图7所示。
[0130] 分类精度如表1所示：
[0131]
[0132] 从表1中可知树种分类的总体精度为81.05％，kappa系数为0.73，相比于使用原始RGB进行分类分别提高了15.37％和0.2。阔叶的分类精度提升15.89％、针叶的分类精度提升2.4％、毛竹的分类精度提升18.11％。结果表明进行多特征优选可以极大提高可见光无人机遥感数据树种分类精度。
[0133] 以上所述仅为本发明的较佳实施例，凡依本发明申请
专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。
