Coupling Coordination Degree between Ecological Environment Quality and Urban Development in Chengdu–Chongqing Economic Circle Based on the Google Earth Engine Platform

Abstract

Rapid urbanization often exerts massive pressure on the resources relied upon by the ecological environment. It is necessary to quickly evaluate the interaction and mutual influence between regional urbanization and the ecological environment. This paper uses the Google Earth Engine (GEE) platform, integrates MODIS and night light remote sensing data sets, and computes the remote sensing-based ecological index (RSEI) and the coupling coordination degree (CCD) to measure the coupling coordination and analyze the spatiotemporal changes in the Chengdu–Chongqing Economic Circle (CCEC) for 2010, 2015, and 2020. Our results demonstrate four key findings. Firstly, the CCD varies spatially; it peaks at the Chengdu and the West Chongqing Plains, decreasing outwards along the mountains, with the lowest degree of coupling in the central, southern, and northern edge areas of the CCEC. Additionally, it has shown a trend of maintaining unchanged first and then increasing, mainly responding to policy decisions. Secondly, the changes between the different coupling levels were almost stable and mainly occurred between adjacent levels. Thirdly, the coupling level of towns spreads outwards from the centers at Chengdu and Chongqing and has an overall upward trend in time. Fourthly, in the most recent year, the coupling types present a distribution pattern of one developing axis connected with two peaks. Specifically, the environment system lagging type aggregates in Chengdu, Chongqing, and their surrounding areas, and the others mainly are economic system lagging type. The high internal coupling type also mainly occurs in the high and low coupling levels. Under this context, constructive suggestions for developmental optimization in the study area were proposed.

1. Introduction

Urban expansion is the most important human social change in the world [1]. Today, the main pattern of urbanization has transformed from the expansion of individual cities to the overall construction of urban agglomerations [2]. However, urbanized areas and living beings are facing huge ecological and environmental pressures due to polluted waste emissions from intense human activities [3]. Indeed, it is more than just a problem for developed countries but also for developing ones. For addressing these global issues, there are a lot of valuable works on the exploration and assessment of urban sustainability [4]. Moreno et al. believed that a “15-Minute City” is an efficient and advanced planning approach in the context of global crisis [5]. It is similar to the compact city, where the design applies compaction strategies. With the Internet of Things (IoT) technology gradually maturing, Belli et al. argued that smart, sustainable cities based on the IoT could improve urban management and government decision [6]. The coupling coordination is one of the critical aspects of sustainability, generally investigating the interactive cooperatingand dynamically evolving relationship among the social economy, eco-environment, and urban land use. Zuo et al. introduced the cloud model to an indicator system’s construction and evaluated, both quantitatively and qualitatively, the status of ecological civilization construction in China [7]. Ariken et al. established a complex framework for evaluating the CCD along the Silk Road Economic Belt from the demographic, economic, social and spatial aspects [8]. Among studies in this field, a coupling coordination degree model can intuitively explain the degree of interaction between two or more systems [9] and is widely used in research on modern urbanization, agricultural precision, and economic diversification [10,11,12]. However, in the previous study of coordination development with the eco-environment, few studies applied remote sensing technology, and most of the existing studies estimate city coupling by only using statistical data from the government or related institutions. While making the results more reliable, data collection and processing require expensive time and labor costs. Meanwhile, due to the differences in knowledge and target, existing methods are strongly affected by subjective factors. The achievements from various subject-based frameworks are also lacking comparability and generality. Fortunately, the emergence of remote sensing technology has dramatically improved the balance between cost and benefit, providing powerful support to experts in related disciplines.

In recent years, remote sensing technology has developed substantially. The advantages of its wide spatial coverage, rapid measurement times, and rich information that can be recorded mean that it has made urban ecological remote sensing studies more convenient. Compared with previous single indicators, such as NDVI, LAI, EVI, and so on, the current indicators are more comprehensive, complex, and systematic. Firozjaei et al. proposed a new Land Surface Ecological Status Composition Index (LSESCI) to distinguish the status of different land uses and covers [13]. Wu et al. proposed a novel remote sensing ecological vulnerability index (RSEVI) to assess ecological vulnerability [14]. The remote sensing-based ecological index (RSEI) model was proposed by Xu Hanqiu [15]. The model is based entirely on data collected from remote sensing technology, which is easy to obtain, objective, and reliable. The RSEI and its improvement have been widely used in current ecological quality assessments and analyses: Tang et al. evaluated the environmental degradation in typical mining regions based on the RSEI and found that the RSEI could better present the ecological quality than the landscape metrics [16]. Yuan et al. studied spatiotemporal change detection by using the RSEI in the Donting Lake Basin and analyzed the link between the RSEI and potential affecting factors [17]. Zheng et al. provided standard deviation-discretized RSEI to partly resolve the problem of instability in time series and cross-regional measurements [18].

Moreover, research has confirmed that night light data are highly correlated with urbanization, such as that by Andreano et al., who estimated the relation between lights and indicators of poverty. Their results demonstrated that using lights is helpful in poverty measurement and mapping [19]. Mellander et al. examined the nexus between night-time light and economic activities by taking residential and industrial samples from Sweden. They found that night-time light might be a better proxy for urbanization [20]. Wang et al. quantified the urbanization level on the Tibetan Plateau by using Luojia 1-01 (LJ1-01) night-time light data. Their results revealed that Lj1-01 has the potential for efficiently evaluating urbanization levels around the world, especially for less-developed regions [21]. So, the use of night light data to characterize urban development levels is feasible.

Furthermore, the rise of cloud-processing systems, for instance, Google Earth Engine (GEE), provides free-of-charge access to EO datasets worldwide. GEE, like other cloud computing platforms, is very popular because it provides efficient methods for storing, accessing, and analyzing datasets on high-performance servers. GEE was launched by Google in 2010 and made freely available remote sensing datasets via its internet-based Python Application Programming Interfaces (APIs) and a JavaScript web-based Interactive Development Environment (IDE) [22]. To that end, recent studies applied GEE to several environmental issues, such as erosion monitoring [23], wildfire mapping [24], land cover change [25], etc.

In 2021, China proposed a guideline for the construction of the CCEC, which emphasized enhancing the high quality of development, making it an important center of growth in the southwest region of China. Therefore, it is necessary to research the mode of evolution to ensure quality development.

With the aim of informing the economy, restricted with the eco-environment, and indicating the directions for the construction of an agglomeration, this paper proposes a framework entirely using remote sensing data to evaluate the coupling coordination in the CCEC. Specifically, this study uses the Google Earth Engine (GEE) platform to calculate the dryness, humidity, greenness, and heat indicators from Moderate Resolution Imaging Spectroradiometer (MODIS) remote sensing images acquired in 2010, 2015, and 2020. The RSEI is constructed using Principal Component Analysis (PCA). The coupling coordination degree model is used to include night light remote sensing data to quantify the degree of interaction between the urban economic system and the environment and to comprehensively evaluate the level of coupling between urbanization and the ecological environment for the region over the last 10 years. The contribution of this paper is a supplementation to the research gap of using remote sensing data to present urban coupling coordination and also offers a reference that is useful for all urban and rural development in the CCEC.

The structure is as follows: Section 2 introduces the study area and data sources. Section 3 shows the RSEI and CCD calculation and the Markov method. The spatiotemporal change analysis and division of the coupling coordination degree results are presented in Section 4. Section 5 discusses the results, limitations, and future research directions. Section 6 makes a concise conclusion based on the findings.

2. Study Area and Data Sources

2.1. Study Area

The Chengdu–Chongqing Economic Circle (27°39′~32°20′ N, 101°55′~109°02′ E) is located in Southwest China (Figure 1). It has a subtropical, humid monsoon climate with high temperatures and rain in the summer and mild and rainy conditions in the winter; it has an annual average temperature of 18 °C and an annual average rainfall of more than 1000 mm.

In the last decade, the population of the CCEC has increased from 110 to 116 million, with 492 people per square kilometer, accounting for about 8.20% of the national population. At the end of 2021, the CCEC achieved a GDP of 7391.92 billion yuan, more than two times that of ten years ago, accounting for 30.8% of the western region. The total area of the CCEC is 185,000 square kilometers, which has created 6.5% of the country’s GDP with 1.9% of the land area. In addition, the employment in each economic sector is 21.9:25.6:52.5, changed from the proportion of 40.2:25.7:34.1 ten years ago. Due to the richness of mineral resources, developed transportation, and talent gathering, there are four trillion-level industries in the CCEC: electronic information, automobile, equipment manufacturing, and consumer goods.

To ensure the integrity of the research area, our research scope included the Chengdu–Chongqing Economic Circle, Shizhu County, and Pengshui in Chongqing (Figure 1). For the same administrative level between the Sichuan and Chongqing cities, we used cities for the geographical research units for the Sichuan Province, and we used districts and counties as the units for Chongqing, making a total of 46 research units (hereinafter referred to as towns for convenience).

2.2. Data Sources

The data sources used in this paper are shown in Table 1. The RSEI includes four ecological elements: greenness, heat, dryness, and humidity. Its pre-processing is as follows:

The greenness is extracted from the Normalized Difference Vegetation Index (NDVI) band of the MOD13A1 V6 image set. The MOD13A1 data are a Level 3 product of MODIS, with the contents of the grid normalized vegetation index and enhanced vegetation index. Based on the 1B data, the product corrects the edge distortion generated by the remote sensor imaging process and provides a Level 3 sine projection product with a spatial resolution of 500 m that divides the whole year into 23 periods every 16 days. The heat component comes from the MOD11A2 V6 image set. The MOD11A2 is synthesized from the daily 1 km surface temperature/emissivity product of MOD11A1. It stores the average value under sunny weather in 8 days, and the projection is a sine curve projection.

The Normalized Difference Built-up and Soil Index (NDBSI) is calculated based on the 8-day periodic synthetic surface reflectance data provided by the MOD09A1 V6 dataset to represent the dryness. MOD09A1 provides an 8-day synthetic data product in bands 1–7 with a 500-m resolution, with the sine curve projection. Each pixel contains the most likely gridded level-2 (L2G) observation values in 8 days. This product has considered the effects of high observation coverage, the low angle of view, cloudless and cloudless shadows, and aerosol concentration.

The humidity component of the tassel cap transformation is also computed from the MOD09A1 V6 dataset based on transformation coefficients (also known as the Kauth–Thomas transformation, used to transform the spectral information of satellite data into spectral indicators) of MODIS. The MOD11A2 V6 and MOD13A1 V6 datasets have undergone an atmospheric correction to eliminate observation errors caused by water, clouds, aerosols, and cloud shadow. An effective mask fog and other interference factors are used to quality control the MOD09A1 V6 dataset to ensure the data’s integrity and image quality. The number of cloud-polluted pixels in the MODIS data is reduced through multi-temporal synthesis technology; however, we selected only the imagery from the summer (June, July, and August) to further minimize the risk of cloud contamination and, so, to enhance the reliability and comparability of the images in different time periods within our 2010–2020 study period.

Pre-processing includes applying a cloud mask, water mask, and re-projection of the data. Firstly, the quality assessment band (StateQA) of the MOD09A1 V6 dataset is used to obtain clear, cloud-free, and cirrus cloud-free high-quality remote sensing data for each image at the specified time. Secondly, the improved normalized water body index (MNDWI) is used to mask any water, to prevent the objectivity of the PCA from being affected by the water body [26]. The rivers in the Chengdu–Chongqing area are small, and the MODIS spatial resolution is medium, so the threshold value is set to 0. The 1000 m resolution MOD11A2 V6 dataset is re-projected to the SR-ORG:6974 projection coordinate system, and the nearest neighbor method is used to resample the resolution to 500 m. This ensures that the coordinates and resolution of all the data products are consistent.

A Prolonged Artificial Nighttime-light Dataset of China is a yearly synthetic data product. It is produced by the Night-Time Light convolutional Long Short-Term Memory (NTLSTM) network, which is a better socioeconomic indicator for build-up areas, population, and the GDP. Thus, this dataset is also re-projected to match the MODIS dataset and applied normalization in the cropped study area [27].

Table 1. Data description.

Sources	Dataset	Spatial Resolution	Description
Google Earth Engine	The National Aeronautics and Space Administration Digital Elevation Model (NASADEM)	30 m	Reprocessed Shuttle Radar Topography Mission (SRTM) data that offers improved accuracy over the original terrain data.
	MOD09A1 V6	500 m	Terra MODIS surface reflectance data.
	MOD11A2 V6	1000 m	8-day average synthetic surface temperature.
	MOD13A1 V6	500 m	The Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) vegetation index products.
A Big Earth Data Platform for Three Poles	A Prolonged Artificial Nighttime-Light Dataset of China (1984–2020) [28]	1000 m	Annual artificial night light data from China.

Table 1. Data description.

Sources	Dataset	Spatial Resolution	Description
Google Earth Engine	The National Aeronautics and Space Administration Digital Elevation Model (NASADEM)	30 m	Reprocessed Shuttle Radar Topography Mission (SRTM) data that offers improved accuracy over the original terrain data.
	MOD09A1 V6	500 m	Terra MODIS surface reflectance data.
	MOD11A2 V6	1000 m	8-day average synthetic surface temperature.
	MOD13A1 V6	500 m	The Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) vegetation index products.
A Big Earth Data Platform for Three Poles	A Prolonged Artificial Nighttime-Light Dataset of China (1984–2020) [28]	1000 m	Annual artificial night light data from China.

3. Methods

This study used the GEE platform and MODIS products to calculate the RSEI to represent the urban ecological environment quality and used the night light data to represent the urban development level. Using the coupling coordination model, the degree of coupling between the ecological environment quality and the economic development level over the past 10 years was calculated, and the spatial-temporal evolution and characteristics of the coupling are analyzed and discussed (Figure 2).

3.1. RSEI Construction

The NDVI has a linear and positive correlation with the vegetation distribution density. The NDVI is a good indicator of vegetation coverage and can, therefore, be used to quantify the regional vegetation growth status and to judge the regional ecological environment quality and degree of urban greening. The greenness component needed for this study was obtained from the NDVI data in the MOD13A1 V6 dataset.

Land Surface Temperature (LST) is a key parameter for studying the urban environment. The composition and structure of an urban surface are complex and spatially heterogeneous, including man-made structures, water surfaces, and green space. The “LST_Day_1km” band, provided as part of the MOD11A2 V6 dataset, were selected as the heat component for our study.

The tassel cap transformation compresses multi-band spectral information to extract a moisture component. This component is a helpful indicator of the soil moisture, to which it is closely related, and which it is, therefore, taken to represent. Humidity was calculated from the MOD09A1 V6 surface reflectance data, as shown in Equation (1):

$$WET=p_{1}Blue+p_{2}Green+p_{3}Red+p_{4}NIR-p_{5}SWIR1-p_{6}SWIR2$$

(1)

In Equation (1), Blue, Green, Red, NIR, ${\mathrm{SWIR}}_1$, and ${\mathrm{SWIR}}_2$ represent the blue light, green light, red light, near-infrared, the first shortwave infrared, and the second shortwave infrared bands of the MODIS sensor; $p_i$ (i = 1, 2, 3, …, 6) represents the empirical coefficients of the corresponding wave bands when calculating the humidity components for the tassel cap transformation, which are 0.1147, 0.2489, 0.2408, 0.3132, 0.3122, 0.6416, and 0.5087, respectively [29].

Where the impervious surfaces of houses, roads, and parking lots have replaced natural ecological surfaces, a significant urban heat island effect has resulted [30]. The dryness component in this study used the MOD09A1 V6 surface reflectance data to calculate the drought index, as shown in Equations (2)–(4) [31]:

$$SI=\frac{\left(SWI{R}_1+Red\right)-(Blue+NIR)}{\left(SWI{R}_1+Red\right)+(Blue+NIR)}$$

(2)

$$IBI=\frac{\frac{2SWI{R}_2}{\left(SWI{R}_1+NIR\right)}-\left[\frac{NIR}{Red+NIR}+\frac{Green}{SWI{R}_1+Green}\right]}{\frac{2SWI{R}_2}{\left(SWI{R}_1+NIR\right)}+\left[\frac{NIR}{Red+NIR}+\frac{Green}{SWI{R}_1+Green}\right]}$$

(3)

$$NDBSI=\left(SI+IBI\right)/2$$

(4)

In Equations (2)–(4), Blue, Green, Red, NIR, ${\mathrm{SWIR}}_1$, and ${\mathrm{SWIR}}_2$ represent the blue light, green light, red light, near-infrared, the first shortwave infrared, and the second shortwave infrared bands of the MODIS sensor.

The above components are calculated based on MODIS remote sensing data and remote sensing products. The median synthesis method, shown in Equation (5), was used to effectively fill and replace areas of the image that were affected by cloud and fog interference. Under this method, the affected areas were replaced with the median summer values for each component.

$$p=median\left[{\left\{p_i| p_i\ne null\right\}}_{i=1}^n\right]$$

(5)

In Equation (5), n represents the pixel value of the image that meets the filtering conditions, and p is the final pixel value at the same location after median synthesis.

The range standardization method was used to normalize each indicator so as to eliminate the dimensional impact between different indicators (considering unit inconsistencies between the different components and indicators). PCA was used to integrate all indicator information, reducing the redundancy between them. The cumulative contribution rate of the filtering components was more than 75%. Finally, the contribution rates of the different areas were taken as their weights, and the RSEI was calculated and normalized. The specific formula is described in Equations (6) and (7):

$$RSE{I}_0=PCA\left[f\left(NDBSI,WET,LST,NDVI\right)\right]$$

(6)

$$RSEI=\frac{RSE{I}_0-RSE{I}_{0\_min}}{\left(RSE{I}_{0\_max}-RSE{I}_{0\_min}\right)}$$

(7)

The RSEI is a value between 0 and 1. A higher value is interpreted as better regional ecological quality.

3.2. Coupling Coordination Model

We used the concept of capacity coupling to identify the degree of coupling for different towns in different years. This allowed us to build a capacity coupling coefficient model for ecological quality and level of urban development [32]. This is expressed in Equations (8)–(10):

$$C={\left[\frac{f\dot{\left(U\right)}\cdot g\left(E\right)}{{\left(f\left(U\right)+g\left(E\right)\right)}^2}\right]}^{\frac{1}{2}}$$

(8)

$$T=\alpha f\left(U\right)+\beta g\left(E\right)$$

(9)

$$D={\left(C\cdot T\right)}^{\frac{1}{2}}$$

(10)

In Equations (8)–(10), $\mathrm{f}\left(\mathrm{U}\right)$ is the indicator for the level of urban development, $\mathrm{g}\left(\mathrm{E}\right)$ is the indicator for urban ecological environment quality, and C indicates the degree of coupling between urban development and the ecological environment quality. A high value for C demonstrates a robust agreement between the internal elements of the system and means that the system will tend to a new orderly structure. T represents the coordination index for the urban development level (α) and ecological environment (β). The sub-system weight was set to 0.5 to indicate the same importance, and D represents the overall coupling coordination level of the system.

3.3. Markov Chain

A Markov chain is a type of Markov process that can be discretized in time and space. It is often used to represent processes of mutual transfer and change between different states in time series analysis. The transfer probability of each state is related only to its previous state. In this paper, a Markov transfer matrix was used to describe the structural transformation process for each of the research units. The evolutionary pattern was then analyzed, and reasonable countermeasures and suggestions were made. The calculation is given by Equation (11) [33]:

$$P_{ij}\left(E_i\to E_j\right)=\frac{n_{ij}}{n_i},$$

(11)

where $E_i$ and $E_j$ represent the i-th and j-th coupling coordination level, and $n_{ij}$ and $n_i$ represent the number of pixels transferring from the i-th to the j-th and i-th levels.

4. Results

4.1. Spatial Distribution Pattern Analysis

We calculated the CCD for 2010, 2015, and 2020 and obtained the spatial-temporal distribution pattern in Figure 3.

We divided the CCD into five equally-spaced levels, building on previous research [34,35] and our results. The levels are serious imbalance (I) ($\in =$ 0–0.2), moderate imbalance (II) ($\in =$ 0.2–0.4), near imbalance (III) ($\in =$ 0.4–0.6), moderate coordination (IV) ($\in =$ 0.6–0.8), and high coordination (V) ($\in =$ 0.8–1.0). The average degree of the coupling coordination for each town was then calculated (Figure 4).

Figure 3 shows that the CCD presents a diagonal layout with two centers at the Chengdu and Chongqing Plains. As the distance from the center increases, the size and quality of the coupling continually reduce. The coverage of levels 2 and 3 gradually expands outward over time. However, the Longquan Mountains and the eastern Sichuan parallel ridge-and-valley (Figure 1) are strong barriers, interfering with the connection and expansion between the cores and other towns.

Specifically, the central region and the southern and northern fringe regions of the study area have the lowest CCD; that is, level I. The areas of level V are concentrated in each central area of towns. The areas with a level IV or III often surround or are close to the centers of towns, while the areas of level II are mainly distributed in the suburban areas that radiate from the centers of towns. To sum up, the regions at levels V, IV, and III have prominent agglomeration characteristics in the space, with significant regional differences.

At the town scale (Figure 4), the coupling coordination is highest in the central urban area of Chongqing, followed by Chengdu and the towns around the central urban area of Chongqing. As the distance between Chengdu and the central urban area of Chongqing increases, the level of the coupling coordination gradually decreases. From 2010 to 2020, the level of the coupling and coordination of the towns obviously rose, primarily changing from level I to II and from level II to III. The changes in level mainly occurred in the central and border areas of Chengdu and Chongqing. By 2020, the towns of Luzhou, Fengdu, Shizhu, and Pengshui still maintained a level I.

4.2. Level Change Analysis

The descriptive statistical data in the study period were calculated (

Table 2. Spatial-temporal evolution statistics of ecological quality, night light brightness, and degree of coupling coordination.

Indicator/Level of Coupling Coordination	2010	2015	2020
I (%)	55.74	48.86	50.02
II (%)	15.50	20.87	15.92
III (%)	25.04	25.16	28.12
IV (%)	2.78	4.60	4.96
V (%)	0.93	0.50	0.97
Average RSEI	0.74	0.62	0.65
Average night light	305.38	426.86	439.24
Average coupling coordination	0.26	0.26	0.28

).

The results in Figure 3 and Table 2 demonstrate that the average coupling coordination showed a temporal trend of being stable first and then increasing ($0.26\to 0.26\to 0.28$). In particular, the central urban areas of Chengdu and Chongqing changed from level V in 2010 to level IV and were surrounded by areas at level V in 2015. Subsequently, the total area with coordination coupling at levels IV and V increased slightly in 2020. Meanwhile, the average RSEI has a trend of decreasing by 16.2% first and then increasing by 4.8% when the night light intensity kept rising, with a total increase of 43.8%.

From 2010 to 2015, the top two net conversion rates were levels I and II (6.88% and 5.37%, respectively), while becoming levels II and III from 2015 to 2020 (4.95% and 2.96%, respectively). There is clearly an inflection point in 2015.

According to Equation (11) and the principle of division given in Section 4.1, the transformation matrix of the level of CCD was calculated (

Table 3. Markov state transition probability matrices for coupling coordination levels.

	Ratio	I	II	III	IV	V
I	52.30	0.820	0.151	0.029	0.000	0.000
II	18.19	0.331	0.417	0.251	0.001	0.000
III	25.10	0.021	0.115	0.802	0.060	0.001
IV	3.69	0.000	0.000	0.123	0.776	0.101
V	0.72	0.000	0.000	0.005	0.524	0.471

The bold represents the highest transition probability of each level.

).

From 2010 to 2020, the coordination level transition for the whole study area presented a double-peak structure in levels I and II, with area ratios of 52.3% and 25.1%, respectively. The proportion of the level IV and V areas was less than 5%.

The evolution of the coordination coupling includes the following features:

(1)

The elements on the main diagonal of Table 3 are larger than the off-diagonal elements, indicating that the coupling coordination level has stability and dependence;

(2)

The off-diagonal elements in Table 3 decrease rapidly with distance from the main diagonal, indicating that the level changes for the coupling coordination are mainly active in the adjacent levels, and there is limited skipping between the levels. The weighted probability of a rise between the adjacent levels near the main diagonal is 0.143, which is higher than the probability of a level drop, which is 0.097;

(3)

Level rise mainly occurs between levels I and II, and the probability of transition from II to III is the highest, with a probability of 0.251. Decreases between the levels mainly occurred under two marginal conditions, namely, levels II and V, reaching the values of 0.331 and 0.524, respectively, in Table 3. According to Figure 3, this indicates that rapid development of the urban economy was concentrated in undeveloped regions and in highly developed regions, which caused the varying change of the levels in these regions.

4.3. Coupling Coordination Type Division

Figure 5 shows that in 2010, about 41% of towns had outliers, but the number of outliers and towns with outliers decreased significantly over time, and only 35% of towns had a small number of outliers in 2020. In particular, Pengshui, Yubei, Bishan, and Tongliang had both upper and lower abnormal points at the same time in the study periods, which means that their internal coupling coordination degree varied greatly. In addition, in 2010, about 74% of towns had a median coupling degree lower than the average value, which decreased to 57% in 2020, indicating that the coupling degree of all the towns showed an upward trend. According to the average value of all the towns, Suining and Tongnan exceeded the average value during the study period, while Guang’an, Liangping, and Jiangjin witnessed a significant decline.

To provide compatible development strategies for different towns, this paper uses the classification principles below and takes the CCD in 2020 as an example.

First, we used a natural breakpoint method to divide the CCD into five first-level classes and then into three second-level classes, according to the difference between the status of two sub-systems (economic and environmental). We then divide these into two third-level classes according to the relative heights of the mean and median lines in the box chart. If the median line is higher than the mean line, it means that more than half of the areas are higher than the average level, and the internal level of the CCD for the town is high. In contrast, it means the low internal level of the CCD. We divided the 46 towns into the 11 types shown in

Table 4. Classifications used for coupling types in 2020.

Code (Town Number)	Code (Town Number)
V-ENV-H (1)	III-B-L (1)
IV-ENV-H (2)	II-B-L (3)
IV-ENV-L (2)	II-B-H (9)
III-ENV-L (1)	II-ECO-L (15)
III-B-H (6)	II-ECO-H (3)
	I-ECO-L (3)

ENV: environment system lagging, B: systems balance, ECO: economic system lagging, H: high internal coupling coordination degree for the whole town, L: low internal coupling coordination degree for the whole town.

.

Table 4 shows that the proportion of towns with balanced sub-systems in 2020 reached 41.3%. For 58.7% of the towns, one system was lagging in development during the process of urbanization. The type of II-ECO-L has the largest number, accounting for about 33% of the total. They also cover a wide area (58% of the total area) and belong to the areas in urgent need of further planning in agglomeration. Therefore, more attention should be paid to the construction of these areas.

Figure 6 shows that the highest and lowest coupling types are almost all distributed in the area of Chongqing, including the types of V-ENV-H, IV-ENV-H, I-ECO-L, etc. It shows that there is an imbalance in the development between the towns of Chongqing, which highlights the importance of further optimization. Moreover, Chengdu and Zigong are the only remaining sub-high coupling coordination types, which are the III-ENV-L and II-B-L types, respectively. The II-B-H type is mainly located in the radiation overlapping area of Chengdu and Chongqing. Overall, the coupling types presented a spatial pattern of two cores and one development axis. The core area is mainly of the environmental lag type, and the area of the axis linking the cores is the status of the systematic balance. The economic lag type is mainly dominant in the edge regions of the study area. The distribution of the H and L types does not depend on the coupling level and sub-system status and is generally distributed in each coupling type.

5. Discussion

The CCEC is a vital part of the Yangtze River Economic Belt, and there are apparent differences between the economic development and environmental governance capabilities here and in those in the middle and lower urban agglomeration [36]. It is located inland and has a typical basin and hilly terrain. The efforts to reform and open up are relatively weak here, and the transformation and upgrading of the domestic industrial structure are relatively difficult [37]. Thus, how to solve the problem of studying the coordinated development of regional integration, ensuring the local ecology, and accelerating the development speed is a challenging problem.

This study found that the overall CCD of the CCEC continued to rise, but 2015 was a turning point due to policy factors, which is similar to Wan et al.’s conclusion [38]. China issued guiding policies on the construction of an urban agglomeration in the Chengdu–Chongqing region in 2011, 2016, and 2020, respectively, resulting in rapid socio-economic development and lagging development of the ecological environment in the early stage (2010–2015). In the subsequent development stage (2015–2020), the government realized the importance of ecological civilization construction, strengthened the protection of the ecological environment, and advocated the construction of green and sustainable development cities; hence, the overall coupling level in the follow-up rebounded.

In addition, from the perspective of the whole research period, the coupling level fluctuations mainly occurred in the lower and higher coupling levels. They are exactly the regions that need to be developed and the highly developed regions. The former is mainly due to the relatively backward economic system, while the latter is mainly due to the relatively backward ecological environment. The government should formulate environment-friendly development strategies, make rational use of its own resources, and maintain the mutual coordination of the urban system structure and functions, realizing the comprehensive and balanced development of the city.

Finally, this paper used the proposed analytical framework to classify the coupling coordination degree in 2020. The results show that this is basically consistent with the governmental planning concept and layout of the CCEC [39]. The government should actively play the driving role in Chengdu and Chongqing and, first of all, drive the coordinated development of adjacent areas in the space. In particular, in Chongqing as a whole, its internal coupling types are pretty different, which also showed in the study of Yang et al. [40], and the problem of inharmonious development in various regions is evident. The government should make an overall plan for these regions. At the same time, the mountains blocking the exchange and integration of the two cores and the towns affected simultaneously by the cores are relatively low couplings. The government needs to improve the level of urban development on the development axis, further strengthen the relationship between the two cores, expand the influence of the entire urban agglomeration, and radiate and drive the backward towns on the edge.

Although the conclusions of this paper are supported by other empirical studies, research on the CCEC is scarce. The applicability of the methods proposed in this paper still needs further verification, whether on multiple time scales or spatial scales. In addition, the type of remote sensing data and the number of indicators used in this paper are relatively limited. For example, the AOD elements can be added to characterize the air quality [41], and the POI data are good factors for describing the urban spatial features and functional structure [42]. At the same time, the interaction between the multiple systems is complex and nonlinear, which needs to study the mechanism of coupling [43]. The factors of land use, energy consumption, and social equity also affect the sustainable development of the cities [44]. The relevance of the proposed indicator in this paper and the above variables needs to be further explored. Therefore, in future research, it is necessary to improve and optimize the effectiveness, richness, and rationality of the indicators and further study and analyze the driving mechanism and coupling mechanism behind the coordinated development of the cities.

6. Conclusions

Based on remote sensing data, this study quantitatively analyzed spatiotemporal patterns of coupling between urbanization and the ecological environment in the CCEC during the last decade. In addition, we have given constructive suggestions for the coupling types in 2020. The main conclusions of this paper are as follows: (1) No matter the pixel scale or the urban scale, the CCD has a significant spatial aggregation, and its spatial distribution presents a double-core ring structure. The farther away from Chengdu and Chongqing, the more pronounced the decline of the coupling level and coupling size. Throughout the past decade, the poles’ radiation effect has had a positive impact on the surrounding areas, and the overall CCD has shown an upward trend. The rising regions include Nanchong, Tongnan, Ziyang, Yaan, and Yibin. In combination with the transition matrix, the level change mainly occurred in the coupling levels I, II, and V, which are the areas that need to be focused on in the process of urban development. (2) According to the distribution characteristics of the coupling types in 2020, the areas with higher coupling levels were the Chengdu and central urban areas of Chongqing, which are dominated by the type of environmental lag and need to take into account their own environmental protection while playing a leading role. The region with a medium coupling level was the corridor area connecting Chengdu and Chongqing, which is mainly based on balanced development. The low coupling level region was mainly distributed in the central and marginal regions with an economic lag. Internal coupling reflects the degree of harmonious development within the region whilst without an obvious distribution trend. In particular, the spatial differentiation of Chongqing is significant, and this region should be improved in combination with the local advantages and the current coupling state. The methods and conclusions of this paper contribute to the existing literature and help enhance the evolution of the CCEC towards coordination and harmony.