A mountainous terrain modeling method based on multi-point cloud data fusion
By using a multi-point cloud data fusion method and employing techniques such as binocular cameras, unscented Kalman filtering, and GAN algorithms, a high-precision 3D model of mountainous terrain is generated. This solves the problems of low cost and poor visualization effect in existing technologies and is applicable to fields such as UAV surveying and architectural surveying.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- SOUTH CHINA AGRICULTURAL UNIVERSITY
- Filing Date
- 2023-06-27
- Publication Date
- 2026-05-19
AI Technical Summary
Existing technologies lack low-cost and visually appealing methods for modeling mountain terrain, making it difficult to directly obtain high-precision mountain terrain data, and vegetation cover has a significant impact.
A multi-point cloud data fusion method is adopted. Reference point cloud data is acquired through a stereo camera. Data cleaning and registration are performed by combining unscented Kalman filtering, GAN algorithm and ICP algorithm. GAN algorithm is used to fill in data gaps. Finally, a 3D model is generated by piecewise linear interpolation and Savitzky-Golay filtering.
It generates high-precision 3D mountain models in variable environments, reduces errors, and improves data accuracy. It is applicable to fields such as UAV surveying and architectural surveying, providing an accurate 3D data foundation.
Smart Images

Figure CN116758234B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of 3D modeling technology, and more specifically to a method for modeling mountainous terrain based on multi-point cloud data fusion. Background Technology
[0002] Topographic modeling is a technique that uses computer technology to digitally model the Earth's surface topography. Its main purpose is to create digital topographic scenes for applications such as topographic analysis and visualization. The data acquisition process for topographic modeling typically involves multiple methods, including surveying instruments, aerial photography, and even satellite remote sensing. DEM data acquisition can directly obtain elevation information from ground data sources, converting the elevation values into model vertices, and then using interpolation algorithms to transform the data into a 3D model. However, such data acquisition requires high-precision instruments, which are expensive and may be cumbersome and inconvenient. Remote sensing technology is easily affected by vegetation cover, making it difficult to directly obtain mountain topography data. Therefore, to meet the needs of real-world production, it is necessary to design a low-cost yet accurate topographic mapping and modeling method to provide data support for topographic analysis, mountain development planning, planting plans, irrigation design, and visualization of topographic features and impacts.
[0003] Currently, there is a lack of low-cost and visually appealing modeling methods for mountainous terrain. Modeling methods that rely on direct interpolation for output are insufficient to provide a reference for real-world applications. We need terrain modeling methods that have a higher degree of fit to reality and better visualization.
[0004] Therefore, improving the realism of terrain modeling and its visualization effects are problems that urgently need to be solved by those skilled in the art. Summary of the Invention
[0005] In view of this, the present invention provides a mountain terrain modeling method based on multi-point cloud data fusion, which improves the degree of fitting with reality, provides better visualization effects, and provides data support for terrain analysis, mountain development planning, planting plans, irrigation design, visualization of terrain features and terrain impacts.
[0006] To achieve the above objectives, the present invention adopts the following technical solution:
[0007] A method for modeling mountainous terrain based on multi-point cloud data fusion includes the following steps:
[0008] Step 1: Collect depth photos of the mountainous area and generate reference point cloud data B;
[0009] Select appropriate angles and calibration points, use a binocular camera to acquire depth photos of the mountain from four directions, and generate reference point cloud data B based on these photos; reference point cloud data B mainly includes the coordinate information of the point cloud converted from depth information;
[0010] Step 2: Collect multiple sets of experimental data from the same area and perform unscented Kalman filtering to obtain experimental point cloud data A1;
[0011] Multiple sets of experimental data samples were acquired from the same mountain using sensors, and the original sensor data were subjected to unscented Kalman filtering to reduce noise, resulting in experimental point cloud data A1. Experimental point cloud data A1 is sensor data, mainly including three-dimensional coordinates, velocity, and acceleration.
[0012] Step 3: Clean the experimental point cloud data A1 to obtain cleaned point cloud data A2. Based on the cleaned point cloud data A2, construct an outlier data identification model using the GAN algorithm. Use the outlier data identification model to mark the outliers in the cleaned point cloud data A2. Perform an inverse transformation on the outliers to obtain the outlier point cloud data.
[0013] Step 4: Filter out valid data from outlier point cloud data and obtain the axial coordinate information of valid data;
[0014] Based on outlier point cloud data A Outliers The distribution in a two-dimensional coordinate system is used to filter out valid data from one axis with distorted data but valid coordinates on other axes, thereby obtaining outlier point cloud data A. Outliers Axial coordinate information of valid data in the middle Outliers_useful ;
[0015] Step 5: Filter out the outlier point cloud data in the experimental point cloud data A1 to obtain filtered point cloud data A4. Use the ICP algorithm to register the reference point cloud data B and the filtered point cloud data A4 to obtain registered point cloud data AB.
[0016] Filtering GANs 1.0 Labeled outlier point cloud data A Outliers The filtered point cloud data A4 is obtained. The reference point cloud data B generated from the depth image is registered with the cleaned filtered point cloud data A4 using the ICP algorithm to obtain the registered point cloud data AB.
[0017] Step 6: Based on the axial coordinate information of the registered point cloud data AB and the effective data, construct a blank data repair model using the GAN algorithm. Use the blank data repair model to correct the filtered point cloud data A4 to obtain the filled point cloud. Then, fuse the filled point cloud and the filtered point cloud data A4 to obtain the terrain point cloud A5.
[0018] The registered point cloud data AB is used as a constraint and the effective axial coordinate information A from the outlier data. Outliers_useful As input, a GAN trained using the GAN algorithm to fill in blank data is generated. 2.0 The model utilizes a pre-trained GAN.2.0 Generator 2.0 The missing parts of the experimental data after cleaning are generated, and the generated point cloud is A. Generator2.0 Discriminator 2.0 Supervise the generated point cloud, A Generator2.0 The terrain point cloud A5 is obtained by fusing it with A4;
[0019] Step 7: Apply piecewise linear interpolation to interpolate the terrain point cloud A5, and obtain the 3D terrain model after filtering.
[0020] Preferably, in step 1, a mountain landmark is selected, and a binocular camera calibrated using the Zhang Zhengyou calibration method is used to acquire mountain depth photos in four directions from the mountain landmark. The binocular camera is calibrated using the Zhang Zhengyou calibration method, and pixel coordinates in the photos are extracted by collecting photos of different poses using a calibration board composed of two-dimensional squares. The intrinsic parameter matrix is output through the homography matrix, and the parameters are optimized using the maximum likelihood estimation method. The distortion is removed using normalized coordinates to obtain reference point cloud data B.
[0021] The requirements for selecting a mountain landmark are as follows:
[0022] 1. Select points with distinctive features. This includes specific landmarks on the mountain, such as mountain peaks, rocks, and trees, or prominent terrain features, such as ridges, valleys, and slopes.
[0023] 2. Even Distribution: Calibration points should be evenly distributed across the four directions of the mountain. This ensures that accurate depth information can be obtained from photos taken from different locations.
[0024] 3. Consider distance and altitude. Calibration points should be located at different distances and altitudes between the camera and the mountain to help the camera measure depth more accurately and reduce errors caused by changes in viewing angle.
[0025] 4. Artificial Structures. Calibration points can be artificial structures within the mountain, such as buildings, fences, roads, etc.
[0026] 5. Safety and Accessibility. Calibration points should be located in safe and accessible locations for camera calibration and depth measurement.
[0027] 6. Consider environmental factors. Lighting and weather conditions affect camera performance and the quality of depth images. Adjust acquisition strategies and data processing methods based on these factors to ensure that the obtained depth information is accurate and reliable.
[0028] Preferably, the mountainous multi-group experimental data samples obtained in step 2 are mainly obtained by collecting laser distance to the ground and GPS positioning through sensors to obtain latitude and longitude coordinate output. The latitude and longitude coordinate output is a TXT format point cloud file containing coordinates in the x, y, and z directions, which is used as experimental data. The experimental data is subjected to unscented Kalman filtering to obtain experimental point cloud data A1.
[0029] Preferably, step 2 employs unscented Kalman filtering to address the nonlinearity issues arising during sensor sampling, primarily including:
[0030] Step 21: Determine state variables and observation variables: When processing point cloud data, state variables can be selected as information such as the x, y, and z coordinates, velocity, and acceleration of the point cloud. Assume the state at time t is x... t The observed variables are the experimental data.
[0031] Step 22: Initialize state variables and state covariance matrix: Before starting filtering, the state variables and state covariance matrix need to be initialized. Usually, you can choose to set the state variables to the x, y, z coordinates of the point cloud, and set the state covariance matrix to a large diagonal matrix.
[0032] The initial values are estimated based on empirical data and measurement errors, or the predicted state covariance matrix output after the previous unscented Kalman filter is used as the initial state matrix for the next unscented Kalman filter.
[0033] Step 23: Select the sigma point for the unscented Kalman filter:
[0034] Step 231: Calculate the predicted state transition matrix F t and the predicted state covariance matrix P t The formula is as follows:
[0035]
[0036]
[0037] Where Q is the process noise covariance matrix, and P is Gaussian white noise; t-1 Let f represent the predicted state covariance matrix at time t-1; T represents the transpose; 0 indicates that the process noise is a zero vector; f is the nonlinear dynamic function of the system.
[0038] Step 232: Perform unscented Kalman filtering on the state variables to obtain 2n+1 state sigma points, including one center point and 2n extension points, for subsequent state estimation, as shown in the following formula:
[0039]
[0040]
[0041]
[0042] Where λ is the adjustment parameter of the unscented Kalman filter algorithm, which is used to control the distribution of sigma points. It is usually taken as 3-n, where n is the dimension of the defined state vector. P represents the square root of the i-th eigenvalue of the matrix; t Let X be the state at time t. t The predicted state covariance matrix, Pt, is obtained by predicting the state transition matrix F. t It is obtained by combining the covariance matrix of the predicted state at the previous time step;
[0043] At time t, the predicted state transition matrix is F t The predicted covariance matrix is P. t The process noise covariance matrix Q follows a Gaussian distribution; at time t, Here, T stands for transpose operation; taking the first round of unscented Kalman filtering as an example, P0 is the initial Kalman filter, i.e., the initial state covariance matrix in step 22. The next time step after P0 is P1, then the prediction covariance matrix calculated in step 23 is... Therefore, the predicted covariance matrix in step 23 is obtained by predicting the current state and the system model, and is the predicted state covariance matrix considering state transition and process noise.
[0044] Step 24: State Estimation: Perform a nonlinear transformation on the sigma points to achieve state transition, obtaining 2n+1 predicted state sigma points, represented as:
[0045]
[0046] in, Let f represent the predicted state sigma point, i.e., the sigma point after nonlinear transformation; f is the nonlinear dynamic function of the system, since state X t It mainly consists of kinematic physical quantities, and this dynamic function can be obtained based on Newton's laws of motion;
[0047] Step 25: Calculate the new predicted state mean and predicted state covariance matrix: Based on the predicted state sigma points and their corresponding weights, calculate the new predicted state mean. and the predicted state covariance matrix P t∣t-1 ; indicates
[0048]
[0049]
[0050] The predicted state covariance matrix is the predicted state covariance matrix P at time t-1. t-1 The calculations show the uncertainty of the system state.
[0051] Because the sigma points are chosen to approximate the density of the system state, these selected points, after being advanced to the next time step through the state transition function, yield the predicted state sigma points. The predicted state sigma points represent: considering system dynamics and process noise W... t After that, the distribution of the system state. However, even if the selected sigma points are advanced to the next time step, and even if system factors and noise are considered, the predicted state sigma is still not enough to directly obtain the predicted state mean and predicted state covariance matrix. Therefore, after generating sigma points and advancing them from time t to time t+1 through the state transition function, it is necessary to recalculate the new predicted state mean and predicted state covariance by combining the predicted state sigma points and weights.
[0052] Step 26: Calculate the predicted observation mean and covariance matrix: Using the predicted state mean and predicted state covariance matrix, calculate the predicted observation mean. and the predicted observation covariance matrix S t , is represented as:
[0053]
[0054]
[0055] in, It is a prediction of the observed sigma point. and These are the weighting coefficients, and Q is the process noise covariance matrix;
[0056] The predicted observation covariance matrix is the predicted observation covariance matrix at time t obtained by using the predicted state covariance matrix at time t-1. It reflects the uncertainty of the prediction of the observation results in the absence of new observation data, and takes into account the uncertainty of previous state predictions and observation models.
[0057] Step 27: Calculate the Kalman gain: Based on the predicted observation mean and predicted observation covariance matrix, and the observed variables (i.e., the actual observed values), calculate the Kalman gain K. t That is, the cross-covariance matrix C between the predicted state and the predicted observation. t With the predicted observation covariance matrix S t The quotient; represented as:
[0058]
[0059]
[0060] C t This is the cross-covariance matrix between the predicted state and the predicted observations; Predict the state sigma point; It is the predicted state mean; It is to predict the observed sigma point; It is a prediction of the observed mean; These are the weighting coefficients.
[0061] Step 28: Update the new predicted state mean and predicted state covariance matrix: Use Kalman gain to update the new predicted state mean and the predicted state covariance matrix from step 23, and then use the new observation data to update the new predicted state covariance matrix.
[0062] Step 29: Iteration: Repeat steps 22-28, initializing the state variables and state covariance matrix according to the updated new predicted state mean and predicted state covariance matrix, until the convergence condition is met or the preset number of iterations is reached. Output the estimated state variables and covariance matrix, as well as the filtered point cloud data, i.e., the experimental point cloud data A1. The obtained predicted observation mean is the filtered experimental point cloud data.
[0063] Preferably, the specific implementation process of step 3 is as follows:
[0064] Step 31: Perform radius filtering coarse cleaning on the experimental point cloud data A1 to obtain the cleaned point cloud data A2;
[0065] Step 32: Perform plane fitting on the cleaned point cloud data A2, and filter the fitted plane according to the set threshold ΔX to obtain the key point cloud A3;
[0066] Select experimental point cloud data that are less than the set threshold ΔX, that is, select points with small deviations near the fitting plane and mark them as key point cloud A3;
[0067] Step 33: Normalize the cleaned point cloud data A2 and key point cloud A3 respectively, and convert them into voxel data A3 using a voxelization tool. 21 And voxel data A 31 ;
[0068] Step 34: Based on voxel data A 31 Constructing an outlier data identification model using GAN algorithm 1.0 ;
[0069] Outlier Data Identification Model GAN 1.0This refers to an adversarial neural network trained using a GAN network, which includes a Generator (denoted as Generator). 1.0 ) and Discriminator (denoted as Discriminator) 1.0 );
[0070] Step 35: Use outlier data to identify the GAN model 1.0 Discriminator in 1.0 For voxel data A 21 Supervise the process, identify and label outlier data, and inversely transform it into outlier point cloud data A containing x, y, and z coordinates. Outliers ;
[0071] Combine Generator 1.0 It can generate pseudo point clouds while Discriminator 1.0 The Discriminator is trained under supervision. 1.0 This can be used for cleaned point cloud data A2, using Discriminator. 1.0 It serves as a classifier to identify potential outliers and real data in A2's point cloud data.
[0072] Preferably, in step 4, the outlier point cloud data A Outliers The point cloud is transformed from two-dimensional data into point distributions on the xy-axis, xz-axis, and yz-axis planes. A discrete point cloud that does not overlap on only one axis plane and maintains good continuity is considered valid data, characterized by data distortion in one direction but still valid in the other two. The axis coordinate information of this valid data is then backed up to A. Outliers_useful For example: Suppose that a point cloud a does not coincide with any other point on any of the axes, and the straight-line distance between point cloud a and its neighboring point clouds on the xy-axis plane is less than a set outlier threshold, while the straight-line distance between point cloud a and its neighboring point clouds on the other two xz-axis planes and yz-axis plane is greater than or equal to the set outlier threshold, then the coordinates of point cloud a on the x and y axes are valid, but the data of point cloud a on the z-axis is distorted.
[0073] Preferably, the specific implementation process of step 5 is as follows:
[0074] Step 51: Delete the outlier point cloud data in the experimental point cloud data A1 to obtain the filtered point cloud data A4;
[0075] Step 52: Select the initial transformation matrix T0 based on the calibration points marked in the mountainous area, and pre-align the filtered point cloud data A4 with the reference point cloud data B;
[0076] Step 53: Use the KD-tree method to process each point i in the reference point cloud data B.a The corresponding point cloud i was found in the filtered point cloud data A4. b ;
[0077] Step 54: Minimize the point cloud i using the least squares method. a And Point Cloud i b The Euclidean distance error between them is used to obtain the transformation matrix T1;
[0078] Step 55: Refer to the point cloud data B, where each point cloud corresponds to a transformation matrix, and select the optimal transformation matrix T. n ;
[0079] Step 56: After registration, the filtered point cloud data A4 is transformed by the optimal transformation matrix T. n The transformation process yields a point cloud that is as aligned as possible with the reference point cloud data B, denoted as the registered point cloud data AB. The transformation matrix adjusts the spatial position and orientation of point cloud A4, making it as aligned as possible with point cloud B. During this process, the coordinates of each point in the original point cloud A4 change after rotation, translation, and other transformations. Therefore, the registered point cloud A4 (AB) differs from the original point cloud A4 in spatial position. Although the goal of registration is to align point cloud A with point cloud B as closely as possible, because A and B may have different structures (i.e., their point distribution and shape may differ), even with optimal registration, the registered point clouds A and B may not be completely identical.
[0080] Preferably, in step 6, since the binocular camera selects four-directional shooting of the mountain to mainly obtain the mountain's shape, it is less affected by the mountain's vegetation. Therefore, the reference point cloud data B can be used as a reference scale to provide a second GAN algorithm training. The GAN network trained in the second training is a blank data repair model GAN. 2.0 This also includes the Generator (denoted as Generator). 2.0 ) and Discriminator (denoted as Discriminator) 2.0 Discriminator for blank data imputation model 2.0 With the constraints of reference point cloud data B, the generator for repairing blank data in the model can be accelerated. 2.0 The convergence of Generato 2.0 It can be used for predictive point generation, Discriminator 2.0 Then, supervision is performed, and the axial coordinate information A of the valid data is transferred. Outliers_useful Make input, Generato 2.0 Predicting data A that was originally distorted in a certain direction Generator As filling point cloud A Generator2.0To fill the gaps in the experimental point cloud data A1 caused by data cleaning, A Generator2.0 The terrain point cloud A5 is obtained by fusing it with the filtered point cloud data A4.
[0081] Preferably, piecewise cubic sampling interpolation is used. Piecewise linear sampling can effectively improve accuracy, and this sampling method can effectively avoid Runge's phenomenon. The specific implementation process of step 7 is as follows:
[0082] Step 71: Sort the terrain point cloud A5 according to the size of its x-coordinate;
[0083] Step 72: Divide the terrain point cloud A5 into several segments based on the distance between adjacent points in the sorted terrain point cloud A5 and the set distance threshold; each segment contains multiple adjacent points, and the x-coordinates of these points do not differ by more than the distance threshold.
[0084] Step 73: For each point cloud segment, perform interpolation using the cubic linear interpolation method;
[0085] For each segment, a cubic linear interpolation is performed. For each point P(i) in the point cloud of each segment, a cubic function is fitted using the two preceding points P(i-1) and P(i-2) and the two following points P(i+1) and P(i+2), a total of five points. This cubic function passes through the three points P(i-1), P(i), and P(i+1), and simultaneously satisfies the continuity of the first derivative at P(i-2) and P(i+2). This ensures the continuity of the interpolation curve between segments and also allows the interpolation curve to better adapt to the local characteristics of each segment, reflecting the original terrain characteristics of the terrain point cloud A5 as much as possible.
[0086] Step 74: Connect the interpolation results of each point cloud segment to obtain the entire curve, i.e., the interpolation curve.
[0087] Preferably, the Savitzky-Golay filter is used to fit the interpolation curve. That is, the linear least squares method is used to fit a continuous subset of adjacent data points in the interpolation curve with a low-order polynomial. This maintains accuracy while smoothing the point cloud data, thereby improving the accuracy and stability of the terrain model and obtaining a three-dimensional terrain model.
[0088] As can be seen from the above technical solution, compared with the prior art, the present invention discloses a mountain terrain modeling method based on multi-point cloud data fusion. Under variable environments, it obtains a three-dimensional model with mountain scale information, which can realize the digitization of mountain terrain and provide data support and guidance for mountain development and mountain orchard planting planning. It starts from two key directions: improving the modeling accuracy of discontinuous discrete point cloud and avoiding model distortion caused by the loss of key point cloud after filtering. It ensures that a relatively accurate three-dimensional model can still be obtained with as few sampling times as possible, even when the sampling process is variable and the vegetation coverage is large and UAV operation is not applicable. Specifically, through steps such as registration, data imputation, and curve fitting, errors are effectively reduced and data accuracy is improved, making the results in fields such as UAV mapping and architectural mapping more accurate and reducing the difficulty of subsequent work. The processed point cloud data can better reflect three-dimensional objects in the real world, aiding in subsequent analysis and applications, achieving accurate representation and understanding of terrain, and can be used in fields such as urban planning and disaster assessment. The processing flow uses diverse point cloud data, applicable to multiple fields such as UAV mapping, architectural mapping, autonomous driving, and robot guidance. The organic combination of multiple point clouds provides a more accurate three-dimensional data foundation for related technologies. In this invention, multiple algorithms complement each other and work together to realize a complete three-dimensional point cloud data processing flow. Each algorithm plays a crucial role in the entire processing flow and is indispensable. Attached Figure Description
[0089] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.
[0090] Figure 1 This is a schematic diagram of the modeling process provided by the present invention;
[0091] Figure 2 This is a schematic diagram of the unscented Kalman filter algorithm provided by the present invention. Detailed Implementation
[0092] 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.
[0093] This invention discloses a method for mountain terrain modeling based on multi-point cloud data fusion. The method involves acquiring an initial depth image of the mountain, calibrating initial points and the horizontal plane, fusing data collected from multiple sensors to obtain spatial location information data samples, performing data fusion using an unscented Kalman filter algorithm, training a GAN algorithm to identify and remove discrete point clouds, registering key point clouds based on the acquired initial image, generating and filling in data gaps caused by the removed point clouds using the GAN algorithm, interpolating and meshing the sample data, and smoothing the model using Savitzky-Golay filtering to establish a three-dimensional model of the mountain.
[0094] Example 1
[0095] A method for modeling mountainous terrain based on multi-point cloud data fusion, such as Figure 1 The diagram shows the modeling method flowchart, which includes the following steps:
[0096] S1: Select appropriate angles and calibration points, use a binocular camera to acquire depth photos of the mountains in four directions, and generate reference point cloud data B based on these photos;
[0097] S2: Multiple sets of experimental data samples of the same mountain are acquired through sensors, and the original sensor data are subjected to unscented Kalman filtering to reduce noise, resulting in experimental point cloud data A1.
[0098] S3: Radius filtering is applied to the experimental point cloud data A1. After coarse cleaning by radius filtering, the resulting experimental point cloud data A2 is subjected to plane fitting. Points on the plane below a set value ΔX are selected as key point clouds A3. A2 and A3 are normalized and further converted into voxel data A2 using a voxelization tool. 21 And the voxel data A3 of key point cloud 31 Combining the GAN algorithm with voxel data A 31 Use the input as the basis to build a GAN for training outlier data identification. 1.0 Model, trained GAN 1.0 Discriminator in the algorithm 1.0 For voxel data A 21 Supervise the identification and labeling of outlier data, and inversely transform them into outlier point cloud data A containing x, y, and z coordinates. Outliers ;
[0099] S4: Based on outlier data A Outliers In a two-dimensional coordinate system, the data is filtered out by identifying data with distorted coordinates on one axis while valid coordinates on other axes, thereby obtaining valid axis coordinate information A from the outlier data. Outliersuseful ;
[0100] S5: Filter out outlier data labeled GAN1.0 Outliers The cleaned point cloud data A4 is obtained. The reference point cloud data B generated from the depth image is registered with the cleaned experimental point cloud data A4 using the ICP algorithm to obtain the registered point cloud data (AB).
[0101] S6: Use the registered point cloud data (AB) as a constraint and the valid axial coordinate information A from the outlier data. Outliers_useful As input, a GAN trained using the GAN algorithm to fill in blank data is generated. 2.0 The model utilizes a pre-trained GAN. 2.0 Generator 2.0 The missing parts of the experimental data after cleaning are generated, and the generated point cloud is A. Generator2.0 Discriminator 2.0 Supervise the generated point cloud, A Generator2.0 The terrain point cloud A5 is obtained by fusing it with A4;
[0102] S7: The piecewise linear interpolation method is applied to interpolate the terrain point cloud A5, and then the Savitzky-Golay filter is used to smooth the interpolated 3D model to improve the accuracy and stability of the terrain model.
[0103] Example 2
[0104] In one specific embodiment, based on the modeling method flow of Embodiment 1, further,
[0105] S1: Obtain depth photos of the mountain. Prioritize the four directions of the mountain and take photos of each direction with an angle of not less than 45 degrees to capture the entire mountain as much as possible. This will provide reference point cloud data for subsequent repair of blank areas caused by noise filtering.
[0106] When acquiring depth images, Zhang Zhengyou's calibration method is used for calibration. This requires taking more than 10 photos of the printed checkerboard from different directions. By changing the orientation of the checkerboard multiple times, the image is captured to enrich the coordinate information. The intrinsic parameter matrix is output through the homography matrix, and camera intrinsic parameters such as focal length, principal point, and distortion coefficient are calculated. The maximum likelihood estimation method is used for parameter optimization, and then the image coordinates are converted into world coordinates. However, this calibration method only considers radial distortion and does not consider tangential distortion. Therefore, OpenCV's iterative distortion correction algorithm is used to normalize and correct the distortion of the point coordinates.
[0107] After obtaining the depth image, it is converted into point cloud information for reference. The steps are as follows: Take two depth images of the same location of the mountain from different angles.
[0108] S11: Read the depth image; use the OpenCV library function cv2.imread() to read the depth image;
[0109] S12: Align the two depth images;
[0110] S13: Select one of the two images as a reference image for matching;
[0111] S14: Use the cv2.matchTemplate() function to find the region in another image that is most similar to the reference image;
[0112] S15: Use the cv2.minMaxLoc() function to find the coordinates of the best matching region;
[0113] S16: Use the cv2.findHomography() function to calculate the transformation matrix between the two images to align them;
[0114] S17: Convert the depth image to a point cloud; convert the depth value of each pixel into 3D coordinates in the camera coordinate system to obtain a point cloud; define camera parameters, including the camera intrinsic matrix and distortion parameters. Define the transformation matrix between the depth camera's coordinate system and the world coordinate system; this transformation matrix is usually determined by the camera's position and attitude, and the camera's attitude information can be obtained using the sensor's IMU (Inertial Measurement Unit) or other external sensors.
[0115] S18: Convert the depth value of each pixel to 3D coordinates in the camera coordinate system. Assuming the depth camera resolution is W*H, for each pixel (i,j), it can be converted to 3D coordinates (x,y,z) in the camera coordinate system using the following formula:
[0116]
[0117]
[0118] Where cx and cy are the principal point coordinates in the camera intrinsic parameter matrix, fx and fy are the focal lengths in the x-axis and y-axis directions in the camera intrinsic parameter matrix, respectively, and the z coordinate value is the depth value of pixel (i,j).
[0119] S19: Convert the 3D coordinates in the camera coordinate system to the 3D coordinates in the world coordinate system to obtain the reference point cloud coordinate information and perform Gaussian filtering. This set of point cloud data will be referred to as reference point cloud data B below.
[0120] S2: The experimenters carried sensor modules to collect data on the mountain. Before the data collection, the sensor modules needed to be reset and calibrated using a level and ruler, and the initial point needed to be located. The collected information included the laser distance to the ground and the latitude and longitude coordinates obtained from GPS positioning. The calculation of the distance to the ground required the collection of attitude and acceleration information. The attitude of the sensor modules was corrected to improve the accuracy of the laser distance. Combined with the experimenters' displacement trajectory, the x, y, and z coordinates were calculated, and the point cloud coordinates were output.
[0121] The changes in the state variables output by the sensor module and the relationship between the observed variables and the state variables are not linear. Unscented Kalman filtering is used to estimate the state of the system with nonlinear state and observation equations. The process is as follows: Figure 2 As shown;
[0122] Unscented Kalman filtering can be used to process point cloud data containing time series and x, y, z coordinate information. These point cloud data can be used to estimate the state variables such as position, velocity, and acceleration of multi-sensor devices during observation. As mentioned above, the problem of unscented Kalman filtering in processing point clouds can be modeled as a state space model, where the state variables include position, velocity, and acceleration, and the observation variables are point cloud data.
[0123] Specifically, the unscented Kalman filter algorithm consists of two parts: prediction and update. The specific definitions and calculation processes are as follows:
[0124] S21: Define variables:
[0125] Suppose the state at time t is:
[0126] x t =[p x ,p y ,p z ,v x ,v y ,v z ,a x ,a y ,a z ] T
[0127] Where p x ,p y ,p z These represent the positions in the x, y, and z directions at the time of observation, respectively, v x v y x z Let a represent the velocities in the x, y, and z directions at the time of observation. x a y a z These represent the accelerations in the x, y, and z directions at the time of observation, respectively.
[0128] The observed variable is point cloud data, that is:
[0129] z t =[x1,y1,z1,x2,y2,z2,…,x n ,y n ,z n ] T
[0130] Where n represents the number of points in the point cloud data, x i y i z i Represent the x, y, and z coordinates of the i-th point;
[0131] S22: The state-space model is represented as:
[0132] x t+1 =fx t ,w t
[0133] z t =h(x t ,v t )
[0134] Where f and h are the state transition function and the observation function, respectively, and w t and v t This includes process noise and observation noise.
[0135] Because the sensor has a high sampling frequency, and considering the actual sampling process, the motion between two adjacent point clouds can be regarded as linear motion. Based on Newton's laws of motion, the state transition function can be defined as follows:
[0136]
[0137] Observation function:
[0138] z t =hx t =[x t 0 3|n-1 ]+v t
[0139] Among them 0 3|n-1 The zero vector represents the n-1 points, meaning that points other than the first point do not affect the state variables;
[0140] To further explain why points other than the first point do not affect the state variables, the Jacobian matrix of the observation function with respect to the state is:
[0141]
[0142] Each row corresponds to the partial derivative of the observation function with respect to the state variable at a given point. The Jacobian matrix in Kalman filtering describes how the observation depends on the state variable; that is, it represents the degree to which the observation is affected by the state variable through the partial derivative of the observation function with respect to the state variable. More specifically, it reflects how much the predicted observation will change if the state variable changes by a single point. The Jacobian matrix is introduced here to explain "why points other than the first point do not affect the state variable." The state covariance matrix describes the uncertainty of the state variable. This matrix reflects the uncertainty in each direction of the defined dimensions, thereby allowing the calculation of the confidence interval for the state.
[0143] The Jacobian matrix and covariance matrix are two important matrices in state estimation, measuring different metrics. The Jacobian matrix measures the dependence of observations on the state, while the covariance matrix measures the uncertainty of state estimation.
[0144] In unscented Kalman filtering, the observation function of a point cloud is usually obtained by superimposing the observation functions of each point in the point cloud. If each point in the point cloud affects the estimation of the state variables, then the superimposed observation function will also affect the estimation of the state variables.
[0145] In the application, it is assumed that the initial first point affects the estimation of the state variables, while points other than the first point do not affect the estimation of the state variables. Therefore, all points in the point cloud except the first point are ignored. This allows us to obtain an observation function that only considers the first point, thus avoiding useless computation. Consequently, the partial derivatives for points other than the first point are zero, i.e.:
[0146]
[0147] achievable
[0148]
[0149] As shown in the Jacobian matrix above, the observation function of the point cloud does not affect the estimation of the state variables except for the first point because its derivative with respect to the Jacobian matrix of the state vector is 0.
[0150] Defining process noise and observation noise: Both process noise and observation noise are assumed to be Gaussian white noise, i.e.:
[0151] w t ~N0,Q
[0152] v t ~N(0,R)
[0153] Here, Q and R are the covariance matrices of process noise and observation noise, respectively. Assuming that they follow a Gaussian distribution, the mean can be regarded as zero. However, a mean of 0 does not mean that the noise is 0 at any time, nor does it mean that the noise change is 0.
[0154] S23: Prediction Steps:
[0155] Based on the current state x t and process noise w t Calculate the predicted state transition matrix F t and the predicted state covariance matrix P t :
[0156]
[0157]
[0158] Where 0 indicates that the process noise is a zero vector, and Q is the process noise covariance matrix;
[0159] By sampling the state variables using Kalman filtering at sigma points, we obtain 2n+1 sigma points:
[0160]
[0161]
[0162]
[0163] Where λ is an adjustment parameter in the unscented Kalman filter algorithm, used to control the distribution of sigma points, typically taking a value of 3-n, where n is the dimension of the state vector. P represents the square root of the i-th eigenvalue of the matrix. t To predict in state x t The covariance matrix.
[0164] S24: Perform state transitions on the sigma points to obtain 2n+1 predicted state sigma points:
[0165]
[0166] The sigma point represents the predicted state, i.e., the sigma point after nonlinear transformation. f is the nonlinear dynamic function of the system, since x is defined... t It mainly consists of kinematic physical quantities, and this dynamic function can be obtained based on Newton's laws of motion;
[0167] S25: Calculate the mean of the predicted state and predicted state covariance P t∣t-1 :
[0168]
[0169]
[0170] in, and Here, Q represents the weighting coefficients, and Q is the process noise covariance matrix.
[0171] and The specific calculations are as follows:
[0172]
[0173]
[0174]
[0175] Where α and β are parameters of the unscented Kalman filter algorithm, typically with a value of 0.1;
[0176] S26: Update steps:
[0177] S261: Sampling of the sigma points observed in the predicted state yields 2n+1 predicted observation sigma points:
[0178]
[0179] It predicts the observed sigma point, where h is the observation function. The sigma point represents the predicted state, i.e., the sigma point after the state transition. t To observe the noise, v t Size depends primarily on sensor accuracy;
[0180] S262: Calculate the predicted observation mean and its predicted observation covariance S t :
[0181]
[0182] in, The predicted mean is obtained by weighted summation of the predicted observation sigma points. As weight;
[0183]
[0184] S tThe predicted observation covariance is obtained by calculating the outer product of the differences between the predicted observation sigma point and the predicted observation value, and then performing a weighted sum. R represents the weights, and R is the observation noise covariance matrix.
[0185] S263: Calculate the covariance matrix C between the predicted state and the predicted observations. t (i.e., cross covariance):
[0186]
[0187] C t The cross covariance between the predicted state and the predicted observations; This represents the predicted state sigma point; It is the predicted state mean; It is to predict the observed sigma point; It is a prediction of the observed mean. As weight;
[0188] S264: Calculate Kalman gain:
[0189]
[0190] Among them, C t S is the covariance matrix between the predicted state and the predicted observation. t It predicts the observed covariance;
[0191] S265: Based on the observed values, calculate the corrected state value, i.e., the updated state covariance P. t|t :
[0192]
[0193] P t|t It is the updated state covariance, P t|t-1 It is the predicted state covariance, K t It is the Kalman gain, S t It is the predicted observation covariance.
[0194] S266: Calculate the updated state mean:
[0195]
[0196] To predict the state mean, Z t K represents the actual observed value. t For Kalman gain, To predict the observed mean.
[0197] S267: After data correction, continue with the prediction and update steps until the complete set of point cloud data has been processed.
[0198] The above describes the point cloud processing procedure using the unscented Kalman filter algorithm.
[0199] S3: After the filtered 3D point cloud data A1 is cleaned by radius filtering, point cloud data A2 is obtained. A2 is then subjected to plane fitting. A threshold ΔX is set based on the distance from the point cloud to the plane, and points with small deviations are extracted as key point clouds A3. A2 and A3 are normalized and further converted into voxel data A of the experimental point cloud data A2. 21 And the voxel data A3 of key point cloud 31 This can be achieved using commonly used voxelization tools;
[0200] Voxel data is used to generate voxel maps, which can be achieved by using view rendering or point cloud to depth map methods, thus representing 3D point cloud data on a 2D plane; voxel data is preprocessed, such as by resizing and scaling, so that it can be used as input to the discriminator of the GAN algorithm;
[0201] The Generative Adversarial Network (GAN) algorithm was used to identify outliers. Because GAN algorithms have a Generator and a Discriminator, the first model built using the GAN algorithm was named GAN. 1.0 Contains Generator 1.0 and Discriminator 1.0 ;
[0202] The GAN algorithm's Discriminator can perform binary classification on the input data, namely, real data and pseudo data. Therefore, the preprocessed voxel data A 31 The input is fed into the GAN algorithm for the first training, resulting in the corresponding GAN. 1.0 Generator 1.0 and Discriminator 1.0 The voxel data A2 of the experimental point cloud data A2 21 As input, using Discriminator 1.0 By performing identification supervision, a label can be obtained, which can be marked as real data or pseudo data. Voxel data marked as pseudo data are regarded as outlier data and can be marked as outliers.
[0203] The Generator can estimate the data gaps in a mountain area caused by the cleaned-up data after inputting a specific vector.
[0204] S4: Discriminator 1.0 After identifying and labeling outlier data, the neighboring points of each outlier point cloud are also labeled as nearby outliers. These points, along with the outlier data, are then converted into point clouds and output as A. Outliers Because the acquisition of the x, y, and z coordinates of the point cloud relies on a multi-sensor module, and the values obtained by each sensor are independent, the x, y, and z coordinates of outliers will exhibit some distortion. Further screening can be performed to retain as much accurate information as possible. Outliers Projecting the data onto three two-dimensional planes—the xy-axis plane, the xz-axis plane, and the yz-axis plane—the following possibilities exist for the labeled outlier data on these three projected planes:
[0205] If the data in the xy-axis plane does not coincide with other point clouds and the straight-line distance to the adjacent point clouds is less than the set outlier threshold, but the difference between the data and the adjacent point clouds in the xz-axis plane and yz-axis plane exceeds the set outlier threshold, or is seriously deviated from most point clouds, it can be said that the position information of the marked outlier data in the horizontal xy-axis plane is valid, but the height information represented by the z-axis in the vertical height is distorted.
[0206] If, on the xz-axis plane, the outlier data of the outlier point cloud is not more than the set outlier threshold compared to the neighboring point clouds, but on the xy-axis plane, the outlier point cloud deviates from the neighboring point clouds more than the set outlier threshold, or overlaps with other point clouds, or is significantly deviated from most point clouds;
[0207] If the value of the outlier point cloud on the z-axis does not exceed the set threshold with the neighboring point clouds on the yz-axis plane, it can be said that the height information represented by the z-axis in the vertical height of the outlier point cloud is valid, while the position information of the x-axis or y-axis in the horizontal direction is distorted.
[0208] If the outlier point cloud deviates significantly from the referenced adjacent point cloud and most of the point cloud on all three two-dimensional axis planes, then the marked outlier data is completely distorted and can be directly removed; the axis plane coordinate information A of the valid data in the outlier data is obtained. Outliers_useful ;
[0209] S5: Filter out GANs 1.0 Labeled outlier data A Outliers The cleaned point cloud data A4 is obtained;
[0210] The reference point cloud data B and point cloud data A4 obtained from the depth image are correlated by selecting mountain features;
[0211] The ICP algorithm is used to register two sets of point clouds: an initial transformation matrix T0 is selected based on the correspondence, and the KD-tree method is used to register each point cloud i in the reference point cloud data. a The corresponding i was found in the experimental point cloud data. b Based on the least squares method to minimize the point cloud i a i b The transformation matrix T1 is calculated using the Euclidean distance error between the points. By continuously finding corresponding points and calculating new transformation matrices, the optimal solution is obtained, and the final transformation matrix T is output. n After registration, the cleaned point cloud data A4 is transformed by the optimal transformation matrix T. n The transformation yields a point cloud that is as aligned as possible with the reference point cloud data B, denoted as the registered point cloud data AB.
[0212] S6: Register the valid axial coordinate information A from the point cloud data AB and the outlier data. Outliers_useful As input, perform a second GAN algorithm training to build the GAN model. 2.0 This generates the second generation of Generator and Discriminator, namely the Generator. 2.0 and Discriminator 2.0 The Generator obtained at this time 2.0 The generated pseudo-point cloud, because the reference point cloud data is added as a constraint, is closer to the real data.
[0213] A Outliers_useful Input GAN 2.0 Then, by the Generator 2.0 For A Outliers_useful Fill-in generation, by Discriminator 2.0 Supervised judgment is performed to complete the missing parts as much as possible, resulting in modeling point cloud data; for example: A Outliers_useful A certain point cloud coordinate (x) n ,y n ,z n ), the z-mark n It is a severe outlier, x n y n If it is valid, then z will be generated by Generator 2.0. n1 The new z-axis coordinate of this point is denoted as (x n ,y n ,z n1 ). A Generator2.0 The terrain point cloud A5 is obtained by fusing it with A4;
[0214] S7: Perform piecewise cubic linear interpolation on the terrain point cloud data A5 that needs to be modeled and visualized. Assume the coordinates of the i-th point in the point cloud are (x... i ,y i ,z i ):
[0215] Sort the point cloud according to the x-coordinate. For the k-th segment of the point cloud, assume its starting point and ending point are labeled as i. k and j k Where ik <= jk, and a threshold d is defined such that if Then we consider point i k and point j k In the same paragraph;
[0216] Specifically, for the k-th point cloud segment, the interpolation function is:
[0217]
[0218] Where a k b k c k d k These are the four parameters that need to be fitted; in order to satisfy the interpolation function in The first derivative is continuous at , so the interpolation function needs to be divided into two quadratic functions, respectively at . and Connection at the point;
[0219] For the k-th point cloud segment, we need to solve for a. k b k c k d k These four parameters can be viewed as a linear least squares problem. Substituting the sample points P(i-2), P(i-1), P(i), P(i+1), and P(i+2) into the interpolation function, we obtain the following system of equations:
[0220]
[0221]
[0222]
[0223]
[0224]
[0225] Then a can be solved using the least squares method. k b k c k d k ;
[0226] By concatenating the interpolation functions of each point cloud segment, a cubic linear interpolation function for the entire point cloud is obtained. The new x value is then interpolated using the interpolation function to obtain the corresponding y and z values. Piecewise cubic linear interpolation can be performed separately in the x, y, and z dimensions to obtain three consecutive cubic interpolation functions, thus obtaining a three-dimensional continuous curve of the entire point cloud.
[0227] Furthermore, the `signal` module in the SciPy library is used to call the Savitzky-Golay filter to smooth the curve. Filter parameters need to be set, including window size and polynomial fitting order. Then, the `savgol_filter` function is called. After the function returns the smoothed data, it is visualized or used for subsequent analysis and verification, thus obtaining a noise-removed, trend-preserving, and distortion-free 3D model of the mountainous terrain.
[0228] Example 3:
[0229] The method of this invention can be applied to fields such as UAV surveying and architectural surveying. In practical applications, such as UAV laser surveying, the data acquired by the sensor includes the ground distance measured by laser and the latitude and longitude coordinates obtained by GPS positioning, outputting coordinates in the x, y, and z directions. The method includes the following steps:
[0230] (1) Data acquisition and preprocessing:
[0231] The drone cruises and collects raw point cloud data and acquires four-directional depth images of the target mountain.
[0232] Combining unscented Kalman filtering and training a GAN 1.0 The network performs preprocessing operations such as denoising, outlier classification, and registration on the collected raw point cloud data. During UAV mapping, sensor data may be affected by environmental interference, such as wind and rain. Therefore, in the data preprocessing stage, it is necessary to identify outliers and remove noise points.
[0233] By combining depth images, the ICP algorithm is used to achieve accurate registration between point clouds, reduce data errors, and provide higher-precision data input for subsequent processing.
[0234] (2) Data imputation: For the missing parts in the point cloud data, the trained GAN is used. 2.0 The model performs data completion. During UAV mapping, occlusions or sensor blind spots may occur, resulting in incomplete point cloud data. A generator produces point cloud data for the missing parts, while a discriminator identifies the authenticity of the generated data, ensuring high consistency between the generated and original data.
[0235] (3) Curve Fitting: The filled point cloud data is subjected to piecewise cubic linear interpolation. First, the point cloud data is divided into several segments, then cubic linear interpolation is performed on each segment, and finally the interpolation functions of each segment are spliced together to obtain a continuous curve of the entire point cloud. Since the UAV flight control is rigorous and orderly, the data acquisition process is more stable than land operation. This step can be performed in the three dimensions of x, y, and z to obtain three continuous cubic interpolation functions, further improving the integrity and continuity of the point cloud data.
[0236] (4) Smoothing: The Savitzky-Golay filter is used to smooth the curve. In practical applications, the terrain may have abrupt changes. To better reflect the actual terrain, smoothing is required. Set the filter parameters, including window size and polynomial fitting order. The smoothed data can be used for visualization or subsequent analysis and verification to obtain a noise-removed, trend-preserving, and distortion-free 3D model.
[0237] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the apparatus disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple; relevant parts can be referred to the method section.
[0238] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.
Claims
1. A method for modeling mountain terrain based on multi-point cloud data fusion, characterized in that, Includes the following steps: Step 1: Collect mountain depth photos of fixed mountain landmarks and generate reference point cloud data B; Step 2: Collect multiple sets of experimental data from the same area and perform unscented Kalman filtering to obtain experimental point cloud data A1; Step 3: Clean the experimental point cloud data A1 to obtain cleaned point cloud data A2. Based on the cleaned point cloud data A2, construct an outlier data identification model using the GAN algorithm. Use the outlier data identification model to mark the outliers in the cleaned point cloud data A2. Perform an inverse transformation on the outliers to obtain the outlier point cloud data. Step 4: Filter out valid data from outlier point cloud data and obtain the axial coordinate information of valid data; Step 5: Filter out the outlier point cloud data in the experimental point cloud data A1 to obtain filtered point cloud data A4. Use the ICP algorithm to register the reference point cloud data B and the filtered point cloud data A4 to obtain registered point cloud data AB. Step 6: Based on the axial coordinate information of the registered point cloud data AB and the effective data, construct a blank data repair model using the GAN algorithm. Use the blank data repair model to correct the filtered point cloud data A4 to obtain the filled point cloud. Then, fuse the filled point cloud and the filtered point cloud data A4 to obtain the terrain point cloud A5. Step 7: Apply piecewise linear interpolation to interpolate the terrain point cloud A5, and obtain the 3D terrain model after filtering.
2. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, In step 1, a mountain landmark is selected, and a binocular camera calibrated by Zhang Zhengyou's calibration method is used to acquire mountain depth photos in four directions from the mountain landmark. A calibration board composed of two-dimensional squares is used to collect photos of different poses to extract pixel coordinates from the photos. The intrinsic parameter matrix is output through the homography matrix, and the parameters are optimized using the maximum likelihood estimation method. The distortion is removed using normalized coordinates to obtain reference point cloud data B.
3. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, In step 2, the laser-to-ground distance and GPS positioning are collected by sensors to obtain a TXT format point cloud file containing latitude and longitude coordinates in the x, y, and z directions as experimental data. Unscented Kalman filtering is then applied to the experimental data to obtain experimental point cloud data A1. The specific process of unscented Kalman filtering during the data acquisition is as follows: Step 21: Determine the state variables and observed variables; Step 22: Initialize the state variables and the state covariance matrix to obtain the initialized state variables and the initialized covariance matrix; Step 23: Calculate the predicted state transition matrix based on the initialized state variables and the initialized state covariance matrix; Calculate the predicted state covariance matrix based on the predicted state transition matrix and the process noise covariance matrix; Unscented Kalman filtering sigma point sampling is performed on the state variables, the predicted state transition matrix, and the predicted state covariance matrix to obtain the state sigma points; Step 24: Use a nonlinear dynamic function to perform state transition on the state sigma point to obtain the predicted state sigma point; Step 25: Calculate the new predicted state mean and the new predicted state covariance matrix based on the predicted state sigma points and their corresponding weights. Step 26: Using the new predicted state mean and the new predicted state covariance matrix, calculate the predicted observation mean and the predicted observation covariance matrix; Step 27: Calculate the Kalman gain based on the predicted observation mean, the predicted observation covariance matrix, and the observed variables; Step 28: Update the new predicted state mean and the predicted state covariance matrix in step 23 using Kalman gain, and return to step 22 to initialize the state variables and state covariance matrix respectively based on the updated new predicted state mean and predicted state covariance matrix, until the convergence condition is met or the preset number of iterations is reached.
4. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 3, characterized in that, The specific process of step 23 is as follows: Step 231: Calculate the predicted state transition matrix and the predicted state covariance matrix based on the initialized state variables, nonlinear dynamic function, and process noise covariance matrix; Step 232: Sampling the sigma points of the state variables by performing unscented Kalman filtering based on the predicted state covariance matrix to obtain the state sigma points.
5. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, The specific implementation process of step 3 is as follows: Step 31: Perform radius filtering coarse cleaning on the experimental point cloud data A1 to obtain the cleaned point cloud data A2; Step 32: Perform plane fitting on the cleaned point cloud data A2, and filter the fitted plane according to the set threshold ΔX to obtain the key point cloud A3; Step 33: Normalize the cleaned point cloud data A2 and key point cloud A3 respectively, and convert them into voxel data A. 21 And voxel data A 31 ; Step 34: Based on voxel data A 31 Construct an outlier data identification model using GAN algorithm; Step 35: Use the Discriminator in the outlier identification model to analyze voxel data A 21 The system performs supervision, identifies and labels outlier data, and then inversely transforms it into outlier point cloud data containing x, y, and z coordinates.
6. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, In step 4, the outlier point cloud data is converted into a two-dimensional distribution of points on the xy-axis plane, xz-axis plane, and yz-axis plane. Points that do not overlap and remain continuous on only one axis plane are identified as valid data, which are distorted in one direction but remain valid in the other two directions. The axis coordinate information of the valid data is obtained.
7. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, The specific implementation process of step 5 is as follows: Step 51: Delete the outlier point cloud data in the experimental point cloud data A1 to obtain the filtered point cloud data A4; Step 52: Select the initial transformation matrix T0 based on the mountain landmarks, and pre-align the filtered point cloud data A4 with the reference point cloud data B; Step 53: Use the KD-tree method to process each point i in the reference point cloud data B. a The corresponding point cloud i was found in the filtered point cloud data A4. b ; Step 54: Minimize the point cloud i using the least squares method. a And Point Cloud i b The Euclidean distance error between them is used to obtain the transformation matrix T1; Step 55: Refer to the point cloud data B, where each point cloud corresponds to a transformation matrix, and select the optimal transformation matrix T. n ; Step 56: Filter point cloud data A4 and apply the optimal transformation matrix T n The transformation yields the point cloud aligned with the reference point cloud data B, denoted as the registration point cloud data AB.
8. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, The blank data inpainting model constructed in step 6 includes a Generator and a Discriminator. The Discriminator of the blank data inpainting model is subject to constraints and supervision based on the reference point cloud data B, while the Generator of the blank data inpainting model generates prediction points to obtain the inpainted point cloud.
9. The method for mountain terrain modeling based on multi-point cloud data fusion according to claim 1, characterized in that, Step 7 uses piecewise cubic sampling interpolation to interpolate the terrain point cloud A5. The specific implementation process is as follows: Step 71: Sort the terrain point cloud A5 according to the size of its x-coordinate; Step 72: Divide the terrain point cloud A5 into several segments based on the distance between adjacent points in the sorted terrain point cloud A5 and the set distance threshold; Each point cloud segment contains multiple adjacent point clouds, and the difference in x-coordinates between adjacent point clouds in each segment is less than or equal to a distance threshold. Step 73: For each point cloud segment, perform interpolation using the cubic linear interpolation method; Step 74: Connect the interpolation results of each point cloud segment to obtain the interpolation curve.
10. A method for modeling mountain terrain based on multi-point cloud data fusion according to claim 1, characterized in that, The interpolation results are filtered using a Savitzky-Golay filter to obtain a 3D terrain model.