A Parameter Optimized Method for InVEST Model in Sub-Pixel Scale Integrating Machine Learning Algorithm and Vegetation–Impervious Surface–Soil Model

Abstract

The InVEST model, with its ability to perform spatial visualization and quantification, is an important tool for mapping ecosystem services. However, the spatial accuracy and simulating performance of the model are deeply influenced by the land use parameter, which often relies on the accuracy of land use/cover data. To address this issue, we propose a novel method for optimizing the land use parameter of the InVEST model based on the vegetation–impervious surface–soil (V–I–S) model and a machine learning algorithm. The optimized model is called Sub-InVEST, and it improves the performance of assessing ecosystem services on a sub-pixel scale. The conceptual steps are (i) extracting the V–I–S fraction of remote sensing images based on the spectral unmixing method; (ii) determining the mapping relationship of the V–I–S fraction between land use/cover type using a machine learning algorithm and field observation data; (iii) inputting the V–I–S fraction into the original model instead of the land use/cover parameter of the InVEST model. To evaluate the performance and spatial accuracy of the Sub-InVEST model, we employed the habitat quality module of InVEST and multi-source remote sensing data, which were applied to acquire Sub-InVEST and estimate the habitat quality of central Guangzhou city from 2000 to 2020 with the help of the LSMA and ISODATA methods. The experimental results showed that the Sub-InVEST model is robust in assessing ecosystem services in sets of complex ground scenes. The spatial distribution of the habitat quality of both models revealed a consistent increasing trend from the southwest to the northeast. Meanwhile, linear regression analyses observed a robust correlation and consistent linear trends, with R² values of 0.41, 0.35, 0.42, 0.39, and 0.47 for the years 2000, 2005, 2010, 2015, and 2020, respectively. Compared with the original model, Sub-InVEST had a more favorable performance in estimating habitat quality in central Guangzhou. The spatial depictions and numerical distribution of the results of the Sub-InVSET model manifest greater detail and better concordance with remote sensing imagery and show a more seamless density curve and a substantially enhanced probability distribution across interval ranges.

1. Introduction

Ecosystem services (ESs) refer to the material and non-material benefits that directly or indirectly arise from socio-ecological systems [1,2]. ESs are crucial for reflecting the multitude of benefits ecosystems offer to humanity and their capacity to deliver these services [3,4,5]. Assessing and mapping ESs has a crucial significance for increasing our awareness of how ES provision affects the interactions of socio-ecological systems, as well as for the application of sustainable ecological management practices [6,7]. However, due to the complex linkages between ecosystem processes, functions, and human well-being, this remains a significant challenge. There is a need for the development and improvement of feasible methods and approaches for assessing and mapping ESs [8].

To match the demand for mapping and assessing ESs, ecosystem services modeling tools were developed and applied on local or global scales, and they abruptly acquired increasing attention [9]. Commonly utilized modeling tools include the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) suite of models [10], Artificial Intelligence for Ecosystem Services (ARIES) [11], Social Values for Ecosystem Services (SolVES) [12], the Intelligent Urban Ecosystem Management System (IUEMS) [13], etc. Notably, InVEST models have seen widespread and popular application due to their multifunctionality in assessing various ESs, user-friendly operation, and concise data sources [14]. Numerous studies have validated their performance and feasibility in assessing and mapping Ess; examples include a study of the North–East Green Belt region that quantified carbon storage, potential pollination, flood mitigation, and several ESs based on the InVEST model [15]; a study on quantifying terrestrial and aquatic habitat quality in the Llobregat River basin [16]; and a study on evaluating and predicting carbon storage in the Sariska Tiger Reserve [17].

The land use parameter is the most crucial component of the InVEST model. It serves as the foundation for estimating and simulating ESs at the landscape scale and necessitates the adjustment and setting of other parameters in the InVEST model according to different land use types [18]. The data sources for land use parameters in the InVEST model include field data and remotely sensed land use classification data. On the one hand, national land surveys are typically conducted once every ten years, yielding highly accurate land use mapping data. However, the long data-updating cycle and substantial consumption of manpower and material resources hinder continuous observation [19]. On the other hand, remotely sensed land use classification data can quickly and broadly capture regional land use status but are constrained by the spatiotemporal resolution of the remote sensing data and the land use classification methods [20]. This implies that the timeliness and classification scale of the input LULC data directly influence the accuracy and usability of the InVEST model’s results. Therefore, optimizing the land use parameters of the InVEST model remains an area that requires innovation, especially with regard to improving the model’s data inputs.

In this study, we proposed an optimization method for the land use parameter of the InVEST model which employs the vegetation–impervious surface–soil (V–I–S) model and a machine learning algorithm. To capture the component of land surface information at the sub-pixel scale, we introduced the V–I–S model to obtain V–I–S fractions from remote sensing images. The mapping relationship between LULC types and V–I–S fractions was then defined using a machine learning algorithm. Subsequently, the land use data within the InVEST model were substituted with V–I–S fractions for the purpose of developing the optimized model, namely the Sub-InVEST model. A comparative assessment was conducted by evaluating Sub-InVEST against the traditional InVEST model and remote sensing imagery. The main purposes of this study were to (i) propose a land use parameter optimization method for the InVEST model and to (ii) evaluate the spatial consistency of and numerical variation in the InVEST model and the optimized model.

2. Related Work

2.1. Vegetation–Impervious Surface–Soil Model

The V–I–S model is a classical conceptual framework proposed by Ridd [21] to express urban surface component information and urban morphology characteristics. The V–I–S model considers the base material of the urban surface to be a linear combination of green vegetation (V), impervious surfaces (I), and soils (S), excluding water bodies. This model abstracts urban surface characterization through V–I–S fractions, providing an objective and quantitative method for determining the biophysical parameters of urban ecosystems [22].

The V–I–S model has been widely applied to the quantitative characterization of land surfaces, with its feasibility and applicability validated in diverse regions globally. For example, Weng and Lu [23] obtained V–I–S fractions in Indianapolis, USA, using Landsat TM/ETM+ imagery and utilized these fractions to quantitatively characterize the composition of urban land use. Madhavan et al. [24] used the V–I–S model to quantitatively reveal land class trends in the Bangkok region, Thailand. Aina et al. [25] analyzed land use changes and spatial patterns in the Riyadh region, Saudi Arabia, across different years by employing the V–I–S model alongside the landscape index and urban expansion index. These studies demonstrate that V–I–S fractions effectively describe land use compositional information at a sub-pixel scale, providing a deeper and more precise expression of the material composition within a given image.

2.2. Machine Learning Algorithm for Ecosystem Service Assessment

In recent years, machine learning (ML) algorithms have been increasingly applied to ecosystem service assessments. By integrating large socio-ecological datasets with linear or nonlinear algorithms, ML provides convenient tools and realistic models for mapping ESs [9]. ML is crucial for inferring relevant indicators of ESs, particularly when data sources are insufficient, and it paves the way for understanding interactions or variations among multiple ESs [26,27].

Several studies have employed ML algorithms in ecosystem service assessments, demonstrating significant improvements in the modeling and mapping of these services. For instance, Sun et al. [28] investigated the interactions of multiple ESs and various social-ecological factors using a Bayesian belief network. Sanderman et al. [29] developed an ML model to map the distribution of carbon density in mangrove forests. Kundu et al. [30] assessed the ecosystem services value of wetlands by extracting spatial patterns of wetlands using a support vector machine. Alqadhi et al. [4] used an ensemble machine learning algorithm to develop a landslide prediction model, which was then used to assess the potential loss of ecosystem service value. These studies indicate that ML algorithms can be effectively applied to determine the mapping relationship between V–I–S fractions and LULC types, thereby optimizing the InVEST model for mapping ESs.

3. Methodology

3.1. Conceptual Steps for Model Parameter Optimization

The V–I–S model has a clear physical meaning, which indicates the proportion of the three endmembers in each mixed pixel. Based on the definition of the V–I–S model, any mixed pixel of optical remote sensing imagery can be characterized by V–I–S fraction combination. When the V–I–S fraction combination of one mixed pixel is identical or similar to that of another mixed pixel, the component of the land surface in these pixels should be defined in the numerical interval of the V–I–S fraction combination (Figure 1).

Therefore, an optimized method was developed at a sub-pixel scale, centered on the V–I–S model, with a key focus on utilizing V–I–S fractions to enhance the InVEST model (Figure 2). The first step involves constructing a mapping relationship between LULC types and V–I–S fractions. V–I–S fraction data were obtained from remote sensing imagery, and V–I–S fraction rule sets with similar characteristics were categorized. Then, the mapping relationship between LULC types and V–I–S fractions was defined. In the second step, the optimized model was proposed based on the mapping relationship between V–I–S fractions and LULC types. V–I–S fraction data are used to replace the land use data required by the InVEST model, resulting in the optimized model Sub-InVEST. In the third step, a comparative assessment of Sub-InVEST and InVEST was conducted. Evaluation indexes were selected from visual interpretation and numerical variation, through which the model performance was assessed.

3.2. Methods for Operationalizing Model Optimization

Based on the conceptual steps, we adopted the LSMA method and ISODATA clustering algorithm to execute the model-optimization process. The main operationalizing steps involved were as follows: (i) we applied LSMA to extract V–I–S fraction data from remote sensing datasets, (ii) we established the mapping relationship between LULC and the V–I–S fraction based on the ISODATA, and (iii) we acquired the Sub-InVEST model based on the V–I–S fraction and evaluated the accuracy assessment.

3.2.1. Linear Spectral Mixture Analysis (LSMA) Method for Extracting V–I–S Fractions

The linear spectral mixture analysis (LSMA) method can estimate the endmember abundance at the sub-pixel scale [31]. In this study, we applied LSMA and the V–I–S model to extract the fraction of vegetation (V), impervious surface (including high albedo: H and low albedo: L), and soil (S) in the mixed pixels of Landsat images. Among them, the fraction of impervious surface (I) is the sum of high-albedo impervious surface (H) and low-albedo impervious surface (L). The Equation of LSMA is as follows:

$$R_i=\sum _{k=1}^n{f}_k{R}_{i,k}+e{R}_i$$

(1)

where R_i is the reflectance of band i (i = 1, 2, 3, …, m) in mixed pixels, f_k is the fraction of endmember k (k = 1, 2, 3,…, n), R_i,k is the reflectance of band i in the endmember k, and eR_i is the residual error of LSMA. In addition, this Equation must satisfy the following conditions: $\sum _{k=1}^n{f}_k=1$ and f_k ≥ 0.

3.2.2. Machine Learning Algorithm for Determining the Mapping Relationship of V–I–S Fraction and LULC

To determine the mapping relationship of the V–I–S fractions and land use/cover types, we assumed that the V–I–S fraction combination, X_i, of LULC type i and the V–I–S fraction rule sets Y_j showed the mapping relationship; the following conditions must be satisfied: (i) the V–I–S fraction combination of LULC type i of the sample point x∈X_i and the class j for the V–I–S fraction rule set Y_j are a continuous closed interval; (ii) the numerical range of the V–I–S fraction combination of LULC type i is a subset of the V–I–S fraction rule set, where $X_i\in Y_j$; (iii) the V–I–S fraction combination, X_i, of LULC type i corresponds to the uniquely determined V–I–S fraction rule set Y_j.

Considering the aforementioned issue, the following steps determined the mapping relationship of the V–I–S fraction and LULC: (i) establishing the V–I–S fraction rule sets based on a ML algorithm and (ii) taking a vote to determine the mapping relationship between the V–I–S fraction rule sets and LULC types based on sufficient field observation records.

Among them, the clustering algorithm has the advantage of segmenting the similar or different subsets of each sample point in their feature space based on Euclidean or correlation distance, which can identify the similar or same characteristics or expressions for different clusters [32]. In the first step, we selected the iterative self-organizing data analysis techniques algorithm (ISODATA) to establish the V–I–S fraction rule sets, which is an unsupervised learning algorithm that can automatically estimate different clusters by separating the feature space of V, I, and S, using Euclidean distance to measure the similarity of the V, I, and S fractions. In the next step, we first recorded the range of the V–I–S fraction of all field observations and counted the frequency of each observation record falling into the V–I–S fraction rule set, which could determine the mapping relationship of the V–I–S fraction and LULC according to the most frequent V–I–S fraction rule set. The details of the ISODATA used to determine the mapping relationship of V–I–S fraction and LULC types are listed below.

Step 1: Input the V–I–S fraction dataset and pre-setting the parameters. The pre-set parameters include the maximum number of cluster centers N_max, the minimum number of samples that can be a cluster ${\theta }_N$, the standard deviation of the samples from their cluster ${\theta }_S$, the minimum distance between two center clusters ${\theta }_C$, the maximum number of cluster centers that can be merged L, and the maximum number of iterations I.

Step 2: The V–I–S fraction samples are divided into the nearest cluster based on the minimum distance, and then the average distance between each cluster center, ${\overline{D}}_j$, and the cluster centers, $Z_j$, are corrected and calculated. It is assumed that if $D_j=min\left\{‖x-z_i‖,i=1,2,\dots ,N_c\right\}$, then $x\in S_j$.

$$Z_j={\displaystyle \frac{1}{N_{j\in S_j}}}\sum X,j=1,2,\dots ,N_c$$

(2)

$${\overline{D}}_j={\displaystyle \frac{1}{N_{j\in S_j}}}\sum ‖X-Z_j‖,j=1,2,\dots ,N_c$$

(3)

where $N_c$ is the number of initial cluster centers; $N_c$ can or cannot be equal to $N_{max}$, and $N_j$ is the number of samples of $S_j$.

Step 3: According to the pre-set parameters, the classified clusters are split and merged, and a new cluster center is recalculated. The inter-cluster distances between all of the cluster centers can be calculated by using the following Equation:

$$D_{ij}=‖Z_i-Z_j‖,i=1,2,\dots ,N_c-1;j=i+1,\dots ,N_c$$

(4)

Step 4: Iterative calculation is carried out repeatedly to judge whether the clustering results meet the requirements until the clustering results that meet the requirements of separability are obtained.

3.2.3. Evaluation Metrics for an Accuracy Comparison Assessment

We assumed that if the Sub-InVEST model shows similar or the same simulation results as the original model, the optimized method is valid and feasible for assessing ESs. To validate and quantify the performance of the Sub-InVEST model, we adopted two main evaluation metrics: numerical difference and visual interpretation. For numerical difference assessment, we employed the linear fitting method to quantify the linear trend in the Sub-InVEST model compared to the original model. This allowed us to measure the degree of similarity between the simulation results of the two models. In addition, visual interpretation was utilized as an evaluation metric. This involved dividing the study area into several sample regions and comparing the model results with remote sensing data. By assessing the degree of accuracy between the ground truth scenes and simulation results, we could determine the effectiveness of the Sub-InVEST model in accurately depicting ESs.

These evaluation metrics collectively provided insights into the performance and efficacy of the Sub-InVEST model, helping to ascertain its validity and feasibility for assessing ESs.

4. Experiments

4.1. Experimental Area and Data Preparation

In this study, we took central Guangzhou as the study area to examine the performance of the Sub-InVEST model, which was designed as a comparative study of the optimized model and the original model. The location of central Guangzhou is shown in Figure 3, which includes the Districts of Tianhe, Haizhu, Yuexiu, Liwan, Huangpu, and Baiyun, with a total area of 1559.59 km². Meanwhile, we selected the habitat quality model, which is a widely used application module for InVEST, to compare the sequence generation results under the Sub-InVEST and InVEST models. To determine the degree of accuracy, the ground truth scenes, and the simulation results, we combined the remote sensing data to participate in a performance comparison.

The InVEST habitat quality model is a GIS-based model that can calculate the regional habitat quality index based on LULC maps, a raster layer of the threat factor. The parameters of the original model included a land use parameter, threat factors, the weight of each threat factor, the maximum distance of each threat factor, the sensitivity score of habitat types to threat factors, and the habitat suitability score. In this study, we took the V–I–S fraction generated by Landsat series images to replace the original LULC maps of the land use/cover parameter, which can be used to acquire the Sub-InVEST model to optimize the original InVEST model.

To examine and verify the performance of the Sub-InVEST model, the experiment data included remote sensing datasets, LULC datasets, and threat factors datasets.

(1)

Remote sensing datasets and preprocesses

The remote sensing datasets included Landsat datasets and Google Earth images. The Landsat datasets were used to extract the V–I–S fraction, and Google Earth images were applied to verify the accuracy of the V–I–S fraction and the performance of Sub-InVEST. The Landsat series datasets were Landsat Collection 2 Level-2 Tier 1 Surface reflectance products, which included Landsat 5 TM (Thematic Mapper) data from 1999, 2004, and 2009, and Landsat 8 Operational Land Imager (OLI) data from 2015 to 2020 (the cloud cover of each image is less than 5%); the specific descriptions are shown in

Table 1. The descriptions of the Landsat series satellite dataset.

Year	Datasets	Acquired Time	Band	Path/Row
2000	Landsat 5 TM	9 December 1999	Band 1–5, Band 7	122/044
2005	Landsat 5 TM	21 January 2004	Band 1–5, Band 7	122/044
2010	Landsat 5 TM	2 January 2009	Band 1–5, Band 7	122/044
2015	Landsat 8 OLI	18 October 2015	Band 1–7	122/044
2020	Landsat 8 OLI	18 February 2020	Band 1–7	122/044

. The datasets were derived from the US geological survey (https://earthexplorer.usgs.gov/, accessed on 14 July 2022), with a 30 m spatial resolution. Landsat Collection 2 Level-2 Tier 1 Surface reflectance products have completed data-preprocessing processes, such as geometric correction, radiometric calibration, and atmospheric correction, and their DN values range from 1 to 65455, which can acquire surface reflectance by calculating 0.0000275 × DN − 0.2. All remote sensing images were transferred to Albers equal area projection, and the outlier values were removed. The Google Earth images were applied to extract the regions of V, I, and S, respectively, at the same time as the Landsat dataset, which compared the proportion of V, I, and S in each pixel with the V–I–S fraction of Landsat images. In addition, the Google Earth images were used to compare the sequence simulation results of the optimized the original models for the evaluation metrics of visual interpretation.

(2)

Land use/cover data and preprocesses

To match the spatial resolution of Landsat datasets, we applied land use/cover data generated by Landsat series satellite remote sensing datasets with the same or the closed data source of the V–I–S fraction. The land use/cover data were China land/use cover change (CNLULC) data, which were acquired from the Resources and Environmental Sciences Data Center, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn, accessed on 14 July 2022), with a fine spatial resolution of 30 m. The CNLULC data had a comparatively high comprehensive accuracy, exceeding 90%, which was generated from Landsat series images based on a supervised classification method. We selected the CNLULC data for the years 2000, 2005, 2010, 2015, and 2020 to calculate habitat quality by using the original InVEST model, which was used to compare the simulation results of the optimized model. In this study, the eight categories of LULC types were divided, including cropland, forest, grassland, shrubland, water bodies, wetland, built-up areas, and unused land.

(3)

Threat factors datasets and preprocesses

The threat factors mainly considered the threats of anthropogenic disturbance to each habitat type. We referred to the InVEST user’s guide [18] and the relevant literature [33,34] to determine artificial surfaces as the threat factors, which include the road network, built-up areas, unused land, cropland, residential areas, and industrial activities. Among them, the built-up areas and unused land were selected from CNLULC. The road network datasets, the residential areas, and industrial activities were obtained from OpenStreetMap (OSM) (https://download.geofabrik.de/, accessed on 20 January 2021), obtained from the attribution of the land use layer and road layer of OSM data.

The parameter settings of the threat factors included the weight of each threat factor, the maximum distance of each threat factor, and the sensitivity score of habitat types to threat factors. In this study, the threat factor parameters were given by the InVEST user’s guide [18], the relevant literature [16,33,35,36], and our previous study [37], which are shown in

Table 2. The parameter settings of the threat factors.

Threat Factors	Weight	Maximum Distance	Decay Type
Unused land	0.2	3	linear
Built-up areas	1	10	exponential
Cropland	0.68	8	linear
Railway	0.9	9	exponential
Trunk road	1	10	exponential
Primary road	1	8	linear
Secondary road	0.75	5	linear
Industrial activities	1	12	exponential
Residential area	0.5	5	exponential

and

Table 3. The sensitivity of the habitat types between threat factors.

Land Use/Cover Type	Threat Factors
	Cropland	Built-Up Areas	Unused Land	Railway	Trunk Road	Primary Road	Secondary Road	Industrial Activities	Residential Areas
Cropland	0	0.4	0.1	0.35	0.35	0.3	0.2	0.6	0.1
Forest	0.3	0.8	0.2	0.75	0.75	0.7	0.6	0.8	0.8
Grassland	0.35	0.6	0.1	0.7	0.7	0.5	0.35	0.7	0.6
Shrubland	0.35	0.6	0.1	0.7	0.7	0.5	0.35	0.7	0.6
Wetland	0.3	0.85	0.3	0.8	0.8	0.75	0.65	0.8	0.8
Water bodies	0.9	0.9	0.5	0.5	0.5	0.45	0.3	0.9	0.7
Built-up areas	0	0	0.3	0.6	0.6	0.5	0.5	0.2	0.1
Unused land	0	0.5	0	0.1	0.1	0.1	0.1	0.2	0.2

.

4.2. Habitat Quality Model Optimization Based on the V–I–S Model

In this study, a model optimization method introducing sub-pixel information was designed to replace the land use data of the original model with the V–I–S fraction and obtain the optimized Sub-InVEST model based on the V–I–S model and a ML algorithm. We selected the calculation of habitat quality as an example, and the optimization test experiment was carried out using the habitat quality model. Therefore, the specific calculation formula of the habitat quality model after optimizing land use parameters based on the V–I–S fraction is as follows:

$$vis{Q}_{xj}={H^{\prime }}_j(1-({\displaystyle \frac{vis{D}_{xj}^z}{vis{D}_{xj}^z+k^z}}))$$

(5)

where $vis{Q}_{xj}$ is the habitat quality of pixel x, ${H^{\prime }}_j$ is the habitat suitability score of type j in the V–I–S fraction rule sets, H’_j∈[0,1], z is a scaling factor to reflect spatial heterogeneity, k is a half-saturation constant, and $vis{D}_{xj}$ is the total threat level in pixel x with type j of the V–I–S fraction rule sets; the Equation is as follows:

$$vis{D}_{xj}=\sum _{r=1}^R\sum _{y=1}^{Y_r}({\displaystyle \frac{w_r}{{\sum }_{r=1}^R{w}_r}})r_y{i}_{rxy}{\beta }_x{S^{\prime }}_{jr}$$

(6)

where $R$ is the number of threat factors (r = 1, 2, 3, …, $R$), $Y_r$ is the pixel number of threat factor r (y = 1, 2, 3, …, $Y_r$), $w_r$ is the weight of threat factor r, $w_r$ ∈ [0,1], $r_y$ is the level of delay in threat factor r, ${\beta }_x$ is the level of accessibility in pixel x, ${S^{\prime }}_{jr}$ is the sensitivity of type j of the V–I–S fraction rule sets to the r type of a threat, ${S^{\prime }}_{jr}$ ∈ [0,1], and $i_{rxy}$ is the distance degradation function of the r type of a threat in pixel y on a natural habitat in pixel x.

It should be mentioned that the habitat suitability score ${H^{\prime }}_j$ and the sensitivity parameter ${S^{\prime }}_{jr}$ of class j regarding the V–I–S fraction rule set for threat factor r in the above formula need to be adjusted and corrected based on the mapping relationship of the land use/cover type and the V–I–S fraction rule set; the other parameters are the same as in the original model. The specific steps for correcting the values of ${H^{\prime }}_j$ and ${S^{\prime }}_{jr}$ are as follows: (i) list the land use/cover types corresponding to the V–I–S fraction rule set of Class j, L₁, L₂, …, L_n; (ii) the habitat suitability scores were set for the first-to-n surface types H₁, H₂, …, H_n, and the sensitivity parameters S₁, S₂, …, S_n; (iii) the mean value of H₁, H₂, …, H_n was set as the habitat suitability score ${H^{\prime }}_j$; (iv) the mean value of S₁, S₂ …, S_n sets the sensitivity parameter ${S^{\prime }}_{jr}$ of class j for the V–I–S fraction rule set for stress factor r.

4.3. Results and Analyses

4.3.1. V–I–S Fraction Results Based on the LSMA Method

To optimize the InVEST model through unmixing the V–I–S fraction, we utilized LSMA and Landsat series imagery to extract the V, I, and S fraction layers in central Guangzhou from 2000 to 2020. Additionally, we randomly selected 100 ground truth validation samples to assess the accuracy of the V–I–S fractions. This assessment involved comparing the root mean square error (RMSE) and mean absolute error (MAE) calculated from the V–I–S proportions obtained from Google Earth imagery with those derived from Landsat series imagery for the corresponding years. The accuracy assessment results indicate that the V–I–S fraction layers obtained through LSMA and Landsat series imagery exhibit a comparatively high accuracy in central Guangzhou from 2000 to 2020, as evidenced by a range of suitable indicators. Specifically, the RMSE ranges from 0 to 0.3, and the MAE ranges from 0 to 0.25.

Figure 4 displays the V, I, and S fraction layers generated through spectral unmixing with LSMA and Landsat series imagery. A value closer to 1 indicates a higher proportion of the endmember in the pixel, while a lower value indicates a lower proportion. From a spatiotemporal perspective, the distribution of the V, I, and S fractions in central Guangzhou exhibits subtle variations on an inter-annual scale. Moreover, the spatial patterns of the three endmembers are notably distinct. Pixels with a high proportion of V endmembers are predominantly concentrated in the eastern region of the study area, while the proportion of the I endmember gradually decreases from southwest to northeast. The overall proportion of the S endmember is relatively small, with pixels exhibiting a scattered distribution.

4.3.2. The Mapping Relationship of the LULC and V–I–S Fractions

To establish the mapping relationship between the V–I–S fractions and LULC types, we initially employed the ISODATA clustering algorithm to segment the V–I–S fraction rule sets based on the feature space comprising the V, I, and S fractions. The value range of each rule set was evaluated using a confidence interval of 5% to 95%, ultimately yielding the V–I–S fraction rule sets.

As depicted in

Table 4. The rule sets of the vegetation–impervious surface–soil fraction clustering.

No.	The Numerical Intervals of V–I–S Fractions Rule Sets
	2000	2005	2010	2015	2020
1	V ≤ 0.442	V ≤ 0.399	V ≤ 0.428	V ≤ 0.359	V ≤ 0.285
	I ≤ 0.561	I ≤ 0.592	I ≤ 0.545	I ≤ 0.702	I ≤ 0.729
	S ≤ 0.674	S ≤ 0.503	S ≤ 0.628	S ≤ 0.494	S ≤ 0.488
2	V ≤ 0.424	V ≤ 0.352	V ≤ 0.447	V ≤ 0.388	V ≤ 0.426
	I ≤ 0.506	I ≤ 0.530	I ≤ 0.596	I ≤ 0.667	I ≤ 0.577
	0.204 ≤ S ≤ 0.894	0.302 ≤ S ≤ 0.761	0.20 ≤ S ≤ 0.945	0.110 ≤ S ≤ 0.870	0.203 ≤ S ≤ 0.820
3	V ≤ 0.564	V ≤ 0.646	V ≤ 0.547	V ≤ 0.508	V ≤ 0.438
	0.02 ≤ I ≤ 0.612	0.02 ≤ I ≤ 0.475	0.02 ≤ I ≤ 0.655	0.180 ≤ I ≤ 0.780	0.263 ≤ I ≤ 0.816
	S ≤ 0.643	S ≤ 0.537	S ≤ 0.631	S ≤ 0.553	S ≤ 0.514
4	V ≤ 0.584	V ≤ 0.588	V ≤ 0.580	V ≤ 0.612	V ≤ 0.591
	0.294 ≤ I ≤ 0.773	0.161 ≤ I ≤ 0.745	0.278 ≤ I ≤ 0.784	0.286 ≤ I ≤ 0.820	0.329 ≤ I ≤ 0.831
	S ≤ 0.282	S ≤ 0.349	S ≤ 0.278	S ≤ 0.239	S ≤ 0.189
5	V ≥ 0.404	V ≥ 0.475	V ≥ 0.400	V ≥ 0.455	V ≥ 0.447
	I ≤ 0.498	I ≤ 0.463	I ≤ 0.490	I ≤ 0.439	I ≤ 0.482
	S ≤ 0.325	S ≤ 0.241	S ≤ 0.361	S ≤ 0.158	S ≤ 0.122
6	V ≤ 0.275	V ≤ 0.253	V ≤ 0.251	V ≤ 0.176	V ≤ 0.247
	I ≥ 0.408	I ≥ 0.502	I ≥ 0.541	I ≥ 0.729	I ≥ 0.631
	S ≤ 0.305	S ≤ 0.212	S ≤ 0.345	S ≤ 0.133	S ≤ 0.177

, six distinct V–I–S fraction rule sets were identified, each exhibiting significant differences from 2000 to 2020. These rule sets effectively convey information regarding the three endmembers (V, I, and S) obtained from Landsat series images at the sub-pixel level. Notably, there is a relatively consistent range within each rule set on an inter-annual scale, with stable overall numerical distributions. Under V–I–S fraction rule set 1, the V fraction maintains an average range of 0 to 0.383, the I fraction ranges around an upper-bound average value of 0.626, and the S fraction mainly falls within the interval range of 0 to 0.557. Conversely, under V–I–S fraction rule set 2, the S fraction exhibits a mean range of 0.204 to 0.858, indicating a stricter numerical range compared to rule set 1. The V fraction maintains an average range of 0 to 0.407, while the I fraction ranges from 0 to 0.575. Under V–I–S fraction rule set 3, the mean ranges are V $\le$ 0.541, 0.101 $\le$ I $\le$ 0.716, and S $\le$ 0.576, with the average upper bounds of the V and S fractions being less than 0.6. The upper and lower bounds of the I fraction fluctuate between 0.2 and 0.6, respectively. Similarly, under V–I–S fraction rule set 4, the lower bounds of the V and S fractions are 0, with upper bounds of 0.591 and 0.2674, respectively, exhibiting a narrow numerical range, while the I fraction ranges from 0.270 to 0.791. Under V–I–S fraction rule set 5, the mean ranges are V $\ge$ 0.436, I $\le$ 0.474, and S $\le$ 0.241, with the V fraction having relatively high values, the S fraction maintaining a small interval span, and the I fraction taking low-value intervals as the limiting condition. Finally, under V–I–S fraction rule set 6, the V fraction is strictly limited to an interval range around the upper bound average value of 0.240, the S fraction ranges from 0 to 0.234, and the I fraction ranges from 0.562 to 1, indicating a notable bias towards the I endmember. These findings elucidate the nuanced relationships between V–I–S fractions and LULC types, providing valuable insights for habitat quality assessments and land management strategies.

In this study, we established the mapping relationship between the V–I–S fractions and LULC types by determining the highest frequency of land use/cover samples observed in Google Earth imagery, which covered all LULC sample types. As shown in

Table 5. The mapping relationship of V–I–S fraction rule sets and land use/cover types.

The V–I–S Fraction Rule Sets	Land Use/Cover Types
	2000	2005	2010	2015	2020
1	-	-	-	-	-
2	unused land	unused land	unused land	unused land	-
3	grassland	-	shrubland–grassland	-	unused land
4	Wetland–cropland	wetland–shrubland– grassland–cropland	wetland	cropland	grassland-shrubland
5	shrubland–forest	forest	forest–cropland	wetland–shrubland– forest–grassland	forest–wetland– cropland
6	built-up areas	built-up areas	built-up areas	built-up areas	built-up areas

, V–I–S fraction rule sets 2 and 6 exhibit relatively stable correlations with unused land and built-up areas, respectively, from 2000 to 2020, remaining consistent over time. The LULC types corresponding to rule set 3 vary on an inter-annual scale. For instance, in 2000, grassland types correspond to rule set 3, while in 2010, it corresponds to shrubland and grassland types, and in 2020, it corresponds to unused land. Rule set 4 predominantly corresponds to vegetation types, encompassing wetland and cropland types in 2000, multiple vegetation types (including wetland, shrubland, grassland, and cropland) in 2005, only wetland in 2010, only cropland in 2015, and grassland and shrubland in 2020. Similarly, rule set 5, akin to rule set 4, primarily corresponds to vegetation types. From 2000 to 2020, forest types exhibit stable correspondence with rule set 5, including shrubland and forest in 2000, only forest in 2005, forest and cropland in 2010, and four surface types in 2015 (wetland, shrubland, grassland, and forest), which further includes forest, wetland, and cropland in 2020. Any samples lacking an assigned LULC type are automatically classified as water bodies. These findings elucidate the dynamic relationships between V–I–S fractions and LULC types over time, providing valuable insights into land use and land cover dynamics for effective land management strategies.

4.3.3. The Habitat Quality Results Based on the Sub-InVEST Model and InVEST Model

In this study, the habitat quality model was utilized to evaluate the performance of the Sub-InVEST model, leveraging the mapping relationship between the V–I–S fraction and LULC types. For a comparative analysis with the original model, we selected central Guangzhou as the testing area and obtained habitat quality data spanning from 2000 to 2020 using both the Sub-InVEST model and the original InVEST model (Figure 5). The spatial distribution of habitat quality, as depicted by both models, reveals a consistently increasing trend from the southwest to the northeast. From the perspective of spatial evolution, the habitat quality results obtained from both the InVEST and Sub-InVEST models were entirely observed to have a decreasing trend, particularly in the southwestern areas of central Guangzhou. Notably, the spatial characteristics observed in the Sub-InVEST model exhibit greater detail, with smoother transitions from low to high values compared to the results obtained from the traditional InVEST model. The latter displays a blocky spatial pattern, particularly evident in the low-value range.

To gain further insight into the numerical disparities in habitat quality derived from the Sub-InVEST and InVEST models, we present the numerical distribution ranges and averages of the habitat quality values obtained from both models in Figure 5b. Although the average values of habitat quality between the models are similar, there are notable differences in the numerical distributions. Specifically, the distribution density curve of habitat quality in the Sub-InVEST model exhibits a smoother profile compared to the InVEST model. Moreover, the probability of values within each range is significantly augmented in the Sub-InVEST results, while the habitat quality outcomes from the InVEST model appear more concentrated and segmented. For instance, there is a conspicuous peak in the distribution density within the low-value range (0.2–0.4) in the results of the traditional model, suggesting that habitat quality values are predominantly concentrated within this range. In contrast, although the distribution density within the same low-value range is also high in the Sub-InVEST results, the increasing trend is more gradual. This indicates that while the habitat quality value distribution of the optimized model reflects the numerical characteristics of the traditional model, it offers a more nuanced and spread-out representation. Overall, the spatial and numerical distribution characteristics of habitat quality derived from both models suggest that the Sub-InVEST model is well-suited for application in central Guangzhou. Utilizing the V–I–S fraction to replace land use data proves to be a viable method for optimizing the InVEST model in this context.

4.3.4. A Comparison of the Sub-InVEST Model and InVEST Model

To assess the efficacy of both the Sub-InVEST and InVEST models, we generated 100 random sample points and delineated seven sample regions (3000 m × 3000 m) within the study area, which were purposed to (i) assess the optimization impact of the models employing Landsat imagery or very-high-resolution imagery from Google Earth for analysis and validation; and (ii) quantitatively evaluate the numerical variances between the models utilizing the linear regression analysis method. The spatial distribution of the sample points and regions is illustrated in Figure 6.

Figure 7 displays the linear regression outcomes concerning habitat quality as evaluated by both the Sub-InVEST and InVEST models. From 2000 to 2020, the habitat quality results obtained from both the optimized and traditional models exhibit a notably consistent upward trend. This consistency suggests that the Sub-InVEST model and its conventional counterpart yield comparable effects in assessing habitat quality, thereby affirming the efficacy of the Sub-InVEST model in this regard. The substitution of land use data with V–I–S abundance further validates the feasibility and effectiveness of this approach. The regression analyses for the years 2000, 2005, 2010, 2015, and 2020 reveal coefficient of determination (R²) values between the traditional model’s habitat quality outcomes and those of the Sub-InVEST model as 0.41, 0.35, 0.42, 0.39, and 0.47, respectively. These values indicate a stronger alignment between the pre- and post-optimization models in 2000, 2010, and 2020, suggesting improved concordance with habitat quality outcomes. Additionally, the positive slopes observed in the linear regression models for each year signify a positive correlation between the results of the traditional model and those of the optimized model. In summary, the linear regression analyses spanning from 2000 to 2020 underscore a robust correlation and consistent linear trends between the pre- and post-optimization model outcomes, affirming their mutual explanatory power.

Figure 8 illustrates seven sample regions utilized for analyzing and validating the optimization effects of the Sub-InVEST and InVEST models utilizing Landsat imagery. These sample regions encompass diverse surface composition scenarios; for instance, sample region 1 represents a peri-urban transition zone, sample region 4 denotes a heavily vegetated region, and sample region 6 signifies a densely populated urban area. The habitat quality assessments vary significantly according to these distinct surface composition scenarios. In this study, we present the habitat quality results for the years 2000 and 2020 at the sample region scale. The Sub-InVEST model yields highly consistent results with Landsat imagery, providing appropriate detail. Conversely, while the traditional model’s results generally align with the overall spatial patterns observed in Landsat imagery, some land surface features are poorly represented or indiscernible. This suggests that Sub-InVEST offers a superior and more refined portrayal of details, which aligns more closely with remote sensing imagery.

For 2000, the Sub-InVEST habitat quality results for sample regions 1–7 reveal land surface scenarios that the traditional model failed to identify and evaluate. For instance, sample region 6 exhibits a patchy distribution of habitat quality values in the traditional model’s results, whereas the optimized model discerns various habitat quality values for complex land use types such as vegetation, wetlands, and bare areas within dense urban environments. Similarly, for 2020, the Sub-InVEST habitat quality results provide a more accurate assessment. Notably, the low habitat-quality values identified by the traditional model in sample region 5 appear fragmented and interconnected, with crude handling of composition within these zones. In contrast, the low-value zones identified by the optimized habitat quality results align more closely with the feature composition observed in Landsat imagery, offering enhanced accuracy.

5. Discussion

5.1. The Importance of Optimizing the Parameters of the InVEST Model

Ecosystem service models serve as the primary approach for quantitatively evaluating ESs, with notable examples including the InVEST model, the ARIES model, and the SolVES model. Among these, the InVEST model is widely favored in ESs assessment due to its concise formulation, ease of access to input parameters, and versatile simulation capabilities [18].

The accuracy of input parameters in the InVEST model significantly impacts the validity and precision of output results. Therefore, when applied in diverse regions, it is crucial to set model parameters locally. Of these parameters, the land use parameter holds particular importance, typically derived from land use data. Adjustments to certain model parameters are made based on these input land use/cover data. There were several studies that used remotely sensed land use classification data to assess ES in various spatial resolutions, such as 30 m [38,39], 300 m [40], 1 km [41,42], etc.

However, the lengthy update intervals of land use data and the coarse spatial resolution often limited the timeliness and accuracy of ecosystem service assessments using the InVEST model. In this study, we proposed replacing the land use data of the InVEST model with V–I–S fraction data at the sub-pixel scale to enhance the model’s timeliness and accuracy. By comparing the spatial and numerical distributions of the two models before and after optimization and employing visual analysis and numerical indicators to evaluate the optimization effect, we conducted a comprehensive optimization process. This process encompasses improvements to land use parameters, optimization testing, and the evaluation of the resultant effects.

The land use parameter optimization tests of the InVEST model involve extracting V–I–S fraction data and selecting surface type samples based on Landsat series remote sensing imagery and Google Earth imagery. A mapping relationship is established between land surface types and V–I–S fraction data, which helps to evaluate the feasibility of the Sub-InVEST model by substituting V–I–S fraction data and LULC data on a regional scale. It is applicable to applications entailing long-term time series and the high-accuracy assessment of ES, which can facilitate the direct linkage of satellite remote sensing datasets and enable the convenient selection of multi-source remote sensing data for specific scales and time intervals. For instance, high- and moderate-resolution remote sensing data, such as those from Landsat, Sentinel-1, and Sentinel-2, can be utilized at the city or county scale, while low-resolution remote sensing data can be employed on a large scale, such as urban agglomerations.

5.2. Mapping Relationship Between LULC and the V–I–S Fraction

The V–I–S model, a conceptual framework describing land surface components in remote sensing, has evolved with advancements in remote sensing technology, particularly in its integration with LSMA to unmix remote sensing images from various sources, yielding V–I–S component information at the sub-pixel scale [43,44]. Currently, the Landsat series images serve as the primary data source for extracting V–I–S fractions [45,46]. In this study, the V fraction, I fraction, and S fraction were extracted based on Landsat series images and the V-H-L-S model from 2000 to 2020, with evaluation metrics such as RMSE and MAE falling within reasonable ranges. These results indicate that V–I–S fraction data possess adequate accuracy for optimizing the InVEST model based on the mapping relationship between V–I–S fractions and LULC types. The spatial patterns of the V, I, and S fractions exhibited significant spatial differentiation. The areas with high V fractions were mainly concentrated in the eastern region of Guangzhou. The areas with high I fraction showed a decreasing trend from southwest to northeast in Guangzhou, and the areas of the S fraction were relatively small. There was a high-consistency pattern with other studies on extracting impervious surface areas in Guangzhou [47], the Guangdong-Hong Kong-Macao Greater Bay Area [48], and the Pearl River Delta [49].

The ISODATA clustering algorithm is employed to conduct clustering analysis in the feature space comprising three surface components (V, I, and S). Through unsupervised learning, different combinations of V–I–S fractions are classified, searching for rule sets with similar features in the V–I–S fraction feature space. Focusing on central Guangzhou as the test area, this study identifies six distinct V–I–S fraction rule sets, defining specific ranges for each rule set based on V–I–S fraction values within the 5–95% confidence interval. By combining V–I–S fractions from various selected land surface type samples, the V–I–S fraction rule set is determined to establish the mapping relationship between V–I–S fractions and LULC types. The results of ISODATA clustering also demonstrate the feasibility of dividing and defining V–I–S fractions in the V–I–S feature space.

Among these rule sets, rule set 1 is dominated by low-value intersections, automatically categorized as water bodies due to the absence of corresponding land surface types. Rule set 2 represents mixed pixels primarily characterized by a combination of S proportion averaging above 0.2 and smaller proportions of V and I, corresponding to unused land. Rule set 3 narrows the range restriction of the I fraction while expanding the V and S fractions, corresponding to land use/cover types of grassland, shrubland, and unused land. Rule set 4 exhibits a V fraction similar to rule set 3 but with a less restricted S fraction, corresponding to land surface types, including cropland, grassland, and wetland. Rule set 5 is dominated by the V fraction, with high proportions of V and I and restricted S proportions, corresponding to land use/cover types, including forest, wetland, and shrub. Rule set 6 provides a clearer definition, with the I proportion restricted to high-value areas and minimal proportions of S and V. These findings suggest that combinations of grassland, shrub, and cropland fractions span a wide range, primarily concentrated in the V–I–S fraction rule sets 3 and 4 with some mixed states.

5.3. Benefits of Optimized Land Use Parameters for the InVEST Model Based on V–I–S Fractions

The utilization of the V–I–S fraction offers a promising avenue for optimizing the InVEST model by replacing conventional land use data. This study leverages the habitat quality model of InVEST to assess its performance in central Guangzhou both before and after optimization. The results demonstrate that the Sub-InVEST model yields habitat quality assessments with enhanced spatial detail and smoother numerical distributions. Furthermore, the comparison involves 100 random samples and seven random sample regions, evaluating numerical accuracy and visual effects to assess parameter optimization effectiveness.

Numerical variations in habitat quality before and after optimization were evaluated using linear regression models and fitting parameters for 100 random points. The results indicate a significant positive correlation with a reasonable linear trend between the optimized and traditional models from 2000 to 2020, suggesting comparable performance in assessing habitat quality. Additionally, seven random sample regions, representing various urban and vegetated landscapes, were selected for further comparison. The optimized Sub-InVEST model exhibited notably finer detail portrayal processing compared to the InVEST model.

Detailed descriptions and comparisons of habitat quality results for sample regions 1, 3, and 4 in 2020, along with Landsat and Google Earth images, are provided in Figure 9. The optimized Sub-InVEST model demonstrated better alignment with Landsat and Google Earth imagery, particularly in depicting land use/cover features, compared to the InVEST model. Notably, sample regions 1, 3, and 4 represent suburban transition areas, dense urban areas, and lush vegetated regions, respectively. The Sub-InVEST model effectively evaluated these diverse scenarios, accurately capturing surface textures and scenes depicted in Landsat and Google Earth images.

In summary, the approach of replacing land use data with V–I–S fraction data for optimizing land use parameters in the InVEST model proves to be both feasible and effective. The Sub-InVEST model enhances assessment precision and demonstrates strong applicability across various complex surface composition scenarios, thus providing more accurate ecosystem service assessment results.

5.4. Limitations and Future Outlook

This study attempts to address the issue of the InVEST model’s land use parameters being directly influenced by land use data, thereby impacting the accuracy of ecosystem service assessments. By integrating V–I–S fraction data into the land use parameters and substituting the requisite land use data with V–I–S fraction data, this study endeavors to augment the timeliness and precision of the InVEST model. Among them, the Landsat series images were utilized in this study to extract the V–I–S fraction from 2000 to 2020 with a 30 m spatial resolution since the Landsat images are the most frequently used multispectral remote sensing data in V–I–S fraction extraction with a low cost [50,51]. Although it can provide long-term series and appropriate spatial resolution data sources for testing the Sub-InVEST model for ES mapping in this study, more extensive datasets are still needed to evaluate the performance. In the next study, we will attempt to test different performances of the Sub-InVEST model by applying a series of remote sensing datasets with higher or lower spatial resolutions, such as Sentinel-1/2 data, GaoFen satellite data, and MODIS data. In addition, it is a valuable endeavor to extract the V–I–S fraction from the fusion of optical and synthetic aperture radar to improve the spatiotemporal resolution.

Central Guangzhou serves as the test area for conducting optimization experiments with the habitat quality model. The results demonstrate the feasibility and efficacy of land use parameter optimization within the InVEST model. Moreover, the study area, characterized by typical urban ecosystem types with complex surface compositions and ecological processes, provides a robust testbed for validating the feasibility and effectiveness of the model optimization method. The optimized model exhibits finer assessment accuracy, particularly in dense urban areas with extensive construction land. In future research, it would be valuable to expand the application of the Sub-InVEST model to diverse geographic units, various scales, and multiple ecosystem types, thereby further strengthening its applicability and generalizability.

6. Conclusions

In this study, we devised an optimization method for the land use parameter of the InVEST model at the sub-pixel scale, leveraging the V–I–S model and a ML algorithm. The resulting optimized model, Sub-InVEST, was derived by establishing a mapping relationship between LULC types and V–I–S fraction data, thus replacing conventional land use data with V–I–S fraction data. To achieve this, we utilized the LSMA method and ISODATA to extract V–I–S fraction data and delineate the mapping relationship between V–I–S fractions and LULC types, drawing from Landsat series imagery and Google Earth imagery.

We conducted a comprehensive evaluation of the optimization effect of the model, both visually and numerically, employing the habitat quality model as the optimization testbed. The evaluation encompassed a comparison of model performance before and after optimization in central Guangzhou as the experimental area. Our findings reveal a high degree of spatial consistency between the optimized model and the traditional model, characterized by a prevailing trend from southwest to northeast with minimal inter-annual variation. Notably, the optimized model effectively delineates land use/cover scenes that elude identification and evaluation using the traditional model in each sample region. Moreover, the spatial descriptions of habitat quality produced by the optimized model exhibit greater detail and alignment with remote sensing imagery. In terms of numerical distribution, the optimized model’s results demonstrate a smoother density curve and significantly improved probability distribution across interval ranges compared to the traditional model. These results underscore the feasibility and efficacy of replacing land use data with V–I–S fraction data for land use parameter optimization in the InVEST model. Furthermore, the optimized Sub-InVEST model offers enhanced accuracy and versatility in ecosystem service assessment, applicable across a broad spectrum of complex surface composition scenarios.