Land cover change detection method based on multi-temporal remote sensing images

Abstract

The application discloses a surface cover change detection method based on multi-temporal remote sensing images and relates to the technical field of natural resource law enforcement and ecological protection supervision; the method comprises the following steps: obtaining monthly multi-temporal remote sensing images covering 12-24 months, basic geographic information data and ecological function regionalization data and performing pretreatment; performing per-pixel time series decomposition on the standardized images, extracting change starting time, duration, change rate and stage spectral characteristics, and constructing change behavior fingerprints of each pixel; generating to-be-detected change map spots through spatial clustering; performing similarity matching on the change behavior fingerprints of the map spots and a pre-trained standard behavior fingerprint library, identifying surface cover change types; and finally, calculating the severity of illegal activities and evaluating ecological hazards, and generating a dynamic law enforcement report. The application breaks through the static thinking of traditional double-temporal change detection, can automatically identify multiple illegal activity types and evaluate ecological hazards.

Description

Technical Field

[0001] This invention relates to the field of natural resource law enforcement and ecological protection supervision technology, and in particular to a method for detecting land cover change based on multi-temporal remote sensing images. Background Technology

[0002] With the rapid development of remote sensing technology, satellite remote sensing imagery has become an important data source for natural resource law enforcement and ecological protection supervision. Currently, the mainstream methods for detecting land cover change are mainly based on dual-temporal remote sensing imagery. By comparing the differences in spectral characteristics between images at two different time points, the methods output binary results showing changes or no changes, or simple land cover conversion results.

[0003] However, existing technologies have the following significant drawbacks: First, relying solely on static spectral features cannot distinguish between land cover changes that are spectrally similar but different in nature. For example, legal agricultural rotation and illegal deforestation have highly similar spectral features in a single temporal image, leading to a very high misjudgment rate. Second, they cannot automatically identify the specific type of change. For different illegal activities such as illegal mining, deforestation, land reclamation from lakes, and illegal construction, a large amount of manual visual interpretation and on-site verification is required, resulting in low enforcement efficiency. Third, they cannot determine the legality and harmfulness of changes, nor can they assess the development stage and trend of illegal activities, leading to delayed enforcement, often only being discovered after serious ecological damage has been caused. Fourth, for gradual illegal activities such as illegal mining and deforestation, dual-temporal detection methods are unable to capture their slow change processes and cannot achieve early warning.

[0004] Although some change detection methods based on time-series remote sensing images have emerged in recent years, most of these methods only focus on whether changes have occurred or not, without modeling the unique temporal development process of different illegal activities. They cannot achieve intelligent identification of change types and automatic assessment of ecological hazards, and still cannot meet the actual needs of natural resource law enforcement and ecological protection supervision. Summary of the Invention

[0005] This invention overcomes the shortcomings of the prior art and provides a method for detecting land cover change based on multi-temporal remote sensing images.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a method for detecting land cover change based on multi-temporal remote sensing images, comprising the following steps:

[0007] S1. Acquire monthly multi-temporal remote sensing images, basic geographic information data, and ecological functional zoning data covering 12 to 24 months. Perform radiometric correction, geometric correction, and atmospheric correction preprocessing on the monthly multi-temporal remote sensing images to obtain a standardized multi-temporal remote sensing image dataset.

[0008] S2. Perform pixel-by-pixel temporal decomposition on the standardized multi-temporal remote sensing image dataset, extract the change start time, change duration, change rate and multiple stage spectral features corresponding to each pixel, and construct the change behavior fingerprint of each pixel based on the change start time, the change duration, the change rate and the multiple stage spectral features.

[0009] S3. Spatial clustering of the change behavior fingerprints of all pixels in the area to be detected to generate multiple change patches to be detected;

[0010] S4. Perform similarity matching between the change behavior fingerprint of the change patch to be detected and the standard behavior fingerprint database to identify the land cover change type corresponding to the change patch to be detected;

[0011] S5. Calculate the severity of illegal activities based on the duration of changes and spatial impact range of the detected change patches, assess the ecological hazards in conjunction with the ecological functional zoning data, and generate a dynamic law enforcement report that includes the start time, development stage, change trend, and hazard warning of the illegal activities.

[0012] In a preferred embodiment of the present invention, in S2, the step of performing pixel-by-pixel temporal decomposition on the standardized multi-temporal remote sensing image dataset includes:

[0013] S201. The STL time series decomposition algorithm is used to decompose the spectral time series of each pixel into a trend term, a seasonal term and a residual term. The abrupt change point is detected in the residual term to determine the start time of the change.

[0014] S202. Calculate the slope difference of the trend term before and after the abrupt change point to obtain the rate of change;

[0015] S203. Calculate the duration of the change by counting the time interval from the start time of the change to when the spectral characteristics tend to stabilize.

[0016] S204. Extract multiple stage spectral features during the change process according to a preset time interval. The stage spectral features include red band reflectance, near-infrared band reflectance, normalized vegetation index, and normalized water index.

[0017] In a preferred embodiment of the present invention, in S2, the step of constructing the change behavior fingerprint of each pixel includes:

[0018] S211. Convert the change start time into a relative timestamp relative to the start time of the monthly multi-temporal remote sensing image dataset;

[0019] S212. Normalize the duration of the change and the rate of change.

[0020] S213. Arrange the multiple stage spectral features in chronological order to form a spectral feature sequence;

[0021] S214. The relative timestamp, the normalized duration of change, the normalized rate of change, and the spectral feature sequence are concatenated into a one-dimensional feature vector, which serves as the change behavior fingerprint of the pixel.

[0022] In a preferred embodiment of the present invention, step S3, generating multiple change patches to be detected, includes:

[0023] S301. The DBSCAN density clustering algorithm is used to perform spatial clustering on the change behavior fingerprint of all pixels, and spatial distance threshold and feature similarity threshold are set.

[0024] S302. Merge adjacent pixels whose spatial distance is less than the spatial distance threshold and whose change behavior fingerprint similarity is greater than the feature similarity threshold into the same change patch to be detected;

[0025] S303. Remove noise patches with an area smaller than the preset minimum patch area threshold to obtain the final multiple change patches to be detected.

[0026] In a preferred embodiment of the present invention, the steps for constructing the standard behavioral fingerprint database include:

[0027] S401. Obtain a historical multi-temporal remote sensing image dataset containing multiple known land cover change types;

[0028] S402. Perform pixel-by-pixel temporal decomposition on the historical multi-temporal remote sensing image dataset and extract the standard change behavior fingerprint corresponding to each known change type.

[0029] S403. The standard variation behavior fingerprint is subjected to feature enhancement and cluster optimization using a self-supervised pre-training method of contrastive learning;

[0030] S404. Label the corresponding land cover change type for each cluster center to form a standard behavioral fingerprint database that includes types such as illegal mining, deforestation, land reclamation around lakes, and illegal construction.

[0031] In a preferred embodiment of the present invention, in step S4, the step of identifying the land cover change type corresponding to the change patch to be detected includes:

[0032] S411. Calculate the cosine similarity between the average change behavior fingerprint of the change patch to be detected and each standard change behavior fingerprint in the standard behavior fingerprint library;

[0033] S412. Select the land cover change type corresponding to the standard change behavior fingerprint with the highest cosine similarity and greater than the preset similarity threshold as the recognition result of the change patch to be detected.

[0034] S413. If all cosine similarities are less than the preset similarity threshold, then the detected change patch is marked as a suspected change patch.

[0035] In a preferred embodiment of the present invention, in step S5, the step of calculating the severity of the illegal activity includes:

[0036] S501. Obtain the spatial influence area and duration of change of the detected change patch;

[0037] S502. Based on the preset area weighting coefficient and duration weighting coefficient, calculate the severity index of illegal activities using a weighted summation method;

[0038] S503. According to the preset severity grading standard, the severity index of the illegal activity is divided into four levels: minor, moderate, relatively serious, and serious.

[0039] In a preferred embodiment of the present invention, step S5, which involves assessing the ecological hazards using the ecological functional zoning data, includes:

[0040] S511. Obtain the ecological protection red line range and ecological function zone level from the ecological function zoning data;

[0041] S512. Calculate the percentage of the overlapping area between the detected change patch and the ecological protection red line range;

[0042] S513. Based on the overlapping area ratio and the ecological functional zone level, the ecological hazard index is calculated using a graded assignment method.

[0043] S514. According to the preset hazard classification standard, the ecological hazard index is divided into four levels: low, medium, high and extremely high.

[0044] In a preferred embodiment of the present invention, the step of generating a dynamic law enforcement report in S5 includes:

[0045] S521. Extract the geographic location information, land cover change type, severity level of illegal activities, ecological hazard level, change start time, change duration and change trend of the detected change patch;

[0046] S522. Associate the information with the administrative division information in the basic geographic information data;

[0047] S523. Generate a dynamic law enforcement report that includes text descriptions, before-and-after image comparisons, time-series curves of the change process, and spatial distribution maps, and export it as a PDF file according to a preset format.

[0048] In a preferred embodiment of the present invention, the land cover change detection method further includes:

[0049] S6. Regularly acquire the latest monthly remote sensing images and continuously track and monitor the identified illegal activity patches;

[0050] S7. Update the change behavior fingerprint, change duration, and spatial impact range of the illegal activity patch;

[0051] S8. Recalculate the severity of illegal activities and their ecological harm;

[0052] S9. If illegal activities continue to develop or the degree of harm escalates, a hazard warning message will be automatically generated and sent to the law enforcement terminal.

[0053] This invention addresses the shortcomings of the prior art and has the following beneficial effects:

[0054] (1) By performing pixel-by-pixel temporal decomposition on multi-temporal remote sensing images, the starting time, duration, rate of change and stage spectral features of the change are extracted, and a feature vector that can characterize the entire process of land cover change is constructed. This breaks through the limitation of traditional dual-temporal change detection relying solely on static spectral features, solves the industry problem of difficulty in distinguishing changes with similar static spectral features but different behavioral patterns, and can accurately distinguish easily confused change types such as legal agricultural production and illegal deforestation and land reclamation.

[0055] (2) By pre-training through self-supervised learning, a standard behavioral fingerprint database containing multiple types of violations is constructed. By matching the change behavior fingerprint of the patch to be detected with the standard database, the type of violation can be automatically identified. This realizes an integrated process of land cover change detection, intelligent identification of change type and ecological hazard assessment, without the need for a large amount of manual visual interpretation, which significantly improves law enforcement efficiency.

[0056] (3) By using monthly multi-temporal remote sensing images covering 12 to 24 months, it is possible to capture the complete development process of progressive illegal activities such as illegal mining and deforestation, and to detect the signs of illegal activities several months in advance, thus achieving early warning.

[0057] (4) The present invention can calculate the severity of illegal activities based on the duration of change and the scope of impact, and assess the ecological hazards by combining ecological functional zoning data, and generate a dynamic law enforcement report that includes the start time, development stage, change trend and hazard warning of illegal activities, providing law enforcement personnel with a comprehensive and accurate basis for decision-making. Attached Figure Description

[0058] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0059] Figure 1 This is a flowchart of a preferred embodiment of the present invention; Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

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

[0062] The land cover change detection method based on multi-temporal remote sensing imagery overcomes the limitations of traditional two-temporal change detection, which relies solely on static spectral features, by constructing a temporal behavioral feature model of land cover change. It can automatically identify various types of illegal activities such as illegal mining, deforestation, land reclamation from lakes, and illegal construction, and assess their ecological hazards.

[0063] like Figure 1 As shown, the land cover change detection method based on multi-temporal remote sensing imagery includes the following steps:

[0064] S1. Acquire monthly multi-temporal remote sensing images, basic geographic information data, and ecological functional zoning data covering 12 to 24 months.

[0065] The monthly multi-temporal remote sensing images use optical satellite images with a spatial resolution of 10 to 30 meters, including but not limited to Sentinel-2 and Landsat-8 / 9 satellite images.

[0066] The basic geographic information data includes administrative division boundary data, road network data, water system data, and land use status data.

[0067] The ecological functional zoning data includes data on the scope of ecological protection red lines, data on the levels of ecological functional zones, and data on the distribution of ecologically sensitive areas.

[0068] The monthly multi-temporal remote sensing images were preprocessed sequentially using cloud detection, cloud removal, radiometric correction, geometric correction, and atmospheric correction. Cloud detection employed the SCL (Scene Classification Layer) cloud mask product built into the Sentinel-2 imagery to identify cloud-covered and cloud-shadowed pixels. Cloud removal used linear interpolation to fill missing values ​​in areas with cloud coverage less than 30% in a single image. For single images with cloud coverage greater than or equal to 30%, the weighted average of cloudless images from adjacent months was used for replacement, with the weights inversely proportional to the time interval. Radiometric correction converted the image's DN values ​​to atmospheric top reflectance using sensor calibration coefficients. Geometric correction used a quadratic polynomial transformation model to register the images to the WGS84 coordinate system, with registration errors controlled within 0.5 pixels. Atmospheric correction used a 6S radiative transfer model to remove the effects of atmospheric scattering and absorption, ultimately yielding a standardized multi-temporal remote sensing image dataset.

[0069] S2. Perform pixel-by-pixel temporal decomposition on the standardized multi-temporal remote sensing image dataset, extract the change start time, change duration, change rate and multiple stage spectral features corresponding to each pixel, and construct the change behavior fingerprint of each pixel based on the above features.

[0070] Specifically, the STL time series decomposition algorithm is used to decompose the spectral time series Y(t) of each pixel into a trend term T(t), a seasonal term S(t), and a residual term R(t), the mathematical expression of which is:

[0071] ;

[0072] Where t is the timestamp, ranging from 1 to N, and N is the total number of images in the standardized multi-temporal remote sensing image dataset. The seasonal cycle length of the STL time series decomposition is set to 12 months, corresponding to the interannual seasonal variation of the monthly images. Pettitt mutation point detection is performed on the residual term R(t), and the statistic U(t) is calculated:

[0073]

[0074] Where sgn(・) is the sign function, sgn(x)=1 when x>0, sgn(x)=0 when x=0, and sgn(x)=-1 when x<0. The time corresponding to the maximum value The significance level of the Pettitt mutation point detection is used as the start time of the change. Set to 0.05, when the maximum The corresponding p value is less than At that time, it is determined that the pixel has undergone a significant change. The rate of change v is obtained by calculating the difference in the slope of the trend term before and after the abrupt change point, and its mathematical expression is:

[0075] ;

[0076] in, The start time of the change The trend item for the previous three consecutive months The slope of the linear regression, The start time of the change The trend item continued for three consecutive months. The slope of the linear regression.

[0077] Statistics from the start time of the change The duration of the change is obtained from the time interval until the spectral characteristics tend to stabilize. The criterion for determining that the spectral characteristics tend to be stable is the residual terms within three consecutive time steps. The absolute values ​​are all less than the preset threshold. , A value of 0.01 is used, which can be adjusted within the range of 0.005-0.02 depending on the image noise level. If the change is still ongoing during detection (the spectral characteristics have not reached a stable state), the duration of the change, d, is calculated as the sum of the current detection time and the start time of the change. The difference.

[0078] According to the preset time interval Extract spectral features at multiple stages during the change process. The timeframe is set to 1 month, and the number of phases M is set to 6, meaning that spectral features are extracted for 6 consecutive months after the start of the change. The phased spectral features include red band reflectance, near-infrared band reflectance, normalized difference vegetation index (NDVI), and normalized difference water index (NDWI), wherein:

[0079] ;

[0080] ;

[0081] In the formula, For red band reflectivity, For green band reflectivity, This refers to the reflectivity in the near-infrared band.

[0082] Based on the features extracted above, a change behavior fingerprint is constructed for each pixel. First, the start time of the change is... Converted to the start time relative to the monthly multi-temporal remote sensing image dataset relative timestamp Then, the duration of the change was determined. and the rate of change The mathematical expression for min-max normalization is:

[0083] ;

[0084] ;

[0085] in, and These are the minimum and maximum durations of all pixel changes, respectively. and These are the minimum and maximum values ​​of the rate of change for all pixels, respectively.

[0086] The multiple stage spectral features are arranged in chronological order to form a spectral feature sequence. ,in The number of stage spectral features, each It is a four-dimensional vector containing red band reflectance, near-infrared band reflectance, NDVI, and NDWI.

[0087] Spectral feature sequence The unfolding method involves assembling the four-dimensional features of each time point in chronological order, and then using the relative timestamps. Duration of normalized changes Normalized rate of change The spectral feature sequence F is concatenated into a one-dimensional feature vector.

[0088] Then the relative timestamp Duration of normalized changes Normalized rate of change and the spectral feature sequence The features are concatenated into a one-dimensional feature vector, which serves as the fingerprint of the pixel's change behavior. The change behavior fingerprint FB of the pixel has a dimension of 1+1+1+6×4=27.

[0089] S3. Spatial clustering of the change behavior fingerprints of all pixels in the area to be detected to generate multiple change patches to be detected.

[0090] The DBSCAN density clustering algorithm is used to spatially cluster the change behavior fingerprints of all pixels, and a spatial distance threshold is set. and feature similarity threshold For any two adjacent pixels and If its Euclidean spatial distance is less than the spatial distance threshold Furthermore, the cosine similarity of its changing behavior fingerprint is greater than the feature similarity threshold. Then the pixel and They are grouped into the same cluster. The formula for calculating the cosine similarity is:

[0091] ;

[0092] in, The dot product of two fingerprint vectors representing changes in behavior. and These are the magnitudes of the two change behavior fingerprint vectors, respectively. All pixels within the same cluster are merged into a single change patch to be detected, and patches with areas smaller than a preset minimum patch area threshold are removed. Noise patches, The area is set to 500 square meters, but can be adjusted within the range of 300 to 1000 square meters according to law enforcement needs, resulting in multiple final change patches to be detected.

[0093] S4. Perform similarity matching between the change behavior fingerprint of the change patch to be detected and the standard behavior fingerprint database to identify the land cover change type corresponding to the change patch to be detected.

[0094] The standard behavioral fingerprint database is pre-built through the following steps:

[0095] Obtain a historical multi-temporal remote sensing image dataset containing various known land cover change types. The historical multi-temporal remote sensing image dataset includes manually labeled samples of change types such as illegal mining, deforestation, land reclamation from lakes, illegal construction, and legal agricultural production, with a total sample size of no less than 10,000.

[0096] The historical multi-temporal remote sensing image dataset is decomposed pixel by pixel to extract the standard change behavior fingerprint corresponding to each known change type.

[0097] A self-supervised pre-training method based on contrastive learning is used to enhance features and optimize clustering of the standard behavioral fingerprint. The feature extractor is trained using the InfoNCE loss function, the mathematical expression of which is:

[0098] ;

[0099] in, and These are two feature vectors obtained by applying different data augmentations to the same standard variation behavior fingerprint. For the feature vector of other standard variation behavior fingerprints, For temperature coefficient, The batch size was set to 100 training epochs. The AdamW optimizer was used, with an initial learning rate of 0.001 and a cosine annealing learning rate decay strategy. Specific data augmentation methods included time shifting (±1 month), Gaussian noise addition (standard deviation of 0.01), and feature scaling (scaling factor of 0.9~1.1).

[0100] The K-Means clustering algorithm was used to cluster the trained feature vectors. The number of clusters was determined by the elbow rule and was ultimately set to 5, corresponding to illegal mining, deforestation for cultivation, land reclamation from lakes, illegal construction, and legal agricultural production. Each cluster center was labeled with its corresponding land cover change type, forming a standard behavioral fingerprint database.

[0101] When performing similarity matching, the average value of the change behavior fingerprint of all pixels within the change patch to be detected is first calculated to obtain the average change behavior fingerprint of the change patch to be detected. .

[0102] Calculate the average change behavior fingerprint Each standard variation behavior fingerprint in the standard behavior fingerprint library Cosine similarity between .

[0103] Select the one with the highest cosine similarity that is greater than the preset similarity threshold. The type of land cover change corresponding to the standard change behavior fingerprint is used as the identification result of the change patch to be detected.

[0104] If all cosine similarities are less than the preset similarity threshold If the detected change patch is marked as a suspected change patch, it will be available for further manual verification.

[0105] S5. Calculate the severity of illegal activities based on the duration of changes and spatial impact range of the detected change patches, assess the ecological hazards in conjunction with the ecological functional zoning data, and generate a dynamic law enforcement report that includes the start time, development stage, change trend, and hazard warning of the illegal activities.

[0106] Specifically, the spatial influence area A and the duration d of the change in the detected change patch are obtained, and then a preset area weighting coefficient is applied. and duration weighting coefficient The severity index S of illegal activities is calculated using a weighted summation method:

[0107] ;

[0108] in, The preset maximum reference area, The preset maximum reference duration, The weighting is determined based on the experience in natural resource law enforcement practice that the area of ​​illegal activity has a greater impact on ecological damage than its duration.

[0109] According to the preset severity grading standard, the severity index S of the illegal activity is divided into minor. ,generally Heavy and serious Four levels.

[0110] Ecological hazards are assessed by combining the aforementioned ecological functional zoning data, and the ecological protection red line range and ecological functional zone level are obtained from the ecological functional zoning data. The ecological functional zone level is divided into three levels: general ecological functional zone, important ecological functional zone, and core ecological functional zone, corresponding to level coefficients of 1, 2, and 3, respectively. This level coefficient is determined based on the degree of ecological importance of different ecological functional zones in the "National Ecological Functional Zoning".

[0111] Calculate the percentage of the overlapping area between the detected change patch and the ecological protection red line. ,in The area of ​​overlap between the detected change patch and the ecological protection red line is denoted as .

[0112] According to the overlapping area ratio and the ecological functional zone level coefficient The ecological hazard index H was calculated using a graded assignment method:

[0113] ;

[0114] According to the preset hazard classification standard, the ecological hazard index H is divided into low hazard index H and high hazard index H. ,middle ,high and extremely high Four levels.

[0115] When generating a dynamic law enforcement report, the geographical location information, land cover change type, severity level of illegal activities, ecological hazard level, change start time, change duration, and change trend of the detected change patch are extracted.

[0116] The trend of change is determined by analyzing the trend items after the start time of the change. Linear regression prediction was performed to obtain the information. This information was then correlated with the administrative division information in the basic geographic information data to determine the province, city, county, and township-level administrative regions where the illegal activity occurred. Finally, a dynamic enforcement report was generated, including textual descriptions, before-and-after image comparisons, a time-series curve of the change process, and a spatial distribution map, and exported as a PDF file according to a preset format. The preset format PDF report includes the following core sections: basic information on the illegal activity, analysis of the change process, severity and harm assessment, enforcement recommendations, and attachments (image comparisons, time-series curves, and spatial distribution maps).

[0117] Furthermore, the method of this invention also includes a continuous tracking and monitoring step. The latest monthly remote sensing images are acquired periodically to continuously track and monitor identified illegal activity patches. The change behavior fingerprint, duration of change, and spatial impact range of the illegal activity patches are updated, and the severity and ecological harm of the illegal activity are recalculated. If any of the following conditions are met, it is determined that the illegal activity is continuing to develop or the degree of harm is escalating, and a hazard warning message is automatically generated and sent to the law enforcement terminal:

[0118] (1) The area of ​​illegal activity zones has expanded for three consecutive months;

[0119] (2) The severity level of the illegal activity has been increased by 1 level or more;

[0120] (3) The ecological hazard level is upgraded by 1 level or above.

[0121] Example

[0122] This embodiment takes a 100 km² ecological protection red line area in a certain county of a certain province as the detection object, and uses Sentinel-2 satellite imagery (spatial resolution of 10 meters) from January 2024 to December 2025 for 24 months as the data source. The basic geographic information data adopts the data from the 2023 National Geographic Census, and the ecological function zoning data adopts the results of the county's ecological protection red line delineation in 2022.

[0123] Step S1: Preprocess 24 Sentinel-2 images. After cloud removal, 22 valid images are obtained. The geometric correction registration error is 0.3 pixels. Finally, a standardized multi-temporal remote sensing image dataset is obtained, containing 10000×10000 pixels.

[0124] Step S2: Perform STL temporal decomposition (12-month seasonal cycle) on each pixel, and use Pettitt mutation point detection ( A total of 1256 pixels showed significant changes. Taking one pixel involved in illegal mining as an example, the start time of the change was detected. For June 2024, the rate of change The change lasted for d = 8 months (and was still ongoing as of February 2025). Six stage spectral features were extracted to construct a 27-dimensional change behavior fingerprint:

[0125] FB=[5,0.33,0.6,0.15,0.22,0.35,0.18,0.12,0.18,0.42,0.15,0.09,0.21,0.1 5,0.52,0.12,0.07,0.25,0.12,0.58,0.10,0.05,0.28,0.10,0.62,0.08,0.03].

[0126] Step S3: Perform spatial clustering using the DBSCAN algorithm ( =30 meters =0.8, =5), after removing noise patches with an area of ​​less than 500 square meters, a total of 32 change patches to be detected were generated, with the largest patch area being 12,500 square meters and the smallest patch area being 600 square meters.

[0127] Step S4: Perform similarity matching between the average change behavior fingerprint of the 32 change patches to be detected and the standard behavior fingerprint database. The study identified 3 illegal mining patches, 8 deforestation and land reclamation patches, 5 illegal construction patches, 14 legal agricultural production patches, and 2 suspected altered patches. The cosine similarity between the patches containing the illegal mining pixels and the standard illegal mining fingerprint was 0.92, thus identifying them as illegal mining types.

[0128] Step S5: Calculate the severity index of the illegal mining activity of this illegal mining plot: S = 0.7 × (12500 / 1000000) + 0.3 × (8 / 24) = 0.00875 + 0.1 = 0.10875, which is classified as minor; this plot is entirely located within the ecological protection red line area. =1.0), located in the core ecological functional area ( =3), the ecological hazard index H=3×(0.5×1.0+0.5×0.10875)=3×0.554375=1.663, which belongs to the extremely high level.

[0129] Generate a dynamic enforcement report, which includes the geographical location of the patch (a certain county, a certain township, a certain village), the start time of the change (June 2024), the development stage (active period), the change trend (the area will continue to expand in the next 3 months), the severity (minor), the ecological harm (extremely high), and the enforcement recommendation (immediate on-site verification and cessation).

[0130] Based on the preferred embodiments of the present invention described above, those skilled in the art can make various changes and modifications without departing from the inventive concept. The technical scope of this invention is not limited to the contents of the specification, but must be determined according to the scope of the claims.

Claims

1. A land cover change detection method based on multi-temporal remote sensing images, characterized in that, Includes the following steps: S1. Acquire monthly multi-temporal remote sensing images, basic geographic information data, and ecological functional zoning data covering 12 to 24 months. Perform radiometric correction, geometric correction, and atmospheric correction preprocessing on the monthly multi-temporal remote sensing images to obtain a standardized multi-temporal remote sensing image dataset. S2. Perform pixel-by-pixel temporal decomposition on the standardized multi-temporal remote sensing image dataset, extract the change start time, change duration, change rate and multiple stage spectral features corresponding to each pixel, and construct the change behavior fingerprint of each pixel based on the change start time, the change duration, the change rate and the multiple stage spectral features. S3. Spatial clustering of the change behavior fingerprints of all pixels in the area to be detected to generate multiple change patches to be detected; S4. Perform similarity matching between the change behavior fingerprint of the change patch to be detected and the standard behavior fingerprint database to identify the land cover change type corresponding to the change patch to be detected; S5. Calculate the severity of illegal activities based on the duration of changes and spatial impact range of the detected change patches, assess the ecological hazards in conjunction with the ecological functional zoning data, and generate a dynamic law enforcement report that includes the start time, development stage, change trend, and hazard warning of the illegal activities.

2. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S2, the step of performing pixel-by-pixel temporal decomposition on the standardized multi-temporal remote sensing image dataset includes: S201. The STL time series decomposition algorithm is used to decompose the spectral time series of each pixel into a trend term, a seasonal term and a residual term. The abrupt change point is detected in the residual term to determine the start time of the change. S202. Calculate the slope difference of the trend term before and after the abrupt change point to obtain the rate of change; S203. Calculate the duration of the change by counting the time interval from the start time of the change to when the spectral characteristics tend to stabilize. S204. Extract multiple stage spectral features during the change process according to a preset time interval. The stage spectral features include red band reflectance, near-infrared band reflectance, normalized vegetation index, and normalized water index.

3. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S2, the steps for constructing the change behavior fingerprint for each cell include: S211. Convert the change start time into a relative timestamp relative to the start time of the monthly multi-temporal remote sensing image dataset; S212. Normalize the duration of the change and the rate of change. S213. Arrange the multiple stage spectral features in chronological order to form a spectral feature sequence; S214. The relative timestamp, the normalized duration of change, the normalized rate of change, and the spectral feature sequence are concatenated into a one-dimensional feature vector, which serves as the change behavior fingerprint of the pixel.

4. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S3, the steps for generating multiple change patches to be detected include: S301. The DBSCAN density clustering algorithm is used to perform spatial clustering on the change behavior fingerprint of all pixels, and spatial distance threshold and feature similarity threshold are set. S302. Merge adjacent pixels whose spatial distance is less than the spatial distance threshold and whose change behavior fingerprint similarity is greater than the feature similarity threshold into the same change patch to be detected; S303. Remove noise patches with an area smaller than the preset minimum patch area threshold to obtain the final multiple change patches to be detected.

5. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: The steps for constructing the standard behavioral fingerprint database include: S401. Obtain a historical multi-temporal remote sensing image dataset containing multiple known land cover change types; S402. Perform pixel-by-pixel temporal decomposition on the historical multi-temporal remote sensing image dataset and extract the standard change behavior fingerprint corresponding to each known change type. S403. The standard variation behavior fingerprint is subjected to feature enhancement and cluster optimization using a self-supervised pre-training method of contrastive learning; S404. Label the corresponding land cover change type for each cluster center to form a standard behavioral fingerprint database that includes types such as illegal mining, deforestation, land reclamation around lakes, and illegal construction.

6. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S4, the step of identifying the land cover change type corresponding to the change patch to be detected includes: S411. Calculate the cosine similarity between the average change behavior fingerprint of the change patch to be detected and each standard change behavior fingerprint in the standard behavior fingerprint library; S412. Select the land cover change type corresponding to the standard change behavior fingerprint with the highest cosine similarity and greater than the preset similarity threshold as the recognition result of the change patch to be detected. S413. If all cosine similarities are less than the preset similarity threshold, then the detected change patch is marked as a suspected change patch.

7. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S5, the steps for calculating the severity of illegal activities include: S501. Obtain the spatial influence area and duration of change of the detected change patch; S502. Based on the preset area weighting coefficient and duration weighting coefficient, calculate the severity index of illegal activities using a weighted summation method; S503. According to the preset severity grading standard, the severity index of the illegal activity is divided into four levels: minor, moderate, relatively serious, and serious.

8. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S5, the steps for assessing ecological hazards based on the aforementioned ecological functional zoning data include: S511. Obtain the ecological protection red line range and ecological function zone level from the ecological function zoning data; S512. Calculate the percentage of the overlapping area between the detected change patch and the ecological protection red line range; S513. Based on the overlapping area ratio and the ecological functional zone level, the ecological hazard index is calculated using a graded assignment method. S514. According to the preset hazard classification standard, the ecological hazard index is divided into four levels: low, medium, high and extremely high.

9. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: In S5, the steps for generating dynamic law enforcement reports include: S521. Extract the geographic location information, land cover change type, severity level of illegal activities, ecological hazard level, change start time, change duration and change trend of the detected change patch; S522. Associate the information with the administrative division information in the basic geographic information data; S523. Generate a dynamic law enforcement report that includes text descriptions, before-and-after image comparisons, time-series curves of the change process, and spatial distribution maps.

10. The land cover change detection method based on multi-temporal remote sensing imagery according to claim 1, characterized in that: The land cover change detection method also includes: S6. Regularly acquire the latest monthly remote sensing images and continuously track and monitor the identified illegal activity patches; S7. Update the change behavior fingerprint, change duration, and spatial impact range of the illegal activity patch; S8. Recalculate the severity of illegal activities and their ecological harm; S9. If illegal activities continue to develop or the degree of harm escalates, a hazard warning message will be automatically generated and sent to the law enforcement terminal.

Images