A Geomorphological Regionalization using the Upscaled DEM: the Beijing-Tianjin-Hebei Area, China Case Study

Characterizing geomorphological patterns based on digital elevation models (DEMs) has become a basic focus of current geomorphology. A new DEM upscaling method based on the high-accuracy surface modelling method (HASM-US method) has been developed to improve the accuracy of current models and the subjectivity of macroscopic geomorphological patterns. The topographic variables of elevation (EL), slope (SL), aspect (AS), relief amplitude (RA), surface incision (SI), surface roughness (SR), and profile curvature (PC) with a spatial resolution of 1 km × 1 km in the Beijing-Tianjin-Hebei (BTH) area of China have been obtained by using the HASM-US method combined with the principal component analysis (PCA) method in terms of the elevation data of the SRTM-4 DEM, meteorological station location information, and field measurements with a GPS receiver. A geomorphological regionalization pattern has been developed to quantitatively classify the geomorphological types in the BTH area by combining the seven topographic factors of EL, SL, AS, RA, SI, SR, and PC that have significant spatial variation. The results show that the upscaling accuracy of elevation (mean difference only −2.32 m) with the HASM-US method is higher than that with the bilinear interpolation method and nearest neighbour interpolation method. The geomorphologic distribution in the BTH area includes 11 types: low plain, low tableland, low hill, low basin, middle plain, middle hill, low mountain with low RA values, low mountain with medium RA values, middle mountain with low RA values, middle mountain with medium RA values, and middle mountain with high RA values. The low plain is the dominant geomorphological type that covers 40.58% of the whole BTH area. The geomorphological distribution shows the different significant characteristics: the elevation rapidly decreases from the Taihang Mountains to the eastern area, gradually decreases from the Yanshan Mountains to the southern area, and first increases and then decreases from the Bashang Plateau to the southeastern area in the whole BTH area.

www.nature.com/scientificreports www.nature.com/scientificreports/ can be classified on the basis of slope, relief, and profile types using statistical windows 1,[13][14][15] . It remains the case, however, that topographic factors will not be exactly the same in all applications; Dragut 16 used a suite of four parameters, elevation, profile and plan curvature, as well as slope gradient, while Liu utilized six parameters: relief amplitude, surface incision, surface roughness, and the elevation variance coefficient, as well as the mean slope and elevation 7 . In contrast, Hu 7 selected just three parameters: relief amplitude, surface incision, and the terrain position parameter (TPI). Large-scale geomorphological patterns can mostly be extracted from SRTM DEM data 4,17,18 ; indeed, a great deal of attention has been afforded to these data, as this application boasts near-global coverage (i.e., between 56°S and 60°N), has relatively high spatial resolution and is free of charge. SRTM DEM data outputs have therefore been applied widely 9,19 ; in terms of regional-scale digital geomorphology, SRTM data have commonly been upscaled to obtain macroscopic geomorphological information 20 . The accuracy of DEM data is mainly controlled by the measurement accuracy, resolution of DEM data, the slope of a pixel, and modelling error of DEM model 20 , which depends to a large extent on sample data reliability 21 and therefore means that simple resampling methods can significantly increase errors 19 . The traditional DEM models in this area suffers from a series of inherent shortcomings, including subjective parameter selection for depicting macrogeographical patterns as well as a greater level of error and the loss of detailed information from DEM data subsequent to upscaling processing. Therefore, on the same condition of the measurement accuracy, terrain complex level and resolution of DEM data, how to develop a suitable upscaling method to improve the DEM modelling accuracy is important in the geomorphological Regionalization.
The Beijing-Tianjin-Hebei (BTH) region is located between 113°04′-119°53′E and 36°01′-42°37′N within China (Fig. 1). Administrative zoning in this area includes the municipalities of Beijing and Tianjin as well as Hebei Province, encompasses a total area of 21.60 × 10 4 km 2 22 and is the political, cultural, and economic centre of China. This region also includes the Taihang and Yanshan Mountains as well as the North China Plain, is located within the transitional zone between the second and third steps and thus contains both complex and diverse geological structures and types. The BTH region, therefore, has important political and economic status and occupies a very important physical geographical region within China. There has also been no research to date that specifically addresses the macroscopic geomorphological features of this area. The main objective of this paper is therefore to propose a new upscaling method and to reveal the basic characteristics and spatial patterns of geomorphology within the BTH region. These results not only reveal the macroscopic landform features of www.nature.com/scientificreports www.nature.com/scientificreports/ this region but also provide a scientific basis for the sustainable utilization of regional resources, urban and rural development, and industrial and agricultural layout.

Data and Methods
Data sources. The data sources used in this study include three categories. The first was SRTM-4 grid data (CGIAR-CSI, http://srtm.csi.cgiar.org/) with a spatial resolution of 3 arcsec (approx. 90 m 2 at the Equator); multiple sets of these data were downloaded and projection transformed (as in the continuous operation reference station system, CORS). The second category of data comprised coordinate data from 102 meteorological stations at elevations less than 50 m within the BTH region but with large positional error eliminated. The third category comprised the elevation data of 458 field points with information on altitude and geographic coordinates observed using GPS-based from CORS.

HASM-US method.
A method for high accuracy surface modelling (HASM) [23][24][25][26] has been developed since 1986 to integrate the extrinsic and intrinsic properties and find solutions for the error and multi-scale problems that are prevalent when each type of information is used separately 27 . This approach uses these inputs as its drivers as well as locally accurate information (e.g., ground observation and/or sampling data) as optimum control constraints. As this modelling approach is based on differential geometry and optimal control theories, it has been widely applied across numerous fields, including in Earth surface modelling DEMs 19,[26][27][28] .
Applying the HASM algorithm, a HASM-upscaling (HASM-US) algorithm for DEM upscaling simulation was developed in this study.
Grids with different spatial resolutions within the same region, Ω H ( ) and Ω h ( ), are expressed as follows: Thus, assuming the grid resolution is .., the sample set is П, and the HASM algorithm based on the equality constraint is ∏ HASM h ( , ), as follows:  Fig. 2(a)). Thus, during the upscaling process, a coarse-resolution DTM retains its overall trend and encompasses additional details by optimizing the control point data. www.nature.com/scientificreports www.nature.com/scientificreports/ A projection was formed between coarse and fine grids throughout the process of spatial grid data upscaling. The coarse grid comprises the following 16 fine grids ( Fig. 2(b)): The fine grid can therefore be simplified as follows: A simple arithmetic averaging method was then adopted to reduce computation, and a value of X I J ( , ) H was obtained for more accurate simplification using the inverse distance interpolation method. The upscaling process of each coarse grid can therefore be similarly simplified; thicalculation process can be defined as a grid projection operator G h H , as follows: As X H is the trend surface, it is possible to calculate the residual of the sample point set П and the trend surface, which enables HASM simulation to be performed on residuals. A grid residual surface was therefore obtained by superimposing the trend surface, X I J ( , ) H , to obtain the final DTM. This upscaling algorithm was defined as HASM-US and was formulated as follows: Data processing. The ellipsoidal heights (He) derived from GPS measurements were converted to orthometric heights (Ho) of EGM96 used in the SRTM: Ho = He -N. Here, N is the geoid height. During the upscaling simulation process of SRTM data, the elevation value of measuring point by GPS instrument and meteorological station were used as the basic value, and the differences between the corresponding pixel value of SRTM data / upscaled DEM data and the basic value were calculated and used to validate the modelling errors of upscaled SRTM data with HASM-US and other methods. The elevation differences (Table 1) were calculated between 560 pixel values of SRTM sample data and the basic values, which correspond to 458 pixels of field measuring points and 102 pixels of meteorological stations, respectively. The calculated results show that the maximum positive difference and maximum negative difference in this analysis were 48.61 m and −38.41 m, respectively, encompassing a mean difference of 1.81 m, a root mean square (RMS) difference of 10.75 m, and a mean absolute difference of 7.65 m. The SRTM elevation values tend to be higher because measurement points are mainly distributed within the northern and western mountainous areas where RA is large; for this reason, an difference statistical SRTM analysis was performed using meteorological station elevation data. This additional analysis revealed that the overall difference is not large, especially within the flat plain area where the main meteorological stations are located. The data in Fig. 3 reveal the distribution of DEM grid value distributions at each sample point; within the northern and western mountainous areas of the BTH region, the SRTM data include relatively high elevation values, while relatively low elevation values are seen in the northeastern mountainous, central, and southeastern plains. DEM data is one of the key parameters for many eco-environmental models. Except the DEM data, most eco-environmental models usually involve many other parameters such as vegetation, soil, climate and human activity data In fact, even though there is a high-resolution DEM data, if the high-resolution data of other parameters needed by the eco-environment model is very difficult to obtain, the high-resolution DEM data is often need to be upscaled to ensure union all parameters in the eco-environmental model 27 . Therefore, to overcome problems of missing detailed information and the loss of accuracy when directly resampling the DEM, SRTM grid data were scaled up using the high-accuracy surface modelling (HASM) method 19 . Combined with meteorological station coordinates and field-measured GPS data, DEM data for 1 km by 1 km cells were obtained; combined with meteorological station coordinates and field-measured GPS data, a series of parameters (i.e., elevation, slope, aspect, RA, SI, SR, profile curvature) were calculated using ArcGIS software. Applying the mean change-point analysis method 29 , it is possible to conclude that a grid size of 5 km by 5 km comprises the best window size for www.nature.com/scientificreports www.nature.com/scientificreports/ extracting BTH region RA from 1 km-resolution DEM data. Principal component analysis (PCA) was then used to screen dominant topographic factors to classify geomorphological types within the study area.
Classification method of geomorphological types. Terrain analysis is a prerequisite for geomorphological regionalization. Due to the complicated landforms in the study area, diverse common parameters should be selected to characterize the topographic features and geomorphological regionalization. According to   30 , Hu (2015) 7 , and Li (2013) 31,9 , the parameters include elevation (EL), slope (SL), aspect (SA), relief amplitude (RA), surface incision (SI), surface roughness (SR) and profile curvature (PC) (Fig. 4).
Geomorphological studies carried out based on DEMs have led to a number of different correlations between parameters. Thus, to determine which parameters are most important, principal component analysis (PCA) was carried out (Tables 2-4). The cumulative contribution rate of the first two principal components (Pc1, Pc2) can reach up to 97.60%, while the highest correlation of the first Pc is elevation (0.950) and the second is RA (0.877). Elevation and RA were therefore used as parameters for classifying geomorphological types within the study area.
The geomorphological regionalization scheme (Table 5) of the BTH area was approved by integrating with the actual situation, EL and RA parameters in the whole BTH area. The landform of the BTH area was classified into 11 geomorphological types: low plain, low tableland, low hill, low basin, middle plain, middle hill, low mountain with low RA values, low mountain with medium RA values, middle mountain with low RA values, middle mountain with medium RA values, and middle mountain with high RA values.

Results
Accuracy validation. To verify the advantages and disadvantages of the upscaling algorithm used here combined with sample information and common resampling methods, the methods of cross-validation, bilinear interpolation, nearest neighbour interpolation, and HASM-US were all used to classify geomorphological types within the BTH study area. A total of 560 sampling points were randomly divided into ten components, of which nine     www.nature.com/scientificreports www.nature.com/scientificreports/ the corresponding differences of the results from the other two kinds of conventional resampling are approximately −6.90 m, 28 m, and 15 m, respectively. The accuracy of the coarse-resolution DEM can be effectively improved by introducing sample data and performing DEM upscaling using the HASM-US algorithm. As shown in Tables 6-8, the mean difference is greater than zero, and it is therefore the case that the elevation values of the upscaling results in each method are slightly higher than the measured values.    www.nature.com/scientificreports www.nature.com/scientificreports/ In addition to cross validation, this study also qualitatively evaluated the effects of DEM upscaling simulations using cross-section line analysis. Thus, two elevation section lines were extracted in the northern and western regions of the BTH area (Fig. 5). Section lines of the relief were revealed as profiles by SRTM 4 DEM data, including the 1 km spatial resolution DEM profile obtained by bilinear interpolation, the nearest neighbour interpolation method, and the HASM-US method. The simulation results of HASM-US can therefore better capture extreme DEM points on the basis of section lines.
Geomorphological Regionalization in the BtH Area. Geomorphological types within low-altitude areas (<1,000 m) include six categories (Table 9, Fig. 6). The low plain area within BTH comprises approximately 8.98 × 10 4 km 2 , 41.75% of the total area. The eastern and southern regions extend to the southeastern province boundary as well as to the coastline of Hebei Province, west to the eastern fringes of the Taihang Mountains, and north to the southern fringes of the Yanshan Mountains. On the western and northern sides of this region, the extent of this line is consistent with numerous highways and national highways and encompasses multiple towns    www.nature.com/scientificreports www.nature.com/scientificreports/ and cities. The low-altitude tableland area within the BTH area encompasses 6,773 km 2 , 3.32% of the study area. These low hills are widely distributed and encompass a total area of approximately 3.20 × 10 4 km 2 , 14.89% of the total area. The spatial pattern of this land type is a belt: the frontal mountain transition zone to the east of the Taihang Mountains as well as the transition zone to the south of the frontal Yanshan Mountains encompass areas of approximately 1.55 × 10 4 km 2 and 1.65 × 10 4 km 2 , respectively. Low basins within the region encompass an area of 4,273 km 2 , 1.99% of the total area. The low mountain area with low RA values is 2.52 × 10 4 km 2 , 11.70% of the total area, and includes the Taihang Mountain area south of Xiaowutai Mountain. This land area comprises a discontinuous band and is more continuous in the northeast, accounting for a large proportion of the Yanshan Mountain area. Low mountains with medium RA areas encompass approximately 9,884 km 2 , 4.59% of the total area.
The middle plain and middle hill regions comprise the main body of the Bashang Plateau. The first of these two regions comprises an area of 5,497 km 2 (2.55%), while the latter region covers an area of 9,481 km 2 (4.57%). The middle hill regions within the BTH area are mainly located in places such as the transition zone of Damaqun Mountain as well as on the middle plain. The middle mountain area with low RA values also comprises 2.21 × 10 4 km 2 , 10.28% of the area. The middle mountain area within the BTH area that has medium RA values comprises 9,085 km 2 , 4.22% of the total area. This kind of land is widely distributed but scattered, including on the eastern side of the Taihang Mountains, the eastern part of Xiaowutai Mountain, and in the middle-western region of the Yanshan Mountains. The area of middle mountain land with high RA values comprises approximately 263 km 2 , 0.12% of the total area. Geomorphological features Analysis in the BtH Area. Geomorphological patterns within the BTH region reveal a basic framework that comprises more low plains in the southeast as well as more mountainous areas in the west and north. The southern section of the western part of the study area is rapidly decreasing in altitude from the Taihang Mountains to the east onto low plains. The profile curve for the BTH area is therefore www.nature.com/scientificreports www.nature.com/scientificreports/ jagged and is characterized by an obvious downward trend. Low-altitude hills in the mountain fronts are wide; mountain basins and high-amplitude mountains interlace with one another and are distributed in the northern section of the western part of the BTH area. Subsequent to the highest Xiaowutai Mountain peak, landforms rapidly decrease to a low plain, while the low hilly area in front of the mountain is narrow. Moving from the middle plains and hilly areas on Bashang in the northwest to the southeast, Damaqun Mountain and other areas increase in amplitude before gradually falling to low plains, while wide valleys develop between mountains. The relief of the Yanshan Mountain district in the north generally decreases from north to south; the profile line across this area is extremely tortuous and indicates strong tectonic movement across this area. The terrain is also relatively fragmented within this region because of river erosion.

Discussion
DeM Upscaling methods. Upscaling conversion of DEM data is necessary before they are used to join the simulation and analysis with other data at different scales. The methods of upscaling DEM data mainly include classical resampling methods, terrain simplification based on the LOD model 32 , the Douglas Peucker algorithm 33 , the point spread function 34 , the filtering method 35 , wavelet transform 36 , fractal 37 , etc., which are mainly used for simplification of the terrain information of scattered points; some are mainly used for DEM scale conversion of TIN format, but the most widely used method is resampling.
The general resampling method (e.g., the nearest neighbour sampling and bilinear interpolation methods), as a grid data processing method, is the most classical method of upscaling transformation of grid data, which is supported by mainstream GIS software and widely used in remote sensing image processing and DTM spatial analysis 38,39 . The accuracy of the DEM will decrease with increasing grid resolution when the DEM data are upscaled, and different resampling methods will have different effects on the accuracy of upscaled DTM data 40 .
The measurement accuracy, resolution of DEM data, terrain relief, pixel size, and modelling error of DEM model directly affect the accuracy of DEM data 20,41 . On the same condition of the measurement accuracy, terrain complex level and resolution of DEM data, the modelling error of DEM is mainly decided by the adopted upscaling model. However, in order to reduce the amount of data and maintain more terrain details in the upscaled DEM data, some elevation samples are usually added to participate in the upscaling transformation of the DEM in the classical upscaling methods [42][43][44] . Moreover, the resampling method, wavelet transform, fractal and other classical methods can't deal with upscaling the DEM data when some new sample data are joined for simulation. The HASM-US method developed in this paper can not only involve these new sample data in the upscaling process of DEM but also decrease the accuracy loss of the upscaled DEM data and retain more terrain details in the original DEM data Geomorphological patterns in the BtH Area. Terrain is one of the most fundamental components of the natural environment 31,45 . Geomorphological regionalization is key to the study of the spatial differentiation of geographical phenomena 31 . The similarities and differences of the combination of landform types and the causes of landforms are the basis of geomorphological zones at all levels 31 . There are various methods for the classification of landform types, but they have similar characteristics in general. Based on remote sensing and GIS methods, Boccao (2001) 46 used the parameters of relief amplitude, slope steepness and dominant lithology to classify major landforms. Gao (2019) 47 made geomorphological divisions of China on the basis of regional differences and the formation of landforms.   48 divided the landform zones of the North China Plain according to morphology and formation. Cheng (2017) 49 characterized the geomorphological zones on the northwestern edge of the Qinghai-Tibet Plateau with topographic profiles and gradients based on digital elevation model (DEM) data. With remote sensing images and SRTM-DEM and other multisource data, intelligent regionalization can be realized using the methods of geographical grid and systematic cluster analysis 50 . Therefore, with the development of remote sensing and geographic information technology, rapid regionalization can be achieved based on digital terrain analysis (DTA) of DEMs. Regarding the geomorphological regionalization of the BTH area, Zhao and Cheng (2016) first divided the landform types into plains and mountains based on remote sensing images, DEMs and historical geomorphological maps and then calculated the relief amplitude of the SRTM-DEM and classified it as plains, tablelands, hills, and mountains with low, medium and high relief amplitudes 51 . In this paper, 2 primary parameters were selected from 7 topographic parameters based on the DTA of more accurate DEMs and PCA and then classified as the quantitative criterion for geomorphological regionalization, in which the subjectivity can be reduced to the greatest extent. Compared to Zhao's program, the grading of elevation and subdivisions of plains and hills produced 11 types; the diverse geomorphological regions in our study area can characterize the topographic features in detail.

conclusions
Multisource elevation data were upscaled in this analysis to obtain DEM grid data with a resolution of 1 km based on the HASM-US method; seven geomorphological parameters, including EL, SL, SA, RA, SI, SR, and PC, were calculated using ArcGIS software. Dominant geomorphological factors were extracted via PCA for geomorphological classification. This analysis leads to a number of key conclusions.
First, cross validation and surface profile analyses reveal that HASM-US enables higher accuracy and better reflects extreme values within DEMs in terms of upscaling simulation results when compared to those of either bilinear interpolation or nearest neighbour interpolation methods. HASM-US is therefore a more efficient method for upscaling DEMs.
The application of PCA reveals that elevation and RA provide representative factors for characterizing regional geomorphological features. Thus, combining these two, the BTH study area can be divided into 11 geomorphological types: low plains, low tablelands, low hills, low basins, low altitude mountains with low RA values, low altitude mountains with medium RA values, middle plains, middle hills, middle altitude mountains with low RA www.nature.com/scientificreports www.nature.com/scientificreports/ values, middle altitude mountains with medium RA values, and middle altitude mountains with high RA values. Low plains are the major land type within the BTH region and account for 40.58% of the total area. Trends in amplitude changes vary greatly along east-west, north-south, and northwest-southeast transect directions.