梁 策 齐 华,2 卢建康 黄华平

(1.西南交通大学信息科学与技术学院,四川成都 610031；2.西南交通大学土木工程学院,四川成都 610031；3.中铁二院工程集团有限责任公司测绘分院,四川成都 610031)

道路横断面是道路走向法线方向地形变化趋势的反映,一般通过野外实测得到。但野外测量任务量大、周期长。随着机载LiDAR技术的推广应用,自动提取高精度横断面成为可能。

移动曲面拟合法具有计算方法简单、分块灵活等优点,但地面采样点较稀或分布不均匀是影响其精度的主要因素[1]。机载LiDAR测量得到的地面采样点数据具有海量、高密度和高精度等特性,若提供给用户使用的地面采样点在平面上呈格网状或接近格网状分布(见图1),对于移动曲面拟合法刚好可以扬长避短。

1 内插函数模型选择

移动曲面拟合法是以内插点为中心,确定一个邻域范围,根据落在邻域范围内的全部或部分采样点个数选择拟合函数,从而展铺成一张数学曲面计算内插点的高程值[2]。

离散状分布的地面采样点适宜用移动曲面拟合法内插高程[3]。由采样点个数n选择拟合函数,当n≥8时,函数选用z=a0+a1x+a2y+a3xy+a4x2+a5y2；当n=6或7时,舍去xy项；当n=4或5时,舍去平方项；当n=3时,仅用线性项[3-4]。格网状分布的地面采样点,常采用双线性多项式内插,即选取4个最靠近内插点的采样点,函数选用z=a0+a1x+a2y+a3xy[3]。针对呈格网状或接近格网状分布的高密度地面采样点,试验表明宜采用的函数为

当采样点个数大于待求系数个数时,式(1)中系数ai可根据误差方程式V=MX-Z,按间接平差最小二乘原理求解,有VTPV=min[5],权值取采样点到内插点距离的倒数。

2 规则格网索引LiDAR地面采样点

首先,利用规则格网对LiDAR地面采样点建立索引(见图1)[4]。本文中规则格网尺寸选用20 m ×20 m。内插时用内插点平面坐标定位格网,并从相应格网中选择参与拟合的采样点,进而完成内插计算。

图1 采样点分布及规则格网索引

3 利用分级扇形区域选择采样点

移动曲面拟合法选择采样点一般考虑三个因素。范围：即选用多大邻域范围内的采样点；点数：即选择多少采样点参加计算[2]；分布[3-4]：即采样点的分布是否包围内插点。

3.1 自适应确定邻域半径

设第一级邻域含有6个采样点；以式(2)计算邻域半径R1；设第二级邻域8个方位各增加一个点,即含有14个采样点,以式(3)计算邻域半径R2；第三级半径取1.25倍格网尺寸,即R3=25 m。

式(2)(3)[4]中,A为格网面积,n为内插点所在格网包含的采样点个数,k是邻域含有的采样点个数,n1为内插点周围共9个格网(即3×3格网)包含的采样点个数。

针对覆盖区域地面采样点总体稀稠不均的情况,本方法将根据内插点所处的局部采样点密度确定邻域范围,准确度相对会提高。

3.2 分级扇形区域

以内插点为中心,设采样点基于内插点的坐标偏移量为 Δx、Δy,建立8 方位划分(见图2),以式(4)计算方位编号q。

设采样点到内插点的距离为d,建立3个分级,以式(5)计算级别编号 p。

图2 8方位划分

3个分级和8方位构成24个分级扇形区域,[p,q]是扇形区域下标(见图3)。

图3 24个分级扇形区域

分级扇形区域收集采样点的方法：每个方位只保留一个距内插点最近的采样点；根据采样点与内插点的距离d和坐标偏移量 Δx、Δy,计算扇形区域下标,并将采样点收集到该区域。

3.3 判断采样点包围内插点

设 θi(i=0～7)用于标记各个方位是否含有采样点,含有则为真(如图4所示)。

表1 采样点空间分布及判断条件

满足表1分布条件的采样点必然能包围内插点。

3.4 选择格网并遍历

图4 采样点空间分布

格网位置是固定的,而内插点是移动的。随着内插点位置的变化,邻域可能覆盖到内插点所在格网的8个相邻格网。选择采样点的简单办法是先遍历这9个格网内的所有采样点,对属于邻域的采样点再选择参与拟合的采样点,这将增加额外的计算量。

本文采用的方法是先遍历内插点所在格网中的采样点,用分级扇形区域收集采样点；若p=0级中的采样点不包围内插点,则从8个相邻格网中选择格网并遍历,从而收集采样点。

选择格网的方法及实现：

(1)数组下标从0到7对应8个方位编号(见图2)和8个相邻格网的编号(见图5(a))；

(2)每个数组成员对应一个格网,且存储其索引下标,初始状态全为F；

(3)依据表2的对应关系,将缺少采样点的方位相应数组成员的状态改为T；

(4)数组中状态为T的成员,其对应的格网即选定的格网。

表2 方位与数组下标对应关系

例如当3、4、5、6方位中已经收集到采样点时,只需遍历0 、1 、2 、3、7 格网,(见图 5)。

3.5 选择参与拟合的采样点

用分级扇形区域收集采样点后,从其中选择采样点(即参考样点)进行拟合的步骤为：

(1)取p=0级中的采样点,若采样点包围内插点,则跳至(4)；

(2)增加p=1级中的采样点,若采样点包围内插点,则跳至(4)；

(3)增加p=2级中的采样点,若采样点不包围内插点,则退出；

(4)取到的采样点为参与拟合的采样点。

图5 选择格网

有效邻域范围R即参与拟合的采样点到内插点的最大距离。

4 道路横断面内插

4.1 内插中两种特殊情况处理

若采样点距内插点不足0.2 m,以采样点高程作为内插高程值[6]。

若矩阵不可逆或采样点无法包围内插点,则标注内插失败,该内插地形点不记入横断面。

4.2 自定义等间距内插

地形特征点距中桩的距离是未知的,为了能得到比较准确的地形特征点,内插间距需要尽量小些,用户可以根据实际地形情况自行设置内插间距l。

实现方法为：已知中桩高程及横断面的中桩和边桩平面坐标；从边桩开始,每次向中桩靠近l距离内插地形点的高程；首个内插成功的地形点作为实际边桩,记载内插成功的地形点,直到中桩为止。地形点只标记其与中桩的高差和平距(左侧为负,右侧为正)；中桩视为高差和平距均为零的地形点。

4.3 D-P算法提取地形特征点

等间距内插得到的地形点未必都是地形特征点,需要进一步提取,从而在一定高程精度约束下降低数据冗余。目前最为稳定和实用的是Douglas-Peucker(简称 D-P)法[7]。

算法思想：给定阈值 δ；连接曲线首、末两点,依次计算中间点到该线段的距离；找出最大距离点,并判断该距离是否小于 δ,若是,则舍去该曲线的所有中间点；否则,保留该点,并以该点为界,将曲线分为两部分,对这两部分重复进行上述操作[7]。以中桩为界分别对横断面左右两侧进行D-P算法处理。

5 试验

本试验数据是丘陵区1136 m ×370 m范围内,呈格网状及接近格网状分布的365 974个LiDAR地面采样点,平均密度为0.87个/m2。以呈离散状分布的534个实测地面采样点作为内插对象,以实测采样点的平面坐标为内插点。内插高程与实测高程比较后统计精度如表3所示。选择参与拟合的采样点时,分级扇形区域增加到的级别,其出现次数统计如表4所示。

表3 误差及参与拟合的平均采样点个数

表4 分级扇形区域选择采样点增加到的级别使用次数

以该区域采样点建立TIN,并提取出73个横断面,共计7 956个地形点。为了与TIN中内插到的地形点做比较,利用73个横断面中的平距换算地形点平面坐标,应用移动曲面拟合法内插,计算内插值与TIN中对应地形点的误差,统计精度如表5所示。

表5 拟合法内插值与TIN中地形点的误差统计

从表3和表5可以看出,内插精度优于0.1 m；从表4可以看出,自适应确定邻域半径,绝大部分情况下R1和R2两级邻域可以提供满足分布要求的采样点,R3可以作为补充。与TIN中提取的横断面比较(见图6),可以看出移动曲面拟合法的结果在反映地形趋势方面是一致的。经D-P算法提取地形特征点后,地形趋势得以保留,冗余数据量显著减少。

图6 横断面

6 结论

利用分级扇形区域选择参与拟合的采样点,可以在确定内插邻域中兼顾范围、采样点个数及其空间分布三方面影响精度的因素,这对于保证内插的可靠性及精度是有利的；此外,对3×3规则格网有选择性的遍历,缩小了搜索空间,从而提高了选择采样点的时间效率。本文的方法对于从呈格网状及接近格网状分布的地面采样点中提取道路横断面具有一定的参考价值。

[1]兰 燕,王明华,等.逐点内插法建立DEM的研究[J].测绘科学,2009(1)：214-216

[2]李志林,朱 庆.数字高程模型[M].武汉：武汉大学出版社,2005

[3]汤国安,刘学军,闾国年.数字高程模型及地学分析的原理与方法[M].北京：科学出版社,2005

[4]柯正谊,何建邦,池天河.数字地面模型[M].北京：中国科学技术出版社,1993

[5]武汉大学测绘学院测量平差学科组.误差理论与测量平差基础[M].武汉：武汉大学出版社,2003

[6]潘 攀,王光霞,等.DEM建模方法的研究与实践[J].测绘科学技术学报,2007(2)：57-60

[7]贾利峰,齐 华.矢量曲线压缩算法与实现[J].铁道勘察,2005(2)：20-21