An ecological quality evaluation method based on an improved index set

Abstract

The application discloses an ecological quality evaluation method based on an improved index set, comprising remote sensing image acquisition and preprocessing, improved ecological index set construction, spatial point set construction, ecological index standardization and geographical weighted principal component analysis, post-processing and result output; the application expands the traditional two-dimensional dryness (SI+IBI) / 2 into three-dimensional evidence (SI+IBI+EI) / 3, and additionally adds NDTI as the fifth dimension, so that the model can distinguish between hardening type dryness and ploughing bareness / soil disturbance type dryness, and the bare disturbance signal is no longer pressed into the error term of SI in the county dominated by farmland. kNDVI alleviates the saturation of NDVI through nonlinear stretching, so that the greenness expression of the high coverage area and the bare land grassland transition zone is more consistent.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing image processing and ecological environment assessment technology, specifically to an ecological quality assessment method based on an improved index set. Background Technology

[0002] The Remote Sensing Ecological Index (RSEI) incorporates four representative indicators—Greenness Index (NDVI), Humidity Index (WET), Dryness Index (NDBSI), and Heat Index (LST)—into a single framework. Principal Component Analysis (PCA) is used to compress multidimensional information into a single principal component (PC1), thereby achieving comprehensive quantification and visualization of regional ecological quality. This index requires no manual weighting and is entirely driven by the covariance structure of the data itself, offering advantages such as objectivity and repeatability. Currently, RSEI and its various improved forms have been widely validated across different scales and landscape types, including watersheds, cities, arid regions, and mining areas, becoming one of the mainstream methods in the field of remote sensing ecological assessment.

[0003] However, the existing RSEI indicator system still has two shortcomings. First, the traditional dryness index NDBSI usually uses the calculation method of (SI+IBI) / 2, where the spectral responses of SI (bare soil index) and IBI (building index) are mutually contaminated in bright bare topsoil, building corners, and dry vegetation areas, leading to distortion of the dryness signal. Second, the traditional greenness index NDVI tends to saturate in high-coverage vegetation areas and cannot sensitively reflect changes in the actual vegetation condition. To address these issues, this invention introduces the Normalized Farm Disturbance Index (NDTI) into the dryness index to more accurately capture the ecological degradation pressure represented by exposed tillage, sparse stubble, and dry topsoil in counties dominated by arable land, rather than simply being indirectly accommodated by SI or IBI; at the same time, the greenness index can adopt alternative indices with better saturation characteristics, such as EVI, as appropriate, thereby constructing an improved indicator set with greater ecological explanatory power.

[0004] In addition to the indicator system, traditional RSEI employs Global Principal Component Analysis (PCA), which calculates a fixed set of loading coefficients based on all pixels within the study area. This implicitly assumes that the main control mechanism of ecological quality is identical across all regions, and that ecological degradation or improvement is driven by the same factors with the same contribution rate regardless of geographical location. This assumption of spatial stationarity often leads to distortion of local ecological assessment results in heterogeneous regions (especially counties with complex topography, diverse land use, or significant ecological gradients). While recent studies have attempted to address this through methods such as zoning correction and local entropy weighting, these either still require empirical zoning or lose the explanatory power of the original covariance structure of PCA. Therefore, this invention introduces Geographically Weighted Principal Component Analysis (GWPCA): a local neighborhood is established based on spatial distance for each sample point, and principal component loadings are estimated within each local window, allowing for differentiated indicator contribution structures at different locations. This method, while fully preserving the covariance interpretation logic of PCA, directly characterizes the spatial non-stationarity of ecological driving factors. Combined with the aforementioned improved indicator set, it achieves a more adaptive and accurate ecological quality assessment. Summary of the Invention

[0005] The purpose of this invention is to provide an ecological quality assessment method based on an improved index set to solve the problems existing in the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: an ecological quality assessment method based on an improved index set, comprising the following steps:

[0007] S1: Acquire multispectral remote sensing images and perform radiometric calibration, atmospheric correction and geometric correction preprocessing;

[0008] S2: Based on the aforementioned multispectral remote sensing imagery, calculate an ecological indicator dataset comprising five indicators: greenness, humidity, dryness, heat, and tillage index; wherein, the greenness indicator uses kNDVI, and the calculation formula is as follows:

[0009]

[0010] In the formula, For near-infrared reflectivity, Reflectivity in the red light band;

[0011] The humidity index uses the tasseled cap transformation (WET), and the calculation formula is as follows:

[0012]

[0013] In the formula, These are the blue, green, red, near-infrared, SWIR1, and SWIR2 bands, respectively.

[0014] The dryness index used is the modified normalized dryness index (NDBSI), calculated using the following formula:

[0015]

[0016] The formula for calculating the bare soil index SI is as follows:

[0017]

[0018] The formula for calculating the Building Index (IBI) is as follows:

[0019]

[0020] The formula for calculating the visible contrast index (EI) is:

[0021]

[0022] The thermal index uses land surface temperature (LST), and the calculation formula is as follows:

[0023]

[0024] In the formula, Brightness temperature obtained by inversion of remote sensing images in the thermal infrared band;

[0025] The tillage index uses the Normalized Tillage Disturbance Index (NDTI), calculated using the following formula:

[0026]

[0027] In the formula, and The reflectance values ​​are SWIR1 and SWIR2 bands, respectively; a five-dimensional ecological index set is constructed based on this, including five indicators: kNDVI, WET, NDBSI, LST, and NDTI.

[0028] S3: Constructing a spatial point set: The ecological indicator dataset is sampled according to a regular grid, and the planar coordinates of each sampling point are extracted. The longitude and latitude are converted to projected coordinates to form a multidimensional dataset with coordinates.

[0029] S4: Standardize and unify the direction of the five-dimensional ecological index set: keep the positive indices kNDVI and WET at their original values, and multiply the negative indices NDBSI, LST, and NDTI by -1 to achieve positive orientation; use Z-score standardization to scale each index to a mean of 0 and a standard deviation of 1; use the geographic weighted principal component analysis (GWPCA) model to construct local neighborhoods with adaptive bandwidth, calculate the local covariance matrix and loading vector, and obtain the first principal component PC1 score for each sampling point; ensure that PC1 and kNDVI are aligned through anchoring operations;

[0030] S5: Interpolate the PC1 score to a raster image, and use Min-Max normalization to scale the value range to [0,1] to obtain the ecological quality index raster image; according to the ecological index value range [0,1], divide it into five levels of poor, poor, moderate, good and excellent with an interval of 0.2, and output the ecological quality spatial distribution map.

[0031] Furthermore, in step S2, kNDVI adopts a hyperbolic tangent function form to transform the nonlinear saturation characteristics of NDVI into a continuous and smooth spatial response, thereby alleviating the exponential saturation problem in high-density vegetation areas.

[0032] Furthermore, in step S2, EI is used to enhance the evidence constraint of dryness index on bright bare topsoil and artificial bright surface, making up for the shortcomings of traditional NDBSI which only relies on SI and IBI; NDTI, as the fifth independent channel, is used to distinguish between hardened drying and cultivated bare drying.

[0033] Furthermore, in step S4, the formula for calculating the Z-score standardization is as follows:

[0034]

[0035] In the formula, Let be the original value of the i-th indicator. and These represent the mean and standard deviation of the indicator across all sampling points, respectively; the positiveification of negative indicators is performed before standardization.

[0036] Furthermore, in step S4, the adaptive bandwidth of GWPCA adopts the nearest neighbor method, and the local neighborhood of each sampling point consists of its 200 nearest points in space; the local load vector is obtained by eigenvalue decomposition of the local covariance matrix, and the first principal component PC1 corresponds to the direction of the largest eigenvalue.

[0037] Furthermore, in step S4, the determination condition for the PC1 anchoring operation is: calculate the correlation coefficient between the PC1 sequence and the kNDVI sequence; if the correlation coefficient is negative, multiply PC1 by -1; to ensure that the final PC1 value is positively correlated with kNDVI.

[0038] Furthermore, in step S5, the interpolation from PC1 to the grid adopts the inverse distance weighted IDW method, with a power exponent of 2, a search radius of 3000m, and a maximum number of neighboring points of 12.

[0039] Furthermore, in step S1, the spatial resolution of the multispectral remote sensing image is better than or equal to 30m.

[0040] Furthermore, in step S1, when using multispectral surface reflectance products, the formula for converting the original DN value into surface reflectance is as follows:

[0041]

[0042] In the formula, Let be the surface reflectance of the i-th optical band. This is the original DN value for this band; the formula for converting the original DN value of the thermal infrared band to brightness temperature is:

[0043]

[0044] In the formula, For brightness temperature, This represents the original DN value for the thermal infrared band.

[0045] Compared with the prior art, the beneficial effects of the present invention are:

[0046] This invention expands the traditional two-dimensional dryness index (SI+IBI) / 2 to a three-dimensional evidence index (SI+IBI+EI) / 3, and adds NDTI as a fifth dimension, enabling the model to distinguish between hardened drying and bare / topsoil disturbed drying. In counties with predominantly arable land, the bare soil disturbance signal is no longer pushed into the error term of SI. kNDVI alleviates NDVI saturation through nonlinear stretching, making the greenness expression more consistent between high-coverage areas and bare grassland transition zones. GWPCA replaces global PCA, with the feature vector changing with location, and ±1 inversion is eliminated point-by-point through kNDVI anchoring. It has a high degree of automation and is suitable for large-scale applications. The PC1 load is entirely driven by local covariance data, eliminating the need for manual determination of which indicator is more important, making it suitable for cross-regional replication. Attached Figure Description

[0047] Figure 1 This is a flowchart illustrating the overall technical process of the method of the present invention.

[0048] Figure 2 This is the dataset of ecological indicators calculated in this invention. Figure 2 (a) is the kNDVI indicator chart, (b) is the WET indicator chart, (c) is the NDBSI indicator chart, (d) is the LST indicator chart, and (e) is the NDTI indicator chart.

[0049] Figure 3 This is a comparison chart of the results of the present invention (MRSEI) and the results of the conventional RSEI in the study area 1;

[0050] Figure 4 This is a comparison chart of the results of the present invention (MRSEI) and the results of the conventional RSEI for the study area 2. Detailed Implementation

[0051] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. 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 of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0052] Please see Figure 1-4 This invention provides an ecological quality assessment method based on an improved index set, comprising the following steps:

[0053] S1: Remote Sensing Image Acquisition and Preprocessing

[0054] The quality of remote sensing images directly determines the accuracy of subsequent ecological index calculations and evaluations. Therefore, the core objective of this step is to acquire high-quality multispectral remote sensing images covering the study area and eliminate radiometric errors, geometric errors, and atmospheric interference through standardized preprocessing to obtain radiometrically consistent and geometrically accurate surface reflectance and surface temperature data.

[0055] This invention supports multiple mainstream multispectral remote sensing data sources and requires the spatial resolution of the images to be better than or equal to 30m in order to ensure that the types of land features and ecological characteristics at the county scale can be clearly identified.

[0056] In this embodiment, publicly available standardized products are selected, and system-level radiometric calibration, atmospheric path radiometric correction and aerosol correction have been completed. Standardized surface reflectance data are directly output, eliminating the need for users to run complex atmospheric correction models such as 6S or FLAASH, which greatly improves the efficiency and consistency of data preprocessing. This product includes the following core bands: SR_B2 (blue light band, wavelength 0.45-0.51μm, spatial resolution 30m), SR_B3 (green light band, wavelength 0.53-0.59μm, spatial resolution 30m), SR_B4 (red light band, wavelength 0.64-0.67μm, spatial resolution 30m), SR_B5 (near-infrared band, wavelength 0.85-0.88μm, spatial resolution 30m), SR_B6 (SWIR1 band, wavelength 1.57-1.65μm, spatial resolution 30m), SR_B7 (SWIR2 band, wavelength 2.11-2.29μm, spatial resolution 30m), as well as the thermal infrared band ST_B10 (wavelength 10.6-11.2μm, spatial resolution 30m) and the pixel quality evaluation band QA_PIXEL.

[0057] The specific steps of data preprocessing are as follows:

[0058] S1.1: DN value physical quantity conversion

[0059] The raw data for the instance is stored as a 16-bit integer DN value, which needs to be converted into physically meaningful surface reflectance and brightness temperature using the officially provided scaling factors. For the optical bands (SR_B2 to SR_B7), the formula for calculating surface reflectance is:

[0060]

[0061] In the formula, Let be the surface reflectance of the i-th optical band, with a value ranging from [-0.2, 1.6]. This represents the original DN value for this band, ranging from [0, 65535]. The coefficient 0.0000275 in this formula is the gain coefficient, and -0.2 is the offset coefficient. For the thermal infrared band (ST_B10), the brightness temperature is calculated using the following formula:

[0062]

[0063] In the formula, The brightness temperature of the top layer of the atmosphere is expressed in Kelvin (K). This represents the original DN value in the thermal infrared band, with a value range of [0, 65535].

[0064] S1.2: Low-quality pixel mask

[0065] By utilizing the bit marker information of the QA_PIXEL band, low-quality pixels such as clouds, cloud shadows, snow, and filler pixels in the image are removed to avoid interference from these pixels in subsequent ecological indicator calculations. Each bit of the QA_PIXEL band represents a different land feature or quality information, where:

[0066] Bit0: Fill pixel, a value of 1 indicates that the pixel is a fill value and has no valid data;

[0067] Bit3: Cloud, a value of 1 indicates that the pixel is covered by a cloud;

[0068] Bit4: Cloud shadow, a value of 1 indicates that the pixel is in the shadow of a cloud;

[0069] Bit5: Snow, a value of 1 indicates that the pixel is covered by snow.

[0070] In this embodiment, Bit0, Bit3, and Bit4 are set as masks, that is, all filling pixels, cloud pixels, and cloud image pixels are removed; at the same time, the surface reflectance data is truncated, and reflectance values ​​less than 0 or greater than 1 are set as invalid values ​​to eliminate the influence of outliers.

[0071] S1.3: Water Mask

[0072] Because the ecological characteristics of aquatic areas differ significantly from those of terrestrial areas, and the spectral properties of water can lead to anomalies in indicators such as greenness and dryness, it is necessary to perform separate masking on aquatic areas. This embodiment uses the modified Normalized Difference Water Index (MNDWI) for water body extraction; the calculation formula is as follows:

[0073]

[0074] In the formula, For green light band reflectivity, This refers to the reflectivity of the SWIR1 band. By setting an appropriate threshold (usually 0), pixels with MNDWI greater than 0 are identified as water bodies and masked.

[0075] S1.4: Projection Transformation and Cropping

[0076] The preprocessed multispectral and land surface temperature images were uniformly converted to a projection coordinate system suitable for the study area, such as the Gauss-Kruger projection or UTM projection commonly used in China, to ensure spatial coordinate consistency. Then, the images were cropped using the vector boundaries of the study area to obtain standardized multispectral and land surface temperature images containing only the study area, which were used for subsequent ecological index calculations.

[0077] S2: Construction of an Improved Ecological Indicator Set

[0078] This step, based on preprocessed multispectral and surface temperature images, calculates an improved ecological index set comprising five indicators: greenness, humidity, dryness, heat, and tillage index. Compared to the four indicators of the traditional RSEI, this invention's index set features targeted improvements in greenness and dryness, and adds a new NDTI index specifically characterizing tillage disturbance, enabling a more comprehensive and accurate reflection of the ecological status at the county scale. For ease of subsequent calculations, the reflectance of each band is first uniformly named:

[0079] make (Blue light band reflectivity) (Reflectivity in the green light band) (Reflectivity in the red light band) (Near-infrared reflectivity) (SWIR1 band reflectivity). (SWIR2 band reflectivity).

[0080] S2.1: Calculation of greenness index kNDVI

[0081] Greenness is a core indicator reflecting vegetation growth and coverage, directly related to the primary productivity and ecosystem service functions of a regional ecosystem. The traditional formula for calculating the Greenness Index (NDVI) is:

[0082]

[0083] The NDVI value ranges from -1 to 1. Generally, the higher the vegetation cover, the larger the NDVI value. However, NDVI has two inherent drawbacks: first, in high-coverage vegetation areas (NDVI > 0.7), a significant saturation phenomenon occurs, meaning that as vegetation biomass increases, the growth rate of NDVI slows down significantly, failing to sensitively capture subtle changes in vegetation growth in high-coverage areas; second, in low-coverage vegetation areas, NDVI is easily affected by soil background reflection, leading to deviations in greenness representation. To address these issues, this invention uses kNDVI instead of traditional NDVI as the greenness indicator, calculated using the following formula:

[0084]

[0085] In the formula, tanh is the hyperbolic tangent function, and its mathematical expression is:

[0086]

[0087] The hyperbolic tangent function has the following characteristics: when x is small, tanh(x) ≈ x, and the function is approximately linear; as x increases, tanh(x) gradually approaches 1, but the growth rate is slower than that of a linear function. By performing a hyperbolic tangent transformation on NDVI and then squaring it, kNDVI can transform the nonlinear saturation characteristics of NDVI into a continuous and smooth spatial response: in low-coverage vegetation areas, kNDVI is approximately linearly correlated with NDVI, maintaining sensitivity to soil background; in high-coverage vegetation areas, the response curve of kNDVI is steeper, effectively alleviating the saturation problem and improving the ability to identify vegetation changes in high-coverage areas. Simultaneously, the squaring operation can amplify differences in greenness, making the contrast between areas with different vegetation cover more obvious.

[0088] S2.2: Humidity index WET calculation

[0089] Humidity is an important indicator reflecting the surface moisture content, including soil moisture, vegetation water content, and water bodies. This invention uses the tasseled cap transform (WET) component as the humidity index. The tasseled cap transform is a linear transformation method based on the statistical characteristics of multispectral images, which can convert the original multispectral bands into components with clear physical meaning. Among them, the humidity component can effectively reflect the surface moisture status. The formula for calculating WET is:

[0090]

[0091] The coefficients in this formula refer to the tasseled cap transformation coefficients, proposed by Crist and Cicone in 1984, and have been widely used in remote sensing monitoring of surface moisture. A higher WET value indicates higher surface moisture content and generally better ecological quality; conversely, a lower WET value indicates drier surface and poorer ecological quality.

[0092] S2.3: Calculation of the Improved Dryness Index NDBSI

[0093] Dryness index is a core indicator reflecting the degree of surface exposure and drying status, including built-up land, bare soil, and withered vegetation. The traditional dryness index NDBSI uses the calculation method of (SI + IBI) / 2, where SI is the bare soil index and IBI is the building index. As mentioned in the background section, the traditional NDBSI suffers from the problem of mutual contamination between the SI and IBI spectra, leading to distortion of the dryness signal.

[0094] To address this issue, this invention improves upon the traditional NDBSI by introducing the Visible Contrast Index (EI) to construct a three-dimensional evidence dryness index. The calculation formula is as follows:

[0095]

[0096] The calculation methods for each sub-index are as follows:

[0097] (1) Bare Soil Index (SI)

[0098] The SI (Soil Exposure Index) is primarily constructed based on the spectral differences between the shortwave infrared and red light bands, effectively reflecting the degree of bare soil exposure. The calculation formula is as follows:

[0099]

[0100] The SI value ranges from [-1, 1]. The SI value is usually higher in bare soil areas and lower in vegetation and water areas.

[0101] (2) Building Index (IBI)

[0102] IBI effectively extracts information about artificial building land use by combining the ratios of multiple bands. The calculation formula is as follows:

[0103]

[0104] The IBI value ranges from [-1, 1]. The IBI value is usually higher for built-up areas and lower for vegetated and water areas.

[0105] (3) Visible light contrast index EI

[0106] EI is a newly added dryness index in this invention, constructed based on the spectral differences between the red and blue light bands, and its calculation formula is as follows:

[0107]

[0108] The EI value ranges from -1 to 1. Bright bare soil and artificially bright surfaces (such as cement pavements and white building roofs) have significantly higher reflectance in the red light band than in the blue light band, resulting in higher EI values. Conversely, vegetation, water bodies, and dark soils show less difference in reflectance between the red and blue light bands, leading to lower EI values. Introducing EI enhances the ability of the dryness index to distinguish between bright bare soil and artificially bright surfaces, effectively compensating for the shortcomings of traditional NDBSI which relies solely on SI and IBI, and reducing the distortion of the dryness signal caused by spectral contamination.

[0109] The improved NDBSI integrates the sensitivity of SI to bare soil, IBI to built-up land, and EI to bright bare soil and artificial bright surfaces, thus reflecting the degree of surface drying more comprehensively and accurately. A higher NDBSI value indicates a drier surface, greater exposure, and poorer ecological quality; conversely, a lower NDBSI value indicates better vegetation cover, higher moisture content, and better ecological quality.

[0110] S2.4: Calculation of the Heat Index (LST)

[0111] Heat is a core indicator reflecting the state of the Earth's surface thermal environment, usually expressed as surface temperature (LST). Surface temperature not only directly affects plant growth, development, and transpiration, but is also closely related to ecological processes such as the urban heat island effect and regional climate regulation.

[0112] The surface temperature product used in this invention has undergone atmospheric and surface emissivity correction using a split-window algorithm, enabling the direct output of high-precision surface temperature data. The calculation formula is as follows:

[0113]

[0114] In the formula, The value of LST is the brightness temperature (K) obtained from the thermal infrared band inversion, and LST is the surface temperature (°C). The higher the LST value, the worse the surface thermal environment and the worse the ecological quality; conversely, the lower the LST value, the better the surface thermal environment and the better the ecological quality.

[0115] S2.5: Calculation of the Tillage Disturbance Index (NDTI)

[0116] Farmland disturbance is a unique ecological degradation pressure in agricultural counties. Agricultural activities such as tilling, sowing, and harvesting lead to reduced surface vegetation cover, exposed topsoil, and damage to soil structure, thereby causing a decline in ecological quality. The traditional RSEI index system lacks an independent indicator specifically characterizing farmland disturbance, resulting in the inability of agricultural ecological quality assessment results to accurately reflect the impact of farming activities.

[0117] To address this issue, this invention introduces the Normalized Tillage Disturbance Index (NDTI) as a fifth-dimensional ecological indicator, specifically designed to capture tillage disturbance information of cultivated land. The NDTI is constructed based on the spectral differences between the SWIR1 and SWIR2 bands, and its calculation formula is as follows:

[0118]

[0119] In the formula, and The reflectance values ​​are SWIR1 and SWIR2 bands, respectively. The NDTI value ranges from [-1, 1]. Bare farmland after tilling has significantly higher NDTI values ​​in the SWIR1 band than in the SWIR2 band due to its dry topsoil and low organic matter content. Conversely, vegetated farmland, water bodies, and built-up areas show smaller differences in SWIR1 and SWIR2 band reflectance, resulting in lower NDTI values. Introducing NDTI effectively distinguishes between drying caused by tillage disturbance and drying caused by built-up area expansion, significantly improving the relevance and explanatory power of agricultural county-level ecological quality assessment.

[0120] S3: Spatial Point Set Construction

[0121] Geographically weighted principal component analysis (GWPCA) is a computationally intensive spatial statistical method. If we directly perform pixel-by-pixel calculations on all pixels of an image, the computational load will be enormous, especially for large images at the county scale, which may lead to excessively long calculation times or even failure to complete the task.

[0122] To address this issue, this invention employs a regular grid sparse sampling method, which significantly reduces the number of samples involved in GWPCA calculations while ensuring uniform spatial coverage and evaluation accuracy. The specific steps for constructing the spatial point set are as follows:

[0123] Determine the sampling step size: Use Δ=30m as the sampling step size to generate a regular grid within the vector boundary of the study area. This step size is consistent with the image spatial resolution, ensuring a sufficiently high spatial density of sampling points and avoiding the omission of important ecological spatial variation information.

[0124] Selecting valid sampling points: Only the center points of the grids falling within the effective surface mask are retained as sampling points, while grid points covered by masks such as clouds, cloud shadows, and water bodies are removed, forming a sampling point set. (i=1,2,...,n), where n is the number of valid sampling points.

[0125] Coordinate transformation and index extraction: Extract the original geographic coordinates (longitude and latitude) of each sampling point and convert them into projected coordinates consistent with the preprocessed image; then extract the values ​​of the five ecological indicators (kNDVI, WET, NDBSI, LST, NDTI) calculated in step S2 to the corresponding sampling points to form a multidimensional dataset with coordinates. Each row of this dataset represents a sampling point, and each column represents an ecological indicator or coordinate information.

[0126] S4: Ecological Indicator Standardization and Geographically Weighted Principal Component Analysis

[0127] This step is the core of the invention. By standardizing and unifying the direction of ecological indicators, the problems of dimensional differences and inconsistencies in direction are eliminated. Then, the geographically weighted principal component analysis method is used to calculate the first principal component score of each sampling point, realizing adaptive local integration of multidimensional ecological information. Finally, the PC1 directional anchoring operation is used to solve the problem of arbitrary PCA principal axis direction and ensure the spatial consistency of evaluation results.

[0128] S4.1: Unification and Standardization of Indicator Direction

[0129] Since different ecological indicators have different dimensions, value ranges and ecological significance, it is necessary to unify and standardize the indicators before conducting principal component analysis in order to eliminate the influence of dimensional differences on the analysis results.

[0130] First, the direction of the indicators was unified: based on the correlation between each indicator and ecological quality, the indicators were divided into positive and negative indicators. Among them, kNDVI and WET are positive indicators, and the larger their values ​​are, the better the ecological quality. Therefore, their original values ​​are kept unchanged. NDBSI, LST, and NDTI are negative indicators, and the larger their values ​​are, the worse the ecological quality. Therefore, they are multiplied by -1 to make them positive, so that the direction of all indicators is unified as "the larger the value, the better the ecological quality".

[0131] Then Z-score standardization is performed: Z-score standardization can scale all indicators to the same scale with a mean of 0 and a standard deviation of 1. The calculation formula is:

[0132]

[0133] In the formula, This is the original value of the i-th indicator (the direction has been unified). and These represent the mean and standard deviation of the index across all sampling points, respectively. These are the standardized index values. The standardized indexes have a mean of 0 and a standard deviation of 1, eliminating dimensional differences and allowing them to be directly used for principal component analysis. It is particularly important to note that the positive conversion of negative indices must be completed before standardization; otherwise, the standardization results will be directionally incorrect, affecting the accuracy of subsequent principal component analysis.

[0134] S4.2: Geographically Weighted Principal Component Analysis (GWPCA)

[0135] Traditional global principal component analysis (PCA) calculates a fixed set of loading coefficients based on the covariance matrix of all samples within the study area, which fails to reflect the spatial non-stationarity of ecological driving factors. Geographically weighted principal component analysis (GWPCA), however, introduces a spatial weighting function to construct a local neighborhood centered on each sample point. Principal component loadings are calculated separately within each local neighborhood, allowing for differentiated contribution structures of indicators at different locations, thus accurately characterizing the spatial variability of key ecological factors.

[0136] The specific calculation steps for GWPCA are as follows:

[0137] Determining the Adaptive Bandwidth: Bandwidth is a key parameter in GWPCA, determining the size of the local neighborhood. This invention employs an adaptive nearest neighbor method to determine the bandwidth, meaning that the local neighborhood of each sampling point consists of its 200 spatially nearest points. The advantage of adaptive bandwidth is that in densely populated regions, the spatial range of the local neighborhood is smaller; in sparsely populated regions, the spatial range of the local neighborhood is larger. This ensures that each local neighborhood has sufficient samples for estimating the covariance matrix, while avoiding excessively large bandwidth that could lead to excessive spatial smoothing.

[0138] Calculate spatial weights: for each sampling point Each sample point within its local neighborhood Each is assigned a spatial weight, the magnitude of which is related to... arrive The spatial weights are inversely proportional to the spatial distance. This invention uses the bisquare kernel function to calculate the spatial weights, and the calculation formula is as follows:

[0139] In the formula, is the spatial weight coefficient, representing For the center point Influence weight; Sampling points and The Euclidean distance between them Sampling points The local bandwidth (i.e., the distance to the 200th nearest neighbor). When ≤ The time weight is calculated according to the above formula, when The time weight is 0. The bisquare kernel function is a continuously decreasing kernel function that gives greater weight to sample points closer to the center and smaller weight to sample points farther from the center. Sample points exceeding the bandwidth have a weight of 0, which conforms to the spatial decay law of geographical phenomena.

[0140] Constructing locally weighted decentralized data: for each sampling point By using spatial weights to weight standardized indicator data within a local neighborhood and then decentralizing the data, a locally weighted decentralized data matrix is ​​obtained. .

[0141] Calculating the local covariance matrix: Based on the locally weighted decentralized data matrix Calculate the local covariance matrix The calculation formula is:

[0142]

[0143] In the formula, This represents the number of samples within the local neighborhood (200 in this example). This is a local space weight diagonal matrix, whose diagonal elements are... .

[0144] Eigenvalue decomposition and local load vector calculation: for the local covariance matrix Eigenvalue decomposition is performed to obtain eigenvalues ​​and their corresponding eigenvectors. The eigenvalues ​​are then arranged in descending order; the eigenvector corresponding to the largest eigenvalue is the local loading vector of the first principal component. It reflects the sampling point Within a local neighborhood, the magnitude of the contribution of each ecological indicator to ecological quality.

[0145] Calculate PC1 score: sampling points The first principal component PC1 score is obtained by the dot product of the standardized index vector and the local load vector, and the calculation formula is as follows:

[0146]

[0147] In the formula, Sampling points Standardized index vector, This is the PC1 score for this sampling point. The PC1 score integrates information from five ecological indicators within the local neighborhood and can reflect the overall ecological quality of that location.

[0148] S4.3: PC1 directional anchoring

[0149] Principal component analysis (PCA) has an inherent problem: the direction of the principal axis is arbitrary, meaning that multiplying an eigenvector by -1 still results in an eigenvector corresponding to the same eigenvalue. In global PCA, the principal axis direction is typically determined by examining the correlation between the first principal component and the greenness index; if the correlation coefficient is negative, PC1 is multiplied by -1. However, in GWPCA, the local loading vector of each sampling point is calculated independently. This can lead to situations where the PC1 direction of some sampling points is opposite to that of other sampling points, resulting in spatial inconsistencies in the evaluation results. Specifically, two adjacent sampling points might exhibit anomalies where one has a very high PC1 value and the other a very low PC1 value.

[0150] To address this issue, this invention designs a PC1 orientation anchoring mechanism, using kNDVI as a reference to unify the PC1 orientation of all sampling points. The specific operation is as follows:

[0151] Calculate the Pearson correlation coefficient r between the PC1 sequence and the kNDVI sequence at all sampling points.

[0152] If r < 0, then multiply the PC1 score of all sampling points by -1; if r ≥ 0, then keep the PC1 score unchanged.

[0153] After anchoring, the PC1 score is positively correlated with kNDVI, meaning that a higher PC1 value indicates better ecological quality and a lower PC1 value indicates worse ecological quality, thus ensuring the spatial consistency and clarity of the physical meaning of the evaluation results.

[0154] S5: Post-processing and Result Output

[0155] This step converts the PC1 scores of discrete sampling points into continuous ecological quality index raster images, performs normalization and classification processing, and finally outputs a spatial distribution map of ecological quality.

[0156] S5.1: Spatial Interpolation

[0157] Since the sampling points are discrete, the PC1 score needs to be converted into a continuous grid surface using spatial interpolation. This invention employs Inverse Distance Weighted (IDW) interpolation, which is based on the spatial interpolation principle of "the closer the distance, the greater the weight." This method is simple to calculate, produces smooth interpolation results, and effectively reflects the spatial continuity of ecological quality. The formula for IDW interpolation is as follows:

[0158]

[0159] In the formula, Points to be interpolated PC1 score, Sampling points PC1 score, Sampling points Relative to the interpolation point The weights are denoted by n, where n is the number of neighboring sampling points involved in the interpolation. The calculation formula is:

[0160]

[0161] In the formula, Points to be interpolated to sampling point The distance is p, where p is the power exponent. In this embodiment, the parameters for IDW interpolation are set as follows:

[0162] Power exponent p=2: This is the most commonly used power exponent in IDW interpolation. It can balance the rate of distance decay, so that the interpolation result is neither too smooth nor too fragmented.

[0163] Maximum search radius maxdist=3000m: Limits the maximum distance between sampling points participating in interpolation, avoiding unreasonable influence of sampling points that are too far apart on the interpolation results.

[0164] Maximum number of nearest neighbors nmax=12: Each point to be interpolated uses at most 12 of its nearest sampling points for interpolation, which improves computational efficiency while ensuring interpolation accuracy.

[0165] By using IDW interpolation, the PC1 scores of discrete sampling points are converted into a continuous raster image with the same spatial resolution as the original remote sensing image.

[0166] S5.2: Min-Max Normalization

[0167] Since the range of PC1 scores varies across different study regions and years, to facilitate cross-sectional comparisons, the interpolated PC1 raster needs to be Min-Max normalized, linearly stretching its value range to the [0,1] interval. The normalization formula is as follows:

[0168]

[0169] In the formula, and These are the minimum and maximum values ​​of the PC1 raster in the current study area, respectively. This is the normalized ecological quality index. After normalization, the closer the RSEI value is to 1, the better the ecological quality; the closer it is to 0, the worse the ecological quality.

[0170] S5.3: Ecological Quality Classification and Output

[0171] To more intuitively illustrate the spatial distribution of ecological quality, based on the RSEI value range [0,1], ecological quality is divided into five levels with an interval of 0.2:

[0172] Difference: RSEI∈[0,0.2)

[0173] Poor: RSEI∈[0.2,0.4)

[0174] Medium: RSEI∈[0.4,0.6)

[0175] Good: RSEI∈[0.6,0.8)

[0176] Excellent: RSEI∈[0.8,1]

[0177] Finally, the classified raster is cropped using the vector boundary of the study area to obtain the final spatial distribution map of ecological quality, and output as commonly used raster data formats such as TIFF, which facilitates subsequent spatial analysis and application.

[0178] Example:

[0179] This embodiment selects two typical agricultural counties as research areas to conduct ecological quality assessment application verification. Research areas 1 and 2 are both located in the Songliao Plain in Northeast China, which is an important commodity grain production base in my country. The counties are mainly composed of arable land, but also have various land use types such as forest land, grassland, and construction land. The ecological driving factors have obvious spatial heterogeneity, making them suitable for verifying the applicability and effectiveness of the method of this invention in agricultural counties.

[0180] This embodiment uses surface reflectance products from the 2022 growing season (July-August) as the data source and conducts an ecological quality assessment according to steps S1 to S5 above. The specific implementation details are as follows:

[0181] Remote sensing image preprocessing: Products covering the study area were retrieved, and DN values ​​were converted into surface reflectance and brightness temperature using official calibration coefficients; clouds, cloud shadows, and fill pixels were removed using QA_PIXEL band masking; water bodies were extracted using the MNDWI index and masked; the images were converted to WGS84UTM51N projection and cropped using the vector boundary of the study area to obtain standardized multispectral images and surface temperature images.

[0182] Calculation of the improved ecological index set: According to the formula in step S2, calculate the five ecological indices kNDVI, WET, improved NDBSI, LST and NDTI respectively, and obtain the raster images of each index.

[0183] Spatial point set construction: A regular grid is generated with a step size of 30m, valid sampling points are selected, and the coordinates of the sampling points and the values ​​of five ecological indicators are extracted to form a multidimensional dataset with coordinates.

[0184] GWPCA Integration and Anchoring: Adaptive bandwidth bw=200NN, bisquare kernel function; directional unification and Z-score standardization of the five indicators; calculation of local PC1 scores for each sampling point; anchoring PC1 direction with kNDVI as reference to ensure that the higher the PC1 value, the better the ecological quality.

[0185] Post-processing and output: IDW interpolation was used to convert discrete PC1 scores into 30m resolution raster images; Min-Max normalization was performed to obtain the RSEI index; ecological quality was divided into five levels, and a spatial distribution map of ecological quality was output.

[0186] By comparing the evaluation results of the method of this invention with those of the traditional RSEI method, it can be found that: the ecological quality spatial distribution boundary generated by the method of this invention better conforms to the actual land cover pattern, and can clearly identify the differences in ecological quality among different land cover types such as cultivated land, forest land, and construction land; the transition zone of ecological quality is more natural and is not smoothed out by a single main axis across the county; the local ecological coupling structure is clearer, and can accurately reflect the spatial impact of different ecological driving factors such as cultivation disturbance and urban expansion. In contrast, the evaluation results of the traditional RSEI method have obvious problems of spatial smoothing transition, the boundaries of different land cover types are blurred, local ecological differences are masked, and it cannot accurately reflect the impact of cultivation disturbance in agricultural areas.

[0187] In summary, this invention effectively addresses a series of inherent defects in the traditional RSEI method by constructing an improved set of ecological indicators and introducing a geographically weighted principal component analysis method. It significantly improves the accuracy, spatial consistency, and ecological interpretability of ecological quality assessment at the county scale, providing strong technical support for county-level ecological environmental protection, ecological restoration planning, and land space management.

[0188] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and scope of the present invention should be included within the protection scope of the present invention.

Claims

1. An ecological quality assessment method based on an improved index set, characterized in that: Includes the following steps: S1: Acquire multispectral remote sensing images and perform radiometric calibration, atmospheric correction and geometric correction preprocessing; S2: Based on the aforementioned multispectral remote sensing imagery, calculate an ecological indicator dataset comprising five indicators: greenness, humidity, dryness, heat, and tillage index; wherein, the greenness indicator uses kNDVI, and the calculation formula is as follows: In the formula, For near-infrared reflectivity, Reflectivity in the red light band; The humidity index uses the tasseled cap transformation (WET), and the calculation formula is as follows: In the formula, These are the blue, green, red, near-infrared, SWIR1, and SWIR2 bands, respectively. The dryness index used is the modified normalized dryness index (NDBSI), calculated using the following formula: The formula for calculating the bare soil index SI is as follows: The formula for calculating the Building Index (IBI) is as follows: The formula for calculating the Visible Contrast Index (EI) is: The thermal index uses land surface temperature (LST), and the calculation formula is as follows: In the formula, Brightness temperature obtained by inversion of remote sensing images in the thermal infrared band; The tillage index uses the Normalized Tillage Disturbance Index (NDTI), calculated using the following formula: In the formula, and The reflectivity of SWIR1 band and SWIR2 band are respectively. This leads to the construction of a five-dimensional ecological indicator set, including kNDVI, WET, NDBSI, LST, and NDTI. S3: Constructing a spatial point set: The ecological indicator dataset is sampled according to a regular grid, and the planar coordinates of each sampling point are extracted. The longitude and latitude are converted to projected coordinates to form a multidimensional dataset with coordinates. S4: Standardize the five-dimensional ecological index set constructed in step S2: keep the positive indices kNDVI and WET at their original values, and multiply the negative indices NDBSI, LST, and NDTI by -1 to make them positive; use Z-score standardization to scale each index to a mean of 0 and a standard deviation of 1; use the geographic weighted principal component analysis (GWPCA) model to construct local neighborhoods with adaptive bandwidth, calculate the local covariance matrix and loading vector, and obtain the first principal component PC1 score for each sampling point; ensure that PC1 and kNDVI are aligned through anchoring operations; S5: Interpolate the PC1 score to a raster image, and use Min-Max normalization to scale the value range to [0,1] to obtain the ecological quality index raster image; according to the ecological index value range [0,1], divide it into five levels of poor, poor, moderate, good and excellent with an interval of 0.2, and output the ecological quality spatial distribution map.

2. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S2, kNDVI undergoes a nonlinear transformation of NDVI using the tanh hyperbolic tangent function.

3. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S4, the formula for calculating the Z-score standardization is as follows: In the formula, Let be the original value of the i-th indicator. and These represent the mean and standard deviation of the indicator across all sampling points, respectively; the positiveification of negative indicators is performed before standardization.

4. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S4, the adaptive bandwidth of GWPCA adopts the nearest neighbor method, and the local neighborhood of each sampling point consists of its 200 nearest points in space; the local load vector is obtained by eigenvalue decomposition of the local covariance matrix, and the first principal component PC1 corresponds to the direction of the largest eigenvalue.

5. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S4, the determination condition for the PC1 anchoring operation is: calculate the correlation coefficient between the PC1 sequence and the kNDVI sequence; if the correlation coefficient is negative, then multiply PC1 by -1. Ensure that the final PC1 value is positively correlated with kNDVI.

6. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S5, the interpolation from PC1 to the grid uses the inverse distance weighted IDW method, with a power exponent of 2, a search radius of 3000m, and a maximum number of neighboring points of 12.

7. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S1, the spatial resolution of the multispectral remote sensing image should be better than or equal to 30m.

8. The ecological quality assessment method based on an improved index set according to claim 1, characterized in that: In step S1, when using surface reflectance data from multispectral satellite remote sensing imagery, the formula for converting the original DN value into surface reflectance is as follows: In the formula, Let be the surface reflectance of the i-th optical band. This is the original DN value for this band; the formula for converting the original DN value of the thermal infrared band to brightness temperature is: In the formula, For brightness temperature, This represents the original DN value for the thermal infrared band.

Images