Expansion of Urban Impervious Surfaces in Xining City Based on GEE and Landsat Time Series Data

Urban expansion is often studied in large cities such as Beijing, Shanghai, and Guangzhou, while scant attention is paid to smaller cities such as Xining. However, Xining is the largest city on the Tibetan Plateau, and an important city in China’s “Belt and Road Initiative”. As its economy and society develops, Xining will play an increasingly important role in connecting the central and western regions. In order to quantify the impacts of rapid urbanization, it is extremely important to collect data on the time and space variations of impervious surfaces. As such, we collected Landsat long-term sequence data about Xining City from 1987–2019 using the random forest method, and then optimized the feature parameters to obtain the dataset. Our results demonstrated that the overall accuracy of land use classification in Xining city is 83.4% and that the urban impervious surface accuracy is 89.5%. Additionally, the overall accuracy improved by 2.4% after optimizing the characteristic parameters, while the urban impervious surface accuracy is 92.8%. In 27 of the 33 years we studied, the classification accuracy of impervious surfaces exceeded 90%. After correcting for the temporal consistency check, the accuracy of impervious surfaces improved by 2% compared to the original sequence. We analyzed the change of impervious surfaces in Xining based on the results of the final dataset and found that the impervious surface area of Xining increased from 55 km2 in 1987 to 334 km2 in 2019. Xining is a typical semi-open river valley city which shares spatial and temporal characteristics with other urban centers. The spatial and temporal characteristics of the expansion of urban spaces in the main urban area of Xining are obvious and are primarily spread around the central area toward tree branch shaped road, which help other cities located in river valleys better understand how urbanization progresses.


I. INTRODUCTION
Cities are home to most humans and economic activity, making them the focus of much academic research [1]. Urbanization involves changes in population, the economy, and land use [2]. According to a UN report, China and India will be the main sources of urban population growth in the next 30 years [3]. China's urbanization has been accelerating since the beginning of 1980s, with a large number of people relocating from rural areas to urban areas [4]- [8].
The associate editor coordinating the review of this manuscript and approving it for publication was Chaker Larabi .
Rapid economic growth and the gradual transformation of agricultural lands to urban cities has rapidly increased the population of China's cities [9]. This pace of urbanization greatly impacts the environment and climate on different scales: the urban heat island effect makes Xining significantly warmer than surrounding rural areas [10], [11], local and regional precipitation levels change, water quality deteriorates [12]- [14], increases in urban impervious surfaces induce water runoff, and the loss of agricultural land negatively impacts food production [15]- [17]. In order to quantify the impact of rapid urbanization, detailed information on the temporal and spatial changes of impervious surfaces is VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ important, making it necessary to study the expansion of urban impervious surfaces [18].
In the past few decades, remote sensing has been proven to be a reliable and effective method of monitoring changes to impervious surfaces [19]. Remote sensing data from multiscale platforms provides different angles for the monitoring of impervious surfaces [20]. Open access to remote sensing data provides new opportunities for collaboration [21]. Landsat Thematic Mapper / Enhanced Thematic Mapper / Operational Land Imager data (Landsat TM/ETM+/OLI ) in multispectral remote sensing images has been free since 2008. It is considered promising data for building a timespace model of impervious surfaces, due to its superior spatial resolution and spectral coverage and dataset stretching back to 1982 [22]- [24]. As the availability of high spatial resolution satellite images has increased, products providing a spatial resolution of 30m of several global and regional land cover products can be used directly [25]- [28]. However, none of them update frequently enough to monitor long-term changes in land cover. Additionally, there have previously been no quick and effective mapping technologies and tools, but in recent years, high-performance computing and more efficient methods of analysis have been integrated with powerful cloud computing resources [27], [29]. The emergence of these powerful cloud computing resources, such as Amazon Web Services, NASA Earth Exchange, Microsoft Azure, and Google Cloud, has conferred benefits on large caches of geographical data [30]. For example, Google Earth Engine (GEE) is an access-based cloud platform, which has a large number of satellite images and geospatial data sets, facilitating the development of algorithms and visualizations with reasonable processing times [31]- [33]. In addition to computing and storage capabilities, many well-known machine learning algorithms have also made significant advances [34].
Previous methods of analysis such as the decision tree model [35], [36], regression modeling [37], [38], and machine learning methods have been developed from single-phase Landsat satellite data to analyze maps of urban land cover and monitor associated changes [39], [40]. However, most of these studies only use Landsat data from limited periods [22], [41]- [46], which makes it difficult to analyze urban expansion dynamics in short periods of time. Many researchers currently use multi-temporal remote sensing images to analyze the spatial and temporal patterns of urban expansion. Change detection technology has also been applied to these multi-temporal images. For example, Li et al. used Landsat time-series data to map the average annual frequency of impervious surfaces in Beijing from 1984 to 2013 [1]. These researchers proposed a time consistency check algorithm, which included time filtering and logical reasoning, and subsequently resulted in a more accurate classification of impervious surfaces. Chai et al. proposed a new method to analyze urban land expansion, using Tianjin from 1990-2014 to analyze the spatio-temporal dynamics of urban expansion [18]. Shi et al. proposed an uncertainty-based spatial-temporal consistency (USTC) model to improve the accuracy of the classification of impervious surfaces across a long period of time, applying this model to Landsat images from 1987-2016 to obtain a map of the annual impervious surfaces of the city of Wuhan [24]. Song et al. used the estimated annual impervious surface coverage in Landsat time series data to obtain the scale of urban expansion in the Washington D.C.-Baltimore metropolitan area [47]. Time smoothing technology was used to remove noise in the time-series of urban land coverage. Fu et al. used Landsat time series images and an optimized Crulist tree model to analyze the ISP sequence [48]. They monitored and analyzed the changes of urban expansion and deurbanization in Guangzhou from 2000 to 2010 and determined trends of annual and long-term urban growth. Using long-term continuous annual average frequency captures more detailed information about rapid urbanization than using a single year or multiple time node images [49]. The existing TC or USTC algorithm addresses the spatial and temporal irrationality of impervious surface spots, however, this kind of algorithm is only used to modify and optimize the original results to obtain satisfactory results. Better results will be obtained if the classification accuracy improves in the early stages of classification and space-time optimization processing.
In this study, we performed research on urban expansion using the open cloud platform to obtain data and tools in Xining, a high altitude city with complex terrain. Specifically, we focused on the following three tasks: (1) Using Google Earth Engine to obtain high-quality images of long time series, method of sample data sets, and the choice of classifier; (2) Adding spectral information, texture, terrain, and climatic factors to the classification process. Using feature parameter optimization steps to establish feature optimization and time consistency algorithms to more accurately classify data and obtain a long-term sequence of high-quality surface data from Xining; (3) Analyze the development of impervious surfaces in Xining from 1987 to 2019.

A. STUDY AREA
Xining is the capital of Qinghai Province and is located between 36 • 12 27 ∼37 • 30 09 N and 100 • 47 56 ∼101 • 56 49 E in the Huangshui valley of the Tibetan Plateau. The Huangshui Valley runs east to west through the city, while the Beichuan River flows into the valley from the north, and the Nanchuan River flows into the valley from the south. It is a narrow city running east to west and is situated in a valley city [50], [51]. The highest and lowest altitudes are 4849m and 2102m, respectively, while the urban area has an average altitude of 2295m [52]. Xining's primary industries are agriculture and animal husbandry, and its culture is a fusion of Chinese culture and Tibetan culture. It is the gateway to the Tibetan Plateau, and is the political, economic, cultural, scientific, technological, transportation, and medical center of Qinghai Province. As shown in Figure 1, the scope of this study includes not only Xining   Municipal District, but also its surrounding parts. The Landsat uses a global reference system whose two-dimensional coordinates are identified by path and row. The study area boundary is located within one scene (path/row:132/035).

B. DATASETS 1) LANDSAT DATA
The data in this study comes from Landsat 5 TM, Land L and sat 7 ETM+, and Landsat 8 OLI images, all of which were obtained from the GEE platform. The data was processed into L1T products, which, after geometric and atmospheric corrections, are surface reflectance products. Figure 2 shows 377 scenes of Landsat images, whose cloud coverage is less than 10%. Two periods were considered in the imaging selection: the lush period of plant growth (June to October) and the withering period of plant growth (November to March of the following year). The withering period was used to calculate the Normalized Difference Vegetation Index (NDVI) of the area in this period and can be used as one a future characteristic parameter during classification.

1) TECHNICAL PROCESS
This study analyzes the annual changes in urban land usage over a long period of time using the gee platform. It is divided into three stages: the first stage uses Google Earth to collect samples of different years and select appropriate images; the second stage selects spectrum, texture, terrain, and climate features to identify their construction, and use the random forest method to classify Landsat data by year; the third stage performs post-classification processing and confirms the time consistency of 33 years of urban land series data, filters out noise caused by classification errors, obtains data related to annual land use and land cover, confirms the accuracy of our results, and performs space-time analysis on the dataset of annual impervious surfaces.

2) OTHER AUXILIARY DATA a: SAMPLE COLLECTION
The quality of samples directly determines the quality of the classification results, making it very important to collect high-quality and relatively stable samples. This study is a high-resolution study of long-term sequences, meaning the strategy for obtaining samples varied from year to year. From 1987-2000, the samples were selected from Landsat images by visual interpretation since there was no high-resolution image. Samples from 2001-2019 were selected from Google Earth Pro by manual visual interpretation. Based on the 2016 image, the two blue sample are type of urban. By 2015, the urban plot on the left becomes cropland, then delete the original urban plots and add cropland plots; in 2014, the urban plot on the right also becomes cropland, delete the original urban samples and add cultivated land samples. By analogy, according to the actual changes of the ground features, point by point comparison and year by year modifi-cation, so as to complete the establishment of all types of sample libraries in each year. This method ensured the stability and continuity of samples across the 33-year time period ( Figure 4). Finally, 80% of the samples were used as training samples and 20% were used as test samples [55]. The classification system of the study area was divided into 6 categories, based on the National Remote Sensing Monitoring Land Use/land cover classification system: city, cropland, grassland, forest, water body, and bare land. The collected training sample information is shown in Table 1.

b: METHOD OF FEATURE OPTIMIZATION
The selection of feature variables is an important step in remote sensing classification and can be improved by using a variety of feature variables and their combinations. Too many feature vectors in the model will also cause data redundancy, and not all features vectors all play a positive role in improving accuracy, so we must optimize the selection of feature vectors [56]- [59]. We optimized the steps used in this study according to the following: Step 1: The random function is used to perform 10 different spatial distributions on the training samples, and then the accuracy of the classification results of these 10 times is compared, and the spatial distribution with the highest accuracy is selected as the final sample point.
Step 2: Selection of texture features and optimal window: obtain texture features from GEE by GLCM (Gley-Level Co-ourrence Matrix) calculations, set texture feature window from 1-10 and conduct 10 experiments to find the best texture feature and the optimal window size; Step 3: Climate factor selection: analyze all climatic factors and select the one with the highest classification accuracy for subsequent classification; Step 4: Sorting selection: classification engineering demonstrated that the order of characteristic parameters impacted classification results. Considering the calculation resources and time cost, only the latest five characteristics were permutated, resulting in 120 different combinations. Finally, the combination with the highest classification accuracy was selected for the final classification parameters ( Table 2).

3) POST-CLASSIFICATION PROCESSING AND TEMPORAL CONSISTENCY CHECK
Using RF to classify remote sensing images from year to year, we identified fragments and ''pepper-salt phenomena'' in the results. Mode filtering in the ArcGIS software was used to perform fragmented speckle removal and post-classification merged clustering. We also obtained a long-time (33 years) series map of urban land classification. However, there could be classification errors in these results for the following reasons: there were low-quality Landsat images in earlier years, some inappropriate samples could have been selected in earlier years, and uncertainty during the classification process, such as ''same object with different spectrums'' and ''different objects have the same spectrum''. As such, we preprocessed the classification results. To test whether the obtained results conformed to the objective laws of urban development, we had to distinguish occurrences of change from classification errors across the long-term dataset. We obtained a more reliable time-series result by using the temporal consistency check method proposed by Dr. Li Xuecao and modifying it locally in order to improve the initial urban land classification results.
The process consists of two steps: temporal filtering and logical reasoning. The temporal filtering formula is shown in Figure 1 indicates a more confident label L i in the temporal sequence, whereas a lower value may reflect an erroneously classified pixel in that year. We set a threshold value of 0.5 to distinguish whether the current category needs to be converted. If it exceeds 0.5, the current category label belongs to the dominant category in the time domain, with the unchanged category; if it is less than 0.5, it will convert the current category label to the opposite. Through this process, individual errors can be corrected. The output of the model is a series of continuous cities (1) and non-cities (0).
In this formula, L i is the label of the target year and L j is its temporal neighborhood. Besides, Con() is an identifying function returning 1 when L j = L i , otherwise 0. Its main purpose to remove the isolated urban patches in the time dimension. For example, 1 is a city and 0 is a non-city. When [0 1 0] appears, the unreasonable situation of this isolated year is modified to[0 0 0], and it will also be modified to [1 1 1] for [1 0 1]. The purpose of this step is only for isolated urban patches. The following will be performed on a long term sequence logical judgment.
The logical judgment mainly includes a presupposition, that is, under normal circumstances, the development of cities is irreversible. Then the urban area in 2019 is the largest,  and there should be no urban distribution beyond 2019 in previous years. After that, there will be two cases. For a pixel time series, if 1987 is a city, and 2019 is also a city, then the entire time series from 1987 to 2019 will be modified to be as city. This situation is rare, the main corresponding actual situation is in the old part of the city, urban architecture has been in existence from 1987 to 2019. The second case, for a pixel time series, 1987 is not a city, but it is a city from a certain year (Xyear), then it is a city from Xyear→2019, and Xyear is a sudden change from non-city to city. This year accounted for most of the year, that is new urban construction.

1) ANALYSIS OF 33-YEARS PRELIMINARY CLASSIFICATION RESULTS
The overall classification accuracy in the article is obtained by calculating the confusion matrix from the verification sample. Many studies have demonstrated that adding texture features can improve classification accuracy [60], [61]. The green line refers to classification accuracy with additional parameters such as Normalized Difference Water Index (NDWI), Radio Vegetation Index (RVI), NDVI, Enhanced Vegetation Index (EVI), and Normalized Difference Building Index (NDBI), the red line refers to classification accuracy after selection using texture features, optimal windows, climatic factors, and feature parameters ( Figure 6). The red line is typically higher than the green line, meaning the accuracy of classification results following feature optimization is higher (by 2.4% higher on average) than that of classification results without optimization. The one exception is in 1991, when the accuracy of the latter was slightly higher than that of the former by 0.44%. In 2011, the difference between them was largest, reaching 8.98%. While determining the final feature elements in this article (see 2.3.2.2 for the process), there were cases when the overall classification accuracy of these schemes was almost the same. When this situation occurred, we gave priority to the classification accuracy of the single category of impervious surfaces in order to select the final feature parameter optimization scheme. The blue line in Figure 6 displays the overall accuracy of the city. The accuracy of other years exceeds 87%, except for 1989 and 1993, where the respective classification accuracy of impervious surfaces was 82.35% and 83.33% The accuracy of 27 years exceeded 90%. This provides a solid basis for subsequent processing and analysis..

2) ANALYSIS OF 33-YEARS PRELIMINARY CLASSIFICATION RESULTS
The accuracy of the modified classification results must be confirmed to test the effects of the temporal consistency algorithm. Figure 7 shows the overall accuracy and impervious surface accuracy results from 1987-2019, both before and after the time consistency treatment. Compared to the original classification results, applying time consistency slightly improved the overall accuracy. In Figure 7, the detection range of impervious surfaces following time consistency processing is 84%-96%, which is about 2% higher than the initial sequence. Therefore, the time consistency algorithm helps solve inconsistent long-term sequence classifications.  For example, during the initial classification of impervious water surfaces, some overestimates or underestimates will be corrected after applying the time consistency algorithm. The accuracy report has improved the accuracy of the classification of impervious surfaces.

B. DISPLAY OF CLASSIFICATION RESULTS
Comparing the results of classification with and without feature optimization (Figure 8), column A is the summer image synthesized by standard false color, column B is the preliminary classification result, and column C is the classification result after feature optimization. To compared the first row, the image of feature optimization show that the water area is increased and the spatial connectivity is better; the second row, the optimized results reduce the expansion of the urban area; the third row, misclassification of forest features have been supplemented.
Additional analysis of classification results following the temporal consistency check (Figure 9) demonstrated that deficiencies (such as data loss or image quality) in the  initial classification can be mostly corrected. For example, the overestimation of urban land use caused by the confusion of bare land in the initial image can be corrected ( Figure 9A), while underestimation, missing, or misclassification of land use caused by poor image quality can also be corrected ( Figure 9B). As such, time-series analysis can improve long-term serial urban remote sensing classification. VOLUME 8, 2020  Figure 10 displays the impervious surface map of the study area from 1987 to 2019. Terrain determines the spatial pattern and direction of urban expansion [62]. The most fundamental difference between river valley cities and typical cities is the river valley's restrictions on urban expansion [63]. Cheng et al. (2003) have researched the expansion of mountainous urban spaces, and found that rivers and mountains prevent cities from growing on a large scale [64]. Typically, their structural forms are star-shaped, fan-shaped, or irregularly clustered, while most other cities form a radial ring. Therefore, the direction, pattern, and land structure of urban expansion in river valleys are restricted by geographical geomorphology. The urban expansion model of Xining reflects a branch shaped spatial structure. The geographical location of Xining has significantly influenced the urban expansion of Xining along the river valley from 1987 to 2019. Figure 11 shows the land use conversion map, the annual impervious area change histogram, and the six types of land use change percentage pie charts for the study area from 1987 to 2019. The impermeable area of the study area kept increasing throughout the period. The land-use conversion map of the study area demonstrates that cropland is the primary source of land converted to urban areas, followed by grassland. Among them, the impermeable area of the study area has maintained growth throughout the period. The impervious surface area increased from 54.67km 2 in 1987 to 334.44km 2 in 2019, with average annual growth of 8.5km 2 .

C. SPATIAL AND TEMPORAL DISTRIBUTION OF URBAN IMPERVIOUS SURFACE
Statistics from the increase of the area in the next two years showed that impervious surfaces increased significantly in 1993, 2003, and 2013 ( Figure 11). We analyzed four periods to better understand the annual growth rate of Xining: 1987Xining: to 1992Xining: , 1993Xining: to 2002Xining: , 2003 to 2012, and 2013 to 2019.
• From 1987 to 1993, the area of impervious surfaces was 54.67km 2 , accounting for 3.3% of the total area. The cropland and grassland respectively accounted for 42% and 44.7% of the total area. By 1993, the impervious area was 70.09km 2 , with an average annual growth rate of 4.23%. At this stage, the growth of urban areas was slow, and impervious surfaces in rural areas were scattered with no significant changes.
• From 1994 to 2003, Xining experienced slow urban growth, with an average annual growth rate of 5.34%. Compared with the first period (1987-1993), the annual growth rate was slightly higher, while the growth rate of impervious surfaces remained relatively stable. The area of impervious surfaces increased from 71.48km 2 in 1994 to 114.22km 2 in 2003. Due to the acceleration of industrialization and urbanization, Xining experienced a period of slow urban growth.
• From 2004 to 2013, the annual growth rate was much higher than in the previous two periods. Xining experienced a period of rapid urban expansion, increasing from 117.68 km 2 in 2004 to 228.48 km 2 in 2013, for an average annual growth rate of 7.65%. The large increase was related to the government's 2006 decision to accelerate the development of the Haihu New District and the Nanchuan Industrial Park in the northwest [63]. • From 2014 to 2019, the annual average growth rate was 6.28%, making for slower and more stable urban development. The urban area increased from 246.66km 2 in 2014 to 334.44km 2 in 2019 and is currently 20.4%. In 2013, as the industrial structure was adjusted and upgraded, Xining entered a stage of stable urbanization where the urban fringe was the primary area for urban expansion. Figure 12 displays the expansion of impervious surfaces in eight directions during the four stages (1987-1993, 1993-2003, 2003-2013, and 2013-2019). The expanded area of impervious surfaces is calculated by subtracting the start time from the end time of each period.

D. EXPANSION OF IMPERVIOUS SURFACE IN DIFFERENT DIRECTIONS
From 1987 to 1993, the fastest-growing direction of impervious surfaces was the seventh direction (45 • Northwest), followed by the third direction (45 • Southeast). These two directions form the fastest-growing directions, accounting for 55.33% of the total impervious surface area. Since its implementation in 1983, the Xining City Master Plan (1981-2000 version) has played an important role in guiding urban development and construction. Public service facilities and infrastructure have been established to adapt to Xining's economic structure under China's Western Development strategy and accommodate future social and economic development and construction. From 1993 to 2003, the two fastest-growing directions were still the seventh direction (45 • Northwest), accounting for 25.92% of the total growth of impermeable surfaces and the third direction (45 • Southeast), however, the growth rate of these two directions was slightly lower than in the previous period. Direction five direction (45 • Southwest) expanded faster than the previous period, increasing from 8.89% to 17.51%. Compared with the first stage, the impervious surface expansion across the eight directions fluctuated. For example, the proportion of direction three decreased by 7.36% between the first period and the second period, but the proportion of direction five increased by 8.62% between the first period and the second period. Since 1999, the most significant impact on land coverage in Xining was a policy of returning cropland to forests and grasslands and barren mountain greening projects. The former converted cropland to green land, and the latter increased the vegetation coverage of Xining. The social economy primarily focused on the old city reconstruction plan, meaning Xining City developed from the city center outward to surrounding areas. The primary locations of urban expansion changed considerably from 2003 to 2013. The fastest growth was still observed in the seventh direction (45 • northwest), which accounted for 29.8% of the growth area. The sixth direction (45 • southsouthwest) was the fastest growing direction during this period, accounting for 26.58% of the growth area of the entire impervious surface (an increase of 15.11% over the second stage). Northwest expansion is mainly related to greenway construction and development of the Haihu New Area, which accelerated high-end construction in the service industry. Southwest expansion is mainly due to the expansion of the Chengnan New Town. During 2013-2019, the seventh direction was once again the site of the fastest growth. It is dominated by the Haihu New District, a center for entertainment, culture, education, and sports. In the fourth period, impervious surfaces grew evenly in all eight directions compared to the three periods. The expansion of the eighth direction increased compared with the previous three periods, mainly due to the expansion of North District's high-tech industrial park, the construction of a logistics park, the transformation of the old city, and the construction of a university town, greatly stimulating the expansion of impervious surfaces. In recent years, the construction of the Chengnan New Town and the Haihu New Area had the greatest impact on the Xining's expansion. The development of these two areas will result in multiple centers of commerce in the future.
We then analyzed the direction of expansion of impervious surfaces. The formation and evolution of Xining's landforms trends NW-NWW, which is typical of a branched semi-open city in a valley [25]. The river flows from west to east through the urban area of Xining. The elevation is generally high in the west and low in the east, high in the north and south, low in the middle, and restricted by mountains in the north and south. Urban expansion is limited by the direction and topography of the river valley. A flat river valley lies eastwest, making it suitable for urban construction. Xining's urban space is organized for ''internal life and external production'', meaning that the construction of Xining's living and services are centered around four branch shaped road intersections, which are the most important centers of public service and living in Xining. The expansion of these controls mainly revolves around the central area and spreads to the four directions. At the same time, various industrial facilities were constructed around the urban area, far from the living spaces [65]. From 2001 to 2020, the Xining city master plan was to develop eastward and southward in the nearand medium-term, and to develop westward and northward in the long-term. The East is the most densely populated area in Qinghai Province and the primary route to Lanzhou. Situated to the west of Xining is the Qinghai Tibet railway, which leads to Tibet [66]. From 1987 to 2016, urban expansion trended southeast and northwest, with uneven expansion of impervious surfaces. Additionally, urbanization is accelerating as both the economy and population increase. The Haihu New Area (a center of entertainment, culture, education, and sports) and the southern part of the city (a center of tourism, refined resource processing, ecological agricultural parks, and high-quality residential communities) were rapidly connected and share a contiguous border. Xichuan New City (a center of finance, business, and creative research and development) was rapidly built. The Xining Economic and Technological Development Zone, the Dongchuan Industrial Park and Beichuan Biological Park, which are mainly based on high-tech industries and trade and logistics, have been connected to the central urban area. This model of urban expansion reflects Xining's unique spatial structure, while the government's administrative policy is reflected in the scale of urban expansion. As such, Xining's natural geographical location and administrative orientation have influenced the urban expansion of Xining [65].

E. OVERALL DEVELOPMENT TREND
We produced a map visualizing the distribution of Xining's urban expansion from 1987 to 2019 ( Figure 13). The gradient from red to blue displays the annual changes of impervious surfaces. Red areas are urban lands that developed early, while blue areas are newly developed urban land. The trends of Xining's expansion over the past 33 years can be summarized as follows: (1) Urban development along land that is already developed, (2) Urban development around branch shaped road, (3) Spatial expansion of Xining's central urban area is affected by the terrain. The basic trend is expansion from the center to the northwest, southwest, and southeast. Northeast development was slow. The areas B, C, and D, shown in Figure 13, are schematic diagrams of zoning, which indicate the dynamic expansion of the city in detail. A large amount of land that was cropland in 1987 was urbanized by 2019. Impervious surfaces have expanded from the city center to the suburbs, which is consistent with the spatial and temporal changes of urban expansion obtained from the mapping.
Datasets of annual impervious surfaces over the course of 33 years provide an opportunity to better understand the spatial and temporal patterns of Xining's expansion. Highresolution products provide detailed information on land cover changes in developing regions where urbanization rapidly occurs (1-2 years), compared to products monitoring urbanization across longer time intervals (5-6 years). Rapid urbanization has made it necessary to map impervious surfaces and their progression in order to monitor the expansion of the city.

IV. CONCLUSION
This study is based on Landsat long-term sequence data from 1987-2019 and used random forest machine learning on the GEE platform to identify impervious surfaces in the city of Xining. We used feature optimization to improve classification accuracy, improving the overall accuracy of the area by 2.4%, The overall accuracy was 85.81%. We then used the space-time optimization algorithm to filter errors from the results, improving the accuracy. The final accuracy of impervious surfaces in most years exceeded 90%. Our results determined that the impervious surface area of Xining increased from 55km 2 in 1987 to 334km 2 in 2019. Xining is a typical semi-open city located in a valley, with spatial and temporal features characteristic of urban spaces, primarily developed from the central area to branch shaped road. These results provide effective information for the modeling and future planning of Xining's urban growth. It is recommended to apply feature optimization and time consistency algorithms when drawing impervious surfaces repeatedly.