融合无人机LiDAR和高分影像的路面平整度监测方法与流程
本发明涉及公路路面平整度遥感监测技术，是利用无人机激光雷达(lightdetectionandranging,lidar)和高分辨率遥感技术针对公路路面进行大范围、高效率、低成本的路面平整度监测与评估。
背景技术：
机载激光雷达(airbornelidar)是一种主动式对地观测系统，集全球导航卫星系统、惯性导航、激光测距及计算机处理技术于一身，能够直接获取高精度、高分辨率的数字地面模型以及地物的三维空间信息，具有传统摄影测量方法无法取代的优越性，并在地形测绘、道路勘测设计、植被参数反演、灾害监测、城市规划及三维建模等领域得到了较为广泛的应用。公路作为一种现代化的运输通道在社会经济中发挥着越来越重要的作用，对沿线的物流、资源开发、招商引资、产业结构调整、横向经济联合等起到积极推动作用，因而其质量好坏直接影响到地区和国家的经济发展。文献[1](sunl,zhangzm,ruthj.modelingindirectstatisticsofsurfaceroughness[j].journaloftransportationengineering-asce,2001,127(2):105-111)中，1960年美国国家公路与交通运输行政联合会开展的道路试验表明：约95％的路面服务性能取决于道路表面的平整度。路面平整度作为评价路面质量的主要指标之一，也是执行公路路面养护管理的重要决策依据。据《2015年交通运输行业发展统计公报》统计，2015年中国公路总里程达到457.73万公里，公路养护里程446.56万公里，占公路总里程97.6％，且呈逐年快速增长趋势。随着我国公路建设的飞速发展，如何实现大范围、高效率、低成本的公路路面健康状况监测就成为一个亟待解决的问题。
平整度监测是常用的公路路面健康状况监测与评价的重要指标之一。平整度评价指标的测量有3米直尺测定最大间隙(h)、国际平整度指数(internationalroughnessindex,iri)、行使质量指数(runningqualityindex,rqi)、颠簸累积值(bumpcumulativevalue,vbi)、平整度标准差(σ)、功率谱密度(powerspectraldensity,psd)等。现有传统的路面平整度监测方法有：3m直尺、连续式平整度仪、手推式断面仪、车载式颠簸累积仪以及近年出现的激光平整度检测仪(包括集成激光平整度仪的道路检测车)，其各自的特点概括于表1。这些方法的共性之一就是只能依靠地面测量获得为数很少的道路纵剖面信息，对道路平整度的描述也较粗略，而且无法快速获取道路表面的三维几何信息，使得管理部门对道路整体状况的了解非常有限。
表1.常用的路面平整度监测方法
低空无人机lidar技术的出现不仅为获取道路表面三维信息提供了可能，而且该技术较少受地形限制、高精度、高效率、低成本(尤其是大范围监测)、可同时获取地面遥感影像的优势更是为路面平整度监测提供了一种新的有效解决方案。目前利用lidar技术进行路面平整度监测的相关案例和研究少见。文献[2](china,olsenmj.evaluationoftechnologiesforroadprofilecapture,analysis,andevaluation[j].journalofsurveyingengineering,2014,141(1):4014011)记载了利用地面激光雷达技术对高速公路进行了平整度监测，但是，该方法并未对路面平整度的空间可视化方法进行尝试，也未对点云的分类方法和路面表面模型的构建方法进行探索，该技术的另外一个致命点就是只能在某个固定点位置附近处获取小范围(半径<50m)的点云数据。如果要获取大范围的路面信息则需大量的点云拼接、坐标纠正、点云配准、坐标转换等一系列复杂、误差不断积累的后处理过程，不仅流程复杂而且非常低效。文献[3](郭姣.基于车载激光的道路平整度检测系统研究[d].首都师范大学,2013)记载了利用车载激光雷达技术分析道路剖面平整度特征并定位道路损坏超标位置，但是，该方法同样没有对点云的分类方法进行探索，而是借助传统的点云处理平台来完成分类工作。此外，该方法仅仅依靠单个数值进行表达，未能对路面质量的空间分布特征进行探索，并要求检测车行驶速度稳定、gps信号稳定，而这在某些建筑物遮挡、崎岖不平的复杂路段难以应用。
综上所述，现有传统的公路路面平整度监测方法大多存在效率低下、自动化程度低，无法进行大范围全路段监测，并且只能获得为数很少的道路纵剖面信息，不能得到道路表面的三维结构信息。现有的激光雷达技术的平整度监测方法则存在测量范围小、作业效率低、缺乏空间可视化表达、数据源单一的缺陷，并且会受到地形条件限制。然而公路路面平整度监测迫切需要探索新的技术途径，提高大范围路面平整度监测效率与准确性。
技术实现要素：
针对现有公路路面平整度监测方法中的不足，本发明提出一种融合无人机lidar和高分辨率影像的路面平整度监测方法及其实施方案。首先，对无人机lidar系统获取的点云数据和高分影像数据进行预处理，然后从高分影像数据提取光谱特征，对点云数据提取多尺度几何特征。将所得到的几何特征和光谱特征作为点云分类的特征变量，并对这些特征变量进行去冗余和降维操作，在特征级别上对两类数据进行融合，利用机器学习中的随机森林分类算法进行点云分类；再进一步通过滤波操作和粗差纠正得到比较精准的道路点云数据集。然后利用三维插值法构建高精度的精细路面数字表面模型(digitalsurfacemodel,dsm)。基于iri指数原理计算路面平整度，再利用gis进行空间可视化表达与平整度质量评价。相比于以往的路面平整度监测技术只能用单个iri值来表达整个路段信息而言，该方法则可给出某条路不同路段不同车道iri值的空间分布，提高了路面养护管理工作决策的科学性，有效提高了无人机遥感技术在公路养护管理应用中的实用性，可广泛应用于城市道路及其他低等级道路的路面平整度监测需要。
本发明的原理是：无人机lidar遥感系统获取的点云数据中包含公路及周边地物的三维坐标信息，通过选择合适的点云去噪手段可有效去除噪声干扰，而同期获得的高分辨率影像数据富含地物的光谱信息且非常易于目视判读。点云数据则因为点集分散而不易直观判读，为此将高分影像与点云数据进行融合就能同时将光谱信息与空间几何信息进行集成，为点云分类提供更多的可利用信息。目前，针对高分影像的地物分类存在一个挑战，就是对“异物同谱”的地物目标无法区分，而点云数据的空间几何信息则可以对此起到很好的补充。同样的，点云数据自身由于点集非常离散而导致在一些边缘交界处以及几何结构特征非常相似的地物之间出现混淆，而高分影像则能对此进行很好的辅助作用。基于此，本发明一方面从高分影像提取典型地物类别的光谱特征，另一方面从点云上提取典型地物类别的几何特征，通过构建多尺度几何特征来强化地物类型的空间形态差异，然后采用高效且抗噪声强的机器学习分类算法对点云进行分类，继而对分类后的点云进行粗差纠正和类别滤波即可得到较精准的道路点云数据，再结合点云内插算法构建精细的道路dsm。在此基础上，从道路表面选取一系列纵剖面曲线，分段计算国际平整度值并进行空间可视化表达，最后对平整度值进行分等定级，从而实现对道路平整度质量的评价。本发明方法经过精度评价证明可用来进行大范围公路路面的快速自动化平整度质量状况监测。
本发明提供的技术方案：
一种融合无人机lidar与高分影像的路面平整度监测方法，包括：数据获取与预处理、特征提取与融合、点云分类、路面dsm模型的构建及路面平整度评价等技术环节。
其中，数据获取与预处理的目的是对低空无人机lidar系统获得的点云数据进行去噪，对同期获取的高分影像进行配准与拼接，得到实验区的正射影像数据，为数据融合做准备。数据融合的目的是将点云数据的几何空间信息与高分影像数据的光谱信息进行集成，从而使得点云数据同时具备空间几何信息和光谱信息。点云分类的目的是为了得到道路点云数据，数据融合的过程为点云分类提供了可利用信息，对实验区进行典型地物类别划分，构建一系列几何特征参数来表达地物类型的差异，并结合光谱特征参数就构成了分类器的特征变量，接着使用随机森林算法实现对点云数据的分类，得到路面点云数据。路面dsm模型的构建可以为路面三维空间结构的表达提供基础，而且dsm模型的每个栅格单元存储的都是高程信息，便于平整度计算，也为平整度的空间可视化表达提供了条件。路面平整度评价的目的就是通过选取一定的平整度指标对路面的平整度状况进行量化描述，并建立分等定级体系对实验区的道路状况进行客观评价，从而对路面当前的整体状况给出前瞻性分析和决策。
本发明提供的路面平整度监测方法具体包括以下步骤：
1)利用低空无人机lidar系统(包括：无人机平台、lidar扫描仪、多光谱相机、姿态定位单元、自动驾驶仪等单元的集成)，获取实验区的lidar点云数据和高分辨率影像数据；
2)对获取的lidar点云数据转换为当前通用的点云数据文件格式(*.las)，然后对点云数据和无人机拍摄的高分辨率影像数据进行预处理，包括点云的去噪、影像的拼接与配准处理及正射影像数据的生成；
3)对步骤2)中得到的正射影像的波段进行分析并提取反映地物差异的光谱特征参数，连同高分辨率影像原来的波段，构成包含光谱信息的多个光谱特征集；
4)对得到的光谱特征数据与点云数据进行特征融合，并以点云为空间参考基准进行地理空间配准；
5)对实验区的主要地物类型进行划分，并在步骤4)得到的点云数据基础上经过抽稀得到抽稀的点云数据，然后对抽稀的点云数据选取一系列的几何特征参数进行计算，并通过改变空间分析的尺度得到包含多尺度几何特征参数和光谱特征参数的点云数据，接着从中选取部分相应地物类型的点云数据作为后续点云分类的训练样本，将剩余的抽稀点云数据作为待分数据。点云的抽稀过程减少了几何特征参数计算过程中的运算量，同时也降低了分类所需处理的数据量，大大提高了数据处理效率；
6)步骤5)中得到的特征变量通常会很多，不可避免的会出现某些变量之间存在较强的相关性，导致特征变量冗余。为了找到关键特征变量，从而对数据起到降维作用，进而简化分类器构建过程的复杂度，在特征级别上，融合多光谱影像的光谱信息与激光点云数据的空间几何信息，构成新的点云分类特征集；
7)利用随机森林模型基于步骤6)的特征集对步骤5)的训练样本进行模型训练，通过不断的调整模型参数取得满意的分类效果后，应用于待分类样本点云数据进行分类，然后将样本点云的类别属性按照空间最近邻插值原理赋给经去噪后的点云数据，从而完成全部点云的分类。分类结果可能还包含不属于路面的点云，此时要对分类后的点云数据进行人工检查和粗差纠正，并进一步通过类别属性过滤即得到道路点云数据；
8)对步骤7)得到的道路点云数据的高程信息，采用不规则三角网模拟路面起伏变化，然后通过自然邻域插值法构建规则格网的高精度精细道路路面的dsm；
9)从步骤8)中得到的dsm模型数据中选取一系列的纵剖面线，运行iri指数模型计算各剖面线对应的iri值并经过栅格化处理，得到路面的iri值的空间分布图；
10)对步骤9)中得到的iri值，依据道路平整度的分等定级标准进行分类，从而得到待测路段对应的平整度质量评价结果。
针对上述公路路面平整度监测方法，进一步地，步骤2)中涉及的lidar点云数据，具体为满足以下4个条件的点云数据：
a)坐标系：点云数据的坐标系为wgs84坐标系(经度、纬度)或者与当地相匹配的地方平面投影坐标系(x，y)；
b)属性信息：点云数据属性信息除经度、纬度外，还应包括高程、扫描角、反射强度；
c)点云密度：点云密度在400点/平方米以上(即激光脚点半径在50mm以内)；
针对上述平整度监测方法，进一步地，对步骤2)中涉及的数据预处理过程，主要是指以下4种：lidar点云的噪声点过滤、lidar点云的高程异常点过滤、lidar点云的扫描角过滤与高分辨率影像的拼接与配准。针对上述平整度监测方法，进一步地，对步骤2)中涉及的高分辨率影像数据，是指市面上性能良好的多光谱数码相机所拍摄的高清数码照片，并且无人机能提供相匹配的姿态与定位(pos)信息文件。针对上述平整度监测方法，进一步地，对步骤3)中涉及的光谱特征参数，为以下3类：
a)高分影像的红、绿、蓝波段的灰度值。
b)高分影像的红、绿、蓝波段的灰度值之标准差；
c)激光点云的反射强度值(lidar扫描仪发射的激光脉冲也属于电磁波的一种，同样具有波长信息)。
针对上述平整度监测方法，进一步地，对步骤4)中涉及的特征融合，具体是指将点云数据与高分辨率影像数据的光谱特征参数进行融合，并且融合的过程是以点云数据的空间参考为基准的，采用仿射逆变换模型将影像上每个像元的光谱值匹配到点云数据坐标系；
针对上述平整度监测方法，进一步地，对步骤5)中涉及的几何特征参数，为以下3类：
a)局部粗糙度(localdegreeofroughness,ldr)。该特征变量的具体含义是指点云中的某个参考点到该参考点所在的某个空间尺度下的邻域点所形成的最佳拟合平面的距离；
b)局部维度特征(localdimensionfeature,ldf)。文献[4](brodun,lagued.3dterrestriallidardataclassificationofcomplexnaturalscenesusingamulti-scaledimensionalitycriterion:applicationsingeomorphology[j].isprsjournalofphotogrammetryandremotesensing,2012,68:121-134)记载了局部维度特征ldf，该特征变量的具体含义是指点云中的某个参考点与该参考点的邻域点之间在空间上是以一维、二维还是三维分布的特征度量值，用一个数值来描述它分属于这三种维度的大小程度，该值就是局部维度特征参数；
c)局部高程差(localheightdifference,lhd)。该特征变量的具体含义是指点云中某个参考点与该参考点在某个空间尺度下的邻域点所构成的全体点集中，到达该点集所形成的最佳拟合面的最大与最小距离之差，即为局部高程差，起伏越大的平面该值也越大。
针对上述平整度监测方法，进一步地，对步骤5)中涉及的多尺度几何特征，其具体含义是指，在计算某个几何特征值时其空间分析尺度有多个(固定某个空间分析尺度会得到一个对应的空间几何特征参数)，不同地物类型在不同的空间分析尺度下其几何特征值会有不同。
针对上述平整度监测方法，进一步地，对步骤5)中涉及的点云数据类别，主要包括道路点云(含交通标线)和低矮植被(农作物和草地)、土壤、树木、移动目标、电力线(包括电线杆)这6大类。
针对上述平整度监测方法，进一步地，对步骤5)中涉及的点云抽稀过程，是在空间中按照一定的空间间隔进行采样，该过程要保证点云在空间的分布相对于原始点云而言是均匀的(而不是零乱的分布)。
针对上述平整度监测方法，进一步地，对步骤6)中涉及的特征选择，主要包括特征子集的搜索和特征子集的评价两个过程。特征子集的搜索采用的是贪心策略，即：初始的最优子集为空，所有待评估的特征变量则组成备选特征集，从备选特征集中筛选解释能力最强的特征变量加入到最优特征子集中，一次循环只添加一个变量；某一特征变量添加结束后需要从备选特征集中剔除该特征变量，重新下一次循环。最优特征子集的评价标准：新加入的特征变量与已选中的特征变量之间的相关性较低。
针对上述平整度监测方法，进一步地，对步骤8)中涉及的dsm模型的分辨率要尽可能高，并且其高程精度应控制在15mm以内。具体的水平和高程精度评估准则为：从现场的路面中选择一些具有明显特征的目标物(如坑槽)并记录其对应的gps坐标，通过测量其长宽高属性，然后在高分辨率影像中找到对应的目标物，利用影像与lidar点云的融合，易于在dsm模型中找到该目标的对应位置，进而可利用长度测量工具交互式地测定该目标的长宽高，再统计分析完成对dsm模型的水平和高程的精度评估。
针对上述平整度监测方法，进一步地，对步骤9)中涉及的道路纵剖面线，其间距越小则对应的路面平整度空间分布越精细。iri指数模型的理论依据是四分之一车模型，在计算过程中需要严格控制采样间隔(50mm以内)和测量单元的长度(10m或20m)。最佳的采样间隔和测量单元长度可通过蒙特卡罗模拟方法定量分析在高程误差一定时iri的误差随采样间隔和测量单元长度的变化规律来确定。
针对上述平整度监测方法，进一步地，步骤9)中涉及的栅格化是指将路面纵剖面线的iri值按照最邻近原则赋给周围的栅格单元，从而得到待测路段的平整度空间分布图。
针对上述平整度监测方法，进一步地，通过在待测路段选择一系列包含不同平整度质量等级的路面样本作为参考数据，可对本方法的评价结果进行精度评定，具体步骤如下：
1)根据设定的平整度质量分级标准，同期在待测路段选择一系列(包含平整度质量较好与较差，覆盖健康与病害类型，每种路面选择一定数量的地面实测样本)的路面样本作为参考数据；
2)对参考数据选取外观、表面粗糙度、车辙深度、病害面积、损坏程度这5项指标使用专家评分法对其进行平整度质量评价，并作为参考值；
3)运用本方法同样的对此路面实测样本进行平整度值进行计算，依据平整度质量分级标准，得到本方法的路面平整度质量的评价结果；
4)对平整度质量结果的评价可以类比于遥感影像分类结果的评价，而混淆矩阵是用来评价遥感影像分类精度的最常用形式，因此本方法也采用混淆矩阵对路面平整度质量结果进行评价。引进精度评价指标总体精度(overallaccuracy,oa)和kappa系数对结果进行分析，得到本方法的精度。
与现有技术相比，本发明的有益效果是：
本发明基于低空无人机平台获取的lidar数据和高分辨率遥感影像，首次提出融合影像的光谱特征与lidar点云的多尺度几何特征构建点云分类特征集的方法，解决当前路面点云与非路面点云分类的难题；再利用路面点云构建精细三维路面数字表面模型(dsm)，通过引入国际平整度指数的计算原理，快速计算获得公路路面的平整度参数，实现路面质量的快速监测与评估。本发明首次提出了一种融合低空无人机lidar和高分影像的公路路面快速高效的平整度监测方法，有效解决了现有路面平整度监测方法中效率低下、自动化程度较低，不适合大范围全路段监测，且只能获得为数很少的道路纵剖面平整度信息，不能得到道路三维结构信息等不足。本方法可快速掌握大范围的路面平整度状况，从而为路面养护维修决策提供科学依据。同时，本发明设计的点云数据处理流程简单易实现，效果显著。另外，本发明提供的平整度质量空间分布图能全面反映公路路面的平整度状况，相比于以往的平整度监测技术只能依靠单个数值来表达整个路段的平整度信息而言，有力地促进了道路养护管理的信息化与精细化发展，提高了路面养护管理工作决策的科学性。
本发明融合影像和lidar进行路面平整度监测，与现有技术相比，一是采用了低空无人机lidar数据源，并首次提出融合光谱特征信息与多尺度几何特征信息来优化道路点云分类，在此基础上构建精细的dsm，快速提取具有空间分布的路面平整度信息；二是数据处理过程简洁易实现，使用的数学模型易于理解，效果显著，所有操作均可在个人普通电脑上完成，对硬件的要求不高；三是在提供可行性的操作流程上，进一步给出了优化数据处理的方法，并通过实验完成验证，具有很好的可重复性。
附图说明
图1为本发明提供的路面平整度监测方法的流程图
图2为局部粗糙度计算方法示意图。
图3为局部维度特征的简化表达示意图；
其中，(a)为一维空间分布特征，(b)为二维空间分布特征，(c)为三维空间分布特征。
图4为局部高程差计算示意图。
图5为本方法中进行点云插值的自然邻域插值算法示意图。
图6为本方法中使用的国际平整度指数(iri)的算法示意图；
其中，zs表示簧载质量的垂直位移，zu表示非簧载质量的垂直位移，ms表示簧载质量的大小，mu表示非簧载质量的大小，ks表示连接簧载质量和非簧载质量的弹簧的劲度系数，cs表示连接簧载质量和非簧载质量的线性阻尼系数，ks表示连接轮胎的弹簧的精度系数，可模拟轮胎所产生的减震效果，y(x)表示路面高程，b代表的是轮胎对地面的包容长度(轮胎与地面接触的部分)。
图7为本发明实施例中地面采集的不同平整度特征的路面样本照片
其中，(a1)坑槽，特征为石料暴露且呈大块状分布，形状近似为圆形或者椭圆形，平整度质量评价为极差(failed)；(a2)坑槽，特征为石料暴露但面积相对较小，形状近似圆形或椭圆形，平整度质量评价为极差(failed)；(b1)塌陷，特征为条带状分布且深度较大(石料没有暴露)、覆盖面积大，平整度质量评价为极差(failed)；(b2)塌陷，特征为条带状分布且深度较小(石料没有暴露)、覆盖面积小，平整度质量评价为差(poor)或极差(failed)，具体取决于沉陷深度以及覆盖面积；(c)裂缝(包括细裂缝和粗裂缝)，特征为狭长状分布、裂缝间隙在10mm以内，平整度质量评价为中等(fair)或差(poor)，具体取决于裂缝宽度和面积；(d)麻面，特征为路面非常粗糙或者有许多细微的凹坑，平整度质量评价为中等(fair)或差(poor)，具体取决于表面破坏的程度和面积；(e)沙石覆盖，特征为路面出现块状破损且被沙石覆盖，平整度质量评价为差(poor)；(f)平整度一般的旧路面，特征为路面纹理较粗糙但整体表面流畅，平整度质量评价为中等(fair)；(g)平整度较好的旧路面，特征为路面色泽较为光亮、纹理致密、表面非常流畅，平整度质量评价为好(good)。
图8为本发明实施例中实验区内获取并经过格式转换lidar点云数据。
图9为本发明实施例经影像拼接与配准处理的实验区多光谱正射影像。
图10为本发明实施例中融合影像波段信息的lidar点云数据(图中显示的是蓝波段)。
图11为对本发明实施例中使用5个光谱信息和36个多尺度几何特征信息进行特征选择后所保留的23个关键特征变量及训练样本在这些特征变量下的属性值分布特点；
其中，(a)不同地物的反射强度、r-g-b波段灰度标准差、r-g-b波段灰度值等5个光谱特征；(b)不同地物的局部粗糙度(ldr)与局部维度特征(ldf)；(c)不同地物的局部维度特征(ldf)和局部高程差(lhd)；(d)不同地物的局部高程差(lhd)。括弧里的长度数字代表空间分析尺度大小，ldf_1表示一维的局部维度特征，ldf_1表示二维的局部维度特征
图12为本发明实施例中对路面点云插值生成tin表达的dsm。
图13为本发明实施例中利用自然邻域插值算法得到的道路dsm。
图14为本发明实施例利用iri指数模型计算得到的路面平整度空间分布图。
图15为本发明实施例对平整度值进行重分类后得到的路面平整度质量空间分布图。
图16为本发明实施例中不同路面质量类型的空间分布统计图。
图17为评估本发明实施案例所建立的平整度质量评价结果的混淆矩阵示意图。
具体实施方式
下面结合附图，并通过具体实例进一步描述本发明的具体实施过程，但不以任何方式限制本发明的适用范围。
图1为本发明提出的融合无人机lidar和高分影像的路面平整度监测方法的流程框图，包括以下步骤：
步骤1：利用无人机lidar系统同时获取待测路段lidar点云数据和高分辨率影像数据；
步骤2：由于不同的lidar系统生产厂商提供的点云数据格式不同(生产厂商大都有一套自定义的存储格式)，为了能在一些通用软件上进行处理，因此需要对lidar点云数据进行格式转换，在输出时注意投影坐标系的选择以及保留的字段属性；
步骤3：对步骤2中得到的通用点云文件(*.las)进行去噪处理，初步去除干扰噪声信息，从而得到去燥后的点云数据。具体的去噪处理方式包括：
a)噪声点过滤。主要是对lidar点云中由于激光脉冲的折射或者多路径效应产生的异常点进行剔除；
b)高程异常点过滤。主要是剔除局部高程异常点(粗差)，具体的判断方法为：设置一定的搜索半径阈值r，如果某点与该点所在半径r内的邻域点的高程平均值相差超过3倍的标准差，则可视为高程异常点并加以剔除；
c)扫描角过滤。由于在无人机航迹规划时，为了能全面获取路段的表面信息，无人机的飞行路线应该在路段正上方，因此利用点云的扫描角属性可以对数据进行过滤从而减少数据处理量，节省不必要的数据处理时间。假设某一点到飞机正下方的水平距离为d，飞机的飞行高度为h，则该点的扫描角θ定义为式1：
θ＝arctan(d/h)(式1)
根据道路的边界范围可以确定道路的最大扫描角θmax，因此设置合适的扫描角阈值就可将道路区域分离。考虑到航线有一定的偏差，实验中需要尽可能将扫描角阈值设置稍大一些；
步骤4：利用无人机影像拼接软件(如本实例中使用的由瑞士pix4d公司开发的商用软件pix4dmapper)，对待测路段内无人机拍摄的高分辨率影像数据进行影像拼接和配准处理，可得到待测路段的正射遥感影像；
步骤5：对高分辨率影像进行光谱特征提取，本发明用到的光谱特征参数包括以下4个：红、绿、蓝波段的灰度值以及这三个波段的标准差，标准差的数学模型为式2：
其中，std.代表标准差，gi(i＝1,2,3)是红、绿、蓝三波段的灰度值。
步骤6：对包含多波段高分影像与点云数据进行特征融合，以点云数据的地理空间参考为基准，通过仿射变换将影像的像素坐标映射到该像素中心点的地理坐标，从而将影像上对应的光谱信息匹配到点云数据。仿射变换的数学模型为式3：
但在实际操作中，由于影像的像素数量远远大于点云的数据量，而且两者也不完全重合，因此为了减少计算量，本发明采用的是从地理坐标转换到图像坐标的过程，即仿射反变换，表示为式4：
式3、式4中，g即是仿射变换矩阵，ai,bi,ci(i＝1,2)是仿射变换系数，mapx、mapy分别是点云的x坐标、y坐标，row、col分别是影像的行、列值。geotiff或其他带有地理坐标信息的影像文件的元数据中一般就有仿射系数信息，另外也可通过影像上的4个边界点坐标及其对应的行列数(像素坐标)信息通过构建方程来对仿射系数进行求解而得到。
步骤7：对点云数据构建八叉树进行存储，以八叉树格网单元的长度为最小空间间隔对点云进行抽稀，得到抽稀的点云数据。然后对抽稀的点云数据进行几何特征参数的计算，本发明中用到的几何特征参数包括以下三个：
a)局部粗糙度(localdegreeofroughness,ldr)，如图2所示，假设整个三维点云点集为c＝{pi|pi＝(xi,yi,zi),i＝1,2,3,...,n}，当前计算点为p＝(x,y,z)∈c，三维空间中其半径为r的邻域点集为p＝{pj|||pj-p||≤r,pj≠p,j＝1,2,3,...,n}，设平面t:ax+by+cz+d＝0是点集q＝{p∪p}的最优拟合平面，那么ldr的数学形式可以表述为式5：
式5中的a,b,c,d分别代表平面方程的基本参数，x,y,z分别代表当前搜索点的x,y,z坐标。本发明中的最优拟合平面t是指点集中所有点投射到x-y平面上时，使得z值误差平方和s为最小的平面，具体的数学形式为式6：
式6中xi,yi,zi分别代表第i个点的x,y,z坐标，其他参数含义同式5；
b)局部维度特征(localdimensionfeature,ldf)，如图3所示，三维点云点集和邻域点集的构造类似于局部粗糙度，假设当前搜索点为p＝(x,y,z)∈c，令点集q＝{p∪p}＝[xyz]，首先对q进行主成分变换，得到矩阵q的三个主成分系数μ1,μ2,μ3(μ1≥μ2≥μ3)，通过式7进一步对这三个主成分系数进行归一化：
其中，λ1、λ2、λ3分别对应于当前搜索点在该邻域内服从一维、二维、三维空间分布的程度，也即一维、二维、三维特征值，由于λ1+λ2+λ3＝1，因此只需要前两个特征参数即可表达当前搜索点的局部维度特征，这样可以提高点云数据的几何特征计算效率；
c)局部高程差(localheightdifference,lhd)，如图4所示，其中最优拟合平面的定义与局部粗糙度计算中的最优拟合平面定义一致，假设当前搜索点为p＝(x,y,z)∈c，三维空间中其半径为r的邻域点集为p＝{pj|||pj-p||≤r,pj≠p,j＝1,2,3,...,n}，点集q＝{p∪p}是包含当前搜索点的邻域点集，那么lhd的数学形式可以表述为式8：
式9中
实验过程中为了加快计算搜索点的局部高程差的速度，对q进行svd分解可以快速得到平面t的4个关键参数a,b,c,d，然后寻找距离平面t最大的点和距离平面t最小的点进行距离求差，从而得到搜索点的局部高程差。
通过不断改变搜索邻域的半径r的大小，并依次计算各个几何特征(ldr、ldf和lhd)在不同搜索空间半径r下的对应值，可以得到关于这三个几何特征的一组列向量，这组列向量就是对地物目标(搜索点所在的物体)的多尺度几何特征参数。
步骤8：对高分辨率影像与点云数据进行融合后，此时的点云数据包含了大量的特征变量(几何特征参数和光谱特征参数)，形成一个高维的矩阵(每个点本身就是一个3维向量)，本身这些特征变量之间也可能存在相关性，导致信息会产生一定的冗余，去除冗余信息的过程就是特征选择的过程。本发明使用的特征选择策略分为两个过程：特征子集搜索和特征子集评价。特征子集搜索的策略为前向搜索，而特征子集的评价则主要是借助相关系数来评估子集内部的相关性程度。具体的实现过程如下：
1)首先建立两个集合形式的数据结构，假设分别为p和q，p代表候选子集，q代表最优子集，类别属性变量为y，除类别属性外的其他属性为x＝{x1,x2,x3,...,xn}。将除y外的其他特征变量添加到候选子集中，设置子集内部的相关性阈值为δ；
2)如果
式10中，xi、y分别代表第i个特征、分类属性，
3)将待选特征变量xc与q内的每个特征变量进行相关性评价，如果xc与q内的任意一个特征变量的相关性均小于δ，则将xc添加到q内，同时从p内删除该特征变量，否则直接从p内删除该特征变量，继续步骤2)；
4)输出q，结束特征子集搜索。
特别说明的是，如果特征变量个数不是很多(比如30个以内)该步骤也可省略。
步骤9：获得抽稀后的点云数据后，交互选择其中的一部分作为训练样本，训练样本要求覆盖所有的地物类型，并且训练样本不能过于集中，适当分散在点云所在空间内，每种地物类型的训练样本数量不宜太少，否则分类器无法构建分类模型。随机森林模型具有对数据集的适应能力强(能同时处理离散型、数值型数据)、训练速度快、抗噪声能力强、不容易陷入过拟合、模型参数少、分类效果好等优点而成为机器学习领域内最受欢迎的算法之一，本发明借鉴其优点来对同时融合多尺度几何特征和光谱特征的点云数据进行地物类别分类。
该算法的基本思想为：假设输入样本数据的个数为n(每个样本对应一行)，特征变量的个数为m(每个特征变量对应一列，即原始样本数据大小为n*m的矩阵)，初始时随机有放回的从n个样本中选择n个样本作为一棵树(随机森林中每一棵树都是二叉树)的训练集，该树的根节点就是训练集；随机地从m个特征中选择m个特征变量，然后从这m个特征变量中按照纯度最小的原则选择特征变量进行分裂，从而可将该树分为两个节点，这两个节点各自包含训练子集的一部分，如此递归的进行划分，直至该树无法继续分裂(到达二叉树的最大递归深度)或者这m个特征使用完毕为止。重复刚才构建树的过程可得到含有任意棵树的森林(随机森林)，并且每棵树的训练样本都不完全一样，使用的特征变量也不完全一样。在预测类别阶段，当输入一个待分类样本时，森林中的每棵树对该样本都会得到一个决策类别，对决策类别进行统计便可得到待分类样本在各个预测类别下的分布频率，输出其中频率最大的类别作为其预测类别，从而完成分类。
本发明中使用的纯度度量指标为基尼系数，采用文献[5](陈华舟,陈福,石凯,等.基于随机森林的鱼粉蛋白近红外定量分析[j].农业机械学报,2015(05):233-238)记载的计算方法，其数学模型为：
1)假设训练样本集包含p个类别，第i个类别的样本数比例为pi，则其基尼系数gini(p)值为：
对该训练样本集按照某个特征变量划分为m个子集，第i个子集的样本个数为ni，该训练样本集的样本总数为n，则划分后子集的基尼系数ginisplit(p)为：
步骤9：对抽稀后的点云进行分类，依据类别属性对其进行过滤可获得路面点云数据，进一步通过相关的点云内插算法就可以得到待测路段的dsm数据。为减少高程噪声对表面模型效果及后续平整度计算的影响，本方法先利用不规则三角网(tin)生成三维地形数据，然后采用自然邻域插值法构建道路的dsm模型。文献[6](sibsonr.abriefdescriptionofnaturalneighbourinterpolation[j].interpretingmultivariatedata,1981:21-36)记载的自然邻域插值法通过式13计算得到待插值点的高程值：
式13中的wi代表邻域内第i个高程点的权重系数，zi则是邻域内第i个高程点的高程，i代表邻域内第i个高程点，z为待插值点的高程值。该方法确定邻近点的原则和距离无关，而是与待插值点的泰森(voronoi)多边形是否有交集来确定邻近关系，而且邻近点的权重则是待插值点的泰森多边形与邻近点的泰森多边形的公共面积所占比例(如图5)。需要注意的是：待插值点的泰森多边形构建过程中考虑的是所有点，而邻近点的泰森多边形构建过程中则将待插值点排除在外。该插值方法的基本特点是它具有局部性，仅使用待插值点周围的样本点子集，权重系数保证了插值得到的结果在所使用的样本值范围之内。其优点在于不会出现外推趋势且不会生成输入样本点尚未表示的山峰、凹地、山脊或山谷，而且得到的表面在除输入样本位置之外的其他所有位置均是平滑的；
步骤10：利用得到的道路dsm选取一系列一定间距(如0.2m)的纵剖面线，然后利用国际平整度指数iri模型来表达实验路段的平整度值。本发明具体实施中，采用文献[7](sayersmw,gillespietd,queirozca.theinternationalroadroughnessexperiment[j].1986)记载的方法计算得到路段平整度值。iri指数的核心原理为四分之一车模型(如图6)，当四分之一车辆以一定的速度沿道路纵剖面行驶时，在路面坡度的激励作用下，系统将产生振动，计算其以规定速度(如80km/h)行驶一定距离后(如1km)悬挂系统的累积位移值即为iri，单位为m/km。为求解该悬挂系统的相对位移，建立二阶振动微分方程如式14～式15：
其中，y是路面高程，zs和zu分别代表簧载质量位移和非簧载质量位移，μ、c、k1、k2是系数，且依据文献[7]记载，有：
其中，ms表示簧载质量的大小，mu表示非簧载质量的大小，ks表示连接簧载质量和非簧载质量的弹簧的劲度系数，cs表示连接簧载质量和非簧载质量的线性阻尼系数，ks表示连接轮胎的弹簧的精度系数。
构造状态变量
其中：
式16是一个非齐次线性微分方程，其对应的状态方程为：
其中，
pr＝a-1(st-i)b,(式19)
式18、式19中的i是单位矩阵，结合
文献[8](王建锋,宋宏勋,马荣贵.路面平整度评价指标iri的影响因素[j].重庆交通大学学报(自然科学版),2012,31(06):1145-1148)记载了上述计算iri的方法，但使用的数据为传统地面测量数据。本实例中为实现iri指数的计算过程，首先需要将dsm数据进行等距离间隔采样处理得到路面高程点序列y＝{y0,y1,...,yn}，并将其作为模型的输入数据。然后依据矩阵a和b计算状态矩阵st和pr，然后对y进行一阶差分运算得到
步骤11：步骤6的计算结果仅仅是得到了单个剖面线的iri值序列，为了能模拟整个路段表面的平整度值，可对道路dsm上的每个栅格单元按照最邻近原则分配对应的iri值，这种处理方式的误差大小主要取决于纵剖面线的间距大小，间距越小则误差也越小。经过上述过程便可得到整个待测路段的平整度空间分布。
步骤12：参照如下标准(如表2)对得到的路面iri值进行分类，即完成对路面平整度质量的评价。该分级设置标准需要结合我国《公路沥青路面养护技术规范》的路面养护质量的相关要求来设置，本方法是结合《公路沥青路面养护技术规范》(jtj073.2-2001)规范来确定的。通过统计待测路段各个质量等级的像元数目可计算其对应的分布比例，并结合道路养护维修的相关规范，从而对该路段需要何种维护措施做出科学决策。
表2路面平整度质量分级设置标准
为了对本方法的可靠性进行评价，采用下述步骤对该方法路面平整度质量监测结果进行精度评价。
1)根据设定的平整度质量分级标准，在待测路段选择一系列(包含平整度质量较好与较差，包含健康与病害类型，每种路面选择一定数量的地面实测样本数据)的路面样本作为参考数据；
2)对参考数据选取外观、表面粗糙度、车辙深度、病害面积、损坏程度这5项指标使用专家评分法对其进行平整度质量评价，并作为参考值；
3)运用本方法同样的对此路面样本进行平整度值进行计算，依据平整度质量分级标准，得到本方法的对应评价结果；
4)对平整度质量结果的精度评价可以类比于遥感影像分类结果的精度评价，而混淆矩阵是用来评价遥感影像分类精度的最常用形式。因此本方法也采用混淆矩阵对路面平整度质量结果进行评价。借助分类精度评价指标总体精度(oa)和kappa系数对结果进行分析，得到本方法的精度，具体的精度评价指标计算方法为：
式21、式22中，n是样本总数，nii代表本方法的评价结果与参考评价结果均为i等级的样本数量，ni.则表示参考评价结果为i等级的样本总数，n.i则代表本方法的评价结果为i等级的样本总数。
以下以某市一条县级公路为例说明本方法具体实施及精度评价过程。
(1)数据获取
该县级公路所在地理位置为北纬44°24′47″、东经85°53′47″附近，道路宽度约为8m。研究人员于2016年6月23日借助瑞士生产的scoutb1-100无人直升机搭载的激光扫描仪系统(rigelvux-1lr)在实验区正上方进行实验数据采集，该激光扫描仪的基本参数信息如表3所示(来源：http://www.riegl.com/products/newriegl-vux-1-series/newriegl-vux-1lr/)。为了使得获取的点云数据密度足够高，实验中无人机飞行的距离路面高度设置为30m、飞行速度为5m/s、扫描角为110°、扫描频率为550khz，最后得到的激光脚点密度为300-600pts，扫描幅宽为85.7m，数据获取当天晴朗无(少)云。
表3激光扫描仪系统的关键参数信息
在获取lidar点云数据的同一天，研究人员在飞行区域内进行地物特征点与路面样本数据的采集。地物特征点数据主要是沿道路周边采集，通过选取一些关键的地物目标(如车身、标志线、坑槽、大块裂缝)进行几何尺寸度量，利用实时差分gps设备进行定位并进行笔录存档；然后在飞行区域区道路表面上选取不同类型的路面样本(包含平整度质量较好与较差，覆盖健康与病害类型)进行取样拍照并作特征描述与分类统计，所有的数据采集工作都进行位置记录与拍照保存。本研究总共采集了10个关键地物特征点，9种不同类型的路面样本(样本点总数52，每种路面样本大约为5-7个)。路面样本类型包括：坑槽(分大块状和小块状)、塌陷(分石料完全暴露和石料未暴露)、裂缝(大裂缝和小裂缝)、麻面、沙石覆盖、路面平整度一般及平整度良好，如图7所示。
(2)数据处理
数据的具体处理步骤如图1所示，具体的步骤执行过程如下。
第一步：由rigelvux-1lr激光扫描仪所获取的lidar点云数据包含的属性有：三维坐标、扫描角、回波数量、回波次数、激光强度。原始点云文件并非通用的点云文件格式，借助rigel公司提供的配套软件对点云数据进行格式转换得到通用的.las格式，如图8所示；
第二步：对上述lidar点云数据进行去噪处理，实验过程中先通过目视识别初步将噪声点进行剔除，然后设置搜索半径为0.2m对局部高程异常点进行剔除，接着根据道路宽度和无人机飞行高度计算得到的路面点云的扫描角范围在±8°内，考虑到航线有一定的偏差，实验过程中设置扫描角阈值为±24°对点云数据进行过滤；
第三步：利用瑞士pix4d公司开发的商用无人机数据处理软件pix4dmapper对待测路段内无人机拍摄的高分辨率影像数据进行影像拼接和配准处理后可得到待测路段的正射影像图(图9)。然后对该正射影像进行特征提取，计算红绿蓝三波段的灰度值标准差作为新的特征波段，并与红绿蓝三波段组合，得到含有光谱信息的4个特征波段；
第四步：采用仿射变换模型将点云中的点与影像中的像素进行匹配，然后将影像上的波段值融合到点云数据中，结果如图10所示(以蓝波段信息为例)；
第五步：对点云数据与高分辨率影像数据进行特征融合后，需要对该融合数据进行多尺度几何特征参数计算。由于点云数据量非常大，为节省时间，同时又不影响结果的可靠性，本例先对点云数据的存储结构进行优化，即构建八叉树存储结构，设置一个比较适宜的数值(本例中使用的是0.3m)作为空间存储单元长度进行重采样，得到抽稀后的点云数据；
第六步：以上述点云数据的每个点为中心，设置搜索半径在[0.2m,1.0m]内且以0.1m为间隔单位进行搜索，搜索的邻域点则是去噪后的点云数据中与当前搜索点距离小于等于搜索半径的点集，然后分别计算当前尺度下的局部粗糙度(ldr)、局部维度特征(ldf)、局部高程差(lhd)。随着空间搜索尺度的调整，即可得到不同尺度几何特征信息的点云数据；
第七步：对步骤六得到的点云数据通过目视选择出覆盖所有地物类型的部分点云作为训练样本，然后通过特征选择对特征变量进行降维(如图11)，容易看出特征选择后多尺度几何特征参数lhd相对于其他特征变量而言对各类地物的区分能力最强。进一步借助随机森林模型进行分类器构建。本例在构建过程中设置树的个数为200，树的深度为40。其次利用所构建的随机森林模型对将抽稀后的点云数据进行分类，从而完成抽稀点云的分类过程；
第八步：由于到目前为止只完成了对抽稀的点云进行了分类，需要采取一定的方法对去噪后的全部点云进行分类。由于空间距离越近的两个物体其属性越相似，因此本例按照空间最近邻原则给去噪后的点云数据的类别属性赋以距离该点最近的已分类点的类别属性值；
第九步：对分类后的点云数据进行类别筛选，得到道路点云数据，再通过不规则三角网构建道路表面tin模型(图12)，进一步利用自然邻域插值得到道路表面的dsm模型，为了满足iri指数计算需要同时兼顾数据处理需要，栅格分辨率设为50mm，结果如图13；
第十步：在arcgis10中通过交互矢量化方式，在道路内部沿道路方向数字化生成间隔为0.125m的一系列纵剖面线，然后利用iri指数模型计算这些剖面线上的iri值。为了反映各路段纵向上的差别，在计算iri时需设置合适的计算单元长度，本实例设置为10m。为了可视化效果再将计算的iri值进行栅格化得到整个路段的平整度空间分布结果，如图14；
第十一步：依据本方法设定的平整度质量分级设置标准，对得到的路面平整度空间分布图进行重分类，得到路面的平整度质量分布图，如图15所示，通过对不同质量等级的栅格单元进行统计分析，可确定该路段的各质量等级分布比例统计图，如图16所示；
(3)结果分析
总体上来看，纵向方面的特点为道路南段的质量评价结果主要以“fair”为主，少部分区域是“failed”；道路中部的质量评价最差，评价结果主要为“poor”或“failed”；而道路北段的质量评价结果则以“fair”或“failed”居多，这与平整度的空间分布趋势一致。横向方面的特点则是中心线附近的质量等级要比两侧车行道低，以“failed”为主，而行车道附近的评价结果则主要以“fair”或“poor”为主，道路边界地带的路面质量等级则相对较低，为“poor”或“failed”，这部分区域由于路面沉陷、砂土混合等原因导致表面粗糙不平。结合各质量等级分布图来看(图15)，可知该路段整体的平整度质量以“fair”为主，其分布比例达41.1％，质量为“poor”或“failed”的路段分布比例居其次，各为27.4％和24.9％，而质量为“good”的路段所占比例最小，仅6.6％。而且质量为“poor”或者“failed”的分布比例总和占整个路段的一半以上，这就表明该路段整体的质量偏差，需要做出重大维修。
(4)精度评价
以实地采集的52个路面样本为分析数据，通过选取外观、表面粗糙度、车辙深度、病害面积、损坏程度这5项指标使用专家评分法对其进行平整度质量评价，并将其作为参考值。然后对这些路面样本依据实地测量的gps点坐标并辅助以高分辨率影像数据在dsm模型上进行地图定位，接着根据路面平整度空间质量空间分布图对这些路面样本的平整度质量等级进行识别，再对所有路面样本的识别结果进行统计，并由此得到用于进行精度评定的混淆矩阵(如图17)。最后利用式21和式22分别计算得到了本方法的总体分类精度为75％，kappa系数为0.65，表明该方法的整体评价结果与实际结果大体相同，且一致性程度较高。
需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。
- 还没有人留言评论。精彩留言会获得点赞！