CN103914618A - 一种基于电子海图的海水深度等值线生成方法 - Google Patents

Abstract

本发明属于海洋地理信息系统领域，具体涉及一种基于电子海图的海水深度等值线生成方法。根据航海作业区域选择用于水深数据提取的电子海图；从电子海图中提取水深数据，生成离散数据源；利用克里金插值法对离散数据源进行插值，生成网格水深数据；利用克里金插值算法逐一对网格G_ij(x_i,y_j)处的水深h_ij进行插值；利用关联表法生成水深等值线；等值线平滑和标注；水深等值线与电子海图叠加显示。本发明以电子海图中的水深数据为离散数据源，生成海水深度等值线，并与电子海图进行叠加显示，可以为航海作业人员提供直观、形象、可靠的海水深度信息。

Description

一种基于电子海图的海水深度等值线生成方法

技术领域

本发明属于海洋地理信息系统领域，具体涉及一种基于电子海图的海水深度等值线生成方法。

背景技术

海水深度是船舶航行中必须考虑的海洋环境要素之一，直接影响船舶的航行安全和作业任务的顺利实施。单纯的水深数据难以直观、形象地呈现海水深度的分布情况和变化规律，而如果将单一的水深数据转换为直观可视的等值线，则可以清楚地了解海水深度的空间分布，方便相关人员开展航海作业任务。因此，研究海水深度等值线生成方法，实现海水深度的连续表达和可视化显示，具有十分重要的意义。海水深度数据是生成海水深度等值线的基础，Shapefile格式电子海图是我国发行的标准电子海图，是根据海洋测绘绘制而成的数字产品，可以反映真实的海洋地理信息。电子海图中包含了大量的水深测绘数据，可以作为制作海水深度等值线的离散数据源。因此，基于电子海图数据生成海水深度等值线具有重要的意义和实用价值。

发明内容

本发明的目的在于提供一种为航海作业人员提供直观、形象、可靠的海水深度信息的基于电子海图的海水深度等值线生成方法。

本发明的目的是这样实现的：

（1）根据航海作业区域选择用于水深数据提取的电子海图：

（1.1）确定航海作业区域：

航海作业区域Ω的经纬度范围为（BLat～ELat），（BLon～ELon），水深数据提取海区Ω₁的范围为（BLat-temp～ELat+temp），（BLon-temp～ELon+temp），其中temp＞0，单位为度：

$\mathrm{temp}=\frac{\mathrm{Max}(\mathrm{fabs}(\mathrm{Blon}-\mathrm{Elon}),\mathrm{fabs}(\mathrm{Blat}-\mathrm{Elat}))}{5};$

（1.2）选择电子海图：

逐一判断海图数据库中的电子海图，如果水深数据提取海区Ω₁与当前电子海图的覆盖区域有交集，将电子海图的名字存入待提取海图列表MapList，直到所有比例尺的海图均判断完毕为止；

（2）从电子海图中提取水深数据，生成离散数据源：

（2.1）从MapList中选择一幅尚未进行水深提取的电子海图；

（2.2）打开选择的电子海图文件，找到水深层所在的位置；

（2.3）读取水深层内的离散水深数据，判断当前水深是否在海区Ω₁内，若水深所在位置在Ω₁内，则将水深信息经度、纬度和水深值存储到Db.dat文件中，重复操作直到此海图中的所有水深数据均进行判断；

（2.4）判断MapList中是否还有尚未进行水深提取的电子海图，若有则从中选择一幅电子海图并执行步骤（2.2），否则说明水深数据提取完毕，Db.dat存储的即是进行插值的离散数据源；

（3）利用克里金插值法对离散数据源进行插值，生成网格水深数据：

（3.1）将航海作业区域Ω划分为M×N个网格，其中M为沿经度方向的网格数，N为沿纬度方向的网格数，G_ij(x_i,y_j)表示第i行，第j列处的网格，x_i为网格中心的经度，y_j为网格中心的纬度，网格边长l为

$l=\mathrm{Max}(\frac{\mathrm{fabs}(\mathrm{Blon}-\mathrm{Elon})}{N},\frac{\mathrm{fabs}(\mathrm{Blat}-\mathrm{Elat})}{N}),$

保持l≤0.05度；

（3.2）利用克里金插值算法逐一对网格G_ij(x_i,y_j)处的水深h_ij进行插值：取以(x_i,y_j)为中心，边长为L的正方形区域内的水深数据为参考水深，其中L=2ωl；ω＞4，保持参数水深的个数W大于3；l为网格边长，h^k为第k个参考水深的值，

（4）利用关联表法生成水深等值线：

（4.1）建立关联表：

（4.1.1）节点链表：

用节点链表存储所有节点信息,经度、纬度和深度值：

（4.1.2）边链表：

边链表存储每条边所包含的两个节点信息和每条边所属的两个邻面信息：

（4.1.3）面链表：

面链表存储它所包含的四条边在边链表中的索引号，其数据结构设计如下：

（4.1.4）索引号：

从网格左上角开始，按照从左至右，从上至下的顺序逐个进行标记,对于处于网格区域边界处的网格边，对于位于网格区域边界的网格边，把边链表中的下侧网格面的索引号标记为-1；

（4.2）等值线追踪：

将航海作业区域Ω内的水深范围按照间隔进行划分，计算出每一条要追踪的等值线属性值，航海作业区域Ω的最大水深为H_max，最小水深为H_min，等值线属性间隔为Δh，第i条等值线的属性值H_i＝H_min+Δh×(i-1)；利用双向追踪法对于任一属性值的等值线进行追踪：

（4.2.1）建立搜索标记数组

边搜索标记数组中的元素与边链表中的边一一对应，数组中所有元素初始化为FALSE，每追踪到一个等值点，搜索标记数组中的对应元素设为TRUE；

（4.2.2）追踪分支等值线

查询边搜索标记数组，以没有被搜索过的边作为等值线起点，开始双向追踪分支等值线，建立分支等值线的等值点链表并把等值点信息存入文件：

（4.2.2.1）初始化前向交点链表ForArray和后向交点链表BackArray，分别用来存储前向追踪和后向追踪过程中所得到的等值点信息；

（4.2.2.2）扫描边链表，找到等值点，作为分支等值线的线头，把等值点存入前向交点链表ForArray中；

（4.2.2.3）以步骤（4.2.2.2）中找到的等值点为起点，进行等值线前向追踪，当追踪到达区域边界时，即通过边链表检索到的面的索引值为-1时，或者在网格内找不到等值线出口时，等值线追踪结束；

（4.2.2.4）以等值线线头为起点，进行等值线后向追踪，当追踪到达区域边界时，即通过边链表检索到的面的索引值为-1时，或者在网格内找不到等值线出口时，等值线追踪结束；

（4.2.2.5）建立分支等值线所有等值点链表，根据前向交点表ForArray和后向交点表BackArray，建立包含分支等值线所有等值点的等值点链表AllArray；

（4.2.3）遍历边搜索标记数组，判断是否已经把边链表中的所有边搜索过，如果边链表中还有边没有被搜索过，则执行步骤（4.2.2），继续搜索等值线，如果边链表中的所有边都被搜索过，则属性值的所有等值线搜索完毕；

（5）等值线平滑和标注：

采用三次B样条算法对等值线进行平滑处理；

从等值点文件中以一定间隔读取等值点数据，并将等值点数据以等值线标注结构体的形式存入等值线标记数组，根据等值线标记数组中存储的等值点坐标确定要写标注的矩形区域，判断矩形区域是否在航海作业区域Ω外，如果不在Ω外，则可以标注等值线，等值线标注结构体为：

（6）水深等值线与电子海图叠加显示：

从等值线文件中读取数值等值线上的数据点，将等值点的经纬度坐标转换为海图坐标系下的图上坐标，在电子海图上将各等值点连接成线即可得到该数值下的等值线，采用步骤（5）生成的等值线标记数组对绘制好的等值线进行标注。

本发明的有益效果在于：

本发明以电子海图中的水深数据为离散数据源，生成海水深度等值线，并与电子海图进行叠加显示，可以为航海作业人员提供直观、形象、可靠的海水深度信息。

附图说明

图1为本发明所述基于电子海图的海深等值线生成流程图。

图2为本发明所述任一属性值的所有等值线追踪流程图。

图3为本发明所述双向追踪某一分支等值线流程图。

具体实施方式

下面结合附图对本发明作进一步描述。

本发明公开了一种基于电子海图的海水深度等值线生成方法，该方法的具体内容包括根据航行区域进行电子海图选择，从电子海图中提取离散水深数据并存入水深文件，基于离散水深数据利用克里金插值算法生成网格水深数据，根据网格数据建立节点链表、边链表和面链表来表示节点、边、面之间的拓扑关系并利用关联表法追踪水深等值线，等值线平滑和标注以及将水深等值线与电子海图叠加显示几个步骤。本发明以电子海图中的水深数据为离散数据源，生成海水深度等值线，并与电子海图进行叠加显示，可以为航海作业人员提供直观、形象、可靠的海水深度信息。

本发明针对航海作业任务的需求，提出了一种基于电子海图的海水深度等值线生成方法。本发明中提供的方法区别于现有方法的显著特征在于：从电子海图中提取水深数据，生成水深离散数据源。然后利用克里金插值法进行水深插值，构建水深网格数据。在此基础上，基于关联表生成水深等值线。最后实现水深等值线与电子海图的叠加，为航海作业人员提供直观、形象、可靠的海水深度分布信息。

本发明提供一种基于电子海图的海水深度等值线生成方法，具体包括以下几个步骤：

步骤一：根据航海作业区域选择用于水深数据提取的电子海图。

Shapefile格式电子海图是以比例尺进行分级管理的，每一比例尺下均由多幅海图，且每幅电子海图都有固定的覆盖范围，每一位置可能有多幅电子海图覆盖。因此，为了提取出尽可能多的水深数据，必须找到能够覆盖航海作业区域的所有电子海图。

步骤二：从电子海图中提取水深数据，生成离散数据源。

Shapefile格式电子海图中的要素分层进行存储，水深数据是点状要素,存储在单独的水深图层中。提取水深要素时,首先遍历选定海图的所有要素,当遍历到与指定区域有交集的水深图层时,读取该水深点并将其存入相应的水深数据文件中。依次按照上述方法从所有步骤一选取的电子海图中读取水深数据并存入水深文件。

步骤三：利用克里金插值法对离散数据源进行插值，生成网格水深数据。

Shapefile格式电子海图中的水深数据是离散分布的，而等值线的生成需要网格数据，因此必须通过插值将离散水深数据转换为网格水深数据。插值方法有很多，本专利针对水深数据所表现出的随机性和结构性的地理空间分布特征，采用克里金插值法进行水深插值运算。

步骤四：利用关联表法生成水深等值线。

根据网格数据建立节点链表、边链表和面链表来表示节点、边、面之间的拓扑关系，然后按一定间隔改变属性值，依次追踪各个属性值对应的等值线，直到所有等值线均追踪完毕。

步骤五：等值线平滑和标注。

直接把等值点依次连接起来，生成的等值线不光滑，因此需要对等值线进行平滑处理。为了方便用户的查看，每条等值线的属性值必须在等线附近标注出来。

步骤六：水深等值线与电子海图叠加显示。

电子海图已经成为航海作业人员常用的工具之一，能够为航海作业信息的显示提供基础背景。以电子海图为背景实现水深等值线的叠加绘制，能够将水深等值线直观形象地呈现给航海作业人员，有助于用户了解水深的分布和变化情况。

本发明提出一种基于电子海图的海水深度等值线生成方法，具体流程如图1所示，包括以下几个步骤：

步骤一：根据航海作业任务选择用于水深数据提取的电子海图。

设航海作业区域Ω的经纬度范围为（BLat～ELat），（～ELon），为了保证水深插值的完整性，水深数据提取海区Ω₁的范围应该为（BLat-temp～ELat+temp），（BLon-temp～ELon+temp）。其中temp＞0，单位为度，计算公式如下：

$\mathrm{temp}=\frac{\mathrm{Max}(\mathrm{fabs}(\mathrm{Blon}-\mathrm{Elon}),\mathrm{fabs}(\mathrm{Blat}-\mathrm{Elat}))}{5}$

本专利采用是Shapefile格式的电子海图作为数据源。该格式电子以比例尺进行分类，不同比例尺的海图所含要素的数目不同。为了获得尽可能多的水深数据，必须从海图数据库中将所有相关的电子海图选出。海图选择的具体方法如下：逐一判断海图数据库中的电子海图，如果水深数据提取海区Ω₁与当前电子海图的覆盖区域有交集，则说明该电子海图中含有所需的水深数据，将该电子海图的名字存入待提取海图列表MapList，重复上述方法，直到所有比例尺的海图均判断完为止。MapList存储的即是所有可供水深数据提取的电子海图。

步骤二：从电子海图中提取水深数据，生成离散数据源。

水深数据存储在电子海图的水深层，该步骤需要从MapList存储的所有电子海图中提取海区Ω₁内的所有水深，并将提取的水深数据依次存入基础数据文件Db.dat。因此具体方法如下：

步骤2.1从MapList中选择一幅尚未进行水深提取的电子海图；

步骤2.2打开选择的电子海图文件，找到水深层所在的位置；

步骤2.3读取水深层内的离散水深数据，判断当前水深是否在海区Ω₁内，若该水深所在位置在Ω₁内，则将该水深信息（包括经度、纬度和水深值）存储到Db.dat文件中。重复该操作直到此海图中的所有水深数据均进行判断。

步骤2.4判断MapList中是否还有尚未进行水深提取的电子海图，若有则从中选择一幅电子海图并执行步骤2.2，否则说明水深数据提取完毕，Db.dat存储的即是进行插值的离散数据源。

步骤三：利用克里金插值法对离散数据源进行插值，生成网格水深数据。

步骤二得到的水深数据是离散分布的，为了生成等值线必须通过插值生成网格水深数据。克里金法作为一种高精度的插值方法被广泛地应用在很多实际领域。该算法充分利用空间数据集自身的变异特性来确定待估计点的有效关联的区域范围，然后通过选取该区域范围内的点集计算待估计点的具体属性。因此，针对水深数据所表现的随机性和结构性的地理空间分布特征，使用克里金法可以满足要求。克里金插值原理为：假设某一样品分布变量为h(x_i)，空间上不同的位置点之间存在着相互的联系，这种相互的联系可以通过一系列的权值表示。克里金插值的关键计算规则是利用待估计点周围的n个已知点的h(x_i)数值来估算待估计点的h(x₀)数值，其表达式为：其中，λ_i是变量h(x_i)的权系数。根据克里金法的无偏估计规则，变量权系数λ_i(i＝1,2...n),满足公式：具体插值方法如下：

步骤3.1将航海作业区域Ω划分为M×N个网格，其中M为沿经度方向的网格数，N为沿纬度方向的网格数。G_ij(x_i,y_j)表示第i行，第j列处的网格，xi为该网格中心的经度，y_j为该网格中心的纬度，网格边长l的计算公式如下为 $l=\mathrm{Max}(\frac{\mathrm{fabs}(\mathrm{Blon}-\mathrm{Elon})}{N},\frac{\mathrm{fabs}(\mathrm{Blat}-\mathrm{Elat})}{N}).$ M、N的大小可以根据需要进行设置，以保证l≤0.05度为宜。

步骤3.2利用克里金插值算法逐一对网格G_ij(x_i,y_j)处的水深h_ij进行插值。具体实现方法如下：取以(x_i,y_j)为中心，边长为L的正方形区域内的水深数据为参考水深，其中L=2ωl；ω＞4，具体取值以保证参数水深的个数W大于3为宜；l为网格边长。h_ij的计算公式为其中，h^k为第k个参考水深的值。λ_k的计算公式为

步骤四：利用关联表法生成水深等值线。

步骤三获得的网格水深数据中矩形网格各节点间以边互相连接，每四条边又组成了一个网格面，除边界外，每条边均为两个网格面所共有。因此，可以建立节点链表、边链表和面链表来表示点、边、面之间的拓扑关系。然后基于关联表进行等值线追踪和搜索。具体步骤如下：

步骤4.1关联表的建立

（1）节点链表的设计

将步骤三得到的每个网格视为一个节点，用节点链表存储所有节点信息（经度、纬度和深度值），对每个节点标记唯一的索引号，节点链表数据结构设计如下：

（2）边链表的设计

边链表存储每条边所包含的两个节点信息和每条边所属的两个邻面信息，并对每条网格边标记唯一的索引号。边链表的数据结构设计如下：

在上述数据结构中，对于横向网格边，index_pt1，index_pt2分别表示该边所包含的左端节点和右端节点的索引号，index_regon1，index_regon2分别表示该边所属的上侧网格面和下侧网格面的索引号。对于纵向网格边，index_pt1，index_pt2分别表示该边所包含的上端节点和下端节点的索引号；index_regon1，index_regon2分别表示该边所属的左侧网格面和右侧网格面的索引号。同时，因为计算等值点时要区分横向网格边和纵向网格边，所以还需要在边链表中增加一个量sign用来标记网格边的种类，其中sign＝1表示该边为横向网格边，sign＝0表示该边为纵向网格边。

（3）面链表的设计

面链表存储它所包含的四条边在边链表中的索引号，并对每个网格面标记唯一的索引号，其数据结构设计如下：

（4）索引号的确定

为了便于计算，必须按照一定的次序来确定点、边、面的索引号。对于点、边、面的编号，均从网格左上角开始，按照从左至右，从上至下的顺序逐个进行标记。以保证在已知边链表中某条边的索引号时，可以计算出它所对应的点和面的索引号。已知某个面的索引号时，可以计算出其所对应的四条边的索引号。另外，对于处于网格区域边界处的网格边，由于它们只被一个面所有，记录它们所属的网格面信息时要作特殊处理。具体方法为：对于位于网格区域底边边界的网格边，把它们在边链表中的下侧网格面的索引号标记为-1。利用同样的方法处理位于网格区域最上边、最左边、最右边的网格边所属网格面的编号。

步骤4.2等值线追踪

等值线追踪的目的是按一定次序得到某属性值的所有等值点信息，这就要求等值线追踪时要扫描边链表中的每一条边，并且每条边只能被扫描一次，才能保证等值线搜索的正确性和完整性，本专利采用边搜索边利用标记数组进行标记已被扫描过的边的方法来达到上述效果。

在进行等值线追踪前首先将航海作业区域Ω内的水深范围按照一定的间隔进行划分，计算出每一条要追踪的等值线属性值。设航海作业区域Ω的最大水深为H_max，最小水深为H_min，等值线属性间隔为Δh，那么第i条等值线的属性值H_i可以通过下式计算出：

H_i＝H_min+Δh×(i-1)

利用双向追踪法对于任一属性值的等值线进行追踪，具体实现步骤如下：

步骤4.2.1建立搜索标记数组

针对某属性值等值线的追踪建立唯一的边搜索标记数组，该数组中的元素和边链表中的边一一对应，用于标记边链表中索引号与该数组下标相同的边是否被搜索过。数组中所有元素初始化为FALSE，每追踪到一个等值点，就把该搜索标记数组中的对应元素设为TRUE。当该属性值等值线所对应的边搜索标记数组中的所有元素值均为TRUE时，则表明边链表中的所有边都已经被搜索过，该属性值的等值线搜索完成。

步骤4.2.2追踪某一条分支等值线

查询边搜索标记数组，任何一条没有被搜索过的边都可以作为等值线起点，并从该点开始双向追踪该分支等值线，建立该分支等值线的等值点链表并把等值点信息存入文件。具体方法如图3所示，即：

a.初始化前向交点链表和后向交点链表

首先要初始化前向交点链表ForArray和后向交点链表BackArray，分别用来存储前向追踪和后向追踪过程中所得到的等值点信息。

b.寻找分支等值线线头

扫描边链表，找到一个等值点，作为分支等值线的线头，并把该等值点存入前向交点链表ForArray中。

c.前向追踪等值线

以步骤b中找到的等值点为起点，进行等值线前向追踪，当追踪到达区域边界时，即通过边链表检索到的面的索引值为-1时，或者在网格内找不到等值线出口时，该方向上等值线追踪结束。

d.后向追踪等值线

同样以等值线线头为起点，进行等值线后向追踪，当追踪到达区域边界时，即通过边链表检索到的面的索引值为-1时，或者在网格内找不到等值线出口时，该方向上等值线追踪结束。

e.建立分支等值线所有等值点链表

根据前向交点表ForArray和后向交点表BackArray，建立包含该条分支等值线所有等值点的等值点链表AllArray。具体分为如下两种情况：

当BackArray为空链表时，说明前向追踪已经完成了该条分支等值线的追踪，也就是说，等值线追踪线尾与线头在同一网格，该条分支等值为闭曲线。此时，要在ForArray链表末尾添加一个节点，该节点信息与链表的初始节点信息完全相同，即把闭曲线的起点重复记录一次当作闭曲线终点，这样得到的链表即为所求的等值点链表AllArray。

当BackArray不为空链表时，说明该分支等值线为开曲线，前向追踪和后向追踪都到达区域边界。此时，首先要把前向交点链表ForArray逆向，然后连接逆向后的ForArray和后向交点链表BackArray，得到的新链表就是所求的等值点链表AllArray。

步骤4.2.3判断是否完成该属性值的所有分支等值线的绘制

遍历边搜索标记数组，判断是否已经把边链表中的所有边搜索过。如果边链表中还有边没有被搜索过，则执行步骤4.2.2，继续搜索等值线。如果边链表中的所有边都已经被搜索过，则该属性值的所有等值线搜索完毕。

步骤五：等值线平滑和标注

直接把步骤四获得的等值点依次连接起来，生成的等值线不光滑，因此需要对等值线进行平滑处理，本专利采用三次B样条算法对等值线进行平滑处理。

另外，为了方便用户查看，每条等值线的属性值必须标注出来。通过读取记录等值点的数据文件，建立等值线标注数组，为标注等值线做好准备工作。具体方法如下：从等值点文件中以一定间隔读取等值点数据，并将等值点数据以等值线标注结构体的形式存入等值线标记数组。然后根据等值线标记数组中存储的等值点坐标确定要写标注的矩形区域。判断该矩形区域是否在航海作业区域Ω外，如果不在Ω外，则可以标注等值线。某条等值线的标注数目分闭曲线和开曲线区别对待，如果是闭曲线，只标注一次，如果是开曲线，则在不同的位置最多标注3次。等值线标注结构体设计如下：

步骤六：水深等值线与电子海图叠加显示。

为了方便用户的使用，需要将水深等值线直观形象地呈现给航海作业人员。电子海图能够为航海作业信息的显示提供基础背景，以电子海图为背景实现水深等值线的绘制，能够方便用户了解水深的分布和变化情况。具体方法为：首先，从等值线文件中读取某数值等值线上的数据点；然后将等值点的经纬度坐标转换为海图坐标系下的图上坐标，具体转化方法由电子海图显示系统提供；在电子海图上将各个等值点连接成线即可得到该数值下的等值线；采用步骤五生成的等值线标记数组对绘制好的等值线进行标注。重复以上操作即可完成所有等值线与海图的叠加显示。

Claims (1)

1.一种基于电子海图的海水深度等值线生成方法，其特征在于：

（1）根据航海作业区域选择用于水深数据提取的电子海图：

（1.1）确定航海作业区域：

航海作业区域Ω的经纬度范围为（BLat～ELat），（BLon～ELon），水深数据提取海区Ω₁的范围为（BLat-temp～ELat+temp），（BLon-temp～ELon+temp），其中temp＞0，单位为度：

$\mathrm{temp}=\frac{\mathrm{Max}(\mathrm{fabs}(\mathrm{Blon}-\mathrm{Elon}),\mathrm{fabs}(\mathrm{Blat}-\mathrm{Elat}))}{5};$

（1.2）选择电子海图：

逐一判断海图数据库中的电子海图，如果水深数据提取海区Ω₁与当前电子海图的覆盖区域有交集，将电子海图的名字存入待提取海图列表MapList，直到所有比例尺的海图均判断完毕为止；

（2）从电子海图中提取水深数据，生成离散数据源：

（2.1）从MapList中选择一幅尚未进行水深提取的电子海图；

（2.2）打开选择的电子海图文件，找到水深层所在的位置；

（2.3）读取水深层内的离散水深数据，判断当前水深是否在海区Ω₁内，若水深所在位置在Ω₁内，则将水深信息经度、纬度和水深值存储到Db.dat文件中，重复操作直到此海图中的所有水深数据均进行判断；

（2.4）判断MapList中是否还有尚未进行水深提取的电子海图，若有则从中选择一幅电子海图并执行步骤（2.2），否则说明水深数据提取完毕，Db.dat存储的即是进行插值的离散数据源；

（3）利用克里金插值法对离散数据源进行插值，生成网格水深数据：

（3.1）将航海作业区域Ω划分为M×N个网格，其中M为沿经度方向的网格数，N为沿纬度方向的网格数，G_ij(x_i,y_j)表示第i行，第j列处的网格，x_i为网格中心的经度，y_j为网格中心的纬度，网格边长l为

$l=\mathrm{Max}(\frac{\mathrm{fabs}(\mathrm{Blon}-\mathrm{Elon})}{N},\frac{\mathrm{fabs}(\mathrm{Blat}-\mathrm{Elat})}{N}),$

保持l≤0.05度；

（3.2）利用克里金插值算法逐一对网格G_ij(x_i,y_j)处的水深h_ij进行插值：取以(x_i,y_j)为中心，边长为L的正方形区域内的水深数据为参考水深，其中L=2ωl；ω＞4，保持参数水深的个数W大于3；l为网格边长，hk为第k个参考水深的值，

（4）利用关联表法生成水深等值线：

（4.1）建立关联表：

（4.1.1）节点链表：

用节点链表存储所有节点信息,经度、纬度和深度值：

（4.1.2）边链表：

边链表存储每条边所包含的两个节点信息和每条边所属的两个邻面信息：

（4.1.3）面链表：

面链表存储它所包含的四条边在边链表中的索引号，其数据结构设计如下：

（4.1.4）索引号：

从网格左上角开始，按照从左至右，从上至下的顺序逐个进行标记,对于处于网格区域边界处的网格边，对于位于网格区域边界的网格边，把边链表中的下侧网格面的索引号标记为-1；

（4.2）等值线追踪：

将航海作业区域Ω内的水深范围按照间隔进行划分，计算出每一条要追踪的等值线属性值，航海作业区域Ω的最大水深为H_max，最小水深为H_min，等值线属性间隔为Δh，第i条等值线的属性值H_i＝H_min+Δh×(i-1)；利用双向追踪法对于任一属性值的等值线进行追踪：

（4.2.1）建立搜索标记数组

边搜索标记数组中的元素与边链表中的边一一对应，数组中所有元素初始化为FALSE，每追踪到一个等值点，搜索标记数组中的对应元素设为TRUE；

（4.2.2）追踪分支等值线

查询边搜索标记数组，以没有被搜索过的边作为等值线起点，开始双向追踪分支等值线，建立分支等值线的等值点链表并把等值点信息存入文件：

（4.2.2.1）初始化前向交点链表ForArray和后向交点链表BackArray，分别用来存储前向追踪和后向追踪过程中所得到的等值点信息；

（4.2.2.2）扫描边链表，找到等值点，作为分支等值线的线头，把等值点存入前向交点链表ForArray中；

（4.2.2.3）以步骤（4.2.2.2）中找到的等值点为起点，进行等值线前向追踪，当追踪到达区域边界时，即通过边链表检索到的面的索引值为-1时，或者在网格内找不到等值线出口时，等值线追踪结束；

（4.2.2.4）以等值线线头为起点，进行等值线后向追踪，当追踪到达区域边界时，即通过边链表检索到的面的索引值为-1时，或者在网格内找不到等值线出口时，等值线追踪结束；

（4.2.2.5）建立分支等值线所有等值点链表，根据前向交点表ForArray和后向交点表BackArray，建立包含分支等值线所有等值点的等值点链表AllArray；

（4.2.3）遍历边搜索标记数组，判断是否已经把边链表中的所有边搜索过，如果边链表中还有边没有被搜索过，则执行步骤（4.2.2），继续搜索等值线，如果边链表中的所有边都被搜索过，则属性值的所有等值线搜索完毕；

（5）等值线平滑和标注：

采用三次B样条算法对等值线进行平滑处理；

从等值点文件中以一定间隔读取等值点数据，并将等值点数据以等值线标注结构体的形式存入等值线标记数组，根据等值线标记数组中存储的等值点坐标确定要写标注的矩形区域，判断矩形区域是否在航海作业区域Ω外，如果不在Ω外，则可以标注等值线，等值线标注结构体为：

（6）水深等值线与电子海图叠加显示：

从等值线文件中读取数值等值线上的数据点，将等值点的经纬度坐标转换为海图坐标系下的图上坐标，在电子海图上将各等值点连接成线即可得到该数值下的等值线，采用步骤（5）生成的等值线标记数组对绘制好的等值线进行标注。