Land reclamation engineering data analysis and optimization method considering terrain constraints

Abstract

The application relates to the field of land consolidation and engineering calculation technology, and discloses a land consolidation engineering data analysis optimization method considering terrain constraints, which comprises the following steps: obtaining an original digital elevation model matrix and engineering constraint parameters, extracting bottom layer differential calculation features to calculate grid microscopic geometric features; based on discrete Morse theory, terrain feature topology is denoised and constructed, and a critical earthwork volume constant is calibrated; an optimal transport model with entropy regularization is constructed for dynamic optimization, real-time monitoring of earthwork transfer amount is carried out, and when a set proportion is reached, local terrain feature topology reconstruction and boundary stitching are triggered, and a global optimal allocation matrix is output; finally, the matrix is reversely mapped to generate a vector flow field, and a limited Laplace grid smoothing operation is carried out to output construction-level engineering data. The application effectively combines mechanical operation constraints and a topological dynamic reconstruction mechanism, and improves the on-site executability of the earthwork allocation scheme.

Description

Technical Field

[0001] This invention relates to the field of land consolidation and engineering calculation technology, specifically to a method for analyzing and optimizing land consolidation engineering data while taking into account topographic constraints. Background Technology

[0002] In land consolidation and site leveling projects, earthwork allocation is a core element affecting construction costs and project timelines. With the application of digital elevation models (DEMs), most existing earthwork allocation algorithms are based on grid partitioning and linear programming models to calculate earthwork transfer volumes and allocation paths. However, traditional algorithms have significant limitations when dealing with complex natural terrain. Existing technologies typically use spatial distance as the sole criterion for calculating allocation weights, failing to effectively consider the physical constraints faced by construction machinery during on-site operations, such as maximum climbing limits. This calculation mode, which relies solely on distance, often results in allocation paths that are unenforceable in areas with undulating terrain.

[0003] Meanwhile, during earthwork transfer calculations, the actual elevation of the ground surface dynamically changes as the cutting and filling operations progress, thereby altering the original topological connectivity. Most existing allocation algorithms are based on static terrain assumptions and lack a topology adaptive reconstruction mechanism that adapts to terrain deformation. When attempting to introduce phased local elevation updates during the calculation process, errors often arise due to mismatches between computational grid variations and solution matrix sizes. To address the dimension alignment problem, some existing solutions employ a global parameter reset and iterative approach from scratch. This approach wastes early computational parameters and convergence states, reducing the algorithm's computational efficiency when processing high-resolution engineering data.

[0004] Furthermore, discretized grid allocation calculations are prone to creating step-like elevation abrupt changes at the boundaries. Existing grid smoothing algorithms lack rigid constraints on the engineering design elevation baseline when eliminating these abrupt changes, which can easily lead to the smoothed terrain deviating from the original engineering design benchmark. Therefore, a data analysis and optimization method is needed that can take into account mechanical and physical operational constraints and has the ability to dynamically reconstruct the topology and perform parameter hot-start, so as to output data that meets the requirements of on-site construction. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a land consolidation engineering data analysis and optimization method that takes into account terrain constraints. This method solves the problem that existing land consolidation earthwork allocation algorithms fail to effectively integrate the physical operational constraints of engineering machinery in complex terrain, resulting in planned earthwork allocation paths being infeasible in actual construction. Furthermore, when earthwork transfer causes dynamic changes in terrain elevation, existing algorithms lack an adaptive topological reconstruction mechanism, easily leading to calculation errors such as mismatched solution matrix dimensions and wasted computational resources for early iteration parameters. Consequently, their output results cannot be directly converted into engineering data to guide actual construction operations.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for data analysis and optimization of land consolidation projects that takes into account topographic constraints, comprising the following steps: Obtain the original digital elevation model matrix and engineering constraint parameters of the target area, and use the finite difference method to extract the low-level differential calculation features of the original digital elevation model matrix to calculate the micro-geometric features of the grid. Based on the discrete Morse theory of continuous homology, feature extraction and denoising are performed on the original digital elevation model matrix to construct a topographic feature map. The critical earthwork volume constant is calibrated by surface integral algorithm. The critical earthwork volume constant refers to the value that reflects the maximum amount of engineering work required to completely level a single local micro-topography, using the boundary of the local micro-topography controlled by the critical point retained in the topographic feature map as the integration solution domain. The difference between the actual surface elevation in the domain and the local target flat field elevation set to eliminate the topographic feature of the critical point is used for surface integral. Calculate the physical resistance weights of the terrain feature topology map based on the engineering constraint parameters, construct the physical constraint cost matrix, and initialize the global supply and demand distribution and iteration variables; An optimal transport model with entropy regularization is constructed, and a local dynamic optimal transport calculation with a hot-start mapping mechanism is performed to monitor the earthwork transfer volume. When the cumulative earthwork transfer volume reaches the set proportion of the critical earthwork volume constant, the reconstruction of the local terrain feature topology map and boundary stitching are triggered. The local micro-topography region controlled by the node with the change is extracted, and the first-order elevation gradient field is used as the constraint benchmark to map the cumulative earthwork transfer volume into the elevation change increment of each discrete grid according to the gradient inverse weight distribution to update the actual surface elevation. The local terrain feature topology sub-map is regenerated for the local region after the elevation update. According to the Euclidean distance nearest neighbor principle of spatial coordinates, the undirected edge pointers of the topological streamlines suspended at the edge of the sub-map are rebound and connected to the breakpoint streamline nodes at the boundary of the unupdated region to complete the boundary stitching. After the boundary stitching is completed, for the new nodes generated by the local update, extract the dual scale vectors of all unchanged nodes that are first-order adjacent to the topological streamline. Use the modified inverse distance weighted interpolation formula with the introduction of a preset minimum positive smoothing constant to assign initial values, and start the optimal transmission iteration from the updated and reorganized set of dual scale vectors until the preset convergence condition is reached, and output the global optimal allocation matrix. The globally optimal allocation matrix is ​​transmitted and reverse-mapped to generate a directed earthwork allocation vector flow field, and a restricted Laplace mesh smoothing operation is performed to output construction-level engineering data.

[0007] Preferably, the step of obtaining the original digital elevation model matrix and engineering constraint parameters of the target area, and calculating the micro-geometric features of the grid, includes: obtaining the engineering constraint parameters, which include the set of polygon boundaries of the engineering red line, the maximum allowable longitudinal slope, the maximum single cut-and-fill depth, and the physical operation error threshold; for the effective pixels in the original digital elevation model matrix located inside the set of polygon boundaries of the engineering red line and after filtering out missing outliers, calculating the partial derivatives in the horizontal and vertical directions using the finite difference method, and synthesizing the first-order elevation gradient field and geometric slope; constructing a second-order Hessian matrix using the elevation information of adjacent nodes, solving for eigenvalues ​​to obtain the mean curvature and Gaussian curvature, and combining the first-order elevation gradient field and geometric slope to form the micro-geometric features of the grid.

[0008] Preferably, the step of extracting features and denoising the original digital elevation model matrix based on the discrete Morse theory of continuous homology to construct a terrain feature topology map includes: extracting local maxima, local minima, and saddle points of the terrain based on the micro-geometric features of the grid; calculating the elevation difference between each extreme point and its associated saddle point; merging extreme points with elevation differences lower than the physical operation error threshold and their associated topological streamlines into adjacent deep catchment areas or slope-ridge associated areas to perform topology denoising; retaining the effective critical points and streamlines after denoising to form the terrain feature topology map, and defining the retained local maxima as earthwork supply source nodes and the retained local minima as earthwork demand sink nodes.

[0009] Preferably, the step of calibrating the critical earthwork volume constant using the surface integral algorithm includes: using the local micro-topography boundary controlled by each retained critical point in the topographic feature map as the integration solution domain, wherein the boundary of the integration solution domain is determined by a closed surface region enclosed by the topographic streamlines connecting the critical point and its adjacent associated saddle points; performing surface integral using the difference between the actual surface elevation within the integration solution domain and the local target flat field elevation set to eliminate the topographic features of the critical point, to obtain the critical earthwork volume constant reflecting the maximum amount of engineering work required to completely flatten a single local micro-topography.

[0010] Preferably, the step of calculating the physical resistance weights of the terrain feature topology map based on the engineering constraint parameters and constructing the physical constraint cost matrix includes: extracting the spatial physical paths corresponding to the topologically connected edges in the terrain feature topology map in the original digital elevation model matrix; combining the physical horizontal distance of the spatial physical paths with the spatial average longitudinal slope, and setting a piecewise nonlinear penalty function based on the maximum permissible longitudinal slope: when the spatial average longitudinal slope is less than or equal to the maximum permissible longitudinal slope, the physical resistance weights are calculated by combining the physical horizontal distance with the square of the ratio of the terrain sensitivity penalty coefficient to the slope; when the spatial average longitudinal slope exceeds the maximum permissible longitudinal slope, the physical resistance weights are forcibly mapped to a maximum penalty constant representing the physically inaccessible state, so as to cut off infeasible earthwork allocation paths; and the physical constraint cost matrix is ​​composed of the set of physical resistance weights between all nodes.

[0011] Preferably, the step of initializing the global supply and demand distribution and dual scale vector includes: performing discrete volume summation based on a three-dimensional mesh to calculate the solid volume of each earthwork supply source node region, and multiplying it by the engineering earthwork loosening or compaction coefficient to obtain the effective virtual volume, and constructing a global supply distribution vector according to the node index order; calculating the volume gap where the actual ground elevation in the control area of ​​each earthwork demand sink node is lower than the local target leveling elevation, and generating a global demand distribution vector; generating a first dual scale vector with the same dimension as the global supply distribution vector, and a second dual scale vector with the same dimension as the global demand distribution vector, and setting the initial values ​​of all internal elements to a constant 1.

[0012] Preferably, the step of triggering local terrain feature topology map reconstruction and boundary stitching includes: extracting the local micro-topography region controlled by the changed node; using the first-order elevation gradient field as a constraint benchmark, mapping the accumulated earthwork transfer volume into the elevation change increment of each discrete grid according to the gradient inverse weight distribution, and correcting the actual surface elevation; regenerating the local terrain feature topology sub-map for the local region after elevation update; searching for the suspended topological streamlines at the edge of the local terrain feature topology sub-map; and rebinding its undirected edge pointers to the breakpoint streamline nodes at the boundary of the unupdated region according to the Euclidean distance nearest neighbor principle of spatial coordinates, thereby completing the boundary stitching.

[0013] Preferably, after the boundary stitching is completed, the step of performing local dynamic optimal transport calculation with a hot-start mapping mechanism includes: re-extracting the physical resistance weights of the mutated connected edges and reconstructing the local physical constraint cost matrix of the updated region; recalculating the supply and demand vectors of the changed nodes to align with the global mathematical dimension; performing hot-start mapping of the dual scale vector space: for the original nodes that have not changed, directly inheriting the dual scale vector of the previous state before the process suspension; for the new nodes generated by the local update, locating all the unchanged nodes that are first-order adjacent to them on the stitched topological streamline, extracting the dual scale vectors of the adjacent nodes, and assigning initial values ​​using the modified inverse distance weighted interpolation formula that introduces a preset minimum positive smoothing constant; resuming the optimal transport iteration starting from the updated and recombined set of dual scale vectors until the supply and demand edge residuals meet the preset approximation accuracy conditions or the internal loop reaches the preset maximum safe iteration number, and outputting the global optimal allocation matrix.

[0014] Preferably, the step of transmitting and reverse mapping the global optimal allocation matrix to generate a directed earthwork allocation vector flow field includes: extracting the sequence of non-zero elements in the global optimal allocation matrix; taking the three-dimensional physical coordinates of the earthwork supply source node as the starting point, the three-dimensional physical coordinates of the corresponding earthwork demand sink node as the ending point, and the connected topological streamlines as the spatial skeleton, assigning the values ​​of the sequence of non-zero elements to the spatial skeleton as physical transfer quantities; strictly calibrating the direction of the spatial transfer vector from the source node to the sink node, and reverse mapping back to the three-dimensional geometric space to generate the directed earthwork allocation vector flow field.

[0015] Preferably, the step of performing restricted Laplace mesh smoothing calculation and outputting construction-level engineering data includes: calculating the spatial weight parameters of each adjacent node using the inverse distance weighting method, performing normalization processing before substituting them into the iterative calculation, and calculating the standard smoothing intermediate elevation; using the engineering design elevation as the base, performing real-time calculations based on the product of the elevations of adjacent nodes and the maximum allowable longitudinal slope over the span of the grid physical resolution, and deducing and determining the dynamic upper and lower limits of the node elevations; using the dynamic upper and lower limits of the node elevations to rigidly truncate the standard smoothing intermediate elevations, and outputting the updated grid elevations; extracting the updated grid elevations and plane coordinate parameters to render and generate a three-dimensional grid elevation map, summarizing the discrete volumes of each grid region to compile and generate a cut and fill detail table, and combining the directed earthwork allocation vector flow field decoding to output the engineering machinery travel path matrix, and encapsulating the above data into the construction-level engineering data for output presentation.

[0016] This invention provides a data analysis and optimization method for land consolidation projects that takes into account topographic constraints. It has the following beneficial effects: 1. This invention constructs a terrain feature topology map by extracting the differential features of the elevation model, and sets a piecewise nonlinear penalty function based on engineering constraint parameters such as the maximum allowable longitudinal slope when constructing the physical constraint cost matrix. It assigns a large penalty constant to physical paths that exceed the mechanical climbing limit. This processing method directly cuts off physically infeasible spatial connections at the bottom data input stage, ensuring that the globally optimal allocation matrix output by calculation can objectively reflect the actual operating capability of engineering machinery in real complex terrain, and improves the executability of the allocation scheme.

[0017] 2. This invention monitors the earthwork transfer volume in real time during the optimal transmission iteration process. When the local terrain deformation reaches a set proportion, it triggers topology reconstruction and boundary stitching. It also uses a modified inverse distance weighted interpolation to perform spatial hot-start mapping of the dual scale vector. This dynamic reconstruction mechanism solves the problem of mismatch in the dimensions of the computational matrix caused by terrain changes due to construction filling and cutting. At the same time, the hot-start mapping allows the grid nodes after local mutation to directly inherit the previous computational state, avoiding global reset calculation, effectively saving computational resources and ensuring the convergence speed of the algorithm.

[0018] 3. This invention introduces a restricted Laplace grid smoothing operation, which is limited by the engineering design elevation and the maximum allowable longitudinal slope, after generating a directed earthwork allocation vector flow field through transmission reverse mapping. This operation is rigidly truncated by dynamically setting the upper and lower limits of the node elevation, which not only eliminates the elevation abrupt changes caused by discrete grid block allocation at the boundary, but also prevents the terrain elevation from violating the engineering design benchmark due to unrestricted smoothing. Thus, the theoretical calculation data is accurately converted into a three-dimensional grid elevation map and engineering machinery driving path matrix that can be directly used for on-site guidance. Attached Figure Description

[0019] Figure 1 This is a flowchart of the method of the present invention; Figure 2 A detailed flowchart of the target region data initialization and micro-geometric feature extraction of this invention; Figure 3 This is a detailed flowchart of the construction of the discrete Morse topology graph and the calibration of the critical volume in this invention; Figure 4 This is a detailed flowchart of the physical constraint cost matrix construction and global initialization of the present invention; Figure 5 This is a detailed flowchart of the local dynamic optimal transmission and hot start mapping calculation of the present invention; Figure 6 This is a detailed flowchart of the transmission reverse mapping and constrained grid smoothing output of the present invention; Figure 7 This is a comparison chart of the core performance indicators of the traditional algorithm of this invention. Detailed Implementation

[0020] The technical solutions in 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.

[0021] See attached document Figure 1 , Figure 1 This is a flowchart of a method according to an embodiment of the present invention. The present invention provides a method for data analysis and optimization of land consolidation projects that takes into account topographic constraints, comprising the following steps: S1. Obtain the original digital elevation model matrix and engineering constraint parameters of the target area, and calculate the micro-geometric features of the grid. In this embodiment, the computing device reads the original digital elevation model matrix of the target area through the data interface, and simultaneously receives the set of polygonal boundaries of the engineering red line, the maximum allowable longitudinal slope (which is usually limited by the physical climbing performance of the mechanical tracks and can be set between 15% and 30%), the maximum single excavation and filling depth, and the physical operation error threshold. Specifically, the physical operation error threshold is not a subjective assumption, but is calibrated based on the minimum operating height difference accuracy of engineering equipment such as bulldozers, for example, 0.1 meters. The processor uses the finite difference method to calculate the first-order elevation gradient field and the second-order Hessian matrix for the effective pixels in the original digital elevation model matrix, and then obtains the geometric slope and curvature data of the grid nodes by calculating the matrix eigenvalues, thus completing the conversion of physical spatial data into underlying computational features. S2, based on the discrete Morse theory of continuous cohomology, performs feature extraction and denoising on the original digital elevation model matrix, constructs a topological map of terrain features, and calibrates the critical earthwork volume constant. Since the original digital elevation model matrix often contains subtle numerical fluctuations, the system is prone to generating pseudo-extremes after extracting local maxima, local minima, and saddle points from the geometric feature data. To eliminate such redundant nodes, the system introduces a Morse hidden-face elimination mechanism. The processor calculates the elevation difference between each extreme point and its associated saddle point, merging extreme points below the physical operation error threshold and their associated streamlines into adjacent catchment areas. After denoising, the system constructs... A topological map of terrain features composed of critical points and streamlines is constructed. For the retained nodes, the critical earthwork volume constant required to change the local bottom-level topology is further calculated by surface integral. The critical earthwork volume constant specifically refers to the physical value obtained by surface integral of the difference between the actual surface elevation in the domain and the local target flat field elevation set to eliminate the topological features of the critical point, using the closed surface region enclosed by the topological streamlines connecting the critical point and its adjacent associated saddle points as the integration solution domain. The critical earthwork volume constant physically represents the maximum amount of engineering work required to completely level a single local micro-topography. S3 calculates the physical resistance weights of the terrain feature topology map based on engineering constraint parameters, constructs the physical constraint cost matrix, and initializes the global supply and demand distribution and iteration variables; the system traverses the edges in the terrain feature topology map and extracts the corresponding spatial physical paths in the original digital elevation model matrix; combining the physical horizontal distance of the path and the spatial average longitudinal slope, a penalty function is set according to the maximum allowable longitudinal slope, and the initial physical constraint cost matrix is ​​generated through nonlinear mapping calculation.

[0022] As a preferred approach, when the spatial average longitudinal slope of the path exceeds the maximum allowable longitudinal slope, the path weight tends to infinity, thereby directly trunculating infeasible paths at the spatial accessibility level; simultaneously, the processor calculates the total available earthwork and the total fill required in each node area, generates a global supply distribution vector and a demand distribution vector, and initializes the dual scale vector for the subsequent optimal transmission algorithm, with the initial value set to an all-1 vector. S4 executes local dynamic optimal transport calculation with a hot-start mapping mechanism, monitors earthwork transfer volume and triggers local terrain feature topology reconstruction and boundary stitching, outputting a globally optimal allocation matrix; the system constructs an optimal transport model with entropy regularization and solves iteratively through alternating projection; during the iteration process, the processor accumulates the earthwork transfer volume of each node in real time; when the accumulated earthwork transfer volume of a specific node reaches a set proportion of the critical earthwork volume constant, the set proportion is preferably between 0.85 and 0.95 to reserve a calculation buffer, it indicates that the macroscopic geometric features of the terrain at that location are about to undergo a fundamental reversal, at which point the system suspends the current transport process; the processor, based on the transferred earthwork transfer volume, ... Using the first-order elevation gradient field as a constraint benchmark, the accumulated earthwork transfer volume is mapped to the elevation change increment of each discrete grid according to the gradient inverse weight distribution. This updates the elevation of the local grid, and feature extraction and denoising are re-executed to generate a local terrain feature topology subgraph. The boundary of this local terrain feature topology subgraph is then stitched with the global terrain feature topology region that has not undergone elevation change. Specifically, based on the Euclidean distance nearest neighbor principle of spatial coordinates, the undirected edge pointers of the topological streamlines suspended at the edge of the subgraph are reconnected to the breakpoint streamline nodes at the boundary of the unupdated region. After the local terrain feature topology is reconstructed, to ensure the mathematical convergence of the algorithm, the system performs a hot-start mapping of the dual scale vector space. For unchanged nodes, the dual scale vector of the previous state is directly inherited. For new nodes generated by local updates, the system uses the dual scale vector of the first-order adjacent unchanged nodes on the topological flow line to perform a modified inverse distance weighted interpolation operation with a preset minimum positive smoothing constant to assign initial values. The system uses the updated set of dual scale vectors as the starting point to recover the optimal transmission iteration until the supply and demand edge residuals meet the accuracy conditions or reach the preset maximum safe iteration number, and outputs a globally optimal allocation matrix containing the allocation amount between nodes. S5 transmits and reverse-maps the globally optimal allocation matrix, performs constrained mesh smoothing, and outputs construction-level engineering data. Since the globally optimal allocation matrix only reflects the macroscopic mathematical flow direction between topological nodes, the system extracts the sequence of non-zero elements from the globally optimal allocation matrix, combines it with the source-sink topological edge relationship, and reverse-maps it back to the three-dimensional geometric space to generate a directed earthwork allocation vector flow field. To eliminate abrupt elevation jumps in the mesh generated by theoretical calculations, the processor performs constrained Laplace mesh smoothing within each mesh region. This constrained Laplace mesh smoothing is subject to hard constraints of engineering design elevation and maximum allowable longitudinal slope to prevent excessive smoothing from causing the terrain to deviate from the design baseline. After smoothing, the system outputs a three-dimensional mesh elevation map, cut and fill details, and engineering machinery travel path matrix to an external storage medium or display terminal.

[0023] See attached document Figure 2 , Figure 2This is a detailed flowchart of target area data initialization and micro-geometric feature extraction according to an embodiment of the present invention.

[0024] In step S1, which involves obtaining the original digital elevation model matrix and engineering constraint parameters of the target area and calculating the macroscopic generalization of the microscopic geometric features of the grid, this embodiment further refines it into the following sub-steps to clarify the specific implementation logic of transforming physical space data into underlying computational features.

[0025] S101, in this embodiment, the computing device loads the original digital elevation model matrix of the target area through the hardware data interface, and simultaneously receives vector coordinate data representing the set of polygonal boundaries of the engineering red line. As for the specific underlying parsing mechanism of the computing device reading multi-source spatial data files and performing coordinate system registration through the data interface, those skilled in the art can call the standard open-source environment of geographic information system to implement it. Its data reading process is a well-known technology in the field and will not be described in detail here. While loading spatial data, the system receives engineering constraint parameters input from the outside. As a preferred method, the engineering constraint parameters specifically include the maximum allowable longitudinal slope, the maximum single cut and fill depth, and the physical operation error threshold. To accurately reflect the real-world land reclamation environment, the maximum permissible longitudinal slope is limited by the physical climbing performance limits of the tracks of bulldozers and excavators at the construction site, with a range of 15% to 30%. Meanwhile, the physical operation error threshold is set to 0.1 meters based on the minimum operating height difference accuracy calibration of specific engineering equipment. The system establishes a two-dimensional Cartesian coordinate system in memory, mapping the original digital elevation model matrix into a set of regular grid nodes, recording the planar coordinates of any node in the grid and its corresponding initial ground elevation, thus establishing the physical boundary conditions for subsequent algorithm calculations.

[0026] In step S102, to extract the local spatial variation trend of the terrain, the processor performs partial derivative calculations on the effective pixels within the original digital elevation model matrix. Here, effective pixels specifically refer to actual elevation data points located within the engineering redline polygon boundary set, after filtering out missing outliers, to avoid numerical calculations getting stuck in dead zones. Since the original digital elevation model matrix is ​​a discrete set of grid points, the processor uses the finite difference method to calculate the first-order elevation gradient field. The system traverses each internal grid node within the engineering redline polygon boundary set, using the central difference algorithm to calculate the ratio of the elevation difference between horizontally adjacent nodes to the horizontal grid's physical resolution, and the ratio of the elevation difference between vertically adjacent nodes to the vertical grid's physical resolution, respectively, to obtain the horizontal and vertical first-order partial derivatives of the grid node. The processor synthesizes these two orthogonal first-order partial derivatives to obtain the first-order elevation gradient field of the node, and simultaneously calculates the geometric slope in pure scalar dimension based on the vector magnitude of this first-order elevation gradient field, thus intuitively reflecting the degree of tilt of the ground surface. S103, based on the first-order elevation gradient field, the processor solves for the second-order differential components reflecting the unevenness of the terrain surface. For each grid node in the original digital elevation model matrix, the processor, based on finite difference calculation logic, uses the elevation information of adjacent nodes to calculate their second-order partial derivatives in the horizontal, vertical, and intersecting directions, thereby constructing a second-order Hessian matrix characterizing the local terrain curvature. After constructing the second-order Hessian matrix, the processor obtains two independent real eigenvalues ​​by solving the characteristic equation of the second-order Hessian matrix. The system adds these two real eigenvalues ​​and takes the arithmetic mean to obtain the average curvature of the grid node. Simultaneously, it multiplies these two real eigenvalues ​​to obtain the Gaussian curvature of the grid node, thus reflecting the terrain curvature. From a theoretical perspective, the mean curvature reflects the overall undulation of the terrain surface, while the positive or negative attribute of Gaussian curvature constitutes the underlying physical criterion for determining whether a local terrain belongs to a mountain peak, depression, or saddle surface. The calculation process of this curvature and the finite difference operation are both mature and well-known numerical analysis techniques, and no additional formula expansion is required here. Based on the first-order elevation gradient field, geometric slope, and spatial curvature data containing mean curvature and Gaussian curvature obtained by the finite difference method, the grid micro-geometric features are formed together. At this point, the computing device has completed the transformation of the original digital elevation model matrix from three-dimensional physical coordinates to underlying differential calculation features. The output underlying differential calculation features constitute the basic input data for triggering subsequent feature dimensionality reduction operations using discrete Morse theory.

[0027] See attached document Figure 3 , Figure 3 This is a detailed flowchart of the construction of a discrete Morse topology graph and the calibration of the critical volume according to an embodiment of the present invention.

[0028] In step S2, which involves feature extraction and denoising of the original digital elevation model matrix based on discrete Morse theory with continuous cohomology, constructing a topological map of terrain features and calibrating the critical earthwork volume constant, this embodiment further refines it into the following sub-steps to clarify the underlying mapping relationship between abstract topological mathematics and actual engineering physical quantities.

[0029] S201, In this embodiment, the system extracts local maxima, local minima, and saddle points of the terrain based on the micro-geometric features of the grid obtained in the aforementioned steps. Specifically, based on Discrete Morse theory, if the elevation of a grid node in the original digital elevation model matrix is ​​greater than the elevations of all its first-order adjacent nodes, the processor determines it as a local maximum; if it is less than the elevations of all its first-order adjacent nodes, it is determined as a local minimum. For the traversal algorithm for extracting critical points on the discrete grid, those skilled in the art can use conventional neighborhood elevation comparison logic to implement it. The basic search process is a well-known technology in the field and will not be described in detail here. To establish a mapping relationship between pure mathematical graph theory and actual physical construction, the system defines the extracted local maxima as the source node of earthwork supply and the local minima as the sink node of earthwork demand. The elevation gradient descent path connecting the nodes is the topological streamline, which represents the spatial operation path of the engineering machinery for earthwork transfer.

[0030] S202. Because the original digital elevation model matrix often contains subtle numerical fluctuations, it is prone to producing pseudo-extreme values. To eliminate such redundant nodes and avoid meaningless iterations in the subsequent algorithm optimization process, the system introduces a Morse hidden point elimination mechanism. The processor calculates the elevation difference between each extreme point and its associated saddle point. Its mathematical calculation logic is defined as follows: ; In the formula, This represents the elevation value of the detected extreme point. This represents the elevation value of the saddle point that is closely adjacent to the extreme point in terms of topological connectivity; the system will calculate the elevation difference. A rigorous comparison is performed with the physical operation error threshold loaded in the aforementioned steps; when the elevation difference... When the value falls below the physical operation error threshold, the processor determines that the extreme point is a numerical artifact of non-real terrain undulation. The system then performs topology denoising, merging the extreme point below the error threshold and its associated topological streamlines into adjacent deep catchment areas (i.e., topological control regions with local minimum points as sink nodes) or hillside and ridge association areas (i.e., topological control regions with local maximum points as source nodes) based on the maximum or minimum values ​​of their extreme attributes. After denoising, the system formally constructs a terrain feature topology map composed of effective critical points and streamlines. This terrain feature topology map is equivalent to a Morse-Smail complex in its underlying data structure, providing a basic spatial topological relationship for subsequent earthwork allocation calculations.

[0031] S203, for the effective critical points retained after denoising, the system further calculates the critical earthwork volume constant required to change the local bottom-level topology. As a preferred approach, the processor uses the local micro-topography boundaries controlled by each retained critical point as the integration solution domain. The maximum amount of engineering work required to level the region is calculated using a surface integral algorithm; in this embodiment, the integral solution domain... The specific boundary is determined by the closed surface region enclosed by the topological streamlines connecting the critical point and its adjacent associated saddle points, and its integral calculation formula is as follows: ; In the formula, and Denote the solution domain of the integral, respectively. The horizontal and vertical coordinates of the two-dimensional plane within; This represents the current actual ground elevation at this plane coordinate location; The local target flat-field elevation is set to eliminate the topological features of the critical point; in the actual physical state simulation, this local target flat-field elevation... The value is taken as the elevation of the saddle point that has a spatial connection with the critical point being solved; from the perspective of the topological cancellation principle of discrete Morse theory, leveling the elevation of a local peak or depression to the elevation level of its associated saddle point is the critical condition for the complete disappearance of the local geometric feature; this critical earthwork volume constant. In a physical sense, it establishes the maximum amount of engineering work required to completely flatten a single local micro-topography; by transforming the abstract topological deformation into a specific earthwork volume value, this parameter provides a quantitative physical criterion for triggering local topological reconstruction by the subsequent dynamic optimal transmission algorithm, thereby supporting the various functional features of dynamic monitoring and reconstruction in the claims.

[0032] See attached document Figure 4 , Figure 4 This is a detailed flowchart of the construction and global initialization of the physical constraint cost matrix according to an embodiment of the present invention.

[0033] In step S3, which involves calculating the physical resistance weights of the terrain feature topology map based on engineering constraint parameters, constructing the physical constraint cost matrix, and initializing the macroscopic generalization of global supply and demand distribution and iteration variables, this embodiment further refines it into the following sub-steps to clarify the specific implementation logic of transforming actual engineering machinery physical constraints into underlying graph theory mathematical weights.

[0034] S301, the system traverses each topologically connected edge in the denoised terrain feature topology map and extracts its corresponding spatial physical path in the original digital elevation model matrix. For tracking the shortest spatial path between two points in the discrete grid model and calculating the physical horizontal distance, those skilled in the art can use the conventional Bressenham line algorithm or Dijkstra path search logic. The extraction process of the path horizontal distance is a well-known technology in the field and will not be described in detail here.

[0035] After obtaining the physical horizontal distance of the topological edge, the processor calculates the spatial average longitudinal slope of the path by combining the elevation difference between the two grid nodes. From the basic physical laws of engineering machinery operation, the energy consumption and driving resistance of mechanical equipment in climbing mode do not increase linearly, but rather exhibit a significant nonlinear increasing trend with the increase of the slope. To realistically reflect the physical driving constraints of engineering machinery in complex terrain in the underlying mathematical model, the system sets a piecewise nonlinear penalty function based on the maximum allowable longitudinal slope loaded in the aforementioned steps, and thereby constructs a physical constraint cost matrix. In this example, for any source node... With sink node The physical resistance weight of the topological edge between them The mathematical computational logic is defined as follows: ; In the formula, Indicates the source node With sink node The physical horizontal distance between the spatial physical paths; This represents the average longitudinal slope of the physical path in space. This indicates the maximum permissible longitudinal slope set due to limitations in the physical climbing performance of the mechanical tracks; A terrain sensitivity penalty coefficient is used to calibrate the degree of impact of terrain undulation on mechanical energy consumption; as a preferred method, this terrain sensitivity penalty coefficient... The value is set to a range of 1.5 to 3.0, and its specific value is determined by looking up a table based on the engineering soil characteristics and the rated power of the machinery in the input system. For example, for the working combination of high-cohesive soil and low-power bulldozer, the processor will adaptively extract a higher penalty coefficient close to 3.0 to amplify the resistance sensing. This represents the maximum penalty constant for physically inaccessible states. It should be noted that when the spatial average longitudinal slope of the path... Exceeding the maximum permissible longitudinal slope At that time, if the physical resistance weight is directly applied in the underlying code... Assigning a strict mathematical positive infinity can easily lead to floating-point overflows or invalid numbers (NaN) in subsequent matrix iteration operations, creating a computational dead zone. Therefore, the system forcibly maps this state to a maximum penalty constant in the memory computation domain. The constant The value of needs to be much greater than the sum of the physical resistance weights of all normal connected paths in the matrix, for example, set to the power of 10 to the sixth power. This directly cuts off physically infeasible earthwork allocation paths at the spatial accessibility level. This mapping mechanism eliminates the possibility of the optimal transmission model generating incorrectly planned streamlines in areas with extreme elevation differences from the algorithm's bottom layer. It effectively avoids subsequent deductions that output operation instructions that violate engineering common sense. The set of physical resistance weights between all nodes in the terrain feature topology map together form the global physical constraint cost matrix that supports subsequent allocation operations.

[0036] S302, after establishing the physical resistance boundary for spatial allocation, the system simultaneously evaluates the earthwork capacity of each work node under the initial terrain condition. The processor calculates the total available earthwork volume of each source node region and the total fill volume required for each sink node region in the terrain feature topology map. Specifically, for any local micro-topography region controlled by a source node, the system performs discrete volume summation based on a three-dimensional mesh, accumulating the volume of the entity whose actual surface elevation exceeds the local target flat elevation, and uses this as the total available earthwork volume of that source node. Simultaneously, to conform to the incompressibility of earthwork in a real physical environment, the system introduces… The system uses the loosening or compaction coefficient of the earthwork in the project to perform physical volume conversion on the total available earthwork of the source nodes. That is, it multiplies the solid volume calculated by the source node area by the loosening or compaction coefficient to obtain the effective virtual volume that can be used for filling allocation, which is the final total available earthwork. By traversing all extracted source nodes and performing the above volume summation and conversion, the system arranges the total available earthwork of all source nodes in a fixed node index order to form a global supply distribution vector. Similarly, the system calculates the volume gap between the actual ground elevation and the local target leveling elevation within the control area of ​​each sink node, and generates a global demand distribution vector.

[0037] After completing the numerical transformation of the capacity space distribution, a neutral computational benchmark is constructed to support the subsequent dynamic optimal transmission optimization algorithm. The system assigns values ​​to the internal iterative variables. In this embodiment, the system sets the dual scale vector applied to the optimal transmission alternating projection solution logic as the initial state to be activated. The processor generates a first dual scale vector with the same dimension as the global supply distribution vector, and at the same time generates a second dual scale vector with the same dimension as the global demand distribution vector. The system strictly sets the initial values ​​of all elements in these two dual scale vectors to a constant of 1. This global initialization strategy of all-1 vectors can ensure the numerical stability and mathematical convergence of the first multiplication operation of the subsequent alternating projection iteration algorithm without introducing artificial spatial bias.

[0038] Thus, the numerically calibrated physical constraint cost matrix, the global supply and demand distribution vector reflecting the local construction volume, and the dual scale vector that has completed memory initialization together provide an initial numerical environment with a complete error prevention mechanism for the dynamic transmission simulation that triggers real-time transfer volume monitoring and local topology reconstruction in the next stage.

[0039] See attached document Figure 5 , Figure 5 This is a detailed flowchart of the local dynamic optimal transmission and hot start mapping calculation according to an embodiment of the present invention.

[0040] For step S4, which involves performing local dynamic optimal transport calculation with a hot-start mapping mechanism, monitoring earthwork transfer volume, triggering local terrain feature topology reconstruction and boundary stitching, and outputting a macro-general summary of the global optimal allocation matrix, this embodiment further refines it into the following sub-steps to fully disclose the underlying mathematical closed loop of the dynamic deduction algorithm.

[0041] S401, In this embodiment, in order to efficiently calculate the earthwork transfer scheme in the node group of the terrain feature topology map, the system constructs an optimal transmission model with entropy regularization. Considering that the traditional linear programming optimal transmission model is prone to the curse of dimensionality when solving high-dimensional discrete grids, the system introduces an information entropy regularization term into the optimal transmission objective function. From the underlying mathematical principle, this mechanism can transform the original linear programming problem into a strictly convex optimization problem, thereby smoothing the solution space and allowing the system to use matrix scaling technology to accelerate algorithm convergence. Specifically, the processor takes the global supply distribution vector, the global demand distribution vector, and the physical constraint cost matrix as inputs, calculates the global cumulative value of the product of the physical resistance weight and the earthwork allocation amount between the corresponding nodes, and subtracts the information entropy of the allocation matrix controlled by the regularization coefficient, thereby establishing the objective function to be optimized. As a preferred method, the value of the regularization coefficient is adaptively determined according to the numerical magnitude of the physical constraint cost matrix, usually taking a value of one-hundredth to one-tenth of the median of the elements inside the physical constraint cost matrix, so as to achieve a reasonable balance between iterative convergence speed and physical allocation accuracy.

[0042] After establishing the objective function, the processor uses the standard alternating projection algorithm to perform internal iterative approximation between the supply and demand sides. The processor will continue to execute this internal alternating projection until the computational residuals of the current supply and demand network state converge, and output the local allocation matrix of the current stage. For the underlying scalar multiplication and matrix update logic of this alternating projection process, those skilled in the art can configure it according to the standard nonnegative matrix scaling theory. Its basic iterative process is a well-known technology in the field and will not be described in detail here.

[0043] In S402, during the multi-stage optimization process, the processor accumulates the historical earthwork transfer volume of each specific node in real time. Specifically, after each stage of global convergence and output of the allocation matrix, the system sums the allocation elements involving specific source or sink nodes according to the transfer share specified in the matrix. The system compares and monitors the accumulated actual earthwork transfer volume with the critical earthwork volume constant associated with the node. When the accumulated earthwork transfer volume of a specific node reaches the set proportion of its critical earthwork volume constant, it indicates that the solid volume of the micro-topography at that location is close to the leveling threshold, and the macroscopic geometric features are about to undergo a fundamental reversal. To prevent the computational grid from collapsing due to the instantaneous collapse of the underlying topology during subsequent transmission, the system triggers an error-prevention suspension mechanism, immediately pausing the optimal transmission calculation for the next stage and retaining the dual scale vector state data of each node.

[0044] S403, after triggering the suspension mechanism, the processor updates the elevation of the local grid in 3D space based on the amount of earthwork transferred. In actual operation, the system extracts the local micro-topographic region controlled by the node that has changed. To prevent the original micro-first and second derivative continuity from being instantly destroyed due to simple volume averaging, the system uses the original first-order elevation gradient field as a constraint benchmark and maps the accumulated earthwork transfer amount to the elevation change increment of each discrete grid according to the gradient inverse weight distribution. Subsequently, the system synchronously corrects the actual surface elevation of the corresponding grid according to the physical operation logic of reducing the elevation of the source node and increasing the elevation of the sink node, so as to ensure that the deformed grid still has the continuity of conforming to the natural landform. After the elevation update is completed, the system only re-executes micro-geometric feature extraction and discrete Morse denoising for the local area where the elevation has changed, generates a local topographic feature sub-map, and stitches the boundary of this local topographic feature sub-map with the global topographic feature map area that has not changed in elevation.

[0045] S404. After reconstructing the local terrain feature topology map, due to local variations in the number of nodes and spatial connectivity, directly continuing the previous iteration would lead to matrix dimension errors. Therefore, the system first re-extracts the physical resistance weights of newly added or mutated connected edges based on the stitched latest topology map, reconstructs the local physical constraint cost matrix of the updated region, and recalculates the remaining available earthwork and fill volume of the changed nodes to align the latest mathematical dimensions of the global supply and demand vectors. After dimension alignment, to ensure the mathematical convergence of the algorithm and improve computational efficiency, the system performs a hot-start mapping of the dual scale vector space. For the original nodes that have not changed, the system instructs their direct inheritance process to suspend the dual scale vector of the previous state. For the new nodes generated by the local update, the system uses spatial correlation to assign them reasonable initial values. The system uses the updated and recombined tensor dimension and dual scale vector set as the starting point to resume the optimal transmission iteration of alternating projection. The processor continuously monitors the global supply and demand edge residuals until the supply and demand edge residuals meet the preset approximation accuracy conditions or the internal loop reaches the preset maximum safe iteration count, and outputs the globally optimal allocation matrix.

[0046] For new nodes generated by local updates, the system assigns them reasonable initial values ​​using spatial correlation. Specifically, the processor locates all unchanged nodes that are first-order adjacent to the new node on the stitched topological streamline, extracts the dual scale vectors of these adjacent nodes, and performs inverse distance weighted interpolation. Considering the numerical robustness of actual computing devices when handling extremely small distances, and to completely eliminate the program dead zone caused by spatially overlapping nodes resulting in a zero denominator, the system adopts a modified inverse distance weighted interpolation formula, the mathematical definition of which is as follows: ; In the formula, This represents the initial value of the dual scale vector to be assigned to the new node; This represents the total number of nodes whose first-order adjacency remains unchanged. Indicates the first The dual scaling vector of a first-order adjacent node with unchanged adjacency; This represents the topological path distance between the new node and its neighboring node; This represents a preset minimum positive smoothing constant. As a preferred approach, the value of this minimum positive smoothing constant is set to 10 to the power of negative 6.

[0047] The system starts with the updated and reorganized set of dual scale vectors, resumes the optimal transmission iteration of alternating projection, and the processor continuously monitors the supply and demand edge residuals on the global supply and demand sides until the supply and demand edge residuals meet the preset approximation accuracy condition or the internal loop reaches the preset maximum safe iteration number. At this point, the algorithm is determined to have reached the convergence state. As an optimal method, the value range of the approximation accuracy condition is set between 10 to the power of -3 and 10 to the power of -5, and the maximum safe iteration number is set between 500 and 2,000. The system stops iterating and outputs the global optimal allocation matrix containing the final calculated allocation amount between all nodes to the external computing memory.

[0048] See attached document Figure 6 , Figure 6 This is a detailed flowchart of transmission reverse mapping and constrained grid smoothing output according to an embodiment of the present invention.

[0049] In step S5, which involves transmitting and reverse mapping the global optimal allocation matrix, performing constrained mesh smoothing operations, and outputting a macroscopic summary of construction-level engineering data, this embodiment further refines it into the following sub-steps to clarify the implementation details of restoring the underlying mathematical abstraction to the real three-dimensional physical space.

[0050] S501, in this embodiment, the system reads the global optimal allocation matrix output by the aforementioned steps. Since the global optimal allocation matrix only reflects the macroscopic mathematical flow between the dimensionality-reduced topological nodes at the mathematical calculation level, the processor needs to restore it to the physical operation path that guides mechanical construction. The system extracts the sequence of non-zero elements in the global optimal allocation matrix. Each non-zero element represents a substantial earthwork transfer instruction in a physical sense. The processor combines these non-zero element sequences with the source-sink topological edge relationship and maps them back to the three-dimensional geometric space. Specifically, the system takes the three-dimensional physical coordinates of the earthwork supply source node as the starting point, the three-dimensional physical coordinates of the corresponding earthwork demand sink node as the ending point, and the connected topological streamlines as the spatial skeleton. The system assigns the values ​​in the global optimal allocation matrix to the spatial skeleton as physical transfer quantities and strictly marks the direction of the spatial transfer vector from the source node to the sink node. Thus, a directed earthwork allocation vector flow field is generated on the original digital elevation model matrix. This directed earthwork allocation vector flow field eliminates the ambiguity of the pure mathematical matrix in spatial position and constitutes the underlying trajectory data for subsequent guidance of the on-site driving and bulldozing operations of construction machinery.

[0051] S502, after generating the directed earthwork allocation vector flow field, since the underlying optimal transport algorithm is based on discrete topological nodes for block earthwork allocation, theoretical calculations are prone to causing discontinuous elevation steps at the boundaries of adjacent grids. To eliminate such discontinuities in the physical grid caused by numerical calculations, the processor performs restricted Laplace smoothing operations within each grid region. From the perspective of underlying geometric processing principles, Laplace smoothing, by introducing the weighted average elevation of adjacent nodes, can effectively diffuse and homogenize locally excessively concentrated earthwork volumes. However, ordinary numerical smoothing will uncontrollably smooth out terrain undulations, thereby causing the original engineering design benchmarks to fail. Therefore, this embodiment introduces hard constraints in the smoothing iteration model, that is, the smoothing operation must be strictly limited by the engineering design elevation and the maximum allowable longitudinal slope to prevent excessive smoothing from causing the terrain to violate the design bottom line. In this embodiment, the specific mathematical update logic of the restricted Laplace smoothing operation is defined as follows: ; In the formula, Indicates the standard smoothed intermediate elevation; This indicates the actual ground elevation of the current grid node; Indicates the smoothing relaxation coefficient; This represents the total number of neighboring nodes of the given grid node; This represents the spatial weight parameter of each adjacent node; in this embodiment, this spatial weight parameter... The inverse distance weighting method is used for calculation. Its value is inversely proportional to the Euclidean distance from the adjacent node to the current grid node, and strict normalization is performed before being substituted into the iterative calculation (i.e., ensuring...). To maintain the conservation of physical elevation, Indicates the elevation of adjacent nodes.

[0052] ; In the formula, This represents the updated mesh elevation output after constraint truncation. and These represent the dynamic upper and lower limits of the node elevation, respectively.

[0053] Specifically, the aforementioned dynamic upper and lower limits are based on the engineering design elevation initially imported by the system, and are calculated in real time by combining the product of the elevation of adjacent nodes and the maximum allowable longitudinal slope over the span of the grid's physical resolution. This ensures, mathematically, that the smoothed elevation difference between adjacent grids will never exceed the physical ramping limit of the engineering equipment. As a preferred method, this smoothing relaxation coefficient... The value range is set between 0.1 and 0.5 to ensure the numerical stability of the iterative update process. This restricted Laplacian grid smoothing operation eliminates the risks of slope instability and cut-and-fill boundary violations caused by conventional image smoothing logic in terrain processing by rigidly truncating the intermediate elevation of the standard smoothing algorithm from the bottom layer of the algorithm.

[0054] S503, after completing the above-mentioned underlying spatial smoothing process, the derivation data of the underlying mathematical model has been completely translated into entity parameters that conform to the standards of real physical construction operations. The system performs formatted compilation and encapsulation operations on the above data; the processor extracts and updates the grid elevation and plane coordinate parameters, and renders and generates a three-dimensional grid elevation map for intuitively displaying the results of terrain transformation; simultaneously, the system performs discrete volume summarization on the fill volume and cut volume of each grid area, and compiles and generates a detailed cut and fill table containing specific work volume indicators; combined with the directional earthwork allocation vector flow field generated in the previous steps, the system decodes the data channel and outputs the engineering machinery travel path matrix carrying the earthwork transfer volume and direction commands.

[0055] Finally, the system confirms the compiled 3D mesh elevation map, cut and fill details table, and engineering machinery travel path matrix as the final construction-level engineering data, and packages and outputs them to external storage media, or directly transmits them to a display terminal for visualization. As for the construction of the 3D rendering pipeline of the graphical interface and the underlying communication interface calling mechanism of the external storage media, those skilled in the art can use general computer graphics toolkits and standard file input / output protocols for environment configuration. The data rendering and disk storage process is a well-known technology in the field and will not be described in detail here.

[0056] Specific application examples are described. In this embodiment, three representative unleveled hilly plots are selected as target test areas. The computing device imports the original digital elevation model matrix of each area through the data interface. In the stage of extracting the micro-geometric features of the grid, the processor uses the finite difference method to calculate the first-order elevation gradient field and the second-order Hessian matrix. Based on the discrete Morse theory, the processor calculates the elevation difference between each extreme point and its associated saddle point, and performs topology denoising on extreme points that are lower than the set physical operation error threshold to extract effective source nodes and sink nodes.

[0057] During the calculation process, the system traverses the topological connected edges and assigns a large penalty constant to spatial paths exceeding the physical slope limit based on the preset maximum allowable longitudinal slope, thus constructing a physical constraint cost matrix. In the optimal transmission alternating projection iteration, the processor monitors the earthwork transfer volume of each node in real time. When the local cumulative transfer volume reaches the set proportion of the corresponding critical earthwork volume constant, the system suspends the current transmission process, uses the gradient inverse proportional weight distribution to map elevation changes, and triggers local topology reconstruction and boundary stitching. Subsequently, the system performs inverse distance weighted interpolation hot-start mapping of the dual scale vector, resumes iteration until convergence, outputs the global optimal allocation matrix, and maps it inversely to a directed earthwork allocation vector flow field. Finally, it performs Laplace grid smoothing operation subject to hard constraints of engineering design elevation and maximum allowable longitudinal slope, and outputs construction-level engineering data.

[0058] To verify the computational performance of the technical solution of the present invention under different terrain complexities, a traditional unconstrained linear programming allocation algorithm was used as a control group, and the specific operating indicators of each test area were recorded and statistically analyzed.

[0059] Table 1 Comparison of experimental data between the traditional allocation algorithm and the proposed technical solution in different test areas. Basic grid dimensions (horizontal × vertical) 125×140 215×185 88×95 Maximum permissible longitudinal slope setting value 22% 18% 26% Extract the number of topology nodes (source nodes / sink nodes) 17 / 14 31 / 28 11 / 13 Control group: Infeasibility rate of allocation path 19.34% 27.62% 15.11% This plan: Allocation path infeasibility rate 0.00% 0.00% 0.00% Control group: Global convergence time 16.27s 34.51s 9.44s This solution: Global convergence time 4.82s 8.19s 2.53s Control group: Number of grid cells crossing elevation boundaries 412 893 207 This plan: Number of grid cells exceeding elevation boundaries 0 0 0 Control group: Theoretical total earthwork transfer volume 92451.3m³ 154320.8m³ 45211.7m³ This plan: Actual total earthwork allocation 95103.8m³ 161874.2m³ 46835.4m³ Based on the data in Table 1 and Figure 7 The land consolidation engineering data analysis and optimization method provided by this invention, which takes into account topographic constraints, demonstrates definite technical effects in three core dimensions: allocation feasibility, calculation efficiency, and engineering benchmark compliance.

[0060] First, regarding the infeasibility rate of the allocated paths, the control group showed that 19.34%, 27.62%, and 15.11% of the paths in regions A, B, and C exceeded the mechanical climbing limit, respectively. In contrast, the infeasibility rate of the allocated paths in this scheme was reduced to 0.00% in all test areas. The underlying technical logic lies in the fact that when constructing the physical constraint cost matrix, this scheme uses a piecewise nonlinear penalty function to forcibly map the path weights exceeding the maximum allowable longitudinal slope (e.g., 18% in region B) to a maximum constant. Table 1 shows that the actual total earthwork allocated by this scheme increased by a reasonable margin compared to the control group; for example, in region A, the volume increased from 92451.3 m³. 3 Increased to 95103.8m 3 The volume increment reflects the mathematical cost of the algorithm's spatial detour to avoid physically inaccessible steep slopes, confirming that the physical resistance weight model established in this scheme can change the optimal transmission flow direction at the bottom layer and output a dispatching scheme that conforms to the actual operating capabilities of tracked machinery on site.

[0061] Secondly, regarding the global convergence time of the algorithm, the computational efficiency of this scheme is significantly higher than that of the control group. Taking region B with the largest grid dimension and the most complex terrain as an example, the global convergence time is reduced from 34.51 seconds to 8.19 seconds. The control group needs to reset global parameters and re-optimize when dealing with terrain deformation due to the lack of a topology reconstruction mechanism. This scheme monitors the cut and fill status in real time by calibrating the critical earthwork volume constant, performs boundary stitching when there are local grid elevation variations, and uses inverse distance weighted interpolation to perform hot-start mapping on the dual scale vector of unchanged nodes. The data shows that this hot-start mapping mechanism allows the system to avoid clearing the global iteration variable to zero when the matrix dimension changes due to the addition or removal of topology nodes, directly inherits the convergence state of the previous calculation, and eliminates redundant internal alternating projection loops.

[0062] Finally, the control group generated 412, 893, and 207 elevation boundary violations in each region during the smoothing output stage, respectively. This indicates that conventional grid smoothing disrupted the original engineering design elevation when smoothing out abrupt changes caused by discrete adjustments. In contrast, the number of elevation boundary violations in this scheme remained at 0 in all regions. This test result corresponds to the technical characteristics of the restricted Laplace grid smoothing operation in this scheme, namely, the algorithm calculates the dynamic upper and lower limits of node elevations based on the engineering design elevation in each smoothing iteration, and performs rigid boundary truncation on the intermediate elevation of the standard smoothing. The data verifies that this rigid truncation mechanism eliminates the risk of terrain deviation from the design baseline due to excessive numerical smoothing, ensuring that the final output construction-level engineering data is strictly controlled within the boundaries of engineering parameters.

[0063] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for data analysis and optimization of land consolidation projects that takes into account topographic constraints, characterized in that, Includes the following steps: Obtain the original digital elevation model matrix and engineering constraint parameters of the target area, and use the finite difference method to extract the low-level differential calculation features of the original digital elevation model matrix to calculate the micro-geometric features of the grid. Based on the discrete Morse theory of continuous homology, feature extraction and denoising are performed on the original digital elevation model matrix to construct a terrain feature topology map. Using the control area of ​​the retained critical point as the integration domain, surface integral is performed on the difference between the actual elevation and the target flat field elevation to obtain the critical earthwork volume constant corresponding to each critical point. Calculate the physical resistance weights of the terrain feature topology map based on the engineering constraint parameters, construct the physical constraint cost matrix, and initialize the global supply and demand distribution and iteration variables; Construct an optimal transport model with entropy regularization and monitor the cumulative earthwork transfer volume of nodes. When the volume reaches the set proportion of the corresponding critical earthwork volume constant, constrained by the first-order elevation gradient field, the transfer volume is allocated as grid elevation increments according to the inverse gradient ratio. Reconstruct the local topological subgraph and stitch the suspended streamlines of the subgraph with the discontinuous streamlines of the unupdated region according to the Euclidean nearest neighbor. Unchanged nodes inherit the original dual scale vector. New nodes are initialized with inverse distance weighted interpolation using the dual scale vectors of adjacent unchanged nodes and then iterated until convergence. Output the global optimal allocation matrix. The globally optimal allocation matrix is ​​transmitted and reverse-mapped to generate a directed earthwork allocation vector flow field, and a restricted Laplace mesh smoothing operation is performed to output construction-level engineering data.

2. The land consolidation engineering data analysis and optimization method considering topographic constraints according to claim 1, characterized in that, The steps of obtaining the original digital elevation model matrix and engineering constraint parameters of the target area, and calculating the micro-geometric features of the grid, include: Obtain the engineering constraint parameters, which include the set of polygonal boundaries of the engineering red line, the maximum allowable longitudinal slope, the maximum single cut and fill depth, and the physical operation error threshold. For the effective pixels in the original digital elevation model matrix located inside the polygon boundary set of the project red line and after filtering out missing outliers, the partial derivatives in the horizontal and vertical directions are calculated using the finite difference method to synthesize the first-order elevation gradient field and geometric slope. A second-order Hessian matrix is ​​constructed using the elevation information of adjacent nodes. The eigenvalues ​​are solved to obtain the average curvature and Gaussian curvature. Combined with the first-order elevation gradient field and geometric slope, the micro-geometric features of the grid are formed.

3. The land consolidation engineering data analysis and optimization method considering terrain constraints according to claim 2, characterized in that, The steps of extracting features and denoising the original digital elevation model matrix based on the discrete Morse theory of continuous homology to construct a terrain feature topology map include: Based on the micro-geometric features of the grid, local maxima, local minima, and saddle points of the terrain are extracted; Calculate the elevation difference between each extreme point and the saddle point associated with the extreme point, and merge the extreme points with elevation differences lower than the physical operation error threshold and the topological streamlines associated with the extreme points into the adjacent deep water catchment area or the area associated with the hillside and ridge to perform topological noise reduction. The effective critical points and streamlines after denoising are retained to form the topological map of the terrain features, and the retained local maxima are defined as earthwork supply source nodes, and the retained local minima are defined as earthwork demand sink nodes.

4. The land consolidation engineering data analysis and optimization method considering topographic constraints according to claim 1, characterized in that, The step of obtaining the critical earthwork volume constant corresponding to each critical point through surface integral includes: The local micro-topography boundary controlled by each retained critical point in the topographic feature map is used as the integral solution domain, wherein the boundary of the integral solution domain is determined by the closed surface region enclosed by the topological streamlines connecting the critical point and the adjacent associated saddle point of the critical point. By performing surface integrals on the difference between the actual surface elevation within the integral solution domain and the local target flat field elevation set to eliminate the topological features of the critical point, the critical earthwork volume constant reflecting the maximum amount of engineering work required to completely flatten a single local micro-topography is obtained.

5. The land consolidation engineering data analysis and optimization method considering terrain constraints according to claim 2, characterized in that, The step of calculating the physical resistance weights of the terrain feature topology map based on the engineering constraint parameters and constructing the physical constraint cost matrix includes: Extract the spatial physical paths of the topologically connected edges in the terrain feature topology map corresponding to the original digital elevation model matrix; Combining the physical horizontal distance of the spatial physical path with the average longitudinal slope, a piecewise nonlinear penalty function is set based on the maximum permissible longitudinal slope: When the average longitudinal slope of the space is less than or equal to the maximum allowable longitudinal slope, the physical resistance weight is calculated by combining the physical horizontal distance with the square of the terrain sensitivity penalty coefficient and the slope ratio. When the average longitudinal slope of the space exceeds the maximum allowable longitudinal slope, the physical resistance weight is forcibly mapped to a maximum penalty constant representing the physical inaccessibility state, so as to cut off the infeasible earthwork allocation path. The physical constraint cost matrix is ​​composed of the set of physical resistance weights between all nodes.

6. The land consolidation engineering data analysis and optimization method considering terrain constraints according to claim 4, characterized in that, The steps for initializing the global supply and demand distribution and the dual scale vector include: Based on the three-dimensional mesh, a discrete volume summation operation is performed to calculate the solid volume of each earthwork supply source node area. The volume is then multiplied by the loose or compacted earthwork coefficient to obtain the effective virtual volume. The volume is then constructed into a global supply distribution vector according to the node index order. Calculate the volume gap between the actual ground elevation and the local target leveling elevation within the control area of ​​each earthwork demand sink node, and generate a global demand distribution vector; Generate a first dual scale vector with the same dimension as the global supply distribution vector, and a second dual scale vector with the same dimension as the global demand distribution vector, and set the initial values ​​of all internal elements to a constant 1.

7. The land consolidation engineering data analysis and optimization method considering topographic constraints according to claim 1, characterized in that, The steps of constraining with a first-order elevation gradient field, distributing the transfer amount as grid elevation increments according to the inverse ratio of the gradient, reconstructing the local topological subgraph, and stitching the suspended streamlines of the subgraph with the discontinuities of the unupdated region according to the Euclidean nearest neighbor method include: Extract the local micro-topographic region controlled by the node where the change occurred, and use the first-order elevation gradient field as a constraint benchmark to map the accumulated earthwork transfer volume into the elevation change increment of each discrete grid according to the gradient inverse weight distribution, thereby correcting the actual surface elevation. For the local area after elevation update, a local terrain feature topology sub-graph is regenerated. The suspended topological streamlines at the edge of the local terrain feature topology sub-graph are searched. Based on the Euclidean nearest neighbor principle of spatial coordinates, the undirected edge pointers of the suspended topological streamlines are re-bound and connected to the breakpoint streamline nodes at the boundary of the unupdated area to complete the boundary stitching.

8. The land consolidation engineering data analysis and optimization method considering terrain constraints according to claim 7, characterized in that, The steps for the unchanged nodes to inherit the original dual scale vector, the new nodes to be initialized with inverse distance weighted interpolation using the dual scale vectors of adjacent unchanged nodes, the iteration to convergence, and the output of the globally optimal allocation matrix include: The physical resistance weights of the mutated connected edges are re-extracted and the local physical constraint cost matrix of the updated region is reconstructed. The supply and demand vectors of the changed nodes are recalculated to align with the global mathematical dimension. Perform a warm-start mapping of the dual scale vector space: For existing nodes that have not changed, directly inherit the dual scale vector of the previous state before the process was suspended. For the new node generated by the local update, locate all the unchanged nodes that are first-order adjacent to the new node on the stitched topological streamline, extract the dual scale vector of the adjacent node, and assign initial values ​​using the modified inverse distance weighted interpolation formula that introduces a preset minimum positive smoothing constant. Starting with the updated and recombined set of dual scale vectors, the optimal transmission iteration is restored until the supply and demand edge residuals meet the preset approximation accuracy conditions or the internal loop reaches the preset maximum safe iteration number, and the globally optimal allocation matrix is ​​output.

9. The land consolidation engineering data analysis and optimization method considering terrain constraints according to claim 3, characterized in that, The step of transmitting and reverse mapping the global optimal allocation matrix to generate a directed earthwork allocation vector flow field includes: Extract the sequence of non-zero elements from the global optimal allocation matrix; take the three-dimensional physical coordinates of the earthwork supply source node as the starting point, the three-dimensional physical coordinates of the corresponding earthwork demand sink node as the ending point, and the connected topological streamlines as the spatial skeleton, and assign the values ​​of the sequence of non-zero elements to the spatial skeleton as physical transfer quantities. The spatial transfer vector direction from the source node to the sink node is strictly calibrated, and then mapped back to the three-dimensional geometric space to generate the directional earthwork allocation vector flow field.

10. The land consolidation engineering data analysis and optimization method considering terrain constraints according to claim 2, characterized in that, The steps of performing restricted Laplacian mesh smoothing calculations and outputting construction-level engineering data include: The spatial weight parameters of each adjacent node are calculated using the inverse distance weighting method. Normalization is performed before the parameters are substituted into the iterative calculation, and the standard smooth intermediate elevation is calculated. Based on the engineering design elevation, and combined with the product of the elevation of adjacent nodes and the maximum allowable longitudinal slope over the grid physical resolution span, real-time calculations are performed to deduce and determine the dynamic upper and lower limits of node elevations. The standard smooth intermediate elevation is rigidly truncated using the dynamic upper and lower limits of the node elevation, and the updated grid elevation is output. The updated grid elevation and plane coordinate parameters are extracted and rendered to generate a three-dimensional grid elevation map. The discrete volumes of each grid region are summarized and compiled to generate a detailed cut and fill table. Combined with the directed earthwork allocation vector flow field decoding, the engineering machinery driving path matrix is ​​output. The above data is encapsulated into the construction-level engineering data and presented.

Images