A Spatial Zoning Model of Municipal Administrative Areas Based on Major Function-Oriented Zones

: The multi planning contradictions and maladjustment of spatial management that evolved from an unbalanced development are becoming increasingly serious in China; therefore, it is imperative to know how to coordinate spatial planning. The scientific division of spatial unit attributes and the identification of ecological – agricultural production – living spaces (EPLs) have become the key scientific problems of this period. The existing urban structure zoning model, the urban growth boundary model, and the function-oriented zoning model show three main shortcomings: the difficult data acquisition and the large limitation of the application because of their high requirements of the related parameters, strong subjective weight setting, and lack of comprehensive zoning in the EPLs practice. This paper builds a quantitative and easily operated model for the spatial zoning of municipal administrative areas (MAA) based on the existing spatial equilibrium model for regional development. Three representative cities in China were modeled empirically with this model and the results were compared with those obtained by the existing models. We investigated the image consistency of this model to the three existing models where the Kappa values were 85.9%, 88.2%, and 85.2%, respectively, with an average of 86.4%. This showed that the model could reduce the data limitation and expand the scope of the application while ensuring the accuracy of the model’s analysis results. Meanwhile, the EPLs zoning was clearer, which made the spatial plan coordination more efficient and scientific. Overall, this model could not only solve the problem of the EPLs zoning and spatial plan coordination in MAA, but also guide urban land use planning from two dimensions of space and time and effectively promote the coordination and sustainable development of spatial planning.


Introduction
Spatial zoning is a public policy to achieve strategic guidance, rigid management, and control.The question of how to coordinate spatial zoning has become a global topic [1].Meanwhile, the content and scale of spatial zoning are distinctly different from the national conditions and development.Overall, the spatial zoning of developed countries such as in Europe, the United States and Japan, is relatively mature and complete [2][3][4], while that of developing countries such as China is relatively backward and unsystematic.In terms of the content of spatial zoning, developed countries have turned from material construction planning to social development planning, and the social and environmental factors involved have been given much more attention [5].From the scope of spatial zoning, more attention has been paid to wide range spatial zoning.For example, the regional development planning of Holland, Britain, Germany, and Hungary take the whole country as the planning object [2,[6][7][8].While undertaking a comprehensive survey of the spatial zoning of developed countries, we found that they had three important characteristics: (1) This spatial zoning was relatively complete [2][3][4] and the level was clear [1].Taking Japan as an example, the national comprehensive development plan has been compiled since the 1960s and a round of planning is compiled every ten years.At present, it has established a perfect planning and management system and formed spatial zoning for the full coverage of the country to wide areas, and then to city-townvillage, which was adapted to the Japanese administrative system [9].(2) Spatial zoning has changed from a traditional concept with economic growth as the primary goal to the comprehensive and sustainable development of the economic, social, and ecological environments [10][11][12][13][14][15][16][17].(3) Much attention has been paid to the major existing problems of these countries such as the industrial structure adjustment after economic growth, resource and environment pressures, the imbalance of regional and intra-regional development, energy shortages, global climate change, and so on [18][19][20][21][22][23][24].In all, the spatial zoning of developed countries has made positive contributions to the local development in several ways such as a balance of regional development, control of land use, environmental protection, and so on [4].In comparison, China is a developing country with a high population density and the land ownership of China is land nationalization.Meanwhile, Chinese urbanization is developing at a high speed [25,26].All of these natures make Chinese spatial zoning unique and very different from that of developed countries.The study of Chinese spatial zoning can not only be used to supplement and refine the existing spatial zoning of western developed countries, but can also provide strong support and a basis for the local spatial planning (zoning) of all levels of governments in China.
Municipal administrative areas (MAAs) are a key component of administrative zoning systems and an important field of zoning adjustment in China.Spatial zoning based on MAAs plays an important role in administrative systems and promotes the implementation of space control policies.However, in actual urban spatial zoning, plans are often contradictory and the effect of management overlap is particularly serious.Meanwhile, the planning system is hierarchical, and strategies where the government value the city but neglect the rural areas, and extensive land use have caused great difficulties for governments at all levels [27].To overcome these problems, the construction of a unified spatial system to identify ecological-agricultural production-living spaces (EPLs) [28,29] is crucial for spatial plan coordination [30].Here, the EPLs include the ecological space mainly used for ecological services and ecosystem maintenance (such as woodlands, grasslands, and so on), the agricultural production space that is mainly for agricultural production (such as basic farmlands, general farmlands for agricultural production, etc.), and the living space, which mainly undertakes the construction and economic development of the urban and rural areas (such as the urban built-up area, urban planning area and development zones, and rural living land).
As for the empirical research into spatial zoning reported previously, most do not use a pattern of equilibrium development [31], which gives priority to the development of MAAs.Objectively speaking, this pattern helps society achieve a higher economic growth [32,33].However, there are growing gaps of economic, societal, and environmental factors between urban and rural MAAs, resulting in problems such as the deterioration of urban environmental quality, house and traffic crowding, and an increase in urban crime rates.Meanwhile, the amount of agricultural and ecological land has decreased, and the ecological environment has been damaged [34][35][36][37].Recently, with the proposal of a more scientific outlook on development, decreasing the gaps between urban and rural areas and ensuring the equilibrium development of the economy, society, and the environment in MAAs under limited resources have become the main principles of urban spatial zoning [38].Similarly, developed countries such as Germany, France, the United States, and the United Kingdom have encountered the same difficulties in the process of industrial development.They have also adopted strategies of equilibrium development to balance the development of land space and environmental protection [39][40][41][42].
There has been a great deal of research on MAA spatial zoning models, and three main models have been established: the urban structure zoning model, the urban growth boundary model, and the function-oriented zoning model.Of them, the urban structure zoning model has been widely used for early MAA spatial zoning.Using this model, the mono-centric model and multi-centric model have been used to study the influence of city sprawl in Seoul and the Mediterranean during different periods [43][44][45][46][47][48].For the urban structure zoning model, its mathematical expression has advantages in exploring the influence of spatial dynamic change factors and their relationships.With the intensification of urban expansion to agriculture and ecology, empirical research on urban growth boundaries have been studied using many methods such as GIS and cellular automaton, and the urban growth boundary model became the mainstream model for the zoning of regional living space [49][50][51][52][53][54][55].Generally, the urban growth boundary model realizes the boundary simulation of urban living space in different periods using different zoning methods and establishes an evaluation system of subjective weight settings.Therefore, it is useful in guiding urban expansion and land use.With the increasing chaos of urban spatial functions, the function-oriented zoning model has become increasingly useful.This model is divided into two categories.The first is the Major Function-Oriented zoning model such as the spatial equilibrium model for regional development, which is the most representative model [56].This model is based on the main influencing factors in the economy, society, environment, and politics, and the model uses the method of setting a subjective weight to zone the major function of the region, providing the basis for regional function orientation and a future development direction.It has been widely used in land development [57][58][59][60][61] in many of the municipal-county administrative areas of China, and the correctness of the model has been verified empirically.The second is the EPLs zoning model.In recent years, with the proposal of multi-plan coordination, growing attention has been paid to a number of EPLs zoning models using subjective weight settings based on function zones.For example, Liu and colleagues explored an evaluation method to classify and evaluate EPLs based on the classifications of land use [62].This evaluation system is a reference for the classification planning and optimization of EPLs.The function-oriented zoning model has aided in the construction of a set of subjective evaluation systems, and due to its practicality, it has been widely applied in the field of spatial zoning.
Generally, the urban structure zoning model studies the factors that affect the urban form and the relationship between them based on strong economic logic derivation.It has been widely used in urban spatial zoning and urban sprawl.However, the factors it studies belong to the economic field and its combination with urban planning is not strong, and many factors cannot be quantified, so it cannot effectively solve the practical problems of spatial plan coordination.The key purpose of the urban growth boundary model is to construct an evaluation system of subjective weight settings and to simulate the urban development boundary during different periods using different models.Given the existence of subjective weight settings, the zoning of the spatial boundary obtained with this model is relatively subjective.As a result, it cannot be used effectively to guide the spatial plan coordination.The function-oriented zoning model can be used to effectively classify optimal development zones, key development zones, prohibited zones, limited development zones, and EPLs in a city, and gives a qualitative evaluation of urban functional attributes from the view of development.Nevertheless, in this model, there is a strong subjective influence of subjective weight settings on zoning [62].
In all, the existing MAA spatial zoning cannot effectively guide the implementation of MAA planning, which leads to various intensifying contradictions between the planning content and spatial governance.Meanwhile, the land-use planning indicator of MAA cannot be used directly, construction projects are difficult to perform, and the management of ecological space is disordered.These problems have led to the resistance of social and economic development in MAAs.In order to solve these problems, it is of utmost importance to establish a scientific and effective MAA spatial zoning model (M-MSZ) for the zoning of MAAs.In this study, based on the spatial equilibrium model for regional development, we constructed an M-MSZ under the framework of spatial equilibrium theory [63][64][65].The obtained model was optimized by a statistical analysis and compared with the existing spatial zoning model (including the urban structure zoning model, urban growth boundary model, and major function-oriented zoning model) under the same factors.It was found that the model showed advantages related to its quantitative, objective, scientific mechanism, and easily accessible data sources.Moreover, it provided an efficient and unified EPLs zoning map for MAA spatial plan coordination, and it showed a higher cost-effectiveness, and was more efficient and precise in the simulation of urban land use planning.
It is necessary to explain that in this paper, the major function-oriented zones' delimitation of MAA was used to coordinate with the M-MSZ to formulate a zoning standard of MAA.EPLs zoning was carried out based on this zoning standard, as shown in Figure 1.Meanwhile, the M-MSZ realized the EPLs spatial zoning of MAAs by constructing an evaluation system of spatial development suitability.Therefore, the M-MSZ can not only carry out the delimitation of EPLs, but can also realize the zoning of spatial development suitability.In the following part of this paper, the applications in both EPLs' delimitation and spatial development suitability of MAA were conducted.

Construction, Analysis, and Optimization of M-MSZ
For the zoning of MAAs, we first need to consider the spatial development suitability.Starting with the selection of factors that affect the development of the MAA space, we constructed the initial model through the mathematical modeling of economics, then used a correlation analysis and regression analysis to determine the effective variables.After this, regression analysis was used to illustrate the explaining the degree of variables to dependent variables.

Selection of Factors Influencing the Suitability for Spatial Development
How to select important representative factors from a large number of factors that influence the spatial function of MAA was the focus of building the M-MSZ.In light of the existing problems of MAA spatial zoning in China such as intensified planning contradictions, the difficult implementation of construction projects, disordered ecological space management and control, etc., we selected the influencing factors of the M-MSZ through statistics found in the literature, expert advice, and a questionnaire survey of MAA participants regarding the leading factors that affect Chinese MAA spatial zoning.The questionnaire surveys are shown in Section 1 of the Supporting Information.
The literature selection was carried out by combining forward retrieval with reverse retrieval.For forward retrieval, we studied the literature published in the Web of Science, KCI-Korean Journal Database, Russian Science Citation Index, MEDLINE, Derwent Innovations Index, SciELO Citation Index, and CNKI from 1950-2018.(1) The topics or keywords included spatial equilibrium, spatial zoning, urban spatial growth boundary zoning, urban structure model, major function-oriented zones delimitation, planning coordination, multi-plan coordination, three-plan coordination, spatial planning, sustainable development model, land use, and land expansion.(2) These studies were published in journals on sociology, architectural science and engineering, environmental science and resource utilization, macroeconomic management, and sustainable development.(3) The types of literature included articles, meetings, reviews, patents, books, and letters.(4) The research domains included the social sciences and science technology.(5) The research directions included geography, urban studies, sociology, operational research, and management science.(6) The literature selection was chosen first in the order of greatest relevance, then the frequency of citation, and the publication date third.For the reverse retrieval, we investigated the cited references of each study that was found by the forward retrieval and found the relevant literature.Finally, a combination of the forward retrieval and reverse retrieval were used to eliminate any duplicated studies.A literature statistical survey was done from the statistics of relevant articles including 334 Chinese articles that were published during 1981-2018 (journals, conferences, and theses) and 424 SCI and SSCI papers that were published from 1961-2018 (journals and conferences).The operation of the literature surveys and a list of part of the literature are shown in Section 2 of the Supporting Information.We filtered the factors according to the occurrence frequency, so only the factors with a frequency of more than 60% were chosen for the M-MSZ, as shown in Table 1.Land renovation, administrative rank, and interregional interaction are not universal (factor occurrence frequency of less than 60%) and were not used as factors in the initial model.For expert advice, each group of subjects contained five experts, and the total number of experts was 50.The average score of each class of five experts was scored according to the 50 valid questionnaires.The opinion of each of them was equally important.The results are shown in Table 2 and include experts working in Chinese universities and experts in engineering planning and design in China.The experts working in universities majored in five fields including land planning, urban and regional planning, urban planning, urban design, and environmental protection and planning.The experts in engineering planning and design also majored in these five fields.It was found that land renovation, administrative rank, and interregional interaction received the lowest scores in the expert survey, which meant that they were less important than the other factors.Herein, only factors with a score of more than 60% were chosen for the model.Meanwhile, from the perspective of human-land relations, the scale of MAA living space also depends on the choice of MAA participants for the living space.A total of 25,341 valid questionnaires investigating the living space choices were made from 25 MAAs with different scales, as shown in Table 3.The 25 MAAs included Beijing, Shanghai, Harbin, Changchun, Taiyuan, Shijiazhuang, Jinan, Kunming, Lanzhou, Urumqi, Langfang, Qigihar, Jilin, Yichang, Lijiang, Baotou, Yinchuan, Sansha, Siping, Jingdezhen, Zunyi, Songyuan, Guangan, Huanan, and Gongzhuling.The questionnaire used the influencing factors of the related research as the main reference, and it was designed in the form of open answers with a structure of a single-factor scoring protocol.These factors with a score of more than 60 were considered to be the most important by the surveyed MAA participants when choosing a living space.As a result, 11 factors were chosen as the initial influencing factors of the MAA spatial development model, these being: economic situation, location advantage, transport superiority degree, population quantity, topography and terrain, land quantity, environmental capacity, ecosystem vulnerability, utilization of water resources, natural disasters, and urban structure.

Construction of M-MSZ
Through the selection of the above factors, we attempt to build an M-MSZ based on the spatial equilibrium model for regional development.Generally, it is believed that the realization of regional spatial equilibrium requires a reasonable flow of all kinds of resources (economic, social, and ecological) among the regions [56].Therefore, it is necessary to derive the inevitable relationship between the indexes of the MAA comprehensive development state to create a reasonable flow.The index reflects the utility of MAA spatial development.Here, we used D to represent the spatial development suitability of the MAA.Due to the existence of spatial equilibrium, the difference in D values for different regions can be generated.If we set the population quantity of region Ri and region Rj as Ni and Nj, respectively, then the spatial equilibrium model for regional development is given by Equation (1):   [56] (1) We assume that the response of MAA participants to land development is the same and the capital elements have complete liquidity.The transportation costs for each commodity from different enterprises to each MAA participant should be paid.The spatial distribution of the main bodies in the MAA can be used to decide the MAA spatial function [66].The construction process of the M-MSZ can be formed as follows: ( ) First, we simply considered the economic factors and divided the MAA commodities into two groups: land commodities and all other commodities.Then, the maximum utility of the MAA development was given as the price function of three variables including the economic situation of the MAA, and the price of land and other commodities [67] as shown in Equation ( 2), where D is the utility of spatial development suitability in MAA, Y is the economic situation, G is the price of other commodities, R is the price of land commodities, and μ is a constant indicating the share of the expenditure of other commodities on the total consumption of MAA commodities.
Given the existence of iceberg transport costs and the Marshallian demand for other commodities during transactions [67][68][69], the price of other commodities is a function of three variables.These three variables encompass the preference of MAA participants for other commodities [70,71], their own location, and the number of other commodities available [72], as shown in Equation ( 3), where t is the transaction cost of other commodities per unit distance; m is the type of the other commodity; σ represents the substitution elasticity between two kinds of MAA commodities; P1 is the preference of MAA participants to other commodities; y + r0 is the location advantage of MAA participants; y represents the external location advantage (i.e., the distance between the research area and the urban center of the surrounding cities); and r0 is the internal location advantage (i.e., the distance between the research area and the MAA center).
In this paper, we assumed that the size of the urban space was determined by the number of urban production enterprises [73], which is related to the urban structure.Given the existence of iceberg transport costs, we can obtain Equation (4).M(r) is an urban structure with a radius of r (the distance from the MAA center to the research boundary).Then, Equations ( 3) and ( 4) can be entered into Equation ( 2) to obtain Equation (5).
As there is a Marshallian demand during the transaction of land commodities [68], we used λ = (1-μ)/μ and obtained Equation ( 6), where L is the total land commodities.As is already known, there is an inverse relationship between L and population agglomeration [74].When the function of population agglomeration is integrated, we can obtain the population quantity of the area, N, as shown in Equation (7).Then, the population quantity is a function of N1 (population agglomeration) and L, as shown in Equation ( 8).Generally, it is believed that there are two factors affecting the price of land commodities in MAAs: the location advantage (y + r0) and transport superiority degree (J) [75,76].However, from the aspect of spatial equilibrium, MAA participants also consider the factors of resource difference, environment, and externality when choosing a living space.Therefore, the preference variable of MAA participants to land commodities can be introduced as P. According to questionnaires regarding the main consideration of MAA participants from 25 MAAs with different scales when choosing the living space, P is predominantly attributed to five aspects: topographic and terrain, natural disaster, utilization of water resources, environmental capacity, and ecosystem vulnerability.In this study, the development utility function focused on the comparison of MAA land development suitability instead of the utility calculation of other commodities, so the preference of MAA participants to other commodities was not considered.If we assumed that MAA participants had the same preference for other commodities, then p1 = (σ−1)/σ and Equation ( 9) can be obtained, where X1 is the share of preference of MAA participants to the price of land commodities, X2 is the share of location advantage in the price of land commodities, and X3 is the share of transport superiority degree in the price of land commodities.When Equations ( 8) and ( 9) are entered into Equation ( 5), we can obtain Equation (10), which is the initial form of the M-MSZ built in this paper.

Statistical Analysis and Optimization of M-MSZ
The model was analyzed and optimized by using a statistical method, and the final model was determined on the basis of optimization.In this paper, 420,436 random valid samples from the three MAAs above-mentioned in Heilongjiang Province were selected.The land sample was taken from the part of the three MAAs and was merged and decomposed according to the above 13 factors, and then 420,436 vector landscapes with different sizes and factors were obtained.Correlation and regression analyses were carried out.The correlation analysis was based on the nonlinear correlation, and the results of curve correlation coefficients and significance are shown in Table 4.This showed that the Sig value of all of the coefficients was compared with 0.05 and 0.01, and it can be seen that only the difference between the coefficients of x6 and x13 were not significant.The regression analysis was done according to the multivariable nonlinear regression, and the regression analysis parameters and the significance are shown in Table 4.By comparing F with Fa (k-1, n-k) of various parameters, it was found that only x6 and x13 showed no significant influence on the interpretation of the dependent variable, y.All of the other parameters had a significant influence on y, which was similar to the results of correlation analyses.Based on the above analysis, this paper eliminated x6 and x13 because of their poor relevance and explanation to y.After the elimination, the correlation between the other factors was not strong, as shown in Table 5, which proved that each of the selected factors could be regarded as an independent factor.Therefore, we removed x6 and x13.Then, the value of L was mainly determined by the available land volume, L1, and the value of M was mainly determined by M1, or the pending development zones.Through the above analysis, the model was adjusted to give Equation (11). (

Numerical Simulation Analysis of M-MSZ
Based on the above models and Equation (10), the numerical simulation was carried out under the support of MATLAB to obtain the contribution rate of various parameters in the model.As we know, the values of the coefficients such as μ, t, x1, x2, x3, and σ will not affect the evaluation of spatial development suitability in the same MAA.Therefore, this paper used the coefficient values reported by Fujita et al., which were taken as 0.5, 0.005, 0.45, 0.2, 0.35, and 2, respectively [77][78][79].
These coefficient values have also been cited by many other works [80][81][82].Meanwhile, it was assumed that MAA participants had the same aspirations in the choice of living space.The simulation results are shown in Figure 2. It can be seen that the attributes of land quantity, population quantity, and economic situation in the model are higher than any other aspects.Therefore, this suggests that developers will not pay much attention to the constraints of other factors and develop the region as much as possible when the land quantity is sufficient, the population size is big, and the economic situation is favorable.When the regional land quantity, population quantity, and economic situation are unfavorable, developers will consider all of the factors that influence MAA development and make a comprehensive judgment of whether the researched zone is suitable for development.Compared with land quantity, population quantity, and economic situation, the contribution rates of location advantage, transport superiority degree, and urban structure to the model were low.Meanwhile, the five parameters of P including topographic and terrain, natural disaster, utilization of water resources, environmental capacity, and ecosystem vulnerability showed the lowest attribution to the model.Through the numerical simulation analysis, it was found that the contribution rate of all kinds of influence factors to the M-MSZ conformed to the basic situation of current MAAs zoning in China [57][58][59][60][61][62].

Validation and Application of the M-MSZ
This model cannot only zone EPLs in MAA, but can also make the MAA spatial zoning more reasonable due to the introduction of an evaluation system of spatial development suitability.In addition to the simulation of the application in EPLs zoning, the comparison of M-MSZ and three other existing models in the application of spatial development suitability were also made to verify the rationality and superiority of the M-MSZ.

Introduction to the Research Area
The Heilongjiang Province is located in the northeast of China, and its development has been relatively slow when compared with the developed provinces.The spatial zoning method of the Heilongjiang Province is relatively subjective and traditional, but the new urbanization of China requires the rapid, steady, and sustainable development of MAA.The mismatch between the spatial zoning method and the urban development goal is a common problem in most developing countries.To reduce the negative influence of the subjective factors and the constraint of urban zoning on urban development, it is necessary to decrease the contribution rate of the subjective factors in the evaluation method of spatial zoning.This paper selected three representatives of MAAs in the Heilongjiang Province according to the type of major function, these were MAAs with the major function of development, agricultural production, and ecological protection, respectively, as shown in Figure 3.The MAAs with the major function of development (MAA-D) were located in the southwest of the Heilongjiang Province.In 2015, the urban population accounted for half and the urban area accounted for about 1.2% of the total area of the MAA, and the GDP of the second and third industries accounted for 85% of the total GDP.The MAA with the major function of agricultural production (MAA-AP) was located in the northeast of the Heilongjiang Province.In 2017, the agricultural population accounted for half of the total area of the MAA, and the GDP of the first industry accounted for about half of the total GDP.Cultivated land accounted for more than 60% of the total MAA and the forest coverage rate was about 20%.The MAA with the major function of ecological protection (MAA-EP) was located in the southeast of the Heilongjiang Province.In 2015, the agricultural population accounted for about half of the total area of the MAA, the GDP of the first industry accounted for a quarter of the total GDP, and the area of cultivated land accounted for about 13% of the total MAA area and the forest coverage rate reached 70%.The basic information of the research areas is shown in Table 6.

Validation of the M-MSZ on Spatial Development Suitability
The representative spatial zoning model can be divided into three categories: the urban structure zoning model, urban growth boundary model, and function-oriented zoning model, respectively [43,49,61].Here, we took a city in northeastern Heilongjiang as an example and compared the M-MSZ and these three existing models in space development suitability to verify the rationality and superiority of the M-MSZ.
First, the evaluation of space development suitability was done by using the M-MSZ and the urban structure zoning model, respectively.Different from the urban growth boundary model and the function-oriented zoning model, both the M-MSZ and the urban structure zoning model can be used to grade the spatial development suitability of living space.In order to verify the grading function of the M-MSZ for spatial development suitability, this paper compared the application results of the M-MSZ with the urban structure zoning model in the grading of spatial development suitability.Here, the divided living spaces included two levels: the first level is the most suitable living space, and the second level is less suitable living space.The results are shown in Figure 4.It can be seen when the living space is concerned, that the consistency of the Kappa value of the M-MSZ to the urban structure zoning model is great.The Kappa value of the first level and the second level of living space were 86.4% and 85.3%, respectively, and the average Kappa value was 85.9%.We believe that the discrepancy was mainly caused by the difficult factors assigned in the urban structure zoning model.Furthermore, we note that the application of the urban structure zoning model was not as wide as the M-MSZ because it can only be used to zone the living space when the latter can be used to zone the EPLs of MAA.

Model Application on EPLs Zoning
The development mode and cultural style are diversified in cities, which leads to differences in their comprehensive development level per capita (Di) in different regions.As noted earlier, the status of economic development (Di1), social development (Di2), and ecological environment (Di3) are closely related to the functional positioning of MAAs.That is, the MAA for economic development usually reaches a high value for the index of economic development, while the MAA for ecological protection usually reaches a high value for the index of the ecological environment [56].
In the study of MAA zoning, we need to determine the best functional positioning of the MAA to give the most favorable functional advantages of the MAA as well as to narrow the gap between the Di.
Therefore, this paper introduced the major function of MAA into the M-MSZ.The index system of major function orientation proposed by Hasbagen [83] was selected to determine the F value (comprehensive classification index), then the major function of the MAA was limited.This evaluation system has been widely applied in previous MAA studies [84][85][86][87], and it has been shown to effectively guide the major function orientation of MAAs.The MAA was divided into three groups according to the major functions: the MAA space mainly for development, the MAA space mainly for agricultural production, and the MAA space mainly for ecological protection.Then, the table of MAAs zoning standards and the classification table of urban built-up areas, transitional zones, and limited development zones are determined according to the actual situation of the research zone and how the MAAs space is zoned.

b) a)
From the perspective of spatial equilibrium, we should follow the principle of giving priority to the major function of the MAA and take into account the principles of protection and utilization.We divided the MAAs space according to the 2017 edition of the "Classification of City Land Use and Standard of Land Use Planning and Construction Land" [88], the 2014 edition of "A Guide to the Classification of Village Planning Land, Housing Construction Department" [89], the "Present Land Use Classification" (GBT 21010-2017) [90], the 2017 edition of "Technical Guide of Red Line for National Ecological Protection and Ecological Function Delineating" [91], and the relevant standards of forestry, agriculture, transportation and other departments.The standards are shown in Tables 7  and 8.

MAA-D Living Space
The urban built-up area; the type I transitional zone that is adjacent to the urban built-up area of level one, two or three; sand barrier, stacking material, other artificial dig land, saline-alkali land, clay land, gravel land, rock land in the type II transitional zone that is adjacent to the urban built-up area of level one or two; type III transitional zone of level one, two or three.Open pits with large destruction in type III transition areas of level one, two or three.Other zones in type III transition areas.

Agricultural Production Space
Protective zone of basic farmland; the type I transitional zone that is adjacent to the urban built-up area of level four or five; the type I transitional zone that is not adjacent to the urban built-up area.Other zones in type III transitional zones.

Ecological Space
Limited development zones except for protection areas of basic farmland; the land in type II transition zones used for producing natural ecological forest and grass, as well the upland and paddy field with slope of bigger than 25 o except the basic farmland; sand barrier, other artificial dig land, saline-alkali land, clay land, gravel land, rock land in type II transitional zones adjacent to urban built-up area of level three, four or five.Type II transitional zone that is not adjacent to the urban built-up area; Open pits with large destruction in type III transition areas of level four or five.

Living Space
The urban built-up area; the type I transitional zone that is adjacent to the urban built-up area of level one; sand barrier, stacking material, other artificial dig land, saline-alkali land, clay land, gravel land, rock land in the type II transitional zone that is adjacent to the urban built-up area of level one or two.Other zones in type III transitional zone.

Agricultural Production Space
Protective zone of basic farmland; the type I transitional zone that is adjacent to the urban built-up area of level two, three, four or five; Paddy fields and uplands of the type II transitional zone except basic farmland; the type I transitional zone that is not adjacent to the urban built-up area; the type III transitional zone with level one, two or three.Open pits with large destruction in type III transition areas

Ecological Space
Limited development zones except for protection areas of basic farmland; the land used for producing natural ecological forest and grass in type II transition zone; sand barriers, other artificial dig land, saline-alkali land, clay land, gravel land, rock land in type II transitional zone adjacent to urban built-up area of level three, four or five.Type II transitional zone that is not adjacent to the urban built-up area; Open pits with large destruction in type III transition areas.

Living Space
The urban built-up area; the type I transitional zone that is adjacent to the urban built-up area of level one; sand barriers, stacking material, other artificial dig land, saline-alkali land, clay land, gravel land, rock land in the type II transitional zone that is adjacent to the urban built-up area of level one or two.

Agricultural Production Space
Protective zone of basic farmland; the type I transitional zone that is adjacent to the urban built-up area of level four or five; the type I transitional zone that is not adjacent to the urban built-up area.

Ecological Space
Limited development zones except for protection areas of basic farmland; the type I transitional zone that is adjacent to the urban built-up area of level two or three; the land in type II transition zone used for producing natural ecological forest and grass, as well the upland and paddy field with slope of bigger than 25 o except the basic farmland; sand barrier, other artificial dig land, saline-alkali land, clay land, gravel land, rock land in type II transitional zone Type I The arable land, garden land, and other agricultural production lands that with a slope smaller than 25 o except the protective zone of basic farmland.

Type II
A woodland and grassland dominated by producing natural ecological forest and grass.The upland and paddy fields with slope bigger than 25 o except the basic farmland.Sand barriers, other artificial dig lands, saline-alkali land, clay land, gravel land, rock land.Type III Open pits with large destruction and other areas.

Other protective zones
Natural reverse Scenic area Forest parks Geo-parks World cultural and natural heritage Land for water and water conservancy facilities Wetland Protection zone for drinking water source Zone that is not suitable for the development of large-scale industrialization and urbanization Using the M-MSZ, the development suitability of the urban built-up area, the control zone (the basic farmland, protected zones, and prohibited and limited development zones), and the transitional zone (except for the urban built-up area and the control zone) of the MAAs located in Heilongjiang Province were calculated.Then, the MAA was divided into levels one to five according to the calculation results.The first level was mostly for development; whereas, basic farmland, protected, and prohibited development zones were level five.
If F > 0, it is suitable for development.If F < 0, it is not suitable for development.There are areas where F < 0 can be divided into two parts: agricultural zones and ecological zones.When the proportion of farmland is greater than that of ecological land, it is defined as an agricultural area.
Otherwise, it is an ecological region [83].When F = 0, this is the threshold for developing protective function zones.According to the actual situation in the Heilongjiang Province, this paper defined it as a zone suitable for development.
According to the actual situation of various MAAs in Heilongjiang, 11,546 valid questionnaires were obtained from experts and MAA participants.For the expert advice, each group of subjects contained three experts, and the total number of experts was 30.The average score of each class of five experts was scored according to the 30 valid questionnaires.The opinion of each expert was equally important and included experts working in Chinese universities and experts in engineering planning and design in China.The experts working in universities majored in five fields including land planning, urban and regional planning, urban planning, urban design, and environmental protection and planning.The experts in engineering planning and design also majored in these five fields.MAA participants included residents, the government, developers, and farmers.Experts and MAA participants were investigated with different consulting methods.The opinions of each type of research subject are equally important and the weight ratio of the opinions of experts to MAA participants was 7:3.Next, all of the opinions were combined and the assignment of a neural network was carried out.Then, the weights of the various factors in the P value were set by using the analysis method of the neural network from the decisions of the expert and MAA participants [93] as follows: natural disasters, 0.15, topography and terrain, 0.15, utilization of water resources, 0.3, environmental capacity, 0.2, and ecosystem vulnerability, 0.2.
By using the proposed M-MSZ, the development suitability of these three MAAs was leveled, and the major function of the MAA was determined using the index system of major function orientation.Then, a MAA spatial zoning map showing different major functions was generated using the development suitability level, the major function of the MAA, the zoning standard table of the MAA, and the classification table of urban built-up areas, transitional zones, and limited development zones as shown in Figures 7-9.In the MAA-D, the living space, agricultural production, and ecological protection accounted for 8.2%, 41.37%, and 50.43% of the total area of the MAA, respectively.Compared with the actual value, living space increased by 5.7%, while agricultural space and ecological space decreased by 4.8% and 0.9%, respectively.The estimated rate of the population growth was 5‰ per year and the construction land per capita was 120 m 2 (the following two research areas are based on this standard), and the growth of living space was estimated to reach the spatial zoning boundary in 2070.As for the MAA-AP, the living space, agricultural production, and ecological protection accounted for 0.73%, 51.2%, and 48.07% of the total area of MAA, respectively.Compared to the actual value, living space and agricultural space increased by 4.8% and 0.9%, respectively.Furthermore, ecological space decreased by 13.07%.The growth rate of living space was estimated to reach the spatial zoning boundary in 2069.When living space reaches this boundary, agricultural production remains in surplus, except in the case of self-sufficiency.Where the MAA-EP is concerned, living space, agricultural production, and ecological protection accounted for 0.89%, 16.38%, and 82.73% of the total area of the MAA, respectively.Compared to the actual value, the living space and ecological space increased by 0.22% and 11.9%, respectively, and agricultural space decreased by 12.12%.As a result, the growth rate of living space was estimated to reach the spatial zoning boundary in 2072, as summarized in Table 9.

Summary
In summary, the M-MSZ could not only simulate the EPLs zoning, but also classify the EPLs from the space development suitability.Meanwhile, this model was highly consistent with the three existing spatial zoning models.The main reason for the discrepancy lay in the strong subjective weight setting and the difficult factor assignment. There, we can conclude that the M-MSZ was more efficient and cheaper than the existing spatial zoning models in EPLs zoning and spatial planning coordination.

Discussion and Conclusions
As we know, unbalanced development is the root of the increasing gap between urban and rural areas, which resists the urban and rural symbiosis in most cities and regions of the world [94,95].Scientific and rational MAA spatial zoning will help to solve many imbalance problems of MAA development such as urban and rural development and protection, land use, etc.Recently, the question of how to build a scientific and reasonable MAA spatial zoning model and coordinate spatial planning has received much attention.Currently, several spatial zoning models have been reported such as the urban structure zoning model, urban growth boundary model, and major functionoriented zoning model, and so on, and they have been successfully used to zone different functional zones and to guide the planning of MAA.Particularly for western developed countries, these models have been applied widely.The experience of them has provided strong support and reference for Chinese spatial zoning.However, it cannot be used directly for the spatial zoning of China for several reasons.First, each country chooses the model suitable for their development and makes the model work in their national conditions.Hence, different national conditions will inevitably create differences in the content and structure of spatial zoning, as shown in Figure 10.Here, we selected several representative developed countries that have shown advanced research in spatial planning to compare with China [2][3][4]6,7].As for China, this is a country with a unique national condition (public ownership of land, vast expanses of land, and a huge population) that has been undergoing the high-speed development of urbanization [96].Therefore, we cannot mechanically apply the spatial zoning models of western developed countries to China.Second, the assignment of factors of the existing models is often difficult and their subjective weight settings are strong, which makes them ineffective to guide the EPLs zoning and spatial plan coordination of China.In this paper, we launched theoretical and empirical research on the spatial zoning model based on Chinese national conditions and development stages.Existing spatial zoning models often show shortcomings such as difficult data acquisition and limited application due to their high requirements of the related parameters, insufficient objectivity because of their strong subjective weight settings, and lack of comprehensive zoning on EPLs.In this study, we attempted to construct a quantitative and easily operated MAA spatial zoning model (M-MSZ) by using the organic integration of economic logic deduction and urban planning practice under the framework of a spatial equilibrium theory, based on the planning experience on the balanced development of ecology, humanity, protection, and development.The M-MSZ and its construction method provide a new perspective and technological support to alleviate and even solve the extensive economic behaviors such as the inefficient operation of land resources (the speed of land urbanization is faster than the speed of population urbanization) and the deterioration of the ecological environment in China and other developing countries with rapid economic growth.It also demonstrates a new feasible demonstration for China to promote the implementation of new urbanization strategies in the new period and new normal situation [99].First, the selection of influence factors for the model was carried out from the literature statistics, the questionnaire survey of MAA participants (residents, government, developers, and farmers) about the leading factors of choosing a living space, and expert advice, which made the model factors more accurate.Second, during the construction of the model, the iceberg transport costs and Marshallian demand theory were introduced, and the development suitability and major functions of MAA were taken into synthetic consideration, which ensured the rationality and feasibility of the zoning model.Then, the model was optimized by statistical analysis and the rationality of the model was verified by the numerical simulation analysis, which guaranteed the scientific application of the model construction.Finally, typical MAAs were used as samples to complete the comparison and induction with the existing spatial zoning models, and the scientificity and ease of operation of the M-MSZ were verified again.The results of optimization, verification, and application are shown as follows.
During the construction of the model, the M-MSZ combined the advantages of the three existing spatial zoning models to reduce their disadvantages.For example, four factors that represent a high contribution rate and easy collection in the urban structure zoning model were retained, these being: economic situation, population agglomeration, location advantage, and transport superiority, respectively.At the same time, the total land commodities, urban structure, and the preference variable of MAA participants to land commodities were optimized, which made these three factors easier to acquire.As a result, the M-MSZ improved the practical application of the urban structure zoning model and its application was much easier to extend.On the other hand, the M-MSZ reduced the subjectivity of the urban growth boundary model and the major function-oriented zoning model (including the spatial equilibrium model for regional development and the existing EPLs zoning model) while maintaining their high practicality.The ratio of quantitative parameters was enhanced significantly by using the method of economic derivation.Many quantitative parameters such as the economic situation, population agglomeration, location advantage, transport superiority, available land volume, the pending development zones, and the preference variable of MAA participants to land commodities were combined in the model to make it more accurate and more objective in the mechanism expression.After comparison, it can be seen that the M-MSZ showed that it was more efficient, comprehensive, cost-effective, and widely used in the coordination of spatial planning and the EPLs zoning than the three existing spatial zoning models.Moreover, this model could also effectively guide the hierarchical scale zoning of the urban spatial system as well as the zoning of land use functions.
The correlation analysis showed that the relevance of unavailable land and the pending development zones to the suitability of urban spatial development was very small.From the comparison of the F values and the Fa (k-1, n-k) of each factor after regression analysis, all factors were F > Fa (k-1, n-k) except for the unavailable land and pending development zones, which was consistent with the results of the correlation analysis.After eliminating the unavailable land and pending development zones, the rest of the 11 factors including the economic situation, location advantage, transport superiority degree, population quantity, topographic and terrain, land quantity, environmental capacity, ecosystem vulnerability, utilization of water resources, natural disaster, and urban structure could be used to effectively explain the spatial development suitability.
The results of the numerical simulation showed that the economic development level, population quantity, and land resource were the main factors that affected the suitability of urban spatial development, while location advantage, transport superiority degree, and urban structure had a smaller influence on the suitability of urban spatial development.The five factors in the preference of MAA participants to land commodity including topographic and terrain, natural disaster, utilization of water resources, environmental capacity, and ecosystem vulnerability had the smallest influence on the suitability of the urban spatial development.Importantly, the accuracy of the M-MSZ was confirmed by the numerical simulation results.
The comparison of the M-MSZ and the three existing models (urban structure zoning model, urban growth boundary model, and major function-oriented zoning model) [43,49,61] showed that the Kappa values of the consistency test were 85.9%, 88.2%, and 85.2%, respectively, with an average of 86.4%.We believe that the M-MSZ efficiently inherited and extended the advantages and interpretations of the existing zoning modes.Moreover, they avoided the limitations of the application range of the existing models and looked like an organic integration of the three existing spatial zoning models.In addition, the EPLs zoning of the M-MSZ was clearer than the three existing models.Therefore, the M-MSZ could be used as a supplement and extension to the research of current existing spatial zoning.From the perspective of the application, this model can not only solve the problem of the EPLs zoning and spatial plan coordination in MAA, but can also effectively guide urban land use planning from two dimensions of space and time and effectively promote the coordination and sustainable development of spatial planning.When the application scope is concerned, the method and construction concept of M-MSZ are more practical to the spatial planning of regions including developing countries or underdeveloped areas (which occupy more than 70% of the world′s land area and the total population), and regions that are in rapid urbanization [100,101].
In conclusion, the M-MSZ can be used to solve the problems of urban and rural co-ordination and regional coordination such as delimiting EPLs of urban and rural areas in MAAs, determining land use functions of MAAs, evaluating the suitability of land exploitation in MAAs, and formulating the timetable of land development in MAAs.However, the M-MSZ is limited in its contribution to various transmission flows such as in the planning of city traffic, population, economy, and tourism [102].Although the M-MSZ is based on the mathematical derivation and comprehensive determination, there are still some subjective judgments such as the expert survey, the questionnaire survey of the factors affecting the choice of living space and its examination, and requires further verification of many more samples and regions.This is a potential area for our future research.

Figure 1 .
Figure 1.The scheme of research outline for Spatial Zoning of MAA basing on the M-MSZ and index system of major function orientation.

Figure 2 .
Figure 2. The curves of D versus r for different factors.(a) L1, the land quantity; (b) N1, the population agglomeration; (c) Y, economic situation; (d) y, the external location advantage; (e) r0, the internal location advantage; (f) J, transport superiority degree; (g) P, the preference variable of MAA participants to land commodities; (h) M1, urban structure.

Figure 3 .
Figure 3.The location of the research area.

Figure 4 .Figure 5 .Figure 6 .
Figure 4.The evaluation of spatial development suitability with (a) the M-MSZ and (b) the urban structure zoning model.

Figure 7 .
Figure 7.The spatial zoning of the MAA-D using M-MSZ.

Figure 8 .
Figure 8.The spatial zoning of the MAA-EP using M-MSZ.

Figure 9 .
Figure 9.The spatial zoning of the MAA-AP using M-MSZ.

Table 1 .
The occurrence frequency of each factor in the literature statistics.

Table 2 .
The expert scoring table on the factor selection (full score is 100).
a: economic situation; b: population quantity; c: land quantity; d: transport superiority degree; e: environmental capacity; f: ecosystem vulnerability; g: utilization of water resources; h: natural disaster; i: location advantage; j: topographic and terrain; k: interregional interaction; l: administrative rank; m: land renovation; n: urban structure; o: urban size.

Table 3 .
The questionnaires about the main consideration of MAA participants from 25 MAAs with different scales when they are choosing the living space.
a: economic situation; b: population quantity; c: land quantity; d: transport superiority degree; e: environmental capacity; f: ecosystem vulnerability; g: utilization of water resources; h: natural disaster; i: location advantage; j: topographic and terrain; k: urban structure; l: urban size.

Table 4 .
The correlation analysis with the R 2 , SIG, F, and Fa parameters.
Note: * is significantly correlated at the 0.05 level; # is significantly correlated at the 0.01 level; R Square = Residual Sum of Squares/Variance Sum of Squares; F = (Explained Sum of Squares/number of independent variable)/(Variance Sum of Squares/degree of freedom); degree of freedom = n-k-1; n: number of samples; k: the number of independent variables.

Table 5 .
The curve correlation coefficients and significance of the influence factors and variables.Note: * is significantly correlated at the 0.05 level; # is significantly correlated at the 0.01 level.y is the spatial development suitability; x1 is the economic situation; x2 is the location advantage; x3 is the transport superiority degree; x4 is the population agglomeration; x5 is the topographic and terrain; x6 is the unavailable land; x7 is the environmental capacity; x8 is the ecosystem vulnerability; x9 is the available land; x10 is the utilization of water resources; x11 is natural disaster; x12 is the number of the pending development zones; and x13 is the number of developed areas.

Table 6 .
The basic information of the research areas.

Table 9 .
The index of various space in MAAs after application with the M-MSZ.