Remote sensing detection and verification of disappeared reservoirs along the Grand Canal of China

Abstract

Digital Earth is an information-rich expression of the real Earth and is a new way of understanding the Earth in the twenty-first century. Archeology has found great potentialities in Digital Earth, strongly increasing its development and its interdisciplinary experimentations. With the rapid development of remote sensing, the importance of applying Earth observation technology to archeological research has caught a great amount of attention worldwide. The objective of this study is to apply Radarsat-1 and Landsat imagery, historical maps, and aerial photographs in the detection and verification of disappeared reservoirs built 500 years ago as part of the Grand Canal of China. The results of the study show that the potential sites present distinct variations in soil moisture, biomass, and biological vigor, characterized by weak backscattering in Radarsat imagery; high vegetation indices in Landsat imagery; and distinct anomalies in a Landsat principle components image. Historical maps and aerial photographs also verified the existence of the sites of disappeared reservoirs. Since the sites have so far not been extensively explored, the results should represent a significant contribution to the understanding of the Grand Canal of China and its auxiliary facilities.

Keywords:

• remote sensing
• digital archeology
• Grand Canal
• detection
• anomaly

1. Introduction

Archeology has a long history of searching the Earth's surface for minute traces of human impact created centuries or millennia prior, in the hope of piecing together the history and development of our species (Osicki 2005). Since the beginning of the twentieth century, aerial photography has been used in archeology primarily to view features on the Earth's surface that are difficult if not impossible to visualize from the ground (De Laet et al. 2007). With the launch of the first Landsat satellite in 1972, satellite remote sensing also became accessible to the archeological community (Clark et al. 1998). The use of remotely sensed images and aerial photographs has become an essential part of archeological landscape studies (Scollar et al. 1990, Bescoby 2006). Such images often reveal spatial or connective relationships between extant archeological remains and other morphological features relating to the use and organization of the landscape in the past (Bescoby 2006).

Remote sensing technologies have provided archeologists a powerful tool set with which to further analyze past human lifeways. Satellite and aerial sensors, individually or in conjunction, have allowed for a rapid, accurate, and nondestructive means to prospect for unidentified sites, as well as to identify and analyze characteristics within and around sites (Zhao 2004). The integrated use of remote sensing technology and methods, besides having important scientific implications, can contribute significantly to the preservation, monitoring, and management of cultural resources (Masini et al. 2008).

This article presents preliminary results obtained from the application of Radarsat-1 and Landsat imagery to detect, identify, and locate disappeared reservoirs built to support the Grand Canal of China over 500 years ago, and verified potential sites using historical surveying maps and aerial photographs. The results are essential for decision-making in the Grand Canal's preservation and management. This study suggests that the complementary use of different imagery, especially combined use of optical and synthetic aperture radar (SAR) data imagery, could provide satisfactory results.

2. History and background

The Grand Canal of China is the longest ancient canal in the world. The oldest part of the canal, Han Gou, was built in 486 BC during the Spring and Autumn Period (722 to 481 BC). Completed in the Sui dynasty (581–618), the Grand Canal laid the foundation for the prosperity of the Tang and subsequent dynasties (Doar 2007). During the Ming (1368–1644) and Qing (1644–1911) dynasties, the Grand Canal served as the main artery between northern and southern China and was essential for the transport of grain to Beijing (Figure 1), the main artery of China's historical development and the symbol of great ecological civilization in China (Mao et al. 2008).

The first written account of the Grand Canal came from China's sixth-century geographer Li Dao Yuan in his ‘Shuijingzhu’ (Commentaries on the Book of Rivers, the most ancient literature on rivers in China). Around 1280, Marco Polo was appointed by Kublai Khan to the governorship of the important city of Yangzhou, which was born with the birth of the Grand Canal. He described the canal as ‘a wide and deep channel dug between stream and stream, between lake and lake, forming as it were a great river on which large vessels can ply.’ In 1937, the American natural historian Willard Price made a visit to the Grand Canal and published his article in National Geographic Magazine (Price 1937). In his article, the status and cultures of the Grand Canal were clearly shown to the western world for the first time and aroused great interest around the world. In 2008, several municipalities and provinces in China submitted an application to have the canal added to UNESCO's World Heritage List. A decision is expected by 2014.

Figure 1.  Modern course of the Grand Canal of China. The red line denotes the course of the Grand Canal of China, while the blue lines denote the course of the major natural rivers.

The elevation of the canal bed varies (Figure 2). It reaches 29.5 m at Beijing and 38.5 m at its summit at Nanwang. In order to guarantee water transportation, auxiliary projects were undertaken at Nanwang in the Ming Dynasty, which used various facilities to adjust the water flow of the canal so that transport could be controlled. During the projects, five reservoirs, namely Anshan, Nanwang, Mata, Shushan, and Machang, were dug to store sufficient water in the rainy season and to feed the Grand Canal in the dry season. Because of various factors, the canal languished for more than one and a half centuries (Deng and Guo 2009, Deng et al. 2010). Many of the canal sections fell into disrepair, and some parts were returned to flat fields, including the reservoirs.

Figure 2.  Elevations of cities along the canal. The elevations were extracted from an ASTER Global Digital Elevation Model, produced with 30 m postings, and formatted in 1×1 degree tiles as GeoTIFF files.

3. Rational basis

Archeological sites do not exhibit spectral signatures that can be used for generic detection purposes. Rather, it is hypothesized that archeological residues produce localized contrasts in the landscape matrix which can be detected using an appropriate sensor under appropriate conditions (Beck 2007). Ancient archeological remains will affect their surrounding landscapes in different ways, whether through altering the surrounding soils, affecting vegetation, or absorbing moisture at a greater or lesser rate (Parcak 2009).

The primary aim of archeological remote sensing is to find contrasts between archeological features and the natural background. When archeological deposits or features possess physical properties different from those of the surrounding matrix, a distinction may be noticed. Contrasts are referred to as ‘anomalies’ until they can be identified, a task that may require excavation (Kvamme 2005). Detection attempts must be conducted with a clear understanding of the archeological and environmental parameters of the region under investigation. This need is amplified in the investigation of disappeared reservoirs, which in general leave a subtle imprint on the landscape. Direct detection need not be the only or even first option; with available remote sensing imagery, environmental indicators connected to archeological sites are more easily identified.

3.1. SAR detection

Imaging radar is a powerful tool for archeology. This technology's unique ability to penetrate underneath desert sands and to unveil ancient topography has astonished experts in and out of the community of archeological remote sensing (El-Baz 1997). Radar transmissions have the ability to penetrate certain objects and are therefore useful where surface features are covered by dense vegetation and can, to a certain extent, detect subsurface features (Hall 2010). While reflectance measurements in the visible and infrared (IR) regions are primarily related to molecular resonances of surface materials, the backscatter of a radar wave is largely determined by the physical properties of surface objects such as surface roughness and the dielectric constant (Henderson and Xia 1997). Moisture in the landscape influences the backscattering coefficient through changes in the dielectric constant of landscape materials. The dielectric constant is a measure of the ability of a substance to conduct electrical energy – an important variable determining the response of a substance that is illuminated with microwave energy. The backscatter of a radar wave can also be influenced by the geometric configuration of targets. This behavior is caused by objects classified as corner reflectors, which often are in fact corner-shaped features (such as the corners of buildings and the alleyways between them in a dense urban landscape), but are also formed by other objects of complex shape. Corner reflectors are common in urban areas due to the abundance of concrete, masonry, and metal surfaces constructed in complex angular shapes.

When a ditch or a reservoir is dug and refilled, the local soil-profile gets disturbed. In most cases, the refilled material is different from the undisturbed soil. The differences are twofold. Firstly, the material itself is different. In most cases, Refilled material is not so compact, and it contains more humus components, making it look darker. Secondly, the grain size distribution differs from the undisturbed soil. Refillings mostly have a larger number of smaller grain sizes. Both cases lead to better water storage and a smaller soil roughness, which backscatter less radar waves than the surrounding soil and can be detected very well from radar imagery.

3.2. Landsat detection

One common source of remote sensing imagery used for archeology is Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper Plus (ETM +). The space-borne TM/ETM+ sensors are capable of returning seven/eight channels of digital data back to a receiver on Earth: three in the red, green, and blue visible spectrum; three in the reflective IR spectrum; one in the thermal IR spectrum; and one panchromatic band (for ETM+ only). Band 1 (with a central wavelength of 0.49 µm), referred to here as TM1, is used for coastal water studies; TM2 (0.56 µm) is used for crop identification and vegetation stage studies; TM3 and TM4 (0.66 and 0.83 µm, respectively) are used to calculate vegetation indexes, such as the Normalized Difference Vegetation Index (NDVI); TM5 and TM7 (1.65 and 2.22 µm, respectively) can be used for clouds, ice, snow, and geological formation discrimination; and finally TM6 (with an effective wavelength of 11.457 µm) is used for land surface temperature retrieval.

Direct detection of archeological sites such as disappeared reservoirs is implausible both because of the unknown nature of the remains, and because there is no detailed location of the remains in our study site. In the case of disappeared reservoirs, however, refilled archeological features can generally store more water, and the plants growing over them will have more water available. This will make them grow higher and ripen later than the plants around them. The plants over them will stay green for a longer period and will have a darker green, producing a tonal contrast. Then, the use of Landsat for archeological prospection in our study site becomes a search for anomalies. Discrete locations that evidence more biomass and biological vigor than their surroundings are anomalous and may represent areas of disappeared reservoirs.

4. Materials

In this study, Radarsat-1 satellite imagery was acquired for the region in 1997. The Radarsat SAR operates at the same C-band frequency of 5.3 GHz with a polarization of horizontal transmit–horizontal receive. Data were obtained in ScanSAR Wide-beam B mode, which has a swath width of 450 km over incidence angles from 20 to 47 degrees, and a 100 m pixel size in the imagery presented here.

We have used Landsat Multispectral Scanner (MSS) data from 1975, 1978, and 1979 together with Landsat TM data from 1987, and Landsat ETM+ data from 2001, 2002, and 2004. Details on satellite, sensor, track and frame, and acquisition date are given in Table 1. Global Land Cover Facility (a NASA-funded member of the Earth Science Information Partnership at the University of Maryland) provided the data in a geo-corrected format. The selection of the satellite images was based on the availability of high-quality satellite imagery and minimal cloud coverage.

Table 1. Satellite data collected in the project.

Satellite	Sensor or image mode	Track/frame	Date	Spatial resolution
Radarsat-1	ScanSAR	Mosaic	17 August 1997	100 m (resampled to 50 m)
Landsat-2	MSS	132/35	21 August 1975	80 m (resampled to 57 m)
Landsat-3	MSS	132/35	3 June 1978	80 m (resampled to 57 m)
Landsat-2	MSS	131/35	19 May 1979	80 m (resampled to 57 m)
Landsat-5	TM	122/35	14 May 1987	28.5 m
Landsat-5	TM	122/35	19 September 1987	28.5 m
Landsat-7	ETM +	122/35	12 May 2001	28.5 m
Landsat-7	ETM +	122/35	31 May 2002	28.5 m
Landsat-7	ETM +	122/35	9 September 2004	28.5 m

5. Methods and results

The SAR images were obtained in mosaic format. The SAR mosaic imagery was imported into ERDAS IMAGINE software version 8.7 for processing and analysis. The images were filtered using a 3×3 median filter to reduce noise. Different filters and window sizes from ERDAS IMAGINE radar texture classification tools were investigated and compared before making a decision on the most suitable filter and window combination. The median filter was chosen because it preserves edges while smoothing the data (Bourgeau-Chavez et al. 2005). The images were then georeferenced to a common base map. A Landsat TM image from 31 May 2002 was used in geometric correction and registration of the SAR images.

Figure 3 shows a dark area in the middle, indicating low backscattering for that area. It shows stronger backscattering than the water body and weaker backscattering than the land between the settlements. As mentioned in the previous section, this phenomenon can be attributed to the following: surface roughness, dielectric constant, and geometric configuration of targets. As for the geometric configuration of targets, this can be excluded. Surface roughness can also be excluded in that it is not likely that distinct variations of surface roughness occur in this study site. So, the only possible reason for this phenomenon is a low dielectric constant in that area compared with the surrounding areas.

Figure 3.  Subset from the Radarsat-1 mosaic imagery of the study site. The red ellipse denotes the area with a low backscattering coefficient.

In order to delineate this area with a low backscattering coefficient, the point density was calculated using ArcGIS Spatial Analyst Toolbox. To create comparable densities at different locations, one usually chooses comparable distance weighting functions around each point. Mathematically, these weights are functions K(s, t) of the vector displacement (s, t) of any point relative to the location (x, y), where the density is computed. A population (or other value) f(x+s, y+t) located at (x+s, y+t) will receive weight K(s, t) and thereby contribute a term f(x+s, y+t)*K(s, t) to the density estimate at (x, y). The density at (x, y) is the sum of all such terms (the dimensions of K() therefore are 1/area). For this to be a true density, the integral of K(s, t) over all possible (s, t) must be in accord. Spatial Analyst Toolbox supports only two kernels. One is the unweighted circular kernel K(s, t) = 1/(Pi*R ²) for values (s, t) inside the circle of radius R (and equal to zero otherwise). The other is a quartic approximation to a Gaussian kernel. The Gaussian kernel is K(s, t) = exp(−(s ²+t ²)/(2C ²))/(2*Pi*C ²). The value of C determines the effective range of K; points (s, t) beyond a distance of about 3*C barely contribute to K. K() is essentially flat beyond 3*C. The process for density mapping of the study site is as follows. First, the boundaries of the settlements and towns were extracted by density slicing and morphological filtering in ENVI 4.7, and then the file type of resulting boundaries was converted from evf to shp to import to ArcGIS. Then a cluster analysis was conducted to the boundaries of the settlements and towns to refine the results, and the point density was calculated using ArcGIS Spatial Analyst Toolbox with an unweighted circular kernel. And finally, the point density value was classified to five intervals and mapped, as shown in Figure 4. From this point density map, it is easy to see that an area with low density lies in the middle of the map. This area with low point density can depict the soil marks shown in Figure 4, in that the spatial distribution of the low density area was almost identical to the dark zones found in Figure 3.

Figure 4.  Point density map of the study site. The red ellipse denotes the area of low point density.

The false-color composite using bands 7, 4, and 2 of the Landsat ETM+ image from 2001 depicts the study site, shown in Figure 5 (left). This combination provides a ‘natural’ rendition. Healthy vegetation will be a bright green and can saturate in seasons of heavy growth; grasslands will appear green; pink areas represent barren soil; oranges and browns represent sparsely vegetated areas and dry vegetation; and water will be blue. This false-color image shows that there are distinct variations of biomass and biological vigor in the vegetated areas in our study site. Based on Landsat remote sensing data, biomass and biological vigor variations which may represent areas of disappeared reservoirs could be suitably identified by exploiting vegetation indices that are spectral combinations of different bands. Such indices are quantitative measures, based on vegetation spectral properties, which attempt to measure biomass or vegetative vigor. Vegetation indices are mainly derived from reflectance data from the red and near infrared (NIR) bands. The simplest form of vegetation index is the ratio between two digital values from the red and NIR spectral bands. The most widely used index is the well-known NDVI obtained using the following formula:

Figure 5.  Subset of the Landsat ETM + image of the study site (left) and NDVI map (right). The red ellipse denotes a strong vegetation area, covering almost the same area as the low backscattering coefficient area shown in Figure 3.

where Red and NIR stand for the digital number (DN) of the IR and NIR channels, respectively. NDVI is indicative of plant photosynthetic activity and has been found to be related to the green leaf area index and the fraction of photosynthetically active radiation absorbed by vegetation. High (low) values of the vegetation index identify pixels covered by substantial proportions of healthy (disease or stressed) vegetation. It is expected that biomass and biological vigor variations created by vegetation patterns should be emphasized by using NDVI. The NDVI of the study area is shown in Figure 5 (right).

In order to enhance the low contrast between archeological features and the surrounding areas, a principal component analysis (PCA) was conducted. PCA has proved to be of value in the analysis of multispectral, remotely sensed data. It is a linear transformation which decorrelates multivariate data by translating and/or rotating the axes of the original feature space, so that the data can be represented without correlation in a new component space (Richards and Jia 2006, Masinia et al. 2009). The PCA transforms the input multispectral bands into new components that should be able to make the identification of distinct features and surface types easier. The major portion of the variance in a multispectral data set is associated with homogeneous areas, whereas localized surface anomalies will be enhanced in later components, which contain less of the total dataset variance. This is the reason why latter components may represent information variance for a small area (Lasaponara and Masini 2007). The transformation of the raw remote sensing data using PCA can result in new PC images that may be more interpretable than the original data. In our work, bands 1–5 and 7 of Landsat ETM+ from 31 May 2001 were employed in PCA, and the output was six components. The principle component 2 (PC2), PC3, PC4, and PC5 images are shown in Figure 6. PC1 and PC6 contain a high proportion of correlated information from the homogeneous areas and noise and are not necessary to show here. PC2 highlights subtle variations contained in the multispectral bands in denoted areas, while PC3 highlights more, PC4 the most, and PC5 only subtly (Figure 6). Also, it is interesting to find that PC4 shows almost identical anomalous areas of possible disappeared reservoirs in our study site.

Figure 6.  Landsat ETM+ PCA map: upper left, PC2; upper right, PC3; lower left, PC4; lower right, PC5. The red ellipse denotes greater contrast than the surrounding area.

6. Verification

An excellent source of remote sensing information for archeologists is archival aerial photography and historical maps (the older the better). Archeologists have been using aerial prospecting and photography to locate buried structures, roads, and other features of ancient landscapes. In this study, however, the historical map and aerial photographs are used for verification of the detection results by Radarsat-1 and Landsat imagery. The old surveying map and aerial photograph used here were produced a century ago and half a century ago, respectively, and have been manually and digitally processed (Table 2).

Table 2. Historical materials used in verification of the detection results.

Data	Date	Resolution
Surveying map	1907	16.95 m
Aerial photograph	1954	16 m

The old surveying map was scanned and compiled at a scale of 1:200,000 at 400 dpi, with each sheet covering an area of 1°×40′, or over 6000 km². The map was georeferenced by corner coordinates in the original conical projection, projected into the Universal Transverse Mercator projection used for this study, and produced as digital images with a pixel size of 16.95 ground meters (about one-tenth the size of the spatial error expected for maps of this scale under modern map accuracy standards). While the accuracy of these prints must be corroborated by other sources, this surveying map represents the earliest systematic, georeferenced data covering the study region. Features about the disappeared reservoirs have been located and mapped, shown in Figure 7 and 8 It can be seen from Figure 7 that in our study site there were three individual reservoirs, namely Shushan, Nanwang, and Mata. Among the three, Shushan reservoir was the biggest and deepest, which explains why in the resulting detection images, such as Figure 3, Figure 5 (right), and Figure 6 (lower left), the lower left of the center shows anomalies more clearly. Figure 8 shows the existence of the two disappeared reservoirs, Nanwang and Mata. That is to say, as late as the 1950s, the Nanwang reservoir was still well maintained, and the Mata reservoir was shrinking and was about to return to field.

Figure 7.  The historic map above shows the location of the disappeared reservoirs. The red ellipse denotes the existence of three reservoirs, Shushan, Mata and Nanwang.

Figure 8.  The aerial photo above shows the location of the disappeared reservoirs. It shows that in the 1950s the Nanwang and Mata reservoirs still remained, especially the Nanwang reservoir.

7. Conclusions

Ancient archeological remains will affect their surrounding landscapes in different ways, whether through altering the surrounding soils, affecting vegetation, or absorbing moisture at a greater or lesser rate. The primary aim of archeological remote sensing is to find contrasts between archeological features and the natural background. In our work, Radarsat-1 and Landsat imagery were employed in detecting disappeared reservoirs through soil moisture variations and biomass and biological vigor variations linked to the presence of underlying ancient remains. Through the analysis of remotely sensed SAR and optical imagery of archeological sites, we can gain a more acute understanding of how varying intensities in images may relate to the sites of disappeared reservoirs. In order to delineate soil moisture variations, a point density analysis was conducted, and NDVI and PCA were conducted to delineate the areas of biomass and biological vigor variations. The detection result was verified by an historical map and aerial photographs. The study showed that SAR imagery is very useful for detecting buried archeological remains, especially when integrated with optical imagery for further analysis. Some interesting detection capabilities were observed in our work. We found that Radarsat-1 is more accurate for detecting archeological remains, especially with further analysis such as density slicing, point density, and reclassifying, as expected due to its sensitivity to soil moisture difference in the target area. Radarsat-1 shows more power in discovering and locating possible disappeared reservoirs in that the information it contained such as soil moisture and settlement distribution were more closely related to archeological elements. For the Landsat detection of archeological remains, the latter components of PCA, such as component 4, may enhance localized surface anomalies because a minor portion of the variance is associated with heterogeneous areas, usually indicating the presence of archeological remains. The outcome of this work enabled us to better understand the spatial distribution and historical extent of the disappeared reservoirs and potential archeological remains in the surrounding areas. This is of great importance to preservation and protection of the Grand Canal in that the disappeared reservoirs were inalienable parts of the Grand Canal in history.

Although remote sensing techniques may prove, in certain instances, an efficient means to obtain information on feature distribution, cost and efficiency must be weighed against the reliability and completeness of results. Remote sensing detection of archeological sites is an aid to discovery and interpretation that allows the field archeologist to focus research areas and goals, and to understand regional human–environment interactions. Traditional methods, however, including field expeditions, are necessary to ground-truth and precisely identify activity noted from a remote platform. Remote sensing techniques do not allow for the characterization of features and must, therefore, be combined with a program of selective excavation or exposure of features. They are a complement to, not a substitute for, subsurface testing. Therefore, the research will continue in several directions. Additional and better data sources, for example, high-resolution multispectral (or even hyperspectral) satellite images, multipolarization, or multifrequency radar images, and Lidar for better elevation data in test sites will be employed. Also, the analyses and processing may be improved with additional archeological knowledge derived from future research, including the methods of predictive modeling. Finally, in-depth site campaigns in conjunction with archeological knowledge are most needed.

Notes on contributors

Dr. Biao Deng is a research assistant at the Institute of Remote Sensing Applications, Chinese Academy of Sciences. He holds B.C. surveying and mapping degree from the Central South University and PhD in Cartography and GIS from the Graduate University of Chinese Academy of Sciences. His research interests lie in remote sensing applications in archeology.

Prof. Yueping Nie is the Director of the Department of Remote Sensing in archeology of Institute of Remote Sensing Applications, Chinese Academy of Sciences. He holds B.C. hydrogeology and engineering geology degree from Guizhou University of Technology and PhD in Cartography and Geography Information System from Nanjing University.

Prof. Huadong Guo is the Director-General of the Center for Earth Observation and Digital Earth, Chinese Academy of Sciences. He is Secretary General of the International Society for Digital Earth (ISDE), Editor-in-Chief of International Journal of Digital Earth, and president of the Committee on Data for Science and Technology (CODATA). He is currently a Principal Scientist of National Basic Research Project on ‘Research on the Mechanism and Methodology for Sensitive Factors of Global Change Using Space Technologies,’ and was Principle Investigator for more than 20 key national projects or programs in China, also PI for over 10 international remote sensing projects. He has published over 200 papers and 13 books. He is the recipient of three national awards on science and technology, and nine awards by the Chinese Academy of Sciences.

Prof. Changlin Wang, obtained his PhD degree from Leicester University, UK in 1995, and did post-doctoral research at Institute of Remote Sensing Applications, CAS. He was an associate Professor at this institute with main researches on radar remote sensing and remote sensing for archeology. He took a professorship position at the Center for Earth for Earth Observation and Digital Earth, CAS in 2007 and now is secretariat Director of the International Society for Digital Earth.

Dr. Shenglin Lei is an associate professor at the Center for Remote Sensing and Aeronautic Photographic Archeology of National Museum of China. He holds master's degree in History and Archeology from the Peking University. His research interests lie in application of spatial information technology to archeology. He has published several monographs and reports.

Dr. Rong Li is a research assistant at the Center for Remote Sensing and Aeronautic Photographic Archeology of National Museum of China. He holds PhD in History and Archeology from the Peking University. His research interests lie in application of spatial information technology to archeology. He has published several monographs and reports.

Acknowledgements

This work is part of interdisciplinary research programs (grant number: 41001272 and 2010BAK67B07) under the auspices of the National Science Foundation of China and Ministry of Science and Technology of China, respectively.