地面激光扫描仪测量地形的方法与流程

本发明涉及及测绘科学与技术、计算机图形学及视觉技术等领域，尤其是一种基于卫星定位服务和迭代最近点算法测量地形的方法。

背景技术：
目前大比例尺地形测绘仍以全站仪采集细部为主导的全野外数字化测图模式，地面激光扫描仪(以下简称扫描仪)进行地形测量仍在研究阶段。扫描仪是测绘领域的高新技术，获取的数据由点云和影像组成，不仅记录了地物地貌的坐标数据和尺寸信息，更能自动记录其拓扑信息、纹理信息。与传统测量手段相比，地面激光扫描仪具有不用照准部、扫描作业自动化、数据记录自动化、获取的数据信息丰富等特点。虽然地面激光扫描仪单测站采集数据精度高、速度快，但是要获取完整的地形点云数据，则须采用多站扫描，并将多站点云拼接，然后转换到大地坐标系下。现行地面激光扫描仪测量地形的作业方法，首先采用全站仪布设导线、测量标靶，然后扫描地形点云，扫描标靶，基于标靶进行内业测站间拼接和坐标转换，从而得到大地坐标系下完整的点云数据。造成该作业过程复杂、繁琐的主要原因有四个。(1)标靶：布设标靶，测量标靶，扫描标靶回收标靶，内业提取标靶等一系列针对标靶的操作，使得每测站耗时估计增加5分钟。(2)全站仪：采用全站仪布设导线，然后测量标靶，在导线点上布设扫描测站，使得每测站平均增加至少5分钟。(3)对中整平：有些扫描仪作业要求对中整平，使得每测站耗时增加2分钟。(4)三脚架：采用三脚架固定仪器，测站转站时，为保护扫描仪需关机，下一站需开机并初始化，使得作业时间增加至少3分钟。上述原因造成地面激光扫描仪作业效率低，拼接后的点云精度不能保证，影响了地面激光扫描仪在地形测量、工程测量等领域的推广应用。中国专利公告号为CN102393183A的专利案公开了一种基于控制网和球标靶的地面激光扫描仪扫描作业 和拼接转换的方法，提高了测站点云的精度，但是在作业效率方面无改善。中国专利公告号为CN202994081U的专利案针对上述第一个因素(标靶)公开了一种定位定向基座，使得外业扫描不需要布设标靶，减轻了部分外业工作量，但是扫描仪仍需对中整平，且需要全站仪配合测量基座左右臂上的测点，作业效率没有明显提高。

技术实现要素：
本发明所要解决的技术问题是提供一种地面激光扫描仪测量地形的方法，作业效率高。为解决上述技术问题，本发明的技术方案是：一种地面激光扫描仪测量地形的方法，包括步骤：步骤一：将GPS天线与地面激光扫描仪同轴连接；步骤二：将地面激光扫描仪安装在车辆顶部；步骤三：外业扫描作业时，地面激光扫描仪实施360度扫描；步骤四：基于卫星定位服务综合系统采用GPS载波相位动态实时差分法同步测量地面激光扫描仪的站心大地坐标；步骤五：沿道路、河流对较大测区进行分块；步骤六：采用测站间公共地物点进行测块粗配准；步骤七：采用点到切平面的迭代最近点算法进行测块内多站点云精确拼接；步骤八：利用站心大地坐标将测块点云转换到大地坐标下；步骤九：将三维点云分割为地面点和非地面点，利用地面点生成等高线；步骤十：基于三维点云通过人机交互碎部采集和点云切片方法测制地物地貌线划图；步骤十一：外业调绘、修补测，对地物地貌线划图编辑整理得到地形测绘成果。本发明为地面激光扫描仪在大比例尺地形图测绘方面的推广应用打下基础，将促进该测绘仪器在城市测量、工程测量等领域的应用，对于古建筑保护、数字城市、逆向工程等领域也具有重要应用价值。本发明克服现有技术影响作业效率的四因素，本发明基于连续运行卫星定位服务 综合系统(CORS，ContinuousOperationalReferenceSystem)和迭代最近点算法(ICP，IterativeClosestPoint)，利用地面激光扫描仪快速测量地形。作为改进，所述步骤一中，应采用连接装置将GPS天线和地面激光扫描仪同轴连接，保证两仪器对中偏差在mm级。因两仪器中心高差较小(一般为20cm左右)，5度以内的倾斜造成的两仪器中心水平投影偏差不足2cm，高度偏差不足1cm，则扫描时不严格整平(倾斜5度范围)对测量精度影响不大。作为改进，采用小型轿车作为载体，在汽车顶部设置方便仪器安装和拆卸的支架，测量时，将仪器安装在支架上。作为改进，所述步骤三中，外业扫描作业时保证车辆静止，地面激光扫描仪实施360度扫描。测站转站时，车辆匀速低速移动，减低扫描仪受到的振动，不需要关闭扫描仪。为了提高测站之间的拼接精度，增加测站之间的联系强度，测站不易距离过远，同时要综合考虑点云密度、扫描作业时间、地物分布。一般地面激光扫描仪都具备蓝牙或wifi功能无线连接功能，可以在车内或车外遥控操作数据采集。作为改进，所述步骤四中，在能够满足RTK作业条件的区域，选择平整的地方停车扫描，同步按照《全球定位系统实时动态(RTK)测量技术规范》，基于CORS采用RTK测量站心大地坐标，因RTK测量时间和扫描时间基本相当，所以不增加测量作业时间。作为改进，所述步骤五中，当测区作业面积较大时，为了减小测站间配准的累积误差，测区应沿道路、河流划分为较小测块，按照测块分别进行点云拼接和坐标转换。作为改进，所述步骤六中，以测块为单元，采用粗拼接的方法完成整个测块内各扫描站的拼接，粗拼接利用相邻两测站间公共地物点(地物角点、尖锐特征点)计算坐标转换矩阵，计算方法采用七参数坐标变换法。由于点云不存在扭曲和缩放，所以点云坐标转换为刚体变换，缩放因子为1，其他六参数包括三个角度转换量(α、β、γ表示沿X、Y、Z轴的旋转角)、三个坐标平移量(tx、ty、tz)，点p为目标点，q为源坐标系点，转换公式如下：p＝Rq+T式1式中R表示旋转矩阵，T表示平移矩阵。设两个测站点云集合P＝{pi},Q＝{qi},i＝1,2，…,N，以式2为目标函数采用最小二乘法计算得到R和T的最优解；作为改进，所述步骤七中，点到切平面的ICP算法是Chen和Medioni等人在1992年提出的。为了解决算法效率问题，提高算法精确度，首先对点云按下列步骤预处理：1)对测站点云包围盒的空间按某初始边长均匀划分为立方块空间；2)遍历每一个方块，将块空间内的点云采用最小二乘法拟合成平面；3)若方块拟合平面的标准偏差小于阈值，则对方块内的点云执行重心化，简化为一个点，记录重心坐标和所拟合平面的法向量；4)否则，若立方块内的点个数大于某设定值，且块边长大于规定最小边长，则将该立方块空间继续均匀细分为八个小立方块，重复步骤2；5)全部方块处理完毕，产生了带有平面法向量的新点集。对新点集按照式3所示的目标函数执行点到切平面的ICP配准算法。设两测站新点集P’＝{pi,Npi}，Q’＝{qi,Nqi},i＝1,2，…,N，则目标函数为式中，R’旋转矩阵，T’为平移矩阵，D(p,n)为点p到法向量n对应平面的距离。作为改进，所述步骤八中，在测块四角和中心附近，从GPS测量的站心大地坐标中选择5至6个点，按照步骤七所述的坐标转换方法将精确拼接后的点云整体转换到大地坐标系下，利用剩余的已测大地坐标的站心点作为检查点，检查测块点云精度。作为改进，所述步骤九，采用滤波算法将三维点云分割为地面点和非地面点，对地面点采用Delaunay三角剖分，生成三角网数字地面模型，在三角网中内插等值点，追踪等值线，拟合光滑曲线，生成等高线。作为改进，所述步骤十中，参照地面激光扫描仪外业拍摄的照片，在点云三维场景中人机交互式采集地物碎部点，线面要素采集时，对点云按设定高程切片，导入到测图系统描绘线划图。作为改进，所述步骤十一中，对于内业点云不能识别的地物，进行外业调绘，采用全站仪和钢尺、测距仪结合对点云被遮挡区域进行补测；根据补测调绘结果在地形测绘软件中进一步编辑整理，得到最终地形测绘成果。本发明与现有技术相比所带来的有益效果是：(1)效率高：本发明不采用标靶、全站仪，不要求扫描仪对中整平，仪器安装在车辆上，测站转站不必关闭扫描仪，从而避免了上述四因素的影响，每测站总耗时仅3-5分钟，使得采用地面激光扫描进行大比例尺地形测量的作业效率提高了数倍；(2)精度高：由于采用了ICP匹配算法，减小了测站拼接的累积误差，又因采用测区分块和测块点云整体坐标转换的方法，进一步将误差降低。经检测，该方法获取的点云平面精度优于5cm，高程优于10cm，基于该点云测制的地形图满足1:500地形图的精度要求；本发明在效率和精度上的提升，为地面激光扫描仪在大比例尺地形测量上的推广应用奠定了基础。本发明不适用公共地物点难以选取的树木密集区。在不便行车的地方，要将仪器固定在三脚架上，测站转站需开关仪器，使得作业时间每测站增加约3分钟。附图说明图1为本发明的流程图。图2为点云包围盒空间划分示意图。图3为点云切片效果图。图4为案例测区航空摄影正射影像图及分块示意图。图5为测块1坐标转换选择的已测站心点分布示意图。图6为由切片点云到地形线划图的成图过程的效果图。图7为本专利测制的居民地区域1:500地形图。图8为本专利测制的苗圃和厂房区域1:500地形图。图9为本专利测制的农田区域1:500地形图。具体实施方式下面结合说明书附图对本发明作进一步说明。如图1所示，一种基于CORS和ICP算法的地面激光扫描仪快速测量地形的方法，包含以下步骤：步骤一：采用RieGLVZ400地面激光扫描仪，扫描仪自带有外置相机，以及连接装置使GPS天线与扫描仪(含相机)同轴连接；连接后保证两仪器对中偏差在mm级，因两仪器中心高差较小(一般为20cm左右)，5度以内的倾斜造成的两仪器中心水平投影偏差不足2cm，高度偏差不足1cm，则扫描时不严格整平(倾斜5度范围)对测量精度影响不大。步骤二：采用小型轿车作为载体，在汽车顶部设置方便仪器安装和拆卸的支架，测量时，将仪器安装在支架上。步骤三：选择如图4所示典型地带作为案例实施区，外业扫描作业时，车辆熄火静止，实施360度扫描，扫描期间车内人员不要有较大动作，以免影响点云质量；外业采集点云密度设为4～6cm(距离仪器100m处的点间距)，扫描测站间的距离不宜大于50m。步骤四：采用了广州市已有的GZCORS网络基准站，坐标系采用广州地方坐标系；在能够满足RTK作业条件的区域，选择较平整的地方停车扫描，同步采用GPS-RTK测量站心大地坐标，一测站内网络RTK测量时间约3-5分钟，与扫描时间相当。步骤五：内业处理时，如图4所示沿道路、河流、围墙将测区按线条进行分块；当测区作业面积较大时，为了减小测站间配准的累积误差，测区应沿道路、河流划分为较小测块，按照测块分别进行点云拼接和坐标转换。步骤六：各分块内测站配准，先选择一个视野开阔的测站，作为固定的基准测站，其他测站两两依次配准，粗配准选择不少于4个公共地物点计算转换矩阵，在居民地、厂区选择房屋角点、路灯顶点，在农田、水域采用高压线塔、电线杆、棚子角点等。粗拼接利用相邻两测站间公共地物点(地物角点、尖锐特征点)计算坐标转换矩阵。计算方法采用七参数坐标变换法。由于点云不存在扭曲和缩放，所以点云坐标转换为刚体变换，缩放因子为1，其他六参数包括三个角度转换量(α、β、γ表示沿X、Y、Z轴的旋转角)、三个坐标平移量(tx、ty、tz)。点p为目标点，q为源坐标系点，转换公式如下：p＝Rq+T式1式中R表示旋转矩阵，T表示平移矩阵。设两个测站点云集合P＝{pi},Q＝{qi},i＝1,2，…,N，以式2为目标函数采用最小二乘法计算得到R和T的最优解。步骤七：首先按照初始边长为1m对各测站内点云进行预处理，设定方块平面拟合标准偏差阈值为2cm，方块内最少点个数设为100，立方体块最小边设为20cm，预处理后的点云采用点到切平面的ICP算法进行测块内多站自动精确拼接。点到切平面的ICP算法是Chen和Medioni等人在1992年提出的，为了解决算法效率问题，提高算法精确度，首先对点云按下列步骤预处理：1)对测站点云包围盒的空间按某初始边长均匀划分为如图2所示的立方块空间；2)遍历每一个方块，将块空间内的点云采用最小二乘法拟合成平面；3)若方块拟合平面的标准偏差小于阈值，则对方块内的点云执行重心化，简化为一个点，记录重心坐标和所拟合平面的法向量；4)否则，若立方块内的点个数大于某设定值，且块边长大于规定最小边长，则将该立方块空间继续均匀细分为八个小立方块，重复步骤2；5)全部方块处理完毕，产生了带有平面法向量的新点集。对新点集按照式3所示的目标函数执行点到切平面的ICP配准算法，设两测站新点集P’＝{pi,Npi}，Q’＝{qi,Nqi},i＝1,2，…,N，则目标函数为式中，R’旋转矩阵，T’为平移矩阵，D(p,n)为点p到法向量n对应平面的距离。步骤八：如图5所示，在某测块周边和中心共选择了5个已测大地坐标的站心点，将测块点云转换到大地坐标下；转换后内符合精度见下表1。表1内符合精度计算表单位m对点云的精度检查，选择本测块内的剩余已测站心点RTK测量的坐标和转换得到坐标进行比较，点云成果的精度计算见表2；表2点云精度检测计算表单位m步骤九：采用滤波算法将三维点云分割为地面点和非地面点，对地面点采用Delaunay三角剖分，生成三角网数字地面模型，在三角网中内插等值点，追踪等值线，拟合光滑曲线，生成等高线，利用地面点基于等值线追踪法生成等高线。步骤十：参照扫描仪所拍照片，在三维点云中，通过人机交互进行碎部采集，绘制线划图；对于线面地形要素，如图3所示对点云按设定高程切片，如图6所示基于切片点云描绘线划图。步骤十一：对于内业点云不能识别的地物，进行外业调绘。采用全站仪和钢尺、测距仪结合对点云被遮挡区域进行补测。本案例在EPS2008地形图测绘软件对地形线划图编辑整理，得到最终地形成果，图7、图8、图9为本专利测量的三测块1:500地形图。对地形测绘成果采用现行全站仪测量方式进行点位精度检查，计算见表3。对地物点的间距精度采用外业钢尺或测距仪测量方式进行检查，计算结果见表4。采用全站仪测量的高程与点云测的地形图的高程进行比较，计算结果见表5；表3地物点坐标较差统计表单位m表4地物间距较差统计表单位m表4“图号”为按照广州市地形图图幅分幅标准的图幅编号，表中检测位置由字母加数字构成，表示表房屋类型和层数，比如“A3”即A类房屋三层，w表示围墙。表5点云测制的地形图高程精度检测单位m图7、图8、图9为本发明案例的三测块1:500地形图对应图4所示块1、块2、块3，合计面积约0.15平方公里，折合图幅3幅，采用本发明以RieGLVZ400扫描仪为例，外业测量消耗1.5组日，扫描了110站，测站间平均距离约30m，外业调绘与补测1组日，内业数据处理和制图约6工日。表63测块详细情况表测块名面积(km2)主要地物类型块10.03居民地、水域块20.05苗圃、厂区块30.07苗圃、农田合计0.15与采用现行扫描仪测量方式相比，外业效率极大提高，工作强度显著降低；与传统全野外测量模式比较，内业效率略有降低，各测量模式工作量详细情况见表7。表7以3幅1:500地形图为例各测量方式工作量比较表注：上述组日为测量作业小组(约3人、主要仪器各一台)工作8小时。