Daniel Borbor
email: d_borbor@ciise.concordia.ca

Lingyu Wang
email: wang@ciise.concordia.ca

Sushil Jajodia
email: jajodia@gmu.edu

Anoop Singhal
email: anoop.singhal@nist.gov

Diversifying Network Services under Cost Constraints for Better Resilience against Unknown Attacks

Diversity as a security mechanism has received revived interest recently due to its potential for improving the resilience of software and networks against unknown attacks. Recent work show diversity can be modeled and quantified as a security metric at the network level. However, such an effort does not directly provide a solution for improving the network diversity, and existing network hardening approaches are largely limited to handling previously known vulnerabilities by disabling existing services. In this paper, we take the first step towards an automated approach to diversifying network services under various cost constraints in order to improve the network's resilience against unknown attacks. Specifically, we provide a model of network services and formulate the diversification requirements as an optimization problem. We devise optimization and heuristic algorithms for efficiently diversifying relatively large networks under different cost constraints. We also evaluate our approach through simulations.

Introduction

Many critical infrastructures, governmental and military organizations, and enterprises have become increasingly dependent on networked computer systems today. Such mission critical computer networks must be protected against not only known attacks, but also potential zero day attacks exploiting unknown vulnerabilities. However, while traditional solutions, such as firewalls, vulnerability scanners, and IDSs, are relatively successful in dealing with known attacks, they are less effective against zero day attacks.

To this end, diversity has previously been considered for a security mechanism for hardening software systems against unknown vulnerabilities, and it has received a revived interest recently due to its potential for improving networks' resilience against known attacks. In particular, a recent work shows diversity can be modeled and quantified as a security metric at the network level [START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Zhang | Network diversity: A security metric for evaluating the resilience of networks against zero-day attacks[END_REF]. However, the work does not directly provide a systematic solution for improving the network diversity under given cost constraints, which can be a challenging task for large and complex networks. On the other hand, existing efforts on network hardening (a detailed review of related work will be given later in Section 2) are largely limited to handling previously known vulnerabilities by disabling existing services.

In this paper, we propose an automated approach to diversifying network services under various cost constraints in order to improve the network's resilience against unknown attacks. Specifically, we devise a model of network services and their different instances by extending the resource graph model; such a model allows us to formulate the diversification requirements and cost constraints as an optimization problem; we apply optimization techniques to solve the formulated problems, and design heuristic algorithms to more efficiently handle larger networks. We evaluate our approach through simulations in order to study the effect of optimization parameters on accuracy and running time, and the effectiveness of optimization for different types of networks. In summary, the main contribution of this paper is twofold:

-To the best of our knowledge, this is the first effort on formulating the problem of network service diversification for improving the resilience of networks, which enables the application of existing optimization techniques and also provides a practical application for existing diversity metrics [START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Zhang | Network diversity: A security metric for evaluating the resilience of networks against zero-day attacks[END_REF]. -As evidenced by the simulation results, the optimization and heuristic algorithms provide a relatively accurate and efficient solution for diversifying network services while considering various cost constraints. By focusing on zero day attacks, our work provides a complementary solution to existing network hardening approaches that focus on fixing known vulnerabilities.

The remainder of this paper is organized as follows: The rest of this section first builds the motivation through a running example. Section 2 reviews related work. In Section 3, we present the model and formulate the optimization problem, and in Section 4 we discuss the methodology and show case studies. Section 5 shows simulation results and Section 6 concludes the paper.

Motivating Example

We present a motivating example to demonstrate that diversifying network services can be a tedious and error-prone task if done manually, even if the considered network is of a small size. Figure 1 shows a hypothetical network, which is roughly based on the virtual penetration lab described in [START_REF]Penetration testing virtual labs[END_REF]. Despite its relatively small scale, it mimics a typical enterprise network, e.g., with DMZ, Web server behind firewall accessible from public Internet, and a private management network protected by another firewall. Specifically, the network consists of four hosts running one or more services allowing accesses from other hosts. We assume the two firewalls and other host-based mechanisms (e.g., personal firewalls or iptables) together enforce the connectivity described inside the connectivity table shown in the figure. We consider attackers on external hosts (represented as h0) attempting to compromise the database server (h4), and we assume the network is secured against known vulnerabilities (we exclude exploits and conditions that involve the firewalls).

To measure the network's resilience against unknown zero day attacks, we consider the k-zero day safety metric [START_REF] Wang | k-zero day safety: A network security metric for measuring the risk of unknown vulnerabilities. Dependable and Secure Computing[END_REF] (which will be referred to as k0d from now on for simplicity), which basically counts how many distinct zero day vulnerabilities must exist and be exploited before an attacker may reach the goal. For simplicity, although the attacker may follow many paths to compromise h4, here we only consider the Web servers as the initial targets. We can observe that there must exist at least two distinct zero-day vulnerabilities, one for the Apache servers (h2, h3, and h4) and one for the IIS server1 (h1), and the attacker must exploit both in order to compromise h4. Finally, we assume the administrator has the option of replacing those Web servers with either an NGINX 1.9 or a Litespeed 5.0.14 Web server and each replacement will incur a given installation/maintenance cost (we will discuss the cost model in more details later in Section 3). Based on above assumptions, consider the following scenarios:

-Scenario 1: The administrator aims to render the network as resilient as possible to zero-day attacks (which means to maximize the aforementioned k0d metric). -Scenario 2: He/she aims at the same goal as in above Scenario 1, but under the constraint that the overall diversification cost must be less than a given budget. -Scenario 3: He/she aims at the same goal as in above Scenario 2, but under an additional constraint that at most two Web servers may be replaced. -Scenario 4: He/she aims at the same goal as in above Scenario 3, but under an additional constraint that replacing the IIS web server in the DMZ (h4) should be given a higher priority.

Clearly, many more use cases may exist in practice than those listed above, and the solution may not always be straightforward even for such a small network. For example, while the administrator can easily increase the k0d metric value to 4 under Scenario 1 (by having four different Web servers), the optimal solution in other scenarios will critically depend on the specific cost constraints and given budgets. Considering that the attacker may also follow other paths to attack (e.g., starting with SMTP, instead of Web, on h1), the problem becomes even more complicated. This shows the need for an automated approach, which will be the subject matter of the remainder of this paper.

In general, the security of networks may be qualitatively modeled using attack trees [START_REF] Dewri | Optimal security hardening using multi-objective optimization on attack tree models of networks[END_REF][START_REF] Dewri | Optimal security hardening on attack tree models of networks: a cost-benefit analysis[END_REF][START_REF] Ray | Using attack trees to identify malicious attacks from authorized insiders[END_REF] or attack graphs [START_REF] Paul Ammann | Scalable, graph-based network vulnerability analysis[END_REF][START_REF] Sheyner | Automated generation and analysis of attack graphs[END_REF]. A majority of existing quantitative models of network security focus on known attacks [START_REF] Wang | Measuring the overall security of network configurations using attack graphs[END_REF][START_REF] Albanese | Time-efficient and cost-effective network hardening using attack graphs[END_REF], while few works have tackled zero day attacks [START_REF] Wang | k-zero day safety: Measuring the security risk of networks against unknown attacks[END_REF][START_REF] Wang | k-zero day safety: A network security metric for measuring the risk of unknown vulnerabilities. Dependable and Secure Computing[END_REF][START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Zhang | Network diversity: A security metric for evaluating the resilience of networks against zero-day attacks[END_REF] which are usually considered unmeasurable due to the uncertainties involved [START_REF] Mchugh | Quality of protection: measuring the unmeasurable?[END_REF].

Early works on network hardening typically rely on qualitative models while improving the security of a network [START_REF] Sheyner | Automated generation and analysis of attack graphs[END_REF][START_REF] Wang | Minimum-cost network hardening using attack graphs[END_REF][START_REF] Wang | Network Hardening: An Automated Approach to Improving Network Security[END_REF]. Those works secure a network by breaking all the attack paths that an attacker can follow to compromise an asset, either in the middle of the paths or at the beginning (disabling initial conditions). Also, those works do not consider the implications when dealing with budget constraints nor include cost assignments, and tend to leave that as a separate task for the network administrators. While more recent works [START_REF] Albanese | Time-efficient and cost-effective network hardening using attack graphs[END_REF][START_REF] Yigit | Cost-aware network hardening with limited budget using compact attack graphs[END_REF] generally provide a cost model to deal with budget constraints, one of the first attempts to systematically address this issue is by Gupta et al. [START_REF] Gupta | Matching information security vulnerabilities to organizational security profiles: a genetic algorithm approach[END_REF]. The authors employed genetic algorithms to solve the problem of choosing the best set of security hardening options while reducing costs. Dewri et a. [START_REF] Dewri | Optimal security hardening using multi-objective optimization on attack tree models of networks[END_REF] build on top of Gupta's work to address the network hardening problem using a more systematic approach. They start by analyzing the problem as a single objective optimization problem and then consider multiple objectives at the same time. Their work consider the damage of compromising any node in the cost model in order to determine the most cost-effective hardening solution. Later on, in [START_REF] Dewri | Optimal security hardening on attack tree models of networks: a cost-benefit analysis[END_REF] and in [START_REF] Wang | Exploring attack graph for costbenefit security hardening: A probabilistic approach[END_REF], the authors extrapolate the network hardening optimization problem as vulnerability analysis with cost/benefit assessment, and risk assessment respectively. In [START_REF] Poolsappasit | Dynamic security risk management using bayesian attack graphs. Dependable and Secure Computing[END_REF] Poolsappasit et al. extend Dewri's model to also take into account dynamic conditions (conditions that may change or emerge while the model is running) by using Bayesian attack graphs in order to consider the likelihood of an attack. Unlike our work, most existing work on network hardening are limited to known vulnerabilities and focus on disabling existing services.

There exist a rich literature on employing diversity for security purposes. The idea of using design diversity for tolerating faults has been investigated for a long time, such as the N-version programming approach [START_REF] Avizienis | On the implementation of n-version programming for software fault tolerance during execution[END_REF], and similar ideas have been employed for preventing security attacks, such as the N-Variant system [START_REF] Cox | N-variant systems: a secretless framework for security through diversity[END_REF], and the behavioral distance approach [START_REF] Gao | Behavioral distance measurement using hidden markov models[END_REF]. In addition to design diversity and generated diversity, recent work employ opportunistic diversity which already exists among different software systems. For example, the practicality of employing OS diversity for intrusion tolerance is evaluated in [START_REF] Garcia | Os diversity for intrusion tolerance: Myth or reality? In Dependable Systems & Networks (DSN)[END_REF]. More recently, the authors in [START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Zhang | Network diversity: A security metric for evaluating the resilience of networks against zero-day attacks[END_REF] adapted biodiversity metrics to networks and lift the diversity metrics to the network level. While those works on diversity provide motivation and useful models, they do not directly provide a systematic solution for improving diversity, which is the topic of this paper.

Model

We first introduce the extended resource graph model to capture network services and their relationships, then we present the diversity control and cost model, followed by problem formulation.

Extended Resource Graph

The first challenge is to model different resources, such as services (e.g., Web servers) that can be remotely accessed over the network, different instances of each resource (e.g., Apache and IIS), and the causal relationships existing among resources (e.g., a host is only reachable after an attacker gains a privilege to another host). For this purpose, we will extend the concept of resource graph introduced in [START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Zhang | Network diversity: A security metric for evaluating the resilience of networks against zero-day attacks[END_REF], which is syntactically equivalent to attack graphs, but models network resources instead of known vulnerabilities as in the latter.

Specifically, we will define an extended resource graph by introducing the notion of Service Instance to indicate which instance (e.g., Apache) of a particular service (e.g., Web server) is being used on a host. Like the original resource graph, we only consider services that can be remotely accessed. The extended resource graph of the running example is shown in Figure 2 and detailed below.

Fig. 2. The example network's resource graph

In Figure 2, each pair shown in plaintext is a security-related condition (e.g., connectivity source, destination or privilege privilege, host ). Each exploit node (oval) is a tuple that consists of a service running on a destination host, the source host, and the destination host (e.g., the tuple http, 1, 2 indicates a potential zero day vulnerability in the http service on host 2, which is exploitable from host 1). The small one-column table beside each exploit indicates the current service instance using a highlighted integer (e.g., 1 means Apache and 2 means IIS) and other potential instances in lighter text. The self-explanatory edges point from pre-conditions to an exploit (e.g., from 1, 2 and http, 2 to http, 1, 2 ), and from the exploit to its post-conditions (e.g., from http, 1, 2 to user, 2 ).

A design choice here is whether to associate the service instance concept with a condition indicating the service (e.g., http, 2 ) or the corresponding exploits (e.g., http, 1, 2 ). While it is more straightforward to have the service instance defined as the property of a condition, which can then be inherited by the corresponding exploits, we have opted to define this property as a label for the exploit nodes in the graph, be-cause this will make it easier to check the number of distinct services along a path, as we will see later. One complication then is that we must ensure all exploits with the same service and destination host (e.g., http, 1, 2 and http, 3, 2 ) to be associated with the same service instance. Definitions 1 and 2 formally introduce these concepts.

Definition 1 (Service Pool and Service Instance). Denote S the set of all services and Z the set of integers, for each service s ∈ S, the function sp(.) : S → Z gives the service pool of s which represent all available instances of that service.

, h d | h s ∈ H, h d ∈ H, s ∈ serv(h d )},
and R r ⊆ C ×E and R i ⊆ E ×C be the collection of pre and post-conditions in C. We call the labeled directed graph, G(E ∪ C, R r ∪ R i ), v the extended resource graph.

Diversity control and cost model

We employ the notion of diversity control as a model for diversifying one or more services in the resource graph. Since we represent the service instance using integers, it will be straightforward to regard each pair of service and destination host on which the service is running as an optimization variable, and formulate diversity control vectors using those variables as follows. We note that the number of optimization variables present in a network will depend on the number of conditions indicating services, instead of the number of exploits (since many exploits may share the same service instance, and hence the optimization variable). Since we only consider remotely accessible services in the extended resource graph model, we would expect in practice the number of optimization variables to grow linearly in the size of network (i.e., the number of hosts). Definition 3 (Optimization Variable and Diversity Control). Given an extended resource graph G, v , ∀e ∈ E, v(e) is an optimization variable. A diversity control vector is the integer valued vector V = (v(e 1 ), v(e 2 ), ..., v(e |E| ).

Changing the value of an optimization variable has an associated diversification cost and the collection of such costs is given in a diversity cost matrix in a self-explanatory manner. We assume the values of cost are assigned by security experts or network administrators. Like in most existing works (e.g., [START_REF] Dewri | Optimal security hardening using multi-objective optimization on attack tree models of networks[END_REF]), we believe an administrator can estimate the diversification costs based on monetary, temporal, and scalability criteria like i) installation cost, ii) operation cost, iii) training cost, iv) system downtime cost and, v) incompatibility cost. We define the diversity cost, diversity cost matrix, and the total diversity cost. Definition 4 (Diversification Cost and Diversity Cost Matrix). Given s ∈ S and sp(s), the cost to diversify a service by changing its service instance to another inside the service pool is called the diversification cost. The collection of all the costs constraints associated with changing services in S are given as a diversity cost matrix DCM in which the element at i th row and j th column indicates the diversification cost of changing the i th service instance to be the j th service instance. Let v s (e i ) be the service associated with the optimization variable v(e i ) and V 0 the initial service instance values for each of the exploits in the network. The total diversification cost, Q d , given by the diversity vector V is obtained by

Q d = |E| i=1 DCM vs(ei) (V 0 (i), V (i))
We note that the above definition of diversification cost between each pair of service instances has some advantages. For example, in practice we can easily imagine cases where the cost is not symmetric, i.e., changing one service instance to another (e.g., from Apache to IIS) carries a cost that is not necessarily the same as the cost of changing it back (from IIS to Apache). Our approach of using a collection of two-dimensional matrices allows us to account for cases like this. Also, the concept can be used to specify many different types of cost constraints, which we will examine in the coming section. For example, an administrator who wants to restrict the total cost to diversify all servers running the http service can do so by simply formulating the cost as the addition of all the optimization variables corresponding to http.

Problem formulation

As demonstrated in Section 1.1, the k0d metric is defined as the minimum number of distinct resources on a single path in the resource graph [START_REF] Wang | k-zero day safety: A network security metric for measuring the risk of unknown vulnerabilities. Dependable and Secure Computing[END_REF]. For example, a closer look at Figure 2 shows that the k0d value for our example network is 1. That is, an attacker needs only one zero-day vulnerability (in http service instance 1) to compromise this network. The dashed line in Figure 2 depicts the shortest path that provides this metric value.

The k0d value can be increased by changing the service instances as long as we respect the available budget of cost. For example, consider a total budget of 78 units, and assume the costs to diversify the http service from service instance 1 to 2, 3 or 4 be 78, 12, and 34 units, respectively. We can see that changing http, 2, 3 from instance 1 to 2 would respect the budget, as well as increasing the k0d value of the network to be 2. We may also see that this is not the optimal solution, since we could also replace http, 2, 3 and http, 3, 4 with instances 3 and 4, respectively, increasing k0d to 3 and still respecting the budget. In the following, we formally formulate this as an optimization problem.

Problem 1 (k0d Optimization Problem). Given an extended resource graph G, v , find a diversity control vector V which maximizes min(k0d( G(V ), v )) subject to the constraint Q ≤ B, where B is the availble budget and Q is the total diversification cost as given in Definition 4.

Since our problem formulation is based on an extended version of the resource graph, which is syntactically equivalent to attack graphs, many existing tools developed for the latter (e.g., the tool in [START_REF] Jajodia | Topological analysis of network attack vulnerability[END_REF] has seen many real applications to enterprise networks) may be easily extended to generate the extended resource graphs we need as inputs. Additionally, our problem formulation assumes a very general budget B and cost Q, which allows us to account for different types of budgets and cost constraints that an administrator might encounter in practice, as will be explained in the following section.

Methodology

This section details the optimization and heuristic algorithms used for solving the formulated diversification problem and describes a few case studies.

Genetic Algorithm Optimization

The genetic algorithm (GA) is a simple and robust search method and optimization technique inspired by the mechanisms of natural selection. We employ GA for our automated optimization approach because it requires little information to search effectively in a large search space in contrast to other optimization methods (e.g., the mixed integer programming [START_REF] Md Azamathulla | Comparison between genetic algorithm and linear programming approach for real time operation[END_REF]). It also provides a simple way to encode candidate solutions to the problem [START_REF] Kalyanmoy | An efficient constraint handling method for genetic algorithms[END_REF]. While inspired by [START_REF] Dewri | Optimal security hardening using multi-objective optimization on attack tree models of networks[END_REF], we focus on service diversification and not on disabling services.

The extended resource graph is the input to our automated optimization algorithm where the function to be optimized (fitness function) is k0d defined on the resource graph (later we will discuss cases where directly evaluating k0d is computationally infeasible). One important point to consider when optimizing the k0d function on the extended resource graph is that, for each generation of the GA, the graph's labels will dynamically change. This in turn will change the value of k0d, since the shortest path may have changed with each successive generation of the GA. Our optimization tool takes this into consideration. We also note one limitation here is that the optimization does not provide a priority if there are more than one shortest path that provide the optimized k0d since the optimization only aims at maximizing the minimum k0d.

The constraints are defined as a set of inequalities in the form of q ≤ b, where q represents one or more constraint conditions and b represents one or more budgets. These constraint conditions can be overall constraints (e.g., the total diversity cost Q d ) or specific constraints to address certain requirements or priorities while diversifying services (e.g., the cost to diversify http services should be less than 80% of the cost to diversify ssh). Those constraints are specified using the diversity control matrix.

The number of independent variables used by the GA (genes) are the optimization variables given by the extended resource graph. For our particular network hardening problem, the GA will be dealing with integer variables representing the selection of the service instances. Because v(e) is defined as an integer, the optimization variables need to be given a minimum value and a maximum value. This range is determined by the number of instances provided in the service pool of each service. The initial service instance for each of the services is given by the extended resource graph while the final diversity control vector V is obtained after running the GA.

The population size that we defined for our tool was set to be at least the value of optimization variables (more details will be provided in the coming section). This way we ensure the individuals in each population span the search space. We ensure the population diversity by testing with different settings in genetic operations (like crossover and mutation). In the following, we discuss several test cases to demonstrate how the optimization works under different types of constraints. For all the test cases, we have used the following algorithm parameters: population size = 100, number of generations = 150, crossover probability = 0.8, and mutation probability = 0.2.

Test case A: Q d ≤ 124 units with individual constraints per service. We start with the simple case of one overall budget constraint (Q d ≤ 124). The solution provided by the GA is V = [3, 2, 1, 4, 1, 1, 1] (represented by label column a in Figure 3). The associated costs for V (1), V (2), and V (4) are 12, 78, and 34, respectively, and the test network's k0d metric becomes 4 while keeping Q d within the budget (Q d ≤ 124). On the other hand, if we assign individual budgets per services, while maintaining overall budget Q d ≤ 124, the optimization results will be quite different. In this case, assume the budget to diversify the http services cannot exceed 100 units (q http ≤ 100); for f tp, it cannot exceed 3 units (q f tp ≤ 3); for ssh, it cannot exceed 39 units (q ssh ≤ 39); finally, for smtp, it cannot exceed 50 units (q smtp ≤ 50). The solution provided by the GA is a V vector where V (1) = 2 and V (2) = 3, with a cost of 78 and 12 units, respectively. The value of the k0d metric rises to 3 with Q d = 90. This total diversification cost satisfies both the overal budget constraint and each of the individual constraints per service.

From this test case, we can see that even with the minimum requeried budget to maximize the k0d metric, additional budget constraints might not allow to achieve the maximum k0d possible. We can see the result of running the GA for this test case in label column b in Figure 3.

Test case B: Q d ≤ 124 units while q http + q ssh ≤ 100. While test case 1 shows how individual cost constraints can affect the k0d metric optimization, in practice not all services may be of concern and some may have negligible cost. This test case models such a scenario by assigning a combined budget restriction for only the http and ssh services, i.e., the cost incurred by diversifying these two services should not exceed 100 units.

The solution provided by the GA is V = [3, 4, 3, 1, 1, 3, 2] (lable column a in Figure 4). Since V (1) to V (3) deal with the http service, we can see that the total incurred cost for http is q http =12+34+12=58 units. Because V (6) and V (7) are the only optimization variables that deal with the f tp and ssh services respectively, we can see that q f tp = 8, and q ssh = 40. The value of the k0d metric rises from 1 to 3 by incurring a total cost of Q d = 106 units. The combined http/ssh budget constraint of 100 units is also satisfied since q http + q ssh = 98 units. Test case C: Q d ≤ 124 units while q http ≤ 0.8 • q ssh . This final case deals with scenarios where some services might have a higher priority over others. The constraint in this test case is that the total incurred cost while diversifying the http service should not exceed 80% of what is incurred by diversifying the ssh service.

The solution provided by the GA is V =[3,1,3,1,1,1,4] (see column b in Figure 4). Here V (1) and V (3) have changed from service instance 1 to 3, while V (7) have changed from service instance 1 to 4. The incurred cost for the http service is q http = 12+12=24 units and for the ssh service is q ssh = 34 units. While the value of the k0d metric only rises from 1 to 2, the budget constraints are satisfied.

As seen from the above test cases, our model and problem formulation makes it relatively easy to apply any standard optimization techniques, such as the GA, to optimize the k0d metric through diversity while dealing with different budget constraints.

Fig. 5. A Heuristic algorithm for calculating m-shortests paths

The algorithm starts by topologically sorting the graph (line 1) and proceeds to go through each one of the nodes on the resource graph collection of attack paths, as set of exploits σ(), that reach that particular node. The main loop cycles through each unprocessed node. If a node is an initial conditions, the algorithm assumes that the node itself is the only path to it and it marks it as processed (lines 6-8). For each exploit e, all of its preconditions are placed in a set (line 10). The collection of attack paths α(e) is constructed from the attack paths of those preconditions (lines 10 and 11). In a similar way, σ (ov(e)) is constructed with the function ov() which, aside of using the exploits includes value of element of the diversity control vector that supervises that exploit.

If there are more than m paths to that node, the algorithm will use the function U nique to first look for unique combinations of service and service instance in α (ov(e)). Then, the algorithm creates a dictionary structure where the key is a path from α(e) and the value is the number of unique service/service instance combinations given by each one of the respective paths in α (ov(e)). The function ShortestM () selects the top m keys whose values are the smallest and returns the m paths with the minimum number of distinct combination of services and service instances (line 13). If there are less than m paths, it will return all of the paths (line 15). After this, it marks the node as processed (line 16). The process is similar when going through each one of the intermediate conditions (lines [START_REF] Sheyner | Automated generation and analysis of attack graphs[END_REF][START_REF] Wang | Network Hardening: An Automated Approach to Improving Network Security[END_REF][START_REF] Wang | k-zero day safety: A network security metric for measuring the risk of unknown vulnerabilities. Dependable and Secure Computing[END_REF][START_REF] Wang | k-zero day safety: Measuring the security risk of networks against unknown attacks[END_REF][START_REF] Wang | Minimum-cost network hardening using attack graphs[END_REF][START_REF] Wang | Measuring the overall security of network configurations using attack graphs[END_REF][START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Wang | Exploring attack graph for costbenefit security hardening: A probabilistic approach[END_REF]. Finally, the algorithm returns the collection of m paths that can reach the goal condition c g . It is worth noting that the algorithm does not make any distinction in whether or not a particular path has a higher priority over another when they share the same number of unique service/service instance combinations.

Simulations

In this section, we show simulation results. All simulations are performed using a computer equipped with a 3.0 GHz CPU and 8GB RAM in the Python 2.7.10 environment under Ubuntu 12.04 LTS and MATLAB 2015a's GA toolbox. To generate a large number of resource graphs for simulations, we first construct a small number of seed graphs based on real networks and then generate larger graphs from those seed graphs by injecting new hosts and assigning resources in a random but realistic fashion (e.g., the number of pre-conditions of each exploit is varied within a small range since real world exploits usually have a constant number of pre-conditions). The resource graphs were used as the input for the optimization toolbox where the objective function is to maximize the minimum k0d value subject to budget constraints. In all the simulations, we employ the heuristic algorithm described in section 4.2. Figure 6 shows that the processing time increases almost linearly as we increase the number of optimization variables or the parameter m of the heuristic algorithm. The results show that the algorithm is relatively scalable with a linear processing time. On the other hand, the accuracy of the results is also an important issue to be considered. Here the accuracy refers to the approximation ratio between the result obtained using the heuristic algorithm and that of the brute force algorithm (i.e., simply enumerating and searching all the paths while assuming all services and service instances are different). For the simulations depicted in Figure 7, we settled for 50 iterations per graph per m-paths. The diversity control vector provided by the GA is used to calculate the accuracy. From the results, we can see that when m is greater or equal to 4 the approximation ratio reaches an acceptable level. For the following simulations, we have settled with an m value of 6 and 100 generations.

Our simulations also showed that (detailed simulation results are omitted here due to page limitations), when no budget constraints are in effect, using the GA with a crossover probability of 80%, a mutation rate of 20%, and setting the number of generations to 50 will be sufficient to obtain good results. However, this is no longer the case when dealing with budget constraints. We have noticed that, by decreasing the crossover probability (and consequently increasing the mutation rate), we can reach a viable solution with less generations. We have therefore settled with a crossover probability of 40% which provides us with a fast (with less generations) way to converge to viable solutions. Additionally, our experiences also show that, when dealing with a diversity control vector (also known as a chromosome in the GA) of less than 100 variables (genes in the GA), the population size could be equal to the amount of variables in the diversity control vector; when dealing with a bigger number, the population size should be at least twice the amount of variables. The number of minimum service instances available

Average gain

Fig. 9. The average gain vs the number of nodes.

Figure 8 shows the results when the diversity control vector has different numbers of sevice instances to take from (i.e., different sizes of the service pools). In this simulation, we have picked graphs with a relative high difference in the length of the shortest path before and after all services are diversified using the algorithm (the maximum k0d value is 16 and the minimum 3). We can see an increasing gain in the k0d value after optimization, when more service instances are available. However, this trend begins to stall after a certain number [START_REF] Mchugh | Quality of protection: measuring the unmeasurable?[END_REF]. From this observation it can be inferred that the number of available service instances will affect the difference between the maximum k0d value possible and the minimum k0d, but such an effect also depends on the size of the network (or the extended resource graph) and increasing the number of available service instances does not always help.

In Figure 9, we analyze the average gain in the optimized results for different sizes of graphs. In this figure, we can see that we have a good enough gain for graphs with a relatively high amount of nodes. As expected, as we increase the size of the graphs, the gain will decrease if we keep the same optimization parameters. For those simulations, we have used a population size of 300, 50 generations, and a crossover fraction of 50%. It is interesting to note that the decrease in gain is very close to being linear.

Figure 10 and Figure 11 show the optimization results on different shapes of resource graphs. While it may be difficult to exactly define the depth of a resource graph, we have relied on the relative distance, i.e., the difference of the shortest path before and after all services are diversified. There is a relative linear increase in the gain as we increase the relative distance in the shortest path. While this does not provide an accurate description of the graph's shape, it does provide an idea of how much our algorithm can increase the minimum k0d for graphs with different depths, as shown in Figure 10.

The k0d difference on the shortest path The number of directly reachable exploits Finally, in Figure 11, we can see the effect of the network's degree of exposure, which is defined as the number of exploits that are directly reachable by the attacker from the external host h0. As we increase the degree of exposure, the gain in optimization decreases in almost a linear way. That is, there will less room for diversification if the network is more exposed.

Conclusions

In this paper, we have formulated service diversity as an optimization problem and proposed an automated diversity-based network hardening approach against zero-day attacks. This automated approach used a heuristic algorithm that helped to manage the complexity of computing the k0d value as well as limiting the time for optimization to an acceptable level. We have shown some sample cost constraints while our model and problem formulation would allow for other practical scenarios to be specified and optimized. We have tested the scalability and accuracy of the proposed algorithms through simulation results, and we have also discussed how the gain in the k0d value will be affected by the number of available service instances in the service pools and different sizes and shapes of the resource graphs.

We discuss several aspects of the proposed automated optimization technique where additional improvements and evaluations can be done.

-While this paper focuses on diversifying services, a natural future step is to integrate this approach with other network hardening options, such as addition or removal of services, or relocating hosts or services (e.g., firewalls). -This study has relied on a simplified model by assuming all service instances to be completely different from each another and all service instances are equally likely to be exploited. A possible future research direction would be to model the degree of difference (or similarity) between the different types of service instances. -We have assumed an abstract cost model in this paper and an important direction is to elaborate the model from different aspects of potential cost for diversifying network resources. -We will also consider other optimization algorithms, in addition to GA, to compare and potentially use them in hybrid optimization schemes when searching for more efficient and effective solutions to our problem. -This study relies on a static network configuration. A future research direction would be to consider a dynamic network model in which both attackers and defenders may cause changes in the network.

Fig. 1 .

 1 Fig. 1. Example network.

Fig. 3 .

 3 Fig. 3. Test case A: general and individual budget constraints.

Fig. 4 .

 4 Fig. 4. Test case B and test case C.Test case C: Q d ≤ 124 units while q http ≤ 0.8 • q ssh . This final case deals with scenarios where some services might have a higher priority over others. The constraint in this test case is that the total incurred cost while diversifying the http service should not exceed 80% of what is incurred by diversifying the ssh service.The solution provided by the GA is V =[3,1,3,1,1,1,4] (see column b in Figure4). Here V (1) and V (3) have changed from service instance 1 to 3, while V (7) have changed from service instance 1 to 4. The incurred cost for the http service is q http = 12+12=24 units and for the ssh service is q ssh = 34 units. While the value of the k0d metric only rises from 1 to 2, the budget constraints are satisfied.As seen from the above test cases, our model and problem formulation makes it relatively easy to apply any standard optimization techniques, such as the GA, to optimize the k0d metric through diversity while dealing with different budget constraints.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Processing time.

Fig. 10 .

 10 Fig. 10. Average gain based on relative distance of shortest path.

Fig. 11 .

 11 Fig. 11. Average gain based on directly reachable vulnerabilities.

  Definition 2 (Extended Resource Graph). Given a network composed of a set of hosts H, a set of services S, with the service mapping serv(.) : H → 2 S , the collection of service pools SP = {sp(s) | s ∈ S}, and the labelling function v(.) : E → SP , which satisfies ∀h s ∈ S ∀h s ∈ S, v( s, h s , h d ) = v( s, h s , h d ) (meaning all exploits with common service and destination host must be associated with the same service instance, as explained earlier). Let E be the set of zero day exploits { s, h s

If different software are considered likely to share common vulnerabilities, a similaritysensitive diversity metric may be needed[START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF]

, 26].

Acknowledgements. This work was partially supported by the National Institute of Standards and Technology grant 60NANB15D091, by the National Science Foundation grant IIP-1266147, and by the Natural Sciences and Engineering Research Council of Canada under Discovery Grant N01035.

All the test cases described above rely on the assumption that all the attack paths are readily available. However, this is not always the case in practice. Due to the well known complexity that resource graphs have inherited from attack graphs due to their common syntax [START_REF] Wang | Modeling network diversity for evaluating the robustness of networks against zeroday attacks[END_REF][START_REF] Zhang | Network diversity: A security metric for evaluating the resilience of networks against zero-day attacks[END_REF], it is usually computationally infeasible to enumerate all the available attack paths in a resource graph for large networks. Therefore, we design a heuristic algorithm to reduce the search complexity when calculating and optimizing the k0d metric by only storing the m-shortest paths at each step, as depicted in Figure 5 and detailed below.

Procedure Heuristic m-shortest Input: Extended resource graph G, v , goal condition cg, number of paths m, diversified diversity control vector, D Output: σ(cg) Method:

1. Let vlistbe any topological sort of G 5. While all vlist elements are unprocessed 6.

If c ∈ C I and c is unprocessed 7.

Let σ(c) = c 8.

Mark c as processed 9.

Else if e ∈ E (e is not processed) and (∀c ∈ C)((c, e) ∈ Rr ⇒ c is processed) 10.

Let {c ∈ C : (c, e) ∈ Rr} = {c1, c2, . . . , cn} 11.

Let α(e) = a1 ∪ a2 . . . ∪ e : ai ∈ σ(ci),

If n > m 13.

Let

σ(e) = a1 ∪ a2 . . . ∪ e : ai ∈ σ(ci), 1 ≤ i ≤ m 16.

Mark e as processed 17.

Else (c s.t. (e, c) ∈ Ri and c is unprocessed) [START_REF] Wang | Network Hardening: An Automated Approach to Improving Network Security[END_REF].

If (∀e ∈ E)((e , c) ∈ Ri ⇒ e is processed) 19.

Let


Christof Ferreira Torres

Rolando Trujillo-Rasua
email: rolando.trujillo@uni.lu

The Fréchet/Manhattan distance and the trajectory anonymisation problem

Mobile communication has grown quickly in the last two decades. Connections can be wirelessly established from almost any habitable place in the earth, leading to a plethora of connection-based tracking mechanisms, such as GPS, GSM, RFID, etc. Trajectories representing the movement of people are consequently being gathered and analysed in a daily basis. However, a trajectory may contain sensitive and private information, which raises the problem of whether spatio-temporal data can be published in a private manner. In this article, we introduce a novel distance measure for trajectories that captures both aspect of the microaggregation process, namely clustering and obfuscation. Based on this distance measure we propose a trajectory anonymisation heuristic method ensuring that each trajectory is indistinguishable from k -1 other trajectories. The proposed distance measure is loosely based on the Fréchet distance, yet it can be computed efficiently in quadratic time complexity. Empirical studies on synthetic trajectories show that our anonymisation approach improves previous work in terms of utility without sacrificing privacy.

Introduction

Not long ago, visual identification was the only mean to collect spatio-temporal data from people. Nowadays this task is far easier since there is no need of direct human intervention for monitoring and tracking. Instead, surveillance cameras, social networks, credit card transactions, and many other worldwide adopted technologies and services, automatically collect this type of data. Today's pervasiveness of location-aware devices like mobile phones and GPS receivers helps even further companies and governments to easily collect huge amount of information about people's movement.

Analysing and mining this type of information, also known as trajectories, might reveal new trends and previously unknown knowledge to be used in traffic, sustainable mobility management, urban planning and supply chain management. By doing so, resources might be optimised and business and government decisions can be solid and well-founded. In this sense, both companies and citizens profit directly from the publication and analysis of databases of trajectories.

Despite of all these benefits, there are obvious threats to people's privacy if their movement data are published in a way which allows re-identification of the person behind a trajectory. Just considering the locations visited by a trajectory, it may reveal sensitive information about users like religious, political, or sexual preferences. The privacy threat grows when the time information exposes user's habits that may be used for unauthorized advertisement and user profiling.

A tentative solution to achieve anonymity is de-identification by means of removing identifying attributes of individuals. However, this is often insufficient to preserve privacy due to other type of attributes called quasi-identifiers, which are non-identifying attributes that together with external information might uniquely identify the individual behind a record. Unfortunately, in the case of spatio-temporal data, every location can be regarded as a quasi-identifier [START_REF] Yarovoy | Anonymizing moving objects: how to hide a mob in a crowd[END_REF]. Therefore, just knowing some locations visited by an individual could be enough to identify his trajectory in a database. As an example, let's consider a GPS application recording trajectories of citizens. Daily routine indicates that an early morning trajectory is likely to begin at the user's home and end at the user's workplace. This simple assumption might be enough to accurately re-identify a user's trajectory.

The above problem has been addressed relying on k-anonymity [START_REF] Samarati | Protecting privacy when disclosing information: k-anonymity and its enforcement through generalization and suppression[END_REF][START_REF] Sweeney | k-anonymity: a model for protecting privacy[END_REF], a widely used privacy notion. A set S is said to satisfy k-anonymity if each combination of quasi-identifier attribute values is shared by at least k records in S. Therefore, considering that all identifying attributes have been removed, kanonymity ensures that no anonymised record can be correctly linked to an individual with probability higher than 1/k. In microdata, the set of quasi-identifiers is typically considered small and known in advance. In spatio-temporal data, however, a similar assumption can hardly hold; any location can be regarded as a quasi-identifier. As a result, anonymisation methods aimed at achieving k-anonymity on microdata cannot be directly applied on spatio-temporal data and vice versa.

Contributions.

In this article we propose a distance measure for trajectories specially suited for clustering and obfuscation. The distance is loosely based on the Fréchet distance [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF], yet it is efficiently computable. The novel construction has significant advantages: (i) it can deal with non-overlapping trajectories, (ii) it outputs, in addition to a distance value, a set of matching points that are exploited later in the obfuscation process, and (iii) it considers the shape of the trajectories due to the very nature of the Fréchet distance. We use the proposed distance measure as the basis of a trajectory anonymisation technique that releases datasets satisfying k-anonymity, regardless of the adversary knowledge. We show, through experiments on synthetic spatio-temporal data, that our approach outperforms previous comparable work in terms of utility.

Outline of the paper. This paper is structured as follows. Section 2 next provides related work. Section 3 introduces a novel distance measure based on the Fréchet distance and the Manhattan norm. A microaggregation-based method for trajectory anonymisation is proposed in Section 4, which is empirically evaluated in Section 5. Finally, Section 6 draws conclusions and future work.

Trajectory k-anonymity is aimed at hiding a single trajectory into a crowd of at least k -1 other trajectories. The idea is that every trajectory in the published dataset be indistinguishable from k -1 other trajectories and, as a consequence, an adversary cannot identify the individual behind a trajectory with probability higher than 1/k.

An approach to achieve k-anonymity is by means of suppression of attribute values, which is generally used on discrete and/or semantic data where perturbation methods are not well suited. One of the first suppression-based methods for trajectory anonymization is due to Terrovitis and Mamoulis [START_REF] Terrovitis | Privacy preservation in the publication of trajectories[END_REF]. They consider trajectories to be sequences of addresses taken from an address domain P. The adversary controls subsets of addresses of P, and thus his knowledge is represented as projections of original trajectories over the addresses in P that are in the adversary's knowledge. A greedy algorithm aimed at guaranteeing that no address unknown by the adversary can be linked with an user with probability higher than a given threshold is proposed in [START_REF] Terrovitis | Privacy preservation in the publication of trajectories[END_REF]. The main problem with this approach is that dealing with all possible adversary's knowledge becomes harder than the original k-anonymity problem, which is already known to be NP-Hard [START_REF] Meyerson | On the complexity of optimal k-anonymity[END_REF]. There exist other suppression-based methods in the literature, e.g., [START_REF] Chen | Privacy-preserving trajectory data publishing by local suppression[END_REF]. However, they target privacy notions different to k-anonymity.

Like Terrovitis and Mamoulis in [START_REF] Terrovitis | Privacy preservation in the publication of trajectories[END_REF], Yarovoy et al. also consider an adversary controlling a subset of user's locations or quasi-identifiers [START_REF] Yarovoy | Anonymizing moving objects: how to hide a mob in a crowd[END_REF], with the distinction that such a subset may differ for different users. Trajectory k-anonymity is defined in terms of a bipartite attack graph relating original trajectories with the anonymised trajectories. The authors propose to create anonymised groups through generalisation with respect to the joint set of quasi-identifiers from the users within the group. K-anonymity is thus achieved by creating anonymised groups such that the bipartite attack graph is symmetric and the degree of each vertex representing an anonymised trajectory is at least k. It is worth remarking that the privacy model considered in this article is different, as any user's location is regarded as a quasi-identifier.

Another generalisation-based approach was proposed by Monreale et. al [START_REF] Monreale | Movement data anonymity through generalization[END_REF]. As in [START_REF] Terrovitis | Privacy preservation in the publication of trajectories[END_REF], they ignore the time information. Therefore, k-anonymity is achieved if the generalisation of every original trajectory is a sub-trajectory of the generalisation of k -1 other trajectories. In order to preserve the utility of the original dataset, a Voronoi tessellation of the geographical area is created so that each location is transformed into the Voronoi cell that contains it. Utility is measured by simply comparing clustering results.

In [START_REF] Domingo-Ferrer | Privacypreserving publication of trajectories using microaggregation[END_REF][START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF], Domingo-Ferrer et al. propose a different approach based on microagrregation and permutation rather than on generalisation. First, they introduce a novel distance measure that consider both spatial and temporal aspects of trajectories. The distance measure is flexible enough to be used either for spatio-temporal data or time series. Based on this distance measure, the authors propose to create clusters of trajectories so as to minimise the intra-cluster distance. Within a given cluster, locations are randomly swapped with other k -1 unswapped close locations. Locations that cannot be swapped are removed and so are the trajectories without swapped locations.

Abul, Bonchi, and Nanni [START_REF] Abul | Never walk alone: uncertainty for anonymity in moving objects databases[END_REF][START_REF] Abul | Anonymization of moving objects databases by clustering and perturbation[END_REF] proposed two trajectory anonymisation methods: Never Walk Alone (NWA) and Wait For Me (W4M). Both are partially based on microggregation [START_REF] Domingo-Ferrer | Practical data-oriented microaggregation for statistical disclosure control[END_REF]. The microggregation technique works as follows. The dataset of trajectories is partitioned into several clusters of size at least k and at most 2k -1. To do so, NWA relies on the Euclidean distance while W4M uses on the edit distance on real sequences (EDR) [START_REF] Lei Chen | Robust and fast similarity search for moving object trajectories[END_REF]. Trajectories within a cluster are perturbed by using space translation. The claimed privacy of these proposals has proven to be flawed [START_REF] Trujillo-Rasua | On the privacy offered by (k, δ)-anonymity[END_REF], though.

In [START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF][START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF], Nergiz et al. consider a trajectory to be a sequence of square geographical areas where a user moves randomly within a given time frame. For clustering, the authors use the log cost metric that balances the spatial and temporal distortion with user-provided weights. Since the log cost metric is based on point matching, the anonymisation process is directly inferred from the clustering process, which improves efficiency.

Recently, Gao et al. proposed a privacy-preserving technique that does not target trajectory k-anonymity directly, as most previous work do, but a tradeoff between privacy and utility [START_REF] Gao | Balancing trajectory privacy and data utility using a personalized anonymization model[END_REF]. Privacy is measured in terms of anonymity sets that are created based on a similarity measure that takes the angles and directions of the trajectories into account. Utility relies on the classical Euclidean distance.

In the literature we can find a variety of distance measures for trajectories and time-series. Vlachos et al. proposed two distance measures based on the Longest Common Subsequence problem (LCSS) [START_REF] Vlachos | Robust similarity measures for mobile object trajectories[END_REF]. The first one matches only points that are within a given spatio-temporal region. Unmatched points are discarded and taken as outliers. This criterion for outliers detection is smoothed in their second distance measure by using a weighted matching function that considers the distance between points. Another distance measure that has been designed to cope with noise is the Edit Distance on Real sequences (EDR) [START_REF] Lei Chen | Robust and fast similarity search for moving object trajectories[END_REF]. The problem is that it requires a fixed and global distance threshold that defines whether a location is too far from another location. A survey on distance measures for trajectory clustering can be found in [START_REF] Zhang | Comparison of similarity measures for trajectory clustering in outdoor surveillance scenes[END_REF].

A distance measure for trajectory microaggregation

We consider trajectories describing the movements of objects on the surface of the earth. Even though a movement is assumed to be continuous, it is typically described by a finite polyline. Formally, a trajectory is defined as a sequence of time-stamped locations τ

= 1 • • • n such that i .t < i+1 .t ∀i ∈ {1, • • • , n -1}
where .t, .x, and .y, denote the time, latitude, and longitude of the location , respectively. In general, trajectories can be recorded at different and irregular sampling rates, are not noise-free, and the velocity between two consecutive locations is assumed to be constant. A collection of trajectories is called a spatio-temporal database. For large databases, the size of a trajectory is considered to be significantly smaller than the size of the database. The choice of the distance measure is critical in microaggregation. It influences the way trajectories are clustered and usually it also impacts on the anonymisation process. There exist different factors that characterise a trajectory distance measure. For example, a distance measure may consider only trajectories within a given timespan, or look for spatial similarity regardless of direction and sampling rate, or take into account trajectory's features such as speed and angle.

The distance measure we propose in this article is loosely based on the Fréchet distance [START_REF] Fréchet | Sur quelques points du calcul functionnel [On some points of functional calculus[END_REF]. The Fréchet distance, also known as the dog-leash distance, assumes that a person walks over one trajectory and his dog over the other trajectory. Both may travel at independent but positive speed. The Fréchet distance outputs the minimum-length leash required for that person to walk his dog. Intuitively, the shorter the leash the closer the two curves.

Alt and Godau proposed in 1995 an algorithm to compute the Fréchet distance for two polylines [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF] with computational complexity O(pq log(p+q)) where p and q are the size of the polylines. To the best of our knowledge, this computational complexity has not been improved significantly without making assumptions on the curves. We thus consider variations of the Fréchet distance such as the Coupling distance [START_REF] Eiter | Computing Discrete Fréchet Distance[END_REF] and the Dynamic Time Warping (DTW) distance [START_REF] Yi | Efficient retrieval of similar time sequences under time warping[END_REF], which are significantly simpler and runs in O(pq) time complexity.

We say that a sequence

L = (u a1 , v b1 ) • • • (u an , v bn ) is a coupling between two trajectories U = u 1 • • • u p and V = v 1 • • • v q if the
following conditions are satisfied:

a 1 = 1 and b 1 = 1 a n = p and b n = q -For every i ∈ {1, • • • , n -1} it holds that a i+1 = a i or a i+1 = a i + 1, and

b i+1 = b i or b i+1 = b i + 1
A coupling can be seen as sequence of matching points, as defined in the Edit Distance on Real sequences (EDR) [START_REF] Lei Chen | Robust and fast similarity search for moving object trajectories[END_REF]. The difference, however, is that a coupling respects the order of the locations and also ensures that all points are considered.

Definition 1 (Coupling distance). Let U = u 1 • • • u p and V = v 1 • • • v q
be two trajectories and let L be the set of all couplings between U and V . Let . denote a norm on L. The coupling distance is defined as follows:

coupling dist(U, V ) = min{ L |L ∈ L}
The coupling distance can be computed by a simple dynamic algorithm. The norm that directly relates to the original discrete Fréchet distance is the Infinite norm. Given L = (u a1 , v b1 ) • • • (u an , v bn ), the Infinite norm L ∞ is the longest distance between a pair of linked locations in L, i.e., L ∞ = max i∈{1,••• ,n} d(u ai , v bi ). Another relevant norm, which we use in this article, is the Manhattan norm, defined as

L 1 = i∈{1,••• ,n} d(u ai , v bi ).
Using the Infinite norm in the coupling distance has a clear interpretation in microaggregation of trajectories, that is, the longest distance that ought to be covered in order to spatially translate a trajectory into another one. However, accounting for the longest distance may lead to non-robust behaviors, because small variations in the trajectories can cause large variations in the distance function. For this reason, we propose to use the infinity norm to compute the optimal coupling between trajectories, yet we consider the average Manhattan norm to represent the actual distance between them. We claim that the average Manhattan norm approximates better the required distortion to microaggregate trajectories. Formally, the distance measure used in the present article is defined as follows.

Definition 2 (Fréchet/Manhattan coupling distance). Let U = u 1 • • • u p and V = v 1 • • • v q be
two trajectories and let L be the set of all couplings between U and V . Let L ⊆ L such that for every l ∈ L it holds that l ∞ is minimum amongst the couplings in L. The average coupling distance is defined as:

min l∈L 1 |l| l 1
Computing the Fréchet/Manhattan distance is a bit more elaborated than computing the coupling distance. Nevertheless, it can still be computed in O(pq) time complexity as shown by Algorithm 1. Given two trajectories U = u 1 • • • u p and V = v 1 • • • v q , we create a matrix I of size p × q where we store the optimal coupling with respect to the Infinite norm. Such computation is performed by the standard dynamic approach proposed in [START_REF] Eiter | Computing Discrete Fréchet Distance[END_REF]. In order to determine the Fréchet/Manhattan distance, we consider another matrix M where we store the optimal coupling distance with respect to the Manhattan norm among those optimal couplings with respect to the Infinite norm. To do so, we need to find where those optimal couplings with respect to the Infinite norm come from. Let us analyse what is the impact of having the pair (u x , v y ) in an optimal coupling l. First, we should notice that if (u

x , v y ) ∈ l then l ∞ ≥ d(u x , v y ). Indeed, if d(u x , v y ) < min{[x-1, y], [x-1, y -1], [x, y -1]} then l ∞ = min{[x-1, y], [x- 1, y -1], [x, y -1]}, otherwise l ∞ = d(u x , v y ).
We thus store in a set C all pairs that lead to an optimal coupling with respect to the Infinite norm amongst the pairs

{[x -1, y], [x -1, y -1], [x, y -1]}. Finally, M [x, y] is computed as min{M [x, y]/L[x, y]|[x, y] ∈ C} where L[x,
y] is the size of the optimal coupling with respect the Manhattan norm for the subtrajectories

u 1 • • • u x and v 1 • • • v y .

A microaggregation-based approach

The anonymisation method proposed in this article is based on k-microaggregation, that is, a process whereby clusters of at least k homogeneous trajectories are anonymised independently. A usual homogeneity criterion is the sum of squared

I[i, 1] = max{I[i -1, 1], d(ui, v1)} 7: M [i, 1] = M [i -1, 1] + d(ui, v1) 8:
L[i, 1] = i 9: end for 10: for j = 2 to q do 11:

I[1, j] = max{I[1, j -1], d(u1, vj)} 12: M [1, j] = M [1, j -1] + d(u1, vj) 13:
L[1, j] = j 14: end for 15: for i = 2 to p do 16:

for j = 2 to q do 17:

Let R be the set end for 30: end for 31: return M [p, q]/L [p, q] pairwise distances between trajectories within a cluster (intra-cluster distance). Hence, an optimal microaggregation can be intuitively defined as the one maximising the within-groups homogeneity.

{[i -1, j], [i -1, j -1], [i, j -1]} 18: Let C ⊆ R such that for every [x, y] ∈ C it holds that I[x, y] ≤ d(ui, vj) 19: if C is empty then 20: Let [a, b] ∈ R such that for every [x, y] ∈ R it holds that I[a, b] ≤ I[x,
The optimal microaggregation problem for multivariate points, like trajectories, has proven to be NP-hard [START_REF] Oganian | On the complexity of optimal microaggregation for statistical disclosure control[END_REF]. That justifies the use of heuristics in microaggregation-based approaches for trajectory anonymisation [START_REF] Trujillo-Rasua | Privacy in spatio-temporal databases: A microaggregation-based approach[END_REF]. An additional challenge to be addressed is that distance measures between trajectories tend to be computationally expensive. This implies that computing all pairwise distances between trajectories in a large database is may not be feasible.

Below, we detail the two main components of our microaggregation-based approach, namely the proposed heuristic for trajectory clustering and the obfuscation technique.

Clustering

We use a greedy approach to address the k-microaggregation problem explained above. Each cluster is represented by a pivot trajectory, and contains k -1 other trajectories that are close to the pivot trajectory. In other words, we consider as homogeneity criterion the sum of squared distances between the pivot trajectory and the other trajectories in the cluster. Once a cluster C is created from a pool D of trajectories, all trajectories in C are removed from D. As shown by Algorithm 2, this process is repeated until D contains less than k trajectories.

Algorithm 2 Trajectory clustering

Require: D = {τ1, . . . , τN } a set of trajectories; a distance measure d : D × D → R; a natural number δ representing the number of clusters generated at each iteration; and an anonymization parameter k 1: Let C be an empty set of clusters of trajectories 2: while |D| ≥ k do 3:

Let τ 1 be a random trajectory in D 4:

Let τ δ be the farthest trajectory to τ 1 with respect to d 5:

Let

τ 1 , τ 2 , • • • , τ δ be δ trajectories in D that minimizes the sum of squares i∈[1..δ-1] d(τ i+1 , τ i ) 2 6:
Let C0 be an empty set of trajectories and d0 = ∞ 7:

for all i = 1 to δ do 8:

Create the cluster of trajectories Ci containing τ i and the closest k -1 trajectories to τ i 9:

Compute di = τ ∈C i d(τ i , τ ) 2 10:

if di < d0 then 11: C0 = Ci and d0 = di 12:

end if 13:

end for 14:

C = C ∪ {C0} 15:
Remove all trajectories in C0 from D 16: end while 17: return C

In our approach, depicted in Algorithm 2, finding the optimal set of clusters is equivalent to finding the optimal sequence of pivot trajectories. The most effective greedy solution to this problem is to choose the best cluster amongst the |D| clusters that can be created considering each trajectory in D a pivot trajectory. However, that requires the computation of all pairwise distances between the trajectories in D. As a trade-off, given a natural number δ |D|, we choose a random trajectory τ 1 in D and find the sequence τ 1 , τ 2 , • • • τ δ such that: i) τ δ is the farthest trajectory to τ 1 and ii) the sum of squares i∈[1..δ-1] d(τ i+1 , τ i ) 2 is minimum. We thus choose the best cluster amongst the δ clusters that can be built considering either τ 1 , or τ 2 , • • • , or τ δ , as pivot. Note that, if δ = |D| then we actually find the optimal set of clusters.

Obfuscation technique

Our privacy-preserving method for the publication of trajectories is based on the clustering technique and the Fréchet/Manhattan coupling distance described above. Even though the coupling distance deals well with trajectories recorded at different sampling rates, the lower is the sampling rate the better it approximates the classical Fréchet distance. We thus use linear interpolation to decrease and homogenise the sampling rate of two trajectories as follows. . An analogous procedure is used to increase the sampling rate of U with respect to V . Note that, the trajectories U and V resulting from re-sampling U and V , respectively, have equal size.

Let U = u 1 • • • u p and V = v 1 • • • v q be two trajectories. For every i ∈ {1, • • • , p},
We use the Fréchet/Manhattan distance in Algorithm 2 to partition a collection of trajectories

{τ 1 , • • • , τ N } into a set of homogeneous clusters {C 1 , • • • , C m }.
For every i ∈ {1, • • • , m}, let X be the pivot trajectory in the cluster C i as considered in Algorithm 2. For each Y ∈ C i (X = Y ), let (u a1 , u b1 ) • • • (u an , v bn ) be the optimal coupling between X and Y with respect to the Fréchet/Manhattan distance, where U = u 1 • • • u p and V = v 1 • • • v q are the re-sampling of X and Y , respectively. For each j ∈ {1, • • • , n} and if u aj ∈ X, i.e., if u aj is an original location of X rather than an interpolated location added during the re-sampling procedure, we add to the set S(u ai ) the location v bi . Once this process is finished for all trajectories in C i , we consider, for every location x ∈ X, the set S(x) containing those locations from other trajectories in C i that formed a pair with x in an optimal coupling. We always include x into S(x) whenever S(x) is not empty. The anonymised trajectory for the cluster C i will be that formed by the average locations obtained from the sets {S(x)|x ∈ X, S(x) = ∅}. A pseudo-code description of this procedure is given in Algorithm 3.

Privacy analysis

Several notions of trajectory k-anonymity exist. For example, in [START_REF] Monreale | Movement data anonymity through generalization[END_REF][START_REF] Terrovitis | Privacy preservation in the publication of trajectories[END_REF], the adversary ignores the time dimension. In [START_REF] Abul | Never walk alone: uncertainty for anonymity in moving objects databases[END_REF][START_REF] Abul | Anonymization of moving objects databases by clustering and perturbation[END_REF], an adversary is considered unable to distinguish two locations if their distance is below a predefined threshold. In [START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF], the model is defined considering that original locations must be preserved, which means that random spatial distortion is disallowed.

In this article we consider trajectory k-anonymity as a property of the anonymised dataset regardless the adversary capabilities. Our notion of k-anonymity is indeed similar to that presented in [START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF][START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF] for generalised trajectories.

Algorithm 3 Trajectory anonymization algorithm

Require: {τ1, . . . , τN } a collection of original trajectories; a number δ to be used in the clustering process; the Fréchet/Manhattan distance d; an anonymisation parameter k 1: Use the clustering technique defined by Algorithm 2 on input {τ1, . . . , τN }, the distance measure d, δ, and k, to obtain a set of clusters {C1, • • • , Cm} 2: Let D be an empty set of trajectories 3: for i = 1 to m do 4:

Let X be the pivot trajectory in Ci as defined in Algorithm 2 5:

Let S(x) be an empty set for every x ∈ X 6:

for Y ∈ Ci and X = Y do 7:

Let (ua if ua j ∈ X then 10:

1 , v b 1 ) • • • (ua n , v bn )
S(ua j ) = S(ua j ) ∪ {v b j } 11:

end if 12:

end for 13:

end for 14:

Let τ be an empty trajectory 15:

for x ∈ X and S(x) = ∅ do 16: Proof. The proof trivially follows from the fact that Algorithm 3 produces k equal trajectories for each cluster (see Step 19 in Algorithm 3).

S(x) = S(x) ∪

Empirical evaluation

As the privacy-preserving anonimisation technique introduced in this article replaces a cluster of k close, but potentially different, trajectories by k identical trajectories, it is of paramount importance to evaluate utility loss in this method. Next, we introduce spatial-range queries as a measure of utility loss. We finally compare our anonymisation approach with other state-of-the-art privacy preserving techniques.

Trajectory analysis and utility measures

There exist a plethora of trajectory analysis techniques developed within the Geographic Information Science and Data Mining fields. These techniques may look for movement patterns such as flocking, leadership, commuting, and encounter, or may be aimed at answering basic queries such as nearest neighbor or range queries.

In this article we mainly focus on queries that are used for aggregate statistics. This queries are typically measurable, and thus they can be defined as functions on the domain of all spatio-temporal databases ranging over a metric space. Let D be the universe of all possible collections of trajectories and let (M, d) be a metric space. A spatio-temporal query Q is formally defined as a function Q : D → M . Examples of measurable queries are traffic density, travel time, peak hours, amongst many others.

Measurable queries can be naturally used to define utility measures for anonymization techniques as follows. Let D ∈ D be an original spatio-temporal database and D ∈ D its anonymized version. Given a measurable query Q : D → M , we measure utility loss by the formula d(Q(D), Q(D )). The closer this measure to zero the better D approximates D with respect to Q.

A well-known type of measurable query in trajectory analysis is spatiotemporal range queries, which were introduced by Trajcevski et al. in [START_REF] Trajcevski | Managing uncertainty in moving objects databases[END_REF] in 2004. In particular, we consider the two following queries.

-Sometime Definitely Inside(T , R, t b , t e ) is true if and only if there exists a time t ∈ [t b , t e ] at which trajectory T is inside region R. -Always Definitely Inside(T , R, t b , t e ) is true if and only if at every time t ∈ [t b , t e ], trajectory T is inside region R.

At a first sight, it may seem that the query Always Definitely Inside(AI) is stronger than Sometime Definitely Inside(SI). However, with the later we can formulate questions at a local level like: how many users pass through the Grand Place in Belgium?, whilst with AI the shape of trajectories becomes more relevant and might be useful for questions like: how many users take the toll highway placed between Barcelona and Tarragona cities?

Other important points to be remarked are the area of R and the time interval [t b , t e ]. Both provide flexibility when dealing with uncertain or perturbed trajectories. Asking for trajectories passing through a single location at a given time-stamp is meaningless in this type of imprecise data. The size of the area and the time interval should not be too large either, though.

Similarly to [START_REF] Abul | Never walk alone: uncertainty for anonymity in moving objects databases[END_REF][START_REF] Abul | Anonymization of moving objects databases by clustering and perturbation[END_REF][START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF], we used both queries to define a distortion metric of the anonymised dataset T with respect to original dataset T . The idea is to define a large set of queries according to some distribution of regions and time intervals. The same set of queries is applied to both datasets T and T and the number of trajectories satisfying SI and AI are counted as shows the following SQL style code.

-Query Q 1 (T , R, t b , t e ):

SELECT COUNT (*) FROM T WHERE SI(T .traj, R, t b , t e ) -Query Q 2 (T , R, t b , t e ):

SELECT COUNT (*) FROM T WHERE AI(T .traj, R, t b , t e )

Two different range query distortions SID(T , T ) and AID(T , T ) are defined by using the accumulative queries Q 1 and Q 2 , respectively.

-SID(T , T ) = 1 Both metrics SID and AID are bounded by 0 and 1. The minimum is achieved when Q i (T , R, t b , t e ) = Q i (T , R, t b , t e ), and the maximum if Q i (T , R, t b , t e ) = 0, where i ∈ {1, 2}. Therefore, the lower the range query distortion the lower the utility loss of the anonymised dataset.

Implementation details of the considered methods

We compare our method with the generalisation-based and permutation-based approach proposed in [START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF] and [START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF], respectively. The generalisation-based method relies on a distance threshold, which allows the Log-cost distance measure to discard outlier locations. Because in this section we only consider noiseless synthetic data, we have set up such a distance threshold to its maximum value. The permutation-based method, instead, discard outlier locations during the obfuscation process by considering both a distance and a time threshold. Again, we set up both thresholds to their maximum values so as to avoid outlier removal in a noiseless dataset. The permutation-based method considered in this article is the one named SwapLocations in [START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF].

Results on synthetic trajectories

We compare our anonymisation method with other approaches by using a synthetic dataset generated with Brinkhoff's framework [START_REF] Brinkhoff | A framework for generating network-based moving objects[END_REF], which is used often to evaluate privacy-preserving approaches. Synthetic data generated with Brinkhoff's generator have the advantage of being easily transferable and reproducible. We thus provide next the parameters used to generate the dataset of trajectories considered in our experiments.

The generation parameters over the map of Oldenburg were: 6 moving object classes and 3 external object classes; 5 moving objects and 3 external object generated per time-stamp; the maximum lifespan of a trajectory was set up to 1, 000 time-stamps; speed 10; and report probability 1, 000. This resulted in 5, 000 synthetic trajectories provided by Brinkhoff's generator [START_REF] Brinkhoff | A framework for generating network-based moving objects[END_REF], which contain a total of 492, 105 locations in the German city of Oldenburg and 98.421 locations per trajectory in average.

In order to generate spatial-range queries, we considered regions whose radius randomly distributes over the interval of natural numbers [0, 500]. The maximum of this interval is a small fraction of the average length of each trajectory, which is 7284. Remark that the smaller the spatial interval the tighter is the spatial-range query and the harder become for an anonymisation technique to apply spatial distortion without bringing down utility. We respect to the time dimension we considered different time intervals [0, 0], [0, 300], [0, 600], [0, 1800], [0, 3600]. For a given time interval [0, t], we generate a spatial-range query by choosing a random interval [t b , t e ] such that 0 ≤ t b ≤ t e ≤ t.

We generated for each time interval 100, 000 spatial-range queries of both types: Q 1 and Q 2 . Armed with these set of queries, we computed the range query distortions SID and AID of the anonymised data sets provided by three different anonymisation methods: the Generalisation-based approach [START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF], the Permutation-based approach [START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF], and our method. Each anonymisation method provided three different datasets satisfying k-anonymity with k ∈ {2, 4, 8}. The results are depicted in Figure 1.

It can be seen from Figure 1 that our method performs better than the approaches proposed in [START_REF] Mehmet Ercan Nergiz | Towards trajectory anonymization: a generalization-based approach[END_REF][START_REF] Domingo-Ferrer | Microaggregation-and permutation-based anonymization of movement data[END_REF] for every cluster size and every time interval. The improvement in terms of utility increases as the offered privacy increases. For k = 2, our method is just slightly better than the generalisation-based approach, while for k ∈ {4, 8} our method performs significantly better. This means that our technique clusters and anonymises trajectories more efficiently.

Figure 1 also shows that more research on trajectory anonymisation techniques ought to be conducted. The ideal range query distortion is zero, and none of the three considered techniques gets close to this optimal value. This issue can be overcome by considering larger datasets of original trajectories. Intuitively, the larger the dataset the easier is to find clusters with low intra-cluster distance. Other solution approach consists in removing outlier trajectories, that is, trajectories that cannot be clustered with other k -1 trajectories without dramatically increasing the intra-cluster distance. The study and evaluation of these solution approaches, as well as reporting on results over real-life datasets, are left as future work.

Conclusions

In this article we have introduced a novel distance measure for trajectories, which is well suited for both clustering and anoymisation. The proposed distance measure resembles to other types of coupling distance measures, such as the Fréchet distance, with the particularity that the Infinite norm and the Manhattan norm are considered together. To demonstrate the suitability of our distance measure, we presented a trajectory-anonymisation heuristic method that creates cluster with low intra-cluster distance and satisfies trajectory k-anonymity. Empirical results show that our method offers better utility than other state-of-the-art methods, such as the generalisation-based and permutation-based approaches.

Algorithm 1

 1 Average coupling distance Require: Two trajectories U = u1 • • • up and V = v1 • • • vq 1: Let I, M , L be three matrices of size p × q. Intuitively, I and M represent the Infinity and Manhattan norms, respectively, while L is the length of the optimal coupling 2: Let d represent the Euclidean distance. 3: I[1, 1] = M [1, 1] = d(u1, v1) 4: L[1, 1] = 1 5: for i = 2 to p do 6:

y] 21 :

 21 I[i, j] = I[a, b] 22: Add to C every element [x, y] in R such that I[x, y] = I[a, b] 23: else 24: I[i, j] = d(ui, vj) 25: end if 26: Let [x, y] ∈ C such that M [x, y]/L[x, y] ≤ M [a, b]/L[a, b] for every [a, b] ∈ C 27: M [i, j] = M [x, y] + d(ui, vj) 28: L[i, j] = L[x, y] + 1 29:

  be the optimal coupling between X and Y with respect to the Fréchet/Manhattan distance d, where U = u1 • • • up and V = v1 • • • vq are the re-sampling of X and Y , respectively 8: j = 1 to n do 9:

{x} 17 :} 20 :

 1720 Add to τ the average location formed by the locations in S(x) 18: end for 19: D = D ∪ {τ, . . . , τ k end for 21: return D Definition 3 (Trajectory k-anonymity). Let D be a collection of trajectories. D meets trajectory k-anonymity if every trajectory in D is equal to other k -1 trajectories in D . Theorem 1. Let D be a collection of original trajectories and D the output of Algorithm 3 on input D. D satisfies trajectory k-anonymity.

  |ξ|∀<R,t b ,te>∈ξ|Q1(T ,R,t b ,te)-Q1(T ,R,t b ,te)| max (Q1(T ,R,t b ,te),Q1(T ,R,t b ,te))where ξ is a large set of SI queries.-AID(T , T ) = 1 |ξ| ∀<R,t b ,te>∈ξ |Q2(T ,R,t b ,te)-Q2(T ,R,t b ,te)| max (Q2(T ,R,t b ,te),Q2(T ,R,t b ,te))where ξ is a large set of AI queries.

Fig. 1 .

 1 Fig. 1. Six charts showing the range query distortions of three different anonymisation methods. Charts on the left depict the SI query distortion (SID), while charts on the right show the AI query distortion (AID).

  we insert in V by using linear interpolation a new point at time v 1 .t +

	(vq.t-v1.t)(ui.t-u1.t)
	up.t-u1.t

Future work will be directed towards reaching optimal range-query distortion values.