*Article* **Thermal Environment Effects of Built-Up Land Expansion in Shijiazhuang**

**Ling Qin 1,2,3, Han Liu 3, Guofei Shang 1,2,3, Huicai Yang 2,3,\* and Haiming Yan 2,3**


**Abstract:** Exploring the thermal environment effects of built-up land expansion can lay a firm foundation for urban planning and design. This study revealed the spatiotemporal dynamic characteristics of built-up land and heat island center points in Shijiazhuang using land-use/land-cover data and land surface temperature (LST) products from 1996 to 2019, and the response mechanism between the percentage of built-up land (PLAND) and LST with the grid sampling method and statistical analysis. Results indicated that heat islands are mainly clustered in the downtown, built-up areas of counties and the Hutuo River Basin. The spatiotemporal shift direction of the center point of the urban heat island (UHI) and built-up land in the whole study area varied due to the eco-environmental transformation of the Hutuo River Basin. In areas far from the Hutuo River Basin, the center points of UHI and built-up land were shifted in a similar direction. There is a remarkable linear correlation between the PLAND and LST, the correlation coefficient of which was higher than 0.7 during the study period. Areas with PLAND > 60% are urban regions with stronger heat island effects, and areas with PLAND < 55% are villages and towns where the temperature raised more slowly.

**Keywords:** surface thermal environment; standard deviation ellipse; Shijiazhuang

#### **1. Introduction**

The expansion of built-up land along with the acceleration of industrialization and urbanization has led to a series of ecological problems, among which the thermal environment effects were most closely related to human life and directly affected individuals' normal life [1]. The concept of the urban heat island (UHI) has attracted the attention of many scholars around the world [2]. The acquisition of the land surface temperature (LST) based on remote sensing images provides firm data support for exploring the UHI [3,4], and spatial analysis technology has provided powerful tools for effectively revealing thermal environment effects of land-use/land-cover change along with urbanization [5–7].

A number of previous studies have shown that land-use/land-over change has greatly affected the spatial pattern of the thermal environment [8,9]. Most scholars have generally focused on characterizing the evolution of thermal environment effects of land-use/landcover change with the distribution condition of landscape pattern indices [3,8], but some previous studies also indicated that not all landscape pattern indices are applicable and the surface thermal environment is predominantly associated with the percentage of built-up land (PLAND). In fact, the percentage of vegetation and the PLAND were both highly correlated with surface temperature [9–13]. Besides this, scholars have conducted studies related to thermal environment effects of the land-use/land-cover change with various

**Citation:** Qin, L.; Liu, H.; Shang, G.; Yang, H.; Yan, H. Thermal Environment Effects of Built-Up Land Expansion in Shijiazhuang. *Land* **2022**, *11*, 968. https:// doi.org/10.3390/land11070968

Academic Editor: Luca Salvati

Received: 31 May 2022 Accepted: 20 June 2022 Published: 24 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

methods in recent decades [14,15]. The correlation analysis of the surface temperature and remote sensing index have been widely used, and Autocorrelation using Getis-Ord Gi\* was also widely used to explore the relationship between the built-up land and LST [14–16]. However, the above-mentioned studies generally provide limited information on the response mechanism between the PLAND and LST, and it is of particular significance to deeply understand the impacts of built-up land expansion on the surface thermal environment in fast-growing cities [17–20]. For example, Chen et al. suggested that the PLAND = 35% split the heating effects of the PLAND on LST by grid sampling and statistics approach in Wu'an [18].

Shijiazhuang, as the capital of Hebei Province of China, has a very high urbanization rate, with 70.18% of the resident population in the urban area in 2020, and the heat island effect due to urban expansion is also increasingly prominent, so it is extremely meaningful to explore the response mechanism between the PLAND and LST in this highly urbanized region. This study adopted the grid sampling statistics approach and took the downtown, two districts, and two counties with fast economic development in Shijiazhuang as the study area, combining the standard deviation ellipse (SDE), spatial autocorrelation, and correlation analysis method to explore the effects of the PLAND on LST, aiming to provide a firm theoretical foundation for the urban planning and design of Shijiazhuang.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Shijiazhuang is located at the eastern piedmont of the Taihang Mountains, and the eastward airflow is blocked by the mountains and then sinks, causing the temperature in the area at the piedmont of the mountains to rise. Shijiazhuang, as the capital of Hebei Province, has experienced rapid urbanization since the 1990s, with large-scale urban land expansion and migration of people to the center of the city for better public services and job opportunities. The study area in this study covers the downtown (Xinhua, Changan, Qiaoxi, Yuhua), two districts (Gaocheng, Luancheng), and two counties (Wuji, Zhengding) (Figure 1) in the east and central parts of Shijiazhuang, with plain terrain outside the mountain ranges, which jointly account for 47.66% of the total population of Shijiazhuang.

**Figure 1.** Location of the study area.

#### *2.2. Data Sources and Processing*

The major data sources used in this study include: (1) Landsat series images from the Google Earth Engine (GEE) platform (https://earthengine.google.com/, accessed on 31 March 2022). Images covering the extent of the study area were selected in six periods (1996, 2001, 2007, 2011, 2015, and 2019). The data parameters are shown in Table 1. (2) Landuse/land-cover data from the China Land Use/Cover Dataset (CLUDs) by Professor Huang Xin's team at Wuhan University. The spatial resolution is 30 m [21], and the data


were obtained from the PIE-Engine platform (https://engine.piesat.cn/, accessed on 31

**Table 1.** Remote sense statistic parameters.

March 2022).

The land-use and land-cover data were programmed to be downloaded locally with the application programming interface of the PIE Engine platform. Then this study reclassified the land-use and land-cover data into four categories, that is, cropland, built-up land, water, and other land (Figure 2), with the help of the reclassification tool in the ArcGIS software. The statistics of the built-up land area for each year were also summarized based on the ArcGIS software (Figure 3).

**Figure 2.** Land-use distribution in the study area from 1996 to 2019.

**Figure 3.** Time series statistic diagram of built-up land area.

This study selected the Landsat images in six less cloudy periods as the foundation for the temperature inversion and acquired 12 scenes of LST single-band images in six periods in total by running the open-source code. This study then mosaicked these remote sensing images of each period together with the help of a seamless mosaicking tool in the ENVI software, and thereafter the data of the study area were extracted based on the study area boundary data (Figure 4). This study thereafter used the practical single channel (PSC) method to acquire the LST data with the help of the GEE cloud platform and the Landsat series satellite surface temperature data inversion framework [21,22]. The overall bias of the obtained temperature products was 0.17 K according to the cross-validation with MODIS surface temperature products, indicating high accuracy.

**Figure 4.** Spatial pattern of the land surface temperature (LST) in the study area from 1996 to 2019.

#### *2.3. Calculation of the Percentage of Built-Up Land*

This study explored the urban thermal environment based on the landscape pattern index, which is a quantitative expression to describe the spatial pattern of the landscape [8,9]. The landscape pattern index is one of the current research hotspots to investigate the urban thermal environment, and previous studies generally suggested that the PLAND index has the best fitting effect on the thermal environment by investigating the correlation between each landscape pattern index and LST [20]. This study specifically described the spatial and temporal distribution pattern of built-up land in the study area with the PLAND, which was calculated based on the moving window approach of the FRAGSTATS software as follows.

$$PLAND = \left(\sum\_{j=1}^{n} a\_{ij} / A\right) \times 100\tag{1}$$

where the PLAND is the percentage of built-up land, which ranges between 0 and 100, *aij* represents the area of built-up land patches, *A* indicates the background area. When the value of the index is close to 0, it indicates that the PLAND is decreasing in a single calculation window, and when the value of a single window is 100, it means that the window consists of built-up land only.

#### *2.4. Thermal Environment Grading Classification*

Climatic differences over the years have significant impacts on the true LST, such differences need to be eliminated when analyzing the evolution of the thermal environment, and the equation for standardization is as follows.

$$T\_{nor} = \frac{T\_i - T\_{min}}{T\_{max} - T\_{min}} \tag{2}$$

where *Tnor* denotes the normalized temperature of the *i th* grid, *Ti* indicates the temperature of the *i th* grid, *Tmax* indicates the maximum temperature in the study area, *Tmin* indicates the minimum temperature in the study area. The mean-standard deviation method was used to classify the thermal environment in the study area, and the LST was divided into five thermal classes: low, sub-low, medium, sub-high temperature, and high temperature (Table 2). By grading the LST, this study generated a map of the thermal class distribution of each period, which makes the thermal environment of each period comparable and lays a firm foundation for the in-depth investigation of the spatial and temporal evolution of cold and heat islands in the study area.

**Table 2.** Land surface temperature (LST) grading criteria.


Note: *Tnor* denotes the normalized temperature of the *i th* grid, *Ta* indicates the average temperature in the study area for each year, and *sd* represents the standard deviation of the LST.

#### *2.5. Establishment of the Standard Deviation Ellipse*

This study explored the evolution trends of the built-up land and heat island using the SDE, a classical method widely used for analyzing the directional characteristics of spatial distribution. In this study, the high-temperature region and the sub-high-temperature region were divided into heat island areas with the spatial statistics module of the ArcGIS software, the direction publishing tool was used to obtain the standard deviation of the heat island and built-up land. The expansion characteristics of the heat island and built-up land were explored by connecting the center points and SDE of the heat island and built-up land. The comparative analysis of the center point shift direction was carried out to lay the foundation for the following quantitative study of the relationship between built-up land expansion and the surrounding thermal environment.

#### *2.6. Spatial Autocorrelation Analysis and Correlation Analysis*

A grid scale of 990 × 990 m was selected based on the Fishnet tool of the ArcGIS software to establish spatial grids to investigate the response mechanism between the PLAND and LST at the same spatial scale [18,19]. The sample grid was collected by combining the PLAND and LST raster maps for each year. It is notable that the thermal field in the Hutuo River Basin may disturb the surrounding thermal environment. Therefore, the interfering areas were excluded from the grid sampling, and the sample grid covered all the areas where the PLAND existed in the study area.

The spatial autocorrelation analysis of the distribution of the PLAND and LST is a prerequisite for exploring the correlation between the PLAND and LST. This study explored the aggregation characteristics of the spatial distribution of built-up land and LST based on Moran's I, which can effectively detect the spatial autocorrelation of geographic entities. This study estimated the Moran's I of built-up land and LST with the spatial autocorrelation tool in the ArcGIS software as follows.

$$\text{Moran\'s I} = \frac{\text{n}}{\sum\_{i=1}^{n} \sum\_{j=1}^{n} w\_{ij}} \times \frac{\sum\_{i=1}^{n} \sum\_{j=1}^{n} w\_{ij} \left(X\_i - \overline{X}\right) \left(X\_j - \overline{X}\right)}{\sum\_{i=1}^{n} \left(X\_i - \overline{X}\right)^2} \tag{3}$$

where *n* is the number of samples, *wij* is the spatial weight matrix, *Xi* indicates the temperature value and the PLAND value of each period.

This study further explored the correlation between the PLAND and the corresponding LST using the Pearson correlation analysis as follows.

$$\tau = \frac{\sum\_{i=1}^{n} (\mathbf{x}\_i - \overline{\mathbf{x}}) \left(\mathbf{y}\_j - \overline{\mathbf{y}}\right)}{\sqrt{\sum\_{i=1}^{n} (\mathbf{x}\_i - \overline{\mathbf{x}})^2} \sqrt{\sum\_{i=1}^{n} (\mathbf{y}\_i - \overline{\mathbf{y}})^2}} \tag{4}$$

where *r* represents the correlation between the two sets of variables, *x* represents the arithmetic mean of the LST of the *i th* grid, *yi* is the percentage of built-up land in the *i th* grid, *y* is the arithmetic mean of the PLAND in the *i th* grid.

#### **3. Results**

#### *3.1. Characteristics of the Distribution of Thermal Environment*

The results indicated the spatial pattern of the thermal environment in six districts and two counties of Shijiazhuang has changed dramatically along with accelerated urbanization and industrialization (Figure 5). In this study, the sub-high temperature and high-temperature zones were categorized as the heat island, the medium-temperature zones were categorized as the normal zone, and the low-temperature zone and sub-low temperature zone were categorized as the cold island. The spatial and temporal expansion of the thermal environment in downtown Shijiazhuang was the most remarkable, indicating the strong heat island effect. Besides this, the heat island range in the main urban area of Shijiazhuang tended to expand to the east as time went by, followed by a large expansion of the heat island in the main urban area to the southeast, which merged with the heat island patches in Luancheng District and reflected the development orientation towards Luancheng District. In addition, the Hutuo River Basin has been maintained in a high-temperature zone state from 1996 to 2011 due to its poor ecological environment and bad geological condition within the dry riverbed.

Cold islands are primarily located in rural and rural farmland far from the builtup areas of districts and towns, which are more widely distributed in the northeast and southwest of the downtown. Meanwhile, the area of cold islands decreased gradually with the expansion of cities and tow. In particular, cold islands appeared in the downtown section of the Hutu River Basin after 2011, which expanded to the east over time due to the ecological restoration project of the Hutuo River Basin, reaching its maximum extent in 2019.

**Figure 5.** Distribution of heat levels from 1996 to 2019.

It is notable that the proportion of sub-high temperature in the heat island area of downtown in 2011 is larger than that of high temperature (Figures 4 and 5), which is primarily attributed to meteorological issues such as the cloud content. In fact, it is not surprising that the heat island area in 2011 is larger than that in 2007 since areas with sub-high temperature and high temperature were together classified as the heat island in this study, and this phenomenon does not affect the accuracy of these results. After 2011, ecological restoration projects such as the "Inflow of the eastern ring water system", "Coordination of the South-North water transfer to divert river water" and "Ecological recharge of the Gang Huang reservoir to the Hutuo River" have been gradually carried out. The ecological environment of the Hutuo River basin has improved significantly since then, which has gradually changed from a hot spot area to a cold spot area. The cropland in the northwest of Zhengding County also showed a band of high-temperature areas. In fact, the cropland in this part is in the old Magnetic River, which was converted to cropland after the river dried up, forming an area with temperatures higher than the surrounding cropland. Some relative authorities have proposed to reduce the urban heat island effect by building the ecological corridor of the Old Magnetic River and the spatial structure of the ecological corridor of the South-North Water Diversion in the "Urban and Rural Master Plan of Zhengding County, Hebei Province (2014–2030)". In fact, the heat island effect of the Old Magnetic River has weakened in 2019 after a period of implementation, but it is still in the sub-high temperature area. Additionally, built-up areas in districts and counties always had higher temperature levels than the area of townships in these districts and counties.

#### *3.2. Spatial Dynamics of Urban Heat Island and Built-Up Land*

Figures 6 and 7 depict the schematic diagrams of the UHI and the built-up land SDE and center points in the whole study area from 1996 to 2019. The results suggest the center point migration trajectory of the heat island and built-up land within the study area was well coupled between 2001 and 2011, both of which tend to expand towards the northeast direction. However, the heat island center point leaned towards the northwest direction before 2001 and turned to the southwest direction after 2011, which was contrary to the build-up land during these periods. The more in-depth analysis results suggested there were strong heating effects in the Hutuo River Basin and the old channel of the Cihe River in Zhengding County, except for the heat island in the built-up area. In fact, the expansion of built-up land was very slow in the downtown before 2001, while the Hutuo River Basin and the old channel of the Cihe River had an increasingly powerful heat island effect. In other words, the heat island effect in the built-up area of the downtown is not as strong as that in the Hutuo River Basin and the old channel of the Cihe River in the northwest part of the study area, which made the heat island center point more prone to the northwest direction before 2001. By contrast, the spatial pattern of the heat island was consistent with the trend of the center point of built-up land expansion from 2001 to 2011. The heat island effect of built-up land expansion was more significant than that of the thermal environment effect in the Hutuo River Basin and the old channel of the Cihe River during this period. The built-up land expansion played a dominant role in influencing the spatial variation of the heat island, and therefore the center point of the heat island followed the evolution direction of the center point of built-up land. The center point of the heat island in 2011 was within the Gaocheng District on the north side of the Hutuo River. This is primarily due to the improvement of the eco-environment of the Hutuo River Basin, where the heat island effect was weakened and the cold island effect was enhanced, making the center point of the heat island shift to the south bank of Hutuo River. If the cold island effect in this region gets intensified as time goes by, the center point of the heat island may further shift towards the southwest direction.

**Figure 6.** Evolution of center points and standard deviation ellipse (SDE) of the urban heat island (UHI) of the study area from 1996 to 2019.

The results mentioned above suggested the ecological environmental transformation of the Hutuo River Basin had a significant influence on the evolution of the center point of the heat island along with built-up land expansion. It is, therefore, necessary to carry out more in-depth analyses of the areas strongly and weakly influenced by the thermal environment effects in the Hutuo River Basin. This study has therefore taken the downtown area and Yuhua District as an example to reveal the influence of the Hutuo River Basin on the shift characteristics of the heat island center. Figure 8 portrayed the schematic diagrams of the UHI, built-up land SDE, and center point in the downtown area from 1996 to 2019. The heat island center point in the downtown area tended to move in the same direction as the center point of built-up land in general from 2001 to 2011, while the heat island center point shifted in the direction opposite to the center point of built-up land from 2011 to 2019, which was generally consistent with the results of the whole study area.

**Figure 7.** Evolution of center points and standard deviation ellipse (SDE) of built-up land of the study area from 1996 to 2019.

**Figure 8.** Evolution of center points and standard deviation ellipse (SDE) of the built-up land and urban heat island (UHI) of the downtown from 1996 to 2019.

Figure 9 depicted the schematic diagram of the UHI and the SDE and center point of built-up land in the Yuhua District from 1996 to 2019. Both the built-up land expansion and the center point of the heat island migrated towards the southwest direction in Yuhua District, indicating the center shift of the heat island and built-up land kept consistent in areas less affected by the thermal environment effect in the Hutuo River Basin. These results further confirmed the thermal environment effect in the Hutuo River Basing, suggesting that the ecological restoration projects in the Hutuo River Basin since 2014 have effectively improved the regional eco-environment quality and significantly influenced the surrounding thermal environment. In summary, the heat island effect in the Hutuo River Basin affected the temperature of the surrounding area, causing the center points of built-up land and the heat island to shift away from the same direction within a certain spatial scope.

**Figure 9.** Evolution of center points and standard deviation ellipse (SDE) of the built-up land and urban heat island (UHI) of Yuhua District from 1996 to 2019.

#### *3.3. Results of Spatial Autocorrelation Analysis and Correlation Analysis*

The spatial autocorrelation analysis results showed that the *p*-values of the PLAND and LST were zero in all years, which passed the 95% confidence test. Besides this, the Z scores of both the PLAND and LST were at a high level (Table 3), indicating the stronger spatial autocorrelation of both the PLAND and LST. In addition, the Moran's I values in all years were above zero, and the Moran's I of the LST were higher than that of the PLAND in each period, indicating the positive spatial correlation of the surface temperature was stronger than that of the build-up. In particular, the spatial autocorrelation of PLAND and LST was both the strongest in 2015.

Figure 10 shows that there was a significant linear correlation between the PLAND and LST of the study area, which passed the significance test and proved the objective credibility of the results. Besides this, Figure 10 indicates that the samples are heavily clustered in the interval of the PLAND values of 0–60%. In addition, the Pearson correlation coefficients from 1996 to 2019 were close to 1, reflecting a high correlation level (Table 4).


**Table 3.** Moran's I of the percentage of built-up land (PLAND) and land surface temperature (LST) from 1996 to 2019.

**Figure 10.** Correlation analysis results of the percentage of built-up land (PLAND) and land surface temperature (LST) from 1996 to 2019.


**Table 4.** Correlation between the percentage of built-up land (PLAND) and land surface temperature (LST) from 1996 to 2019.

*3.4. Response Mechanism of the Percentage of Built-Up Land (PLAND) and Land Surface Temperature (LST)*

Figure 11 shows the response mechanism between the PLAND and LST. When the PLAND value is below 20–25%, the number of the sample grids decreased year by year, and when it is above this range, the number of the sample grids increased year by year. In particular, the number of sample grids in the 25–30% range tended to be the same, and the difference in the number of sample grids gradually increased as the value increased, indicating that the low-density built-up land in each grid gradually was transformed into high-density built-up land, and there was significant built-up land expansion.

**Figure 11.** Statistics of sample grids of different percentages of built-up land (PLAND) types.

This study categorized the grid in the study area according to the PLAND value, based on which the average temperature changes in the PLAND categories were revealed (Figure 12). The results showed that the temperature changed most slightly in the build-up land with the PLAND of 55–60%, and the 55–60% interval can be used to divide the sample grids of the whole study area into two categories. Specifically, for areas with the PLAND below the 55–60% range, the temperature increase tended to increase and then decrease as the PLAND increased. In particular, the maximum temperature was in the areas with the PLAND below 55–60%, and increases occurred in the areas with the PLAND of 25–30%. By contrast, the temperature increases in the areas with the PLAND above 55–60% showed first an increasing and then a decreasing trend as the PLAND increased. In particular, the maximum temperature increases in the areas with the PLAND above 55–60% occurred in the areas with the PLAND at 65–70%. The average surface temperature increase tended to slow down in the areas with the PLAND exceeding 70%, where the temperature still showed an increasing trend.

**Figure 12.** Statistics of the average temperature of different percentage of built-up land (PLAND) types.

This study further classified the sample grids into two categories based on the PLAND of 0–60% and 60–100% and carried out clustering analysis to further investigate the influence of the PLAND on the thermal environment. The LISA clustering diagram suggested that the high-high clusters were mainly in the built-up areas of the downtown and counties (Figure 13). The grids of high-high clusters gradually increased as time went by, reflecting the increasing heat island effects due to the built-up land expansion. Meanwhile, the low-high clusters were scattered at the edge of the built-up areas of the downtown, and the high-low clusters were distributed in the suburban areas of each district and county. All these results can also provide a firm basis for the boundary division of built-up areas.

**Figure 13.** LISA cluster map of impacts of the percentage of built-up land (PLAND) on land surface temperature (LST) from 1996 to 2011.

#### *3.5. Discussion*

This study revealed the thermal environmental effects of the built-up land expansion based on long time series data, the results of which can provide valuable references for urban development and ecological environment improvement, but it is still necessary to carry out some more in-depth research. Firstly, the study period of this study lasted 24 years, but the time interval is not uniform, and it is necessary to carry out more in-depth case studies in more sample cities. Besides this, it is necessary to link the research results with the relevant urban planning in order to effectively improve the thermal environment in urban areas. For example, the results of this study suggested it is necessary for the relevant departments to expand the green area of urban areas as vastly as possible and accelerate the ecological restoration projects of the Hutuo River Basin in order to enhance its cold island effect and thus reduce the adverse thermal environment effects in the surrounding areas. Additionally, it is also necessary for relevant departments to strengthen the control of urban size in order to avoid further deterioration of the thermal environment resulting from the uncontrolled built-up land expansion according to the results of this study.

#### **4. Conclusions**

This study explored the thermal environment effects of built-up land expansion in Shijiazhuang from 1996 to 2019 based on land-use data and temperature products and revealed the response mechanism between the PLAND and LST in six districts and two counties of Shijiazhuang. The results showed that: (1) The heat island in the study area was mainly distributed in the downtown, built-up areas of each district and county, and the Hutuo River Basin. However, ecological restoration projects have remarkably improved the regional ecological environment in the Hutuo River Basin, where the heat island was weakened and significantly influenced the evolution of the center point shift of the heat island in the study area. (2) There was a significant correlation between the PLAND and the LST, with the correlation coefficients all above 0.7 during the study period, indicating the significant role of built-up land expansion in affecting the thermal environment. (3) PLAND = 55–60% separated the surface thermal environmental effects of built-up land expansion in urban and rural areas, which reflects the heat island effect of the urban built-up land expansion and the cold island effect of rural cropland. The areas with PLAND > 60% are generally urban regions with stronger heat island effects, while areas with PLAND < 55% are generally villages and towns where the temperature raised more slowly. These results of this study can provide a valuable reference for urban development and ecological environment improvement in Shijiazhuang and other cities, but it is necessary to link the research results with the relevant urban planning in order to effectively improve the thermal environment in urban areas.

**Author Contributions:** L.Q. and H.L.: Investigation, data curation, writing—original draft preparation, writing—review and editing, funding acquisition; G.S. and H.Y. (Haiming Yan): conceptualization, methodology, supervision, project administration; H.Y. (Huicai Yang): software, validation, visualization. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the Social Science Foundation of Hebei Province (HB17YJ026).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

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

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

