A method for processing and predicting multi-source remote sensing yield data in the early season of crops

Abstract

The invention discloses a method for processing and predicting multi-source remote sensing yield data in the early season of crops, wherein the data processing method comprises the following steps: preprocessing multispectral images and laser radar data to obtain a plurality of multispectral feature maps and laser radar feature maps corresponding to a target area; merging all the multispectral feature maps and laser radar feature maps into an original data matrix, performing non-negative matrix decomposition of bivariate parameters on the original data matrix to obtain a new feature weight matrix; extracting the original feature sequence number with the largest proportion of each new feature from the new feature weight matrix, screening the original data matrix to obtain the new data matrix, thereby completing the data processing.

Description

Technical Field

[0001] The invention relates to the technical field of agricultural remote sensing image processing, and in particular to a method for processing and predicting multi-source remote sensing yield data in the early season of crops. Background Art

[0002] Traditional crop yield prediction methods mainly include agronomic yield estimation, statistical yield estimation, and meteorological statistics, etc., but these methods are often only applicable to small-area yield estimation. Since the 1970s, remote sensing technology has been widely used in agricultural production due to its macroscopic, dynamic and rapid characteristics, providing an effective way to monitor crop growth and predict grain yield.

[0003] In recent years, on the one hand, sensors that obtain spectral and other structural information have gradually become smaller and more integrable, and on the other hand, the payload and stability of drone platforms equipped with these sensors have also made great progress. These changes have made the combined and synchronous data acquisition of multiple sensors increasingly one of the mainstream remote sensing detection solutions. The combination of multi-modal sensor data and its application in yield prediction can combine the advantages of rich spectral, spatial and structural information contained in different sensor systems to obtain more comprehensive and detailed current crop information, which helps to improve the prediction accuracy of yield prediction.

[0004] The existing crop remote sensing yield estimation technology can be roughly divided into three steps: estimating the crop planting area by analyzing remote sensing image data; extracting the corresponding crop vegetation index by analyzing remote sensing image data to monitor the crop growth status; constructing a yield estimation model based on vegetation index, crop yield, and other meteorological and agronomic parameters, and then obtaining the total yield through further calculation. When constructing a yield estimation model using multi-source remote sensing data features and crop yield, in order to avoid redundant features or invalid features that reduce the prediction performance of the yield estimation model, the PCA algorithm is often used to reduce the dimension of the data, and then the multivariate linear regression method is used to construct a relationship function for the remote sensing data and yield data after dimensionality reduction to obtain a simple yield prediction model. However, the obtained model has poor interpretability and insufficient prediction accuracy. Summary of the invention

[0005] In order to solve the problems of the above existing technologies, the present invention provides a method for processing and predicting multi-source remote sensing yield data in the early stage of a crop season.

[0006] In order to achieve the above-mentioned purpose of the present invention, the technical scheme adopted is as follows:

[0007] A method for processing multi-source remote sensing yield data of crop season in the early stage, comprising the following steps:

[0008] Preprocess the multispectral image and lidar data to obtain multiple multispectral feature maps and lidar feature maps corresponding to the target area;

[0009] All multispectral feature maps and lidar feature maps are combined into an original data matrix, and the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix;

[0010] The original feature number that accounts for the largest proportion of each new feature is extracted from the new feature weight matrix, and the new data matrix is ​​obtained by screening from the original data matrix, thereby completing the data processing.

[0011] Preferably, the multispectral image is preprocessed to obtain a multispectral feature map, and the specific steps are as follows:

[0012] The multiple channel images of the multispectral image are processed using N1 different vegetation indices that can reflect the growth status of crops to obtain N1 vegetation index images;

[0013] Taking the size of the biomass sampling area at the final yield harvest as the standard, N1 vegetation index maps were gridded synchronously;

[0014] Then, N2 different derivatives are extracted in each grid unit of each vegetation index map to form N1×N2 multispectral feature maps, where N1 and N2 are both positive integers.

[0015] Furthermore, the vegetation index includes N1=10 items, namely normalized vegetation index, normalized difference red edge index, green light normalized difference vegetation index, enhanced vegetation index, modified anthocyanin content index, optimized soil conditioning vegetation index, simplified canopy chlorophyll content index, transformed chlorophyll absorption and reflectance index, triangular green index, and visible atmospheric resistance index.

[0016] Furthermore, the derivatives include N2=7 items in total, including maximum value, minimum value, average value, standard value, lower quartile value, median value, and upper quartile value.

[0017] Preferably, the laser radar image is preprocessed to obtain a laser radar feature map, and the specific steps are as follows:

[0018] Taking the size of the biomass sampling area at the time of final yield harvest as the standard, the required vegetation data is extracted in each regional unit to generate multiple lidar feature maps.

[0019] Furthermore, the vegetation data include maximum height, average height, mean square height, height standard deviation, height skewness, height kurtosis, 5th to 95th height percentiles, scores from 5th to 95th from the ground to the height percentiles, and the proportion of the four height points.

[0020] Preferably, the multispectral feature map and the laser radar feature map are respectively arranged into row vectors X in the same position order. 1,n , n represents the number of samples in the feature map;

[0021] Then, the row vectors of the multispectral feature map and the laser radar feature map, a total of m feature map transformations, are merged to obtain the original data matrix X that combines the multispectral feature map and the laser radar feature map m,n .

[0022] Furthermore, the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix, as follows:

[0023] According to the original data matrix X m,n The dimension size and the number of new features set are initialized to construct a bivariate parameter matrix Q k,m and P m,k , where m is the number of feature maps, n is the number of samples in the feature map, and k is the number of new features;

[0024] In order to make the original data matrix X m,n ≈P m,k Q k,m X m,n , will solve Q k,m and P m,k The problem is transformed into a nonlinear optimization problem:

[0025]

[0026] The alternating least squares method is used to randomly initialize the parameter p, and then use p and formula (1) to calculate the parameter q. Then, use the parameter q and formula (1) to solve the parameter p. This process is repeated until it converges to the local minimum.

[0027] By continuously optimizing P m,k and Q k,m , get the best new feature matrix Q k,m X m,n , and the new feature weight matrix Q k,m .

[0028] Furthermore, a new data matrix is ​​obtained by screening the original data matrix, as follows:

[0029] Initialize an empty new data matrix Y k,n ;

[0030] For the new feature weight matrix Q k,m Traverse each row to find the data with the largest value, and record the column coordinate of the data as index;

[0031] In the original data matrix Xm,n Retrieve the data of the index row from the matrix and merge it into the new data matrix Y k,n middle;

[0032] After k traversals, the new data matrix Y after key feature extraction is obtained k,n .

[0033] A method for predicting crop yield data from multi-source remote sensing in the early season, comprising the following steps:

[0034] The new data matrix Y is obtained by the method for processing multi-source remote sensing yield data in the early season of crops. k,n ;

[0035] By collecting the yield data Z of the biomass sampling area at the corresponding final yield harvest 1，n , the new data matrix Y k,n With production data Z 1，n Input deep neural network for regression training, construct the potential relationship function between remote sensing data features and final yield, obtain the yield prediction model, and retain the input feature parameters of the yield prediction model;

[0036] The final yield prediction model is used to predict crop yield.

[0037] The beneficial effects of the present invention are as follows:

[0038] The present invention targets multi-spectral and laser radar remote sensing data of two different modes that can provide rich spectral information and crop three-dimensional structure information, performs pre-processing and extracts derivatives within the same grid division unit, and mines more prediction factors that have potential benefits for yield prediction. Then, the data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to screen out the more critical prediction factors, while eliminating the redundancy between the original features and making the new features obtained by dimensionality reduction have practical significance. Finally, a neural network is used to perform nonlinear regression fitting on the prediction factors and yield results, and the obtained yield prediction model can further improve the prediction accuracy of yield on the basis of the linear regression model. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a flow chart of the data processing method described in the present invention.

[0040] Figure 2 It is a schematic diagram of the conversion of feature map into a data matrix.

[0041] Figure 3 It is a schematic diagram of obtaining a new data matrix from a new feature weight matrix.

[0042] Figure 4 It is a schematic diagram of the prediction method of the present invention.

[0043] Figure 5 The present invention is used to obtain the sugarcane yield prediction model effect diagram.

[0044] Figure 6 This is the effect diagram of the traditional sugarcane yield prediction model. DETAILED DESCRIPTION

[0045] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0046] Example 1

[0047] like Figure 1 As shown, a method for processing multi-source remote sensing yield data of crop season in the early stage comprises the following steps:

[0048] Preprocess the multispectral image and lidar data to obtain multiple multispectral feature maps and lidar feature maps corresponding to the target area;

[0049] All multispectral feature maps and lidar feature maps are combined into an original data matrix, and the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix;

[0050] The original feature number that accounts for the largest proportion of each new feature is extracted from the new feature weight matrix, and the new data matrix is ​​obtained by screening from the original data matrix, thereby completing the data processing.

[0051] In a specific embodiment, the multispectral image is preprocessed to obtain a multispectral feature map, and the specific steps are as follows:

[0052] The multiple channel images of the multispectral image are processed using N1 different vegetation indices that can reflect the growth status of crops to obtain N1 vegetation index images (VIs);

[0053] Taking the sugarcane remote sensing dataset as an example, there are a total of N1=10 vegetation indices, namely, normalized difference vegetation index (NDVI), normalized difference red edge index (NDRE), green normalized difference vegetation index (GNDVI), enhanced vegetation index (EVI), modified anthocyanin content index (MACI), optimized soil conditioned vegetation index (OSAVI), simplified canopy chlorophyll content index (SCCCI), transformed chlorophyll absorption and reflectance index (TCARI), triangular green index (TGI), and visible atmospheric resistance index (VARI).

[0054] The calculation formula of the vegetation index is as follows:

[0055] NDVI = (NIR-R) / (NIR+R)

[0056] NDRE=(NIR-RE) / (NIR+RE)

[0057] GNDVI=(NIR-G) / (NIR+G)

[0058] EVI=2.5(NIR-R) / (NIR+6R-7.5B+1)

[0059] MACI=NIR / G

[0060] OSAVI=(1+0.16)(NIR-R) / (NIR+R+0.16)

[0061] SCCCI=NDRE / NDVI

[0062] TCARI=3[RE-R-0.2(RE / G)(RE / R)] / OSAVI

[0063] TGI=-0.5[(668-475)(RG)-(668-560)(RB)]

[0064] VARI=(G*R) / (G+RB)

[0065] In the formula, RE represents the red edge channel value of the sensor, NIR represents the near infrared channel value of the sensor, R represents the standard red channel value of the sensor, G represents the standard green channel value of the sensor, and B represents the standard blue channel value of the sensor;

[0066] The 10 vegetation index maps were gridded synchronously based on the size of the biomass sampling area at the final yield harvest. For example, if the multispectral image resolution is 0.5 and the designated area is 2m*2m, each grid unit contains 16 pixels.

[0067] Then, N2=7 different derivatives are extracted in each grid unit of each vegetation index map to form N1×N2 multispectral feature maps, where N1 and N2 are both positive integers. For example, the maximum value of NDVI in each unit is taken to form the max NDVI feature map. There are 7 derivatives, namely, the maximum value max, the minimum value min, the average value avg, the standard value std, the lower quartile value p25, the median value p50, and the upper quartile value p75. Using derivatives as multispectral data features, each VIs map generates 7 multispectral feature maps, and a total of 70 multispectral feature maps are obtained.

[0068] In a specific embodiment, the laser radar image is preprocessed to obtain a laser radar feature map, and the specific steps are as follows:

[0069] Taking the size of the biomass sampling area at the time of final yield harvest as the standard, LasTools and other software were used to extract the required vegetation data in each regional unit to generate multiple lidar feature maps.

[0070] The vegetation data include maximum height max, average height avg, average square height qav, height standard deviation std, height skewness ske, height kurtosis kur, 5th to 95th height percentiles p05~p95, 5th to 95th scores from the ground to the height percentiles b05~b95, and the proportion of four height points d00~d03. A total of 48 prediction variables generate 48 lidar feature maps.

[0071] In a specific embodiment, the multispectral feature map and the laser radar feature map are respectively arranged into row vectors X in the same position order. 1,n , n represents the number of samples in the feature map;

[0072] Then, the row vectors of the multispectral feature map and the laser radar feature map, a total of m feature map transformations, are merged to obtain the original data matrix X that combines the multispectral feature map and the laser radar feature map m,n .

[0073] like Figure 2 As shown in the figure, taking the sugarcane prediction experiment as an example, the above extraction has 70 multispectral feature maps and 48 lidar feature maps, that is, a total of 118 feature maps, and there are 56 random biomass sampling points (key blocks with 2m*2m units and recording final yield data), so putting the 118 features of these 56 blocks together and converting them into row vectors will get 118 row vectors in the format of [1, 56], and finally merge them into an original data matrix X[118, 56], where the number of rows 118 represents the number of features and the number of columns 56 represents the number of samples.

[0074] In a specific embodiment, the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix, as follows:

[0075] According to the original data matrix X m,n =[118,56] dimension size and the number of new features set to initialize the bivariate parameter matrix Q k,m =[15,118] and P m,k =[118, 15], where m is the number of feature maps, n is the number of samples in the feature map, and k=5 is the number of new features;

[0076] In order to make the original data matrix X m,n ≈P m,k Q k,m X m,n(or X[118,56]≈P[118,15]*Q[15,118]*X[118,56]), which will solve for Q k,m and P m,k The problem is transformed into a nonlinear optimization problem:

[0077]

[0078] The alternating least squares method is used to randomly initialize the parameter p, use p and formula (1) to calculate the parameter q, and then use the parameter q and formula (1) to solve the parameter p. This process is repeated until it converges to the local minimum.

[0079] By continuously optimizing P m,k and Q k,m , get the best new feature matrix Q k,m X m,n , and the new feature weight matrix Q k,m In this embodiment, Q k,m X m,n The dimension of this matrix is ​​[15, 56]), Q k,m is Q[15,118].

[0080] In a specific embodiment, the new data matrix is ​​obtained by screening from the original data matrix as follows:

[0081] Initialize an empty new data matrix Y k,n ; This embodiment Y k,n is [15, 56].

[0082] For the new feature weight matrix Q k,m Traverse each row of , find the data with the largest value, and record the column coordinate of the data as index; in this embodiment, Q k,m It is [15, 118], with a total of 15 rows.

[0083] In the original data matrix X m,n Retrieve the data of the index row from the matrix and merge it into the new data matrix Y k,n middle;

[0084] After k = 15 traversals, the new data matrix Y after key feature extraction is obtained. k,n In this embodiment, the new data matrix Y k,n In this embodiment, 15 features are selected, and a new data matrix of 56 sample data is obtained, as shown in FIG. Figure 3 shown.

[0085] In this embodiment, all sugarcane planting areas are first divided into 2m*2m grids, and then drones are used to carry sensors to collect multispectral images and lidar data that can accurately locate the grid positions in the early stages of sugarcane growth. The yield values ​​of each yield sampling area are recorded at the end of the harvest (that is, the number of kg of sugarcane in the black area is recorded).

[0086] The multispectral image and lidar data are stored in the form of pixels (with a resolution less than 2m*2m).

[0087] Example 2

[0088] like Figure 4 As shown, a method for predicting crop yield from multi-source remote sensing in the early season includes the following steps:

[0089] Preprocess the multispectral image and lidar data to obtain multiple multispectral feature maps and lidar feature maps corresponding to the target area;

[0090] All multispectral feature maps and lidar feature maps are combined into an original data matrix, and the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix;

[0091] Extract the original feature number that accounts for the largest proportion of each new feature from the new feature weight matrix, and filter the new data matrix Y from the original data matrix k,n , thus completing the data processing.

[0092] By collecting the yield data Z of the biomass sampling area at the corresponding final yield harvest 1，n , the new data matrix Y k,n With production data Z 1，n Input deep neural network for regression training, construct the potential relationship function between remote sensing data features and final yield, obtain the yield prediction model, and retain the input feature parameters of the yield prediction model;

[0093] The final yield prediction model is used to predict crop yield.

[0094] The new data matrix Y[15, 56] and the yield data Z[1, 56] at the same grid position in step 1 are sent to the deep neural network for nonlinear regression fitting training to construct the potential relationship function between the remote sensing data characteristics and the final yield. This can obtain a yield prediction model with higher prediction accuracy while retaining the actual significance of the model input feature parameters and improving the interpretability of the yield prediction method.

[0095] In this embodiment, the data is preprocessed into a feature map in grid units rather than the original feature map in pixel units in order to keep it on the same scale as the production data. The feature values ​​in grid units can correspond one-to-one to the production values ​​at the same grid position, thereby constructing the relationship between input and output.

[0096] In this embodiment, the data used is farmland remote sensing data, such as Figure 5 As shown in FIG. 1 , it is the yield prediction model performance obtained for sugarcane by the present invention. Figure 6 As shown, it is a schematic diagram of the fitting effect comparison of the yield prediction model performance obtained by the traditional method. Therefore, according to the comparison, it can be seen that the yield prediction model obtained by the present invention has a better effect.

[0097] Example 3

[0098] A computer device comprises a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the following steps are implemented:

[0099] Preprocess the multispectral image and lidar data to obtain multiple multispectral feature maps and lidar feature maps corresponding to the target area;

[0100] All multispectral feature maps and lidar feature maps are combined into an original data matrix, and the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix;

[0101] The original feature number that accounts for the largest proportion of each new feature is extracted from the new feature weight matrix, and the new data matrix is ​​obtained by screening from the original data matrix, thereby completing the data processing.

[0102] Obviously, the above embodiments of the present invention are only examples for clearly explaining the present invention, and are not intended to limit the implementation methods of the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention shall be included in the protection scope of the claims of the present invention.

Claims

1. A method for processing multi-source remote sensing yield data in the early season of crops, characterized in that: The following steps are involved: Preprocess the multispectral image and lidar data to obtain multiple multispectral feature maps and lidar feature maps corresponding to the target area; All multispectral feature maps and lidar feature maps are combined into an original data matrix, and the original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix; The original data matrix is ​​subjected to non-negative matrix decomposition of bivariate parameters to obtain a new feature weight matrix, as follows: According to the original data matrix X m,n The dimension size and the number of new features set are initialized to construct a bivariate parameter matrix Q k,m and P m,k , where m is the number of feature maps, n is the number of samples in the feature map, and k is the number of new features; In order to make the original data matrix X m,n ≈P m,k Q k,m X m,n , will solve Q k,m and P m,k The problem is transformed into a nonlinear optimization problem: The alternating least squares method is used to randomly initialize the parameter p, and then use p and formula (1) to calculate the parameter q. Then, use the parameter q and formula (1) to solve the parameter p. This process is repeated until it converges to the local minimum. By continuously optimizing P m,k and Q k,m , get the best new feature matrix Q k,m X m,n , and the new feature weight matrix Q k,m ; The original feature number that accounts for the largest proportion of each new feature is extracted from the new feature weight matrix, and the new data matrix is ​​obtained by screening from the original data matrix, thereby completing the data processing.

2. The method for processing crop season early multi-source remote sensing yield data according to claim 1, characterized in that: Preprocess the multispectral image to obtain the multispectral feature map. The specific steps are as follows: The multiple channel images of the multispectral image are processed using N1 different vegetation indices that can reflect the growth status of crops to obtain N1 vegetation index images; Taking the size of the biomass sampling area at the final yield harvest as the standard, N1 vegetation index maps were gridded synchronously; Then, N2 different derivatives are extracted in each grid unit of each vegetation index map to form N1×N2 multispectral feature maps, where N1 and N2 are both positive integers.

3. The method for processing crop season early multi-source remote sensing yield data according to claim 2, characterized in that: The vegetation index includes N1=10 items in total, namely normalized vegetation index, normalized difference red edge index, green light normalized difference vegetation index, enhanced vegetation index, modified anthocyanin content index, optimized soil conditioning vegetation index, simplified canopy chlorophyll content index, transformed chlorophyll absorption and reflectance index, triangular green index, and visible atmospheric resistance index.

4. The method for processing crop season early multi-source remote sensing yield data according to claim 2, characterized in that: The derivatives include N2=7 items in total, including maximum value, minimum value, average value, standard value, lower quartile value, median value, and upper quartile value.

5. The method for processing crop season early multi-source remote sensing yield data according to claim 1, characterized in that: The laser radar image is preprocessed to obtain the laser radar feature map. The specific steps are as follows: Taking the size of the biomass sampling area at the time of final yield harvest as the standard, the required vegetation data is extracted in each regional unit to generate multiple lidar feature maps.

6. The method for processing crop season early multi-source remote sensing yield data according to claim 5, characterized in that: The vegetation data include maximum height, average height, mean square height, height standard deviation, height skewness, height kurtosis, 5th to 95th height percentiles, scores from the 5th to 95th from the ground to the height percentiles, and the proportion of the four height points.

7. The method for processing crop season early multi-source remote sensing yield data according to claim 1, characterized in that: Set the multispectral feature map and the lidar feature map into row vectors X in the same position order. 1,n , n represents the number of samples in the feature map; Then, the row vectors of the multispectral feature map and the laser radar feature map, a total of m feature map transformations, are merged to obtain the original data matrix X that combines the multispectral feature map and the laser radar feature map m,n .

8. The method for processing multi-source remote sensing yield data of crops in early season according to claim 1, characterized in that: The new data matrix is ​​obtained by screening from the original data matrix, as follows: Initialize an empty new data matrix Y k,n ; For the new feature weight matrix Q k,m Traverse each row to find the data with the largest value, and record the column coordinate of the data as index; In the original data matrix X m,n Retrieve the data of the index row from the matrix and merge it into the new data matrix Y k,n middle; After k traversals, the new data matrix Y after key feature extraction is obtained k,n .

9. A method for predicting crop yield data from multi-source remote sensing in the early season, characterized by: The following steps are involved: The new data matrix Y is obtained by the method for processing crop season early multi-source remote sensing yield data as described in any one of claims 1 to 8. k,n ; By collecting the yield data Z of the biomass sampling area at the corresponding final yield harvest 1，n , the new data matrix Y k,n With production data Z 1，n Input deep neural network for regression training, construct the potential relationship function between remote sensing data features and final yield, obtain the yield prediction model, and retain the input feature parameters of the yield prediction model; The final yield prediction model is used to predict crop yield.

Images