The study area of this study is Hantai District which is located in the center of the Hanzhong Basin in the southwest of Shaanxi Province in China, with the Hanjiang River to the south and the Qinling Mountains to the north. The terrain of Hantai District is high in the north and low in the south, and the north belongs to the Qinling Mountains, with an elevation of 700–2000 m, accounting for 34% of the total area; the central part is a hilly area, with an elevation of 541–700 m, accounting for 28% of the total area; and the south is the Hanjiang alluvial plain, accounting for 38% of the total area (Wang L. et al., 2021). Hantai District is the largest commodity distribution center in southern Shaanxi and the core area of the Qinling and Bashan mountains, which has important economic and ecological value. Since the urban construction in Hantai District is mainly concentrated in the south, the study area is composed of the administrative boundaries of nine towns in the south (Figure 1).

This study utilized five types of data: GF-2 imagery, NPP-VIIRS imagery, POI data, WorldPop data, and administrative boundaries. The GF-2 satellite is the first civil optical RS satellite with a spatial resolution better than 1 m independently developed by China, which achieves sub-meter spatial resolution and multi-spectral comprehensive RS data acquisition (Zhou et al., 2021). It has the characteristics of high positioning accuracy, high spatial resolution, and high time resolution. Four GF-2 images covering the study area in July 2020 were used.

The data of POI comprises geographical objects that can be abstracted as points, especially some geographical entities closely related to people’s lives (Wenhao and Tinghua, 2015). In this study, the POI map data was divided into 18 service categories: automobile, motorcycle, catering, shopping, life, sports and leisure, healthcare, accommodation, scenic spots, commercial residences, government institutions and social organizations, science and education, transportation facilities, financial and insurance, company enterprises, road ancillary facilities, place name address information, and public facilities. Given road ancillary facilities and place names usually represent non-substantive information, such as road section identification, administrative place names, natural place names, and traffic place names, they were not selected.

The NPP-VIIRS data were obtained from the National Geophysical Data Center (NGDC) of the United States and were detected by the Suomi NPP satellite with a visible infrared imaging radiation instrument (Shi et al., 2014). The monthly data of lighting products in July 2020 were selected with a spatial resolution of 500 m.

The WorldPop project, hosted by the Department of Geography and the Institute of New Pathogens at the University of Florida, aims to provide spatial population data sets for Central and South America, Africa, and Asia to support development, disaster response, and health applications (Tatem, 2017). Table 1 lists the specific information of all the data used in this study.

SVM is a classification algorithm based on the VC dimension theory of statistical learning theory and the structural risk minimization principle (Cortes and Vapnik, 1995). Compared with traditional classification ideas based on neural networks or statistics, SVM controls the complexity of the model through the number of vectors and does not need to reduce the feature variables through dimensionality reduction processing to control the complexity of the model (Thaseen and Kumar, 2017). Therefore, in the classification process, the SVM classifier will not lose the feature information of ground objects and reduces the occurrences of some over-fitting phenomena.

The basic principle of SVM is to transform the original feature vector into a high-dimensional feature space and to solve the optimal classification hyperplane in a high-dimensional space. The classification accuracy depends on the distance between the hyperplane and the boundary plane on both sides of the hyperplane. The larger the distance, the higher the accuracy of the classifier and the smaller the error. Given a sample set ( x i , y i ), i is 1, 2, 3.N, x i ϵ R n , y i ϵ − 1 , 1 , the classification line is H = ω ∙ x + b = 0 . ω is the normal vector, x is the training sample, and b is the offset. The optimal classification line needs to have the maximum classification interval, that is, the maximum distance of 2 ω between the two types of samples, equivalent to finding the minimum value of the function. φ ω = 1 2 ω 2 (1)

The classification of land use belongs to the linear inseparable problem, so it is necessary to introduce the relaxation variable ξ i and the penalty parameter c into the above expression, which can then be expressed as: min ⁡ ⁡ 1 2 ω 2 + c ∑ i = 1 l ξ i + ξ i * , ξ i ≥ 0   i = 1 , 2 , … N s . t . ω x i + b − y i ≤ ε + ξ i y i − ω x i + b ≤ ε + ξ i * ξ i , ξ i * ≥ 0 (2) where ε is the insensitive loss coefficient, and its value affects the number of support vectors; ξ i and ξ i * are relaxation variables, indicating the extent to which the sample deviates from the ε-insensitive region.

For optimization problems, it needs to be transformed into the dual problem of the Lagrangian operator, and the final expression of SVM is: f x = ∑ i = 1 n α i − α i * k x i , x + b (3) where k x i , x is the kernel function, α and α * are Lagrange multipliers.

The object-oriented classification method breaks through the limitation of the traditional classification method which takes a single pixel as the basic classification and processing unit, classifies the image from the object level, and reduces the loss rate of semantic information contained in the traditional pixel-based classification method (Stumpf and Kerle, 2011). The object-oriented SVM classification method takes into account the advantages of object-oriented multi-scale segmentation and SVM. Firstly, multi-scale segmentation is carried out according to the properties of the object region on the image, which not only considers the spectral information of the image, but also adds features such as texture, geometric shape, and spatial topological relations, and then uses training samples for SVM classification. The object-oriented SVM classification method has obvious advantages in accuracy, generalization, and high-dimensional data processing, and thus has been widely used in RS image classification (Li et al., 2010; Lin et al., 2018; Shen et al., 2020). To obtain the land use types in the study area, this research utilized the object-oriented SVM classification method to classify remote sensing images.

The landscape disorder degree can indicate the degree of fragmentation and dispersion of urban landscapes, as well as reflect the heterogeneity and homogeneity of landscape space (Dong et al., 2010). The higher the heterogeneity of land use patches per unit area, the greater the degree of landscape disorder. Urban and rural areas usually have a single type of land, mostly for construction or agricultural use, and the degree of their landscape disorder is low. Since the urban fringe is an active expansion zone between the urban landscape and agricultural hinterland, it has various types of land use and a high degree of landscape disorder (Scott et al., 2013; Yang et al., 2017; Sharaf El Din, 2020). Therefore, the scope of the urban fringe can be determined by the differences in landscape disorder. The formula for the landscape disorder degree is as follows: W = − ∑ n = 1 N X n ⁡ ln X n ( Vranken et al ., 2015 ) (4)

Where W is the value of the degree of landscape disorder, X n represents the ratio of a certain type of land in the unit area, and n is the number of land use type patches per unit area.

The KDE calculation is often used to evaluate the density value of the neighborhood of point or line elements, and to simulate the spatial distribution of elements, which is widely used in geospatial analyses (Lin et al., 2021). Its main principle is that the estimated density of elements decreases with an increase in distance within a certain bandwidth. The kernel density at the center of the feature is the highest, and the kernel density at the edge of the bandwidth is 0. As KDE obeys the law of spatial correlation, the closer the distance is, the greater the correlation is, hence the POI data also conform to this law. The kernel density formula is as follows: λ s = ∑ l = 1 n 1 π r 2 φ d l s r (5) where λ s is the estimated density at region s , r is the bandwidth set by the kernel density function, n is the total number of elements involved in the calculation, d l s is the distance between POI points l and s , and φ is the weight of the distance.

Geodetectors are a group of statistical methods used to detect spatial diversity and explain the driving force behind it, including factor, interaction, risk, and ecological detectors (Wang et al., 2010). The main principle of the geodetector is to assume that the study area is divided into several sub-regions; if the sum of the variance of the sub-region is less than the total variance of the region, there is a spatial difference; if the spatial distribution of the two variables tends to be consistent, there is a statistical correlation between the two variables. Geodetectors can evaluate spatial differentiation, detect explanatory factors, and analyze the interaction between variables, and therefore have been widely used in nature, environmental science, human health, and other fields (Wang et al., 2010; Cao et al., 2013; Liu et al., 2020). Based on the factor detector, this study obtained the degree of landscape disorder, POI kernel density, and night light intensity to explain the spatial distribution of population.

Of the detector components, the factor detector detects the spatial differentiation of attribute Y and the explanatory power of factor X to attribute Y, which is measured by q value. Given population density is closely related to urban development (Liu et al., 2012; Yue et al., 2013; Li et al., 2019), this study determined the weight of each factor according to the explanatory power of landscape disorder degree, POI kernel density, and night light intensity on population spatial distribution. The value q can be defined as: q = 1 − ∑ h = 1 L N h σ h 2 N σ 2 = 1 − S S W S S T   h = 1,2 … S S W = ∑ h = 1 L N h σ h 2 , S S T = N σ 2 (6) where L is the stratification of variable Y or factor X, that is, classification or partition ;   N h and N are the numbers of units in layer h and the whole region, respectively; σ h 2 and σ 2 are the variances of layer h and Y values of the whole region, respectively; and S S W and S S T are the sums of the variance of the layer and the variance of the whole region, respectively.

Taking the landscape disorder degree, POI kernel density, and night light intensity as urban characteristic factors, this study proposes the FEM-SMU model (Figure 2). The spatial distribution of population density in an urban area is reflective of the differences between the urban center, urban fringe, and rural areas (Simon, 2008). The explanatory power of individual factors (landscape disorder degree, POI kernel density, night light intensity) on population distribution can reflect their explanatory power on the differences in urban center, urban fringe, and rural characteristics. Therefore, the relative weights of the various factors can be established based on the ratio of the explanatory power of a single factor to the sum of the explanatory power of all factors. The implementation process of the model is divided into the following three steps.

The first step was to calculate the degree of landscape disorder, POI kernel density, and night light intensity. For the calculation of landscape disorder degree: First, Envi 5.3 software was used to perform atmospheric correction, fusion, mosaic, and cropping operations on GF-2 images in the study area, then the object-oriented SVM method was used for land use classification, and finally, Eq. 4 was used to calculate the degree of landscape disorder. For the calculation of POI kernel density: First, the POI data were filtered and reprojected using ArcGIS 10.3 software, and then the kernel density tool was used to evaluate the kernel density. For the calculation of night light intensity: First, the NPP-VIIRS data of the study area were reprojected, denoised, and cropped using ArcGIS 10.3, and then the DN value of pixels was obtained to represent the night light intensity. It should be noted that the degree of landscape disorder, POI kernel density, and night light intensity should be normalized and resampled to a resolution of 10 m after calculation to ensure accurate analysis.

The second step was to build a comprehensive index. First, the Wordpop data were reprojected and cropped using ArcGIS 10.3. Second, the explanatory power of the landscape disorder degree, POI kernel density, and night light intensity on population spatial distribution was determined based on the factor detector. Finally, the comprehensive index was established after the weight was determined according to the ratio of the explanatory power degree of a single factor to all explanatory power degrees.

The third step was to extract the fringe of small- and medium-sized urban areas. First, the comprehensive index grid was divided into three categories by using the Natural Breaks Classification (NBC) method, which represents the urban core, urban fringe, and rural hinterland, and then the classified results were transformed into vectors and smoothed to obtain the range of the urban fringe. It is worth mentioning that this paper determines the NBC method as the classification method for the composite index through comparison of the accuracy of equal interval method, NBC method, and standard deviation method, and reference to literature (Xu et al., 2013; Xiao et al., 2020).

Based on the preprocessed GF-2 imagery of the study area, the land use was divided into four categories: vegetation (cultivated land, forest land, and grassland), construction land, water, and unused land, using the object-oriented SVM classification method (Figure 3A). It can be seen from the figure that the construction land was mainly distributed in the middle and south of the study area, the water was distributed along the southern boundary, the unused land was mainly distributed in the north, and the vegetation was mainly distributed in the northwest, northeast, and southeast. According to statistics, the vegetation area of the study area was 80.34 km², the construction land area was 60.83 km², the water area was 6.95 km², and the unused land area was 3.31 km².

To highlight the circular structure of landscape disorder in the study area, a grid of 100 × 100 m was constructed as the scale unit of spatial calculation. ArcGIS 10.3 was applied to calculate the area ratios of vegetation, construction land, water, and unused land in the cell grid, and the degree of landscape disorder in the study area was calculated using Eq. 4; (Figure 3B). As can be seen from the figure, the landscape structure characteristics of the urban core area were prominent, the landscape disorder degree was low, and there were concentrated and contiguous low-value areas. After repeated experiments, a threshold value of less than 0.46 was used as the symbol to identify the urban core area. However, if there was a large range of green spaces in the urban core area, the degree of landscape disorder was higher, and the landscape disorder degree threshold method could not identify the complete urban core area. The difference in the degree of landscape disorder between the urban fringe and rural hinterland was not obvious, and the degree of landscape disorder was higher. We consider that, compared with larger cities, small- and medium-sized cities have fewer populations and smaller villages; further, the study area is located in the Hanzhong Plain, where villages are relatively concentrated and cultivated land is relatively scattered. These circumstances resulted in a higher degree of landscape disorder between the urban fringe and rural hinterland. Although the landscape disorder degree threshold method is currently widely used in the identification of metropolitan fringes (Huang et al., 2016; Wang Y. et al., 2021; Liu, 2021), in this study, the common shortcomings of this as a single-factor threshold method, which include discontinuity, lack of detail, and poor versatility, were amplified obviously in the process of urban fringe identification in the study area. Therefore, we conclude that the landscape disorder degree threshold method is not suitable for the identification of small- and medium-sized urban fringes.

The setting of the bandwidth in KDE has an important influence on the results. After referring to previous research results (Heidenreich et al., 2013; Lin et al., 2021), the distances of 500 m, 1000 m, and 1,500 m were selected for KDE analysis (Figures 4A–C). It can be seen from the figure that with the shortest bandwidth of 500 m, the KDE results were fragmented and discontinuous, and the overall distribution of urban POI was not obvious. On the other end of the spectrum, with a bandwidth of 1,500 m, the local characteristics of the overall distribution trend of urban POI were difficult to show, and the details were insufficient. Finally, with the medium bandwidth of 1,000 m, the KDE results had good stability and an obvious overall distribution, which can meet the analysis needs of the urban fringe in the study area.

Taking this selected bandwidth, Figure 4D shows the result of dividing the POI kernel density with a bandwidth of 1,000 m into three categories by the NBC method. As can be seen from Figures 4B, D, the POI kernel density showed an obvious circular structure distribution. The urban core area had a large area of continuous high-density area, the density of the urban fringe was low, and the density of most rural areas was close to 0. Compared with the landscape disorder degree threshold method, the results of the POI kernel density breakpoint analysis method can clearly and completely identify the urban core area; however, both of them had poor recognition results for the urban fringe area. Moreover, the POI kernel density breakpoint analysis method identifies well-developed villages with a large number of POIs as urban fringes. As the urban fringe was in the stage of development, with less POI data and slow updating speed, this led to obvious errors between the results and the actual urban fringe. Different from larger cities, the POI data integrity of small- and medium-sized urban fringe was lower and the update speed was relatively slow, so we conclude that the POI kernel density breakpoint analysis method had too much difficulty in accurately extracting the urban fringe of the study area.

Finally, the urban fringe of the study area was extracted using the FEM-SMU model proposed in this study. First, the degree of landscape disorder, POI kernel density, and night light intensity were calculated. Then, the weight of each factor was determined by combining the Geodetector and WorldPop data (landscape disorder degree: 0.10, POI kernel density: 0.51, night light intensity: 0.39) to construct a comprehensive index. Finally, the urban fringe was identified by the NBC method, and the results were post-processed.

Figure 5A shows the result of dividing the comprehensive index into three categories by the NBC method, and Figure 5B shows the final result of the FEM-SMU model. As can be seen from the figure, the FEM-SMU model could accurately and completely identify the urban fringe of the study area. The urban fringe was mainly concentrated in the north and east, with an area of approximately 39 km². Compared with the two single-factor extraction methods, the performance of the FEM-SMU model was a great improvement, especially for the results of the outer boundary of the urban fringe. The overall pattern of the inner boundary of the urban fringe (urban core area) extracted by the FEM-SMU model and POI kernel density breakpoint analysis was more consistent, and further, the former was more detailed. The difference between the results of the two methods was mainly concentrated in the southwest in the Hanjiang New Area which is still in the stage of rapid development with rapidly changing landscape patterns, and the POI data updated slowly. Therefore, there is a partial error between the inner boundary of the urban fringe extracted by the POI kernel density breakpoint analysis and the actual boundary. Nevertheless, the FEM-SMU model focuses on the comprehensive performance of the regional landscape disorder degree, POI kernel density, and night light intensity, and is less dependent on the performance of a single factor, thus it can more accurately identify the urban fringe.

As the landscape disorder degree threshold method makes it difficult to identify the urban fringe of the study area, this study only evaluates the accuracy of the POI kernel density breakpoint analysis method and FEM-SMU model. Through detailed visual analysis in Section 4, we found that the FEM-SMU model was significantly stronger than the POI kernel density breakpoint analysis method in detail and integrity. To further evaluate the extraction accuracy of different methods, field verification and landscape pattern index evaluations were used in this study. A total of 100 sample points were evenly selected along the road around the urban fringe (Figure 5B), and the accuracy of the extraction results was analyzed by field verification (Table 2). It can be seen from the table that the overall accuracy of the FEM-SMU model was significantly higher than that of POI kernel density breakpoint analysis method, reaching 98%. The POI kernel density breakpoint analysis method had a relatively large number of errors, with an overall accuracy of only 67%.

As another validation method, the landscape pattern index is often used to evaluate the extraction accuracy of urban fringes (Peng et al., 2016). After referring to previous research results (Vizzari and Sigura, 2015; Yang et al., 2021), patch density (PD) and the Shannon diversity index (SHDI) were selected to evaluate the accuracy of the two methods at the grade and landscape levels. The PD value indicates the degree of landscape fragmentation, and SHDI represents the degree of richness and complexity of landscape types. Generally speaking, the PD and SHDI in urban fringes should be higher, while those in urban centers and rural areas should be lower. Table 3 shows the PD and SHDI values of the two methods in different regions calculated by the Fragstats 4.2 software. It can be seen from the table that the PD and SHDI values of the FEM-SMU model were significantly higher than those of the POI kernel density breakpoint analysis method in the urban fringe. This shows that the degree of landscape fragmentation, complexity, and diversity in the urban fringe extracted by the former was higher than that of the latter, which can better reflect the characteristics of the urban fringe. In the rural hinterland, the PD and SHDI values of the FEM-SMU model were significantly lower than those of the POI kernel density breakpoint analysis method. However, the PD and SHDI values are approximate in the urban core areas because of the relatively small difference between the boundaries extracted by the two methods. In conclusion, the two evaluation methods show that the FEM-SMU model has higher accuracy and can accurately extract the urban fringe of the study area.

To further verify the applicability of the model proposed in this study in different regions and different types of small- and medium-sized cities, urban fringe identification was completed for Shangzhou District in Shangluo City and Hanbin District in Ankang City (Figure 6). Shangzhou District is a typical banded urban structure, which is severely limited by resource conditions. As shown in Figure 6A, the extraction result of the POI kernel density breakpoint analysis method is incomplete, especially in the southeast. The reason for this is likely that the southeast part of Shangzhou District is mainly an industrial park, and the number of POIs is small and scattered, which is not enough to support the identification of the urban fringe. Hanbin District is separated from the middle by the Hanjiang River, with the old urban area in the southeast and the new urban area in the northwest, which is a multi-center urban structure. As can be seen from Figure 6B, the extraction result of the POI nuclear density breakpoint analysis was poor, especially in the old urban area of the southeast. The main reason for this is likely that the single-factor method requires high data quality when extracting the fringe of multi-center cities, while the development of the old urban area is relatively backward, the population distribution is concentrated, and the POI data are not complete. As a better option, the FEM-SMU model identifies the urban fringe according to the comprehensive performance differences of regional landscape disorder degree, POI kernel density, and night light intensity combined, and is less dependent on the performance of a single factor, thus it can adapt to different regions and different types of small- and medium-sized cities.