Scale effect of coastal landscape pattern stability and driving forces: a case study of Guangdong Province, China

1. Introduction

Coastal zones are key transition regions between marine and land environments, supporting 40% of the world’s species (Wetlands International, 2021). They are also sites of high carbon accumulation, with rates as high as (210 ± 20) g/m² (Costanza et al., 1997). The ecological functions of coastal wetlands include climate regulation, water purification, and wave protection (Halder et al., 2021; Sandilyan and Kathiresan, 2015), and they form a natural barrier against flooding and storm surge erosion in estuarine areas. Coastal zones support approximately 60% of the global population and two thirds of world’s large and medium-sized cities (Syvitski et al., 2009). However, with the advent of global warming and the intensification of anthropogenic disturbance over the past century, coastal zones globally have become threatened by degradation or are already seriously degraded (Tessler et al., 2015). Only 15.5% of coastal regions worldwide remain intact, with the remainder severely impacted by fishing, agricultural development, and urbanization (Williams et al., 2022). Coastal zones are threatened by excessive reclamation, erosion retreat, functional degradation, and resource decay (Kirwan et al., 2013; Jiang et al., 2015), which in turn pose a serious threat to the regional ecological balance. With the increasingly prominent degradation of coastal ecosystems (Sun et al., 2017), there has been a growing need for studies on coastal landscape pattern stability (Zhao et al., 2012; Liu and Ma, 2017). Studies conducted thus far have formed the scientific basis for ensuring the safety of estuarine and coastal levees as well as species protection. They have also contributed to the formulation of the trinity (quantity, quality, and ecology) protection policy for coastal wetlands.

Coastal landscape patterns are inherently spatially correlated and scale-dependent (Wu, 2007; Chi et al., 2019), exhibiting a ubiquitous spatial heterogeneity that is critical for both their structure and function across various scales (Wu, 2004). Different landscape patterns can be revealed by different spatial scales of observation (Hanke et al., 2014). The strength of the relationship between land use and coastal water conditions depends upon the spatial scale (Zhou et al., 2023). Even within the same region, landscape features at different scales can be highly heterogeneous (Wickham et al., 1997; Cui et al., 2016). Relevant studies include observing the effects of changing the grain on observed spatial patterns (Turner et al., 1989); examining how pattern indices change with scale in real landscape of different kinds (Wu et al., 2002); emphasizing the significance of multiscale analysis of landscape composition and configuration at both class and landscape levels (Wu, 2004). Moreover, other studies have shown that the stability of the coastal landscape is characterized by spatiotemporal variation, with specific stability characteristics manifesting at different spatiotemporal scales (Luo et al., 2004; Xie et al., 2005). However, research into the scale-dependent effects on landscape patterns, especially investigations into the optimal scale threshold for stability of coastal landscape patterns at the provincial level, remains relatively sparse (Oehri et al., 2017). Therefore, studies of the evolution of coastal landscape stability patterns and their mechanisms require multi-scale hierarchical analyses as a prerequisite (Cushman and McGarigal, 2002; Lv et al., 2007), both to gain an in-depth understanding of ecological processes in coastal zones and as a key premise for correctly understanding pattern-process mechanisms.

Recent studies of coastal zones have mainly focused on landscape structure changes and environmental effects (Hauser et al., 2017), biomass estimation (Ghosh et al., 2021), and ecological service function and risk assessment (Trégarot et al., 2021; Wang et al., 2021), whereas few studies have investigated coastal landscape pattern stability and its evolutionary mechanisms (Wu and Zhang, 2021). Xu et al. (2018) reported that the area of extremely stable and extremely unstable landscape in the Beijing–Tianjin–Hebei urban agglomeration continued to expand during 1980 to 2010. Tang et al. (2019) pointed out that over past four decades, the coastline stability index of Liaoning Province had gradually decreased by 0.89 due to the combined effects of natural and human influences. Moreover, Tian et al. (2021) reported that human land use behaviors could notably affect the coastal stability along the southern coast of Laizhou Bay. However, most of the previous studies on landscape pattern stability have been focused on small-scale area (i.e., cities, urban agglomeration and estuaries), and they qualitatively or semi-quantitatively explored the factors influencing landscape pattern stability changes. In addition, Chen et al. (2023) suggested that aquaculture land, arable land and hydraulic engineering land had negative effects on the stability of coastal wetlands in Jiangsu Province, while salt industry land had positive effects. Therefore, comprehensively understanding the dynamic processes of coastal landscape pattern stability in medium to large scales and quantifying their driving mechanisms are of great significance for the management and protection of coastal ecosystems.

Guangdong Province has 4 114 km of mainland coastline, accounting for approximately 1/4 of the length of China’s mainland coastline (People’s Government of Guangdong Province, 2017). Guangdong Province is regarded as a pioneer region of China’s economic reform, and is undergoing rapid urbanization (Cheng et al., 2022). Thus, the coastal landscape of Guangdong Province is among the most dynamic regions globally; its Zhujiang River (Pearl River) Estuary is the most concentrated area of economic and human activities in China and possibly worldwide (Wu et al., 2018). Due to large-scale reclamation that has accompanied urbanization, the tidal flat area of the Lingdingyang Estuary was reduced by about 200 km² during 1955 and 2015 (Chen et al., 2020). Rapid urbanization has also induced buildup in rural areas that have expanded irregularly, causing fragmentation of the coastal landscape of the Hanjiang Delta (Qiu et al., 2022). The ecological risk to the landscape in the coastal zone of Guangdong Province has gradually increased over the past 15 years (Cheng et al., 2022), with climate change and human activities as the major factors driving the shrinkage of coastal wetlands in the Zhujiang River Delta (Zhou et al., 2019). However, the coastal landscape pattern stability evolution process in Guangdong Province influenced by change in natural factors and rapid urbanization is still not well understood, and its quantitative driving mechanism remains unclear. Thus, more relevant studies need to be further explored.

This study explored coastal landscape pattern stability and its underlying mechanisms in Guangdong Province against a background of global climate change and frequent anthropogenic disturbances in the context of China’s strategic goals of integrated land and sea spatial planning and coastal zone sustainability. The main goals of this study were to further explore the scale effect of coastal landscape pattern stability to determine the optimal scale for analyses, to analyze the spatiotemporal evolution of coastal landscape pattern stability at the optimal analytical scale, and to quantitatively explore the mechanisms driving landscape pattern stability evolution. Our findings will provide a better understanding of the long-term dynamics of coastal landscape evolution and scientific support for the high-quality development of the coastal zone in Guangdong Province.

2. Study area and materials

2.1 Study area

Guangdong Province (20°09′–25°31′N, 109°45′–117°20′E), located in the southernmost region of Chinese mainland (Fig. 1a), covers an area of 1.8 × 10⁵ km². Since 1989, it has become the most developed province in China (Guo et al., 2022). The mainland coastline of Guangdong Province is the longest in China, at 4.1 × 10³ km (https://www.gd.gov.cn/gkmlpt/content/0/146/post_146486.html#7), and includes extensive beaches and numerous other sea landscapes (Chen et al., 2023). The climate of Guangdong Province is subtropical marine monsoon, and thus hot and rainy in summer and warm and humid in winter. The average annual temperature is 21.2–22.5℃, and average annual rainfall ranges from 1700 mm to 2500 mm. The terrain of Guangdong Province is generally high in the north and low in the south, and mainly consists of mountains, hills, and alluvial plains.

Figure  1.  Location of the coastal zone in Guangdong Province, China. ZJ: Zhanjiang City; MM: Maoming City; YJ: Yanjiang City; JM: Jiangmen City; ZH: Zhuhai City; ZS: Zhongshan City; GZ: Guangzhou City; DG: Dongguan City; SZ: Shenzhen City; HZ: Huizhou City; SW: Shanwei City; JY: Jieyang City; ST: Shantou City; CZ: Chaozhou City. Drawing from Minstry of Natural Resources of China (http://bzdt.eh.mnr.gov.cn/). Approval number: GS(2020)no.4619.

DownLoad: Full-Size Img PowerPoint

The coastal zone is divided into 41 country districts that account for 27.67% of the total mainland area of Guangdong Province, covering an area of 4.98 × 10⁴ km², and is home to 47.63% of the permanent population of the province, which contributed 56.41% of the province’s total gross domestic province in 2021 (Guangdong Provincial Bureau of Statistics and Survey Office of the National Bureau of Statistics in Guangdong, 2022). Over the past few decades, rapid urbanization has induced a high population concentration along the coast (Zhu et al., 2018) and disorderly urban land expansion (Lin et al., 2020), which together with intensive tideland reclamation (Chen et al., 2020) have led to highly complicated coastal landscape patterns. Coupled with frequent natural disasters such as typhoons, floods, and droughts (Guo et al., 2022), the trend in coastal landscape evolution is highly uncertain.

2.2 Materials

Recent developments in Earth observation technologies have provided archives of satellite images. Landsat is the longest-running set of operational platforms collecting images of the Earth’s surface. Landsat-4, -5, -7, and -8, launched between 1982 and 2013, operated the thematic mapper (TM), enhanced thematic mapper (ETM), enhanced thematic mapper plus (ETM+), and operational land imager (OLI), respectively. These sensors have detected the changes that occurred since the 1980s, providing multiple remote sensing images (119–124 paths, 43–46 rows) from Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI, which are available in the historic satellite data archive (HSDA) within Google Earth Engine (GEE). Among other applications, these images have been used to monitor the dynamic evolution of coastal landscape pattern stability since the 1980s. Annual image availability for the coastal zone of Guangdong Province is currently approximately 15 to 16 images, with a 30-m spatial resolution and 16-d revisit time (Table 1). In this study, 125 images were selected for processing based on the world geodetic system 1984 (National Geospatial-Intelligence Agency) geographical coordinate system and universal transverse mercator projection coordinate system. In general, the selection of images was based on the following principles: (1) there was minimal cloud cover (less than 5%); (2) most of the acquired images were taken at low tide levels (ranging from –0.28 m to 0.21 m), so that all mangrove canopies were exposed above the seawater; (3) images were selected from June to September in each typical year, so that all vegetation was in the growing season.

Table  1.  Remote sensing data used in this study

Year	Data	Sensor	Number of scenes	Month of acquisition	Cloud cover	Spatial resolution and revisit time	Path and row
1985	Landsat-5 TM/ETM+	thematic mapper/enhanced thematic mapper plus	16	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
1990	Landsat-5 TM/ETM+	thematic mapper/enhanced thematic mapper plus	16	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
1995	Landsat-5 TM/ETM+	thematic mapper/enhanced thematic mapper plus	15	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
2000	Landsat-5 TM/ETM+	thematic mapper/enhanced thematic mapper plus	16	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
2005	Landsat-5 TM/ETM+	thematic mapper/enhanced thematic mapper plus	16	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
2010	Landsat-5 TM/ETM+	thematic mapper/enhanced thematic mapper plus	16	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
2015	Landsat-8 OLI	operational land imager	15	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
2020	Landsat-8 OLI	operational land imager	15	from June to September in each typical year	<5%	30 m × 30 m, 16 d	paths: 119–124, rows: 43–46
Note: Total number of scenes is 125.

 | Show Table

DownLoad: CSV

3. Methodology

3.1 Image processing

GEE is a cloud-based geoprocessing platform that provides vast amounts of publicly available geospatial datasets. It offers a multi-petabyte catalog of satellite imagery, including Landsat imagery, and facilitates fast data access and parallel processing for rapid analysis (Gorelick et al., 2017; Kennedy et al., 2018). In this study, we used atmospherically corrected United States geological survey landsat surface reflectance products to retrieve the distributions of land use types from 1985 to 2020 based on the GEE platform. We used the random forest algorithm to classify ground objects, as it is a computationally efficient and highly effective supervised learning algorithm. The theory of random forest model is to build a group of decision trees to classify the land cover according to the input data. This method works by picking features from the data randomly and using them to make the best possible divisions within the trees. The success of these divisions is measured by something called the Gini coefficient. Essentially, the lower the Gini coefficient, the better the division is at accurately classifying the land (Chen et al, 2023), as

where n is the number of categories generated by the spilt, j is category j generated by the spilt, k is the number of the types of land use cover, i means the land use cover i and p is the frequency of occurrence of land use cover in a category.

Refer to previous studies in this study area (Li et al., 2010; Chen et al., 2023b; Zhang and Lang, 2022) and considering the uniqueness and complexity of coastal landscape (Phillips, 2018; Chen et al., 2024), the prepared images were classified into six land cover categories: impervious water surface, water, forest, farmland, mangrove forest, and other. The above processing involved five steps as follows.

(1) The image dataset with few clouds in the study area was constructed. The GEE platform was used to select Landsat-5 (TM), Landsat-7 (ETM+) and Landsat-8 (OLI) images in each typical year (1985, 1990, 1995, 2000, 2005, 2010, 2015, and 2020). The CloudScoreRange was then set to 5 and images meeting set conditions were spliced to generate images with less clouds in each typical year for the study area.

(2) The region of interest (ROI) for coastal landscape classification was built. Google Earth Pro was used to randomly generate 1 800 sample points in the image of the study area for each year (300 sample points for each land cover type). These sample points were imported into the GEE application programming interface (API) for further analysis. The samples were then randomly divided into training and validation datasets at a 7:3 ratio.

(3) A classification decision tree was further created and landscape types were rapidly extracted. Some previous studies have shown that the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), mangrove vegetation index (MVI), and normalized difference water index (NDWI) were helpful and widely applied to extract complex landscape information in coastal zones (Baloloy et al., 2020; Karaman, 2021). The detailed calculation method and physical significance of four spectral indices were consulted in relevant studies (Baloloy et al., 2020; Zheng et al., 2021). Additionally, based on previous related studies (Baloloy et al., 2020; Karaman, 2021; Davis et al., 2023), and considering the uniqueness of the coastal areas in Guangdong Province, a sensitivity analysis was conducted on four indices to determine the optimal threshold values for the indices in the study area (Fig. 2). Firstly, the NDWI was used to identify the water and non-water areas and the threshold was set to 0 (Fig. 2a); secondly, the NDVI and EVI were applied to extract the vegetation area from the non-water area, and the threshold values were set to 0.05 (Fig. 2a), respectively; thirdly, the MVI was used to extracted the mangrove forest from vegetation area, and the optimal minimum threshold was set to 4.5 (Fig. 2b). Finally, these decision tree rules were input into the random forest classifier for fast classification.

Figure  2.  Sensitivity analysis of four indices to classification accuracy. NDWI: normalized difference water index; NDVI: normalized difference vegetation index; EVI: enhanced vegetation index; MVI: mangrove vegetation index.

DownLoad: Full-Size Img PowerPoint

(4) Post-classification processing was applied to the fast classification results. A manual visual interpretation approach was used to correct the image classification results for mangroves within a 10-km buffer zone extending inland from the coastline. This procedure effectively corrected misclassifications between forest and mangroves resulting from “spectral confusion”.

(5) A classification confusion matrix was constructed for accuracy evaluation. The user’s accuracy, total accuracy and Kappa coefficients for the landscape classification results in typical year were obtained in Table 2. Among them, the total accuracy ranged from 83.06% to 90.17%. Moreover, the average Kappa coefficient exceeded 0.86 (Table 2). Therefore, the accuracy of the classification was sufficient to analyze changes in landscape pattern stability along the coast of Guangdong Province and the influencing factors.

Table  2.  User’s accuracy, total accuracy and Kappa coefficients for the landscape classification result

Calendar (year)	User’s accuracy/%						Total accuracy/%	Kappa coefficient
	Impervious water surface	Water	Forest	Farmland	Mangrove forest	Other
1985	83.33	85.67	85.00	82.67	81.33	80.33	83.06	0.826
1990	84.67	87.00	86.00	83.33	82.33	82.00	84.22	0.842
1995	86.00	88.33	87.00	85.00	84.00	82.33	85.41	0.851
2000	87.00	89.33	88.67	87.00	84.67	82.67	86.56	0.854
2005	88.33	90.67	90.00	87.00	86.67	84.67	87.89	0.865
2010	89.33	91.33	90.67	88.33	87.00	85.33	88.67	0.876
2015	91.00	92.33	91.33	89.00	87.67	86.00	89.56	0.888
2020	91.67	93.00	92.00	89.33	88.33	86.67	90.17	0.898

 | Show Table

DownLoad: CSV

3.2 Landscape index selection

The number of landscape types, area and shape of patches, degrees of aggregation and segmentation, and landscape diversity were used to select six commonly used landscape metrics: patch density (PD), which measures the number of patches per km²; total edge contrast index (TECI), which equals the sum of the lengths (m) of each edge segment in the landscape multiplied by the corresponding contrast weight, divided by the total length (m) of edge in the landscape, multiplied by 100 (to convert to a percentage); mean patch size (MPS), which measures the average area (hm²) of all patches; mean proximity index (MPI), which measures the proximity and fragmentation of each landscape type; aggregation index (AI), which measures the degree of aggregation for each patch type; contagion (CONTAG), an information theory-based index that measures the spatial aggregation of each patch type; and Shannon’s diversity index (SHDI), which measures patch type diversity. Detailed descriptions of the equations used to calculate each index and their ecological significance are provided in the help files of the Fragstats 4.2 software and a related study (You et al., 2006). Generally, the higher PD and TECI value indicates a greater degree of landscape fragmentation. For MSP, MPI, AI, and CONTA, a more clustered landscape indicates a smaller degree of fragmentation. PD is used to compare the degrees of fragmentation of various landscape types at various levels, but is applicable only for horizontal comparisons (Wang et al., 2014), and such comparisons require additional indicators such as area and proximity. A higher SHDI value indicates a more even distribution of patch types in the landscape.

Generally, the dramatic change of landscape metrics can indicate the scale sensitivity of coastal pattern, especially the inflection point. In addition, Moran Ⅰ index was used to quantitatively explore the spatial autocorrelation of coastal landscape pattern stability. The larger the Moran Ⅰ index value, the more obvious the spatial correlation and the more comprehensive the regional information. Therefore, the inflection point of landscape metrics curve and the peak value of Moran Ⅰ index were the important principles of optimal scale selection in this study.

3.3 Coastal landscape pattern stability model

Landscape metrics can highly concentrate landscape pattern information and quantitatively reflect the structural composition and spatial configuration of the landscape (Wu, 2007). Therefore, given the complexity of the coastal landscape, in this study PD, CONTAG, and TECI were used to construct a model for measuring landscape pattern stability (S) (Xu et al., 2018) as follows:

where S is the landscape pattern stability metric, C is the contagion metric, P is the patch density metric, and T is the total edge contrast index. The larger the S value, the higher the landscape pattern stability. On the contrary, the landscape pattern stability is lower.

3.4 Factor detector

The geodetector spatial variance analysis method (Wang and Xu, 2017) was used in this study to calculate spatial variance and interpret its mechanisms. We imported the GD package into RStudio (R Core Team, Vienna, Austria) to calculate the q values of each continuous variable in different hierarchical forms and at different numerical intervals. The factor and interaction detector modules were selected to explore the dominant drivers of variation in coastal landscape pattern stability (Wang et al., 2021).

The factor detection module is mainly applied to determine the explanatory power of each factor driving the spatial differentiation of landscape pattern stability, expressed by q (Wang et al., 2010), as follows (Wang et al., 2022b):

where h is the identification partition, h = 1, 2, …, L, where L is the number of partitions; $ {N}_{h} $ and N are the numbers of units in layer h and the whole region, respectively; $ {\sigma }_{h}^{2} $ and $ {\sigma }^{2} $ are the variances of the coastal landscape pattern stability values of layer h and the whole region, respectively; and SSW and SST are the sums of intralayer variance and total region variance, respectively. The q value ranges from 0 to 1, where a larger q value indicates stronger explanatory power; thus, q = 1 indicates complete spatial differentiation for independent factors, and q = 0 indicates no spatial variance in the study area and therefore random spatial distribution of coastal landscape pattern stability.

The interaction detection module is mainly applied to identify the interaction between various potential driving factors in determining the spatial differentiation of coastal landscape pattern stability (Song et al., 2020; Wang et al., 2022b). The various types of interaction detection are described in Table 3.

Table  3.  Types of potential factors driving landscape pattern stability

Criteria	Interaction type
q(X1 ∩ X2) < Min(q(X1), q(X2))	nonlinear weakening
Min(q(X1), q(X2)) < q(X1 ∩ X2) < Max(q( X1), q( X2))	single-factor nonlinear attenuation
q(X1 ∩ X2) > Max(q(X1), q(X2))	two-factor interaction enhancement
q(X1 ∩ X2) = q(X1) + q(X2)	independence
q(X1 ∩ X2) > q(X1) + q(X2)	nonlinear enhancement

 | Show Table

DownLoad: CSV

In this study, 41 country districts were used as detection areas. Based on previous studies and data accessibility, the variation of coastal landscape pattern stability in each township from 1985 to 2020 was defined as Y. The eight factors comprised four socioeconomic factors (X₁, change in population density; X₂, change in agricultural output; X₃, change in gross industrial output; and X₄, change in road network density) and four natural environment factors [X₅, digital elevation model (DEM); X₆, slope; X₇, average annual temperature change; and X₈, average annual precipitation change].

4. Results

4.1 Selection of the optimal landscape pattern stability scale for the coastal zone of Guangdong Province

The coastal zone of Guangdong Province has undergone rapid urbanization over the past several decades. Therefore, we explored the scale-dependent characteristics of coastal landscape patterns at landscape and class levels, and then measured landscape fragmentation during urbanization at the optimal scale and landscape level.

4.1.1   Landscape sensitivity effects reflected changes in spatial resolution

The highly complex landscape patterns in coastal zones are closely related to their processes, and they change with scale. These multi-scale pattern changes reflect the hierarchical characteristics of ecosystems. The ability to measure the scale-dependent characteristics of coastal landscape patterns is a prerequisite for the quantitative detection of the mechanisms driving dynamic landscape pattern evolution. Therefore, to explore the scale threshold of landscape patterns at the landscape and class levels, seven landscape metrics and six land use types defined in 2020 were selected and 11 spatial resolutions (or grain sizes) were defined: 30 m × 30 m, 60 m × 60 m, 90 m × 90 m, 120 m × 120 m, 150 m × 150 m, 180 m × 180 m, 210 m × 210 m, 240 m × 240 m, 360 m × 360 m, 480 m × 480 m, and 600 m × 600 m. Due to the basic pixel size of the Landsat image data being 30 m × 30 m, to ensure the continuity and integrity of the ground object information as much as possible, the 11 spatial resolutions (or grain sizes) were all multiples of the basic pixel size.

At the landscape level, PD, TECI, AI and CONTAG metrics decreased slightly at first and then sharply with increasing spatial scale (Figs 3a and b). However, both MPS and MPI metrics increased slightly at first and then significantly with increasing spatial scale (Fig. 3a). In addition, the SHDI metric fluctuated up at first, and then fluctuated notably down, and reached a clear peak of 1.351 at 240 m × 240 m spatial scale (Fig. 3b). In general, seven landscape-level pattern metrics response curves indicated that high sensitivity of the coastal landscape pattern to a changing spatial scale.

Figure  3.  Seven landscape-level pattern metrics response curves to progressively increasing grain size. PD: patch density; MPS: mean patch size; TECI: total edge contrast index, MPI: mean proximity index, AI: aggregation index, CONTAG: contagion; SHDI: Shannon’s diversity index; p < 0.01.

DownLoad: Full-Size Img PowerPoint

In order to further explore the scale sensitivity of coastal landscape pattern at the class level, PD and AI metrics were mainly selected to measure the response curves for six land use types to progressively increasing grain size. As shown in Figs 3c and d, both PD and AI decreased significantly with increasing spatial scale, suggesting high sensitivity of the landscape pattern to a changing spatial scale. Specifically, the PD of farm land deceased substantially, from 1.55 at a scale of 30 m × 30 m to 0.28 at 240 m × 240 m, and decreased slightly with the continuous increase of the spatial scale (Fig. 3c). A similar phenomenon occurred in forest and construction land use types. PD for water body, mangrove, and other types changed slightly at multi-scale levels. In both mangrove and other types, AI decreased from 80.73% and 77.80% at a spatial scale of 30 m × 30 m to 31.92% and 33.99% at 240 m × 240 m, respectively, and decreased slightly with the continuous increase of the spatial scale (Fig. 3d). Similar trends were observed for construction, water body, forest, and farm land use types.

A global autocorrelation model was applied to quantitatively investigate the spatial aggregation of PD and AI at different spatial scales, with significance evaluated using Moran Ⅰ, the Z score, and the p value. As shown in Fig. 3e, as the spatial scale increased from 30 m × 30 m to 600 m × 600 m, Moran Ⅰ decreased slightly, then increased, and finally decreased dramatically for both PD and AI. Specifically, PD and AI showed stronger spatial aggregation at a spatial scale of 240 m × 240 m than at the other spatial scales, with Moran Ⅰ index values of 0.624 and 0.559, respectively (Fig. 3e). The Z scores of PD and AI at a spatial scale of 240 m × 240 m were 2 155.88 and 2 022.41, respectively, whereas p < 0.01 for all spatial scales.

Generally, seven landscape metrics changed considerably in response to changing spatial scales, particularly at a spatial scale of 240 m × 240 m; Moran Ⅰ was also highest at this spatial scale, indicating that 240 m × 240 m is the optimal spatial scale threshold for analyzing coastal zone landscape stability patterns.

4.1.2   Changes in the overall landscape pattern during urbanization

With accelerating urbanization in the coastal zone of Guangdong Province, PD was stable and relatively low from 1985 to 2010, but increased greatly from 2010 to 2015, suggesting a dramatic increase in landscape fragmentation during this period. After 2015, PD decreased slightly (Fig. 4a). By contrast, MPS decreased gradually during urbanization, especially from 2010 to 2015, and then showed a mild rise from 2015 to 2020 (Fig. 4a), which was consistent with a more fragmented landscape after 2010. TECI generally decreased at first and then increased with a turning point around 2000 (Fig. 4a). SHDI increased gradually from 1985 to 2000, reached a clear peak of 1.56 in 2000 (Fig. 4b), and then generally decreased. Thus, as urbanization unfolded, the area of a particular land use type decreased and the landscape gradually became fragmented.

Figure  4.  Spatiotemporal patterns of urbanization in the coastal zone of Guangdong Province as described by seven landscape-level pattern metrics at the optimal analysis scale. PD: patch density; TECI: total edge contrast index; MPS: mean patch size; MPI: mean proximity index; AI: aggregation index; CONTAG: contagion; SHDI: Shannon’s diversity index.

DownLoad: Full-Size Img PowerPoint

CONTAG measures the degree of patch clumping and pertains to both the compositional and the configurational aspects of the landscape pattern. Overall, CONTAG increased slightly, with the greatest landscape fragmentation occurring in 2010 (Fig. 4b). The AI trend showed that the degree of aggregation decreased constantly from 2010, reaching a minimum in 2015 (Fig. 4b). MPI generally remained stable over time, with no significant changes. Thus, urbanization in the coastal zone of Guangdong Province generally resulted in an increasingly rapid increase in PD and large decreases in MPS, AI, and CONTAG, indicating increasing landscape fragmentation and shape complexity.

4.2 Evolution of coastal landscape pattern stability

4.2.1   Spatiotemporal patterns in coastal landscape pattern stability

As the optimal scale threshold for coastal landscape pattern analysis in Guangdong Province was 240 m × 240 m, we used this scale to examine spatiotemporal patterns in coastal landscape pattern stability, calculated using Eq. (1), for Guangdong Province from 1985 to 2020.

Generally, the mean landscape pattern stability of Guangdong Province increased slightly from 1985 to 2005, but decreased significantly from 2005 to 2020 (Fig. 5), with values 12.68% and 15.33% higher in the eastern coastal zone (ECZ) and middle coastal zone (MCZ) than in the western coastal zone (WCZ) of Guangdong Province (Fig. 5i). Spatially, landscape pattern stability was higher in areas with intact patches and clear boundaries, and lower in transition zones, where multiple land types interleaved. Regions of high stability were mainly the large water bodies of the Zhujiang River Estuary and the large forests in the eastern MCZ. With the increase in urban built-up areas, the artificial boundary became increasingly obvious, and the landscape pattern stability of the central urban area was gradually enhanced (Xu et al., 2018). By contrast, regions of low stability were mainly detected in the northern MCZ and most of the WCZ. Over time, landscape pattern stability exhibited significant stage characteristics. From 1985 to 2000, we observed large increases in the mean landscape pattern stability (S) of the ECZ, WCZ, and MCZ, by 26.88%, 24.05%, and 18.22%, respectively (Fig. 5i), followed by slight changes in all three regions from 2000 to 2005. From 2005 to 2020, the S values of the ECZ and MCZ decreased significantly, by 41.68% and 34.80%, respectively, whereas that of the WCZ decreased from 0.79 to 0.51 from 2010 to 2020, representing a difference of 35.79% (Fig. 5i).

Figure  5.  Spatiotemporal patterns of coastal landscape stability in Guangdong Province from 1985 to 2020. WCZ: western coastal zone of Guangdong Province; MCZ: middle coastal zone of Guangdong Province; ECZ: eastern coastal zone of Guangdong Province.

DownLoad: Full-Size Img PowerPoint

4.2.2   Variation in coastal landscape pattern stability

The evolution of spatial and temporal landscape pattern stability in the coastal zone of Guangdong Province differed significantly over the study period (Fig. 6). From 1985 to 1990, landscape pattern stability was maintained or increased (Fig. 6a), with a large increase in the area of stability observed in the northern regions of the MCZ (Panyu and Zhongshan) and ECZ (Chaozhou) and the western WCZ (Zhanjiang), and a large decrease at the junction of forest, grass, farmland, and city land use types in the western MCZ and WCZ. The evolution of landscape pattern stability from 1990 to 1995 was dominated by a slight decrease in stability (Figs 6b and h), but with a significant decrease in the area of stability mainly occurring in the urban fringe area (northern MCZ and southern WCZ). From 1995 to 2000, an increase in the area of stability was mainly observed in the eastern (Huizhou), western WCZ (Zhanjiang), and central ECZ (Shanwei), and a decrease in the area of stability in the western MCZ (Zhongshan and Jiangmen) (Fig. 6c). From 2000 to 2005, an increase in the area of stability was mainly observed in the northern MCZ (Dongguan and Zhongshan), eastern ECZ (Chaozhou), and western WCZ (Zhanjiang), whereas stability decreased across a wide area covering the eastern MCZ, western ECZ, and eastern WCZ (Fig. 6d).

Figure  6.  Changes in coastal landscape pattern stability in Guangdong Province from 1985 to 2020. WCZ: western coastal zone of Guangdong Province; MCZ: middle coastal zone of Guangdong Province; ECZ: eastern coastal zone of Guangdong Province. Landscape stability decreased: the annual mean stability metric belowed –0.05; landscape stability increased: the annual mean stability metric begonded 0.05; landscape relative stability: the annual mean stability metric ranged from –0.05 to 0.05.

DownLoad: Full-Size Img PowerPoint

Generally, the evolution of landscape pattern stability in the coastal zone in Guangdong Province was characterized by a decrease after 2005, including a 41.24% decrease in area from 2010 to 2015 (Fig. 6h). However, there was a significant increase in the area of stability in the eastern MCZ and ECZ from 2005 to 2010, the western WCZ from 2010 to 2015, and the northern ECZ from 2015 to 2020. Due to continuous expansion of urban construction and the fragmentation of forest, grass, and cultivated land, the degree of land patch aggregation decreased, the degree of spreading increased, the boundary became more complex, and landscape pattern stability gradually decreased.

4.3 Mechanisms driving coastal landscape pattern stability

According to our results and a preliminary study of the coastal zone in the Zhujiang River Estuary (Chen et al., 2023), X₁ (change in population density), X₅ (DEM), and X₆ (slope) were optimally divided into 14, 13, and 11 categories, respectively, using the geometric interval method. X₂ (change in agricultural output), X₃ (change in gross industrial output), X₄ (change in road network density), X₇ (average annual temperature change), and X₈ (average annual precipitation change) were divided into 9–15 categories using the natural interval method.

The relative importance of each factor was detected using the geodetector method based on optimal parameters and expressed as a q value (Wang and Xu, 2017; Wu et al., 2023). Among the eight factors driving landscape pattern stability evolution, X₂, X₃, X₈, and X₇ were ranked as the top four factors, with average q values of 0.203, 0.191, 0.107, and 0.087, respectively (Fig. 7). A comparison of the q values showed that each of these four factors changed uniquely during the seven periods. From 1985 to 1995, q values were low, and the change in the agricultural output value played a critical role in coastal landscape pattern stability, with q values of 0.096 during 1985–1990 and 0.158 during 1990–1995 (Figs 7a and b). From 1995 to 2010, the explanatory power of X₂ increased gradually, and from 2010 to 2015, the explanatory power of this factor and the change in X₃ increased, with q values of 0.404 and 0.417, respectively (Fig. 7f), whereas the explanatory power of the changes in X₇ and X₈ decreased. From 2015 to 2020, the explanatory power of almost all factors strongly decreased, with q < 0.1 for most factors (Fig. 7g). Generally, the regularity of the explanatory power of the various driving factors on the evolution of coastal landscape pattern stability was poor, indicating unique and complex short- and long-term links between variations in natural and anthropogenic factors and landscape pattern stability in Guangdong Province, and most likely other major global coastal zone systems.

Figure  7.  q values of influential factors (X₁−X₈) from 1985 to 2020. X₁: Change in population density; X₂: change in the agricultural output value; X₃: change in the gross industrial output value; X₄: change in the road network density; X₅: DEM, digital elevation model; X₆: slope; X₇: average annual temperature change; X₈: average annual precipitation change. p < 0.01.

DownLoad: Full-Size Img PowerPoint

Essentially, landscape pattern stability evolution in Guangdong Province reflected the rapid urbanization of China’s coastal zones, subsuming interactions between natural drivers and anthropogenic interferences. The geodetector method was applied to quantitatively identify the mutual explanatory power of paired factors with respect to the change in coastal landscape pattern stability; the results are presented in Fig. 8. Among the evaluated factors, the explanatory power of paired factors was greater than that of single factors, with significantly more non-linear enhancements (152 groups) than two-factor enhancements (44 groups). Therefore, the coastal landscape pattern stability of Guangdong Province was likely the result of a combination of factors. From 1985 to 2000, the non-linear reinforcement effect among factors X₄, X₅, and X₆ was weak, whereas that among X₂, X₃, and other factors was relatively strong. For example, the q values of X₂ ∩ X₇ in 1985–1990, X₂ ∩ X₇ in 1990–1995, and X₃ ∩ X₇ in 1995–2000 were 0.334, 0.378, and 0.360, respectively (Figs 8a and c), suggesting that combined changes in agricultural and industrial development, temperature, and precipitation drove the variation in coastal landscape pattern stability during those periods. However, the interactions among X₂, X₃, and other factors resulted in two-factor enhancements and were slightly strengthened from 2000 to 2010, with the q values of X₃ ∩ X₈ in 2000–2005 and X₂ ∩ X₇ in 2005–2010 reaching 0.398 and 0.381, respectively (Figs 8d and e). The q values of X₇ ∩ X₈ during the same periods were 0.417 and 0.405 (Figs 8d and e), respectively, indicating that climate change was the key factor driving variation in landscape pattern stability. From 2010 to 2015, most of the interactions among X₂, X₃, and other factors changed from two-factor enhancements to non-linear reinforcement and were further strengthened, with the q values of X₃ ∩ X₈ and X₃ ∩ X₇ reaching 0.604 and 0.601 (Fig. 8f), respectively. The q values of X₁ ∩ X₂, X₁ ∩ X₃, X₂ ∩ X₄, and X₃ ∩ X₄ were all higher than 0.5. From 2015 to 2020, the q values of any two factors strongly declined and the explanatory power was weaker than that during any previous period.

Figure  8.  Interaction detection results of potential influencing factors (X₁−X₈). X₁: change in population density; X₂: change in the agricultural output value; X₃: change in the gross industrial output value; X₄: change in the road network density; X₅: DEM, digital elevation model; X₆: slope; X₇: average annual temperature change; X₈: average annual precipitation change.

DownLoad: Full-Size Img PowerPoint

In conclusion, anthropogenic factors had stronger explanatory power than natural factors in explaining the variation in coastal landscape pattern stability in Guangdong Province. However, our geodetector findings suggested that anthropogenic and natural factors were mutually dependent and that the influence of multifactor interactions (i.e., changes in agricultural and industrial output, average annual temperature, and average annual precipitation) on the spatial differentiation of the variation in landscape pattern stability consisted of double-factor and nonlinear reinforcements rather than a simple linear superposition.

5. Discussion

5.1 Complex evolution of coastal landscape pattern stability under rapid urbanization

Urbanization is the most extreme form of land transformation (Wu et al., 2011), as it results in the direct loss and fragmentation of natural habitats and profoundly affects biodiversity and ecosystem processes (Grimm et al., 2008). In this study, we explored the complex evolution of the coastal landscape in response to rapid urbanization by introducing a landscape metric analysis of coastal landscape pattern stability in Guangdong Province. From 1985 to 2000, coastal landscape pattern stability increased gradually, especially in areas of rapid urbanization (i.e., the MCZ) (Figs 5a–d). This trend was mainly due to the change in urban development from a point-like expansion to planar development (Zhang et al., 2022), with the pattern of newly built-up areas around the original building land maintaining relatively intact landscape patches and thus improving landscape pattern stability (Xu et al., 2018). The landscape pattern stability of the Guangdong coastal zone reached a maximum from 2000 to 2005 (Fig. 5i), with the most dramatic land cover changes occurring during this period, mainly in the form of transitions from cropland to artificial built-up areas (Liu and Ma, 2017; Lin et al., 2020). However, after 2005, the landscape pattern stability of the Guangdong coastal zone declined dramatically, mainly due to disorderly township expansion, which eroded both forests and cultivated land, resulting in increased fragmentation of land patches (Fig. 3a) and increasingly complex boundaries of various land use types (Xu et al., 2018). Serra et al. (2014) suggested both positive and negative impacts of urbanization on landscape pattern stability. For example, positive effects are observed when urban intensification is conducive to reducing anthropogenic intervention in surrounding forests, water systems, and cultivated land (Serra et al., 2014; Xu et al., 2018), thus maintaining relatively intact regional landscape patches and enhancing landscape pattern stability. Negative effects mainly occur in the disorderly expansion of towns and villages and through coastalization (Serra et al., 2014), causing marginalization, landscape fragmentation, and habitat loss in coastal hinterlands (Rico-Amoros et al., 2009). Cantos et al. (2010) found that urban expansion in coastal regions in the form of consolidated towns and new residential areas increased the flood risk. The short- and long-term links between the speed, scale, and type of urban expansion and landscape pattern stability in major coastal zones worldwide are both unique and complex.

5.2 Effect of scale change on coastal landscape pattern stability

Landscape pattern stability is rooted in landscape spatial heterogeneity, which in turn stems from variations in spatial dependence (Wu, 2004). Because landscape metrics vary with spatial magnitude and granularity, their use requires clear indication of the scale of analysis (Wu, 2007). Previous studies have shown that remote sensing data with different spatial resolutions may lead to differences in spatial analysis results for the same area (Wickham et al., 1997; Gong et al., 2006). The plasticity of area units has also had a significant impact in studies involving spatial data in ecology, because ecosystems and their patterns exist at multiple scales or have a hierarchical structure (Wu et al., 2011), and patterns and processes can vary at different scales (Rao et al., 2021; Wu et al., 2023). Therefore, this study mainly explored the determination of optimal analysis scale for landscape pattern stability in Section 4.1, which established an important foundation for the subsequent quantitative analysis of the evolution of landscape pattern stability patterns and their driving mechanisms.

Coastal zones have highly complex landscape heterogeneity, the degree of which varies with scale (Chen and Dong, 2003; Gong et al., 2006). The stability of the coastal landscape is based on spatiotemporal variations, with different landscape characteristics manifesting at different scales (Luo et al., 2004; Xie et al., 2005). In their study of Haikou, China, Liu et al. (2018) suggested that the appropriate landscape grain size for landscape pattern stability analysis based on land use was 70 m × 70 m, whereas another study evaluating landscape pattern stability in the Beijing–Tianjin–Hebei region applied a spatial grain size of 20 km × 20 km (Xu et al., 2018). Coastal landscape pattern stability is more complex at the provincial scale than at the city or urban agglomeration scale. As Guangdong Province has a long, highly varied coastline compared to other Chinese coastlines (Zhao, 2017), the first step in the analysis of its coastal landscape pattern stability was to identify the optimal grain size. As shown in Fig. 3e, we found that Moran Ⅰ index and the Z score were much higher at a spatial scale of 240 m × 240 m than at other spatial scales, indicating that this spatial scale had the largest impact on coastal landscape pattern stability, which was consistent with the findings of other studies (Xu et al., 2020, 2021; Hu et al., 2023). Because the computation of a landscape metric at an improper resolution can lead to ecological fallacies in correlation and regression analyses (Wu, 2004), it is critically important to identify an appropriate level of detail for the ecological processes of interest. Moreover, the MPS metric declined from 132.90 hm² in 2005 to 84.21 hm² in 2015, a decrease of 36.65% (Fig. 4a). On the contrary, the average landscape pattern stability index of the ECZ, MCZ, and WCZ decreased by 35.95%, 29.18%, and 25.17% (Fig. 5i), respectively. This indicated that the patch size played a crucial role in maintaining landscape pattern stability. As several studies have shown that aggregation methods can affect the scaling relations of landscape metrics (Turner et al., 2001; Wu et al., 2002; Wu, 2004), it is of interest to compare the effects of such methods on the characteristics of coastal landscape pattern stability.

In addition, the study outcomes allowed specific guidelines for efficient utilization and integrated management of coastal zone resources to be formulated. (1) The optimal spatial scale of 240 m × 240 m for coastal landscape pattern stability in the present study should be applied to guide the delineation of coastal wetland reserve areas and community parks. (2) The urban-rural transition zones with significant changes in landscape pattern stability identified in this study should draw more attention of multiple government departments and should be developed and managed according to local conditions in the future.

5.3 Complexity and uncertainty of driving force detection

Coastal zones are ecological transition zones where land and marine ecosystems interact to form a rich geographical unit and distinctive landscape (Yuan and Chang, 2021). However, because of their complexity and marginality, coastal zones are also ecologically fragile areas that present significant challenges (Zheng, 2023). Landscape patterns in coastal regions reflect the interactions of natural and anthropogenic interferences that develop over long periods of time (Yeh and Huang, 2009). In this study, eight factors potentially influencing coastal landscape pattern stability were quantitatively examined using optimal parameters according to geodetector analysis. As shown in Fig. 7, socioeconomic factors (i.e., X₂ and X₃) had stronger exploratory power than natural factors in changing the coastal landscape pattern stability of Guangdong Province over the past several decades. A combination of q statistics showed that the explanatory power of paired factors was greater than that of single driving factors. For example, the q values of X₃ (change in gross industrial output) ∩ X₈ (change in average annual precipitation) and X₃ ∩ X₇ (change in average annual temperature) were 0.604 and 0.601, respectively (Fig. 8f). These dominant drivers of landscape pattern stability variations are consistent with those of previous studies (Jiang et al., 2021; Li et al., 2022; Wang et al., 2023). However, the additional impacts of mariculture, tidal flat reclamation, and storm surges on coastal landscape pattern stability (Liu and Ma, 2017; Tian et al., 2021) could lead to complexities and uncertainties in the analysis of coastal landscape pattern stability.

Governmental policy ultimately determines coastal urbanization and is the leading force in coastal landscape evolution. In 2001, the “Outline of the Tenth Five-Year Plan for National Economic and Social Development’ of Guangdong Province” stated that urbanization should be accelerated and that the urbanization level of the Zhujiang River Delta region should reach about 60% by 2005 (People’s Government of Guangdong Province, 2001). In the 2013 “Twelfth Five-Year Plan for Urbanization Development of Guangdong Province”, an annual increase in the urbanization rate of the province by 0.8% was targeted, with a 70% increase to be reached by 2015 (People’s Government of Guangdong Province, 2011). These governmental policies accelerated coastal zone urbanization, with the resulting expansion of built-up areas significantly eroding green space, water bodies, and farmland, leading to a sharp decline in coastal landscape pattern stability in Guangdong Province after 2005 (Fig. 5). Despite the key role of governmental policies in the development of coastal zones, quantifying their impact on the evolution of coastal landscape pattern stability remains difficult, such that it was not considered in this study and is a potential limitation.

5.4 Future studies

Although many studies have shown that rapid urbanization leads to a decline in coastal landscape pattern stability (Tian et al., 2021; Li et al., 2022), this study demonstrated that the relationship between urbanization and coastal landscape pattern stability is not a simple linear negative correlation, but rather a nonlinear relationship that is closely related to the expansion rate, stage, and form of coastal urbanization. Studies on the quantitative relationship between urbanization and coastal landscape pattern stability will improve the governance and sustainable development of coastal zones. Further discussion and in-depth study of quantifiable policy factors as well as improved approaches toward analyzing the structural system of driving factors are needed to understand the development direction of coastal zones and trends in the evolution of landscape pattern stability under the guidance of governmental policies.

6. Conclusions

In this study, landscape metrics were calculated for the coastal zone of Guangdong Province at 11 different scales. From the dual perspective of spatial heterogeneity and spatial autocorrelation, the optimal scale for coastal landscape pattern stability in Guangdong Province was selected and then used in an overlay analysis to explore the spatiotemporal evolution of coastal landscape pattern stability. Geodetector analysis at an optimal scale based on optimal parameters revealed the forces driving coastal landscape pattern stability. In Guangdong Province, the MCZ was the major distribution area, with high and low values of landscape pattern stability. In recent decades, the annual mean stability of the coastal landscape in Guangdong Province increased slightly, and then decreased substantially with rapid urbanization, with a transition period occurring around 2005. The most significant changes in landscape pattern stability were mainly observed in the transition zone of rural urban expansion.

The response characteristics of the coastal landscape pattern suggested that a spatial scale of 240 m × 240 m was optimal landscape pattern stability analysis in Guangdong Province. Based on this scale, the optimal discretization method and classification categories of each driving factor were obtained using geodetector analysis, a practical statistical method that has provided reliable scientific data support for other coastal zones worldwide subjected to similar stresses.

The mechanisms driving coastal landscape pattern stability differed markedly over time. The use of q statistics revealed that the explanatory power of paired factors was greater than that of single driving factors, with significantly greater non-linear reinforcement than double-factor reinforcement. Changes in agricultural and industrial output, average annual temperature, and average annual precipitation were the primary determinants of coastal landscape pattern stability changes in Guangdong Province, with the q values of X₃ ∩ X₈ and X₃ ∩ X₇ reaching 0.604 and 0.601, respectively, from 2010 to 2015. In the future, additional efforts will be needed to guide the sustainable high-quality development of human-land relationships and the regional economy in the coastal zone of Guangdong Province in accordance with local conditions.