Geographic prior information-based crop remote sensing image cell automatic segmentation method

Abstract

The invention relates to the technical field of crop monitoring, in particular to an automatic cell segmentation method for a crop remote sensing image based on geographic prior information. Comprising the steps that all images are determined to be in the same coordinate system through geographical registration, point cloud file conversion and ground specific point matching, and image registration is achieved; performing automatic segmentation or semi-automatic segmentation on the images determined in the same coordinate system to obtain a cell segmentation map of the multi-source crop remote sensing image; performing semi-automatic segmentation on an unplanted bare soil region and a region of crops with unknown growth vigor; aiming at a crop area with known growth vigor, carrying out batch full-automatic segmentation processing on the multi-source sensor image by adopting an improved Ground-SAM segmentation large model, and outputting a segmentation result; and carrying out optimization processing on the segmentation result, and finally outputting an optimized cell segmentation map of the multi-source crop remote sensing image. The method has the advantages of accurate registration, efficient segmentation, multi-sensor compatibility and automatic processing, and remote sensing monitoring is improved.

Description

Technical Field

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

[0002] With the continuous development of agricultural informatization and precision agriculture, remote sensing technology has been increasingly widely used in the agricultural field, especially in field monitoring and analysis. There are significant differences between images obtained by different sensors in traditional UAV remote sensing, and it is impossible to use the same image processing method for accurate segmentation. For example, in the patent application with the publication number CN118334522A, publication date July 12, 2024, and invention title "Method and device for automatically extracting plots of high-throughput field crops based on low-altitude remote sensing", an image processing method is used to segment the land, and then the picture is cropped based on the estimated value, with certain errors. Geographical location information is not considered, and it is impossible to batch segment all the shots in one experiment using the same coordinates. Another example is the patent application with the publication number CN116469019A, publication date July 21, 2023, and invention title "Method for automatically segmenting field plots based on fully convolutional neural network", which also does not consider the positioning of geographical location information and only uses the network for batch segmentation. For large pictures in the field, a huge amount of computing is usually required for segmenting each individual plot.

[0003] The common problems in the prior art include the following aspects: (1) Sensor differences: Images obtained by different sensors (such as RGB cameras, thermal infrared cameras, hyperspectral cameras, lidar, etc.) have significant differences in spectral characteristics, resolution, and data format, and it is difficult for traditional methods to be applicable to the segmentation of images from multiple sensors simultaneously. (2) Inefficiency of manual segmentation: In agricultural remote sensing detection, the planning and segmentation of field plots are essential steps. Currently, many studies and practical applications still rely on manual segmentation or use software such as arcgis for plot planning and then batch segmentation. However, due to the mismatch of geographical location information, the manual matching and segmentation process is time-consuming and laborious, and prone to errors. (3) Processing of images with different time series: In long-term monitoring, it is necessary to segment and analyze remote sensing images with different time series. Existing methods cannot effectively handle the differences between images with different time series and are difficult to achieve automation and batch processing. Summary of the Invention

[0004] 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 object of the present invention is to provide a method for automatically segmenting plots of crop remote sensing images based on geographic prior information, which specifically includes the following steps: S1. Through georegistration, point cloud file conversion, and matching of specific ground points, all images are determined in the same coordinate system to achieve image registration. S2. Automatically or semi-automatically segment the images determined in the same coordinate system to obtain the sub-region segmentation map of the multi-source crop remote sensing images. Specifically, it includes the following steps: S21. Semi-automatically segment the unplanted bare soil areas and areas with unknown crop growth. S22. For the areas with known crop growth, use the improved Grounded-SAM segmentation model to perform batch full-automatic segmentation processing on the multi-source sensor images and output the segmentation results. S23. Optimize the segmentation results and finally output the optimized sub-region segmentation map of the multi-source crop remote sensing images.

[0006] Preferably, in step S21, for the unplanted bare soil areas, manually segment the sub-region boundaries; for the areas with unknown crop growth, first select the UAV images of the crop sub-regions that have grown to the maximum projected area for manual segmentation, and then, based on the sub-region boundary range obtained through manual segmentation, batch crop the other images and then perform sub-region image segmentation.

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

[0008] Preferably, the prompt words in step S222 are plants. The images output in step S225 include the confidence levels or visualization contours of the segmented regions.

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

[0010] Preferably, step S23 specifically includes the following sub-steps: S231. For the connected regions in the independent image segmented and extracted in step S224, traverse each pixel row by row, calculate the number of consecutive 1s upward at each position to form a height matrix height, and determine the shape and size of the connected regions of each plot; S232. For each row of the height array, use the monotonic stack algorithm to solve the maximum rectangle area of the current row; the specific process is as follows: S2321. Initialize the stack: The stack is used to save the indexes of the columns with increasing heights. The initial state of the stack is empty, and the heights corresponding to the bottom to the top of the stack increase in turn; S2322. Traverse each element of the current row: When the stack is not empty and the current height is less than the height of the top of the stack, pop the top element top of the stack and calculate the rectangle area with the height height[top]; S2323. Process the remaining elements in the stack: After the traversal is completed, perform step S2322 on the remaining elements in the stack to calculate the possible areas; S2324. Update the global maximum area: After processing each row, record the maximum rectangle area of the current row; traverse all rows and take the maximum value among all rows as the global maximum area; S233. For each connected region of the segmented independent plot image, calculate its maximum circumscribed rectangle and output the maximum growth projection area of each plot.

[0011] Preferably, the height matrix formula is as follows: ; In the formula, i represents the row index in the image, and j represents the column index in the image; The rectangle area formula is as follows: ; Where h is the popped height, i is the current index, s top-1 is the new top index of the stack after popping; push the current index onto the stack.

[0012] Preferably, step S1 specifically includes the following steps: S11. Extract the UAV flight data, and the flight data includes GNSS observation files, navigation data files, event record files, binary log files, and marker files; S12. Assign the geographical location information collected from the flight data extracted in step S11 to the image to achieve the georegistration of the image; 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 the precise registration of at least two UAV remote sensing images.

[0013] Preferably, step S12 specifically includes the following sub-steps: S121. Time synchronization calibration: Ensure that the shooting time of the image is consistent with the time in the GPS trajectory log 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 shooting time as the interpolation time point, and adopt cubic spline interpolation or linear interpolation to output the precise position (X, Y, Z) at each image shooting moment, that is, the WGS84 coordinates; S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S122 into the CGCS2000 coordinate system; S124. Use the EGM2008 model for elevation correction to convert the altitude measured by GPS into orthometric height; The said 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 in the RGB three channels; if there is no color information, it is thrown as an exception; S132. Calculate the size of the grid: Calculate the number of columns and rows of the grid according to the spatial range of the point cloud data and the preset resolution; S133. Binned statistics: Use the binned_statistic_2d function to distribute the data of each RGB channel into the calculated grid; S134. RGB merging: Merge the three color channels after binning 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 it is necessary to scale the color values from 16 bits to 8 bits according to actual needs; if scale color = True, then scale the color values from 16 bits to 8 bits by dividing the value of each color channel by 65535 and multiplying by 255; update the scaled color values to the RGB stack; S136. Generate a GeoTIFF file: Determine the metadata of the GeoTIFF file; Open the GeoTIFF file using the rasterio.open function and pass in the metadata; Write the three color channels of the RGB stack to the corresponding bands of the GeoTIFF file. The specific sub-steps of step S14 include: S141. Identify and select specific fixed ground points in the image, re-assign fixed coordinate points to the selected specific ground points, and use the fixed coordinate points as the 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. Use the cross-correlation method to automatically detect and match the corresponding feature points in the two images. S144. Geometric transformation calculation: Adopt a polynomial transformation model, measure the positions of the image pixels, register them with the corresponding geographical positions, and fit the polynomial transformation coefficients using the paired points; After the fitting is completed, use the obtained polynomial transformation coefficients to perform geometric correction on the entire image. S145. Apply the cubic convolution interpolation method to resample the image to be registered to adapt to the spatial resolution of the reference image. S146. Use the multiple dynamic link display function of ENVI to compare the reference image and the registered image and evaluate the registration accuracy.

[0014] Preferably, the calculation formula for compensating the image time in step S121 is as follows: t corrected = t image + Δt; Where: t corrected is the calibrated image capture time; t image is the original image capture time; Δt is the system clock deviation compensation amount. When using cubic spline interpolation in step S122, for the trajectory segment index j, the interpolation formula is: ; Where: t represents the time point where interpolation is required; j represents the trajectory segment index; t j represents the starting time point of the trajectory segment index j ; S j ( t ) represents the position function value of the trajectory segment index j at time t ; represents the trajectory segment indexj 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 = HGPS - N geoid ; In the formula: H ortho represents the orthometric height; H GPS represents the altitude measured by GPS; N geoid represents the geoid height anomaly calculated by model interpolation.

[0015] Compared with the prior art, the present invention can achieve the following beneficial effects: 1. Precise registration and batch segmentation: By fusing the flight parameters of the unmanned aerial vehicle and geospatial information, the registration of images taken at different times and by different sensors is realized, so that all remote sensing images are kept in the same coordinate system, thereby improving the segmentation accuracy and supporting batch segmentation; 2. Reducing the amount of computation and improving the calculation efficiency: First, use geographic information for cropping, and then perform cell segmentation, avoiding complex operations on the entire large image and reducing the calculation cost; 3. Applicable to a variety of remote sensing sensors: Support various data sources such as RGB, hyperspectral, multispectral, lidar, thermal infrared, etc., and solve the problem that data from different sensors cannot be uniformly processed; 4. Compatible with automatic and semi-automatic segmentation: For crops with known growth conditions, use an empirical model for fully automatic segmentation; for crops with unknown growth conditions, first determine the maximum growth projection area and then perform batch cropping to improve flexibility; 5. Improved segmentation method: Adopt the Grounded-SAM model combined with DINO-X object detection to achieve accurate target recognition and improve the segmentation accuracy; 6. Efficient processing of multi-temporal images: Through the calculation of the maximum growth projection area, batch segmentation of full-temporal remote sensing images is realized, reducing manual intervention and improving the automation degree of remote sensing monitoring. Description of the Drawings

[0016] Figure 1 is a visual display effect diagram of two flight files for registration in a method for automatic cell segmentation of crop remote sensing images based on geographic prior information provided by an embodiment of the present invention.

[0017] Figure 2 is the image segmentation effect after automatic cell segmentation provided by an embodiment of the present invention.

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

[0019] Figure 4is the automatic segmentation effect of the visible light image provided by the embodiment of the present invention. Detailed implementation mode

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

[0021] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below 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 to the present invention.

[0022] The present invention provides a method for automatically segmenting plots of crop remote sensing images based on geographical prior information, which specifically includes the following steps: S1. Through georegistration, point cloud file conversion, and ground specific point matching, all images are determined in the same coordinate system to achieve image registration; specifically including the following steps: S11. Extract the flight data of the unmanned aerial vehicle, and the flight data includes GNSS observation files (.obs), navigation data files (NAV files), event record files (EVENT files), binary log files (BIN files), and marker files (MRK files); specifically as follows: GNSS observation file (.obs): A standard data file that mainly stores the raw measurement data from the global navigation satellite system; 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; Navigation data file (NAV file): Stores information about the flight path, GPS coordinates, altitude, and other relevant data, which can be used for geographical reference of images, flight analysis, or troubleshooting; Event record file (EVENT file): Records various events that occur during the flight, such as when the unmanned aerial vehicle takes off, when an image is captured, or when a certain waypoint is reached; can help understand the progress of the flight and identify any problems that may occur during the mission; Binary log file (BIN file): Contains detailed information about the flight of the unmanned aerial vehicle, including telemetry data, sensor readings, and system information, which can be used for in-depth analysis and troubleshooting of the performance of the unmanned aerial vehicle; Marker file (MRK file): Records the relevant information of the photo-taking points, where: Column 1: Serial number of the photo-taking point, i.e., the serial number of the photo record information stored in this folder; Column 2: Seconds within GPS week, indicating the photo-taking moment, expressed in the form of seconds within GPS week (TOW); Column 3: GPS week, indicating the GPS week of the photo-taking moment; Column 4: Compensation value in the north direction, in millimeters, with north being positive; Column 5: Compensation value in the east direction, in millimeters, with east being positive; Column 6: Compensation value in the elevation direction, in millimeters, with down being positive; Column 7: Longitude after compensation; Column 8: Latitude after compensation; Column 9: Ellipsoidal height; Column 10: Positioning standard deviation (north); Column 11: Positioning standard deviation (east); Column 12: Positioning standard deviation (elevation); Column 13: Positioning status.

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

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

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

[0026] S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S122 to the CGCS2000 coordinate system (meeting Chinese requirements); use the seven-parameter Bursa model for coordinate transformation: ; In the formula: [X CGCS , Y CGCS , Z CGCS represents the coordinates in the 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, representing the scale difference between the two coordinate systems; represents the rotation matrix, , , are the rotation parameters, representing the rotation angles around the X, Y, and Z axes respectively.

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

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

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

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

[0031] S134. RGB Merging: 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 red, green, and blue channels, facilitating the generation of color images; The specific operations are as follows: Initialize the RGB stack: Create a three-dimensional array (rgb stack ), with a shape of (rows, cols, 3), where rows and cols are the number of rows and columns of the grid respectively, and 3 represents the three RGB color channels. Fill the RGB stack: Fill the binned red, green, and blue channel data into the first three positions of the third dimension of rgb stack , for example, fill the red channel data into rgb stack [:,:,0], the green channel data into rgb stack [:,:,1], and the blue channel data into rgb stack [:,:,2]. Ensure that the data types of all channels are consistent before merging, usually unsigned 8-bit integers (uint8), representing the color value range from 0 to 255.

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

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

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

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

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

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

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

[0039] The convolution kernel calculates weights based on the pixel intervals from the position of the target pixel; the Catmull-Rom interpolation kernel is adopted, which is a smooth cubic function that can effectively reduce the jagged effect and blurring phenomenon during the interpolation process; pixels closer to the target pixel will obtain higher weights, while pixels farther away will obtain lower weights; S1453. Generate the 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 there is no obvious jagged effect or blurring phenomenon; If necessary, adjust the parameters of the convolution kernel to optimize the result; S1454. Output the final image: Output the resampled image as the final result for subsequent image analysis or applications; The bicubic convolution interpolation method uses more pixel value information, so it can usually produce higher-quality results when enlarging images; However, its computational complexity is also relatively high.

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

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

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

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

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

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

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

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

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

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

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

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

[0052] For different types of images, if the background complexity is not high, etc., not limited to the improved Grounded-SAM segmentation large model, the NDVI, EXGR or Otsu algorithm can also be used to optimize the vegetation and soil segmentation results to obtain the small cell segmentation map of the multi-source crop remote sensing image; specifically as follows: For RGB images, hyperspectral images, multispectral images or lidar images, use the improved EXGR segmentation algorithm for segmentation; Set a threshold according to the EXGR value, and determine the pixel area greater than the threshold as vegetation, and the rest as non-vegetation areas; The segmentation formula is as follows: EXGR = 3G - 2.1R - B; Among them, G, R, and B respectively represent the intensity values of the green, red, and blue channels of each pixel in the image; For hyperspectral images or multispectral images, the NDVI algorithm can also be used for segmentation; The NDVI calculation formula is as follows: ; In the formula, NDVI represents the Normalized Difference Vegetation Index, and NIR and Red respectively represent the intensity values of the near-infrared and red light channels of each pixel in the image; By setting a reasonable threshold, the vegetation and soil areas can be accurately distinguished; the area with an NDVI value greater than 0.7 is determined as vegetation, and the NDVI value of the soil is usually 0.1 - 0.25.

[0053] For the thermal infrared image, since the thermal infrared image usually reflects the temperature in the imaging area, the Otsu algorithm based on threshold is used to segment the thermal infrared image; the temperature is divided into two categories, the higher temperature is for vegetation and the lower temperature is for soil; specifically including: (1) Calculate the gray-level histogram p i ; assume the image has a total of L gray levels (such as 0 - 255), and p i represents the frequency of the gray level i appearing in the image, satisfying ; (2) Calculate the class probabilities: select a threshold t, and divide the pixels into two categories: Background class C0: gray levels 0 - t, with a probability of ; Foreground class C1: gray levels t + 1 - L - 1, with a probability of ; (3) Calculate the class means and the overall mean: Background mean: ; Foreground mean: ; In the formula: μ 0( t ) and μ 1( t ) respectively represent the means of the background and foreground at the threshold t ; ω 0 and ω 1 are respectively the weights of the background and foreground, that is, the probabilities of background-class pixels; i is the index of the gray level; Calculate the overall mean: ; (4) Calculate the between-class variance: ; It can be simplified to: ; (5) Optimize the threshold: traverse all possible thresholds t, and calculate σ 2 (t); select the 2 that makes σ (t) the largest as the segmentation threshold; ; (6) According to Generate a binary image, determine the areas with high temperature as vegetation, and the areas with low temperature as soil: ; Compare the results of automatic segmentation and semi-automatic segmentation to verify the accuracy and consistency of the segmentation.

[0054] The key technical points of the present invention include the following aspects: 1. Design of the cropping method: By batch cropping the registered images, multi-source remote sensing image segmentation can be performed by manually or automatically setting the coordinate system once; 2. Automatic segmentation based on deep learning: Adopt the Grounded-SAM large model, input crop keywords for automatic object detection and segmentation; 3. Batch segmentation based on projected area: Through the maximum projected area calculation method, batch plot segmentation is performed on multi-temporal images of the same crop; 4. Optimize the calculation efficiency: Use the monotonic stack algorithm to solve the maximum rectangle area and improve the calculation speed of plot segmentation; 5. Improve the traditional algorithm: For RGB and hyperspectral / multispectral images: Use the EXGR (3G - 2.1R - B) algorithm and the improved NDVI for segmentation; For thermal infrared images: Use the Otsu adaptive threshold method for vegetation and soil classification.

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

[0056] Figure 3 The semi-automatic segmentation effect is shown, that is, after image correction, first obtain the UTM coordinates of the vertices of each plot, and use the UTM coordinates 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 is shown, including the automatic recognition of plots in visible light images, as well as the complete process and results of generating a mask to extract vegetation images.

[0057] It should be understood that various forms of processes shown above can be reordered, steps can be added or deleted. For example, the steps recorded in the disclosure of the present invention can be executed in parallel, sequentially, or in a different order, as long as the desired results of the technical solution disclosed in the present invention can be achieved, and no limitations are imposed herein.

[0058] The above specific embodiments do not constitute a limitation to the protection scope of the present invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions and improvements made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for automatically segmenting plots in remote sensing images of crops based on geographical prior information, characterized in that: Specifically, it includes the following steps: S1. Through georegistration, point cloud file conversion, and matching of specific ground points, all images are determined in the same coordinate system to achieve image registration; S2. Automatically or semi-automatically segment the images determined in the same coordinate system to obtain the sub-region segmentation map of the multi-source crop remote sensing images; specifically, it includes the following steps: S21. Semi-automatically segment the unplanted bare soil areas and areas with unknown crop growth; S22. For the areas of known crop growth, use the improved Grounded-SAM segmentation large model to perform batch full-automatic segmentation processing on the multi-source sensor images and output the segmentation results; S23. Optimize the segmentation results and finally output the optimized sub-region segmentation map of the multi-source crop remote sensing images.

2. The automatic sub - area segmentation method of crop remote sensing images based on geographical prior information according to claim 1, characterized in that: In the step S21, for the unplanted bare soil areas, manually segment the sub-region boundaries; for the areas with unknown crop growth, first select the UAV captured images of the crop sub-regions that have grown to the maximum projected area for manual segmentation, and then, based on the sub-region boundary range obtained through manual segmentation, batch crop the other images and then perform sub-region image segmentation.

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

4. The automatic sub - area segmentation method of crop remote sensing images based on geographic prior information according to claim 3, characterized in that: The prompt word in the step S222 is "plant"; The image output in the step S225 includes the confidence level or visualization contour of the segmented region.

5. The automatic sub - area segmentation method of crop remote sensing images based on geographical prior information according to claim 3, characterized in that: The optimization processing in the step S23 specifically includes: constructing a height matrix to determine the shape of the vegetation area, using the monotonic stack algorithm to calculate the maximum rectangle area row by row, and calculating the maximum growth projected area of each sub-region based on the empirical model, and finally outputting the optimized sub-region segmentation map of the multi-source crop remote sensing images.

6. The automatic sub - area segmentation method of crop remote sensing images based on geographic prior information according to claim 5, characterized in that: The step S23 specifically includes the following sub-steps: S231. For the connected components in the independent images segmented and extracted in step S224, traverse each pixel row by row, calculate the number of consecutive 1s upward at each position to form the height matrix height, and determine the shape and size of the connected components of each sub-region; S232. For each row of the height array, use the monotonic stack algorithm to solve the maximum rectangle area of the current row; the specific process is as follows: S2321. Initialize the stack: The stack is used to save the indices 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. S2322. Traverse each element in the current row: When the stack is not empty and the current height is less than the height at the top of the stack, pop the top element top of the stack, and calculate the area of the rectangle with height height[top]. S2323. Process the remaining elements in the stack: After the traversal, perform step S2322 on the remaining elements in the stack to calculate the possible areas. S2324. Update the global maximum area: After processing each row, record the maximum rectangle area of that row; traverse all rows and take the maximum value among all rows as the global maximum area. S233. For each separated independent small cell image connected component, calculate its maximum circumscribed rectangle, and output the maximum growth projection area of each small cell.

7. The automatic sub - area segmentation method of crop remote sensing images based on geographical prior information according to claim 6, characterized in that: The height matrix formula is as follows: ； In the formula, i represents the row index in the image, and j represents the column index in the image. The rectangle area formula is as follows: ； where h is the popped height, i is the current index, and s top-1 is the new stack top index after popping; push the current index onto the stack.

8. The automatic sub - area segmentation method of crop remote sensing images based on geographical prior information according to claim 1, characterized in that: The specific steps of step S1 are as follows: S11. Extract the UAV flight data, and the flight data includes GNSS observation files, navigation data files, event record files, binary log files, and marker files. S12. Assign the geographic location information collected in the flight data extracted in step S11 to the image to achieve the georegistration of the image. 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 the precise registration of at least two UAV remote sensing images.

9. The method for automatically segmenting plots of crop remote sensing images based on geographic prior information according to claim 8, wherein: The specific steps of step S12 are as follows: S121. Time synchronization calibration: Ensure that the shooting time of the image is consistent with the time in the GPS trajectory log 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 shooting time as the interpolation time point, and adopt cubic spline interpolation or linear interpolation to output the precise position (X, Y, Z), that is, the WGS84 coordinates, at each image shooting moment. S123. Convert the WGS84 coordinates in the GPS trajectory obtained in step S122 into the CGCS2000 coordinate system. S124. Use the EGM2008 model for elevation correction to convert the altitude measured by GPS into orthometric height. The specific steps of step S13 are as follows: S131. Read the LAS file and extract the point cloud data; check whether the LAS file contains color information for the RGB three channels; if there is no color information, it is thrown as an exception. S132. Calculate the size of the grid: According to the spatial range of the point cloud data and the preset resolution, calculate the number of columns and rows of the grid. S133. Binned statistics: Use the binned_statistic_2d function to allocate the data of each RGB channel into the calculated grid. S134. RGB Merging: Merge the three color channels after binning into a three-dimensional RGB stack, where each pixel contains the color values of the red, green, and blue channels, facilitating the generation of color images; S135. Color Scaling: Determine whether it is necessary to scale the color value from 16 bits to 8 bits according to actual needs; if scale color = True, scale the color value from 16 bits to 8 bits by dividing the value of each color channel by 65535 and multiplying by 255; update the scaled color value to the RGB stack; S136. Generating 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 specific fixed ground points in the image, reassign fixed coordinate points to the selected specific ground points, and use the fixed coordinate points as the reference points for image registration; 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. Use the cross-correlation method to automatically detect and match the corresponding feature points in the two images; S144. Geometric transformation calculation: Adopt a polynomial transformation model, measure the positions of the image pixels, register them with the corresponding geographical positions, and use the paired points to fit the polynomial transformation coefficients; After the fitting is completed, use the obtained polynomial transformation coefficients to perform geometric correction on the entire image; S145. Apply 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 function of ENVI to compare the reference image and the registered image and evaluate the registration accuracy.

10. The automatic sub - area segmentation method of crop remote sensing images based on geographical prior information according to claim 9, 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 calibrated image capture time; t image is the original capture time of the image; Δt is the system clock deviation compensation amount; When using cubic spline interpolation in step S122, for the trajectory segment index j, the interpolation formula is: ； Wherein: t represents the time point at which interpolation is required; j represents the trajectory segment index; t j represents the trajectory segment index j of the starting time point; S j ( t ) represents the trajectory segment index j at time t of the position function value; represents the trajectory segment index j at time t of the second derivative of the position function; represents the trajectory segment index j +1 at time t j of the second derivative of the position function; coefficient a j 、b j 、c j 、d j are solved by boundary conditions; When using linear interpolation, the linear interpolation formula is as follows: ； In the formula: X ( t ) represents the X - coordinate value at time t; X {k} represents at time t {k} the X - coordinate value; X {k+1} represents at time t {k+1} the X - coordinate value; t represents the time point at which the X - coordinate value needs to be estimated; t {k} and t {k+1} represent two time points with known X - coordinate values; The Y and Z coordinate formulas are the same as the X coordinate formula: ； ； Step S123 uses the seven-parameter Bursa model for 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, representing the scale difference between the two coordinate systems; represents the rotation matrix, , , are the rotation parameters, representing the rotation angles around the X, Y, and Z axes respectively; The elevation correction formula in step S124 is: H ortho = H GPS - N geoid ； In the formula: H ortho represents the orthometric height; H GPS represents the altitude measured by GPS; N geoid represents the geoid height anomaly calculated by model interpolation.

Images