Assessment of the urban habitat quality service functions and their drivers based on the fusion module of graph attention network and residual network

ABSTRACT

Land use/cover change is a major cause of ecological degradation. Reliable LUCC data are essential for evaluating habitat quality. The current method of surface cover classification based on the convolutional neural networks (CNNs) is usually a local spatial operation using a regular convolutional kernel, which ignores the correlation between adjacent image elements. This paper proposes a combination network with two branches, branch 1 uses the K-nearest neighbor clustering algorithm to construct superpixels and then uses the data transformation module to construct a graph attention network (GAT); branch 2 constructs the CNN using attention and residual modules to obtain the spatial and higher-order semantic information of the images. Finally, the features are fused using weighted fusion, and a classification map with less point noise and greater consistency with the real surface coverage is obtained. The classification results of this network are better than those of the other competitive methods. In addition, the urbanization of Sanya has resulted in significant habitat degradation. A good fit (${\mathit{R}}^2$ in 2020 = 0.639) between habitat quality (HQ) and natural and socioeconomic factors was observed in Sanya. Natural factors are more relevant to HQ than socioeconomic factors and vary spatially.

KEYWORDS:

• Residual network
• graph attention network
• super-pixel segmentation
• habitat quality
• driving force analysis

1. Introduction

Habitat quality (HQ) reflects the strengths and weaknesses of natural resources and environments on which humans depend. It is a notable symbol of ecosystem health and indicates the strength of ecosystem services and the ability of a species to provide a living environment (Fellman et al. 2015; Hillard, Nielsen, and Groninger 2017; Newbold et al. 2015). As an important ecosystem service, it has a very important position, and knowledge of the spatial and temporal mechanisms of HQ can help regulate human behavior in the use of modified resources and the environment and achieve regional sustainable development (Yin et al. 2020). LUCC is the most direct signal for characterizing the impact of human activities on natural ecosystems on the Earth’s land surface and is the link between human socio-economic activities and natural ecological processes (Dai et al. 2019; Terrado et al. 2016). Land-use change can directly alter the distribution of biological habitats and resources, which in turn affect ecosystem processes. Therefore, studying the effects of land-use change on HQ can provide a theoretical basis for the construction of regional ecological security patterns and sustainable use of land resources (Otto et al. 2016; Sala et al. 2000; Song et al. 2020). In recent years, China has experienced rapid urbanization, with the urbanization rate of its resident population increasing from 36.22% in 2000 to 63.89% in 2020. Rapid urbanization has considerably changed the integrity, distribution pattern, and function of the original habitats in the region, and the service capacity and health status of regional ecosystems have been greatly impacted, thereby affecting regional ecological security (Huang et al. 2018; Ke et al. 2018). Exploring regional HQ assessments and analyzing their intrinsic change patterns are important for urban ecological security. Scholars have explored the spatial and temporal distribution of HQ, driving factors, and prediction models from different perspectives and scales, combined them with concepts such as landscape patterns, land-use changes, urban expansion, and human activity intensity, and have achieved multi-scale and multi-level research results (Chu, Lu, and Sun 2022; Song et al. 2020). Satir et al. studied the effects of rapid urbanization and rising transport infrastructure on ecological service functions, showing that ecological service functions near transport facilities were substantially reduced, which was related to the intensity of human activities (Satir et al. 2023). In most studies, it declined, whereas the higher intensity of human activities significantly changed the regional landscape pattern and affected the ecological processes of biodiversity. Therefore, it is crucial to investigate the effects of urbanization on it (Wei et al. 2022; Wei et al. 2023).

In addition, natural factors (altitude, NDVI, and climate) strongly influence HQ variability (Zhang et al. 2022; Zheng and Li 2022). Therefore, the identification of factors affecting HQ can provide essential information for maintaining biodiversity (Ahrends et al. 2015). Although previous studies have focused on spatiotemporal evolutionary features and the factors affecting them, further research on this topic is required. Several methods have been applied to explore the factors affecting the HQ, such as correlation analysis (Yohannes et al. 2021), ordinary least squares (OLS) models, geographically weighted regression (GWR) models (Zhu et al. 2020), and GeoDetector (Xie and Zhang 2023). GeoDetector can only explore the interpretive force of single- or multi-factor interactions of each factor, whereas GWR models can explore the correlation of impact factors on their spatial location; however, there are still fewer experiments using GWR models to explore impact factors and apply them to HQ. The exploration of the impact factors is crucial for sustainability.

Several scholars have used InVEST, SolVES, landscape ecological risk, and HIS models to assess habitat quality (Qing et al. 2021; Seoane et al. 2006; Zhu et al. 2020). The InVEST model stands out among these models and has the advantages of convenient access to data, few required parameters, strong spatial representation and visualization ability, accurate analysis, and simple data processing (Berta Aneseyee et al. 2020; Sun et al. 2021; Wang and Cheng 2022; Wu, Sun, and Fan 2021). Studying the influence of land use change on ecosystems is important for promoting ecology, improving ecosystem management, guiding urban planning and construction, and raising citizen awareness of ecological conservation (Fu et al. 2022; Li et al. 2021). Considering the input data of InVEST, all data are studied based on land use data, and it is important to obtain high-accuracy land use data; however, when scholars used InVEST to study regional habitat quality, most of them were based on finished data, which are generally large-scale, and the accuracy is still unsatisfactory. Yu et al. showed that previous research has largely underestimated the impact of LUCC in ecosystem assessment and also emphasized the importance of reliable LUCC data in ecosystem service function assessment (Yu et al. 2022). Ersoy Mirici, Satir, and Berberoglu (2020) used a hybrid classification approach for land use cover classification that takes full advantage of almost every machine learning classification, using landscape metrics to determine the habitat quality of the area (Syrbe and Walz 2012). The advantage of using this method is that the habitat quality of each land class can be rapidly counted, whereas using the InVEST model makes it easier to visualize the habitat quality of the whole region; however, the disadvantage of the landscape metrics method is that more metrics have to be analyzed. Wang et al. used the method of deep learning to produce LUCC data and used it as input for the InVEST model (Wang et al. 2023a). Based on this idea, we first used deep learning to produce LUCC data which were subsequently used to assess HQ.

Remote sensing images have been widely used in environmental assessment, crop extraction, disaster assessment, and prediction because of their richness in information and have received wide attention from scholars because of their classification with important applications (Gyaneshwar and Nidamanuri 2021; Wang et al. 2022a; Xie et al. 2022; Zhang et al. 2020). Initially, some simple linear classifiers in machine learning were used for classification, such as random forests and K-nearest neighbor algorithms, which were widely used because of their simplicity and the good classification accuracy they could achieve in general (Belgiu and Drăguţ 2016; Xue et al. 2014; Xue, Du, and Su 2014; Zhang and Zhou 2007). However, the above methods are based on handcrafted spectral space features, which rely heavily on rich expertise, and the features are generally shallow, severely limiting classification accuracy. To address these shortcomings, deep-learning technology has been extensively utilized for feature classification tasks by relying on powerful self-learning capabilities and feature representations (Sapkota, Sharma, and Mann 2022; Tulapurkar, Banerjee, and Buddhiraju 2023).

Deep learning techniques can automatically extract features, effectively avoiding classification errors caused by manual involvement in the feature design and threshold setting. With the development of classification networks, graph neural networks have received increasing attention owing to their ability to handle arbitrarily structured data (Kipf and Welling 2016) because CNN can utilize the spectral and local spatial features of image pixels but ignore the global contextual information and the interaction between spatiotemporal information. In addition, there is a large variability between pixels of the same or different categories in the region, which can lead to misclassification of boundary pixels in the classification process (Guo et al. 2020). CNN is designed for Euclidean data to deal with regular spatial structures, and thus cannot easily capture the internal connections of different land covers in images (Liu et al. 2022). To solve these problems, classification networks of GNNs have been designed. GNNs can learn adaptive kernel parameters based on the specific distribution of land cover and perform flexible convolution in arbitrary irregular land cover regions (Xu et al. 2022). Thus, the GNN was used to learn the correlations between different land cover types and model their spatial topology on graphs (Wan et al. 2020). Dong et al. (2022) designed a fusion network combining a GAT network and a depth-separable convolutional network, which was able to combine and fuse the features extracted from two different convolutions and obtained better classification results; Liu et al. (2023) designed a self-supervised classifier using graph convolutional networks, which can produce accurate class images despite the small sample size; with the development of CNN, several excellent models have emerged, such as DFInet, which addresses the challenges present in the wetland classification task. Its module enhances the utilization of complementary information from multispectral and hyperspectral data and preserves the differential features of the target (Gao et al. 2022). Good quality raw data can improve the accuracy of the classification task, and CSMFormer solves both problems by allowing not only data fusion, but also feature classification (Gao et al. 2023a). ACL-CNN is mainly used to discover the complementary advantages of multimodal data and suppress noise through PSM before performing feature classification, outperforming DFInet on specific datasets (Gao et al. 2023b). RSRNet proposes modal enhancement and semantic enhancement to address the accumulation of errors in the model after multiple rounds of learning. A status replay strategy used on the classifier can make the decision boundary more stable, in addition, a module combining multiple sources of information to jointly optimize the parameters is proposed, and RSRNet is better than the above competitive models in some cases (Wang et al. 2023b). The innovation of MGSNet is the use of target background information to aid decision-making, and its well-designed network, which not only improves the distinguishability between target-backgrounds, but also improves the expression of similarity between targets, which makes full use of the texture features of the image (Wang et al. 2023c). The structural information of multi-source data is not well handled in feature extraction, which leads to underutilization of the complementary information of multi-source data, weakening the classification accuracy. SOT-NET is used to increase the correlation between different types of data, which consist of three branches to optimize the information of multi-source data in the transmission. Experiments have shown that it still exhibits a good semantic classification ability with a small number of samples (Zhang et al. 2023a). ADGAN uses a single classifier, which solves the problem of GAN classification being self-contradictory due to unbalanced data samples. And an AdapDrop module is carefully designed, which can generate adaptive shapes and performs well in dealing with objects of different shapes and sizes (Wang et al. 2021). MDA-NET solves the problem of collaborative integration of multi-source data and cross-scene classification in domain adaptive tasks with its established adapters for different modalities, and the information extracted from different modalities is aggregated by mutual aid classifiers, which provide better results than other competing DA methods (Zhang et al. 2023b). These advanced networks have more or less residual connections; residual connections have gained increasing attention because they can fuse low-level features with high-level features. Kong et al. (2021) used residual networks for scene classification because the network learns more features, thus reducing misclassification. The main contributions of this study are as follows:

1. To fully fuse the spectral and spatial features of images, a two-branch graph convolutional fusion network is proposed, where the GNN branch can fuse the node information using the topology of the graph. The graph nodes are the superpixel segmentation units in the image, and the edges are the relationships between neighboring nodes.

2. The CNN branch passes the processed pixel-level spectral features to the network, in combination with an attention module, to capture spatial and higher-order semantic information. The CNN branch is constructed using a residual network to combine the features initially extracted by the network with high-level features.

3. The algorithm combines the spectral differences between features and the spatial relationships between image elements to significantly reduce misclassification of features and point noise in the classification map. The method proposed in this study can obtain more discriminative feature representations and improve the accuracy of ground cover classification.

4. The InVEST model simulated HQ in Sanya City, and the MGWR model investigated the factors influencing HQ.

2. Materials and methods

2.1. Research area

Sanya City is located at 18°09'34″∼18°37'27″$\mathrm{N}$, 108°56'30″∼109° 48'28″$\mathrm{E}$ (Figure 1), at the southernmost tip of Hainan Island, which is a tropical monsoon maritime climate zone. The topography is high in the north and low in the south. It is north of the mountains and south of the sea and city forms a ‘mountain-river-sea’ pattern. The east–west length is 91.6 $\mathrm{k}\mathrm{m}$, and the north–south length is 51.75 $\mathrm{k}\mathrm{m}$. The total land area is 1921.51 $\mathrm{k}{\mathrm{m}}^2$, and the sea area is 3226 $\mathrm{k}{\mathrm{m}}^2$. It comprises of four districts under the jurisdictions of Yazhou, Tianya, Jiyang, and Haitang. The coastline starts from Fuwan in Fujikiao in the east and ends at Jiao Tou Bay in Meishan in the west, with a length of 209.1 $\mathrm{k}\mathrm{m}$. There are 19 harbors and bays, 40 islands of various sizes, rich wildlife resources of more than 300 species, and various types of vegetation (Lei et al. 2022). As a transport hub in southern Hainan and an important port along the southeastern coast under China’s opening-up policy, Hainan Island is an important ecological functional area in China and a global biodiversity hotspot. The urbanization level of Sanya City on Hainan Island has been increasing, and plenty of arable and unused land has been occupied by construction. Currently, the tourism industry is well developed in Sanya Bay, but the erosion of the coastline is becoming increasingly serious due to the lack of management. Land use is mainly tourism development, which has had an immeasurable impact on the land and ecological environment. Long-term time series of the temporal–spatial evolution patterns of HQ need to be urgently investigated.

Figure 1. Map of Sanya administrative region.

2.2. Data sources and processing

The remote sensing images from 2000 to 2020 are Landsat series data, which should be selected for the same period of the year as much as possible. However, owing to the constraints of both data quality and quantity, we selected data with similar time and high image quality. The images were subsequently radiometrically calibrated, atmospherically corrected, and cropped using the administrative vector map of Sanya City, and to quantitatively analyze the habitat quality, we uniformly used the administrative boundary of Sanya City in 2020. The land types were divided into six categories: arable land, forest land, water bodies, grassland, construction land, and unused land. The interpretation of remote sensing images is shown in Table 1 (Wang et al. 2022b). Night-time lighting data for 2000 were obtained from Chen et al. (2023), and all other data in this experiment were obtained from the Resource and Environmental Science and Data Centre of the Chinese Academy of Sciences (https://www.resdc.cn/).

Table 1. Land use classification and remote sensing interpretation signs in Sanya City.

Land Use Type	Interpretive Marker
Cropland	Mainly green and dark green, with regular shapes and concentrated distribution, and generally distributed around buildings
Forest	Distributed near relatively undulating mountains with rough shapes and blurred boundaries, it is widely distributed in the north and southeast of the study area, making it easy to distinguish
Grassland	Mainly green and dark green, with Sporadic distribution
Water Land	The geometric features are obvious, the shape is partially flat and straight, showing a linear distribution, and it becomes blue or black under certain
Built land	In true color images, it is mainly bright white, with a relatively scattered distribution, and presents black purple under certain band combinations.
Unused land	Concentrated in mountainous and overgrown areas, sporadic distribution, generally showing brown color on true color images

2.3. Methods

The proposed remote sensing image classification network consists of two branches: branch 1 is composed of a GAT network to extract superpixel-level features of remote sensing images; branch 2 is composed of the residual network and attention module to extract pixel-level features of images. Remote sensing images are preprocessed using a convolutional network before entering the two branch networks to remove redundant information (Figure 2).

Figure 2. The framework of the FGR model.

2.3.1. Data cleaning module

The data cleaning module consists of two 1 × 1 convolutional layers, which coarsely process the remote sensing data, remove some noise and redundant information, and improve the subsequent network’s ability to recognize the features, while compressing the data and initially extracting some features, which can reduce the computational cost of the whole model. For Landsat data with a small number of spectra, the module can appropriately increase the data dimensions to utilize the original band information fully. The specific implementation method is expressed by the following equation: (1) ${\mathit{X}}_{h,w,n}^{\ast }=\sigma \left({\sum }_{i,j,b}{\mathit{K}}_{i,j,b,n}{\mathit{X}}_{h,w,b}+{\mathit{B}}_{i,j,b,n}\right)$(1) where ${\mathit{X}}_{h,w,n}^{\ast }$ is the output feature mapping, and the input Landsat image features are ${\mathit{X}}_{h,w,b}$ h, w, and b denote the length, width, and number of bands, respectively. ${\mathit{K}}_{i,j,b,n}$ is the convolution kernel of the input feature map in row i, column j, and channel b, and n denotes the number of kernels. ${\mathit{B}}_{i,j,b,n}$ is the bias, $\sigma (\cdot )$ denotes the leaky ReLU layer. The input features pass through a 1 × 1 convolutional layer, and the number of channels changes from b to n. This achieves either a dimensionality reduction or increase in the data. In this experiment, b = 7 and n = 128.

2.3.2. Pre-activated residuals module

With the development of CNN, residual connections can directly propagate the results of shallow stacked layers into deep stacked layers. A series of residual networks emerged based on residual connections. Figure 3 indicates that the activation function (ReLU and BN) is placed before the convolution layer. This concept creates a new residual unit compared to the traditional concept of post-activation, which significantly improves the performance of ResNet. The complexity of the model is reduced by the pre-activation residual unit module, which allows the model to converge faster. It is a nonlinear activation function whose expression and derivative from (ReLU’) are shown in Eqs. (2) and (3). (2) $ReLU\left(x\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left(x,0\right)$(2) (3) $ReL{U}^{\mathrm{\prime }}\left(x\right)=\{\begin{array}{c}1,x⩾0\\ 0,x<0\end{array}$(3) where x denotes the eigenvalue of the input and $\mathrm{m}\mathrm{a}\mathrm{x}\left(x,0\right)$ denotes the maximum of the two. As shown in Eq. (3), when the learning rate of the network is large, the back-propagation process generates a large gradient in the ReLU layer, which makes the bias and weights updated in this layer too small, resulting in the output of this layer being zero in the next round of forward propagation; the relevant parameters are no longer updated, eventually leading to the overall data distribution leaning toward the upper and lower limits of the ReLU function and the deactivation of neurons. Thus, the BN layer is used to forcibly convert the data before the input ReLU layer into a standard normal distribution, as shown in Eq. (4): (4) ${\mathit{x}}_{ϵ}={\displaystyle \frac{\mathit{x}-E\left[\mathit{x}\right]}{\sqrt{V_{ar}\left[\mathit{x}\right]}}}$(4) (5) $y=\gamma {\mathit{x}}_{ϵ}+\text{β}$(5)

Figure 3. Residual structure.

where$x$ is the matrix, $\mathrm{E}\left[\mathit{x}\right]$ is the mean of the smallest batch, and ${\mathrm{V}}_{\mathrm{a}\mathrm{r}}\left[\mathit{x}\right]$ is the variance of the smallest batch. However, forcing the data distribution in this manner will result in some features not being learned. Therefore, we introduce $\gamma$ and $\text{β}$, while making the model learn these two parameters adaptively by back-propagation to obtain the output$y$ of the BN layer.

2.3.3. Data transformation and construction of superpixel-based graph attention networks

The GNN can process non-Euclidean data and learn both node feature and structural information. It enhances the interaction between image elements in the neighborhood, removes irrelevant and redundant information, and improves the characterization ability of the feature information. In addition, neighborhood information helps obtain global information when it is passed between layers. The GAT is an attention mechanism in graph convolution that assigns different weights to nodes, causing the network to focus more on important nodes.

Although we can use image pixels as nodes in the graph, the GNN requires an input adjacency matrix, which results in a high computational cost. To reduce the computational cost and complexity of the algorithm, we performed superpixel segmentation of the raw multispectral data using simple linear iterative clustering (SLIC), which uses the K-means algorithm to divide the raw remote sensing image into multiple superpixels that are spatially connected and have similar internal spectral characteristics, further reflecting the spatial structure of the image (Achanta et al. 2012). Through this clustering segmentation operation, the image is transformed into topologically structured graph data ($\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{h}=\left\{\mathrm{V},\mathrm{E}\right\}$), where V and E represent nodes and edges, respectively, and graph attention allows the neural network to learn useful graph representations by selectively focusing on nodes.

For example, we partition the remote sensing data into Z regions (Z = 7 in Figure 4), which can be expressed as $S=\{S_i{\}}_{\dot{i}=1}^Z$, where S represents the set of superpixels, ${\mathcal{S}}_i=\{x_j^i{\}}_{j=1}^{N_i}$ represents the i-th superpixel, $x_j^i$ represents the j-th pixel in ${\mathcal{S}}_i$, and $N_i$ represents the total number of pixels in ${\mathcal{S}}_i$.

Figure 4. Generation of an association matrix and graph structure, where $x_i$ is the i-th pixel in the flattened input image, and $v_j$ is the centroid (mean) of the j-th superpixel $S_j$.

Because the size of the partitioned regions is not the same after superpixel segmentation, resulting in different numbers of pixels in different regions, it is impossible to apply the segmented data directly to the branch of the GNN network. We designed a data transformation module to allow the features to be propagated at the pixel and superpixel levels by flattening the original features after segmentation and by directly flattening. From this, we set $O\in R^{HW\ast Z}$ $\left(HW=\sum _{i=1}^Z{N}_i=n\right)$ to be the correlation matrix between the two, where Z is the number of superpixels, and the matrix $O_{i,j}$ can be obtained according to the following equation: (6) ${\mathit{O}}_{i,j}=\{\begin{array}{cc}1,& \mathrm{i}\mathrm{f}{\overline{\mathit{X}}}_i\in S_j\\ 0,& otherwise\end{array}\overline{\mathit{X}}=Flatten\left(\mathit{X}\right)$(6) where Flatten indicates that the data features are spread into a one-dimensional vector. $O_{i,j}$ denotes the value of O at position (i, j). ${\overline{X}}_i$ denotes the ith pixel in $\overline{X}$. After obtaining the correlation matrix, the transformation is performed using the following equation: (7) $\mathrm{V}=\mathrm{E}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{r}\left(\mathit{X};\mathit{O}\right)={\overline{\mathit{O}}}^T\mathrm{F}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\left(\mathit{X}\right)$(7) (8) $\tilde{\mathit{X}}=Decoder\left(V;\mathit{O}\right)=Reshape\left(\mathit{O}V\right)$(8) where $\overline{O}$ denotes O normalized by the column. $V$ denotes the node consisting of superpixels, and $Reshape(\cdot )$ means that the spatial dimension is being reassembled. $\tilde{X}$ denotes the features transformed back into a grid.

After converting the pixel-level features into super-pixel-level features, the graph attention module is applied to process them. The structure of the GAT is shown in Figure 5. The central node has four neighboring nodes, in this paper, only first-order neighboring nodes are considered, and the current node is represented by aggregating neighbor information. Specifically, to obtain adequate network representation, we first perform a linear transformation of the input node features $\left(V=\left[v_1,v_2,\cdots ,v_Z\right]\right)$ using a weight matrix $W\in R^{B\times C}$, where C is the number of channels of v, and then compute the attention weight coefficients of the input vectors using a self-attentive mechanism a: (9) $e_{ij}=a\left(\mathit{W}{v}_i,\mathit{W}{v}_j\right)$(9) where $W$ denotes the weight parameter matrix of the nodes transformed from the input to the output feature dimension. $e_{ij}$ represents the attention weight of node j with respect to node i. Next, Softmax normalization is performed to obtain the final attention weight coefficient ${\alpha }_{ij}$ which is shown below: (10) ${\alpha }_{ij}=softmax\left(e_{ij}\right)={\displaystyle \frac{exp\left(e_{ij}\right)}{{\sum }_{k\in N_i}exp\left(e_{ik}\right)}}$(10)

The self-attention mechanism can be parameterized by a weight vector and nonlinearized using the activation function; thus, the expression for the attention weights can be further derived as (11) ${\alpha }_{ij}={\displaystyle \frac{exp\left(L\left(a^T\left[\mathit{W}{v}_i\parallel \mathit{W}{v}_j\right]\right)\right)}{{\sum }_{k\in N_i}exp\left(L\left(a^T\left[\mathit{W}{v}_i\parallel \mathit{W}{v}_j\right]\right)\right)}}$(11) where $a^T$ is the transpose of the attention weight vector; L denotes the LeakyReLU activation function; and $\parallel$ is the splicing operation.

Figure 5. Graph attention layer.

After obtaining the attention weight coefficients, a weighted sum is obtained by summing the central node. The output feature vector $v_i^{\mathrm{\prime }}$ of (12) $v_i^{\mathrm{\prime }}=\sigma \left({\sum }_{j\in N_i}{\alpha }_{ij}\left(\mathit{W}{v}_i\right)\right)$(12) where $v_i^{\mathrm{\prime }}$ is the new feature vector of the node. The graph attention layer uses a multiheaded self-attentive mechanism to capture different information, and the multiheaded self-attentive mechanism obtains more comprehensive information by calculating the K sets of attention independently. (13) $v_i^{\mathrm{\prime }}={\parallel }_{k=1}^K\sigma \left(\sum _{j\in N_i}{\alpha }_{ij}^k\left({\mathit{W}}^k{v}_j\right)\right)$(13) where K is the number of attention heads; ${\alpha }_{ij}^k$ and $W^k{v}_j$ are the attention coefficients and weight parameter matrices of the k-th group of self-attentive mechanisms, respectively.

2.3.4. Attention module

We were inspired to construct a dual attention mechanism consisting of spatial and channel attention mechanisms suitable for remote sensing images (Zhang et al. 2021), where the graph attention layer in branch one learns the feature relationships between superpixels, and a well-designed spatial attention module (SAM) using three parallel convolutions with a small receptive field is used to learn the features of small objects, as shown in Figure 6. The SAM for a given feature vector is $A\in {\mathbb{R}}^{C^{\ast }{H}^{\ast }W}$. The SAM generates two new feature vectors $Q\in {\mathbb{R}}^{C^{\ast }{H}^{\ast }W}$ and $K\in {\mathbb{R}}^{C^{\ast }{H}^{\ast }W}$ using two independent convolution operations, respectively. Matrix multiplication of Q and K is performed, and the normalized spatial attention weight matrix $U_{ij}\in {\mathbb{R}}^{N^{\ast }N}$ is obtained by applying the Softmax function. (14) ${\mathit{U}}_{ij}=softmax({\mathit{Q}}_i\cdot {\mathit{K}}_j)={\displaystyle \frac{exp({\mathit{Q}}_i\cdot {\mathit{K}}_j^T)}{{\sum }_{i=1}^N{\sum }_{j=1}^{N}exp({\mathit{Q}}_i\cdot {\mathit{K}}_j^T)}}$(14) where N denotes the feature map size, i.e. N = H*W; $Q_i$ and $K_j$ are the feature representations of the points i and j in Q and K, respectively; $U_{ij}$ is the similarity between points i and j, and the higher the similarity, the larger the value, which reflects the correlation between the points to some extent. The new feature vector $V\in {\mathbb{R}}^{C^{\ast }{H}^{\ast }W}$ is obtained by applying the convolution operation to the feature vector A, and multiplying it with the attention weight $U_{ij}$ and finally with the scale parameter α. It is then summed element-wise with the input feature A to obtain the enhanced output feature F. (15) $\mathit{F}=\alpha {\sum }_{j=1}^N\left({\mathit{U}}_{ij}{\mathit{V}}_j\right)+{\mathit{A}}_i$(15) where α is a learnable scale parameter with an initial value of 0.

Figure 6. The structure of the SAM module.

As shown in Figure 7, the channel attention module (CAM) first performs max pooling (MaxPool) and average pooling (AvgPool) on the input feature vector A in the spatial dimension, so that it retains only the channel dimension, and the other dimensions are 1 to obtain a one-dimensional vector. The outputs are then summed element-wise after a shared weight MLP, and a sigmoid activation function is passed to obtain the channel attention weights. Finally, the weight and input feature vector A are multiplied to obtain feature vector $F$ after the channel attention mechanism. (16) $\mathit{F}=\sigma \left(\begin{array}{c}MLP\left(MaxPool\left(\mathit{A}\right)\right)\\ +MLP\left(AvgPool\left(\mathit{A}\right)\right)\end{array}\right)\mathit{A}$(16)

Figure 7. The structure of the CAM module.

2.3.5. Estimation of habitat quality

The InVEST model has proven to be a reliable tool (Benra et al. 2021; Luan, Huang, and Zheng 2023). To run the habitat quality model, the following five datasets are generally required as input: current LUCC data, threat factor element table, threat source raster layer dataset, threat source sensitivity table, half-saturation, and K constant values. Land cover data were obtained according to the classification of the deep learning module, the threat source raster layer dataset was produced using ArcGIS software, and the threat factor element table and threat source sensitivity table were obtained from previous studies.

The habitat quality of land cover or habitat type j at raster cell x ($Q_{xj}$) is (17) $Q_{xj}=H_j\times \left(1-{\displaystyle \frac{D_{xj}^z}{D_{xj}^z+k^2}}\right)$(17) where $H_j$ is the habitat suitability of land-cover type j, $D_{xj}$ is the habitat degradation of raster cells x in land type j, k is the half-saturation constant, and Z is the default model parameter. The more sensitive the habitat type is to the threat factor, the more it is affected by the threat factor, and the more degraded the habitat. The degree of habitat degradation ($D_{xj}$) is calculated by the formula (18) $D_{xj}={\sum }_{r=1}^R{\sum }_{y=1}^{Y_r}\left({\displaystyle \frac{W_r}{{\sum }_{r=1}^R{W}_r}}\right)\times r_y\times i_{rxy}\times {\text{β}}_x{S}_{jr}$(18) where $i_{rxy}$ are 2 types: (19) $i_{rxy}=\{\begin{array}{cc}1-{\displaystyle \frac{d_{xy}}{d_{rmax}},}& \mathrm{L}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{y}\\ exp\left({\displaystyle \frac{-2.99d_{xy}}{d_{rmax}}}\right),& \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{y}\end{array}$(19) where$d_{xy}$ is the linear distance between raster cell x (habitat) and raster cell y (stressor), and $d_{rmax}$ is the maximum influence range of stressor r. The required parameters were set in conjunction with the recommended values of the InVEST model, literature (He, Huang, and Li 2017; Rahimi, Malekmohammadi, and Yavari 2020; Terrado et al. 2016), and expert recommendations.

2.3.6. Multi-scale geo-weighted regression models

In the traditional global regression model, the relationship pattern between variables is explored based on the spatial homogeneity of the observations. However, the factors influencing HQ vary with the geographical location of the observation sites, which shows non-smoothness of the spatial data. To avoid distortion of the estimation results owing to the assumption of homogeneity, the MGWR model was used to investigate the response of HQ to climate. Compared with the classical GWR model, the MGWR model adds spatially smooth variables to the GWR model. The expressions are (20) $y_i={\text{β}}_o(u_i,v_i)+{\sum }_{j=1}^k{\text{β}}_{bwj}(u_i,v_i)x_{ij}+{ϵ}_i$(20) where $\left(u_{\mathrm{i}},v_i\right)$ are the geospatial coordinates of sample point i (i = 1, 2, … , n); $y_i$ and $x_{ij}$ are the measured values of the dependent variable y and the j-th explanatory variable x at spatial position $\left(u_i,v_i\right)$, respectively; ${\text{β}}_o\left(u_i,v_i\right)$ is the location i intercept; ${\text{β}}_{bwj}\left(u_i,v_i\right)$ is the local regression coefficient of j variables at location i. ${ϵ}_i$ is the error term of the model. In this study, $y_i$ denotes habitat quality and $x_{ij}$ denotes impact factors.

3. Results

3.1. Setting of experimental parameters

Considering that the kernel size is directly related to the computational cost, we set the kernel size of the convolutional layer of the first and second residual units to (3, 3), (3, 5), (5, 3), and (5, 5), respectively. As shown in Table 2, the effect of the kernel size on the results inspired us to use (3, 3), where the three regions were arbitrarily selected areas of Sanya City. The effect of weighted fusion coefficients η of CNN and GAT results on the results is shown in Figure 8. We set η to 0.1 because all three regions have high accuracy at this point. The framework configuration of the FGR net is presented in Table 3.

Figure 8. Classification ability of FGR-net with different weighted fusion coefficients η.

Table 2. Classification performance of FGR-net with different kernel size k on different datasets. (3, 3) denotes the convolutional kernel of the convolutional layer of the first residual unit and the second residual unit is 3 × 3.

	(3,3)	(3,5)	(5,3)	(5,5)
Region I	98.61	97.94	95.58	93.43
Region II	99.71	99.43	99.39	99.58
Region III	99.58	99.06	98.54	99.09

Table 3. Layer configuration of the FGR-net. The number of hidden layer neurons and the number of GAT heads are 30 and 3, respectively; 0.1 is the weighted fusion coefficient, which represents the final result as 0.1 × CNN_result + 0.9 × GAT_result.

	K and η	Branch one	Branch two
conv	1 × 1@128	–	–
	1 × 1@128	–	–
GAT/P-ResNet(1)	-	30,3	3 × 3@128
FC/P-ResNet(2)	-	64	3 × 3@64
Weighted feature fusion	0.1	-	-

3.2. Pre-activated residual network

We designed a weak classifier (containing only one residual unit) to explore the improvement in model performance due to pre-activation. The Indian Pines dataset was used for this experimental test, and it can be observed from Figure 9 that as the epoch increases, the pre-activation-improved residual network loss decreases faster, and the OA improves faster, suggesting a higher application value owing to the faster fitting speed.

Figure 9. Performance comparison of pre-activated residual networks with conventional residual networks.

3.3. Comparison of classification performance with other competing methods

To compare the classification performance of the experimental model, we compared the model of this experiment with machine learning methods (SVM) and deep learning methods (2d-cnn, ResNet, 3d-cnn, GCN), and selected a typical area in Sanya City (Figure 10), which includes four land types: water, cultivated land, forest land, and building land The user’s accuracy of each land type, overall accuracy (OA), and kappa coefficient were compared separately.

Figure 10. Comparison of the accuracy of FGR-net with other classical models. (a) True color image. (b) Ground truth. (c) SVM. (d) 2D-CNN. (e) ResNet. (f) 3D-CNN. (g) GCN. (h) FGR-net.

From the classification results of typical features in Sanya City, we can see (Table 4) that although the samples for training in this experiment are all manually labeled, and the labeling on roads and some linear objects may be unbalanced with other samples, their classification results are still relatively satisfactory. Because of the need for follow-up studies, and selected the best result from ten experiments. we have added standard deviations to the table to demonstrate the stability of the model. The FGR-net is less stable than SVM, but has good robustness compared to other competing methods. Furthermore, SVM and 2D-CNN show the worst classification effect, in which the classification effect of the SVM algorithm is even higher than that of 2D-CNN. The overall accuracy of the SVM algorithm is also higher than that of the 3D-CNN algorithm, which shows that the SVM algorithm has its own unique superiority, but there is a large amount of point noise in its classification results. The CNN network is stronger than the SVM algorithm in suppressing noise; among them, the overall accuracy and Kappa coefficient of the GCN and ResNet algorithms are high, and their classification effects on water bodies are better than those of the SVM, 2D-CNN, and 3D-CNN. However, the accuracy of the FGR-net algorithm proposed in this paper is superior to that of the GCN, and the overall accuracy and Kappa coefficient are the highest; the network does not increase the computational cost, and its training speed and testing speed are the same as those of other networks.

Table 4. Comparison precision results obtained by various models on the dataset (Standard deviation in parentheses).

	SVM	2D-CNN	ResNet	3D-CNN	GCN	FGR-net
Water Land	63.49(0)	74.53(0.68)	92.94(0.92)	90.23(0.71)	95.31(1.07)	96.61(0.21)
Cropland	83.76(0)	80.68(2.26)	98.16(20.69)	94.99(11.8)	98.56(4.99)	98.98(2.69)
Forest	98.43(0)	96.89(0.25)	97.45(0.98)	99.34(0.23)	99.90(0.65)	100.00(0.29)
Built Land	95.83(0)	82.79(0.83)	98.46(4.73)	97.49(1.89)	99.47(2.94)	99.64(0.92)
OA (%)	94.70(0)	89.16(0.42)	98.09(1.10)	91.47(0.52)	98.31(0.59)	98.81(0.13)
Kappa × 100	91.43(0)	90.34(0.67)	96.58(1.74)	94.74(0.85)	98.32(0.96)	99.43(0.20)

The classification results are shown in Figures 10 and 11. The 2D-CNN shows a large number of noise points, the ResNet algorithm reduces the noise, the GCN network obtains smoother results, and the ensuing problem is the misclassification of some small features. The FGR model combines the advantages of the two networks, and by assigning different weights to the CNN and GNN branches, the classification results are guaranteed in a smoother image. The classification result is more consistent with the real features, and there is almost no point noise in the result, indicating that is is an excellent classifier.

Figure 11. Ablation visualization results of the FGR model. (a) True color remote sensing image; (b) Visualization of the selected typical areas by branch one; (c) Visualization of the selected typical areas by branch two; (d) Visualization of the selected typical areas by FGR-net (η = 0.1).

To compare the fusion effect of the network in this study on the branch 1 and branch 2 more intuitively, we compared the results obtained from each process of the network in Figure 11. After the weighted fusion of the two network branches, the classification effect depends on the branch one network for large objects, which can maintain a relatively smooth effect. The complexity of branch two network for small object extraction can also be preserved, as shown in the orange oval marker in the figure, which is a road under the true color image. The result of using the branch one GAT for classification is shown in Figure 11b. Almost no road information is detected. The result of using the branch two residual network for classification is shown in Figure 11c, the road information can be extracted more completely, and after fusing the two branches, a smooth and more complete result is obtained. In addition, as shown in the green rounded rectangle marker, the small building in Figure 11(b), the result is too smooth, and the classification of small buildings in Figure 11(c) produces more noise. The fusion of the two is more consistent with the real ground land use type, and the overall accuracy (%) and Kappa (×100) coefficients are above 95 using the network proposed in this paper to classify Landsat images of Sanya City (2000 and 2020). The use of this network by manually plotting samples is operationally convenient and provides a highly accurate land use type classification result, which provides reliable data for subsequent habitat quality assessment.

A comparison of classification results in a zoomed-in typical region is shown in Figure 12, where the SVM algorithm exhibits a large number of noise points; 2D-CNN and Resnet perform better but show some obvious misclassification of forest and built-up land; GCN generates a very smooth result, with obvious misclassification of built-up land and arable land; and the results of the FGR model are closest to the real results.

Figure 12. Ground truth and classification maps obtained by different methods on a typical region.

3.4. The spatiotemporal pattern changes in habitat quality

Combining the high precision LUCC we obtained using FGR-net with the previously collected data, we obtained the HQ of Sanya City according to the InVEST model, and the HQ indices were divided into five intervals and classified into five classes of Ⅰ, Ⅱ, Ⅲ, Ⅳ, and V accordingly (Table 5).

Table 5. HQ grading in Sanya City from 2000 to 2020.

The level of HQ	Grade	Score interval	Proportion in 2000	Proportion in 2020
Low	Ⅰ	0–0.2	0.89%	2.66%
Relatively low	Ⅱ	0.2–0.4	23.90%	29.37%
Medium	Ⅲ	0.4–0.6	0.01%	0.01%
Relatively high	Ⅳ	0.6–0.8	1.21%	1.15%
High	V	0.8–1	73.99%	66.81%

From the time scale, the overall HQ in Sanya shows a declining trend from 2000 to 2020, and the mean value of HQ reduces from 0.82 (in 2000) to 0.76 (in 2020), with a decline rate of 6.88%. Class V and Class II are the dominant types of HQ in Sanya, with Class V HQ areas account for approximately 70% (Table 5), mainly due to the good natural environment in Sanya, emphasis on ecological protection, and large woodland area. The study area had the lowest proportion of habitat quality belonging to class III, and the proportion of the area of class I and class II (i.e. low-grade habitat quality areas) increased, while the proportion of the area of class IV and class V (i.e. high-grade habitat quality areas) decreased between 2000 and 2020, with class V decreasing the most, by 7.18%. This indicates that the high-grade habitat is continuously transforming into a low-grade habitat, and Sanya City has potential habitat degradation risk. The reason for the decrease in Class V areas and the increase in Class II areas is that the accelerated urbanization of Sanya City has led to considerable forest land, grassland, and farmland being occupied by construction land on the one hand. On the other hand, deforestation and reclamation still exist in some areas, as well as reclamation around reservoirs, ponds, and beaches. Class III and Class IV areas are basically unchanged because the policy of returning farmland to forest (grass) and field to lake (wet) is promoted.

The spatial distribution of HQ in Sanya is shown in Figure 13. Figure 14 shows the changes in Sanya's HQ. It is obvious that the areas with decreased HQ are significantly more than those with improved HQ and are concentrated in the southern coastal area. The Sanya City Land Use Master Plan (2006–2020) clearly highlighted the strategy of land use (TPGSC 2010): to give priority to guaranteeing the land for the construction of Sanya central town, and to meet the land for the key projects of coastal tourism, new industry, cultural industry, port and transport in the east and west flanks areas of Haitang Bay, Sanya Bay and Yacheng Bay; to adjust the industrial structure with a high starting point, and to prompt highly efficient land-use from the coast to the inner land step by step. This resulted in the rapid development of the tourism industry and the industry in the coastal areas, and the emergence of a large amount of land for building purposes. Sanya's HQ has a serious degradation trend. However, there are a small number of areas in the north of Sanya City where the quality of habitats has been improved, which is because these areas have actively responded to the Sanya City Land Use Master Plan (2006–2020), which stabilised the area of cropland during the planning period, increased the area of forest and other agricultural land, and implemented the afforestation of barren slopes. The closure of mountains to forests and the return of land to forests to a moderate extent in the areas of high-altitude mountains and areas along the rivers and reservoirs is promoted to maintain a good ecological environment.

Figure 13. Spatial distribution of HQ in Sanya City in 2000 and 2020.

Figure 14. Changes in HQ in Sanya from 2000 to 2020.

3.5. Correlation analysis of habitat quality in sanya on various factors

There were some differences in the strength and direction of the effects of the changes in natural and economic factors on the spatial patterns of HQ. Correlation analysis was performed using the SPSS software, and the correlation coefficients between the HQ and each influencing factor passed the significance level test of 0.01 in both years (Figure 15). HQ was significantly negatively correlated with TEM, POP, NTL, and GDP, whereas the remaining indices were significantly positively correlated. Among the socioeconomic indicators, POP was significantly negatively correlated with HQ in 2000 (correlation coefficient of −0.424), and the negative correlation was stronger than that of NTL (correlation coefficient of −0.376) and GDP (correlation coefficient of −0.416); NTL had the strongest negative correlation with HQ in 2020. Among the natural factor indices, NDVI had the highest correlation coefficient with HQ.

Figure 15. The correlation between HQ and influencing factors in Sanya City. PRE is the annual precipitation; TEM is the annual average temperature; POP represents population density; NTL is nighttime lighting; DEM is the elevation; GDP is the gross domestic product; NDVI is the normalized differential vegetation index; and the size of the circle denotes the magnitude of the correlation coefficient.

The study area was divided into 100 × 100 grid cells using the equally spaced sampling method, and 4032 sample points were obtained after removing error points. Habitat quality indices in the grid were used as the explanatory variables, and natural factors and socio-economic factors were used as explanatory variables. From the results, it is clear that the geographically weighted regression model fits well (Table 7).

To compare the impact of each factor on HQ, ordinary least squares (OLS) was used to quantify the influence of natural and social factors on regional HQ. The results showed that the variance inflation factor values of each variable were less than 7.5 (Table 6), indicating that there was no redundancy between the selected variables. To facilitate a comparative analysis, the core factors affecting habitat quality in both periods were selected based on the findings of the OLS model, and GWR was applied (Ji et al. 2022) with the MGWR model to calculate the relevant indicators (Table 7). Compared with the classical GWR model, the AICc value of the MGWR model was smaller, and the R² values of the goodness of fit of the model in 2000 and 2020 were 0.552 and 0.639, respectively, which were higher than those of the GWR models. This indicates that the MGWR model has better fitting results, better explains the local variations in the influencing factors, and reduces the spatial autocorrelation of the model residuals.

Table 6. Results of variance inflation factor.

Variables	DEM	TEM	PRE	NDVI	GDP	POP	NTL
2000VIF	4.781	5.028	1.180	1.845	4.143	3.998	1.175
2020VIF	6.660	6.043	1.934	2.297	4.383	3.265	2.301

Table 7. Indicators of the model diagnosis.

Model	AICc	R^2	Adjusted R^2
2000gwr	2390.236	0.536	0.508
2000mgwr	2379.051	0.552	0.514
2020gwr	2313.309	0.612	0.583
2020mgwr	2296.709	0.639	0.596

The MGWR model was used to analyze the spatial heterogeneity of the core factors affecting habitat quality. The socioeconomic indicators, NTL and POP, showed significant differences in HQ in different regions (Figure 16).

Figure 16. Spatial distribution pattern of correlation coefficients between HQ and socio-economic factors, night-time lighting (NTL) and population density (POP), in Sanya.

The 2000 NTL had a significant negative effect on HQ, with the proportion of negative values reaching 98.43% of the regression coefficient and only 17 positive cases in a significant quantity of red areas in the figure. The negative impact of NTL in 2020 weakened and the positive impact strengthened in some regions, with the proportion of negative regression coefficients decreasing and the proportion of positive values increasing, with small differences in the values taken by the regions.

The negative effect of POP on HQ was dominated by the negative effect, which increased in 2020 and exhibited a negative effect in all regions. POP had a large negative driving effect on habitat quality in 2000, with the proportion of negative regression coefficients reaching 70.66%; such effects gradually increased in size and regional extent by 2020. The negative impact of POP on HQ is greater than that of NTL, as shown by the fact that the overall correlation coefficient in 2000 is smaller than that of NTL, and the regional POP correlation coefficients are all around −0.7 in 2020. In contrast, the negative impact of NTL varies greatly from region to region, and there are a few positive impact regions.

Regional economic development affects regional HQ to a certain extent, and with the increase in regional income and GDP per capita, people's demand for environmental quality increases, and the government has more funds to invest in environmental management and protection. NTL has a certain enhancement effect on habitat quality in some regions (Zhao et al. 2022). The main manifestation is that NTL and HQ show a negative correlation in the northern part of Sanya city, which accounts for 80% of the total area. However, the central and southern coastal areas of Sanya City showed a weak positive correlation between NTL and HQ (red and orange areas are shown in Fig). Most of the northern part of Sanya City is forested, and the negative correlation coefficient indicates that increased human activities in forested areas will decrease habitat quality. The southern coastal area shows a slightly positive correlation coefficient, which indicates that the increased intensity of human activities in the built-up areas of the southern coast will not decrease habitat quality.

The correlation between the HQ and natural factors showed regional variability (Figure 17). In 2000–2020, the natural factors TEM and PRE had positive and negative ‘bidirectional’ effects on HQ in different regions. In 2000, in all regions, PRE showed a negative correlation with habitat quality, but the absolute value of the correlation coefficient was small, indicating that the negative effect was weak, and the grid of the positive effect increased with time. In 2020, in all regions, PRE and habitat quality showed a positive correlation, with a correlation coefficient of approximately 0.3, and the maximum positive area of the regression coefficient with HQ was mainly distributed in the northern mountainous areas where there was relatively little anthropogenic disturbance. In terms of TEM, the negative effects in 2000 were mainly distributed along the southern coast and its surrounding areas, and the number of grids with negative effects further shrank over time, from 60% (2000) to 0 (2020). Most areas along the river cities showed strong positive effects, and the distribution pattern of correlation coefficients changed from north–south graded to east–west graded. DEM and NDVI had positive effects on HQ, and the positive correlation between DEM and HQ has increased over time. However, the overall distribution pattern of the correlation coefficients has not changed significantly (although the correlation coefficient between DEM and HQ has a negative value in 2020, the overall proportion does not exceed 8%), and the positive effect of NDVI on habitat quality has weakened, showing an overall positive effect of NDVI on HQ. The positive influence of NDVI on HQ weakened, as shown by the decrease in the overall correlation coefficient, and its distribution pattern did not change significantly, showing a roughly east–west graded distribution pattern. The regression coefficients were not significantly different between groups.

Figure 17. Spatial distribution pattern of correlation coefficients between HQ and natural factor indicators in Sanya.

In terms of natural factors in 2020, NDVI, DEM (excluding less than 8% negative values), TM, PRE, and HQ correlation coefficients are positive in all areas of Sanya city, where positive correlation is graded by regional gradients and TM and PRE have similar spatial distribution patterns (all graded east–west and gradually decreasing from east to west). Compared to social factors (nighttime lighting data), the impact of natural factors on HQ shows a large spatial variation. However, the mean absolute values of the regression coefficients indicate that the natural factors have a stronger association with HQ, and the strength of the association between the four natural factors and HQ is from large to small (NDVI > DEM > TM > PRE). This result is consistent with the findings of Xie and Zhang (2023). Xie detected differences in the explanatory power (q-values) of various factors affecting HQ changes using a geographic probe; the q-values were in the order of land use type > NDVI > slope > temperature > rainfall > population density.

4. Discussion

In this paper, we constructed a two-branch classification network based on Landsat images of Sanya City and designed several intuitive comparison experiments to verify the superiority of the FGR-Net, which proved that the model could well combine the features of CNN and GNN. After fusing the two features, the classification results have less point-like noise than CNN and higher classification accuracy than GNN, and its classification results are more in line with the real ground cover, which provides good data for the subsequent experiments. The findings of Gao et al. (2021) supported this result. They designed a CNN branch using channel attention and a GCN branch, and finally spliced the features of the two branches to obtain good classification results, which proved the complementarity of the CNN and GNN. The CNN branch of FGR-net uses a well-designed attention network and two improved residual units. Although the improvement step is simple, it results in faster convergence of the model. In contrast to some advanced CNN classification networks (Gao et al. 2022; Gao et al. 2023a; Zhang et al. 2023a; Zhang et al. 2023b), such as DFInet, FGR-net uses an attention network to generate weights for channel and spatial information, and data information is focused on when it has a high weight in both space and channel. In addition, it uses GNN branching to suppress noise and ultimately obtains a more realistic result, which is difficult to achieve using a CNN alone. The proposed model also has some limitations, the FGR-net needs to set too many parameters, such as the number of hyperpixels after performing hyperpixel segmentation, and the weighted fusion coefficient of the two branch networks, if these two numbers are not set reasonably, it may lead to the GCN branch obtaining an overly smooth image, and this reveals that future work should reduce the parameter settings, so that the model extracts and fuses the features in an adaptive way; it has a high requirement for the original labels to be correct, and correctly labeling remote sensing images will cost a lot of human and material resources, and in order to label the existing images, the researchers have done a lot of work, but due to a variety of reasons, the labeled images are not necessarily all correct, call this labeling error as category noise, SVM shows a certain resistance to category noise, but when the category noise reaches a certain amount, the accuracy of SVM will also decline, however, the ability of FGR-net to resist category noise is poor, but FGR-net performs well in small sample tasks.

The spatial distribution pattern of HQ is influenced by physical geographic factors and socioeconomic activities, and land use change is considered a significant factor in biodiversity change. The results showed that when Sanya City was divided into three regions from north to south, the overall habitat quality of Sanya City showed a spatial distribution pattern of high in the upper and lower regions and low in the middle region, with a clear hierarchical structure that was highly compatible with the topographic characteristics of the city. There is a correlation between HQ distribution and topographic conditions, and natural elements and socio-economic factors also affect the spatial pattern of regional HQ. Generally speaking, areas with higher topography have less human activity, whereas areas with lower topography have intense human activities. In the context of this thesis, the correlation of HQ with topographic and climatic factors was explored using the MGWR model. The findings indicated that HQ was positively correlated with topographic factors, which is similar to the study by Jia et al., who concluded that habitat quality has a significant topographic gradient effect (Jia et al. 2022). Assessing the habitat quality of the three major watersheds in Hainan Island, the results showed that the overall HQ of the three major watersheds showed a spatial distribution pattern of high in the upstream source area and low in the middle and downstream areas, with an obvious hierarchical structure. The distribution of HQ was correlated with topographic conditions, and the low-value areas were mostly areas with lower elevations. First, because the complex topography of Sanya city directly affects the spatial differences of water and heat and other climatic conditions, resulting in the differences of terrestrial vegetation zonal distribution and the diversity it exhibits, the concentrated distribution of woodlands in the upper and lower higher elevation areas supports the excellent habitat quality and biodiversity level. In contrast, in the lower elevation areas of Sanya city, especially in the coastal areas, population density, land use intensity, and other intensive socio-economic activities have led to the gradual expansion of construction land in Sanya from coastal areas like inland. This has posed a great threat to habitat quality and led to the loss and degradation of local habitats, and with the future demand for natural ecological space for infrastructure construction and regional development in Sanya, there will still be considerable growth in the amount of construction land in Sanya in the future, leading to a more severe test of the ecological environment.

The findings of the correlation analysis revealed that POP and NTL had a significant negative effect on HQ, while NDVI and DEM were positively correlated with HQ. The spatial autocorrelation of HQ was considered, and the spatial heterogeneity of natural factors and socioeconomic indices on HQ was comprehensively explored based on a combination of OLS and MGWR models, which improved the scientific validity of the study to a certain extent. Natural factors and socio-economic indices have positive and negative ‘bidirectional’ effects on HQ. The spatial distribution of the regression coefficients showed that the effects of POP and NTL on habitat quality were significantly spatially heterogeneous, which is consistent with the findings of related studies (Zhu et al. 2020). The effect of POP on HQ was negative, with a greater degree of negative effect than that of NTL. The negative effect increased with time, while the positive effect of NTL on HQ gradually increased in some areas. Urban expansion within the framework of rapid urbanization changes regional patterns, destroys ecosystem integrity, exacerbates environmental pollution, and indirectly leads to a decline in HQ (Wang et al. 2020). Some studies have further demonstrated that urban sprawl due to rapid urbanization has significantly altered ecological landscape patterns in densely populated areas of the region, posing a serious threat to habitat quality (Liang et al. 2022; Wei et al. 2022; Wei et al. 2023). Similar studies have focused on areas with high-intensity human activities, such as metropolitan areas and population concentrations (Wang et al. 2020; Wei et al. 2023). The impact of NTL on HQ exhibited different correlations and intensities at spatial and temporal scales, and the high values of negative NTL on habitat quality were mainly located in the northern part of Sanya City, which may be related to local vegetation types, with land use types of mainly water or forest grassland. In some areas, owing to the location advantage, a large number of industries have been relocated, towns have been developed significantly owing to the influence of tourism and other development industries, arable land and construction land have encroached on the ecological land, causing great disturbance to the ecology of the area. It is urgent to strengthen the ecological and environmental protection and management.

NTL characterizes the indirect disturbance effect of urban socioeconomic development and the intensity of human activities on ecosystems, and high NTL and high habitat quality are not contradictory to each other. When cities develop to a certain level, ecological and environmental protection will attract significant attention from local governments, thus alleviating habitat pressure (Sallustio et al. 2017; Wang et al. 2022c). Urban-related policies and development strategies can affect habitat quality to a certain extent, and some areas can achieve synergistic economic and ecological development (Tang et al. 2021). NTL and habitat quality have evolved over time. As time evolves, NTL and habitat quality show positive correlation effects in some areas, with higher regression coefficients in local areas such as the central and southern regions, which is because some areas were originally land types, mainly wasteland/farmland, etc., with low habitat quality indices. In recent years, influenced by the project of returning farmland to the forest (grass) and ecological protection policies, the government planted trees on a large scale to enhance green spaces and maintain the ecosystem in the region. The strength and correlation properties of natural factors (DEM, PRE, TEM, and NDVI) on habitat quality differed somewhat at the spatial and temporal scales, which is more consistent with the findings of related studies (Yan et al. 2018). The same natural factor indices show different correlation properties and strengths with habitat quality depending on the region, and ecological regulation in the process of habitat quality optimization should focus on regional differences in the correlation analysis between natural factor indices and habitat quality (Li et al. 2021). Therefore, attention should be paid to natural environmental protection in socioeconomic development and construction activities to optimize production, living, and ecological space and balance the relationship between urbanization and the ecological environment. The scale of construction land should be reasonably controlled, green ecological spaces should be constructed in construction land-intensive areas, ecological corridors should be built, ecological security patterns should be optimized, and ecological environmental protection should be increased to achieve green and high-quality development.

5. Conclusion

1. The proposed two-branch model (FGR) has a high classification accuracy on the dataset used and is superior to traditional machine learning algorithms, CNN networks, and GNN networks in terms of its ability to accurately identify ground cover types.

2. In the habitat quality change, V HQ area and II HQ area are the main types, and most of the high-grade habitat quality areas are located in forest grassland. Most of the low-grade habitat quality areas are located in cultivated land for construction, and there is a tendency that V habitat area is transformed to II habitat quality area with time.

3. Over the past 20 years, the overall HQ of Sanya City has been high and shows obvious spatial and temporal heterogeneity, with a slight degradation trend in time and a spatial evolution pattern with the built-up areas of large and medium-sized cities and small towns as low habitat quality areas, followed by a gradual increase in agricultural production areas and ecological functional areas.

4. The MGWR model fits better than the classical GWR model and reveals different relationships between various urbanization parameters and habitat quality. The HQ was significantly influenced by rapid urbanization and natural factors, and there was significant spatial heterogeneity in the influence of different factors. Natural factor indicators (PRE, TEM, DEM, and NDVI) had negative or positive effects on the HQ in different regions. In future regional economic development, more attention should be paid to coordination with ecological protection to maximize comprehensive social, economic, and ecological benefits.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

This study analyzed publicly available datasets, which can be found here: Geospatial Data Cloud (https://www.gscloud.cn/) and the Resource and Environmental Science and Data Centre of the Chinese Academy of Sciences (https://www.resdc.cn/).

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (NSFC) [grant number 42341206], Chunhui Program Cooperative Research Project of Chinese Ministry of Education [HZKY20220279], Henan Provincial Science and Technology Research Project (232102211019, 222102210131), the Key Research Project Fund of Institution of Higher Education in Henan Province (23A520029), and Japan Society for the Promotion of Science (JSPS) KAKENHI [grant number 20K12146].