Spatiotemporal Evolution and Coupling Pattern Analysis of Urbanization and Ecological Environmental Quality of the Chinese Loess Plateau

Abstract

Understanding the interactive coupling mechanism between urbanization and eco-environmental quality is crucial to achieve the goal of urban sustainable development. The Chinese Loess Plateau (CLP) was taken as the research object, and the city nighttime light index (CNLI) and remote sensing ecological index with local adaptability (LARSEI) were constructed based on the data from the Defense Meteorological Satellite Program’s Operational Linescan System (DMSP/OLS), National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP/VIIRS), and Moderate Resolution Imaging Spectroradiometer (MODIS). Then, trend analysis, standard deviation ellipse (SDE), coupling degree (C), and coupling coordination degree (CCD) models were used to determine the spatiotemporal variation of urbanization and eco-environmental quality and its coupling relationship. The results show that: (1) the urbanization level of the CLP showed a trend of continuous improvement from 2000 to 2019. A significant increasing trend was found from the CNLI (slope_CNLI = 0.0030 yr⁻¹, p < 0.01), and its value rose from 0.07 in 2000 to 0.14 in 2019. In terms of spatial distribution, a multi-core distribution pattern with provincial capital cities as the core was presented in the CLP. The cities expanded at different degrees and presented a gradual concentrated expansion towards the southeast on the whole. (2) The eco-environmental quality in the CLP greatly increased during 2000 to 2019. An area with an increasing trend in the remote sensing ecological index with local adaptability (LARSEI) accounted for 58.82% and was mainly concentrated in the west and central part of the CLP. (3) The C and CCD between urbanization and eco-environmental quality in the CLP presented a trend of significant increase during 2000 to 2019 (slope_C = 0.0051 yr⁻¹, p < 0.01; slope_CCD = 0.0040 yr⁻¹, p < 0.01). The cities with a higher coupling degree were mainly located in the southeastern and northern parts of the CLP, while those with a higher coordination degree were scattered in the marginal parts of the CLP. The research results can provide suggestions for decision-making to achieve high-quality coordinated development of the cities in the CLP.

1. Introduction

The process of urbanization expresses the process in the type of society, from a traditional rural model, dominated by agriculture, to a modern urban society, characterized by the prominence of non-agricultural sectors, such as industry and services [1,2]. In parallel, with the rapid development of industrialization and the continuous growth of the national economy, the urbanization rate of China has increased from 36.09% in 2000 to 60.60% in 2019 [3,4,5]. Rapid urbanization promotes social and economic development; although it brings a variety of ecological problems, such as resource consumption, water shortages, and environmental pollution [6,7]. Therefore, it is of crucial importance to identify the spatiotemporal evolution and coupling relationship between urbanization and eco-environmental quality.

Several scholars have studied the spatiotemporal evolution of urbanization and eco-environmental quality, and their relationship. For example, Wang et al. [8] structured a comprehensive evaluation system of urbanization and eco-environmental quality in the Beijing-Tianjin-Hebei region based on statistical data, and further affirmed the existence of an interactive stress effect. Feng and Li [9] used statistical data to analyze the whole process of interaction of urbanization and eco-environmental quality of the Qinghai–Tibet Plateau. Gong et al. [10] analyzed the coordination relationship between urbanization and eco-environmental quality in the Guizhou Province based on statistical data and by adopting the entropy weight method and the coupling effect model. The abovementioned studies mostly used statistical data to evaluate the spatiotemporal evolution of urbanization and eco-environmental quality at different scales (i.e., province, city, and county). However, it is considerably difficult to analyze this phenomenon at finer scales (such as at pixel scale) [10,11,12].

Thanks to the development of space technology, some scholars constructed normalized difference building index (NDBI) and soil index (SI), or normalized difference vegetation index (NDVI), fractional vegetation coverage (FVC), enhanced vegetation index (EVI), ecological index (EI) and remote sensing ecological index (RSEI) to analyze the spatiotemporal evolution of urbanization and eco-environmental quality [13,14,15,16,17]. The proposal of comprehensive ecological assessment indices, such as the EI and RSEI, allowed us to offset the difficulty to analyze changes in the eco-environmental state using statistical data [17,18]. Not only this, but it also allowed us to fill a current gap: the difficulty to perform a multi-dimensional integrated inversion of eco-environmental quality. Although the EI index has the advantages of comprehensiveness and integrity, its constituent indices are difficult to obtain for several regions, which limits its further application [19,20,21]. The RSEI is constructed through principal component analysis (PCA) on the basis of greenness, dryness, humidity, and heat. This index has been widely used in the evaluation of regional eco-environmental quality, thanks to its relative comprehensiveness and effectiveness [22,23,24,25]. For example, Zhang et al. [26] used the Google Earth Engine (GEE) to construct an RSEI index, and evaluated the seasonal variations in eco-environmental quality in the Yangtze River Basin. Nie et al. [27] constructed an RSEI and evaluated the evolution, clustering situation, and influencing factors of eco-environmental quality in coal mining areas. In some studies, the RSEI considered the whole research area as a research object; thus making the evaluation results not as significant by performing an equalization treatment of the different environmental variables. Therefore, it is necessary to set an appropriate sliding window to calculate the weight of each indicator separately in order to fulfill the requirement of the RSEI for spatial adaptability [19,28,29].

In addition, some scholars used nighttime light data, such as DMSP/OLS and NPP/VIIRS, to analyze the spatiotemporal evolution of urbanization. Such data can not only better reflect the level of human activities and social economic development, but also compensate for the lack of statistical panel data. Consequently, they are extensively applied in urbanization research [30,31,32]. As for research on the relationship between urbanization and eco-environmental quality, the environmental Kuznets curve (EKC) [33,34], the coupled human and natural cube (CHNC) [35], the pressure–state–response (PSR) model [36], the CCD model [10], the local and tele-coupling theory [37], and the urbanization and eco-environment coupling coil theory [38] have commonly been applied. Specifically, the CCD and EKC are the two most commonly used models. However, in relation to developing countries, the results obtained by using the EKC are not optimal, possibly because of the independent consideration of the urbanization system and the eco-environmental system [39]. Moreover, the CCD is more suitable for China, because it focuses on the harmonious, symbiotic relationship between multiple systems [40].

The Chinese Loess Plateau (CLP), located in the middle reaches of the Yellow River, is the main area for national ecological security and for the implementation of high-quality development in the whole basin. Since the implementation of the western development strategy, urbanization of this region increased rapidly, and the five urban agglomerations of Guanzhong, Jinzhong, Lanxi, Hubao eyu, and Ningxia Coastal Yellow River have developed. Compared to the more developed eastern regions, urbanization in this region is still in a state of rapid development, and the overall urbanization level is relatively low. Moreover, the eco-environmental damage caused by rapid urbanization is particularly serious. Although some researchers have analyzed urbanization and eco-environmental effects in the CLP based on nighttime light data and statistical data, there is currently a lack of analysis of the spatiotemporal evolution of urbanization and eco-environmental quality, and of their coupling relationship at city and pixel levels, from multi-scale perspectives.

Thus, taking the CLP as an example, this study focused on the following: (1) it analyzed the temporal and spatial variations of the urbanization trend in the CLP from 2000 to 2019, based on nighttime light data, by constructing the CNLI index; (2) it evaluated the temporal and spatial evolution situation of eco-environmental quality in the CLP, based on the NDVI, normalized difference building soil index (NDBSI), WET, and land surface temperature (LST), by using the PCA model to calculate the remote sensing ecological index with local adaptability (LARSEI) index; and (3) it revealed the coupling coordination relationship between urbanization and eco-environmental quality in terms of the CCD model.

The contributions of this study are as follows: first, an evaluation index of the urbanization level from the perspective of nighttime light was constructed to analyze the development trend of urbanization level. Second, the remote sensing ecological index with local adaptability (LARSEI) index was constructed to analyze the changes of ecological environment quality. Finally, the coupling characteristics between the urbanization level and eco-environmental quality were analyzed across the urban scale and whole CLP. The results of this study provide some theoretical bases and scientific guidance for regional high-quality development.

The rest of this article is organized as follows. Section 2 mainly describes the general situation of the study area and the data sources used in the study. Section 3 introduces the data preprocessing, the main methods used in the research, and the process of index system construction. Section 4 focuses on the analysis of the results of this study, including three parts. The first part mainly analyzes the spatiotemporal variation trend of urbanization in the CLP, the second part focuses on the spatiotemporal evolution trend of eco-environmental quality in the CLP, and the coupling coordination between urbanization and eco-environmental quality is evaluated in the third part. Section 5 mainly discusses the applicability of the CLP of the indicators, the reasons for the above results, and suggestions. Finally, in Section 6, we summarize the main conclusions of this article.

2. Study Area and Materials

2.1. Study Area

The CLP, located in the central and northern parts of China (33°39′~41°17′ N and 100°52′~114°31′ E), is one of the four major plateaus in China, with an area of 6.35 × 10⁵ km² and an average elevation of 1407 m (Figure 1a). As an important bridge for the economic development of the eastern and western regions of China, the CLP borders Wushaoling in the west, the Taihang Mountain in the east, the Great Wall in the north, and Qinling in the south. It spans up to 1000 km from east to west and 850 km from north to south. The CLP includes all parts of the Ningxia Autonomous Region and Shanxi Province and parts of the Qinghai Province, Gansu Province, Shaanxi Province, Inner Mongolia Autonomous Region, and Henan Province (Figure 1c). The CLP suffers from a variety of human and geographical factors, including complex and diverse landcover types (Figure 1b), fragile eco-environment, significant differences in social and economic structure, and the rapid development of urban agglomerations, such as Guanzhong and Jinzhong. The ecological system of the CLP is at great risk [41].

2.2. Materials

The data product of nighttime light is an important tool to analyze human activities including the urbanization process [42]. The DMSP/OLS annual stable nighttime light data products, whose time span is from 1992 to 2013, and the NPP/VIIRS average monthly nighttime light radiation data products from 2012 to 2019, are developed by the National Geophysical Data Center of the National Oceanic and Atmospheric Administration (NOAA/NGDC; https:// www.ngdc.noaa.gov/ngdc.html; accessed on 25 November 2021). The spatial resolution of the DMSP/OLS is 1 × 1 km, and the radiation resolution is 6 bits. The pixel value DN (digital number), which refers to the relative brightness value without on-board radiometric calibration, represents nighttime light intensity. Although the influence of short-term lights, such as lightning and fishing boat lights, is eliminated in the data, the different sequences of images are not comparable [43,44]. The spatial resolution of the NPP/VIIRS is about 750 × 750 m, and the radiation resolution is 12 bits. The pixel value DN is the absolute brightness value after radiometric calibration, and the image data of different time periods are comparable. However, the influence of short-term lights, such as firelight and fishing boat lights are not eliminated [43,44].

The statistical panel data come from the Statistical Yearbook, the Statistical Yearbook of Chinese Cities, and the Statistical Bulletin of each province. The time range of the data is from 2001 to 2020, the temporal resolution is one year, and the spatial resolution is prefecture-level city.

The optical remote sensing data required for the construction of the LARSEI index include the MOD09A1, MOD13A3, and MOD11A2 data products, which come from the United States Geological Survey (USGS; https://www.usgs.gov/; accessed on 10 December 2021). The spatial resolution of MOD09A1 is 500 × 500 m, and the temporal resolution is 8 days. The spatial resolution of both MOD11A2 and MOD13A3 is 1 × 1 km, and their temporal resolution is 8 days and one month, respectively.

3. Methods

3.1. Construction of the Urbanization Index Based on Nighttime Light Data

3.1.1. Pre-Processing of Nighttime Light Data

Since nighttime light data were derived from different sensors, each with different parameters, they had to be jointly matched and corrected to obtain a long-term series [43,44]. The specific correction steps performed are as follows: (1) correction of outliers: because the nighttime light intensity is greater than zero, and the nighttime light intensity in Shanghai is the maximum, the threshold consists of two values. The outliers of all images were processed with this threshold method. (2) Spatial resolution correction: because the spatial resolution of the DMSP/OLS and NPP/VIIRS is about 1 × 1 km and 750 × 750 m, the NPP/VIIRS data were resampled to 1 × 1 km by upscaling, in order to solve the problem of inconsistency of resolutions between DMSP/OLS and NPP/VIIRS data. (3) Saturation correction: taking F16_2006 data, which are after performing saturation and radiometric correction, as the reference image and Jixi City of Heilongjiang Province as the reference area, the invariant target region method was used to perform exponential, logarithmic, and polynomial fitting of the annual DMSP/OLS images, so as to select the optimal function model and parameters (polynomial equations) as the correction function of DMSP/OLS images for saturation correction. (4) Intra-year fusion: due to the coexistence of multiple sensors in certain years, the mean value method was adopted to synthesize the images derived from different sensors. (5) Inter-annual correction: after performing saturation correction and intra-year fusion, as a problem of pixel incomparability between different years still existed, the inter-annual correction function (1) was used for correction. (6) Combined correction of DMSP/OLS and NPP/VIIRS: the invariant target region method was used to perform function fitting for 2013, and the optimal fitting function was selected as the correction function to correct NPP/VIIRS images. (7) Inter-annual recalibration: so as to ensure better flow and readability in different years after correction, inter-annual recalibration was used to obtain long-term series corrected nighttime light data.

$$D{N}_{(i,j)}=\{\begin{array}{ll}0& \hspace{1em}D{N}_{(i+1,j)}=0\\ D{N}_{(i-1,j)}& \hspace{1em}D{N}_{(i+1,j)}>0\&D{N}_{(i-1,j)}>D{N}_{(i,j)}\\ D{N}_{(i,j)}& \hspace{1em}other\end{array}$$

(1)

where DN_{(i, j)} represents the jth pixel in the image of the ith year; DN_{(i−1, j)} represents the jth pixel in the image of year (i−1)th; and DN_{(i+1, j)} represents the jth pixel in the image of year (i+1)th.

3.1.2. Urbanization Index of Nighttime Light

Based on the assumption that nighttime light intensity represents the social development process and the intensity of human activities, the CNLI was constructed to measure the degree of urbanization [45,46]. The calculation formula is as follows:

$$\mathrm{CNLI}=\frac{{\displaystyle \sum _{i=k}^{n}D{N}_i\times N_i}}{{\displaystyle \sum _{i=1}^{n}D{N}_{\max }\times N_i}}$$

(2)

where CNLI represents the city nighttime light intensity of the CLP; k represents the minimum DN value corresponding to the built-up area; DN_i represents the DN value corresponding to the ith gray level; N_i represents the number of pixels of the ith gray level; and DN_max represents the maximum of the gray level.

3.2. Construction of the Urbanization Indicators Based on Panel Data

3.2.1. Composition of the Index System

In order to verify the applicability of the CNLI to the CLP, this study further evaluated in a comprehensive way the urbanization process from four criteria layers, namely population urbanization, economic urbanization, spatial urbanization, and social urbanization [3]. Different methods of extreme value normalization were adopted to eliminate the differences in layer dimension in accordance with the attributes that influence the degree of urbanization development [47]. Then, the entropy and weight were calculated by adopting the entropy weight method [4]. The results are shown in

Table 1. Table of index system of urbanization.

Target Layer	Criterion Layer	Index Layer	Attribute	Entropy	Weight
Urbanization	Population Urbanization	Proportion of the urban population	+	0.952	0.043
		The density of the population	+	0.892	0.098
		Proportion of the non-agricultural employees	+	0.966	0.031
	Economic Urbanization	Proportion of the non-agricultural industry in output value	+	0.907	0.084
		Public budget revenue	+	0.885	0.104
		Per capita GDP	+	0.911	0.080
	Spatial Urbanization	Built-up area	+	0.953	0.042
		Per capita green area	+	0.869	0.118
		Per capita road area	+	0.890	0.099
	Social Urbanization	Number of beds in health care institutions per 10,000 people	+	0.874	0.114
		Total sales of goods	+	0.851	0.135
		Number of college students per 10,000 people	+	0.941	0.054

.

3.2.2. Urbanization Indicators of the Statistical Panel Data

The urbanization development index (UDI) was constructed on the basis of the four criteria layers mentioned in Section 3.2.1 to evaluate the applicability of the CNLI index to the CLP [48]. The calculation formula employed is as follows:

$${\mathrm{UDI}}_i={\displaystyle \sum _{k=1}^{12}{\omega }_k \times X_{(i,k)}}$$

(3)

where UDI_i is the urbanization development index in the ith year; ${\omega }_k$ represents the weight of the kth index; and $X_{(i,k)}$ represents the standardized value of the kth index in the ith year.

3.2.3. Standard Deviation Ellipse (SDE)

Clearly understanding the spatial morphological characteristics of each city is conducive to comprehend the development degree and direction of urbanization [49]. Therefore, the SDE was used to intuitively analyze the spatial expansion form of urbanization. The ellipse coverage represents the area of the cities with the highest urbanization, the larger the coverage area, the larger the spatial spread of urbanization, and the smaller the area, the higher the concentration of urbanization. The ellipse centroid was calculated to catch the urbanization change trend of the whole region, and the azimuth angle combined with the centroid position could be used to analyze the urbanization expansion direction of the whole region [50,51]. The specific calculation formulas are as follows:

$$\overline{X}=\frac{{\displaystyle \sum _{i=1}^n{\omega }_i{x}_i}}{{\displaystyle \sum _{i=1}^n{\omega }_i}},\overline{Y}=\frac{{\displaystyle \sum _{i=1}^n{\omega }_i{y}_i}}{{\displaystyle \sum _{i=1}^n{\omega }_i}}$$

(4)

$${\mathrm{SDE}}_x=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{(x_i-\overline{x})}^2}}{n}},{\mathrm{SDE}}_y=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{(y_i-\overline{y})}^2}}{n}}$$

(5)

$$\tan \theta =\frac{({\displaystyle \sum _{i=1}^n{\tilde{x}}_i{}^2-}{\displaystyle \sum _{i=1}^n{\tilde{y}}_i{}^2})+\sqrt{{({\displaystyle \sum _{i=1}^n{\tilde{x}}_i{}^2-}{\displaystyle \sum _{i=1}^n{\tilde{y}}_i{}^2})}^2+4({\displaystyle \sum _{i=1}^n{\tilde{x}}_i{\tilde{y}}_i{)}^2}}}{2{\displaystyle \sum _{i=1}^n{\tilde{x}}_i{\tilde{y}}_i}}$$

(6)

$${\sigma }_x=\sqrt{2\frac{{\displaystyle \sum _{i=1}^n{({\tilde{x}}_i\cos \theta -{\tilde{y}}_i\sin \theta )}^2}}{n}},{\sigma }_y=\sqrt{2\frac{{\displaystyle \sum _{i=1}^n{({\tilde{x}}_i\sin \theta +{\tilde{y}}_i\cos \theta )}^2}}{n}}$$

(7)

$${(\frac{x}{{\mathrm{SDE}}_x})}^2+{(\frac{y}{{\mathrm{SDE}}_y})}^2=S$$

(8)

where (x_i, y_i) are the geographic coordinates of the ith element; ${\omega }_i$ is the weight of the ith element; and $(\overline{X},\overline{Y})$ are the centroid coordinates of the SDE. $(\overline{x},\overline{y})$ are the arithmetic mean centers of (x, y), and $(\widetilde{x_i},\widetilde{y_i})$ are the difference between (x, y) and $(\overline{x},\overline{y})$. SDE_x and SDE_y indicate the long and short axes of the SDE, respectively. $\theta$ is the azimuth angle of the standard deviation ellipse, in which clockwise rotation in the due north direction is positive. ${\sigma }_x$ and ${\sigma }_y$ are the standard deviations of the x and y axes, respectively, and S is the standard deviation ellipse.

3.3. Construction of Remote Sensing Ecological Index with Local Adaptability (LARSEI)

Since greenness (NDVI) [24], humidity (WET) [51], heat (LST) [25], and dryness (NDBSI) [27] were required to construct the evaluation index of eco-environmental quality, in this study, the surface reflectance data of MOD09A1 and the finished product data, such as MOD13A3 and MOD11A2, were firstly mosaicked and clipped, the outliers were processed, and the annual values were synthesized. Meanwhile, in order to ensure the consistency of the constructed index resolution, WET and NDBSI were resampled to 1 × 1 km. Then, by applying the PCA, these data were further calculated to obtain the RSEI index. Since the RSEI index could not fulfill the requirement of local adaptability, the method developed by Zhu et al. [19] was referred to improve the RSEI. Different from RSEI, which took the whole image as a research unit, the remote sensing ecological index with local adaptability (LARSEI) took each sliding window as a research unit. All pixels in the sliding window were used to calculate the weight of each ecological index by the PCA method, and assigned it to the center pixel (i, j) of the sliding window (Figure 2). Therefore, PC1 was obtained by weighted summation of all ecological indexes and then the LARSEI images were obtained by the formula in

Table 2. Data sources and calculation formulas required for the construction and validation of LARSEI.

Indicators	Spatial Resolution	Time Resolution	Sources	Product or Formula	Functions
NDVI	1 km	Month	USGS	MOD13A3	To calculate RSEI and LARSEI
LST	1 km	8 days	USGS	MOD11A2	To calculate RSEI and LARSEI
Surface reflectance	500 m	8 days	USGS	MOD09A1	To calculate WET and NDBSI
WET	500 m	8 days	Surface reflectance	$\begin{array}{l}\mathrm{WET}=0.1509\times {\rho }_{Blue}+0.1973\times {\rho }_{Green}+0.3279\times {\rho }_{Red}\\ +0.3406\times {\rho }_{NIR}-0.7112\times {\rho }_{SWIR1}-0.4572\times {\rho }_{SWIR2}\end{array}$	To calculate RSEI and LARSEI
NDBSI	500 m	8 days	Surface reflectance	$\begin{array}{c}BSI=\frac{({\rho }_{SWIR1}+{\rho }_{Red})-({\rho }_{NIR}+{\rho }_{Blue})}{({\rho }_{SWIR1}+{\rho }_{Red})+({\rho }_{NIR}+{\rho }_{Blue})}\\ BI=\frac{\frac{2{\rho }_{SWIR1}}{{\rho }_{SWIR1}+{\rho }_{NIR}}-(\frac{{\rho }_{NIR}}{{\rho }_{NIR}+{\rho }_{Red}}+\frac{{\rho }_{Green}}{{\rho }_{Green}+{\rho }_{SWIR1}})}{\frac{2{\rho }_{SWIR1}}{{\rho }_{SWIR1}+{\rho }_{NIR}}+(\frac{{\rho }_{NIR}}{{\rho }_{NIR}+{\rho }_{Red}}+\frac{{\rho }_{Green}}{{\rho }_{Green}+{\rho }_{SWIR1}})}\\ \mathrm{NDBSI}=\frac{BSI+BI}{2}\end{array}$	To calculate RSEI and LARSEI
RSEI-2	1 km	Annual	Xu et al. [51]	$\begin{array}{l}PC1=PCA\left(\mathrm{NDNV},\mathrm{NDBSI},\mathrm{LST},\mathrm{WET},\mathrm{AI}\right)\\ \mathrm{RSEI}\text{-}2=\frac{PC1-PC{1}_{\min }}{PC1-PC{1}_{\max }}\end{array}$	To compare the applicability of the ecological environmental index of the CLP
RSEI	1 km	8 days	NDVI NDBSI LST WET	$\begin{array}{c}PC1=PCA(\mathrm{NDVI},\mathrm{MNDWI},\mathrm{LST},\mathrm{NDBSI})\\ \mathrm{RSE}{\mathrm{I}}_0=1-PC1\\ \mathrm{RSEI}=\frac{\mathrm{RSE}{\mathrm{I}}_0-\mathrm{RSE}{\mathrm{I}}_0{}_{\min }}{\mathrm{RSE}{\mathrm{I}}_{0\max }-\mathrm{RSE}{\mathrm{I}}_{0\min }}\end{array}$	To compare the applicability of the ecological environmental index of the CLP
LARSEI	1 km	Month	NDVI NDBSI LST WET	$\begin{array}{c}\mathrm{LARSE}{\mathrm{I}}_0=1-PC1\\ \mathrm{LARSEI}=\frac{\mathrm{LARSE}{\mathrm{I}}_0-\mathrm{LARSE}{\mathrm{I}}_0{}_{\min }}{\mathrm{LARSE}{\mathrm{I}}_{0\max }-\mathrm{LARSE}{\mathrm{I}}_{0\min }}\end{array}$	To assess the dynamic changes of eco-environmental quality

. In this way, the eco-environmental assessment index (LARSEI), whose temporal and spatial resolution are one month and 1 × 1 km, was constructed to invert the changes in eco-environmental quality in the CLP. The specific data sources and calculation formulas are shown in Table 2.

3.4. Coupling Degree (C) and Coupling Coordination Degree (CCD) Models

In the process of economic development, urbanization is bound to impact the eco-environment. The understanding of the relationship between them is crucial to grasp the goal of coordinated development of the human progress and eco-environmental protection [23]. Therefore, coupling degree (C) and coupling coordination degree (CCD) were adopted to analyze the coupling degree and coordination degree between the urbanization development level and eco-environmental quality in the CLP. By following Xu et al. [51], the C and CCD were classified as shown in

Table 3. The classified criterions of C and CCD.

Coupling Degree (C)		Coupling and Coordination Degree (CCD)
Type	Range of C	Type	Range of CCD
Primary coupling	(0.0, 0.2)	Extreme disorder	(0.0, 0.2)
Severe confrontation	(0.2, 0.3)	Severe disorder	(0.2, 0.3)
Moderate confrontation	(0.3, 0.4)	Moderate disorder	(0.3, 0.4)
Slight confrontation	(0.4, 0.5)	Slight disorder	(0.4, 0.5)
Primary adjustment	(0.5, 0.6)	Primary coordination	(0.5, 0.6)
Middle adjustment	(0.6, 0.7)	Middle coordination	(0.6, 0.7)
High adjustment	(0.7, 0.8)	High coordination	(0.7, 0.8)
Extreme coupling	(0.8, 1.0)	Extreme coordination	(0.8, 1.0)

. The specific formula employed is as follows:

$$\mathrm{C}={(\frac{\mathrm{U}\times \mathrm{E}}{{(\mathrm{U}+\mathrm{E})}^2})}^{\frac{1}{2}}$$

(9)

$$\mathrm{T}=\alpha \mathrm{U}+(1-\alpha )\mathrm{E}$$

(10)

$$\mathrm{CCD}=\sqrt{\mathrm{C}\times \mathrm{T}}$$

(11)

where U represents the urbanization development index; E represents the eco-environmental quality index; and T represents the comprehensive coordination index. Since urbanization and ecological environment are of common importance, $\alpha$ was set at 0.5 [51].

4. Results

4.1. Spatiotemporal Evolution Pattern of Urbanization

4.1.1. Temporal Changes in Urbanization

The CNLI value in the CLP increased significantly in the study period, from 0.07 in 2000 to 0.14 in 2019, with an increase rate of 0.0030 yr⁻¹ (p < 0.01), as shown in Figure 3a. This indicates a rapid increase of urbanization in the study region. The CNLI value followed a trend of rapid growth for all cities investigated. The CNLI values for XA, YIC, HD, and ORDOS showed the largest increase, reaching 0.0089 yr⁻¹ (p < 0.01), 0.0079 yr⁻¹ (p < 0.01), 0.005 yr⁻¹ (p < 0.01), and 0.0047 yr⁻¹ (p < 0.01), respectively. In parallel, the CNLI values for JC, LF, DT, and LL showed the smallest increase, equaling to 0.0012 yr⁻¹ (p < 0.01), 0.0009 yr⁻¹ (p < 0.01), 0.0006 yr⁻¹ (p < 0.05), and 0.0005 yr⁻¹, respectively (Figure 3b). During this process, the growth rate of the CNLI slowed down during 2003–2005 and 2007–2008, due to the impact of SARS and the global financial crisis, and the volatility characteristics were similar to Song et al. [52].

4.1.2. Intensity of Urban Expansion

Figure 4 shows the expansion trend of the core urban areas of the CLP from 2000 to 2019. It was found that urbanization in the CLP presented a multi-core distribution pattern, with each city expanding rapidly at different speeds. In the northern part of the CLP, the urban expansion degree of BT and HET reached 214.34% and 41.65%, respectively. The maximum occurred in the period of 2000–2003, reaching 104.09% and 16.60%, respectively. In the eastern part of the CLP, the urban expansion rate of TY was 92.09%, with the largest expansion rate (34.46%) being recorded from 2003 to 2007. In the southern part of the CLP, the cities of XA and LY showed an increase of 274.70% and 145.68%, respectively, and the largest increase reaching 46.14% and 28.93% were recorded in 2011–2015 and 2007–2011, respectively. In the western part of the CLP, LZ, YIC, and XN had an expansion rate of 135.73%, 296.98%, and 60.66%, with the largest growth rates of 54.99%, 62.50%, and 34.79%, respectively, during 2000–2003 and 2011–2015.

4.1.3. Direction of Urban Expansion

Figure 5 shows the deviation trajectory of the SDE centroid of urban expansion in the CLP from 2000 to 2019. During 2000–2019, the centroid azimuth angle of the CLP changed from 83.33° to 88.37°, showing an overall trend of concentrated expansion towards the southeast. From 2000 to 2010, the centroid azimuth angle of urban expansion in the CLP changed from 83.33° to 86.52°, and the cities of the area showed an overall trend of eastward expansion (Figure 5a,c). During 2010–2015, the centroid azimuth angle of the CLP changed little, and the direction of urban expansion began to turn back towards the southwest in combination with location elements. From 2015 to 2019, the centroid azimuth angle of the CLP changed from 86.27° to 88.37°, showing a trend of gradual expansion towards the southeast. At the same time, the coverage area of the SDE fluctuated from 4.62 × 10⁵ km² in 2000 to 4.40 × 10⁵ km² in 2019, with a decrease rate of 0.0107 km² yr⁻¹, showing a trend of gradual concentration to the southeast (Figure 5b).

4.2. Spatiotemporal Evolution Pattern of Eco-Environmental Quality

4.2.1. Temporal Changes in Eco-Environmental Quality

Figure 6 illustrates the temporal variation trend of LARSEI and its components in the CLP from 2000 to 2019. The LARSEI value increased from 0.31 in 2000 to 0.38 in 2019, with an increase rate of 0.0026 yr⁻¹ (p < 0.01), indicating a significant improvement in the eco-environmental quality of the CLP. In relation to the various factors considered, NDVI exhibited a trend of significant increase (p < 0.01), while WET represented a trend of significant decrease (p < 0.05), and both LST and NDBSI changed insignificantly.

4.2.2. Spatial Changes of Eco-Environmental Quality

For further understanding the changes in eco-environmental quality in different areas of the CLP, the spatiotemporal evolution trends of LARSEI were analyzed at both city and pixel scale (Figure 7). In the field of city scale, the eco-environmental quality in XN, XA, and BJ were better, with LARSEI reaching 0.51, 0.50, and 0.49, respectively. The eco-environmental quality of WH and ORDOS relatively worsened, with LARSEI reaching 0.20 and 0.22, respectively. In the majority of cities, such as TY, YQ, and CZ, LARSEI increased at first and then decreased following an inverted U-shaped trend (Figure 7a). As for the pixel scale, LARSEI increased in 89.07% of the CLP area, with 58.82% of the areas showing a significant increase trend (p < 0.05), mainly in the western and central regions of the CLP. The LARSEI showed a significant decrease trend in only 2.32% of the areas, mainly in the southeast of the CLP (Figure 7b,c). These results indicated that LARSEI had generally an upward trend in the CLP, indicating a significant improvement of the eco-environmental quality of the whole CLP. At the same time, it had a downward trend in the areas with rapid urban expansion. Therefore, more attention should be paid to the influence of urbanization on the ecological environment.

4.3. C and CCD Distribution Changes between Urbanization and Eco-Environmental Quality

Figure 8 shows the spatiotemporal changes of the C and CCD between urbanization and eco-environment in the CLP during 2000 to 2019. As can be seen from the Figure 8a, the C in the CLP changed from a stage of high adjustment to a stage of extreme coupling, and its value increased from 0.76 in 2000 to 0.83 in 2019, at a rate of 0.0051 yr⁻¹ (p < 0.01). The CCD in the CLP changed from a stage of moderate disorder to a stage of slight disorder, and its value increased from 0.36 in 2000 to 0.45 in 2019, at a rate of 0.0040 yr⁻¹ (p < 0.01) (Figure 8b). These indicate that coupling degree and coordination degree in the CLP constantly improved.

In terms of spatial distribution of C, Figure 8c described the time variation trend of C in different cities of the CLP. There were almost no cities in the stages of primary coupling and confrontation in the CLP in 2000. HD was the only city in the primary adjustment stage, and 13 cities were in the middle adjustment stage, mainly located in the southern part of the CLP. At the same time, 11 cities were in the high adjustment stage, mainly in the east, north, and southeast of the CLP, and 16 cities stayed in the extreme coupling stage, mainly in the north and west of the CLP. With the vigorous development of urbanization, the C between urbanization and eco-environmental quality gradually increased. By 2019, all cities in the primary and middle adjustment stages passed into the high adjustment stage or above. Among them, the cities in the extreme coupling stage were mainly located in the central part and north of the CLP (Figure 8d).

As for the CCD, Figure 8e described the time variation trend of CCD in different cities of the CLP, and there were no cities in the stage of extreme disorder in the CLP in 2000. There was only one city in the stage of severe disorder, namely ZW, located in the northwestern part of the CLP. In total, 24 cities, mainly located in the central part and west of the CLP, were in the moderate disorder stage, and 16 cities were in the slight disorder stage, mainly located in the eastern part of the CLP, among which XA, TY, and XN had the highest CCD. After the promulgation by the government of the strategy for the symbiotic development of environment and urbanization, the CCD began to increase, especially in the central part and west of the CLP. By 2019, the cities in the severe disorder stage passed into a moderate disorder stage or below, and there were only five cities in a moderate disorder stage, mainly located in the southwestern part of the CLP. The number of cities in the slight disorder stage was up to 24, mainly located in the central part and east of the CLP, and there were 11 cities in the stage of primary coordination, mainly in the southeast of the CLP. Meanwhile, it was worth noting that XA was the only city in the middle coordination stage (Figure 8f). The results show that all cities in the CLP belonged to the category of city with advanced eco-environmental development (Figure 3b and Figure 7a). Although the coupling coordination relationship between urbanization and eco-environment in the CLP gradually improved and the contradiction relationship between urbanization and ecological development gradually decreased, the relevant departments of the cities in the CLP should still take corresponding measures to achieve a coordinated coexistence between urban development and ecological environment.

5. Discussion

5.1. Applicability Analysis of CNLI to the CLP

Because of its capacity to clearly and directly describe the intensity of human behavior activities and the level of social and economic development, nighttime light data have been widely used to estimate the various social factors, such as the inversion of urbanization level [53], the simulation of population density [42], the estimation of urban land expansion [54], the evaluation of urban carbon emissions [55,56] and power consumption [57]. In terms of evaluation of urbanization level indicators, several studies have been conducted based on nighttime light data. For example, based on statistical data and by using NPP/VIIRS, Xie et al. [58] found a strong linear relationship between urbanization and nighttime light intensity in the process of inversion of the urbanization level of 64 cities in the United States. Wang et al. [49] used nighttime light data to construct an average nighttime light intensity index to evaluate the urbanization level of Beijing-Tianjin-Hebei at county and city levels, and found that the average nighttime light intensity had a strong linear relationship with the urbanization index constructed based on statistical data.

Accordingly, based on nighttime light data, in this study a long-term series of CNLI was constructed to evaluate the urbanization development level of the CLP. The results show that urbanization in the CLP followed a trend of significant increase (Figure 3a), which may be the result of the country’s efforts to develop urbanization with the advancement of industrialization. These results are similar to the findings of Song et al. [41] and Xu et al. [46]. As for urban expansion intensity, due to the national reform of the urbanization development system and the promulgation of new urbanization policies, the expansion degree of urbanization was different across the cities in the CLP, which is in line with the expectations of General Secretary Xi Jinping for the construction of new urbanization. In relation to the direction of urban expansion, the CLP showed a general trend of concentrated expansion toward the southeast, which may be the result of the rapid development of the eastern coastal cities. This result fits well with the findings of Wang et al. [49] in their study of the migration direction of the Beijing-Tianjin-Hebei urban agglomeration. At the same time, this result is different from the conclusions of Lu et al. [50], who concluded that urbanization tended to be scattered when calculating the coverage area by using the SDE, which may be caused by the fact that coastal cities were taken as the research object.

To further verify the applicability of CNLI to the CLP, the correlation between CNLI and UDI was analyzed based on statistical data. As can be seen from Figure 9a, the correlation between CNLI and UDI was high (R² = 0.90, p < 0.01). From the perspective of each city (Figure 9b), the correlation coefficient between CNLI and UDI was relatively better on the whole, with an average Pearson coefficient of 0.82. Only the cities of DT, LL, DX, and TS had lower correlation coefficients, mainly because they are resource-dominated cities and in recent years they experienced a process of transformation [59]. The correlation coefficient between CNLI and UDI in the other cities was higher, indicating that CNLI can better reflect the urbanization level in the CLP. This finding fits well with the results of Zhang et al. [60] in their study of analyzing the coupling state of urbanization and socio-economic development in the Yellow River Basin.

5.2. Applicability Analysis of LARSEI to the CLP

The rapid development of urbanization has generated serious damage to the ecological environment. The issue of how to quickly and accurately assess the state of ecological environment has become crucial for its restoration. In this study, LARSEI was constructed to evaluate changes in eco-environmental quality in the CLP based on greenness, dryness, heat, and humidity. The results indicate that in 58.82% of the area in the CLP the LARSEI showed a significant increase trend, mainly distributed in the central part and west of the CLP, primarily because these cities are the key areas of ecological engineering construction. The regions with a decrease in LARSEI value were mainly located in the Guanzhong and Jinzhong urban agglomerations in the southeastern part of the CLP (Figure 7b), mainly due to the rapid expansion of urbanization [14]. These results are in line with Sun et al. [61] and contrast with the findings of Xu et al. [51], mainly due to the addition of the abundance index (AI) in the RSEI and to diverse standardization methods.

To further verify the applicability of LARSEI to the CLP, RSEI, RSEI-2, and LARSEI were used for correlation analysis with NDVI, LST, WET, and NDBSI, respectively, to measure the applicability of eco-environmental assessment. The results, presented in the Figure 10a, show that RSEI, RSEI-2, and LARSEI had the highest correlation coefficients with NDVI and the lowest correlation coefficients with LST. The reasons may be caused by the increase of the transpiration and the decrease of the respiration under the higher surface temperatures [62]. These results did not fit well with the findings of Shan et al. [22], the reason being the difference between the study regions. For the single index correlation coefficient, RSEI had the highest correlation coefficient with NDVI and NDBSI (0.99 and 0.85, respectively), and LARSEI had the highest correlation coefficients with LST and WET (0.65 and 0.89, respectively). Owing to their own advantages of the indices, their applicability to the CLP cannot be significantly compared. Therefore, the $\overline{Cp}$ index [22] was introduced to comprehensively evaluate the relevance of LARSEI, RSEI, and RSEI-2 to these four indicators. As can be seen from Figure 10b, the C_p index of LARSEI was the highest (1.12), followed by that of RSEI (1.10), while the C_p index of RSEI-2 was the lowest (0.68). This confirms that LARSEI could be used for eco-environmental assessment of the CLP.

6. Conclusions

Based on DMSP/OLS, NPP/VIIRS, and MODIS data, the CNLI and LARSEI indices were constructed to analyze the development trends of urbanization and eco-environmental quality, and the coupling coordination degree model was combined to distinguish the development status of cities in the CLP. The results show that: urbanization in the CLP increased rapidly from 2000 to 2019. Urbanization in the cities of XA, ORDOS, HD, and YIC developed more rapidly. In terms of spatial distribution, a multi-core expansion pattern was present in the CLP, and the expansion intensity gradually decreased from the provincial capital cities to the periphery. Moreover, the cities are expanding at different rates, and a trend of concentrated expansion towards to the southeast was found from all of the cities in the whole CLP. Meanwhile, the eco-environmental quality of the CLP has significantly improved. LARSEI showed an upward trend, and the eco-environmental quality of the cities of XN, XA, and BJ was relatively better. As for the spatial distribution, LARSEI presented an upward trend in most parts of the CLP, and a significant increase in half of the area, mainly corresponding to the central part and west of the CLP. In a few parts of the CLP, LARSEI significantly decreased, mainly corresponding to the southeastern part of the CLP. It was found that the C and CCD between urbanization and eco-environmental quality of the CLP improved significantly progressively, respectively, by analyzing the relationship between urbanization and ecological environment. The cities with a higher C were mainly located in the north, central part, and west of the CLP, while the cities with a greater CCD were mainly clustered around provincial capitals.

However, this research also presents some shortcomings. Because the nighttime light data were derived from various sensors, the problem of how to match data from different sources is debatable. Moreover, because the resolution employed represents the degree to which geographic objects could be identified, the relatively low resolution of the nighttime light data may hinder the good reflection of urban information at a small scale. Thus, it is imperative to improve the resolution of the nighttime light data by combining them with high-resolution optical images. Moreover, the nighttime light data have a lot to do with consumer habits in different areas and government attitudes to energy use, but these factors have not been considered in this study. In addition, the difference of natural environment base is the key element to evaluate the quality of cities and eco-environment, and this should be taken into account in future studies. Furthermore, the impact of urbanization on other social systems and the possible solutions still remained to be studied. In future research, a profound analysis, including the coupling mechanism and influence factors between urbanization and other life systems, should be performed from each respect, in order to shape a balanced and sustainable development between humans and the ecological environment.