High-fidelity terrain surface interpolation method and system for large-scale point cloud data

Abstract

The invention relates to the technical field of computer information processing, and discloses a high-fidelity terrain surface interpolation method and system for large-scale point cloud data, and the method comprises the steps: carrying out the multi-level cleaning of original point cloud data, and generating a cleaned point cloud data set; performing adaptive tile grid division based on the cleaned data to generate a tile set covering the target area; in each tile, taking the to-be-interpolated grid point as a query point, constructing and optimizing a local Gaussian process model, and obtaining a real elevation predicted value of the query point; processing all query points in the tiles in parallel to obtain a tile elevation interpolation result set; performing weighted average on each tile result to generate a global seamless terrain surface; according to the method, the problems of poor expandability and surface discontinuity caused by local interpolation in large-scale point cloud data processing are effectively solved, and high-precision and high-fidelity terrain surface construction is realized.

Description

Technical Field

[0001] This invention relates to the field of computer information processing technology, and more specifically, to a high-fidelity terrain surface interpolation method and system for large-scale point cloud data. Background Technology

[0002] Efficient processing and high-precision modeling of large-scale 3D point cloud data are core components of building digital elevation models (DEMs) and digital surface models (DSMs). With the development of technologies such as LiDAR and UAV photogrammetry, the scale of point cloud data has increased from millions of points to billions or even trillions of points. This has placed three demands on technologies such as spatial indexing, distributed computing, and interpolation modeling: "high scalability, high interpolation accuracy, and high surface continuity." Although existing related technologies have explored aspects such as tile processing and data distribution, a complete solution adapted to large-scale point cloud elevation interpolation has not yet been formed.

[0003] In the prior art, Chinese patent application CN119645842A discloses an automated testing method, apparatus, device, and medium for map operation response. It proposes to establish an index relationship between tile map location and serial number, simulate map operation, and compare the actual and target tile serial numbers to determine the test results. Its core is to use tile indexing to achieve accurate verification of map operation response. It only focuses on map test scenarios and does not involve spatial neighborhood query and elevation interpolation logic of point cloud data. Chinese patent application CN116226299A discloses a distributed data processing method, apparatus, device, and system. It divides map source data into tile grids, extracts adjacent edge data of target tiles and expands them to adjacent tiles to achieve deduplication, and then updates the map through distributed computing. Its technical focus is on distributed deduplication of tile data and map update, without considering the statistical distribution characteristics of point cloud elevation data and the boundary continuity of interpolation results.

[0004] Existing technologies cannot simultaneously solve several core problems when adapting to high-fidelity terrain surface interpolation scenarios using large-scale point cloud data, such as "scalability of large-scale point cloud processing, interpolation accuracy of non-Gaussian data, and continuity of tile boundaries," making it difficult to support the construction of high-fidelity terrain surfaces. Specifically, the tile index only serves map operation testing and cannot provide spatial index support for "quickly locating local neighboring points" for point cloud interpolation. Point cloud interpolation requires building a model based on local neighborhoods. The lack of this support will increase the neighborhood query time for billions of point clouds from microseconds to seconds, resulting in a loss of real-time processing capabilities. Adjacent edge data processing is only used for data deduplication and no fusion mechanism for tile boundary elevation values ​​is designed. If the tile division method is directly used, it will cause elevation abrupt changes of several meters or even tens of meters at the boundary of adjacent tile interpolation results (i.e., "tile effect"), making the terrain surface visually present a distinct facet, which cannot meet the requirements of terrain continuity in scenarios such as hydrological analysis and road design. No adaptation scheme is designed for the non-Gaussian distribution characteristics of point cloud data (such as data concentration in plain areas and long-tail distribution of data in mountainous areas). If the traditional interpolation model based on Gaussian assumptions is used, it will lead to a large fitting deviation of mountain micro-topographic features and cannot quantify the interpolation uncertainty, further reducing the fidelity of the terrain surface. Summary of the Invention

[0005] To overcome the aforementioned shortcomings of existing technologies, this invention provides a high-fidelity terrain surface interpolation method and system for large-scale point cloud data. Through a complete workflow design of "multi-level point cloud cleaning - K-dimensional tree neighborhood search - local Gaussian process modeling (including quantile-quantile normalization) - tile weighted fusion," it not only solves the "tile effect" caused by existing tile technology but also adapts to the non-Gaussian distribution characteristics of terrain elevation data, significantly improving the interpolation accuracy and terrain surface continuity of large-scale point cloud data. The parallelized tile processing architecture distributes single-tile computation tasks to independent computing units, greatly reducing the processing time of point clouds, while ensuring the continuity of the global terrain surface through a weighted fusion mechanism.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] High-fidelity terrain surface interpolation methods for large-scale point cloud data include:

[0008] Input the original point cloud dataset, perform multi-level cleaning on the original point cloud dataset to generate a cleaned point cloud dataset; based on the cleaned point cloud dataset, perform tile mesh division to generate a set of tiles covering the entire target area, and generate a set of grid points to be interpolated in each tile of the tile set;

[0009] Using the interpolation grid points in the set of interpolation grid points as query points, a local Gaussian process model is constructed and optimized. Based on the optimized local Gaussian process model, the predicted true elevation of the query point in its respective tile is obtained.

[0010] All query points within a tile are processed in parallel to obtain the true elevation prediction values ​​of all query points within the tile. Based on the true elevation prediction values ​​of all query points within the tile, a set of tile elevation interpolation results is obtained. The weighted average of each set of tile elevation interpolation results is then performed to obtain a global seamless terrain surface.

[0011] Each point in the original point cloud dataset contains three-dimensional coordinate information, which refers to the horizontal position coordinates X and Y and the elevation value Z.

[0012] The multi-level cleaning method includes:

[0013] Construct a hash mapping table based on horizontal position, and divide the 3D coordinate points in the original point cloud dataset into multiple horizontal position groups according to the hash mapping table; each horizontal position group has a single elevation value or multiple elevation values.

[0014] Statistical analysis was performed on horizontal location groups containing multiple elevation values, and outliers were removed to obtain a set of valid elevation values.

[0015] Aggregate the valid elevation values ​​within each horizontal location group to generate a cleaned point cloud dataset with a unique elevation value for each horizontal location.

[0016] The method for dividing the horizontal position groups includes:

[0017] Traverse the original point cloud dataset, extract the horizontal position coordinates (X,Y) of each point in the original point cloud dataset to generate a unique key value, and store the corresponding elevation value Z in the list under the key value to form a horizontal position group.

[0018] The method for constructing and optimizing a local Gaussian process model includes:

[0019] Obtain the local neighborhood point set of the query point, construct the original elevation value set based on the local neighborhood point set, transform the original elevation value set into a normal elevation value set through quantile-quantile normalization transformation, and construct and store the transformation mapping table at the same time;

[0020] Based on the normalized set of elevation values, a local Gaussian process model based on the Matérn kernel function is constructed. The hyperparameters of the Matérn kernel function are optimized by maximum likelihood estimation to obtain the optimized local Gaussian process model.

[0021] The method for obtaining the local neighborhood point set of the query point includes:

[0022] Use a median-based recursive partitioning strategy to process the cleaned point cloud dataset and construct a K-dimensional tree spatial index structure;

[0023] Use the K-dimensional tree spatial index structure to perform an adaptive hybrid neighborhood search on any query point to obtain the local neighborhood point set of the query point.

[0024] The method of using the K-dimensional tree spatial index structure to perform an adaptive hybrid neighborhood search on any query point includes:

[0025] Determine the initial radius and density threshold. Using the K-dimensional tree spatial index structure, perform an initial radius search centered on the query point and count the number of points Nr within the initial radius;

[0026] The density threshold includes a density threshold lower limit kmin and a density threshold upper limit kmax;

[0027] When kmin ≤ Nr ≤ kmax, take all Nr points within the initial radius r' as neighborhood points to form the local neighborhood point set of the query point;

[0028] When Nr < kmin, obtain the kmin points closest to the query point as the neighborhood points of the query point and construct the local neighborhood point set of the query point;

[0029] When Nr > kmax, select the kmax points closest to the query point as the neighborhood points of the query point to form the local neighborhood point set of the query point.

[0030] The method of transforming the original elevation value set into a normalized elevation value set includes:

[0031] Sort the original elevation values in the original elevation value set in ascending order to obtain an ordered elevation value sequence;

[0032] Calculate the empirical cumulative probability for each original elevation value in the ordered elevation value sequence and apply the inverse cumulative distribution function of the standard normal distribution to generate a normalized elevation value set.

[0033] The method of optimizing the hyperparameters of the Matérn kernel function by maximum likelihood estimation includes:

[0034] Determine the objective function for hyperparameter optimization. The objective function is the log marginal likelihood function, which includes three terms. The first term is the fitting error term, the second term is the regularization term of the initial kernel matrix, and the third term is a constant term; the fitting error term is constructed based on the normalized elevation value set;

[0035] Based on the objective function, use the gradient ascent algorithm to iteratively optimize the hyperparameters to obtain an optimized hyperparameter set.

[0036] The method for obtaining the predicted true elevation value of the query point in its respective tile area includes:

[0037] Based on the optimized local Gaussian process model, the mean normalized elevation prediction value of the query point is predicted.

[0038] The normalized elevation prediction mean of the query point is transformed by inverse quantile transformation through the stored transformation mapping table to obtain the true elevation prediction value of the query point in the tile to which it belongs.

[0039] A high-fidelity terrain surface interpolation system for large-scale point cloud data, which implements the above-mentioned high-fidelity terrain surface interpolation method for large-scale point cloud data, the system comprising:

[0040] Point cloud data cleaning module: Input the raw point cloud dataset, perform multi-level cleaning on the raw point cloud dataset, and generate a cleaned point cloud dataset;

[0041] Tile partitioning module: Used to partition the tile grid based on the cleaned point cloud dataset, generate a set of tiles covering the entire target area, and generate a set of grid points to be interpolated within each tile in the tile set;

[0042] Elevation prediction module: The interpolation grid points in the set of interpolation grid points are used as query points. A local Gaussian process model is constructed and optimized. Based on the optimized local Gaussian process model, the true elevation prediction value of the query point in the tile to which it belongs is obtained.

[0043] Terrain surface generation module: It is used to process all query points within a tile in parallel, obtain the true elevation prediction value of all query points within the tile, obtain the tile elevation interpolation result set based on the true elevation prediction value of all query points within the tile, and perform a weighted average on the tile elevation interpolation result set to obtain the global seamless terrain surface.

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

[0045] This invention effectively removes data redundancy, noise, and elevation value conflicts at the same horizontal position by performing multi-level cleaning on the original point cloud dataset, solving the problem of poor quality of the original point cloud data. This provides a high-quality and consistent data foundation for subsequent interpolation, ensuring the reliability of elevation prediction from the source. Based on the tile grid division of the cleaned data, large-scale point cloud data can be reasonably decomposed into several independent processing units. Combined with the parallel processing of all query points within the tile, the processing efficiency of massive point cloud data is greatly improved, solving the problems of weak scalability and computational time consumption when traditional interpolation methods process large-scale data. By weighted averaging of the elevation interpolation results of each tile, the boundary discontinuity phenomenon that is prone to occur in local interpolation can be effectively eliminated. At the same time, combined with the construction and optimization of the local Gaussian process model, the accuracy of the elevation prediction of the query point is ensured. Ultimately, it achieves the generation of a high-fidelity global seamless terrain surface while efficiently processing large-scale point cloud data, meeting the core requirements of data processing efficiency, interpolation accuracy, and terrain surface continuity in the construction of digital elevation models or digital surface models. Attached Figure Description

[0046] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. 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.

[0047] Figure 1 This is a flowchart of a high-fidelity terrain surface interpolation method for large-scale point cloud data provided in an embodiment of the present invention;

[0048] Figure 2 This is a flowchart of a method for multi-level cleaning of raw point cloud datasets provided in an embodiment of the present invention;

[0049] Figure 3 This is a schematic diagram of an overlapping area of ​​adjacent tiles provided in an embodiment of the present invention;

[0050] Figure 4 This is a schematic diagram illustrating the principle of obtaining a local neighborhood point set of a query point according to an embodiment of the present invention.

[0051] Figure 5 A flowchart of a method for transforming an original set of elevation values ​​into a normalized set of elevation values ​​is provided in an embodiment of the present invention;

[0052] Figure 6 A scatter plot comparing the interpolation results provided in this embodiment of the invention with the actual elevation;

[0053] Figure 7The original data heatmap, interpolation result heatmap, and error heatmap provided for embodiments of the present invention;

[0054] Figure 8 The histogram of prediction error distribution provided in the embodiments of the present invention;

[0055] Figure 9 The image shows the application effect of the interpolation method of the present invention in a mountainous terrain scene according to an embodiment of the present invention.

[0056] Figure 10 The illustration shows the application effect of the interpolation method of the present invention in a plain and hilly terrain scene provided in the embodiment of the present invention.

[0057] Figure 11 This is a functional block diagram of a high-fidelity terrain surface interpolation system for large-scale point cloud data provided in an embodiment of the present invention. Detailed Implementation

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

[0059] Example 1

[0060] Please see Figure 1 As shown, this embodiment provides a high-fidelity terrain surface interpolation method for large-scale point cloud data, including:

[0061] Step S10: Input the original point cloud dataset, perform multi-level cleaning on the original point cloud dataset to generate a cleaned point cloud dataset; based on the cleaned point cloud dataset, perform tile grid division to generate a tile set covering the entire target area, and generate a set of grid points to be interpolated in each tile of the tile set.

[0062] Step S10, through a progressive process of "multi-level cleaning - spatial index construction - adaptive tile partitioning", solves the problems of "poor data quality, slow neighborhood query, and unbalanced computational allocation" in the processing of massive point cloud data, laying the foundation for subsequent high-precision interpolation and seamless integration.

[0063] Further, step S10 includes:

[0064] Step S11: Input the original point cloud dataset, perform multi-level cleaning on the original point cloud dataset, and generate a cleaned point cloud dataset.

[0065] Further, see Figure 2 Step S11 includes:

[0066] Step S111: Construct a hash mapping table based on horizontal position, and divide the 3D coordinate points in the original point cloud dataset into multiple horizontal position groups according to the hash mapping table; each horizontal position group has a single elevation value or multiple elevation values.

[0067] Step S112: Perform statistical analysis on horizontal position groups containing multiple elevation values ​​and remove outliers to obtain a set of valid elevation values;

[0068] Step S113: Aggregate the set of valid elevation values ​​within each horizontal position group to generate a cleaned point cloud dataset with a unique elevation value for each horizontal position.

[0069] Specifically, the original point cloud dataset is used to construct a continuous digital terrain surface. It consists of massive point cloud data acquired through methods such as LiDAR, UAV photogrammetry, or InSAR. Each point in the original point cloud dataset contains three-dimensional coordinate information, namely horizontal position coordinates X and Y and an elevation value Z. In step S11, the hash map table is a key-value pair storage structure built by using the horizontal position coordinates (X, Y) of the three-dimensional coordinate points as keys and the corresponding elevation value Z as values. Its core function is to achieve rapid aggregation of points with the same horizontal position, solving the problem of drastic efficiency decline with increasing data volume during traditional traversal grouping. When constructing the hash map table, the original point cloud dataset is traversed, the horizontal position coordinates (X, Y) of each point in the original point cloud dataset are extracted to generate a unique key value, and the corresponding elevation value Z is stored in the list under the key value, forming horizontal position groupings. The reason for using a hash map table is that its insertion and lookup operations have a time complexity close to constant level. Even with billions of point clouds, it can still quickly complete grouping, significantly improving processing efficiency compared to array traversal grouping. For example, the original point cloud dataset contains several points, where the coordinates (X1,Y1) correspond to Z values ​​Z1-1, Z1-2, and Z1-3, and (X2,Y2) corresponds to Z value Z2-1. Through a hash mapping table, the three Z values ​​of (X1,Y1) can be directly grouped into the same group, and the single Z value of (X2,Y2) can be grouped independently, without the need to compare all coordinates point by point.

[0070] If a horizontal location group contains only a single elevation value (i.e., a single point), that elevation value is directly retained as the representative elevation value for that horizontal location. For statistical outlier removal in horizontal location groups containing multiple elevation values, the median Z-mean is used. medianThe median is chosen as a combined standard with three times the standard deviation σ, rather than the traditional mean standard, because the mean is easily affected by extreme outliers, leading to a shift in the statistical baseline. The median, however, is minimally affected by extreme values ​​and more stably reflects the central tendency of elevation values ​​within a group. In the calculation, the elevation values ​​Z (i.e., Z-values) within a group containing multiple elevation values ​​are first sorted in ascending order, and the value at the median position is taken as the median Z. median Then, the standard deviation σ is calculated using the standard deviation formula. The standard deviation formula measures the difference between each Z value and Z0. median The degree of deviation is used to quantify the data dispersion. An outlier detection threshold of 3 times the standard deviation is set because of the normal distribution characteristic; the elevation values ​​corresponding to normal measurement errors are mostly concentrated in the Z-axis. median Points outside the ±3σ range are likely noise points. After removing outliers, the remaining valid elevation values ​​(i.e., valid Z-values) constitute the set of valid elevation values. Aggregating the set of valid elevation values ​​within each horizontal position group refers to aggregating the remaining valid elevation values ​​Z-values. x Calculate the arithmetic mean as the representative elevation value Z of the horizontal position corresponding to the horizontal position grouping. d The arithmetic mean is used to average and integrate multiple reliable measurements, i.e., effective elevation values, reducing the impact of single measurement errors. This ensures that each horizontal location retains only a unique and reliable elevation value, i.e., a representative elevation value. All horizontal locations are grouped together to form a unique horizontal elevation. The location coordinates (X,Y) and the representative elevation value of the group together constitute the cleaned point cloud dataset.

[0071] The multi-level cleaning process in step S11 addresses three core issues in the original point cloud dataset: spatial grouping and deduplication eliminates duplicate points with the same (X,Y) coordinates, reducing data redundancy and significantly decreasing the amount of data processed subsequently; statistical anomaly removal eliminates noise points caused by multipath reflection and atmospheric disturbances, improving data purity; and effective data aggregation resolves conflicts between multiple elevation values ​​at the same location, ensuring data consistency. Without step S11, noise and duplicate points in the original point cloud dataset would directly enter subsequent stages, leading to redundancy in the K-dimensional tree index nodes, decreased query efficiency, and severe fitting bias in the local Gaussian process model due to unreliable input data, compromising the accuracy of the entire solution.

[0072] Step S12: The cleaned point cloud dataset is processed using a median-based recursive segmentation strategy to construct a K-dimensional tree space index structure.

[0073] The K-dimensional tree spatial index is a hierarchical spatial partitioning structure built on two-dimensional horizontal coordinates (X,Y). The value of K is 2, corresponding to the two dimensions of the horizontal coordinates. Its core function is to transform unordered scattered points in a cleaned point cloud dataset into structured data that can be retrieved quickly, thus overcoming the efficiency bottleneck of traditional brute-force search. The K-dimensional tree spatial index structure is constructed using the cleaned point cloud dataset as its original data foundation, specifically employing a recursive partitioning strategy based on the median: at the beginning of construction, the root node directly contains all three-dimensional coordinate points in the cleaned point cloud dataset, meaning that the initial data of the root node comes entirely from the cleaned point cloud dataset, serving as a complete representation of the cleaned point cloud dataset. Starting from the root node, alternately select X and Y dimensions for segmentation: For the current node (initially the root node, subsequently child nodes), all points it contains come from a subset of the cleaned point cloud dataset. First, calculate the median coordinates of all points within the current node along the selected dimension. Using this median as the boundary, divide the point set within the current node into left and right subsets. Points with coordinate values ​​less than or equal to the median are assigned to the left subtree, and points with coordinate values ​​greater than the median are assigned to the right subtree. Repeat this segmentation process, with each segmentation involving a subset of the cleaned point cloud dataset (subsets derived from the parent node), until the termination condition is met. The termination condition is: if the number of points in the current node is less than a preset node count threshold, for example, 10 points, then stop segmentation and treat it as a leaf node; when the tree depth reaches log2(N)... x When N is reached, forced termination is performed to prevent the tree from degenerating into a linked list in extreme cases, thus ensuring the balance of the tree. x This represents the total number of points in the cleaned point cloud dataset. Throughout the construction process, the cleaned point cloud dataset serves as the "source" of the K-tree: the initial data of the root node and the segmentation objects of each layer's child nodes are all points or subsets of the cleaned point cloud dataset. This ensures that the K-tree index uses cleaned, high-quality data, laying the foundation for efficient neighborhood queries in the future.

[0074] The method described above for constructing a K-dimensional tree spatial index structure ensures that the number of nodes in the left and right subtrees of the K-dimensional tree is approximately equal. During a query, the search range can be quickly narrowed down by comparing dimensions, reducing the query time complexity from O(N) to O(N). x The value decreases to O(log2(N)). x O(N) x ) represents the execution time of the algorithm versus N. x It exhibits a linear growth relationship, O(log2(N) x )) indicates the execution time of the algorithm versus "N". xThe logarithm to base 2 exhibits a logarithmic growth relationship. This optimization is crucial for processing massive amounts of data. For example, a single query of a billion-level point cloud dataset requires traversing all points using brute force, taking seconds or even minutes. However, using a K-tree index can shorten this time to microseconds, meeting the needs of real-time interpolation. The cleaned data, because each (X,Y) coordinate contains only a unique Z value, avoids point conflicts at the same spatial location during K-tree construction, making the index structure more stable. Simultaneously, the reduced data volume in step S11 lowers the difficulty of building the K-tree and storage overhead, improving index construction efficiency. Without step S12, the subsequent neighborhood search in step S21 would not achieve rapid location. Even with an adaptive neighborhood strategy, the long query time would compromise real-time performance, rendering the scalability of the solution ineffective.

[0075] Step S13: Calculate the spatial extent and global point cloud density of the cleaned point cloud dataset. Based on the spatial extent and global point cloud density of the cleaned point cloud dataset, adaptively divide the tile grid to generate a set of tiles covering the entire target area. Generate a set of grid points to be interpolated within each tile.

[0076] The core of step S13 is to achieve adaptive tile mesh generation and interpolation point preset based on the spatial characteristics of the cleaned point cloud data and the hardware computing capabilities. First, the minimum bounding rectangle boundary (Xmin, Xmax, Ymin, Ymax) of the cleaned point cloud is calculated. Xmin refers to the minimum horizontal position coordinate X of all points in the cleaned point cloud dataset, representing the leftmost boundary of the target terrain area in the X direction; Xmax refers to the maximum X coordinate of all points, representing the rightmost boundary in the X direction; Ymin refers to the minimum Y coordinate of all points, representing the bottommost boundary in the Y direction; Ymax refers to the maximum Y coordinate of all points, representing the topmost boundary in the Y direction. These four values ​​together constitute the minimum bounding rectangle boundary of the cleaned point cloud data, accurately defining the spatial range of the target area and avoiding incomplete coverage or redundancy in subsequent tile generation.

[0077] The total number of points N in the cleaned point cloud dataset x The global point cloud density ρ is calculated as the ratio of the area covered by the point cloud to the total area covered by the point cloud. total The total area covered by the point cloud is (Xmax-Xmin)×(Ymax-Ymin); the target number of points N for each tile is set based on the number of CPU cores P available for computing resources. w =N x / (4×P), the 4-fold factor is set to reserve computational redundancy and avoid excessive load on a single core. ρ totalA larger value indicates more points per unit area, providing a density basis for calculating the tile reference size and avoiding point imbalance caused by mismatch between tile size and point cloud density. The number of CPU cores P available for computing resources refers to the number of cores with independent computing capabilities in the central processing unit that can be called upon when executing this method. It directly reflects the hardware's parallel processing potential, such as the number of cores in a server-grade CPU or the number of cores in a desktop CPU. Its value needs to be determined based on the actual hardware configuration, providing a basis for matching the number of tiles with computing resources.

[0078] Determine the number of target points N for the tile w :N w =N x / (4×P), where the 4-fold coefficient is set based on the following: overlapping tile areas will generate data redundancy, and there are additional overheads such as cache swapping when a single CPU core processes tiles. Reserving 3 times the redundancy (total coefficient 4) allows the actual data volume processed by each core to approach its optimal capacity, avoiding core overload and timeout due to too many tile points, or core idle waste due to too few points. This solves the problem of traditional fixed target tile numbers not adapting to changes in hardware resources. Target tile number N w It is the expected number of point cloud data contained in a single tile. Its core function is to control the point cloud scale of a single tile to adapt to hardware computing resources (such as the number of CPU cores), so as to avoid overload caused by too many points processed by a single computing unit (core) or idle resources caused by too few points. In essence, it is a "control index of single-tile computing volume".

[0079] According to N w and ρ total Calculate the width W and height H of the tile using the following formula: ,in, The formula uses a square root; its logic is as follows: the target number of points on the tile is divided by the global point cloud density to obtain the required horizontal area of ​​the tile. Taking the square root yields the side length of a square tile. This ensures that tiles are small and have a controllable number of points in dense point cloud areas, while being large and have a sufficient number of points in sparse areas. This solves the problem of point imbalance in areas with uneven density caused by traditional fixed tile sizes. For example, ρ total When it is large, As the result decreases, W and H also decrease, maintaining the number of points within the tile close to the target value; ρ total When the tile size is smaller, the tile size is increased while still ensuring the number of dots matches. Then, the overlap area width M is set, initially set to 20% of W, such as... Figure 3 As shown, the shaded area between adjacent tiles 1 and 2 is the overlapping area, with a width of M. The overlapping area between adjacent tiles needs to be wide enough to ensure the transition space during fusion, avoid discontinuity caused by missing boundary data, and at the same time, the 20% ratio will not cause too much data redundancy due to excessive overlap, thus solving the boundary problem of traditional non-overlapping or improper overlap ratio.

[0080] Based on the determined tile width W, tile height H, and overlapping area width M, a tile grid is created. All tile grids form a set of tiles covering the entire target area. In the X direction, starting from Xmin, tile centers are set every (WM). The same logic applies to the Y direction. The logic is that the spacing between adjacent centers ensures that the overlap width between the right boundary of the previous tile and the left boundary of the next tile is exactly M, guaranteeing coverage of the entire target area without omissions, thus solving the coverage problem of traditional arrangement methods. Finally, a set of interpolation grid points is generated within each tile according to a preset interpolation resolution. The interpolation resolution is determined based on the terrain accuracy requirements; high resolution results in denser grid points, while low resolution results in sparser grid points. The generation method involves generating points sequentially within the X and Y ranges of the tile according to the resolution, clearly defining the interpolation target location and solving the problem of disordered computation in traditional methods. This invention, based on hardware and point cloud density partitioning, fully utilizes multi-core resources, avoiding overload or idleness; the overlapping area provides a transition for subsequent S32 fusion. The adaptive mechanism is robust to sudden changes in local point cloud density, avoiding a surge in the number of local tile points. Simultaneously, the preset grid points allow for pre-planning of memory, reducing memory fluctuations during computation and improving system stability.

[0081] Step S10 addresses the core challenges of massive point cloud preprocessing through a three-tiered collaborative approach: multi-level cleaning provides high-quality data input for subsequent steps, K-dimensional tree indexing offers an efficient tool for neighborhood queries, and adaptive tile partitioning provides a reasonable framework for parallel computing. The cleaned data improves the efficiency of K-dimensional tree index construction and querying; the index's efficiency allows for rapid neighborhood searches within tiles; and the parallelism of tile partitioning further amplifies the efficiency gains from index optimization, forming a positive cycle of "data purification - index acceleration - parallel computing." This preprocessing system fundamentally overcomes the bottlenecks of traditional methods in massive data processing—low efficiency, poor accuracy, and weak scalability—enabling efficient and accurate implementation of subsequent local Gaussian process modeling and global fusion.

[0082] Step S20: Take the interpolation grid points in the set of interpolation grid points as query points, construct and optimize the local Gaussian process model, and obtain the true elevation prediction value of the query point in the tile it belongs to based on the optimized local Gaussian process model.

[0083] The core of step S20 is to address the problems of fitting bias caused by the non-Gaussian distribution of real terrain elevation data, unstable neighborhood quality caused by uneven point cloud density, and meaningless output of Gaussian process model by using a progressive process of "adaptive neighborhood construction - non-Gaussian data normalization - local Gaussian process optimization - elevation prediction and inverse transformation" for the grid points to be interpolated (query points).

[0084] Further, step S20 includes:

[0085] Step S21: Take the interpolation grid points in the set of interpolation grid points as query points, and use the K-dimensional tree space index structure to perform adaptive hybrid neighborhood search on any query point to obtain the local neighborhood point set of the query point.

[0086] The method for obtaining the local neighborhood point set of the query point includes: determining the initial radius and density threshold, using a K-dimensional tree space index structure, performing an initial radius search with the query point as the center, counting the number of points within the initial radius, selecting an appropriate neighborhood construction strategy based on the relationship between the number of points within the initial radius and the density threshold, and constructing the local neighborhood point set of the query point.

[0087] Specifically, the core of step S21 is to achieve adaptive hybrid neighborhood search of query points based on the density features of the cleaned point cloud and the K-dimensional tree index, so as to solve the problem of unstable number of neighborhood points caused by uneven spatial distribution of point cloud density, and provide a reliable local neighborhood point set for subsequent construction of local Gaussian process model.

[0088] The initial radius r' refers to the radius of the initial circular search range centered on the query point. Its value needs to be adapted to the global point cloud density to ensure that a sufficient number of neighborhood points can be obtained in the average density region. The density threshold refers to the critical range for judging whether the number of neighborhood points is suitable for the modeling requirements. It includes the lower limit of the density threshold kmin and the upper limit of the density threshold kmat. kmin is the minimum number of neighborhood points required to build a local Gaussian process model, and kmat is the maximum number of neighborhood points to control the computational complexity, so as to avoid excessive matrix operation time due to too many points. The setting of both needs to be based on the trade-off between model complexity and computational efficiency. For example, kmin can be set to 15-20 and kmat can be set to 40-60. The specific values ​​need to be adjusted according to the actual modeling accuracy requirements and hardware computing capabilities.

[0089] After determining the initial radius and density threshold, the constructed K-dimensional tree spatial index structure is invoked. A circular region search with radius r' is performed, centered on the query point q(xq,yq). The hierarchical spatial partitioning structure of the K-dimensional tree can quickly locate all points within this region, avoiding the inefficiency of traditional brute-force search that traverses all points. After the search is complete, the number of points Nr within the region is calculated. Subsequently, based on the relationship between Nr and kmin and kmat, a neighborhood construction strategy is dynamically selected. Please refer to [link to relevant documentation]. Figure 4As shown, when kmin ≤ Nr ≤ kmax, it indicates that the point cloud density in the current area is moderate. The number of neighboring points searched can meet the requirements of model modeling without causing excessive computational complexity. Therefore, directly take all Nr points within the initial radius r' as neighboring points to form the local neighboring point set Nq of the query point q. The advantage of this strategy is to make full use of the density characteristics of the current area without additional screening, ensuring the spatial integrity of neighboring points. When Nr < kmin, it indicates that the point cloud in the current area is too sparse, and insufficient points can be obtained by searching according to r'. If the points in this area are continued to be used to build the model, it will lead to unreliable estimation of model parameters. Therefore, abandon the radius constraint and use K-nearest neighbor search instead. Quickly obtain the kmin points closest to the query point q as the neighboring points of the query point q through the K-dimensional tree, and build the local neighboring point set Nq of the query point q. This strategy forcibly ensures the number of neighboring points, avoiding fitting deviation of the model due to insufficient data. If the number of available point clouds in the entire domain is less than kmin, then take all available neighboring points. When Nr > kmax, it indicates that the point cloud in the current area is too dense. Too many neighboring points will increase the complexity of kernel matrix calculation and inversion, and the spatial correlation between distant points and the query point is relatively low, which has limited improvement on the modeling accuracy. Therefore, retain the spatial constraint of the radius r', and only select the kmax points closest to the query point q as the neighboring points of the query point q to form the local neighboring point set Nq of the query point q. This strategy improves the model fitting accuracy while controlling the computational amount by screening high-correlation points. The switching of the three dynamic strategies makes the present invention robust to sparse areas, dense areas, and moderate areas, avoiding the problems of modeling failure in sparse areas and computational overload in dense areas of traditional methods.

[0090] It should be noted that K in the K-dimensional tree spatial index structure is the number of spatial dimensions, based on the spatial coordinate dimensions of the point cloud. Since only two-dimensional horizontal coordinates (X, Y) are used, K = 2 in the K-dimensional tree spatial index structure; K in the K-nearest neighbor search refers to the "number of nearest neighbor points" to be obtained, based on the modeling requirements of the local Gaussian process model, that is, the lower density threshold kmin, to ensure that the model has sufficient data support.

[0091] The core problem addressed by step S21 is the inability of fixed neighborhood parameters in existing technologies to adapt to the uneven density space of point clouds. It achieves stable neighborhood point set quality through a "global density-based initial radius + local point number dynamic adjustment strategy." Step S21 relies on the K-dimensional tree spatial index structure from step S12 for efficient search. Without this structure, radius search and K-nearest neighbor search cannot be completed quickly, leading to excessively long neighborhood acquisition times. Step S21 provides a reliable local neighborhood point set Nq for step S22. If the number of points in Nq is unstable or contains low-relevance points, the normalization transformation data in step S22 will be unreliable, thus affecting the accuracy of subsequent models. The combination of K-nearest neighbor search and radius constraints ensures both the number of points and filters for highly relevant points, improving the model's ability to capture local terrain features. The strategy of filtering the nearest point can automatically exclude abnormal residual points in the local neighborhood, such as noise points that were not completely removed in step S11. Because abnormal points are usually far away from the query point, they are more likely to be filtered out, which indirectly improves the purity of the neighborhood point set. The dynamic strategy switching does not require manual intervention, enabling the method to have end-to-end processing capabilities and adapt to point cloud data of different terrain types.

[0092] Step S22: Extract the representative elevation values ​​of all neighboring points in the local neighborhood point set of the query point to form the original elevation value set. Transform the original elevation value set into a normalized elevation value set through quantile-quantile normalization transformation, and at the same time construct and store the transformation mapping table.

[0093] Further, see Figure 5 Step S22 includes:

[0094] Step S221: Extract the representative elevation value of each neighboring point from the local neighboring point set of the query point as the original elevation value to form an original elevation value set. Sort the original elevation values ​​in the original elevation value set in ascending order to obtain an ordered elevation value sequence.

[0095] Step S222: Calculate the empirical cumulative probability for each original elevation value in the ordered elevation value sequence, and apply the inverse cumulative distribution function of the standard normal distribution to generate a normalized elevation value set;

[0096] Step S223: Record the mapping relationship between the ordered elevation value sequence and the normalized elevation value set, and construct and store the transformation mapping table.

[0097] Specifically, the core of step S22 is to convert the non-Gaussian distribution elevation values ​​of the local neighborhood point set into a standard normal distribution through quantile-quantile normalization transformation, thereby solving the problem that the Gaussian process model's assumptions about data distribution do not match the actual terrain data. At the same time, a transformation mapping table is constructed to provide a basis for the subsequent inverse transformation.

[0098] From the local neighborhood point set Nq, the representative elevation value of each neighborhood point is extracted as the original elevation value. All extracted original elevation values ​​constitute the original elevation value set. The extraction operation must ensure that no point in Nq is missed to avoid data loss that could lead to subsequent transformation deviations. Then, all original elevation values ​​in the original elevation value set are sorted in ascending order to obtain an ordered elevation value sequence Zq={z (1) ,z (2) ,…,z (nq)}, where z (1) ≤z (2) ≤…≤z (nq) , i.e. z (1) It is the smallest original elevation value in the ordered elevation value sequence, z (2) This is the second smallest original elevation value, and so on, z (nq) As the last element of the ordered elevation value sequence, it corresponds to the largest original elevation value after sorting, where nq is the number of points in the local neighborhood point set Nq. The purpose of sorting is to assign a unique position index i' (i' ranges from 1 to nq) to each original elevation value, ensuring the orderliness and uniqueness of subsequent empirical cumulative probability calculations. Without sorting, the disordered positions of the original elevation values ​​would lead to errors in the cumulative probability calculation, thus causing the normalization transformation to fail.

[0099] Methods for calculating the empirical cumulative probability for each original elevation value in an ordered elevation value sequence include: for each original elevation value z in the ordered elevation value sequence... (i') Calculate the empirical cumulative probability p (i') Using formula p (i') =(i'-0.5) / nq, where "i'-0.5" is the continuity correction term. This is because the original elevation values ​​are discrete data, while the standard normal distribution is continuous. Continuity correction reduces the deviation of discrete data during the transformation process, making the transformed normalized elevation values ​​closer to the continuous standard normal distribution, thus avoiding distribution distortion caused by discreteness. The basis for taking the correction term as 0.5 is: for a neighborhood point set of nq≥15 (in this example, kmin≥15, therefore nq≥15), 0.5 conforms to the conventional choice of the statistical "midpoint correction method," balancing the deviation between discrete samples and continuous distributions; if nq<10 (small sample), it can be adjusted to 0.3~0.7, requiring dynamic adaptation based on the sample size. Each p... (i') Substituting the inverse cumulative distribution function Φ of the standard normal distribution -1 The corresponding normalized elevation value z' is calculated. (i') =Φ -1 (p (i') ), all z' (i') This forms a normalized set of elevation values ​​Z'q. Φ is chosen. -1The reason is that the theoretical basis of the Gaussian process model is based on the assumption that the data follows a normal distribution. This function can be used to represent p from any distribution. (i') Mapping to values ​​from a standard normal distribution ensures that the normalized elevation value set satisfies the model assumptions, thus resolving the fitting bias problem caused by directly using non-Gaussian data in traditional methods. Φ -1 The implementation can be achieved using commonly used numerical calculation methods in this field, such as Newton's iteration method and Chebyshev's approximation method. Those skilled in the art can choose the appropriate implementation method according to the requirements of calculation accuracy and efficiency.

[0100] Construct a transformation mapping table Mq for query point q, and map each original elevation value z in the ordered elevation value sequence Zq. (i') Its corresponding normalized elevation value z' (i') Form a mapping pair, i.e.: Mq={(z (1) ,z' (1) ),(z (2) ,z' (2) ),…,(z (nq) ,z' (nq) The process involves storing Mq locally or in a cache. The purpose of storing Mq is to provide a basis for the inverse quantile transformation in subsequent step S25. Furthermore, Mq must correspond one-to-one with each query point q because the distribution of the local neighborhood point set Nq differs for different query points, and the corresponding normalization transformation rules also differ. A dedicated Mq ensures the accuracy of the inverse transformation and avoids the precision loss caused by a globally unified mapping table.

[0101] The core problem addressed in step S22 is the contradiction between the non-Gaussian distribution of real terrain elevation data and the normal distribution assumption of the Gaussian process model. Through a technical approach of "sorting-cumulative probability calculation-normalization transformation-mapping table storage," the non-Gaussian data is converted to normal data, while simultaneously providing support for the inverse transformation. Step S22 relies on the high-quality local neighborhood point set Nq provided in step S21. If Nq contains noisy points or has insufficient points, the original elevation value set will have an abnormal distribution, resulting in poor normalization transformation. Step S22 provides the local Gaussian process model in step S23 with input data that conforms to the assumptions, namely, the normalized elevation value set, ensuring the model's fitting accuracy and the reliability of uncertainty quantification. The inverse transformation in step S25 depends entirely on Mq. If Mq is incorrectly constructed or not stored, the true elevation prediction value cannot be obtained.

[0102] Ascending order in the ordered elevation value sequence ensures the orderliness of cumulative probability calculation, avoiding transformation bias caused by positional confusion; continuity correction improves the fit between normalized data and the standard normal distribution, reducing fitting errors caused by data distribution mismatches; Φ -1This method achieves the transformation from arbitrary distributions to a standard normal distribution, enabling Gaussian process models to operate under optimal assumptions. The construction and storage of Mq ensures the accuracy of the inverse transformation, allowing prediction results to revert to the original elevation scale, thus possessing practical application value. Continuity correction not only improves normalization accuracy but also indirectly suppresses the influence of residual outliers in local neighborhoods. Outliers, after sorting, are typically located at the ends of the sequence, and their corresponding p... (i') After correction, the mapped z' (i') It will not deviate excessively from the reasonable range of the standard normal distribution, reducing the interference of outliers on the model; the unique Mq for each query point allows the transformation rules to adapt to the local data distribution. Compared with a globally uniform transformation, it can more accurately preserve the subtle features of the local terrain and improve the subsequent model's ability to fit terrain details.

[0103] Step S23: Based on the normalized elevation value set, construct a local Gaussian process model based on the Matérn kernel function, and optimize the hyperparameters of the Matérn kernel function through maximum likelihood estimation to obtain the optimized local Gaussian process model.

[0104] Further, step S23 includes:

[0105] Step S231: Based on the normalized elevation value set, select the Matérn kernel function and initialize the hyperparameters of the Matérn kernel function to construct the initial local Gaussian process model and the initial kernel matrix of the initial local Gaussian process model.

[0106] Step S232: Using the maximum likelihood estimation principle, the hyperparameters of the Matérn kernel function are iteratively optimized through the gradient ascent algorithm to obtain the optimized hyperparameter set and the corresponding optimized kernel matrix. The optimized hyperparameter set and the corresponding optimized kernel matrix together constitute the optimized local Gaussian process model.

[0107] Specifically, the core of step S23 is to construct and optimize a local Gaussian process model containing the Matérn kernel function based on the normalized elevation value set generated in step S22. This solves the problem that traditional fixed-parameter kernel functions cannot adapt to the spatial heterogeneity of terrain (the difference in spatial correlation between plains and mountains), and at the same time provides accurate model support that conforms to local terrain characteristics for subsequent elevation prediction.

[0108] The Matérn kernel function is a covariance function that characterizes the correlation between two points in space. The smoothness parameter ν of the Matérn kernel function determines the continuity level of the terrain surface. We set the smoothness parameter ν = 2.5 for the Matérn kernel function. The reason for choosing the Matérn kernel function is that traditional radial basis functions (RBF) can only characterize infinitely smooth spatial correlations, while real terrain is relatively smooth in plains areas and has dramatic undulations in mountainous areas, resulting in significant differences in smoothness. The Matérn kernel function allows for flexible adjustment of smoothness through the value of ν. When ν = 2.5, it ensures the second-order continuity of the terrain surface (avoiding sharp edges or abrupt changes in interpolation results) while also adapting to the subtle undulations of complex terrains such as mountains, thus solving the fitting distortion problem caused by the fixed smoothness of traditional kernel functions. Subsequently, the hyperparameter set θ of the Matérn kernel function is initialized. The hyperparameter set θ includes the signal variance σ. f 2 Length scale L, noise variance σ n 2 , where σ f 2 The overall fluctuation range of terrain elevation is quantified, L controls the decay rate of terrain spatial correlation, and σ n 2 The allowable noise level during model fitting, along with the signal variance σ, collectively determine the kernel function's ability to characterize local terrain features. f 2 The initial value is set to the variance of the normalized elevation value set. Logically, the variance of the normalized elevation value set directly reflects the fluctuation range of elevation data within the local neighborhood, and is related to σ. f 2 The physical meaning of "controlling the amplitude of terrain undulation" is highly consistent. Matching the initial value with data features can reduce the number of subsequent optimization iterations and improve optimization efficiency. The initial value of the length scale L is set to 0.1 × the spatial span of the local neighborhood point set Nq. The spatial span refers to first calculating the coordinate difference of the local neighborhood point set Nq in the X direction (i.e., the difference between the maximum and minimum X coordinates of all points in Nq) and the coordinate difference in the Y direction (i.e., the difference between the maximum and minimum Y coordinates of all points in Nq), and then taking the larger of these two differences. The basis for setting the coefficient to 0.1 is that L needs to be adapted to the spatial scale of the local terrain. This coefficient can ensure that the initial L can initially characterize the local spatial correlation, avoiding L being too large, which would lead to an overly coarse characterization of the correlation (failing to capture the subtle undulations of the mountains) or too small, which would lead to an overly fine characterization (overfitting noise). Noise variance σ n 2 The initial value is set to a small preset value, for example, 0.01. The reason is that the normalization transformation in step S22 has reduced the noise interference in the original data. The small initial value can avoid the excessive influence of noise on the model fitting and ensure that the model learns the true spatial correlation of the terrain first.

[0109] The initial kernel matrix Kq of the initial local Gaussian process model is constructed. Kq is an nq×nq square matrix composed of the covariances between all points in the local neighborhood point set. It is used to quantify the correlation strength between samples and is the core matrix for Gaussian process regression calculation. The method for constructing the initial kernel matrix Kq includes: for any two points xi and xj in the local neighborhood point set Nq, calculating the Euclidean distance d between the two points. ij Substituting these values ​​into the Matérn kernel function formula (when ν=2.5), we obtain the covariance k(xi,xj) between the two points; the element in the i-th row and j-th column of the initial kernel matrix Kq is Kq[i,j]=k(xi,xj)+δ ij ×σ n 2 Where i and j are index variables of neighborhood points in the local neighborhood point set Nq, i ranges from 1 to nq, j ranges from 1 to nq, and δ ij For the Kronecker function (δ when i=j) ij =1, when i≠j, δ ij =0), add σ n 2 The reason is to avoid the kernel matrix becoming singular (determinant is 0) due to excessive correlation between samples, to ensure the feasibility of subsequent matrix inversion operations, and to solve the problem of computational failure caused by the noise-free terms of the traditional kernel matrix.

[0110] The implementation process of step S232 is as follows: First, determine the objective function for hyperparameter optimization, which is the logarithmic marginal likelihood function, and the formula is:

[0111] ;

[0112] in:

[0113] : The symbol for the probability distribution function, representing "the probability of the target variable under given conditions";

[0114] : The set of normalized elevation values ​​(in column vector form) of a local neighborhood point set;

[0115] Given the two premises of "the two-dimensional horizontal coordinates Xq of the local neighborhood point set" and "the set of hyperparameters θ of the Matérn kernel function", the probability of "observing the set of normalized elevation values ​​Z′q corresponding to the current local neighborhood point set"; where, the symbol... The conditional probability separator indicates the condition under Xq and θ.

[0116] : A set of two-dimensional horizontal coordinates (in matrix form) of a local neighborhood point set;

[0117] : The transpose of (row vector form);

[0118] : The initial kernel matrix (covariance matrix) constructed based on the Matérn kernel function.

[0119] : The inverse matrix;

[0120] : is the logarithm of the determinant of the initial kernel matrix Kq.

[0121] The logic of this formula is: the first term The fitting error term measures the error between the model's predicted values ​​and the observed data (Z'q). A smaller value indicates a smaller fitting error. In Gaussian process regression, the "fitting error between the model's predicted values ​​and the observed values" is quantified by the "quadratic form of the residuals," where... It is the observed normalized elevation value (true value). It is the inverse of the initial kernel matrix, which can be understood as a "spatial correlation weight matrix," assigning weights to the residuals at different locations that are adapted to the local terrain correlation; the negative sign ( )and These are the inherent coefficients of the Gaussian probability density function, ensuring that "the smaller the error, the larger the function value".

[0122] Second item The regularization term of the initial kernel matrix prevents the model from overfitting to local noise during hyperparameter optimization, ensuring the model's ability to depict the true spatial relevance of the terrain and avoiding a decrease in generalization ability due to overfitting; the determinant of the kernel matrix Reflecting the "spatial correlation strength of local data", if If the value is too large, it indicates that the spatial correlation of the data is too strong, and the model is prone to "underfitting"; if If the value is too small, it indicates that the spatial correlation of the data is too weak, and the model is prone to "overfitting"; adding Multiplying by a negative sign balances the need for "small fitting error" and "strong model generalization ability," preventing the model from overfitting in pursuit of minimizing fitting error. (Third term) This is a constant term that does not affect the direction of hyperparameter optimization; it only ensures the mathematical integrity of the logarithmic marginal likelihood function, making it strictly conform to the marginal likelihood formula of the Gaussian distribution. A larger overall function value indicates a stronger interpretability of the current hyperparameter combination for local terrain data, thus addressing the deficiency of traditional empirical hyperparameter setting without quantitative evaluation standards.

[0123] The core function of the objective function for hyperparameter optimization is to provide a quantitative evaluation standard for the hyperparameter optimization of the local Gaussian process model. Traditional methods often set hyperparameters empirically, lacking objective evaluation criteria. This formula, however, transforms the quality of hyperparameters into the magnitude of the logarithmic marginal likelihood function value. The larger the function value, the better the fit of the set of hyperparameters to the local terrain data, thus providing a clear optimization objective for the gradient ascent algorithm (i.e., updating hyperparameters in the direction of increasing function value). This formula relies on the quantile-quantile normalization transformation in step S22, and is constructed based on the normalized elevation value set Z'q, enabling the Gaussian process model to adapt to the original non-Gaussian model. Topographic data avoids fitting bias caused by inconsistent data distribution; by optimizing the hyperparameter set θ, the model can automatically adapt to the spatial correlation characteristics of different terrains (e.g., L is larger in plains and smaller in mountainous areas), solving the problem that traditional fixed-parameter kernel functions cannot adapt to the spatial heterogeneity of terrain; the regularization term in the formula can effectively avoid the model overfitting local noise, achieving a balance between fitting accuracy and predictive ability for unknown areas; logarithmic transformation avoids the numerical overflow problem of the original marginal likelihood function (exponential form), and gradient calculation can be efficiently derived through the chain rule, which can adapt to the real-time processing needs of massive point cloud data.

[0124] Based on the objective function, the gradient ascent algorithm is used to iteratively optimize the hyperparameters: the first step is to calculate the log marginal likelihood function for each hyperparameter (σ). f 2 L, σ n 2 The gradient of the function is calculated, with the gradient direction being the direction in which the function value increases. The gradient expression is derived using the chain rule to ensure that each update moves closer to the optimal hyperparameters. The second step updates the hyperparameters along the gradient direction. The update step size is set according to computational resources and convergence speed, using an adaptive step size strategy, such as adjusting according to the gradient magnitude. When the gradient magnitude is large, the step size is reduced to avoid oscillation; when the gradient magnitude is small, the step size is increased to accelerate convergence. The third step determines whether convergence has occurred. The convergence condition is that the difference between the logarithmic marginal likelihood function values ​​of two adjacent iterations is less than a preset difference threshold, which can be set to 1×10 for example. -6 Or, the number of iterations reaches a preset maximum, for example, 100 times. After convergence, the optimized hyperparameter set θ'=[σ'] is obtained. f 2 , L', σ' n 2 ], where σ' f 2 Let L' be the optimized signal variance, L' be the optimized length scale, and σ' be the variable. n 2 The optimized noise variance is calculated, and the kernel matrix is ​​recalculated based on θ' to obtain the optimized kernel matrix K'q.

[0125] Step S23 addresses the core issue of the inability of fixed-parameter kernel functions in existing technologies to adapt to the spatial heterogeneity of terrain. It employs a technique of "adaptive kernel function selection (Matérn kernel function ν=2.5) + quantization of the objective function (logarithmic marginal likelihood function) + gradient ascent optimization" to achieve accurate model adaptation to local terrain features. Quantization optimization of the logarithmic marginal likelihood function ensures optimal hyperparameters, significantly improving accuracy compared to empirical settings. Iterative updates of the gradient ascent algorithm avoid getting trapped in local optima, ensuring global model adaptability. The optimized length scale L' indirectly reflects the complexity of local terrain; a large L' value indicates gentle terrain (e.g., plains) with slow spatial correlation decay, while a small L' value indicates complex terrain (e.g., mountains) with rapid spatial correlation decay. This information can serve as an auxiliary basis for terrain classification, expanding the terrain analysis function of the method. Optimizing the determinant of the kernel matrix K'q quantifies the correlation strength of local data; a large determinant value indicates strong data correlation and low model prediction uncertainty, providing additional evidence for judging the reliability of subsequent prediction results. If step S23 is missing, the subsequent step S24 will only be able to use a kernel function with fixed parameters, which cannot adapt to the terrain features of different regions. This will result in overly smooth predictions in plain areas (losing subtle undulations) and large prediction deviations in mountainous areas (failing to capture dramatic undulations). At the same time, the lack of hyperparameter optimization will prevent the kernel matrix from reflecting the true correlation of local data, which will lead to a significant decrease in the accuracy of elevation prediction and make it impossible to achieve the "high fidelity" goal of the entire scheme.

[0126] Step S24: Based on the optimized local Gaussian process model, predict the mean normalized elevation of the query point;

[0127] The core of step S24 is based on the optimized local Gaussian process model obtained in step S23. It takes the two-dimensional horizontal coordinates of the query point as input and outputs the normalized elevation prediction mean and variance. This addresses the problem that traditional interpolation methods only provide point estimates and lack uncertainty quantification, providing a basis for judging the reliability of the prediction results. The implementation process is as follows: First, obtain the two-dimensional horizontal coordinates (xq, yq) of the query point q, the optimized hyperparameter set θ', and the optimized kernel matrix K'q. Calculate the covariance vector kq*, where kq* is a 1×nq vector composed of the covariances of the query point q and all points in the local neighborhood point set Nq. A large value for the i-th element kq*(q, xi) in the covariance vector kq* indicates a strong correlation between the query point q and the neighborhood point xi, and that the neighborhood point xi has a significant impact on the elevation prediction of the query point q, thus solving the problem that traditional fixed-parameter covariance calculations cannot reflect local correlation differences.

[0128] Then, the prior variance kq** is calculated, where kq** represents the covariance of the query point q itself, and the formula for kq** is kq** = σ'. f 2 +σ' n 2The logic is as follows: the prior variance must include the fluctuations of the signal itself (σ'). f 2 ) and the noise allowed by the model (σ') n 2 This involves initially quantizing the uncertainty range of the query point elevation to avoid estimation biases caused by considering only the signal or noise. For example, if σ' f 2 =0.8、σ' n 2 =0.02, then kq**=0.82, indicating that the initial fluctuation range of the elevation of the query point is dominated by the signal and the noise has little impact.

[0129] Based on Gaussian process regression theory, the formula for calculating the posterior distribution: the normalized elevation prediction mean μq of query point q is μq=(kq*). T ×(K'q) -1 ×Z'q, where Z'q is the set of normalized elevation values. The logic of this formula is: (K'q) -1 The weight coefficients of neighboring points are calculated using ×Z'q, and then multiplied by the covariance vector kq* (association strength) to obtain the predicted elevation value of the query point. This weight allocation reflects the differences in influence among neighboring points; neighboring points with strong associations have higher weights, resulting in higher accuracy compared to traditional equal-weight interpolation. The normalized elevation prediction variance σq of the query point q... 2 The formula is σq 2 =kq**-(kq*) T ×(K'q) -1 ×kq*, the logic is: use the correlation strength between the query point and its neighboring points to correct the prior variance. If the value of kq* is large (strong correlation), then the correction term is large, σq 2 A smaller normality indicates lower uncertainty, and vice versa, overcoming the limitation of traditional interpolation in quantifying uncertainty. Finally, the 95% confidence interval is calculated: Lower confidence bound normality = μq - 1.96 × The upper bound of confidence is normality = μq + 1.96 × The value of 1.96 represents the 95th percentile of the standard normal distribution, derived from the "3σ principle" of the normal distribution. This ensures that the true normalized elevation of the query point falls within this interval with a 95% probability, providing a clear boundary for subsequent reliability assessment.

[0130] Step S25: The normalized elevation prediction mean of the query point is transformed by inverse quantile transformation through the stored transformation mapping table to obtain the true elevation prediction value of the query point in the tile to which it belongs.

[0131] The core of step S25 is to use the transformation mapping table Mq constructed in step S22 to inversely transform the normalized prediction result obtained in step S24 into the true value of the original elevation space, thereby solving the problem that the normalized spatial result has no practical topographic meaning (cannot directly reflect the true elevation) and meeting the practical application requirements of topographic interpolation.

[0132] The implementation process of step S25 is as follows: First, query the normalized elevation prediction mean μq of point q, and calculate the cumulative probability pq=Φ(μq), where Φ is the cumulative distribution function of the standard normal distribution. The purpose is to map the mean of the normalized space to a probability value in the interval [0,1]. This probability value is consistent with the interval of the empirical cumulative probability calculated in step S22, providing an index basis for finding Mq and avoiding mapping failure due to spatial mismatch. The cumulative probability pq represents the probability position of the normalized elevation prediction mean μq in the standard normal distribution.

[0133] Then, the index j' is searched in the transformation mapping table Mq: the empirical cumulative probability p in Mq is traversed. (i') Find the one that satisfies p (j') ≤pq≤p (j'+1) j', where p (j') p represents the empirical cumulative probability of the j'-th mapping pair (j'+1) This represents the cumulative empirical probability of the (j'+1)th mapping pair. The logic of this search is: pq falls on p (j') With p (j'+1) The interval between these values ​​indicates that the actual elevation value of the query point is between the corresponding original elevation value z. (j') With z (j'+1) This ensures the uniqueness and accuracy of the mapping, resolving the multi-valued mapping problem that may exist in traditional global inverse transforms. For example, if p in Mq... (10) =0.5625, p (11) =0.5938, and pq=0.5793, then j'=10, indicating that the true elevation is between z and q. (10) With z (11) Between. The predicted true elevation value μq* is calculated using linear interpolation, with the following formula:

[0134]

[0135] in, and These are the original elevation values ​​of the j' and j'+1th mapping pairs in the transformation mapping table Mq, respectively. This formula calculates the cumulative probability pq in... and The ratio between ,exist and Interpolation is performed between them to avoid accuracy loss caused by the discreteness of Mq, thus solving the problem of traditional discrete mapping (directly taking the nearest value). The error problem. Similarly, perform the above inverse transformation on the lower confidence bound normal and upper confidence bound normal of the confidence interval in step S24: first calculate the corresponding cumulative probability pq. 下界 and pq 上界 Then, by finding Mq and linear interpolation, the lower and upper bounds of the true elevation confidence are obtained, making the quantification of uncertainty in the original elevation space practically meaningful.

[0136] S24 addresses the deficiency of traditional interpolation methods in lacking uncertainty quantification, providing a basis for reliable prediction results. S25 addresses the deficiency of normalized spatial results lacking practical topographic significance, transforming the prediction results into directly applicable real elevation values. The synergy between the preceding and following steps is manifested in the following ways: S24 relies on the optimized local Gaussian process model of step S23; if the model parameters are not optimized, the accuracy of the predicted mean and variance will be low, failing to reflect the true terrain features. S25 relies on the transformation mapping table Mq of step S22 and the normalization result of step S24; without Mq, the inverse transformation cannot be achieved, and without the normalization result, there is no transformation object. The real elevation prediction value output by S25 provides the basic data for tile fusion in step S30; without it, there is no input for fusion, and global terrain surface generation cannot be achieved.

[0137] The uncertainty quantification in step S24 can guide subsequent data supplementation, σq 2 Large areas (such as the edges of mountainous regions with complex terrain) indicate low prediction reliability, requiring increased point cloud data acquisition, which is more efficient than blind acquisition using traditional interpolation. The posterior variance correction mechanism ensures accurate uncertainty estimation and avoids the coarseness of prior variance (considering only signal and noise without incorporating neighborhood correlation). The linear interpolation in step S25 guarantees the accuracy of the true elevation, significantly reducing errors compared to traditional discrete mapping, especially in areas where the original elevation values ​​are sparsely distributed. The one-to-one association of Mq ensures the uniqueness of the inverse transformation, avoiding accuracy loss caused by a globally unified mapping table (ignoring local data distribution differences). The true elevation confidence interval provides a basis for risk assessment in terrain engineering applications such as earthwork calculation and road design. A smaller interval indicates lower engineering risk, while a larger interval requires careful adjustment of the design scheme.

[0138] Step S30: Process all query points within the tile in parallel to obtain the true elevation prediction value of all query points within the tile. Obtain the tile elevation interpolation result set based on the true elevation prediction value of all query points within the tile. Perform a weighted average on each tile elevation interpolation result set to obtain the global seamless terrain surface.

[0139] The core of step S30 is to solve the problems of low processing efficiency of massive query points and discontinuous local interpolation boundaries by using the process of "parallel processing of query points within tiles - dynamic adjustment of overlapping areas - weighted fusion of tile results", and to integrate the local elevation prediction results of step S20 into a global seamless terrain surface.

[0140] Further, step S30 includes:

[0141] Step S31: Perform parallel processing on all query points within the tile to obtain the true elevation prediction value of all query points within the tile, and obtain the tile elevation interpolation result set based on the true elevation prediction value of all query points within the tile.

[0142] The core of step S31 is based on the tile division in step S13. Each tile is assigned to an independent parallel computing unit, and the local modeling process of steps S21-S25 is executed synchronously to generate a set of tile elevation interpolation results. The width of the overlapping area is dynamically adjusted by collecting and optimizing the length scale, which solves the problems of low efficiency of traditional serial processing and poor fusion quality caused by insufficient overlapping area.

[0143] The specific implementation process of step S31 includes: obtaining the tile set T generated in step S13 and the number of CPU cores P of available computing resources. Let the number of tiles be Nt. Take the CPU cores as computing units and allocate computing tasks according to the ratio of the number of tiles to P: if Nt≤P, each tile is allocated to an independent computing unit; if Nt>P, the tiles are evenly allocated to each computing unit in a round-robin manner to ensure that the load of each unit is balanced and to avoid overload caused by too many tiles in a certain unit, while a certain unit is idle.

[0144] Each computing unit performs the following operations on the assigned tile: traverse all interpolation grid points within the tile, sequentially execute the adaptive hybrid neighborhood search in step S21, the quantile-quantile normalization transformation in step S22, the local Gaussian process optimization in step S23, the elevation prediction in step S24, and the inverse transformation in step S25 to obtain the true elevation prediction values ​​of all query points within the tile. The true elevation prediction values ​​of all query points within the tile are then summarized into the elevation interpolation result set for that tile. At the same time, the optimized length scale of that tile output by each computing unit is recorded.

[0145] Calculate the arithmetic mean of the optimized length scales of the tiles output by all computing units to obtain the global average length scale L*. Update the overlapping region width M based on the global average length scale L*. Calculate the optimized overlapping region width M* = max(3 × L*, M) based on L* and M. If the difference between M* and the initially set overlapping region width M is greater than or equal to a preset deviation threshold (for example, the deviation threshold can be 10% of the initially set overlapping region width M), it is determined that the initial overlapping region is insufficient, and the tile division in step S13 and the parallel processing in step S31 are re-executed to ensure sufficient transition space during subsequent fusion; if the difference is less than the deviation threshold, the current tile division result and overlapping region width M are retained.

[0146] Traditional methods perform local modeling serially for all query points, resulting in processing time that increases linearly with the number of query points, failing to meet the demands of massive data volumes. This invention, based on tile partitioning in step S13, decomposes the task into Nt independent subtasks, significantly reducing parallel processing time and resolving the efficiency bottleneck. Traditional methods use fixed overlapping regions, which cannot adapt to variations in terrain complexity. This invention dynamically adjusts M using L*, with 3×L* set based on the fact that the terrain correlation decay range is approximately three times the length scale. The overlapping region must cover this range to ensure a smooth transition between adjacent tile results, resolving the instability in fusion quality caused by fixed overlap. Dynamically updating the overlapping region ensures sufficient transition space for regions with different terrain complexities, avoiding abrupt boundary changes during subsequent fusion in step S32 and providing more data support for smooth fusion.

[0147] Step S32: Identify all tiles belonging to each query point, calculate the fusion weight of all tiles at the query point, and perform a weighted average of the elevation interpolation results of each tile according to the fusion weight to obtain a global seamless terrain surface.

[0148] The core of step S32 is to eliminate the discontinuity of adjacent tile boundaries and generate a globally seamless terrain surface through the process of "identifying covering tiles - calculating normalized distance - smoothing weight allocation - weighted average fusion," thus solving the problem of the "tile effect" (abrupt boundary) in traditional local interpolation. The implementation method is as follows: traverse all query points and identify the set of covering tiles Sq covering each query point, where Sq represents the set of covering tiles covering query point q. The method for identifying the set of covering tiles Sq is as follows: for each tile in the tile set T, determine whether the coordinates of query point q fall within the spatial range of the tile. If they do, the tile is defined as a covering tile and included in Sq. Calculate the normalized distance (u) from query point q to the center of each covering tile in Sq. i* ,v i* ), where i* is the index of the covering tile in the set of covering tiles Sq, u i*Let v be the normalized distance from query point q in the X direction to the center of the i*th covering tile. i* Let q be the normalized distance in the Y direction from the query point to the center of the i*th overlay tile; then calculate the fusion weight w of each overlay tile at the query point q. i* (q):

[0149] Normalized distance u i* and v i* Substitute into the smooth decay function , if u i* If S(u) > 1, then S(u) i* If v = 0; i* If S(v) > 1, then S(v) i* )=0; the weight formula is: w i* (q)=S(u i* )×S(v i* ); where, smooth decay function The cubic polynomial form is:

[0150]

[0151] Where t is the normalized distance from the query point to the center of the covering tile, this function ensures that when t=0... That is, the query point is located at the center of the tile, and has the highest weight; hour The query point is at the tile boundary and has a weight of 0. Meanwhile, derivative exist Continuous in time: ; and This achieves a smooth transition of weights. The elevation interpolation results of all covered tiles in Sq are aggregated, and a weighted average is performed by fusion weights to obtain the final elevation interpolation result; after traversing all query points to complete the fusion, a global seamless terrain surface is obtained.

[0152] Traditional methods employ linear decay functions, such as S(t) = 1 - t, where the derivative is -1 when t ∈ [0,1]. This discontinuity leads to abrupt weight changes at tile boundaries, resulting in elevation abrupt changes. This invention uses a cubic polynomial function to ensure that both the function value and the derivative are continuous when t ∈ [0,1], and the weight smoothly decreases from 1 to 0, solving the boundary abruptness problem. Traditional methods use simple averaging, failing to consider the intensity of tile influence on the query point. For example, tiles far from the query point contribute the same amount as those nearby, leading to distortion in the fusion result. This invention uses weighted fusion to ensure that tiles near the query point contribute more and tiles far away contribute less, solving the distortion problem caused by uneven influence intensity. The smooth decay function ensures continuous weight transition at adjacent tile boundaries, solving the core problem of "discontinuous local interpolation boundaries." Weighted averaging integrates the interpolation results of multiple tiles, rather than a single tile, making the final elevation interpolation result closer to the actual terrain elevation.

[0153] Example 2

[0154] This embodiment, based on Embodiment 1, provides experimental verification of the high-fidelity terrain surface interpolation method for large-scale point cloud data in Embodiment 1, as detailed below:

[0155] To verify the computational efficiency and high-fidelity effect of the high-fidelity terrain surface interpolation method for large-scale point cloud data described in Embodiment 1 of this invention in practical applications, this experiment uses a point cloud dataset containing three million points. Before the experiment, a multi-level cleaning process was carried out to remove redundant and duplicate points, measurement noise points, and elevation conflict points at the same horizontal position from the original point cloud data, ensuring the accuracy and consistency of the input data and laying a reliable data foundation for subsequent interpolation experiments. During the experiment, the hardware environment (number of CPU cores, memory capacity) and software configuration of all interpolation methods were kept consistent to eliminate the interference of external condition differences on the experimental results.

[0156] First, based on the cleaned three million point cloud data, the IDW interpolation method, the Kriging interpolation method, and the full-process interpolation method of "data terrain reading-processing-regional-local interpolation-global smoothing processing" described in Embodiment 1 of this invention were used to perform interpolation calculations on the same data to be interpolated. The complete time taken by the three methods from data loading to generating the final interpolation result was recorded simultaneously. At the same time, the interpolation calculation time of a single point to be interpolated was additionally recorded for the method of this invention. As shown in Table 1, the calculation time for the IDW interpolation method is approximately 60 minutes, and for the Kriging interpolation method it is approximately 80 minutes. In contrast, the complete process of the interpolation method of this invention takes only about 4 minutes, representing a speed improvement of more than 15 times compared to IDW interpolation and more than 20 times compared to Kriging interpolation. Furthermore, the interpolation calculation time for a single point to be interpolated (query point) is approximately 0.02 seconds, enabling rapid response to the elevation prediction needs of a single point. This efficiency advantage provides crucial support for real-time interactive operations of subsequent 3D intelligent route selection products on large-scale terrain data, such as real-time adjustment of route selection paths and real-time preview of terrain cutting effects. It solves the problem that traditional interpolation methods, due to their excessive time consumption, cannot meet real-time interactive requirements.

[0157] Table 1. Comparison of Interpolation Method Time Consumption

[0158] Interference method time consuming IDW interpolation method 60 minutes Kriging interpolation method 80 minutes The interpolation method of the present invention 4 minutes

[0159] To verify the high-fidelity characteristics of the method of this invention, experiments were conducted from three dimensions: consistency between the interpolation results and the true values, the reproducibility of terrain features, and the error distribution. First, a scatter plot comparing the interpolation results and the true elevations was generated, as shown below. Figure 6 As shown, this graph uses the true elevation value (z_true) as the horizontal axis and the interpolated result (i.e., the calculated elevation z_pred) as the vertical axis. The red dashed line is the y=x reference line, representing that the interpolation result is completely consistent with the true value; from Figure 6 As can be visually observed, the calculation results of the interpolation method of this invention are closely distributed near the red line, with only a very few points showing slight deviations. This indicates that the deviation between the interpolation results and the true elevation values ​​is extremely small, preliminarily verifying the high-fidelity characteristics of the method of this invention. Secondly, the experiment simultaneously generated heatmaps of the original data, interpolation results, and errors, such as... Figure 7As shown, the spatial coordinate ranges of the three types of heatmaps are kept consistent to ensure the effectiveness of the comparison. The original data heatmap clearly presents the elevation distribution characteristics of the experimental point cloud data, including the spatial location and range of high elevation areas (darker color areas) and low elevation areas (lighter color areas). The color distribution pattern of the interpolated heatmap is highly consistent with that of the original data heatmap. Especially in areas with complex terrain undulations, the interpolated heatmap accurately reproduces the elevation change trend of the original terrain, without obvious loss or distortion of terrain features. In the error heatmap, the error values ​​of most areas are within a very small range (the color is close to the unbiased reference color), with only slight errors in local edge areas. This further proves that the method of this invention can accurately reproduce the surface features of the original terrain, providing reliable interpolation data for subsequent terrain analysis and engineering design.

[0160] The experiment also performed statistical analysis on the prediction error of the interpolation method of the present invention, generating a histogram of the prediction error distribution, such as... Figure 8 As shown in the figure, the zero error line is the baseline of the error distribution, the vertical axis is the frequency (error frequency), and the horizontal axis is the error. Statistical calculations show that the error distribution conforms to the characteristics of a Gaussian distribution, with a mean (μ) of 0.20 and a standard deviation (σ) of 14.47. The histogram shows that the vast majority of error values ​​are concentrated in the range of -50 to 50, and are symmetrically distributed around the zero error line. Only a very small number of error values ​​exceed the range of -100 to 100. This result indicates that the interpolation error of the method of the present invention is relatively small and uniformly distributed, with no systematic deviation. It can meet the stringent requirements for the accuracy of elevation prediction in scenarios such as road engineering design and terrain modeling, and avoid safety hazards or economic problems in engineering design schemes due to excessive errors.

[0161] Furthermore, the experiments also verified the adaptability of the method of this invention in different terrain scenarios, including mountainous terrain and plain / hilly terrain. Figures 9-10 As shown. In mountainous terrain scenarios, after interpolation using the method of this invention, the system can clearly present the complex terrain undulations of mountainous areas, including steep slopes, valleys, and other micro-topographic details. Furthermore, the subsequent longitudinal and cross-sectional data of the route generated based on the interpolation results can accurately reflect changes in terrain elevation, providing precise data support for route slope optimization and roadbed excavation calculation. In plain and hilly terrain scenarios, the interpolation results can smoothly reproduce the gentle terrain of plains and the gentle slope characteristics of hilly areas, without the "tile effect" (abrupt changes in boundary elevation) that is prone to occur in traditional local interpolation methods. Moreover, the interpolation results have a high degree of matching with the actual terrain, which can accurately guide the layout of the surrounding landscape and the selection of transportation facilities, verifying the adaptability and stability of the method of this invention under different terrain types.

[0162] This experiment, through comparison of computational efficiency, verification of high-fidelity characteristics, and testing of adaptability to different terrain scenarios, fully demonstrates that the interpolation method described in Embodiment 1 of this invention can achieve efficient and high-fidelity terrain surface interpolation on large-scale point cloud data. Compared with existing IDW and Kriging interpolation methods, it has significant efficiency and accuracy advantages, and can adapt to various terrain scenarios in road engineering design. It provides high-quality terrain data support for downstream applications such as 3D intelligent route selection, roadbed design, and landscape planning, and solves the core problems of poor scalability, low interpolation fidelity, and weak terrain adaptability of existing technologies in large-scale point cloud processing.

[0163] Example 3

[0164] This embodiment, based on Embodiment 1, provides a high-fidelity terrain surface interpolation system for large-scale point cloud data, such as... Figure 11 As shown, it includes:

[0165] Point cloud data cleaning module: Input the raw point cloud dataset, perform multi-level cleaning on the raw point cloud dataset, and generate a cleaned point cloud dataset;

[0166] Tile partitioning module: Used to partition the tile grid based on the cleaned point cloud dataset, generate a set of tiles covering the entire target area, and generate a set of grid points to be interpolated within each tile in the tile set;

[0167] Elevation prediction module: The interpolation grid points in the set of interpolation grid points are used as query points. A local Gaussian process model is constructed and optimized. Based on the optimized local Gaussian process model, the true elevation prediction value of the query point in the tile to which it belongs is obtained.

[0168] Terrain surface generation module: It is used to process all query points within a tile in parallel, obtain the true elevation prediction value of all query points within the tile, obtain the tile elevation interpolation result set based on the true elevation prediction value of all query points within the tile, and perform a weighted average on the tile elevation interpolation result set to obtain the global seamless terrain surface.

[0169] Furthermore, in the point cloud data cleaning module, the method for multi-level cleaning of the original point cloud dataset includes:

[0170] Step S111: Construct a hash mapping table based on horizontal position, and divide the 3D coordinate points in the original point cloud dataset into multiple horizontal position groups according to the hash mapping table; each horizontal position group has a single elevation value or multiple elevation values.

[0171] Step S112: Perform statistical analysis on horizontal position groups containing multiple elevation values ​​and remove outliers to obtain a set of valid elevation values;

[0172] Step S113: Aggregate the set of valid elevation values ​​within each horizontal position group to generate a cleaned point cloud dataset with a unique elevation value for each horizontal position.

[0173] Furthermore, in the elevation prediction module, the method for constructing and optimizing a local Gaussian process model, and obtaining the predicted true elevation of the query point in its respective tile based on the optimized local Gaussian process model, includes:

[0174] Step S21: Take the interpolation grid points in the set of interpolation grid points as query points, and use the K-dimensional tree space index structure to perform adaptive hybrid neighborhood search on any query point to obtain the local neighborhood point set of the query point.

[0175] Step S22: Extract the representative elevation values ​​of all neighboring points in the local neighborhood point set of the query point to form the original elevation value set. Transform the original elevation value set into a normalized elevation value set through quantile-quantile normalization transformation, and at the same time construct and store the transformation mapping table.

[0176] Step S23: Based on the normalized elevation value set, construct a local Gaussian process model based on the Matérn kernel function, and optimize the hyperparameters of the Matérn kernel function through maximum likelihood estimation to obtain the optimized local Gaussian process model.

[0177] Step S24: Based on the optimized local Gaussian process model, predict the mean normalized elevation of the query point;

[0178] Step S25: The normalized elevation prediction mean of the query point is transformed by inverse quantile transformation through the stored transformation mapping table to obtain the true elevation prediction value of the query point in the tile to which it belongs.

[0179] Further, step S23 includes:

[0180] Step S231: Based on the normalized elevation value set, select the Matérn kernel function and initialize the hyperparameters of the Matérn kernel function to construct the initial local Gaussian process model and the initial kernel matrix of the initial local Gaussian process model.

[0181] Step S232: Using the maximum likelihood estimation principle, the hyperparameters of the Matérn kernel function are iteratively optimized through the gradient ascent algorithm to obtain the optimized hyperparameter set and the corresponding optimized kernel matrix. The optimized hyperparameter set and the corresponding optimized kernel matrix together constitute the optimized local Gaussian process model.

[0182] The methods and systems of this application may be implemented in many ways. For example, they may be implemented by software, hardware, firmware, or any combination of software, hardware, and firmware. The above-described order of steps for the method is for illustrative purposes only, and the steps of the method of this application are not limited to the order specifically described above, unless otherwise specifically stated.

[0183] In addition, the parts of the technical solutions provided in the embodiments of this application that are consistent with the implementation principles of the corresponding technical solutions in the prior art have not been described in detail, so as to avoid excessive elaboration.

[0184] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the invention. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A high-fidelity terrain surface interpolation method for large-scale point cloud data, characterized in that, The method includes: Input the original point cloud dataset, perform multi-level cleaning on the original point cloud dataset to generate a cleaned point cloud dataset; based on the cleaned point cloud dataset, perform tile mesh division to generate a set of tiles covering the entire target area, and generate a set of grid points to be interpolated in each tile of the tile set; Using the interpolation grid points in the set of interpolation grid points as query points, a local Gaussian process model is constructed and optimized. Based on the optimized local Gaussian process model, the predicted true elevation of the query point in its respective tile is obtained. All query points within a tile are processed in parallel to obtain the true elevation prediction values ​​of all query points within the tile. Based on the true elevation prediction values ​​of all query points within the tile, a set of tile elevation interpolation results is obtained. The weighted average of each set of tile elevation interpolation results is then performed to obtain a global seamless terrain surface.

2. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 1, characterized in that, Each point in the original point cloud dataset contains three-dimensional coordinate information, which refers to the horizontal position coordinates X and Y and the elevation value Z. The multi-level cleaning method includes: Construct a hash mapping table based on horizontal position, and divide the 3D coordinate points in the original point cloud dataset into multiple horizontal position groups according to the hash mapping table; each horizontal position group has a single elevation value or multiple elevation values. Statistical analysis was performed on horizontal location groups containing multiple elevation values, and outliers were removed to obtain a set of valid elevation values. Aggregate the valid elevation values ​​within each horizontal location group to generate a cleaned point cloud dataset with a unique elevation value for each horizontal location.

3. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 2, characterized in that, The method for dividing the horizontal position groups includes: Traverse the original point cloud dataset, extract the horizontal position coordinates (X,Y) of each point in the original point cloud dataset to generate a unique key value, and store the corresponding elevation value Z in the list under the key value to form a horizontal position group.

4. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 3, characterized in that, The method for constructing and optimizing a local Gaussian process model includes: Obtain the local neighborhood point set of the query point, construct the original elevation value set based on the local neighborhood point set, transform the original elevation value set into a normal elevation value set through quantile-quantile normalization transformation, and construct and store the transformation mapping table at the same time; Based on the normalized set of elevation values, a local Gaussian process model based on the Matérn kernel function is constructed. The hyperparameters of the Matérn kernel function are optimized by maximum likelihood estimation to obtain the optimized local Gaussian process model.

5. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 4, characterized in that, The method for obtaining the local neighborhood point set of the query point includes: A recursive segmentation strategy based on the median is used to process the cleaned point cloud dataset, and a K-dimensional tree spatial index structure is constructed. An adaptive hybrid neighborhood search is performed on any query point using a K-dimensional tree spatial index structure to obtain the local neighborhood point set of the query point.

6. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 5, characterized in that, The method for performing adaptive hybrid neighborhood search for any query point using a K-dimensional tree space index structure includes: Determine the initial radius and density threshold, use a K-dimensional tree space index structure to perform an initial radius search with the query point as the center, and count the number of points Nr within the initial radius; The density threshold includes a lower limit of density threshold kmin and an upper limit of density threshold kmat; When \(k_{min}\leq N_r\leq k_{max}\), all \(N_r\) points within the initial radius \(r'\) are used as neighborhood points to form the local neighborhood point set of the query point; When \(N_r < k_{min}\), the \(k_{min}\) points closest to the query point are obtained as the neighborhood points of the query point, and the local neighborhood point set of the query point is constructed; When \(N_r>k_{max}\), the \(k_{max}\) points closest to the query point are selected as the neighborhood points of the query point to form the local neighborhood point set of the query point.

7. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 6, characterized in that, The method for transforming the original elevation value set into a normalized elevation value set includes: Sort the original elevation values in the original elevation value set in ascending order to obtain an ordered elevation value sequence; Calculate the empirical cumulative probability for each original elevation value in the ordered elevation value sequence, and apply the inverse cumulative distribution function of the standard normal distribution to generate a normalized elevation value set.

8. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 7, characterized in that, The method for optimizing the hyperparameters of the Matérn kernel function by maximum likelihood estimation includes: Determine the objective function for hyperparameter optimization. The objective function is the log marginal likelihood function, which includes three terms. The first term is the fitting error term, the second term is the regularization term of the initial kernel matrix, and the third term is a constant term; the fitting error term is constructed based on the normalized elevation value set; Based on the objective function, use the gradient ascent algorithm to iteratively optimize the hyperparameters to obtain an optimized hyperparameter set.

9. The high-fidelity terrain surface interpolation method for large-scale point cloud data according to claim 8, characterized in that, The method for obtaining the predicted true elevation value of the query point in the tile it belongs to includes: Based on the optimized local Gaussian process model, predict the normalized elevation prediction mean of the query point; Perform inverse quantile transformation on the normalized elevation prediction mean of the query point through the stored transformation mapping table to obtain the predicted true elevation value of the query point in the tile it belongs to.

10. A high-fidelity terrain surface interpolation system for large-scale point cloud data, used to implement the high-fidelity terrain surface interpolation method for large-scale point cloud data as described in any one of claims 1-9, characterized in that, The system includes: Point cloud data cleaning module: Input the original point cloud data set, perform multi-level cleaning on the original point cloud data set, and generate a cleaned point cloud data set; Tile division module: Used to perform tile grid division according to the cleaned point cloud data set, generate a tile set covering the entire target area, and generate a set of interpolation grid points to be interpolated within each tile in the tile set; Elevation value prediction module: Use the interpolation grid points to be interpolated in the set of interpolation grid points to be interpolated as query points, construct and optimize the local Gaussian process model, and based on the optimized local Gaussian process model, obtain the predicted true elevation value of the query point in the tile it belongs to; Terrain surface generation module: Used to process all query points within the tile in parallel, obtain the predicted true elevation values of all query points within the tile, and obtain a tile elevation interpolation result set according to the predicted true elevation values of all query points within the tile; Perform weighted averaging on each tile elevation interpolation result set to obtain a global seamless terrain surface.

Images