Crop canopy height estimation method based on laser radar and multispectral fusion

Abstract

The invention discloses a crop canopy height estimation method based on laser radar and multispectral fusion, and the method comprises the steps: scanning a crop field through the laser radar, and obtaining the field point cloud data; aerial photography is carried out on a crop field through an airborne multispectral camera to obtain field multispectral image data, and after data image splicing, crop region extraction and land parcel segmentation, two point cloud feature parameters with the strongest correlation with crops are extracted from point cloud data to obtain a point cloud feature parameter; the method comprises the following steps of: extracting a vegetation index with the strongest correlation with the crop height from multispectral image data, solving through multiple linear regression, and establishing a height estimation model; according to the invention, the multispectral vegetation index and the laser radar data are combined, so that the crop canopy height monitoring precision is improved, the crop growth monitoring precision and reliability can be improved, the crop growth condition can be comprehensively understood, and an agricultural decision maker can determine the crop canopy height based on the crop canopy structure data. And a more accurate, timely and effective farmland management decision can be made.

Description

Technical Field

[0001] The present invention relates to the field of crop canopy height estimation methods, and in particular to a crop canopy height estimation method based on laser radar and multispectral fusion. Background Art

[0002] LiDAR point cloud features and spectral vegetation indices are commonly used to monitor crop canopy height. For short crops, the short round-trip time of the radar signal can significantly affect the accuracy of detection results. When monitoring short crops, radar data can contain invalid or noisy points due to crop characteristics and environmental factors. Furthermore, when the collected point cloud density is low, the crop point cloud data may not be accurately captured. When crops reach saturation height or form a multi-layered canopy structure, vegetation indices are less sensitive to small changes in crop height, limiting their application in certain situations. Therefore, models based on point cloud features are susceptible to point density distribution and noise, while vegetation indices can become saturated due to high vegetation density, affecting canopy height monitoring. Combining point cloud features with spectral vegetation indices to construct a canopy height estimation model complements these two data types. This approach effectively addresses the inaccuracy of crop canopy height estimation using either point cloud features or vegetation indices alone, providing a more precise method for remote sensing phenotypic monitoring of short, densely planted crops. Summary of the Invention

[0003] The purpose of this invention is to provide a crop canopy height estimation method based on lidar and multispectral fusion. The technical solution adopted by the present invention to solve the above technical problems is: a crop canopy height estimation method based on lidar and multispectral fusion, comprising the following steps: Step 1: Scan the crop field with a lidar to obtain field cloud data; use an airborne multispectral camera to take aerial photos of the crop field to obtain field multispectral image data; Step 2: Perform POS calculation on the field point cloud data to achieve point cloud data stitching; perform radiometric correction on the field multispectral image data to achieve multispectral image stitching; Image control points are extracted from the spliced ​​field point cloud data, and the visible light image in the spliced ​​field multispectral image is matched with the point cloud image through the puncture point processing method to obtain the field visible light image registered with the field point cloud data; Using the registered visible light image of the field as the reference image, the reference image is spatially aligned with the multispectral image of the field to be calibrated to obtain a multispectral image of the field that is spatially registered with the field cloud data. Step 3: Classify the field point cloud data obtained by step 2 using the RGB vegetation index threshold method, extract the crop area, and obtain the crop point cloud data. Step 4: Cluster the crop point cloud data using the K-means algorithm to divide the crop field into multiple regions. Then, use the height threshold to segment the primary clustered point cloud data into upper and lower layers, retaining the upper layer data as the canopy point cloud data. Perform a secondary clustering of the canopy point cloud data using the DBSCAN algorithm, and draw the bounding box of the clustered vegetation point cloud to obtain the contour point cloud. Obtain the segmentation distance value of each crop field from the contour point cloud, and segment the crop point cloud data into plots based on the segmentation distance value. Step 5: Select multiple vegetation indices and set threshold ranges for each of them. Classify the pixel points of the field multispectral image registered in step 2 according to the threshold ranges of the multiple vegetation indices, extract the crop area from it, and obtain crop multispectral image data. Then, segment the crop multispectral image data into plots according to the segmentation distance value in step 4. Step 6: extracting point cloud feature parameters highly correlated with the crop from the segmented crop point cloud data, and extracting vegetation indices highly correlated with the crop from the segmented crop multispectral image data; Calculate the Pearson correlation coefficient between each point cloud feature parameter and the actual height of the crop, and select the two point cloud feature parameters with the strongest correlation with the actual height of the crop; Calculate the Pearson correlation coefficient between each vegetation index and the actual height of the crop, and select the vegetation index with the strongest correlation with the actual height of the crop; Step 7: Input all the data of the two point cloud feature parameters and one vegetation index selected in Step 5 into the Python software, and also input the actual crop height data into the Python software. The data set of all the data is divided into a modeling set and a validation set. The two point cloud feature parameters and one vegetation index are used as the three input variables of the multiple linear regression, and the crop height is used as the target variable of the multiple linear regression. The multiple linear regression solution is solved using the data of the modeling set to establish a height estimation model; The accuracy of the height estimation model is verified by the data of the validation set. When the accuracy is not up to standard, the modeling set and validation set are re-divided and the multivariate linear regression is performed again until the accuracy is qualified. The crop canopy height estimation model is obtained; by establishing crop canopy height estimation models for different crops and different production stages, the crop canopy height is estimated.

[0004] Preferably, in step 2, an iterative linear regression rotation method is used to perform angle rotation correction on the registered field multispectral image.

[0005] Preferably, the point cloud feature parameters include roughness, mean flatness, density, mean height difference, standard deviation of height difference, skewness, and three-dimensional volume, and the vegetation indices include normalized difference vegetation index, enhanced vegetation index, soil-adjusted vegetation index, optimized soil-adjusted vegetation index, weighted difference red and infrared vegetation index, and chlorophyll absorption ratio index.

[0006] Preferably, in step five, three vegetation indices, namely the normalized vegetation index NDVI, enhanced vegetation index EVI, and normalized difference water index NDWI, are selected, and the threshold ranges are 0.4 < NDVI ≤ 1, 0.4 < EVI ≤ 1, and -0.5 < NDWI ≤ 1; then the pixel points that meet at least two threshold ranges are classified as vegetation.

[0007] Preferably, in step six, the data volume ratio of the modeling set to the validation set is 7:3.

[0008] Preferably, in step six, the method for solving multiple linear regression is to represent the multiple linear regression model as: y = β 0 + β 1x1 + β 2x2 + β 3x3; In the above formula, y is the target variable, x1, x2, and x3 are input variables, β 0 is the intercept term, β [[ID=二十五]]1, β 2, β 3 are regression coefficients; The objective function to be minimized is: ; Then, the optimal solutions of the regression coefficients β 1, β 2, β 3 are derived respectively by the least squares method, so that the residual sum of squares between the predicted value and the true value y is minimized.

[0009] According to the above technical solution, the beneficial effect of the present invention is: The present invention combines multi-spectral vegetation indices with lidar data to improve the accuracy of crop canopy height monitoring. By using point cloud and spectral data, the canopy height of low-growing and dense crops such as winter wheat can be monitored based on lidar and multi-spectral fusion data. The three-dimensional spatial information of the point cloud and the physiological state information of the vegetation index are mutually fused, which can not only improve the accuracy and reliability of crop growth monitoring, but also enable a more comprehensive understanding of the crop growth status, enabling agricultural decision-makers to make more accurate, timely, and effective farmland management decisions based on the crop canopy structure data, thereby improving crop yield and quality. Detailed implementation manners

[0010] This embodiment provides a crop canopy height estimation method based on lidar and multispectral fusion, including the following steps: Step 1: Scan the crop field with lidar to obtain field cloud data; use an airborne multispectral camera to take aerial photos of the crop field to obtain field multispectral image data.

[0011] This example uses a multi-rotor drone AA10 laser aerial survey system (Huace Navigation Technology Co., Ltd., Shanghai) to obtain wheat canopy lidar point cloud data and orthophoto data. A DJI Phantom 4 RTK drone acquires five-band (blue, green, red, red-edge, and near-infrared) multispectral image data.

[0012] Step 2: Perform POS calculation on the field point cloud data to achieve point cloud data stitching; perform radiation correction on the field multispectral image data to achieve multispectral image stitching.

[0013] In this embodiment, the supporting software Copre2 is used to perform POS solution on the lidar point cloud data to achieve point cloud data stitching; the multispectral data is processed using DJITerra software, and standard reflectivity correction plates of 25%, 50%, and 75% are selected for radiation correction to achieve multispectral image stitching.

[0014] Image control points are extracted from the spliced ​​field point cloud data. The visible light image in the spliced ​​field multispectral image is matched with the point cloud image through the puncture point processing method to obtain the field visible light image that is aligned with the field point cloud data.

[0015] In this example, four white markers were placed on the ground at the boundaries of the plots, arranged in specific configurations such as "□" and "X" shapes. These marker frames were used to extract image control points (ICPs) from the LiDAR point cloud data. Four ICPs were extracted for each marker, for a total of 16 ICPs, ensuring uniform spatial coverage. The 16 selected ICPs (XYZ coordinates) were matched to the corresponding ground locations in the imagery. Aerial triangulation was performed using DJI Terra software. Then, a point-processing method was used to match the multispectral imagery with the point cloud, resulting in a visible light image that was precisely aligned with the LiDAR point cloud data.

[0016] The registered visible light image of the field is used as the reference image, and the reference image is spatially aligned with the multispectral image of the field to be calibrated to obtain a multispectral image of the field that is spatially registered with the field point cloud data. Specifically, the inter-image registration function of ENVI software is used to perform fine spatial alignment on the reference image and the image to be calibrated.

[0017] Due to the geographical location of the experimental field itself, the obtained remote sensing images deviate from the due north-south axis, which is not conducive to the subsequent extraction of the characteristics of the experimental field plots. Therefore, in this embodiment, the iterative linear regression rotation method is used for angle rotation correction. Linear regression analysis is performed using the X and Y coordinates of the point cloud data points to obtain the best-fit curve, and then the inclination angle of the curve relative to the X-axis is further calculated. The accuracy of the angle measurement is ensured through multiple iterative processes.

[0018] Step 3: Classify the stitched field point cloud data through the RGB vegetation index threshold method, extract the crop area therefrom, and obtain the crop point cloud data. The RGB vegetation index threshold method is based on the spectral reflection characteristics of the red (Red), green (Green), and blue (Blue) color channels in the original point cloud data. Utilizing the absorption characteristics of plants, which strongly absorb red light, weakly absorb green light, and the absorption of blue light is between the two, the vegetation and non-vegetation are distinguished by evaluating the reflection intensities of these color channels.

[0019] Step 4: Perform a primary clustering on the crop point cloud data through the K-means algorithm, divide the crop field into multiple regions, and then perform upper and lower layer segmentation on the point cloud data of the primary clustering through a height threshold, retaining the upper layer data as the canopy point cloud data; perform a secondary clustering on the canopy point cloud data through the DBSCAN algorithm, and draw the bounding box of the clustered vegetation point cloud to obtain the contour point cloud; obtain the segmentation distance value of each crop field from the contour point cloud, and perform plot segmentation on the crop point cloud data according to the segmentation distance value.

[0020] Step 5: Select multiple vegetation indices and set threshold ranges respectively. Classify the pixel points of the multi-spectral image of the field registered in Step 2 according to the threshold ranges of the multiple vegetation indices, extract the crop area therefrom, and obtain the crop multi-spectral image data. Then perform plot segmentation on the crop multi-spectral image data according to the segmentation distance value in Step 4. The same segmentation distance can ensure the consistency of the multi-spectral image and the plot segmentation result in the spatial scale.

[0021] In this embodiment, three vegetation indices, namely the normalized difference vegetation index NDVI, the enhanced vegetation index EVI, and the normalized difference water index NDWI, are selected, and the threshold ranges are 0.4 < NDVI ≤ 1, 0.4 < EVI ≤ 1, -0.5 < NDWI ≤ 1; then the pixel points that satisfy at least two threshold ranges are classified as vegetation.

[0022] Step 6: Extract the point cloud feature parameters related to the crop height from the segmented crop point cloud data, and extract the vegetation indices related to the crop height from the segmented crop multi-spectral image data.

[0023] The point cloud feature parameters related to crop height selected in this embodiment include: roughness (Roughness), mean planeness (MP), density (Density, Den), mean height difference (MHD), standard deviation of height difference (SDHD), skewness (Skew), and 3D volume (3DVol) for modeling.

[0024] The vegetation indices highly correlated with crops selected in this embodiment include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Soil Adjusted Vegetation Index (SAVI), Optimized Soil Adjusted Vegetation Index (OSAVI), Weighted Difference in Red and Infrared Vegetation Indices (WDRVI), Modified Chlorophyll Absorption Ratio Index (MCARI), and other parameters.

[0025] The Pearson correlation coefficient (Pearson's correlation coefficient) was then used as an evaluation metric to analyze the correlation between point cloud feature parameters, vegetation indices, and measured heights. The Pearson correlation coefficient is a statistic used to measure the degree of linear correlation between two variables. The absolute value of the correlation coefficient indicates the strength of the correlation, while the sign indicates the direction of the correlation. The Pearson's correlation coefficient (r), commonly expressed as r, ranges from -1 to +1. Values ​​close to +1 indicate a strong positive correlation, values ​​close to -1 indicate a strong negative correlation, and values ​​close to 0 indicate no linear correlation.

[0026] The Pearson correlation coefficient between each point cloud feature parameter and the actual height of the crop is calculated respectively, and the two point cloud feature parameters with the strongest correlation with the actual height of the crop are selected.

[0027] The Pearson correlation coefficient between each vegetation index and the actual height of the crop was calculated respectively, and the vegetation index with the strongest correlation with the actual height of the crop was selected.

[0028] The results of selecting point cloud feature parameters and vegetation indices for wheat at different growth stages in this embodiment are shown in Table 1.

[0029] Table 1 Results of point cloud characteristic parameters and vegetation index selection at different growth stages of wheat In Table 1, Rough is roughness, Den is density, SDHD is standard deviation of height difference, MP is mean flatness, 3DVol is three-dimensional volume, Skew is skewness, and MHD is mean height difference.

[0030] In Table 1, MCARI is the modified chlorophyll absorption ratio index, and WDRVI is the weighted difference red and infrared vegetation index.

[0031] Step 7: Input all data from the two point cloud feature parameters and one vegetation index selected in Step 5 into Python software, along with the actual crop height data. Divide the resulting dataset into a modeling set and a validation set, with the ratio of 7:3 between the two sets. Use the two point cloud feature parameters and one vegetation index as the three input variables for a multiple linear regression, and crop height as the target variable. Perform a multiple linear regression on the modeling set data to establish a height estimation model.

[0032] The accuracy of the height estimation model is verified using the validation set data. If the accuracy fails, the modeling set and validation set are repartitioned and the multivariate linear regression solution is repeated until the accuracy meets the requirements. This results in a crop canopy height estimation model. Crop canopy height estimation models are established for different crops and production stages to estimate crop canopy height.

[0033] The solution method of multiple linear regression is to express the multiple linear regression model as: y= β 0+ β 1x1+ β 2x2+ β 3x3; In the above formula, y is the target variable, x1, x2, x3 are input variables, β 0 is the intercept term, β 1. β 2. β 3 is the regression coefficient; The minimized objective function is: ; Then the regression coefficients were derived by the least squares method. β 1. β 2. β The optimal solution of 3 makes the predicted value The sum of squares of the residuals and the true value y is minimized: ; Where X is the n×(p+1) design matrix, y is the n×1 target vector, is the estimated coefficient vector.

[0034] The wheat canopy height estimation model obtained by fusing point cloud features and vegetation index in this embodiment is shown in Table 2.

[0035] Table 2 Wheat canopy height estimation model based on integration of point cloud features and vegetation index In addition, a comparative example is given. In this example, only two point cloud feature parameters are selected, and no vegetation index is selected. A wheat canopy height estimation model based on point cloud feature parameters is established, as shown in Table 3.

[0036] Table 3 Wheat canopy height estimation model based on point cloud feature parameters The model evaluation was comprehensively analyzed using statistical indicators such as the coefficient of determination and root mean square error. 2 This reflects the extent to which the model explains the variability in the observed data. The closer its value is to 1, the stronger the model's explanatory power. The root mean square error (RMSE) measures the deviation between the predicted value and the actual observed value. Lower RMSE values ​​indicate higher model prediction accuracy. In addition to the aforementioned comparison examples, several other comparisons were also set up. Some of these comparisons varied the estimation model establishment method, using random forests and partial least squares regression methods for solution. Other comparisons varied the number of point cloud features and vegetation indices selected, selecting a larger number of point cloud features and vegetation indices for multivariate linear regression. Using these different comparisons, multiple different canopy height estimation models were established.

[0037] The determination coefficient and root mean square error of the estimation results of the multiple linear regression, random forest and partial least squares regression models at different growth stages were compared and analyzed, and the determination coefficient R was selected. 2Analysis of models with higher values ​​(the closer to 1, the better) and their corresponding combinations shows that the wheat canopy height estimation model solved by the multivariate linear regression method in this embodiment, which uses two point cloud features and a vegetation index, has the most accurate estimation accuracy. Its accuracy is higher than that of models established by other methods and also higher than the comparison model using only two point cloud features. In terms of estimation accuracy and ease of implementation, it outperforms height estimation models that use more point cloud features or vegetation indices. This present invention integrates the three-dimensional spatial information of point clouds with the physiological status information of vegetation indices, not only improving the accuracy and reliability of crop growth monitoring but also enabling a more comprehensive understanding of crop growth conditions, optimizing agricultural management decisions such as irrigation and fertilization, thereby increasing crop yield and quality.

Claims

1. A crop canopy height estimation method based on laser radar and multi-spectral fusion, characterized in that: The following steps are involved: Step 1: Scan the crop field with a laser radar to obtain field cloud data; take aerial photos of the crop field with an airborne multispectral camera to obtain field multispectral image data; Step 2: Perform POS calculation on the field point cloud data to achieve point cloud data stitching; perform radiation correction on the field multispectral image data to achieve multispectral image stitching; Image control points are extracted from the spliced ​​field point cloud data, and the visible light image in the spliced ​​field multispectral image is matched with the point cloud image by the puncture point processing method to obtain the field visible light image registered with the field point cloud data; The registered visible light image of the field is used as the reference image, and the reference image is spatially aligned with the field multispectral image to be calibrated to obtain the field multispectral image that is spatially registered with the field point cloud data. Step 3: Classify the field point cloud data spliced ​​in step 2 by using the RGB vegetation index threshold method, extract the crop area from it, and obtain the crop point cloud data; Step 4: Perform a clustering of the crop point cloud data using the K-means algorithm to divide the crop field into multiple areas. Then, use the height threshold to segment the point cloud data of the first clustering into upper and lower layers, and retain the upper layer data as the canopy point cloud data. Perform a secondary clustering of the canopy point cloud data using the DBSCAN algorithm, and draw the clustered vegetation point cloud boundary box to obtain the contour point cloud. Obtain the segmentation distance value of each crop field from the contour point cloud, and segment the crop point cloud data into plots according to the segmentation distance value. Step 5: Select multiple vegetation indices and set threshold ranges respectively, classify the pixel points of the field multispectral image registered in step 2 according to the threshold ranges of the multiple vegetation indices, extract the crop area therefrom, obtain the crop multispectral image data, and then segment the crop multispectral image data according to the segmentation distance value of step 4; Step 6: extracting point cloud feature parameters highly correlated with the crop from the segmented crop point cloud data, and extracting vegetation indices highly correlated with the crop from the segmented crop multispectral image data; Calculate the Pearson correlation coefficient between each point cloud feature parameter and the actual height of the crop, and select the two point cloud feature parameters with the strongest correlation with the actual height of the crop; Calculate the Pearson correlation coefficient between each vegetation index and the actual height of the crop, and select the vegetation index with the strongest correlation with the actual height of the crop; Step 7: Input all the data of the two point cloud feature parameters and one vegetation index selected in step 5 into the Python software, and also input the actual crop height data into the Python software, divide the data set of all the data into a modeling set and a validation set, use the two point cloud feature parameters and one vegetation index as the three input variables of multiple linear regression, use the crop height as the target variable of multiple linear regression, perform multiple linear regression on the data of the modeling set, and establish a height estimation model; Verify the accuracy of the height estimation model using the data of the validation set. When the accuracy is unqualified, re-divide the modeling set and the validation set and perform multiple linear regression solution again until the accuracy is qualified, then the crop canopy height estimation model is obtained; estimate the crop canopy height by establishing crop canopy height estimation models for different crops and different production stages.

2. The crop canopy height estimation method based on laser radar and multi-spectral fusion according to claim 1, characterized in that: In step two, the iterative linear regression rotation method is used to correct the angle rotation of the registered field multispectral image.

3. The crop canopy height estimation method based on laser radar and multi-spectral fusion according to claim 1, characterized in that: The point cloud feature parameters include roughness, average flatness, density, average height difference, standard deviation of height difference, skewness and three-dimensional volume, and the vegetation indices include normalized difference vegetation index, enhanced vegetation index, soil-adjusted vegetation index, optimized soil-adjusted vegetation index, weighted difference red and infrared vegetation index, and chlorophyll absorption ratio index.

4. The crop canopy height estimation method based on laser radar and multi-spectral fusion according to claim 3 is characterized by: In step five, select three vegetation indices, namely the normalized vegetation index NDVI, the enhanced vegetation index EVI, and the normalized difference water index NDWI. The threshold ranges are 0.4 < NDVI ≤ 1, 0.4 < EVI ≤ 1, and -0.5 < NDWI ≤ 1; then divide the pixel points that meet at least two threshold ranges into vegetation.

5. The crop canopy height estimation method based on laser radar and multi-spectral fusion according to claim 1, characterized in that: In step six, the data volume ratio of the modeling set to the validation set is 7:

3.

6. The crop canopy height estimation method based on laser radar and multi-spectral fusion according to claim 1, characterized in that: In step six, the multiple linear regression solution method is to represent the multiple linear regression model as: y= β 0+ β 1x1+ β 2x2+ β 3x3; In the above formula, y is the target variable, x1, x2, x3 are input variables, β 0 is the intercept term, β 1. β 2. β 3 is the regression coefficient; The minimized objective function is: ； Then the regression coefficients were derived by the least squares method. β 1. β 2. β The optimal solution of 3 makes the predicted value The sum of squares of the residuals with the true value y is minimized.

Images