本发明属于地理信息技术领域，涉及一种gis矢量线要素数据的表达方法，尤其是一种矢量线要素多尺度bézier曲线分段拟合方法。
背景技术：
在科学和工程领域，定量观测数据的表达问题十分重要，经常一个合适的函数方程可以为问题的解决提供极大的方便。作为科学和工程相结合的测绘地理信息领域，测绘地理数据中线状地理实体诸如河流、区域边界、海岸线、等高线等，通常采用离散采样点的坐标串(即多边线)来表达。多边线的表达形式是对地理线要素的过度采样表达，是一种使用冗余的坐标点来对地理要素的曲线形态进行的模拟表达。多边线的表达形式简单直观，在数据处理简单性，以及可视化绘制便捷性上有一定的优势，但在这两个方面也存在不足。一方面多边线是一种冗余的数据存储，大量数据处理需要进行较多的冗余计算，计算资源消耗大；另一方面，多边线难以满足地图要素的无极缩放可视化的要求，线要素点串的采样精度决定了最终的可视化绘制的详细程度，不能随设备环境的改变而自适应的改变。
目前，用几何函数表达地理实体在工程实和践理论研究中得到推崇。随着以计算机辅助设计(cad)技术为基础的计算机辅助地图制图(cac)技术的普及，曲线几何造型工具(如：bézier、b样条曲线等)常常被用来表达线状地物。有研究提出了数图统一的地理实体基元表达模型，建立了14种基元模型与bézier曲线之间的严格数学关系，提出了地理实体用曲线表达的理论及相关实践应用。在gis数据表达上也有研究使用不同的几何基元对点串表达的gis线目标进行分段的局部拟合，得到线目标的分段表达模型，用于降低数据的采用频率，达到数据压缩的作用，也有利于数据的计算处理。参数曲线表达线要素在gis领域已有初步的应用和研究，但是考虑gis线要素整体的参数曲线表达方法，并顾及要素的多尺度特性问题，已有研究和相关技术方案还存在欠缺。因而，可以看到在地图制图的工程实践中，用分段的参数曲线表达线状地理要素的方法由来已久，但在gis系统中将线要素按照分段参数曲线的方式进行表达和分析的技术方法还存在不足。
技术实现要素：
针对以上现有技术问题，本发明设计了一种基于blg树结构对点串进行层次剖分，对不同层次的分段点串进行三次bézier曲线分段拟合，建立曲线的多尺度层次分段bézier曲线参数表达的方法。
本发明技术方案提供一种矢量线要素多尺度bézier曲线分段拟合方法，包括以下步骤：
步骤1，坐标点串表达的gis矢量线要素数据通过线化简二叉树结构层次化剖分，得到多尺度的点串曲线分段点，通过分段点的局部特征计算分段点处的切线方向，用于曲线拟合时各个分段bézier曲线之间几何连续性的保持；
步骤2，对于步骤1所得每一个分段中离散的、顺次衔接的点集
步骤3，三次bézier曲线的8元组参数按照多尺度的分段剖分进行组织，得到多尺度的bézier曲线参数表达模型。
而且，步骤1的实现过程如下，
步骤1.1，通过douglas-peucker算法，点串表达的线要素被递归分割成不同层次的线段存储在二叉树中，得到线化简二叉树结构，通过该结构得到曲线多层次分段点；
步骤1.2，计算连续三个分段点连线之间构成的一个弯曲结构的角平分线，计算角平分线的法向量，得到曲线局部结构的切线方向；
步骤1.3，对于曲线的首尾点的切线，通过计算首尾局部连续多个点的趋势线作为切线方向。
而且，步骤2的实现过程如下，
步骤2.1，设给定分段点串中各点对于待拟合的三次bézier曲线段的参数t＝{tj|j＝1,2,…,m}，以数据点的累积弧长与总弧长的商给出一个参数的初始值，公式如下：
其中，
步骤2.2，在步骤1.2和步骤1.3所得切线约束条件下，使用初始参数，通过最小二乘法进行三次bézier曲线的拟合计算，得到曲线初始的拟合模型如下，使得代价函数s最小，
其中，q3(tj)为参数tj对应三次bézier曲线上的点；
步骤2.3，利用每个点与拟合曲线的误差函数进行牛顿迭代计算得到不断优化的拟合模型，每个点与曲线的逼近度函数如下：
f(tj)＝[q3(tj)-pj]·[q3(tj)-pj]t.
迭代计算公式为：
其中，
当满足迭代的终止条件时停止迭代。
而且，迭代的终止条件是迭代前后两次拟合误差ε小于相应的预设阈值，针对每一段拟合曲线，使用
其中，e为拟合误差，e2为拟合误差的平方；
前后两次迭代误差之差ε求取如下：
ε＝el-el-1
其中，el为第l次迭代拟合误差，el-1为第l-1迭代拟合误差。
或者，迭代的终止条件是达到预定的迭代次数。
而且，步骤3的实现过程如下，
步骤3.1，将步骤2中计算得到的各个分段点串的三次bézier曲线的控制点表达转换成8元组的参数表达；
步骤3.2，将多层次的拟合曲线段按照线化简二叉树结构进行组织表达，树节点存储8元组参数以及曲线段的拟合精度值，通过检索不同的拟合精度的曲线段，得到不同的曲线表达结果。
本发明通过建立点串线要素数据的多层次bézier曲线表达结构，改变了传统的点串采样的冗余表达法，利用参数函数模型进行表达，一方面降低了数据存储的体积，另一方面有利于复杂的地理几何计算的模型化，计算效率高，同时也能实现曲线不同尺度，不同设备条件下的自适应可视化表达。
附图说明
图1是本发明的实施例的流程图。
图2是本发明的实施例点串多尺度分割示意图，其中图2a是点串的blg树结构douglas-peucker算法分割示意图，图2b是点串blg树结构组织示意图。
图3是本发明的实施例曲线分段点的切线计算示意图，其中图3a是连续分段点局部结构切线估算示意图，图3b是分段点串切线方向约束拟合示意图。
图4是本发明的实施例分段拟合迭代结果示意图，其中图4a是1次迭代拟合曲线结果示意图，图4b是50次迭代拟合曲线结果示意图，图4c是100次迭代拟合曲线结果示意图，图4d是500次迭代拟合曲线结果示意图。
图5是本发明的实施例多层次bézier曲线树结构表达示意图，其中图5a为曲线分段的示意图，图5b为其存储的树形结构图。
图6是本发明的实施例精度约束下曲线多尺度bézier曲线表达结果图，图6a是最大图面拟合误差小于50mm的分段bézier曲线表达结果图，图6b是最大图面拟合误差小于30mm的分段bézier曲线表达结果图，图6c是最大图面拟合误差小于20mm的分段bézier曲线表达结果图，图6d是最大图面拟合误差小于10mm的分段bézier曲线表达结果图，图6e是最大图面拟合误差小于5mm的分段bézier曲线表达结果图，图6f是最大图面拟合误差小于2mm的分段bézier曲线表达结果图。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
本发明考虑到，作为gis线状要素表达的一个可选方案，分段拟合参数函数表达的线数据模型具备一些有潜在的优势。首先，点串表达是对曲线的过度采样，具有较大的数据冗余，而参数化表达是对曲线的模型化拟合，有数据存储体积小，计算效率高的优势；其次，参数模型表达的几何实体在数值计算效率和计算精度方面的优势有利于gis中的几何分析和计算；最后，点串表达方式在可视化图形输出时受限于采样精度，参数模型表达可以根据可视化条件、尺度要求以及图形显示精度自适应的输出较优的可视化效果。将点串表达的gis线要素转换成参数函数模型化的表达是一个曲线拟合问题，目标是通过尽量少的分段函数获得尽量高精度的逼近。同时，gis中的线状要素具有多尺度特征，在建立线要素分段曲线表达模型时也需要顾及尺度效应。
如图1所示，本发明实施例提出的一种矢量线要素多尺度bézier曲线分段拟合方法，包括以下步骤：
步骤1：首先对地理信息系统(gis)中坐标点串采样表达的线要素进行多尺度的分段剖分，分别估算分段点的切线方向，包括对坐标点串表达的gis矢量线要素数据通过blg树(binarylinegeneralization-tree，线要素化简二叉树)结构层次化剖分，得到多尺度的曲线分段点，通过分段点的局部特征计算分段点处的切线方向，用于曲线拟合时各个分段bézier曲线之间几何连续性的保持；
实施例的步骤1采用以下子步骤实现：
步骤1.1：通过douglas-peucker算法，点串表达的线要素被递归分割成不同层次的线段存储在二叉树中，得到线化简二叉树结构，通过该结构得到曲线多层次分段点：如图2，通过douglas-peucker算法，点串表达的线要素被递归分割成不同层次的线段存储在二叉树中，即blg树结构，通过该结构得到曲线多层次分段点；
步骤1.2：计算连续三个分段点连线之间构成的一个弯曲结构的角平分线，计算角平分线的法向量，得到曲线局部结构的切线方向：如图3左所示分段点切线计算模型，s1,s2,s3分别是连续的三个分段点，它们控制着一定尺度下的一个地理弯曲特征，可以使用∠s1s2s3的角平分线作为曲线沿分段点s2的法向
步骤1.3：对于曲线首尾的分段点的切向则通过起始和末尾连续3个点的趋势线(直线拟合)来估算。
步骤2：对分段的坐标点串中各点进行三次bézier曲线参数估计，进行最小二乘拟合，并利用牛顿法不断迭代修正各点的曲线参数，得到分段点串的最佳三次bézier曲线拟合，包括对于步骤1中每一个分段中包含m个离散的、顺次衔接的点pj，构成点集
实施例的步骤2采用以下子步骤实现：
步骤2.1：给定分段点串中各点对于待拟合的三次bézier曲线段的参数t＝{tj|j＝1,2,…,m}，以数据点的累积长度与总长度的商给出一个参数的初始值
其中，
步骤2.2：在步骤1.2和步骤1.3中所得切线约束条件下，使用初始参数，通过最小二乘法进行三次bézier曲线的拟合计算，得到曲线的初始模型，计算模型为：
其中，q3(tj)为参数tj对应三次bézier曲线上的点。
即已知曲线的首尾端控制点p1,pm，模型含有两个待求解的参数，即求两个内控制点c1,c2(如图3)，使得代价函数s最小，也就是需要满足方程以下对c1,c2的偏导数为0
方程求解得：
其中
为了简化结果表达，令变量：
那么内控制点结果为：
步骤2.3：如图4所示迭代过程，其中图4a、图4b、图4c、图4d分别是1次，50次，100次，500次迭代拟合曲线结果示意图，修正参数t，需要计算每一个点pj与对应曲线上参数为tj的点q3(tj)的逼近度，其评价函数f(tj)为：
f(tj)＝[q3(tj)-pj]·[q3(tj)-pj]t
迭代公式为下式，l为迭代次数：
其中，
多次迭代更新参数值可以计算新的内控制点，直至拟合精度达到一定的要求或者达到预定的迭代次数即可，例如在实际的执行中控制迭代的次数，如l≤200。实施例中，迭代的终止点是迭代前后两次拟合误差的差值ε小于相应的预设阈值。针对每一段拟合曲线，使用点集
其中，e为拟合误差，e2为拟合误差的平方。
那么每一次迭代拟合过程中迭代前后两次拟合误差的差值ε求取如下：
ε＝el-el-1
其中，el为第l次迭代拟合误差，el-1为第l-1迭代拟合误差。
步骤3：将三次bézier曲线的8元组参数按照多尺度的分段剖分进行组织，得到多尺度的bézier曲线参数表达模型。
实施例的步骤3采用以下子步骤实现：
步骤3.1：将步骤2中计算得到的各个分段点串的三次bézier曲线的控制点表达转换成8元组的参数表达，将多层次的拟合曲线段按照blg树结构进行组织表达，树节点存储一个三次bézier曲线段的参数表达形式的8元组参数值(即曲线横、纵坐标的三次多项式函数各自的四个参数值)以及其拟合精度值，其示意结构如图5，图5a为曲线分段的示意图，图5b为其存储的树形结构图，例如，最初的粗略的拟合为al曲线段，al还可以进一步分段成ad和dl，拟合的中间部分，如ac段进一步分割成ab，bc段，依次递归分割直至曲线完全表达，树结构中每个节点存储的是三次bézier曲线的8元组参数以及其拟合精度值。
步骤3.2：通过检索不同的拟合精度的曲线段，得到不同的曲线表达结果，如图6，图6a、图6b、图6c、图6d、图6e、图6f分别是最大图面拟合误差小50mm，30mm，20mm，10mm，5mm，2mm的分段bézier曲线表达结果图，即分别是随着比例尺变大后的拟合结果。
具体实施时，可采用计算机软件技术实现以上流程的自动运行。
应当理解的是本说明书未进行详细阐述的部分均属于现有技术。
应当理解的是，上述实施例的描述较为详细，并不能因此而认为是对发明专利保护范围的限制，本领域的技术人员在本发明的启发下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明请求的保护范围应以所附权利要求为准。