2.1.2. Spatial Scope

The Sanhe region is composed of the Hedong region, Henan region and Henei region. Before the Qin and Han dynasties, the Sanhe region was only used as a generalized geographical location [12], and whether there was a county or not in the pre-Qin period does not have clear evidence in academic research. Until the Qin and Han Dynasties, there was a clear administrative boundary. In the Qin Dynasty, the region was named Hedong County [13], Hanei County and Sanchuan County, chronologically, and during the Western Han Dynasty, it was named Hedong County, Hanei County and Henan County, chronologically [14]. Compared with the Qin Dynasty, the area of the Sanhe region in the Western Han Dynasty had expand. In this study, the research boundary is determined by the Sanhe region map in the Qin and Han Dynasties depicted in the "Atlas of Chinese History and the Records of the Sanhe Region in the pre-Qin Period" written by Tan Qixiang (Figure 1), which roughly overlaps the city area of Yuncheng City and Linfen City in Shanxi Province (namely the ancient Hedong region), Luoyang City and Zhengzhou City

in Henan Province (namely the ancient Henan region), Jiyuan City, Jiaozuo City, Xinxiang City, Anyang City and Hebi City (namely the ancient Henei region) [15]. Xinxiang City, Anyang City and Hebi City (namely the ancient Henei region) [15].

*Water* **2022**, *14*, x FOR PEER REVIEW 3 of 17

Western Han Dynasty had expand. In this study, the research boundary is determined by the Sanhe region map in the Qin and Han Dynasties depicted in the "Atlas of Chinese History and the Records of the Sanhe Region in the pre-Qin Period" written by Tan Qixiang (Figure 1), which roughly overlaps the city area of Yuncheng City and Linfen City in Shanxi Province (namely the ancient Hedong region), Luoyang City and Zhengzhou City in Henan Province (namely the ancient Henan region), Jiyuan City, Jiaozuo City,

**Figure 1.** Location of the Sanhe region in the Yellow River Basin.

#### **Figure 1.** Location of the Sanhe region in the Yellow River Basin. *2.2. Data*

recorded.

research results.

*2.2. Data*  The sites in this study are mainly chosen from the list of "China national key cultural relics protection units", provincial cultural relic protection units published on the website of the Shanxi Provincial Cultural Relics Bureau (http://wwj.shanxi.gov.cn/, accessed on 20 May 2021), and municipal cultural relic protection units, county-level cultural relic protection units, and unclassified cultural relic protection units in the Chinese Cultural Relics Atlas—Shanxi volume; and the Henan Provincial Cultural Relics Bureau website (http://wwj.henan.gov.cn/, accessed on 20 January 2022) announced the national key cul-The sites in this study are mainly chosen from the list of "China national key cultural relics protection units", provincial cultural relic protection units published on the website of the Shanxi Provincial Cultural Relics Bureau (http://wwj.shanxi.gov.cn/, accessed on 20 May 2021), and municipal cultural relic protection units, county-level cultural relic protection units, and unclassified cultural relic protection units in the Chinese Cultural Relics Atlas—Shanxi volume; and the Henan Provincial Cultural Relics Bureau website (http://wwj.henan.gov.cn/, accessed on 20 January 2022) announced the national key cultural relics protection units in Henan Province, as well as the volumes on the immovable cultural relic lists of Zhengzhou, Luoyang, Anyang, Hebi, Xinxiang, Jiyuan, and Jiaozuo from the third Henan Provincial Cultural Relics survey. Among them, a total of 478 Paleolithic sites, a total of 2137 Neolithic sites, and a total of 2318 Xia–Shang-Zhou sites are recorded.

tural relics protection units in Henan Province, as well as the volumes on the immovable cultural relic lists of Zhengzhou, Luoyang, Anyang , Hebi, Xinxiang, Jiyuan, and Jiaozuo from the third Henan Provincial Cultural Relics survey. Among them, a total of 478 Paleolithic sites, a total of 2137 Neolithic sites, and a total of 2318 Xia–Shang-Zhou sites are The DEM elevation data was obtained from Google Maps. In addition, the slope, aspect, and topographic relief data were also based on the DEM elevation data and calculated by GIS software. The data used in the research is raster elevation data. As the research area is very large, high-precision data cannot be calculated in GIS software, so DEM elevation data with a 9.5 m precision was selected. In fact, for a large research area, existing research generally uses 30 m precision data provided by government websites, and its accuracy is

The DEM elevation data was obtained from Google Maps. In addition, the slope, as-

lated by GIS software. The data used in the research is raster elevation data. As the research area is very large, high-precision data cannot be calculated in GIS software, so DEM elevation data with a 9.5 m precision was selected. In fact, for a large research area, existing research generally uses 30 m precision data provided by government websites, and its accuracy is difficult to guarantee. Therefore, the highest-precision data on the research area that could be operated by computers was selected to improve the accuracy of the difficult to guarantee. Therefore, the highest-precision data on the research area that could be operated by computers was selected to improve the accuracy of the research results.

#### *2.3. Methods*

### 2.3.1. Adjacent Index

The spatial distribution types of the sites in the study area are determined through adjacent index analysis of the sites. The adjacent index analysis results will return five values: average observation distance, theoretical average distance, nearest neighbor index, *Z* score and *p* value. *p* value represents probability; when the *p* value is lower, the observed spatial pattern is unlikely to be generated in a random process (small probability). The *Z* score is a multiple of the standard deviation. The higher the the *Z* value, the greater the degree of aggregation. The classification rules of the spatial distribution types of sites are as follows: if the adjacent index is less than 1, the spatial distribution type belongs to the aggregation distribution type; if the nearest neighbor index is equal to 1, the spatial distribution type is random; if the adjacent index is greater than 1, it belongs to uniform distribution type. The calculation formula is as follows:

$$R = \frac{\overline{r}\_1}{r\_E} = 2\sqrt{Dr\_1} \tag{1}$$

In the formula, *r*<sup>1</sup> refers to the average value of the distance *r*<sup>1</sup> between the nearest neighbor points, and *r<sup>E</sup>* refers to the average distance of the nearest neighbor points in the random distribution pattern. *D* refers to the density of all points, and finally the adjacent index *R* is calculated [16,17].

## 2.3.2. Kernel Density Estimation

The kernel density estimation method can reflect the spatial distribution and aggregation characteristics of site points. The larger the kernel density, the denser the distribution of sites. In addition, the kernel density map of geodiversity elements and the distribution of cultural sites were superimposed, and the correlation between them was discussed. The calculation formula is as follows:

$$f\_{\boldsymbol{\mu}}(\boldsymbol{x}) = \frac{1}{\text{nh}} \sum\_{i=1}^{n} k(\frac{\boldsymbol{x} - \boldsymbol{X}\_{i}}{h}) \tag{2}$$

In the formula, the kernel density is estimated as the probability that the density function *f* is at the point *x*, where *k*( *x*−*X<sup>i</sup> h* ) is the kernel function, *h* represents a search radius greater than 0, and *x* − *X<sup>i</sup>* means the distance (km) between the estimated point *x* and the event *X<sup>i</sup>* [18,19].

#### 2.3.3. Standard Deviation Ellipse

The distribution direction of the sites in each period is determined by the standard deviation ellipse. The long axis of the ellipse represents the main direction of the site distribution, the short axis represents the secondary direction of the site distribution, and the inclination angle of the ellipse represents the distribution trend of the sites. By comparing the distribution trend and direction of the sites in each period, the evolution characteristics of their spatial and temporal distribution can be obtained. The calculation formula is as follows:

$$\begin{array}{l} SDE\_X = \sqrt{\frac{\sum\_{i=1}^n (x\_i - \overline{x})^2}{n}}\\ SDE\_y = \sqrt{\frac{\sum\_{i=1}^n (y\_i - \overline{y})^2}{n}} \end{array} \tag{3}$$

In the formula, *x<sup>i</sup>* and *y<sup>i</sup>* are the coordinates of points, *n* is the total number of samples, and {*x*,*y*} is the average center of all points.

$$\begin{array}{c} \tan \theta = \frac{A + B}{C} \\ A = \left( \sum\_{i=1}^{n} \widetilde{\mathbf{x}}\_{i}^{2} - \sum\_{i=1}^{n} \widetilde{y}\_{i}^{2} \right) \\ B = \sqrt{\left( \sum\_{i=1}^{n} \widetilde{\mathbf{x}}\_{i}^{2} - \sum\_{i=1}^{n} \widetilde{y}\_{i}^{2} \right)^{2} + 4 \left( \sum\_{i=1}^{n} \widetilde{\mathbf{x}}\_{i} \widetilde{y}\_{i} \right)^{2}} \\ \mathbf{C} = 2 \sum\_{i=1}^{n} \widetilde{\mathbf{x}}\_{i} \widetilde{y}\_{i} \end{array} \tag{4}$$

In the formula, tan *<sup>θ</sup>* is the tangent of the ellipse rotation angle, and *<sup>x</sup>*e*<sup>i</sup>* , *y*e*i* is the deviation of *xy* coordinates from the average center.

$$\begin{aligned} \sigma\_{\mathcal{X}} &= \sqrt{\frac{\sum\_{i=1}^{n} \left(\tilde{\mathbf{x}}\_{i} \cos \theta - \tilde{\mathbf{y}}\_{i} \sin \theta\right)^{2}}{n}} \\ \sigma\_{\mathcal{Y}} &= \sqrt{\frac{\sum\_{i=1}^{n} \left(\tilde{\mathbf{x}}\_{i} \sin \theta - \tilde{\mathbf{y}}\_{i} \cos \theta\right)^{2}}{n}} \end{aligned} \tag{5}$$

In the formula, *σ<sup>x</sup>* is the length of the long axis of the ellipse, and *σ<sup>y</sup>* is the length of the short axis of the ellipse [20].

#### 2.3.4. Analysis of Geographical Factors


curves of sites with different slopes in three periods to explore the impact of slopes on the distribution of early sites.

