Simulation of an Urban-Rural Spatial Structure on the Basis of Green Infrastructure Assessment: The Case of Harbin, China

Due to their long-term dual structures and rapid urbanization, cities and villages in developing countries are undergoing the challenges of urban-rural integration and ecological security. This study aims to determine the pattern of urban-rural spatial structures under the circumstances of ecological security in the future to promote the integrated, coordinated, green, and sustainable development of urban-rural spaces. Using a quantitative evaluation method, the logistic-CA model, the LCP (least cost path) model, and a classification of ecological importance, this study constructed an integrated simulation model based on a green infrastructure assessment and applied the model to simulate and predict the urban-rural spatial structure of the Harbin city territory (Harbin) in 2035. The results indicate that the urban-rural hierarchical scale structure of Harbin comprises a central city, sub-central city, central town, major town, common town, central village, and general village. The urban-rural traffic network structure forms a pattern of “radiation + grid”, with Harbin city at the center of the structure. The urban-rural land use zoning structure consists of eco-spaces, agricultural spaces, and construction spaces. It can be concluded that in 2035, the urban-rural spatial structure of Harbin will show an increasing development tendency, where single-center, medium, and small cities in will Harbin develop, and traffic systems above the county level will also improve.


Introduction
The cities and villages of developed countries, like Britain, France, the USA, and Japan, have broadly evolved into a state of balance [1,2]. However, influenced by the extreme differences between urban and rural areas [3,4], the contradiction of urban-rural development in China remains prominent. Urban and rural areas are equally essential and should be combined organically [5]. Therefore, how to intensively make use of urban-rural resources and spaces, rationally optimize the overall layout of urban and rural areas, protect and rectify the urban-rural ecological environment, and coordinate the relationship between urban and rural development to promote the sustainable and healthy development of urban-rural society and economy, is one of main tasks confronting China at the moment. The urban-rural spatial structure, a regional reflection of humans' social, economic, and cultural activities under specific circumstances, illustrates the spatial dichotomy between urban and rural areas. In addition, with the rapid development of urbanization, urban sprawl has brought great challenges to urban and rural eco-environments [6][7][8], especially in China [9,10]. How to protect and maintain ecosystem security during the process of urban growth, to resolve the contradiction between Land 2019, 8, 196 3 of 21 This study offers two contributions to knowledge: (1) This study attempts to construct an integrated simulation model to simulate urban-rural spatial structures; and (2) this study offers an analysis of Harbin's urban-rural spatial structure on a regional scale from an ecological perspective. Previous research on urban spatial simulations have focused on central Harbin city [31] or Harbin city [32]. The scope of this study will be the administrative area of Harbin, that is the Harbin city territory. This territory contains Harbin city and nine counties (or county-level cities). Research on the spatial structure of Harbin aims to reveal the spatial expansion and evolution of the construction land in central Harbin city [33][34][35] without considering rural areas and ecological problems. Additionally, some studies have placed Harbin city into a bigger scale and broader context, analyzing the urban agglomeration's spatial structure evolution [36] and its patterns [37].
The main purpose of the study was to understand the patterns of the urban-rural spatial structure under ecological security via the constructed simulation model to propose optimization and implementation strategies in the next step. This study also aims to accelerate harmonious urban and rural development through overall consideration of urban-rural settlements, traffic systems, and land use zoning, thus further promoting the achievement of the final target of urban-rural integration. This is a significant and key aspect of Chinese national spatial planning.

Study Area and Data
Harbin city territory (i.e., the administrative area of Harbin), lies in the northeastern part of China ( Figure 1) and is the capital city of Heilongjiang Province, which contains Harbin city and nine counties (or county-level cities). It has a total administrative area of 53,100 km 2 , a city area of 10,198 km 2 , and a central city area of 4187 km 2 . The total population in the administrative area was above 9.61 million by the end of 2015 (data from the 'Harbin Statistics Yearbook 2016'). The urban-rural spatial structure for 2015 is shown in Figure 2. The Harbin urban-rural residential system comprises the central city (1), the sub-central city (2), the central town (14), the major town (58), the common town (111), the central village (1689), and the general village. The traffic network structure consists of a railway network and a highway network. The density of the traffic network above the county level is 0.39 km/km 2 . This study offers two contributions to knowledge: (1) This study attempts to construct an integrated simulation model to simulate urban-rural spatial structures; and (2) this study offers an analysis of Harbin's urban-rural spatial structure on a regional scale from an ecological perspective. Previous research on urban spatial simulations have focused on central Harbin city [31] or Harbin city [32]. The scope of this study will be the administrative area of Harbin, that is the Harbin city territory. This territory contains Harbin city and nine counties (or county-level cities). Research on the spatial structure of Harbin aims to reveal the spatial expansion and evolution of the construction land in central Harbin city [33][34][35] without considering rural areas and ecological problems. Additionally, some studies have placed Harbin city into a bigger scale and broader context, analyzing the urban agglomeration's spatial structure evolution [36] and its patterns [37].
The main purpose of the study was to understand the patterns of the urban-rural spatial structure under ecological security via the constructed simulation model to propose optimization and implementation strategies in the next step. This study also aims to accelerate harmonious urban and rural development through overall consideration of urban-rural settlements, traffic systems, and land use zoning, thus further promoting the achievement of the final target of urban-rural integration. This is a significant and key aspect of Chinese national spatial planning.

Study Area and Data
Harbin city territory (i.e., the administrative area of Harbin), lies in the northeastern part of China ( Figure 1) and is the capital city of Heilongjiang Province, which contains Harbin city and nine counties (or county-level cities). It has a total administrative area of 53,100 km 2 , a city area of 10,198 km 2 , and a central city area of 4187 km 2 . The total population in the administrative area was above 9.61 million by the end of 2015 (data from the 'Harbin Statistics Yearbook 2016'). The urban-rural spatial structure for 2015 is shown in Figure 2. The Harbin urban-rural residential system comprises the central city (1), the sub-central city (2), the central town (14), the major town (58), the common town (111), the central village (1689), and the general village. The traffic network structure consists of a railway network and a highway network. The density of the traffic network above the county level is 0.39 km/km 2 .   The urban-rural development of Harbin city territory has a long history dating back to the first appearance of the five Ancient Cities during the Jin Dynasty (1115-1234 BC). Compared with other cities in Heilongjiang Province, the level of urban-rural integration in Harbin city is high [38], but compared to the whole of China, its level is relatively low [39]. Furthermore, the urban-rural integration between Harbin city and the other counties shows a development imbalance [40]. In order to enhance the overall urban-rural integration, this study, taking the Harbin city territory as the object of the study, simulates the urban-rural spatial structure of ecological priority to develop optimization strategies. A flowchart of the study structure is shown in Figure 3. For convenience, the abbreviation "Harbin" hereafter represents "Harbin city territory".
Land use data used in this study were offered by Beijing Digital Technology Co. Ltd. Digital Elevation Model (DEM) and Moderate-resolution Imaging Spectroradiometer (MODIS) data were downloaded separately from the Geospatial Data Cloud [41] and National Aeronautics and Space Administration (NASA) [42]. Traffic, demographic, and Gross Domestic Product (GDP) data were taken from the 'Harbin Yearbook', 'Harbin Statistics Yearbook', and relevant government departments. All data were processed using ArcGIS.

Quantitative Evaluation Method
Green infrastructure refers to a network of natural and semi-natural green and blue areas that can maintain and improve ecosystem services, which comprises different kinds of green vegetation The urban-rural development of Harbin city territory has a long history dating back to the first appearance of the five Ancient Cities during the Jin Dynasty (1115-1234 BC). Compared with other cities in Heilongjiang Province, the level of urban-rural integration in Harbin city is high [38], but compared to the whole of China, its level is relatively low [39]. Furthermore, the urban-rural integration between Harbin city and the other counties shows a development imbalance [40]. In order to enhance the overall urban-rural integration, this study, taking the Harbin city territory as the object of the study, simulates the urban-rural spatial structure of ecological priority to develop optimization strategies. A flowchart of the study structure is shown in Figure 3. For convenience, the abbreviation "Harbin" hereafter represents "Harbin city territory". The urban-rural development of Harbin city territory has a long history dating back to the first appearance of the five Ancient Cities during the Jin Dynasty (1115-1234 BC). Compared with other cities in Heilongjiang Province, the level of urban-rural integration in Harbin city is high [38], but compared to the whole of China, its level is relatively low [39]. Furthermore, the urban-rural integration between Harbin city and the other counties shows a development imbalance [40]. In order to enhance the overall urban-rural integration, this study, taking the Harbin city territory as the object of the study, simulates the urban-rural spatial structure of ecological priority to develop optimization strategies. A flowchart of the study structure is shown in Figure 3. For convenience, the abbreviation "Harbin" hereafter represents "Harbin city territory".
Land use data used in this study were offered by Beijing Digital Technology Co. Ltd. Digital Elevation Model (DEM) and Moderate-resolution Imaging Spectroradiometer (MODIS) data were downloaded separately from the Geospatial Data Cloud [41] and National Aeronautics and Space Administration (NASA) [42]. Traffic, demographic, and Gross Domestic Product (GDP) data were taken from the 'Harbin Yearbook', 'Harbin Statistics Yearbook', and relevant government departments. All data were processed using ArcGIS.

Quantitative Evaluation Method
Green infrastructure refers to a network of natural and semi-natural green and blue areas that can maintain and improve ecosystem services, which comprises different kinds of green vegetation Land use data used in this study were offered by Beijing Digital Technology Co. Ltd. Digital Elevation Model (DEM) and Moderate-resolution Imaging Spectroradiometer (MODIS) data were downloaded separately from the Geospatial Data Cloud [41] and National Aeronautics and Space Administration (NASA) [42]. Traffic, demographic, and Gross Domestic Product (GDP) data were taken from the 'Harbin Yearbook', 'Harbin Statistics Yearbook', and relevant government departments. All data were processed using ArcGIS.

Quantitative Evaluation Method
Green infrastructure refers to a network of natural and semi-natural green and blue areas that can maintain and improve ecosystem services, which comprises different kinds of green vegetation and water bodies [20][21][22]. Assessing green infrastructure with a quantitative evaluation method aims to understand the distribution of key ecological resources and sequence the ecological importance of these resources. The results will be used as restrictive conditions in the integrated simulation model. Therefore, an evaluation index system was built from resource protection and ecological services. According to the data availability, the scientific principle of index selection, and previous research [43][44][45][46][47], this study selected 12 indices to assess Harbin's green infrastructure. Specific indicators can be seen in Table 1. The distance to each resource refers to the distance of the resource boundaries from a certain point, which can be measured by the buffer tool in ArcGIS. Generally, the closer the distance to each resource, the higher ecological sensitivity the point has during the process of land development. Thus, the ecological importance of this point is also higher. The distance of the buffer in this study was determined by multiple empirical studies [47][48][49]. The indices (e.g., the capability of support, the ability for temperature reduction, etc.) were measured by the product of the equivalents of the ecosystem service value supplied by the per unit area of the ecosystem [50] and the areas of various green infrastructures. The weights of the indicators were determined by expert marking and literature research on the network construction, planning, and evaluation of green infrastructure [45][46][47]51,52]. Table 1. The evaluation index system of green infrastructure 1 .

Categories Factors Indices (X ij ) Weights
Resources protection (X 1 ) Drinking water sources The distance to water sources (X 11 ) 0.10 Great parks (>= 100 km 2 ) The distance to parks (X 12 ) 0.07 Nature reserves The distance to nature reserves (X 13  The expert marking method refers to consulting experts in the fields of landscape, ecology, and planning in the form of a questionnaire to obtain the importance value (X ij ) of each index. Table 2 is an example of the expert marking table. A new score (X ij ) is reassigned to the average score (X ij ) of 10 experts according to Table 3. This new value is considered to be the final importance value for each indicator. Then, the weight (w ij ) of each indicator can be calculated by Equation (1). Meanwhile, these 10 experts are invited to score the relative importance between the resource protection factors and the ecological service factors to obtain the combined weight (w i ). Consequently, the final weight (w ij ) of each indicator is the result of Equation (2). The overlay process is fulfilled in ArcGIS.
where n is the number of indicators in the i category.

Prediction of the Urban-Rural Construction Land Scale
In general, the increase in construction land originates from the requirements for residential, public service, and industrial land due to population growth [53,54]. Therefore, this study built a linear regression equation between population, GDP, and construction land, but the results showed that this equation could not pass the DW (Durbin-Watson) test. According to the 'Harbin Statistics Yearbook 2018', the population of Harbin has experienced a downward trend since 2013. As a result, this study supposes that industrial production brought about the expansion of construction land, so the correlation between construction land and per capita GDP was analyzed and the correlation coefficient reached 0.9978. Meanwhile, the per capita GDP increased over time, and the two factors showed a linear correlation. Thus, using data from 2000 to 2016, first, the per capita GDP for the year 2035 was projected by linear regression Equation (3), and then the scale of the construction land was calculated in accordance with linear regression Equation (4). In addition, the linear regression equation of the construction land and time was also established [55] (Equation (5)). All equations passed the R test and the DW test.
where y 1 is GDP; y 2 and y 3 are the scale of urban-rural construction land; and t is time (year). Finally, the average of y 2 and y 3 was determined as the scale of urban-rural construction land in 2035.

Logistic-Celluar Aotomata Model
Geographical CA is a commonly used model in simulating urban land development [56,57] and comprises four basic parameters: cells, state, neighborhood, and transformation rules. The state of each cell at time t + 1 is decided by the state of the cell itself and its surroundings and the transformation rules at time t. This can be expressed by Equation (6).
where S t+1 ij and S t ij are the land use states of cell ij at time t + 1 and t, respectively; N t ij is the state set of cells around cell ij; and R represents the transformation rule. With the development of the CA model, the constrained urban CA model is able to mimic changes in the urban patterns of different developing targets or scenarios [25,27,58]. For instance, if ecological conditions are added to a CA model, then this where con S t ij expresses a constrained condition set that decides whether the state of cell ij at t moment can be transformed; its value is 0~1.
According to Equation (7), adopting an 8-neighborhood rule and transformation rule obtained by the logistic regression equation method [59], taking the results of green infrastructure assessment as an ecological constrained condition, this study simulates the evolution of urban-rural construction land. On the basis of the simulation results, the size of the urban and rural construction land was computed. In accordance with the computed size, urban and rural settlements in a specific region can be divided into the central city, sub-central city, central town, major town, common town, etc., which were combined into the hierarchical size structure [29,60]. Generally, the determination of a city's categories is based on population size, which was decided by the size of the construction land in this study due to the limitations of data acquisition. Therefore, the "size" in the hierarchical size structure of this study was the size of the construction land.

Least Cost Path Model
The LCP model reflects the cost consumed or the work done when organisms move from the sources to the destinations. As it features a concise data structure, a fast algorithm, and intuitive results, LCP is widely used to study route simulations and planning problems [61][62][63][64][65]. This study aims to simulate the optimal travel path between urban-rural settlements on the basis of ecological protection. Consequently, an LCP model was selected. In order to apply the LCP model, three elements (sources, resistance factors, and resistance) need to be determined in advance. Sources are the origins of species dispersal and maintenance, which are usually composed of protected species or existing habitats. In this study, sources refer to urban-rural settlements that support the agglomeration and diffusion of the population and economic elements. The resistance factors vary based on the research objectives. Construction suitability, ecological condition, and land use types constitute the resistance factors of this study. The resistance values of the factors are different in terms of the subjects, which are usually decided upon by experts and empirical values [49,63], but these values only represent the relative resistance for organisms to migrate and disperse. Taking the setting of the resistance value for land use types as an example, in the field of ecology, the resistance of construction land is normally the highest, while the resistances of natural reserves and forests are relatively lower in terms of the migration of species, which is the opposite in this study because the resistance value of this LCP model refers to the cost of the movement of residents or material flows. Finally, the resistance factors and their resistance values were overlaid to obtain the resistance surface, which represents the difficulty of flowing.

Classification of Ecological Importance
The purpose of dividing land by its level of ecological importance is to provide instruction for making land protection and development measures through understanding the geographical distribution of land with its eco-value and, ultimately, to realize urban-rural sustainable development. In this study, the boundaries and ranks of ecological spaces can be determined according to the classification of ecological importance, which is the basis of dividing urban-rural functional spaces and determining controlling measures.
This study estimates and sequences the ecological importance of urban-rural land on the basis of green infrastructure assessment and adopts a natural break method to rank importance. Based on the obtained results of ecological importance zoning, coupled with the scope of the urban-rural construction land and farmland, the urban-rural land use zoning structure can be determined.

Harbin Green Infrastructure Assessment
According to the definition of green infrastructure used in this study, Harbin's green infrastructure consists of several land use types including farmland, forest, grassland, water bodies, and other land [66]. Applying the established evaluation index system to the green infrastructure assessment of Harbin in 2015, Figure 4 describes the sequencing results of the ecological importance of Harbin's green infrastructure. In general, the green infrastructure with high ecological importance is distributed in the middle region of Harbin, while those with lower ecological importance are primarily concentrated in Harbin city. Natural reserves, great parks, etc. that need to be protected, gain the highest value, which means that strategies of priority protection should be applied to these areas. It should be noted that the purple parts in Figure 4 are urban-rural construction areas, which do not belong to the green infrastructure. Consequently, in the subsequent study on urban-rural spatial structure simulation, areas with higher ecological importance are regarded as limited development areas, while urban-rural construction areas are considered suitable development areas, where the value of importance is 0. This would benefit eco-system security and sustainable land development in Harbin. According to the definition of green infrastructure used in this study, Harbin's green infrastructure consists of several land use types including farmland, forest, grassland, water bodies, and other land [66]. Applying the established evaluation index system to the green infrastructure assessment of Harbin in 2015, Figure 4 describes the sequencing results of the ecological importance of Harbin's green infrastructure. In general, the green infrastructure with high ecological importance is distributed in the middle region of Harbin, while those with lower ecological importance are primarily concentrated in Harbin city. Natural reserves, great parks, etc. that need to be protected, gain the highest value, which means that strategies of priority protection should be applied to these areas. It should be noted that the purple parts in Figure 4 are urban-rural construction areas, which do not belong to the green infrastructure. Consequently, in the subsequent study on urban-rural spatial structure simulation, areas with higher ecological importance are regarded as limited development areas, while urban-rural construction areas are considered suitable development areas, where the value of importance is 0. This would benefit eco-system security and sustainable land development in Harbin.

Hierarchical Size Structure Extraction
This study adopted a land scale instead of a population scale to discuss the urban-rural hierarchical size structure of Harbin in 2035. Therefore, a logistic-CA model was used to simulate the expansion of urban-rural construction land. On the basis of this simulation, the urban-rural hierarchical size structure can be extracted.

Parameter Settings and Accuracy Test
Before applying a logistic-CA model, it is necessary to set parameters and examine the accuracy of the model in advance. The resolution of the raster data used in the model in this study was 90 m x 90 m. The parameter settings of the model are as follows: • Spatial variables of the transformation rule (R) based on logistic regression: distance to the highway, distance to the main road, distance to the urban center, distance to the town center, per capita GDP, population density, elevation, and slope;

Hierarchical Size Structure Extraction
This study adopted a land scale instead of a population scale to discuss the urban-rural hierarchical size structure of Harbin in 2035. Therefore, a logistic-CA model was used to simulate the expansion of urban-rural construction land. On the basis of this simulation, the urban-rural hierarchical size structure can be extracted.

Parameter Settings and Accuracy Test
Before applying a logistic-CA model, it is necessary to set parameters and examine the accuracy of the model in advance. The resolution of the raster data used in the model in this study was 90 m x 90 m. The parameter settings of the model are as follows: • Spatial variables of the transformation rule (R) based on logistic regression: distance to the highway, distance to the main road, distance to the urban center, distance to the town center, per capita GDP, population density, elevation, and slope; • Neighborhood rule (N t ij ): a neighborhood space of 3 x 3 grid cells; • Constrained condition (con S t ij ): the results of green infrastructure assessment. Using the land use data of Harbin from 2005, 2013, and 2015, the transformation rule and the accuracy of the logistic-CA model was obtained and tested. In SPSS, taking two land use maps as the inputs of the dependent variable, and the distance to the highway (X 1 ), distance to the main road (X 2 ), distance to the urban center (X 3 ), distance to the town center (X 4 ), per capita GDP (X 5 ), population density (X 6 ), elevation (X 7 ), and slope (X 8 ) as the independent variables, the logistic regression Equation (8) was established: The slope and elevation factors have the most obvious influence on the development of urban-rural construction land. Comparing the simulation results of 2015 (Figure 5a) with the real urban-rural construction land development in 2015 (Figure 5b), and adopting a pixel by pixel comparison method, the results indicate that the overall accuracy value of the model reached 96.88%. The Kappa simulation, which refers to the coefficient of agreement between the simulated land use transitions and the observed land use transitions [67], was 0.82. In addition, the Moran I index, which can reflect the extent of the spatial concentration or dispersion [68], was also calculated ( Table 4). The results show that the simulated Moran I index value was lower than the real value of 2015 as well as the value of 2005. Thus, it can be seen that the land use pattern of Harbin tends to disperse evolution. On the whole, the simulated results from the model were close to the observed development of construction land. Consequently, this model can be used to predict land use changes in the future.  (Figure 5b), and adopting a pixel by pixel comparison method, the results indicate that the overall accuracy value of the model reached 96.88%. The Kappa simulation, which refers to the coefficient of agreement between the simulated land use transitions and the observed land use transitions [67], was 0.82. In addition, the Moran I index, which can reflect the extent of the spatial concentration or dispersion [68], was also calculated ( Table 4). The results show that the simulated Moran I index value was lower than the real value of 2015 as well as the value of 2005. Thus, it can be seen that the land use pattern of Harbin tends to disperse evolution. On the whole, the simulated results from the model were close to the observed development of construction land. Consequently, this model can be used to predict land use changes in the future.   (4) and (5), is 2341.06 km 2 and 2320.75 km 2 , respectively. Therefore, the mean of 2330.91 km 2 is considered as the urban-rural construction land scale for 2035.
Based on the predicted scale of 2035, the results of the green infrastructure assessment, and the tested logistic-CA model, the land use change of Harbin in 2035 was obtained ( Figure 6). It was found that the Moran I index for 2035 declines, which means that the land use pattern in Harbin is more dispersed. Appropriate land use adjustment measures should be taken to enhance the intensive   (4) and (5), is 2341.06 km 2 and 2320.75 km 2 , respectively. Therefore, the mean of 2330.91 km 2 is considered as the urban-rural construction land scale for 2035.
Based on the predicted scale of 2035, the results of the green infrastructure assessment, and the tested logistic-CA model, the land use change of Harbin in 2035 was obtained ( Figure 6). It was found that the Moran I index for 2035 declines, which means that the land use pattern in Harbin is more dispersed. Appropriate land use adjustment measures should be taken to enhance the intensive utilization of urban-rural land. Subsequently, the urban-rural hierarchical size structure for Harbin in 2035 was extracted based on the land use of 2035. As shown in Figure 7, Harbin's urban-rural residential system comprises a central city (1), sub-central city (5), central town (12), major town (72), common town (94), central village (2301), and general village (28,665). Out of these categories, because Harbin city has a land area of more than 400 km 2 , it was classed as the central city; cities with land areas between 10 and 20 km 2 including Wuchang, Shangzhi, Yanshou, Tonghe, and Binzhou would be classified as sub-central cities; the land areas of Yilan, Bayan, Mulan, Fangzheng, Dalianhe, Yimianpo, Weihe, Yabuli, Lalin, Binxi, Xinglong, and Shanhe towns range from 5 to 10 km 2 , which would make them central towns; towns with a scale of construction land from 1 to 5 km 2 are major towns; other towns whose areas are smaller than 1 km 2 are common towns; villages with a land area between 0.2 and 1 km 2 are considered central villages; other villages are general villages. By comparing the hierarchical size structure for 2035 with that of 2015, the number of sub-central cities and major towns increases, but there is a decrease in the number of central towns, which occurs due to the upgrading of the Wuchang, Yanshou, and Binzhou towns, which indicates that the medium and small cities of Harbin will develop. in 2035 was extracted based on the land use of 2035. As shown in Figure 7, Harbin's urban-rural residential system comprises a central city (1), sub-central city (5), central town (12), major town (72), common town (94), central village (2301), and general village (28,665). Out of these categories, because Harbin city has a land area of more than 400 km 2 , it was classed as the central city; cities with land areas between 10 and 20 km 2 including Wuchang, Shangzhi, Yanshou, Tonghe, and Binzhou would be classified as sub-central cities; the land areas of Yilan, Bayan, Mulan, Fangzheng, Dalianhe, Yimianpo, Weihe, Yabuli, Lalin, Binxi, Xinglong, and Shanhe towns range from 5 to 10 km 2 , which would make them central towns; towns with a scale of construction land from 1 to 5 km 2 are major towns; other towns whose areas are smaller than 1 km 2 are common towns; villages with a land area between 0.2 and 1 km 2 are considered central villages; other villages are general villages. By comparing the hierarchical size structure for 2035 with that of 2015, the number of sub-central cities and major towns increases, but there is a decrease in the number of central towns, which occurs due to the upgrading of the Wuchang, Yanshou, and Binzhou towns, which indicates that the medium and small cities of Harbin will develop.

Traffic Network Structure Simulation
With the help of ArcGIS, based on the extracted hierarchical size structure, an LCP model was applied to simulate potential traffic routes. Harbin's urban-rural traffic network system, based on the green infrastructure assessment used in this study, was obtained via the extraction of key traffic corridors and integration with current and forecasted traffic.

Simulation of Potential Traffic Routes
The simulation of potential traffic routes can be divided into three steps: source identification, building a resistant surface, and traffic network construction. In this study, urban-rural settlements were regarded as sources. The building of a resistant surface was implemented by using an overlay tool in ArcGIS, and the resistance factors and their resistance can be seen in Table 5. To set the range, the range of the green infrastructure assessment results was divided according to a manual classification method in ArcGIS, and the range of the slope was determined based on the suitability of the engineering construction [69]. The setting of resistance references was done according to previous studies [47,49,63]. The resistance value is a relative value, which only expresses the relative resistance of the landscape by a numerical value. Finally, with the aid of cost distance and cost path tools, the potential urban-rural traffic network (Figure 8) of Harbin for 2035, based on green infrastructure assessment, was mapped.

Construction of Traffic Network
According to the importance of the highway, the grading of roads can be divided into five grades: motorway, national road, provincial road, county road, and rural road. On the basis of the 'National Highway Network Planning', and the 'Heilongjiang Provincial Highway Network Planning', 'Comprehensive Transportation Planning of Harbin during the 13 th Five-Year Plan' proposed that the layout of the motorway system is "two rings, nine axes, and two links" (Figure 10a), the layout of the common trunk highway system is "one ring, thirteen axes, three vertical lines, and eleven links" (Figure 10b). In accordance with this transportation planning, coupled with Figures 8 and 9, and the existing roads, Harbin's urban-rural traffic network structure in 2035, based on a green infrastructure assessment, was constructed by adjusting the grade of some simulated roads. As shown in Figure 11, this structure consists of a railway network and a highway network. In this structure, the highway network contains a motorway system, a common trunk highway system composed of national roads, provincial roads, and county roads, and a rural highway system. The density of the urban-rural traffic network excluding the rural highway system will have increased by 0.08 km/km 2 since 2015, when the density was 0.39 km/km 2 . This increase shows that traffic systems above the county level have improved.

Extraction of Key Traffic Corridors
In accordance with the grade of the urban-rural settlements connected by traffic routes, key traffic corridors were extracted. Settlements with a higher grade in their hierarchical size structure are more important, and the traffic routes that link those settlements are more crucial. As a result, this study considered the potential routes that connected the central city, sub-central cities, and central towns as key traffic paths, and those that link major towns are deemed vice-key traffic paths (Figure 9).

Construction of Traffic Network
According to the importance of the highway, the grading of roads can be divided into five grades: motorway, national road, provincial road, county road, and rural road. On the basis of the 'National Highway Network Planning', and the 'Heilongjiang Provincial Highway Network Planning', 'Comprehensive Transportation Planning of Harbin during the 13 th Five-Year Plan' proposed that the layout of the motorway system is "two rings, nine axes, and two links" (Figure 10a), the layout of the common trunk highway system is "one ring, thirteen axes, three vertical lines, and eleven links" (Figure 10b). In accordance with this transportation planning, coupled with Figures 8 and 9, and the existing roads, Harbin's urban-rural traffic network structure in 2035, based on a green infrastructure assessment, was constructed by adjusting the grade of some simulated roads. As shown in Figure 11, this structure consists of a railway network and a highway network. In this structure, the highway network contains a motorway system, a common trunk highway system composed of national roads,

Construction of Traffic Network
According to the importance of the highway, the grading of roads can be divided into five grades: motorway, national road, provincial road, county road, and rural road. On the basis of the 'National Highway Network Planning', and the 'Heilongjiang Provincial Highway Network Planning', 'Comprehensive Transportation Planning of Harbin during the 13 th Five-Year Plan' proposed that the layout of the motorway system is "two rings, nine axes, and two links" (Figure 10a), the layout of the common trunk highway system is "one ring, thirteen axes, three vertical lines, and eleven links" (Figure 10b). In accordance with this transportation planning, coupled with Figures 8 and 9, and the existing roads, Harbin's urban-rural traffic network structure in 2035, based on a green infrastructure assessment, was constructed by adjusting the grade of some simulated roads. As shown in Figure 11, this structure consists of a railway network and a highway network. In this structure, the highway network contains a motorway system, a common trunk highway system composed of national roads, provincial roads, and county roads, and a rural highway system. The density of the urban-rural traffic network excluding the rural highway system will have increased by 0.08 km/km 2 since 2015, when the density was 0.39 km/km 2 . This increase shows that traffic systems above the county level have improved.

Land Use Zoning Structure Determination
In urban-rural spatial planning, the division of urban-rural functional spaces on the basis of ecological importance zoning is an essential way to put the idea of ecological priority into practice. This study divides the urban-rural spaces of Harbin into urban-rural construction spaces, ecological spaces, and agricultural spaces based on the ecological importance zoning obtained via a classification of ecological importance. In this way, the land use zoning structure is determined.

Ecological Importance Zoning
According to the sequencing of ecological importance from the green infrastructure assessment, and taking a natural break classification approach [70,71] (natural break values in this approach is

Land Use Zoning Structure Determination
In urban-rural spatial planning, the division of urban-rural functional spaces on the basis of ecological importance zoning is an essential way to put the idea of ecological priority into practice. This study divides the urban-rural spaces of Harbin into urban-rural construction spaces, ecological spaces, and agricultural spaces based on the ecological importance zoning obtained via a classification of ecological importance. In this way, the land use zoning structure is determined.

Ecological Importance Zoning
According to the sequencing of ecological importance from the green infrastructure assessment, and taking a natural break classification approach [70,71] (natural break values in this approach is obtained by the Jenk optimization method of statistics, which can minimize the sum of internal

Land Use Zoning Structure Determination
In urban-rural spatial planning, the division of urban-rural functional spaces on the basis of ecological importance zoning is an essential way to put the idea of ecological priority into practice. This study divides the urban-rural spaces of Harbin into urban-rural construction spaces, ecological spaces, and agricultural spaces based on the ecological importance zoning obtained via a classification of ecological importance. In this way, the land use zoning structure is determined.

Ecological Importance Zoning
According to the sequencing of ecological importance from the green infrastructure assessment, and taking a natural break classification approach [70,71] (natural break values in this approach is obtained by the Jenk optimization method of statistics, which can minimize the sum of internal variances at all levels [72]) in GIS, Harbin (Figure 12) can be divided into extremely important (5.335-7.28), very important (3.758-5.335), relatively important (2.826-3.758), generally important (1.135-2.826), and not important (0-1.135) zones. Extremely important and very important zones are regarded as excluded from development areas, where protection and restoration behavior are conducted. Relatively important and generally important zones are taken as restricted areas, in which measures of protection priority will be taken. Not important zones contain all the construction land that can be developed and can be classified as areas suitable for construction.
In accordance with the theory of urban-rural integration and the definition of territorial functional space, this study classifies the functional spaces of Harbin as urban-rural construction spaces, eco-spaces, and agricultural spaces ( Figure 13). Eco-spaces are divided into control regions of ecological red lines and other eco-spaces. The boundary lines of excluded from development areas with high eco-values are ecological red lines. Regions within the scope of the ecological red lines are classified as control regions that should be strictly protected to keep the urban-rural eco-system safe. The scopes of urban-rural construction land, as determined by the land use simulation for 2035, are identified as the urban-rural construction spaces of 2035. Comparing the current ecological importance zoning picture with the predicted land use change picture of 2035, agricultural spaces and other eco-spaces located in restricted areas were defined. In this study, agricultural spaces refer to agricultural production spaces, but do not include rural living spaces. Other eco-spaces mainly comprise grasslands and other types of land. The percentages of each kind of space in the total land area are as follows: eco-space, 51.4% (which controls regions of ecological red lines comprising 46.7%, while other eco-spaces comprise 4.7%); urban-rural space, 4.4%; and agricultural space, 44.0%. Construction spaces and agricultural spaces are primarily distributed in the smoother western areas, while eco-spaces are mostly concentrated in the middle area, where the terrain is higher, and the surface is covered by forests and rivers.

Delimitation of Urban-Rural Functional Spaces
In accordance with the theory of urban-rural integration and the definition of territorial functional space, this study classifies the functional spaces of Harbin as urban-rural construction spaces, eco-spaces, and agricultural spaces ( Figure 13). Eco-spaces are divided into control regions of ecological red lines and other eco-spaces. The boundary lines of excluded from development areas with high eco-values are ecological red lines. Regions within the scope of the ecological red lines are classified as control regions that should be strictly protected to keep the urban-rural eco-system safe. The scopes of urban-rural construction land, as determined by the land use simulation for 2035, are identified as the urban-rural construction spaces of 2035. Comparing the current ecological importance zoning picture with the predicted land use change picture of 2035, agricultural spaces and other eco-spaces located in restricted areas were defined. In this study, agricultural spaces refer to agricultural production spaces, but do not include rural living spaces. Other eco-spaces mainly comprise grasslands and other types of land. The percentages of each kind of space in the total land area are as follows: eco-space, 51.4% (which controls regions of ecological red lines comprising 46.7%, while other eco-spaces comprise 4.7%); urban-rural space, 4.4%; and agricultural space, 44.0%. Construction spaces and agricultural spaces are primarily distributed in the smoother western areas, while eco-spaces are mostly concentrated in the middle area, where the terrain is higher, and the surface is covered by forests and rivers.

Integration of Urban-rural Spatial Structures
The obtained urban-rural hierarchical size structure, traffic network structure, and land use zoning structure were integrated to obtain the overall urban-rural spatial structure of Harbin in 2035 based on the green infrastructure assessment. As shown in Figure 14, generally, the structure features a single center (Harbin city), a core distribution, a point-axis, and main cities and towns distributed in series along a trunk network that is made up of highways and railways. Measuring the extent of the single center by the urban primacy index showed that the centrality of Harbin city is predicted to keep increasing from 25.74 (2015) to 28 (2035). The urban-rural hierarchical size structure can be described as a seven-level settlement system composed of a central city, a sub-central city, a central town, a major town, a common town, a central village, and a general village. The urban-rural traffic network structure comprises railways and highways of different grades and forms a pattern of "radiation + grid", with Harbin city as the center of the structure. The urban-rural land use zoning structure was classified into eco-spaces, agricultural spaces, and construction spaces.

Integration of Urban-Rural Spatial Structures
The obtained urban-rural hierarchical size structure, traffic network structure, and land use zoning structure were integrated to obtain the overall urban-rural spatial structure of Harbin in 2035 based on the green infrastructure assessment. As shown in Figure 14, generally, the structure features a single center (Harbin city), a core distribution, a point-axis, and main cities and towns distributed in series along a trunk network that is made up of highways and railways. Measuring the extent of the single center by the urban primacy index showed that the centrality of Harbin city is predicted to keep increasing from 25.74 (2015) to 28 (2035). The urban-rural hierarchical size structure can be described as a seven-level settlement system composed of a central city, a sub-central city, a central town, a major town, a common town, a central village, and a general village. The urban-rural traffic network structure comprises railways and highways of different grades and forms a pattern of "radiation + grid", with Harbin city as the center of the structure. The urban-rural land use zoning structure was classified into eco-spaces, agricultural spaces, and construction spaces.

Integration of Urban-rural Spatial Structures
The obtained urban-rural hierarchical size structure, traffic network structure, and land use zoning structure were integrated to obtain the overall urban-rural spatial structure of Harbin in 2035 based on the green infrastructure assessment. As shown in Figure 14, generally, the structure features a single center (Harbin city), a core distribution, a point-axis, and main cities and towns distributed in series along a trunk network that is made up of highways and railways. Measuring the extent of the single center by the urban primacy index showed that the centrality of Harbin city is predicted to keep increasing from 25.74 (2015) to 28 (2035). The urban-rural hierarchical size structure can be described as a seven-level settlement system composed of a central city, a sub-central city, a central town, a major town, a common town, a central village, and a general village. The urban-rural traffic network structure comprises railways and highways of different grades and forms a pattern of "radiation + grid", with Harbin city as the center of the structure. The urban-rural land use zoning structure was classified into eco-spaces, agricultural spaces, and construction spaces.

Discussion
This study simulated and predicted the overall arrangement of the urban-rural spatial structure under the basis of ecological constraints in order to maintain the lowest effect on ecosystem security when urban-rural development is in progress. Controlling urban growth is a common approach for promoting urban-rural sustainable development [73,74]. Green infrastructures in this study were assessed in terms of resource protection and ecological services, instead of only from the perspective of ecological services [51,75]. The results showed that key areas, like natural reserves, gained a higher value of importance. This would be especially helpful for the establishment of subsequent development strategies.
This study considered settlements, traffic systems, and land use zoning as the three determinants of urban-rural spatial structures. In the process of hierarchical size structure extraction, although the Moran I index, which can reflect the extent of spatial concentration or dispersion, is lower for 2035 than for 2015, medium and small cities of Harbin will be further developed by 2035. These results are in accord with the national strategy in China. The rapid development of medium and small cities could have positive effects on the development of surrounding rural areas [76], which is of great significance for realizing balanced regional development and urban-rural integration, especially for developing countries. Meanwhile, the urbanization progress of the central city (Harbin city) is also continuing, whose speed is faster than that of medium and small cities. Consequently, the urban primacy index increased, which indicates that the centrality of the central city will still increase. The traffic network supports urban-rural material flow and is the necessary basis for urban-rural integration. Many studies have predicted traffic flow and explained spatial-temporal correlations [77,78], or forecasted traffic volume through a significant number of surveys [79], but their aims were not to explore the spatial layout of the traffic network. According to the simulation result of this study, the traffic network structure of Harbin in 2035 tends to be more complex. The grade and quality of the roads will also be enhanced. Drawing on the experience of urban-rural integration in developed countries [80,81], the current service capacity including the grade and quality of highways and railways will be improved to connect all settlements through the traffic network. It should be noted that because the Harbin railway development will mainly be fulfilled by upgrading existing railway facilities, the layout of the railway in this study was maintained the same as the current situation. This could be further studied when more data are available. The land in this study was zoned in terms of living, productive, and ecological functions responding to national spatial planning in China. Finally, the overlaid result showed that the urban-rural spatial structure still featured a significant single center. This is currently a common phenomenon for the development of regions in China. A review of urbanization in developed countries has found that the tendency to develop a single center could be an inevitable process preceding urban-rural integration [82].
In accordance with the equivalents of the ecosystem service values supplied by the per unit area of the ecosystem [50], water bodies can offer a superior water resources supply service, hydrological adjustment service, biodiversity support service, and cultural service; forests can provide a better gas regulation service, climate regulation service, and biodiversity support service; while agricultural land can offer a better food supply service. The application of the integrated model in this study can lead to policy recommendations that may support more sustainable development in Harbin. For example, through the policy of prohibiting development activities in water bodies and forest areas and the restricting misuse of agricultural land, multiple ecological services including water resource supply services, hydrological adjustment services, biodiversity support services, and so on, can be preserved to support sustainable development in Harbin.
Previous research on spatial structure only simulated land use structures from the perspective of the surface [83][84][85], while the integrated model developed in this study comprehensively predicts the urban-rural spatial structure in the future from a point-line-surface perspective. In addition, the model of this study can determine the sequence of land development in advance on the basis of the results of a green infrastructure assessment. Then, the urban-rural spatial structure can be simulated. Therefore, the integrated model in this study can be considered as a model of sustainable development. The simulation results based on this model could effectively protect important ecological resources. To further explain the strength of the integrated model in this study, a comparison between the ecological constraint scenario that we developed and the scenario of urbanization as usual should be undertaken. This work is being done at present. Furthermore, to realize the results of the simulation in this study, some optimization measurements must be undertaken, for example, how to promote the fast development of medium and small settlements in Harbin, how to improve the service capacity of the traffic network in Harbin, and how to control and manage developing activities in ecological, agricultural, and construction spaces for Harbin. All of these questions may be addressed in future research in this area.
This study focused on problems at the macro scale and considers urban areas and rural areas in a comprehensive way. To some extent, by giving priority to the protection of areas with higher ecological importance to restrict the location of construction land development, the problem of mixed urban-rural land uses (which entails more detailed and complicated issues such as fragmentation, leap-frog development, and low-density urbanization) can be modified. However, this study is not concerned about problems of large-scale development within the city. For these problems, some restoration and compulsory measures will be necessary. This will also be studied more deeply in the future.

Conclusions
From an ecological perspective, this study proposed an integrated simulation model to simulate the urban-rural spatial structure of Harbin for 2035, the result of which is expected to further promote the sustainable development of urban-rural integration through the limitation of green infrastructure on land use transitions induced by urbanization. This study, using a quantitative evaluation method, the logistic-CA model, the LCP model, and a classification of ecological importance, simulated the hierarchical size structure, traffic network structure, and land use zoning structure, and integrated them into an urban-rural spatial structure for Harbin. The conclusions based on the simulated results from this integrated model are as follows: • The integrated simulation model based on green infrastructure assessment in this study can be used to simulate urban-rural spatial structures.

•
According to the simulated results of this study, although land use patterns in Harbin in 2035 will be more dispersed, the medium and small cities of Harbin will become developed.

•
The simulated results indicate that in 2035, the traffic network of Harbin will still consist of a railway network and highway network, and the traffic systems above the county level will be improved.

•
For Harbin, based on the simulated results, most of the land in 2035 will be covered by ecological land, followed by agricultural land, while construction spaces will be primarily distributed in the smoother western area.

•
Through the integrated simulation model, the urban-rural spatial structure of Harbin shows an increasing development tendency for a single center. The main cities and towns are distributed in series along with a trunk network.