A cotton recognition method

Abstract

This invention discloses a cotton identification method, comprising the following steps: extracting Sentinel-2, MODIS, and CHIRPS time-series data of cotton growing season in the target area; constructing pixel-by-pixel confidence coefficients and confidence observation marker masks based on the CHIRPS time-series data to identify low-reliability observations affected by rainfall; using the MODIS normalized vegetation index time-series to guide the reconstruction of the original spectral index time-series of the Sentinel-2 path, and mapping the MODIS path's normalized vegetation index time-series to the value domain space of the reconstructed spectral index time-series of the Sentinel-2 path to obtain a dual-source spectral index time-series; extracting key phenological windows for cotton, calculating the dual-path cumulative spectral phenological index, and performing weighted fusion and residual correction to obtain an enhanced cotton spectral phenological index; finally, inputting it into a classification and identification model to output the spatial distribution results of cotton in the target area. This invention can reduce the impact of cloud and rain interference on early cotton identification and leverage the complementary advantages of multi-source remote sensing data in terms of spatial detail and temporal continuity.

Description

Technical Field

[0001] This invention belongs to the field of agricultural remote sensing and information technology, and more specifically, relates to a cotton identification method. Background Technology

[0002] As an important economic crop, cotton's planting area and spatial distribution information are crucial for cotton production management, yield assessment, and regional agricultural resource allocation. During the cotton's entire growth cycle, the early growing season, from seedling emergence to budding, is a critical stage where vegetative and reproductive growth occur simultaneously. The spatial distribution information during this period directly affects the timeliness and accuracy of subsequent field management measures such as topping, irrigation, and fertilization. Therefore, how to obtain timely and accurate cotton planting area distribution results in the early stages of cotton growth is a critical problem that urgently needs to be solved in the field of agricultural remote sensing monitoring. In recent years, early crop mapping methods based on remote sensing imagery have been applied to crops such as potatoes, winter wheat, and rice, indicating that using remote sensing temporal features for early growing season crop identification has good application potential.

[0003] However, compared to some food crops, early remote sensing mapping of cotton remains significantly challenging. Existing cotton identification methods largely rely on significant spectral features at maturity or boll opening stages. For example, they utilize the white canopy of cotton during boll opening to construct identification indicators such as the White Boll Index (WBI) and the Cotton Boll Opening Index (COBI). While these methods can achieve identification in the later stages of cotton growth, the identification time is significantly delayed, making it difficult to support early field management decisions. Some pre-harvest identification methods, although capable of early cotton identification, generally suffer from issues such as threshold dependence on experience, insufficient cross-regional migration ability, weak inter-year stability, and the inability of fixed phenological windows to adapt to large-scale spatial heterogeneity. With the development of remote sensing technology, multi-source remote sensing data fusion, machine learning, and deep learning methods are increasingly being applied to cotton mapping, improving classification accuracy by fusing multi-temporal, multi-band, or multi-source remote sensing features. However, in the early stages of cotton growth, vegetation cover is low, spectral signals are relatively weak, soil background has a significant impact, and frequent cloud and rain weather can lead to missing remote sensing observations, increased noise, and instability in the temporal characteristics of key phenological periods, making it difficult to guarantee the accuracy and stability of early cotton identification.

[0004] Existing patent application CN117132897A proposes an automatic cotton sample extraction method based on Landsat and MODIS. This method uses Landsat remote sensing images to extract farmland information, and constructs cotton sample point screening rules based on cotton phenological characteristics using MODIS-NDVI time series data. Then, it uses the screened cotton sample points to train a classifier to obtain long-term spatial distribution data of cotton. Existing patent application CN112395914A proposes a plot crop identification method that integrates remote sensing image time series and texture features. This method constructs NDVI time series curves through spectral normalization of multi-source, multi-temporal remote sensing images, and combines plot texture features and a classifier to achieve crop identification. Although the aforementioned existing technologies involve multi-source remote sensing data, MODIS-NDVI time series, utilization of phenological features, inter-sensor spectral normalization, and crop classification, they are still mainly limited to sample screening, time series curve construction, spectral normalization, or classification feature fusion. They have not yet explicitly identified and reconstructed low-confidence observations caused by frequent cloud and rain in the early stages of cotton planting. Therefore, it is difficult to stably express the early spectral phenological features of cotton in the early cotton mapping scenario with insufficient effective observations and strong spatial heterogeneity. Summary of the Invention

[0005] This invention provides a cotton identification method to overcome the problems of existing cotton early remote sensing identification methods being susceptible to cloud and rain weather and the failure of multi-source remote sensing data to fully utilize the complementary advantages of spatial details and temporal continuity.

[0006] The primary objective of this invention is to solve the aforementioned technical problems. The technical solution of this invention is as follows: This invention provides a cotton identification method, comprising the following steps: S1: Extract time series data corresponding to high spatial resolution optical remote sensing data, high temporal resolution remote sensing data and precipitation monitoring data during the cotton growing season in the target area, respectively. S2: Construct pixel-by-pixel confidence coefficients and confidence observation marker masks based on precipitation monitoring time series data; S3: Combining the confidence coefficient and the confidence observation marker mask, the normalized vegetation index time series corresponding to the high temporal resolution remote sensing data is used to guide the reconstruction of the original spectral index time series of the high spatial resolution optical remote sensing data path; and the normalized vegetation index time series of the high temporal resolution remote sensing data path is mapped to the value domain space of the reconstructed spectral index time series of the high spatial resolution optical remote sensing data path to obtain the dual-source spectral index time series; S4: Extract key phenological windows for cotton using the reconstructed spectral index time series of high spatial resolution optical remote sensing data paths; S5: Using the dual-source spectral index time series and the key phenological window of cotton, the cumulative spectral phenological index of the two paths is calculated respectively, and the cumulative spectral phenological index of the two paths is weighted and fused and the residual is corrected to obtain the enhanced cotton spectral phenological index. S6: Input the enhanced cotton spectral phenological index into the classification and recognition model to classify and identify cotton, and output the spatial distribution results of cotton in the target area.

[0007] Compared with the prior art, the beneficial effects of the technical solution of the present invention are: This invention extracts time-series data corresponding to high spatial resolution optical remote sensing data, high temporal resolution remote sensing data, and precipitation monitoring data during the cotton growing season in the target area. Based on the precipitation monitoring time-series data, it constructs pixel-by-pixel confidence coefficients and confidence observation marker masks, allowing different synthesis periods with varying degrees of rainfall interference to have different confidence expressions. This reduces the interference of low-reliability observations on the early spectral time series of cotton before spectral index time series reconstruction. Simultaneously, this invention combines the confidence coefficients and confidence observation marker masks, utilizing the normalized vegetation index time series corresponding to the high temporal resolution remote sensing data to guide the reconstruction of the original spectral index time series of the high spatial resolution optical remote sensing data path. The normalized vegetation index time series under the high temporal resolution remote sensing data path is then mapped to the value domain space of the reconstructed spectral index time series of the high spatial resolution optical remote sensing data path, forming a dual-source spectral index time series. This approach, while preserving the spatial details of the high spatial resolution optical remote sensing data, introduces the temporal continuity advantage of the high temporal resolution remote sensing data, alleviating the problems of discontinuous and unstable spectral index time series in the early cotton growing season caused by cloud and rain interference and insufficient effective observations. Attached Figure Description

[0008] To make the objectives and technical solutions of this invention clearer, the following drawings are provided and described: Figure 1 A flowchart of the method provided in an embodiment of the present invention; Figure 2 Example spatial distribution diagram of classification results in 2020 for three typical cotton-growing areas: Aksu, Kashgar, and Changji, provided for embodiments of the present invention; Figure 3 The results of long-term supervised classification mapping for three study areas from 2019 to 2024 provided for embodiments of the present invention. Detailed Implementation

[0009] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0010] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0011] Example 1: This invention provides a cotton identification method, such as... Figure 1 The diagram shown is a flowchart of a cotton identification method provided by the present invention. The specific steps are as follows: S1: Extract the time series data corresponding to the high spatial resolution optical remote sensing data (Sentinel-2), high temporal resolution remote sensing data (MODIS), and precipitation monitoring data (CHIRPS) during the cotton growing season in the target area, respectively; S2: Based on the precipitation monitoring time series data, the cumulative rainfall and the maximum daily rainfall for each synthesis period are calculated pixel by pixel year by year. The annual adaptive elimination threshold is determined according to the preset physical lower limit and the cumulative rainfall quantile. The rainstorm detection threshold is determined according to the daily heavy rainfall event. Then, the confidence coefficient and confidence observation mark mask are constructed pixel by pixel. S3: Combining the confidence coefficient and the confidence observation marker mask, the original spectral index time series of the Sentinel-2 path is reconstructed using the Normalized Difference Vegetation Index (MODIS NDVI) time series corresponding to the high temporal resolution remote sensing data, resulting in a reconstructed spectral index time series. Under the MODIS path, NDVI is used as a spectral index substitute time series, and its value range is mapped to the value range space of the reconstructed spectral index time series of the Sentinel-2 path using linear regression, resulting in an aligned spectral index time series. The reconstructed spectral index time series and the aligned spectral index time series are used as a dual-source spectral index time series with unified physical and dimensional properties. S4: Using the temporal sequence of reconstructed spectral indices from the Sentinel-2 path, key phenological windows for cotton are extracted pixel by pixel in an adaptive manner. S5: Using the dual-source spectral index time series and the key phenological window of cotton, the cumulative spectral phenological index (CSP) of the two paths is calculated respectively, and the cumulative spectral phenological index of the two paths is weighted and fused and the residual is corrected to obtain the enhanced cotton spectral-phenological index (ECSPI). S6: Input ECSPI into the classification and recognition model to classify and recognize cotton, and output the spatial distribution results of cotton in the target area.

[0012] This invention extracts time-series data corresponding to Sentinel-2, MODIS, and CHIRPS data during the cotton growing season in the target area, respectively. Based on the CHIRPS time-series data, it constructs pixel-by-pixel confidence coefficients and confidence observation marker masks, allowing synthesis periods with varying degrees of rainfall interference to have different confidence expressions. This reduces the interference of low-reliability observations on the early spectral time series of cotton before spectral index time-series reconstruction. Simultaneously, this invention combines the confidence coefficients and confidence observation marker masks, utilizing the normalized vegetation index time series corresponding to the MODIS data to guide the reconstruction of the original spectral index time series of the Sentinel-2 path. It then maps the normalized vegetation index time series under the MODIS path to the spectral index space, forming a dual-source spectral index time series. Thus, while preserving the spatial details of the Sentinel-2 data, it introduces the temporal continuity of the MODIS data. This invention offers several advantages, alleviating the problems of discontinuous and unstable spectral index timing caused by cloud and rain interference and insufficient effective observation in the early stages of cotton growth. Furthermore, this invention utilizes the reconstructed spectral index timing of the Sentinel-2 path to extract key phenological windows for cotton, rather than using fixed phenological windows. This allows for adaptation to differences in growth rhythms among different pixels within the target area, mitigating the difficulty of fixed phenological windows adapting to large-area spatial heterogeneity. Moreover, this invention uses dual-source spectral index timing and key phenological windows to calculate dual-path cumulative spectral phenological indices and performs weighted fusion to obtain enhanced cotton spectral phenological indices. These indices integrate spatial details, temporal continuity, and spectral variation characteristics within key phenological windows, improving the stability and applicability of early cotton identification features. Finally, through a classification and identification model, the spatial distribution results of cotton in the target area are output, providing spatial information support for early cotton identification and subsequent field management decisions in the target area.

[0013] Example 2: This embodiment further supplements the description of the cotton identification method proposed in Embodiment 1.

[0014] The Sentinel-2 data is Sentinel-2 MSI L2A data, the high temporal resolution remote sensing data is MOD13Q1 V6.1 data, and the precipitation monitoring data is CHIRPS Daily data; The extraction of time series data in step S1 specifically includes the following steps: The Sentinel-2 MSI L2A data was processed with dual cloud masking using the QA60 and SCL bands to filter low cloud cover images. The low spatial resolution bands in the filtered Sentinel-2 MSI L2A data were resampled to the spatial resolution corresponding to the high spatial resolution bands. The effective observations during the cotton growing season were synthesized using a preset time interval as the synthesis period to obtain the Sentinel-2 time series data. The MOD13Q1 V6.1 data was processed to remove cloud and snow pixels using the SummaryQA band; the MOD13Q1 V6.1 data after removing cloud and snow pixels was upscaled to the spatial resolution grid corresponding to the high spatial resolution optical remote sensing time series data through bilinear interpolation, and aligned with the time nodes of the Sentinel-2 time series data to obtain MODIS time series data. The time range data corresponding to the cotton growing season in the CHIRPS Daily precipitation data are retained, and the data are accumulated and synthesized according to the same time synthesis period as the Sentinel-2 time series data to obtain the CHIRPS time series data.

[0015] The construction of the confidence coefficient and the confidence observation marker mask in step S2 includes the following steps: S2.1: Based on CHIRPS time series data, calculate the cumulative rainfall and maximum daily rainfall during each synthetic period of the cotton growing season using daily precipitation data; S2.2: Statistical characteristics of cumulative rainfall calculated pixel by pixel year by year, wherein the statistical characteristics include at least quartiles. , , and ; (Set to 10mm in this embodiment) The annual adaptive rejection threshold is determined by a dynamic threshold based on the statistical characteristics of cumulative rainfall, as shown in the following expression: = max( , ) At the same time, the rainstorm detection threshold based on single-day heavy rainfall events was determined. In this embodiment, it is set to 20mm.

[0016] S2.3: The confidence coefficient for the t-th synthesis period. Calculated using the following piecewise function:

[0017] In the formula, To calculate the cumulative rainfall, the confidence coefficients for each composite period are... As confidence markers for each pixel, a confidence observation marker mask with reliability gradient information is generated, wherein... = 0 indicates completely reliable. = 0.95 indicates completely unreliable, while intermediate values ​​reflect different levels of observation reliability.

[0018] This invention extracts time-series data corresponding to high spatial resolution optical remote sensing data, high temporal resolution remote sensing data, and precipitation monitoring data during the cotton growing season in the target area. Based on the precipitation monitoring time-series data, it constructs pixel-by-pixel confidence coefficients and confidence observation marker masks. This allows for reliability labeling of observation periods affected by rainfall before subsequent spectral index time-series reconstruction. Therefore, this invention does not directly process all remote sensing observation data uniformly. Instead, it first uses precipitation monitoring information to constrain the observation quality, resulting in different confidence expressions for synthesis periods with varying degrees of rainfall interference. This reduces the interference of low-reliability observations on the early spectral time series of cotton, improving the reliability of the basic data for subsequent time-series reconstruction and phenological feature extraction.

[0019] The guided reconstruction of the original spectral index time series of the high spatial resolution optical remote sensing data path in step S3 includes the following steps: S3.1: Constructing the original spectral index time series of the Sentinel-2 path The expression is as follows:

[0020] In the formula, , , These represent the surface reflectance at time t in the red-edge 2 band, red-edge 3 band, and near-infrared band, respectively. S3.2: Using the confidence coefficient and the confidence observation marker mask, determine the low-confidence observation period and the high-confidence observation period in the high spatial resolution optical remote sensing data path. For the t-th composite period, if the cumulative rainfall during the composite period... Greater than or equal to the annual adaptive removal threshold Or, the maximum daily rainfall during the synthesis period is greater than or equal to the rainstorm detection threshold. The composite period is defined as the low-confidence observation period; if the cumulative rainfall during the composite period is... Less than the annual adaptive rejection threshold Furthermore, the maximum daily rainfall was less than the rainstorm detection threshold. The synthesis period is defined as the high-confidence observation period; S3.3: For each low-confidence observation period, search for the nearest high-confidence observation period as an anchor point within a preset number (6 in this embodiment) of composite periods before and after it. Guided values ​​are calculated for three cases: available anchor points on both sides, available anchor points on only one side, and no reliable anchor points. The guide weight of NDVI is set according to the extent to which it exceeds a threshold; that is, the more severe the rainfall interference, the more it relies on the more stable curve shape of MODIS NDVI. S3.3.1: If both anchor points are available, a hybrid interpolation strategy combining MODIS NDVI morphological weights and time-linear weights is adopted. The expression for the hybrid weight w(t) is as follows:

[0021] In the formula, , and Represent the current synthesis period t to be reconstructed and the forward reliable anchor point, respectively. and backward trusted anchors Normalized vegetation index values ​​corresponding to high temporal resolution remote sensing data; , and Represent the current synthesis period t to be reconstructed and the forward reliable anchor point, respectively. and backward trusted anchors The corresponding day sequence, , As weighting coefficients, this embodiment uses 0.6 and 0.4 respectively, taking into account the stability of MODIS NDVI phenological information and time linear interpolation, and then obtaining the guiding value through linear interpolation between the two anchor points. The expression is as follows:

[0022] S3.3.2: If only one anchor point is available, extrapolation is performed using MODIS NDVI ratio scaling. It is assumed that the spectral index SI and NDVI change approximately proportionally during the same growth stage, and the guiding value... The expression is as follows:

[0023] In the formula, For single-sided available anchor points corresponding to the synthesis period; to prevent extreme values, the ratio scaling factor... Limit it to the range [0.5, 2.0]; S3.3.3: If no reliable anchor point is available, no guiding value is calculated for the synthesis period, and its missing state is maintained. The information of the synthesis period is compensated by the adaptive weighted fusion mechanism of Sentinel-2 path and MODIS path in subsequent step S5, so as to avoid the introduction of false signals due to missing observations or low confidence data. S3.4: Utilizing guiding values , and confidence coefficient The reconstructed spectral index time series of high spatial resolution optical remote sensing data paths, which takes into account both spatial detail and phenological continuity, was calculated. The expression is as follows: .

[0025] Step S3, which maps the normalized vegetation index time series of high temporal resolution remote sensing data paths to their value range, includes the following steps: The normalized vegetation index (NDVI) time series was obtained using the high temporal resolution remote sensing time series data. Since the MODIS sensor lacks a red edge band, the spectral index cannot be directly calculated. Moreover, the NDVI can effectively reflect the growth status of the vegetation canopy and is similar to the spectral index in a physical sense. Therefore, this time series was used as a substitute time series for the spectral index of the MODIS path. Using high-confidence observation samples of the same pixel and the same time phase, pixel-level linear regression is performed by replacing the time series of the reconstructed spectral index of the high spatial resolution optical remote sensing data path with the spectral index of the MODIS path, and the first regression coefficient a and the second regression coefficient b of the corresponding pixel are obtained. Using the first regression coefficient 'a' and the second regression coefficient 'b', a linear mapping is performed on the spectral index substitution time series of the MODIS path to obtain the range-aligned spectral index time series of the MODIS path. The calculation formula is:

[0026] In the formula, This represents the spectral index substitution value of the MODIS path at time t. This represents the spectral index value of the MODIS path after range alignment at time t.

[0027] This invention combines confidence coefficients and confidence observation marker masks. Under a high spatial resolution optical remote sensing data path, it utilizes the normalized vegetation index (NVI) time series corresponding to high temporal resolution remote sensing data to guide the reconstruction of the original spectral index time series of the high spatial resolution optical remote sensing data path, resulting in a reconstructed spectral index time series. Simultaneously, under a high temporal resolution remote sensing data path, the NVI time series is used as a substitute time series for the spectral index, and its value range is mapped to the spectral index space, resulting in an aligned spectral index time series. Through the above processing, this invention can preserve the spatial details of high spatial resolution optical remote sensing data while introducing the temporal continuity advantage of high temporal resolution remote sensing data, forming a dual-source spectral index time series with relatively unified physical meaning and dimensions. This alleviates the problem of discontinuous and unstable spectral index time series in cotton early stages due to cloud and rain interference and insufficient effective observation.

[0028] Extracting the key phenological window for cotton in step S4 includes the following steps: S4.1: The reconstructed spectral index time series of the Sentinel-2 path is filtered and smoothed using the second-order Savitzky-Golay smoothing method to suppress residual noise while preserving phenological transition characteristics, thus obtaining the smoothed time series. S4.2: Based on the smoothed time sequence, detect the day sequence corresponding to the maximum value of the spectral index during the cotton growing season pixel by pixel, and determine the day sequence corresponding to the maximum value of the spectral index as the termination point of the key phenological window for cotton. ; S4.3: At the termination point Within the previously preset sowing to seedling stage time range (set to DOY∈[90,180] in this embodiment), the smoothed time series is used to calculate the second-order difference, and the date corresponding to the maximum value of the second-order difference is determined as the starting point of the key phenological window for cotton. ; S4.4: For the starting point and termination point Apply validity constraints when the termination point With the starting point When the time interval between the two is less than p (3 in this embodiment) composite periods, pixels that exceed the preset range are reverted to the values ​​obtained in advance from cotton samples within the study area. , Mean.

[0029] The method for constructing the enhanced cotton spectral phenological index in step S5 includes the following steps: S5.1: Utilizing the timing of dual-source spectral indexes and the starting point of key phenological windows for cotton. and termination point The cumulative spectral phenological index of the high spatial resolution optical remote sensing data path was calculated respectively. Cumulative spectral phenological index of high temporal resolution remote sensing data path Calculated according to the following expression:

[0030]

[0031]

[0032]

[0033] In the formula, , These represent the reconstructed spectral index time series of the high spatial resolution optical remote sensing data path and the spectral index time series of the high temporal resolution remote sensing data path, respectively, at time t. , This represents the growth rate within the corresponding spectral index time-series phenological window. , Indicates the end point The corresponding spectral index value, Indicates the starting point The corresponding spectral index value, where n represents the number of synthesis periods within the phenological window; S5.2: with The distribution of is used as the target value range for . Z-score normalization and range alignment are performed to obtain the cumulative spectral phenological index after range alignment. ; S5.3: Construct pixel-wise adaptive fusion weights based on the temporal gap severity score, local spatial consistency score, and dual-path consistency score within the phenological window. ,include: S5.3.1: with to Based on the number of valid Sentinel-2 observations within the window, a temporal gap severity score is obtained in the interval [0, 1] after transformation using the Sigmoid function. The expression is as follows:

[0034] In the formula, To normalize the observation density, k is the steepness parameter of the Sigmoid function (k = 8 in this embodiment). The inflection point of the Sigmoid function (in this embodiment) = 0.4).

[0035] S5.3.2: Evaluate each cell The degree of deviation of the value from the mean of its 3×3 neighborhood is used to obtain the local spatial consistency score. The expression is as follows:

[0036] In the formula, x and y are the horizontal and vertical coordinates, respectively. and Within a 3×3 neighborhood The mean and standard deviation.

[0037] S5.3.3: Calculate the value of each cell within the window. , The difference between the mean of the two paths and the mean of its neighborhood is recorded as 1 if they have the same sign, and 0 otherwise. Then, the mean of the two paths is taken in the neighborhood to obtain the bipath consistency score. This is used to examine whether two data sources exhibit a consistent spatial variation trend within a local area.

[0038] S5.3.4: The temporal gap severity score, local spatial consistency score, and dual-path consistency score are weighted and summed according to a preset ratio to obtain the pixel-by-pixel adaptive fusion weight. The expression is as follows:

[0039] In the formula, , , The weighting coefficients for the three indicators are 0.70, 0.15, and 0.15 respectively in this embodiment. The time series gap is given the highest weight as the core indicator, and local spatial consistency and dual-path consistency are used as auxiliary constraints.

[0040] S5.4: Utilizing the aforementioned adaptive fusion weights right Cumulative spectral phenological index of high temporal resolution remote sensing data path after value range alignment Perform weighted fusion to obtain the weighted fusion result. The expression is as follows:

[0041] S5.5: Perform residual bias correction on the weighted fusion results, and calculate the fusion results with... The mean residual Δ between the two is used to subtract nΔ (n is 0.5 in this embodiment) from the fusion result to eliminate systematic bias, thus obtaining the enhanced cotton spectral phenology index ECSPI.

[0042] This invention further utilizes the reconstructed spectral index time series of high spatial resolution optical remote sensing data paths to extract key phenological windows for cotton. Using the dual-source spectral index time series and the key phenological windows, dual-path cumulative spectral phenological indices are calculated separately. These dual-path cumulative spectral phenological indices are then weighted and fused to obtain an enhanced cotton spectral phenological index. Since the key phenological windows are extracted based on the reconstructed spectral index time series, rather than using a fixed phenological window, it can adapt to the differences in growth rhythms among different pixels within the target area, alleviating the problem that fixed phenological windows are difficult to adapt to large-area spatial heterogeneity. Simultaneously, the enhanced cotton spectral phenological index can integrate the spatial details of high spatial resolution optical remote sensing data, the temporal continuity of high temporal resolution remote sensing data, and the spectral variation characteristics within the key phenological windows, improving the stability and applicability of early cotton identification features.

[0043] If the classification and recognition model described in step S6 is an unsupervised classification model, obtaining the unsupervised classification result of cotton includes the following steps: The enhanced cotton spectral phenology index was input into the Otsu-Sauvola adaptive threshold segmentation model. The global optimal segmentation threshold was determined in the entire study area by the Otsu algorithm, and the local threshold was calculated pixel by pixel in the preset sliding window based on the Sauvola algorithm. The final segmentation threshold is determined based on the global threshold and the local threshold. The ECSPI is then binarized using the final segmentation threshold to obtain the cotton / non-cotton binarization result. A circular structuring element with a radius of 1 pixel is applied to the binarized result to perform a closing operation, filling the plot holes and smoothing the boundaries, thus obtaining the cotton unsupervised classification result.

[0044] If the classification and identification model described in step S6 is a supervised classification model, obtaining the cotton supervised classification result includes the following steps: Using the enhanced cotton spectral phenology index (ECSPI) as the core feature, a simplified feature set was constructed by combining the time-series statistics of NDVI and phenological parameters. A machine learning classifier is trained using the simplified feature set; Input the simplified feature set corresponding to the region to be identified into the trained machine learning classifier, and output the cotton supervised classification result.

[0045] This invention inputs an enhanced cotton spectral phenological index into a classification and identification model to classify and identify cotton, outputting the spatial distribution results of cotton in the target area. Since the enhanced cotton spectral phenological index is obtained based on confidence constraints, time-series reconstruction of dual-source spectral indices, extraction of key phenological windows, and fusion of dual-path cumulative spectral phenological indices, it can reduce the impact of low-reliability observations and limitations of single data sources on the classification results, making the cotton spatial distribution results more suitable for early identification of cotton in the target area and subsequent field management decisions.

[0046] After outputting the spatial distribution results of cotton in the target area in step S6, the method further includes accuracy evaluation and feedback correction of the spatial distribution results of cotton in the target area, specifically including the following steps: (1) Using a test set divided by a spatial block strategy, the classification accuracy of the cotton spatial distribution results in the target area is evaluated. The classification accuracy evaluation index includes overall accuracy OA, Kappa coefficient and F1-score. (2) Apply the cotton identification method to remote sensing data from multiple years to obtain the spatial distribution results of cotton in multiple years, and evaluate the multi-year stability of the cotton identification method based on the coefficient of variation (CV) of the interannual classification area. (3) When the classification accuracy evaluation index or the multi-year stability evaluation result does not reach the preset threshold, the relevant parameters in steps S2 to S6 are corrected by feedback, and steps S2 to S6 are re-executed based on the corrected parameters. The feedback correction includes at least one of the following: a) Adjust the local sliding window size, Sauvola correction coefficients, and fusion strategy of global and local thresholds in the Otsu-Sauvola adaptive threshold segmentation model; b) Adjust the preset physical lower limit or cumulative rainfall quantile level in the annual adaptive rejection threshold, and adjust the rainstorm detection threshold; c) Adjust the half-window size, second-order difference search range, minimum growth cycle constraint, and starting point of the second-order Savitzky-Golay smoothing method filter. and termination point The effective scope; d) Adjust the weight coefficients corresponding to the temporal gap severity score, local spatial consistency score, and dual-path consistency score in the adaptive fusion weights, as well as the steepness parameter and inflection point parameter of the Sigmoid function.

[0047] Example 3: This embodiment provides a cotton identification method, the specific steps of which are as follows: Step 1: Temporal extraction of multi-source remote sensing data To construct a stable temporal feature set covering the entire cotton growing season, this embodiment takes three typical cotton-growing areas in Xinjiang—Aksu, Kashgar, and Changji—as the research areas and acquires all available Sentinel-2 MSI L2A data and MODIS MOD13Q1 V6.1 data for the growing season from 2019 to 2024 on the Google Earth Engine (GEE) cloud platform.

[0048] The specific process of temporal extraction from Sentinel-2 data is as follows: First, a dual cloud mask is constructed using the QA60 and SCL bands to screen low-cloud-coverage images during the cotton growing season. Second, the 20m bands (B5-B7, B8A, B11, B12) are resampled to a 10m spatial resolution using nearest neighbor resampling. Finally, with a 10-day composite period, the median of all valid observations within the growing season is used to generate a time-series composite image, thereby minimizing noise while ensuring temporal consistency.

[0049] The specific process of MODIS data time series extraction is as follows: The MOD13Q1 V6.1 product is selected. After removing cloud and snow pixels based on the SummaryQA band, the data is upscaled to a 10m grid through bilinear interpolation and aligned with the Sentinel-2 time node as a time anchor reference for subsequent spatiotemporal fusion.

[0050] The specific process of extracting CHIRPS data time series is as follows: retain the time range data that matches the cotton growing season, and synthesize the daily precipitation of CHIRPS to the same 10-day synthesis cycle as Sentinel-2.

[0051] Step 2: Construct confidence coefficients and confidence observation marker masks To address the sensitivity of spectral indices to rainfall, a CHIRPS daily precipitation product with a spatial resolution of 0.05° is introduced to construct a quality mask for rainfall observations. The specific implementation steps are as follows: First, calculate the cumulative rainfall for each period and the maximum daily rainfall; Subsequently, the quartiles of the cumulative rainfall for all composite periods were statistically analyzed pixel by pixel year by year. / / / Set annual adaptive rejection threshold For max(10mm, The threshold for detecting heavy rainfall is set at 10 mm, where 10 mm is the lower physical limit, and a threshold for detecting heavy rainfall exceeding 20 mm in a single day is also set. Accurately identify pixels affected by heavy rainfall; For the t-th synthesis period, the corresponding confidence coefficient Calculated using the following piecewise function:

[0052] In the formula, To calculate the cumulative rainfall, the confidence coefficients for each composite period are... As confidence markers for each pixel, a confidence observation marker mask with reliability gradient information is generated, wherein... = 0 indicates complete reliability. = 0.95 indicates completely unreliable, while intermediate values ​​reflect different levels of observation reliability.

[0053] Step 3: Time series construction of dual-source spectral index To obtain spectral time series with high spatiotemporal continuity and to accurately locate key phenological windows for cotton, this step constructs and reconstructs the spectral indices of two independent paths, Sentinel-2 and MODIS, based on the low-confidence observation marker mask.

[0054] (1) Sentinel-2 path SI calculation and MODIS-guided reconstruction Sentinel-2 is equipped with three bands that continuously cover the red-edge to near-infrared region: Red Edge 2 (RE2), Red Edge 3 (RE3), and Near Infrared (NIR). These three bands are used together to form the Spectral Index (SI), calculated using the following formula:

[0055] The low-confidence and high-confidence observation periods in the Sentinel-2 path are determined using confidence coefficients and confidence observation marker masks. High-confidence observation periods refer to observation periods where cumulative rainfall is less than the adaptive threshold and daily maximum rainfall is less than the rainstorm threshold. Low-confidence observations do not meet the above conditions. The guiding weight of NDVI is set according to the extent of exceeding the threshold. That is, the more severe the rainfall interference, the more it relies on the more stable curve shape of MODIS NDVI.

[0056] Next, we need to process the data marked as low confidence. The values ​​are reconstructed. This embodiment employs a gradual MODIS-guided reconstruction strategy: using continuous confidence coefficients... Characterize the degree of disbelief in disbelief observation periods to achieve a smooth transition from complete trust in S2 to complete reliance on MODIS guiding values:

[0057] in The degree to which rainfall exceeds a threshold is determined: when rainfall is less than the calculated threshold, = 0, fully trust the original S2 observations; when rainfall exceeds the threshold, = 0.1-0.9, gradually changing according to the excess ratio; when a heavy rainfall event occurs with a maximum daily rainfall of ≥20 mm, = 0.95, almost entirely dependent on MODIS bootloader. The calculation method is as follows: For each synthesis period marked as low confidence, the nearest confidence observation within the range of 6 synthesis periods before and after it is searched as anchor points. Interpolation or extrapolation is performed in three cases: two-sided anchor points are available, only one-sided anchor points are available, and no reliable anchor points are available.

[0058] (2) MODIS path SI calculation and linear regression reconstruction Since the MODIS sensor lacks the red-edge band of Sentinel-2, SI cannot be directly calculated. Therefore, in this embodiment, MODIS NDVI is used as a substitute variable, denoted as:

[0059] because and The dimensions and numerical ranges differ significantly. This embodiment uses pixel-level linear regression to... Mapped to time series reconstruction The range space of the values ​​is adjusted to ensure dimensional consistency between the two. Finally, the regression coefficients are applied to the MODIS NDVI time series to obtain the reconstructed time series with aligned ranges. :

[0060] Step 4: Pixel-by-pixel automated extraction of key phenological windows Based on reconstruction and In this embodiment, the starting point of the key phenological window is further specified. and end point The detection method has been changed to a pixel-by-pixel adaptive approach. The specific implementation involves the following steps: a. SG filtering smoothing: Apply a second-order Savitzky-Golay filter to the applied time sequence.

[0061] b. Detection: Defined as the day sequence (DOY) corresponding to the maximum value of SI during the entire growing season, with the effective range set in DOY∈[160,260].

[0062] c. Detection: Defined as Previously, the DOY corresponding to the second-order difference maximum value in the interval [90, 180].

[0063] d. Validity constraint: requires - ≥3 synthesis periods, otherwise regress to the mean of the study area.

[0064] Step 5: Construction of Enhanced Cotton Spectral Phenological Index (ECSPI) (1) CSP Index Calculation Timing of acquiring dual-source spectral index , After the pixel-by-pixel phenological window, the CSP indexes for the Sentinel-2 and MODIS paths are calculated using the formula, denoted as... and :

[0065]

[0066] (2) Adaptive fusion weight construction Fusion weights of the Sentinel-2 path It is composed of a weighted combination of three components: a. Time gap severity score (70% weighting): based on to The number of valid S2 observations within the window is the basic indicator, which is converted into a gap score in the [0,1] interval by the Sigmoid function.

[0067] b. Local Spatial Consistency Score (15% weight): Evaluates each cell The degree of deviation of the value from the mean of its 3×3 neighborhood.

[0068] c. S2-MODIS Consistency Score (15% weight): Examines whether S2 and MODIS show consistent spatial variation trends in local regions.

[0069] Finally, the three components are weighted and summed according to a preset ratio to obtain the fusion weight:

[0070] (3) Range alignment and weighted fusion because and There may be systematic biases; this embodiment uses Z-score normalization for value range alignment. The distribution of is used as the target value range. A linear transformation is applied to achieve the same mean and variance. Subsequently, weighted fusion is performed using pixel-wise adaptive weights.

[0071] (4) Residual deviation correction After fusion, residual bias correction is still required. The calculation results are then compared with... The mean residual Δ between the two values ​​is used to subtract 0.5Δ from the fusion result to eliminate systematic bias. The fusion index after the above processing is the enhanced crop spectral phenology index (ECSPI) proposed in this invention.

[0072] Step 6: Early identification and accuracy evaluation of cotton To comprehensively evaluate the effectiveness of the method of the present invention, a variety of classification methods were designed for comparative experiments.

[0073] (1) Unsupervised classification: ECSPI + Otsu-Sauvola After obtaining ECSPI feature images, the Otsu-Sauvola two-layer adaptive thresholding method is used to achieve automatic cotton identification. First, the OTSU algorithm is used to determine the global threshold. Secondly, based on the Sauvola algorithm, the local threshold t(x,y) is calculated pixel by pixel within a sliding window of r×r pixels; the maximum value of the local threshold and the global threshold is taken as the final segmentation threshold T(x,y)=max[t(x,y), The final segmentation threshold is determined based on global and local thresholds. This final threshold is then used to perform binarization segmentation on ECSPI, yielding cotton / non-cotton binarized results. Finally, a closing operation is performed on the binarized results using a circular structuring element with a radius of 1 pixel to improve the result quality.

[0074] (2) Supervised classification: Random Forest ECSPI simplified feature classification To verify ECSPI's ability to condense cotton information, a 12-dimensional simplified feature set was constructed, including the ECSPI index value, ECSPI intermediate values, and so on. , , ρ, The data includes four time-series statistics and three phenological parameters of NDVI. A simplified feature set is combined with a random forest classifier to classify cotton-growing areas.

[0075] (3) Accuracy Evaluation The test set was partitioned using a spatial segmentation strategy, and overall accuracy (OA), Kappa coefficient, and F1-score were selected as evaluation metrics. Furthermore, the method was applied to six years from 2019 to 2024, and the stability of the method was assessed using the coefficient of variation (CV) of the inter-year classification area.

[0076] The verification results of the embodiments of the present invention in three typical cotton-growing areas, Aksu, Kashgar, and Changji, are as follows: Figure 2 , Figure 3 As shown in Tables 1 to 3.

[0077] like Figure 2The image shows a spatial distribution comparison of classification results for three typical cotton-growing areas in 2020. From left to right, the image represents the three sub-regions of Aksu, Kashgar, and Changji; from top to bottom, it shows the original remote sensing image, reference data, unsupervised classification results, and supervised classification results. Specifically, the first column represents the Aksu cotton-growing area. Figure 2 (a) is the original remote sensing image of Aksu region. Figure 2 (d) is the corresponding ground reference data (yellow area). Figure 2 (g) and Figure 2 (j) shows the unsupervised and supervised classification results for this region (green area). The second column represents the Kashgar cotton region. Figure 2 (b) is the original remote sensing image of the Kashgar region. Figure 2 (e) is the corresponding ground reference data. Figure 2 (h) and Figure 2 (k) represents the unsupervised and supervised classification extraction results for this region, respectively. The third column represents the Changji cotton region. Figure 2 (c) is the original remote sensing image of the Changji area. Figure 2 (f) represents the corresponding ground reference data. Figure 2 (i) and Figure 2 (l) are the unsupervised and supervised classification extraction results for this region, respectively.

[0078] like Figure 3 The image shows the results of long-term supervised classification mapping from 2019 to 2024. Specifically, Figure 3 (a) ~ (c) are the results of the supervision and classification mapping of cotton-growing areas in Aksu, Kashgar and Changji in 2019 (green areas). Figure 3 (d) ~ (f) Figure 3 (g) ~ (i), Figure 3 (j) ~ (l), Figure 3 (m) ~ (o) and Figure 3 (p) ~ (r) correspond to the supervision classification mapping results of the above three cotton-growing areas in 2020, 2021, 2022, 2023 and 2024, respectively.

[0079] Table 1 summarizes the accuracy evaluation results of supervised and unsupervised classification in three study areas from 2019 to 2021 for the embodiments of the present invention.

[0080] Table 1

[0081] As can be seen from the statistical results in Table 1, the supervised classification method used in this embodiment of the invention exhibits significantly superior mapping accuracy in all study areas and in different years. Specifically, in the Aksu, Kashgar, and Changji study areas, the overall accuracy (OA) and F1 score of supervised classification remained consistently above 96%, and the Kappa coefficient was above 92%. In contrast, the OA of unsupervised classification mainly ranged between 83% and 92%, while the Kappa coefficient fluctuated between 66% and 84%.

[0082] Data comparison shows that the supervised classification method of this invention not only significantly outperforms unsupervised classification in overall accuracy, but also exhibits extremely high stability and robustness across different years and geographical regions, enabling the extraction of long-term cotton planting areas with high accuracy and high reliability.

[0083] Table 2 compares the available time windows of the ECSPI proposed in this invention with those of the existing mainstream methods WBI and COBI in the three major cotton-growing areas of Xinjiang.

[0084] Table 2

[0085] The results showed that ECSPI completed mapping approximately 60–70 days earlier than WBI and COBI, demonstrating a significant time advantage.

[0086] Table 3 shows the interannual variation of cotton planting area in each study area from 2019 to 2024, including the mean, standard deviation, and coefficient of variation (CV), to assess the multi-year stability of the method.

[0087] Table 3

[0088] Statistical results show that the cotton planting area extracted in this embodiment of the invention exhibits extremely high interannual stability and robustness over long time series. As shown in Table 3, the coefficient of variation (CV) of cotton planting area in the three study areas remained at a low level (less than 10%) between 2019 and 2024. Among them, the interannual fluctuation in the Changji cotton area was the most stable, with a CV of only 1.9% over six years and a stable average area of ​​approximately 10360.5. The CVs in the Aksu and Kashgar cotton areas were also controlled at 5.5% and 8.8%, respectively.

[0089] The low coefficient of variation and standard deviation fully demonstrate that the method of the present invention can effectively suppress noise interference caused by interannual phenological differences, illumination conditions and atmospheric environmental changes when conducting long-term remote sensing monitoring, and exhibits excellent multi-year mapping stability and time-series robustness.

[0090] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art can make other variations or modifications based on the above description. It is neither necessary nor possible to exhaustively describe all embodiments here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the claims of the present invention.

Claims

1. A cotton recognition method, characterized in that, Includes the following steps: S1: Extract time series data corresponding to high spatial resolution optical remote sensing data, high temporal resolution remote sensing data and precipitation monitoring data during the cotton growing season in the target area, respectively. S2: Construct pixel-by-pixel confidence coefficients and confidence observation marker masks based on precipitation monitoring time series data; S3: Combining confidence coefficients and confidence observation marker masks, the original spectral index time series of high spatial resolution optical remote sensing data paths is reconstructed using the normalized vegetation index time series corresponding to high temporal resolution remote sensing data. The normalized vegetation index time series of the high temporal resolution remote sensing data path is mapped to the value domain space of the reconstructed spectral index time series of the high spatial resolution optical remote sensing data path to obtain the dual-source spectral index time series. S4: Extract key phenological windows for cotton using the reconstructed spectral index time series of high spatial resolution optical remote sensing data paths; S5: Using the dual-source spectral index time series and the key phenological window of cotton, the cumulative spectral phenological index of the two paths is calculated respectively, and the cumulative spectral phenological index of the two paths is weighted and fused and the residual is corrected to obtain the enhanced cotton spectral phenological index. S6: Input the enhanced cotton spectral phenological index into the classification and recognition model to classify and identify cotton, and output the spatial distribution results of cotton in the target area.

2. The cotton identification method according to claim 1, characterized in that, The extraction of time-series data in step S1 specifically includes the following steps: High spatial resolution optical remote sensing data were processed with dual cloud masking using QA60 and SCL bands to filter low cloud cover images. The low spatial resolution bands in the filtered high spatial resolution optical remote sensing data were resampled to the spatial resolution corresponding to the high spatial resolution bands. The effective observations during the cotton growing season were synthesized using a preset time interval as the synthesis period to obtain high spatial resolution optical remote sensing time series data. Cloud and snow pixel removal was performed on high temporal resolution remote sensing data using the SummaryQA band. The high temporal resolution remote sensing data after cloud and snow pixel removal was then upscaled to the spatial resolution grid corresponding to the high spatial resolution optical remote sensing time series data through bilinear interpolation and aligned with the time nodes of the high spatial resolution optical remote sensing time series data to obtain high temporal resolution remote sensing time series data. The time range data corresponding to the cotton growing season in the daily precipitation data of precipitation monitoring data is retained, and the data is accumulated and synthesized according to the same time synthesis period as the high spatial resolution optical remote sensing time series data to obtain precipitation monitoring time series data.

3. The cotton identification method according to claim 1, characterized in that, The construction of the confidence coefficient and confidence observation marker mask in step S2 includes the following steps: S2.1: Based on precipitation monitoring time series data, calculate the cumulative rainfall and maximum daily rainfall during each synthetic period of the cotton growing season using daily precipitation data; S2.2: Statistical characteristics of cumulative rainfall calculated pixel by pixel year by year, wherein the statistical characteristics include at least quartiles. , , and ; The annual adaptive removal threshold is determined by a dynamic threshold based on the statistical characteristics of cumulative rainfall, as shown in the following expression: = max( , ) At the same time, the rainstorm detection threshold based on single-day heavy rainfall events was determined. ; S2.3: For the t-th synthesis period, the corresponding confidence coefficient Calculated using the following piecewise function: In the formula, To calculate the cumulative rainfall, the confidence coefficients for each composite period are... As confidence markers for each pixel, a confidence observation marker mask is generated.

4. The cotton identification method according to claim 3, characterized in that, Step S3 involves guided reconstruction of the time series of the original spectral index of the high spatial resolution optical remote sensing data path, including the following steps: S3.1: Constructing the original spectral index time series of high spatial resolution optical remote sensing data paths The expression is as follows: In the formula, , , These represent the surface reflectance at time t in the red-edge 2 band, red-edge 3 band, and near-infrared band, respectively. S3.2: Using the confidence coefficient and the confidence observation marker mask, determine the low-confidence observation period and the high-confidence observation period in the high spatial resolution optical remote sensing data path. For the t-th composite period, if the cumulative rainfall during the composite period... Greater than or equal to the annual adaptive removal threshold Or, the maximum daily rainfall during the synthesis period is greater than or equal to the rainstorm detection threshold. The composite period is defined as the low-confidence observation period; if the cumulative rainfall during the composite period is... Less than the annual adaptive rejection threshold Furthermore, the maximum daily rainfall was less than the rainstorm detection threshold. The synthesis period is defined as the high-confidence observation period; S3.3: For each low-confidence observation period, search for the nearest high-confidence observation period as an anchor point within a preset number of composite periods before and after it. Guided value calculation is performed in three cases: when both sides of the anchor point are available, when only one side of the anchor point is available, and when no reliable anchor point is available. S3.3.1: If both anchor points are available, the expression for the mixed weight w(t) is as follows: In the formula, , and Represent the current synthesis period t to be reconstructed and the forward reliable anchor point, respectively. and backward trusted anchors Normalized vegetation index values ​​corresponding to high temporal resolution remote sensing data; , and Represent the current synthesis period t to be reconstructed and the forward reliable anchor point, respectively. and backward trusted anchors The corresponding day sequence, , These are the weighting coefficients; The guiding value is then obtained through linear interpolation between the two anchor points. The expression is as follows: S3.3.2: If only one side of the anchor point is available, guide value The expression is as follows: In the formula, For each available anchor point on one side, the corresponding synthesis period is specified. S3.3.3: If no trusted anchor point is available, no guiding value is calculated for this synthesis period, and the guiding value is left missing. S3.4: Utilizing guiding values , and confidence coefficient The reconstructed spectral index time series of high spatial resolution optical remote sensing data paths was calculated. The expression is as follows: 。 5. The cotton identification method according to claim 1, characterized in that, Step S3 involves mapping the time series of the normalized vegetation index for high temporal resolution remote sensing data paths to their value range, including the following steps: Normalized vegetation index time series is obtained using high temporal resolution remote sensing time series data, and the normalized vegetation index time series is used as the spectral index replacement time series of the high temporal resolution remote sensing data path. By using the reconstructed spectral index time series of high spatial resolution optical remote sensing data paths and the spectral index of high temporal resolution remote sensing data paths to replace the high-confidence observation samples of the same pixel and the same time phase in the time series, a pixel-level linear regression model is established to obtain the first regression coefficient a and the second regression coefficient b of the corresponding pixel. Using the first regression coefficient 'a' and the second regression coefficient 'b', a linear mapping is performed on the spectral index substitution time series of the high temporal resolution remote sensing data path to obtain the spectral index time series of the high temporal resolution remote sensing data path after value range alignment, as shown in the following expression: In the formula, The spectral index substitution value represents the high temporal resolution remote sensing data path corresponding to time t. This represents the spectral index value of the high temporal resolution remote sensing data path corresponding to time t after value range alignment.

6. The cotton identification method according to claim 1, characterized in that, Step S4 involves extracting the key phenological window for cotton, including the following steps: S1: The reconstructed spectral index time series of high spatial resolution optical remote sensing data paths is filtered and smoothed using the second-order Savitsky-Gory smoothing method to suppress residual noise while preserving phenological transition characteristics, resulting in a smoothed time series. S2: Based on the smoothed time sequence, the day sequence corresponding to the maximum value of the spectral index during the cotton growing season is detected pixel by pixel, and the day sequence corresponding to the maximum value of the spectral index is determined as the termination point of the key phenological window for cotton. ; S3: In Within the previously preset sowing to seedling stage time range, the smoothed time series is used to calculate the second-order difference, and the date corresponding to the maximum value of the second-order difference is determined as the starting point of the key phenological window for cotton. ; S4: Yes and Apply validity constraints when and When the time interval between samples is less than p synthesis periods, pixels exceeding the preset range are reverted to values ​​pre-calculated from cotton samples within the study area. , Mean.

7. The cotton identification method according to claim 1, characterized in that, The method for constructing the enhanced cotton spectral phenological index described in step S5 includes the following steps: S1: Utilizing the timing of dual-source spectral indexes and the starting point of key phenological windows for cotton. and termination point The cumulative spectral phenological index of the high spatial resolution optical remote sensing data path was calculated respectively. Cumulative spectral phenological index of high temporal resolution remote sensing data path Calculated according to the following expression: In the formula, , These represent the reconstructed spectral index time series of the high spatial resolution optical remote sensing data path and the spectral index time series of the high temporal resolution remote sensing data path, respectively, at time t. , This represents the growth rate within the corresponding spectral index time-series phenological window. , Indicates the end point The corresponding spectral index value, Indicates the starting point The corresponding spectral index value, where n represents the number of synthesis periods within the phenological window; S2: with The distribution of is used as the target value range for . Z-score normalization and range alignment are performed to obtain the cumulative spectral phenological index after range alignment. ; S3: Construct pixel-wise adaptive fusion weights based on the temporal gap severity score, local spatial consistency score, and dual-path consistency score within the phenological window. ,include: S3.1: with to Based on the number of effective observations of high spatial resolution optical remote sensing data within the window, a temporal gap severity score is obtained after Sigmoid function transformation. The expression is as follows: In the formula, To normalize the observation density, k is the steepness parameter of the Sigmoid function. This is the inflection point of the Sigmoid function; S3.2: Evaluate each cell The degree of deviation of the value from the mean of its 3×3 neighborhood is used to obtain the local spatial consistency score. The expression is as follows: In the formula, x and y are the abscissa and ordinate, respectively. and Within a 3×3 neighborhood The mean and standard deviation; S3.3: Calculate the value of each cell within the window. , The difference between the mean of the two paths and the mean of its neighborhood is recorded as 1 if they have the same sign, and 0 otherwise. Then, the mean of the two paths is taken in the neighborhood to obtain the bipath consistency score. ; S3.4: The temporal gap severity score, local spatial consistency score, and dual-path consistency score are weighted and summed according to a preset ratio to obtain the pixel-by-pixel adaptive fusion weight. The expression is as follows: In the formula, , , These are the weighting coefficients for the three indicators; S4: Utilizing the adaptive fusion weights right Cumulative spectral phenological index of high temporal resolution remote sensing data path after value range alignment Perform weighted fusion to obtain the weighted fusion result. The expression is as follows: S5: Perform residual bias correction on the weighted fusion results, and calculate the fusion results with... The mean residual Δ between the two results is used to subtract nΔ from the fusion result to eliminate systematic bias, thus obtaining the enhanced cotton spectral phenological index.

8. The cotton identification method according to claim 1, characterized in that, The classification and recognition model mentioned in step S6 is an unsupervised classification model. Obtaining the unsupervised classification result of cotton includes the following steps: The enhanced cotton spectral phenological index is input into the Otsu-Sauvola adaptive threshold segmentation model. The global optimal segmentation threshold is determined in the entire study area by the Otsu algorithm, and the local threshold is calculated pixel by pixel within the preset sliding window based on the Sauvola algorithm. The final segmentation threshold is determined based on the global threshold and the local threshold. The ECSPI is then binarized using the final segmentation threshold to obtain the cotton / non-cotton binarization result. A circular structuring element with a radius of 1 pixel is applied to the binarized result to perform a closing operation, filling the plot holes and smoothing the boundaries, thus obtaining the cotton unsupervised classification result.

9. The cotton identification method according to claim 1, characterized in that, The classification and recognition model mentioned in step S6 is a supervised classification model. Obtaining the supervised classification result of cotton includes the following steps: Using the enhanced cotton spectral phenological index as the core feature, a simplified feature set was constructed by combining the normalized vegetation index time series. A machine learning classifier is trained using the simplified feature set; Input the simplified feature set corresponding to the region to be identified into the trained machine learning classifier, and output the cotton supervised classification result.

10. The cotton identification method according to claim 1, characterized in that, After outputting the spatial distribution results of cotton in the target area in step S6, the method further includes accuracy evaluation and feedback correction of the spatial distribution results of cotton in the target area, specifically including the following steps: A test set based on a spatial block partitioning strategy is used to evaluate the classification accuracy of the cotton spatial distribution results in the target area. The classification accuracy evaluation indicators include overall accuracy OA, Kappa coefficient and F1-score. The cotton identification method was applied to remote sensing data from multiple years to obtain the spatial distribution results of cotton in multiple years, and the multi-year stability of the cotton identification method was evaluated based on the coefficient of variation (CV) of the interannual classification area. When the classification accuracy evaluation index or the multi-year stability evaluation result does not reach the preset threshold, the relevant parameters in steps S2 to S6 are corrected by feedback, and steps S2 to S6 are re-executed based on the corrected parameters. The feedback correction includes at least one of the following: Adjust the local sliding window size, Sauvola correction coefficients, and fusion strategy of global and local thresholds in the Otsu-Sauvola adaptive threshold segmentation model; Adjust the preset physical lower limit or cumulative rainfall quantile level in the annual adaptive rejection threshold, and adjust the rainstorm detection threshold; Adjust the half-window size of the second-order Savitzky-Gorye smoothing filter, the second-order difference search range, the minimum growth cycle constraint, and the starting point. and termination point The effective scope; Adjust the weight coefficients corresponding to the temporal gap severity score, local spatial consistency score, and dual-path consistency score in the adaptive fusion weights, as well as the steepness parameter and inflection point parameter of the Sigmoid function.

Images