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

Abstract

A crop canopy height estimation method based on laser radar and multispectral fusion, through laser radar scanning of crop fields, point cloud data of the fields are obtained; through aerial photography of the crop fields by an airborne multispectral camera, multispectral image data of the fields are obtained, after data image splicing, crop area extraction and plot segmentation, two strongest point cloud feature parameters related to crop height are extracted from the point cloud data, one strongest vegetation index related to crop height is extracted from the multispectral image data, a height estimation model is established through multivariate linear regression solution. The present application combines multispectral vegetation index and laser radar data to improve the precision of crop canopy height monitoring, can improve the accuracy and reliability of crop growth monitoring, and can comprehensively understand the growth status of crops, so that agricultural decision makers can make more accurate, timely and effective farmland management decisions based on crop canopy structure data.

Description

Technical Field

[0001] This invention relates to the field of crop canopy height estimation methods, and more particularly to a crop canopy height estimation method based on lidar and multispectral fusion. Background Technology

[0002] LiDAR point cloud features and spectral vegetation indices are commonly used to monitor crop canopy height. However, for low-growing crops, the short round-trip signal of the radar can significantly impact the accuracy of the detection results. When monitoring and studying low-growing crops, radar data may contain invalid or noisy points due to crop characteristics and environmental factors. Furthermore, when the point cloud density is low, it may be impossible to accurately capture crop point cloud data. When crops reach height saturation or form multi-layered canopy structures, vegetation indices are not sensitive enough to small changes in crop height, limiting their application in specific situations. Therefore, point cloud feature-based models are susceptible to point density distribution and noise, while vegetation indices can become saturated due to high vegetation density, affecting canopy height monitoring. By integrating point cloud features and spectral vegetation indices to construct a canopy height estimation model, their complementary advantages can be effectively addressed, resolving the insufficient accuracy of estimating crop canopy height using either point cloud features or vegetation indices alone. This provides a more accurate processing method for remote sensing monitoring of phenotypic information of low-growing, densely planted crops. Summary of the Invention

[0003] The purpose of this invention is to provide a method for estimating crop canopy height based on lidar and multispectral fusion.

[0004] The technical solution adopted by this invention to solve the above-mentioned technical problems is: a method for estimating crop canopy height based on lidar and multispectral fusion, comprising the following steps:

[0005] Step 1: Scan the crop field with lidar to obtain field ground cloud data; take aerial photos of the crop field with an airborne multispectral camera to obtain multispectral image data of the field.

[0006] 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.

[0007] Image control points are extracted from the stitched field point cloud data. The visible light image in the stitched field multispectral image is matched with the point cloud image using the point processing method to obtain the field visible light image registered with the field point cloud data.

[0008] Using the registered field visible light image as the reference image, the reference image and the field multispectral image to be calibrated are spatially aligned to obtain the field multispectral image that is spatially registered with the field ground cloud data.

[0009] Step 3: Classify the field point cloud data after step 2 by using the RGB vegetation index threshold method, extract the crop area, and obtain the crop point cloud data.

[0010] Step 4: Perform a first clustering of the crop point cloud data using the K-means algorithm to divide the crop field into multiple regions. Then, use a height threshold to segment the first clustered point cloud data into upper and lower layers, retaining the upper layer data as canopy point cloud data. Perform a second 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 outline point cloud. Obtain the segmentation distance value of each crop field from the outline point cloud, and segment the crop point cloud data into plots based on the segmentation distance value.

[0011] Step 5: Select multiple vegetation indices and set threshold ranges for each. Based on the threshold ranges of the multiple vegetation indices, classify the pixels of the field multispectral image registered in Step 2, extract the crop region, obtain crop multispectral image data, and then perform plot segmentation on the crop multispectral image data according to the segmentation distance value in Step 4.

[0012] Step 6: Extract point cloud feature parameters related to crop height from the segmented crop point cloud data, and extract vegetation indices related to crop height from the segmented crop multispectral image data;

[0013] 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 that are most strongly correlated with the actual height of the crop.

[0014] Calculate the Pearson correlation coefficient between each vegetation index and the actual crop height, and select the vegetation index with the strongest correlation to the actual crop height.

[0015] Step 7: Input all the data of the two point cloud feature parameters and one vegetation index selected in Step 6 into the Python software, and also input the actual crop height data into the Python software. Divide the dataset of all 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, and use the crop height as the target variable of multiple linear regression. Solve the multiple linear regression using the data of the modeling set to establish a height estimation model.

[0016] The accuracy of the height estimation model is validated using data from the validation set. If the accuracy is not up to standard, the modeling set and validation set are re-divided and multiple linear regression is performed again until the accuracy is up to standard, thus obtaining the crop canopy height estimation model. By establishing crop canopy height estimation models for different crops and different production stages, the crop canopy height is estimated.

[0017] Preferably, in step two, the iterative linear regression rotation method is used to perform angular rotation correction on the registered multi-spectral images of the field.

[0018] 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.

[0019] Preferably, in step five, three vegetation indices, namely the normalized difference 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.

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

[0021] Preferably, in step seven, the method for solving multiple linear regression is to represent the multiple linear regression model as:

[0022] y = β0 + β1x1 + β2x2 + β3x3;

[0023] In the above formula, y is the target variable, x1, x2, and x3 are input variables, β0 is the intercept term, and β1, β2, and β3 are regression coefficients;

[0024] The objective function to be minimized is:

[0025] ;

[0026] Then, the optimal solutions of the regression coefficients β1, β2, and β3 are respectively derived by the least squares method to minimize the sum of the squared residuals of the predicted values and the true values.

[0027] According to the above technical solution, the beneficial effects of the present invention are:

[0028] 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 the data of lidar and multi-spectral fusion. 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 growth status of crops, 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

[0029] This embodiment provides a method for estimating crop canopy height based on lidar and multispectral fusion, including the following steps:

[0030] Step 1: Scan the crop field with lidar to obtain field ground cloud data; take aerial photos of the crop field with an airborne multispectral camera to obtain multispectral image data of the field.

[0031] This embodiment uses the AA10 multi-rotor UAV laser aerial survey system (Huace Navigation Technology Co., Ltd., Shanghai) to acquire wheat canopy lidar point cloud data and orthophoto data. A DJI Phantom 4 RTK is used to acquire five-band (blue, green, red, red-edge, and near-infrared) multispectral image data.

[0032] 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.

[0033] In this embodiment, the lidar point cloud data is processed using the accompanying software Copre2 to perform POS calculation and achieve point cloud data stitching; the multispectral data is processed using DJI Terra software, and radiometric correction is performed using standard reflectivity correction plates of 25%, 50%, and 75% to achieve multispectral image stitching.

[0034] Image control points are extracted from the stitched field point cloud data. The visible light image in the stitched field multispectral image is matched with the point cloud image using the point piercing processing method to obtain the field visible light image registered with the field point cloud data.

[0035] In this embodiment, four white markers are placed at the boundary of the land parcel, arranged in specific layouts such as "□" and "X" shapes on the ground. These marker frames are used to extract control points (VCS) from the LiDAR point cloud data. Four VCS are extracted from each marker, for a total of 16 VCS, ensuring uniform spatial coverage. The selected 16 VCS (XYZ coordinates) are matched with their corresponding ground locations on the image using DJI Terra software for aerial triangulation. Then, a point-splitting method is used to match the multispectral image with the point cloud image, obtaining a visible light image precisely registered with the LiDAR point cloud data.

[0036] Using the registered field visible light image as the reference image, the reference image and the field multispectral image to be calibrated are spatially aligned to obtain a field multispectral image that is spatially registered with the field ground cloud data. Specifically, the image registration function of ENVI software is used to perform fine spatial alignment processing between the reference image and the image to be calibrated.

[0037] Due to the geographical location of the experimental field itself, the obtained remote sensing images deviate from the true 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 fitting 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 iteration processes.

[0038] Step 3: Classify the spliced point cloud data of the field through the RGB vegetation index threshold method, extract the crop area from it, 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. By utilizing the absorption characteristics of plants, which absorb strongly in the red light, weakly in the 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.

[0039] Step 4: Perform a primary clustering on the crop point cloud data through the K-means algorithm to divide the crop field into multiple regions. Then, perform upper and lower layer segmentation on the point cloud data of the primary clustering through a height threshold, and retain the upper layer data as the canopy point cloud data. Perform a secondary clustering on the canopy point cloud data through the DBSCAN algorithm, draw the boundary box of the clustered vegetation point cloud, and 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.

[0040] 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 from it, 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.

[0041] 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 meet at least two threshold ranges are classified as vegetation.

[0042] 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.

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

[0044] The vegetation indices selected in this embodiment that are correlated with crop height 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), and Modified Chlorophyll Absorption Ratio Index (MCARI).

[0045] Then, the Pearson correlation coefficient was used as an evaluation index to analyze the correlation between point cloud feature parameters and vegetation indices and measured height. 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 represents the strength of the correlation, while the sign indicates the direction of the correlation. The Pearson correlation coefficient, usually denoted as r, ranges from -1 to +1. A value close to +1 indicates a strong positive correlation, a value close to -1 indicates a strong negative correlation, and a value close to 0 indicates no linear correlation.

[0046] 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 that have the strongest correlation with the actual height of the crop.

[0047] Calculate the Pearson correlation coefficient between each vegetation index and the actual crop height, and select the vegetation index with the strongest correlation to the actual crop height.

[0048] 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.

[0049] Table 1. Selection results of point cloud feature parameters and vegetation indices for wheat at different growth stages.

[0050]

[0051] In Table 1, Rough represents roughness, Den represents density, SDHD represents standard deviation of height difference, MP represents average flatness, 3DVol represents three-dimensional volume, Skew represents skewness, and MHD represents average height difference.

[0052] In Table 1, MCARI is the improved chlorophyll uptake ratio index, and WDRVI is the weighted differential red light and infrared vegetation index.

[0053] Step 7: Input all data from the two point cloud feature parameters and one vegetation index selected in Step 6 into Python software, and also input the actual crop height data into Python software. Divide the dataset into a modeling set and a validation set, with a data volume ratio of 7:3. Use the two point cloud feature parameters and one vegetation index as the three input variables for multiple linear regression, and crop height as the target variable. Solve the multiple linear regression using the data from the modeling set to establish a height estimation model.

[0054] The accuracy of the height estimation model was validated using data from the validation set. If the accuracy was unsatisfactory, the modeling set and validation set were re-split, and multiple linear regression was performed again until the accuracy was satisfactory, thus obtaining the crop canopy height estimation model. By establishing crop canopy height estimation models for different crops and different production stages, crop canopy height was estimated.

[0055] The solution method for multiple linear regression is as follows: the multiple linear regression model is expressed as:

[0056] y = β0 + β1x1 + β2x2 + β3x3;

[0057] In the above formula, y is the target variable, x1, x2, and x3 are the input variables, β0 is the intercept term, and β1, β2, and β3 are the regression coefficients;

[0058] The objective function to be minimized is:

[0059] ;

[0060] Then, the optimal solutions for the regression coefficients β1, β2, and β3 are derived using the least squares method, minimizing the sum of squared residuals between the predicted and actual values:

[0061] ;

[0062] Where X is an n×(p+1) design matrix, and y is an n×1 target vector. It is the estimated coefficient vector.

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

[0064] Table 2. Wheat canopy height estimation model integrating point cloud features and vegetation index.

[0065]

[0066] In addition, a comparative example is given, in which only two point cloud feature parameters are selected and vegetation index is not selected. A wheat canopy height estimation model based on point cloud feature parameters is established, as shown in Table 3.

[0067] Table 3. Wheat canopy height estimation model based on point cloud feature parameters

[0068]

[0069] A comprehensive analysis of the model evaluation was conducted using statistical indicators such as the coefficient of determination and root mean square error. The coefficient of determination R0 was... 2 The root mean square error (RMSE) reflects the model's ability to explain the variability of observed data; the closer the RMSE value is to 1, the stronger the model's explanatory power. The RMSE measures the magnitude of the deviation between predicted and observed values; a lower RMSE value indicates higher accuracy in the model's predictions. In addition to the above comparisons, several other comparisons were set up. Some comparisons changed the method of establishing the estimation model, using random forest and partial least squares regression methods. Other comparisons changed the number of point cloud features and vegetation indices selected, using more point cloud features and vegetation indices for multiple linear regression. Multiple different canopy height estimation models were established using these different comparisons.

[0070] The determination coefficients and root mean square errors of 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. 2The analysis results show that the wheat canopy height estimation model obtained by the multiple linear regression method in this embodiment, which uses two point cloud features and one vegetation index, has the highest estimation accuracy. Its accuracy is higher than models established by other methods, and also higher than the estimation model in the comparative example that only uses two point cloud features. Furthermore, it is superior to height estimation models that use more point cloud features or vegetation indices in terms of estimation accuracy and ease of implementation. This invention integrates the three-dimensional spatial information of point clouds with the physiological state information of vegetation indices, which not only improves the accuracy and reliability of crop growth monitoring, but also enables a more comprehensive understanding of crop growth status, optimizing agricultural management decisions such as irrigation and fertilization, thereby improving crop yield and quality.

Claims

1. A crop canopy height estimation method based on laser radar and multispectral fusion, characterized in that, The method comprises the following steps: Step one, scanning the crop field by laser radar to obtain the field point cloud data; taking aerial photography of the crop field by the airborne multispectral camera to obtain the field multispectral image data; Step two, POS solving is performed on the field point cloud data to realize point cloud data splicing; radiation correction is performed on the field multispectral image data to realize multispectral image splicing; Extracting the image control points from the spliced field point cloud data, and matching the visible light image in the spliced field multispectral image with the point cloud image through the pin-point processing method to obtain the field visible light image registered with the field point cloud data; Taking the registered field visible light image as the reference image, performing spatial alignment processing on the reference image and the field multispectral image to be calibrated to obtain the field multispectral image registered with the field point cloud data in the spatial position; Step three, classifying the field point cloud data spliced in step two through the RGB vegetation index threshold method to extract the crop region and obtain the crop point cloud data; Step four, performing primary clustering on the crop point cloud data through the K-means algorithm, dividing the crop field into multiple regions, then performing upper and lower layer segmentation on the primary clustered point cloud data through the height threshold, retaining the upper layer data as the canopy point cloud data; performing secondary clustering on the canopy point cloud data through the DBSCAN algorithm, and drawing the boundary box of the clustered vegetation point cloud to obtain the contour point cloud; obtaining the segmentation distance value of each crop field from the contour point cloud, and performing land segmentation on the crop point cloud data according to the segmentation distance value; Step five, selecting multiple vegetation indexes and setting the threshold range respectively, classifying the pixel points of the field multispectral image registered in step two according to the threshold range of the multiple vegetation indexes, extracting the crop region to obtain the crop multispectral image data, and then performing land segmentation on the crop multispectral image data according to the segmentation distance value in step four; Step six, extracting the point cloud feature parameters related to the crop height from the segmented crop point cloud data, and extracting the vegetation indexes related to the crop height from the segmented crop multispectral image data; Calculating the Pearson correlation coefficient of each point cloud feature parameter and the actual height of the crop, and selecting the two point cloud feature parameters with the strongest correlation with the actual height of the crop; Calculating the Pearson correlation coefficient of each vegetation index and the actual height of the crop, and selecting the vegetation index with the strongest correlation with the actual height of the crop; Step seven, inputting all the data of the two point cloud feature parameters and the vegetation index selected in step six into the python software, and inputting the actual height data of the crop into the python software, dividing the data set of all the data into a modeling set and a verification set, taking the two point cloud feature parameters and the vegetation index as three input variables of the multiple linear regression, taking the crop height as the target variable of the multiple linear regression, performing multiple linear regression solving through the data of the modeling set, and establishing a height estimation model. The accuracy of the height estimation model is verified by the data of the verification set, and when the accuracy is unqualified, the modeling set and the verification set are re-divided and multiple linear regression is solved again until the accuracy is qualified, that is, the crop canopy height estimation model is obtained; the crop canopy height estimation model of different crops and different production stages is established to estimate the crop canopy height.

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

3. The crop canopy height estimation method based on laser radar and multispectral fusion according to claim 1, characterized in that: The point cloud feature parameters include roughness, average flatness, density, average height difference, height difference standard deviation, skewness and three-dimensional volume, and the vegetation indexes 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 multispectral fusion according to claim 3, characterized in that: In step five, three vegetation indexes, normalized vegetation index NDVI, enhanced vegetation index EVI and normalized difference water index NDWI, are selected, and the threshold range is 0.4<NDVI≤1, 0.4<EVI≤1 and-0.5<NDWI≤1; then the pixel points meeting at least two threshold ranges are divided into vegetation.

5. The crop canopy height estimation method based on laser radar and multispectral fusion according to claim 1, characterized in that: In step seven, the data amount ratio of the modeling set and the verification set is 7:

3.

6. The crop canopy height estimation method based on laser radar and multispectral fusion according to claim 1, characterized in that: In step seven, the multiple linear regression solving method 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 and x3 are input variables, β0 is the intercept term, β1, β2 and β3 are regression coefficients; The minimum objective function is: ； Then the optimal solutions of the regression coefficients β1, β2 and β3 are derived by the least square method respectively, so that the residual sum of squares of the predicted value and the true value is minimum.

Images