Farmland irrigation process remote sensing dynamic monitoring method supported by spatio-temporal data fusion

Abstract

The invention discloses a farmland irrigation process remote sensing dynamic monitoring method supported by spatio-temporal data fusion, which comprises the following steps: acquiring a high time-low spatial resolution image sequence and a low time-high spatial resolution image sequence which are unified by spatio-temporal reference in a monitoring period, and calculating a fusion weight based on spatio-temporal correlation of the two image sequences, generating a corresponding day-by-day high-temporal-spatial-resolution fusion image sequence based on the day-by-day low-spatial-resolution image; calculating a spectral index sensitive to moisture change based on the fused image sequence, constructing a time sequence characteristic set, slicing the time sequence characteristic set by using a sliding window, inputting a pre-constructed irrigation identification model, analyzing time sequence evolution characteristics in the window, judging whether an irrigation event occurs at a space-time position corresponding to the window, and obtaining a dynamic distribution diagram of the irrigation area. According to the method, the problems that the irrigation abrupt change signal is smoothed and the identification precision is interfered by fusion noise are effectively solved, and fine dynamic monitoring of the irrigation event of the fragmented farmland area is realized.

Description

Technical Field

[0001] This invention relates to the field of agricultural monitoring technology, and in particular to a remote sensing dynamic monitoring method for farmland irrigation process supported by spatiotemporal data fusion. Background Technology

[0002] Dynamic monitoring of farmland irrigation events is a crucial foundation for achieving refined management of agricultural water resources. Because irrigation events are typically short-lived and spatially fragmented, a single remote sensing data source cannot simultaneously meet the dual requirements of high temporal resolution (capturing diurnal water changes) and high spatial resolution (distinguishing between small plots). Utilizing multi-source remote sensing data for spatiotemporal fusion to construct continuous image sequences with high spatiotemporal resolution, and combining this with time-series analysis techniques to identify irrigation behavior, is currently one of the mainstream technical approaches in this field.

[0003] Current mainstream technical solutions typically employ spatiotemporal fusion algorithms, such as ESTARFM (Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model), which utilizes temporal variation information from low spatial resolution images like MODIS (Moderate Resolution Imaging Spectroradiometer) to predict missing temporal data from high spatial resolution images like Sentinel-2, generating synthetic image sequences. Deep learning models, such as Bi-LSTM (Bidirectional Long Short-Term Memory), are then used to extract temporal features of spectral indices sensitive to water conditions to determine irrigation events. Under ideal conditions, these methods can effectively fill data gaps and restore the overall trend of crop growing season changes mapped in temporal image data.

[0004] However, existing techniques still have shortcomings in handling abrupt changes in irrigation events, mainly manifested in the smoothing effect of spatiotemporal fusion on abrupt signals and the sensitivity of recognition models to fusion noise. Traditional fusion algorithms are usually based on linear variation or spectral similarity assumptions, tending to treat rapid changes in reflectance caused by non-rainfalling factors, such as sudden changes in wetness caused by irrigation, as observation noise or registration errors and smoothing them out, resulting in the neglect of signals characterizing irrigation features. Affected by cloud shadows and other factors, the fused image sequence inevitably contains local quality uncertainties, and existing irrigation recognition models lack a mechanism to perceive the quality of input data, treating low-confidence interpolation errors the same as real surface signals, leading to misjudgments under complex meteorological conditions. Summary of the Invention

[0005] The purpose of this invention is to provide a remote sensing dynamic monitoring method for farmland irrigation process supported by spatiotemporal data fusion, in order to solve at least one of the aforementioned problems existing in the prior art.

[0006] According to one aspect of this application, a remote sensing dynamic monitoring method for farmland irrigation process supported by spatiotemporal data fusion includes:

[0007] Acquire high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences of the target area during the monitoring period, and unify the spatiotemporal reference of the two;

[0008] The fusion weight is calculated based on the spatiotemporal correlation of high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences, and a daily high spatiotemporal resolution fused image sequence covering the monitoring period is generated accordingly.

[0009] Based on the calculation of spectral indices sensitive to moisture changes using daily high spatiotemporal resolution fused image sequences, a temporal feature set is constructed.

[0010] The time-series feature set is sliced ​​based on a sliding window. The sliced ​​time-series features are then input into a pre-built irrigation identification model. The time-series evolution features within the window are analyzed to determine whether an irrigation event has occurred at the spatiotemporal location corresponding to the window, thus obtaining a dynamic distribution map of the irrigation area.

[0011] Beneficial effects: This invention effectively solves the problems of smoothing irrigation mutation signals and reducing noise interference in the identification accuracy of traditional methods, and realizes fine dynamic monitoring of irrigation events in fragmented farmland areas. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the overall process of a remote sensing dynamic monitoring method for farmland irrigation process supported by spatiotemporal data fusion, provided in an embodiment of this application.

[0013] Figure 2 This is a schematic diagram of the image preprocessing process provided in the embodiments of this application.

[0014] Figure 3 This is a schematic diagram of the irrigation identification process using a bidirectional long short-term memory network provided in the embodiments of this application. Detailed Implementation

[0015] Example 1 details the overall process of a remote sensing dynamic monitoring method for farmland irrigation process supported by spatiotemporal data fusion, such as... Figure 1As shown, this invention addresses the technical challenges of existing technologies where a single remote sensing data source cannot simultaneously meet the requirements for high spatial and temporal resolution, and where traditional monitoring methods struggle to accurately capture dynamic irrigation events. By combining multi-source data fusion with temporal deep learning pattern recognition, it achieves refined and dynamic monitoring of farmland irrigation activities.

[0016] The system environment used in this embodiment can be a local workstation based on a high-performance computing cluster or a cloud computing platform. Taking a certain irrigation district as an example, this area is mainly used for rice cultivation, with a high degree of land fragmentation, and irrigation behavior is affected by both human management and rainfall. It has typical characteristics of southern irrigation districts and is suitable as an application scenario for this method.

[0017] Step 101: Obtain high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences of the target area during the monitoring period, and unify the spatiotemporal reference of the high temporal-low spatial resolution image sequences and the low temporal-high spatial resolution image sequences.

[0018] In this embodiment, the target area refers to the geographic space to be monitored for irrigation, usually defined by regional divisions or irrigation district boundaries. The monitoring period refers to a continuous time period covering the main crop growth cycle, such as June to November during the rice growing season. High temporal-low spatial resolution image sequences refer to remote sensing image data with short revisit periods, such as daily revisits, but low spatial resolution, such as 250m, 500m, or 1000m. This type of data can capture rapid changes on the Earth's surface at the daily scale, but it is difficult to distinguish small plots. Common low spatial resolution image data sources include MODIS series products. Low temporal-high spatial resolution image sequences refer to remote sensing image data with high spatial resolution but limited by satellite revisit periods or the influence of cloud and rain weather, resulting in discontinuous temporal distribution and long intervals. This type of data can better show some plot boundaries and texture details, but it is difficult to provide continuous temporal information. Common high spatial resolution image data sources include Sentinel-2 or GF-1 / 2 / 6, etc.

[0019] Specifically, the unification of spatiotemporal references is a prerequisite for data fusion. This process includes unifying geographic coordinate systems, transforming spatial projections, and registering pixel locations. For example, all image data is projected uniformly onto a UTM (Universal Transverse Mercator) coordinate system, such as UTM Zone 50N, ensuring that low-spatial-resolution image pixels and high-spatial-resolution image pixels have a clear spatial nesting or correspondence. Furthermore, the acquisition dates of the images are documented to ensure a unified temporal reference, facilitating subsequent time-series analysis. Through this step, a spatiotemporal reference for multi-source data collaboration is constructed, laying the foundation for subsequent data fusion.

[0020] Step 102: Calculate the fusion weight based on the spatiotemporal correlation of high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences, and use the fusion weight to generate daily high spatiotemporal resolution fused image sequences covering the monitoring period.

[0021] This step utilizes the temporal variation information provided by low-spatial-resolution (LSRP) imagery to fill in the gaps in key temporal phases of high-SRP imagery. Simultaneously, it leverages the spatial texture information provided by high-SRP imagery to enhance the spatial details of the low-SRP imagery, generating a fused image that combines high spatial and temporal resolution. Spatiotemporal correlation refers to the inherent consistency or correlation of spectral variation patterns of ground features at different times and spatial locations. For example, crop phenological information within the same plot should exhibit relatively similar trends in both low-SRP and high-SRP imagery.

[0022] In practice, this step typically involves establishing a mapping relationship between the known time (reference time) of the high spatial resolution image and the time to be predicted (target time). It requires searching for neighboring pixels with spectral similarity to the target pixel within a local window; calculating the fusion weights of similar pixels based on factors such as spectral similarity, spatial distance, and time interval. A larger weight indicates a greater contribution of the similar pixel to the predicted value of the target pixel. The calculated weights are then used to weight and sum the reflectance changes of the similar pixels to deduce the high spatial resolution reflectance at the target time. This process can be based on the ESTARFM algorithm or the improved uncertainty-enhanced spatiotemporal adaptive reflectance fusion model IU-ESTARFM algorithm described in subsequent embodiments of this invention. The generated daily high spatiotemporal resolution fused image sequence will serve as the basic data source for subsequent feature extraction.

[0023] Step 103: Calculate spectral indices sensitive to moisture changes based on daily high spatiotemporal resolution fused image sequences, and construct a temporal feature set.

[0024] Spectral indices are indicators constructed through mathematical operations between spectral bands that can enhance predetermined land cover features, such as vegetation, water bodies, and soil moisture. In this embodiment, the focus is on spectral indices sensitive to changes in moisture, as the most direct physical manifestation of irrigation activity is a sharp increase in soil moisture content or the appearance of a water layer on the field surface. The temporal feature set refers to the numerical sequence of the aforementioned spectral indices over time during the monitoring period.

[0025] Specifically, the Revised Vertical Drought Index (RPDI) and the Normalized Difference Water Index (NDWI) can be calculated. Taking RPDI as an example, soil moisture is retrieved by utilizing the distribution characteristics of red and near-infrared bands near the soil line, demonstrating high sensitivity to soil moisture changes. For each pixel in the fused image sequence, its corresponding spectral index value is calculated daily, forming a continuous time-series curve. This curve records the moisture change trajectory of that pixel throughout the growing season, incorporating the superposition of multiple signals such as crop growth, natural rainfall, and artificial irrigation. Constructing a time-series feature set aims to transform high-dimensional image data into one-dimensional signals that more directly reflect physical processes, reducing the complexity of subsequent recognition models and improving the signal-to-noise ratio of the features.

[0026] Step 104: Slice the temporal feature set based on a sliding window, input the sliced ​​temporal features into a pre-built irrigation identification model, analyze the temporal evolution features within the window, determine whether an irrigation event has occurred at the spatiotemporal location corresponding to the window, and obtain a dynamic distribution map of the irrigation area.

[0027] In this step, the sliding window technique is used to capture the dynamic characteristics of irrigation events within a local timeframe. Since irrigation is typically a short process lasting from several hours to several days, and the irrigation time varies across different plots, a simple global classification of the entire time series curve is not feasible. The sliding window, by setting a fixed window length and step size (e.g., 10, 15, or 20 days, and a step size of 1 day), slides across the time series curve day by day, extracting segments of the local sequence, i.e., slices.

[0028] The irrigation identification model is a pre-trained algorithm model. Its input is the aforementioned time-series slice, and its output is the probability or category label of irrigation occurring at the center or end of the window. This model can be a rule-based threshold judgment model, or a machine learning-based model such as random forest or support vector machine, or a deep learning classification model such as Long Short-Term Memory (LSTM) or Bi-LSTM. The model's task is to analyze the temporal evolution characteristics within the window, especially identifying irrigation-specific signal patterns, such as the V-shaped feature of a sharp drop in spectral index within a short period (dry crops: becoming wet; rice and other crops: water layer in the field) followed by a slow recovery (dry crops: becoming dry; rice and other crops: no water layer in the field), and distinguishing it from the wetting characteristics caused by natural rainfall. By classifying all pixels across the entire region and time period window by window, all spatiotemporal points identified as irrigated can be mapped back to geographic space, generating a daily dynamic distribution map of irrigated area. This distribution map not only reflects the spatial location of irrigation but also records the specific time of irrigation, providing refined dynamic monitoring data for water resource management.

[0029] Example 2 describes data preprocessing and temporal feature construction, detailing the specific selection of multi-source remote sensing data, band alignment, cloud removal and restoration, and the calculation method of the core spectral index RPDI. It also supplements the preferred implementation methods for auxiliary index calculation and temporal smoothing processing, addressing the technical problems of spatiotemporal mismatch in multi-source data and the inability of a single spectral feature to comprehensively characterize the irrigation process. The image preprocessing workflow is as follows: Figure 2 As shown.

[0030] Step 201: Obtain MODIS surface reflectance products as a high temporal-low spatial resolution image sequence, and obtain Sentinel-2 images as a low temporal-high spatial resolution image sequence.

[0031] In this embodiment, for remote sensing image sequences with high temporal resolution and low spatial resolution, MODIS' MCD43A4 V6.1 surface reflectance product is preferred. This product, based on multi-day observation data, uses a two-way reflectance distribution function (BRDF) model to normalize surface reflectance, generating reflectance equivalent to that under zenith observation conditions, thus effectively eliminating the influence of differences in observation angles. This product provides data output daily with a spatial resolution of 500 meters, enabling relatively continuous characterization of the dynamic changes in surface reflectance characteristics. Data acquisition can be performed in batches through the GEE cloud platform (a cloud-based geospatial analysis platform), with filtering criteria including the geometric boundaries of the study area and the time range from May 25th to November 5th, 2025, resulting in 165 daily images.

[0032] In this embodiment, for low temporal-high spatial resolution image sequences, Level-2A data from Sentinel-2 MSI (Sentinel-2 MultiSpectral Image) is preferred. Level-2A data has undergone radiometric calibration and atmospheric correction, providing surface reflectance values ​​that can be directly used for quantitative remote sensing analysis. This data has a spatial resolution of up to 10 meters, a single-satellite revisit period of approximately 10 days, and a two-satellite revisit period of approximately 5 days. During acquisition, in addition to limiting the temporal and spatial range, a cloud cover threshold must be set. Preferably, the cloud cover is set to be less than or equal to 10% to ensure the acquisition of high-quality cloudless or low-cloud imagery as a high-resolution benchmark for spatiotemporal fusion.

[0033] Step 202: Select the band intersection of the high temporal-low spatial resolution image sequence and the low temporal-high spatial resolution image sequence, and resample the high temporal-low spatial resolution image sequence to the same spatial resolution as the low temporal-high spatial resolution image sequence.

[0034] In this embodiment, due to the differences in sensor band settings between MODIS and Sentinel-2, common bands are extracted to ensure the consistency of physical meaning in subsequent spatiotemporal fusion. Specifically, six bands common to both types of images are selected: Blue, Green, Red, NIR, SWIR1, and SWIR2. In Sentinel-2, the center wavelengths corresponding to these bands are approximately 490nm, 560nm, 665nm, 842nm, 1610nm, and 2190nm, respectively.

[0035] Specifically, to achieve pixel alignment, both types of images are first projected onto the same coordinate system, such as UTM zone 50N. For the MCD43A4 V6.1 image with a spatial resolution of 500 meters and the non-10-meter image bands in Sentinel-2, bilinear interpolation or the nearest neighbor method is used to resample them to 10 meters, and the pixel grids of the MCD43A4 V6.1 image and the Sentinel-2 image are strictly aligned. Although this step does not increase the actual information content of the low spatial resolution image, it constructs a unified computational grid, so that each 10-meter pixel corresponds to a low spatial resolution reflectance value, providing a data foundation for subsequently establishing the high-resolution-low-resolution mapping relationship.

[0036] Step 203: Perform cloud masking on the low temporal-high spatial resolution image sequence, and use a spatiotemporal interpolation method based on the combination of K-nearest neighbor and inverse distance weights to fill in the invalid value regions covered by the mask.

[0037] In this embodiment, although a cloud cover threshold is set during the data acquisition phase, local thin clouds or cloud shadows may still exist in the selected Sentinel-2 images. Cloud masks generated using the QA60 or SCL bands provided with the Sentinel-2 data can identify the pixel locations of clouds and cloud shadows. The reflectance values ​​at these locations are then marked as invalid.

[0038] Specifically, to ensure the integrity of the reference image, these invalid value regions must be repaired. This embodiment preferably employs a weighted interpolation method combining K-nearest neighbor (KNN) and inverse distance weighted (IDW). K valid approximate pixels are searched within the spatial neighborhood of the invalid pixel; for example, K can be 5 to 10. The Euclidean distance between these valid pixels and the target invalid pixel is calculated. The reciprocal of the distance is used as a weight to perform a weighted average of the reflectance of the valid pixels, filling in the invalid values. For cases where large-area cloud cover leads to poor spatial interpolation results, temporal co-interpolation can be performed by combining nearby images, such as pixels at the same location in adjacent time phases, to ensure that each Sentinel-2 image used as input is spatially continuous within the study area.

[0039] Step 204, the spectral indices include the revised vertical drought index (RPDI), and the construction of the time-series feature set is specifically based on the calculation of the daily RPDI sequence using the following formula:

[0040] ;

[0041] Where RPDI is the Recorrected Vertical Aridity Index, Red is the red band reflectance, and m t The soil line slope is corrected based on multi-temporal data, i.e., the soil line slope obtained by linearly fitting the red light and near-infrared reflectance of bare soil pixels in the current temporal image. NIR represents near-infrared reflectance, and F... v Red is the vegetation cover calculated based on the normalized vegetation index. v NIR represents the reflectance of vegetation in the red light band, specifically the reflectance of pixels with total vegetation cover in the red light band. v This represents the reflectance of vegetation in the near-infrared band, specifically the reflectance of pixels with full vegetation cover in the near-infrared band.

[0042] In this embodiment, the RPDI index is designed to eliminate the interference of vegetation background in mixed pixels on soil moisture retrieval. The parameter m... t This was obtained by constructing a Red-NIR spectral feature space. Specifically, bare soil pixel samples from the study area were selected, and their Red and NIR reflectances were linearly fitted. The slope of the fitted line is the soil line slope M, while m... t This is the value of M after correcting for multi-temporal variations. Parameter F v The vegetation cover within a pixel can be estimated using the Normalized Difference Vegetation Index (NDVI) and a pixel-based dichotomy model, i.e.:

[0043] F v =(NDVI-NDVI soil ) / (NDVI veg -NDVI soil );

[0044] Among them, NDVI soil and NDVI veg The NDVI thresholds are for bare soil and full vegetation cover, respectively.

[0045] The smaller the RPDI value calculated by this formula, the higher the soil moisture, i.e., the wetter the soil; the larger the RPDI value, the drier the soil. During irrigation, soil moisture increases sharply, and the RPDI will show a significant downward trend.

[0046] Specifically, the multi-temporal correction method for soil line slope involves selecting, for example, 500-1000 bare soil sample pixels with NDVI less than 0.15 from the study area at each time phase t, and linearly fitting their Red and NIR reflectance to obtain the soil line slope M for that time phase. t ;m t That is, M calculated on that day. t The value is adjusted to account for the seasonal variations in soil moisture, solar altitude angle, and other factors on the soil line. In phases where bare soil samples are insufficient, the M value from adjacent time phases can be used. t Perform linear interpolation.

[0047] In some optional implementations, to capture irrigation signals more comprehensively, the constructed time-series feature set may also include the Enhanced Vegetation Index (EVI) and the Water Index (WI2015) as auxiliary features. The EVI is calculated using the following formula:

[0048] EVI=2.5*(NIR-Red) / (NIR+6*Red-7.5*Blue+1);

[0049] Where Blue represents the reflectance in the blue light band. EVI is used to help determine the growth stage of crops, for example, determining the heading stage of rice by the EVI peak value. The formula for calculating WI2015 is:

[0050] WI2015=1.7204+171*Green+3*Red-70*NIR-45*SWIR1-71*SWIR2;

[0051] Wherein, Green represents the reflectance in the green light band, NIR represents the reflectance in the near-infrared band, SWIR1 represents the reflectance in the short-wave infrared band 1, and SWIR2 represents the reflectance in the short-wave infrared band 2. WI2015 is highly sensitive to surface water bodies and can be used to assist in identifying visible water signals at the beginning of irrigation.

[0052] Furthermore, considering the unavoidable high-frequency random noise in remote sensing observations, this embodiment also includes a step of smoothing the time series features for the generated daily RPDI time series data. The Savitzky-Golay (SG) filter, i.e., the least squares polynomial smoothing algorithm, is preferably used. The filter window size is set to 7 days, i.e., taking 3 days before and after the target date, and the polynomial fitting order is 2. This processing can effectively remove spike noise from the time series curve while better preserving the V-shaped abrupt change characteristics caused by irrigation, providing high-quality input data for subsequent irrigation identification models.

[0053] Specifically, the multi-temporal correction method for soil line slope is as follows: at each target time phase t, bare soil pixels with a Normalized Difference Vegetation Index (NDVI) less than 0.15 are selected from the study area image as samples, with the sample size typically ranging from 500 to 1000. Linear regression fitting is then performed on the red band reflectance (Red) and near-infrared band reflectance (NIR) of the bare soil pixels, with the fitting equation being NIR = M. t ×Red+b b b b Let M be the intercept term, where the fitted slope is M. t This is the slope of the soil line at that time phase, in m. t The direct value is M t .

[0054] The physical basis of this correction method is that the slope of the soil line is affected by factors such as soil type, soil moisture, and solar altitude angle. These factors change with the seasons and need to be updated sequentially. When the number of bare soil sample pixels in a certain time phase is insufficient, for example, less than 100, linear interpolation can be performed using the M value of the nearest time phase.

[0055] m t =M t1 ×(t2-t) / (t2-t1)+M t2 ×(t-t1) / (t2-t1);

[0056] Where t1 and t2 are the most recent effective soil line calculation dates, M t1 and M t2 These are the effective soil line slopes obtained from regression of bare soil samples at dates t1 and t2, respectively.

[0057] Example 3 describes the basic spatiotemporal fusion and plot constraint search strategy, and elaborates on the spatiotemporal data fusion process based on the weighted mechanism. In particular, it explains how to achieve high-precision image reconstruction in the highly heterogeneous farmland fragmentation area through optimized similar pixel search rules and multidimensional distance metrics, and solve the problem of spatial detail loss caused by the mixed pixel effect of low spatial resolution images.

[0058] Step 301: Construct a search window centered on the target pixel to be predicted, and filter similar pixels within the search window whose spectral features differ from the target pixel by less than a preset threshold.

[0059] In this embodiment, the target pixel refers to the pixel at prediction time t. p The pixel location (x, y) for reconstructing high-resolution reflectance values ​​needs to be determined. The search window is a rectangular area geometrically centered on the target pixel. Its size is typically set to an integer multiple of the low spatial resolution pixel size, such as 50 x 50 x 10-meter spatial resolution pixels, or a 500m x 500m range, corresponding to one MODIS pixel size. Alternatively, a larger window, such as 1500m x 1500m, can be used in more heterogeneous regions. The search window is constructed to determine a local spatial neighborhood, assuming that features within this neighborhood have similar spectral variation patterns to the target pixel.

[0060] Specifically, traditional screening methods rely solely on spectral thresholds, i.e., calculating the difference between candidate pixels and target pixels within a window at a reference time t. k Spectral differences (at times when high-resolution images are available) are considered similar pixels if the difference is less than a preset threshold. However, in the hilly irrigated areas of southern China, the land is fragmented and the planting structure is complex. Simple spectral screening can easily introduce pixels from different plots and different crops, i.e., heterogeneous pixels, leading to crosstalk or blurring in the fusion results.

[0061] This embodiment provides a preferred filtering strategy based on land parcel mask constraints. This strategy introduces land parcel mask data as prior knowledge. The land parcel mask data is a raster image where the value of each cell represents the farmland parcel ID it belongs to. During filtering, the same-parcel constraint is prioritized: only cells within the search window that have the same parcel ID as the target cell are selected as the candidate set. For example, if the target cell is located within a paddy field parcel with ID 101, the system only searches for similar cells within the region of ID 101.

[0062] When the number of similar pixels within the same plot that meet the spectral threshold is less than the preset minimum sample number N min When there are 10 (e.g.), a hierarchical expansion strategy is adopted:

[0063] The first level of expansion maintains the plot constraint and relaxes the spectral threshold, for example, by increasing the spectral difference threshold from 0.02 to 0.05, and re-screens within the same plot.

[0064] The second level of expansion, if the first level of expansion is still insufficient, removes the plot constraint and expands the search range to spatially adjacent plots, prioritizing adjacent plots with the same crop type, and supplementing the search with pixels whose spectral characteristics are closest from the adjacent plots.

[0065] In the third stage of processing, if there are still fewer than 5 similar pixels after two stages of expansion, the target pixel is marked as low quality and spatial fusion is not performed. Instead, linear interpolation in the time dimension is used to fill the gap, that is, linear interpolation is performed using the reflectance of the pixel at the time of the nearest effective high spatial resolution image before and after it.

[0066] It should be noted that the time interpolation strategy without extension described in Example 8 is suitable for application scenarios that require high purity of plot boundaries but relatively relaxed requirements for time accuracy, such as annual irrigation volume statistics; while the hierarchical extension strategy described in this example is suitable for dynamic monitoring scenarios that require precise capture of irrigation start dates.

[0067] The following three methods are available for obtaining land parcel mask data:

[0068] The vector data is rasterized to obtain vector cadastral maps or land use status maps provided by relevant departments. These maps contain the vector boundaries and unique numbers of each cultivated land parcel. Geographic Information System (GIS) tools are then used to rasterize the vector boundaries, generating a raster mask map with the same spatial resolution as the high spatial resolution imagery (e.g., 10 meters), where the value of each cell corresponds to the parcel's number.

[0069] Image segmentation was employed, using high spatial resolution remote sensing imagery such as Sentinel-2 or GF-2, and an object-oriented multi-scale segmentation algorithm was used to automatically extract parcel boundaries. Segmentation parameters were adjusted according to parcel fragmentation, for example, setting the segmentation scale to 50-100, the shape factor to 0.3, and the compactness to 0.5. The segmentation results were manually verified and then converted into a raster mask.

[0070] Deep learning semantic segmentation is performed by training semantic segmentation models such as U-Net (a U-shaped convolutional neural network) or DeepLab (a deep learning-based semantic segmentation model) using labeled plot boundary samples to extract plot boundaries from the images of the study area end-to-end.

[0071] This embodiment preferably uses vector data rasterization to ensure that the plot boundaries are consistent with the actual agricultural management units, which facilitates the subsequent integration of irrigation monitoring results with the farmland management system.

[0072] Step 302: For each similar pixel, calculate the spectral similarity distance, spatial distance, and temporal distance between it and the target pixel.

[0073] In this embodiment, to quantify the contribution of each similar pixel to the predicted value of the target pixel, the distance needs to be calculated from three dimensions: spectral, spatial, and temporal. The smaller the distance, the more similar or closer the similar pixel is to the target pixel, and the higher its reference value.

[0074] Specifically, spectral similarity distance reflects the spectral differences between pixels. For the i-th similar pixel, its spectral similarity distance S i Based on reference time t k It is calculated from high-resolution images. The calculation formula is:

[0075] S i =∑(∣F i,tk,b -F target,tk,b ∣) / B;

[0076] Among them, S i Let F be the spectral similarity distance of the i-th similar pixel, where ∑ represents the summation operation. i,tk,b For the i-th similar pixel at time t k The reflectivity of band b, F target,tk,b For the target pixel at time t k The reflectivity of band b, where B is the total number of bands involved in the calculation.

[0077] Spatial distance reflects the geometrical differences between pixels. Pixels closer to the target pixel generally have more similar geographical environments. The calculation formula is:

[0078] D i =sqrt((x i -x target ) 2 +(y i -y target ) 2 );

[0079] Among them, D i Let be the spatial distance between the i-th similar pixels, and sqrt denote the square root operation. i ,y i Let (x) be the coordinates of the i-th similar pixel. target ,y target The coordinates of the target pixel are shown below. For ease of subsequent calculations, the spatial distance is usually normalized to between 0 and 1, for example, by dividing it by half the diagonal length of the search window.

[0080] The temporal distance reflects the time span between the baseline image and the prediction image. When making predictions using two consecutive baseline images, t1 and t2, the temporal distance is used to measure which baseline image is closer to the prediction time t. p The calculation formula is as follows:

[0081] T k =∣t p -t k |;

[0082] Among them, T k Based on reference time tk The calculated time distance, t p t is the number of days to predict the date. k This refers to the number of days from the base image date. For example, if the forecast date is day 200, and the base dates are day 190 and day 210 respectively, then the reference value of the two is theoretically balanced.

[0083] Step 303: Based on spectral similarity distance, spatial distance and temporal distance, determine the weight contribution capability of each similar pixel. The weight contribution capability is used to calculate the normalized weight to generate the fused image.

[0084] In this embodiment, this step integrates the distances of the three dimensions mentioned above to form a comprehensive index that can characterize the credibility of similar pixels. Typically, the three distances are independent of each other, and therefore are integrated through linear combination or product. For example, a comprehensive distance index C can be constructed. i This makes C i =S i +a*D i +b'*T k Here, a and b' are coefficients that adjust the weights for different distances. The smaller the overall distance, the stronger the weight contribution of the similar pixels. This step establishes the basic physical principle that the more similar, closer, and shorter the time, the greater the weight.

[0085] Step 304: Utilizing the relationship between the changes in the high temporal-low spatial resolution image sequence and the low temporal-high spatial resolution image sequence at the reference time and the prediction time, and based on the correspondence between the reflectance of the high spatial resolution image and the reflectance of the low spatial resolution image at the reference time for similar pixels, calculate the conversion coefficient for each similar pixel.

[0086] In this embodiment, considering that land cover is not static, the reflectance relationship between low-spatial-resolution and high-spatial-resolution images may change linearly or non-linearly over time. Conversion coefficient V i This is used to capture this relationship, describing how changes in low-spatial-resolution pixels are mapped to changes in high-spatial-resolution pixels.

[0087] In practice, for each similar pixel i, its position at reference time t is extracted. k and prediction time t p Low spatial resolution image reflectance C i,tk and C i,tp and its time at reference time t k High spatial resolution image reflectance F i,tkBy assuming a linear relationship between changes in high-spatial-resolution pixels and changes in low-spatial-resolution pixels within a local window, the least squares method can be used to perform linear regression on the set of similar pixels, and the regression coefficients can be obtained as the transformation coefficients V. i This coefficient corrects for systematic biases caused by phenological changes or differences in lighting conditions, enabling the prediction model to incorporate not only spatial interpolation but also temporal evolution information.

[0088] Conversion factor V i The calculation formula is:

[0089] V i =(F i,tk -F i,tk_mean ) / (C i,tk -C i,tk_mean +ε);

[0090] Among them, F i,tk_mean C is the mean reflectance of a set of similar pixels in a high spatial resolution image at a reference time. i,tk_mean To represent the mean reflectance of the low spatial resolution image, ε is a small constant to prevent division by zero, preferably 0.0001. This coefficient describes the change in the high spatial resolution image corresponding to a unit change in the low spatial resolution image, and is used to correct the scale transformation error caused by the heterogeneity of ground features, thus achieving scale transformation.

[0091] Step 305: Calculate the comprehensive distance using the reciprocal or negative exponential relationship between spectral similarity distance, spatial distance, and temporal distance, and normalize the comprehensive distance to obtain the normalized weight.

[0092] In this embodiment, inverse distance weighting is typically used to transform the distance metric into normalized weights for weighted summation. The unnormalized weights R... i The calculation formula can be expressed as:

[0093] R i =(1 / (S i *D i *T k ));

[0094] Alternatively, the exponential form can be used when considering the logarithmic scale:

[0095] R i =exp(-S i / h s )*exp(-D i / h d )*exp(-T k / h t );

[0096] Among them, h sh d h t These are the attenuation parameters for the spectrum, space, and time, respectively.

[0097] The weights of all N similar pixels within the search window are normalized to ensure that the sum of all weights is 1. The normalized weights W... i The calculation formula is:

[0098] W i =R i / ∑(R j );

[0099] Among them, W i For the normalized weight of the i-th similar pixel, ∑(R j ) indicates that the unnormalized weights of all similar pixels j within the window are summed.

[0100] This step ensures that the prediction results are unbiased in terms of radiation magnitude.

[0101] Step 306: Based on the normalized weights and conversion coefficients, the reflectance of all similar pixels is weighted and summed to obtain the fused reflectance value of the target pixel at the monitoring time.

[0102] In this embodiment, this is the final output step of the spatiotemporal fusion algorithm. Based on the aforementioned calculated conversion coefficient V... i and normalized weight W i Combined with the high spatial resolution image reflectance F at the reference time i,tk And the change C of the low spatial resolution image at two time points i,tp -C i,tk The predicted value of the target pixel is calculated.

[0103] The specific calculation formula is as follows:

[0104] F target,tp =F target,tk +∑(W i *V i *(C i,tp -C i,tk ));

[0105] Among them, F target,tp For the target pixel at prediction time t p The fused reflectivity value, F target,tk For the target pixel at reference time t k The known high-resolution reflectance value, C i,tp Let C be the low spatial resolution image reflectance of the location corresponding to the i-th similar pixel at the prediction time. i,tk The reflectance of its low spatial resolution image at the reference time.

[0106] When there are two reference times t1 and t2 (where t1 < t p < t2), the predicted values F target,t1 and F target,t2 based on t1 and t2 are calculated respectively. The final fusion result adopts the time distance T k for weighted average:

[0107] F target,tp =(T2 × F target,t1 + T1 × F target,t2 ) / (T1 + T2);

[0108] where, T1 = |t p - t1| is the time interval from the prediction time to the first reference time (in days), and T2 = |t2 - t p | is the time interval from the second reference time to the prediction time.

[0109] The physical meaning of this formula is that the reference image closer to the prediction time has a greater contribution weight to the prediction result, ensuring a smooth transition of the fusion result in the time dimension.

[0110] For example, if t1 is the 190th day, t p is the 195th day, and t2 is the 200th day, then T1 = 5, T2 = 5, and the weights of the two reference images are each 50%; if t p is the 192nd day, then T1 = 2, T2 = 8, the weight of the predicted value based on t1 is 80%, and the weight of the predicted value based on t2 is 20%.

[0111] This formula makes the prediction result closer to the reference image with a shorter time distance, conforming to the assumption of temporal continuity of surface changes. The dual-reference fusion strategy can better handle the situation where the prediction time is between two observation days and smooth the surface change trend.

[0112] Through the above steps, a high spatio-temporal resolution fusion image of the entire map can be generated pixel by pixel and band by band.

[0113] Example 4 describes an improved IU-ESTARFM spatio-temporal fusion algorithm, and elaborates on an improved spatio-temporal fusion algorithm IU-ESTARFM that introduces a mutation protection mechanism and confidence output. It solves the technical problem that the basic ESTARFM algorithm often smooths the short-term surface wetting mutation caused by irrigation as noise or error, resulting in the elongation or weakening of the V-shaped feature. At the same time, it fills the gap in the lack of reliability measurement of the traditional fusion algorithm and provides a key quality gating signal for the subsequent recognition model.

[0114] Generating daily high spatiotemporal resolution fused image sequences covering the monitoring period, specifically including calculating fusion confidence, includes the following steps:

[0115] Step 401: Obtain daily rainfall data for the monitoring period.

[0116] In this embodiment, the rainfall data P t This is a key auxiliary variable for determining the causes of surface wetting. Its source can be measured data from meteorological stations or gridded precipitation products at global or regional scales. If station data is used, a raster precipitation map covering the study area needs to be generated using kriging interpolation or inverse distance-weighted interpolation. The spatial resolution of the precipitation data does not need to be very high, but the temporal resolution must be daily to match the monitoring period.

[0117] Step 402: Calculate the daily moisture index based on the shortwave infrared and near-infrared bands of the high temporal-low spatial resolution image sequence, and calculate the intensity of the sudden change in moisture on the target date based on the difference in moisture index between adjacent dates.

[0118] In this embodiment, the high temporal-low spatial resolution MCD43A4 V6.1 is used to capture changes in surface moisture. Specifically, the target date t is calculated. p The intensity of the wet change abruptly A tp The formula is:

[0119] A tp =max(0,J tp -J tp_prev );

[0120] Among them, A tp J represents the intensity of the wetted change on the target date. tp For target date t p Low spatial resolution moisture index, J tp_prev The day before the target date (t) p -1) is the low spatial resolution moisture index, and max indicates the operation of taking the maximum value.

[0121] The physical meaning of this formula is that it only considers the case where the moisture index increases significantly (i.e., becomes wetter). If the index decreases or remains unchanged (i.e., becomes drier or stable), then A... tp The value is 0. For example, if the value of J was 0.2 the previous day and the value of J was 0.5 on the current day, then A... tp =0.3, indicating a strong wetting mutation.

[0122] Step 403: Calculate the mutation protection factor based on rainfall data and the intensity of wet change mutation. The mutation protection factor is used to indicate the degree of wet change caused by non-rainfall.

[0123] In this embodiment, this is the step to distinguish between natural rainfall-induced wetting and artificial irrigation-induced wetting. If surface wetting is detected, i.e., A... tp >0, and the daily rainfall P tp If the rainfall is significant, the increased wetness is likely caused by widespread rainfall, in which case no specific protection is needed, and the original spatiotemporal integration logic should be maintained. Conversely, if there is little or no rain, but the surface is noticeably wet, it is likely due to irrigation, in which case the mutation protection factor E needs to be introduced. tp This is to amplify the mutation signal and prevent it from being smoothed out by surrounding unirrigated pixels or historical imagery.

[0124] The mutation protection factor is calculated based on rainfall data and the intensity of wet change abruptly, according to the following rules:

[0125] If the rainfall data for the target date is greater than the preset rainfall threshold, or the intensity of the wet change mutation is zero, then the mutation protection factor will be set to the baseline value of 1.

[0126] If the rainfall data for the target date is less than or equal to the rainfall threshold, and the wet change mutation intensity is greater than zero, then the baseline value is amplified using the wet change mutation intensity to obtain the mutation protection factor. The greater the wet change mutation intensity, the larger the value of the mutation protection factor.

[0127] This embodiment provides E tp The specific calculation formula is as follows:

[0128] E tp =1+ρ*(A tp / (mean Atp +θ))*I(P tp <=P th );

[0129] Among them, E tp ρ is the mutation protection factor; ρ is the mutation gain coefficient, used to control the amplification factor, preferably ranging from 1.0 to 3.0, for example, set to 2.0; mean Atp θ is the average of the humidity intensity of all pixels within the region, used for normalization; θ is the abrupt saturation constant, used to prevent the denominator from being too small, preferably set to 0.05; P tp The target daily rainfall; P th The preset rainfall threshold is set to 5mm; I(...) is an indicator function that takes the value 1 when the condition in parentheses is met (i.e., non-rainfall-dominated), and 0 otherwise.

[0130] For example, suppose the rainfall threshold P th =5mm, no rain that day, i.e. P tp =0, and a pixel was detected to have become intensely wet. tp =0.4, regional average humidity meanAtp =0.1. Setting ρ=2.0 and θ=0.05, the mutation protection factor for this pixel is calculated as follows:

[0131] E tp =1+2.0*(0.4 / (0.1+0.05))*1=1+2.0*2.67=6.34.

[0132] It can be seen that the value of this factor is much greater than 1, which will significantly amplify the competitiveness of the pixel weight in subsequent steps.

[0133] The selection criteria for the above parameters are explained below:

[0134] Rainfall threshold P th The threshold of 5 mm is based on agricultural meteorological research, where daily rainfall less than 5 mm is generally defined as ineffective or trace rainfall. This type of rainfall has little impact on soil surface moisture content and is unlikely to produce significant changes in reflectance in remote sensing images. When daily rainfall exceeds 5 mm, surface wetness can be reasonably attributed to rainfall; conversely, the possibility of artificial irrigation needs to be considered. This threshold can be adjusted according to the climatic characteristics and soil type of the study area; for example, it can be appropriately lowered to 3 mm in sandy soil areas and appropriately increased to 8 mm in clay soil areas.

[0135] The mutation gain coefficient ρ controls the mutation protection factor E tp The amplification factor of the original weights. If ρ is too small, such as ρ < 1.0, the amplification effect is insufficient, and abrupt signals may still be smoothed; if ρ is too large, such as ρ > 3.0, it may cause a few abnormal pixels to excessively dominate the fusion result. Experiments show that when ρ is 2.0, a good balance is achieved between preserving abrupt signals and maintaining fusion stability.

[0136] The abrupt saturation constant θ is to prevent the average regional wetting intensity from changing when the mean saturation intensity is reached. Atp When the value is small, division by zero errors or numerical instability may occur. The value of θ=0.05 corresponds to 5% of the change in reflectivity (normalized to 0-1). This order of magnitude ensures numerical stability without excessively suppressing small but real humidity changes.

[0137] The sliding window length of L=20 days is set based on the following considerations: the single irrigation cycle for major crops such as rice is typically 7 to 14 days, and a 20-day window can completely cover the V-shaped cycle of irrigation-water retention-consumption. If the window is too short, such as L<15 days, the V-shaped features may be incomplete, increasing the false negative rate; if the window is too long, such as L>30 days, multiple irrigations may occur within a single window, increasing the difficulty of model learning. Experimental verification shows that L=20 days achieves an optimal balance between detection accuracy and computational efficiency.

[0138] Step 404: Use the mutation protection factor to correct the initial weights determined based on spectral similarity distance, spatial distance and temporal distance to obtain the corrected normalized weights.

[0139] In this embodiment, the calculated mutation protection factor E tp This is injected into the weight generation process. In Example 3, the initial integrated distance or weight is mainly determined by spectral, spatial, and temporal distances. In this step, E is used... tp The initial weights are multiplicatively corrected.

[0140] Specifically, the corrected unnormalized weight R v_new The calculation formula is:

[0141] R v_new =R v_old *E tp ;

[0142] Among them, R v_old The weights are the original weights calculated solely based on 3D distance; that is, the smaller the distance, the greater the weight. R v_new This is the corrected weight. Continuing with the above example, since E... tp =6.34, the weight of this pixel has been amplified by more than 6 times. This indicates that during weighted fusion, the system highly trusts the information of low spatial resolution image pixels that show non-rainfall abrupt changes and forcibly injects their changing trends into the fusion result, preserving the V-shaped left-wing descent feature caused by irrigation and avoiding averaging by surrounding unwetted pixels.

[0143] The corrected weights are then normalized.

[0144] W v =R v_new / ∑(R k_new );

[0145] Among them, W v For the final normalized weights used for fusion, ∑(R k_new ) is the sum of the corrected weights of all similar pixels within the search window.

[0146] Step 405: The reflectance of similar pixels is weighted and summed using the corrected normalized weights to obtain the fused image sequence. The fusion confidence is calculated based on the distribution characteristics of the corrected normalized weights. The fusion confidence is used as the gating input for the subsequent irrigation identification model.

[0147] In this embodiment, in addition to generating the fused image, the fusion confidence score q is also output. tThis metric is used to quantify the reliability of the current pixel fusion result. If there are very similar pixels with high weights within the search window, i.e., the weight distribution is concentrated, the fusion result is reliable; if all similar pixels have roughly the same weight, i.e., the weight distribution is uniform and chaotic, it indicates that no good reference pixels can be found, and the fusion result has greater uncertainty.

[0148] The fusion confidence score is calculated based on the distribution characteristics of the corrected normalized weights, using the following formula:

[0149] H tp =-∑ v∈Ω (W v *ln(W v ));

[0150] q tp =1-(H tp / ln(N)), 0≤q tp ≤1;

[0151] Among them, H tp W represents the information entropy of the weight distribution. v Ω is the corrected normalized weight of similar pixels v; ln is the natural logarithm function; N is the total number of pixels in the similar pixel set; q tp The fusion confidence level for the target date, ranging from 0 to 1.

[0152] For example, if there are N=50 similar pixels in the window.

[0153] Scenario A, high confidence: For one pixel, W1 = 0.9, while the remaining 49 pixels have lower values. In this case, the entropy H... tp Approaching 0, q tp Approaching 1. This usually occurs when an ideal plot of land with the same quality of pixels is found.

[0154] Scenario B, low confidence: all 50 pixels have equal weights, each 0.02. In this case, the entropy H... tp The maximum value ln(50) is reached, causing q to... tp =1-1=0. This usually occurs in areas with cloud cover, highly fragmented terrain features, or poor data quality.

[0155] The final output of daily q tp The raster image will serve as the gating input for the subsequent confidence-gated and event-aligned bidirectional long short-term memory network model CGEA-BiLSTM, informing the model which days' data are reliable and which days' data are unreliable.

[0156] Example 5 details a method for constructing and applying an irrigation identification model based on a Bidirectional Long Short-Term Memory (Bi-LSTM) network. Leveraging the advantages of deep learning technology in handling time-series dependencies, this method replaces traditional identification methods based on simple thresholds or static classifiers. It can effectively capture the dynamic evolution characteristics of irrigation events on the RPDI time-series curve. The irrigation identification process is as follows: Figure 3 As shown.

[0157] Step 501: The time-series feature sequences within the sliding window are normalized and then mapped into high-dimensional feature vectors.

[0158] In this embodiment, a sliding window is a time-series data preprocessing technique used to divide a continuous long sequence into fixed-length local subsequences, facilitating the model's capture of local event features. Specifically, the sliding window length L is set to 20 days, with a step size of 1 day. For each pixel within the study area, subsequences of length 20 are sequentially extracted from its complete RPDI time series from May 25th to November 5th. For example, the first window covers days 1 to 20, the second window covers days 2 to 21, and so on.

[0159] Specifically, in this embodiment, the time-series feature set only includes the Recorrected Vertical Drought Index (RPDI) sequence. The pre-built irrigation identification model is to directly use the smoothed RPDI sequence as input and use a bidirectional long short-term memory network to capture the V-shaped change features of the RPDI sequence before and after irrigation.

[0160] RPDI directly reflects the relative change of soil moisture with respect to vegetation background, making it the strongest carrier of irrigation signals. Before inputting it into the model, RPDI data typically needs to be normalized, for example, through Z-score normalization or Min-Max normalization, mapping it to a distribution with a mean of 0 and a variance of 1, or an interval between 0 and 1, to accelerate the convergence of the neural network. After mapping, each sliding window is converted into a tensor of dimension (L,1), i.e., (20,1), where L represents the time step and 1 represents the feature dimension, specifically RPDI. This tensor is the initial high-dimensional feature vector fed into the neural network.

[0161] In some alternative implementations, the length L of the sliding window can be adjusted according to local irrigation practices, provided computational resources allow. For example, in areas with longer irrigation durations, L can be set to 30 days; while in areas with higher irrigation frequency, L can be set to 15 days. Furthermore, in addition to Z-score normalization, decimal scaling normalization can be used, which involves dividing the data by 10 to the power of k to make its absolute value less than 1.

[0162] Step 502: Use a bidirectional long short-term memory network to extract the dependencies of the temporal feature sequences in the forward and backward time dimensions, respectively, and generate the hidden state sequence.

[0163] In this embodiment, the Bi-LSTM bidirectional long short-term memory network consists of forward LSTM layers and backward LSTM layers. Standard recurrent neural networks (RNNs) are prone to gradient vanishing or gradient exploding problems during long-sequence training, while LSTMs effectively solve this problem by introducing gating mechanisms (forget gate, input gate, and output gate), enabling them to remember key information for a long time. The bidirectional Bi-RNN structure further allows the model to see not only past information but also future information at the current moment, which is particularly important for judging irrigation events, as the determination of irrigation often depends on the upward trend after irrigation, i.e., the right wing of a V-shaped recovery.

[0164] Specifically, the computation process inside an LSTM cell can be described by the following linear formula:

[0165] f t =σ(W f *x t +U f *h t-1 +b f );

[0166] i t =σ(W i *x t +U i *h t-1 +b i );

[0167] o t =σ(W o *x t +U o *h t-1 +b o );

[0168] C t =f t *C t-1 +i t *tanh(W c *x t +U c *h t-1 +b c );

[0169] h t =o t *tanh(C t );

[0170] Where, x tThe input at time t is the RPDI value for that day, h t-1 h represents the hidden state of the previous time step. t f is the hidden state at the current moment. t i t o t C represents the activation values ​​of the forget gate, input gate, and output gate, respectively. t In cellular state, C t-1 Let σ represent the cell state at the previous time step, σ be the sigmoid activation function, tanh be the hyperbolic tangent activation function, and W be the cell state at the previous time step. f W i W o W c U is the weight matrix input to each gate / candidate memory. f U i U o U c Let b be the weight matrix from the hidden state to each gate / candidate memory. f b i b o b c This is the bias vector for each gate / candidate memory, where * indicates matrix multiplication or element-wise multiplication.

[0171] In this embodiment, the preferred model configuration is a two-layer stacked Bi-LSTM structure. Each Bi-LSTM layer has 64 hidden units. For each time step t, the forward LSTM outputs the hidden state h. f,t The hidden state h is then output to the LSTM. b,t The system concatenates the two vectors to obtain the combined hidden state h at the current time step. t =[h f,t ,h b,t For a window of length L=20, the model ultimately generates a hidden state sequence of length 20:

[0172] H={h1,h2,...,h 20}, where each h t The dimension is 128, or 64+64. To prevent overfitting, a Dropout layer is introduced between the two Bi-LSTM layers, with a dropout rate set to 0.5.

[0173] Step 503: Perform feature aggregation on the hidden state sequence and output the irrigation probability corresponding to the sliding window through a fully connected layer.

[0174] In this embodiment, feature convergence is used to compress variable-length or serialized hidden state information into a fixed-length window-level feature vector for classification. The simplest convergence method is to directly take the hidden state h_L of the last time step, or to take the average of the hidden states of all time steps. In the basic scheme of this embodiment, average pooling or max pooling strategies are preferred, that is, the sequence H is averaged or maximized in the time dimension to obtain the window feature vector Z.

[0175] Specifically, the fully connected layer acts as a classifier, mapping the window feature vector Z to the class space of the samples. The calculation formula is:

[0176] y=σ(W fc *Z+b fc );

[0177] Where Z is the converged feature vector, and W fc Let b be the weight matrix of the fully connected layer. fc σ is the bias, σ is the Sigmoid function, and y is the output predicted probability value, which ranges from 0 to 1.

[0178] In this embodiment, the probability value y represents the likelihood of irrigation occurring on the date corresponding to the center or end of the sliding window. Through training, this model can accurately capture the V-shaped feature of the RPDI sequence. Specifically, when the input RPDI sequence exhibits a pattern of first a sharp decrease (irrigation), maintaining a low value (water retention), and then a slow increase (depletion), the temporal feature Z extracted by Bi-LSTM will be closer to the irrigation class center in the feature space, causing the output probability y to approach 1. In practical applications, a discrimination threshold is typically set, for example, 0.5. When y > 0.5, it is determined that an irrigation event has occurred in the window.

[0179] In some optional implementations, if the specific start date of irrigation needs to be output, not only can the end of the window be classified, but a many-to-many output mode can also be adopted, that is, outputting a probability value at each time step to achieve daily marking of the irrigation process. Furthermore, the Adam optimizer is preferably used during model training, with an initial learning rate of 0.001, a batch size of 32, and cross-entropy as the loss function.

[0180] Example 6 details a bidirectional long short-term memory network (CGEA-BiLSTM) model with confidence gating and event alignment, addressing the technical problem of misjudgment in traditional Bi-LSTM when facing local noise from spatiotemporal fusion, cloud shadow interpolation errors, and rainfall interference. By deeply coupling the fusion confidence and rainfall data into the recognition network, adaptive suppression of external interference factors such as rainfall and accurate alignment of the spatiotemporal location of irrigation events are achieved.

[0181] Step 601: Obtain the daily RPDI sequence, fusion confidence sequence, and rainfall data within the sliding window, and construct the daily composite input vector.

[0182] In this embodiment, the input layer is no longer a single-dimensional RPDI, but is expanded into a multi-dimensional composite vector containing feature terms, quality terms, and environmental terms. Specifically, for each time step t, t=1,...,20, the system obtains three components: the Re-corrected Vertical Aridity Index (RPDI). t ; Fusion confidence q t ; and rainfall P t .

[0183] Specifically, the three components are concatenated to construct a composite input vector X. t For a sliding window of length L, the dimension of the input tensor changes from (L,1) to (L,3). Here, RPDI... t It carries physical signals of changes in surface moisture; q t It carries prior knowledge (values ​​from 0 to 1) regarding whether the RPDI data is reliable at that moment; P t It carries meteorological background information. This design allows the neural network not only to see what changes have occurred, but also to determine whether these changes are real and whether they were caused by rainfall.

[0184] Step 602: Before inputting the composite input vector into the bidirectional long short-term memory network, the input features are gated and weighted using the fusion confidence to suppress the feature responses at low confidence time steps.

[0185] In this embodiment, to address the common noise problem in remote sensing time-series data, this invention does not simply remove low-quality data, which would disrupt temporal continuity. Instead, a soft gating mechanism is designed. This mechanism is located before the Bi-LSTM layer, as part of a preprocessing layer or embedding layer. The gating coefficient g is used... t The input features are weighted, and when the data quality is poor or there is strong interference from rainfall, the amplitude of the input at that moment is forcibly reduced to be close to 0, thereby blocking the propagation of noise into the hidden state of the LSTM.

[0186] First, calculate the first difference of the RPDI sequence to capture the trend of change:

[0187] ΔR t =RPDI t -RPDI t-1 ;

[0188] Where, ΔR t This represents the change on day t relative to the previous day.

[0189] A rainfall suppression term is constructed using rainfall data, and the gating coefficient is calculated by combining the fusion confidence level. The gating coefficient g is then constructed. t This coefficient takes into account both fusion quality and rainfall interference. The calculation formula is:

[0190] g t =q t *(1-min(1,P t / P th ));

[0191] Among them, g t q is the gating coefficient, ranging from 0 to 1; t For fusion confidence; P t P represents the rainfall for that day. th t represents the rainfall threshold, for example, 5mm; min represents the minimum value operation; and t represents the time.

[0192] This formula means:

[0193] Case 1, if q t If the value is close to 0, such as in the cloud shadow region or the fusion failure region, then regardless of P... t How, g t All values ​​approach 0, and the model automatically ignores the input at that moment.

[0194] Case 2, if P t Greater than P th On days with heavy rainfall, (1-1) = 0, g t A value of 0 indicates that even if the RPDI decreases significantly, i.e., becomes wetter, the model will still be affected by g. t The inhibitory effect of rainfall is ignored, effectively eliminating false irrigation signals caused by rainfall.

[0195] Scenario 3: Only when q t High (data reliable) and no rain (P) t When =0), g t Only when it approaches 1 does the variation in RPDI fully enter the network. By X t With g t Perform element-wise multiplication to generate a weighted input vector X. t 'Input Bi-LSTM.'

[0196] Step 603: When performing feature aggregation on the hidden state sequence, an attention mechanism is introduced. The event alignment score calculated based on the RPDI decrease magnitude and fusion confidence is used to calculate the attention weight, so that the model focuses on the non-rainfall wet change abrupt moment with high confidence.

[0197] In this embodiment, to address the phase drift problem that may arise from a fixed sliding window—that is, irrigation may occur on the 3rd day or the 18th day of the window—this invention proposes an attention mechanism. Unlike traditional attention mechanisms that rely solely on the self-similarity of hidden states, this scheme utilizes the left-wing descent event score s. t This is to explicitly guide the network to focus on important time nodes.

[0198] The left-wing descent event score is obtained by weighting the descent magnitude of the first-order difference based on the gating coefficient. t The calculation formula is:

[0199] s t =g t *max(0,-1*ΔR t );

[0200] Where, -1*ΔR t This indicates the magnitude of the decrease in RPDI, i.e., the degree of humidification; max(0,...) ensures that only the downward trend (humidification) is considered, and if RPDI increases (drierness), the score is 0; g t The aforementioned gating coefficient. This formula defines a relatively ideal candidate irrigation time: it must be a time when the data is reliable (q). t High), non-rainfalling (P) t Low), and undergoes drastic wetting (ΔR) t (Large negative values). s t The larger the value, the greater the likelihood that irrigation will occur at that moment.

[0201] Further utilize s t Calculate attention weight a t The calculation formula is:

[0202] a t =exp(μ*s t ) / ∑(exp(μ*s k ));

[0203] Among them, a t is the normalized attention weight at time step t; μ is the attention temperature coefficient, used to control the degree of attention focus, preferably ranging from 1.0 to 5.0, for example, set to 2.0. The larger μ is, the more concentrated the weight is on the moment with the highest score; exp is the exponential function; ∑ represents the summation over all time steps k within the window.

[0204] Using a t The hidden state sequence H={h1,...,h} output by Bi-LSTM 20 We perform weighted aggregation to obtain the window-level feature vector z. u :

[0205] z u =∑(a t *h t );

[0206] Among them, h t Let be the hidden state vector at time t.

[0207] Through the above mechanism, regardless of where irrigation occurs within the window, the attention weight a t The system will automatically display the moment of mutation, resulting in the aggregated feature z. u It always includes the most obvious irrigation behavior features, achieving a soft alignment that makes the model more robust to the location of the irrigation start date, no longer relying on the assumption that irrigation must occur in the center of the window.

[0208] For example, suppose the parameter for day t=10 within a sliding window is: fusion confidence q 10 =0.85, rainfall P 10 =0mm, rainfall threshold P th =5mm, RPDI9=0.3, RPDI 10 =0.1. Therefore:

[0209] First-order difference: ΔR 10 =0.1-0.3=-0.2, indicating that it becomes wetter;

[0210] Gating coefficient: g 10 =0.85×(1-min(1,0 / 5))=0.85×1=0.85;

[0211] Left-wing descent event score: s 10 =0.85×max(0,-(-0.2))=0.85×0.2=0.17;

[0212] The higher event score will cause the attention mechanism to give greater weight to day 10 in subsequent calculations, allowing the model to focus on the moment when irrigation occurs.

[0213] In contrast, if day t=12 is a rainy day, P 12 =15mm>5mm, even if RPDI decreases, its gating factor g 12 =0.85×(1-1)=0, event score s 12 =0, this moment will be automatically ignored by the model, eliminating false irrigation signals caused by rainfall.

[0214] In some alternative implementations, to enhance the model's adaptability to different crop growth stages, the temperature coefficient μ can be set as a learnable parameter, allowing the network to automatically adjust its alignment sensitivity during training. Furthermore, in addition to using Softmax (i.e., exp / ∑) for normalization, a sigmoid function can be used to score each time step independently to adapt to complex scenarios where multiple irrigations (such as intermittent irrigation) may occur within a window.

[0215] Example 7 details the offline construction process of the irrigation identification model, addressing how to prevent overfitting of the model to low-quality data when cloud shadow noise, interpolation errors, and rainfall interference samples are unavoidable in the training data, thus training a more robust CGEA-BiLSTM model. By introducing a loss weighting mechanism based on fusion confidence, reinforcement learning on high-quality samples and automatic suppression of low-quality samples are achieved.

[0216] Step 701: Construct a training sample set containing a set of historical time-series features and corresponding irrigation labels.

[0217] In this embodiment, the quality of the training sample set determines the model's recognition accuracy. The construction process consists of three stages: sample selection, feature extraction, and label creation. Based on field survey data within the study area, such as farmland irrigation records from June to November 2025, typical plots are selected as ground truth samples. Simultaneously, to expand the sample size, reservoir water release records and hydrological monitoring station data are combined, and typical V-shaped features from remote sensing images are used to assist in labeling some samples. Specifically, for each sample point, its complete time-series data within the monitoring period is extracted, including the RPDI sequence and the fused confidence sequence q. t And rainfall data P t .

[0218] Samples are generated using the same sliding window strategy as in the prediction phase. The window length L is set to 20 days, and the step size is 1 day. For each sliding window, the input feature X is a matrix of dimension (20,3). The corresponding label Y is a binary category, 0 or 1.

[0219] For each sliding window, the time correspondence of its label Y is defined as follows, with the label determination position set to the end date t of the window. end If the on-site records or auxiliary judgments are displayed in the time window [t] start ,t end [Inside, especially with t] end If at least one irrigation event occurs within a 3-day period before and after the center, the window is labeled Y=1; otherwise, Y=0. When outputting the dynamic distribution map of irrigated area, windows determined to be irrigated (i.e., with an output probability greater than the threshold) are mapped to the window's last date t. endThe spatial location of the pixel is used as the irrigation detection result for that pixel on that date.

[0220] The reason for using the end of the window as the label location is that irrigation decisions require observing a complete 'V'-shaped trend in RPDI from its decline to its rebound, and this trend only fully manifests several days after irrigation. Mapping the label to the end of the window ensures that the model has already identified sufficient subsequent information when making the decision.

[0221] Using the above method, a dataset of 400 sample pairs (X,Y) was constructed. To ensure the model's generalization ability, the dataset is typically randomly divided into a training set and a validation set in a 7:3 ratio.

[0222] Irrigation labels can be created from the following three sources.

[0223] Source 1: Field survey data. During the monitoring period, irrigation records were obtained by visiting irrigation district management departments and typical farmers, including irrigation dates, irrigated plot numbers, and irrigation volume. This type of data has the highest accuracy, but its coverage is limited.

[0224] The second source is water conservancy facility monitoring data. This involves using flow monitoring data from irrigation canal head gates, reservoir water release records, and pump station operation records to infer the time frame and beneficiary areas of irrigation. This type of data has a wide coverage, but there is some uncertainty in its spatiotemporal accuracy.

[0225] Source 3: Remote sensing-assisted annotation. Experienced remote sensing interpreters manually interpret and annotate typical spectral variation features on the original Sentinel-2 imagery, such as the V-shaped drop in RPDI and clear water signals in visible light images. To ensure annotation quality, each sample is independently annotated by at least two interpreters, and only samples with consistent annotations are retained.

[0226] When constructing the training sample set, data from Source 1 is used first as the positive sample (Y=1), followed by data from Sources 2 and 3 to supplement it. The selection principle for negative samples (Y=0) is to choose a spatiotemporal window with no rainfall or irrigation records during the monitoring period, and where the RPDI curve is stable without a significant decline. The ratio of positive to negative samples is controlled between 1:2 and 1:3 to alleviate the class imbalance problem.

[0227] Step 702: For each training sample corresponding to a sliding window, calculate the gating coefficient g at each time step within the window. t The statistical mean is used to obtain the loss weight ω for this training sample. k .

[0228] In this embodiment, not all training samples have equal importance or confidence. Traditional training methods typically assign equal weights (i.e., a weight of 1) to all samples, causing the model to forcibly fit dirty samples generated by cloud shadow interpolation errors or rainfall noise, thus undermining the model's learning of real irrigation features. This invention proposes a confidence-based sample weighting strategy. For the k-th training sample, i.e., the k-th sliding window, the gating coefficient g at each time step t within the window is first calculated. t .

[0229] Specifically, the loss weight ω k The calculation formula is:

[0230] ω k =(1 / L)*∑(g t );

[0231] Where, ω k Let L be the loss weight for the k-th sample, and L be the window length, for example, 20, ∑(g t The formula sums the gating coefficients over all time steps within the window. This formula calculates the average data quality level within the window.

[0232] For example, for a high-quality sample of irrigation during sunny days, the fusion confidence q of most time steps is... t High and no rain, g t Approaching 1, the final calculated ω k It is also close to 1. Conversely, for a sample that is severely affected by continuous cloud and rain coverage and is mainly generated by interpolation, its g t Generally low, the calculated ω k It may be only 0.2 or lower. ω k This quantifies the influence of the sample during the training process.

[0233] Step 703: The predicted loss of the training sample is weighted using the loss weights, and the model parameters are updated using the backpropagation algorithm based on the weighted total loss to reduce the impact of low-confidence samples on model training.

[0234] In this embodiment, model training aims to minimize the difference between the predicted result and the true label. Binary cross-entropy is used as the basic loss function. Loss weights ω are introduced. k Then, the weighted loss l of the k-th sample k The definition is as follows, and the calculation formula is:

[0235] l k =-1*ω k *(y k *ln(p k )+(1-y k )*ln(1-pk ));

[0236] Among them, l k For the weighted loss of the k-th sample, y k For real labels, 0 or 1, p k Let ln be the predicted probability output by the model, and ln be the natural logarithm.

[0237] Each training batch contains M samples, for example, M=32. The total loss of the system is L. total The average weighted loss of all samples within this batch:

[0238] L total =(1 / M)*∑(l k );

[0239] Where ∑ represents the summation over all samples k in the batch.

[0240] Calculate L using the backpropagation algorithm total The model parameters, such as the weight matrix and bias gradient of the LSTM, are updated using the optimizer.

[0241] Specifically, the Adam optimizer is preferred, with an initial learning rate set to 0.001. During training, the model will preferentially learn ω. k Larger, high-confidence samples generate larger gradients and contribute more significantly to parameter updates. However, for ω... k Even with smaller, lower-quality samples, the loss from model prediction errors will be reduced by ω. k This narrowing of the model prevents it from deviating from the correct decision boundary in order to accommodate noisy data. This mechanism improves the model's robustness under complex weather conditions.

[0242] In some alternative implementations, in addition to taking the mean, the loss weight ω k You can also retrieve g within the window t The minimum or median value can be used to implement stricter quality control. For example, if a minimum value strategy is used, the weight of the entire sample will be downgraded if the data quality is poor on any one day within the window, which is suitable for applications with high requirements for data purity. In addition, to prevent the weights from being too small and causing training instability, a basic weight threshold can be set, such as 0.1, to ensure that all samples contribute at least slightly.

[0243] Once the irrigation recognition model is trained, it can be deployed in the following environments for inference applications.

[0244] Batch processing is deployed in the cloud, and the trained model is deployed on the cloud computing platform to perform batch inference on the fused image sequence of the entire irrigation area to generate a daily dynamic distribution map of irrigation area.

[0245] Deployed on a local workstation, the model is exported to the open neural network exchange format ONNX or TensorFlowLite, and deployed on a local workstation equipped with a GPU. After the user uploads image data of a specified area, the system automatically completes the fusion and recognition, and promptly returns the irrigation detection results.

[0246] The inference process is as follows: load the pre-trained model weights; perform normalization preprocessing on the input RPDI sequence (and confidence sequence, rainfall data, if CGEA-BiLSTM is used); input the model by sliding window slices; the model outputs the irrigation probability of each window; mark windows with probabilities greater than a threshold (e.g., 0.5) as irrigation events; map the detection results back to geographic coordinates to generate an irrigation distribution map.

[0247] Example 8 provides an alternative implementation scheme for application scenarios lacking refined meteorological data, such as daily rainfall, or where computational resources are limited, making it difficult to run attention mechanisms. This scheme does not introduce the mutation protection factor of IU-ESTARFM or the gating alignment mechanism of CGEA-BiLSTM. It optimizes the basic fusion process through plot mask constraints and utilizes basic Bi-LSTM for irrigation identification. This effectively solves the problems of mixed pixel interference in low spatial resolution images and the low accuracy of single threshold methods, and is a basic implementation of this invention.

[0248] Step 801: Construct a search window centered on the target pixel to be predicted, and filter similar pixels within the search window whose spectral features differ from the target pixel by less than a preset threshold.

[0249] In this embodiment, as an optimization scheme for heterogeneous and fragmented farmland areas, this step focuses on implementing a plot masking constraint strategy. Unlike general search, this embodiment enforces strict same-plot search logic to minimize spectral confusion across plots.

[0250] Specifically, first, a vector plot map of the study area is loaded and rasterized to match the spatial resolution of the high spatial resolution image, for example, a 10-meter plot mask M. For a target pixel at location (x,y), its plot ID is identified. target After constructing a search window of size w*w, for example, 51*51, the system only retains the elements within the window that satisfy M(i,j)=ID. target The pixels are used as the candidate set.

[0251] The spectral Euclidean distance between candidate pixels and target pixels on the reference image is calculated, and the N pixels with the smallest distance are selected as the final similar pixels. If there are fewer than N candidate pixels in the same plot, the system will not expand outward, but will directly use the weighted average of all valid pixels in the same plot as the prediction result, or mark the pixel as low quality and fill it using linear interpolation in the time dimension.

[0252] While this strategy may make some boundary cells difficult to predict, it ensures the purity of the predicted cells within the field and avoids the erroneous fusion of spectral features from nearby roads or woodlands into the farmland.

[0253] Step 802: Calculate the conversion coefficient for each similar pixel by utilizing the relationship between the high temporal-low spatial resolution image sequence and the low temporal-high spatial resolution image sequence at different times.

[0254] In this embodiment, the conversion coefficient calculation logic of standard ESTARFM is adopted. It is assumed that the reflectance of ground objects changes linearly over a short period of time.

[0255] For each selected similar pixel k, its value at reference time t is used. k and prediction time t p Low spatial resolution image reflectance to (C k,tk C k,tp ) and the low spatial resolution image reflectance of the target pixel (C target,tk C target,tp Linear regression analysis was performed.

[0256] The calculation formula is: C k,tp =a k *C k,tk +b k ;

[0257] Among them, a k b is the gain coefficient. k This is the offset. Applying this linear relationship to high spatial resolution scales, we assume the transformation coefficient V... k Approximately equal to a k This step ensures that the fusion model can capture the overall brightness fluctuations caused by phenological changes.

[0258] In some embodiments, the process may further involve calculating global conversion coefficients using the low spatial resolution image reflectance variation relationship of all similar pixels within the search window at the reference time and the prediction time. Specifically, for all N similar pixels within the window, the low spatial resolution image reflectance pairs (C... k,tk C k,tp For k=1,...,N, a linear regression is performed, and the fitted model is:

[0259] C k,tp =a*C k,tk +b a ;

[0260] Where 'a' is the common slope fitted together from all similar pixels, i.e., the uniform gain coefficient of the entire search window, and 'b' is the common slope fitted together from all similar pixels. a It is the common intercept fitted together from all similar pixels, that is, the uniform offset of the entire search window. The slope 'a' obtained from the regression is used as the uniform transformation coefficient within this search window and applied to all similar pixels.

[0261] Step 803: Calculate the comprehensive distance using the reciprocal or negative exponential relationship between spectral similarity distance, spatial distance, and temporal distance, and normalize the comprehensive distance to obtain the normalized weight.

[0262] In this embodiment, the weights are calculated based on physical distance, without introducing a mutation protection factor. This is to adapt to scenarios where there is no rainfall data input. The comprehensive distance D... k The calculation formula is:

[0263] D k =1.0-(S k / S max )*(D space_k / D space_max )*(T k / T max );

[0264] Among them, S k For spectral differences, S max D represents the maximum spectral difference among all candidate pixels within the search window. space_k D is the spatial Euclidean distance (the pixel or geographic distance between the target pixel and the k-th similar pixel). space_max T is the maximum spatial radius of the search window, such as half the side length of the window. k T represents the number of days in the time interval. max This is the maximum effective time window set.

[0265] Alternatively, a more general inverse distance form can be used:

[0266] W k_raw =(1 / S k )*(1 / D space_k )*(1 / T k );

[0267] Among them, W k_raw The original weight (unnormalized) of the k-th similar pixel.

[0268] Then normalization was performed, W k =Wk_raw / ∑(W i_raw ); where W k The final weight is represented by ∑, which indicates the summation over all similar pixels within the window.

[0269] Step 804: Use a sliding window to slice the time-series feature set, input the pre-built irrigation identification model, analyze the time-series evolution features within the window, determine whether an irrigation event has occurred at the spatiotemporal location corresponding to the window, and obtain a dynamic distribution map of the irrigation area.

[0270] In this embodiment, as an alternative, the irrigation identification model adopts a standard two-layer Bi-LSTM network, namely a bidirectional long short-term memory network Bi-LSTM structure, and the input feature is only a univariate RPDI sequence, that is, the feature dimension is 1. Although this model lacks the gating mechanism of CGEA-BiLSTM, in the data preprocessing stage, the influence of noise can be partially offset by more stringent smoothing filtering, such as increasing the window of SG filtering to 9 days.

[0271] In practice, the model's output layer is connected to a Sigmoid activation function, outputting probability P. irr To reduce false positives, this embodiment employs a dual-threshold decision-making strategy instead of a single threshold: setting a high threshold Th... high =0.8 and low threshold Th low =0.4. If P irr >Th high If it is 0.4, it is directly determined as irrigation; <P irr If the value is less than 0.8, the system further checks whether the minimum RPDI value within the window is below a certain absolute threshold, such as -0.2. If it is, the system is considered irrigation; otherwise, it is considered non-irrigation. This post-processing logic, as a simple supplement to hard rules, compensates to some extent for the uncertainty of deep learning models on fuzzy samples.

[0272] Example 9: Through specific experimental data and comparative analysis, the beneficial effects of the technical method of the present invention compared with the prior art are verified. The data in this example comes from actual observation experiments in an irrigation area in 2025.

[0273] To objectively evaluate the performance of the spatiotemporal fusion algorithm and the irrigation identification model, the following quantitative evaluation indicators were selected in this embodiment.

[0274] For evaluating the accuracy of spatiotemporal fusion, the root mean square error (RMSE), correlation coefficient (R), and structural similarity index (SSIM) are used. The formula for calculating RMSE is:

[0275] RMSE=sqrt(∑((y pred,i -y true,i ) 2 ) / N);

[0276] Among them, y pred,i To fuse the predicted reflectance of the i-th pixel in the image, y true,i The RMSE is the observed reflectance of the Sentinel-2 real image at the corresponding time, and N is the total number of pixels. The smaller the RMSE, the smaller the fusion error.

[0277] The correlation coefficient R measures the degree of linear correlation between predicted and actual values, and its value ranges from -1 to 1. The closer the value is to 1, the stronger the linear correlation. The calculation formula is as follows:

[0278] ;

[0279] in, To predict the local mean of reflectivity, This represents the local mean of the observed reflectance.

[0280] The accuracy of irrigation identification is evaluated using overall accuracy (OA) and F1 score (F1-score). The formula for calculating OA is:

[0281] OA = (TP + TN) / (TP + TN + FP + FN);

[0282] Among them, TP represents true positives, which is the number of correctly identified irrigated samples; TN represents true negatives, which is the number of correctly identified non-irrigated samples; FP represents false positives, which are false alarms; and FN represents false negatives, which are false misses.

[0283] The formula for calculating the F1 score is:

[0284] F1=2*(Precision*Recall) / (Precision+Recall);

[0285] Where Precision is the precision rate and Recall is the recall rate, Precision = TP / (TP+FP) and Recall = TP / (TP+FN).

[0286] The F1 score, which takes into account both precision and recall, is a key indicator for evaluating the performance of imbalanced class classification.

[0287] The experiment selected a time window within the irrigation district that included significant irrigation events, such as July 20 to August 10. During this period, a paddy field underwent flood irrigation.

[0288] In fused image sequences generated using the traditional ESTARFM algorithm, although the RPDI value decreases on the day irrigation occurs, such as July 25th, the decrease is gradual, and the recovery process is smoothly stretched, exhibiting a broad U-shaped characteristic. This is because ESTARFM relies excessively on linear interpolation of surrounding unirrigated pixels and historical images.

[0289] In the fused image sequence generated by the IU-ESTARFM algorithm of this invention, the RPDI value on July 25th showed a precipitous drop and remained low the following day, exhibiting a sharp V-shaped feature, which is consistent with the changing trend of the actual Sentinel-2 observation data.

[0290] The experiment constructed a test set containing 120 samples, covering various meteorological conditions such as sunny, cloudy, and rainy days. The basic Bi-LSTM model performed well on sunny day samples, but it was prone to generating false alarms on days with severe cloud shadow interference or intermittent rainfall, mistaking cloud shadows for irrigation, or having difficulty distinguishing between natural wetting and artificial irrigation on rainy days.

[0291] The improved CGEA-BiLSTM model introduces a fusion confidence level q. t Gating was implemented to automatically suppress the weight of low-quality data; rainfall data P was used as a reference. t Masking is applied to effectively eliminate rainfall interference; event scoring is used. t Attention alignment was performed to accurately pinpoint the irrigation start date. The accuracy of the identification results was improved to 87.5%.

[0292] The method proposed in this invention solves the signal fidelity problem in data fusion through mutation protection and addresses the noise propagation problem in the identification process through confidence level cascading. This method can not only generate higher-quality, high-spatiotemporal-resolution data, but also achieve high-precision daily dynamic monitoring of irrigation under complex weather conditions, demonstrating practical value and promising prospects for wider application.

[0293] This application introduces a mutation protection mechanism in the fusion stage. By calculating the wetting intensity of low spatial resolution images and combining it with rainfall data, a mutation protection factor is constructed. This factor can identify anomalous wetting caused by non-rainfall and forcibly amplifies the contribution of this mutation signal during weight generation. It preserves the V-shaped feature of the sharp drop in surface reflectivity when irrigation occurs, avoiding the smoothing of this as noise by traditional algorithms, and ensuring the fidelity of the irrigation signal in spatiotemporal reconstruction. This solves the problem of smoothing irrigation mutation signals in spatiotemporal fusion.

[0294] This application establishes a confidence-gating and event alignment mechanism. The fusion process outputs a fusion confidence score based on weighted entropy. The recognition model uses this confidence score and rainfall data to construct gating coefficients, automatically suppressing the feature responses of low-quality interpolated data or cloud shadow noise, thus blocking error propagation. Simultaneously, a physics-driven attention mechanism is used to focus the model on abrupt changes with high confidence. This reduces the misclassification rate under complex meteorological conditions and improves attribution accuracy. It also addresses the problem of recognition models being susceptible to fusion noise interference.

[0295] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A remote sensing dynamic monitoring method for farmland irrigation process supported by spatiotemporal data fusion, characterized in that, include: Acquire high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences of the target area during the monitoring period, and unify the spatiotemporal reference of the two; The fusion weight is calculated based on the spatiotemporal correlation of high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences, and a daily high spatiotemporal resolution fused image sequence covering the monitoring period is generated accordingly. Based on the calculation of spectral indices sensitive to moisture changes using daily high spatiotemporal resolution fused image sequences, a temporal feature set is constructed. The time-series feature set is sliced ​​based on a sliding window. The sliced ​​time-series features are then input into a pre-built irrigation identification model. The time-series evolution features within the window are analyzed to determine whether an irrigation event has occurred at the spatiotemporal location corresponding to the window, thus obtaining a dynamic distribution map of the irrigation area.

2. The method according to claim 1, characterized in that, Acquire high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences of the target area during the monitoring period, including: MODIS surface reflectance products were acquired as a high temporal-low spatial resolution image sequence, and Sentinel-2 images were acquired as a low temporal-high spatial resolution image sequence. The band intersection of a high temporal-low spatial resolution image sequence and a low temporal-high spatial resolution image sequence is selected, and the high temporal-low spatial resolution image sequence is resampled to the same spatial resolution as the low temporal-high spatial resolution image sequence. Cloud masking is applied to low temporal-high spatial resolution image sequences, and spatiotemporal interpolation methods are used to fill in the invalid value regions covered by the mask.

3. The method according to claim 1, characterized in that, Spectral indices include the Recorrected Vertical Aridity Index (RPDI). The construction of the time-series feature set is specifically based on the calculation of daily RPDI sequences using the following formula: ； Where Red represents the reflectivity in the red light band, m t The slope of the soil line is corrected for multi-temporal data; NIR represents the near-infrared reflectance; F... v For vegetation cover, Red v NIR represents the reflectance of vegetation in the red light band. v The reflectance of vegetation in the near-infrared band.

4. The method according to claim 1, characterized in that, The fusion weights are calculated based on the spatiotemporal correlation between high-temporal-low-spatial-resolution image sequences and low-temporal-high-spatial-resolution image sequences, including: A search window is constructed centered on the target pixel to be predicted, and similar pixels whose spectral features differ from the target pixel by less than a preset threshold are selected within the search window. For each similar pixel, calculate the spectral similarity distance, spatial distance, and temporal distance between it and the target pixel. Based on spectral similarity distance, spatial distance, and temporal distance, the weight contribution capability of each similar pixel is determined and used to calculate normalized weights to generate fused images.

5. The method according to claim 4, characterized in that, Generate daily high spatiotemporal resolution fused image sequences covering the monitoring period, including: By utilizing the relationship between the changes of high temporal-low spatial resolution image sequences and low temporal-high spatial resolution image sequences at the reference time and the prediction time, the conversion coefficient is calculated for each similar pixel. The comprehensive distance is calculated using the reciprocal or negative exponential relationship between spectral similarity distance, spatial distance, and temporal distance, and then normalized to obtain the normalized weight. Based on normalized weights and conversion coefficients, the reflectance of all similar pixels is weighted and summed to obtain the fused reflectance value of the target pixel at the monitoring time.

6. The method according to claim 4, characterized in that, Generating daily high spatiotemporal resolution fused image sequences covering the monitoring period also includes calculating the fusion confidence level: Obtain daily rainfall data during the monitoring period; The daily moisture index is calculated based on high temporal-low spatial resolution image sequences, and the intensity of the abrupt change in moisture intensity on the target date is calculated based on the difference in moisture index between adjacent dates. A mutation protection factor is calculated based on rainfall data and the intensity of abrupt changes in humidity to indicate the degree of humidity caused by non-rainfalling factors. The initial weights determined based on spectral similarity distance, spatial distance, and temporal distance are corrected using a mutation protection factor to obtain the corrected normalized weights. The reflectance of similar pixels is weighted and summed using the corrected normalized weights to obtain the fused image sequence, and the fusion confidence is calculated based on the distribution characteristics of the corrected normalized weights.

7. The method according to claim 6, characterized in that, The mutation protection factor is calculated based on rainfall data and the intensity of wet change abruptly, according to the following rules: If the rainfall data for the target date is greater than the preset rainfall threshold, or the intensity of the wet change mutation is zero, then the mutation protection factor will be set to the baseline value. If the rainfall data for the target date is less than or equal to the rainfall threshold, and the wet change mutation intensity is greater than zero, then the baseline value is amplified using the wet change mutation intensity to obtain the mutation protection factor. The greater the wet change mutation intensity, the larger the value of the mutation protection factor.

8. The method according to claim 6, characterized in that, The fusion confidence score is calculated based on the distribution characteristics of the corrected normalized weights, using the following formula: A tp =-∑ v∈Ω (W v *ln(W v )); q tp =1-(H tp / ln(N)),0≤q tp ≤1; Among them, H tp Let q be the information entropy of the weight distribution. tp W represents the fusion confidence level for the target date. v Let Ω be the corrected normalized weight of similar pixel v, ln be the set of similar pixels, N be the natural logarithm, and N be the total number of pixels in the set of similar pixels.

9. The method according to claim 1, characterized in that, The irrigation identification model employs a bidirectional long short-term memory network structure, slicing the temporal feature set based on a sliding window. The input to the pre-built irrigation identification model includes: The time-series feature sequences within the sliding window are normalized and then mapped into high-dimensional feature vectors. By using a bidirectional long short-term memory network, the dependencies of temporal feature sequences in the forward and backward time dimensions are extracted to generate hidden state sequences. The hidden state sequence is feature-gathered, and the irrigation probability corresponding to the sliding window is output through a fully connected layer. The time-series feature set includes the Recorrected Vertical Drought Index (RPDI) sequence; The pre-built irrigation identification model is input by directly using the smoothed RPDI sequence as input, and a bidirectional long short-term memory network is used to capture the V-shaped change features of the RPDI sequence before and after irrigation.

10. The method according to claim 6, characterized in that, The irrigation identification model is a confidence-gated and event-aligned model. The input to the pre-built irrigation identification model includes: Obtain the daily RPDI sequence, fused confidence sequence, and rainfall data within the sliding window, and construct the daily composite input vector; Before inputting the composite input vector into the bidirectional long short-term memory network, the input features are gated and weighted using the fusion confidence to suppress the feature responses at low confidence time steps; When performing feature aggregation on the hidden state sequence, an attention mechanism is introduced. The event alignment score, calculated based on the RPDI decrease magnitude and fusion confidence, is used to calculate the attention weight, so that the model focuses on the non-rainfall wet change abrupt moment with high confidence.

Images