Time Series Analysis of Land Cover Change Using Remotely Sensed and Multisource Urban Data Based on Machine Learning: A Case Study of Shenzhen, China from 1979 to 2022

: Shenzhen has experienced rapid urbanization since the establishment of the Special Economic Zone in 1978. However, it is rare to witness high-speed urbanization in Shenzhen. It is important to study the LUCC progress in Shenzhen (regarding refusing multisource data), which can provide a reference for governments to solve the problems of land resource shortages and urban expansion spaces. In this paper, nine Landsat images were used to retrieve land cover maps in Shen-zhen, China, from 1979 to 2022. The classiﬁcation method is based on support vector machines with assistance from visual interpretation. The results show that the urban area increased by 756.84 km 2 , the vegetation area decreased by 546.27 km 2 , the water area decreased by 132.95 km 2 , and the bare area decreased by 77.62 km 2 in the last 43 years of our research region. Urban sprawl starts from the Luohu district, then propagates to Futian, Nanshan, and Yantian districts, and ﬁnally expands to other outlying districts (Baoan, Longgang, Guangming, Dapeng, and Pingshan). The spatial–temporal characteristics and the impact factors of urbanization were further analyzed. The visualization of land cover changes based on a complex network approach reveals that the velocity of urban expansion is growing. The coastline distributions were retrieved from nine observation times from 1979 to 2022; the results show that the west coastline changed more dramatically than the east and most of the east coastline remained stable, except for the parts near Yantian port and Mirs Bay, which experienced some changes. The impact factors of coastline changes are further discussed. Through a correlation analysis using urban data, such as transportation and socioeconomic factors, it was found that elevation and roads have strong constraints on the spatial patterns of a city’s expansion. There is exponential decay in the urban land increase against the distance to the roads, implying that trafﬁc factors greatly determine urban land expansion. The turning point of the exponential decay is a distance of around 150 m. Time and population are highly correlated with land use development, indicating that urban land grows linearly with time and the population, which are important driving forces of urban land development. Compared with secondary and tertiary industries, the primary industry is less related to urban land use in Shenzhen.


Introduction
The land is a resource that we rely on for survival in the world.Land use/cover is the result of the interaction between human activities and nature processes [1,2].Land network approach to analyze land use changes from an overall perspective.A complex network approach was widely used in various fields, such as brain networks, social networks, transportation networks, regional innovation networks, disease transmission networks, social relationship networks, and industrial cluster networks [76][77][78].In LUCC fields, the complex network approach can be a very appropriate tool for analyzing the LUCC process, such as the structural characteristics of land use and the relationships between the individual land type behaviors [62].Moreover, a complex network with weights and directions was built to distinguish the key land use types and different conversion relationships of the LUCC progress [79].Meanwhile, the correlations between urban expansion and some natural and socioeconomic factors were investigated.

Study Area and Data
Our study area is Shenzhen city in Guangdong province of southern China (Figure 1), located on the east coast of the Pearl River Delta, bordering Hong Kong to the south, Dongguan to the north, and Huizhou to the northeast.Shenzhen is a global center for research, technology, manufacturing, commerce, finance, tourism, and transportation, and the Port of Shenzhen is the fourth busiest container port in the world [80][81][82].Shenzhen features a subtropical monsoon climate, and the weather is relatively moderate with an annual average temperature of 23.0 • C. The northwestern terrain is lower than the southeastern part.Most areas are covered by low hills.The coastal areas in Baoan and Nanshan districts compose a part of the Pearl River delta plain.The research region covers most of Shenzhen city and some surrounding waters.
cover maps from nine time points from 1979 to 2022.We took advantage of the comp network approach to analyze land use changes from an overall perspective.A comp network approach was widely used in various fields, such as brain networks, soc networks, transportation networks, regional innovation networks, disease transmissi networks, social relationship networks, and industrial cluster networks [76][77][78].In LU fields, the complex network approach can be a very appropriate tool for analyzing LUCC process, such as the structural characteristics of land use and the relationsh between the individual land type behaviors [62].Moreover, a complex network w weights and directions was built to distinguish the key land use types and differ conversion relationships of the LUCC progress [79].Meanwhile, the correlations betwe urban expansion and some natural and socioeconomic factors were investigated.

Study Area and Data
Our study area is Shenzhen city in Guangdong province of southern China (Figu 1), located on the east coast of the Pearl River Delta, bordering Hong Kong to the sou Dongguan to the north, and Huizhou to the northeast.Shenzhen is a global center research, technology, manufacturing, commerce, finance, tourism, and transportati and the Port of Shenzhen is the fourth busiest container port in the world [80-8 Shenzhen features a subtropical monsoon climate, and the weather is relatively moder with an annual average temperature of 23.0 °C.The northwestern terrain is lower than southeastern part.Most areas are covered by low hills.The coastal areas in Baoan a Nanshan districts compose a part of the Pearl River delta plain.The research region cov most of Shenzhen city and some surrounding waters.After rapid development in Shenzhen in recent decades, its economic aggregate has increased from 196 million in 1979 to 3066 billion in 2021, ranking third among Chinese cities.During this period, the permanent population has grown rapidly from 314,000 to 17.68 million in 2021 [82], ranking sixth among Chinese cities.However, the land area of Shenzhen ranked 330th among Chinese cities. Shenzhen has the highest population density in China, with 7173 people per square kilometer [82].The high population density has brought about heavy pressure on the environment.
The data used in our study include the Landsat MSS/TM/ETM/OLI images (Table 1) [83], Shenzhen socioeconomic data [82], and Shenzhen road network data.We used 9 Landsat images, covering 1979 to 2022.In general, the remote sensing images at the start or end of the growing season are optimal, when it is relatively easy to distinguish land type, especially in southern China [84,85].Therefore, most of the Landsat images without clouds were selected around November of that year.The Landsat images were already preprocessed by geometric correction.

Classification System
We grouped the land cover in Shenzhen into four classes (vegetation, urban land, water, and bare land).The spectral characteristics of the TM image are shown in Figure 2.After rapid development in Shenzhen in recent decades, its economic aggregate has increased from 196 million in 1979 to 3066 billion in 2021, ranking third among Chinese cities.During this period, the permanent population has grown rapidly from 314,000 to 17.68 million in 2021 [82], ranking sixth among Chinese cities.However, the land area of Shenzhen ranked 330th among Chinese cities. Shenzhen has the highest population density in China, with 7173 people per square kilometer [82].The high population density has brought about heavy pressure on the environment.
The data used in our study include the Landsat MSS/TM/ETM/OLI images (Table 1) [83], Shenzhen socioeconomic data [82], and Shenzhen road network data.We used 9 Landsat images, covering 1979 to 2022.In general, the remote sensing images at the start or end of the growing season are optimal, when it is relatively easy to distinguish land type, especially in southern China [84,85].Therefore, most of the Landsat images without clouds were selected around November of that year.The Landsat images were already preprocessed by geometric correction.

Classification System
We grouped the land cover in Shenzhen into four classes (vegetation, urban land, water, and bare land).The spectral characteristics of the TM image are shown in Figure 2.  We can see that the four land classes differ greatly from each other in the seven spectra bands.For vegetation, the near-infrared band (band 5) has the largest DN value among all bands.The difference between it and the visible bands (bands 1-4) is very significant.The spectrum of urban land is not dominated by the near-infrared band.Water can absorb electromagnetic waves and, therefore, the values in most bands are lower than other classes, except the values in some bands (bands 1-4) are higher than that of vegetation.Bare land has larger values in most bands than others, except for bands 1-2.In the standard false-color image, vegetation, urban land, water, and bare land appear to be red, celadon, black-blue, and white, respectively.So, they can be distinguished clearly in Landsat images.In some studies, vegetation is further divided into subclasses, such as forest land, wetland, orchard, and grassland [86].We also compare the spectra of wetland and grassland in the Landsat TM image and find that their characteristics are very close (Figure 2).Additionally, the Landsat MSS sensor only has four bands, which makes it more difficult to distinguish them.In order to ensure the reliability of classification results, we decided to not divide vegetation into subclasses.

Land Use Classification
With the development of remote sensing techniques, there are already a number of classification methods, such as decision tree [87], k-nearest neighbor [88], maximum likelihood [89], and artificial neural network [90].Most of these methods are based on pixelbased spectral features and ignore the characteristics of high-scoring remote sensing images (i.e., the rich textures and shapes).A support vector machine (SVM) is a machine learning method based on the statistic studying theory, which has satisfactory performance in a situation with limited training samples.A SVM is a supervised non-parametric statistical learning method, which can solve some classification problems existing in other methods, such as small samples, devilishly-learning, big dimensions, and the local minimum [91].
Previous studies have compared SVMs with maximum likelihood and artificial neural networks in the classification of Landsat ETM+ images [92][93][94].The results revealed that SVMs outperformed other methods in terms of better land-use classification precision.We also tested other methods in our study.After visual interpretations of their results, we found that SVMs can provide results closer to real situations.So, we adopted the SVM method to classify the Landsat images in our study.SVM algorithms use a set of mathematical functions defined as kernels.The function of the kernel is to take data as input and transform them into the desired forms.Different types of kernel functions are used in different SVM algorithms.Since different kernel functions have different characteristics, the classification effects caused by them will also be different.The common kernel functions include linear, polynomial, sigmoid, and Gaussian.Linear kernel: Polynomial kernel: Remote Sens. 2022, 14, x FOR PEER REVIEW 5 of 19 We can see that the four land classes differ greatly from each other in the seven spectra bands.For vegetation, the near-infrared band (band 5) has the largest DN value among all bands.The difference between it and the visible bands (bands 1-4) is very significant.The spectrum of urban land is not dominated by the near-infrared band.Water can absorb electromagnetic waves and, therefore, the values in most bands are lower than other classes, except the values in some bands (bands 1-4) are higher than that of vegetation.Bare land has larger values in most bands than others, except for bands 1-2.In the standard false-color image, vegetation, urban land, water, and bare land appear to be red, celadon, black-blue, and white, respectively.So, they can be distinguished clearly in Landsat images.In some studies, vegetation is further divided into subclasses, such as forest land, wetland, orchard, and grassland [86].We also compare the spectra of wetland and grassland in the Landsat TM image and find that their characteristics are very close (Figure 2).Additionally, the Landsat MSS sensor only has four bands, which makes it more difficult to distinguish them.In order to ensure the reliability of classification results, we decided to not divide vegetation into subclasses.

Land Use Classification
With the development of remote sensing techniques, there are already a number of classification methods, such as decision tree [87], k-nearest neighbor [88], maximum likelihood [89], and artificial neural network [90].Most of these methods are based on pixel-based spectral features and ignore the characteristics of high-scoring remote sensing images (i.e., the rich textures and shapes).A support vector machine (SVM) is a machine learning method based on the statistic studying theory, which has satisfactory performance in a situation with limited training samples.A SVM is a supervised nonparametric statistical learning method, which can solve some classification problems existing in other methods, such as small samples, devilishly-learning, big dimensions, and the local minimum [91].
Previous studies have compared SVMs with maximum likelihood and artificial neural networks in the classification of Landsat ETM+ images [92][93][94].The results revealed that SVMs outperformed other methods in terms of better land-use classification precision.We also tested other methods in our study.After visual interpretations of their results, we found that SVMs can provide results closer to real situations.So, we adopted the SVM method to classify the Landsat images in our study.SVM algorithms use a set of mathematical functions defined as kernels.The function of the kernel is to take data as input and transform them into the desired forms.Different types of kernel functions are used in different SVM algorithms.Since different kernel functions have different characteristics, the classification effects caused by them will also be different.The common kernel functions include linear, polynomial, sigmoid, and Gaussian.
Linear kernel: Polynomial kernel: where  is the degree of the polynomial, ϒ is set as the reciprocal of the band number, and  is the offset of the polynomial kernel.Sigmoid kernel: ,  = ℎ (  + ) where  is the slope of the sigmoid kernel,  is the offset of the sigmoid kernel.Gaussian kernel: We can see that the four land classes differ greatly from each other in the se spectra bands.For vegetation, the near-infrared band (band 5) has the largest DN v among all bands.The difference between it and the visible bands (bands 1-4) is v significant.The spectrum of urban land is not dominated by the near-infrared band.W can absorb electromagnetic waves and, therefore, the values in most bands are lower t other classes, except the values in some bands (bands 1-4) are higher than tha vegetation.Bare land has larger values in most bands than others, except for bands In the standard false-color image, vegetation, urban land, water, and bare land appea be red, celadon, black-blue, and white, respectively.So, they can be distinguished cle in Landsat images.In some studies, vegetation is further divided into subclasses, suc forest land, wetland, orchard, and grassland [86].We also compare the spectra of wetl and grassland in the Landsat TM image and find that their characteristics are very c (Figure 2).Additionally, the Landsat MSS sensor only has four bands, which mak more difficult to distinguish them.In order to ensure the reliability of classification resu we decided to not divide vegetation into subclasses.

Land Use Classification
With the development of remote sensing techniques, there are already a numbe classification methods, such as decision tree [87], k-nearest neighbor [88], maxim likelihood [89], and artificial neural network [90].Most of these methods are based pixel-based spectral features and ignore the characteristics of high-scoring remote sen images (i.e., the rich textures and shapes).A support vector machine (SVM) is a mach learning method based on the statistic studying theory, which has satisfac performance in a situation with limited training samples.A SVM is a supervised n parametric statistical learning method, which can solve some classification probl existing in other methods, such as small samples, devilishly-learning, big dimensions, the local minimum [91].
Previous studies have compared SVMs with maximum likelihood and artif neural networks in the classification of Landsat ETM+ images [92][93][94].The results revea that SVMs outperformed other methods in terms of better land-use classifica precision.We also tested other methods in our study.After visual interpretations of t results, we found that SVMs can provide results closer to real situations.So, we adop the SVM method to classify the Landsat images in our study.SVM algorithms use a se mathematical functions defined as kernels.The function of the kernel is to take dat input and transform them into the desired forms.Different types of kernel functions used in different SVM algorithms.Since different kernel functions have diffe characteristics, the classification effects caused by them will also be different.The comm kernel functions include linear, polynomial, sigmoid, and Gaussian.
Linear kernel: where  is the degree of the polynomial, ϒ is set as the reciprocal of the band num and  is the offset of the polynomial kernel.Sigmoid kernel:   ,  = ℎ (  + ) where  is the slope of the sigmoid kernel,  is the offset of the sigmoid kernel.Gaussian kernel: where d is the degree of the polynomial, Remote Sens. 2022, 14, x FOR PEER REVIEW 5 of 19 We can see that the four land classes differ greatly from each other in the seven spectra bands.For vegetation, the near-infrared band (band 5) has the largest DN value among all bands.The difference between it and the visible bands (bands 1-4) is very significant.The spectrum of urban land is not dominated by the near-infrared band.Water can absorb electromagnetic waves and, therefore, the values in most bands are lower than other classes, except the values in some bands (bands 1-4) are higher than that of vegetation.Bare land has larger values in most bands than others, except for bands 1-2.In the standard false-color image, vegetation, urban land, water, and bare land appear to be red, celadon, black-blue, and white, respectively.So, they can be distinguished clearly in Landsat images.In some studies, vegetation is further divided into subclasses, such as forest land, wetland, orchard, and grassland [86].We also compare the spectra of wetland and grassland in the Landsat TM image and find that their characteristics are very close (Figure 2).Additionally, the Landsat MSS sensor only has four bands, which makes it more difficult to distinguish them.In order to ensure the reliability of classification results, we decided to not divide vegetation into subclasses.

Land Use Classification
With the development of remote sensing techniques, there are already a number of classification methods, such as decision tree [87], k-nearest neighbor [88], maximum likelihood [89], and artificial neural network [90].Most of these methods are based on pixel-based spectral features and ignore the characteristics of high-scoring remote sensing images (i.e., the rich textures and shapes).A support vector machine (SVM) is a machine learning method based on the statistic studying theory, which has satisfactory performance in a situation with limited training samples.A SVM is a supervised nonparametric statistical learning method, which can solve some classification problems existing in other methods, such as small samples, devilishly-learning, big dimensions, and the local minimum [91].
Previous studies have compared SVMs with maximum likelihood and artificial neural networks in the classification of Landsat ETM+ images [92][93][94].The results revealed that SVMs outperformed other methods in terms of better land-use classification precision.We also tested other methods in our study.After visual interpretations of their results, we found that SVMs can provide results closer to real situations.So, we adopted the SVM method to classify the Landsat images in our study.SVM algorithms use a set of mathematical functions defined as kernels.The function of the kernel is to take data as input and transform them into the desired forms.Different types of kernel functions are used in different SVM algorithms.Since different kernel functions have different characteristics, the classification effects caused by them will also be different.The common kernel functions include linear, polynomial, sigmoid, and Gaussian.
Linear kernel: Polynomial kernel: where  is the degree of the polynomial, ϒ is set as the reciprocal of the band number, and  is the offset of the polynomial kernel.Sigmoid kernel:   ,  = ℎ (  + ) where  is the slope of the sigmoid kernel,  is the offset of the sigmoid kernel.Gaussian kernel: is set as the reciprocal of the band number, and α is the offset of the polynomial kernel.Sigmoid kernel: where s is the slope of the sigmoid kernel, α is the offset of the sigmoid kernel.Gaussian kernel: We can see that the four land classes differ greatly from each other in the seven spectra bands.For vegetation, the near-infrared band (band 5) has the largest DN value among all bands.The difference between it and the visible bands (bands 1-4) is very significant.The spectrum of urban land is not dominated by the near-infrared band.Water can absorb electromagnetic waves and, therefore, the values in most bands are lower than other classes, except the values in some bands (bands 1-4) are higher than that of vegetation.Bare land has larger values in most bands than others, except for bands 1-2.In the standard false-color image, vegetation, urban land, water, and bare land appear to be red, celadon, black-blue, and white, respectively.So, they can be distinguished clearly in Landsat images.In some studies, vegetation is further divided into subclasses, such as forest land, wetland, orchard, and grassland [86].We also compare the spectra of wetland and grassland in the Landsat TM image and find that their characteristics are very close (Figure 2).Additionally, the Landsat MSS sensor only has four bands, which makes it more difficult to distinguish them.In order to ensure the reliability of classification results, we decided to not divide vegetation into subclasses.

Land Use Classification
With the development of remote sensing techniques, there are already a number of classification methods, such as decision tree [87], k-nearest neighbor [88], maximum likelihood [89], and artificial neural network [90].Most of these methods are based on pixel-based spectral features and ignore the characteristics of high-scoring remote sensing images (i.e., the rich textures and shapes).A support vector machine (SVM) is a machine learning method based on the statistic studying theory, which has satisfactory performance in a situation with limited training samples.A SVM is a supervised nonparametric statistical learning method, which can solve some classification problems existing in other methods, such as small samples, devilishly-learning, big dimensions, and the local minimum [91].
Previous studies have compared SVMs with maximum likelihood and artificial neural networks in the classification of Landsat ETM+ images [92][93][94].The results revealed that SVMs outperformed other methods in terms of better land-use classification precision.We also tested other methods in our study.After visual interpretations of their results, we found that SVMs can provide results closer to real situations.So, we adopted the SVM method to classify the Landsat images in our study.SVM algorithms use a set of mathematical functions defined as kernels.The function of the kernel is to take data as input and transform them into the desired forms.Different types of kernel functions are used in different SVM algorithms.Since different kernel functions have different characteristics, the classification effects caused by them will also be different.The common kernel functions include linear, polynomial, sigmoid, and Gaussian.
Linear kernel: Polynomial kernel: where  is the degree of the polynomial, ϒ is set as the reciprocal of the band number, and  is the offset of the polynomial kernel.Sigmoid kernel: where  is the slope of the sigmoid kernel,  is the offset of the sigmoid kernel.Gaussian kernel: We can see that the four land classes differ greatly from each other in the sev spectra bands.For vegetation, the near-infrared band (band 5) has the largest DN val among all bands.The difference between it and the visible bands (bands 1-4) is ve significant.The spectrum of urban land is not dominated by the near-infrared band.Wat can absorb electromagnetic waves and, therefore, the values in most bands are lower th other classes, except the values in some bands (bands 1-4) are higher than that vegetation.Bare land has larger values in most bands than others, except for bands 1-In the standard false-color image, vegetation, urban land, water, and bare land appear be red, celadon, black-blue, and white, respectively.So, they can be distinguished clear in Landsat images.In some studies, vegetation is further divided into subclasses, such forest land, wetland, orchard, and grassland [86].We also compare the spectra of wetlan and grassland in the Landsat TM image and find that their characteristics are very clo (Figure 2).Additionally, the Landsat MSS sensor only has four bands, which makes more difficult to distinguish them.In order to ensure the reliability of classification resul we decided to not divide vegetation into subclasses.

Land Use Classification
With the development of remote sensing techniques, there are already a number classification methods, such as decision tree [87], k-nearest neighbor [88], maximu likelihood [89], and artificial neural network [90].Most of these methods are based pixel-based spectral features and ignore the characteristics of high-scoring remote sensi images (i.e., the rich textures and shapes).A support vector machine (SVM) is a machi learning method based on the statistic studying theory, which has satisfacto performance in a situation with limited training samples.A SVM is a supervised no parametric statistical learning method, which can solve some classification problem existing in other methods, such as small samples, devilishly-learning, big dimensions, an the local minimum [91].
Previous studies have compared SVMs with maximum likelihood and artific neural networks in the classification of Landsat ETM+ images [92][93][94].The results reveal that SVMs outperformed other methods in terms of better land-use classificati precision.We also tested other methods in our study.After visual interpretations of the results, we found that SVMs can provide results closer to real situations.So, we adopt the SVM method to classify the Landsat images in our study.SVM algorithms use a set mathematical functions defined as kernels.The function of the kernel is to take data input and transform them into the desired forms.Different types of kernel functions a used in different SVM algorithms.Since different kernel functions have differe characteristics, the classification effects caused by them will also be different.The comm kernel functions include linear, polynomial, sigmoid, and Gaussian.where  is the slope of the sigmoid kernel,  is the offset of the sigmoid kernel.
Previous studies have compared the different kernel functions used in SVMs [95][96][97].The results show that the SVM model with the Gaussian kernel function has the highest performance.We also tested these kernel functions in our study.After a visual interpretation of their results, we found that the Gaussian kernel can give results closer to a real situation.Based on the above analysis, we adopted the SVM model with Gaussian kernel to classify the Landsat images in our study.
In the Landsat image series of Shenzhen, the one captured in 1990 was influenced greatly by the clouds over the southern coastal areas.The clouds were often mistakenly classified as rare land even though the land cover under clouds might belong to vegetation, water, or urban land.The automated classification method cannot give a satisfactory solution to such a problem.In order to ensure classification accuracy, manual intervention was implemented on the classification.We corrected the affected areas using visual interpretation, which is based on the assumptions of spatial and temporal continuity for landscapes.The classes for cloud-covered areas were then inferred from their surrounding landscapes or the classes for them at different times.

Complex Network of LUCC Progress
A complex network with weights and directions was composed of two parts: nodes and connection edges.The nodes of the network in the LUCC progress were the land use types, and the connection edges of the network were the conversion relationships between land use types.The transfer matrix of the LUCC progress can directly reflect the conversion relationships between land use types in different periods.The transfer matrix of the LUCC progress is shown in Equation (5): where i is the land use type before the transfer, j is the land use type after the transfer, and L ij is the transfer value from land use type i to j.In the transfer matrix, every column reflects the transition values from other land use types to land use type j, and every row reflects the transition values from land use type i to other land use types.

Correlation Analysis Using Multisource Urban Data
Furthermore, we conducted a correlation analysis between urban land data and multisource urban data, such as transportation data and socioeconomic data, respectively.The transportation data mainly included the road network map data of Shenzhen city.The socioeconomic data mainly included population, gross domestic product (GDP), the primary industry, secondary industry, and tertiary industry in nine corresponding years.The socioeconomic data mainly came from Shenzhen's statistical yearbook [82].In order to find the relationship between the road network and urban land change, we overlapped the road network map and the urban land expansion map.Meanwhile, through the correlation analysis between urban land and socioeconomic data, we speculated the main impacting factors and driving forces for urbanization in Shenzhen.

Results
As described above, we classified the Landsat images based on SVMs with assistance from visual interpretation.We obtained nine land cover maps from 1979 to 2022 (Figure 3).They showed that Shenzhen experienced a rapid urbanization process and that urban land expanded significantly during this period.Then, we conducted a statistical analysis of the land cover change (Table 2).We found that urban land increased by 743.3% from 101.82 km 2 in 1979 to 858.66 km 2 in 2022, and there was an acceleration in the rate of urban land growth from 1979 to 2003; in particular, the rate of urban land growth averaged 62.95 km 2 per year in 2000-2003.However, the growth fell in 2005 and slightly upgraded in 2008.After that, the growth fell again in 2017 and increased in 2022.This reveals that rapid urbanization at the early stage excessively used land resources.When urban land achieves a certain level, the urbanization speed will gradually decrease.
Remote Sens. 2022, 14, x FOR PEER REVIEW 7 of 19 of the land cover change (Table 2).We found that urban land increased by 743.3% from 101.82 km 2 in 1979 to 858.66 km 2 in 2022, and there was an acceleration in the rate of urban land growth from 1979 to 2003; in particular, the rate of urban land growth averaged 62.95 km 2 per year in 2000-2003.However, the growth fell in 2005 and slightly upgraded in 2008.After that, the growth fell again in 2017 and increased in 2022.This reveals that rapid urbanization at the early stage excessively used land resources.When urban land achieves a certain level, the urbanization speed will gradually decrease.In contrast to the increase in urban land, the vegetation area decreased yearly but at a fluctuant rate.The vegetation area decreased by 48.28% from 1131.36 km 2 in 1979 to  In contrast to the increase in urban land, the vegetation area decreased yearly but at a fluctuant rate.The vegetation area decreased by 48.28% from 1131.36 km 2 in 1979 to 585.09 km 2 in 2022.We found that there was an acceleration in the rate of the vegetation area decline from 1979 to 2003, in particular, the vegetation dropped by 10.19% from 2000 to 2003 with a loss of 36.49km 2 per year.This shows that urbanization in Shenzhen is implemented at a cost of vegetation reduction.
The water area types in our study region included reservoir, river, and sea.Table 2 shows the water area decreased by 24.44% from 543.88 km 2 in 1979 to 410.93 km 2 in 2022.The water areas exhibited a decreasing trend for most years in this period, with the exception of 1990, 2005, and 2022.During this period, the water area dropped by 13.52% from 2000 to 2003 with a loss of 30.26 km 2 per year.The sea areas decreased due to the reclamation activity, but the reservoir areas increased because of the artificial expansion that supported the increasing population in Shenzhen.This suggests that reclamation activity is the main factor that affects the water area in Shenzhen.In contrast to the other three land types, bare land showed no significant trend at this time.Bare land presented a decreasing trend in this period, with the exception of 1990, 2003, 2008, and 2022.We found that bare land increased by 114.86% from 1979 to 1990 with a gain of 13.57km 2 per year.However, bare land decreased by 52.86% from 1990 to 2000 with a loss of 14.76 km 2 per year.
When investigating the spatial characteristics of urbanization based on Figures 3 and 4, we found that most urban expansions from 1979 to 1990 in Shenzhen radiated from the Luohu port (the first border crossing between Shenzhen and Hong Kong), and then gradually propagated outside.Some slight independent expansions can be found around Shenzhen airport, Yantian port, and Shenzhen Bay.This reveals that the early development of Shenzhen was strongly relevant to its interaction with Hong Kong.The facilities construction (e.g., airport and port) related to trade and transit also contributed to the early urbanization characteristics in Shenzhen.Until 2005, urban land expanded to the west of the Luohu district and reached Futian and Nanshan districts.
At the same time, the urban land in coastal areas of Baoan district (Fuyong and Baoan centers), and some inland areas of Longgang district (Longhua and Longgang centers) increased.This suggests that urbanization during the ten years had fully permeated all central districts in Shenzhen, gradually transferring to outlying districts.The Tiegang reservoir and Xili reservoir experienced significant expansion, implying that the requirements for freshwater continued to increase during the Shenzhen development.Until 2013, the increase in urban land mainly occurred in the Baoan and Longgang districts.The land cover changes in the central districts were represented by the transformation from the sea area to urban land on Qianhai Bay and Houhai Bay, implying that urban land use in central districts had reached a critical point.Until 2022, the increase in urban land mainly occurred in Baoan, Longgang, and Pingshan districts; the land use changes of this period are still dominated by urbanization in outlying districts and reclamation in central districts, which can be regarded as a continuation of the previous process.reclamation in central districts, which can be regarded as a continuation of the previous process.

Complex Network Analysis on Land Cover Change
The statistics of land cover change on a single class are not enough to represent the shift between different classes for the whole land cover system.Most studies use the transfer matrix to describe the characteristics of land cover changes among classes, but it is still not intuitive and easily understood.The complex network approach is an important method for the study of sociology and informatics [98], which has also recently been introduced to analyze land cover changes [99].It can visually reveal the relationship between different land cover classes from an overall perspective.We took advantage of this method to describe the shift between different land cover classes in Shenzhen.Each node in the network represents a land use type.The line denotes the shift between two land types, and the arrow indicates the transforming direction.The transformation is represented by the thickness of the line and arrow.Figure 5 presents the complex network visualization of the land cover transformation between the nine observation time points.It can be clearly seen that the transformation from vegetation to bare land dominates the land cover changes during 1979-1990.Because bare land acts as the medium for turning vegetation into urban land, the urbanization of Shenzhen in this period can be characterized as the preparation for upcoming urban construction.
From 1990 to 2000, the transformation from bare land and vegetation to urban land dominated the land cover change.This suggests that bare land prepared before was transformed into urban land.During the ten years, many areas covered by vegetation quickly experienced two transformations, namely vegetation to bare land and bare land to urban land.This implies that construction velocity was much faster than the decade before.In the two periods (2000-2003 and 2003-2005), the land cover changes were only dominated by that transformation from vegetation to urban.The bare land medium during urban construction could hardly be witnessed, even though the time intervals

Discussions 4.1. Complex Network Analysis on Land Cover Change
The statistics of land cover change on a single class are not enough to represent the shift between different classes for the whole land cover system.Most studies use the transfer matrix to describe the characteristics of land cover changes among classes, but it is still not intuitive and easily understood.The complex network approach is an important method for the study of sociology and informatics [98], which has also recently been introduced to analyze land cover changes [99].It can visually reveal the relationship between different land cover classes from an overall perspective.We took advantage of this method to describe the shift between different land cover classes in Shenzhen.Each node in the network represents a land use type.The line denotes the shift between two land types, and the arrow indicates the transforming direction.The transformation is represented by the thickness of the line and arrow.Figure 5 presents the complex network visualization of the land cover transformation between the nine observation time points.It can be clearly seen that the transformation from vegetation to bare land dominates the land cover changes during 1979-1990.Because bare land acts as the medium for turning vegetation into urban land, the urbanization of Shenzhen in this period can be characterized as the preparation for upcoming urban construction.
From 1990 to 2000, the transformation from bare land and vegetation to urban land dominated the land cover change.This suggests that bare land prepared before was transformed into urban land.During the ten years, many areas covered by vegetation quickly experienced two transformations, namely vegetation to bare land and bare land to urban land.This implies that construction velocity was much faster than the decade before.In the two periods (2000-2003 and 2003-2005), the land cover changes were only dominated by that transformation from vegetation to urban.The bare land medium during urban construction could hardly be witnessed, even though the time intervals were shorter than the two periods before.This reveals that the transformation from vegetation to urban became faster.
over time with the reduction of available land resources.From 2013 to 2017, the land cover change was only dominated by the transformation from vegetation to urban.During this period, the transformation from bare land to urban land was slower, implying that the bare land prepared before was transformed into urban land.From 2017 to 2022, it can be clearly seen that the transformation from vegetation to bare land dominated the land cover change, and the transformation from vegetation to urban was slightly faster than the 2013-2017 period, implying that the velocity of urban expansion is growing over time.The rapid urbanization of Shenzhen is also accompanied by the dramatic change in its coastline.After the visual interpretation of the Landsat images, we retrieved the coastline distributions from the nine observation times from 1979 to 2022 (Figure 6).During that period, Shenzhen started its rapid urbanization and industrialization process.With the rapid expansion of the urban construction land scale, the contradiction between man and land became increasingly prominent.Reclamation and land-building projects became the main means to solve the problem of land resource supply exhaustion in coastal zones, resulting in drastic changes in the shoreline in some areas.As shown in Figure 6, the west coastline obviously changed more dramatically than the east coastline, i.e., it kept extending outside during the observing period.However, most of the east coastline remained stable, except for the parts near Yantian port and Mirs Bay, which experienced some changes.Such asymmetry along the coastline resulted from its natural characteristics.
The west coastline is located to the east of the Pearl River Estuary.Most are silty coasts, providing favorable conditions for reclamation engineering.However, the east At the same time, the transformation from vegetation to urban was significantly slower than the decade before, implying that the velocity of urban expansion is decreasing over time with the reduction of available land resources.From 2013 to 2017, the land cover change was only dominated by the transformation from vegetation to urban.During this period, the transformation from bare land to urban land was slower, implying that the bare land prepared before was transformed into urban land.From 2017 to 2022, it can be clearly seen that the transformation from vegetation to bare land dominated the land cover change, and the transformation from vegetation to urban was slightly faster than the 2013-2017 period, implying that the velocity of urban expansion is growing over time.
The rapid urbanization of Shenzhen is also accompanied by the dramatic change in its coastline.After the visual interpretation of the Landsat images, we retrieved the coastline distributions from the nine observation times from 1979 to 2022 (Figure 6).During that period, Shenzhen started its rapid urbanization and industrialization process.With the rapid expansion of the urban construction land scale, the contradiction between man and land became increasingly prominent.Reclamation and land-building projects became the main means to solve the problem of land resource supply exhaustion in coastal zones, resulting in drastic changes in the shoreline in some areas.As shown in Figure 6, the west coastline obviously changed more dramatically than the east coastline, i.e., it kept extending outside during the observing period.However, most of the east coastline remained stable, except for the parts near Yantian port and Mirs Bay, which experienced some changes.Such asymmetry along the coastline resulted from its natural characteristics.
coastline is mainly composed of rocky coast, which is difficult for construction.The purpose of reclamation on the east coastline mainly involves the Yantian port construction.Conversely, the functions of the west coastline are more complex.On the Shajin-Xinan coastline, the reclamations mainly involve the mariculture zone expansion.On the central Qianhai Bay, the reclaimed lands are planned for various uses, including industry development, public facilities building, and port construction.The Qianhai Shenzhen-Hong Kong Service Industry Cooperation Zone is located in the reclaimed land of Qianhai Bay.Shenzhen Bay lies on the east of Qianhai Bay.The reclaimed lands here are mainly used for building public facilities.Mangrove Park, Shenzhen Bay Park, and Overseas Chinese Town are all located there.
Based on the above analysis, we created a 3 km buffer zone of the coastlines as shown in Figure 6.The complex network approach was applied to analyze the land cover transformation in the buffer zone (Figure 7).It is shown that an obvious transformation from water to urban land occurred during the periods (2000-2003 and 2005-2008), but the water along the coastline transforming to other land types significantly decreased from 2008 to 2022.It can be clearly seen that the transformation from urban land to water occurred during 2003-2005.This reveals that the reservoir areas increased because of the artificial expansion to support the increasing population during this period.Corresponding to that trend is the water area, which increased from 451.11 km 2 in 2003 to 464.13 km 2 in 2005, as shown in Table 2.
The vegetation along the coastline transformed into other land types significantly in the early stages of Shenzhen's development.Such a transformation gradually decreased over time.Meanwhile, the transformation from bare land to urban land became slower, implying that the bare land prepared before was transformed into urban land.The bare land along the coastline transforming to other land types significantly decreased from 2008 to 2022.Increased urban land through reclamation activity played an important role in promoting regional economic development, but this change had negative impacts on the ecological environment.In the 2003-2022 period, other land types changed to vegetation.This reveals that Shenzhen government is paying increasing attention to the protection of the coastal zone ecosystem.The west coastline is located to the east of the Pearl River Estuary.Most are silty coasts, providing favorable conditions for reclamation engineering.However, the east coastline is mainly composed of rocky coast, which is difficult for construction.The purpose of reclamation on the east coastline mainly involves the Yantian port construction.Conversely, the functions of the west coastline are more complex.On the Shajin-Xinan coastline, the reclamations mainly involve the mariculture zone expansion.On the central Qianhai Bay, the reclaimed lands are planned for various uses, including industry development, public facilities building, and port construction.
The Qianhai Shenzhen-Hong Kong Service Industry Cooperation Zone is located in the reclaimed land of Qianhai Bay.Shenzhen Bay lies on the east of Qianhai Bay.The reclaimed lands here are mainly used for building public facilities.Mangrove Park, Shenzhen Bay Park, and Overseas Chinese Town are all located there.
Based on the above analysis, we created a 3 km buffer zone of the coastlines as shown in Figure 6.The complex network approach was applied to analyze the land cover transformation in the buffer zone (Figure 7).It is shown that an obvious transformation from water to urban land occurred during the periods (2000-2003 and 2005-2008), but the water along the coastline transforming to other land types significantly decreased from 2008 to 2022.It can be clearly seen that the transformation from urban land to water occurred during 2003-2005.This reveals that the reservoir areas increased because of the artificial expansion to support the increasing population during this period.Corresponding to that trend is the water area, which increased from 451.11 km 2 in 2003 to 464.13 km 2 in 2005, as shown in Table 2.
The vegetation along the coastline transformed into other land types significantly in the early stages of Shenzhen's development.Such a transformation gradually decreased over time.Meanwhile, the transformation from bare land to urban land became slower, implying that the bare land prepared before was transformed into urban land.The bare land along the coastline transforming to other land types significantly decreased from 2008 to 2022.Increased urban land through reclamation activity played an important role in promoting regional economic development, but this change had negative impacts on the ecological environment.In the 2003-2022 period, other land types changed to vegetation.This reveals that Shenzhen government is paying increasing attention to the protection of the coastal zone ecosystem.

Subsection Traffic Factors in Urban Land Change
The traffic network exhibits obvious constraints on urban land use.When overlapping the road network map and the urban land expansion map (1979-2022), an obvious coupling can be found (Figure 8a).The urban land use increase was mainly located around the traffic network.We applied a histogram analysis on the transportation distance (Figure 8b), and found that

Subsection Traffic Factors in Urban Land Change
The traffic network exhibits obvious constraints on urban land use.When overlapping the road network map and the urban land expansion map (1979-2022), an obvious coupling can be found (Figure 8a).
The urban land use increase was mainly located around the traffic network.We applied a histogram analysis on the transportation distance (Figure 8b), and found that there was an exponential decay of urban land increase against the distance to the roads, implying that traffic factors greatly determine the urban land expansion.The turning point of the exponential decay was a distance of around 150 m.This suggests that the traffic effects on urban land expansion reduce sharply with a transportation distance beyond 150 m.

Socioeconomic Factors in Urban Land Change
Urban land use is also affected by socioeconomic factors, such as time, population, gross domestic product (GDP), primary industry, secondary industry, and tertiary industry.We implemented a correlation analysis between these factors and our retrieved urban land areas, respectively.The analysis was based on the socioeconomic data from Shenzhen's statistical yearbook [82] in the nine corresponding years.
It can be seen from Figure 9 that time and population were highly related to land use development, indicating that urban land grew linearly with time and that the population factor is an important driving force in urban land development.Meanwhile, we can see that the weakest relevance existed in the primary industry.A negative correlation even appeared between 2003 and 2013.This suggests that the land use in Shenzhen became less related to the primary industry.Compared with the primary industry, the secondary and tertiary industries have stronger correlations with urban land use in Shenzhen.

Subsection Traffic Factors in Urban Land Change
The traffic network exhibits obvious constraints on urban land use.When overlapping the road network map and the urban land expansion map (1979-2022), an obvious coupling can be found (Figure 8a).The urban land use increase was mainly located around the traffic network.We applied a histogram analysis on the transportation distance (Figure 8b), and found that there was an exponential decay of urban land increase against the distance to the roads, implying that traffic factors greatly determine the urban land expansion.The turning point of the exponential decay was a distance of around 150 m.This suggests that the traffic effects on urban land expansion reduce sharply with a transportation distance beyond 150 m.

Socioeconomic Factors in Urban Land Change
Urban land use is also affected by socioeconomic factors, such as time, population, gross domestic product (GDP), primary industry, secondary industry, and tertiary industry.We implemented a correlation analysis between these factors and our retrieved urban land areas, respectively.The analysis was based on the socioeconomic data from Shenzhen's statistical yearbook [82] in the nine corresponding years.
It can be seen from Figure 9 that time and population were highly related to land use development, indicating that urban land grew linearly with time and that the population factor is an important driving force in urban land development.Meanwhile, we can see that the weakest relevance existed in the primary industry.A negative correlation even appeared between 2003 and 2013.This suggests that the land use in Shenzhen became less related to the primary industry.Compared with the primary industry, the secondary and tertiary industries have stronger correlations with urban land use in Shenzhen.

Conclusions
The LUCC statistics on a single class are not enough to represent the shift between different classes for the whole land cover system.Most studies explored the LUCC progress from certain land type view changes, ignoring the integrality of the land use system and the correlation between land use and other factors.Moreover, most studies used the transfer matrix to describe the characteristics of land cover changes among classes, but it is still not intuitive or easily understood.In this paper, we studied the land cover and its changes in Shenzhen from 1979 to 2022, using Landsat time series images.The spatial-temporal characteristics and the impact factors of urbanization were further analyzed.The main conclusions are as follows: First, Shenzhen experienced a rapid urbanization process from 1979 and 2022.In the study area, the cumulative increase in urban land was about 756.84 km 2 , the reduction of vegetation was about 546.27 km 2 , the reduction of water was about 132.95 km 2 , and the reduction of bare land was about 77.62 km 2 .The urbanization started in the Luohu district, gradually radiated to the urban districts (Futian, Nanshan, and Yantian districts), and finally reached the outlying districts (Baoan, Longgang, Guangming, Dapeng, and Pingshan).Rapid urbanization at an early stage results in excessive use of land resources; however, when urban land achieves a certain level, the urbanization speed will gradually decrease with time.
Second, the complex network visualization reveals that the main changes in land cover are different over time.The transformation from vegetation to bare land is hard to witness in the late stage.As bare land is the medium for transforming vegetation into urban land, it means that the velocity of vegetation-urban land transformation keeps increasing.We can conclude that the urbanization of Shenzhen continuously speeds up with time.Meanwhile, the coastline distributions show that the west coastline changed more dramatically than the east coastline; most of the east coastline remained stable, except for the parts near Yantian port and Mirs Bay, which experienced some changes.This reveals that an obvious transformation from water to urban land occurred during the early period, but the water along the coastline transforming to other land types significantly decreased in recent years.
Finally, the road network greatly impacted the land cover changes from 1979 to 2022.A histogram analysis was applied to the transportation distance, and we found an exponential decay of urban land increase against the distance to the roads.The turning point of the exponential decay was a distance of around 150 m.This suggests that the traffic effect on urban land expansion reduces sharply with a transportation distance beyond 150 m.By investigating the socioeconomic data, we found that urban land grows linearly with time and the population is highly correlated with urban area changes.The primary industry, compared with the secondary and tertiary industries, is less related to urban land use in Shenzhen.These socioeconomic factors can be considered the main driving forces for urbanization in Shenzhen.In the future, under the current economic policies, the population and economy of Shenzhen will continue to increase.However, the increase will be accompanied by a decrease in vegetation and water resources.According to this present study, the LUCC changes in Shenzhen will gradually enter a stable period.Thus, sustainable land use policies should be formulated dynamically according to the socioeconomic realities of Shenzhen.Although each city has a different development stage and faces different problems, the research method of this paper can provide a reference for other cities.

Figure 1 .
Figure 1.The TM false color image in the research region.

Figure 1 .
Figure 1.The TM false color image in the research region.

Figure 2 .
Figure 2. Spectral characteristics of the Landsat TM image for different land types.

Figure 2 .
Figure 2. Spectral characteristics of the Landsat TM image for different land types.
Linear kernel:   ,  =   ( Polynomial kernel:   ,  = (ϒ  + ) , ϒ > 0 ( where  is the degree of the polynomial, ϒ is set as the reciprocal of the band numb and  is the offset of the polynomial kernel.Sigmoid kernel:   ,  = ℎ (  + ) (

Figure 3 .
Figure 3.The land cover maps in Shenzhen from 1979 to 2022.

Figure 3 .
Figure 3.The land cover maps in Shenzhen from 1979 to 2022.

Figure 4 .
Figure 4.The administrative map of Shenzhen city.

Figure 4 .
Figure 4.The administrative map of Shenzhen city.

Figure 5 .
Figure 5.The complex network visualization of the land cover transformation in Shenzhen from 1979 to 2022.

Figure 5 .
Figure 5.The complex network visualization of the land cover transformation in Shenzhen from 1979 to 2022.In the two periods(2005-2008 and 2008-2013), the transformation from bare land and vegetation to urban land dominated the land cover changes.This implies that the construction velocity was slightly slower during the 2008-2013 period than the 2008-2013 period.At the same time, the transformation from vegetation to urban was significantly slower than the decade before, implying that the velocity of urban expansion is decreasing over time with the reduction of available land resources.From 2013 to 2017, the land cover change was only dominated by the transformation from vegetation to urban.During this period, the transformation from bare land to urban land was slower, implying that the bare land prepared before was transformed into urban land.From 2017 to 2022, it can be clearly seen that the transformation from vegetation to bare land dominated the land cover change, and the transformation from vegetation to urban was slightly faster than the 2013-2017 period, implying that the velocity of urban expansion is growing over time.The rapid urbanization of Shenzhen is also accompanied by the dramatic change in its coastline.After the visual interpretation of the Landsat images, we retrieved the coastline distributions from the nine observation times from 1979 to 2022 (Figure6).During that period, Shenzhen started its rapid urbanization and industrialization process.With the rapid expansion of the urban construction land scale, the contradiction between man and land became increasingly prominent.Reclamation and land-building projects became the main means to solve the problem of land resource supply exhaustion in coastal zones, resulting in drastic changes in the shoreline in some areas.As shown in Figure6, the west coastline obviously changed more dramatically than the east coastline, i.e., it kept extending outside during the observing period.However, most of the east coastline remained stable, except for the parts near Yantian port and Mirs Bay, which experienced some changes.Such asymmetry along the coastline resulted from its natural characteristics.

Figure 6 .
Figure 6.The coastline changes in Shenzhen from 1979 to 2022.The shadow area denotes the 3 km buffer zone of the coastline.

Figure 6 .
Figure 6.The coastline changes in Shenzhen from 1979 to 2022.The shadow area denotes the 3 km buffer zone of the coastline.

Figure 7 .
Figure 7.The complex network visualization of the land cover transformation in the 3 km buffer zone of the Shenzhen coastline from 1979 to 2022.

Figure 8 .
Figure 8.The relationship between transportation and urban land increase from 1979 to 2022.(a) The road network overlaid on the map of urban land expansion.(b) The transportation distance histograms for urban land increase.

Figure 7 .
Figure 7.The complex network visualization of the land cover transformation in the 3 km buffer zone of the Shenzhen coastline from 1979 to 2022.

Figure 7 .
Figure 7.The complex network visualization of the land cover transformation in the 3 km buffer zone of the Shenzhen coastline from 1979 to 2022.

Figure 8 .
Figure 8.The relationship between transportation and urban land increase from 1979 to 2022.(a) The road network overlaid on the map of urban land expansion.(b) The transportation distance histograms for urban land increase.

Figure 8 .
Figure 8.The relationship between transportation and urban land increase from 1979 to 2022.(a) The road network overlaid on the map of urban land expansion.(b) The transportation distance histograms for urban land increase.

Figure 9 .
Figure 9.The correlation analysis between socioeconomic factors and the urban land area.

Table 1 .
Parameters of the used Landsat images.

Table 1 .
Parameters of the used Landsat images.

Table 2 .
Statistics on the land cover of Shenzhen at different times.