Canopy morphology 3D feature analysis method based on lightweight UAV and point cloud singular value decomposition

Abstract

A canopy morphology 3D feature analysis method based on lightweight drones and point cloud singular value decomposition belongs to the field of plant morphology analysis. The method is as follows: a drone-mounted RGB camera is combined with the SfM‑MVS algorithm to perform high-precision three-dimensional reconstruction of the soybean canopy structure. At the same time, a label-based method introduces the singular value decomposition (SVD) method for the point cloud coordinate matrix, and 9 canopy morphology 3D features at multiple scales are established as a new interpretable feature at the crop plot level. The results show that the global canopy 3D features are significantly correlated with soybean yield, lodging, canopy height, and vegetation index (p <= 0.05). 3D feature extraction based on local range can extract more detailed structural information of the canopy surface, which helps to characterize the differences in more growth stages of soybeans.

Description

Technical Field

[0001] The present invention belongs to the field of plant morphological analysis, and specifically relates to a canopy morphological 3D feature analysis method based on a lightweight unmanned aerial vehicle and point cloud singular value decomposition. Background Art

[0002] Crop breeding often requires surveying hundreds or even thousands of breeding materials. Traditional field surveys are time-consuming, labor-intensive, and costly, while also limiting the traits they cover. Unmanned aerial vehicle (UAV) platforms, with their high throughput and adaptability, are increasingly being used to study soybean (Glycine max (L.) Merr.) breeding phenotyping. The rapid development of UAV proximity sensing platforms and information technology has enabled effective crop phenotyping and aided breeding decision-making (Zhou et al., 2022; Ren et al., 2023). By equipping with sensors such as RGB, multispectral, hyperspectral, and lidar, phenotypic characteristics such as vegetation index, canopy cover, texture, and grain can be extracted (Mardanisamani et al., 2019; Yue et al., 2019). These characteristics can be used to estimate biomass (Ogawa et al., 2021; Che et al., 2022), yield (Fei et al., 2023), nitrogen content (Ding et al., 2022), and chlorophyll content (Hu et al., 2010), as well as for early diagnosis of nutrients (Wiwart et al., 2009) and identification of diseases at the plot scale (Cao et al., 2018).

[0003] Compared to two-dimensional spectral data, which suffers from resolution limitations and occlusion issues, three-dimensional spatial data provides more comprehensive structural information about crops (Herrero-Huerta et al., 2020; Sun et al., 2024). Methods for acquiring high-precision three-dimensional crop data primarily include Light Detection and Ranging (LiDAR) and photogrammetry. LiDAR has been successfully used in crop phenotyping, but its high cost and limited portability limit its application in agriculture (Paulus, 2019; Xiao et al., 2023). Photogrammetry, equipped with a consumer-grade high-resolution RGB camera, effectively generates dense point clouds and color information from sparse point clouds by fusing structure-from-motion (SFM) and multi-view stereo (MVS) (Westoby et al., 2012; Wu, 2013). Due to terrain restrictions and crop occlusion (Xiao et al., 2024), traditional straight-line oblique photography cannot capture a sufficient number of viewpoints in complex field environments (Sun et al., 2024). UAV-based CCO imaging, by adjusting flight altitude and shooting angle according to specific crop conditions, can cost-effectively obtain more comprehensive surface information and higher-quality point cloud data (Xiao et al., 2023; Zhu et al., 2023; Sun et al., 2024).

[0004] 3D data can effectively extract spatial structural attributes of the canopy (Jimenez-Berni et al., 2018; Deery et al., 2020; Herrero-Huerta et al., 2020). Currently, there is an urgent need for rapid and automated 3D analysis methods to efficiently extract crop phenotypes, link them to genotypes, and improve breeding efficiency (Duan et al., 2016). Canopy structural features can be obtained from 3D point clouds, such as the 3D voxel index (3DVI) and 3D silhouette coefficient (3DPI) (Jimenez-Berni et al., 2018; Deery et al., 2020), and the percentage of plant height at different heights (Luo et al., 2021). Although their effectiveness in estimating parameters such as crop biomass has been demonstrated, 3D canopy features remain insufficient for discerning trait differences across multiple genotypes when faced with hundreds or even thousands of crop breeding materials (Li et al., 2021; Ren et al., 2023). Point cloud feature description methods are categorized by feature extraction strategy into statistical feature-based and deep learning-based methods (Han et al., 2018). Deep learning-based methods preprocess raw 3D point cloud data into structured, fixed-size representations that can be efficiently processed by different network architectures (Qi et al., 2017; Zhou and Tuzel, 2017). However, the features derived from these methods struggle to interpret trait differences between crop varieties and require data annotation (Xiao et al., 2024). Statistical feature-based methods can intuitively represent the characteristics of crop canopy point cloud data, making the feature extraction process easier to understand and interpret.

[0005] Methods based on statistical features are mainly divided into label-based methods and histogram-based methods (Han et al., 2018). Histogram-based methods simply sum the geometric statistics of each region to form a histogram, effectively describing local features. However, the bin width of the histogram requires careful adjustment (Rusu et al., 2009). When the data distribution is uneven, histogram-based methods can lose some global information (Flint et al., 2007). Label-based methods, on the other hand, encode geometric information within the local area of ​​feature points to reflect the characteristics of the local area (Pauly et al., 2003; West et al., 2004). They have been widely used in fields such as urban scene classification (Weinmann et al., 2015; Hackel et al., 2016) and cotton boll identification and yield estimation (Sun et al., 2020; Xiao et al., 2024). This method has high classification accuracy in pattern recognition, but in the field of crop phenotyping, the extraction of local information ignores the overall three-dimensional structural information of the crop canopy, and local calculations require a large amount of memory space for storage.

[0006] In summary, there is an urgent need for rapid, automated 3D analysis methods to efficiently extract crop phenotypic traits, correlate them with genotypes, and optimize breeding strategies. However, for hundreds of breeding materials, characterizing 3D canopy traits in field plots is lacking, making it difficult to effectively identify trait differences between multiple genotypes. Summary of the Invention

[0007] The purpose of the present invention is to solve the problem of the lack of point cloud feature indicators in the field of agricultural point cloud applications, and to provide a canopy morphological 3D feature analysis method based on lightweight drones and point cloud singular value decomposition. The method combines cross-circle (CCO) flight line oblique photography with the SFM-MVS algorithm to obtain high-quality, low-cost canopy 3D point clouds of different varieties of crops. The label-based method of the present invention introduces the singular value decomposition (SVD) method for the point cloud matrix, and establishes 9 canopy morphological 3D features at global and local scales. The feature extraction process is automated, and the overall features of the canopy are efficiently extracted at the global scale, and information on the surface details of the canopy is obtained at the local scale. In addition, the significance test of the features was carried out under different planting densities, lodging grades and yield grades. At the same time, the 3D features of the canopy are also used in phenotypic estimation, such as yield estimation and lodging discrimination, to test its effectiveness.

[0008] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0009] A canopy morphology 3D feature analysis method based on a lightweight drone and point cloud singular value decomposition comprises the following steps:

[0010] Step 1: Field data collection:

[0011] During field data collection, canopy height, bottom pod height, number of main stem nodes, number of branches, number of pods per plant, number of grains per plant, 100-grain weight, and grain yield were measured for crops of different plot varieties. Agronomic knowledge was used to accurately determine the degree of crop lodging. Crops in each plot were harvested, threshed, and air-dried individually, and the crop yield was weighed and measured when its moisture content fell below 13%.

[0012] Step 2: UAV multi-sensor image acquisition:

[0013] UAVs can be equipped with high-resolution RGB to conduct aerial surveys of crops at different stages;

[0014] Step 3: Raw data preprocessing:

[0015] After the drone aerial operation is completed, the captured data is geometrically corrected and stitched to generate a digital elevation model and a digital orthophoto; the geometric correction and stitching processing specifically includes image alignment, GCP import and dense point cloud generation;

[0016] Step 4: Extracting canopy spectrum, point cloud structure, and texture information

[0017] Step 41: Extracting canopy spectral and texture information: Using high-resolution orthophotos, the CIVE index is established to segment plants and soil, thereby extracting soybean canopy coverage. A gray-level co-occurrence matrix is ​​used to extract RGB texture feature information.

[0018] Step 42: Extract global and local scale features based on soybean 3D canopy point cloud;

[0019] Step 5: Based on the global and local scale characteristics of the canopy point cloud calculated in Step 4, construct a time series analysis box plot based on yield grade and lodging grade, and perform a significance test. In addition, calculate the correlation between characteristics of different growth stages to characterize the differences within the characteristics.

[0020] Step 6: In phenotypic estimation, a machine learning model was constructed. A regression model was used for yield estimation, and a classification model was used for lodging discrimination. The maximum number of models was set to 10, the maximum run time was set to 60 seconds, and 5-fold cross-validation was used to improve the stability and reliability of the evaluation.

[0021] Furthermore, in step 2, the drone adopts two flight modes:

[0022] (1) When the flight altitude is 15 m, the shooting mode is oblique photography, the conventional route is adopted, and the image overlap rate is 85%;

[0023] (2) The flight altitude is 5 m, a cross-circling flight route is adopted, and the overlap rate between images is 50%.

[0024] Furthermore, in step three, for the images obtained by CCO route photography, the multi-view images captured along the CCO path are stitched, and the multi-view stereo and motion structure algorithms are used to reconstruct the point cloud of the soybean canopy; the image data collected along the CCO route are aligned to generate a dense point cloud for subsequent modeling; a vector layer is established, and a single plot is segmented to obtain the RGB orthophoto map, plant height and point cloud data of the plot; during the point cloud data processing, outliers and noise points are removed, and then the RANSAC method is used for ground area fitting, retaining high-precision point cloud data related to the soybean canopy.

[0025] Furthermore, in step 42, in the preprocessing part, the soybean three-dimensional canopy point cloud is screened based on color information and ground height information, and the ground area point cloud is eliminated; then, the noise points in the point cloud are eliminated using the remove radius outlier method of Open3D, where the neighborhood radius is set to 0.1m and the minimum number of neighborhood points is set to 20; finally, the point cloud is normalized by the center value.

[0026] Furthermore, in step 42, in the global scale canopy 3D feature extraction, the three-dimensional coordinates of the global point cloud are constructed into a matrix NP, and the eigenvectors and eigenvalues ​​of the matrix NP are obtained through singular value decomposition.

[0027] Furthermore, in step 42, in the local scale canopy 3D feature extraction, by traversing each point cloud, using the K nearest neighbor search area of ​​20, 50, 100, 200, and 500 point clouds of ROI, the small-scale point clouds of different levels are subjected to SVD to obtain eigenvalues ​​and eigenvectors, thereby obtaining , traverse each point cloud in turn and calculate the local canopy 3D features, and finally obtain the average value .

[0028] Compared with the prior art, the present invention has the following advantages:

[0029] This study uses a lightweight drone combined with a CCO flight path to reconstruct crop canopy point clouds with high precision. Furthermore, the SVD method for point cloud coordinate matrices was introduced to establish nine 3D canopy morphological features at different scales and evaluate their potential for phenotypic estimation.

[0030] 1. The global canopy 3D features are significantly correlated with crop yield, lodging, canopy height and vegetation index (p <= 0.05). When the planting density, lodging grade and yield grade are used as the classification basis, the characteristic entropy ( ), features and ( ) and the total variance ( ) index showed significant differences between the initial pod stage and the full pod stage (p<=0.05), and the sphericity ( ), curvature( ) and anisotropy ( ) index showed significant differences in the first grain to full grain early maturity stages (p<=0.01).

[0031] 2. Global canopy 3D features primarily characterize yield differences across crop yield grades from the initial grain stage to the full grain stage. Localized 3D feature extraction can extract more detailed structural information, helping to characterize crop variability across multiple growth stages and improving phenotypic estimation capabilities.

[0032] 3. Compared with the fusion of only one type of features, multimodal feature fusion can effectively improve the accuracy of phenotypic estimation. When all features are fused, the yield estimation accuracy is the best (RMSE = 362.10 kg ha -1 , RMAE=278.31 kg ha -1 , rRMSE%=14.31, r=0.63). The accuracy of distinguishing four types of lodging reached the best (Accuracy=0.575, F1-score=0.569), and the accuracy of distinguishing two types of lodging reached the best (Accuracy=0.901, F1-score=0.892).

[0033] In summary, 3D canopy characterization at different scales has great potential and could help accelerate crop breeding. Furthermore, future work will require the collection of data from multiple years and time points, and the use of data representative of varietal or genotype differences, to improve the applicability and stability of 3D canopy characterization. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 Schematic diagram of the test site's geographical location and site design; (a) Test area in Heihe City, Heilongjiang Province, China; (b) Orthophoto generated by stitching RGB images; (c) 3D crop canopy RGB images acquired using the cross-surround CCO flight method using a DJI Phantom 4 RTK; (d) 3D crop canopy reconstruction using the SFM-MVS algorithm, followed by removal of the ground point cloud and cropping to obtain a single plot point cloud.

[0035] Figure 2This is a flow chart for the calculation of 3D features of canopy morphology. The content is divided into two sections: global scale and local scale feature extraction. 1. In the global scale calculation, directly input all point clouds, and then use SVD to obtain the eigenvalues ​​( ) is used to calculate global 3D features. 2. In local scale calculations, by searching for nearby point clouds, K-nearest neighbors are used to search for different numbers of point clouds around each point cloud for feature calculation. By traversing all point clouds, the mean value represents the local 3D feature.

[0036] Figure 3 Figure 2 is the correlation matrix of point cloud structural features at different periods. The characteristic factors are yield, lodging grade, global canopy morphology 3D features, voxel ratio, plant height, and vegetation indices (EXR, EXG, EXGR). Red ellipses indicate positive correlations, and blue ellipses indicate negative correlations. Indicates a significant correlation at P <= 0.05. (a) indicates the initial flowering stage of the crop; (b) indicates the full flowering stage of the crop; (c) indicates the initial pod formation stage of the crop; (d) indicates the full pod formation stage of the crop; (e) indicates the initial grain formation stage of the crop; and (f) indicates the full grain and initial maturity stage of the crop.

[0037] Figure 4 The box plots of the global canopy 3D characteristics at different growth stages under different planting densities. (a) Linear ( ); (b) planarity ( ); (c) sphericity ( ); (d) characteristic entropy ( ); (e) characteristics and ( ); (f) curvature ( ); (g) total variance ( ); (h) anisotropy ( ); (i) trace ( D1 represents low-density planting (green); D2 represents medium-density planting (orange); and D3 represents high-density planting (red). Indicates p<=0.05, there is a significant difference. Indicates p<=0.01, indicating a very significant difference. Indicates p <= 0.001, indicating a highly significant difference. The purple line between the box plots represents the connecting line of the mean.

[0038] Figure 5 The box plots of the global canopy 3D characteristics at different growth stages under different lodging levels. (a) Linear ( ); (b) planarity ( ); (c) sphericity ( ); (d) characteristic entropy ( ); (e) characteristics and ( ); (f) curvature ( ); (g) total variance ( ); (h) anisotropy ( ); (i) trace ( ). Among them, L1 indicates no lodging (green), L3 indicates mild lodging (yellow), L5 indicates moderate lodging (orange), and L7 indicates severe lodging (red). Indicates p<=0.05, there is a significant difference. Indicates p<=0.01, indicating a very significant difference. Indicates p <= 0.001, indicating a highly significant difference. The purple line between the box plots represents the connecting line of the mean.

[0039] Figure 6 The box plots of the global canopy 3D characteristics at different growth stages under different yield levels are shown. The yield level is divided into four intervals according to 25%, 50%, 75%, and 100%. (a) Linear ( ); (b) planarity ( ); (c) sphericity ( ); (d) characteristic entropy ( ); (e) characteristics and ( ); (f) curvature ( ); (g) total variance ( ); (h) anisotropy ( ); (i) trace ( ). Among them, Y1 represents lower yield (red), Y2 represents low-medium yield (orange), Y3 represents medium-high yield (yellow), and Y4 represents high yield (green). Indicates p<=0.05, there is a significant difference. Indicates p<=0.01, indicating a very significant difference. Indicates p <= 0.001, indicating a highly significant difference. The purple line between the box plots represents the connecting line of the mean.

[0040] Figure 7 Boxplots of linear metrics at different scales, divided by production level and period. (a) Local scale with 20 points; (b) Local scale with 50 points; (c) Local scale with 100 points; (d) Local scale with 200 points; (e) Local scale with 500 points; (f) Global scale with all points included. Indicates p<=0.05, there is a significant difference. Indicates p<=0.01, indicating a very significant difference. Indicates p <= 0.001, indicating a highly significant difference. The purple line between the box plots represents the connecting line of the mean.

[0041] Figure 8 Boxplots of planarity metrics at different scales, divided by production level and period. (a) Local scale with 20 points; (b) Local scale with 50 points; (c) Local scale with 100 points; (d) Local scale with 200 points; (e) Local scale with 500 points; (f) Global scale with all points included. Indicates p<=0.05, there is a significant difference. Indicates p<=0.01, indicating a very significant difference. Indicates p <= 0.001, indicating a highly significant difference. The purple line between the box plots represents the connecting line of the mean.

[0042] Figure 9 Boxplots of sphericity at different scales, divided by production level and period. (a) Local scale with 20 points; (b) Local scale with 50 points; (c) Local scale with 100 points; (d) Local scale with 200 points; (e) Local scale with 500 points; (f) Global scale with all points included. Indicates p<=0.05, there is a significant difference. Indicates p<=0.01, indicating a very significant difference. Indicates p <= 0.001, indicating a highly significant difference. The purple line between the box plots represents the connecting line of the mean.

[0043] Figure 10 Boxplots of the characteristic entropy index at different scales, divided by production level and period. (a) The number of local point clouds is 20; (b) The number of local point clouds is 50; (c) The number of local point clouds is 100; (d) The number of local point clouds is 200; (e) The number of local point clouds is 500; (f) All point clouds at the global scale are included in the calculation.

[0044] Figure 11 Boxplots are calculated for features and metrics at different scales, divided by production level and period. (a) Local scale point cloud with 20 points; (b) Local scale point cloud with 50 points; (c) Local scale point cloud with 100 points; (d) Local scale point cloud with 200 points; (e) Local scale point cloud with 500 points; (f) Global scale point cloud with all points included.

[0045] Figure 12Boxplots of curvature indices calculated for different scales, divided by production level and period. (a) Local scale point cloud number is 20; (b) Local scale point cloud number is 50; (c) Local scale point cloud number is 100; (d) Local scale point cloud number is 200; (e) Local scale point cloud number is 500; (f) All point clouds in the global scale are included in the calculation.

[0046] Figure 13 Boxplots of the total variance at different scales, divided by production level and period. (a) 20 local point clouds; (b) 50 local point clouds; (c) 100 local point clouds; (d) 200 local point clouds; (e) 500 local point clouds; (f) all global point clouds included.

[0047] Figure 14 Boxplots of anisotropy indices calculated for different scales, divided by production level and period. (a) Local scale with 20 point clouds; (b) Local scale with 50 point clouds; (c) Local scale with 100 point clouds; (d) Local scale with 200 point clouds; (e) Local scale with 500 point clouds; (f) Global scale with all point clouds participating in the calculation.

[0048] Figure 15 Boxplots of trace metrics calculated for different scales, divided by production level and period. (a) Local scale with 20 point clouds; (b) Local scale with 50 point clouds; (c) Local scale with 100 point clouds; (d) Local scale with 200 point clouds; (e) Local scale with 500 point clouds; (f) Global scale with all point clouds participating in the calculation.

[0049] Figure 16 Yield estimation results for different feature combinations at different crop growth stages. Results obtained using H2O-AutoML with 5-fold validation. ALL indicates that all features were used for modeling. (a) RMSE; (b) MAE; (c) rRMSE%; (d) Pearson's.

[0050] Figure 17 Results of different feature combinations for distinguishing four types of lodging at different crop growth stages. Results obtained using H2O-AutoML with 5-fold validation. ALL indicates that all features were used for modeling. (a) Accuracy; (b) Precision; (c) Recall; (d) F1-score.

[0051] Figure 18Results of different feature combinations for different lodging classifications at different crop growth stages. Results obtained using H2O-AutoML with 5-fold validation. ALL indicates that all features were used for modeling. (a) Accuracy; (b) Precision; (c) Recall; (d) F1-score. DETAILED DESCRIPTION

[0052] The technical solution of the present invention is further described below with reference to the accompanying drawings, but is not limited thereto. Any modification or equivalent replacement of the technical solution of the present invention that does not depart from the spirit and scope of the technical solution of the present invention should be included in the scope of protection of the present invention.

[0053] The present invention uses an RGB lens mounted on a drone combined with the SfM-MVS algorithm to perform high-precision three-dimensional reconstruction of the soybean canopy structure. At the same time, a label-based method introduces the singular value decomposition (SVD) method for the point cloud coordinate matrix, and establishes 9 canopy morphological 3D features at multiple scales as a new interpretable feature at the crop plot level. The results show that the global canopy 3D features are significantly correlated with soybean yield, lodging, canopy height and vegetation index (p<=0.05). 3D feature extraction based on the local range can extract more detailed structural information of the canopy surface, which helps to characterize the differences in more growth periods of soybeans. Finally, when the vegetation index, texture and canopy 3D features at different scales are fused, the yield estimation accuracy using H2O-AutoML is the best (RMSE=362.10 kg ha -1 , rRMSE%=14.31), and both the four-category and two-category lodging classification accuracy reached the highest level (4-category: Accuracy=0.575, F1-score=0.569, 2-category: Accuracy=0.901, F1-score=0.892). In summary, the high-precision 3D reconstruction and canopy feature extraction technology developed in this study based on a UAV platform demonstrates great potential and provides rapid and reliable data support for breeding research.

[0054] The present invention provides a canopy morphology 3D feature analysis method based on a lightweight unmanned aerial vehicle and point cloud singular value decomposition. Taking soybean as an example, the method includes the following steps:

[0055] 1. Test location and experimental setup

[0056] The test site is located in Heihe City, Heilongjiang Province, China (47°42′-51°03′ N, 124°45′-129°18′ E). Heihe City is located at the eastern end of the Greater Khingan Range and the northern end of the Lesser Khingan Range ( Figure 1a). It has a cold temperate continental monsoon climate, with warm springs and strong winds, and abundant rainfall in summer. Figure 1 b) The crops are all soybeans, and the ridge is 3.1 m long, 0.6 m wide, and covers an area of ​​approximately 1.9 m 2 The combination of three ridges is the area of ​​a plot for one variety, the spacing between ridges is about 0.3 m, and the area of ​​the plot is about 5.9 m 2 .

[0057] A total of 186 soybean varieties were planted in the study area, with a total of 558 planting plots, of which the low-density (D1) planting density was 3.0×10 5 seed hm -1 The planting density of medium density (D2) is 3.5×10 5 seed hm -1 The high density (D3) is 4.2×10 5 seed ha -1 All soybean varieties were planted on May 14, 2023, sown on May 30, 2023, and harvested on October 19, 2023. Field management was carried out according to conventional practices.

[0058] 2. Data Collection

[0059] 2.1 Field Data Collection

[0060] A quantitative evaluation of soybean canopy height, bottom pod height, number of main stem nodes, number of branches, number of pods per plant, number of grains per plant, 100-grain weight, grain yield, and lodging grade was conducted for 558 plots. Lodging grade was determined by visual inspection and categorized into four levels: no lodging (grade 1), mild lodging (grade 3), moderate lodging (grade 5), and severe lodging (grade 7). Each plot was harvested, threshed, and dried individually. Soybeans were weighed when their moisture content was below 13%. Yield was expressed in kg ha. -1 Table 1 shows the field statistics.

[0061] Table 1 Descriptive statistics of field measurement parameters

[0062]

[0063] 2.2. UAV Data Collection Methods

[0064] The drone survey was conducted in multiple stages on July 1, 2023 (first flowering), July 10 (full flowering), July 20 (first pod), July 28 (full pod), August 14 (first grain), and September 2 (full grain and first ripening). The drone model was DJI Phantom 4RTK, and the RGB image resolution obtained was 4608×3456 pixels. The drone used two flight modes: (1) When operating at an altitude of 15 meters, oblique photography was used with 85% forward and lateral overlap for image capture. The spatial resolution of the RGB image was 0.38 cm / pixel. (2) In order to obtain high-precision canopy point cloud data ( Figure 1 c) The drone was flown at an altitude of 5 m, using a cross-circle (CCO) pattern (Xiao et al., 2023). The overlap ratio was 50%, and the camera tilt angle was 45°. All drone imagery was collected during noon (11:00–14:00) in clear, windless weather.

[0065] 2.2 Raw Data Preprocessing

[0066] After completing the UAV orthophoto shooting mission, the captured data was geometrically corrected and stitched using Agisoft PhotoScan Professional software (Agisoft LLC, St. Petersburg, Russia). This includes image alignment, GCP import, and dense point cloud generation. Finally, a digital elevation model and a digital orthophoto were generated. For the images acquired by CCO route photography, Agisoft PhotoScan Professional was used in the study to stitch the multi-view images captured along the CCO path, and the multi-view stereo and motion structure (SFM-MVS) algorithm was used to reconstruct the point cloud of the soybean canopy. Agisoft PhotoScan Professional software performed image alignment on the image data collected along the CCO route to generate a dense point cloud ( Figure 1 d) for subsequent modeling. QGIS 3.24 software was used to create a vector layer and segment the individual plots ( Figure 1 d), obtain the RGB orthophoto map, plant height and point cloud data of the plot. In the process of point cloud data processing ( Figure 1 d) First, outliers and noise points are removed, and then the ground area is fitted using the RANSAC method (Yang et al., 2022). This allows the ground and soybean plants in the data to be separated. Ultimately, high-precision point cloud data related to the soybean canopy is retained.

[0067] 3. Modeling methods and phenotypic information extraction

[0068] 3.1. Extraction of canopy spectrum and texture information

[0069] The orthophoto image constructed by the drone consists of three bands: red, green and blue. Based on these three bands, the following Figure 2 The vegetation index shown is used for yield estimation and lodging identification. Furthermore, the CIVE index was developed using high-resolution orthophotos to segment plants and soil, thereby extracting soybean canopy cover (CC). Furthermore, texture information has been shown to be effective for yield and lodging monitoring in high-throughput field plant phenotyping using drones (Nichol and Sarker, 2011). Therefore, the present invention uses the gray-level co-occurrence matrix (GLCM) (Haralick and Shanmugam, 1973) to extract RGB texture feature information. After converting the RGB image to grayscale, homogeneity (HO), contrast (CO), dissimilarity (DI), characteristic entropy (EN), second-order moment (SE), and correlation (CO) are extracted, as shown in Table 2.

[0070] Table 2 Extraction of vegetation index and texture features

[0071]

[0072] 3.2 Extraction of canopy point cloud structure information

[0073] The canopy point cloud structure information established by the present invention mainly includes canopy morphological 3D features, canopy height and voxel ratio features. The point cloud surface normal vector is an important attribute of the geometric surface. In the present invention, the problem of estimating the surface normal is converted into calculating the eigenvectors and eigenvalues ​​of the covariance matrix. In the present invention, global scale and local scale feature extraction is constructed based on the soybean three-dimensional canopy point cloud, such as Figure 2 As shown in Algorithm 1. In the preprocessing phase, the soybean canopy 3D point cloud is filtered based on color and ground height information, removing ground-level points. Next, noise points are removed using the Open3D noise removal method, with a neighborhood radius set to 0.1m and a minimum number of neighborhood points set to 20. Finally, the point cloud is normalized using its center value.

[0074] In global scale canopy 3D feature extraction ( Figure 2 Global scale), the three-dimensional coordinates of the global point cloud are constructed into a matrix NP, and the eigenvectors and eigenvalues ​​of the matrix NP are obtained through singular value decomposition (SVD). First, the covariance matrix of the global point cloud is calculated , and then find The eigenvalue of and eigenvectors , as shown in Equations 1 and 2. Calculate 3D feature linearity ( ), planarity ( ), sphericity ( ), characteristic entropy ( ), features and ( ), curvature( ), total variance ( ), anisotropy ( ),trace( )(Weinmann et al., 2015), the calculation method is shown in Table 3.

[0075]

[0076]

[0077] In the local scale canopy 3D feature extraction ( Figure 2 Local scale) is different from global scale in that it takes into account the local characteristic information of soybean canopy. By traversing each point cloud, 20, 50, 100, 200, and 500 point clouds of ROI (full name first) within the K nearest neighbor search area are used to perform SVD on small-scale point clouds of different levels to obtain eigenvalues. and eigenvectors , thus obtaining Traverse each point cloud in turn and calculate the local canopy 3D features, and finally obtain the average value .

[0078] The plant height and point cloud voxel ratio of soybeans at different growth stages are also important information. The present invention obtains high-precision canopy point clouds of soybeans through the CCO route, and subtracts the ground point cloud elevation from the canopy point cloud elevation to obtain canopy height information. The canopy height information extraction is divided into the point cloud height average, maximum, and cumulative height 95% quantile. The point cloud voxel ratio is an indicator used to evaluate the sparsity of point cloud data. The voxel ratio is usually defined as the ratio of the number of actual points in the point cloud to the number of non-empty (containing points) voxels in the voxels after point cloud voxelization (Liu et al., 2021). In the study, the sizes of the voxelized grids were set to (0.05, 0.10, 0.20, 0.30, 0.40, 0.50) for calculation.

[0079] Table 3 3D feature extraction methods

[0080]

[0081]

[0082] 3.3 Model building method and environment

[0083] In the present invention, the development language is Python 3.7.5, the GPU model is Tesla V100, the CPU is 2 Cores, and the memory is 32GB. The 3D features of the canopy point cloud are calculated using Open3d 3.7 and Numpy 1.25.0 libraries. In phenotypic estimation, a machine learning model was constructed using H2O-AutoML 3.46 (LeDell and Poirier, 2020). A regression model was used for yield estimation, and a classification model was used for lodging discrimination. The maximum number of models was set to 10, and the maximum running time was set to 60s. In order to reduce the randomness of phenotypic estimation and improve the generalization ability of the model, 5-fold cross-validation was used in the study to improve the stability and reliability of the evaluation.

[0084] 3.4 Statistical methods

[0085] 3.4.1 Regression Indicators

[0086] The root mean square error (RMSE), root mean square error (RMAE), relative RMSE (rRMSE%), and correlation coefficient r are used to evaluate the estimation accuracy of soybean yield, etc., as shown in formula 3-6. and These are the measured and estimated values ​​of soybean yield, respectively. is the average value of the measured soybean yield, and n is the total number of samples in the test set.

[0087]

[0088]

[0089]

[0090]

[0091] 3.4.2 Lodging classification indicators

[0092] Accuracy, precision, recall, and F1-score are used to evaluate the accuracy of crop lodging and measure the training quality of the model, as shown in Equations 7-10. TP is the number of true positive samples, FP is the number of false positive samples, and FN is the number of false negative samples.

[0093]

[0094] Precision=

[0095] Recall=

[0096]

[0097] 3.4.3 ANOVA test

[0098] To explore the performance and applicability of canopy 3D features across different planting densities, lodging classes, and yield classes, a one-way analysis of variance (ANOVA) was conducted using the HSD Tukey test (α = 0.001) to examine significant differences in canopy 3D features across planting densities, lodging classes, and yield classes. Yield classes were determined by categorizing measured data into four classes, each consisting of the 25th, 50th, 75th, and 100th percentiles of yield. 0-25% was defined as low yield (Y1), 25-50% as low-medium yield (Y2), 50%-75% as medium-high yield (Y3), and 75-100% as high yield (Y4).

[0099] 4. Results

[0100] 4.1 Multi-period correlation of different characteristic factors

[0101] At different growth stages, yield, lodging grade, global canopy morphology 3D features, voxel ratio, plant height and vegetation index (EXR, EXG, EXGR) were used to construct the correlation matrix, e.g. Figure 3 As shown, The representative calculation of Pearson's calculation showed a significant correlation with p <= 0.05. In the analysis of different growth stages, the point cloud structure characteristics showed a certain correlation with yield, lodging grade and vegetation index. ) as an example, in the initial flowering to full grain and initial maturity stage, it showed a significant negative correlation with yield (p <= 0.05). The higher the value, the lower the yield. In the initial flowering to full flowering stage, since lodging has not yet occurred, there is no significant correlation. In the initial grain to full grain early maturity period, soybean lodging has occurred, so it shows a significant positive correlation with the lodging grade (Pearson's initial grain: 0.40; full grain early maturity: 0.44; p <= 0.05). This shows that the higher the lodging grade, the more serious the lodging phenomenon. In addition, during the initial flowering to full grain and initial ripening stage, There is a significant positive correlation with canopy height (H95%, Hmean, Hmax) (Pearson's: 0.08-0.49, p<=0.05). There was also a significant correlation with VIs (Pearson's: 0.08-0.38, p<=0.05), and in the pod-full grain stage, The correlation between plant height and yield was weak (Pearson's: -0.02-0.15, p<=0.05). The point cloud voxel ratio feature showed a positive correlation during the initial flowering to full flowering period, but the correlation was weaker during the initial pod to full grain and early maturity period. The point cloud-extracted canopy plant height feature (H95%, Hmean, Hmax) showed a significant positive correlation with yield during the initial flowering to full pod period (Pearson's: 0.12-0.36, p<=0.05). During the initial grain to full grain and early maturity period, plant height showed a weak negative correlation with yield (Pearson's: -0.28 to -0.09, p<=0.05), as early-maturing soybean varieties experienced leaf shedding and lodging. However, plant height showed a significant positive correlation with lodging grade (Pearson's: 0.33-0.55, p<=0.05).

[0102] 4.2. Global canopy 3D feature time series results under different planting densities

[0103] The global canopy 3D feature time series results were evaluated at different planting densities ( Figure 4 For linear ( ) and planarity ( )index( Figure 4 ab), in the initial flowering to full flowering stage, D1, D2 and D3 showed significant differences (p <= 0.05). In the early growth stage of soybean, the higher the planting density, The values ​​are relatively high, while For spherical ( ), curvature( ) and anisotropy ( )index( Figure 4 cfh), significant differences were found in the initial pod, initial grain, and full grain early maturity stages. Especially in the full grain early maturity stage, there were extremely significant differences among D1, D2, and D3 (p<=0.001). As the planting density increased, the point cloud features of soybeans under D3 showed higher divergence and curvature. 、 The values ​​increased relatively The value is relatively reduced. For the characteristic entropy ( ), features and ( ) and the total variance ( )index( Figure 4 deg), mainly in the initial pod-initial grain stage, there was a significant difference (p<=0.001). The higher the planting density, There is an increase, and It shows a decrease. )index( Figure 4i) There were significant differences in the pod-early grain stage (p <= 0.01). The higher the planting density, The lower.

[0104] 4.3. Global canopy 3D feature time series results at different lodging levels

[0105] The global canopy 3D feature time series results were evaluated at different lodging levels ( Figure 5 For linear ( ) and planarity ( )index( Figure 5 ab), the differences in different lodging grades were small during the initial flowering to full grain and initial maturity period. ), curvature( ) and anisotropy ( )index( Figure 5 cfh), in the initial grain-full grain early maturity stage, there are extremely significant differences among L1, L3, L5 and L7 (p<=0.001). As the soybean lodging becomes more serious, the point cloud features show higher divergence and curvature. 、 The values ​​increased significantly, and The value is relatively reduced. ), features and ( ) and the total variance ( )index( Figure 5 deg), mainly showing significant differences in the initial pod-full pod stage (p<=0.05). The more severe the soybean lodging, There is an increase, and It shows a decrease. )index( Figure 5 i) There were significant differences in different lodging grades during the pod-filling period (p <= 0.01). The higher the planting density, The lower.

[0106] 4.4. Global canopy 3D feature time series results under different yield levels

[0107] The global canopy 3D feature time series results were evaluated at different yield levels ( Figure 6 For linear ( ) and planarity ( )index( Figure 6 ab), the differences in different lodging grades were small during the initial flowering to full grain and initial maturity period. ), curvature( ) and anisotropy ( )index( Figure 6cfh), in the initial grain-full grain stage, there were extremely significant differences between Y1 and Y2, Y3, and Y4 (p<=0.001). The lower the soybean plot yield, the higher the divergence, curvature, and anisotropy of the canopy 3D characteristics. 、 The values ​​increased significantly, and The value is relatively reduced. ), features and ( ) and the total variance ( )index( Figure 6 deg), mainly showing significant differences in the initial grain stage (p<=0.01). The higher the yield of the soybean plot, It shows a gradual increase, and and It shows a decrease. )index( Figure 6 i) There was a significant difference between Y1 and Y4 at the initial grain stage (p <= 0.001). The higher the plot yield, the The lower the value.

[0108] 4.5. Time series results of local extraction of 3D features of different canopies

[0109] The 3D characteristics of the canopy at different scales were evaluated. The yield grade was used as the classification, and the significant differences in the 3D characteristics of the canopy at different stages were statistically analyzed. Figure 7-8 ), although the feature differences at the global scale are not significant ( Figure 7 f), but in the primary particle stage, when the local range ROI=(50, 100, 200) ( Figure 7 ad, Figure 8 be), can obtain more detailed characteristics of soybeans, so there are significant differences (p <= 0.01). Yield grade Y1 and (Y2, Y3, Y4) show significant differences. As the yield grade increases, and Shows an upward trend. )index( Figure 15 ), the difference between local features and global features is small.

[0110] For spherical ( ), curvature( ) and anisotropy ( )index( Figure 9 , Figure 12 , Figure 14), with significant differences appearing at different scales. At the global scale, significant differences were observed between Y1 and other yield grades during the first kernel to full kernel and early maturity stages (p <= 0.001). At ROIs (20, 50, 100), significant differences were observed between the low-yield grade and Y4 during the first kernel stage (p <= 0.01). At ROIs (500), more detailed information was captured within the local area, resulting in significant differences across the first pod to full kernel and early maturity stages (p <= 0.05). Index in global scope and local scope ( Figure 9 af), as the production level continues to increase, All showed a downward trend.

[0111] For the characteristic entropy ( ), features and ( ) and the total variance ( )index( Figure 10 , Figure 11 , Figure 13 ), and The differences between different yield grades are small under local feature extraction. However, When ROI=(20,50), local feature extraction shows more significant differences, such as the difference between Y4 and other yield grades. At the same time, local feature extraction and global feature extraction show an opposite upward trend ( Figure 13 ), although at the global scale, the higher the yield level, The lower the value, the higher the yield level will be because of the more detailed information obtained from local extraction. The higher the state.

[0112] 4.6 Results of Yield Estimation and Lodging Identification Using Different Strategies

[0113] The accuracy of yield estimation of different soybean trait combination strategies was evaluated at different growth stages ( Figure 16 At the pod-filling stage, when only 3D (global) canopy features were used, the yield estimation RMSE was 428.05 kg ha. -1 When the 3D (global + local) features were combined with local information, the yield estimation RMSE was reduced by 5.16 kg ha -1 When VIs and 3D (global + local) combined features are used, the yield estimation RMSE is 376.10 kg ha -1 , RMAE=285.63 kg ha -1 When different features were fused, the yield estimation accuracy was the best (RMSE = 362.10 kg ha -1, RMAE = 278.31 kg ha -1 , rRMSE% =14.31, r = 0.63).

[0114] The results of the four categories of lodging discrimination with different feature combinations at different growth stages are as follows: Figure 17 As the soybean growth period progresses, lodging becomes increasingly pronounced, and the F1-score for distinguishing the four types of lodging continues to improve. At the full-grain, early-ripening stage, the F1-score for four-category discrimination using VIs was 0.438. This improvement was 0.024 when using 3D canopy (global) features, and 0.092 when using 3D (global + local) features. This indicates that combining 3D canopy features can effectively improve the accuracy of lodging discrimination. The combination of VIs and 3D (global + local) features achieved the best accuracy for distinguishing the four types of lodging (Accuracy = 0.575, Precision = 0.571, Recall = 0.575, F1-score = 0.569).

[0115] The results of different feature combinations for distinguishing the two types of lodging at different growth stages are as follows: Figure 18 As shown in Figure 2, the accuracy of two-category lodging discrimination gradually improves from the initial flowering stage to the full grain and initial maturity stage. In the full grain and initial maturity stage, when using only VIs for lodging discrimination, the F1-score is 0.855, which is 0.021 higher than when using 3D (global + local) features. When all features are combined, the accuracy reaches the highest level (Accuracy = 0.901, Precision = 0.903, Recall = 0.901, F1-score = 0.892).

[0116] 5.1 Global 3D Feature Interpretation of Canopy Morphology

[0117] Although 3D features of canopy morphology have been named (Pauly et al., 2003; West et al., 2004; Weinmann et al., 2015), their implications for crop phenotyping still require further evaluation. Figure 3 ), linear( ), planarity ( ) and traces ( ) index showed a significant correlation. The correlation test was conducted with the lodging grade in the initial flowering, initial pod and full pod stages (p<=0.05), but there was no significant correlation in the initial grain-full grain early maturity stage, indicating that , , It helps to characterize the variability of early crop lodging. In addition, (Dittrich et al., 2017) showed that linear ( ) and planarity ( ) is susceptible to noise, which makes the classification task susceptible.

[0118] For spherical ( ), curvature( ) and anisotropy ( ) index, in the correlation test, and showed a significant positive correlation (Pearson's = 0.98-1.00, P <= 0.05), and There was a significant negative correlation (Pearson's = -1.00 - -0.98, P <= 0.05). In the first grain to full grain early maturity period, the soybean lodging condition was more obvious, so it showed a more significant positive correlation with the lodging grade ( Figure 3 ef). Through the significance test of different lodging grades and yield grades ( Figure 5-6 ), indicating that the stronger the soybean lodging grade, the stronger the surface curvature and structural divergence of the point cloud. , The value increases, The value decreases. As the soybean yield level increases, the curvature and divergence of the point cloud surface decrease. Therefore, 、 as well as It has the potential to be used to describe soybean lodging traits and to estimate yield. Figure 4 efh), the increase in planting density will also lead to an increase in the curvature and divergence of the point cloud structure.

[0119] For the characteristic entropy ( ), features and ( ) and the total variance ( ) index, in the correlation test, and( , ) showed a significant negative correlation (Pearson's = -1.00 - -0.90, P <= 0.05). It is defined by Shannon entropy (Shannon, 1948), which measures the complexity of the point cloud information it contains (Weinmann et al., 2015; Dittrich et al., 2017). In the soybean pod-to-early stage, there is a significant correlation with yield and lodging ( Figure 3de). At the same time, the significance analysis of different planting densities was conducted ( Figure 4 ), D1 and D3 showed significant differences. The lower the value, the more significant the structure of the point cloud data is under low density, while under high density, the point cloud structure is more complex and tends to be chaotic. In the significance test of different lodging levels, the soybeans that did not lodge at the initial pod-full pod stage The value is lower, showing significant differences from L5 and L7 (p <= 0.001), indicating that for soybean plots without lodging, the point cloud structure is stronger, while for plots with moderate and severe lodging, the point cloud structure is more chaotic and disordered. Therefore, it is consistent with the study of (Dittrich et al., 2017) that the curvature ( ) and characteristic entropy ( ) features are more suitable for classification tasks.

[0120] 5.2 Effect of scale on canopy 3D characteristics

[0121] In technologies such as photogrammetry and computer vision, scale selection is crucial for automated analysis of canopy 3D point clouds (Dittrich et al., 2017). (Pauly et al., 2003) Neighborhood scale selection relies on covariance features of local surface variation. (Demantké et al., 2012) Covariance features based on linearity, planarity, and divergence are used for dimensional scale selection. (Weinmann et al., 2015) Covariance features based on feature entropy are used for eigen-feature entropy-based scale selection. The latter two methods consider the scale corresponding to the minimum Shannon entropy as the optimal scale. However, these scale selection methods are primarily used for semantic classification tasks of large-scale point clouds (Weinmann et al., 2015; Hackel et al., 2016). In multi-variety soybean phenotyping, especially at the plot scale, canopy 3D feature calculation is limited to the individual plots. Using a global point cloud to calculate canopy 3D features can effectively estimate yield and identify lodging ( Figure 16-18 ), and characterize the differences in different planting densities, lodging grades, and yield grades ( Figure 4-6 In addition, compared with the global scale, 3D feature extraction at the local scale can obtain more detailed structural information (Weinmann et al., 2015; Hackel et al., 2016), which is helpful for characterizing the differences in other reproductive stages ( Figure 7-15 For example, linear features are not very useful at the global scale ( Figure 9 f) showed no difference, while when ROI=(20, 50, 100, 200) ( Figure 9 ad), there were significant differences in yield levels at the primary grain stage (p <= 0.01). Figure 9 f), mainly manifested in significant differences in yield levels from the initial grain to the full grain initial maturity stage. However, at the local scale, when ROI=500 ( Figure 9 e) Significant differences in yield were found in different grades from the initial pod stage to the full grain and initial maturity stage.

Claims

1. A 3D feature analysis method of canopy morphology based on lightweight UAV and point cloud singular value decomposition, characterized by The method comprises the following steps: Step 1: Field data collection: During field data collection, canopy height, bottom pod height, number of main stem nodes, number of branches, number of pods per plant, number of grains per plant, 100-grain weight, and grain yield were measured for crops of different plot varieties. Agronomic knowledge was used to accurately determine the degree of crop lodging. Crops in each plot were harvested, threshed, and air-dried individually, and the crop yield was weighed and measured when its moisture content fell below 13%. Step 2: UAV multi-sensor image acquisition: UAVs can be equipped with high-resolution RGB to conduct aerial surveys of crops at different stages; Step 3: Raw data preprocessing: After the drone's aerial operation is completed, the captured data is geometrically corrected and stitched to generate a digital elevation model and a digital orthophoto. The geometric correction and stitching processing specifically includes image alignment, GCP import, and dense point cloud generation. For images acquired by CCO flight line photography, multi-view images captured along the CCO path are stitched, and multi-view stereo and motion structure algorithms are used to reconstruct the soybean canopy point cloud. Image alignment processing is performed on the image data collected along the CCO route to generate a dense point cloud for subsequent modeling. A vector layer is established to segment individual plots to obtain the RGB orthophoto map, plant height, and point cloud data of the plot. During the point cloud data processing, outliers and noise points are removed, and then the RANSAC method is used for ground area fitting, retaining high-precision point cloud data related to the soybean canopy. Step 4: Extracting canopy spectrum, point cloud structure, and texture information Step 41: Extracting canopy spectral and texture information: Using high-resolution orthophotos, the CIVE index is established to segment plants and soil, thereby extracting soybean canopy coverage. A gray-level co-occurrence matrix is ​​used to extract RGB texture feature information. Step 42: Extract global and local scale features based on soybean 3D canopy point cloud; Step 5: Based on the global and local scale characteristics of the canopy point cloud calculated in Step 4, construct a time series analysis box plot based on yield grade and lodging grade, and perform a significance test. In addition, calculate the correlation between characteristics of different growth stages to characterize the differences within the characteristics. Step 6: In phenotypic estimation, a machine learning model was constructed. A regression model was used for yield estimation, and a classification model was used for lodging discrimination. The maximum number of models was set to 10, the maximum run time was set to 60 seconds, and 5-fold cross-validation was used to improve the stability and reliability of the evaluation.

2. The canopy morphology 3D feature analysis method based on a lightweight UAV and point cloud singular value decomposition according to claim 1 is characterized by: In step 2, the drone uses two flight modes: (1) When the flight altitude is 15 m, the shooting mode is oblique photography, the conventional route is adopted, and the image overlap rate is 85%; (2) The flight altitude is 5 m, a cross-circling flight route is adopted, and the overlap rate between images is 50%.

3. The canopy morphology 3D feature analysis method based on a lightweight UAV and point cloud singular value decomposition according to claim 1 is characterized by: Step 42: In the preprocessing part, the soybean canopy 3D point cloud is screened based on color information and ground height information to remove the ground area point cloud; then, the noise points in the point cloud are removed using the Open3D remove radius outlier method, where the neighborhood radius is set to 0.1m and the minimum number of neighborhood points is set to 20; finally, the point cloud is normalized by the center value.

4. The canopy morphology 3D feature analysis method based on a lightweight UAV and point cloud singular value decomposition according to claim 1 is characterized by: Step 42: In the global scale canopy 3D feature extraction, the three-dimensional coordinates of the global point cloud are constructed into a matrix NP, and the eigenvectors and eigenvalues ​​of the matrix NP are obtained through singular value decomposition.

5. The canopy morphology 3D feature analysis method based on a lightweight UAV and point cloud singular value decomposition according to claim 1 is characterized by: Step 42: In the local scale canopy 3D feature extraction, by traversing each point cloud, using the K nearest neighbor search area of ​​20, 50, 100, 200, and 500 point clouds of ROI, the small-scale point clouds of different levels are subjected to SVD to obtain eigenvalues ​​and eigenvectors, thus obtaining , traverse each point cloud in turn and calculate the local canopy 3D features, and finally obtain the average value .

Images