本发明公开了一种基于地形采样点位置的海底地形数据建立方法，属于地理信息系统空间领域。所述方法首先从电子海图中提取海底地形采样点数据，对海底地形采样数据进行四叉树空间索引，查找插值点所在四叉树空间索引中位置；根据插值点所处的位置得到合适的插值方法；最后对已经进行插值的次数和需要插值的总次数进行比较，如果小于需要插值的总次数，则继续插值过程，否则插值过程结束。本发明根据海底地形采样数据边界附近插值点无法获取连成环形自然邻点的特点，采用了自然邻点和加权平均法插值混合插值的方式，充分考虑到了插值方法的特点，有效减少了地形插值误差。
1.基于地形采样点位置的海底地形数据建立方法，其特征在于：包括如下步骤：步骤一：从电子海图中提取海底地形采样点数据；
步骤2.1、确定插值区域的范围，遍历所有海底地形采样点数据，依次比较它们的坐标，得到地形采样点的范围，即xmin,xmax,ymin,ymax，设地形采样点坐标为(xi,yi,zi)，i＝1,2,3…n，n为采样点个数，那么xmin＝min(x1,x2,x3,…,xn)，xmax＝max(x1,x2,x3,…,xn)，ymin＝min(y1,y2,y3,…,yn)，ymax＝max(y1,y2,y3,…,yn)，建立根节点，将插值区域内所有地形采样点移入根节点，若采样点个数n≤NP，则执行步骤2.6结束；否则执行步骤2.2；NP表示叶子节点内海底地形采样点个数的阈值；
步骤2.2、若节点内的海底地形采样数据n＞NP,将节点的最小外包矩形进行4等分，生成4个子区，每个子区作为一个孩子节点，共四个孩子节点Si(i＝1,2,3,4)；
步骤2.3、依次将步骤2.2中被划分的节点中的海底地形采样点划分的四个范围分别移入它的四个孩子节点中，并计算落入各孩子节点Si内的地形采样点个数ni，并保存它的深度LEVEL，最后将父节点作为中间节点存储；
步骤2.4、依次检查四个孩子节点Si内的地形采样点个数ni,若ni＞NP，则递归执行步骤2.2～步骤2.4；若ni≤NP，执行步骤2.5；
步骤2.5、若某个孩子节点满足ni≤NP，将该孩子节点作为叶子节点存储；
步骤三：查找插值点所在四叉树空间索引中位置，具体为：从四叉树的根节点开始，孩子节点被分成NE、NW、SE、SW四个节点，依次判断插值点在哪个孩子节点的MBRi内，若该孩子节点为中间节点，则继续判断下去，一直找到叶子节点，即得到该插值点位于整个海底地形采样区域中的位置；
步骤四：根据插值点所处的位置得到合适的插值方法；如果插值点位于海底地形采样区域的边界附近，即所在的叶子节点的最小外包矩形和整个海底地形采样区域边界相交，则采用加权平均法；如果插值点不位于边界点附近，则采用自然邻点法插值；
加权平均法插值采用搜索圆来寻找附近的数据，从而求得插值点的高程值z；公式如下：其中，n表示搜索圆内海底地形采样点的个数，pi为第i个海底地形采样点的权重，zi为第i个海底地形采样点的高程值；
在这里我们使用快速寻找自然邻点的方法，以插值点为中心，取边长为2r的正方形虚框内的节点作为插值点X的初始影响区域；首先找到离插值点X最近的节点P1；以插值点X和节点P1的连线X-P1作为初始边，在X-P1左侧的所有节点中按照最大角准则(Delaunay准则)找到一点,使 和 夹角最大，按照上述的方法，找到其他节点；然后使用Laplace插值函数；设X为待插值点，它的测试函数可以定义为： 其中：N为待插值点X的自然邻点的数目；φi(X)为插值点X第i个自然邻点的形函数，它表示了第i个自然邻点所占的权重， 表示表示第i个自然邻点对插值点X的影响权值，u(Pi)为第i个自然邻点Pi的属性值，u(X)为插值点X的属性估计值；Laplace插值形函数定义为：式中：φi(X)满足： si(X)表示X与节点Pi相关联的Voronoi边的长度，表示区间的 测度，用 表示；hi(X)表示X到自然邻点Pi的距离，用d(X,Pi)表示；
步骤五：已经进行插值的次数和需要插值的总次数进行比较；如果小于需要插值的总次数，则继续插值过程，否则插值过程结束。
2.根据权利要求1所述的基于地形采样点位置的海底地形数据建立方法，其特征在于：步骤一具体为：步骤1.1、读取电子海图头文件；
步骤1.2、从电子海图要素数据层中读出水深层信息，水深层信息以海底地形数据的形式存储，包含位置信息和高程信息；
步骤1.3、坐标系的转换：从电子海图中读出的位置信息是相对于海图左下角的建立的海图坐标，然后将海图坐标转换成地理坐标，将地理坐标转换成经纬度；
步骤1.4、将步骤1.2读出的水深层信息存储在文本文件中。
3.根据权利要求1所述的基于地形采样点位置的海底地形数据建立方法，其特征在于：阈值NP的值取经验值NP＝5。
4.根据权利要求1所述的基于地形采样点位置的海底地形数据建立方法，其特征在于：所述的四叉树节点数据结构如下:(1)、根节点：(＜n,MBR＞)，其中，n为整个海底地形采样点的个数，MBR为整个海底地形区域的范围，包括上下左右边界范围；
(2)、孩子节点：(＜ni,MBRi＞)，其中，ni为孩子节点对应范围内的海底地形采样点的个数，MBRi为该孩子节点中地形采样点的最小包围矩形，其中，i＝1,2,3,4；
(LEVEL,＜S1,MBR1＞,＜S2,MBR2＞,＜S3,MBR3＞,＜S4,MBR4＞)，其中，LEVEL为四叉树节点层次值，＜Si,MBRi＞(i＝1,2,3,4)为索引项，Si为分别指向四个孩子节点的指针，MBRi为中间节点各个子节点的最小包围矩形；
(4)、叶子节点：(LEVEL,COUNT＜T1,MBR1＞,＜T2,MBR2＞,…＜Tn,MBRn＞)其中，LEVEL为四叉树节点层次值，COUNT为叶子节点中海底地形采样点的数目，n表示叶子节点的个数；＜Ti,MBRi＞为地形采样点数据信息，Ti为每个采样点的经纬度信息和水深值信息，MBRi为该地形采样点的最小包围矩形，其中i＝1,2,…n。
5.根据权利要求1所述的基于地形采样点位置的海底地形数据建立方法，其特征在于：所述的加权平均法插值采用搜索圆来寻找附近的数据，从而求得插值点的高程值z，公式如下：其中，n表示搜索圆内海底地形采样点的个数，pi为第i个海底地形采样点的权重，zi为第i个海底地形采样点的高程值。
基于地形采样点位置的海底地形数据建立方法
技术领域
[0001] 本发明属于地理信息系统空间领域，涉及海底地形的拟合或模拟，特别涉及一种基于地形采样点位置的海底地形数据建立方法。
背景技术
[0002] 数字高程模型是描述地表起伏形态特征的空间数据模型，是由地面规则格网点的高程值构成的矩阵，形成栅格结构数据集。数字高程模型建立过程的关键环节是格网点上高程的计算，在数学上属于数值分析中的插值问题。任意一种插值方法都是基于原始地形起伏变化的连续性即邻近的数据点间的相关性来求解待定点的高程。地形内插方法还没有统一的分类标准，从数据分布规律来讲地形插值有基于规则分布的插值方法、基于不规则分布的插值方法；按插值点分布范围，分为整体内插、局部内插和逐点内插。
[0003] 整体内插是在整个区域用一个数学函数来表达地形曲面，要求采样点的个数大于或等于多项式的待定系数，缺点是不易得到稳定的数值解，且求解速度较慢。
[0004] 局部内插是将区域按一定方法进行分块，对每块区域根据地形曲面特征单独进行曲面拟合，分块插值简化了地形曲面的形态，使得每块都可用不同曲面进行表达，缺点是如何进行分块和不能保证各块曲面之间的连续性。局部内插方法有最小二乘配置法、克里金法、样条函数法、多层曲面叠加法，局部插值方法有线性插值和双线性插值。
[0005] 逐点插值是以插值点为中心，确定一个邻域范围，用落在邻域范围内的采样点计算内插点的高程。逐点内插与局部内插的区别是，局部内插的分块范围一经确定，凡是落在该范围内的点都要进行曲面拟合，而逐点内插法邻域的范围大小、形状、采样点的个数都随内插点变动，一套数据只能用来进行一个内插点的计算。逐点内插法由于内插效率较高而成为目前数字高程模型生产中常采用的方法。逐点插值包含了加权平均法、自然邻点插值法等。加权平均法采用加权平均插值法思想，可以使用确切的或者圆滑的方式插值，它简单易行，可以为变量值变化很大的数据集提供一个合理的插值结果，不会产生无意义的插值结果。Braun在Nature上发表的文章《A numerical method for solving partial differential equations on highly irregular evolving grids》和Sukumer发表的文章《Natural neighbour Galerkin methods》分别在1995年或者1998年利用自然领域插值法构造形函数，提出了自然单元法，这是一种真正的无网格法，它采用了基于Voronoi图的Sibson或者Laplace插值方法全域构造函数。自然单元法求解精度高、计算量小、可准确施加边界条件。Sibson函数需要构造二阶Voronoi图，而Laplace插值利用Voronoi边的距离构造插值基函数，使得Laplace插值计算比Sibson插值要大为简化，在三维空间尤为突出。Sibson插值在凸区域的边界是线性精确地，但是对凹区域的边界，插值是不精确的。然而Laplace插值没有这个限制，因此我们一般采用Laplace插值。根据海底地形的特点，整体内海底地形分布联系并不很大，但是在局部内，水下地形联系紧密，单独个体的水深点受到周围邻近水深点的影响比较大。样条函数法适合非常光滑的表面，适合于根据很密的点内插等值线，然而它难以对误差进行估计，计算量也较大。Kriging法是一个误差很小的插值方法，它考虑到了空间采样点所提供的各种信息，然而它计算量很大，插值速度慢，实时性不好，所以我们选择Laplace插值用于海底地形的插值。
[0006] 自然邻点是在Voronoi图中与这个节点拥有公共边的其他节点。使用Laplace插值，最关键的是求出自然邻接点，然后计算出自然邻点的形函数（即得到该点所占的权重），最终求得插值点的高程值，在计算中，我们要求自然邻点构成环形，这样才能计算出形函数。田峰敏的《基于先验地形数据处理的水下潜器地形辅助导航方法研究》论文（哈尔滨工程大学硕士论文2007年1月）中使用了海底地形采样点的Laplace插值，具有插值域稳定、形函数满足kroneker条件、高平滑度、满足分解及线性连续条件等优点，但是它只在整个地形采样点的中间部分进行插值，没有考虑到内插点在边界附近的插值。然而对离散点进行Delaunay三角网剖分时候，位于凸壳上的点的自然邻点是不能构成环形的，这就是说位于离散点边界附近的插差点通过Delaunay三角网剖分得到的自然邻点就不能构成环形，这样就会导致无法用Laplace插值。
发明内容
[0007] 针对现有的Laplace插值在海底地形采样点区域边界处无法插值，导致插值后形成的地形的边界处附近地形不光滑，影响了生成的地形的连续性、真实性。发明提出一种针对从电子海图中得到海底地形采样数据边界处无法获取连成环形自然邻点，从而导致无法采用Laplace插值的特点，采用了先用四叉树(quad-tree)对海底地形采样区域进行分块索引，近似得到海底地形采样点边界点，从而选择合适的插值方法。如果插值点位于海底地形采样点区域边界附近，我们采用加权平均法插值，如果插值点不位于边界点附近，我们采用自然邻点法插值。整个大区域地形采用速度较快，误差较小的Laplace插值，而边界处使用加权平均法去弥补其在海底地形边界处无法插值的弊端，因此混合插值使得位于整个海底地形区域的插值点都找到合适的插值方法。采用混合插值的方式充分考虑到了实际插值中遇到的问题，充分考虑到了自然邻点插值和加权平均法插值的优势，得到了一种快速性好，误差较低的插值方法，应用该插值方法生成的地形数据具有边界处连续性好、真实自然的特点。
[0008] 本发明公开了一种基于地形采样点所在位置的海底地形混合插值方法，主要包括了以下步骤：
[0009] 步骤一：从电子海图中提取海底地形采样点数据；
[0010] 步骤二：对海底地形采样数据进行四叉树空间索引；
[0011] 步骤三：查找插值点所在四叉树空间索引中位置；
[0012] 步骤四：根据插值点所处的位置得到合适的插值方法；
[0013] 步骤五：已经进行插值的次数和需要插值的总次数进行比较。如果小于需要插值的总次数，则继续插值过程，否则插值过程结束。
[0014] 本发明的优点在于：
[0015] (1)本发明采用四叉树索引的方式，得到海底地形采样点所在位置，近似得到位于海底地形采样点边界处的点，为我们选择合适的插值提供了一种方法。
[0016] (2)本发明根据海底地形采样数据边界附近插值点无法获取连成环形自然邻点的特点，采用了自然邻点和加权平均法插值混合插值的方式。当插值点位于海底地形采样点边界附近，我们采用加权平均法，如果插值点不位于边界附近，我们采用自然邻点法插值。这样充分考虑到了插值方法的特点，有效减少了地形插值误差。
附图说明
[0017] 图1：本发明提出的基于海底地形采样点位置的海底地形混合插值方法流程图；
[0018] 图2：本发明提出的海底地形采样点四叉树索引后近似得到海底地形采样区域边界点方法示意图；
[0019] 图3：本发明设计的插值点的快速定位方法示意图。
具体实施方式
[0020] 下面将结合附图对本发明进行详细说明。
[0021] 本发明提出一种基于海底地形采样点位置的海底地形混合插值方法，如图1所示流程，具体包括了以下几个步骤：
[0022] 步骤一：从电子海图中提取海底地形数据。
[0023] 从电子海图中读出海底地形数据的方法如下：
[0024] 步骤1.1、读取电子海图头文件，电子海图头文件包含了整个电子海图的基本信息，包含了图号、比例尺、投影方式、深度单位、索引区首地址、要素数目、数据区首地址、数据区大小等基本信息。
[0025] 步骤1.2、从电子海图要素数据层中读出水深层信息，水深层信息以海底地形数据的形式存储，包含位置信息和高程信息。
[0026] 步骤1.3、坐标系的转换。从电子海图中读出的位置信息是相对于海图左下角的建立的海图坐标（单位：0.01毫米），然后转换成地理坐标（单位：米），将地理坐标转换成经纬度（单位：弧度）。
[0027] 步骤1.4、将步骤2读出的水深层信息存储在文本文件中，以便使用时读取。
[0028] 步骤二：对海底地形采样点数据进行四叉树空间索引。
[0029] 本发明中首先定义四叉树节点数据结构如下:
[0030] (1)、根节点：(＜n,MBR＞)，其中，n为整个海底地形采样点的个数，MBR为整个海底地形区域的范围，包括上下左右边界范围。
[0031] (2)、孩子节点：(＜ni,MBRi＞)，其中，ni为孩子节点对应范围内的海底地形采样点的个数，MBRi为该孩子节点中地形采样点的最小包围矩形，其中，i＝1,2,3,4。
[0032] (3)、中间节点（子树个数不为0的节点）：
[0033] (LEVEL,＜S1,MBR1＞,＜S2,MBR2＞,＜S3,MBR3＞,＜S4,MBR4＞)，其中，LEVEL为四叉树节点层次值，＜Si,MBRi＞(i＝1,2,3,4)为索引项，Si为分别指向四个孩子节点的指针，MBRi为中间节点各个子节点的最小包围矩形。
[0034] (4)、叶子节点（子树个数为0的节点）：(LEVEL,COUNT＜T1,MBR1＞,＜T2,MBR2＞,...＜Tn,MBRn＞)其中，LEVEL为四叉树节点层次值，COUNT为叶子节点中海底地形采样点的数目，n表示叶子节点的个数。＜Ti,MBRi＞为地形采样点数据信息，Ti为每个采样点的经纬度信息和水深值信息，MBRi为该地形采样点的最小包围矩形，其中i＝1,2,...n。
[0035] 如图2所示，建立四叉树空间索引的方法如下：
[0036] 步骤2.1、确定插值区域的范围，遍历所有海底地形采样点数据，依次比较它们的坐标，得到地形采样点的范围，即xmin,xmax,ymin,ymax，设地形采样点坐标为(xi,yi,zi)，i＝1,2,3...n，n为采样点个数，那么xmin=min(x1,x2,x3,…,xn)，xmax=max(x1,x2,x3,…,xn)，ymin=min(y1,y2,y3,…,yn)，ymax=max(y1,y2,y3,…,yn)，建立根节点，将插值区域内所有地形采样点移入根节点，若采样点个数n≤NP，则执行步骤2.6结束；否则执行步骤2.2。
[0037] 步骤2.2、若节点内（若步骤2第一次执行，此节点为步骤2.1中的根节点）的海底地形采样数据n＞NP,将节点的最小外包矩形进行4等分，生成4个子区，每个子区作为一个孩子节点，共四个孩子节点Si(i＝1,2,3,4)；
[0038] 步骤2.3、依次将步骤2.2中被划分的节点中的海底地形采样点划分的四个范围分别移入它的四个孩子节点中，并计算落入各孩子节点Si内的地形采样点个数ni，并保存它的深度LEVEL，最后将父节点作为中间节点存储；
[0039] 步骤2.4、依次检查四个孩子节点Si内的地形采样点个数ni,若ni＞NP，则递归执行步骤2.2~步骤2.3。
[0040] 步骤2.5、若某个孩子节点满足ni≤NP，将该孩子节点作为叶子节点存储。
[0041] 步骤2.6、结束。
[0042] NP表示叶子节点内海底地形采样点个数的阈值。合理选择阈值NP的值来平衡四叉树的构建时间和检索时间，NP越大，四叉树的深度越小，构建的索引时间越短，但是节点内搜索地形采样点的速率降低；反之NP越小，四叉树的深度越大，索引构建时间加长并增加存储空间，但节点内的搜索地形采样点的搜索效率提高，一般下我们取经验值NP＝5。
[0043] 步骤三：查找插值点所在四叉树空间索引中位置。
[0044] 如图3所示，查找插值点所在四叉树空间索引中位置的方法，具体如下：
[0045] 从四叉树的根节点开始，孩子节点被分成NE、NW、SE、SW四个节点，依次判断插值点在哪个孩子节点的MBRi内，若该孩子节点为中间节点，则继续判断下去，一直找到叶子节点，即得到该插值点位于整个海底地形采样区域中的位置。
[0046] 步骤四：根据插值点所处的位置得到合适的插值方法。
[0047] 根据插值点所在叶子节点的MBRi，我们便得到该插值点在整个海底地形采样区域的位置。如果插值点位于海底地形采样区域的边界附近，即所在的叶子节点的最小外包矩形和整个海底地形采样区域边界相交，我们采用加权平均法；如果插值点不位于边界点附近，我们采用自然邻点法插值。整个插值方法充分考虑到实际插值中遇到的问题，考虑到了自然邻点插值和加权平均法各自的优点，实现了速度较快误差较小的插值。将进行插值的插值次数加1。
[0048] 步骤4.1加权平均法插值方法；
[0049] 加权平均法插值采用搜索圆来寻找附近的数据，从而求得插值点的高程值z。公式如下：
[0050]
[0051] 其中，n表示搜索圆内海底地形采样点的个数，pi为第i个海底地形采样点的权重，zi为第i个海底地形采样点的高程值。
[0052] 步骤4.2自然邻点插值方法；
[0053] 在这里我们使用快速寻找自然邻点的方法，以插值点为中心，取边长为2r的正方形虚框内的节点作为插值点X的初始影响区域。首先找到离插值点X最近的节点P1；以插值点X和节点P1的连线X-P1作为初始边，在X-P1左侧的所有节点中按照最大角准则（Delaunay准则）找到一点,使 和 夹角最大，按照上述的方法，找到其他节点。然后使用Laplace插值函数。设X为待插值点，它的测试函数可以定义为：其中：N为待插
[0054] 值点X的自然邻点的数目；φi(X)为插值点X第i个自然邻点的形函数，它表示了第i个自然邻点所占的权重， 表示表示第i个自然邻点对插值点X的影响权值，u(Pi)为第i个自然邻点Pi的属性值，u(X)为插值点X的属性估计值；Laplace插值形函数定义为：
[0055]
[0056] 式中：φi(X)满足： si(X)表示X与节点Pi相关联的Voronoi边的长度，表示区间的 测度，用 表示；hi(X)表示X到自然邻点Pi的距离，用
d(X,Pi)表示
[0057] 步骤五：已经进行插值的次数和需要插值的总次数进行比较。如果小于需要插值的总次数，则继续插值过程，否则插值过程结束。
[0058] 上述步骤结束后得到基于地形采样点位置的海底地形数据，利用所述建立方法得到的海底地形数据为规则的网格数据，规则网格数据便于存储和索引，这为以后的三维可视化建模、地形分析提供了方便。应用该插值方法生成的地形数据具有边界处连续性好、真实自然的特点。
