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

Abstract

The invention provides a remote sensing AI intelligent interpretation-based cultivated land change pattern spot intelligent identification method and system, and the method comprises the steps: firstly obtaining a multi-temporal remote sensing image data set of a target region, comprehensively recording the cultivated land conditions of the target region in different time periods, carrying out the preprocessing operation of the obtained multi-temporal remote sensing image data set, and carrying out the recognition of the cultivated land change pattern spots. Then feature extraction operation is carried out on the preprocessed multi-temporal remote sensing image data set, a spectral feature sequence and a texture feature sequence of a target area are generated respectively, farmland features are described in detail from different angles, and finally, a time sequence analysis method is combined to analyze the farmland features. And performing change detection processing on the generated spectral feature sequence and texture feature sequence to generate a cultivated land change pattern spot data set of the target area, thereby effectively improving accuracy and efficiency of cultivated land change pattern spot identification.

Description

Technical Field

[0001] The present invention relates to the field of artificial intelligence technology, and in particular to a method and system for intelligently identifying cultivated land change patterns based on remote sensing AI intelligent interpretation. Background Art

[0002] In today's agricultural resource management and land resource monitoring, accurate and efficient monitoring of cultivated land changes is crucial. Relevant technologies rely on accurate information on cultivated land changes to rationally plan grain planting areas and adjust agricultural production strategies, thereby ensuring a stable food supply to meet the needs of a growing population. At the same time, sound land resource planning relies on a clear understanding of cultivated land changes. By understanding the increase, decrease, and distribution of cultivated land, we can scientifically arrange the ratio of various types of construction land to agricultural land, thereby improving the efficiency of land resource utilization. Furthermore, cultivated land changes are closely related to ecological and environmental protection. Irrational cultivated land development and utilization can lead to ecological problems such as soil erosion and soil degradation. Therefore, accurate monitoring of cultivated land changes can help to implement appropriate ecological protection measures.

[0003] There are two main traditional methods for monitoring cultivated land changes. One is manual field surveys. This method requires significant manpower, material resources, and time. Investigators must personally visit each cultivated field site to conduct on-site surveys, measurements, and records. Furthermore, the scope of manual field surveys is often limited by manpower and time, making it difficult to conduct comprehensive and rapid monitoring of large areas. Furthermore, survey results are subject to subjective judgment by investigators. Different investigators may reach different conclusions due to differences in experience, knowledge, and judgment criteria, making it difficult to ensure the accuracy and reliability of monitoring results.

[0004] Another traditional method is simple visual interpretation of remote sensing imagery. While this method utilizes remote sensing imagery to obtain cultivated land information, it requires a high level of expertise and practical experience to accurately identify cultivated land and its changes from remote sensing imagery. Furthermore, visual interpretation is a tedious and time-consuming process, requiring interpreters to examine remote sensing imagery piece by piece and region by region, making it highly inefficient. Furthermore, visual interpretation is prone to omissions and misjudgments. Even subtle cultivated land changes are difficult to accurately identify due to factors such as image resolution and visual fatigue. Furthermore, traditional methods often fail to fully exploit the time series information inherent in multi-temporal remote sensing imagery data. Cultivated land changes are a dynamic process, with inherent connections and patterns of change between image data at different time points. Traditional methods fail to effectively utilize this information, resulting in insufficient dynamic monitoring capabilities for cultivated land changes and an inability to accurately and timely grasp their real-time status. Summary of the Invention

[0005] In view of the above-mentioned problems, in combination with the first aspect of the present invention, an embodiment of the present invention provides a method for intelligently identifying cultivated land change patches based on remote sensing AI intelligent interpretation, the method comprising: Acquire a multi-temporal remote sensing image data set of a target area, wherein the multi-temporal remote sensing image data set includes remote sensing image data collected at different time points; Performing a preprocessing operation on the multi-temporal remote sensing image data set to obtain a preprocessed multi-temporal remote sensing image data set; Performing a feature extraction operation on the preprocessed multi-temporal remote sensing image data set to generate a spectral feature sequence and a texture feature sequence of the target area; The spectral feature sequence and the texture feature sequence are subjected to change detection processing in combination with time series analysis to generate a cultivated land change patch data set of the target area.

[0006] On the other hand, an embodiment of the present invention also provides an intelligent identification system for cultivated land change patterns 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 codes, and the processor is used to execute the programs, instructions or codes in the machine-readable storage medium to implement the above method.

[0007] Based on the above aspects, the embodiments of the present invention significantly improve the efficiency and accuracy of cultivated land change monitoring. By obtaining a multi-temporal remote sensing image data set containing data collected at different time nodes in the target area, covering the cultivated land information status of the target area at different time stages, and performing preprocessing operations on the obtained multi-temporal remote sensing image data set, it is possible to effectively remove noise, errors and other interference factors in the data, improve the quality and availability of the data, and the preprocessed data is more in line with the requirements of subsequent feature extraction and analysis, ensuring the accuracy and reliability of the entire recognition process. Feature extraction operations are performed on the preprocessed multi-temporal remote sensing image data set to generate spectral feature sequences and texture feature sequences of the target area, and a comprehensive and detailed description of the cultivated land characteristics of the target area is performed from different dimensions. Spectral features can reflect the reflectance characteristics of cultivated land in different bands, while texture features can reflect the structure and texture information of the cultivated land surface. The combination of the two can more accurately depict the characteristics of cultivated land and improve the ability to identify and distinguish cultivated land. Combined with time series analysis, the generated spectral and texture feature sequences are processed for change detection, fully utilizing the temporal dimension of multi-temporal data to keenly capture subtle changes in cultivated land between different time points. By analyzing time series data, various types of changes, such as additions, reductions, and changes in land use, can be accurately identified, generating a dataset of cultivated land change patches in 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 limited dynamic monitoring capabilities. This method can provide more scientific, accurate, and timely information on cultivated land changes for agricultural resource management and land planning departments. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] Figure 1 It is a schematic diagram of the execution flow of the method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation provided by an embodiment of the present invention.

[0009] Figure 2 This is a schematic diagram of exemplary hardware and software components of the intelligent identification system for cultivated land change patterns based on remote sensing AI intelligent interpretation provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0010] The present invention will be described in detail below with reference to the accompanying drawings. Figure 1 This is a flow chart of a method for intelligently identifying cultivated land change patterns based on remote sensing AI intelligent interpretation provided by an embodiment of the present invention. The following is a detailed introduction to the method for intelligently identifying cultivated land change patterns based on remote sensing AI intelligent interpretation.

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

[0012] In this embodiment, to achieve intelligent identification of cultivated land change patterns in a target area, a multi-temporal remote sensing image dataset of the target area is first acquired. This multi-temporal remote sensing image dataset relies on remote sensing sensors carried by remote sensing satellites, aircraft, and other vehicles, and data is collected from the target area at different time points. For example, given a target area A, to fully understand the dynamic changes in cultivated land there, remote sensing image data must be acquired at multiple different time points. The selection of these time points should fully consider the characteristic changes of cultivated land across different growth cycles and seasons, such as collecting data at key stages such as the sowing period, growing period, and harvest period.

[0013] During data collection, remote sensing sensors record the reflectance information of the target area at different wavelengths, forming multi-band remote sensing image data. Each time point in the remote sensing image data can be considered an independent dataset, containing the spectral and spatial information of the objects within the target area at that moment. These datasets collected at different time points are combined to form a multi-temporal remote sensing image dataset. Each data subset in this multi-temporal remote sensing image dataset corresponds to a specific point in time.

[0014] Step S120: performing a preprocessing operation on the multi-temporal remote sensing image data set to obtain a preprocessed multi-temporal remote sensing image data set.

[0015] After acquiring a multi-temporal remote sensing image data set, these raw data may have radiometric distortion, geometric distortion, noise interference and other problems, which will affect the subsequent feature extraction and change detection results. Therefore, they need to be preprocessed to improve the quality and usability of the data. The preprocessing operation mainly includes radiometric correction, geometric correction, registration and denoising steps, as follows: Step S121: calling a radiation correction algorithm to perform radiation correction processing on each scene of remote sensing image data in the multi-temporal remote sensing image data set.

[0016] During the acquisition of remote sensing image data, factors such as atmospheric scattering and absorption, as well as the sensor's inherent response characteristics, can lead to deviations between the image's radiometric value and the actual reflectivity of the ground. To eliminate this deviation, a radiometric correction algorithm is required to perform radiometric correction on each remote sensing image.

[0017] Assume that for a scene I in a multi-temporal remote sensing image dataset, its original radiometric value is R. The core concept of the radiometric correction algorithm 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 the ground object. Common radiometric correction algorithms include those based on radiation transfer models and those based on ground control points.

[0018] Taking the correction method based on the radiation transfer model as an example, this correction method needs to consider factors such as the optical properties of the atmosphere, the solar altitude angle, and the sensor's observation angle. First, based on the optical parameters of the atmosphere, such as atmospheric transmittance and atmospheric scattering coefficient, an atmospheric radiation transfer model is established. Then, this model is used to calculate the impact of the atmosphere on the radiation value of the remote sensing image, and then the original radiation value is corrected. Specifically, for each pixel in image I, its corrected radiation value R' can be calculated using the following steps: The first step is to obtain the atmospheric optical parameters at the pixel location, including atmospheric transmittance τ, atmospheric scattering coefficient σ, etc. These parameters can be obtained through atmospheric sounding data or meteorological data.

[0019] In the second step, the propagation path length L of solar radiation in the atmosphere is calculated based on the solar altitude angle θ and the sensor's observation angle φ.

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

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

[0022] By performing such radiation correction processing on each image in the multi-temporal remote sensing image data set, the radiation amounts of image data at different time points can be made comparable.

[0023] Step S122: geometrically correcting the radiometrically corrected remote sensing image data to eliminate image spatial position deviations by matching ground control points, wherein the ground control points are automatically extracted based on high-precision topographic maps and remote sensing image feature corner points.

[0024] Remote sensing image data may still have geometric distortion after radiometric correction. This is caused by factors such as sensor attitude changes, terrain undulations, and the curvature of the earth. To eliminate these geometric distortions, the image needs to be geometrically corrected.

[0025] The key to geometric correction is to find ground control points (GCPs). These control points are points on the image with clear geographical locations. 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 the spatial position deviation of the image.

[0026] In this embodiment, ground control points are automatically extracted based on high-precision topographic maps and remote sensing image feature corner points. The specific steps are as follows: First, points with distinct features, such as road intersections and building corners, are extracted from a high-precision topographic map and used as control points on the topographic map. At the same time, feature extraction algorithms, such as the Harris corner detection algorithm, are used to extract feature corners from the remote sensing image.

[0027] The extracted corner points of the remote sensing image are then matched with the control points on the topographic map. This matching process can use a matching method based on feature descriptors, such as SIFT (Scale-Invariant Feature Transform) or SURF (Speeded Up Robust Features). The best matching point pairs are found by calculating the similarity between the feature descriptors of the corner points and the control points.

[0028] Assume that the set of feature corner points extracted from the remote sensing image is C1, and the set of control points extracted from the topographic map is C2. For each corner point c1 in C1, calculate the similarity s between its feature descriptor and that of each control point c2 in C2. Select the control point with the greatest similarity as the matching point, forming a matching point pair (c1, c2).

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

[0030] Through geometric correction processing, remote sensing image data at different time nodes can be made consistent in spatial position, providing accurate spatial reference for subsequent alignment and change detection.

[0031] Step S123: performing a registration operation on the geometrically corrected remote sensing image data to unify the remote sensing image data of different time nodes into the same geographic coordinate system. The registration operation uses a feature matching algorithm to align the image edge with the contour of the object.

[0032] Although the remote sensing image data has been corrected in spatial position after geometric correction, there may still be certain deviations in the images at different time nodes, and alignment operations are required to unify them into the same geographic coordinate system.

[0033] The main purpose of the registration operation is to align remote sensing image data at different time points so that the same ground objects in the image overlap in spatial position. In this embodiment, a feature matching algorithm is used to achieve image registration. The specific steps are as follows: 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, a feature extraction algorithm, such as the ORB (Oriented FAST and Rotated BRIEF) algorithm, is used to extract feature points from the image.

[0034] Assume that 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 the similarity s between its feature descriptor and that of each feature point f_base in F_base. Select the feature point with the greatest similarity as the matching point, forming a matching point pair (f_reg, f_base).

[0035] To improve matching accuracy, the RANSAC (Random Sample Consensus) algorithm can also be used to screen matching point pairs and remove incorrectly matched point pairs. The RANSAC algorithm randomly selects a portion of matching point pairs, calculates a geometric transformation model, and then verifies the other matching point pairs based on the model, retaining the point pairs that meet the model.

[0036] Based on the screened matching point pairs, the geometric transformation matrix T between the image I_reg and the reference image I_base is calculated. The geometric transformation matrix T can be an affine transformation matrix or a perspective transformation matrix, which describes the spatial transformation relationship between the image I_reg and the reference image I_base.

[0037] 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 the 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).

[0038] Through the registration operation, remote sensing image data at different time points can be unified in the geographic coordinate system, which facilitates subsequent comparison and analysis of the images.

[0039] Step S124: performing denoising processing on the registered remote sensing image data to obtain a pre-processed multi-temporal remote sensing image data set.

[0040] The remote sensing image data after registration may still contain noise, which will affect the subsequent feature extraction and change detection results, so it needs to be denoised.

[0041] Common denoising methods include mean filtering, median filtering, Gaussian filtering, etc. Taking median filtering as an example, median filtering is a nonlinear 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 the pixel.

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

[0043] The specific steps are as follows: The first step is to determine the neighborhood range N for each pixel P in the image I.

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

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

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

[0047] By performing such median filtering on each pixel in the registered remote sensing image data, the noise in the image can be effectively removed, and a preprocessed multi-temporal remote sensing image data set can be obtained.

[0048] Step S130: performing a feature extraction operation on the pre-processed multi-temporal remote sensing image data set to generate a spectral feature sequence and a texture feature sequence of the target area.

[0049] The preprocessed multi-temporal remote sensing image dataset contains rich information about the target area. However, in order to more accurately identify cultivated land change patches, useful features need to be extracted from it. This step mainly extracts spectral feature sequences and texture feature sequences. The specific steps are as follows: Step S131: performing band synthesis processing on the pre-processed remote sensing image data of each scene to generate multispectral image data including visible light and near-infrared bands. The band synthesis processing dynamically selects band combinations according to the vegetation index calculation requirements.

[0050] In a 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 both visible and near-infrared bands.

[0051] Band compositing requires dynamic selection of band combinations based on the vegetation index calculation requirements. Different vegetation indices have different band selection requirements. 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.

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

[0053] The specific band synthesis process can be achieved through the following steps: 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.

[0054] The second step is to combine the data of each band in Bselect according to the same pixel position 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 the pixel in Imulti.

[0055] Through band synthesis processing, multispectral image data containing visible light and near-infrared bands can be generated, providing a basis for subsequent spectral feature and texture feature extraction.

[0056] Step S132: Calculating the mean reflectance of each pixel in different bands based on the multispectral image data to generate the spectral feature sequence of the target area. The mean reflectance is normalized to eliminate the influence of different lighting conditions.

[0057] After obtaining multispectral image data containing visible light and near-infrared bands, in order to generate the spectral feature sequence of the target area, it is necessary to calculate the mean reflectance of each pixel in different bands and eliminate the influence of different lighting conditions through normalization. This process involves a series of detailed steps, as follows: Step S1321: The land cover types of the target area are divided according to a preset land feature classification system. The land cover types are automatically labeled using a supervised classification algorithm combined with training samples. The land cover types include at least three types: cultivated land, forest land, and water area.

[0058] The preset land feature classification system is constructed based on the spectral characteristics, spatial distribution characteristics and other related attributes of the land features. In order to accurately classify the land cover types in the target area, it is necessary to use a supervised classification algorithm and automatically label the training samples.

[0059] The first step is to select training samples. These samples need to cover a wide range of typical landform types within the target area. For cultivated land, pixels representing different crops and growth stages should be selected. For woodlands, 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 should be representative and accurately reflect the characteristics of each type of landform.

[0060] Take the maximum likelihood classification algorithm, for example. This is a commonly used supervised classification algorithm, whose core concept is based on the Bayesian criterion in probability theory. For each pixel in a multispectral image, the algorithm calculates the probability of that pixel belonging to each land cover type based on the training samples, and then classifies it into the category with the highest probability.

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

[0062] For example, let the pixel's eigenvector be A, the eigenmean vector of a particular land cover type be B, and the covariance matrix be C. The process for calculating the Mahalanobis distance D from pixel A to that land cover type is: first calculate the difference vector E between A and B, then multiply the transpose of E by the inverse matrix of C, and then 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 the pixel.

[0063] In this way, all pixels in the target area are classified to complete the automatic labeling of land cover types, and at least three types of land cover are divided: cultivated land, forest land, and water area.

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

[0065] After completing the land cover classification, for each land cover type, the mean reflectance of the corresponding pixels in the visible and near-infrared bands is calculated. 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.

[0066] The basic principle of the region growing algorithm is to start from one or more seed pixels and gradually merge pixels that have similar characteristics to the seed pixels and are spatially adjacent to them into the same region.

[0067] The first step is to select seed pixels. For each land cover type, some pixels can be randomly selected as seed pixels, or representative pixels can be selected as seed pixels based on their characteristic values ​​(such as reflectance).

[0068] Then, starting with the seed pixel, the pixels within its neighborhood are examined. Typically, an 8-neighborhood (i.e., the 8 pixels adjacent to the seed pixel) is used for this examination. For a pixel within the neighborhood, if it belongs to the same land cover type as the seed pixel and the difference between its reflectance in the visible and near-infrared bands and the reflectance of the seed pixel is within a preset threshold, the neighboring pixel is merged into the region where the seed pixel resides.

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

[0070] Repeat this process and expand the area continuously until there are no neighboring pixels that meet the conditions. In this way, the aggregated similar pixel area is obtained.

[0071] For each aggregated area, the sum of the reflectances of all pixels in the area in the visible light band and the near-infrared band is calculated, and then divided by the number of pixels in the area to obtain the mean reflectance of the area in the visible light band and the near-infrared band.

[0072] Step S1323: constructing a spectral change curve for each land cover type based on the reflectivity mean values ​​at different time nodes, wherein the spectral change curve is fitted with the reflectivity estimation values ​​of the missing time nodes using a cubic spline interpolation method.

[0073] After obtaining the mean reflectance values ​​of the visible and near-infrared bands for each land cover type at different time points, we need to construct a spectral change curve for each land cover type. The spectral change curve can intuitively reflect the change in reflectance of the land cover type over time.

[0074] Since there may be missing data at certain time nodes during the actual data collection process, in order to make the spectral change curve more continuous and accurate, the cubic spline interpolation method is used to fit the reflectance estimation value of the missing time nodes.

[0075] Cubic spline interpolation is a piecewise polynomial interpolation method that approximates the original data by constructing a cubic polynomial on each small interval. Specifically, for a given time node and its corresponding mean reflectivity, the interval between adjacent time nodes is divided into several small intervals. In each small interval, a cubic polynomial is constructed such that the polynomial is equal to the known mean reflectivity at the endpoints of the small interval and has a second-order continuous derivative over the entire interval.

[0076] For example, given that the mean reflectivity values ​​corresponding to time nodes t1, t2, and t3 are R1, R2, and R3, respectively, a cubic polynomial P1(t) = a1*t^3 + b1*t^2 + c1*t + d1 is constructed for the small interval t1 to t2, where a1, b1, c1, and d1 are the coefficients to be determined. These coefficients are determined by requiring P1(t1) = R1 and P1(t2) = R2, and by ensuring that the first and second derivatives of P1(t) at t1 and t2 are equal to the first and second derivatives of the polynomial for the adjacent small interval at the corresponding endpoints.

[0077] By performing this process on all small intervals, the cubic spline interpolation function for the entire time interval is obtained. For the reflectivity estimate of the missing time node, the time node is substituted into the corresponding cubic spline interpolation function for calculation.

[0078] In this way, the spectral change curves of each land cover type are constructed, making the reflectance data more continuous and complete in the temporal dimension.

[0079] Step S1324: Based on the spectral change curve, identify the spectral abnormal fluctuation period of the cultivated land type pixel, use the mutation point detection algorithm to locate the time node where the reflectance change amplitude within the set time interval is greater than the set amplitude, generate the time series change mark in the spectral feature sequence, perform time intersection operation on the spectral abnormal fluctuation period and the cultivated land change pattern appearance period, and retain the pattern whose overlapping period exceeds the preset threshold.

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

[0081] The basic idea of ​​the mutation point detection algorithm is to find the time point when the reflectivity changes significantly by comparing the reflectivity changes at adjacent time nodes or within a set time interval.

[0082] For example, set a time interval of T and calculate the reflectance difference between each time node and the previous T time nodes for the spectral change curve of the cultivated land type pixel. If the reflectance difference at a time node is greater than the set amplitude threshold, the time node is considered to be a mutation point.

[0083] After the mutation points are determined, the time periods between adjacent mutation points are defined as spectral abnormal fluctuation periods. These spectral abnormal fluctuation periods are marked as time series change marks and added to the spectral feature sequence.

[0084] Then, a temporal intersection operation is performed on the periods of spectral anomaly fluctuation and the periods of cropland change patches. The temporal intersection operation involves finding the portion of time that overlaps between each period of spectral anomaly fluctuation and each period of cropland change patches. The proportion of this overlapping time in the total number of periods of spectral anomaly fluctuation and cropland change patches is then calculated.

[0085] Only when the proportion of overlapping periods in the periods of spectral abnormal fluctuations or the periods of cultivated land change patches exceeds a preset threshold will the corresponding cultivated land change patches be retained. This can filter out the true cultivated land change patches related to spectral abnormal fluctuations and improve the accuracy of cultivated land change patch identification.

[0086] Step S13241: extracting historical reflectance data of cultivated land type pixels in the growing season and non-growing season, performing time dimension normalization processing on the historical reflectance data, and generating a normalized reflectance baseline interval for each time node.

[0087] To more accurately identify periods of abnormal spectral fluctuations in cultivated land pixels, it is necessary to extract historical reflectance data for cultivated land pixels during the 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 in both normal growing and non-growing states.

[0088] When extracting historical reflectance data, it is important to ensure its accuracy and completeness. Reflectance data for cultivated land pixels can be filtered from historical multispectral imagery and categorized by growing and non-growing seasons.

[0089] The extracted historical reflectance data are normalized in the time dimension to eliminate the influence of factors such as lighting conditions and atmospheric conditions between different time nodes and make the data comparable.

[0090] Time dimension normalization can be performed using the minimum-maximum normalization method. For each time node's reflectance data, find the minimum and maximum values ​​within the dataset. For a reflectance value R, its normalized reflectance value R' is calculated by first calculating the difference between R and the minimum value, then dividing it by the difference between the maximum and minimum values.

[0091] By performing this normalization process on the reflectivity data at all time points, normalized reflectivity data is obtained. Based on this normalized reflectivity data, the reflectivity mean and standard deviation for each time point are calculated. The normalized reflectivity baseline interval for each time point is determined based on the set standard deviation multiple, centered around the reflectivity mean. For example, the range of the mean plus or minus two standard deviations can be used as the baseline interval.

[0092] Step S13242: Obtain the multispectral image reflectance data of the current time node, and use the same normalization processing method as the historical reflectance data to generate a normalized reflectance vector of the current time node.

[0093] Obtain the multispectral image reflectance data of the current time node. These data reflect the reflectance of the cultivated land at the current moment.

[0094] The multispectral image reflectance data at the current time node is normalized using the same normalization method used for historical reflectance data. This involves finding the minimum and maximum values ​​in the reflectance dataset at the current time node and calculating the normalized value for each reflectance value using the minimum-maximum normalization method.

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

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

[0097] The standardized Euclidean distance between the normalized reflectance vector of the current time node and the baseline interval of the corresponding time node is calculated to generate the spectral deviation index of each time node.

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

[0099] The process for calculating the standardized Euclidean distance is as follows: First, calculate the difference vector D between the normalized reflectivity 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', to obtain the standardized Euclidean distance.

[0100] This normalized Euclidean distance is used as the spectral deviation index for that time node. The larger the spectral deviation index, the greater the difference between the reflectance at the current time node and the baseline interval, and the possibility of spectral anomaly.

[0101] Step S13244: Perform weighted averaging processing on the spectrum deviation indicators of consecutive time nodes based on the sliding window mechanism to generate a dynamic deviation threshold curve.

[0102] Based on the sliding window mechanism, the spectral deviation indicators of consecutive time nodes are weighted averaged to generate a dynamic deviation threshold curve.

[0103] The sliding window mechanism is to set a fixed-size window on the time series, and the window moves gradually on the time series. The spectral deviation index within each window is processed using a weighted average method.

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

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

[0106] Connecting these weighted average spectral deviation index values ​​forms a dynamic deviation threshold curve. This curve can dynamically adjust the threshold for judging spectral anomalies based on time, improving the accuracy of identifying abnormal spectral fluctuations.

[0107] Step S13245: Identify the time period in which the spectrum deviation index continuously exceeds the dynamic deviation threshold curve, and extract its starting time node and ending time node as candidate abnormal fluctuation time periods.

[0108] After the dynamic deviation threshold curve is obtained, a time period in which the spectrum deviation index continuously exceeds the dynamic deviation threshold curve is identified.

[0109] Starting from the starting point of the time series, the spectral deviation index of each time node is checked one by one. If the spectral deviation index of a time node exceeds the threshold corresponding to the dynamic deviation threshold curve, and the spectral deviation index of several subsequent time nodes continues to exceed the threshold, the period of continuous exceeding the threshold is recorded.

[0110] The start and end time nodes of this period are extracted and used as candidate abnormal fluctuation periods. These candidate abnormal fluctuation periods indicate that the reflectivity of the cultivated land has continuously deviated from the normal range for a period of time, which may be an abnormal situation.

[0111] Step S13246: Match the candidate abnormal fluctuation period with the time series of meteorological disaster event records by timestamp, and calculate the overlapping time ratio.

[0112] The candidate abnormal fluctuation periods are timestamp-matched with the time series of meteorological disaster event records in order to exclude the spectral abnormal fluctuations caused by natural factors such as meteorological disasters and find out the spectral abnormal fluctuations that are actually caused by cultivated land changes.

[0113] The time series of meteorological disaster event records contains the time information of meteorological disasters in the target area. For each candidate abnormal fluctuation period, check whether it overlaps with the time interval of each disaster event in the time series of meteorological disaster event records.

[0114] The process for calculating the overlap time percentage is as follows: for a candidate abnormal fluctuation period and a meteorological disaster time interval, find the length of their overlap. Then, divide the overlap time length by the length of the candidate abnormal fluctuation period and the length of the meteorological disaster time interval, respectively, to obtain two overlap time percentages.

[0115] Step S13247: According to the overlapping time ratio and the continuous growth rate of the spectral deviation index, a pseudo-abnormal fluctuation discrimination function is established to eliminate the candidate time periods caused by meteorological interference and output the spectral abnormal fluctuation time periods corresponding to the real cultivated land changes.

[0116] Based on the overlap time ratio and the sustained growth rate of the spectral deviation index, a pseudo-abnormal fluctuation discriminant function was established. The purpose of the pseudo-abnormal fluctuation discriminant function is to determine whether the candidate abnormal fluctuation period is caused by meteorological interference or real farmland changes.

[0117] The sustained growth rate of the spectral deviation index can be obtained by calculating the difference between the spectral deviation indicators of adjacent time nodes and dividing it by the time interval.

[0118] The pseudo-anomalous fluctuation discrimination function can be a function that comprehensively considers the overlap time percentage and the sustained growth rate of the spectral deviation index. For example, a threshold combination can be set. When the overlap time percentage exceeds a certain threshold and the sustained growth rate of the spectral deviation index falls below another threshold, the candidate abnormal fluctuation period is considered to be caused by meteorological interference and is eliminated.

[0119] By performing such judgment and screening on all candidate abnormal fluctuation periods, the candidate periods caused by meteorological interference are eliminated, and finally the spectral abnormal fluctuation periods corresponding to the real cultivated land changes are output.

[0120] Therefore, based on the multispectral image data, the mean reflectance of each pixel in different bands can be accurately calculated to generate the spectral feature sequence of the target area, and effectively eliminate the influence of differences in lighting conditions. At the same time, the spectral abnormal fluctuation period of the cultivated land type pixels can be accurately identified, providing reliable spectral feature information for the subsequent identification of cultivated land change spots.

[0121] 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 the object, and the homogeneity index is used to quantify the uniformity of the gray level distribution in the local area.

[0122] In addition to spectral features, texture features are also important for identifying cultivated land change patches. This step extracts the gray-level co-occurrence matrix from the multispectral image data and calculates the contrast and homogeneity index of each pixel based on this matrix to generate a texture feature sequence for the target area.

[0123] First, the multispectral image data Imulti is converted into a grayscale image Igrey. The weighted average method can be used to convert the multiple band data 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 by the following formula: G(x, y)=w1*R1(x, y)+w2*R2(x, y)+...+wm*Rm(x, y) Among them, w1, w2, ..., wm are the weight coefficients of each band, and they satisfy w1+w2+...+wm=1.

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

[0125] Next, the contrast and homogeneity index of each pixel are calculated based on the gray level co-occurrence matrix. The contrast index C is used to characterize the edge sharpness of the object, and its calculation formula is: C=∑(i,j)[(ij)^2*GLCM(i,j;d,θ)] The homogeneity index H is used to quantify the uniformity of grayscale distribution in a local area, and its calculation formula is: H=∑(i,j)[GLCM(i,j;d,θ) / (1+(ij)^2)] Finally, the contrast and homogeneity indices of all pixels in the target area are arranged in the 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.

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

[0127] Step S134: performing a spatial alignment operation on the spectral feature sequence and the texture feature sequence at the same time node to generate a spectral feature sequence and a texture feature sequence of the target area.

[0128] After generating the spectral feature sequence and texture feature sequence respectively, in order to ensure their spatial consistency, the spectral feature sequence and texture feature sequence of the same time node need to be spatially aligned.

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

[0130] Assume that the spectral feature sequence at the same time node 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, the new feature sequence F = [F1, F2, ..., Fk] is generated, where Fi = [Si, Ti]. That is, the feature of each pixel is the concatenation of its spectral and texture features.

[0131] Through the spatial alignment operation, the spectral features and texture features can be organically combined to generate the spectral feature sequence and texture feature sequence of the target area, providing more comprehensive feature information for subsequent change detection.

[0132] Step S140: performing change detection processing on the spectral feature sequence and the texture feature sequence in combination with time series analysis to generate a cultivated land change patch data set of the target area.

[0133] After obtaining the spectral feature sequence and texture feature sequence of the target area, in order to identify the change patterns of cultivated land, it is necessary to combine time series analysis to perform change detection on these two feature sequences. This process includes the following sub-steps: Step S141: constructing a time series analysis model, inputting the spectral feature sequence and the texture feature sequence into the time series analysis model in chronological order, and the time series analysis model adopts a sliding window mechanism to dynamically capture feature evolution trends.

[0134] First, we need to build a time series analysis model. This time series analysis model analyzes the spectral feature sequence and texture feature sequence to capture their changing trends over time. Here, a sliding window mechanism is used, which can dynamically observe the changes in features in the time dimension.

[0135] Assume that the spectral feature sequence is represented by S and the texture feature sequence is represented by T. They are both feature sets arranged in time sequence. The size of the sliding window is represented by w, and the step size of each movement of the window in the time series is represented by s.

[0136] When building a time series analysis model, you can use models suitable for processing time series data, such as long short-term memory (LSTM) networks or gated recurrent units (GRUs). Taking the LSTM as an example, it consists of an input layer, hidden layers, and an output layer. The input layer receives data on spectral and texture feature sequences. The LSTM units in the hidden layer memorize the long-term dependencies in the time series, and the output layer outputs predictions of feature trends.

[0137] 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 features and texture features of that time point and the w-1 time points before it are combined into an input vector and 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)].

[0138] The sliding window mechanism allows the model to dynamically capture the evolving trends of features. As the window moves across the time series, the model continuously updates its analysis of feature changes.

[0139] Step S142: Calculate the spectral feature difference and texture feature variation between adjacent time nodes in the time series analysis model. The spectral feature difference measures the amplitude of the band reflectance change by Euclidean distance, and the texture feature variation quantifies the texture structure change rate based on the time derivative of the contrast index.

[0140] In the time series analysis model, it is necessary to calculate the spectral feature differences and texture feature variations between adjacent time nodes to determine whether the cultivated land has changed.

[0141] For spectral feature difference, the Euclidean distance is used to measure the magnitude of the 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 calculated by calculating the Euclidean distance between these two vectors. The specific calculation process is to first calculate the sum of the squares of the differences between the corresponding elements of the two vectors and then take the square root of this 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 square of (s1(t+1) - s1(t)) plus the square of (s2(t+1) - s2(t)) and so on up to (sn(t+1) - sn(t)), and then take the square root of this sum.

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

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

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

[0145] Since the value ranges of spectral feature difference and texture feature variability may be different, they need to be standardized separately for the convenience of comparison and unified processing.

[0146] The Z-score standardization method can be used for normalization. For the spectral feature difference Ds, its mean μs and standard deviation σs are first calculated. Then, for each spectral feature difference value Ds(i), its normalized difference index NDs(i) can be obtained by dividing (Ds(i)-μs) by σs. For the texture feature variability Dt, its mean μt and standard deviation σt are similarly calculated. The normalized variation index NDt(i) of each texture feature variability value Dt(i) can be obtained by dividing (Dt(i)-μt) by σt.

[0147] 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 characteristics. The pixel is marked as a candidate change area, which can preliminarily screen out pixels where cultivated land changes may occur.

[0148] Step S144: performing spatial clustering processing on all candidate change areas, and using a density peak clustering algorithm to merge adjacent pixels to form an initial change patch. 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.

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

[0150] The algorithm first calculates the number of adjacent pixels for 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 a radius of r is drawn around it. The number of candidate pixels within this neighborhood is counted, and this number is the local density ρ(p) of that pixel.

[0151] Then, the minimum distance δ(p) between each pixel and the pixel with a larger local density is calculated. By comparing the local density ρ(p) and the minimum distance δ(p) of all pixels, the pixels with larger local density and larger minimum distance are found. These pixels are the initial cluster centers.

[0152] For each initial cluster center, the pixels in its neighborhood are merged into the cluster where the cluster center is located to form the initial change patch. The specific merging process is to determine whether a pixel is in the neighborhood of an initial cluster center. If so, the pixel is marked as part of the cluster.

[0153] Step S145: Perform morphological filtering on the initial change patches, perform erosion operations to eliminate isolated noise points and dilation operations to smooth patch boundaries in sequence, and generate the cultivated land change patch data set. The size of the structural element of the morphological filter is adaptively selected according to the average area of ​​the patches.

[0154] The initial change 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 change patches.

[0155] Morphological filtering involves erosion and dilation. The purpose of the erosion operation is to eliminate isolated noise points. An appropriate structuring element is selected, and its size can be adaptively chosen based on the average area of ​​the patch. Assuming the structuring element is represented by SE, for each pixel in the initial change patch, if its neighborhood completely matches the SE, the pixel is retained; otherwise, it is deleted. This operation can remove some isolated, small noise points.

[0156] The purpose of the dilation operation is to smooth the edges of the patch. After the erosion operation, the dilation operation is performed on the patch. Again, using the structuring element SE, for each pixel in the patch, all pixels in its neighborhood that intersect with the structuring element SE are marked as part of the patch. This dilation operation can smooth the edges of the patch.

[0157] After the erosion and dilation operations, the obtained patches are the final cultivated land change patch data set.

[0158] Based on the above examples, the generated cultivated land change patch data set can be subsequently processed or applied, such as for result display and data storage. In this example, it is assumed that the generated cultivated land change patch data set is stored in a database for subsequent query and analysis. The relevant information for each patch in the cultivated land change patch data set, such as location, area, change time, etc., is stored in the corresponding table of the database according to a predefined data format. This facilitates further statistics and analysis of cultivated land changes.

[0159] Furthermore, the method may further include: Step S210: extracting spatial distribution characteristics and temporal change frequency characteristics from the cultivated land change patch data set, wherein the spatial distribution characteristics are quantified by kernel density estimation method to determine the degree of patch aggregation, and the temporal change frequency characteristics are calculated based on the frequency and duration of patch occurrence.

[0160] In order to further optimize the results of change detection, it is necessary to extract the spatial distribution characteristics and temporal change frequency characteristics of the cultivated land change patch data set.

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

[0162] The temporal variation frequency feature is calculated based on the frequency and duration of each patch. The number of times each patch appears in the entire time series is counted, recorded as the frequency N, and the duration T of each patch is also recorded. The frequency N is divided by the duration T to obtain a temporal variation frequency index F, which reflects the frequency of patch changes.

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

[0164] Based on the extracted spatial distribution characteristics, the geographic weight parameter of the change detection threshold can be adjusted to improve the accuracy of change detection. The specific steps are as follows: Step S221: Divide the target area into a plurality of geographic grid units.

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

[0166] Step S222: Counting the degree of clustering of spots in each geographic grid unit, wherein the degree of clustering of spots is calculated by comprehensively calculating the number of spots in a unit area and the average spot area.

[0167] For each geographic grid cell, count the number of patches N and their total area A. Then calculate the number of patches per unit area (N / A) and the average area (A / N). Taking the number of patches per unit area and the average patch area into account, we can obtain a patch aggregation index (C). For example, we can take a weighted sum of the number of patches per unit area and the average patch area to obtain the patch aggregation index (C) = α*(N / A) + β*(A / N), where α and β are weight coefficients, and α+β=1.

[0168] Step S223: establishing a regression relationship model between the altitude and slope data of the geographic grid unit and the density of the change patch. The regression relationship model uses a ridge regression algorithm to deal with the multicollinearity problem.

[0169] Geographic data, such as altitude H and slope S, were collected for each geographic grid cell. The density of change patches, D, within each grid cell was calculated (the density of change patches can be calculated by dividing the number of patches by the grid cell area). A regression model was developed to predict the density of change patches, D, using altitude H and slope S. Because altitude and slope data may contain multicollinearity, a ridge regression algorithm was used. Ridge regression, based on ordinary least squares, incorporates a regularization term to reduce the variance of parameter estimates and improve model stability.

[0170] Step S224: predicting theoretical change densities under different geographical conditions based on the regression relationship model, where the theoretical change density reflects the natural constraint effect of topographic factors on cultivated land changes.

[0171] Using the established regression model, we input altitude and slope data for different geographic grid cells to predict the theoretical change density D' for each grid cell. This theoretical change density reflects the natural constraints of topography on cultivated land change, meaning that the likelihood of cultivated land change varies under different topographic conditions.

[0172] 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 unit 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.

[0173] 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 means that the actual change density of the geographic grid cell is higher than the theoretical change density, and there may be some changes in cultivated land caused by non-topographic factors. Reduce the geographic weight parameter of the change detection threshold of the area according to the first dynamic mapping coefficient k1, so that these changes can be more easily detected. If the residual value R is negative, it means that the actual change density of the geographic grid cell is lower than the theoretical change density. Increase the geographic weight parameter of the change detection threshold of the area according to the second dynamic mapping coefficient k2 to reduce the possibility of false detection.

[0174] 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 area to capture the change events, and extend the time window length in the stable area. When the time window length exceeds the preset maximum window, reset the starting point of the time window to the current time node minus the fixed period.

[0175] According to the extracted time-varying frequency characteristics, the time window length parameters of the time series analysis model can be optimized.

[0176] For the changing areas, that is, the areas with high frequency of temporal changes, shorten the time window length. Because in these areas, cultivated land changes are more frequent, shorter time windows can capture change events more timely. For example, shorten the time window length from the original w to w1 (w1 <w)。

[0177] For stable regions, where the frequency of temporal changes is low, extend the time window length. A longer time window can smooth out small fluctuations and more accurately reflect the characteristics of stable regions. For example, extend the time window length from the original w to w2 (where w2 > w).

[0178] When the time window length exceeds the preset maximum window length Wmax, in order to avoid information redundancy caused by the time window being too long, the starting point of the time window is reset to the current time node minus a fixed period T. This can ensure that the time window is within a reasonable range while continuing to pay attention to changes in cultivated land.

[0179] For data collection involved in the above process, if privacy-sensitive data is involved, differential privacy technology is used to protect privacy and prevent data leakage. The basic idea of ​​differential privacy technology is to add a certain amount of noise to the data, so that the difference in query results before and after the noise is added 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 based on the sensitivity of the data and the privacy budget. This method can effectively prevent the leakage of privacy-sensitive data while ensuring data availability.

[0180] Figure 2 A schematic diagram illustrates exemplary hardware and software components of a system 100 for intelligently identifying cultivated land change patterns based on remote sensing AI intelligent interpretation, which can implement the concepts of this application, as provided in some embodiments of this application. For example, the processor 120 can be used in the system 100 for intelligently identifying cultivated land change patterns based on remote sensing AI intelligent interpretation, and can be used to perform the functions described in this application.

[0181] The system 100 for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation can be a general-purpose server or a special-purpose server, both of which can be used to implement the intelligent identification method of cultivated land change patterns 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.

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

[0183] For ease of explanation, only one processor is described in the intelligent identification system 100 for cultivated land change patterns based on remote sensing AI intelligent interpretation. However, it should be noted that the intelligent identification system 100 for cultivated land change patterns based on remote sensing AI intelligent interpretation in this application may also include multiple processors, so the steps performed by one processor described in this application may also be performed jointly or individually by multiple processors. For example, if the processor of the intelligent identification system 100 for cultivated land change patterns based on remote sensing AI intelligent interpretation executes step A and step B, it should be understood that step A and step B may also be executed jointly by two different processors or individually in one processor. For example, the first processor executes step A, the second processor executes step B, or the first processor and the second processor execute steps A and B together.

[0184] In addition, an embodiment of the present invention also provides a readable storage medium, which has computer-executable instructions preset in the readable storage medium. When the processor executes the computer-executable instructions, the above-mentioned intelligent identification method of cultivated land change patterns based on remote sensing AI intelligent interpretation is implemented.

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

Claims

1. A method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation, characterized in that: The method comprises: Acquire a multi-temporal remote sensing image data set of a target area, wherein the multi-temporal remote sensing image data set includes remote sensing image data collected at different time points; Performing a preprocessing operation on the multi-temporal remote sensing image data set to obtain a preprocessed multi-temporal remote sensing image data set; Performing a feature extraction operation on the preprocessed multi-temporal remote sensing image data set to generate a spectral feature sequence and a texture feature sequence of the target area; The spectral feature sequence and the texture feature sequence are subjected to change detection processing in combination with time series analysis to generate a cultivated land change patch data set of the target area.

2. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 1 is characterized in that: The preprocessing operation is performed on the multi-temporal remote sensing image data set to obtain a preprocessed multi-temporal remote sensing image data set, including: Calling a radiation correction algorithm to perform radiation correction processing on each scene of remote sensing image data in the multi-temporal remote sensing image data set; Performing geometric correction on the radiometrically corrected remote sensing image data to eliminate image spatial position deviations by matching ground control points, wherein the ground control points are automatically extracted based on high-precision topographic maps and remote sensing image feature corner points; Performing a registration operation on the geometrically corrected remote sensing image data to unify the remote sensing image data of different time nodes into the same geographic coordinate system, wherein the registration operation uses a feature matching algorithm to align image edges with ground feature contours; The registered remote sensing image data is subjected to denoising processing to obtain a pre-processed multi-temporal remote sensing image data set.

3. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 1 is characterized in that: The performing of a feature extraction operation on the pre-processed multi-temporal remote sensing image data set to generate a spectral feature sequence and a texture feature sequence of the target area includes: Perform band synthesis processing on the pre-processed remote sensing image data of each scene to generate multispectral image data containing visible light and near-infrared bands. The band synthesis processing dynamically selects band combinations according to the requirements of vegetation index calculation; Calculating the mean reflectivity of each pixel in different bands based on the multispectral image data to generate the spectral feature sequence of the target area, wherein the mean reflectivity is normalized to eliminate the influence of different lighting conditions; Extracting the gray level co-occurrence matrix of the multispectral image data, calculating the contrast and homogeneity index of each pixel based on the gray level co-occurrence matrix, and generating the texture feature sequence of the target area, wherein the contrast index is used to characterize the edge sharpness of the object, and the homogeneity index is used to quantify the uniformity of the gray level distribution in the local area; A spatial alignment operation is performed on the spectral feature sequence and the texture feature sequence at the same time node to generate the spectral feature sequence and texture feature sequence of the target area.

4. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 3 is characterized in that: The step of calculating the reflectance mean of each pixel in different bands based on the multispectral image data to generate the spectral feature sequence of the target area includes: Classifying the land cover types of the target area according to a preset land feature classification system, wherein the land cover types are automatically labeled using a supervised classification algorithm combined with training samples, and the land cover types include at least three types: cultivated land, forest land, and water area; For each land cover type, the mean reflectance of the corresponding pixel in the visible light band and the near-infrared band is calculated respectively. In the process of calculating the mean reflectance, spatially adjacent pixels of the same type are aggregated using a region growing algorithm. Constructing a spectral change curve for each land cover type based on the mean reflectivity values ​​at different time nodes, wherein the spectral change curve is fitted with the reflectivity estimation values ​​of the missing time nodes using a cubic spline interpolation method; Based on the spectral change curve, the spectral abnormal fluctuation period of the cultivated land type pixel is identified, and the mutation point detection algorithm is used to locate the time node when the reflectance change amplitude in the set time interval is greater than the set amplitude, and the time series change mark in the spectral feature sequence is generated. The spectral abnormal fluctuation period and the cultivated land change pattern appearance period are subjected to time intersection operation, and the patterns whose overlapping period exceeds the preset threshold are retained.

5. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 4 is characterized in that: The identifying of the abnormal spectral fluctuation period of the cultivated land type pixel based on the spectral change curve includes: Extracting historical reflectance data of cultivated land type pixels in the growing season and non-growing season, normalizing the historical reflectance data in the time dimension, and generating normalized reflectance baseline intervals at each time node; Obtaining the multispectral image reflectance data of the current time node, and using the same normalization processing method as the historical reflectance data to generate a normalized reflectance vector of the current time node; Calculating the normalized Euclidean distance between the normalized reflectance vector and the baseline interval of the corresponding time node to generate a spectral deviation index for each time node; Based on the sliding window mechanism, the spectral deviation index of consecutive time nodes is weighted averaged to generate a dynamic deviation threshold curve; Identify the period in which the spectral deviation index continuously exceeds the dynamic deviation threshold curve, and extract the start time node and the end time node as candidate abnormal fluctuation periods; Match the candidate abnormal fluctuation period with the time series of meteorological disaster event records by timestamp, and calculate the overlapping time ratio; According to the overlapping time proportion and the continuous growth rate of the spectral deviation index, a pseudo-abnormal fluctuation discriminant function is established to eliminate the candidate time periods caused by meteorological interference and output the spectral abnormal fluctuation time periods corresponding to the real cultivated land changes.

6. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 1 is characterized in that: The combining of time series analysis to perform change detection processing on the spectral feature sequence and the texture feature sequence to generate a cultivated land change patch data set of the target area includes: Constructing a time series analysis model, inputting the spectral feature sequence and the texture feature sequence into the time series analysis model in chronological order, and the time series analysis model adopts a sliding window mechanism to dynamically capture feature evolution trends; Calculating the spectral feature difference and texture feature variation between adjacent time nodes in the time series analysis model, wherein the spectral feature difference measures the amplitude of the band reflectance change by Euclidean distance, and the texture feature variation quantifies the rate of change of the texture structure based on the time derivative of the contrast index; The spectral feature difference and texture feature variation are respectively normalized to generate a normalized difference index and a normalized variation index. When both the normalized difference index and the normalized variation index exceed a preset change detection threshold, the corresponding pixel is marked as a candidate change area. Perform spatial clustering on all candidate change areas, and use a density peak clustering algorithm to merge adjacent pixels to form an initial change map. The density peak clustering algorithm determines the initial cluster center by calculating the number of adjacent pixels of each pixel within a preset neighborhood. The initial change spots are subjected to morphological filtering, and the corrosion operation is performed in sequence to eliminate isolated noise points and the dilation operation is performed to smooth the spot boundaries to generate the cultivated land change spot data set. The size of the structural element of the morphological filtering is adaptively selected according to the average area of ​​the spots.

7. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 6 is characterized in that: The method further comprises: Extracting spatial distribution characteristics and temporal change frequency characteristics from the cultivated land change patch data set, wherein the spatial distribution characteristics are quantified by the kernel density estimation method to determine the degree of patch aggregation, and the temporal change frequency characteristics are calculated based on the frequency and duration of patch occurrence; adjusting a geographic weight parameter of the change detection threshold according to the spatial distribution characteristics; The time window length parameter of the time series analysis model is optimized according to the time change frequency characteristics, the time window length is shortened in the changing area to capture the change events, and the time window length is extended in the stable area. 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 the fixed period.

8. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 7 is characterized in that: The adjusting the geographic weight parameter of the change detection threshold according to the spatial distribution characteristics includes: Dividing the target area into a plurality of geographic grid cells; Counting the degree of clustering of spots within each geographic grid unit, where the degree of clustering is calculated by combining the number of spots per unit area with the average spot area; Establishing a regression relationship model between the altitude and slope data of geographic grid cells and the density of change patches, wherein the regression relationship model uses a ridge regression algorithm to deal with multicollinearity problems; Predicting theoretical change densities under different geographical conditions based on the regression relationship model, wherein the theoretical change density reflects the natural constraint effect of topographic factors on cultivated land changes; 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 unit 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.

9. The method for intelligent identification of cultivated land change patterns based on remote sensing AI intelligent interpretation according to claim 6 is characterized in that: The spatial clustering process is performed on all candidate change regions, and a density peak clustering algorithm is used to merge adjacent pixels to form an initial change patch, including: A density peak clustering algorithm is used to perform spatial density analysis on the candidate change area, identify the spatial clustering area of ​​the candidate change points through the density peak clustering algorithm, and calculate the candidate point density of the clustering area to which each pixel belongs; Based on the eight-neighborhood connectivity criterion, the spatial clustering pixels identified by the density peak clustering algorithm are marked as connected domains. The density peak clustering algorithm divides the clustering area according to a preset density threshold. The connected domain marking recursively merges adjacent pixels through the seed filling algorithm. The spatial clustering pixels are pixels whose candidate point density is greater than the set density. Calculate the geometric center coordinates and minimum bounding rectangle of each marked connected domain, where the geometric center coordinates are determined based on the weighted average of all pixel coordinates in the connected domain, and the minimum bounding rectangle is obtained using a rotating caliper algorithm. Merging connected domains with overlapping bounding rectangles and recalculating geometric properties, wherein the merging operation is triggered according to the ratio of the intersection area of ​​the rectangles to the original area; The merged connected domains are screened by area, and connected domains with an area smaller than the minimum patch area threshold are eliminated. The minimum patch area threshold is dynamically adjusted according to the landscape fragmentation index calculated based on the ratio of the number of regional cultivated land patches to the total cultivated land area and the monitoring accuracy requirements, and the remaining connected domains are retained as the initial change patches.

10. An intelligent identification system for cultivated land change patterns based on remote sensing AI intelligent interpretation, characterized in that: It includes a processor and a memory, the memory is connected to the processor, the memory is used to store programs, instructions or codes, and the processor is used to execute the programs, instructions or codes in the memory to implement the intelligent identification method of cultivated land change patterns based on remote sensing AI intelligent interpretation as described in any one of claims 1 to 9.

Images