Intelligent identification method and system for cultivated land change map based on remote sensing AI intelligent interpretation

Abstract

The application provides a cultivated land change graph spot intelligent recognition method and system based on remote sensing AI intelligent interpretation. First, a multi-temporal remote sensing image data set of a target region is acquired to comprehensively record the cultivated land conditions of the target region at different time periods. Then, the acquired multi-temporal remote sensing image data set is preprocessed. Then, the preprocessed multi-temporal remote sensing image data set is subjected to feature extraction operation to generate a spectral feature sequence and a texture feature sequence of the target region, respectively, to describe the cultivated land features from different angles in detail. Finally, in combination with a time series analysis method, the generated spectral feature sequence and texture feature sequence are subjected to change detection processing to generate a cultivated land change graph spot data set of the target region, thereby effectively improving the accuracy and efficiency of cultivated land change graph spot recognition.

Description

Technical Field

[0001] This invention relates to the field of artificial intelligence technology, and more specifically, to a method and system for intelligent identification of farmland change patches based on remote sensing AI intelligent interpretation. Background Technology

[0002] In today's agricultural resource management and land resource monitoring fields, accurate and efficient monitoring of arable land changes is crucial. Related technologies rely on accurate information about arable land changes to rationally plan grain planting areas and adjust agricultural production strategies, thereby ensuring a stable food supply to meet the growing population's demands. Simultaneously, rational land resource planning also depends on a clear understanding of arable land changes. By understanding the increase, decrease, and distribution of arable land, the proportion of various types of construction land and agricultural land can be scientifically allocated, improving land resource utilization efficiency. Furthermore, arable land changes are closely related to ecological environmental protection. Unreasonable arable land development and utilization can lead to ecological problems such as soil erosion and degradation; therefore, accurate monitoring of arable land changes helps in taking appropriate ecological protection measures.

[0003] Traditional methods for monitoring farmland change mainly fall into two categories: manual field surveys and manual surveys. The former requires significant investment of manpower, resources, and time. Surveyors must personally visit each farmland site to conduct on-site investigations, measurements, and records. Furthermore, the scope of manual field surveys is often limited by manpower and time constraints, making it difficult to comprehensively and rapidly monitor large areas. In addition, the survey results are susceptible to the influence of the surveyors' subjective judgment; different surveyors may reach different conclusions due to differences in experience, knowledge level, and judgment standards, making it difficult to guarantee the accuracy and reliability of the monitoring results.

[0004] Another traditional method is simple visual interpretation of remote sensing images. While this method utilizes remote sensing imagery to obtain farmland information, it demands a high level of expertise from professionals. Interpreters need extensive professional knowledge and practical experience to accurately identify farmland and its changes from the remote sensing images. Moreover, visual interpretation is a tedious and time-consuming process, requiring interpreters to observe remote sensing images frame by frame and area by area, resulting in very low efficiency. Furthermore, visual interpretation is prone to omissions and misjudgments. For subtle changes in farmland, factors such as image resolution and visual fatigue make accurate identification difficult. Additionally, traditional methods often fail to fully extract the temporal series information contained in multi-temporal remote sensing image data. Farmland change is a dynamic process; there are inherent connections and patterns of change between image data from different time points. Traditional methods cannot effectively utilize this information, leading to insufficient dynamic monitoring capabilities for farmland changes and an inability to grasp real-time changes in farmland in a timely and accurate manner. Summary of the Invention

[0005] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for intelligent identification of farmland change patches based on remote sensing AI intelligent interpretation, the method comprising:

[0006] Acquire a multi-temporal remote sensing image data set of the target area, wherein the multi-temporal remote sensing image data set includes remote sensing image data collected at different time nodes;

[0007] The multi-temporal remote sensing image dataset is preprocessed to obtain a preprocessed multi-temporal remote sensing image dataset.

[0008] Feature extraction is performed on the preprocessed multi-temporal remote sensing image dataset to generate spectral feature sequences and texture feature sequences of the target region.

[0009] By combining time series analysis, the spectral feature sequence and the texture feature sequence are subjected to change detection processing to generate a set of cultivated land change patch data for the target area.

[0010] In another aspect, embodiments of the present invention also provide an intelligent identification system for farmland change patches based on remote sensing AI intelligent interpretation, including a processor and a machine-readable storage medium. The machine-readable storage medium is connected to the processor, the machine-readable storage medium is used to store programs, instructions or code, and the processor is used to execute the programs, instructions or code in the machine-readable storage medium to implement the above-mentioned method.

[0011] Based on the above, the embodiments of the present invention significantly improve the efficiency and accuracy of farmland change monitoring. By acquiring a multi-temporal remote sensing image dataset containing data collected at different time points in the target area, covering the farmland information status of the target area at different time stages, and performing preprocessing operations on the acquired multi-temporal remote sensing image dataset, noise, errors, and other interference factors in the data can be effectively removed, improving the quality and usability of the data. The preprocessed data is more in line with the requirements of subsequent feature extraction and analysis, ensuring the accuracy and reliability of the entire identification process. Feature extraction operations are performed on the preprocessed multi-temporal remote sensing image dataset to generate spectral feature sequences and texture feature sequences of the target area, providing a comprehensive and detailed description of the farmland characteristics of the target area from different dimensions. Spectral features can reflect the reflectance characteristics of farmland in different bands, while texture features can reflect the structure and texture information of the farmland surface. The combination of the two can more accurately depict the characteristics of farmland and improve the ability to identify and distinguish farmland. By combining time series analysis with the generated spectral and textural feature sequences for change detection, this method fully utilizes the temporal dimension information in multi-temporal data to keenly capture subtle changes in arable land at different time points. Through analysis of time series data, various types of changes in arable land, such as additions, reductions, and changes in land use, can be accurately identified, generating a dataset of arable land change patches for the target area. This change detection method based on multi-dimensional features and time series analysis effectively overcomes the inefficiencies of traditional methods in monitoring large areas, their inability to identify subtle changes, and their limitations in dynamic monitoring capabilities. Therefore, it can provide agricultural resource management departments and land planning departments with more scientific, accurate, and timely information on arable land changes. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the execution flow of the intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation provided in an embodiment of the present invention.

[0013] Figure 2 This is a schematic diagram of exemplary hardware and software components of the intelligent identification system for farmland change patches based on remote sensing AI intelligent interpretation provided in an embodiment of the present invention. Detailed Implementation

[0014] The present invention will now be described in detail with reference to the accompanying drawings. Figure 1 This is a flowchart illustrating a method for intelligent identification of farmland change patches based on remote sensing AI intelligent interpretation, according to an embodiment of the present invention. The following is a detailed description of this method for intelligent identification of farmland change patches based on remote sensing AI intelligent interpretation.

[0015] Step S110: Obtain a multi-temporal remote sensing image data set of the target area, wherein the multi-temporal remote sensing image data set includes remote sensing image data collected at different time nodes.

[0016] In this embodiment, to achieve intelligent identification of farmland change patches in the target area, a multi-temporal remote sensing image dataset of the target area is first acquired. The construction of this dataset relies on remote sensing sensors mounted on remote sensing satellites, aircraft, etc., to collect data on the target area at different time points. For example, given a target area A, to comprehensively understand the dynamic changes in its farmland, remote sensing image data needs to be acquired at multiple different time points. The selection of these time points should fully consider the characteristic changes of farmland in different growth cycles and seasons, such as data collection during key stages like the sowing, growth, and harvesting periods of crops.

[0017] During data acquisition, remote sensing sensors record the reflectance information of the target area in different bands, forming multi-band remote sensing image data. Each time point in time acquired can be considered an independent dataset, containing the spectral and spatial information of ground features within the target area at that moment. Combining these datasets from different time points forms a multi-temporal remote sensing image dataset. Each subset of data in this multi-temporal remote sensing image dataset corresponds to a specific time point.

[0018] Step S120: Perform preprocessing operations on the multi-temporal remote sensing image data set to obtain a preprocessed multi-temporal remote sensing image data set.

[0019] After acquiring multi-temporal remote sensing image datasets, these raw data may contain issues such as radiometric distortion, geometric aberration, and noise interference, which can affect subsequent feature extraction and change detection results. Therefore, preprocessing is necessary to improve data quality and usability. Preprocessing mainly includes radiometric correction, geometric correction, registration, and denoising steps, as detailed below:

[0020] Step S121: Call the radiometric correction algorithm to perform radiometric correction processing on each remote sensing image data in the multi-temporal remote sensing image data set.

[0021] During the acquisition of remote sensing image data, factors such as atmospheric scattering and absorption, as well as the sensor's own response characteristics, can cause a discrepancy between the image's radiance and the actual reflectance of ground objects. To eliminate this discrepancy, a radiometric correction algorithm is needed to perform radiometric correction processing on each remote sensing image.

[0022] Suppose that for a scene I in a multi-temporal remote sensing image dataset, its original radiometric value is R. The core idea of ​​radiometric correction algorithms is to establish a correction model to convert the original radiometric value R into a corrected radiometric value R', so that the corrected radiometric value can more accurately reflect the true reflectivity of ground objects. Commonly used radiometric correction algorithms include correction methods based on radiative transfer models and correction methods based on ground control points.

[0023] Taking the correction method based on the radiative transfer model as an example, this method needs to consider factors such as the optical properties of the atmosphere, the solar altitude angle, and the sensor's observation angle. First, an atmospheric radiative transfer model is established based on atmospheric optical parameters, such as atmospheric transmittance and atmospheric scattering coefficient. Then, the influence of the atmosphere on the radiance values ​​of the remote sensing image is calculated using this model, and the original radiance values ​​are corrected accordingly. Specifically, for each pixel in image I, its corrected radiance value R' can be calculated through the following steps:

[0024] The first step is to obtain the atmospheric optical parameters of the location of the pixel, including atmospheric transmittance τ and atmospheric scattering coefficient σ. These parameters can be obtained through atmospheric sounding data or meteorological data.

[0025] The second step is to calculate the propagation path length L of solar radiation in the atmosphere based on the solar altitude angle θ and the observation angle φ of the sensor.

[0026] The third step is to use the atmospheric radiative transfer model to calculate the influence of the atmosphere on the radiative value of the pixel and obtain the atmospheric correction term ΔR.

[0027] The fourth step is to subtract the atmospheric correction term ΔR from the original radiation value R to obtain the corrected radiation value R', i.e., R' = R - ΔR.

[0028] By performing this radiometric correction on each image in a multi-temporal remote sensing image dataset, the radiometric data of image data at different time points can be made comparable.

[0029] Step S122: Perform geometric correction processing on the radiometrically corrected remote sensing image data, and eliminate spatial position deviation of the image by matching ground control points. The ground control points are automatically extracted based on the feature corner points of the high-precision topographic map and remote sensing image.

[0030] Radiometrically corrected remote sensing image data may still contain geometric distortions due to factors such as sensor attitude variations, terrain undulations, and the curvature of the Earth. To eliminate these distortions, geometric correction processing is required.

[0031] The key to geometric correction is finding ground control points (GCPs), which are points with clear geographical locations on the image. By matching the control points on the image with the corresponding control points on the high-precision topographic map, a geometric transformation relationship between the image and the topographic map can be established, thereby eliminating spatial positional deviations in the image.

[0032] In this embodiment, the ground control points are automatically extracted based on the feature corner points of high-precision topographic maps and remote sensing images. The specific steps are as follows:

[0033] First, points with distinctive features, such as road intersections and building corners, are extracted from high-precision topographic maps and used as control points on the topographic map. Simultaneously, feature extraction algorithms, such as the Harris corner detection algorithm, are used to extract characteristic corner points from the remote sensing imagery.

[0034] Then, the extracted remote sensing image feature corner points are matched with control points on the topographic map. The matching process can employ feature descriptor-based methods, such as SIFT (Scale Invariant Feature Transform) descriptors or SURF (Accelerated Robust Feature Transform) descriptors. By calculating the similarity between the feature descriptors of the feature corner points and the control points, the best-matching point pairs are found.

[0035] Suppose that the set of feature corner points extracted from remote sensing imagery is C1, and the set of control points extracted from topographic map is C2. For each corner point c1 in C1, calculate its feature descriptor similarity s with each control point c2 in C2. Select the control point with the highest similarity as the matching point to form a matching point pair (c1, c2).

[0036] Finally, based on the matching point pairs, a geometric transformation model, such as an affine transformation model or a polynomial transformation model, is established between the imagery and the topographic map. This geometric transformation model is then used to perform geometric correction on the remote sensing imagery, mapping each pixel in the imagery to the correct geographical location.

[0037] Geometric correction can ensure that remote sensing image data from different time points are consistent in spatial location, providing an accurate spatial reference for subsequent registration and change detection.

[0038] Step S123: Perform a registration operation on the geometrically corrected remote sensing image data to unify the remote sensing image data from different time points into the same geographic coordinate system. The registration operation uses a feature matching algorithm to align the image edges with the contours of ground features.

[0039] Although the remote sensing image data has been corrected in spatial location after geometric correction, there may still be some deviations in the images at different time points. Registration is required to unify them into the same geographic coordinate system.

[0040] The main purpose of registration is to align remote sensing image data from different time points, ensuring that identical ground features in the images coincide spatially. This embodiment uses a feature matching algorithm to achieve image registration, and the specific steps are as follows:

[0041] First, a reference image I_base is selected, usually an image with good quality and high resolution. Then, for other images I_reg that need to be registered, feature extraction algorithms, such as ORB (Oriented Fast and Rotated BRIEF), are used to extract feature points from the images.

[0042] Assume the set of feature points extracted from the reference image I_base is F_base, and the set of feature points extracted from the image to be registered I_reg is F_reg. For each feature point f_reg in F_reg, calculate its similarity s with the feature descriptor of each feature point f_base in F_base. Select the feature point with the highest similarity as the matching point to form a matching point pair (f_reg, f_base).

[0043] To improve matching accuracy, the RANSAC (Random Sample Consensus) algorithm can be used to filter matching point pairs and remove mismatched pairs. The RANSAC algorithm randomly selects a subset of matching point pairs, calculates a geometric transformation model, and then uses this model to verify other matching point pairs, retaining those that conform to the model.

[0044] Based on the selected matching point pairs, calculate the geometric transformation matrix T between image I_reg and reference image I_base. The geometric transformation matrix T can be an affine transformation matrix or a perspective transformation matrix, which describes the spatial transformation relationship from image I_reg to reference image I_base.

[0045] Finally, the image I_reg is transformed using the geometric transformation matrix T and registered to the geographic coordinate system of the reference image I_base. Specifically, for each pixel (x, y) in image I_reg, its corresponding position (x', y') in the reference image coordinate system is calculated using the geometric transformation formula (x', y') = T*(x, y).

[0046] Registration can unify remote sensing image data from different time points in the geographic coordinate system, making it easier to compare and analyze the images later.

[0047] Step S124: Denoise the registered remote sensing image data to obtain a preprocessed multi-temporal remote sensing image data set.

[0048] Registered remote sensing image data may still contain noise, which can affect subsequent feature extraction and change detection results, so it is necessary to perform noise reduction processing.

[0049] Commonly used noise reduction methods include mean filtering, median filtering, and Gaussian filtering. Taking median filtering as an example, median filtering is a non-linear filtering method that removes noise by sorting the pixel values ​​in the neighborhood of each pixel in the image and then taking the median value as the new value of that pixel.

[0050] Suppose that for a registered remote sensing image I, a pixel P needs to be processed using median filtering. First, determine the neighborhood range of this pixel, such as a 3x3 neighborhood or a 5x5 neighborhood. Then, sort all pixel values ​​within the neighborhood and take the median value as the new value for pixel P.

[0051] The specific steps are as follows:

[0052] The first step is to determine the neighborhood range N for each pixel P in image I.

[0053] The second step is to store all pixel values ​​within the neighborhood N in an array.

[0054] The third step is to sort the pixel values ​​in the array.

[0055] The fourth step is to take the middle value of the sorted array as the new value of cell P.

[0056] By performing this median filtering process on each pixel in the registered remote sensing image data, noise in the image can be effectively removed, resulting in a preprocessed multi-temporal remote sensing image data set.

[0057] Step S130: Perform feature extraction on the preprocessed multi-temporal remote sensing image data set to generate the spectral feature sequence and texture feature sequence of the target area.

[0058] The preprocessed multi-temporal remote sensing image dataset contains rich information about the target area. However, in order to more accurately identify farmland change patches, it is necessary to extract useful features from it. This step mainly extracts spectral feature sequences and texture feature sequences, and the specific steps are as follows:

[0059] Step S131: Perform band synthesis processing on each preprocessed remote sensing image data to generate multispectral image data containing visible light and near-infrared bands. The band synthesis processing dynamically selects the band combination according to the vegetation index calculation requirements.

[0060] In the preprocessed multi-temporal remote sensing image dataset, each image typically contains data from multiple bands. To facilitate subsequent vegetation index calculation and feature extraction, these bands need to be synthesized to generate multispectral image data containing visible and near-infrared bands.

[0061] Band synthesis processing requires dynamic selection of band combinations based on the needs of vegetation index calculation. Different vegetation indices have different requirements for band selection. For example, the Normalized Difference Vegetation Index (NDVI) requires near-infrared and red bands, while the Enhanced Vegetation Index (EVI) requires blue, red, and near-infrared bands.

[0062] Suppose that for a preprocessed remote sensing image I, its band set is B = {B1, B2, ..., Bn}. Based on the vegetation index calculation requirements, appropriate bands are selected for synthesis. For example, if NDVI calculation is needed, the near-infrared band Bnir and the red band Bred are selected, and the data from these two bands are combined to generate a new multispectral image data Imulti.

[0063] The specific band synthesis process can be achieved through the following steps:

[0064] The first step is to determine the band set Bselect that needs to be synthesized, and select appropriate bands from the original band set B according to the vegetation index calculation requirements.

[0065] The second step is to combine the data from each band in Bselect according to the same pixel location to form a new multispectral image data Imulti. For example, for each pixel (x, y) in image I, the pixel values ​​of the corresponding band in Bselect are combined together as the multispectral value of that pixel in Imulti.

[0066] By using band synthesis processing, multispectral image data containing visible and near-infrared bands can be generated, providing a foundation for subsequent extraction of spectral and texture features.

[0067] Step S132: Calculate the average reflectance of each pixel in different bands based on the multispectral image data, and generate the spectral feature sequence of the target area. The average reflectance is normalized to eliminate the influence of differences in illumination conditions.

[0068] After obtaining multispectral image data containing visible and near-infrared bands, in order to generate the spectral feature sequence of the target area, it is necessary to calculate the average reflectance of each pixel in different bands and then perform normalization to eliminate the influence of differences in illumination conditions. This process involves a series of detailed steps, as follows:

[0069] Step S1321: Divide the land cover type of the target area according to the preset land cover classification system. The land cover type is automatically labeled by supervised classification algorithm combined with training samples. The land cover type includes at least three types: cultivated land, forest land and water area.

[0070] The pre-defined land cover classification system is constructed based on the spectral characteristics, spatial distribution characteristics, and other relevant attributes of land cover. To achieve accurate classification of land cover types in the target area, a supervised classification algorithm is used, combined with automatic labeling of training samples.

[0071] The first step is the selection of training samples. Training samples need to cover various typical land cover types within the target area. For cultivated land, pixels representing different crops and growth stages should be selected; for forest land, pixels representing different tree species and ages should be considered; for water bodies, pixels representing different water qualities and sizes should be included. The selected training samples must be representative and accurately reflect the characteristics of various land cover types.

[0072] Taking the maximum likelihood classification algorithm as an example, this is a commonly used supervised classification algorithm whose core idea is based on the Bayesian criterion in probability theory. For each pixel in a multispectral image, the algorithm calculates the probability of the pixel belonging to each land cover type based on the training samples, and then classifies it into the category with the highest probability.

[0073] Specifically, for a pixel in a multispectral image, its reflectance values ​​in each band constitute a feature vector. The algorithm calculates the feature mean vector and covariance matrix for each land cover type based on the training samples. For a pixel to be classified, the algorithm calculates its Mahalanobis distance to the feature mean vectors of each land cover type. Mahalanobis distance takes into account the covariance structure of the data and can more accurately measure the similarity between the pixel and the features of various land cover types. By comparing the Mahalanobis distance of the pixel to each land cover type, it is assigned to the land cover type with the closest distance.

[0074] For example, let the feature vector of a pixel be A, the feature mean vector of a certain land cover type be B, and the covariance matrix be C. The process of calculating the Mahalanobis distance D from pixel A to the land cover type is as follows: First, calculate the difference vector E between A and B; then, calculate the product of the transpose of E and the inverse of C; finally, multiply the result by E to obtain the Mahalanobis distance D. This calculation is performed for all land cover types, and the type with the smallest Mahalanobis distance is selected as the classification result for that pixel.

[0075] In this way, all pixels within the target area are classified, thereby completing the automatic labeling of land cover types, and at least three types are identified: cultivated land, forest land, and water area.

[0076] Step S1322: For each land cover type, calculate the average reflectance of the corresponding pixel in the visible light band and near-infrared band, wherein, in the process of calculating the average reflectance, spatially adjacent pixels of the same type are aggregated by a region growing algorithm.

[0077] After classifying land cover types, the average reflectance of the corresponding pixels in the visible and near-infrared bands is calculated for each land cover type. To improve the accuracy and reliability of the calculation results, a region growing algorithm is used to aggregate spatially adjacent pixels of the same type.

[0078] The basic principle of the region growing algorithm is to start with one or more seed pixels and gradually merge pixels with similar characteristics and spatial adjacency into the same region.

[0079] The first step is the selection of seed pixels. For each land cover type, some pixels can be randomly selected as seed pixels, or representative pixels can be selected based on the pixel's characteristic values ​​(such as reflectance).

[0080] Then, starting from the seed pixel, the pixels in its neighborhood are checked. Typically, an 8-neighborhood (i.e., the 8 pixels adjacent to the seed pixel) is used for this check. For a pixel within the neighborhood, if it belongs to the same land cover type as the seed pixel, and the difference in reflectance between the seed pixel and the seed pixel in the visible and near-infrared bands is within a preset threshold range, then that neighboring pixel is merged into the area containing the seed pixel.

[0081] For example, a seed pixel of the cultivated land type has a reflectance of R1 in the visible light band and R2 in the near-infrared band. For a pixel within its 8-neighborhood, its reflectance is r1 in the visible light band and r2 in the near-infrared band. If |r1-R1| is less than the reflectance difference threshold in the visible light band, and |r2-R2| is less than the reflectance difference threshold in the near-infrared band, and this neighboring pixel is also classified as cultivated land, then this neighboring pixel is merged into the region containing the seed pixel.

[0082] Repeat this process, continuously expanding the region, until there are no more neighboring pixels that meet the criteria. This yields the aggregated region of pixels of the same type.

[0083] For each aggregated region, calculate the sum of the reflectance of all pixels in the visible and near-infrared bands, and then divide it by the number of pixels in the region to obtain the average reflectance of the region in the visible and near-infrared bands.

[0084] Step S1323: Construct spectral variation curves for each land cover type based on the mean reflectance values ​​at different time points. The spectral variation curves are fitted with the estimated reflectance values ​​for missing time points using cubic spline interpolation.

[0085] After obtaining the average reflectance values ​​of each land cover type in the visible and near-infrared bands at different time points, spectral variation curves for each land cover type should be constructed. These spectral variation curves can intuitively reflect the changes in reflectance of land cover types over time.

[0086] Since there may be missing data at certain time points during the actual data acquisition process, cubic spline interpolation is used to fit the reflectance estimates of the missing time points in order to make the spectral change curve more continuous and accurate.

[0087] Cubic spline interpolation is a piecewise polynomial interpolation method that approximates the original data by constructing a cubic polynomial over each small interval. Specifically, for a known time point and its corresponding mean reflectance, the time points are divided into several small intervals. In each small interval, a cubic polynomial is constructed such that it equals the known mean reflectance at the endpoints of the intervals and has a continuous second derivative over the entire interval.

[0088] For example, given the mean reflectance values ​​at time points t1, t2, and t3 as R1, R2, and R3 respectively, a cubic polynomial P1(t) = a1*t^3 + b1*t^2 + c1*t + d1 is constructed over the small interval from t1 to t2, where a1, b1, c1, and d1 are the coefficients to be determined. These coefficients are determined by requiring P1(t1) = R1, P1(t2) = R2, and that the first and second derivatives of P1(t) at t1 and t2 are equal to the first and second derivatives of the polynomials between adjacent cells at their respective endpoints.

[0089] By performing this process on all the smaller intervals, we obtain the cubic spline interpolation function for the entire time interval. For the reflectance estimate of a missing time node, the result is obtained by substituting that time node into the corresponding cubic spline interpolation function.

[0090] In this way, spectral variation curves for each land cover type are constructed, making the reflectance data more continuous and complete in the time dimension.

[0091] Step S1324: Identify the spectral abnormal fluctuation period of cultivated land type pixels based on the spectral change curve, use the mutation point detection algorithm to locate the time node where the reflectance changes by a set time interval greater than a set amplitude, generate the temporal change marker in the spectral feature sequence, perform time intersection operation on the spectral abnormal fluctuation period and the cultivated land change patch occurrence period, and retain the patches whose overlapping period exceeds a preset threshold.

[0092] Based on the constructed spectral variation curves for each land cover type, the focus is on identifying periods of abnormal spectral fluctuations in cultivated land type pixels. A mutation point detection algorithm is used to locate time points where the reflectance changes by a greater than a set range within a set time interval.

[0093] The basic idea of ​​the mutation point detection algorithm is to find the time point where the reflectance changes significantly by comparing the reflectance changes between adjacent time points or within a set time interval.

[0094] For example, given a time interval T, for the spectral variation curve of cultivated land type pixels, the reflectance difference between each time node and the previous T time nodes is calculated. If the reflectance difference at a certain time node is greater than a set amplitude threshold, that time node is considered a point of abrupt change.

[0095] After identifying the abrupt change points, the time intervals between adjacent abrupt change points are defined as spectral anomalous fluctuation periods. These spectral anomalous fluctuation periods are marked as time-series change markers and added to the spectral feature sequence.

[0096] Then, a temporal intersection operation is performed between the periods of spectral anomaly fluctuations and the periods of farmland change patch occurrences. The process involves identifying the overlapping time portion for each period of spectral anomaly fluctuation and each period of farmland change patch occurrences. The proportion of this overlapping time within each period of spectral anomaly fluctuation and farmland change patch occurrence is then calculated.

[0097] Only when the proportion of overlapping time periods in the periods of spectral anomalies or the periods in which farmland change patches appear exceeds a preset threshold are the corresponding farmland change patches retained. This allows for the screening of real farmland change patches related to spectral anomalies, improving the accuracy of farmland change patch identification.

[0098] Step S13241: Extract historical reflectance data of cultivated land type pixels during the growing season and non-growing season, perform time dimension normalization processing on the historical reflectance data, and generate normalized reflectance baseline intervals for each time node.

[0099] To more accurately identify periods of spectral anomalies in cultivated land type pixels, it is necessary to extract historical reflectance data of cultivated land type pixels during both growing and non-growing seasons. This historical reflectance data records the reflectance of cultivated land at different times and can reflect the reflectance characteristics of cultivated land under normal growing and non-growing conditions.

[0100] When extracting historical reflectance data, it is essential to ensure the accuracy and completeness of the data. Reflectance data for cultivated land type pixels can be filtered from historical multispectral imagery and then categorized according to growing season and non-growing season.

[0101] The extracted historical reflectance data is normalized over time to eliminate the influence of factors such as lighting conditions and atmospheric conditions at different time points, so as to make the data comparable.

[0102] Time-dimensional normalization can be achieved using the min-max normalization method. For reflectance data at each time point, the minimum and maximum values ​​in the dataset are identified. For a reflectance value R, the normalized reflectance value R' is calculated as follows: first, calculate the difference between R and the minimum value, then divide by the difference between the maximum and minimum values.

[0103] By normalizing the reflectance data across all time points, normalized reflectance data is obtained. Based on this normalized reflectance data, the mean and standard deviation of reflectance are calculated for each time point. Using the mean reflectance as the center, a baseline interval for the normalized reflectance at each time point is determined according to a predetermined multiple of the standard deviation. For example, the baseline interval could be the range of the mean plus or minus twice the standard deviation.

[0104] Step S13242: Obtain the multispectral image reflectance data at the current time node, and generate the normalized reflectance vector at the current time node using the same normalization processing method as the historical reflectance data.

[0105] Obtain multispectral image reflectance data for the current time point, which reflects the reflectance of farmland at the current moment.

[0106] The same normalization process used for historical reflectance data was applied to the multispectral image reflectance data at the current time point. Specifically, the minimum and maximum values ​​in the reflectance dataset at the current time point were identified, and the normalized value for each reflectance value was calculated using the minimum-maximum normalization method.

[0107] The normalized reflectance values ​​for each band at the current time point are combined to form the normalized reflectance vector for the current time point. This vector can comprehensively reflect the spectral characteristics of the cultivated land at the current time point.

[0108] Step S13243: Calculate the normalized Euclidean distance between the normalized reflectance vector and the baseline interval of the corresponding time node, and generate the spectral deviation index for each time node.

[0109] Calculate the normalized Euclidean distance between the normalized reflectance vector at the current time node and the baseline interval at the corresponding time node to generate the spectral deviation index for each time node.

[0110] Standardized Euclidean distance is a distance metric that takes into account the characteristics of data distribution. For the normalized reflectance vector V at the current time point and the baseline interval at the corresponding time point, the baseline interval can be represented by the baseline mean vector M and the standard deviation vector S.

[0111] The process of calculating the standardized Euclidean distance is as follows: First, calculate the difference vector D between the normalized reflectance vector V and the baseline mean vector M. Then, divide each element of the difference vector D by the corresponding element of the standard deviation vector S to obtain the standardized difference vector D'. Finally, calculate the Euclidean distance of the standardized difference vector D', which is the square root of the sum of the squares of the elements of the standardized difference vector D'.

[0112] This standardized Euclidean distance serves as an indicator of spectral deviation at that time point. A larger spectral deviation indicator indicates a greater difference between the reflectance at the current time point and the baseline range, potentially suggesting spectral anomalies.

[0113] Step S13244: Perform weighted averaging of the spectral deviation index at consecutive time nodes based on the sliding window mechanism to generate a dynamic deviation threshold curve.

[0114] The spectral deviation index at consecutive time points is weighted and averaged using a sliding window mechanism to generate a dynamic deviation threshold curve.

[0115] The sliding window mechanism involves setting a fixed-size window on the time series, which moves gradually across the time series. For the spectral deviation index within each window, a weighted average method is used for processing.

[0116] The weights in a weighted average can be determined based on the distance between time nodes within the window and the center time node of the window. Generally, time nodes closer to the center time node have a larger weight. For example, a Gaussian weighting function can be used to determine the weights; the closer a time node is to the center time node, the larger its weight value, and the weight values ​​follow a Gaussian distribution.

[0117] For each window position, a weighted average of the spectral deviation index within the window is calculated. As the window moves over time, the weighted average is continuously updated to obtain a series of weighted average spectral deviation index values.

[0118] Connecting these weighted average spectral deviation index values ​​forms a dynamic deviation threshold curve. This curve dynamically adjusts the threshold for judging spectral anomalies based on changes over time, improving the accuracy of identifying spectral anomaly fluctuations.

[0119] Step S13245: Identify the time periods in which the spectral deviation index continuously exceeds the dynamic deviation threshold curve, and extract its start time node and end time node as candidate abnormal fluctuation periods.

[0120] After obtaining the dynamic deviation threshold curve, identify the period during which the spectral deviation index continuously exceeds the dynamic deviation threshold curve.

[0121] Starting from the beginning of the time series, the spectral deviation index of each time point is examined one by one. If the spectral deviation index of a certain time point exceeds the threshold corresponding to the dynamic deviation threshold curve, and the spectral deviation index of several subsequent consecutive time points continues to exceed the threshold, then this period of continuous exceedance of the threshold is recorded.

[0122] The start and end times of this period are extracted and used as candidate periods of abnormal fluctuation. These candidate periods of abnormal fluctuation indicate that the reflectance of cultivated land has been continuously deviating from the normal range over a period of time, which may indicate an anomaly.

[0123] Step S13246: Match the candidate abnormal fluctuation period with the time series of meteorological disaster events by timestamp, and calculate the percentage of overlapping time.

[0124] The purpose of matching the candidate abnormal fluctuation periods with the time series of meteorological disaster events with timestamps is to eliminate spectral abnormal fluctuations caused by natural factors such as meteorological disasters and to identify spectral abnormal fluctuations that are truly caused by changes in cultivated land.

[0125] The time series of meteorological disaster event records contains information on the timing of meteorological disasters occurring within the target area. For each candidate abnormal fluctuation period, it is checked whether it overlaps with the time intervals of each disaster event in the time series of meteorological disaster event records.

[0126] The process of calculating the overlap time ratio is as follows: For a candidate anomalous fluctuation period and a meteorological disaster event time interval, find the length of their overlap. Then, divide the overlap time length by the length of the candidate anomalous fluctuation period and the length of the meteorological disaster event time interval, respectively, to obtain the two overlap time ratios.

[0127] Step S13247: Based on the continuous growth rate of the overlap time ratio and the spectral deviation index, establish a pseudo-anomaly fluctuation discrimination function, eliminate candidate periods caused by meteorological interference, and output the spectral anomaly fluctuation period corresponding to the actual changes in cultivated land.

[0128] Based on the overlap time ratio and the continuous growth rate of the spectral deviation index, a pseudo-anomaly fluctuation discrimination function is established. The purpose of the pseudo-anomaly fluctuation discrimination function is to determine whether the candidate anomaly fluctuation period is caused by meteorological disturbances or by actual changes in cultivated land.

[0129] The sustained growth rate of the spectral deviation index can be obtained by calculating the difference in the spectral deviation index between adjacent time points and dividing by the time interval.

[0130] The pseudo-anomaly fluctuation discrimination function can be a function that comprehensively considers the overlap time ratio and the continuous growth rate of the spectral deviation index. For example, a threshold combination can be set. When the overlap time ratio exceeds a certain threshold and the continuous growth rate of the spectral deviation index is lower than another threshold, the candidate anomaly fluctuation period is considered to be caused by meteorological interference and is removed.

[0131] By judging and screening all candidate abnormal fluctuation periods, eliminating candidate periods caused by meteorological interference, the spectral abnormal fluctuation periods corresponding to the actual changes in cultivated land are finally output.

[0132] Therefore, the average reflectance of each pixel in different bands can be accurately calculated based on multispectral image data, generating a spectral feature sequence of the target area, effectively eliminating the influence of differences in lighting conditions, and accurately identifying the periods of abnormal spectral fluctuations in cultivated land type pixels, providing reliable spectral feature information for subsequent cultivated land change patch identification.

[0133] Step S133: Extract the gray-level co-occurrence matrix of the multispectral image data, calculate the contrast and homogeneity index of each pixel based on the gray-level co-occurrence matrix, and generate the texture feature sequence of the target area. The contrast index is used to characterize the edge sharpness of ground features, and the homogeneity index is used to quantify the uniformity of gray-level distribution in the local area.

[0134] Besides spectral features, texture features are also important information for identifying changes in cultivated land. This step extracts the gray-level co-occurrence matrix from the multispectral image data and calculates the contrast and homogeneity indices of each pixel based on this matrix to generate a texture feature sequence for the target area.

[0135] First, the multispectral image data Imulti is converted into a grayscale image Igrey. A weighted averaging method can be used to convert multiple bands of the multispectral image into grayscale values. For example, for each pixel (x, y) in the multispectral image, its grayscale value G(x, y) can be calculated using the following formula:

[0136] G(x, y)=w1*R1(x, y)+w2*R2(x, y)+...+wm*Rm(x, y)

[0137] Where w1, w2, ..., wm are the weight coefficients of each band, and satisfy w1 + w2 + ... + wm = 1.

[0138] Then, for each pixel in the grayscale image Igrey, its gray-level co-occurrence matrix is ​​calculated. The gray-level co-occurrence matrix is ​​a matrix that describes the spatial distribution of gray levels in an image, reflecting the adjacency relationship between different gray levels. Assuming the gray level range of the grayscale image is [0, L-1], for a given distance d and angle θ, the gray-level co-occurrence matrix GLCM(i, j; d, θ) represents the co-occurrence frequency of the pixel with gray value i and the pixel with gray value j at distance d and angle θ.

[0139] Next, the contrast and homogeneity indices of each pixel are calculated based on the gray-level co-occurrence matrix. The contrast index C is used to characterize the sharpness of ground feature edges, and its calculation formula is as follows:

[0140] C=∑(i,j)[(ij)^2*GLCM(i,j;d,θ)]

[0141] The homogeneity index H is used to quantify the uniformity of gray-level distribution in a local area, and its calculation formula is as follows:

[0142] H=∑(i,j)[GLCM(i,j;d,θ) / (1+(ij)^2)]

[0143] Finally, the contrast and homogeneity indices of all pixels in the target area are arranged in a set order to form a texture feature sequence T=[C1, H1, C2, H2, ..., Ck, Hk], where Ci and Hi are the contrast and homogeneity indices of the i-th pixel, respectively.

[0144] By extracting the gray-level co-occurrence matrix of multispectral image data and calculating contrast and homogeneity indices, a texture feature sequence of the target area can be generated, providing important texture information for the identification of farmland change patches.

[0145] Step S134: Spatial alignment operation is performed on the spectral feature sequence and the texture feature sequence at the same time point to generate the spectral feature sequence and texture feature sequence of the target region.

[0146] After generating spectral feature sequences and texture feature sequences respectively, in order to ensure their spatial consistency, it is necessary to perform spatial alignment operations on the spectral feature sequences and texture feature sequences at the same time point.

[0147] Since both the spectral feature sequence and the texture feature sequence are generated from the same multispectral image data, their pixel positions are one-to-one. Therefore, spatial alignment can be achieved simply by combining the spectral and texture features of the same pixel.

[0148] Suppose the spectral feature sequence at the same time point is S=[S1, S2, ..., Sk], and the texture feature sequence is T=[T1, T2, ..., Tk], where Si and Ti are the spectral and texture features of the i-th pixel, respectively. After spatial alignment, a new feature sequence F=[F1, F2, ..., Fk] is generated, where Fi=[Si, Ti], meaning that the feature of each pixel is composed of its spectral and texture features.

[0149] Spatial alignment operations can organically combine spectral and texture features to generate spectral and texture feature sequences of the target region, providing more comprehensive feature information for subsequent change detection.

[0150] Step S140: Combine time series analysis to perform change detection processing on the spectral feature sequence and the texture feature sequence to generate a set of cultivated land change patch data for the target area.

[0151] After obtaining the spectral and texture feature sequences of the target area, in order to identify the changed patches of cultivated land, it is necessary to perform change detection processing on these two feature sequences by combining time series analysis. This process includes the following sub-steps:

[0152] Step S141: Construct a time series analysis model. Input the spectral feature sequence and the texture feature sequence into the time series analysis model in chronological order. The time series analysis model uses a sliding window mechanism to dynamically capture the feature evolution trend.

[0153] First, a time series analysis model needs to be constructed. This model analyzes spectral and texture feature sequences to capture their trends over time. A sliding window mechanism is used here, which allows for dynamic observation of feature changes over time.

[0154] Suppose that the spectral feature sequence is represented by S and the texture feature sequence is represented by T, both of which are feature sets arranged in chronological order. The size of the sliding window is represented by w, and the step size of the window moving across the time series is represented by s.

[0155] When building time series analysis models, several models suitable for processing time series data can be used, such as Long Short-Term Memory (LSTM) networks or Gated Recurrent Units (GRUs). Taking LSTM as an example, it has an input layer, hidden layers, and an output layer. The input layer receives data from spectral feature sequences and texture feature sequences, the LSTM units in the hidden layers can memorize long-term dependency information in the time series, and the output layer outputs the prediction results of feature change trends.

[0156] The spectral feature sequence S and the texture feature sequence T are input into the time series analysis model in chronological order. For each time point t, the spectral and texture features of that time point and the w-1 time points before it are combined into an input vector, which is then input into the LSTM model. For example, at time point t, the input vector can be represented as [S(t-w+1), T(t-w+1), S(t-w+2), T(t-w+2), ..., S(t), T(t)].

[0157] By using a sliding window mechanism, the model can dynamically capture the evolution trend of features. As the window moves over time, the model continuously updates its analysis of feature changes.

[0158] Step S142: Calculate the spectral feature difference and texture feature variability between adjacent time nodes in the time series analysis model. The spectral feature difference is measured by Euclidean distance to measure the change in band reflectance, and the texture feature variability is quantified based on the time derivative of the contrast index to quantify the rate of change of texture structure.

[0159] In time series analysis models, it is necessary to calculate the spectral feature differences and texture feature variability between adjacent time points in order to determine whether farmland has changed.

[0160] For spectral feature difference, Euclidean distance is used to measure the magnitude of change in band reflectance. Assuming the spectral feature vectors at time points t and t+1 are S(t) and S(t+1) respectively, the spectral feature difference Ds can be obtained by calculating the Euclidean distance between these two vectors. Specifically, the calculation process involves first calculating the sum of squared differences between corresponding elements of the two vectors, and then taking the square root of that sum. For example, if S(t) = [s1(t), s2(t), ..., sn(t)] and S(t+1) = [s1(t+1), s2(t+1), ..., sn(t+1)], then Ds is the sum of the squares of (s1(t+1) - s1(t)) and (s2(t+1) - s2(t)) up to (sn(t+1) - sn(t)), and then taking the square root of that sum.

[0161] For texture feature variability, the rate of change of texture structure is quantified based on the time derivative of the contrast index. Assuming the contrast indices of the texture features at time points t and t+1 are C(t) and C(t+1) respectively, the texture feature variability Dt can be obtained by calculating the ratio of (C(t+1) - C(t)) to the time interval. The time interval is typically one time unit (in the context of this embodiment), so the texture feature variability Dt is approximately equal to C(t+1) - C(t).

[0162] By calculating the spectral feature difference and texture feature variability, the changes in features between adjacent time points can be quantified, providing a basis for subsequent change detection.

[0163] Step S143: Standardize the spectral feature difference degree and texture feature variation degree respectively to generate a normalized difference index and a normalized variation index. When both the normalized difference index and the normalized variation index exceed the preset change detection threshold, mark the corresponding pixel as a change candidate region.

[0164] Since the ranges of spectral feature difference and texture feature variability may differ, they need to be standardized separately for easier comparison and unified processing.

[0165] Standardization can be achieved using the Z-score standardization method. For the spectral feature variability Ds, first calculate its mean μs and standard deviation σs. Then, for each spectral feature variability value Ds(i), its normalized variability index NDs(i) can be obtained by dividing (Ds(i) - μs) by σs. For the texture feature variability Dt, similarly calculate its mean μt and standard deviation σt. The normalized variability index NDt(i) for each texture feature variability value Dt(i) can be obtained by dividing (Dt(i) - μt) by σt.

[0166] A change detection threshold Th is set. When both the normalized difference index NDs(i) and the normalized variation index NDt(i) exceed the change detection threshold Th, it indicates that the pixel has undergone significant changes in both spectral and texture features. The pixel is then marked as a candidate region for change, which can initially screen out pixels that may have undergone changes in cultivated land.

[0167] Step S144: Perform spatial clustering on all candidate change regions, and use the density peak clustering algorithm to merge adjacent pixels to form initial change patches. The density peak clustering algorithm determines the initial cluster center by calculating the number of adjacent pixels of each pixel within a preset neighborhood range.

[0168] After obtaining the candidate regions of change, in order to merge adjacent changed pixels into meaningful changed patches, these candidate regions need to be spatially clustered. Here, the density peak clustering algorithm is used.

[0169] The algorithm first calculates the number of neighboring pixels of each pixel within a preset neighborhood. Assuming the preset neighborhood is represented by a radius r, for each candidate pixel p, a circular neighborhood with radius r is drawn centered on p, and the number of candidate pixels within this neighborhood is counted. This number of candidate pixels is the local density ρ(p) of that pixel.

[0170] Then, the minimum distance δ(p) between each cell and the cell with a higher local density is calculated. By comparing the local density ρ(p) and the minimum distance δ(p) of all cells, the cells with both larger local density and larger minimum distance are identified, and these cells are the initial cluster centers.

[0171] For each initial cluster center, the cells in its neighborhood are merged into the cluster to which that cluster center belongs, forming the initial changed patch. Specifically, the merging process involves determining whether a cell is in the neighborhood of an initial cluster center; if so, the cell is marked as part of that cluster.

[0172] Step S145: Perform morphological filtering on the initial changed patches, sequentially perform erosion operation to eliminate isolated noise points and dilation operation to smooth patch boundaries, and generate the farmland changed patch data set. The size of the structuring element of the morphological filter is adaptively selected according to the average area of ​​the patches.

[0173] The initial variation patches obtained after spatial clustering may contain some isolated noise points, and the boundaries of the patches may not be smooth enough. Therefore, morphological filtering is required for the initial variation patches.

[0174] Morphological filtering includes erosion and dilation operations. The purpose of erosion is to eliminate isolated noise points. A suitable structuring element is selected, the size of which is adaptively chosen based on the average area of ​​the patch. Assuming the structuring element is represented by SE, for each pixel in the initial changed patch, if the pixel's neighborhood perfectly matches the structuring element SE, the pixel is retained; otherwise, it is deleted. Through this operation, some isolated small noise points can be removed.

[0175] The purpose of the dilation operation is to smooth the boundaries of the patch. After the erosion operation, a dilation operation is performed on the patch. Using the same structuring element SE, for each cell in the patch, all neighboring cells that intersect with the structuring element SE are marked as part of the patch. Through the dilation operation, the boundaries of the patch can be made smoother.

[0176] After the erosion and expansion operations, the resulting patches are the final set of farmland change patch data.

[0177] Based on the above embodiments, the generated farmland change map data set can be further processed or applied, such as for result display and data storage. In this embodiment, it is assumed that the generated farmland change map data set is stored in a database for subsequent querying and analysis. The relevant information of each map patch in the farmland change map data set, such as location, area, and change time, is stored in the corresponding table of the database according to a set data format. This facilitates further statistical analysis of farmland changes.

[0178] Furthermore, the method may also include:

[0179] Step S210: Extract the spatial distribution features and temporal frequency features from the farmland change patch data set. The spatial distribution features are quantified by kernel density estimation, and the temporal frequency features are calculated based on the frequency of patch occurrence and duration.

[0180] To further optimize the change detection results, it is necessary to extract the spatial distribution features and temporal frequency features from the farmland change patch dataset.

[0181] For spatial distribution characteristics, kernel density estimation is used to quantify the clustering degree of patches. The basic idea of ​​kernel density estimation is to calculate the density of patches within a defined range around the center location of each patch, using a kernel function as the weight. Assuming the kernel function is represented by K and the bandwidth by h, the kernel density estimate f(x) for a location x can be obtained by weighted summation over the center locations xi of all patches. The weight is the kernel function K((x-xi) / h) divided by the bandwidth h raised to the power of n (n is the spatial dimension, usually 2). By calculating the kernel density estimate for each location, a spatial distribution density map of the patches can be obtained, thereby quantifying the clustering degree of the patches.

[0182] For the time-varying frequency characteristic, it is calculated based on the occurrence frequency and duration of the patch. The number of times each patch appears in the entire time series is counted, denoted as the occurrence frequency N, and the duration T of each patch is also recorded. A time-varying frequency index F can be obtained by dividing the occurrence frequency N by the duration T, which reflects the frequency of patch changes.

[0183] Step S220: Adjust the geographic weight parameter of the change detection threshold according to the spatial distribution characteristics.

[0184] Based on the extracted spatial distribution features, the geographic weight parameters of the change detection threshold can be adjusted to improve the accuracy of change detection. The specific steps are as follows:

[0185] Step S221: Divide the target area into multiple geographic grid units.

[0186] The target area is divided into multiple geographic grid units of the same size. Each grid unit can be regarded as an independent geographic region. The purpose of this division is to take into account the characteristic differences of different geographic regions in more detail.

[0187] Step S222: Calculate the degree of clustering of map features within each geographic grid cell. The degree of clustering is calculated by combining the number of map features per unit area with the average area of ​​map features.

[0188] For each geographic grid cell, the number of polygons N and the total area A are counted. Then, the number of polygons per unit area N / A is calculated, along with the average area A / N. By combining the number of polygons per unit area and the average area, a polygon clustering index C can be obtained. For example, the number of polygons per unit area and the average area can be weighted and summed to obtain the polygon clustering index C = α*(N / A) + β*(A / N), where α and β are weighting coefficients, and α + β = 1.

[0189] Step S223: Establish a regression model for the elevation, slope data and change patch density of geographic grid units. The regression model uses the ridge regression algorithm to handle the multicollinearity problem.

[0190] Geographic data such as elevation H and slope S are collected for each geographic grid cell. Simultaneously, the density of changeable patches D within each grid cell is calculated (the density of changeable patches can be obtained by dividing the number of patches by the area of ​​the grid cell). A regression model is established to predict the density of changeable patches D using elevation H and slope S. Due to the potential for multicollinearity in elevation and slope data, ridge regression is employed to address this issue. The ridge regression algorithm, based on ordinary least squares, incorporates a regularization term to reduce the variance of parameter estimates and improve model stability.

[0191] Step S224: Predict the theoretical change density under different geographical conditions based on the regression relationship model. The theoretical change density reflects the natural constraint effect of topographic factors on farmland change.

[0192] Using the established regression model, and inputting elevation and slope data for different geographic grid units, the theoretical change density D' of each geographic grid unit is predicted. This theoretical change density reflects the natural constraint effect of topographic factors on farmland change, meaning that the probability of farmland change varies under different topographic conditions.

[0193] Step S225: Calculate the residual value between the actual change density and the theoretical change density, dynamically adjust the geographic weight parameter of the corresponding geographic grid cell based on the residual value, reduce the geographic weight parameter of the change detection threshold of the positive residual area according to the first dynamic mapping coefficient, and increase the geographic weight parameter of the change detection threshold of the negative residual area according to the second dynamic mapping coefficient.

[0194] For each geographic grid cell, calculate the residual value R = D - D' between the actual change density D and the theoretical change density D'. If the residual value R is positive, it indicates that the actual change density of the geographic grid cell is higher than the theoretical change density, and there may be some farmland changes caused by non-topographic factors. The geographic weight parameter of the change detection threshold for this area is reduced according to the first dynamic mapping coefficient k1, making it easier to detect these changes. If the residual value R is negative, it indicates that the actual change density of the geographic grid cell is lower than the theoretical change density. The geographic weight parameter of the change detection threshold for this area is increased according to the second dynamic mapping coefficient k2 to reduce the possibility of false detections.

[0195] Step S230: Optimize the time window length parameter of the time series analysis model according to the time change frequency characteristics. Shorten the time window length in the changing region to capture changing events, and extend the time window length in the stable region. When the time window length exceeds the preset maximum window, reset the starting point of the time window to the current time node minus a fixed period.

[0196] Based on the extracted time-varying frequency characteristics, the time window length parameter of the time series analysis model can be optimized.

[0197] For areas of high frequency of change, such as farmland, the time window length should be shortened. This is because farmland changes frequently in these areas, and a shorter time window allows for more timely detection of these events. For example, the time window length could be shortened from the original w to w1 (w1... <w)。

[0198] For stable regions, i.e., regions with low frequency of time variation, the time window length can be extended. A longer time window can smooth out some small fluctuations and more accurately reflect the characteristics of the stable region. For example, the time window length can be extended from the original w to w2 (w2>w).

[0199] When the time window length exceeds the preset maximum window length Wmax, in order to avoid information redundancy caused by an excessively long time window, the starting point of the time window is reset to the current time node minus a fixed period T. This ensures that the time window is within a reasonable range, while also allowing continuous monitoring of changes in cultivated land.

[0200] For data collection involved in the above process, if privacy-sensitive data is involved, differential privacy technology is used for privacy protection and leakage prevention. The basic idea of ​​differential privacy technology is to add a certain amount of noise to the data so that the difference between the query results before and after adding noise is not obvious, thereby protecting the privacy of the data. Specifically, for the collected remote sensing image data and related feature data, Laplace noise is added to the sensitive information in the data before storage and processing. The parameters of the Laplace noise are determined according to the sensitivity of the data and the privacy budget. In this way, the leakage of privacy-sensitive data can be effectively prevented while ensuring data availability.

[0201] Figure 2 The illustration shows exemplary hardware and software components of an intelligent farmland change map patch identification system 100 based on remote sensing AI intelligent interpretation, which can implement the ideas of this application, according to some embodiments of this application. For example, a processor 120 can be used in the intelligent farmland change map patch identification system 100 based on remote sensing AI intelligent interpretation and to perform the functions in this application.

[0202] The intelligent identification system 100 for farmland change patches based on remote sensing AI intelligent interpretation can be a general-purpose server or a special-purpose server. Both can be used to implement the intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation of this application. Although only one server is shown in this application, for convenience, the functions described in this application can be implemented in a distributed manner on multiple similar platforms to balance the processing load.

[0203] For example, the intelligent identification system 100 for farmland change maps based on remote sensing AI intelligent interpretation may include a network port 110 connected to a network, one or more processors 120 for executing program instructions, a communication bus 130, and various forms of storage media 140, such as a disk, ROM, or RAM, or any combination thereof. Exemplarily, the intelligent identification system 100 for farmland change maps based on remote sensing AI intelligent interpretation may also include program instructions stored in ROM, RAM, or other types of non-transitory storage media, or any combination thereof. The methods of this application can be implemented according to these program instructions. The intelligent identification system 100 for farmland change maps based on remote sensing AI intelligent interpretation also includes an I / O interface 150 between the computer and other input / output devices.

[0204] For ease of explanation, only one processor is described in the intelligent farmland change map identification system 100 based on remote sensing AI intelligent interpretation. However, it should be noted that the intelligent farmland change map identification system 100 based on remote sensing AI intelligent interpretation in this application may also include multiple processors. Therefore, the steps performed by one processor described in this application may also be performed jointly by multiple processors or individually. For example, if the processor of the intelligent farmland change map identification system 100 based on remote sensing AI intelligent interpretation performs steps A and B, it should be understood that steps A and B may also be performed jointly by two different processors or individually by one processor. For example, the first processor performs step A, the second processor performs step B, or the first processor and the second processor jointly perform steps A and B.

[0205] Furthermore, this embodiment of the invention also provides a readable storage medium, wherein computer-executable instructions are preset in the readable storage medium, and when the processor executes the computer-executable instructions, the above-mentioned intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation is realized.

[0206] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for intelligent identification of farmland change patches based on remote sensing AI intelligent interpretation, characterized in that, The method includes: Acquire a multi-temporal remote sensing image data set of the target area, wherein the multi-temporal remote sensing image data set includes remote sensing image data collected at different time nodes; The multi-temporal remote sensing image dataset is preprocessed to obtain a preprocessed multi-temporal remote sensing image dataset. Feature extraction is performed on the preprocessed multi-temporal remote sensing image dataset to generate spectral feature sequences and texture feature sequences of the target region. By combining time series analysis, the spectral feature sequence and the texture feature sequence are processed to detect changes, thereby generating a set of farmland change patch data for the target area. The step of performing feature extraction on the preprocessed multi-temporal remote sensing image dataset to generate spectral feature sequences and texture feature sequences for the target region includes: Each pre-processed remote sensing image data is subjected to band synthesis processing to generate multispectral image data containing visible light and near-infrared bands. The band synthesis processing dynamically selects the band combination according to the vegetation index calculation requirements. The average reflectance of each pixel in different bands is calculated based on the multispectral image data to generate the spectral feature sequence of the target area. The average reflectance is normalized to eliminate the influence of differences in lighting conditions. Extract the gray-level co-occurrence matrix of the multispectral image data, calculate the contrast and homogeneity index of each pixel based on the gray-level co-occurrence matrix, and generate the texture feature sequence of the target area. The contrast index is used to characterize the sharpness of the ground object edge, and the homogeneity index is used to quantify the uniformity of gray-level distribution in the local area. Spatially align the spectral feature sequence and the texture feature sequence at the same time point to generate the spectral feature sequence and texture feature sequence of the target region.

2. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 1, characterized in that, The preprocessing operation on the multi-temporal remote sensing image dataset to obtain a preprocessed multi-temporal remote sensing image dataset includes: The radiometric correction algorithm is invoked to perform radiometric correction processing on each scene of remote sensing image data in the multi-temporal remote sensing image dataset. The radiometrically corrected remote sensing image data is subjected to geometric correction processing, and the spatial position deviation of the image is eliminated by ground control point matching. The ground control points are automatically extracted based on the feature corner points of high-precision topographic maps and remote sensing images. The geometrically corrected remote sensing image data is registered to unify remote sensing image data from different time points into the same geographic coordinate system. The registration operation uses a feature matching algorithm to align image edges with ground feature outlines. The registered remote sensing image data is then denoised to obtain a preprocessed multi-temporal remote sensing image data set.

3. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 1, characterized in that, The step of calculating the average reflectance of each pixel in different bands based on the multispectral image data to generate the spectral feature sequence of the target area includes: The land cover types of the target area are classified according to a preset land cover classification system. The land cover types are automatically labeled by a supervised classification algorithm combined with training samples. The land cover types include at least three types: cultivated land, forest land, and water area. For each land cover type, the average reflectance of the corresponding pixel in the visible light band and near-infrared band is calculated. In the process of calculating the average reflectance, spatially adjacent pixels of the same type are aggregated by a region growing algorithm. Spectral variation curves for each land cover type are constructed based on the mean reflectance values ​​at different time points. The spectral variation curves are fitted with the reflectance estimates for missing time points using cubic spline interpolation. Based on the spectral change curve, the period of spectral abnormal fluctuation of cultivated land type pixels is identified. The mutation point detection algorithm is used to locate the time node where the change in reflectance is greater than the set range in a set time interval. The temporal change marker in the spectral feature sequence is generated. The time intersection operation is performed between the period of spectral abnormal fluctuation and the period of occurrence of cultivated land change patches, and the patches with overlapping periods exceeding the preset threshold are retained.

4. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 3, characterized in that, The period of spectral anomaly fluctuation in the arable land type pixels identified based on the spectral change curve includes: Historical reflectance data of cultivated land type pixels during the growing season and non-growing season are extracted, and the historical reflectance data is normalized in the time dimension to generate normalized reflectance baseline intervals for each time node. Obtain the multispectral image reflectance data at the current time point, and use the same normalization processing method as the historical reflectance data to generate the normalized reflectance vector at the current time point; Calculate the normalized Euclidean distance between the normalized reflectance vector and the baseline interval of the corresponding time node to generate the spectral deviation index for each time node; The spectral deviation index at consecutive time points is weighted and averaged using a sliding window mechanism to generate a dynamic deviation threshold curve. Identify the time periods in which the spectral deviation index continuously exceeds the dynamic deviation threshold curve, and extract the start and end time nodes as candidate abnormal fluctuation periods; The candidate abnormal fluctuation periods are timestamped with the time series of meteorological disaster events, and the percentage of overlapping time is calculated. Based on the continuous growth rate of the overlap time ratio and the spectral deviation index, a pseudo-anomaly fluctuation discrimination function is established to eliminate candidate periods caused by meteorological interference and output the spectral anomaly fluctuation period corresponding to the actual changes in cultivated land.

5. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 1, characterized in that, The process of combining time series analysis to perform change detection processing on the spectral feature sequence and the texture feature sequence generates a set of farmland change patch data for the target area, including: A time series analysis model is constructed by inputting the spectral feature sequence and the texture feature sequence into the time series analysis model in chronological order. The time series analysis model uses a sliding window mechanism to dynamically capture the feature evolution trend. In the time series analysis model, the spectral feature difference and texture feature variability between adjacent time nodes are calculated. The spectral feature difference is measured by Euclidean distance to measure the change in band reflectance, and the texture feature variability is quantified by the time derivative of the contrast index to quantify the rate of change of texture structure. The spectral feature difference degree and texture feature variation degree are standardized respectively to generate a normalized difference index and a normalized variation index. When both the normalized difference index and the normalized variation index exceed the preset change detection threshold, the corresponding pixel is marked as a change candidate region. Spatial clustering is performed on all candidate change regions, and the density peak clustering algorithm is used to merge adjacent pixels to form initial change patches. The density peak clustering algorithm determines the initial cluster center by calculating the number of adjacent pixels of each pixel within a preset neighborhood range. The initial changed patches are subjected to morphological filtering, and erosion operation is performed sequentially to eliminate isolated noise points and dilation operation to smooth the patch boundaries, generating the farmland changed patch data set. The size of the structuring element of the morphological filter is adaptively selected according to the average area of ​​the patches.

6. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 5, characterized in that, The method further includes: The spatial distribution features and temporal frequency features of the cultivated land change patch data set are extracted. The spatial distribution features are quantified by kernel density estimation, and the temporal frequency features are calculated based on the frequency and duration of patch occurrence. The geographic weight parameter of the change detection threshold is adjusted according to the spatial distribution characteristics; The time window length parameter of the time series analysis model is optimized based on the time change frequency characteristics. The time window length is shortened in the changing region to capture changing events, and the time window length is extended in the stable region. When the time window length exceeds the preset maximum window, the starting point of the time window is reset to the current time node minus a fixed period.

7. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 6, characterized in that, The geographical weight parameter for adjusting the change detection threshold based on the spatial distribution characteristics includes: The target area is divided into multiple geographic grid units; The degree of clustering of map features within each geographic grid cell is statistically analyzed. The degree of clustering is calculated by combining the number of map features per unit area with the average area of ​​map features. A regression model is established between the elevation and slope data of geographic grid units and the density of changing patches. The regression model uses the ridge regression algorithm to handle the multicollinearity problem. The theoretical change density is predicted under different geographical conditions based on the regression model, and the theoretical change density reflects the natural constraint effect of topographic factors on farmland change. Calculate the residual value between the actual change density and the theoretical change density, and dynamically adjust the geographic weight parameter of the corresponding geographic grid unit based on the residual value. Decrease the geographic weight parameter of the change detection threshold of the positive residual area according to the first dynamic mapping coefficient, and increase the geographic weight parameter of the change detection threshold of the negative residual area according to the second dynamic mapping coefficient.

8. The intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation according to claim 5, characterized in that, The step of spatial clustering all candidate change regions and merging adjacent pixels using the density peak clustering algorithm to form initial change patches includes: The density peak clustering algorithm is used to perform spatial density analysis on the change candidate region. The density peak clustering algorithm is used to identify the spatial clustering region of the change candidate points and calculate the candidate point density of the clustering region to which each pixel belongs. Based on the eight-neighbor connectivity criterion, spatially clustered pixels identified by the density peak clustering algorithm are labeled with connected components. The density peak clustering algorithm divides clustered regions according to a preset density threshold. The connected component labeling recursively merges adjacent pixels through a seed filling algorithm. The spatially clustered pixels are pixels whose candidate point density is greater than a set density. Calculate the geometric center coordinates and minimum bounding rectangle range of the connected component for each marker. The geometric center coordinates are determined based on the weighted average of the coordinates of all cells in the connected component, and the minimum bounding rectangle is obtained by the rotating caliper algorithm. Merge connected components with overlapping bounding rectangles and recalculate their geometric properties. The merge operation is triggered based on the ratio of the area of ​​the rectangle intersection to the original area. The merged connected domains are filtered by area, and connected domains with an area smaller than the minimum patch area threshold are removed. The minimum patch area threshold is dynamically adjusted based on the landscape fragmentation index calculated based on the ratio of the number of cultivated land patches in the region to the total cultivated land area and the monitoring accuracy requirements. The remaining connected domains are retained as initial changed patches.

9. A system for intelligent identification of farmland change patches based on remote sensing AI intelligent interpretation, characterized in that, The method includes a processor and a memory, the memory and the processor being connected. The memory is used to store programs, instructions or code, and the processor is used to execute the programs, instructions or code in the memory to implement the intelligent identification method for farmland change patches based on remote sensing AI intelligent interpretation as described in any one of claims 1-8.

Images