Crop planting area intelligent extraction method and system based on multispectral remote sensing

Abstract

The invention provides a crop planting area intelligent extraction method and system based on multispectral remote sensing, and relates to the field of remote sensing image processing, and the method comprises the steps: obtaining and correcting a multiband remote sensing image; calculating a vegetation index and constructing a crop growth characterization index to extract canopy features; obtaining a multi-period feature map, calculating a spatial distribution entropy, and establishing an evaluation function to determine a time sequence fusion weight; a spatial constraint function is established based on the spectral distance, and a segmentation criterion is constructed by combining the spatial constraint function with time sequence features for classification iteration. According to the method, the accuracy and the discrimination degree of crop planting area extraction are improved, and crop identification requirements in a complex agricultural environment are met.

Description

Technical Field

[0001] This invention relates to remote sensing image processing technology, and more particularly to a method and system for intelligent extraction of crop planting areas based on multispectral remote sensing. Background Technology

[0002] Multispectral remote sensing technology is an important tool for extracting crop planting areas in modern agricultural management. Current remote sensing identification methods mainly rely on vegetation indices and spectral features for crop region division. However, due to the dynamic changes in spectral features during crop growth, remote sensing features from a single period are insufficient to accurately represent the crop's growth status. Existing methods typically use fixed weights or simple arithmetic averages when fusing features from multiple periods, failing to fully utilize the feature contributions of different growth stages, resulting in feature fusion results that do not accurately reflect the crop's growth process. Furthermore, spatial correlations between pixels are often ignored during region segmentation, making classification results susceptible to noise interference and resulting in inaccurate segmentation boundaries. Therefore, a smart identification method is needed that can effectively utilize multi-period remote sensing data to accurately extract crop planting areas. Summary of the Invention

[0003] This invention provides a method and system for intelligent extraction of crop planting areas based on multispectral remote sensing, which can solve the problems in the prior art.

[0004] A first aspect of this invention provides a method for intelligent extraction of crop planting areas based on multispectral remote sensing, comprising:

[0005] Acquire remote sensing image data in multiple bands and perform atmospheric and geometric corrections to obtain corrected remote sensing image data;

[0006] The vegetation coverage index and vegetation growth vitality index are calculated based on the corrected remote sensing image data. The band combination coefficient is calculated based on the vegetation coverage index and vegetation growth vitality index. A new crop growth characterization index is constructed using the band combination coefficient. The crop canopy structure features are extracted based on the crop growth characterization index to generate a crop feature map.

[0007] The crop feature maps of multiple periods within the crop growth cycle are obtained. The spatial distribution entropy of pixel values ​​in the crop feature maps of each period is calculated. A growth status evaluation function is established. The feature contribution of each period is calculated based on the growth status evaluation function as the temporal fusion weight. The crop feature maps of multiple periods are weighted and combined to obtain temporal feature data.

[0008] The spectral distance between pixels is calculated from the temporal feature data. A spatial constraint function is established based on the spectral distance. The spatial constraint function is combined with the temporal feature data to construct a region segmentation criterion. The pixels are then classified and iteratively processed to output the location information of the crop planting area and the crop type information.

[0009] In one alternative embodiment,

[0010] Multiple bands of remote sensing image data were acquired and atmospheric and geometric corrections were performed to obtain the corrected remote sensing image data, including:

[0011] Remote sensing image data in the visible and infrared bands is acquired using a multispectral sensor. The remote sensing image data includes spectral and spatial information of ground object reflection.

[0012] Atmospheric parameters of the region corresponding to the remote sensing image data are collected, atmospheric influence features are extracted based on the atmospheric parameters, and atmospheric correction is performed on the remote sensing image data based on the atmospheric influence features to obtain spectral reflectance data after eliminating atmospheric influence.

[0013] Obtain uniformly distributed ground control point information within the geographic area corresponding to the remote sensing image data. The ground control point information includes the geographic coordinates of the control points and the corresponding image features. Construct a coordinate mapping relationship based on the ground control point information. Use the coordinate mapping relationship to perform geometric correction on the spectral reflectance data to generate initial corrected data with geographic coordinates.

[0014] The spectral response change between adjacent bands and the spatial position deviation between adjacent pixels are calculated in the initial correction data. Abnormal regions are identified based on the spectral response change and spatial position deviation. The abnormal regions are then corrected to obtain the corrected remote sensing image data.

[0015] In one alternative embodiment,

[0016] The vegetation cover index and vegetation growth vitality index are calculated based on the corrected remote sensing image data. The band combination coefficients are then calculated based on these indices, including:

[0017] Collect spectral data of farmland soil in the crop planting area, extract soil background features based on the spectral data, and establish a feature library containing the soil background features;

[0018] Each pixel in the corrected remote sensing image data is matched with the soil background features in the feature library. The corresponding soil type is identified based on the matching results. The spectral response of the pixel is corrected based on the soil type to obtain spectral data that eliminates the influence of soil background and thus obtain vegetation spectral data.

[0019] Red edge location features and chlorophyll absorption features are extracted from the vegetation spectral data. The vegetation coverage index is calculated based on the red edge location features and chlorophyll absorption features. At the same time, the vegetation growth vitality index is calculated using the first derivative features of the vegetation spectral data.

[0020] The vegetation coverage index and vegetation growth vitality index are decomposed into multiple scales respectively. By calculating the complementary information of the vegetation coverage index and vegetation growth vitality index at each scale, the information weight of each scale is determined. The information weight is then normalized to obtain the initial band combination coefficient.

[0021] Canopy spectral data of crops at different growth stages are collected, and the variation of the canopy spectral data with the growth stage is analyzed. The initial band combination coefficients are dynamically corrected according to the variation pattern to obtain the band combination coefficients.

[0022] In one alternative embodiment,

[0023] New crop growth characterization indices are constructed using band combination coefficients. Crop canopy structure features are then extracted based on these indices, generating crop feature maps, including:

[0024] The spatial distribution gradients of vegetation coverage index and growth vitality index are calculated separately. A regional segmentation threshold is established, and the crop planting area is segmented using the regional segmentation threshold to obtain multiple sub-regions.

[0025] Calculate the statistical characteristics of vegetation coverage index and growth vitality index in each sub-region to obtain the spectral distribution characteristics of the sub-region. Based on the spectral distribution characteristics, perform regional correction on the band combination coefficient to generate regional band combination coefficients.

[0026] The vegetation coverage index and growth vitality index are weighted and fused based on the regional band combination coefficients to construct a crop growth characterization index.

[0027] Wavelet decomposition is performed on the crop growth characterization index to obtain multi-scale wavelet coefficients. The energy distribution of wavelet coefficients at each scale is calculated, and an inter-scale correlation metric is constructed. The combined weights of the coefficients at each scale are determined based on the inter-scale correlation metric. The multi-scale wavelet coefficients are reconstructed using the combined weights to obtain the crop canopy structure characteristics.

[0028] The vertical height distribution and horizontal density distribution are extracted based on the characteristics of the crop canopy structure.

[0029] The height and density distributions of multiple time phases are obtained, the distribution changes between time phases are calculated, and spatiotemporal interpolation is performed based on the distribution changes to generate a crop feature map.

[0030] In one alternative embodiment,

[0031] The spatial distribution entropy of pixel values ​​in crop feature maps at different periods is calculated, a growth state evaluation function is established, and the feature contribution of each period is calculated as the temporal fusion weight based on the growth state evaluation function. The crop feature maps from multiple periods are then weighted and combined to obtain the temporal feature data, including:

[0032] Perform multi-scale decomposition on the crop feature map of each period to obtain feature components at different scales, calculate the pixel value distribution probability of each scale feature component, and construct the spatial distribution entropy based on the pixel value distribution probability.

[0033] Analyze the changing trends of spatial distribution entropy at different scales in adjacent periods, extract the principal direction vectors of spatial distribution entropy at different scales, calculate the orthogonal projection metric between scales based on the principal direction vectors, and fuse the spatial distribution entropy at each scale according to the orthogonal projection metric to obtain the quantitative value of the growth state of the period.

[0034] Calculate the spatial autocorrelation coefficient of local regions in the crop feature map, construct a growth status evaluation function based on the spatial autocorrelation coefficient, input the quantitative value of the growth status of the period into the growth status evaluation function, and calculate the feature contribution of the crop feature map of each period.

[0035] The gradient direction field of crop feature maps from adjacent periods is extracted, and the gradient direction field is used to construct a temporal variation tensor to dynamically adjust the feature contribution and generate temporal fusion weights.

[0036] The crop feature maps from multiple periods are weighted and combined based on the temporal fusion weight to obtain a fused feature map. The anisotropy index of the fused feature map is calculated to construct a regional smoothing constraint. The fused feature map is then subjected to temporal analysis based on the regional smoothing constraint to generate temporal feature data.

[0037] In one alternative embodiment,

[0038] Analyze the changing trends of spatial distribution entropy at various scales in adjacent periods, extract the principal direction vectors of spatial distribution entropy at different scales, calculate the orthogonal projection metric between scales based on the principal direction vectors, and fuse the spatial distribution entropy at each scale according to the orthogonal projection metric to obtain the quantitative value of the growth state of each period, including:

[0039] Multi-scale decomposition of spatial distribution entropy in adjacent periods yields spatial distribution entropy at different scales. Gradient field of spatial distribution entropy at each scale is calculated. Local direction information is obtained by smoothing the gradient field. Direction frequency histogram is generated based on the local direction information.

[0040] The initial principal direction is obtained by detecting the peak position of the direction frequency histogram. The direction jump is eliminated based on the cyclic consistency constraint to obtain a continuous principal direction. The continuous principal direction is organized into a principal direction vector at each scale.

[0041] Construct a projection transformation matrix between adjacent scales, map the first scale principal direction vector to the second scale space through the projection transformation matrix, calculate the orthogonal components of the mapped vector and the second scale principal direction vector, construct an orthogonal projection metric based on the orthogonal components, and calculate the initial fusion weights of the spatial distribution entropy of each scale based on the orthogonal projection metric.

[0042] The initial fusion features are obtained by weighting and combining the spatial distribution entropy at each scale using the initial fusion weights; the mutual information between the initial fusion features and the spatial distribution entropy at each scale is calculated, and the fusion weights are iteratively optimized based on the mutual information until convergence, thus obtaining the final fusion weights.

[0043] The final fusion weights are applied to the weighted fusion of spatial distribution entropy at each scale to generate quantitative values ​​of the growth state during the period.

[0044] In one alternative embodiment,

[0045] The spectral distance between pixels is calculated from the temporal feature data. A spatial constraint function is established based on the spectral distance. The spatial constraint function is combined with the temporal feature data to construct a region segmentation criterion. The pixels are then classified and iteratively processed. The output includes crop planting area location information and crop type information.

[0046] Extract the feature changes of adjacent pixels in the temporal feature data at each growth stage, construct a temporal difference weight matrix, and use the temporal difference weight matrix to calculate the spectral distance between pixels.

[0047] An adaptive kernel function is constructed based on the spectral distance. The spatial correlation strength between pixels is calculated based on the adaptive kernel function to generate a spatial constraint function. The spatial constraint function is combined with temporal feature data to construct a region segmentation criterion. The region segmentation criterion includes inter-class feature differences and region spatial constraints.

[0048] Classification samples are extracted from time-series feature data and initial feature centers are calculated. The criterion value of each pixel is calculated according to the region segmentation criterion, and the pixel is assigned to the category with the smallest criterion value. The initial feature centers of each category are updated based on the classification results, and the parameters of the adaptive kernel function are adjusted using the updated feature centers.

[0049] The region segmentation criterion gradient of the classification boundary pixels is recalculated based on the adjusted adaptive kernel function. For boundary pixels that exceed the preset gradient range, the category assignment is re-executed. The iteration is repeated until the region segmentation converges, and the location information of the crop planting area and the crop type information are output.

[0050] A second aspect of this invention provides an intelligent extraction system for crop planting areas based on multispectral remote sensing, comprising:

[0051] The first unit is used to acquire remote sensing image data of multiple bands and perform atmospheric and geometric corrections to obtain corrected remote sensing image data.

[0052] The second unit is used to calculate the vegetation coverage index and vegetation growth vitality index based on the corrected remote sensing image data, calculate the band combination coefficient based on the vegetation coverage index and vegetation growth vitality index, construct a new crop growth characterization index using the band combination coefficient, extract crop canopy structure features based on the crop growth characterization index, and generate a crop feature map.

[0053] The third unit is used to obtain crop feature maps of multiple periods within the crop growth cycle, calculate the spatial distribution entropy of pixel values ​​in the crop feature maps of each period, establish a growth status evaluation function, calculate the feature contribution of each period as the temporal fusion weight based on the growth status evaluation function, and perform weighted combination of crop feature maps of multiple periods to obtain temporal feature data.

[0054] The fourth unit is used to calculate the spectral distance between pixels in the temporal feature data, establish a spatial constraint function based on the spectral distance, combine the spatial constraint function with the temporal feature data to construct a region segmentation criterion, perform iterative classification calculations on pixels, and output the location information of crop planting areas and crop type information.

[0055] A third aspect of the present invention provides an electronic device, comprising:

[0056] processor;

[0057] Memory used to store processor-executable instructions;

[0058] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0059] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0060] In this embodiment, the intelligent extraction method for crop planting areas based on multispectral remote sensing utilizes multispectral remote sensing data to construct specific crop growth characterization indices. This effectively extracts crop canopy structure features, reduces soil background and weed interference, and improves the accuracy of crop identification. By analyzing the spatial distribution entropy of crop feature maps from multiple periods, a growth status evaluation function is constructed to calculate the feature contribution of each period. This achieves adaptive fusion of multi-period data, fully utilizing the temporal feature differences of crops at different growth stages and enhancing the distinguishability of crops from other land features. When classifying temporal feature data, a spatial constraint function based on spectral distance is introduced, considering the spatial correlation between pixels. This effectively reduces noise interference during the classification process, improves the continuity and completeness of crop planting area identification, and achieves accurate extraction of crop planting area location and crop type. Attached Figure Description

[0061] Figure 1 This is a flowchart illustrating the intelligent extraction method for crop planting areas based on multispectral remote sensing, as described in an embodiment of the present invention.

[0062] Figure 2 This is a flowchart of the crop time-series image feature fusion and classification method according to an embodiment of the present invention. Detailed Implementation

[0063] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0064] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0065] Figure 1 This is a flowchart illustrating the intelligent extraction method for crop planting areas based on multispectral remote sensing, as described in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0066] Acquire remote sensing image data in multiple bands and perform atmospheric and geometric corrections to obtain corrected remote sensing image data;

[0067] The vegetation coverage index and vegetation growth vitality index are calculated based on the corrected remote sensing image data. The band combination coefficient is calculated based on the vegetation coverage index and vegetation growth vitality index. A new crop growth characterization index is constructed using the band combination coefficient. The crop canopy structure features are extracted based on the crop growth characterization index to generate a crop feature map.

[0068] The crop feature maps of multiple periods within the crop growth cycle are obtained. The spatial distribution entropy of pixel values ​​in the crop feature maps of each period is calculated. A growth status evaluation function is established. The feature contribution of each period is calculated based on the growth status evaluation function as the temporal fusion weight. The crop feature maps of multiple periods are weighted and combined to obtain temporal feature data.

[0069] The spectral distance between pixels is calculated from the temporal feature data. A spatial constraint function is established based on the spectral distance. The spatial constraint function is combined with the temporal feature data to construct a region segmentation criterion. The pixels are then classified and iteratively processed to output the location information of the crop planting area and the crop type information.

[0070] In one optional implementation, remote sensing image data of multiple bands is acquired and atmospheric and geometric corrections are performed to obtain corrected remote sensing image data, including:

[0071] Remote sensing image data in the visible and infrared bands is acquired using a multispectral sensor. The remote sensing image data includes spectral and spatial information of ground object reflection.

[0072] Atmospheric parameters of the region corresponding to the remote sensing image data are collected, atmospheric influence features are extracted based on the atmospheric parameters, and atmospheric correction is performed on the remote sensing image data based on the atmospheric influence features to obtain spectral reflectance data after eliminating atmospheric influence.

[0073] Obtain uniformly distributed ground control point information within the geographic area corresponding to the remote sensing image data. The ground control point information includes the geographic coordinates of the control points and the corresponding image features. Construct a coordinate mapping relationship based on the ground control point information. Use the coordinate mapping relationship to perform geometric correction on the spectral reflectance data to generate initial corrected data with geographic coordinates.

[0074] The spectral response change between adjacent bands and the spatial position deviation between adjacent pixels are calculated in the initial correction data. Abnormal regions are identified based on the spectral response change and spatial position deviation. The abnormal regions are then corrected to obtain the corrected remote sensing image data.

[0075] In this embodiment, a multispectral sensor is used to acquire remote sensing image data in the visible light and infrared bands. Multispectral sensors can typically acquire image data simultaneously in multiple bands, including blue (0.45-0.52 μm), green (0.52-0.60 μm), red (0.63-0.69 μm), near-infrared (0.76-0.90 μm), and mid-infrared (1.55-1.75 μm). The acquired remote sensing image data contains spectral information about ground features, such as the reflectance characteristics of different ground features in various bands, as well as spatial information, such as the shape, size, texture, and spatial distribution of ground features.

[0076] While acquiring remote sensing image data, atmospheric parameters for the corresponding areas are also collected. Specifically, atmospheric parameters such as temperature, air pressure, humidity, and visibility can be collected in real time through ground meteorological stations, or atmospheric optical thickness can be measured using a solar photometer, or aerosol optical thickness and particle size distribution can be measured using aerosol monitoring equipment. For example, in the study area where data was acquired on May 15, 2023, the collected atmospheric parameters included: temperature 25°C, relative humidity 65%, visibility 15 km, atmospheric optical thickness 0.35, and urban aerosol type. Atmospheric influence characteristics, including atmospheric scattering and absorption effects, are extracted based on these atmospheric parameters.

[0077] To address atmospheric scattering effects, Rayleigh and Mie scattering coefficients are calculated based on atmospheric parameters. For example, the Rayleigh scattering coefficient in the blue light band (0.45-0.52 μm) is 0.018, and the Mie scattering coefficient is 0.035. To address atmospheric absorption effects, atmospheric transmittance for each band is calculated based on the content of gases such as water vapor, ozone, and carbon dioxide in the atmosphere. For example, the atmospheric transmittance in the near-infrared band (0.76-0.90 μm) is 0.82. Based on these atmospheric influence characteristics, atmospheric correction is performed on the remote sensing image data. Specifically, a lookup table-based atmospheric correction method is used. By consulting a pre-defined table of correspondence between atmospheric parameters and surface reflectance, the grayscale values ​​in the original remote sensing image data are converted into surface reflectance. For example, a pixel with an original grayscale value of 85 in the red light band has a surface reflectance of 0.16 after atmospheric correction. Through atmospheric correction, the influence of atmospheric scattering and absorption on the remote sensing image data is eliminated, resulting in spectral reflectance data free from atmospheric effects.

[0078] Ground control point information was obtained within the geographic area corresponding to the remote sensing image data. Twenty-five prominent landmarks were selected as ground control points within the study area using the Global Positioning System (GPS), ensuring their uniform distribution. The geographic coordinates (longitude and latitude) of each control point were recorded. Simultaneously, the pixel locations corresponding to these control points were located in the remote sensing image, and their image features, such as pixel coordinates (row number, column number) and grayscale values, were extracted. For example, the geographic coordinates of the first control point are (116.3°E, 40.0°N), corresponding to image coordinates (row 1205, column 3568); the geographic coordinates of the second control point are (116.5°E, 40.2°N), corresponding to image coordinates (row 2350, column 4125). Based on this ground control point information, a coordinate mapping relationship was constructed, and a quadratic polynomial transformation model was used to establish the correspondence between the image coordinate system and the geographic coordinate system.

[0079] Geometric correction of spectral reflectance data is performed using the established coordinate mapping relationship, mapping pixels in the original image to a new geographic coordinate system through resampling. Specifically, bilinear interpolation resampling is used to calculate the gray value of the target pixel. For example, if the gray values ​​of the four original pixels surrounding the target pixel are 0.18, 0.20, 0.19, and 0.21, respectively, after bilinear interpolation, the gray value of the target pixel is 0.195. Through geometric correction, initial corrected data with standard geographic coordinates is generated, achieving accurate spatial correspondence of the image.

[0080] The initial calibration data undergoes quality assessment and further calibration. The spectral response changes between adjacent bands in the initial calibration data are calculated to detect spectral anomalies. For example, the difference image between the red and near-infrared bands is calculated; standard vegetation areas should exhibit large positive values ​​(e.g., 0.4-0.6). If this value is abnormal in a certain area (e.g., less than 0.1), a spectral anomaly may exist. Simultaneously, the spatial positional deviation between adjacent pixels is calculated to detect geometric correction anomalies. For example, the image gradient is calculated through edge detection; if the gradient value in a certain area is abnormal (e.g., greater than 30), there may be inaccurate geometric correction. Based on the anomalies identified by the spectral response changes and spatial positional deviations, targeted correction processing is applied. For spectral anomalies, correction is performed using a spectral matching method, comparing the spectral characteristics of the anomaly area with those of neighboring normal areas, and adjusting the spectral characteristics of the anomaly area using a weighted average method. For example, if the near-infrared reflectance of the anomaly area is 0.25, the average reflectance of the surrounding normal area is 0.45, and the adjusted reflectance is 0.42. For geometrically anomalous areas, local image registration is used for correction. A feature-point-based registration algorithm is employed to precisely adjust the geometric position of the anomalous area. For example, if the positional deviation of an anomalous area is 3.5 pixels, after local registration, the deviation is reduced to 0.8 pixels. Through these correction processes, high-quality corrected remote sensing image data is obtained, providing an accurate data foundation for subsequent remote sensing applications.

[0081] In this embodiment, comprehensive image data is captured by acquiring multi-band visible and infrared data and combining it with the spectral and spatial information of ground features. Atmospheric influence features are extracted using regional atmospheric parameters, and atmospheric correction is applied to the remote sensing image, significantly eliminating the interference of atmospheric scattering and absorption on spectral reflectance and improving the accuracy of reflectance data. By acquiring uniformly distributed ground control point information and establishing precise coordinate mapping relationships, geometric correction is performed on the spectral reflectance data, ensuring a high degree of matching between the image and geographical location, thus providing an accurate spatial reference for subsequent crop area identification. Furthermore, the spectral response changes and pixel spatial deviations between adjacent bands are calculated to identify and correct abnormal areas, effectively eliminating or correcting noise and error regions, and improving the overall consistency and reliability of the image. These technical effects collectively ensure the accuracy and stability of subsequent crop growth status analysis and intelligent extraction based on multi-temporal and multi-band data.

[0082] In one optional implementation, the vegetation coverage index and vegetation growth vitality index are calculated based on the corrected remote sensing image data, and the band combination coefficients are calculated based on the vegetation coverage index and vegetation growth vitality index, including:

[0083] Collect spectral data of farmland soil in the crop planting area, extract soil background features based on the spectral data, and establish a feature library containing the soil background features;

[0084] Each pixel in the corrected remote sensing image data is matched with the soil background features in the feature library. The corresponding soil type is identified based on the matching results. The spectral response of the pixel is corrected based on the soil type to obtain spectral data that eliminates the influence of soil background and thus obtain vegetation spectral data.

[0085] Red edge location features and chlorophyll absorption features are extracted from the vegetation spectral data. The vegetation coverage index is calculated based on the red edge location features and chlorophyll absorption features. At the same time, the vegetation growth vitality index is calculated using the first derivative features of the vegetation spectral data.

[0086] The vegetation coverage index and vegetation growth vitality index are decomposed into multiple scales respectively. By calculating the complementary information of the vegetation coverage index and vegetation growth vitality index at each scale, the information weight of each scale is determined. The information weight is then normalized to obtain the initial band combination coefficient.

[0087] Canopy spectral data of crops at different growth stages are collected, and the variation of the canopy spectral data with the growth stage is analyzed. The initial band combination coefficients are dynamically corrected according to the variation pattern to obtain the band combination coefficients.

[0088] In one embodiment, to construct a soil background feature library, representative farmland plots are selected within a specific crop planting area, and spectral data of different types of farmland soil are collected using a portable spectrometer. In practice, continuous spectral data with a wavelength range of 400 nm to 2500 nm are collected for each sampling point. Each sampling point is measured five times, and the average value is taken to reduce random errors. The collected raw soil spectral data is preprocessed, including noise removal, smoothing, and standardization. Key features are extracted from the processed soil spectral data using principal component analysis, retaining the top few principal components with a cumulative contribution rate reaching 95% as soil background features. These feature vectors, along with the corresponding soil type information, are stored in a feature library to form a feature library containing soil background features.

[0089] The corrected remote sensing image data undergoes pixel-by-pixel processing to eliminate the influence of soil background on vegetation spectra. For each pixel, its spectral features are extracted and matched with soil background features stored in the feature library. The minimum distance classification method is used to determine the best-matching soil type. For example, for a pixel with an extracted spectral feature vector of [0.23, 0.35, 0.41, 0.54, 0.67], the Euclidean distance between this vector and the feature vectors of each soil type in the feature library is calculated, and the soil type with the smallest distance is found to be "sandy loam". Based on the determined soil type, the spectral response of the pixel is corrected using the soil-adjusted vegetation index method. Taking a practical example, for a pixel identified as having a sandy loam background, its original spectrum has a reflectance of 0.42 in the near-infrared band and 0.18 in the red band. After correction using the soil line equation, the adjusted near-infrared reflectance is 0.39 and the red band reflectance is 0.15, thus eliminating the influence of soil background and obtaining more accurate vegetation spectral data.

[0090] Red-edge location features and chlorophyll absorption features were extracted from the vegetation spectral data after removing the influence of soil background. Red-edge location features refer to the locations where vegetation reflectance increases sharply in the 680nm to 750nm wavelength range. The wavelength corresponding to the maximum first derivative was determined by interpolating and fitting the spectral curve within this range. For example, if the maximum first derivative of the processed vegetation spectrum for a certain pixel occurs at 724nm in the red-edge region, this wavelength is the red-edge location feature value. Chlorophyll absorption features were determined by calculating the absorption depth and area of ​​the spectral curve in the 550nm to 680nm range. Based on the extracted red-edge location features and chlorophyll absorption features, a vegetation coverage index was calculated. Simultaneously, a vegetation growth vitality index was calculated using the first derivative features of the vegetation spectral data in the 400nm to 2500nm range. For example, the calculated vegetation coverage index for a certain pixel was 0.76, and the vegetation growth vitality index was 0.65. These two indices reflect the degree of crop coverage and growth status, respectively.

[0091] The calculated vegetation cover index and vegetation growth vitality index were decomposed into multiple scales using wavelet transform to divide the indices into four scale levels. Complementary information content was calculated for the vegetation cover index and vegetation growth vitality index at each scale to measure the amount of independent information contained in the two indices at that scale. Taking a certain region as an example, the complementary information contents at the four scales were 0.87, 0.65, 0.42, and 0.31, respectively. Based on these values, the information weights at each scale were calculated. After normalizing the information weights, the initial band combination coefficients for the four scales were obtained as 0.39, 0.29, 0.19, and 0.13, respectively. These coefficients reflect the contribution of information at different scales to the final vegetation monitoring.

[0092] To adapt to the changing spectral characteristics of crops at different growth stages, canopy spectral data were collected during key stages of the crop growth season (such as seedling emergence, jointing, heading, and maturity). A hyperspectral imager was used to image the sample plots vertically downwards, acquiring hyperspectral data in the wavelength range of 400 nm to 1000 nm. By analyzing the variation patterns of canopy spectra at different growth stages, a correspondence between growth stages and spectral characteristics was established. For example, during the jointing stage of wheat, near-infrared reflectance increased significantly, while red reflectance decreased, and the red edge shifted approximately 15 nm towards longer wavelengths. Based on these patterns, the initial band combination coefficients were dynamically corrected. Specifically, for wheat at the jointing stage, the weights related to the near-infrared band in the initial band combination coefficients were increased by 10%, and the weights related to the red band were decreased by 5%, resulting in band combination coefficients [0.43, 0.28, 0.18, 0.11] suitable for the current growth stage. This dynamic correction mechanism allows vegetation monitoring indicators to better adapt to the changing characteristics of crops at different growth stages, improving monitoring accuracy.

[0093] In this embodiment, by establishing a soil background feature database and performing soil type matching and correction on remote sensing image pixels, the interference of soil background on vegetation spectra is effectively eliminated, improving the purity and accuracy of vegetation spectral data. By extracting red edge positions and chlorophyll absorption features, vegetation coverage index and growth vitality index are accurately calculated, enhancing the sensitivity and characterization ability of vegetation growth status. Multi-scale decomposition and complementary information analysis enable multi-level information fusion of vegetation indices, rationally allocating information weights at different scales, and improving the scientificity and adaptability of band combination coefficients. By combining the dynamic changes in canopy spectra at different crop growth stages, the band combination coefficients are dynamically corrected, allowing the band weights to adaptively adjust with changes in crop growth, thus more accurately reflecting crop growth status. This significantly improves the accuracy and timeliness of vegetation index calculation, providing a more reliable spectral basis and data support for subsequent crop growth monitoring and regional extraction.

[0094] In one optional implementation, a new crop growth characterization index is constructed using band combination coefficients. Based on this index, crop canopy structure features are extracted to generate a crop feature map, including:

[0095] The spatial distribution gradients of vegetation coverage index and growth vitality index are calculated separately. A regional segmentation threshold is established, and the crop planting area is segmented using the regional segmentation threshold to obtain multiple sub-regions.

[0096] Calculate the statistical characteristics of vegetation coverage index and growth vitality index in each sub-region to obtain the spectral distribution characteristics of the sub-region. Based on the spectral distribution characteristics, perform regional correction on the band combination coefficient to generate regional band combination coefficients.

[0097] The vegetation coverage index and growth vitality index are weighted and fused based on the regional band combination coefficients to construct a crop growth characterization index.

[0098] Wavelet decomposition is performed on the crop growth characterization index to obtain multi-scale wavelet coefficients. The energy distribution of wavelet coefficients at each scale is calculated, and an inter-scale correlation metric is constructed. The combined weights of the coefficients at each scale are determined based on the inter-scale correlation metric. The multi-scale wavelet coefficients are reconstructed using the combined weights to obtain the crop canopy structure characteristics.

[0099] The vertical height distribution and horizontal density distribution are extracted based on the characteristics of the crop canopy structure.

[0100] The height and density distributions of multiple time phases are obtained, the distribution changes between time phases are calculated, and spatiotemporal interpolation is performed based on the distribution changes to generate a crop feature map.

[0101] In one optional embodiment, multispectral remote sensing image data of the crop planting area is first acquired, including red, near-infrared, and green bands. Based on this band data, a vegetation cover index is calculated, such as the Normalized Difference Vegetation Index (NDVI) = (Near-infrared band - Red band) / (Near-infrared band + Red band). Simultaneously, a growth vigor index, such as the Enhanced Vegetation Index (EVI), is calculated. This index incorporates blue band data for atmospheric correction to reduce soil background interference. During the summer maize growing season in a certain experimental area, NDVI values ​​ranged from 0.25 to 0.95, and EVI values ​​ranged from 0.15 to 0.75.

[0102] The spatial distribution gradients of vegetation cover index (NDVI) and growth vigor index (EVI) were calculated. These gradients were obtained by calculating the rates of change of the indices in the horizontal and vertical directions. For NDVI, the gradient rate of change in the horizontal direction was ±0.02 / pixel, and in the vertical direction, it was ±0.025 / pixel. For EVI, the gradient rate of change in the horizontal direction was ±0.015 / pixel, and in the vertical direction, it was ±0.02 / pixel. Based on these gradient values, a region segmentation threshold was set. For example, when the absolute value of the gradient of NDVI was greater than 0.02 and the absolute value of the gradient of EVI was greater than 0.015, that location was considered the region boundary. The crop planting area was segmented using the region segmentation threshold to obtain multiple sub-regions; a total of 12 sub-regions were obtained in the experimental area.

[0103] For each subregion, the statistical characteristics of the vegetation cover index and growth vigor index were calculated, including the mean, standard deviation, maximum, minimum, and median. For example, in subregion 1, the NDVI mean was 0.76, the standard deviation was 0.08, the maximum was 0.92, the minimum was 0.58, and the median was 0.78; the EVI mean was 0.52, the standard deviation was 0.06, the maximum was 0.67, the minimum was 0.38, and the median was 0.53. These statistical characteristics constitute the spectral distribution characteristics of the subregion.

[0104] The band combination coefficients are regionally corrected based on spectral distribution characteristics to generate regional band combination coefficients. The correction process considers the mean and standard deviation of the index and adjusts the band weights accordingly. For example, the original weight of the red band is 0.6 and the weight of the near-infrared band is 0.4. In sub-region 1, based on the ratio of the mean NDVI to the mean EVI (0.76 / 0.52≈1.46), the weight of the red band is adjusted to 0.65 and the weight of the near-infrared band is adjusted to 0.35.

[0105] A crop growth characterization index is constructed by weighting and fusing the vegetation cover index and the growth vigor index using regional band combination coefficients. The specific fusion method is: Growth Characterization Index = Regional Band Combination Coefficient [Red Light] × NDVI + Regional Band Combination Coefficient [Near Infrared] × EVI. In sub-region 1, the growth characterization index = 0.65 × 0.76 + 0.35 × 0.52 = 0.68.

[0106] Crop growth characterization indices were decomposed using wavelet decomposition at three levels, employing the db4 wavelet basis function to obtain multi-scale wavelet coefficients. Each level yielded an approximation coefficient matrix and three detail coefficient matrices, corresponding to detail information in the horizontal, vertical, and diagonal directions, respectively. The energy distribution of the wavelet coefficients at each scale was calculated, with energy values ​​derived from the sum of squared coefficients. In the experimental area data, the energy distribution of the first level decomposition was as follows: approximation coefficients accounted for 85%, horizontal detail coefficients for 6%, vertical detail coefficients for 7%, and diagonal detail coefficients for 2%.

[0107] An inter-scale correlation metric is constructed by calculating the correlation between coefficients at different scales. For example, the correlation coefficient for the horizontal detail coefficients of the first and second layers is 0.82, the correlation coefficient for the vertical detail coefficients is 0.78, and the correlation coefficient for the diagonal detail coefficients is 0.65. The combined weights of the coefficients at each scale are determined based on the inter-scale correlation metric, with the weight values ​​proportional to the energy distribution and correlation. In this embodiment, the weight for the first layer is 0.5, the weight for the second layer is 0.3, and the weight for the third layer is 0.2.

[0108] The crop canopy structure features are obtained by reconstructing multi-scale wavelet coefficients using combined weights. The reconstruction process involves weighted fusion of approximation and detail coefficients from each layer, followed by inverse wavelet transform. The reconstructed canopy structure features contain multi-scale information about the crop's growth status.

[0109] The vertical height distribution and horizontal density distribution of the crop canopy were extracted based on the characteristics of the canopy structure. The height distribution was represented by the rate of change of the reconstructed signal in the vertical direction, and the density distribution was represented by the cumulative value of the reconstructed signal in the horizontal direction. In this example, the maize canopy height distribution ranged from 65 cm to 240 cm, with an average height of 182 cm; the density distribution ranged from 35 plants / m² to 75 plants / m², with an average density of 52 plants / m².

[0110] Data on height and density distributions across multiple time phases were acquired, including the five key growth stages: emergence, jointing, heading, grain filling, and maturity. The distribution changes between these phases were calculated; for example, from the jointing stage to the heading stage, the average height increased by 75 cm, and the density increased by 8 plants per square meter. Based on these distribution changes, spatiotemporal interpolation was performed using cubic spline interpolation to generate continuous crop characteristic maps, including height and density change maps. These maps clearly demonstrate the spatiotemporal dynamics of crop growth and development.

[0111] By calculating the spatial gradients of vegetation cover index and growth vigor index, and utilizing regional segmentation thresholds, detailed division of crop planting areas is achieved, enhancing spectral consistency and spatial coherence within the regions. For the spectral distribution characteristics of different sub-regions, regional band combination coefficients are corrected and weighted fusion is performed, improving the regional adaptability and accuracy of crop growth characterization indices. Wavelet decomposition is used to obtain multi-scale information, and weights are rationally allocated through inter-scale correlation metrics to achieve multi-scale reconstruction of canopy structure features, effectively capturing the spatial complexity of the crop canopy. Furthermore, vertical height distribution and horizontal density distribution are extracted from canopy structure features to reflect the three-dimensional growth status of crops. By analyzing height and density changes across multiple time phases, spatiotemporal interpolation techniques are used to smooth and fill in the data, enabling continuous monitoring and refined characterization of crop growth dynamics.

[0112] like Figure 2 The diagram illustrates the process of fusion and classification of crop time-series image features in this embodiment.

[0113] In one optional implementation, the spatial distribution entropy of pixel values ​​in the crop feature maps of each period is calculated, a growth state evaluation function is established, and the feature contribution of each period is calculated as the temporal fusion weight based on the growth state evaluation function. The crop feature maps of multiple periods are then weighted and combined to obtain the temporal feature data, including:

[0114] Perform multi-scale decomposition on the crop feature map of each period to obtain feature components at different scales, calculate the pixel value distribution probability of each scale feature component, and construct the spatial distribution entropy based on the pixel value distribution probability.

[0115] Analyze the changing trends of spatial distribution entropy at different scales in adjacent periods, extract the principal direction vectors of spatial distribution entropy at different scales, calculate the orthogonal projection metric between scales based on the principal direction vectors, and fuse the spatial distribution entropy at each scale according to the orthogonal projection metric to obtain the quantitative value of the growth state of the period.

[0116] Calculate the spatial autocorrelation coefficient of local regions in the crop feature map, construct a growth status evaluation function based on the spatial autocorrelation coefficient, input the quantitative value of the growth status of the period into the growth status evaluation function, and calculate the feature contribution of the crop feature map of each period.

[0117] The gradient direction field of crop feature maps from adjacent periods is extracted, and the gradient direction field is used to construct a temporal variation tensor to dynamically adjust the feature contribution and generate temporal fusion weights.

[0118] The crop feature maps from multiple periods are weighted and combined based on the temporal fusion weight to obtain a fused feature map. The anisotropy index of the fused feature map is calculated to construct a regional smoothing constraint. The fused feature map is then subjected to temporal analysis based on the regional smoothing constraint to generate temporal feature data.

[0119] For example, firstly, crop feature maps at multiple growth stages are acquired. These feature maps can be multispectral images acquired by satellites or drones. In this embodiment, feature maps of wheat from sowing to maturity in a certain area are collected, with a resolution of 10 meters and a coverage area of ​​100 hectares. Multi-scale decomposition is performed on the crop feature map of each stage. Haar wavelet transform is used to decompose the feature map into four scale levels to obtain feature components at different scales. For the feature map of the first growth stage, the first-scale feature component reflects the field boundary information, the second-scale feature component reflects the large-scale farmland structure, the third-scale feature component reflects the changes in medium-sized growth areas, and the fourth-scale feature component reflects the subtle crop growth status. Taking the first growth stage as an example, the feature map is divided into a 256×256 pixel grid. A pixel value histogram is calculated for each grid, and the pixel value range [0, 255] is divided into 32 intervals. The number of pixels in each interval is counted and normalized to obtain the pixel value distribution probability. Based on these probability distributions, the spatial distribution entropy of the feature components at each scale is calculated. For example, the spatial distribution entropy at the first scale is 4.82, at the second scale it is 5.36, at the third scale it is 6.14, and at the fourth scale it is 3.89.

[0120] The changing trends of spatial distribution entropy at various scales in adjacent periods were analyzed. Comparing the spatial distribution entropy of the first and second growth periods, the changes were found to be +0.56 at the first scale, +0.78 at the second scale, +1.25 at the third scale, and +0.43 at the fourth scale. These changes constitute a four-dimensional vector. Principal component analysis was used to extract the principal direction vector, which was found to be [0.23, 0.32, 0.51, 0.18]. For the comparison between the second and third growth periods, the principal direction vector was [0.19, 0.37, 0.45, 0.24]. The orthogonal projection metric between these two principal direction vectors was calculated, yielding a value of 0.89, indicating that the direction of growth state changes was basically consistent. Based on the orthogonal projection metric, the spatial distribution entropy at each scale was weighted and fused, resulting in a state quantification value of 5.21 for the first growth period, 6.37 for the second period, 7.18 for the third period, 7.85 for the fourth period, and 7.49 for the fifth period.

[0121] The spatial autocorrelation coefficient of local regions in the crop feature map was calculated. The feature map was divided into 50×50 pixel local regions, and the Moran index was calculated as the spatial autocorrelation coefficient for each region. Taking the third growth stage as an example, the calculated average spatial autocorrelation coefficient was 0.78, indicating that the pixel values ​​have strong spatial clustering and the crop growth is relatively uniform. A growth status evaluation function was constructed based on the spatial autocorrelation coefficient. This function combines the spatial autocorrelation coefficient with the quantitative value of the growth status of each stage, in the form of a weighted product of the two. For the first growth stage, the spatial autocorrelation coefficient was 0.52, the evaluation function result was 2.71, and the corresponding feature contribution was 0.15. For the second growth stage, the spatial autocorrelation coefficient was 0.67, the evaluation function result was 4.27, and the corresponding feature contribution was 0.23. For the third growth stage, the evaluation function result was 5.60, and the feature contribution was 0.30. For the fourth growth stage, the evaluation function result was 4.71, and the feature contribution was 0.25. For the fifth growth stage, the evaluation function result is 1.35, and the feature contribution is 0.07.

[0122] Gradient orientation fields were extracted from crop feature maps of adjacent periods. The Sobel operator was used to calculate the horizontal and vertical gradients of the feature maps, and then the gradient direction and intensity for each pixel were calculated. Taking the comparison between the second and third growth periods as an example, the consistency of the gradient orientation fields was calculated, revealing that 73% of the pixels had similar gradient directions, indicating strong directional consistency in crop growth changes. A temporal variation tensor was constructed using the gradient orientation fields. The tensor's principal eigenvalue was 5.34, secondary eigenvalue was 2.18, and eigenvalue ratio was 2.45, indicating a clear principal directionality in the changes. Based on the eigenvalue structure of the temporal variation tensor, the feature contribution was dynamically adjusted, increasing the feature contribution of the third growth period by 0.05, resulting in the adjusted temporal fusion weights: 0.13 for the first period, 0.21 for the second period, 0.35 for the third period, 0.24 for the fourth period, and 0.07 for the fifth period.

[0123] A weighted combination of crop feature maps from multiple periods is performed based on temporal fusion weights. For each pixel position in the feature map, a weighted average of the feature values ​​from five periods is calculated to generate a fused feature map. The anisotropy index of the fused feature map is calculated, and the feature map is divided into 20×20 pixel analysis windows. Within each window, the directional distribution statistics of pixel values ​​are calculated, resulting in an average anisotropy index of 0.67, indicating that the fused feature map has a clear directional structure. Based on the anisotropy index, a regional smoothing constraint is constructed. Regions with high anisotropy indices maintain clear edges, while regions with low anisotropy indices are moderately smoothed. The smoothing radius increases as the anisotropy index decreases, with a maximum smoothing radius of 5 pixels. After applying the regional smoothing constraint to the fused feature map, the rate of change and acceleration over time are calculated to generate temporal feature data characterizing the dynamic characteristics of crop growth, including growth rate, growth stability, and growth cycle features. The final generated temporal feature data shows that the crops in the test area reach their optimal growth state in the third growth period, with the maximum growth rate occurring between the second and third periods, and an average growth rate of 0.032 units / day.

[0124] By calculating the spatial distribution entropy of pixel values ​​in crop feature maps, the growth status of crops at different time points is quantified, effectively reflecting temporal variation characteristics. Multi-scale decomposition is used to obtain spatial information at different scales, and combined with principal direction vectors and orthogonal projection measures, the fusion of multi-scale spatial distribution entropy is achieved, improving the ability to quantify growth status with fine granularity. A growth status evaluation function is constructed using spatial autocorrelation coefficients, making the calculation of feature contribution more spatially consistent and scientific. The introduction of gradient direction fields and temporal variation tensors to dynamically adjust the fusion weights enhances the ability to capture the inherent changing trends of temporal data. Regional smoothing constraints are constructed based on anisotropy indices, improving the spatial continuity and stability of the fused feature maps and further optimizing the temporal analysis results.

[0125] In one optional implementation, the variation trend of spatial distribution entropy at each scale in adjacent periods is analyzed, the principal direction vector of spatial distribution entropy at different scales is extracted, the orthogonal projection metric between scales is calculated based on the principal direction vector, and the spatial distribution entropy at each scale is fused according to the orthogonal projection metric to obtain the quantitative value of the period growth state, including:

[0126] Multi-scale decomposition of spatial distribution entropy in adjacent periods yields spatial distribution entropy at different scales. Gradient field of spatial distribution entropy at each scale is calculated. Local direction information is obtained by smoothing the gradient field. Direction frequency histogram is generated based on the local direction information.

[0127] The initial principal direction is obtained by detecting the peak position of the direction frequency histogram. The direction jump is eliminated based on the cyclic consistency constraint to obtain a continuous principal direction. The continuous principal direction is organized into a principal direction vector at each scale.

[0128] Construct a projection transformation matrix between adjacent scales, map the first scale principal direction vector to the second scale space through the projection transformation matrix, calculate the orthogonal components of the mapped vector and the second scale principal direction vector, construct an orthogonal projection metric based on the orthogonal components, and calculate the initial fusion weights of the spatial distribution entropy of each scale based on the orthogonal projection metric.

[0129] The initial fusion features are obtained by weighting and combining the spatial distribution entropy at each scale using the initial fusion weights; the mutual information between the initial fusion features and the spatial distribution entropy at each scale is calculated, and the fusion weights are iteratively optimized based on the mutual information until convergence, thus obtaining the final fusion weights.

[0130] The final fusion weights are applied to the weighted fusion of spatial distribution entropy at each scale to generate quantitative values ​​of the growth state during the period.

[0131] In this embodiment, when performing multi-scale decomposition of the spatial distribution entropy of adjacent periods, wavelet transform is used to decompose the spatial distribution entropy into multiple scales. For example, the spatial distribution entropy can be decomposed into four scales, corresponding to frequency ranges of 0-0.125Hz, 0.125-0.25Hz, 0.25-0.5Hz, and 0.5-1Hz, respectively. At each scale, the spatial distribution entropy is represented as a two-dimensional matrix, reflecting the complexity of the spatial distribution at that scale. For the first-scale spatial distribution entropy matrix E1(i,j) of the first period and the first-scale spatial distribution entropy matrix E2(i,j) of the second period, the changing trend can be analyzed by calculating their difference matrix ΔE(i,j) = E2(i,j) - E1(i,j).

[0132] When calculating the gradient field of the spatial distribution entropy at each scale, gradient operators in the horizontal and vertical directions are used. For the spatial distribution entropy matrix Ek at scale k, the horizontal gradient Gx and vertical gradient Gy are obtained by calculating the differences between adjacent elements in the horizontal and vertical directions, respectively. To reduce the influence of noise, a Gaussian filter is applied to smooth the gradient field, with the filter kernel size set to 5×5 and the standard deviation set to 1.2. After smoothing, the local orientation information θ(i,j) can be calculated using the arctangent function of the gradient, representing the local orientation at each pixel location.

[0133] When generating a direction frequency histogram based on local direction information, the direction angle range from 0° to 360° is divided into 36 intervals, each with a width of 10°. The direction frequency histogram H(θ) is formed by counting the number of pixels within each direction interval. For example, when the local direction is mainly distributed within the range of 45°±15° at a certain scale, the histogram will show a significant peak in this interval. When detecting the peak positions of the direction frequency histogram to obtain the initial principal direction, a peak detection algorithm is used to identify the local maximum values ​​in the histogram. A threshold of 30% of the maximum value of the histogram is set, and peaks below this threshold are filtered out. The initial principal direction is represented by the angle value α; for example, the initial principal directions detected at a certain scale are 43°, 126°, and 215°.

[0134] When eliminating directional jumps to continuous principal directions based on cyclic consistency constraints, the continuity of principal directions between adjacent periods is considered. A directional jump is identified when the difference between the principal directions of adjacent periods exceeds 90°. The direction value is adjusted by adding or subtracting 180° to ensure a smooth change in principal direction over time. For example, if the principal direction of the first period is 30° and the initial principal direction of the second period is 210°, the adjusted principal direction of the second period is 30°, maintaining directional continuity. When organizing continuous principal directions into principal direction vectors at various scales, for each scale of each period, the most significant principal direction (corresponding to the highest peak in the histogram) is selected to form the principal direction vector v. For example, the principal direction vector of period t at four scales can be represented as v. t = [42°, 45°, 47°, 43°], which indicates that the main directions at different scales during this period have a certain degree of consistency.

[0135] When constructing the projection transformation matrix between adjacent scales, the rotation matrix is ​​calculated using the principal direction angle. For scale k and scale k+1, the projection transformation matrix T passes through the principal direction angle difference Δθ = θ. k+1 -θ k Calculated. For example, when θ k =45°, θ k+1 When the angle is 50°, the projection transformation matrix represents a 5° rotation transformation from scale k to scale k+1.

[0136] The principal direction vector of the first scale is mapped to the second scale space through a projection transformation matrix. When calculating the orthogonal component between the mapped vector and the principal direction vector of the second scale, the principal direction angle is first converted into a unit vector form. For example, the unit vector corresponding to the angle 45° is [0.7071, 0.7071]. The first-scale unit vector is then mapped to the second-scale space through the projection transformation matrix to obtain the mapped vector. Then, the magnitude of the orthogonal component between the mapped vector and the second-scale principal direction unit vector is calculated; the smaller this value, the more consistent the principal directions of the two scales are.

[0137] An orthogonal projection metric is constructed based on orthogonal components. The orthogonal component values ​​are normalized to the range [0, 1] to obtain the orthogonal projection metric OPM(k, k+1) between scale k and scale k+1. For example, when the calculated orthogonal component value is 0.087, the orthogonal projection metric is 0.913, indicating that the principal directions of the two scales are highly consistent. When calculating the initial fusion weights of the spatial distribution entropy at each scale based on the orthogonal projection metric, the orthogonal projection metric is regarded as a similarity metric between adjacent scales. The initial fusion weight w is proportional to the orthogonal projection metric. For example, when the orthogonal projection metrics between the four scales are 0.92, 0.88, and 0.85, the initial weights can be set to w = [0.28, 0.26, 0.24, 0.22], and the sum of the weights is 1.

[0138] The initial fusion features are obtained by weighting the spatial distribution entropy at each scale using initial fusion weights. For the four scale spatial distribution entropies E1, E2, E3, and E4 at time t, the initial fusion feature F = 0.28×E1 + 0.26×E2 + 0.24×E3 + 0.22×E4.

[0139] When calculating the mutual information between the initial fusion features and the spatial distribution entropy at each scale, the probability distribution is estimated using a histogram method. The value ranges of the fusion features and the spatial distribution entropy at each scale are divided into 32 equally wide intervals, and the frequency of each interval is counted as a probability estimate. Mutual information represents the statistical correlation between the fusion features and the spatial distribution entropy at each scale; a larger value indicates a stronger correlation. When iteratively optimizing the fusion weights based on the mutual information, the gradient ascent method is used. In each iteration, the weight values ​​are adjusted according to the gradient of the mutual information, with a step size of 0.01. Convergence is determined when the weight change between two consecutive iterations is less than a preset threshold of 0.0001. For example, after 25 iterations, the final optimized weights are w = [0.32, 0.28, 0.23, 0.17].

[0140] The final fusion weights are applied to the weighted fusion of spatial distribution entropy at each scale to generate a quantitative value of the growth state for each period: Q = 0.32×E1 + 0.28×E2 + 0.23×E3 + 0.17×E4. This value reflects the overall growth state of the current period and can be used for subsequent time series analysis and prediction. In practical applications, the quantitative values ​​of the growth state for 10 consecutive periods may be [0.412, 0.438, 0.455, 0.487, 0.521, 0.563, 0.592, 0.608, 0.615, 0.621], showing a continuously rising growth trend.

[0141] By analyzing the changing trends of spatial distribution entropy through multi-scale decomposition, the principal direction vectors at different scales are accurately extracted, effectively capturing the directional characteristics of spatial structures. Utilizing directional frequency histograms and cyclic consistency constraints, the continuity and smoothness of the principal directions are achieved, avoiding errors caused by directional jumps. By constructing inter-scale projection transformation matrices and calculating the orthogonal projection metric of the principal direction vectors, a scientific basis for weight allocation is provided for spatial distribution entropy fusion. The weighted combination of initial fusion weights, combined with iterative optimization of mutual information, ensures that the fusion process takes into account the complementarity and correlation of information at each scale, achieving dynamic adaptive adjustment of weights. Finally, the spatial distribution entropy at each scale is weighted and fused to generate accurate quantitative values ​​of the growth state at different periods. This significantly improves the accuracy and stability of multi-scale spatial information fusion, enhances the ability to depict the details of crop growth states, and provides effective data support and analytical tools for accurately assessing crop growth processes.

[0142] In one optional implementation, the spectral distance between pixels is calculated from the temporal feature data, a spatial constraint function is established based on the spectral distance, the spatial constraint function is combined with the temporal feature data to construct a region segmentation criterion, and the pixels are classified and iteratively processed to output crop planting area location information and crop type information, including:

[0143] Extract the feature changes of adjacent pixels in the temporal feature data at each growth stage, construct a temporal difference weight matrix, and use the temporal difference weight matrix to calculate the spectral distance between pixels.

[0144] An adaptive kernel function is constructed based on the spectral distance. The spatial correlation strength between pixels is calculated based on the adaptive kernel function to generate a spatial constraint function. The spatial constraint function is combined with temporal feature data to construct a region segmentation criterion. The region segmentation criterion includes inter-class feature differences and region spatial constraints.

[0145] Classification samples are extracted from time-series feature data and initial feature centers are calculated. The criterion value of each pixel is calculated according to the region segmentation criterion, and the pixel is assigned to the category with the smallest criterion value. The initial feature centers of each category are updated based on the classification results, and the parameters of the adaptive kernel function are adjusted using the updated feature centers.

[0146] The region segmentation criterion gradient of the classification boundary pixels is recalculated based on the adjusted adaptive kernel function. For boundary pixels that exceed the preset gradient range, the category assignment is re-executed. The iteration is repeated until the region segmentation converges, and the location information of the crop planting area and the crop type information are output.

[0147] In this implementation, when extracting the characteristic changes of adjacent pixels in the temporal feature data at each growth stage, multi-temporal remote sensing image data of the crop growth cycle is first acquired, including biophysical parameters such as Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Chlorophyll Index (CI) at different stages. For a 4×4 pixel sample area, the NDVI value of each pixel at six growth stages (e.g., seedling stage, tillering stage, jointing stage, heading stage, grain filling stage, and maturity stage) is extracted. Taking the pixel at position (2,2) as the center, the NDVI difference of the eight adjacent pixels at each stage is calculated to generate a temporal difference matrix. For example, the NDVI values ​​of the center pixel in six periods are [0.25, 0.42, 0.68, 0.85, 0.76, 0.45], and the NDVI values ​​of its right-side neighboring pixels are [0.26, 0.44, 0.69, 0.83, 0.74, 0.44]. The calculated temporal difference values ​​are [0.01, 0.02, 0.01, -0.02, -0.02, -0.01]. The temporal difference values ​​are calculated for each of the eight neighboring pixels, and the difference weights are determined by combining these values ​​to construct a temporal difference weight matrix. Periods with smaller differences are assigned higher weights (e.g., 0.25), and periods with larger differences are assigned lower weights (e.g., 0.10), with a total weight of 1. The spectral distance between pixels is calculated using this weight matrix. For example, if the NDVI values ​​of two pixels differ by [0.02, 0.03, 0.05, 0.04, 0.03, 0.02] at six different times, the corresponding weights are [0.20, 0.15, 0.10, 0.15, 0.20, 0.20], and the calculated weighted spectral distance is 0.0305.

[0148] When constructing the adaptive kernel function based on the calculated spectral distance, the Gaussian kernel function is chosen as the basic form, and the kernel function parameters are related to the spectral and spatial distances between pixels. For two pixels with a spectral distance of 0.0305 and a spatial distance of 1.414 (diagonally adjacent), the initial kernel width parameter is set to 5.0, and the calculated spatial association strength is 0.962. Similar calculations are performed on all pixel pairs in the image to generate a spatial association strength matrix. Based on the spatial association strength matrix, a spatial constraint function is constructed, which represents the degree of spatial consistency of pixels within a region. In a certain 4×4 region, the spatial association strengths between the center pixel (2,2) and the surrounding pixels are [0.985, 0.975, 0.962, 0.978, 0.960, 0.958, 0.950, 0.945], and the calculated spatial constraint is 0.964. The spatial constraint function is combined with temporal feature data to construct a region segmentation criterion, which includes two parts: inter-class feature difference and region spatial constraint. The inter-class feature difference reflects the degree of difference between pixel features and the center of each class, while the region spatial constraint reflects the degree of spatial clustering of pixels within a region. These two parts are combined with a weight ratio of 0.65:0.35 to form the final region segmentation criterion.

[0149] In the iterative classification process, classification samples are first extracted from the temporal feature data. For example, for a region containing three crops—rice, corn, and soybean—ten pixels of known categories are selected as training samples for each crop. The average temporal feature value of each category is calculated as the initial feature center. For example, the initial feature centers for rice are [0.22, 0.45, 0.72, 0.88, 0.75, 0.40], for corn are [0.18, 0.35, 0.65, 0.82, 0.70, 0.32], and for soybeans are [0.20, 0.38, 0.58, 0.75, 0.68, 0.35]. For each pixel in the image, its region segmentation criterion value relative to each category center is calculated. For example, a pixel with temporal features of [0.21, 0.43, 0.70, 0.86, 0.73, 0.38] has criterion values ​​of 0.042, 0.087, and 0.065 for the rice, corn, and soybean categories, respectively. Therefore, this pixel is assigned to the rice category, which has the smallest criterion value. After the initial classification of all pixels, the feature centers for each category are recalculated. For example, the average temporal features of all pixels assigned to the rice category are [0.23, 0.46, 0.71, 0.87, 0.74, 0.41], which will be used as the updated feature centers for the rice category.

[0150] Based on the updated category feature centers, adjust the parameters of the adaptive kernel function. For boundary regions with significant feature changes, decrease the kernel width parameter to 4.0 to enhance boundary discrimination; for regions with gradual feature changes, increase the kernel width parameter to 6.0 to improve region consistency. Recalculate the region segmentation criterion gradient of the classification boundary pixels using the adjusted kernel function. A boundary pixel is defined as a pixel within its 8-neighborhood that belongs to a different category. Calculate the criterion gradient for each boundary pixel. For example, if a boundary pixel has a criterion value of 0.055 in its current category (rice) and a criterion value of 0.060 in its adjacent possible category (soybean), the criterion gradient is 0.005. Set a gradient threshold of 0.01. Since the criterion gradient of this pixel is less than the threshold, its current category remains unchanged. For boundary pixels with criterion gradients greater than the threshold, reassign them to the category with the smallest criterion value.

[0151] Repeat the above steps of feature center update, kernel function parameter adjustment, and boundary pixel reclassification until the rate of change of classification results between two consecutive iterations is less than 1%, or the maximum number of iterations (20) is reached. The final output classification results include crop planting area location information (represented by pixel coordinates) and crop type information (rice, corn, soybean, etc.). For example, pixel (100, 150) in the image is classified as rice, and pixel (200, 250) is classified as corn. These results can be used to generate crop distribution maps and area statistics tables.

[0152] By calculating the feature changes of adjacent pixels in temporal feature data, a temporal difference weight matrix is ​​constructed to effectively quantify the spectral distance between pixels, enhancing the perception of temporal dynamic changes between pixels. An adaptive kernel function is designed based on spectral distance, and a spatial constraint function is generated by combining spatial correlation strength, achieving a reasonable modeling of spatial continuity and proximity relationships. The spatial constraint function is fused with temporal feature data to construct a region segmentation criterion that includes inter-class differences and spatial constraints, improving classification accuracy and boundary clarity. Dynamic adaptive optimization of the classification process is achieved by iteratively updating the feature centers of the classification samples and adjusting the parameters of the adaptive kernel function. Gradient analysis and dynamic redistribution are performed on the classification boundary pixels to effectively avoid the accumulation of classification errors and ensure the stability and convergence of region division. The overall scheme significantly improves the accuracy and robustness of crop planting area and crop type identification, providing efficient and reliable spatial-temporal classification technology support for precision agriculture management.

[0153] A second aspect of this invention provides an intelligent extraction system for crop planting areas based on multispectral remote sensing, the system comprising:

[0154] The first unit is used to acquire remote sensing image data of multiple bands and perform atmospheric and geometric corrections to obtain corrected remote sensing image data.

[0155] The second unit is used to calculate the vegetation coverage index and vegetation growth vitality index based on the corrected remote sensing image data, calculate the band combination coefficient based on the vegetation coverage index and vegetation growth vitality index, construct a new crop growth characterization index using the band combination coefficient, extract crop canopy structure features based on the crop growth characterization index, and generate a crop feature map.

[0156] The third unit is used to obtain crop feature maps of multiple periods within the crop growth cycle, calculate the spatial distribution entropy of pixel values ​​in the crop feature maps of each period, establish a growth status evaluation function, calculate the feature contribution of each period as the temporal fusion weight based on the growth status evaluation function, and perform weighted combination of crop feature maps of multiple periods to obtain temporal feature data.

[0157] The fourth unit is used to calculate the spectral distance between pixels in the temporal feature data, establish a spatial constraint function based on the spectral distance, combine the spatial constraint function with the temporal feature data to construct a region segmentation criterion, perform iterative classification calculations on pixels, and output the location information of crop planting areas and crop type information.

[0158] A third aspect of the present invention provides an electronic device, comprising:

[0159] processor;

[0160] Memory used to store processor-executable instructions;

[0161] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0162] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0163] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0164] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for intelligent extraction of crop planting areas based on multispectral remote sensing, characterized in that, include: Acquire remote sensing image data in multiple bands and perform atmospheric and geometric corrections to obtain corrected remote sensing image data; The vegetation coverage index and vegetation growth vitality index are calculated based on the corrected remote sensing image data. The band combination coefficient is calculated based on the vegetation coverage index and vegetation growth vitality index. A new crop growth characterization index is constructed using the band combination coefficient. The crop canopy structure features are extracted based on the crop growth characterization index to generate a crop feature map. The crop feature maps of multiple periods within the crop growth cycle are obtained. The spatial distribution entropy of pixel values ​​in the crop feature maps of each period is calculated. A growth status evaluation function is established. The feature contribution of each period is calculated based on the growth status evaluation function as the temporal fusion weight. The crop feature maps of multiple periods are weighted and combined to obtain temporal feature data. The spectral distance between pixels is calculated from the temporal feature data. A spatial constraint function is established based on the spectral distance. The spatial constraint function is combined with the temporal feature data to construct a region segmentation criterion. The pixels are then classified and iteratively processed to output the location information of the crop planting area and the crop type information.

2. The method according to claim 1, characterized in that, Multiple bands of remote sensing image data were acquired and atmospheric and geometric corrections were performed to obtain the corrected remote sensing image data, including: Remote sensing image data in the visible and infrared bands is acquired using a multispectral sensor. The remote sensing image data includes spectral and spatial information of ground object reflection. Atmospheric parameters of the region corresponding to the remote sensing image data are collected, atmospheric influence features are extracted based on the atmospheric parameters, and atmospheric correction is performed on the remote sensing image data based on the atmospheric influence features to obtain spectral reflectance data after eliminating atmospheric influence. Obtain uniformly distributed ground control point information within the geographic area corresponding to the remote sensing image data. The ground control point information includes the geographic coordinates of the control points and the corresponding image features. Construct a coordinate mapping relationship based on the ground control point information. Use the coordinate mapping relationship to perform geometric correction on the spectral reflectance data to generate initial corrected data with geographic coordinates. The spectral response change between adjacent bands and the spatial position deviation between adjacent pixels are calculated in the initial correction data. Abnormal regions are identified based on the spectral response change and spatial position deviation. The abnormal regions are then corrected to obtain the corrected remote sensing image data.

3. The method according to claim 1, characterized in that, The vegetation cover index and vegetation growth vitality index are calculated based on the corrected remote sensing image data. The band combination coefficients are then calculated based on these indices, including: Collect spectral data of farmland soil in the crop planting area, extract soil background features based on the spectral data, and establish a feature library containing the soil background features; Each pixel in the corrected remote sensing image data is matched with the soil background features in the feature library. The corresponding soil type is identified based on the matching results. The spectral response of the pixel is corrected based on the soil type to obtain spectral data that eliminates the influence of soil background and thus obtain vegetation spectral data. Red edge location features and chlorophyll absorption features are extracted from the vegetation spectral data. The vegetation coverage index is calculated based on the red edge location features and chlorophyll absorption features. At the same time, the vegetation growth vitality index is calculated using the first derivative features of the vegetation spectral data. The vegetation coverage index and vegetation growth vitality index are decomposed into multiple scales respectively. By calculating the complementary information of the vegetation coverage index and vegetation growth vitality index at each scale, the information weight of each scale is determined. The information weight is then normalized to obtain the initial band combination coefficient. Canopy spectral data of crops at different growth stages are collected, and the variation of the canopy spectral data with the growth stage is analyzed. The initial band combination coefficients are dynamically corrected according to the variation pattern to obtain the band combination coefficients.

4. The method according to claim 1, characterized in that, New crop growth characterization indices are constructed using band combination coefficients. Crop canopy structure features are then extracted based on these indices, generating crop feature maps, including: The spatial distribution gradients of vegetation coverage index and growth vitality index are calculated separately. A regional segmentation threshold is established, and the crop planting area is segmented using the regional segmentation threshold to obtain multiple sub-regions. Calculate the statistical characteristics of vegetation coverage index and growth vitality index in each sub-region to obtain the spectral distribution characteristics of the sub-region. Based on the spectral distribution characteristics, perform regional correction on the band combination coefficient to generate regional band combination coefficients. The vegetation coverage index and growth vitality index are weighted and fused based on the regional band combination coefficients to construct a crop growth characterization index. Wavelet decomposition is performed on the crop growth characterization index to obtain multi-scale wavelet coefficients. The energy distribution of wavelet coefficients at each scale is calculated, and an inter-scale correlation metric is constructed. The combined weights of the coefficients at each scale are determined based on the inter-scale correlation metric. The multi-scale wavelet coefficients are reconstructed using the combined weights to obtain the crop canopy structure characteristics. The vertical height distribution and horizontal density distribution are extracted based on the characteristics of the crop canopy structure. The height and density distributions of multiple time phases are obtained, the distribution changes between time phases are calculated, and spatiotemporal interpolation is performed based on the distribution changes to generate a crop feature map.

5. The method according to claim 1, characterized in that, The spatial distribution entropy of pixel values ​​in crop feature maps at different periods is calculated, a growth state evaluation function is established, and the feature contribution of each period is calculated as the temporal fusion weight based on the growth state evaluation function. The crop feature maps from multiple periods are then weighted and combined to obtain the temporal feature data, including: Perform multi-scale decomposition on the crop feature map of each period to obtain feature components at different scales, calculate the pixel value distribution probability of each scale feature component, and construct the spatial distribution entropy based on the pixel value distribution probability. Analyze the changing trends of spatial distribution entropy at different scales in adjacent periods, extract the principal direction vectors of spatial distribution entropy at different scales, calculate the orthogonal projection metric between scales based on the principal direction vectors, and fuse the spatial distribution entropy at each scale according to the orthogonal projection metric to obtain the quantitative value of the growth state of the period. Calculate the spatial autocorrelation coefficient of local regions in the crop feature map, construct a growth status evaluation function based on the spatial autocorrelation coefficient, input the quantitative value of the growth status of the period into the growth status evaluation function, and calculate the feature contribution of the crop feature map of each period. The gradient direction field of crop feature maps from adjacent periods is extracted, and the gradient direction field is used to construct a temporal variation tensor to dynamically adjust the feature contribution and generate temporal fusion weights. The crop feature maps from multiple periods are weighted and combined based on the temporal fusion weight to obtain a fused feature map. The anisotropy index of the fused feature map is calculated to construct a regional smoothing constraint. The fused feature map is then subjected to temporal analysis based on the regional smoothing constraint to generate temporal feature data.

6. The method according to claim 5, characterized in that, Analyze the changing trends of spatial distribution entropy at various scales in adjacent periods, extract the principal direction vectors of spatial distribution entropy at different scales, calculate the orthogonal projection metric between scales based on the principal direction vectors, and fuse the spatial distribution entropy at each scale according to the orthogonal projection metric to obtain the quantitative value of the growth state of each period, including: Multi-scale decomposition of spatial distribution entropy in adjacent periods yields spatial distribution entropy at different scales. Gradient field of spatial distribution entropy at each scale is calculated. Local direction information is obtained by smoothing the gradient field. Direction frequency histogram is generated based on the local direction information. The initial principal direction is obtained by detecting the peak position of the direction frequency histogram. The direction jump is eliminated based on the cyclic consistency constraint to obtain a continuous principal direction. The continuous principal direction is organized into a principal direction vector at each scale. Construct a projection transformation matrix between adjacent scales, map the first scale principal direction vector to the second scale space through the projection transformation matrix, calculate the orthogonal components of the mapped vector and the second scale principal direction vector, construct an orthogonal projection metric based on the orthogonal components, and calculate the initial fusion weights of the spatial distribution entropy of each scale based on the orthogonal projection metric. The initial fusion features are obtained by weighting and combining the spatial distribution entropy at each scale using the initial fusion weights; the mutual information between the initial fusion features and the spatial distribution entropy at each scale is calculated, and the fusion weights are iteratively optimized based on the mutual information until convergence, thus obtaining the final fusion weights. The final fusion weights are applied to the weighted fusion of spatial distribution entropy at each scale to generate quantitative values ​​of the growth state during the period.

7. The method according to claim 1, characterized in that, The spectral distance between pixels is calculated from the temporal feature data. A spatial constraint function is established based on the spectral distance. The spatial constraint function is combined with the temporal feature data to construct a region segmentation criterion. The pixels are then classified and iteratively processed. The output includes crop planting area location information and crop type information. Extract the feature changes of adjacent pixels in the temporal feature data at each growth stage, construct a temporal difference weight matrix, and use the temporal difference weight matrix to calculate the spectral distance between pixels. An adaptive kernel function is constructed based on the spectral distance. The spatial correlation strength between pixels is calculated based on the adaptive kernel function to generate a spatial constraint function. The spatial constraint function is combined with temporal feature data to construct a region segmentation criterion. The region segmentation criterion includes inter-class feature differences and region spatial constraints. Classification samples are extracted from time-series feature data and initial feature centers are calculated. The criterion value of each pixel is calculated according to the region segmentation criterion, and the pixel is assigned to the category with the smallest criterion value. The initial feature centers of each category are updated based on the classification results, and the parameters of the adaptive kernel function are adjusted using the updated feature centers. The region segmentation criterion gradient of the classification boundary pixels is recalculated based on the adjusted adaptive kernel function. For boundary pixels that exceed the preset gradient range, the category assignment is re-executed. The iteration is repeated until the region segmentation converges, and the location information of the crop planting area and the crop type information are output.

8. A smart extraction system for crop planting areas based on multispectral remote sensing, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to acquire remote sensing image data of multiple bands and perform atmospheric and geometric corrections to obtain corrected remote sensing image data. The second unit is used to calculate the vegetation coverage index and vegetation growth vitality index based on the corrected remote sensing image data, calculate the band combination coefficient based on the vegetation coverage index and vegetation growth vitality index, construct a new crop growth characterization index using the band combination coefficient, extract crop canopy structure features based on the crop growth characterization index, and generate a crop feature map. The third unit is used to obtain crop feature maps of multiple periods within the crop growth cycle, calculate the spatial distribution entropy of pixel values ​​in the crop feature maps of each period, establish a growth status evaluation function, calculate the feature contribution of each period as the temporal fusion weight based on the growth status evaluation function, and perform weighted combination of crop feature maps of multiple periods to obtain temporal feature data. The fourth unit is used to calculate the spectral distance between pixels in the temporal feature data, establish a spatial constraint function based on the spectral distance, combine the spatial constraint function with the temporal feature data to construct a region segmentation criterion, perform iterative classification calculations on pixels, and output the location information of crop planting areas and crop type information.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images