Fine identification method for crop types of cultivated land based on time-series remote sensing images

Abstract

The application discloses a farmland crop type fine identification method based on time sequence remote sensing images and relates to the technical field of farmland crop type identification. Time sequence remote sensing images and geological environment monitoring data of a whole growth cycle of a region to be identified are acquired, and a multi-modal heterogeneous input tensor is constructed. Through a double-flow channel of an orthogonal decoupling convolution coding network, crop biomass feature components and geological background feature components are separated and extracted. The components are taken as joint input to a differentiable logic reasoning identification model. A global utility function with the maximum classification confidence, feature orthogonality and physical constraint satisfaction degree as the target is defined, a gradient ascent method is used to iteratively optimize the logic node parameters of the model, and the optimal decision path distribution is output until convergence. In combination with the safety margin of the physical constraint boundary, a crop-environment adaptation analysis index is calculated, and high-precision and interpretable crop type identification is realized. The application improves the identification accuracy and decision transparency through feature decoupling and physical logic fusion.

Description

Technical Field

[0001] This invention relates to the field of arable land crop type identification technology, and more specifically, to a method for refined identification of arable land crop types based on time-series remote sensing imagery. Background Technology

[0002] In the field of agricultural remote sensing monitoring technology, the automated identification of crop types and their distribution using remote sensing image data is a crucial technical means to achieve digital and precision agriculture. Current mainstream technologies typically rely on multispectral time-series remote sensing image data, utilizing the differentiated spectral reflectance characteristics of crops at different growth stages, such as sowing, growth, and maturity, to construct time-series curves for classification. However, current technologies still face significant bottlenecks in practical applications. On the one hand, the spectral characteristics of crops in remote sensing images are often deeply intertwined with the characteristic signals of environmental factors such as soil background, topography, and water and fertilizer conditions. The phenomena of identical spectral features and dissimilar spectral features are extremely common, making it difficult for purely data-driven models to effectively extract the pure growth information of the crop itself from the mixed signals, thus limiting classification accuracy. On the other hand, the utilization of prior knowledge in crop phenology, soil science, and geography is severely insufficient. The model training process relies solely on statistical data patterns, and its decision-making logic is often a black box, easily leading to erroneous judgments that violate basic natural common sense, such as identifying rice on a dry hillside. Summary of the Invention

[0003] The purpose of this invention is to provide a refined identification method for crop types in cultivated land based on time-series remote sensing images. It separates crop and environmental features through an orthogonal decoupling network and innovatively incorporates physical logic constraints into model training. This not only improves the identification accuracy but also makes the decision-making process transparent, ultimately producing highly reliable identification results that are quantifiable and physically interpretable.

[0004] This invention is achieved through the following technical solution:

[0005] A method for refined identification of crop types in cultivated land based on time-series remote sensing imagery, comprising the following steps:

[0006] Acquire time-series remote sensing image sets and geological environment monitoring data of the area to be identified throughout its entire growth cycle, and construct a multimodal heterogeneous input tensor through spatiotemporal registration and dimensional mapping;

[0007] The multimodal heterogeneous input tensor is input into an orthogonal decoupled convolutional coding network, processed through a dual-stream convolutional channel, and crop biomass feature components and geological background feature components are extracted and output.

[0008] The crop biomass characteristic components and geological background characteristic components are used as joint inputs and injected into the differentiable logic reasoning recognition model, wherein the differentiable logic reasoning recognition model has built-in geological-phenological physical constraint boundaries.

[0009] Define a global utility maximization objective function, with the joint optimization objectives of maximizing classification confidence, maximizing feature component orthogonality, and maximizing physical constraint satisfaction. Perform gradient ascent iterative solution on the logic node parameters of the differentiable logic reasoning recognition model until the objective function converges and output the optimal decision path distribution.

[0010] Based on the optimal decision path distribution output by the differentiable logic reasoning identification model, and combined with the safety margin of the physical constraint boundary, a quantitative crop-environment adaptation analysis index is calculated to complete the refined identification of arable land crop types.

[0011] Optionally, the construction of the multimodal heterogeneous input tensor specifically involves:

[0012] Based on Gaussian process regression, interpolation is performed on geological environment monitoring data to generate a continuous geological grid field with spatial resolution consistent with the time-series remote sensing images of the entire growth cycle.

[0013] The geological grid field is replicated and expanded in the time dimension based on the time axis broadcasting mechanism, so that it is aligned with the frame number of the full growth cycle time-series remote sensing image set.

[0014] The aligned geological environment monitoring data is used as an additional channel and deeply stitched with the full-cycle time-series remote sensing image channel to form a dimension of [missing information]. Multimodal heterogeneous input tensors;

[0015] in, This represents the dimension of the multimodal heterogeneous input tensor. This indicates the number of spectral band channels in a remotely sensed image. This indicates the number of characteristic channels in geological environment monitoring data; Indicates the frame length of the time series; This represents the spatial height in pixels of the input image; This represents the spatial width of the input image in pixels.

[0016] Optionally, the orthogonal decoupled convolutional coding network includes two parallel encoder branches, namely: a biomass encoder. and background encoder These correspond to the output crop biomass characteristic components. and geological background feature components ;

[0017] In the objective function of maximizing global utility, an orthogonality metric is introduced. As the maximization sub-item, its calculation formula is:

[0018]

[0019] in, This represents the extracted dynamic feature vector of crop biomass; This represents the extracted static feature vector of the geological background; A smoothing constant is set to prevent the denominator from being zero.

[0020] Optionally, the differentiable logic reasoning recognition model is constructed using a probabilistic soft routing mechanism, specifically as follows:

[0021] The differentiable logic reasoning recognition model consists of multiple cascaded soft decision nodes, where each soft decision node performs routing probability calculations.

[0022] After decoupling Input decision node, node The output is defined as the routing probability. :

[0023]

[0024] in, This represents the Sigmoid activation function; Represents a node Learnable weight vectors; Represents a node The bias term;

[0025] For any leaf node Its final path probability Defined as:

[0026]

[0027] in, Represents the distance from the root node to the leaf node. The set of path nodes; As a direction indicator variable, when the path is to the left When the path turns right .

[0028] Optionally, the global utility maximization objective function The calculation formula is as follows:

[0029]

[0030] in, This represents the total optimization utility value; Indicates based on the final classification probability distribution The log-likelihood utility function; Indicates an orthogonality metric; Indicates physical and logical compliance metrics; The weight hyperparameters representing orthogonality constraints; The weight hyperparameters represent the physical logic constraints.

[0031] Optionally, the physical logic compliance indicators The specific calculation formula is as follows:

[0032]

[0033] in, Indicates the total number of physical constraints; This represents the input feature vector, i.e. and The combination; Indicates the first The theoretically feasible hyperplane defined by a physical constraint; Representing the eigenvector To the feasible hyperplane The Euclidean distance; This represents the bandwidth parameter that controls the width of the potential energy trap.

[0034] Optionally, the logic node parameters of the differentiable logic reasoning recognition model are solved using gradient ascent iterative methods, specifically as follows:

[0035] Constructing the Lagrange function :

[0036]

[0037] in, This represents the set of all learnable parameters in the network; Represent the set of Lagrange multipliers; Represents the total number of inequality constraints; Indicates the first A physical constraint function, when When indicates that the constraint is satisfied, This indicates a violation of the constraint;

[0038] During the iteration process, the following updates are performed synchronously:

[0039]

[0040]

[0041] in, The learning rate represents the network parameters; Indicates the ascending step size of the Lagrange multiplier; Indicates the parameter The gradient operator.

[0042] Optionally, the quantitative crop-environment adaptation resolution index The calculation formula is:

[0043]

[0044] in, This represents the final analytical index; Represents the set of all leaf nodes in a decision tree; This indicates that the sample falls into the first... The path probability of each leaf node; Indicates the first The category prediction confidence vector of each leaf node; This represents the normalized mapping function; This represents the geometric distance of the feature vector to the nearest physical forbidden boundary, used to characterize the safety margin.

[0045] Optionally, after completing the refined identification of crop types on arable land, the system further includes logical conflict backtracking and correction logic, which specifically includes:

[0046] When calculated Below the preset confidence threshold At that time, it was determined to be a "logically conflicting sample";

[0047] Activate the inverse gradient attention mechanism and calculate Gradient plot of the input tensor in the time dimension :

[0048] according to The peak positioning of the key time frame that causes logical conflict is masked or time-series interpolation repaired, and then the recognition is re-executed.

[0049] Optionally, the calculation Gradient plot of the input tensor in the time dimension The specific calculation formula is as follows:

[0050]

[0051] in, Indicates the input tensor in time frame Data slices.

[0052] The technical solution of the present invention has at least the following advantages and beneficial effects:

[0053] This invention, by constructing an orthogonal decoupled coding network, can forcibly separate biomass features characterizing crop growth dynamics from static background features characterizing the geological environment from mixed multimodal input data, fundamentally solving the feature confusion problem and obtaining a pure crop signal for identification. Furthermore, this invention creatively transforms agronomical prior knowledge such as crop phenological patterns, growth cycles, and geological environment compatibility into differentiable physical logic constraints, deeply integrating them into the model optimization process. This ensures that the model's inference results are not only statistically optimal but also physically logically sound, significantly improving the model's robustness and accuracy. In addition, this invention employs a differentiable logic inference model, making the decision path clear and transparent, and ultimately outputting a quantitative crop-environment fit resolution index. This index integrates statistical confidence and physical logic conformity, providing users with an unprecedented, physically meaningful measure of the reliability of classification results. Attached Figure Description

[0054] Figure 1 This is a flowchart illustrating the method for refined identification of crop types in cultivated land based on time-series remote sensing images provided by the present invention.

[0055] Figure 2 This is a schematic diagram illustrating the principle of the refined identification system for crop types in cultivated land based on time-series remote sensing images provided by the present invention. Detailed Implementation

[0056] 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. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0057] As the most important agricultural production resource, the accurate identification of crop types and growth conditions on arable land is a crucial data foundation for ensuring food security, optimizing agricultural subsidies, guiding agricultural production, and forecasting yields. Traditional crop identification methods heavily rely on manual field surveys, which are not only time-consuming and labor-intensive but also difficult to implement large-scale, high-time dynamic monitoring. With the development of remote sensing technology, crop classification using image data acquired by satellites or drones has become mainstream. However, existing technologies still face many challenges:

[0058] First, the phenomena of identical but different spectra and different spectra of identical objects are common. Different crops may exhibit extremely similar spectral characteristics at certain growth stages, while the spectral characteristics of the same crop may also show huge variations due to differences in planting density, water and fertilizer management, and soil background. In particular, the soil background, with its color, humidity, and mineral composition, greatly interferes with remote sensing signals, resulting in a deep overlap between crop characteristics and environmental background characteristics.

[0059] Second, there is insufficient utilization of temporal information and prior knowledge. Crop growth is a continuous and dynamic process, and its phenological rhythms (when to germinate, when to head, and when to mature) are the core basis for distinguishing different crops. Existing methods mostly rely on single-phase images, or the use of multi-phase images is relatively simplistic, failing to fully explore the phenological patterns contained in the time-series curves. At the same time, a great deal of agronomic knowledge and common geological knowledge, such as "rice is mostly planted in low-lying, clayey soil areas" and "the growth cycle of corn is usually around 120 days," are often ignored in existing purely data-driven models, leading to erroneous identification results that may violate basic physical principles.

[0060] Third, the "black box" nature of the model's decision-making process and the insufficient quantitative analytical capability of the results pose challenges. While deep learning models have achieved significant success in classification accuracy, their decision-making processes are often difficult to interpret, failing to explain to users why a particular plot of land is identified as a specific crop. Attribution analysis is also challenging when the model makes mistakes. Furthermore, traditional classification results typically provide only a category label, lacking a quantitative assessment of the reliability of the classification result. In other words, they cannot express the difference between "this plot is identified as wheat, and its characteristics closely match the standard wheat growth pattern" and "this plot is identified as wheat, but some of its characteristics are on the fuzzy edge."

[0061] Therefore, there is an urgent need to propose a new technical solution that can decouple the interference between crop characteristics and environmental background from a mechanistic perspective, deeply integrate the temporal phenological patterns and the physical constraints of the geological environment, and provide quantitative identification results with clear physical meaning and credibility.

[0062] Reference Figure 1 This embodiment details a method for refined identification of crop types in cultivated land based on time-series remote sensing imagery. This embodiment uses a major agricultural production area as the region to be identified, with the goal of refined identification of three main crops: wheat, corn, and soybeans.

[0063] Step 1: Acquisition, preprocessing and input tensor construction of multimodal heterogeneous data.

[0064] The goal of this step is to construct a high-dimensional, aligned, and information-complete data cube as input for subsequent deep learning models.

[0065] Temporal Remote Sensing Image Set: This embodiment selects the Sentinel-2 MSI sensor released by the European Space Agency (ESA) as the primary data source. This data has a spatial resolution of 10 meters and a revisit period of 5 days, effectively capturing key phenological changes in crops. The acquisition timeframe covers the complete growth cycle of the target area from March 2022 to October 2022, ensuring coverage of wheat's greening, heading, and maturity stages, as well as the sowing, growth, and maturity processes of summer maize and summer soybeans. The acquired data is L2A level product, which has undergone preliminary atmospheric correction. The downloaded image bands include blue, green, red, near-infrared, and three red-edge bands, totaling seven bands.

[0066] Geological environment monitoring data include:

[0067] Soil type data: derived from 1:500,000 scale soil map vector data generated from the Second National Soil Survey. This data contains classification information of soils within the study area, such as sandy soil, loam, and clay.

[0068] Soil physicochemical properties data: Soil samples were obtained from 100 sampling points by grid sampling of the study area, and their physicochemical indicators such as soil organic matter content, pH value, and soil permeability coefficient were measured in the laboratory.

[0069] Digital Elevation Model (DEM) Data: NASA-released SRTM (Shuttle Radar Topography Mission) data with a spatial resolution of 30 meters was used to extract topographic factors such as elevation, slope, and aspect of the land parcel.

[0070] Data preprocessing of the temporal remote sensing image set was performed as follows: First, cloud, cloud shadow, and snow cover masks were applied to all temporal images using the quality assessment bands provided by Sentinel-2 imagery, generating a high-quality cloud-free pixel observation sequence. Then, all bands were resampled to unify the spatial resolution to 10 meters. Next, commonly used vegetation indices such as the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EDI), and Normalized Difference Water Index (NDDI) were calculated as new feature channels. Finally, images from 20 valid temporal phases acquired within one growth cycle were time-stacked.

[0071] The geological environment data was rasterized and registered. Specifically, the vector-format soil map was rasterized to generate a soil type raster map with the same coordinate system and 10-meter resolution as the remote sensing image. Using the 100 soil physicochemical property sampling points obtained in step 1, spatial interpolation was performed using the Gaussian process regression method. Gaussian process regression is a non-parametric Bayesian method that can generate continuous predictive surfaces with uncertain estimates based on a small number of sampling points. This method generated a continuous raster field of soil organic matter, pH, and permeability coefficients at a 10-meter resolution covering the entire study area. The 30-meter resolution DEM data was resampled to 10 meters using bilinear interpolation, and slope and aspect were calculated.

[0072] Constructing multimodal heterogeneous input tensors is crucial for achieving multi-source data fusion. This embodiment fuses preprocessed static geological environment data with dynamic time-series remote sensing data.

[0073] This embodiment treats the static geological environment raster data of seven channels—soil type, organic matter, pH value, permeability coefficient, altitude, slope, and aspect—as a single dimension. The tensor of, where, and These represent the height and width in pixels of the study area, respectively. Using a time-axis broadcasting mechanism, this static tensor is copied and expanded 20 times along the time dimension, resulting in its dimension becoming... This aligns with the number of frames in the time-series remote sensing image set (20 in this embodiment).

[0074] The temporal tensor (dimension 1) of the preprocessed remote sensing image (containing 7 original bands and 3 vegetation indices, totaling 10 channels) is used to generate the tensor of the image. ), and perform deep splicing with the geological environment tensor after broadcast in the feature channel dimension.

[0075] Ultimately, a dimension is constructed. ,Right now The input tensor is a multimodal heterogeneous tensor. Each spatial location in this tensor contains a spectral time series from March to October, and each frame carries invariant geological background information for that location. This tensor serves as the sole data entry point for the subsequent recognition model.

[0076] Step 2: Feature extraction based on orthogonal decoupled convolutional coding networks.

[0077] The core objective of this step is to effectively separate the crop growth information and geological background information mixed in the input tensor, so as to obtain pure crop features for classification and environmental features for physical constraints.

[0078] This embodiment designs a dual-stream parallel convolutional encoder network. This network receives the multimodal heterogeneous input tensor constructed in step one, specifically including:

[0079] Biomass encoder This branch is specifically designed to extract dynamic temporal features related to crop growth rhythms. Its network structure employs a 3D convolutional neural network (3D-CNN). The 3D convolutional kernels can slide simultaneously across the temporal and spatial dimensions, thus directly learning dynamic patterns within the spatiotemporal cube. Specifically, the encoder consists of three consecutive 3D convolutional layers, a batch normalization layer, and a ReLU activation function layer, finally outputting a fixed-length feature vector through a global average pooling layer. . It is designed to characterize the intrinsic patterns of crop phenology from sowing to maturity.

[0080] Background Encoder This branch is specifically designed for extracting static geological background features. Its network structure employs a standard 2D convolutional neural network (2D-CNN), such as a simplified version of ResNet. Since the geological background information is constant over time, this encoder processes only the first time frame (or the average of all time frames) of the input tensor. The network consists of several 2D convolutional residual blocks and also outputs a fixed-length feature vector through a global average pooling layer. . It is designed to characterize the inherent properties of a plot of land, such as soil type, fertility, and topography.

[0081] To ensure and True independence, that is It does not contain geological information. The model does not include crop information, and orthogonality constraints are introduced into its optimization objective. Specifically, this is achieved by defining an orthogonality metric. To achieve this, the calculation formula is as follows:

[0082]

[0083] in, This represents the extracted dynamic feature vector of crop biomass; This represents the extracted static feature vector of the geological background; A smoothing constant is set to prevent the denominator from being zero.

[0084] In the specific application of this embodiment, the physical meaning of this formula is that the absolute value of the dot product of two vectors divided by the product of their respective lengths yields the absolute value of the cosine of the angle between them. When two vectors are completely orthogonal (perpendicular), the angle is 90 degrees, and the cosine value is 0. The cosine value reaches its maximum value of 1. When two vectors are parallel or collinear, the absolute value of the cosine is 1. The minimum value of 0 is reached. During model training, this is achieved by maximizing... The gradient descent algorithm will backpropagate to adjust... and The network weights are designed to force the feature representations they learn to be mutually perpendicular in the feature space. This mechanism fundamentally forces the network to decouple time-varying crop growth information from time-invariant geological background information, separating them into mutually orthogonal subspaces.

[0085] Step 3: Construction and optimization of the identification model based on differentiable logic reasoning.

[0086] This step replaces the traditional black-box classifier with an innovative, differentiable logical model that incorporates physical constraints.

[0087] In this embodiment, the construction process of the differentiable logic reasoning recognition model is as follows:

[0088] The differentiable logic reasoning recognition model borrows the idea of ​​decision trees, but makes its decision-making process fully differentiable by "softening" it, so that it can be jointly trained end-to-end with the aforementioned encoding network. The model consists of multiple cascaded soft decision nodes.

[0089] Traditional decision trees partition samples into left or right subtrees at each node based on a hard threshold (e.g., "NDVI > 0.5"). The soft decision node in this embodiment differs; it does not perform hard partitioning but instead calculates a routing probability. This decouples the previous step to obtain a pure crop biomass feature vector. Input is fed into the decision node. For any internal node in the model... Its output is defined as the routing probability. :

[0090]

[0091] in, This represents the Sigmoid activation function, whose output ranges between (0,1); Represents a node The learnable weight vector, whose dimension is... same; Represents a node The bias term; The calculated result, for example, 0.8, is interpreted as meaning that the current sample has an 80% probability of belonging to the logic of the left branch and a 20% probability of belonging to the logic of the right branch. The information flow of the sample will be weighted according to this probability and will flow to both the left and right subtrees simultaneously, rather than being an either-or choice.

[0092] After passing through multiple layers of soft routing, the sample information flow eventually reaches the leaf nodes of the model. Each leaf node stores the probability distribution of the categories. For any leaf node... The total probability that a sample eventually "reaches" that leaf node. It is the product of the routing probabilities of all internal nodes traversed along the path from the root node to the leaf node. For example, the path probability of a path that goes left first and then right is... Ultimately, the output of the entire model is a weighted average of the probabilities of all leaf nodes, forming the final classification probability distribution for the three categories of wheat, corn, and soybeans. .

[0093] The global utility maximization objective function defined in this embodiment is the core driving force for the entire model training. This embodiment constructs an objective that maximizes the utility function, more intuitively reflecting the process of finding the optimal solution. Objective function It consists of three weighted parts:

[0094]

[0095] in, This represents the total optimization utility value; Indicates based on the final classification probability distribution The log-likelihood utility function; Indicates an orthogonality metric; Indicates physical and logical compliance metrics; The weight hyperparameter representing the orthogonality constraint is set to [value] in this embodiment. ; This represents the weight hyperparameter of the physical logic constraints. In this embodiment, it is set to... .

[0096] log-likelihood effect Used to measure the accuracy of model predictions. It is calculated as the inverse of the classification cross-entropy. Specifically, for a training sample with the true label "wheat," if the model outputs a probability distribution... In this model, if the probability of wheat is 0.9, then the utility value of this term is high; conversely, if the probability of wheat is low, then the utility value is also low. Maximizing this term drives the model to learn how to classify correctly.

[0097] Orthogonality metrics This refers to the orthogonality index defined in step two. Incorporating it into the total utility function ensures that the model, while pursuing classification accuracy, must also consider the task of feature decoupling.

[0098] Physical and logical compliance indicators This is the key innovation that distinguishes this invention from purely data-driven models. This metric is used to measure whether the model's predictions conform to predefined agronomic and geological common sense. Its construction employs a potential energy trap function form:

[0099]

[0100] in, Indicates the total number of physical constraints; This represents the input feature vector, i.e. and The combination; Indicates the first The theoretically feasible hyperplane defined by a physical constraint; Representing the eigenvector To the feasible hyperplane The Euclidean distance; The bandwidth parameter represents the control of the potential energy trap width; in this embodiment, it is set to 1.5.

[0101] In this embodiment, when the feature vector When it falls exactly in the center of the physically feasible region Approaching 0 It approaches its maximum value of 1. As the feature vector deviates further and further from the feasible region,... It will rapidly decrease to 0 in a Gaussian decay manner. Maximizing this term means that the model must push the feature representation of the samples into a region that conforms to physical laws.

[0102] In this embodiment, the following are defined: A specific physical constraint:

[0103] Constraint 1 is a geological matching constraint: Wheat and corn are dryland crops, requiring a certain level of soil permeability and unsuitable for growth in clay soils prone to prolonged flooding. Soybeans, on the other hand, are more adaptable to various soil types. Therefore, a constraint is defined: if the model predicts wheat or corn, then the corresponding geological background feature components... The decoded soil permeability coefficient must be greater than cm / s. A prediction that violates this rule will result in... Increase reduce.

[0104] Constraint 2 is a phenological monotonicity constraint: the NDVI value of wheat should show a monotonically increasing trend from the time it turns green in spring until it heads. Therefore, a constraint is defined: if the model predicts wheat, then its corresponding crop biomass characteristic components... The decoded time-series NDVI curve must have a positive slope in any sub-interval during April and May. An abnormal decline violates this rule. Increase reduce.

[0105] Constraint 3 is a growth cycle constraint: the growth cycle of summer maize is typically 90-120 days. Therefore, a constraint is defined: if the model's prediction is maize, then its corresponding... The complete cycle length of the decoded NDVI time series curve, from trough to peak and back to trough, must be within a time window of 80-130 days. Significantly exceeding or falling short of this range constitutes a violation of this rule. Increase reduce.

[0106] In this embodiment, the optimization process of the model is to find the optimal solution. The largest network parameters (including the weights of the two encoders and the weights of all soft decision nodes) With bias The process employs an optimization algorithm based on the dual ascent strategy of the augmented Lagrange multiplier method.

[0107] Specifically, this algorithm can be simply understood as a penalized hill-climbing process. The aforementioned physical constraints (such as a permeability coefficient greater than a threshold) can be expressed as inequalities. Construct the Lagrange function in the form of... ,in, It is the collection of all network parameters. These are Lagrange multipliers. In each iteration, the model parameters... Along make Update a small step in the direction of increase (gradient ascent). Simultaneously, check if the current sample violates any physical constraints. If violated Then the corresponding Lagrange multipliers Let's increase it a little. This increase... In the next iteration, the penalty for violating this constraint will be increased, acting like a dynamically adjusted "fine," forcing the model parameters to be updated back to the physically feasible region. In this way, classification accuracy is maximized while strictly adhering to all pre-defined physical logic.

[0108] Step 4: Quantitative analysis index calculation and refined identification are completed.

[0109] After the model is trained, for a new land parcel to be identified, the final identification result is no longer a simple category label, but a more informative quantitative index.

[0110] In this embodiment, the crop-environment adaptability index aims to quantify the comprehensive reliability of the identification results at both statistical and physical-logical levels. Its calculation formula is as follows:

[0111]

[0112] In this specific application, the first term on the right side of the equals sign is the purely statistical classification confidence score output by the model. Wherein, Represents the set of all leaf nodes in a decision tree; This indicates that the sample falls into the first... The path probability of each leaf node; Indicates the first The category prediction confidence vector for each leaf node (a vector representing the probability of wheat, corn, and soybean) is calculated as the weighted average probability of the plot belonging to each category.

[0113] The second term on the right side of the equation is the physical safety margin coefficient. Among them, It is the feature vector of the land parcel ( and This refers to the geometric distance from the nearest physical prohibition boundary. For example, if a plot of land is identified as wheat, we would calculate how far its soil permeability characteristics are from the prohibition boundary for "too low permeability". The greater this distance, the closer it is to the center of the safe zone. It is a normalization function that maps this distance to the interval (0,1). The larger the distance, the better. The closer the value is to 1.

[0114] final, It is the product of these two items. A higher one. A value (e.g., 0.95) means that the plot is not only statistically highly probable to contain a particular crop (e.g., wheat), but its biophysical characteristics (time series curves, soil matching, etc.) are also perfectly located in the central region of standard wheat, making it a very "standard" and "healthy" sample, and the identification result is highly reliable. A medium-sized... A value (e.g., 0.6) may mean that although it is statistically likely to be wheat, a certain physical characteristic, such as soil pH, is close to the tolerance boundary of wheat, making it a "marginal" sample, and the result is subject to some uncertainty.

[0115] Furthermore, this embodiment also includes a post-processing self-checking and correction mechanism. A confidence threshold is set, for example... When the highest calculated value of a certain plot of land... If the value is still below this threshold, the system will identify it as a logical conflict sample. The typical characteristic of such a sample is that its spectral time series characteristics are statistically very similar to crop A (e.g., corn), but its geological environment forcibly restricts it from being crop A (e.g., located in a high-altitude, cold mountainous area, which does not meet the accumulated temperature requirements for corn growth).

[0116] At this point, the system activates the inverse gradient attention mechanism. Specifically, it calculates the physical logic compliance index. The gradient map is obtained by taking the absolute value of the partial derivative of the original multimodal input tensor with respect to each frame in the time dimension. .

[0117]

[0118] in, Indicates the input tensor in time frame Data slices. Gradient values. The larger the time frame, the greater the "contribution" of the image data in that frame to the physical logic conflict. In practice, this usually corresponds to thin clouds that were not perfectly removed, cloud shadow pollution, or abnormal agricultural activities (such as irrigation or harvesting outside of the normal season). After locating this critical time frame, the system automatically performs masking processing on it and uses linear interpolation with the image data from the two preceding and following time frames to generate a repaired "reasonable" data. Then, the repaired data is re-inputted into the already trained model, and the recognition process is executed again. This backtracking correction process significantly improves the model's robustness to complex and dirty data and the final recognition accuracy.

[0119] Reference Figure 2 As shown, a refined crop identification system for cultivated land implementing the above method can have its hardware deployed on a high-performance computing server cluster or cloud platform. Logically, the system includes:

[0120] Heterogeneous Data Stream Tensorization Module: This module is responsible for establishing interfaces with various remote sensing data archives and geographic information databases, automatically pulling and preprocessing data. Its core is a GPU-based parallel computing subsystem used to efficiently perform spatial interpolation, rasterization, and time-axis broadcasting operations on geological data, ultimately constructing a multimodal heterogeneous input tensor that meets the required format.

[0121] Orthogonal Feature Extraction Engine: This module solidifies the two-stream convolutional network structure described in step two. During the inference phase, it receives the output of the tensor quantization module and, through one forward propagation, directly calculates and outputs mutually orthogonal crop biomass feature components. Geological background characteristic components .

[0122] Differentiable logic solver: During the training phase, this module is the core optimizer, with a built-in automatic differentiation framework responsible for constructing the global utility maximization objective function. It then executes a complex Lagrange dual ascent algorithm to update the parameters of the entire network. During the inference phase, it receives feature components, performs probabilistic routing through a trained soft decision tree, and outputs the final classification probability distribution.

[0123] Resolution Index Calculation Unit: This module is the final output of the system. It receives the classification probability from the solver and the feature components from the feature extraction engine, and, combined with a pre-set physical constraint boundary database, calculates the crop-environment adaptation resolution index for each plot in real time. Furthermore, it includes control logic for logical conflict detection and backtracking correction. Ultimately, this unit generates a visualized crop classification map, along with corresponding... Credibility graphs provide intuitive and reliable data support for agricultural management decisions.

[0124] In summary, this embodiment constructs a multimodal heterogeneous input, designs an orthogonal decoupling network to separate crop and environmental features, and innovatively integrates physical constraints into a differentiable logic model for maximizing utility optimization. Finally, it uses a quantitative analytical index to measure the reliability of the results, forming a complete, clear, feasible, and significantly superior refined identification scheme for cultivated land crops compared to existing technologies.

[0125] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for fine identification of crop types in cultivated land based on time-series remote sensing images, characterized in that, The steps of the method comprise: Obtain the full growth cycle time series remote sensing image set and the geological environment monitoring data of the to-be-identified region, construct a multi-modal heterogeneous input tensor through space-time registration and dimension mapping; Input the multi-modal heterogeneous input tensor into an orthogonal decoupling convolutional coding network, process through a double-flow convolutional channel, extract and output crop biomass feature components and geological background feature components; Take the crop biomass feature components and the geological background feature components as joint inputs, inject them into a differentiable logic reasoning identification model, wherein the differentiable logic reasoning identification model is built-in with a geological-phenology physical constraint boundary; Define a global utility maximization objective function to maximize the classification confidence, maximize the orthogonality of the feature components, and maximize the satisfaction degree of the physical constraint as the joint optimization goal, perform gradient ascent iterative solution on the logic node parameters of the differentiable logic reasoning identification model until the objective function converges, and output the optimal decision path distribution; Based on the optimal decision path distribution output by the differentiable logic reasoning identification model, combined with the safety margin of the physical constraint boundary, a quantitative crop-environment adaptation analysis index is calculated to complete the fine identification of the cultivated crop type; The orthogonal decoupled convolutional encoding network comprises two parallel encoder branches, which are respectively a biomass encoder and a background encoder , corresponding to output crop biomass feature components and geological background feature components respectively; In the global utility maximization objective function, the orthogonality measure index is introduced As the maximization sub-item, its calculation formula is: wherein, denotes the extracted crop biomass dynamic feature vector; denotes the extracted geology background static feature vector; a smoothing constant set to prevent the denominator from being zero; The differentiable logic reasoning identification model is constructed by using a probability soft routing mechanism, specifically: The differentiable logic reasoning identification model is composed of multiple cascaded soft decision nodes, wherein each soft decision node performs routing probability calculation; Decoupled Input decision node, node Output of the node is defined as routing probability : wherein, denotes a Sigmoid activation function; denotes a learnable weight vector for a node ; denotes a bias term for a node ; For any leaf node its final path probability is defined as: wherein, represents a set of path nodes from the root node to the leaf node ; is a direction indication variable, when the path goes left , when the path goes right . 2.The method of claim 1, wherein the method further comprises: determining a time interval between the first time point and the second time point; and determining a time interval between the second time point and the third time point. The construction of the multi-modal heterogeneous input tensor is specifically: Interpolate the geological environment monitoring data based on Gaussian process regression to generate a continuous geological grid field consistent with the spatial resolution of the full growth cycle time series remote sensing image; Based on the time axis broadcast mechanism, the geological grid field is replicated and expanded in the time dimension to align with the frame number of the full growth cycle time series remote sensing image set; The aligned geological environment monitoring data is taken as an additional channel, and is deeply spliced with a whole growth cycle time sequence remote sensing image channel to form a multi-modal heterogeneous input tensor with a dimension of ​ wherein, denotes the dimension of the multi-modal heterogeneous input tensor, denotes the spectral band channel number of the remote sensing image; denotes the feature channel number of the geological environment monitoring data; denotes the frame length of the time series; denotes the spatial height pixel number of the input image; denotes the spatial width pixel number of the input image. 3.The method of claim 2, wherein, The global utility maximization objective function The formula is: wherein, represents the total optimization utility value; represents a log-likelihood utility function based on the final classification probability distribution represents an orthogonality measure indicator; represents a physical logic compliance measure indicator; represents a weight hyperparameter for the orthogonality constraint; represents a weight hyperparameter for the physical logic constraint.​ 4. The method according to claim 3, wherein, The physical logic compliance indicator The specific calculation formula is: wherein, denotes the total number of physical constraints; denotes the input feature vector, i.e. in combination with ; denotes the theoretical feasible region hyperplane defined by the th physical constraint; denotes the Euclidean distance of the feature vector to the feasible region hyperplane ; denotes the bandwidth parameter controlling the potential well width.

5. The method according to claim 4, wherein, The gradient ascent iterative solution of the logic node parameters of the differentiable logic reasoning identification model is specifically: Constructing a Lagrangian function : wherein, denotes the set of all learnable parameters in the network; denotes the set of Lagrange multipliers; denotes the total number of inequality constraints; denotes the th physical constraint function, denotes satisfaction of the constraint when denotes violation of the constraint when denotes violation of the constraint when During the iteration process, the following updates are performed synchronously: wherein, denotes a learning rate of a network parameter; denotes an ascent step size of a Lagrange multiplier; denotes a gradient operator on a parameter . 6.The method of claim 5, wherein the method further comprises: The quantified crop-environment adaptation analysis index The calculation formula is: wherein, represents the final resolution index; represents the set of all leaf nodes of the decision tree; represents the path probability of the sample falling into the th leaf node; represents the class prediction confidence vector of the th leaf node; represents the normalization mapping function; represents the geometric distance of the feature vector to the nearest physical forbidden boundary, used to characterize the safety margin.

7. The method according to claim 6, wherein, After completing the fine identification of the cultivated crop type, it further includes a logic conflict backtracking correction logic, which is specifically: When the calculated is lower than a preset trust threshold , it is determined as a "logic conflict sample"; activate the reverse gradient attention mechanism, compute gradient graph on the input tensor time dimension : According to the peak positioning leads to a critical time frame of logical conflicts, which is masked or time-interpolated repaired, and the recognition is re-executed. 8.The method of claim 7, wherein, The calculation Gradient graph for input tensor time dimension The specific calculation formula is: wherein, represents a data slice of the input tensor at a time frame t.

Images