Automatic crop segmentation method based on multi-source remote sensing images

Abstract

The present invention relates to the field of crop monitoring technology, and in particular to a method for automatic crop segmentation based on multi-source remote sensing images. 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; performing full-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; specifically comprising: semi-automatic segmentation of unplanted bare soil areas and areas with unknown crop growth; automatic segmentation using traditional algorithms for crop areas with known growth; and batch fully automatic segmentation processing of the images segmented in the previous step using an improved Grounded‑SAM segmentation model, outputting the segmentation results. The advantages are: achieving accurate image registration, multi-sensor compatibility, improving accuracy and efficiency through improved algorithm batch segmentation processing, and reducing manual intervention.

Description

Technical Field

[0001] The present invention relates to the technical field of crop monitoring, and in particular to an automatic crop segmentation method based on multi-source remote sensing images. Background Art

[0002] With the continuous advancement of agricultural modernization, the demand for accurate monitoring and analysis of crop growth is growing. Unmanned aerial vehicle (UAV) remote sensing, as a flexible and efficient remote sensing method, can be equipped with a variety of sensors, such as RGB cameras, hyperspectral / multispectral cameras, and thermal infrared cameras, to obtain a wealth of crop-related information. However, when using multi-source remote sensing imagery for crop segmentation, there are significant differences between the images obtained by different sensors using traditional UAV remote sensing. RGB images primarily reflect crop color information, while hyperspectral / multispectral images cover a wider spectral range and provide richer spectral features. Thermal infrared images focus on temperature information. These different image types differ significantly in spectral characteristics, resolution, and imaging principles, making it difficult to accurately segment them using the same image processing methods. For example, commonly used color thresholding-based segmentation methods for RGB images may not perform well on hyperspectral / multispectral images due to their more complex spectral information. Simple color thresholding makes it difficult to accurately distinguish crops from background. Thermal infrared image segmentation, on the other hand, requires consideration of temperature distribution characteristics, making traditional spectral feature-based segmentation algorithms often inapplicable.

[0003] In the actual process of segmenting remote sensing images, specialized remote sensing software is required to first segment the ground and vegetation in the image, and then manually crop the image to complete the plot division. This process relies on manual accuracy and is time-consuming and labor-intensive. Manual cropping is easily influenced by subjective factors, resulting in inconsistent segmentation results and reducing the efficiency and accuracy of crop monitoring. Especially when dealing with large areas of farmland and a large number of remote sensing images, the difficulty and workload of manual operation increase exponentially, seriously affecting the efficiency and scientificity of agricultural production decision-making. Existing plant phenotypic remote sensing research still relies mainly on traditional thresholding methods, exponential boosting methods (such as NDVI and ExG), or supervised learning-based classification models for the segmentation of vegetation and non-vegetation areas based on visible light or multispectral imagery. However, these methods have numerous practical challenges: first, they require a large number of manually annotated samples for model training, making them difficult to generalize to unknown areas; second, they have poor adaptability to scene variations (such as different soil backgrounds, obstructions, and water pipe interference) and limited robustness; and third, they are expensive to deploy and inefficient for large-scale remote sensing image processing tasks.

[0004] Therefore, how to develop a method that can effectively process the differences in multi-source remote sensing images and achieve automatic and accurate crop segmentation has become a technical problem that needs to be solved urgently in the current field of agricultural remote sensing. Summary of the Invention

[0005] In order to solve the above problems, the present invention provides a method for automatic crop segmentation based on multi-source remote sensing images.

[0006] The present invention aims to provide a method for automatic crop segmentation based on multi-source remote sensing images, which specifically comprises the following steps:

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

[0008] S2. Perform fully automatic or semi-automatic segmentation on the images in the same coordinate system to obtain a plot segmentation map of the multi-source crop remote sensing image. This includes semi-automatic segmentation of unplanted bare soil areas and areas with unknown crop growth; and automatic segmentation of areas with known crop growth using traditional algorithms.

[0009] S3. Use the improved Grounded-SAM segmentation model to perform batch fully automatic segmentation processing on the images segmented in step S2 and output the segmentation results.

[0010] Preferably, the step S3 specifically includes the following sub-steps:

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

[0012] S32. Using the improved Grounded-SAM segmentation model, input a prompt word or click a prompt. Grounding DINO receives the image and text prompt and outputs a bounding box in the image that is semantically related to the text to determine the location of the plant.

[0013] S33. Use the bounding box output by Grounding DINO as a hint to further refine each bounding box region and generate an accurate segmentation mask;

[0014] S34. Cropping the multi-source sensor image according to the boundary range set by the empirical model; applying the segmentation mask generated in step S33 to the cropped image for image segmentation to extract independent images;

[0015] S35. Overlay the segmentation result with the original image to generate a new image, or output a segmentation mask image; the new image or mask image includes the confidence level or visual outline of the segmented area.

[0016] Preferably, step S2 specifically includes:

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

[0018] S22. For crop areas with known growth conditions, perform batch, fully automatic segmentation processing on multi-source sensor images using traditional algorithms and output segmentation results.

[0019] The traditional algorithm is NDVI algorithm, improved EXGR algorithm or K-Means segmentation algorithm;

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

[0021] Preferably, in step S21, the boundaries of the plots are manually segmented for the unplanted bare soil areas; for the areas with unknown growing crops, the drone-photographed images of the crop plots that have grown to the maximum projected area are first selected for manual segmentation, and then the other images are batch cropped based on the range of the plot boundaries obtained by manual segmentation, and then the plot images are segmented.

[0022] Preferably, step S22 specifically includes:

[0023] S221. For RGB images or LiDAR point cloud images, use the improved EXGR segmentation algorithm for segmentation.

[0024] S222. For hyperspectral or multispectral images, use the improved EXGR segmentation algorithm or NDVI algorithm for segmentation. The NDVI calculation formula is as follows:

[0025] ;

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

[0027] S223. For thermal infrared images, use the K-Means segmentation algorithm to effectively divide different areas in the image.

[0028] The improved EXGR segmentation algorithm sets a threshold according to the EXGR value, and determines the pixel area greater than the threshold as vegetation, and the rest as non-vegetation area; the segmentation formula is as follows:

[0029] EXGR=3G-2.1RB;

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

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

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

[0033] S231. For the connected domain in the independent image segmented and extracted in step S223, 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;

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

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

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

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

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

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

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

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

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

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

[0044] S14. Use specific ground points for image matching to achieve accurate registration of two or more UAV remote sensing images.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[0061] 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. Once fitted, the resulting polynomial coefficients are used to perform geometric correction on the entire image.

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

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

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

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

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

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

[0068] ;

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

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

[0071] ;

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

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

[0074] ;

[0075] ;

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

[0077] ;

[0078] Where: [X CGCS ,YCGCS ,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;

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

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

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

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

[0083] 1. Improve segmentation accuracy: Through image registration technology, ensure that multi-source sensor images remain in the same coordinate system, improve segmentation accuracy, and avoid errors caused by sensor differences.

[0084] 2. Improve processing efficiency: Combining automatic segmentation and semi-automatic segmentation methods reduces dependence on manual operations and significantly improves the efficiency of remote sensing image processing.

[0085] 3. Applicable to multiple sensors: This method is not only applicable to the segmentation of single sensor images such as RGB, hyperspectral, thermal infrared, and lidar, but also supports the unified processing of multi-source remote sensing data, thereby improving the scope of application.

[0086] 4. Enhanced robustness: Using improved segmentation algorithms (such as the improved EXGR algorithm, K-Means algorithm, and SAM model) can adapt to different lighting conditions, environmental changes, and different crop types.

[0087] 5. Achieve efficient batch segmentation: Intelligently segment field plots based on an automatic detection model and combine it with GNSS coordinate information to enable batch and automated plot segmentation, improving large-scale data processing capabilities, reducing computational complexity, and improving computational time.

[0088] 6. Optimize crop growth monitoring: By calculating the circumscribed rectangle of the projected area, automatic segmentation of full-time remote sensing data is achieved, facilitating long-term crop growth monitoring and analysis.

[0089] 7. Reduce human intervention: Compared with the traditional method that relies on manual cropping, this method reduces human involvement, reduces human errors, and improves the consistency and stability of segmentation. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0091] Figure 2 This is a diagram showing the effect of segmentation using the improved EXGR segmentation algorithm provided by an embodiment of the present invention.

[0092] Figure 3 This is a front view effect diagram of the laser radar point cloud image segmentation provided according to an embodiment of the present invention.

[0093] Figure 4 This is a hyperspectral image segmentation effect diagram provided according to an embodiment of the present invention.

[0094] Figure 5 This is a segmentation effect diagram of a thermal infrared image provided according to an embodiment of the present invention.

[0095] Figure 6 4 is a schematic diagram of the process framework of Grounding DINO provided according to an embodiment of the present invention.

[0096] Figure 7 It is a segmentation flow chart of the Grounded-SAM segmentation model provided according to an embodiment of the present invention.

[0097] Figure 8 This is a diagram showing the effect of segmentation using an improved Grounded-SAM segmentation model provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0099] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0100] The present invention provides a method for automatic crop segmentation based on multi-source remote sensing images, which specifically includes the following steps:

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

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

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

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

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

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

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

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

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

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

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

[0112] ;

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

[0114] The specific steps include the following:

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

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

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

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

[0119] ;

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

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

[0122] ;

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

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

[0125] ;

[0126] ;

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

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

[0129] ;

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

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

[0132] Hortho = H GPS - N geoid ;

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

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

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

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

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

[0138] ;

[0139] ;

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

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

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

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

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

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

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

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

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

[0149] The linear scaling formula is as follows:

[0150] ;

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

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

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

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

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

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

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

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

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

[0160] ;

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

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

[0163] The algorithm logic is as follows:

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

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

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

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

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

[0169] ;

[0170] ;

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

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

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

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

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

[0176] The specific processing process is as follows:

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

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

[0179] ;

[0180] 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 jis 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 .

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

[0182] S1453. Generate a resampled image: Repeat steps S1451 and S1452 for each pixel in the image to generate the entire resampled image; check the quality of the resampled image to ensure that there is no obvious aliasing or blurring; if necessary, adjust the parameters of the convolution kernel to optimize the result;

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

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

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

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

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

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

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

[0190] S2. Performing fully automatic or semi-automatic segmentation on the images determined to be in the same coordinate system to obtain a plot segmentation map of the multi-source crop remote sensing image;

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

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

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

[0194] The specific steps include:

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

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

[0197] S22. For crop areas with known growth, optimize the vegetation and soil segmentation results using NDVI, improved EXGR, or K-Means segmentation algorithms to obtain a plot segmentation map for the multi-source crop remote sensing image.

[0198] For different types of images, if the background complexity is not high, the following algorithm can be used:

[0199] S221. For RGB images or LiDAR point cloud images, use the improved EXGR segmentation algorithm for segmentation. Set a threshold based on the EXGR value and define pixel areas with a value greater than the threshold as vegetation, while the rest as non-vegetation areas. The segmentation formula is as follows:

[0200] EXGR=3G-2.1RB;

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

[0202] Since the point cloud information has echo interference and missing information in the Z-axis direction, it is necessary to first perform clustering segmentation to remove noise, and then perform image subtraction on the bare soil image to subtract the land information contained in the EXGR segmented point cloud to ensure accuracy. The final front view effect of the lidar point cloud image segmentation is shown in Figure 3 .

[0203] S222. For hyperspectral or multispectral images, use the improved EXGR segmentation algorithm or NDVI algorithm for segmentation. The NDVI calculation formula is as follows:

[0204] ;

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

[0206] By setting a reasonable threshold, vegetation and soil areas can be accurately distinguished; areas with NDVI values ​​greater than 0.7 are identified as vegetation, and soil NDVI values ​​are usually 0.1~0.25. Figure 4 .

[0207] S223. For thermal infrared images, since the cell images are grayscale images, the K-Means segmentation algorithm is used to effectively divide the image into different regions. The basic process of the K-Means algorithm is as follows:

[0208] Initialization phase: Determine the constant K, which is the number of target cluster categories (pre-set according to the specific image content and segmentation requirements); randomly initialize K cluster centers, which are selected in the image grayscale value space as the initial category representatives;

[0209] Calculate the similarity between each sample and each cluster center, and divide the samples into the most similar categories; similarity is usually measured by a distance metric (such as Euclidean distance), that is, calculating the square root of the sum of the squares of the gray value differences;

[0210] Calculate the mean of all sample features divided into each category and use the mean as the new cluster center of the category;

[0211] Repeat the above iterative calculation process until the cluster center no longer changes; output the final cluster center and the category to which each sample belongs, and complete the segmentation of the thermal infrared image. Figure 5 .

[0212] S3. Use the improved Grounded-SAM segmentation model to perform batch automatic segmentation processing on the images segmented in step S2, and output the segmentation results; specifically, the following sub-steps are included:

[0213] S31. 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.

[0214] S32. Using the improved Grounded-SAM segmentation model, input a prompt word (e.g., "soybean leaf") or a point prompt (clicking one or more points in the image). Grounding DINO receives the image and text prompt and outputs a bounding box in the image that is semantically related to the text to determine the location of the plant.

[0215] Grounding DINO is an advanced zero-shot target detection model that combines the Transformer-based detector DINO with Grounded Pre-Training. Grounding DINO has zero-shot learning capabilities and can automatically identify any semantic target in the image through natural language prompts (Text Prompt) without additional training; SAM provides general, high-resolution image segmentation capabilities, which can achieve pixel-level target boundary division under the guidance of detection results. Its advantages lie in zero-shot generalization capabilities, multimodal interactivity, and efficient and accurate segmentation performance, which overcomes the traditional methods' reliance on data annotation, scene customization, and complex training. The purpose of this step is to guide the model to detect and segment target crops through natural language prompts, as well as to quickly locate the target crop area in the image. The process framework diagram of GroundingDINO is shown below. Figure 6 shown.

[0216] S33. Using the bounding box output by Grounding DINO as a hint, each bounding box region is further refined to generate an accurate segmentation mask.

[0217] S34. Crop the multi-source sensor image according to the boundary range set by the empirical model; apply the segmentation mask generated in step S33 to the cropped image to perform image segmentation and extract independent images.

[0218] S35. Overlay the segmentation result with the original image to generate a new image, or output a segmentation mask image, which includes the confidence level or visual outline of the segmented area.

[0219] The segmentation process of Grounded-SAM segmentation large model is shown in Figure 7First, an image is input into the model; the image is processed by the image encoder to generate an image embedding (a representation of image features); the image embedding is further processed by the convolutional layer to generate a mask; the prompt encoder processes the input prompt information, which can be in the form of points, boxes, or text, to guide the model for more accurate segmentation; the mask decoder receives the output of the convolutional layer and combines it with the output of the prompt encoder to generate multiple segmentation results; the model outputs multiple images (for example, images 1, 2, and 3) and corresponding valid masks, which represent the segmented areas in the image.

[0220] The effect of using improved Grounded-SAM segmentation large model segmentation is shown in Figure 8 .

[0221] Brief Description of the Invention: This invention utilizes a multimodal joint detection and segmentation system combining the Grounding DINO open-set object detector and the Segment Anything Model (SAM), specifically designed for plant phenotypic structure recognition and region extraction tasks in remote sensing imagery. Grounding DINO possesses zero-shot learning capabilities, automatically identifying any semantic target in an image through natural language prompts (Text Prompts) without the need for additional training. SAM provides general, high-resolution image segmentation capabilities, enabling pixel-level object boundary delineation guided by detection results. Its advantages lie in its zero-shot generalization, multimodal interactivity, and efficient and accurate segmentation performance, overcoming the reliance of traditional methods on data annotation, scene customization, and complex training. In particular, in large-scale remote sensing applications, it can rapidly adapt to new tasks (such as drought monitoring, land parcel segmentation, and texture segmentation), significantly reducing deployment costs.

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

[0223] 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 crop segmentation based on multi-source remote sensing images, characterized by: The specific steps include: S1. Align all images to the same coordinate system through georeferencing, point cloud file conversion, and ground point matching to achieve image registration. S2. Perform fully automatic or semi-automatic segmentation on the images in the same coordinate system to obtain a plot segmentation map for the multi-source crop remote sensing image. This includes: semi-automatic segmentation of unplanted bare soil areas and areas with unknown crop growth; and automatic segmentation of crop areas with known growth using NDVI, improved EXGR, or K-Means segmentation algorithms. This includes: S21. semi-automatically segmenting unplanted bare soil areas and areas with unknown crop growth; S22. For crop areas with known growth conditions, use the NDVI algorithm, the improved EXGR algorithm, or the K-Means segmentation algorithm to perform batch, fully automatic segmentation processing on multi-source sensor images and output the segmentation results. S23. Optimizing the segmentation results, and finally outputting an optimized cell segmentation map of the multi-source crop remote sensing image; S3. Use the improved Grounded-SAM segmentation model to perform batch automatic segmentation processing on the images segmented in step S2, and output the segmentation results; specifically, the following sub-steps are included: S31. Determine empirical model parameters: Set the empirical model based on known crop planting parameters and growth characteristics; S32. Using the improved Grounded-SAM segmentation model, input a prompt word or click a prompt; GroundingDINO receives the image and text prompt, outputs a bounding box in the image that is semantically related to the text, and determines the location of the plant; S33. Use the bounding box output by Grounding DINO as a hint to further refine each bounding box region and generate an accurate segmentation mask; S34. Cropping the multi-source sensor image according to the boundary range set by the empirical model; applying the segmentation mask generated in step S33 to the cropped image for image segmentation to extract independent images; S35. Overlay the segmentation result with the original image to generate a new image, or output a segmentation mask image; the new image or mask image includes the confidence level or visual outline of the segmented area.

2. The method for automatic crop segmentation based on multi-source remote sensing images 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 crop segmentation based on multi-source remote sensing images according to claim 2, characterized in that: The step S22 specifically includes: S221. For RGB images or LiDAR point cloud images, use the improved EXGR segmentation algorithm for segmentation. S222. For hyperspectral or multispectral images, use the improved EXGR segmentation algorithm or NDVI algorithm for segmentation. The NDVI calculation formula is as follows: ； 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; S223. For thermal infrared images, use the K-Means segmentation algorithm to perform image segmentation and effectively divide different areas in the image. The improved EXGR segmentation algorithm sets a threshold according to the EXGR value, and determines the pixel area greater than the threshold as vegetation, and the rest as non-vegetation area; the segmentation formula is as follows: EXGR=3G-2.1RB; Among them, G, R, and B represent the intensity values ​​of the green, red, and blue channels of each pixel in the image, respectively.

4. The method for automatic crop segmentation based on multi-source remote sensing images according to claim 1, characterized in that: The optimization process 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.

5. The method for automatic crop segmentation based on multi-source remote sensing images according to claim 4, characterized in that: The step S23 specifically includes the following sub-steps: S231. For the connected domain in the segmented independent image, 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.

6. The method for automatic crop segmentation based on multi-source remote sensing images 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 two or more UAV remote sensing images.

7. The method for automatic crop segmentation based on multi-source remote sensing images 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 is consistent with the GPS track log 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 S202 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 their corresponding geographic locations. The polynomial transformation coefficients are then fitted using the paired points. Once fitted, the resulting polynomial 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 crop segmentation based on multi-source remote sensing images 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 at which interpolation is required; j Indicates the trajectory segment index; t j Indicates the trajectory segment index j The starting time point; S j ( t ) represents the trajectory segment index j In time t The position function value of ; Indicates the trajectory segment index j In time t The second derivative of the position function; Indicates the trajectory segment index j +1 in time t j The second derivative of the position function; coefficient a j 、b j 、c j d j Solve by boundary conditions; 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