Road planning method between upper reservoir and lower reservoir of pumped storage power station based on improved A* algorithm

Abstract

The application discloses a kind of based on the improved A* algorithm upper and lower reservoir road planning method of pumped storage power station, including S1, the three-dimensional terrain point cloud data of pumped storage power station target area is rasterized, and digital elevation model DEM is generated;S2, identify natural obstacle area in target area, generate impassable area matrix map;S3, construct heuristic function considering distance and slope, enhanced cost function, and multilayer neighborhood search strategy, path search is carried out based on improved A* algorithm, and initial path is constructed;S4, the initial path is optimized, and the final path is obtained;S5, the final path after optimization is superimposed to original terrain grid map and visualized, and the key design index of path is calculated and displayed, to assist designer to carry out scheme comparison and engineering decision-making.The method is high in intelligent degree, and strong in engineering adaptability, and is suitable for the rapid route selection and optimization design of upper and lower reservoir road of complex mountain pumped storage power station.

Description

Technical Field

[0001] This invention belongs to the field of computer-aided design and terrain path planning technology, specifically relating to a method for planning roads between the upper and lower reservoirs of a pumped storage power station based on an improved A* algorithm. Background Technology

[0002] Pumped storage power stations are usually located in mountainous areas with complex terrain. The design of the road between the upper and lower reservoirs faces many challenges, such as large topographic relief, steep slopes, and complex geological conditions.

[0003] Current mountain road design still relies heavily on engineers' experience and traditional CAD software, resulting in low efficiency and accuracy dependent on subjective judgment. Some studies have attempted to use automated path search algorithms, such as A* and Dijkstra's algorithm, for path planning. However, these algorithms often prioritize shortest distance as their optimization objective, neglecting the complex constraints of mountain road design, such as maximum / average slope, turning radius, and geological safety. Consequently, the generated paths are difficult to apply directly in practical engineering projects.

[0004] In recent years, the development of UAV technology has made it possible to acquire high-precision terrain point cloud data, providing rich data support for road planning. The introduction of point cloud data provides a foundation for automatic route selection, but due to the dramatic changes in mountainous terrain, uneven distribution of point density, and data gaps in some areas, direct use can easily lead to path distortion. Furthermore, the traditional A* algorithm suffers from local optima and path non-smoothness in grid-level path search, further limiting its engineering applicability. Summary of the Invention

[0005] The purpose of this invention is to provide a road planning method between the upper and lower reservoirs of a pumped storage power station based on an improved A* algorithm, which solves the technical problems of poor adaptability of road selection to complex mountainous terrain, weak integration of design specifications, and lack of feasibility of path results in the existing technology.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following solution:

[0007] A road planning method between the upper and lower reservoirs of a pumped storage power station based on an improved A* algorithm includes the following steps:

[0008] Step S1: Rasterize the 3D topographic point cloud data of the target area of ​​the pumped storage power station to generate a digital elevation model (DEM). By converting complex terrain into structured data (raster matrix) that can be processed by algorithms, it facilitates subsequent digital analysis; the DEM provides accurate elevation information, laying the foundation for slope calculation, obstacle identification, etc., and ensuring the quantitative accuracy of terrain constraints.

[0009] Step S2: Based on the elevation model (DEM), identify natural obstacle areas corresponding to rivers, valleys, and steep slopes in the target area, generating an inaccessible area matrix map to provide obstacle avoidance information for path search. By identifying natural obstacles, paths avoid traversing high-risk or restricted areas, reducing engineering risks and compliance hazards; the matrix representation allows the algorithm to quickly identify "restricted areas," reducing invalid searches and improving planning efficiency.

[0010] Step S3: Construct a heuristic function that comprehensively considers distance and slope, an enhanced cost function, and a multi-layer neighborhood search strategy. Based on the improved A* algorithm, perform path search to construct an initial path. By designing a heuristic function that integrates distance (shortening mileage) and slope (reducing travel difficulty), the "steep path" problem caused by the traditional A* algorithm focusing only on the shortest distance is avoided, improving path safety and engineering feasibility. By designing a multi-layer neighborhood search strategy (such as prioritizing the expansion of low-cost areas), computational redundancy is reduced, enabling the rapid generation of an initial path that balances efficiency and quality in complex terrain.

[0011] Step S4: Optimize the initial path. First, remove redundant turning nodes. Then, use a third-order Bézier curve to interpolate and smooth the path, resulting in the final path. Removing redundant nodes makes the path simpler, reduces unnecessary turns, and lowers construction difficulty and earthwork volume. Furthermore, using Bézier curves for smoothing avoids abrupt path changes, reduces vehicle bumps, and improves operational safety and comfort.

[0012] Step S5: Path Result Visualization and Key Parameter Extraction: The optimized final path is overlaid onto the original terrain raster map for visualization, and key design indicators of the path are calculated to assist designers in comparing different options and making engineering decisions. Visualizing the planned path overlaid with the terrain helps designers quickly understand the spatial layout of the scheme, reducing communication costs. By quantifying and displaying indicators such as length, slope, and earthwork volume, subjective judgment biases are avoided, supporting scientific comparison of multiple schemes and accelerating decision-making efficiency.

[0013] Further optimization, step S1 specifically includes:

[0014] Step S1.1, Point Cloud Boundary Range Extraction and Raster Parameter Definition: Extract the measured 3D point cloud data P = {(x...} of the target area of ​​the pumped storage power station. i ,y i ,z i )|i=1,2,…,Q}, using a fixed grid size G, the point cloud data of the entire terrain area is spatially discretized, and the total number of grid rows and columns is calculated based on the point cloud boundary according to formula (1).

[0015]

[0016] In the above formula, where (xi ,y i (x) represents planar coordinates, z i n represents the elevation value, Q represents the total number of point clouds; x n is the total number of columns in the grid elevation matrix. y The grid represents the total number of rows in the grid elevation matrix; G represents the side length of the grid cell, x max This represents the maximum value of the point cloud data of the target region at the boundary in the x-direction, x min This represents the minimum value of the point cloud data of the target region at the boundary in the x-direction, and y-direction... max This represents the maximum value of the point cloud data of the target region at the boundary in the y-direction. min This represents the minimum value of the point cloud data of the target region at the boundary in the y direction. This indicates the rounding up operation.

[0017] Step S1.2: Map the point cloud to the raster index: According to formula (2), for each point (x i ,y i ,z i ) Calculate its raster index in the x and y directions respectively, and map each point to the corresponding raster index (I). i J i That is, the point cloud is assigned to the corresponding raster cells, and each raster cell (I, J) is associated with a set of points;

[0018]

[0019] In the above formula, x i Let y be the X-coordinate of the i-th point. i Let I be the Y-coordinate of the i-th point; k I represents the column index of the grid cell. i ∈[0, n x ); J i J represents the row index of the grid cell. i ∈[0, n y ), This indicates the floor function.

[0020] Step S1.3: Generate the elevation matrix model: Calculate the elevation of the points falling in each grid cell (I, J) to form the corresponding elevation set Z. I,J ={z i |(I i =I,J i =J)}.

[0021] If a grid cell contains at least one point, then the representative elevation H(I, J) of that grid cell is:

[0022]

[0023] In the above formula, Z I,J H(I,J) represents the set of elevation values ​​for all points within a raster cell; H(I,J) represents the representative elevation of the raster cell; |Z I,J | represents the number of elements in the elevation set.

[0024] If there are no points in the grid cell, its elevation is marked as missing, and a representative elevation is repaired by interpolation.

[0025] Finally, the representative elevations of all grid elements constitute a two-dimensional elevation matrix (DEM), denoted as...

[0026] Further optimization involves step S1.3, where for raster cells lacking elevation values, a Gaussian weighted quadratic polynomial is used to fit their representative elevation. Specifically, this includes:

[0027] Step S1.3.1: For each raster cell (I, J) with missing elevation values, construct a neighborhood window with its geometric center (x0, y0) as the target point, and search for points (x0, y0) in the neighboring valid raster cells. j ,y j ,z j Let n points be found in total, forming a neighborhood point set. n is a positive integer greater than or equal to 6; a valid grid cell refers to a grid cell containing at least one point.

[0028] Step S1.3.2, Solving for the coefficients A of the quadratic polynomial, specifically:

[0029] Let the quadratic polynomial be of the form F(x,y)=a0+a1x+a2y+a3x 2 +a4xy+a5y 2 ;

[0030] Based on the neighborhood point set formed in step S1.3.1 Construct the design matrix X, weight matrix W, and neighborhood point set elevation vector Z, and then apply the normal equation (X... T WX))A=X T WZ solves for the coefficient A.

[0031] Where Z = [z1, z2, ..., z n ] T A = [a0, a1, a2, a3, a4, a5] T ,

[0032]

[0033]

[0034] In the above formula, R is the smoothing coefficient.

[0035] Step S1.3.3: Substitute the coefficient A obtained in step S1.3.2 into the quadratic polynomial, and then substitute the target point (x0, y0) into the fitting function to calculate the interpolated elevation value z0 of the grid cell.

[0036]

[0037] Further optimization involves using the following methods to identify natural obstacles, specifically including:

[0038] Step S2.1, Slope Identification: Calculate the slope s of each grid cell (I,J) according to formula (4). I,J The elevation values ​​of this grid cell and its adjacent grid cells are used to calculate the following:

[0039]

[0040] In the above formula, d I,J→M,N This represents the horizontal distance between the current raster cell (I,J) and its neighboring raster cells (M,N), where M∈[0,n]. x ), N∈[0, n y );

[0041] For each grid cell that is adjacent to its four surrounding grid cells, i.e., a four-neighborhood, let d. I,J→M,N =G; For each grid cell that is adjacent to its eight surrounding grid cells, i.e., the eight-neighborhood, take...

[0042] If the slope is s I,J ≥s th If so, the grid unit is marked as impassable, i.e., an obstacle area;

[0043] If the slope is s I,J <s th If so, then mark the grid unit as passable;

[0044] s th For the set slope threshold, such as s th =0.3, meaning the slope is ≥30%.

[0045] Step S2.2, River Identification: For a certain continuous region Ω, if the standard deviation of its internal elevation satisfies formula (5), it is determined to be a candidate area for water bodies;

[0046]

[0047] Then, the connected components of the raster graph are used to analyze the candidate water bodies. Areas that extend in a strip shape and connect the upper and lower reservoirs of the pumped storage power station are identified as river channels, i.e., obstacle areas, used to exclude non-area water bodies such as lakes. Among them, continuous flat areas are candidate water bodies, and strip-shaped extensions are linear characteristics of river channels.

[0048] In the above formula, σ th The set elevation standard deviation threshold, such as σ th =0.5m; Ω represents a continuous region, referring to a set of spatially adjacent raster cells in a DEM, used for local terrain analysis, such as detecting flat water bodies; |Ω| represents the number of raster cells in the continuous region; This represents the average elevation value of a continuous area.

[0049] Step S2.3: Map the identified obstacle areas onto a two-dimensional grid map to generate an obstacle matrix map;

[0050] If a grid cell (I,J) is located in a passable area, then its corresponding map(I,J) = 1;

[0051] If a grid cell is located in a non-passable area, then its corresponding map(I,J) = 0.

[0052] Further optimization involves incorporating restricted areas as additional barrier areas into the map. These restricted areas include, but are not limited to, construction control zones and protection zones.

[0053] Further optimization involves step S3, which performs path search based on an improved A* algorithm, specifically including...

[0054] Step S3.1: Construct the mapping between grid cells and nodes: Define a node m = (I, J, g(m), h(m), f(m), parent), where each grid cell corresponds to a unique node. The node inherits the slope and passability attribute information of the corresponding grid cell. Here, (I, J) represents the row and column index of the grid cell; f(m) represents the cost function value, g(m) is the actual cost, and h(m) is the heuristic function; parent represents the parent node pointer, which is used for path backtracking.

[0055] Step S3.2: Construct a heuristic function h(m) that comprehensively considers the adaptability of distance and slope:

[0056] h(m)=α·d norm (m)+β·Δs norm (m) (6)

[0057] In the above formula, This represents the normalized Euclidean distance from the current node m to the target node, which is the Euclidean distance d(m,goal) from the current node m to the target node divided by the maximum Euclidean distance D from the starting point to the target node. max ; This indicates the slope deviates from the normalized value; Δs(m) = |s(m) - s opt | represents the node slope s (m) and the recommended design slope s opt The degree of deviation, s max The maximum allowable slope is represented by s(m); s(m) represents the slope of the current node, which is obtained by equation (7), that is, the elevation difference between the current node m and its parent node p divided by the horizontal distance.

[0058]

[0059] α and β represent the distance weight and slope weight, respectively, and their sum is 1. These values ​​are adjustable according to engineering requirements; the recommended initial values ​​are α = 0.6 and β = 0.4. Considering actual path variations and terrain complexity, the weights α and β are set to be dynamically adjustable, as follows:

[0060] α = 1 - λ·ρ(m)

[0061] β=λ·ρ(m)

[0062] In the formula, λ∈[0,1] is the weight adjustment coefficient; ρ(m) represents the path local complexity index, which is taken as the average slope change rate in this invention. When the terrain is relatively flat, distance is emphasized; when the slope changes drastically, the focus is increased on the slope.

[0063] Step S3.3: Introduce the slope penalty cost function and construct the enhanced cost function g(m), whose expression is shown in equation (8):

[0064] g(m)=g(p)+d(m,p)+γ·q(m) (8)

[0065] In the above formula, g(m) represents the total cost of the parent node; d(m,p) represents the actual distance between the current node m and the parent node; q(m) represents the slope penalty function, defined as:

[0066]

[0067] In the formula, γ is the penalty factor adjustment parameter; s opt The desired slope value is usually taken as 5%; s max To ensure maximum permissible slope, pumped storage power station projects are typically limited to 9%.

[0068] Step S3.4: For complex path planning regions, a multi-level neighborhood search strategy is adopted, specifically: Define the multi-level neighborhood node set N of the current node m. r(I,J):

[0069] N r (I,J)={(M,N)||MI|≤r,|NJ|≤r}(10)

[0070] In the above formula, r represents the neighborhood layer number, which is set to 1 to 10 to flexibly adjust the search range and reduce the search depth without sacrificing accuracy. (I,J) is the raster cell index corresponding to the current node m.

[0071] Step S3.5: Based on steps S3.1-S3.4, perform path search, specifically including:

[0072] Step S3.5.1: Determine the starting node and the target node, create an open list to store nodes to be explored, create a closed list to store nodes that have been explored, and set a parent node pointer for each node to track the path.

[0073] Step S3.5.2: Add the starting point to the open list and set its total cost to 0. Then, repeat the following steps until the target is found or the open list is empty:

[0074] Step S3.5.2.1: Select the node m with the smallest total cost f(m) in the open list as the current node, f(m) = g(m) + h(m), and determine whether the current node is the target:

[0075] If so, end the search;

[0076] Otherwise, remove the current node m from the open list and add it to the closed list;

[0077] Step S3.5.2.2: Traverse the multi-level neighborhood set N of the current node m. r Collect all child nodes in (I,J) and determine whether each child node satisfies the preset constraints; if not, delete the child node; if satisfied, continue to the next step.

[0078] Step S3.5.2.3: Calculate the total generation value of the child nodes that meet the conditions, and store these child nodes in the open list in ascending order. Update the open list, and then start the judgment again from step S3.5.2.1.

[0079] Step S3.5.3: Starting from the target, backtrack along the parent node pointer to the starting point to construct a complete initial path.

[0080] Further optimization involves employing a bidirectional A* search algorithm in step S3.5, specifically by setting a starting point P0 and a target P0. t Simultaneously from the starting point P0 and the target P tA path search is performed towards the center using the improved A* search algorithm until the two search paths meet. The forward and reverse paths are then obtained, and the meeting node is denoted as p. z The forward and reverse paths are merged to obtain the initial path {P0,…,P}. z ,…,P t}, ensure the passability of the corresponding grid unit Pz, i.e., map[I z J z ] = 1.

[0081] Further optimization, in step S4, involves optimizing the initial path, specifically including:

[0082] Step S4.1: For the generated initial path node sequence {P0, P2, ..., P...} t The sliding window method is applied to optimize redundant nodes; specifically, the initial path points {P0, P2, ..., P} are traversed. t For each path node P k Calculate vector and The included angle θ k , k∈(0,t).

[0083] If θ k If ≤θ0, then delete the path node P. k At this point, the node is a redundant turning point, and it needs to be extracted to reduce construction complexity, such as the number of fill sections. If θ k If θ > 0, then retain the path node P. k As a curve control node, it ensures that the path curvature conforms to the specifications, such as the minimum half.

[0084] Arrange the starting node, the target node, and all control nodes in sequence to construct the critical path point set {P0′, P1′, ..., P′}. t′}

[0085] Step S4.2: Interpolate the critical path points selected in step S4.1 using a third-order Bézier curve to eliminate the inflection points of the broken line, thereby generating a continuous and smooth path curve; specifically, starting from the first point P0′, sequentially select P′... l ,P′ l+1 ,P′ l+2 ,P′ l+3 As a group, construct the third-order Bézier curve as shown in equation (10) to obtain the corresponding interpolation nodes. l∈[0,3,6,…3n'], where n' is a positive integer.

[0086]

[0087] In the above formula, For parameters, Control the interpolation progress from the starting point to the target;

[0088] The starting node, target node, all control nodes, and corresponding interpolation nodes are arranged in sequence to form a secondary optimization path.

[0089] Step S4.3: Based on the obstacle matrix generated in step S2, check each node of the secondary optimized path to determine whether it crosses an impassable area. If there are non-compliant points, perform local backtracking search or forced avoidance processing, and finally obtain the final path.

[0090] Further optimization, step S5 specifically includes:

[0091] Step S5.1: Based on the final path node set {P0″, P″1, ..., P″} t} Calculate the key design parameters using the following formulas:

[0092] 1) Total path length L total :

[0093]

[0094] 2) Average slope S avg :

[0095]

[0096] 3) Maximum slope S max :

[0097]

[0098] In the above formula, u∈[0, t], node P u The three-dimensional coordinates of "″ are (x u ,y u ,z u );

[0099] Step S5.2: Compare the above key parameters with the set road design specifications to determine whether they meet the design specifications and mark them accordingly;

[0100] Step S5.3: When optimizing parameters to obtain multiple alternative routes, score and rank them according to key parameters to assist designers in making the best decision. Specifically:

[0101] Score = χ1L total +χ2S avg +χ3S max ; χ1+χ2+χ3=1, where χ1, χ2, and χ3 represent the weights corresponding to the key parameters respectively.

[0102] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0103] 1. The method described in this invention can efficiently generate upper and lower library connection paths that conform to engineering specifications under complex mountainous terrain conditions. By using 3D point cloud data to generate a digital elevation model (DEM) through rasterization, it can realistically reproduce the characteristics of complex mountainous terrain, adapting to problems such as changes in point density and abrupt changes in terrain elevation, providing accurate data support for road alignment. Through spatial analysis methods such as slope analysis, river identification, and low-lying area extraction, an inaccessible area matrix is ​​constructed, effectively avoiding unfavorable terrain obstacles and improving the safety and feasibility of road alignment.

[0104] 2. The method described in this invention introduces a joint cost function of distance and slope, a multi-layer neighborhood search mechanism, and a slope penalty weight control strategy on the basis of the traditional A* algorithm, which significantly improves the efficiency and rationality of path planning, and is particularly suitable for the design of mountain roads where the average slope and turning radius need to be controlled.

[0105] 3. This invention realizes intelligent processing of the entire process from terrain data processing and obstacle analysis to intelligent path generation and result analysis, which significantly improves the efficiency and quality of road design for pumped storage power stations. It is especially suitable for large-scale energy infrastructure projects with complex mountainous terrain and tight design cycles. Attached Figure Description

[0106] Figure 1 This is an overall flowchart of the road planning method between the upper and lower reservoirs of a pumped storage power station based on the improved A* algorithm described in this invention.

[0107] Figure 2 This is a schematic diagram of the three-layer neighborhood grid division used in this embodiment of the invention;

[0108] Figure 3 This is a flowchart of the improved bidirectional A* algorithm search in an embodiment of the present invention;

[0109] Figure 4 This is a schematic diagram of path planning control points and search routes for a local grid area in an embodiment of the present invention;

[0110] Figure 5 This is a path comparison diagram between the automatic route selection result and the manual route selection result in an embodiment of the present invention. Detailed Implementation

[0111] The present invention will now be described in detail with reference to specific embodiments to enable those skilled in the art to better understand and implement the invention. However, it should be understood that the scope of protection of the present invention is not limited to the specific embodiments listed. For those skilled in the art, any equivalent transformations, substitutions, adjustments, or extensions made without departing from the essential spirit and scope defined by the claims should be considered to fall within the scope of protection of the present invention.

[0112] like Figure 1 As shown, a smart planning method for the upper and lower reservoir roads of a pumped storage power station based on point cloud and improved A* algorithm is presented:

[0113] Step S1: Rasterize the 3D topographic point cloud data of the target area of ​​the pumped storage power station to generate a digital elevation model (DEM), specifically including:

[0114] Step S1.1, Point Cloud Boundary Range Extraction and Raster Parameter Definition: Extract the measured 3D point cloud data P = {(x...} of the target area of ​​the pumped storage power station. i ,y i ,z i )|i=1,2,…,Q}, using a fixed grid size G, the point cloud data of the entire terrain area is spatially discretized, and the total number of grid rows and columns is calculated based on the point cloud boundary according to formula (1);

[0115]

[0116] In the above formula, where (x i ,y i (x) represents planar coordinates, z i n represents the elevation value, Q represents the total number of point clouds; x n is the total number of columns in the grid elevation matrix. y The grid represents the total number of rows in the grid elevation matrix; G represents the side length of the grid cell, x max This represents the maximum value of the point cloud data of the target region at the boundary in the x-direction, x min This represents the minimum value of the point cloud data of the target region at the boundary in the x-direction, and y-direction... max This represents the maximum value of the point cloud data of the target region at the boundary in the y-direction. min This represents the minimum value of the point cloud data of the target region at the boundary in the y direction. This indicates the rounding up operation.

[0117] Step S1.2: Map the point cloud to the raster index: According to formula (2), calculate the x-axis of each point. i ,y i ,z i Map each point to its corresponding raster index (I) in both the x and y directions. i Ji That is, the point cloud is assigned to the corresponding raster cells, and each raster cell (I, J) is associated with a set of points;

[0118]

[0119] In the above formula, x i Let y be the X-coordinate of the i-th point. i Let I be the Y-coordinate of the i-th point; k I represents the column index of the grid cell. i ∈[0, n x ); J i J represents the row index of the grid cell. i ∈[0, n y ), This indicates the floor function.

[0120] Step S1.3: Generate Elevation Matrix Model (EDM): Calculate the elevation of the points falling within each grid cell (I, J) to form the corresponding elevation set Z. I,J ={z i |(I i =I,J i =J)};If there is at least one point in the grid cell, then the representative elevation H(I,J) of the grid cell is as shown in equation (3):

[0121]

[0122] In the above formula, Z I,J H(I,J) represents the set of elevation values ​​for all points within a raster cell; H(I,J) represents the representative elevation of the raster cell; |ZI,J| represents the number of elements in the elevation set.

[0123] If there are no points in the grid cell, its elevation is marked as missing, such as NaN or 0.

[0124] In this embodiment, for the grid cells of elevation markers, the representative elevation is calculated using a Gaussian weighted quadratic polynomial fitting method, specifically including:

[0125] Step S1.3.1: For each raster cell (I, J) with missing elevation values, construct a neighborhood window with its geometric center (x0, y0) as the target point, and search for points (x0, y0) in the neighboring valid raster cells. j ,y j ,z j Let n points be found in total, forming a neighborhood point set. n is a positive integer greater than or equal to 6; a valid grid cell refers to a grid cell containing at least one point.

[0126] Step S1.3.2, Solving for the coefficients A of the quadratic polynomial, specifically:

[0127] Let the quadratic polynomial be of the form F(x,y)=a0+a1x+a2y+a3x 2 +a4xy+a5y 2 ;

[0128] Based on the neighborhood point set formed in step S1.3.1 Construct the design matrix X, weight matrix W, and neighborhood point set elevation vector Z, and then apply the normal equation (X... T WX))A=X T WZ solves for the coefficient A.

[0129] Where Z = [z1, z2, ..., z n ] T A = [a0, a1, a2, a3, a4, a5] T ,

[0130] W = diag(w1, w2, ..., w n ),

[0131]

[0132] In the above formula, R is the smoothing coefficient;

[0133] Step S1.3.3: Substitute the coefficient A obtained in step S1.3.2 into the quadratic polynomial, and then substitute the target point (x0, y0) into the fitting function to calculate the interpolated elevation value z0 of the grid cell.

[0134]

[0135] Finally, the representative elevations of all grid elements are obtained to form an elevation matrix DEM, denoted as .

[0136] Step S2: Based on the elevation model (DEM), identify natural obstacle areas in the target region and generate an impassable area matrix map. Specifically, this includes:

[0137] Step S2.1, Slope Identification: Since roads cannot have negative slopes starting from the starting point, grids with elevations lower than the starting point elevation are set as impassable. The slope s of each grid cell (I,J) is calculated according to formula (4). I,J The elevation values ​​of this grid cell and its adjacent grid cells are used to calculate the following:

[0138]

[0139] If the slope is s I,J≥s th If the slope s is not clear, then the grid unit is marked as impassable, i.e., an obstacle area; if the slope s is not clear, then the grid unit is marked as impassable, i.e., an obstacle area. I,J <s th If so, the grid unit is marked as passable. th In this embodiment, s is set as the slope threshold. th =0.3, meaning the slope is ≥30%.

[0140] In the above formula, d I,J→M,N This represents the horizontal distance between the current raster cell (I,J) and its neighboring raster cells (M,N), where M∈[0,n]. x ), N∈[0, n y For the four-neighborhood, take d. I,J→M,N =G; for the eight-neighborhood, take

[0141] Step S2.2, River Identification: For a certain continuous region Ω, if the standard deviation of its internal elevation satisfies formula (5), it is determined to be a candidate area for water bodies;

[0142]

[0143] Then, the connected components of the raster graph are used to analyze the candidate water areas. Areas that extend in a strip and connect the upper and lower reservoirs of the pumped storage power station are identified as river channels, which are also impassable areas, i.e., obstacle areas.

[0144] In the above formula, σ th The set elevation standard deviation threshold; Ω represents the continuous region, referring to the set of spatially adjacent raster cells in the DEM; |Ω| represents the number of raster cells in the continuous region; This represents the average elevation value of a continuous area.

[0145] Step S2.3: Map the identified obstacle areas onto a two-dimensional grid map to generate an obstacle matrix map; if a grid cell (I,J) is located in a passable area, then its corresponding map(I,J) = 1; if a grid cell is located in an impassable area, then its corresponding map(I,J) = 0.

[0146] The restricted areas will be included in the map as additional barrier areas. These restricted areas include, but are not limited to, construction control zones and protection zones.

[0147] Step S3: Construct a heuristic function that comprehensively considers distance and slope, an enhanced cost function, and a multi-layer neighborhood search strategy. Based on the improved A* algorithm, perform path search and construct an initial path. Specifically, this includes...

[0148] Step S3.1: Construct the mapping between grid cells and nodes: Define a node m = (I, J, g(m), h(m), f(m), parent), where each grid cell corresponds to a unique node. The node inherits the slope and passability attribute information of the corresponding grid cell. Here, (I, J) represents the row and column index of the grid cell; f(m) represents the cost function value, g(m) is the actual cost, and h(m) is the heuristic function; parent represents the parent node pointer, which is used for path backtracking.

[0149] Step S3.2: Construct a heuristic function h(m) that comprehensively considers the adaptability of distance and slope, as shown in equation (6):

[0150] h(m)=α·d norm (m)+β·Δs norm (m)(6)

[0151] In the above formula, This represents the normalized Euclidean distance from the current node m to the target node, which is the Euclidean distance d(m,goal) from the current node m to the target node divided by the maximum Euclidean distance D from the starting point to the target node. max ; This indicates the slope deviates from the normalized value; Δs(m) = |s(m) - s opt | represents the node slope s (m) and the recommended design slope s opt The degree of deviation, in this embodiment, s opt Set to 5%; s max The maximum allowable slope is denoted as 9% in this embodiment; s(m) represents the slope of the current node, obtained by equation (7), which is the difference in elevation between the current node m and its parent node p divided by the horizontal distance.

[0152]

[0153] α and β represent distance weight and slope weight, respectively, and their sum is 1. They are adjustable according to engineering requirements, and the recommended initial values ​​are α = 0.6 and β = 0.4.

[0154] In this embodiment, a dynamically adjustable weighting mechanism is introduced when constructing the heuristic function. Specifically, when the distance between the current node and the target node is less than 20% of the total path length, and the average slope tends to 5% of the design value of the target area, the system automatically shifts the weight towards the distance term, prioritizing the selection of the shortest path for rapid convergence and avoiding unnecessary detours at the path terminal.

[0155] Step S3.3: Introduce the slope penalty cost function and construct the enhanced cost function g(m), whose expression is shown in equation (8):

[0156] g(m)=g(p)+d(m,p)+γ·q(m) (8)

[0157] In the above formula, g(m) represents the total cost of the current node; d(m,p) represents the actual distance between the current node m and its parent node; q(m) represents the slope penalty function, defined as:

[0158]

[0159] In the above formula, γ is the penalty factor adjustment parameter; s opt The desired slope value is usually taken as 5%; s max To ensure maximum permissible slope, pumped storage power station projects are typically limited to 9%.

[0160] Step S3.4: For complex path planning areas, such as mountainous terrain, conventional 4-neighbor or 8-neighbor expansion is prone to getting stuck in local optima, especially in narrow areas such as canyons and steep cliffs where it is difficult to escape the current search range. Therefore, this invention adopts a multi-layer neighborhood search strategy, specifically: defining a multi-layer neighborhood node set N for the current node m. r (I,J):

[0161] N r (I,J)={(M,N)||MI|≤r,|NJ|≤r}(10)

[0162] In the above formula, r represents the number of neighborhood layers, set to 1 to 10; (M,N) is the grid cell index corresponding to the neighborhood node; (I,J) is the grid cell index corresponding to the current node m.

[0163] In this example, a three-layer neighborhood search strategy is adopted, such as Figure 2 As shown:

[0164] The first layer consists of 8 neighborhoods, adjacent to the grid cell where the current node is located; the second layer expands to 16 neighborhoods, jumping 2 grid cells from the current node to this layer; the third layer expands to 24 neighborhoods, jumping 3 grid cells from the current node to this layer.

[0165] During each expanded search, all 48 nodes in the three-layer neighborhood are included in the candidate set, and their cost functions g(m) are uniformly evaluated and ranked. This strategy improves the global path search capability, effectively avoids getting stuck in local dead ends, and can flexibly bypass large-area obstacles or areas with abrupt terrain changes.

[0166] Step S3.5: Based on steps S3.1-S3.4, perform a path search. For example... Figure 3 As shown, it specifically includes:

[0167] Step S3.5.1: Determine the starting node and the target node, create an open list to store the nodes to be explored, create a closed list to store the nodes that have been explored, and set a parent node pointer for each node to track the path;

[0168] Step S3.5.2: Add the starting node to the open list and set its total cost to 0. Check if the initial node and the target node overlap. If they overlap, end the search directly. If they do not overlap, repeat the following steps until the target is found or the open list is empty:

[0169] Step S3.5.2.1: Select the node m with the smallest total cost f(m) in the open list as the current node, f(m) = g(m) + h(m), and determine whether the current node is the target:

[0170] If so, end the search;

[0171] Otherwise, remove the current node m from the open list and add it to the closed list;

[0172] Step S3.5.2.2: Traverse the multi-level neighborhood set N of the current node m. r Extract all child nodes in (I,J) and determine whether each child node satisfies the preset constraints (such as whether it is in an obstacle area); if not, delete the child node; if it satisfies, continue to the next step.

[0173] Step S3.5.2.3: Calculate the total generation value of the child nodes that meet the conditions, and store these child nodes in the open list in ascending order. Update the open list, and then start the judgment again from step S3.5.2.1.

[0174] Step S3.5.3: Starting from the target, backtrack along the parent node pointer to the starting point to construct a complete initial path.

[0175] In this embodiment, a bidirectional A* search algorithm is used, specifically: a starting point P0 (as in the lower database) and an ending point P (as in the upper database) are set. t Given the location coordinates in a two-dimensional raster map, execute the A* search algorithm simultaneously from the starting point P0 and the ending point P0. t Proceed with a path search towards the center until two search paths meet, thus obtaining the forward path and the reverse path respectively. The meeting node is denoted as p. z The forward and reverse paths are merged to obtain the initial path {P0,…,P}. z ,…,P t}, and ensure the passability of the corresponding grid unit Pz, i.e., map[I z J z ] = 1, such as Figure 4As shown in the figure, the blue dashed line represents the initial path formed in the two-dimensional raster coordinate map.

[0176] In the iterative process of the forward search, each iteration checks whether the pre-set number of steps has been reached. If the step limit has been reached, the current forward node is used as the target point for the reverse search, and the reverse search phase begins; otherwise, the forward search continues. Similarly, in the iterative process of the reverse search, each iteration checks whether the pre-set number of steps has been reached. If the step limit has been reached, the current reverse node is used as the target point for the forward search, and the forward search phase begins; otherwise, the forward search continues.

[0177] Step S4: Optimize the initial path. First, remove redundant turning nodes, then use a third-order Bézier curve to interpolate and smooth the path, resulting in the final path. Specifically, this includes:

[0178] Step S4.1: For the generated initial path node sequence {P0, P2, ..., P...} t The sliding window method is applied to optimize redundant nodes; specifically, the initial path points {P0, P2, ..., P} are traversed. t For each path node P k Calculate vector and The included angle θ k , k∈(0,t), and compared with the set included angle threshold θ0;

[0179] If θ k If ≤θ0, then delete the path node P. k This is to reduce redundant inflection points.

[0180] If θ k If θ > 0, then retain the path node P. k , serving as a path curve control node.

[0181] Arrange the starting node, the target node, and all control nodes in sequence to form the critical path node set {P0′, P1′, ..., P′}. t′}

[0182] Step S4.2: Interpolate the critical path nodes selected in step S4.1 using a third-order Bézier curve to eliminate the inflection points of the broken line, thereby generating a continuous and smooth path curve; specifically, starting from the initial node P0′, sequentially select P... l ′,P′ l+1 ,P′ l+2 ,P′ l+3 Grouping four consecutive nodes together, construct a third-order Bézier curve as shown in equation (11), and obtain the corresponding interpolation nodes. l∈[0,3,6,…3n'], where n' is a positive integer, to achieve path segment interpolation and smooth transition;

[0183]

[0184] In the above formula, For parameters, Control the interpolation progress from the starting point to the target.

[0185] The starting node, target node, all control nodes, and corresponding interpolation nodes are arranged in sequence to form a secondary optimization path.

[0186] Step S4.3: Based on the obstacle matrix generated in step S2, verify each node of the secondary optimized path to determine whether it traverses impassable areas. If non-compliant points exist, perform local backtracking or forced avoidance processing to finally obtain the final path {P0″, P″1, ..., P″}. t}

[0187] Step S5: Overlay the optimized final path onto the original terrain raster map for visualization, calculate and display the key design indicators of the path, and assist designers in comparing different options and making engineering decisions. Specifically, this includes:

[0188] Step S5.1: Based on the final path node set {P0″, P″1, ..., P″} t The key design metrics are calculated using the following formula:

[0189] 1) Total path length L total :

[0190]

[0191] 2) Average slope S avg :

[0192]

[0193] 3) Maximum slope S max :

[0194]

[0195] In the above formula, u∈[0, t], node P u The three-dimensional coordinates of "″ are (x u ,y u ,z u );

[0196] Step S5.2: Compare the above key design indicators with the established road design specifications to determine whether they meet the design specifications, and mark and display the determination results;

[0197] Step S5.3: When multiple comparison paths are obtained by optimizing the parameters, they are scored and ranked according to key design indicators to assist designers in making the best decision. Specifically:

[0198] Score = χ1L total +χ2S avg +χ3S max ; χ1+χ2+χ3=1, where χ1, χ2, and χ3 represent the weights corresponding to the key design indicators respectively.

[0199] like Figure 5 The image shown is a diagram illustrating the results of automatic and manual route selection in this embodiment.

[0200] Based on the above-described preferred embodiments of the present invention, and through the foregoing description, those skilled in the art can make various changes and modifications without departing from the inventive concept. The technical scope of this invention is not limited to the contents of the specification, but must be determined according to the scope of the claims.

Claims

1. A road planning method between upper and lower reservoirs of a pumped storage power station based on an improved A* algorithm, characterized in that, Includes the following steps: Step S1: Rasterize the three-dimensional topographic point cloud data of the target area of ​​the pumped storage power station to generate a digital elevation model (DEM). Step S2: Based on the elevation model DEM, identify natural obstacle areas in the target area and generate an impassable area matrix map; Step S3: Construct a heuristic function that comprehensively considers distance and slope, an enhanced cost function, and a multi-layer neighborhood search strategy; perform path search based on the improved A* algorithm; and construct an initial path. Specifically, it includes: Step S3.1: Construct grid cell and node mapping: Define node m=(I,J,g(m),h(m),f(m),parent), each grid cell corresponds to a unique node, and the node inherits the slope and passability attribute information of the corresponding grid cell; Where (I,J) represents the row and column indices of the raster cell; f(m) represents the cost function value, g(m) is the actual cost, and h(m) is the heuristic function; parent represents the parent node pointer, used for path backtracking; Step S3.2: Construct a heuristic function that comprehensively considers the adaptability of distance and slope. For example, equation (6): （6） In the above formula, , represents the normalized Euclidean distance from the current node m to the target node, which is the Euclidean distance from the current node m to the target node. Divide by the maximum Euclidean distance from the starting point to the target ; This indicates that the slope deviates from the normalized value; Indicates the slope of the node relative to the recommended slope The degree of deviation, Indicates the maximum permissible slope; The slope of the current node is represented by equation (7), which is the difference in elevation between the current node m and its parent node p divided by the horizontal distance. （7） These represent the distance weight and the slope weight, respectively, and their sum is 1. Step S3.3: Introduce the slope penalty cost function and construct the enhanced cost function. Its expression is as shown in equation (8): （8） In the above formula, This represents the total cost of the current node; This represents the actual distance between the current node m and its parent node; The slope penalty function is defined as follows: （9） In the above formula, This is the parameter for adjusting the penalty factor; The desired slope value is usually taken as 5%; The maximum permissible slope; Step S3.4: For complex path planning regions, a multi-level neighborhood search strategy is adopted, specifically: defining the multi-level neighborhood node set of the current node m. : （10） In the above formula, r represents the number of neighborhood layers, set to 1 to 10; (M,N) is the grid cell index corresponding to the neighborhood node; (I,J) is the grid cell index corresponding to the current node m; Step S4: Optimize the initial path. First, remove redundant turning nodes, and then use a third-order Bézier curve to interpolate and smooth the path to obtain the final path. Step S5: Overlay the optimized final path onto the original terrain raster map for visualization, calculate and display the key design indicators of the path, and assist designers in comparing schemes and making engineering decisions.

2. The method according to claim 1, characterized in that, Step S1 specifically includes: Step S1.1, Point Cloud Boundary Range Extraction and Raster Parameter Definition: Extract the measured 3D point cloud data P={(x...} of the target area of ​​the pumped storage power station. i ,y i ,z i )|i=1,2,…,Q}, using a fixed grid size G, the point cloud data of the entire terrain area is spatially discretized, and the total number of grid rows and columns is calculated based on the point cloud boundary according to formula (1); ， （1） In the above formula, where (x i ,y i (x) represents planar coordinates, z i n represents the elevation value, Q represents the total number of point clouds; x n is the total number of columns in the grid elevation matrix. y The grid represents the total number of rows in the grid elevation matrix; G represents the side length of the grid cell, x max This represents the maximum value of the point cloud data of the target region at the boundary in the x-direction, x min This represents the minimum value of the point cloud data of the target region at the boundary in the x-direction, and y-direction... max This represents the maximum value of the point cloud data of the target region at the boundary in the y-direction. min This represents the minimum value of the point cloud data of the target region in the y-direction boundary, and ⌈ ⌉ represents the rounding up operation; Step S1.2: Map the point cloud to the raster index: According to formula (2), calculate the index of each point ( Map each point to its corresponding raster index (I) in both the x and y directions. i J i That is, the point cloud is assigned to the corresponding raster cells, and each raster cell (I, J) is associated with a set of points; ， （2） In the above formula, x i Let y be the X-coordinate of the i-th point. i Let I be the Y-coordinate of the i-th point; k I represents the column index of the grid cell. i ∈[0, n x ); J i J represents the row index of the grid cell. i ∈[0, n y ), ⌊⌋ represents the floor operation; Step S1.3: Generate Elevation Matrix Model (EDM): Calculate the elevation of the points falling within each raster cell (I, J) to form the corresponding elevation set. If there is at least one point in the grid cell, then the representative elevation H(I, J) of the grid cell is as shown in equation (3): （3） In the above formula, Z I,J H(I,J) represents the set of elevation values ​​for all points within a raster cell; H(I,J) represents the representative elevation of the raster cell; |ZI,J| represents the number of elements in the elevation set. If there are no points in the grid cell, its elevation is marked as missing and repaired by interpolation to represent the elevation. Finally, the representative elevations of all grid elements constitute an elevation matrix (DEM), denoted as . .

3. The method according to claim 2, characterized in that, In step S1.3, for raster cells with missing elevation values, the representative elevation is calculated using a Gaussian weighted quadratic polynomial fitting method, specifically including: Step S1.3.1: For each raster cell (I, J) with missing elevation values, take its geometric center... Using the target point as a reference, a neighborhood window is constructed to search for points in the nearest valid raster cells. Set up a total search to find a total of Each point forms a neighborhood point set. The term "valid grid cell" refers to a grid cell containing at least one point. Step S1.3.2, Solving for the coefficients A of the quadratic polynomial, specifically: Let the quadratic polynomial form be ; Based on the neighborhood point set formed in step S1.3.1 Construct the design matrix X, weight matrix W, and neighborhood point set elevation vector Z, and then use the normal equation (X) T WX))A=X T WZ solves for coefficient A; Where Z=[z1,z2,…,z] n ] T A = [a0, a1, a2, a3, a4, a5] T , W=diag(w1,w2,…,w n ), ; ； In the above formula, R is the smoothing coefficient; Step S1.3.3: Substitute the coefficient A obtained in step S1.3.2 into the quadratic polynomial, and then set the target point... Substitute the values ​​into the fitting function to calculate the interpolated elevation value of the raster cell. ; 。 4. The method according to claim 3, characterized in that, In step S2, natural obstacles are identified using the following method, specifically including: Step S2.1, Slope Identification: Calculate the slope of each grid cell (I,J) according to formula (4). The elevation values ​​of this grid cell and its adjacent grid cells are used to calculate the following: (4)； In the above formula, This represents the horizontal distance between the current raster cell (I,J) and its neighboring raster cells (M,N), where M∈[0,n]. x ), N∈[0, n y ); For the four neighboring regions, take For the eight neighboring domains, take ; If the slope If the grid unit is marked as impassable, i.e., an obstacle zone, then the grid unit is marked as impassable. The set slope threshold; If the slope If so, then mark the grid unit as passable; Step S2.2, River Identification: For a certain continuous area If the standard deviation of its internal elevation satisfies formula (5), then it is determined to be a candidate area for water body; （5） Then, the connected components of the raster graph are used to analyze the candidate water bodies. Areas that extend in a strip and connect the upper and lower reservoirs of the pumped storage power station are identified as river channels, i.e., obstacle areas. In the above formula, The set elevation standard deviation threshold; Ω represents the continuous region, referring to the set of spatially adjacent raster cells in the DEM; |Ω| represents the number of raster cells in the continuous region; This represents the average elevation value of a continuous area; Step S2.3: Map the identified obstacle areas onto a two-dimensional grid map to generate an impassable area matrix map; If a grid cell (I,J) is located in a passable area, then its corresponding map(I,J)=1; If a grid cell is located in a non-passable area, then its corresponding map(I,J)=0.

5. The method according to claim 4, characterized in that, In step S2.3, the restricted area is included in the map as an additional barrier area. The restricted area includes, but is not limited to, the construction control area and the protection area.

6. The method according to claim 5, characterized in that, Step S3, which performs path search based on the improved A* algorithm, also includes: Step S3.5: Based on steps S3.1-S3.4, perform path search, specifically including: Step S3.5.1: Determine the starting node and the target node, create an open list to store the nodes to be explored, create a closed list to store the nodes that have been explored, and set a parent node pointer for each node to track the path; Step S3.5.2: Add the starting node to the open list and set its total cost to 0. Then, repeat the following steps until the target is found or the open list is empty: Step S3.5.2.1: Select the node m with the smallest total cost f(m) in the open list as the current node, f(m) = g(m) + h(m), and determine whether the current node is the target: If so, end the search; Otherwise, remove the current node m from the open list and add it to the closed list; Step S3.5.2.2: Traverse the multi-level neighborhood set of the current node m. The algorithm identifies all child nodes and determines whether each child node satisfies the preset constraints. If not, the child node is deleted; otherwise, the algorithm continues to the next step. Step S3.5.2.3: Calculate the total generation value of the child nodes that meet the conditions, and store these child nodes in the open list in ascending order. Update the open list, and then start the judgment again from step S3.5.2.

1. Step S3.5.3: Starting from the target, backtrack along the parent node pointer to the starting point to construct a complete initial path.

7. The method according to claim 6, characterized in that, In step S3.5, a bidirectional A* search algorithm is used, specifically: setting the starting point P0 and the target P0. t The improved A* search algorithm is adopted to simultaneously search from the starting point P0 and the target P. t Proceed with a path search towards the center until two search paths meet, thus obtaining the forward path and the reverse path respectively. The meeting node is denoted as p. z The forward and reverse paths are merged to obtain the initial path {P0,…,P}. z ,…,P t }, and ensure the passability of the corresponding grid unit Pz, i.e., map[I z J z ]=1.

8. The method according to claim 7, characterized in that, In step S4, the initial path is optimized, specifically including: Step S4.1: For the generated initial path node sequence {P0, P2, ..., P...} t The sliding window method is applied to optimize redundant nodes; specifically, the initial path points {P0, P2, ..., P} are traversed. t For each path node P k Calculate vector and The included angle , k∈(0,t), and within the set included angle threshold Compare; like Then delete the path node P. k ; like > Then retain the path node P. k , serving as a path curve control node; Arrange the starting node, the target node, and all control nodes in sequence to form the critical path node set. ; Step S4.2: Interpolate the critical path nodes selected in Step S4.1 using a third-order Bézier curve to eliminate the inflection points of the broken line, thereby generating a continuous and smooth path curve; specifically, starting from the starting node... Start by taking... As a group, construct the third-order Bézier curve as shown in equation (11) to obtain the corresponding interpolation nodes. l∈[0,3,6,…3n'], where n' is a positive integer; (11) In the above formula, For parameters, Control the interpolation progress from the starting point to the target; Arrange the starting node, target node, all control nodes, and corresponding interpolation nodes in sequence to form a secondary optimization path; Step S4.3: Based on the obstacle matrix generated in step S2, verify each node of the secondary optimized path to determine whether it traverses impassable areas. If non-compliant points exist, perform local backtracking or forced avoidance to obtain the final path. .

9. The method according to claim 8, characterized in that, Step S5 specifically includes: Step S5.1: Based on the final path node set The key design metrics are calculated using the following formula: 1) Total path length : （12） 2) Average slope : （13） 3) Maximum slope : （14） In the above formula, u∈[0, t], node The three-dimensional coordinates are (x u ,y u ,z u ); Step S5.2: Compare the above key design indicators with the established road design specifications to determine whether they meet the design specifications, and mark and display the determination results; Step S5.3: When multiple comparison paths are obtained by optimizing the parameters, they are scored and ranked according to key design indicators to assist designers in making the best decision. Specifically: ; , This indicates the weights corresponding to each key design metric.

Images