A method for analyzing coupled evolution of land use and ecological environment

Abstract

The application provides a land use and ecological environment coupling evolution analysis method, relates to the geographic information processing field, and specifically comprises the following steps: acquiring multi-temporal remote sensing images and auxiliary data of a research area; pre-processing the multi-temporal remote sensing images to obtain a ground reflectivity image set; calculating at least one ecological index based on the pre-processed images to generate ecological index time series data, and performing interannual trend analysis on the ecological index; performing land use classification on the research area through a classification model based on the pre-processed images and the ecological index to obtain annual land use data; extracting land use change characteristics according to the annual land use data and calculating the coupling coordination degree between land use and ecological environment; and combining driving factor data and the coupling coordination degree to obtain a predicted land evolution result. The technical scheme of the application overcomes the problem in the prior art that long-time-series multi-source remote sensing data cannot be processed and land evolution is not accurately described.

Description

Technical Field

[0001] This invention relates to the field of geographic information processing technology, and specifically to a method for analyzing the coupled evolution of land use and ecological environment. Background Technology

[0002] Describing the spatiotemporal evolution patterns of land use and the ecological environment is beneficial for implementing protection measures for ecosystems.

[0003] Existing technologies typically rely on remote sensing imagery for single land use classifications or ecological index calculations, but they suffer from the following shortcomings: First, the quality and consistency of multi-source data are difficult to guarantee, with cloud cover and sensor attenuation causing radiometric shifts in time-series data; second, traditional local computing models are inefficient when processing long-term, large-scale data and are prone to memory overflows; third, there is a lack of full-chain coupling analysis of land use change and ecological environment response, making it particularly difficult to reveal the spatial heterogeneity of driving factors in different regions; and fourth, general classification systems are not optimized for the complex land types of coastal wetlands, resulting in low accuracy for categories such as wetlands and saline-alkali land.

[0004] Therefore, there is a need for a method that can efficiently process long-term multi-source remote sensing data and accurately characterize the land use and ecological environment coupling evolution analysis of coastal wetland land use evolution. Summary of the Invention

[0005] The main objective of this invention is to provide a method for coupled evolution analysis of land use and ecological environment, in order to solve the problem that existing technologies cannot process long-term multi-source remote sensing data and thus cannot accurately characterize land evolution.

[0006] To achieve the above objectives, this invention provides a method for analyzing the coupled evolution of land use and the ecological environment, specifically including the following steps: S1. Acquire multi-temporal remote sensing images of the study area and auxiliary data, including DEM data and driving factor data.

[0007] S2, preprocesses multi-temporal remote sensing images to obtain a set of surface reflectance images.

[0008] S3. Calculate at least one ecological index based on the preprocessed image, generate time series data of the ecological index, and perform interannual trend analysis on the ecological index.

[0009] S4, based on preprocessed images and ecological indices, classifies land use in the study area using a classification model to obtain annual land use data.

[0010] S5 extracts land use change characteristics based on annual land use data and calculates the coupling coordination degree between land use and the ecological environment.

[0011] S6 combines driving factor data and coupling coordination degree to identify the spatial differentiation characteristics of the influence of driving factors on coupling coordination degree, and obtain the predicted land evolution results.

[0012] Furthermore, the specific steps include the following: S2.1 Radiometric calibration is performed on multi-temporal remote sensing images to convert the original DN values ​​into apparent reflectance.

[0013] S2.2, The Sen2Cor algorithm is used to perform atmospheric correction on the radiometrically calibrated image to obtain the true surface reflectance.

[0014] S2.3, based on the image quality band, performs cloud, shadow and snow pixel masking on the atmospherically corrected image to retain the effective observed pixels.

[0015] S2.4, based on the boundary vector of the study area, the masked images are cropped and stitched together to form a set of surface reflectance images covering the study area.

[0016] Furthermore, the ecological indices include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Free Vegetation Capacity (FVC), Normalized Difference Water Index (NDWI), Improved Normalized Difference Water Index (MNDWI), Salinity Index (SI), Normalized Difference Building Index (NDSI), and Normalized Difference Building and Soil Index (NDBSI). The formulas for calculating the ecological indices are as follows: ; ; ; ; ; ; ; ; in, For near-infrared reflectivity, For red light band reflectivity, For blue light band reflectivity, For green light band reflectivity, For shortwave infrared reflectivity, For mid-infrared reflectivity, and These represent the minimum and maximum NDVI values ​​within the study area, respectively.

[0017] Furthermore, step S3 specifically includes the following steps: S3.1, the Mann-Kendall nonparametric test is used to determine the significance of the interannual variation trend of each ecological index at the pixel level.

[0018] S3.2, calculate the pixel-level change rate of each ecological index by combining Sen's Slope.

[0019] S3.3 generates a distribution map of the significance of ecological index trends and a distribution map of the rate of change.

[0020] Furthermore, step S4 specifically includes the following steps: S4.1 Construct a land use classification system applicable to the study area. The classification system includes: cultivated land, forest land, grassland, wetland, construction land, unused land, and water bodies.

[0021] S4.2 Select training samples for each type of land cover, with a sample size of not less than N and a uniform spatial distribution.

[0022] S4.3 uses the spectral bands and ecological indices of the preprocessed image as input features to train a random forest classification model.

[0023] S4.4 Optimize the parameters of the random forest classification model through cross-validation.

[0024] S4.5. The optimized model is applied to the preprocessed images year by year to obtain the land use classification raster map year by year.

[0025] Furthermore, step S5 specifically includes the following steps: S5.1 uses a zonal statistical method to calculate the mean, standard deviation, and coefficient of variation of ecological indicators for each land use type.

[0026] S5.2 uses Pearson correlation coefficient to perform correlation analysis, quantifying the association strength between land use type and ecological indicators. The calculation formula is as follows: ; in, The correlation coefficient, These represent land use type codes and ecological indicator values, respectively, where m is the sample size. for The average value, for The average value.

[0027] S5.3, the land use type distribution map and the ecological indicator spatial distribution map are rasterized to generate a composite layer of land use and ecological indicators. The expression is: ; in, For the coupled layer pixel values, For land use type pixel values, For the first The pixel value of each ecological indicator This is a spatial superposition operation.

[0028] S5.4, using Getis-Ord Hotspot or coldspot analysis is conducted to identify high-value and low-value clusters of ecological indicators. The calculation formula is as follows: ; in, This is a spatial clustering statistic. This is the spatial weight matrix. For pixel index values, S is the mean of the indicator, and S is the standard deviation of the indicator.

[0029] Furthermore, step S5 also includes the following steps: S5.5, Greenness NDVI, Humidity Wet, Dryness NDBSI, and Heat LST are selected as evaluation factors. Each evaluation factor is normalized to the [0,1] interval. Positive factors are: ; Negative factors are adopted as follows: ; Where NI is the standardized value, I is the original value, and I max I min These are the maximum and minimum values ​​of the factor, respectively.

[0030] S5.6, perform principal component analysis on the standardized factors, using the first principal component PC1 as the initial value of RSEI. The calculation formula is as follows: ; in, , The first principal component loading coefficient is denoted as .

[0031] Normalizing PC1 yields the final RSEI index, calculated using the following formula: ; The RSEI ranges from 0 to 1, with a higher value indicating better regional ecological quality.

[0032] S5.7 performs overlay analysis on the annual land use classification map, calculates the land use transfer matrix, and clarifies the transfer-in area, transfer-out area, and conversion direction of various land features; based on the transfer matrix, it calculates the land use dynamic degree and expansion intensity index to identify hotspot areas of land use change.

[0033] S5.8, construct evaluation indicators for the land use system and the ecological environment system. The land use system evaluation indicators include the area proportion of various land features, land use dynamics, and the area of ​​key transformation types. The ecological environment system evaluation indicators include the mean vegetation index, mean salinity index, mean water body index, and remote sensing ecological index (RSEI). Range standardization is performed on each evaluation indicator. Positive indicators adopt: ;

[0034] Negative indicators are adopted as follows: ; in, The standardized values ​​of the evaluation indicators The original values ​​of the evaluation indicators, For the first The global maximum value of the evaluation index For the first The global minimum value of each evaluation indicator.

[0035] The weights, information entropy, and difference coefficients of the evaluation indicators are calculated using the entropy weight method to obtain the weights of each evaluation indicator. ; ; ; ; in, For the first Evaluation indicator number The weight of each evaluation unit For the first Information entropy of the evaluation indicators For the first The coefficient of difference of the evaluation indicators For the first The entropy weight of each evaluation indicator The total number of evaluation units, The total number of evaluation indicators, It is the natural logarithm function.

[0036] Furthermore, step S5 also includes the following steps: S5.9, calculate the comprehensive evaluation value of the land use system respectively. Comprehensive evaluation value of ecological environment system The formula is: ; ; in, This is the comprehensive evaluation value of the land use system. This is a comprehensive evaluation value for the ecological environment system. For the land use system The weight of each evaluation indicator For the ecological environment system The weight of each evaluation indicator For the land use system Standardized values ​​of the evaluation indicators For the ecological environment system Standardized values ​​of the evaluation indicators This represents the total number of indicators in the land use system. This represents the total number of indicators for the ecological and environmental system.

[0037] S5.10, Calculate the coupling degree C between the land use system and the ecological environment system according to the coupling degree formula: ; Where C is the coupling degree, which takes the value [0,1]. The larger the value, the higher the coupling degree.

[0038] S5.11, Calculate the overall development level T: ; Where T represents the comprehensive development level of land use and ecological environment. , These are the weighting coefficients;

[0039] S5.12, Calculate the coupling coordination degree D, and classify the coupling coordination level according to the D value: ; Where D is the coupling coordination degree, with a value of [0,1].

[0040] Furthermore, step S6 specifically includes the following steps: S6.1 Quantitative analysis of regional landscape spatial pattern is conducted by selecting patch area, patch shape index, fractal dimension, patch density, edge density, and aggregation index.

[0041] S6.2 uses Getis-Ord Spatial statistics are used to detect hotspots or cold spots, identifying high-value clusters, low-value clusters, and insignificantly stable areas of land use or ecological environment change. ; in, This is a spatial clustering statistic for coupling coordination degree. For the coupling coordination degree of pixel q, The global mean of the coupling coordination degree. This represents the global standard deviation of the coupling coordination degree.

[0042] like >0 indicates that pixel p is a high-coordination hotspot area; <0 indicates that pixel p is a low-coordination cold spot area. , indicating that pixel p is a non-significantly stable region.

[0043] S6.3 Select driving factors, standardize the driving factors and perform multicollinearity tests to eliminate redundant factors.

[0044] S6.4. Using the driving factors selected in step S6.3 as input, random forest, XGBoost and multiple linear regression algorithms are used to construct driving models in high coordination hot spots, low coordination cold spots and non-significant stable areas respectively. The root mean square error, coefficient of determination and mean absolute error are used to compare and evaluate the models. The optimal model is used to calculate the predicted value of coupling coordination degree.

[0045] Section S6.5 introduces the SHAP value interpretation method to perform interpretability analysis on the optimal model, quantifying the contribution, direction, and nonlinear relationship of each driving factor. The formula is: ; in, This is the predicted value for coupling coordination. The average predicted value for all samples is given by m, where m is the total number of driving factors. is the SHAP value of the j-th driving factor. A positive value indicates a positive contribution, and a negative value indicates a negative contribution. The larger the absolute value, the stronger the contribution.

[0046] Furthermore, the driving factor data in step S6 include: temperature, precipitation, altitude, slope, aspect, soil type, and soil organic matter.

[0047] The present invention has the following beneficial effects: (1) By adopting a specific combination of “GEE cloud platform + Sen2Cor algorithm + QA60 quality band mask”, more than 1,200 Sentinel-2 images were processed by cloud parallel computing, the calculation cycle was shortened by more than 70%, and the effective pixel ratio was ≥85%. This solved the contradiction between efficiency and accuracy in long time series data processing. It was not a simple superposition, but produced an unexpected technical effect of “high efficiency + high fidelity”.

[0048] (2) A full-chain analysis system of “classification change → single factor response → coupling coordination → spatial driving” was constructed. For the first time in coastal wetland research, the land use transfer matrix, ecological index response, coupling coordination degree and geographical weighted regression were organically combined. This broke through the limitations of the traditional single analysis model and could clearly identify the differences in ecological response and dominant driving factors in different regions and different transformation types, filling the technical gap.

[0049] (3) By integrating the coupling coordination degree with the GWR model and using the coupling coordination degree as the dependent variable for spatial regression, the spatial heterogeneous influence of driving factors such as climate and population on the coordination state of the "land use-ecological environment" system was revealed, which deepened the understanding of the coupling mechanism of human-land relationship.

[0050] (4) A seven-category exclusive classification system including "coastal wetlands / inland wetlands" was constructed for a certain region. Combined with high-resolution image samples and random forest optimization, the wetland classification accuracy was improved to more than 88%, which solved the problem of insufficient accuracy of the general classification system in complex wetland areas. Attached Figure Description

[0051] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings: Figure 1 A flowchart of a method for analyzing the coupled evolution of land use and ecological environment according to Embodiment 1 of the present invention is shown.

[0052] Figure 2 The diagram shows the structure of a system for analyzing the coupled evolution of land use and ecological environment according to Embodiment 2 of the present invention.

[0053] Figure 3 A heat map of land use transfer matrix analysis according to Embodiment 1 of the present invention is shown.

[0054] Figure 4 A stacked diagram showing the land use structure proportion statistics of Embodiment 1 of the present invention is shown.

[0055] Figure 5 The results show the comparison of land use identification accuracy in a certain area using the method provided by the present invention. Detailed Implementation

[0056] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. 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.

[0057] Example 1 like Figure 1 The method for analyzing the coupled evolution of land use and ecological environment shown includes the following steps: S1. Acquire multi-temporal remote sensing images and auxiliary data of the study area. The auxiliary data includes DEM data and driving factor data. The multi-temporal remote sensing images are Sentinel-2 images of each growing season from 2016 to 2024.

[0058] S2 uses the Google Earth Engine cloud platform to preprocess multi-temporal remote sensing images to obtain a set of surface reflectance images.

[0059] S3. Calculate at least one ecological index based on the preprocessed image, generate time series data of the ecological index, and perform interannual trend analysis on the ecological index.

[0060] S4, based on preprocessed images and ecological indices, classifies land use in the study area using a classification model to obtain annual land use data.

[0061] S5 extracts land use change characteristics based on annual land use data and calculates the coupling coordination degree between land use and the ecological environment.

[0062] S6 combines driving factor data and coupling coordination degree to identify the spatial differentiation characteristics of the influence of driving factors on coupling coordination degree, and obtain the predicted land evolution results.

[0063] Specifically, it includes the following steps: S2.1 Radiometric calibration is performed on multi-temporal remote sensing images to convert the original DN values ​​into apparent reflectance.

[0064] S2.2, The Sen2Cor algorithm is used to perform atmospheric correction on the radiometrically calibrated image to obtain the true surface reflectance.

[0065] S2.3, based on the image quality band, performs cloud, shadow and snow pixel masking on the atmospherically corrected image, retaining the effective observed pixels, the image quality band is the QA60 band of Sentinel-2.

[0066] S2.4, based on the boundary vector of the study area, the masked images are cropped and stitched together to form a set of surface reflectance images covering the study area.

[0067] Specifically, the ecological indices include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Free Vegetation Capacity (FVC), Normalized Difference Water Index (NDWI), Improved Normalized Difference Water Index (MNDWI), Salinity Index (SI), Normalized Difference Building Index (NDSI), and Normalized Difference Building and Soil Index (NDBSI). The formulas for calculating the ecological indices are as follows: ; ; ; ; ; ; ; ; in, For near-infrared reflectivity, For red light band reflectivity, For blue light band reflectivity, For green light band reflectivity, For shortwave infrared reflectivity, For mid-infrared reflectivity, and These represent the minimum and maximum NDVI values ​​within the study area, respectively.

[0068] Specifically, step S3 includes the following steps: S3.1, the Mann-Kendall nonparametric test is used to determine the significance of the interannual variation trend of each ecological index at the pixel level.

[0069] S3.2, calculate the pixel-level change rate of each ecological index by combining Sen's Slope.

[0070] S3.3 generates a distribution map of the significance of ecological index trends and a distribution map of the rate of change.

[0071] Specifically, step S4 includes the following steps: S4.1 Construct a land use classification system applicable to the study area. The classification system includes: cultivated land, forest land, grassland, wetland, construction land, unused land, and water bodies.

[0072] S4.2 Select training samples for each type of land cover, with a sample size of not less than N and a uniform spatial distribution; in this embodiment, N is 500; the selection of training samples combines high spatial resolution imagery and field survey data.

[0073] S4.3 uses the spectral bands and ecological indices of the preprocessed image as input features to train a random forest classification model.

[0074] S4.4 Optimize the parameters of the random forest classification model through cross-validation; ensure that the overall classification accuracy of the random forest model is not less than 85% and the Kappa coefficient is not less than 0.8.

[0075] S4.5. The optimized model is applied to the preprocessed images year by year to obtain the land use classification raster map year by year.

[0076] Specifically, step S5 includes the following steps: Using the annual land use type map of the study area as the spatial basis, key ecological indicators such as vegetation cover, water index, soil moisture, and biomass were selected to carry out statistical coupling analysis, calculate the mean, standard deviation, and coefficient of variation, and conduct correlation analysis to quantify the correlation strength between land use type and ecological indicators.

[0077] S5.1 uses a regional statistical method to calculate the mean, standard deviation, and coefficient of variation of ecological indicators for each land use type. The calculation formulas include: ; ; ; in, Let be the mean value of a certain ecological indicator under the i-th land use category. Standard deviation The coefficient of variation is 1. Let be the index value of the j-th pixel within the i-th land use category, and n be the total number of pixels.

[0078] S5.2 uses Pearson correlation coefficient to perform correlation analysis, quantifying the association strength between land use type and ecological indicators. The calculation formula is as follows: ; in, The correlation coefficient, These represent land use type codes and ecological indicator values, respectively, where m is the sample size. for The average value, for The average value.

[0079] Using a spatial explicit approach, a composite layer of "land use-ecological indicators" is generated by overlaying GIS raster data. Hotspot or cold spot analysis is then performed to identify areas where ecological indicators are clustered with high or low values, thus depicting the spatial pattern of land use ecological effects.

[0080] S5.3, the land use type distribution map and the ecological indicator spatial distribution map are rasterized to generate a composite layer of land use and ecological indicators. The expression is: ; in, For the coupled layer pixel values, For land use type pixel values, For the first The pixel value of each ecological indicator This is a spatial superposition operation.

[0081] S5.4, using Getis-Ord Hotspot or coldspot analysis is conducted to identify high-value and low-value clusters of ecological indicators. The calculation formula is as follows: ; in, This is a spatial clustering statistic. This is the spatial weight matrix. For pixel index values, S is the mean of the indicator, and S is the standard deviation of the indicator.

[0082] Specifically, step S5 also includes the following steps: Greenness index (NDVI), humidity index (Wet), dryness index (NDBSI), and heat index (LST) were selected as evaluation factors. After data standardization, principal component analysis, and normalization, a remote sensing ecological index (RSEI) was constructed to achieve quantitative and spatially visualized evaluation of regional ecological quality.

[0083] S5.5, Greenness NDVI, Humidity Wet, Dryness NDBSI, and Heat LST are selected as evaluation factors. Each evaluation factor is normalized to the [0,1] interval. Positive factors are: ; Negative factors are adopted as follows: ; Where NI is the standardized value, I is the original value, and I max I min These are the maximum and minimum values ​​of the factor, respectively.

[0084] S5.6, perform principal component analysis on the standardized factors, using the first principal component PC1 as the initial value of RSEI. The calculation formula is as follows: ; in, , The first principal component loading coefficient is denoted as .

[0085] Normalizing PC1 yields the final RSEI index, calculated using the following formula: ; The RSEI ranges from 0 to 1, with a higher value indicating better regional ecological quality.

[0086] S5.7 extracts land use change characteristics based on annual land use data, and clarifies the patterns of land use transformation and hotspot areas.

[0087] Overlay analysis of annual land use classification maps is performed to calculate the land use transfer matrix, clarifying the transfer-in area, transfer-out area, and conversion direction of various land features; based on the transfer matrix, land use dynamics and expansion intensity indices are calculated to identify hotspots of land use change.

[0088] S5.8, construct evaluation indicators for the land use system and the ecological environment system. The land use system evaluation indicators include the area proportion of various land features, land use dynamics, and the area of ​​key transformation types. The ecological environment system evaluation indicators include the mean vegetation index, mean salinity index, mean water body index, and remote sensing ecological index (RSEI). Range standardization is performed on each evaluation indicator. Positive indicators adopt: ; Negative indicators are adopted as follows: ; in, The standardized values ​​of the evaluation indicators The original values ​​of the evaluation indicators, For the first The global maximum value of the evaluation index For the first The global minimum value of each evaluation indicator.

[0089] The weights, information entropy, and difference coefficients of the evaluation indicators are calculated using the entropy weight method to obtain the weights of each evaluation indicator. ; ; ; ; in, For the first Evaluation indicator number The weight of each evaluation unit For the first Information entropy of the evaluation indicators For the first The coefficient of difference of the evaluation indicators For the first The entropy weight of each evaluation indicator The total number of evaluation units, The total number of evaluation indicators, It is the natural logarithm function.

[0090] Specifically, step S5 also includes the following steps: S5.9, calculate the comprehensive evaluation value of the land use system respectively. Comprehensive evaluation value of ecological environment system The formula is: ; ; in, This is the comprehensive evaluation value of the land use system. This is a comprehensive evaluation value for the ecological environment system. For the land use system The weight of each evaluation indicator For the ecological environment system The weight of each evaluation indicator For the land use system Standardized values ​​of the evaluation indicators For the ecological environment system Standardized values ​​of the evaluation indicators This represents the total number of indicators in the land use system. This represents the total number of indicators for the ecological and environmental system.

[0091] S5.10, Calculate the coupling degree C between the land use system and the ecological environment system according to the coupling degree formula: ; Where C is the coupling degree, which takes the value [0,1]. The larger the value, the higher the coupling degree.

[0092] S5.11, Calculate the overall development level T: ; Where T represents the comprehensive development level of land use and ecological environment. , The weighting coefficients satisfy α+β=1, and in this invention, α=β=0.5 are used.

[0093] S5.12, Calculate the coupling coordination degree D, and classify the coupling coordination level according to the D value: ; Where D represents the coupling coordination degree, taking values ​​[0,1], with larger values ​​indicating higher coordination. Coupling coordination levels are classified based on the D value: 0≤D<0.4 indicates imbalance and decline, 0.4≤D<0.6 indicates transition and harmony, and 0.6≤D≤1 indicates improved coordination.

[0094] Specifically, step S6 includes the following steps: S6.1 Quantitative analysis of regional landscape spatial pattern is conducted by selecting patch area, patch shape index, fractal dimension, patch density, edge density, and aggregation index.

[0095] S6.2 uses Getis-Ord Spatial statistics are used to detect hotspots or cold spots, identifying high-value clusters, low-value clusters, and insignificantly stable areas of land use or ecological environment change. ; in, This is a spatial clustering statistic for coupling coordination degree. For the coupling coordination degree of pixel q, The global mean of the coupling coordination degree. This represents the global standard deviation of the coupling coordination degree.

[0096] like >0 indicates that pixel p is a high-coordination hotspot area; <0 indicates that pixel p is a low-coordination cold spot area. , indicating that pixel p is a non-significantly stable region.

[0097] S6.3 Select driving factors, standardize the driving factors and perform multicollinearity tests to eliminate redundant factors.

[0098] Driver factor Min-Max normalization (mapped to [0,1]): ; Z-score normalization (mapped to [-1,1]): ; in, The standardized value The original value, The factor mean The standard deviation of the factor is 1.

[0099] S6.4. Using the driving factors selected in step S6.3 as input, random forest, XGBoost and multiple linear regression algorithms are used to construct driving models in high coordination hot spots, low coordination cold spots and non-significant stable areas respectively. The root mean square error, coefficient of determination and mean absolute error are used to compare and evaluate the models. The optimal model is used to calculate the predicted value of coupling coordination degree.

[0100] Based on the spatial partitioning results of high-coordination hotspot areas, low-coordination coldspot areas, and insignificantly stable areas, driving factor datasets and true values ​​of coupling coordination degrees are extracted for each region. Using driving factors as independent variables and coupling coordination degree as dependent variables, regression prediction models are constructed for each region using three algorithms: random forest, XGBoost, and multiple linear regression. The performance of each algorithm is compared by means of root mean square error, coefficient of determination, and mean absolute error. The optimal algorithm and corresponding region driving model are selected for each region, forming three regions-specific driving model systems for coupling coordination degree prediction and subsequent driving mechanism analysis.

[0101] Section S6.5 introduces the SHAP value interpretation method to perform interpretability analysis on the optimal model, quantifying the contribution, direction, and nonlinear relationship of each driving factor. The formula is: ; in, This is the predicted value for coupling coordination. The mean predicted value (baseline value) for all samples is given, where m is the total number of driving factors. is the SHAP value of the j-th driving factor. A positive value indicates a positive contribution, and a negative value indicates a negative contribution. The larger the absolute value, the stronger the contribution.

[0102] Specifically, the driving factor data in step S6 includes: temperature, precipitation, altitude, slope, aspect, soil type, and soil organic matter.

[0103] In this embodiment, the land use and ecological environment coupled evolution analysis method is implemented collaboratively in the Google Earth Engine cloud platform and the Python environment. The batch processing of remote sensing images is performed on the Google Earth Engine platform, while statistical modeling and spatial analysis are performed in the Python environment.

[0104] The research process of this invention follows a closed-loop logic of "multi-source data integration → standardized preprocessing → index and trend quantification → dynamic land use analysis → ecological coupling and driver attribution → verification output." Figure 2 As shown, taking a certain region as the study area and spanning the time period from 2016 to 2024, the specific implementation process of the land use and ecological environment coupled evolution analysis method of this application for the land use and ecological environment coupled evolution analysis of a certain region is as follows.

[0105] Phase 1: Multi-source data input and classification / organization phase To address the research needs of a specific region, systematically collect and archive multi-dimensional data, clarifying the purpose and format requirements of each data point: (1) Remote sensing image data, including: Sentinel-2 multi-temporal imagery: growing season (May-October) images from 2016 to 2024, with a spatial resolution of 10m, including 13 bands such as red, green, blue, and near-infrared, used for ecological index calculation and land use classification; High-resolution satellite imagery (GF-1 / 2): 2m spatial resolution, used to assist in selecting training samples for land use classification and improve sample representativeness.

[0106] (2) Basic spatial data, including: Geographic boundary vector of a certain region: administrative boundaries, wetland reserve boundaries, etc., used to define the research scope; Historical land use data: Land use classification maps before 2016, serving as reference samples for the classification model.

[0107] (3) Auxiliary attribute data, including: Driving factor data: temperature, precipitation, altitude, slope, aspect, soil type, soil organic matter, etc.

[0108] Phase Two: Data Preprocessing Phase Based on the characteristics of Sentinel-2 images, standardization preprocessing was performed to eliminate noise interference. Specific steps included:

[0109] (1) Radiation calibration: The Sentinel-2 radiometric calibration tool of GEE was used to convert the raw DN values ​​of the image into apparent reflectance, thus unifying the radiometric dimensions of the data.

[0110] (2) Atmospheric correction: The Sen2Cor algorithm is used to eliminate the effects of atmospheric scattering and absorption, converting apparent reflectance into true surface reflectance and improving the spectral accuracy of the image.

[0111] (3) Cloud / Shadow / Snow Mask: Based on the QA60 quality band of Sentinel-2, pixels in cloud, cloud shadow, and snow-covered areas are identified and masked, retaining only valid observation pixels without cloud contamination.

[0112] (4) Image cropping and splicing: Based on the boundary vector of a certain region, the images covering the study area are cropped, and multiple images covering the study area are seamlessly stitched together to form a complete image dataset of the study area.

[0113] Phase Three: Ecological Index Calculation and Time-Series Trend Analysis Based on the preprocessed Sentinel-2 images, core ecological indices were calculated in batches, and their interannual variation trends were analyzed:

[0114] (1) Calculation of ecological index: ①Vegetation index: NDVI: Normalized Difference Vegetation Index, the most classic and commonly used vegetation index. It is used to detect vegetation growth status, vegetation cover, and eliminate some radiation errors. The value range is [-1, 1], with positive values ​​indicating vegetation cover and higher values ​​indicating more abundant vegetation; negative values ​​indicate non-vegetated surfaces, such as water bodies, clouds, and snow. It is calculated using the reflectance of the near-infrared (NIR) and red (Red) bands, as shown in the following formula: ; Wherein, NIR represents the reflectance in the near-infrared band, and Red represents the reflectance in the red band.

[0115] EVI: An enhancement of the vegetation index over NDVI. It corrects for atmospheric aerosol and soil background interference by introducing a blue band, making it particularly suitable for high biomass areas. Its dynamic range and sensitivity are better than NDVI. It takes into account the effects of atmospheric and soil background, as shown in the following formula: , Wherein, Blue represents the reflectivity of the blue light band, NIR represents the near-infrared band, and Red represents the red light band.

[0116] FVC (Greenness Value): Vegetation cover is an important biophysical parameter that directly estimates the percentage of area occupied by the vertical projection of vegetation canopy within a pixel. In the provided code, FVC is derived from NDVI using an empirical linear model. This method assumes that the minimum NDVI value represents bare soil and the maximum value represents complete vegetation cover. FVC transforms vegetation "greenness" information (NDVI) into a more intuitive physical quantity and is widely used in eco-hydrological models, soil erosion assessments, urban greening surveys, and climate change research. It is a key indicator for quantifying surface vegetation structure. The calculation uses a pixel-based bipartite model, as shown in the following formula: ; in, and These represent the minimum and maximum NDVI values ​​within the study area, respectively.

[0117] ② Water body index: NDWI (Normalized Difference Water Index) identifies water bodies by utilizing their high reflectivity in the green light band and extremely low reflectivity in the near-infrared band due to strong absorption. It is effective in identifying large areas of open water. However, a major drawback of NDWI is its susceptibility to confusion with built-up areas, as built-up areas also have relatively high reflectivity in the green light band, potentially leading to misclassification of "false water bodies" in urban areas. Furthermore, it is less sensitive to water bodies with vegetation cover, as the near-infrared reflection from vegetation weakens the water body signal. The formula for calculating reflectivity in the green and near-infrared bands is as follows: ; NIR stands for near-infrared band, and Green stands for green band.

[0118] MNDWI: The Improved Normalized Water Index is a significant optimization addressing the shortcomings of NDWI. It replaces the near-infrared band with the short-wave infrared band because both water and vegetation have low reflectivity in the short-wave infrared, while buildings and dry bare soil have very high reflectivity. This substitution allows MNDWI to effectively suppress background noise from urban areas and soil, thus more accurately delineating water body boundaries, especially in complex areas where urban environments intersect with water bodies. MNDWI performs exceptionally well in monitoring small rivers, reservoirs, coastline changes, and eutrophication, making it one of the preferred indices for water body information extraction. MNDWI is an improvement on NDWI, using the short-wave infrared band (SWIR1) instead of the near-infrared band, as shown in the formula below: ; Among them, SWIR1 is the shortwave infrared band, and Green is the green band;

[0119] ③ Environmental indices: SI (Salinity Index) is a specific index used to indicate the degree of soil salinization. Saline soils often have unique spectral characteristics; for example, the formation of a salt crust on the surface leads to increased reflectance in specific wavelengths (such as shortwave infrared). SI amplifies this signal by calculating the ratio of shortwave infrared to near-infrared. A higher index value generally indicates higher surface salinity and worse soil conditions. It plays an important role in surveying saline-alkali land in arid and semi-arid regions, monitoring secondary salinization of farmland, and assessing ecological degradation. The index is calculated using reflectance in the near-infrared (NIR) and shortwave infrared (SWIR1) bands, as shown in the following formula: ; Among them, SWIR1 is the shortwave infrared band, and NIR is the near-infrared band.

[0120] NDSI (Normalized Difference Building Index / Salinity Index) is a multi-purpose index. It utilizes the fact that buildings and impermeable layers have higher reflectance in the shortwave infrared than in the near-infrared, effectively distinguishing urban built-up areas from natural surfaces. Simultaneously, because bare land and certain saline-alkali lands have similar spectral characteristics to built-up areas, NDSI is also frequently used as a bare land index to track urban sprawl, desertification, and bare surfaces with extremely low vegetation cover. It is calculated using reflectance in the visible light band (red band) and the mid-infrared band (MIR), as shown in the following formula: ; Among them, Red represents the red light band, and MIR represents the mid-infrared band.

[0121] NDBSI (Normalized Difference Building and Soil Index) is a comprehensive index that combines multiple bands of information sensitive to buildings and bare land. Its design aims to more comprehensively characterize surface conditions under urban environments and human disturbance. A higher NDBSI value indicates more pronounced building or bare land characteristics within a pixel. This index is a key component in constructing the Remote Sensing Ecological Index (RSEI). In the comprehensive evaluation of regional ecological environment quality, it quantifies the distribution of buildings and bare land representing "dry" stress, complementing and contrasting with vegetation cover (representing "green") and humidity (representing "wet"). The formula is as follows: ; Among them, SWIR1 is the shortwave infrared band, Red is the red light band, NIR is the near-infrared band, and Blue is the blue light band.

[0122] (2) Time series trend analysis: ① The Mann-Kendall nonparametric test was used to determine whether the interannual variation trend of each ecological index was significant (P<0.05 was considered significant). ② Combining Sen's Slope calculation, the rate of change of the index (such as the annual growth / decrease rate of NDVI) is quantified, and finally, trend significance plots and rate of change plots of each index are generated.

[0123] Phase Four: Land Use Classification and Dynamic Change Monitoring Phase Based on remote sensing imagery and ecological indices, annual classification and change analysis of land use were completed.

[0124] (1) Classification system and sample selection: ① Construct a classification system specific to a certain region: divide it into 7 categories: cultivated land, forest land, grassland, wetland (coastal wetland / inland wetland), construction land, unused land, and water bodies; ② Sample selection: Combining GF-1 / 2 high-resolution imagery with field survey data, select at least 500 sample points for each class in GEE (to ensure uniform spatial distribution of samples).

[0125] (2) Classification model construction and optimization: ① A classification model is constructed using the random forest algorithm, with input features such as spectral bands and ecological indices from Sentinel-2 images; ② Optimize model parameters (such as the number of decision trees) through 5-fold cross-validation to ensure that the overall model accuracy is ≥85% and the Kappa coefficient is ≥0.8.

[0126] (3) Annual land use classification: The optimized model was applied to preprocessed images from 2016 to 2024 to generate annual land use classification raster maps.

[0127] (4) Change detection and analysis: Calculate the land use transfer matrix: By overlaying the classification maps year by year, the direction and area of ​​transfer in / out for each type of land use are clarified.

[0128] (5) Extracting change indicators: Calculate land use dynamics (representing the rate of change) and expansion intensity index (representing the degree of expansion of types such as construction land) to identify hotspots of change (such as wetland shrinkage areas and construction land expansion areas).

[0129] Phase Five: Land Use-Ecological Environment Coupling Analysis (1) Single-factor response analysis: For different land use conversion types (such as "wetland → farmland" and "farmland → construction land"), the changes in indices such as NDVI, SI, and NDWI within the conversion area are statistically analyzed, and the direction of the impact (improvement / degradation) of each type of conversion on ecological elements is analyzed.

[0130] (2) Calculation of coupling coordination degree: The application process of the coupling coordination degree model is optimized as follows: First, based on the regional characteristics of a certain area, a coupled evaluation index system of "land use-ecological environment" is constructed. The land use system focuses on the core of type and change, and selects the proportion of land use type (the area proportion of key types such as cultivated land, wetland, and construction land), land use dynamics (representing the annual change rate of type), and core transformation area of ​​the transfer matrix (such as the area of ​​major change paths such as "wetland → cultivated land" and "cultivated land → construction land") as evaluation indicators. The ecological environment system still selects the annual average value of NDVI, the annual average value of SI, and the annual average value of NDWI as core indicators.

[0131] Next, the indicators are distinguished by their attributes. Positive indicators such as land use dynamics (positive within a reasonable range) and annual average NDVI, as well as negative indicators such as annual average SI and construction land ratio, are processed using corresponding standardized formulas to eliminate the influence of dimensional and numerical differences.

[0132] Subsequently, the weights of each indicator were objectively determined using the entropy weight method to avoid subjective assignment bias, and the comprehensive evaluation value of the land use system (U1) and the comprehensive evaluation value of the ecological environment system (U2) were calculated. ; ; In the formula, U1 is the comprehensive evaluation value of the land use system, and U2 is the comprehensive evaluation value of the ecological environment system. Let j be the weight of the j-th indicator in the land use system. Let k be the weight of the indicator in the ecological environment system. Let j be the standardized value of the j-th indicator in the land use system. Let k be the standardized value of the indicator in the ecological environment system. This represents the total number of indicators in the land use system. This represents the total number of indicators for the ecological and environmental system.

[0133] Then, the interaction strength (C value, range 0-1) of the two systems is calculated using the coupling degree formula: Introducing the comprehensive development level (T value, α=β=0.5, taking into account the importance of both systems): , Where α and β are weighting coefficients, and α+β=1.

[0134] The coupling coordination degree is derived. The coupling coordination level is determined based on the D value.

[0135] Finally, based on the degree of coupling coordination, the coupling coordination state is divided into 3 categories and 8 levels: 0≤D<0.4 is the imbalance and decline category, 0.4≤D<0.6 is the transitional and harmonious category, and 0.6≤D≤1 is the coordination and improvement category. The matching status and changing trend of land use type, change characteristics and ecological environment in a certain region from 2016 to 2024 are accurately quantified, providing quantitative support for regional land use optimization and ecological protection.

[0136] Phase Six: Driver Factor Analysis Phase (1) Driving factors include: temperature, precipitation, altitude, slope, aspect, soil type, and soil organic matter; (2) Factor analysis and modeling: ① Preprocessing: Pearson correlation coefficient was used to screen for factors without multicollinearity, and principal component analysis was used to reduce the dimensionality of multiple factors to a few principal components; ② Spatial modeling: The geographically weighted regression (GWR) model is used to quantify the contribution of each driving factor in different regions and identify the dominant driving factors in the study area (e.g., natural factors dominate in coastal areas and human factors dominate in urban areas).

[0137] This invention integrates GEE parallel computing, Sen2Cor correction, and QA60 masking, significantly improving the efficiency and accuracy of long-term data processing. It pioneered a full-chain analysis of "classification change → single-factor response → coupling coordination → spatial driving" and revealed the spatial heterogeneous influence of driving factors, providing scientific support for the territorial spatial planning of coastal wetland areas.

[0138] like Figure 3 As shown, the land use transfer matrix heatmap transforms the numerical land use transfer matrix into an intuitive visual representation through color gradients, enhancing the efficiency of data feature recognition. A two-color gradient scheme is used, with the change from light to dark colors corresponding to an increase in the number of transfer pixels; the darker the color, the larger the scale of the transfer path. Observing the overall color distribution, the color blocks in the diagonal positions are generally darker, reflecting the stable dominance of each land use type. In the off-diagonal areas, the transfer blocks from cultivated land to water areas and from cultivated land to construction land are the darkest, confirming the importance of these two transformation paths. The colors of the forest land rows and columns are generally lighter, consistent with its low transformation rate. The transfer matrix heatmap can quickly identify the main directions of land use type changes, thereby determining the main land use type directions identified by the model.

[0139] like Figure 4 As shown, the chart compares the percentage composition of six land use types—arable land, forest, grassland, water area, unused land, and water area—between 2016 and 2024. By comparing the data from the two years, the changes in the proportion of each land use type in the total area are clearly visible, such as the potential increase in construction land and the potential decrease in arable land or forest. This chart helps to identify the main directions of regional land use transformation from a macro-structural perspective, such as the expansion of construction land due to urbanization or the restoration of forests and grasslands brought about by ecological protection policies. This structural analysis is fundamental to understanding the evolution of regional human-land relationships and provides background information for subsequent changes in spatial patterns and landscape indices.

[0140] like Figure 5The final comparison results show that the method of this invention performs excellently in terms of land use identification accuracy, ecological index inversion accuracy, and driving factor analysis reliability. All evaluation indicators are higher than those of single algorithms and traditional statistical analysis methods. Among them, the overall accuracy of land use classification reaches a high level, and key indicators such as the significance test accuracy of ecological trend analysis and the explanatory power of driving factors are also significantly improved. This fully demonstrates that the technical system constructed in this study is suitable for the ecological environment research characteristics of a certain wetland.

[0141] Example 2 like Figure 2 As shown, an embodiment of the present invention provides a system for coupled evolution analysis of land use and ecological environment. Utilizing the method provided in Embodiment 1, the system includes: The data acquisition module is used to acquire multi-temporal remote sensing images of the study area and auxiliary data, including driving factor data. The preprocessing module is used to preprocess multi-temporal remote sensing images to obtain a set of surface reflectance images; The index calculation module is used to calculate at least one ecological index based on the preprocessed image, generate time series data of the ecological index, and perform interannual trend analysis on the ecological index. The land use classification module is used to classify land use in the study area based on preprocessed images and ecological indices, and obtain annual land use data. The coupling analysis module is used to extract land use change characteristics based on annual land use data and calculate the coupling coordination degree between land use and the ecological environment. The driver identification module is used to combine driver factor data and coupling coordination degree to identify the spatial differentiation characteristics of the influence of driver factors on coupling coordination degree.

[0142] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for coupled evolution analysis of land use and ecological environment, characterized in that, Specifically, the steps include the following: S1. Acquire multi-temporal remote sensing images of the study area and auxiliary data, including DEM data and driving factor data. S2, preprocessing multi-temporal remote sensing images to obtain a set of surface reflectance images; S3. Calculate at least one ecological index based on the preprocessed image, generate time series data of the ecological index, and perform interannual trend analysis on the ecological index. S4. Based on the preprocessed images and ecological indices, land use is classified in the study area using a classification model to obtain annual land use data. S5. Based on annual land use data, extract land use change characteristics and calculate the coupling coordination degree between land use and the ecological environment. S6. By combining driving factor data and coupling coordination degree, the spatial differentiation characteristics of the influence of driving factors on coupling coordination degree are identified, and the predicted land evolution results are obtained. Step S5 specifically includes the following steps: S5.1, using the zoning statistical method, calculate the mean, standard deviation, and coefficient of variation of ecological indicators under each land use type; S5.2 uses Pearson correlation coefficient to perform correlation analysis, quantifying the association strength between land use type and ecological indicators. The calculation formula is as follows: ； in, The correlation coefficient, These represent land use type codes and ecological indicator values, respectively, where m is the sample size. for The average value, for The average value; S5.3, the land use type distribution map and the ecological indicator spatial distribution map are rasterized to generate a composite layer of land use and ecological indicators. The expression is: ； in, For the coupled layer pixel values, For land use type pixel values, For the first The pixel value of each ecological indicator This is a spatial superposition operation; S5.4, using Getis-Ord Hotspot or coldspot analysis is conducted to identify high-value and low-value clusters of ecological indicators. The calculation formula is as follows: ； in, This is a spatial clustering statistic. This is the spatial weight matrix. For pixel index values, S is the mean of the indicator, and S is the standard deviation of the indicator. Step S5 also includes the following steps: S5.5, Greenness NDVI, Humidity Wet, Dryness NDBSI, and Heat LST are selected as evaluation factors. Each evaluation factor is normalized to the [0,1] interval. Positive factors are: ； Negative factors are adopted as follows: ； Where NI is the standardized value, I is the original value, and I max I min These are the maximum and minimum values ​​of the factor, respectively; S5.6, perform principal component analysis on the standardized factors, using the first principal component PC1 as the initial value of RSEI. The calculation formula is as follows: ； in, , The first principal component loading coefficient; Normalizing PC1 yields the final RSEI index, calculated using the following formula: ； The RSEI ranges from 0 to 1, with a higher value indicating better regional ecological quality. S5.7 performs overlay analysis on the annual land use classification map, calculates the land use transfer matrix, and clarifies the transfer-in area, transfer-out area, and conversion direction of various land features; based on the transfer matrix, it calculates the land use dynamic degree and expansion intensity index to identify hotspot areas of land use change; S5.8, construct evaluation indicators for the land use system and the ecological environment system. The land use system evaluation indicators include the area proportion of various land features, land use dynamics, and the area of ​​key transformation types. The ecological environment system evaluation indicators include the mean vegetation index, mean salinity index, mean water body index, and remote sensing ecological index (RSEI). Range standardization is performed on each evaluation indicator. Positive indicators adopt: ； Negative indicators are adopted as follows: ； in, The standardized values ​​of the evaluation indicators The original values ​​of the evaluation indicators, For the first The global maximum value of the evaluation index For the first The global minimum value of the evaluation index; The weights, information entropy, and difference coefficients of the evaluation indicators are calculated using the entropy weight method to obtain the weights of each evaluation indicator. ； ； ； ； in, For the first Evaluation indicator number The weight of each evaluation unit For the first Information entropy of the evaluation indicators For the first The coefficient of difference of the evaluation indicators For the first The entropy weight of each evaluation indicator The total number of evaluation units, The total number of evaluation indicators, It is the natural logarithm function; Step S5 also includes the following steps: S5.9, calculate the comprehensive evaluation value of the land use system respectively. Comprehensive evaluation value of ecological environment system The formula is: ； ； in, This is the comprehensive evaluation value of the land use system. This is a comprehensive evaluation value for the ecological environment system. For the land use system The weight of each evaluation indicator For the ecological environment system The weight of each evaluation indicator For the land use system Standardized values ​​of the evaluation indicators For the ecological environment system Standardized values ​​of the evaluation indicators This represents the total number of indicators in the land use system. This represents the total number of indicators for the ecological and environmental system. S5.10, Calculate the coupling degree C between the land use system and the ecological environment system according to the coupling degree formula: ； Where C represents the coupling degree, taking values ​​[0,1], with larger values ​​indicating higher coupling. S5.11, Calculate the overall development level T: ； Where T represents the comprehensive development level of land use and ecological environment. , These are the weighting coefficients; S5.12, Calculate the coupling coordination degree D, and classify the coupling coordination level according to the D value: ； Where D is the coupling coordination degree, with a value of [0,1]; Step S6 specifically includes the following steps: S6.1 Select patch area, patch shape index, fractal dimension, patch density, edge density, and aggregation index to conduct a quantitative analysis of the regional landscape spatial pattern; S6.2 uses Getis-Ord Spatial statistics are used to detect hotspots or cold spots, identifying high-value clusters, low-value clusters, and insignificantly stable areas of land use or ecological environment change. ； in, This is a spatial clustering statistic for coupling coordination degree. For the coupling coordination degree of pixel q, The global mean of the coupling coordination degree. The global standard deviation of the coupling coordination degree; like >0 indicates that pixel p is a high-coordination hotspot area; <0 indicates that pixel p is a low-coordination cold spot area. , indicating that pixel p is a non-significantly stable region; S6.3 Select driving factors, standardize the driving factors and perform multicollinearity test to eliminate redundant factors; S6.

4. Using the driving factors selected in step S6.3 as input, random forest, XGBoost and multiple linear regression algorithms are used to construct driving models in high coordination hot spots, low coordination cold spots and non-significant stable areas respectively. The root mean square error, coefficient of determination and mean absolute error are used to compare and evaluate the models. The optimal model is used to calculate the predicted value of coupling coordination degree. S6.5 introduces the SHAP value interpretation method to perform interpretability analysis on the optimal model, quantifying the contribution, direction, and nonlinear relationship of each driving factor; the formula is: ； in, This is the predicted value for coupling coordination. The average predicted value for all samples is given by m, where m is the total number of driving factors. Let be the SHAP value of the j-th driving factor. A positive value indicates a positive contribution, and a negative value indicates a negative contribution. The larger the absolute value, the stronger the contribution. The driving factor data in step S6 include: temperature, precipitation, altitude, slope, aspect, soil type, and soil organic matter.

2. The method for coupled evolution analysis of land use and ecological environment according to claim 1, characterized in that, Specifically, the steps include the following: S2.1 Radiometric calibration of multi-temporal remote sensing images is performed to convert the original DN values ​​into apparent reflectance; S2.2, The Sen2Cor algorithm is used to perform atmospheric correction on the radiometrically calibrated image to obtain the true surface reflectance; S2.3, based on the image quality band, performs cloud, shadow and snow pixel masking on the atmospherically corrected image to retain effective observation pixels; S2.4, based on the boundary vector of the study area, the masked images are cropped and stitched together to form a set of surface reflectance images covering the study area.

3. The method for coupled evolution analysis of land use and ecological environment according to claim 1, characterized in that, Ecological indices include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Free Vegetation Capacity (FVC), Normalized Difference Water Index (NDWI), Improved Normalized Difference Water Index (MNDWI), Salinity Index (SI), Normalized Difference Building Index (NDSI), and Normalized Difference Building and Soil Index (NDBSI). The formulas for calculating these ecological indices are as follows: ； ； ； ； ； ； ； ； in, For near-infrared reflectivity, Reflectivity in the red light band For blue light band reflectivity, For green light band reflectivity, For shortwave infrared reflectance, For mid-infrared reflectivity, and These represent the minimum and maximum NDVI values ​​within the study area, respectively.

4. The method for coupled evolution analysis of land use and ecological environment according to claim 1, characterized in that, Step S3 specifically includes the following steps: S3.1, the Mann-Kendall nonparametric test was used to determine the significance of the interannual variation trend of each ecological index at the pixel level; S3.2, calculate the pixel-level change rate of each ecological index by combining Sen's Slope slope; S3.3 generates a distribution map of the significance of ecological index trends and a distribution map of the rate of change.

5. The method for coupled evolution analysis of land use and ecological environment according to claim 1, characterized in that, Step S4 specifically includes the following steps: S4.1 Construct a land use classification system applicable to the study area. The classification system includes: cultivated land, forest land, grassland, wetland, construction land, unused land, and water bodies. S4.2 Select training samples for each type of land cover, with a sample size of not less than N and a uniform spatial distribution; S4.3, using the spectral bands and ecological indices of the preprocessed image as input features, train a random forest classification model; S4.4, optimize the parameters of the random forest classification model through cross-validation; S4.5 applies the optimized model to the preprocessed images year by year to obtain the land use classification raster map year by year.

Images