Multi-source high-resolution remote sensing image intelligent identification method for cultivated land change

Abstract

The application discloses a multi-source high-resolution remote sensing image cultivated land change intelligent identification method, relates to the technical field of remote sensing image data processing, and comprises the following steps: acquiring multi-temporal optical remote sensing images and synthetic aperture radar images and registering, constructing optical phenology sequences and radar phenology sequences pixel by pixel, calculating the cross-modal coupling degree between the two sequences, determining the cultivated land change region according to the coupling breakage degree exceeding the change criterion threshold, identifying the cultivated land change type according to the coupling breakage mode, inversely calculating the change occurrence time window based on the time inflection point of the coupling breakage degree, and feeding back the coupling breakage mode corresponding to the change type to calibrate the change criterion threshold. The application takes the breakage of the cross-modal coupling relationship instead of the difference amplitude in the mode as the change criterion, effectively suppresses the false changes caused by clouds, fog and crop rotation, and reduces the false positive rate of cultivated land change detection.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image data processing technology, specifically to a method for intelligent identification of farmland changes in multi-source high-resolution remote sensing images. Background Technology

[0002] Farmland protection is a core foundation of the national food security strategy. Refined, round-the-clock, and dynamic monitoring of farmland is a prerequisite for the timely detection of typical changes such as farmland abandonment and non-grain conversion. Traditional farmland supervision methods mainly rely on manual on-site inspections by grassroots natural resources departments and periodic cadastral map updates. These methods suffer from limited coverage, long monitoring cycles, high labor costs, and difficulty in timely detection of illegal changes, making it difficult to meet the timeliness requirements of large-scale dynamic farmland supervision.

[0003] With the development of remote sensing Earth observation technology, change detection methods based on remote sensing images have been widely applied to the dynamic monitoring of cultivated land. However, cultivated land is a special type of land cover with strong periodic apparent changes: the same cultivated land will experience multiple phenological stages such as sowing, emergence, greening, heading, maturity, and harvest within a single farming year. Its spectral response and radar backscattering in remote sensing images undergo dramatic and regular recurring changes with crop growth. The state of bare soil after harvest and bare land after construction occupation are easily confused in a single temporal image. At the same time, optical remote sensing images are affected by cloud cover, atmospheric conditions, and differences in the angle of sunlight, resulting in radiation differences between different temporal images that are unrelated to the actual changes in the land surface. These factors make change detection in cultivated land scenarios face a core technical problem that differs from general land cover change detection: how to distinguish the actual changes in cultivated land use from the aforementioned pseudo-changes under real multi-source remote sensing observation conditions where cloud cover, light differences, crop rotation, and phenological cycles coexist, thereby reducing the false alarm rate of cultivated land change detection.

[0004] Chinese invention patent application CN112541904A discloses an unsupervised remote sensing image change detection method, storage medium, and computing device. This method constructs a multi-scale graph convolutional neural network (MCN) to extract spatial and spectral features from two-temporal images and jointly calculate initial pseudo-labels. The concatenated two-temporal images are then input into the MCN for training to generate a two-channel difference map. A metric learning module is used to update the pseudo-labels and train the two-channel difference map. Finally, a binary change map is obtained by comparing the two channels. However, this scheme uses the difference map extracted between two-temporal images as the change criterion, essentially using the intensity of intramodal difference as a measure of change. In the context of cultivated land, the strongest differences between consecutive temporal phases come precisely from crop phenological cycles and bare soil after harvest. However, changes in actual land use may not manifest as a single maximum value of difference intensity. Therefore, this scheme struggles to distinguish between pseudo-changes caused by crop phenology and crop rotation and changes in the actual use of cultivated land, leading to numerous false alarms in cultivated land dynamic monitoring.

[0005] US Patent Application Publication No. US20250118054A1 discloses a method and system for remote sensing change detection using a segmented all-encompassing model. This method encodes previous and subsequent event images using the frozen features of the segmented all-encompassing model. After processing by a change modeler and a change cue, the segmented all-encompassing model decodes the data to obtain a change mask. However, this approach also uses the feature differences between the two time phases as the change criterion, focusing only on instantaneous difference detection between the two time phases. It does not utilize the temporal evolution of ground features and cannot characterize the periodic structure of farmland phenology. Therefore, it is equally difficult to suppress false changes in farmland scenarios where crop phenological cycles and light differences coexist.

[0006] Another type of method performs phenological decomposition on long-term remote sensing images of a single modality, fitting the intra-annual periodic trajectory of crops with a harmonic model, and identifying changes through residual or breakpoint detection in order to suppress spurious changes caused by seasonal differences. However, this type of method relies on densely and regularly sampled optical time series, and its phenological fitting accuracy drops sharply under real-world conditions where cloud pollution leads to sparse and irregular effective observations. More importantly, this type of method only uses time series of a single modality. When crop rotation occurs on an entire cultivated land, its phenological trajectory undergoes an overall phase shift, which is misjudged as a breakpoint by this type of method, thus still failing to distinguish between crop rotation and changes in the actual use of cultivated land.

[0007] In summary, existing methods for detecting farmland change, whether based on dual-temporal differences or single-modal time-series decomposition, rely on the magnitude of differences within the same modality as their change criteria. This makes it difficult to distinguish between true and false changes in farmland use under conditions of cloud cover, sunlight, crop rotation, and phenological cycles, resulting in a persistently high false alarm rate. Therefore, there is an urgent need for an intelligent method for identifying farmland change from multi-source high-resolution remote sensing images that can overcome the limitations of intra-modal difference criteria at the fundamental level. Summary of the Invention

[0008] To address the core technical bottleneck of existing technologies for detecting changes in cultivated land through remote sensing, which relies on intramodal difference amplitude as the change criterion and struggles to distinguish between true and false changes in cultivated land use under conditions of cloud cover, light differences, crop rotation, and phenological cycles, leading to excessively high false alarm rates, this invention provides an intelligent identification method for cultivated land changes from multi-source high-resolution remote sensing images. By treating the cross-modal phenological coupling degree between two phenological trajectories—optical remote sensing and synthetic aperture radar—as a physical invariant of the living, managed cultivated land, and using the degree of breakage of this coupling relationship rather than the difference amplitude of any single mode as the change criterion, this method achieves high-precision intelligent identification of changes in the true use of cultivated land from the perspective of cross-modal phenological coupling mechanisms without relying on densely regular sampling of optical time series.

[0009] The technical solution of this invention is as follows:

[0010] A method for intelligent identification of farmland changes from multi-source high-resolution remote sensing imagery includes the following steps:

[0011] Step S1: Acquire multi-temporal optical remote sensing images and multi-temporal synthetic aperture radar images covering the target area, register the optical remote sensing images and the synthetic aperture radar images to a unified coordinate system, and construct, pixel by pixel, an optical phenological sequence reflecting the change of canopy spectrum over time and a radar phenological sequence reflecting the change of backscattering over time.

[0012] Step S2: Calculate the cross-modal phenological coupling degree between the optical phenological sequence and the radar phenological sequence as phenological progresses, and obtain the cross-modal phenological coupling degree sequence. Pixels whose coupling breakage degree of the cross-modal phenological coupling degree sequence exceeds the change criterion threshold are identified as cultivated land change areas. The coupling breakage degree is the amount of attenuation of the cross-modal phenological coupling degree relative to its historical stable value.

[0013] Step S3: For the cultivated land change area, identify the cultivated land change type based on the coupling rupture mode presented by the cross-modal phenological coupling degree sequence. The cultivated land change type includes cultivated land converted to construction land, cultivated land converted to facility agricultural land, cultivated land converted to orchard and forest land, cultivated land abandoned and fallow, and cultivated land replanting.

[0014] Step S4: Determine the time window of change occurrence based on the time inflection point of the abrupt change in coupling breakdown degree in the cross-modal phenological coupling degree sequence, and output the farmland change area, the farmland change type and the change occurrence time window as farmland change identification results; feed back the coupling breakdown mode corresponding to the farmland change type to the change criterion threshold for calibration.

[0015] The beneficial effects of this invention are as follows:

[0016] First, by using the degree of breakage of the cross-modal phenological coupling between the optical phenological sequence and the radar phenological sequence as the criterion for change, rather than using the difference amplitude within any single mode, the spurious changes caused by clouds and illumination are effectively suppressed. The mechanism is as follows: the optical phenological trajectory reflects the canopy chlorophyll content and biomass, while radar backscattering reflects plant structure, water content, and surface roughness. Both are driven by the same agricultural management calendar and have a stable physical coupling. Cloud cover and illumination differences only disrupt the observation values ​​of the optical single mode, not the all-weather observation of the synthetic aperture radar, nor the existence of a coupling relationship between the two trajectories. Therefore, when using the degree of coupling breakage as the criterion, spurious changes in cloud cover will not cause coupling attenuation and will be naturally filtered out. Compared to using intra-modal difference maps as the criterion and including illumination and phenological differences in the change intensity, this scheme eliminates the influence of single-mode radiation disturbances on change detection from the source of the criterion.

[0017] Second, by using time-varying mutual information between the optical phenological sequence and the radar phenological sequence to characterize the cross-modal phenological coupling degree, zero false alarms for crop rotation across the entire cultivated land are achieved. The mechanism is as follows: crop rotation causes a synchronized overall phase shift between the optical and radar phenological trajectories, without disrupting the statistical dependency structure between them. The time-varying mutual information measures precisely this statistical dependency structure between the two trajectories, remaining unchanged by the synchronized shift of their rhythms and only sensitive to disruptions in the dependency structure itself. Therefore, although crop rotation changes the absolute phase of the phenological trajectory, it does not change the coupling degree and is not considered a change. Compared to single-modal phenological time series decomposition methods that misjudge the phase shift caused by crop rotation as a breakpoint, this scheme eliminates crop rotation false alarms by measuring the invariance of the measurement itself.

[0018] Third, by placing the calculation of the cross-modal phenological coupling degree on the phenological phase axis established by the all-weather regular sampling of the radar phenological sequence, and weighting the coupling degree according to the optical effective observation density, robust coupling degree estimation is achieved under the condition of sparse and irregular optical effective observation. The mechanism is as follows: Synthetic Aperture Radar provides regular all-weather sampling unaffected by clouds and fog. Establishing the phenological phase axis with it as a time reference can avoid the instability of phenological fitting caused by sparse optical observation. Furthermore, aligning the two sequences in the phenological phase domain rather than the calendar time domain eliminates spurious coupling breaks caused by the non-overlapping of the two modal sampling grids.

[0019] Fourth, by formalizing the coupling breakdown modes corresponding to different types of farmland change into discriminative constraints corresponding to physical mechanisms and feeding them back to the change criterion threshold for calibration, adaptive optimization of the change detection criterion is achieved while identifying change types, forming a closed-loop collaboration between detection and identification. The mechanism is as follows: different change types, such as farmland conversion to construction land, farmland abandonment, and farmland recultivation, physically correspond to different breakdown modes, such as permanent zeroing of coupling, slow decay, and re-increase from low to high. Feeding back the coupling breakdown degree corresponding to the identified breakdown modes to the change criterion threshold allows the threshold to adaptively tighten or loosen according to the dominant change type in the region, avoiding missed detections and false detections in areas where different change types are mixed.

[0020] Fifth, the synergistic effect of the aforementioned cross-modal coupling rupture criterion, time-varying mutual information measurement, phenological phase axis alignment, and rupture mode feedback exceeds the sum of their individual effects. Specifically, the cross-modal coupling rupture criterion eliminates false changes in cloud illumination; the time-varying mutual information measurement further eliminates false changes in crop rotation based on this criterion; phenological phase axis alignment provides a computational benchmark for the former two that still holds true under sparse and irregular sampling; and the rupture mode feedback uses the identification results to tighten the criterion in reverse, thus enhancing the accuracy of identifying real changes and the ability to suppress false changes. At the same time, the decay inflection point of the time-varying mutual information sequence, designed to suppress false alarms in crop rotation, also encodes the specific time of change occurrence, allowing the time window of change occurrence to be inverted without additional modules. This serves the enforcement and traceability of farmland abandonment and non-grain conversion, demonstrating the synergistic gain of a single technical feature solving multiple technical problems simultaneously. Attached Figure Description

[0021] Figure 1 This is an overall flowchart of the intelligent identification method for farmland changes in multi-source high-resolution remote sensing images according to the present invention;

[0022] Figure 2 This is a flowchart of step S1 of the present invention: multi-source image registration and phenological sequence construction;

[0023] Figure 3This is a flowchart of step S2 of the present invention, which involves calculating the cross-modal phenological coupling degree and determining the cultivated land change area.

[0024] Figure 4 This is a flowchart of step S3 of the present invention: identification of cultivated land change types and threshold feedback calibration;

[0025] Figure 5 This is a flowchart of the time window inversion of step S4 of the present invention. Detailed Implementation

[0026] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment takes a grain production functional area in southern China as the target area, with an area of ​​approximately 1200 km². The area is mainly composed of paddy fields and dry land, with some facility agriculture and scattered construction land. It exhibits various surface dynamics such as farmland conversion to construction land, farmland abandonment, and crop rotation, providing a typical scenario for verifying the farmland change identification capability of the present invention. The processing flow of this embodiment is implemented on a cloud-based remote sensing computing platform, and the entire process is as follows: Figure 1 As shown, it includes steps S1 to S4.

[0027] Step S1: Acquisition, registration, and phenological sequence construction of multi-source, multi-temporal images. For example... Figure 2 As shown, step S1 acquires multi-temporal optical remote sensing images and multi-temporal synthetic aperture radar (SAR) images covering the target area. The optical remote sensing images in this embodiment include images from domestic Gaofen series satellites, commercial high-resolution satellites, UAV aerial images, and Sentinel series satellites. The spatial resolution of the domestic Gaofen series satellite images is 2m to 8m, the spatial resolution of the commercial high-resolution satellite images is 0.5m, the spatial resolution of the UAV aerial images is 0.1m, and the spatial resolution of the Sentinel-2 optical satellite images is 10m. The SAR images use C-band images from the Sentinel-1 satellite, including both vertical transmission and vertical reception polarization modes, with a spatial resolution of 10m and a revisit period of 12 days. This embodiment acquires no fewer than 15 optical remote sensing images and no fewer than 30 SAR images of the target area within a single cropping year to cover the complete phenological cycle of crop growth.

[0028] After acquiring multi-source images, the optical remote sensing images and the synthetic aperture radar images are registered to a unified coordinate system. In this embodiment, the National Geodetic Coordinate System 2000 is used as the unified coordinate system. Images with different spatial resolutions are resampled and aligned, and all images are resampled to a unified 10m grid. Bilinear interpolation is used for resampling, and the root mean square error of registration is controlled within 0.5 grids to ensure the spatial consistency of multi-source images at the pixel scale.

[0029] Before constructing the optical phenological sequence, cloud detection and cloud masking are performed on the optical remote sensing images based on the quality assessment band, and atmospheric correction is applied. In this embodiment, the quality assessment band of Sentinel-2 imagery is used to identify clouds and cirrus clouds pixel by pixel in each optical image, removing cloud-contaminated pixels and retaining only clear-sky observations. Atmospheric correction is then performed on the retained valid observations to obtain the surface reflectance. Before constructing the radar phenological sequence, speckle noise suppression is achieved by applying local window Lee filtering to the synthetic aperture radar imagery. In this embodiment, a 3×3 local window Lee filter is used, and the filter intensity is dynamically adjusted through local mean and local variance to suppress speckle noise while preserving the edges and textures of ground features.

[0030] After the above preprocessing, an optical phenological sequence reflecting the change of canopy spectrum over time and a radar phenological sequence reflecting the change of backscattering over time are constructed pixel by pixel. In this embodiment, for each pixel, the optical phenological sequence is constructed using the normalized vegetation index over time, and the radar phenological sequence is constructed using the backscattering coefficient of vertical transmission and vertical reception polarization over time. To fully characterize the seasonal backscattering rhythm of crops, this embodiment further organizes synthetic aperture radar images within a growing year by quarter, calculates the backscattering ratios of the first and third quarters and the second and fourth quarters respectively, and uses the above-mentioned quarterly backscattering ratios as additional rhythmic components of the radar phenological sequence, so that the radar phenological sequence reflects instantaneous backscattering while explicitly characterizing the differences in the structural and water content rhythms of crops in different growing seasons. This embodiment found that the backscatter ratio of actively cultivated paddy fields and dry land exhibits a stable and regular rhythmic characteristic in different seasons. However, changes in land use such as construction occupation and abandonment can significantly alter the rhythmic structure of the backscatter ratio in that season. Therefore, the introduction of the quarterly backscatter ratio enhances the sensitivity of the radar phenological sequence to the dynamics of cultivated land.

[0031] In constructing the optical phenological sequence, in addition to the normalized vegetation index (NVI), this embodiment also extracts vegetation indices such as enhanced vegetation index and soil-adjusted vegetation index, water indices such as normalized water index and improved normalized water index, and building and bare soil indices such as normalized building index and bare soil index on a pixel-by-pixel basis. These are then integrated with elevation and slope information extracted from the digital elevation model to form a multi-dimensional optical phenological characterization. This embodiment utilizes a random forest classifier to systematically analyze the importance of each band and index, and selects core features by ranking them according to their feature contribution, thus reducing the computational burden caused by redundant features while ensuring the separability of land features. The core features in the above-mentioned multi-dimensional optical phenological characterization are used to assist in the identification of cultivated land change types in step S3, while the calculation of the cross-modal phenological coupling degree is mainly based on the optical phenological sequence composed of the normalized vegetation index to ensure the clarity of the physical meaning of the coupling degree calculation. Both the optical phenological sequence and the radar phenological sequence are indexed by calendar time. Due to the lack of effective optical observations during periods of cloud pollution, the optical phenological sequence is sparse and irregular on the time axis, while the radar phenological sequence is composed of all-weather observations and is regular and dense on the time axis. The optical phenological sequence and the radar phenological sequence are the inputs for the cross-modal phenological coupling degree calculation in step S2.

[0032] Step S2: Calculation of cross-modal phenological coupling degree and determination of cultivated land change areas. For example... Figure 3 As shown, step S2 calculates the cross-modal phenological coupling degree between the optical phenological sequence and the radar phenological sequence as phenological progresses. Considering that the optical phenological sequence is sparse and irregular while the radar phenological sequence is regular and dense, directly pairing and calculating the coupling relationship between the two at each calendar time would result in spurious coupling breaks unrelated to actual surface changes due to the non-overlapping sampling times of the two sequences. Therefore, this embodiment establishes a phenological phase axis based on the all-weather regular sampling of the radar phenological sequence and projects the effective observations of the optical phenological sequence onto the phenological phase axis.

[0033] First, the phenological phase of each pixel is determined. The phenological phase maps calendar time to a normalized position within the crop growth cycle, and its calculation formula is as follows:

[0034] ,in: For calendar time The corresponding phenological phase is a scalar, with a range of values ​​of . , dimensionless, is calculated by this formula and represents the normalized position of that moment in the crop growth cycle; The time is a calendar time, a scalar value in d, and is determined by the image imaging time. The starting time of crop phenology is a scalar quantity in d, determined by the inflection point of the backscattering recovery of the radar phenological sequence, corresponding to the crop greening period. Let be the length of the crop growth cycle, and be a scalar quantity with a range of values. The unit is d, which is determined by the cropping system of the dominant crop in the target area. In this example, it is taken as 365 days. The modulo operator represents... right Take the remainder. Then set the phenological phase axis. The effective observations of the optical phenological sequence are projected onto each phase grid point at equal intervals, based on the phenological phase. The projection uses phase-domain weighted interpolation, and the calculation formula is as follows:

[0035] ,in: For the optical phenological sequence in the 1st The projection value at each phase grid point is a scalar, and its value ranges from 1 to 2. , dimensionless, is calculated by this formula and characterizes the canopy spectral state at this phenological phase; For the first The phenological phase at each phase grid point is a scalar, with a value range of [value range missing]. , dimensionless, obtained by dividing the phase axis at equal intervals; The total number of phase grid points is a positive integer. In this embodiment, it is 48, which corresponds to dividing the crop growth cycle into 48 phase intervals. The phase grid point number has a value range of 1. to Integers; For the optical phenological sequence in the 1st One optically effective observation time The value of is a scalar, dimensionless, and is given by the optical phenological sequence constructed in step S1; The total number of effective optical observations is a positive integer, determined by the number of clear-sky observations retained after the cloud masking in step S1. The optically valid sequence number has a value range of [value range missing]. to Integers; For the first The first optically effective observation pair The projection weight of each phase grid point is a scalar with a value range of [value missing]. , dimensionless, is calculated using this formula; Phase grid points Phase with optically effective observation The annular phase distance between them is a scalar with a range of values ​​of 1. Dimensionless, defined as ,in For the minimum value operator, To be an absolute value operator, the definition of the annular distance makes the phase axes connect end to end to adapt to the cyclical characteristics of the phenological cycle; Let be the phase bandwidth, and be a scalar with a range of values. Dimensionless, in this embodiment, we take This controls the smoothness of the projection phase. If the value is too large, it will over-smooth and cover up phenological details. If the value is too small, it will create holes in sparse areas of optical observation. It is an exponential function with the natural constant as its base; Indicates to from to Summing each item individually.

[0036] During the projection process, the optically effective observation density projected onto the phenological phase axis is simultaneously counted as an indicator of cloud pollution level, and its calculation formula is as follows:

[0037] ,in: For the first The effective optical observation density at each phase grid point is a scalar quantity with a range of values ​​of . , dimensionless, is calculated by this formula, and characterizes the sufficiency of clear-sky optical observations at this phenological phase. The smaller the value, the more severe the cloud pollution and the sparser the optical observations at this phase. , , , , , The definition is the same as the formula mentioned above.

[0038] After phase domain alignment, the cross-modal phenological coupling degree between the optical phenological sequence and the radar phenological sequence is quantified. Considering that crop rotation occurs across the entire cultivated land, the two phenological trajectories undergo synchronous overall phase shift. If a coupling metric based on linear correlation is used, this synchronous shift might be misjudged as a coupling breakdown. Therefore, this embodiment uses time-varying mutual information between the optical phenological sequence and the radar phenological sequence to characterize the cross-modal phenological coupling degree. (The text then abruptly shifts to a different topic, mentioning a "first" and "second" phenological sequence, but without further context, a precise translation is impossible.) Centered on a phase grid point, with a width of [missing information]... Within a sliding phase window of 1 phase grid points, the joint distribution of the projected optical phenological values ​​and radar phenological values ​​is statistically analyzed, and their time-varying mutual information is calculated. The calculation formula is as follows: ,in: For the first The cross-modal phenological coupling degree at each phase grid point is a scalar with a value range of [value missing]. The unit is nat, which is calculated by this formula and represents the statistical dependence strength between the optical phenological sequence and the radar phenological sequence at this phenological phase. The larger the value, the tighter the coupling between the two phenological modes. The total number of bins for the value is a positive integer. In this embodiment, it is 8. The optical phenological values ​​and radar phenological values ​​are discretized into 8 bins with equal probability. This refers to the bin number for optical phenological values, with a value range of [value range missing]. to Integers; This refers to the bin number for radar phenological data, with a value range of [value range missing]. to Integers; In order to be in the first Within the sliding phase window centered on the first phase grid point, the optical phenological value falls within the first... The sub-box and the radar phenological value fall within the first... The joint probability of the bins is a scalar, and its value ranges from 1 to 2. , dimensionless, obtained from the frequency statistics of samples within the window; For optical phenological values ​​to fall into the first The marginal probability of each bin is a scalar, and its value ranges from 1 to 2. Dimensionless, by Calculated; For radar phenological values ​​to fall into the first The marginal probability of each bin is a scalar, and its value ranges from 1 to 2. Dimensionless, by The half-width of the sliding phase window is calculated. The value is a positive integer; in this embodiment, it is 4, which corresponds to a window covering 9 phase grid points. It is the natural logarithm function; Indicates to and Each from to The double summation. Since mutual information measures the statistical dependency structure between two variables, when the optical phenological trajectory and the radar phenological trajectory undergo synchronous overall phase shift, the paired samples within the window remain unchanged. Unchanged, therefore The synchronous translation remains unchanged, thus crop rotation does not cause a decrease in coupling.

[0039] As an alternative measure of the cross-modal phenological coupling degree, this embodiment can also measure the cross-modal phenological coupling degree based on the phase consistency between the key phenological periods and the backscattering rhythm phase. The key phenological periods are the phases of several characteristic phenological nodes representing the crop growth process extracted from the optical phenological sequence on the phenological phase axis, including the greening period, heading period, and maturity period. The phases of each key phenological period are located by the extreme and peak points of the growth rate of the optical phenological sequence on the phenological phase axis. The backscattering rhythm phase is the phase of the backscattering periodic rhythm characteristic points corresponding to the above-mentioned key phenological periods extracted from the radar phenological sequence on the phenological phase axis, located by the rising and falling inflection points of the backscattering of the radar phenological sequence as phenological progresses. On active cultivated land driven by the same agricultural management calendar, there is a stable phase correspondence between the key phenological periods determined by the crop canopy spectrum and the backscattering rhythm phase determined by plant structure and water content. Therefore, the phase consistency between the two can be used as a measure of the cross-modal phenological coupling degree. Specifically, the average annular phase distance between each phenological key period phase and its corresponding backscattering rhythm phase is denoted as the average annular phase distance. The cross-modal phenological coupling degree is 1 minus twice the average annular phase distance. The definition of the annular phase distance is the same as that in the aforementioned phase domain projection, and its value ranges from 0 to 0.5. Therefore, twice the average annular phase distance ranges from 0 to 1. The cross-modal phenological coupling degree is correspondingly between 0 and 1 and is a dimensionless quantity. The closer the value is to 1, the more consistent the phenological key period and the backscattering rhythm phase are, and the tighter the cross-modal phenological coupling is. When cloud cover or light differences cause only partial loss or anomalies in optically effective observations, the phases of each key phenological period recovered after phase-domain weighted projection remain essentially unchanged, while the backscattering rhythm phase is kept stable by all-weather radar observations. The phase consistency between the two does not decrease, and the cross-modal phenological coupling degree does not decay. This has the same physical effect as the time-varying mutual information-based metric in suppressing cloud cover and light spurious changes. When the actual use of cultivated land changes, the key phenological periods of crops disappear or migrate significantly, and the backscattering rhythm phase changes inconsistently with the change in land cover. The phase consistency between the two decreases significantly, and the cross-modal phenological coupling degree decays accordingly, triggering the determination of cultivated land change areas. In this embodiment, the aforementioned time-varying mutual information-based metric is preferentially used in areas with sufficient optically effective observations. In areas with sparse optically effective observations or high computational efficiency requirements, the aforementioned phase consistency-based metric can be used. Both can serve as implementation methods for the cross-modal phenological coupling degree.

[0040] The cross-modal phenological coupling degrees of each phase grid point are arranged in phase order and weighted according to the optical effective observation density to obtain the cross-modal phenological coupling degree sequence. The physical significance of weighting according to optical effective observation density here is that the coupling degree estimation for phase grid points with sparse optical observations has greater uncertainty. Assigning lower weights avoids unreliable estimations of overall pollution during cloud-polluted periods. Furthermore, the all-weather observations of synthetic aperture radar ensure that even for phase grid points where optical observations are completely absent, the radar phenological sequence still provides a continuous rhythmic reference, maintaining the temporal integrity of the cross-modal phenological coupling degree sequence. This design allows this scheme to robustly estimate coupling degrees in cloudy and rainy areas with low optical effective observation density, without relying on dense and regularly sampled optical time series as in single-modal phenological time series decomposition methods. Subsequently, the coupling breakage degree of the cross-modal phenological coupling degree sequence is calculated. The coupling breakage degree is the amount of decay of the cross-modal phenological coupling degree relative to its historical stable value, and its calculation formula is: ,in: For the first The coupling breakage degree at each phase grid point is a scalar, with a value range of [value missing]. , dimensionless, is calculated by this formula, and the closer the value is to 1, the more complete the coupling breakdown; The historical stable value of the cross-modal phenological coupling degree is a scalar with the unit nat. It is determined by the median of the cross-modal phenological coupling degree over multiple stable historical years and represents the coupling baseline of the pixel when there is no change in its use. The definition is the same as the aforementioned formula; Let be the pixel-level coupling breakage degree, and let be a scalar with a value range of . The dimensionless value is obtained by weighting the coupling breakage degree of each phase grid point with the effective optical observation density as the weight in this formula, so that the phase grid points with sufficient and high reliability of optical observation contribute more to the pixel-level determination. The definition is the same as the aforementioned formula; Indicates to from to Summation. When pixel-level coupling breakdown... Exceeding the change criterion threshold At that time, the pixel was identified as an area of ​​farmland change, among which Let be the threshold value for the change criterion, and let be a scalar value with a range of values. The value is dimensionless, and the initial value in this embodiment is 0.45. Feedback calibration is performed in step S3 based on the coupling breakdown mode corresponding to the type of cultivated land change. For pixels where only the optical phenological sequence deviates from its historical phenological trajectory, but the phase of the radar phenological sequence's backscattering rhythm does not shift accordingly, the cross-modal phenological coupling degree does not decrease. No more than Therefore, it is judged as a false change and excluded from the cultivated land change area. Step S3: Cultivated land change type identification and change criterion threshold feedback calibration. Figure 4 As shown, step S3 identifies the type of farmland change in the farmland change area determined in step S2 based on the coupling rupture pattern presented by the cross-modal phenological coupling degree sequence. This embodiment found that different types of farmland change physically correspond to different coupling breakdown modes: when farmland is converted to construction land, the surface is hardened and covered, crop phenology completely disappears, corresponding to the cross-modal phenological coupling degree permanently returning to zero and the backscattering of the radar phenological sequence increasing stepwise due to strong reflection from the hard surface; when farmland is abandoned, the surface is gradually covered by weeds, corresponding to the slow decay of the cross-modal phenological coupling degree and the presence of low amplitude and weak periodicity; when farmland is replanted, the surface recovers from a non-cultivated state to regular cultivation, corresponding to the cross-modal phenological coupling degree rising again from below the change criterion threshold; when farmland is converted to facility agriculture land, the surface is covered by greenhouse film, corresponding to the suppression of optical phenology and the appearance of regular mirror-like features in radar backscattering; when farmland is converted to orchards and woodlands, annual crops change to perennial vegetation, corresponding to the phenological cycle changing from a single peak within the year to a slow change over perennials. The above-mentioned types of rupture modes are physically distinct: construction and occupation cause the two phenological trajectories to disappear completely at the same time; abandonment causes the two trajectories to slowly degenerate into weak cycles; recultivation causes the two trajectories to strengthen again from weak; facility agricultural land causes the optics to be blocked by a thin film and the radar to appear as a mirror; and orchards and woodlands cause the cycle scale to shift from an intra-annual scale to a multi-year scale. Therefore, based on the coupled rupture modes, the change types can be finely distinguished while the change is detected. Different change types correspond to differentiated farmland supervision and disposal strategies, which can assist natural resource management departments in carrying out classified law enforcement and precise governance.

[0041] To improve the accuracy of change type identification by utilizing the aforementioned physical laws, this embodiment formalizes the coupling rupture mode into a discriminative constraint corresponding to the physical mechanism. First, rupture mode feature vectors are extracted from the cross-modal phenological coupling degree sequence and the radar phenological sequence, and their structure is as follows:

[0042] ,in: The rupture mode feature vector is a four-dimensional column vector, dimensionless, constructed by this formula, and characterizes the physical mode of coupling rupture. Let be the final value of the fracture, and be a scalar with a range of values. Dimensionless, the coupling breakage degree of the last segment of the cross-modal phenological coupling degree sequence is taken to characterize the degree of completeness of the breakage; Let be the fracture decay rate, and be a scalar quantity, in units of . The slope of the attenuation segment of the cross-modal phenological coupling degree sequence is fitted to characterize the rate of rupture. The backscattering step is a scalar quantity in dB, obtained from the difference in the mean backscattering before and after the change in the radar phenology sequence, and characterizes the degree of surface hardening. Let be the coupling recovery ratio, and be a scalar with a value range of . Dimensionless, obtained from the recovery ratio of the last segment of the cross-modal phenological coupling degree sequence relative to its valley value, characterizing the degree of recovery after replanting; superscript This represents the vector transpose. To eliminate the dimensional differences among the components, this embodiment... Each component is standardized to zero mean and unit variance before being used for discrimination. Each component is a dimensionless quantity.

[0043] The feature vectors of the rupture modes are then matched with physical templates for each type of farmland change using Mahalanobis distance constraints to identify the types of farmland change. The calculation formula is as follows:

[0044] ,in: To identify the types of changes in cultivated land, a category label is assigned, with a value that is one of the following: cultivated land converted to construction land, cultivated land converted to facility agricultural land, cultivated land converted to orchard or forest land, cultivated land abandoned or fallow, or cultivated land replanted. The label is dimensionless and is obtained by taking the category that minimizes the objective function in this formula. For candidate types of farmland change, iterate through the five types mentioned above; The definition is the same as the aforementioned formula; For the first The mean vector of the physical template of the rupture mode for the type of farmland change is a four-dimensional column vector, dimensionless, determined by the mean of the feature vector of the rupture mode of the known samples of this type, and reflects the typical physical rupture mode of this type. For the first The covariance matrix of the rupture mode for the type of farmland change is a 4×4 matrix, dimensionless, and is determined by the covariance of the eigenvectors of the rupture mode for known samples of this type. for The inverse matrix; For the first The physical prior probability of a type of farmland change is a scalar, and its value range is... , dimensionless, determined by the proportion of historical change types in the target area; It is the natural logarithm function; This indicates taking the candidate type that minimizes the subsequent expression. ; superscript This represents the vector transpose. By introducing Mahalanobis distance and physical prior discriminative constraints, the identification of change types is upgraded from purely data-driven to gray-box discrimination constrained by physical rupture modes, thus improving the ability to distinguish between types with similar physical meanings.

[0045] When identifying the types of farmland change, auxiliary discriminant features are also integrated to further improve accuracy. These auxiliary discriminant features include vegetation indices, water indices, and bare soil indices extracted from the optical remote sensing imagery, as well as elevation and slope extracted from the digital elevation model. The vegetation indices include the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EDI), and the Soil-Adjusted Vegetation Index (SDI). In this embodiment, a pre-trained classification model is used to integrate the fracture mode feature vector and the auxiliary discriminant features to output the types of farmland change. The classification model extracts spatial texture features from a convolutional neural network and then inputs them into a random forest classifier for high-dimensional feature classification. The convolutional neural network uses three layers of convolution and global average pooling, with an input image patch size of 9×9, and outputs a 64-dimensional spatial feature vector. To reduce the interference of farmland change area samples on the training of the classification model, the training samples are taken from long-term stable coupling regions where the cross-modal phenological coupling degree remains stable over multiple years. This concentrates the training samples in areas with stable land types and reliable labeling, improving the stability and reliability of the classification model.

[0046] After completing the change type identification, the prior coupling breakdown degree corresponding to the identified coupling breakdown mode is fed back to the change criterion threshold for adaptive calibration. The calculation formula is as follows:

[0047] ,in: The threshold value for the calibrated change criterion is denoted as , which is a scalar with a range of values ​​of . , dimensionless, is calculated by this formula and is used to determine the area of ​​farmland change in the next iteration; The threshold value for the change criterion before calibration is denoted as , which is a scalar and has a range of values ​​of . , dimensionless, is the threshold value for the change criterion in the current round, with an initial value of 0.45; The learning rate is the feedback rate, which is a scalar and has a range of values. The value is dimensionless; in this embodiment, it is set to 0.1 to control the magnitude of the threshold adjustment. A value that is too large will cause the threshold to oscillate, while a value that is too small will cause slow convergence. To identify type The corresponding prior coupling failure degree is a scalar, with a value range of [value range missing]. Dimensionless, determined by the typical coupling breakdown degree that this type should physically achieve. For example, the breakdown degree prior for converting farmland to construction land is close to 1, and the breakdown degree prior for abandoning farmland is a moderate value. The definition is the same as the aforementioned formula. This feedback adaptively tightens or relaxes the change criterion threshold according to the physical rupture characteristics of the dominant change type in the region, thereby forming a closed-loop synergy between change detection and change type identification.

[0048] Step S4: Inversion of the time window of change occurrence and output of results. For example... Figure 5As shown, step S4 determines the time window of change occurrence based on the time inflection point of the abrupt change in coupling breakdown degree in the cross-modal phenological coupling degree sequence. The time-varying mutual information constitutes a time-varying mutual information sequence along the phenological phase axis. Within the cultivated land change area, the time-varying mutual information sequence exhibits monotonically decaying and has an inflection point with an extremely high decay rate. This embodiment found that the time-varying mutual information sequence, originally designed to suppress false alarms about crop rotation, actually encodes the time of the actual change in cultivated land at its decay inflection point. Therefore, the time of change occurrence can be inverted from the coupling degree sequence without the need for an additional module.

[0049] First, calculate the decay rate of the time-varying mutual information sequence and locate the inflection point where the decay rate reaches its maximum. The calculation formula is as follows:

[0050] ,in: For the first The coupling breakage attenuation rate at each phase grid point is a scalar quantity in units of . The formula is obtained by central difference calculation and characterizes the rate of change of coupling breakage at this phase; and The first With the The coupling breakage degree at each phase grid point is defined as in the aforementioned formula; The phase interval between adjacent phase grid points is denoted by , which is a scalar and has a range of values ​​of . Dimensionless, in this embodiment, we take ; Let be the phase at the inflection point where the decay rate is extremely high, and let be a scalar with a value range of . , dimensionless, is obtained by taking the phase grid point that maximizes the decay rate from this formula, corresponding to the phenological phase in which the changes in cultivated land are most drastic; This indicates taking the phase grid point that maximizes the subsequent expression. The corresponding phase.

[0051] The inflection point phase is then inverted to a calendar time as the start time of the change window, and the phase corresponding to the coupling breakdown attenuation to a preset lower limit is inverted to the termination time. The calculation formula is as follows: ,in: Let be the starting time of the time window in which the change occurs, and be a scalar in d. The inflection point phase is obtained by inverting the calendar time using this formula. and The definition is the same as the aforementioned formula; The definition is the same as the aforementioned formula; The period number is a non-negative integer, determined by the number of phenological cycles spanned by the observation period, so that the inversion time falls within the actual observation period; Let be the termination phase of the time window in which the change occurs, and let be a scalar with a range of values. , dimensionless, is obtained by taking the phase grid point at which the coupling breakage degree first reaches the lower limit of complete breakage from this formula; Let be the complete fracture tolerance, and let be a scalar with a range of values. Dimensionless, in this embodiment we take 0.05. That is, the lower limit for determining whether the coupling breakage degree has reached complete breakage; This indicates retrieving the minimum element from the set; This indicates that the subsequent conditions are met. Termination phase. Press and The same method is used to reverse the process to the termination time.

[0052] Finally, the farmland change area, the farmland change type, and the change occurrence time window are output as farmland change identification results. This embodiment automatically marks short-term rapidly expanding areas, continuously coupled broken areas, and abnormal land occupation areas, and generates change risk levels based on the rate of change area growth, duration of change, and spatial expansion trend. This provides decision support for application scenarios such as the Ministry of Natural Resources' farmland protection satellite enforcement, dynamic monitoring of permanent basic farmland, monitoring of grain production functional zones by agricultural and rural departments, and evaluation of the effectiveness of high-standard farmland construction.

[0053] Steps S1 to S4 in this embodiment constitute a deeply coupled closed-loop collaborative process. The optical phenological sequence and the radar phenological sequence constructed in step S1 serve as inputs for the cross-modal phenological coupling degree calculation in step S2. The cross-modal phenological coupling degree sequence obtained in step S2 serves as both the basis for identifying farmland change types in step S3 and the basis for retrieving the change occurrence time window in step S4. The prior coupling breakdown mode corresponding to the farmland change type identified in step S3 is fed back to the change criterion threshold in step S2 for adaptive calibration, enabling subsequent rounds of change detection to dynamically optimize with the physical characteristics of the dominant change type. The input coupling between the preceding and following steps and the feedback calibration of the preceding parameters by the subsequent steps distinguish this scheme from the processing paradigm of existing change detection methods where detection and identification are mutually separated, achieving mutual enhancement of change detection accuracy and change type identification accuracy.

[0054] To verify the effectiveness of this invention, experiments were conducted in the aforementioned target area. The experimental data consisted of Sentinel-2 optical images and Sentinel-1 synthetic aperture radar images of the target area over a single farming year, supplemented by domestically produced high-resolution satellite imagery and UAV aerial images for ground truth verification. A total of 286 ground truth samples of farmland changes were manually verified and marked on-site, including 112 instances of farmland converted to construction land, 74 instances of farmland abandoned or fallow, 58 instances of farmland recultivated and replanted, 26 instances of farmland converted to facility agriculture land, and 16 instances of farmland converted to orchards or forest land.

[0055] In terms of change detection performance, the accuracy of this invention for detecting changes in cultivated land areas is 0.93, recall is 0.89, and F1 score is 0.91. As a control, a method using intramodal difference criteria has an F1 score of 0.74 on the same data and generates a large number of false alarms during the bare soil period after crop harvest and during cloud pollution periods, with a false alarm rate of 21%. This invention, by using cross-modal phenological coupling breakdown as the criterion, reduces the false alarm rate of the aforementioned spurious changes to 6%. In sample areas where crop rotation occurs across the entire cultivated land, the control method using linear correlation coupling metrics has a false alarm rate of 18%, while this invention, using time-varying mutual information metrics, reduces the false alarm rate for crop rotation to near zero.

[0056] In terms of change type recognition performance, after integrating physical discrimination constraints of rupture modes and auxiliary discrimination features, the overall accuracy of farmland change type recognition is 0.90, which is about 7 percentage points higher than the scheme that only uses auxiliary discrimination features without introducing physical constraints of rupture modes. Among them, the recognition accuracy of farmland converted to construction land is the highest, reaching 0.95. This is because the physical rupture mode corresponding to this type, which has a permanent zero cross-modal phenological coupling degree and a step increase in backscattering, is the most significant and has the greatest distinguishability from other types. Due to the transitional stage in phenological appearance, the recognition accuracy of farmland abandonment and farmland replanting is 0.88 and 0.86, respectively. However, after introducing the coupling recovery ratio as a rupture mode feature, the recognition accuracy of farmland replanting is improved by about 9 percentage points compared with the scheme without the coupling recovery ratio. In terms of change criterion threshold feedback calibration, after three rounds of iterative calibration, the overall F1 score of the target area is improved by about 4 percentage points compared with the fixed threshold scheme, indicating that the closed-loop collaboration between detection and recognition effectively improves the change detection accuracy. Regarding the inversion of the change occurrence time window, compared with the change occurrence period verified manually, the average error of the start time of the change occurrence time window inverted by the present invention is 11 days, which meets the accuracy requirements for locating the change period in the law enforcement traceability of abandoned farmland and non-grain farmland.

[0057] The above experimental results show that the present invention effectively solves the technical problem of excessively high false alarm rate in existing cultivated land change detection under the coexistence of cloud and fog, light and crop rotation and phenological cycle conditions by using cross-modal phenological coupling breakage as a change criterion, eliminating false alarms of crop rotation with time-varying mutual information, robustly estimating the coupling degree on the phenological phase axis and calibrating the change criterion threshold with breakage mode feedback. It achieves high-precision intelligent identification of changes in the true use of cultivated land.

[0058] It should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for intelligent identification of farmland changes from multi-source high-resolution remote sensing images, characterized in that, Includes the following steps: Acquire multi-temporal optical remote sensing images and multi-temporal synthetic aperture radar images covering the target area, register the optical remote sensing images and the synthetic aperture radar images to a unified coordinate system, and construct an optical phenological sequence reflecting the change of canopy spectrum over time and a radar phenological sequence reflecting the change of backscattering over time for each pixel. The cross-modal phenological coupling degree between the optical phenological sequence and the radar phenological sequence, which changes with the progression of phenology, is calculated to obtain a cross-modal phenological coupling degree sequence. The cross-modal phenological coupling degree characterizes the statistical dependence strength between the optical phenological sequence and the radar phenological sequence as phenology progresses. It is obtained by calculating the time-varying mutual information between the optical phenological sequence and the radar phenological sequence, or by calculating the phase consistency between the key phenological period extracted from the optical phenological sequence and the backscattering rhythm phase extracted from the radar phenological sequence. Pixels whose coupling breakage degree of the cross-modal phenological coupling degree sequence exceeds the change criterion threshold are identified as areas of farmland change. The coupling breakage degree is the amount of attenuation of the cross-modal phenological coupling degree relative to its historical stable value. For the areas of farmland change, the types of farmland change are identified based on the coupling breakdown patterns presented by the cross-modal phenological coupling degree sequence. The types of farmland change include farmland converted to construction land, farmland converted to facility agricultural land, farmland converted to orchards and forest land, farmland abandoned and fallow, and farmland replanting. The time window for the change is determined based on the time inflection point of the abrupt change in the coupling breakdown degree in the cross-modal phenological coupling degree sequence, and the farmland change area, the farmland change type, and the change occurrence time window are output as farmland change identification results; the coupling breakdown mode corresponding to the farmland change type is fed back to the change criterion threshold for calibration.

2. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 1, characterized in that, The calculation of the cross-modal phenological coupling degree between the optical phenological sequence and the radar phenological sequence as phenological progresses includes: extracting key phenological periods based on the optical phenological sequence, extracting backscattering rhythm phases based on the radar phenological sequence, and measuring the cross-modal phenological coupling degree by the phase consistency between the key phenological periods and the backscattering rhythm phases; for pixels where only the optical phenological sequence deviates from its historical phenological trajectory but the backscattering rhythm phase does not shift accordingly, they are determined to be pseudo-changes and excluded from the cultivated land change area.

3. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 2, characterized in that, The calculation of the cross-modal phenological coupling degree between the optical phenological sequence and the radar phenological sequence as phenological progresses further includes: when the effective observations of the optical remote sensing image are sparse due to cloud pollution and the sampling time does not coincide with the synthetic aperture radar image, establishing a phenological phase axis based on the all-weather regular sampling of the radar phenological sequence, projecting the effective observations of the optical phenological sequence onto the phenological phase axis, and then calculating the cross-modal phenological coupling degree; statistically analyzing the density of effective optical observations projected onto the phenological phase axis as an indicator of cloud pollution level, and weighting the cross-modal phenological coupling degree based on the density of effective optical observations.

4. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 3, characterized in that, When quantifying the cross-modal phenological coupling degree on the phenological phase axis, the time-varying mutual information between the optical phenological sequence and the radar phenological sequence is used to characterize the cross-modal phenological coupling degree. The time-varying mutual information remains unchanged in the overall translation of the rhythmic synchronization between the optical phenological sequence and the radar phenological sequence, and is only sensitive to the disruption of the statistical dependency structure between the optical phenological sequence and the radar phenological sequence, so that the rhythmic synchronization translation caused by crop rotation is not identified as the cultivated land change area.

5. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 4, characterized in that, The method of determining the change occurrence time window based on the time inflection point of the abrupt change in coupling breakdown degree in the cross-modal phenological coupling degree sequence includes: the time-varying mutual information forms a time-varying mutual information sequence along the phenological phase axis, and the time-varying mutual information sequence in the cultivated land change area exhibits monotonically decaying and has an inflection point with a very large decay rate; the phenological phase corresponding to the inflection point is inverted into a calendar time as the starting time of the change occurrence time window, and the phenological phase corresponding to the decay of the time-varying mutual information sequence to a preset lower limit is inverted into the ending time of the change occurrence time window.

6. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 5, characterized in that, The step of identifying farmland change types based on the coupling breakdown patterns presented by the cross-modal phenological coupling degree sequence includes: formalizing the coupling breakdown patterns into discriminant constraints corresponding to physical mechanisms: farmland conversion to construction land corresponds to the cross-modal phenological coupling degree permanently returning to zero and the backscattering of the radar phenological sequence increasing by a step; farmland abandonment corresponds to the cross-modal phenological coupling degree slowly decaying and retaining low amplitude and weak periodicity; farmland replanting corresponds to the cross-modal phenological coupling degree rising again from below the change criterion threshold; identifying the farmland change type based on the discriminant constraints, and feeding back the prior coupling breakdown degree corresponding to the identified coupling breakdown patterns to the change criterion threshold for adaptive calibration.

7. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 1, characterized in that, The multi-temporal optical remote sensing images and the multi-temporal synthetic aperture radar images include domestic Gaofen series satellite images, commercial high-resolution satellite images, UAV aerial images, and Sentinel series satellite images; the registration of the optical remote sensing images and the synthetic aperture radar images to a unified coordinate system includes resampling and aligning images with different spatial resolutions.

8. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 1, characterized in that, Before constructing the optical phenological sequence, cloud detection and cloud masking are performed on the optical remote sensing image based on the quality assessment band, and atmospheric correction is applied; before constructing the radar phenological sequence, speckle noise is suppressed on the synthetic aperture radar image using local window Lee filtering.

9. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 1, characterized in that, When identifying the types of farmland changes, auxiliary discriminant features are also integrated. These auxiliary discriminant features include vegetation indices, water indices, and bare soil indices extracted from the optical remote sensing images, as well as elevation and slope extracted from the digital elevation model. The vegetation indices include normalized vegetation index, enhanced vegetation index, and soil-adjusted vegetation index.

10. The intelligent identification method for farmland changes in multi-source high-resolution remote sensing imagery according to claim 1, characterized in that, The classification model is pre-trained to identify the types of farmland change. The training samples of the classification model are taken from long-term stable coupling regions where the cross-modal phenological coupling degree remains stable over multiple years, so as to reduce the interference of farmland change region samples on the training of the classification model.

Images