Automatic calculation and extraction method of crop phenotypic parameters based on multi-source remote sensing images

Abstract

The present invention relates to the field of crop monitoring technology, and more particularly to a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images. The method includes aligning all images to the same coordinate system through georeferencing, point cloud file conversion, and ground-specific point matching to achieve image registration; automatically or semi-automatically segmenting the images aligned to the same coordinate system to obtain a multi-source crop remote sensing image plot segmentation map; and using an algorithm to extract plant phenotypic, physical, and chemical parameters from the multi-source crop remote sensing image plot segmentation map. Plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy cover, and vegetation projected area; physical parameters include canopy mean temperature, canopy temperature standard deviation, and canopy temperature coefficient of variation; and chemical parameters include soil and vegetation chemical elements such as nitrogen, phosphorus, potassium, calcium, and magnesium. The method has the advantage of enabling large-scale data processing and high-precision regional prediction, enhancing crop growth monitoring capabilities.

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:

[0006] S1. Align all images to the same coordinate system through georeferencing, point cloud file conversion, and ground point matching to achieve image registration.

[0007] S2. Automatically segment or semi-automatically segment the images determined in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image;

[0008] S3. extracting plant phenotypic parameters, physical parameters, and chemical parameters from the plot segmentation map of the multi-source crop remote sensing image;

[0009] Plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy cover and vegetation projected area; physical parameters are canopy temperature characteristics, which include canopy mean temperature, canopy temperature standard deviation and canopy temperature variation coefficient; chemical parameters include soil chemical elements and vegetation chemical elements. Soil chemical elements include nitrogen, phosphorus, potassium, calcium, magnesium and organic matter, and vegetation chemical elements include nitrogen, phosphorus, potassium, calcium and magnesium.

[0010] Preferably, in step S3, the vegetation index is calculated by using the band value of the image, and the vegetation index includes the normalized difference vegetation index, the enhanced vegetation index, and the optimized soil adjusted vegetation index;

[0011] Calculating the plant height using a random sampling consensus algorithm;

[0012] Calculate surface area using Delaunay triangulation or Poisson reconstruction algorithms;

[0013] Use Gauss's theorem to calculate the volume by converting volume integral to surface integral;

[0014] Project the 3D vegetation point cloud onto the horizontal plane and rasterize it. Count the ratio of the number of covered grids to the total number of valid grids. Then interpolate and fill in the number of covered grids or scale up the number of covered grids to obtain the canopy coverage.

[0015] The vegetation projection area is obtained by counting the number of pixels belonging to vegetation in the segmented vegetation image and multiplying it by the actual area corresponding to each pixel.

[0016] Preferably, in step S3, the K-means clustering algorithm is used to further segment the plot segmentation map of the multi-source crop remote sensing image, and the canopy average temperature value, canopy temperature standard deviation and canopy temperature variation coefficient are extracted from the segmented canopy area;

[0017] The formulas for the canopy temperature standard deviation and the canopy temperature coefficient of variation are as follows:

[0018] ;

[0019] ;

[0020] Where: represents the average canopy temperature, It is expressed as the temperature of soybean pixels in the thermal infrared image of the UAV, is the number of pixels;

[0021] The formula of K-means clustering algorithm is:

[0022] ;

[0023] 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 A data point in μ i Is a cluster The center point of is a data point and its cluster center μ i The square of the Euclidean distance between them.

[0024] 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:

[0025] S331. Data cleaning: Use the validate_data function for data loading, header row detection, interactive validation, data cleaning, error reporting, and saving corrected data.

[0026] Data cleaning specifically includes: first, performing real content cleaning, cleaning the last element of each row of data, removing spaces, non-digits, and decimal characters, retaining only digits and decimal points. If the real content after cleaning is empty, an exception is thrown. Then, performing spectral data cleaning, traversing the other values ​​of each row of spectral data, replacing commas with dots, and removing illegal characters. If the data of a certain band is invalid, an exception is thrown.

[0027] Error reporting specifically includes: error handling. If an exception occurs during the cleaning process, the error information of the row is recorded in the errors list;

[0028] S332. Data loading and preprocessing: Use optimizer.preprocess for preprocessing, including outlier handling, data augmentation, and data normalization.

[0029] S333. Model Training: Soil chemical element content is inverted using the ridge regression model; vegetation chemical element content is inverted using the random forest and XGBoost models. The best performing model is selected for prediction.

[0030] 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 use the scoring criteria R 2 :

[0031] ;

[0032] Where: y true is the true value, y pred is the predicted value, is the mean of the true values;

[0033] S335. Prediction: Use the trained model to make predictions about new hyperspectral images.

[0034] Preferably, step S332 specifically includes:

[0035] Outlier handling: Use the median and median absolute deviation to identify and remove outliers. The calculation formula is:

[0036] ;

[0037] Then remove more than 3 times MAD outliers;

[0038] Data augmentation: The dataset is augmented by adding noise to generate augmented data; the noise comes from a normal distribution with a standard deviation of 0.05;

[0039] Use StandardScaler to standardize the data so that the mean is 0 and the standard deviation is 1.

[0040] Preferably, the goal of the ridge regression model in step S333 is to minimize the following loss function:

[0041] ;

[0042] in: is the least squares loss, i.e. the error in fitting the regression model; is the L2 regularization term; λ is the regularization parameter; β is the model parameter vector; For the i The feature vector of each sample; β jRepresents the parameter vector β No. j elements; is the number of features;

[0043] The loss function formula of the XGBoost model is as follows:

[0044] ;

[0045] in: is the true label of the sample; is the model's prediction; It is The prediction function of the tree model; is the regularization term of model complexity;

[0046] The performance of the random forest model and the XGBoost model were compared using cross-validation and test set evaluation methods, and the model with better performance was selected for the final prediction.

[0047] Preferably, step S334 further includes drawing a learning curve, feature importance, residual analysis, and a scatter plot comparing true values ​​and predicted values; and analyzing the performance of the model through visualization; specifically as follows:

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

[0049] 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;

[0050] Residual analysis: calculate the difference between the true value and the predicted value and plot the residual distribution;

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

[0052] Preferably, step S335 specifically includes:

[0053] Load model and standardization objects: load the trained model and standardizer; load the chemical element content range saved during training;

[0054] 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;

[0055] 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;

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

[0057] Preferably, step S2 specifically includes:

[0058] S21. semi-automatically segmenting unplanted bare soil areas and areas with unknown crop growth;

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

[0060] 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:

[0061] EXGR=3G-2.1RB;

[0062] Among them, G, R, and B represent the intensity values ​​of the green, red, and blue channels of each pixel in the image respectively;

[0063] S23. Optimize the segmentation results and finally output an optimized cell segmentation map of the multi-source crop remote sensing image.

[0064] Preferably, in step S22, using the improved Grounded-SAM segmentation model specifically includes the following sub-steps:

[0065] S221. Determine empirical model parameters: Set the empirical model based on known crop planting parameters and growth characteristics;

[0066] S222. Use the improved Grounded-SAM segmentation model and input the prompt word; receive the image and text prompt through DINO-X, output the bounding box in the image that is semantically related to the text, and determine the location of the plant;

[0067] S223. Using the bounding box output by DINO-X as a hint, further refine each bounding box region to generate an accurate segmentation mask;

[0068] S224. Cropping the multi-source sensor image according to the boundary range set by the empirical model; applying the segmentation mask generated in step S223 to the cropped image to perform cell image segmentation and extract an independent image of each cell;

[0069] S225. Overlay the segmentation result with the original image to generate a new image, or output a segmentation mask image.

[0070] Compared with the prior art, the present invention can achieve the following beneficial effects:

[0071] 1. Improve the efficiency of crop phenotyping: This invention adopts automated image segmentation and phenotypic parameter calculation methods, which reduces the dependence on manual measurement and greatly improves the efficiency of field crop phenotyping.

[0072] 2. Applicable to multi-source remote sensing data: This method can be applied to multiple remote sensing data such as RGB, multispectral, hyperspectral, thermal infrared, lidar, etc., and can extract rich crop growth parameters.

[0073] 3. Improve data accuracy and consistency: Through automated segmentation and calculation, human errors are avoided, data consistency is ensured, and measurement accuracy is improved.

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

[0075] 5. Enhance crop growth monitoring capabilities: By extracting multiple vegetation indices, canopy structures, physical parameters, and chemical parameters, crop growth status can be accurately monitored, providing a scientific basis for precision agriculture.

[0076] 6. Improve model generalization capabilities: Use cross-validation, regularization, hyperparameter optimization and other techniques to improve the generalization capabilities of the model and make it applicable to different regions and different crop types.

[0077] 7. Visual analysis supports decision-making: Through 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.

[0078] 8. Achieve high-precision regional predictions: Based on a trained machine learning model, hyperspectral images are processed block by block and predicted pixel by pixel. Combined with standardization and inverse normalization operations, the true values ​​of chemical elements such as calcium are effectively restored. Abnormal outputs are controlled through predicted value clipping, improving the accuracy and reliability of regional inversion and providing data support for crop nutrient monitoring and soil improvement. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 This is a visualization effect diagram of two registered flight files provided according to an embodiment of the present invention.

[0080] Figure 2 This is an effect diagram of the nitrogen content of the soil chemical element provided by an embodiment of the present invention.

[0081] Figure 3 This is an effect diagram of the soil chemical element phosphorus content provided by an embodiment of the present invention.

[0082] Figure 4 This is an effect diagram of the soil chemical element potassium content provided by an embodiment of the present invention.

[0083] Figure 5 This is an effect diagram of the soil chemical element calcium content provided by an embodiment of the present invention.

[0084] Figure 6 This is an effect diagram of the soil chemical element magnesium content provided by an embodiment of the present invention.

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

[0086] Figure 8 These are the results of random forest and XGBoost model analysis of the nitrogen content of chemical elements in vegetation provided by an embodiment of the present invention.

[0087] Figure 9 These are the results of random forest and XGBoost model analysis of the phosphorus content of the chemical element in vegetation provided by an embodiment of the present invention.

[0088] Figure 10 These are the results of random forest and XGBoost model analysis of the potassium content of the chemical element in vegetation provided by an embodiment of the present invention.

[0089] Figure 11 These are the results of random forest and XGBoost model analysis of the calcium content of the chemical element in vegetation provided by an embodiment of the present invention.

[0090] Figure 12 These are the results of random forest and XGBoost model analysis of magnesium content in vegetation according to an embodiment of the present invention.

[0091] Figure 13 3 is a predicted nitrogen and phosphorus content distribution diagram provided according to an embodiment of the present invention.

[0092] Figure 14 1 is a diagram of the predicted potassium, calcium, and magnesium content distribution provided according to an embodiment of the present invention. DETAILED DESCRIPTION

[0093] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. In the following description, identical modules are denoted by the same reference numerals. Where identical reference numerals are used, their names and functions are also identical. Therefore, their detailed description will not be repeated. To make the objectives, technical solutions, and advantages of the present invention more clearly understood, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0094] Example 1

[0095] 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:

[0096] S1. Align all images to the same coordinate system through georeferencing, point cloud file conversion, and ground point matching. This involves the following steps:

[0097] S11. Extract UAV flight data, which includes GNSS observation files (.obs), navigation data files (NAV files), event log files (EVENT files), binary log files (BIN files), and marker files (MRK files). The details are as follows:

[0098] GNSS observation file (.obs): A standard data file that primarily stores raw measurement data from the Global Navigation Satellite System. It follows 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.

[0099] Navigation data files (NAV files): Store information about the flight path, GPS coordinates, altitude, and other relevant data that can be used for georeferencing of images, flight analysis, or troubleshooting;

[0100] Event log files (EVENT files): Record various events that occurred during the flight, such as when the drone took off, when it captured an image, or when it reached a waypoint. They can help understand the progress of the flight and identify any problems that may have occurred during the mission.

[0101] Binary log files (BIN files): Contain detailed information about the drone's flight, including telemetry data, sensor readings, and system information, which can be used for in-depth analysis and troubleshooting of the drone's performance;

[0102] Marker file (MRK file): records the relevant information of the photo point, including:

[0103] Column 1: Photo point serial number, i.e. the serial number of the photo record information stored in this folder; Column 2: GPS seconds of the week, indicating the time when the photo was taken, expressed in GPS seconds of the week (TOW); Column 3: GPS week, indicating the GPS week at the time of the photo; Column 4: Compensation value in the north direction, in millimeters, with north as positive; Column 5: Compensation value in the east direction, in millimeters, with east as positive; Column 6: Compensation value in the elevation direction, in millimeters, with down as positive; Column 7: Longitude after compensation; Column 8: Latitude after compensation; Column 9: Ellipsoid height; Column 10: Positioning standard deviation (north); Column 11: Positioning standard deviation (east); Column 12: Positioning standard deviation (elevation); Column 13: Positioning status.

[0104] S12. Assigning the geographic location information collected from the flight data extracted in step S11 to the image to achieve image georegistration;

[0105] 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 calibration.

[0106] In order to realize the registration process, it is necessary to use the flight data file collected during the UAV flight. The file contains key data such as geographic location information, flight attitude information, and timestamp. Through this data, the spatial trajectory of the photography center can be solved, and then the georegistration of the image can be realized. The core formula for solving the spatial trajectory of the photography center is based on the collinearity equation and space-time interpolation. The formula is as follows:

[0107] ;

[0108] Where: λ represents the proportional 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 coordinates of the ground point; (X0, Y0, Z0) represents the coordinates of the photography center (GPS position to be solved); R represents the exterior orientation rotation matrix composed of the drone attitude angles (ω, φ, κ); f represents the camera focal length; (x, y) represents the image point coordinates;

[0109] The specific steps include the following:

[0110] S121. Time synchronization calibration: Ensure that the image capture time and the GPS track log time are consistent in the same time system (UTC or local time zone). If there is a system clock deviation, the image time needs to be compensated:

[0111] t corrected =t image +Δt;

[0112] Where: t corrected is the image capture time after calibration; t image is the original capture time of the image (extracted from the metadata of the image file, i.e., the DateTimeOriginal field in the EXIF ​​information); Δt is the system clock deviation compensation, which can be obtained by capturing a reference point of known position at the same time and inversely calculating the time offset.

[0113] S122. Trajectory interpolation algorithm: using the calibrated image capture time t corrected As the interpolation time point, cubic spline interpolation is used to output the exact position (X, Y, Z) of each image capture moment, that is, the WGS84 coordinates; for the trajectory segment index j, the interpolation formula is:

[0114] ;

[0115] in: t Indicates the time point at which interpolation is required; j Indicates the trajectory segment index; t j Indicates the trajectory segment index j The starting time point; S j ( t ) represents the trajectory segment index j In time t The position function value of ; Indicates the trajectory segment index j In time t The second derivative of the position function; Indicates the trajectory segment index j +1 in time t j The second derivative of the position function; coefficient a j 、b j 、c j d j Solve by boundary conditions;

[0116] The simplified linear interpolation formula is as follows:

[0117] ;

[0118] Where: X ( t ) represents the X coordinate value at time t; X {k} Indicates time t {k} The X coordinate value at time ; X {k+1} Indicates time t {k+1} The X coordinate value at time ; t Indicates the time point at which the X coordinate value needs to be estimated; t {k} and t {k+1} Represents two time points with known X coordinate values;

[0119] The Y and Z coordinate formulas are the same as the X coordinate formulas:

[0120] ;

[0121] ;

[0122] Brief Principle: A drone's GPS trajectory is typically recorded at a frequency of 1 Hz, while a camera's shutter speed can reach multiple frames per second. This means that the sampling frequency of GPS position information is lower than the image capture frequency. To obtain the precise position (X, Y, Z) at the time of each image capture, the GPS trajectory must be interpolated. The cubic spline interpolation used in this step is suitable for high-precision aerial surveys and can ensure the continuity of position and velocity.

[0123] S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S202 to the CGCS2000 coordinate system (which complies with Chinese requirements); use the seven-parameter Bursa model for coordinate conversion:

[0124] ;

[0125] Where: [X CGCS ,Y CGCS ,Z CGCS ] represents the coordinates in the converted CGCS2000 coordinate system; [Δ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, which represents the scale difference between the two coordinate systems; represents the rotation matrix, 、 、 are rotation parameters, representing the rotation angles around the X, Y, and Z axes respectively.

[0126] S124. Use the Earth Gravitational Model 2008 (EGM2008) to perform elevation corrections, converting GPS-measured altitudes (based on the ellipsoid) to orthometric altitudes (based on the geoid). The elevation correction formula is:

[0127] H ortho = H GPS - N geoid ;

[0128] Where: H ortho Indicates positive height; H GPS Indicates the altitude measured by GPS; N geoid Represents the geoid elevation anomaly calculated through model interpolation.

[0129] S13. Image data generation: Convert the LAS point cloud file with color information into a GeoTIFF image file with RGB bands. This specifically includes the following steps:

[0130] S131. Read the LAS file and extract the point cloud data; check whether the LAS file contains color information (RGB channels); if there is no color information, throw an exception to ensure the validity of subsequent operations;

[0131] Specifically, the laspy.read method is used to read the input LAS file; the point cloud data includes spatial coordinates x, y and color channel data red, green, blue.

[0132] S132. Calculate the grid size: Based on the spatial extent of the point cloud data (minimum and maximum coordinate values) 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:

[0133] ;

[0134] ;

[0135] Where: max x ,min x ,max y ,min y They are the maximum and minimum coordinate values ​​of the point cloud data respectively; resolution is the resolution;

[0136] The resolution of the grid determines the actual spatial range represented by each pixel. The resolution is the side length of each pixel in space, usually in meters. Based on the given resolution, the width and height of the grid can be calculated.

[0137] S133. Binned statistics: Use the binned_statistic_2d function to assign the data of each RGB channel to the calculated grid;

[0138] binned_statistic_2d is a SciPy function used to perform two-dimensional binning statistics; the formula is as follows: statistic=statistic method (data in grid cells);

[0139] Among them, x and y are the spatial coordinates of the point cloud data, data is the corresponding RGB color channel data; the bin edge is specified by the bins parameter, here is x edges and y edges ; statistic specifies the calculation method used for statistics, such as mean (mean), max (maximum), min (minimum), etc.;

[0140] The present invention uses mean as the default statistical method, that is, the color values ​​in each grid cell are averaged to generate a mean grid.

[0141] S134.RGB Merge: Merge the three binned color channels (red, green, blue) into a three-dimensional RGB stack (RGB stack ), where each pixel contains the color values ​​of the three channels of red, green and blue, making it easy to generate color images;

[0142] The specific operations are 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 RGB stack: fill the binned red, green, and blue channel data into RGB respectively. stack The first three positions of the third dimension, such as the red channel data filled into rgb stack [:,:,0], green channel data is filled into rgb stack [:,:,1], blue channel data is filled into rgb stack [:,:,2]. Before merging, ensure that all channels have the same data type, usually an unsigned 8-bit integer (uint8), representing a color value range of 0 to 255.

[0143] S135. Color scaling: Determine whether to scale the color value from 16 bits to 8 bits based on actual needs; if scale color =True, the color value is scaled from 16 bits to 8 bits by dividing the value of each color channel by 65535 (the maximum value of 16 bits) and multiplying it by 255 (the maximum value of 8 bits); the scaled color value is updated to the RGB stack, ensuring that the data type of all color channels is uint8; this step determines whether scaling is required based on the actual situation;

[0144] The linear scaling formula is as follows:

[0145] ;

[0146] Where: scaled color Represents the scaled color value; color represents the original color value, which comes from the color information in the LAS file.

[0147] S136. Generate 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;

[0148] A GeoTIFF file contains three bands (RGB), with the appropriate data type set for each band; and is georeferenced using a given geographic transformation (Affine) and coordinate system (crs); the metadata of a GeoTIFF file includes band information, data type, geographic transformation, and coordinate system.

[0149] S14. Use specific ground points for image matching to achieve accurate registration of two or more UAV remote sensing images. The specific steps are as follows:

[0150] S141. Identify and select fixed ground specific points in the image, reassign fixed coordinate points to the selected ground specific points, and use the fixed coordinate points as reference points for image registration;

[0151] Ground specific points can be selected from fixed points in the field, such as houses, roads, water pipes, etc. in the image; fixed coordinate points are reassigned to the ground specific points so that two or more images can be matched, so that all remote sensing images of the detection area are kept in the same coordinate system to achieve accurate registration;

[0152] 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. Pixel coordinates, arbitrary coordinates without projection information, or pseudo coordinates are not acceptable. There are no strict constraints on the image to be registered, but if coordinate information is missing, manually select at least three points of the same name.

[0153] The Image Registration Workflow tool in ENVI software can automatically, quickly, and accurately register two images with different geometric positions. This automated, accurate, and rapid image registration workflow integrates complex parameter setting steps into a unified panel, enabling rapid and accurate automatic registration of images with little or no manual intervention. Data includes a pair of data with the same resolution but different imaging times, and a pair of panchromatic and multispectral data with different resolutions and the same imaging time. The registration tool supports the following data formats: ENVI, TIFF, NITF, JPEG2000, JPEG, Esri® raster layer, Geodatabase raster, and Web Services.

[0154] S143. Feature Extraction and Matching: Automatically detect and match two images using the cross correlation method ( I 1 and I 2) The corresponding feature points in the calculation formula are as follows:

[0155] ;

[0156] Where: CC ( I 1, I 2) Representing an image I 1 and I 2; M and N represent the number of rows and columns of the image respectively; I 1 ( i , j )and I 2 ( i , j ) represent images respectively I 1 and I 2 in position ( i , j ); cross-correlation measures the similarity of two images by calculating the product of their pixel values.

[0157] Cross Correlation (CC) is a common method for measuring the similarity between two images. It works by calculating the similarity of pixel intensities between the images to determine whether or how the images are aligned. Cross Correlation is generally used for images with similar morphology or features (for example, two images taken by the same sensor or of the same type, such as optical images versus optical images).

[0158] The algorithm logic is as follows:

[0159] (1) Calculate the matching window: select the image I 1 A fixed area (called a window) on the image I 2. Scan the same size area;

[0160] (2) Sliding window calculation 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;

[0161] (3) Maximum value position: Finally, the position with the largest cross-correlation value is selected, which is the estimated position of the best alignment of the two images.

[0162] S144. Geometric Transformation Calculation: Using a polynomial transformation model, the positions of image pixels are measured and aligned with their corresponding geographic locations. The polynomial transformation coefficients are then fitted using the paired points. After fitting, the entire image is geometrically corrected using the obtained polynomial coefficients. After correction, the accuracy of the geometric correction is verified by comparing the corrected image with the reference image.

[0163] Polynomial transformation usually uses polynomial functions to fit the position transformation of image pixels. The most common form is two-dimensional polynomial transformation, which is expressed as follows:

[0164] ;

[0165] ;

[0166] Where: X and Y are the row and column coordinates of the original image pixels respectively; X’ and Y’ are the row and column coordinates of the corrected image pixels, a i and b i (i=0,1,2,3,4,5,…,n) are the polynomial coefficients;

[0167] The geometric correction step maps the positions of image pixels to their real geographic locations in the real world, achieving precise alignment of the images.

[0168] In this step, for the geometric correction of the lidar point cloud data, the orthographic projection of the lidar point cloud data is first extracted, and then the orthographic projection image is aligned with the original lidar point cloud image; the orthographic projection converts the three-dimensional model into an orthographic angle projection image, which retains the geographic location information; because the two images have the same viewing angle at the orthographic angle and the surface features at this angle are consistent, the alignment process mainly focuses on assigning accurate geographic location information to point cloud data at different heights to ensure the consistency of the point cloud data and the orthographic projection image in geographic space.

[0169] S145. Resampling: Apply cubic convolution to resample the image to match the spatial resolution of the reference image.

[0170] Cubic convolution interpolation is a high-order interpolation method that uses more pixel values ​​around the target point to calculate the interpolation result, thereby producing a smoother image. It performs weighted calculations by selecting multiple pixels around the target pixel (usually a 4×4 area containing 16 pixels). Each neighboring pixel is assigned a different weight value based on its relative position to the target pixel to achieve neighborhood pixel weighting.

[0171] The specific processing process is as follows:

[0172] S1451. Determine a resampling area and collect neighboring pixel values: Select an area (4×4 area) around the target pixel, which contains multiple pixels for calculating the interpolated value; collect the values ​​of 16 pixels in the selected area as the neighboring pixels of the target pixel, and use them to calculate the interpolated value of the target pixel;

[0173] S1452. Calculate weights: Use a cubic convolution kernel to calculate weights for each neighboring pixel. Based on the calculated weights, weight the neighboring pixel values. Ultimately, the interpolated value of the target pixel (x', y') is the weighted average of the neighboring pixel values, expressed as:

[0174] ;

[0175] Where: I ( x ', y ') indicates that the image is at position ( x ', y ')'s pixel value; I ( x + i , y + j ) is the value of the neighborhood pixel; i andj is the relative position offset within the neighborhood; is the horizontal distance x’ With neighboring points x + i The weight between is the longitudinal distance y’ With neighboring points y + j The weight between .

[0176] The convolution kernel calculates weights based on the distance between pixels and the target pixel. It uses the Catmull-Rom interpolation kernel, which is a smooth cubic function that effectively reduces jagged effects and blurring during interpolation. Pixels closer to the target pixel receive higher weights, while pixels farther away receive lower weights.

[0177] S1453. Generate a 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 that there is no obvious aliasing or blurring; if necessary, adjust the parameters of the convolution kernel to optimize the result;

[0178] S1454. Output the final image: Output the resampled image as the final result, ready for subsequent image analysis or application;

[0179] Cubic convolution interpolation uses more pixel value information and therefore usually produces higher quality results when enlarging images; however, its computational complexity is also relatively high.

[0180] S146. Use the multiple dynamic link display feature of ENVI to compare the baseline image and the registered image to evaluate the registration accuracy.

[0181] The above method can be used to register images acquired by different sensors (such as RGB, multispectral, hyperspectral, and lidar);

[0182] For example, RGB images are three-dimensional, multispectral images are five- or more-dimensional, and hyperspectral images are multi-dimensional. LiDAR images are also three-dimensional, and these differ from sensor to sensor. RGB, hyperspectral, and multispectral images all have RGB or pseudo-RGB bands, but pixel accuracy varies between sensors. As mentioned above, images from the same sensor across different missions also present issues. Therefore, it is necessary to address the issues of varying image resolution, geographic location information, and image formats. The aforementioned methods are used to assign geographic information to RGB, multispectral, hyperspectral, and thermal infrared images. For the registration of geographic location information between RGB images and point clouds, it is necessary to first extract the orthographic projection of the lidar point cloud data and then register the orthographic projection image with the original lidar point cloud image; the orthographic projection converts the three-dimensional model into an orthographic angle projection image, which retains the geographic location information; since the two images have the same viewing angle at the orthographic angle and the surface features at this angle are consistent, the registration process mainly focuses on assigning accurate geographic location information to point cloud data at different heights to ensure the consistency of the point cloud data and the orthographic projection image in geographic space.

[0183] Figure 1 The following figure shows a visualization of two flight files (including deviations in flight turning points and altitude). The left image is a three-dimensional Flight Trajectories (3D) diagram, showing the spatial trajectories of the two flights (Flight 1 and Flight 2), with the horizontal axis representing the UTM X coordinate, the vertical axis representing the UTM Y coordinate, and the vertical axis representing the flight altitude. The right image is a two-dimensional Flight Trajectories (2D) diagram, which removes the altitude dimension to more clearly illustrate the ground coverage paths of the two flights. As can be seen from the figure, for flights with the same planned route, there are differences in altitude and position during actual execution.

[0184] For the registration of thermal infrared images of different dates, the hyperspectral images are registered with thermal infrared images to achieve minimal error.

[0185] S2. Automatically segment or semi-automatically segment the images determined in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image;

[0186] Specifically, unplanted bare soil areas and areas with unknown crop growth are semi-automatically segmented; areas with known crop growth are automatically segmented;

[0187] The semi-automatic segmentation method is as follows: For unplanted bare soil areas, the plot boundaries are manually segmented. Since the maximum growth range of some crops with unknown growth cannot be directly determined from remote sensing images, crop plots with known maximum growth ranges need to be manually segmented. Specifically, drone-captured images of crop plots that have grown to their maximum projected area are selected for manual segmentation. Since the maximum projected area is known, it is used as a basis for batch cropping other images and then performing plot segmentation.

[0188] The fully automatic segmentation method is as follows: For crops with known growth, such as soybeans, 6 to 8 plants are planted in a 2×1.5 meter plot. According to the empirical model, the length and width of their maximum growth area will not exceed 2 meters. Based on the empirical model, batch segmentation of multi-source sensor images can be achieved by matching image coordinates with geographic location information (GNSS coordinate system) (i.e., step S1).

[0189] The specific steps include:

[0190] S21. semi-automatically segmenting unplanted bare soil areas and areas with unknown crop growth;

[0191] For unplanted bare soil areas, the cell boundaries were manually segmented. For areas with unknown crop growth, drone images of crop plots with the largest projected area were first selected for manual segmentation. The cell boundaries obtained through manual segmentation were then used as a benchmark for batch cropping of other images before performing cell segmentation.

[0192] S22. For crop areas with known growth conditions, perform batch, fully automatic segmentation processing on multi-source sensor images using the improved Grounded-SAM segmentation model and output segmentation results. This specifically includes the following sub-steps:

[0193] S221. Determine empirical model parameters: Set the empirical model based on known crop planting parameters and growth characteristics; this step provides prior knowledge for subsequent automated segmentation and reduces unnecessary calculations;

[0194] For example, if 6 to 8 soybean plants are planted in a 2×1.5 meter plot, according to the empirical model, the length and width of the maximum growth area will not exceed 2 meters.

[0195] S222. Using the improved Grounded-SAM segmentation model, input the prompt word "plant"; receive the image and text prompt via DINO-X, output a bounding box in the image that is semantically related to the text, and determine the location of the plant;

[0196] DINO-X is an open-source object detector capable of detecting relevant objects in images based on arbitrary free-form text prompts. Since this large model performs segmentation based on prompt words, simply inputting the word "plant" will cause the model to detect and segment crops within the image plot based on this vocabulary. The cropped plot image contains only land and crops, making segmentation extremely simple. This step aims to guide the model through natural language prompts to detect and segment the target crop, and to quickly locate the target crop area within the image.

[0197] S223. Using the bounding box output by DINO-X as a hint, each bounding box region is further refined to generate an accurate segmentation mask.

[0198] 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 to perform cell image segmentation and extract an independent image for each cell.

[0199] S225. The segmentation result is superimposed on the original image to generate a new image, or a segmentation mask image is output, which includes the confidence or visual outline of the segmented area;

[0200] In a specific embodiment, the RGB image of the field plot is segmented using the Grounded-SAM model to extract the area of ​​the soybean plant; the segmentation mask (for example, green represents vegetation and red represents soil) is superimposed on the original RGB image to generate a new image, in which the vegetation area of ​​the original image is covered with green and the soil area is covered with red.

[0201] S23. Optimize the segmentation results, including constructing a height matrix to determine the shape of the vegetation area, calculating the maximum rectangular area row by row using a monotone stack algorithm, and calculating the maximum growth projection area of ​​each plot based on an empirical model. Finally, output an optimized plot segmentation map of the multi-source crop remote sensing image. This specifically includes the following sub-steps:

[0202] S231. Construct a height matrix: For the connected domain (foreground is 1, background is 0) in the independent image (binary image) segmented and extracted in step S224, traverse each pixel row by row and calculate the number of consecutive 1s at each position upward to form a height matrix height. Determine the shape and size of the connected domain (vegetation area) of each cell. The height matrix formula is as follows:

[0203] ;

[0204] Where i represents the row index in the image, and j represents the column index in the image;

[0205] S232. Calculate the maximum rectangular area row by row: For each row of the height array, use the monotone stack algorithm to find the maximum rectangular area of ​​the current row. The specific process is as follows:

[0206] Initialize the stack: The stack is used to store the indexes of columns with increasing heights. The stack is initially empty, and the heights from the bottom to the top of the stack increase in sequence.

[0207] Traverse each element of the current row: When the stack is not empty and the current height is less than the top height of the stack, pop the top element of the stack and calculate the area of ​​the rectangle with the height height[top] as the height;

[0208] ;

[0209] Among them, h is the height of the pop-up, i is the current index, s top-1 The new top index of the stack after popping; push the current index into the stack;

[0210] Process the remaining elements in the stack: After the traversal is completed, perform similar operations on the remaining elements in the stack to calculate the possible area;

[0211] Update the global maximum area: After processing each row, record the maximum rectangular area of ​​the row; traverse all rows and take the maximum value of all rows as the global maximum area (indicating the maximum rectangular area of ​​the vegetation area in all rows).

[0212] The simplified formula is as follows:

[0213] Height matrix: accumulate the height of consecutive 1s in each column;

[0214] Maximum area of ​​a single line: Use the monotone stack to find the maximum width, area = height × width;

[0215] Global maximum area: take the maximum value among all rows;

[0216] S233. Calculate the maximum bounding rectangle based on the empirical model: For each connected domain of the segmented independent cell image, calculate its maximum bounding rectangle; this step can determine the maximum growth projection area of ​​each cell, providing a basis for subsequent analysis.

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

[0218] S3. Extracting plant phenotypic parameters, physical parameters, and chemical parameters from the plot segmentation map of the multi-source crop remote sensing image; specifically comprising the following sub-steps:

[0219] S31. Extracting plant phenotypic parameters;

[0220] Plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy cover, and vegetation projected area;

[0221] (1) Calculation of vegetation index: The vegetation index reflects the crop growth status through multispectral data. The vegetation index is calculated based on the band values ​​of the image. Vegetation indices include the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI), the Optimized Soil Adjusted Vegetation Index (OSAVI), etc., and are calculated using existing formulas.

[0222] (2) Calculation of plant height: The segmented crop point cloud data is processed using the Random Sample Consensus (RANSAC) algorithm. 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);

[0223] 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);

[0224] The random sampling consensus algorithm is an iterative method for estimating the parameters of a mathematical model from a data set containing outliers. It includes:

[0225] Initialization parameters: set the maximum number of iterations T; set the inlier threshold to determine whether a point belongs to the model; set the minimum number of inliers to determine whether the model is valid;

[0226] Random sampling: Randomly select some samples from all samples (data sets) as a sample subset; these sample points are used to estimate the initial model;

[0227] Model estimation: Fitting a mathematical model using a subset of samples; for example, if the goal is to fit a plane, the points in the sample subset are used to calculate the parameters of the plane.

[0228] 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, the point is marked as an "inlier"; otherwise, it is marked as an "outlier";

[0229] Evaluate the model: count the number of inliers; if the number of inliers is greater than or equal to the minimum number of inliers, the model is considered valid; otherwise, discard the model and resample;

[0230] 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, the number of inliers of the current model is recorded, and the model with the most inliers is retained as the best model;

[0231] Output the best model: At the end of the algorithm, the model with the most inliers is output as the final result; if necessary, all inliers can be used to refit the model to improve the accuracy of the model.

[0232] (3) Calculate surface area: Use Delaunay triangulation or Poisson reconstruction algorithm to reconstruct the point cloud data into a triangular mesh; the generated triangular mesh consists of multiple triangles, which are recorded as a triangle set {T j}; For each triangle T j , whose vertices are P1, P2 and P3 respectively. The area of ​​the triangle is calculated as follows:

[0233] ;

[0234] Where, Area( T j ) represents the area of ​​a triangle, ∥∥ represents the modulus of a vector;

[0235] The surface area is calculated as follows:

[0236] ;

[0237] Where, A surface Represents the surface area of ​​the crop.

[0238] (4) Calculate volume: Fill the holes in the vegetation point cloud to ensure surface closure; use Gauss's theorem (divergence theorem) to convert the volume integral into the surface integral. The formula is as follows:

[0239] ;

[0240] In the formula: N is the total number of triangles; c i is the centroid coordinate of the i-th patch, c i =(V1+V2+V3) / 3, V1, V2, V3 represent the coordinate vectors of the three vertices of the i-th triangle patch, that is, the positions of the three points that make up the triangle; n i is the unit normal vector of the patch (needs to point outward); A i Represents the area of ​​the patch;

[0241] This method solves the volume indirectly by calculating the surface properties of the triangular mesh.

[0242] (5) Calculation of canopy cover: Canopy cover (CC) is an indicator that estimates the coverage area ratio by analyzing the projection distribution of vegetation point clouds on the horizontal plane. The calculation steps are as follows:

[0243] Point cloud projection and rasterization, including:

[0244] Define grid resolution: select the grid size as required (e.g. 0.1m×0.1m);

[0245] Projection and statistics: Project the 3D point cloud vertically onto the horizontal plane (XY plane) and perform statistics on each grid cell: if there is at least one vegetation point in each grid cell, it is marked as "covered" (value 1); otherwise, it is marked as "uncovered" (value 0);

[0246] Generate a two-dimensional grid map: The projection result can be expressed as a binary matrix M∈{0,1}m×n, where m×n is the total number of grids;

[0247] Canopy coverage calculation; the calculation formula is:

[0248] ;

[0249] Where: the number of covered grids is the number of grids with a value of 1 in the matrix M; the total number of valid grids is the number of all grids in the region of interest (ROI) (excluding areas with no data).

[0250] Dealing with non-closed errors: Non-closed point clouds may underestimate coverage due to missing data, and the following compensation methods are required:

[0251] Interpolation filling: Use the coverage status of neighboring grids to interpolate the hole area, such as nearest neighbor interpolation or Kriging interpolation;

[0252] Statistical modeling: Assuming that the point cloud missing area has the same coverage density as the surrounding area, the number of coverage grids is increased proportionally:

[0253] .

[0254] (6) Calculation of vegetation projection area: The vegetation projection area refers to the projection area of ​​vegetation on the horizontal plane, which is calculated by analyzing the segmented vegetation image. The calculation steps are as follows:

[0255] Traverse the segmented vegetation image and count the number of pixels P belonging to vegetation. Determine the actual area A (ratio) corresponding to each pixel based on the image resolution and shooting scale, for example, 1 pixel = A square meters. Calculate the vegetation projection area:

[0256] Vegetation projected area = P × A.

[0257] S32. Extracting physical parameters: The physical parameters are canopy temperature characteristics, including the canopy mean temperature (Tc), canopy temperature standard deviation (CTSD), and canopy temperature coefficient of variation (CTSV). Because thermal infrared images have low resolution and only provide single-band grayscale information, it is difficult to distinguish between plants and soil by simply visually extracting canopy pixel values. Therefore, a K-means clustering algorithm is used to further segment the plots of multi-source crop remote sensing images, and temperature characteristics are extracted from the segmented canopy regions.

[0258] The formulas for canopy temperature standard deviation (CTSD) and canopy temperature coefficient of variation (CTSV) are as follows:

[0259] ;

[0260] ;

[0261] Where: represents the average canopy temperature, It is expressed as the temperature of soybean pixels in the thermal infrared image of the UAV, is the number of pixels.

[0262] The core of the K-means clustering algorithm is to minimize the sum of the distances between data points within a cluster and their cluster centers, also known as the Within-Cluster Sum of Squares (WCSS). The algorithm iterates until a stopping condition is met, such as the change in cluster center being less than a certain threshold, or until a preset maximum number of iterations is reached. The core formula is:

[0263] ;

[0264] Where: 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 A data point in μ i Is a cluster The center point of the cluster (i.e., the mean of all points in the cluster); is a data point and its cluster center μ i The square of the Euclidean distance between them.

[0265] S33. Chemical parameter extraction: Chemical parameters include soil chemical elements (nitrogen, phosphorus, potassium, calcium, magnesium, and organic matter) and vegetation chemical elements (nitrogen, phosphorus, potassium, calcium, and magnesium). Chemical parameter extraction includes data cleaning, data loading and preprocessing, model training, model tuning and training evaluation, and prediction.

[0266] S331. Data cleaning: Use the validate_data function, which mainly includes data loading, header row detection, interactive validation, data cleaning, error reporting, and saving corrected data;

[0267] Data loading: loading raw spectrum data and chemical element content data;

[0268] Header row detection: Check whether the header row of the data file meets expectations;

[0269] Interactive verification: Verify data integrity through an interactive interface;

[0270] Data cleaning: First, perform real content cleaning. Clean the last element (real content) of each row of data, remove spaces, and use the regular expression re.sub(r"[^\d.]", "", k_val) to remove non-digits and decimal characters, retaining only digits and decimal points. If the real content after cleaning is empty, an exception is thrown. Then, perform spectral data cleaning. Iterate over the other values ​​of each row of data (spectral data), replace commas (,) with dots (.), and use the regular expression re.sub(r"[^\d.Ee+-]", "", val_str) to remove illegal characters (retaining only digits, exponential symbols, etc.). If a band of data is invalid (cannot be converted to a floating number), an exception is thrown.

[0271] Error reporting: Performs error handling. If an exception occurs during the cleaning process (such as incorrect or missing data format), the error information of the row will be recorded in the errors list, including the cause of the error.

[0272] Save the corrected data: Save the cleaned data as a new file.

[0273] S332. Data loading and preprocessing: Use optimizer.preprocess for preprocessing, including outlier handling, data augmentation, and data normalization.

[0274] Outlier processing: Use the median and median absolute deviation (MAD) to identify and remove outliers. The calculation formula is:

[0275] ;

[0276] Outliers with more than 3 times the MAD are then removed.

[0277] Data augmentation: The dataset is augmented by adding noise. The noise comes from a normal distribution with a standard deviation of 0.05. The augmented data is generated:

[0278] ;

[0279] Data standardization: Use StandardScaler to standardize the data so that the mean is 0 and the standard deviation is 1.

[0280] S333. Model training:

[0281] S3331. Inversion of soil chemical element content using the Ridge Regression model;

[0282] Ridge regression is a variant of linear regression that addresses the multicollinearity problem that may be encountered 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 cause the variance of the regression coefficients to be too large, 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:

[0283] ;

[0284] in: is the least squares loss, that is, the error of the regression model fitting; is the L2 regularization term; λ is the regularization parameter, which controls the influence of the regularization term on the model. λ The larger it is, the stronger the regularization effect; β is the model parameter vector; For the i The feature vector of each sample; β j Represents the parameter vector β No. j elements; is the number of features;

[0285] S3332. Invert vegetation chemical element content using Random Forest and XGBoost models, selecting the best performing model for prediction.

[0286] Random Forest is an ensemble learning method that trains multiple decision trees and aggregates the results to make predictions. The basic unit of a random forest is the decision tree. Multiple subsets are extracted from the original training dataset with replacement (each subset may have duplicate samples). For each subset, a decision tree is constructed. When constructing each tree, the number of features considered at each split is limited, thereby introducing more randomness and reducing overfitting. The prediction of each tree can be expressed as:

[0287] ;

[0288] in, f i is the i-th decision tree, x is the input feature; the final prediction is the average or vote of the prediction results of all trees:

[0289] ;

[0290] Where T is the number of decision trees.

[0291] XGBoost is a powerful ensemble learning method based on and optimized for the gradient boosting algorithm. It improves prediction accuracy by progressively training multiple weak learners (typically decision trees) and weightedly aggregating the results. The core idea of ​​gradient boosting is that for each iteration, the gradient between the current model's prediction and the true value (i.e., the negative gradient of the loss function) is calculated. A new model is then trained based on this gradient information to "correct" the errors of the previous model.

[0292] XGBoost has made many optimizations based on the traditional gradient boosting method: the first is regularization. XGBoost adds a penalty term for the complexity of the tree (such as the number and weight of leaves) in the loss function, regularizes the loss function, makes the model more generalizable, and prevents overfitting. The second is pruning. XGBoost uses post-pruning (pruning) in the tree generation process. By setting a maximum depth limit, it avoids excessive growth in traditional gradient boosting. Then there is column sampling. XGBoost uses column sampling technology to randomly select feature subsets during training to reduce the amount of calculation and the risk of overfitting. Finally, there is parallelization. XGBoost improves training speed through parallel calculations (especially column-based splitting of data). The loss function formula is as follows:

[0293] ;

[0294] in: is the true label of the sample; is the model's prediction; It is The prediction function of the tree model; It is a regularization term for the model complexity and is usually defined as a function of the leaf size and leaf weights of the tree.

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

[0296] 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 use R as the scoring standard. 2 :

[0297] ;

[0298] Where: y true is the true value, y pred is the predicted value, is the mean of the true values;

[0299] Use cross-validation to evaluate the stability of the model: perform 5-fold cross-validation on the training set through cross_val_score and calculate the R of each fold 2 , and output its mean (cv_mean) and standard deviation (cv_std); use the test set for final evaluation and calculate the R of the model 2 And the root mean square error (RMSE):

[0300] ;

[0301] Draw learning curves, feature importance, residual analysis, and scatter plots comparing true values ​​and predicted values; analyze the performance of the model through visualization; the details are as follows:

[0302] 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;

[0303] 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;

[0304] Residual analysis: Calculate the residuals (the difference between the true value and the predicted value) and plot the residual distribution; the residual distribution should be close to normal distribution. Deviation from normal distribution may indicate problems with the model;

[0305] Scatter plot of true values ​​versus predicted values: Use a scatter plot to compare the true values ​​with the predicted values; ideally, the points should be concentrated near the diagonal line; the diagonal line represents y=x, that is, the predicted value is equal to the true value.

[0306] Based on the above principles, hyperspectral images are combined with real data, and ridge regression is used to invert the effect map of soil chemical element content, and random forest and XGBoost are used to invert the effect map of vegetation chemical element content.

[0307] S335. Prediction: Use the trained model to make predictions about new hyperspectral images; the method is as follows:

[0308] Load model and standardization objects: load the trained model and normalizer; load the chemical element content range (minimum and maximum values) saved during training;

[0309] 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;

[0310] 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;

[0311] Generate prediction map: Write the predicted value of each image block into a new image file to realize the conversion of hyperspectral data into chemical element content map, ensuring that the prediction result is both reasonable and has spatial resolution.

[0312] Example 2

[0313] This embodiment provides a method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images, wherein step S1 is the same as in embodiment 1;

[0314] S2. Automatically segment or semi-automatically segment the images determined in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image;

[0315] Specifically, unplanted bare soil areas and areas with unknown crop growth are semi-automatically segmented; areas with known crop growth are automatically segmented;

[0316] The semi-automatic segmentation method is as follows: For unplanted bare soil areas, the plot boundaries are manually segmented. Since the maximum growth range of some crops with unknown growth cannot be directly determined from remote sensing images, crop plots with known maximum growth ranges need to be manually segmented. Specifically, drone-captured images of crop plots that have grown to their maximum projected area are selected for manual segmentation. Since the maximum projected area is known, it is used as a basis for batch cropping other images and then performing plot segmentation.

[0317] The fully automatic segmentation method is as follows: For crops with known growth, such as soybeans, 6 to 8 plants are planted in a 2×1.5 meter plot. According to the empirical model, the length and width of their maximum growth area will not exceed 2 meters. Based on the empirical model, batch segmentation of multi-source sensor images can be achieved by matching image coordinates with geographic location information (GNSS coordinate system) (i.e., step S1).

[0318] The specific steps include:

[0319] S21. Semi-automatically segment unplanted bare soil areas and areas with unknown crop growth. For unplanted bare soil areas, manually segment the plot boundaries. For areas with unknown crop growth, manually segment the drone-captured crop plots that have grown to their maximum projected area. Then, based on the manually segmented plot boundaries, batch crop the remaining images and perform plot segmentation.

[0320] S22. For crop areas with known growth, optimize the vegetation and soil segmentation results using NDVI, modified EXGR, or the Otsu algorithm to obtain a plot segmentation map for multi-source crop remote sensing images. For different image types, if the background complexity is not high, the following algorithms can be used:

[0321] S221. For RGB images, hyperspectral images, multispectral images, or lidar images, use the improved EXGR segmentation algorithm for segmentation. A threshold is set based on the EXGR value, and pixel areas with pixels greater than the threshold are defined as vegetation, while the rest are defined as non-vegetation areas. The segmentation formula is as follows:

[0322] EXGR=3G-2.1RB;

[0323] Among them, G, R, and B represent the intensity values ​​of the green, red, and blue channels of each pixel in the image, respectively.

[0324] S222. For hyperspectral or multispectral images, NDVI algorithm can also be used for segmentation. The NDVI calculation formula is as follows:

[0325] ;

[0326] Where NDVI represents the normalized vegetation index, NIR and Red represent the intensity values ​​of the near-infrared and red channels of each pixel in the image, respectively;

[0327] By setting reasonable thresholds, vegetation and soil areas can be accurately distinguished; areas with NDVI values ​​greater than 0.7 are identified as vegetation, and the NDVI value of soil is usually 0.1~0.25.

[0328] S223. For thermal infrared images, since they typically reflect the temperature within the imaging area, a threshold-based binary classification (Otsu) algorithm is used to segment the thermal infrared images. The temperature is divided into two categories: high temperature vegetation and low temperature soil. Specifically, the following are performed:

[0329] S2231. Calculate the grayscale histogram of the image p i ; Assume that the image has L gray levels (such as 0~255), p i Represents the frequency of gray level i in the image, satisfying ;

[0330] S2232. Calculating Class Probability: Selecting a Threshold t , the pixels are divided into two categories:

[0331] Background class C0: gray level 0~t, probability is ;

[0332] Foreground class C1: gray level t+1~L-1, probability is ;

[0333] S2233. Calculate class mean and overall mean:

[0334] Background mean: ;

[0335] Outlook mean: ;

[0336] Where: μ 0( t )and μ 1( t ) represent the threshold values t The mean of the background and foreground; ω 0 and ω 1 is the weight of background and foreground, i.e. the probability of background class pixels; i is the index of gray level;

[0337] Calculate the population mean: ;

[0338] S2234. Calculate the between-class variance: ;

[0339] can be simplified to: ;

[0340] S2235. Optimize threshold: traverse all possible thresholds t and calculate σ 2 (t); choose σ 2 (t) Maximum As the segmentation threshold;

[0341] ;

[0342] S2236. According to Generate a binary map, identifying areas with high temperatures as vegetation and areas with low temperatures as soil:

[0343] .

[0344] S3. Extracting plant phenotypic parameters, physical parameters, and chemical parameters from the plot segmentation map of the multi-source crop remote sensing image; specifically comprising the following sub-steps:

[0345] S31 extract plant phenotypic parameters; plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy cover and vegetation projected area; vegetation index, surface area, volume, canopy cover and vegetation projected area with Example 1;

[0346] Calculating plant height: Using the Random Sample Consensus (RANSAC) algorithm, the crop point cloud data segmented by the improved EXGR segmentation algorithm in step S221 is processed, and the processing steps are the same as those in Example 1;

[0347] Step S32 and step S33 are the same as those in Example 1.

[0348] The crop phenotypic parameters were automatically calculated and extracted according to the methods of Examples 1 and 2. The results are as follows: Figure 2-Figure 7 , is the effect diagram of soil chemical element content, from Figure 2-Figure 7 The following are the effect diagrams of nitrogen (N), phosphorus (P), potassium (K), calcium (Ca), magnesium (Mg) and organic matter (Organic Matter). A, B, C, and D in the figure are the learning curve, feature importance, residual distribution, and scatter plots comparing true values ​​and predicted values.

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

[0350] The final prediction effect diagram is shown in Figure 13-14 , Figure 13 To predict the distribution of nitrogen and phosphorus content, Figure 14 is the predicted distribution of potassium, calcium and magnesium contents.

[0351] The key technical features and advantages of this invention lie in the following aspects: Phenotypic Parameter Extraction: Utilizing image analysis technology to calculate crop vegetation indices, canopy structure parameters, and physical and chemical parameters, this improves crop growth monitoring capabilities. Data Preprocessing and Optimization: Through methods such as data cleaning, outlier detection, and standardization, data quality is improved, noise interference is reduced, and model prediction accuracy is enhanced. Machine Learning and Deep Learning Modeling: Combining machine learning methods such as Random Forest and XGBoost, crop phenotypic parameters are predicted, improving model accuracy and stability. Hyperspectral Inversion Analysis: Based on hyperspectral data and combined with ground-truth data, crop nutrient content (such as nitrogen, phosphorus, potassium, calcium, and magnesium) is inverted to achieve high-precision crop health diagnosis. Model Tuning and Automated Analysis: Utilizing methods such as GridSearchCV hyperparameter optimization, cross-validation, and feature selection, this improves model generalization and predictive performance. Visualization Analysis and Decision Support: Data visualization tools (such as learning curves, feature importance, residual analysis, and scatter plots) enhance the interpretability of results and assist in agricultural production management and scientific research analysis. Prediction accuracy control and regional inversion: For pixel-level inversion within the hyperspectral image area, image block prediction, invalid value elimination, numerical standardization and inverse normalization recovery methods are used, combined with the measured calcium content range to tailor the predicted values, improve the accuracy and stability of regional predictions, ensure that the prediction results are real and usable, and provide technical support for the construction of soil nutrient spatial distribution maps and precision agriculture.

[0352] It should be understood that the various forms of processes shown above can be used to reorder, add or delete steps. For example, the steps disclosed in the present invention can be executed in parallel, sequentially or in different orders, as long as the expected results of the technical solutions disclosed in the present invention can be achieved, and this document does not limit them here. The above specific implementation methods do not constitute limitations on the scope of protection of the present invention. It should be understood by those skilled in the art 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 principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images, characterized by: The specific steps include: S1. Align all images to the same coordinate system through georeferencing, point cloud file conversion, and ground point matching to achieve image registration. S2. Automatically segment or semi-automatically segment the images determined in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image; S3. Extract plant phenotypic, physical, and chemical parameters from plot segmentation maps of multi-source crop remote sensing images. Plant phenotypic parameters include vegetation index, plant height, surface area, volume, canopy cover, and vegetation projected area. Physical parameters include canopy temperature characteristics, including canopy mean temperature, canopy temperature standard deviation, and canopy temperature coefficient of variation. Chemical parameters include soil chemical elements (nitrogen, phosphorus, potassium, calcium, magnesium, and organic matter) and vegetation chemical elements (nitrogen, phosphorus, potassium, calcium, and magnesium). Chemical parameter extraction includes data cleaning, data loading and preprocessing, model training, model tuning and training evaluation, and prediction. Specifically, it includes: S331. Data cleaning: Use the validate_data function for data loading, header row detection, interactive validation, data cleaning, error reporting, and saving corrected data. S332. Data loading and preprocessing: Use optimizer.preprocess for preprocessing, including outlier handling, data augmentation, and data normalization. S333. Model Training: Soil chemical element content is inverted using the ridge regression model; vegetation chemical element content is inverted using the random forest and XGBoost models. The best performing model is selected 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 use the scoring criteria R 2 : ； Where: 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 make predictions about new hyperspectral images; specifically: 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.

2. The method for automatically calculating and extracting 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 using the band values ​​of the image, where the vegetation index includes the normalized difference vegetation index, the enhanced vegetation index, and the optimized soil adjusted vegetation index; Calculating the plant height using a random sampling consensus algorithm; Calculate surface area using Delaunay triangulation or Poisson reconstruction algorithms; Use Gauss's theorem to calculate the volume by converting volume integral to surface integral; Project the 3D vegetation point cloud onto the horizontal plane and rasterize it. Count the ratio of the number of covered grids to the total number of valid grids. Then interpolate and fill in the number of covered grids or scale up the number of covered grids to obtain the canopy coverage. The vegetation projection area is obtained by counting the number of pixels belonging to vegetation in the segmented vegetation image and multiplying it by the actual area corresponding to each pixel.

3. The method for automatically calculating and extracting 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 plot segmentation map of the multi-source crop remote sensing image, and the canopy mean 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 coefficient of variation are as follows: ； ； Where: represents the average canopy temperature, It is expressed as the temperature of soybean pixels in the thermal infrared image of the UAV, is the number of pixels; The formula of 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 A data point in μ i Is a cluster The center point of is a data point and its cluster center μ i The square of the Euclidean distance between them.

4. The method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images according to claim 1, characterized in that: In step S3, data cleaning specifically includes: first, performing real content cleaning, cleaning the last element of each row of data, removing spaces, non-digits and decimal characters, and retaining only digits and decimal points. If the real content after cleaning is empty, an exception is thrown; then performing spectral data cleaning, traversing other values ​​of each row of spectral data, replacing commas with dots, and removing illegal characters. If a band of data is invalid, an exception is thrown; The error report specifically includes: error handling. If an exception occurs during the cleaning process, the error information of the row is recorded in the errors list.

5. The 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 S332 specifically includes: Outlier handling: Use the median and median absolute deviation to identify and remove outliers. The calculation formula is: ； Then remove more than 3 times MAD outliers; Data augmentation: The dataset is augmented by adding noise to generate augmented data; the noise comes from a normal distribution with a standard deviation of 0.05; Data standardization: Use StandardScaler to standardize the data so that the mean is 0 and the standard deviation is 1.

6. The method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images according to claim 4, characterized in that: The goal of the ridge regression model in step S333 is to minimize the following loss function: ； in: is the least squares loss, i.e. the error in fitting the regression model; is the L2 regularization term; λ is the regularization parameter; β is the model parameter vector; For the i The feature vector of each sample; β j Represents the parameter vector β No. j elements; is the number of features; The loss function formula of the XGBoost model is as follows: ； in: is the true label of the sample; is the model's prediction; It is The prediction function of the tree model; is the regularization term of model complexity; The performance of the random forest model and the XGBoost model were compared using cross-validation and test set evaluation methods, and the model with better performance was selected for the final prediction.

7. The 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 drawing a learning curve, feature importance, residual analysis, and a scatter plot comparing true values ​​and predicted values; and visually analyzing 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.

8. The method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images according to claim 1, characterized in that: The 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.

9. The method for automatically calculating and extracting crop phenotypic parameters based on multi-source remote sensing images according to claim 8, characterized in that: In step S22, the improved Grounded-SAM segmentation model is used to specifically include the following sub-steps: S221. Determine empirical model parameters: Set the empirical model based on known crop planting parameters and growth characteristics; S222. Use the improved Grounded-SAM segmentation model and input the prompt word; receive the image and text prompt through DINO-X, output the bounding box in the image that is semantically related to the text, and determine the location of the plant; S223. Using the bounding box output by DINO-X as a hint, further refine each bounding box region to generate an accurate segmentation mask; S224. Cropping the multi-source sensor image according to the boundary range set by the empirical model; applying the segmentation mask generated in step S223 to the cropped image to perform cell image segmentation and extract an independent image of each cell; S225. Overlay the segmentation result with the original image to generate a new image, or output a segmentation mask image.

Images