Traffic flow calculation method using spherical hexagonal grid vector lines

Abstract

The invention belongs to the technical field of ship traffic, and particularly relates to a traffic flow calculation method using a spherical hexagonal grid vector line. According to the method, a spherical hexagonal grid coding scheme is designed, a spherical hexagonal grid vector line generation algorithm is realized by means of coding addition operation, and a spherical hexagonal grid result corresponding to ship route data is generated on the basis and stored in a database; and the ship traffic flow calculation in any buffer area is realized through code matching and query operation in the database. Experiments prove that the ship traffic volume calculation method provided by the invention has good effects in the aspects of efficiency and accuracy, and has practical application potential.

Description

Technical Field

[0001] The invention belongs to the technical field of ship traffic, and in particular relates to a traffic flow calculation method using spherical hexagonal grid vector lines. Background Art

[0002] In order to improve the capacity and efficiency of maritime transportation, the number of ships has increased significantly, and the types of ships have shown a trend of specialization, large-scale and high-speed development. As a result, the traffic density of important channel waters, ports and waterways has continued to increase, and the navigation environment of waterways has become increasingly complex. Under this situation, research related to maritime transportation, such as ship abnormal behavior detection, ship route prediction, and ship route clustering, has gradually become a hot topic. For these studies, research related to ship routes and ship traffic volume is an important foundation.

[0003] Ship routes refer to the standardized navigation tracks that ships must follow under the influence of factors such as fuel efficiency, marine protected area policies, and security threats. They reflect the ship traffic rules under the combined effects of marine traffic planning policies, ship maneuvering behaviors, and hydrological characteristics, and can be used to characterize real ship trajectories. Ship traffic volume is a basic quantity that characterizes the actual situation of ship traffic in a certain water area. Calculating and analyzing ship traffic volume based on existing ship route data is the basis for research on route design, maritime management, and detection of marine traffic hotspots. At present, the ship traffic calculation methods can be roughly divided into three categories. The first category is to calculate the ship traffic by judging whether each route track point is in a specific water area, but this method may involve a large number of geometric operations in the process of judging the point-surface relationship; the second category is to generate a route grid map based on the square grid where the route track points are located, and obtain the corresponding ship traffic volume of each grid by superimposing the route grid map, but there is an error between the route grid map generated based on the track points and the route grid result, which affects the calculation accuracy of the route traffic volume; the third category is to use the geographic grid model to divide the research sea area into square grids, and determine the route traffic volume by judging whether the ship route intersects with the grid. Although this method introduces the grid idea into the ship traffic calculation, it only divides the target water area into grids without gridding the ship route. It may still involve a lot of geometric calculations when judging the relationship between the route and the grid.

[0004] Compared with square grids and other types of grids, hexagonal grids have more advantages in plane coverage, unit adjacency, etc. After decades of development, the technical theories related to the hexagonal grid system have become mature and gradually put into use. Sahr designed a "Central Place Indexing" (CPI) scheme for the mixed aperture grid system and obtained a US patent authorization (Sahr K. 2012. Central Place Indexing Systems. US Patented. United States Patent Awarded: US2010054550); Robertson developed a data integration and analysis system IDEAS based on the hexagonal grid, which can be used for large-scale data modeling and analysis (ROBERTSON C, CHAUDHURIC, HOJATI M, et al. An integrated environmental analytics system (IDEAS) based on a DGGS [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2020, 162: 214-228); Canada's PYXIS Innovation company developed a spatial data fusion service system PYXIS WorldView Studio based on the ISEA3H grid (Peterson P. 2006. Close-Packed Uniformly Adjacent, Multi-resolution Overlapping Spatial Data Ordering. United States Patent Application: 2006 / 0265197A1); The European Space Agency uses hexagonal grids to store and process remote sensing image data in the "Soil Moisture and Ocean Salinity Mission" of the "Living Planet Project" (Dumedah G, Walker JP, Rudiger C. Can SMOS Data be Used Directly on the 15-km Discrete Global Grid? [J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(5): 2538-2544.).The expression of vector lines in hexagonal grids is the key to whether vector data can be well integrated with the grid system. Many scholars have conducted in-depth research on this. Amiri et al. realized the gridding of vector data in hexagonal images and proposed an effective management scheme for gridded data (Mahdavi-Amiri A, Alderson T, Samavati F. Data Management Possibilities for Aperture 3 Hexagonal Discrete Global Grid Systems[R]. Science, 2016.); Tsatcha et al. proposed a bidirectional hexagonal grid path-finding algorithm and applied it to the automatic planning of the optimal navigation route (Tsatcha, Dieudonné, Saux, éric, Claramunt CA bidirectional path-finding algorithm and data structure for maritime routing[J]. International Journal of Geographical Information Science, 2014, 28(7): 1355-1377.); Liu Yongkui et al. realized the planar hexagonal grid vector line expression algorithm by selecting the grid cells closest to the vector line one by one (L. Yong-Kui. The Generation of. Straight Lines on Hexagonal Grids.[J].Comput.Graph.Forum,1993,12(1):27-32.); Yu Wenshuai et al. studied the accuracy control of vector line expression in hexagonal grid system, which promoted the application of vector line grid (YU Wenshuai, TONG Xiaochong, BEN Jin, XIE Jinhua. The Accuracy Control in the Process of Vector Line Data Drawing in the Hexagon Discrete Global Grid System[J].Journal of Geo-information Science,2015,17(7):804-809https: / / doi.org / 10.3724 / SP.J.1047.2015.00804); Vince established a unified mathematical model for vector line grid expression for multidimensional grids that can be obtained by translating a single Voronoi unit, which can be applied to planar hexagonal grids (Vince A.Discrete Lines and Wandering Paths[J].Siam Journal on Discrete Mathematics,2007,21(3):647-661.); Based on Vince's theory, Du Lingyu et al. proposed a planar hexagonal grid vector line generation algorithm based on dimensionality reduction (Du L,Ma Q,Ben J,et al.Duality and DimensionalityReduction Discrete Line Generation Algorithm for a Triangular Grid[J].International Journal of Geo-Information,2018,7(10):391.DOI:10.3390 / ijgi7100391.).

[0005] In general, there are abundant studies on ship traffic calculation and hexagonal grid vector line generation algorithms, but few studies combine the two to realize ship traffic calculation in arbitrary waters with the help of hexagonal gridding results of ship routes. Summary of the invention

[0006] The present invention aims to solve the problem that the existing ship traffic volume calculation method has a complex calculation process, a large amount of calculation and thus a relatively low efficiency, and provides a traffic flow calculation method using spherical hexagonal grid vector lines. The method of the present invention shows good results in both calculation efficiency and accuracy.

[0007] The present invention specifically adopts the following technical solutions:

[0008] Traffic flow calculation method using spherical hexagonal grid vector lines, including:

[0009] S1: Generate grid data corresponding to the ship route from the compressed ship route data using a spherical hexagonal grid vector line generation algorithm, and store it in a database; the compression is: adopt a route data compression method based on square tiles, divide the water area where the route data is located into square tiles at the same latitude and longitude intervals, and when continuous track points fall in the same tile, only retain one representative track point.

[0010] S1 specifically includes:

[0011] S1.1: compress the ship route data;

[0012] S1.2: storing the compressed ship route data in the form of a table in the database; extracting key fields from the processed ship route data to form entries in the TraceItem table in the database, wherein the key fields include mmsi, lon, lat, and dt_pos_utc fields; using a spherical hexagonal grid vector line generation algorithm to generate grid data corresponding to the ship route from the processed ship route data, and storing the grid data in the form of a table in the database;

[0013] The spherical hexagonal grid vector line generation algorithm includes:

[0014] A: Plane hexagonal grid vector line generation based on vector tracking:

[0015] Based on the vertical distance from the center of the plane hexagonal grid unit to the vector line, the vector line is tracked using the adjacent vector to generate the vector line grid result; the hexagonal grid adopts a hierarchical recursive method to divide the grid unit, and encodes it with {00, 01, 02, 03, 04, 05, 06, 10, 20, 30, 40, 50, 60} as the basic code element. For any vector line with the center of the grid unit as the starting point in the plane hexagonal grid, there are two direction vectors closest to its extension direction, which are the adjacent vectors.

[0016] B: The planar algorithm is extended to the spherical hexagonal grid system using Snyder equal-area projection. For spherical vector lines, three "nearby vectors" whose angles with the sphere are less than 180 degrees are selected for tracking; the three spherical nearby vectors are replaced by the corresponding basic code elements {basecode1, basecode2, basecode3} in {01, 02, 03, 04, 05, 06}, and coded addition operations are used to realize vector tracking in the spherical hexagonal grid.

[0017] Step B specifically includes:

[0018] B-1: Convert the spherical vector line endpoint coordinates (lon1, lat1) (lon2, lat2) to grid codes (Code1, Code2); lon represents longitude, lat represents latitude, and Code represents grid code;

[0019] B-2: Determine the basic code elements {basecode1, basecode2, basecode3} of the corresponding direction, where the basic code elements of the corresponding direction refer to: taking (lon1, lat1) and (lon2, lat2) as the starting points, determining a vector line on the sphere, and then determining three "adjacent vectors" whose spherical angles with the vector line are less than 180 degrees, and the three basic code elements corresponding to the directions pointed by these three "adjacent vectors"; Basecode represents the basic code element;

[0020] B-3: Determine the spherical grid cells one by one using coded addition operations and point-to-surface vertical distances i , output spherical gridding results {Grid1, Grid2, Grid3..., Grid n}; Grid represents a grid unit. The point-to-plane vertical distance refers to the vertical distance d from the center V of the grid unit to the plane determined by the endpoint AB of the vector line and the center O of the sphere. VQ Assume that the radius of the sphere is R, and the vertical distance d from the center of the grid unit V to the plane determined by the endpoint AB of the vector line and the center of the sphere O is VQ The spherical distance d from the cell center to the vector line VN The following conversion relationship exists:

[0021]

[0022] S2: Determine the grid cells in the buffer zone, and obtain the ship traffic volume by matching the grid cell codes in the buffer zone with the grid cell codes in the grid data corresponding to the ship routes in the database. The grid cells in the buffer zone are determined as follows: the buffer zone is determined by the geographical coordinates of the center point and the radius, and the grid cells whose centers fall within the buffer zone are regarded as the grid cells in the buffer zone.

[0023] S2 specifically includes:

[0024] S2.1: In the TraceItem table, according to the buffer center coordinates, radius, and query time interval, filter the ship routes that meet the conditions to form a set T: {trace1, trace2, trace3, ..., trace n}, trace represents the selected ship route;

[0025] S2.2: Determine the grid cell code set C in the buffer zone according to the buffer zone radius: {code1, code2, code3, ..., code m}(m≤19).

[0026] S2.3: trace each route in set T i (i≤n) corresponds to the 18th layer of grid data and each grid code in the set C m (m≤19) to query and match, if trace i (i≤n) The grid data corresponding to the grid contains code m (m≤1) and the corresponding timestamp is in the query time interval, the traffic volume L is accumulated. After all queries are completed, L is the traffic volume calculation result.

[0027] The method of the invention can be applied to the calculation of ship traffic volume in international energy transportation.

[0028] The beneficial effects of the present invention are:

[0029] 1. The present invention combines the excellent characteristics of hexagonal grids, such as hierarchical decomposition, unique grid unit coding, and calculable coding, with the calculation of ship traffic. On the basis of designing a spherical hexagonal grid coding scheme, the present invention proposes a spherical hexagonal grid vector line generation algorithm based on vector tracking, and fully utilizes the advantage of the discreteness of the grid being suitable for computer storage, pre-stores the gridding results of ship routes in the database, and uses the grid coding hierarchical relationship to perform data query and matching, thereby realizing a traffic flow calculation method using spherical hexagonal grid vector lines. The method of the present invention converts the spherical geometric operations originally required for calculating ship traffic into simple database query operations by exchanging storage space for query time, thereby improving the computational efficiency of ship traffic calculation. After experimental verification, the traffic flow calculation method using spherical hexagonal grid vector lines proposed by the present invention has shown good results in both efficiency and accuracy, and has practical application potential.

[0030] 2. The vector tracking-based spherical hexagonal grid vector line generation algorithm of the present invention proposes the concepts of "direction vector" and "adjacent vector", takes the vertical distance from the grid unit center to the vector line as the basis, and uses the "adjacent vector" to track the vector line, thereby generating a vector line grid result. And the algorithm is extended to the sphere with the help of grid system code addition operation.

[0031] 3. The strategy for determining hexagonal grid units in the buffer zone of the present invention: first determine the nth layer of grids whose grid unit side length is just less than the regional radius, and regard the n+1th layer of grids as the applicable grids based on this. Determine the 19 grid units in the unit where the regional center is located and its surrounding adjacent units in the n+1th layer of grids through proximity query and coded addition operations, and determine the grid units in the region based on the spherical distance from the center of each grid unit to the regional center. This method adaptively determines the grid level and the corresponding level grid units in the region according to the regional radius while taking into account both accuracy and computational efficiency.

[0032] 4. Traffic flow calculation using spherical hexagonal grid vector lines of the present invention: a storage structure of ship route data in a database is designed, and the characteristics of the spherical hexagonal grid system that can be hierarchically decomposed and the grid unit code is unique are utilized. Combined with the spherical hexagonal grid vector line generation algorithm, the ship routes are pre-gridded in the hexagonal grid system and stored in the database, and the ship traffic volume is further calculated through database query operations based on conditions such as location, time, and grid code. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 This is a schematic diagram of the grid division and coding of the HLST scheme.

[0034] Figure 2 Schematic diagram of encoding addition operation for the HLST scheme.

[0035] Figure 3 This is a schematic diagram of the conversion between geographic coordinates (latitude and longitude) and grid codes.

[0036] Figure 4 Schematic diagram of the algorithm for generating vector lines in a hexagonal grid based on vector tracking. (a) Schematic diagram of “nearby vectors”; (b) Schematic diagram of vector tracking.

[0037] Figure 5 Schematic diagram of Snyder projection deformation.

[0038] Figure 6 This is a schematic diagram of the "nearby vector" on a sphere.

[0039] Figure 7 This is a schematic diagram of the relationship between the spherical distance and the perpendicular distance between the point and the plane.

[0040] Figure 8 Schematic diagram of the algorithm for generating vector lines in a planar hexagonal grid based on vector tracing.

[0041] Fig. 9 Schematic diagram of the spherical hexagonal grid vector line generation algorithm based on vector tracking. (a) The fifth-layer grid vector line generation effect; (b) The sixth-layer grid vector line generation effect.

[0042] Fig.10 This is the flow chart for entering ship route data into the database.

[0043] Fig.11 Diagram of the strategy for determining grid cells within the buffer zone.

[0044] Fig.12 Schematic diagram of extreme radius.

[0045] Fig.13 Schematic diagram of grid cells in the buffer zone under extreme radius conditions.

[0046] Fig.14 Flowchart for vessel traffic calculation.

[0047] Fig.15 This is a schematic diagram of the experimental area and ship route distribution.

[0048] Fig.16 Calculate time statistics for different data volume traffic.

[0049] Fig.17 This is a schematic diagram of intersection analysis.

[0050] Fig.18 Calculate average deviation rate statistics for traffic volume.

[0051] Fig.19 Schematic diagram of the non-overlapping area between the circular buffer and the grid cells. DETAILED DESCRIPTION

[0052] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution in the embodiment of the present invention will be clearly described below in conjunction with the drawings in the embodiment of the present invention. Obviously, the described embodiment is a part of the embodiment of the present invention, not all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0053] Terminology explanation:

[0054] The hexagonal grid is a divisible multi-level grid data model with hexagons as basic units.

[0055] Ship routes refer to the standardized navigation paths that ships must follow, influenced by factors such as fuel efficiency, marine protected area policies, and security threats.

[0056] Ship traffic volume refers to the number of all ship routes passing through a certain water area within a specified time period. It is a basic quantity that characterizes the actual situation of ship traffic in a certain water area.

[0057] Example 1

[0058] 1. Design of spherical hexagonal grid coding scheme

[0059] The planar "Hexagon Lattice Symmetry Tree" (HLST) scheme adopted by the present invention is based on a four-hole hexagonal grid (i.e., the area of ​​adjacent hierarchical grid cells is four times the relationship), and the grid cells are divided by hierarchical recursion, and {00, 01, 02, 03, 04, 05, 06, 10, 20, 30, 40, 50, 60} is used as the basic code element for encoding. Figure 1 As shown, the n+2th (n≥1)th layer grid cells and corresponding codes are jointly determined by the nth layer and the n+1th layer grid cells.

[0060] The HLST scheme implements coded addition operations, coded multiplication operations, proximity relationship queries, hierarchical relationship queries, and conversion between coded and plane Cartesian coordinates. Among them, the coded addition operation is defined based on the parallelogram rule, and its geometric meaning corresponds to the grid unit translation operation, such as Figure 2 shown.

[0061] The planar "Hexagon Lattice Symmetry Tree" scheme (HLST) is extended to the sphere using Snyder projection and regular icosahedron, and the mutual conversion between the longitude and latitude of any point on the sphere and the grid unit code is realized by means of projection forward and reverse calculation and coding and plane Cartesian coordinate conversion, such as Figure 3 shown.

[0062] 2. Spherical hexagonal grid vector line generation algorithm based on vector tracking

[0063] (1) Plane hexagonal grid vector line generation algorithm based on vector tracking

[0064] The adjacent consistency of the hexagonal grid ensures that the grid cells are adjacent to each other. Six "direction vectors" can be formed by connecting any grid cell and its adjacent cell center. For any vector line with the grid cell center as the starting point in the plane hexagonal grid, there will always be two "direction vectors" closest to its extension direction, which are called "adjacent vectors", such as Figure 4(a) As shown. Based on the vertical distance from the center of the grid unit to the vector line, the vector line can be tracked using the "near vector" to generate the vector line grid result, as shown in Figure 4 (b) is shown. This algorithm is called “planar hexagonal grid vector line generation algorithm based on vector tracking”.

[0065] (2) Extension of the algorithm to the spherical hexagonal grid system

[0066] like Figure 5 As shown in the figure, the Snyder projection is severely deformed in some areas, which may cause distortion of the regular grid on the plane after being projected onto the sphere, causing the grid originally arranged in a fixed direction to deviate. In this case, the spherical vector line may not be correctly tracked using only two "near vectors" according to the plane algorithm. Therefore, for the spherical vector line, three "near vectors" with an angle less than 180 degrees with the sphere are selected for tracking, such as Figure 6 shown.

[0067] Projection deformation causes the length of the spherical connection line between the centers of adjacent grids to be inconsistent, and it is impossible to track the spherical vector line with a fixed "near vector" as in the plane algorithm. The "near vector" tracking in the plane algorithm is equivalent to the translation of the grid unit, and the geometric meaning of the grid code addition operation is the grid unit translation operation. The translation operation realized by the code addition operation only retains the topological relationship between the grid units, but shields the concept of the distance between the unit centers, and the basic code elements {01, 02, 03, 04, 05, 06} just correspond to the extension direction of the six adjacent grid units. For example, the grid unit coded as Code performs the code addition operation Code+01, which is equivalent to translating it by one unit in the direction corresponding to the basic code element 01. With this feature, the three "near vectors" on the spherical surface can be replaced by the corresponding basic code elements {basecode1, basecode2, basecode3} in {01, 02, 03, 04, 05, 06}, so as to realize the vector tracking in the spherical hexagonal grid with the help of code addition operation.

[0068] like Figure 7 As shown, let the radius of the sphere be R, and the vertical distance d from the center of the grid unit V to the plane determined by the endpoint AB of the vector line and the center of the sphere O is VQ The spherical distance d from the cell center to the vector line VN The following conversion relationship exists:

[0069]

[0070] From the above relationship, we can see that the spherical distance d VN The perpendicular distance d to the point VQ There is a positive correlation between them, then d VQ Instead of dVN As a condition for determining the spherical grid unit, the computational complexity is reduced and the computational efficiency is improved.

[0071] In summary, the algorithm flow of spherical hexagonal grid vector line generation based on vector tracking is as follows: Figure 8 Selecting some regional administrative boundaries as vector line data, the algorithm is verified in the 5th and 6th grid layers. The results are shown in Fig. 9 shown.

[0072] 3. Ship traffic calculation based on spherical hexagonal grid vector line generation algorithm

[0073] (1) Ship route data processing and database storage structure

[0074] The situation where data is too dense often causes computational redundancy and thus affects efficiency, so the ship route data needs to be compressed. The present invention adopts a route data compression method based on square tiles (Varlamis, Iraklis Kontopoulos, Ioannis Tserpes, Konstantinos Etemad, Mohammad Soares, Amilcar Matwin, Stan. Building navigation networks from multi-vessel trajectory data [J]. Geoinformatica: An international journal of advances of computer science for geographic, 2021, 25 (1).), and divides the waters where the route data is located into square tiles at the same longitude and latitude intervals. When continuous track points fall in the same tile, only one representative track point is retained. This method is particularly effective for processing dense data in nearshore and stopover areas. The tile interval used in this embodiment is 0.01°, which has been proven to compress the volume of ship route data by about 5 to 10 times.

[0075] The processed ship route data is stored in the database (using SQLITE database), and the mmsi, lon, lat, dt_pos_utc (timestamp of each track point in utc time in timestamp form) fields are extracted from the data. The above fields are screened and the maximum value is extracted to form the TraceItem table in the database. Each row in the TraceItem table stores a set of basic information of ship route data. The table structure is shown in Table 1. The basic information of each set of route data can be used to achieve preliminary screening of ship route data. When the data is stored in the TraceItem table, the processed ship route data is directly stored in the database in the form of a table. The table name is named in the form of "mmsi_starttime_endtime" (such as 305530000_1487173125_1487666170), and stored in the origin_data_tablename field of the corresponding entry of this route data in the TraceItem table for original data query. At the same time, the mmsi_starttime_endtime corresponding to each ship route is different, so mmsi_starttime_endtime can also be used as a unique identifier for each ship route.

[0076] Table 1. TraceItem table structure

[0077]

[0078]

[0079] According to the lon, lat, and dt_pos_utc fields in the processed ship route data, the grid data corresponding to the ship route is generated using the spherical hexagonal grid vector line generation algorithm. The higher the grid level, the smaller the grid unit, and the higher the accuracy of the vector line fitting. This embodiment generates the 18th layer of grid data corresponding to the ship route, and the average spherical side length of the grid unit is about 20 meters, which can meet the accuracy requirements of actual applications. The ship route data is composed of sequentially arranged vector lines, and only the endpoints of each vector line have geographic coordinates and timestamp attributes. The gridded result generated thereby has no time information except for the two grid units at the endpoints. Although the projection may cause a large deformation of the grid units in the local spherical area, most of the spherical grids basically maintain their original form. Therefore, the gridded result of the ship route can be regarded as an approximately regular arrangement, and the timestamp information of other grid units is obtained by averaging the timestamps of the endpoint grid units. In this way, it is possible to ensure that each grid unit has time information, refine the time granularity of the ship route, and thus improve the accuracy of ship traffic calculation.

[0080] The 18th layer gridding result corresponding to the ship route is stored in the database in the form of a table. The table name is named in the form of "mmsi_starttime_endtime_18" (such as 305530000_1487173125_1487666170_18). The table structure is shown in Table 2. The above steps can complete the operation of storing the ship route data. The specific process is as follows: Fig.10 shown.

[0081] Table 2: Ship route grid data table structure

[0082] Field Name describe id Data item index Level Grid level Code Grid cell coding Center_lat Grid cell center latitude Center_lon Grid cell center longitude Timestamp Grid cell UTC timestamp

[0083] (2) Determination of grid cells in the buffer zone and calculation of ship traffic volume

[0084] The circular buffer used in the present invention is jointly determined by the geographical coordinates of the center point and the radius, and the grid unit whose center falls within the buffer zone is regarded as the buffer zone inclusion unit (i.e., the grid unit within the buffer zone). The calculation result of the ship traffic volume is obtained by matching the grid unit code within the buffer zone with the grid unit code in the grid data table corresponding to the route in the database. The higher the grid level, the more grid units are contained in the buffer zone, the higher the fitting accuracy of the grid unit for the circular buffer zone, and the higher the accuracy of traffic volume calculation, but the calculation and matching amount will also increase significantly. Therefore, it is crucial to adaptively determine the grid level and the corresponding level grid units in the buffer zone according to the buffer zone radius while taking into account both accuracy and computational efficiency.

[0085] like Fig.11 As shown, the strategy adopted in this embodiment is to first determine the nth layer of grids whose side length is just less than the radius of the buffer zone, and regard the n+1th layer of grids as a suitable grid based on this. The cell where the buffer zone center is located in the n+1th layer of grids and its surrounding adjacent cells, a total of 19 grid cells, are determined through proximity query and coded addition operations. The grid cells in the buffer zone are determined based on the spherical distance from the center of each grid cell to the center of the buffer zone. The grid cells whose spherical distance from the center of the grid cell to the center of the buffer zone is less than the radius of the circular buffer zone are grid cells in the circular buffer zone.

[0086] The buffer zone radius range with the n+1th layer as the applicable grid level is between the radius of the nth layer grid unit and the radius of the n-1th layer grid, such as Fig.12 As shown. The two extreme radius cases are further discussed, as shown in Fig.13 As shown in the figure, in the case of extreme radius, when the center of the buffer zone is at special positions such as the center of the grid unit, the vertex of the grid unit, and an edge of the grid unit, it can still ensure that the coverage area of ​​the grid unit in the buffer zone overlaps with most of the buffer zone, and the non-overlapping area will not be too concentrated, and the effect is relatively ideal.

[0087] The HLST scheme uses {00, 01, 02, 03, 04, 05, 06, 10, 20, 30, 40, 50, 60} as the basic code element for encoding. In the actual algorithm implementation, the basic code element can be replaced by the hexadecimal number {0, 1, 2, 3, 4, 5, 6, A, B, C, D, E, F}, then the n-th layer grid unit corresponds to the n-bit code. The grid coding scheme is designed according to the grid hierarchy, so it can reflect the hierarchical relationship of the grid system, and the code is The nth grid unit can be encoded by its first m (m < n) bits The corresponding m-th grid unit is divided. Based on this, it is determined that the buffer contains the ω (w≤18)th grid code w After that, you can check whether there is a code in the grid data table corresponding to each route in the database. ω The code at the beginning of the buffer zone is used to determine whether the ship route passes through the buffer zone.

[0088] In summary, after the ship route data is stored, given the geographical coordinates of the buffer center point, radius and query time interval, the geographical coordinates of the buffer center point and the query time interval can be used to preliminarily filter out the qualified ship route set T according to the max_lon, min_lon, max_lat, min_lat, starttime, and endtime fields in the TraceItem table: {trace1trace2, trace3, ....trace u}, where each ship route is identified by mmsi_starttime_endtime. Further, the applicable hierarchical grid unit code set C in the area is determined according to the buffer radius: {code1, code2, code3, ..., code m}(m≤19), trace each route in set T i (i≤n) corresponds to the 18th layer grid data (stored in the mmsi_starttime_endtime_18 table corresponding to each route) and each grid code code in set C m (m≤19) to query and match, if trace i (i≤n) The grid data corresponding to the grid contains code m (m≤19) and the corresponding timestamp is in the query time interval, the traffic volume L is accumulated. After all queries are completed, L is the traffic volume calculation result. The specific process is as follows Fig.14 shown.

[0089] Example 2 Application Example

[0090] The ship route data used in this embodiment is obtained by using the route extraction algorithm based on the starting and ending ports of the tanker AIS data in February 2024 (reference to (Guo X, Wang N, Ren Y, et al. Ship trajectory segmentation by movement states while addressing uncertainty and sparsity[J]. Ocean Engineering, 2024, 312. DOI: 10.1016 / j.oceaneng.2024.119218.) Ship trajectory segmentation by movement states while addressing uncertainty and sparsity[J]. Ocean Engineering, 2024, 312. DOI: 10.1016 / j.oceaneng.2024.119218.) The experimental area and route distribution are as follows: Fig.15 As shown, this area contains a total of about 8,300 route data. The spherical hexagonal grid vector line generation algorithm based on vector tracking is written in C++ language, the development tool is VisualStudio 2022, and the code is compiled as Release version; ship route data processing, warehousing and ship traffic volume calculation are written in Python language, and the development tool is PyCharm Community Edition 2023. The program is tested on a desktop computer (hardware configuration: Intel (R) Core (TM) i7-7700 CPU @ 3.60GHz, 16G RAM; operating system: Windows 10 Enterprise Edition). The experiment verifies the method of the present invention from the two aspects of ship traffic calculation efficiency and accuracy.

[0091] When the buffer and query time intervals are the same, the total amount of data is different, and the number of routes initially screened will also be different, which may affect the calculation time of ship traffic. The query time interval is set to 15 days from 00:00:00 on February 1, 2024 to 00:00:00 on February 16, 2024. The buffer radius is 100 km. Eight points are selected as the buffer center points near the straits and port areas where ship routes are relatively dense in the experimental area. The geographic coordinates are Loc1 (5.57°W, 35.94°N), Loc2 (2.86°W, 51.43°N), Loc3 (43.39°E, 12.59°N), Loc4 (36.13°E, 36.69°N), Loc5 (55.79°E, 26.36°N), Loc6 (102.43°E, 1.83°N), Loc7 (29.05°E, 41.09°N), and Loc8 (12.96°E, 36.83°N). The point distribution is as follows: Fig.15 The average ship traffic calculation time of the above 8 locations under different data totals is recorded, and the results are as follows Fig.16 shown.

[0092] Under the premise that the total amount of data remains unchanged, different buffer radius will result in different corresponding grid levels and the number of grid cells contained, and the change of query time interval will lead to different numbers of ship routes obtained by preliminary screening, both of which may also affect the traffic volume calculation time. The above 8 locations are still selected, with 00:00:00 on February 1, 2024 as the starting time, and the average ship traffic volume calculation time under different buffer radius and query time interval is counted. The results are shown in Table 3.

[0093] Table 3. Traffic volume calculation time statistics for different radius and query time intervals

[0094]

[0095] Depend on Fig.16 As shown in the results of Table 3, in general, the calculation time of ship traffic volume is positively correlated with the total amount of data, query time interval and radius. This is mainly because when the query conditions remain unchanged, the increase in the total amount of data and the increase in the query time interval will lead to an increase in the number of routes initially screened, so the database query and matching time will also increase. As for the buffer radius, although the buffer radius is not used in the process of preliminary screening of routes, when determining the grid units in the buffer zone, the buffer radius is needed as a basis for judgment. The increase in the buffer radius may cause an increase in the number of grid units in the buffer zone, which will also cause an increase in the database query and matching time. At the same time, it can be seen from the above results that the method proposed by the present invention has a high efficiency, and with the increase in the total amount of data, the query time interval and the buffer radius, the calculation time of ship traffic volume maintains a stable upward trend, and there will be no surge in calculation time, which basically meets the actual application needs.

[0096] Furthermore, the accuracy of the method of the present invention is verified. Intersection analysis is a method for calculating the geometric intersection of input elements. It can be used to accurately obtain the ship routes passing through the circular buffer zone. Most of the current mainstream commercial geographic information system software has built-in intersection analysis modules to enhance its spatial analysis function, such as Fig.17 As shown. Therefore, the intersection analysis results can be used as a reference standard to verify the accuracy of the method of the present invention. The above 8 locations are still selected, and the query time interval is set to 30 days from 00:00:00 on February 1, 2024 to 00:00:00 on March 2, 2024. The traffic volume calculation results of the method in this paper and the corresponding intersection analysis results are recorded under different buffer zone radii, as shown in Table 4. At the same time, the absolute value of the difference between the traffic volume calculation result of the method of the present invention and the intersection analysis result and the ratio of the traffic volume calculation result of the method of the present invention are defined as the deviation rate, which is used to characterize the accuracy of the method of the present invention. The lower the deviation rate, the higher the accuracy of the traffic volume calculation. The average accuracy of the 8 locations under different buffer zone radii is shown in the figure. Fig.18 shown.

[0097] Table 4. Comparison of traffic volume calculation results

[0098]

[0099] From the results in Table 4, it can be seen that in most cases, there is a difference between the traffic volume calculation results of the method of the present invention and the accurate value. The main reason for the difference is that, under the condition of taking into account both efficiency and accuracy, it is impossible to achieve high-precision fitting of the circular buffer zone with a certain number of discrete grid units. There must be non-overlapping areas between the grid units and the circular buffer zone, such as Fig.19 As shown in the red area, the routes that only pass through non-overlapping areas are the reason for the difference in traffic volume. In fact, using discrete grids instead of continuous vector data for calculation and analysis will inevitably lead to a loss in calculation accuracy, which is unavoidable. Fig.18 It can be seen that for the 8 areas selected in the experiment, the highest deviation rate is 14.34%, the lowest is 2.00%, and the overall average deviation rate is about 10.30%, which can basically meet the requirements of practical applications such as hot spot area determination, ship route planning, and ship route trend analysis.

[0100] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A traffic flow calculation method using spherical hexagonal grid vector lines, characterized in that: include: S1: Generate grid data corresponding to the ship route from the compressed ship route data using a spherical hexagonal grid vector line generation algorithm, and store it in a database; The spherical hexagonal grid vector line generation algorithm includes: A: Plane hexagonal grid vector line generation based on vector tracking: Based on the vertical distance from the center of the plane hexagonal grid unit to the vector line, the vector line is tracked using the adjacent vector to generate a vector line grid result; the hexagonal grid is divided into grid units in a hierarchical recursive manner, and {00, 01, 02, 03, 04, 05, 06, 10, 20, 30, 40, 50, 60} is used as the basic code element for encoding; B: The planar algorithm is extended to the spherical hexagonal grid system using Snyder equal-area projection. For the spherical vector line, three adjacent vectors whose angles with the spherical surface are less than 180 degrees are selected for tracking; the three adjacent spherical vectors are replaced by the corresponding basic code elements {basecode1, basecode2, basecode3} in {01, 02, 03, 04, 05, 06}, and the coded addition operation is used to realize the vector tracking in the spherical hexagonal grid; S2: Determine the grid cells in the buffer zone, and obtain the ship traffic volume by matching the grid cell codes in the buffer zone with the grid cell codes in the grid data corresponding to the ship routes in the database.

2. The method according to claim 1, characterized in that In S1, the compression is: using a route data compression method based on square tiles, dividing the waters where the route data is located into square tiles at the same longitude and latitude intervals, and when continuous track points fall in the same tile, only one representative track point is retained.

3. The method according to claim 1, characterized in that S1 includes: S1.1: compress the ship route data; S1.2: Store the compressed ship route data in the database in the form of a table; extract the key fields in the compressed ship route data to form entries and store them in the TraceItem table in the database; use the spherical hexagonal grid vector line generation algorithm to generate grid data corresponding to the ship route from the compressed ship route data, and store it in the database in the form of a table.

4. The method according to claim 3, characterized in that The key fields include mmsi, lon, lat, and dt_pos_utc fields.

5. The method according to claim 1, characterized in that Step B includes: B-1: Convert the spherical vector line endpoint coordinates (lon1, lat1) (lon2, lat2) to grid codes (Code1, Code2); lon represents longitude, lat represents latitude, and Code represents grid code; B-2: Determine the basic code element {basecode1, basecode2, basecode3} of the corresponding direction; Basecode represents the basic code element; B-3: Determine the spherical grid cells one by one using coded addition operations and point-to-surface vertical distances i , output spherical grid result {Grid1,Grid2,Grid 3。。。。。。 ,Grid n }; Grid represents a grid unit; the vertical distance between a point and a plane refers to the vertical distance d from the center V of the grid unit to the plane determined by the endpoint AB of the vector line and the center O of the sphere VQ .

6. The method according to claim 5, characterized in that Assume that the radius of the sphere is R, and the vertical distance d from the center of the grid unit V to the plane determined by the endpoint AB of the vector line and the center of the sphere O is VQ The spherical distance d from the cell center to the vector line VN The following conversion relationship exists:

7. The method according to claim 3, characterized in that S2 include: S2.1: In the TraceItem table, according to the buffer center coordinates, radius, and query time interval, select the ship routes that meet the conditions to form a set T: {tace1, tace2, tace 3…… ,tace n }; trace indicates the selected ship route; S2.2: Determine the grid cell code set C in the buffer zone according to the buffer zone radius: {code1, code2, code 3…… ,code m }, m≤19; S2.3: Set each route in set T to i Corresponding to the 18th layer of grid data and each grid code in set C m Perform query and match, if tace i The corresponding grid data contains code m The grid cells at the beginning and whose corresponding timestamps are in the query time interval, the traffic volume L is accumulated. After all queries are completed, L is the traffic volume calculation result; Among them, i≤n; m≤19.

8. The method according to claim 1, characterized in that The buffer zone is determined by the geographical coordinates of the center point and the radius, and the grid cell whose center falls within the buffer zone is the grid cell within the buffer zone.

Images