Route optimization method based on particle swarm optimization and geographic information

Abstract

This invention discloses a route optimization method based on particle swarm optimization (PSO) and geographic information, specifically relating to the field of ship route optimization technology. First, multi-source geographic information of the target shipping area is collected and preprocessed to obtain a standard voyage operation dataset, and a two-layer control constraint model is constructed. A path search algorithm is used to pre-calculate the forward and backward operating cost matrices. Iterative optimization and verification are then performed using a PSO algorithm incorporating an operating path traction term, outputting the globally optimal first waypoint sequence. The globally optimal first waypoint sequence is smoothed using B-spline interpolation and then verified a second time, outputting the verified optimal route. This invention constructs a bidirectional feasible region to define the targeted search boundary of the PSO algorithm and innovatively introduces an operating path traction term into the velocity update formula, integrating geographic information throughout the process. This solves the problems of premature convergence and local optima caused by a lack of geographic guidance.

Description

Technical Field

[0001] This invention relates to the field of ship route optimization technology, specifically a route optimization method based on particle swarm optimization and geographic information. Background Technology

[0002] With the large-scale development of the shipping industry, the optimization and scheduling of ship routes is a core part of enterprise operation and management. In response to the route planning needs in complex geographical environments, it is necessary to strictly adhere to various operational constraints such as terrain obstacles, no-fly and no-navigation zones, weather conditions, and transit costs to achieve ship route operation management. At present, for the ship route optimization problem, the mainstream technologies mostly adopt intelligent optimization algorithms combined with geographic information to achieve the implementation path. Among them, particle swarm optimization algorithm is widely used in various path planning scenarios.

[0003] However, when traditional particle swarm optimization algorithms are applied to ship route optimization, the particle population is usually randomly initialized across the entire area. Most existing route optimization schemes based on particle swarm optimization use geographical information as a single constraint, simply distinguishing between passable and impassable areas. They mostly perform routine optimizations on the algorithm's own parameters, and simply discard or randomly reselect particles that fly out of the feasible area or touch hard constraints during iteration, thereby achieving the purpose of ship route optimization.

[0004] Existing route optimization schemes still have some limitations. For example, they do not deeply integrate operational geographic information into the algorithm iteration process, lack operational path guidance, and are prone to premature convergence and getting trapped in local optima. They cannot quantify the degree of matching between routes and geographic environments, and the final output routes may deviate from operational constraints and have excessively high costs, making it difficult to meet the high-precision and highly practical route optimization needs in complex geographic environments. In summary, existing route optimization technologies based on particle swarm optimization algorithms, due to defects such as coarse utilization of geographic information and insufficient coupling between algorithms and geographic constraints, cannot achieve route planning that balances compliance with hard operational constraints and optimal soft operational costs in complex geographic environments. There is an urgent need for a route optimization method that uses hierarchical adaptation of geographic information and targeted guidance from particle swarm search to solve many of the shortcomings of existing technologies. Summary of the Invention

[0005] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a route optimization method based on particle swarm optimization and geographic information to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a route optimization method based on particle swarm optimization and geographic information, comprising:

[0007] S1: Collect multi-source geographic information of the target voyage's operating area and preprocess it to obtain the voyage's standard operating dataset;

[0008] S2: Based on the standard dataset of voyage operations, construct a two-layer voyage operation control and constraint model, including a hard operation constraint layer and a soft operation cost layer;

[0009] S3: Based on the voyage start point, end point and hard operational constraints layer, the forward operational cost matrix and the backward operational cost matrix are pre-calculated using a path search algorithm and then fused to generate a bidirectional feasible region;

[0010] S4: Based on the bidirectional feasible region, initialize the particle swarm population, construct the fitness function based on the soft operating cost layer, and perform iterative optimization and verification by introducing the operating path traction term of the particle swarm algorithm, and output the globally optimal first waypoint sequence;

[0011] S5: The globally optimal first waypoint sequence is smoothed using B-spline interpolation. The smoothing result is then used for secondary verification, and the verified optimal route is output.

[0012] The technical effects and advantages of this invention are as follows:

[0013] 1. This invention delineates the target search boundary of the particle swarm by constructing a two-layer voyage operation management and control constraint model, and innovatively introduces an operation path traction term into the velocity update formula, so that the particles are not only guided by individual and global optima during the iteration process, but also continuously guided by geographical constraints, integrating geographical information throughout the process, and solving the problems of premature convergence and local optima caused by lack of geographical guidance.

[0014] 2. This invention constructs a soft operating cost layer that integrates multi-source geographic information and builds a geographic adaptation cost function based on this layer, thereby realizing a quantitative assessment of the degree of matching between flight routes and geographic environment. By incorporating "compliance" and "economy" into a unified optimization framework through the fitness function, it solves the problem of excessively high flight route costs caused by the inability to quantify the degree of adaptation in existing technologies.

[0015] 3. This invention achieves deep coupling between the algorithm search space and geographical constraints by using a bidirectional feasible region selected through bidirectional cost matrix fusion. This avoids the traditional coarse global search, improves the algorithm's convergence efficiency, and ensures that the final output route resolves the core contradiction of the disconnect between the algorithm and geographical constraints. It perfectly adapts to the engineering requirements for high-precision and high-practicability routes in complex geographical environments. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the overall process of the present invention.

[0017] Figure 2 This is a schematic diagram of the method flow of the present invention.

[0018] Figure 3This is a schematic diagram of the global optimal first waypoint sequence output process of the present invention. Detailed Implementation

[0019] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] Please see Figure 1 As shown, the present invention provides a route optimization system based on particle swarm optimization and geographic information, including a shipping operation data acquisition module, an operation constraint hierarchical modeling module, a bidirectional feasible domain construction module, a particle swarm iterative optimization module, and a route optimization output module.

[0021] The shipping operation data acquisition module is connected to the operation constraint hierarchical modeling module, the bidirectional feasible domain construction module is connected to both the operation constraint hierarchical modeling module and the particle swarm iterative optimization module, and the route optimization output module is connected to the particle swarm iterative optimization module.

[0022] Shipping operation data acquisition module: Collects multi-source geographic information of the target voyage's operating area and preprocesses it to obtain a standard dataset for voyage operations;

[0023] Operational constraint hierarchical modeling module: Based on the standard dataset of voyage operations, a two-layer voyage operation control constraint model is constructed, including a hard operational constraint layer and a soft operational cost layer;

[0024] Bidirectional feasible region construction module: Based on the voyage start point, end point and hard operational constraint layers, the forward operational cost matrix and the backward operational cost matrix are pre-calculated using a path search algorithm and then fused to generate a bidirectional feasible region;

[0025] Particle swarm iterative optimization module: Based on the bidirectional feasible region, the particle swarm population is initialized, the fitness function is constructed based on the soft operating cost layer, and the particle swarm algorithm with the introduction of the operating path traction term is used for iterative optimization and verification, and the globally optimal first waypoint sequence is output.

[0026] The route optimization output module uses B-spline interpolation to smooth the globally optimal first waypoint sequence, performs secondary verification based on the smoothing results, and outputs the verified optimal route.

[0027] Please see Figure 2As shown, the route optimization method based on particle swarm optimization and geographic information includes: S1: Collecting multi-source geographic information of the target voyage's operating area and preprocessing it to obtain a standard voyage operation dataset; S2: Constructing a two-layer voyage operation control constraint model based on the standard voyage operation dataset, including a hard operation constraint layer and a soft operation cost layer; S3: Pre-calculating the forward operation cost matrix and the backward operation cost matrix using a path search algorithm based on the voyage's starting point, ending point, and hard operation constraint layer, and fusing them to generate a bidirectional feasible region; S4: Initializing the particle swarm population based on the bidirectional feasible region, constructing a fitness function based on the soft operation cost layer, and iteratively optimizing and verifying it using a particle swarm optimization algorithm with an introduced operation path traction term, outputting the globally optimal first waypoint sequence; S5: Smoothing the globally optimal first waypoint sequence using B-spline interpolation, performing secondary verification based on the smoothing results, and outputting the verified optimal route.

[0028] S1: Collect multi-source geographic information of the target voyage's operating area and preprocess it to obtain a standard voyage operation dataset, including:

[0029] In this embodiment, it should be specifically explained that the target voyage operating area refers to the entire sea area / waterway where the ship will sail. It refers to the geographical space area surrounded by the starting point and ending point of the route that needs to be geographic information modeled and route optimization planned, including but not limited to airspace, sea area, and land transportation area.

[0030] S1.1: First, define the original multi-source geographic dataset D. raw It includes at least the following data sets: water depth topography data set, coastline data set, obstacle data set (such as reefs, shipwrecks, drilling platforms, etc.), ocean current data set (such as current velocity, current direction, etc.), wind field data set (such as wind speed, wind direction, etc.), wave data set (such as wave height, wave direction, period, etc.), sea ice data set (such as ice density, ice thickness, ice perimeter, etc.), and meteorological data set (such as air pressure, air temperature, visibility, etc.). Each data set D k ={(p i ,t i ,v i | i=1,2,...,N k}, N k D represents the total number of data points contained in the k-th data type. k Each element in the vector represents the spatial location vector of the i-th data point, the timestamp of the i-th data point, and the observation vector of the i-th data point; the spatial location vector p i =(x i ,y i ,z i ), (x i ,y i) represents the planar coordinates (such as latitude and longitude) in the original coordinate system, z i The depth (unit: meters, water depth is negative, elevation is positive), observation vector v i For example: ocean current data v i =(u i ,v i ), u i and v i These are the eastward and northward velocity components (unit: m / s), wind field data: v i =(ws i ,wd i ), ws i and wd i These represent wind speed (m / s) and wind direction (degrees), respectively; coordinate system one, converting the spatial position vector to p coordinates. i z p i z =T coord (p i ,src i ,p tar ), T coord () represents the coordinate transformation operator, p i src i and p tar These are the original spatial coordinates of the i-th data point before transformation, the original coordinate system parameters (original coordinate system parameters are standard terms in the surveying field, including reference ellipsoid parameters (such as semi-major axis, flattening), projection method (such as Mercator projection, UTM projection), central meridian, eastward offset, northward offset, etc. These parameters can usually be obtained from the data source metadata or by looking up a table according to the data source region. A coordinate transformation method known in the field, such as the Bursa-Wolf model or raster offset correction method, is used to transform the data point from the original coordinate system to the target coordinate system (such as WGS-84)). Time synchronization and interpolation are also considered. Let the target location be p, the target time be t, and a spatial neighborhood N(p) be defined with p as the center and rr (e.g., a radius of 10 kilometers) as the radius. A time window (e.g., ±1 hour) is defined with t as the center and Δt as the half-width. A set of known data points SS={(p) that simultaneously satisfies the spatial neighborhood and time window conditions is selected. i ,t i ,v i The estimated value v of the target point is calculated using a spatiotemporal weighted interpolation function. , In the data point set SS, obtain the variance σ of the data point location p. p 2 variance σ over time t 2 ,||ppi || 2 The square of the Euclidean distance between the positions, |tt i | represents the absolute time difference, and δ represents the minimum value (e.g., 10). -6 To avoid denoising to zero, filtering algorithms are used to denoise terrain, meteorological, and sensor data; outliers are determined using the 3σ criterion or quartile method, and the mean of adjacent valid values ​​is used to replace outliers; based on the target coordinate system, the target voyage operating area is divided into M×N regular grids g according to a preset resolution Δx×Δy (e.g., 0.01° for near-shore routes and 0.05° for ocean routes in the geodetic coordinate system). mn (Each raster has a unique row and column number (m, n), corresponding to a fixed geographic range), resulting in a set of regular rasters G. For each raster, all data points falling within that cell are collected, and data points of the same type are aggregated, for example, using an arithmetic mean or median filter, to obtain the representative attribute value of that raster for that type. By traversing all types of geographic information (water depth, ocean currents, wind fields, etc.), each raster has a corresponding attribute value, resulting in a geographic information dataset D. ba Let m and n be the raster row index and column index, respectively; let the coordinates of the lower left corner of the target coordinate system be (x0, y0) and the coordinates of the upper right corner be (x0, y0). max ,y max ), grid g mn The spatial range is [x0+(n-1)Δx,x0+nΔx]×[y0+(m-1)Δy,y0+mΔy], n=⌊(x-x0) / Δx⌋+1, m=⌊(y-y0) / Δy⌋+1, (x,y) are the coordinates in the target coordinate system, N=⌈(x max -x0) / Δx⌉+1,M=⌈(y max -y0) / Δy⌉+1;

[0031] S1.2: Based on geographic information dataset D ba Define a multi-source cost dataset C, where C = {C 1 C 2 ,...,C n1}, where C is the set of cost factors that need to be normalized, C n1 The cost factors to be considered for the n1th voyage, where n1 is the number of different cost factors, and each cost factor C... j In grid g mn The original value on is C mn j Then, by using the maximum and minimum values ​​of the j-th cost factor in all grids, C mn j Normalize to C mn,no j The normalized cost dataset C is obtained.no ={C no j |j=1,2,...,n1},C no j ={C mn,no j |m=1,...,M,n=1,...,N};Integrate geographic information dataset D ba and C no This yields the standard dataset for voyage operations.

[0032] This embodiment requires specific explanation that the cost factor includes at least the ocean current resistance cost C. cu Wind farm cost factor C wi Wave cost factor C wa Sea ice cost factor C ci And the cost factor C of fuel consumption fu In a preferred embodiment of the present invention, the plurality of basic cost factors and their calculation methods are as follows: Ocean current resistance cost factor C cu According to formula C cu =1-max(0,cosθ) is used for calculation, where θ is the angle between the heading and the current direction. When the heading and the current direction are in the same direction, θ=0°, cosθ=1, C cu =0, resistance cost is 0; when the heading is perpendicular / opposite to the flow direction, cosθ≤0, C cu =1, resulting in the highest drag cost; wind farm cost factor C wi Characterizing the effect of wind on ship navigation, headwind increases drag, while tailwind decreases drag. C wi =max(0,WS×cosθ1 / WS max ), WS, cosθ1 and WS max These are wind speed, the angle between wind speed and heading, and the maximum possible wind speed (e.g., 50 knots); wave cost factor C. wa Based on the significant wave height H s With the maximum possible wave height H s,max The square of the ratio (e.g., 15 meters) is used to determine the cost factor C of sea ice. ci Characterizing the impact of sea ice on navigation, higher ice density results in greater navigation risk and higher costs. ci =IC / 100×(1+IT / IT max IC represents ice density (ice coverage), with a value ranging from 0% to 100%. IT and IT max Ice thickness and maximum ice thickness, respectively; energy cost factor C fuBased on ship load and speed estimates, this is used to characterize the fuel consumption cost per unit distance of a ship under different operating conditions. Those skilled in the art will understand that ship energy consumption is related to various factors such as speed, load, hull fouling, and main engine efficiency. It can be estimated using well-known ship energy consumption models. As a specific implementation method, a simplified energy consumption model based on the naval coefficient method can be used: P E =Δ 2 / 3 ×V 3 / CC, C fu =P E ×SFOC×Pr fu / V,P E Where Δ is the ship's effective power, V is its displacement, CC is the naval coefficient, SFOC is the main engine fuel consumption rate, and Pr is the ship's effective power. fu For fuel prices; as another implementation method, a machine learning model trained on real ship data, such as random forest or gradient boosting tree, can be used to predict the fuel consumption rate at different speeds and loads based on historical navigation data, and the energy cost factor can be calculated by selecting any of the above methods or other energy consumption estimation methods known in the art.

[0033] S2: Based on the standard voyage operation dataset, construct a two-layer voyage operation control and constraint model, including a hard operation constraint layer and a soft operation cost layer, including:

[0034] S2.1: Based on geographic information dataset D ba Get the raster g mn Hard constraint value h mn If the grid g mn h satisfies any one type of hard constraint condition mn =0 indicates that the grid is not passable; otherwise, h mn =1 indicates that the grid is passable. A binary matrix hard operational constraint layer H is constructed, H=[h mn ] M×N h mn ∈{0,1}; The hard constraints include at least one or more of the following: land constraints (land areas, islands, peninsulas, etc.), shoal constraints (areas with water depth less than the safe navigation threshold for ships, which can be dynamically set according to the ship's draft), obstruction constraints (fixed obstacles such as reefs, shipwrecks, drilling platforms, and offshore wind power facilities), military restricted area constraints (military exercise areas, weapons testing areas, and permanently military-controlled waters announced by various countries), marine protected area constraints (prohibited navigation areas such as marine nature reserves, ecological red line areas, and coral reef protected areas), and special control area constraints (temporary restricted navigation areas such as high-risk pirate areas, epidemic isolation areas, and port quarantine anchorages).

[0035] This embodiment requires specific explanation of the hard constraints: First, multi-source geographic hard constraint data is collected, including global water depth data, no-navigation zones in official electronic charts (ENC), navigational obstruction layers, navigation warnings issued by maritime authorities of various countries, and land and shoal distribution interpreted from satellite remote sensing imagery. This multi-source geographic hard constraint data is then classified according to the constraint types required for the route and mapped to a unified raster coordinate system. For each constraint type c, a binary mask matrix H=[h mn ] M×N Land constraint data comes from the GSHHG coastline database, satellite imagery, etc.; shoal constraint data comes from GEBCO depth data, ENC charts, etc.; navigational obstruction constraint data comes from the ENC chart navigational obstruction layer; military restricted area data comes from national navigation notices, IMO reports, etc.; protected area data comes from announcements by national marine management departments, etc.; temporary control area data comes from real-time navigation warning services (International Maritime Organization (IMO), International Maritime Bureau (IMB) Piracy Reporting Center, national maritime safety notices, etc.). Constraint data from different sources and formats are mapped to a unified raster coordinate system: First, the coordinate reference of all data is unified to the target raster coordinate system (e.g., WGS-84); then, vector data is rasterized: for point-like navigational obstructions (e.g., reefs, shipwrecks, etc.), the row and column number of the target raster is calculated based on its coordinates (x, y): m = ⌊(yy) min ) / Δy⌋+1,n=⌊(xx min ) / Δx⌋+1, sets the mask of the corresponding grid to 1, (x min ,y min Let (x, y) be the minimum value in the x and y directions of the raster, and (Δx, Δy) be the resolution of the raster. For area-based no-navigation zones (such as military restricted zones and nature reserves), a polygon rasterization algorithm is used to determine whether the center point of each raster is located inside the polygon. If it is inside, the raster is marked as 1. This can be implemented using the ray method or the scanline algorithm. Secondly, raster data resampling: For geographic data in existing raster format (such as water depth networks), if its resolution or coordinate system is inconsistent with the target raster, resampling is performed. Let the value of the source raster at point (x, y) be v(x, y), then the value of the target raster g is... mn (The coordinates of its center point are (x) mn ,y mn Estimate the value v at that point. mnThe algorithm can take the value of the source raster point closest to the target point (dividing the source raster space into a Voronoi diagram (Thieson polygon) so that the target point falls within the resolution range of a certain source raster point), use linear interpolation (fitting the target point onto the line connecting two source raster points), or use bilinear interpolation (fitting the target point into the rectangle enclosed by four source raster points). Finally, for text data parsing and mapping: for constraint information in text formats such as navigation warnings, the described geographic coordinate range is first parsed. If it is a structured format (such as GeoJSON), the coordinates are directly extracted; if it is a natural language description, the latitude and longitude range is extracted through rule matching. The parsed coordinate range is treated as a temporary isometric feature and rasterized using vector data rasterization methods.

[0036] S2.2: Based on the normalized cost dataset C no Each cost element C mn,no j Weighted synthesis is performed to obtain the raster g. mn soft cost value s mn (Value range [0,1]) ω j The weight of the j-th cost factor (assigned directly based on the experience of the route planner, captain, and operations management personnel according to voyage requirements; for example, on a certain container liner route, the operations team determines: distance cost weight is 0.4, ocean current resistance weight is 0.2, weather risk weight is 0.1, and energy cost weight is 0.3; or the weights are obtained through the analytic hierarchy process (AHP) or entropy weight method), ω j ≥0, all weights sum to 1, construct a soft operating cost layer S, S=[s mn ] M×N Finally, the two-layer voyage operation control constraint model MM is obtained, where MM=(H,S,G), and G is the set of rule grids.

[0037] S3: Based on the voyage's origin, destination, and hard operational constraints layers, a path search algorithm is used to pre-calculate the forward and backward operational cost matrices, which are then fused to generate a bidirectional feasible region, including:

[0038] S3.1: Define the raster g in the hard operational constraint layer H mn Let the neighborhood set N(m,n) be defined, where N(m,n) = {(m',n')||m-m'|≤1,|n-n'|≤1,(m',n')≠(m,n)} (i.e., eight-neighbor connectivity, including top, bottom, left, right, and diagonal directions), and m' and n' be the row and column numbers of the candidate neighbor grids, respectively; define the grid from grid g mn Move to adjacent grid g ab The single-step operating cost function d(g) mn ,g ab If the hard constraint value h of adjacent grid cellsab =0, then d(g) mn ,g ab If h = ∞, ab =1, neighborhood soft cost fusion value s_avg=(s mn +s ab If it is a four-neighbor movement (|ma| + |nb| = 1), then d(g) / 2, mn ,g ab If d(g) = 1 × s_avg, and it is a diagonal neighborhood move (|ma| = 1 and |nb| = 1), then d(g) = 1 × s_avg. mn ,g ab )=2 1 / 2 ×s_avg;

[0039] S3.2: Employ the path search algorithm (Dijkstra's algorithm) to obtain the path from the rasterization starting point st of the flight path to the raster g. mn Minimum cumulative cost f mn f , π is the value from st to g mn A feasible path, f(st,g) mn ) for all from st to g mn The set of feasible paths, where |π| is the number of grid cells traversed by path π. t Let g be the t-th cell on path π, where the boundary condition is: if g mn =st, then f mn f =0, if the hard constraint value h mn =0, then f mn f =∞, construct the forward operating cost matrix F for =[f mn f ] M×N Similarly, define the distance from the rasterized endpoint EN of the flight path to the raster g. mn Minimum cumulative cost f mn b The boundary condition is: if g mn =EN, then f mn b =0, if h mn =0, then f mn b =∞, construct the backward operation cost matrix F bac =[f mn b ] M×N For each grid g mn Calculate the two-way cost V mn V mn =fmn f +f mn b V mn The smaller the value, the better the path the grid is on. Simultaneously, the maximum bidirectional cost V of all accessible grids is calculated. mn,max and minimum value V mn,min , , Finally, a bidirectional feasible region R is generated by setting a bidirectional cost threshold th, where R = [r]. mn ] M×N r mn ∈{0,1}, r mn For grid g mn The bidirectional feasible region label value, if h mn =1 and V mn ≤th, r mn =1 indicates that the grid g mn It belongs to a bidirectional feasible region, where particles can move freely; otherwise, r mn =0 indicates that the grid g mn Particles are prohibited from entering regions that are not within the bidirectional feasible region; the bidirectional cost threshold th can be an absolute threshold or a relative threshold, where the absolute threshold method is: th = V mn,min +α×(V mn,max -V mn,min ), where α is the threshold coefficient, α∈[0,1]. α=0 means only the grid cells on the shortest path are retained (most stringent). For open coastal routes, α=0.2; for complex island / reef routes, α=0.5; for ocean-going / transoceanic routes, α=0.1. α=1 means all passable grid cells are retained (most lenient). Relative thresholding method: th=β×V mn,min β is a multiplier coefficient, β≥1, β=1.2 indicates that the cost of retaining a grid does not exceed 1.2 times the minimum cost, and α and β can be obtained by training based on historical flight route data;

[0040] This embodiment specifically describes the use of Dijkstra's algorithm, well-known in the art, for searching under the constraints of a hard operational layer H. The difference between this search process and the standard Dijkstra algorithm lies in the following: the neighborhood of a node is defined as an eight-neighbor connected region, and the edge weight is set to 1 or 2 depending on the direction of movement. 1 / 2 And only if the hard constraint value g of the target node ab Expansion is only allowed when =1, or A can be used. * An algorithm can replace Dijkstra's algorithm by designing a suitable heuristic function (such as Euclidean distance) to guide the search direction. * When the heuristic function of the algorithm satisfies the acceptability (i.e., the heuristic value does not exceed the actual minimum cost from the current node to the destination), the resulting cost matrix is ​​equivalent to that of Dijkstra's algorithm.

[0041] Please see Figure 3 As shown, S4: Based on the bidirectional feasible region, initialize the particle swarm population, construct the fitness function based on the soft operating cost layer, and perform iterative optimization and verification using a particle swarm algorithm that incorporates an operating path traction term. Output the globally optimal first waypoint sequence, including:

[0042] S4.1: First, define the particle population size as K1 (e.g., 50-80). Each particle L contains a starting point, an ending point, and M1-2 intermediate waypoints, totaling M1 waypoints (M1 waypoints can be determined empirically based on the route type, such as M1=20 for ocean routes, M1=15 for coastal routes, and M1=10 for port entry / exit segments. These empirical values ​​can be adjusted based on the characteristics of the route and planning requirements. The M1-2 intermediate waypoints are generated within the bidirectional feasible region R). M1≥4, and the coordinates of each waypoint are p. L,m1 For the origin (st) and destination (EN), first convert their geographic coordinates to waypoint coordinates (p). L,1 and p L,M1 Then, based on the bidirectional feasible region R, the particle position X is initialized. L 0 X L 0 =[p L,1 0 ,p L,2 0 ,...,p L,M1-1 0 ,p L,M1 0 ], where intermediate waypoints satisfy p L,m1 0 ∈R, m1=2,3,...,M1-1, p L,m1 0 To initialize the position of the Lth particle at its m1-th waypoint; finally, initialize the particle velocity V. L 0 V L 0 =[0,v L,2 0 ,...,v L,M1-1 0 ,0],v L 0 =(v x,L 0 ,v y,L 0 Let v be the velocity component of the waypoint (in normalized position units / iteration step), where the velocity component satisfies: v x,L,m1 0 ∈[-v max ,v max ], vy,L,m1 0 ∈[-v max ,v max The starting and ending velocities are both 0, v max The maximum velocity threshold is typically set to 1% to 5% of the feasible region's side length.

[0043] S4.2: In the T-th iteration, the velocity V of the L-th particle at the m1-th waypoint. L,m1 T Combined with the operational path traction item G r (X L,m1 T Updated to V L,m1 T+1 V L,m1 T+1 =ω1×V L,m1 T +c1r1(P best,L,m1 -X L,m1 T )+c2r2(G best,m1 -X L,m1 T )+G r (X L,m1 T ), where ω1 is the inertia weight (controlling the tendency of the particle to maintain its original velocity; the inertia weight adopts a linear decreasing strategy, with an initial value set to 0.9, and iterates to the maximum number of iterations T). max When the time decays to 0.4, ω1(T) = 0.9 - [(0.9 - 0.4) / T] max [×T, which enables the algorithm to have a strong global search capability in the early stage to explore the feasible region), c1 and c2 are the individual optimal and global optimal learning factors (value range [1,2.5], such as c1=c2=1.5), r1 and r2 are random numbers in the interval [0,1], X L,m1 T Let P be the position of the L-th particle and the m1-th waypoint in the T-th iteration. best,L,m1 G represents the optimal position of the Lth particle at the m1th waypoint. best,m1 G represents the globally optimal position at the m1-th waypoint (the position with the best fitness found by all particles in all iterations). r (X L,m1 T ) is the operational path guiding item. X L,m1 T p is the m1th intermediate waypoint of the Lth particle in the Tth iteration. sp e For the shortest path P spThe target towing point coordinates are scaled to match the current waypoint position (grid coordinates are first converted to the target coordinate system, then normalized to the same coordinates as the waypoint position). The scaling rule is: e = ⌊m1 / (M1-1)×E⌋. Let P be the path with the minimum sum of forward and backward operating costs within the bidirectional feasible region. sp =[p sp 1 ,p sp 2 ,...,,p sp E E is the number of grid cells included in the shortest path, η is the traction coefficient, η∈[0.01,0.1], δ is the minimum value, G is the minimum value when m1=1 and m1=M1. r (X L,m1 T The value is 0 (no traction at the start / end point), and the start and end point positions remain unchanged; only intermediate waypoints participate in iterative optimization. The update velocity V of the Lth particle is obtained by traversing all waypoint positions and velocity updates. L T+1 The position update formula for the Lth particle is: X L T+1 =X L T +V L T+1 X L T Let c1r1 be the position vector of all waypoints of the Lth particle in the Tth iteration; where c1r1(P best,L,m1 -X L,m1 T ) and c2r2(G best,m1 -X L,m1 T () represents the change in position within a unit iteration step;

[0044] In this embodiment, it should be specifically noted that during the iterative optimization phase of the particle swarm optimization algorithm, to improve algorithm stability and convergence efficiency, the waypoint coordinates are the normalized coordinates of the waypoint positions. The coordinates (x, y) in the target coordinate system are converted to raster row and column numbers (m... d ,n d ), in the grid (m d ,n d Within the grid, the offset of point (x, y) relative to the bottom left corner (x0, y0) is (ξ_x, ξ_y), where ξ_x = [x - (x0 + (n...]]. d -1)Δx)] / Δx,ξ_y=[y-(y0+(m d -1)∆y)] / ∆y, then the waypoint position coordinates p=(x_no d ,y_no d )=((nd -1+ξ_x) / N,(m d -1+ξ_y) / M).

[0045] S4.3: For the new position X L T+1 For each waypoint in the equation, hard constraint verification and feasible region verification are performed. The hard constraint verification rule ISV is: if waypoint p L,m1 T+1 The grid g mn Hard constraint value h mn =1, then ISV is True, otherwise False. The feasible region check verifies ISI as follows: if waypoint p L,m1 T+1 The grid g mn The bidirectional feasible region label value r mn =1, then ISI is True; otherwise, it is False. If ISV is False or ISI is False, then waypoint position correction is triggered: p L,m1 T+1 =Correct(p L,m1 T ,v L,m1 T+1 The corrected particle position X is obtained by traversing all waypoints, where R and H are the bidirectional feasible region and hard operational constraint layers, respectively. The corrected particle position X is obtained by traversing all waypoints. L T+1,corr ;

[0046] This embodiment requires specific explanation regarding the situation where the particle, after being updated, flies out of the bidirectional feasible region R or the hard constraint value h of its grid cell. mn When p = 0, the position correction mechanism is triggered. As a specific implementation, the position correction function Correct can be corrected using the well-known backtracking + reflection method, for example, by pulling the particle back to the intersection of its historical trajectory and the boundary of the bidirectional feasible region, and adjusting its velocity direction, i.e., p L,m1 T+1 =p L,m1 T +γ×v L,m1 T+1 γ∈[0,1] is the backoff coefficient, and the maximum γ that satisfies feasibility is found through binary search or linear search: γ=max{γ∈[0,1]|ISV(p L,m1 T +γ×v L,m1 T+1 =True∩ISI(p L,m1 T +γ×v L,m1T+1 =True}; Simultaneously adjust the velocity direction: decompose the velocity into a normal component v_n and a tangential component v_t of the bidirectional feasible region boundary, reverse the normal component and multiply it by the bounce coefficient v L,m1 T+1 =-β1×v_n+v_t, where β1∈[0,1] is the bounce coefficient. The bounce coefficient β1 is determined based on the particle's out-of-bounds distance. When the out-of-bounds distance is ≤1 / 2 of the feasible region grid side length, β1=0.8; when the out-of-bounds distance is >1 / 2 of the feasible region grid side length, β1=0.5, so that the particle gradually converges into the feasible region after bounce. If the corrected particle position still does not meet the feasible region constraint, the backtracking method + bounce method is repeated until the waypoint falls into the feasible region. The maximum number of consecutive corrections is set (e.g., 3 times). If it is still not feasible, the waypoint is reinitialized.

[0047] S4.4: Connect any two adjacent waypoints (p L,m1 ,p L,m1+1 Defined as segment I (L,m1)_(L,m1+1) For segment I (L,m1)_(L,m1+1) From p L,m1 to p L,m1+1 The line segment rasterization algorithm (Bresenham line segment algorithm or DDA algorithm, both well-known techniques in the field of computer graphics) is used to obtain the line segment I. (L,m1)_(L,m1+1) The average soft cost s of all traversed graticules. mn I Euclidean distance ∆d between the flight segment and the flight segment I Calculate flight segment I (L,m1)_(L,m1+1) Total cost CC I , ,∆d I =||p L,m1 -p L,m1+1 ||2; Finally, the total fitness of the Lth particle is obtained. L , pr(X) L T Let pr(X) be the penalty term for the violation at the position of the Lth particle. If all waypoints of the particle satisfy the bidirectional feasible region constraint and the hard constraint, then pr(X) L T If pr(X) = 0, then if there is an infeasible violation point, then pr(X) = 0. L T )=λ1×CC avg ×N vi , where N is the penalty multiplier (recommended value is 2-5). vi The number of infeasible violations, CC avg This represents the average total cost across all segments of the particle's flight path.

[0048] S4.5: If the position of the Lth particle is XL T+1,corr Total fitness fit(X) L T+1,corr The optimal individual position P of the Lth particle best,L Total fitness fit(P) best,L ), then P best,L Updated to X L T+1,corr Otherwise, keep P best,L Unchanged; G best The globally optimal position G best The updated formula is: G best =argmin L=1,...,K1 fit(P best,L ), K1 is the particle population size; finally, when |fit(G best T+1 )-fit(G best T | < convergence threshold th1 (value range 10) -3 ~10 -5 (or the number of iterations T > the maximum number of iterations T) max If the iteration reaches 100 iterations, then stop iterating and output the globally optimal particle G. best This includes the coordinates of the grid center point and waypoint positions of the starting point, ending point, and the M1-2 optimal waypoints located in between;

[0049] S4.6: The globally optimal particle G best All waypoint position coordinates (x_no d ,y_no d Performing inverse normalization yields the raster row and column numbers as (m) d ,n d )(n d =min(N,⌊x_no d ×N⌋)+1,m d =min(M,⌊y_no d ×M⌋+1), then obtain the grid offset (ξ_x,ξ_y), and then obtain the waypoint coordinates (X,Y) in the target coordinate system (if the target coordinate system is the geodetic coordinate system, then X is longitude and Y is latitude), X=x0+(n d -1+ξ_x)×Δx,Y=y0+(m d -1+ξ_y)×Δy, (x0,y0) are the coordinates of the lower left corner of the target coordinate system, Δx and Δy are the grid resolution, and finally the globally optimal particle G after inverse normalization is obtained. best z This is the globally optimal first waypoint sequence;

[0050] S5: Smooth the globally optimal first waypoint sequence using B-spline interpolation, perform secondary verification based on the smoothing result, and output the verified optimal route, including:

[0051] S5.1: Using the globally optimal first waypoint sequence as the control points of a cubic B-spline curve, the number of curve segments is M1-3 (the number of control points is M1). Based on the control points, define the parametric equation C of the J-th curve segment. J (tt), controlled by point p J+1 p J+2 p J+3 and p J+4 Define t ∈ [0,1] as a local parameter, such as M1=5, the number of curve segments is 2, the 0th segment is defined by w1, w2, w3 and w4, the 1st segment is defined by w2, w3, w4 and w5, C J (tt) = 1 / 6[(1-tt)] 3 p J+1 +(3tt 3 -6tt 2 +4)p J+2 +(-3tt 3 +3tt 2 +3tt+1)p J+3 +tt 3 ×p J+4 According to the total length L of the route tot Given the desired sampling interval ∆s (e.g., 1.11 nautical miles), determine the total number of sampling points Q required for the route. tot Q tot =max(Q min ,⌈L tot / ∆s⌉), Q min Set the minimum number of sampling points (e.g., 100); calculate the arc length L of each curve segment. J (Arc length can be calculated using numerical integration or discrete myopia calculation), J=0,1,...,M1-4, based on the arc length accounting for L tot The total number of sampling points is allocated to each curve segment according to the ratio, resulting in the number of sampling points Z for the Jth curve segment. J Z J =max(2,⌊L J / L tot ×Q tot +0.5⌋) (at least 2 to ensure the endpoints are included); then for the J-th curve segment, take Z on the local parameter tt. J At each of the three equally spaced points, the local parameter value tt at each interval point is obtained. J,jj Calculate the corresponding sampling point q J,jj =C J (tt J,jj ), tt J,jj=jj / (Z J -1), jj=0,1,...,Z J -1; Finally, the sampling points of all curve segments are concatenated in sequence, and overlapping points between adjacent segments are removed (i.e., the last point of each segment is the same as the first point of the next segment), to obtain the complete smooth waypoint sequence G. best h , Q act Let Q be the number of sampling points for all curve segments. act ≠Q tot According to Q act With Q tot The absolute deviation |∆Q| is used to sort the curve segments in descending order of arc length. If Q act <Q tot Add one sampling point to the first |∆Q| longest curve segments in sequence. If Q act >Q tot Then, decrease the number of sampling points by 1 for the first |∆Q| longest curve segments, so that Q act =Q tot ;

[0052] S5.2: For each sampling point q jj Perform point feasibility verification, verifying the hard constraint verification rule ISV of the grid where it is located. If the verification result is satisfied, it is True; otherwise, it is False. For each pair of adjacent sampling points (q... jj ,q jj+1 Define segment II jj_jj+1 And perform segment feasibility verification, using a segment rasterization algorithm to obtain segment II. jj_jj+1 For all raster sequences G1 passed through, if any raster g mn ∈G1 and hard constraint value h mn If the value is 1, the feasibility verification result is True; otherwise, it is False. For infeasible points or flight segments with a verification result of False, the infeasible sampling points are projected to the nearest passable grid center to obtain the verification-corrected sampling point q. jj co q jj co =argmin cc∈C_va ||q jj -cc||2, C_va is the set of center points cc of all passable grids. After completing the single-point projection correction, the feasibility verification of the segment is re-executed for all segments composed of adjacent sampling points. If there are still infeasible segments, then segment II is... jj_jj+1Waypoint encryption is performed by adding sampling points at both ends of infeasible areas and projecting them to the nearest passable area until all segments meet the hard constraints, ensuring the continuity and compliance of the smooth route throughout. Finally, a waypoint sequence G with secondary verification is obtained. best hy This is the verified optimal route;

[0053] Secondly: The accompanying drawings of the embodiments disclosed in this invention only involve the structures involved in the embodiments disclosed in this invention. Other structures can refer to the general design. In the absence of conflict, the same embodiment and different embodiments of this invention can be combined with each other.

[0054] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A route optimization method based on particle swarm optimization and geographic information, characterized by: include: S1: Collect multi-source geographic information of the target voyage's operating area and preprocess it to obtain the voyage's standard operating dataset; S2: Based on the standard dataset of voyage operations, construct a two-layer voyage operation control and constraint model, including a hard operation constraint layer and a soft operation cost layer; The hard operational constraint layer: based on geographic information dataset D ba Get the raster g mn Hard constraint value h mn If the grid g mn h satisfies any one type of hard constraint condition mn =0 indicates that the grid is not passable; otherwise, h mn =1 indicates that the grid is passable. A binary matrix hard operational constraint layer H is constructed, H=[h mn ] M×N h mn ∈{0,1}; The soft operating cost layer: based on the normalized cost dataset C no Each cost element C mn,no j Weighted synthesis is performed to obtain the raster g. mn soft cost value s mn Construct a soft operating cost layer S, S=[s mn ] M×N Finally, the two-layer voyage operation control constraint model MM is obtained, where MM=(H,S,G), and G is the set of rule grids. S3: Based on the voyage start point, end point and hard operational constraints layer, the forward operational cost matrix and the backward operational cost matrix are pre-calculated using a path search algorithm and then fused to generate a bidirectional feasible region; The bidirectional feasible region includes: S3.1: defining the grid g in the hard operational constraint layer H. mn The set of neighborhoods N(m,n); defined from the grid g mn Move to adjacent grid g ab The single-step operating cost function d(g) mn ,g ab ); S3.2: Employ a path search algorithm to obtain the path from the starting point st of the flight path to the grid g. mn Minimum cumulative cost f mn f The boundary condition is: if g mn =st, then f mn f =0, if the hard constraint value h mn =0, then f mn f =∞, construct the forward operating cost matrix F for =[f mn f ] M×N Similarly, define the distance from the route endpoint EN to the grid g. mn Minimum cumulative cost f mn b The boundary condition is: if g mn =EN, then f mn b =0, if h mn =0, then f mn b =∞, construct the backward operating cost matrix F bac =[f mn b ] M×N For each grid g mn Calculate the two-way cost V mn V mn =f mn f +f mn b Simultaneously calculate the maximum bidirectional cost V of all accessible grid cells. mn,max and minimum value V mn,min Finally, a bidirectional feasible region R is generated by setting a bidirectional cost threshold th. S4: Based on the bidirectional feasible region, initialize the particle swarm population, construct the fitness function based on the soft operating cost layer, and perform iterative optimization and verification by introducing the operating path traction term of the particle swarm algorithm, and output the globally optimal first waypoint sequence; The particle swarm optimization algorithm that introduces the operational path traction term: In the T-th iteration, the velocity V of the L-th particle at the m1-th waypoint... L T Combined with the operational path traction item G r (X L,m1 T Updated to V L,m1 T+1 ; By iterating through all waypoints and velocity updates of the Lth particle, the update velocity V of the Lth particle is obtained. L T+1 The position update formula for the Lth particle is: X L T+1 =X L T +V L T+1 X L T Let L be the position vector of all waypoints of the Lth particle in the Tth iteration; For the new location X L T+1 For each waypoint in the equation, hard constraint verification and feasible region verification are performed. The hard constraint verification rule ISV is: if waypoint w L,m1 T+1 The grid g mn Hard constraint value h mn =1, then ISV is True, otherwise it is False. The feasible region check verifies ISI as follows: if waypoint w L,m1 T+1 The grid g mn The bidirectional feasible region label value r mn =1, then ISI is True; otherwise, it is False. If ISV is False or ISI is False, then waypoint position correction is triggered: w L,m1 T+1 =Correct(w L,m1 T ,v L,m1 T+1 The function `Correct()` is the position correction function, where `R` and `H` represent the bidirectional feasible region and the hard operational constraint layer, respectively. Finally, the corrected particle position `X` is obtained by traversing all waypoints. L T+1,corr ; Any two adjacent waypoints (p) L,m1 ,p L,m1+1 ) is defined as segment I, and for segment I, from p L,m1 to p L,m1+1 A line segment rasterization algorithm is used to obtain all the grids traversed by the flight segment, and the average soft cost s of all grids is combined. mn I Euclidean distance △d between the flight segment and the flight segment I Calculate the total cost CC of flight segment I. I Finally, the total fitness of the Lth particle is obtained. L ; S5: The globally optimal first waypoint sequence is smoothed using B-spline interpolation. The smoothing result is then used for secondary verification, and the verified optimal route is output.

2. The route optimization method based on particle swarm optimization and geographic information according to claim 1, characterized in that: The implementation of S4 includes: S4.1: First, define the particle population size as K1. Each particle L contains a total of M1 waypoints, including a starting point, an ending point, and M1-2 intermediate points. Each waypoint m1 includes the coordinates w of the grid center point. L,m1 and waypoint position coordinates p L,m1 Then, based on the bidirectional feasible region R, the particle position X is initialized. L 0 Finally, initialize the particle velocity V. L 0 .

3. The route optimization method based on particle swarm optimization and geographic information according to claim 1, characterized in that: The implementation of S4 also includes: S4.5: If the position of the Lth particle is X L T+1,corr Total fitness fit(X) L T+1,corr The optimal individual position P of the Lth particle best,L Total fitness fit(P) best,L ), then P best,L Updated to X L T+1,corr Otherwise, keep P best,L Unchanged; G best The globally optimal position G best The updated formula is: G best =argmin L=1,...,K1 fit(P best,L ), K1 is the particle population size; finally, when |fit(G best T+1 )-fit(G best T | <convergence threshold th1 or iteration count T> maximum iteration count T max If the iteration stops, the globally optimal particle G is output. best This includes the coordinates of the grid center point and waypoint positions of the starting point, ending point, and the M1-2 optimal waypoints located in between; S4.6: The globally optimal particle G best All waypoint position coordinates (x_no d ,y_no d Performing inverse normalization yields the raster row and column numbers as (m) d ,n d Then, the offset within the grid (ξ_x, ξ_y) is obtained, followed by the waypoint coordinates (X, Y) in the target coordinate system. Finally, the globally optimal particle G after inverse normalization is obtained. best z This is the globally optimal first waypoint sequence.

4. The route optimization method based on particle swarm optimization and geographic information according to claim 1, characterized in that: The implementation of S5 includes: S5.1: Using the globally optimal first waypoint sequence as the control points of a cubic B-spline curve, the number of curve segments is M1-3, and the parametric equation C of the J-th curve segment is defined based on the control points. J (tt); based on the total length of the route L tot Given the desired sampling interval Δs, determine the total number of sampling points Q required for the flight path. tot ; Calculate the arc length L of each curve segment. J J=0,1,...,M1-4, according to the arc length as a percentage of L tot The total number of sampling points is allocated to each curve segment according to the ratio, resulting in the number of sampling points Z for the Jth curve segment. J Next, for the J-th curve segment, take Z on the local parameter tt. J At each of the three equally spaced points, the local parameter value tt at each interval point is obtained. J,jj Calculate the corresponding sampling point q J,jj Finally, the sampling points of all curve segments are concatenated sequentially, and overlapping points between adjacent segments are removed to obtain the complete smooth waypoint sequence G. best h ; S5.2: For each sampling point q jj Perform point feasibility verification, verifying the hard constraint verification rule ISV of the grid where it is located. If the verification result is satisfied, it is True; otherwise, it is False. For each team's adjacent sampling points (q jj ,q jj+1 To verify the feasibility of the segment, a line segment rasterization algorithm is used to obtain all raster sequences G1 traversed by the segment. If any raster g mn ∈G1 and hard constraint value h mn If the value is 1, the feasibility verification result is True; otherwise, it is False. For infeasible points or flight segments with a verification result of False, the infeasible sampling points are projected to the nearest passable grid center to obtain the verification-corrected sampling point q. jj co After completing the single-point projection correction, the feasibility verification of the segments composed of all adjacent sampling points is re-executed. If there are still infeasible segments, waypoints are densified for those segments, and new sampling points are added at both ends of the infeasible areas and projected to the nearest passable area. This process continues until all segments meet the hard constraints, and finally, the waypoint sequence G of the secondary verification is obtained. best hy This is the optimal route after verification.

Images