Urban engineering planning method based on unmanned aerial vehicle remote sensing surveying and mapping

Abstract

The invention discloses an urban engineering planning method based on unmanned aerial vehicle remote sensing surveying and mapping, which relates to the technical field of unmanned aerial vehicle remote sensing surveying and mapping, and comprises the steps of mapping a gradient magnitude into an exponential growth difficulty value, forming a terrain traffic difficulty field, calculating an Euclidean distance from each grid point to a nearest obstacle grid based on an obstacle binary image, and calculating a reciprocal. The terrain traffic difficulty field is defined as an obstacle exclusion potential field, a steady-state diffusion equation is constructed in combination with the terrain traffic difficulty field and the obstacle exclusion potential field, the steady-state diffusion equation is solved by adopting a Gaussian-Seidel iteration method, and a dynamic cost surface is obtained; and generating an optimal path sequence, and constructing a three-dimensional engineering entity model based on the optimal path sequence. According to the method, through combination of NDVI vegetation identification and DSM morphological operation, the extraction accuracy of urban ground object obstacles is enhanced, and through combination of a terrain difficulty field mapped by a slope index and an obstacle exclusion potential field, the continuity and physical rationality of cost field construction are improved.

Description

Technical Field

[0001] This invention relates to the field of unmanned aerial vehicle (UAV) remote sensing mapping technology, and in particular to an urban engineering planning method based on UAV remote sensing mapping. Background Technology

[0002] With the rapid advancement of urbanization, urban engineering planning has placed higher demands on the accuracy, real-time nature, and multi-dimensional integration capabilities of spatial information. Traditional urban planning surveying and mapping relies mainly on ground surveying, single-view aerial photogrammetry, or satellite remote sensing, which have inherent limitations in terms of spatial resolution, the impact of ground cover occlusion, and dynamic updating capabilities.

[0003] The existing technical approach to combining UAV mapping with urban planning still has many shortcomings. The preprocessing and fusion of multi-source remote sensing data still rely on independent processes. LiDAR point clouds, oblique photogrammetric images, and multispectral images often need to be processed step by step. The lack of a unified spatial alignment framework leads to limited DSM construction quality and insufficient discrimination accuracy for elevation-sensitive areas such as buildings and vegetation. Existing path planning technologies mostly use discrete graph search (such as A*, Dijkstra) or cost field methods based on empirical weights, lacking continuous mathematical modeling of terrain gradients, obstacle potential fields, and multi-source ground feature characteristics, making it difficult for the generated paths to reflect engineering feasibility. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides an urban engineering planning method based on UAV remote sensing mapping, which solves the problem that the preprocessing and fusion of multi-source remote sensing data still rely on independent processes. LiDAR point clouds, oblique photogrammetric images, and multispectral images often need to be processed step by step, and there is a lack of a unified spatial alignment framework, which limits the quality of DSM construction and is insufficient in the discrimination accuracy of elevation-sensitive areas such as buildings and vegetation. Existing path planning technologies mostly use discrete graph search or cost field methods based on empirical weights, which lack continuous mathematical modeling of terrain gradients, obstacle potential fields, and multi-source ground feature features, making it difficult for the generated path to reflect the feasibility of engineering.

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

[0007] In a first aspect, the present invention provides an urban engineering planning method based on unmanned aerial vehicle (UAV) remote sensing mapping, which includes,

[0008] Collect and preprocess multi-source data, including multispectral images, oblique photogrammetric images, and lidar point cloud data;

[0009] Generate image point clouds, spatially register the lidar point clouds with the image point clouds, and merge them to obtain fused point clouds. Perform regular meshing on the fused point clouds to generate digital surface models. Perform binarization processing on multispectral images to obtain vegetation binary maps. Perform morphological opening operations on the digital surface models to extract elevation values ​​from the smoothed surface models and generate building binary maps. Perform a logical "OR" operation on the vegetation binary maps and building binary maps to generate obstacle binary maps.

[0010] Based on the digital surface model, the gradient of each grid point is calculated, and the gradient magnitude is mapped to an exponentially increasing difficulty value to form a terrain access difficulty field. Based on the obstacle binary map, the Euclidean distance from each grid point to the nearest obstacle grid is calculated, and the reciprocal is calculated and defined as the obstacle repulsion potential field. Combining the terrain access difficulty field and the obstacle repulsion potential field, a steady-state diffusion equation is constructed. The Gauss-Seidel iterative method is used to solve the steady-state diffusion equation to obtain the dynamic cost surface.

[0011] The dynamic cost surface is converted into a comprehensive travel velocity field of each point in space. The equation is constructed and solved using the fast travel method to obtain the minimum arrival time. The negative gradient direction of the minimum arrival time field is calculated and converted into the next path point to generate the optimal path sequence. A three-dimensional engineering entity model is constructed based on the optimal path sequence, and the data of the three-dimensional engineering entity model is packaged.

[0012] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the steps of binarizing multispectral images to obtain a vegetation binary map, performing morphological opening operations on a digital surface model to extract elevation values ​​from the smoothed surface model, generating a building binary map, and performing a logical OR operation between the vegetation binary map and the building binary map to generate an obstacle binary map include:

[0013] Jointly perform motion reconstruction and multi-view stereo matching on multispectral and oblique photogrammetric images to generate image point clouds;

[0014] Spatial registration is performed between the lidar point cloud and the image point cloud. The registered lidar point cloud and the image point cloud are then merged to obtain a fused point cloud.

[0015] The fused point cloud is subjected to rule-based meshing, and the elevation value of each grid cell is calculated using the inverse distance weighted interpolation method to generate a digital surface model.

[0016] Normalized vegetation index (NVI) maps are calculated based on multispectral images. Regional thresholds are set, and the NVI maps are binarized to obtain binary vegetation maps.

[0017] Morphological opening operations are performed on the digital surface model to obtain a smoothed surface model. Elevation values ​​are extracted from the smoothed surface model. Fixed windows are set, and the standard deviation of elevation within each fixed window and the difference between the center elevation and the lowest elevation of the window are calculated with each grid as the center. Height threshold and standard deviation threshold are set to generate a binary map of the building.

[0018] A binary image of obstacles is generated by performing a logical OR operation between the binary image of vegetation and the binary image of buildings.

[0019] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the following steps are included: Based on a digital surface model, the gradient of each grid point is calculated, and the gradient magnitude is mapped to an exponentially increasing difficulty value to form a terrain accessibility difficulty field. Based on a binary obstacle map, the Euclidean distance from each grid point to the nearest obstacle grid is calculated, and its reciprocal is calculated and defined as the obstacle repulsion potential field, including:

[0020] Based on the digital surface model, the Sobel operator is used to calculate the gradient of each grid point on x and y, the magnitude of the gradient is defined as the slope, and the slope is mapped to an exponentially increasing difficulty value to form a terrain accessibility field.

[0021] Based on the binary obstacle map, calculate the Euclidean distance from each grid point to the nearest obstacle grid, and calculate the reciprocal, which is defined as the obstacle repulsive potential field.

[0022] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the following steps are taken: A steady-state diffusion equation is constructed by combining the terrain accessibility field and the obstacle repulsion potential field; the Gauss-Seidel iterative method is used to solve the steady-state diffusion equation to obtain the dynamic cost surface, including:

[0023] By combining the terrain accessibility field and the obstacle repulsion potential field, a steady-state diffusion equation is constructed and solved using the Gauss-Seidel iterative method. A stopping threshold is set using the residual criterion. The iteration stops when the L2 norm of the residual vector of two consecutive iterations is less than the stopping threshold, thus obtaining the dynamic cost surface.

[0024] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the following steps are included: converting the dynamic cost surface into a comprehensive travel velocity field for each point in space, constructing a process function equation, solving the process function equation using the fast travel method to obtain the minimum arrival time, calculating the negative gradient direction of the minimum arrival time field, converting it into the next path point, and generating an optimal path sequence.

[0025] The dynamic cost surface is converted into a comprehensive travel velocity field at various points in space. An equation is constructed, and the equation is solved using the fast travel method to obtain the minimum arrival time.

[0026] Initialize the path sequence, calculate the negative gradient direction of the minimum arrival time field, and convert it into the next path point;

[0027] When the Euclidean distance between a newly generated path point and the initial path point is less than one backtracking step, the initial path point is added to the sequence as the last point, the backtracking terminates, and the path point coordinate sequence is obtained. The sequence is then reversed and defined as the optimal path sequence.

[0028] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the step of constructing a three-dimensional engineering entity model based on the optimal path sequence includes:

[0029] Connect the optimal path sequence to form the initial centerline geometric entity, and perform cubic spline smoothing to obtain the final path centerline;

[0030] Using the path centerline as the axis, perform Euclidean buffer analysis to generate a fixed-distance buffer zone, which is defined as the engineering land area.

[0031] The binary map of the obstacle is spatially overlaid with the area of ​​the engineering land area to calculate the area encroached by the obstacle.

[0032] Design station points are set up along the path centerline at fixed pile intervals. The elevation values ​​of the design station points are extracted from the digital surface model, and Kriging interpolation is used to generate a continuous design elevation surface.

[0033] The two-dimensional boundary polygon of the engineering site is stretched vertically to the design elevation surface to generate a closed three-dimensional engineering solid model.

[0034] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the step of packaging the three-dimensional engineering entity model into data includes:

[0035] The 3D engineering solid model is triangulated, and the output file is automatically named according to the rule of "Project ID_Deliverable Type_Timestamp".

[0036] As a preferred embodiment of the urban engineering planning method based on UAV remote sensing mapping described in this invention, the step of collecting multi-source data and performing preprocessing includes:

[0037] Once the planning task is triggered, drones are dispatched to survey the target planning area, collect multi-source data of the city, and perform noise reduction and standardization processing.

[0038] In a second aspect, the present invention provides a computer device including a memory and a processor, wherein the memory stores a computer program, wherein when the computer program is executed by the processor, it implements any step of the urban engineering planning method based on UAV remote sensing mapping as described in the first aspect of the present invention.

[0039] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein: when the computer program is executed by a processor, it implements any step of the urban engineering planning method based on UAV remote sensing mapping as described in the first aspect of the present invention.

[0040] The beneficial effects of this invention are as follows: By combining NDVI vegetation recognition with DSM morphological operations, this invention enhances the accuracy of urban terrain obstacle extraction; by combining the terrain difficulty field mapped by the slope index with the obstacle repulsion potential field, it improves the continuity and physical rationality of the cost field construction; and by combining the solution of the steady-state diffusion equation with the fast travel method, it improves the convergence stability and geometric smoothness of cost propagation and optimal path calculation, making the planned path significantly superior to existing technologies in terms of engineering feasibility. Attached Figure Description

[0041] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0042] Figure 1 This is a flowchart of the urban engineering planning method based on UAV remote sensing mapping in Example 1.

[0043] Figure 2 This is a schematic diagram of the urban engineering planning method based on UAV remote sensing mapping in Example 1. Detailed Implementation

[0044] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0045] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0046] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0047] Example 1, referring to Figure 1 and Figure 2 This is the first embodiment of the present invention, which provides an urban engineering planning method based on UAV remote sensing mapping, including the following steps:

[0048] S1. Collect and preprocess multi-source data, including multispectral images, oblique photogrammetric images, and lidar point cloud data;

[0049] Specifically, this involves collecting and preprocessing multi-source data, including:

[0050] Once the planning task is triggered, a drone equipped with a lidar (such as DJI L1), a multispectral camera (such as RedEdge-MX), and a five-lens oblique photography camera (such as SW-5) is dispatched to survey the target planning area (based on the planning area setting), collect multi-source data of the city, and perform noise reduction and standardization processing.

[0051] The joint collection of multi-source data significantly improves the spatiotemporal consistency of different data sources. The preprocessed data has the ability to be transferred across scales and regions and can be directly applied to different urban engineering planning scenarios. Data standardization significantly reduces the algorithm's dependence on human intervention, which is conducive to the formation of an end-to-end UAV surveying-modeling-planning technical system.

[0052] S2. Generate image point cloud, spatially register the lidar point cloud with the image point cloud, and merge them to obtain a fused point cloud. Perform regular meshing on the fused point cloud to generate a digital surface model. Perform binarization processing on the multispectral image to obtain a vegetation binary map. Perform morphological opening operation on the digital surface model to extract the elevation values ​​in the smoothed surface model and generate a building binary map. Perform a logical "OR" operation on the vegetation binary map and the building binary map to generate an obstacle binary map.

[0053] Specifically, the multispectral image is binarized to obtain a vegetation binary map; morphological opening is performed on the digital surface model to extract elevation values ​​from the smoothed surface model, generating a building binary map; and a logical OR operation is performed between the vegetation binary map and the building binary map to generate an obstacle binary map, including:

[0054] The multispectral and oblique photogrammetric images are jointly processed to perform Structure from Motion (SfM) and Multi-View Stereo (MVS) matching to generate image point clouds.

[0055] SfM uses feature matching between images and camera pose optimization to reconstruct sparse point clouds. MVS generates dense point clouds based on this, and its three-dimensional coordinates are strictly aligned with the original image POS data (the camera's position (latitude, longitude, altitude, and other geographic coordinate information) and attitude (such as heading angle, pitch angle, and roll angle) when capturing remote sensing images, usually provided by GNSS (Global Navigation Satellite System) and IMU (Inertial Measurement Unit).

[0056] Spatial registration is performed between the lidar point cloud and the image point cloud. The geographic coordinates of the two are used as the initial alignment basis. The Iterative Closest Point (ICP) algorithm is used to solve the optimal rigid transformation, so as to minimize the distance error between the two point clouds. The registered lidar point cloud and the image point cloud are then merged to obtain the fused point cloud.

[0057] The fused point cloud is subjected to regular meshing, and the elevation value of each mesh cell (0.5 m × 0.5 m) is calculated using the Inverse Distance Weighting (IDW) method to generate a digital surface model. The formula is as follows:

[0058] ,

[0059] in For digital surface models, For grid points, To merge the elevation values ​​of the k-th point in the point cloud, is the Euclidean distance from the grid point to the horizontal coordinate of the k-th point cloud point, and n is the number of nearest neighbors;

[0060] The normalized vegetation index map is calculated based on multispectral imagery, using the following formula:

[0061] ,

[0062] in This is a normalized vegetation index (NDI) map, representing the NDI at grid point (x, y). Let (x, y) be the reflectance value of the grid point (x, y) in the near-infrared band. The reflectance value of grid point (x,y) in the red light band;

[0063] By setting a regional threshold, the normalized vegetation index map is binarized to obtain a binary vegetation map. The formula is as follows:

[0064] ,

[0065] in This is a binary image of vegetation. The Otsu threshold method is used to set the region threshold.

[0066] Morphological opening operations (erosion followed by dilation) are performed on the digital surface model to filter out small noises (such as cars and small shrubs) and obtain a smoothed surface model.

[0067] Elevation values ​​are extracted from the smoothed surface model. An empirical window method is used to set fixed windows (e.g., 5×5). Centered on each grid, the standard deviation of elevation within each fixed window and the difference between the window's center elevation and its lowest elevation are calculated. A height threshold (typical minimum building height) and a standard deviation threshold (to distinguish between flat roofs and natural slopes) are set to generate a binary image of the building. The formula is as follows:

[0068] ,

[0069] in A binary image of the building. The elevation difference between grid points (x, y) The altitude threshold is set based on the "Classification Standard for Low-Altitude Flight Obstacles". Let (x, y) be the standard deviation of the elevation of the grid point. The standard deviation threshold (calibrated using the statistical method of surface roughness of typical ground features to ensure that flat roofs are preserved and natural slopes are excluded) is recommended for low-altitude obstacle avoidance modeling in the "Technical Guidelines for Obstacle Extraction in UAV Photogrammetry" (Trial, Ministry of Natural Resources, 2023).

[0070] Perform a logical OR operation between the vegetation binary map and the building binary map to generate an obstacle binary map. An obstacle binary map value of 1 indicates that there are insurmountable obstacles such as vegetation or buildings at that location, while an obstacle binary map value of 0 indicates that the area is passable.

[0071] LiDAR maintains penetration in vegetation and low-texture areas, while image point clouds perform better on building facades and areas with rich texture. The fusion of the two achieves a global point cloud with high density and high geometric quality. Compared with radial basis function or simple averaging, IDW is more sensitive to local elevation gradients. NDVI can effectively eliminate reflectivity deviations caused by changes in illumination and exposure, ensuring that vegetation areas still have stable discriminability under different environmental conditions. Elevation difference and local standard deviation constitute a dual-threshold judgment system, which can eliminate natural slopes and stepped terrain, and enhance the geometric consistency of building recognition. Obstacle binary maps provide basic semantic constraints for distance transformation, potential energy modeling, and path planning.

[0072] S3. Based on the digital surface model, calculate the gradient of each grid point, map the gradient magnitude to an exponentially increasing difficulty value to form a terrain access difficulty field. Based on the obstacle binary map, calculate the Euclidean distance from each grid point to the nearest obstacle grid and calculate the reciprocal, which is defined as the obstacle repulsion potential field. Combine the terrain access difficulty field and the obstacle repulsion potential field to construct a steady-state diffusion equation. Use the Gauss-Seidel iterative method to solve the steady-state diffusion equation to obtain the dynamic cost surface.

[0073] Specifically, based on the digital surface model, the gradient of each grid point is calculated, and the gradient magnitude is mapped to an exponentially increasing difficulty value to form a terrain traversal difficulty field. Based on the obstacle binary map, the Euclidean distance from each grid point to the nearest obstacle grid is calculated, and its reciprocal is calculated and defined as the obstacle repulsion potential field, including:

[0074] Based on the digital surface model, the Sobel operator is used to calculate the gradient of each grid point in the east-west x direction and the north-south y direction, i.e. the elevation change rate. The magnitude of the gradient in the two directions is defined as the slope.

[0075] Mapping slope to exponentially increasing difficulty values ​​creates a terrain traversal difficulty field, the formula is:

[0076] ,

[0077] in This represents the terrain traversal difficulty field, specifically the difficulty value at the grid point (x, y). Let (x, y) be the slope at grid point (x, y). To achieve the maximum allowable slope, the climb capability of typical multirotor aircraft is based on the "Civil Unmanned Aerial Vehicle Operation Performance Standards";

[0078] Based on the binary obstacle map, a linear-time Euclidean distance transformation algorithm (such as the Felzenszwalb algorithm) is used to calculate the Euclidean distance from each grid point to the nearest obstacle grid, and the reciprocal is calculated and defined as the obstacle repulsive potential field, with the formula as follows:

[0079] ,

[0080] ,

[0081] in Let (x, y) be the Euclidean distance from the grid point (x, y) to the nearest obstacle. Let B be the coordinates of the pixels with a value of 1 in the binary image of the obstacle, and let B be the binary image of the obstacle. To prevent infinite values ​​from appearing on obstacles, a small positive number needs to be added as a smoothing constant. This constant is set to half the grid resolution. Let be the obstacle repulsion potential at grid point (x,y).

[0082] Compared to traditional threshold-based piecewise functions, the exponential difficulty field exhibits stronger smoothness, avoiding abrupt changes in the cost field and providing support for the analytical and numerical stability of subsequent diffusion equations. The exponential function allows for a sharp increase in cost, enhancing the planning's ability to identify critical engineering structures such as steep slopes and building platform edges. In engineering design where smoothness is paramount, the exponential difficulty field naturally prevents planned paths from crossing steep slopes, thus reducing earthwork and improving construction economy. As grid points approach obstacles, the potential energy increases rapidly, naturally guiding the planned path away from buildings, vegetation, and inaccessible areas, achieving a continuous expression of obstacle avoidance capabilities. This invention's potential field expresses the "radius of influence" of obstacles on the surrounding area, preventing paths from hugging narrow passages and improving safety.

[0083] Furthermore, by combining the terrain accessibility field and the obstacle repulsion potential field, a steady-state diffusion equation is constructed. The Gauss-Seidel iterative method is used to solve the steady-state diffusion equation, yielding the dynamic cost surface, including:

[0084] Combining the terrain accessibility field and the obstacle repulsion potential field, a steady-state diffusion equation (elliptic partial differential equation) is constructed, as follows:

[0085] ,

[0086] ,

[0087] in and For gradient operator and divergence operator, The diffusion term represents the flow of cost R from the high-value range to the low-value range. For attenuation term, For source terms, Let be the value of the dynamic cost surface to be determined at the grid point (x,y). and The source strength coefficients of the terrain difficulty field and the obstacle repulsion potential field are set by the reciprocals of the mean values ​​of the terrain traversal difficulty field and the obstacle repulsion potential field, respectively. The attenuation coefficient is set using experimental calibration methods (such as test indicators).

[0088] On the boundary of the computational domain, a homogeneous Neumann boundary condition is applied, as shown in the formula:

[0089] ,

[0090] in For the planning area The boundary, Let be the unit outward normal vector at the point (x, y) on the boundary. The directional derivative of the cost field R at the boundary point along the outward normal direction;

[0091] The steady-state diffusion equation is discretized on a regular grid using the finite difference method, and then solved using the Gauss-Seidel iterative method, as shown in the following formula:

[0092] ,

[0093] ,

[0094] Where A is the coefficient matrix, R is the cost field vector to be determined, and b is the source term vector. Let be the approximate cost at grid point (i,j) after the k-th iteration, where i and j are the indices of the regular grid in the row and column directions, respectively. For the unknowns in a discrete linear system of equations The diagonal coefficient, For each term in the discrete linear system of equations, the term on the right-hand side (source term) corresponding to equation (i,j) is denoted as . It is the set of neighboring indices that are directly related to grid point (i,j) in the discrete equation (including (i±1,j) and (i,j±1)). For the off-diagonal coefficients corresponding to the neighboring point (p,q), To take the value of the neighboring points in the k-th iteration when summing (using Gauss-Seidel iteration, some neighboring points may have been updated to the value of the k+1th iteration);

[0095] Using the residual criterion to set a stopping threshold, the iteration stops when the L2 norm of the residual vector of two consecutive iterations is less than the stopping threshold, thus obtaining the dynamic cost surface.

[0096] Unlike traditional cost superposition, this invention generates a natural "cost density concentration" phenomenon through diffusion suppression, which is closer to the energy distribution on real terrain. The Gauss-Seidel method converges faster and has lower computational latency than the Jacobi method, making it suitable for real-time planning needs of large-scale grids (such as millions of nodes). The grid structure is naturally matched with finite difference, making it easy for engineers to understand. L2 norm monitoring ensures that the final static cost field reaches a stable state, avoiding premature stopping or excessive iteration. The cost concentration area can correspond to potentially unstable slopes, weathering troughs, and other geological structures that require special attention.

[0097] S4. Convert the dynamic cost surface into a comprehensive travel velocity field of each point in space, construct the equation, solve the equation using the fast travel method to obtain the minimum arrival time, calculate the negative gradient direction of the minimum arrival time field, convert it into the next path point, generate the optimal path sequence, construct a three-dimensional engineering entity model based on the optimal path sequence, and package the data of the three-dimensional engineering entity model.

[0098] Specifically, the dynamic cost surface is converted into a comprehensive travel velocity field at each point in space, an equation is constructed, and the equation is solved using the fast travel method to obtain the minimum arrival time. The negative gradient direction of the minimum arrival time field is calculated and converted into the next path point, generating an optimal path sequence, including:

[0099] The dynamic cost surface is converted into a comprehensive travel velocity field at various points in space, as shown in the formula:

[0100] ,

[0101] in The combined velocity field at grid point (x,y) The benchmark traffic speed is determined according to the design speed value of the corresponding road grade in the "Code for Design of Urban Road Engineering" (CJJ 37-2012);

[0102] A procedural equation is constructed, which satisfies the wavefront propagation process. This equation describes the physical relationship between the gradient magnitude and velocity at the wavefront arrival time, and the formula is as follows:

[0103] ,

[0104] in Let be the minimum arrival time from the starting point to the grid point (x, y). Let be the gradient magnitude of the time field. This equation is a classic equation in geometric optics describing wavefront propagation.

[0105] The Fast Marching Method (FMM) is used to solve the equation of the process function to obtain the minimum arrival time;

[0106] Mark all grid points as "far point" with their arrival time initialized to infinity. Mark the grid point where the starting point is located as "accepting point" with its arrival time set to 0. For each of the four neighboring grid points of the starting point, mark its state as "test point" and estimate its initial arrival time using the first-order upwind difference formula. Insert all "test points" into a min-heap priority queue according to their currently estimated initial arrival time. Pop the "test point" with the smallest current arrival time from the priority queue, denoted as Q, and change the state of point Q to "accepting point". For each of the four neighboring points P of Q, if P's state is "accepting point", skip it. For a neighboring point P, calculate a candidate value of U(P) based on the arrival times of its surrounding neighbors that are already "accepting points". Consider only the neighbors with the state of "accepting point", and construct a quadratic equation about U, the formula is:

[0107] ,

[0108] in This is the minimum arrival time of the neighbor whose state is 'accepting point' among point P's two direct neighbors (left and right grid points) in the x-direction. If neither neighbor is a 'accepting point', this item is not included in the calculation. Let P be the minimum arrival time of the neighbor whose state is 'accepting point' among the two direct neighbors (upper and lower grid points) of point P in the y direction. Let the combined travel speed be at the neighboring point p. Let U be the spacing of the regular grid, and U be the candidate arrival time of the neighboring point P to be determined.

[0109] Solving the quadratic equation yields the candidate arrival times. ;

[0110] If the calculated candidate arrival time is less than the currently stored U value of point P, then update the U value of P to the candidate arrival time. If the original state of point P was "far point", then change its state to "test point" and set it to... Insert into the priority queue. If the state of point P is already "test point", update its U value in the priority queue and adjust the heap structure until the state of the destination becomes "accept point" or the priority queue is empty, and finally obtain the minimum arrival time field.

[0111] Initialize the path sequence, calculate the negative gradient direction of the minimum arrival time field, and convert it to the next path point. The formula is:

[0112] ,

[0113] ,

[0114] in It is a unit direction vector. Let K be the two-dimensional spatial coordinates of the k-th path point. for, The backtracking step size is set to the grid resolution. Let be the unit vector pointing in the direction of the fastest time decrease at the k-th path point. and For the minimum arrival time field U at the path point The approximate values ​​of the partial derivatives along the x and y directions are calculated using the one-sided difference method;

[0115] When the Euclidean distance between a newly generated path point and the initial path point is less than one backtracking step, the initial path point is added to the sequence as the last point, the backtracking terminates, and the path point coordinate sequence is obtained. The sequence is then reversed and defined as the optimal path sequence.

[0116] The greater the cost and the lower the speed, the slower the propagation speed in unfavorable areas (slopes, dense obstacles), thus providing a physical basis for the concept of "minimum time" in the subsequent equation. The fast-moving method is logically similar to a "monotonically outward expanding wavefront". It only needs a single scan from the inside to the outside to obtain the minimum arrival time field of the entire domain, avoiding the high complexity of traditional Dijkstra's operation in irregular grids. By constructing a local time equation using only "accepted" neighbors through first-order upwind difference, the propagation direction always follows the wavefront's direction of movement, ensuring the correctness and physical rationality of the minimum time solution and avoiding pseudo-solutions of "backflow propagation". Unlike the shortest path based on grid neighborhood, gradient backtracking operates in continuous space. Even if the grid is coarse, the obtained path still appears as a smooth curve. Only a small amount of geometric smoothing is needed to meet the engineering requirements for curve radius and compliance.

[0117] Furthermore, a three-dimensional engineering entity model is constructed based on the optimal path sequence, including:

[0118] Connect the optimal path sequence to form the initial centerline geometry, and then perform cubic spline smoothing to obtain the final path centerline, as shown in the formula:

[0119] ,

[0120] Where c represents the initial centerline geometry, N is the number of path segments, and represents the index of the last point in the optimal path sequence. For connection and The straight line segment;

[0121] Using the path centerline as the axis, perform Euclidean buffer analysis to generate a fixed-distance buffer zone, defined as the engineering land area surface, using the following formula:

[0122] ,

[0123] in The area of ​​the project site is defined as W, which is the standard cross-sectional width of the project, determined according to the engineering design specifications, and o is a point on the center line of the path.

[0124] The obstacle binary map is spatially overlaid with the engineering land area to calculate the obstacle encroachment area. The formula is as follows:

[0125] ,

[0126] in The area occupied by the obstacle;

[0127] Design station points are laid out along the centerline of the path at fixed station intervals (e.g., 10 meters, which can be adjusted to 5–20 meters according to the accuracy requirements of the project). The elevation values ​​of the design station points are extracted from the digital surface model, and Kriging interpolation is used to generate a continuous design elevation surface.

[0128] The two-dimensional boundary polygon of the engineering land area is stretched along the vertical direction (Z-axis) to the design elevation surface to generate a closed three-dimensional engineering solid model.

[0129] The stretching operation is achieved through triangulation.

[0130] After the centerline is smoothed by three splines, the curvature changes continuously, which can meet the hard indicators of minimum curve radius and line smoothness in engineering projects such as roads, utility tunnels, and tunnels. The smooth line reduces sharp turns and angles of intersection, which is conducive to controlling superelevation and reducing abrupt changes in cross slope. At the same time, it also reduces stress concentration caused by geometric abrupt changes in structural stress. By setting the buffer radius to half the standard cross section width, the land area that meets the specifications for road red line width or utility tunnel width can be obtained. The amount of demolition and the degree of ecological impact can be predicted in the planning stage, avoiding the need to make major adjustments to the route when large-scale encroachment on sensitive areas is discovered in the later design stage.

[0131] Furthermore, the data of the 3D engineering solid model is packaged, including:

[0132] The path centerline and the engineering land area are saved in GeoJSON and ESRI Shapefile formats respectively. The area occupied by obstacles is stored in a structured form (such as CSV or JSON format) and associated with the vector data (path centerline and engineering land area). The 3D engineering entity model is triangulated. All output files are automatically named according to the rule of "Project ID_Result Type_Timestamp".

[0133] This embodiment also provides a computer device applicable to urban engineering planning methods based on UAV remote sensing mapping, including: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to realize the urban engineering planning method based on UAV remote sensing mapping proposed in the above embodiment.

[0134] The computer device can be a terminal, comprising a processor, memory, communication interface, display screen, and input devices connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, carrier networks, NFC (Near Field Communication), or other technologies. The display screen can be an LCD screen or an e-ink screen. The input devices can be a touch layer covering the display screen, buttons, a trackball, or a touchpad on the computer device's casing, or an external keyboard, touchpad, or mouse.

[0135] This embodiment also provides a storage medium storing a computer program, which, when executed by a processor, implements the urban engineering planning method based on UAV remote sensing mapping as proposed in the above embodiments. The storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read Only Memory (EPROM), Programmable Red-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0136] In summary, this invention enhances the accuracy of urban obstacle extraction by combining NDVI vegetation recognition with DSM morphological operations. It improves the continuity and physical rationality of the cost field construction by combining the terrain difficulty field mapped by the slope index with the obstacle repulsion potential field. Furthermore, it improves the convergence stability and geometric smoothness of cost propagation and optimal path calculation by solving the steady-state diffusion equation and combining it with the fast travel method. As a result, the planned path is significantly superior to existing technologies in terms of engineering feasibility.

[0137] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for urban engineering planning based on UAV remote sensing mapping, characterized in that: include, Collect and preprocess multi-source data, including multispectral images, oblique photogrammetric images, and lidar point cloud data; Generate image point clouds, spatially register the lidar point clouds with the image point clouds, and merge them to obtain fused point clouds. Perform regular meshing on the fused point clouds to generate digital surface models. Perform binarization processing on multispectral images to obtain vegetation binary maps. Perform morphological opening operations on the digital surface models to extract elevation values ​​from the smoothed surface models and generate building binary maps. Perform a logical "OR" operation between the vegetation binary maps and the building binary maps to generate obstacle binary maps. Based on the digital surface model, the gradient of each grid point is calculated, and the gradient magnitude is mapped to an exponentially increasing difficulty value to form a terrain access difficulty field. Based on the obstacle binary map, the Euclidean distance from each grid point to the nearest obstacle grid is calculated, and the reciprocal is calculated and defined as the obstacle repulsion potential field. Combining the terrain access difficulty field and the obstacle repulsion potential field, a steady-state diffusion equation is constructed. The Gauss-Seidel iterative method is used to solve the steady-state diffusion equation to obtain the dynamic cost surface. The dynamic cost surface is converted into a comprehensive travel velocity field of each point in space. The equation is constructed and solved using the fast travel method to obtain the minimum arrival time. The negative gradient direction of the minimum arrival time field is calculated and converted into the next path point to generate the optimal path sequence. A three-dimensional engineering entity model is constructed based on the optimal path sequence, and the data of the three-dimensional engineering entity model is packaged.

2. The urban engineering planning method based on UAV remote sensing mapping as described in claim 1, characterized in that: The process involves binarizing the multispectral image to obtain a vegetation binary map, performing morphological opening operations on the digital surface model to extract elevation values ​​from the smoothed surface model, generating a building binary map, and performing a logical OR operation between the vegetation binary map and the building binary map to generate an obstacle binary map. Jointly perform motion reconstruction and multi-view stereo matching on multispectral and oblique photogrammetric images to generate image point clouds; Spatial registration is performed between the lidar point cloud and the image point cloud. The registered lidar point cloud and the image point cloud are then merged to obtain a fused point cloud. The fused point cloud is subjected to rule-based meshing, and the elevation value of each grid cell is calculated using the inverse distance weighted interpolation method to generate a digital surface model. Normalized vegetation index (NVI) maps are calculated based on multispectral images. Regional thresholds are set, and the NVI maps are binarized to obtain binary vegetation maps. Morphological opening operations are performed on the digital surface model to obtain a smoothed surface model. Elevation values ​​are extracted from the smoothed surface model. Fixed windows are set, and the standard deviation of elevation within each fixed window and the difference between the center elevation and the lowest elevation of the window are calculated with each grid as the center. Height threshold and standard deviation threshold are set to generate a binary map of the building. A binary image of obstacles is generated by performing a logical OR operation between the binary image of vegetation and the binary image of buildings.

3. The urban engineering planning method based on UAV remote sensing mapping as described in claim 2, characterized in that: The process, based on a digital surface model, calculates the gradient at each grid point, maps the gradient magnitude to an exponentially increasing difficulty value, forming a terrain traversal difficulty field. Based on the obstacle binary map, it calculates the Euclidean distance from each grid point to the nearest obstacle grid and calculates its reciprocal, defining it as the obstacle repulsion potential field, including: Based on the digital surface model, the Sobel operator is used to calculate the gradient of each grid point on x and y, the magnitude of the gradient is defined as the slope, and the slope is mapped to an exponentially increasing difficulty value to form a terrain accessibility field. Based on the binary obstacle map, calculate the Euclidean distance from each grid point to the nearest obstacle grid, and calculate the reciprocal, which is defined as the obstacle repulsive potential field.

4. The urban engineering planning method based on UAV remote sensing mapping as described in claim 3, characterized in that: By combining the terrain accessibility field and the obstacle repulsion potential field, a steady-state diffusion equation is constructed. The Gauss-Seidel iterative method is used to solve the steady-state diffusion equation, resulting in a dynamic cost surface, including: By combining the terrain accessibility field and the obstacle repulsion potential field, a steady-state diffusion equation is constructed and solved using the Gauss-Seidel iterative method. A stopping threshold is set using the residual criterion. The iteration stops when the L2 norm of the residual vector of two consecutive iterations is less than the stopping threshold, thus obtaining the dynamic cost surface.

5. The urban engineering planning method based on UAV remote sensing mapping as described in claim 4, characterized in that: The process involves converting the dynamic cost surface into a comprehensive travel velocity field for each point in space, constructing an equation, solving the equation using the fast travel method to obtain the minimum arrival time, calculating the negative gradient direction of the minimum arrival time field, converting it into the next path point, and generating an optimal path sequence, including: The dynamic cost surface is converted into a comprehensive travel velocity field at various points in space. An equation is constructed, and the equation is solved using the fast travel method to obtain the minimum arrival time. Initialize the path sequence, calculate the negative gradient direction of the minimum arrival time field, and convert it into the next path point; When the Euclidean distance between a newly generated path point and the initial path point is less than one backtracking step, the initial path point is added to the sequence as the last point, the backtracking terminates, and the path point coordinate sequence is obtained. The sequence is then reversed and defined as the optimal path sequence.

6. The urban engineering planning method based on UAV remote sensing mapping as described in claim 5, characterized in that: The construction of a three-dimensional engineering entity model based on the optimal path sequence includes: Connect the optimal path sequence to form the initial centerline geometric entity, and perform cubic spline smoothing to obtain the final path centerline; Using the path centerline as the axis, perform Euclidean buffer analysis to generate a fixed-distance buffer zone, which is defined as the engineering land area. The binary map of the obstacle is spatially overlaid with the area of ​​the engineering land area to calculate the area encroached by the obstacle. Design station points are set up along the path centerline at fixed pile intervals. The elevation values ​​of the design station points are extracted from the digital surface model, and Kriging interpolation is used to generate a continuous design elevation surface. The two-dimensional boundary polygon of the engineering site is stretched vertically to the design elevation surface to generate a closed three-dimensional engineering solid model.

7. The urban engineering planning method based on UAV remote sensing mapping as described in claim 6, characterized in that: The process of packaging data from a 3D engineering entity model includes: The 3D engineering solid model is triangulated, and the output file is automatically named according to the rule "Project ID_Deliverable Type_Timestamp".

8. The urban engineering planning method based on UAV remote sensing mapping as described in claim 1, characterized in that: The process of collecting and preprocessing multi-source data includes: Once the planning task is triggered, drones are dispatched to survey the target planning area, collect multi-source data of the city, and perform noise reduction and standardization processing.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, it implements the steps of the urban engineering planning method based on UAV remote sensing mapping as described in any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by the processor, it implements the steps of the urban engineering planning method based on UAV remote sensing mapping as described in any one of claims 1 to 8.

Images