Group canopy configuration index construction and estimation method

Abstract

The invention discloses a canopy leaf vertical distribution index construction method, and relates to the field of crop monitoring, and the method comprises the steps: S1, manually measuring the plant height (PH) of each cell of corn, dividing the plant height into five layers according to the PH, and obtaining LAI values of the whole canopy and leaves of each layer; s2, fitting cumulative LAI values of different leaf layers by using a Logistic equation; s3, calculating a first-order derivative of the fitted curve to obtain a canopy LAD distribution curve; and S4, extracting the kurtosis (Kurtosis) index and the skewness (Skewness) index of the LAD curve. Meanwhile, a rapid and accurate estimation method for canopy configuration indexes is provided. The method has the advantages that a novel canopy leaf vertical distribution index is constructed according to accumulated LAI values of different heights; and fusing unmanned aerial vehicle multispectral, RGB and LIDAR data, estimating canopy LAI, plant height and constructed indexes, and realizing rapid identification of the corn plant type.

Description

Technical Field

[0001] The present invention relates to the field of crop monitoring, and in particular to a method for constructing and estimating a group canopy configuration index. Background Art

[0002] In the summary of corn planting, how to effectively tap the photosynthetic potential of corn to increase yield is one of the current research focuses. Studies have shown that dense planting can optimize the canopy structure of the corn population and more effectively coordinate the relationship between individuals and groups, thereby significantly increasing the yield per unit area of ​​corn. Constructing an "ideal plant type" with a reasonable canopy leaf distribution can achieve the goal of increasing yield by raising the upper limit of planting density. The "ideal plant type" characteristics of corn are characterized by a canopy structure with compact upper leaves and relatively stretched middle and lower leaves. However, the lack of relevant plant type indicators is one of the main challenges hindering plant type research.

[0003] The vertical distribution of maize leaves is uneven, and the vertical distribution of leaves is an important determinant of light, water and nutrient distribution. Compared with the leaf area index (LAI), the leaf area density (LAD) quantifies the leaf area per unit volume and is an effective indicator for characterizing the vertical distribution of canopy leaves. However, the unified study of LAD at different heights is inefficient and cumbersome, and further quantification of LAD distribution is needed. In addition, compared with the determination of LAI, the direct determination of LAD is more cumbersome, requiring the leaf area at multiple heights of the canopy to be counted separately, and only limited sampling point data can be obtained, which is not only time-consuming and labor-intensive, but also destructive. Therefore, it is not suitable for long-term monitoring of the dynamic growth of leaf area in space and time.

[0004] Laser radar (LIDAR) is an active sensor that accurately obtains the height information of the surface and objects by emitting laser pulses and measuring their reflection time. Due to its unique three-dimensional observation capability, LIDAR is widely used to accurately quantify the horizontal and vertical structure of corn canopies. Oblique photography technology uses the camera to tilt and shoot at different angles, combines the principles of photogrammetry and computer vision, and uses technologies such as feature point matching and bundle adjustment to generate a visualized three-dimensional scene. The point cloud generated by oblique photography has low accuracy and serious distortion, occlusion and noise problems, but the point cloud density and refinement are high; the accuracy of LIDAR point cloud is high, and the point cloud generated by LIDAR is higher in accuracy and more complete, but the density is low and lacks texture information. The fusion of oblique photography technology and LIDAR technology, through the complementary advantages of different sensors, organically combines visual information with precise geometric information, and can capture the structural information of vegetation more accurately and comprehensively. In view of this, the present invention proposes a method for constructing and estimating a group canopy configuration index. Summary of the invention

[0005] The purpose of the present invention is to solve the problems in the prior art and to propose a method for constructing and estimating a group canopy configuration index.

[0006] A group canopy configuration index, the construction method of which comprises the following steps:

[0007] S1. Manually measure the plant height pH of each corn plot, divide it into 5 layers according to plant height, and obtain the LAI value of the entire canopy and leaves in each layer;

[0008] S2. Use the Logistic equation to fit the cumulative LAI values ​​of different leaf layers. The Logistic equation formula is:

[0009]

[0010] Among them, H is the relative height, M is the theoretical maximum LAI value, k is the initial value parameter, and b is the growth rate parameter;

[0011] S3. Calculate the first-order derivative of the fitting curve to obtain the canopy LAD distribution curve. The first-order derivative calculation formula is:

[0012]

[0013] Where, d is the difference;

[0014] S4. Extract the kurtosis and skewness index of the LAD curve. The calculation formula of the kurtosis and skewness is:

[0015]

[0016] Where μ represents the mean value of LAD and σ represents the standard deviation of LAD.

[0017] In the construction of the above-mentioned group canopy configuration index, in step S1, multiple corn plants with uniform growth are selected for height measurement, and the leaf area of ​​each layer is measured using a photographic leaf area meter in the laboratory, and the LAI is calculated based on the planting density and the plot area.

[0018] A method for estimating a group canopy configuration index comprises the following steps:

[0019] S1. Obtain multispectral, RGB images and LIDAR point cloud data of the target area;

[0020] S2. Perform multispectral reconstruction to generate a multispectral raster image of the entire target area;

[0021] S3. Perform LIDAR point cloud reconstruction to generate a three-dimensional point cloud and digital surface model DSM of the entire target area;

[0022] S4. reconstruct the RGB image to generate a dense point cloud of the entire target area;

[0023] S5. Fusing the LIDAR point cloud and the dense point cloud generated by the RGB image;

[0024] S6. Calculate 15 vegetation indices based on multispectral raster images;

[0025] S7. Extracting the corn canopy height model CHM according to the DSM digital surface model image, the calculation formula of the CHM is:

[0026] CHM=DSM-DEM

[0027] Among them, DEM is digital elevation model;

[0028] S8. Extract the point cloud density and geometric neighborhood features of the fused point cloud, using two density representation forms: surface density and volume density. The calculation formula of the surface density is:

[0029] SD=N / (π·R 2 )

[0030] The calculation formula of the volume density is:

[0031]

[0032] Among them, N is the number of neighbors, R is the neighborhood radius, and R is set to 0.05;

[0033] S9. Based on vegetation index, point cloud density and point cloud geometric neighborhood characteristics, the estimation models of LAI, Kurtosis and Skewness are constructed, and the CHM model is used to estimate PH;

[0034] S10. Using the determination coefficient R 2 , root mean square error RMSE and mean absolute error MAE to evaluate the estimation performance of the model;

[0035] The R 2 , RMSE, and MAE calculation formulas are:

[0036]

[0037] Where N is the number of samples, y i and are the measured parameters and the predicted parameter values, respectively. is the average value of the measured parameters.

[0038] In the above-mentioned estimation method of group canopy configuration indicators, in step S1, the flight altitude of the drone multispectral, RGB and LIDAR is 12 meters, the heading overlap rate and the lateral overlap rate are 80% and 70% respectively, the image acquisition method is oblique photography, and before performing the multispectral data acquisition task, standard reflectors with reflectivities of 25%, 50% and 75% are manually photographed, the number of LIDAR echoes is 3 echoes, and the scanning method is repeated scanning.

[0039] In the above-mentioned method for constructing and estimating group canopy configuration indicators, in step S2, multispectral reconstruction is performed using DJI Terra software. Before multispectral reconstruction, calibration plate data is imported for radiation correction, and the “Band Combination” tool in ENVI software is used to construct a multispectral raster image.

[0040] In the above-mentioned method for constructing and estimating group canopy configuration indicators, in step S3, DJI Terra software is used to reconstruct LIDAR point cloud, the point cloud density is set to 100%, accuracy optimization is turned on, and the generated point cloud data is imported into ArcGIS software to construct DSM.

[0041] In the above-mentioned group canopy configuration index construction and estimation method, in step S4, Agisoft Metashape Professional software is used for RGB reconstruction, which mainly includes four steps: loading source data, photo alignment, camera calibration, and generating dense point cloud.

[0042] In the above-mentioned group canopy configuration index construction and estimation method, in step S5, real-time dynamic positioning RTK is used to establish ground control points gcp in the same coordinate system to achieve geographic alignment of LIDAR point cloud and dense point cloud generated by RGB, Context Capture software is used to fuse the two data sets, and LiDAR360 software is used to denoise, filter and crop the point cloud.

[0043] Compared with the existing technology, the advantages of the present invention are: the present invention constructs a new canopy leaf vertical distribution index according to the cumulative LAI values ​​at different heights; then fuses the drone multispectral, RGB and LIDAR data to estimate the canopy LAI, plant height and constructed indicators, thereby realizing rapid identification of corn plant type. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 This is a flow chart of a method for constructing and estimating a group canopy configuration index proposed in the present invention.

[0045] Figure 2 This is the canopy LAI vertical distribution feature extraction diagram in the present invention.

[0046] Figure 3 This is a graph of the LAI, kurtosis, skewness and PH estimation results based on drone data in the present invention.

[0047] Figure 4 This is a canopy characteristic diagram of different yield levels under the density of D4 (7.5×104 plants / hectare) in the present invention. DETAILED DESCRIPTION

[0048] Reference Figure 1-4 In order to enable those skilled in the art to better understand the technical solution of the present invention, the present invention is further described in detail below in conjunction with specific embodiments.

[0049] The present invention discloses a group canopy configuration index, and a construction method thereof comprises the following steps:

[0050] S1. Select 10 corn plants with the same growth in each plot and measure their height. Then, divide them into 5 layers according to the pH value of corn and collect leaves from each layer separately ( Figure 1 The leaf area of ​​each layer was measured in the laboratory using a photographic leaf area meter, and the LAI was calculated based on the planting density and plot area. Data were collected twice during the entire growth period: silking stage (R1) and wax maturity stage (R5).

[0051] S2. Use the "Logistic" equation to fit the cumulative LAI values ​​of different leaf layers. The Logistic equation formula is:

[0052]

[0053] Among them, H is the relative height, M is the theoretical maximum LAI value, k is the initial value parameter, and b is the growth rate parameter.

[0054] S3. Calculate the first-order derivative of the fitting curve to obtain the canopy LAD distribution curve. The first-order derivative calculation formula is:

[0055]

[0056] Where d is the difference.

[0057] S4. Extract the kurtosis and skewness indicators of the LAD curve. The greater the kurtosis, the greater the density of leaves at a certain height in the canopy. Similarly, the greater the skewness, the closer the densest position of leaves is to the lower layer. The calculation formulas for the kurtosis and skewness are:

[0058]

[0059] Where μ represents the mean value of LAD and σ represents the standard deviation of LAD.

[0060] S5. Use DJI Mavic 3 multispectral to obtain multispectral and RGB images of the target area. The flight time is from 10:00 to 14:00 in clear weather without wind and cloud. The flight altitude of the drone is 12 meters, the heading overlap rate and the side overlap rate are 80% and 70% respectively, and the image acquisition method is oblique photography. Before performing the multispectral data acquisition task, manually shoot standard reflectors with reflectivity of 25%, 50% and 75%. Use DJI M300RTK equipped with Zenmuse L1 laser radar to obtain LIDAR point cloud data of the target area, the number of echoes is 3 echoes, and the scanning method is repeated scanning. The flight altitude, overlap rate, and acquisition method are consistent with multispectral.

[0061] S6. Use DJI Terra software to perform multispectral reconstruction. Before multispectral reconstruction, import calibration plate data for radiation correction. Use the "Band Combination" tool in ENVI software to construct a multispectral raster image.

[0062] S7. Use DJI Terra software to reconstruct LIDAR point cloud. Set the point cloud density to 100% and turn on accuracy optimization. Import the generated point cloud data into ArcGIS software to construct DSM.

[0063] S8. RGB reconstruction was performed using Agisoft Metashape Professional software, which mainly included four steps: loading source data, photo alignment, camera calibration, and generating dense point clouds.

[0064] S9. Use real-time kinematic positioning (RTK) to establish ground control points (GCPs) in the same coordinate system to achieve geo-registration of LIDAR point cloud and dense point cloud generated by RGB. Use Context Capture software to fuse the two datasets. Use LiDAR360 software to denoise, filter and crop the point cloud.

[0065] S10. Calculate 15 vegetation indices based on the multispectral raster image. The calculation formula of the vegetation index is:

[0066]

[0067] NIR, RedE, R, and G represent the reflectance of near-infrared, red-edge, red, and green bands, respectively.

[0068] S11. Extract corn canopy height model (CHM) based on DSM image. The calculation formula of CHM is:

[0069] CHM=DSM-DEM

[0070] Among them, DEM is a digital elevation model.

[0071] S12. Extract the point cloud density and geometric neighborhood features of the fused point cloud. Use two density representation forms: surface density and volume density. The calculation formula of the surface density is:

[0072] SD=N / (π·R 2 )

[0073] The calculation formula of the volume density is:

[0074]

[0075] Among them, N is the number of neighbors, R is the neighborhood radius, and R is set to 0.05.

[0076] The point cloud geometric neighborhood features are:

[0077]

[0078] S13. Based on vegetation index, point cloud density and point cloud geometric neighborhood characteristics, the estimation models of LAI, Kurtosis and Skewness are constructed; and the CHM model is used to estimate PH.

[0079] S14. Using the coefficient of determination (R 2 ), root mean square error (RMSE) and mean absolute error (MAE) to evaluate the estimation performance of the model. 2 , RMSE, and MAE calculation formulas are:

[0080]

[0081]

[0082] Where N is the number of samples, y i and are the measured parameters and the predicted parameter values, respectively. is the average value of the measured parameters.

[0083] Figure 2 a shows the average of the cumulative LAI values ​​of all cells, and the initial LAI value is set to 0. Figure 2 b shows the fitted R 2 The distribution of R1 is 2 The range was 0.984 to 0.999, with an average of 0.993; the R 2 The range is 0.970 to 0.992, with an average value of 0.991, showing a very high goodness of fit. After differentiating the fitting curve, the canopy LAD distribution curve is obtained, as shown in Figure 2c. The peak of the LAD curve represents the area with the densest leaf area distribution inside the canopy. As can be seen from the figure, this leaf layer is located near the lower two-fifths of the plant, which is consistent with the actual situation. The area under the LAD curve represents the total LAI value. The larger the area under the curve, the larger the LAI value. Finally, the kurtosis and skewness characteristics of different curves were calculated, and their distribution is shown in Figure 2 d. The kurtosis and skewness of R1 are slightly different, with standard deviations of 0.075 and 0.059, respectively; while the kurtosis and skewness of R5 are significantly different, with standard deviations of 0.268 and 0.155, respectively.

[0084] Figure 3 The estimation results of the four canopy configuration indices measured manually are shown in Figure 2. The ten-fold cross validation results show that all parameters have achieved high estimation accuracy, and R 2 All of them were higher than 0.79, with the highest reaching 0.929. 2 In the R1 and R5 periods, the predicted values ​​were 0.881 and 0.899, respectively, and the RMSE and MAE were both within 0.25 and 0.2. However, when the LAI was large, the predicted values ​​were often larger than the measured values. This may be because the canopy was more lush when the LAI was large, causing the UAV data such as vegetation index to reach saturation, which was also the reason for the R1 error in estimating LAI. 2 One of the reasons why it is lower than R5 period. Benefiting from the high precision of LIDAR active sensor, the PH estimation shows a higher accuracy, R 2 The estimation accuracy of kurtosis and skewness was consistent in the two growth stages. The estimation accuracy of R5 stage was higher than that of R1 stage. 2 They are 0.881 and 0.901 respectively.

[0085] There was no significant difference in LAI between different yield levels. The average LAI values ​​of low-yield, medium-yield, and high-yield corn groups were 5.153, 5.292, and 5.114 in R1, and 3.417, 3.668, and 3.619 in R5, respectively. The difference between R1 and R5 was significantly different in low-yield and high-yield corn (p<0.001), indicating that high-yield varieties of corn have stronger stay-green properties and can produce more photosynthetic products.

[0086] Theoretically, a larger kurtosis indicates that the density of leaves at a certain height is higher, which may be detrimental to the photosynthesis of the lower leaves; a smaller kurtosis may limit the projected area, indicating a lower LAI. The larger the skewness, the closer the height of the maximum LAD is to the root, indicating that the middle and upper leaves can obtain more light. Therefore, a larger kurtosis combined with a larger skewness has a higher high-yield potential and is an important quantitative indicator of the ideal plant type of corn. In the present invention, the kurtosis and skewness of the high-yield group were significantly higher than those of the low-yield group (p<0.001), which is consistent with the theory.

[0087] Plant height has a decisive influence on the light energy interception capacity and light energy utilization rate of the corn canopy. Generally, a larger plant height is more likely to achieve a high yield. At a lower planting density, the yield is significantly positively correlated with plant height (p<0.001), while at a high density, the correlation between plant height and yield is not significant. With the increase of planting density, plant height usually increases, but when it exceeds the appropriate density, plant height will decrease. In addition, with the increase of planting density, the plant height and ear height will increase, but the stem diameter will gradually become thinner, the lodging resistance will decrease, and the empty stalk rate will increase, which will ultimately affect the corn yield.

[0088] Please refer to Figure 4 ,D4 density under different yield levels of canopy structure characteristics change diagram,The test results show that the embodiment of the present invention proposes a group canopy configuration index that is close to the actual situation and can reflect the actual field corn canopy situation. In addition, the use of drones can quickly and accurately estimate the four canopy configuration indicators, improving the efficiency of plant type identification.

[0089] It is known from common technical knowledge that the present invention can be implemented by other embodiments that do not deviate from its spirit or essential features. Therefore, the above disclosed embodiments are only illustrative in all respects and are not exclusive. All changes within the scope of the present invention or within the scope equivalent to the present invention are included in the present invention.

Claims

1. A group canopy configuration index, characterized in that: The construction method includes the following steps: S1. Manually measure the plant height pH of each corn plot, divide it into 5 layers according to the plant height, and obtain the leaf area index LAI value of the entire canopy and leaves in each layer; S2. Use the Logistic equation to fit the cumulative LAI values ​​of different leaf layers. The Logistic equation formula is: Among them, H is the relative height, M is the theoretical maximum LAI value, k is the initial value parameter, and b is the growth rate parameter; S3. Calculate the first-order derivative of the fitting curve to obtain the canopy leaf area density LAD distribution curve. The first-order derivative calculation formula is: Where, d is the difference; S4. Extract the kurtosis and skewness index of the LAD curve. The calculation formula of the kurtosis and skewness is: Where μ represents the mean value of LAD and σ represents the standard deviation of LAD.

2. A method for constructing and estimating a group canopy configuration index according to claim 1, characterized in that: In step S1, multiple corn plants with uniform growth are selected for height measurement, and the leaf area of ​​each layer is measured using a photographic leaf area meter in a laboratory, and the LAI is calculated based on the planting density and the plot area.

3. A method for estimating a group canopy configuration index, characterized in that: The configuration index of claim 1 is estimated, the steps comprising: S1. Obtain multispectral, RGB images and LIDAR point cloud data of the target area; S2. Perform multispectral reconstruction to generate a multispectral raster image of the entire target area; S3. Perform LIDAR point cloud reconstruction to generate a three-dimensional point cloud and digital surface model DSM of the entire target area; S4. reconstruct the RGB image to generate a dense point cloud of the entire target area; S5. Fusing the LIDAR point cloud and the dense point cloud generated by the RGB image; S6. Calculate 15 vegetation indices based on multispectral raster images; S7. Extracting the corn canopy height model CHM according to the DSM digital surface model image, the calculation formula of the CHM is: CHM=DSM-DEM Among them, DEM is digital elevation model; S8. Extract the point cloud density and geometric neighborhood features of the fused point cloud, using two density representation forms: surface density and volume density. The calculation formula of the surface density is: SD=N / (π·R 2 ) The calculation formula of the volume density is: Among them, N is the number of neighbors, R is the neighborhood radius, and R is set to 0.05; S9. Based on vegetation index, point cloud density and point cloud geometric neighborhood characteristics, the estimation models of leaf area index LAI, Kurtosis and Skewness are constructed, and the CHM model is used to estimate PH; S10. Using the determination coefficient R 2 , root mean square error RMSE and mean absolute error MAE to evaluate the estimation performance of the model; The R 2 , RMSE, and MAE calculation formulas are: Where N is the number of samples, y i and are the measured parameters and the predicted parameter values, respectively. is the average value of the measured parameters.

4. A method for constructing and estimating a group canopy configuration index according to claim 3, characterized in that: In step S1, the flight altitude of the UAV multispectral, RGB and LIDAR is 12 meters, the heading overlap rate and the lateral overlap rate are 80% and 70% respectively, the image acquisition method is oblique photography, and before performing the multispectral data acquisition task, standard reflectors with reflectivity of 25%, 50% and 75% are manually photographed. The number of LIDAR echoes is 3 echoes, and the scanning mode is repeated scanning.

5. The method for constructing and estimating a group canopy configuration index according to claim 3, characterized in that: In step S2, multispectral reconstruction is performed using software. Before multispectral reconstruction, the calibration plate data is imported for radiation correction, and the "Band Combination" tool in ENVI software is used to construct a multispectral raster image.

6. A method for estimating group canopy configuration indexes according to claim 3, characterized in that: In step S3, the LIDAR point cloud is reconstructed using software, the point cloud density is set to 100%, the accuracy optimization is turned on, and the generated point cloud data is imported into ArcGIS software to construct the DSM.

7. A method for constructing and estimating a group canopy configuration index according to claim 3, characterized in that: In step S4, software is used to perform RGB reconstruction, which mainly includes four steps: loading source data, photo alignment, camera calibration, and generating dense point cloud.

8. A method for constructing and estimating a group canopy configuration index according to claim 3, characterized in that: In step S5, real-time dynamic positioning RTK is used to establish ground control points gcp in the same coordinate system to achieve geo-registration of the LIDAR point cloud and the dense point cloud generated by RGB, and the two data sets are fused using software, and the point cloud is denoised, filtered and cropped using software.

Images