A method for constructing a high-precision digital water depth model for nautical charts with isobath constraints

Abstract

The present invention discloses a method for constructing a high-precision digital water depth model of a nautical chart with isobath constraints, and belongs to the field of marine surveying and mapping technology. The present invention mainly targets nautical chart data with relatively sparse water depth data, and combines elements such as isobaths, coastlines, and obstructions on the nautical chart to construct a high-precision digital water depth model that conforms to the actual expression intent of the nautical chart. Compared with the existing method, the present invention effectively overcomes the problem of obvious deviation between the depth value of the local area and the expression intent of the nautical chart due to the dense and sparse water depth, and significantly improves the quality of the digital water depth model.

Description

Technical Field

[0001] The invention belongs to the technical field of ocean surveying and mapping, and relates to a method for constructing a high-precision digital water depth model of a nautical chart with isobath constraints. Background Art

[0002] Digital Depth Model (DDM) is a digital expression of the seabed surface morphology through limited water depth data. An important application of DDM is to ensure the safe navigation of ships at sea. Especially with the promotion and application of digital nautical charts, DDM, as an information model, can serve ship navigation more accurately and reliably compared with the traditional manual operation mode.

[0003] In recent years, many scholars have explored interpolation methods for various sources of depth data (including multi-beam depth data, single-beam depth data, chart depth data, etc.) under different circumstances. Among them, some methods are directly borrowed from DEM, and the interpolation of elevation values ​​in DEM is introduced into the depth interpolation of DDM; some methods are based on the DEM elevation interpolation method, and some improvements are made to the interpolation method based on the characteristics of the depth data itself, making it more suitable for DDM interpolation under specific conditions. All of the above methods have achieved good construction results under their specific conditions. However, it should be noted that most of the above methods are aimed at high-density data with dense depths, while nautical chart data with sparse depth density are less considered.

[0004] Nautical charts often use the coordination of water depth and isobaths to express the seabed topography. In the actual process of nautical chart making, in order to make the map expression beautiful and clear and easy to read, the water depth of some element-dense areas is deleted according to the specification requirements, and the isobaths that are close to each other are selected, interrupted, merged, etc. Therefore, the isobaths and water depth data on the nautical chart are actually "simple" data after cartographic synthesis. When this type of nautical chart is used as the data source and the DDM is directly constructed using the above method, the construction result will deviate from the actual nautical chart expression intention, and even the error of "expanding depth" will occur. For this reason, there is an urgent need for a model that is specifically used for constructing DDM for nautical charts based on nautical chart data to improve the accuracy of DDM and expand the application scope of nautical charts. Summary of the invention

[0005] In order to solve the deviation problem that is easy to occur in the current DDM construction method when constructing nautical charts, the present invention provides a high-precision DDM construction method for nautical charts based on isobath constraints.

[0006] The technical solution adopted by the present invention to achieve the above-mentioned purpose is:

[0007] A method for constructing a high-precision digital water depth model of a nautical chart with isobath constraints, comprising the following steps:

[0008] S1: Input nautical chart data, and extract the coastline, isobath (including beach line), special line (chart outline, mandatory connection line, additional buffer zone, etc.) and water depth point and other elements according to the nautical chart data;

[0009] S2: According to the relationship between the isobath area (the isobath area in the present invention refers to the polygonal area surrounded by isobaths, shorelines, beach lines, etc.), isobaths and water depth points, the depth range value of the isobath area is automatically analyzed and extracted (the depth range of the isobath area in the present invention refers to the depth value range of the polygonal area surrounded by isobaths, shorelines, beach lines, etc.);

[0010] S3: Make necessary preliminary preparations for the data based on the input coastline, isobaths (including beach lines), water depth points and other elements;

[0011] S4: Construct a constrained Delaunay triangulation network and classify triangles according to the water depth of the triangle vertices and the properties of the three sides, and automatically identify triangles that cannot be directly interpolated;

[0012] S5: Construct a ternary tree model. According to the triangle type obtained in step S4, the triangles that cannot be directly interpolated are systematically organized using a ternary tree model to identify the area to be interpolated, laying the foundation for subsequent interpolation transfer;

[0013] S6: Use the "source tracing" algorithm to transfer the interpolation value to the interpolation area. Based on the above ternary tree model, use the "source tracing" algorithm to transfer the interpolation value to the interpolation area, and record the transferred interpolation value on the triangle edge;

[0014] S7: Convert and output the constructed high-precision digital water depth model. Based on the above interpolated triangles, the triangles are further split, and each triangle area is interpolated according to the current mainstream interpolation method, and the interpolated high-precision grid digital water depth model is output.

[0015] Furthermore, in step S2, the method for automatically analyzing and extracting the depth range value of the equal depth zone is:

[0016] S2.1: According to the attributes of the isobath polygon in the nautical chart, if the attribute is land or dry beach area, the upper and lower limits of its depth range are directly set to 0m; if the attribute is shallow area, sub-shallow area or deep area, the depth range value is further determined based on the isobath and water depth information, and then go to step S2.2;

[0017] S2.2: Based on the relationship between the isobath polygon and the isobaths, all the isobaths (including shorelines and beach lines) that constitute the isobath are extracted. To increase the reliability of the algorithm, the isobaths with a length less than 0.25 cm (this value is the empirical value obtained after multiple tests) are directly discarded, and then the maximum depth value Depth of the above isobaths is extracted. max ;

[0018] S2.3: According to Table 1, obtain the above maximum depth value Depth max The depth values ​​of the isobaths of the upper and lower layers (respectively denoted as D + , D ++ , D - , D -- For example, if the current scale is 1:250,000 and the current maximum water depth is 30m, then the depth values ​​D of the two upper and lower levels mentioned above are required. + , D ++ , D - , D -- Then, all normal water depth points (excluding dangerous water depth) in the current depth zone are extracted and compared with the depth values ​​of each level mentioned above one by one. - , D - ~Depth max 、Depth max ~D + , D + ~∞), and calculate their percentages, which are denoted as ρ D- , D-Depth , DepthD+ , D+ .

[0019] Table 1 Basic depth contours on maps of different scales

[0020]

[0021] S2.4: Based on the percentage calculated in step S2.3, make the following judgment:

[0022] like The depth range of the equal depth zone is set to D -- ~Depth max ;

[0023] On the contrary The depth range of the equal depth zone is set to D - ~Depth max ;

[0024] On the contrary Then set the depth range of the equal depth area to Depthmax ~D ++ ;

[0025] On the contrary Then set the depth range of the equal depth area to Depth max ~D + .

[0026] Furthermore, in step S3, the main steps of data preparation are as follows:

[0027] S3.1: According to the input nautical chart data, traverse the equal depth polygons in the data to obtain the depth range value;

[0028] S3.2: Extract the islands, reefs, water depth and other elements in the current bathymetric zone, and compare these elements with the depth range of the bathymetric zone one by one. If the depth is less than the lower limit of the bathymetric zone, mark the element as "to be processed", otherwise no processing is performed;

[0029] S3.3: Traverse all "elements to be processed", take the element point as the center, and construct a buffer zone (the present invention refers to these buffer zones as additional buffer zones) with a smaller distance (distance value is 0.01 cm) as the radius, and set the depth of the buffer zone boundary line to the lower limit of the isobath zone to prevent the interpolation value from deviating from the actual expression intention of the nautical chart due to the element depth being shallower than the lower limit of the isobath zone.

[0030] Furthermore, in step S4, the triangle classification and identification method of triangles that cannot be directly interpolated are constructed as follows:

[0031] S4.1: Construct a constrained Delaunay triangulation based on the coastline, isobaths (including beach lines), special lines (outline lines, mandatory connection lines, additional buffer zones, etc.), and water depth points, and record the triangle vertex numbers;

[0032] S4.2: Based on the coordinates of the centroid of each triangle, obtain the isobath polygon in the surface layer where the triangle is located, and based on the upper and lower limit depth values ​​of the isobath polygon, determine whether the depth values ​​of the three vertices of the triangle are less than the lower limit of the depth range of the isobath; if so, record the depth value at the vertex as the lower limit of the depth range of the isobath;

[0033] S4.3: After the constrained Delaunay triangulation is completed, in order to effectively identify the areas that cannot be directly interpolated and facilitate the subsequent construction of the ternary tree model, the triangles are divided into 11 types of triangles (such as the attached triangles) based on the triangle vertex attributes (point islands are also treated as curve points), the corrected depth value (whether the depth is equal), and the three-side attributes. Figure 4 shown);

[0034] S4.4: According to the triangle type classification method of the present invention (as shown in the attached Figure 4 After the constrained Delaunay triangulation is constructed, the area on the map that cannot be directly interpolated should be the attached Figure 4 The areas where the V-type triangle, VI-type triangle, VIII-type triangle, X-type triangle, and XI-type triangle are located are shown. For ease of description, the non-constrained edges of the above triangles with the same vertex depth value are called "expandable edges".

[0035] Furthermore, in step S5, in order to systematically organize various types of triangles, the present invention adopts a ternary tree model to organize the triangle data that cannot be directly interpolated in step S4, and stores the parent node Father, three child nodes Son, all triangles allTri in the node area, and whether the node can be directly processed IsPossible in the ternary tree node.

[0036] Since the XI-type triangles among the expandable triangles obtained in step S4 have the most expandable situations and are the most complex, this type of triangle is processed first, and the idea of ​​permutation and combination is adopted to process them one by one in the order of XI-type, X-type, V-type, and VI-type triangles to ensure that no conflict occurs during the triangle expansion process. For the convenience of subsequent processing, the present invention marks all triangles as "false", and if they are processed, they are set to "true". At the same time, the corresponding interpolation value subDepths is set at the midpoint of the three sides of the triangle, and the initial value is -1.

[0037] The basic steps of constructing the ternary tree are as follows:

[0038] a. First, traverse the incoming triangles and obtain the unprocessed type XI triangles. If such type XI triangles exist, store them and go to step c. If such type of triangles do not exist, go to step b.

[0039] b. Traverse all triangles and obtain unprocessed X-type, V-type, and VI-type triangles. If such triangles exist, store them and go to step e. If such triangles do not exist, it indicates that all triangles to be processed have been organized and the ternary tree model has been constructed. Output the root nodes of all tree models.

[0040] c. Create the root node T of the ternary tree root , store the XI-type triangle passed in step a in the root node, set its processing flag to "true", and start expanding to the expandable edges on its three sides. For example, if the triangle passed in is a XI-type triangle m0, store it directly in the created root node T root , and create three child nodes for the root node. rootThe last triangle m0 is stored in allTri of the three child nodes, and the expansion direction is specified for it. The root node is set as the father node of the above three child nodes, and the three child nodes are added one by one. Pass in step d;

[0041] d. Based on the incoming child nodes and its specified expansion direction, and stores the adjacent triangles in the specified expansion direction in the child node In the process, set its processing flag to "true", and then judge the child node The type of the last triangle. If it is a type VIII triangle, the adjacent triangles in the extension direction are stored in the child node. Continue to judge The type of the last triangle, until the last triangle is not a type VIII triangle; if the last triangle is a type XI triangle, it is a child node Create two child nodes and store the last triangle in the child nodes In the two child nodes of , specify the expansion direction, set the child node For the parent node (Father) of the above two child nodes, repeat step d; if the last triangle is a V, VI, or X type triangle, stop the expansion operation, and go to step a after all child nodes have stopped expanding;

[0042] e. Create the root node T of the ternary tree root , store the triangle passed in step b into the root node T root In the example above, set its processing flag to "true" and start expanding according to the triangle type. The specific method is to get the root node T root The adjacent triangles in the expandable direction of the last triangle are stored in node T root In the process, set its processing flag to "true" and judge the root node T root The type of the last triangle, if it is a type VIII triangle, then according to the expansion direction, its adjacent triangles are stored in the root node T root In step 2, continue to determine the type of the last triangle until the last triangle is not a type VIII triangle; if the last triangle is a type V, VI, or X triangle, stop expanding and go to step b.

[0043] Furthermore, the basic method of the "tracing source" algorithm in step S6 is as follows:

[0044] S6.1: Input the root nodes in turn and determine whether the root node has child nodes. If not, it indicates that the node is a single node and can be directly interpolated and transferred. Interpolation and transfer processing is directly performed on the node; if child nodes exist, these child nodes are passed to step S6.2;

[0045] S6.2: Input the child nodes in sequence, and judge whether the node can be interpolated (whether there is a difference between the depth values ​​of the first and last triangles in the node and the expandable edge). If interpolation can be performed directly, perform the interpolation operation and return the value "true"; if interpolation cannot be performed, determine whether the child node has child nodes. If so, pass these child nodes to step S6.2: recursion; if not, return the value "false"; after the child node judgment is completed, go to step S6.3;

[0046] S6.3: According to the return value of the child node in step S6.2, if the return value is "false", then directly return "false", otherwise, it indicates that there is a node for interpolation transfer in the child nodes of the node, so that one side of the end XI type triangle of this node is interpolated and transferred, and there is an additional water depth point (the water depth point is set at the midpoint of the expandable edge of the triangle during the interpolation transfer process, and its water depth value is calculated by interpolation); according to the additional water depth point, the node is interpolated and transferred, and all child nodes under this node (except the child nodes for interpolation transfer) are interpolated and transferred, and then the value "true" is returned. Until the node does not have a parent node, it indicates that all nodes under the root node have been interpolated and transferred, and then go to step S6.1.

[0047] The basic process of interpolation transfer is as follows: Taking the X-type and V-type triangles as the head and tail, the central axis is extracted using the constructed triangulated network as shown in the attached figure. Figure 6 middle As shown, due to the curve point P L7 and water depth point P D It is not equal to the value on the surrounding isobaths, so the depth value of each expandable edge midpoint can be calculated by linear interpolation. The calculation formula is as follows:

[0048]

[0049] Among them, Depth Psubi Represents the midpoint P of the expandable edge subi The depth value of the edge is referred to as the additional depth value of the edge in the present invention. L7 P subi | represents point P on the central axis L7 Solstice P subi The distance, |P L7 P D | represents point P on the central axis L7 Solstice P D distance.

[0050] Furthermore, the triangle splitting and interpolation method in step S7 is as follows:

[0051] First, split the triangle according to the number of additional points on the triangle:

[0052] If only one side has an additional point: For this type of triangle, directly connect the additional point and its corresponding vertex, directly divide it into two triangles, and use the current mainstream interpolation algorithm (there are many methods of this type and no requirements are made) to calculate the internal depth values ​​of the two triangles respectively.

[0053] If there are additional points on both sides: For triangles with two additional point depth values, first directly connect the two additional points, then calculate the distance between each additional point and its corresponding vertex, and connect the smaller distance, thereby dividing the triangle into three small triangles; finally, use the current mainstream interpolation algorithm to calculate the internal depth value of each triangle.

[0054] If there are additional points on all three sides: For a triangle with 3 additional points, connect the additional points in sequence to divide the triangle into 4 small triangles, and calculate the internal depth value of each triangle.

[0055] Beneficial effects of the present invention: The present invention aims at the current digital water depth model construction method, which is mainly aimed at water depth data with high density (such as multi-beam data). When the data is changed to a nautical chart, the existing method will appear to be contrary to the actual intention of the nautical chart. The present invention mainly aims at nautical chart data with relatively sparse water depth data, and combines the elements such as isobaths, coastlines and obstacles on the nautical chart to construct a high-precision digital water depth model that conforms to the actual intention of the nautical chart. Compared with the existing method, the present invention effectively overcomes the problem of obvious deviation between the depth value of the local area and the intention of the nautical chart due to the dense and sparse water depth, and significantly improves the quality of the digital water depth model. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 It is a general flow chart of the method for constructing a high-precision digital water depth model of a nautical chart with equal depth zone constraints of the present invention;

[0057] Figure 2 It is a flow chart of the automatic analysis and extraction method of depth range of isobath area of ​​the present invention;

[0058] Figure 3 It is the main flow chart of data preprocessing of the present invention;

[0059] Figure 4 It is a schematic diagram of constructing triangle classification in the present invention;

[0060] Figure 5 It is the main flow chart of the ternary tree construction of the present invention;

[0061] Figure 6 It is a schematic diagram of the basic process of interpolation transfer of the present invention;

[0062] Figure 7 It is the main flow chart of the interpolation transfer of the tracing algorithm of the present invention. DETAILED DESCRIPTION

[0063] The specific implementation modes of the present invention will be further described below in conjunction with the accompanying drawings.

[0064] The method for constructing a high-precision digital water depth model of a nautical chart with equal depth zone constraints proposed in the present invention is implemented by using a computer to construct a high-precision digital water depth model of a nautical chart. Assuming that there is an existing nautical chart A, the method of the present invention is used to construct a high-precision digital water depth model, including the following steps: Figure 1 As shown:

[0065] S1: Input nautical chart data, extract coastlines, isobaths (including beach lines), special lines (chart outlines, mandatory connecting lines, additional buffer zones, etc.), and water depth points and store them in memory. For the hypothetical nautical chart A, the set of extracted elements is called S(L);

[0066] S2: Based on the relationship between the isobaths (the isobaths refer to the polygonal area surrounded by isobaths, shorelines, beachlines, etc.), isobaths and water depth points in S(L), the Figure 2 The method shown automatically analyzes and extracts the depth range values ​​of the equal depth zones;

[0067] S3: Basis Figure 3 The method shown performs necessary preliminary preparation on the input feature set S(L) to obtain the feature set S′(L);

[0068] S4: Construct a constrained Delaunay triangulation for the feature set S′(L) to obtain a triangle set S(T), and use the water depth of the triangle vertices and the properties of the three sides to Figure 4 The triangles are classified in the manner shown;

[0069] S5: Based on the triangle type obtained in step S4, use Figure 5 The ternary tree (the ternary tree set S(Tree)) organization method shown is systematically organized to identify the area to be interpolated, laying the foundation for subsequent interpolation transfer;

[0070] S6: Based on the above ternary tree model, adopt Figure 7 The "tracing back" algorithm transfers the interpolation value to the interpolation area, and records the transferred interpolation value on the triangle edge, and obtains the triangle set Sr(T) after the interpolation transfer;

[0071] S7: Based on the above interpolated triangle set Sr(T), interpolate each triangle area according to the mainstream interpolation method, and output the interpolated high-precision nautical chart DDM.

[0072] The present invention is described by way of embodiments, and those skilled in the art will appreciate that various changes or equivalent substitutions may be made to these features and embodiments without departing from the spirit and scope of the present invention. In addition, under the teachings of the present invention, these features and embodiments may be modified to adapt to specific circumstances and materials without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application belong to the protection scope of the present invention.

Claims

1. A method for constructing a high-precision digital water depth model for a nautical chart with isobath constraints, characterized in that: The method comprises the following steps: S1: Input nautical chart data, and extract elements including coastline, isobath, special line and water depth point according to the nautical chart data; the special line includes outline line, mandatory connection line and additional buffer zone; S2: Automatically analyze and extract the depth range of the isobath area based on the relationship between the isobath area, isobath line and water depth point; S3: Prepare the data based on the input elements including shoreline, depth contour and water depth points; S4: Construct a constrained Delaunay triangulation network and classify triangles according to the water depth of the triangle vertices and the properties of the three sides, and automatically identify triangles that cannot be directly interpolated; S5: According to the triangle types classified in step S4, the triangles that cannot be directly interpolated are systematically organized using a ternary tree model to identify the area to be interpolated, laying a foundation for subsequent interpolation transfer; S6: Based on the ternary tree model in step S5, the "tracing back" algorithm is used to transfer the interpolation value to the interpolation area, and the transferred interpolation value is recorded on the triangle edge; S7: Based on the above interpolated triangles, further split the triangles, and interpolate each triangle area according to the current mainstream interpolation method, and output the interpolated high-precision grid digital water depth model; In step S2, the method for automatically analyzing and extracting the depth range value of the equal depth zone is as follows: S2.1: According to the attributes of the isobath polygon in the nautical chart, if the attribute is land or dry beach area, the upper and lower limits of its depth range are directly set to 0m; if the attribute is shallow area, sub-shallow area or deep area, the depth range value is further determined according to the isobath and water depth information, and then go to step S2.2; S2.2: Based on the relationship between the depth contour polygon and the depth contour line, extract all the depth contour lines that constitute the depth contour area, directly discard the depth contour lines with a length less than 0.25 cm, and then extract the maximum depth value Depth of the above depth contour lines. max ; S2.3: Obtain the above maximum depth value Depth according to the basic isobaths shown on the maps of different scales. max The depth values ​​of the depth contours of the upper and lower layers are recorded as D + , D ++ , D - , D -- ; Then extract all the normal water depth points in the current depth equalization area, compare them with the depth values ​​of each layer one by one, and count the depth values ​​of each layer 0~D respectively. - , D - ~Depth max 、Depth max ~D + , D + The number of water depths in the range of ∞ and their percentages are calculated and recorded as The basic isobaths shown on maps of different scales are as follows: When the scale is greater than 1:200,000, the basic isobaths are 0, 2, 5, 10, 20, 30, 50, 100, 200, in meters; When the scale is 1:200,000 to 1:490,000, the basic isobaths are 0, 5, 10, 20, 30, 50, 100, 200, in meters; When the scale is 1:500,000 to 1:1,000,000, the basic isobaths are 0, 5, 10, 20, 30, 50, 100, 200, 500, in meters; When the scale is 1:1,000,000, the basic isobaths are 0, 10, 20, 30, 50, 100, 200, 500, 1000, 2000, in meters; When the scale is less than 1:1,000,000, the basic isobaths are 0, 10, 50, 200, 500, 1000, 2000, 4000, 6000, in meters; S2.4: Based on the percentage calculated in step S2.3, make the following judgment: like The depth range of the equal depth zone is set to D -- ~Depth max ; On the contrary The depth range of the equal depth zone is set to D - ~Depth max ; On the contrary Then set the depth range of the equal depth area to Depth max ~D ++ ; On the contrary Then set the depth range of the equal depth area to Depth max ~D + ; In step S3, the steps of data preparatory are as follows: S3.1: According to the input nautical chart data, traverse the equal depth polygons in the data to obtain the depth range value; S3.2: Extract the elements including point islands, reefs and water depth in the current bathymetric area, and compare these elements with the depth range of the bathymetric area one by one. If the depth is less than the lower limit of the bathymetric area, mark the element as "to be processed", otherwise no processing is performed; S3.3: Traverse all "unprocessed elements", build a buffer with the element point as the center and a distance value of 0.01cm as the radius, and set the depth of the buffer boundary line to the lower limit of the bathymetric zone to prevent the interpolation value from deviating from the actual expression intention of the nautical chart due to the element depth being shallower than the lower limit of the bathymetric zone; In step S4, the triangle classification and identification method of triangles that cannot be directly interpolated are constructed as follows: S4.1: Construct a constrained Delaunay triangulation based on the elements including shorelines, isobaths, special lines, and water depth points, and record the triangle vertex numbers; S4.2: Based on the coordinates of the centroid of each triangle, obtain the isobath polygon in the surface layer where the triangle is located, and based on the upper and lower limit depth values ​​of the isobath polygon, determine whether the depth values ​​of the three vertices of the triangle are less than the lower limit of the depth range of the isobath; if so, record the depth value at the vertex as the lower limit of the depth range of the isobath; S4.3: After the constrained Delaunay triangulation is constructed, in order to effectively identify the areas that cannot be directly interpolated and facilitate the subsequent construction of the ternary tree model, the triangles are divided into 11 categories according to the triangle vertex attributes, the corrected depth values ​​and the three-edge attributes; S4.4: The areas on the graph that cannot be directly interpolated should be the areas where type V triangles, type VI triangles, type VIII triangles, type X triangles, and type XI triangles are located. For ease of description, the non-constrained edges of the above triangles with the same vertex depth value are called "extendable edges"; In step S5, in order to systematically organize triangles of various types, a ternary tree model is used to organize the triangle data that cannot be directly interpolated in step S4, and the parent node Father, three child nodes Son, all triangles allTri in the node area, and whether the node can be directly processed IsPossible are stored in the ternary tree node; First, process the XI triangles, and use the idea of ​​permutation and combination to process them one by one in the order of XI, X, V, and VI triangles to ensure that there is no conflict in the triangle expansion process; for the convenience of subsequent processing, all triangles are marked as "false", and if they are processed, they are set to "true"; at the same time, the corresponding interpolation value subDepths is set at the midpoint of the three sides of the triangle, and the initial value is -1; The basic steps of constructing the ternary tree are as follows: a. First, traverse the incoming triangles and obtain the unprocessed type XI triangles. If such type XI triangles exist, store them and go to step c. If such type of triangles do not exist, go to step b. b. Traverse all triangles and obtain unprocessed X-type, V-type, and VI-type triangles. If such triangles exist, store them and go to step e. If such triangles do not exist, it indicates that all triangles to be processed have been organized and the ternary tree model has been constructed. Output the root nodes of all tree models. c. Create the root node T of the ternary tree root , store the XI-type triangle passed in step a in the root node, set its processing flag to "true", and start expanding to the expandable edges on its three sides; if the passed-in triangle is a XI-type triangle m0, store it directly in the created root node T root , and create three child nodes for the root node. root The last triangle m0 is stored in allTri of the three child nodes, and the expansion direction is specified for it. The root node is set as the parent node Father of the above three child nodes, and the three child nodes are added one by one. Pass in step d; d. Based on the incoming child nodes and its specified expansion direction, and stores the adjacent triangles in the specified expansion direction in the child node In the process, set its processing flag to "true", and then judge the child node The type of the last triangle. If it is a type VIII triangle, the adjacent triangles in the extension direction are stored in the child node. Continue to judge The type of the last triangle, until the last triangle is not a type VIII triangle; if the last triangle is a type XI triangle, it is a child node Create two child nodes and store the last triangle in the child nodes In the two child nodes of , specify the expansion direction, set the child node Repeat step d for the parent node Father of the above two child nodes; if the last triangle is a V, VI, or X type triangle, stop the expansion operation, and go to step a after all child nodes have stopped expanding; e. Create the root node T of the ternary tree root , store the triangle passed in step b into the root node T root In the example, set its processing flag to "true" and start expanding according to the triangle type; the specific method is to obtain the root node T root The adjacent triangles in the expandable direction of the last triangle are stored in node T root In the process, set its processing flag to "true" and judge the root node T root The type of the last triangle, if it is a type VIII triangle, then according to the expansion direction, its adjacent triangles are stored in the root node T root In step 2, continue to determine the type of the last triangle until the last triangle is not a type VIII triangle; if the last triangle is a type V, VI, or X triangle, stop expanding and go to step b.

2. The method for constructing a high-precision digital water depth model of a nautical chart with isobath constraints according to claim 1, characterized in that: In step S6, the "tracing back" algorithm is as follows: S6.1: Input the root nodes in turn and determine whether the root node has child nodes. If not, it indicates that the node is a single node and the node is directly interpolated and transferred. If child nodes exist, these child nodes are passed to step S6.2; S6.2: Input the child nodes in sequence, and determine whether interpolation transfer is possible based on whether there is a difference between the depth values ​​of the first and last triangles in the node and the expandable edge. If interpolation transfer is possible directly, perform the interpolation transfer operation and return the value "true"; If interpolation transfer is not possible, determine whether the child node has child nodes. If so, pass these child nodes to step S6.2: recursion; if not, return the value "false"; after the child node determination is completed, go to step S6.3; S6.3: Based on the return value of the child node in step S6.2, if the return value is "false", then directly return "false", otherwise, it indicates that there is a node for interpolation transfer among the child nodes of the node, so that one side of the end XI type triangle of this node is interpolated and transferred, and there is an additional water depth point; based on the additional water depth point, this node is interpolated and transferred, and all child nodes under this node except the child nodes that have been interpolated and transferred are interpolated and transferred, and then the value "true" is returned; until the node has no parent node, it indicates that all nodes under the root node have been interpolated and transferred, and then go to step S6.

1.

3. The method for constructing a high-precision digital water depth model of a nautical chart with isobath constraints according to claim 1, characterized in that: In step S7, the triangle splitting and interpolation method is as follows: Split triangles based on the number of additional points on them: If there is an additional point on only one side: directly connect the additional point and its corresponding vertex, directly divide it into two triangles, and use the current mainstream interpolation algorithm to calculate the internal depth values ​​of the two triangles respectively; If there are additional points on both sides: first directly connect the two additional points, then calculate the distance between each additional point and its corresponding vertex, and connect the smaller distance, so as to divide the triangle into three small triangles; finally, the current mainstream interpolation algorithm is also used to calculate the internal depth value of each triangle; If there are additional points on all three sides: connect the additional points in sequence to divide the triangle into 4 small triangles, and calculate the internal depth value of each triangle.

4. The method for constructing a high-precision digital water depth model of a nautical chart with isobath constraints according to claim 2, characterized in that: The interpolation transfer process is as follows: Taking the X-type and V-type triangles as the head and tail as an example, the central axis is extracted using the constructed triangulated network Since the curve point P L7 and water depth point P D The depth values ​​of the midpoints of the expandable edges are not equal to the values ​​on the surrounding isobaths, so linear interpolation is used to calculate the depth values ​​of the midpoints of the expandable edges. The calculation formula is as follows: in, Represents the midpoint P of the expandable edge subi The depth value of |P L7 P subi | represents point P on the central axis L7 Solstice P subi The distance, |P L7 P D | represents point P on the central axis L7 Solstice P D distance.

Images