Multi-terrain path intelligent planning method fusing high-resolution satellite remote sensing data

Abstract

The application discloses a multi-terrain path intelligent planning method fusing high-resolution satellite remote sensing data, and belongs to the technical field of remote sensing data analysis and path planning. The method acquires target mountainous area high-resolution satellite remote sensing data, extracts a digital surface model and an orthographic image; terrain analysis is performed on the digital surface model to generate a slope graph, a slope direction graph and a curvature graph; ground cover types are classified and identified according to the orthographic image; the slope graph, the slope direction graph, the curvature graph and the ground cover types are fused to construct a multi-terrain resistance grid graph; a graph cut algorithm is used to convert path planning into a graph cut problem, a graph cut energy function is constructed on the passing resistance grid graph, a minimum cumulative resistance path is solved as an initial passing path; morphological inflation is performed on the initial passing path to generate a path buffer zone, and local smoothing processing is performed on the path buffer zone by using the digital surface model to obtain a final passing path.

Description

Technical Field

[0001] This invention relates to the field of remote sensing data analysis and path planning technology, specifically to a multi-terrain path intelligent planning method that integrates high-resolution satellite remote sensing data. Background Technology

[0002] Route planning in mountainous areas with complex terrain and surface environments presents a significant challenge in remote sensing applications and spatial analysis. Traditional route planning methods typically rely on digital elevation models (DEMs) to extract slope or single terrain factors to construct a route cost map, then use graph search algorithms to generate routes connecting the starting and ending points. These approaches neglect the rich land cover information contained in high-resolution remote sensing data, making it difficult to impose differentiated traffic resistance constraints on different underlying surface types such as forests, bare rock, glaciers, and snow cover. Furthermore, they fail to comprehensively represent the combined effects of slope, aspect, surface curvature, and land cover type, resulting in routes with poor traversability and high safety risks under actual mountainous conditions. In terms of path solving, existing methods mostly adopt point-by-point expansion search strategies such as A* or Dijkstra. These algorithms are inefficient in handling a large number of high-resistance or impassable areas in the traffic resistance grid. The generated paths tend to adhere to the edges of obstacles, resulting in rigid and tortuous paths. Furthermore, they lack re-analysis and adaptive smoothing of the local terrain of the path and do not utilize high-precision digital surface models to make local adjustments to the path at the micro-topographic level. They are unable to effectively avoid unfavorable micro-topographic features such as steep slopes and grooves on both sides of the path. Summary of the Invention

[0003] This invention provides a multi-terrain path intelligent planning method that integrates high-resolution satellite remote sensing data to solve the problems of insufficient collaborative resistance modeling of multiple factors of terrain and land cover in existing technologies and the lack of local terrain adaptive optimization in path generation results.

[0004] The objective of this invention can be achieved through the following technical solutions:

[0005] This invention discloses a multi-terrain path intelligent planning method that integrates high-resolution satellite remote sensing data, comprising:

[0006] High-resolution satellite remote sensing data of the target mountainous area is acquired, and a digital surface model and orthophotos are extracted from it. Topographic analysis is performed on the digital surface model to generate slope, aspect, and curvature maps; land cover types are identified based on the orthophotos. The slope, aspect, curvature, and land cover types are fused to construct a multi-topographic resistance raster map. A graph cut algorithm is used to transform the path planning problem into a graph cut problem. Starting and target points are set on the traffic resistance raster map, and a graph cut energy function is constructed at the cost of traffic resistance. The path with the minimum cumulative resistance from the starting point to the target point is obtained as the initial traffic path. Morphological dilation is applied to the initial traffic path to generate a path buffer. Local smoothing is performed within the buffer using the digital surface model to obtain the final traffic path.

[0007] As a technical solution of this invention, when extracting the digital surface model, a stereo image pair of the target mountain area acquired by a high-resolution satellite is received. A semi-global dense matching algorithm is used to perform pixel-level matching of the forward-looking and backward-looking images, calculate the disparity value of corresponding pixels, and generate an initial digital surface model. After median filtering to remove isolated noise points, the final digital surface model is obtained. Simultaneously, multispectral satellite images are received, and surface reflectance images are obtained after radiometric calibration and atmospheric correction. Orthorectification is then performed using the digital surface model to eliminate geometric deformation caused by terrain, resulting in an orthorectified image. This process ensures the geometric consistency and radiometric quality of the basic topographic data and image data from the source, providing a reliable data foundation for subsequent accurate analysis.

[0008] In a preferred implementation, during terrain analysis, each grid cell in the digital surface model is traversed, and the elevation values ​​of nine grid cells within its nine-square neighborhood window are extracted. A third-order finite difference operator is used to calculate the elevation change rates in the east-west and north-south directions, and based on this, the slope and aspect values ​​of the central grid cell are calculated and filled into the corresponding positions on the slope and aspect maps, respectively. A second-order difference operation is performed on the slope map to obtain the curvature map. This method can finely characterize micro-topographic features, accurately reflect the steepness of mountains, orientation distribution, and surface undulations, providing high-precision terrain parameters for traffic resistance modeling.

[0009] When identifying land cover types, spectral feature vectors containing reflectance in blue, green, red, and near-infrared bands are extracted from orthophotos for each pixel. Gray-level co-occurrence matrix analysis is then performed on local neighborhoods to extract texture features such as contrast, correlation, and homogeneity. The spectral and texture features are concatenated and fed into a pre-trained radial basis function (RBF) support vector machine classifier. The output classifier calculates the probability that a pixel belongs to a predefined type, such as coniferous forest, broadleaf forest, alpine meadow, bare rock, glacier, snow cover, or water body. The type with the highest probability is selected as the identification result. This classification strategy integrates multispectral and spatial texture information, effectively improving the discrimination accuracy of complex mountainous land cover features and avoiding confusion caused by single spectral features.

[0010] When constructing a multi-topographic resistance raster map, the initial passability resistance value of each raster cell is composed of four parts: slope value, aspect deviation value, curvature correction value, and baseline resistance value for land cover type. For raster cells with slopes exceeding a preset upper limit, their resistance is set to infinity to indicate impassability. For raster cells within the limit, the angle between the slope aspect and the prevailing wind direction of the region is calculated, and an additional resistance term is added when the angle exceeds a preset angle threshold. For concave areas with negative curvature, the resistance is multiplied by a concave resistance amplification factor greater than one. A corresponding baseline resistance value is assigned to each land cover type, and the baseline resistance for water bodies is set to infinity. Through multi-factor fusion, this raster map comprehensively reflects the steepness of slope, headwind influence, risk of subsidence in concave terrain, and the ease of passage through different underlying surfaces, ensuring that path planning fully considers real-world traffic constraints and risk factors. Preferably, the threshold for the angle between the slope aspect and the prevailing wind direction is set as a configurable parameter between 45 and 90 degrees to adapt to path selection preferences under different wind conditions.

[0011] In the graph cut solution stage, each grid cell in the traffic resistance raster graph is constructed as a graph node, and the spatial adjacency relationship between adjacent grid cells is used as an edge. A binary label variable is assigned to each graph node, and an energy function consisting of a data term and a smoothing term is constructed. The data term makes nodes with lower resistance more likely to be marked as traffic paths, while the smoothing term imposes a higher penalty for dissimilarity when adjacent nodes have large resistance differences or large spatial distances. The starting point is forcibly connected to the source point, and the target point is forcibly connected to the sink point, with the capacity of the corresponding edges set to infinity. After solving for the minimum cut, the sequence of grid cells marked as traffic paths and connected to the starting and ending points is extracted as the initial traffic path. This scheme transforms global path optimization into a graph cut problem, efficiently finding continuous paths that satisfy the minimum cumulative resistance, avoiding the problem of traditional greedy search getting trapped in local optima. A further preferred approach is to use an augmenting path and residual network iterative method to solve for the minimum cut. Starting from the source node, a hierarchical graph is constructed and augmenting paths are continuously searched. The residual capacity is updated until the sink node can no longer be reached. The edges between the source node set and the sink node set are the minimum cuts, and then the connected paths are extracted. This process ensures the global optimality of the solution.

[0012] To improve planning quality over a larger spatial scale, a multi-scale iterative optimization strategy can be adopted: Gaussian pyramid downsampling is performed on the original traffic resistance raster map to generate multiple scale layers; starting from the coarsest scale, graph cut solutions are performed to obtain the initial path region, which is then upsampled and mapped to a finer-scale layer using bilinear interpolation as prior information to constrain graph cut solutions. This process is repeated layer by layer until the solution is completed on the original resolution layer, extracting the initial traffic path. This strategy utilizes the coarse scale to quickly determine the approximate path direction, while the fine scale performs fine adjustments under prior constraints, significantly reducing computation time on high-resolution maps while maintaining path continuity and global optimum characteristics.

[0013] Based on the initial travel path, when generating the path buffer, the initial path point set is obtained, and the local terrain undulation of each point is calculated according to the digital surface model. The radius of the circular expansion structure element is dynamically determined based on the terrain undulation; the larger the undulation, the larger the expansion radius, ranging from one to ten grid cells. The expansion operation is performed on all points, and the covered areas are merged to form the path buffer. This dynamic expansion mechanism enables a wider buffer search band in undulating and unstable terrain areas, providing sufficient terrain information for subsequent smoothing, while avoiding unnecessary expansion of the buffer range in flat areas. Preferably, the positive correlation between terrain undulation and expansion radius is defined using a piecewise linear function. When the terrain undulation is less than a first threshold, the radius takes its minimum value; when it is greater than a second threshold, the radius takes its maximum value; and when it is in between, it increases linearly, allowing the expansion range to flexibly match the terrain complexity.

[0014] When performing local smoothing within the buffer zone, the elevation point cloud of the buffer zone is extracted. For each path point, its normal direction on the digital surface model is calculated, and a one-dimensional search is performed within the buffer zone along the normal direction to find the raster cell with the lowest elevation value as a candidate replacement point. All candidate replacement points are connected in the order of the initial path to obtain a preliminary smoothed path. Then, cubic spline interpolation is used to generate a continuous curve and convert it back to a raster cell sequence to obtain the final travel path. The one-dimensional search range is limited to the search interval formed by extending a preset number of steps on both sides of the normal direction with the current path point as the center. This smoothing process guides the path to the point of lowest local elevation in the normal direction without leaving the safe buffer zone, effectively avoiding protruding micro-terrain obstacles and reducing the amount of undulation in actual travel. Cubic spline interpolation ensures the continuity and smoothness of the path curve, improving the feasibility of hiking or trekking.

[0015] This invention utilizes deep interpretation of high-resolution satellite remote sensing data to fuse fine-grained terrain factors and land cover attributes at multiple levels, constructing a raster map that accurately represents various terrain resistances. By combining global map cut optimization and adaptive buffer local smoothing, it generates safe, efficient, and realistic travel paths. The entire method is highly automated, adaptable to the complex and varied geographical environments of mountainous areas, and provides scientifically reliable route support for field exploration, emergency search and rescue, and border patrols.

[0016] The beneficial effects of this invention are:

[0017] The path planning problem is transformed into a graph cut problem. On a traffic resistance raster map constructed by integrating slope, aspect, curvature, and land cover type, a graph cut energy function containing data and smoothing terms is built. The data term is set according to the traffic resistance value of each raster cell; the lower the resistance value, the lower the cost of classifying the raster cell as a path region. The smoothing term defines a penalty cost for inconsistent labeling based on the difference in traffic resistance and spatial distance between adjacent raster cells; the more significant the difference in resistance or the larger the spatial span of adjacent raster cells, the heavier the penalty for assigning different labels. The minimum cut of this graph cut model is solved, extracting low-resistance connected regions connected to the starting and target points from the partitioning of the source and sink sets as the initial traffic path. This global optimization mechanism can simultaneously balance traffic resistance and path smoothness throughout the entire raster space, avoiding paths that are too close to the edges of high-resistance obstacles or have redundant turns in narrow passages, outputting paths with low cumulative resistance and smooth morphology under complex resistance distributions. When performing morphological dilation on the initial travel path, the radius of the circular dilation structural element is dynamically determined based on the terrain undulation at the location of each path point. Terrain undulation is defined as the difference between the highest and lowest elevations within a preset window centered on that point. The more pronounced the terrain undulation, the larger the dilation radius, with a lower limit of one raster cell and an upper limit of ten raster cells. The path buffer generated by the dilation captures micro-terrain information related to travel around the path. During local smoothing within the buffer, a one-dimensional search is performed along the normal direction of each path point on the digital surface model. The raster cell with the lowest elevation value within the search interval on both sides of the normal direction is selected as a candidate replacement point. After connecting all candidate replacement points in the original path order, cubic spline interpolation is used to generate a continuous final travel path. This process allows the path to actively migrate towards locally low-lying and flat areas while preserving the globally optimal orientation, avoiding unfavorable terrain such as cross slopes and rock embankments. This reduces the longitudinal profile undulation of the path and significantly improves path continuity and traversability. Attached Figure Description

[0018] The invention will now be further described with reference to the accompanying drawings.

[0019] Figure 1 This is a flowchart of a multi-terrain path intelligent planning method that integrates high-resolution satellite remote sensing data;

[0020] Figure 2 This is a flowchart of digital surface model and orthophoto generation and terrain analysis;

[0021] Figure 3 This is a schematic diagram illustrating the process of constructing a multi-terrain resistance raster map;

[0022] Figure 4 This is a flowchart of path smoothing processing based on terrain-adaptive morphological dilation. Detailed Implementation

[0023] 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.

[0024] See Figure 1 This invention provides a multi-terrain path intelligent planning method that integrates high-resolution satellite remote sensing data. The overall implementation scheme of this method is as follows:

[0025] High-resolution satellite remote sensing data of the target mountainous area is acquired, and digital surface models and orthophotos are extracted from the high-resolution satellite remote sensing data. Topographic analysis is performed on the digital surface models to generate slope maps, aspect maps, and curvature maps. Land cover types are identified based on the orthophotos. The slope maps, aspect maps, curvature maps, and land cover types are fused to construct a multi-topographic resistance raster map. A graph cut algorithm is used to transform the path planning problem into a graph cut problem. Starting and target points are set on the traffic resistance raster map, and a graph cut energy function is constructed at the cost of traffic resistance. The path with the minimum cumulative resistance from the starting point to the target point is obtained as the initial traffic path. Morphological dilation is applied to the initial traffic path to generate a path buffer. Local smoothing is performed within the buffer using the digital surface model to obtain the final traffic path. This method's multi-topographic path intelligent planning refers to the automatic generation and optimization of traffic paths in the specific geographical environment of mountainous areas, comprehensively considering the combined effects of multiple topographic factors and land cover types.

[0026] In specific implementation, please refer to Figure 2 The study acquires high-resolution satellite remote sensing data of the target mountainous area, from which digital surface models and orthophotos are extracted. High-resolution satellite remote sensing data is acquired by receiving stereo image pairs of the target mountainous area from high-resolution satellites. These stereo image pairs contain forward-looking and backward-looking images of the same area acquired from different observation angles. A semi-global dense matching algorithm is used to perform pixel-level matching on the forward-looking and backward-looking images. The algorithm calculates the aggregated matching cost value for each pixel along multiple directions and determines the disparity value of each corresponding pixel by minimizing the energy function.

[0027] The initial digital surface model is subjected to median filtering. The median filtering process uses a sliding window to traverse each grid cell in the initial digital surface model. The elevation value of the central grid cell in the sliding window is replaced with the median of the elevation values ​​of all grid cells in the sliding window to filter out isolated noise points and obtain the final digital surface model.

[0028] Multispectral satellite imagery of the target mountainous area was received simultaneously. This imagery includes data from multiple bands, including blue, green, red, and near-infrared. Radiometric calibration was performed on the multispectral satellite imagery, converting image pixel grayscale values ​​to apparent radiance at the top of the atmosphere. The calibration coefficients used in the radiometric calibration process were obtained from the satellite imagery metadata. Atmospheric correction was then applied to the radiometrically calibrated imagery. Atmospheric correction employed a correction method based on a radiative transfer model to eliminate the influence of atmospheric absorption and scattering on the surface reflectance signal, resulting in a surface reflectance image.

[0029] Orthorectification of surface reflectance images is performed using a digital surface model as an elevation reference. Based on the ground point elevation information provided by the digital surface model and the orbital attitude parameters during satellite imaging, a collinearity equation model is employed to geometrically resample each pixel in the surface reflectance image, eliminating geometric distortions caused by terrain and generating an orthorectified image.

[0030] Topographic analysis is performed on the digital surface model to generate slope, aspect, and curvature maps. Each raster cell in the digital surface model is traversed, and a 3x3 grid neighborhood window is extracted centered on the currently traversed raster cell. The elevation values ​​of the nine raster cells within this window are recorded. A third-order finite difference operator is used to calculate the east-west and north-south elevation change rates of the central raster cell. The east-west elevation change rate is obtained by dividing the difference in elevation values ​​between the left and right sides of the central column by twice the raster spacing, and the north-south elevation change rate is obtained by dividing the difference in elevation values ​​between the top and bottom sides of the central row by twice the raster spacing. Specifically, the elevation value of the central raster cell within the 3x3 neighborhood window is recorded as... The elevation values ​​of the grid cells located to the left, right, top, and bottom of the central grid cell in the nine-grid neighborhood window are respectively denoted as... , , , The grid spacing is denoted as The rate of change of elevation in the east-west direction Represented as North-South Elevation Change Rate Represented as .

[0031] Based on the east-west elevation change rate and the rate of change of elevation in the north-south direction Calculate the slope value of the central grid cell. The calculation formula is:

[0032]

[0033] in, The slope value is expressed in radians. Represents the arctangent function. This represents the square root operation. The calculated slope value is then used. Fill in the corresponding positions on the slope map.

[0034] Based on the east-west elevation change rate and the rate of change of elevation in the north-south direction Calculate the slope aspect value of the central grid cell. The calculation formula is and according to and The positive and negative relationship will Convert to an angle range of 0 to 360 degrees. Calculate the slope aspect value. Fill in the corresponding positions on the aspect diagram.

[0035] A second difference operation is performed on the slope map, using the same third-order finite difference operator as used in calculating the slope map to calculate the slope change rate of each grid cell, resulting in a curvature map. The second difference operation iterates through each grid cell in the slope map, extracting the slope value within a 3x3 grid neighborhood window, repeating the calculation process of the third-order finite difference operator, and filling the obtained slope change rate into the corresponding position in the curvature map.

[0036] In practice, land cover types are identified based on orthophotos. A spectral feature vector is extracted from each pixel of the orthophoto image. This vector contains four dimensions: blue band reflectance, green band reflectance, red band reflectance, and near-infrared band reflectance. The blue band reflectance, green band reflectance, red band reflectance, and near-infrared band reflectance are directly read from the corresponding spectral band data of the orthophoto image.

[0037] Gray-level co-occurrence matrix (GLCM) analysis is performed on the local neighborhood window of each pixel. The local neighborhood window, centered on the currently analyzed pixel, covers a pre-defined rectangular region with a width and height of 21 pixels each. Within the local neighborhood window, the gray values ​​of the near-infrared band of the orthophoto image are quantized into 32 gray levels. The GLCM of all pixel pairs within the window is calculated with a horizontal offset of 1 pixel and a vertical offset of 0 pixels. Contrast, correlation, and homogeneity features are extracted from the GLCM. Contrast features... The calculation formula is:

[0038]

[0039] in, Represents the contrast feature value. This represents the total number of gray levels. The value is 32. This represents the row index of the gray-level co-occurrence matrix. The column index of the gray-level co-occurrence matrix. Indicates row index With column index The square of the difference This indicates that the gray-level co-occurrence matrix is ​​located at the th position. Line number The element values ​​of the column. The calculation of correlation and homogeneity eigenvalues ​​is performed based on the conventional definition of the gray-level co-occurrence matrix.

[0040] The spectral feature vector, contrast feature value, correlation feature value, and homogeneity feature value are concatenated to form the joint feature input vector of the currently analyzed pixel. The concatenation operation is performed in the following order: blue band reflectance, green band reflectance, red band reflectance, near-infrared band reflectance, contrast feature value, correlation feature value, and homogeneity feature value, generating a seven-dimensional joint feature input vector.

[0041] The joint feature input vector is fed into a pre-trained Support Vector Machine (SVM) classifier. The core architecture of the SVM classifier consists of an input layer, a kernel function mapping layer, and a decision output layer. The input layer receives a seven-dimensional joint feature input vector, the kernel function mapping layer uses a radial basis function to map the input vector to a high-dimensional feature space, and the decision output layer constructs multiple binary classification decision functions based on the support vectors and Lagrange multipliers. The training process of the Support Vector Machine (SVM) classifier is based on remote sensing data of sample areas with labeled land cover types. The training process includes the following steps: In the remote sensing data of the sample areas, for each predefined land cover type, several sample areas are manually selected, and the joint feature input vector of all pixels in each sample area is extracted. Each joint feature input vector is then labeled with the corresponding land cover type label. The labeled set of joint feature input vectors is randomly divided into a training set and a validation set. The number of samples in the training set accounts for 70% of the total number of samples in the joint feature input vector set, and the number of samples in the validation set accounts for 30% of the total number of samples in the joint feature input vector set. Multiple binary SVMs are combined using a one-to-one strategy. For a classification task containing seven predefined land cover types, twenty-one binary SVMs are combined. Each binary SVM uses a radial basis function (RBF) kernel function, which takes the form of: ,in, and These represent two input vectors, Represents kernel parameters, This represents an exponential function with base e. Let represent the Euclidean norm of the vector; set the penalty parameter when training each binary support vector machine. Select values ​​from the set {0.1, 1, 10, 100} to set the kernel parameters. Values ​​are taken from the set {0.01, 0.1, 1, 10}, depending on the penalty parameter. and kernel parameters For each combination, solve the dual optimization problem of the support vector machine on the training set to obtain the support vectors and their corresponding Lagrange multipliers, and evaluate the classification accuracy on the validation set; select the penalty parameter that yields the highest classification accuracy on the validation set. and kernel parameters The combined parameters are used as the final model parameters. The system is then retrained on the set of joint feature input vectors of the remote sensing data from the entire labeled sample area to obtain the final support vector machine classifier.

[0042] The Support Vector Machine (SVM) classifier outputs a probability vector for each pixel belonging to a predefined land cover type. The probability vector is calculated using a pairwise coupling method. The output decision values ​​of the twenty-one binary SVM classes are converted into twenty-one posterior probability estimates using Platt scaling. These twenty-one posterior probability estimates are then fused using the pairwise coupling method to generate a seven-dimensional probability vector. Each dimension of the probability vector corresponds to a predefined land cover type, including coniferous forest, broadleaf forest, alpine meadow, bare rock, glacier, snow cover, and water body. The predefined land cover type with the highest probability value in the probability vector is selected as the recognition result for that pixel.

[0043] In specific implementation, please refer to Figure 3 A multi-topographic resistance raster map is constructed by integrating slope maps, aspect maps, curvature maps, and land cover types. The construction process of the traffic resistance raster map begins with calculating the initial traffic resistance value for each raster cell. The initial traffic resistance value is composed of four parts: slope value, aspect deviation value, curvature correction value, and land cover type baseline resistance value. The superposition calculation follows the formula below:

[0044]

[0045] in, Indicates the coordinates of the traffic resistance grid. The initial passage resistance value of the grid cell at that location, This indicates the coordinates read from the slope map. The slope value of the grid cell at that location, expressed in degrees. Indicates the location at coordinates The slope deviation value of the grid cell at that location. Indicates the location at coordinates The curvature correction value of the grid cell at that location. Indicates the location at coordinates The reference resistance value for the land cover type corresponding to the grid cell at that location.

[0046] Slope deviation value The calculation method is as follows: read the grid cells from the slope aspect map. Slope value at the location , slope value Compared with the pre-statistical regional prevailing wind angle value Take the difference and get the absolute value, then multiply the absolute value by the slope deviation coefficient. Obtain the slope deviation value Slope aspect deviation coefficient The setting is 1.0, based on the principle that every 1 degree deviation of the slope from the prevailing wind direction corresponds to 1 unit of initial traffic resistance. Regional prevailing wind direction angle value. It is the azimuth angle of the prevailing wind direction, which is statistically obtained from the annual wind direction observation data of the meteorological stations in the target mountainous area.

[0047] Curvature correction value The calculation method is as follows: read the grid cells from the curvature map. curvature value at curvature value Representing terrain in grid cells The degree of unevenness and bending strength at the point will affect the curvature value. The absolute value multiplied by the curvature correction factor Obtain curvature correction value Curvature correction factor The value is set to 10.0, based on the curvature value in the curvature plot. The order of magnitude is usually Each meter, after being multiplied by 10.0, has a curvature correction value on the same order of magnitude as the slope value.

[0048] Land cover type benchmark resistance value The baseline resistance values ​​were assigned based on the land cover type classification results. The baseline resistance value for coniferous forest cover was set at 50, for broadleaf forest cover at 40, for alpine meadow cover at 20, for bare rock cover at 30, for glacier cover at 60, for snow cover at 70, and for water cover at infinity. These baseline resistance values ​​were set based on prior knowledge of the degree of obstruction to pedestrian movement caused by different land cover types. The infinity baseline resistance value for water cover indicates that walking through water bodies is impossible.

[0049] After obtaining the initial traffic resistance values ​​for all grid cells, conditional checks and resistance value adjustments are performed on each grid cell. This is done for the slope value. For grid cells exceeding the preset slope limit, the passage resistance value of the grid cell will be directly set to infinity. The slope limit is set to 45 degrees, based on the fact that people cannot maintain a stable standing posture under slope conditions exceeding 45 degrees, posing a risk of slipping and falling, and are therefore considered impassable areas.

[0050] Regarding slope value For grid cells with a slope not exceeding the upper limit of 45 degrees, the slope aspect value of the grid cell is read from the aspect map. The angle between the slope aspect of the grid cell and the pre-statistically calculated prevailing wind direction of the region is calculated. When the angle exceeds the pre-set angle threshold of 45 degrees, an additional resistance term is added to the calculated initial passage resistance value. The angle threshold of 45 degrees is based on mountain meteorological statistics. When the angle between the slope aspect and the prevailing wind direction of the region exceeds 45 degrees, people on the hillside will experience significant lateral wind force or wind resistance, requiring additional penalties for passage resistance. The additional resistance term is set to a fixed value of 500. This fixed value of 500 is based on the principle that the additional resistance value should be significantly greater than the regular resistance component to highlight the tendency of strong winds to influence path selection.

[0051] Based on the curvature map, concave grid cells with negative curvature values ​​are identified. A negative curvature value indicates that the terrain is concave at that grid cell. For the identified concave grid cells, the adjusted travel resistance value of the grid cell (i.e., the resistance value after slope upper limit check and slope aspect angle addition) is multiplied by a concave resistance amplification factor. Concave drag amplification factor The value is set to 1.5 because concave terrain areas are prone to snow accumulation, water retention, and poor air circulation, resulting in significantly increased traffic resistance compared to convex or flat terrain. A magnification factor of 1.5 reasonably reflects the additional difficulty of passage in concave terrain. The traffic resistance values ​​of grid cells in non-concave areas remain unchanged. After adjusting all the above grid cells, the final multi-terrain resistance raster map is generated.

[0052] In practice, the graph cut algorithm is used to transform the path planning problem into a graph cut problem. The starting point and the target point are set on the traffic resistance grid map, and a graph cut energy function with traffic resistance as the cost is constructed. The path with the minimum cumulative resistance from the starting point to the target point is obtained as the initial traffic path.

[0053] Each grid cell in the traffic resistance grid diagram is treated as a graph node in the graph cut model, and the spatial adjacency between adjacent grid cells is treated as an edge in the graph cut model. The spatial adjacency uses an eight-neighborhood connection method, where each grid cell establishes edge connections with eight adjacent grid cells in the top, bottom, left, right, and four diagonal directions. Each graph node is assigned a binary label variable, which takes the value of either a first label or a second label. The first label indicates that the grid cell belongs to the traffic path region, and the second label indicates that the grid cell does not belong to the traffic path region.

[0054] Construct the graph cut energy function, which has the following form:

[0055]

[0056] in, Indicates the field of the label The corresponding total energy value, Let V be a vector consisting of the label variables of all graph nodes, and let V represent the set of graph nodes. Represents the set of graph nodes Any graph node in the graph, Represents graph nodes binary labeled variables, The value can be either the first label or the second label. Represents graph nodes Assigned a label The cost of data items at that time Represents the weighting coefficients of the smoothing term. Denotes the set of edges. Represents the set of edges Connecting Graph Nodes Graph Nodes One of the edges, and These represent the two graph nodes at either end of the edge. Represents graph nodes binary labeled variables, Indicates adjacent graph nodes Graph Nodes They were assigned labels. and The smoothing term penalty cost.

[0057] Smoothing term weight coefficient The setting is 2.0, which is based on balancing the magnitude of the data item cost and the smoothing item penalty cost, so that the path planning result both respects the magnitude of traffic resistance and maintains the spatial continuity of the path.

[0058] Data item cost The value is set based on the passage resistance value of each grid cell. Grid cells with lower passage resistance values ​​have lower data item costs when assigned the first label. When a graph node... When assigned the first label, the cost of the data item Set as the passage resistance value of the grid cell corresponding to this graph node; when the graph node When assigned a second label, the cost of the data item The passage resistance value of the corresponding grid cell for this node is set to a preset scaling factor. The product of the products, with a preset proportionality coefficient. The value is set to 0.3, which is based on the principle of ensuring that the penalty cost of a path not selected is lower than the cost baseline of a path selected, thus preventing excessive path expansion.

[0059] Smoothing term penalty cost The settings are based on the difference in passage resistance and spatial distance between adjacent grid cells. and When the labels are the same, It equals 0. When and When the labels are different, The value is inversely proportional to the difference in traffic resistance and inversely proportional to the spatial distance, which means that the penalty cost of assigning different labels to adjacent grid cells with greater resistance differences or greater spatial distance is lower, and the path boundary is more likely to cut along the location with significant resistance differences or greater distance.

[0060] Force the graph node corresponding to the starting point to be connected to the source point, and set the capacity of the connecting edge to infinity. Force the graph node corresponding to the target point to be connected to the sink point, and set the capacity of the connecting edge to infinity. The source point is the super source defined in the graph cut model, and the sink point is the super sink defined in the graph cut model.

[0061] Find the minimum cut in the graph cut model. Initialize the residual network by setting the current residual capacity of all directed edges in the graph cut model to the original capacity of the corresponding directed edges. The original capacity of the directed edges is set by connecting the graph nodes. Graph Nodes An undirected edge corresponds to two directed edges in opposite directions in the residual network, and the original capacity of each directed edge is equal to... The corresponding values ​​under different label combinations; the original capacity of the directed edge connecting the source node and the graph node is set according to the data item cost; the original capacity of the directed edge connecting the graph node and the sink node is set according to the data item cost. Set the source node as the currently active node and set the sink node as the terminated node.

[0062] A breadth-first search algorithm is used to search for all reachable nodes in the residual network starting from the source node. All visited nodes are marked, and each visited node is assigned a distance value representing the number of directed edges traversed from the source node to that node. A hierarchical graph is constructed. Only nodes with a distance value of [missing value] are retained in the hierarchical graph. The node points to a distance value of The directed edges of the nodes.

[0063] Starting from the source node, find an augmenting path from the source node to the sink node along the hierarchical graph. Each directed edge on the augmenting path satisfies the condition of increasing distance value. For each directed edge on the augmenting path, record the current remaining capacity of the directed edge, and determine the bottleneck capacity of the augmenting path. The bottleneck capacity is the minimum current remaining capacity among all directed edges on the augmenting path.

[0064] Subtract the bottleneck capacity from the current remaining capacity of each directed edge along the augmenting path, and add the bottleneck capacity to the current remaining capacity of the reverse edge of each directed edge, and update the residual network.

[0065] Repeat the steps of breadth-first search to build a hierarchical graph, find augmenting paths, and update the residual network until it is impossible to reach the sink from the source in the residual network. At this point, all nodes reachable from the source form the source set, and the remaining nodes form the sink set. The set of edges between the source set and the sink set is the minimum cut.

[0066] Extract all graph nodes marked with the first label from the graph node labels corresponding to the minimum cut. From these first-labeled nodes, select those belonging to the same connected component as the starting node. The connected component determination is based on the eight-neighbor spatial adjacency relationship between the graph nodes. Connect the raster cells corresponding to the selected nodes according to their spatial adjacency relationships to form a path, and output this path as the initial travel path.

[0067] In transforming the path planning problem into a graph cut problem using a graph cut algorithm, a multi-scale iterative optimization strategy is employed for graph cut solution. Gaussian pyramid downsampling is performed on the original resolution traffic resistance raster map to generate multiple traffic resistance raster maps at different scales. The Gaussian pyramid has three layers: the bottom layer is the original resolution layer, the middle layers are the layers after the first downsampling, and the top layer is the layer after the second downsampling. The resolution of each layer at a higher scale is half the square of the resolution of the layer at the next lower scale. During downsampling, bilinear interpolation is used to calculate the traffic resistance value of each downsampled raster cell; that is, the traffic resistance value of the downsampled raster cell is obtained by the weighted average of the traffic resistance values ​​of four adjacent raster cells within the corresponding spatial range in the original layer.

[0068] The graph cut solution is performed starting from the coarsest-scale layer, which is the top layer of the pyramid. After constructing the graph cut energy function and solving for the minimum cut on the coarsest-scale layer, the set of first-labeled nodes corresponding to the minimum cut obtained on the coarsest-scale layer is used as the initial path region under the coarsest-scale layer.

[0069] The initial path region at the coarse scale is upsampled and mapped to the adjacent finer-scale layer using bilinear interpolation. Bilinear interpolation upsampling assigns the label information of one raster cell in the coarse scale to the corresponding four subdivision raster cells in the finer-scale layer. The mapped region is then used as the initial label prior information for graph cuts in the finer-scale layer, and graph cut solutions with prior constraints are performed on the finer-scale layer. The prior constraints are reflected by modifying the data item cost; for graph nodes marked by the prior information as favoring the first label, the data item cost is adjusted. The original value is multiplied by an adjustment factor less than 1, with a value of 0.5. This adjustment factor is set to maintain the original labeling tendency of nodes in the prior region. For graph nodes that are marked by prior information as tending towards the second label, the data item cost is... Multiply the original value by an adjustment factor of 0.5.

[0070] The process is refined layer by layer, from the coarsest scale layer to the intermediate scale layer, and then to the original resolution layer. Each layer performs graph cut solutions under the prior constraints provided by the previous layer. After completing the graph cut solution on the original resolution layer, the set of the first labeled nodes in the graph cut results obtained on the original resolution layer is taken as the final graph cut solution, and the initial travel path is extracted from the final graph cut solution.

[0071] In specific implementation, please refer to Figure 4 The initial travel path is morphologically dilated to generate a path buffer. The position coordinates of all raster cells occupied by the initial travel path are obtained, and the set of all position coordinates is used as the initial path point set. Based on the digital surface model, the terrain relief of the local area where each path point in the initial path point set is located is calculated. The terrain relief is defined as the difference between the highest and lowest elevation values ​​within a pre-defined window centered on that path point. The pre-defined window size is set to 5 raster cells multiplied by 5 raster cells, with the center of the window being the raster cell containing the currently calculated path point, and the total number of raster cells covered by the window is 25. The elevation values ​​of all raster cells within the window are traversed, and the highest and lowest elevation values ​​are recorded. The difference between the highest and lowest elevation values ​​is used as the terrain relief of that path point.

[0072] For each pathpoint in the initial pathpoint set, the radius of the corresponding expansion structural element is dynamically determined based on the terrain relief of that pathpoint. There is a positive correlation between terrain relief and expansion radius; the greater the terrain relief, the larger the corresponding expansion radius. The minimum expansion radius is set to 1 grid cell, and the maximum expansion radius is set to 10 grid cells. The positive correlation between terrain relief and expansion radius is defined using a piecewise linear function, the specific form of which is:

[0073]

[0074] in, This represents the expansion radius, expressed in units of grid cells. The terrain relief of the waypoints is expressed in meters. Indicates the minimum expansion radius. The value is 1 grid cell. Indicates the maximum expansion radius. The value is 10 grid cells. Indicates the first threshold. The value is 10 meters. This represents the second threshold. The value is 50 meters. First threshold. The setting of 10 meters is based on the fact that when the terrain undulation is less than 10 meters, the terrain is relatively flat, and the path does not require a large buffer on both sides to meet the local adjustment needs of the path. Only a minimum expansion radius of 1 grid cell is needed. Second threshold The value of 50 meters is chosen because when the terrain undulation exceeds 50 meters, the terrain is drastically uneven, and the potential feasible areas on both sides of the path are widely distributed, requiring buffer coverage with a maximum expansion radius of 10 grid cells. When the terrain undulation... In the range of 10 meters to 50 meters, the expansion radius With the undulation of the terrain As the expansion radius increases linearly, the change in expansion radius is in a fixed proportion to the change in topographic relief.

[0075] A circular dilatation structuring element is used to perform a morphological dilatation operation on each path point in the initial path point set. The radius of the circular dilatation structuring element is the dilatation radius determined by the terrain relief. A circular expansion structural element is defined as a circle with a center and a radius of [missing information]. The path buffer is a collection of all raster cells. The morphological dilation operation places the center of the circular dilation structuring element at the center of the raster cell containing the current path point, marking all raster cells covered by the circular dilation structuring element as the path buffer region. The union of the raster cells covered by the morphological dilation operation on all path points is then used to generate a binary mask for the path buffer. In the binary mask, raster cells with a value of 1 (the first value) belong to the path buffer, while raster cells with a value of 0 (the second value) do not belong to the path buffer.

[0076] Within the buffer, local smoothing is performed using a digital surface model to obtain the final travel path. All raster cells with the first value are extracted from the binary mask of the path buffer. The elevation values ​​of the corresponding positions of the extracted raster cells on the digital surface model are used to construct the buffer elevation point cloud. Each point in the buffer elevation point cloud has three-dimensional spatial coordinates, which include the row and column position of the raster cell and its elevation value on the digital surface model.

[0077] For each pathpoint on the initial path, the normal direction of that pathpoint on the digital surface model is calculated. The normal direction is calculated as follows: a three-row, three-column local window is extracted from the digital surface model centered on the current pathpoint. The elevation values ​​of nine grid cells within this local window are used to determine a local fitting plane through least-squares plane fitting. The projection direction of the unit normal vector of this local fitting plane onto the horizontal plane is the normal direction of that pathpoint. A one-dimensional search is performed along the normal direction within the buffer elevation point cloud. The search range is limited to a search interval extending from the current pathpoint along both sides of the normal direction by a preset step size. The preset step size is set to 5 grid cells, based on the consideration that the maximum allowable offset during local path smoothing adjustments, without significantly altering the overall path orientation, is 5 grid cells, which can cover the optimal location of the local terrain. Within the search interval, all grid cells in the buffer elevation point cloud are traversed, and the elevation values ​​of each grid cell are compared one by one. The grid cell with the lowest elevation value in the normal direction is selected as the candidate replacement point for the current pathpoint.

[0078] Connect all candidate replacement points corresponding to all path points in sequence according to the order of the path points on the initial travel path. The connection method is to connect adjacent candidate replacement points with straight line segments to generate a preliminary smooth path.

[0079] A cubic spline interpolation algorithm is used to refine the initial smoothed path. The algorithm uses the spatial coordinates of all candidate replacement points as control nodes, constructing a cubic polynomial curve segment between every two adjacent control nodes. This segment satisfies the conditions of continuous function values, continuous first derivatives, and continuous second derivatives at the control nodes. Natural boundary conditions are specified for the entire curve, meaning the second derivatives at both ends are zero. The interpolation refinement step size is set to 0.5 grid cells, generating a continuous curve path passing through each candidate replacement point. This continuous curve path is then converted back to a grid cell sequence by traversing each interpolation point on the path, rounding the interpolation point coordinates to the nearest grid cell center, and sequentially removing duplicate grid cells to obtain the final travel path.

[0080] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A multi-terrain path intelligent planning method integrating high-resolution satellite remote sensing data, characterized in that, include: Acquire high-resolution satellite remote sensing data of the target mountainous area and extract digital surface models and orthophotos from it; Perform terrain analysis on the digital surface model to generate slope maps, aspect maps, and curvature maps; Identify land cover types based on orthophotos; A multi-topographic resistance raster map is constructed by integrating slope map, aspect map, curvature map and land cover type; The path planning problem is transformed into a graph cut problem using a graph cut algorithm. A starting point and a target point are set on a traffic resistance grid map. A graph cut energy function is constructed at the cost of traffic resistance. The path with the minimum cumulative resistance from the starting point to the target point is obtained as the initial traffic path. The initial travel path is morphologically dilated to generate a path buffer. Within the buffer, a digital surface model is used for local smoothing to obtain the final travel path.

2. The intelligent multi-terrain path planning method fusing high-resolution satellite remote sensing data according to claim 1, characterized in that, The steps of performing terrain analysis on the digital surface model to generate slope maps, aspect maps, and curvature maps specifically include: Traverse each grid cell in the digital surface model, extract a 3x3 grid neighborhood window centered on that grid cell, and record the elevation values ​​of the nine grid cells within the neighborhood window. The east-west and north-south elevation change rates of the central grid cells are calculated using a third-order finite difference operator. The east-west elevation change rate is obtained by dividing the elevation difference between the left and right sides of the central column by twice the grid spacing, and the north-south elevation change rate is obtained by dividing the elevation difference between the top and bottom sides of the central row by twice the grid spacing. The slope value of the central grid cell is calculated based on the east-west elevation change rate and the north-south elevation change rate. This slope value is the angle value corresponding to the square root of the sum of the squares of the east-west elevation change rate and the squares of the north-south elevation change rate. The calculated slope value is then filled into the corresponding position on the slope map. The slope aspect value of the central grid cell is calculated based on the east-west elevation change rate and the north-south elevation change rate. This slope aspect value is the arctangent angle of the ratio of the east-west elevation change rate to the north-south elevation change rate. The calculated slope aspect value is then filled into the corresponding position on the slope aspect diagram. A second difference operation is performed on the slope map, and the slope change rate of each grid cell is calculated using the same third-order finite difference operator as used to calculate the slope map, thus obtaining the curvature map.

3. The intelligent multi-terrain path planning method fusion method based on high-resolution satellite remote sensing data according to claim 1, characterized in that, The step of classifying and identifying land cover types based on orthophotos specifically includes: Extract the spectral feature vector of each pixel from the orthophoto image. The spectral feature vector includes four dimensions: blue band reflectance, green band reflectance, red band reflectance, and near-infrared band reflectance. Gray-level co-occurrence matrix analysis is performed on the local neighborhood window of each pixel, and texture feature values ​​are extracted from the gray-level co-occurrence matrix. These texture feature values ​​include contrast, correlation, and homogeneity. The spectral feature vector and texture feature value are concatenated to form the joint feature input vector of the pixel; The joint feature input vector is fed into a pre-trained support vector machine classifier, which uses a radial basis kernel function. Its training process is based on remote sensing data of sample areas with labeled land cover types. The support vector machine classifier outputs a probability vector for each pixel to belong to a predefined land cover type. The land cover type with the highest probability is selected as the recognition result for that pixel. The predefined land cover types include coniferous forest, broad-leaved forest, alpine meadow, bare rock, glacier, snow cover, and water body.

4. The intelligent multi-terrain path planning method fusing high-resolution satellite remote sensing data according to claim 1, characterized in that, The step of fusing slope maps, aspect maps, curvature maps, and land cover types to construct a multi-topographic resistance raster map specifically includes: The initial passage resistance value of each grid cell in the passage resistance raster map is calculated. This initial passage resistance value is composed of four parts: slope value, aspect deviation value, curvature correction value, and land cover type baseline resistance value. For grid cells whose slope value exceeds the preset upper limit of slope, set their initial passage resistance value to infinity; For grid cells whose slope values ​​do not exceed the upper limit of the slope, the angle between the slope direction of the grid cell and the pre-statistical prevailing wind direction of the area is calculated according to the aspect diagram. When the angle exceeds the pre-set angle threshold, an additional resistance term is added to the initial passage resistance value. Based on the curvature diagram, identify concave grid cells with negative curvature, and multiply the initial passage resistance value of such grid cells by a concave resistance amplification factor greater than one. Based on the classification results of land cover types, a corresponding baseline resistance value is assigned to each type of land cover, and this baseline resistance value is used as the basic component of the initial passage resistance value. The baseline resistance value for water cover types is set to infinity.

5. The multi-terrain path intelligent planning method fusion of high-resolution satellite remote sensing data according to claim 4, characterized in that, The angle threshold between the slope aspect and the pre-statistical prevailing wind direction in the region is set to a configurable parameter between 45 degrees and 90 degrees.

6. The intelligent multi-terrain path planning method fusing high-resolution satellite remote sensing data according to claim 1, characterized in that, The step of transforming the path planning problem into a graph cut problem using a graph cut algorithm, setting a starting point and a target point on a traffic resistance grid map, constructing a graph cut energy function at the cost of traffic resistance, and solving for the minimum cumulative resistance path from the starting point to the target point as the initial traffic path, specifically includes: Each grid cell in the traffic resistance grid diagram is treated as a graph node in the graph cut model, and the spatial adjacency between adjacent grid cells is treated as an edge in the graph cut model. Set a binary label variable for each graph node. The value of the binary label variable can be either a first label or a second label. The first label indicates that the grid cell belongs to the passageway area, and the second label indicates that the grid cell does not belong to the passageway area. Construct a graph cut energy function, which includes a data term and a smoothing term. The data term is set according to the passage resistance value of each grid cell. The smaller the passage resistance value, the lower the cost of assigning the first label to the grid cell. The smoothing term is set according to the passage resistance difference and spatial distance between adjacent grid cells. The larger the resistance difference or the larger the spatial distance between adjacent grid cells, the higher the penalty cost of assigning different labels to adjacent grid cells. Force the graph node corresponding to the starting point to be connected to the source point, and set the capacity of the edge connecting the source point to this graph node to infinity. Force the graph node corresponding to the target point to be connected to the sink point, and set the capacity of the edge connecting this graph node to the sink point to infinity. To find the minimum cut in the graph cut model, the path formed by all grid cells marked with the first label and connected to the starting point and the target point in the graph node label division corresponding to the minimum cut is taken as the initial travel path.

7. The intelligent multi-terrain path planning method fusing high-resolution satellite remote sensing data according to claim 1, characterized in that, The step of generating a path buffer by morphological dilation of the initial travel path specifically includes: Obtain the position coordinates of all grid cells occupied by the initial travel path to form the initial path point set; The terrain relief of the local area where each point in the initial path point set is located is calculated based on the digital surface model. The terrain relief is defined as the difference between the highest and lowest elevation values ​​within a pre-defined window centered on that point. For each initial path point, the radius of the corresponding expansion structure element is dynamically determined based on its terrain undulation. The terrain undulation is positively correlated with the expansion radius. The minimum value of the expansion radius is set to one grid cell, and the maximum value is set to ten grid cells. A circular dilation structuring element is used to perform a morphological dilation operation on each point in the initial path point set. The grid cells covered by all dilation operations are then combined to generate a binary mask for the path buffer. In this mask, grid cells with the first value belong to the path buffer, while grid cells with the second value do not belong to the path buffer.

8. The intelligent multi-terrain path planning method fusion method based on high-resolution satellite remote sensing data according to claim 7, characterized in that, The positive correlation between the terrain undulation and the expansion radius is defined by a piecewise linear function. When the terrain undulation is less than the first threshold, the expansion radius takes the minimum value. When it is greater than the second threshold, the expansion radius takes the maximum value. When it is in between, the expansion radius increases linearly with the terrain undulation.

9. The intelligent multi-terrain path planning method fusion method based on high-resolution satellite remote sensing data according to claim 7, characterized in that, The step of performing local smoothing processing within the buffer using a digital surface model to obtain the final travel path specifically includes: Extract all raster cells belonging to the path buffer from the binary mask of the path buffer to form the buffer elevation point cloud; For each path point on the initial travel path, calculate the normal direction of the path point on the digital surface model, and perform a one-dimensional search in the buffer elevation point cloud along the normal direction to find the raster cell with the lowest elevation value in the normal direction as the candidate replacement point of the path point. Connect all candidate replacement points of all path points in the order of the initial travel path to generate a preliminary smooth path; The initial smooth path is interpolated and encrypted using a cubic spline interpolation algorithm to generate a continuous curved path that passes through each candidate replacement point. This continuous curved path is then converted back into a raster cell sequence to obtain the final travel path.

10. The intelligent multi-terrain path planning method fusing high-resolution satellite remote sensing data according to claim 9, characterized in that, When the one-dimensional search finds the grid cell with the lowest elevation value in the normal direction, the search range is limited to a search interval formed by extending a preset number of steps on both sides of the normal direction with the current path point as the center.

Images