A 3D land parcel mapping method based on the fusion of GIS and digital twins

Abstract

This invention relates to the field of 3D point cloud construction technology, and more specifically, to a 3D land parcel rendering method based on the fusion of GIS and digital twins. The method includes: Step 1: Generating original digital twin 3D point cloud data of a target land parcel region using GIS data; Step 2: Performing 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; Step 3: Determining the location strategy for generating new points in each iteration; defining a fractal dimension enhancement operator and outputting an elevation adjustment value; Step 4: Applying the fractal dimension enhancement operator and adjusting the elevation of all points according to the elevation adjustment value to complete the 3D land parcel rendering. This invention significantly improves the realism and detail level of the 3D land parcel model, solving the problems of insufficient detail and lack of natural complexity in traditional methods.

Description

Technical Field

[0001] This invention belongs to the field of 3D point cloud construction technology, specifically involving a 3D land parcel drawing method based on the fusion of GIS and digital twins. Background Technology

[0002] In modern societal development, accurate three-dimensional digital representation of landforms is of paramount importance and has a wide range of applications. Whether serving urban planning and management, farmland analysis in precision agriculture, soil and water conservation and environmental monitoring, civil engineering design and earthwork calculation, natural disaster simulation and risk assessment, or meeting the demands for realistic scenes in virtual reality, game development, and film production, high-precision, highly realistic three-dimensional landform models are indispensable foundational data support. With the rapid development of computer graphics, remote sensing technology, and geographic information science, the technical means for acquiring and constructing three-dimensional landform models are becoming increasingly diverse. However, new challenges also arise, particularly in efficiently generating specific landform models that both conform to macroscopic geographical facts and showcase the complex textures and details of natural landforms at a microscopic scale.

[0003] Models generated from traditional GIS data (especially low-to-medium resolution DEMs or contour lines) often lack detail and realism. Standard DEM data has limited resolution (for example, common SRTM data is 30 meters or 90 meters, and even national-level basic DEMs are mostly 5 meters or 10 meters), failing to capture subtle surface undulations, gullies, field ridges, and other geomorphic details. These details are crucial for applications such as hydrological analysis, soil erosion simulation, and precision agriculture management. Surfaces generated based on contour line interpolation are prone to unnatural smooth transitions in flat areas or areas with sparse contour lines, making it difficult to accurately represent the true curvature changes of the terrain, and may even result in "terraced" interpolation artifacts. These models may meet accuracy requirements on a macroscopic level, but when viewed up close or requiring detailed analysis, their overly smooth or regular appearance is far removed from the complexity of natural terrain. While photogrammetric methods can generate high-resolution texture models, they face challenges in capturing pure ground geometry. Especially in vegetated areas, SfM / MVS technology struggles to penetrate dense foliage to obtain accurate ground elevation. The generated models are often just the surface of the vegetation canopy or contain significant vegetation noise, requiring complex post-processing to attempt surface information extraction, and accuracy is difficult to guarantee. Furthermore, this technology is highly dependent on the texture features of the imagery, and matching difficulties easily arise on surfaces with simple or repetitive textures (such as deserts, snowfields, and large areas of grassland), leading to reconstruction failure or decreased accuracy. Changes in lighting conditions and the presence of shadows also negatively impact the geometric accuracy and texture consistency of the reconstruction results. Even if dense point clouds can be generated, subsequent surface reconstruction processes may sacrifice some realistic micro-topographic details in pursuit of smoothness. Summary of the Invention

[0004] The main objective of this invention is to provide a 3D land parcel rendering method based on the fusion of GIS and digital twins. This method can effectively utilize existing GIS data and, by introducing fractal geometry principles and an adaptive iterative enhancement mechanism, intelligently generate rich microscopic details that far exceed the resolution of the original data and conform to the laws of natural landforms. This significantly improves the realism and detail of the 3D land parcel model, solving the problems of insufficient detail and lack of natural complexity in traditional methods. As a result, it can construct digital twin models with higher fidelity, serving various applications that require detailed terrain representation.

[0005] To solve the above problems, the technical solution of the present invention is implemented as follows:

[0006] A method for drawing 3D land parcels based on the fusion of GIS and digital twins, the method comprising:

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

[0008] Step 2: Perform 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 and the 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 iteration number and the local average fractal exponent around the point;

[0010] Step 4: Apply the fractal dimension-enhancing 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 is processed. Perform K-nearest neighbor search or radius neighborhood search to determine each point neighborhood point set ; For point X-axis coordinate; For point Y-axis coordinate; For point Z-axis coordinate; The total number of points; for each point and its neighborhood point set Analyzing elevation at different scales The relationship between the magnitude of change and the horizontal distance; calculate the relationship for each point. Local topographic fractal index Local topographic fractal index The higher the value, the more complex the local terrain.

[0012] Furthermore, in step 2, each point is calculated. Local topographic fractal index for:

[0013] ;

[0014] in, For point neighborhood points The Z-axis coordinate represents the point. Elevation value; For point The number of neighboring points; To achieve the minimum elevation resolution, avoid zero or negative values ​​within the logarithm; For point The average distance to its neighboring points; Define the radius used when performing a neighborhood search; The weighting coefficient for adjusting the influence of the volume-area ratio is set between 0.1 and 0.5. For point and its neighborhood point set The volume of the formed three-dimensional convex hull; For point and its neighborhood point set exist The area of ​​the projected region on the plane; The total projected area of ​​the original digital twin 3D point cloud data.

[0015] The 3D land parcel mapping method based on the fusion of GIS and digital twins of this invention has the following beneficial effects: Firstly, this method effectively utilizes widely available and relatively low-cost geographic information system data, such as digital elevation models or existing 3D point clouds, as the foundation for constructing digital twins. This reduces the dependence on expensive, high-precision, high-density raw data, improving the accessibility and economy of the technology. Secondly, by introducing the calculation and analysis of the local terrain fractal index of the raw point cloud data, this method can quantitatively "perceive" and "understand" the inherent complexity and irregularity of the land surface at different locations. This represents a deep mining of the original geometric information, surpassing the limitations of traditional methods that only focus on coordinate accuracy. Based on this understanding of local complexity, the core fractal dimensionality enhancement operator can perform intelligent and adaptive elevation adjustments: in areas with complex and rapidly changing terrain, more significant adjustments are applied to enhance its features; in areas with gentle and simple terrain, stability is maintained or minor adjustments are made. This differentiated processing approach makes the added details more targeted, avoiding the unnatural feeling caused by uniform or random perturbations of the entire model. More importantly, this method, through iterative generation of new points combined with fractal dimensionality enhancement, 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 3D land parcel model. The generated model is no longer a simple geometry that is too smooth or angular, but rather presents a complex, multi-scale, self-similar morphology similar to the Earth's surface in the real world. Furthermore, by considering various factors such as the iteration process, boundary effects, overall parcel scale, and slope during the adjustment process, the global consistency and rationality of the generated details are ensured. The introduction of a neighborhood smoothing mechanism and an iteration weight decay strategy further guarantees the stable convergence of the iteration process and the quality of the final model. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the method flow for 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 Implementation

[0017] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0018] Example 1, Reference Figure 1 A method for drawing three-dimensional land parcels based on the fusion of GIS and digital twins, the method comprising:

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

[0020] When implementing a 3D land parcel mapping method based on the fusion of 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 of the target land parcel area requiring 3D mapping, existing in the form of a 3D point cloud. This process transforms abstract, potentially dispersed geospatial information into a concrete, centralized 3D geometric representation, laying a solid data foundation for subsequent terrain detail enhancement and refined mapping based on fractal theory.

[0021] Consider a specific application scenario: A detailed 3D site mapping is needed for an approximately eight-hectare mountainous area in the suburbs of a city, for ecological restoration planning and visualization simulation. The site has significant topographical variations, including gentle slopes, steep embankments, and streams. Traditional 2D maps or rough 3D models cannot meet the accuracy requirements of the planning and design. Therefore, the first step in implementing this 3D site mapping method is to collect various types of GIS data related to the eight-hectare site. Data collection sources may be diverse, including but not limited to official databases of local natural resources departments, past surveying project results, commercial data providers, or self-acquisition through modern methods such as UAV aerial surveying. For this mountainous site, the types of data that may need to be acquired mainly include: precise site boundary data, which is usually stored in vector format (such as Shapefile), clearly defining the legal or physical boundaries of the site, and whose coordinates may come from cadastral surveys or high-precision Global Navigation Satellite System (GNSS) measurements; elevation data, which is crucial for generating the 3D shape, and can come from digital elevation models (DEMs) or digital surface models (DSMs). These models store the elevation information of each point on the surface in raster form, and their resolution determines the level of detail of the initial model. For example, five-meter resolution DEM data covering the area may be acquired; in addition, high-precision LiDAR scan data is also an excellent choice, as it directly records dense 3D coordinate points on the surface in point cloud form, providing richer terrain details. Assuming that LiDAR data with an average point spacing of 0.5 meters is acquired for this project; at the same time, orthophoto maps (DOMs) of the area may also be collected, although they are mainly used for texture mapping or reference, but can also help determine the accuracy of terrain features and boundaries in the initial stage.

[0022] After acquiring the data, the next crucial step is data preprocessing and integration. Since the data may originate from different periods, using different technologies, and employing different coordinate systems, standardization is essential. It's necessary to check the consistency of the Coordinate Reference System (CRS) for all data. If inconsistencies exist, coordinate transformation must be performed using professional GIS software (such as ArcGIS or QGIS) to unify it to the standard coordinate system required for the project, such as the national 2000 coordinate system or a local independent coordinate system, ensuring a consistent elevation datum. For parcel boundary vector data, careful verification is required to ensure correct topology, closed boundaries, and the absence of self-intersections or other errors. For DEM or DSM raster data, denoising, missing value filling, or resampling to match specific resolution requirements may be necessary. For LiDAR point cloud data, filtering and classification are required, such as distinguishing between ground points, vegetation points, and building points. For terrain mapping, typically only ground point data is used to accurately reflect the true surface undulations. In this case of an eight-hectare mountainous area, the processing steps might include: converting the shapefile 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 the absence of abnormal elevation values; classifying the 0.5-meter point spacing LiDAR data, extracting the ground point cloud, and removing any possible noise points.

[0023] Once the data is ready, the core step of generating the original digital twin 3D 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, a raster-to-point method can be adopted. The GIS software will read the land parcel boundary vector file and determine which DEM raster cells fall completely or partially within the land parcel area. For each raster cell within the area, the plane coordinates (X,Y) of its center point are taken, and the raster 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 (eighty-thousand-square-meter) land parcel, approximately three thousand two hundred initial points (80,000 / (55)) can be generated, which constitute 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 inherently in point cloud format, the only need to utilize the land parcel boundary vector file as a clipping tool is to perform spatial queries or clipping operations in GIS software to accurately extract all ground points falling within the eight-hectare land parcel. Since LiDAR data density is much higher than DEM, for example, with a point spacing of 0.5 meters, the eight-hectare land parcel may contain as many as 320,000 (80,000 / (0.50.5)) or even more ground points, forming a higher-density initial point cloud that better reflects micro-topographic features. In this case, considering the high requirements for terrain detail, it was decided to use the filtered LiDAR ground point cloud as the basis. Using the spatial selection function of the GIS software, based on the land parcel boundary Shapefile, approximately 300,000 ground points belonging to the eight-hectare land parcel were selected from the vast regional LiDAR dataset.

[0024] The generated point cloud data still requires final checking and processing. It's necessary to confirm the correctness of the coordinate units for all points (e.g., uniformly in meters), and to identify any outliers caused by data processing errors (such as coordinates outside a reasonable range, extreme elevation anomalies, etc.), which should be removed or corrected. The coverage area of ​​the point cloud must be checked to ensure it perfectly matches the site boundaries. Finally, these filtered, processed, and transformed point sets, representing the initial three-dimensional morphology of the target site, are exported and stored in a standard point cloud file format (such as LAS, LAZ, XYZ, or CSV). This file contains the three-dimensional coordinates of tens of thousands or even millions of points, each point being a digital sample of the site surface at a specific location. For example, the final result might be an XYZ text file containing approximately 300,000 points, with each line containing the X, Y, and Z coordinates of one point, accurate to the centimeter level. This completes the process of generating the original digital twin three-dimensional point cloud data of the target site area based on GIS data. This output point cloud file, the "original digital twin three-dimensional point cloud data," is a preliminary but crucial digital mirror of the eight-hectare mountainous terrain. While it may not be rich enough in detail, especially in areas with low data source resolution or gentle terrain changes, it provides an indispensable geometric framework and data foundation for subsequent iterative refinement using fractal dimension-enhancing operators to generate more realistic and detailed final 3D plot models.

[0025] Step 2: Perform 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 and the 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 iteration number and the local average fractal exponent around the point;

[0027] Step 4: Apply the fractal dimension-enhancing 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 perform in-depth local geometric analysis on the raw point cloud data, quantifying the terrain complexity at each location. Specifically, this requires traversing every 3D point in the raw point cloud; in this eight-hectare mountain case, this means processing each of the approximately 300,000 ground points extracted from LiDAR data. For any given point, its surrounding local environment, i.e., its neighborhood set, needs to be determined. This can be achieved using two commonly used spatial search strategies: K-nearest neighbor search or radius neighborhood search. K-nearest neighbor search finds a fixed number of spatially nearest neighbors for the current point, for example, setting it to find the fifteen nearest neighbors. The advantage of this method is that it can ensure 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 finds all points within a preset radius that are horizontally less than the current point, for example, setting the radius to two meters. This method defines a fixed neighborhood spatial extent, which better reflects local point density information, but may find few or no neighbors in sparse areas. In this case, assuming the LiDAR point cloud has a high and relatively uniform density, a neighborhood search strategy with a two-meter radius is chosen. Then, for each original point, the computer will use an efficient spatial indexing structure (such as a kd-tree or octree) to quickly retrieve all other points within a two-meter horizontal distance, forming the neighborhood point set of that point.

[0029] After determining the neighborhood, the next crucial step is to perform geometric feature analysis on each point and its neighboring point set, focusing on how elevation changes with horizontal distance. For a central point and all its found neighbors, 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 tends to be gentle (small elevation change) or steep (large elevation change) within a small local area at that point, and the degree of this change. Essentially, this involves 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 index—the local terrain fractal index—is finally calculated for each point. This index aims to capture and quantify the inherent complexity or irregularity of the terrain around the point with a numerical value, which is a core feature of natural landforms. A higher fractal index typically indicates a more fragmented, rugged, and detailed local terrain. For example, in the case study, points on steep slopes, exposed rock areas, or stream edges in mountainous terrain would have relatively high fractal indices, possibly between 1.8 and 2.2. Conversely, if a point is located on a relatively flat and open plateau or gentle slope with gradual elevation changes, its fractal index would be lower, possibly ranging from 1.1 to 1.4. After this step, the original 300,000 point cloud data points are no longer just a collection of coordinates; each point is assigned an attribute value describing the complexity of its local environment. This point cloud with fractal index annotations becomes the intelligent foundation for subsequent refined operations, demonstrating the advancement of digital twins from simple geometric replication to deep feature understanding.

[0030] Step three then sets the rules and tools for subsequent iterative enhancement. First, it's necessary to determine how and where to generate new points to increase model density during iteration. A common and intuitive strategy is the "midpoint generation strategy," where new points are placed at the midpoint of the lines connecting existing neighboring point pairs. For example, arbitrarily select a pair of neighboring points (within a two-meter radius) in the current point set, calculate the average of their 3D coordinates, and obtain the spatial location of a new point. Using this strategy, new points naturally fill the gaps between existing data points, especially in areas with relatively low sampling density or requiring further refinement. This strategy aims to improve the continuity and detail of the model by densifying the interior of the existing structure. Another core aspect of step three is defining the "fractal dimensionality enhancement operator." This operator embodies the intelligence of the entire method; essentially, it's a mathematical function or a set of rules used to calculate how much the elevation of each point (including original and newly generated points) should be adjusted in each iteration. This adjustment is not arbitrary, but rather influenced by a combination of factors, aiming to simulate the formation process of natural terrain, making complex areas more complex, and keeping flat areas flat or creating subtle natural undulations.

[0031] When calculating the elevation adjustment value, this operator mainly relies on the following input information: First, the local average fractal index around the point, which is the most crucial input. It directly links the terrain complexity calculated in step two with the elevation adjustment. Generally, areas with higher local average fractal indices will receive larger elevation adjustments (which may be positive or negative) to enhance their terrain features. Second, the number of iterations. The operator may be designed to generate larger adjustments in early iterations to shape the macroscopic form, while gradually decreasing the adjustments in later iterations to refine the details. This design contributes to the stability of the algorithm. The system is designed to ensure convergence. Furthermore, the macroscopic characteristics of the entire site may also be considered, such as the maximum elevation difference and average slope. This global information can be used to scale the adjustment value to match the overall topographic scale of the site. Additionally, boundary effects need to be considered. Points near the site boundary may have their elevation adjustments moderately suppressed or smoothed to avoid unnatural sharp abrupt changes or breaks at the model edges, ensuring 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 site, to highlight the uniqueness of local features. This fractal dimensionality enhancement operator ultimately outputs a specific elevation adjustment value, such as +0.05 meters or -0.01 meters.

[0032] This value will be applied in step four to guide the dynamic evolution of the point cloud elevation. After defining the new point generation strategy and the fractal dimension enhancement operator, the final iterative rendering stage, i.e., step four, begins. This stage is the core of the entire method, gradually evolving the original point cloud with fractal index annotations into the final refined 3D land parcel model by repeatedly applying the rules defined in step three. Initialization is performed first, using the original point cloud generated in step one (approximately 300,000 points with fractal indices in this case) as the starting point for the zeroth iteration. Then, a total number of iterations is set, for example, five iterations based on experience or the requirements for model detail. Next, the iterative loop begins. 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, calculating 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. The second step involves calculating or updating the local average fractal index for all points participating in this iteration (including the original point and all newly generated potential point locations). For the original point, the value calculated in step two can be used, or it can be recalculated based on the slightly changed neighborhood. For newly generated midpoints, their initial fractal index can be estimated by interpolating their parent point or based on their neighborhood in the current point set. The third step applies the fractal dimension-increasing operator defined in step three to each original point and each newly generated location point. The operator is input with the point's location, its local average fractal index, and the current iteration number (k=0) to calculate the corresponding elevation adjustment value.

[0033] The fourth step is to update the existing elevations. For each point in the original point set, its current elevation value is added to the calculated elevation adjustment value to obtain the new elevation of that point at the end of the first iteration. More refined adjustment mechanisms may be applied here, such as introducing an iteration weight factor to make the adjustment larger in early iterations and smaller in later iterations; or introducing a smoothing term for neighborhood elevation differences to ensure that the adjustment is based not only on its own complexity but also on the relative elevations of its neighbors, thus playing a smoothing and coordinating role; or even considering morphological factors such as the overall perimeter and volume ratio of the plot to correct the adjustment. The fifth step is to determine the elevation of the new point. For each newly generated midpoint location, an initial elevation needs to be estimated first. This can usually be obtained by simply averaging the elevation values ​​of its two parent points (after the elevation update), or by using a more complex local interpolation method. Then, the elevation adjustment value calculated for this new point location is added to this initial estimated elevation to obtain the final three-dimensional coordinates of the new point. The sixth step is to integrate the point set. All original points that have undergone elevation updates are merged with all newly created points whose complete 3D coordinates have been calculated, forming a new point cloud with a larger number of points and higher density. For example, after the first iteration, the number of points in the point cloud may 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 merging the point set) is repeated a preset total number of times (e.g., five times).

[0034] In each new iteration, neighborhood search, fractal index calculation, elevation adjustment, and new point addition are performed based on the results of the previous round. As the iteration progresses, the point cloud density continuously increases, especially in areas with complex terrain. The fractal dimensionality enhancement operator continues to operate, potentially making ridges more prominent, ravines deeper, or generating subtle, natural-looking undulations in flat areas. Due to the effect of the iteration weight factor or the design of the operator itself, the adjustment range gradually decreases, causing the model to stabilize and the details to be refined. After completing all the preset number of iterations (e.g., after five iterations), the final point cloud is the final result of this 3D land parcel rendering method based on the fusion of GIS and digital twins and using fractal principles for iterative enhancement. For this eight-hectare mountainous land parcel case, the final point cloud may contain more than one million points. Its geometric shape not only accurately reflects the macroscopic terrain captured by the original LiDAR data, but more importantly, through the fractal dimensionality enhancement process, it intelligently generates microscopic details that conform to the laws of natural landforms, making features such as steep slopes, streams, and terraces more vivid and realistic.

[0035] Furthermore, in step 2, the original digital twin 3D point cloud data is processed. Perform K-nearest neighbor search or radius neighborhood search to determine each point neighborhood point set ; For point X-axis coordinate; For point Y-axis coordinate; For point Z-axis coordinate; The total number of points; for each point and its neighborhood point set Analyzing elevation at different scales The relationship between the magnitude of change and the horizontal distance; calculate the relationship for each point. Local topographic fractal index Local topographic fractal index The higher the value, the more complex the local terrain.

[0036] For example, a search radius of two meters can be set. Then, for a given point, all other points horizontally within two meters are considered its neighbors. This method defines a neighborhood with a clear physical spatial scale, better reflecting local point density and the actual terrain influence range. However, its drawback is that in very sparse point cloud areas, a point may have few or no neighbors, while in dense areas it may have many. The specific strategy used depends on the characteristics of the point cloud data and the application requirements. In practice, to efficiently complete neighborhood searches, especially for point cloud data containing hundreds of thousands or even millions of points, advanced spatial indexing data structures, such as kd-trees or octrees, are typically used. These data structures organize the point cloud data, allowing the search for neighboring points to quickly locate spatial regions that may contain neighbors, rather than performing inefficient global pairwise distance comparisons, thus greatly improving computational efficiency.

[0037] For an eight-hectare mountainous plot, assuming a neighborhood search with a two-meter radius is chosen, the computer will perform a neighborhood query once for each of the 300,000 points, efficiently finding all neighboring points within a two-meter horizontal distance and storing their identifiers or coordinates to form the neighborhood set of that point. This process generates a corresponding neighbor list for each point in the original point cloud. Once the neighborhood set for each point is determined, the core analysis in step two begins: for each point and its determined neighborhood set, the relationship between the magnitude of elevation change and horizontal distance is analyzed at different scales. The "different scales" here are implicit in the definition of the neighborhood; the neighborhood itself defines a local analytical scale. The purpose of the analysis is to understand how the vertical variation of the terrain (elevation change) changes with the extension of horizontal distance within this local area. Specifically, for the center point and each of its neighboring points, the absolute difference between their horizontal distance (based on their X and Y coordinates) and their elevation values ​​(Z coordinate) is calculated. By examining these two quantities (horizontal distance and elevation difference) of all neighboring points within this neighborhood, one can understand the topographical features around the central point. For example, if a large elevation difference occurs within a small horizontal distance, it indicates that the terrain is steep or very uneven; conversely, if even points near the neighborhood boundary (e.g., within a radius of about two meters) have a small elevation difference relative to the central point, it indicates that the terrain is relatively flat or has a gentle slope. This analysis essentially probes the degree of vertical undulation in the local terrain and the regularity of its variation with spatial distance.

[0038] Based on the detailed analysis of the relationship between elevation variation and horizontal distance, and combined with other possible local geometric features, the ultimate goal is to calculate a single numerical value for each original point—the local topographic fractal index. This index is a quantitative expression of the local topographic complexity at that point. It does not directly measure slope or curvature, but rather attempts to capture the irregularity, fragmentation, or self-similarity characteristics of the terrain at a local scale—all manifestations of the fractal features of natural landforms. Generally, a higher local topographic fractal index indicates a more complex, rugged, and detailed local topography. For example, in areas with numerous small undulations, steep slopes, gullies, or exposed rock, the fractal index of a point will be significantly higher. Conversely, on a smooth plane or a slope with a uniform gradient, the fractal index will be relatively lower. The calculation method for the fractal index often borrows the idea of ​​estimating the fractal dimension in fractal geometry. It examines how the average elevation difference or cumulative elevation difference within the neighborhood changes with the average horizontal distance or neighborhood scale in logarithmic coordinates (i.e., the so-called scale behavior). This rate of change, or logarithmic slope, is closely related to the fractal index. The calculation process may also integrate other information, such as the ratio of the volume of the convex hull formed by the neighborhood point set in three-dimensional space to its projected area on a two-dimensional horizontal plane. This ratio reflects the fullness of the three-dimensional morphology of a local area, further aiding in the assessment of terrain complexity. Furthermore, adjustment factors or reference values ​​may be introduced during the calculation, such as the total projected area of ​​the entire plot or the minimum elevation resolution, to ensure the robustness and comparability of the calculation. Returning to the example of the eight-hectare mountainous area, a point located on the steep bank of a stream has both higher points on the bank and lower points on the riverbed within its two-meter neighborhood. The elevation difference changes drastically over a short distance, and the calculated local terrain fractal index may be as high as 2.1. Conversely, a point located in the center of a gentle grassland has very similar elevation values ​​for all points in its neighborhood, with minimal elevation difference over distance; its fractal index may only be around 1.2. By completing this series of neighborhood searches, geometric analyses, and fractal index calculations for all 300,000 points, the original point cloud data is given a new dimension. Now, each point not only has its spatial location information (X, Y, Z), but also carries a fractal index describing the complexity of its local environment. This point cloud dataset with fractal index attributes is like drawing a detailed "terrain complexity map" for the original digital twin. It reveals the differences in geomorphic features in different areas within the plot, providing crucial and quantitative basis for adaptively adding details and adjusting elevations based on local complexity in subsequent steps.

[0039] Furthermore, in step 2, each point is calculated. Local topographic fractal index for:

[0040] ;

[0041] in, For point neighborhood points The Z-axis coordinate represents the point. Elevation value; For point The number of neighboring points; To achieve the minimum elevation resolution, avoid zero or negative values ​​within the logarithm; For point The average distance to its neighboring points; Define the radius used when performing a neighborhood search; The weighting coefficient for adjusting the influence of the volume-area ratio is set between 0.1 and 0.5. For point and its neighborhood point set The volume of the formed three-dimensional convex hull; For point and its neighborhood point set exist The area of ​​the projected region on the plane; 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 variations in terrain relative to horizontal distance, similar to methods commonly used in fractal geometry to characterize surface roughness. Specifically, the numerator of this part involves a logarithmic term, which contains points. Its neighboring regions All neighboring points The average of the absolute values ​​of the elevation differences between them. This average elevation difference, i.e. ,in It is the number of neighboring points. and These are the elevation values ​​of the center point and neighboring points, respectively, which intuitively reflect the elevation of the area formed by... This refers to the average vertical undulation of the terrain within this defined local area. For example, in the previously mentioned case of an eight-hectare mountainous plot, if the point... If a point is located on a flat plateau, then its neighboring points within a two-meter radius are... elevation and If they are very close, the average elevation difference will be very small. But if the points are... Located on a steep stream bank, its neighborhood may include points higher up on the bank and lower down below, resulting in a significantly larger average elevation difference. Taking the logarithm of this average elevation difference in the formula serves two purposes: firstly, to compress the numerical range, and secondly, because in fractal geometry analysis, many scale relationships exhibit power-law forms, while in logarithmic coordinates they become linear relationships. Furthermore, a factor is subtracted from the numerator. Item, here This represents a preset minimum elevation resolution. This parameter is introduced to avoid problems in logarithmic calculations when the terrain is extremely flat and the average elevation difference approaches zero (logarithms of zero or negative numbers are undefined). It essentially sets a lower limit for sensitivity to terrain changes, ensuring effective calculation results even in very flat areas. Looking at the denominator in the first part, it also contains a logarithmic term, involving... and .in It is a point To all its neighboring points The average horizontal distance. This value reflects the average distribution scale of neighboring points relative to the center point in the horizontal direction within a defined neighborhood. If neighboring points are closely surrounding the center point, then... Smaller; if neighboring points are more dispersed (but still within the neighborhood radius), then Relatively large. This is the radius value used when performing a neighborhood search (specifically when using a radius-based neighborhood search strategy). The denominator is... In essence This measure is the logarithmic ratio of the average distribution scale of points within a neighborhood to the maximum search scale. Dividing the numerator (a logarithmic measure reflecting the change in average elevation) by the denominator (a logarithmic measure reflecting the change in horizontal scale), the first term constructs a concept similar to the slope in fractal dimension estimation. It quantifies the rate at which the average vertical variability of local terrain changes with horizontal distance at a logarithmic scale. A large ratio indicates significant elevation changes even over small horizontal distances, suggesting steep or highly irregular terrain, i.e., high local complexity. Conversely, a small ratio indicates relatively flat terrain and low complexity.

[0043] The second part of the formula supplements the description of local complexity from another perspective, focusing on the three-dimensional morphological features of local point cloud clusters, namely the relationship between volume and area. This term is derived from... Composition. Among them, It refers to the point from the center. and all its neighboring points The volume of the three-dimensional convex hull formed by these points. A three-dimensional convex hull can be imagined as the smallest convex shape formed by wrapping these points with an elastic membrane. This volume... It measures the "volume" or "thickness" occupied by this set of local points in three-dimensional space. These are the same set of points ( and The area covered by the projection of a point onto a two-dimensional horizontal (XY) plane represents the "territory" of that set of points in the horizontal direction. Therefore, the ratio... This provides a measure of the "height-to-thickness ratio" of local terrain. A large ratio indicates that these points are widely distributed vertically relative to their horizontal coverage, potentially forming a localized protrusion (such as a hill or rock) or depression (such as a pit), suggesting significant three-dimensional structural complexity. A small ratio indicates that these points are roughly distributed within a relatively flat area, close to a plane or gentle slope. In the example of an eight-hectare mountain, if... Located atop an isolated small mound, its neighboring points are distributed across the slope of the mound, forming a convex hull volume. Relative to its horizontal projected area It will be relatively large, resulting in a higher ratio. And if... In flat regions, where neighboring points are all near the same plane, the convex hull volume will be very small, and this ratio will be close to zero. This item also includes... ,in This is the total projected area of ​​the entire study plot (eight hectares). Multiplying its logarithm by the volume-area ratio is to link the local three-dimensional morphological measurement with the macroscopic scale of the entire plot, performing a degree of normalization or scale adjustment to make this morphological complexity measure comparable across plots of different sizes. Finally, This is a weighting coefficient, with a value limited to between 0.1 and 0.5. This coefficient adjusts the weighting of the second part (the 3D morphological complexity represented by the volume-to-area ratio) relative to the first part (roughness or scale effect complexity) in the final result. The magnitude of the contribution in the calculation. Users can adjust this according to the specific application scenario and the emphasis on terrain complexity. The value can be set to a smaller value if you are more concerned with the subtle undulations and textures of the terrain surface. If more attention is paid to local three-dimensional protrusions or depressions, the size can be appropriately increased. The value of .

[0044] Furthermore, point Average distance from its neighboring points for:

[0045] ;

[0046] in, For neighborhood points X-axis coordinate; For neighborhood points The Y-axis coordinate.

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

[0048] First, regarding points Average distance from its neighboring points The purpose of the calculation is to obtain a point that can represent the point. A single numerical value representing the average distribution range or scale of surrounding neighboring points on a horizontal plane. The previously discussed local topographic fractal index It plays a crucial role in the calculation formula, especially in the denominator, and is used to calculate the neighborhood search radius. The level scale upon which the characterization analysis is based. Calculation. The process begins with a point neighborhood point set The determination of (this was completed in the preceding operation of step two, for example, by obtaining it through a neighborhood search with a radius of two meters) (Number of neighboring points). For the set of neighboring points... Each neighboring point in Both require calculation of its distance from the center point. The horizontal distance between the two points. This horizontal distance is calculated using the coordinates of the two points on a two-dimensional plane (i.e., the X-axis coordinates). and Y-axis coordinates It is calculated using the standard Euclidean distance formula; specifically, it involves calculating... and The square of the difference plus and The square of the difference is taken, and then the square root is calculated. This calculation completely ignores the elevation difference between the two points. We only consider their distances on the horizontal projection plane. This calculation process will affect... All of them Neighboring points Proceed one by one to obtain Each horizontal distance is calculated individually. Finally, all these calculated horizontal distance values ​​are summed and divided by the total number of neighboring points. And that's how you get some points. Average distance from its neighboring points This arithmetic mean Summarized at the point Surrounded by Within this defined local microenvironment, the approximate average interval or dispersion of neighboring points in the horizontal direction. For example, in the case of dealing with an eight-hectare mountainous plot, if a point... Located in an area where LiDAR data collection is very dense, its two-meter radius neighborhood may contain dozens of neighboring points, and most of these neighboring points are distributed at a distance of... Within one meter, then the calculated The value might be relatively small, such as 0.8 meters. Conversely, if the point... In areas with sparse point clouds, such as the edges of plots or shaded areas, only a few neighboring points may be found within a two-meter radius, and these points may be relatively far apart, close to the boundary of the two-meter radius. Therefore, the calculated... The value could be quite large, for example, 1.6 meters. The value dynamically reflects the actual horizontal scale characteristics of the local neighborhood of each point, and is used in the fractal index. In the calculation, it helps to make the measure of complexity adaptable to changes in local point cloud density, thereby more accurately assessing the roughness of the terrain at the corresponding scale.

[0049] Secondly, regarding the new point generation strategy determined in step three, this relates to how to add new points to the existing point cloud to improve the model's level of detail during the subsequent iterations in step four. This method explicitly employs a specific strategy: new points are generated at the midpoint of the lines connecting existing neighboring points. This strategy defines the spatial insertion method for new information. In each iteration, the algorithm identifies pairs of mutually neighboring points in the current point cloud (the neighbor relationships can follow the neighborhood relationships defined in step two, such as points within two meters being neighbors, or based on other connection relationships as shown in the figure). For each identified pair of neighboring points, assuming their coordinates are... and The algorithm will calculate the midpoint of the line connecting them. The horizontal coordinates of the new point... It is set to the average of the horizontal coordinates of the two neighboring points, that is... and This calculation determines the position of the new point on the two-dimensional plane. As for the elevation of the new point... In step four, a more complex mechanism is used to determine the location of the new point. Typically, an initial estimate is obtained by interpolation based on the elevation of neighboring points, and then a fractal dimensionality-enhancing operator is applied to adjust it based on factors such as the local complexity of that location. However, the key to step three lies in determining the "location" where the new point is generated—the midpoint of the neighboring connection. This midpoint generation strategy is chosen for several reasons: it is an intuitive and easy-to-implement point cloud encryption method that effectively inserts new information nodes between existing data points, thereby improving the density and coverage uniformity of the point cloud, especially in areas where the original point spacing was large. Since the new point is generated based on existing neighbor relationships, this strategy tends to refine areas with existing data structures, helping to maintain and enhance the continuity of the terrain representation. It also has a degree of adaptability, as dense point cloud areas naturally have more neighbor pairs, theoretically potentially generating more new points, although control may be needed in practice to avoid over-encryption. For the example of an eight-hectare mountainous area, 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 location will be generated with horizontal coordinates of... , Throughout the iteration process, tens of thousands of such new point locations will be calculated, forming the basis for the next round of point cloud.

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

[0051] ;

[0052] The elevation adjustment value is:

[0053] ;

[0054] in, To the new point Place, No. The elevation adjustment value generated in the next iteration; and These are the X-axis and Y-axis coordinates of the new point, respectively. This represents the current iteration number; For new points The local average fractal index at that location; The scaling factor is the basic scaling factor for the dimension-adding effect, and its value ranges from 0.5 to 0.9; This is the mean of all local terrain fractal indices; The power exponent of the fractal exponent influence controls the sensitivity of complexity to elevation adjustments, with a value ranging from 0.4 to 0.6. and These represent the maximum and minimum elevation values ​​of the original digital twin 3D point cloud data, respectively. The diagonal length of the original digital twin 3D point cloud data; The maximum number of iterations is preset. For new points The horizontal distance to the boundary of the nearest raw digital twin 3D point cloud data; This is the boundary influence attenuation factor, which controls the degree of weakening of the fractal characteristics in the boundary region, and its value ranges from 0.6 to 1; The overall slope influencing factor of the plot ranges from 0.4 to 0.8. This represents the average slope value of the original digital twin 3D point cloud data.

[0055] Specifically, this fractal dimension-enhancing operator is composed of the multiplication of six main factors. Each factor adjusts the final elevation value from a specific perspective. Modulation is performed to enable the adjustment process to respond to various information such as local terrain features, overall geomorphological background, iterative process stage, spatial boundary conditions, and global slope. The first factor is... This is the basic scaling factor, known as the dimensionality-enhancing effect, and its value is typically set between 0.5 and 0.9. This factor... It functions like a global "intensity" controller, setting the basic adjustment range for the entire fractal dimensionality enhancement process. The larger the value, the more calculated the elevation adjustment values ​​will be. The larger the absolute value of , the more significant and prominent the generated additional terrain features will be; conversely, the smaller the value, the more prominent the generated additional terrain features will be. A value of [value] will result in a relatively gentle adjustment, producing more subtle details. Choosing the appropriate [value]... The value is crucial for controlling the overall visual effect and level of detail of the final model.

[0056] The second factor is This is directly related to the complexity of the terrain. Here... In the current iteration step, for the point where the adjustment value is being calculated. The locally average fractal index, which quantifies the topographic complexity of the region immediately surrounding the point, is calculated in step two or dynamically updated during iteration. This represents the average local topographic fractal index of all points on the entire site, signifying the average complexity level of the site as a benchmark. Therefore, This difference reflects the degree to which the local complexity of the current point deviates from the average level of the entire site. A positive difference indicates that the area where the point is located is more complex than the average; a negative difference indicates that it is simpler. This difference is then applied to a power factor. (Its value typically ranges from 0.4 to 0.6). This index The sensitivity of elevation adjustments to complexity deviations was controlled, and because its value is less than 1, it introduces a non-linear response relationship, which may make moderate complexity deviations have a relatively more significant impact than small or very large deviations (depending on the specific situation). The sign of the value and the difference). The core function of this factor is to achieve adaptive adjustment: those regions identified as having significantly higher complexity than the average ( Much larger The regions with lower-than-average complexity will receive larger adjustments (which may be positive or negative, depending on other factors and the specific implementation, but the magnitude will be greater), thereby further enhancing their complex characteristics; while those regions with lower-than-average complexity will receive relatively smaller adjustments, which will help maintain the smoothness of the terrain or make minor adjustments.

[0057] The third factor is It incorporates macro-scale information about the entire land parcel. Among them... and These are the maximum and minimum elevation values ​​of the entire plot recorded in the original digital twin 3D point cloud data, and their difference is... It represents the total vertical undulation range of the land parcel. This is the diagonal length of the bounding box of the original point cloud data, which can be considered as a characteristic horizontal dimension of the land parcel. Therefore, this ratio Overall, it reflects the ratio of the vertical scale to the horizontal scale of the site, and can be regarded as a measure of the overall slope or the degree of topographic relief of the site. Incorporating this factor into the operator means that the calculated elevation adjustment value... The scale will be adjusted based on the macroscopic topographic features of the site. For a site with significant elevation differences and large topographic relief (such as the eight-hectare mountainous area with steep valleys in the previous case), this ratio will be larger, allowing for relatively large elevation adjustments; while for a flat site with minimal overall elevation differences, this ratio will be smaller, correspondingly limiting the extent of elevation adjustments. This ensures that the generated topographic details are consistent in size with the overall natural scale of the site.

[0058] The fourth factor is It introduces the time dimension of the iterative process. Here It is the current iteration number (counting from 0 or 1). This is the pre-set total number of iterations. With... From 0 to ,ratio From 0 to 1, then It changes from 0 to The square of the sine function The values ​​within this range start from 0 (when...) or (when) changes to 1 (when) The elevation adjustment factor (TED) is set at 0 and then returns to 0. This means that the TED will result in smaller elevation adjustments at the beginning and end of the iteration process, and larger adjustments in the middle stage. This design strategy helps to achieve a smooth evolution process: small adjustments in the early stage can avoid causing too drastic impacts on the original terrain, allowing the model to gradually adapt to changes; the largest adjustment in the middle stage focuses on shaping the main fractal features; and the adjustment decreases again in the final stage, which helps the model converge and fine-tune details, avoiding unnecessary oscillations in the final stage.

[0059] The fifth factor is It is used to handle boundary effects. Refers to the current calculation point The nearest horizontal distance to the original data boundary of the entire plot (usually the convex hull of the point cloud or the defined region boundary). This is a parameter called the boundary effect decay factor (ranging from 0.6 to 1), which controls the decay rate or range of the boundary effect. (Exponential decay function) The characteristic is that when point When the distance is very far from the boundary, It is very large, the exponential part is a very large negative number, and the value of the entire factor approaches 0; when point When very close to the boundary, The value of the factor approaches 0, the exponential part approaches 0, and the overall factor value approaches 0. Strictly follow the formula. It appears that the factor has the highest value (close to 1) near the boundary and the lowest value (close to 0) far from the boundary. This seems to enhance rather than weaken the boundary characteristics. If the goal is to "weaken the boundary effect," then this factor should result in a larger value near the boundary. The adjustment range has decreased. This means that when... When the value is small, this factor should approach 0 or be less than 1, while when... When the elevation is large (i.e., the point is located inside the plot), the factor should approach 1 to allow for full adjustment. In the case of an eight-hectare mountainous area, regardless of the specific function form, the design intent is to ensure that the calculated elevation adjustment value is consistent at points on the edge of the plot, such as within five meters of the boundary line. This factor significantly weakens the effect, thus preventing the generation of unnatural artificial structures at the model boundaries and ensuring a smooth transition of refinement effects within the model to the boundaries. The sixth factor is... It takes into account the overall average slope of the plot. It is the average slope value of the entire plot calculated from the original point cloud data. This is a slope influence factor (ranging from 0.4 to 0.8) used to adjust the degree of influence of the average slope on the adjustment range. The form of this factor indicates the overall average slope of the plot. When the value is large, the value of this factor will be greater than 1, thus slightly increasing the elevation adjustment range across the entire region. Generally steeper terrain is better able to naturally accommodate more undulating details. For generally flatter plots, The factor is relatively small, close to 1, and has little impact on the adjustment range.

[0060] Furthermore, step 4 specifically includes: initializing the point set as the original digital twin 3D point cloud data. ;conduct The nth iteration, where the nth The process of each iteration includes: determining the location of the new point based on the location strategy for generating the new point in each iteration. ; Calculate these new points and the first The point set of the next iteration Local average fractal index of all points ;right Each point in Calculate the elevation adjustment value And update the elevation ; For the first The point of the next iteration The Z-axis coordinate represents the point. Elevation; for each new point Through the Interpolate from nearby points to calculate the initial elevation estimate. Then, by applying the elevation adjustment value, the Z-axis coordinates of the new point are obtained: ; New point Add to the point set to form the first The point set of the next iteration ; To the new point Place, No. The elevation adjustment value generated in the next iteration.

[0061] Further adjustments The elevation is:

[0062] ;

[0063] in, For the first The iteration weight factor for the next iteration; For the point Apply fractal dimension-increasing operators; For point The local average fractal index at that location; For point In the The neighborhood point set at the next iteration; For neighborhood points In the Elevation of the next iteration; For point and In the The 3D Euclidean distance at the next iteration is equal to ; The smoothing scaling factor for the influence of neighborhood elevation differences ranges from 0.6 to 0.8. The perimeter of the target plot area; This is an estimated volume of the target plot area.

[0064] The original digital twin 3D point cloud data generated from GIS data in step one The set of points designated for the zeroth iteration is denoted as... Regarding the eight-hectare mountainous land parcel case that has been under discussion, this This initially consisted of a set of approximately 300,000 LiDAR ground points. Additionally, a total number of iterations needed to be pre-defined. This value determines the depth and computational cost of the refinement process; for example, it can be set... The iterations will be performed five times. After setting the initial state and the total number of iterations, the core iterative loop will begin, which will execute... Second-rate.

[0065] In each iteration, for example, the first iteration... The next iteration ( The value ranges from 0 to Each of these steps involves a series of closely linked operations. First, based on the new point generation strategy determined in step three (e.g., generating a new point at the midpoint of a line connecting existing neighboring points), the positions of the batch of points to be newly generated in this iteration are determined. This process will be based on the current set of points. This is used to identify neighbor relationships and calculate midpoints, thereby generating a large number of potential new point 2D coordinates. In the first iteration ( When it is based on the original This is used to generate these new locations. Next, a crucial preparatory step is to calculate or update the local average fractal index. This step needs to cover all points participating in this iteration, not just the current set of points. Each existing vertex in This also includes all newly identified potential point locations. Calculate or estimate these points. This is to ensure that subsequent elevation adjustments accurately reflect the latest local terrain complexity assessment at each location. For new points, its It could be done by interpolating the index of its "parent" point (the pair of neighboring points that generated it) or based on its position in the system. The expected neighborhood is used to estimate.

[0066] Subsequently, the iterative process enters the core stage of elevation adjustment, which consists of two parts: updating existing points and calculating new points. For Each existing vertex in ,in It is in the first At the start of the next iteration, an elevation adjustment value needs to be calculated. This adjustment value is obtained by applying the fractal dimension-enhancing operator defined in step three and inputting the location of the point. Current iteration number and the local average fractal index at that point Obtained from parameters such as... Calculated Then, it is added to the current elevation to obtain the new elevation of that point at the end of this iteration. This update process will be applied to At all points, the height of the entire point cloud is adjusted based on its local complexity and other control factors. For each newly generated location... Then it is necessary to first determine its initial elevation estimate. This is usually done by adjusting its position in... The nearest points in the data (such as the two parent points that generated it) are used to update the elevation. This is accomplished through interpolation; the simplest method is to take the average elevation of the parent point. This yields the initial elevation. Then, using the same fractal dimension-increasing operator, a value specifically for this new location was calculated. Elevation adjustment value Then, this adjustment value is added to the initial estimate to obtain the final elevation of the new point. In this way, each new point obtains a complete [database / structure]. The three-dimensional coordinates became a valid new vertex. .

[0067] The final step is consolidation. This involves integrating all the original points that have undergone elevation updates (their current elevations are...). ) and all newly created points with calculated elevations When combined, they form the first... The point set of the next iteration Compared to New point set It contains more points, and the elevations of these points are adjusted according to fractal dimensionality enhancement rules, theoretically enabling a more refined depiction of the terrain. For example, in the case of an eight-hectare plot, after the first iteration, the point set might... The number of points increased to 300,000. There are 450,000 points, and the elevations of these points have varied according to factors such as their local complexity. Then, as the next iteration (the... (This process is repeated) with input (times), from generating new point locations to integrating the point set. This loop continues until the preset time is completed. The iteration continues. As the number of iterations increases, the point cloud density increases, and the terrain details gradually become more apparent and richer under the action of the fractal operator. At the same time, due to the design of the iteration decay factor and other factors that may be included in the operator, the adjustment process tends to stabilize.

[0068] Furthermore, this method also provides a more refined elevation adjustment calculation formula, specifically designed for updating existing vertices. elevation This formula introduces more adjustment terms on top of the basic adjustment value, aiming to achieve a smoother adjustment effect that better conforms to the overall shape. Its core idea is that the new elevation... Equal to old elevation In addition, a comprehensive adjustment term that has undergone multiple modulations is added. This comprehensive adjustment term first includes an iterative weighting factor. This factor is typically designed to vary with the number of iterations. The increase and decrease (e.g.) This means that in the early stages of the iteration A value close to 1 indicates a strong adjustment effect, which helps to quickly shape the main features; however, in the later stages of iteration... As the value approaches zero, the adjustment effect weakens, which helps the model converge and stabilize. Next are the basic output values ​​of the fractal dimension-increasing operator. It is still based on points Position, current iteration number and its local average fractal index The calculated core adjustment driver. However, this output value is not used directly, but is modulated by a term reflecting the elevation differences in the neighborhood. This modulation term is a fractional weighted average. Its numerator is the calculated point... Rather than in the All neighboring points in the next iteration Elevation difference between And multiply each elevation difference by a factor based on the three-dimensional Euclidean distance between them. Exponentially decaying weights Then sum all these weighted elevation differences. Here It is in the Next iteration time point and The complete three-dimensional spatial distance between them, and It is a smoothing scaling factor (ranging from 0.6 to 0.8) that controls the rate at which the weights decay with distance. The denominator is the sum of all these exponentially decaying weights, used for normalization. The result of the entire score calculation is... The term is a weighted average relative elevation compared 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, the 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 trends of its surrounding neighbors, helping 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. .in It is the total perimeter of the target plot area. It is its estimated total volume. This ratio can be viewed as a feature length scale derived from the volume. This may reflect the complexity or extensibility of the site's shape (e.g., a larger perimeter relative to volume may indicate a more irregular or elongated shape). Adding 1 to this ratio and using it as a multiplier means that for sites with more complex overall shapes, the elevation adjustment may be slightly amplified. This refined formula introduces iterative weights. The effects of the basic fractal dimensionality enhancement operator are more comprehensively modulated by the neighborhood elevation difference smoothing term and the global shape factor, aiming to generate a 3D terrain model that is both rich in detail and overall harmonious and smooth. In summary, step four... The iterative loop systematically performs a series of operations, including new point generation, complexity assessment, elevation updating (possibly using basic or refined formulas), and point set integration. It is a dynamic evolutionary process that transforms the initially static digital twin point cloud constructed from GIS data into a unified digital twin. It is gradually transformed into a final 3D point cloud with more points, richer details, and more realistic shapes. This completed the task of high-precision 3D land parcel mapping based on the fusion of GIS and digital twins.

[0069] Furthermore, the first The iteration weight factor of the next iteration for:

[0070] .

[0071] Let the iteration number be... Start counting from 0 and continue until... End (Total progress) (nth iteration). At the beginning of the iteration process, i.e. when At that time, weighting factor The value is The result is 1. This means that in the first iteration, the elevation adjustment calculated using other terms (such as fractal dimensionality increase operators, neighborhood elevation difference terms, etc.) will be applied to the point's elevation update in its entirety. As the iteration progresses, The value gradually increases, the molecule The value of decreases linearly accordingly, while the denominator The weighting factor remains unchanged. The value of also decreases linearly. For example, if halfway through the iteration, Approximately equal to ,So The value is approximately At this point, the elevation adjustment range is reduced to approximately half. When the iteration is nearing its end, for example, in the final iteration... At that time, weighting factor The value is This is a positive number much less than 1 (unless...). The magnitude of the elevation adjustment is very small (in itself), indicating that the range of elevation adjustments has been significantly compressed in the final stage. If we assume the iterative process conceptually proceeds to the [number missing]th [stage missing]... Step (i.e., complete all) After the next iteration), then Will equal to This means that the adjustment effect has completely stopped. Therefore, this weighting factor... Its function is to linearly decrease the "strength" of elevation adjustment from the maximum value (100%) at the beginning of the iteration to the minimum value (close to 0%) at the end of the iteration.

[0072] Introducing this weighting factor that decreases with the number of iterations The principles and purposes of this approach are multifaceted. Firstly, it helps control the stability of the iterative 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, allowing for relatively large elevation adjustments (…). Approaching 1) helps to quickly shape the terrain's framework and main undulations based on fractal index analysis results. However, if such a large adjustment range is maintained indefinitely, subsequent iterations may introduce excessive perturbations as point cloud density and detail increase, leading to model oscillations or difficulty in converging to a stable state. By gradually reducing... This design ensures that in later iterations, when the model is already quite refined, adjustments become more gentle, focusing primarily on fine-tuning and local optimization, thus making it easier to reach a stable and well-coordinated final state. Secondly, this design aligns with the characteristics of many natural processes or optimization algorithms, namely the transition from coarse to fine adjustments. This can be likened to an artist creating a sculpture: initially, the excess is drastically removed (corresponding to early iterations, high precision, and fine-grained optimization). Then gradually switch to smaller tools for fine-tuning (corresponding to later iterations, low-end). This strategy helps to gradually add and refine details while maintaining the overall structure.

[0073] Taking an eight-hectare mountainous plot as an example, let's assume a total number of iterations is set. For 5 iterations, the iteration number is... From 0 to 4. Then in the first iteration ( )hour, The elevation adjustment fully takes effect and may produce significant elevation changes in complex areas such as streams and steep slopes, based on the initially calculated fractal index. In the second iteration ( )hour, All calculated adjustments are multiplied by 0.8 before application, resulting in a weaker adjustment magnitude compared to the first iteration. In the third iteration (… ), The adjustment range was further reduced to 60%. Fourth iteration ( ), The adjustment range was reduced to 40%. In the last iteration ( )hour, The adjustment range is only 20% of the initial potential, mainly making very minor corrections, which helps the entire model achieve harmony and stability at the level of detail.

[0074] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for drawing three-dimensional land parcels based on the fusion of GIS and digital twins, characterized in that, The method includes: Step 1: Generate the original digital twin 3D point cloud data of the target plot area based on 3D point cloud using GIS data; Step 2: Perform 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 and the 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 iteration number and the local average fractal exponent around the point; Step 4: Apply fractal dimension enhancement operators to generate new points based on the original digital twin 3D point cloud data through multiple iterations, and adjust the elevation of all points according to the elevation adjustment value to complete the 3D plot drawing; In step 3, the fractal dimension-increasing operator for: The elevation adjustment value is: in, To the new point Place, No. The elevation adjustment value generated in the next iteration; and These are the X-axis and Y-axis coordinates of the new point, respectively. This represents the current iteration number; For new points The local average fractal index at that location; The scaling factor is the basic scaling factor for the dimension-adding effect, and its value ranges from 0.5 to 0.9; This is the mean of all local terrain fractal indices; The power exponent of the fractal exponent influence controls the sensitivity of complexity to elevation adjustments, with a value ranging from 0.4 to 0.

6. and These represent the maximum and minimum elevation values ​​of the original digital twin 3D point cloud data, respectively. The diagonal length of the original digital twin 3D point cloud data; The maximum number of iterations is preset. For new points The horizontal distance to the boundary of the nearest raw digital twin 3D point cloud data; This is the boundary influence attenuation factor, which controls the degree of weakening of the fractal characteristics in the boundary region, and its value ranges from 0.6 to 1; The overall slope influencing factor of the plot ranges from 0.4 to 0.

8. This represents the average slope value of the original digital twin 3D point cloud data.

2. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twins as described in claim 1, characterized in that, In step 2, the original digital twin 3D point cloud data is processed. Perform K-nearest neighbor search or radius neighborhood search to determine each point neighborhood point set ; For point X-axis coordinate; For point Y-axis coordinate; For point Z-axis coordinate; The total number of points; for each point and its neighborhood point set Analyzing elevation at different scales The relationship between the magnitude of change and the horizontal distance; calculate the relationship for each point. Local topographic fractal index Local topographic fractal index The higher the value, the more complex the local terrain.

3. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twins as described in claim 1, characterized in that, In step 2, calculate each point Local topographic fractal index for: in, For point neighborhood points The Z-axis coordinate represents the point. Elevation value; For point Z-axis coordinate; For point The number of neighboring points; To achieve the minimum elevation resolution, avoid zero or negative values ​​within the logarithm; For point The average distance to its neighboring points; Define the radius used when performing a neighborhood search; The weighting coefficient for adjusting the influence of the volume-area ratio is set between 0.1 and 0.

5. For point and its neighborhood point set The volume of the formed three-dimensional convex hull; For point and its neighborhood point set exist The area of ​​the projected region on the plane; 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 twins as described in claim 3, characterized in that, point Average distance from its neighboring points for: in, For neighborhood points X-axis coordinate; For neighborhood points The Y-axis coordinate.

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

6. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twins as described in claim 5, characterized in that, Step 4 specifically includes: initializing the point set as the original digital twin 3D point cloud data. ;conduct The nth iteration, where the nth The process of each iteration includes: determining the location of the new point based on the location strategy for generating the new point in each iteration. ; Calculate the new point and the first The point set of the next iteration The local average fractal index of all points; for Each point in Calculate the elevation adjustment value And update the elevation ; For the first The point of the next iteration The Z-axis coordinate represents the point. Elevation; for each new point Through the Interpolate from nearby points to calculate the initial elevation estimate. Then, by applying the elevation adjustment value, the Z-axis coordinates of the new point are obtained: ; New point Add to the point set to form the first The point set of the next iteration ; To the new point Place, No. The elevation adjustment value generated in the next iteration.

7. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twins as described in claim 6, characterized in that, Adjustment points The elevation is: in, For the first The iteration weight factor for the next iteration; For the point Apply fractal dimension-increasing operators; For point The local average fractal index at that location; For point In the The neighborhood point set at the next iteration; For neighborhood points In the Elevation of the next iteration; For point and In the The 3D Euclidean distance at the next iteration is equal to ; The smoothing scaling factor for the influence of neighborhood elevation differences ranges from 0.6 to 0.

8. The perimeter of the target plot area; This is an estimated volume of the target plot area.

8. The three-dimensional land parcel drawing method based on the fusion of GIS and digital twins as described in claim 7, characterized in that, No. The iteration weight factor of the next iteration for: 。

Images