Three-dimensional plot drawing method based on GIS and digital twinborn fusion

Abstract

The invention relates to the technical field of three-dimensional point cloud construction, and further relates to a three-dimensional plot drawing method based on GIS and digital twinborn fusion. The method comprises the following steps: step 1, generating original digital twin three-dimensional point cloud data of a target plot area based on three-dimensional point cloud through GIS (Geographic Information System) data; 2, performing K-neighbor search or radius neighborhood search on the original digital twin three-dimensional point cloud data, and determining a neighborhood point set of each point; 3, determining a position strategy of new point generation in each iteration; defining a fractal dimension-increasing operator, and outputting an elevation adjustment value; and 4, applying a fractal dimension-increasing operator, and adjusting the elevations of all the points according to the elevation adjustment values to complete three-dimensional plot drawing. The reality sense and the detail level of the three-dimensional plot model are remarkably improved, and the problems that a traditional method is insufficient in details and lacks natural complexity are solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of three-dimensional point cloud construction, and specifically relates to a three-dimensional land parcel drawing method based on the fusion of GIS and digital twins. Background Art

[0002] In the development of modern society, accurate three-dimensional digital representation of the earth's surface morphology is of great significance and has widespread application demand. Whether serving urban planning and management, farmland plot analysis for precision agriculture, soil and water conservation and environmental monitoring, civil engineering design and earthwork calculation, natural disaster simulation and risk assessment, or meeting the demand for realistic scenes in fields such as virtual reality, game development, and film and television production, high-precision, highly realistic three-dimensional land parcel models are indispensable basic data support. With the rapid development of computer graphics, remote sensing technology, and geographic information science, the technical means of acquiring and constructing three-dimensional surface models are becoming increasingly diverse, but at the same time, new challenges are also emerging, especially in how to efficiently generate specific land parcel models that both conform to macroscopic geographical facts and display the complex textures and details of natural landforms at the microscale.

[0003] Models generated based on traditional GIS data (especially medium- and low-resolution DEMs or contours) often lack detail and realism. Standard DEM data has a limited resolution (for example, common SRTM data is 30 or 90 meters, and even national-level basic DEMs are mostly 5 or 10 meters). This fails to capture subtle surface undulations, gullies, ridges, scours, and other topographic details, which are crucial for applications such as hydrological analysis, soil erosion simulation, and precision agricultural management. Surfaces generated based on contour interpolation are prone to unnaturally smooth transitions in flat areas or areas with sparse contours, making it difficult to accurately represent the true curvature of the terrain and even creating interpolation artifacts resembling "terracing." These models may meet accuracy requirements at a macroscopic level, but when observed up close or requiring detailed analysis, their overly smooth or regular appearance deviates significantly from the complexity of the natural landscape. While photogrammetry methods can generate high-resolution textured models, capturing the pure ground geometry presents challenges. Especially in areas covered with vegetation, SfM / MVS technology has difficulty penetrating dense branches and leaves to obtain the true ground elevation. The generated model is often the surface of the vegetation canopy or contains a large amount of vegetation noise. Complex post-processing is required to try to extract surface information, and the accuracy is difficult to guarantee. At the same time, this technology is highly dependent on the texture features of the image. It is easy to encounter matching difficulties on surfaces with single or repeated textures (such as deserts, snowfields, and large areas of grassland), resulting in reconstruction failure or decreased accuracy. Changes in lighting conditions and the presence of shadows will also have an adverse effect on the geometric accuracy and texture consistency of the reconstruction results. Even if a dense point cloud can be generated, the subsequent surface reconstruction process may sacrifice some of the real micro-topography details in the pursuit of smoothness. Summary of the Invention

[0004] The main purpose of this invention is to provide a three-dimensional land parcel rendering method based on the integration of GIS and digital twins. It can effectively utilize the existing GIS data foundation, and by introducing fractal geometry principles and adaptive iterative enhancement mechanisms, it can intelligently generate rich microscopic details that far exceed the resolution of the original data and conform to the laws of natural landforms, significantly improving the realism and detail level of the three-dimensional land parcel model, and solving the problems of insufficient details and lack of natural complexity in traditional methods, thereby being able to construct a higher-fidelity digital twin model to serve various applications that require fine terrain expression.

[0005] In order to solve the above problems, the technical solution of the present invention is achieved as follows:

[0006] A three-dimensional land parcel drawing method based on the integration of GIS and digital twins, the method comprising:

[0007] Step 1: Generate original digital twin 3D point cloud data of the target land area based on 3D point cloud using GIS data;

[0008] Step 2: Perform a K-nearest neighbor search or radius neighborhood search on the original digital twin 3D point cloud data to determine the neighborhood point set of each point. For each point and its neighborhood point set, analyze the relationship between the elevation change amplitude and horizontal distance at different scales and calculate the local terrain fractal index of each point.

[0009] Step 3: Determine the location strategy for generating new points in each iteration; define a fractal dimension-increasing operator that outputs an elevation adjustment value based on the current number of iterations and the local average fractal index around the point;

[0010] Step 4: Apply the fractal dimension increase operator and generate new points based on the original digital twin 3D point cloud data through multiple iterations. Adjust the elevation of all points according to the elevation adjustment value to complete the 3D plot drawing.

[0011] Furthermore, in step 2, the original digital twin 3D point cloud data P = {p i =(x i ,y i ,z i )|i=1,...,N} perform K-nearest neighbor search or radius neighborhood search to determine each point p i Neighborhood point set x i For point p i The X-axis coordinate of y i For point p i The Y-axis coordinate of z i For point p i The Z-axis coordinate; N is the total number of points; for each point p iand its neighborhood Analyze the relationship between the change of elevation z and horizontal distance at different scales; calculate the p of each point i The local terrain fractal index Ψ i , local terrain fractal index Ψ i The higher it is, the more complex the local terrain is.

[0012] Furthermore, in step 2, calculate each point p i The local terrain fractal index Ψ i for:

[0013]

[0014] Among them, z j For point p i The neighborhood point p of j The Z-axis coordinate of point p j Elevation value; N i For point p i The number of neighboring points of min is the minimum elevation resolution, avoiding zero or negative values inside the logarithm; For point p i The average horizontal distance with its neighboring points; R neigh is the radius used to define the neighborhood search; λ is the weight coefficient for adjusting the volume-area ratio, ranging from 0.1 to 0.5; V hull( i) is point p i and its neighboring points The volume of the three-dimensional convex hull formed; A proj (i) is point p i and its neighboring points The area of the projection area on the XY plane; A plot is the total projected area of the original digital twin 3D point cloud data.

[0015] The present invention's 3D land parcel rendering method, based on the integration of GIS and digital twins, has the following beneficial effects: It effectively utilizes widely available and relatively low-cost geographic information system data, such as digital elevation models (DEMs) or existing 3D point clouds, as the foundation for building digital twins. This reduces reliance on expensive, high-precision, high-density raw data, improving the accessibility and affordability of the technology. Secondly, by incorporating the calculation and analysis of local terrain fractal indices from the raw point cloud data, the method can quantitatively "perceive" and "understand" the inherent complexity and irregularities of the surface at different locations. This deep exploration of the original geometric information transcends the limitations of traditional methods that focus solely on coordinate accuracy. Based on this understanding of local complexity, the core fractal dimension-increasing operator enables intelligent and adaptive elevation adjustments: in areas with complex and drastically changing terrain, more significant adjustments are applied to enhance their features; in areas with flatter, simpler terrain, the terrain remains stable or undergoes minor adjustments. This differentiated processing allows for more targeted added details and avoids the unnatural effect of uniform or random perturbations across the entire model. More importantly, by iteratively generating new points and combining them with fractal dimensionality adjustment, this method can generate rich details and textures that far exceed the resolution of the original data and conform to the fractal characteristics of natural landforms, greatly enhancing the visual realism and immersion of the final three-dimensional land model. The generated model is no longer a simple geometric body that is overly smooth or angular, but instead presents a multi-scale, self-similar complex morphology similar to the surface of the real world. In addition, by considering multiple factors such as the iterative process, boundary effects, and the overall scale and slope of the land in the adjustment process, the global coordination and rationality of the generated details are ensured. The introduction of a neighborhood smoothing mechanism and an iterative weight decay strategy further ensures the stable convergence of the iterative process and the quality of the final model. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A schematic diagram of the method flow of a three-dimensional land parcel drawing method based on the fusion of GIS and digital twins provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0017] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0018] Example 1, reference Figure 1 : A three-dimensional land parcel drawing method based on the integration of GIS and digital twins, the method comprising:

[0019] Step 1: Generate original digital twin 3D point cloud data of the target land area based on 3D point cloud using GIS data;

[0020] When implementing a 3D land parcel rendering method based on the integration of a geographic information system (GIS) and digital twins, the initial step is crucial. Its core task is to utilize existing GIS data to construct a basic, original digital twin model in the form of a 3D point cloud for the target land parcel area to be 3D rendered. This process transforms abstract and potentially dispersed geospatial information into a concrete, centralized 3D geometric representation, laying a solid data foundation for subsequent fractal-based terrain detail enhancement and refined rendering.

[0021] Consider a specific application scenario: a mountainous area of approximately eight hectares located on the outskirts of a city needs to be mapped in a detailed three-dimensional plot for the purpose of ecological restoration planning and visualization simulation. The terrain of this plot is highly undulating, with gentle slopes, steep slopes, and small streams. Traditional two-dimensional maps or crude three-dimensional models cannot meet the accuracy requirements of planning and design. At this point, the first step of the three-dimensional plot mapping method is initiated. First, various types of GIS data related to the eight-hectare plot need to be collected. Data collection may come from various sources, including but not limited to official databases of local natural resources departments, results of past surveying and mapping projects, commercial data providers, or self-acquisition through modern methods such as drone aerial surveys. For this mountainous parcel, the types of data that may need to be obtained include: precise parcel boundary data, which is usually stored in a vector format (such as Shapefile) and clearly defines the legal or physical boundaries of the parcel. Its coordinates may be derived from cadastral surveys or high-precision Global Navigation Satellite System (GNSS) measurements; elevation data, which is key to generating three-dimensional form, can be derived from digital elevation models (DEMs) or digital surface models (DSMs). These models store elevation information of each surface point in raster form. Their resolution determines the fineness of the initial model. For example, five-meter resolution DEM data covering the area may be obtained. In addition, high-precision laser detection and ranging (LiDAR) scanning data is also an excellent choice. It directly records dense three-dimensional coordinate points on the surface in the form of a point cloud, which can provide richer terrain details. Assume that LiDAR data with an average point spacing of 0.5 meters is obtained for this project. At the same time, an orthophoto map (DOM) of the area may also be collected. Although it is mainly used for texture mapping or reference, it can also assist in determining terrain features and boundary accuracy in the initial stages.

[0022] After acquiring the data, the next key step is data preprocessing and integration. Since the data may come from different periods, different technical means, and different coordinate systems, they must be standardized. It is necessary to check whether the coordinate reference system (CRS) of all data is consistent. If inconsistent, the coordinates must be converted through professional GIS software (such as ArcGIS or QGIS) to unify them into the standard coordinate system required by the project, such as the national 2000 coordinate system or a local independent coordinate system, and ensure that the elevation benchmark is unified. For the plot boundary vector data, it is necessary to carefully check to ensure that its topological structure is correct, the boundaries are closed, and there are no errors such as self-intersection. For DEM or DSM raster data, denoising may be required, missing values may be filled, or resampling may be performed as needed to match specific resolution requirements. For LiDAR point cloud data, filtering and classification are required, for example, to distinguish ground points, vegetation points, building points, etc. For terrain mapping, only ground point data is usually used to accurately reflect the actual surface undulations. In the case of this eight-hectare mountainous area, the processing process may include: converting the shapefile file of the plot boundary from the local coordinate system to the project's unified National 2000 coordinate system; checking the five-meter resolution DEM data to confirm that there are no abnormal elevation values; and classifying the 0.5-meter point spacing LiDAR data to extract the ground point cloud and remove any possible noise points.

[0023] After the data is ready, the core step of generating the original digital twin three-dimensional point cloud data begins. The essence of this step is to combine the two-dimensional boundary information and elevation information to generate a series of points with three-dimensional coordinates (X, Y, Z). The implementation method varies depending on the main elevation data source selected. If DEM data is used, the raster conversion method can be used. The GIS software will read the plot boundary vector file and determine which DEM grid cells fall completely or partially within the plot area. For each grid cell within the range, the plane coordinates (X, Y) of its center point are taken, and the grid value corresponding to the cell is read as its elevation value (Z), thereby generating a point. For example, using a five-meter resolution DEM, for an eight-hectare (80,000 square meter) plot, theoretically about 3,200 initial points (80,000 / (55)) can be generated, which form a regular grid-like initial point cloud. If processed LiDAR ground point cloud data is used, the process is relatively straightforward. Because LiDAR data is natively in point cloud format, it's simple to use the parcel boundary vector file as a clipping tool, performing spatial queries or clipping operations in the GIS software to accurately extract all ground points falling within the eight-hectare parcel. Because LiDAR data density is much higher than a DEM, with a point spacing of, for example, 0.5 meters, an eight-hectare parcel may contain up to 320,000 (80,000 / (0.50.5)) or even more ground points, forming an initial point cloud with a higher density and better representation of microtopography. In this case, given the high level of topographic detail required, the decision was made to use a filtered LiDAR ground point cloud as the basis. Using the spatial selection function of the GIS software and based on the parcel boundary Shapefile, approximately 300,000 ground points belonging to the eight-hectare parcel were filtered from the massive regional LiDAR dataset.

[0024] The generated point cloud data requires final inspection and organization. It's necessary to confirm that all points have correct coordinate units (for example, meters), and to remove or correct any outliers caused by data processing errors (such as coordinates outside the reasonable range or extremely abnormal elevation values). The point cloud's coverage is then checked to ensure it fully aligns with the parcel boundaries. Finally, the resulting point set representing the initial 3D shape of the target parcel, after screening, processing, and conversion, is exported and stored in a standard point cloud file format (such as LAS, LAZ, XYZ, or CSV). This file contains the 3D coordinates of tens of thousands or even millions of points, each representing a digital sample of a specific location on the parcel's surface. For example, the resulting file might be an XYZ text file containing approximately 300,000 points, with each line containing the X, Y, and Z coordinates of a point, accurate to the centimeter level. This completes the process of generating the original digital twin 3D point cloud data for the target parcel area based on GIS data. This output point cloud file, known as the "original digital twin 3D point cloud data," provides a preliminary but crucial digital representation of the eight-hectare mountainous surface. Although it may not be rich in details, especially in areas with low data source resolution or gentle terrain changes, it provides an indispensable geometric framework and data foundation for the subsequent iterative refinement using fractal dimensionality increase operators to generate a more realistic and detailed final three-dimensional land model.

[0025] Step 2: Perform a K-nearest neighbor search or radius neighborhood search on the original digital twin 3D point cloud data to determine the neighborhood point set of each point. For each point and its neighborhood point set, analyze the relationship between the elevation change amplitude and horizontal distance at different scales and calculate the local terrain fractal index of each point.

[0026] Step 3: Determine the location strategy for generating new points in each iteration; define a fractal dimension-increasing operator that outputs an elevation adjustment value based on the current number of iterations and the local average fractal index around the point;

[0027] Step 4: Apply the fractal dimension increase operator and generate new points based on the original digital twin 3D point cloud data through multiple iterations. Adjust the elevation of all points according to the elevation adjustment value to complete the 3D plot drawing.

[0028] The core task of step two is to conduct in-depth local geometric analysis of the raw point cloud data to quantify the terrain complexity at each location. Specifically, this requires iterating over every 3D point in the raw point cloud. For this eight-hectare mountainous area, this involves individually processing the approximately 300,000 ground points extracted from the LiDAR data. For each point, its local surroundings, namely its set of neighboring points, must be determined. This can be achieved using two common spatial search strategies: K-nearest neighbor search or radius neighborhood search. K-nearest neighbor search finds a fixed number of spatially closest neighbors for the current point, for example, the fifteen nearest points. This approach has the advantage of ensuring a sufficient number of neighbors for analysis even in areas with uneven point cloud density, but the actual spatial extent of the neighborhood may vary. Another strategy is radius neighborhood search, which searches for all points within a preset horizontal distance from the current point, for example, a radius of two meters. This approach defines a fixed neighborhood range and better reflects local point density, but may find few or no neighbors in sparse areas. Assume that in this case, given the high and relatively uniform density of the LiDAR point cloud, a two-meter radius neighborhood search strategy is chosen. Then, for each original point, the computer will use an efficient spatial index structure (such as a kd-tree or octree) to quickly retrieve all other points within a two-meter horizontal distance, forming the point's neighborhood point set.

[0029] After defining the neighborhood, the next key step is to analyze the geometric characteristics of each point and its neighboring points. The core objective is to explore how elevation varies with horizontal distance. For a central point and all its neighboring points, the absolute value of the elevation difference between them and their respective horizontal distances are calculated. By analyzing these data pairs, we can understand whether the terrain is flat (small elevation change) or steep (large elevation change) within the local area of the point, as well as the severity of this change. This is essentially observing the characteristics of terrain undulation at different scales (determined by the neighborhood range). Based on this analysis of the relationship between elevation change and horizontal distance, combined with other possible local geometric metrics, such as the ratio of the convex hull volume of the three-dimensional point cluster formed by the point and its neighbors to its projected area on the horizontal plane, a key numerical indicator is calculated for each point: the local terrain fractal index. This index aims to capture and quantify the inherent complexity or irregularity of the terrain around the point with a single value, a core characteristic of natural landforms. The higher the fractal index, the more broken, rugged, and detailed the terrain in the local area where the point is located. For example, in the case of steep slopes of mountainous land, exposed rock areas, or stream edges, the calculated fractal index of points in these places will be relatively high, possibly between 1.8 and 2.2. On the contrary, if a point is located on a relatively flat and open plateau or gentle slope, and the elevation changes around it are gentle, the calculated fractal index will be lower, possibly in the range of 1.1 to 1.4. After completing this step, the original 300,000 point cloud data is not just a collection of coordinates. Each point is assigned an attribute value that describes the complexity of its local environment. This point cloud annotated with the fractal index becomes the intelligent foundation for subsequent refined operations, reflecting the digital twin's progress from simple geometric replication to deep feature understanding.

[0030] Next, we proceed to step three, which involves setting the rules and tools for the subsequent iterative enhancement process. First, we need to determine how and where new points will be generated to increase the density of the model during the iteration process. A common and intuitive strategy is the "midpoint generation strategy," where new points are placed at the midpoint of the line connecting existing pairs of neighboring points. For example, a pair of neighboring points (within a two-meter radius) is randomly selected from the current point set, and their three-dimensional coordinates are averaged to obtain the spatial location of the new point. This strategy allows new points to naturally fill in the gaps between existing data points, particularly in areas where sampling density was relatively low or where further refinement is needed. This strategy aims to improve the continuity and detail of the model by densifying the interior of existing structures. Another key aspect of step three is the definition of the "fractal dimensionality increase operator." This operator embodies the intelligence of the entire method. It is essentially a mathematical function or set of rules that calculates the appropriate adjustment to be made to the elevation of each point (both the original and the newly generated points) in each iteration. This adjustment is not random, but is influenced by a combination of factors. It aims to simulate the formation process of natural terrain, making complex areas more complex and flat areas remain flat or produce subtle natural undulations.

[0031] When the operator calculates the elevation adjustment value, it mainly relies on the following input information: First, the local average fractal index around the point, which is the most critical input. It directly links the terrain complexity calculated in step 2 with the elevation adjustment. Generally, the area with a higher local average fractal index will have a larger elevation adjustment (which may be positive or negative) to enhance its terrain characteristics; second, the number of current iterations. The operator may be designed to produce a larger adjustment amount in the early iterations to shape the macro form, and gradually reduce the adjustment amount in the later iterations to fine-tune the details. This design helps the stability of the algorithm. Furthermore, the macroscopic characteristics of the entire plot may also be taken into account, such as the maximum elevation difference and average slope of the entire plot. This global information can be used to scale the adjustment value to match the overall landform scale of the plot. Furthermore, boundary effects must be considered. Elevation adjustments at points near the plot boundary may be moderately suppressed or smoothed to avoid unnatural sharp mutations or breaks at the model edge and ensure a smooth transition with the surrounding environment. Finally, the magnitude of the adjustment value may also be related to the degree of deviation of the fractal index at that point from the average fractal index of the entire plot, thereby highlighting the uniqueness of local features. This fractal dimensionality increase operator ultimately outputs a specific elevation adjustment value, such as positive 0.05 meters or negative 0.01 meters.

[0032] This value will be used in step 4 to guide the dynamic evolution of the point cloud elevation. After defining the new point generation strategy and fractal dimension increase operator, the final iterative rendering phase, step 4, begins. This phase is the core of the entire method. By repeatedly applying the rules defined in step 3, the original point cloud annotated with fractal indices is gradually evolved into the final refined 3D land parcel model. Initialization is performed, using the original point cloud generated in step 1 (approximately 300,000 points with fractal indices in this case) as the starting point for iteration 0. A total number of iterations is then set, such as five based on experience or the required model detail. Next, an iterative loop is entered. Taking the first iteration (k = 0) as an example, the process is roughly as follows: First, according to the preset new point generation strategy (midpoint generation), all neighboring point pairs in the current point set are traversed to calculate the horizontal coordinates (X, Y) of a large number of potential new points. For example, based on 300,000 points, hundreds of thousands of potential new point locations may be generated. In the second step, it is necessary to calculate or update the local average fractal index for all points participating in this iteration (including the original point and all newly generated potential point positions). For the original point, the value calculated in step 2 can be used or recalculated based on the current slightly changed neighborhood; for the newly generated midpoint, its initial fractal index can be estimated by interpolating its parent point or based on its neighborhood in the current point set. In the third step, for each original point and each newly generated position point, the fractal dimension increase operator defined in step 3 is applied. The position of the point, its local average fractal index, the current number of iterations (k = 0), and other information are input into the operator to calculate the corresponding elevation adjustment value.

[0033] Step 4: Update the existing elevation. For each point in the original point set, add the calculated elevation adjustment to its current elevation to obtain the point's new elevation at the end of the first iteration. More sophisticated adjustment mechanisms may be applied here, such as introducing an iteration weight factor to make adjustments larger in early iterations and smaller in later ones; or introducing a smoothing term for neighborhood elevation differences, so that adjustments are based not only on the complexity of the point itself but also on the relative elevations of surrounding neighbors, providing a degree of smoothing and coordination; or even taking into account morphological factors such as the overall perimeter-to-volume ratio of the plot to modify the adjustments. Step 5: Determine the elevation of the new point. For each newly generated midpoint location, an initial elevation must first be estimated. This can usually be obtained by simply averaging the elevations of its two parent points (after the elevations have been updated) or using more complex local interpolation methods. The calculated elevation adjustment for the new point location is then added to this initial estimate to obtain the final 3D coordinates of the new point. Step 6: Integrate the point set. All original points with updated elevations and all newly created points with fully calculated 3D coordinates are merged together to form a new point cloud with a larger number of points and a higher density. For example, after the first iteration, the number of points in the point cloud might increase from 300,000 to 450,000. This new point cloud serves as the input for the next iteration (k = 1). This complete iterative process (from generating new point locations to integrating the point cloud) is repeated a predetermined number of times (e.g., five).

[0034] Each new iteration builds upon the results of the previous iteration by performing neighborhood searches, calculating fractal indices, adjusting elevations, and adding new points. As the iterations progress, the point cloud density increases. Particularly in areas with complex terrain, the fractal dimensionality multiplication operator continues to operate, potentially making ridges more prominent and ravines deeper, or generating subtle, natural-looking undulations in flat areas. Due to the effects of the iterative weighting factors or the operator's inherent design, the magnitude of these adjustments gradually decreases, stabilizing the model and refining details. After completing the preset number of iterations (e.g., five), the resulting point cloud represents the final result of this fractal-based, iteratively enhanced 3D land parcel rendering method. For this eight-hectare mountainous parcel, the resulting point cloud may contain over one million points. Its geometry not only accurately reflects the macroscopic topography captured by the original LiDAR data but, more importantly, through the fractal dimensionality multiplication process, intelligently generates microscopic details consistent with natural landform patterns, rendering features such as sloping slopes, streams, and terraces more vivid and realistic.

[0035] Furthermore, in step 2, the original digital twin 3D point cloud data P = {p i =(x i ,y i ,z i)|i=1,...,N} perform K-nearest neighbor search or radius neighborhood search to determine each point p i Neighborhood point set x i For point p i The X-axis coordinate of y i For point p i The Y-axis coordinate of z i For point p i The Z-axis coordinate; N is the total number of points; for each point p i and its neighborhood Analyze the relationship between the change of elevation z and horizontal distance at different scales; calculate the p of each point i The local terrain fractal index Ψ i , local terrain fractal index Ψ i The higher it is, the more complex the local terrain is.

[0036] For example, the search radius can be set to two meters. Then, for a point, all other points with a horizontal distance of less than two meters from it are regarded as its neighbors. The neighborhood defined by this method has a clear physical spatial scale, which can better reflect the local point density and the actual terrain influence range, but the disadvantage is that in areas where the point cloud is very sparse, a point may have few or even no neighbors, while in dense areas there may be a large number of neighbors. The specific strategy to be adopted needs to be selected according to the characteristics of the point cloud data and application requirements. In actual operation, in order to efficiently complete the neighborhood search, especially for point cloud data containing hundreds of thousands or even millions of points, advanced spatial index data structures such as kd trees or octrees are usually used. These data structures can organize point cloud data so that the process of finding neighboring points does not have to perform inefficient global pairwise distance comparisons, but can quickly locate the spatial area that may contain neighbors, greatly improving computational efficiency.

[0037] For an eight-hectare mountainous plot, assuming a two-meter radius neighborhood search is chosen, the computer will perform a neighborhood query for each of the 300,000 points, efficiently identifying all neighboring points within a two-meter horizontal distance. The computer then stores the identifiers or coordinates of these neighboring points to form the point's neighborhood point set. This process generates a corresponding neighbor list for each point in the original point cloud. Once the neighborhood point set for each point is determined, the core analysis in step two begins: for each point and its identified neighborhood point set, the relationship between elevation change and horizontal distance is analyzed at different scales. The "different scales" here are implicit in the definition of the neighborhood, which itself defines the scale range for the local analysis. The goal of this analysis is to understand how the vertical variation (elevation change) of the terrain changes with horizontal distance within this local area. Specifically, for the central point and each of its neighboring points, the horizontal distance (based on their X and Y coordinates) and the absolute difference in their elevation values (Z coordinates) are calculated. By examining these two quantities (horizontal distance and elevation difference) for all neighboring points within this neighborhood, we can understand the terrain characteristics surrounding the central point. For example, if a large elevation difference appears within a small horizontal distance, it indicates that the terrain is steep or very uneven. Conversely, if the elevation difference relative to the central point is still small even near the neighborhood boundary (for example, within a radius of nearly two meters), it indicates that the terrain is relatively flat or the slope changes gently. This analysis essentially detects the vertical undulation of the local terrain and its regularity in changing with spatial distance.

[0038] Based on the detailed analysis of the relationship between elevation variation and horizontal distance, combined with possible other local geometric features, the ultimate goal is to calculate a single numerical value for each origin point—the local terrain fractal index. This index quantifies the complexity of the local terrain at that point. Rather than directly measuring slope or curvature, it attempts to capture the irregularity, fragmentation, or self-similarity of the terrain at the local scale—all manifestations of fractal characteristics in natural landforms. Generally speaking, a higher local terrain fractal index indicates a more complex, rugged, and detailed local terrain at that point. For example, a point's fractal index will be significantly higher in areas with numerous small undulations, steep slopes, gullies, or exposed rock. Conversely, a point on a smooth surface or a slope with a uniform gradient will have a relatively low fractal index. The calculation of the fractal index often draws on the concept of estimating fractal dimension in fractal geometry. It examines how the average or cumulative elevation difference within a neighborhood changes with horizontal horizontal distance or neighborhood scale (so-called scaling behavior) in logarithmic coordinates. This rate of change, or logarithmic slope, is closely related to the fractal index. The calculation may also incorporate other information, such as the ratio of the convex hull volume formed by the neighborhood point set in three-dimensional space to its projected area on a two-dimensional horizontal plane. This ratio can reflect the 3D morphological richness of the local area and further assist in determining terrain complexity. Furthermore, adjustment factors or reference values, such as the total projected area of the entire plot or the minimum elevation resolution, may be incorporated into the calculation to ensure robustness and comparability. Returning to the eight-hectare mountain example, a point located on a steep stream bank has a two-meter neighborhood encompassing both high bank points and low riverbed points. The elevation difference varies dramatically over a short distance, resulting in a calculated local terrain fractal index as high as 2.1. On the other hand, a point located in the center of a flat grassland has a neighborhood where all points have very similar elevation values, with little variation in elevation difference over distance. Therefore, the fractal index may be only around 1.2. By performing this series of neighborhood searches, geometric analysis, and fractal index calculations on all 300,000 points, the original point cloud data is given a new dimension. Each point now carries not only its spatial location information (X, Y, Z) but also a fractal index describing the complexity of its local environment. This point cloud dataset, with its fractal index attribute, acts like a detailed "terrain complexity map" for the original digital twin. It reveals the differences in topographical features across different regions within the plot, providing a critical, quantitative basis for adaptively adding details and adjusting elevations based on local complexity in subsequent steps.

[0039] Furthermore, in step 2, calculate each point p i The local terrain fractal index Ψ i for:

[0040]

[0041] Among them, z j For point p i The neighborhood point p of j The Z-axis coordinate of point p j Elevation value; N i For point p i The number of neighboring points of min is the minimum elevation resolution, avoiding zero or negative values inside the logarithm; For point p i The average horizontal distance to its neighboring points; r neigh is the radius used to define the neighborhood search; λ is the weight coefficient for adjusting the volume-area ratio, ranging from 0.1 to 0.5; V hull (i) is point p i and its neighboring points The volume of the three-dimensional convex hull formed; A proj (i) is point p i and its neighboring points The area of the projection area on the XY plane; A plot is the total projected area of the original digital twin 3D point cloud data.

[0042] The first part focuses on describing the scale effect of local vertical changes in terrain relative to horizontal distance, similar to the method commonly used in fractal geometry to characterize surface roughness. Specifically, the numerator of this part involves a logarithmic term, whose interior is the point p i All neighboring points p in its neighborhood N(i) j The average value of the absolute value of the elevation difference between the two. This average elevation difference, where N i is the number of neighborhood points, z i and z j are the elevation values of the center point and the neighboring points, respectively, which intuitively reflect the average vertical fluctuation amplitude of the terrain within the local small area defined by N(i). For example, in the eight-hectare mountain plot case mentioned earlier, if point p i Located on a flat terrace, then the neighboring points p within a radius of two meters around it j The elevation z j With z i are very close, the average elevation difference will be very small. i On the banks of steep streams, the neighborhood may contain points higher up the bank and lower down the bank, so the average elevation difference will be significantly increased. The formula takes the logarithm of this average elevation difference, on the one hand to compress the numerical range, and on the other hand because in fractal geometry analysis, many scale relationships are presented in the form of power laws, which appear as linear relationships in logarithmic coordinates. In addition, a log(δmin ) term, where δ min Represents a preset minimum elevation resolution. This parameter is introduced to avoid problems in logarithmic calculation when the terrain is extremely flat and the average elevation difference approaches zero (logarithmic zero or negative is undefined). It is equivalent to setting a lower limit of sensitivity to terrain changes, ensuring that valid calculation results can be obtained even in very flat areas. Looking at the denominator of the first part, it also contains logarithmic terms, involving and log(R neigh ).in It's point p i The average horizontal distance to all its neighboring points N(i). This value reflects the average distribution scale of neighboring points relative to the center point in the horizontal direction within a certain neighborhood. If the neighboring points are closely around the center point, then Smaller; if the neighbor points are more dispersed (but still within the neighborhood radius), then Larger. neigh It is the radius value used when performing neighborhood search (especially when using the radius neighborhood search strategy). In essence It measures the logarithmic ratio of the average distribution scale of points in the neighborhood to the maximum search scale. Dividing the numerator (the logarithmic measure reflecting the average elevation change) by the denominator (the logarithmic measure reflecting the horizontal scale change), the first term constructs a concept similar to the slope in fractal dimension estimation. It quantifies the rate at which the average vertical variation of local terrain changes with horizontal distance on a logarithmic scale. If this ratio is large, it means that the elevation has changed significantly even within a small range of horizontal distance changes, indicating steepness or high irregularity of the terrain, that is, higher local complexity. Conversely, if the ratio is small, it means that the terrain is relatively flat and less complex.

[0043] The second part of the formula supplements the description of local complexity from another perspective. It focuses on the three-dimensional morphological characteristics of the local point cloud cluster, that is, the relationship between volume and area. This term is given by ·log(A plot ) is composed of. Among them, V hull (i) refers to the central point p i The volume of the three-dimensional convex hull formed by and all its neighboring points N(i). The three-dimensional convex hull can be imagined as the smallest convex body formed by wrapping these points with an elastic membrane. This volume V hull (i) Measures the volume or thickness of this set of local points in three-dimensional space. proj (i) is the same set of points (p iThe area covered by the projection of N(i) and N(i)) on the two-dimensional horizontal (XY) plane represents the "occupancy range" of this set of points in the horizontal direction. Therefore, the ratio Provides a measure of the "height-to-thickness ratio" of the local terrain. If this ratio is large, it means that the points are widely distributed in the vertical direction relative to their horizontal coverage, possibly forming a local protrusion (such as a hillock, rock) or depression (such as a pothole), suggesting significant three-dimensional structural complexity. If the ratio is small, it means that the points are roughly distributed in a relatively flat area, close to a plane or gentle slope. In the example of eight hectares of mountain, if p i Located on the top of an isolated small mound, its neighborhood points are distributed on the slope of the mound, forming a convex hull volume V hull (i) relative to its horizontal projection area A proj (i) will be larger, resulting in a higher ratio. i In a flat area, the neighboring points are all near the same plane, the convex hull volume will be very small, and this ratio is close to zero. This term also includes log(A plot ), where A plot is the total projected area of the entire study plot (eight hectares). The purpose of multiplying its logarithm by the volume-to-area ratio is to relate the local three-dimensional morphological measure to the macroscopic scale of the entire plot, to perform a certain degree of normalization or scale adjustment, so that this morphological complexity measure has a certain degree of comparability between plots of different sizes. Finally, λ is a weight coefficient whose value range is limited to between 0.1 and 0.5. The role of this coefficient is to adjust the second part (the three-dimensional morphological complexity represented by the volume-to-area ratio) relative to the first part (roughness or scale effect complexity) in the final Ψ i The contribution of the λ function to the calculation. Users can adjust the value of λ based on the specific application scenario and the emphasis on terrain complexity. For example, if you are more concerned with the subtle undulations and textures of the terrain surface, you can set a smaller λ; if you are more concerned with local structural features such as three-dimensional protrusions or depressions, you can increase the λ value appropriately.

[0044] Further, click p i The average horizontal distance to its neighbors for:

[0045]

[0046] Among them, x j is the neighborhood point p j The X-axis coordinate of y j is the neighborhood point p j The Y-axis coordinate of .

[0047] Furthermore, in step 3, the location strategy for generating new points in each iteration is: the new point is generated at the midpoint of the line connecting the existing neighborhood points.

[0048] First, about point p i The average horizontal distance to its neighbors The purpose of calculation is to obtain a value that can represent point p i A single value representing the average distribution range or scale of neighboring points on the horizontal plane. In the local terrain fractal index Ψ discussed previously i plays a key role in the calculation formula, especially in the denominator, which is used to calculate the neighborhood search radius R. neigh Together they characterize the horizontal scale on which the analysis is based. The process starts with the point p i Determination of the neighborhood point set N(i) (this has been completed in the previous operation of step 2, for example, N is obtained by searching the neighborhood with a radius of two meters i neighbor points). For each neighbor point p in the neighborhood point set N(i) j , both need to be calculated between it and the center point p i The horizontal distance between the two points. This horizontal distance is calculated using the coordinates of the two points on the two-dimensional plane (i.e. the X-axis coordinate x i ,x j and the Y-axis coordinate y i ,y j ) is calculated using the standard Euclidean distance formula, specifically by calculating x i with x j The square of the difference plus y i with y j This calculation completely ignores the difference in elevation z between the two points. i ,z j , only focusing on their distance on the horizontal projection plane. This calculation process will be done for all N in N(i) i neighbor point p j Do this one by one and get N i Finally, add up all these calculated horizontal distance values and divide by the total number of neighbor points N i , I got some points i The average horizontal distance to its neighbors This arithmetic mean Summarizes the point p i The approximate average spacing or dispersion of neighboring points in the horizontal direction within the local microenvironment defined by N(i). For example, in the case of an eight-hectare mountain plot, if a point p iIn an area where LiDAR data is collected very densely, there may be dozens of neighboring points within a two-meter radius, and most of these neighboring points are distributed at a distance of p. i If the area is within one meter, then the calculated The value may be relatively small, such as 0.8 meters. On the contrary, if point p i In areas where the point cloud is relatively sparse, such as the edge of a plot or an obstructed area, there may be only a few neighboring points within a two-meter radius, and these points may be distributed relatively far away, close to the boundary of the two-meter radius, so the calculated The value may be larger, such as 1.6 meters. The value dynamically reflects the actual horizontal scale characteristics of the local neighborhood of each point, which is used as the fractal index Ψ i This helps make the complexity measure adapt to the changes in local point cloud density, thereby more accurately evaluating the roughness of the terrain at the corresponding scale.

[0049] Secondly, regarding the strategy for generating new points determined in step three, this involves how to add new points to the existing point cloud in the subsequent iterative process of step four to improve the level of detail of the model. This method explicitly adopts a specific strategy: new points are generated at the midpoint of the line connecting the existing neighborhood points. This strategy defines how new information is inserted into space. During each iteration, the algorithm identifies pairs of points in the current point cloud that are neighbors (the neighbor relationship here can follow the neighborhood relationship defined in step two, for example, points within two meters are neighbors, or based on other connection relationships as shown in the figure). For each pair of identified neighbor points, assuming that their coordinates are (x1, y1, z1) and (x2, y2, z2) respectively, the algorithm will calculate the midpoint of their line. The horizontal coordinate of the new point (x new ,y new ) is set to the average of the horizontal coordinates of the two neighboring points, that is, x new =(x1+x2) / 2 and y new =(y1+y2) / 2. This calculation determines the position of the new point on the two-dimensional plane. As for the height z of the new point new, it is determined in step four through a more complex mechanism. Typically, an initial estimate is obtained by interpolation based on the elevations of neighboring points, and then a fractal dimensionality increase operator is applied to adjust the estimate based on factors such as the local complexity of the location. However, the key to step three lies in determining the "location" where the new point is generated—the midpoint of the neighbor lines. This midpoint generation strategy was chosen for its advantages: it is an intuitive and easy-to-implement method for densifying point clouds, effectively inserting new information nodes between existing data points, thereby improving the density and coverage uniformity of the point cloud, especially in areas with large spacing between points. Because new points are generated based on existing neighbor relationships, this strategy tends to refine areas where data structure already exists, helping to maintain and enhance the continuity of the terrain. It also has a certain degree of adaptability, as densely populated areas naturally have more neighbor pairs, theoretically resulting in more new points, although this may require control in implementation to avoid over-densification. For the example of eight hectares of mountain land, suppose that in a certain iteration, point A with coordinates (150.2, 310.5, 88.4) and point B with coordinates (150.8, 310.1, 88.9) are identified as neighbors (the horizontal distance between them is less than two meters). Then, according to the midpoint strategy, a new point position will be generated with horizontal coordinates x = (150.2 + 150.8) / 2 = 150.5, y

[0050] = (310.5 + 310.1) / 2 = 310.3. During the entire iterative process, thousands of such new point positions will be calculated, forming the basis for the next round of point cloud.

[0051] Furthermore, in step 3, the fractal dimension-increasing operator for:

[0052]

[0053] The elevation adjustment values are:

[0054]

[0055] Where Δz (k+1) (x, y) is the elevation adjustment value generated by the k+1th iteration at the new point (x, y); x and y are the X-axis coordinates and Y-axis coordinates of the new point respectively; k is the current iteration number; is the local average fractal index at the new point (x, y); α is the basic scaling factor of the dimensionality increase effect, ranging from 0.5 to 0.9; Ψ base is the mean of all local terrain fractal indices; p is the power exponent of the fractal index, which controls the sensitivity of complexity to elevation adjustment and has a value range of 0.4 to 0.6; Z max and Z minare the maximum and minimum elevation values of the original digital twin 3D point cloud data; L diag is the diagonal length of the original digital twin 3D point cloud data; K max is the preset maximum number of iterations; d boundaty (x, y) is the horizontal distance from the new point (x, y) to the nearest boundary of the original digital twin 3D point cloud data; σ b is the boundary influence attenuation factor, which controls the weakening degree of the fractal characteristics of the boundary area, and its value range is 0.6 to 1; η is the overall slope influence factor of the plot, and its value range is 0.4 to 0.8; S plot is the average slope value of the original digital twin 3D point cloud data.

[0056] Specifically, the composition of this fractal dimensionality increase operator is obtained by multiplying six main factors. Each factor modulates the final elevation adjustment value Δz from a specific angle, so that the adjustment process can respond to a variety of information such as local terrain characteristics, overall landform background, iterative process stage, spatial boundary conditions, and global slope. The first factor is α, which is called the basic scaling factor of the dimensionality increase effect. Its value range is usually set between 0.5 and 0.9. This coefficient α is equivalent to a global "intensity" controller, which sets the basic adjustment amplitude of the entire fractal dimensionality increase process. The larger the α value, the larger the absolute value of all calculated elevation adjustment values Δz tends to be, and the additional terrain features generated will be more significant and prominent; conversely, a smaller α value will result in a relatively mild adjustment and more delicate details. Choosing the right α value is crucial to controlling the overall visual effect and level of detail of the final model.

[0057] The second factor is It is directly related to the complexity of the terrain. It is the local average fractal index evaluated at the point (x, y) where the adjustment value is being calculated in the current iteration step. It quantifies the terrain complexity of the area immediately surrounding the point. This index itself is calculated in step 2 or dynamically updated during the iteration. base is the average value of the local terrain fractal index of all points in the entire plot, representing the average complexity level of the entire plot, as a comparison benchmark. This difference reflects the degree of deviation of the local complexity of the current point from the overall average level of the plot. If the difference is positive, it means that the area where the point is located is more complex than the average level; if it is negative, it means that it is simpler. This difference is then subjected to a power exponent p (its value range is usually between 0.4 and 0.6). This exponent p controls the sensitivity of the elevation adjustment to complexity deviation, and because its value is less than 1, it introduces a nonlinear response relationship, which may make moderate complexity deviations have a relatively more significant effect than smaller deviations or very large deviations (depending on the p value and the sign of the difference). The core role of this factor is to achieve the adaptability of the adjustment: those areas that are identified as significantly more complex than the average ( Much larger than Ψ base ), will receive a larger absolute adjustment (which may be positive or negative, depending on other factors and the specific implementation, but the magnitude will be larger), thereby further enhancing its complexity; while those areas with complexity below average will have relatively smaller adjustments, which help maintain the smoothness of the terrain or make slight adjustments.

[0058] The third factor is It introduces the macro-scale information of the entire plot. max and Z min are the maximum and minimum elevation values of the entire plot recorded in the original digital twin 3D point cloud data, and their difference Z max ―Z min Represents the total vertical relief of the plot. diag is the diagonal length of the bounding box outside the original point cloud data, which can be regarded as a characteristic horizontal size of the plot. Therefore, this ratio Overall, it reflects the ratio of the vertical scale of the plot to its horizontal scale, and can be regarded as a measure of the overall slope of the plot or the severity of the terrain undulation. Including this factor in the operator means that the calculated elevation adjustment value Δz will be scaled according to the macro-geomorphological characteristics of the plot itself. For a plot with a large height difference and large terrain undulation (such as the eight-hectare mountain with a steep valley in the previous case), this ratio will be larger, allowing a relatively large elevation adjustment; while for a plot with flat terrain and a small overall height difference, this ratio will be smaller, which will correspondingly limit the amplitude of the elevation adjustment. This ensures that the size of the generated terrain details is coordinated with the overall natural scale of the plot.

[0059] The fourth factor is It introduces the time dimension of the iterative process. Here k is the current iteration number (counting from 0 or 1), K max is the total number of iterations set in advance. As k changes from 0 to K max ,ratio From 0 to 1, then It changes from 0 to π. The square of the sine function sin 2 (·) The value in this interval is from 0 (when k=0 or k=K max ) changes to 1 (when k=K max / 2), then back to 0. This means that this factor causes the elevation adjustment to be smaller at the beginning and end of the iteration process, while reaching its maximum in the middle. This design strategy helps achieve a smooth evolutionary process: small adjustments in the early stages avoid excessive impact on the original terrain, allowing the model to gradually adapt to the changes; the largest adjustments in the middle stages focus on shaping the main fractal features; and a further reduction in adjustments at the end facilitates model convergence and fine-tuning of details, avoiding unnecessary oscillations in the final stages.

[0060] The fifth factor is Used to handle boundary effects. boundary (x, y) refers to the nearest horizontal distance from the current calculation point (x, y) to the original data boundary of the entire plot (usually the convex hull of the point cloud or the defined region boundary). b Is a parameter called the boundary effect attenuation factor (range 0.6 to 1), which controls the decay rate or range of the boundary effect. Exponential decay function The characteristic is that when the point (x, y) is very far from the boundary, d boundary Very large, the exponential part is a very large negative number, and the value of the entire factor approaches 0; when the point (x, y) is very close to the boundary, d boundary Approaches 0, the exponential part approaches 0, and the value of the entire factor approaches exp(0) = 1. Strictly follow the formula From the above, we can see that it will be the largest (close to 1) near the boundary and the smallest (close to 0) far from the boundary. This seems to enhance the boundary characteristics rather than weaken them. If we assume that its purpose is to "reduce the influence of the boundary", then this factor should reduce the Δz adjustment amplitude of the points close to the boundary. This means that when d boundary When d is small, this factor should be close to 0 or less than 1, and when d boundary When ηS is large (i.e., the point is inside the plot), the factor should approach 1, allowing for full adjustment. In the case of the eight-hectare mountain area, regardless of the specific function form, the design intention is to ensure that at the edge of the plot, such as points within five meters of the property line, the calculated elevation adjustment value Δz will be significantly weakened by this factor, thereby avoiding the generation of unnatural artificial structures at the model boundary and ensuring that the refinement effect inside the model can smoothly transition to the boundary. The sixth factor is (1+ηS plot ), which takes into account the average slope of the entire plot. plotis the average slope value of the entire plot calculated on the original point cloud data. η is a slope influence factor (ranging from 0.4 to 0.8), which is used to adjust the influence of the average slope on the adjustment amplitude. The form of this factor indicates that when the overall average slope of the plot S plot When it is large, the value of this factor will be greater than 1, which will slightly increase the elevation adjustment amplitude Δz of the entire area. The overall steeper terrain can more naturally carry the details of the larger undulations. For the overall flatter plots, S plot Small, this factor is close to 1, and has little effect on the adjustment range.

[0061] Furthermore, step 4 specifically includes: initializing the point set to the original digital twin 3D point cloud data P (0) =P; proceed to K max The process of the kth iteration includes: determining the location of the new point (x′) according to the location strategy of the new point generated in each iteration new ,y′ new ); Calculate these new points and the point set P of the kth iteration (k) The local average fractal index of all points P (k) Each point in Calculate the elevation adjustment value Δz (k+1) (x i ,y i ) and update the elevation is the point of the kth iteration The Z-axis coordinate of the point The elevation of ; for each new point (x ′ new ,y′ new ), by P (k) Interpolate the neighboring points and calculate the initial elevation estimate z′ init , then apply the elevation adjustment value to get the Z coordinate of the new point: z′ new =z′ init +Δz (k+1) (x′ new ,y′ new ); the new point p′ new =(x′ new ,y′ new ,z′ new ) is added to the point set to form the point set P of the k+1th iteration (k+1) ;Δz (k +1) (x i ,y i ) is the new point (x i ,y i), the elevation adjustment value generated by the k+1th iteration.

[0062] Furthermore, adjust point p i The elevation is:

[0063]

[0064] Among them, w k is the iteration weight factor of the kth iteration; For point p i Apply fractal dimension-increasing operators; For point p i The local average fractal index at ; For point p i The neighborhood point set at the kth iteration; is the neighborhood point p j Elevation at the kth iteration; For point p i and p j The three-dimensional Euclidean distance at the kth iteration is equal to σ smooth is the smoothing scale factor affected by the neighborhood elevation difference, ranging from 0.6 to 0.8; C plot is the perimeter of the target plot area; V plot The estimated volume of the target plot area.

[0065] The original digital twin 3D point cloud data P generated by GIS data in step 1 is designated as the point set of the zeroth iteration, denoted as P (0) For the eight-hectare mountain plot case that has been discussed, this P (0) It is a set of about 300,000 LiDAR ground points. At the same time, a total number of iterations K needs to be set in advance. max , this value determines the depth and computational complexity of the refinement process. For example, you can set K max After setting the initial state and the total number of iterations, the core iteration loop is entered, which will execute K max Second-rate.

[0066] At each iteration, for example, the kth iteration (the value of k ranges from 0 to K max ―1), which includes a series of closely connected operations. First, according to the new point generation position strategy determined in step 3 (for example, generating a new point at the midpoint of the line connecting the existing neighborhood points), the positions of the new points that need to be generated in this iteration (x′) are determined. new ,y′ new ). This process will be based on the current point set P (k)To identify neighbor relationships and calculate midpoints, a large number of potential new point two-dimensional coordinates are generated. In the first iteration (k = 0), based on the original P (0) To generate these new positions. Then, a key preparation is to calculate or update the local average fractal index This step needs to cover all points participating in this iteration, including not only the current point set P (k) Every existing vertex in Also includes all newly determined potential point positions (x′ new ,y′ new ). Calculate or estimate these points This is to ensure that subsequent elevation adjustments accurately reflect the latest local terrain complexity assessment at each location. Possibly by interpolating the indices of its "parent" point (the pair of neighboring points that generated it) or based on its position in P (k) , is estimated from the expected neighborhood in .

[0067] Then, the iterative process enters the core link of elevation adjustment, which is divided into two parts: updating existing points and calculating new points. (k) Every existing vertex in in is its elevation at the beginning of the kth iteration, and an elevation adjustment value Δz needs to be calculated (k+1) (x i ,y i This adjustment value is obtained by applying the fractal dimension-increasing operator defined in step 3 and inputting the position of the point (x i ,y i ), the current number of iterations k and the local average fractal index of the point Calculate Δz (k+1) (x i ,y i ), add it to the current elevation to get the new elevation of the point at the end of this iteration This update process will be applied to P (k) For each newly generated position (x′ new ,y′ new ), then it is necessary to first determine its initial elevation estimate z′ init This is usually done by (k) The neighboring points in (such as the two parent points that generated it, use their updated elevation z (k+1) ) is interpolated to complete the process. The simplest method is to take the average of the parent point elevations. Get the initial elevation z′ initThen, the fractal dimension-increasing operator is also applied to calculate the value of the new position (x′ new ,y′ new ) elevation adjustment value Δz (k+1) (x′ new ,y′ new This adjustment is then added to the initial estimate to give the final elevation z′ of the new point. new =z′ init +Δz (k+1) (x′ new ,y′ new ). In this way, each new point obtains the complete (x′ new ,y′ new ,z′ new ) three-dimensional coordinates, becoming a valid new vertex p′ new .

[0068] The final step is integration. All the original points with updated elevations (now their elevations are ) and all newly created and calculated elevation points p′ new Merged together, we form the point set P of the k+1th iteration (k+1) Compared with P (k) , the new point set P (k+1) It contains more points, and the elevation of the points is also adjusted according to the fractal dimension increase rule, which can theoretically describe the terrain more finely. For example, in the case of an eight-hectare plot, after the first iteration, the point set may change from P to (0) The 300,000 points increased to P (1) 450,000 points, and the elevations of these points have changed according to factors such as their local complexity. (k+1) Then, as the input for the next iteration (k+1), the whole process from generating new point positions to integrating point sets is repeated. This cycle will continue until the preset K max As the number of iterations increases, the point cloud density continues to increase, and the terrain details gradually appear and become richer under the action of the fractal operator. At the same time, due to the design of the iterative attenuation factor that may be included in the operator, the adjustment process will tend to be stable.

[0069] In addition, this method also provides a more refined elevation adjustment calculation formula specifically for updating the existing vertex p i Elevation This formula introduces more adjustment items based on the basic adjustment value, aiming to obtain a smoother and more consistent adjustment effect. The core idea is that the new elevation Equal to old elevation Add a multi-modulated comprehensive adjustment term. This comprehensive adjustment term first includes an iterative weight factor w k This factor is usually designed to decrease as the number of iterations k increases (e.g. w k =(K max ―k) / K max ), which means that at the beginning of the iteration w k Close to 1, the adjustment effect is strong, which helps to quickly shape the main features; in the later stage of iteration, w k Close to 0, the adjustment effect is weakened, which helps the model converge and stabilize. Next is the basic output value of the fractal dimension increase operator It is still based on point p i The position, the current iteration number k and its local average fractal index The calculated core adjustment driving force. However, the output value of this is not used directly, but is modulated by a term that reflects the difference in neighborhood elevation. This modulator term is a weighted average in the form of a fraction. Its numerator is the calculation point p i All its neighbor points p at the kth iteration j ∈N (k) (i) Elevation difference between And multiply each elevation difference by a three-dimensional Euclidean distance between them Exponential decay weight All these weighted elevation differences are then summed. is the point p at the kth iteration i and p j The full three-dimensional space distance between smooth is a smoothing scale factor (ranging from 0.6 to 0.8) that controls how quickly the weight decays with distance. The denominator is the sum of all these exponentially decaying weights, used for normalization. The result of the entire fraction is p i A weighted average relative elevation relative to its neighbors. This term introduces a local smoothing or consensus mechanism: if a point's neighbors are generally higher than it, even if the operator itself indicates a downward adjustment, a positive value of this term may partially offset or mitigate this adjustment, and vice versa. It makes the adjustment of a point depend not only on its own complexity, but also on the elevation trend of the surrounding neighboring points, which helps to avoid overly isolated or abrupt adjustment results. Finally, the entire adjustment term is multiplied by a factor related to the overall shape of the plot Among them C plot is the total perimeter of the target plot area, V plot is its estimated total volume. It can be regarded as a characteristic length scale derived from the volume. This may reflect the shape complexity or ductility of the plot (for example, a larger perimeter relative to the volume may mean a more irregular or slender shape). Adding 1 to this ratio and using it as a multiplier means that for plots with more complex overall shapes, the magnitude of the elevation adjustment may be slightly exaggerated. This refined formula is further refined by introducing the iterative weight w k , neighborhood elevation difference smoothing term and global shape factor, the effect of the basic fractal dimensionality increase operator is more comprehensively modulated, aiming to generate a three-dimensional terrain model that is both rich in details and overall coordinated and smooth. In summary, step 4 uses K max The iterative cycle systematically performs a series of operations such as new point generation, complexity evaluation, elevation update (may use basic or refined formulas) and point set integration. It is a dynamic evolution process that transforms the static digital twin point cloud P originally constructed by GIS data into a (0) , gradually transformed into a final 3D point cloud with more points, richer details and more realistic shapes This completed the task of high-precision three-dimensional land plot drawing based on the fusion of GIS and digital twins.

[0070] Furthermore, the iteration weight factor w of the kth iteration is k for:

[0071] w k =(K max ―k) / K max .

[0072] Assume that the iteration number k starts counting from 0 to K max ―1 ends (a total of K max At the initial stage of the iteration process, that is, when k = 0, the value of the weight factor w0 is (K max ―0) / K max , the result is equal to 1. This means that in the first iteration, the elevation adjustment calculated by other terms (such as fractal dimension increase operator, neighborhood elevation difference term, etc.) will be applied to the elevation update of the point with its full amplitude. As the iteration proceeds, the value of k gradually increases, and the numerator (K max ―k) decreases linearly accordingly, and the denominator K max remains unchanged, so the weight factor w k The value of also decreases linearly. For example, halfway through the iteration, if k is approximately equal to K max / 2, then w k The value of is approximately (K max ―K max / 2) / K max = 0.5, the elevation adjustment is reduced to about half. When the iteration is nearing the end, for example, in the last iteration, k = K max-1, the weight factor The value of (K max ―(K max ―1)) / K max =1 / K max This is a positive number much smaller than 1 (unless K max itself is very small), indicating that in the final stage, the amplitude of the elevation adjustment has been significantly compressed. max Step (i.e. complete all K max After iterations), then will be equal to (K max ―K max ) / K max = 0, which means that the adjustment function stops completely. Therefore, the weight factor w k The role of is to make the "strength" of the elevation adjustment decay linearly from the maximum value (100%) at the beginning of the iteration to the minimum value (close to 0%) at the end of the iteration.

[0073] Introducing this weight factor w that decreases with the number of iterations k The principle and purpose of w is multifaceted. First, it helps to control the stability of the iteration process. In the early stages of iteration, the point cloud is relatively sparse and the main features of the terrain have not yet fully emerged. At this time, a large elevation adjustment (w k Close to 1) helps to quickly shape the skeleton and main undulations of the terrain based on the fractal index analysis results. However, if such a large adjustment range is always maintained, as the point cloud density increases and the details become richer, subsequent iterations may introduce excessive disturbances, causing the model to oscillate or have difficulty converging to a stable state. By gradually reducing w k , which ensures that in the later stages of iteration, when the model is already relatively refined, adjustments become more moderate, focusing on fine-tuning and local optimization, making it easier to reach a stable and detailed final state. Secondly, this design also conforms to the characteristics of many natural processes or optimization algorithms, that is, transitioning from coarse adjustments to fine adjustments. It can be compared to the process of an artist creating a sculpture: starting with a drastic removal of the redundant parts (corresponding to the early iterations, high w k ), and then gradually switch to small tools for fine-tuning (corresponding to later iterations, low w k ). This strategy helps you gradually add and refine details while maintaining the overall structure.

[0074] Taking the case of eight hectares of mountain land as an example, assuming that the total number of iterations K is set maxThe number of iterations is 5, and the iteration number k is from 0 to 4. Then in the first iteration (k = 0), w0 = (5-0) / 5 = 1.0, and the elevation adjustment is fully effective. It may produce significant elevation changes in complex areas such as streams and steep slopes according to the initially calculated fractal index. In the second iteration (k = 1), w1

[0075] = (5-1) / 5 = 0.8. All calculated adjustments are multiplied by 0.8 before application, resulting in a smaller adjustment than the first iteration. By the third iteration (k = 2), w2 = (5-2) / 5 = 0.6, further reducing the adjustment to 60%. By the fourth iteration (k = 3), w3 = (5-3) / 5 = 0.4, reducing the adjustment to 40%. By the final iteration (k = 4), w4 = (5-4) / 5 = 0.2, representing only 20% of the initial potential. This primarily involves very subtle corrections, helping to achieve harmony and stability at the detailed level of the entire model.

[0076] As described above, the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that the technical solutions described in the above embodiments can still be modified, or some of the technical features thereof can be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A three-dimensional plot drawing method based on the integration of GIS and digital twins, characterized by: The method comprises: Step 1: Generate original digital twin 3D point cloud data of the target land area based on 3D point cloud using GIS data; Step 2: Perform a K-nearest neighbor search or radius neighborhood search on the original digital twin 3D point cloud data to determine the neighborhood point set of each point. For each point and its neighborhood point set, analyze the relationship between the elevation change amplitude and horizontal distance at different scales and calculate the local terrain fractal index of each point. Step 3: Determine the location strategy for generating new points in each iteration; define a fractal dimension-increasing operator that outputs an elevation adjustment value based on the current number of iterations and the local average fractal index around the point; Step 4: Apply the fractal dimension increase operator and generate new points based on the original digital twin 3D point cloud data through multiple iterations. Adjust the elevation of all points according to the elevation adjustment value to complete the 3D plot drawing.

2. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 1, characterized in that: In step 2, the original digital twin 3D point cloud data P = {p i =(x i ,y i ,z i )|i=1,...,N} perform K-nearest neighbor search or radius neighborhood search to determine each point p i Neighborhood point set x i For point p i The X-axis coordinate of y i For point p i Y-axis coordinate; z i For point p i The Z-axis coordinate; N is the total number of points; for each point p i and its neighborhood Analyze the relationship between the change of elevation z and horizontal distance at different scales; calculate the p of each point i The local terrain fractal index Ψ i , local terrain fractal index Ψ i The higher it is, the more complex the local terrain is.

3. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 1, characterized in that: In step 2, calculate each point p i The local terrain fractal index Ψ i for: Among them, z j For point p i The neighborhood point p of the sum j The Z-axis coordinate of point p j Elevation value; N i For point p i The number of neighboring points of min is the minimum elevation resolution, avoiding zero or negative values inside the logarithm; For point p i The average horizontal distance with its neighboring points; R neigh is the radius used to define the neighborhood search; λ is the weight coefficient for adjusting the volume-area ratio, ranging from 0.1 to 0.5; V hull (i) is point p i and its neighboring points The volume of the three-dimensional convex hull; A proj (i) is point p i and its neighboring points The area of the projection area on the XY plane; A plot is the total projected area of the original digital twin 3D point cloud data.

4. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 3 is characterized in that: Click p i The average horizontal distance to its neighbors for: Among them, x j is the neighborhood point p j The X-axis coordinate of y j is the neighborhood point p j The Y-axis coordinate of .

5. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 4 is characterized in that: In step 3, the location strategy for generating new points in each iteration is: the new point is generated at the midpoint of the line connecting the existing neighborhood points.

6. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 5, characterized in that: In step 3, the fractal dimension-increasing operator for: The elevation adjustment values are: Where Δz (k+1) (x, y) is the elevation adjustment value generated by the k+1th iteration at the new point (x, y); x and y are the X-axis coordinates and Y-axis coordinates of the new point respectively; k is the current iteration number; is the local average fractal index at the new point (x, y); α is the basic scaling factor of the dimensionality increase effect, ranging from 0.5 to 0.9; Ψ base is the mean of all local terrain fractal indices; p is the power exponent of the fractal index, which controls the sensitivity of complexity to elevation adjustment and has a value range of 0.4 to 0.6; Z max and Z min are the maximum and minimum elevation values of the original digital twin 3D point cloud data respectively; l diag is the diagonal length of the original digital twin 3D point cloud data; K max is the preset maximum number of iterations; d boundary (x, y) is the horizontal distance from the new point (x, y) to the nearest boundary of the original digital twin 3D point cloud data; σ b is the boundary influence attenuation factor, which controls the weakening degree of the fractal characteristics of the boundary area, and its value range is 0.6 to 1; η is the overall slope influence factor of the plot, and its value range is 0.4 to 0.8; S plot is the average slope value of the original digital twin 3D point cloud data.

7. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 6, characterized in that: Step 4 specifically includes: initializing the point set to the original digital twin 3D point cloud data P (0) =P; proceed to K max The process of the kth iteration includes: determining the location of the new point (x′) according to the location strategy of the new point generated in each iteration new ,y′ new ); Calculate the new point and the point set P of the kth iteration (l) The local average fractal index of all points P (k) Each point in Calculate the elevation adjustment value Δz (k+1) (x i ,y i ) and update the elevation is the point of the kth iteration The Z-axis coordinate of the point The height of each new point (x′ new ,y′ new ), through P (k) Interpolate the neighboring points and calculate the initial elevation estimate z′ init , then apply the elevation adjustment value to get the Z coordinate of the new point: z′ new =z′ init +Δz (k+1) (x′ new ,y′ new ); the new point p′ new =(x′ new ,y′ new ,z′ new ) is added to the point set to form the point set P of the k+1th iteration (k+1) ;Δz (k+1) (x i ,y i ) is the new point (x i ,y i ), the elevation adjustment value generated by the k+1th iteration.

8. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 7, characterized in that: Adjustment point p i The elevation is: Among them, w k is the iteration weight factor of the kth iteration; For point p i Apply fractal dimension-increasing operators; For point p i The local average fractal index at ; For point p i The neighborhood point set at the kth iteration; is the neighborhood point p j Elevation at the kth iteration; For point p i and p j The three-dimensional Euclidean distance at the kth iteration is equal to σ smooth is the smoothing scale factor affected by the neighborhood elevation difference, ranging from 0.6 to 0.8; C plot is the perimeter of the target plot area; V plot The estimated volume of the target plot area.

9. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twin according to claim 8, characterized in that: Iteration weight factor w for the kth iteration k for: w k =(K max ―k) / K max 。

Images