Incremental 3D Reconstruction Method and System for Unmanned Aerial Vehicles in Complex Mountainous Environments

Abstract

This invention discloses an incremental 3D reconstruction method and system for unmanned aerial vehicles (UAVs) in complex mountainous environments, belonging to the field of 3D reconstruction technology. The method includes: acquiring multi-view monocular images; establishing feature correspondences between images; optimizing newly added images to construct a sparse 3D point cloud; initializing each 3D point in the sparse 3D point cloud as a 3D Gaussian primitive; based on a comparison of the normal change rate and a preset steepness threshold, making the Gaussian ellipsoid anisotropic or isotropic distribution; projecting and transforming the 3D Gaussian primitives and the 3D covariance matrix to generate 2D Gaussian distributed primitives, and weighted accumulating them to generate a 2D image; constructing a joint loss function to back-optimize the parameters of the 3D Gaussian primitives to obtain a static 3D Gaussian model; and generating 3D reconstruction results based on the static 3D Gaussian model. This invention enhances the robustness of feature detection and matching in complex mountainous environments by adaptively adjusting the scale parameters based on the slope angle and mean curvature.

Description

Technical Field

[0001] This invention belongs to the field of 3D reconstruction technology, and in particular relates to an incremental 3D reconstruction method and system for UAVs in complex mountainous environments. Background Technology

[0002] Currently, UAV 3D reconstruction technology mainly relies on traditional structured light, LiDAR point clouds, and multi-view geometry (SfM / MVS) algorithms. While these methods can achieve high accuracy in regular scenes or urban environments, they still have significant limitations in complex mountainous environments.

[0003] (1) Poor robustness of image matching: In mountain scenes, due to the presence of texture repetition (rock walls), large scale changes, drastic lighting changes, and many shadows, the stability of feature points is poor, making it difficult to match images.

[0004] (2) Severe pose drift in large-scale environments: The flight path in mountainous areas is long and the field of view is wide. Traditional SfM is prone to error accumulation when incrementally adding images, resulting in overall drift or deformation of the three-dimensional point cloud.

[0005] (3) Sparse point clouds are insufficient to express complex terrain: Conventional SfM sparse point clouds have limited density and are difficult to express the shape and steep texture of the mountain surface in detail.

[0006] (4) Lack of depth perception of complex terrain: Traditional MVS often fails in mountainous scenes with long distances, high elevation differences and severe obstruction, and is prone to large-area depth loss.

[0007] (5) Incremental processing is inefficient: When faced with thousands of aerial images, traditional global BA (Bundle Adjustment) is extremely time-consuming and difficult to use in real-time or near real-time in engineering applications. Summary of the Invention

[0008] The purpose of this invention is to overcome the shortcomings of the prior art and provide an incremental 3D reconstruction method and system for UAVs in complex mountainous environments. By constructing a scale-space adjustment factor for terrain perception and adaptively adjusting the scale parameters based on the slope angle and average curvature, the problem of unstable image matching in complex mountainous environments is solved.

[0009] The objective of this invention is achieved through the following technical solution:

[0010] The first aspect of this invention discloses an incremental 3D reconstruction method for unmanned aerial vehicles (UAVs) in complex mountainous environments, comprising:

[0011] Control the drone to fly along a preset route and collect multi-view monocular images and corresponding attitude and position information;

[0012] Extract feature points from the images and match them to establish feature correspondences between images;

[0013] The camera pose is estimated based on the pose and position information corresponding to the first frame image; for each new image, the new image and its co-view associated images are jointly optimized by local bundle adjustment to update the camera pose and 3D point coordinates; a sparse 3D point cloud and camera pose sequence are constructed based on the updated camera pose and 3D point coordinates.

[0014] Each 3D point is initialized as a 3D Gaussian primitive based on a sparse 3D point cloud. The 3D Gaussian primitive includes Gaussian center position parameters, rotation parameters, scaling parameters, color spherical harmonic function coefficients, and opacity parameters.

[0015] Calculate the rate of change of the normal in the neighborhood of each 3D point in the sparse 3D point cloud; if it is greater than the preset steepness threshold, adjust the scaling matrix to make the Gaussian ellipsoid lengthen along the slope normal and shorten along the tangent, and align the principal axis of the rotation matrix with the slope direction to form an anisotropic distribution that is lengthened along the slope; otherwise, set the scaling matrix to a multiple of the identity matrix to make the Gaussian ellipsoid be isotropic.

[0016] The three-dimensional Gaussian primitives are projected onto the two-dimensional image plane using the view transformation matrix to obtain the two-dimensional center point. The three-dimensional Gaussian covariance matrix is ​​projected onto the Jacobian matrix to obtain the two-dimensional covariance. The two constitute the two-dimensional Gaussian distribution primitives. The two-dimensional Gaussian distribution primitives are weighted and accumulated based on opacity to render and generate a two-dimensional image.

[0017] A joint loss function for the two-dimensional image and the real image is constructed, and the parameters of the three-dimensional Gaussian unit are optimized in reverse until convergence, thus obtaining a static three-dimensional Gaussian model.

[0018] Incremental 3D reconstruction results are generated based on the static 3D Gaussian model.

[0019] Furthermore, feature points are extracted from the images and matched to establish feature correspondences between images, including:

[0020] The slope angle and average curvature of the terrain surface are calculated. The slope angle is obtained by calculating the first partial derivative of the elevation function, and the average curvature is obtained by calculating the second partial derivative of the elevation function.

[0021] The adjusted scale parameter is determined based on the slope angle and the average curvature. The adjusted scale parameter increases when the slope angle increases and when the absolute value of the average curvature increases.

[0022] The adjusted scale parameters are substituted into the Gaussian kernel function to generate an adaptive scale space. Local extreme points are detected in the adaptive scale space to obtain candidate feature points.

[0023] Candidate feature points extracted from different images are matched across images to obtain initial matching pairs;

[0024] Based on epipolar geometry constraints, erroneous matches in the initial matching pairs are eliminated through essential matrix verification, and feature correspondences between images are established, which contain multiple sets of matching point pairs.

[0025] Furthermore, feature points are extracted from the images and matched to establish feature correspondences between images, including:

[0026] A semantic segmentation network is used to perform pixel-level terrain classification on the image, generating semantic masks for vegetation areas, rock areas, and bare soil areas.

[0027] The slope angle and average curvature of the terrain surface are calculated. The slope angle is obtained by calculating the first partial derivative of the elevation function, and the average curvature is obtained by calculating the second partial derivative of the elevation function.

[0028] The foundation dimensional parameters are determined based on the slope angle and the average curvature; the foundation dimensional parameters increase as the slope angle increases and as the absolute value of the average curvature increases.

[0029] Adjusted scale parameters are generated based on the basic scale parameters: for vegetated areas, the basic scale parameters are weighted by superimposing the normalized vegetation index and used as the adjusted scale parameters; for rocky areas, the basic scale parameters are weighted by increasing the curvature weight coefficient and used as the adjusted scale parameters; for bare soil areas, the basic scale parameters are directly used as the adjusted scale parameters.

[0030] The adjusted scale parameters are substituted into the Gaussian kernel function to generate an adaptive scale space, and a regional adaptive Gaussian scale space is generated according to the semantic mask.

[0031] Local extreme points are detected in the regional adaptive Gaussian scale space, and the detected extreme points are used as candidate feature points: For vegetated areas, the scale space image is weighted based on the normalized vegetation index (NVI) to detect local extreme points, and extreme points with the NVI below a preset threshold are filtered out; For rocky areas, the scale space image is stretched in grayscale to detect local extreme points, and extreme points located in areas spanning slopes are removed; For bare soil areas, local extreme points are detected directly.

[0032] Based on the semantic mask, feature descriptors are generated for the candidate feature points: for vegetated areas, SIFT descriptors are generated based on the normalized vegetation index and optical flow information between adjacent frames is fused; for rocky areas, SIFT-AID descriptors are generated; for bare soil areas, standard SIFT descriptors are generated.

[0033] Candidate feature points extracted from different images are matched across images to obtain initial matching pairs;

[0034] The initial matching pairs are filtered according to the semantic mask: for initial matching pairs in vegetated areas, initial matching pairs that do not meet the dynamic optical flow matching constraints are removed; for initial matching pairs in rock areas, initial matching pairs that do not meet the local normal consistency test are removed; for initial matching pairs in bare soil areas, all initial matching pairs are retained.

[0035] Based on epipolar geometry constraints, the distance from each matching point pair to the epipolar line is calculated using the essential matrix. According to the semantic mask, for matching pairs in vegetated areas, pairs are discarded if the distance is greater than a first preset threshold; for matching pairs in rocky areas, pairs are discarded if the distance is greater than a second preset threshold; for matching pairs in bare soil areas, pairs are discarded if the distance is greater than a third preset threshold. The retained matching pairs are used as feature correspondences between images, and these feature correspondences contain multiple sets of matching point pairs. The second preset threshold is greater than the first preset threshold, and the first preset threshold is greater than the third preset threshold.

[0036] Furthermore, for each new image, local bundle adjustment is used to jointly optimize the new image and its associated co-view images, updating the camera pose and 3D point coordinates, including:

[0037] For each new image, the camera pose of the new image is solved using the PnP algorithm based on the matching point pairs between the new image and the reconstructed image.

[0038] Based on the obtained camera pose, the unreconstructed matching point pairs between the current image and the reconstructed image are triangulated to generate the coordinates of the newly added 3D points.

[0039] The newly added 3D point coordinates are merged with the existing 3D point coordinates to update the 3D point cloud;

[0040] Local bundle adjustment is performed, and the newly added image and the preset number of associated images with the highest co-view degree are used to form a local window. The camera pose and 3D point coordinates of all images in the local window are optimized by minimizing the reprojection error.

[0041] Furthermore, the rate of change of the normal vector within the neighborhood of each 3D point in the sparse 3D point cloud is calculated, including:

[0042] For each three-dimensional point in the sparse three-dimensional point cloud, the spatial neighborhood of the three-dimensional point is formed by selecting the adjacent three-dimensional points within a preset radius.

[0043] By performing plane fitting on all three-dimensional points in the spatial neighborhood using principal component analysis or least squares method, a local plane and its normal vector are obtained.

[0044] The rate of change of the normal vector of a three-dimensional point is calculated based on the difference between the normal vector and the normal vectors of each adjacent three-dimensional point in the spatial neighborhood.

[0045] Furthermore, a two-dimensional Gaussian primitive is projected onto a two-dimensional image plane using a view transformation matrix to obtain a two-dimensional center point. A two-dimensional covariance matrix is ​​then projected onto the three-dimensional Gaussian covariance matrix using a Jacobian matrix. These two components constitute a two-dimensional Gaussian distribution primitive. An opacity-based weighted accumulation is then performed on the two-dimensional Gaussian distribution primitive to render and generate a two-dimensional image, including:

[0046] For each 3D Gaussian primitive, multiply its center position coordinates by the view transformation matrix to obtain the projected coordinates of the 3D Gaussian primitive on the 2D image plane, which serve as the 2D center point.

[0047] The three-dimensional Gaussian covariance matrix of the three-dimensional Gaussian element is multiplied sequentially with the view transformation matrix and the Jacobian matrix to obtain the two-dimensional covariance after projection.

[0048] Based on the two-dimensional center point and the two-dimensional covariance, construct the two-dimensional Gaussian distribution primitive of the Gaussian primitive on the two-dimensional image plane;

[0049] For each pixel on the two-dimensional image plane, determine all two-dimensional Gaussian distribution primitives covering that pixel and sort them in ascending order of depth value;

[0050] Following the sorted order, the color value and opacity of each two-dimensional Gaussian distribution primitive are successively fused to obtain the final color value of the pixel, thus generating a two-dimensional image.

[0051] Furthermore, a joint loss function for the 2D image and the real image is constructed, and the parameters of the 3D Gaussian primitives are optimized inversely until convergence, including:

[0052] Construct a joint loss function that includes the L1 loss function and the structural similarity loss function;

[0053] An adaptive moment estimation optimizer is used to backpropagate the gradient of the loss function to each parameter of the Gaussian unit through a differentiable rendering pipeline, and iteratively update the parameters until the loss function converges.

[0054] Furthermore, incremental 3D reconstruction results are generated based on the static 3D Gaussian model, including:

[0055] A dynamic Gaussian mapping network is constructed. Based on the static three-dimensional Gaussian model, the four-dimensional spatiotemporal coordinates are decomposed and mapped into six two-dimensional feature planes. The six two-dimensional feature planes include three pure spatial planes and three spatiotemporal planes.

[0056] A lightweight multilayer perceptron decoder is used to extract and fuse features from each two-dimensional feature plane, and output the displacement, rotation, and scaling offsets of each three-dimensional Gaussian element in the time dimension.

[0057] The attention mechanism of the dynamic Gaussian deformation network is guided by terrain features, which include at least one of slope, aspect, and vegetation cover. The vegetation cover is obtained by semantic segmentation of the color and texture features of aerial images. The loss weight guided by terrain semantics makes the dynamic deformation modeling focus on dense vegetation areas and cliff edge areas.

[0058] Dynamic temporal smoothing regularization constraints and terrain consistency constraints are applied. Temporal smoothing is achieved by applying Laplace filtering to the eigenvalues ​​of the spatiotemporal plane, and terrain consistency constraints are achieved through regularization terms. In vegetated areas, the magnitudes of displacement offset, rotation offset, and scaling offset are constrained in relation to the terrain slope and vegetation type, and constraints are applied using a pre-statistical vegetation sway model. In rocky areas, displacement offset, rotation offset, and scaling offset are forced to be zero to maintain static conditions. In terrain transition areas, bilinear interpolation is used to smoothly transition deformation parameters between dynamic and static regions, constructing a dynamic three-dimensional Gaussian model.

[0059] A second aspect of the present invention discloses an incremental 3D reconstruction system for unmanned aerial vehicles (UAVs) in complex mountainous environments, comprising:

[0060] The image acquisition module is used to simultaneously acquire multi-view monocular images and corresponding attitude and position information when the UAV is flying along a preset route;

[0061] The feature matching module is used to extract feature points from images and perform matching to establish feature correspondences between images.

[0062] The 3D point cloud construction module is used to estimate the camera pose based on the pose and position information corresponding to the first frame image. For each new image, the new image and its co-view associated images are jointly optimized by local bundle adjustment to update the camera pose and 3D point coordinates. Based on the updated camera pose and 3D point coordinates, a sparse 3D point cloud and camera pose sequence are constructed.

[0063] The 3D Gaussian primitive generation module is used to initialize each 3D point as a 3D Gaussian primitive based on a sparse 3D point cloud. The 3D Gaussian primitive includes Gaussian center position parameters, rotation parameters, scaling parameters, color spherical harmonic function coefficients, and opacity parameters.

[0064] The normal change rate calculation module is used to calculate the normal change rate in the neighborhood of each 3D point in a sparse 3D point cloud. If it is greater than the preset steepness threshold, the scaling matrix is ​​adjusted to make the Gaussian ellipsoid lengthen along the slope normal and shorten along the tangent, and the principal axis of the rotation matrix is ​​aligned with the slope direction to form an anisotropic distribution that is lengthened along the slope. Otherwise, the scaling matrix is ​​set to a multiple of the identity matrix to make the Gaussian ellipsoid be isotropic.

[0065] The 2D image generation module is used to project the 3D Gaussian primitives onto the 2D image plane through the view transformation matrix to obtain the 2D center point, and to project the 3D Gaussian covariance matrix using the Jacobian matrix to obtain the 2D covariance. The two constitute the 2D Gaussian distribution primitives. The 2D Gaussian distribution primitives are weighted and accumulated based on opacity to render and generate a 2D image.

[0066] The static model generation module is used to construct a joint loss function between the two-dimensional image and the real image, and to back-optimize the parameters of the three-dimensional Gaussian unit until convergence, thus obtaining a static three-dimensional Gaussian model.

[0067] The 3D reconstruction module is used to generate incremental 3D reconstruction results based on a static 3D Gaussian model.

[0068] The beneficial effects of this invention are:

[0069] (1) This invention can complete the three-dimensional reconstruction of complex mountainous environments by relying only on visible light image data captured by UAVs. It does not require expensive hardware equipment such as lidar and structured light, which reduces the system hardware cost and platform load requirements. It is suitable for small and medium-sized UAV platforms and has stronger engineering practicality and promotion value.

[0070] (2) This invention constructs a scale space adjustment factor for terrain perception, and adaptively adjusts the Gaussian blur scale parameter based on the slope angle and average curvature of the terrain surface. In areas with steep slopes and large curvature, the scale parameter is increased to retain more feature points, and in flat areas, the scale parameter is decreased to reduce feature redundancy. At the same time, combined with the illumination invariant descriptor, it effectively enhances the robustness of feature detection and matching in complex mountainous environments such as weak texture vegetation areas, shadowed areas, and strong light rock areas, and greatly reduces the matching failure rate.

[0071] (3) In the incremental camera pose solution process, the present invention introduces a joint optimization mechanism of PnP pose estimation and local bundle adjustment (BA). By forming a local window with the current newly added image and the preset number of associated images with the highest co-view degree, the camera pose and three-dimensional point coordinates are optimized at the same time within the window. This can effectively suppress the drift problem caused by the gradual accumulation of errors during long-haul flight and improve the geometric accuracy of reconstruction.

[0072] (4) This invention uses a terrain-guided Gaussian initialization strategy to dynamically adjust the anisotropic distribution of the Gaussian ellipsoid according to the rate of change of the normal in the neighborhood of each 3D point in the sparse 3D point cloud: in steep terrain areas, the Gaussian ellipsoid is stretched along the slope direction, while in flat areas it is kept isotropic, thus more accurately approximating the complex terrain surface; combined with differentiable rendering and joint loss function inverse optimization, it effectively fills the large-area depth loss and structural void problems that are easy to occur in traditional SfM / MVS in areas such as vegetation occlusion and steep cliffs, and improves the geometric integrity of the reconstruction model;

[0073] (5) The present invention adopts a splash-based three-dimensional Gaussian differentiable rendering method, which projects three-dimensional Gaussian primitives onto a two-dimensional image plane through view transformation, and uses the alpha mixing formula to perform pixel-by-pixel fusion rendering, avoiding a large number of light sampling calculations in the neural radiation field, and has a faster convergence speed and higher rendering efficiency, which is particularly suitable for rapid modeling tasks in large-scale mountainous areas.

[0074] (6) This invention models the color of three-dimensional Gaussian primitives using spherical harmonic functions and combines alpha blending rendering with opacity parameters, so that the reconstructed model can present realistic and continuous lighting effects from different perspectives, with higher visual fidelity, thereby accurately reproducing the texture and light and shadow changes of complex mountain surfaces.

[0075] (7) This invention introduces a dynamic Gaussian mapping network on the basis of static three-dimensional Gaussian modeling, decomposes the four-dimensional spatiotemporal coordinates into six two-dimensional feature planes (three pure spatial planes and three spatiotemporal planes), and outputs the displacement, rotation and scaling offsets of each Gaussian element in the time dimension through a lightweight multilayer perceptron decoder; at the same time, it uses terrain features (slope, aspect and vegetation coverage) to guide the attention mechanism of the dynamic network, and constrains the deformation amplitude by using dynamic temporal smoothing regularization and terrain consistency constraints, associating the vegetation swing model in the vegetation area to constrain the deformation amplitude, and forcing the deformation parameter to be zero in the rock area, thereby achieving effective decoupling between the dynamic and static areas, making the dynamic modeling conform to the physical laws of mountains, and improving the physical interpretability of dynamic deformation;

[0076] (8) This invention forms a complete closed-loop system from data acquisition, incremental reconstruction, terrain-guided initialization, differentiable rendering optimization to dynamic modeling. It has good automated processing capabilities and can be widely applied to complex mountainous scenarios such as natural resource surveys, disaster monitoring, landslide assessment, power line inspection, and emergency mapping. It can also be further extended to application scenarios such as multi-UAV collaborative reconstruction, real-time updates and digital twins, and has significant engineering practical value. Attached Figure Description

[0077] Figure 1 This is a flowchart of an incremental 3D reconstruction method for UAVs in this invention;

[0078] Figure 2 This is a schematic diagram of the Gauss Pyramid.

[0079] Figure 3 This is a schematic diagram of the spatial distribution of spherical harmonic functions of different orders.

[0080] Figure 4 This is a schematic diagram of a 3D sand table of point clouds;

[0081] Figure 5 This is a schematic diagram of the 3D sand table after sharpening and smoothing.

[0082] Figure 6 This is a block diagram of one component of the UAV incremental 3D reconstruction system in this invention. Detailed Implementation

[0083] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. 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.

[0084] See Figures 1-6 This embodiment provides an incremental 3D reconstruction method and system for UAVs in complex mountainous environments:

[0085] The first aspect of this embodiment discloses an incremental 3D reconstruction method for unmanned aerial vehicles (UAVs) in complex mountainous environments, such as... Figure 1 As shown, the UAV incremental 3D reconstruction method includes steps S100 to S800.

[0086] Step S100. Control the UAV to fly along the preset route and collect multi-view monocular images and corresponding attitude and position information.

[0087] In this embodiment, adaptive route planning is performed based on the terrain undulation of the mountainous area to be reconstructed. The preset route must meet at least the following conditions: the forward overlap rate is not less than 60%, preferably 70% to 85%; the lateral overlap rate is not less than 50%, preferably 60% to 75%; in areas with large slopes or severe terrain undulations, the route density is increased or multi-altitude layered routes are adopted; the route view coverage includes a top view and at least one oblique view (such as 15° to 45°) to enhance the three-dimensional structure observation capability.

[0088] The multi-view monocular image sequence must meet at least the following conditions: the number of images is not less than 100, preferably not less than 500; any surface area is covered by at least 3 or more images from different perspectives; there is a parallax variation of not less than a preset threshold between adjacent images to ensure triangulation stability; the temporal continuity is maintained during image acquisition, and the time interval between adjacent frames is less than a preset threshold (e.g., 1 second) to reduce the impact of dynamic disturbances.

[0089] In areas with steep terrain, by adjusting the flight attitude or introducing side-view compensation routes, multi-directional cross-view coverage can be created in local areas to improve the observability of cliffs or obstructed areas.

[0090] The UAV attitude information is acquired by an inertial measurement unit (IMU), and the position information is acquired by a positioning module (GPS / RTK), forming a time-synchronized image-pose data pair.

[0091] Step S200. Extract feature points from the image and perform matching to establish feature correspondence between images.

[0092] In some embodiments of this example, step S200 specifically includes steps S210 to S250.

[0093] Step S210. Calculate the slope angle and average curvature of the terrain surface. The slope angle is calculated using the first partial derivative of the elevation function, and the average curvature is calculated using the second partial derivative of the elevation function.

[0094] Partial derivatives can be obtained through convolution using the difference of Gaussians operator or the Sobel operator to ensure computational accuracy and efficiency.

[0095] The slope angle of the terrain surface It reflects the slope of the terrain surface, and is derived from the modulus of the first partial derivative and the rate of change of elevation. Its calculation formula is as follows:

[0096]

[0097] In the formula, and The first partial derivative of the elevation function; Represents the terrain surface in pixels The greater the rate of change of elevation at a given location, the greater the slope angle. The larger the value, the steeper the terrain, consistent with the design logic of subsequent scale adjustment factors.

[0098] The average curvature of the terrain surface It reflects the curvature of the terrain surface, combining the average effect of the two principal curvatures. It is simple to calculate and effectively characterizes the undulations of the terrain. Its calculation formula is:

[0099]

[0100] In the formula, and The first partial derivative of the elevation function; , and The second partial derivative of the elevation function is given. In the formula for calculating the average curvature, the numerator reflects the combined effect of the second derivative of the terrain surface, while the denominator is used for normalization to avoid the influence of the elevation change rate on the curvature calculation result, ensuring that the curvature value can stably represent the degree of curvature of the terrain. The larger the absolute value, the more severe the terrain undulation.

[0101] Step S220. Determine the adjusted scale parameter based on the slope angle and the average curvature. The adjusted scale parameter increases when the slope angle increases and when the absolute value of the average curvature increases.

[0102] Based on the average curvature of the terrain surface, a curvature-driven scale adjustment factor is constructed. The calculation formula for the curvature-driven scale adjustment factor is as follows:

[0103]

[0104] In the formula, A curvature-driven scaling factor; It is the absolute value of the average curvature, which can avoid the influence of positive and negative curvature on the direction of scale adjustment; This is the curvature weighting coefficient, used to balance the influence of curvature on scale adjustment. Its value ranges from 0.1 to 0.5 and can be adaptively adjusted according to the degree of undulation of specific terrain data.

[0105] Based on the slope angle of the terrain surface, a slope-driven scale adjustment factor is constructed. The calculation formula for the slope-driven scale adjustment factor is as follows:

[0106]

[0107] In the formula, A slope-driven scaling factor; It reflects the steepness of the slope; the larger the slope angle, the larger its tangent value, and the more significant the scale adjustment range. The slope weighting coefficient has a value range of 0.05 to 0.2. It is used to control the sensitivity of the slope to scale adjustment and to avoid abnormal scale increase due to excessive slope.

[0108] The formula for calculating the adjusted scale parameter is as follows:

[0109]

[0110] In the formula, The adjusted scale parameters; It is a preset base scale, serving as the benchmark scale for flat terrain areas.

[0111] Step S230. Substitute the adjusted scale parameters into the Gaussian kernel function to generate an adaptive scale space, and detect local extreme points in the adaptive scale space to obtain candidate feature points.

[0112] The formula for calculating the scale-space image with adjustment factor is:

[0113]

[0114] In the formula, The original terrain image (or terrain point cloud projection image); Represents a two-dimensional convolution operation; For A two-dimensional Gaussian kernel function with standard deviation is used to achieve image smoothing at different scales, thereby generating an adaptive scale space that fits the terrain features, providing a reliable scale basis for subsequent tasks such as terrain feature extraction and matching.

[0115] The process of generating the adaptive scale space is as follows: the adjusted scale parameters are introduced into the Gaussian kernel function, and the adaptive scale space is generated through Gaussian convolution.

[0116] In some embodiments of this example, the feature point extraction uses the Scale Invariant Feature Transform (SIFT) algorithm, specifically including steps S231 to S236.

[0117] Step S231. Constructing the scale space of the image: Apply Gaussian kernel functions with different blur scale factors to the original image for convolution operations to generate scale space images. These scale space images form a Gaussian pyramid, where each scale space image represents an image version at a different scale, such as... Figure 2 As shown.

[0118] The Gaussian blur kernel function is:

[0119]

[0120] In the formula, For fuzzy scaling factor, For pixel coordinates, is the Gaussian blur kernel function at a given scale.

[0121] Step S232. Construct Gaussian difference scale space: Based on the difference in Gaussian (DoG) function, differ images of adjacent scales in the Gaussian pyramid to obtain a difference in Gaussian image.

[0122] The difference-of-Gaussian (DoG) function is:

[0123]

[0124] In the formula, Indicated in scale The difference-of-Gaussian image below; Represents the pixel coordinates in the image; Indicates the current scale parameter; is the scaling factor between adjacent scale layers, and ; This indicates that the original image has a standard deviation of 1. The scale-space image obtained after Gaussian kernel convolution; This represents the scale-space image corresponding to an adjacent higher-scale layer.

[0125] Scale space image It can be represented as:

[0126]

[0127] In the formula, The standard deviation is The two-dimensional Gaussian kernel function, For the original image, This represents the convolution operation.

[0128] Step S233. Detect scale space extreme points: In the Gaussian difference scale space, each pixel is compared with its 8 neighboring pixels at the same scale and 9 pixels at each of the adjacent scales above and below. If the pixel is the maximum or minimum value, it is taken as a candidate feature point.

[0129] Step S234. Precisely locate feature points: refine the candidate feature points at the sub-pixel level by Taylor series expansion, and remove feature points with contrast below the preset threshold and unstable edge response.

[0130] Step S235. Assign a primary direction to the feature point: Using the feature point as the center, select all pixels within the neighborhood window and calculate the gradient magnitude and direction for each pixel. Construct a gradient direction histogram within the region surrounding the feature point, dividing the gradient directions into 36 intervals, each 10 degrees. During the statistical analysis, Gaussian weighting is applied based on the pixel's distance from the feature point. The direction with the largest peak in the histogram is selected as the primary direction for that feature point; if the magnitude of other directions reaches more than 80% of the magnitude of the primary peak, then that direction is designated as an auxiliary direction, meaning the feature point is assigned multiple directions.

[0131] Step S236. Generate feature descriptor: Rotate the neighborhood of the feature point to the main direction to eliminate the effect of rotation. Take a 16×16 pixel window in the rotated neighborhood and divide it into 4×4 sub-blocks. Calculate the gradient histogram of 8 directions in each sub-block to form a 128-dimensional feature vector. Normalize this vector to obtain the final feature descriptor.

[0132] Step S240. Perform cross-image matching on the candidate feature points extracted from different images to obtain initial matching pairs.

[0133] In this embodiment, for feature point matching between adjacent images, the nearest neighbor and second nearest neighbor distance ratio criterion is adopted. For each feature point, the Euclidean distance between its feature descriptor and all feature point descriptors in the other image is calculated, and the minimum distance d1 and the second minimum distance d2 are taken to calculate the distance ratio. If the NNDR is less than a preset threshold (e.g., 0.7), then the two feature points are considered to be successfully matched, and an initial matching pair is obtained.

[0134] Step S250. Based on epipolar geometry constraints, erroneous matches in the initial matching pairs are eliminated by verifying the essential matrix, and feature correspondences between images are established, wherein the feature correspondences contain multiple sets of matching point pairs.

[0135] The essential matrix is:

[0136]

[0137] In the formula, E is the essential matrix; x1 and x2 are the matching points in the two perspectives, respectively.

[0138] In this embodiment, the random sampling consensus algorithm is used to iteratively estimate the essential matrix. The distance from the matching point pair to the epipolar line is calculated based on the epipolar geometric constraints. Matching point pairs with a distance greater than a preset threshold are identified as outliers and removed, thereby obtaining a reliable feature correspondence between images.

[0139] In some embodiments of this example, step S200 specifically includes steps K201 to K210.

[0140] Step K201. Use a semantic segmentation network to perform pixel-level terrain classification on the image and generate semantic masks for vegetation areas, rock areas and bare soil areas.

[0141] Step K202. Calculate the slope angle and average curvature of the terrain surface. The slope angle is calculated using the first partial derivative of the elevation function, and the average curvature is calculated using the second partial derivative of the elevation function.

[0142] Step K203. Determine the basic scale parameters based on the slope angle and the average curvature; the basic scale parameters increase as the slope angle increases, and the basic scale parameters increase as the absolute value of the average curvature increases.

[0143] Step K204. Generate adjusted scale parameters based on the basic scale parameters: For vegetated areas, the basic scale parameters are weighted by superimposing the normalized vegetation index and used as the adjusted scale parameters; for rocky areas, the basic scale parameters are weighted by increasing the curvature weight coefficient and used as the adjusted scale parameters; for bare soil areas, the basic scale parameters are directly used as the adjusted scale parameters.

[0144] Specifically, in vegetated areas, the basic scale parameters are combined with the Normalized Difference Vegetation Index (NDVI) using a product-stacking method to obtain scale parameters applicable to vegetated areas. The weighting formula is as follows:

[0145]

[0146] in, These are the scale parameters after vegetation area adjustment; Basic scale parameters; This is a vegetation index weighting coefficient, ranging from 0.1 to 0.5, used to control the intensity of the NDVI's influence on scale. is the normalized vegetation index at pixel (x,y), normalized to the [0,1] interval.

[0147] In some embodiments of this example, the curvature weighting coefficient is increased when any of the following conditions are met. The current region is a rocky area; the local rate of change of curvature exceeds a preset threshold; and the local elevation gradient changes drastically. The curvature weighting coefficient can be expressed as:

[0148]

[0149] in, For position Adaptive curvature weighting coefficient at the location; Basic weighting coefficients; Intensity of curvature change (such as curvature gradient or variance); This is a region indicator function, which takes the value 1 when the region is rocky and 0 otherwise. , This is the adjustment coefficient.

[0150] In these implementations, differentiated scale adjustment strategies are applied to the three types of terrain regions based on semantic masks: For vegetated regions, a normalized vegetation index weighting is superimposed on the dual-drive scale adjustment, appropriately reducing the base scale in densely vegetated areas to enhance the detection capability of fine-scale feature points and effectively suppress interference caused by repetitive vegetation textures; for rocky regions, the curvature weight coefficient is increased to enlarge the Gaussian blur scale, suppressing noise in rocky regions while preserving the weak texture features of rock edges; for bare soil regions, the original dual-drive scale adjustment formula is directly used with default scale parameters to ensure the stability of feature extraction. Through the above regional scale adjustments, a regional adaptive Gaussian scale space is generated, completing the optimization and enhancement of the original "constructing an adaptive scale space" step.

[0151] Specifically, in this embodiment, the calculation method for the Normalized Difference Vegetation Index (NDVI) is as follows: A lightweight UAV multispectral camera adapted for mountain aerial photography scenarios is used. This camera must support simultaneous imaging in the red band (wavelength range 620–670 nm) and the near-infrared band (wavelength range 800–850 nm), and have a timestamp synchronization function with the UAV's main camera (used to acquire RGB images required for feature extraction). This avoids pixel position shifts caused by acquisition time differences, ensuring accurate correspondence between NDVI calculation and the semantic mask and feature extraction positions of the vegetation area. When performing mountain aerial photography tasks, the UAV, equipped with a multispectral camera, flies synchronously with the main camera, acquiring image data according to a preset aerial flight path (flight altitude 80–120 m, meeting the needs of mountain terrain reconstruction). The multispectral camera simultaneously acquires red band and near-infrared band images, while the main camera acquires RGB aerial images. During the acquisition process, the UAV's GPS position information and attitude parameters (heading angle, pitch angle, roll angle) are recorded simultaneously for subsequent image registration and correction. The acquired red band and near-infrared band images are preprocessed to remove noise and interference during the acquisition process and ensure the accuracy of the reflectance values. The specific preprocessing steps include: (1) eliminating the radiation deviation caused by the difference in response of the multispectral camera and the change in light intensity, converting the original image gray value into the actual reflectance value, and ensuring that the band reflectance values ​​acquired in different areas and at different times are comparable; (2) combining the UAV GPS position information and attitude parameters, performing geometric distortion correction on the band images, eliminating the image offset caused by flight attitude fluctuations and terrain undulations, so that the red band image, near-infrared band image and the RGB image and semantic mask acquired by the main camera can be accurately aligned at the pixel level; (3) using the feature matching algorithm, registering the preprocessed red band image, near-infrared band image and the corresponding RGB image, ensuring that the red reflectance value (Red) and near-infrared reflectance value (NIR) corresponding to the same pixel position are accurately matched, providing a basis for NDVI pixel-by-pixel calculation. After preprocessing, the red reflectance value corresponding to each pixel (x,y) is extracted from the red band image and denoted as Red(x,y); the near-infrared reflectance value corresponding to the same pixel (x,y) is extracted from the near-infrared band image and denoted as NIR(x,y). The extracted reflectance values ​​for both bands are normalized to the [0,1] interval for subsequent NDVI calculation. The formula for calculating NDVI is:

[0152]

[0153] In the formula, NDVI represents the normalized vegetation index, NIR represents the near-infrared reflectance value, and Red represents the red light reflectance value.

[0154] Step K205. Substitute the adjusted scale parameters into the Gaussian kernel function to generate an adaptive scale space, and generate a regional adaptive Gaussian scale space according to the semantic mask.

[0155] Step K206. Detect local extreme points in the regional adaptive Gaussian scale space and use the detected extreme points as candidate feature points: For vegetated areas, the scale space image is weighted based on the normalized vegetation index to detect local extreme points, and extreme points with a normalized vegetation index lower than a preset threshold are filtered out; For rocky areas, the scale space image is stretched in grayscale to detect local extreme points, and extreme points located in the slope crossing area are removed; For bare soil areas, local extreme points are detected directly.

[0156] Specifically, the calculation formula for the grayscale value stretching is as follows:

[0157]

[0158] In the formula, For the original image in pixels The grayscale value at that location; The grayscale value after stretching; , These are the minimum and maximum gray values ​​in the current image or local region; The target grayscale range (usually 0 and 255).

[0159] Specifically, the method for determining extreme points located in the slope-crossing region is as follows: for each candidate extreme point, calculate its neighborhood local surface normal vector n from the currently reconstructed sparse point cloud or depth map; for the matched feature point pair (p1, p2), obtain its corresponding local normal vector. , And calculate the angle between the normal vectors. If the included angle If the angle is greater than a preset angle threshold (30° in this embodiment), the feature point is determined to be located on different slopes and belongs to a cross-slope extreme point, and is therefore removed; if the included angle is greater than the .... If the value is less than or equal to the threshold, it is considered a valid feature point on the same slope and is retained.

[0160] In these implementations, regional adaptive Gaussian scale space is used, combined with semantic masks, to achieve regional extreme point detection and filtering: For vegetated areas, vegetation texture features are enhanced by combining vegetation indices (such as NDVI), and extreme points in the vegetation canopy area are enriched and completed. Invalid extreme points with low vegetation response are removed by using an NDVI threshold, solving the problem of extreme point sparsity caused by repetitive vegetation textures; For rocky areas, DoG image grayscale enhancement is performed on weak texture areas to increase the number of extreme point detections. At the same time, false extreme points across slopes are removed by local normal consistency pre-screening to avoid interference from invalid feature points in subsequent matching; For bare soil areas, a standard extreme point detection and precise localization strategy is adopted, using the Taylor expansion precise localization method in step S200 to ensure the accuracy of feature point detection. Finally, a regional accurate feature point set is generated, completing the optimization and enhancement of the original "detection of local extreme points" step.

[0161] Step K207. Generate feature descriptors for the candidate feature points based on the semantic mask: for vegetated areas, generate SIFT descriptors based on the normalized vegetation index and fuse optical flow information between adjacent frames; for rock areas, generate SIFT-AID descriptors; for bare soil areas, generate standard SIFT descriptors.

[0162] Specifically, the SIFT-AID descriptor generation method is as follows: Image feature points are detected using the SIFT algorithm to determine their position, scale, and main direction. A fixed-size image block is cropped centered on the feature point according to the scale and main direction, and this image block is rotated and scaled to obtain an affine-normalized image block. Subsequently, the normalized image block is input into a preset affine-invariant feature network, which outputs a high-dimensional real-valued feature vector through multi-layer convolution and pooling operations. Each element of this high-dimensional real-valued feature vector is symbolized, with components greater than zero recorded as 1 and components not greater than zero recorded as 0, thus obtaining a binarized SIFT-AID descriptor. During feature matching, the Hamming distance is obtained by calculating the number of bits differing at the same position between two SIFT-AID descriptors, and the similarity of the matches is measured by the Hamming distance.

[0163] In these implementations, during the feature descriptor generation stage, differentiated feature descriptor generation strategies are adopted based on the region categories defined by the semantic mask: For vegetated regions, the SIFT descriptor generation process is performed based on the NDVI-enhanced image, and optical flow information between adjacent frames is superimposed for weighted processing to add temporal displacement features to feature points, generating vegetation-enhanced SIFT descriptors and solving the descriptor instability problem caused by dynamic swaying of vegetation; For rocky regions, the original SIFT descriptor is replaced with a SIFT-AID illumination-invariant descriptor, and the influence of sudden illumination changes on the descriptor is eliminated by converting the grayscale of the feature point's neighborhood into a relative gradient difference, while the descriptor structure is optimized by combining local normal information to improve the robustness of feature description in rocky regions; For bare soil regions, the standard SIFT feature descriptor generation method is directly adopted without additional modification, adapting to the characteristics of uniform illumination and weak texture in bare soil regions. Through the above differentiated processing, customized feature descriptors for different regions are generated, completing the optimization and enhancement of the original "feature descriptor generation" stage.

[0164] The method for calculating the local normal information is as follows: taking the three-dimensional point corresponding to the feature point as the center, select the neighborhood three-dimensional point set to construct the covariance matrix, perform eigenvalue decomposition on it, take the eigenvector corresponding to the smallest eigenvalue as the local normal vector, and normalize it to the unit normal vector. At the same time, unify the normal vector to point to the viewpoint direction, and finally obtain the local normal information of the feature point.

[0165] Step K208. Perform cross-image matching on the candidate feature points extracted from different images to obtain initial matching pairs.

[0166] Step K209. Filter the initial matching pairs according to the semantic mask: For the initial matching pairs in the vegetation area, remove the initial matching pairs that do not meet the dynamic optical flow matching constraints; for the initial matching pairs in the rock area, remove the initial matching pairs that do not meet the local normal consistency test; for the initial matching pairs in the bare soil area, retain all the initial matching pairs.

[0167] Specifically, the method for determining whether the initial matching pair satisfies the dynamic optical flow matching constraint is as follows: the dense optical flow field is calculated from the aerial images of adjacent frames, and the optical flow displacement vector of the current feature point in the reference frame is obtained. in These represent the predicted optical flow displacements of pixels in the horizontal and vertical directions, respectively. Based on the initial matching pairs, the actual matching displacement vectors of feature points between the current frame and the reference frame are obtained. Calculate the angle between the direction of the optical flow predicted displacement and the actual matching displacement. Calculate the relative deviation of displacement amplitude. The initial matching pair is considered to satisfy the optical flow dynamic matching constraint and is a valid match only if both of the following conditions are met simultaneously; otherwise, it is considered a false match and is discarded: the directional angle is not greater than a preset angle threshold, which is 15° in this embodiment; and the relative deviation of the displacement amplitude is not greater than a preset proportional threshold, which is 0.2 in this embodiment.

[0168] Specifically, the method for determining whether the initial matching pair satisfies the local normal consistency test is as follows: during feature matching, a local normal consistency test is performed on the initial matching point pair. For any matching point pair... Calculate the unit normal vector of each corresponding spatial point. and And calculate the included normal angle. when If the matching point pair satisfies the local normal consistency constraint, then it is determined that the pair satisfies the constraint; otherwise, it is determined to be a mismatched point pair across the slope and is discarded. This is the normal consistency threshold, used to measure the allowable deviation range between the surface normals corresponding to two feature points. When the angle between two normal vectors is less than this threshold, they are considered to be on the same or similar slope structure.

[0169] Step K210. Based on the epipolar geometric constraints, calculate the distance from each matching point pair to the epipolar line using the essential matrix; based on the semantic mask, for matching pairs in vegetated areas, discard them if the distance is greater than a first preset threshold; for matching pairs in rocky areas, discard them if the distance is greater than a second preset threshold; for matching pairs in bare soil areas, discard them if the distance is greater than a third preset threshold; the retained matching pairs are used as feature correspondences between images, and the feature correspondences contain multiple sets of matching point pairs.

[0170] In this embodiment, the second preset threshold is greater than the first preset threshold, and the first preset threshold is greater than the third preset threshold. For example, the first preset threshold (vegetated area) is 0.8 to 1.5 pixels, preferably 1.2 pixels; the second preset threshold (rock area) is 1.5 to 2.5 pixels, preferably 2.0 pixels; and the third preset threshold (bare soil area) is 0.5 to 1.0 pixels, preferably 0.8 pixels.

[0171] These implementations introduce regional constraints introduced by semantic segmentation to achieve regional mismatch elimination: For vegetated areas, based on the original nearest neighbor distance ratio matching + essential matrix geometric verification, dynamic optical flow matching constraints are superimposed. By calculating the optical flow information of adjacent frames, dynamically mismatched point pairs with excessive deviations between the displacement direction of feature points and optical flow are eliminated, solving the matching instability problem caused by vegetation swaying; For rocky areas, based on the original matching and geometric verification, a local normal consistency check is superimposed. By calculating the local normal vector of the rocky area, cross-slope mismatched point pairs with normal vector angles exceeding a preset threshold are eliminated, solving the matching false alarm problem caused by weak rock texture and sudden changes in illumination; For bare soil areas, the original standard feature matching and geometric verification strategy in step S200 is used, and the distance ratio criterion is adopted to achieve mismatch elimination. Finally, regionally accurate matching point pairs are obtained, completing the optimization and enhancement of the original "feature matching and geometric verification" step.

[0172] In these implementations, the terrain classification results obtained from semantic segmentation are integrated throughout the entire process of scale space construction, feature point detection, feature description, and feature matching. This achieves regional adaptive optimization under a unified feature extraction framework, specifically addressing the pain points of feature extraction in different mountainous terrain regions. It significantly improves the robustness of feature matching in complex mountainous environments, providing a reliable guarantee for improving the accuracy of subsequent 3D reconstruction. Furthermore, this collaborative scheme demonstrates clear technical integration and engineering feasibility. This method, by introducing a terrain semantic-driven mechanism into a unified feature extraction framework, achieves a transformation from "general feature extraction" to "terrain-aware feature extraction." Experiments show that this method improves the matching success rate by approximately 37% in areas with weak texture and vegetation, and reduces the mismatch rate by 42% in shadowy rocky areas. Regarding overall reconstruction accuracy, the average reprojection error of the sparse point cloud is only 0.43 pixels, the camera pose drift rate is reduced to 5.8%, and the geometric integrity rate of the dense model reaches 96.2%. Compared with existing methods, it effectively solves core problems such as unstable matching, severe pose drift, and incomplete reconstruction in mountainous, weak-textured, and shadow-occluded scenes.

[0173] Step S300. Estimate the camera pose based on the pose and position information corresponding to the first frame image; for each new image, perform joint optimization on the new image and its co-view associated images by local bundle adjustment, and update the camera pose and 3D point coordinates; construct a sparse 3D point cloud and camera pose sequence based on the updated camera pose and 3D point coordinates.

[0174] In some embodiments of this example, step S300 specifically includes steps S310 to S380.

[0175] Step S310. Obtain the feature correspondence between images, wherein the feature correspondence includes multiple sets of matching point pairs.

[0176] Step S320. For the first image, estimate its camera pose based on the synchronously recorded UAV position and attitude information, or by using the homography matrix decomposition method.

[0177] Step S330. For each new image, the camera pose of the new image is solved using the PnP algorithm based on the matching point pairs between the new image and the reconstructed image.

[0178] Step S340. Based on the solved camera pose, triangulate the unreconstructed matching point pairs between the current image and the reconstructed image to generate new 3D point coordinates.

[0179] Step S350. Merge the newly added 3D point coordinates with the existing 3D point coordinates to update the 3D point cloud.

[0180] Step S360. Perform local bundle adjustment, and construct a local window by associating the newly added image with the preset number of images with the highest co-view degree. Optimize the camera pose and 3D point coordinates of all images in the local window by minimizing the reprojection error.

[0181] Step S370. Repeat steps S330 to S360 until all images have been processed.

[0182] Step S380. Construct a sparse 3D point cloud and camera pose sequence based on the updated camera pose and 3D point coordinates.

[0183] Step S400. Initialize each 3D point as a 3D Gaussian primitive based on the sparse 3D point cloud. The 3D Gaussian primitive includes Gaussian center position parameters, rotation parameters, scaling parameters, color sphere harmonic function coefficients, and opacity parameters.

[0184] Specifically, the Gaussian center position parameter μ is initialized to the spatial coordinates of the 3D point; the rotation parameter is initialized to the identity matrix; the scaling parameter is initialized to an isotropic matrix, for example, its scale factor is 0.5 times the average distance between the 3D point and its neighboring 3D points; the opacity parameter α is initialized to a set value, for example, 0.5; and the color spherical harmonic function coefficients are obtained by mapping the image pixel color values ​​corresponding to the 3D point.

[0185] The geometry of each 3D Gaussian element is described by the covariance matrix Σ. The general form of the covariance matrix Σ is:

[0186]

[0187] In the formula, Indicates the first A three-dimensional Gaussian element in spatial position The response value or probability density value at that location; For input multidimensional sample vectors; For the first The mean vector of Gaussian elements (i.e., the center position of Gaussian). This is the covariance matrix, used to describe the shape and orientation of Gaussian elements; It is the inverse of the covariance matrix; This is the transpose of the deviation vector.

[0188] To facilitate the subsequent optimization process of learning the pose and shape of the Gaussian ellipsoid separately, the covariance matrix is ​​decomposed into a rotation matrix. and scaling matrix The product form ( Represents a special orthogonal group in three dimensions; Indicates the first A rotation matrix of three-dimensional Gaussian elements is used to describe the orientation of the Gaussian ellipsoid in three-dimensional space; Represents a set of real number elements Matrix space; Indicates the first The scaling matrix of three-dimensional Gaussian elements (used to describe the scale of the Gaussian ellipsoid along the three principal axes) can be represented as:

[0189]

[0190] In the formula, Indicates the first A rotation matrix of three-dimensional Gaussian elements is used to describe the orientation of the Gaussian ellipsoid in three-dimensional space; Represents the transpose of the rotation matrix; Indicates the first A scaling matrix of three-dimensional Gaussian elements is used to describe the scale of the Gaussian ellipsoid along the three principal axes.

[0191] Under normal circumstances, For diagonal matrices:

[0192]

[0193] In the formula, They represent the first The scaling of a Gaussian element along three principal axes.

[0194] Based on the above parameters, the probability density function of the three-dimensional Gaussian element in space is defined as:

[0195]

[0196] In the formula, Representing multidimensional random variables The probability density value; The mean vector representing the distribution determines the "central location" of the multivariate normal distribution.

[0197] Colors are represented using spherical harmonic functions, the general expression of which is:

[0198]

[0199] In the formula, the spherical harmonic function is usually expressed as: ; n represents the order of the spherical harmonic function, and ; Let represent the degree of the spherical harmonic function, and ; Represents the polar angle in spherical coordinates; Represents the azimuth angle in spherical coordinates; This represents an order of n and a degree of n. The association of Legendre polynomials; Represents the natural constant; Represents the imaginary unit; This represents a complex exponential term related to the azimuth angle; It represents pi (π).

[0200] Depend on Figure 3 It is known that low-order spherical harmonic functions correspond to relatively smooth, low-frequency directional changes, and can be used to represent slow changes in overall illumination or color. As the order of the spherical harmonic function increases, the positive and negative lobe distribution of the spherical harmonic basis functions on the sphere becomes more complex, enabling the expression of higher-frequency and more refined view-dependent color changes. This embodiment uses spherical harmonic functions to model the color of three-dimensional Gaussian primitives, allowing each Gaussian primitive to present different color responses under different viewing directions, thereby improving the ability to express the view-dependent appearance of areas such as rocks, vegetation, and shadows in complex mountain scenes.

[0201] Step S500. Calculate the rate of change of the normal in the neighborhood of each 3D point in the sparse 3D point cloud; if it is greater than the preset steepness threshold, adjust the scaling matrix to make the Gaussian ellipsoid lengthen along the slope normal and shorten along the tangent, and align the principal axis of the rotation matrix with the slope direction to form an anisotropic distribution that is lengthened along the slope; otherwise, set the scaling matrix to a multiple of the identity matrix to make the Gaussian ellipsoid be isotropic.

[0202] Specifically, the rate of change of the normal vector within the neighborhood of each 3D point in the sparse 3D point cloud is calculated. When the rate of change of the normal vector is greater than a preset steepness threshold, the point is determined to be located in a steep terrain region. The scaling matrix is ​​adjusted so that it has a first scale along the slope normal direction and a second scale along the slope tangential direction, with the first scale being greater than the second scale. The rotation matrix is ​​then used to align the principal axis of the Gaussian ellipsoid with the principal direction of the slope, forming an anisotropic distribution elongated along the slope. When the rate of change of the normal vector is less than or equal to the preset steepness threshold, the point is determined to be located in a non-steep terrain region. The scaling matrix is ​​set to a multiple of the identity matrix, so that the Gaussian ellipsoid exhibits an isotropic distribution.

[0203] In some embodiments of this example, calculating the rate of change of the normal vector within the neighborhood of each 3D point in the sparse 3D point cloud includes: for each 3D point in the sparse 3D point cloud, selecting adjacent 3D points within a preset radius to form a spatial neighborhood of the 3D point; performing plane fitting on all 3D points within the spatial neighborhood using principal component analysis or least squares method to obtain a local plane and the normal vector of the local plane; and calculating the rate of change of the normal vector of the 3D point based on the difference between the normal vector and the normal vectors of each adjacent 3D point within the spatial neighborhood.

[0204] The preset steepness threshold is a pre-set constant.

[0205] The preset radius range of the three-dimensional points is set by the density of the sparse three-dimensional point cloud, for example, by taking 2 to 3 times the average spacing of the point cloud.

[0206] In this embodiment, the initial scaling parameters and rotation matrix of the Gaussian distribution are dynamically adjusted according to the local normal variation rate of the sparse three-dimensional point cloud and the degree of terrain undulation, so that the Gaussian ellipsoid has stronger directionality in steep terrain areas (stretched along the slope) and tends to be isotropic in flat areas, thereby more accurately approximating complex terrain surfaces.

[0207] Step S600. Project the three-dimensional Gaussian primitives onto the two-dimensional image plane using the view transformation matrix to obtain the two-dimensional center point. Project the three-dimensional Gaussian covariance matrix using the Jacobian matrix to obtain the two-dimensional covariance. The two constitute the two-dimensional Gaussian distribution primitives. Perform opacity-based weighted accumulation on the two-dimensional Gaussian distribution primitives and render to generate a two-dimensional image.

[0208] Specifically, the three-dimensional Gaussian primitives are projected onto a two-dimensional image plane using a view transformation matrix to obtain a two-dimensional center point, and the three-dimensional Gaussian covariance matrix is ​​transformed by a Jacobian matrix to obtain a two-dimensional covariance. The two-dimensional center point and the two-dimensional covariance constitute a two-dimensional Gaussian distribution primitive. The two-dimensional Gaussian distribution primitive is then weighted and accumulated based on opacity on the two-dimensional image plane to render and generate a two-dimensional image.

[0209] In some embodiments of this example, step S600 specifically includes steps S610 to S650.

[0210] Step S610. For each 3D Gaussian element, multiply its center position coordinates by the view transformation matrix to obtain the projection coordinates of the 3D Gaussian element on the 2D image plane, which is used as the 2D center point.

[0211] Step S620. Perform multiplication operations on the three-dimensional Gaussian covariance matrix of the three-dimensional Gaussian element with the view transformation matrix and the Jacobian matrix in sequence to obtain the two-dimensional covariance after projection.

[0212] The two-dimensional covariance matrix is:

[0213]

[0214] In the formula, J is the Jacobian matrix of the affine approximation of the projection transformation, and W represents the view transformation matrix.

[0215]

[0216]

[0217] In the formula, , , These represent scale, rotation, and displacement, respectively.

[0218] Step S630. Based on the two-dimensional center point and the two-dimensional covariance, construct the two-dimensional Gaussian distribution primitive of the Gaussian primitive on the two-dimensional image plane.

[0219] Step S640. For each pixel on the two-dimensional image plane, determine all two-dimensional Gaussian distribution primitives covering the pixel and sort them in ascending order of depth value to achieve a fusion order from near to far.

[0220] Step S650. Following the sorted order, the color value and opacity of each two-dimensional Gaussian distribution primitive are successively fused to obtain the final color value of the pixel, thereby generating a two-dimensional image.

[0221] The fusion formula is:

[0222]

[0223] In the formula, This represents the final color value of the target pixel, and N represents the N ordered points that can influence i pixels within a certain range. This represents the color of the i-th pixel. This represents the opacity of the current point i.

[0224] Step S700. Construct a joint loss function between the rendered 2D image and the real image, and back-optimize the parameters of the 3D Gaussian unit until the convergence condition is met to obtain a static 3D Gaussian model.

[0225] The two-dimensional image is the image generated in step S600, and the real image is the image acquired in step S100.

[0226] In some embodiments of this example, step S700 specifically includes steps S710 to S720.

[0227] Step S710. Construct a joint loss function that includes the L1 loss function and the structural similarity loss function.

[0228] The joint loss function is:

[0229]

[0230] In the formula, L is the joint loss function, and L1 is the L1 loss function. Let λ be the structural similarity loss function, and λ be the weight coefficient.

[0231] Step S720. The adaptive moment estimator optimizer is used to backpropagate the gradient of the loss function to each parameter of the Gaussian unit through the differentiable rendering pipeline, and the parameters are iteratively updated until the loss function converges.

[0232] Step S800. Generate incremental 3D reconstruction results based on the static 3D Gaussian model.

[0233] In some embodiments of this example, step S800 specifically includes steps S810 to S840.

[0234] Step S810. Construct a dynamic Gaussian mapping network. Based on the static three-dimensional Gaussian model, decompose and map the four-dimensional spatiotemporal coordinates into six two-dimensional feature planes. The six two-dimensional feature planes include three pure spatial planes and three spatiotemporal planes.

[0235] The tensors of each characteristic plane are defined as follows:

[0236]

[0237] in, Represents a set of three purely spatial feature planes, specifically including , respectively corresponding Three spatial coordinate planes; Represents three sets of spatiotemporal feature planes, specifically including , respectively corresponding Three spacetime coordinate planes.

[0238] Step S820. Extract and fuse features from each two-dimensional feature plane using a lightweight multilayer perceptron decoder, and output the displacement, rotation, and scaling offsets of each three-dimensional Gaussian element in the time dimension.

[0239] The parameter expression for the dynamic Gaussian element is:

[0240]

[0241]

[0242] In the formula, Represents dynamic Gaussian elements; Represents a static three-dimensional Gaussian element; This represents the offset of the Gaussian meta-parameters output by the dynamic Gaussian mapping network; This represents the dynamically updated center position parameter; This represents the dynamically updated rotation parameters; This represents the dynamically updated scaling parameters; and These represent opacity and color, respectively.

[0243] Static 3D Gaussian elements It can be represented as:

[0244]

[0245] In the formula, This represents the center position parameter of a static three-dimensional Gaussian element; The rotation parameters represent the static three-dimensional Gaussian elements; The scaling parameter represents the static 3D Gaussian primitive; This represents the opacity parameter; This represents the color parameter.

[0246] dynamic offset It can be represented as:

[0247]

[0248] In the formula, Indicates the offset of the center position; Indicates the rotation offset; This indicates the scaling offset.

[0249] Therefore, dynamic Gaussian elements It can be represented as:

[0250]

[0251] In the formula, This represents the center position parameter of a static three-dimensional Gaussian element; The rotation parameters represent the static three-dimensional Gaussian elements; The scaling parameter represents the static 3D Gaussian primitive; This represents the opacity parameter; Indicates color parameters; Indicates the offset of the center position; Indicates the rotation offset; This indicates the scaling offset.

[0252] Step S830. Use terrain features to guide the attention mechanism of the dynamic Gaussian deformation network. The terrain features include at least one of slope, aspect, and vegetation cover. The vegetation cover is obtained by semantic segmentation of the color and texture features of the aerial image. The loss weight guided by terrain semantics makes the dynamic deformation modeling focus on dense vegetation areas and cliff edge areas.

[0253] In other embodiments of this example, the dynamic Gaussian mapping network employs a terrain-driven dynamic activation mechanism. A dynamic importance weight map is generated based on the terrain slope map, roughness map, and vegetation cover density map obtained from static reconstruction. In densely vegetated areas, the dynamic network samples at a higher temporal resolution, allowing displacement offsets of up to 0.5 meters; in rocky areas, the dynamic sampling frequency is reduced by 80%, and rigid constraints are applied. The dynamic modeling results are adjusted for static Gaussian covariance via a feedback optimization module: vegetation sway statistics are used to correct the anisotropy of the static model, making the reconstructed vegetation area more consistent with natural morphology.

[0254] Step S840. Apply dynamic temporal smoothing regularization constraints and terrain consistency constraints. Temporal smoothing is achieved by applying Laplace filtering to the eigenvalues ​​of the spatiotemporal plane, and terrain consistency constraints are achieved through regularization terms. In vegetated areas, the magnitudes of displacement offset, rotation offset, and scaling offset are correlated with the terrain slope and vegetation type, and constraints are applied using a pre-statistical vegetation sway model. In rocky areas, displacement offset, rotation offset, and scaling offset are forced to be zero to maintain static conditions. In terrain transition areas, deformation parameters are smoothly transitioned between dynamic and static areas through bilinear interpolation to construct a dynamic three-dimensional Gaussian model.

[0255] Dynamic temporal smoothing regularization constraints are implemented using Laplace filtering, and their expression is:

[0256]

[0257] In the formula, This represents a dynamic temporal smoothing regularization term; Represents a set of spatiotemporal characteristic planes; Represents the set of spatiotemporal planes participating in temporal smoothing constraints; This represents the number of planes in the set of spatiotemporal planes; Represents one of the spacetime planes; Indicates a spatial location index in the spacetime plane; Indicates a time index; This represents the sampling resolution of the spatiotemporal plane in a single spatial dimension. This represents the total number of spatial sampling locations on each spatiotemporal plane; Representing the spacetime plane In spatial location and time Eigenvalues ​​at; and Representing the same spacetime plane Spatial location The feature values ​​of the adjacent previous and adjacent next time steps; This represents the square of the L2 norm.

[0258] Furthermore, set the spatiotemporal plane constraint terms:

[0259]

[0260] In the formula, Represents the spatiotemporal plane constraint regularization term; Indicates the first A spatiotemporal characteristic plane; Indicates and Tensors with the same dimensions and all element values; It represents the first norm, used to constrain the value distribution of the spatiotemporal characteristic plane.

[0261] Figure 4 This diagram illustrates one possible result of incremental 3D reconstruction using the method described in this embodiment. The result fully preserves the natural undulations of the mountainous terrain (such as ridges, valleys, gullies, terraces, cliffs, and other micro-topographical units) and accurately delineates areas of abrupt terrain change (such as steep slopes, rock edges, and boundaries of man-made structures). Semantic classification labels (such as vegetated areas, rocky areas, bare soil areas, and roads) can be overlaid on the reconstruction results. The reconstruction results can be directly applied to subsequent processing such as mesh reconstruction, 3D Gaussian modeling, and volume calculation.

[0262] Figure 5This diagram illustrates the effect of edge enhancement and smoothing on the 3D reconstruction results obtained by the method in this embodiment. This processing strengthens terrain edge features (such as cliff edges, rock edges, and road boundaries), making the model outline clearer and more three-dimensional, avoiding a blurry, "blurred" appearance. Simultaneously, it eliminates high-frequency noise (such as point cloud jitter and reconstruction artifacts), resulting in natural transitions between slopes and gentle slopes, avoiding jagged or uneven rough surfaces. Ultimately, it yields a 3D model with clear outlines, a smooth surface, and rich details, improving visual presentation quality while preserving the true undulating features of the mountains.

[0263] In one example, a real-world aerial photography area in a mountainous region was selected as the test scenario. The flight altitude was 120m, the forward overlap rate was 80%, and the lateral overlap rate was 70%, resulting in the acquisition of 2100 frames of images. The UAV flight platform was equipped with an image acquisition module, an inertial measurement unit (IMU), and a GPS / RTK positioning module, simultaneously acquiring multi-view monocular image sequences and corresponding attitude and position information.

[0264] Simulation results show that: the average reprojection error of sparse point clouds is less than 0.43 pixels; the geometric integrity of dense models exceeds 96%; stable reconstruction is still possible in steep slopes, woodlands, and deep shadow areas; camera pose drift is reduced by about 41% compared to the traditional SfM method; in the dynamic modeling part, the efficiency of dynamic region modeling is improved by about 35% and the physical consistency score of dynamic deformation is improved by 28% compared to the general 4DGS method, and the overall modeling efficiency is improved by about 2.3 times compared to the traditional MVS method.

[0265] The second aspect of this embodiment discloses an incremental 3D reconstruction system for unmanned aerial vehicles (UAVs) in complex mountainous environments, such as... Figure 6 As shown, the UAV incremental 3D reconstruction system includes an image acquisition module, a feature matching module, a 3D point cloud construction module, a 3D Gaussian primitive generation module, a normal change rate calculation module, a 2D image generation module, a static model generation module, and a 3D reconstruction module.

[0266] The image acquisition module is used to simultaneously acquire multi-view monocular images and corresponding attitude and position information when the UAV is executing a preset flight path.

[0267] The feature matching module is used to extract feature points from the image and perform matching to establish feature correspondence between images.

[0268] The 3D point cloud construction module is used to estimate the camera pose based on the pose and position information corresponding to the first frame image; for each new image, the new image and its co-view associated images are jointly optimized by local bundle adjustment to update the camera pose and 3D point coordinates; and a sparse 3D point cloud and camera pose sequence are constructed based on the updated camera pose and 3D point coordinates.

[0269] The 3D Gaussian primitive generation module is used to initialize each 3D point as a 3D Gaussian primitive based on a sparse 3D point cloud. The 3D Gaussian primitive includes Gaussian center position parameters, rotation parameters, scaling parameters, color sphere harmonic function coefficients, and opacity parameters.

[0270] The normal change rate calculation module is used to calculate the normal change rate in the neighborhood of each 3D point in the sparse 3D point cloud. If it is greater than the preset steepness threshold, the scaling matrix is ​​adjusted to make the Gaussian ellipsoid lengthen along the slope normal and shorten along the tangent, and the principal axis of the rotation matrix is ​​aligned with the slope direction to form an anisotropic distribution that is lengthened along the slope. Otherwise, the scaling matrix is ​​set to a multiple of the identity matrix to make the Gaussian ellipsoid be isotropic.

[0271] The two-dimensional image generation module is used to project three-dimensional Gaussian primitives onto a two-dimensional image plane using a view transformation matrix to obtain a two-dimensional center point, and to project a three-dimensional Gaussian covariance matrix using a Jacobian matrix to obtain a two-dimensional covariance. The two constitute two-dimensional Gaussian distribution primitives. The two-dimensional Gaussian distribution primitives are then weighted and accumulated based on opacity to render and generate a two-dimensional image.

[0272] The static model generation module is used to construct a joint loss function between the two-dimensional image and the real image, and to inversely optimize the parameters of the three-dimensional Gaussian unit until convergence, thereby obtaining a static three-dimensional Gaussian model.

[0273] The 3D reconstruction module is used to generate incremental 3D reconstruction results based on a static 3D Gaussian model.

[0274] It should be noted that those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working process of the system or device described above can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here.

[0275] The above description is merely a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the forms disclosed herein and should not be construed as excluding other embodiments. It can be used in various other combinations, modifications, and environments, and can be altered within the scope of the concept described herein through the above teachings or related technologies or knowledge. Modifications and variations made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the protection scope of the appended claims.

Claims

1. A method for incremental 3D reconstruction of unmanned aerial vehicles in complex mountainous environments, characterized by, The method comprises the following steps: controlling the unmanned aerial vehicle to fly along a preset flight path, and collecting multi-view monocular images and corresponding attitude and position information; extracting feature points of the images and performing matching to establish a feature correspondence relationship between the images; estimating the camera pose of the first image based on the attitude and position information corresponding to the first image; for each new image, performing joint optimization on the new image and its co-visibility associated images through local bundle adjustment to update the camera pose and three-dimensional point coordinates; and constructing a sparse three-dimensional point cloud and a camera pose sequence based on the updated camera pose and three-dimensional point coordinates; initializing each three-dimensional point in the sparse three-dimensional point cloud as a three-dimensional Gaussian primitive, wherein the three-dimensional Gaussian primitive comprises a Gaussian center position parameter, a rotation parameter, a scaling parameter, a color spherical harmonic function coefficient and an opacity parameter; calculating the normal variation rate in the neighborhood of each three-dimensional point in the sparse three-dimensional point cloud; if the normal variation rate is greater than a preset steepness threshold, adjusting the scaling matrix to make the Gaussian ellipsoid elongate along the slope normal and shorten in the tangential direction, and aligning the principal axis of the rotation matrix with the slope direction to form an anisotropic distribution along the slope; otherwise, setting the scaling matrix as a multiple of the unit matrix to make the Gaussian ellipsoid have an isotropic distribution; projecting the three-dimensional Gaussian primitive to a two-dimensional image plane through a view transformation matrix to obtain a two-dimensional center point, and projecting the three-dimensional Gaussian covariance matrix through a Jacobian matrix to obtain a two-dimensional covariance, which together form a two-dimensional Gaussian distribution primitive; performing weighted accumulation on the two-dimensional Gaussian distribution primitive based on the opacity to render a two-dimensional image; constructing a joint loss function of the two-dimensional image and the real image, and reversely optimizing the parameters of the three-dimensional Gaussian primitive to convergence to obtain a static three-dimensional Gaussian model; and generating an incremental three-dimensional reconstruction result based on the static three-dimensional Gaussian model.

2. The method of claim 1, wherein, extracting feature points of the images and performing matching to establish a feature correspondence relationship between the images, comprising: calculating the slope angle and the average curvature of the terrain surface, wherein the slope angle is calculated through the first-order partial derivative of the elevation function, and the average curvature is calculated through the second-order partial derivative of the elevation function; determining an adjusted scale parameter based on the slope angle and the average curvature, wherein the adjusted scale parameter increases when the slope angle increases and increases when the absolute value of the average curvature increases; substituting the adjusted scale parameter into a Gaussian kernel function to generate an adaptive scale space, detecting local extreme points in the adaptive scale space, and obtaining candidate feature points; performing cross-image matching on the candidate feature points extracted in different images to obtain initial matching pairs; verifying and removing false matches in the initial matching pairs through an essential matrix based on epipolar geometry constraints to establish a feature correspondence relationship between the images, wherein the feature correspondence relationship comprises a plurality of matching point pairs.

3. The method of claim 1, wherein, extracting feature points of the images and performing matching to establish a feature correspondence relationship between the images, comprising: performing pixel-level terrain classification on the images using a semantic segmentation network to generate semantic masks of vegetation regions, rock regions and bare soil regions; calculating the slope angle and the average curvature of the terrain surface, wherein the slope angle is calculated through the first-order partial derivative of the elevation function, and the average curvature is calculated through the second-order partial derivative of the elevation function; The foundation dimensional parameters are determined based on the slope angle and the average curvature; the foundation dimensional parameters increase as the slope angle increases and as the absolute value of the average curvature increases. Adjusted scale parameters are generated based on the basic scale parameters: for vegetated areas, the basic scale parameters are weighted by superimposing the normalized vegetation index and used as the adjusted scale parameters; for rocky areas, the basic scale parameters are weighted by increasing the curvature weight coefficient and used as the adjusted scale parameters; for bare soil areas, the basic scale parameters are directly used as the adjusted scale parameters. The adjusted scale parameters are substituted into the Gaussian kernel function to generate an adaptive scale space, and a regional adaptive Gaussian scale space is generated according to the semantic mask. Local extreme points are detected in the regional adaptive Gaussian scale space, and the detected extreme points are used as candidate feature points: For vegetated areas, the scale space image is weighted based on the normalized vegetation index (NVI) to detect local extreme points, and extreme points with the NVI below a preset threshold are filtered out; For rocky areas, the scale space image is stretched in grayscale to detect local extreme points, and extreme points located in areas spanning slopes are removed; For bare soil areas, local extreme points are detected directly. Based on the semantic mask, feature descriptors are generated for the candidate feature points: for vegetated areas, SIFT descriptors are generated based on the normalized vegetation index and optical flow information between adjacent frames is fused; for rocky areas, SIFT-AID descriptors are generated; for bare soil areas, standard SIFT descriptors are generated. Candidate feature points extracted from different images are matched across images to obtain initial matching pairs; The initial matching pairs are filtered according to the semantic mask: for initial matching pairs in vegetated areas, initial matching pairs that do not meet the dynamic optical flow matching constraints are removed; for initial matching pairs in rock areas, initial matching pairs that do not meet the local normal consistency test are removed; for initial matching pairs in bare soil areas, all initial matching pairs are retained. Based on epipolar geometry constraints, the distance from each matching point pair to the epipolar line is calculated using the essential matrix. According to the semantic mask, for matching pairs in vegetated areas, pairs are discarded if the distance is greater than a first preset threshold; for matching pairs in rocky areas, pairs are discarded if the distance is greater than a second preset threshold; for matching pairs in bare soil areas, pairs are discarded if the distance is greater than a third preset threshold. The retained matching pairs are used as feature correspondences between images, and these feature correspondences contain multiple sets of matching point pairs. The second preset threshold is greater than the first preset threshold, and the first preset threshold is greater than the third preset threshold.

4. The method of claim 1, wherein, For each new image, local bundle adjustment is used to jointly optimize the new image and its associated co-view images, updating the camera pose and 3D point coordinates, including: For each new image, the camera pose of the new image is solved using the PnP algorithm based on the matching point pairs between the new image and the reconstructed image. Based on the obtained camera pose, the unreconstructed matching point pairs between the current image and the reconstructed image are triangulated to generate the coordinates of the newly added 3D points. The newly added 3D point coordinates are merged with the existing 3D point coordinates to update the 3D point cloud; Local bundle adjustment is performed, and the newly added image and the preset number of associated images with the highest co-view degree are used to form a local window. The camera pose and 3D point coordinates of all images in the local window are optimized by minimizing the reprojection error.

5. The method of claim 1, wherein, Calculate the rate of change of the normal vector within the neighborhood of each 3D point in a sparse 3D point cloud, including: For each three-dimensional point in the sparse three-dimensional point cloud, the spatial neighborhood of the three-dimensional point is formed by selecting the adjacent three-dimensional points within a preset radius. By performing plane fitting on all three-dimensional points in the spatial neighborhood using principal component analysis or least squares method, a local plane and its normal vector are obtained. The rate of change of the normal vector of a three-dimensional point is calculated based on the difference between the normal vector and the normal vectors of each adjacent three-dimensional point in the spatial neighborhood.

6. The method of claim 1, wherein, The three-dimensional Gaussian primitives are projected onto the two-dimensional image plane using a view transformation matrix to obtain the two-dimensional center point. The three-dimensional Gaussian covariance matrix is ​​then projected using a Jacobian matrix to obtain the two-dimensional covariance. These two components constitute the two-dimensional Gaussian distribution primitives. An opacity-based weighted accumulation is performed on these two-dimensional Gaussian distribution primitives, and a two-dimensional image is rendered, including: For each 3D Gaussian primitive, multiply its center position coordinates by the view transformation matrix to obtain the projected coordinates of the 3D Gaussian primitive on the 2D image plane, which serve as the 2D center point. The three-dimensional Gaussian covariance matrix of the three-dimensional Gaussian element is multiplied sequentially with the view transformation matrix and the Jacobian matrix to obtain the two-dimensional covariance after projection. Based on the two-dimensional center point and the two-dimensional covariance, construct the two-dimensional Gaussian distribution primitive of the Gaussian primitive on the two-dimensional image plane; For each pixel on the two-dimensional image plane, determine all two-dimensional Gaussian distribution primitives covering that pixel and sort them in ascending order of depth value; Following the sorted order, the color value and opacity of each two-dimensional Gaussian distribution primitive are successively fused to obtain the final color value of the pixel, thus generating a two-dimensional image.

7. The method of claim 1, wherein, Construct a joint loss function for the 2D image and the real image, and inversely optimize the parameters of the 3D Gaussian trigonometric function until convergence, including: Construct a joint loss function that includes the L1 loss function and the structural similarity loss function; An adaptive moment estimation optimizer is used to backpropagate the gradient of the loss function to each parameter of the Gaussian unit through a differentiable rendering pipeline, and iteratively update the parameters until the loss function converges.

8. The method of claim 1, wherein, Incremental 3D reconstruction results are generated based on the static 3D Gaussian model, including: A dynamic Gaussian mapping network is constructed. Based on the static three-dimensional Gaussian model, the four-dimensional spatiotemporal coordinates are decomposed and mapped into six two-dimensional feature planes. The six two-dimensional feature planes include three pure spatial planes and three spatiotemporal planes. A lightweight multilayer perceptron decoder is used to extract and fuse features from each two-dimensional feature plane, and output the displacement, rotation, and scaling offsets of each three-dimensional Gaussian element in the time dimension. The attention mechanism of the dynamic Gaussian deformation network is guided by terrain features, which include at least one of slope, aspect, and vegetation cover. The vegetation cover is obtained by semantic segmentation of the color and texture features of aerial images. The loss weight guided by terrain semantics makes the dynamic deformation modeling focus on dense vegetation areas and cliff edge areas. Dynamic temporal smoothing regularization constraints and terrain consistency constraints are applied. Temporal smoothing is achieved by applying Laplace filtering to the eigenvalues ​​of the spatiotemporal plane, and terrain consistency constraints are achieved through regularization terms. In vegetated areas, the magnitudes of displacement offset, rotation offset, and scaling offset are constrained in relation to the terrain slope and vegetation type, and constraints are applied using a pre-statistical vegetation sway model. In rocky areas, displacement offset, rotation offset, and scaling offset are forced to be zero to maintain static conditions. In terrain transition areas, bilinear interpolation is used to smoothly transition deformation parameters between dynamic and static regions, constructing a dynamic three-dimensional Gaussian model.

9. A UAV incremental 3D reconstruction system for complex mountainous environments, characterized in that, include: The image acquisition module is used to simultaneously acquire multi-view monocular images and corresponding attitude and position information when the UAV is flying along a preset route; The feature matching module is used to extract feature points from images and perform matching to establish feature correspondences between images. The 3D point cloud construction module is used to estimate the camera pose based on the pose and position information corresponding to the first frame image. For each new image, the new image and its co-view associated images are jointly optimized by local bundle adjustment to update the camera pose and 3D point coordinates. Based on the updated camera pose and 3D point coordinates, a sparse 3D point cloud and camera pose sequence are constructed. The 3D Gaussian primitive generation module is used to initialize each 3D point as a 3D Gaussian primitive based on a sparse 3D point cloud. The 3D Gaussian primitive includes Gaussian center position parameters, rotation parameters, scaling parameters, color spherical harmonic function coefficients, and opacity parameters. The normal change rate calculation module is used to calculate the normal change rate in the neighborhood of each 3D point in a sparse 3D point cloud. If it is greater than the preset steepness threshold, the scaling matrix is ​​adjusted to make the Gaussian ellipsoid lengthen along the slope normal and shorten along the tangent, and the principal axis of the rotation matrix is ​​aligned with the slope direction to form an anisotropic distribution that is lengthened along the slope. Otherwise, the scaling matrix is ​​set to a multiple of the identity matrix to make the Gaussian ellipsoid be isotropic. The 2D image generation module is used to project the 3D Gaussian primitives onto the 2D image plane through the view transformation matrix to obtain the 2D center point, and to project the 3D Gaussian covariance matrix using the Jacobian matrix to obtain the 2D covariance. The two constitute the 2D Gaussian distribution primitives. The 2D Gaussian distribution primitives are weighted and accumulated based on opacity to render and generate a 2D image. The static model generation module is used to construct a joint loss function between the two-dimensional image and the real image, and to back-optimize the parameters of the three-dimensional Gaussian unit until convergence, thus obtaining a static three-dimensional Gaussian model. The 3D reconstruction module is used to generate incremental 3D reconstruction results based on a static 3D Gaussian model.

Images