CN102521882A - 基于离散高程和自适应混合加权得到海床地形数据的方法 - Google Patents

Abstract

本发明提出一种基于离散高程和自适应混合加权的海床地形数据处理方法。将海床地形的水平面划分成均匀网格，网格顶点为控制点，控制点距水平面的高度为读入数字海图的高程数据，这些控制点的高度保存在三维数组中，根据高程点，首先使用加权内插和加权外插的方法计算基本的子点参考数据，然后使用混合加权动态校正方法动态推算出新子点，最后按照自适应LOD插值关联度表格模型进行动态调度，建立一个经过所有该LOD层所有点的曲面，生成海床地形，最后在曲面上贴上纹理就得到逼真的海床地形。

Description

基于离散高程和自适应混合加权得到海床地形数据的方法

技术领域

本发明涉及一种基于离散高程和自适应混合加权得到海床地形数据的方法，属于虚拟现实技术领域。

背景技术

大规模地形三维图像数据生成是虚拟战场技术研究中一项重要内容，目前生成地形数据的方法主要有：

1.基于分形算法创建地形，其特点是多边形数量完全可控、没有拓扑错误，但是迭代计算次数较多，地形规模很大时计算开销可观，而且难以生成真实的地形；

2.根据BMP格式的等高线图生成三维地形的方法，其优点是可生成真实地形，但是各高程间网格点的三角生成算法较复杂；

3.利用电子海图数据生成海底地形，其特点是所得地形数据一般比较粗糙，由于精度高且数据点密集的海床数字地图的获取难度较大，需要对粗粒度海床地形生成细粒度海床地形的进一步处理，但单一的加权外推算法会造成地形失真。

国内外流行的Vega三维图像开发软件中，除了基本模块外，已经开发出作为工具的若干个扩展模块，包括传感器、云彩、海洋等模块，缺少对海洋洋面以下空间环境进行开发的水下空间模块，水下模块的开发中要解决的任务之一是海床地形生成技术，即根据给定的电子海图高程离散数据生成起伏的海床地形。由于水声测量方法的局限性，一般海图存在着相邻数据间隔距离大，据此连成的曲面不平滑和生成的图像逼真度差，而已有的其它地形建立方法存在着算法复杂、计算开销大的缺陷。

发明内容

为了避免现有技术的不足之处，本发明提出一种基于离散高程和自适应混合加权得到海床地形数据的方法。

一种基于离散高程和自适应混合加权得到海床地形数据的方法，其特征在于步骤如下：

步骤1：根据数字海图数据，将欲生成的海床地形区域按数据间隔大小均匀划分为网格，其中的网格顶点为控制点P；

步骤2：将网格按地形图像数据要求，以某个网格顶点为子点，分三个过程求取该子点的高程数据；

步骤a：通过子点周边第一层所有控制点的高程数据的加权内插算法生成制点网格的x行和y列围成的区域内的某子点高度数据z_xy：

$z_{\mathrm{xy}}=\frac{\underset{i=x}{\overset{x+2}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+1}{\mathrm{\Σ}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+2,y+1)}{\mathrm{\Σ}}}{k}_{(i,j)}}+\mathrm{\ϵ},$ 其中：ε为-1～1之间的随机数

利用线性关系求出子点距离每个控制点的距离r_(i，j)，取1/r_(i，j)为加权系数k_(i，j)(i，j∈Z⁺)。

步骤b：通过子点同经度或者纬度上两层内相邻两控制点的高程数据的加权外推算法生成控制点网格的x行和y列围成的区域内的某子点高度数据z′_xy：

$z_{\mathrm{xy}}^{\′}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\′}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\′\′}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}})$

其中：

z_xy为控制点网格的x行和y列围成的区域内的某子点高度；

x、y分别为控制点网格的行、列序号，

n为控制点总行数，m为控制点总列数；

r_xy表示子点在水平面的投影点Q_xy与子点相邻右控制点在水平面的投影点Q_i，j之间的距离；

所述 $z_v^{\′}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2})$

所述 $z_v^{\′\′}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1})$

i，j均为自然数；

步骤c：通过综合步骤a和步骤b两种方法得到子点的某子点高度进行混合加权动态校正方法生成新子点的高程数据z_x，混合加权动态校正公式为：

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\′}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\′}|}{2}\mathrm{\ζ}$

其中：ζ为-1～1内随机数；

步骤3：以步骤2生成的新子点与控制点归于同一新的LOD层级，以该LOD层级重复循环步骤2，得到下一层级子点形成下一层级的LOD层级；重复循环次数小于等于17。

在步骤3完成后，根据视距、插值点数及LOD层级关系建立17级的自适应多级LOD插值关联度表格模型，利用OpenGL将步骤3得到的各子点的高程数据和原始控制点依次连接，生成三角形面片，得到由各三角形组成的海床地形曲面；

所述自适应多级LOD插值关联度表格模型为：

超宏观视距	宏观视距	微观视距	超微观视距
				10000	1000	60	3
5000	500	20	1
				4000	200	10	0.5
2000	100	5	0.3
				1～4级LOD	5～8级LOD	9～12级LOD	0.1

其中：超微观视距分为5级，精细模型特增加一层LOD，层级数对应插值点数，

变换步长为5000米、2000米、1000米、500米、300米、100米、40米、10米、5米、2米、1米、0.5米、0.2米。

在步骤2所述的步骤c中，控制点网格的x行和y列围成的区域内的某子点高度采用如下公式生成：

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\′}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\′}|}{2}\mathrm{\ζ},$ 其中：ζ为-1～1内随机数，

$z_{\mathrm{xy}}=\frac{\underset{i=x}{\overset{x+1}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+2}{\mathrm{\Σ}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+1,y+2)}{\mathrm{\Σ}}}{k}_{(i,j)}}+\mathrm{\ϵ},$ ε为-1～1之间的随机数

$z_{\mathrm{xy}}^{\′}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\′}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\′\′}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}}),$

$z_v^{\′}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2}),$ $z_v^{\′\′}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1}).$

在海床地形曲面上添加纹理，生成海床地形三维图像数据。

本发明提出的一种基于离散高程和自适应混合加权的海床地形数据处理方法。将海床地形的水平面划分成均匀网格，网格顶点为控制点，控制点距水平面的高度为读入数字海图的高程数据，这些控制点的高度保存在三维数组中，根据高程点，首先使用加权内插和加权外插的方法计算基本的子点参考数据，然后使用混合加权动态校正方法动态推算出新子点，最后按照自适应LOD插值关联度表格模型进行动态调度，建立一个经过所有该LOD层所有点的曲面，生成海床地形，最后在曲面上贴上纹理就得到逼真的海床地形。

本发明方法，可获取符合海床三维图像绘制要求的各高程数据，且计算速度快可满足实时性要求，从而解决了电子海图中高程数据过少且不连续难以生成真实地形的难题，可以根据给定的电子海图高程离散数据生成起伏的海床地形。克服了一般海图数据存在着相邻数据间隔距离大，据此连成的曲面不平滑和生成的图像逼真度差的问题。利用本方法生成的海床地形三维图像逼真、实用，而且算法简单、计算开销小。

附图说明

图1：海床地形数据处理方法流程图

图2：一维的插值示意图

图3：二维平面内部插值示意图

图4：一维加权外推模型分析

图5：二维平面外部推算插值示意图

表1：自适应多级LOD插值关联度表格

图6：4×4控制点的均匀网格控制点高程数据

图7：混合加权后新生的16×16控制点的均匀网格控制点高程数据

具体实施方式

现结合实施例、附图对本发明作进一步描述：

将海床地形的水平面划分成均匀网格，网格顶点为控制点，控制点距水平面的高度为读入数字海图的高程数据，这些控制点的高度保存在三维数组中，根据高程点，使用混合内外加权校正及自适应LOD插值的方法，建立一个经过所有点的曲面，生成海床地形，最后在曲面上贴上纹理就得到逼真的海床地形。海床地形数据处理方法流程图见图1。

由海床的自然属性可知，其地表是起伏变化的。假设起伏较平缓，则相邻控制点之间的高度变化是连续的。即子点的高度由其周围控制点的高度决定，下面分析建立及调度子点高程的模型。

方法的基本原理：在视域范围内，大地形由粗粒度模型生成细粒度模型的方法有两种：内部插值与外部推算。每一种方法都能生成地形子点，但是其不具有自校正子点能力，现对两种方法产生的子点进行混合加权来校正子点，生成符合多级梯度场的最佳子点，并根据视距通过多层LOD实现自适应分级调度。

混合加权动态推算方法的物理意义：两控制点原始真实数据的变化趋势反映在其连线上的子点，也反映在连线以外的子点，且子点偏离控制点连线越远，受其变化趋势影响越小，子点与最邻近参考点连线的方向与梯度变化方向夹角越大，受其变化趋势影响越小。

(1)内部插值方法：

一般的地形模型都是等间距的网格结构，从侧面看，对两个三维坐标数据点的插值就演变成一维的几何关系，其插值示意图如图2。这里新生成的子点P_x的经纬度信息可以从插值个数及相邻两点得到，其高度信息可以通过加权得出。随机数可以上下浮动。

(其中，分别取k₁、k₂为

ε为-1～1之间的随机数)    (1)

以上这种方法虽然理论上正确，但是失真太大，主要原因是没有考虑相邻的其他点对其影响。从二维的平面分析，其插值示意图如图3。这里新生成的子点Q_x的经纬度信息可以从插值个数及相邻两点得到，其高度信息可以通过加权得出。Z_(i，j)表示第i行第j列的高度值。则通用的计算公式如下：

$z_{\mathrm{xy}}=\frac{\underset{i=x}{\overset{x+2}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+1}{\mathrm{\Σ}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+2,y+1)}{\mathrm{\Σ}}}{k}_{(i,j)}}+\mathrm{\ϵ},$ 其中：ε位于x轴上的子点，

利用线性关系求出子点距离每个控制点的距离r_(i，j)，取1/r_(i，j)为加权系数k_(i，j)(i，j∈Z⁺)，ε为-1～1之间的随机数。

(2)外部推算方法：

设空间任意两点P₁和P₂，在基准平面(平行于水平面)的投影为Q₁和Q₂，相对于基准面的高度分别为Z₁和Z₂，如图4所示。现在推算Q₁和Q₂两点区间外任一子点Q_x的高度，它距离Q₁、Q₂点分别为r₁、r₂。首先考虑最简单的线性外推的方法，则得到图中的P_V点，但这样形成的海底轮廓线为直线，不符合真实海底的以曲面起伏波动的特点。因为子点的高度值在其相邻的控制点的高度值附近波动，且距离越近的控制点，对子点的影响越大，故可采用对子点附近的控制点高度作加权平均的方法来推算子点的高度，距离越远，权重越小。首先选取距离的平方的倒数值为加权系数，对控制点进行加权平均后，假设得到高程点为

其高度的计算公式为：

$z_x^{\′}=\frac{\frac{z_1}{r_1^2}+\frac{z_2}{r_2^2}}{\frac{1}{r_1^2}+\frac{1}{r_2^2}}---\left(3\right)$

上式整理可得：

$z_x^{\′}=\frac{z_1\·r_2^2+z_2\·r_1^2}{r_1^2+r_2^2}=\frac{z_1\·r_2^2+(z_2\·r_1^2+z_2\·r_2^2)-z_2\·r_2^2}{r_1^2+r_2^2}=\frac{z_1-z_2}{r_1^2+r_2^2}{\·r}_2^2+z_2---\left(4\right)$

可知由P₁点到P₂点的高度场梯度变化是降低的，但对上式分析可知，

点总是高过P₂点，当r₂逐渐增大时，z′_x相对于z₂的增量部分的变化趋势分析如下：

$\underset{r_2\→\∞}{\lim }\frac{z_1-z_2}{r_1^2+r_2^2}\·r_2^2=(z_1-z_2)\·\underset{r_2\→\∞}{\lim }\frac{1}{{\left(\frac{r_1}{r_2}\right)}^2+1}=(z_1-z_2)\·\underset{r_2\→\∞}{\lim }\frac{1}{{\left(\frac{r_2+\stackrel{\‾}{Q_1{Q}_2}}{r_2}\right)}^2+1}$

$=\frac{z_1-z_2}{2}---\left(5\right)$

式中

为Q₁和Q₂两点间的距离，

点的变化曲线如图4中曲线段1所示。

显然，推算公式(3)不符合控制点的梯度变化趋势，且r₂越大，误差越大。考虑到线性推算得到的P_V点高度小于P₂点高度，故可将P_V点亦加权引入到上面的推算公式中，其中加权系数选为子点与较近的控制点的距离r₂的平方的倒数，下面证明，由此式得到的子点P″_x小于P₂，符合控制点变化趋势。期望子点变化曲线如图4中曲线段2所示。修正后的推算公式为：

$z_x^{\′\′}=\frac{\frac{z_1}{r_1^2}+\frac{z_2}{r_2^2}+\frac{z_v}{r_2^2}}{\frac{1}{r_1^2}+\frac{1}{r_2^2}+\frac{1}{r_2^2}}---\left(6\right)$

同理，当Z₁＜Z₂时，以及沿P₂和P₁点向左推算时，类似地可以推出以上结论，修正公式总是合理的。

由

点的变化曲线和期望曲线可知，随着r₂的增加，其误差加大，这种变化趋势在修正后的(6)式中依然存在，因此可加大线性推算点P_V的权重，只取距离一次方的倒数为权重系数，加权平均公式如(7)式所示，求得的P_x点的高度z_x会比z₂更小。

$z_x=\frac{\frac{z_1}{r_1^2}+\frac{z_2}{r_2^2}+\frac{z_v}{r_2}}{\frac{1}{r_1^2}+\frac{1}{r_2^2}+\frac{1}{r_2}}---\left(7\right)$

需要说明的是，对于生成子点Q_x，y做加权外推需要做六次计算，外推参考点分别是：Q_i+1，j+1Q_i+2，j+1、Q_i+1，jQ_i+2，j、Q_i+3，j Q_i+2，j、Q_i+3，j+1Q_i+2，j+1、Q_i+2，j-1Q_i+2，j、Q_i+2，j+2Q_i+2，j+1。但是前四次加权外推对子点影响不大，因为两参考点外推方向也即梯度变化方向只能说明其生成子点在x轴方向上的变化趋势，不能表明其在y轴上的变化情况，且新生子点必须在沿y轴方向两相邻参考点Q_i+2，jQ_i+2，j+1内。所以只能从Q_i+2，j-1到Q_i+2，j、从Q_i+2，j+2到Q_i+2，j+1做两次加权外推计算。然后对两次外推做均值即可得到Q_x，y。二维平面加权的外部推算插值示意图见图5。

$z_{\mathrm{xy}}^{\′}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\′}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$ (8)

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\′\′}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}})$

其中：

z_xy为控制点网格的x行和y列围成的区域内的某子点高度，

x、y分别为控制点网格的行、列序号，

n为控制点总行数，m为控制点总列数，

r_xy表示Q_xy(子点在水平面的投影点)与Q_i+2，j+1(子点相邻右控制点在水平面的投影点)之间的距离；

$z_v^{\′}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2})$ 和 $z_v^{\′\′}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1})$ (i，j都是自然数)

(3)混合加权动态校正方法：

内部插值方法主要根据其相邻的参考点得出子点的高程，其优点是所得子点具有一定的可信度，缺点是没有考虑外侧参考点对其的影响。外部推算方法则是根据外侧到相邻参考点的梯度变化进行推算的，其子点的梯度变化符合地形变化规律，但是失去了其它相邻参考点对子点的影响。综合两种方法，可以考虑这两个原因，把外侧到参考点的梯度变化和相邻高权重参考点对子点的影响通过加权的方法，进行二次拟合，两个子点相互校正，得出最新子点满足真实地形数据特征。两种方法得到的子点做差，设值为a，然后设[-1，1]为新随机变量ζ的范围。取两个子点高程的均值，则可求出新子点的校正位置。

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\′}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\′}|}{2}\mathrm{\ζ}$ (ζ为-1～1内随机数)(9)

(4)自适应分级调度策略：

为了实现自适应功能，根据视觉效果需要对视距与地形网格内插值点个数结合多层LOD调度进行关联，常规的LOD细分方法为每层按照2ⁿ-1个点进行插值，其跳变幅化较大。此时可以按照自然数顺序使用逐级***子点对多层LOD进行预处理，LOD层级数与插值点数一一对应，共设定17级LOD，超过17级的更多插值可以使用云计算进行处理。自适应多级LOD插值关联度表格模型如表1。

视距就是指视点到地形模型的最短距离，当视距超过10000时可视化效果很差，此时的地形模型可以用一个平面代替，从10000到0米的视距可划分为超宏观、宏观、微观与超微观四种视距，每一种视距都以递减形式分别设定多级LOD层的变换步长，递减步长分别为5000米、2000米、1000米、500米、300米、100米、40米、10米、5米、2米、1米、0.5米、0.2米等。

表1自适应多级LOD插值关联度表格

由此地形模型，我们得到了较细的地形网格划分，而且应用式(10)非常容易编程计算。只要利用OpenGL的绘制三角形片的函数将相邻的网格点依次相连，便得到了起伏的海床地形。

例如从数字海图中读取4×4控制点的均匀网格控制点高程数据。并存于数组AN中，其海床地形如图6所示。再使用自适应混合加权动态推算地形生成算法，设视距为2000米，此时自动调用关联度模型，动态调度并***子点成为4级LOD地形，也即每一个单元格内每一级***一个子点，共***四个子点，数组变为16×16个控制点，保存在A_S中。其海床地形如图7所示。

$A_N=\left(\begin{array}{cccc}2.911& 2.707& 2.461& 2.240\\ 2.473& 3.333& 3.320& 3.100\\ 2.735& 2.761& 2.961& 2.387\\ 2.864& 2.651& 2.372& 2.456\end{array}\right),$ $A_S=\left(\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ A_{41}& A_{42}& A_{43}& A_{44}\end{array}\right)$

$A_{11}=\left(\begin{array}{cccc}2.911& 2.599& 2.547& 2.549\\ 2.980& 3.185& 3.030& 3.239\\ 3.190& 3.037& 3.284& 3.024\\ 3.258& 3.019& 2.819& 3.097\end{array}\right),$ $A_{12}=\left(\begin{array}{cccc}2.537& 2.707& 2.973& 2.554\\ 2.833& 2.853& 3.217& 2.807\\ 3.109& 3.263& 3.325& 3.059\\ 3.185& 2.903& 3.169& 3.309\end{array}\right)$

$A_{13}=\left(8\begin{array}{cccc}2.834& 2.698& 2.461& 2.269\\ 3.136& 3.272& 2.841& 2.708\\ 3.180& 3.037& 3.014& 2.624\\ 3.089& 3.325& 2.918& 2.499\end{array}\right),$ $A_{14}=\left(\begin{array}{cccc}2.367& 2.319& 2.198& 2.240\\ 2.354& 2.324& 2.466& 2.761\\ 2.554& 2.795& 2.771& 2.992\\ 2.961& 3.376& 3.111& 3.214\end{array}\right)$

$A_{21}=\left(\begin{array}{cccc}2.903& 2.814& 2.803& 3.083\\ 2.473& 2.527& 2.748& 2.758\\ 3.008& 2.899& 2.926& 3.188\\ 3.199& 2.979& 2.844& 2.934\end{array}\right),$ $A_{22}=\left(\begin{array}{cccc}2.800& 3.025& 3.374& 3.100\\ 2.975& 3.333& 3.510& 3.258\\ 3.181& 3.346& 3.172& 3.103\\ 3.306& 3.293& 3.418& 3.353\end{array}\right)$

$A_{23}=\left(\begin{array}{cccc}2.853& 2.686& 2.672& 3.061\\ 3.234& 3.195& 3.320& 3.150\\ 3.102& 3.121& 3.063& 3.326\\ 2.999& 2.576& 2.516& 2.830\end{array}\right),$ $A_{24}=\left(\begin{array}{cccc}3.241& 3.367& 3.174& 3.231\\ 3.081& 3.096& 3.028& 3.100\\ 3.105& 2.740& 2.578& 2.821\\ 2.833& 2.942& 2.642& 2.504\end{array}\right)$

$A_{31}=\left(\begin{array}{cccc}2.907& 2.644& 2.873& 3.219\\ 2.560& 2.663& 2.827& 2.871\\ 2.735& 2.717& 2.959& 3.296\\ 3.114& 3.060& 3.203& 2.983\end{array}\right),$ $A_{32}=\left(\begin{array}{cccc}3.442& 2.959& 3.128& 3.483\\ 3.169& 2.962& 2.810& 3.146\\ 3.105& 2.761& 2.683& 2.653\\ 3.264& 3.198& 3.209& 2.873\end{array}\right)$

$A_{33}=\left(\begin{array}{cccc}3.344& 2.933& 2.666& 2.396\\ 3.388& 2.967& 2.880& 2.602\\ 2.860& 3.044& 2.961& 2.569\\ 2.574& 2.757& 2.815& 2.737\end{array}\right),$ $A_{34}=\left(\begin{array}{cccc}2.536& 2.743& 2.748& 2.403\\ 2.391& 2.754& 2.625& 2.426\\ 2.435& 2.236& 2.572& 2.387\\ 2.449& 2.286& 2.329& 2.395\end{array}\right)$

$A_{41}=\left(\begin{array}{cccc}3.184& 3.805& 2.967& 2.738\\ 2.995& 2.776& 2.944& 2.934\\ 2.976& 2.682& 2.607& 2.614\\ 2.864& 2.752& 2.418& 2.326\end{array}\right),$ $A_{42}=\left(\begin{array}{cccc}2.796& 2.759& 3.098& 3.116\\ 2.579& 2.723& 3.060& 3.073\\ 2.487& 2.747& 3.806& 3.162\\ 2.651& 2.737& 2.895& 2.856\end{array}\right)$

$A_{43}=\left(\begin{array}{cccc}2.684& 2.400& 2.561& 2.839\\ 2.954& 2.608& 2.311& 2.712\\ 3.056& 2.720& 2.292& 2.477\\ 2.758& 2.770& 2.372& 2.227\end{array}\right),$ $A_{44}=\left(\begin{array}{cccc}2.647& 2.394& 2.240& 2.125\\ 2.564& 2.721& 2.618& 2.583\\ 2.325& 2.502& 2.484& 2.513\\ 2.158& 2.444& 2.406& 2.456\end{array}\right)$

由本实施例可见：本发明方法可获取符合海床三维图像绘制要求的各高程数据，且计算速度快可满足实时性要求，从而解决了电子海图中高程数据过少且不连续难以生成真实地形的难题，可以根据给定的电子海图高程离散数据生成起伏的海床地形。克服了一般海图数据存在着相邻数据间隔距离大，据此连成的曲面不平滑和生成的图像逼真度差的问题。利用本方法生成的海床地形三维图像逼真、实用，而且算法简单、计算开销小。

Claims (4)

1.一种基于离散高程和自适应混合加权得到海床地形数据的方法，其特征在于步骤如下：

步骤1：根据数字海图数据，将欲生成的海床地形区域按数据间隔大小均匀划分为网格，其中的网格顶点为控制点P；

步骤2：将网格按地形图像数据要求，以某个网格顶点为子点，分三个过程求取该子点的高程数据；

步骤a：通过子点周边第一层所有控制点的高程数据的加权内插算法生成制点网格的x行和y列围成的区域内的某子点高度数据z_xy：

$z_{\mathrm{xy}}=\frac{\underset{i=x}{\overset{x+2}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+1}{\mathrm{\Σ}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+2,y+1)}{\mathrm{\Σ}}}{k}_{(i,j)}}+\mathrm{\ϵ},$ 其中：ε为-1～1之间的随机数

利用线性关系求出子点距离每个控制点的距离r_(i，j)，取1/r_(i，j)为加权系数k_(i，j)(i，j∈Z⁺)。

步骤b：通过子点同经度或者纬度上两层内相邻两控制点的高程数据的加权外推算法生成控制点网格的x行和y列围成的区域内的某子点高度数据z′_xy：

$z_{\mathrm{xy}}^{\′}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\′}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\′\′}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}})$

其中：

z_xy为控制点网格的x行和y列围成的区域内的某子点高度；

x、y分别为控制点网格的行、列序号，

n为控制点总行数，m为控制点总列数；

r_xy表示子点在水平面的投影点Q_xy与子点相邻右控制点在水平面的投影点Q_i，j之间的距离；

所述 $z_v^{\′}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2})$

所述 $z_v^{\′\′}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1})$

i，j均为自然数；

步骤c：通过综合步骤a和步骤b两种方法得到子点的某子点高度进行混合加权动态校正方法生成新子点的高程数据z_x，混合加权动态校正公式为：

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\′}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\′}|}{2}\mathrm{\ζ}$

其中：ζ为-1～1内随机数；

步骤3：以步骤2生成的新子点与控制点归于同一新的LOD层级，以该LOD层级重复循环步骤2，得到下一层级子点形成下一层级的LOD层级；重复循环次数小于等于17。

2.根据权利要求1所述的基于离散高程和自适应混合加权的海床地形数据处理方法，其特征在于：在步骤3完成后，根据视距、插值点数及LOD层级关系建立17级的自适应多级LOD插值关联度表格模型，利用OpenGL将步骤3得到的各子点的高程数据和原始控制点依次连接，生成三角形面片，得到由各三角形组成的海床地形曲面；

所述自适应多级LOD插值关联度表格模型为：

  超宏观视距   宏观视距   微观视距   超微观视距   10000   1000   60   3   5000   500   20   1   4000   200   10   0.5   2000   100   5   0.3   1～4级LOD   5～8级LOD   9～12级LOD   0.1

其中：超微观视距分为5级，精细模型特增加一层LOD，层级数对应插值点数，

变换步长为5000米、2000米、1000米、500米、300米、100米、40米、10米、5米、2米、1米、0.5米、0.2米。

3.根据权利要求1所述的基于离散高程和自适应混合加权的海床地形数据处理方法，其特征在于：在步骤2所述的步骤c中，控制点网格的x行和y列围成的区域内的某子点高度采用如下公式生成：

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\′}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\′}|}{2}\mathrm{\ζ},$ 其中：ζ为-1～1内随机数，

$z_{\mathrm{xy}}=\frac{\underset{i=x}{\overset{x+1}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+2}{\mathrm{\Σ}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+1,y+2)}{\mathrm{\Σ}}}{k}_{(i,j)}}+\mathrm{\ϵ},$ ε为-1～1之间的随机数

$z_{\mathrm{xy}}^{\′}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\′}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\′\′}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}}),$

$z_v^{\′}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2}),$ $z_v^{\′\′}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1}).$

4.根据权利要求2所述的基于离散高程和自适应混合加权的海床地形数据处理方法，其特征在于：在海床地形曲面上添加纹理，生成海床地形三维图像数据。

Images