Introduction

The level of urbanization in China has accelerated rapidly since the twenty-first century. The urbanization rate of permanent residents has increased from 36.2% in 2000 to 64.72% in 2021. Rapid urbanization has brought enormous economic and social benefits, but has also exacerbated the ecological crises, such as biodiversity loss1. In the face of ecological problems caused by rapid urbanization, conferences such as the Future Earth Plan and the 2030 Agenda for Sustainable Development emphasize the coordination of urbanization with ecological civilization, along with the promotion of the sustainable development of biodiversity2. The United Nations has identified 17 Sustainable Development Goals (SDGs)3, that provide specific criteria for assessing environmental, economic and social sustainable development. The quantification of SDGs has become an important basis for countries and regions to comprehensively evaluate the progress of their sustainable development.

With the acceleration of urbanization, human beings are causing biodiversity destruction at an unsustainable rate4,5,6. Biodiversity is the sum of ecological complex formed by organisms and their environment, as well as various ecological processes related to it7. It usually includes genetic, species, ecosystem, and landscape diversity. Currently, the scope of biodiversity evaluation in China is at the county, city, provincial, and national scales, among which the spatial distribution pattern at the county scale is a prerequisite for implementing regional biodiversity conservation measures8,9,10. Traditional biodiversity studies were dominated by expensive and laborious field investigation methods, including distribution and quantitative surveys of specific species11,12. With the continuous development of remote-sensing technology, more high-resolution data sources have been used to construct and monitor indicators of regional biodiversity13,14. The biodiversity monitoring methods can be roughly divided into two categories—the first uses a single remote-sensing feature to predict species richness and landscape diversity15,16,17, and the second indirectly evaluates biodiversity by establishing a comprehensive indicator system, which usually includes multilevel indicators, such as species, ecosystems, and landscapes. For example, Liu et al.10,18, Peng et al.19, Yang et al.13, Li et al.20, and Deng et al.21 selected multiple environmental indicators, such as habitat quality, species distribution, ecosystem status, and landscape structure, to build a comprehensive evaluation system based on remote sensing technology for the dynamic monitoring of county biodiversity. In conclusion, conducting a comprehensive evaluation of biodiversity is feasible at small and medium scales based on remote-sensing technology and ecological models in areas where fieldwork is lacking. In addition, the United Nations has established specific targets for biodiversity evaluation. SDG 15.1.2 measures the “proportion of important sites for terrestrial and freshwater biodiversity that are covered by protected areas, by ecosystem type,” focusing on the identification and extraction of key biodiversity areas (KBAs). However, the SDG framework has several limitations in terms of scale, data, and calculation methods. The application of SDG 15.1.2, therefore, needs to be localized according to the situation in different study areas22,23,24.

Urbanization is closely related to the ecosystem. Existing urbanization research generally adopts single and comprehensive index methods to quantify the urbanization levels, where the comprehensive method can reflect the overall urbanization levels more comprehensively25. Human activities have exerted a profound impact on the changes in ecosystem and urban spaces with the development of the social economy. Human activity intensity of the land surface (HAILS) represents the degree of human utilization, transformation, and exploitation of the land surface in a certain region. From the perspective of land use, HAILS is calculated based on the change in natural land properties caused by social and economic activities of humans26 and is widely used in the Loess Plateau and Yellow River Delta27,28. Currently, the study of the interaction between urbanization and ecosystem is of significant interest among scholars29,30,31,32. Ecosystems with high biodiversity are better at resisting environmental changes. However, there are relatively few studies on the coupling coordination between biodiversity and urbanization at the county scale are also limited. According to the system coupling theory, two or more subsystems can produce a coupling phenomenon through interaction. The coupling degree reflects the coordination degree of the interaction between subsystems33. The coupling coordination degree (CCD) model has obvious advantages in analyzing the relationship between complex systems because it can better reflect the structure and function of complex systems34. The relative development degree29 is typically used to examine the degree of adaptability and synchronization of two systems in their interactive developmental processes. The introduction of grey relational degree model35 to identify the key influencing factors of the composite system plays a crucial and effective role in improving the ability of the overall coordination level. Summarizing, biodiversity and urbanization are multidimensional systems. The interaction between them constitutes a more complex system, in which various elements constitute several complex connections that affect its development. Therefore, the study of coupling mechanism based on multi-source data between biodiversity and urbanization is necessary.

Although the scholars have conducted various studies on biodiversity, urbanization and their relationship, some shortcomings are still present. In terms of the evaluation indicator, previous biodiversity studies mainly consider single species indicator and rarely involve ecosystem and landscape conditions. The urbanization indicator, on the other hand, takes less account of the impact of human activities. In terms of the study scale, most studies focus on major strategic areas in China. Very few studies are available on the lower Yellow River at the county scale. In terms of the data type, the construction of indicator system is primarily based on statistical data and lacks multiple source data types. In terms of the study content, most scholars have only studied a certain theory of biodiversity or urbanization, and seldom attempt to visualize them from the perspective of remote sensing. The interaction mechanism and its influencing factors are usually not discussed. In terms of study significance, previous studies have lacked effective suggestions considering the sustainable coordinating development of multiple influencing factors.

The areas along the Yellow River of Shandong province (AYRSP) play an irreplaceable position and role in major national strategies such as the Belt and Road Initiative and the ecological protection and high-quality development in the Yellow River basin. Additionally, the AYRSP boasts rich terrestrial and marine ecosystems, nurturing abundant biodiversity resources, which provide valuable samples for biodiversity research at home and abroad. Impacted by resource constraints, environmental changes and rapid urbanization, the vulnerability of ecosystems is increasing, and biodiversity conservation is under threat36. Therefore, the demand of coordinated development between urban expansion and biodiversity protection in the AYRSP has become increasingly urgent. In 2021, the outline of Ecological Protection and High-quality Development in the Yellow River Basin clearly emphasized the ecological protection and management of downstream wetlands, so as to maintain the biodiversity of the lower Yellow River. In 2022, Shandong issued an outline of Ecological Protection and High-quality Development in AYRSP. This outline emphasized that the Shandong section of the Yellow River is the main area for ecological protection in the lower Yellow River region. Currently, there are many studies on ecosystems and urbanization in the AYRSP both domestically and internationally2,29,32, but there are few studies on the complex interaction between biodiversity and urbanization from the empirical aspect. Quantitatively studying the interactive coupling relationship between urbanization and biodiversity in the AYRSP has strategic significance for promoting the innovation of sustainable development mode and achieving the coordination between economic development and ecological protection within the Yellow River basin and even globally. It is a crucial issue in coordinating human-land relationships and promoting sustainable development. Taken the AYRSP as the study area, this study evaluated the spatio-temporal heterogeneity and coupling relationship between biodiversity and urbanization under the land-use change based on multi-source data and multiple models. In terms of system coupling, the study provides reference suggestions for improving regional coordination level which helps to promote sustainable development of biodiversity protection, ecosystem stability, and new urbanization in the lower Yellow River region. Meanwhile, it also provides reference for the development of other provinces along the Yellow River.

Study area and data sources

Study area

The AYRSP is in the lower reaches of the Yellow River. The Yellow River Basin is an important ecological defense area in China and is strategically important for economic development and ecological security. As the only province located east of the nine provinces along the Yellow River, Shandong Province plays an important role in its protection strategies. The Yellow River flows from Dongming County in Heze City in the Shandong Province into the Bohai Sea in Kenli County in Dongying City. It flows 628 km from west to east through nine cities and 25 counties in Shandong (Fig. 1). The terrain in most areas is relatively gentle, and the Mount Tai is primarily located in the middle and south of the AYRSP. In recent years, with rapid urbanization, human activities have increasingly intervened this region. Consequently, biodiversity is threatened, which has a serious impact on ecological security. Moreover, biodiversity protection is an important aspect of ecological protection and high-quality development in the Yellow River. Shandong issued the Biodiversity Protection Regulations to strengthen biodiversity protection and revised and implemented a provincial biodiversity protection strategy and action plan in the country.

Fig. 1
figure 1

Location of the research area. Generated in the ArcGIS 10.6 software (www.esri.com).

Full size image

Data sources

According to the standard of Monitoring Contents and Indicators of Geographical National Conditions (GQJC 03-2020), the actual situation of Shandong Province, and the research demand, the geographic national conditions monitoring land use data37 (resolution: 2.5 m) in 2015, 2017, 2019, and 2021 were reclassified into 15 categories (Fig. 2). It was used to calculate the habitat quality index (HQI), percentage of habitat area (\({S}_{p}\)), landscape indicators, and land urbanization indicators.

Fig. 2
figure 2

Geographic national conditions monitoring data of the AYRSP in 2021. Generated in the ArcGIS 10.6 software (www.esri.com).

Full size image

HQI was calculated by using the Habitat Quality module of Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model. Based on the Google Earth Engine platform (https://code.earthengine.google.com/), this study selected Landsat 8 images (resolution: 30 m) of vegetation growing seasons (May–October) in 2015, 2017, 2019, and 2021 after cloud removal, format conversion, data splicing, cropping, and calculation to obtain the normalized difference vegetation index (NDVI) and wetness (WET). Shannon diversity index (SHDI) and landscape shape index (LSI) were calculated by Fragstats 4. The projected coordinate system of all the above remote-sensing data was set to CGCS2000_3_Degree_GK_CM_117E.

The statistical data in society and economy were mainly obtained from statistical yearbooks and statistical bulletins of various cities in Shandong province. In order to ensure the reliability of the data, this study supplemented the small amount of missing data based on the data of adjacent years.

Materials and methods

Framework

This study was divided into three sections (Fig. 3). In the first part, this study selected six indicators from species, ecosystem, and landscape diversity to construct a remote sensing system based on the InVEST model. Then, the level of biodiversity sustainability was quantified by the proportion of KBAS. In the second part, five indicators were selected from three dimensions to build comprehensive urbanization system based on panel data. In the third part, the coupling coordination relationship between biodiversity and urbanization was quantified based on the CCD model. The lead-lag type was quantified between the two in different developmental processes based on the relative development degree model. Furthermore, the key influencing factors in the coupling coordination system were revealed based on the grey relation model.

Fig. 3
figure 3

Analysis framework of this study.

Full size image

Evaluation of biodiversity

Species diversity

HQI

Habitat quality determines the ability of the natural environment to provide for the suitable survival of individuals and sustainable development of the populations. It is positively correlated with biodiversity abundance, and is often used as an important index to assess biodiversity38,39. The habitat quality module of InVEST model can characterize the response degree of different habitats to threat factors and the interactions between threat factors40. HQI has some reference value in areas where species distribution information is poor or where there is a mixture of habitat types coexisting. To a certain extent, HQI is a representation of biodiversity. The calculation process is as follows9,10,40:

$${Q}_{xj}={H}_{j}[1-(\frac{{D}_{xj}^{z}}{{D}_{xj}^{z}+{k}^{z}})]$$
(1)

where \({Q}_{xj}\) is the habitat quality of land type \(j\) in grid cell \(x\), \({H}_{j}\) is the habitat suitability of land type \(j\), \(k\) is half-saturation and constant, and \(z\) is the model default parameter which can be set to 2.5. By referring to the recommended values of the model, current biodiversity situation of the lower Yellow River and related study results41,42, the value of the correlation factor of the threat source are determined. Threat factors of related data in the AYRSP are shown as Table 1, as well as the data of the habitat suitability and the sensitivity of threat factors are shown as Table 2.

$${D}_{xj}=\sum_{r=1}^{R}\sum_{y=1}^{{Y}_{r}}\left(\frac{{w}_{r}}{\sum_{r=1}^{R}{w}_{r}}\right){r}_{y}{i}_{rxy}{\beta }_{x}{S}_{jr}$$
(2)

where \({D}_{xj}\) denotes the habitat degradation degree of grid \(x\) in landscape \(j\), \(R\) is the number of threat factors, \({w}_{r}\) is the weight of threat factor \(r\), \({Y}_{r}\) is the total number of grid cells for threat factor \(r\), \({r}_{y}\) is the intensity of threat factor,\({\beta }_{x}\) represents the anti-interference level of the habitat,\({S}_{jr}\) is the magnitude of sensitivity of the landscape \(j\) to threat factors, \({i}_{rxy}\) means the threat level of threat factor to habitat.

Table 1 Weights for threat factors of the study area.
Full size table
Table 2 Habitat suitability and its relative sensitivity to different threat factors.
Full size table
NDVI

NDVI is a vegetation index that not only reflects the photosynthetic activity of vegetation, but also predicts the richness and diversity of species43. Some scholars believe that there is a significant correlation between vegetation indices such as NDVI and plant diversity44. The formulation is as follows:

$$NDVI=({\rho }_{NIR}-{\rho }_{red})/({\rho }_{NIR}+{\rho }_{red})$$
(3)

where \({\rho }_{NIR}, {\rho }_{red}\) represent the reflectance of near-infrared and red bands of Landsat 8 images, respectively.

Ecosystem diversity

WET

WET obtained through the tasseled-cap transformations has been widely used in ecosystem monitoring. WET reflects the moisture in vegetation, water, and soil, which is the basic condition for the reproduction and growth of organisms, and also reflects the richness of biodiversity in an ecosystem to a certain extent45. The larger the WET is, the more suitable it is for the survival and reproduction of species. It is as shown:

$$WET=0.1511{b}_{2}+0.1973{b}_{3}+0.3283{b}_{4}+0.3407{b}_{5}-0.7117{b}_{6}-0.4459{b}_{7}$$
(4)

where \({b}_{2}\)-\({b}_{7}\) represent the reflectance of bands 2–7 of Landsat 8 /OLI data, respectively.

\({S}_{p}\)

Since the nineteenth century, scientists have recognized the relationship between the number of species and the habitat area. Later, Fu et al.46 pointed out that the habitat patch area not only affects the number of species, but also affects the productivity level and the distribution of energy and nutrients in the ecosystem. Theoretically, there is a positive correlation between habitat area and biological species, and species richness is also correlated with ecosystem type. Therefore, \({S}_{p}\) can better characterize the structure and type of ecosystem. In addition, some studies13,18,47 have used \({S}_{p}\) as a measure of ecosystem diversity to evaluate biodiversity. \({S}_{p}\) is calculated as follows:

$${S}_{p}=\frac{{A}_{e}}{{A}_{l}}$$
(5)

where \({S}_{p}\) denotes the percentage of habitat area; \({A}_{e}\) is the habitat area including forest, grassland, garden, and water body; \({A}_{l}\) is the total ecosystem area.

Landscape diversity

SHDI

SHDI focuses on the diversity of patch types, and can accurately capture the uneven distribution of patch types in the landscape. The higher the degree of landscape fragmentation, the higher the SHDI, and the richer the regional biodiversity18,20,21. SHDI is calculated as follows:

$$SHDI=-\sum_{i=1}^{m}{P}_{i}\text{ln}({P}_{i})$$
(6)

where \(m\) is the number of landscape types,\({P}_{i}\) is the proportion of area occupied by landscape type \(i\).

LSI

LSI focuses on the complexity of patch shape in the landscape, and can comprehensively reflect the shape heterogeneity and connectivity of the landscape. Some studies have found that LSI is significantly positively correlated with regional biodiversity48. The larger the LSI, the higher the complexity of the overall distribution of each landscape, and the richer the biodiversity. LSI is calculated by the following formula:

$$LSI=\frac{0.25E}{\sqrt{A}}$$
(7)

where \(E\) is the total length of all patch boundaries in the landscape, A is the total area of the landscape.

Construction of biodiversity indicator

The above indicators are positive, so the standardized calculation is as follows:

$$X_{ij}^{\prime } = \frac{{X_{ij} - min \left\{ {X_{ij} } \right\}}}{{max \left\{ {X_{ij} } \right\} - min \left\{ {X_{ij} } \right\}}}$$
(8)

where \({X}_{ij}^{\prime}\) is the value of the evaluation indicator after normalization, \({X}_{ij}\) is the value of the evaluation indicator before normalization, \(min\left\{{X}_{ij}\right\}\) is the minimum value in the evaluation indicator before normalization, \(max\left\{{X}_{ij}\right\}\) is the maximum value in the evaluation indicator before normalization.

The analytic hierarchy process (AHP)49 is used to determine the weights. This study made reference to relevant studies10,18,47,50 and fully took into account the influencing factors of biodiversity in the AYRSP to determine the relative importance of indicators. Then, the judgment matrix is constructed, and the eigenvector corresponding to the maximum eigenvalue of the matrix is calculated and normalized to obtain the weight value of each indicator (Table 3). At last, the consistency test is conducted by:

$$CI=\frac{LA-N}{N-1}$$
(9)
$$CR=\frac{CI}{RI}$$
(10)

where \(CI\) is the consistency index, \(LA\) is the maximum eigenvalue of the judgment matrix, \(N\) is the order of the judgment matrix, N = 6; \(CR\) represents the proportion of random consistency, \(RI\) is the average random consistency index. After calculation, CR = 0.06 < 0.1. The judgment matrix passed the consistency test on the whole, so the weight value is reasonable and reliable.

Table 3 Composition and weight of each biodiversity indicator.
Full size table

Finally, the biodiversity indicator is calculated as follows:

$$BI=HQI\times {\beta }_{1}+NDVI\times {\beta }_{2}+{S}_{p}\times {\beta }_{3}+WET\times {\beta }_{4}+SHDI\times {\beta }_{5}+LSI\times {\beta }_{6}$$
(11)

where \({\beta }_{i}\) is the weight coefficient corresponding to each index.

Evaluation of sustainable development level of biodiversity

The quantitative assessment of SDG 15.1.2 mainly focuses on the identification and extraction of KBAs to provide support for biodiversity conservation and sustainable development. The 25 counties on the mainstream of the Yellow River in the Shandong province are the core areas in the protection of Yellow River basin. Accordingly, the AYRSP was regarded as a key policy protection area. We referred to the method of Liu et al.10,18 to evaluate the sustainable development level of biodiversity by taking the proportion of KBAs to the land area. SDG 15.1.2 (\({X}_{s}\)) is calculated as follows:

$${X}_{s}=\frac{{A}_{i}}{A}$$
(12)

where \({A}_{i}\) is the area of high-level biodiversity area (BI ≥ 0.6), which is considered as KBAs; A is the area of the core protection areas. In order to facilitate the comparison of the results, the Sustainable Development Solutions Network (SDSN) delimited the classification threshold for each index value (Table 4).

Table 4 The classification threshold for SDG 15.1.2
Full size table

Evaluation of urbanization level

HAILS

Selecting HAILS to represent the development of land urbanization can reflect the comprehensive effects of human, social, and economic activities on land resources. HAILS is calculated based on the construction land equivalent sum of different land use types in the region26.

$${HAILS}_{i}=\frac{{S}_{CLE-i}}{{S}_{i}}\times 100\text{\%}$$
(13)
$${S}_{CLE-i}=\sum_{j=1}^{m}({SL}_{ij}\times {CI}_{ij})$$
(14)

where \({HAILS}_{i}\) is the human activity intensity of the land surface in area \(i\) (%), \({S}_{CLE-i}\) is the equivalent area of construction land in region \(i\), \({S}_{i}\) is the land area of the region \(i\), \({SL}_{ij}\) is the area of the land type \(j\) in the region \(i\), \({CI}_{ij}\) is the construction land coefficient discounted for the land type \(j\) in the region \(i\), \(m\) is the number of land types in the region \(i\).

Construction of comprehensive urbanization indicator

This study selected the per capita GDP, shares of secondary and tertiary industries, urbanization rate of permanent residents, proportion of built-up areas in total areas, and HAILS from three dimensions to construct an urbanization indicator system using entropy method (Table 5). Additionally, it is necessary to standardize the indicators before constructing the indicator system. The calculation process of urbanization indicators is as follows29:

$${P}_{ikj}={y}_{ikj}/{\sum }_{i=1}^{m}{\sum }_{k=1}^{n}{y}_{ikj}$$
(15)
$${e}_{j}=-\frac{1}{\text{ln}(m\times n)}{\sum }_{i=1}^{m}{\sum }_{k=1}^{n}{P}_{ikj}ln{P}_{ikj}$$
(16)
$${W}_{j}=\frac{1-{e}_{j}}{n-{\sum }_{j=1}^{z}{e}_{j}}$$
(17)
$$UI=\sum_{j=1}^{z}{W}_{j}{Y}_{ij}$$
(18)

where \({P}_{ikj}\) is the proportion of the indicator \(j\) of the county \(k\) in year \(i\), \({e}_{j}\) is the entropy value of indicator \(j\), \({W}_{j}\) is the weight of the indicator \(j\) in urbanization indicator system, \({Y}_{ij}\) is the normalized value of the indicator \(j\) in the region \(i\) in this system, \(z\) is the number of indicators in this system.

Table 5 Evaluation system for comprehensive urbanization level.
Full size table

Analysis of the coupling relationship between biodiversity and urbanization

Coupling coordination degree model

CCD model is used to evaluate the relationship between biodiversity and urbanization of the 25 counties in the AYRSP.

$$C=2\sqrt{\frac{B\times U}{{(B+U)}^{2}}}$$
(19)
$$T=aB+bU$$
(20)
$$D=\sqrt{C\times T}$$
(21)

where \(C\) represents the coupling degree between biodiversity and urbanization level, \(B\) is the evaluation value of BI, \(U\) is the evaluation value of UI, \(T\) is the evaluation value of BI and UI, \(D\) is the coupling coordination degree, which is in the range of [0, 1]. Two systems are of equal importance, thus, a = b = 0.5. Referring to some studies29,31, the classification of CCD is shown in Table 6.

Table 6 Classification of coupling coordination degree.
Full size table

Relative development degree model

In order to further determine the lead or lag relationship between urbanization and biodiversity on the county scale, the study calculated their relative development degree. Based on the results of other studies29,31, the relative development types are categorized as follows (Table 7).

Table 7 Classification criteria of relative development types.
Full size table

Grey relation degree model

In order to study the influencing factors of coupling coordination development between biodiversity and urbanization on the county scale, the grey relation model is used to measure the correlation degree of the indicators in two systems to the coupling coordination degree. This study took the sequence of the coupling coordination degree between biodiversity and urbanization as the reference sequence, and took each indicator’s value of two systems as the comparison sequence. The formula is as follows35:

$${\xi }_{i}(k)=\frac{\underset{i}{\text{min}}\underset{k}{\text{min}}\left|{x}_{0}\left(k\right)-{x}_{i}\left(k\right)\right|+\rho \cdot \underset{i}{\text{max}}\underset{k}{\text{max}}|{x}_{0}\left(k\right)-{x}_{i}\left(k\right)|}{\left|{x}_{0}\left(k\right)-{x}_{i}\left(k\right)\right|+\rho \cdot \underset{i}{\text{max}}\underset{k}{\text{max}}|{x}_{0}\left(k\right)-{x}_{i}\left(k\right)|}$$
(22)
$${\xi }_{i}=\frac{1}{n}\sum_{k=1}^{n}{\xi }_{i}(k)$$
(23)

where \({\xi }_{i}(k)\) represents the grey relation coefficients between the degree of coupling coordination and the indicators of the two systems in year \(k\); \({x}_{0}\left(k\right)\) is the normalized value of the coupling coordination degree of two systems in the year \(k\); \({x}_{i}\left(k\right)\) is the normalized value of the internal indicators of two systems in year \(k\) (\(i\)=1,2…,7); ρ denotes the coefficient of variation, which usually takes the value of 0.5; \(n\)=4; \({\xi }_{i}\) represents the grey relation degree between the coupling coordination degree and the influencing indicators of two systems. The more \({\xi }_{i}\) tends to 1, indicating that the greater the influence of the indicator on the coupling coordination degree. Referring to the related research35, the classification of the grey relation degree is shown as Table 8.

Table 8 Classification criteria of the grey relation degree.
Full size table

Results

Evaluation of biodiversity

Spatio-temporal change in BI

The classification criteria of BI primarily refers to the standard for the assessment of regional biodiversity55 and similar studies50, and was combined with the physical and natural geography of the AYRSP. The biodiversity status was then divided into four ratings (Table 9).

Table 9 Biodiversity indicator rating criteria.
Full size table

The BIs in the AYRSP were 0.522, 0.506, 0.504, and 0.509 in 2015, 2017, 2019, and 2021, respectively, with an average BI value of 0.510; therefore, a decreasing trend was observed. Among them, the fastest declining trend in BI was observed in 2015–2017, with steady development after 2017, and a trend of recovery in 2021. Due to the influence of climate, the average precipitation in Shandong amounted to 988.2 mm in 2021, an increase by 53.2% compared to previous years. Therefore, vegetation, water, and soil moisture increased sharply compared with previous years. Simultaneously, the flow of rainwater transported soil nutrients to the subsoil. This cyclic effect maintains the balance of the ecosystem and plays an important role in maintaining biodiversity.

In terms of spatial distribution (Fig. 4), the spatial differentiation of biodiversity in the AYRSP was obvious from 2015 to 2021, presenting a pattern of “low in the western, and high in mid-southern and eastern regions”. According to statistics, Area I accounts for 24.31–27.56%, most of which is concentrated in Mount Tai and parts of the Yellow River Delta. Mount Tai and the Yellow River Delta are the KBAs in the lower Yellow River region. These areas are rich in arbor forests and shrubs, densely covered with river networks and diverse ecosystem types, and are rarely disturbed by human activities. Most are protected areas or forestry management areas, with the distribution of rare plants and animals, such as the Ciconia boyciana, Mergus squamatus, Grus leucogeranus, Nipponia nippon, Glycine soja, and Glycyrrhiza uralensis. Area II accounts for 27.29–31.48% of the AYRSP and these areas are primarily dominated by shrubs, other forests, grasslands, and water bodies, with high habitat quality and rich biodiversity. Area III accounts for 23.28–29.19%, which is primarily dominated by large areas of farmland and roads; therefore, landscapes are fragmented, and human interference is frequent. Area IV includes primarily concentrated built-up areas of urban and bare land, where the intensity of human activities is high, ecosystem type is fragile, and habitat quality is poor, accounting for 17.68–20.65% of the total area.

Fig. 4
figure 4

Spatio-temporal distribution patterns of biodiversity in the AYRSP on the grid. Generated in the ArcGIS 10.6 software (www.esri.com).

Full size image

The change of BI on the county scale is shown in Fig. 5. The five counties with the highest BI were Changqing (0.633), Licheng (0.632), Pingyin (0.578), Zhangqiu (0.565) and Hekou (0.530). The large mountainous areas of these counties and rich forest resources provide a good basis for biodiversity. Hekou is located north of the Yellow River Delta and has a dense water network, abundant wetland resources, a large area of high-quality grassland, and abundant nature reserves. The five counties with the lowest BI values were Mudan (0.469), Juancheng (0.461), Dongming (0.455), Yuncheng (0.454) and Yanggu (0.439). These counties are located in western Shandong, which is dominated by cropland and construction land with a single land cover. Additionally, their ecology is fragile and biodiversity development is limited.

Fig. 5
figure 5

Evolution of biodiversity in the AYRSP on the county scale. (a–d) Spatio-temporal distribution patterns of biodiversity in the AYRSP on the county scale; (e) the situation of the average value of BI in the AYRSP.

Full size image

Evaluation of sustainable development level of biodiversity

The KBAs contribute significantly to biodiversity conservation and sustainability. They are additionally important for evaluating SDG 15.1.2. The score was calculated according to the proportion of KBAs—27.0 points in 2015, 24.4 points in 2017, 24.3 points in 2019, and 25.9 points in 2021. The comprehensive scoring was 25.4 points, which was within the orange grade. This indicates poor implementation of SDG 15.1.2 in the AYRSP from 2015 to 2021. There are still some challenges to reaching the optimal level.

The distribution of KBAs in the 25 counties along the Yellow River of Shandong from 2015 to 2021 is presented in Fig. 6. Licheng (70.1 points) and Changqing (60.1 points) were at the green levels. Six counties were at the yellow level. The largest number of counties were at the orange level. Only Dongming (7.3 points) and Yanggu (8.2 points) were found to be at the red level, which indicates the worst implementation of SDG 15.1.2. Dongming has a small habitat area because of limited natural conditions. Moreover, the economy is relatively backward and infrastructure is inadequate. Dongming should, therefore, continuously strengthen source prevention measures to improve environmental quality and create favorable conditions for the survival of organisms. Yanggu is also limited by natural conditions. Industrial development has caused significant damage to the ecological environment. The environmental monitoring of Yanggu should, therefore, be strengthened to improve the quality of its ecological environment. Other counties should additionally make full use of the land and plant woodlands or grasslands. They should strive to expand green spaces in the urban areas, country, and habitat areas to make biodiversity more sustainable.

Fig. 6
figure 6

Spatio-temporal distribution of KBAs in the AYRSP on the county scale. (a) Overall change in trend of KBAs in 2015–2021; (b) spatial distribution patterns of KBAs in the AYRSP on the county scale.

Full size image

Evaluation of comprehensive urbanization

Referring to the existing studies31 and the actual situation in the AYRSP, the grading criteria of comprehensive urbanization indicator are determined as—lowest-level urbanization area (0.16 > UI ≥ 0), low-level urbanization area (0.24 > UI ≥ 0.16), medium-level urbanization area (0.32 > UI ≥ 0.24), high-level urbanization area (0.4 > UI ≥ 0.32), and higher-level urbanization area (1 ≥ UI ≥ 0.4). This study explored the characteristics of spatiotemporal changes of comprehensive urbanization level in the AYRSP (Fig. 7).

Fig. 7
figure 7

Spatio-temporal distribution patterns of comprehensive urbanization in the AYRSP on the county scale.

Full size image

The comprehensive urbanization level in the AYRSP continued to improve and develop during 2015–2021. The UI increased from 0.294 to 0.364, with an average of 0.328, and a growth rate of 24%. The degree of comprehensive urbanization of the 25 counties improved at the county scale. In terms of overall growth rate, the counties with the highest values were Juancheng (60.80%) and Dongming (60.51%), whereas the counties with a low growth rate of comprehensive urbanization were Dongping (2.08%) and Kenli (2.96%). Although the level of urbanization varies greatly, the speed of urbanization development of some high-level urbanization areas has been lower than that of medium- and low-level urbanization areas in recent years, which have shown rapid development.

In terms of the spatial distribution, the pattern observed was “high in the east and middle, low in the west.” Jinan and Dongying City were the “first echelons” in the development of urbanization in the AYRSP. In addition, the area of medium-level urbanization in the 25 counties along the Yellow River of Shandong continuously expanded, the area of low- and lowest-level urbanization gradually shrunk.

Coupling system and its influencing factors between biodiversity and comprehensive urbanization

Spatio-temporal change in the coupling coordination degree between biodiversity and comprehensive urbanization

The coupling degree of the 25 counties exceeded 0.80 during 2015–2021. The CCD values were 0.612, 0.627, 0.621, and 0.643 in 2015, 2017, 2019, and 2021, respectively, indicating the moderate-coordination phase. In terms of the growth rate, the CCD increased by 5.2%. In short, the coupling coordination level between biodiversity and urbanization showed an upward trend.

In terms of spatial distribution (Fig. 8), the eastern coastal and central regions had higher values compared to the southern and western regions of Shandong. Spatial differences in the CCD were caused by differences in economic development and habitats. First, the CCD of Jinan and Dongying City represented the highest levels, among which Huaiyin entered a high-coupling phase in 2019. Second, the number of counties in the moderate-coordination phase gradually increased to 15. Third, most counties in the west presented the lowest CCD. Although the counties in the west had removed the uncoordinated phase, they were still in the lowest position among the 25 counties along the Yellow River, particularly Liaocheng and Heze.

Fig. 8
figure 8

Spatio-temporal change in the coupling coordination between biodiversity and comprehensive urbanization in the AYRSP on the county scale.

Full size image

Spatio-temporal change in the relative development degree between biodiversity and comprehensive urbanization

The spatio-temporal distribution of the relative development degree of biodiversity and comprehensive urbanization is presented in Fig. 9. Only Huaiyin and Tianqiao presented lag types of BI before 2019. In 2019, Dongying entered the stage of lag types of BI. In 2021, Licheng and Bincheng entered the stage of BI-UI synchronous development type. In terms of quantity change, the number of lag types of UI decreased from 22 to 20, the number of BI-UI synchronous development types increased from 1 to 2, and the number of lag types of BI increased from 2 to 3. Although Huaiyin and Tianqiao have the high-level urbanization and CCD, they are the old urban areas of Jinan. These regions have highly concentrated populations, insufficient resource-carrying capacities, and prominent environmental problems that restrict the development of biodiversity. Other counties should, therefore, learn from the experience of urbanization construction and concentrate on improving the quality of urbanization development, while being aware of the trend of biodiversity degradation.

Fig. 9
figure 9

Spatio-temporal change in the relative development of biodiversity and comprehensive urbanization in the AYRSP on the county scale.

Full size image

Influencing factors of coupling coordination system

To study the influencing factors of the coupling coordination system between biodiversity and comprehensive urbanization in different regions, the grey relation degree model was used to evaluate the correlation between each indicator and the CCD, and ranking was given (Table 10). In 2015–2021, the grey relation degree between the influencing indicators and the overall CCD was above 0.6 in most counties along the Yellow River of Shandong, demonstrating high relevancy. In the biodiversity system, the grey relation degree of \({S}_{p}\) and HQI was the highest. This indicated that the resource endowment of ecosystems and habitats creates favorable space and carrier for the construction of urbanization, which can promote the coordinated development between biodiversity and urbanization. In the urbanization system, the HAILS and shares of secondary and tertiary industries presented the highest degree of grey relation. This meant that human activities and economic development have a profound impact on biodiversity and urbanization. In order to achieve sustainable biodiversity and urbanization conservation, a series of measures can be taken to reduce the negative impact of human activities. According to the ranking, the three major factors affecting the coupling coordination development of the two systems were \({S}_{p}\), HAILS, and shares of secondary and tertiary industries. They had a correlation coefficient of more than 0.9, which was the highest relevancy. In general, the grey relation degree of indicators in the urbanization system was higher than that of indicators in the biodiversity system, indicating that urbanization had a greater impact on the coupling coordination system. This means that a large amount of funds is used to improve and protect the ecological environment and biodiversity with the gradual development of urbanization, which can promote the coordination and sustainability of biodiversity and urbanization to a large extent.

Table 10 Grey relation ranking between each indicator and the CCD.
Full size table

Discussion

Biodiversity monitoring and conservation

Biodiversity monitoring is the basis for implementing biodiversity conservation and restoration. The construction of indicator system is the pivotal focus. Three principles have been considered in the selection of indicators—the typicality of model and calculation method, the interconnections among different levels of biodiversity, and the applicability for the fine scale of county. Based on this, this study selected 6 indicators that are closely related to biodiversity to construct biodiversity assessment methods. Compared with that of Peng et al.19, this study extended the evaluation of landscape diversity into the system, which is essential for maintaining biodiversity and ecological balance. Due to the difference of regional ___location, biological condition, data acquisition, and other factors, diverse indicators appeared in different studies56. For example, Li et al.20 constructed an evaluation indicator system of biodiversity including SHDI, HQI, NDVI, NPP, NDVI, and other indicators to dynamically monitor biodiversity in Maduo County, Qinghai province. However, the impact of the shape complexity in landscape patch was ignored, which was solved by LSI in this study. By quantifying the diversity and complexity of landscape patterns, SHDI and LSI provided an important basis for assessing biodiversity status, predicting the trend of biodiversity change, and formulating strategies of biodiversity conservation. However, the impact of landscape diversity on biodiversity is complex. On the one hand, complex and diverse landscapes may mean more habitats and effects of ecological edge, which contribute to increasing the species diversity and habitat, especially in areas strongly disturbed by human activities. On the other hand, the landscape fragmentation may lead to the fragmentation and loss of habitats, posing a threat to the survival and migration of species. Therefore, it is necessary to fully consider the impact of landscape fragmentation when protecting and restoring biodiversity and take corresponding measures to reduce its negative impact on species survival and migration. Compared with indicators of Yang et al.13, this study not only added LSI, but also supplemented the influence factor of WET. Suitable humidity conditions can not only maintain the balance of the ecosystem, but also promote the interaction and symbiosis between organisms. Therefore, these indicators selected in this study are scientific with mechanism basis and representativeness; the same number of indicators at each biodiversity level ensure the fairness of the construction of biodiversity indicators; meanwhile, the indicator system based on remote sensing guarantees the easy data acquisition, simple calculation, and good operability. In addition, the Chinese government attaches great importance to the implementation of the 2030 Sustainable Development Goals, which was listed as a national strategy in the 13th Five-Year Plan. The study innovatively identified KBAs in the AYRSP and quantified the result of SDG 15.1.2. Therefore, the evaluation method of biodiversity in the study is conducive to the local evaluation of SDG 15.1.2 to advance the reality of the SDGs.

Compared with the previous studies on biodiversity and habitat quality, this study used geographic national conditions monitoring data as the basic data source, which had the characteristics of high spatial resolution, strong spatiotemporal continuity, and a good current situation57. According to relevant studies41,42 and the present situation of the AYRSP, 15 landscape types were classified with more details. Different forest types differ greatly in structure and have different effects on evapotranspiration. Therefore, the subdivision of the main landscape types, such as forests and grasslands, could better reflect the differences in the impact of different vegetation types on the habitat, which was more advantageous for understanding the habitat suitable for biological survival or potential biological distribution. Subdividing construction land into four categories better demonstrated the impact of different building type and structure on habitat quality. In summary, the 15 landscape types were effective in elaborating local differences in habitat quality and biodiversity in the AYRSP. According to the relevant research58, there are higher species richness in the Dongping County and some counties of Dongying (Dongying, Hekou, Kenli, and Lijin) and Jinan City (Pingyin, Changqing, Huaiyin, Licheng, and Zhangqiu). The BI values calculated in this study were basically consistent with the previous research and actual situation of biodiversity in Shandong Province. Additionally, the study found that human activities such as the expansion of built-up areas and loss of natural habitats have threatened the development of biodiversity. Therefore, comprehensive measures are needed to mitigate the negative impact of human activities.

Analysis of urbanization level

Compared with previous urbanization studies, this study introduced the concept and index of HAILS for land urbanization. First, an urban ecosystem is a complex integrated system. In the process of urban development, the greater the HAILS, the more economic and social activities are present, and higher is the urbanization. Second, human activities are also the main driving factors of biodiversity degradation. From the perspective of social economy and land structure development, studies on the coupling coordination relationship between biodiversity and urbanization with the introduction of human activity intensity factors can better protect biodiversity. Additionally, the combination of remote-sensing and statistical data provides new ideas for the construction of urban systems. The study was found that the level of comprehensive urbanization continuously improved, which is consistent with the results of relevant studies29,31. Furthermore, the study found that the development speed of urbanization in high-urbanization areas was lower than that in medium- and low-level urbanization areas. Although the development speed of medium- and low-level urbanization areas is accelerating, some problems are still present, such as imbalances between urban infrastructure construction and population urbanization level, fierce contradictions between the deteriorating ecological environment and sustainable urbanization construction, and insufficient impetus for sustainable development. Such areas should accordingly change their development mode by simply pursuing urbanization speed and focusing on improving the quality and efficiency of urbanization.

The interaction and coupling effects of biodiversity and urbanization

In exploring the complex and critical relationship between biodiversity and urbanization, it is crucial to choose the suitable study methodology. Currently, the studies on the interaction between ecology and urbanization are various. Some scholars have explored the correlation between urbanization level and ecosystem services by using mathematical, statistical or spatial analysis models based on a single factor59,60, but such studies often lack in-depth analysis of the spatial differentiation of coupling coordination. In contrast, some other scholars, such as Sun et al.42 used CCD model to estimate the interactive coupling relationship between habitat quality and urbanization level in the lower Yellow River from 2000 to 2020, and found that although the two systems were in the basic coordination stage, the coupling coordination degree were constantly improving, which was also roughly similar to the results of this study. On this basis, this study innovatively introduced the relative development degree model to deeply analyze the leading-lagging relationship between biodiversity and urbanization in different time periods, providing a new perspective for understanding the dynamic changes between the two. At the same time, in view of the problem that previous studies paid little attention to the influencing factors of the coupling coordination system2,29,32,42, this study identified the key factors affecting the coupling coordination system by adopting the grey relational degree model with strong applicability, simple calculation, and strong interpretation, so as to make the study of the coupling system more complete and reasonable.

During the study period, the CCD between biodiversity and urbanization presented an upward trend, which indicated that the interaction between them in the AYRSP was very close. With the development of the social economy, people are paying more attention to protecting the ecosystem of the Yellow River Basin and promoting urbanization by improving resource utilization efficiency and changing lifestyles. In terms of the results of the relative development degree, the AYRSP still has a long way to go in the degree and quality of urban development. The study of the coupling relationship between biodiversity and urbanization in the AYRSP is a reminder that the coordination between biodiversity and economic society to improve the resilience of local ecosystems should be prioritized. Land use and urban forest construction should be planned as a whole to reduce the risk of species extinction, and ecosystem or landscape methods should be applied to reduce the pressure of land on forests and other natural ecosystems to avoid further habitat loss and create conditions to restore more habitats54.

Previous studies29,54,61 have seldom explored the influence of internal factors on the entire coupling coordination degree system. This study used the grey relation model to explore the main driving factors in the development of coupling coordination. First, the degree of correlation between the landscape diversity index and the coupling coordination system was minimal. Some studies62 have found that the relationship between biodiversity and the landscape shape index is inconsistent, which may be related to landscape types and the sampling scale of different study areas. It is, therefore, necessary to further explore the different effects of landscape type and sampling scale on biodiversity and coupling coordination systems. Second, the urbanization indicators presented a higher grey relation than the biodiversity indicators. In the future, the focus should be on improving the quality of urbanization to promote coordination between biodiversity and urbanization. Through the process of environmental assessment and large-scale zoning covering biological diversity in infrastructure investment planning and development, such as transportation system design and management, the consideration of biological diversity is reflected to avoid its vulnerability. Measures to protect ecological connectivity, such as underground passages and green infrastructure, promote the transition to a more sustainable urbanization model63.

In summary, the study aimed to establish an effective modeling method of coupling system to reveal the characteristics and driving mechanism of coupling coordination between biodiversity and urbanization, so as to provide theoretical support for promoting biodiversity conservation and urbanization development. By identifying key areas where biodiversity is rich or threatened, the system of natural reserve can be optimized to more effectively formulate ecological restoration plans and conserve biodiversity. Through scientific planning of urban spatial layout, the negative impact of urbanization on biodiversity can be reduced, and the coordinated development of urbanization and ecological protection can be realized. Additionally, the study can also provide useful ideas for similar research in other geographical regions. However, the specific coupling mechanism and influencing factors need to vary from region to region. For example, biodiversity in the upper reaches of the Yellow River are mainly mountainous ecosystems, which may be more dependent on natural factors such as topography, climate, and soil type; while urbanization may be more significantly influenced by factors such as topographic constraints, transport conditions, and economic activities. Therefore, when applying the conclusions of this study to other regions of the Yellow River Basin, regional particularities need to be fully considered, and in-depth analysis and verification should be carried out in combination with local actual conditions. In the future, the research scope can be actively expanded to provide a more comprehensive and in-depth scientific basis for the sustainable development of the Yellow River Basin through comparative analysis of the coupling mechanism of biodiversity and urbanization in different regions.

Limitations

Although remote-sensing data have the advantages of wide coverage, diversity, and easy processing, field survey data should still be combined as validation and supplementation if conditions permit, especially the species number and distribution data, to make the monitoring results of remote sensing more accurate and reliable. In the evaluation of the sustainable development of biodiversity, the use of remote-sensing method to identify KBAs and establish databases is an effective means of addressing the limitations of the calculation of the existing SDG 15.1.2. However, there are still some shortcomings. SDG indicators originally corresponded to the national scale, and the use of the “red or green grading” method to indicate the good or bad results of indicators also corresponds to the standards at the large-scale level. Blindly connecting to international standards may not be suitable for the actual situation of small regions. The issue of classification standards after the quantification of biodiversity indicators should, therefore, be explored more deeply. Currently, studies on the coordination between biodiversity and urbanization at the county scale in China are still in the exploratory stage. Biodiversity and urbanization systems in different regions are affected by geographical ___location and policies. Deficiencies inevitably occur during the construction of indicator systems. For example, relevant factors were not fully considered because of the influence of indicator acquisition. In the future, we will focus on the internal links between indicators and try to ensure that each indicator has its own unique importance. According to the evaluation results and the actual situation, we will also regularly evaluate and adjust the indicator system, and constantly optimize the construction of the indicator system. The selection of research scale may have a certain impact on the results of indicator construction. In the future, data at different scales can be combined for comprehensive analysis, which can reveal more levels of ecological processes and patterns. For example, it is possible to capture major ecological features at the macro scale and focus on the dynamic changes of specific ecological patches at the micro scale. This study only examined the spatial differentiation law of coupling coordination in 2015–2021. The research period was relatively short; this should, therefore, be extended in future studies. Furthermore, the specific driving mechanisms and path analyses should be further investigated. In addition, the degree of influence of the same influencing factor on the coupling coordination degree of different counties was not discussed; future work could, therefore, take into account the combination of the geographically and temporal weighted regression (GTWR) model and GeoDetector to investigate the influencing mechanism, which helps discover the interaction of the significant single factor and dual factors, and exploring the spatiotemporal heterogeneity of different factors.

Conclusions

Based on multi-source data, this study used CCD, relative development degree, and grey relation degree models to identify the coupling coordination stage, lead-lag relationship, and the main factors influencing the coupling system between biodiversity and urbanization in the AYRSP from 2015 to 2021. Then, the interactive mechanism and collaborative development path of biodiversity and urbanization were deeply discussed. It is of great significance for monitoring and evaluating the dynamic development of biodiversity and urbanization and the implementation of the Sustainable Development Goals. The conclusions were as follows.

  1. 1.

    From 2015 to 2021, the average BI value was 0.510, which represented the biodiversity level was high in local areas of AYRSP. However, the comprehensive BI presented fluctuating and declining trends. The overall spatial pattern observed was “low in the west and high in the mid-southern and eastern regions,” in which Mount Tai and parts of the Yellow River Delta are the core regions of the KBAs. At the county scale, the biodiversity of Changqing and Licheng was determined to be the most abundant. From the perspective of the implementation of SDG 15.1.2, the implementation was poor, with a score of 25.4, which remained a challenge from the optimal green level classified by the SDSN organization. The capacity for the sustainable development of biodiversity, therefore, needs to be strengthened.

  2. 2.

    From 2015 to 2021, the comprehensive urbanization level in the AYRSP presented an upward trend, with a growth rate of 24%, demonstrating a distribution trend of “high in the east and middle, and low in the west.” Jinan and Dongying City were the “first echelons” in the development of urbanization in the AYRSP. Areas with moderate urbanization were expanding, whereas those with low urbanization were shrinking.

  3. 3.

    The coupling coordination level between biodiversity and urbanization in the AYRSP presented an upward trend from 2015 to 2021, belonging to a moderate-coordination phase. The spatial difference in CCD was obvious, being higher in the coastal and central regions, but lower in southern and western Shandong. According to the results of the relative degree of development, Huaiyin, Tianqiao, and Dongying were lagging types of biodiversity. Licheng and Bingcheng were BI-UI synchronous development types. Other counties were mostly lagging in the urbanization type.

  4. 4.

    The grey relation degree between most of the influencing indicators and the overall CCD was above 0.6 from 2015 to 2021, demonstrating high relevancy. According to the results of the grey relation degree sorting, the three major influencing factors affecting the coupling coordination development of the two systems were the percentage of habitat area, human activity intensity of the land surface, and shares of secondary and tertiary industries. Overall, the urbanization demonstrated a greater impact on the coupling coordination level.