基于金字塔的海量激光雷达点云存储方法

摘要

本发明涉及测绘数据处理技术领域，尤其涉及基于金字塔的海量机载激光雷达点云存储方法。通过建立一种渐进加密金字塔式的数据存储方式，在金字塔的各层通过抽稀存储不同密度的瓦片式点云对海量点云数据进行空间索引。基于金字塔的点云数据管理方式可以实现对海量点云高效有序的存储，实际数据加工操作可根据需求自适应加载或修改相应的数据瓦片，从而大幅降低制作过程中所消耗的计算资源，并直接提高基于激光雷达点云的空间数字产品制作效率。

说明书

技术领域

本发明涉及测绘数据处理技术领域，尤其涉及基于金字塔的海量机载激光雷达点云存储 方法。

背景技术

激光雷达测量(Light Detecting and Ranging，LiDAR)技术源自1970年美国航天局(NASA) 的研发，20世纪80年代得到迅速发展，20世纪90年代前后，将其与GPS和惯性导航系统 集成应用，它可以提供精确的地球空间信息。我国应用激光雷达测量技术始于20世纪90年 代，现在广泛应用于测绘、建筑、市政、电力、军事等各项国家基础建设领域。

激光雷达测量是一种主动式探测技术，它主要分为激光发射端和接收端两个模块。发射 端产生并向地物发射一束激光脉冲，在地物上产生反射，而反射为接收器接收。通过计算发 射到被反射回的传播时间，即可计算得激光传播的精确距离(光速*时间/2)。结合激光器的 位置及姿态，即可准确地计算出每个激光束被反射点的三维坐标。激光雷达可放置于多种平 台之上。其中，安放于飞机上在飞行过程中对地面进行连续激光雷达探测的系统成为机载激 光雷达系统。激光雷达产生的空间数据通称“点云”。激光雷达点云数据以文件形式存放，称 为“点云文件”。点云文件记录的信息通常包括点的三维坐标、回波次数、反射率等属性。

随着机载激光雷达硬件系统的不断成熟，机载激光雷达扫描的密度与精度日益提高。目 前，主流商业机载激光雷达的发射脉冲频率可以达到数十万Hz。飞行高度距地面1000米的 正常作业中，地表扫描线间距可以达到亚米级。这意味着每平方公里地表测量区域记录的激 光点个数可达2000万(单次条带)，而一次包含多个条带的飞行作业往往可获取数十亿个激 光点。如果不对海量激光点云进行合理的空间索引，点云的浏览、特征提取、数字产品生成 等操作均无法正常有效进行。

现有的激光雷达点云数据处理商业软件对于海量点云的存储主要是直接将原始点云根据 二维分布范围存储为一个个大的文件。按照国际通用格式las格式存储，一套包含10亿激光 点的点云数据文件大小在数十GB数量级，而每个单独瓦片往往在数百MB。由于只是简单 地将点云根据其覆盖范围分文件存储，没有考虑具体操作对数据密度的不同需求，导致每次 请求一个局部区域的数据时，均需将该区域所在文件的所有数据读出，I/O冗余操作严重。例 如，主流商业软件在制作数字高程模型(DEM)时首先将原始点云文件全部读入，然后在滤 波、内插阶段再根据要制作DEM的分辨率作为点间距，抽样选取真正所需的激光点进行计 算。假设使用密度为d米(如1米)的原始点云生产分辨率为D米(如90米)的DEM产品， 则数据读取冗余度最大可达(D-d)/d(如89)倍。

发明内容

本发明的目的在于解决现有的机载激光雷达点云数据处理商业软件在海量机载激光雷达 点云数据存储上的缺陷，提供一种有效的基于金字塔的海量激光点云数据存储方法。

本发明的技术方案采用基于金字塔分层瓦片的渐进加密式数据存储，包括以下主要步骤：

A、金字塔建立；

B、点云数据读取；

C、点云数据写回。

进一步地，步骤A中金字塔建立细化为：

A1 根据原始点云确定整片激光点云的二维边界范围，并对所有激光点进行统一 编号；

A2 用户指定第一层(底层)单位瓦片的长度(Tile_Height)、宽度(Tile_Width)、抽 稀因子(Thin_Factor，整数)，并根据此配置计算出所需金字塔的各项结构指数：金 字塔总层数N(整数)、各层点云平均密度Dn、每个瓦片的二维边界及瓦片索引 坐标(N，i，j)。

A3 通过抽稀方法逐层生成并分配第1至第N层点云数据。常用的点云抽稀方法 包括：每n(2＜＝n＜＝点总个数)个点中采样第1个点、每n(2＜＝n＜＝点总个数)个 点中随机采样1个点、基于kd树的最小二维距离采样、基于kd树的最小三维距 离采样。对于每一层金字塔，按照本层瓦片边界将点云分配至相应的瓦片并存储 为瓦片文件；

A4 删除各层中没有分配到任何激光点的空瓦片；

进一步地，步骤A2所需金字塔的各项结构指数的配置计算方法为：

A21 结合原始点云总二维边界的长度和宽度和第一层单位瓦片的长度、 宽度，可计算出第一层金字塔瓦片的总行数M和总列数N，也可计算出 第一层金字塔每个瓦片的二维范围；瓦片索引横坐标按照从左向右顺序， 纵坐标按照从上到下顺序排列；

A22 第二层单位瓦片的长度与宽度为第一层的Thin_Factor倍。以该层单 位瓦片长宽计算出该层瓦片的总行、总列数、每个瓦片的二维范围。实 际上，第二层每个瓦片的二维范围与第一层的Thin_Factor*Thin_Factor 个瓦片重合；

A23 以此方法继续生成更高层数的瓦片结构，直至整个原始数据都在同 一个瓦片内部，则金字塔最高层N确定。从底层算起，每一层的点云平 均密度为前一层的1/Thin_Factor。

进一步地，步骤B中点云数据读取方法细化为：

B1 根据数据请求的分辨率，确定一个最佳的点云金字塔层数n，即从最高层开 始第一个平均点密度高于请求分辨率的金字塔层；

B2 根据数据请求的二维范围，确定第n层金字塔中需要读取的瓦片索引范围； 进一步地，步骤C中点云数据写回方法细化为：

当第n层某瓦片(n，i，j)数据被修改后，同步修改金字塔其他层中对应瓦片数据。

其中，第n+1层到顶层的对应范围瓦片数据重新由该层数据渐进抽稀生成；第1 层到第n-1层的对应范围瓦片数据根据激光点号进行对应修改。

本发明提供了一种基于金字塔的海量机载激光雷达点云数据存储方法。该方法能够有效 地存储海量激光点云，可以实现点云数据的快速高效读入与写回，尤其适合于海量机载激光 点云数据的即时显示、快速处理、以及数字空间产品制作。如目前的主流商业软件浏览点云 时仅能加载某一固定密度的整片点云。庞大的数据量导致用户改变视角时，系统无法按照用 户的键盘鼠标输入进行同步屏幕刷新，直接影响了人机交互效率，并间接影响了数字空间产 品的生产效率。而基于本方法的点云浏览可根据用户视角高度及范围，仅加载显示最优金字 塔层中视角范围内的若干瓦片，对任意大点云数据均可达到即时交互式浏览效果。此外，本 方法也大幅提高了数字产品生产的效率。实验数据证明，使用密度为1米的20平方公里点云 数据制作10米分辨率的数字高程模型，基于传统固定数据块存储方法的生产过程中，数据读 入约需5分钟，滤波分类约需17分钟，格网内插约需4分钟，处理时间(不含人工操作)共 计26分钟；而使用基于本存储方法的生产过程中，数据读入约需0.5分钟，滤波分类约需5 分钟，格网内插约需4分钟，处理时间(不含人工操作)共计9.5分钟，效率较传统方法提 高了近3倍。

附图说明

图1是金字塔第1、2层结构示意图；

图2是某瓦片配置参数文件内容图；

图3是金字塔各层数据存储结构关系图；

图4是金字塔数据写回示意图；

具体实施方案

下文将参照附图对本发明的具体实施方案进行更详细的说明：

A 金字塔建立：

A1 循环原始点云文件，将每一个激光点赋予唯一编号，并确定整个原始点云所在区域的 二维边界范围，即点云坐标X和Y的最大和最小值(X_Min，X_Max，Y_Min，Y_Max)。

A2 用户指定第一层(最底层)单位瓦片的大小和相邻层之间的抽稀因子。单位瓦片的大 小即每个单位瓦片的高度(Tile_Height)和宽度(Tile_Width)，抽稀因子(Factor_Value)确定相 邻层单位瓦片的大小比例以及抽稀倍数。结合原始点云二维边界范围以及用户指定参数， 可计算出金字塔所有结构配置参数，即金字塔总层数、各层平均点云密度、各层瓦片索引 坐标及瓦片二维边界范围，并保存金字塔配置文件。具体步骤包括：

A21 结合原始点云总二维边界的长度和宽度和第一层单位瓦片的长度、宽度，可计算 出第一层金字塔瓦片的总行数M和总列数N，也可计算出第一层金字塔每个瓦片的二 维范围；瓦片索引横坐标按照从左向右顺序，纵坐标按照从上到下顺序排列.金字塔 各层瓦片索引坐标及瓦片二维边界范围计算公式参见表1。

表1 金字塔层各瓦片二维边界范围计算公式

(“/”表示取整，i、j表示行列坐标，k表示金字塔层数)

A22 第二层单位瓦片的长度与宽度为第一层的Thin_Factor倍。以该层单位瓦片长宽 计算出该层瓦片的总行、总列数、每个瓦片的二维范围。实际上，第二层每个瓦片的 二维范围与第一层的Thin_Factor*Thin_Factor个瓦片重合。为了便于本领域技术人员 参考理解，提供了金字塔第1、2层结构示意图(参见附图1)。

A23 以此方法继续生成更高层数的瓦片结构，直至整个原始数据都在同一个瓦片内部， 则金字塔最高层N确定。从底层算起，每一层的点云平均密度为前一层的 1/Thin_Factor。公式(1)给出了点云金字塔的总层数(Floor_Number)、原始点云的边界 范围、单位瓦片的长宽和抽稀因子这四者之间的换算关系。

$\mathrm{Floor}\_\mathrm{Number}=\left(\begin{array}{c}\mathrm{int}\left({\log }_{\mathrm{Factor}\_\mathrm{Value}}\left(\frac{X\_\mathrm{Max}-X\_\mathrm{Min}}{\mathrm{Tile}\_\mathrm{Width}}\right)\right)+1\\ \mathrm{int}\left({\log }_{\mathrm{Factor}\_\mathrm{Value}}\left(\frac{Y\_\mathrm{Max}-Y\_\mathrm{Min}}{\mathrm{Tile}\_\mathrm{Height}}\right)\right)+1\end{array}\right)---\left(1\right)$

举例，第3层的某瓦片L3(i＝3，j＝5)大小为Tile_Height、Tile_Width，则第4层的单位 瓦片大小为：Tile_Height*Factor_Value、Tile_Width*Factor_Value，第4层内点云平均间 距为第3层的1/Factor_Value。第3层该瓦片配置参数文件内容见附图2。

A3 根据A2确定的各层瓦片二维边界范围，从底层开始逐层抽稀分配原始数据至金字塔 各层各瓦片。金字塔底层数据直接来自原始数据，而第i+1(i＞0)层的数据由第i层数据 抽稀生成。具体步骤为：

A31 遍历原始点云中每个点，根据该层瓦片边界判断该点属于哪个瓦片，将点号、点 坐标及点属性存入该瓦片对应的数据文件中。第一层的行列数参见公式(2)。

$\left(\begin{array}{c}\mathrm{Row}\_\mathrm{Count}=((Y\_\mathrm{Max}-Y\_\mathrm{Min})/\mathrm{Tile}\_\mathrm{Height})+1\\ \mathrm{Column}\_\mathrm{Count}=((X\_\mathrm{Max}-X\_\mathrm{Min})/\mathrm{Tile}\_\mathrm{Width})+1\end{array}\right)---\left(2\right)$

A32 为第二层金字塔瓦片分配数据。由于第一层中已经进行过边界判断，第二层分配 点时只需判断相对应(同样边界范围)的第一层数据瓦片内中的点是否满足条件即可。 附图3给出了金字塔相邻层之间的点云抽稀关系。图中第二层L2(1，1)存储的数据由第 一层瓦片L1(1，1)，L1(1，2)，L1(2，1)，L1(2，2)中的全部点云数据抽稀而成，分配L2(1，1) 数据时只需抽稀L1中该四个瓦片即可。点云抽稀的原则可自定义，常用的抽稀方法 包括：每隔n点选取一个点、每隔n点随机选取一个点、基于kd树的最小二维距离， 基于kd数的最小三维距离等。存储第三层数据时，以第二层数据为参考，以此类推， 直至存储到最高层数据。假定Factor_Value＝2，抽稀方法为每隔一个点选取一个点存 储，则第二层金字塔中的行列数参见公式(3)。

$\left(\begin{array}{c}\mathrm{Row}\_\mathrm{Count}=((Y\_\mathrm{Max}-Y\_\mathrm{Min})/\mathrm{Tile}\_\mathrm{Height}/\mathrm{Factor}\_\mathrm{Value})+1\\ \mathrm{Column}\_\mathrm{Count}=((X\_\mathrm{Max}-X\_\mathrm{Min})/\mathrm{Tile}\_\mathrm{Width}/\mathrm{Factor}\_\mathrm{Value})+1\end{array}\right)---\left(3\right)$

A33 以同样的方法分配其余各层瓦片数据，直至金字塔顶层。

B 点云数据读取：

B1 根据数据请求的分辨率，确定一个最佳的点云金字塔层数n。该分辨率可以为要生产 的空间数据产品分辨率，如要生成5米的数字高程模型。也可根据当前用户浏览的视角高度 动态确定，如视点高度1000米对应1米点云浏览分辨率，高度5000米对应5米点云浏览分 辨率。假设请求分辨率为r米，则最佳的点云金字塔层应当为从最上层数起，第一个点云平 均密度大于r的金字塔层。

B2 根据数据请求的范围，确定第n层金字塔中需要读取的瓦片索引范围。首先计算数据 请求范围的二维边界；其次，将此边界与瓦片结构进行重叠，确定覆盖区域的瓦片，即点云 数据应当被读取返回的瓦片；再次，对于完全重叠的瓦片，读取全部数据，对于部分重叠的 瓦片，选取落于请求范围多边形内的激光点。

C 点云数据写回：

当第n层某瓦片(n，i，j)点云被删除、添加、或修改属性后，同步修改金字塔其他层中对应 瓦片数据。对于第n+1层至第N层，直接根据修改后的第n层瓦片点云进行渐进抽稀，并替 换原有瓦片点云；对于第1到第n-1层，根据第n层瓦片点云的点编号修改对应瓦片内的对 应点。

为了便于本领域技术人员参考理解，提供了点云金字塔数据写回示意图，参见附图4： 图中第一层金字塔中L1(i，j)(其中i[1，4]，j[5，8])中点云数据被删除，则第二层与第三层中 相应的数据也被删除；第二层中L2(4，1)和L2(4，2)中全部点及(3，1)和L2(3，2)中部分点的颜色 属性被改变，第一层和第三层的相应的点属性也被改变，其中Factor_Value＝2，抽稀方法为每 隔一个点选取一个点存储。