Urban Sprawl and Changes in Land-Use Efﬁciency in the Beijing–Tianjin–Hebei Region, China from 2000 to 2020: A Spatiotemporal Analysis Using Earth Observation Data

: Sustainable development in urban areas is at the core of the implementation of the UN 2030 Agenda and the Sustainable Development Goals (SDG). Analysis of SDG indicator 11.3.1—Land-use efﬁciency based on functional urban boundaries—provides a globally harmonized avenue for tracking changes in urban settlements in different areas. In this study, a methodology was developed to map built-up areas using time-series of Landsat imagery on the Google Earth Engine cloud platform. By fusing the mapping results with four available land-cover products—GlobeLand30, GHS-Built, GAIA and GLC_FCS-2020—a new built-up area product (BTH_BU) was generated for the Beijing–Tianjin–Hebei (BTH) region, China for the time period 2000–2020. Using the BTH_BU product, functional urban boundaries were created, and changes in the size of the urban areas and their form were analyzed for the 13 cities in the BTH region from 2000 to 2020. Finally, the spatiotemporal dynamics of SDG 11.3.1 indicators were analyzed for these cities. The results showed that the urban built-up area could be extracted effectively using the BTH_BU method, giving an overall accuracy and kappa coefﬁcient of 0.93 and 0.85, respectively. The overall ratio of the land consumption rate to population growth rate (LCRPGR) in the BTH region ﬂuctuated from 1.142 in 2000–2005 to 0.946 in 2005–2010, 2.232 in 2010–2015 and 1.538 in 2015–2020. Diverged changing trends of LCRPGR values in cities with different population sizes in the study area. Apart from the megacities of Beijing and Tianjin, after 2010, the LCRPGR values were greater than 2 in all the cities in the region. The cities classed as either small or very small had the highest LCRPGR values; however, some of these cities, such as Chengde and Hengshui, experienced population loss in 2005–2010. To mitigate the negative impacts of low-density sprawl on environment and resources, local decision makers should optimize the utilization of land resources and improve land-use efﬁciency in cities, especially in the small cities in the BTH region.


Introduction
According to the United Nations, at present, more than 4 billion people live in urban areas worldwide, and this number continues to rise [1]. The global urbanization trend lead to excessive land use change and expansion of urban land [2]. The physical growth of urban areas is often disproportionate in relation to population growth which results in inefficient However, the accuracy of Earth observation products can vary between different regions, which can lead to uncertainties in evaluation results. Efforts have been made to measure LUE changes using satellite data at the local city scale. Using Landsat and SPOT images from 1996, 2001-2002 and 2011, the SDG 11.3.1 indicator was evaluated for four South African cities, and it was demonstrated that the population growth rate and rate of spatial expansion of small and second-tier cities in Africa was faster than that of large cities [20]. Ghazaryan et al. developed an automated classification approach and monitored urban sprawl and densification processes in western Germany using SDG Indicator 11.3.1 [21]. Spatiotemporal changes in urban land-use efficiency were analyzed and urban land lease policy gaps were emphasized to help ensure sustainable urban growth in fastgrowing cities in Ethiopia [22].
Since the reform and opening up which began in 1978, China's urbanization rate has increased from 17.9% in 1978 to 60.6% in 2019 [23]. This urbanization process is constantly accelerating and greatly promotes the nation' s social and economic development. One results of the country's industrialization and urbanization over the past four decades has been the gradual formation of large urban agglomerations or megalopolises [24,25]. The Beijing-Tianjin-Hebei (BTH) region is the largest such urban agglomeration in northern China. Monitoring and analyzing the dynamic changes in urban space is necessary to understand the urban development processes in this region. Therefore, the objectives of this study include: (1) to develop an approach to extract urban built-up areas by using Earth observation data; (2) to analyze urban expansion processes based on remote sensing products of built-up areas; and (3) to assess the spatial and temporal changes in urban landuse efficiency in the cities of the Beijing-Tianjin-Hebei urban agglomeration, China.

Beijing-Tianjin-Hebei Region
The Beijing-Tianjin-Hebei (BTH) region lies between the latitudes 36°03′ N and 42°40′ N and the longitudes 113°27′ E and 119°50′ E in northern China ( Figure 1). It consists of Beijing and Tianjin municipalities, and 11 prefecture-level cities in Hebei Province (Zhangjiakou, Chengde, Qinhuangdao, Tangshan, Cangzhou, Hengshui, Langfang, Baoding, Shijiazhuang, Xingtai and Handan). Among them, Beijing and Tianjin are megacities with a population of more than 10 million. The area of this region comprises 218,000 km², accounting for 2.3% of China's land area. The total population of the study area in 2020 is 113 million, accounting for 7.8% of China's total population.  Since 2014, the coordinated development of the BTH region has become a major national strategic policy of China. As the center of economic development in the north of the country, the BTH region plays a significant role in China's political and economic development [8]. To realize the coordinated development of the region and build a worldclass urban agglomeration, the urbanization of the region should be developed in a planned and sustainable way, which requires timely monitoring of the processes and status of urban growth.

Datasets
The datasets used in this study included Landsat satellite images, topographic data, population data and administrative boundaries. Details of these datasets are described in the following sections.

Landsat Imagery
The Google Earth Engine (GEE), which can process satellite and other Earth observation data, is a cloud-based computing platform [26]. The platform stores nearly 40 years of publicly available global remote sensing images, including Landsat, Sentinel, MODIS and DMSP/OLS satellite imagery [27]. In this study, Landsat 5 tier-1 surface reflectance data from 2000, 2005 and 2010 and Landsat 8 data from 2015 and 2020 provided by the Google Earth Engine were used as the main data source for built-up area extraction. These datasets were atmospherically corrected using the GEE platform; this processing included cloud, shadow, water and snow masks.

Built-Up Area Products
Four urban built-up land or land-cover products were used in this study, including GlobeLand30, GAIA, GHS-Built and GLC_FCS-2020. These all have a spatial resolution of 30 m. The characteristics of these datasets are shown in Table 1.
The GAIA product is a global high-spatial-resolution (30 m) annual artificial surface dynamic data product (1985-2018) released in 2019. GAIA is produced using the Google Earth Engine platform together with long-time series of Landsat images (nearly 1.5 million scenes) as the main data source and auxiliary data such as night-light data and Sentinel-1 SAR data [28]. The GAIA dataset reveals the different rates of change in artificial impervious surface cover in major countries and regions of the world.
The Global Human Settlement Layer (GHSL) was developed by the Joint Research Centre (JRC) and contains spatial information about human settlements at a global scale [29]. The GHSL uses the core method of symbolic machine learning (SML) for data extraction and uses Landsat images from the past 40 years as the main data source. The GHSL is an open and free dataset. The GHS-Built product was used in this study. This product consists of multi-temporal data extracted from Landsat imagery (consisting of the GLS1975, GLS1990 and GLS2000 datasets and ad-hoc Landsat 8 data for 2013/2014). There are three types of land use classified in this product: built-up areas, non-built-up areas and water bodies.
The GLC_FCS30-2020 data is a global 30 m land-cover product with a fine classification system for the year 2020 [30]. The GLC_FCS30-2020 dataset is produced based on timeseries of Landsat surface reflectance data (from 2019-2020), Sentinel-1 SAR data, DEM data, and auxiliary datasets. These data reflect the land cover types found globally (except for Antarctica) at a spatial resolution of 30 m in 2020.
The GlobeLand30 datasets from 2000, 2010 and 2020 were used in this study. These data include 10 land-cover types such as artificial land and bare land and cover the land within the range 80 • S-80 • N. The GlobeLand30 data were mainly classified using Landsat TM\ETM+ images and China Environmental Disaster Reduction Satellite (HJ-1) data based on a comprehensive pixel-object-knowledge classification method [31].
The data set used in this study are built-up area mapping results from these products. In the GAIA, GHS-Built, GLC_FCS30-2020, and GlobeLand30 products, the built-up areas are defined as impervious, built-up, impervious surfaces, and artificial surfaces, respectively. Hence, the four products are reclassified into two land cover types: non-built-up areas and built-up areas. The non-built-up area includes all other land cover types except built-up area, such as water, cropland, and forest. The specific reclassification schemes for each land cover product are shown in Table 2. In order to fuse theses datasets with the mapping results we derived using Google Earth Engine, the four built-up area products were re-projected to the same geo-reference system consisting of the WGS-84 coordinate system and UTM projection and resampled to a spatial resolution of 30 m.

Population Data
The WorldPop project, which is funded by the Bill & Melinda Gates Foundation, aims to offer spatial population datasets for Central and South America, Africa and Asia to support development, disaster response and health applications. The WorldPop data constitute an open-access gridded population data set (https://www.worldpop.org/ (10 April 2021)), produced by using a combination of machine learning and asymmetric modeling [32]. In this study, the population datasets for China for 2000, 2005, 2010, 2015 and 2020, which have a spatial resolution of 100 m and are based on the WGS84 geographic coordinate system, were used. The WorldPop population totals for different countries have been adjusted to match the corresponding official United Nations population estimates obtained from the Population Division of the Department of Economic and Social Affairs of the United Nations Secretariat.

Ancillary Datasets
The SRTM (Shuttle Radar Topography Mission) DEM is a digital elevation terrain model published by NASA and the National Bureau of Surveying and Mapping (NIMA). It has a spatial resolution of 30 m [33] and can be accessed directly in GEE. The administrative boundaries, including municipal district boundaries, provincial boundaries and municipal boundaries, can be obtained from National Geomatics Center of China (http://www.ngcc. cn/ngcc// (6 March 2021)).

Built-Up Area Extraction
Three spectral indices-the normalized difference built-up index (NDBI), normalized difference vegetation index (NDVI) and normalized difference water index (NDWI)were calculated from time-series of Landsat images [34][35][36]. Annual composite images were created from the optical, near-infrared, NDVI, NDBI and NDWI bands using the 'reduce' method by provided by GEE. The reduce method reduces a collection of images to an image by performing operation such as summation, median, etc. pixel by pixel. Four reduce functions-min, max, mean and median-were applied to improve the classification accuracy [37]. Sample points were selected by visual interpretation using the GEE platform. A stratified random sampling method was implemented with mapped land cover classes as strata [38]. The numbers of the sample points were determined based on the areal proportion of each land cover type. 80% of the sample points were randomly taken as training samples and 20% as validation sample points. Elevation and slope data calculated from the SRTM DEM were used as supplementary data; image classification was performed using the random forest algorithm(RFA) [39]. The number of classification trees in the random forest classifier is 500, and the remaining parameters are default. The confusion matrix between the validation samples and the classified products was calculated and the overall accuracy (OA) and kappa coefficients were obtained [40]. For the dichotomy of builtup area products (built-up area and non-built-up area) in this study, OA is the accuracy of the overall pixel being correctly classified, which can directly reflect the proportion of the correct classification. Kappa coefficient is a statistical index that can measure the balance of the correct classification of different categories. Omission error refers to the pixels belonging to the built-up area but classified as non-built-up area. Commission error refers to the pixels that is not built-up but is classified as built-up area. We tested different strategies for splitting the training and validation samples. If the classification result meets the requirement of OA > 0.85 and kappa > 0.8 [41,42], it was exported. The RFA was run iteratively and classification result with the highest OA was taken as the final output. The resolution of classification result was 30 m, and the spatial reference system was WGS-84 coordinate system and UTM projection.
The voting method is one of the most commonly used multi-classifier combination methods [43,44] in which the category with the most votes is taken as the final classification result. When multiple categories have the same number of votes, one of them is randomly selected as the final result. Based on the existing products for classifying urban built-up land and the classification results obtained from the RFA, a multi-temporal urban built-up area dataset was generated using the majority voting method. After the voting method had been applied, a temporal filtering method was used to ensure that the amount of built-up land in the product exhibited a growing trend [45]. A majority filter with a 3 × 3 window size was used to remove isolated pixels to obtain the final built-up area products with 30 m resolution, BTH_BU 2000, BTH_BU 2005, BTH_BU 2010, BTH_BU 2015 and BTH_BU 2020 [36]. The data-processing workflow is shown in Figure 2.
tested different strategies for splitting the training and validation samples. If the classification result meets the requirement of OA > 0.85 and kappa > 0.8 [41,42], it was exported. The RFA was run iteratively and classification result with the highest OA was taken as the final output. The resolution of classification result was 30 m, and the spatial reference system was WGS-84 coordinate system and UTM projection.
The voting method is one of the most commonly used multi-classifier combination methods [43,44] in which the category with the most votes is taken as the final classification result. When multiple categories have the same number of votes, one of them is randomly selected as the final result. Based on the existing products for classifying urban built-up land and the classification results obtained from the RFA, a multi-temporal urban built-up area dataset was generated using the majority voting method. After the voting method had been applied, a temporal filtering method was used to ensure that the amount of built-up land in the product exhibited a growing trend [45]. A majority filter with a 3 × 3 window size was used to remove isolated pixels to obtain the final built-up area products with 30 m resolution, BTH_BU 2000, BTH_BU 2005, BTH_BU 2010, BTH_BU 2015 and BTH_BU 2020 [36]. The data-processing workflow is shown in Figure 2.

Accuracy Assessment
In order to confirm the feasibility of the BTH_BU product for urban expansion analysis across time, it is necessary to evaluate its thematic accuracy. In this study, confusion

Accuracy Assessment
In order to confirm the feasibility of the BTH_BU product for urban expansion analysis across time, it is necessary to evaluate its thematic accuracy. In this study, confusion matrix was used to evaluate and compare the accuracy of the built-up area products [46][47][48]. The comparative analysis can also provide useful information for applications of these products in urban expansion and other fields. The confusion matrix compares the consistency of the land cover categories in the reference data and the data to be verified at a specific location, and then creates a two-dimensional table that compares the two (i.e., a confusion matrix). Based on this matrix, several different indices, including the overall accuracy (OA), kappa coefficient, omission error (OE) and misclassification error (CE), which can measure the accuracy of the product, can be calculated [49]. The equations used to calculate these indices are: Remote Sens. 2021, 13, 2850 where TP = true positive, TN = true negative, FP = false positive, FN = false negative and P o = OA. In this study, a true positive result indicated that the built-up area pixels had been correctly extracted. A true negative indicated that image pixels considered to be non-built-up areas had been classified as not built-up areas. A false positive meant that image pixels that did not correspond to built-up areas had mistakenly been extracted, and a false negative meant that image pixels that did correspond to built-up areas had not been extracted.
A total of 300 points including 150 classed as built-up land and 150 classed as non-builtup land were randomly created for each city using ArcGIS 10.5 software (Esri, Redlands, CA, USA). High-resolution satellite images from Google Earth obtained in the same year as the Landsat data were used to identify the land-cover types corresponding to these random points manually. An accuracy evaluation database was then built and used to calculate the error matrix and accuracy assessment indicators.

Functional City Boundary
Determining the extent of the built-up area is a prerequisite for measuring and comparing SDG 11.3.1. Instead of the administratively defined boundaries, the dynamic and functional urban boundary proposed by UN-Habitat-the actual area within which urban growth occurs over a defined period of time-was adopted in our study [10]. The functional urban boundary was obtained using several steps (Figure 3). The BTH_BU products are reclassified based on the density of built-up area per square kilometer. When the built-up area density is greater than 50%, it is defined as urban. If the built-up area density is greater than 0.25 and less than 0.5, it is defined as suburban. If the built-up area density is less than 0.25, it is considered as rural. The urban and suburban pixels were then used to determine the fringe open spaces, which are defined as areas within 100 m of urban and suburban pixels. A polygon was created by merging urban, sub-urban and the fringe open space pixels. Finally, buffer that extended the area of the polygon by 25% was created around each feature and the maximum extent of the buffer area was taken to be the functional urban boundary [7].

Change in Urban Built-Up Area
After obtaining the functional urban boundaries of the 13 cities in the study area in 2000, 2005, 2010, 2015 and 2020,the annual expansion area (AEA, km 2 ) was used to quantify the amount and rate of urban expansion in our study [17]. The AE is defined as follows:

Change in Urban Built-Up Area
After obtaining the functional urban boundaries of the 13 cities in the study area in 2000, 2005, 2010, 2015 and 2020, the annual expansion area (AEA, km 2 ) was used to quantify the amount and rate of urban expansion in our study [17]. The AE is defined as follows: where A end represents the urban built-up area at the end of a certain period, A start is the built-up area at the beginning of a period and d is the number of years between the start and end of the period. The AEA can directly measure the annual change of built-up area and reveal the difference of the expansion amount of a city in different stages.

Change in Urban Form
The changes in urban form were analyzed by calculating compactness and aggregation indices which are spatial metrics. The compactness (C) of a city's boundaries is an indicator that can reflect the form of the city's space. It is defined as: where C is the compactness of the city-that is, the compactness of the urban functional boundary. A is the area of the city and P denotes the length of the city perimeter. The greater the compactness value, the more compact the city is. A circle is the most compact shape: the compactness of a space in which each part is highly compressed is 1. In contrast, less compact shapes tend to be narrow or long. Also, the more complex a shape, the smaller its compactness. The aggregation index (AI) of a city is defined as: where AI denotes the degree of aggregation and g ii denotes the number of adjacent similar patches of the same landscape type. The AI is calculated based on the length of the common boundary between pixels corresponding to the same type of patch. When there is no common boundary between any of these pixels, the degree of aggregation is the lowest; when the common boundary between all pixels corresponding to the same type of patch reaches its maximum value, the aggregation index is the highest. In this study, the compactness and aggregation index for each built-up area, as defined by the functional urban boundaries, were calculated for each time period.

Derivation of LCR, PGR and LCRPGR
The SDG 11.3.1 indicators give the ratio of the land consumption rate (LCR) to the population growth rate (PGR). The PGR refers to the rate of population change in a city/urban area over a given period (usually a year). The LCR is defined as the rate of appropriation of land by urban development. The ratio of land consumption rate to population growth rate(LCRPGR) is calculated as the ratio of LCR to PGR. SDG 11.3.1 indicators are calculated using the following formulae: where Pop t is the total population of the urban area in the past/initial year, Pop (t+n) is the total population in the current/final year, Urb t is the total extent of the urban area (in km 2 ) for the past/initial year, Urb (t+n) is the sum of the urban areas in the current/final year and y is the number of years between the two measurement periods. Ideally, the ratio of land consumption to the population growth rate would be equal to 1, implying that the amount of urban land is increasing at the same rate as the population. Based on LCRPGR values, cities can be categorized into 5 types. These are shown in Table 3. Table 3. City classification based on the ratio of land consumption rate to population growth rate values.

LCRPGR Value Meaning
LCRPGR < −1 the rate of population decline is greater than the rate of built-up area expansion −1 < LCRPGR < 0 the rate of population decline is less than the rate of built-up area expansion 0 < LCRPGR < 1 the rate of population growth is greater than the rate of built-up area expansion 1 < LCRPGR < 2 the rate of built-up area expansion is 1-2 times the rate of population growth LCRPGR > 2 the rate of built-up area expansion is greater than 2 times the rate of population growth

Accuracy Assessment
Accuracy evaluation indices were obtained for each of these built-up area products for each study period (Table 4). Among the built-up area products, BTH_BU products had the highest average OA and Kappa coefficients and the lowest average OE and CE (Table 4). It shows that the proposed method can effectively improve the extraction accuracy of built-up area.  Figure 4 shows the OA, kappa coefficient, OE and CE of the five remote sensing products for each city. The classification accuracy of BTH_BU products in 11 out of 13 cities was the highest. BTH_BU and GlobeLand30 have the smallest variability in classification accuracy, followed by GLC_FCS and GHS-Built. GAIA has the largest variability in classification accuracy between the different cities.  Figure 4 shows the OA, kappa coefficient, OE and CE of the five remote sensing products for each city. The classification accuracy of BTH_BU products in 11 out of 13 cities was the highest. BTH_BU and GlobeLand30 have the smallest variability in classification accuracy, followed by GLC_FCS and GHS-Built. GAIA has the largest variability in classification accuracy between the different cities.  Figure 5 shows the spatial pattern of urban expansion in the cities in the Beijing-Tianjin-Hebei region. Different types of urban spatial growth, infilling, edge-expansion, and leapfrogging, can be observed in these cities. The scattered built-up patches in suburban and rural areas was gradually connected with the central urban area, and thus forming a larger urban core in large cities.

Changes in Urban Built-Up Area
The annual expansion area (AEA) of urban built-up area showed obvious regional differences ( Figure 6). From 2010 to 2015, the AEA of the Beijing-Tianjin-Hebei built-up area is much higher than that in other periods. The dramatic changes in these values may indicate that a fusion between suburban patches and the core urban district or explosive growth resulting from the sprawl of built-up areas occurred during this period.  Figure 5 shows the spatial pattern of urban expansion in the cities in the Beijing-Tianjin-Hebei region. Different types of urban spatial growth, infilling, edge-expansion, and leapfrogging, can be observed in these cities. The scattered built-up patches in suburban and rural areas was gradually connected with the central urban area, and thus forming a larger urban core in large cities.  The annual expansion area (AEA) of urban built-up area showed obvious regional differences ( Figure 6). From 2010 to 2015, the AEA of the Beijing-Tianjin-Hebei built-up area is much higher than that in other periods. The dramatic changes in these values may indicate that a fusion between suburban patches and the core urban district or explosive growth resulting from the sprawl of built-up areas occurred during this period.

Changes in Urban Form
The changes in the compactness (C) and aggregation index (AI) for each city are shown in Figure 7. The compactness of all the cities in the Beijing-Tianjin-Hebei region shows a downward trend, while for the aggregation index (AI) there is an overall upward trend from 2000 to 2020 (Figure 7). The increasing trend in AI and the decreasing trend in C for all of the cities indicate that the built-up areas in all of these cities are gradually

Changes in Urban Form
The changes in the compactness (C) and aggregation index (AI) for each city are shown in Figure 7. The compactness of all the cities in the Beijing-Tianjin-Hebei region shows a downward trend, while for the aggregation index (AI) there is an overall upward trend from 2000 to 2020 (Figure 7). The increasing trend in AI and the decreasing trend in C for all of the cities indicate that the built-up areas in all of these cities are gradually amalgamating and the spatial forms of the functional urban areas are becoming increasingly complex [50]. amalgamating and the spatial forms of the functional urban areas are becoming increasingly complex [50].

Variations in LCR, PGR and LCRPGR
The process of urbanization in the Beijing-Tianjin-Hebei region from 2000 to 2020 was analyzed by calculating the LCR, PGR and LCRPGR values ( Table 5). The LCR, PGR and LCRPGR for the region all followed the same trend from 2000 to 2020: a decrease, followed by an increase, followed by another decrease. The decrease and increase in the LCRPGR during the period 2000-2015 agrees with the observed trend for all Chinese cities [51]. As a result of development in the early decades of the 21st century, there was also a consistent outward spread of the urban area in the Beijing-Tianjin-Hebei urban agglomeration. Before 2010, the rate of urban land growth was less than or similar to urbanization rate of the population. After Beijing hosted the Olympic Games in 2008, economic devel-  The process of urbanization in the Beijing-Tianjin-Hebei region from 2000 to 2020 was analyzed by calculating the LCR, PGR and LCRPGR values ( Table 5). The LCR, PGR and LCRPGR for the region all followed the same trend from 2000 to 2020: a decrease, followed by an increase, followed by another decrease. The decrease and increase in the LCRPGR during the period 2000-2015 agrees with the observed trend for all Chinese cities [51]. As a result of development in the early decades of the 21st century, there was also a consistent outward spread of the urban area in the Beijing-Tianjin-Hebei urban agglomeration. Before 2010, the rate of urban land growth was less than or similar to urbanization rate of the population. After Beijing hosted the Olympic Games in 2008, economic development led to rapid urban expansion in the Beijing-Tianjin-Hebei region. From 2010 to 2015, due to the acceleration in urban expansion and population growth, the LCR increased to 0.154, the PGR increased to 0.069 and the LCRPGR increased to 2.232. The urbanization rate of land was 2.232 times of the urbanization rate of population. In 2014, the National New Urbanization Plan (2014-2020) was issued by the national government and put forward a coordinated development strategy for the Beijing-Tianjin-Hebei region. As a guideline for China's urbanization, the plan emphasized the importance of achieving a more environmentally and sustainable form of urbanization. During the period 2015-2020, the rate of increase in urban land and population urbanization slowed down, the LCR decreased to 0.048, the PGR decreased to 0.031 and the LCRPGR decreased to 1.538. The temporal changes in the LCR, PGR and LCRPGR values of each city in the Beijing-Tianjin-Hebei region were analyzed. In the economically developed areas such as Beijing and Tianjin, the population growth and urban expansion were more rapid than in most cities in Hebei province. From 2000 to 2010 and from 2015 to 2020, although Beijing, as the capital of China, attracted a great inward flow of people, the rate of expansion of the built-up area was less than the population growth rate and the LCRPGR values were less than 1. Tianjin, the second largest city in the BTH region, had a higher LCRPGR value of between 1.0 and 1.5 (Figure 8). From 2010 to 2015, the built-up areas in all the cities in the Beijing-Tianjin-Hebei region expanded rapidly and their LCRPGR values were all greater than 1.5. The LCRPGR values for Langfang, Cangzhou, Baoding, Zhangjiakou and Tangshan have been increasing since 2000. By 2020, the LCRPGR values of these cities were all greater than 2, indicating that the cities had entered a phase of extensive development. In Qinhuangdao, the LCRPGR was greater than 2 during the period 2000-2020. Population growth in Hengshui has been relatively slow. The population growth here was negative from 2005 to 2010 and the absolute value of the LCRPGR from 2000 to 2020 was greater than 3. This indicates that a balance was maintained between the rates of urban expansion and population growth in Beijing and Tianjin. However, development in the cities in Hebei was unbalanced, with most of these cities underwent a phase of extensive development.
opment. In Qinhuangdao, the LCRPGR was greater than 2 during the period 2000-2020. Population growth in Hengshui has been relatively slow. The population growth here was negative from 2005 to 2010 and the absolute value of the LCRPGR from 2000 to 2020 was greater than 3. This indicates that a balance was maintained between the rates of urban expansion and population growth in Beijing and Tianjin. However, development in the cities in Hebei was unbalanced, with most of these cities underwent a phase of extensive development.

Differences in LCRPGR by City Type
We divided the 13 cities in the study into five categories based on population: very small cities, small cities, medium cities, large cities, and megacities ( Table 6). The LCRPGR values for cities with different population sizes are shown in Figure 9. Cities of different sizes undergone diverged urbanization processes during the study period ( Figure 9). The LCRPGR values of megacities and large cities in 2015-2020 were lower than that before 2015, while the LCRPGR value of medium-sized cities, small cities and very small cities in 2015-2020 is higher than that before 2015. The slowed down

Differences in LCRPGR by City Type
We divided the 13 cities in the study into five categories based on population: very small cities, small cities, medium cities, large cities, and megacities ( Table 6). The LCRPGR values for cities with different population sizes are shown in Figure 9.

Urban Growth Analysis
In this study, a new built-up area product (BTH_BU) was generated for the Beijing- Cities of different sizes undergone diverged urbanization processes during the study period ( Figure 9). The LCRPGR values of megacities and large cities in 2015-2020 were lower than that before 2015, while the LCRPGR value of medium-sized cities, small cities and very small cities in 2015-2020 is higher than that before 2015. The slowed down sprawling development of megacities and large cities since 2015 might be attributed to the implementation of the new urbanization policy. However, the land use efficiency remains very low in medium, small and very small cities after 2015.

Urban Growth Analysis
In this study, a new built-up area product (BTH_BU) was generated for the Beijing-Tianjin-Hebei (BTH) region by fusing the mapping results with four available land-cover products-GlobeLand30, GHS-Built, GAIA and GLC_FCS-2020. The accuracy of different built-up area or land cover products varies across cities. The accuracy of built-up area mapping mainly depends on the data sources, sample samples and classification methods employed. BTH_BU and GAIA rely mostly on the spectral features of Landsat data and a pixel-based classifier, whereas GHS-Built and GlobeLand30 also employ textural and contextual information about settlements that is obtained from remote sensing images. The GAIA and GLC_FCS-2020 uses both optical and Sentinel-1 SAR data. Our results showed that the integrative use of optical and SAR data improved the identification of buildings in less-dense urban areas. The sensitivity of SAR backscatter to building structures can help distinguish bare lands or sparse vegetated lands from built-up areas which are often confused in multispectral images [52]. By fusing the urban mapping results with four available products with a majority voting method, the BTH_BU takes advantage of different data products, and thus improved the accuracy of the final classification results.
Using the BTH_BU product and urban land use efficiency indicator, our results revealed diverged changing trends of LCRPGR values in cities with different population sizes in the study area. The extensive sprawling growth patterns in medium and small sized cities of this region suggest that economical and efficient use of land resources was neglected in urban planning and management strategies in these cities [17]. As the third largest urban agglomeration in China, the urban sprawl trajectory of cities with different sizes in Beijing-Tianjin-Hebei region can provide implications for other urban agglomerations [9]. Moreover, the urban mapping method and spatio-temporal urban sprawl analysis workflow using Earth observation data and geospatial technology can be applied to other cities and urban agglomerations.

Uncertainties and Limitations
Although the urban sprawl was characterized with the LCRPGR indicator effectively in our study, there are some limitations which should be considered. Multi-temporal Landsat imagery with 30 m spatial resolution was used as the main data source to extract urban land in this study, which limits the accuracy of built-up area extraction. The application of high resolution satellite data such as Sentinel-1/2, SPOT-5/6/7 and WorldView-2/3 may provide more detailed and accurate information on urban built-up areas [53,54]. For built-up area data fusion, different methods such as entropy weighting and performance weighting and the Dempster-Shafer theory can be applied in future studies [43]. When the LCRPGR value less than 0 is obtained, it should be interpreted with cautions. If PGR is less than 0 and LCR is greater than 0, the LCRPGR represents urban expansion and population loss. If PGR is greater than 0 and LCR is less than 0, it represents the decrease of built-up area and the increase of population in cities. Moreover, the indicators of SDG11.3.1 only use land and population as indicators, which may not fully characterize the urbanization process. Other environmental and economic indicators in SDGs can be used together with LCRPGR to perform a more comprehensive assessment of the urban development trajectories and trends.

Conclusions
In this study, a land-cover classification method was developed using a random forest classifier and Google Earth Engine cloud platform. The built-up area was extracted using time-series of Landsat imagery, a DEM and other ancillary data of the Beijing-Tianjin-Hebei region for the years 2000, 2005, 2010, 2015 and 2020. The data were also fused with four existing land-classification products-GlobeLand30, GHS-Built, GAIA and GLC_FCS-2020 to generate the built-up products BTH_BU 2000, BTH_BU 2005, BTH_BU 2010, BTH_BU 2015 and BTH_BU 2020. An accuracy assessment produced an overall accuracy of 0.93 and a kappa coefficient of 0.85 for the BTH_BU products. The built-up areas in all of these cities are aggregating and the spatial forms of the functional urban areas are becoming more complex. The spatiotemporal dynamics of the SDG 11.3.1 land-use efficiency indicators were monitored using BTH_BU and population data. The results showed that, for the BTH region overall, the LCRPGR values were close to 1 from 2000 to 2010 but rose to 2 or higher after 2010. Diverged changing trends of LCRPGR values in cities with different population sizes in the study area. The rates of land and population urbanization were found to be relatively balanced in the megacities of Beijing and Tianjin. Except for the megacities, the LCRPGR values were greater than 2 after 2010. The small and very small cities had the highest LCRPGR values after 2015; however, some of these cities, such as Chengde and Hengshui, experienced population loss. To mitigate the negative impacts of low-density sprawl on the environment and land resources, local decision makers should optimize the utilization of land resources and improve land-use efficiency in cities, especially in the small cities in the Beijing-Tianjin-Hebei region.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.