Using Remote Sensing Data to Study the Coupling Relationship between Urbanization and Eco-Environment Change: A Case Study in the Guangdong-Hong Kong-Macao Greater Bay Area

Promoting the coordinated development of urbanization and eco-environment is essential for the development of urban agglomerations. As an emerging economic growth pole in China, the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) has become a research hotspot in recent years. However, relevant studies in this area have been largely constrained by the incomparability of statistical data between the inland part of the GBA and the two Special Administrative Regions (Hong Kong and Macao) of the GBA. This study used nighttime light data and the normalized difference vegetation index (NDVI) to evaluate the level of urbanization and eco-environment in the GBA, respectively. Then we adopted the gravity center model to analyze the overall coupling situation between urbanization and eco-environment change. We also adopted a coordination index to determine the spatial differentiation characteristics of this coupling in the GBA from 2000 to 2018. The results show that (1) the spatial pattern of urbanization and eco-environment both show an approximately circular structure, and the change rates of these two variables show significant spatiotemporal differentiation; (2) on the whole, urban development and eco-environment construction became more coupled in the GBA during 2000–2018, as indicated by the continuously decreasing distance between the gravity centers of urbanization and eco-environment; and (3) as for the spatial differentiation characteristic of this coupling, the GBA was generally dominated by slightly uncoupled units, while the spatial distribution of different coupling types transformed from a circular structure to a relatively random form.


Introduction
In the past few decades, China has experienced unprecedented urbanization and economic development. According to a blue paper report released by Morgan Stanley in 2019, titled "The Rise of China's Supercities: New Era of Urbanization" [1], the country's urbanization rate will reach 75% by 2030, driven by city clusters, smart cities, and agricultural modernization (the growth of China's agricultural labor productivity over the next decade is expected to free up more rural population for further urbanization). This new era of urbanization (which is referred to as "Urbanization 2.0") is expected to increase China's urban population by 220 million by 2030, half of which are predicted to settle in China's five largest city clusters, namely the Yangtze River Delta, the Jing-Jin-Ji Area, the Guangdong-Hong Kong-Macao Greater Bay Area, the Mid-Yangtze River Area, and the Chengdu-Chongqing Area. However, population aggregation in big cities is always accompanied by for integrating DMSP/OLS and NPP/VIIRS data to evaluate the city light dynamics in Syria during the civil war (2011-present), which effectively extended the temporal coverage of NTL data.
The normalized difference vegetation index (NDVI) is another widely used type of remote sensing data in the study of socio-economic activities and environmental change. It reflects the vegetation growth state and vegetation coverage, which can be used to effectively detect spatiotemporal patterns of above-ground vegetation dynamics [38,39]. In particular, the annual variation of the NDVI has been extensively used to assess the impacts of human activities on the urban eco-environment [40].
With the incorporation of GBA into China's national strategies, this region has been recognized as a hotspot for future research. This study used NTL data and the NDVI to solve the problem of the incomparability of statistics between the inland part of the GBA and the Hong Kong-Macao area. The mean nighttime light (MNL) and the yearly maximum value of the NDVI were used to represent the level of urbanization and the level of eco-environment, respectively. Then, we adopted the gravity center model to analyze the overall coupling situation between urbanization and eco-environment and the coordination index to analyze the spatial differentiation characteristics of this coupling. This study was intended to (1) reveal the spatiotemporal patterns of urbanization and eco-environment in the GBA; (2) identify the overall coupling situation between urbanization and eco-environment change in the GBA; and (3) analyze the spatial differentiation characteristics of this coupling during different periods in the GBA. The results of this study are of important guiding significance to promote the coordinated development of urbanization and eco-environment in the GBA. Figure 1 is a flowchart showing the methodology of this study.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 21 the city light dynamics in Syria during the civil war (2011-present), which effectively extended the temporal coverage of NTL data. The normalized difference vegetation index (NDVI) is another widely used type of remote sensing data in the study of socio-economic activities and environmental change. It reflects the vegetation growth state and vegetation coverage, which can be used to effectively detect spatiotemporal patterns of above-ground vegetation dynamics [38,39]. In particular, the annual variation of the NDVI has been extensively used to assess the impacts of human activities on the urban eco-environment [40].
With the incorporation of GBA into China's national strategies, this region has been recognized as a hotspot for future research. This study used NTL data and the NDVI to solve the problem of the incomparability of statistics between the inland part of the GBA and the Hong Kong-Macao area. The mean nighttime light (MNL) and the yearly maximum value of the NDVI were used to represent the level of urbanization and the level of eco-environment, respectively. Then, we adopted the gravity center model to analyze the overall coupling situation between urbanization and eco-environment and the coordination index to analyze the spatial differentiation characteristics of this coupling. This study was intended to (1) reveal the spatiotemporal patterns of urbanization and eco-environment in the GBA; (2) identify the overall coupling situation between urbanization and eco-environment change in the GBA; and (3) analyze the spatial differentiation characteristics of this coupling during different periods in the GBA. The results of this study are of important guiding significance to promote the coordinated development of urbanization and eco-environment in the GBA. Figure 1 is a flowchart showing the methodology of this study. Nighttime

Study Area
The GBA is located at 111 • 12 -115 • 35 E and 21 • 25 -24 • 30 N and consists of nine cities in Guangdong Province (Guangzhou, Shenzhen, Foshan, Dongguan, Zhongshan, Zhuhai, Huizhou, Jiangmen, and Zhaoqing) and two Special Administrative Regions (Hong Kong and Macao) ( Figure 2). Its total area is 56,000 km 2 , and its population had reached 71.51 million by the end of 2018 [41].
As one of the most urbanized city agglomerations in the world, this region is expected to promote the cooperative development of the Pan-Pearl River Delta and construct a new open economic system for further integration into the world economy [16]. Up to now, two core areas-the Guangzhou-Foshan core area and the Hong Kong-Shenzhen core area, have emerged in the GBA, which will become new potential growth poles in the regional integration of the GBA. On the other hand, the boundary effect has reduced remarkably and border areas are growing rapidly in the GBA [42], especially the border regions of Shenzhen-Dongguan, Shenzhen-Huizhou, Dongguan-Huizhou, Guangzhou-Dongguan, and Guangzhou-Foshan. Under this development trend, the GBA will play a crucial role in the implementation of the national regional development strategy and the sustainable development of Hong Kong and Macao.
As one of the most urbanized city agglomerations in the world, this region is expected to promote the cooperative development of the Pan-Pearl River Delta and construct a new open economic system for further integration into the world economy [16]. Up to now, two core areasthe Guangzhou-Foshan core area and the Hong Kong-Shenzhen core area, have emerged in the GBA, which will become new potential growth poles in the regional integration of the GBA. On the other hand, the boundary effect has reduced remarkably and border areas are growing rapidly in the GBA [42], especially the border regions of Shenzhen-Dongguan, Shenzhen-Huizhou, Dongguan-Huizhou, Guangzhou-Dongguan, and Guangzhou-Foshan. Under this development trend, the GBA will play a crucial role in the implementation of the national regional development strategy and the sustainable development of Hong Kong and Macao.

Data Sources and Pre-Processing
The data used in this study included remote sensing data and administrative border data. The remote sensing data included NTL data and MODIS (Moderate-Resolution Imaging Spectroradiometer) NDVI data. The administrative border data included the boundary of the GBA

Data Sources and Pre-Processing
The data used in this study included remote sensing data and administrative border data. The remote sensing data included NTL data and MODIS (Moderate-Resolution Imaging Spectroradiometer) NDVI data. The administrative border data included the boundary of the GBA and the boundaries of the 620 township (street) units in the GBA. Both the NTL data and NDVI data were extracted according to the boundary of the GBA. DMSP/OLS data are provided in 30-arc-second grids, with only one band containing a greyscale value. The visible band has 6-bit quantization, producing a digital number (DN) value ranging from 0 to 63 for each pixel [43]. These data can be used to detect artificial lights, such as those from human settlements, fires, gas flares, and fishing boats, as well as lightning and the aurora [36,44]. Currently, the products of DMSP/OLS include three kinds of annual average data: average visible, stable lights, and cloud-free coverage. The night stable light is an annual cloud-free composite of the average DN value for the detected lights (e.g., lights from cities, towns, and other sites with persistent lighting) with ephemeral lights and background noise removed [45]. Therefore, we chose stable light data in our study. The data were not available subsequent to February 2014 [37]. In 2011, NASA (National Aeronautics and Space Administration, Washington, DC, USA) and NOAA (National Oceanic and Atmospheric Administration, Silver Spring, MD, USA) launched the Suomi National Polar Partnership (NPP) satellite carrying the first Visible Infrared Imaging Radiometer Suite (VIIRS) instrument. NPP/VIIRS data, which have a spatial resolution of 15 arc-seconds, include two data types: VIIRS Cloud Mask Configuration (VCMCFG) and VIIRS Cloud Mask Stray-light Corrected Configuration (VCMSLCFG). NPP/VIIRS data offer substantial improvements in spatial resolution, radiometric calibration, and usable dynamic range compared to the DMSP/OLS data [46].
This study used data from the years 2000, 2006, 2012, and 2018. We chose nighttime stable light of DMSP/OLS and VCMCFG of NPP/VIIRS in December, both of which were obtained from NGDC (National Genomics Data Center; https://ngdc.noaa.gov) ( Table 1). In order to match the spatial resolution of the two datasets, the DMSP/OLS data and NPP/VIIRS data were both resampled to a pixel size of 1km and projected into the Lambert Conformal Conic projection with reference to the WGS84 coordinate system. The negative radiation values in NPP/VIIRS images were set to 0.

MODIS NDVI Data
A MODIS data product (MOD13A2 data) in HDF format was acquired for 2000, 2006, 2012, and 2018 from the NASA Earth Observing System (EOS) data gateway (https://ladsweb.modaps.eosdis. nasa.gov/search/). This product is provided every 16 days at a 1-km spatial resolution. The strip number of the study area is h28v06. The NDVI data were projected into the Lambert Conformal Conic projection with reference to the WGS84 coordinate system and converted to TIFF format using the Modis Reprojection Tool (MRT).

Evaluation Indices of Urbanization and the Eco-Environment
We chose the MNL as the evaluation index of the level of urbanization. The MNL integrates information on the built-up area and the intensity of socio-economic activities, which reflects the average intensity of human activities in the study area [47]. The MNL was calculated as follows: where k represents the number of pixels in the analytical unit, DN i represents the DN value of the ith pixel. The processing of NTL data included calibration and integration, which will be further discussed in Section 3.2.
For the evaluation of the eco-environment, the yearly maximum values of the NDVI were extracted to reduce the influence of the atmosphere, aerosols, and clouds [48,49] (Figure S1). Then the mean value of the NDVI in each analytical unit was calculated to represent the level of eco-environment.

Calibration of DMSP/OLS Data
Since DMSP/OLS sensors have no onboard radiometric calibration for the visible band [50,51], the acquired data differ from each other in terms of radiometric performance. In this section, we followed the method proposed by Cao et al. [52] to achieve radiometric and saturation calibration of the DMSP/OLS data.
Firstly, we used the power function for intercalibration (Formula (2)): where DN orig(n,i) represents the original DN value of the ith pixel of the DMSP/OLS image in the nth year, DN cal(n,i) represents the intercalibrated DN value, and a and b are intercalibration parameters, which are listed in Table 2. Furthermore, in order to make full use of the DMSP/OLS data derived from multiple sensors in the same year, it was necessary to perform continuity calibration using Formula (3) if there were two images in the same year.
where DN 1 cal(n,i) and DN 2 cal(n,i) represent the intercalibrated DN values of the two NTL images, and DN con(n,i) represents the DN value after continuity calibration.
Finally, considering the accelerating urbanization process in the GBA since the start of the 21st century, further continuity calibration steps were performed using Formula (4) to ensure that the DN value of the same pixel keeps increasing over time: DN con(n+6,i) = 0 DN con(n−6,i) , DN con(n+6,i) < 0 and DN con(n−6,i) < DN con(n,i) where DN con(n,i) represents the DN value after the first step of continuity calibration, and DN (n,i) represents the DN value after further continuity calibration. The above calibration processes effectively solved the problem of pixel saturation in the DMSP/OLS data and greatly improved their continuity and comparability [44].

Integration of Nighttime Light Data
This study referred to the integration method proposed by Li et al. [37], which was used to simulate DMSP/OLS data from NPP/VIIRS data using a power function based on the temporal overlap between the two datasets. The integration process included the following three steps: (1) For each DN value, all pixels in the DMSP/OLS image from 2012 and the corresponding pixels in the NPP/VIIRS image for the same year were extracted. Then, the average radiation value of the extracted NPP/VIIRS pixels was calculated (Table S1, Supplementary Material), and a scatterplot was generated with the average radiation value of NPP/VIIRS data as the independent variable (the unit is nW/cm 2 /s) and the DN value of DMSP/OLS data as the dependent variable. We selected the pixels whose DN value was equal to or less than 50 to reduce the influence of pixel saturation [37]. The fitted curve for 2012 is shown in Figure 3. The integration function is given in Formula (5): where radiance represents the radiation value of the NPP/VIIRS data, and DN represents the DN value simulated from the NPP/VIIRS data. (2) Formula (5) was used to simulate DMSP/OLS data for 2018. Then a Gaussian low-pass filter with a window size of 3 pixels was used to smooth the simulated image. (3) Finally, the continuity calibration was performed again using Formula (4) to ensure the continuity of the whole dataset.
where ( , ) represents the DN value after the first step of continuity calibration, and ( , ) represents the DN value after further continuity calibration. The above calibration processes effectively solved the problem of pixel saturation in the DMSP/OLS data and greatly improved their continuity and comparability [44].

Integration of Nighttime Light Data
This study referred to the integration method proposed by Li et al. [37], which was used to simulate DMSP/OLS data from NPP/VIIRS data using a power function based on the temporal overlap between the two datasets. The integration process included the following three steps: (1) For each DN value, all pixels in the DMSP/OLS image from 2012 and the corresponding pixels in the NPP/VIIRS image for the same year were extracted. Then, the average radiation value of the extracted NPP/VIIRS pixels was calculated (Table S1, Supplementary Material), and a scatterplot was generated with the average radiation value of NPP/VIIRS data as the independent variable (the unit is nW/cm 2 /s) and the DN value of DMSP/OLS data as the dependent variable. We selected the pixels whose DN value was equal to or less than 50 to reduce the influence of pixel saturation [37]. The fitted curve for 2012 is shown in Figure 3. The integration function is given in Formula (5): where radiance represents the radiation value of the NPP/VIIRS data, and DN represents the DN value simulated from the NPP/VIIRS data.
(2) Formula (5) was used to simulate DMSP/OLS data for 2018. Then a Gaussian low-pass filter with a window size of 3 pixels was used to smooth the simulated image. (3) Finally, the continuity calibration was performed again using Formula (4) to ensure the continuity of the whole dataset.

The Gravity Center Model
The gravity center model [9] was adopted to analyze the overall coupling situation of the GBA. The calculation formulas for the gravity centers of urbanization and eco-environment are shown in Formulas (6) and (7), respectively: where C U and C E represent the gravity centers of urbanization and eco-environment, respectively, U i and E i represent the value of MNL and NDVI of the ith township (street) unit, respectively, and C(x i , y i ) are the coordinates of the gravity center of the ith township (street) unit. The gravity center model consists of two parts: the overlapping of gravity centers and the consistency of their movement tracks.

Overlapping of Gravity Centers
The overlapping of gravity centers was calculated by the distance between the two gravity centers, with a shorter distance implying a higher level of coordination. The distance was calculated using Formula (8): where x U (x E ) and y U (y E ) represent the x-coordinate and y-coordinate of the gravity center of urbanization (eco-environment), respectively.

The Consistency of Movements
The consistency of movements was defined as the cosine of the angle (θ) between the movement tracks of the two gravity centers, which was calculated using Formula (9) according to the cosine law: where ∆ x U (∆ x E ) and ∆ y U (yx E ) represent the variations of the x-coordinate and y-coordinate of the gravity center of urbanization (eco-environment), respectively, during the specific period. The range of C was [−1,1]; C = 1 implies consistent directions of the two movement tracks, and C = −1 implies opposite directions.

The Coordination Index
Aside from the gravity center model, we adopted the coordination index [7] to evaluate the spatial differentiation characteristics of the coupling relationship between urbanization and eco-environment in the GBA. The coordination index was defined as follows: where C UE represents the coordination index of urbanization and eco-environment change, and UR and ER represent the average annual change rates of the MNL and NDVI, respectively. The range of C UE was [0,1]. When UR was equal to ER, C UE was 1, and urbanization and eco-environment change were maximally coupled; when the absolute value of UR was equal to ER but of opposite sign, C UE was 0, and urbanization and eco-environment change were minimally uncoupled. The classification of C UE and the relationship between UR and ER are shown in Table 3 (adapted from [2,7], Table S2). Table 3. Classification of the coordination index of urbanization and eco-environment change (C UE ) and the relationship between the average annual change rate of the mean nighttime light (MNL) (UR) and the average annual change rate of the normalized difference vegetation index (NDVI) (ER).

Analysis of the Pattern of Urbanization
The MNL was calculated for each township (street) unit using Formula (1). Figure 4 shows the spatial and temporal variations of the urbanization level in the GBA. Between 2000 and 2018, the urbanization pattern mainly presented the following two characteristics.
where CUE represents the coordination index of urbanization and eco-environment change, and UR and ER represent the average annual change rates of the MNL and NDVI, respectively. The range of CUE was [0,1]. When UR was equal to ER, CUE was 1, and urbanization and eco-environment change were maximally coupled; when the absolute value of UR was equal to ER but of opposite sign, CUE was 0, and urbanization and eco-environment change were minimally uncoupled. The classification of CUE and the relationship between UR and ER are shown in Table 3 (adapted from [2,7], Table S2).

Analysis of the Pattern of Urbanization
The MNL was calculated for each township (street) unit using Formula (1). Figure 4 shows the spatial and temporal variations of the urbanization level in the GBA. Between 2000 and 2018, the urbanization pattern mainly presented the following two characteristics.

A Circular Structure
On the whole, there was a large spatial differentiation of the urbanization level in the GBA (Figure 4). The high-value units were concentrated in the central area, such as the border of Guangzhou and Foshan, northern Dongguan, and most parts of Shenzhen and Hong Kong, while the low-value units were mainly distributed at the periphery, which showed a circular structure following the gradient of urbanization level. This significant imbalance was caused by several factors. Firstly, the central area is conveniently connected to other cities of the Pearl River Delta by the efficient transportation system, while the mountainous western and eastern regions are separated from central cities and mainly function as ecological barriers. Secondly, the central area has attracted a number of foreign enterprises to invest since China's reform and opening-up and has become the first production base in the "pre-shop and post-factory" model [53]. In recent years, the industrial structure in cities in the central GBA has been continuously optimized and a high value-added modern industrial system has been gradually formed; however, peripheral cities are still in the early stages of industrialization. As a result, the gap between the center and the periphery of the GBA has further widened (Figure 4d).

Rise in Urban Border Areas
In recent years, the urban border areas took advantage of the cheap land and convenient transportation to undertake core industry diffusion and industrial transfer. Economic activities in the urban border areas and the cross-city cooperation between towns and streets were rather frequent [54]. As shown in Figure 4, it was found that, between 2000 and 2018, the urbanization level of border areas such as Guangzhou-Foshan, Guangzhou-Dongguan, Shenzhen-Dongguan, Shenzhen-Huizhou, and Dongguan-Huizhou significantly increased, which was closely related to the driving effects of core cities such as Guangzhou and Shenzhen.
Additionally, the change rates of the urbanization level were also found to vary greatly among different units during the three periods ( Figure 5). Most units developed stably during the research period, while western Zhaoqing and eastern Huizhou developed significantly during the first and the third period, especially from 2012 to 2018. The periphery area remained underdeveloped for a long time due to its mountainous terrain ( Figure S3). However, with the large-scale land development policy enacted since the start of the 21st century, most cities in the GBA, especially those in less developed regions, began to exploit their underdeveloped areas through land reclamation and expropriation (Figure 5a). By 2006, the concepts of "ecological cities" and "circular economy" were adopted in most cities in China, including those in the GBA. Therefore, urban construction entered into a relatively steady development period (Figure 5b). After 2012, Zhaoqing and Huizhou continued to exploit idle land to expand the border of the urban core area, which resulted in the rapid expansion of urban land in the peripheries of the two cities (Figure 5c).

Analysis of the Pattern of Eco-Environment
The mean value of the yearly maximum value of the NDVI was calculated for each township (street) unit. Similar to the spatial pattern of urbanization (Section 4.1), the spatial pattern of the ecoenvironment of the township (street) units showed a circular structure ( Figure 6). Areas with a high level of eco-environment were mainly concentrated in the periphery of the GBA, such as Zhaoqing, northern Huizhou, and western Jiangmen; these areas are dominated by mountains and hills and perform key ecological functions in the cities. On the contrary, in urban areas in the central GBA, the

Analysis of the Pattern of Eco-Environment
The mean value of the yearly maximum value of the NDVI was calculated for each township (street) unit. Similar to the spatial pattern of urbanization (Section 4.1), the spatial pattern of the eco-environment of the township (street) units showed a circular structure ( Figure 6). Areas with a high level of eco-environment were mainly concentrated in the periphery of the GBA, such as Zhaoqing, northern Huizhou, and western Jiangmen; these areas are dominated by mountains and hills and perform key ecological functions in the cities. On the contrary, in urban areas in the central GBA, the eco-environment was greatly degraded between 2000 and 2018 as a consequence of urban land exploitation (Figure 6d), and therefore the level of eco-environment was generally low in these areas compared to peripheral areas. The change rates of the eco-environment also showed significant spatial and temporal differentiation (Figure 7). From 2000 to 2006, except for the central district of Guangzhou, the ecoenvironment in most central cities showed a significantly deteriorating trend (Figure 7a). During this period, the central district of Guangzhou exhibited a unique pattern of stably increasing ecoenvironment levels while being surrounded by units with the deteriorating eco-environment. This phenomenon was closely related to policies involving population redistribution and industrial relocation outward from the central district [55]. After 2006, the level of eco-environment improved greatly in central areas, which corresponded to the prevalence of ecological city construction in China (Figure 7b,c). On the other hand, in most peripheral areas, the level of eco-environment increased steadily during the research period, except for central and western Zhaoqing, which showed a slight drop during 2006-2012 (Figure 7b).

Overall Coupling Situation
This study adopted the gravity center model to evaluate the overall coupling situation between urbanization and eco-environment change in the GBA. The overlapping of gravity centers and the consistency of movements were calculated as the distance between the two gravity centers (D) and the cosine of the angle (θ) between their movement tracks (C), respectively ( Table 4). The two movement tracks are shown in Figure 8.

Overall Coupling Situation
This study adopted the gravity center model to evaluate the overall coupling situation between urbanization and eco-environment change in the GBA. The overlapping of gravity centers and the consistency of movements were calculated as the distance between the two gravity centers (D) and the cosine of the angle (θ) between their movement tracks (C), respectively ( Table 4). The two movement tracks are shown in Figure 8. The locations of the gravity centers of urbanization and eco-environment both changed significantly during the research period. On the one hand, in southern Guangzhou (Huangge town of Nansha district), the gravity center of urbanization moved from the southeast to the northwest with an average annual rate of 0. better-protected eco-environment of the northwestern area, while the unusual phenomenon from 2006 to 2012 might reflect the remarkable effect of ecological restoration and recovery in the central area resulting from environmental conservation policies during this period. On the whole, the variation of both gravity centers was largest during 2012-2018, indicating the rapid northwestward expansion of urban development and green space construction in the GBA. Therefore, during the research period, Zhaoqing continued to strengthen its status as an ecological barrier while maintaining sufficient momentum of urban development.
In general, during the research period, the distance between the two gravity centers decreased continuously with an increasing rate, suggesting that the urbanization process and eco-environment change became increasingly coupled in the whole GBA. From 2000 to 2006, the distance between the two gravity centers decreased from 47.39 to 46.12 km, and the movement tracks were almost parallel, with a C of 0.87. Therefore, both the "hotspots" of urbanization development and eco-environment construction were concentrated in the northwestern area during this period. During 2006-2012, the two gravity centers became closer, however, the negative value of C implied opposite movement directions of the gravity centers. The gravity center of urbanization continued to move northwestward, while the movement track of the gravity center of the eco-environment turned eastward, reflecting different spatial patterns of urbanization and eco-environment change. This difference was partly due to the different foci of relevant policies. From 2012 to 2018, the distance decreased further at an average annual rate of 0.42 km/a, and the movement track of the gravity center of the eco-environment turned northwestward again with a significantly increased C of 0.85. On the whole, both the variation of the distance between the two gravity centers and the consistency of their movements were largest between 2012 and 2018, when urban development and environmental construction became more coupled and spatial patterns of the two processes were generally consistent. This can be attributed to the incorporation of the GBA into the national strategy and relevant policies targeting the area as a whole during this period. On the whole, the variation of both gravity centers was largest during 2012-2018, indicating the rapid northwestward expansion of urban development and green space construction in the GBA. Therefore, during the research period, Zhaoqing continued to strengthen its status as an ecological barrier while maintaining sufficient momentum of urban development.
In general, during the research period, the distance between the two gravity centers decreased continuously with an increasing rate, suggesting that the urbanization process and eco-environment change became increasingly coupled in the whole GBA. From 2000 to 2006, the distance between the two gravity centers decreased from 47.39 to 46.12 km, and the movement tracks were almost parallel, with a C of 0.87. Therefore, both the "hotspots" of urbanization development and eco-environment construction were concentrated in the northwestern area during this period. During 2006-2012, the two gravity centers became closer, however, the negative value of C implied opposite movement directions of the gravity centers. The gravity center of urbanization continued to move northwestward, while the movement track of the gravity center of the eco-environment turned eastward, reflecting different spatial patterns of urbanization and eco-environment change. This difference was partly due to the different foci of relevant policies. From 2012 to 2018, the distance decreased further at an average annual rate of 0.42 km/a, and the movement track of the gravity center of the eco-environment turned northwestward again with a significantly increased C of 0.85. On the whole, both the variation of the distance between the two gravity centers and the consistency of their movements were largest between 2012 and 2018, when urban development and environmental construction became more coupled and spatial patterns of the two processes were generally consistent. This can be attributed to the incorporation of the GBA into the national strategy and relevant policies targeting the area as a whole during this period.

Spatial Differentiation Characteristics of Coupling
This study adopted the coordination index to analyze the spatial differentiation characteristics of coupling between urbanization and eco-environment change in the GBA. The 620 township (street) units were classified into five categories-namely superiorly coupled, barely coupled, slightly uncoupled, moderately uncoupled, and seriously uncoupled-according to the value of the coordination index.  Figure 10.
On the whole, the GBA was dominated by slightly uncoupled units, while the number of moderately and seriously uncoupled units decreased over time ( Figure 11). A comparison between Figures 9 and 10 showed that units in the stage of uncoupled development were mostly generated by higher UR vs. ER (UR − ER > 0.1, eco-environment construction lag behind urbanization process), while coupled units showed a relatively synchronous development of urbanization and eco-environment (0 ≤ |UR − ER| ≤ 0.1).
From 2000 to 2006, the spatial distribution of different coupling types generally showed a circular structure, transitioning from seriously uncoupled to moderately uncoupled to slightly uncoupled to barely coupled to superiorly coupled units when moving from the center to the periphery (Figure 9a), which approximately paralleled the spatial pattern of urbanization and eco-environment during this period (Figures 4 and 6). In western Zhaoqing, western Jiangmen, and northern Guangzhou, urban development lagged behind other cities for a long time. In these areas, the original land-use types were mostly retained and the carrying capacity of the eco-environment was relatively high. Therefore, in these regions, the level of urbanization and eco-environment both increased steadily, and the urban expansion and eco-environment construction were more coordinated than in other areas. However, in central areas such as southern Guangzhou, Foshan, and Dongguan, land-development activities starting at the end of the 20th century largely reduced the original natural space, leading to environmental deterioration and the continuing unbalanced development between urban expansion and environmental protection. The only exception was the central district of Guangzhou, where a few coupled units were surrounded by a number of moderately and seriously uncoupled units. On the one hand, the process of urbanization in this area nearly reached a plateau and the construction of green space was greatly constrained by the widespread impervious ground. Therefore, the change rates of urbanization and eco-environment in this district were both limited and basically synchronous. On the other hand, however, as a consequence of the limited space for urban sprawl in the central district of Guangzhou, surrounding areas were forced to expand construction land rapidly. Therefore, in these areas, the change rates of the urbanization level significantly exceeded the change rates of the eco-environment level (Figure 10a    From 2006 to 2012, the percentages of superiorly coupled, barely coupled, and slightly uncoupled units increased respectively from 13%, 17%, and 45% to 17%, 22%, and 49%, while the percentage of moderately and seriously uncoupled units decreased respectively from 16% and 9% to 9% and 3%. On the whole, most seriously and moderately uncoupled units during 2000-2006 were converted to coupled and slightly uncoupled units (Figure 9b). In this period, serious environmental problems caused by previous large-scale urban development aroused common public concerns regarding ecological protection and urban landscape planning in central urban areas. With the proposal of the sustainable development strategy and ecological civilization construction, the restoration and conservation of urban green space effectively protected the environment from continuous destruction by the expansion of construction land (Figure 10b). Therefore the process of urbanization and eco-environment change became more coupled during this period. However, note that a number of moderately uncoupled units appeared in central Zhaoqing, which were converted from barely coupled and slightly uncoupled units. With the process of urbanization, the expansion From 2006 to 2012, the percentages of superiorly coupled, barely coupled, and slightly uncoupled units increased respectively from 13%, 17%, and 45% to 17%, 22%, and 49%, while the percentage of moderately and seriously uncoupled units decreased respectively from 16% and 9% to 9% and 3%. On the whole, most seriously and moderately uncoupled units during 2000-2006 were converted to coupled and slightly uncoupled units (Figure 9b). In this period, serious environmental problems caused by previous large-scale urban development aroused common public concerns regarding ecological protection and urban landscape planning in central urban areas. With the proposal of the sustainable development strategy and ecological civilization construction, the restoration and conservation of urban green space effectively protected the environment from continuous destruction by the expansion of construction land (Figure 10b). Therefore the process of urbanization and eco-environment change became more coupled during this period. However, note that a number of moderately uncoupled units appeared in central Zhaoqing, which were converted from barely coupled and slightly uncoupled units. With the process of urbanization, the expansion of urban construction land started to encroach the original natural green space, resulting in the somewhat uncoupled development of urbanization and eco-environment construction in this area.  From 2012 to 2018, the percentage of slightly uncoupled units increased from 49% to 69%, while the percentage of other coupling types all decreased slightly. Most areas were dominated by slightly uncoupled development units, which were mainly converted from superiorly and barely coupled units. Moreover, some seriously uncoupled units appeared in Macao and the border area of Shenzhen and Hong Kong (Figure 9c). After 2012, the GBA gradually became a key growth pole in China's implementations of innovation-driven development and commitment to reform and opening-up, and the integration of urban agglomerations in the GBA became the basic goal from the perspective of regional planning. However, due to the immature mechanism of cooperation in cross-border areas, ecological protection was overlooked, and thus the environment greatly deteriorated during this period, especially in the previously underdeveloped Yuen Long District in northern Hong Kong. In peripheral areas such as Zhaoqing, northern Huizhou, and western Jiangmen, urban construction From 2012 to 2018, the percentage of slightly uncoupled units increased from 49% to 69%, while the percentage of other coupling types all decreased slightly. Most areas were dominated by slightly uncoupled development units, which were mainly converted from superiorly and barely coupled units. Moreover, some seriously uncoupled units appeared in Macao and the border area of Shenzhen and Hong Kong (Figure 9c). After 2012, the GBA gradually became a key growth pole in China's implementations of innovation-driven development and commitment to reform and opening-up, and the integration of urban agglomerations in the GBA became the basic goal from the perspective of regional planning. However, due to the immature mechanism of cooperation in cross-border areas, ecological protection was overlooked, and thus the environment greatly deteriorated during this period, especially in the previously underdeveloped Yuen Long District in northern Hong Kong. In peripheral areas such as Zhaoqing, northern Huizhou, and western Jiangmen, urban construction developed rapidly, and ecological protection projects were promoted steadily (Figure 10c). As a result, urbanization and eco-environment change remained coupled or slightly uncoupled during this period, and the level of coupling decreased slightly.

Conclusions
As a hotspot area in China's national strategies, the GBA has attracted public attention from various research fields in recent years. Recent research on the GBA has mainly focused on economic advancement and policy-making. Taking the GBA as the example, we adopted a gravity center model and the coordination index to analyze the overall coupling situation between urbanization and eco-environment and the spatial differentiation characteristics of this coupling, respectively. The results showed that the process of urbanization and eco-environment change in the GBA became increasingly coupled between 2000 and 2018, as indicated by the decreasing distance between the two gravity centers which both generally moved northwestward. Specifically, the spatial pattern of coupling transformed from a circular structure to a more or less random form, which corresponded to the adjustment of policies and development paths in some cities in the GBA.
Remote sensing data have been widely recognized as an effective tool in research on urbanization and eco-environment, especially in the case of cross-regional studies where the utility of socio-economic statistics was always hampered by the problem of incomparability. This research adopted nighttime light data (DMSP/OLS data and NPP/VIIRS data) and MODIS NDVI data to evaluate the level of urbanization (MNL) and eco-environment (NDVI), respectively. The spatial patterns of urbanization and eco-environment both presented a circular structure, and their change rates showed significant spatial and temporal differentiation. These findings correspond to the empirical evidence presented in previous studies of the urbanization process and eco-environment change in the GBA [13,56,57], and in turn, corroborate the applicability of nighttime light data and MODIS NDVI for quantitatively evaluating the level of urbanization and eco-environment as well as their spatiotemporal differentiation.
The analysis of the overall coupling situation suggested that, between 2000 and 2018, urban development and eco-environmental construction became more coupled in the whole area, but the spatial and temporal variations of the coordination index showed some contrasting patterns. From 2000 to 2006, the spatial distribution of different coupling types was mostly a legacy of policy orientation and land-development activities from the last century. During this period, most central areas were transformed into modern cities accompanied by environmental deterioration, while peripheral areas remained almost undeveloped. Severe environmental issues began to raise public attention and many developed areas began to emphasize the coordinated development of urban construction and environmental protection, which significantly mitigated the previous adverse outcomes. Therefore, the number of coupled units increased during 2006-2012.
However, from 2012 to 2018, most moderately and superiorly coupled units were substituted by slightly uncoupled units, which contradicts the overall trend of strong coupling between urbanization and eco-environment. We suggested that this represents a "false-coupled" phenomenon, where the number of coupled units decreased in the whole region, while the overall coupling situation showed an improving trend instead. The gravity center model treated the region as a whole and just reflected the average coupling situation during a specific period. On the other hand, the coordination index quantitatively evaluated the coupling stage of each analytical unit and further revealed the spatial differentiation characteristics of the overall coupling situation. This "false-coupled" phenomenon, in fact, is probably misleading when achieving the overall regional development goal is taken as the only golden rule in the policy-making process. Therefore, understanding this phenomenon is essential for guiding sustainable urban development in the GBA. When formulating regional development plans, local governments should pay more attention to some critical areas (for example, areas with a fragile environment, areas with a low level of urbanization that lags behind most other areas, and areas with uncoordinated development of urbanization and eco-environment) and develop a differentiation strategy targeting different areas to achieve coordinated and integrated planning for the GBA.  Figure S3: The elevation map of the GBA, Table S1: The DN value of DMSP/OLS and the average radiation value of corresponding pixels of NPP/VIIRS in the images in 2012, Table S2: The classification of coupling types in References [1,5].