Polar region ship navigation planning method based on sea ice prediction and multi-stage path optimization

Abstract

The invention discloses a polar region ship navigation planning method based on sea ice prediction and multi-stage path optimization, and the method comprises the steps: obtaining a navigation risk map based on a to-be-navigated ship according to the updated single-day sea ice data in a time group; the method comprises the steps of obtaining a candidate path set based on a K shortest path algorithm, determining final U candidate paths according to path similarity of a result set of the candidate paths, further obtaining comprehensive scores of the paths in the result set, and finally determining a final driving path of a to-be-sailed ship according to the comprehensive score of each final candidate path. And completing the path planning of the polar region ship navigation. According to the method, the time, the risk and the path smoothness under the distance suboptimal constraint are taken as optimization targets, and similarity calculation is performed on the paths in the result set of the obtained candidate paths, so that diversified risk avoiding scheme paths can be provided. According to the method, the prediction time range is sliced, so that the high-frequency environment updating requirement in the prior art is met.

Description

Technical Field

[0001] The present invention relates to the field of polar navigation safety technology, and in particular to a polar ship navigation planning method based on sea ice prediction and multi-stage path optimization. Background Art

[0002] With global warming and the accelerated melting of Arctic sea ice, the commercial navigation windows of the Northeast Passage (NSR) and the Northwest Passage (NWP) are extended year by year, and research on polar navigation technology has gradually focused on dynamic risk modeling and path planning methods.

[0003] Existing technologies include the Polar Operation Limit Assessment Risk Indexing System (POLARIS), which provides a theoretical framework for risk classification by quantifying the match between ice conditions and ship performance; and generates the shortest route by simulating sea ice data, with distance as the optimization target.

[0004] Although existing technologies have made some progress, the following key bottlenecks still exist:

[0005] The POLARIS system is not deeply coupled with dynamic path planning, and risk assessment results cannot guide navigation decisions in real time. Furthermore, existing path planning technologies mostly focus on minimizing the cost of a single path, lacking attention to controlling the diversity of the path set and ensuring the smoothness of the candidate path segment structure. Furthermore, the route planning process does not incorporate auxiliary information from multiple source layers (such as user priority zones and ecological zones), resulting in a lack of human-machine collaboration. These factors, in the complex polar environment, can easily lead to large local fluctuations in the execution path, poor environmental response, and overfitting of routes, limiting the system's generalization and safety performance. Summary of the Invention

[0006] The present invention discloses a polar ship navigation planning method based on sea ice prediction and multi-stage path optimization to overcome the above technical problems.

[0007] In order to achieve the above object, the technical solution of the present invention is:

[0008] A polar ship navigation planning method based on sea ice prediction and multi-stage path optimization includes the following steps:

[0009] S1: Obtain historical sea ice data for the ice area, and obtain predicted daily sea ice data for the next P days based on the ConvLSTM model; wherein the sea ice data includes sea ice thickness data and sea ice density data;

[0010] S2: Time-slicing the next P days to obtain multiple time groups consisting of N days, so as to update the single-day sea ice data in the time groups according to the predicted sea ice data for the next P days;

[0011] S3: Establishing a grid environment map of the ice area to obtain a comprehensive navigation risk value of each grid cell for the ship to be sailed based on the single-day sea ice data in the updated time group, and generating a navigation risk map corresponding to the single-day sea ice data in the updated time group to obtain a navigation risk map for the ship to be sailed in the time group;

[0012] S4: According to the navigation risk map of the ships to be sailed in the time group, based on the K shortest path algorithm, obtain the route from the starting point v0 to the end point v M A set of U candidate paths arranged in ascending order of path distance; to obtain the path node sequence v0,v1,…,v of the uth,u=1,2,…U candidate paths m ,…,v M}; and add the uth, u=1, 2, …Uth candidate path to the result set of candidate paths; where v m Indicates the mth path node of the candidate path; m represents the index number of the path node;

[0013] S5: According to the Yen algorithm, obtain the first path node v1 to the end point v of the u, u=1, 2, ... U candidate paths. M The shortest path of , to obtain the updated path of the u-th candidate path;

[0014] If the updated path is different from the path in the result set, then the updated path is added to the result set of the candidate paths;

[0015] S6: If the number of paths in the candidate path result set is less than the set threshold, then the mth path node v of the uth, u=1, 2, ... Uth candidate paths are sequentially selected. m , m=2,3,…, repeat S5 until the number of paths in the result set of the candidate paths reaches the set threshold;

[0016] S7: Obtain the similarity of the paths in the result set of the candidate paths; if there is a τth path and a τth path in the result set of the candidate paths, If the similarity of the paths is less than the set dissimilarity threshold δ, the τth path or the τth path will be randomly The path is removed from the result set of candidate paths and based on the m+1th path node v of the candidate path m+1 , re-execute S5;

[0017] S8: When the result set of candidate paths does not contain the τth path and the When the similarity of the paths is less than the set dissimilarity threshold δ, the final set of candidate paths is obtained;

[0018] S9: Obtaining the comprehensive path scores of the candidate paths in the final candidate path result set to determine the final driving path of the ship to be navigated, so as to complete the path planning for polar ship navigation.

[0019] Furthermore, the method for obtaining the comprehensive path score of the candidate path is as follows:

[0020] S91: Get the total cost of the path;

[0021] S92: Get the structural smoothness penalty of the path;

[0022] S93: Obtain a comprehensive score of the path according to the total cost of the path and the structural smoothness penalty term of the path.

[0023] Furthermore, the formula used to obtain the total cost of the path is as follows:

[0024] C path =α1T total +α2C total

[0025]

[0026] Where: C path represents the total cost of the path; α1 represents the weight of the cumulative navigation time of the path; α2 represents the weight of the cumulative cost value of the path; T total Indicates the cumulative navigation time of the path; C total Indicates the cumulative cost of the path; m indicates the index number of the path node; M indicates the total number of path nodes; Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 Duration; C m Represents the total cost of the mth path node;

[0027] The formula for obtaining the total cost of the mth path node is as follows:

[0028]

[0029] Where: C m represents the total cost of the mth path node; n represents the index number of the layer; N represents the total number of layers; Represents the cost value at the mth path node on the nth layer; γ n Represents the weight of the nth layer; RIO m represents the comprehensive navigation risk value of the mth path node;

[0030] Get the ship from the mth path node v m Travel to the m+1th path node v m+1 The formula used for the duration is as follows:

[0031]

[0032] Where: Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 duration of use; Represents the mth path node v m and the m+1th path node v m+1 The path distance between them; Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 Speed m Indicates that the ship is at the mth path node v m Speed m+1 Indicates that the ship is at the m+1th path node v m+1 The speed of the ship.

[0033] Furthermore, the formula used to obtain the structural smoothness penalty term of the path is as follows:

[0034] P angle =∑ m ω m max(0,Δθ m -θ max )

[0035] Where: P angle represents the path structure smoothness penalty term; ω m represents the weight of the mth path; Δθ m represents the angle between the mth path segment and the m+1th path segment; θ max Indicates the steering angle threshold;

[0036] in,

[0037]

[0038] Where: v → m represents the direction vector of the mth path; ||v → m|| represents v → The magnitude of the m vector.

[0039] Furthermore, the formula used to obtain the comprehensive score of the path is as follows:

[0040] S total =β·C path +λ·P angle

[0041] Where: S total Indicates the comprehensive score of the path; C path is the total proxy value of the path; P angle is the penalty term for path structure smoothness; β represents the weight of the total path cost; λ represents the weight of the penalty term for path structure smoothness.

[0042] Furthermore, the method used to update the single-day sea ice data within the time group is as follows:

[0043] If, within the time group, only the single-day sea ice concentration data of the nth day exists as the maximum single-day sea ice concentration data within the time group, then each single-day sea ice data within the time group is updated to the sea ice data of the nth day;

[0044] If within the time group, the single-day sea ice concentration data on day n1 and day n are both the largest single-day sea ice concentration data within the time group, then

[0045] If the single-day sea ice thickness data of the nth day is greater than the single-day sea ice thickness data of the n1th day, then each single-day sea ice data in the time group is updated to the sea ice data of the nth day;

[0046] If the daily sea ice thickness data for day n is equal to the daily sea ice thickness data for day n1, then if day n is later than day n1, each daily sea ice data in the time group is updated to the sea ice data for day n; where n and n1 both represent the index of the day in the time group.

[0047] Furthermore, the formula used to obtain the similarity of the paths in the result set of candidate paths is as follows:

[0048]

[0049] in: represents the τth path and the Similarity of paths; G(p τ ) represents the set of path nodes covered by the τth path; Indicates the The set of path nodes covered by the paths; |·| represents the absolute value; τ, Both are index numbers of the paths in the result set of candidate paths.

[0050] Beneficial Effects: The present invention proposes a polar ship navigation planning method based on sea ice prediction and multi-stage path optimization. Using a ConvLSTM model, the method obtains predicted sea ice data for the next P days. Furthermore, based on the updated single-day sea ice data within a time group, a navigation risk map for the vessel to be sailed is generated. Using a K-shortest path algorithm, a set of candidate paths is obtained. A final set of U candidate paths is determined based on the path similarity of the resulting set of candidate paths. A comprehensive score is then obtained for each of the paths in the resulting set. Finally, based on the comprehensive score of each final candidate path, the final route for the vessel to be sailed is determined, completing path planning for polar ship navigation. This method addresses the prior art problem of using distance as a single optimization objective, ignoring time constraints and the cumulative risk effect, which can easily lead to locally high-risk paths. By optimizing time, risk, and path smoothness under suboptimal distance constraints, the method calculates similarity for the resulting set of candidate paths to provide diverse risk-avoidance solutions. Furthermore, the method addresses the need for frequent environmental updates in the prior art by slicing the predicted time range and reducing complexity through phased local optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative labor.

[0052] Figure 1 This is a flow chart of the polar ship navigation planning method based on sea ice prediction and multi-stage path optimization of the present invention;

[0053] Figure 2 Schematic diagram of the technical route in an embodiment of the present invention. DETAILED DESCRIPTION

[0054] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0055] This embodiment introduces a polar ship navigation planning method based on sea ice prediction and multi-stage path optimization, including the following steps: Figure 1 and Figure 2 As shown:

[0056] S1: Obtain historical sea ice data of the ice zone (Arctic region) and obtain predicted sea ice data for the next P days (30 days) based on the ConvLSTM model; wherein the sea ice data includes sea ice thickness data and sea ice density data;

[0057] Specifically, in this embodiment, the ice area in the Arctic region is taken as the research object, and future sea ice data is predicted based on historical sea ice data, specifically including sea ice density data and sea ice thickness data. The historical sea ice density and sea ice thickness data are from the CMEMS platform, with a spatial coverage range of 52°N~85°N, -180°~179.96°. The time range is 765 days of daily data from August 1 to December 31 each year from 2020 to 2024, with a spatial resolution of 3×3km. Among them, 612 days of data from August 1 to December 31 each year from 2020 to 2023 are used for model training, and 153 days of data in 2024 are used for model testing and effect verification. In addition, because there can be no singular values ​​in the convolution calculation, the Nan value of land is assigned to 0 during model training, that is, it is considered that there is no sea ice here. Based on the ConvLSTM model, the data from the first to the fifth day is used to predict the data from the sixth to the eighth day. Therefore, the prediction results for a total of 148 days from August 6 to December 31 of the same year can be obtained through calculation.

[0058] S2: Time-slicing the next P days to obtain multiple time groups consisting of N days, so as to update the single-day sea ice data in the time groups according to the predicted sea ice data for the next P days;

[0059] Specifically, the time in the time group for slicing is determined based on the set update time of the navigation risk map. In this embodiment, the predicted 30 days are sliced ​​into 3-day segments to obtain multiple 3-day time groups, so as to determine the single-day sea ice data within the time group based on the predicted sea ice data for the next 30 days.

[0060] Preferably, the method used to update the single-day sea ice data in the time group is as follows:

[0061] If, within the time group, only the single-day sea ice concentration data of the nth day exists as the maximum single-day sea ice concentration data within the time group, then each single-day sea ice data within the time group is updated to the sea ice data of the nth day;

[0062] If within the time group, the single-day sea ice concentration data on day n1 and day n are both the largest single-day sea ice concentration data within the time group, then

[0063] If the single-day sea ice thickness data of the nth day is greater than the single-day sea ice thickness data of the n1th day, then each single-day sea ice data in the time group is updated to the sea ice data of the nth day;

[0064] If the daily sea ice thickness data for day n is equal to the daily sea ice thickness data for day n1, then if day n is later than day n1, each daily sea ice data in the time group is updated to the sea ice data for day n; where n and n1 both represent the index of the day in the time group.

[0065] Specifically, in this embodiment, the single-day sea ice data with the highest sea ice density within the group is first selected. If there are multiple single-day sea ice data with the same and maximum sea ice density, the single-day sea ice data with the largest sea ice thickness is selected among the single days with the same sea ice density. If the sea ice density and sea ice thickness of all single days are the same, the single-day sea ice data of the last day of each group is selected by default.

[0066] In this embodiment, the predicted future 30 days are sliced ​​into 3-day segments, where the sea ice concentrations of the three days in a certain time group are C1, C2, and C3, and the sea ice thicknesses are T1, T2, and T3. If only the sea ice concentration C1 of the first day is the largest, the sea ice data of each day in the time group is replaced by the sea ice data of the first day. If the sea ice concentration C1 of the first day and the sea ice concentration C2 of the second day are equal and both are the maximum sea ice concentrations in the time group, the sea ice thicknesses of the first day and the second day are compared. If the sea ice thickness of the second day is greater than that of the first day, the sea ice data of each day in the time group is replaced by the sea ice data of the second day. If, at this time, the sea ice thickness of the second day is also equal to that of the first day, then since the second day is later than the first day, the sea ice data coverage of each day in the time group is replaced by the sea ice data of the second day.

[0067] S3: Establish a grid environment map of the Arctic region to obtain the comprehensive navigation risk value of the ship to be sailed in each grid cell based on the single-day sea ice data in the updated time group and the POLARIS rule, and generate a navigation risk map corresponding to the single-day sea ice data in the updated time group to obtain a navigation risk map based on the ship to be sailed in the time group;

[0068] Specifically, this embodiment introduces the Polar Operational Restrictions Assessment Risk Index System (POLARIS), which analyzes sea ice density and thickness data and combines them with ship ice class to construct a grid environment map with a navigation risk index.

[0069] Specifically, the Polar Operational Restrictions Assessment Risk Index System (POLARIS) first determines an Ice Navigation Risk Index (RIV) for each type of ship based on the ship's ice class and ice type (as shown in Table 1). This value indicates the severity of the impact of sea ice on the ship. The method used to determine the risk index is conventional for those skilled in the art and will not be described in detail here. A RIV < 0 indicates a high navigation risk, making the area unsuitable for navigation by that type of ship (i.e., an unnavigable area (obstacle)). A RIV > 0 indicates that the area is navigable, and higher values ​​indicate lower navigation risk.

[0070] Table 1 Definition of ice types

[0071]

[0072] The formula for calculating the integrated navigation risk (RIO) is:

[0073]

[0074] Where: RIO T Indicates the comprehensive navigation risk value of the ship in the T type sea ice cover area; T represents the ice type of the sea ice cover area; C T Indicates the density of T-type sea ice in the grid cell, with a value between [0,1]; The ice breaking level is ice class Risk value of ships sailing in T type sea ice cover area; ice class The icebreaking level of the ship to be sailed.

[0075] The pixel colors of the navigation area are divided according to the navigation risk levels corresponding to different ice classes of ships in different sea ice conditions, and the navigation risk levels are divided into:

[0076] (1) High-risk area: RIO < -10;

[0077] (2) Medium risk area: -10≤RIO<0;

[0078] (3) Low-risk area: RIO ≥ 0.

[0079] High-risk areas are displayed as red (255, 51, 0), medium-risk areas as yellow (255, 255, 0), and low-risk areas as green (146, 208, 80). In addition, areas without sea ice data, i.e., land areas, are displayed as black (0, 0, 0). The generated risk map varies depending on the ship's ice class. Ships cannot sail in high-risk areas or land areas, and are only allowed to sail in medium- and low-risk areas, i.e., yellow and green areas.

[0080] Table 2 Risk Index Value (RIV) of Ice Navigation for Various Types of Ships in the POLARIS System

[0081]

[0082] This embodiment combines the POLARIS rules to map the predicted data into multi-dimensional risk values ​​(ice thickness, density, and time decay risk), constructs a spatiotemporal coupled rasterized navigation risk map, and realizes dynamic quantitative assessment of complex ice zone threats.

[0083] S4: According to the navigation risk map of the ships to be sailed in the time group, based on the K shortest path algorithm, obtain the route from the starting point v0 to the end point v M U candidate sets arranged in ascending order of path distance; to obtain the path node sequence {v0,v1,…,v m ,…,v M}; and add the uth, u=1, 2, …Uth candidate path to the result set of candidate paths; where v m represents the mth path node of the candidate path; m represents the index number of the path node; among them, the ship starts from the mth path node v m Travel to the m+1th path node v m+1 The time is less than the time of each time group.

[0084] Specifically, based on the K shortest paths algorithm and according to the navigation risk map of the ship to be sailed, U candidate paths are generated and arranged in ascending order of path distance to form a candidate path set.

[0085] Specifically, a preliminary set of paths is generated based on the existing K shortest path algorithm, including:

[0086] First, use the A* algorithm to determine the initial shortest distance route, where the cost function of the A* algorithm is defined as:

[0087] f(n)=g(n)+h(n)

[0088] Where: g(n) is the cost function, which represents the cumulative cost from the initial path node to the current path node; h(n) is the heuristic function. This embodiment uses the Euclidean distance as the heuristic function, which represents the estimated cost from the current path node to the target path node; f(n) is the estimated cost from the initial path node through the current path node to the target path node.

[0089] S5: According to the Yen algorithm, obtain the first path node v1 to the end point v of the u, u=1, 2, ... U candidate paths. M The shortest path of , to obtain the updated path of the u-th candidate path;

[0090] If the updated path is different from the path in the result set, the updated path is added to the result set of candidate paths;

[0091] Specifically, when the first path node v1 to the end point v of the u-th candidate path is obtained M The shortest path of the u-th candidate path is spliced ​​with the sub-path between the starting point v0 and the first path node v1 of the navigation route to form a complete updated path of the u-th candidate path. If the updated path is different from the candidate in the result set, the updated path is added to the result set of the u-th candidate path;

[0092] S6: If the number of paths in the candidate path result set is less than the set threshold, then the mth path node v of the uth, u=1, 2, ... Uth candidate paths are sequentially selected. m , m=2,3,…, repeat S5 until the number of paths in the result set of the candidate paths reaches the set threshold;

[0093] S7: Get the similarity of the paths in the result set of candidate paths; if there is a τth path and a τth path in the result set of candidate paths, If the similarity of the paths is less than the set dissimilarity threshold δ, the τth path or the τth path will be randomly The path is removed from the result set of candidate paths and S5 is re-executed;

[0094] Specifically, after S6 is executed, based on the mth path node v m The number of paths in the result set of candidate paths reaches the set threshold. At this time, the path diversity maintenance mechanism is introduced to calculate the similarity of the paths in the result set of candidate paths. If there is a τth path and a If the similarity of the paths is less than the set dissimilarity threshold δ, the τth path or the τth path will be randomly The path is removed from the result set of candidate paths. Then, based on the m+1th path node v m+1 Repeat S5.

[0095] Preferably, the formula used to obtain the similarity of the paths in the result set of candidate paths is as follows:

[0096]

[0097] in: represents the τth path and the Similarity of paths; G(p τ ) represents the set of path nodes covered by the τth path; Indicates the The set of path nodes covered by the paths; |·| represents the absolute value; τ, The index number of the path in the result set of candidate paths; path τ represents the τth path; Indicates the path.

[0098] Specifically, this embodiment sets the dissimilarity threshold δ∈[0.2,0.4], if Consider the τth path and the If the paths have repeated structures, any redundant paths with similarities below the dissimilarity threshold are eliminated, and a set of paths with significant structural differences is output for subsequent selection. This embodiment uses a path dissimilarity screening mechanism to eliminate structurally similar paths from the candidate path set, retaining only paths with significant differences. This creates a diversified path set, increases user selection space and risk redundancy, and expands the selection space for resulting paths.

[0099] S8: When the result set of candidate paths does not contain the τth path and the When the similarity of the paths is less than the set dissimilarity threshold δ, the result set of candidate paths at this time is the final result set of candidate paths;

[0100] S9: Obtaining the comprehensive path scores of the candidate paths in the final candidate path result set to determine the final driving path of the ship to be navigated, so as to complete the path planning for polar ship navigation.

[0101] Specifically, this embodiment selects the candidate path with the smallest comprehensive score from the candidate path result set as the final travel path for the ship to be sailed, so as to optimize time, risk, and path smoothness under the suboptimal distance constraint, and complete the path planning for polar ship navigation.

[0102] Preferably, the method for obtaining the comprehensive path score of the candidate path is as follows:

[0103] S91: Get the total cost of the path;

[0104] C path =α1T total +α2C total

[0105]

[0106] Where: C path represents the total cost of the path; α1 represents the weight of the cumulative navigation time of the path; α2 represents the weight of the cumulative cost value of the path; T total Indicates the cumulative navigation time of the path; C total Indicates the cumulative cost of the path; m indicates the index number of the path node; M indicates the total number of path nodes; Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 Duration; C m Represents the total cost of the mth path node;

[0107] Specifically, when m=0, it represents the starting point v0 of the navigation route, and when m=M, it represents the end point v M ; Weights α1, α2∈[0,1], satisfying α1+α2=1.

[0108] Specifically, this embodiment introduces multi-dimensional auxiliary layer information into the cost calculation of path nodes, specifically including: allowing users to load several custom layers, such as ecologically sensitive areas, user preference areas, etc. The layer data uses a format compatible with the geographic information system (GeoTIFF), which can correspond to each node on the path raster map and use a numerical field to identify the preference level or restriction level of the area described by the node. The layer data is manually annotated by the user or imported from an external platform; each layer is assigned a weight γ n ∈R, R represents the set of real numbers; where γ n A positive value indicates encouragement, γ n A negative value indicates avoidance, and the node total cost correction function is as follows:

[0109]

[0110] Where: C m represents the total cost of the mth path node; n represents the index number of the layer; N represents the total number of layers; Represents the cost value at the mth path node on the nth layer; γ n Represents the weight of the nth layer; RIO m represents the comprehensive navigation risk value of the mth path node;

[0111] Get the ship from the mth path node v m Travel to the m+1th path node v m+1 The formula used for the duration is as follows:

[0112]

[0113] Where: Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 duration of use; Represents the mth path node v m and the m+1th path node v m+1 The path distance between them; Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 Speed, where the speed between two adjacent nodes is the minimum of the two speeds to ensure navigation safety; Speed m Indicates that the ship is at the mth path node v m Speed m+1 Indicates that the ship is at the m+1th path node v m+1 Speed ​​at 1 hour;

[0114] In this embodiment, by This constraint limits the total flight time from the subpath's starting point to its end point to less than or equal to 72 hours (the time of each time group);

[0115] Specifically, by traversing the navigation risk graph, each path node in the ship's path is found according to the path node risk level. m , where the mapping relationship between the sea ice risk value and the navigation speed associated with each node in the rasterized risk map is as follows:

[0116]

[0117] Specifically, in this embodiment, if the total time of the path nodes in a sub-path plus the last segment of the path exceeds 72 hours, the remaining time (Remain_Time) from the end of the previous segment to 72 hours is calculated, and then the maximum distance that can be navigated at the corresponding speed within this time period (Remain_Distance) is calculated to obtain the final stop point of the sub-path.

[0118]

[0119] Where: Indicates the starting point node of the path_i subpath; The end node of the path_i-1th subpath is set as the starting point of each subpath, and the end point of the previous subpath is set. When the end point of the subpath coincides with the global end point, the algorithm terminates. This ensures that the navigation time of each subpath does not exceed 72 hours, that is, the number of days in a time group, to ensure that the ship can complete the journey from the starting point to the end point when the navigation risk map is updated.

[0120]

[0121] Where: Remain Time Indicates the remaining allocatable time after deducting the used time from the maximum allowed time budget (72 hours); Remain_Distance indicates the remaining allocatable time after Remain TimeThe maximum remaining distance that can be traveled based on the speed of the current path segment;

[0122] Specifically, this embodiment employs a multi-stage path generation strategy (72-hour sub-path iteration) and a multi-criteria optimization method under suboptimal distance constraints. Using a K-shortest path algorithm, a candidate path set is generated, screening global routes with short time and manageable risks. This avoids the local suboptimal nature of single-objective optimization and achieves a global balance between long-range safety and efficiency. By using three-day data slices and dynamically correcting route deviations, the vessel can continuously adapt to the rapid evolution of sea ice during a 30-day voyage, reducing the risk of entrapment and enhancing autonomous navigation capabilities in complex ice areas.

[0123] This embodiment also uses a layer-guided path generation module to integrate multi-source auxiliary layer information (such as ecologically sensitive areas and user-preferred areas) into path cost calculation, enabling path structure control tailored to user preferences, ecological sensitivity, and regional authority. This enhances the environmental responsiveness of the generated path.

[0124] S92: Get the structural smoothness penalty of the path;

[0125] Specifically, the present embodiment detects the continuity of path segments and outputs the optimal path, specifically including: for each path in the candidate path, sequentially extracting the direction vectors v of its adjacent path segments; → m,v → m+1, calculate the angle Δθ m , if Δθ m >θ max , it is considered as a jump segment on the path, and in the comprehensive scoring of the path, a structural smoothness penalty term P is added. angle :

[0126] P angle =∑ m ω m max(0,Δθ m -θ max )

[0127] Where: P angle represents the path structure smoothness penalty term; ω m represents the weight of the mth path, ω m ∈[0.5,2.0]; Δθ m represents the angle between the mth segment path and the m+1th segment path; θ max Indicates the steering angle threshold;

[0128] Specifically, in this embodiment, the subpath from the starting point v0 of the navigation route to the first path node v1 is the first subpath, the subpath from the first path node v1 of the navigation route to the second path node v2 is the second subpath, and so on.

[0129] in,

[0130]

[0131] Where: v → m represents the direction vector of the mth path; ||v → m|| represents v → The modulus of the m vector;

[0132] This embodiment introduces a path segment splicing smoothness detection mechanism, which evaluates the continuity of inter-segment connections through vector angles to prevent path jumps and reversals, thereby improving the navigability of the path and the robustness of the path generation structure.

[0133] S93: Combine the total generation value and the structural penalty to obtain the comprehensive score of the path;

[0134] S total =β·C path +λ·P angle

[0135] Where: S total Indicates the comprehensive score of the path; C path is the total cost of the path (including time and risk); β represents the weight of the total cost of the path; P angle represents the path structure smoothness penalty term; λ represents the weight of the path structure smoothness penalty term; where β, λ∈[0,1], and β+λ=1;

[0136] Specifically, this embodiment is based on S total The scoring ultimately selects the optimal path with both smoothness and efficiency as the execution plan.

[0137] Specifically, in this embodiment, among the U candidate paths found, for each candidate path, a time constraint is used to limit the maximum sailing time of each sub-path on the current candidate path to within 72 hours: based on the rasterized risk map, the distance between adjacent nodes on the sub-path is obtained, and the navigation speed of the ship is determined according to the sea ice conditions at the node, thereby determining the navigation time of the ship, ensuring that the navigation time used for each sub-path remains within 72 hours.

[0138] This embodiment uses speed-risk mapping and a multi-source auxiliary layer weighting mechanism to calculate the comprehensive cost of each path and construct an initial set of candidate paths. Based on this set of candidate paths, the path similarity between any pair of paths is calculated. A difference threshold is set to filter out structurally diverse paths, thus constructing a candidate set for path diversity control.

[0139] This embodiment solves the problem in the prior art of using distance as a single optimization objective, ignoring time constraints and the cumulative effect of risk, which easily leads to local high-risk paths. This embodiment uses time, risk, and path smoothness under the suboptimal distance constraint as optimization objectives, and calculates similarity for the paths in the result set of obtained candidate paths to provide diversified risk avoidance solution paths. This embodiment solves the high-frequency environment update requirements in the prior art by slicing the predicted time range and reducing complexity through phased local optimization. At the same time, this embodiment uses the A* algorithm to generate the initial path, which can solve the problem of low search efficiency in high-resolution grid maps due to reliance on Dijkstra initialization, thereby improving the efficiency of path planning.

[0140] 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 above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A polar ship navigation planning method based on sea ice prediction and multi-stage path optimization, characterized in that: The steps include: S1: Obtain historical sea ice data for the ice area, and obtain predicted daily sea ice data for the next P days based on the ConvLSTM model; wherein the sea ice data includes sea ice thickness data and sea ice density data; S2: Time-slicing the next P days to obtain multiple time groups consisting of N days, so as to update the single-day sea ice data in the time groups according to the predicted sea ice data for the next P days; S3: Establishing a grid environment map of the ice area, obtaining a comprehensive navigation risk value of each grid cell for the ship to be sailed based on the single-day sea ice data in the updated time group, generating a navigation risk map corresponding to the single-day sea ice data in the updated time group, and obtaining a navigation risk map for the ship to be sailed in the time group; S4: According to the navigation risk map of the ships to be sailed in the time group, based on the K shortest path algorithm, obtain the route from the starting point v0 to the end point v M A set of U candidate paths arranged in ascending order of path distance; to obtain the path node sequence v0,v1,…,v of the uth,u=1,2,…U candidate paths m ,…,v M }; and add the uth, u=1, 2, …Uth candidate path to the result set of candidate paths; where v m Indicates the mth path node of the candidate path; m represents the index number of the path node; S5: According to the Yen algorithm, obtain the first path node v1 to the end point v of the u, u=1, 2, ... U candidate paths. M The shortest path to obtain the updated path of the u-th candidate path; If the updated path is different from the path in the result set, then the updated path is added to the result set of the candidate paths; S6: If the number of paths in the candidate path result set is less than the set threshold, then the mth path node v of the uth, u=1, 2, ... Uth candidate paths are sequentially selected. m , m=2,3,…, repeat S5 until the number of paths in the result set of the candidate paths reaches the set threshold; S7: Obtain the similarity of the paths in the result set of the candidate paths; if there is a τth path and a τth path in the result set of the candidate paths, If the similarity of the paths is less than the set dissimilarity threshold δ, the τth path or the τth path will be randomly The path is removed from the result set of candidate paths and based on the m+1th path node v of the candidate path m+1 , re-execute S5; S8: When the result set of candidate paths does not contain the τth path and the When the similarity of the paths is less than the set dissimilarity threshold δ, the final set of candidate paths is obtained; S9: Obtaining the comprehensive path scores of the candidate paths in the final candidate path result set to determine the final driving path of the ship to be navigated, so as to complete the path planning for polar ship navigation.

2. The polar ship navigation planning method based on sea ice prediction and multi-stage path optimization according to claim 1 is characterized in that: The method for obtaining the comprehensive path score of a candidate path is as follows: S91: Get the total cost of the path; S92: Get the structural smoothness penalty of the path; S93: Obtain a comprehensive score of the path according to the total cost of the path and the structural smoothness penalty term of the path.

3. The polar ship navigation planning method based on sea ice prediction and multi-stage path optimization according to claim 2 is characterized in that: The formula used to obtain the total cost of the path is as follows: C path =α1T total +α2C total Where: C path represents the total cost of the path; α1 represents the weight of the cumulative navigation time of the path; α2 represents the weight of the cumulative cost value of the path; T total Indicates the cumulative navigation time of the path; C total Indicates the cumulative cost of the path; m indicates the index number of the path node; M indicates the total number of path nodes; Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 Duration; C m Represents the total cost of the mth path node; The formula for obtaining the total cost of the mth path node is as follows: Where: C m represents the total cost of the mth path node; n represents the index number of the layer; N represents the total number of layers; Represents the cost value at the mth path node on the nth layer; γ n Represents the weight of the nth layer; RIO m represents the comprehensive navigation risk value of the mth path node; Get the ship from the mth path node v m Travel to the m+1th path node v m+1 The formula used for the duration is as follows: Where: Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 duration of use; Represents the mth path node v m and the m+1th path node v m+1 The path distance between them; Indicates that the ship starts from the mth path node v m Travel to the m+1th path node v m+1 Speed m Indicates that the ship is at the mth path node v m Speed m+1 Indicates that the ship is at the m+1th path node v m+1 The speed of the ship.

4. The polar ship navigation planning method based on sea ice prediction and multi-stage path optimization according to claim 2 is characterized in that: The formula used to obtain the structural smoothness penalty term of the path is as follows: P angle =∑ m oh m ·max(0,Δθ m -θ max ) Where: P angle represents the path structure smoothness penalty term; ω m represents the weight of the mth path; Δθ m represents the angle between the mth path segment and the m+1th path segment; θ max Indicates the steering angle threshold; in, Where: v → m represents the direction vector of the mth path; ||v → m|| represents v → The magnitude of the m vector.

5. The polar ship navigation planning method based on sea ice prediction and multi-stage path optimization according to claim 2 is characterized in that: The formula used to obtain the comprehensive score of the path is as follows: S total =β·C path +λ·P angle Where: S total Indicates the comprehensive score of the path; C path is the total proxy value of the path; P angle is the penalty term for path structure smoothness; β represents the weight of the total path cost; λ represents the weight of the path structure smoothness penalty term.

6. The polar ship navigation planning method based on sea ice prediction and multi-stage path optimization according to claim 1 is characterized in that: The method used to update the daily sea ice data within the time group is as follows: If, within the time group, only the single-day sea ice concentration data of the nth day exists as the maximum single-day sea ice concentration data within the time group, then each single-day sea ice data within the time group is updated to the sea ice data of the nth day; If within the time group, the single-day sea ice concentration data on day n1 and day n are both the largest single-day sea ice concentration data within the time group, then If the single-day sea ice thickness data of the nth day is greater than the single-day sea ice thickness data of the n1th day, then each single-day sea ice data in the time group is updated to the sea ice data of the nth day; If the daily sea ice thickness data for day n is equal to the daily sea ice thickness data for day n1, then if day n is later than day n1, each daily sea ice data in the time group is updated to the sea ice data for day n; where n and n1 both represent the index of the day in the time group.

7. The polar ship navigation planning method based on sea ice prediction and multi-stage path optimization according to claim 1 is characterized in that: The formula used to obtain the similarity of the paths in the result set of candidate paths is as follows: in: represents the τth path and the Similarity of paths; G(p τ ) represents the set of path nodes covered by the τth path; Indicates the The set of path nodes covered by the paths; |·| represents the absolute value; τ, Both are index numbers of the paths in the result set of candidate paths.

Images