Spatiotemporal Distributions of Multiple Land Use Functions and Their Coupling Coordination Degree in the Yangtze River Delta Urban Agglomeration, China

Abstract

In order to gain a comprehensive understanding of land system changes and regional sustainable development, it is crucial to explore the spatiotemporal distributions of multiple land use functions (LUFs). Therefore, herein, considering the Yangtze River Delta urban agglomeration (YRDUA) as the research object, we constructed an evaluation system based on the perspective of production–living–ecological (PLE) land functions. Furthermore, the coupling coordination model, kernel density curve, trend surface analysis, and spatial autocorrelation model were used to examine the spatial and temporal characteristics of LUFs and their coupling coordination and analyze the spatial clustering of the coupling coordination effect in the urban agglomeration from 2005 to 2020. The following results were obtained: The level of PLE functions and its coupling coordination degree in the YRDUA have been improved in the study period, and the distribution of high value areas of different functions is diverse. In terms of the spatial distribution of the coupling coordination degree, the high value areas of production function (PF)–living function (LF) is distributed in “clusters”, the PF–ecological function (EF) demonstrated a single-polarized development structure, and the LF-EF exhibited a multi-core structure. The coupling coordination of the LUFs demonstrated a “high in the east and low in the west” pattern in the east–west direction and an inverted “U” shape in the north–south direction. Moreover, both high-value and low-value areas exhibited a clustering phenomenon, with an evident spatial positive correlation. The results of this study can serve as a reference for the integrated socioeconomic development of the Yangtze River Delta region and the sustainable development of urban land use.

1. Introduction

Land use and cover change are crucial drivers of global change, with significant implications for ecosystems, global biogeochemistry, climate change, and human vulnerability [1,2,3]. In the context of economic globalization, human demand for land is ever increasing [4,5]. Therefore, a growing number of countries are facing the issue of scarce land resources and compact spaces required for human activities [6]. In this process, economic growth, social development, and changes in human activities all contribute to the multiple functions of land in a region and exhibit significant spatial differences and temporal variability [7]. Multifunctional land use (MLU), which can be described as integrating various functions on the same territorial units, has proved to be an effective way to coordinate limited land resources and multiple human demands, especially in densely populated countries and regions [8,9,10].

Land use functions (LUFs) are defined as the private and public goods and services provided by different land uses that summarize a region’s key economic, environmental, and social issues [11]. The concept of LUFs was first proposed in the European Union’s (EU) Sixth Framework SENSOR Project [12]. The term “multifunctionality” was coined by the Organisation for Economic Co-operation and Development (OECD) and the EU in their theoretical considerations on agricultural policy reforms [13]. With the advancement of “multifunctionality” studies, researchers have increasingly recognized that agriculture is not the only sector with multifunctional characteristics [14]. With the intensification of global land planning research, the concept of multifunctionality has been increasingly associated with land use [15]. The concept of MLU has evolved with time and has become the frontier of international sustainable research [16,17]. Owing to the variations among different research objects and disciplines, the classification methods also vary [18,19]. Based on the regional system of human–earth relationships, the land use system often includes three fundamental functions, i.e., the production function (PF), living function (LF), and ecological function (EF) [20]. PF is regarded as the driving force, EF as the foundation, and LF as the link of LUFs. The three aspects of LUFs are closely related and mutually transformed and, in combination, determine the overall efficiency of land use [21]. Therefore, in recent years, many countries have undertaken enormous efforts to improve land use efficiency to alleviate the conflict between people and land. However, in most countries, the conflict between the functions within the land use system remains entrenched [22,23]. To this end, the 18th National Congress of the Communist Party of China proposed to build a livable national land space with “intensive and efficient production space, moderate and livable living space, and beautiful ecological space” [24]. Moreover, the need to build a regional economic layout and territorial spatial system with complementary advantages and high-quality development was further emphasized at the party’s 20th National Congress. This indicates that coordinating the relationship between MLUs and exploring the synergy among different regions from the perspective of production, living, and ecology (PLE) can better provide certain policy guidance to promote people’s livelihoods and achieve sustainable regional development [25,26].

Research on LUFs in China started in the 1980s and is still at a preliminary stage. In general, the research framework and methodological system of LUFs have been initially formed. However, there is no consensus on the connotation of LUFs nor has a scientific and unified evaluation system been formed. Therefore, how to determine the connotation of LUFs from the perspective of system, establish an evaluation system based on the actual land use in the region, and validate the evaluation system is a key point of future relevant research. In addition, the impact of the coordination of LUFs on regional development has attracted the attention of scholars, and the research direction is gradually evolving toward the coordination of multifunctional relationships. The existing research mainly focuses on the “point” study of developed cities or megacities. Then, it is important to gradually strengthen the development from “point” to “surface” and pay attention to the spatial and temporal evolution characteristics of the coupled coordination relationship and synergistic relationship between different functions of land use in urban clusters or specific regions for territorial spatial planning.

Coupling, which originates from physics and was gradually introduced into the field of geography, is a phenomenon in which two or more subsystems interact and influence each other via miscellaneous interactions in the whole physical system. Coordination, including the good interaction between the various departments in a system, can explain the sustainable development of the system. A complex nonlinear coupling relationship usually exists between the elements or subsystems. The coordinated development of a system refers to the harmonious development of its subsystems and the relationship between them [27]. Therefore, the economic significance of obtaining the variable of the coupling coordination degree is that it can help us better understand the output effects of the multiple systems when they work together and develop together. This model has been widely used by researchers to measure the coupled coordination of land use systems with urbanization [28], livelihood efficiency [29], population growth, and housing supply [30] to provide a theoretical basis for urban development planning. Given that there exists a coupling relationship among the PLE functions with mutual promotion and coordination, the coupling coordination degree model may prove effective for the quantitative measurement of both the coupling and coordination relationships among production, living, and ecological space functions [31]. However, the application of this model to measure the internal relationships of LUF systems from this perspective remains poorly explored. Moreover, most existing studies on the evolutionary patterns of the coupled coordination of LUFs lack a consideration of the coupling coordination differences within the study area. In this direction, kernel density estimation and trend surface analysis will allow further analysis. The kernel density estimation curve can portray the distribution pattern, location, and extension of the random variables, so it can be used to examine the time-series evolutionary variation characteristics of the coupling coordination degree of LUFs. In addition, the trend surface can be used to fit a mathematical surface to illustrate the spatial trends and distribution patterns of the observed values of geographic elements over a large spatial span, which will then be used to further analyze the overall spatial divergence trends of the coupling coordination within the study area [32].

As one of the three major urban agglomerations in China, the Yangtze River Delta urban agglomeration (YRDUA) encompasses some of the most important modern service industries and advanced manufacturing industries in the world. Since the reform and opening up of China, the complex land use types of the YRDUA make the temporal and spatial evolution of the PLE functions more acute, which in turn leads to significantly enhanced regional changes in natural, economic, and social factors in the region. Therefore, it is an urgent problem to accurately and comprehensively reveal the interaction between the PF, LF, and EF of the YRDUA and analyze its spatial and temporal evolution characteristics in order to realize the coordinated development of the PLE functions. All in all, in this study, considering the municipal level as the research unit, we aimed to explore the coordination characteristics of each LUF from 2005 to 2020 using the coupling coordination degree model, kernel density curve, and trend surface analysis. Furthermore, we applied the spatial correlation model to analyze the spatial agglomeration characteristics of the study area. Finally, we put forward some policy implications based on research conclusions. This study will provide a scientific basis for sustainable development by considering land management planning and land resource use in urban agglomerations.

2. Materials and Methods

2.1. Study Area

The YRDUA, which is located in the “T” zone along the river coast of China (115°46′–123°25′ E, 23°46′–29°20′ N), encompasses the Jiangsu, Zhejiang, and Anhui provinces and Shanghai and includes 26 cities in total [33]. It predominantly exhibits a subtropical monsoon climate, with complex ecosystem types and diverse land cover (Figure 1). As the most developed urbanized area in China, the YRDUA plays a crucial strategic role in China’s overall modernization drive and opening-up pattern. However, the rapid economic development in urban agglomerations is also accompanied by outstanding problems such as strained land supply and demand, large regional disparities, and severe environmental pollution. Therefore, the Communist Party of China’s 20th National Congress report clearly states the need to continuously promote the development of the Yangtze River Economic Belt and the integrated development of the Yangtze River Delta. The orderly and coordinated development of PLE functions, faced with the contradiction between land supply and demand and the optimization of land use structure, is a pressing issue of land spatial planning.

2.2. Data Sources

In this study, the economic, ecological, demographic, and social development data of 26 cities in the Yangtze River Delta from 2005 to 2020 were collected from the Statistical Yearbook, Water Resources Bulletin, National Economic and Social Development Bulletin, and other relevant government work reports of each province and city in the corresponding year. Missing data for various areas for some of the years were obtained by linear interpolation of recent years.

2.3. Calculation Methods and Analysis Models

2.3.1. Production–Living–Ecological Land Use Function Classification System

The evaluation of LUFs is fundamental for research on and the analysis of multifunctional utilization and activities of land. Based on the PLE functions of land space and considering the characteristics of economic and local environmental development and the level of living conditions of residents in the study area, this study divides the LUFs into PF, LF, and EF. The indicator system is the basis of multifunctional evaluation of land use, and the selection of indicators has an important impact on the evaluation results. Therefore, the evaluation indexes especially emphasize the relevant elements reflecting the characteristics of the three dimensions of PF, LF, and EF. At present, the LUF index system is being enriched and developed with the depth of research. According to the principles of scientificity, data availability, and accuracy, and referring to the existing research results [12,34,35], the corresponding secondary functions and specific indicators were selected to construct a multifunctional evaluation system for land use (

Table 1. Assessment of the index system for production–living–ecological land use functions.

Land Use Functions	Sub-Functions	Indicators	Units	Weight
Production Function (0.422)	Economic development	GDP per capita	RMB per area−1	0.172
		Density of economy	RMB per capita−1	0.436
		Industrial structure	%	0.084
	Agricultural production	Grain yield per area	Kg ha−1	0.116
		Land reclamation rate	%	0.192
Living Function (0.361)	Life bearing	Construction land area	10,000 km2	0.378
		Road network density	%	0.298
	Social security	Number of beds per 10,000 people in hospital	10,000 beds	0.142
		Disposable income of urban residents	RMB	0.160
	Employment support	Urban registered unemployment rate	%	0.022
Ecological Function (0.217)	Ecological conservation	Park green space per capita	m2	0.112
		Total water resources per area	m3	0.271
	Resource supply	Cultivated area per capita	m2	0.533
	Environmental purification	Discharge of industrial wastewater	ton/km2	0.064
		Emissions of industrial sulfur dioxide	ton/km2	0.020

). First, PF considers the production of various products or services by diverse social production activities, taking land as the carrier and object of human labor. It primarily includes economic development and agricultural PF. The GDP per capita, density of economy, and industrial structure represent the economic development function, whereas the agricultural PF is characterized by grain yield per area and land reclamation rate. Second, LF refers to the ability of the land to provide all kinds of life guarantees, such as education, medical treatment, health, and transportation, in the process of human survival and development, including three sub-functions: life bearing, social security, and employment support functions. The life bearing function is characterized by construction land area and road network density. The social security function is the basis for maintaining people’s living needs, mainly in terms of the number of hospital beds per 10,000 people and the density of the road network. The employment support function is characterized by the urban registered unemployment rate. Finally, EF supports the development of human economy and society and limits the excessive use of resources. It is the basis and support of the multifunctional system of land and primarily involves the functions of ecological preservation, resource supply, and environmental purification. The ecological conservation function is expressed by park green space per capita and total water resources per area. Cultivated area per capita is used to characterize the land resource supply function, and industrial wastewater discharge and industrial sulfur dioxide emissions are utilized to express the environmental purification function.

2.3.2. Entropy Method

To eliminate the subjectivity in the process of weight determination, we herein adopted the method of combining subjective and objective weights to calculate the weight of the index system. We applied the expert scoring method to determine the weight of the first-level functions and the entropy method to determine the weight of each specific index. Based on the opinions collected from five experts, including two from the field of sociology, two from the field of land resources management, and one from the field of economics, we assigned weights of 0.422, 0.361, and 0.217 to the PF, LF, and EF, respectively. Entropy was originally introduced into information theory by Shannon. The entropy method relies on the physical characteristics of indicators to objectively assign weights and has been widely used in data analysis [36]. The specific steps of this method are as follows:

First, as the data should be standardized and dimensionless, this study used the polar difference method to standardize the original data [37]:

Positive indicator:

$$X_{ij}=\frac{x_{ij}-\min (x_{ij})}{\max (x_{ij})-\min (x_{ij})}$$

(1)

Negative indicator:

$$X_{ij}=\frac{\max (x_{ij})-x_{ij}}{\max (x_{ij})-\min (x_{ij})}$$

(2)

where x_ij denotes the value of indicator i in year j, X_ij denotes the standard value of x_ij, and max(x_ij) and min(x_ij) are the maximum and minimum values of indicator i in all years, respectively. Thus, the index values ranged from 0 to 1.

Ratio and entropy calculation:

$$Y_{ij}=\frac{X_{ij}}{{\displaystyle \sum X_{ij}}}$$

(3)

$$e_j=-\frac{1}{{\ln }_j}{\displaystyle \sum (Y_{ij}\times \ln Y_{ij})}$$

(4)

where Y_ij denotes the normalized data matrix, and e_j denotes the indicator entropy value.

Weight calculation:

$$W_i=\frac{1-e_j}{{\displaystyle {\sum }_{j=1}^n(1-e_j)}},i=1,2,3,\dots ,n$$

(5)

where W_i denotes the objective weight value of the index i.

Comprehensive score calculation:

$$S_{ij}={\displaystyle \sum _{i=1}^n{w}_j\ast X_{ij}}$$

(6)

where S_ij denotes the comprehensive score of the index i.

2.3.3. Coupling Coordination Model

Coupling coordination is used to measure the degree of association, perception, and dependence between independent systems and is being widely applied to measure the degree of interaction between systems. Accordingly, in this study, we constructed a coupling coordination degree model to calculate the coupling coordination degree index for multiple LUFs among the cities in the YRDUA. The equations are as follows [38]:

$$C=3\times {\left[\frac{P_i\times L_i\times E_i}{(P_i+L_i+E_i)}\right]}^{\frac{1}{3}}$$

(7)

$$T=\alpha P_i+\beta L_i+\chi E_i$$

(8)

$$D={(C\times T)}^{\frac{1}{2}}$$

(9)

where C denotes the coupling degree of PF, LF, and EF, and the value range is [0, 1]. Further, the value of C is determined by the evaluation value of PF, LF, and EF. P_i, L_i, and E_i denote the comprehensive evaluation values of PF, LF, and EF, respectively. α, β, and χ represent the weights of these three functions in the coupling coordination model, respectively, where α + β + χ = 1. Based on the expert opinions and the weight calculation of each specific index in the relevant literature, this study hypothesized that the three functions are equally significant in the coupling coordination model. Therefore, the mean value of all the undetermined coefficients is 1/3. D denotes the coupling coordination index of urban LUFs, and t represents the coupling degree.

To further analyze the degree of interaction between two functions of multiple LUFs of the YRDUA, we constructed a coupled coordination model between the two systems. The equations are as follows [39]:

$$C_1=2\times {\left[\frac{P_i\times L_i}{(P_i+L_i)}\right]}^{\frac{1}{2}},T_1={\alpha }_{1}P+{\beta }_{1}L,D_1={(C_1\times T_1)}^{\frac{1}{2}}$$

(10)

$$C_2=2\times {\left[\frac{P_i\times E_i}{(P_i+E_i)}\right]}^{\frac{1}{2}},T_1={\alpha }_{2}P+{\beta }_{2}E,D_2={(C_2\times T_2)}^{\frac{1}{2}}$$

(11)

$$C_3=2\times {\left[\frac{L_i\times E_i}{(L_i+E_i)}\right]}^{\frac{1}{2}},T_1={\alpha }_{3}L+{\beta }_{3}E,D_3={(C_3\times T_3)}^{\frac{1}{2}}$$

(12)

where C₁, T₁, and D₁ denote the coupling, coordination, and coupling coordination degrees between PF and LF, respectively; C₂, T₂, and D₂ denote those between PF and EF, respectively; C₃, T₃, and D₃ denote those between LF and EF, respectively. Furthermore, α₁ and β₁, α₂ and β₂, and α₃ and β₃ refer to the weights of PF, LF, and EF, respectively, in the coupled coordination model; i.e., they are the coefficients to be determined. In general, α₁ + β₁ = 1, α₂ + β₂ = 1, and α₃ + β₃ = 1. Thus, the mean value of all the determined coefficients is taken as 0.5.

According to the actual situation of the study area and the existing research results, the coupling coordination degree is graded using the non-equidistant method. The specific grading criteria are shown in

Table 2. Classification standard of the coupling coordination degree.

Coupling Coordination Degree	Level	Coupling Coordination Degree	Level
0 < D ≤ 0.200	High conflict	0.501 < D ≤ 0.600	Barely coordinate
0.201 < D ≤ 0.300	Intermediate conflict	0601 < D ≤ 0.700	Low coordination
0.301 < D ≤ 0.400	Low conflict	0.701 < D ≤ 0.800	Intermediate coordination
0.401 < D ≤ 0.500	Slight conflict	0801 < D ≤ 1	High coordination

.

2.3.4. Nonparametric Kernel Density Estimation Curve

Kernel density estimation, as a nonparametric estimation method, is not subject to any assumptions and is used to study the characteristics of data distribution, utilizing the data sample itself. It offers the advantages of intuitive expression, conceptual simplicity, and ease of calculation and is often used to analyze dynamic sample distribution. Different types of kernel functions include Gaussian, triangular, and quadratic kernel functions. Herein, we selected the Gaussian kernel function commonly used in theoretical circles to analyze the dynamic changes in spatial distribution of coupling coordination degree of multiple LUFs in the YRDUA. The formula is as follows [40]:

$$f(x)=\frac{1}{Nh}{\displaystyle \sum _{i=1}^{N}K(\frac{x-x_i}{h}})$$

(13)

$$K(\frac{x-x_i}{h})=\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{1}{2}{(\frac{x-x_i}{h})}^2\right]$$

(14)

where f(x) refers to the kernel density function; N denotes the number of observations; K(-) denotes the kernel function; and h denotes the bandwidth, which is used to control the degree of smoothing. The optimal bandwidth expression, h = 1.06Se × N − 1/5 (where Se is the standard deviation of random observations), proposed by Silverman, is used to determine the value of h.

2.3.5. Trend Surface Analysis

The trend surface analysis fits the values of spatial sampling points defined by mathematical functions using global polynomials. Herein, it was used to express the spatial divergence trend of the coupling coordination degree of multiple LUFs in the YRDUA and is represented as follows [41]:

$$z_i(x_i,y_i)=T_i(x_i,y_i)+{\epsilon }_i$$

(15)

where (x_i,y_i) denotes the geographical coordinate, and ε_i denotes the residual, i.e., the deviation of the true value from the fitted value.

2.3.6. Spatial Autocorrelation

Spatial autocorrelation detects the dependence of the spatial distribution pattern of things and phenomena; determines the diffusivity, polarization, or randomness of the spatial distribution; and, thus, reveals the interaction mechanism between the attribute values of the research object. In this study, we used global autocorrelation (Moran’s I) to measure the overall spatial autocorrelation degree of the coupling coordination degree of multiple LUFs in the YRDUA. The calculation formula can be expressed as follows [42]:

$$I(d)=\frac{n{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{\omega }_{ij}(x_i-\overline{x})(x_j-\overline{x})}}}{{\displaystyle \sum _{j=1}^n(x_i-\overline{x}){\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{\omega }_{ij}}}}}$$

(16)

where x_i and x_j denote the values of variable x at adjacent paired spatial points in regions i and j, respectively; $\overline{x}$ denotes the mean value; ω_ij denotes the economic gap between provinces (i.e., the inverse of the mean difference between the actual per capita GDP of regions i and j during the sample period); n refers to the total number of spatial points. The value of Moran’s I coefficient is [−1, 1]. A value less than 0 indicates a negative correlation, a value equal to or close to 0 indicates no correlation, and a value greater than 0 indicates a positive correlation.

3. Result Analysis

3.1. Spatial Patterns of Production–Living–Ecological Land Use Functions

By constructing the multifunctional evaluation system of land use, we obtained the evaluation results of multiple LUFs in the YRDUA. Using the natural breakpoint method in ArcGIS, four key years (2005, 2010, 2015, and 2020) were selected. The evaluation results were divided into five levels: high, moderately high, moderate, moderately low, and low, and accordingly, the spatial distribution map of each function was drawn (Figure 2).

3.1.1. Evolution of Production, Living, and Ecological Functions

Figure 2 shows that the spatial distribution of the development level of each LUF in the YRDUA varies. In particular, the PF level gradually evolved from “high in the north and low in the south” in 2005 to “high in the east and low in the west” in 2020. Hefei and Chizhou exhibited gradual development from high to moderate and low levels. The PF level of most cities remained moderate, indicating that it has not changed significantly.

The spatial differentiation pattern of the development level of the LF demonstrated no significant change, with a slight local adjustment, showing a pattern of “high in the east and low in the west”. However, with the advancement of time, its functional differences exhibited a trend of first increasing and then decreasing. The LF level significantly improved from 2010 to 2020, and the LF of most cities reached the moderate level and above. In this urban agglomeration, Shanghai, Nanjing, Hangzhou, and Shaoxing demonstrated a high level, with the low-value areas being distributed in the surrounding cities in the west of the YRDUA, including Anqing, Chizhou, and Xuancheng.

The EF level demonstrated evident differentiation during the study period. In general, the high-value areas are mainly located at the western periphery. In 2005, the high and moderately high value areas of the EF were scattered at the edge of the YRDUA, such as in Yancheng, Chuzhou, Chizhou, Jinhua, Taizhou, Ningbo, and Zhoushan, and by 2015, they were concentrated in the southern region of the YRDUA, such as in Chizhou, Xuancheng, Hangzhou, and Jinhua. In 2020, the EF was dominated by medium and moderately low levels, exhibiting a trend of “high in the west and low in the east”. The EF level of Anqing and Chizhou was the highest. The closer the area to the eastern region, the lower the function level. Accordingly, Shanghai, Suzhou, and Taizhou exhibited a low EF level.

The value of the LUFs in the YRDUA increased during the study period, and the functional level differentiation was evident. Moreover, the regional difference showed a trend of first increasing and then decreasing with time. Although the overall level remained primarily beyond moderate, the moderately high level gradually evolved, demonstrating a scattered distribution in 2005 and a centralized distribution in 2020. In the urban agglomeration, Shanghai and Xuancheng consistently demonstrated the highest and lowest levels, respectively, and Chizhou exhibited an evolution from the highest to the lowest level. Furthermore, at the end of the study period, nine cities reached the moderately high levels. The LUFs of the central north–south part of the YRDUA exhibited fast development. Additionally, Hangzhou and Shaoxing, located in the south, also exhibited moderately high levels of LUF development.

3.1.2. Evolution of Sub-Functions of Land Use

This study further analyzed the secondary functions of land use and expressed their changes in four key years using radar charts. Figure 3 shows that there exist large differences in the development sub-function levels of land use in the YRDUA. In particular, the economic development and social security functions demonstrated an evident growth trend. Moreover, the development rate of the ecological conservation function exhibited a trend of first increasing and then decreasing. The development of the life bearing function remained well balanced, and the development of the employment support and environmental purification functions has been relatively slow. The agricultural PF exhibited a trend of first increasing and then decreasing, whereas the resource supply function continued to decline.

3.2. Coupling Coordination Relationships among Production–Living–Ecological Land Use Functions

We prepared a line chart according to the evaluation value calculated using the coupling coordination model and the time development trend (Figure 4). Moreover, based on the above-mentioned classification criteria and applying ArcGIS10.8 software, we generated a spatiotemporal distribution map of the coupling coordination degree between two functions and the total functions in the study area for years 2005, 2010, 2015, and 2020 (Figure 5).

3.2.1. Spatiotemporal Distribution of Coupling Coordination Degree between Two Functions

Figure 4 and Figure 5 demonstrate that the coupling coordination degree of PF–LF increased year by year, demonstrating a rapid rate of development. During the study period, the coupling coordination degree of PF–LF could be characterized by a pattern of “high in the east and low in the west”. In particular, the PF–LF coupling coordination in most cities remained in a low conflict state at the early stage, and Shanghai, which exhibited the highest coupling coordination level, demonstrated scarce coordination. By the end of the study period, the PL–LF coupling coordination levels demonstrated scarce coordination and slight conflict in most cities of this urban agglomeration. In Anqing, the PL–LF coupling coordination level evolved from highly to slightly conflicting. Moreover, Shanghai reached a high coordination level, whereas Nanjing, Wuxi, Suzhou, Hangzhou, and Shaoxing reached a low coordination level.

The PF–EF coupling coordination degree presented a zigzag upward trend, with the pace of development among cities being relatively consistent. In particular, in 2005, it showed a spatial distribution pattern of “high in the north and low in the south,” with the southern and northern regions being primarily dominated by low conflict levels. By 2020, the differences in the levels of the regions decreased, and overall, the urban agglomerations demonstrated a scarce coordination level. Only Chizhou, Xuancheng, Jinhua, and Taizhou exhibited slight conflict. In general, the PF–EF coupling coordination degree presented a single-polarized development structure centered on Shanghai.

The LF–EF coupling coordination degree exhibited a fluctuating upward trend, with a slow development rate. In terms of spatial distribution, the high-value areas presented a multi-core structure, with Hefei, Nanjing, Shanghai, Hangzhou, and Nanjing at the core. From 2005 to 2015, some cities developed rapidly and reached scarce and low coordination levels, with eight cities still exhibiting slightly conflicting levels. By 2020, the difference in the levels of the cities decreased, with low coordination level cities being scattered in the middle of the urban agglomeration. By this time, 21 cities exhibited a scarce coordination level.

3.2.2. Spatiotemporal Distribution of Coupling Coordination Degree of Land Use Functions

To further explore the urban variation characteristics of the LUF coupling coordination, Stata16.0 software was used to plot kernel density estimation curves for an in-depth analysis (Figure 6).

The curve exhibited a multi-wave peak in 2005, with the coupling coordination being evident at 0.25, 0.35, and 0.47. Moreover, the peak at 0.35 was significantly higher than the peaks at 0.25 and 0.47, implying a certain gap in the development of various parts of the YRDUA. However, the peaks in 2010 and 2015 decreased, and the wave width gradually increased. Compared with that in 2015, the coupling coordination degree in 2020 only increased slightly, and all of them showed a single-peak phenomenon in the last four years, indicating the increasing gap of LUF coupling coordination degree among cities. Simultaneously, the curve for 2020 showed a right-trailing phenomenon around 0.6, implying that the coupling coordination level in a small number of cities reached a higher level.

In terms of spatial distribution, the level of coupling coordination degree also exhibited some differences. In particular, the overall level in 2005 remained low, and most cities exhibited a low conflict level. Moreover, the number of cities exhibiting slight conflict increased between 2005 and 2015. Shanghai, belonging to the cities with low coordination, demonstrated the highest level of coordination, and only Anqing still presented moderate coordination. In 2020, Shanghai further evolved into a moderately coordinated city, with Nanjing, Changzhou, Zhejiang, and Shaoxing exhibiting low coordination levels and most cities exhibiting scarce coordination levels. Only Chizhou, Xuancheng, Tongling, and Taizhou demonstrated a low level of coupling coordination and slight conflict.

The trend surface analysis using ArcGIS software was conducted to explore the spatial trend changes of the LUF coupling coordination, the results of which are shown in Figure 7. As shown in the figure, the east–west trend line slopes to the upper right, and its slope tends to first increase and then decrease, indicating that the LUF coupling coordination degree in the eastern region of the YRDUA was higher than that in the western region. The variation in development among different cities exhibited a trend of first increasing and then decreasing. As time advanced, the trend surface gradually changed from a horizontal line to an inverted “U” shape in the north–south direction. This suggests that interregional development differences gradually emerged with the rapid development of the central region, and the LUF coupling coordination degree was significantly higher in the central region of the YRDUA than that in the southern and northern regions.

3.2.3. Spatial Correlation Pattern of the Coupling Coordination of Land Use Functions

To further reveal the spatial agglomeration characteristics, the Moran’s I value of the multifunctional coupling coordination degree of land use between different cities in the Yangtze River Delta from 2005 to 2020 was calculated using the spatial autocorrelation model (

Table 3. The global Moran’s I of coupling coordination degree of LUFs.

Year	Moran’s I	Z Value	Year	Moran’s I	Z Value
2005	0.220 **	2.953	2013	0.179 ***	5.605
2006	0.180 **	2.480	2014	0.460 ***	5.432
2007	0.249 ***	3.243	2015	0.474 ***	5.584
2008	0.292 ***	3.695	2016	0.477 ***	5.853
2009	0.332 ***	4.060	2017	0.480 ***	5.621
2010	0.284 ***	3.494	2018	0.432 ***	5.159
2011	0.462 ***	5.434	2019	0.411 ***	4.980
2012	0.479 ***	5.630	2020	0.303 ***	3.813

Note: ** p < 0.05; *** p < 0.01.

). The Moran’s I scatter plot (Figure 8) for years 2005, 2010, 2015, and 2020 was drawn to explore spatial correlation of the coupling coordination.

It is evident from Table 3 that the Moran’s I values of the coupling coordination degree were all positive from 2005 to 2020, and both the P and Z values passed the significance test, indicating that the coupling coordination degree between the various functions of land use showed a significant spatial positive correlation. Combined with the scattered distribution of Moran’s I, most of the observed values in the four years were clustered in the first and third quadrants. This means that the multifunctional coupling coordination of land use in the YRDUA exhibited certain spatial agglomeration, with certain agglomeration phenomena in high- and low-value areas.

4. Discussion

Currently, research on the development, utilization, protection, and governance of national land resources is emerging from the perspective of LUFs [43]. The quantification and understanding of the spatial and temporal dynamic evolution characteristics of global land use multifunctionality are becoming increasingly important for strengthening land resource management [44]. The YRDUA is a crucial international gateway in the Asia-Pacific region, and via an in-depth analysis of its land use multifunctionality and coupling coordination degree characteristics, the present study can provide a reference for the coordinated development of LUF at the urban agglomeration scale. This study serves as a guidance for the coordinated development of multiple LUFs on the scale of urban agglomerations.

Optimizing the land space is key to building a society with well-developed production, happy and affluent lives, and a friendly ecological environment and strengthening the interaction of the PLE functions. From this study, it is evident that the high-value areas for PF and LF, such as Shanghai, Suzhou, Wuxi, and Nanjing, demonstrate relatively low values for EF. On the contrary, the PF and LF of Anqing, Chizhou, Chuzhou, and other regions on the edge of the Yangtze River Delta exhibit low levels. However, these cities are high-value areas for the EF. These results are consistent with those of Shan et al. [45]. From the evolution characteristics of the sub-functions of land use, it is evident that the resource supply function in the EF of the entire YRDUA gradually decreased. This could be attributed to the extraordinary period of economic growth experienced by China since its economic reforms in 1978 because the economic recovery of most countries is dependent on the exploitation of natural resources [46]. The gradual acceleration of industrialization, trade, and urbanization stimulated the rapid development of China’s economy. During this period, the PF and LF of coastal cities, provincial capitals, and resource-based cities such as Shanghai, Nanjing, and Suzhou improved owing to their superior economic location, regional policies, administrative resources, and development foundation [47]. Simultaneously, the ecological land was damaged to a certain extent, leading to a negative impact on EFs [48]. Accordingly, the environment continues to deteriorate, which is consistent with the findings of Liu et al. [49]. Although this impact may not be visible in the short term [50], it is undeniable that the ecological problems caused by industrial development pose a significant challenge to sustainable development [51]. Thus, China has proposed an ecological protection “red line,” with continuous emphasis on the need to maintain and develop the land use EF.

The destruction of the ecological environment also affects the comprehensive development of cities for a long time. In this study, it is evident that the coupling coordination degree of PF–EF and LF–EF in the YRDUA has rapidly increased in recent years. Ren revealed that the climate change policies aimed at improving the area and quality of ecological lands are more conducive to the coupling development of the climate–economy nexus than those focusing on limiting carbon emissions [52], further indicating that the development of land use EFs can better stimulate the development of other fields. Moreover, we found that the level of coupling coordination degree of PF–LF in the YRDUA remains relatively stable with the evolution of time. However, the spatial distribution is more complicated. In 2020, five different levels of coupling coordination degree were noted, indicating its variation between different regions. Shanghai exhibited a high coordination level, and most of the fringe cities to the north and south demonstrated a moderate conflict level, with Chizhou and Xuancheng being at high conflict levels. The population and economic carrying capacity of these areas are low, which may be due to the low population and economic carrying capacity in these areas, and the uneven development of production and living space needs further coordination. Moreover, in terms of the coupling coordination degree of LUFs, the difference between the north and south directions was not evident at the beginning of the study period. Over time, this difference became more pronounced. By 2020, the coupling coordination level in the central region was significantly higher than that in the north and south regions, with the high-value areas being “grouped” in Shanghai, Nanjing, Wuxi, and their surrounding cities.

Natural resources in different regions can complement each other to a certain extent, and, as the resource elements of a region interact with its neighbors, a variety of linkages arise between these regions [53]. After examining the spatial correlation of the multifunctional coupling coordination degree of land use in the YRDUA, we found that it has a certain spatial agglomeration effect, exhibiting a trinuclear structure, with Shanghai, Nanjing, and Hangzhou at the apex. This pattern is different from that of the other two economic agglomerations in China: Beijing–Tianjin–Hebei and the Pearl River Delta. In a study of PLE functions of land use in Beijing and Tianjin, Yu et al. found that the coupling coordination degree presents an outward decreasing dual-core structure centered on Beijing and Shijiazhuang [47]. Moreover, Han et al. [54], in a study of the coupling and coordinated development of “economy–society–ecology–population-land” urbanization in the Pearl River Delta metropolitan area, found that the high-value areas are concentrated in the central and southern regions of the Pearl River Delta. The reason for the relatively fragmented distribution of the YRDUA high-value areas is mainly due to the fact that there are many fast-developing cities, such as Suzhou, Wuxi, Changzhou, Shaoxing, etc. These cities, which are scattered within the city cluster, are less dependent on the core cities and can effectively improve the coupling coordination level of LUFs with their own advantages. Additionally, Yu et al., in a study on the spatial correlation of land use carbon emissions in the Yangtze River Delta region from 1995 to 2018, found that Shanghai and Wuxi have been playing a dominant role in the spatial correlation network and have a large land use carbon radiation range [55]. The spatial evolution of the coupling coordination degree of LUFs in the YRDUA revealed that the coupling coordination among some cities in the downstream region of the Yangtze River basin, such as Xuancheng and Chizhou, is at a low level, which is consistent with the results of Lin et al. [56]. This could be attributed to the fact that these regions have replaced the relatively low-end industries in Jiangsu, Zhejiang, and Shanghai. Therefore, considering this, the urban agglomeration in the future should design a differentiated optimization and coordination path of urban PLE functions based on the spatial difference characteristics of the coupling coordination level and its evolution law. Simultaneously, the radiation-driving ability of regions with high values of coupling coordination should be impressed upon, and the support for low-value regions should be increased to better promote the integration and high-quality development of the Yangtze River Delta.

5. Conclusions and Policy Implications

The study of LUFs is an important proposition in land science as well as geography research. It has important theoretical and applied values for improving land use efficiency and promoting regional sustainable development. Taking the YRDUA as the research object, in this study, we constructed the evaluation index system from the perspective of “production–living–ecological” functions, and used the coupling coordination model to explore the coordination status of the PF, LF, and EF during 2005–2020. The kernel density estimation and trend surface analysis are used to further investigate the spatial and temporal variation of the coupling coordination. Finally, we also used the autocorrelation model to explore the spatial correlation characteristics of the coupling coordination degree of the YRDUA. The following primary conclusions were drawn:

(1)

Levels of LUFs in the YRDUA improved from 2005 to 2020, with the spatial distribution of each function level varying significantly. Among them, the high-value areas of the PF, LF, EF, and other LUFs are concentrated in the eastern, central, southwestern, and southern regions of the YRDUA, respectively. Overall, the EF level is low in the areas with high PF and LF levels.

(2)

During the study period, the coupling coordination degree of LUFs in the YRDUA showed an upward trend to varying degrees, and the spatial and temporal differentiation was evident. In particular, the coupling coordination degree of PF–LF progressed slowly, and the high-value areas are clustered around Shanghai, Nanjing, Wuxi, and their surrounding cities. Moreover, the PF–EF coupling coordination degree exhibited a fluctuation upward, with the spatial distribution showing a single-polarized development structure centered on Shanghai. The LF–EF coupling coordination degree rose slowly, with the high-value areas demonstrating a multi-core structure.

(3)

The coupling coordination degree of LUFs increased steadily, fluctuating in a few years, and the difference of coordination degree between cities increased. Overall, the change of coupling coordination degree in the north–south direction of the urban agglomeration was evidently higher than that in the east–west direction, showing an inverted “U” shape in the north–south direction and a pattern of “high in the east and low in the west” in the east–west direction. The development of marginal areas is relatively lagging, and the high-value areas exhibited a three-core structure, with Shanghai, Nanjing, and Hangzhou at the vertices.

(4)

The coupling coordination degree of LUFs in the YRDUA showed evident spatial positive correlation, with a certain agglomeration phenomenon in both high- and low-value areas.

Based on the above conclusions, this paper puts forward several policy implications:

(1)

Based on a full understanding of the interrelationship of the “PLE” functions, the relevant government departments should drive the optimization of the “PLE” functions, by taking the advantageous functions as the traction, and pay attention to the synergy with other functions in the process in order to improve the overall coordination of land functions in the region. In addition, government departments can further improve the level of coupling and coordination by combining the LUF coupling and coordination index with performance evaluations.

(2)

In developing the PF, the relevant departments should change the economic development mode to promote an intensive, efficient, and green low-carbon economy. This will realize the pulling effect of the PF to the LF. In addition, the proportion of the agricultural production function and non-agricultural production function in the process of land use should be coordinated to avoid the negative impact of excessive non-agricultural production on residents, which can also meet people’s demand for food. In turn, it can achieve the promotion of the PF to the LF. Finally, the relevant departments can gradually improve the living environment of residents by strengthening the construction of urban green parks and other methods. At the same time, the relevant departments can guide people to choose an environmentally friendly daily lifestyle, which can promote the coordination of the LF and EF.

(3)

In our research results, it can be proven that there is a radiating effect of a central city on the development of the surrounding cities in terms of the coupled coordination of LUFs. Then, by strengthening the position of the central city, the relevant departments can stimulate higher economic radiation and a trickle-down effect, which can drive the coordinated development of the surrounding towns. In addition, relevant departments should also develop and implement territorial spatial planning policies according to local conditions in order to reduce the disparity between regions.

6. Study Limitations

Although this study combines the coupling coordination degree of the YRDUA’s multiple LUFs as a whole and the coupling coordination degree between two functions for comparative analysis, it overcomes the research limitations of single measurement of the coupling coordination degree of LUFs and provides more realistic guidance for clarifying the optimization and coordination of the three functions of land use in urban agglomerations. However, there are still some limitations in this paper. Firstly, in the selection of indicators, due to the limitation of data availability, the function expression is not complete. Secondly, due to the serious lag of most county-level statistical data, this study only discusses the LUFs of the YRDUA with cities as the research unit. In the future, with the improvement of the statistical system, we will consider including county-level units into the sample to further enrich the research results. Finally, based on the analysis of the spatial and temporal evolution of LUFs and their coupling coordination degree, this paper still lacks an analysis of the influencing factors of their coupling coordination relationship, which is also the focus of our future research. However, the above problems do not affect the significance of the references cited in this paper.