Automatic surveying and mapping system and method for urban and rural planning complex terrain based on image analysis

Abstract

The invention belongs to the technical field of surveying and mapping remote sensing, and discloses an automatic surveying and mapping system and method for urban and rural planning complex terrains based on image analysis. The method comprises the following steps of: acquiring multi-view heterogeneous image data of a target surveying and mapping area, performing layered decoupling processing and multi-view registration on an image, identifying spectral texture features and geometric three-dimensional features of a terrain, and performing complexity function partitioning on the surveying and mapping area according to the spectral texture features and the geometric three-dimensional features; the regional surveying and mapping demand degree is calculated by analyzing a terrain local complexity index and spatial distribution dispersion characteristics, so that a personalized sampling parameter range is screened out; constructing an optimized three-dimensional terrain model by fusing a terrain feature map, a point cloud reconstruction result and a geometric constraint attribute, predicting surveying and mapping precision, and performing comparative analysis on the predicted surveying and mapping precision and a back projection result to obtain a confidence evaluation result; the method is applied to a personalized urban and rural planning surveying and mapping scheme; according to the invention, accurate identification and personalized automatic surveying and mapping of the complex terrain are realized.

Description

Technical Field

[0001] This invention relates to the field of surveying and remote sensing technology, and more specifically, to an automatic surveying system and method for complex terrain in urban and rural planning based on image analysis. Background Technology

[0002] With the rapid development of remote sensing technology and computer vision, multi-source heterogeneous data fusion mapping technology has entered a stage of refined application. Currently, technologies such as high-resolution satellite remote sensing, UAV oblique photography, and lidar scanning can acquire terrain data with centimeter-level accuracy, multispectral imaging can achieve vegetation penetration analysis, 3D reconstruction algorithms support complex terrain modeling, and real-time positioning technology has achieved sub-meter-level spatial positioning accuracy. These technologies have demonstrated great value in fields such as urban planning, infrastructure construction, and environmental monitoring, but they also face challenges in mapping accuracy under complex terrain environments.

[0003] However, existing technologies generally lack end-to-end collaborative processing mechanisms for multi-source heterogeneous data, making it difficult to guarantee the mapping accuracy of complex terrain areas. Specifically: 1. Traditional surveying methods rely solely on the geometric features of a single data source, neglecting the complete processing chain of terrain data acquisition and model reconstruction.

[0004] 2. Existing technologies struggle to identify key stages of error accumulation, cannot determine the specific sources and propagation paths of accuracy loss, and are unable to achieve targeted accuracy optimization.

[0005] 3. Existing technologies are particularly vulnerable in multimodal data fusion scenarios. When optical images, radar data, and point cloud information need to be processed collaboratively to build a complete terrain model, the lack of cross-modal consistency verification mechanisms makes the fusion results prone to geometric inconsistencies and accuracy mismatches.

[0006] In view of this, the present invention proposes an automatic mapping system and method for complex terrain in urban and rural planning based on image analysis to solve the above problems. Summary of the Invention

[0007] To overcome the aforementioned deficiencies of the prior art and to achieve the above objectives, the present invention provides the following technical solution: An automatic mapping system for complex terrain in urban and rural planning based on image analysis includes: The data acquisition module is used to acquire multi-view heterogeneous image data of the target mapping area. By extracting the features of the terrain anchor points corresponding to each view in the multi-view heterogeneous images, and combining the spatial geometric constraint relationship of the terrain anchor points with the view occlusion compensation coefficient, multi-view registration is performed to obtain multi-view image data, which includes vertical aerial images, oblique photography images and satellite remote sensing images. The feature extraction module performs layered decoupling on the surveyed area based on the spectral texture and geometric 3D features of multi-view image data. The layered decoupling includes penetration separation of vegetation-covered areas, depth inference of shadowed areas, gradient correction of steep slope areas, and texture enhancement of exposed surface areas to obtain a terrain feature map. The terrain zoning module is used to obtain the local complexity index of terrain based on the spatial distribution dispersion and gradient mutation frequency of terrain features at each level within the terrain feature map; and to divide the survey area based on the clustering distribution characteristics of the local complexity index of terrain to obtain the complexity zoning results. The model building module is used to generate an initial terrain point cloud based on the complex regional results, and to reconstruct a three-dimensional terrain model by identifying hollow areas and noise outliers in the initial terrain point cloud and combining the curvature propagation characteristics and topological continuity constraints of the neighboring terrain. The model optimization module is used to identify areas with insufficient accuracy in the 3D terrain model by analyzing the back projection consistency of multi-view image data and the elevation deviation of terrain control points. Based on the elevation deviation analysis results, it obtains the spatial propagation path and cumulative effect of the deviation, dynamically adjusts the layered extraction parameters and reconstruction strategy, and obtains the final mapping model and its confidence distribution map.

[0008] Furthermore, the construction process of the terrain feature atlas includes: Based on satellite remote sensing images, vegetation-covered areas are penetrated and separated. By analyzing the spectral differences between the near-infrared reflectance and the visible light band of vegetation-covered areas in satellite remote sensing images, a vegetation density distribution map is constructed. Based on the correlation model between vegetation density and the angle of incidence, the elevation shift below the vegetation is estimated, and the surface elevation after vegetation penetration is obtained. Based on oblique photogrammetry images, depth reasoning is performed on the shadow occlusion area to extract the geometric contour of the shadow boundary and the solar elevation angle information, reconstruct the three-dimensional spatial location of the shadow projection source, and combine the continuity of the neighborhood terrain gradient and the directional consistency of the surface texture of the shadow occlusion area to use conditional generative adversarial network to reason about the hidden terrain features of the shadow occlusion area. Gradient correction of steep slope areas is performed based on oblique photography images. By calculating the angle between the local normal vector of the terrain surface and the imaging optical axis, the terrain gradient distortion area is identified, and the true slope of the steep slope area is restored based on the difference in gradient observation values ​​of the terrain gradient distortion area in images from different viewpoints. Texture enhancement of exposed surface layers based on vertical aerial images: Multi-scale wavelet decomposition of vertical aerial images is performed to separate the high-frequency detail components and low-frequency background components of surface texture. Adaptive gain adjustment of high-frequency detail components is performed based on the local consistency of texture direction gradient and edge response intensity to obtain texture-enhanced surface feature map. The surface elevation estimate, hidden terrain features, gradient correction results, and surface feature map are aligned to obtain the terrain feature map.

[0009] Furthermore, the process of obtaining the complexity partitioning results includes: The local complexity of the terrain feature map is quantified to generate a terrain local complexity index; Spatial autocorrelation analysis is performed on the local terrain complexity index. By calculating the Moran index and local spatial outliers of the complexity index of adjacent regions, high-value and low-value clustering regions of complexity are identified. Based on the spatial connectivity and area threshold of the clustering regions, preliminary partition boundaries are determined. The preliminary partition boundary is refined by extracting terrain fault lines and gradient abrupt change zones in the boundary neighborhood and performing morphological optimization of the boundary using the watershed algorithm to generate a refined partition boundary. Based on the fine partition boundaries, the target mapping area is divided into boundaries to obtain complexity partitioning results. The complexity partitioning results include flat baseline areas, regular undulating areas, complex steep slope areas, and extreme terrain areas, and each partition is assigned a mapping priority weight.

[0010] Furthermore, the reconstruction process of the three-dimensional terrain model includes: The complexity partitioning results are applied to a partitioned adaptive point cloud generation strategy, and an initial terrain point cloud is generated by stereo matching of multi-view image data. Hole detection and noise identification are performed on the initial terrain point cloud. By calculating the local density distribution and K-nearest neighbor distance variance of the initial terrain point cloud, low-density hole regions and high-variance noise outliers in the initial terrain point cloud are identified. Based on the area size and shape complexity of the holes, the holes are classified into regular holes and irregular holes. The regular voids are repaired using a surface interpolation method based on radial basis functions, while the irregular voids are repaired using an iterative repair method based on neighborhood terrain curvature propagation to gradually fill the void region until the surface smoothness constraint is met. The initial terrain point cloud after repair is reconstructed by triangulation and the curvature of the mesh surface is optimized. At the same time, non-manifold edges and self-intersecting surfaces in the mesh are detected, and a local mesh reconstruction strategy is used to eliminate topological anomalies and generate a three-dimensional terrain model.

[0011] Furthermore, the process of obtaining the final mapping model and its confidence distribution map includes: The three-dimensional terrain model is back-projected onto the image planes of each viewpoint, the Hausdorff distance between the projected contour and the actual image contour is calculated, and the difference areas between the model and the actual terrain are identified by comparing the structural similarity index between the virtual viewpoint image generated by the model and the real image. Local accuracy assessment is performed on the discrepancy areas by setting up virtual terrain control points within the discrepancy areas and calculating the root mean square error between the model elevation and the actual elevation of the terrain control points. The accuracy deficiency level of the discrepancy areas is determined, and the causes of the error are analyzed based on the spatial distribution pattern of the accuracy deficiency levels. Based on the error cause analysis results, an error spatial propagation model is constructed. By tracing the propagation path of the error from the initial point cloud to the 3D terrain model and calculating the error amplification coefficient of each processing stage along the path, the key stages of error accumulation are identified, and parameters are backtracked and adjusted for the key stages. The adjusted parameters are iteratively optimized by constructing a nonlinear optimization problem with mapping accuracy as the objective function and hierarchical extraction parameters and reconstruction strategy parameters as optimization variables. Particle swarm optimization algorithm is used for global optimization, and the final mapping model and elevation confidence distribution maps of each region are output.

[0012] Furthermore, the process of penetrating and separating vegetation-covered areas includes: Spectral analysis was performed on the vegetation cover area in the multi-view image data. By calculating the combined features of the normalized vegetation index and the enhanced vegetation index, a preliminary vegetation cover classification map was generated. Morphological closing operation was performed on the classification map to obtain the vegetation cover classification map. The vegetation coverage classification map is quantified by density analysis. By analyzing the texture roughness and shadow ratio within the vegetation coverage area and combining the spectral mixing decomposition results of the vegetation, a vegetation density map is generated. Based on the vegetation density map, calculate the visibility probability of the ground surface under different vegetation densities; Based on the visibility probability, the vegetation-covered area is divided into high-visibility areas and low-visibility areas; Direct elevation observations are used for high visibility areas, while neighboring topographic trend extrapolation is used for low visibility areas. By weighted fusion of direct observations and extrapolated values, a surface elevation estimate after vegetation penetration is generated, and the elevation confidence level of each pixel is labeled.

[0013] Furthermore, the implementation process of local complexity metric includes: Multi-scale gradient calculation is performed on the terrain feature map to generate multi-scale gradient feature vectors; Perform curvature tensor analysis on the terrain feature map to generate curvature change entropy values; The terrain feature map is evaluated for directional heterogeneity by dividing the sliding window into eight fan-shaped sub-regions and calculating the average gradient direction of each sub-region. Based on the circumferential variance of the eight directional vectors, a directional heterogeneity index is generated. The norm, curvature change entropy, and directional heterogeneity index of the multi-scale gradient feature vector are normalized, and the comprehensive features are extracted by principal component analysis to determine the threshold for complexity classification.

[0014] Furthermore, the implementation process of the iterative repair method based on neighborhood terrain curvature propagation for the irregular voids includes: Boundary features are extracted from the irregular cavity. The three-dimensional coordinate sequence of the cavity boundary point cloud is identified, and curve fitting is performed on the cavity boundary points to calculate the curvature tensor and its principal direction at each boundary point. Based on the curvature tensor of the cavity boundary points, an initial field for curvature propagation is constructed by establishing a regular sampling grid inside the cavity and assigning initial curvature values ​​to the grid nodes. The initial field of curvature propagation is iteratively propagated and updated by calculating the curvature Laplacian operator of the grid nodes in each iteration and performing anisotropic diffusion along the principal curvature direction; Based on the curvature field after iterative convergence, a complete surface for the hole region is generated by a curvature-driven surface reconstruction method. The surface smoothness of the complete surface is evaluated. The total curvature variation of the surface is calculated. If the total curvature variation exceeds the preset variation threshold, the propagation parameters are backtracked and the iteration is repeated until a repair result that meets the smoothness constraint is generated.

[0015] Furthermore, the process of backtracking and adjusting parameters in key processes includes: For areas with insufficient accuracy and low level of difference, the error source is traced. The processing records of the area in the three-dimensional model reconstruction, point cloud generation, feature extraction and image registration stages are traced in reverse. Intermediate result data of the area are extracted at each stage to establish a complete data flow graph from the original image to the three-dimensional terrain model. An error sensitivity analysis is performed on each processing node in the data flow graph, and the error propagation coefficient of each node is calculated. The error propagation coefficient is the ratio of the change in output error to the change in input error. Based on the error propagation coefficient, an error accumulation path tree is constructed. By calculating the cumulative error amplification factor of all possible paths from the initial input to the final output, and identifying the critical propagation path with the cumulative amplification factor exceeding the threshold, the node with the largest error propagation coefficient on the critical propagation path is determined as the critical link of error accumulation. For the critical link of error accumulation, a gradient descent optimization framework is constructed with the final model accuracy as the target and the critical link parameters as the optimization variables. The critical link parameters are iteratively adjusted until the global error converges, and the optimized parameter configuration and its corresponding error suppression effect analysis are output.

[0016] An automatic mapping method for complex terrain in urban and rural planning based on image analysis includes: Step S1: Obtain multi-view heterogeneous image data of the target mapping area. Extract the features of the terrain anchor points corresponding to each view image in the multi-view heterogeneous image, and combine the spatial geometric constraint relationship of the terrain anchor points with the view occlusion compensation coefficient to perform multi-view registration to obtain multi-view image data. The multi-view image data includes vertical aerial images, oblique photography images and satellite remote sensing images. Step S2: Based on the spectral texture and geometric 3D features of multi-view image data, the survey area is decoupled in layers. The decoupling includes penetration separation of vegetation-covered areas, depth inference of shadowed areas, gradient correction of steep slope areas, and texture enhancement of bare surface areas to obtain a terrain feature map. Step S3: Obtain the local complexity index of the terrain based on the spatial distribution dispersion and gradient mutation frequency of terrain features at each level within the terrain feature map; and divide the survey area according to the clustering distribution characteristics of the local complexity index to obtain the complexity partitioning results. Step S4: Generate an initial terrain point cloud based on the complex regionalization results, and reconstruct a three-dimensional terrain model by identifying hollow areas and noise outliers in the initial terrain point cloud, combined with the curvature propagation characteristics of the neighboring terrain and topological continuity constraints. Step S5: Based on the 3D terrain model, identify areas with insufficient accuracy in the 3D terrain model by analyzing the back projection consistency of multi-view image data and the elevation deviation of terrain control points. Based on the elevation deviation analysis results, obtain the spatial propagation path and cumulative effect of the deviation, dynamically adjust the layered extraction parameters and reconstruction strategy, and obtain the final mapping model and its confidence distribution map.

[0017] The technical effects and advantages of the image analysis-based automatic mapping system and method for complex terrain in urban and rural planning according to this invention are as follows: This invention establishes a full-link collaborative processing mechanism for multi-source heterogeneous data, effectively restoring complete surveying and mapping information of complex terrain. This enhanced accuracy improves the reliability of urban and rural planning and design, enabling precise acquisition of the true terrain of complex areas with vegetation cover, shadow occlusion, and steep slope deformation, thus optimizing planning schemes and clarifying design parameters, providing accurate data support for infrastructure construction and land development. This invention helps build a surveying accuracy assurance defense line, preventing engineering risks and design defects caused by inaccurate terrain data, and improving the technical level and credibility of the surveying and mapping industry. The cross-modal consistency verification mechanism of this invention effectively integrates carefully registered multi-view data sources, maintaining the accuracy and consistency of surveying and mapping results. This invention can promptly detect and correct error accumulation behavior during the surveying and mapping process. By reconstructing a trustworthy surveying and mapping link for complex terrain, this invention provides a technical foundation for building an accurate and reliable digital terrain environment. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the image analysis-based automatic mapping system for complex terrain in urban and rural planning according to the present invention. Figure 2 This is a schematic diagram of the automatic mapping method for complex terrain in urban and rural planning based on image analysis according to the present invention. Detailed Implementation

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

[0020] Example 1 Please see Figure 1 As shown in this embodiment, the automatic mapping system for complex terrain in urban and rural planning based on image analysis includes: The data acquisition module is used to acquire multi-view heterogeneous image data of the target mapping area under different observation conditions. The multi-view heterogeneous image data includes key information such as vertical aerial images, oblique photogrammetry images, and satellite remote sensing images, which are acquired in real time through a multi-sensor data acquisition interface. Vertical aerial images record the orthophoto geometric information and surface texture details of the mapping area, oblique photogrammetry images provide the three-dimensional geometric structure and side-view angle features of the terrain, and satellite remote sensing images reflect the large-scale surface cover distribution and spectral response characteristics. By extracting the spatial features of the terrain anchor points corresponding to each view image, and combining the geometric constraints of the terrain anchor points with the view occlusion compensation coefficient, multi-view registration is performed to ensure the integrity and spatial consistency of the data acquisition.

[0021] The feature extraction module performs layered decoupling processing on the surveyed area based on the spectral texture and geometric 3D features of multi-view image data. This layered decoupling process includes penetration separation of vegetation-covered areas, depth inference of shaded areas, gradient correction of steep slope areas, and texture enhancement of exposed surface areas. Vegetation penetration separation identifies vegetation cover density by analyzing band differences in multispectral data and estimates the true surface elevation beneath the vegetation. Shade depth inference reconstructs the 3D structural information of shaded areas using geometric projection relationships. Gradient correction geometrically corrects perspective distortion in steep slopes to restore the true slope. Texture enhancement extracts details and optimizes contrast in exposed surface areas at multiple scales. These layered processing results are then aligned and fused to generate a topographic feature map containing multi-level topographic information, providing rich feature descriptions for subsequent analysis.

[0022] The terrain zoning module quantifies the terrain complexity distribution characteristics of the surveyed area by analyzing the spatial distribution dispersion and gradient abrupt change frequency of terrain features at each level in the terrain feature atlas. This module first calculates the spatial variability index of terrain features, extracts the statistical characteristics of gradient changes and the entropy characteristics of curvature distribution, and generates a local terrain complexity index. Then, it identifies the clustering patterns of the complexity index through spatial clustering analysis, and determines reasonable zoning boundaries by combining terrain continuity constraints and regional area thresholds. Based on the complexity distribution characteristics, the surveyed area is divided into flat baseline areas, regular undulating areas, complex steep slope areas, and extreme terrain areas, forming a hierarchical complexity zoning result. Corresponding surveying strategy parameters are assigned to each zoning to ensure that appropriate processing methods are used for areas of different complexity.

[0023] The model building module is used to generate a 3D terrain model based on the complexity partitioning results using a partitioning adaptive strategy. This module first determines the corresponding point cloud sampling density based on the complexity level of each partition, and extracts the initial terrain point cloud from multi-view image data using a stereo matching algorithm. Then, it performs a quality assessment on the initial point cloud, identifying void regions and noise outliers, and selects appropriate repair strategies based on the geometric characteristics of the defect types. Regular voids are filled using surface interpolation, while irregular voids are iteratively repaired using the curvature propagation characteristics of the neighboring terrain, ensuring that the repair results meet the requirements for smoothness and continuity of the terrain surface. Finally, through triangular mesh reconstruction and topology optimization, a 3D terrain model with a complete geometric structure is generated, providing basic data for accuracy verification and model optimization.

[0024] The model optimization module is used to verify the accuracy of the 3D terrain model. It iteratively optimizes the model through error propagation analysis and parameter backtracking adjustment. This module identifies areas with insufficient accuracy in the 3D terrain model by checking the consistency of backprojection from multi-view image data and analyzing the elevation deviation of terrain control points. It quantifies the error distribution characteristics and confidence levels of each area. Based on the error analysis results, it constructs an error spatial propagation model, tracing the cumulative path of errors from data acquisition to model reconstruction, and identifying key links and sensitive parameters that amplify errors. For key links, a parameter backtracking optimization strategy is adopted. By constructing a nonlinear optimization problem with mapping accuracy as the objective function, it dynamically adjusts the hierarchical extraction parameters and reconstruction strategy parameters to achieve global optimization. Finally, it outputs the final mapping model after multiple rounds of iterative optimization and its corresponding confidence distribution map, ensuring the accuracy and reliability of the mapping results.

[0025] The modules mentioned above are connected via wired and / or wireless means to enable data transmission and collaborative processing between them.

[0026] In embodiments of the present invention, the detailed implementation steps of multi-view registration include: Terrain anchor point features are extracted from multi-view heterogeneous images to generate stable feature descriptors for each viewpoint. Terrain anchor point feature extraction identifies and describes key terrain points with significant geometric features in the images, providing a reliable matching benchmark for subsequent registration. The feature extraction process employs a fusion algorithm of Scale Invariant Feature Transform (SIFT) and Histogram of Oriented Gradients (HOG), identifying candidate anchor points through multi-scale spatial extremum detection, such as peaks and valleys of terrain undulations, rock edges, and feature boundaries. Specifically, a Gaussian pyramid is first constructed for each viewpoint image, and local extrema are detected at different scale levels as candidate anchor points. Then, low-quality points with weak edge responses are eliminated through principal curvature analysis, retaining stable anchor points with good positioning accuracy. Finally, a 128-dimensional feature descriptor is constructed for each anchor point, encoding the gradient distribution and texture information of the anchor point's neighborhood.

[0027] Based on the 3D coordinate information of terrain anchor points, spatial geometric constraints between anchor points are constructed, and the reliability weights of these constraints are calculated. Spatial geometric constraints are the core mechanism for ensuring registration accuracy. By utilizing the inherent geometric properties of the terrain to constrain the registration process, erroneous matching and accumulated errors are avoided. The constraint construction process first calculates the 3D coordinates of each anchor point using the principle of stereo vision, and reconstructs the spatial position of the anchor points using camera calibration parameters and the spatial distance between the corresponding imaging positions in images from different viewpoints. Then, geometric relationships such as distance constraints, angle constraints, and coplanar constraints between anchor points are established to form a set of constraint equations. Finally, based on the degree of constraint violation and the quality of anchor point features, a reliability weight is assigned to each constraint relationship, with the weight range set to (0,1) to reflect the importance of the constraint relationship in the registration optimization. Among them, the distance constraint requires that the same terrain anchor points maintain consistent Euclidean distance in images from different viewpoints, the angle constraint ensures the geometric invariance of the angle relationship formed by the anchor points in each viewpoint, and the coplanar constraint further enhances the registration stability by utilizing the local planar characteristics of the terrain surface.

[0028] To address the imaging occlusion problem from different viewpoints, an occlusion compensation coefficient is calculated, and the anchor point matching weight is dynamically adjusted based on this coefficient. Occlusion compensation is a key technology for solving multi-view registration in complex terrain. By quantifying and compensating for visibility differences between different viewpoints, the robustness and accuracy of registration are improved. The calculation process first analyzes the imaging geometry of each viewpoint, assessing the visibility of anchor points based on camera pose, terrain slope, and local occlusion. Then, it considers the impact of atmospheric scattering, shadow occlusion, and vegetation cover on feature extraction quality, quantifying the differences in imaging conditions between viewpoints. Finally, considering both visibility and imaging quality factors, the occlusion compensation coefficient for each anchor point in different viewpoint combinations is calculated. The formula for calculating the occlusion compensation coefficient is as follows: In the formula, For terrain anchor points From the perspective The occlusion compensation coefficient is below. Due to the difference in viewing angle, The degree of shadow occlusion, As the degree of vegetation shading, , , These are the weight parameters.

[0029] Based on anchor point features, geometric constraints, and occlusion compensation coefficients, multi-view registration is achieved using a graph optimization algorithm, and the spatial registration accuracy evaluation results are output. Multi-view registration is a process of comprehensively utilizing all constraint information to solve for the optimal transformation parameters. A global optimization objective function is constructed, simultaneously satisfying multiple constraints such as feature matching, geometric consistency, and occlusion compensation. The registration process first transforms the anchor point matching problem into a graph optimization problem, constructing a registration graph model with anchor points as nodes and matching relationships as edges. Then, a composite objective function including feature similarity terms, geometric constraint terms, and occlusion compensation terms is established, and the optimal transformation matrix is ​​solved using a gradient descent algorithm. Finally, the consistency of the registration results is verified, reprojection error and geometric constraint violation are calculated, and the registration accuracy and reliability are evaluated. The objective function is in the form of: In the formula, For the overall registration error, For feature similarity error term, For geometric constraint error terms, To compensate for the occlusion error term, , , These correspond to the weighting coefficients. Registration accuracy typically requires a reprojection error of less than 1 pixel and a geometric constraint violation of less than 0.05 to ensure spatial consistency and mapping accuracy requirements for multi-view image data.

[0030] It should be noted that, in one embodiment of the present invention, before constructing the terrain feature map, it is necessary to use remote sensing image classification technology to perform preliminary identification of the terrain type of the target surveying area from multi-view image data, thereby obtaining the vegetation-covered area, shadow-occluded area, steep slope area and bare surface area within the surveying area; wherein, the preliminary identification process of the corresponding terrain type is prior art, and the steps in this application are described in detail.

[0031] In embodiments of the present invention, the detailed implementation steps for penetrating and separating vegetation-covered areas include: Multispectral vegetation identification and analysis was performed on multi-view image data. The composite features of the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) were calculated to generate preliminary vegetation cover classification results. Morphological closing operations were then used to optimize the classification boundary, resulting in a spatially continuous vegetation cover classification map. Multispectral vegetation identification and analysis is a fundamental technology for vegetation information extraction. By utilizing the unique spectral response differences of vegetation in the near-infrared and visible light bands, accurate identification and classification of vegetation areas can be achieved. The identification process first involves radiometric and geometric correction of the multi-view images to eliminate sensor differences and atmospheric effects, obtaining standardized reflectance data. Then, the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) were calculated separately. A composite vegetation index was constructed based on a weighted combination of NDVI and EVI. Statistical analysis was used to determine the classification threshold, dividing pixels into four coverage levels: no vegetation, sparse vegetation, moderate vegetation, and dense vegetation. The preliminary classification results are processed by morphological closing operations. A 3×3 circular structuring element is used to first expand and then erode to fill small cavities in the vegetation area, eliminate noise interference, and generate a vegetation coverage classification map with good spatial continuity.

[0032] A quantitative analysis of vegetation canopy structure is conducted based on vegetation cover classification maps. This analysis integrates texture roughness, shadow distribution characteristics, and spectral mixture decomposition results to calculate pixel-level vegetation canopy density values, generating a continuously distributed vegetation density map. This quantitative analysis of vegetation canopy structure represents a crucial shift from qualitative classification to quantitative assessment, achieving precise quantification of vegetation density through multi-feature fusion. The quantification process first calculates the texture roughness of the vegetation area, using the gray-level co-occurrence matrix method to extract texture parameters such as contrast, correlation, and entropy. High texture roughness indicates a complex vegetation canopy structure and higher density. Next, shadow proportion characteristics are statistically analyzed. Shadow pixels are identified through the luminance component of the HSV color space, and the ratio of shadow area to total area is calculated. A high shadow proportion reflects strong shading and high density. Simultaneously, linear spectral mixture decomposition technology is employed to decompose mixed pixels into the contribution proportions of vegetation endmembers, soil endmembers, and shadow endmembers. The vegetation endmember proportion directly quantifies vegetation cover density, generating a continuous vegetation density map with density values ​​ranging from 0 to 1, providing a precise quantitative basis for subsequent penetration effect assessment.

[0033] Based on vegetation density maps and optical transmission physics models, the probability distribution of surface optical visibility under different vegetation cover densities was calculated, and a quantitative relationship model between visibility probability and vegetation penetration depth was established. Calculating the optical visibility probability is a core step in evaluating vegetation penetration, quantifying the impact of vegetation on surface observation through optical physics principles. The calculation process is based on Beer-Lambert's law of light attenuation, combined with the influence of factors such as vegetation density, leaf area index, leaf tilt angle distribution, and observation geometry on light propagation. First, a three-dimensional optical model of the vegetation canopy was established, dividing the canopy into multiple thin layers, each with specific extinction coefficients and scattering characteristics. Then, the propagation path of light from the top of the canopy to the ground surface was calculated, considering multiple scattering and absorption effects. Finally, the optical visibility probability of the ground surface was calculated by comprehensively considering propagation attenuation. The formula for calculating the visibility probability is: ; In the formula, This represents the probability of visibility on the ground. The extinction coefficient is the vegetation extinction coefficient, and LAI is the leaf area index per unit area. The density of the vegetation canopy; among which, The calculation process for LAI is existing technology and will not be elaborated upon in this application.

[0034] Based on the obtained visibility probabilities, a visibility probability distribution map of the entire region is generated. The higher the probability value, the easier the surface is to be observed, providing a scientific basis for the zoning processing strategy.

[0035] Based on the statistical distribution characteristics of surface visibility probability, the vegetation cover area is adaptively divided into high visibility area and low visibility area by setting the optimal visibility threshold, and corresponding differentiated elevation acquisition and processing strategies are formulated. Adaptive regional division is the basic step to achieve accurate penetration processing. Through reasonable threshold setting and zoning strategy, processing accuracy and computational efficiency are balanced. The partitioning process first performs statistical analysis on visibility probability, calculating the mean, variance, and quantile characteristics of the probability distribution. The optimal partitioning threshold is determined using the Otsu automatic thresholding algorithm. Areas above the optimal threshold are classified as high-visibility areas, indicating sparse vegetation and relatively clear surface features; areas below the threshold are classified as low-visibility areas, indicating dense vegetation and difficult surface observation. Simultaneously, based on spatial connectivity constraints, connected component analysis is used to eliminate isolated areas with excessively small areas, ensuring the spatial rationality of the partitioning results. For high-visibility areas, direct stereo matching is used to obtain elevation observations, offering high accuracy and fast computational efficiency. For low-visibility areas, a combination of neighborhood topographic trend extrapolation and vegetation canopy height compensation is employed. Based on the topographic gradient characteristics of the surrounding visible areas, the elevation of the obscured surface is estimated by interpolation. Simultaneously, elevation offset compensation is performed based on vegetation density and canopy height information, ultimately generating an estimated surface elevation after vegetation penetration. Each pixel is then labeled with a comprehensive confidence score based on visibility probability, neighborhood consistency, and interpolation error.

[0036] In embodiments of the present invention, the process of performing depth inference on the shadow-occluded region includes: Geometric projection reconstruction of shadow-occluded areas is performed based on oblique photogrammetry images. The geometric features of the shadow boundary and solar radiation angle information are extracted to infer the 3D spatial coordinates of the shadow source. Conditional generative adversarial networks (GANs) are then used to infer the hidden terrain structure of the occluded area. Geometric projection reconstruction is the core technology for solving the shadow occlusion problem. By analyzing the geometric projection relationship of the shadow, the spatial relationship between the occluding object and the occluded terrain is inversely deduced. The reconstruction process first segments the shadow area of ​​the oblique photogrammetry image using thresholding based on the HSV color space and morphological processing to accurately extract the shadow boundary contour. Then, the shooting time is extracted using the EXIF ​​information of the image (EXIF information records key information such as shooting time, equipment parameters, and GPS positioning), and the corresponding solar altitude and azimuth angle parameters are calculated. Based on the geometric features of the shadow contour and the solar angle information, a ray tracing algorithm is used to reconstruct the light propagation path from the solar source to the shadow boundary. The 3D coordinates of the shadow source are determined by calculating the intersection points of the light rays and the terrain surface. The formula for calculating the shadow source location is: ; In the formula, The three-dimensional coordinates of the projection source Let t be the coordinates of the shadow boundary point, and t be the distance the light ray travels. The solar altitude angle, The azimuth of the sun; The system represents the unit direction vector of sunlight. Then, it analyzes the gradient continuity and surface texture direction consistency characteristics of the surrounding terrain in the shaded area to establish terrain continuity constraints. A conditional generative adversarial network architecture is adopted, using the visible terrain in the neighborhood as conditional input. The generator network generates hidden terrain features of the shaded area through an encoder-decoder structure, while the discriminator network ensures the authenticity of the generated results and the physical rationality of the terrain. Network training employs joint optimization of adversarial loss and terrain continuity loss, ultimately outputting the hidden terrain features of the shaded area, providing crucial information for complete terrain reconstruction.

[0037] In an embodiment of the present invention, the process of gradient correction for steep slope regions includes: Perspective distortion correction for steep slope areas is performed based on oblique photogrammetry images. The geometric angle between the local normal vector of the terrain surface and the camera's optical axis is calculated to identify severely distorted areas. The true geometric gradient of the steep slope area is then restored through difference analysis of multi-view observation data. Perspective distortion correction obtains accurate terrain gradient parameters by compensating for geometric deformation in oblique photogrammetry. The correction process first calculates the local normal vector of each point on the terrain surface based on a digital elevation model, then calculates the elevation gradient using the central difference method, and finally obtains the normal vector through the cross product of the gradient vectors. Next, the angle between the normal vector and the camera's imaging optical axis is calculated, and terrain gradient distortion areas where the angle exceeds a preset threshold are identified. Within the distortion areas, the differences in gradient observation values ​​of the same steep slope in images from different viewpoints are analyzed to establish a multi-view gradient fusion model.

[0038] The formula for correcting the actual slope is as follows: ; In the formula, To correct the actual slope, To observe the slope, The angle between the normal vector and the optical axis. The weights are used to determine the reliability of the viewpoints. By weighted averaging of the correction results from multiple viewpoints, the observation bias of a single viewpoint is eliminated, and the gradient correction results for steep slope areas are obtained, providing reliable geometric constraints for accurate terrain modeling.

[0039] In an embodiment of the present invention, the process of enhancing the texture of the exposed surface layer includes: Multi-frequency texture enhancement is performed on exposed land areas based on vertical aerial images. Wavelet transform is used to separate multi-scale frequency components of the surface texture, extracting high-frequency detail features and low-frequency background features. Adaptive contrast adjustment is then performed based on the spatial consistency of texture direction. Multi-frequency texture enhancement is an effective means to improve the quality of surface feature extraction, highlighting surface morphological details at different scales through frequency division. The enhancement process uses Daubechies wavelet for three-level decomposition, separating low-frequency background components and high-frequency detail components. The low-frequency background component reflects the overall grayscale distribution and large-scale terrain changes of the surface, while the high-frequency detail component contains detailed information such as surface texture, cracks, and micro-topography. Then, the local consistency index of the texture direction gradient is calculated. The gradient direction is extracted using the Sobel operator, and the circumferential variance of the gradient direction is statistically analyzed within a sliding window to quantify the degree of texture direction consistency. Based on edge response intensity and direction consistency, adaptive gain adjustment coefficients are calculated. The formula for obtaining the adaptive gain adjustment coefficient is as follows: ; In the formula, For adaptive gain coefficients, These are gain control parameters. The directional consistency coefficient. For edge response strength, The maximum edge intensity is obtained by adjusting the gain of the high-frequency components and reconstructing the image through wavelet inverse transform, resulting in a texture-enhanced surface feature map that effectively improves the recognizability and contrast of surface details.

[0040] In an embodiment of the present invention, the process of constructing a terrain feature map includes: This paper describes a comprehensive topographic feature map that integrates vegetation penetration elevation estimation, terrain features hidden in shaded areas, gradient correction results for steep slopes, and enhanced surface texture maps through multi-layer feature alignment and fusion. Multi-layer feature alignment and fusion is a key technology for integrating heterogeneous topographic information, ensuring consistency of features from different sources in both spatial and semantic dimensions. The alignment process first involves spatial registration, correcting spatial biases from different data sources through feature point matching and affine transformation. Then, resolution unification is performed by resampling all feature layers onto a unified spatial grid using bicubic interpolation. Next, feature standardization is applied to eliminate differences in the dimensions and dynamic range of different feature types. Feature fusion employs a weighted stacking strategy, with weights dynamically allocated based on the quality assessment and spatial coverage of each feature layer. The final topographic feature map contains multi-dimensional feature information such as elevation, slope, texture, and confidence level, providing a comprehensive feature data foundation for subsequent topographic complexity analysis and 3D modeling.

[0041] In an embodiment of the present invention, the process of obtaining the complexity partitioning result includes: A sliding window local complexity quantification analysis is performed on the topographic feature map, comprehensively calculating the variability of elevation gradient, the complexity of curvature distribution, and the heterogeneity of topographic orientation to generate a comprehensive complexity index reflecting the complexity of local topographic changes. The local complexity quantification analysis comprehensively assesses the spatial complexity of topographic changes through multi-dimensional feature fusion. The quantification process adopts a multi-scale sliding window strategy, with window sizes set to three levels: 3×3, 5×5, and 7×7, corresponding to fine, medium, and coarse-scale topographic feature analysis, respectively. Within each window, the standard deviation of the elevation gradient is first calculated, and the gradient components in the horizontal and vertical directions are extracted using the Sobel operator. Then, the standard deviation of the gradient magnitude is calculated to reflect the degree of variation in topographic undulation within the window.

[0042] Next, the entropy value of curvature change is calculated. The principal curvature and Gaussian curvature are calculated using a second-order difference operator. The frequency distribution of curvature values ​​is statistically analyzed, and the uncertainty of the curvature distribution is quantified using the Shannon entropy formula. Simultaneously, the directional heterogeneity of terrain features is assessed. The sliding window is divided into eight 45° sector sub-regions, and the average gradient direction of each sub-region is calculated. The dispersion of the direction is quantified using the circumferential variance. The formula for calculating the terrain local complexity index is: In the formula, This represents the local complexity index of the terrain. For the gradient standard deviation, The entropy value is the curvature change. The index of directional heterogeneity. , , The weighting coefficient is used to generate the local complexity index distribution of the entire region through this formula, providing a quantitative basis for subsequent spatial clustering analysis.

[0043] Spatial autocorrelation pattern recognition based on the local topographic complexity index calculates the spatial correlation strength and anomalous distribution characteristics of complexity in adjacent regions, identifies high-value and low-value clustered regions of complexity, and determines the spatial distribution pattern and preliminary zoning boundaries of topographic complexity. Spatial autocorrelation pattern recognition quantifies spatial proximity effects and clustering characteristics using statistical methods. The recognition process first constructs a spatial weight matrix, defines adjacency relationships using the Queen adjacency criterion, and applies inverse proportional weighting to distance. Then, the global Moran index is calculated to assess the spatial autocorrelation strength of the overall complexity distribution. Next, the local Moran index is calculated to identify the local spatial correlation pattern at each location. Through Local Indicative Spatial Association (LISA) analysis, four spatial patterns are identified: high-high clustering (high complexity surrounded by high complexity), low-low clustering (low complexity surrounded by low complexity), high-low anomaly, and low-high anomaly. Combining spatial connectivity analysis of clustered regions and area threshold screening, statistically significant complexity clustered regions are identified, forming the preliminary zoning boundaries of topographic complexity and providing a spatial framework for subsequent boundary optimization.

[0044] The initial zoning boundaries undergo terrain structure feature refinement and optimization. Terrain fault lines and gradient abrupt changes in the boundary neighborhood are extracted, and the watershed algorithm is used for boundary morphological optimization to generate refined zoning boundaries that conform to the natural terrain structure. Boundary refinement and optimization is a crucial step in ensuring the rationality of the zoning boundaries. Through terrain structure analysis and morphological processing, the zoning boundaries are made consistent with the actual terrain features. The optimization process first extracts terrain fault lines within the buffer zone of the initial boundary, using edge detection operators to identify linear features with abrupt elevation changes, such as ridgelines, valley lines, and slope transition lines. Then, gradient abrupt change zones are detected, and areas with significant slope change rates are identified by calculating the second derivative of the gradient. These areas are typically natural terrain boundaries.

[0045] Next, the watershed algorithm is used to perform morphological optimization on the boundary, taking the local complexity index as the terrain height, and identifying the natural watershed boundary by simulating the water flow convergence process. The watershed algorithm finds local minimum points as water catchment points through gradient descent, and then traces the watershed boundary backward from the water catchment point to form a continuous watershed line. Finally, the terrain fault line, gradient abrupt change zone and watershed boundary are weighted and fused to generate a fine partition boundary that highly matches the natural structure of the terrain, ensuring the terrain rationality and spatial continuity of the partitioning results.

[0046] Based on fine-grained zoning boundaries, the target surveying area is hierarchically divided into complexity zones. According to the statistical distribution characteristics of the complexity index, different levels of terrain complexity regions are defined, and the surveying priority weights for each zone are calculated in conjunction with urban and rural planning needs. Hierarchical complexity zoning is the core step in realizing differentiated surveying strategies. Through quantitative grading and priority allocation, it provides zoning guidance for subsequent surveying processing. The zoning process first performs statistical analysis on the complexity index of the entire region, calculating the mean, standard deviation, and quantile characteristics. The quartile method is used to determine the grading thresholds; for example, below the 25th quantile is a flat baseline area, 25%–50% is a regular undulating area, 50%–75% is a complex steep slope area, and above the 75th quantile is an extreme terrain area. Then, combined with the fine-grained zoning boundaries, the continuous complexity distribution is discretized into four levels of regions, each with relatively uniform terrain complexity characteristics. Next, the surveying priority weights for each zone are calculated, comprehensively considering the zone's complexity index, area proportion, and spatial proximity to key urban and rural planning areas. The formula for calculating the surveying priority weights is: ; In the formula, For surveying priority weights, The normalized complexity index, For area percentage, Distance to key areas of urban and rural planning , , For weight parameters, This is the distance attenuation parameter. The priority weight of each partition is calculated using this formula; a higher weight indicates a higher mapping priority and requires a more refined processing strategy. The final result is a hierarchical complexity partitioning system comprising flat baseline areas, regular undulating areas, complex steep slope areas, and extreme terrain areas. Each partition has clear complexity characteristics and mapping priority, providing a scientific basis for subsequent adaptive mapping processing.

[0047] In embodiments of the present invention, the detailed implementation steps of local complexity metric include: Multi-scale gradient feature extraction and analysis is performed on topographic feature maps. The first-order partial differential gradient and second-order partial differential gradient of topographic elevation are calculated in sliding windows at different spatial scales. The extreme value statistical features of gradient amplitude at each scale level are extracted to generate multi-scale gradient feature vectors that characterize the degree of drastic topographic changes. Multi-scale gradient feature extraction and analysis is a fundamental technique for quantifying topographic complexity. By calculating gradients at different scales, the changing characteristics of topography can be comprehensively depicted. The analysis process employs sliding windows of three scales: 3×3, 5×5, and 7×7, corresponding to the local, neighborhood, and regional spatial levels, respectively. Within each scale window, the first-order gradient of elevation is first calculated using the Sobel and Scharr operators to obtain the gradient components in the x and y directions, followed by the calculation of gradient magnitude and direction angle. Next, the second-order gradient is calculated using the Laplacian operator to reflect the intensity of topographic curvature changes. Statistical analysis is performed on the gradient magnitudes at each scale, extracting characteristic parameters such as maximum, minimum, average, standard deviation, and skewness to construct a multidimensional gradient feature vector. This vector norm comprehensively reflects the overall intensity of topographic gradient changes at different scales, providing a gradient dimension feature description for complex quantification.

[0048] Curvature tensor decomposition is performed on topographic feature maps. The principal curvature and Gaussian curvature distributions of the topographic surface are obtained through Hessian matrix eigenvalue decomposition. The probability distribution characteristics of curvature values ​​within a sliding window are statistically analyzed, and the randomness and complexity of the curvature distribution are calculated using information entropy theory. Curvature tensor decomposition is an important method for describing the geometric complexity of topographic surfaces, revealing the concavity and convexity patterns of the terrain through curvature analysis. The calculation process first constructs the Hessian matrix of the topographic elevation function, containing second-order partial derivatives fxx, fyy, and fxy. Then, two principal curvatures, κ1 and κ2, are obtained through eigenvalue decomposition, representing the degree of curvature of the topographic surface in two principal directions, respectively. Next, the Gaussian curvature K = κ1 × κ2 and the mean curvature H = (κ1 + κ2) / 2 are calculated. The Gaussian curvature reflects the inherent curvature properties of the surface, while the mean curvature reflects the overall curvature of the surface. The frequency distributions of the principal curvature, Gaussian curvature, and mean curvature are statistically analyzed within a sliding window, and a probability density function is constructed using an equally spaced binning method.

[0049] A spatial analysis of directional heterogeneity is performed on topographic feature maps. The sliding window is divided into eight fan-shaped sub-regions based on azimuth angle. The dominant direction of the topographic gradient within each sub-region is calculated, and the consistency of topographic relief direction is quantified through circular statistical analysis of the direction vectors. Spatial analysis of directional heterogeneity is a specialized technique for assessing the complexity of topographic change directions, revealing the spatial organization pattern of topographic relief through directional statistics. The analysis process first divides the sliding window into equal-angle fan-shaped regions with the center point as the origin, according to eight directions: 0°, 45°, 90°, 135°, 180°, 225°, 270°, and 315°. Then, the average directional angle of the elevation gradient is calculated within each fan-shaped region, and the circular property of the angle data is handled using the vector averaging method. Finally, the eight directional angles form a set of direction vectors, and the circular variance of the direction vectors is calculated.

[0050] Data standardization was performed on the multi-scale gradient eigenvector norm, curvature change entropy, and directional heterogeneity index. Principal component analysis (PCA) was used to reduce dimensionality and extract comprehensive features, constructing a unified terrain local complexity index and determining complexity grading thresholds. Comprehensive feature extraction is a key step in fusing multi-dimensional complexity features, and a unified complexity evaluation index was constructed using statistical methods. The process first standardizes the three types of features using Z-scores to eliminate the influence of dimensional differences and numerical ranges among different features. Then, a feature covariance matrix is ​​constructed, and principal component analysis is used to extract the main directions of variation. Next, the contribution rate of each principal component is calculated, and principal components with a cumulative contribution rate of over 85% are selected as comprehensive features. Typically, the contribution rate of the first principal component is between 60% and 70%, which can represent most of the terrain complexity information; therefore, the first principal component is selected as the terrain local complexity index. Finally, statistical analysis was performed on the complexity index, and the mean and standard deviation were calculated. The mean ± 0.5 standard deviation, mean ± standard deviation, and mean ± 1.5 standard deviation were used as the classification thresholds to divide the complexity into five levels: very low, low, medium, high, and very high, providing a quantitative complexity evaluation standard for subsequent partitioning.

[0051] In an embodiment of the present invention, the reconstruction process of the three-dimensional terrain model includes: A differentiated point cloud sampling strategy is implemented based on the complexity partitioning results. The corresponding sampling density level is determined according to the terrain complexity of each partition. Spatial coordinate points are extracted from multi-view image data using stereo vision matching technology to generate initial terrain point cloud data covering the entire region. The differentiated point cloud sampling strategy balances modeling accuracy and computational efficiency through adaptive density control. The sampling process employs a four-level sampling scheme based on the terrain complexity zoning results: A sparse sampling density (e.g., 10-20 points / m²) is used in flat baseline areas to accurately describe gentle terrain with fewer sampling points; a standard sampling density (e.g., 30-50 points / m²) is used in areas with moderate undulations to ensure sufficient representation of the geometric features of moderately undulating terrain; a denser sampling density (e.g., 60-100 points / m²) is used in complex steep slope areas to capture detailed changes in steep terrain; and an ultra-dense sampling density (120-200 points / m²) is used in extreme terrain areas to process complex and rugged terrain through ultra-dense sampling. The stereo matching process uses a semi-global matching algorithm, combining the geometric constraints of multi-view images to reconstruct 3D coordinates through feature point matching and disparity calculation. Simultaneously, the redundancy of multi-view information is utilized to improve the accuracy and completeness of the point cloud. Finally, an initial terrain point cloud containing XYZ coordinates, reflection intensity, and matching confidence is generated, providing the raw data foundation for subsequent quality inspection and model reconstruction.

[0052] Quality defect detection and identification are performed on the initial terrain point cloud. Through point cloud density statistical analysis and neighborhood distance variability assessment, sparse missing regions and anomalous noise points in the point cloud data are identified, and the void types are classified according to the geometric morphological characteristics of the missing regions. Quality defect detection and identification is a crucial step in ensuring the integrity of the point cloud. Data defects are discovered and classified using statistical methods and geometric analysis. The detection process first calculates the local density distribution of each point, counts the number of points within a neighborhood of radius r, and generates a density distribution map. Then, the variance of the K nearest neighbor distance is calculated, and the distance distribution characteristics from each point to its K nearest neighbors are analyzed to identify outlier noise points with significantly large distance variances. Low-density void regions are identified using density threshold discrimination; when the local density is below a preset density threshold, it is marked as a sparse region. High-variance noise points are identified through statistical anomaly detection; when the variance of the K nearest neighbor distance exceeds N times the standard deviation, it is marked as a noise point, where N is an integer; typically N=3.

[0053] For the identified void regions, morphological classification is performed based on the geometric characteristics of the void boundaries: the void area and perimeter are calculated, and the roundness index is calculated as follows: Roundness Index = Void Area ÷ (Void Perimeter) 2 The regularity of the shape is evaluated by ÷4π. Holes with a roundness index greater than the preset index threshold are classified as regular holes, while those with a roundness index not greater than the preset index threshold are classified as irregular holes. This classification provides a basis for selecting appropriate repair algorithms and ensures that different types of defects are treated in a targeted manner.

[0054] A classification-based repair strategy is employed for identified point cloud defects. For holes with regular geometric shapes, radial basis function interpolation is used for reconstruction. For complex, irregular holes, a curvature propagation iterative algorithm is used for filling, ensuring the geometric continuity and physical rationality of the repair results with the surrounding terrain. The classification-based repair strategy is a key method for improving repair quality, achieving optimal defect repair results through algorithm matching. For regular hole repair, multivariate quadratic radial basis function (MQ-RBF) is used for surface interpolation to generate a smooth, completed surface. For irregular hole repair, an iterative algorithm based on neighborhood terrain curvature propagation is used. This algorithm first extracts the curvature tensor information of the hole boundary points and calculates the principal curvature direction and curvature value. Then, a regular mesh is established inside the hole, and initial curvature values ​​are assigned to the mesh nodes using an inverse distance weighting method. Next, iterative propagation updates are performed, with each iteration calculating the Laplace diffusion of curvature and propagating anisotropically along the principal curvature direction, while applying boundary constraints to ensure the continuity of the curvature field. Through multiple iterations until convergence, a hole-filling result consistent with the surrounding terrain geometry is generated.

[0055] The repaired terrain point cloud is modeled using triangulation to construct the topological connectivity of the terrain surface. Through mesh quality optimization and topological anomaly detection, a 3D terrain model with good geometric properties is generated. Triangulation modeling is a crucial step in converting discrete point clouds into continuous curved surfaces, generating a high-quality 3D model through topological reconstruction and geometric optimization. The modeling process first uses the Delaunay triangulation algorithm to construct an initial triangular mesh, ensuring the optimality and circularity of the triangles. Then, mesh quality is evaluated by calculating the aspect ratio and interior angle distribution of the triangles, identifying thin, elongated triangles with poor quality. Next, curvature optimization is performed, adjusting vertex positions using the Laplacian smoothing operator. The optimization formula is: In the formula, This represents the positions of the n neighboring vertices; This represents the original position of the m-th vertex; The optimized vertex position, For smoothing parameters, As vertices The neighborhood, The weighting coefficients are used. Simultaneously, mesh topological anomalies are detected, identifying non-manifold edges (edges connecting more than two faces) and self-intersecting surfaces. A local mesh reconstruction strategy is employed to eliminate these anomalies: non-manifold edges are split, and self-intersecting surfaces are re-triangulated. The final result is a topologically correct and geometrically smooth 3D terrain model, containing vertex coordinates, facet connectivity, and normal vector information, providing a complete 3D terrain representation for subsequent accuracy evaluation and application analysis.

[0056] In an embodiment of the present invention, the detailed implementation steps of the iterative repair method for irregular voids based on neighborhood terrain curvature propagation include: This study analyzes the boundary geometric features of irregular cavities, identifying the 3D coordinate point sequence of the cavity boundary. A continuous geometric contour of the boundary is reconstructed using spline curve fitting technology. Tangent and normal vectors of the boundary curve at each point are extracted, and the curvature tensor and principal curvature directions of the boundary points are calculated. The boundary geometric feature analysis provides reliable constraints for subsequent propagation through a precise boundary description. The analysis process first sorts the point sequence of the cavity boundary, using nearest neighbor search and angle discrimination to ensure the continuity of the boundary points. Then, cubic B-spline curves are used to fit the boundary points, generating a smooth and continuous boundary contour. Next, the first derivative (tangent vector) and second derivative of the fitted curve at each boundary point are calculated, and the curve curvature and normal vector are calculated using the Flyner formula. For the 3D boundary curve, the curvature tensor matrix of each point is constructed, and the principal curvatures κ1 and κ2, as well as the corresponding principal direction vectors e1 and e2, are obtained through eigenvalue decomposition.

[0057] The initial conditions for curvature propagation inside the cavity are constructed based on the boundary curvature tensor information. A uniformly distributed computational grid is established within the cavity region, and initial curvature values ​​are assigned to the grid nodes using a distance-weighted interpolation method, forming the initial state field for curvature propagation. The construction of the initial conditions for curvature propagation is a key step in ensuring the rationality of propagation. Scientific initial value setting ensures the convergence of the iterative process and the physical meaning of the results. The construction process first generates a regular rectangular grid inside the cavity, with the grid spacing adaptively determined according to the cavity size. Then, the Euclidean distance from each grid node to all boundary points is calculated, and a distance weight matrix is ​​constructed. Finally, the inverse distance-weighted (IDW) interpolation method is used to calculate the initial curvature values ​​of the grid nodes.

[0058] Iterative diffusion update calculations are performed on the initial field of curvature propagation based on the initial curvature value. In each iteration, the partial differential diffusion equation of curvature is solved, and anisotropic propagation occurs along the principal curvature directions. Simultaneously, fixed boundary constraints are applied to ensure a smooth transition between the propagation process and the boundary curvature. Iterative diffusion update calculation is the core algorithm for curvature propagation, achieving spatial diffusion of curvature information through numerical solution of the partial differential equation. The calculation process uses the finite difference method to discretize the curvature diffusion equation, considering anisotropic diffusion characteristics. The principal direction of the diffusion tensor is consistent with the local principal curvature direction, and the diffusion coefficient is adaptively adjusted according to the gradient of the curvature value.

[0059] In each iteration, the curvature Laplacian operator for the grid nodes is calculated, and a central difference scheme is used for numerical approximation. Then, weighted diffusion is performed along the principal curvature direction, with the diffusion weight inversely proportional to the curvature gradient to ensure a stronger diffusion effect in regions with gentle curvature changes. Boundary constraint terms are introduced, and the Lagrange multiplier method is used to ensure that the curvature values ​​of the boundary nodes remain fixed, preventing boundary information from drifting during iteration. The convergence criterion for iterative updates is based on the maximum change in the curvature field between adjacent iterations; convergence is considered achieved when the change is less than a preset variable threshold.

[0060] Curvature-driven surface reconstruction is performed based on a convergent curvature field. Curvature information is used to constrain the surface geometry, generating a completed surface for void regions. The reconstruction quality is evaluated through surface smoothness checks. Curvature-driven surface reconstruction is the inverse process of recovering a geometric surface from a curvature field. A curvature-constrained surface fitting algorithm generates a completed surface that harmonizes with the surrounding terrain. The reconstruction process employs an energy minimization method based on curvature constraints, constructing a composite energy function that includes surface smoothness and curvature fitting terms. An optimal surface is solved using a variational method, ensuring that the curvature distribution of the reconstructed surface most closely approximates the target curvature field while maintaining overall surface smoothness. After reconstruction, surface quality is evaluated, and the total curvature variation of the surface is calculated as a smoothness index.

[0061] in, Let S be the variation of total curvature, and S be the surface region. The curvature gradient is used. If the total curvature variation exceeds a preset variation threshold, it indicates that the surface exhibits excessive oscillation or is not smooth. In this case, it is necessary to backtrack and adjust the propagation parameters (such as the diffusion coefficient and boundary weights) and re-execute the iteration process until a high-quality repair result that meets the smoothness requirements is generated. The final output completed surface is continuous with the original terrain at the boundary, ensuring the geometric continuity and visual realism of the overall terrain model.

[0062] In an embodiment of the present invention, the detailed implementation steps for obtaining the final mapping model and its confidence distribution map include: Multi-view backprojection accuracy verification of a 3D terrain model is performed. The reconstructed 3D terrain model is projected onto the image planes of each original viewpoint. Through geometric difference measurement and image similarity analysis, regions and error distribution characteristics that fail to meet the model reconstruction quality standards are identified. Multi-view backprojection accuracy verification comprehensively evaluates the reconstruction accuracy and reliability by comparing the differences between model predictions and actual observations. The verification process first involves backprojecting the 3D terrain model onto the corresponding image plane according to the camera parameters and geometric relationships of each viewpoint, generating the model-predicted contour lines and surface textures. Then, terrain contour features are extracted from the actual images, and edge detection algorithms are used to identify terrain boundaries and salient feature lines. Next, the Hausdorff distance between the projected contour and the actual image contour is calculated to quantify the degree of geometric shape matching. Simultaneously, the structural similarity index (SSIM) between the virtual viewpoint image generated by the model and the real image is calculated to evaluate the matching quality of texture and lighting.

[0063] By comprehensively evaluating geometric differences and image similarity, regions where the Hausdorff distance exceeds a preset distance threshold or the SSIM is less than a similarity threshold are identified as areas of difference between the model and the actual terrain, providing target areas for subsequent accuracy analysis.

[0064] Local accuracy quantitative assessment is performed on the identified discrepancy areas. Through the deployment of virtual control points and statistical analysis of elevation deviations, the mapping accuracy level of the discrepancy areas is quantified, and the causal mechanism of insufficient accuracy is analyzed based on the spatial distribution pattern of errors. Local accuracy quantitative assessment is the core method for accurately quantifying model quality. Statistical analysis identifies the severity and distribution pattern of errors. The assessment process first involves deploying virtual terrain control points in the discrepancy areas according to a regular grid. The density of control points is adaptively determined based on the degree of difference, typically 4–16 points per 100 square meters. Then, the true elevation values ​​of the control points are obtained through high-precision GPS measurements or lidar data, serving as the benchmark for accuracy assessment. Next, the elevation deviation between the model elevation and the true elevation of each control point is calculated, and the root mean square error (RMSE) of the deviation is statistically analyzed.

[0065] The accuracy level of the regions with discrepancies is determined based on the root mean square error: excellent, good, average, poor, and very poor. Spatial autocorrelation analysis and cluster analysis are used to identify the spatial distribution patterns of regions with insufficient accuracy. The main causes of insufficient accuracy are analyzed by considering factors such as terrain complexity, data quality, and algorithm parameters.

[0066] Based on error causal analysis, a spatial propagation tracking model of errors in the processing chain is constructed. By reverse analysis of the cumulative path of errors from data acquisition to model reconstruction at each stage, the error amplification effect of each processing link is quantified, and the key links with the greatest impact on the final accuracy are identified. The error spatial propagation tracking model reveals the generation mechanism and propagation law of errors through systematic analysis. The model construction process first establishes a complete processing flowchart from raw image data to a 3D terrain model, including major links such as image preprocessing, feature extraction, multi-view registration, point cloud generation, defect repair, and mesh reconstruction. Then, the processing history of areas with insufficient accuracy is traced, and intermediate result data and parameter settings of these areas in each processing link are extracted. Next, the cumulative effect of errors in each link is analyzed through error propagation theory, and a set of error propagation equations is established. For each processing link, a known small perturbation is applied to the input data, and the degree of influence of the perturbation on the output result is observed, and the error propagation coefficient is calculated. Then, all possible propagation paths from the initial input to the final output are constructed, the cumulative error amplification factor of each path is calculated, and the key propagation paths and key processing links with amplification factors exceeding the critical threshold are identified, providing accurate target positioning for parameter optimization.

[0067] Parameter backtracking optimization is implemented for key identification steps. By constructing a multi-objective optimization function and intelligent optimization algorithm, the configuration of key parameters is iteratively adjusted until the optimal mapping accuracy is achieved, outputting a fully optimized final mapping model and confidence distribution map. Parameter backtracking optimization is the core technology for improving model accuracy, maximizing global accuracy through systematic parameter optimization. The optimization process first constructs a multi-objective optimization function with mapping accuracy as the primary objective and computational efficiency as the secondary objective, using key parameters such as hierarchical extraction parameters, point cloud sampling density, and mesh reconstruction parameters as optimization variables. Then, a particle swarm optimization algorithm is used for global optimization. The particle swarm optimization algorithm has the advantages of strong global search capability and fast convergence speed, making it suitable for handling multi-parameter nonlinear optimization problems. Through continuous iterative optimization, when the accuracy improvement of several consecutive iterations is less than the preset accuracy threshold, it is considered to have converged or reached the maximum number of iterations, and the optimal parameter configuration is output. Finally, a fully optimized 3D terrain model is generated, and a confidence value is calculated for each grid cell. The confidence value comprehensively considers factors such as measurement accuracy, algorithm reliability, and data quality, and generates an elevation confidence distribution map of the entire region, providing a quality assessment basis for the application of surveying and mapping results.

[0068] In this embodiment of the invention, the detailed implementation steps for backtracking and adjusting parameters in key processes include: For low-precision difference areas identified in the accuracy assessment, the error sources in the processing chain are traced backward. Through processing log analysis and intermediate data extraction, the complete data transformation path of this area in the entire surveying and mapping process is reconstructed, and a data flow map from the original multi-view images to the final 3D model is constructed. The reverse tracing of the error sources in the processing chain identifies the error generation nodes and propagation mechanisms through systematic backtracking analysis. The tracing process first extracts the processing records of the insufficiently accurate areas at each stage of image registration, feature extraction, point cloud generation, and 3D reconstruction, including algorithm parameter settings, intermediate calculation results, and quality assessment indicators. Then, according to the time sequence and logical relationship, the complete processing history of the data in this area is reconstructed, and the input data, algorithm operations, and output results of each processing node are marked. Next, the quality change trend of intermediate results at each stage is analyzed. By comparing the data accuracy and completeness before and after processing, the possible links where errors are introduced are initially identified. A complete data flow graph containing node attributes, connection relationships and data flow is established. Nodes represent processing links, edges represent data transmission relationships, and node weights reflect the degree of influence of the link on the final accuracy. This data flow graph provides a structured analytical framework for subsequent sensitivity analysis, ensuring the comprehensiveness and accuracy of error tracing.

[0069] Error sensitivity quantitative analysis is performed on each processing node in the data flow graph. Through controlled perturbation experiments and sensitivity coefficient calculations, the influence of each node on the output accuracy is quantified, and the sensitive nodes and parameters that contribute the most to the final error are identified. Error sensitivity quantitative analysis is the core method for determining key optimization objectives. By quantitatively assessing the impact of parameter changes on the results, it guides the formulation of optimization strategies. The analysis process adopts a single-factor perturbation experiment design, applying small perturbations to the key parameters of each processing node while keeping other parameters unchanged, and observing the magnitude of change in the final model accuracy. Then, the error propagation coefficient of each node is calculated, defined as the ratio of the change in output error to the amount of input perturbation.

[0070] An error propagation path analysis tree is constructed based on the error propagation coefficients of each node. This tree calculates the cumulative amplification effect of errors on all possible paths from the initial data input to the final model output, identifying key error propagation paths where the cumulative amplification effect exceeds a significance threshold. Error propagation path analysis is a crucial method for identifying error amplification mechanisms in systems, revealing the accumulation patterns of errors in complex processing chains through path analysis. The analysis process first converts the data flow graph into a directed acyclic graph (DAG), assigning a weight to each edge based on its corresponding error propagation coefficient. Then, a depth-first search algorithm is used to traverse all possible paths from the starting node to the terminal node, calculating the cumulative error amplification factor for each path. The cumulative factor equals the product of all propagation coefficients on the path. Next, key propagation paths with cumulative amplification factors exceeding a preset threshold are identified; these paths are the main channels for error amplification. Finally, the node with the largest error propagation coefficient on the key propagation path is located, identifying it as a critical link in error accumulation. These critical links are typically located at the beginning or middle of the processing chain, significantly impacting all subsequent steps and serving as key targets for parameter tuning.

[0071] For key stages of error accumulation in the identification process, parameter backpropagation optimization based on gradient descent is implemented. A parameter optimization framework aimed at maximizing model accuracy is constructed, and a systematic reduction in global error is achieved through iterative adjustment of key stage parameters. Parameter backpropagation optimization utilizes intelligent optimization techniques to achieve automatic parameter tuning and systematic error suppression. The optimization process first constructs an optimization problem with the final model's RMSE as the objective function, using adjustable parameters of key stages as optimization variables and establishing parameter constraints to ensure the physical rationality of the parameters. Then, the gradient descent algorithm is used for optimization, with gradient calculation achieved through automatic differentiation technology to avoid the accuracy loss associated with numerical differentiation. The parameter update formula is: In the formula, and These represent the parameter vectors at training iterations t+1 and t, respectively. For learning rate, Let be the objective function. The algorithm employs a gradient operator; and to avoid local optima, adaptive learning rate adjustment and momentum term correction are used to improve the convergence and stability of the optimization. Global error changes are monitored in real time during the optimization process, and convergence is considered achieved when the error improvement over several consecutive iterations is less than a preset accuracy threshold. The optimized parameter configuration is output, including optimal parameter values ​​for each key stage, an evaluation of error suppression effectiveness, and a parameter sensitivity analysis report, providing parameter setting references and error control experience for subsequent similar projects.

[0072] This invention achieves fully automated processing of high-precision 3D mapping under complex terrain conditions through multi-view heterogeneous data acquisition, hierarchical decoupling feature extraction, complexity adaptive partitioning, point cloud defect repair and reconstruction, and error propagation tracing optimization.

[0073] Example 2 Please see Figure 2 As shown, the parts not described in detail in this embodiment are described in Embodiment 1. An automatic mapping method for complex terrain in urban and rural planning based on image analysis is provided. This method is implemented based on the automatic mapping system for complex terrain in urban and rural planning based on image analysis as described in any one of claims 1 to 9. Its features include: Step S1: Obtain multi-view heterogeneous image data of the target mapping area. Extract the features of the terrain anchor points corresponding to each view image in the multi-view heterogeneous image, and combine the spatial geometric constraint relationship of the terrain anchor points with the view occlusion compensation coefficient to perform multi-view registration to obtain multi-view image data. The multi-view image data includes vertical aerial images, oblique photography images and satellite remote sensing images. Step S2: Based on the spectral texture and geometric 3D features of multi-view image data, the survey area is decoupled in layers. The decoupling includes penetration separation of vegetation-covered areas, depth inference of shadowed areas, gradient correction of steep slope areas, and texture enhancement of bare surface areas to obtain a terrain feature map. Step S3: Obtain the local complexity index of the terrain based on the spatial distribution dispersion and gradient mutation frequency of terrain features at each level within the terrain feature map; and divide the survey area according to the clustering distribution characteristics of the local complexity index to obtain the complexity partitioning results. Step S4: Generate an initial terrain point cloud based on the complex regionalization results, and reconstruct a three-dimensional terrain model by identifying hollow areas and noise outliers in the initial terrain point cloud, combined with the curvature propagation characteristics of the neighboring terrain and topological continuity constraints. Step S5: Based on the 3D terrain model, identify areas with insufficient accuracy in the 3D terrain model by analyzing the back projection consistency of multi-view image data and the elevation deviation of terrain control points. Based on the elevation deviation analysis results, obtain the spatial propagation path and cumulative effect of the deviation, dynamically adjust the layered extraction parameters and reconstruction strategy, and obtain the final mapping model and its confidence distribution map.

[0074] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0075] All formulas in this manual are dimensionless and calculated numerically. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

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

Claims

1. An automatic mapping system for complex terrain in urban and rural planning based on image analysis, characterized in that, include: The data acquisition module is used to acquire multi-view heterogeneous image data of the target mapping area. By extracting the features of the terrain anchor points corresponding to each view in the multi-view heterogeneous images, and combining the spatial geometric constraint relationship of the terrain anchor points with the view occlusion compensation coefficient, multi-view registration is performed to obtain multi-view image data, which includes vertical aerial images, oblique photography images and satellite remote sensing images. The feature extraction module performs layered decoupling on the surveyed area based on the spectral texture and geometric 3D features of multi-view image data. The layered decoupling includes penetration separation of vegetation-covered areas, depth inference of shadowed areas, gradient correction of steep slope areas, and texture enhancement of exposed surface areas to obtain a terrain feature map. The terrain zoning module is used to obtain the local complexity index of terrain based on the spatial distribution dispersion and gradient mutation frequency of terrain features at each level within the terrain feature map; and to divide the survey area based on the clustering distribution characteristics of the local complexity index of terrain to obtain the complexity zoning results. The model building module is used to generate an initial terrain point cloud based on the complex regional results, and to reconstruct a three-dimensional terrain model by identifying hollow areas and noise outliers in the initial terrain point cloud and combining the curvature propagation characteristics and topological continuity constraints of the neighboring terrain. The model optimization module is used to identify areas with insufficient accuracy in the 3D terrain model by analyzing the back projection consistency of multi-view image data and the elevation deviation of terrain control points. Based on the elevation deviation analysis results, it obtains the spatial propagation path and cumulative effect of the deviation, dynamically adjusts the layered extraction parameters and reconstruction strategy, and obtains the final mapping model and its confidence distribution map.

2. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 1, characterized in that, The process of constructing the terrain feature map includes: Based on satellite remote sensing images, vegetation-covered areas are penetrated and separated. The difference between the near-infrared reflectance and the visible light spectrum of vegetation-covered areas in satellite remote sensing images is analyzed. Based on the correlation model between vegetation density and the angle of incidence, the elevation shift is predicted, and the surface elevation after vegetation penetration is estimated. Based on oblique photogrammetry images, depth reasoning is performed on the shadow-occluded area to extract the geometric contour of the shadow boundary and the solar elevation angle information, reconstruct the three-dimensional spatial location of the shadow projection source, and combine the continuity of the neighborhood terrain gradient and the directional consistency of the surface texture of the shadow-occluded area to infer the hidden terrain features of the shadow-occluded area. Gradient correction of steep slope areas is performed based on oblique photography images. By calculating the angle between the local normal vector of the terrain surface and the imaging optical axis, the terrain gradient distortion area is identified, and the true slope of the steep slope area is restored based on the difference in gradient observation values ​​of the terrain gradient distortion area in images from different viewpoints. Texture enhancement of exposed surface layers based on vertical aerial images: By performing multi-scale wavelet decomposition on the vertical aerial images, the high-frequency detail components and low-frequency background components of the surface texture are separated, and adaptive gain adjustment is performed on the high-frequency detail components to obtain a texture-enhanced surface feature map. The surface elevation estimate, hidden terrain features, gradient correction results, and surface feature map are aligned to obtain the terrain feature map.

3. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 1, characterized in that, The process of obtaining the complexity partitioning results includes: The local complexity of the terrain feature map is quantified to generate a terrain local complexity index; Spatial autocorrelation analysis is performed on the local terrain complexity index. By calculating the Moran index and local spatial outliers of the complexity index of adjacent regions, high-value and low-value clustering regions of complexity are identified. Based on the spatial connectivity and area threshold of the clustering regions, preliminary partition boundaries are determined. The preliminary partition boundary is refined by extracting terrain fault lines and gradient abrupt change zones in the boundary neighborhood and performing morphological optimization of the boundary using the watershed algorithm to generate a refined partition boundary. Based on the fine partition boundaries, the target mapping area is divided into boundaries to obtain complexity partitioning results. The complexity partitioning results include flat baseline areas, regular undulating areas, complex steep slope areas, and extreme terrain areas, and each partition is assigned a mapping priority weight.

4. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 1, characterized in that, The reconstruction process of a 3D terrain model includes: The complexity partitioning results are applied to a partitioned adaptive point cloud generation strategy, and an initial terrain point cloud is generated by stereo matching of multi-view image data. Hole detection and noise identification are performed on the initial terrain point cloud. By calculating the local density distribution and K-nearest neighbor distance variance of the initial terrain point cloud, low-density hole regions and high-variance noise outliers in the initial terrain point cloud are identified. Based on the area size and shape complexity of the holes, the holes are classified into regular holes and irregular holes. The regular voids are repaired using a surface interpolation method based on radial basis functions, while the irregular voids are repaired using an iterative repair method based on neighborhood terrain curvature propagation to gradually fill the void region until the surface smoothness constraint is met. The initial terrain point cloud after repair is reconstructed by triangulation and the curvature of the mesh surface is optimized. At the same time, non-manifold edges and self-intersecting surfaces in the mesh are detected, and a local mesh reconstruction strategy is used to eliminate topological anomalies and generate a three-dimensional terrain model.

5. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 1, characterized in that, The process of obtaining the final mapping model and its confidence distribution map includes: The three-dimensional terrain model is back-projected onto the image planes of each viewpoint, the Hausdorff distance between the projected contour and the actual image contour is calculated, and the difference areas between the model and the actual terrain are identified by comparing the structural similarity index between the virtual viewpoint image generated by the model and the real image. Local accuracy assessment is performed on the discrepancy areas by setting up virtual terrain control points within the discrepancy areas and calculating the root mean square error between the model elevation and the actual elevation of the terrain control points. The accuracy deficiency level of the discrepancy areas is determined, and the causes of the error are analyzed based on the spatial distribution pattern of the accuracy deficiency levels. Based on the error cause analysis results, an error spatial propagation model is constructed. By tracing the propagation path of the error from the initial point cloud to the 3D terrain model and calculating the error amplification coefficient of each processing stage along the path, the key stages of error accumulation are identified, and parameters are backtracked and adjusted for the key stages. The adjusted parameters are iteratively optimized by constructing a nonlinear optimization problem with mapping accuracy as the objective function and hierarchical extraction parameters and reconstruction strategy parameters as optimization variables. Particle swarm optimization algorithm is used for global optimization, and the final mapping model and elevation confidence distribution maps of each region are output.

6. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 2, characterized in that, The process of penetrating and separating vegetation-covered areas includes: Spectral analysis was performed on the vegetation cover area in the multi-view image data. By calculating the combined features of the normalized vegetation index and the enhanced vegetation index, a preliminary vegetation cover classification map was generated. Morphological closing operation was performed on the classification map to obtain the vegetation cover classification map. The vegetation coverage classification map is quantified by density analysis. By analyzing the texture roughness and shadow ratio within the vegetation coverage area and combining the spectral mixing decomposition results of the vegetation, a vegetation density map is generated. Based on the vegetation density map, calculate the visibility probability of the ground surface under different vegetation densities; Based on the visibility probability, the vegetation-covered area is divided into high-visibility areas and low-visibility areas; Direct elevation observations are used for high visibility areas, while neighboring topographic trend extrapolation is used for low visibility areas. By weighted fusion of direct observations and extrapolated values, a surface elevation estimate after vegetation penetration is generated, and the elevation confidence level of each pixel is labeled.

7. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 3, characterized in that, The implementation process of local complexity metric includes: Multi-scale gradient calculation is performed on the terrain feature map to generate multi-scale gradient feature vectors; Perform curvature tensor analysis on the terrain feature map to generate curvature change entropy values; The terrain feature map is evaluated for directional heterogeneity by dividing the sliding window into eight fan-shaped sub-regions and calculating the average gradient direction of each sub-region. Based on the circumferential variance of the eight directional vectors, a directional heterogeneity index is generated. The norm, curvature change entropy, and directional heterogeneity index of the multi-scale gradient feature vector are normalized, and the comprehensive features are extracted by principal component analysis to determine the threshold for complexity classification.

8. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 4, characterized in that, The implementation process of the iterative repair method based on neighborhood terrain curvature propagation for the aforementioned irregular voids includes: Boundary features are extracted from the irregular cavity. The three-dimensional coordinate sequence of the cavity boundary point cloud is identified, and curve fitting is performed on the cavity boundary points to calculate the curvature tensor and its principal direction at each boundary point. Based on the curvature tensor of the cavity boundary points, an initial field for curvature propagation is constructed by establishing a regular sampling grid inside the cavity and assigning initial curvature values ​​to the grid nodes. The initial field of curvature propagation is iteratively propagated and updated by calculating the curvature Laplacian operator of the grid nodes in each iteration and performing anisotropic diffusion along the principal curvature direction; Based on the curvature field after iterative convergence, a complete surface for the hole region is generated by a curvature-driven surface reconstruction method. The surface smoothness of the complete surface is evaluated. The total curvature variation of the surface is calculated. If the total curvature variation exceeds the preset variation threshold, the propagation parameters are backtracked and the iteration is repeated until a repair result that meets the smoothness constraint is generated.

9. The automatic mapping system for complex terrain in urban and rural planning based on image analysis according to claim 5, characterized in that, The process of backtracking and adjusting parameters in key processes includes: For areas with insufficient accuracy and low level of difference, the error source is traced. The processing records of the area in the three-dimensional model reconstruction, point cloud generation, feature extraction and image registration stages are traced in reverse. Intermediate result data of the area are extracted at each stage to establish a complete data flow graph from the original image to the three-dimensional terrain model. An error sensitivity analysis is performed on each processing node in the data flow graph, and the error propagation coefficient of each node is calculated. The error propagation coefficient is the ratio of the change in output error to the change in input error. Based on the error propagation coefficient, an error accumulation path tree is constructed. By calculating the cumulative error amplification factor of all possible paths from the initial input to the final output, and identifying the critical propagation path with the cumulative amplification factor exceeding the threshold, the node with the largest error propagation coefficient on the critical propagation path is determined as the critical link of error accumulation. For the critical link of error accumulation, a gradient descent optimization framework is constructed with the final model accuracy as the target and the critical link parameters as the optimization variables. The critical link parameters are iteratively adjusted until the global error converges, and the optimized parameter configuration and its corresponding error suppression effect analysis are output.

10. An automatic mapping method for complex terrain in urban and rural planning based on image analysis, implemented based on the automatic mapping system for complex terrain in urban and rural planning based on image analysis as described in any one of claims 1 to 9, characterized in that, include: Step S1: Obtain multi-view heterogeneous image data of the target mapping area. Extract the features of the terrain anchor points corresponding to each view image in the multi-view heterogeneous image, and combine the spatial geometric constraint relationship of the terrain anchor points with the view occlusion compensation coefficient to perform multi-view registration to obtain multi-view image data. The multi-view image data includes vertical aerial images, oblique photography images and satellite remote sensing images. Step S2: Based on the spectral texture and geometric 3D features of multi-view image data, the survey area is decoupled in layers. The decoupling includes penetration separation of vegetation-covered areas, depth inference of shadowed areas, gradient correction of steep slope areas, and texture enhancement of bare surface areas to obtain a terrain feature map. Step S3: Obtain the local complexity index of the terrain based on the spatial distribution dispersion and gradient mutation frequency of terrain features at each level within the terrain feature map; and divide the survey area according to the clustering distribution characteristics of the local complexity index to obtain the complexity partitioning results. Step S4: Generate an initial terrain point cloud based on the complex regionalization results, and reconstruct a three-dimensional terrain model by identifying hollow areas and noise outliers in the initial terrain point cloud, combined with the curvature propagation characteristics of the neighboring terrain and topological continuity constraints. Step S5: Based on the 3D terrain model, identify areas with insufficient accuracy in the 3D terrain model by analyzing the back projection consistency of multi-view image data and the elevation deviation of terrain control points. Based on the elevation deviation analysis results, obtain the spatial propagation path and cumulative effect of the deviation, dynamically adjust the layered extraction parameters and reconstruction strategy, and obtain the final mapping model and its confidence distribution map.

Images