Phenological Observations on Classical Prehistoric Sites in the Middle and Lower Reaches of the Yellow River Based on Landsat NDVI Time Series

Abstract

Buried archeological features show up as crop marks that are mostly visible using high-resolution image data. Such data are costly and restricted to small regions and time domains. However, a time series of freely available medium resolution imagery can be employed to detect crop growth changes to reveal subtle surface marks in large areas. This paper aims to study the classical Chinese settlements of Taosi and Erlitou over large areas using Landsat NDVI time series crop phenology to determine the optimum periods for detection and monitoring of crop anomalies. Burial areas (such as the palace area and the sacrificial area) were selected as the research area while the surrounding empty fields with a low density of ancient features were used as reference regions. Landsat NDVI covering two years’ growth periods of wheat and maize were computed and analyzed using Pearson’s correlation coefficient and Euclidean distance. Similarities or disparities between the burial areas and their empty areas were computed using the Hausdorff distance. Based on the phenology of crop growth, the time series NDVI images of winter wheat and summer maize were generated to analyze crop anomalies in the archeological sites. Results show that the Hausdorff distance was high during the critical stages of water for both crops and that the images of high Hausdorff distance can provide more obvious subsurface archeological information.

1. Introduction

Remote sensing images offer an effective approach for archeological/historical site detection and monitoring by providing dynamic information about the landscape and environment. High spatial resolution and hyperspectral optical images acquired by multiple sensors onboard aerial aircrafts [1], satellites [2] and unmanned aerial vehicles (UAVs) [3] have been widely applied to detect crop and soil marks of ancient sites and remnants above or below the surface [4,5] in numerous archeological investigations worldwide. Ancient human activities have greatly impacted the soil characteristics by altering the natural soil forming processes. The long and strong presence of buried ruins has also had a prolonged interference on the soil material and its properties over time [6]. Features such as the remnants of temples and palaces built by stones and rammed earth, lineal features comprising of walls, paleochannels and ancient roads, emerging or buried, affect crop growth and are displayed as crop anomalies that can be interpreted explicitly from detailed remote sensing images, depending on their bands and acquired time.

Previous research has shown that spectral regions of around 700–800 nm (near infrared and red bands) are suitable for detecting crop marks [7]. The Archaeology Index—computed from narrow band imagery—can enhance crop marks related to buried archeological features [8]. For wide band optical satellite images (such as ASTER and Landsat), Normalized Difference Vegetation Index (NDVI) computed from near infrared and red bands has proven to be an effective way to extract crop anomalies by using a relative spectral response (RSR) filter and linear coefficients for different sensors to derive linear equations in combination with VNIR (Visible and Near-Infrared Region) bands to get the component images of crop marks, vegetation and soil [9].

Crop marks or crop anomalies are marks on the ground surface that have obvious or subtle differences from their surrounding areas of crop growth caused by near-surface archeological features that either enhance (positive) or inhibit (negative) crop growth. These differences can be detected by remote sensing images [10]. The appearance of the crop anomalies and crop marks is accounted for by a combination of different factors and parameters (crop, soil depth, etc.). A difference in soil moisture of the features is considered a key factor as it may enhance or reduce the growth of crops overlying archeological features [11]. Some research focusing on the appearance time of crop marks shows that archeological crop marks in England or the Mediterranean region are more visible at certain seasons every year [12,13]. Unlike high-resolution image data, Landsat imagery has not only the advantage of low cost, but also periodicity, which is useful in temporal research based on multi-temporal data (such as phenological [14,15] and change monitoring [16]). NDVI series images derived from a Landsat time series have been used to monitor the phenological growth patterns for crops, so it can determine the optimum temporal window for crop mark detection [17].

There are many prehistoric settlement sites with splendid cultures well known for their importance in the origins of Chinese civilization, such as Liangzhu, Taosi, Shijiahe, and Erlitou [18]. Although a number of them are believed to be early cities [19], only a few of them have been studied as part of a regional context by considering their unique attributes such as large scale, clear boundaries, few modern structures and an identical cropping system. In this paper, the Yellow River valley was selected with the aim of determining the best growth stages of different crops for shallow burial archeological site detection and monitoring. Taosi [20] and Erlitou [21] were selected as research sites based on extensive fieldwork, historical documents, and remote sensing interpretation, as they were both capitals at the origin of Chinese civilization in the middle and lower reaches of the Yellow River Basin.

This study focuses on shallow burial archeological sites that are less than 200 cm below the surface and the phenology of the crops above the ground. Key areas (palace area, sacrifice area, etc.) that represent classic Chinese prehistoric sites were selected as the study area while their surrounding empty fields with low ancient features density were used as reference regions. Based on the acquisition time of the effective Landsat data and the phenological changes of the crops at the site, the study focused on the critical stages of crop growth with little human impact. For every stage, NDVI was computed and analyzed using the Hausdorff distance to measure the similarity between the burial and empty areas. Based on the theory of set, a high Hausdorff distance value shows less similarity between two sets, therefore leading to a high possibility of exposing abnormal crop growth. A Landsat TM/ETM+ (Thematic Mapper/Enhanced Thematic Mapper) NDVI image series was used to check the interaction between the crops and archeological features. Results show that the Hausdorff distance was high during a critical stage of water for wheat and maize. The Landsat time series has confirmed the presence of crop anomalies arising from buried rammed earth relics within the period with a higher distance. The results indicate that this method can provide useful details about crop growth so that the optimum time window for crop anomaly identification can be obtained. This work can be used as a pre-study for high-cost remote sensing data for archeological applications and field archeological research.

2. Materials and Methods

2.1. Study Area

Taosi is situated on the piedmont alluvial plains of the Linfen Basin at the foot of Ta’ershan Hill, to the east of the Fenhe River, Xiangfen County of Shanxi Linfen (Figure 1b). Currently, it is considered as the biggest settlement with walls and a political center of the Late Neolithic Culture (2500–1900 BC) in the Linfen basin [20]. Research on the Loess-Paleosol Sequence at the Taosi section shows that the paleosol occurs to a depth of about 40–204 cm below the surface and pottery sherds from the Late Taosi Age (2050–1850 BC) accumulated between 60 and 124 cm below the surface [22]. Records show that there are cultural layers that contain remains from the Ming-Qing Period (1368–1840) or Tang-Song Period (618–1279) above the Longshan cultural layer, which creates a 10–15 cm thick hard stratum 10–15 cm below the surface [23]. Continuous walls, concentrated palaces, tombs and cellar holes that were used for underground storage were also found at the Taosi site.

The Erlitou site is considered as the political center of the Yiluo River catchment during the formative period of early Chinese states [24]. It is located to the south of the Luo River, on a triangular upland in Yanshi County, Henan Luoyang (Figure 1c). The lower reaches of both the Yi River and Luo River run through the floodplains and unite 10 km southeast of the Erlitou site. There is sufficient evidence to show that Erlitou was seen as the first capital site and developed into a large political and economic center around 1800 BC [24] in China, while most prehistoric agricultural culture declined around 2000 BC. Physicochemical analysis of the soil profile indicated that the paleosol has cultural and heritage features at depths of 45–65 cm and 105–150 cm underground [6]. Excavation records show that the cultural strata were buried 30–50 cm deep and included treasures of the Erlitou Age (1750–1500 BC) and Erligang of the Shang Dynasty (1600–1300 BC) [25]. A large number of palaces and tombs were also found at the Erlitou site, as was a cast copper region, used for casting bronze.

Both sites have a continental monsoon climate characterized by concentrated rainfall from July–September and a dry season from March–May. Figure 1b,e are high spatial resolution true color images and Corona panchromatic imagery of the Taosi and Erlitou sites. Both sites are located in a rural area beside some villages and the landscape shows no considerable change, except for the development of larger villages from the 1960s–2000s. The consistency of the crop planting patterns and varieties at these sites have provided the basis of this research. Until 2012, the two sites have not been excavated in large scale, but have experienced conservation tillage after the archeological features were discovered. The irrigation practices of both sites are primitive, and mainly rely on a natural water supply. The Taosi site is located at a higher altitude and has a low groundwater level, so relies on rainwater and surface runoff from precipitation and terrain. The Erlitou site is located on the south bank of Luohe, closer to the river, so the water supply for crop growth is reliant on river irrigation and precipitation.

Figure 2 shows the NDVI average curves derived from MODIS (Moderate Resolution Imaging Spectroradiometer) from 2000–2010 for each sample group, indicating planting time and the growth stage. The NDVI data has been smoothed using the Savitzky–Golay filter to reduce noise. Previous research shows that the results can provide a better fit, particularly in pixels where multicropping was present [26]. Winter wheat and summer maize both have fibrous root systems, and the root depth for winter wheat is generally 200 cm around loess regions [27] while the maximum root depth for summer maize is 230 cm [28]. Before the ripening period of summer maize, 80% of the root weight is concentrated within 40 cm of soil underground [29,30] while the filling stage of winter wheat has 90% of the root weight concentrated within 80 cm of soil underground [31]. It is worth noting that in the Taosi group, the summer maize growth was stunted in the research area. The site lies on the piedmont plain of the Ta’ershan Hill where groundwater exists at a depth of between 80 to 120 m [32]. Statistics show there were crop failures in most of these areas during the dry summer of 1997–2002 and 2004–2010. As above-mentioned, the water supply for crops in these areas is largely reliant on rainfall and surface runoff on the ground, which caused the poor growth of summer maize at the Taosi site.

2.2. Data and Samples

2.2.1. Data

A Landsat TM/ETM+ NDVI time series was used in this study. The time distribution of Landsat data across different periods is not uniform due to cloud cover on the images that led to an unequal sampling interval of available Landsat data. Periods of crop changes arising from human farming factors (sowing and seeding, ripening and harvest) that affected the NDVI the most were removed immediately. According to the Landsat data acquired time and the change trend of NDVI data in each period, the study period focused on critical stages of crop growth, and similarity and distance measurements were carried out for data sampling across each stage. For consistency, data were collected following the winter and maize seasons. For Taosi, the study time was from October 2000 (winter wheat sowing) to September 2002 (summer corn end of harvest); for Erlitou, it was from October 1999 (winter wheat planting season) to September 2001 (summer corn end of harvest). Data were collected for two seasons for each site.

Table 1. Landsat data acquisition time.

Site	Crops	Stages (Start Date-End Date)	Acquired date of Landsat
Taosi	Winter wheat	Tillering and Over-wintering (15 November 2000–15 February 2001 and 15 November 2001–15 February 2002)	22 November 2000, 16 December 2000, 24 December 2000, 31 December 2000, 16 January 2001, 17 January 2001, 2 February 2001, 9 February 2001, 10 February 2001, 17 February 2001, 16 November 2001, 17 November 2001, 3 January 2002, 4 January 2002, 20 January 2002, 4 February 2002, 5 February 2002, 11 February 2002.
		Jointing (5 March 2001–4 April 2001 and 5 March 2002–4 April 2002)	13 March 2001, 21 March 2001, 29 March 2001, 30 March 2001, 8 March 2002, 9 March 2002, 17 March 2002, 24 March 2002, 2 April 2002.
		Heading and Flowering (5 April 2001–10 May 2001 and 5 April 2002–10 May 2002)	14 April 2001, 1 May 2001, 8 May 2001, 19 May 2001, 18 April 2002.
	Summer maize	Jointing and Booting (18 June 2001–24 July 2001 and 18 June 2002–24 July 2002)	18 June 2001, 25 June 2001, 4 July 2001, 12 July 2001, 29 June 2002, 15 July 2002.
		Heading and Flowering (25 July 2001–23 August 2001 and 25 July 2002–23 August 2002)	4 August 2001, 5 August 2001, 12 August 2001, 13 August 2001, 30 July 2002.
		Filling (24 August 2001–15 September 2001 and 24 August 2002–15 September 2002)	6 September 2001, 13 September 2001, 24 August 2002, 31 August 2002, 1 September 2002, 3 September 2002, 9 September 2002.
Erlitou	Winter wheat	Three-leaf and Tillering (1 November 1999–4 December 1999 and 1 November 2000–4 December 2000)	12 November 1999, 20 November 1999, 28 November 1999, 29 November 1999, 22 November 2000.
		Over-wintering (5 December 1999–10 February 2000 and 5 December 2000–10 February 2001)	7 December 1999, 14 December 1999, 22 December 1999, 30 December 1999, 24 December 2000, 10 January 2001, 26 January 2001.
		Jointing (10 March 2000–10 April 2000 and 10 March 2001–10 April 2001)	12 March 2000, 19 March 2000, 27 March 2000, 28 March 2000, 4 April 2000, 5 April 2000, 30 March 2001, 31 March 2001, 7 April 2001.
		Heading and Flowering (11 April 2000–11 May 2000 and 11 April 2001–11 May 2001)	12 April 2000, 13 April 2000, 28 April 2000, 6 May 2000, 7 May 2000, 16 April 2001, 24 April 2001, 1 May 2001, 2 May 2001, 10 May 2001.
	Summer maize	Jointing and Booting (15 June 2000–30 July 2000 and 15 June 2001–30 July 2001)	15 June 2000, 16 June 2000, 23 June 2000, 4 July 2000, 9 July 2000, 17 July 2000, 18 June 2001, 27 June 2001, 4 July 2001, 12 July 2001.
		Heading to Filling (31 July 2000–10 September 2000 and 31 July 2001–10 September 2001)	10 August 2000, 19 August 2000, 27 August 2000, 6 August 2001, 21 August 2001, 22 August 2001, 6 September 2001.

shows the acquisition time of the data used. Radiometric and geometric corrections were applied before the NDVI computations. The radiation corrections adopted in this paper include radiation calibration and QUAC (Quick Atmospheric Correction). QUAC is a quick atmospheric correction method of ENVI (The Environment for Visualizing Images) for hyperspectral and multispectral images, and automatically collects spectral information of different substances from images and obtains empirical values.

2.2.2. Samples

Buried remains affect surface crops in different ways based on their densities and properties. Key zones such as a palace or sacrificial area at the site were vectorized as per the site map diagrams and high-resolution images. Empty areas near the sites (with less cultural relics and cultural layers) were selected as reference samples based on archived, archeological documents and past field excavation records for comparative purposes. Pixel values from the NDVI images were extracted as phenological samples based on its acquisition date. Figure 3 shows the locations and attributes of the samples on the high spatial resolution images and

Table 2. The main attribute of samples.

Site	Area (km2)	Rotation Crops	Sample	Area (m2)	Region Attribute	Number of Pixels
Taosi	2.8	Winter wheat/Summer maize	S1	71,018	Sacrificial area	76
			S2	55,010	Palace area	61
			R1	19,000	Blank area	24
			R2	59,461	Blank area	64
Erlitou	2.9	Winter wheat/Summer maize	S	116,900	Palace area	128
			R	115,781	Blank area	137

represents the attributes and the pixel numbers (N) of each sample for each site. Bare areas were removed and the main crop areas were extracted. Both sites had the same crop species and planting modes.

2.3. Method and Application

2.3.1. Similarity Measure of Temporal Series

Time series similarity/dissimilarity is a key research issue in the data mining community [33]. Distance measure is one significant approach to time series clustering and classification, despite the numerous measurements that can be applied to measure the difference among time series. Time series distance can be measured in time (Euclidean distance or correlation based distances); in shape (Dynamic Time Warping (DTW) and Longest Common Sub-Sequence (LCSS) etc.) and in change (structural) (such as HMM (Hidden Markov Model) based distance) for comparison and categorization [34]. All of the above-mentioned methods are generally employed to measure the distance between two time-series. For measuring the similarity of two sets of time series, the Hausdorff distance [35] and the modified Hausdorff (MODH) distance are designed and employed in the clustering and classification of time-series data. In this paper, the Hausdorff distance was employed to measure the similarity between two sample groups.

2.3.2. Hausdorff Distance

Given a metric space $(\mathrm{M},\mathsf{\delta })$, with metric $\mathsf{\delta }$, the Hausdorff distance [36] between two subsets $\mathrm{A},\mathrm{B}\subseteq \mathrm{M}$ is defined as:

$$d_H(A,B)=\max \{d_h(A,B),d_h(B,A)\}$$

(1)

where $d_h(A,B)=\underset{a\in A}{\max }\left(\underset{b\in B}{\min }a-b\right)=\underset{a\in A}{\sup }\underset{b\in B}{\inf }\delta (a,b)$, which measures the largest out of all distances with a ∈ A to set B. This function is not symmetric, $d_h(A,B)\ne d_h(B,A)$.

So the Hausdorff distance,

$$d_H(A,B)=\max \{\underset{a\in A}{\sup }\underset{b\in B}{\inf }(a,b),\underset{b\in B}{\sup }\underset{a\in A}{\inf }(a,b)\}$$

(2)

is clearly symmetric. If the data set is finite and consists of N elements, all distances can be arranged in an N × N matrix δ_ij and Equation (2) reads

$$d_H(A,B)=\max \{\underset{i\in A}{\max }\underset{j\in B}{\min }{\delta }_{ij},\underset{j\in B}{\max }\underset{i\in A}{\min }{\delta }_{ij}\}.$$

(3)

Based on the rigorous definition of the mathematical concept of distance, the Hausdorff distance is an efficient method for the analysis of complex sets with complicated, irregular and fractal structures and stands out more than the other methods. The Hausdorff distance measures the biggest non-matching extent of point sets and can detect the abnormal point of two sets, but yields error as noise.

In this study, we only needed to study the unidirectional Hausdorff distance $d_h(S_{burial},S_{empty})$, which is the max-minimal distance from the burial area to the empty area used to detect crop anomalies in the burial area. Most of the methods mentioned in this paper are used for similarity measurements between two individual time series; however, in this paper similarity is measured between pairs of time series sets. According to the Hausdorff distance definition, the unidirectional Hausdorff distance is the maximum of the minimum distance from one set to another. The purpose of this paper is to find the abnormal growth present in one of the sets by comparing the two sets. The Hausdorff distance is consistent with the purpose of this paper, so the Hausdorff distance was employed.

2.3.3. Application

There are many distance methods that can be used to quantify the similarity or difference between two consecutive time series. Euclidean distance and Pearson’s correlation coefficient related distance [37] were adopted to measure the difference between two time series extracted from two pixels.

The Euclidean distance was obtained by

$$d_{xy}=\sqrt{{\displaystyle \sum _{k=1}{\left(x_k-y_k\right)}^2}}$$

(4)

where x_k and y_k are the observations of samples x and y. The Euclidean distance was calculated directly and cannot be used to compare the phenological periods due to the different number of observations in each phenological period. For comparison, after the Euclidean distance was calculated, all the results in a period were linearly normalized and rescaled to the range of 0–1 and then the Hausdorff distance used for calculation.

The Pearson’s correlation coefficient was defined as a rescaled variance given by

$$d_{xy}=\sqrt{2(1-c_{xy})}$$

(5)

where c_xy is the Pearson’s correlation coefficient of two pixels computed over the investigated time period and is obtained by

$$c_{xy}=\frac{{\displaystyle \sum _{k=1}(x_k-\overline{x})(y_k-\overline{y})}}{\sqrt{{\displaystyle \sum _{k=1}{(x_k-\overline{x})}^2{\displaystyle \sum _{k=1}{(y_k-\overline{y})}^2}}}}$$

(6)

where $\overline{x}$ and $\overline{y}$ are the average of the samples x and y, the distance [38] is bounded by the interval from 0–2, as the two series changes from perfectly correlated $(c_{ij}=+1)$ to inversely correlated $(c_{ij}=-1)$, the distance decreased monotonically with correlation. In this paper, to study the obvious seasonal growth differences of the same cropping system, it was necessary to find the stages of two time series that had weak correlation $(0.3\le c_{ij}<0.6)$, no-correlation $(0<c_{ij}\le 0.3)$ or negative correlation $(c_{ij}<0)$. When the correlation is $(c_{ij}<0.6)$, the distance $d_{ij}>0.89$, so a distance threshold $d_{ij}>0.9$ was set up to detect the different time of the two series. The Hausdorff distance inherits the concept of similarity measures from the distance defined at the level of the individual time series and applies it to the level of time series sets. It is more similar among the series in one set than those belonging to other sets as the max-minimum distance within a set is smaller.

For every stage the Euclidean distance, Pearson Correlation Coefficient and its related distances from a given pixel to another pixel in the same stage was calculated. On the basis of the above-mentioned, unidirectional Hausdorff distance from a burial area to an empty area was calculated, and the stages with higher Hausdorff distance indicates periods of more obvious difference of crops growth. There is much relevant literature for the archeological study sites following extensive excavation and research over the years. Analysis was thus constrained on the archeological sites due to lack of sufficient strata records for the surrounding fields. The Landsat time series NDVI images were displayed to further examine the result in combination with the excavation reports [39,40,41]. Crop growth status was compared with existing archeological field reports on detecting archeological features.

3. Results

3.1. Taosi

Table 3. Distance measurements result for Taosi site.

Crops	Stages	Euclidean Distance		Pearson’s Correlation Distance
		S1	S2	S1	S2
Winter wheat	Tillering and Over-wintering	0.341	0.337	0.408	0.402
	Jointing	0.645	0.274	0.798	0.339
	Heading and Flowering	0.569	0.262	1.130	0.520
Summer maize	Jointing and Booting	0.386	0.269	0.707	0.492
	Heading and Flowering	0.339	0.150	0.501	0.221
	Filling	0.282	0.386	0.504	0.691

shows the distance measurement results for the Taosi site. Of the Euclidean distance results, 83% were smaller than 0.4. The calculation results for a palace area were less than 0.4 with an average of 0.28, indicating that the samples with a significant growth difference between the palace area and empty area was not indicated for each phenological stage. The Hausdorff distances derived from the Pearson’s Correlation distance were significantly higher in the jointing stage and heading and flowering stages of winter wheat than that of the other stages in the sacrificial areas. Both the Hausdorff distances derived from the Euclidean distance and Pearson’s correlation of the sacrificial areas were higher than the palace area except for the filling period of summer maize. The Pearson’s correlation distances were large; however, most of these values were below 0.9 except for the heading and flowering stage of winter wheat, implying a strong correlation between the samples. For the palace area samples, the results in the heading and flowering stage of winter wheat, the jointing and booting period, and the filling period of summer maize was relatively high. The results from the jointing period to the heading and flowering stage of winter wheat in the sacrificial area were significantly higher than those in other periods, and was the same as the Euclidean distance. At the same time, the results from the sacrificial areas were different from the Euclidean distance in the summer maize season, which were relatively high in the jointing stage of summer maize.

Figure 4 and Figure 5 are the Landsat NDVI time series images of 2000–2001 and 2001–2002. The NDVI of this area shows crop anomalies from mid-March (the late jointing stage) to mid-April, and the heading and flowering stages. Prior to this period in the tillering stage, the NDVI revealed no anomalies, but during this period the anomalies caused by archeological features were clearly identified. In the southeast, a small town in a sacrificial area and an ancient astronomical rammed earth observatory building, which was excavated during the 2003–2005 period, show a distinct negative mark. This was the same for the north sacrificial area outside the walled site and the dense storage area to the middle east of the site. The palace area where a large number of rammed earth building foundations were found from 2002–2007 could also be vaguely identified (Figure 8). The last three images show that the NDVI of the whole site was approximately equal to that of the bare soil, except at the gullies and the southeast small town with the ancient astronomical observatory. It is consistent with the unidirectional Hausdorff distance result that in the last images, which were acquired at the filling stages of the summer maize, the palace area also had positive crop anomalies.

3.2. Erlitou

Table 4. Distance measurement results for Erlitou site.

Crops	Stages	Euclidean Distance	Pearson’s Correlation Distance
Winter wheat	Three-leaf and Tillering	0.606	0.603
	Over-wintering	0.481	0.393
	Jointing	0.587	0.571
	Heading and Flowering	0.596	0.540
Summer maize	Jointing and Booting	0.513	0.389
	Heading and Flowering	0.606	1.000

shows the distance measurement results from the Erlitou site. The change of Hausdorff distance derived from Euclidean distance was not large and ranges between 0.48–0.61. However, it can be seen that the over-wintering period of winter wheat and the jointing period of summer maize were significantly lower (smaller than 0.52) than at other periods. Pearson’s correlation distance results also showed the same trend. The heading and flowering stages of summer maize were the only values greater than 0.9, indicating a strong correlation between the growth trends of the samples during most phenophases. The results from the three-leaf and tillering stage, the jointing stage, the heading and flowering stage were all relatively high.

Based on the Landsat NDVI images series from 1999–2000 and 2000–2001 (Figure 6 and Figure 7), from the jointing stage to the flowering stage of winter wheat, palace base No. 1, the southern ditch of cast copper relics region, the sacrificial area, the ancient road and the eastern wall of the imperial palace and palace base No. 5 displayed crop anomalies with remarkably low values. This buried building foundation of rammed earth was confirmed in 1974, 1984, 2003 and 2010. At the heading and flowering stages of the summer maize planting time, the eastern wall of the imperial palace, and palace bases No. 5 and No. 3 can be identified due to their low values. The archeological features stood out gradually as the root length increased. During August to early September, the archeological features became more distinct. It should be noted that at the three-leaf stage in 2000, the crop overlying the eastern wall of the imperial palace shows better vigor compared to its bare surroundings.

4. Discussion

4.1. Negative Crop Anomalies

Winter wheat and summer maize are the main grain crops grown in the middle and lower reaches of the Yellow River where the Taosi and Erlitou archeological sites are located. There are clear crop anomalies during the jointing to flowering stages of the winter wheat planting time for both sites and the heading to filling stages of the summer maize planting time for the Erlitou site. These stages were the critical and peak periods of water for winter wheat and summer maize. Furthermore, the length of root was deep enough to surpass buried archeological features. During the winter wheat planting time, the NDVI values in both groups rose sharply in the fast growth stages and large volumes of water were absorbed for growth. During the heading and flowering stages, the NDVI value stopped rising and dropped from its peak. Additionally, when the winter wheat and summer maize are transforming from vegetative growth to reproductive growth, they are very sensitive to water changes during this period. Figure 2 shows that for both sites, the NDVI curves rose to the first peak at the jointing stage and then began to fall. For the Erlitou site, the NDVI curve began to drop from the second peak at the heading and flowering stages.

The archeological layers contain remains such as palace building bases with rammed earth that restrained the growth of wheat and maize, leading to lower NDVI values and an emergence of crop anomalies. At the Taosi site, the southern small walled site and ancient astronomical observatory had lower changes of NDVI due to the concentrated distribution of buried features and water stress on the roots impairing growth of the winter wheat. The palace area, the storage regions with dense cellar holes and the northern sacrificial area with walls outside the walled town also displayed vaguely in the images during the jointing to flowering stages of winter wheat. Similarly, the NDVI performed well in the identification of archeological layers of the Erlitou site in the jointing and heading stages of the winter wheat, and the heading to filling stages of the summer maize planting season, especially in regions such as the palace area where archeological features were concentrated.

Figure 8 and Figure 9 show the distribution of negative crop anomalies and the main archeological features of Taosi and Erlitou. Obvious negative crop anomalies are seen in the Landsat NDVI images, and detailed distribution about their relative archeological features is shown in the high spatial resolution images. It should be noted that in Figure 8, the NDVI values were affected by the dry gullies at the Taosi site, which are lower than the other parts from the jointing to flowering stages. The Taosi site is located at the foot of the Ta’ershan Hill with an altitude of 500–600 m, a low under groundwater level, and is classified as a continental monsoon climate zone where most of the phenological window was within the dry season from March–May. During the dry season, intense evaporation occurs leading to high soil salinization due to salt accumulation. This resulted in the restrained growth of winter wheat, as indicated by the low values in gully areas where the groundwater level is high.

4.2. Positive Crop Anomalies

Figure 10 shows the distribution of positive crop anomalies and the main archeological features at the Taosi site, where significant positive crop anomalies were seen at the palace area and the southern small town on the Landsat NDVI image. The Google Earth image shows the distribution of the archeological features and crop growth in this region during the planting of the summer maize. Crop anomalies on the Landsat image series of the winter wheat and summer maize show the apparent differences between the gullies and the burial area. Contrary to the winter wheat rapid growth stages, the Landsat TM/ETM+ data illustrated that the summer maize reduced across the whole Taosi site, except at the gullies and the southern small walled site. The NDVI values at the dried gullies showed higher values than at other parts of the site. The southeast small town of worship and sacrifice with the ancient astronomical observatory had remarkably high values. On the images of the summer maize filling stage, the palace area also showed high values, as it did at Erlitou, where the three-leaf stage of winter wheat and other effective phenological stages also showed the opposite crop marks. Figure 11 shows the NDVI image series of the early stages (emergence, three-leaf, and tillering) of the winter wheat planting period in 1999 and 2000 at the Erlitou site. By monitoring crop anomalies, it was found that the growth of abnormal crops was better than that of the surrounding areas at the emergence and three-leaf stages of winter wheat, but that the anomaly gradually disappeared after the tillering stage of winter wheat.

The Taosi site is located along the hill at an altitude of 500–600 m and its low groundwater level contributed to the unusual growth of summer maize. During the summer maize planting time, the dry weather caused the summer crops to perform poorly at most parts of this site except in the gullies, where the groundwater level and soil moisture are high. The buried features slowed down the groundwater flow and acted as water storage that sustained the growth of summer crops unlike other areas across the site. The Erlitou site has a typical continental monsoon climate with rainfall concentrated in the months of July–September. The rainfall in September and October affects soil moisture, which directly affects the sprouts and seedlings of winter wheat. According to the excavation records, the existence of the ancient road stopped the downward flow of groundwater, and is where positive crop anomalies were revealed in the early stage images. This resulted from soil moisture accumulating between the shallow buried features and the surface being higher than the surrounding area. Therefore, during the early winter wheat planting period, some of the crops above the shallow buried features displayed positive marks.

Positive crop anomalies at both sites first appeared at the emergence stage when the budding rate and growth rate are heavily dependent on soil moisture. Buried features with dense texture hinders the downward infiltration of groundwater, resulting in soil moisture between the buried features and the higher topsoil than the surrounding area, which leads to better crop growth than the surrounding area. This phenomenon continued until the three-leaf stage. For the winter wheat grown at the Erlitou site, up to the tillering period, the anomalies gradually disappeared and the crop growing trend within the entire site became more consistent. During the over-wintering period, the growth of winter wheat was homogeneous, thus there were no significant differences. However, for the summer maize at the Taosi site, due to the planting failure of the surrounding crops, the positive anomaly continued throughout the summer maize season. It was evident that, although the image series showed effective results, image quality was affected by weather changes. It also illustrated that for the same crop rotation system, soil aridity and salinity impacted the interaction of the buried features and the crops, which displayed a different performance of crop anomalies.

4.3. Distances

Two Hausdorff distance measurements derived from the Euclidean and Pearson correlation were compared with the Landsat image series and it was found that the results from the Euclidean distance were not significant for the positive crop anomalies at the Taosi site, but significant for the negative crop anomalies. The optimum time window of the crop anomalies during the summer maize could not be identified. The results from the Pearson correlation coefficient were more sensitive to the presence of both anomalies. At the Erlitou site, for both positive and negative crop anomalies, the results of the Euclidean distance and Pearson’s correlation distance measurements were significant and could identify the periods when crop anomalies appeared.

One reason is that the Taosi site—due to the failure of the summer maize—the overall NDVI value of the samples was low, hence the correlation coefficient was more sensitive. Even in areas where positive vegetation markers were present, crop growth was inhibited by drought during the early period and the buried features during the late period with NDVI values ranging from 0–0.3 for all samples during this time. Another reason is that of comparison. After calculating the Euclidean distance between every two samples in a period, the whole Euclidean distance was normalized in the range of 0–1 before the Hausdorff distance was calculated. However, during this process, the size of the Hausdorff distance was affected by the full range of the sample distance, particularly the maximum value of the sample distances.

The statistical evidence showed that the ratio of images with crop anomalies was higher in the periods of high Hausdorff distance. In the study time, 100% of the NDVI images during the jointing period of winter wheat, the jointing and booting stage, the heading and flowering stage, and filling stage of summer maize at the Taosi site and the jointing stage, heading and flowering stage of winter wheat, and the heading and flowering stage of summer maize of Erlitou site showed crop anomalies. Similarly, 80% of the NDVI images of the heading and flowering period of the winter wheat at the Taosi site and 60% of the NDVI images of the three-leaf and tillering period of winter wheat at the Erlitou site also showed crop anomalies. However, in the period of lower Hausdorff distance, only 22% of the NDVI images of the tillering and over-wintering of winter wheat at the Taosi site, and 14% of the NDVI images of the over-wintering of winter wheat and 30% of the NDVI images of the jointing and booting of summer maize at the Erlitou site showed crop anomalies. Furthermore, weather changes that caused different soil moisture content also affected the emergence and area of crop anomalies for the same crop. In some cases, crop anomalies appeared in images of some years within the stages of low Hausdorff distances and no anomalies occurred in other years of this period, while within the stages of high Hausdorff distance, most of the images displayed crops anomalies. Therefore, further study of the data sampling within the phenological stages should be as uniform as possible with the sampling year range extended. In this study, Landsat NDVI time series was used to observe the phenology of crops in the middle and lower reaches of the Yellow River. The jointing and heading stages of winter wheat can be used as the optimum time windows for data selection for remote sensing archeological applications of the middle and lower reaches of the Yellow River. During this time, the weather had less impact on the appearance of crop anomalies and crop anomalies were more clearly identified.

5. Conclusions

In this paper, two classic prehistoric sites, Taosi and Erlitou in the middle and the lower reaches of the Yellow River were studied. The Hausdorff distance method, which is based on the concept of set, was used to measure the similarity of the archeological sites between mainly burial areas and empty areas. The crop anomalies and their phenological time windows were determined by comparing the sampling area with the empty area.

The Hausdorff distance was high between the burial area and their reference areas during the critical stages of water consumption: the three-leaf stages, from the jointing to flowering stages of winter wheat for the Taosi and Erlitou sites, and the heading to filling stages of summer maize for the Erlitou site. During these seasons, large growth differences appeared, which is evidenced by the presence of crop anomalies, thus indicating the possibility of buried relics in such places. Results also show the consistency between the Hausdorff distance and the Landsat image series. In the study period, when the Hausdorff distance calculated using the Pearson correlation coefficient distance was between 0.5–1.2, 60–100% of the images showed crop anomalies, and only 14–30% of the NDVI images at a Hausdorff distance of 0.2–0.5 showed crop anomalies. The Landsat images series also show that the crop marks or anomalies were more distinct at stages of high Hausdorff distances and vague at stages of lower distances.

The Taosi and Erlitou studies showed the effect of archeological features on the growth of crops. Crop marks and anomalies showed in the denser region or large areas of buried rammed earth relics (palace area, sacrificial area). Each site has unique complexities including unique climate and land conditions that have different impacts on the surface crop and therefore affected detection results. Weather changes can disrupt the normal growth of crops as well as affect the quality of the remote sensing images, thus should be considered in the selection of study period for remote sensing archeological applications. The phenology played an important role in the study of the stress impacts that the buried archeological features have on crop growth. In this study, the Pearson correlation coefficient was more sensitive to the presence of crop anomalies while the Euclidean distances were relatively insensitive to low NDVI values. Furthermore, the Euclidean distance posed some difficulties due to differences in the number of observations per period.

In practice, areas with high density of archeological features are considered as study samples while the adjacent empty areas are used as reference samples. The methods mentioned in this paper are used to detect crop anomalies and their phenological time windows. For the large area sites, an appropriate size of the sampling window can be set at the study area for spatial sampling. Finally, each sampling area is compared with the reference region to determine crop anomalies and their phenological time windows. The above-mentioned method can also be used to determine the phenological window period of crop anomalies when the burial area is determined and the empty area is uncertain. This method can be used as a preparatory study for large-area archeological prospection before deploying more accurate remote sensing data and other ground methods. As it is difficult to conduct large-area and continuous timing studies for radar and high-resolution image data, the area and time range of the images must be very precise. This method can also be used to determine the area and time of the required data. Additionally, the Landsat satellite data in this paper has the characteristics of massive archiving and traceability, so the method can also be used in systematic remote sensing archeological investigation.