Characterisation of Terrain Variations of an Underwater Ancient Town in Qiandao Lake

: The underwater ancient town of Chunan is of great importance in archaeology and tourism. Hence, the e ﬃ cient mapping and monitoring of the topographical changes in this town are essential. An attractive choice for the e ﬃ cient mapping of underwater archaeology is the multibeam echo sounder system (MBES). The MBES has several advantages including noncontact survey, high precision, and low cost. In this study, the topographical changes of the ancient town under Qiandao Lake were quantitatively assessed on the basis of time-series MBES data collected in 2002 and 2015. First, the iterative closest point (ICP) algorithm was applied to eliminate the coordinate deviations between two point sets. Second, the robust estimation method was used to analyse the characterisations of the terrain variations of the town on the basis of the di ﬀ erences between the two matched point sets. Obvious topographical changes ranging from − 0.89 m to 0.88 m were observed in a number of local areas in the town. On the global scale, the mean absolute value of the depth change in the town was merely 0.12 m, which indicated a weak global deformation pattern. The experiment proved the e ﬀ ectiveness of applying MBES data to analyse the deformation of the ancient town. The results are beneﬁcial to the study of underwater ancient towns and the development of protection strategies.


Introduction
Underwater archaeology is of great importance in the study of the evolution of society and the culture of human beings over long periods [1][2][3][4][5]. However, underwater cultural materials and heritage sites in a complex and hostile underwater environment are threatened by physical, chemical, and biological factors such as dissolved oxygen, various ions, and water flow [6][7][8]. Moreover, underwater archaeological sites are inadvertently or deliberately damaged due to human activities such as dredging and clearing, deep trawling, and other fishing activities [9]. Hence, protecting these treasures requires effective methods to monitor the status of cultural materials and heritage sites and to repair damaged parts and facilitate in situ stabilisation [10][11][12].
Given the nature of the underwater environment, underwater sites are inaccessible and thus present difficulties in the conduct of field surveys. Early underwater archaeological methods combined diving and land archaeological techniques. Divers were required to be equipped with ventilators and appropriate archaeological equipment when conducting underwater archaeological work [13] where such procedures are labour intensive and time consuming. Moreover, the measurement quality strongly relies on the professional skills of personnel. For instance, when the wreck at Uluburun was discovered, the divers had to conduct up to 22,413 dives to depths between 44 and 61 m [14]. With the development of acoustic and optical technology, devices such as scanning sonar can obtain data of underwater targets with high precision and high resolution [15]. In the Thistlegorm Project, underwater photogrammetry was used in the survey of a shipwreck and proved to be an important tool for acquiring high-quality images, even in an unfavourable underwater environment [16]. However, underwater photogrammetry fails to obtain the high-precision geolocation data of objects and geometric information. Moreover, the efficiency of underwater photogrammetry remains low, as it can only cover several meters along a track. Westley et al. [17] combined single-beam sonar with side-scan sonar to explore a practical and economical way to detect underwater objects. During the measurement, single-beam sonar was used to provide high-precision locations of the seabed, and side-scan sonar was used to provide high-resolution acoustic images of the seabed and target objects. Although single-beam sonar provides an accurate depth measurement, it only covers the nadir region along the track. Thus, its measuring efficiency is restricted. In contrast, the side-scan sonar provides swath coverage during measurement, but the spatial accuracy of the measurement data is rather limited. To improve measurement resolution and efficiency, Doneus et al. [18] conducted a large-scale underwater archaeological survey on an ancient villa in the northern Adriatic coast of Croatia by using airborne laser bathymetry (ALB) technology. The study proved the enormous potential of ALB applied to shallow water archaeology. However, the measurement depth and quality of the terrain data collected by ALB heavily rely on the depth and quality of water.
Amongst various underwater remote sensing approaches, the multibeam echo sounder system (MBES) is the most widely used in bathymetric surveying because of its high resolution, high accuracy, and high efficiency. Unlike single-beam sonar or underwater photogrammetry, the MBES provides wide swath coverage of up to three to four times the depth. As a noncontact remote sensing approach, the MBES provides a safe approach to surveying targets and operators. Therefore, the MBES is considered as a promising technique to survey underwater objects [19,20]. Plets et al. [21] surveyed the north of Joint Irish with the MBES to search the remains of a shipwreck and successfully extracted its geometric information. Ernstsen et al. [22] used the MBES to monitor the evolution of the dunes in the Gradyb tidal inlet channel of the Danish Wadden Sea and verified the feasibility of the MBES in analysing the time-varying topographical changes of underwater targets. Tassetti et al. [23] quantitatively assessed the time-varying changes of the coral reefs in a seabed by comparing multitemporal MBES data.
Quantitatively and efficiently mapping and monitoring underwater heritage sites is essential for providing relevant data for the development of related protection strategies. The current study focused on assessing the time-varying characterisations of an underwater ancient town in Qiandao Lake with MBES data collected in two periods over a 13-year gap. To address the coordinate biases between the two datasets, we used the iterative closest point (ICP) algorithm in calculating the registration parameters. Then, the robust estimation method was used to compute the uncertainty of the bathymetric data and analyse the topographical changes of the town. The remainder of this paper is organised as follows. Section 2 describes the methodology of the proposed algorithm for detecting areas of the underwater ancient town with terrain variations. Section 3 presents the computation results of the proposed method. Sections 4 and 5 present the discussions and conclusions, respectively.

Topographical Change Monitoring of Underwater Ancient Town with Multibeam Echo Sounder System (MBES)
In multibeam operation, the transmitting transducer sends beams downward with a wide opening angle sector. Then, the receiving transducer receives the reflected signals at various angles within the sector and generates beam footprints on the seabed. The location of footprints can be calculated by applying the corresponding mathematical model using relevant measurements including the receiving beam angle, the travel time of the signal from the transducer to the seabed, and other auxiliary information (e.g., sound velocity profile of water column, altitude data and positioning data of the vessel) [24,25]. With the calculated footprints, the digital elevation model (DEM) of the lake including the underwater ancient town can be generated [26].
Underwater topographical changes are generated by various factors such as crustal deformation, water currents, or even shipwrecks [27]. A topographical change monitoring methodology involves calculating and assessing the differences between DEMs obtained by two repeated MBES surveys. The resulting 'DEM of Difference' (DoD) quantifies the changes in terrain elevation with positive values indicating deposition, negative values indicating erosion, and null values indicating an unchanged terrain. Specifically, the two datasets of the same terrain can be represented as {h 1i , i = 1, . . . , n1} and {h 2i , i = 1, . . . , n2}, respectively, and the DoD at the i-th point dH i can be expressed as Theoretically, the DoD value can reflect the bathymetric changes of the terrain. However, measurement errors may affect the accuracy of the quantitative analysis of terrain changes. Errors in MBES measurements stem from a range of factors such as altitude changes of the transducer with the motion of carriers, the residual errors in instruments, and background noise [28]. Therefore, the DoD value could not be simply regarded as indicative of changes in terrain elevation, and the influence of measurement noise must be accounted for. Based on the assumption that the noise in MBES data follows a Gaussian distribution, the probability distribution of the DoD value is expressed as where σ 1 and σ 2 are the standard deviation (SD) of random errors introduced during the measurement. Therefore, topographical changes can be represented by the points that are statistically significant when the absolute value of dH is greater than a certain predefined threshold (e.g., 3σ dH ).

Coordinate Frame Registration
Ideally, the terrain variation in the two periods can be calculated and analysed on the basis of the aforementioned procedure as long as both point sets are correctly georeferenced. However, the deviation in the horizontal and vertical coordinate frames between the two point sets occurs because of inaccurate positioning and variation in the vertical direction [27,29,30]. If not addressed properly, the deviation could seriously affect the results of topographical change analysis. To address the problem, we applied the ICP algorithm to transfer the two point sets acquired by different sensors into the same common coordinate system. Figure 1 shows the basic flow of the ICP algorithm.
where R is a 3D rotation matrix and the T is a translation vector. The ICP algorithm attempts to find the best transformation relationship between point sets that aligns the source point set with the target point set. For this purpose, the following cost function E(R, T) is used: To derive the 3D rotation matrix R and translation vector T, the covariance matrix C is decomposed by using the singular value decomposition: where 1 S ， 2 S are the mean values of the point sets 1 S and 2 S , respectively.
An iterative search is conducted to derive the pre-established correspondence between the source point set and the target point set. At each iteration, a correspondence is established, and the transformation that minimises the mean square metric is computed. The data points are transformed using the computed transformation, and the process is repeated whilst reducing the given cost function E(R, T) until it is less than the given threshold, or until a maximum number of iterations is reached. A given surface can be described as a point set. The two point sets describing the same surface collected by different sensors could be represented in three dimensions as S 1 = {p 11 , p 12 , . . . , p 1n , p 1i ∈ R 3 and S 2 = p 21 , p 22 , . . . , p 2n , p 2i ∈ R 3 , respectively. Let S 1 be the source point set, and S 2 be the target point set; S 1 is not necessarily equal to S 2 . No pre-established correspondence exists between the point sets. One of the point sets generally undergoes certain transformations such as translation and rotation to match with the other point set. The equation can be expressed as where R is a 3D rotation matrix and the T is a translation vector. The ICP algorithm attempts to find the best transformation relationship between point sets that aligns the source point set with the target point set. For this purpose, the following cost function E(R, T) is used: To derive the 3D rotation matrix R and translation vector T, the covariance matrix C is decomposed by using the singular value decomposition: where S 1 , S 2 are the mean values of the point sets S 1 and S 2 , respectively. An iterative search is conducted to derive the pre-established correspondence between the source point set and the target point set. At each iteration, a correspondence is established, and the transformation that minimises the mean square metric is computed. The data points are transformed using the computed transformation, and the process is repeated whilst reducing the given cost function E(R, T) until it is less than the given threshold, or until a maximum number of iterations is reached.

Detection of Terrain Variations
With the ICP algorithm, two different temporal point sets can be registered in the same coordinate system. Then, the dynamic information of the terrain can be calculated on the basis of the SD of dH (σ dH ) in the region. In theory, σ dH can be obtained by [Equation (2)], however, accurately modelling σ dH for bathymetric data, whose uncertainty knowledge is rarely available, is difficult because of the numerous factors during MBES measurement that may complicate the noise model in unfavourable underwater conditions [31]. Therefore, we focused on the local statistical distribution of dH and estimated the σ dH value from a phenomenon point of view. Specifically, the whole study area was divided into several subareas by using a circular sliding window, with each point in the point set as the centre point with a radius of 5 m. Then, the SD of dH was calculated locally in the window.
Two factors are considered in choosing the size of the sliding window. First, given a small window, the resulting variance estimate would show a strong statistical fluctuation due to the limited size of the samples in the window. Second, a large window would weaken the representation of the variance estimate. Thus, a trade-off exists in choosing the size of the sliding window. In our study, the radius of the sliding window was set as 5 m. In each circular window, the SD of local dH was computed as where σ dH k is the SD of the k-th window; dH ki is the dH value of the i-th point in the k-th window; and dH mean is the mean value in the window given as Generally, the aforementioned statistics are efficient for Gaussian-distributed variables. However, for our datasets, the potential existence of topographical changes and outliers in measurements are not consistent with such a Gaussian assumption. If not addressed properly, then the computed variance will be seriously biased by the outliers [32,33]. Hence, the variance of depth difference in the local window in the current work was estimated using the robust estimation procedure. Specifically, the robust counterparts of the mean and the SD statistics were calculated. The empirical formula for the SD of the local region is represented aŝ where dĤ k is the median value in the k-th window.
Given the variance estimate, if the points in the window have not undergone any change, then the distribution of dH values will conform to a standard normal distribution. If the absolute value of a depth difference value exceeds a certain threshold (3σ in this study), then the point will be considered having undergone a significant change. Thus, by traversing all the depth difference values based on Equations (12)- (14), the points can be flagged [24].
where S k is the circular region with a radius of 5 m. α is the confidence level and µ α 2 is the corresponding value when the confidence level is equal to α/2.

Study Area and Data Collection
The study area is the ancient town of Chunan, which is submerged in Qiandao Lake. Located in Hangzhou City, China, Qiandao Lake is a vast inland lake with an area of approximately 567.4 km 2 and an average depth of about 34 m (Figure 2). In September 1959, the town of Lion City in Chunan County, with a long history of over a thousand years, was submerged underwater due to the building of the Xin'anjiang Hydroelectric Plant. Significant underwater terrain variation can be expected due to the damage of the ancient buildings, as shown in Figure 3. Hence, the MBES should be employed for repeated measurements to monitor and evaluate the temporal topographical variations of the town.  (14) where Sk is the circular region with a radius of 5 m. α is the confidence level and 2 α μ is the corresponding value when the confidence level is equal to α /2.

Study Area and Data Collection
The study area is the ancient town of Chunan, which is submerged in Qiandao Lake. Located in Hangzhou City, China, Qiandao Lake is a vast inland lake with an area of approximately 567.4 km² and an average depth of about 34 m (Figure 2). In September 1959, the town of Lion City in Chunan County, with a long history of over a thousand years, was submerged underwater due to the building of the Xin'anjiang Hydroelectric Plant. Significant underwater terrain variation can be expected due to the damage of the ancient buildings, as shown in Figure 3. Hence, the MBES should be employed for repeated measurements to monitor and evaluate the temporal topographical variations of the town.
.  Two datasets collected at different times were used for the analysis in this study. The first dataset was collected in 2002 with an EM3000 system. In the survey, the swath angle was set to 100°, and the beam number and depth resolution were 127 and 1 cm, respectively. The average point density was equal to 1.2 points per m², and the average value of the depth measurements was 31.54 m. The second dataset was collected in 2015 using the SONIC2024 system. In the survey, the swath angle was about 120°, and the beam number and depth resolution were 256 and 1.25 cm, respectively. To ensure that the whole terrain was surveyed, the interval distance between survey lines was set corresponding to the depth to provide overlap in the adjacent lines. The average point resolution was 1.5 m, and the average depth was 30.75 m.
Prior to analysis, the bathymetric data were processed using the commercial software CARIS HIPS/SISP 9.1. The specific procedure was as follows: (1) swath edit (mainly removing the noise points and outliers manually); (2) sound velocity correction; (3) tidal correction; (4) data fusion; (5) calibration of the installation parameters; and (6) data refusion and exporting of bathymetric data. Subsequently, the bathymetric data were converted into the XYZ format. Gridded DEMs with a resolution of 1 m were generated by interpolating the point cloud with the Golden Surfer 16. 0 software, as shown in Figure 4. Visual judgment revealed that the main features around the town including the walls, roads, pond, and private houses were well preserved and visible. The current status of the Lion City is in full compliance with the historical records, indicating an area of about 0.432 km² and wall length of approximately 2.57 km. For easy visualisation and analysis, the town and surrounding areas were divided into four regions: A, B, C, and D. Regions A and B were set as the marginal zones of the outer wall of the town with a relatively flat terrain. Region C was characterised as being covered by mountains with great topographic relief. Region D represented the ancient town with visible outlines of the streets and walls. Two datasets collected at different times were used for the analysis in this study. The first dataset was collected in 2002 with an EM3000 system. In the survey, the swath angle was set to 100 • , and the beam number and depth resolution were 127 and 1 cm, respectively. The average point density was equal to 1.2 points per m 2 , and the average value of the depth measurements was 31.54 m. The second dataset was collected in 2015 using the SONIC2024 system. In the survey, the swath angle was about 120 • , and the beam number and depth resolution were 256 and 1.25 cm, respectively. To ensure that the whole terrain was surveyed, the interval distance between survey lines was set corresponding to the depth to provide overlap in the adjacent lines. The average point resolution was 1.5 m, and the average depth was 30.75 m.
Prior to analysis, the bathymetric data were processed using the commercial software CARIS HIPS/SISP 9.1. The specific procedure was as follows: (1) swath edit (mainly removing the noise points and outliers manually); (2) sound velocity correction; (3) tidal correction; (4) data fusion; (5) calibration of the installation parameters; and (6) data refusion and exporting of bathymetric data. Subsequently, the bathymetric data were converted into the XYZ format. Gridded DEMs with a resolution of 1 m were generated by interpolating the point cloud with the Golden Surfer 16. 0 software, as shown in Figure 4. Visual judgment revealed that the main features around the town including the walls, roads, pond, and private houses were well preserved and visible. The current status of the Lion City is in full compliance with the historical records, indicating an area of about 0.432 km 2 and wall length of approximately 2.57 km. For easy visualisation and analysis, the town and surrounding areas were divided into four regions: A, B, C, and D. Regions A and B were set as the marginal zones of the outer wall of the town with a relatively flat terrain. Region C was characterised as being covered by mountains with great topographic relief. Region D represented the ancient town with visible outlines of the streets and walls.

Registration of Point Sets
The systematic residuals and calibration residuals between two datasets collected by different types of MBES may present differences. Hence, the coordinate frames of the two datasets will be inevitably inconsistent. Therefore, the ICP algorithm was applied in this work to reduce the coordinate deviation. Unless the residual value reaches the threshold, the iteration in the ICP algorithm is terminated to ensure optimal registration. As shown in Figure 5, region C, which was covered by mountains with obvious topographical relief, was extracted to highlight the area with a significant contrast effect from the ICP algorithm. Figure 5 shows that relative to the right image in Figure 5b, the morphological features in the left image in Figure 5a are blurry. Hence, detecting topographical changes was challenging. After being registered by the ICP algorithm, the root mean square (RMS) of DoD dropped from 0.25 m to 0.13 m. The terrain features are marked by red frames in Figure 5b.

Registration of Point Sets
The systematic residuals and calibration residuals between two datasets collected by different types of MBES may present differences. Hence, the coordinate frames of the two datasets will be inevitably inconsistent. Therefore, the ICP algorithm was applied in this work to reduce the coordinate deviation. Unless the residual value reaches the threshold, the iteration in the ICP algorithm is terminated to ensure optimal registration. As shown in Figure 5, region C, which was covered by mountains with obvious topographical relief, was extracted to highlight the area with a significant contrast effect from the ICP algorithm. Figure 5 shows that relative to the right image in Figure 5b, the morphological features in the left image in Figure 5a are blurry. Hence, detecting topographical changes was challenging. After being registered by the ICP algorithm, the root mean square (RMS) of DoD dropped from 0.25 m to 0.13 m. The terrain features are marked by red frames in Figure 5b.

Registration of Point Sets
The systematic residuals and calibration residuals between two datasets collected by different types of MBES may present differences. Hence, the coordinate frames of the two datasets will be inevitably inconsistent. Therefore, the ICP algorithm was applied in this work to reduce the coordinate deviation. Unless the residual value reaches the threshold, the iteration in the ICP algorithm is terminated to ensure optimal registration. As shown in Figure 5, region C, which was covered by mountains with obvious topographical relief, was extracted to highlight the area with a significant contrast effect from the ICP algorithm. Figure 5 shows that relative to the right image in Figure 5b, the morphological features in the left image in Figure 5a are blurry. Hence, detecting topographical changes was challenging. After being registered by the ICP algorithm, the root mean square (RMS) of DoD dropped from 0.25 m to 0.13 m. The terrain features are marked by red frames in Figure 5b.   Figure 6 shows the residual variation during the ICP computation process. The residuals around the whole study area dropped from 1.44 m to 0.21 m. By minimising the residuals, the bias between the two datasets in the X and Y direction were identified as 0.01 and 0.16 m, respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 15 Figure 6 shows the residual variation during the ICP computation process. The residuals around the whole study area dropped from 1.44 m to 0.21 m. By minimising the residuals, the bias between the two datasets in the X and Y direction were identified as 0.01 and 0.16 m, respectively.

Regional Analysis
With two point sets being registered, six different statistics (Table 1) of DoD in the four regions, (i.e., A, B, C and D, Figure 4) were calculated. On the global scale, the temporal variation of terrain is relatively small, as indicated by the mean, the mean absolute value (MAV), and the root mean square (RMS) ( Table 1), while from a local perspective, there is an evident difference in the bathymetric variation patterns among different local regions. Generally, the MAVs and RMS in regions A and B are greater than those in regions C and D. This indicates that the temporal variation in bathymetry is more evident in regions A and B. Specifically, the mean depth change values were −0.01 m for region A, which presented an erosion feature; and 0.04 m in region B, which indicated a deposition pattern. As for regions C and D, the mean depth change values were both negligible, which indicates insignificant temporal bathymetry changes around these regions.
Similar patterns were also revealed by the histograms of the DoD values. As shown in Figure 7, evident bathymetry variations are presented in regions A and B. Specifically, the heavy tail in the left hand side of the histogram indicates a general erosion pattern in region A. On the contrary, region B is shown to have a significant deposition pattern, while for regions C and D, the histograms were generally symmetrical, and are similar to the standard normal distribution.

Regional Analysis
With two point sets being registered, six different statistics (Table 1) of DoD in the four regions, (i.e., A, B, C and D, Figure 4) were calculated. On the global scale, the temporal variation of terrain is relatively small, as indicated by the mean, the mean absolute value (MAV), and the root mean square (RMS) ( Table 1), while from a local perspective, there is an evident difference in the bathymetric variation patterns among different local regions. Generally, the MAVs and RMS in regions A and B are greater than those in regions C and D. This indicates that the temporal variation in bathymetry is more evident in regions A and B. Specifically, the mean depth change values were −0.01 m for region A, which presented an erosion feature; and 0.04 m in region B, which indicated a deposition pattern. As for regions C and D, the mean depth change values were both negligible, which indicates insignificant temporal bathymetry changes around these regions.
Similar patterns were also revealed by the histograms of the DoD values. As shown in Figure 7, evident bathymetry variations are presented in regions A and B. Specifically, the heavy tail in the left hand side of the histogram indicates a general erosion pattern in region A. On the contrary, region B is shown to have a significant deposition pattern, while for regions C and D, the histograms were generally symmetrical, and are similar to the standard normal distribution. On the basis of robust estimation, the specific areas experiencing topographical changes were detected. First, the variances of the DoD in the local sliding window were calculated. Second, the confidence level of the DoD was set to 99.7% to identify the areas with significant bathymetry variations. Figure 8 shows the spatial distribution of the identified topographical variations. In region A, a mix of erosion and light deposition was observed. The deposition was predominant in region B, and an increasing tendency was revealed along the wall. On the contrary, few topographical changes were observed in regions C and D. To better illustrate the bathymetry variation pattern, four representative profiles are shown for their DoD and bathymetry information (Figure 9). On the basis of robust estimation, the specific areas experiencing topographical changes were detected. First, the variances of the DoD in the local sliding window were calculated. Second, the confidence level of the DoD was set to 99.7% to identify the areas with significant bathymetry variations. Figure 8 shows the spatial distribution of the identified topographical variations. In region A, a mix of erosion and light deposition was observed. The deposition was predominant in region B, and an increasing tendency was revealed along the wall. On the contrary, few topographical changes were observed in regions C and D. To better illustrate the bathymetry variation pattern, four representative profiles are shown for their DoD and bathymetry information (Figure 9).

Various Threats to the Ancient Town
The Lion City is an ancient town with a history of over a thousand years. It is a part of a large-scale ancient community that includes several other ancient towns (the Crane City, the Weiping Town, the Port Town, and the Tea Plantation Town). All of the ancient towns have been submerged in the Xin'anjiang Reservoir for more than 60 years. Generally, the terrain in the reservoir environment is often subjected to significant bathymetry variation due to the influence of underwater currents [34,35]. Moreover, human activities in Qiandao Lake such as dredging and fishing activities can also cause an evident variation of underwater topography. All of these factors pose a significant threat to the conservation of the underwater ancient towns. For instance, the walls and architectures built with vulnerable bricks and glued with lime and other materials may collapse due to the pressure generated by the near-bottom current. Thus, monitoring the temporal topography variation of the ancient towns based on repeated surveys can provide critical information for the heritage conservation management.

Validation of Analysis
When comparing two repeated bathymetric datasets, the differences in depth resulting from measurement errors can be approximated by a normal distribution. In contrast, the morphological changes of underwater targets can be regarded as outliers. Thus, the temporal changes of the bathymetry around the ancient town can be identified by using hypothesis testing based on DoD values. The robust estimation revealed slight terrain variations in the study area on the global scale. The mean and MAV of DoD value were 0 and 0.10 m, respectively. This finding is consistent with the conclusion of Zhang that the geological environment of Qiandao Lake is stable and that the corresponding annual crust deformation is around 1.1 cm [36]. On the contrary, large terrain variations were identified in some local regions. The spatial distribution of the DoD (Figure 8) indicates that the locations with evident terrain variations are located mainly outside the town wall in regions A and B. This pattern is presumably caused by the influence of the interaction between the underwater current and the wall. For sites located inside the wall, the cause of the significant terrain variation is probably the architectural deformation. In conclusion, all of these findings suggest that the ancient town located underwater is subjected to significant stress due to the influence of hydrological and chemical factors.

Monitoring Suggestions for Ancient Town
The proposed method provides a practical way to understand the characterisations of the topographical changes in the study area. Combining the data collected by the MBES with multitemporal analysis is a tractable, effective, and labour-saving method to determine the dynamics of the town. However, the uncertainty in MBES datasets may prevent the detection of actual surface changes as these datasets are based merely on dH values containing all the information on topographical changes [37]. However, the current methods of estimating topographical changes based on MBES data including determining the difference between two MBES datasets under a stable feature assumption [38,39] or subtracting crest lines and trough lines [40,41], are unable to completely remove the disturbance of outliers; hence, the assessment of terrain variations remains challenging. In addition, these methods do not consider the uncertainty in their calculation. Therefore, we suggest that before using the robust estimation, the impacts of outliers should be mitigated as much as possible. In future research, data from various sources should be used to enhance the accurate characterisation of the dynamics of ancient towns.

Conservation Suggestion of Ancient Towns
The ancient town of Chunan has been submerged in Qiandao Lake for more than 60 years. Recently, the value of this ancient town in culture and tourism has been recognised by the local government. In 2011, Lion City was listed as a heritage conservation unit of Zhejiang Province. However, shipping, fishing, and dredging are still prevalent in Qiandao Lake. Moreover, the frequent variations in water levels due to the storage and discharge of the Xian Reservoir can generate large underwater currents. Although Section 3 showed that the overall depth variation was insignificant, a significant deformation pattern was identified for the local area inside and outside the town wall. Thus, protecting the underwater ancient town requires a strict and targeted policy to balance economic development and heritage conservation. In this context, repeated MBES surveys can effectively collect information about the state of the underwater architecture. Thus, appropriate policy making and management will require an increased frequency of MBES bathymetry surveys around the area.

Conclusions
Repeated MBES surveys provide valuable information about the topography variation of the ancient town of Chunan. In this study, we combined the ICP algorithm and the robust estimation method to characterise the time-varying topographical features of the ancient town between 2002 and 2015. The computation result showed an insignificant global terrain change pattern. In contrast, large topographical changes were identified in a number of local areas. The experimental results verified the feasibility of applying the MBES to effectively quantify the time-varying topographical changes of the underwater ancient town. Furthermore, the result is of reference value for formulating government policies about the ancient town. In future study, the characterisation of the dynamics of the ancient town can benefit greatly from long-term MBES surveys with an increased frequency.