Three-dimensional reconstruction method for canyon vegetation leaf expansion process based on image fusion

Abstract

The present application relates to the technical field of image processing, and further relates to a canyon vegetation leaf stretching process three-dimensional reconstruction method based on image fusion, which comprises the following steps: step one: integrating a laser radar and an optical camera on a UAV platform to obtain a multi-time point cloud data set and a multi-time optical image set; step two: setting the earliest time point cloud data in the multi-time point cloud data set as a reference point cloud to obtain a multi-time color point cloud set; step three: filtering ground points and performing connectivity clustering on the multi-time color point cloud set to obtain independent plant point clouds; step four: performing single-leaf segmentation and cross-time matching on the multi-time leaf point cloud set, and constructing a canyon vegetation leaf space-time stretching dynamic model after terrain correction. The present application effectively solves technical problems such as insufficient multi-time point cloud registration accuracy, difficulty in extracting plant skeletons, inaccuracy in leaf separation, and stretching measurement errors caused by terrain elevation differences in a complex canyon terrain environment.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a three-dimensional reconstruction method for the leaf extension process of canyon vegetation based on image fusion. Background Technology

[0002] In monitoring vegetation growth dynamics, traditional methods mainly rely on manual observation and ground quadrat surveys, which suffer from low efficiency, limited coverage, and insufficient temporal resolution. The introduction of remote sensing technology has provided a new means for monitoring vegetation growth dynamics. Researchers infer vegetation growth status by analyzing changes in vegetation indices through multi-temporal remote sensing images. However, methods based on two-dimensional remote sensing images struggle to obtain information on three-dimensional structural changes in vegetation and cannot accurately quantify the leaf extension process in three-dimensional space. The emergence of multi-temporal lidar scanning technology has brought a breakthrough to three-dimensional dynamic monitoring of vegetation. By comparing point cloud data from different periods, three-dimensional morphological changes in vegetation can be analyzed. However, existing multi-temporal point cloud registration techniques face severe challenges in complex terrain environments. Canyon terrain is characterized by steep slopes, significant elevation differences, and severe occlusion, leading to decreased GPS positioning accuracy for UAVs in canyon environments and significant deviations in directly georeferenced point cloud data. Traditional point cloud registration methods, such as the iterative nearest-point algorithm, require good initial alignment and are sensitive to initial pose deviations, easily getting trapped in local optima when initial deviations are large. Feature-matching-based registration methods estimate transformation parameters by extracting local geometric features from point clouds and establishing correspondences. However, the large number of occlusions and repetitive structures in canyon environments leads to extremely high mismatch rates in feature matching, with outliers typically accounting for 70% to 90%. Existing outlier removal methods, such as the random sampling consensus algorithm, experience a sharp performance drop when the outlier ratio exceeds 50%, making it difficult to obtain reliable registration results under high outlier ratio conditions.

[0003] Plant skeleton extraction is a crucial step in vegetation structure analysis. Existing skeleton extraction methods mainly include voxelization-based refinement methods, Laplace shrinkage-based methods, and graph theory-based methods. Voxelization-based refinement methods convert point clouds into voxel meshes and then peel away outer voxels layer by layer until the skeleton is obtained; however, this method is sensitive to point cloud density and has low computational efficiency. Laplace shrinkage-based methods shrink the point cloud towards the central axis by solving the Laplace equation, which can maintain a good topological structure, but it is sensitive to parameter settings and prone to breakage when dealing with complex branching structures. Graph theory-based methods obtain the skeleton by constructing an adjacency graph of the point cloud and extracting the shortest path tree, but they are sensitive to noise and have difficulty handling multi-branching structures. Furthermore, most existing methods are designed for single plants, facing the dual challenges of computational efficiency and segmentation accuracy when dealing with densely distributed vegetation communities in canyon environments. Summary of the Invention

[0004] The main objective of this invention is to provide a three-dimensional reconstruction method for the leaf extension process of canyon vegetation based on image fusion. This method effectively solves technical problems such as insufficient registration accuracy of multi-temporal point clouds, difficulty in extracting plant skeletons, inaccurate leaf separation, and extension measurement errors caused by terrain elevation differences in complex canyon terrain environments. It can accurately quantify the leaf extension process of canyon vegetation in three-dimensional space, providing a reliable technical means for monitoring the dynamic growth of vegetation and conducting ecological research in high mountain and canyon areas.

[0005] To solve the above problems, the technical solution of the present invention is implemented as follows: A three-dimensional reconstruction method for the leaf extension process of vegetation in canyons based on image fusion, the method including: Step 1: Integrate LiDAR and optical camera on the UAV platform, achieve time synchronization through hardware trigger signal, and repeatedly scan the target canyon vegetation area on multiple observation dates during the growing season to obtain multi-temporal point cloud dataset and multi-temporal optical image set; Step 2: The earliest point cloud data in the multi-temporal point cloud dataset is designated as the reference point cloud, and the remaining point cloud data from each temporal phase are designated as point clouds to be registered. Feature descriptors are extracted and matched from the reference and point clouds to be registered, forming an initial set of corresponding point pairs. Each corresponding point pair in the initial set is used as a node in the invariant measurement map, and the invariant measurement map is constructed based on distance invariance constraints. A branch-and-bound strategy is used on the invariant measurement map to search for the maximum clique to identify the set of interior points. Based on the set of interior points, a global coarse registration transformation is solved through singular value decomposition, followed by local fine registration to obtain a multi-temporal registration point cloud set. The multi-temporal registration point cloud set is fused with the corresponding temporal optical images to obtain a multi-temporal color point cloud set. Step 3: Perform ground point filtering and connectivity clustering on the multi-temporal color point cloud set to obtain independent plant point clouds; perform iterative shrinkage processing on the independent plant point clouds to obtain a skeleton candidate point set; construct a minimum spanning tree on the skeleton candidate point set to form the plant skeleton; separate the independent plant point clouds into leaf point clouds based on the distance of each point to the plant skeleton to obtain a multi-temporal leaf point cloud set. Step 4: Perform single-leaf segmentation and time-span matching on the multi-temporal leaf point cloud, calculate the extension velocity vector of each leaf, and construct a dynamic model of the spatiotemporal extension of canyon vegetation leaves after terrain correction.

[0006] Furthermore, feature descriptor extraction and matching are performed on the reference point cloud and the point cloud to be registered to form an initial set of corresponding point pairs. This includes: uniformly downsampling the reference point cloud and the current point cloud to be registered to obtain the reference sparse point cloud and the sparse point cloud to be registered; calculating a fast point feature histogram descriptor for each point in the reference sparse point cloud and the sparse point cloud to be registered; comparing the fast point feature histogram descriptor of each point in the sparse point cloud to be registered with the fast point feature histogram descriptor of all points in the reference sparse point cloud using Euclidean distance, and selecting the reference point with the smallest descriptor distance for each point to be registered as the corresponding point to form an initial set of corresponding point pairs.

[0007] Furthermore, an invariant measurement graph is constructed based on distance invariance constraints, including: for any two nodes, calculating the Euclidean distance between the registration point in the first node and the registration point in the second node, and the Euclidean distance between the reference point in the first node and the reference point in the second node; calculating the absolute value of the difference between the two Euclidean distances; when the absolute value of the difference is less than the preset distance invariance tolerance, adding an edge between the two nodes; and traversing all node pairs to complete the construction of the invariant measurement graph.

[0008] Furthermore, a branch-and-bound strategy is employed on the invariant measurement graph to search for the maximum clique to identify the set of interior points. This includes: initializing the current maximum clique as an empty set and initializing the candidate node set to all nodes in the invariant measurement graph; selecting the node with the highest degree from the candidate node set and adding it to the current search path; updating the candidate node set to the intersection of the common neighbor nodes of all nodes in the current search path; if the updated candidate node set is not empty, recursively executing the above selection and update process; if the updated candidate node set is empty and the number of nodes in the current search path is greater than the number of nodes in the current maximum clique, updating the current maximum clique to the current search path; backtracking and selecting the node with the second highest degree from the candidate node set to continue the search; when the sum of the remaining number of nodes in the candidate node set and the number of nodes already in the current search path is less than the number of nodes in the current maximum clique, pruning and abandoning the current branch; after the search is completed, the corresponding point pairs of the nodes contained in the current maximum clique constitute the set of interior points.

[0009] Furthermore, the global coarse registration transformation is solved based on the set of interior points using singular value decomposition, including: extracting the 3D coordinates of all points to be registered from the set of interior points to form a coordinate matrix of the points to be registered, and extracting the 3D coordinates of all corresponding reference points to form a reference point coordinate matrix; calculating the centroid coordinates of the coordinate matrices of the points to be registered and the reference point coordinate matrices respectively, and subtracting the corresponding centroid coordinates from the coordinates of each point in each coordinate matrix to obtain a centroid-free coordinate matrix; calculating the product of the transpose of the centroid-free coordinate matrix of the points to be registered and the centroid-free reference point coordinate matrix to obtain the covariance matrix; performing singular value decomposition on the covariance matrix to obtain a left singular vector matrix and a right singular vector matrix; calculating the product of the right singular vector matrix and the transpose of the left singular vector matrix to obtain a rotation matrix; subtracting the product of the rotation matrix and the centroid coordinates of the points to be registered from the centroid coordinates of the reference points to obtain a translation vector; the rotation matrix and the translation vector together constitute the global coarse registration transformation.

[0010] Furthermore, local fine registration is performed to obtain a multi-temporal registration point cloud, including: using the global coarse registration transformation as the initial value, local fine registration is performed using the point-to-surface iterative nearest point method to obtain the fine registration transformation; the global coarse registration transformation and the fine registration transformation are concatenated to obtain the final registration transformation from the current point cloud to the reference point cloud; the above registration process is repeated for all point clouds to be registered in the multi-temporal point cloud dataset to obtain a multi-temporal registration point cloud in a unified coordinate system; based on the camera calibration parameters, each point in the multi-temporal registration point cloud is projected onto the optical image of the corresponding temporal phase, and each point is assigned a color value of the corresponding pixel to obtain a multi-temporal color point cloud.

[0011] Furthermore, an iterative shrinkage process is performed on the independent plant point cloud to obtain a skeleton candidate point set, including: copying the current independent plant point cloud into a shrunken point cloud, and setting an initial shrinkage step size for each point in the shrunken point cloud; for each current point in the shrunken point cloud, searching for all points whose Euclidean distance to the current point is less than a preset neighborhood radius to form a neighborhood point set; for each neighborhood point in the neighborhood point set, calculating the unit direction vector from the current point to the neighborhood point, and calculating the reciprocal of the Euclidean distance between the current point and the neighborhood point as the contribution of the unit direction vector; summing and normalizing the unit direction vectors corresponding to all neighborhood points in the neighborhood point set according to their respective contributions to obtain the shrinkage direction vector of the current point; moving the current point along the shrinkage direction vector by the current shrinkage step size; determining whether the moving distance is less than a preset convergence threshold, and if it is less, marking it as converged; terminating the iteration when the proportion of the number of converged points to the total number of points exceeds a preset proportion threshold, otherwise reducing the shrinkage step size to half of the current value and continuing the iteration; after the iteration terminates, downsampling the shrunken point cloud to remove overlapping points to obtain a skeleton candidate point set.

[0012] Furthermore, a minimum spanning tree is constructed on the skeleton candidate point set to form the plant skeleton, including: constructing a minimum spanning tree on the skeleton candidate point set, wherein the edge of the minimum spanning tree is a line segment connecting two skeleton candidate points, and the length of the edge is the Euclidean distance between the two skeleton candidate points; in the minimum spanning tree, the node with a recognition degree of one is the endpoint, and the node with a recognition degree greater than two is the branch point; the edge sequence connecting adjacent branch points or between a branch point and an endpoint is merged into a skeleton branch, and all skeleton branches constitute the plant skeleton.

[0013] Furthermore, based on the distance from each point to the plant skeleton, the independent plant point cloud is separated into leaf point clouds, including: for each original point in the independent plant point cloud, the shortest Euclidean distance from the original point to the plant skeleton is calculated; when the shortest Euclidean distance is greater than or equal to a preset skeleton distance threshold, the original point is marked as a leaf point; all leaf points constitute a leaf point cloud; the skeleton extraction and leaf separation process is repeated for all phases in the multi-temporal color point cloud set to obtain a multi-temporal plant skeleton set and a multi-temporal leaf point cloud set.

[0014] Furthermore, single-leaf segmentation and time-series matching are performed on the multi-temporal leaf point cloud set. The extension velocity vector of each leaf is calculated, and after terrain correction, a spatiotemporal extension dynamic model of canyon vegetation leaves is constructed. This includes: performing curvature-based region segmentation on the current-temporal leaf point cloud set in the multi-temporal leaf point cloud set to obtain several single-leaf point clouds; calculating the centroid 3D coordinates and principal direction vector of each single-leaf point cloud; performing matching on the single-leaf point cloud sets of two adjacent temporal phases, using the Euclidean distance between the centroid 3D coordinates and the angle between the principal direction vectors as the matching cost, and using the Hungarian algorithm to obtain the correspondence of the same leaf between adjacent temporal phases; for each successfully matched leaf, calculating the difference between the centroid 3D coordinates of the later temporal phase and the centroid 3D coordinates of the previous temporal phase to obtain the corresponding leaf. The extension displacement vector is divided by the time interval between two time phases to obtain the extension velocity vector. Based on the three-dimensional coordinates of the plant skeletons at each time phase in the multi-time phase plant skeleton set, the local terrain slope direction of the three-dimensional coordinates of the centroid of each leaf is calculated. The extension velocity vector is decomposed into a component along the local terrain slope direction and a component perpendicular to the local terrain slope direction. The component perpendicular to the local terrain slope direction is retained as the terrain-corrected leaf extension velocity vector. The terrain-corrected leaf extension velocity vectors of all leaves are spatially organized according to the three-dimensional coordinates of the centroid of the leaves to construct a three-dimensional extension vector field of plant leaves. The three-dimensional extension vector fields of plant leaves from all observation time phases are integrated to obtain a spatiotemporal extension dynamic model of canyon vegetation leaves covering the growing season.

[0015] The three-dimensional reconstruction method for the leaf extension process of canyon vegetation based on image fusion of the present invention has the following beneficial effects: The present invention integrates a lidar and an optical camera on a UAV platform and uses a hardware trigger signal to achieve time synchronization between the two, ensuring strict alignment of point cloud data and optical images in the time domain. This multi-source sensor collaborative acquisition method fully leverages the complementary advantages of the high ranging accuracy of lidar and the rich texture information of optical cameras, providing a high-quality data foundation for subsequent point cloud colorization and semantic segmentation. The acquisition strategy of repeated scanning on multiple observation days during the growing season can completely capture the entire leaf extension process of vegetation from budding to full leaf stage, meeting the temporal resolution requirements for dynamic monitoring of vegetation growth.

[0016] This invention addresses the high mismatch rate caused by canyon occlusion by innovatively proposing an out-point-resistant registration method based on maximum clique search in invariant measurement graphs. By treating corresponding point pairs as graph nodes and constructing edges based on distance invariance constraints, the interior point identification problem is transformed into a maximum clique search problem in graph theory. A branch-and-bound strategy is employed for maximum clique search, and pruning operations significantly improve search efficiency. This method can accurately identify the true matching set of interior points even with a very high out-point ratio, effectively overcoming the performance degradation of traditional registration methods in canyon environments. Based on the interior point set, a global coarse registration transformation is solved through singular value decomposition, combined with a local fine registration strategy, achieving a significant improvement in registration accuracy from coarse to fine.

[0017] This invention employs an iterative shrinkage method to progressively shrink the point cloud of independent plants to a candidate skeleton point set. By calculating the shrinkage direction vector within the neighborhood of each point and summing and normalizing it according to its contribution, the point cloud naturally converges towards the locally dense branch centerline. A minimum spanning tree is constructed on the candidate skeleton point set to form the plant skeleton, which accurately reflects the topological connection relationship from the trunk to each level of branches. The leaf separation method based on the distance from each point to the plant skeleton is simple and effective, and can accurately separate the leaf point cloud from the branch point cloud. Attached Figure Description

[0018] Figure 1 This is a schematic diagram illustrating the process of constructing an invariant measurement map based on distance invariant constraints, as provided in an embodiment of the present invention. Figure 2 This is a schematic diagram of the process curve for searching the maximum clique on the invariance measurement map using a branch and bound strategy, provided by an embodiment of the present invention. Figure 3 A cross-sectional schematic diagram illustrating the principle of iterative shrinkage and skeleton extraction of plant point clouds provided in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the three-dimensional spatial relationship between plant skeleton extraction and leaf point cloud separation provided in an embodiment of the present invention. Detailed Implementation

[0019] A three-dimensional reconstruction method for the leaf extension process of vegetation in canyons based on image fusion, the method including: Step 1: Integrate LiDAR and optical camera on the UAV platform, achieve time synchronization through hardware trigger signal, and repeatedly scan the target canyon vegetation area on multiple observation dates during the growing season to obtain multi-temporal point cloud dataset and multi-temporal optical image set; Step 2: The earliest point cloud data in the multi-temporal point cloud dataset is designated as the reference point cloud, and the remaining point cloud data from each temporal phase are designated as point clouds to be registered. Feature descriptors are extracted and matched from the reference and point clouds to be registered, forming an initial set of corresponding point pairs. Each corresponding point pair in the initial set is used as a node in the invariant measurement map, and the invariant measurement map is constructed based on distance invariance constraints. A branch-and-bound strategy is used on the invariant measurement map to search for the maximum clique to identify the set of interior points. Based on the set of interior points, a global coarse registration transformation is solved through singular value decomposition, followed by local fine registration to obtain a multi-temporal registration point cloud set. The multi-temporal registration point cloud set is fused with the corresponding temporal optical images to obtain a multi-temporal color point cloud set. Step 3: Perform ground point filtering and connectivity clustering on the multi-temporal color point cloud set to obtain independent plant point clouds; perform iterative shrinkage processing on the independent plant point clouds to obtain a skeleton candidate point set; construct a minimum spanning tree on the skeleton candidate point set to form the plant skeleton; separate the independent plant point clouds into leaf point clouds based on the distance of each point to the plant skeleton to obtain a multi-temporal leaf point cloud set. Step 4: Perform single-leaf segmentation and time-span matching on the multi-temporal leaf point cloud, calculate the extension velocity vector of each leaf, and construct a dynamic model of the spatiotemporal extension of canyon vegetation leaves after terrain correction.

[0020] Accurate 3D reconstruction of the leaf extension process in canyon terrain faces numerous challenges. The shading effect created by the mountainsides on both sides of the canyon results in complex and variable lighting conditions. The multi-layered structure of the vegetation canopy makes it difficult for traditional single sensors to fully capture the spatial distribution information of branches and leaves. Furthermore, data matching between different growth stages suffers from numerous mismatches due to terrain undulations and changes in vegetation morphology. To address these issues, this implementation combines multi-source sensor collaborative acquisition with an anti-outside-point registration algorithm to achieve high-precision 3D reconstruction of the leaf extension process in canyons.

[0021] Regarding the selection of the UAV platform and payload configuration, this implementation uses a hexacoach UAV as the flight platform, with a maximum takeoff mass of 25 kg, a payload capacity of 8 kg, and an endurance of up to 45 minutes under full load conditions. The UAV platform integrates two types of sensors: a solid-state LiDAR and an industrial-grade optical camera. The solid-state LiDAR employs a non-repeating scanning mode, achieving a single-echo ranging accuracy better than 3 cm, a scanning field of view of 70 degrees x 77 degrees, and a point cloud acquisition rate of 240,000 points per second. The optical camera uses a global shutter sensor with 45 million effective pixels, a sensor size of 35.9 mm x 24 mm, a lens focal length of 35 mm, and a single-frame image resolution of 8256 pixels x 5504 pixels.

[0022] The hardware trigger signal time synchronization mechanism is crucial to ensuring strict time alignment between LiDAR point cloud data and optical images. This implementation uses a Field-Programmable Gate Array (FPGA) as the synchronization control core, employing a second pulse signal output from a Global Navigation Satellite System (GNSS) receiver as a unified time reference. The rising edge of the second pulse signal triggers a counter within the FPGA to reset. The counter operates at a clock frequency of 100 MHz, achieving a time resolution of 10 nanoseconds. The FPGA outputs trigger pulses to both the LiDAR and the optical camera according to a preset acquisition cycle, with the time deviation between the two trigger pulses controlled within 100 nanoseconds. Upon receiving the trigger pulse, the LiDAR initiates a point cloud scan, and the optical camera performs exposure acquisition, thus ensuring a strict time correspondence between the point cloud data and the optical image acquired simultaneously. In an optional implementation, if the UAV platform lacks FPGA trigger control capabilities, a software timestamp matching method can be used. During post-processing, the point cloud frame and image frame with the closest timestamps are aligned, achieving a time synchronization accuracy of approximately 10 to 50 milliseconds.

[0023] The planning of observation periods during the growing season directly affects the temporal resolution of leaf extension. This implementation method sets eight observation dates during the vegetation growing season from April to September each year, with an interval of 20 to 25 days between adjacent observation dates. The selection of observation dates must avoid rainy weather and require wind speeds below 5 meters per second to minimize the impact of leaf vibration on the quality of point cloud acquisition. Flight missions for each observation date are conducted between 10:00 AM and 2:00 PM local time, during which the solar altitude angle is high and the lighting conditions at the bottom of the canyon are relatively sufficient. The UAV scans the target canyon vegetation area along a preset route with a route spacing of 30 meters, a flight altitude of 80 meters above the ground, and a flight speed of 5 meters per second. The area covered by a single flight mission is approximately 500 meters by 500 meters. After eight repeated scans during the growing season, a multi-temporal point cloud dataset containing eight phases is obtained, with each phase containing approximately 200 three-dimensional points; simultaneously, a multi-temporal optical image set containing eight phases is also obtained, with each phase containing approximately 200 optical images with a 75% overlap.

[0024] Registration between multi-temporal point cloud data is fundamental to subsequent leaf extension analysis. Due to the obstruction of GPS signals by mountains in canyon environments, the direct geolocation accuracy of UAVs is limited. Furthermore, vegetation morphology changes between observation dates, leading to significant deviations in point cloud overlays based on direct geolocation. Therefore, refined point cloud registration processing is necessary. In this implementation, the earliest point cloud data from the multi-temporal point cloud dataset is set as the reference point cloud, specifically the data acquired during the first observation in April. The earliest point cloud data is chosen as the reference because at this time, the vegetation is in the early budding stage, the canopy structure is relatively simple, and stable features such as ground points and tree trunks are easier to identify, facilitating the establishment of a reliable registration reference framework. The point cloud data from the remaining seven temporal phases are sequentially set as point clouds to be registered, and are registered with the reference point cloud one by one in chronological order.

[0025] Before extracting and matching feature descriptors from the baseline point cloud and the point cloud to be registered, uniform downsampling is first performed to reduce computational complexity. Downsampling employs a voxel grid filtering method, dividing the 3D space into a cubic grid with a side length of 0.5 meters. All points within each grid are replaced by their centroid coordinates, thus downsampling the original point cloud from 200 million points to approximately 8 million points, obtaining the baseline sparse point cloud and the sparse point cloud to be registered. The voxel side length for downsampling needs to be adjusted according to the actual scene; too small a side length results in limited downsampling effectiveness, while too large a side length may lead to the loss of important geometric features. A voxel side length of 0.5 meters strikes a good balance between computational efficiency and feature preservation in canyon vegetation scenes.

[0026] For each point in both the baseline sparse point cloud and the sparse point cloud to be registered, a fast point feature histogram descriptor is calculated. The fast point feature histogram descriptor is a feature representation method based on local geometric relationships, capable of characterizing the spatial angular relationships between each point in the point cloud and its neighbors, and is invariant to rigid transformations. For the current point... First, search in the point cloud for... Centered on, with radius All points within a spherical neighborhood constitute the neighborhood point set. Radius The value is set to 2.5 meters, approximately five times the side length of the voxel, to ensure that the neighborhood contains a sufficient number of points for feature calculation. Next, the current point is calculated. With each neighboring point in the neighborhood point set The three angular features between them. Let the current point be... The normal vector is Neighborhood points The normal vector is ,from point to The unit vector is Define three angular features: the first angular feature Normal vector with vector The angle between them; the second angular feature Normal vector with vector The angle between them; the third angular feature Normal vector With normal vector The angle between the three angle features is calculated. The value ranges of the three angle features are divided into 11 intervals. The number of all neighboring points in the neighborhood point set falling into each interval is counted, forming a 33-dimensional histogram vector, which represents the current point. The fast point feature histogram descriptor is used. The normal vector is calculated using principal component analysis, which decomposes the coordinate matrix of the neighborhood point set into eigenvalues. The eigenvector corresponding to the smallest eigenvalue is the direction of the normal vector.

[0027] After completing the feature descriptor calculation, feature matching is performed to establish the correspondence between the baseline sparse point cloud and the sparse point cloud to be registered. For each point in the sparse point cloud to be registered, its fast point feature histogram descriptor is compared with the Euclidean distance of the fast point feature histogram descriptors of all points in the baseline sparse point cloud. Let the descriptor of the point to be registered be... The descriptor of the reference point is The Euclidean distance between the two is ,in and These represent the first and second descriptor vectors, respectively. Each point has one component. For each point to be registered, the reference point with the smallest descriptor distance is selected as the corresponding point. That is, the point with the smallest Euclidean distance to the descriptor of the current point to be registered is found in the reference sparse point cloud, and these two points are paired. All points in the sparse point cloud to be registered are traversed to form an initial set of corresponding point pairs. Due to the large amount of occlusion and repetitive structures in the canyon environment, the mismatch rate in the initial set of corresponding point pairs is usually as high as 70% to 90%. If all corresponding point pairs are directly used for registration transformation, it will lead to serious registration errors.

[0028] To address the high proportion of mismatches in the initial set of corresponding point pairs caused by canyon occlusion, this implementation employs a maximum clique solution method based on invariant measurement maps to identify true interior points. The Euclidean distance between any two points remains unchanged under rigid transformations; this geometric invariance is the key criterion for distinguishing interior and exterior points. Let the initial set of corresponding point pairs contain... The corresponding point pair, the first Each corresponding point pair is denoted as ,in For the points in the sparse point cloud to be registered, In the baseline sparse point cloud and Matching corresponding points. For any two pairs of corresponding points... and If both are correct matches, then the two points on the side to be registered... and The Euclidean distance between the two reference points should be equal to that between the two reference points. and The Euclidean distance between them. Based on the above invariance principle, an invariant measurement graph is constructed: each corresponding point pair in the initial set of corresponding point pairs is taken as a node in the invariant measurement graph, and the graph contains a total of There are 10 nodes; for any two nodes, calculate the Euclidean distance between the points to be registered in the first node and the points to be registered in the second node. And the Euclidean distance between the reference point in the first node and the reference point in the second node. ,in Represent the Euclidean norm; calculate the absolute value of the difference between the two Euclidean distances mentioned above. When the absolute value of the difference is less than the preset distance invariance tolerance In this case, an undirected edge is added between two nodes to indicate that the two corresponding node pairs are compatible. Distance invariance tolerance. The value of needs to take into account the measurement noise level of the point cloud. In this embodiment, the value is 0.1 meters, which is approximately three times the ranging accuracy of the lidar. Traverse all... Each node completes the construction of the invariant measurement graph.

[0029] refer to Figure 1In canyon vegetation point cloud registration scenarios, due to the influence of occlusion and repetitive structures, there are a large number of mismatched outliers in the initial set of corresponding point pairs, with the proportion of outliers typically reaching 70% to 90%. To identify the true inliers from the initial set of corresponding point pairs, this invention constructs an invariant measurement map and searches for the maximum clique on it. Figure 1 The diagram illustrates an invariant measurement graph structure containing 15 nodes, where each node corresponds to a pair of points in the initial set of corresponding points. The nodes in the graph... to These represent 15 corresponding point pairs, where circular nodes represent the actual matching interior points, including... , , , , , and There are 7 nodes in total. Square nodes represent outliers that are mismatched, including... , , , , , , and There are 8 nodes in total. The 2D coordinates of the nodes in the graph are for visualization purposes only and do not represent the actual spatial locations of the corresponding point pairs. The edges of the invariant measurement graph are added according to the distance invariance constraint. For any two nodes, let the corresponding point pair of the first node be denoted as . The corresponding point pair for the second node is ,in and For the points in the sparse point cloud to be registered, and The corresponding point in the baseline sparse point cloud is the one that matches it.

[0030] Calculate the Euclidean distance between the two points on the side to be registered. and the Euclidean distance between the two points on the reference side. If the absolute value of the difference between two distances Less than the preset distance invariance tolerance If two points are matched, an edge is added between them to indicate that the two corresponding point pairs are compatible. Solid lines in the diagram represent edges that satisfy the distance invariance constraint, mainly connecting interior nodes. This is because the distance between truly matched corresponding point pairs remains constant under rigid transformations, thus interior points are highly likely to satisfy the distance invariance constraint. Dashed lines in the diagram represent edges that do not satisfy the distance invariance constraint but happen to satisfy it due to random factors, mainly appearing between exterior nodes or between interior and exterior points. Figure 1It can be observed that the seven interior nodes form a dense network of connections, with each pair of these seven nodes connected by an edge, constituting a complete subgraph, i.e., a clique in graph theory. The connections between exterior nodes are relatively sparse, and the number of connections between exterior and interior nodes is finite. This structural feature indicates that by searching for the maximum clique on the invariant measurement graph, it is possible to effectively identify mutually compatible maximum subsets of interior nodes. Figure 1 Once the invariant measurement graph shown is constructed, it contains 15 nodes and several edges. The maximum clique formed by the internal nodes contains 7 nodes, corresponding to the 7 true matches in the initial set of corresponding point pairs.

[0031] The maximum clique in an invariant measurement graph corresponds to the largest compatible subset of interior vertices in the initial set of corresponding pairs. A clique, a concept in graph theory, refers to a subset of nodes in a graph where every two nodes are connected by an edge; the maximum clique is the clique containing the largest number of nodes. In an invariant measurement graph, all corresponding pairs within a clique satisfy the distance-invariant constraint, thus indicating a high probability that corresponding pairs within a clique are correctly matched. Finding the maximum clique is a combinatorial optimization problem; this implementation uses a branch-and-bound strategy to search for the maximum clique on the invariant measurement graph.

[0032] The specific execution process of branch and bound search is as follows: Initialize the current maximum clique as an empty set, and initialize the candidate node set to all nodes of the invariant measurement graph. Select the node with the highest degree from the candidate node set and add it to the current search path. The node with the highest degree represents the node connected to the most other nodes. Prioritizing the node with the highest degree can improve search efficiency because a node with a high degree is more likely to belong to a large clique. Update the candidate node set to the intersection of the common neighbor nodes of all nodes in the current search path. A common neighbor node is a node that is adjacent to every node in the search path. Only common neighbor nodes can form a clique with the nodes in the search path. If the updated candidate node set is not empty, the above selection and update process is recursively executed to continue expanding the current search path. If the updated candidate node set is empty, it means that the current search path can no longer be expanded. At this time, the current search path is compared with the current maximum clique. If the number of nodes contained in the current search path is greater than the number of nodes contained in the current maximum clique, the current maximum clique is updated to the current search path. Then, backtrack to the previous level search path, select the node with the second highest degree from the candidate node set to continue the search and explore other possible search branches. During the search process, when the sum of the remaining nodes in the candidate node set and the number of nodes already in the current search path is less than the number of nodes in the current largest clique, even if all nodes in the candidate node set are added to the search path, it will still be impossible to surpass the currently found largest clique. Therefore, pruning is performed, abandoning the current branch, and backtracking to the previous level to continue the search. The pruning strategy can significantly reduce invalid searches and improve algorithm efficiency. After the search is completed, the corresponding point pairs of the nodes contained in the current largest clique constitute the interior point set. In actual tests in canyon vegetation scenarios, the initial set of corresponding point pairs usually contains about 5000 corresponding point pairs, and the interior point set identified after solving the largest clique usually contains 200 to 500 corresponding point pairs, with the interior point ratio being about 4% to 10%, which is in stark contrast to the high mismatch rate of the initial matching.

[0033] refer to Figure 2The horizontal axis represents the number of search iterations, ranging from 0 to 24, for a total of 25 iterations. The vertical axis represents the number of nodes in the current search path, ranging from 0 to 9. The curve shows the trend of the current search path length changing with the number of iterations during the search process, with the filled area below the curve visually indicating the growth process of the search path. Branch and bound search starts from an empty set, initializing the current maximum clique as an empty set, and the candidate node set contains all nodes of the invariant measurement graph. In the first iteration, the node with the highest degree is selected from the candidate node set and added to the current search path, increasing the search path length to 1. Subsequently, the candidate node set is updated to the intersection of the common neighbors of all nodes in the current search path. During the second to fourth iterations, the search path length gradually increases to 2 and 3, indicating that the search algorithm continuously adds compatible nodes to the search path. From the fifth iteration onwards, the search path length remains constant in some iterations because the algorithm needs to backtrack when exploring different search branches. When the search path can no longer be expanded, the algorithm compares the current search path with the current maximum clique. If the number of nodes in the current search path is greater than the number of nodes in the current maximum clique, then the current maximum clique is updated.

[0034] The algorithm then backtracks to the previous search path, selecting the node with the second-largest degree from the candidate node set to continue searching other branches. The core optimization of the branch and bound strategy lies in the pruning operation. During the search process, when the sum of the number of remaining nodes in the candidate node set and the number of nodes already in the current search path is less than the number of nodes in the current largest clique, even if all nodes in the candidate node set are added to the search path, it is impossible to obtain a clique larger than the current largest clique. Therefore, the current branch is abandoned and backtracking is initiated. Figure 2 The curve stabilizes after the 17th iteration, with the search path length remaining at 7, indicating that the algorithm has found the largest clique containing 7 nodes. Subsequent iterations, while continuing to explore other branches, failed to find larger cliques. The horizontal reference line in the figure marks the largest clique size as 7. Figure 1 The clique is formed by the seven interior nodes. The entire search process was completed in 25 iterations, and the final determined maximum clique contains seven nodes. The corresponding pairs of these seven nodes form the interior set, which is used for subsequent global coarse registration transformation. The branch and bound strategy significantly reduces invalid search through pruning operations, greatly improving the algorithm efficiency compared to exhaustive search, allowing the maximum clique solution to be completed in an acceptable time.

[0035] The global coarse registration transformation is solved using singular value decomposition based on the interior point set. The 3D coordinates of all points to be registered are extracted from the interior point set and arranged in rows to form a coordinate matrix of the points to be registered. The matrix size is ,in This represents the number of corresponding point pairs in the set of interior points. Similarly, the three-dimensional coordinates of all corresponding reference points are extracted to construct the reference point coordinate matrix. The matrix size is also Calculate the centroid coordinates of the coordinate matrix of the points to be registered. The centroid coordinates of the reference point coordinate matrix ,in and Represent the coordinates of the point to be registered and the reference point, respectively. Line, i.e., the first The three-dimensional coordinates of the registration point and the reference point are aligned using an internal point matrix. The centroid coordinates of the registration point are subtracted from the coordinates of each point in the registration point coordinate matrix to obtain the centroid-free registration point coordinate matrix. ,in For elements all equal to 1 A column vector. Similarly, subtract the centroid coordinates of the reference point from the coordinates of each point in the reference point coordinate matrix to obtain the centroid-free reference point coordinate matrix. The purpose of centroid removal is to separate rotation and translation, simplifying the solution process. The covariance matrix is ​​obtained by calculating the product of the coordinate matrix of the centroid-required point and the transpose of the centroid-required reference point coordinate matrix. The covariance matrix is A square matrix. Perform singular value decomposition on the covariance matrix. ,in for The left singular vector matrix, for The singular value diagonal matrix, for The right singular vector matrix is ​​obtained by multiplying the right singular vector matrix by the transpose of the left singular vector matrix. The translation vector is obtained by subtracting the product of the rotation matrix and the centroid coordinates of the point to be registered from the centroid coordinates of the reference point. Rotation matrix Translation vector Together, they constitute a global coarse registration transformation, which can transform the point cloud to be registered from its original coordinate system to the coordinate system of the reference point cloud.

[0036] Global coarse registration transformation is based on sparse feature points, and its registration accuracy is limited by the distribution and number of feature points, typically resulting in residual errors ranging from several centimeters to tens of centimeters. To further improve registration accuracy, local fine registration is performed using the global coarse registration transformation as the initial value. Local fine registration employs a point-to-surface iterative nearest-point method. Compared to the traditional point-to-point iterative nearest-point method, the point-to-surface method utilizes local surface information, resulting in faster convergence and higher accuracy. The iterative process of the point-to-surface iterative nearest-point method is as follows: First, the point cloud to be registered is transformed to the reference point cloud coordinate system using the current transformation; for each point in the transformed point cloud, the nearest point in the reference point cloud is searched as the corresponding point, and the normal vector of the corresponding point is obtained; an error function is constructed as the sum of the squares of the projected distances along the normal vector direction between the transformed point to be registered and the corresponding reference point; the transformation parameters are updated through linearization and least squares solutions; the above process is repeated until the change in the transformation parameters is less than a preset threshold or the maximum number of iterations is reached. In this embodiment, the maximum number of iterations is set to 50, and the threshold for transformation parameter changes is set to a rotation angle of less than 0.001 degrees and a translation amount of less than 0.001 meters. After local fine registration is completed, a fine registration transformation is obtained. The global coarse registration transformation is concatenated with the fine registration transformation to obtain the final registration transformation from the current point cloud to the reference point cloud. The final registration transformation is used to perform coordinate transformation on all points in the current point cloud to be registered. The above registration process is repeated for the seven point clouds to be registered in the multi-temporal point cloud dataset to obtain a multi-temporal registration point cloud set in a unified coordinate system.

[0037] In one alternative implementation, when the size of the interior point set is small, causing the global coarse registration transformation to be unstable, a random sampling consensus method can be used to verify and optimize the results of the singular value decomposition solution. Specifically, the transformation parameters are solved by randomly sampling the minimum number of corresponding point pairs from the interior point set, the registration error of all interior points under this transformation is statistically analyzed, and after repeated sampling, the transformation with the smallest registration error is selected as the final result.

[0038] The fusion of multi-temporal registration point clouds with corresponding temporal optical images aims to assign color information to each point in the point cloud, facilitating subsequent semantic segmentation and visual analysis. The fusion process is based on camera calibration parameters, which include intrinsic and extrinsic parameters. Intrinsic parameters describe the camera's internal geometry, including focal length, principal point coordinates, and distortion coefficients, obtained before the flight mission using a checkerboard calibration method. Extrinsic parameters describe the rigid transformation relationship between the lidar coordinate system and the camera coordinate system, including rotation matrices and translation vectors, obtained through joint calibration of the lidar and camera. Let the three-dimensional coordinates of a point in the multi-temporal registration point cloud be... (In the lidar coordinate system), firstly, the coordinates are transformed to the camera coordinate system using extrinsic parameters. ,in and These are the rotation matrix and translation vector from the LiDAR coordinate system to the camera coordinate system, respectively. Then, the 3D points in the camera coordinate system are projected onto the image plane using intrinsic parameters. ,in and These are the components of the focal length along the two axes of the image coordinate system. and Principal point coordinates These are the pixel coordinates of the point cloud points projected onto the image. The pixel coordinates are read from the optical image of the corresponding time phase. The red, green, and blue color values ​​of a point are used to assign color values ​​to the point cloud, completing the color assignment for a single point. For point cloud points that fall outside the image boundary or are occluded and cannot obtain valid color values, valid projections are searched from other overlapping images of the same time phase. All points and all time phases in the multi-temporal registration point cloud set are traversed, and the corresponding pixel color value is assigned to each point, obtaining a multi-temporal color point cloud set. The multi-temporal color point cloud set provides rich geometric and spectral information for subsequent ground point filtering, plant segmentation, and leaf recognition.

[0039] After completing multi-temporal point cloud registration and color fusion, the next step is to extract the skeletal structure of the plants from the multi-temporal color point cloud and separate the leaf point cloud, thereby establishing the leaf extension vector field in three-dimensional space. The three-dimensional point cloud of canyon vegetation contains various elements such as ground, rocks, tree trunks, branches, and leaves. It is necessary to gradually peel off various elements through a systematic processing flow to accurately obtain the spatial distribution and temporal changes of the leaves. Ground point filtering is the first step in vegetation point cloud extraction. The canyon terrain is undulating, with ground elevation changes reaching tens of meters. Traditional filtering methods based on fixed elevation thresholds are difficult to adapt to this complex terrain. This implementation uses a cloth simulation filtering method for ground point filtering. This method simulates a soft cloth falling and covering the inverted point cloud surface after the point cloud is inverted. The point corresponding to the final stable position of the cloth is the ground point. Specifically, the color point cloud of the current temporal phase in the multi-temporal color point cloud is inverted vertically, that is, the elevation coordinates of each point are reversed. A virtual cloth composed of a regular grid is placed above the inverted point cloud, with the grid node spacing set to 1 meter. The fabric gradually falls under the influence of gravity. Each grid node checks for collisions with the inverted point cloud during its descent; if a collision occurs, the node stops falling. Simultaneously, elastic constraints exist between adjacent grid nodes to maintain a certain rigidity of the fabric. This falling and constraint adjustment process is iteratively executed until the fabric shape stabilizes. After stabilization, the vertical distance between each point in the original point cloud and the fabric surface is calculated. Points with a distance less than 0.3 meters are marked as ground points and filtered out; the remaining points constitute the vegetation point cloud. The fabric rigidity parameter affects the filtering effect; a rigidity parameter value of 2 can better adapt to moderate undulations in canyon terrain. In an optional implementation, for areas with extremely complex terrain, a progressive morphological filtering method or a multi-scale surface fitting method can be used to filter out ground points.

[0040] The vegetation point cloud contains all plants within the target area, and connectivity clustering is needed to separate different plants into independent point cloud clusters. Connectivity clustering determines the connectivity between points based on Euclidean distance, grouping spatially close points into the same cluster. For any point in the vegetation point cloud, all neighboring points with an Euclidean distance less than the clustering distance threshold are searched in the point cloud, and these neighboring points are grouped into the same cluster as the current point. The above neighboring point search and merging operation is recursively performed on each point newly added to the cluster until no more new points are added to the cluster. Then, a point is selected from the points that have not yet been assigned to any cluster as the seed point of the new cluster, and the above clustering process is repeated. The choice of the clustering distance threshold needs to balance the completeness and accuracy of plant separation. If the threshold is too small, the same plant will be divided into multiple clusters; if the threshold is too large, adjacent plants will be merged into one cluster. In this embodiment, the clustering distance threshold is set to 0.8 meters, which is slightly larger than the typical spacing between adjacent branches and leaves in canyon vegetation. After clustering, small clusters with fewer than 500 points are filtered out, as these clusters typically correspond to shrub fragments or noise points. Each retained point cloud cluster corresponds to an individual plant, thus obtaining several individual plant point clouds.

[0041] refer to Figure 4 This diagram, constructed in a three-dimensional Cartesian coordinate system, visually reflects the geometric determination process of separating leaf point clouds from complex canopy point clouds. Figure 4 In the diagram, the thick black lines represent the plant's skeletal structure constructed using the minimum spanning tree algorithm. This skeletal structure exhibits a clear hierarchical topology, including a vertically upward-growing main trunk and primary and secondary branches extending outwards from the main trunk. These black skeletal line segments represent... Figure 3 The candidate skeleton points generated in the process are obtained by performing topological connections and path smoothing. They constitute the three-dimensional geometric support network of the plant and serve as a spatial reference benchmark for subsequent distance measurement and point cloud classification. Figure 4 The core logic lies in demonstrating how to distinguish between branch points and leaf points based on the Euclidean distance from a point to the skeleton. The diagram uses a semi-transparent blue grid cylinder to visually construct the classification threshold boundary. This blue grid cylinder surrounds the black skeleton lines, and its radius corresponds to the preset skeleton distance threshold. In this embodiment, the threshold is set to 0.25 meters. The physical significance of this threshold lies in defining the effective influence range of the branches. For any given original lidar point in space... The algorithm calculates its connection to the set of black skeleton line segments. Shortest Euclidean distance ,in Let be any point on the skeleton line segment. Figure 4The classification results are visually illustrated using color coding of the points: points distributed inside or on the surface of the blue grid cylinder are marked as brown scattered points, which represent branch points. This is due to their distance... Points smaller than the threshold of 0.25 meters are identified as attachments to the surface of the trunk or branches and are considered to be non-leaf tissue.

[0042] Conversely, points distributed in the region outside the blue grid cylinder are marked as green dots, representing successfully separated blade points. Figure 4 In the image, the green dot clouds exhibit a distribution pattern that spreads outward from the framework, with the shortest Euclidean distance from the framework all greater than or equal to 0.25 meters. This distance-field-based segmentation strategy effectively utilizes botanical morphological characteristics, namely that leaves typically grow at the ends of branches or extend a certain distance beyond the branch via petioles. Figure 4 Clearly, as the spatial location moves away from the center of the black skeleton, the brown dots gradually transition and disappear, replaced by numerous green leaf dots, with the blue threshold boundary serving as the dividing line, achieving high-precision semantic segmentation. Furthermore, Figure 4 It also demonstrates the complexity of the skeletal branching, including the trunk (coordinates). (from 0 to 2 meters) and bifurcation structure (coordinates) (Within a range greater than 2 meters), the blue threshold boundary adaptively extends along the skeleton, ensuring consistency of the segmentation standard throughout the entire plant space. In this way, Figure 4 This study visually verifies the effectiveness of transforming unstructured massive point cloud data into structured data with clear semantic labels (branches / leaves). The separated green leaf point cloud set will be directly used as the basic input data for the subsequent construction of the leaf extension vector field and spatiotemporal dynamic model.

[0043] For each individual plant point cloud, an iterative shrinkage process is performed to extract the plant's skeletal structure. The plant skeleton is a simplified representation of the plant's three-dimensional morphology, composed of a series of connected line segments, reflecting the topological connections from the trunk to branches at various levels. Skeleton extraction plays a crucial role in distinguishing branches from leaves, as branch points are distributed near the skeleton, while leaf points are distributed at the far end of the skeleton. The iterative shrinkage method gradually shrinks the point cloud towards its local central axis, causing points originally scattered on the surface of branches and leaves to gradually converge near the centerline of the branches, ultimately forming an approximately one-dimensional skeletal structure.

[0044] In the initialization phase of the iterative shrinkage process, the current independent plant point cloud is copied into a shrunken point cloud. Points in the shrunken point cloud will continuously move in subsequent iterations. An initial shrinkage step size is set for each point in the shrunken point cloud, with a value of 0.05 meters, approximately 2 to 3 times the typical leaf thickness. The shrinkage step size determines the distance the point moves in each iteration; a step size that is too large may lead to instability in the shrinkage process, while a step size that is too small will require more iterations to converge.

[0045] The core calculation process of iterative shrinkage is as follows: For each current point in the shrinking point cloud... Search within the contracted point cloud for points whose Euclidean distance to the current point is less than a preset neighborhood radius. All points in the set of points form the neighborhood of the current point. The preset neighborhood radius is used. The value is 0.3 meters, and this value needs to be greater than the shrinkage step size to ensure that the neighborhood contains a sufficient number of points for direction calculation. Let the neighborhood point set contain... The nth neighboring point, the th Let the neighboring points be denoted as For each neighboring point in the neighborhood point set... Calculate the current point Point to neighboring points Direction vector Normalize the direction vector to obtain the unit direction vector. ,in The Euclidean norm of the direction vector, i.e., the current point. With neighboring points The Euclidean distance between the points is calculated. Simultaneously, the reciprocal of the Euclidean distance between the current point and its neighbors is calculated. The contribution of the unit direction vector is considered. Closer neighboring points have a greater influence on the contraction direction. This design causes the point cloud to contract towards locally denser regions, and the branch centerline happens to be the region with the highest point cloud density. The unit direction vectors corresponding to all neighboring points in the neighborhood set are weighted and summed according to their respective contributions to obtain the accumulated vector. Normalize the accumulated vector to obtain the contraction direction vector of the current point. The contraction direction vector points towards the weighted centroid direction of the neighboring points of the current point, i.e., towards the center of the local point cloud. The current point... vector along the contraction direction Move the current contraction step size Get the new position of the current point .

[0046] After updating the position at each point, it is necessary to determine whether the point has converged. Calculate the new position of the current point. Original position of the current point Euclidean distance between If the moving distance is less than a preset convergence threshold, the current point is marked as converged. The preset convergence threshold is 0.005 meters. When the moving distance of a point is less than this threshold, it is considered that the point has reached the local central axis position and no further shrinkage is needed. After performing the above shrinkage movement operation on all points in the shrinking point cloud, the proportion of converged points in the shrinking point cloud to the total number of points is counted. When this proportion exceeds a preset proportion threshold, the iteration is terminated. The preset proportion threshold is 0.95, that is, the iteration stops when more than 95% of the points have converged. If the convergence proportion does not reach the threshold, the shrinkage step size is reduced to half of the current value before continuing to the next round of iteration. The gradual reduction of the shrinkage step size makes the iteration process go from coarse to fine, with large step sizes for rapid shrinkage in the early stage and small step sizes for fine adjustment in the later stage. Usually, after 15 to 25 rounds of iteration, the shrinking point cloud converges to the skeleton shape.

[0047] After the iteration terminates, the points in the shrunken point cloud are highly clustered near the skeleton lines, but there are many overlapping points. To remove overlapping points, the shrunken point cloud is downsampled using a voxel grid downsampling method. The voxel side length is set to 0.02 meters, and the points within each voxel are replaced by their centroids. After downsampling, a skeleton candidate point set is obtained, which typically contains hundreds to thousands of points distributed along the plant skeleton lines.

[0048] A minimum spanning tree is constructed on the candidate skeleton points to form the topological structure of the plant skeleton. The minimum spanning tree is an acyclic connected graph with the minimum total edge length connecting all candidate skeleton points, naturally reflecting the branching structure of the plant from the trunk to each level of branches. The minimum spanning tree is constructed using Prim's algorithm: Initially, the point with the lowest elevation in the candidate skeleton point set is selected as the root node, which typically corresponds to the base of the trunk; the root node is added to the visited set; the Euclidean distance between all points in the visited set and all points in the unvisited set is calculated; the pair of points with the smallest distance is selected, the unvisited point is added to the visited set, and an edge is added between the two points; this process is repeated until all candidate skeleton points have been visited. The edges of the minimum spanning tree are line segments connecting two candidate skeleton points, and the length of the edge is the Euclidean distance between the two candidate skeleton points.

[0049] After constructing the minimum spanning tree, topological analysis is performed to identify the structural elements of the skeleton. In the minimum spanning tree, nodes with a degree of one are called endpoints, where degree refers to the number of edges connected to that node. An endpoint is connected to only one edge and usually corresponds to the end of a branch. Nodes with a degree greater than two are called branching points, which are connected to three or more edges and correspond to the branching positions of branches. Nodes with a degree equal to two are ordinary nodes, located in the middle of a branch. The edge sequences connecting adjacent branching points or between a branching point and an endpoint are merged into skeleton branches, each corresponding to a continuous branch segment. All skeleton branches constitute the plant skeleton. Among the skeleton branches, the longest cumulative length sequence of skeleton branches originating from the root node is identified as the main trunk skeleton, corresponding to the central trunk of the plant. The remaining skeleton branches are branch skeletons, corresponding to the various levels of lateral branches branching from the main trunk.

[0050] The individual plant point cloud is separated into leaf point clouds based on the distance from each point to the plant skeleton. For each original point in the individual plant point cloud, the shortest Euclidean distance from the original point to the plant skeleton is calculated. The plant skeleton consists of multiple skeleton branches, and each skeleton branch consists of multiple edges. The distance from the original point to each edge is calculated, and the minimum value is taken as the shortest Euclidean distance from the original point to the plant skeleton. The method for calculating the distance from a point to a line segment is as follows: Let the two endpoints of the line segment be... and The point whose distance needs to be calculated is First calculate the point Online segment The projection position on the line segment is determined by the Euclidean distance: if the projection position falls inside the line segment, the distance is the Euclidean distance from the point to the projection position; if the projection position falls outside the line segment, the distance is the Euclidean distance from the point to the nearest endpoint of the line segment. When the shortest Euclidean distance is greater than or equal to a preset skeleton distance threshold, the original point is marked as a leaf point; when the shortest Euclidean distance is less than the preset skeleton distance threshold, the original point is marked as a branch point. The preset skeleton distance threshold is 0.15 meters, approximately 1.5 times the average diameter of the branch, which allows points attached to the branch surface to be classified as branches, while correctly separating leaf points far from the skeleton. All leaf points constitute a leaf point cloud.

[0051] In one alternative implementation, the skeleton distance threshold can be adaptively adjusted according to the level of the skeleton branches. The distance threshold for the main trunk skeleton is set to 0.25 meters, the distance threshold for the first-level branch skeleton is set to 0.15 meters, and the distance threshold for the second-level and higher-level branch skeleton is set to 0.08 meters. This adaptive threshold strategy can more accurately adapt to the geometric characteristics of branches of different thicknesses.

[0052] The above-described processes of ground point filtering, connectivity clustering, skeleton extraction, and leaf separation are repeated for all phases in the multi-temporal color point cloud set to obtain a multi-temporal plant skeleton set and a multi-temporal leaf point cloud set. The multi-temporal plant skeleton set contains plant skeleton data from 8 phases, and the multi-temporal leaf point cloud set contains leaf point cloud data from 8 phases.

[0053] refer to Figure 3 This image visually presents the contraction and evolution process of a three-dimensional point cloud on a cross-section perpendicular to the branch axis, using a two-dimensional planar projection. Figure 3 In the initial point cloud of individual plants, light green scattered dots are represented on the outermost layer. These dots simulate the original lidar echo data surrounding the plant branches in space. Due to measurement noise and the roughness of the plant surface, the initial point cloud exhibits a ring-shaped structure with a certain thickness centered on the branch center. Its radial distribution follows a Gaussian distribution with a mean of 1.0 m and a standard deviation of 0.1 m. Figure 3 The core geometric elements of the iterative shrinkage algorithm are clearly marked, with the blue dashed circles representing the current point. The defined local neighborhood range, its radius The radius is set to 0.3 meters. The choice of this local neighborhood radius is crucial; it must be larger than the step size of a single contraction and cover enough neighborhood points to ensure robustness of the direction calculation. Figure 3 The gray arrows illustrate the dynamic mechanism of point cloud contraction during the shown contraction process. For any current point located on the annular band... The algorithm first searches for all neighboring points that fall within the range of the blue dashed circle. Then, calculate from the current point Pointing to each neighboring point unit direction vector To direct the contraction direction towards the local central axis location with higher point cloud density, a weighting mechanism based on the reciprocal of the Euclidean distance is introduced, i.e., weighting. . Figure 3 The gray arrow in the diagram represents the final contraction direction vector obtained through weighted summation. This vector The mathematical expression is This geometric definition ensures that, regardless of the initial point's location on the branch's cross-section, the direction of the force always tends to point towards the centroid of the local neighborhood, thus driving the point cloud to converge towards the geometric center.

[0054] Figure 3 The intermediate states of the contraction process are further illustrated, represented by orange dots located in the inner circle. These orange dots are the contraction direction vectors indicated by the gray arrows from the initial light green dots. The contraction step size was shifted. The position after. Initial contraction step size. The distance was set to 0.05 meters. After multiple iterations of position updates, the point cloud gradually broke free from the constraints of the original surface and moved closer to the central axis. Figure 3 The central position shows the final contraction state, with red cross-shaped markers representing the skeleton candidate point set. These red markers are closely clustered near the origin of the cross-section, indicating that the algorithm has successfully compressed the scattered surface point cloud into an approximately one-dimensional linear structure. At this point, the distance between points is less than the preset convergence threshold of 0.005 meters, satisfying the convergence condition. Figure 3 By comparing the spatial distribution of the initial point (light green), intermediate point (orange), and final skeleton candidate point (red), the effectiveness of the Laplace contraction principle in extracting vegetation skeletons is intuitively verified. That is, while keeping the point cloud topology unchanged, the three-dimensional volume data is collapsed into a one-dimensional geometric skeleton through iterative smoothing, providing a precise geometric benchmark for subsequent topological connections and leaf separation.

[0055] Single-leaf segmentation is performed on multi-temporal leaf point clouds, dividing the leaf point cloud of each temporal phase into several single-leaf point clouds. A leaf point cloud is a set of points from all leaves of the plant; the goal of single-leaf segmentation is to group points belonging to the same leaf together. Since leaves typically present a flat, thin morphology with abrupt curvature changes between adjacent leaves, a curvature-based region segmentation method is used for single-leaf segmentation.

[0056] The curvature-based region segmentation process is as follows: First, the local curvature is calculated for each point in the leaf point cloud. For the current point... Search by A local neighborhood point set is formed by all points within a spherical neighborhood centered at a radius of 0.05 meters. Principal component analysis is performed on the three-dimensional coordinates of the local neighborhood point set to obtain three eigenvalues. ,in , , These are the first, second, and third eigenvalues ​​of the covariance matrix, respectively. The smallest eigenvalue... The corresponding eigenvector is the direction of the normal vector of the local surface. Local curvature is defined as... This ratio reflects the flatness of the local point distribution. The smaller the curvature value, the flatter the local area is, and the larger the curvature value, the more severe the local curvature is.

[0057] After calculating the local curvature of all points, a region growing method is used for segmentation. The point with the smallest curvature in the leaf point cloud is selected as the seed point; this point is located in a flat area of ​​the leaf. The seed point is added to the current region, and its neighbors are searched. If the curvature of a neighboring point is less than a curvature threshold and the angle between the neighboring point's normal vector and the seed point's normal vector is less than an angle threshold, then the neighboring point is added to the current region. The curvature threshold is set to 0.1, and the angle threshold is set to 15 degrees. Neighbor search and judgment are recursively performed for each point in the newly added region until the current region no longer expands. Then, the point with the smallest curvature from the remaining unassigned points is selected as the seed point for a new region, and the region growing process is repeated. After segmentation, each region corresponds to a single leaf point cloud. Small regions with fewer than 50 points are filtered out; these small regions usually correspond to leaf debris or noise. This yields several single leaf point clouds.

[0058] For each individual leaf point cloud, its geometric features are calculated for subsequent time-series matching. The 3D coordinates of the centroid of each individual leaf point cloud are calculated. This is the arithmetic mean of the three-dimensional coordinates of all points in the point cloud of a single leaf, i.e. ,in This represents the number of points in a single leaf point cloud. For the first The three-dimensional coordinates of each point are calculated. The three-dimensional coordinates of the centroid reflect the position of the blade in three-dimensional space. Simultaneously, the principal direction vector of each individual blade point cloud is calculated. Principal component analysis is performed on the three-dimensional coordinates of the individual blade point clouds; the eigenvector corresponding to the largest eigenvalue is the principal direction vector. This reflects the main extension direction of the blade. To ensure the consistency of the main direction vector, the perpendicular component of the main direction vector is defined. If the calculated principal direction vector is positive, then the entire vector is inverted if the vertical component is negative.

[0059] Matching is performed on the point cloud sets of single blades from two adjacent time phases to establish the correspondence between the same blade in different time phases. Let the point cloud set of single blades in the previous time phase contain... Each single blade, the point cloud set of the single blade in the later time phase contains A single blade. The construction scale is... The cost matrix, the cost matrix of the first... Line number The column element represents the first element of the previous time phase. The first single blade and the first in the next time phase The cost of matching a single blade. The matching cost comprehensively considers the Euclidean distance between the three-dimensional coordinates of the centroids and the angle between the principal direction vectors. Let the previous time phase be... The three-dimensional coordinates of the centroid of a single blade are as follows: The main direction vector is ; Later, the prime minister The three-dimensional coordinates of the centroid of a single blade are as follows: The main direction vector is Distance between centroids The Euclidean distance between the two centroids in three-dimensional coordinates. The angle between the principal directions. The angle between the two principal direction vectors is taken as the absolute value because the positive or negative direction of the principal direction vectors does not affect the actual orientation of the blade. Matching cost This is the weighted sum of the distance between the centroids and the angles between the principal directions, where... This is the angle penalty coefficient, with a value of 0.5 meters per radian, used to convert angular dimensions to distance dimensions. If the centroid distance... If the maximum matching distance threshold of 2 meters is exceeded, the cost is set to infinity, indicating that the two single blades cannot be the same blade.

[0060] The Hungarian algorithm is used to find the optimal matching on the cost matrix. The Hungarian algorithm is a classic algorithm for finding minimum-weighted matching in a bipartite graph, capable of finding a one-to-one matching scheme that minimizes the total matching cost in polynomial time. The execution result of the Hungarian algorithm is a set of matching pairs, each containing a single-leaf index from the previous time phase and a single-leaf index from the next time phase. Matching pairs with infinite costs are discarded; these pairs correspond to fallen leaves that existed in the previous time phase but disappeared in the next time phase, or newly formed leaves that appeared in the next time phase but did not exist in the previous time phase. The retained matching pairs form the correspondence of the same leaf between adjacent time phases.

[0061] For each successfully matched blade pair, calculate its extension velocity vector between adjacent time phases. Let the three-dimensional coordinates of the centroid of the same blade in the previous time phase be... The three-dimensional coordinates of the centroid in the later time phase are The difference between the three-dimensional coordinates of the centroid in the subsequent time phase and the three-dimensional coordinates of the centroid in the previous time phase is calculated to obtain the blade's extension displacement vector. ,in , , These represent the components of the extension displacement vector in the east, north, and vertical directions, respectively. Let the time interval between the two time phases be . The time interval is measured in days. Dividing the extension displacement vector by the time interval yields the blade extension velocity vector. The unit is meters per day. The extension velocity vector describes the speed and direction of the leaf's center of mass movement in three-dimensional space, reflecting the leaf's growth and extension state.

[0062] The elevation differences in canyon terrain affect the measurement of the extension velocity vector. On steep slopes, even if the actual extension distance of the blade along the slope direction is the same, the components projected onto the horizontal and vertical directions will differ due to variations in slope, leading to systematic biases in the measurement results. To eliminate measurement errors caused by elevation differences, terrain correction is required for the extension velocity vector.

[0063] The specific process of terrain correction is as follows: Based on the three-dimensional coordinates of the plant skeletons at each time phase in the multi-temporal plant skeleton set, the local terrain slope direction at the location of the three-dimensional coordinates of the centroid of each leaf is calculated. Within the neighborhood of the current leaf centroid location, ground points extracted in the ground point filtering step are collected, with the neighborhood set as a cylindrical region with a radius of 5 meters centered on the leaf centroid. Plane fitting is performed on the ground points within the neighborhood, and the least squares method is used to fit the plane equation. ,in , , These are the components of the plane normal vector. This is a constant term. It normalizes the plane normal vector to a unit normal vector. Define the perpendicular component of the normal vector A positive value means the normal vector points upwards. The direction of the local terrain slope is defined as the direction of the projection of the normal vector onto the horizontal plane, i.e. Slope direction vector The projection of the maximum downward direction of the slope onto the horizontal plane.

[0064] The extension velocity vector of each blade is decomposed into a component along the local terrain slope direction and a component perpendicular to the local terrain slope direction. Let the extension velocity vector be... The local terrain slope direction is Component along the local topographic slope direction ,in This is the dot product of the extension velocity vector and the slope direction vector. The component perpendicular to the local terrain slope direction... The component along the slope direction is significantly affected by topographic relief, including systematic biases caused by the slope; the component perpendicular to the slope direction is less affected by topography and better reflects the true extension state of the leaves. The component perpendicular to the local topographic slope direction should be retained. As the blade extension velocity vector after terrain correction.

[0065] In one alternative implementation, if the number of local ground points is insufficient for reliable plane fitting, the local terrain slope direction can be obtained by interpolation using a digital elevation model. A digital elevation model with a resolution of 1 meter is generated from the ground points gathered from the multi-temporal registration point cloud. The ground elevation at the blade centroid location is obtained through bilinear interpolation, and the local slope direction is calculated by combining the elevations of neighboring grid points.

[0066] The topographically corrected leaf extension velocity vectors of all leaves are spatially organized according to the three-dimensional coordinates of the leaf's centroid to construct a three-dimensional leaf extension vector field. A three-dimensional extension vector field is a spatial field data structure where each location corresponds to a vector value, representing the leaf's extension velocity at that location. Specifically, the three-dimensional space of the plant is divided into a cubic grid with sides of 0.5 meters. Each grid cell stores the average value of all leaf extension velocity vectors falling within that grid. If a grid cell contains no leaves, its extension velocity vector is zero. This yields the three-dimensional leaf extension vector field for a single observation period.

[0067] The above-mentioned single-leaf segmentation, time-series matching, extension velocity vector calculation, and terrain correction processes were repeated for all adjacent observation periods within the growing season to obtain seven sets of three-dimensional leaf extension vector fields at different time intervals. By integrating the three-dimensional leaf extension vector fields from all observation phases, a spatiotemporal dynamic model of canyon vegetation leaf extension covering the growing season was obtained. This model is organized using a four-dimensional data structure, with three spatial dimensions corresponding to east, north, and vertical coordinates, and one temporal dimension corresponding to the observation date. This model can demonstrate the spatial distribution changes and temporal evolution of leaf extension rates from budding to full leaf stage, providing quantitative three-dimensional spatiotemporal data support for canyon vegetation growth monitoring and ecological research.

[0068] In one alternative implementation, the spatiotemporal extension dynamic model can be further extended to a comprehensive model incorporating multiple growth indicators such as the rate of change of leaf area and the rate of change of leaf inclination angle. The leaf area is estimated using the convex hull projection area of ​​a single leaf point cloud, and the leaf inclination angle is calculated using the angle between the principal direction vector and the horizontal plane. These extended indicators, together with the extension velocity vector, constitute a multi-dimensional description of leaf growth characteristics.

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

Claims

1. A three-dimensional reconstruction method for the leaf extension process of canyon vegetation based on image fusion, characterized in that, The method includes: Step 1: Integrate LiDAR and optical camera on the UAV platform, achieve time synchronization through hardware trigger signal, and repeatedly scan the target canyon vegetation area on multiple observation dates during the growing season to obtain multi-temporal point cloud dataset and multi-temporal optical image set; Step 2: The earliest point cloud data in the multi-temporal point cloud dataset is designated as the reference point cloud, and the remaining point cloud data from each temporal phase are designated as point clouds to be registered. Feature descriptors are extracted and matched from the reference and point clouds to be registered, forming an initial set of corresponding point pairs. Each corresponding point pair in the initial set is used as a node in the invariant measurement map, and the invariant measurement map is constructed based on distance invariance constraints. A branch-and-bound strategy is used on the invariant measurement map to search for the maximum clique to identify the set of interior points. Based on the set of interior points, a global coarse registration transformation is solved through singular value decomposition, followed by local fine registration to obtain a multi-temporal registration point cloud set. The multi-temporal registration point cloud set is fused with the corresponding temporal optical images to obtain a multi-temporal color point cloud set. Step 3: Perform ground point filtering and connectivity clustering on the multi-temporal color point cloud set to obtain independent plant point clouds; perform iterative shrinkage processing on the independent plant point clouds to obtain a skeleton candidate point set; construct a minimum spanning tree on the skeleton candidate point set to form the plant skeleton; separate the independent plant point clouds into leaf point clouds based on the distance of each point to the plant skeleton to obtain a multi-temporal leaf point cloud set. Step 4: Perform single-leaf segmentation and time-span matching on the multi-temporal leaf point cloud, calculate the extension velocity vector of each leaf, and construct a dynamic model of the spatiotemporal extension of canyon vegetation leaves after terrain correction.

2. The method according to claim 1, characterized in that, Feature descriptor extraction and matching are performed on the reference point cloud and the point cloud to be registered to form an initial set of corresponding point pairs. This includes: uniformly downsampling the reference point cloud and the current point cloud to be registered to obtain the reference sparse point cloud and the sparse point cloud to be registered; calculating the fast point feature histogram descriptor for each point in the reference sparse point cloud and the sparse point cloud to be registered; comparing the fast point feature histogram descriptor of each point in the sparse point cloud to be registered with the fast point feature histogram descriptor of all points in the reference sparse point cloud using Euclidean distance, and selecting the reference point with the smallest descriptor distance for each point to be registered as the corresponding point to form an initial set of corresponding point pairs.

3. The method according to claim 2, characterized in that, Constructing an invariant measurement graph based on distance invariance constraints includes: for any two nodes, calculating the Euclidean distance between the registration point in the first node and the registration point in the second node, and the Euclidean distance between the reference point in the first node and the reference point in the second node; calculating the absolute value of the difference between the two Euclidean distances; when the absolute value of the difference is less than the preset distance invariance tolerance, adding an edge between the two nodes; and traversing all node pairs to complete the construction of the invariant measurement graph.

4. The method according to claim 3, characterized in that, A branch-and-bound strategy is used on the invariant measurement graph to search for the maximum clique to identify the set of interior points. This includes: initializing the current maximum clique as an empty set and initializing the candidate node set to all nodes in the invariant measurement graph; selecting the node with the highest degree from the candidate node set and adding it to the current search path; updating the candidate node set to the intersection of the common neighbors of all nodes in the current search path; if the updated candidate node set is not empty, recursively executing the above selection and update process; if the updated candidate node set is empty and the number of nodes in the current search path is greater than the number of nodes in the current maximum clique, then updating the current maximum clique to the current search path; backtracking and selecting the node with the second highest degree from the candidate node set to continue the search; when the sum of the remaining number of nodes in the candidate node set and the number of nodes already in the current search path is less than the number of nodes in the current maximum clique, pruning and abandoning the current branch; after the search is completed, the corresponding point pairs of the nodes contained in the current maximum clique constitute the set of interior points.

5. The method according to claim 4, characterized in that, The global coarse registration transformation based on the interior point set is solved by singular value decomposition, including: extracting the 3D coordinates of all points to be registered from the interior point set to form the coordinate matrix of the points to be registered, and extracting the 3D coordinates of all corresponding reference points to form the coordinate matrix of the reference points; calculating the centroid coordinates of the coordinate matrices of the points to be registered and the reference points respectively, and subtracting the corresponding centroid coordinates from the coordinates of each point in each coordinate matrix to obtain the centroid-free coordinate matrix; calculating the product of the transpose of the centroid-free coordinate matrix of the points to be registered and the centroid-free reference point coordinate matrix to obtain the covariance matrix; performing singular value decomposition on the covariance matrix to obtain the left singular vector matrix and the right singular vector matrix; calculating the product of the right singular vector matrix and the transpose of the left singular vector matrix to obtain the rotation matrix; subtracting the product of the rotation matrix and the centroid coordinates of the points to be registered from the centroid coordinates of the reference points to obtain the translation vector; the rotation matrix and the translation vector together constitute the global coarse registration transformation.

6. The method according to claim 5, characterized in that, Next, local fine registration is performed to obtain a multi-temporal registration point cloud, including: using the global coarse registration transformation as the initial value, local fine registration is performed using the point-to-surface iterative nearest point method to obtain the fine registration transformation; the global coarse registration transformation and the fine registration transformation are concatenated to obtain the final registration transformation from the current point cloud to the reference point cloud; the above registration process is repeated for all point clouds to be registered in the multi-temporal point cloud dataset to obtain a multi-temporal registration point cloud in a unified coordinate system; based on the camera calibration parameters, each point in the multi-temporal registration point cloud is projected onto the optical image of the corresponding temporal phase, and each point is assigned a color value of the corresponding pixel to obtain a multi-temporal color point cloud.

7. The method according to claim 1, characterized in that, An iterative shrinkage process is performed on the point cloud of an independent plant to obtain a candidate point set for the skeleton. This process includes: copying the current point cloud of an independent plant into a shrunken point cloud, and setting an initial shrinkage step size for each point in the shrunken point cloud; for each current point in the shrunken point cloud, searching for all points whose Euclidean distance to the current point is less than a preset neighborhood radius to form a neighborhood point set; for each neighborhood point in the neighborhood point set, calculating the unit direction vector from the current point to the neighborhood point, and calculating the reciprocal of the Euclidean distance between the current point and the neighborhood point as the contribution of the unit direction vector; summing and normalizing the unit direction vectors corresponding to all neighborhood points in the neighborhood point set according to their respective contributions to obtain the shrinkage direction vector of the current point; moving the current point along the shrinkage direction vector by the current shrinkage step size; determining whether the moving distance is less than a preset convergence threshold, and marking it as converged if it is less; terminating the iteration when the proportion of the number of converged points to the total number of points exceeds a preset proportion threshold, otherwise reducing the shrinkage step size to half of the current value and continuing the iteration; after the iteration terminates, downsampling the shrunken point cloud to remove overlapping points to obtain a candidate point set for the skeleton.

8. The method according to claim 7, characterized in that, The process of constructing a minimum spanning tree on the candidate skeleton point set to form the plant skeleton includes: constructing a minimum spanning tree on the candidate skeleton point set, wherein the edges of the minimum spanning tree are line segments connecting two candidate skeleton points, and the length of the edge is the Euclidean distance between the two candidate skeleton points; in the minimum spanning tree, nodes with a recognition degree of one are endpoints, and nodes with a recognition degree greater than two are branch points; merging the edge sequences connecting adjacent branch points or between a branch point and an endpoint into skeleton branches, and all skeleton branches constitute the plant skeleton.

9. The method according to claim 8, characterized in that, Based on the distance from each point to the plant skeleton, the independent plant point cloud is separated into leaf point clouds, including: for each original point in the independent plant point cloud, the shortest Euclidean distance from the original point to the plant skeleton is calculated; when the shortest Euclidean distance is greater than or equal to a preset skeleton distance threshold, the original point is marked as a leaf point; all leaf points constitute a leaf point cloud; the skeleton extraction and leaf separation process is repeated for all phases in the multi-temporal color point cloud set to obtain a multi-temporal plant skeleton set and a multi-temporal leaf point cloud set.

10. The method according to claim 9, characterized in that, A dynamic model of the spatiotemporal extension of canyon vegetation leaves is constructed by segmenting and matching single leaves across time periods from multi-temporal leaf point clouds, calculating the extension velocity vector of each leaf, and performing terrain correction. This includes: performing curvature-based region segmentation on the current-temporal leaf point cloud in the multi-temporal leaf point cloud set to obtain several single-leaf point clouds; calculating the centroid 3D coordinates and principal direction vector of each single-leaf point cloud; performing matching on single-leaf point cloud sets of two adjacent time periods, using the Euclidean distance between the centroid 3D coordinates and the angle between the principal direction vectors as matching costs, and employing the Hungarian algorithm to obtain the correspondence of the same leaf between adjacent time periods; for each successfully matched leaf pair, calculating the difference between the centroid 3D coordinates of the later time period and the centroid 3D coordinates of the previous time period to obtain the extension velocity vector. The displacement vector is obtained by dividing the extension displacement vector by the time interval between two time phases to obtain the extension velocity vector. Based on the three-dimensional coordinates of the plant skeletons at each time phase in the multi-time phase plant skeleton set, the local terrain slope direction of the three-dimensional coordinates of the centroid of each leaf is calculated. The extension velocity vector is decomposed into a component along the local terrain slope direction and a component perpendicular to the local terrain slope direction. The component perpendicular to the local terrain slope direction is retained as the terrain-corrected leaf extension velocity vector. The terrain-corrected leaf extension velocity vectors of all leaves are spatially organized according to the three-dimensional coordinates of the centroid of the leaves to construct a three-dimensional extension vector field of plant leaves. The three-dimensional extension vector fields of plant leaves from all observation time phases are integrated to obtain a spatiotemporal extension dynamic model of canyon vegetation leaves covering the growing season.

Images