Coastal Change Patterns from Time Series Clustering of Permanent Laser Scan Data

. Sandy coasts are constantly changing environments governed by complex, interacting processes. Permanent laser scanning is a promising technique to monitor such coastal areas and to support analysis of ge-omorphological deformation processes. This novel technique delivers 3D representations of the coast at hourly temporal and centimetre spatial resolution and allows to observe small scale changes in elevation over extended periods of time. These observations have the potential to improve understanding and modelling of coastal deformation processes. However, to be of use to coastal researchers and coastal management, an efﬁcient way to ﬁnd and extract deformation processes from the large spatio-temporal data set is needed. To enable automated data mining, we extract time series of surface elevation and use unsupervised learning algorithms to derive a partitioning of the observed area according to change patterns. We compare three well known clustering algorithms, k-means, agglomerative clustering and DBSCAN, apply them on the set of time series and identify areas that undergo similar evolution during one month. We test if these algorithms fulﬁl our criteria for suitable clustering on our exemplary data set. The three clustering methods are applied to time series over 30 days extracted from a data set of daily scans covering about two kilometres of coast at Kijkduin, the Netherlands. A small section of the beach, where a pile of sand was accumulated by a bulldozer is used to evaluate the performance of the algorithms against a ground truth. The k-means algorithm and agglomerative clustering deliver similar clusters, and both allow to identify a ﬁxed number of dominant deformation processes in sandy coastal areas, such as sand accumulation by a bulldozer or erosion in the intertidal area. The level of detail found with these algorithms depends on the choice of the number of clusters k . The DBSCAN algorithm ﬁnds clusters for only about 44 % of the area and turns out to be more suitable for the detection of outliers, caused for example by temporary objects on the beach. Our study provides a methodology to efﬁciently mine a spatio-temporal data set for predominant deformation patterns with the associated regions, where they occur.

in terms of standard deviation within the clusters and area of the beach covered by the clustering. We compare and evaluate 90 the resulting clusters using these criteria as a first step towards the development of a method to mine the entire data set from permanent laser scanning for deformation processes.
2 The permanent laser scan data set shown as an example in Figure 3. The test area, which is discussed in Section 3.4, is indicated in red at the end of the northern path leading to the beach.
The data set from permanent laser scanning is acquired within the CoastScan project at a typical urban beach in Kijkduin, the Netherlands, . For the acquisition a Riegl VZ-2000 laser scanner was used to scan over a period of six 95 months from December 2016 to May 2017. The full data set consists of hourly scans of a section of sandy beach and dunes.
For the present study, a subset of the available data is used to develop the methodology. This subset consists of 30 daily scans taken at low tide over a period of one month, January 2017. It covers a section of the beach and dunes in Kijkduin and is displayed in top view in Figure 1. The area contains a path and stairs leading down to the beach, a paved area in front of the dunes, a fenced in dune area and the sandy beach. It is about 950 m long, 250 m wide and the distance from the scanner to the 100 4 https://doi.org/10.5194/esurf-2020-34 Preprint. Discussion started: 19 May 2020 c Author(s) 2020. CC BY 4.0 License. farthest points on the beach is just below 500 m. For the duration of the experiment the scanner was mounted on the roof of a hotel just behind the dunes at a height of about 37 m above sea level (as shown in Figure 2).
The data is extracted from the laser scanner output format and converted into a file that contains xyz-coordinates and spherical coordinates for each point. The data is mapped into a local coordinate system, where the origin in x-and y-direction is at the location of the scanner and the height (z-coordinate) corresponds to height above sea level. Since we are interested in relative 105 changes between consecutive scans, we do not transform the data into a geo-referenced coordinate system for this analysis.
Each point cloud is chosen to be at the time of lowest tide between 18:00 and 06:00, in order to avoid people and dogs on the beach, with the exception of two days where only very few scans were available due to maintenance activities. The data from 9 th of January 2017 is entirely removed from the data set, because of poor visibility due to fog. Additionally all points above 14.5 m elevation are removed to filter out points representing the balcony of the hotel and flag posts along the paths that 110 are interfering with the spherical coordinate grid extraction. In this way also a majority of reflections from particles in the air, birds or raindrops are removed. However, some of these particles might still be present at lower heights close to the beach.  Figure 1) on the path that is assumed to be stable throughout the entire month. Elevation is varying within less than 2 cm.
Since the data is acquired from a fixed and stable position we assume that consecutive scans are aligned. Nevertheless, the orientation of the scanner may change slightly due to strong wind, sudden changes in temperature, or maintenance activities. 5 https://doi.org/10.5194/esurf-2020-34 Preprint. Discussion started: 19 May 2020 c Author(s) 2020. CC BY 4.0 License.
The internal inclination sensor of the scanner measures these shifts while it is scanning and we apply a correction for large 115 deviations (more than 0.01 degrees) from the median orientation.
The remaining error in elevation and range is estimated as the average standard deviation of time series in locations that are assumed to be stable during the entire month. We chose stable surfaces that are part of the paved paths on top of the dunes and leading to the beach in northern and southern direction and derived the mean remaining errors shown in Table 1. The elevation does on average not deviate more then 2 cm from the mean elevation of the respective area and the standard deviation of the 120 range is within 5 cm on the northern path and on top of the dunes, but around 11 cm on the southern path. An example time series from the stable paved area on top of the dunes (at location (x 0 , y 0 ) marked in Figure 1) is shown in Figure 3.

Methods
To derive coastal deformation processes from clusters based on change patterns we follow three steps: Extraction of time series in two different coordinate systems, clustering of time series with three different algorithms, and derivation of geomorpholog-125 ical deformation processes. To cluster time series the definition of a distance between two time series (or the similarity) is not immediately obvious. We discuss two different options (Euclidean distance and correlation) to define distances between time series with different effects on the clustering results. The rest of this section is organized as follows: We focus on time series extraction in subsection 3.1, discuss distance metrics for time series (3.2), introduce three clustering algorithms (3.3) and our evaluation criteria (3.4). The derivation of deformation processes will be discusses with the results (section 4).

Time Series Extraction
Time series are extracted from the PLS data set in two different ways: First by using a grid in Cartesian xy-coordinates and extracting time series in elevation and second by using a grid in spherical coordinates and extracting time series in range. For both methods, we only use grid cells that contain at least one point for each of the scans.
A B

Cartesian Coordinates
Before defining a grid in Cartesian coordinates, we rotate the observation area to make sure that the coastline is parallel to the y-axis. This ensures that the grid covers the entire observation area efficiently and leaves as few empty cells as possible. Then we generate a regular grid (as illustrated in Figure 4) with grid cells of 1 m × 1 m. Time series are generated for each grid cell by taking the median elevation z i for each grid cell and for each time stamp t k . That means, per grid cell with center (x i , y i ) we have a time series with the number of time stamps T = 30. To make the time series dependent on change patterns, rather than the absolute elevation values, we remove the mean elevationz i of each time seriesZ i . This leads to time series

Spherical Coordinates
Using spherical coordinates allows to generate a grid with constant angle increment in horizontal and vertical direction with roughly constant point density per grid cell. For each grid cell j, we derive a time series R j consisting of the median range r j (t k ) per grid cell per time stamp t k :

155
where T = 30, as above. The range r j is defined as the line of sight distance from the laser scanner to the respective point. We choose grid cells of 0.1 • ×0.375 • to ensure that grid cells on the beach (close to the dune foot) cover roughly 1 m 2 , the same as in Cartesian coordinates, in order to make both methods comparable. However, transformed back into Cartesian coordinates, the size of the grid cells (in square-meters) varies with distance from the scanner (see Figure 4).
The point density in a point cloud is generally lower, the farther away a point is from the scanner. This property of our data 160 set is preserved using spherical coordinates and represented in the size of the grid cell, or distance between grid cell centres.

Distance Metrics
We consider two different distance metrics for our analysis.

Euclidean Distance
The most common and obvious choice is the Euclidean distance metric defined as: for two time series Z 0 and Z 1 of length n.

Correlation Distance
Another well known distance measure is correlation distance, defined as one minus the Pearson correlation coefficient (see for example Deza and Deza (2009)). It is a suitable measure of similarity between two time series, when correlation in the data is 170 expected (see Iglesias and Kastner (2013)). Correlation between two time series Z 0 and Z 1 is defined as: withZ being the mean value of time series Z and || · || the Euclidean 2-norm as in Equation (4). We have to note here, that correlation cannot compare simple constant time series (leads to division by zeros) and is therefore not a distance metric in the sense of the definition Deza and Deza (2009).

Comparison
For a comparison of the two distances for some example time series see Figure 5. The example shows that the distance between two time series is not intuitively clear. The use of different distance metrics results in different sorting of distances between the shown pairs of time series. However, when normalizing all time series (subtracting the mean and scaling by the standard deviation) correlation distance and Euclidean distance are equivalent (as shown for example by Deza and Deza (2009)).

180
Both Euclidean distance and correlation are not taking into account the order of the values within each time series. For example, two identical time series that are shifted in time are seen as 'similar' with the correlation distance, but not as similar with the Euclidean distance and would not be considered as identical by either of them (see Figure 5). Additionally neither of the two distance metrics can deal with time series of different lengths or containing gaps.

185
Clustering methods for Time Series can be divided into two categories: feature based and raw data based (see for example Liao (2005)). Feature based methods first extract relevant features to reduce dimensionality (for example using Fourier-or wavelet-transforms) and then form clusters based on these features. They could also be used to deal with gaps in time series.
We focus on the raw data based approach to not define features in advance and to make sure that no information within the data set is lost. We use three different methods: k-means clustering, agglomerative clustering and Density-Based Spatial Clustering 190 of Applications with Noise (DBSCAN). In Figure 6 an illustration of a partitioning of a simple 2D data set is shown for each of the three algorithms. The two clusters that can be distinguished in this example have different variances and are grouped differently by each of the algorithms.
For the implementation of all three algorithms, we make use of the Scikit-learn package in Python (see Pedregosa et al. (2011)).
195 Figure 6. Example of clustering of data with two clusters with different variance: The k-means algorithm separates them, but adds a few points in the middle to the purple cluster instead of the yellow one (A). Agglomerative clustering separates both clusters according to their variances (B) and DBSCAN detects the cluster with low variance and high point density (yellow) and discards all other points as outliers (turquoise) (C).

k-means Clustering
The k-means algorithm was first introduced in 1955 and is still one of the most widely used clustering methods (Jain (2010)).
The algorithm is based on minimizing the sum of all distances between points and centroids over all possible choices of k cluster centroids V = {v 1 , . . . , v k }:

200
with Euclidean distance metric || · ||. After the initial choice of k centroids among all points the following steps are repeated iteratively, until the above sum does not change significantly: 1. Assign each point to the cluster with closest centroid 2. Move centroid to mean of each cluster 3. Calculate sum of distances over all clusters (Equation (6)) 205 Note that minimizing the squared sum of distances over all clusters, coincides with minimizing the squared sum of all within cluster variances. The convergence to a local minimum can be shown for the use of Euclidean distance (see for example Jain (2010)). The convergence is sped up using so-called k-means++ initialization: After the random selection of the first centroid, all following centroids are chosen based on a probability distribution proportional to their squared distance to the already defined centroids. In this way the initial centroids are spread out throughout the data set and the dependence on the random 210 initialization of the cluster centroids is reduced.
There are variations of k-means using alternative distance metrics such as the L 1 -norm (k-medoids, Park and Jun (2009)), however the convergence is not always ensured in these cases. Another issue to take into account when considering alternative distance metrics, is the definition of the cluster centroids as mean of time series, which is not automatically defined for any distance metric. For more information on k-means see Jain (2010), Liao (2005) and the documentation of the Scikit-learn 215 package (Pedregosa et al. (2011)).

Agglomerative Clustering
Agglomerative clustering is one form of hierarchical clustering: It starts with each point in a separate cluster and iteratively merges clusters together until a certain stopping criterion is met. There are different variations of agglomerative clustering using different input parameter and stopping criteria (see for example Liao (2005) or the documentation of the scikit-learn 220 package (Pedregosa et al. (2011))). We choose the minimization of the sum of the within cluster variances using the Euclidean distance metric (Equation (6), where the centroids v j are the mean values of the clusters) for a pre-defined number of clusters k.
The algorithm starts with each point in a separate cluster and iteratively repeats the following steps until k clusters are found: 1. Loop through all combinations of clusters: -Calculate squared sum of distances (Equation (6)) for each combination 2. Keep clusters with minimal squared sum of distances In this way we use agglomerative clustering with a similar approach to the k-means algorithm, the same optimization criterion with the same input parameter and Euclidean distance measure. We therefore expect similar results. However, this agglomerative clustering can easily be adapted to alternative distance measures and could therefore potentially deal with time 230 series of different lengths or containing gaps.

DBSCAN Algorithm
Density-Based Spatial Clustering of Applications with Noise, DBSCAN, is a classical example of clustering based on the maximal allowed distance to neighbouring points that automatically derives the numbers of clusters from the data. It was introduced in 1996 by Ester et al. (1996) and recently revisited by Schubert et al. (2017). The algorithm is based on dividing 235 all points into core points or non-core points that are close to core points but not themselves surrounded by enough points to be counted as core points. The algorithm needs the maximum allowed distance between points within a cluster (ε) and the minimum number of points per cluster (N min ) as input parameters. These two parameters define a core point: If a point has a neighbourhood of N min points at ε distance, it is considered a core point. The algorithm consists of the following steps In this way clusters are formed that truly represent a dense collection of 'similar' points. Since we choose to use correlation 245 as distance metric, each cluster will contain correlated time series in our case. All points that can not be assigned to a close surrounding of a core point, are classified as noise or outliers.

Evaluation Criteria
To determine if an algorithm is suitable, we expect that it fulfils the previously defined criteria: -A majority of the observation area is separated into distinct regions,

250
each cluster shows a change pattern that can be associated with a geomorphic deformation process, and time series contained in each cluster roughly follow the mean change pattern.
In order to establish these criteria, we compare the three clustering algorithms, as well as the two different ways to derive time series, using the following evaluation methods.

255
The clustered data in Cartesian coordinates are visualized in a top view of the observation area, where each point represents the location of a grid cell and its colour the corresponding cluster, which contains the time series in that location. The centre of each grid cell in spherical coordinates is mapped back to the mean Cartesian xy-coordinates, to visualize the clusters of the range time series in a comparable way. Each cluster is associated with its cluster centroid, the mean time series in elevation of all time series in the respective cluster. In this way, cluster centroids are visualized as elevation time series independent of the 260 coordinate system that was used to generate them. This allows for a direct comparison of the cluster centroids between both time series extraction methods. We subsequently derive change processes visually from the entire clustered area. We establish which kind of deformation patterns can be distinguished and if they match, with what we expect in the respective areas (for example gradual erosion in the intertidal area).

265
We use the following criteria to compare the respective clustering and grid generation methods quantitatively: percentage of entire area clustered minimum and maximum within cluster variation percentage of correctly identified change in test area with bulldozer work The percentage of the area that is clustered differs depending on the algorithm. Especially DBSCAN sorts out points that 270 are too far away (i.e. too dissimilar) from others as noise. This will be measured over the entire observation area. The number of all complete time series counts as 100%.
Each cluster has a mean centroid time series and all other time series deviate from that to a certain degree. We calculate the average standard deviation over the entire month per cluster and report on the minimum and maximum value out of all realized clusters.

Test Area
To allow for a comparison of the clusters with a sort of ground truth, we selected a test area at the bottom of the footpath. In this area a pile of sand was accumulated by a bulldozer, after the entrance to the path was covered with lots of sand during the storm, as found by Anders et al. (2019). We chose two time stamps for illustration and show the elevation before the bulldozer activity on 12 January, after the bulldozer activity on 16 January and the difference between the elevations on these two days 280 in Figure 7. The area does not change significantly after this event. Within this test area we classify (manually) each point as 'stable' or 'with significant change' depending on a change in elevation of more than 5 cm (positive or negative). Then we

Results
The results are presented in two parts. First, we compare the different time series extraction methods. Then, we further analyse the clustering algorithms on the time series in elevation. We compare the results on the test area, where a bulldozer created 290 a pile of sand (as indicated in Figure 1) and in terms of percentage of data clustered, average standard deviation within each cluster and physical interpretation of clusters.

Elevation vs. Range Time Series
With the method to extract time series based on Cartesian coordinates we extract around 40000 grid cells that contain complete time series in elevation for the entire month. Using the method to extract time series based on spherical coordinates we obtain 295 around 47000 complete time series in range. In Figure 8 it can clearly be seen that the area covered by both time series extraction methods is different. The data in spherical coordinates covers a slightly longer part of the beach and more points at larger distance to the scanner. However, the clusters do not yield a clear partitioning of the beach according to change patterns. Since the spherical coordinate representation is more native to the scanner and the grid cells generally grow in size by variations in elevation, which is not captured well in spherical coordinates.

K-means
For the k-means algorithm, we choose to use k = 6 clusters. From our visual inspection this leads to good, usable results by partitioning into clusters that are small enough to capture geomorphic changes but not too large, which would make them less 305 informative. With the k-means algorithm, the entire observation area is divided into partitions, which change slightly depending on the random initialization. The standard deviation within each cluster is relatively high, and generally increases with the size of the cluster. Even the cluster with the smallest standard deviation over the entire month, shows a standard deviation of 0.7 m (cluster 5).
On the test area the k-means algorithm correctly classifies about 85% of points into 'stable', 'significant negative change', or 310 'significant positive change'. However, as can be seen in Figure 7, a part of the points with negative change are not identified.

Agglomerative Clustering
The agglomerative clustering algorithm is set up, as the k-means algorithm, to find six clusters. It produces results very similar to the clusters found with the k-means algorithm, as can be seen comparing Figures 9 and 10  2 and 3 from agglomerative clustering correspond roughly to the clusters 3 and 2 from k-means clustering. The ordering of clusters is according to size, so more time series are considered 'noisy' according to k-means, whereas agglomerative clustering assigns more of these time series to the gradually accreting cluster. All other clusters appear to be nearly identical and show similar spatial distributions as well as centroid shapes.
On the test area, the detection of negative and positive changes is more balanced and leads to an overall score of 88 % 320 correctly identified points. Agglomerative clustering clearly separates the path that was cleared by the bulldozer and identifies it as eroding. 19 https://doi.org/10.5194/esurf-2020-34 Preprint. Discussion started: 19 May 2020 c Author(s) 2020. CC BY 4.0 License.

DBSCAN
When we use the DBSCAN algorithm on the same data set, with minimum number of points N min = 30 and maximum distance ε = 0.05, a large part of the time series (55 %) is classified as noise, meaning that they are not very similar (i.e. not correlated, 325 since we use correlation as distance measure) to any of the other time series. However they roughly match the combined areas that are identified as stable and noisy by k-means (clusters 0 and 2). The remaining time series are clustered into six clusters.
The standard deviation within each cluster is generally lower than in the clusters generated with k-means (minimum standard deviations is 0.4 m) without considering the time series that are classified as noise.
The intertidal zone cannot be separated clearly from the 'noise' part of the observation area, nor can we distinguish the stable 330 path area or the upper part of the beach. But, two clusters represent areas, which are relatively stable throughout the month, except for a sudden peak in elevation on one day. These peaks are dominated by a van parking on the path on top of the dunes and people passing by and not caused by actual deformation in the observed area, as shown in Figure 11. On the test area the DBSCAN algorithm performs worse than both other algorithms. In total only 45% of points are correctly classified into 'stable', 'significant negative change', or 'significant positive change'. As stable points we count in this case all 335 points that are classified as noise, because only time series that show coherent change patterns are clustered by the DBSCAN algorithm. This matches with about 41% of points classified as stable in the ground truth data. However only 1.5% of the points with positive significant changes are correctly identified, they are mixed up with a large class of only slightly varying points.
We can see in Figure 7, that a significant part of the stable area is also included in the same cluster.

340
From the clustering using range time series, no clear change processes can be distinguished and the beach cannot be partitioned according to deformation patterns. Considering the clusters found by the k-means algorithm and agglomerative clustering on elevation time series, we can clearly distinguish between time series that represent erosion and gradual accretion, as well as a sudden jump in elevation, caused by bulldozer work. In Figure 12 we show the clusters and associated main cause for deformations. The cluster dominated 345 by erosion is close to the water line and roughly represent the inter-tidal zone of the beach and is likely caused by the effects of 21 https://doi.org/10.5194/esurf-2020-34 Preprint. Discussion started: 19 May 2020 c Author(s) 2020. CC BY 4.0 License.
tides and waves. The slowly accreting area is mostly at the upper beach, close to the dune foot and likely dominated by aeolian sand transport. The most obvious change process, is the sand removed from the entrances of the paths leading to the beach by bulldozer works (cluster 4) and accumulated in a pile of sand (cluster 5). The noisy cluster is spread out through the dune area and probably caused by moving vegetation.

Discussion
We successfully applied the presented methods on a data set from permanent laser scanning and demonstrated the identification of deformation processes from the resulting clusters. Here we discuss our results on time series extraction, distance measures, clustering methods and the choice of their respective input parameters and derivation of change processes.

355
We compared two different methods to extract time series from the PLS data set either in elevation or in range. The time series extraction in range, which is a more native way of using the data, is very sensitive to vertical structures in the data set, or points in the air in-between the observed surface and the scanner. After removing those points and using the median range per grid cell in spherical coordinates, the time series appear to be dominated by noisy fluctuations, which do not vary a lot depending on location. Clear change patterns can therefore not be distinguished with any of our algorithms and the distinction of areas 360 that follow a certain change pattern is not possible. A likely cause of this issue, is that most changes are observed in elevation (z-direction) and not in the direction of the range, which makes them less pronounced in spherical coordinates. An alternative approach could be the use of spherical coordinates for the generation of the grid cell, but extraction of time series in elevation instead of range.
In contrast, the time series extraction in Cartesian coordinates provides promising results. Some data is lost, due to low 365 point density in grid cells that are at large distances of the laser scanner. Besides that, the resulting clusters clearly follow recognizable deformation patterns and the clustering allows to separate regions according to these patterns.

Distance Measures
Possible distance measures for the use in time series clustering are analysed among others by Iglesias and Kastner (2013) and Liao (2005). We use Euclidean distance in combination with the k-means algorithm and agglomerative clustering for our 370 analysis. It has been shown by Keogh and Kasetty (2003) that especially for time series with high dimensions, alternative distance measures rarely outperform Euclidean distance. However, we have to note here, that Euclidean distance is affected by the so called 'curse of dimensionality', which causes a space of time series with many epochs (dimensions) to be difficult to cluster. For more details on this issue see Assent (2012) and Verleysen and François (2005). This could possibly lead to difficulties, when extending these methods to the use of longer time series, but does not appear to degrade results on our 375 current data set. We chose for the use of correlation distance with the DBSCAN algorithm, because correlation in principle represents a more intuitive way of comparing time series (see Figure 5). DBSCAN is based on identifying the clusters of high density, which in our case works better using correlation distance instead of Euclidean distance.
Neither of the the two distance measures analysed here can deal with gaps in the time series. They also do not allow to 380 identify identical elevation patterns that are shifted in time as similar. An alternative distance measure suitable to deal with these issues would be Dynamic Time Warping (DTW), which accounts for similarity in patterns even though they might be shifted in time (Keogh and Ratanamahatana (2005)). An interpolation method to fill gaps in elevation over short time spans based on surrounding data or a feature based clustering method could be other alternatives.

385
The use of k-means clustering on elevation time series from the same data set was demonstrated by Lindenbergh et al. (2019) and has shown promising first result. We follow the same approach and, as a comparison, use agglomerative clustering, with the same optimization criterion, distance metric and input parameter. As expected the results are similar, although agglomerative clustering does not depend on random initialization. It therefore delivers the same result for every run, which is an advantage.
Considering our previously defined criteria:

390
a majority of the observation area is separated into distinct regions, each cluster shows a change pattern that can be associated with a geomorphic deformation process, and time series contained in each cluster roughly follow the mean change pattern, both algorithms are suitable and the differences in the resulting clusters are negligible for our specific data set.
However, the computational effort needed to loop through all possible combinations of merging clusters for agglomerative 395 clustering is considerably higher. Of the three algorithms that were used in this study, agglomerative clustering is the only one that regularly ran into memory errors. This is a disadvantage considering the possible extension of our method to a data set with longer time series.
One of the disadvantages of the k-means algorithm and our configuration of agglomerative clustering, is that the number of clusters has to be defined in advance. Our choice of k = 6 clusters yields promising results, but remains somewhat arbitrary, 400 especially without prior knowledge of the data set. We have considered two different methods to determine a suitable value for k: analysis of the overall sum of variances for different values of k and so-called cluster balance following the approach of Jung et al. (2003). Neither of them resolved the problem satisfactorily and we cannot make a generalized recommendation for the choice of k at this point.
To avoid this issue we also compare both approaches with the use of the DBSCAN algorithm. It is especially suitable to 405 distinguish anomalies and unexpected patterns in data as demonstrated by Çelik et al. (2011) using temperature time series.
To decide, which values are most suitable for the two input parameters of the DBSCAN algorithms we plot the percentage of clustered points and the number of clusters depending on both parameters (see Figure 13). However, this did not lead to a clear indication of an 'optimal' set of parameters. After the trade-off analysis between the number of points in clusters and the number of clusters (not too high, so that the clusters become very small and not too low so that we generate only one big 410 cluster) we chose ε = 0.05 and N min = 30 by visually inspecting the resulting clusters.
An alternative clustering approach for time series based on fuzzy C-means is proposed by Coppi et al. (2010). They develop a method to balance the clustering based on the pattern of time series while keeping an approximate spatial homogeneity of the clusters. This approach was successfully applied to time series from socio-economic indicators and could be adapted for our purpose. It could potentially improve detection of features like sand bars, or bulldozer work, but not distinguish moving 415 vegetation in the dunes as our current approach does.
A similar approach would be to use our clustering results and identified change patterns as input to the region-growing approach of Anders et al. (2020). In this way we could combine advantages of both approaches by making the identification of the corresponding regions for each distinct deformation pattern more exact, without having to define possible deformation patterns in advance.

Derivation of Change Processes
As shown in Figure 12, we identified change processes from the clusters generated by k-means and agglomerative clustering.
Each centroid representing the mean time series of its cluster shows a distinct change pattern (see Figures 9 and 10), which allows to conclude on a predominant deformation. By associating the centroids with the location and spatial spread of the clusters, we can derive the main cause for this deformation. In some cases extra information, or an external source of validation 425 data would be useful to verify the origin of the process. This will be taken into account for future studies. However, the steep rise of the centroid from cluster 5 allows to conclude that the cause of this sudden accretion is not natural. The information found by Anders et al. (2019) for the research on their study, confirms the coinciding bulldozer works.
The DBSCAN algorithm successfully identifies parts of the beach that are dominated by a prominent peak in the time series (caused by a van and a small group of people). Out of the three algorithms that we compare, it is most sensitive to these outliers 430 in the form of people or temporary objects in the data. It was not our goal for this study, to detect people or objects on the beach, but this ability could be a useful application of the DBSCAN algorithm to filter the data for outliers in a pre-processing step.
We compared three different clustering algorithms (k-means, agglomerative clustering and DBSCAN) on a subset of a large 435 time series data set from permanent laser scanning on a sandy urban beach. We successfully separated the observed beach and dune area according to their deformation patterns. Each cluster, described by the mean time series, is associated with a specific process (such as bulldozer work, tidal erosion) or surface property (for example moving vegetation cover).
The most promising results are found using k-means and agglomerative clustering, which both correctly classify between 85 and 88 % of time series in our test area. However, they both need the input of the number of clusters we are looking for 440 and agglomerative clustering is computationally expensive. DBSCAN turned out to be more suitable for the identification of outliers or unnatural occurring changes in elevation due to temporary objects or people in the observed area.
Our key findings are summarized as follows: 1. Both k-means and agglomerative clustering fulfil our criteria for a suitable method to cluster time series from permanent laser scanning.

445
2. Predominant deformation patterns of sandy beaches are detected automatically and without prior knowledge using these methods.
3. Change processes on sandy beaches, which are associated with a specific region and time span, are detected in a spatiotemporal data set from permanent laser scanning with the presented methods.
Our results demonstrate a successful method to mine a spatio-temporal data set from permanent laser scanning for pre-450 dominant change patterns. The method is suitable for the application in an automated processing chain to derive deformation patterns and regions of interest from a large spatio-temporal data set. It allows such a data set to be partitioned in space and time according to specific research questions into phenomena, such as for example the interaction of human activities and natural sand transport during storms, recovery periods after a storm event or the formation of sand banks. The presented methods enable the use of an extensive time series data set from permanent laser scanning to support the research on long-term and small 455 scale processes on sandy beaches and improve analysis and modelling of these processes. In this way we expect to contribute to an improved understanding and managing of these vulnerable coastal areas.