Lidar-based three-dimensional terrain mapping method

Abstract

The application discloses a three-dimensional terrain mapping method based on a laser radar, and relates to the technical field of three-dimensional terrain mapping. The method is used for improving the terrain mapping precision and spatial consistency. Airborne laser radar point cloud data is acquired, and multispectral images and synthetic aperture radar data are synchronously collected. Matching points are recognized and matched through point cloud density and multi-source feature correlation, a point cloud quality evaluation atlas is generated by establishing multi-modal registration mapping, the fusion weight of multi-source data is adjusted according to the quality atlas and the density distribution, the best fusion ratio is determined through feature consistency analysis, and enhanced point cloud data is generated. The elevation variation, slope and curvature features of the enhanced point cloud are extracted, the terrain elements are divided based on the spatial distribution law, and the multi-modal fusion terrain classification map is generated in combination with the multi-source features. Further, the abnormal area of the terrain element boundary is detected, the correlation between the abnormal area and the fusion parameters is analyzed, the abnormal spatial distribution map is generated, and the closed-loop calculation and processing from multi-source data collection to abnormal analysis is realized.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional terrain mapping technology, specifically a three-dimensional terrain mapping method based on lidar. Background Technology

[0002] With the increasing demand for applications such as urban planning, geological surveying, and environmental monitoring, higher requirements are being placed on high-precision, full-coverage, and dynamically updatable 3D terrain data. Traditional single-sensor 3D mapping methods are prone to point cloud holes and measurement noise in complex terrain or areas with high vegetation cover, making it difficult to guarantee spatial consistency and feature integrity.

[0003] Existing multimodal data fusion methods suffer from substantial shortcomings in image and point cloud processing: First, the spatial registration accuracy between point clouds and imagery or radar data is unstable, leading to significant errors in local geometric and textural features in the fused point cloud; second, the multi-source data fusion process typically relies on fixed weights or local experience, making it difficult to achieve adaptive optimization in regions with varying terrain density and complexity; finally, terrain feature extraction and classification often depend on single geometric or spectral features, failing to fully utilize the correlation between multi-source information, thus reducing the accuracy of terrain unit segmentation and anomaly detection. These deficiencies limit the application effectiveness of 3D terrain mapping in detailed analysis and dynamic monitoring. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a three-dimensional terrain mapping method based on lidar, which solves the problems mentioned in the background.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a three-dimensional terrain mapping method based on lidar, comprising the following steps: S1: acquiring raw point cloud data scanned by airborne lidar, and simultaneously collecting reflectivity data from multispectral images and backscattering coefficients from synthetic aperture radar; identifying feature matching points in the point cloud data based on the correlation between point cloud density distribution and multi-source features, establishing a registration mapping relationship for multi-modal data, and generating a point cloud quality assessment map; S2: adjusting the fusion weight configuration of multi-source data based on the sum and density distribution characteristics of the point cloud quality assessment map, and testing the feature table of different sensor data. S3: Determine the optimal fusion ratio of each data source based on the consistency of the data, and generate enhanced point cloud data with spatial consistency characteristics; S4: Extract terrain features from the enhanced point cloud data, including calculating elevation variation features, slope distribution features, and curvature change features. Divide the terrain into units according to the spatial distribution pattern of the features, and establish a correspondence with spectral and backscattering features to obtain a multimodal fused terrain classification map; S5: Based on the category distribution features in the terrain classification map, detect abnormal change areas at the boundaries of terrain units, analyze the correlation between abnormal areas and multi-source data fusion parameters, and generate an abnormal spatial distribution map.

[0006] Furthermore, the specific process of acquiring raw point cloud data from airborne lidar scanning and simultaneously collecting reflectivity data from multispectral images and backscattering coefficients from synthetic aperture radar is as follows: A multi-source data acquisition parameter mapping table is established to record the spatial attitude parameters and environmental parameters of each sensor during acquisition; radiometric calibration and atmospheric correction are performed on the multispectral image data to convert digital quantization values ​​into surface reflectivity data; incident angle correction and terrain radiometric correction are performed on the synthetic aperture radar data to eliminate the influence of terrain undulations on the backscattering coefficient; and multispectral reflectivity data and radar backscattering coefficients are registered to the spatial location of the lidar point cloud through spatial interpolation.

[0007] Furthermore, based on the correlation between point cloud density distribution and multi-source features, the specific process of identifying feature matching points in point cloud data, establishing registration mapping relationships for multimodal data, and generating point cloud quality assessment maps is as follows: Based on point cloud density distribution characteristics, salient feature regions are identified; texture and spectral features of the corresponding regions are extracted from multispectral images; scattering characteristics are extracted from synthetic aperture radar data; a correspondence between point cloud features and multi-source features is established using a feature similarity measurement algorithm; spatial transformation parameters between multi-source data are calculated; based on the spatial transformation parameters, the local geometric error, registration residual, density consistency, and multi-source feature consistency of each point cloud point are mapped to multi-dimensional quality indicators, and a three-dimensional quality assessment grid is constructed; the assessment grid is divided into regions based on density distribution characteristics, the quality score of each region is calculated, and a point cloud quality assessment map is generated.

[0008] Furthermore, based on the density distribution characteristics of the point cloud quality assessment map, the specific process for adjusting the fusion weight configuration of multi-source data is as follows: Based on the multi-dimensional quality indicators in the point cloud quality assessment map, a quality-weight mapping function is established to analyze the contribution of each quality dimension to the data fusion effect; based on the density distribution characteristics, regions with different terrain complexity are identified, a dynamic weight adjustment strategy is adopted for high-density variation regions, and a fixed weight configuration is adopted for low-density uniform regions; through an adaptive weight optimization algorithm, the optimal fusion weight of each sensor data at different spatial locations is calculated, and the weight ratio of regions with registration errors or inconsistent features is reduced, ultimately generating a multi-source data fusion weight configuration scheme.

[0009] Furthermore, by testing the consistency of feature representations among different sensor data, the optimal fusion ratio of each data source is determined, and the specific process for generating enhanced point cloud data with spatial consistency features is as follows: A multimodal feature consistency evaluation framework is constructed, and data consistency is evaluated by calculating the differences in feature representations of different sensor data at the same spatial location; a multi-objective optimization algorithm is used, with the optimization objectives of maximizing feature consistency and minimizing geometric error, and the optimal fusion ratio of each data source is determined through Pareto optimal solution search; based on the optimal fusion ratio, an enhanced point cloud dataset is generated using a weighted fusion algorithm.

[0010] Furthermore, the specific process for extracting terrain features from the enhanced point cloud data, including calculating elevation variation features, slope distribution features, and curvature change features, is as follows: The enhanced point cloud data is spatially gridded, and a local neighborhood point set is constructed within each grid cell. The elevation variation features of each grid cell are calculated, including the local elevation mean, elevation standard deviation, and elevation range. A local terrain surface is fitted using the least squares method, and the angle between the surface normal vector and the vertical direction is calculated as the slope value. The slope distribution of all points within the grid is then statistically analyzed. Curvature change features, including maximum curvature, minimum curvature, and average curvature, are calculated based on the residuals from the quadratic surface fitting, and all feature values ​​are stored in the point cloud attributes.

[0011] Furthermore, the specific process of dividing terrain units according to the spatial distribution pattern of features and establishing a correspondence between spectral and backscattering features to obtain a multimodal fused terrain classification map is as follows: Multi-dimensional terrain features are spatially clustered using a clustering algorithm, and the point cloud is divided into different terrain units based on feature similarity; the mean of multispectral reflectance and the variance of radar backscattering coefficient are calculated for each terrain unit, and a correspondence table between terrain type and multi-source features is established; the terrain geometric features and multi-source features are integrated using a weighted fusion algorithm, the terrain unit boundaries are optimized, and the categories of each terrain unit are labeled according to the fusion results to generate a multimodal fused terrain classification map.

[0012] Furthermore, based on the category distribution characteristics in the terrain classification map, the specific process for detecting abnormal change areas at the boundaries of terrain units is as follows: extract the boundary lines of each category region in the terrain classification map, enhance the boundary features using the Canny edge detection algorithm, and calculate the geometric feature differences such as elevation difference, slope difference, and curvature difference of the point clouds on both sides of the boundary; identify abnormal change points that exceed the normal range by setting an adaptive threshold, verify the abnormal points by combining multispectral and radar features, cluster continuous abnormal points to form abnormal change areas, and record their spatial location and feature attributes.

[0013] Furthermore, the specific process of analyzing the correlation between abnormal regions and multi-source data fusion parameters to generate anomaly spatial distribution maps is as follows: extract the multi-source data fusion parameters corresponding to the abnormal regions, including the weight configuration and registration accuracy parameters of each sensor; calculate the correlation coefficient between the feature values ​​of the abnormal regions and the fusion parameters, analyze the influence of parameter settings on the anomaly detection results, generate anomaly intensity distribution maps through spatial interpolation, and visualize the spatial distribution patterns of abnormal regions through heatmaps, generating thematic maps containing anomaly types, spatial distribution characteristics, and parameter correlation analysis.

[0014] The present invention has the following beneficial effects:

[0015] (1) A three-dimensional terrain mapping method based on lidar is proposed. By establishing a multi-source data acquisition parameter mapping table and performing radiometric and terrain correction on data from different sensors, the spatiotemporal consistency and quality standardization of the data are achieved. Furthermore, a point cloud quality assessment map is constructed using multi-dimensional quality indicators to provide an accurate reference for data fusion. Based on the quality map and point cloud density distribution, an adaptive weight adjustment strategy is adopted to achieve optimal multi-source data fusion under different terrain regions and data quality, thereby generating enhanced point cloud data with high spatial consistency and complete features, improving the accuracy and robustness of terrain mapping.

[0016] (2) A three-dimensional terrain mapping method based on lidar provides a basis for terrain unit division by calculating multi-dimensional terrain features such as elevation variation, slope distribution, and curvature change. The terrain features are deeply fused with multispectral reflectivity and radar backscattering characteristics to generate a multimodal terrain classification map; and combined with boundary detection and anomaly analysis methods, abnormal areas are identified and an association model between abnormal features and data fusion parameters is established, providing closed-loop support for subsequent data optimization and processing, and improving the method's adaptability and practicality.

[0017] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0018] Figure 1 This is a flowchart of the three-dimensional terrain mapping method based on lidar of the present invention. Detailed Implementation

[0019] This application's embodiments solve the problems of insufficient accuracy in multi-source sensor data fusion, uneven point cloud quality, and incomplete extraction of terrain features by using a lidar-based three-dimensional terrain mapping method.

[0020] The overall approach of the scheme in this application embodiment is as follows: acquire airborne lidar point cloud data, simultaneously collect multispectral imagery and synthetic aperture radar data, and perform spatiotemporal correction and standardization processing on the multi-source data; based on point cloud density and multi-source feature identification, establish a multimodal registration mapping relationship, generate a point cloud quality assessment map, and adaptively adjust the fusion weights to form enhanced point cloud data; extract terrain features such as elevation, slope, and curvature from the enhanced point cloud data, divide the terrain into units, and generate a terrain classification map by combining multimodal features; finally, perform boundary anomaly detection and correlation analysis between anomaly features and fusion parameters to generate an anomaly spatial distribution map, achieving closed-loop analysis and optimization.

[0021] Please see Figure 1 This invention provides a technical solution: a three-dimensional terrain mapping method based on lidar, comprising the following steps: S1: acquiring raw point cloud data scanned by an airborne lidar, and simultaneously collecting reflectivity data from multispectral images and backscattering coefficients from synthetic aperture radar; identifying feature matching points in the point cloud data based on the correlation between point cloud density distribution and multi-source features, establishing a registration mapping relationship for multi-modal data, and generating a point cloud quality assessment map; S2: adjusting the fusion weight configuration of multi-source data based on the density distribution characteristics of the point cloud quality assessment map, and testing the feature representation of data from different sensors. S3: Determine the optimal fusion ratio for each data source to generate enhanced point cloud data with spatial consistency characteristics; S4: Extract terrain features from the enhanced point cloud data, including calculating elevation variation features, slope distribution features, and curvature change features. Divide the terrain into units based on the spatial distribution patterns of the features, and establish a correspondence between spectral and backscattering features to obtain a multimodal fused terrain classification map; S5: Based on the category distribution features in the terrain classification map, detect abnormal change areas at the boundaries of terrain units, analyze the correlation between abnormal areas and multi-source data fusion parameters, and generate an abnormal spatial distribution map.

[0022] In this implementation scheme, S1: Acquiring raw point cloud data from airborne lidar scanning refers to using a lidar sensor mounted on a UAV to scan the ground surface along its flight path, obtaining a three-dimensional coordinate point set (i.e., point cloud). Simultaneously acquired multispectral images refer to visible and near-infrared images in different bands, used to obtain spectral information of ground features; the backscattering coefficient of Synthetic Aperture Radar (SAR) reflects the scattering characteristics of the ground surface to radar waves, which can help distinguish different ground surface types. Based on this, by analyzing the correlation between point cloud density distribution (i.e., the distribution of point cloud quantity per unit space) and multi-source features, feature matching points are identified, i.e., the same ground feature feature points corresponding to different data sources. Further establishing the registration mapping relationship for multimodal data refers to aligning the lidar point cloud, multispectral images, and radar data in a unified coordinate system to ensure consistent spatial correspondence between different data sources, ultimately generating a point cloud quality assessment map to quantify the completeness, accuracy, and feature consistency of the point cloud data. S2: Based on the point cloud quality assessment map and point cloud density distribution characteristics, the multi-source data is fused with weights. The so-called fusion weight refers to the proportion of each data source's contribution to the final point cloud when merging LiDAR, imagery, and radar data. The optimal fusion ratio is determined by testing the consistency of feature representation across different sensor data, i.e., comparing the feature similarity of different sensor data at the same spatial location. The enhanced point cloud data generated using this method maintains spatial consistency while integrating the advantages of different sensors, improving the accuracy of subsequent terrain analysis. S3: Based on the enhanced point cloud data, terrain features are extracted, including elevation variation features (describing fluctuations in surface elevation), slope distribution features (describing the degree of surface tilt), and curvature variation features (describing changes in surface curvature). Based on the spatial distribution patterns of these features, the surface is divided into different terrain units, such as hills and plains. Further, by combining the spectral information of multispectral imagery and SAR backscattering characteristics, a correspondence between terrain units and multimodal features is established, resulting in a multimodal fused terrain classification map, where each terrain unit contains both morphological features and spectral and radar feature information. S4: Based on the category distribution features in the terrain classification map, abnormal change areas at the boundaries of terrain units are detected. Anomaly zones refer to spatial areas where significant differences in terrain or feature indicators occur at their boundaries. Further analysis of the correlation between these anomaly zones and the multi-source data fusion parameters is needed; specifically, to examine whether the anomalies are caused by weight allocation or registration errors during data fusion, and to generate an anomaly spatial distribution map. This process can provide a basis for subsequent optimization of fusion strategies or terrain analysis, improving the reliability and traceability of the method.

[0023] Specifically, the process of acquiring raw point cloud data from airborne lidar scanning and simultaneously collecting reflectivity data from multispectral images and backscattering coefficients from synthetic aperture radar is as follows: A multi-source data acquisition parameter mapping table is established to record the spatial attitude parameters and environmental parameters of each sensor during acquisition; radiometric calibration and atmospheric correction are performed on the multispectral image data to convert digital quantization values ​​into surface reflectivity data; incident angle correction and terrain radiometric correction are performed on the synthetic aperture radar data to eliminate the influence of terrain undulations on the backscattering coefficient; and multispectral reflectivity data and radar backscattering coefficients are registered to the spatial location of the lidar point cloud through spatial interpolation.

[0024] In this implementation plan, a multi-source data acquisition parameter mapping table is established. This table is generated by recording the spatial attitude parameters (including airborne platform attitude angles, pitch angles, and yaw angles) and environmental parameters (such as atmospheric transmittance, humidity, and illumination conditions) of each sensor during acquisition. The mapping table is used for subsequent spatiotemporal correction and registration of the multi-source data, ensuring consistency of data from different sensors within a unified reference coordinate system. This step does not involve complex calculations and can be described in words. Multispectral image radiometric calibration and atmospheric correction convert the digital quantized values ​​of the multispectral image... Converted to surface reflectance The process can be represented as: in, The surface reflectance of the m-th pixel in the k-th band; The digital quantization value of the m-th pixel in the k-th band; G k The gain coefficient of the k-th band can be obtained by comparing it with a standard reflector; B k The bias coefficient for the k-th band is obtained through dark current measurement; T k (L m Atmospheric transmittance correction function, dependent on pixel L m The observation path length can be estimated using an atmospheric radiative transfer model. By linearly correcting the digitally quantized values ​​and compensating for atmospheric effects, accurate calculation of surface reflectivity is achieved. Synthetic aperture radar incident angle and topographic radiative correction, radar backscattering coefficient... The correction can be expressed as: in, The corrected backscattering coefficient of the nth radar pixel; Original radar backscattering coefficient; θ n The angle of incidence can be calculated from the platform attitude and terrain elevation; H n : The terrain elevation corresponding to the nth pixel; f h (H nThe topographic radiation correction function is used to eliminate the influence of topographic relief on the scattering coefficient, and can be calculated by establishing a local slope-scattering mapping table. Spatial registration of multi-source data will convert multispectral reflectance... With radar backscattering coefficient Registration to the LiDAR point cloud location is achieved using a spatial interpolation algorithm. Let P be the coordinates of the p-th point in the LiDAR point cloud. p =(X p ,Y p Z p The corresponding multispectral reflectance and radar scattering coefficient are mapped using nearest neighbor interpolation or weighted interpolation as follows: in, The multispectral reflectance of the k-th band corresponding to lidar point p; The radar backscattering coefficient corresponding to point p in the lidar system; The set of nearest neighbors of point p in a multispectral image; The set of nearest neighbor pixels of point p in radar data; w q ,w r Weighting coefficients, used for weighted interpolation, are determined by inverse distance weighting, i.e., w. q =1 / d pq And normalize, d pq Let p be the spatial distance between point p and pixel q. Through the above steps, the spatiotemporal alignment and quality correction of multi-source data are completed, providing basic data for subsequent feature extraction and point cloud fusion.

[0025] Specifically, based on the correlation between point cloud density distribution and multi-source features, the process of identifying feature matching points in point cloud data, establishing registration mapping relationships for multimodal data, and generating point cloud quality assessment maps is as follows: Based on point cloud density distribution characteristics, salient feature regions are identified; texture and spectral features of the corresponding regions are extracted from multispectral images; scattering characteristics are extracted from synthetic aperture radar data; a correspondence between point cloud features and multi-source features is established using a feature similarity measurement algorithm; spatial transformation parameters between multi-source data are calculated; based on the spatial transformation parameters, the local geometric error, registration residual, density consistency, and multi-source feature consistency of each point cloud point are mapped into multi-dimensional quality indicators, and a three-dimensional quality assessment grid is constructed; the assessment grid is divided into regions based on density distribution characteristics, the quality score of each region is calculated, and a point cloud quality assessment map is generated.

[0026] In this implementation scheme, regions with significant features are identified based on the point cloud density distribution. The local density of the lidar point cloud is used to calculate the value of each point C. m Density index ρ m : Where, ρ m Point C m Local point cloud density; At point C m The set of points within the neighborhood; V s : Neighborhood volume; salient regions are defined by setting a density threshold τ p The threshold for identification is determined by adaptive calculation based on the mean and standard deviation of the entire point cloud density distribution. in, Overall average density of point cloud; σ ρ : Standard deviation of point cloud density; k: Empirical coefficient, which can be obtained through cross-validation optimization. Multi-source feature extraction: Extract texture features T from multispectral images corresponding to areas with significant density. n and spectral characteristics S n ; Extract scattering characteristic features Γ from synthetic aperture radar data. n This step involves data feature extraction, achieved through standard image processing and radar scattering analysis. Feature correspondences are established, and spatial transformation parameters are calculated using a feature similarity measurement algorithm, transforming point cloud features (X...) into... m ,Y m Z m ,ρ m ) and multi-source features (T n ,S n ,Γ n Perform matching and calculate the spatial transformation parameter H: Where H: spatial transformation matrix, containing rotation, translation, and scale parameters; F m : Feature vector of point m in point cloud; M m : Corresponds to multi-source feature vectors; G(H,M) m ): The mapping of multi-source feature vectors after transformation matrix H; d(·,·): Feature similarity measurement function, which can use weighted Euclidean distance; The weight coefficients are determined by the feature information entropy method in the similarity measurement, that is, for each feature component f k Calculate the information entropy H(f) k The higher the information entropy of a feature component, the greater its weight, ensuring that key features have a more significant impact on matching. Mapping point cloud quality metrics and constructing a 3D mesh will account for the local geometric error E of each point cloud point. m Registration residual R m Density consistency D m Consistency with multi-source features C m Mapped to a multidimensional quality index Q m Q m =α·E m +β·R m +γ·D m +δ·C mWhere α, β, γ, δ: weights of each quality indicator, determined by calculating the variance contribution rate of each indicator based on principal component analysis (PCA); when constructing the 3D mesh, the point cloud space is divided into uniform voxels, and the quality indicator of each voxel is the weighted average of the indicators of all points within that voxel. Among them, Q v Voxel unit quality indicators; voxel interior point set; w m Point weight coefficients are determined based on point density and feature importance. The evaluation grid region division and quality map generation process divides the 3D grid into regions based on density distribution characteristics, calculates the comprehensive quality score for each region, and generates a point cloud quality evaluation map through color coding or grayscale mapping, facilitating subsequent data fusion and processing.

[0027] Specifically, the process of adjusting the fusion weight configuration of multi-source data based on the density distribution characteristics of the point cloud quality assessment map is as follows: A quality-weight mapping function is established based on the multi-dimensional quality indicators in the point cloud quality assessment map, and the contribution of each quality dimension to the data fusion effect is analyzed; regions with different terrain complexities are identified based on density distribution characteristics, and a dynamic weight adjustment strategy is adopted for high-density variation regions, while a fixed weight configuration is adopted for low-density uniform regions; through an adaptive weight optimization algorithm, the optimal fusion weight of each sensor data at different spatial locations is calculated, and the weight ratio of regions with registration errors or inconsistent features is reduced, ultimately generating a multi-source data fusion weight configuration scheme.

[0028] In this implementation plan, a quality weight mapping function is established based on the multidimensional quality index Q in the point cloud quality assessment map. m (e.g., local geometric errors, registration residuals, density consistency, multi-source feature consistency), establish a mapping function W s =f(Q) m This is used to convert quality metrics into fusion weights. Among them, W s Q: The initial fusion weight of the s-th sensor in a specific grid cell; m : The m-th quality index value; η r The contribution coefficient of the quality indicator is determined based on the normalization of the indicator variance. That is, the larger the indicator variance, the more significant the impact on the fusion weight. Smoothing parameters to prevent the denominator from being zero. Divide the terrain complexity regions based on density distribution characteristics and apply this to the point cloud density distribution ρ. i Spatial analysis is performed, dividing the space into high-density variation regions and low-density uniform regions. A dynamic weight adjustment strategy is applied to the high-density variation regions, while a fixed weight configuration is applied to the low-density uniform regions. The specific steps are as follows: High-density variation regions: Initial weight W... sBased on this, a local adjustment coefficient λ is introduced. u The formula is as follows: W s '=W s ·(1+λ u ·δ ρ ); where W s ': Adjusted fusion weights; λ u : Dynamically adjust coefficients based on the local density change gradient δ ρ =ρ max -ρ min Adaptive determination; δ ρ : Local density variation range. Low-density uniform region: Maintain initial weight W s The adaptive weight optimization algorithm, in order to further consider registration error and feature consistency, optimizes the fusion weights at each spatial location, defining a weighted objective function F(W'). s ):F(W' s )=∑ s W' s ·(1-E s )+κ·C s Among them, E s Local registration error of sensor s; C s : Sensor feature consistency index; κ: Balance coefficient, used to adjust the contribution of error and consistency to the weights, determined by selecting the optimal fusion result through local cross-validation; through constraint condition ∑ s W' s =1 and 0≤W' s The optimal set of fusion weights is obtained by solving for values ​​≤1. A multi-source data fusion weight configuration scheme is generated, and the optimized weights are applied. This is applied to the corresponding spatial grid cells to form the final fusion weight configuration scheme, which is used to generate enhanced point cloud data. This step is the mapping execution process.

[0029] Specifically, the process of determining the optimal fusion ratio of each data source and generating enhanced point cloud data with spatial consistency features by testing the consistency of feature representations of different sensor data is as follows: A multimodal feature consistency evaluation framework is constructed, and data consistency is evaluated by calculating the differences in feature representations of different sensor data at the same spatial location; a multi-objective optimization algorithm is used, with the optimization objectives of maximizing feature consistency and minimizing geometric error, and the optimal fusion ratio of each data source is determined through Pareto optimal solution search; an enhanced point cloud dataset is generated based on the optimal fusion ratio using a weighted fusion algorithm.

[0030] In this implementation plan, a multimodal feature consistency evaluation framework is constructed for each spatial location P. v Point cloud data is used to calculate feature representation vectors from different sensors. (For example, features such as elevation, reflectivity, and backscattering coefficient), through a feature similarity measurement function. Calculate the consistency between sensor g and sensor h: Among them, S gh The spatial positions P of sensors g and h v Consistency of features; Sensors g and h at spatial position P v Eigenvectors on; ||·||2: Euclidean norm calculation difference; The smoothing parameter to prevent the denominator from being zero is determined by the smallest non-zero eigenvalue. A comprehensive consistency index for this spatial location is obtained by averaging across all sensor pairs. Among them, S v Spatial location P v The overall characteristic consistency index; M: total number of sensors. Multi-objective optimization determines the optimal fusion ratio, defining the fusion ratio R for each sensor. u Through a multi-objective optimization algorithm, S is maximized based on feature consistency. v and geometric error minimization E v To optimize the objective: max∑ v S v (R u ),min∑ v E v (R u ); where R u The fusion ratio of sensor u can be initially determined through uniform distribution; E v (R u ): Spatial location P of the fused point cloud v Local geometric errors; determined by Pareto optimal solution search for each sensor. The optimal fusion ratio is determined. An enhanced point cloud dataset is generated based on this determined optimal fusion ratio. The point cloud features at each spatial location are weighted and synthesized using a weighted fusion algorithm: in, Enhanced point cloud at spatial location P v The final feature vector; Sensor u at spatial position P v eigenvectors; The final fusion ratio of sensor u is determined by the results of multi-objective optimization. Fusion ratio The Pareto optimal solution search is used to determine a balance between feature consistency and geometric error for the entire point cloud dataset; the feature consistency index S gh Used to quantify the degree of matching of feature representations from different sensors; local geometric error E v (Ru The distance from a point to a neighboring triangular mesh can be calculated; the weighted fusion method linearly superimposes the feature vector for each spatial point to ensure that data from different sensors contribute in an optimized proportion.

[0031] Specifically, the process of extracting terrain features from enhanced point cloud data, including calculating elevation variation features, slope distribution features, and curvature change features, is as follows: The enhanced point cloud data is spatially gridded, and a local neighborhood point set is constructed within each grid cell. The elevation variation features of each grid cell are calculated, including the local elevation mean, elevation standard deviation, and elevation range. A local terrain surface is fitted using the least squares method, and the angle between the surface normal vector and the vertical direction is calculated as the slope value. The slope distribution of all points within the grid is then statistically analyzed. Curvature change features, including maximum curvature, minimum curvature, and average curvature, are calculated based on the residuals from the quadratic surface fitting, and all feature values ​​are stored in the point cloud attributes.

[0032] In this implementation scheme, spatial meshing and local neighborhood construction divide the enhanced point cloud data into a three-dimensional mesh according to a predefined resolution, with each mesh cell G... x It contains several point sets Q y ={p1,p2,…,p N}; where p z =(X z ,Y z Z z The coordinates of a point in the point cloud are represented by ). This local neighborhood allows for statistical analysis of the terrain features of each grid cell. Elevation variation characteristics are calculated in each grid cell G. x Internal calculated average elevation H m Standard deviation H s With range H r : H r =max(Z) z )-min(Z z ); where Z z : Elevation value of the z-th point; N: Number of points in the grid cell; H m : Local mean elevation, used to reflect the local topographic baseline; H s Local elevation standard deviation, used to describe the degree of elevation dispersion; H r Local elevation range, used to measure the magnitude of elevation change. Slope characteristic calculation for each grid cell G. x The local point set Q within y Fit the quadratic surface using the least squares method: Z = aU 2 +bV 2+cUV+dU+eV+f; where U,V: local plane coordinates; Z: elevation; a,b,c,d,e,f: fitting surface coefficients, determined by the least squares method; and the fitting surface normal vector. and vertical direction Calculate the slope θ using the included angle: in, Local fitting surface normal vector; Vertical unit vector; θ: slope value; the slope values ​​of all points within the grid can be used to statistically analyze slope distribution characteristics, such as maximum slope, minimum slope, and average slope. Curvature variation characteristics are calculated based on the quadratic surface fitting results, calculating the maximum curvature K. max Minimum curvature K min and mean curvature K avg : Where A = 2a, B = c, C = 2b: determined by the coefficients of the quadratic term of the fitted surface; K max Maximum curvature, reflecting the degree of curvature in the steepest local direction; K min Minimum curvature reflects the degree of curvature in the gentlest local direction; K avg Mean curvature, used to describe the overall curvature trend of local terrain. Feature storage stores the elevation variation feature H of each grid cell. m H s H r Slope characteristics θ and curvature characteristics K max ,K min ,K avg The attribute table stored in the enhanced point cloud data is used for subsequent terrain unit division and multimodal fusion analysis. The grid resolution is determined by statistical analysis of the point cloud density, ensuring that each grid cell contains a sufficient number of points on average to guarantee the stability of feature statistics; the least squares fitting coefficients are calculated using standard linear algebra methods; in curvature calculation, if there are singularities in the quadratic surface fitting, the triangular patch method within the local minimum neighborhood is used to estimate the curvature.

[0033] Specifically, the process of dividing terrain units according to the spatial distribution patterns of features and establishing a correspondence between spectral and backscattering features to obtain a multimodal fused terrain classification map is as follows: Multi-dimensional terrain features are spatially clustered using a clustering algorithm, and the point cloud is divided into different terrain units based on feature similarity; the mean of multispectral reflectance and the variance of radar backscattering coefficient are calculated for each terrain unit to establish a correspondence table between terrain type and multi-source features; a weighted fusion algorithm is used to integrate terrain geometric features and multi-source features, optimizing terrain unit boundaries; and the fusion results are used to label the categories of each terrain unit, generating a multimodal fused terrain classification map.

[0034] In this implementation scheme, terrain unit spatial clustering enhances the multi-dimensional terrain feature vector F of point cloud data. m =[H m H s H r ,θ,K max ,K min ,K avg The input clustering algorithm divides the point cloud into U terrain units through spatial clustering: Where, C v S: The point set of the v-th terrain unit; v : The set of point indices within the v-th terrain unit; Point p w Multidimensional terrain feature vectors; The feature mean vector of the v-th unit; clustering uses K-means or density-based spatial clustering algorithms, with the number of clusters U adaptively determined based on point cloud density and feature variance. Statistical analysis of multi-source features is performed on each terrain unit C. v The point set within which the mean multispectral reflectance R is statistically analyzed. v With the variance of radar backscattering coefficient Σ v : in, Point p w Multispectral reflectance; Point p w The radar backscattering coefficient; The average scattering coefficient within the v-th unit; |C v |: Number of unit points; These statistical values ​​are used to establish a correspondence table between terrain units and multi-source features, providing a basis for subsequent category labeling. Multimodal feature weighted fusion is applied to each terrain unit C. v Geometric features F v With multi-source features [R v ,Σ v Weighted fusion: G v =λ v ·F v +(1-λ v )·[R v ,Σ v ]; among them, G v : Feature vector of the fused terrain unit; λ v The fusion weighting coefficient is defined as follows: σ R : Standard deviation of multispectral characteristics; σ B Radar backscattering feature standard deviation; this weighting coefficient is adaptively determined by comparing the stability of each feature dimension to ensure the reliability of the fused features. Terrain unit category labeling is based on the fused feature vector G. vFor each terrain unit, a category is determined, using nearest neighbor classification or decision tree methods to label the terrain unit into different categories L. v This generates a multimodal fusion terrain classification map. Boundary optimization, performed after category labeling, optimizes the category boundaries between adjacent cells by minimizing the boundary inconsistency function: E. boundary =∑ (v,u) ||G v -G u || 2 ·δ(L v ≠L u ); where (v,u): adjacent cell index pair; δ(·): indicator function, which is 1 if the condition is met, and 0 otherwise; boundary optimization can reduce the discontinuity of terrain cell boundaries and improve the spatial coherence of the classification map.

[0035] Specifically, based on the category distribution characteristics in the terrain classification map, the specific process for detecting abnormal change areas at the boundary of terrain units is as follows: extract the boundary lines of each category region in the terrain classification map, enhance the boundary features using the Canny edge detection algorithm, and calculate the differences in geometric features such as elevation difference, slope difference, and curvature difference between the point clouds on both sides of the boundary; identify abnormal change points that exceed the normal range by setting an adaptive threshold, verify the abnormal points by combining multispectral and radar features, cluster continuous abnormal points to form abnormal change areas, and record their spatial location and feature attributes.

[0036] In this implementation scheme, boundary lines are extracted and boundary features are enhanced for each category region C in the terrain classification map. m Extract its boundary line B m The Canny edge detection algorithm is used to enhance the boundary, resulting in an enhanced boundary map ε. m B m Category Area C m The set of boundary points; ε m : Boundary map enhanced by the Canny algorithm; Standard deviation σ of Gaussian filtering in the Canny algorithm g It can be determined through local point cloud density statistics, that is, areas with high point cloud density use a smaller σ. g Preserve boundary details; use a larger σ value for areas with low point cloud density. g Smoothing noise. Calculating the geometric feature differences on both sides of the boundary for boundary points. Construct the set of neighboring points on both sides and Calculate the mean of the geometric features respectively: in, The average elevation of neighboring points on both sides of the boundary point; The average slope of the neighboring points on both sides of the boundary point; Mean curvature of neighboring points on both sides of the boundary point; ΔH n ,ΔS n ,ΔK n Geometric feature differences at boundary points are used to quantify the degree of abnormal changes. Adaptive thresholding identifies outliers by setting an anomaly threshold T. H ,T S ,T K Calculate abnormal indicators: Among them, T H ,T S ,T K Adaptive thresholds for elevation, slope, and curvature; the thresholds are obtained by statistically analyzing the local feature distribution of boundary points within the same category region, for example... in This is the local average elevation. Let be the local elevation standard deviation, and λ be an adjustable amplification factor used to control anomaly sensitivity. The set of anomaly points A = {a} is identified by combining multi-source features to verify the anomaly point pairs. n Extract its corresponding multispectral reflectance R. n and radar backscattering characteristics B n The retention of outliers is verified based on multi-source feature consistency rules. If a point's multi-source features deviate from the normal range, it is confirmed as an outlier; otherwise, it is excluded. The rule can use Mahalanobis distance or standardized difference metric to calculate the feature deviation. Clustering is then used to form anomalous regions. Continuous or adjacent outliers are spatially clustered to generate anomalous change regions Z. k It also records the spatial location and characteristic attributes of each region. Clustering algorithms can employ density-based DBSCAN, with parameters... Neighborhood distance and MinPts (minimum number of points) are adaptively determined by the average spacing of the point cloud and the density of outlier distribution. The output of the spatial distribution results of outliers generates a thematic map containing the boundaries, spatial location, and characteristic attributes of the outlier regions, which can be used for subsequent analysis and multi-source data fusion optimization.

[0037] Specifically, the process of analyzing the correlation between abnormal regions and multi-source data fusion parameters to generate anomaly spatial distribution maps is as follows: extract the multi-source data fusion parameters corresponding to the abnormal regions, including the weight configuration and registration accuracy parameters of each sensor; calculate the correlation coefficient between the feature values ​​of the abnormal regions and the fusion parameters, analyze the influence of parameter settings on the anomaly detection results, generate anomaly intensity distribution maps through spatial interpolation, visualize the spatial distribution patterns of abnormal regions through heatmaps, and generate thematic maps containing anomaly types, spatial distribution characteristics, and parameter correlation analysis.

[0038] In this implementation scheme, the fusion parameters corresponding to the abnormal regions are extracted for the detected abnormal region A. q (Obtained from step S4), extract the fusion weight W of each sensor data in this area. r W s W t And the registration accuracy parameter P x ,P y ,P z :W r : Fusion weights of LiDAR data in abnormal areas; W s : Fusion weights of multispectral image data in anomalous regions; W t : Fusion weights of synthetic aperture radar data in anomalous regions; P x ,P y ,P z The local registration accuracy index of the point cloud in the X, Y, and Z directions is used to quantify the geometric error of the point cloud. The weighting coefficients are determined by the adaptive weight optimization algorithm in the previous stage, i.e., calculated based on the feature consistency and density distribution of each sensor data in the local neighborhood. The correlation coefficient between the anomaly features and the fusion parameters is calculated for the feature values ​​F within the anomaly region. a =[H a ,θ a ,K a ] and fusion parameters [W r W s W t ,P x ,P y ,P z Calculate the correlation coefficient matrix R ab : in, Abnormal region A q The terrain feature vector of the p-th point within the area; Abnormal region A q The fusion parameter vector corresponding to the q-th point; Cov(·): covariance function; Standard deviation of topographic features; The standard deviation of the fusion parameters; the correlation coefficient is used to analyze the influence of each fusion parameter on the anomaly detection results. An anomaly intensity distribution map is generated, which represents the anomaly characteristic intensity S of the anomaly region. a Correlation R with fusion parameters ab Mapping to a 3D mesh space via spatial interpolation: Where, S a (x m ,y m ,z m ): Mesh cell (x m ,y m ,z m The abnormal intensity value of S;ak : Anomaly eigenvalues ​​of neighboring point k; d km : Euclidean distance from neighboring point k to the center of the grid cell; γe: distance attenuation coefficient, determined through cross-validation to ensure smooth spatial interpolation while preserving local differences; N m The set of neighboring points around a grid cell. Visualization of the spatial distribution of anomalies, using heatmaps or 3D visualization, to represent the anomaly intensity value S. a (x m ,y m ,z m The map is mapped to a color gradient, and combined with the correlation coefficient matrix of category labels and fusion parameters, to generate thematic maps: - showing anomaly types; - showing the spatial distribution patterns of anomalies; - showing the degree of influence of fusion parameters on anomalies; this visualization supports subsequent anomaly analysis and parameter optimization.

[0039] In summary, this application has at least the following effects:

[0040] This LiDAR-based 3D terrain mapping method constructs a multi-source data acquisition parameter mapping table and performs targeted radiometric and terrain correction on LiDAR, multispectral imagery, and synthetic aperture radar data. This achieves high-precision registration of multi-source data in both spatial and feature dimensions, providing a reliable quality foundation for subsequent data fusion and terrain analysis. Utilizing point cloud quality assessment maps and point cloud density distribution characteristics, the fusion weights of each sensor's data are adaptively adjusted. The optimal fusion ratio is determined through multimodal feature consistency evaluation, generating enhanced point cloud data with spatial consistency characteristics, effectively improving data fusion accuracy and robustness. Multi-dimensional terrain features, such as elevation variation, slope distribution, and curvature changes, are extracted from the enhanced point cloud data. Terrain units are divided based on the spatial distribution patterns of these features. Furthermore, multimodal correspondences are established by combining multispectral and radar features, achieving refined and spatially coherent terrain classification. By analyzing the differences in boundary geometric features, determining adaptive thresholds, and verifying multi-source features, abnormal change areas at the boundaries of terrain units are identified. A correlation model between abnormal areas and multi-source data fusion parameters is established, generating an anomaly spatial distribution map. This achieves a closed loop of anomaly detection, visualization, and cause analysis, providing decision support for subsequent terrain mapping and data processing.

[0041] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0042] This invention is described with reference to flowchart illustrations and / or block diagrams of systems, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0043] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0044] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0045] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0046] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A three-dimensional terrain mapping method based on lidar, characterized in that, Includes the following steps: S1: Acquire raw point cloud data from airborne lidar scanning, and simultaneously collect reflectivity data from multispectral images and backscattering coefficients from synthetic aperture radar. Based on the correlation between point cloud density distribution and multi-source features, identify feature matching points in the point cloud data, establish registration mapping relationships for multimodal data, and generate point cloud quality assessment maps. S2: Based on the density distribution characteristics of the point cloud quality assessment map, adjust the fusion weight configuration of multi-source data, determine the optimal fusion ratio of each data source by testing the consistency of feature representation of data from different sensors, and generate enhanced point cloud data with spatial consistency characteristics. S3: Extract terrain features from enhanced point cloud data, including calculating elevation variation features, slope distribution features, and curvature change features. Divide terrain units according to the spatial distribution pattern of the features, and establish a correspondence by combining spectral and backscattering features to obtain a multimodal fusion terrain classification map. S4: Based on the category distribution characteristics in the terrain classification map, detect abnormal change areas at the boundaries of terrain units, analyze the correlation between abnormal areas and multi-source data fusion parameters, and generate an abnormal spatial distribution map.

2. The three-dimensional terrain mapping method based on lidar according to claim 1, characterized in that: The specific process of acquiring raw point cloud data from airborne lidar scanning and simultaneously collecting reflectivity data from multispectral images and backscattering coefficients from synthetic aperture radar is as follows: Establish a multi-source data acquisition parameter mapping table to record the spatial attitude parameters and environmental parameters acquired by each sensor. Radiometric calibration and atmospheric correction are performed on multispectral image data to convert digital quantization values ​​into surface reflectance data. Incidence angle correction and topographic radiometric correction are performed on synthetic aperture radar data to eliminate the influence of topographic undulation on the backscattering coefficient. Multispectral reflectance data and radar backscattering coefficient are registered to the spatial location of lidar point cloud through spatial interpolation.

3. The three-dimensional terrain mapping method based on lidar according to claim 2, characterized in that: Based on the correlation between point cloud density distribution and multi-source features, the specific process of identifying feature matching points in point cloud data, establishing registration mapping relationships for multimodal data, and generating point cloud quality assessment maps is as follows: Based on the point cloud density distribution characteristics, we identify regions with significant features, extract texture and spectral features of the corresponding regions from multispectral images, and extract scattering characteristics from synthetic aperture radar data. The correspondence between point cloud features and multi-source features is established through a feature similarity measurement algorithm, and the spatial transformation parameters between multi-source data are calculated. Based on the spatial transformation parameters, the local geometric error, registration residual, density consistency, and multi-source feature consistency of each point cloud point are mapped into multi-dimensional quality indicators to construct a three-dimensional quality assessment grid. The evaluation grid is divided into regions based on density distribution characteristics, the quality score of each region is calculated, and a point cloud quality evaluation map is generated.

4. The three-dimensional terrain mapping method based on lidar according to claim 1, characterized in that: The specific process for adjusting the fusion weight configuration of multi-source data based on the point cloud quality assessment map and density distribution characteristics is as follows: Based on the multidimensional quality indicators in the point cloud quality assessment map, a quality-weight mapping function is established to analyze the contribution of each quality dimension to the data fusion effect. Based on density distribution characteristics, regions with different terrain complexities are identified. A dynamic weight adjustment strategy is adopted for high-density variation regions, while a fixed weight configuration is adopted for low-density uniform regions. Calculate the optimal fusion weights of each sensor's data at different spatial locations, reduce the weight ratio of regions with registration errors or inconsistent features, and generate a multi-source data fusion weight configuration scheme.

5. The three-dimensional terrain mapping method based on lidar according to claim 4, characterized in that: The specific process of generating enhanced point cloud data with spatial consistency features by testing the consistency of feature representations from different sensor data sources and determining the optimal fusion ratio of each data source is as follows: A multimodal feature consistency evaluation framework is constructed to evaluate data consistency by calculating the differences in feature representation of data from different sensors at the same spatial location. The optimal fusion ratio of each data source is determined by using a multi-objective optimization algorithm with the optimization objectives of maximizing feature consistency and minimizing geometric error, and by searching for the Pareto optimal solution. An enhanced point cloud dataset is generated using a weighted fusion algorithm based on the optimal fusion ratio.

6. The three-dimensional terrain mapping method based on lidar according to claim 1, characterized in that: The specific process for extracting terrain features from enhanced point cloud data, including calculating elevation variation features, slope distribution features, and curvature variation features, is as follows: The enhanced point cloud data is spatially gridded, a local neighborhood point set is constructed within each grid cell, and the elevation variation characteristics of each grid cell are calculated, including the local elevation mean, elevation standard deviation, and elevation range. The local terrain surface is fitted using the least squares method. The angle between the surface normal vector and the vertical direction is calculated as the slope value, and the slope distribution of all points within the grid is statistically analyzed. Curvature variation characteristics are calculated based on the residual of quadratic surface fitting, including maximum curvature, minimum curvature and average curvature, and the feature values ​​are stored in point cloud attributes.

7. The three-dimensional terrain mapping method based on lidar according to claim 6, characterized in that: The specific process of dividing terrain units based on the spatial distribution patterns of features, establishing a correspondence between spectral and backscattering features, and obtaining a multimodal fused terrain classification map is as follows: Spatial clustering of multi-dimensional terrain features is performed using clustering algorithms, and point clouds are divided into different terrain units based on feature similarity. For each terrain unit, the mean of its multispectral reflectance and the variance of its radar backscattering coefficient are statistically analyzed, and a table of correspondence between terrain type and multi-source features is established. By integrating topographic geometric features with multi-source features through a weighted fusion algorithm, the boundaries of topographic units are optimized, and the categories of each topographic unit are labeled according to the fusion results to generate a multimodal fused topographic classification map.

8. The three-dimensional terrain mapping method based on lidar according to claim 1, characterized in that: Based on the category distribution characteristics in the terrain classification map, the specific process for detecting abnormal change areas at the boundaries of terrain units is as follows: Extract the boundary lines of each category region in the terrain classification map, enhance the boundary features using the Canny edge detection algorithm, and calculate the geometric differences in the point clouds on both sides of the boundary, such as the elevation difference, slope difference, and curvature difference. By setting an adaptive threshold, abnormal change points that exceed the normal range are identified. The abnormal points are verified by combining multispectral and radar features. Continuous abnormal points are clustered into abnormal change regions, and their spatial location and characteristic attributes are recorded.

9. The three-dimensional terrain mapping method based on lidar according to claim 8, characterized in that: The specific process for analyzing the correlation between abnormal regions and multi-source data fusion parameters to generate an anomaly spatial distribution map is as follows: Extract the multi-source data fusion parameters corresponding to the abnormal regions, including the weight configuration and registration accuracy parameters of each sensor; The correlation coefficient between the feature values ​​of the abnormal region and the fusion parameters is calculated, the influence of parameter settings on the anomaly detection results is analyzed, an anomaly intensity distribution map is generated through spatial interpolation, and the spatial distribution pattern of the abnormal region is visualized through heat map, generating thematic maps that include anomaly type, spatial distribution characteristics and parameter correlation analysis.

Images