Detection of City Integration Processes in Rapidly Urbanizing Areas Based on Remote Sensing Imagery

Since China’s reform and development commenced, in the context of rapid urbanization and coordinated regional development, Chinese cities with a close geographic proximity and social ties have gradually formed an integrated city development model. As a new phenomenon in China’s urbanization process, existing research on China’s integrated cities mainly focuses on typical case studies, and most research has been limited to literature reviews and theoretical analyses. The growing application of remote sensing technology in urbanization research in recent years has provided new opportunities for the analysis of city integration. Therefore, based on multi-spectral Landsat-8 and nighttime light images (SNPP/VIIRS, Suomi National Polar-orbiting Platform/Visible Infrared Imaging Radiometer Suite), this paper selects four of the most representative integrated cities with different backgrounds in China to analyze the land-use conversion, plot light fluctuation, and light gravity center shift in the boundary zone between cities. The results show that (1) Guangfo has the highest level of integration and urban expansion is mainly concentrated in the south-central part of the boundary area; (2) Guanshen’s level of integration is second to Guangfo’s and is mainly concentrated in the west; (3) HuSu’s integration is still in the initial stage and its increase in light intensity lags behind the expansion of building land during the study period; (4) although the light intensity and building land area increased significantly during the study period in Xixian, the overall development level of Xixian still lagged behind coastal cities due to the restriction of its geographical location. Our application results expand the data sources for integrated city research and the obtained results can potentially support decision-making and planning in the process of urban development.


Introduction
Since China's reform and opening up policies commenced in 1978, the country has experienced a rapid process of urbanization, with the urbanization rate rising from 17.9% in 1978 to 60.34% in 2020 [1]. Considering that there is still a gap between China's current urbanization rate and the 80% of developed countries, it is foreseeable that China will continue to further promote the urbanization process in the next few decades. On the one hand, continuous and high-intensity urbanization has become the engine that drives China's rapid economic growth, but on the other hand, problems such as over-concentration of population, tight urban land resources, and ecological environmental impacts that appear in the process of urbanization also constantly restrict and affect the sustainable development of cities [2][3][4]. In order to alleviate the pressure on urban land resources during the urbanization process, the original boundaries of the city are constantly expanding into surrounding With China's increasing urbanization, these four groups of cities are experiencing rapid urban expansion (Table 1), and the economic, social, natural, and ecological ties between them are becoming close and developing towards city integration. As a result of the process of integration of neighboring cities, the built-up areas between cities become interconnected and eventually develop into a whole part, also known as a megacity or urban agglomeration. Therefore, as the core area of the integration process, the boundary area of adjacent cities is an excellent window to observe this urban integration. As shown in Figure 1, we especially focus on the border area among the integrated cities, and construct 1 km, 3 km, and 5 km buffer zones at the borders of neighboring cities as the target areas for remote sensing monitoring analysis. With China's increasing urbanization, these four groups of cities are experiencing rapid urban expansion (Table 1), and the economic, social, natural, and ecological ties between them are becoming close and developing towards city integration. As a result of the process of integration of neighboring cities, the built-up areas between cities become interconnected and eventually develop into a whole part, also known as a megacity or urban agglomeration. Therefore, as the core area of the integration process, the boundary area of adjacent cities is an excellent window to observe this urban integration. As shown in Figure 1, we especially focus on the border area among the integrated cities, and construct 1 km, 3 km, and 5 km buffer zones at the borders of neighboring cities as the target areas for remote sensing monitoring analysis.

Data Source
In this paper, the dataset used in the integrated analysis mainly includes four types: nighttime light data, multi-spectral image data, high-resolution satellite street view data, and administrator boundary datasets.
Land 2020, 9, 378 4 of 15 Considering the scale of the study area, we selected the nighttime light data of the Suomi National Polar Orbiting Partnership/Visible Infrared Imaging Radiometer Suite (SNPP/VIIRS) image. As a new generation of nighttime light datasets, VIIRS has significantly improved spatial resolution and data quality compared to a Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS). SNPP/VIIRS monthly average light data filter out the stray light, lightning, moonlight, and cloud cover, and the spatial resolution is around 500 m. Specifically, in this paper, the SNPP/VIIRS nighttime light data used was the monthly VIIRS Cloud Mask product (VCM) from April 2012 to May 2018, with a total of 74 images.
In addition, Landsat-8 multi-spectral images in 2013, 2015, and 2017 were used to monitor land cover change in the integrated city boundary zone. Since the spatial resolution of the Landsat-8 multi-spectral band is limited to 30 m, it cannot further support the requirements for the identification of typical features. Therefore, we used time-series high-resolution satellite image on the Google Earth platform to identify the typical plots within the boundary zone.
In addition, the administrative boundary data for the study area is from GADM Version 3.6 (https://gadm.org/). We combined the Open Street Map (OSM) data with partial corrections to ensure the accuracy of the study area boundaries.

Annual Nighttime Light Composite
As a means of reducing the interference of nighttime light fluctuations with the results, we synthesized annual nighttime light data based on monthly VCM data according to Equation (1), specifically: where VCM year is an annual synthesis of VIIRS nighttime light images; i is the month of images; VCM i is the nighttime light image of month i; k is the selected month; K is total number of selected months. It should be noted that due to the inevitable cloud cover in some months, the monthly images with abnormally low-value light radiation were excluded from the annual composite products. In addition, to minimize the instability of the single-year VCM data, this paper smoothed the selected light images according to Equation (2): where VCM year A is the mean-smoothed annual synthesized VCM of year A; VCM year a−1 , VCM year a , and VCM year a+1 are annual light images for a total of three years before and after year A, respectively. In this paper, we set year A as 2013, 2015, and 2017, and then calculated the difference between the smoothed VCM to get the nighttime light intensity changes of the integrated cities boundary zones.

Land Surface Feature Extraction
Firstly, this paper performed the necessary mosaic cropping, geometric correction, and atmospheric correction of the acquired Landsat 8 images of the study area for 2013, 2015, and 2017. Based on the time-series Landsat images, the spectral index method was used to extract feature types [42]. Specifically, water bodies were extracted using the Modified Normalized Difference Water Index (MNDWI, Equation (3)), vegetation was extracted using the Normalized Difference Vegetation Index (NDVI, Equation (4)), and land for construction was extracted using the Normalized Difference Building Index (NDBI, Equation (5)). Finally, we used high-resolution Google Earth images to evaluate the accuracy of land where ρ Green , ρ Red , ρ NIR , and ρ SWIR1 are the reflectance of the green band (Band 3), red band (Band 4), near infrared band (Band 5), and short-wave infrared band1 (Band 6) of the Landsat-8 image, respectively.

Gravity Model Analysis of Nighttime Light
With the purpose of measuring the spatially fluctuating characteristics of nighttime light intensity in the integrated city boundary zone, the nighttime light gravity center (NTLGC) model was introduced to reveal the horizontal shifting of nighttime light. The NTLGC was defined as: where G VCM (X, Y) is the gravity center coordinate of nighttime light; DN i is the nighttime light intensity value of the light pixel i; G(x i , y i ) is the center coordinate of the light pixel i. Moreover, the shift distance between the NTLGCs at different periods can indicate the strength and direction of nighttime light migration within the boundary zone. During the study period, the closer the shift distance of the gravity center is, the more stable the light changes in space; the farther the distance is, the stronger the difference of night light changes in space.

Spatial Correlation Analysis
Quantifying the relationship between surface vegetation and impervious surface with nighttime lighting changes will help reveal the impact of integrated cities on the regional environment. Therefore, ordinary least squares (OLS) and geographically weighted regression (GWR) were used to conduct the correlation analysis in this paper.
OLS is a commonly used statistical method for analyzing the interdependence between interpreted and explanatory variables, which can be calculated as: where y i is the dependent variable value of point i; β 0 is the intercept; x ik is the value of the independent variable k at point i; β k is the regression coefficient of the independent variable k; ε i is the residual variable.
Since OLS models may ignore the spatial instability of the relationship between the independent and dependent variables, this paper also introduced the GWR model, which adds geographic locations to the regression parameters and allows for local parameter estimation [39]. The GWR model is as follows: Land 2020, 9, 378 x ik is the impact factors; ε i is the random errors terms.

Land-Use Change in Boundary Zones
During the integration process, the downtown area of neighboring cities continues to expand and the boundaries of the built-up areas approach the administrative boundaries, leading to significant adjustments in the land-use structure characterizing the urban boundary zones. We extracted the land use within the multi-level buffer zones constructed along the boundary lines of four groups of integrated cities, the results of this extraction are shown in Figure 2.
Land 2020, 9, x FOR PEER REVIEW 6 of 16 During the integration process, the downtown area of neighboring cities continues to expand and the boundaries of the built-up areas approach the administrative boundaries, leading to significant adjustments in the land-use structure characterizing the urban boundary zones. We extracted the land use within the multi-level buffer zones constructed along the boundary lines of four groups of integrated cities, the results of this extraction are shown in Figure 2. Based on the classification results, this paper counts the percentage of area within the boundary zone for 1 km, 3 km, and 5 km widths, respectively, as shown in Table 2. During 2013-2017, the surface cover of the above boundary zones of four integrated cities changed significantly. Among them, (1) the building area expanded markedly, with Xixian growing the most and Guanshen growing the least; (2) both vegetation cover and water body were occupied, with the average cover decreasing by 49.59% and 10.61%, respectively. In other words, roughly one-sixth of the new construction land was converted from water cover and one-fifth was converted from vegetation cover. Besides, for the expansion rate of the building area within the boundary zone of different widths, the growth rate of Guangfo and Guanshen was more stable, while the growth rate of Husu and Xixian decreased with the buffer distance increases. The two different growth characteristics above reflect the different stages of integration development: the integration process of Guangfo and Guanshen has entered the stage of overall expansion, while Husu and Xixian are still in the initial stage of local expansion.  Based on the classification results, this paper counts the percentage of area within the boundary zone for 1 km, 3 km, and 5 km widths, respectively, as shown in Table 2. During 2013-2017, the surface cover of the above boundary zones of four integrated cities changed significantly. Among them, (1) the building area expanded markedly, with Xixian growing the most and Guanshen growing the least; (2) both vegetation cover and water body were occupied, with the average cover decreasing by 49.59% and 10.61%, respectively. In other words, roughly one-sixth of the new construction land was converted from water cover and one-fifth was converted from vegetation cover. Besides, for the expansion rate of the building area within the boundary zone of different widths, the growth rate of Guangfo and Guanshen was more stable, while the growth rate of Husu and Xixian decreased with the buffer distance increases. The two different growth characteristics above reflect the different stages of integration development: the integration process of Guangfo and Guanshen has entered the stage of overall expansion, while Husu and Xixian are still in the initial stage of local expansion.

Nighttime Light Intensity Changes
Since nighttime light data are superior in depicting the intensity of human activity (especially in downtown areas), it has been widely used in regional or global-scale urbanization studies. Hence, we carried out the fluctuation analysis of nighttime light intensity in the boundary zone to detect the impact of human activities on the integration process ( Table 3). The results show that from 2012 to 2018, the average light intensity within the boundary zone maintained an upward trend, with an increase of 35%, 44%, 36%, and 73% for Guangfo, Guanshen, Husu, and Xixian, respectively. Although Xixian achieved the highest light growth rate (over 70%), the absolute value of light intensity is still apparently lower than that of other integrated cities. In recent years, although the urbanization of Xixian, which is located in the hinterland of northwest China, has been accelerating, the overall development level there is still lagging behind the eastern and southern coastal cities due to the limitation of its non-coastal geographical location and the phased strategy of reform and development. In addition, the fitted equations and coefficient of determination R 2 for light changes in Table 3 similarly show that the increase in nighttime light intensity is more stable in the integrated coastal cities than in the integrated inland cities. To investigate the triggers of light fluctuations in the boundary area, we selected plots with remarkable changes in light intensity in the study area and combined Google Earth high-resolution satellite images for comparative analysis (Figure 3). Taking Guangfo as an example, the area and distribution density of the light increase area within the boundary zone is larger than the attenuation area. The plots with reduced light intensity in Guangfo's boundary zone were mainly: unoccupied arable land (a), river bridge blocks (b, c, d, and h), and idle land in urban areas (e, f, and g). Plots with increasing light intensity could be divided into transportation land, industrial land, and commercial land according to the type of use. Among them, Plot 1 is mainly composed of three industrial parks and two logistics parks; Plot 2 includes a large interchange, two shopping malls, and a logistics center; Plot 3 is located at the junction of the three districts and includes two industrial parks and a railway distribution center; Plot 4 is located in Foshan and includes two commercial sites; Plot 5 includes two industrial zones; Plot 6 includes one industrial land area and two commercial land areas; Plot 7 is located in Guangzhou, which is an industrial area mainly used for garment making; Plot 8 mainly includes two industrial parks in Foshan; Plot 9 is located on the boundary line, which includes Guangzhou South Station; Plot 10 is located in Guangzhou and contains two industrial zones. The comparison shows that the land use types in the light attenuation area are mainly vegetation and water bodies, while the light increase area is mainly industrial land (industrial parks and logistics parks), transportation land, and some commercial land. Similarly, Guanshen's light increase area is mainly located in industrial parks, and its light attenuation area is reservoirs, some urban villages and gardens, etc.; Husu's light attenuation zone mainly includes water bodies and green land, while the light increase plots mainly are new residential communities, industrial parks, and transportation facilities; Xixian's light attenuation area mainly includes land to be developed and green land within the city, while the light increase plots are mainly used for industrial and traffic facilities. Generally, in the boundary areas of integrated cities, logistics parks and industrial parks, driven by location and land prices, continue to attract concentrations of people and goods, and the intensity of human activity then increases significantly. Conversely, for cultivated and green land within a border zone, a decreasing trend in light intensity is observed, and we speculated that the shading effect caused by the increase in vegetation cover of green land is the possible cause of the decrease in light radiation.

The Shifting of the NTLGC
In this paper, the NTLGC model was introduced to reveal the differences in the direction of expansion between cities during the integration process ( Figure 4). For the integrated city of Guangfo, the geometric center of the boundary zone is located on the side of Foshan, while the center for the NTLGC is located on the side of Guangzhou, which indicates that Guangzhou is occupying the dominant position for the process of integration. During the study period, the NTLGC's shift in direction was "southeast-northwest", and the offset distance increased from 545 m (2013a-2015a) to 695 m (2015a-2017a), reflecting the fact that the NTLGC kept moving towards the Foshan side and the integration process advanced gradually. For the Guanshen city, both the geometric center and the NTLGC were located on the side of Shenzhen. The offset distance of the NTLGC increased from 1699 m (2013a-2015a) to 2086.33 m (2015a-2017a), showing a long-span and incremental characteristic. For Husu city, the geometric center was roughly located on the boundary between Shanghai and Suzhou, while the NTLGC was located in Suzhou. The shift trajectory of NTLGC in Husu presented a Land 2020, 9, 378 9 of 15 "V" structure. Although the gravity center shift direction has changed dramatically, in fact, since the shift distances are kept at about 200 m, the boundary zone's light intensity of Husu was found to be stable during the research period. Compared with the gathering characteristics of the NTLGC in Husu boundary zone, the shift of the NTLGC in Xixian boundary zone was more significant, and the shift direction was maintained as being "northeast-southwest".
comparison shows that the land use types in the light attenuation area are mainly vegetation and water bodies, while the light increase area is mainly industrial land (industrial parks and logistics parks), transportation land, and some commercial land. Similarly, Guanshen's light increase area is mainly located in industrial parks, and its light attenuation area is reservoirs, some urban villages and gardens, etc.; Husu's light attenuation zone mainly includes water bodies and green land, while the light increase plots mainly are new residential communities, industrial parks, and transportation facilities; Xixian's light attenuation area mainly includes land to be developed and green land within the city, while the light increase plots are mainly used for industrial and traffic facilities. Generally, in the boundary areas of integrated cities, logistics parks and industrial parks, driven by location and land prices, continue to attract concentrations of people and goods, and the intensity of human activity then increases significantly. Conversely, for cultivated and green land within a border zone, a decreasing trend in light intensity is observed, and we speculated that the shading effect caused by the increase in vegetation cover of green land is the possible cause of the decrease in light radiation.

Spatial Regression Analysis of Light Intensity Change
Since changes of light intensity in the integrated cite boundary zones are mainly caused by the transformation of construction land and vegetation coverage, we selected the determination construction land and vegetation coverage as the main research objective in the correlation analysis of light intensity and urban surface structure characteristics. We used Peng et al.'s [43] method to construct a grid in the study area as the research unit, and obtained the average values of NDVI (annual synthesis), NDBI (annual synthesis), and VCM of each grid unit.
while the NTLGC was located in Suzhou. The shift trajectory of NTLGC in Husu presented a "V" structure. Although the gravity center shift direction has changed dramatically, in fact, since the shift distances are kept at about 200 m, the boundary zone's light intensity of Husu was found to be stable during the research period. Compared with the gathering characteristics of the NTLGC in Husu boundary zone, the shift of the NTLGC in Xixian boundary zone was more significant, and the shift direction was maintained as being "northeast-southwest".

Spatial Regression Analysis of Light Intensity Change
Since changes of light intensity in the integrated cite boundary zones are mainly caused by the transformation of construction land and vegetation coverage, we selected the determination construction land and vegetation coverage as the main research objective in the correlation analysis  Table 4 shows the OLS regression results of nighttime light intensity with NDBI and NDVI in four groups of integrated cites (5 km). In general, the effects of buildings and vegetation on nighttime lighting intensity varies significantly in different time and cities. From the comparison of four groups of integrated cities, OLI regression results show two most obvious characteristics: (1) the goodness of fit between NDBI and VCM is higher than that between NDVI and VCM, which reflects that construction land occupies a greater explanatory power in light distribution; (2) based on the slope of the linear fitting equation, there is a negative correlation between NDVI and VCM, while NDVI and VCM show a positive correlation. Specifically, from the results of the linear regression between NDBI and VCM, the goodness of fit and regression coefficients of the Guangfo and Guanshen boundary zones are higher than those of Husu and Xixian, which implies that the newly constructed building plots in the Guangfo and Guanshen border area can provide higher light intensity. Moreover, the difference in the function parameters in 2013 and 2017 reflects that the expansion of regional built-up plots in the Guangfo and Guanshen boundary zones have been accompanied by rapid growth in regional human activity intensity during the process of urbanization, while Husu and Xixian show a declining trend. We speculate that the reason for this difference is due to the early start of the integrated development of Guangfo and Guanshen, which caused more population increases and commercial activities to occur along with the expansion of construction land during their development. In order to reveal the spatial instability characteristics of the nighttime light intensity distribution due to changes in buildings and vegetation cover, we performed GWR analysis on the above four groups of urban integrated boundary zones. Compared to OLS, the fit goodness of GWR that considers spatial instability was found to be significantly improved, and the R 2 of the four integrated cities in 2013 and 2017 were 0.837 and 0.805 (Guangfo), 0.753 and 0.731 (Guanshen), 0.620 and 0.660 (Husu), and 0.484 and 0.212 (Xixian), respectively.
In addition, GWR regression analysis also provided the local R 2 and coefficients of independent variables. Figure 5 shows the distribution of local R 2 and NDBI's regression coefficients for the regression analysis of the four groups of integrated cities in 2017. The results showed that the local R 2 of the integrated cities gave evidence of obvious spatial differentiation, including that: (1) the high values of the local R 2 of the Guangfo were concentrated in the south and north; (2) the high values of the local R 2 of the Guanshen were concentrated in the middle of the boundary zone; (3) the local R 2 of the Husu showed obvious "segmented" features, decreasing from the central region to the north and south; and (4) the local R 2 of the Xixian was also segmented, and decreasing from east to west. In addition, as noted earlier, the expansion of building land is the main contributor to regional changes in light levels. Thus, the regression coefficient for building land represents the change in lighting per unit of construction area, which reflects the level of regional integration. In Figure 5, the high values of NDBI's coefficients of Guangfo can be seen to be mainly distributed in the west of Baiyun District, the southeast of Liwan District, and the west of Panyu District, while the low values are mainly concentrated in the east side of Nanhai District and the middle part of Liwan District. For Guanshen, high value coefficients are distributed in the west side of the boundary zone, including Changan Town and Baoan District, while the low values are concentrated in the middle part. The high value coefficients for Husu are mainly located in Qingpu District and Kunshan City, while the low values are mainly located in the southern part of Taicang City and the western part of Jiading District. The distribution of the coefficients for Xixian shows a clear east-west difference, with the coefficients being higher in the east than in the west.
high value coefficients are distributed in the west side of the boundary zone, including Changan Town and Baoan District, while the low values are concentrated in the middle part. The high value coefficients for Husu are mainly located in Qingpu District and Kunshan City, while the low values are mainly located in the southern part of Taicang City and the western part of Jiading District. The distribution of the coefficients for Xixian shows a clear east-west difference, with the coefficients being higher in the east than in the west.

Discussion
With the continuous development of urbanization in China and the rapid expansion of cities to surrounding areas, many built-up areas of neighboring cities are gradually becoming connected in regional space. The increasing demand for coordinated development among cities promotes the trend

Discussion
With the continuous development of urbanization in China and the rapid expansion of cities to surrounding areas, many built-up areas of neighboring cities are gradually becoming connected in regional space. The increasing demand for coordinated development among cities promotes the trend of city integration. As an important mode of regional cooperation, city integration is conducive to narrowing the gap between cities and promoting coordinated development within the region. During the process of urban integration, due to China's distinct regional characteristics and significant differences in socio-economic development levels, the process of urbanization presents a variety of characteristics and stages. Therefore, this paper selected four groups of integrated cities with different backgrounds to carry out monitoring and analysis, which helped to reveal the differences in the process of urban integration. The process of city integration is a rather complex dynamic process, which is reflected in many areas such as society, the economy, transportation, the population, and so on. However, given the advantages of remote sensing in monitoring time-series surface changes, this paper focuses on changes related to surface activity and land use types in integration process monitoring, which can be extracted from remotely sensed images.
The remote sensing-based monitoring results of four groups of representative integrated cities under different natural and social-economic conditions can reflect some characteristics and differences of Chinese cities during the process of integration. Generally, the development level of Guangfo and Guanshen in the integration process are significantly higher than that of Husu and Xixian according to the land-use changes and nighttime light intensity fluctuation in the boundary zone. The integration development level of Guangfo is the most mature and thorough among the four groups of cities. Guangzhou and Foshan are adjacent to each other and belong to the birthplace of Guangfu culture, which has an important history. Since the reform and opening-up, the forces of industrialization and urbanization have accelerated the mobility of resource elements, prompting the continuous development of Fangcun, Huangqi, Yanbu, etc. in the boundary zone, and the city boundary has become increasingly blurred. For the integration process of Guanshen, the most prominent feature is the integration heterogeneity in the East-West direction. In recent years, with the rapid development of Shenzhen, the high land price and the shortage of land supply have accelerated the industrial spillover effect in Shenzhen. With the advantages of geography, industrial support, and urban environment, Dongguan has gradually become the first choice for industrial transfer in Shenzhen. Compared with Guangfo and Guanshen, Husu in the Yangtze River Delta region has a lower level of integrated development than the Pearl River Delta region. As the economic center of China, Shanghai, located in Land 2020, 9,378 13 of 15 the Yangtze River Delta, has a much larger economy than the neighboring city of Suzhou. Since the city area of Shanghai is 3.2 times that of Shenzhen, the shortage of land supply for urban development is not as prominent as in Shenzhen. Therefore, in the process of Husu integration, there exists the problem of inadequate initiative and enthusiasm, which is reflected in the changes of NTLGC in the boundary zone, showing different migration directions and the smallest migration range. For the Xixian, the changes in land-use, light intensity, and the shift of NTLGC in Xixian reflect that integration is mainly concentrated in the middle of the boundary zone. In addition, although China's economy and urbanization have made world-renowned achievements in recent decades, the spatial unevenness of development remains prominent. The eastern cities of mainland China, especially the coastal cities, are clearly superior to the inland cities in terms of location, policies, and environmental conditions, leading to serious unbalanced development between "East" and "West". Therefore, as an inland city, Xixian's economic development lags behind that of developed coastal areas, but the establishment of Xi Xian New District (located between Xi'an and Xianyang built-up areas) has provided new opportunities for its integrated development.

Conclusions
At present, research on city integration is mainly focused on describing the characteristics, formation mechanism, and dynamic mechanism of integration, while relatively little research has been conducted on integration based on remote sensing approach. Therefore, under the background of rapid urbanization in China, this paper chooses four of most representative integrated cities as the research sample, taking the boundary zone as the study area, which is the area with the most drastic surface changes in the integration process. Through the monitoring of changes in land use and nighttime light intensity in the boundary zones over the past five years, the differences in the integrated development of the above cities were revealed. The results show that the process of city integration in China is not only unbalanced in the level of development, but also has significant spatial heterogeneity in its internal development. Specifically, the results show that: (1) Guangfo has the highest level of integration and urban expansion is mainly concentrated in the south-central part of the boundary area; (2) Guanshen's level of integration is second to Guangfo's and is mainly concentrated in the west; (3) HuSu's integration is still in the initial stage and its increase in light intensity lags behind its expansion of building land during the study period; (4) although the light intensity and building land area increased significantly during the study period in Xixian, the overall development level of Xixian still lagged behind coastal cities due to the restriction of geographical location. Considering the aforementioned characteristics of integrated city development, planning decision makers should pay more attention to the rational allocation of resources and ensure that population transfer and infrastructure development are synchronized, which can promote the growth of regional dynamism.
It is worth noting that there are limitations in this study. Due to the limitations of nighttime light data, especially the defects of DMSP/OLS data and its difference with SNPP/VIIRS, both of which limited the width of the time window we studied, it was impossible to analyze the process of urbanization characteristics from a longer time span. Moreover, since policy factors are difficult to quantify, we focused only on the effects of human activities on the land surface during integration and ignored the endogenous factors such as policy and population migration. Therefore, in the following research, we will further improve several points: (1) how to integrate DMSP/OLS and SNPP/VIIRS data to achieve a longer time span of integration process research; (2) develop new data sources, such as POI (Point of Interest) data and land price data, to enhance quantitative research on the impact of policy and demographic changes behind the integration process.