Crop phenotype parameter automatic calculation and extraction method based on multi-source remote sensing image

Abstract

The invention relates to the technical field of crop monitoring, in particular to a crop phenotypic parameter automatic calculation and extraction method based on a multi-source remote sensing image. Comprising the steps that all images are determined to be in the same coordinate system through geographical registration, point cloud file conversion and ground specific point matching, and image registration is achieved; performing automatic or semi-automatic segmentation on the images determined in the same coordinate system to obtain a multi-source crop remote sensing image cell segmentation map; extracting plant phenotypic parameters, physical parameters and chemical parameters from the multisource crop remote sensing image cell segmentation map by using an algorithm; the plant phenotype parameters comprise a vegetation index, a plant height, a surface area, a volume, a canopy coverage degree and a vegetation projection area; the physical parameters comprise a canopy average temperature value, a canopy temperature standard deviation and a canopy temperature variation coefficient; the chemical parameters comprise chemical elements such as nitrogen, phosphorus, potassium, calcium and magnesium in soil and vegetation. The method has the advantages that large-scale data processing and high-precision area prediction are realized, and the crop growth monitoring capability is enhanced.

Description

Technical Field

[0001] The present invention relates to the technical field of crop monitoring, and in particular to a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images. Background Art

[0002] In modern agricultural production and research, accurate acquisition of crop phenotypic parameters is of great significance for crop growth monitoring, yield prediction, pest and disease control, and variety improvement. However, traditional field crop phenotypic detection methods mainly rely on manual operations. Workers need to go deep into the fields and measure many parameters of crops one by one, such as plant height, leaf area, and biomass. This method has many disadvantages: on the one hand, the accuracy of manual detection is difficult to guarantee, and differences in operating habits and experience of different testers will cause fluctuations in measurement results; on the other hand, manual detection is inefficient, consumes a lot of manpower and time, and cannot meet the needs of rapid monitoring of large areas of farmland. In addition, many parameters need to be sampled in the field and brought back to the laboratory for measurement. This process not only increases the complexity of the operation, but may also cause changes in the samples during transportation and processing, introducing additional errors.

[0003] In recent years, with the continuous development of remote sensing technology, its application in the extraction of crop phenotypic parameters has gradually attracted attention. Crop information is primarily obtained through multispectral, hyperspectral, thermal infrared, and lidar technologies. Multispectral remote sensing images can capture crop reflectance information in different wavelengths, providing a basis for estimating chlorophyll content, nitrogen nutrition status, and other indicators. Hyperspectral remote sensing images, with their high spectral resolution, enable more detailed analysis of the chemical composition and physiological status of crop leaves. Thermal infrared remote sensing images can be used to monitor crop transpiration and water stress. LiDAR technology can obtain three-dimensional structural information of crops, such as plant height and canopy structure. The fusion of multi-source remote sensing images is expected to overcome the shortcomings of traditional manual detection methods and achieve efficient and accurate extraction of crop phenotypic parameters. However, the process of extracting crop phenotypic parameters through multi-source remote sensing image fusion still faces some urgent challenges, such as the differences in the spatiotemporal resolution of data from different sensors and the need to improve the accuracy and efficiency of data fusion algorithms. Therefore, it is urgent to develop an efficient and accurate automatic extraction method of crop phenotypic parameters based on multi-source remote sensing images to promote the process of agricultural modernization and improve the level of agricultural production management. Summary of the Invention

[0004] In order to solve the above problems, the present invention provides a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images.

[0005] The present invention aims to provide a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images, which specifically comprises the following steps: S1. Through georegistration, point cloud file conversion, and matching of specific ground points, all images are determined in the same coordinate system to achieve image registration; S2. Automatically or semi-automatically segment the images determined in the same coordinate system to obtain the sub-region segmentation maps of multi-source crop remote sensing images; S3. Extract plant phenotypic parameters, physical parameters, and chemical parameters from the sub-region segmentation maps of multi-source crop remote sensing images; The plant phenotypic parameters include vegetation indices, plant height, surface area, volume, canopy coverage, and vegetation projected area; the physical parameter is the canopy temperature feature, and the canopy temperature feature includes the canopy average temperature value, canopy temperature standard deviation, and canopy temperature coefficient of variation; the chemical parameters include soil chemical elements and vegetation chemical elements. The soil chemical elements include nitrogen, phosphorus, potassium, calcium, magnesium, and organic matter, and the vegetation chemical elements include nitrogen, phosphorus, potassium, calcium, and magnesium.

[0006] Preferably, in step S3, the vegetation indices are calculated through the band values of the images, and the vegetation indices include the normalized difference vegetation index, enhanced vegetation index, and optimized soil-adjusted vegetation index; The random sample consensus algorithm is used to calculate the plant height; The Delaunay triangulation or Poisson reconstruction algorithm is used to calculate the surface area; Using Gauss's theorem, through the conversion of volume integral and surface integral, the volume is calculated; Project the three-dimensional vegetation point cloud onto the horizontal plane for rasterization, count the ratio of the number of covered rasters to the total number of effective rasters, and then perform interpolation filling or proportionally enlarge the number of covered rasters to obtain the canopy coverage; By counting the number of pixels belonging to the vegetation in the segmented vegetation image and then multiplying by the true area corresponding to each pixel, the vegetation projected area is obtained.

[0007] Preferably, in step S3, the K-means clustering algorithm is used to further segment the sub-region segmentation maps of multi-source crop remote sensing images, and the canopy average temperature value, canopy temperature standard deviation, and canopy temperature coefficient of variation are extracted from the segmented canopy region; The formulas for the canopy temperature standard deviation and the canopy temperature coefficient of variation are as follows: ; ; In the formula: represents the canopy average temperature, represents the temperature of soybean pixels in the UAV thermal infrared image, is the number of pixels; The formula for the K-means clustering algorithm is: ; In the formula: is the within-cluster sum of squares; is the number of clusters; is the set of data points in the i-th cluster; is the cluster is a data point in μ i is the center point of the cluster ; is the data point and its cluster center μ i is the square of the Euclidean distance between them.

[0008] Preferably, the extraction of chemical parameters in step S3 includes data cleaning, data loading and preprocessing, model training, model tuning and training evaluation, and prediction; specifically including: S331. Data cleaning: Use the validate_data function for data loading, header row detection, interactive verification, data cleaning, error reporting, and saving the corrected data; Data cleaning specifically includes: First, perform true content cleaning. Clean the last element of each row of data. First, remove spaces, non-numeric and non-decimal point characters, and only retain numbers and decimal points. If the true content after cleaning is empty, throw an exception. Then, perform spectral data cleaning. Traverse the other values of each row of spectral data, replace commas with dots, and remove illegal characters. If a certain band data is invalid, throw an exception; Error reporting specifically includes: Perform error handling. If an exception occurs during the cleaning process, record the error information of that row in the errors list; S332. Data loading and preprocessing: Use optimizer.preprocess for preprocessing, including outlier handling, data augmentation, and data normalization; S333. Model training: Use the ridge regression model for inverting the soil chemical element content; use the random forest and XGBoost models for inverting the vegetation chemical element content, and select the best-performing model for prediction; S334. Model tuning and training evaluation: Use GridSearchCV to perform grid search on the hyperparameters of each model; select the best parameters through cross-validation, and the scoring criterion is selected as R 2 : ; In the formula: y true is the true value, y pred is the predicted value, is the mean of the true values; S335. Prediction: Use the trained model to predict new hyperspectral images.

[0009] Preferably, step S332 specifically includes: Outlier handling: Identify and remove outliers using the median and median absolute deviation. The calculation formula is: ; Then remove outliers that exceed 3 times MAD ; Data augmentation: Expand the dataset by adding noise to generate augmented data. The noise is from a normal distribution, and the noise standard deviation is 0.05. Use StandardScaler to standardize the data so that the mean is 0 and the standard deviation is 1.

[0010] Preferably, the objective of the ridge regression model in step S333 is to minimize the following loss function: ; Where: is the least squares loss, that is, the error of the regression model fitting. is the L2 regularization term. λ is the regularization parameter. β is the model parameter vector. is for the i th sample's feature vector. β j represents the β th j element of the parameter vector is the number of features. The loss function formula of the XGBoost model is as follows: ; Where: is the true label of the sample. is the prediction of the model. is the prediction function of the th tree model. is the regularization term of the model complexity. Use the methods of cross-validation and test set evaluation to compare the performance of the random forest model and the XGBoost model, and select the model with better performance for the final prediction.

[0011] Preferably, step S334 also includes plotting learning curves, feature importance, residual analysis, and scatter plots of true values versus predicted values; visually analyze the performance of the model; specifically as follows: Learning curve: Use learning_curve to plot the performance of the model under different training set sizes; R of the training set and validation set 2 Changes with the number of training samples; the learning curve can help determine whether the model is overfitting or underfitting; Feature Importance: Use model.feature_importances to get the importance score of each feature; feature importance reflects the contribution of each feature to model prediction; Residual analysis: calculate the difference between the true value and the predicted value and plot the residual distribution; Scatter plot of true values versus predicted values: Use a scatter plot to compare the true values with the predicted values; ideally, the points are concentrated near the diagonal line; the diagonal line represents y=x, that is, the predicted value is equal to the true value.

[0012] Preferably, step S335 specifically includes: Load model and standardization objects: load the trained model and standardizer; load the chemical element content range saved during training; Process image data: read the input image in blocks and extract the spectral characteristics of each pixel; normalize the non-zero valid pixels; use the model to predict the normalized chemical element content value; and restore it to the true value through inverse normalization; Limit the prediction results: Use np.clip to limit the prediction results to the range of the actual chemical element content in the training set; the all-zero area is forced to be assigned to 0; Generate prediction map: Write the predicted value of each image block into a new image file, converting the hyperspectral data into a chemical element content map.

[0013] Preferably, step S2 specifically includes: S21. semi-automatically segmenting unplanted bare soil areas and areas with unknown crop growth; S22. For crop areas with known growth conditions, use the improved Grounded-SAM segmentation model or traditional algorithms to perform batch, fully automatic segmentation processing on multi-source sensor images and output the segmentation results. The traditional algorithms include NDVI algorithm, improved EXGR algorithm or Otsu algorithm. The improved EXGR algorithm sets a threshold according to the EXGR value, and defines the pixel area greater than the threshold as vegetation, and the rest as non-vegetation area. The segmentation formula is as follows: EXGR=3G-2.1RB; Among them, G, R, and B represent the intensity values of the green, red, and blue channels of each pixel in the image respectively; S23. Optimize the segmentation results and finally output an optimized cell segmentation map of the multi-source crop remote sensing image.

[0014] Preferably, in step S22, the specific sub-steps of using the improved Grounded-SAM segmentation large model are as follows: S221. Determine the empirical model parameters: Set the empirical model according to the known crop planting parameters and growth characteristics; S222. Use the improved Grounded-SAM segmentation large model and input the prompt words; Receive the image and text prompts through DINO-X, output the bounding boxes in the image that are semantically related to the text, and determine the positions of the plants; S223. Use the bounding boxes output by DINO-X as prompts to further refine each bounding box area and generate accurate segmentation masks; S224. Crop the multi-source sensor images according to the boundary range set by the empirical model; Apply the segmentation masks generated in step S223 to the cropped images for plot image segmentation, and extract the independent images of each plot; S225. Overlay and display the segmentation results with the original image to generate a new image, or output the segmentation mask image.

[0015] Compared with the prior art, the present invention can achieve the following beneficial effects: 1. Improve the efficiency of crop phenotype detection: The present invention adopts an automated image segmentation and phenotype parameter calculation method, reduces the dependence on manual measurement, and greatly improves the efficiency of field crop phenotype detection.

[0016] 2. Applicable to multi-source remote sensing data: This method can be simultaneously applicable to various remote sensing data such as RGB, multi-spectral, hyperspectral, thermal infrared, and lidar, and can extract rich crop growth parameters.

[0017] 3. Improve data accuracy and consistency: Through automated segmentation and calculation, artificial errors are avoided, ensuring data consistency and improving measurement accuracy.

[0018] 4. Support large-scale data processing: Combining machine learning and deep learning models (Random Forest, XGBoost), it can achieve efficient processing of a large number of data, providing support for agricultural big data analysis.

[0019] 5. Enhance the ability to monitor crop growth: By extracting various vegetation indices, canopy structures, physical parameters, and chemical parameters, it can accurately monitor the crop growth status and provide a scientific basis for precision agriculture.

[0020] 6. Improve the generalization ability of the model: Adopting techniques such as cross-validation, regularization, and hyperparameter optimization to improve the generalization ability of the model, making it applicable to different regions and different crop types.

[0021] 7. Visual analysis supports decision-making: By means of visualization methods such as learning curves, feature importance analysis, and residual analysis, the interpretability of data analysis is improved to assist agricultural research and decision-making.

[0022] 8. Achieve high-precision regional prediction: Based on the trained machine learning model, the hyperspectral image is processed in blocks and predicted pixel by pixel. Combining standardization and inverse normalization operations, the true values of chemical element contents such as calcium are effectively restored, and the abnormal output is controlled by predicting value clipping, improving the accuracy and reliability of regional inversion, and providing data support for crop nutrient monitoring and soil improvement. Description of the Drawings

[0023] Figure 1 It is the visualization display effect diagram of the registered two-flight files provided by the embodiment of the present invention.

[0024] Figure 2 It is the effect diagram of the nitrogen content of soil chemical elements provided by the embodiment of the present invention.

[0025] Figure 3 It is the effect diagram of the phosphorus content of soil chemical elements provided by the embodiment of the present invention.

[0026] Figure 4 It is the effect diagram of the potassium content of soil chemical elements provided by the embodiment of the present invention.

[0027] Figure 5 It is the effect diagram of the calcium content of soil chemical elements provided by the embodiment of the present invention.

[0028] Figure 6 It is the effect diagram of the magnesium content of soil chemical elements provided by the embodiment of the present invention.

[0029] Figure 7 It is the effect diagram of the soil organic matter content provided by the embodiment of the present invention.

[0030] Figure 8 It is the analysis results of the random forest and XGBoost models for the nitrogen content of vegetation chemical elements provided by the embodiment of the present invention.

[0031] Figure 9 It is the analysis results of the random forest and XGBoost models for the phosphorus content of vegetation chemical elements provided by the embodiment of the present invention.

[0032] Figure 10 It is the analysis results of the random forest and XGBoost models for the potassium content of vegetation chemical elements provided by the embodiment of the present invention.

[0033] Figure 11It is the analysis results of the random forest and XGBoost models for the calcium content of the chemical elements in vegetation provided by the embodiments of the present invention.

[0034] Figure 12 It is the analysis results of the random forest and XGBoost models for the magnesium content of the chemical elements in vegetation provided by the embodiments of the present invention.

[0035] Figure 13 It is the predicted content distribution maps of nitrogen and phosphorus provided by the embodiments of the present invention.

[0036] Figure 14 It is the predicted content distribution maps of potassium, calcium, and magnesium provided by the embodiments of the present invention. Detailed implementation manners

[0037] In the following, embodiments of the present invention will be described with reference to the accompanying drawings. In the following description, the same modules are denoted by the same reference numerals. In the case of the same reference numerals, their names and functions are also the same. Therefore, their detailed descriptions will not be repeated. In order to make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0038] Embodiment 1 This embodiment provides a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images, which specifically includes the following steps: S1. Through georegistration, point cloud file conversion, and ground specific point matching, all images are determined in the same coordinate system to achieve image registration; specifically including the following steps: S11. Extract the flight data of the unmanned aerial vehicle (UAV). The flight data includes GNSS observation files (.obs), navigation data files (NAV files), event record files (EVENT files), binary log files (BIN files), and marker files (MRK files); specifically as follows: GNSS observation files (.obs): Standard data files that mainly store the original measurement data from the global navigation satellite system; follow a specific structured format to facilitate interoperability between different software and hardware platforms; the file header contains information about the site, receiver, antenna, and other metadata; Navigation data files (NAV files): Store information about the flight path, GPS coordinates, altitude, and other relevant data, which can be used for georeferencing of images, flight analysis, or troubleshooting; Event record files (EVENT files): Record various events that occur during the flight, such as when the UAV takes off, when an image is captured, or when a certain waypoint is reached; it can help understand the progress of the flight and identify any problems that may occur during the mission; Binary log file (BIN file): It contains detailed information about the drone flight, including telemetry data, sensor readings, and system information, and can be used for in-depth analysis and troubleshooting of the drone's performance; Marker file (MRK file): Records relevant information about the photo-taking points, where: Column 1: Serial number of the photo-taking point, i.e., the serial number of the photo record information stored in this folder; Column 2: Seconds within GPS week, indicating the photo-taking moment, expressed in the form of seconds within GPS week (TOW); Column 3: GPS week, indicating the GPS week of the photo-taking moment; Column 4: Compensation value in the north direction, in millimeters, with north being positive; Column 5: Compensation value in the east direction, in millimeters, with east being positive; Column 6: Compensation value in the elevation direction, in millimeters, with down being positive; Column 7: Longitude after compensation; Column 8: Latitude after compensation; Column 9: Ellipsoidal height; Column 10: Positioning standard deviation (north); Column 11: Positioning standard deviation (east); Column 12: Positioning standard deviation (elevation); Column 13: Positioning status.

[0039] S12. Assign the geographical location information collected in the flight data extracted in step S11 to the image to achieve the georegistration of the image; Before synchronizing the image time with the flight data, it is necessary to extract the metadata information in the image file, especially the DateTimeOriginal field in the EXIF information. This field records the original shooting time of the image and is an important basis for time synchronization and calibration.

[0040] To achieve the registration process, it is necessary to use the flight data file collected during the drone flight. The file contains key data such as geographical location information, flight attitude information, and timestamps. Through these data, the spatial trajectory of the projection center can be solved, and then the georegistration of the image can be achieved; to solve the spatial trajectory of the projection center, the core formula is based on the collinearity equation and spatio-temporal interpolation, and the formula is as follows: ; In the formula: λ represents the proportionality factor between the distance from a point on the image to the image center and the actual distance of the point on the ground. It is a dimensionless value used to maintain the consistency of the proportional relationship in the collinearity equation; (X, Y, Z) represents the object space coordinates of the ground point; (X0, Y0, Z0) represents the projection center coordinates (GPS position to be solved); R represents the external orientation rotation matrix composed of the drone attitude angles (ω, φ, κ); f represents the camera focal length; (x, y) represents the image point coordinates; Specifically, it includes the following sub-steps: S121. Time synchronization and calibration: Ensure that the image capture time is consistent with the time in the GPS track log in the same time system (UTC or local time zone); if there is a system clock deviation, the image time needs to be compensated: t corrected =t image +Δt; Where: t corrected is the calibrated image capture time; t image is the original image capture time (extracted from the metadata of the image file, i.e., the DateTimeOriginal field in the EXIF information); Δt is the system clock deviation compensation amount, which can be obtained by photographing a reference point at a known location at the same time and calculating the time offset inversely.

[0041] S122. Trajectory interpolation algorithm: Use the calibrated image capture time t corrected as the interpolation time point, and adopt cubic spline interpolation to output the exact position (X, Y, Z), i.e., WGS84 coordinates, at each image capture moment; for the trajectory segment index j, the interpolation formula is: ; Where: t represents the time point for interpolation; j represents the trajectory segment index; t j represents the starting time point of the trajectory segment index j ; S j ( t ) represents the position function value of the trajectory segment index j at time t ; represents the second derivative of the position function of the trajectory segment index j at time t ; represents the second derivative of the position function of the trajectory segment index j +1 at time t j ; the coefficients a j , b j , c j , d j are solved through boundary conditions; The simplified linear interpolation formula is as follows: ; In the formula: X ( t ) represents the X coordinate value at time t; X {k} represents at timet {k} The X - coordinate value at X {k+1} represents the X - coordinate value at time t {k+1} ; t represents the time point at which the X - coordinate value needs to be estimated; t {k} and t {k+1} represent two time points with known X - coordinate values; The Y and Z coordinates follow the same formula as the X - coordinate: ; ; Principle summary: The GPS trajectory of the UAV is usually recorded at a frequency of 1Hz, while the shutter speed of the camera can reach multiple frames per second, which means that the sampling frequency of the GPS position information is lower than the shooting frequency of the images; in order to obtain the precise position (X, Y, Z) at each image shooting moment, interpolation of the GPS trajectory is required; the cubic spline interpolation in this step is applicable to high - precision aerial surveying and can ensure the continuity of position and speed.

[0042] S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S202 to the CGCS2000 coordinate system (meeting Chinese requirements); use the seven - parameter Bursa model for coordinate transformation: ; In the formula: [X CGCS , Y CGCS , Z CGCS represents the coordinates in the CGCS2000 coordinate system after transformation; [ΔX, ΔY, ΔZ] represents the translation vector between the WGS84 coordinate system and the CGCS2000 coordinate system; 1 + m represents the scale change factor, where m is the scale change parameter, representing the scale difference between the two coordinate systems; represents the rotation matrix, , , are the rotation parameters, representing the rotation angles around the X, Y, and Z axes respectively.

[0043] S124. Use the EGM2008 (Earth Gravitational Model 2008) model for elevation correction to convert the altitude measured by GPS (based on the ellipsoid) to orthometric height (based on the geoid); the elevation correction formula is: H ortho = H GPS - Ngeoid ; In the formula: H ortho represents the orthometric height; H GPS represents the altitude measured by GPS; N geoid represents the geoid height anomaly calculated by model interpolation.

[0044] S13. Generation of image data: Convert the LAS point cloud file with color information into a GeoTIFF image file with RGB bands; specifically including the following steps: S131. Read the LAS file and extract the point cloud data; check whether the LAS file contains color information (RGB three channels); if there is no color information, it is thrown as an exception to ensure the effectiveness of subsequent operations; Specifically, use the laspy.read method to read the input LAS file; the point cloud data includes spatial coordinates x, y and color channel data red, green, blue.

[0045] S132. Calculate the size of the grid: According to the spatial range (minimum and maximum coordinate values) of the point cloud data and the preset resolution (how many meters each pixel represents), calculate the number of columns (cols) and rows (rows) of the grid; the formula is as follows: ; ; In the formula: max x , min x , max y , min y are the maximum and minimum coordinate values of the point cloud data respectively; resolution is the resolution; The resolution of the grid determines the actual spatial range represented by each pixel. The resolution (resolution) is the side length of each pixel in space, usually in meters; according to the given resolution, the width and height of the grid can be calculated.

[0046] S133. Binned statistics: Use the binned_statistic_2d function to allocate the data of each RGB channel into the calculated grid; binned_statistic_2d is a SciPy function used to perform two-dimensional binned statistics; the formula is as follows: statistic = statistic method (data in grid cells); Among them, x and y are the spatial coordinates of the point cloud data, and data is the corresponding RGB color channel data; the bin edges are specified by the bins parameter, which is x here edges and y edges ; statistic specifies the calculation method for statistics, such as mean (average), max (maximum value), min (minimum value), etc.; The present invention uses the mean as the default statistical method, that is, the color values within each grid cell are averaged to generate an average grid.

[0047] S134. RGB merging: Merge the three color channels (red, green, blue) after binning into a three-dimensional RGB stack (rgb stack ), where each pixel contains the color values of the red, green, and blue channels, facilitating the generation of a color image; The specific operation is as follows: Initialize the RGB stack: Create a three-dimensional array (rgb stack ), whose shape is (rows, cols, 3), where rows and cols are the number of rows and columns of the grid respectively, and 3 represents the three RGB color channels. Fill the RGB stack: Fill the binned red, green, and blue channel data into the first three positions of the third dimension of rgb stack respectively. For example, the red channel data is filled into rgb stack [:,:,0], the green channel data is filled into rgb stack [:,:,1], and the blue channel data is filled into rgb stack [:,:,2]. Ensure that the data types of all channels are consistent before merging, usually unsigned 8-bit integers (uint8), representing the color value range from 0 to 255.

[0048] S135. Color scaling: Determine whether it is necessary to scale the color values from 16 bits to 8 bits according to actual needs; if scale color = True, then scale the color values from 16 bits to 8 bits by dividing the value of each color channel by 65535 (the maximum value of 16 bits) and multiplying by 255 (the maximum value of 8 bits); update the scaled color values to the RGB stack to ensure that the data types of all color channels are uint8; this step is selected according to the actual situation whether scaling is required; The linear scaling formula is as follows: ; In the formula: scaled color represents the scaled color value; color represents the original color value, from the color information in the LAS file.

[0049] S136. Generate a GeoTIFF file: Determine the metadata of the GeoTIFF file; Open the GeoTIFF file using the rasterio.open function and pass in the metadata; Write the three color channels of the RGB stack to the corresponding bands of the GeoTIFF file; The GeoTIFF file contains three bands (RGB), set appropriate data types for each band; And perform georeferencing using the given geotransform (Affine) and coordinate system (crs); The metadata of the GeoTIFF file includes band information, data type, geotransform, and coordinate system.

[0050] S14. Use specific ground points for image matching to achieve precise registration of two or more UAV remote sensing images; The specific steps are as follows: S141. Identify and select specific fixed ground points in the image, reassign fixed coordinate points to the selected specific ground points, and use the fixed coordinate points as the reference points for image registration; The specific ground points can be the fixed points in the field, such as houses, roads, water pipes, etc. existing in the image; Reassign fixed coordinate points to the specific ground points so that two or more images can be matched, keeping all the remote sensing images of the detection area in one coordinate system to achieve precise registration; S142. Designate one image as the reference image and the other as the image to be registered; The reference image must include standard map coordinates or RPC information and cannot be pixel coordinates, coordinates without projection information (arbitrary coordinate information), and pseudo coordinates (pseudo); There are no strict constraints on the image to be registered, but if there is no coordinate information, at least 3 homologous points need to be manually selected; Use the Image Registration Workflow tool in ENVI software to automatically, quickly, and accurately complete the registration between two images with geometric position deviations; Image Registration Workflow is an automatic, accurate, and fast image registration workflow that integrates the complicated parameter setting steps into a unified panel and can quickly and accurately achieve automatic registration between images with little or no manual intervention; The data includes a pair of data with the same resolution and different imaging times, and a pair of panchromatic and multispectral data with different resolutions and the same imaging time; The data formats supported by the registration tool are: ENVI, TIFF, NITF, JPEG2000, JPEG, Esri® raster layer, Geodatabase raster, Web Services; S143. Feature extraction and matching: Use the Cross Correlation method to automatically detect and match two images (I 1 and I the corresponding feature points in (2); the calculation formula is as follows: ; In the formula: CC ( I 1, I 2) represents the cross - correlation value between images I 1 and I 2; M and N respectively represent the number of rows and columns of the image; I 1( i , j ) and I 2( i , j ) respectively represent the pixel values of images I 1 and I 2 at the position ([[]] i , j ); cross - correlation measures the similarity between two images by calculating the product of the pixel values of the two images.

[0051] Cross - Correlation (CC) is a common method for measuring the similarity between two images. Its principle is to judge whether the images are aligned or how to align by calculating the similarity of pixel intensities between the images; cross - correlation is generally used for images with similar morphology or features (for example, two images taken by the same type of sensor or the same type of images, such as an optical image pair with an optical image); The algorithm logic is as follows: (1) Calculate the matching window: Select a fixed area (called a window) on image I 1 and scan the same - sized area on image I 2; (2) Slide the window to calculate the correlation: For each position, calculate the cross - correlation value within the window, that is, the weighted sum of the pixel values within the two image windows; (3) Maximum value position: Finally, select the position with the largest cross - correlation value, which is the estimated position of the best alignment of the two images.

[0052] S144. Geometric transformation calculation: Adopt a polynomial transformation model. By measuring the positions of image pixels and registering them with the corresponding geographical positions, use the paired points to fit the coefficients of the polynomial transformation; after fitting, use the obtained polynomial coefficients to perform geometric correction on the entire image; after correction, verify the accuracy of the geometric correction by comparing the corrected image with the reference image; Polynomial transformation usually uses a polynomial function to fit the position transformation of image pixels. The common form is a two - dimensional polynomial transformation, and the expression is as follows: ; ; In the formula: X and Y are respectively the row and column coordinates of the pixels in the original image; X’ and Y’ are respectively the row and column coordinates of the pixels in the corrected image, a i and b i (i = 0, 1, 2, 3, 4, 5, …, n) are polynomial coefficients; The geometric correction step maps the positions of the image pixels to their true geographical locations in the real world, achieving precise alignment of the image.

[0053] In this step, for the geometric correction of lidar point cloud data, first extract the orthographic projection of the lidar point cloud data, and then register the orthographic projection map with the original lidar point cloud image; the orthographic projection transforms the three-dimensional model into an orthographic projection map, which retains the geographical location information; since the two images have the same viewing angle in the orthographic direction and the surface features are consistent at this angle, therefore, the registration process mainly focuses on assigning precise geographical location information to the point cloud data at different heights to ensure the consistency of the point cloud data and the orthographic projection map in the geographical space.

[0054] S145. Resampling: Apply the Cubic Convolution method to resample the image to be registered to adapt to the spatial resolution of the reference image; The Cubic Convolution method is a high-order interpolation method that uses more pixel values around the target point to calculate the interpolation result, thus producing a smoother image; calculate the weighted value by selecting multiple pixels (usually a 4×4 area, containing 16 pixels) around the target pixel; each neighboring pixel is assigned a different weight value according to its relative position to the target pixel to achieve weighted neighboring pixels; The specific processing process is as follows: S1451. Determine the resampling area and collect neighboring pixel values: Select an area (a 4×4 area) around the target pixel, which contains multiple pixel points for calculating the interpolation value; collect the values of 16 pixel points within the selected area as the neighboring pixels of the target pixel for calculating the interpolation value of the target pixel; S1452. Calculate weights: Use the Cubic Convolution Kernel to calculate weights for each neighboring pixel; according to the calculated weights, weight the neighboring pixel values; finally, the interpolation of the target pixel (x’, y’) is the weighted average of the neighboring pixel values, expressed as: ; In the formula: I ( x ’, y ) represents the pixel value of the image at the position ([ x ’, y ); I ( x + i , y + j ) are the values of the neighboring pixels; i and j are the relative position offsets within the neighborhood; is the horizontal distance x’ from the neighborhood point x + i ; is the vertical distance y’ from the neighborhood point y + j ;

[0055] The convolution kernel calculates weights based on the pixel intervals from the target pixel position; the Catmull-Rom interpolation kernel is used, which is a smooth cubic function that can effectively reduce the jagged effect and blurring phenomenon during the interpolation process; pixels closer to the target pixel will receive higher weights, while pixels farther away will receive lower weights. S1453. Generate the resampled image: Repeat steps S1451 and S1452 for each pixel in the image to generate the entire resampled image; check the quality of the resampled image to ensure there are no obvious jagged effects or blurring phenomena; if necessary, adjust the parameters of the convolution kernel to optimize the results. S1454. Output the final image: Output the resampled image as the final result for subsequent image analysis or applications. The bicubic interpolation method uses more pixel value information, so it can usually produce higher-quality results when enlarging images; however, its computational complexity is also relatively high.

[0056] S146. Use the multiple dynamic link display function of ENVI to compare the reference image and the registered image and evaluate the registration accuracy.

[0057] For images acquired by different sensors (such as RGB, multispectral, hyperspectral, lidar), the above method can be used for registration. For example, an RGB image is three-dimensional, a multispectral image is five- or multi-dimensional, a hyperspectral image is a multi-dimensional image, and a lidar image is also a three-dimensional image. Sensors vary from one another. There are RGB or pseudo-RGB bands among RGB, hyperspectral, and multispectral images, but the pixel accuracy differs between different sensors. As mentioned above, there are also problems with images from the same sensor in different flight missions. Therefore, it is necessary to solve the problems of different image resolutions, different image geographical location information, and different image formats. For RGB, multispectral, and hyperspectral images, the methods mentioned above are used to assign geographical information to them. For the registration of geographical location information between RGB images and point clouds, it is necessary to first extract the orthophoto of the lidar point cloud data and then register the orthophoto with the original lidar point cloud image. The orthophoto transforms the three-dimensional model into an orthographic projection image that retains the geographical location information. Since the two images have the same viewing angle at the orthographic angle and the surface features are consistent at this angle, the registration process mainly focuses on assigning accurate geographical location information to point cloud data at different heights to ensure the consistency of the point cloud data and the orthophoto in the geographical space.

[0058] Figure 1 The visualization effect of two flight files is shown (with flight inflection point deviation and altitude deviation). The left figure is a three-dimensional flight trajectory diagram (Flight Trajectories (3D)), which shows the spatial trajectories of two flights (Flight 1 and Flight 2). The horizontal axis is the UTM X coordinate, the vertical axis is the UTM Y coordinate, and the vertical axis is the flight altitude (Altitude). The right figure is a two-dimensional plane trajectory diagram (Flight Trajectories (2D)), which removes the altitude dimension and more clearly shows the ground coverage paths of the two flights. It can be seen from the figure that for flights with the same route plan, there are altitude and position differences in actual execution.

[0059] For the registration of thermal infrared images on different dates, the thermal infrared images are used to register the hyperspectral images to achieve minimal error.

[0060] S2. Automatically or semi-automatically segment the images determined to be in the same coordinate system to obtain the sub-region segmentation map of the multi-source crop remote sensing images; Specifically, semi-automatically segment the unplanted bare soil areas and areas of crops with unknown growth conditions; automatically segment the areas of crops with known growth conditions; The method of semi-automatic segmentation is as follows: For the bare soil areas without crops, manually segment the plot boundaries; Since the maximum growth range of some remote sensing images of crops with unknown growth conditions cannot be directly determined, it is necessary to manually segment the crop plots with known maximum growth ranges. Specifically, select the UAV images of the crop plots that have grown to the maximum projected area for manual segmentation. Since the maximum projected area is known, batch cropping of other images is performed based on this, and then plot image segmentation is carried out; The method of fully automatic segmentation is as follows: Taking crops with known growth conditions as an example, for soybeans, 6 to 8 plants are planted in a plot of 2×1.5 meters. According to the empirical model, the length and width of its maximum growth area will not exceed 2 meters; Based on the empirical model, by matching the image coordinates with the geographical location information (GNSS coordinate system) (i.e., step S1), batch segmentation of multi-source sensor images can be achieved.

[0061] Specifically, it includes the following steps: S21. Perform semi-automatic segmentation on the bare soil areas without crops and the areas of crops with unknown growth conditions; Among them, for the bare soil areas without crops, manually segment the plot boundaries; For the areas of crops with unknown growth conditions, first select the UAV images of the crop plots that have grown to the maximum projected area for manual segmentation, and then, based on the plot boundary range obtained through manual segmentation, batch crop other images, and then perform plot image segmentation; S22. For the areas of crops with known growth conditions, use the improved Grounded-SAM segmentation large model to perform batch full-automatic segmentation processing on multi-source sensor images and output the segmentation results; Specifically, it includes the following sub-steps: S221. Determine the empirical model parameters: Set the empirical model according to the known crop planting parameters and growth characteristics; This step provides prior knowledge for subsequent automatic segmentation and reduces unnecessary calculations; For example, for soybeans, 6 to 8 plants are planted in a plot of 2×1.5 meters. According to the empirical model, the length and width of the maximum growth area will not exceed 2 meters.

[0062] S222. Use the improved Grounded-SAM segmentation large model and input the prompt word "plant"; Receive the image and text prompt through DINO-X, output the bounding boxes in the image that are semantically related to the text, and determine the positions of the plants; DINO-X is an open-set object detector that can detect relevant objects in images according to arbitrary free-form text prompts. Since this large model performs segmentation based on the prompts, only by inputting the prompt "plant", the model will detect and segment the crops in the picture plot based on the vocabulary. The cropped plot image only contains two parts: land and crops, so the segmentation is extremely simple. The purpose of this step is to guide the model to detect and segment the target crops through natural language prompts and quickly locate the target crop area in the image.

[0063] S223. Use the bounding boxes output by DINO-X as prompts to further refine each bounding box region to generate accurate segmentation masks.

[0064] S224. Crop the multi-source sensor images according to the boundary range set by the empirical model; apply the segmentation masks generated in step S223 to the cropped images for plot image segmentation to extract independent images of each plot.

[0065] S225. Overlay the segmentation result with the original image to generate a new image, or output the segmentation mask image, which includes the confidence or visualization contour of the segmentation region; In a specific embodiment, for the RGB image of the field plot, the Grounded-SAM model was used for segmentation to extract the regions of soybean plants; the segmentation mask (for example, green represents vegetation and red represents soil) was overlaid on the original RGB image to generate a new image, where the vegetation region of the original image was covered with green and the soil region was covered with red.

[0066] S23. Optimize the segmentation result, including constructing a height matrix to determine the shape of the vegetation region, using the monotonic stack algorithm to calculate the maximum rectangle area row by row, and calculating the maximum growth projection area of each plot based on the empirical model, and finally output the plot segmentation map of the optimized multi-source crop remote sensing image; specifically, it includes the following sub-steps: S231. Construct a height matrix: For the connected components (foreground is 1, background is 0) in the independent image (binary image) segmented and extracted in step S224, traverse each pixel row by row, calculate the number of consecutive 1s upward at each position to form a height matrix height, and determine the shape and size of the connected components (vegetation regions) of each plot; the height matrix formula is as follows: ; In the formula, i represents the row index in the image, and j represents the column index in the image; S232. Calculate the maximum rectangle area row by row: For each row of the height array, use the monotonic stack algorithm to solve the maximum rectangle area of the current row; the specific process is as follows: Initialize the stack: The stack is used to save the indices of the columns with increasing heights. The initial state of the stack is empty, and the heights corresponding to the bottom to the top of the stack increase in sequence. Traverse each element in the current row: When the stack is not empty and the current height is less than the height at the top of the stack, pop the top element top of the stack, and calculate the area of the rectangle with the height height[top]. ; Among them, h is the popped height, i is the current index, and s top-1 is the index of the new top of the stack after popping; Push the current index onto the stack. Process the remaining elements in the stack: After the traversal, perform a similar operation on the remaining elements in the stack to calculate the possible areas; Update the global maximum area: After processing each row, record the maximum rectangle area of that row; Traverse all rows and take the maximum value among all rows as the global maximum area (representing the maximum rectangle area of the vegetation area in all rows).

[0067] The simplified formula is expressed as follows: Height matrix: Accumulate the heights of consecutive 1s in each column; Maximum area of a single row: Use a monotonic stack to find the maximum width, and the area = height × width; Global maximum area: Take the maximum value among all rows; S233. Calculate the maximum circumscribed rectangle based on an empirical model: For each connected region of the independent small cell image segmented, calculate its maximum circumscribed rectangle; Through this step, the maximum growth projection area of each small cell can be determined, providing a basis for subsequent analysis.

[0068] Compare the results of automatic segmentation and semi-automatic segmentation to verify the accuracy and consistency of the segmentation.

[0069] S3. Extract plant phenotypic parameters, physical parameters, and chemical parameters from the small cell segmentation map of the multi-source crop remote sensing image; Specifically, it includes the following sub-steps: S31. Extract plant phenotypic parameters; Plant phenotypic parameters include vegetation indices, plant height, surface area, volume, canopy coverage, and vegetation projection area; (1) Calculate vegetation indices: Vegetation indices reflect the growth status of crops through multi-spectral data; Calculate vegetation indices through the band values of the image. Vegetation indices include the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Optimized Soil Adjusted Vegetation Index (OSAVI), etc., and calculate using existing formulas.

[0070] (2) Calculate plant height: Use the Random Sample Consensus (RANSAC) algorithm to process the segmented crop point cloud data. After processing, the difference between the maximum and minimum values of all points in the point cloud in the vertical direction (such as the Z-axis) is the plant height: Plant height = max(Z) - min(Z); Where: max(Z) represents the maximum value, min(Z) represents the minimum value, and Z is the vertical coordinate of all points in the point cloud data (such as the Z-axis value); The Random Sample Consensus algorithm is an iterative method for estimating the parameters of a mathematical model from a dataset containing outliers, specifically including: Initialize parameters: Set the maximum number of iterations T; Set the inlier threshold for determining whether a point belongs to the model; Set the minimum number of inliers for determining whether the model is valid; Random sampling: Randomly select a part of the samples from all samples (the dataset) as a sample subset; These sample points are used for the estimation of the initial model; Model estimation: Use the sample subset to fit a mathematical model; For example, if the goal is to fit a plane, the points in the sample subset will be used to calculate the parameters of the plane.

[0071] Inlier detection: Substitute all data points into the fitted model and calculate the deviation of each point from the model; If the deviation is less than the set threshold, mark the point as an "inlier"; Otherwise, mark it as an "outlier"; Evaluate the model: Count the number of inliers; If the number of inliers is greater than or equal to the minimum number of inliers, consider the model valid; Otherwise, discard the model and resample; Iterative optimization: Repeat the above steps (random sampling, model estimation, inlier detection, and model evaluation) until the maximum number of iterations T is reached or the best model that meets the conditions is found; In each iteration, record the number of inliers of the current model and retain the model with the most inliers as the best model; Output the best model: At the end of the algorithm, output the model with the most inliers as the final result; If necessary, the model can be refitted using all inliers to improve the accuracy of the model.

[0072] (3) Calculate the surface area: Use the Delaunay triangulation or Poisson reconstruction algorithm to reconstruct the point cloud data into a triangular mesh; The generated triangular mesh consists of multiple triangles, denoted as the triangle set {T j }; For each triangle T j , its vertices are P1, P2, and P3 respectively, and the area calculation formula of the triangle is as follows: ; In the formula, Area ( Tj )represents the area of a triangle, and ∥∥ represents the modulus of a vector; The surface area calculation formula is as follows: ; In the formula, A surface represents the surface area of the crop.

[0073] (4)Calculate the volume: Fill the holes in the vegetation point cloud to ensure the surface is closed; Use Gauss's theorem (divergence theorem) to convert the volume integral into a surface integral, and the formula is as follows: ; In the formula: In the formula: N is the total number of triangular patches; c i is the centroid coordinate of the i-th patch, c i = (V1 + V2 + V3) / 3, where V1, V2, and V3 represent the coordinate vectors of the three vertices of the i-th triangular patch, that is, the positions of the three points that make up the triangle; n i is the unit normal vector of the patch (which needs to point outward); A i represents the area of the patch; This method indirectly solves the volume by calculating the surface properties of the triangular mesh.

[0074] (5)Calculate the canopy cover: Canopy cover (CC) is an index used to estimate the proportion of the covered area by analyzing the projection distribution of the vegetation point cloud on the horizontal plane; The calculation steps are as follows: Point cloud projection and rasterization, including: Define the raster resolution: Select the raster size according to the requirements (such as 0.1m × 0.1m); Projection and statistics: Vertically project the three-dimensional point cloud onto the horizontal plane (X-Y plane) and perform statistics on each raster cell: If there is at least one vegetation point in each raster cell, it is marked as "covered" (value 1); Otherwise, it is "uncovered" (value 0); Generate a two-dimensional raster map: The projection result can be expressed as a binary matrix M ∈ {0, 1}m×n, where m × n is the total number of rasters; Canopy cover calculation; The calculation formula is: ; In the formula: The number of covered rasters is the number of rasters with a value of 1 in the matrix M; The total number of effective rasters is the total number of rasters within the region of interest (ROI) (areas with no data need to be excluded).

[0075] Handle non-closure errors: Non-closed point clouds may lead to an underestimation of the cover due to missing data, and the following compensation methods need to be adopted: Interpolation filling: Interpolate the coverage status of adjacent grids for the void area, such as using nearest neighbor interpolation or Kriging interpolation; Statistical modeling: Assume that the missing area of the point cloud is consistent with the surrounding coverage density, and scale up the number of covered grids proportionally: .

[0076] (6) Calculate the vegetation projected area: The vegetation projected area refers to the projected area of the vegetation on the horizontal plane, which is calculated by analyzing the segmented vegetation image; the calculation steps are as follows: Traverse the segmented vegetation image and count the number of pixels P that all belong to the vegetation; according to the resolution of the image and the shooting scale, determine the true area A corresponding to each pixel (i.e., the ratio), for example, 1 pixel = A square meters; calculate the vegetation projected area: Vegetation projected area = P × A.

[0077] S32. Extract physical parameters: The physical parameter is the canopy temperature feature, and the canopy temperature feature includes the canopy average temperature value Tc, the canopy temperature standard deviation (CTSD), and the canopy temperature coefficient of variation (CTSV); due to the low resolution of the thermal infrared image and only single-band grayscale information, it is difficult to divide plants and soil by simply visually taking the canopy pixel values. Therefore, the K-means clustering algorithm is used to further segment the plot segmentation map of the multi-source crop remote sensing image, and the temperature feature is extracted from the segmented canopy area; The formulas for the canopy temperature standard deviation (CTSD) and the canopy temperature coefficient of variation (CTSV) are as follows: ; ; In the formula: represents the canopy average temperature, represents the temperature of the soybean pixel in the UAV thermal infrared image, is the number of pixels.

[0078] The core of the K-means clustering algorithm is to minimize the sum of the distances between the data points within the cluster and their cluster centers, which is also called the within-cluster sum of squares (WCSS). The iterative process of the algorithm continues until the stopping condition is met, such as the change in the cluster center is less than a certain threshold, or the preset maximum number of iterations is reached. Its core formula is: ; In the formula: is the within-cluster sum of squares (WCSS); is the number of clusters; is the set of data points in the i-th cluster; is a cluster is a data point in the cluster; μ i is the center point of the cluster (i.e., the mean of all points within the cluster); is the data point and its cluster center μ i is the square of the Euclidean distance between them.

[0079] S33. Extract chemical parameters: Chemical parameters include soil chemical elements and vegetation chemical elements. Soil chemical elements include nitrogen, phosphorus, potassium, calcium, magnesium, and organic matter. Vegetation chemical elements include nitrogen, phosphorus, potassium, calcium, and magnesium. The extraction of chemical parameters includes data cleaning, data loading and preprocessing, model training, model tuning and training evaluation, and prediction; S331. Data cleaning: Use the validate_data function, which mainly includes data loading, header row detection, interactive verification, data cleaning, error reporting, and saving the corrected data; Data loading: Load the original spectral data and chemical element content data; Header row detection: Detect whether the header row of the data file meets the expectations; Interactive verification: Verify the integrity of the data through an interactive interface; Data cleaning: First, perform true content cleaning. Clean the last element (true content) of each row of data. First, remove the spaces, and use the regular expression re.sub(r"[^\d.]", "", k_val) to remove non-numeric and non-decimal point characters, keeping only numbers and decimal points. If the cleaned true content is empty, raise an exception. Then, perform spectral data cleaning. Traverse the other values (spectral data) of each row of data, replace the comma (,) with a point (.), and use the regular expression re.sub(r"[^\d.Ee+-]", "", val_str) to remove illegal characters (keeping only numbers, exponential symbols, etc.). If a certain band data is invalid (cannot be converted to a floating number), raise an exception; Error reporting: Perform error handling. If an exception occurs during the cleaning process (such as incorrect data format or missing data), the error information of that row will be recorded in the errors list, including the reason for the error; Saving the corrected data: Save the cleaned data as a new file.

[0080] S332. Data Loading and Preprocessing: Use optimizer.preprocess for preprocessing, including steps such as outlier handling, data augmentation, and data standardization; Outlier Handling: Use the median and median absolute deviation (MAD) to identify and remove outliers. The calculation formula is: ; Then remove outliers that exceed 3 times the MAD.

[0081] Data Augmentation: Augment the dataset by adding noise. The noise is from a normal distribution, and the standard deviation of the noise is 0.05; Generate augmented data: ; Data Standardization: Use StandardScaler to standardize the data so that the mean is 0 and the standard deviation is 1.

[0082] S333. Model Training: S3331. Use the Ridge Regression model for inverting the content of soil chemical elements; Ridge Regression is a variant of linear regression. It solves the problem of multicollinearity that may occur in ordinary least squares regression (OLS) by adding a regularization term to the loss function. When there is a high correlation between features, OLS may lead to too large a variance of the regression coefficients, thus affecting the stability of the model. Ridge Regression suppresses the magnitude of the regression coefficients through regularization, making the model more stable. The goal of Ridge Regression is to minimize the following loss function: ; Where: is the least squares loss, which is the error of the regression model fitting; is the L2 regularization term; λ is the regularization parameter, and λ controls the influence of the regularization term on the model. λ The larger it is, the stronger the influence of regularization; β is the model parameter vector; is the i th feature vector of the sample; β j represents the β th j element of the parameter vector is the number of features; S3332. Use the Random Forest and XGBoost models for inverting the content of vegetation chemical elements, and select the model with the best performance for prediction; [[ID=...]] Random Forest is an ensemble learning method that makes predictions by training multiple decision trees and aggregating the results. The basic unit of Random Forest is the decision tree; multiple subsets are drawn with replacement from the original training dataset (each subset may repeat some samples); for each subset, a decision tree is constructed; for the construction of each tree, the number of features considered at each split is restricted, thereby introducing more randomness and reducing overfitting. The prediction of each tree can be expressed as: ; where, f i is the i-th decision tree, and x is the input feature; the final prediction is the average or vote of the prediction results of all trees: ; where, T is the number of decision trees.

[0083] XGBoost is a powerful ensemble learning method based on the Gradient Boosting algorithm and optimized on this basis. It improves the prediction accuracy by gradually training multiple weak learners (usually decision trees) and aggregating the results with weights. The core idea of gradient boosting is: for each round of iteration, calculate the gradient between the prediction result of the current model and the true value (i.e., the negative gradient of the loss function), and then train a new model based on the gradient information to "correct" the errors of the previous model.

[0084] XGBoost has made many optimizations on the basis of the traditional gradient boosting method: First is regularization. XGBoost adds a penalty term for the complexity of the tree (such as the number and weights of leaves) to the loss function to regularize the loss function, making the model more generalizable and preventing overfitting. Second is pruning. XGBoost adopts post-pruning during the tree generation process. By setting a maximum depth limit, it avoids overgrowth in traditional gradient boosting. Then is column sampling. XGBoost adopts column sampling technology, which can randomly select a subset of features during training to reduce the computational amount and the risk of overfitting. Finally is parallelization. XGBoost improves the training speed through parallel computing (especially column-wise splitting based on data). The loss function formula is as follows: ; where: is the true label of the sample; is the prediction of the model; is the prediction function of the -th tree model; is the regularization term of the model complexity, usually defined as a function of the number of leaves and leaf weights of the tree.

[0085] Compare the performance of the random forest model and the XGBoost model using methods such as cross-validation and test set evaluation, such as accuracy, recall, F1 score, root mean square error, etc.; select the model with better performance for the final prediction.

[0086] S334. Model Tuning and Training Evaluation: Use GridSearchCV to perform grid search on the hyperparameters of each model; select the best parameters through cross-validation, and the scoring criterion is R 2 : ; where: y true is the true value, y pred is the predicted value, is the mean of the true values; Evaluate the stability of the model using cross-validation: Perform 5-fold cross-validation on the training set through cross_val_score, calculate the R 2 for each fold, and output its mean (cv_mean) and standard deviation (cv_std); use the test set for the final evaluation, and calculate the R 2 and root mean square error (RMSE) of the model: ; Plot the learning curve, feature importance, residual analysis, and scatter plot of the true value vs. the predicted value; analyze the performance of the model through visualization; specifically as follows: Learning curve: Use learning_curve to plot the performance of the model at different training set sizes; the R 2 of the training set and the validation set changes with the number of training samples; the learning curve can help determine whether the model is overfitting or underfitting; Feature importance: Use model.feature_importances to obtain the importance scores of each feature; the feature importance reflects the contribution of each feature to the model prediction; Residual analysis: Calculate the residuals (the differences between the true values and the predicted values) and plot the residual distribution; the residual distribution should be close to a normal distribution, and deviation from the normal distribution may indicate problems with the model; Scatter plot of the true value vs. the predicted value: Compare the true values and the predicted values through a scatter plot; ideally, the points should be concentrated near the diagonal; the diagonal represents y = x, that is, the predicted value is equal to the true value.

[0087] Based on the above principles, use hyperspectral images combined with real data to use ridge regression to invert the effect diagram of soil chemical element content, and use random forest and XGBoost to invert the effect diagram of vegetation chemical element content.

[0088] S335. Prediction: Use the trained model to predict new hyperspectral images; the method is as follows: Load the model and normalization object: Load the trained model and the normalizer; load the range of chemical element contents (minimum and maximum values) saved during training; Process the image data: Read the input image in blocks and extract the spectral features of each pixel; perform normalization processing on non-zero valid pixels; use the model to predict the normalized chemical element content values; restore them to real values through inverse normalization; Limit the prediction results: Use np.clip to limit the prediction results within the range of the actual chemical element contents in the training set; the all-zero regions are forced to be assigned 0; Generate the prediction map: Write the predicted values of each image block into a new image file to realize the conversion of hyperspectral data into a chemical element content map, ensuring that the prediction results are both reasonable and have spatial resolution ability.

[0089] Example 2 This example provides a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images, where step S1 is the same as in Example 1; S2. Automatically or semi-automatically segment the images determined to be in the same coordinate system to obtain the plot segmentation map of the multi-source crop remote sensing images; Specifically, semi-automatically segment the areas of bare soil where no crops are planted and the areas of crops with unknown growth conditions; automatically segment the areas of crops with known growth conditions; The method of semi-automatic segmentation is as follows: For the areas of bare soil where no crops are planted, manually segment the plot boundaries; since it is impossible to directly determine the maximum growth range of some remote sensing images of crops with unknown growth conditions, it is necessary to manually segment the crop plots with known maximum growth ranges. Specifically, select the UAV images of the crop plots that have grown to the maximum projected area for manual segmentation. Since the maximum projected area is known, batch cropping of other images is carried out based on this, and then plot image segmentation is performed; The fully automatic segmentation method is as follows: For crops with known growth conditions, taking soybeans as an example, 6 to 8 plants are planted in a plot of 2×1.5 meters. According to the empirical model, the length and width of its maximum growth area will not exceed 2 meters; based on the empirical model, through matching the image coordinates with the geographical location information (GNSS coordinate system) (i.e., step S1), batch segmentation of multi-source sensor images can be achieved.

[0090] Specifically, it includes the following steps: S21. Semi-automatically segment the unplanted bare soil areas and areas with crops of unknown growth conditions. Among them, for the unplanted bare soil areas, manually segment the plot boundaries. For the areas with crops of unknown growth conditions, first select the UAV images of the crop plots that have grown to the maximum projected area for manual segmentation, then batch crop the other images based on the plot boundary range obtained through manual segmentation, and then perform plot image segmentation. S22. For the areas with crops of known growth conditions, use the NDVI, improved EXGR or Otsu algorithm to optimize the vegetation and soil segmentation results to obtain the plot segmentation map of the multi-source crop remote sensing images. For different types of images, if the background complexity is not high and other environments are met, the following algorithms can be used: S221. For RGB images, hyperspectral images, multispectral images or lidar images, use the improved EXGR segmentation algorithm for segmentation. Set a threshold according to the EXGR value, and determine the pixel areas greater than this threshold as vegetation, and the rest as non-vegetation areas. The segmentation formula is as follows: EXGR = 3G - 2.1R - B; Among them, G, R, and B respectively represent the intensity values of the green, red, and blue channels of each pixel in the image.

[0091] S222. For hyperspectral images or multispectral images, the NDVI algorithm can also be used for segmentation. The NDVI calculation formula is as follows: ; In the formula, NDVI represents the normalized difference vegetation index, and NIR and Red respectively represent the intensity values of the near-infrared and red channels of each pixel in the image; By setting a reasonable threshold, accurately distinguish the vegetation and soil areas; determine the areas with NDVI values greater than 0.7 as vegetation, and the NDVI values of the soil are usually 0.1 - 0.25.

[0092] S223. For thermal infrared images, since thermal infrared images usually reflect the temperature in the imaging area, the threshold-based binary classification (Otsu) algorithm is used to segment the thermal infrared images. The temperature is divided into two categories, the higher temperature is vegetation, and the lower temperature is soil. Specifically, it includes: S2231. Calculate the grayscale histogram of the image p i ; Assume that the image has a total of L grayscale levels (such as 0 - 255), p i represents the frequency of the grayscale level i appearing in the image, satisfying ; S2232. Calculate the class probability: Select the threshold t , and divide the pixels into two categories: Background class C0: gray levels 0 to t, with probability ; Foreground class C1: gray levels t + 1 to L - 1, with probability ; S2233. Calculate the class means and the overall mean: Background mean: ; Foreground mean: ; Where: μ 0( t ) and μ 1( t ) represent the means of the background and foreground respectively at the threshold t ; ω 0 and ω 1 are the weights of the background and foreground respectively, i.e., the probabilities of background class pixels; i is the index of the gray level; Calculate the overall mean: ; S2234. Calculate the between-class variance: ; It can be simplified to: ; S2235. Optimize the threshold: Traverse all possible thresholds t and calculate σ 2 (t); Select the 2 that makes σ the largest as the segmentation threshold; ; S2236. Generate a binary image according to and determine the high-temperature regions as vegetation and the low-temperature regions as soil: .

[0093] S3. Extract plant phenotypic parameters, physical parameters, and chemical parameters from the sub-region segmentation map of the multi-source crop remote sensing image; specifically, it includes the following sub-steps: S31. Extract plant phenotypic parameters; Plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy coverage, and vegetation projected area; The vegetation index, surface area, volume, canopy coverage, and vegetation projected area are the same as in Embodiment 1; Calculate the plant height: Use the Random Sample Consensus (RANSAC) algorithm to process the crop point cloud data segmented by the improved EXGR segmentation algorithm in step S221, and the processing steps are the same as in Embodiment 1; Steps S32 and S33 are the same as in Embodiment 1.

[0094] Automatic calculation and extraction of crop phenotypic parameters were carried out according to the methods of Embodiment 1 and Embodiment 2, and the results are as follows: See Figures 2 - 7 , which is the effect diagram of soil chemical element content. From Figures 2 - 7 in turn are the effect diagrams of nitrogen (N), phosphorus (P), potassium (K), calcium (Ca), magnesium (Mg), and organic matter. In the figure, A, B, C, and D are the learning curve, feature importance, residual distribution, and scatter plot of the comparison between the true value and the predicted value, respectively.

[0095] See Figures 8 - 12 , which is the effect diagram of vegetation chemical element content. From Figures 8 - 12 in turn are the effect diagrams of nitrogen (N), phosphorus (P), potassium (K), calcium (Ca), and magnesium (Mg). In the figure, A, B, C, and D are the learning curve, feature importance, residual distribution, and scatter plot of the comparison between the true value and the predicted value of the random forest, respectively. E, F, G, and H are the learning curve, feature importance, residual distribution, and scatter plot of the comparison between the true value and the predicted value of XGBoost, respectively.

[0096] The final predicted effect diagram is shown in Figures 13 - 14 , Figure 13 which is the content distribution of the predicted nitrogen and phosphorus, Figure 14 and

[0097] The key technical points and advantages of the present invention are as follows: Phenotypic parameter extraction: Using image analysis technology to calculate the vegetation index, canopy structure parameters, physical and chemical parameters of crops, and improving the ability to monitor crop growth. Data preprocessing and optimization: By methods such as data cleaning, outlier detection, and standardization, the data quality is improved, noise interference is reduced, and the prediction accuracy of the model is increased. Machine learning and deep learning modeling: Combining machine learning methods such as Random Forest and XGBoost to predict crop phenotypic parameters and improve the accuracy and stability of the model. Hyperspectral inversion analysis: Based on hyperspectral data and combined with ground measured data, the nutrient content of crops (such as nitrogen, phosphorus, potassium, calcium, magnesium, etc.) is inverted to achieve high-precision crop health diagnosis. Model tuning and automated analysis: Using methods such as GridSearchCV hyperparameter optimization, cross-validation, and feature selection to improve the generalization ability and prediction performance of the model. Visualization analysis and decision support: Through data visualization means (such as learning curve, feature importance, residual analysis, scatter comparison diagram, etc.), the interpretability of the results is improved to assist agricultural production management and scientific research analysis. Prediction accuracy control and regional inversion: For pixel-level inversion within the hyperspectral image region, methods such as image block prediction, invalid value elimination, numerical standardization and inverse normalization recovery are used, and combined with the measured calcium content range for prediction value clipping, to improve the accuracy and stability of regional prediction, ensure that the prediction results are true and usable, and provide technical support for the construction of soil nutrient spatial distribution maps and precision agriculture.

[0098] It should be understood that the various forms of processes shown above can be used, with steps reordered, added or deleted. For example, the steps described in the disclosure of the present invention can be executed in parallel, sequentially or in different orders, as long as the desired results of the technical solution disclosed in the present invention can be achieved, and no limitation is imposed herein. The above specific embodiments do not constitute a limitation on the protection scope of the present invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions and improvements made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. An automatic calculation and extraction method for crop phenotypic parameters based on multi-source remote sensing images, characterized in that: Specifically, it includes the following steps: S1. Through georegistration, point cloud file conversion, and matching of specific ground points, all images are determined in the same coordinate system to achieve image registration; S2. Automatically or semi-automatically segment the images determined in the same coordinate system to obtain the sub-region segmentation map of the multi-source crop remote sensing images; S3. Extract plant phenotypic parameters, physical parameters, and chemical parameters from the sub-region segmentation map of the multi-source crop remote sensing images; The plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy coverage, and vegetation projection area; the physical parameter is the canopy temperature characteristic, and the canopy temperature characteristic includes the average canopy temperature value, canopy temperature standard deviation, and canopy temperature variation coefficient; the chemical parameters include soil chemical elements and vegetation chemical elements. The soil chemical elements include nitrogen, phosphorus, potassium, calcium, magnesium, and organic matter, and the vegetation chemical elements include nitrogen, phosphorus, potassium, calcium, and magnesium.

2. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 1, characterized in that: In step S3, the vegetation index is calculated through the band values of the image, and the vegetation index includes the normalized difference vegetation index, enhanced vegetation index, and optimized soil-adjusted vegetation index; The plant height is calculated using the random sample consensus algorithm; The surface area is calculated using the Delaunay triangulation or Poisson reconstruction algorithm; The volume is calculated using Gauss's theorem through the conversion of volume integral and surface integral; Project the three-dimensional vegetation point cloud onto the horizontal plane for rasterization, count the ratio of the number of covered grids to the total number of effective grids, and then perform interpolation filling or proportionally amplify the number of covered grids to obtain the canopy coverage; The vegetation projection area is obtained by counting the number of pixels belonging to the vegetation in the segmented vegetation image and then multiplying by the true area corresponding to each pixel.

3. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 1, characterized in that: In step S3, the K-means clustering algorithm is used to further segment the sub-region segmentation map of the multi-source crop remote sensing images, and the average canopy temperature value, canopy temperature standard deviation, and canopy temperature variation coefficient are extracted from the segmented canopy area; The formulas for the canopy temperature standard deviation and the canopy temperature variation coefficient are as follows: ； ； In the formula: represents the average canopy temperature, is expressed as the temperature of soybean pixels in the UAV thermal infrared image, is the number of pixels; The formula for the K-means clustering algorithm is: ； Where: is the within-cluster sum of squares; is the number of clusters; is the set of data points in the \(i\)th cluster; is a cluster and \(x_{ij}\) is a data point μ i in the cluster and \(\mu_i\) is the center point of the cluster; where \(x_{ij}\) is a data point and \(\mu_i\) is its cluster center μ i and \(d^2(x_{ij},\mu_i)\) is the squared Euclidean distance between them.

4. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 1, characterized in that: The extraction of chemical parameters in step S3 includes data cleaning, data loading and preprocessing, model training, model tuning and training evaluation, and prediction; specifically including: S331. Data cleaning: Use the validate_data function for data loading, header row detection, interactive verification, data cleaning, error reporting, and saving the corrected data; Data cleaning specifically includes: First, perform real content cleaning. Clean the last element of each row of data. First, remove spaces, non-numeric and non-decimal characters, and only keep numbers and decimals. If the real content after cleaning is empty, an exception is thrown; then perform spectral data cleaning. Traverse the other values of each row of spectral data, replace commas with dots, and remove illegal characters. If a certain band data is invalid, an exception is thrown; Error reporting specifically includes: Perform error handling. If an exception occurs during the cleaning process, record the error information of that row in the errors list; S332. Data Loading and Preprocessing: Use optimizer.preprocess for preprocessing, including outlier handling, data augmentation, and data standardization; S333. Model Training: Use the ridge regression model for soil chemical element content inversion; use the random forest and XGBoost models for vegetation chemical element content inversion, and select the best-performing model for prediction; S334. Model Tuning and Training Evaluation: Perform grid search on the hyperparameters of each model using GridSearchCV; select the best parameters through cross-validation, and the scoring criterion is R 2 : ； Wherein: y true is the true value, y pred is the predicted value, is the mean of the true values; S335. Prediction: Use the trained model to predict new hyperspectral images.

5. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 4, characterized in that: The specific steps of S332 include: Outlier Handling: Use the median and median absolute deviation to identify and remove outliers, and the calculation formula is: ； Then remove the outliers that are more than 3 times MAD ; Data Augmentation: Augment the dataset by adding noise to generate augmented data; the noise comes from a normal distribution, and the noise standard deviation is 0.05; Data Standardization: Use StandardScaler to standardize the data so that the mean is 0 and the standard deviation is 1.

6. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 4, characterized in that: In step S333, the goal of the ridge regression model is to minimize the following loss function: ； Wherein: is the least squares loss, i.e., the error of the regression model fitting; is the L2 regularization term; λ is the regularization parameter; β is the model parameter vector; is the i th feature vector of the sample; β j represents the β th j element of the parameter vector; is the number of features; The loss function formula of the XGBoost model is as follows: ； Wherein: is the true label of the sample; is the prediction of the model; is the prediction function of the is the regularization term of the model complexity; Use the methods of cross-validation and test set evaluation to compare the performance of the random forest model and the XGBoost model, and select the better-performing model for the final prediction.

7. A method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images according to claim 4, characterized in that: The step S334 also includes plotting learning curves, feature importances, residual analysis, and scatter plots comparing true values with predicted values; visually analyze the performance of the model; specifically as follows: Learning curve: Use learning_curve to plot the performance of the model with different training set sizes; R of the training set and the validation set 2 varies with the number of training samples; Learning curves can help determine whether the model is overfitting or underfitting; Feature Importance: Use model.feature_importances to obtain the importance scores of each feature; Feature importance reflects the contribution of each feature to the model prediction; Residual Analysis: Calculate the difference between the true value and the predicted value, and plot the residual distribution; Scatter Plot Comparing True Values with Predicted Values: Compare the true value with the predicted value through a scatter plot; ideally, the points are concentrated near the diagonal; the diagonal represents y = x, that is, the predicted value is equal to the true value.

8. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 4, characterized in that: The specific steps of S335 include: Load the model and standardization object: Load the trained model and the standardizer; load the chemical element content range saved during training; Process the image data: Read the input image in blocks, and extract the spectral features of each pixel; perform standardization processing on non-zero valid pixels; use the model to predict the standardized chemical element content value; restore it to the true value through inverse standardization; Limit the prediction results: Use np.clip to limit the prediction results within the range of the actual chemical element content in the training set; the all-zero region is forced to be assigned 0; Generate the prediction map: Write the predicted values of each image block into a new image file, and convert the hyperspectral data into a chemical element content map.

9. The automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images according to claim 1, characterized in that: The specific steps of S2 include: S21. Semi-automatically segment the unplanted bare soil areas and areas of crops with unknown growth conditions; S22. For areas of crops with known growth conditions, use the improved Grounded-SAM segmentation large model or traditional algorithms to perform batch full-automatic segmentation processing on multi-source sensor images, and output the segmentation results; The traditional algorithms are the NDVI algorithm, the improved EXGR algorithm or the Otsu algorithm; the improved EXGR algorithm sets a threshold according to the EXGR value, and determines the pixel regions greater than the threshold as vegetation, and the rest as non-vegetation regions; the segmentation formula is as follows: EXGR = 3G - 2.1R - B; wherein, G, R, and B respectively represent the intensity values of the green, red, and blue channels of each pixel in the image; S23. Optimize the segmentation result, and finally output the plot segmentation map of the optimized multi-source crop remote sensing image.

10. A method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images according to claim 9, characterized in that: In the step S22, specifically including the following sub-steps for using the improved Grounded-SAM segmentation large model: S221. Determine the empirical model parameters: Set the empirical model according to the known crop planting parameters and growth characteristics; S222. Use the improved Grounded-SAM segmentation large model and input the prompt words; receive the image and text prompts through DINO-X, output the bounding boxes in the image that are semantically related to the text, and determine the positions of the plants; S223. Use the bounding boxes output by DINO-X as prompts to further refine each bounding box region to generate accurate segmentation masks; S224. Crop the multi-source sensor image according to the boundary range set by the empirical model; apply the segmentation mask generated in step S223 to the cropped image for plot image segmentation to extract the independent images of each plot; S225. Overlay and display the segmentation result with the original image to generate a new image, or output the segmentation mask image.

Images