Automatic plot segmentation method for crop remote sensing images based on geographic prior information

Abstract

The present invention relates to the field of crop monitoring technology, and in particular to a method for automatic segmentation of crop remote sensing images based on geographic prior information. The method comprises: determining all images in the same coordinate system through geographic registration, point cloud file conversion, and ground specific point matching to achieve image registration; automatically segmenting or semi-automatically segmenting the images determined in the same coordinate system to obtain a segmentation map of the multi-source crop remote sensing image; semi-automatically segmenting unplanted bare soil areas and areas with unknown crop growth; for crop areas with known growth, using an improved Grounded‑SAM segmentation model to perform batch fully automatic segmentation processing on multi-source sensor images and output segmentation results; optimizing the segmentation results and finally outputting an optimized segmentation map of the multi-source crop remote sensing image. The advantages are: precise registration, efficient segmentation, multi-sensor compatibility, and automated processing to enhance remote sensing monitoring.

Description

Technical Field

[0001] The present invention relates to the technical field of crop monitoring, and in particular to a method for automatically segmenting crop plots based on geographical prior information using remote sensing images. Background Art

[0002] With the continuous development of agricultural informatization and precision agriculture, remote sensing technology is increasingly being used in agriculture, particularly in field monitoring and analysis. Traditional drone remote sensing images captured by different sensors exhibit significant discrepancies, making it impossible to accurately segment them using the same image processing method. For example, patent application number CN118334522A, published on July 12, 2024, and titled "Method and Apparatus for Automatic Extraction of High-Throughput Field Crop Plots Based on Low-Altitude Remote Sensing," uses an image processing method to segment land and then crops the image based on estimated values. This introduces certain errors and fails to consider geographic location information. Furthermore, it is impossible to use the same coordinates for batch segmentation across all images captured in a single experiment. Another example is patent application number CN116469019A, published on July 21, 2023, titled "Method for Automatic Field Plot Segmentation Based on a Fully Convolutional Neural Network," which also fails to consider geographic location information and only uses a network for batch segmentation. Segmenting each individual plot in a large field image typically requires a significant amount of computation.

[0003] Common problems with existing technologies include the following aspects: (1) Sensor differences: Images acquired by different sensors (such as RGB cameras, thermal infrared cameras, hyperspectral cameras, lidar, etc.) have large differences in spectral characteristics, resolution and data format. Traditional methods are difficult to apply to the segmentation of multiple sensor images at the same time. (2) Inefficiency of manual segmentation: In agricultural remote sensing detection, field plot planning and segmentation are essential steps. At present, many studies and practical applications still rely on manual segmentation or use software such as ArcGIS to plan plots and then perform batch segmentation. However, due to the mismatch of geographic location information, the manual matching and segmentation process is time-consuming and labor-intensive, and prone to errors. (3) Processing of images of different time series: In long-term monitoring, remote sensing images of different time series need to be segmented and analyzed. Existing methods cannot effectively handle the differences between images of different time series, making it difficult to achieve automation and batch processing. Summary of the Invention

[0004] In order to solve the above problems, the present invention provides a method for automatically segmenting plots of crop remote sensing images based on geographic prior information.

[0005] The present invention aims to provide a method for automatically segmenting crop plots from remote sensing images based on geographic prior information, 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 to be in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image; specifically comprising the following steps:

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

[0009] S22. For crop areas with known growth conditions, use the improved Grounded-SAM segmentation model to perform batch, fully automatic segmentation processing on multi-source sensor images and output the segmentation results.

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

[0011] Preferably, in step S21, the boundaries of the cells are manually segmented for the unplanted bare soil areas; for the areas with unknown growing crops, the drone-photographed images of the crop cells that have grown to the maximum projected area are first selected for manual segmentation, and then the cell boundary range obtained by manual segmentation is used as a reference to batch crop other images, and then cell image segmentation is performed.

[0012] Preferably, step S22 specifically includes the following sub-steps:

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

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

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

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

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

[0018] Preferably, the prompt word in step S222 is plant;

[0019] The image output in step S225 includes the confidence level or the visual outline of the segmented region.

[0020] Preferably, the optimization processing in step S23 specifically includes: constructing a height matrix to determine the shape of the vegetation area, using the monotone stack algorithm to calculate the maximum rectangular area row by row, and calculating the maximum growth projection area of ​​each cell based on the empirical model, and finally outputting the cell segmentation map of the optimized multi-source crop remote sensing image.

[0021] Preferably, step S23 specifically includes the following sub-steps:

[0022] S231. For the connected domain in the independent image segmented and extracted in step S224, traverse each pixel row by row, calculate the number of consecutive 1s at each position upward, form a height matrix height, and determine the shape and size of the connected domain of each cell;

[0023] S232. 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:

[0024] S2321. 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.

[0025] S2322. Traverse each element in the current row: If the stack is not empty and the current height is less than the top of the stack, pop the top element of the stack and calculate the area of ​​the rectangle with height height[top] as the height.

[0026] S2323. Process the remaining elements in the stack: After the traversal is completed, step S2322 is executed on the remaining elements in the stack to calculate the possible area;

[0027] S2324. 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;

[0028] S233. For each connected domain of the segmented independent cell image, calculate its maximum outer rectangle and output the maximum growth projection area of ​​each cell.

[0029] Preferably, the height matrix formula is as follows:

[0030] ;

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

[0032] The formula for the area of ​​a rectangle is as follows:

[0033] ;

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

[0035] Preferably, step S1 specifically includes the following steps:

[0036] S11. Extract UAV flight data, including GNSS observation files, navigation data files, event log files, binary log files, and marker files;

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

[0038] S13. Convert the LAS point cloud file with color information into a GeoTIFF image file with RGB bands;

[0039] S14. Use specific ground points for image matching to achieve accurate registration of at least two UAV remote sensing images.

[0040] Preferably, step S12 specifically includes the following sub-steps:

[0041] S121. Time synchronization calibration: Ensure that the image capture time and the GPS track log time are consistent in the same time system; if there is a system clock deviation, the image time needs to be compensated;

[0042] S122. Trajectory interpolation algorithm: Use the calibrated image capture time as the interpolation time point, use cubic spline interpolation or linear interpolation, and output the precise position (X, Y, Z) of each image capture time, i.e., WGS84 coordinates;

[0043] S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S122 to the CGCS2000 coordinate system;

[0044] S124. Use the EGM2008 model to perform elevation correction and convert the GPS-measured altitude to orthometric height.

[0045] The step S13 specifically includes the following sub-steps:

[0046] S131. Read the LAS file and extract the point cloud data; check whether the LAS file contains color information for the three RGB channels; if there is no color information, an exception is thrown;

[0047] S132. Calculate the size of the grid: Calculate the number of columns and rows of the grid based on the spatial extent of the point cloud data and the preset resolution;

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

[0049] S134.RGB Merge: Merges the three binned color channels into a three-dimensional RGB stack, where each pixel contains the color values ​​of the red, green, and blue channels, facilitating the generation of a color image.

[0050] 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 and multiplying by 255; the scaled color value is updated to the RGB stack;

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

[0052] The specific sub-steps of step S14 include:

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

[0054] S142. Designate one image as the reference image and the other as the image to be registered; the reference image includes standard map coordinates or RPC information;

[0055] S143. Automatically detect and match corresponding feature points in the two images using a cross-correlation method;

[0056] S144. Geometric Transformation Calculation: Using a polynomial transformation model, the positions of image pixels are measured and aligned with their corresponding geographic locations. Polynomial transformation coefficients are then fitted using the paired points. Once fitted, the resulting polynomial transformation coefficients are used to perform geometric correction on the entire image.

[0057] S145. Applying cubic convolution interpolation to resample the image to be registered to adapt to the spatial resolution of the reference image;

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

[0059] Preferably, the calculation formula for compensating the image time in step S121 is as follows:

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

[0061] Where: t corrected is the image capture time after calibration; t image is the original capture time of the image; Δt is the system clock deviation compensation;

[0062] When cubic spline interpolation is used in step S122, for the trajectory segment index j, the interpolation formula is:

[0063] ;

[0064] in: t Indicates the time point where 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;

[0065] When linear interpolation is used, the linear interpolation formula is as follows:

[0066] ;

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

[0068] The Y and Z coordinates have the same formula as the X coordinate:

[0069] ;

[0070] ;

[0071] The step S123 uses the seven-parameter Bursa model to perform coordinate transformation:

[0072] ;

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

[0074] The elevation correction formula in step S124 is:

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

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

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

[0078] 1. Accurate registration and batch segmentation: Integrating UAV flight parameters and geospatial information, it achieves registration of images from different sensors at different times, keeping all remote sensing images in the same coordinate system, thereby improving segmentation accuracy and supporting batch segmentation.

[0079] 2. Reduce the amount of calculation and improve computing efficiency: First use geographic information to crop, then divide the area, avoiding complex calculations on the entire large image and reducing computing costs;

[0080] 3. Applicable to a variety of remote sensing sensors: Supports multiple data sources such as RGB, hyperspectral, multispectral, lidar, and thermal infrared, solving the problem of incompatible processing of data from different sensors;

[0081] 4. Compatible with automatic and semi-automatic segmentation: For crops with known growth conditions, fully automatic segmentation is performed using empirical models. For crops with unknown growth conditions, the maximum growth projection area is first determined, and then batch cropping is performed, improving flexibility.

[0082] 5. Improved segmentation method: Using the Grounded-SAM model combined with DINO-X object detection to achieve precise target recognition and improve segmentation accuracy;

[0083] 6. Efficient processing of multi-time series images: Through the calculation of the maximum growth projection area, batch segmentation of full-time series remote sensing images can be achieved, reducing manual intervention and improving the degree of automation of remote sensing monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 This is a visualization display effect diagram of two registered flight files in a method for automatic plot segmentation of crop remote sensing images based on geographic prior information provided by an embodiment of the present invention.

[0085] Figure 2 This is the image segmentation effect after the cell is automatically segmented according to the embodiment of the present invention.

[0086] Figure 3 It is the result of semi-automatic segmentation of lidar point cloud files, visible light images, and thermal infrared images collected on different dates provided by an embodiment of the present invention.

[0087] Figure 4 This is the automatic segmentation effect of the visible light image provided by the embodiment of the present invention. DETAILED DESCRIPTION

[0088] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. In the following description, identical modules are denoted by identical reference numerals. In the case of identical reference numerals, their names and functions are also identical. Therefore, their detailed description will not be repeated.

[0089] In order to make the purpose, 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. It should be understood that the specific embodiments described herein are only used to explain the present invention and do not constitute a limitation of the present invention.

[0090] The present invention provides a method for automatically segmenting plots of crop remote sensing images based on geographic prior information, which specifically includes the following steps:

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

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

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

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

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

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

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

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

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

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

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

[0102] ;

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

[0104] The specific steps include the following:

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

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

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

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

[0109] ;

[0110] in: t Indicates the time point where 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;

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

[0112] ;

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

[0114] The Y and Z coordinates are the same as the X coordinate:

[0115] ;

[0116] ;

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

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

[0119] ;

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

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

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

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

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

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

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

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

[0128] ;

[0129] ;

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

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

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

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

[0134] 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.;

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

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

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

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

[0139] The linear scaling formula is as follows:

[0140] ;

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

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

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

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

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

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

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

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

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

[0150] ;

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

[0152] 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).

[0153] The algorithm logic is as follows:

[0154] (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;

[0155] (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;

[0156] (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.

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

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

[0159] ;

[0160] ;

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

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

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

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

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

[0166] The specific processing process is as follows:

[0167] S1451. Determine the resampling area and collect neighborhood pixel values: Select an area (4×4 area) around the target pixel, which contains multiple pixels used to calculate the interpolated value; collect the values ​​of 16 pixels in the selected area to calculate the interpolated value of the target pixel;

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

[0169] ;

[0170] 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 and j 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 .

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

[0172] S1453. Generate a resampled image: Repeat steps S4051 and S4052 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;

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

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

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

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

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

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

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

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

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

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

[0183] 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).

[0184] The specific steps include:

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

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

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

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

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

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

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

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

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

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

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

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

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

[0198] ;

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

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

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

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

[0203] ;

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

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

[0206] 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).

[0207] The simplified formula is as follows:

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

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

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

[0211] 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. Figure 2 Demonstrates the image segmentation effect.

[0212] For different types of images, if the background complexity is not high, we are not limited to the improved Grounded-SAM segmentation model. We can also use NDVI, EXGR or Otsu algorithm to optimize the vegetation and soil segmentation results to obtain the cell segmentation map of multi-source crop remote sensing images; the details are as follows:

[0213] For RGB images, hyperspectral images, multispectral images or lidar images, the improved EXGR segmentation algorithm is used for segmentation. A threshold is set according to the EXGR value, and the pixel area greater than the threshold is determined as vegetation, and the rest is non-vegetation area. The segmentation formula is as follows:

[0214] EXGR=3G-2.1RB;

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

[0216] For hyperspectral images or multispectral images, NDVI algorithm can also be used for segmentation; the NDVI calculation formula is as follows:

[0217] ;

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

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

[0220] For thermal infrared images, since they usually 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 included:

[0221] (1) Calculate the grayscale histogram p of the image 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 ;

[0222] (2) Calculate the category probability: Select the threshold t and divide the pixels into two categories:

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

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

[0225] (3) Calculate the class mean and the overall mean:

[0226] Background mean: ;

[0227] Outlook mean: ;

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

[0229] Calculate the population mean: ;

[0230] (4) Calculate the between-class variance: ;

[0231] can be simplified to: ;

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

[0233] ;

[0234] (6) According to Generate a binary map, identifying areas with high temperatures as vegetation and areas with low temperatures as soil:

[0235] ;

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

[0237] The key technical points of the present invention include the following aspects:

[0238] 1. Cropping method design: By batch cropping the registered images, multi-source remote sensing image segmentation can be performed by manually or automatically setting the coordinate system once;

[0239] 2. Automatic segmentation based on deep learning: Using the Grounded-SAM model, crop keywords are input for automatic target detection and segmentation;

[0240] 3. Batch segmentation based on projection area: Batch segmentation of multiple time series images of the same crop is performed using the maximum projection area calculation method;

[0241] 4. Optimize computational efficiency: Use the monotone stack algorithm to solve the maximum rectangular area and improve the computational speed of cell segmentation;

[0242] 5. Improve traditional algorithms: For RGB and hyperspectral / multispectral images, the EXGR (3G-2.1RB) algorithm and improved NDVI are used for segmentation; for thermal infrared images, the Otsu adaptive threshold method is used for vegetation and soil classification.

[0243] The advantage of the present invention is that high-precision segmentation of field plot images is achieved through image processing using modern deep learning models (Grounded-SAM and DINO-X), which not only improves the accuracy and robustness of the segmentation, but also enhances the adaptability and practicality of the present invention.

[0244] Figure 3 The semi-automatic segmentation effect is demonstrated. That is, after image correction, the UTM coordinates of each cell vertex are first obtained, and the UTM coordinates are used to achieve the segmentation results of visible light, radar point cloud files, and thermal infrared images collected on different dates. Figure 4 The automatic segmentation effect, i.e., the automatic identification of cells in visible light images, and the complete process and results of mask generation and vegetation image extraction are demonstrated.

[0245] It should be understood that the various forms of the processes shown above can be used to reorder, add, or delete steps. For example, the steps described in the present disclosure can be performed in parallel, sequentially, or in a different order, as long as the desired results of the technical solutions disclosed in the present disclosure can be achieved. This is not limited herein.

[0246] The above specific embodiments do not limit the scope of protection of the present invention. Those skilled in the art will appreciate that various modifications, combinations, sub-combinations, and substitutions may be made based on design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention are intended to be included within the scope of protection of the present invention.

Claims

1. A method for automatic plot segmentation of crop remote sensing images based on geographic prior information, 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 to be in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image; specifically comprising the following steps: S21. semi-automatically segmenting unplanted bare soil areas and areas with unknown crop growth; 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: 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. The segmentation result is superimposed on the original image to generate a new image, or the segmentation mask image is output; S23. Optimize the segmentation results and ultimately output an optimized cell segmentation map of the multi-source crop remote sensing image; the optimization process specifically includes: constructing a height matrix to determine the shape of the vegetation area, using a monotone stack algorithm to calculate the maximum rectangular area row by row, and calculating the maximum growth projection area of ​​each cell based on an empirical model, and ultimately outputting an optimized cell segmentation map of the multi-source crop remote sensing image.

2. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 1, characterized in that: In step S21, the cell boundaries are manually segmented for the unplanted bare soil area; for the area with unknown crop growth, the drone-captured images of the crop cells that have grown to the maximum projected area are first selected for manual segmentation, and then the cell boundary range obtained by manual segmentation is used as a reference to batch crop other images, and then the cell image segmentation is performed.

3. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 1, characterized in that: The prompt word in step S222 is plant; The image output in step S225 includes the confidence level or the visual outline of the segmented region.

4. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 1, characterized in that: The step S23 specifically includes the following sub-steps: S231. For the connected domain in the independent image segmented and extracted in step S224, traverse each pixel row by row, calculate the number of consecutive 1s at each position upward, form a height matrix height, and determine the shape and size of the connected domain of each cell; S232. 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: S2321. 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. S2322. Traverse each element in the current row: If the stack is not empty and the current height is less than the top of the stack, pop the top element of the stack and calculate the area of ​​the rectangle with height height[top] as the height. S2323. Process the remaining elements in the stack: After the traversal is completed, step S2322 is executed on the remaining elements in the stack to calculate the possible area; S2324. 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; S233. For each connected domain of the segmented independent cell image, calculate its maximum outer rectangle and output the maximum growth projection area of ​​each cell.

5. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 4, characterized in that: The height matrix formula is as follows: ； Where i represents the row index in the image, and j represents the column index in the image; The formula for the area of ​​a rectangle is as follows: ； Among them, h is the height of the pop-up, i is the current index, s top-1 The new stack top index after popping; push the current index into the stack.

6. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 1, characterized in that: The step S1 specifically includes the following steps: S11. Extract UAV flight data, including GNSS observation files, navigation data files, event log files, binary log files, and marker files; S12. Assign the geographic location information collected from the flight data extracted in step S11 to the image to achieve image georegistration; S13. Convert the LAS point cloud file with color information into a GeoTIFF image file with RGB bands; S14. Use specific ground points for image matching to achieve accurate registration of at least two UAV remote sensing images.

7. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 6, characterized in that: The step S12 specifically includes the following sub-steps: S121. Time synchronization calibration: Ensure that the image capture time and the GPS track log time are consistent in the same time system; if there is a system clock deviation, the image time needs to be compensated; S122. Trajectory interpolation algorithm: Use the calibrated image capture time as the interpolation time point, use cubic spline interpolation or linear interpolation, and output the precise position (X, Y, Z) of each image capture time, i.e., WGS84 coordinates; S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S122 to the CGCS2000 coordinate system; S124. Use the EGM2008 model to perform elevation correction and convert the GPS-measured altitude to orthometric height. The step S13 specifically includes the following sub-steps: S131. Read the LAS file and extract the point cloud data; check whether the LAS file contains color information for the three RGB channels; if there is no color information, an exception is thrown; S132. Calculate the size of the grid: Calculate the number of columns and rows of the grid based on the spatial extent of the point cloud data and the preset resolution; S133. Binned statistics: Use the binned_statistic_2d function to assign the data of each RGB channel to the calculated grid; S134.RGB Merge: Merges the three binned color channels into a three-dimensional RGB stack, where each pixel contains the color values ​​of the red, green, and blue channels, facilitating the generation of a color image. 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 and multiplying by 255; the scaled color value is updated to the RGB stack; 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; The specific sub-steps of step S14 include: 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; S142. Designate one image as the reference image and the other as the image to be registered; the reference image includes standard map coordinates or RPC information; S143. Automatically detect and match corresponding feature points in the two images using a cross-correlation method; S144. Geometric Transformation Calculation: Using a polynomial transformation model, the positions of image pixels are measured and aligned with the corresponding geographic locations, and polynomial transformation coefficients are fitted using paired points. After the fitting is completed, the obtained polynomial transformation coefficients are used to perform geometric correction on the entire image; S145. Applying cubic convolution interpolation to resample the image to be registered to adapt to the spatial resolution of the reference image; S146. Use the multiple dynamic link display feature of ENVI to compare the baseline image and the registered image to evaluate the registration accuracy.

8. The method for automatic plot segmentation of crop remote sensing images based on geographic prior information according to claim 7, characterized in that: The calculation formula for compensating the image time in step S121 is as follows: t corrected =t image +Δt; Where: t corrected is the image capture time after calibration; t image is the original capture time of the image; Δt is the system clock deviation compensation; When cubic spline interpolation is used in step S122, for the trajectory segment index j, the interpolation formula is: ； in: t Indicates the time point where 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; When linear interpolation is used, the linear interpolation formula is as follows: ； 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; The Y and Z coordinates have the same formula as the X coordinate: ； ； The step S123 uses the seven-parameter Bursa model to perform coordinate transformation: ； 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; The elevation correction formula in step S124 is: H ortho = H GPS - N geoid ； 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.

Images